Added support for Gaussian curve fitting
JIRA: MATH-400 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@980938 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
7adf364d60
commit
d9e60241d2
|
|
@ -149,6 +149,11 @@
|
|||
<module name="PackageHtml"/>
|
||||
|
||||
<!-- Setup special comments to suppress specific checks from source files -->
|
||||
<module name="SuppressionCommentFilter">
|
||||
<property name="offCommentFormat" value="CHECKSTYLE\: stop JavadocVariable"/>
|
||||
<property name="onCommentFormat" value="CHECKSTYLE\: resume JavadocVariable"/>
|
||||
<property name="checkFormat" value="JavadocVariable"/>
|
||||
</module>
|
||||
<module name="SuppressionCommentFilter">
|
||||
<property name="offCommentFormat" value="CHECKSTYLE\: stop JavadocMethodCheck"/>
|
||||
<property name="onCommentFormat" value="CHECKSTYLE\: resume JavadocMethodCheck"/>
|
||||
|
|
|
|||
3
pom.xml
3
pom.xml
|
|
@ -181,6 +181,9 @@
|
|||
<contributor>
|
||||
<name>Benjamin McCann</name>
|
||||
</contributor>
|
||||
<contributor>
|
||||
<name>J. Lewis Muir</name>
|
||||
</contributor>
|
||||
<contributor>
|
||||
<name>Fredrik Norin</name>
|
||||
</contributor>
|
||||
|
|
|
|||
|
|
@ -0,0 +1,49 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
package org.apache.commons.math.exception;
|
||||
|
||||
import org.apache.commons.math.util.Localizable;
|
||||
import org.apache.commons.math.util.LocalizedFormats;
|
||||
|
||||
/**
|
||||
* Exception to be thrown when zero is provided where it is not allowed.
|
||||
*
|
||||
* @since 2.2
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class ZeroNotAllowedException extends MathIllegalNumberException {
|
||||
|
||||
/** Serializable version identifier */
|
||||
private static final long serialVersionUID = -1960874856936000015L;
|
||||
|
||||
/**
|
||||
* Construct the exception.
|
||||
*/
|
||||
public ZeroNotAllowedException() {
|
||||
this(null);
|
||||
}
|
||||
|
||||
/**
|
||||
* Construct the exception with a specific context.
|
||||
*
|
||||
* @param specific Specific contexte pattern .
|
||||
*/
|
||||
public ZeroNotAllowedException(Localizable specific) {
|
||||
super(specific, LocalizedFormats.ZERO_NOT_ALLOWED, 0);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -0,0 +1,105 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
|
||||
package org.apache.commons.math.optimization.fitting;
|
||||
|
||||
import java.io.Serializable;
|
||||
|
||||
import org.apache.commons.math.FunctionEvaluationException;
|
||||
import org.apache.commons.math.MathRuntimeException;
|
||||
import org.apache.commons.math.analysis.UnivariateRealFunction;
|
||||
import org.apache.commons.math.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math.exception.ZeroNotAllowedException;
|
||||
import org.apache.commons.math.util.LocalizedFormats;
|
||||
|
||||
/**
|
||||
* The derivative of {@link GaussianFunction}. Specifically:
|
||||
* <p>
|
||||
* <tt>f'(x) = (-b / (d^2)) * (x - c) * exp(-((x - c)^2) / (2*(d^2)))</tt>
|
||||
* <p>
|
||||
* Notation key:
|
||||
* <ul>
|
||||
* <li><tt>x^n</tt>: <tt>x</tt> raised to the power of <tt>n</tt>
|
||||
* <li><tt>exp(x)</tt>: <i>e</i><tt>^x</tt>
|
||||
* </ul>
|
||||
*
|
||||
* @since 2.2
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class GaussianDerivativeFunction implements UnivariateRealFunction, Serializable {
|
||||
|
||||
/** Serializable version identifier. */
|
||||
private static final long serialVersionUID = -6500229089670174766L;
|
||||
|
||||
/** Parameter b of this function. */
|
||||
private final double b;
|
||||
|
||||
/** Parameter c of this function. */
|
||||
private final double c;
|
||||
|
||||
/** Square of the parameter d of this function. */
|
||||
private final double d2;
|
||||
|
||||
/**
|
||||
* Constructs an instance with the specified parameters.
|
||||
*
|
||||
* @param b <tt>b</tt> parameter value
|
||||
* @param c <tt>c</tt> parameter value
|
||||
* @param d <tt>d</tt> parameter value
|
||||
*
|
||||
* @throws IllegalArgumentException if <code>d</code> is 0
|
||||
*/
|
||||
public GaussianDerivativeFunction(double b, double c, double d) {
|
||||
if (d == 0.0) {
|
||||
throw new ZeroNotAllowedException();
|
||||
}
|
||||
this.b = b;
|
||||
this.c = c;
|
||||
this.d2 = d * d;
|
||||
}
|
||||
|
||||
/**
|
||||
* Constructs an instance with the specified parameters.
|
||||
*
|
||||
* @param parameters <tt>b</tt>, <tt>c</tt>, and <tt>d</tt> parameter values
|
||||
*
|
||||
* @throws IllegalArgumentException if <code>parameters</code> is null,
|
||||
* <code>parameters</code> length is not 3, or if
|
||||
* <code>parameters[2]</code> is 0
|
||||
*/
|
||||
public GaussianDerivativeFunction(double[] parameters) {
|
||||
if (parameters == null) {
|
||||
throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.NULL_INPUT_ARRAY);
|
||||
}
|
||||
if (parameters.length != 3) {
|
||||
throw new DimensionMismatchException(3, parameters.length);
|
||||
}
|
||||
if (parameters[2] == 0.0) {
|
||||
throw new ZeroNotAllowedException();
|
||||
}
|
||||
this.b = parameters[0];
|
||||
this.c = parameters[1];
|
||||
this.d2 = parameters[2] * parameters[2];
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double value(double x) throws FunctionEvaluationException {
|
||||
final double xMc = x - c;
|
||||
return (-b / d2) * xMc * Math.exp(-(xMc * xMc) / (2.0 * d2));
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -0,0 +1,121 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
|
||||
package org.apache.commons.math.optimization.fitting;
|
||||
|
||||
import org.apache.commons.math.FunctionEvaluationException;
|
||||
import org.apache.commons.math.optimization.
|
||||
DifferentiableMultivariateVectorialOptimizer;
|
||||
import org.apache.commons.math.optimization.OptimizationException;
|
||||
import org.apache.commons.math.optimization.fitting.CurveFitter;
|
||||
import org.apache.commons.math.optimization.fitting.WeightedObservedPoint;
|
||||
|
||||
/**
|
||||
* Fits points to a Gaussian function (that is, a {@link GaussianFunction}).
|
||||
* <p>
|
||||
* Usage example:
|
||||
* <pre>
|
||||
* GaussianFitter fitter = new GaussianFitter(
|
||||
* new LevenbergMarquardtOptimizer());
|
||||
* fitter.addObservedPoint(4.0254623, 531026.0);
|
||||
* fitter.addObservedPoint(4.03128248, 984167.0);
|
||||
* fitter.addObservedPoint(4.03839603, 1887233.0);
|
||||
* fitter.addObservedPoint(4.04421621, 2687152.0);
|
||||
* fitter.addObservedPoint(4.05132976, 3461228.0);
|
||||
* fitter.addObservedPoint(4.05326982, 3580526.0);
|
||||
* fitter.addObservedPoint(4.05779662, 3439750.0);
|
||||
* fitter.addObservedPoint(4.0636168, 2877648.0);
|
||||
* fitter.addObservedPoint(4.06943698, 2175960.0);
|
||||
* fitter.addObservedPoint(4.07525716, 1447024.0);
|
||||
* fitter.addObservedPoint(4.08237071, 717104.0);
|
||||
* fitter.addObservedPoint(4.08366408, 620014.0);
|
||||
* GaussianFunction fitFunction = fitter.fit();
|
||||
* </pre>
|
||||
*
|
||||
* @see ParametricGaussianFunction
|
||||
* @since 2.2
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class GaussianFitter {
|
||||
|
||||
/** Fitter used for fitting. */
|
||||
private final CurveFitter fitter;
|
||||
|
||||
/**
|
||||
* Constructs an instance using the specified optimizer.
|
||||
*
|
||||
* @param optimizer optimizer to use for the fitting
|
||||
*/
|
||||
public GaussianFitter(DifferentiableMultivariateVectorialOptimizer optimizer) {
|
||||
fitter = new CurveFitter(optimizer);
|
||||
}
|
||||
|
||||
/**
|
||||
* Adds point (<code>x</code>, <code>y</code>) to list of observed points
|
||||
* with a weight of 1.0.
|
||||
*
|
||||
* @param x <tt>x</tt> point value
|
||||
* @param y <tt>y</tt> point value
|
||||
*/
|
||||
public void addObservedPoint(double x, double y) {
|
||||
addObservedPoint(1.0, x, y);
|
||||
}
|
||||
|
||||
/**
|
||||
* Adds point (<code>x</code>, <code>y</code>) to list of observed points
|
||||
* with a weight of <code>weight</code>.
|
||||
*
|
||||
* @param weight weight assigned to point
|
||||
* @param x <tt>x</tt> point value
|
||||
* @param y <tt>y</tt> point value
|
||||
*/
|
||||
public void addObservedPoint(double weight, double x, double y) {
|
||||
fitter.addObservedPoint(weight, x, y);
|
||||
}
|
||||
|
||||
/**
|
||||
* Fits Gaussian function to the observed points.
|
||||
*
|
||||
* @return Gaussian function best fitting the observed points
|
||||
*
|
||||
* @throws FunctionEvaluationException if <code>CurveFitter.fit</code>
|
||||
* throws it
|
||||
* @throws OptimizationException if <code>CurveFitter.fit</code> throws it
|
||||
* @throws IllegalArgumentException if <code>CurveFitter.fit</code> throws
|
||||
* it
|
||||
*
|
||||
* @see CurveFitter
|
||||
*/
|
||||
public GaussianFunction fit()
|
||||
throws FunctionEvaluationException, OptimizationException {
|
||||
return new GaussianFunction(fitter.fit(new ParametricGaussianFunction(),
|
||||
createParametersGuesser(fitter.getObservations()).guess()));
|
||||
}
|
||||
|
||||
/**
|
||||
* Factory method to create a <code>GaussianParametersGuesser</code>
|
||||
* instance initialized with the specified observations.
|
||||
*
|
||||
* @param observations points used to initialize the created
|
||||
* <code>GaussianParametersGuesser</code> instance
|
||||
*
|
||||
* @return new <code>GaussianParametersGuesser</code> instance
|
||||
*/
|
||||
protected GaussianParametersGuesser createParametersGuesser(WeightedObservedPoint[] observations) {
|
||||
return new GaussianParametersGuesser(observations);
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,160 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
|
||||
package org.apache.commons.math.optimization.fitting;
|
||||
|
||||
import java.io.Serializable;
|
||||
|
||||
import org.apache.commons.math.FunctionEvaluationException;
|
||||
import org.apache.commons.math.MathRuntimeException;
|
||||
import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
|
||||
import org.apache.commons.math.analysis.UnivariateRealFunction;
|
||||
import org.apache.commons.math.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math.exception.ZeroNotAllowedException;
|
||||
import org.apache.commons.math.util.LocalizedFormats;
|
||||
|
||||
/**
|
||||
* A Gaussian function. Specifically:
|
||||
* <p>
|
||||
* <tt>f(x) = a + b*exp(-((x - c)^2 / (2*d^2)))</tt>
|
||||
* <p>
|
||||
* Notation key:
|
||||
* <ul>
|
||||
* <li><tt>x^n</tt>: <tt>x</tt> raised to the power of <tt>n</tt>
|
||||
* <li><tt>exp(x)</tt>: <i>e</i><tt>^x</tt>
|
||||
* </ul>
|
||||
* References:
|
||||
* <ul>
|
||||
* <li><a href="http://en.wikipedia.org/wiki/Gaussian_function">Wikipedia:
|
||||
* Gaussian function</a>
|
||||
* </ul>
|
||||
*
|
||||
* @see GaussianDerivativeFunction
|
||||
* @see ParametricGaussianFunction
|
||||
* @since 2.2
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class GaussianFunction implements DifferentiableUnivariateRealFunction, Serializable {
|
||||
|
||||
/** Serializable version identifier. */
|
||||
private static final long serialVersionUID = -3195385616125629512L;
|
||||
|
||||
/** Parameter a of this function. */
|
||||
private final double a;
|
||||
|
||||
/** Parameter b of this function. */
|
||||
private final double b;
|
||||
|
||||
/** Parameter c of this function. */
|
||||
private final double c;
|
||||
|
||||
/** Parameter d of this function. */
|
||||
private final double d;
|
||||
|
||||
/**
|
||||
* Constructs an instance with the specified parameters.
|
||||
*
|
||||
* @param a <tt>a</tt> parameter value
|
||||
* @param b <tt>b</tt> parameter value
|
||||
* @param c <tt>c</tt> parameter value
|
||||
* @param d <tt>d</tt> parameter value
|
||||
*
|
||||
* @throws IllegalArgumentException if <code>d</code> is 0
|
||||
*/
|
||||
public GaussianFunction(double a, double b, double c, double d) {
|
||||
if (d == 0.0) {
|
||||
throw new ZeroNotAllowedException();
|
||||
}
|
||||
this.a = a;
|
||||
this.b = b;
|
||||
this.c = c;
|
||||
this.d = d;
|
||||
}
|
||||
|
||||
/**
|
||||
* Constructs an instance with the specified parameters.
|
||||
*
|
||||
* @param parameters <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and <tt>d</tt>
|
||||
* parameter values
|
||||
*
|
||||
* @throws IllegalArgumentException if <code>parameters</code> is null,
|
||||
* <code>parameters</code> length is not 4, or if
|
||||
* <code>parameters[3]</code> is 0
|
||||
*/
|
||||
public GaussianFunction(double[] parameters) {
|
||||
if (parameters == null) {
|
||||
throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.NULL_INPUT_ARRAY);
|
||||
}
|
||||
if (parameters.length != 4) {
|
||||
throw new DimensionMismatchException(4, parameters.length);
|
||||
}
|
||||
if (parameters[3] == 0.0) {
|
||||
throw new ZeroNotAllowedException();
|
||||
}
|
||||
this.a = parameters[0];
|
||||
this.b = parameters[1];
|
||||
this.c = parameters[2];
|
||||
this.d = parameters[3];
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public UnivariateRealFunction derivative() {
|
||||
return new GaussianDerivativeFunction(b, c, d);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double value(double x) throws FunctionEvaluationException {
|
||||
final double xMc = x - c;
|
||||
return a + b * Math.exp(-xMc * xMc / (2.0 * (d * d)));
|
||||
}
|
||||
|
||||
/**
|
||||
* Gets <tt>a</tt> parameter value.
|
||||
*
|
||||
* @return <tt>a</tt> parameter value
|
||||
*/
|
||||
public double getA() {
|
||||
return a;
|
||||
}
|
||||
|
||||
/**
|
||||
* Gets <tt>b</tt> parameter value.
|
||||
*
|
||||
* @return <tt>b</tt> parameter value
|
||||
*/
|
||||
public double getB() {
|
||||
return b;
|
||||
}
|
||||
|
||||
/**
|
||||
* Gets <tt>c</tt> parameter value.
|
||||
*
|
||||
* @return <tt>c</tt> parameter value
|
||||
*/
|
||||
public double getC() {
|
||||
return c;
|
||||
}
|
||||
|
||||
/**
|
||||
* Gets <tt>d</tt> parameter value.
|
||||
*
|
||||
* @return <tt>d</tt> parameter value
|
||||
*/
|
||||
public double getD() {
|
||||
return d;
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,271 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
|
||||
package org.apache.commons.math.optimization.fitting;
|
||||
|
||||
import java.util.Arrays;
|
||||
import java.util.Comparator;
|
||||
|
||||
import org.apache.commons.math.MathRuntimeException;
|
||||
import org.apache.commons.math.exception.NumberIsTooSmallException;
|
||||
import org.apache.commons.math.exception.OutOfRangeException;
|
||||
import org.apache.commons.math.exception.ZeroNotAllowedException;
|
||||
import org.apache.commons.math.util.LocalizedFormats;
|
||||
|
||||
/**
|
||||
* Guesses the parameters (<tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and <tt>d</tt>)
|
||||
* of a {@link ParametricGaussianFunction} based on the specified observed
|
||||
* points.
|
||||
*
|
||||
* @since 2.2
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class GaussianParametersGuesser {
|
||||
|
||||
/** Observed points. */
|
||||
private final WeightedObservedPoint[] observations;
|
||||
|
||||
/** Resulting guessed parameters. */
|
||||
private double[] parameters;
|
||||
|
||||
/**
|
||||
* Constructs instance with the specified observed points.
|
||||
*
|
||||
* @param observations observed points upon which should base guess
|
||||
*/
|
||||
public GaussianParametersGuesser(WeightedObservedPoint[] observations) {
|
||||
if (observations == null) {
|
||||
throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.NULL_INPUT_ARRAY);
|
||||
}
|
||||
if (observations.length < 3) {
|
||||
throw new NumberIsTooSmallException(observations.length, 3, true);
|
||||
}
|
||||
this.observations = observations.clone();
|
||||
}
|
||||
|
||||
/**
|
||||
* Guesses the parameters based on the observed points.
|
||||
*
|
||||
* @return guessed parameters array <code>{a, b, c, d}</code>
|
||||
*/
|
||||
public double[] guess() {
|
||||
if (parameters == null) {
|
||||
parameters = basicGuess(observations);
|
||||
}
|
||||
return parameters.clone();
|
||||
}
|
||||
|
||||
/**
|
||||
* Guesses the parameters based on the specified observed points.
|
||||
*
|
||||
* @param points observed points upon which should base guess
|
||||
*
|
||||
* @return guessed parameters array <code>{a, b, c, d}</code>
|
||||
*/
|
||||
private double[] basicGuess(WeightedObservedPoint[] points) {
|
||||
Arrays.sort(points, createWeightedObservedPointComparator());
|
||||
double[] params = new double[4];
|
||||
|
||||
int minYIdx = findMinY(points);
|
||||
params[0] = points[minYIdx].getY();
|
||||
|
||||
int maxYIdx = findMaxY(points);
|
||||
params[1] = points[maxYIdx].getY();
|
||||
params[2] = points[maxYIdx].getX();
|
||||
|
||||
double fwhmApprox;
|
||||
try {
|
||||
double halfY = params[0] + ((params[1] - params[0]) / 2.0);
|
||||
double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY);
|
||||
double fwhmX2 = interpolateXAtY(points, maxYIdx, +1, halfY);
|
||||
fwhmApprox = fwhmX2 - fwhmX1;
|
||||
} catch (OutOfRangeException e) {
|
||||
fwhmApprox = points[points.length - 1].getX() - points[0].getX();
|
||||
}
|
||||
params[3] = fwhmApprox / (2.0 * Math.sqrt(2.0 * Math.log(2.0)));
|
||||
|
||||
return params;
|
||||
}
|
||||
|
||||
/**
|
||||
* Finds index of point in specified points with the smallest Y.
|
||||
*
|
||||
* @param points points to search
|
||||
*
|
||||
* @return index in specified points array
|
||||
*/
|
||||
private int findMinY(WeightedObservedPoint[] points) {
|
||||
int minYIdx = 0;
|
||||
for (int i = 1; i < points.length; i++) {
|
||||
if (points[i].getY() < points[minYIdx].getY()) {
|
||||
minYIdx = i;
|
||||
}
|
||||
}
|
||||
return minYIdx;
|
||||
}
|
||||
|
||||
/**
|
||||
* Finds index of point in specified points with the largest Y.
|
||||
*
|
||||
* @param points points to search
|
||||
*
|
||||
* @return index in specified points array
|
||||
*/
|
||||
private int findMaxY(WeightedObservedPoint[] points) {
|
||||
int maxYIdx = 0;
|
||||
for (int i = 1; i < points.length; i++) {
|
||||
if (points[i].getY() > points[maxYIdx].getY()) {
|
||||
maxYIdx = i;
|
||||
}
|
||||
}
|
||||
return maxYIdx;
|
||||
}
|
||||
|
||||
/**
|
||||
* Interpolates using the specified points to determine X at the specified
|
||||
* Y.
|
||||
*
|
||||
* @param points points to use for interpolation
|
||||
* @param startIdx index within points from which to start search for
|
||||
* interpolation bounds points
|
||||
* @param idxStep index step for search for interpolation bounds points
|
||||
* @param y Y value for which X should be determined
|
||||
*
|
||||
* @return value of X at the specified Y
|
||||
*
|
||||
* @throws IllegalArgumentException if idxStep is 0
|
||||
* @throws OutOfRangeException if specified <code>y</code> is not within the
|
||||
* range of the specified <code>points</code>
|
||||
*/
|
||||
private double interpolateXAtY(WeightedObservedPoint[] points,
|
||||
int startIdx, int idxStep, double y) throws OutOfRangeException {
|
||||
if (idxStep == 0) {
|
||||
throw new ZeroNotAllowedException();
|
||||
}
|
||||
WeightedObservedPoint[] twoPoints = getInterpolationPointsForY(points, startIdx, idxStep, y);
|
||||
WeightedObservedPoint pointA = twoPoints[0];
|
||||
WeightedObservedPoint pointB = twoPoints[1];
|
||||
if (pointA.getY() == y) {
|
||||
return pointA.getX();
|
||||
}
|
||||
if (pointB.getY() == y) {
|
||||
return pointB.getX();
|
||||
}
|
||||
return pointA.getX() +
|
||||
(((y - pointA.getY()) * (pointB.getX() - pointA.getX())) / (pointB.getY() - pointA.getY()));
|
||||
}
|
||||
|
||||
/**
|
||||
* Gets the two bounding interpolation points from the specified points
|
||||
* suitable for determining X at the specified Y.
|
||||
*
|
||||
* @param points points to use for interpolation
|
||||
* @param startIdx index within points from which to start search for
|
||||
* interpolation bounds points
|
||||
* @param idxStep index step for search for interpolation bounds points
|
||||
* @param y Y value for which X should be determined
|
||||
*
|
||||
* @return array containing two points suitable for determining X at the
|
||||
* specified Y
|
||||
*
|
||||
* @throws IllegalArgumentException if idxStep is 0
|
||||
* @throws OutOfRangeException if specified <code>y</code> is not within the
|
||||
* range of the specified <code>points</code>
|
||||
*/
|
||||
private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points,
|
||||
int startIdx, int idxStep, double y)
|
||||
throws OutOfRangeException {
|
||||
if (idxStep == 0) {
|
||||
throw new ZeroNotAllowedException();
|
||||
}
|
||||
for (int i = startIdx;
|
||||
(idxStep < 0) ? (i + idxStep >= 0) : (i + idxStep < points.length);
|
||||
i += idxStep) {
|
||||
if (isBetween(y, points[i].getY(), points[i + idxStep].getY())) {
|
||||
return (idxStep < 0) ?
|
||||
new WeightedObservedPoint[] { points[i + idxStep], points[i] } :
|
||||
new WeightedObservedPoint[] { points[i], points[i + idxStep] };
|
||||
}
|
||||
}
|
||||
|
||||
double minY = Double.POSITIVE_INFINITY;
|
||||
double maxY = Double.NEGATIVE_INFINITY;
|
||||
for (final WeightedObservedPoint point : points) {
|
||||
minY = Math.min(minY, point.getY());
|
||||
maxY = Math.max(maxY, point.getY());
|
||||
}
|
||||
throw new OutOfRangeException(y, minY, maxY);
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* Determines whether a value is between two other values.
|
||||
*
|
||||
* @param value value to determine whether is between <code>boundary1</code>
|
||||
* and <code>boundary2</code>
|
||||
* @param boundary1 one end of the range
|
||||
* @param boundary2 other end of the range
|
||||
*
|
||||
* @return true if <code>value</code> is between <code>boundary1</code> and
|
||||
* <code>boundary2</code> (inclusive); false otherwise
|
||||
*/
|
||||
private boolean isBetween(double value, double boundary1, double boundary2) {
|
||||
return (value >= boundary1 && value <= boundary2) ||
|
||||
(value >= boundary2 && value <= boundary1);
|
||||
}
|
||||
|
||||
/**
|
||||
* Factory method creating <code>Comparator</code> for comparing
|
||||
* <code>WeightedObservedPoint</code> instances.
|
||||
*
|
||||
* @return new <code>Comparator</code> instance
|
||||
*/
|
||||
private Comparator<WeightedObservedPoint> createWeightedObservedPointComparator() {
|
||||
return new Comparator<WeightedObservedPoint>() {
|
||||
public int compare(WeightedObservedPoint p1, WeightedObservedPoint p2) {
|
||||
if (p1 == null && p2 == null) {
|
||||
return 0;
|
||||
}
|
||||
if (p1 == null) {
|
||||
return -1;
|
||||
}
|
||||
if (p2 == null) {
|
||||
return 1;
|
||||
}
|
||||
if (p1.getX() < p2.getX()) {
|
||||
return -1;
|
||||
}
|
||||
if (p1.getX() > p2.getX()) {
|
||||
return 1;
|
||||
}
|
||||
if (p1.getY() < p2.getY()) {
|
||||
return -1;
|
||||
}
|
||||
if (p1.getY() > p2.getY()) {
|
||||
return 1;
|
||||
}
|
||||
if (p1.getWeight() < p2.getWeight()) {
|
||||
return -1;
|
||||
}
|
||||
if (p1.getWeight() > p2.getWeight()) {
|
||||
return 1;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
};
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,166 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
|
||||
package org.apache.commons.math.optimization.fitting;
|
||||
|
||||
import java.io.Serializable;
|
||||
|
||||
import org.apache.commons.math.FunctionEvaluationException;
|
||||
import org.apache.commons.math.MathRuntimeException;
|
||||
import org.apache.commons.math.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math.exception.ZeroNotAllowedException;
|
||||
import org.apache.commons.math.optimization.fitting.ParametricRealFunction;
|
||||
import org.apache.commons.math.util.LocalizedFormats;
|
||||
|
||||
/**
|
||||
* A Gaussian function. Specifically:
|
||||
* <p>
|
||||
* <tt>f(x) = a + b*exp(-((x - c)^2 / (2*d^2)))</tt>
|
||||
* <p>
|
||||
* The parameters have the following meaning:
|
||||
* <ul>
|
||||
* <li><tt>a</tt> is a constant offset that shifts <tt>f(x)</tt> up or down
|
||||
* <li><tt>b</tt> is the height of the peak
|
||||
* <li><tt>c</tt> is the position of the center of the peak
|
||||
* <li><tt>d</tt> is related to the FWHM by <tt>FWHM = 2*sqrt(2*ln(2))*d</tt>
|
||||
* </ul>
|
||||
* Notation key:
|
||||
* <ul>
|
||||
* <li><tt>x^n</tt>: <tt>x</tt> raised to the power of <tt>n</tt>
|
||||
* <li><tt>exp(x)</tt>: <i>e</i><tt>^x</tt>
|
||||
* <li><tt>sqrt(x)</tt>: the square root of <tt>x</tt>
|
||||
* <li><tt>ln(x)</tt>: the natural logarithm of <tt>x</tt>
|
||||
* </ul>
|
||||
* References:
|
||||
* <ul>
|
||||
* <li><a href="http://en.wikipedia.org/wiki/Gaussian_function">Wikipedia:
|
||||
* Gaussian function</a>
|
||||
* </ul>
|
||||
*
|
||||
* @since 2.2
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class ParametricGaussianFunction implements ParametricRealFunction, Serializable {
|
||||
|
||||
/** Serializable version Id. */
|
||||
private static final long serialVersionUID = -3875578602503903233L;
|
||||
|
||||
/**
|
||||
* Constructs an instance.
|
||||
*/
|
||||
public ParametricGaussianFunction() {
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes value of function <tt>f(x)</tt> for the specified <tt>x</tt> and
|
||||
* parameters <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and <tt>d</tt>.
|
||||
*
|
||||
* @param x <tt>x</tt> value
|
||||
* @param parameters values of <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and
|
||||
* <tt>d</tt>
|
||||
*
|
||||
* @return value of <tt>f(x)</tt> evaluated at <tt>x</tt> with the specified
|
||||
* parameters
|
||||
*
|
||||
* @throws IllegalArgumentException if <code>parameters</code> is invalid as
|
||||
* determined by {@link #validateParameters(double[])}
|
||||
* @throws FunctionEvaluationException if <code>parameters</code> values are
|
||||
* invalid as determined by {@link #validateParameters(double[])}
|
||||
*/
|
||||
public double value(double x, double[] parameters) throws FunctionEvaluationException {
|
||||
validateParameters(parameters);
|
||||
final double a = parameters[0];
|
||||
final double b = parameters[1];
|
||||
final double c = parameters[2];
|
||||
final double d = parameters[3];
|
||||
final double xMc = x - c;
|
||||
return a + b * Math.exp(-xMc * xMc / (2.0 * (d * d)));
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes the gradient vector for a four variable version of the function
|
||||
* where the parameters, <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and <tt>d</tt>,
|
||||
* are considered the variables, not <tt>x</tt>. That is, instead of
|
||||
* computing the gradient vector for the function <tt>f(x)</tt> (which would
|
||||
* just be the derivative of <tt>f(x)</tt> with respect to <tt>x</tt> since
|
||||
* it's a one-dimensional function), computes the gradient vector for the
|
||||
* function <tt>f(a, b, c, d) = a + b*exp(-((x - c)^2 / (2*d^2)))</tt>
|
||||
* treating the specified <tt>x</tt> as a constant.
|
||||
* <p>
|
||||
* The components of the computed gradient vector are the partial
|
||||
* derivatives of <tt>f(a, b, c, d)</tt> with respect to each variable.
|
||||
* That is, the partial derivative of <tt>f(a, b, c, d)</tt> with respect to
|
||||
* <tt>a</tt>, the partial derivative of <tt>f(a, b, c, d)</tt> with respect
|
||||
* to <tt>b</tt>, the partial derivative of <tt>f(a, b, c, d)</tt> with
|
||||
* respect to <tt>c</tt>, and the partial derivative of <tt>f(a, b, c,
|
||||
* d)</tt> with respect to <tt>d</tt>.
|
||||
*
|
||||
* @param x <tt>x</tt> value to be used as constant in <tt>f(a, b, c,
|
||||
* d)</tt>
|
||||
* @param parameters values of <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and
|
||||
* <tt>d</tt> for computation of gradient vector of <tt>f(a, b, c,
|
||||
* d)</tt>
|
||||
*
|
||||
* @return gradient vector of <tt>f(a, b, c, d)</tt>
|
||||
*
|
||||
* @throws IllegalArgumentException if <code>parameters</code> is invalid as
|
||||
* determined by {@link #validateParameters(double[])}
|
||||
* @throws FunctionEvaluationException if <code>parameters</code> values are
|
||||
* invalid as determined by {@link #validateParameters(double[])}
|
||||
*/
|
||||
public double[] gradient(double x, double[] parameters) throws FunctionEvaluationException {
|
||||
|
||||
validateParameters(parameters);
|
||||
final double b = parameters[1];
|
||||
final double c = parameters[2];
|
||||
final double d = parameters[3];
|
||||
|
||||
final double xMc = x - c;
|
||||
final double d2 = d * d;
|
||||
final double exp = Math.exp(-xMc * xMc / (2 * d2));
|
||||
final double f = b * exp * xMc / d2;
|
||||
|
||||
return new double[] { 1.0, exp, f, f * xMc / d };
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* Validates parameters to ensure they are appropriate for the evaluation of
|
||||
* the <code>value</code> and <code>gradient</code> methods.
|
||||
*
|
||||
* @param parameters values of <tt>a</tt>, <tt>b</tt>, <tt>c</tt>, and
|
||||
* <tt>d</tt>
|
||||
*
|
||||
* @throws IllegalArgumentException if <code>parameters</code> is
|
||||
* <code>null</code> or if <code>parameters</code> does not have
|
||||
* length == 4
|
||||
* @throws FunctionEvaluationException if <code>parameters[3]</code>
|
||||
* (<tt>d</tt>) is 0
|
||||
*/
|
||||
private void validateParameters(double[] parameters) throws FunctionEvaluationException {
|
||||
if (parameters == null) {
|
||||
throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.NULL_INPUT_ARRAY);
|
||||
}
|
||||
if (parameters.length != 4) {
|
||||
throw new DimensionMismatchException(4, parameters.length);
|
||||
}
|
||||
if (parameters[3] == 0.0) {
|
||||
throw new ZeroNotAllowedException();
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -40,6 +40,7 @@ import java.util.ResourceBundle;
|
|||
public enum LocalizedFormats implements Localizable {
|
||||
|
||||
// CHECKSTYLE: stop MultipleVariableDeclarations
|
||||
// CHECKSTYLE: stop JavadocVariable
|
||||
|
||||
ARGUMENT_OUTSIDE_DOMAIN("Argument {0} outside domain [{1} ; {2}]"),
|
||||
ARRAY_SIZES_SHOULD_HAVE_DIFFERENCE_1("array sizes should have difference 1 ({0} != {1} + 1)"),
|
||||
|
|
@ -290,8 +291,10 @@ public enum LocalizedFormats implements Localizable {
|
|||
ZERO_FRACTION_TO_DIVIDE_BY("the fraction to divide by must not be zero: {0}/{1}"),
|
||||
ZERO_NORM("zero norm"),
|
||||
ZERO_NORM_FOR_ROTATION_AXIS("zero norm for rotation axis"),
|
||||
ZERO_NORM_FOR_ROTATION_DEFINING_VECTOR("zero norm for rotation defining vector");
|
||||
ZERO_NORM_FOR_ROTATION_DEFINING_VECTOR("zero norm for rotation defining vector"),
|
||||
ZERO_NOT_ALLOWED("zero not allowed here");
|
||||
|
||||
// CHECKSTYLE: resume JavadocVariable
|
||||
// CHECKSTYLE: resume MultipleVariableDeclarations
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -263,3 +263,4 @@ ZERO_FRACTION_TO_DIVIDE_BY = division par un nombre rationnel nul : {0}/{1}
|
|||
ZERO_NORM = norme nulle
|
||||
ZERO_NORM_FOR_ROTATION_AXIS = norme nulle pour un axe de rotation
|
||||
ZERO_NORM_FOR_ROTATION_DEFINING_VECTOR = norme nulle pour un axe de d\u00e9finition de rotation
|
||||
ZERO_NOT_ALLOWED = la valeur z\u00e9ro n''est pas autoris\u00e9e ici
|
||||
|
|
|
|||
|
|
@ -52,6 +52,9 @@ The <action> type attribute can be add,update,fix,remove.
|
|||
If the output is not quite correct, check for invisible trailing spaces!
|
||||
-->
|
||||
<release version="2.2" date="TBD" description="TBD">
|
||||
<action dev="luc" type="add" issue="MATH-400" due-to="J. Lewis Muir">
|
||||
Added support for Gaussian curve fitting.
|
||||
</action>
|
||||
<action dev="erans" type="update" issue="MATH-397">
|
||||
Modified "AbstractUnivariateRealOptimizer" to make it more similar to
|
||||
"BaseAbstractScalarOptimizer".
|
||||
|
|
|
|||
|
|
@ -0,0 +1,321 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
|
||||
package org.apache.commons.math.optimization.fitting;
|
||||
|
||||
import static org.junit.Assert.assertEquals;
|
||||
|
||||
import org.apache.commons.math.FunctionEvaluationException;
|
||||
import org.apache.commons.math.optimization.OptimizationException;
|
||||
import org.apache.commons.math.optimization.general.
|
||||
LevenbergMarquardtOptimizer;
|
||||
import org.junit.Test;
|
||||
|
||||
/**
|
||||
* Tests {@link GaussianFitter}.
|
||||
*
|
||||
* @since 2.2
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class GaussianFitterTest {
|
||||
/** Good data. */
|
||||
protected static final double[][] DATASET1 = new double[][] {
|
||||
{4.0254623, 531026.0},
|
||||
{4.02804905, 664002.0},
|
||||
{4.02934242, 787079.0},
|
||||
{4.03128248, 984167.0},
|
||||
{4.03386923, 1294546.0},
|
||||
{4.03580929, 1560230.0},
|
||||
{4.03839603, 1887233.0},
|
||||
{4.0396894, 2113240.0},
|
||||
{4.04162946, 2375211.0},
|
||||
{4.04421621, 2687152.0},
|
||||
{4.04550958, 2862644.0},
|
||||
{4.04744964, 3078898.0},
|
||||
{4.05003639, 3327238.0},
|
||||
{4.05132976, 3461228.0},
|
||||
{4.05326982, 3580526.0},
|
||||
{4.05585657, 3576946.0},
|
||||
{4.05779662, 3439750.0},
|
||||
{4.06038337, 3220296.0},
|
||||
{4.06167674, 3070073.0},
|
||||
{4.0636168, 2877648.0},
|
||||
{4.06620355, 2595848.0},
|
||||
{4.06749692, 2390157.0},
|
||||
{4.06943698, 2175960.0},
|
||||
{4.07202373, 1895104.0},
|
||||
{4.0733171, 1687576.0},
|
||||
{4.07525716, 1447024.0},
|
||||
{4.0778439, 1130879.0},
|
||||
{4.07978396, 904900.0},
|
||||
{4.08237071, 717104.0},
|
||||
{4.08366408, 620014.0}
|
||||
};
|
||||
/** Poor data: right of peak not symmetric with left of peak. */
|
||||
protected static final double[][] DATASET2 = new double[][] {
|
||||
{-20.15, 1523.0},
|
||||
{-19.65, 1566.0},
|
||||
{-19.15, 1592.0},
|
||||
{-18.65, 1927.0},
|
||||
{-18.15, 3089.0},
|
||||
{-17.65, 6068.0},
|
||||
{-17.15, 14239.0},
|
||||
{-16.65, 34124.0},
|
||||
{-16.15, 64097.0},
|
||||
{-15.65, 110352.0},
|
||||
{-15.15, 164742.0},
|
||||
{-14.65, 209499.0},
|
||||
{-14.15, 267274.0},
|
||||
{-13.65, 283290.0},
|
||||
{-13.15, 275363.0},
|
||||
{-12.65, 258014.0},
|
||||
{-12.15, 225000.0},
|
||||
{-11.65, 200000.0},
|
||||
{-11.15, 190000.0},
|
||||
{-10.65, 185000.0},
|
||||
{-10.15, 180000.0},
|
||||
{ -9.65, 179000.0},
|
||||
{ -9.15, 178000.0},
|
||||
{ -8.65, 177000.0},
|
||||
{ -8.15, 176000.0},
|
||||
{ -7.65, 175000.0},
|
||||
{ -7.15, 174000.0},
|
||||
{ -6.65, 173000.0},
|
||||
{ -6.15, 172000.0},
|
||||
{ -5.65, 171000.0},
|
||||
{ -5.15, 170000.0}
|
||||
};
|
||||
/** Poor data: long tails. */
|
||||
protected static final double[][] DATASET3 = new double[][] {
|
||||
{-90.15, 1513.0},
|
||||
{-80.15, 1514.0},
|
||||
{-70.15, 1513.0},
|
||||
{-60.15, 1514.0},
|
||||
{-50.15, 1513.0},
|
||||
{-40.15, 1514.0},
|
||||
{-30.15, 1513.0},
|
||||
{-20.15, 1523.0},
|
||||
{-19.65, 1566.0},
|
||||
{-19.15, 1592.0},
|
||||
{-18.65, 1927.0},
|
||||
{-18.15, 3089.0},
|
||||
{-17.65, 6068.0},
|
||||
{-17.15, 14239.0},
|
||||
{-16.65, 34124.0},
|
||||
{-16.15, 64097.0},
|
||||
{-15.65, 110352.0},
|
||||
{-15.15, 164742.0},
|
||||
{-14.65, 209499.0},
|
||||
{-14.15, 267274.0},
|
||||
{-13.65, 283290.0},
|
||||
{-13.15, 275363.0},
|
||||
{-12.65, 258014.0},
|
||||
{-12.15, 214073.0},
|
||||
{-11.65, 182244.0},
|
||||
{-11.15, 136419.0},
|
||||
{-10.65, 97823.0},
|
||||
{-10.15, 58930.0},
|
||||
{ -9.65, 35404.0},
|
||||
{ -9.15, 16120.0},
|
||||
{ -8.65, 9823.0},
|
||||
{ -8.15, 5064.0},
|
||||
{ -7.65, 2575.0},
|
||||
{ -7.15, 1642.0},
|
||||
{ -6.65, 1101.0},
|
||||
{ -6.15, 812.0},
|
||||
{ -5.65, 690.0},
|
||||
{ -5.15, 565.0},
|
||||
{ 5.15, 564.0},
|
||||
{ 15.15, 565.0},
|
||||
{ 25.15, 564.0},
|
||||
{ 35.15, 565.0},
|
||||
{ 45.15, 564.0},
|
||||
{ 55.15, 565.0},
|
||||
{ 65.15, 564.0},
|
||||
{ 75.15, 565.0}
|
||||
};
|
||||
/** Poor data: right of peak is missing. */
|
||||
protected static final double[][] DATASET4 = new double[][] {
|
||||
{-20.15, 1523.0},
|
||||
{-19.65, 1566.0},
|
||||
{-19.15, 1592.0},
|
||||
{-18.65, 1927.0},
|
||||
{-18.15, 3089.0},
|
||||
{-17.65, 6068.0},
|
||||
{-17.15, 14239.0},
|
||||
{-16.65, 34124.0},
|
||||
{-16.15, 64097.0},
|
||||
{-15.65, 110352.0},
|
||||
{-15.15, 164742.0},
|
||||
{-14.65, 209499.0},
|
||||
{-14.15, 267274.0},
|
||||
{-13.65, 283290.0}
|
||||
};
|
||||
/** Good data, but few points. */
|
||||
protected static final double[][] DATASET5 = new double[][] {
|
||||
{4.0254623, 531026.0},
|
||||
{4.03128248, 984167.0},
|
||||
{4.03839603, 1887233.0},
|
||||
{4.04421621, 2687152.0},
|
||||
{4.05132976, 3461228.0},
|
||||
{4.05326982, 3580526.0},
|
||||
{4.05779662, 3439750.0},
|
||||
{4.0636168, 2877648.0},
|
||||
{4.06943698, 2175960.0},
|
||||
{4.07525716, 1447024.0},
|
||||
{4.08237071, 717104.0},
|
||||
{4.08366408, 620014.0}
|
||||
};
|
||||
|
||||
/**
|
||||
* Basic.
|
||||
*
|
||||
* @throws OptimizationException in the event of a test case error
|
||||
* @throws FunctionEvaluationException in the event of a test case error
|
||||
*/
|
||||
@Test
|
||||
public void testFit01()
|
||||
throws OptimizationException, FunctionEvaluationException {
|
||||
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
|
||||
addDatasetToGaussianFitter(DATASET1, fitter);
|
||||
GaussianFunction fitFunction = fitter.fit();
|
||||
assertEquals(99200.86969833552, fitFunction.getA(), 1e-4);
|
||||
assertEquals(3410515.285208688, fitFunction.getB(), 1e-4);
|
||||
assertEquals(4.054928275302832, fitFunction.getC(), 1e-4);
|
||||
assertEquals(0.014609868872574, fitFunction.getD(), 1e-4);
|
||||
}
|
||||
|
||||
/**
|
||||
* Zero points is not enough observed points.
|
||||
*
|
||||
* @throws OptimizationException in the event of a test case error
|
||||
* @throws FunctionEvaluationException in the event of a test case error
|
||||
*/
|
||||
@Test(expected=IllegalArgumentException.class)
|
||||
public void testFit02()
|
||||
throws OptimizationException, FunctionEvaluationException {
|
||||
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
|
||||
fitter.fit();
|
||||
}
|
||||
|
||||
/**
|
||||
* Two points is not enough observed points.
|
||||
*
|
||||
* @throws OptimizationException in the event of a test case error
|
||||
* @throws FunctionEvaluationException in the event of a test case error
|
||||
*/
|
||||
@Test(expected=IllegalArgumentException.class)
|
||||
public void testFit03()
|
||||
throws OptimizationException, FunctionEvaluationException {
|
||||
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
|
||||
addDatasetToGaussianFitter(new double[][] {
|
||||
{4.0254623, 531026.0},
|
||||
{4.02804905, 664002.0}},
|
||||
fitter);
|
||||
fitter.fit();
|
||||
}
|
||||
|
||||
/**
|
||||
* Poor data: right of peak not symmetric with left of peak.
|
||||
*
|
||||
* @throws OptimizationException in the event of a test case error
|
||||
* @throws FunctionEvaluationException in the event of a test case error
|
||||
*/
|
||||
@Test
|
||||
public void testFit04()
|
||||
throws OptimizationException, FunctionEvaluationException {
|
||||
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
|
||||
addDatasetToGaussianFitter(DATASET2, fitter);
|
||||
GaussianFunction fitFunction = fitter.fit();
|
||||
assertEquals(-256534.689445631, fitFunction.getA(), 1e-4);
|
||||
assertEquals(481328.2181530679, fitFunction.getB(), 1e-4);
|
||||
assertEquals(-10.5217226891099, fitFunction.getC(), 1e-4);
|
||||
assertEquals(-7.64248239366800, fitFunction.getD(), 1e-4);
|
||||
}
|
||||
|
||||
/**
|
||||
* Poor data: long tails.
|
||||
*
|
||||
* @throws OptimizationException in the event of a test case error
|
||||
* @throws FunctionEvaluationException in the event of a test case error
|
||||
*/
|
||||
@Test
|
||||
public void testFit05()
|
||||
throws OptimizationException, FunctionEvaluationException {
|
||||
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
|
||||
addDatasetToGaussianFitter(DATASET3, fitter);
|
||||
GaussianFunction fitFunction = fitter.fit();
|
||||
assertEquals(491.6310079258938, fitFunction.getA(), 1e-4);
|
||||
assertEquals(283508.6800413632, fitFunction.getB(), 1e-4);
|
||||
assertEquals(-13.2966857238057, fitFunction.getC(), 1e-4);
|
||||
assertEquals(1.725590356962981, fitFunction.getD(), 1e-4);
|
||||
}
|
||||
|
||||
/**
|
||||
* Poor data: right of peak is missing.
|
||||
*
|
||||
* @throws OptimizationException in the event of a test case error
|
||||
* @throws FunctionEvaluationException in the event of a test case error
|
||||
*/
|
||||
@Test
|
||||
public void testFit06()
|
||||
throws OptimizationException, FunctionEvaluationException {
|
||||
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
|
||||
addDatasetToGaussianFitter(DATASET4, fitter);
|
||||
GaussianFunction fitFunction = fitter.fit();
|
||||
assertEquals(530.3649792355617, fitFunction.getA(), 1e-4);
|
||||
assertEquals(284517.0835567514, fitFunction.getB(), 1e-4);
|
||||
assertEquals(-13.5355534565105, fitFunction.getC(), 1e-4);
|
||||
assertEquals(1.512353018625465, fitFunction.getD(), 1e-4);
|
||||
}
|
||||
|
||||
/**
|
||||
* Basic with smaller dataset.
|
||||
*
|
||||
* @throws OptimizationException in the event of a test case error
|
||||
* @throws FunctionEvaluationException in the event of a test case error
|
||||
*/
|
||||
@Test
|
||||
public void testFit07()
|
||||
throws OptimizationException, FunctionEvaluationException {
|
||||
GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer());
|
||||
addDatasetToGaussianFitter(DATASET5, fitter);
|
||||
GaussianFunction fitFunction = fitter.fit();
|
||||
assertEquals(176748.1400947575, fitFunction.getA(), 1e-4);
|
||||
assertEquals(3361537.018813906, fitFunction.getB(), 1e-4);
|
||||
assertEquals(4.054949992747176, fitFunction.getC(), 1e-4);
|
||||
assertEquals(0.014192380137002, fitFunction.getD(), 1e-4);
|
||||
}
|
||||
|
||||
/**
|
||||
* Adds the specified points to specified <code>GaussianFitter</code>
|
||||
* instance.
|
||||
*
|
||||
* @param points data points where first dimension is a point index and
|
||||
* second dimension is an array of length two representing the point
|
||||
* with the first value corresponding to X and the second value
|
||||
* corresponding to Y
|
||||
* @param fitter fitter to which the points in <code>points</code> should be
|
||||
* added as observed points
|
||||
*/
|
||||
protected static void addDatasetToGaussianFitter(double[][] points,
|
||||
GaussianFitter fitter) {
|
||||
for (int i = 0; i < points.length; i++) {
|
||||
fitter.addObservedPoint(points[i][0], points[i][1]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,158 @@
|
|||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
|
||||
package org.apache.commons.math.optimization.fitting;
|
||||
|
||||
import org.apache.commons.math.FunctionEvaluationException;
|
||||
import org.apache.commons.math.exception.ZeroNotAllowedException;
|
||||
import org.apache.commons.math.optimization.OptimizationException;
|
||||
import org.apache.commons.math.optimization.fitting.CurveFitter;
|
||||
import org.apache.commons.math.optimization.general.
|
||||
LevenbergMarquardtOptimizer;
|
||||
import org.junit.Test;
|
||||
|
||||
import static org.junit.Assert.assertEquals;
|
||||
|
||||
/**
|
||||
* Tests {@link ParametricGaussianFunction}.
|
||||
*
|
||||
* @since 2.2
|
||||
* @version $Revision$ $Date$
|
||||
*/
|
||||
public class ParametricGaussianFunctionTest {
|
||||
/** Dataset 1 used by some test cases. */
|
||||
protected static final double[][] DATASET1 = new double[][] {
|
||||
{4.0254623, 531026.0},
|
||||
{4.02804905, 664002.0},
|
||||
{4.02934242, 787079.0},
|
||||
{4.03128248, 984167.0},
|
||||
{4.03386923, 1294546.0},
|
||||
{4.03580929, 1560230.0},
|
||||
{4.03839603, 1887233.0},
|
||||
{4.0396894, 2113240.0},
|
||||
{4.04162946, 2375211.0},
|
||||
{4.04421621, 2687152.0},
|
||||
{4.04550958, 2862644.0},
|
||||
{4.04744964, 3078898.0},
|
||||
{4.05003639, 3327238.0},
|
||||
{4.05132976, 3461228.0},
|
||||
{4.05326982, 3580526.0},
|
||||
{4.05585657, 3576946.0},
|
||||
{4.05779662, 3439750.0},
|
||||
{4.06038337, 3220296.0},
|
||||
{4.06167674, 3070073.0},
|
||||
{4.0636168, 2877648.0},
|
||||
{4.06620355, 2595848.0},
|
||||
{4.06749692, 2390157.0},
|
||||
{4.06943698, 2175960.0},
|
||||
{4.07202373, 1895104.0},
|
||||
{4.0733171, 1687576.0},
|
||||
{4.07525716, 1447024.0},
|
||||
{4.0778439, 1130879.0},
|
||||
{4.07978396, 904900.0},
|
||||
{4.08237071, 717104.0},
|
||||
{4.08366408, 620014.0}
|
||||
};
|
||||
|
||||
/**
|
||||
* Using not-so-good initial parameters.
|
||||
*
|
||||
* @throws OptimizationException in the event of a test case error
|
||||
* @throws FunctionEvaluationException in the event of a test case error
|
||||
*/
|
||||
@Test
|
||||
public void testFit01()
|
||||
throws OptimizationException, FunctionEvaluationException {
|
||||
CurveFitter fitter = new CurveFitter(new LevenbergMarquardtOptimizer());
|
||||
addDatasetToCurveFitter(DATASET1, fitter);
|
||||
double[] parameters = fitter.fit(new ParametricGaussianFunction(),
|
||||
new double[] {8.64753e3, 3.483323e6, 4.06322, 1.946857e-2});
|
||||
assertEquals(99200.94715858076, parameters[0], 1e-4);
|
||||
assertEquals(3410515.221897707, parameters[1], 1e-4);
|
||||
assertEquals(4.054928275257894, parameters[2], 1e-4);
|
||||
assertEquals(0.014609868499860, parameters[3], 1e-4);
|
||||
}
|
||||
|
||||
/**
|
||||
* Using eye-balled guesses for initial parameters.
|
||||
*
|
||||
* @throws OptimizationException in the event of a test case error
|
||||
* @throws FunctionEvaluationException in the event of a test case error
|
||||
*/
|
||||
@Test
|
||||
public void testFit02()
|
||||
throws OptimizationException, FunctionEvaluationException {
|
||||
CurveFitter fitter = new CurveFitter(new LevenbergMarquardtOptimizer());
|
||||
addDatasetToCurveFitter(DATASET1, fitter);
|
||||
double[] parameters = fitter.fit(new ParametricGaussianFunction(),
|
||||
new double[] {500000.0, 3500000.0, 4.055, 0.025479654});
|
||||
assertEquals(99200.81836264656, parameters[0], 1e-4);
|
||||
assertEquals(3410515.327151986, parameters[1], 1e-4);
|
||||
assertEquals(4.054928275377392, parameters[2], 1e-4);
|
||||
assertEquals(0.014609869119806, parameters[3], 1e-4);
|
||||
}
|
||||
|
||||
/**
|
||||
* The parameters array is null.
|
||||
*
|
||||
* @throws FunctionEvaluationException in the event of a test case error
|
||||
*/
|
||||
@Test(expected=IllegalArgumentException.class)
|
||||
public void testValue01() throws FunctionEvaluationException {
|
||||
ParametricGaussianFunction f = new ParametricGaussianFunction();
|
||||
f.value(0.0, null);
|
||||
}
|
||||
|
||||
/**
|
||||
* The parameters array length is not 4.
|
||||
*
|
||||
* @throws FunctionEvaluationException in the event of a test case error
|
||||
*/
|
||||
@Test(expected=IllegalArgumentException.class)
|
||||
public void testValue02() throws FunctionEvaluationException {
|
||||
ParametricGaussianFunction f = new ParametricGaussianFunction();
|
||||
f.value(0.0, new double[] {0.0, 1.0});
|
||||
}
|
||||
|
||||
/**
|
||||
* The parameters d is 0.
|
||||
*
|
||||
* @throws FunctionEvaluationException in the event of a test case error
|
||||
*/
|
||||
@Test(expected=ZeroNotAllowedException.class)
|
||||
public void testValue03() throws FunctionEvaluationException {
|
||||
ParametricGaussianFunction f = new ParametricGaussianFunction();
|
||||
f.value(0.0, new double[] {0.0, 1.0, 1.0, 0.0});
|
||||
}
|
||||
|
||||
/**
|
||||
* Adds the specified points to specified <code>CurveFitter</code> instance.
|
||||
*
|
||||
* @param points data points where first dimension is a point index and
|
||||
* second dimension is an array of length two representing the point
|
||||
* with the first value corresponding to X and the second value
|
||||
* corresponding to Y
|
||||
* @param fitter fitter to which the points in <code>points</code> should be
|
||||
* added as observed points
|
||||
*/
|
||||
protected static void addDatasetToCurveFitter(double[][] points,
|
||||
CurveFitter fitter) {
|
||||
for (int i = 0; i < points.length; i++) {
|
||||
fitter.addObservedPoint(points[i][0], points[i][1]);
|
||||
}
|
||||
}
|
||||
}
|
||||
Loading…
Reference in New Issue