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    }