    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.io.Serializable;
    021
    022    import org.apache.commons.math.exception.DimensionMismatchException;
    023    import org.apache.commons.math.exception.util.LocalizedFormats;
    024    import org.apache.commons.math.exception.ZeroException;
    025    import org.apache.commons.math.exception.NullArgumentException;
    026    import org.apache.commons.math.optimization.fitting.ParametricRealFunction;
    027
    028    /**
    029     * A Gaussian function.  Specifically:
    030     * <p>
    031     * <tt>f(x) = a + b*exp(-((x - c)^2 / (2*d^2)))</tt>
    032     * <p>
    033     * The parameters have the following meaning:
    034     * <ul>
    035     * <li><tt>a</tt> is a constant offset that shifts <tt>f(x)</tt> up or down
    036     * <li><tt>b</tt> is the height of the peak
    037     * <li><tt>c</tt> is the position of the center of the peak
    038     * <li><tt>d</tt> is related to the FWHM by <tt>FWHM = 2*sqrt(2*ln(2))*d</tt>
    039     * </ul>
    040     * Notation key:
    041     * <ul>
    042     * <li><tt>x^n</tt>: <tt>x</tt> raised to the power of <tt>n</tt>
    043     * <li><tt>exp(x)</tt>: <i>e</i><tt>^x</tt>
    044     * <li><tt>sqrt(x)</tt>: the square root of <tt>x</tt>
    045     * <li><tt>ln(x)</tt>: the natural logarithm of <tt>x</tt>
    046     * </ul>
    047     * References:
    048     * <ul>
    049     * <li><a href="http://en.wikipedia.org/wiki/Gaussian_function">Wikipedia:
    050     *   Gaussian function</a>
    051     * </ul>
    052     *
    053     * @since 2.2
    054     * @version $Revision: 1037327 $ $Date: 2010-11-20 21:57:37 +0100 (sam. 20 nov. 2010) $
    055     */
    056    public class ParametricGaussianFunction implements ParametricRealFunction, Serializable {
    057
    058        /** Serializable version Id. */
    059        private static final long serialVersionUID = -3875578602503903233L;
    060
    061        /**
    062         * Constructs an instance.
    063         */
    064        public ParametricGaussianFunction() {
    065        }
    066
    067        /**
    068         * Computes value of function <tt>f(x)</tt> for the specified <tt>x</tt> and
    069         * parameters <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and <tt>d</tt>.
    070         *
    071         * @param x <tt>x</tt> value
    072         * @param parameters values of <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and
    073         *        <tt>d</tt>
    074         *
    075         * @return value of <tt>f(x)</tt> evaluated at <tt>x</tt> with the specified
    076         *         parameters
    077         *
    078         * @throws IllegalArgumentException if <code>parameters</code> is invalid as
    079         *         determined by {@link #validateParameters(double[])}
    080         * @throws ZeroException if <code>parameters</code> values are
    081         *         invalid as determined by {@link #validateParameters(double[])}
    082         */
    083        public double value(double x, double[] parameters) throws ZeroException {
    084            validateParameters(parameters);
    085            final double a = parameters[0];
    086            final double b = parameters[1];
    087            final double c = parameters[2];
    088            final double d = parameters[3];
    089            final double xMc = x - c;
    090            return a + b * Math.exp(-xMc * xMc / (2.0 * (d * d)));
    091        }
    092
    093        /**
    094         * Computes the gradient vector for a four variable version of the function
    095         * where the parameters, <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and <tt>d</tt>,
    096         * are considered the variables, not <tt>x</tt>.  That is, instead of
    097         * computing the gradient vector for the function <tt>f(x)</tt> (which would
    098         * just be the derivative of <tt>f(x)</tt> with respect to <tt>x</tt> since
    099         * it's a one-dimensional function), computes the gradient vector for the
    100         * function <tt>f(a, b, c, d) = a + b*exp(-((x - c)^2 / (2*d^2)))</tt>
    101         * treating the specified <tt>x</tt> as a constant.
    102         * <p>
    103         * The components of the computed gradient vector are the partial
    104         * derivatives of <tt>f(a, b, c, d)</tt> with respect to each variable.
    105         * That is, the partial derivative of <tt>f(a, b, c, d)</tt> with respect to
    106         * <tt>a</tt>, the partial derivative of <tt>f(a, b, c, d)</tt> with respect
    107         * to <tt>b</tt>, the partial derivative of <tt>f(a, b, c, d)</tt> with
    108         * respect to <tt>c</tt>, and the partial derivative of <tt>f(a, b, c,
    109         * d)</tt> with respect to <tt>d</tt>.
    110         *
    111         * @param x <tt>x</tt> value to be used as constant in <tt>f(a, b, c,
    112         *        d)</tt>
    113         * @param parameters values of <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and
    114         *        <tt>d</tt> for computation of gradient vector of <tt>f(a, b, c,
    115         *        d)</tt>
    116         *
    117         * @return gradient vector of <tt>f(a, b, c, d)</tt>
    118         *
    119         * @throws IllegalArgumentException if <code>parameters</code> is invalid as
    120         *         determined by {@link #validateParameters(double[])}
    121         * @throws ZeroException if <code>parameters</code> values are
    122         *         invalid as determined by {@link #validateParameters(double[])}
    123         */
    124        public double[] gradient(double x, double[] parameters) throws ZeroException {
    125
    126            validateParameters(parameters);
    127            final double b = parameters[1];
    128            final double c = parameters[2];
    129            final double d = parameters[3];
    130
    131            final double xMc  = x - c;
    132            final double d2   = d * d;
    133            final double exp  = Math.exp(-xMc * xMc / (2 * d2));
    134            final double f    = b * exp * xMc / d2;
    135
    136            return new double[] { 1.0, exp, f, f * xMc / d };
    137
    138        }
    139
    140        /**
    141         * Validates parameters to ensure they are appropriate for the evaluation of
    142         * the <code>value</code> and <code>gradient</code> methods.
    143         *
    144         * @param parameters values of <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and
    145         *        <tt>d</tt>
    146         *
    147         * @throws IllegalArgumentException if <code>parameters</code> is
    148         *         <code>null</code> or if <code>parameters</code> does not have
    149         *         length == 4
    150         * @throws ZeroException if <code>parameters[3]</code>
    151         *         (<tt>d</tt>) is 0
    152         */
    153        private void validateParameters(double[] parameters) throws ZeroException {
    154            if (parameters == null) {
    155                throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
    156            }
    157            if (parameters.length != 4) {
    158                throw new DimensionMismatchException(4, parameters.length);
    159            }
    160            if (parameters[3] == 0.0) {
    161                throw new ZeroException();
    162            }
    163        }
    164
    165    }