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 }