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 */
017package org.apache.commons.math4.legacy.fitting;
018
019import java.util.List;
020import java.util.Collection;
021
022import org.apache.commons.math4.legacy.analysis.function.Gaussian;
023import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException;
024import org.apache.commons.math4.legacy.exception.NullArgumentException;
025import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
026import org.apache.commons.math4.legacy.exception.OutOfRangeException;
027import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
028import org.apache.commons.math4.core.jdkmath.JdkMath;
029
030/**
031 * Fits points to a {@link
032 * org.apache.commons.math4.legacy.analysis.function.Gaussian.Parametric Gaussian}
033 * function.
034 * <br>
035 * The {@link #withStartPoint(double[]) initial guess values} must be passed
036 * in the following order:
037 * <ul>
038 *  <li>Normalization</li>
039 *  <li>Mean</li>
040 *  <li>Sigma</li>
041 * </ul>
042 * The optimal values will be returned in the same order.
043 *
044 * <p>
045 * Usage example:
046 * <pre>
047 *   WeightedObservedPoints obs = new WeightedObservedPoints();
048 *   obs.add(4.0254623,  531026.0);
049 *   obs.add(4.03128248, 984167.0);
050 *   obs.add(4.03839603, 1887233.0);
051 *   obs.add(4.04421621, 2687152.0);
052 *   obs.add(4.05132976, 3461228.0);
053 *   obs.add(4.05326982, 3580526.0);
054 *   obs.add(4.05779662, 3439750.0);
055 *   obs.add(4.0636168,  2877648.0);
056 *   obs.add(4.06943698, 2175960.0);
057 *   obs.add(4.07525716, 1447024.0);
058 *   obs.add(4.08237071, 717104.0);
059 *   obs.add(4.08366408, 620014.0);
060 *   double[] parameters = GaussianCurveFitter.create().fit(obs.toList());
061 * </pre>
062 *
063 * @since 3.3
064 */
065public final class GaussianCurveFitter extends SimpleCurveFitter {
066    /** Parametric function to be fitted. */
067    private static final Gaussian.Parametric FUNCTION = new Gaussian.Parametric() {
068            /** {@inheritDoc} */
069            @Override
070            public double value(double x, double ... p) {
071                double v = Double.POSITIVE_INFINITY;
072                try {
073                    v = super.value(x, p);
074                } catch (NotStrictlyPositiveException e) { // NOPMD
075                    // Do nothing.
076                }
077                return v;
078            }
079
080            /** {@inheritDoc} */
081            @Override
082            public double[] gradient(double x, double ... p) {
083                double[] v = { Double.POSITIVE_INFINITY,
084                               Double.POSITIVE_INFINITY,
085                               Double.POSITIVE_INFINITY };
086                try {
087                    v = super.gradient(x, p);
088                } catch (NotStrictlyPositiveException e) { // NOPMD
089                    // Do nothing.
090                }
091                return v;
092            }
093        };
094
095    /**
096     * Constructor used by the factory methods.
097     *
098     * @param initialGuess Initial guess. If set to {@code null}, the initial guess
099     * will be estimated using the {@link ParameterGuesser}.
100     * @param maxIter Maximum number of iterations of the optimization algorithm.
101     */
102    private GaussianCurveFitter(double[] initialGuess,
103                                int maxIter) {
104        super(FUNCTION, initialGuess, new ParameterGuesser(), maxIter);
105    }
106
107    /**
108     * Creates a default curve fitter.
109     * The initial guess for the parameters will be {@link ParameterGuesser}
110     * computed automatically, and the maximum number of iterations of the
111     * optimization algorithm is set to {@link Integer#MAX_VALUE}.
112     *
113     * @return a curve fitter.
114     *
115     * @see #withStartPoint(double[])
116     * @see #withMaxIterations(int)
117     */
118    public static GaussianCurveFitter create() {
119        return new GaussianCurveFitter(null, Integer.MAX_VALUE);
120    }
121
122    /**
123     * Guesses the parameters {@code norm}, {@code mean}, and {@code sigma}
124     * of a {@link org.apache.commons.math4.legacy.analysis.function.Gaussian.Parametric}
125     * based on the specified observed points.
126     */
127    public static class ParameterGuesser extends SimpleCurveFitter.ParameterGuesser {
128        /**
129         * {@inheritDoc}
130         *
131         * @return the guessed parameters, in the following order:
132         * <ul>
133         *  <li>Normalization factor</li>
134         *  <li>Mean</li>
135         *  <li>Standard deviation</li>
136         * </ul>
137         * @throws NullArgumentException if {@code observations} is
138         * {@code null}.
139         * @throws NumberIsTooSmallException if there are less than 3
140         * observations.
141         */
142        @Override
143        public double[] guess(Collection<WeightedObservedPoint> observations) {
144            if (observations == null) {
145                throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
146            }
147            if (observations.size() < 3) {
148                throw new NumberIsTooSmallException(observations.size(), 3, true);
149            }
150
151            final List<WeightedObservedPoint> sorted = sortObservations(observations);
152            return basicGuess(sorted.toArray(new WeightedObservedPoint[0]));
153        }
154
155        /**
156         * Guesses the parameters based on the specified observed points.
157         *
158         * @param points Observed points, sorted.
159         * @return the guessed parameters (normalization factor, mean and
160         * sigma).
161         */
162        private double[] basicGuess(WeightedObservedPoint[] points) {
163            final int maxYIdx = findMaxY(points);
164            final double n = points[maxYIdx].getY();
165
166            double fwhmApprox;
167            try {
168                final double halfY = 0.5 * n;
169                final double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY);
170                final double fwhmX2 = interpolateXAtY(points, maxYIdx, 1, halfY);
171                fwhmApprox = fwhmX2 - fwhmX1;
172            } catch (OutOfRangeException e) {
173                // TODO: Exceptions should not be used for flow control.
174                fwhmApprox = points[points.length - 1].getX() - points[0].getX();
175            }
176            final double s = fwhmApprox / (2 * JdkMath.sqrt(2 * JdkMath.log(2)));
177
178            return new double[] { n, points[maxYIdx].getX(), s };
179        }
180    }
181}