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
018package org.apache.commons.math3.analysis.function;
019
020import java.util.Arrays;
021
022import org.apache.commons.math3.analysis.FunctionUtils;
023import org.apache.commons.math3.analysis.UnivariateFunction;
024import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction;
025import org.apache.commons.math3.analysis.ParametricUnivariateFunction;
026import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
027import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction;
028import org.apache.commons.math3.exception.NotStrictlyPositiveException;
029import org.apache.commons.math3.exception.NullArgumentException;
030import org.apache.commons.math3.exception.DimensionMismatchException;
031import org.apache.commons.math3.util.FastMath;
032import org.apache.commons.math3.util.Precision;
033
034/**
035 * <a href="http://en.wikipedia.org/wiki/Gaussian_function">
036 *  Gaussian</a> function.
037 *
038 * @since 3.0
039 */
040public class Gaussian implements UnivariateDifferentiableFunction, DifferentiableUnivariateFunction {
041    /** Mean. */
042    private final double mean;
043    /** Inverse of the standard deviation. */
044    private final double is;
045    /** Inverse of twice the square of the standard deviation. */
046    private final double i2s2;
047    /** Normalization factor. */
048    private final double norm;
049
050    /**
051     * Gaussian with given normalization factor, mean and standard deviation.
052     *
053     * @param norm Normalization factor.
054     * @param mean Mean.
055     * @param sigma Standard deviation.
056     * @throws NotStrictlyPositiveException if {@code sigma <= 0}.
057     */
058    public Gaussian(double norm,
059                    double mean,
060                    double sigma)
061        throws NotStrictlyPositiveException {
062        if (sigma <= 0) {
063            throw new NotStrictlyPositiveException(sigma);
064        }
065
066        this.norm = norm;
067        this.mean = mean;
068        this.is   = 1 / sigma;
069        this.i2s2 = 0.5 * is * is;
070    }
071
072    /**
073     * Normalized gaussian with given mean and standard deviation.
074     *
075     * @param mean Mean.
076     * @param sigma Standard deviation.
077     * @throws NotStrictlyPositiveException if {@code sigma <= 0}.
078     */
079    public Gaussian(double mean,
080                    double sigma)
081        throws NotStrictlyPositiveException {
082        this(1 / (sigma * FastMath.sqrt(2 * Math.PI)), mean, sigma);
083    }
084
085    /**
086     * Normalized gaussian with zero mean and unit standard deviation.
087     */
088    public Gaussian() {
089        this(0, 1);
090    }
091
092    /** {@inheritDoc} */
093    public double value(double x) {
094        return value(x - mean, norm, i2s2);
095    }
096
097    /** {@inheritDoc}
098     * @deprecated as of 3.1, replaced by {@link #value(DerivativeStructure)}
099     */
100    @Deprecated
101    public UnivariateFunction derivative() {
102        return FunctionUtils.toDifferentiableUnivariateFunction(this).derivative();
103    }
104
105    /**
106     * Parametric function where the input array contains the parameters of
107     * the Gaussian, ordered as follows:
108     * <ul>
109     *  <li>Norm</li>
110     *  <li>Mean</li>
111     *  <li>Standard deviation</li>
112     * </ul>
113     */
114    public static class Parametric implements ParametricUnivariateFunction {
115        /**
116         * Computes the value of the Gaussian at {@code x}.
117         *
118         * @param x Value for which the function must be computed.
119         * @param param Values of norm, mean and standard deviation.
120         * @return the value of the function.
121         * @throws NullArgumentException if {@code param} is {@code null}.
122         * @throws DimensionMismatchException if the size of {@code param} is
123         * not 3.
124         * @throws NotStrictlyPositiveException if {@code param[2]} is negative.
125         */
126        public double value(double x, double ... param)
127            throws NullArgumentException,
128                   DimensionMismatchException,
129                   NotStrictlyPositiveException {
130            validateParameters(param);
131
132            final double diff = x - param[1];
133            final double i2s2 = 1 / (2 * param[2] * param[2]);
134            return Gaussian.value(diff, param[0], i2s2);
135        }
136
137        /**
138         * Computes the value of the gradient at {@code x}.
139         * The components of the gradient vector are the partial
140         * derivatives of the function with respect to each of the
141         * <em>parameters</em> (norm, mean and standard deviation).
142         *
143         * @param x Value at which the gradient must be computed.
144         * @param param Values of norm, mean and standard deviation.
145         * @return the gradient vector at {@code x}.
146         * @throws NullArgumentException if {@code param} is {@code null}.
147         * @throws DimensionMismatchException if the size of {@code param} is
148         * not 3.
149         * @throws NotStrictlyPositiveException if {@code param[2]} is negative.
150         */
151        public double[] gradient(double x, double ... param)
152            throws NullArgumentException,
153                   DimensionMismatchException,
154                   NotStrictlyPositiveException {
155            validateParameters(param);
156
157            final double norm = param[0];
158            final double diff = x - param[1];
159            final double sigma = param[2];
160            final double i2s2 = 1 / (2 * sigma * sigma);
161
162            final double n = Gaussian.value(diff, 1, i2s2);
163            final double m = norm * n * 2 * i2s2 * diff;
164            final double s = m * diff / sigma;
165
166            return new double[] { n, m, s };
167        }
168
169        /**
170         * Validates parameters to ensure they are appropriate for the evaluation of
171         * the {@link #value(double,double[])} and {@link #gradient(double,double[])}
172         * methods.
173         *
174         * @param param Values of norm, mean and standard deviation.
175         * @throws NullArgumentException if {@code param} is {@code null}.
176         * @throws DimensionMismatchException if the size of {@code param} is
177         * not 3.
178         * @throws NotStrictlyPositiveException if {@code param[2]} is negative.
179         */
180        private void validateParameters(double[] param)
181            throws NullArgumentException,
182                   DimensionMismatchException,
183                   NotStrictlyPositiveException {
184            if (param == null) {
185                throw new NullArgumentException();
186            }
187            if (param.length != 3) {
188                throw new DimensionMismatchException(param.length, 3);
189            }
190            if (param[2] <= 0) {
191                throw new NotStrictlyPositiveException(param[2]);
192            }
193        }
194    }
195
196    /**
197     * @param xMinusMean {@code x - mean}.
198     * @param norm Normalization factor.
199     * @param i2s2 Inverse of twice the square of the standard deviation.
200     * @return the value of the Gaussian at {@code x}.
201     */
202    private static double value(double xMinusMean,
203                                double norm,
204                                double i2s2) {
205        return norm * FastMath.exp(-xMinusMean * xMinusMean * i2s2);
206    }
207
208    /** {@inheritDoc}
209     * @since 3.1
210     */
211    public DerivativeStructure value(final DerivativeStructure t)
212        throws DimensionMismatchException {
213
214        final double u = is * (t.getValue() - mean);
215        double[] f = new double[t.getOrder() + 1];
216
217        // the nth order derivative of the Gaussian has the form:
218        // dn(g(x)/dxn = (norm / s^n) P_n(u) exp(-u^2/2) with u=(x-m)/s
219        // where P_n(u) is a degree n polynomial with same parity as n
220        // P_0(u) = 1, P_1(u) = -u, P_2(u) = u^2 - 1, P_3(u) = -u^3 + 3 u...
221        // the general recurrence relation for P_n is:
222        // P_n(u) = P_(n-1)'(u) - u P_(n-1)(u)
223        // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
224        final double[] p = new double[f.length];
225        p[0] = 1;
226        final double u2 = u * u;
227        double coeff = norm * FastMath.exp(-0.5 * u2);
228        if (coeff <= Precision.SAFE_MIN) {
229            Arrays.fill(f, 0.0);
230        } else {
231            f[0] = coeff;
232            for (int n = 1; n < f.length; ++n) {
233
234                // update and evaluate polynomial P_n(x)
235                double v = 0;
236                p[n] = -p[n - 1];
237                for (int k = n; k >= 0; k -= 2) {
238                    v = v * u2 + p[k];
239                    if (k > 2) {
240                        p[k - 2] = (k - 1) * p[k - 1] - p[k - 3];
241                    } else if (k == 2) {
242                        p[0] = p[1];
243                    }
244                }
245                if ((n & 0x1) == 1) {
246                    v *= u;
247                }
248
249                coeff *= is;
250                f[n] = coeff * v;
251
252            }
253        }
254
255        return t.compose(f);
256
257    }
258
259}