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 }