Home | History | Annotate | Download | only in fitting 
      1 /*
      2  * Licensed to the Apache Software Foundation (ASF) under one or more
      3  * contributor license agreements.  See the NOTICE file distributed with
      4  * this work for additional information regarding copyright ownership.
      5  * The ASF licenses this file to You under the Apache License, Version 2.0
      6  * (the "License"); you may not use this file except in compliance with
      7  * the License.  You may obtain a copy of the License at
      8  *
      9  *      http://www.apache.org/licenses/LICENSE-2.0
     10  *
     11  * Unless required by applicable law or agreed to in writing, software
     12  * distributed under the License is distributed on an "AS IS" BASIS,
     13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
     14  * See the License for the specific language governing permissions and
     15  * limitations under the License.
     16  */
     17 
     18 package org.apache.commons.math.optimization.fitting;
     19 
     20 import java.io.Serializable;
     21 
     22 import org.apache.commons.math.exception.DimensionMismatchException;
     23 import org.apache.commons.math.exception.util.LocalizedFormats;
     24 import org.apache.commons.math.exception.ZeroException;
     25 import org.apache.commons.math.exception.NullArgumentException;
     26 import org.apache.commons.math.optimization.fitting.ParametricRealFunction;
     27 
     28 /**
     29  * A Gaussian function.  Specifically:
     30  * <p>
     31  * <tt>f(x) = a + b*exp(-((x - c)^2 / (2*d^2)))</tt>
     32  * <p>
     33  * The parameters have the following meaning:
     34  * <ul>
     35  * <li><tt>a</tt> is a constant offset that shifts <tt>f(x)</tt> up or down
     36  * <li><tt>b</tt> is the height of the peak
     37  * <li><tt>c</tt> is the position of the center of the peak
     38  * <li><tt>d</tt> is related to the FWHM by <tt>FWHM = 2*sqrt(2*ln(2))*d</tt>
     39  * </ul>
     40  * Notation key:
     41  * <ul>
     42  * <li><tt>x^n</tt>: <tt>x</tt> raised to the power of <tt>n</tt>
     43  * <li><tt>exp(x)</tt>: <i>e</i><tt>^x</tt>
     44  * <li><tt>sqrt(x)</tt>: the square root of <tt>x</tt>
     45  * <li><tt>ln(x)</tt>: the natural logarithm of <tt>x</tt>
     46  * </ul>
     47  * References:
     48  * <ul>
     49  * <li><a href="http://en.wikipedia.org/wiki/Gaussian_function">Wikipedia:
     50  *   Gaussian function</a>
     51  * </ul>
     52  *
     53  * @since 2.2
     54  * @version $Revision: 1037327 $ $Date: 2010-11-20 21:57:37 +0100 (sam. 20 nov. 2010) $
     55  */
     56 public class ParametricGaussianFunction implements ParametricRealFunction, Serializable {
     57 
     58     /** Serializable version Id. */
     59     private static final long serialVersionUID = -3875578602503903233L;
     60 
     61     /**
     62      * Constructs an instance.
     63      */
     64     public ParametricGaussianFunction() {
     65     }
     66 
     67     /**
     68      * Computes value of function <tt>f(x)</tt> for the specified <tt>x</tt> and
     69      * parameters <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and <tt>d</tt>.
     70      *
     71      * @param x <tt>x</tt> value
     72      * @param parameters values of <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and
     73      *        <tt>d</tt>
     74      *
     75      * @return value of <tt>f(x)</tt> evaluated at <tt>x</tt> with the specified
     76      *         parameters
     77      *
     78      * @throws IllegalArgumentException if <code>parameters</code> is invalid as
     79      *         determined by {@link #validateParameters(double[])}
     80      * @throws ZeroException if <code>parameters</code> values are
     81      *         invalid as determined by {@link #validateParameters(double[])}
     82      */
     83     public double value(double x, double[] parameters) throws ZeroException {
     84         validateParameters(parameters);
     85         final double a = parameters[0];
     86         final double b = parameters[1];
     87         final double c = parameters[2];
     88         final double d = parameters[3];
     89         final double xMc = x - c;
     90         return a + b * Math.exp(-xMc * xMc / (2.0 * (d * d)));
     91     }
     92 
     93     /**
     94      * Computes the gradient vector for a four variable version of the function
     95      * where the parameters, <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and <tt>d</tt>,
     96      * are considered the variables, not <tt>x</tt>.  That is, instead of
     97      * computing the gradient vector for the function <tt>f(x)</tt> (which would
     98      * just be the derivative of <tt>f(x)</tt> with respect to <tt>x</tt> since
     99      * 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 }
    166