    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    }