001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math.optimization.fitting;
019
020 import java.util.Arrays;
021 import java.util.Comparator;
022
023 import org.apache.commons.math.exception.util.LocalizedFormats;
024 import org.apache.commons.math.exception.NumberIsTooSmallException;
025 import org.apache.commons.math.exception.OutOfRangeException;
026 import org.apache.commons.math.exception.ZeroException;
027 import org.apache.commons.math.exception.NullArgumentException;
028
029 /**
030 * Guesses the parameters ({@code a}, {@code b}, {@code c}, and {@code d})
031 * of a {@link ParametricGaussianFunction} based on the specified observed
032 * points.
033 *
034 * @since 2.2
035 * @version $Revision: 983921 $ $Date: 2010-08-10 12:46:06 +0200 (mar. 10 ao??t 2010) $
036 */
037 public class GaussianParametersGuesser {
038
039 /** Observed points. */
040 private final WeightedObservedPoint[] observations;
041
042 /** Resulting guessed parameters. */
043 private double[] parameters;
044
045 /**
046 * Constructs instance with the specified observed points.
047 *
048 * @param observations observed points upon which should base guess
049 */
050 public GaussianParametersGuesser(WeightedObservedPoint[] observations) {
051 if (observations == null) {
052 throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
053 }
054 if (observations.length < 3) {
055 throw new NumberIsTooSmallException(observations.length, 3, true);
056 }
057 this.observations = observations.clone();
058 }
059
060 /**
061 * Guesses the parameters based on the observed points.
062 *
063 * @return guessed parameters array <code>{a, b, c, d}</code>
064 */
065 public double[] guess() {
066 if (parameters == null) {
067 parameters = basicGuess(observations);
068 }
069 return parameters.clone();
070 }
071
072 /**
073 * Guesses the parameters based on the specified observed points.
074 *
075 * @param points observed points upon which should base guess
076 *
077 * @return guessed parameters array <code>{a, b, c, d}</code>
078 */
079 private double[] basicGuess(WeightedObservedPoint[] points) {
080 Arrays.sort(points, createWeightedObservedPointComparator());
081 double[] params = new double[4];
082
083 int minYIdx = findMinY(points);
084 params[0] = points[minYIdx].getY();
085
086 int maxYIdx = findMaxY(points);
087 params[1] = points[maxYIdx].getY();
088 params[2] = points[maxYIdx].getX();
089
090 double fwhmApprox;
091 try {
092 double halfY = params[0] + ((params[1] - params[0]) / 2.0);
093 double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY);
094 double fwhmX2 = interpolateXAtY(points, maxYIdx, +1, halfY);
095 fwhmApprox = fwhmX2 - fwhmX1;
096 } catch (OutOfRangeException e) {
097 fwhmApprox = points[points.length - 1].getX() - points[0].getX();
098 }
099 params[3] = fwhmApprox / (2.0 * Math.sqrt(2.0 * Math.log(2.0)));
100
101 return params;
102 }
103
104 /**
105 * Finds index of point in specified points with the smallest Y.
106 *
107 * @param points points to search
108 *
109 * @return index in specified points array
110 */
111 private int findMinY(WeightedObservedPoint[] points) {
112 int minYIdx = 0;
113 for (int i = 1; i < points.length; i++) {
114 if (points[i].getY() < points[minYIdx].getY()) {
115 minYIdx = i;
116 }
117 }
118 return minYIdx;
119 }
120
121 /**
122 * Finds index of point in specified points with the largest Y.
123 *
124 * @param points points to search
125 *
126 * @return index in specified points array
127 */
128 private int findMaxY(WeightedObservedPoint[] points) {
129 int maxYIdx = 0;
130 for (int i = 1; i < points.length; i++) {
131 if (points[i].getY() > points[maxYIdx].getY()) {
132 maxYIdx = i;
133 }
134 }
135 return maxYIdx;
136 }
137
138 /**
139 * Interpolates using the specified points to determine X at the specified
140 * Y.
141 *
142 * @param points points to use for interpolation
143 * @param startIdx index within points from which to start search for
144 * interpolation bounds points
145 * @param idxStep index step for search for interpolation bounds points
146 * @param y Y value for which X should be determined
147 *
148 * @return value of X at the specified Y
149 *
150 * @throws IllegalArgumentException if idxStep is 0
151 * @throws OutOfRangeException if specified <code>y</code> is not within the
152 * range of the specified <code>points</code>
153 */
154 private double interpolateXAtY(WeightedObservedPoint[] points,
155 int startIdx, int idxStep, double y) throws OutOfRangeException {
156 if (idxStep == 0) {
157 throw new ZeroException();
158 }
159 WeightedObservedPoint[] twoPoints = getInterpolationPointsForY(points, startIdx, idxStep, y);
160 WeightedObservedPoint pointA = twoPoints[0];
161 WeightedObservedPoint pointB = twoPoints[1];
162 if (pointA.getY() == y) {
163 return pointA.getX();
164 }
165 if (pointB.getY() == y) {
166 return pointB.getX();
167 }
168 return pointA.getX() +
169 (((y - pointA.getY()) * (pointB.getX() - pointA.getX())) / (pointB.getY() - pointA.getY()));
170 }
171
172 /**
173 * Gets the two bounding interpolation points from the specified points
174 * suitable for determining X at the specified Y.
175 *
176 * @param points points to use for interpolation
177 * @param startIdx index within points from which to start search for
178 * interpolation bounds points
179 * @param idxStep index step for search for interpolation bounds points
180 * @param y Y value for which X should be determined
181 *
182 * @return array containing two points suitable for determining X at the
183 * specified Y
184 *
185 * @throws IllegalArgumentException if idxStep is 0
186 * @throws OutOfRangeException if specified <code>y</code> is not within the
187 * range of the specified <code>points</code>
188 */
189 private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points,
190 int startIdx, int idxStep, double y)
191 throws OutOfRangeException {
192 if (idxStep == 0) {
193 throw new ZeroException();
194 }
195 for (int i = startIdx;
196 (idxStep < 0) ? (i + idxStep >= 0) : (i + idxStep < points.length);
197 i += idxStep) {
198 if (isBetween(y, points[i].getY(), points[i + idxStep].getY())) {
199 return (idxStep < 0) ?
200 new WeightedObservedPoint[] { points[i + idxStep], points[i] } :
201 new WeightedObservedPoint[] { points[i], points[i + idxStep] };
202 }
203 }
204
205 double minY = Double.POSITIVE_INFINITY;
206 double maxY = Double.NEGATIVE_INFINITY;
207 for (final WeightedObservedPoint point : points) {
208 minY = Math.min(minY, point.getY());
209 maxY = Math.max(maxY, point.getY());
210 }
211 throw new OutOfRangeException(y, minY, maxY);
212
213 }
214
215 /**
216 * Determines whether a value is between two other values.
217 *
218 * @param value value to determine whether is between <code>boundary1</code>
219 * and <code>boundary2</code>
220 * @param boundary1 one end of the range
221 * @param boundary2 other end of the range
222 *
223 * @return true if <code>value</code> is between <code>boundary1</code> and
224 * <code>boundary2</code> (inclusive); false otherwise
225 */
226 private boolean isBetween(double value, double boundary1, double boundary2) {
227 return (value >= boundary1 && value <= boundary2) ||
228 (value >= boundary2 && value <= boundary1);
229 }
230
231 /**
232 * Factory method creating <code>Comparator</code> for comparing
233 * <code>WeightedObservedPoint</code> instances.
234 *
235 * @return new <code>Comparator</code> instance
236 */
237 private Comparator<WeightedObservedPoint> createWeightedObservedPointComparator() {
238 return new Comparator<WeightedObservedPoint>() {
239 public int compare(WeightedObservedPoint p1, WeightedObservedPoint p2) {
240 if (p1 == null && p2 == null) {
241 return 0;
242 }
243 if (p1 == null) {
244 return -1;
245 }
246 if (p2 == null) {
247 return 1;
248 }
249 if (p1.getX() < p2.getX()) {
250 return -1;
251 }
252 if (p1.getX() > p2.getX()) {
253 return 1;
254 }
255 if (p1.getY() < p2.getY()) {
256 return -1;
257 }
258 if (p1.getY() > p2.getY()) {
259 return 1;
260 }
261 if (p1.getWeight() < p2.getWeight()) {
262 return -1;
263 }
264 if (p1.getWeight() > p2.getWeight()) {
265 return 1;
266 }
267 return 0;
268 }
269 };
270 }
271 }