1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 package org.apache.commons.math.optimization.fitting; 19 20 import java.util.Arrays; 21 import java.util.Comparator; 22 23 import org.apache.commons.math.exception.util.LocalizedFormats; 24 import org.apache.commons.math.exception.NumberIsTooSmallException; 25 import org.apache.commons.math.exception.OutOfRangeException; 26 import org.apache.commons.math.exception.ZeroException; 27 import org.apache.commons.math.exception.NullArgumentException; 28 29 /** 30 * Guesses the parameters ({@code a}, {@code b}, {@code c}, and {@code d}) 31 * of a {@link ParametricGaussianFunction} based on the specified observed 32 * points. 33 * 34 * @since 2.2 35 * @version $Revision: 983921 $ $Date: 2010-08-10 12:46:06 +0200 (mar. 10 aot 2010) $ 36 */ 37 public class GaussianParametersGuesser { 38 39 /** Observed points. */ 40 private final WeightedObservedPoint[] observations; 41 42 /** Resulting guessed parameters. */ 43 private double[] parameters; 44 45 /** 46 * Constructs instance with the specified observed points. 47 * 48 * @param observations observed points upon which should base guess 49 */ 50 public GaussianParametersGuesser(WeightedObservedPoint[] observations) { 51 if (observations == null) { 52 throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY); 53 } 54 if (observations.length < 3) { 55 throw new NumberIsTooSmallException(observations.length, 3, true); 56 } 57 this.observations = observations.clone(); 58 } 59 60 /** 61 * Guesses the parameters based on the observed points. 62 * 63 * @return guessed parameters array <code>{a, b, c, d}</code> 64 */ 65 public double[] guess() { 66 if (parameters == null) { 67 parameters = basicGuess(observations); 68 } 69 return parameters.clone(); 70 } 71 72 /** 73 * Guesses the parameters based on the specified observed points. 74 * 75 * @param points observed points upon which should base guess 76 * 77 * @return guessed parameters array <code>{a, b, c, d}</code> 78 */ 79 private double[] basicGuess(WeightedObservedPoint[] points) { 80 Arrays.sort(points, createWeightedObservedPointComparator()); 81 double[] params = new double[4]; 82 83 int minYIdx = findMinY(points); 84 params[0] = points[minYIdx].getY(); 85 86 int maxYIdx = findMaxY(points); 87 params[1] = points[maxYIdx].getY(); 88 params[2] = points[maxYIdx].getX(); 89 90 double fwhmApprox; 91 try { 92 double halfY = params[0] + ((params[1] - params[0]) / 2.0); 93 double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY); 94 double fwhmX2 = interpolateXAtY(points, maxYIdx, +1, halfY); 95 fwhmApprox = fwhmX2 - fwhmX1; 96 } catch (OutOfRangeException e) { 97 fwhmApprox = points[points.length - 1].getX() - points[0].getX(); 98 } 99 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 } 272
