    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.math3.optimization.fitting;
    019
    020    import java.util.Arrays;
    021    import java.util.Comparator;
    022
    023    import org.apache.commons.math3.analysis.function.Gaussian;
    024    import org.apache.commons.math3.exception.NullArgumentException;
    025    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    026    import org.apache.commons.math3.exception.OutOfRangeException;
    027    import org.apache.commons.math3.exception.ZeroException;
    028    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    029    import org.apache.commons.math3.exception.util.LocalizedFormats;
    030    import org.apache.commons.math3.optimization.DifferentiableMultivariateVectorOptimizer;
    031    import org.apache.commons.math3.util.FastMath;
    032
    033    /**
    034     * Fits points to a {@link
    035     * org.apache.commons.math3.analysis.function.Gaussian.Parametric Gaussian} function.
    036     * <p>
    037     * Usage example:
    038     * <pre>
    039     *   GaussianFitter fitter = new GaussianFitter(
    040     *     new LevenbergMarquardtOptimizer());
    041     *   fitter.addObservedPoint(4.0254623,  531026.0);
    042     *   fitter.addObservedPoint(4.03128248, 984167.0);
    043     *   fitter.addObservedPoint(4.03839603, 1887233.0);
    044     *   fitter.addObservedPoint(4.04421621, 2687152.0);
    045     *   fitter.addObservedPoint(4.05132976, 3461228.0);
    046     *   fitter.addObservedPoint(4.05326982, 3580526.0);
    047     *   fitter.addObservedPoint(4.05779662, 3439750.0);
    048     *   fitter.addObservedPoint(4.0636168,  2877648.0);
    049     *   fitter.addObservedPoint(4.06943698, 2175960.0);
    050     *   fitter.addObservedPoint(4.07525716, 1447024.0);
    051     *   fitter.addObservedPoint(4.08237071, 717104.0);
    052     *   fitter.addObservedPoint(4.08366408, 620014.0);
    053     *   double[] parameters = fitter.fit();
    054     * </pre>
    055     *
    056     * @since 2.2
    057     * @version $Id: GaussianFitter.java 1422230 2012-12-15 12:11:13Z erans $
    058     * @deprecated As of 3.1 (to be removed in 4.0).
    059     */
    060    @Deprecated
    061    public class GaussianFitter extends CurveFitter<Gaussian.Parametric> {
    062        /**
    063         * Constructs an instance using the specified optimizer.
    064         *
    065         * @param optimizer Optimizer to use for the fitting.
    066         */
    067        public GaussianFitter(DifferentiableMultivariateVectorOptimizer optimizer) {
    068            super(optimizer);
    069        }
    070
    071        /**
    072         * Fits a Gaussian function to the observed points.
    073         *
    074         * @param initialGuess First guess values in the following order:
    075         * <ul>
    076         *  <li>Norm</li>
    077         *  <li>Mean</li>
    078         *  <li>Sigma</li>
    079         * </ul>
    080         * @return the parameters of the Gaussian function that best fits the
    081         * observed points (in the same order as above).
    082         * @since 3.0
    083         */
    084        public double[] fit(double[] initialGuess) {
    085            final Gaussian.Parametric f = new Gaussian.Parametric() {
    086                    @Override
    087                    public double value(double x, double ... p) {
    088                        double v = Double.POSITIVE_INFINITY;
    089                        try {
    090                            v = super.value(x, p);
    091                        } catch (NotStrictlyPositiveException e) { // NOPMD
    092                            // Do nothing.
    093                        }
    094                        return v;
    095                    }
    096
    097                    @Override
    098                    public double[] gradient(double x, double ... p) {
    099                        double[] v = { Double.POSITIVE_INFINITY,
    100                                       Double.POSITIVE_INFINITY,
    101                                       Double.POSITIVE_INFINITY };
    102                        try {
    103                            v = super.gradient(x, p);
    104                        } catch (NotStrictlyPositiveException e) { // NOPMD
    105                            // Do nothing.
    106                        }
    107                        return v;
    108                    }
    109                };
    110
    111            return fit(f, initialGuess);
    112        }
    113
    114        /**
    115         * Fits a Gaussian function to the observed points.
    116         *
    117         * @return the parameters of the Gaussian function that best fits the
    118         * observed points (in the same order as above).
    119         */
    120        public double[] fit() {
    121            final double[] guess = (new ParameterGuesser(getObservations())).guess();
    122            return fit(guess);
    123        }
    124
    125        /**
    126         * Guesses the parameters {@code norm}, {@code mean}, and {@code sigma}
    127         * of a {@link org.apache.commons.math3.analysis.function.Gaussian.Parametric}
    128         * based on the specified observed points.
    129         */
    130        public static class ParameterGuesser {
    131            /** Normalization factor. */
    132            private final double norm;
    133            /** Mean. */
    134            private final double mean;
    135            /** Standard deviation. */
    136            private final double sigma;
    137
    138            /**
    139             * Constructs instance with the specified observed points.
    140             *
    141             * @param observations Observed points from which to guess the
    142             * parameters of the Gaussian.
    143             * @throws NullArgumentException if {@code observations} is
    144             * {@code null}.
    145             * @throws NumberIsTooSmallException if there are less than 3
    146             * observations.
    147             */
    148            public ParameterGuesser(WeightedObservedPoint[] observations) {
    149                if (observations == null) {
    150                    throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
    151                }
    152                if (observations.length < 3) {
    153                    throw new NumberIsTooSmallException(observations.length, 3, true);
    154                }
    155
    156                final WeightedObservedPoint[] sorted = sortObservations(observations);
    157                final double[] params = basicGuess(sorted);
    158
    159                norm = params[0];
    160                mean = params[1];
    161                sigma = params[2];
    162            }
    163
    164            /**
    165             * Gets an estimation of the parameters.
    166             *
    167             * @return the guessed parameters, in the following order:
    168             * <ul>
    169             *  <li>Normalization factor</li>
    170             *  <li>Mean</li>
    171             *  <li>Standard deviation</li>
    172             * </ul>
    173             */
    174            public double[] guess() {
    175                return new double[] { norm, mean, sigma };
    176            }
    177
    178            /**
    179             * Sort the observations.
    180             *
    181             * @param unsorted Input observations.
    182             * @return the input observations, sorted.
    183             */
    184            private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) {
    185                final WeightedObservedPoint[] observations = unsorted.clone();
    186                final Comparator<WeightedObservedPoint> cmp
    187                    = new Comparator<WeightedObservedPoint>() {
    188                    public int compare(WeightedObservedPoint p1,
    189                                       WeightedObservedPoint p2) {
    190                        if (p1 == null && p2 == null) {
    191                            return 0;
    192                        }
    193                        if (p1 == null) {
    194                            return -1;
    195                        }
    196                        if (p2 == null) {
    197                            return 1;
    198                        }
    199                        if (p1.getX() < p2.getX()) {
    200                            return -1;
    201                        }
    202                        if (p1.getX() > p2.getX()) {
    203                            return 1;
    204                        }
    205                        if (p1.getY() < p2.getY()) {
    206                            return -1;
    207                        }
    208                        if (p1.getY() > p2.getY()) {
    209                            return 1;
    210                        }
    211                        if (p1.getWeight() < p2.getWeight()) {
    212                            return -1;
    213                        }
    214                        if (p1.getWeight() > p2.getWeight()) {
    215                            return 1;
    216                        }
    217                        return 0;
    218                    }
    219                };
    220
    221                Arrays.sort(observations, cmp);
    222                return observations;
    223            }
    224
    225            /**
    226             * Guesses the parameters based on the specified observed points.
    227             *
    228             * @param points Observed points, sorted.
    229             * @return the guessed parameters (normalization factor, mean and
    230             * sigma).
    231             */
    232            private double[] basicGuess(WeightedObservedPoint[] points) {
    233                final int maxYIdx = findMaxY(points);
    234                final double n = points[maxYIdx].getY();
    235                final double m = points[maxYIdx].getX();
    236
    237                double fwhmApprox;
    238                try {
    239                    final double halfY = n + ((m - n) / 2);
    240                    final double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY);
    241                    final double fwhmX2 = interpolateXAtY(points, maxYIdx, 1, halfY);
    242                    fwhmApprox = fwhmX2 - fwhmX1;
    243                } catch (OutOfRangeException e) {
    244                    // TODO: Exceptions should not be used for flow control.
    245                    fwhmApprox = points[points.length - 1].getX() - points[0].getX();
    246                }
    247                final double s = fwhmApprox / (2 * FastMath.sqrt(2 * FastMath.log(2)));
    248
    249                return new double[] { n, m, s };
    250            }
    251
    252            /**
    253             * Finds index of point in specified points with the largest Y.
    254             *
    255             * @param points Points to search.
    256             * @return the index in specified points array.
    257             */
    258            private int findMaxY(WeightedObservedPoint[] points) {
    259                int maxYIdx = 0;
    260                for (int i = 1; i < points.length; i++) {
    261                    if (points[i].getY() > points[maxYIdx].getY()) {
    262                        maxYIdx = i;
    263                    }
    264                }
    265                return maxYIdx;
    266            }
    267
    268            /**
    269             * Interpolates using the specified points to determine X at the
    270             * specified Y.
    271             *
    272             * @param points Points to use for interpolation.
    273             * @param startIdx Index within points from which to start the search for
    274             * interpolation bounds points.
    275             * @param idxStep Index step for searching interpolation bounds points.
    276             * @param y Y value for which X should be determined.
    277             * @return the value of X for the specified Y.
    278             * @throws ZeroException if {@code idxStep} is 0.
    279             * @throws OutOfRangeException if specified {@code y} is not within the
    280             * range of the specified {@code points}.
    281             */
    282            private double interpolateXAtY(WeightedObservedPoint[] points,
    283                                           int startIdx,
    284                                           int idxStep,
    285                                           double y)
    286                throws OutOfRangeException {
    287                if (idxStep == 0) {
    288                    throw new ZeroException();
    289                }
    290                final WeightedObservedPoint[] twoPoints
    291                    = getInterpolationPointsForY(points, startIdx, idxStep, y);
    292                final WeightedObservedPoint p1 = twoPoints[0];
    293                final WeightedObservedPoint p2 = twoPoints[1];
    294                if (p1.getY() == y) {
    295                    return p1.getX();
    296                }
    297                if (p2.getY() == y) {
    298                    return p2.getX();
    299                }
    300                return p1.getX() + (((y - p1.getY()) * (p2.getX() - p1.getX())) /
    301                                    (p2.getY() - p1.getY()));
    302            }
    303
    304            /**
    305             * Gets the two bounding interpolation points from the specified points
    306             * suitable for determining X at the specified Y.
    307             *
    308             * @param points Points to use for interpolation.
    309             * @param startIdx Index within points from which to start search for
    310             * interpolation bounds points.
    311             * @param idxStep Index step for search for interpolation bounds points.
    312             * @param y Y value for which X should be determined.
    313             * @return the array containing two points suitable for determining X at
    314             * the specified Y.
    315             * @throws ZeroException if {@code idxStep} is 0.
    316             * @throws OutOfRangeException if specified {@code y} is not within the
    317             * range of the specified {@code points}.
    318             */
    319            private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points,
    320                                                                       int startIdx,
    321                                                                       int idxStep,
    322                                                                       double y)
    323                throws OutOfRangeException {
    324                if (idxStep == 0) {
    325                    throw new ZeroException();
    326                }
    327                for (int i = startIdx;
    328                     idxStep < 0 ? i + idxStep >= 0 : i + idxStep < points.length;
    329                     i += idxStep) {
    330                    final WeightedObservedPoint p1 = points[i];
    331                    final WeightedObservedPoint p2 = points[i + idxStep];
    332                    if (isBetween(y, p1.getY(), p2.getY())) {
    333                        if (idxStep < 0) {
    334                            return new WeightedObservedPoint[] { p2, p1 };
    335                        } else {
    336                            return new WeightedObservedPoint[] { p1, p2 };
    337                        }
    338                    }
    339                }
    340
    341                // Boundaries are replaced by dummy values because the raised
    342                // exception is caught and the message never displayed.
    343                // TODO: Exceptions should not be used for flow control.
    344                throw new OutOfRangeException(y,
    345                                              Double.NEGATIVE_INFINITY,
    346                                              Double.POSITIVE_INFINITY);
    347            }
    348
    349            /**
    350             * Determines whether a value is between two other values.
    351             *
    352             * @param value Value to test whether it is between {@code boundary1}
    353             * and {@code boundary2}.
    354             * @param boundary1 One end of the range.
    355             * @param boundary2 Other end of the range.
    356             * @return {@code true} if {@code value} is between {@code boundary1} and
    357             * {@code boundary2} (inclusive), {@code false} otherwise.
    358             */
    359            private boolean isBetween(double value,
    360                                      double boundary1,
    361                                      double boundary2) {
    362                return (value >= boundary1 && value <= boundary2) ||
    363                    (value >= boundary2 && value <= boundary1);
    364            }
    365        }
    366    }