Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

merged curve fitting from mantissa into commons-math 

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@786479 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-06-19 12:36:16 +00:00
parent 6463532544
commit 2a1842f3ef
25 changed files with 1250 additions and 1467 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

191
src/java/org/apache/commons/math/optimization/fitting/CurveFitter.java Normal file
View File

@ -0,0 +1,191 @@
/*
* 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.ArrayList;
import java.util.List;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
import org.apache.commons.math.analysis.MultivariateMatrixFunction;
import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.VectorialPointValuePair;
/** Fitter for parametric univariate real functions y = f(x).
* <p>When a univariate real function y = f(x) does depend on some
* unknown parameters p<sub>0</sub>, p<sub>1</sub> ... p<sub>n-1</sub>,
* this class can be used to find these parameters. It does this
* by <em>fitting</em> the curve so it remains very close to a set of
* observed points (x<sub>0</sub>, y<sub>0</sub>), (x<sub>1</sub>,
* y<sub>1</sub>) ... (x<sub>k-1</sub>, y<sub>k-1</sub>). This fitting
* is done by finding the parameters values that minimizes the objective
* function &sum;(y<sub>i</sub>-f(x<sub>i</sub>))<sup>2</sup>. This is
* really a least squares problem.</p>
* @version $Revision$ $Date$
* @since 2.0
*/
public class CurveFitter {
/** Optimizer to use for the fitting. */
private final DifferentiableMultivariateVectorialOptimizer optimizer;
/** Observed points. */
private final List<WeightedObservedPoint> observations;
/** Simple constructor.
* @param optimizer optimizer to use for the fitting
*/
public CurveFitter(final DifferentiableMultivariateVectorialOptimizer optimizer) {
this.optimizer = optimizer;
observations = new ArrayList<WeightedObservedPoint>();
}
/** Add an observed (x,y) point to the sample with unit weight.
* <p>Calling this method is equivalent to call
* <code>addObservedPoint(1.0, x, y)</code>.</p>
* @param x abscissa of the point
* @param y observed value of the point at x, after fitting we should
* have f(x) as close as possible to this value
* @see #addObservedPoint(double, double, double)
* @see #addObservedPoint(WeightedObservedPoint)
* @see #getObservations()
*/
public void addObservedPoint(double x, double y) {
addObservedPoint(1.0, x, y);
}
/** Add an observed weighted (x,y) point to the sample.
* @param weight weight of the observed point in the fit
* @param x abscissa of the point
* @param y observed value of the point at x, after fitting we should
* have f(x) as close as possible to this value
* @see #addObservedPoint(double, double)
* @see #addObservedPoint(WeightedObservedPoint)
* @see #getObservations()
*/
public void addObservedPoint(double weight, double x, double y) {
observations.add(new WeightedObservedPoint(weight, x, y));
}
/** Add an observed weighted (x,y) point to the sample.
* @param observed observed point to add
* @see #addObservedPoint(double, double)
* @see #addObservedPoint(double, double, double)
* @see #getObservations()
*/
public void addObservedPoint(WeightedObservedPoint observed) {
observations.add(observed);
}
/** Get the observed points.
* @return observed points
* @see #addObservedPoint(double, double)
* @see #addObservedPoint(double, double, double)
* @see #addObservedPoint(WeightedObservedPoint)
*/
public WeightedObservedPoint[] getObservations() {
return observations.toArray(new WeightedObservedPoint[observations.size()]);
}
/** Fit a curve.
* <p>This method compute the coefficients of the curve that best
* fit the sample of weighted pairs previously given through calls
* to the {@link #addWeightedPair addWeightedPair} method.</p>
* @param f parametric function to fit
* @param initialGuess first guess of the function parameters
* @return fitted parameters
* @exception FunctionEvaluationException if the objective function throws one during
* the search
* @exception OptimizationException if the algorithm failed to converge
* @exception IllegalArgumentException if the start point dimension is wrong
*/
public double[] fit(final ParametricRealFunction f,
final double[] initialGuess)
throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
// prepare least squares problem
double[] target = new double[observations.size()];
double[] weights = new double[observations.size()];
int i = 0;
for (WeightedObservedPoint point : observations) {
target[i] = point.getY();
weights[i] = point.getWeight();
++i;
}
// perform the fit
VectorialPointValuePair optimum =
optimizer.optimize(new TheoreticalValuesFunction(f), target, weights, initialGuess);
// extract the coefficients
return optimum.getPointRef();
}
/** Vectorial function computing function theoretical values. */
private class TheoreticalValuesFunction
implements DifferentiableMultivariateVectorialFunction {
/** Function to fit. */
private final ParametricRealFunction f;
/** Simple constructor.
* @param f function to fit.
*/
public TheoreticalValuesFunction(final ParametricRealFunction f) {
this.f = f;
}
/** {@inheritDoc} */
public MultivariateMatrixFunction jacobian() {
return new MultivariateMatrixFunction() {
public double[][] value(double[] point)
throws FunctionEvaluationException, IllegalArgumentException {
final double[][] jacobian = new double[observations.size()][];
int i = 0;
for (WeightedObservedPoint observed : observations) {
jacobian[i++] = f.gradient(observed.getX(), point);
}
return jacobian;
}
};
}
/** {@inheritDoc} */
public double[] value(double[] point)
throws FunctionEvaluationException, IllegalArgumentException {
// compute the residuals
final double[] values = new double[observations.size()];
int i = 0;
for (WeightedObservedPoint observed : observations) {
values[i++] = f.value(observed.getX(), point);
}
return values;
}
}
}

300
src/java/org/apache/commons/math/optimization/fitting/HarmonicCoefficientsGuesser.java Normal file
View File

@ -0,0 +1,300 @@
/*
* 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.optimization.OptimizationException;
/** This class guesses harmonic coefficients from a sample.
* <p>The algorithm used to guess the coefficients is as follows:</p>
* <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a,
* &omega; and &phi; such that f (t) = a cos (&omega; t + &phi;).
* </p>
*
* <p>From the analytical expression, we can compute two primitives :
* <pre>
* If2 (t) = &int; f<sup>2</sup> = a<sup>2</sup> &times; [t + S (t)] / 2
* If'2 (t) = &int; f'<sup>2</sup> = a<sup>2</sup> &omega;<sup>2</sup> &times; [t - S (t)] / 2
* where S (t) = sin (2 (&omega; t + &phi;)) / (2 &omega;)
* </pre>
* </p>
*
* <p>We can remove S between these expressions :
* <pre>
* If'2 (t) = a<sup>2</sup> &omega;<sup>2</sup> t - &omega;<sup>2</sup> If2 (t)
* </pre>
* </p>
*
* <p>The preceding expression shows that If'2 (t) is a linear
* combination of both t and If2 (t): If'2 (t) = A &times; t + B &times; If2 (t)
* </p>
*
* <p>From the primitive, we can deduce the same form for definite
* integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> :
* <pre>
* If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>) = A &times; (t<sub>i</sub> - t<sub>1</sub>) + B &times; (If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>))
* </pre>
* </p>
*
* <p>We can find the coefficients A and B that best fit the sample
* to this linear expression by computing the definite integrals for
* each sample points.
* </p>
*
* <p>For a bilinear expression z (x<sub>i</sub>, y<sub>i</sub>) = A &times; x<sub>i</sub> + B &times; y<sub>i</sub>, the
* coefficients A and B that minimize a least square criterion
* &sum; (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p>
* <pre>
*
* &sum;y<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
* A = ------------------------
* &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
*
* &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub>
* B = ------------------------
* &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
* </pre>
* </p>
*
*
* <p>In fact, we can assume both a and &omega; are positive and
* compute them directly, knowing that A = a<sup>2</sup> &omega;<sup>2</sup> and that
* B = - &omega;<sup>2</sup>. The complete algorithm is therefore:</p>
* <pre>
*
* for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute:
* f (t<sub>i</sub>)
* f' (t<sub>i</sub>) = (f (t<sub>i+1</sub>) - f(t<sub>i-1</sub>)) / (t<sub>i+1</sub> - t<sub>i-1</sub>)
* x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub>
* y<sub>i</sub> = &int; f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
* z<sub>i</sub> = &int; f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
* update the sums &sum;x<sub>i</sub>x<sub>i</sub>, &sum;y<sub>i</sub>y<sub>i</sub>, &sum;x<sub>i</sub>y<sub>i</sub>, &sum;x<sub>i</sub>z<sub>i</sub> and &sum;y<sub>i</sub>z<sub>i</sub>
* end for
*
* |--------------------------
* \ | &sum;y<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
* a = \ | ------------------------
* \| &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
*
*
* |--------------------------
* \ | &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
* &omega; = \ | ------------------------
* \| &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
*
* </pre>
* </p>
* <p>Once we know &omega;, we can compute:
* <pre>
* fc = &omega; f (t) cos (&omega; t) - f' (t) sin (&omega; t)
* fs = &omega; f (t) sin (&omega; t) + f' (t) cos (&omega; t)
* </pre>
* </p>
* <p>It appears that <code>fc = a &omega; cos (&phi;)</code> and
* <code>fs = -a &omega; sin (&phi;)</code>, so we can use these
* expressions to compute &phi;. The best estimate over the sample is
* given by averaging these expressions.
* </p>
* <p>Since integrals and means are involved in the preceding
* estimations, these operations run in O(n) time, where n is the
* number of measurements.</p>
* @version $Revision$ $Date$
* @since 2.0
*/
public class HarmonicCoefficientsGuesser {
/** Sampled observations. */
private final WeightedObservedPoint[] observations;
/** Guessed amplitude. */
private double a;
/** Guessed pulsation &omega;. */
private double omega;
/** Guessed phase &phi;. */
private double phi;
/** Simple constructor.
* @param observations sampled observations
*/
public HarmonicCoefficientsGuesser(WeightedObservedPoint[] observations) {
this.observations = observations.clone();
a = Double.NaN;
omega = Double.NaN;
}
/** Estimate a first guess of the coefficients.
* @exception OptimizationException if the sample is too short or if
* the first guess cannot be computed (when the elements under the
* square roots are negative).
* */
public void guess() throws OptimizationException {
sortObservations();
guessAOmega();
guessPhi();
}
/** Sort the observations with respect to the abscissa.
*/
private void sortObservations() {
// Since the samples are almost always already sorted, this
// method is implemented as an insertion sort that reorders the
// elements in place. Insertion sort is very efficient in this case.
WeightedObservedPoint curr = observations[0];
for (int j = 1; j < observations.length; ++j) {
WeightedObservedPoint prec = curr;
curr = observations[j];
if (curr.getX() < prec.getX()) {
// the current element should be inserted closer to the beginning
int i = j - 1;
WeightedObservedPoint mI = observations[i];
while ((i >= 0) && (curr.getX() < mI.getX())) {
observations[i + 1] = mI;
if (i-- != 0) {
mI = observations[i];
} else {
mI = null;
}
}
observations[i + 1] = curr;
curr = observations[j];
}
}
}
/** Estimate a first guess of the a and &omega; coefficients.
* @exception OptimizationException if the sample is too short or if
* the first guess cannot be computed (when the elements under the
* square roots are negative).
*/
private void guessAOmega() throws OptimizationException {
// initialize the sums for the linear model between the two integrals
double sx2 = 0.0;
double sy2 = 0.0;
double sxy = 0.0;
double sxz = 0.0;
double syz = 0.0;
double currentX = observations[0].getX();
double currentY = observations[0].getY();
double f2Integral = 0;
double fPrime2Integral = 0;
final double startX = currentX;
for (int i = 1; i < observations.length; ++i) {
// one step forward
final double previousX = currentX;
final double previousY = currentY;
currentX = observations[i].getX();
currentY = observations[i].getY();
// update the integrals of f<sup>2</sup> and f'<sup>2</sup>
// considering a linear model for f (and therefore constant f')
final double dx = currentX - previousX;
final double dy = currentY - previousY;
final double f2StepIntegral =
dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
final double fPrime2StepIntegral = dy * dy / dx;
final double x = currentX - startX;
f2Integral += f2StepIntegral;
fPrime2Integral += fPrime2StepIntegral;
sx2 += x * x;
sy2 += f2Integral * f2Integral;
sxy += x * f2Integral;
sxz += x * fPrime2Integral;
syz += f2Integral * fPrime2Integral;
}
// compute the amplitude and pulsation coefficients
double c1 = sy2 * sxz - sxy * syz;
double c2 = sxy * sxz - sx2 * syz;
double c3 = sx2 * sy2 - sxy * sxy;
if ((c1 / c2 < 0.0) || (c2 / c3 < 0.0)) {
throw new OptimizationException("unable to first guess the harmonic coefficients");
}
a = Math.sqrt(c1 / c2);
omega = Math.sqrt(c2 / c3);
}
/** Estimate a first guess of the &phi; coefficient.
*/
private void guessPhi() {
// initialize the means
double fcMean = 0.0;
double fsMean = 0.0;
double currentX = observations[0].getX();
double currentY = observations[0].getY();
for (int i = 1; i < observations.length; ++i) {
// one step forward
final double previousX = currentX;
final double previousY = currentY;
currentX = observations[i].getX();
currentY = observations[i].getY();
final double currentYPrime = (currentY - previousY) / (currentX - previousX);
double omegaX = omega * currentX;
double cosine = Math.cos(omegaX);
double sine = Math.sin(omegaX);
fcMean += omega * currentY * cosine - currentYPrime * sine;
fsMean += omega * currentY * sine + currentYPrime * cosine;
}
phi = Math.atan2(-fsMean, fcMean);
}
/** Get the guessed amplitude a.
* @return guessed amplitude a;
*/
public double getGuessedAmplitude() {
return a;
}
/** Get the guessed pulsation &omega;.
* @return guessed pulsation &omega;
*/
public double getGuessedPulsation() {
return omega;
}
/** Get the guessed phase &phi;.
* @return guessed phase &phi;
*/
public double getGuessedPhase() {
return phi;
}
}

132
src/java/org/apache/commons/math/optimization/fitting/HarmonicFitter.java Normal file
View File

@ -0,0 +1,132 @@
/*
* 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.MathRuntimeException;
import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
import org.apache.commons.math.optimization.OptimizationException;
/** This class implements a curve fitting specialized for sinusoids.
* <p>Harmonic fitting is a very simple case of curve fitting. The
* estimated coefficients are the amplitude a, the pulsation &omega; and
* the phase &phi;: <code>f (t) = a cos (&omega; t + &phi;)</code>. They are
* searched by a least square estimator initialized with a rough guess
* based on integrals.</p>
* @version $Revision$ $Date$
* @since 2.0
*/
public class HarmonicFitter {
/** Fitter for the coefficients. */
private final CurveFitter fitter;
/** Values for amplitude, pulsation &omega; and phase &phi;. */
private double[] parameters;
/** Simple constructor.
* @param optimizer optimizer to use for the fitting
*/
public HarmonicFitter(final DifferentiableMultivariateVectorialOptimizer optimizer) {
this.fitter = new CurveFitter(optimizer);
parameters = null;
}
/** Simple constructor.
* <p>This constructor can be used when a first guess of the
* coefficients is already known.</p>
* @param optimizer optimizer to use for the fitting
* @param initialGuess guessed values for amplitude (index 0),
* pulsation &omega; (index 1) and phase &phi; (index 2)
*/
public HarmonicFitter(final DifferentiableMultivariateVectorialOptimizer optimizer,
final double[] initialGuess) {
this.fitter = new CurveFitter(optimizer);
this.parameters = initialGuess.clone();
}
/** Add an observed weighted (x,y) point to the sample.
* @param weight weight of the observed point in the fit
* @param x abscissa of the point
* @param y observed value of the point at x, after fitting we should
* have P(x) as close as possible to this value
*/
public void addObservedPoint(double weight, double x, double y) {
fitter.addObservedPoint(weight, x, y);
}
/** Fit an harmonic function to the observed points.
* @return harmonic function best fitting the observed points
* @throws OptimizationException if the sample is too short or if
* the first guess cannot be computed
*/
public HarmonicFunction fit() throws OptimizationException {
try {
// shall we compute the first guess of the parameters ourselves ?
if (parameters == null) {
final WeightedObservedPoint[] observations = fitter.getObservations();
if (observations.length < 4) {
throw new OptimizationException("sample contains {0} observed points, at least {1} are required",
observations.length, 4);
}
HarmonicCoefficientsGuesser guesser = new HarmonicCoefficientsGuesser(observations);
guesser.guess();
parameters = new double[] {
guesser.getGuessedAmplitude(),
guesser.getGuessedPulsation(),
guesser.getGuessedPhase()
};
}
double[] fitted = fitter.fit(new ParametricHarmonicFunction(), parameters);
return new HarmonicFunction(fitted[0], fitted[1], fitted[2]);
} catch (FunctionEvaluationException fee) {
// this should never happen
throw MathRuntimeException.createInternalError(fee);
}
}
/** Parametric harmonic function. */
private static class ParametricHarmonicFunction implements ParametricRealFunction {
/** {@inheritDoc} */
public double value(double x, double[] parameters) {
final double a = parameters[0];
final double omega = parameters[1];
final double phi = parameters[2];
return a * Math.cos(omega * x + phi);
}
/** {@inheritDoc} */
public double[] gradient(double x, double[] parameters) {
final double a = parameters[0];
final double omega = parameters[1];
final double phi = parameters[2];
final double alpha = omega * x + phi;
final double cosAlpha = Math.cos(alpha);
final double sinAlpha = Math.sin(alpha);
return new double[] { cosAlpha, -a * x * sinAlpha, -a * sinAlpha };
}
}
}

79
src/java/org/apache/commons/math/optimization/fitting/HarmonicFunction.java Normal file
View File

@ -0,0 +1,79 @@
/*
* 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.analysis.DifferentiableUnivariateRealFunction;
/** Harmonic function of the form <code>f (t) = a cos (&omega; t + &phi;)</code>.
* @version $Revision$ $Date$
* @since 2.0
*/
public class HarmonicFunction implements DifferentiableUnivariateRealFunction {
/** Amplitude a. */
private final double a;
/** Pulsation &omega;. */
private final double omega;
/** Phase &phi;. */
private final double phi;
/** Simple constructor.
* @param a amplitude
* @param omega pulsation
* @param phi phase
*/
public HarmonicFunction(double a, double omega, double phi) {
this.a = a;
this.omega = omega;
this.phi = phi;
}
/** {@inheritDoc} */
public double value(double x) {
return a * Math.cos(omega * x + phi);
}
/** {@inheritDoc} */
public HarmonicFunction derivative() {
return new HarmonicFunction(a * omega, omega, phi + Math.PI / 2);
}
/** Get the amplitude a.
* @return amplitude a;
*/
public double getAmplitude() {
return a;
}
/** Get the pulsation &omega;.
* @return pulsation &omega;
*/
public double getPulsation() {
return omega;
}
/** Get the phase &phi;.
* @return phase &phi;
*/
public double getPhase() {
return phi;
}
}

50
src/java/org/apache/commons/math/optimization/fitting/ParametricRealFunction.java Normal file
View File

@ -0,0 +1,50 @@
/*
* 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;
/**
* An interface representing a real function that depends on one independent
* variable plus some extra parameters.
*
* @version $Revision$ $Date$
*/
public interface ParametricRealFunction {
/**
* Compute the value of the function.
* @param x the point for which the function value should be computed
* @param parameters function parameters
* @return the value
* @throws FunctionEvaluationException if the function evaluation fails
*/
public double value(double x, double[] parameters)
throws FunctionEvaluationException;
/**
* Compute the gradient of the function with respect to its parameters.
* @param x the point for which the function value should be computed
* @param parameters function parameters
* @return the value
* @throws FunctionEvaluationException if the function evaluation fails
*/
public double[] gradient(double x, double[] parameters)
throws FunctionEvaluationException;
}

103
src/java/org/apache/commons/math/optimization/fitting/PolynomialFitter.java Normal file
View File

@ -0,0 +1,103 @@
/*
* 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.MathRuntimeException;
import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
import org.apache.commons.math.optimization.OptimizationException;
/** This class implements a curve fitting specialized for polynomials.
* <p>Polynomial fitting is a very simple case of curve fitting. The
* estimated coefficients are the polynomial coefficients. They are
* searched by a least square estimator.</p>
* @version $Revision$ $Date$
* @since 2.0
*/
public class PolynomialFitter {
/** Fitter for the coefficients. */
private final CurveFitter fitter;
/** Polynomial degree. */
private final int degree;
/** Simple constructor.
* <p>The polynomial fitter built this way are complete polynomials,
* ie. a n-degree polynomial has n+1 coefficients.</p>
* @param degree maximal degree of the polynomial
* @param optimizer optimizer to use for the fitting
*/
public PolynomialFitter(int degree, final DifferentiableMultivariateVectorialOptimizer optimizer) {
this.fitter = new CurveFitter(optimizer);
this.degree = degree;
}
/** Add an observed weighted (x,y) point to the sample.
* @param weight weight of the observed point in the fit
* @param x abscissa of the point
* @param y observed value of the point at x, after fitting we should
* have P(x) as close as possible to this value
*/
public void addObservedPoint(double weight, double x, double y) {
fitter.addObservedPoint(weight, x, y);
}
/** Get the polynomial fitting the weighted (x, y) points.
* @return polynomial function best fitting the observed points
* @exception OptimizationException if the algorithm failed to converge
*/
public PolynomialFunction fit()
throws OptimizationException {
try {
return new PolynomialFunction(fitter.fit(new ParametricPolynomial(), new double[degree + 1]));
} catch (FunctionEvaluationException fee) {
// this should never happen
throw MathRuntimeException.createInternalError(fee);
}
}
/** Dedicated parametric polynomial class. */
private static class ParametricPolynomial implements ParametricRealFunction {
/** {@inheritDoc} */
public double[] gradient(double x, double[] parameters)
throws FunctionEvaluationException {
final double[] gradient = new double[parameters.length];
double xn = 1.0;
for (int i = 0; i < parameters.length; ++i) {
gradient[i] = xn;
xn *= x;
}
return gradient;
}
/** {@inheritDoc} */
public double value(final double x, final double[] parameters) {
double y = 0;
for (int i = parameters.length - 1; i >= 0; --i) {
y = y * x + parameters[i];
}
return y;
}
}
}

75
src/java/org/apache/commons/math/optimization/fitting/WeightedObservedPoint.java Normal file
View File

@ -0,0 +1,75 @@
/*
* 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;
/** This class is a simple container for weighted observed point in
* {@link CurveFitter curve fitting}.
* <p>Instances of this class are guaranteed to be immutable.</p>
* @version $Revision$ $Date$
* @since 2.0
*/
public class WeightedObservedPoint implements Serializable {
/** Serializable version id. */
private static final long serialVersionUID = 5306874947404636157L;
/** Weight of the measurement in the fitting process. */
private final double weight;
/** Abscissa of the point. */
private final double x;
/** Observed value of the function at x. */
private final double y;
/** Simple constructor.
* @param weight weight of the measurement in the fitting process
* @param x abscissa of the measurement
* @param y ordinate of the measurement
*/
public WeightedObservedPoint(final double weight, final double x, final double y) {
this.weight = weight;
this.x = x;
this.y = y;
}
/** Get the weight of the measurement in the fitting process.
* @return weight of the measurement in the fitting process
*/
public double getWeight() {
return weight;
}
/** Get the abscissa of the point.
* @return abscissa of the point
*/
public double getX() {
return x;
}
/** Get the observed value of the function at x.
* @return observed value of the function at x
*/
public double getY() {
return y;
}
}

0
src/mantissa/src/org/spaceroots/mantissa/fitting/package.html → src/java/org/apache/commons/math/optimization/fitting/package.html
View File

222
src/mantissa/src/org/spaceroots/mantissa/fitting/AbstractCurveFitter.java
View File

@ -1,222 +0,0 @@
// 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.spaceroots.mantissa.fitting;
import java.io.Serializable;
import java.util.ArrayList;
import java.util.List;
import org.spaceroots.mantissa.estimation.*;
/** This class is the base class for all curve fitting classes in the package.
* <p>This class handles all common features of curve fitting like the
* sample points handling. It declares two methods ({@link
* #valueAt} and {@link #partial}) which should be implemented by
* sub-classes to define the precise shape of the curve they
* represent.</p>
* @version $Id$
* @author L. Maisonobe
*/
public abstract class AbstractCurveFitter
implements EstimationProblem, Serializable {
/** Simple constructor.
* @param n number of coefficients in the underlying function
* @param estimator estimator to use for the fitting
*/
protected AbstractCurveFitter(int n, Estimator estimator) {
coefficients = new EstimatedParameter[n];
measurements = new ArrayList();
this.estimator = estimator;
}
/** Simple constructor.
* @param coefficients first estimate of the coefficients. A
* reference to this array is hold by the newly created object. Its
* elements will be adjusted during the fitting process and they will
* be set to the adjusted coefficients at the end.
* @param estimator estimator to use for the fitting
*/
protected AbstractCurveFitter(EstimatedParameter[] coefficients,
Estimator estimator) {
this.coefficients = coefficients;
measurements = new ArrayList();
this.estimator = estimator;
}
/** Add a weighted (x,y) pair to the sample.
* @param weight weight of this pair in the fit
* @param x abscissa
* @param y ordinate, we have <code>y = f (x)</code>
*/
public void addWeightedPair(double weight, double x, double y) {
measurements.add(new FitMeasurement(weight, x, y));
}
/** Perform the fitting.
* <p>This method compute the coefficients of the curve that best
* fit the sample of weighted pairs previously given through calls
* to the {@link #addWeightedPair addWeightedPair} method.</p>
* @return coefficients of the curve
* @exception EstimationException if the fitting is not possible
* (for example if the sample has to few independant points)
*/
public double[] fit()
throws EstimationException {
// perform the fit
estimator.estimate(this);
// extract the coefficients
double[] fittedCoefficients = new double[coefficients.length];
for (int i = 0; i < coefficients.length; ++i) {
fittedCoefficients[i] = coefficients[i].getEstimate();
}
return fittedCoefficients;
}
public WeightedMeasurement[] getMeasurements() {
return (WeightedMeasurement[]) measurements.toArray(new FitMeasurement[measurements.size()]);
}
/** Get the unbound parameters of the problem.
* For a curve fitting, none of the function coefficient is bound.
* @return unbound parameters
*/
public EstimatedParameter[] getUnboundParameters() {
return (EstimatedParameter[]) coefficients.clone();
}
/** Get all the parameters of the problem.
* @return parameters
*/
public EstimatedParameter[] getAllParameters() {
return (EstimatedParameter[]) coefficients.clone();
}
/** Utility method to sort the measurements with respect to the abscissa.
* <p>This method is provided as a utility for derived classes. As
* an example, the {@link HarmonicFitter} class needs it in order to
* compute a first guess of the coefficients to initialize the
* estimation algorithm.</p>
*/
protected void sortMeasurements() {
// Since the samples are almost always already sorted, this
// method is implemented as an insertion sort that reorders the
// elements in place. Insertion sort is very efficient in this case.
FitMeasurement curr = (FitMeasurement) measurements.get(0);
for (int j = 1; j < measurements.size (); ++j) {
FitMeasurement prec = curr;
curr = (FitMeasurement) measurements.get(j);
if (curr.x < prec.x) {
// the current element should be inserted closer to the beginning
int i = j - 1;
FitMeasurement mI = (FitMeasurement) measurements.get(i);
while ((i >= 0) && (curr.x < mI.x)) {
measurements.set(i + 1, mI);
if (i-- != 0) {
mI = (FitMeasurement) measurements.get(i);
} else {
mI = null;
}
}
measurements.set(i + 1, curr);
curr = (FitMeasurement) measurements.get(j);
}
}
}
/** Get the value of the function at x according to the current parameters value.
* @param x abscissa at which the theoretical value is requested
* @return theoretical value at x
*/
public abstract double valueAt(double x);
/** Get the derivative of the function at x with respect to parameter p.
* @param x abscissa at which the partial derivative is requested
* @param p parameter with respect to which the derivative is requested
* @return partial derivative
*/
public abstract double partial(double x, EstimatedParameter p);
/** This class represents the fit measurements.
* One measurement is a weighted pair (x, y), where <code>y = f
* (x)</code> is the value of the function at x abscissa. This class
* is an inner class because the methods related to the computation
* of f values and derivative are proveded by the fitter
* implementations.
*/
public class FitMeasurement
extends WeightedMeasurement {
/** Simple constructor.
* @param weight weight of the measurement in the fitting process
* @param x abscissa of the measurement
* @param y ordinate of the measurement
*/
public FitMeasurement(double weight, double x, double y) {
super(weight, y);
this.x = x;
}
/** Get the value of the fitted function at x.
* @return theoretical value at the measurement abscissa
*/
public double getTheoreticalValue() {
return valueAt(x);
}
/** Partial derivative with respect to one function coefficient.
* @param p parameter with respect to which the derivative is requested
* @return partial derivative
*/
public double getPartial(EstimatedParameter p) {
return partial(x, p);
}
/** Abscissa of the measurement. */
public final double x;
private static final long serialVersionUID = -2682582852369995960L;
}
/** Coefficients of the function */
protected EstimatedParameter[] coefficients;
/** Measurements vector */
protected List measurements;
/** Estimator for the fitting problem. */
private Estimator estimator;
}

75
src/mantissa/src/org/spaceroots/mantissa/fitting/F2FP2Iterator.java
View File

@ -1,75 +0,0 @@
// 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.spaceroots.mantissa.fitting;
import java.io.Serializable;
import org.spaceroots.mantissa.functions.FunctionException;
import org.spaceroots.mantissa.functions.ExhaustedSampleException;
import org.spaceroots.mantissa.functions.vectorial.SampledFunctionIterator;
import org.spaceroots.mantissa.functions.vectorial.VectorialValuedPair;
/** This class provides sampled values of the function t -> [f(t)^2, f'(t)^2].
* This class is a helper class used to compute a first guess of the
* harmonic coefficients of a function <code>f (t) = a cos (omega t +
* phi)</code>.
* @see FFPIterator
* @see HarmonicCoefficientsGuesser
* @version $Id$
* @author L. Maisonobe
*/
class F2FP2Iterator
implements SampledFunctionIterator, Serializable {
public F2FP2Iterator(AbstractCurveFitter.FitMeasurement[] measurements) {
ffpIterator = new FFPIterator(measurements);
}
public int getDimension() {
return 2;
}
public boolean hasNext() {
return ffpIterator.hasNext();
}
public VectorialValuedPair nextSamplePoint()
throws ExhaustedSampleException, FunctionException {
// get the raw values from the underlying FFPIterator
VectorialValuedPair point = ffpIterator.nextSamplePoint();
double[] y = point.y;
// square the values
return new VectorialValuedPair(point.x,
new double[] {
y[0] * y[0], y[1] * y[1]
});
}
private FFPIterator ffpIterator;
private static final long serialVersionUID = -8113110433795298072L;
}

100
src/mantissa/src/org/spaceroots/mantissa/fitting/FFPIterator.java
View File

@ -1,100 +0,0 @@
// 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.spaceroots.mantissa.fitting;
import java.io.Serializable;
import org.spaceroots.mantissa.functions.FunctionException;
import org.spaceroots.mantissa.functions.ExhaustedSampleException;
import org.spaceroots.mantissa.functions.vectorial.SampledFunctionIterator;
import org.spaceroots.mantissa.functions.vectorial.VectorialValuedPair;
/** This class provides sampled values of the function t -> [f(t), f'(t)].
* This class is a helper class used to compute a first guess of the
* harmonic coefficients of a function <code>f (t) = a cos (omega t +
* phi)</code>.
* @see F2FP2Iterator
* @see HarmonicCoefficientsGuesser
* @version $Id$
* @author L. Maisonobe
*/
class FFPIterator
implements SampledFunctionIterator, Serializable {
public FFPIterator(AbstractCurveFitter.FitMeasurement[] measurements) {
this.measurements = measurements;
// initialize the points of the raw sample
current = measurements[0];
currentY = current.getMeasuredValue();
next = measurements[1];
nextY = next.getMeasuredValue();
nextIndex = 2;
}
public int getDimension() {
return 2;
}
public boolean hasNext() {
return nextIndex < measurements.length;
}
public VectorialValuedPair nextSamplePoint()
throws ExhaustedSampleException, FunctionException {
if (nextIndex >= measurements.length) {
throw new ExhaustedSampleException(measurements.length);
}
// shift the points
previous = current;
previousY = currentY;
current = next;
currentY = nextY;
next = measurements[nextIndex++];
nextY = next.getMeasuredValue();
// return the two dimensions vector [f(x), f'(x)]
double[] table = new double[2];
table[0] = currentY;
table[1] = (nextY - previousY) / (next.x - previous.x);
return new VectorialValuedPair(current.x, table);
}
private AbstractCurveFitter.FitMeasurement[] measurements;
private int nextIndex;
private AbstractCurveFitter.FitMeasurement previous;
private double previousY;
private AbstractCurveFitter.FitMeasurement current;
private double nextY;
private AbstractCurveFitter.FitMeasurement next;
private double currentY;
private static final long serialVersionUID = -3187229691615380125L;
}

269
src/mantissa/src/org/spaceroots/mantissa/fitting/HarmonicCoefficientsGuesser.java
View File

@ -1,269 +0,0 @@
// 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.spaceroots.mantissa.fitting;
import java.io.Serializable;
import org.spaceroots.mantissa.functions.FunctionException;
import org.spaceroots.mantissa.functions.ExhaustedSampleException;
import org.spaceroots.mantissa.functions.vectorial.SampledFunctionIterator;
import org.spaceroots.mantissa.functions.vectorial.VectorialValuedPair;
import org.spaceroots.mantissa.quadrature.vectorial.EnhancedSimpsonIntegratorSampler;
import org.spaceroots.mantissa.estimation.EstimationException;
/** This class guesses harmonic coefficients from a sample.
* <p>The algorithm used to guess the coefficients is as follows:</p>
* <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a,
* &omega; and &phi; such that f (t) = a cos (&omega; t + &phi;).
* </p>
*
* <p>From the analytical expression, we can compute two primitives :
* <pre>
* If2 (t) = &int; f<sup>2</sup> = a<sup>2</sup> &times; [t + S (t)] / 2
* If'2 (t) = &int; f'<sup>2</sup> = a<sup>2</sup> &omega;<sup>2</sup> &times; [t - S (t)] / 2
* where S (t) = sin (2 (&omega; t + &phi;)) / (2 &omega;)
* </pre>
* </p>
*
* <p>We can remove S between these expressions :
* <pre>
* If'2 (t) = a<sup>2</sup> &omega;<sup>2</sup> t - &omega;<sup>2</sup> If2 (t)
* </pre>
* </p>
*
* <p>The preceding expression shows that If'2 (t) is a linear
* combination of both t and If2 (t): If'2 (t) = A &times; t + B &times; If2 (t)
* </p>
*
* <p>From the primitive, we can deduce the same form for definite
* integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> :
* <pre>
* If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>) = A &times; (t<sub>i</sub> - t<sub>1</sub>) + B &times; (If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>))
* </pre>
* </p>
*
* <p>We can find the coefficients A and B that best fit the sample
* to this linear expression by computing the definite integrals for
* each sample points.
* </p>
*
* <p>For a bilinear expression z (x<sub>i</sub>, y<sub>i</sub>) = A &times; x<sub>i</sub> + B &times; y<sub>i</sub>, the
* coefficients A and B that minimize a least square criterion
* &sum; (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p>
* <pre>
*
* &sum;y<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
* A = ------------------------
* &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
*
* &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub>
* B = ------------------------
* &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
* </pre>
* </p>
*
*
* <p>In fact, we can assume both a and &omega; are positive and
* compute them directly, knowing that A = a<sup>2</sup> &omega;<sup>2</sup> and that
* B = - &omega;<sup>2</sup>. The complete algorithm is therefore:</p>
* <pre>
*
* for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute:
* f (t<sub>i</sub>)
* f' (t<sub>i</sub>) = (f (t<sub>i+1</sub>) - f(t<sub>i-1</sub>)) / (t<sub>i+1</sub> - t<sub>i-1</sub>)
* x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub>
* y<sub>i</sub> = &int; f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
* z<sub>i</sub> = &int; f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
* update the sums &sum;x<sub>i</sub>x<sub>i</sub>, &sum;y<sub>i</sub>y<sub>i</sub>, &sum;x<sub>i</sub>y<sub>i</sub>, &sum;x<sub>i</sub>z<sub>i</sub> and &sum;y<sub>i</sub>z<sub>i</sub>
* end for
*
* |--------------------------
* \ | &sum;y<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
* a = \ | ------------------------
* \| &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
*
*
* |--------------------------
* \ | &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
* &omega; = \ | ------------------------
* \| &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
*
* </pre>
* </p>
* <p>Once we know &omega;, we can compute:
* <pre>
* fc = &omega; f (t) cos (&omega; t) - f' (t) sin (&omega; t)
* fs = &omega; f (t) sin (&omega; t) + f' (t) cos (&omega; t)
* </pre>
* </p>
* <p>It appears that <code>fc = a &omega; cos (&phi;)</code> and
* <code>fs = -a &omega; sin (&phi;)</code>, so we can use these
* expressions to compute &phi;. The best estimate over the sample is
* given by averaging these expressions.
* </p>
* <p>Since integrals and means are involved in the preceding
* estimations, these operations run in O(n) time, where n is the
* number of measurements.</p>
* @version $Id$
* @author L. Maisonobe
*/
public class HarmonicCoefficientsGuesser
implements Serializable{
public HarmonicCoefficientsGuesser(AbstractCurveFitter.FitMeasurement[] measurements) {
this.measurements =
(AbstractCurveFitter.FitMeasurement[]) measurements.clone();
a = Double.NaN;
omega = Double.NaN;
}
/** Estimate a first guess of the coefficients.
* @exception ExhaustedSampleException if the sample is exhausted.
* @exception FunctionException if the integrator throws one.
* @exception EstimationException if the sample is too short or if
* the first guess cannot be computed (when the elements under the
* square roots are negative).
* */
public void guess()
throws ExhaustedSampleException, FunctionException, EstimationException {
guessAOmega();
guessPhi();
}
/** Estimate a first guess of the a and &omega; coefficients.
* @exception ExhaustedSampleException if the sample is exhausted.
* @exception FunctionException if the integrator throws one.
* @exception EstimationException if the sample is too short or if
* the first guess cannot be computed (when the elements under the
* square roots are negative).
*/
private void guessAOmega()
throws ExhaustedSampleException, FunctionException, EstimationException {
// initialize the sums for the linear model between the two integrals
double sx2 = 0.0;
double sy2 = 0.0;
double sxy = 0.0;
double sxz = 0.0;
double syz = 0.0;
// build the integrals sampler
F2FP2Iterator iter = new F2FP2Iterator(measurements);
SampledFunctionIterator sampler =
new EnhancedSimpsonIntegratorSampler(iter);
VectorialValuedPair p0 = sampler.nextSamplePoint();
double p0X = p0.x;
double[] p0Y = p0.y;
// get the points for the linear model
while (sampler.hasNext()) {
VectorialValuedPair point = sampler.nextSamplePoint();
double pX = point.x;
double[] pY = point.y;
double x = pX - p0X;
double y = pY[0] - p0Y[0];
double z = pY[1] - p0Y[1];
sx2 += x * x;
sy2 += y * y;
sxy += x * y;
sxz += x * z;
syz += y * z;
}
// compute the amplitude and pulsation coefficients
double c1 = sy2 * sxz - sxy * syz;
double c2 = sxy * sxz - sx2 * syz;
double c3 = sx2 * sy2 - sxy * sxy;
if ((c1 / c2 < 0.0) || (c2 / c3 < 0.0)) {
throw new EstimationException("unable to guess a first estimate");
}
a = Math.sqrt(c1 / c2);
omega = Math.sqrt(c2 / c3);
}
/** Estimate a first guess of the &phi; coefficient.
* @exception ExhaustedSampleException if the sample is exhausted.
* @exception FunctionException if the sampler throws one.
*/
private void guessPhi()
throws ExhaustedSampleException, FunctionException {
SampledFunctionIterator iter = new FFPIterator(measurements);
// initialize the means
double fcMean = 0.0;
double fsMean = 0.0;
while (iter.hasNext()) {
VectorialValuedPair point = iter.nextSamplePoint();
double omegaX = omega * point.x;
double cosine = Math.cos(omegaX);
double sine = Math.sin(omegaX);
fcMean += omega * point.y[0] * cosine - point.y[1] * sine;
fsMean += omega * point.y[0] * sine + point.y[1] * cosine;
}
phi = Math.atan2(-fsMean, fcMean);
}
public double getOmega() {
return omega;
}
public double getA() {
return a;
}
public double getPhi() {
return phi;
}
private AbstractCurveFitter.FitMeasurement[] measurements;
private double a;
private double omega;
private double phi;
private static final long serialVersionUID = 2400399048702758814L;
}

168
src/mantissa/src/org/spaceroots/mantissa/fitting/HarmonicFitter.java
View File

@ -1,168 +0,0 @@
// 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.spaceroots.mantissa.fitting;
import org.spaceroots.mantissa.estimation.EstimatedParameter;
import org.spaceroots.mantissa.estimation.EstimationException;
import org.spaceroots.mantissa.estimation.Estimator;
import org.spaceroots.mantissa.estimation.GaussNewtonEstimator;
import org.spaceroots.mantissa.functions.ExhaustedSampleException;
import org.spaceroots.mantissa.functions.FunctionException;
/** This class implements a curve fitting specialized for sinusoids.
* <p>Harmonic fitting is a very simple case of curve fitting. The
* estimated coefficients are the amplitude a, the pulsation omega and
* the phase phi: <code>f (t) = a cos (omega t + phi)</code>. They are
* searched by a least square estimator initialized with a rough guess
* based on integrals.</p>
* <p>This class <emph>is by no means optimized</emph>, neither versus
* space nor versus time performance.</p>
* @version $Id$
* @author L. Maisonobe
*/
public class HarmonicFitter
extends AbstractCurveFitter {
/** Simple constructor.
* @param estimator estimator to use for the fitting
*/
public HarmonicFitter(Estimator estimator) {
super(3, estimator);
coefficients[0] = new EstimatedParameter("a", 2.0 * Math.PI);
coefficients[1] = new EstimatedParameter("omega", 0.0);
coefficients[2] = new EstimatedParameter("phi", 0.0);
firstGuessNeeded = true;
}
/**
* Simple constructor.
* <p>This constructor can be used when a first estimate of the
* coefficients is already known.</p>
* @param coefficients first estimate of the coefficients.
* A reference to this array is hold by the newly created
* object. Its elements will be adjusted during the fitting process
* and they will be set to the adjusted coefficients at the end.
* @param estimator estimator to use for the fitting
*/
public HarmonicFitter(EstimatedParameter[] coefficients,
Estimator estimator) {
super(coefficients, estimator);
firstGuessNeeded = false;
}
public double[] fit()
throws EstimationException {
if (firstGuessNeeded) {
if (measurements.size() < 4) {
throw new EstimationException("sample must contain at least {0} points",
new String[] {
Integer.toString(4)
});
}
sortMeasurements();
try {
HarmonicCoefficientsGuesser guesser =
new HarmonicCoefficientsGuesser((FitMeasurement[]) getMeasurements());
guesser.guess();
coefficients[0].setEstimate(guesser.getA());
coefficients[1].setEstimate(guesser.getOmega());
coefficients[2].setEstimate(guesser.getPhi());
} catch(ExhaustedSampleException e) {
throw new EstimationException(e);
} catch(FunctionException e) {
throw new EstimationException(e);
}
firstGuessNeeded = false;
}
return super.fit();
}
/** Get the current amplitude coefficient estimate.
* Get a, where <code>f (t) = a cos (omega t + phi)</code>
* @return current amplitude coefficient estimate
*/
public double getAmplitude() {
return coefficients[0].getEstimate();
}
/** Get the current pulsation coefficient estimate.
* Get omega, where <code>f (t) = a cos (omega t + phi)</code>
* @return current pulsation coefficient estimate
*/
public double getPulsation() {
return coefficients[1].getEstimate();
}
/** Get the current phase coefficient estimate.
* Get phi, where <code>f (t) = a cos (omega t + phi)</code>
* @return current phase coefficient estimate
*/
public double getPhase() {
return coefficients[2].getEstimate();
}
/** Get the value of the function at x according to the current parameters value.
* @param x abscissa at which the theoretical value is requested
* @return theoretical value at x
*/
public double valueAt(double x) {
double a = coefficients[0].getEstimate();
double omega = coefficients[1].getEstimate();
double phi = coefficients[2].getEstimate();
return a * Math.cos(omega * x + phi);
}
/** Get the derivative of the function at x with respect to parameter p.
* @param x abscissa at which the partial derivative is requested
* @param p parameter with respect to which the derivative is requested
* @return partial derivative
*/
public double partial(double x, EstimatedParameter p) {
double a = coefficients[0].getEstimate();
double omega = coefficients[1].getEstimate();
double phi = coefficients[2].getEstimate();
if (p == coefficients[0]) {
return Math.cos(omega * x + phi);
} else if (p == coefficients[1]) {
return -a * x * Math.sin(omega * x + phi);
} else {
return -a * Math.sin(omega * x + phi);
}
}
/** Indicator of the need to compute a first guess of the coefficients. */
private boolean firstGuessNeeded;
private static final long serialVersionUID = -8722683066277473450L;
}

44
src/mantissa/src/org/spaceroots/mantissa/fitting/PolynomialCoefficient.java
View File

@ -1,44 +0,0 @@
// 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.spaceroots.mantissa.fitting;
import org.spaceroots.mantissa.estimation.EstimatedParameter;
/** This class represents a polynomial coefficient.
* <p>Each coefficient is uniquely defined by its degree.</p>
* @see PolynomialFitter
* @version $Id$
* @author L. Maisonobe
*/
public class PolynomialCoefficient
extends EstimatedParameter {
public PolynomialCoefficient(int degree) {
super("a" + degree, 0.0);
this.degree = degree;
}
public final int degree;
private static final long serialVersionUID = 5775845068390259552L;
}

107
src/mantissa/src/org/spaceroots/mantissa/fitting/PolynomialFitter.java
View File

@ -1,107 +0,0 @@
// 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.spaceroots.mantissa.fitting;
import org.spaceroots.mantissa.estimation.EstimatedParameter;
import org.spaceroots.mantissa.estimation.Estimator;
import org.spaceroots.mantissa.estimation.GaussNewtonEstimator;
/** This class implements a curve fitting specialized for polynomials.
* <p>Polynomial fitting is a very simple case of curve fitting. The
* estimated coefficients are the polynom coefficients. They are
* searched by a least square estimator.</p>
* <p>This class <emph>is by no means optimized</emph>, neither in
* space nor in time performance.</p>
* @see PolynomialCoefficient
* @version $Id$
* @author L. Maisonobe
*/
public class PolynomialFitter
extends AbstractCurveFitter {
/** Simple constructor.
* <p>The polynomial fitter built this way are complete polynoms,
* ie. a n-degree polynom has n+1 coefficients. In order to build
* fitter for sparse polynoms (for example <code>a x^20 - b
* x^30</code>, on should first build the coefficients array and
* provide it to {@link
* #PolynomialFitter(PolynomialCoefficient[], int, double, double,
* double)}.</p>
* @param degree maximal degree of the polynom
* @param estimator estimator to use for the fitting
*/
public PolynomialFitter(int degree, Estimator estimator) {
super(degree + 1, estimator);
for (int i = 0; i < coefficients.length; ++i) {
coefficients[i] = new PolynomialCoefficient(i);
}
}
/** Simple constructor.
* <p>This constructor can be used either when a first estimate of
* the coefficients is already known (which is of little interest
* because the fit cost is the same whether a first guess is known or
* not) or when one needs to handle sparse polynoms like <code>a
* x^20 - b x^30</code>.</p>
* @param coefficients first estimate of the coefficients.
* A reference to this array is hold by the newly created
* object. Its elements will be adjusted during the fitting process
* and they will be set to the adjusted coefficients at the end.
* @param estimator estimator to use for the fitting
*/
public PolynomialFitter(PolynomialCoefficient[] coefficients,
Estimator estimator) {
super(coefficients, estimator);
}
/** Get the value of the function at x according to the current parameters value.
* @param x abscissa at which the theoretical value is requested
* @return theoretical value at x
*/
public double valueAt(double x) {
double y = coefficients[coefficients.length - 1].getEstimate();
for (int i = coefficients.length - 2; i >= 0; --i) {
y = y * x + coefficients[i].getEstimate();
}
return y;
}
/** Get the derivative of the function at x with respect to parameter p.
* @param x abscissa at which the partial derivative is requested
* @param p parameter with respect to which the derivative is requested
* @return partial derivative
*/
public double partial(double x, EstimatedParameter p) {
if (p instanceof PolynomialCoefficient) {
return Math.pow(x, ((PolynomialCoefficient) p).degree);
}
throw new RuntimeException("internal error");
}
private static final long serialVersionUID = -744904084649890769L;
}

BIN
src/mantissa/src/org/spaceroots/mantissa/fitting/doc-files/org_spaceroots_mantissa_fitting_classes.png
View File

Binary file not shown.

108
src/mantissa/tests-src/org/spaceroots/mantissa/fitting/AbstractCurveFitterTest.java
View File

@ -1,108 +0,0 @@
// 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.spaceroots.mantissa.fitting;
import java.util.Random;
import junit.framework.*;
import org.spaceroots.mantissa.estimation.EstimatedParameter;
import org.spaceroots.mantissa.estimation.LevenbergMarquardtEstimator;
import org.spaceroots.mantissa.estimation.WeightedMeasurement;
public class AbstractCurveFitterTest
extends TestCase {
public AbstractCurveFitterTest(String name) {
super(name);
fitter = null;
}
public void testAlreadySorted() {
for (double x = 0.0; x < 100.0; x += 1.0) {
fitter.addWeightedPair(1.0, x, 0.0);
}
checkSorted();
}
public void testReversed() {
for (double x = 0.0; x < 100.0; x += 1.0) {
fitter.addWeightedPair(1.0, 100.0 - x, 0.0);
}
checkSorted();
}
public void testRandom() {
Random randomizer = new Random(86757343594l);
for (int i = 0; i < 100; ++i) {
fitter.addWeightedPair(1.0, 10.0 * randomizer.nextDouble(), 0.0);
}
checkSorted();
}
public void checkSorted() {
fitter.doSort();
WeightedMeasurement[] measurements = fitter.getMeasurements();
for (int i = 1; i < measurements.length; ++i) {
AbstractCurveFitter.FitMeasurement m1
= (AbstractCurveFitter.FitMeasurement) measurements[i-1];
AbstractCurveFitter.FitMeasurement m2
= (AbstractCurveFitter.FitMeasurement) measurements[i];
assertTrue(m1.x <= m2.x);
}
}
public static Test suite() {
return new TestSuite(AbstractCurveFitterTest.class);
}
public void setUp() {
fitter = new DummyFitter();
}
public void tearDown() {
fitter = null;
}
private static class DummyFitter
extends AbstractCurveFitter {
public DummyFitter() {
super(10, new LevenbergMarquardtEstimator());
}
public double valueAt(double x) {
return 0.0;
}
public double partial(double x, EstimatedParameter p) {
return 0.0;
}
public void doSort() {
sortMeasurements();
}
private static final long serialVersionUID = 4016396219767783678L;
}
private DummyFitter fitter;
}

35
src/mantissa/tests-src/org/spaceroots/mantissa/fitting/AllTests.java
View File

@ -1,35 +0,0 @@
// 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.spaceroots.mantissa.fitting;
import junit.framework.Test;
import junit.framework.TestSuite;
public class AllTests {
public static Test suite() {
TestSuite suite= new TestSuite("org.spaceroots.mantissa.fitting");
suite.addTest(AbstractCurveFitterTest.suite());
suite.addTest(PolynomialFitterTest.suite());
suite.addTest(HarmonicFitterTest.suite());
return suite;
}
}

175
src/mantissa/tests-src/org/spaceroots/mantissa/fitting/HarmonicFitterTest.java
View File

@ -1,175 +0,0 @@
// 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.spaceroots.mantissa.fitting;
import java.util.Random;
import junit.framework.*;
import org.spaceroots.mantissa.estimation.EstimationException;
import org.spaceroots.mantissa.estimation.LevenbergMarquardtEstimator;
import org.spaceroots.mantissa.estimation.WeightedMeasurement;
public class HarmonicFitterTest
extends TestCase {
public HarmonicFitterTest(String name) {
super(name);
}
public void testNoError()
throws EstimationException {
HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtEstimator());
for (double x = 0.0; x < 1.3; x += 0.01) {
fitter.addWeightedPair(1.0, x, f.valueAt(x));
}
double[] coeffs = fitter.fit();
HarmonicFunction fitted = new HarmonicFunction(coeffs[0],
coeffs[1],
coeffs[2]);
assertTrue(Math.abs(coeffs[0] - f.getA()) < 1.0e-13);
assertTrue(Math.abs(coeffs[1] - f.getOmega()) < 1.0e-13);
assertTrue(Math.abs(coeffs[2] - center(f.getPhi(), coeffs[2])) < 1.0e-13);
for (double x = -1.0; x < 1.0; x += 0.01) {
assertTrue(Math.abs(f.valueAt(x) - fitted.valueAt(x)) < 1.0e-13);
}
}
public void test1PercentError()
throws EstimationException {
Random randomizer = new Random(64925784252l);
HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtEstimator());
for (double x = 0.0; x < 10.0; x += 0.1) {
fitter.addWeightedPair(1.0, x,
f.valueAt(x) + 0.01 * randomizer.nextGaussian());
}
double[] coeffs = fitter.fit();
new HarmonicFunction(coeffs[0], coeffs[1], coeffs[2]);
assertTrue(Math.abs(coeffs[0] - f.getA()) < 7.6e-4);
assertTrue(Math.abs(coeffs[1] - f.getOmega()) < 2.7e-3);
assertTrue(Math.abs(coeffs[2] - center(f.getPhi(), coeffs[2])) < 1.3e-2);
WeightedMeasurement[] measurements = fitter.getMeasurements();
for (int i = 0; i < measurements.length; ++i) {
WeightedMeasurement m = measurements[i];
assertTrue(Math.abs(measurements[i].getMeasuredValue()
- m.getTheoreticalValue()) < 0.04);
}
}
public void testUnsorted()
throws EstimationException {
Random randomizer = new Random(64925784252l);
HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtEstimator());
// build a regularly spaced array of measurements
int size = 100;
double[] xTab = new double[size];
double[] yTab = new double[size];
for (int i = 0; i < size; ++i) {
xTab[i] = 0.1 * i;
yTab[i] = f.valueAt (xTab[i]) + 0.01 * randomizer.nextGaussian();
}
// shake it
for (int i = 0; i < size; ++i) {
int i1 = randomizer.nextInt(size);
int i2 = randomizer.nextInt(size);
double xTmp = xTab[i1];
double yTmp = yTab[i1];
xTab[i1] = xTab[i2];
yTab[i1] = yTab[i2];
xTab[i2] = xTmp;
yTab[i2] = yTmp;
}
// pass it to the fitter
for (int i = 0; i < size; ++i) {
fitter.addWeightedPair(1.0, xTab[i], yTab[i]);
}
double[] coeffs = fitter.fit();
new HarmonicFunction(coeffs[0], coeffs[1], coeffs[2]);
assertTrue(Math.abs(coeffs[0] - f.getA()) < 7.6e-4);
assertTrue(Math.abs(coeffs[1] - f.getOmega()) < 3.5e-3);
assertTrue(Math.abs(coeffs[2] - center(f.getPhi(), coeffs[2])) < 1.5e-2);
WeightedMeasurement[] measurements = fitter.getMeasurements();
for (int i = 0; i < measurements.length; ++i) {
WeightedMeasurement m = measurements[i];
assertTrue(Math.abs(m.getMeasuredValue() - m.getTheoreticalValue())
< 0.04);
}
}
public static Test suite() {
return new TestSuite(HarmonicFitterTest.class);
}
/** Center an angle with respect to another one. */
private static double center(double a, double ref) {
double twoPi = Math.PI + Math.PI;
return a - twoPi * Math.floor((a + Math.PI - ref) / twoPi);
}
private static class HarmonicFunction {
public HarmonicFunction(double a, double omega, double phi) {
this.a = a;
this.omega = omega;
this.phi = phi;
}
public double valueAt(double x) {
return a * Math.cos(omega * x + phi);
}
public double getA() {
return a;
}
public double getOmega() {
return omega;
}
public double getPhi() {
return phi;
}
private double a;
private double omega;
private double phi;
}
}

164
src/mantissa/tests-src/org/spaceroots/mantissa/fitting/PolynomialFitterTest.java
View File

@ -1,164 +0,0 @@
// 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.spaceroots.mantissa.fitting;
import java.util.Random;
import junit.framework.*;
import org.spaceroots.mantissa.estimation.EstimationException;
import org.spaceroots.mantissa.estimation.Estimator;
import org.spaceroots.mantissa.estimation.GaussNewtonEstimator;
import org.spaceroots.mantissa.estimation.LevenbergMarquardtEstimator;
public class PolynomialFitterTest
extends TestCase {
public PolynomialFitterTest(String name) {
super(name);
}
public void testNoError()
throws EstimationException {
Random randomizer = new Random(64925784252l);
for (int degree = 0; degree < 10; ++degree) {
Polynom p = new Polynom(degree);
for (int i = 0; i <= degree; ++i) {
p.initCoeff (i, randomizer.nextGaussian());
}
PolynomialFitter fitter =
new PolynomialFitter(degree, new LevenbergMarquardtEstimator());
for (int i = 0; i <= degree; ++i) {
fitter.addWeightedPair(1.0, i, p.valueAt(i));
}
Polynom fitted = new Polynom(fitter.fit());
for (double x = -1.0; x < 1.0; x += 0.01) {
double error = Math.abs(p.valueAt(x) - fitted.valueAt(x))
/ (1.0 + Math.abs(p.valueAt(x)));
assertTrue(Math.abs(error) < 1.0e-5);
}
}
}
public void testSmallError()
throws EstimationException {
Random randomizer = new Random(53882150042l);
for (int degree = 0; degree < 10; ++degree) {
Polynom p = new Polynom(degree);
for (int i = 0; i <= degree; ++i) {
p.initCoeff(i, randomizer.nextGaussian());
}
PolynomialFitter fitter =
new PolynomialFitter(degree, new LevenbergMarquardtEstimator());
for (double x = -1.0; x < 1.0; x += 0.01) {
fitter.addWeightedPair(1.0, x,
p.valueAt(x) + 0.1 * randomizer.nextGaussian());
}
Polynom fitted = new Polynom(fitter.fit());
for (double x = -1.0; x < 1.0; x += 0.01) {
double error = Math.abs(p.valueAt(x) - fitted.valueAt(x))
/ (1.0 + Math.abs(p.valueAt(x)));
assertTrue(Math.abs(error) < 0.1);
}
}
}
public void testRedundantSolvable() {
// Levenberg-Marquardt should handle redundant information gracefully
checkUnsolvableProblem(new LevenbergMarquardtEstimator(), true);
}
public void testRedundantUnsolvable() {
// Gauss-Newton should not be able to solve redundant information
checkUnsolvableProblem(new GaussNewtonEstimator(10, 1.0e-7, 1.0e-7,
1.0e-10),
false);
}
private void checkUnsolvableProblem(Estimator estimator,
boolean solvable) {
Random randomizer = new Random(1248788532l);
for (int degree = 0; degree < 10; ++degree) {
Polynom p = new Polynom(degree);
for (int i = 0; i <= degree; ++i) {
p.initCoeff(i, randomizer.nextGaussian());
}
PolynomialFitter fitter = new PolynomialFitter(degree, estimator);
// reusing the same point over and over again does not bring
// information, the problem cannot be solved in this case for
// degrees greater than 1 (but one point is sufficient for
// degree 0)
for (double x = -1.0; x < 1.0; x += 0.01) {
fitter.addWeightedPair(1.0, 0.0, p.valueAt(0.0));
}
try {
fitter.fit();
assertTrue(solvable || (degree == 0));
} catch(EstimationException e) {
assertTrue((! solvable) && (degree > 0));
}
}
}
public static Test suite() {
return new TestSuite(PolynomialFitterTest.class);
}
private static class Polynom {
public Polynom(int degree) {
coeffs = new double[degree + 1];
for (int i = 0; i < coeffs.length; ++i) {
coeffs[i] = 0.0;
}
}
public Polynom(double[]coeffs) {
this.coeffs = coeffs;
}
public void initCoeff(int i, double c) {
coeffs[i] = c;
}
public double valueAt(double x) {
double y = coeffs[coeffs.length - 1];
for (int i = coeffs.length - 2; i >= 0; --i) {
y = y * x + coeffs[i];
}
return y;
}
private double[] coeffs;
}
}

3
src/site/xdoc/changes.xml
View File

@ -39,6 +39,9 @@ The <action> type attribute can be add,update,fix,remove.
</properties>
<body>
<release version="2.0" date="TBD" description="TBD">
<action dev="luc" type="add" >
Added curve fitting with a general case and two specific cases (polynomial and harmonic).
</action>
<action dev="psteitz" type="update" issue="MATH-276" due-to="Mark Anderson">
Optimized Complex isNaN(), isInfinite() by moving computation to constructor.
</action>

1
src/site/xdoc/userguide/index.xml
View File

@ -120,6 +120,7 @@
<li><a href="analysis.html#a12.3_Linear_Programming">12.3 Linear Programming</a></li>
<li><a href="optimization.html#a12.4_Direct_Methods">12.4 Direct Methods</a></li>
<li><a href="analysis.html#a12.5_General_Case">12.5 General Case</a></li>
<li><a href="analysis.html#a12.6_Curve_Fitting">12.6 Curve Fitting</a></li>
</ul></li>
<li><a href="ode.html">13. Ordinary Differential Equations Integration</a>
<ul>

50
src/site/xdoc/userguide/optimization.xml
View File

@ -40,6 +40,7 @@
using direct search methods (i.e. not using derivatives),</li>
<li>the general package handles multivariate scalar or vector functions
using derivatives.</li>
<li>the fitting package handles curve fitting by univariate real functions</li>
</ul>
</p>
<p>
@ -212,6 +213,55 @@
solver).
</p>
</subsection>
<subsection name="12.6 Curve Fitting" href="fitting">
<p>
The fitting package deals with curve fitting for univariate real functions.
When a univariate real function y = f(x) does depend on some unknown parameters
p<sub>0</sub>, p<sub>1</sub> ... p<sub>n-1</sub>, curve fitting can be used to
find these parameters. It does this by <em>fitting</em> the curve so it remains
very close to a set of observed points (x<sub>0</sub>, y<sub>0</sub>),
(x<sub>1</sub>, y<sub>1</sub>) ... (x<sub>k-1</sub>, y<sub>k-1</sub>). This
fitting is done by finding the parameters values that minimizes the objective
function sum(y<sub>i</sub>-f(x<sub>i</sub>))<sup>2</sup>. This is really a least
squares problem.
</p>
<p>
For all provided curve fitters, the operating principle is the same. Users must first
create an instance of the fitter, then add the observed points and once the complete
sample of observed points has been added they must call the <code>fit</code> method
which will compute the parameters that best fit the sample. A weight is associated
with each observed point, this allows to take into account uncertainty on some points
when they come from loosy measurements for example. If no such information exist and
all points should be treated the same, it is safe to put 1.0 as the weight for all points.
</p>
<p>
The <a
href="../apidocs/org/apache/commons/math/optimization/fitting/CurveFitter.html">
CurveFitter</a> class provides curve fitting for general curves. Users must
provide their own implementation of the curve template as a class implementing
the <a
href="../apidocs/org/apache/commons/math/optimization/fitting/ParametricRealFunction.html">
ParametricRealFunction</a> interface and they must provide the initial guess of the
parameters. The more specialized <a
href="../apidocs/org/apache/commons/math/optimization/fitting/PolynomialFitter.html">
PolynomialFitter</a> and <a
href="../apidocs/org/apache/commons/math/optimization/fitting/HarmonicFitter.html">
HarmonicFitter</a> classes require neither an implementation of the parametric real function
not an initial guess as they are able to compute them by themselves.
</p>
<p>
An example of fitting a polynomial is given here:
</p>
<source>PolynomialFitter fitter = new PolynomialFitter(degree, new LevenbergMarquardtOptimizer());
fitter.addObservedPoint(-1.00, 2.021170021833143);
fitter.addObservedPoint(-0.99 2.221135431136975);
fitter.addObservedPoint(-0.98 2.09985277659314);
fitter.addObservedPoint(-0.97 2.0211192647627025);
// lots of lines ommitted
fitter.addObservedPoint( 0.99, -2.4345814727089854);
PolynomialFunction fitted = fitter.fit();
</source>
</subsection>
</section>
</body>
</document>

132
src/test/org/apache/commons/math/optimization/fitting/HarmonicFitterTest.java Normal file
View File

@ -0,0 +1,132 @@
// 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 static org.junit.Assert.assertTrue;
import java.util.Random;
import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
import org.apache.commons.math.util.MathUtils;
import org.junit.Test;
public class HarmonicFitterTest {
@Test
public void testNoError() throws OptimizationException {
HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtOptimizer());
for (double x = 0.0; x < 1.3; x += 0.01) {
fitter.addObservedPoint(1.0, x, f.value(x));
}
HarmonicFunction fitted = fitter.fit();
assertEquals(f.getAmplitude(), fitted.getAmplitude(), 1.0e-13);
assertEquals(f.getPulsation(), fitted.getPulsation(), 1.0e-13);
assertEquals(f.getPhase(), MathUtils.normalizeAngle(fitted.getPhase(), f.getPhase()), 1.0e-13);
for (double x = -1.0; x < 1.0; x += 0.01) {
assertTrue(Math.abs(f.value(x) - fitted.value(x)) < 1.0e-13);
}
}
@Test
public void test1PercentError() throws OptimizationException {
Random randomizer = new Random(64925784252l);
HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtOptimizer());
for (double x = 0.0; x < 10.0; x += 0.1) {
fitter.addObservedPoint(1.0, x,
f.value(x) + 0.01 * randomizer.nextGaussian());
}
HarmonicFunction fitted = fitter.fit();
assertEquals(f.getAmplitude(), fitted.getAmplitude(), 7.6e-4);
assertEquals(f.getPulsation(), fitted.getPulsation(), 2.7e-3);
assertEquals(f.getPhase(), MathUtils.normalizeAngle(fitted.getPhase(), f.getPhase()), 1.3e-2);
}
@Test
public void testInitialGuess() throws OptimizationException {
Random randomizer = new Random(45314242l);
HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtOptimizer(), new double[] { 0.15, 3.6, 4.5 });
for (double x = 0.0; x < 10.0; x += 0.1) {
fitter.addObservedPoint(1.0, x,
f.value(x) + 0.01 * randomizer.nextGaussian());
}
HarmonicFunction fitted = fitter.fit();
assertEquals(f.getAmplitude(), fitted.getAmplitude(), 1.2e-3);
assertEquals(f.getPulsation(), fitted.getPulsation(), 3.3e-3);
assertEquals(f.getPhase(), MathUtils.normalizeAngle(fitted.getPhase(), f.getPhase()), 1.7e-2);
}
@Test
public void testUnsorted() throws OptimizationException {
Random randomizer = new Random(64925784252l);
HarmonicFunction f = new HarmonicFunction(0.2, 3.4, 4.1);
HarmonicFitter fitter =
new HarmonicFitter(new LevenbergMarquardtOptimizer());
// build a regularly spaced array of measurements
int size = 100;
double[] xTab = new double[size];
double[] yTab = new double[size];
for (int i = 0; i < size; ++i) {
xTab[i] = 0.1 * i;
yTab[i] = f.value(xTab[i]) + 0.01 * randomizer.nextGaussian();
}
// shake it
for (int i = 0; i < size; ++i) {
int i1 = randomizer.nextInt(size);
int i2 = randomizer.nextInt(size);
double xTmp = xTab[i1];
double yTmp = yTab[i1];
xTab[i1] = xTab[i2];
yTab[i1] = yTab[i2];
xTab[i2] = xTmp;
yTab[i2] = yTmp;
}
// pass it to the fitter
for (int i = 0; i < size; ++i) {
fitter.addObservedPoint(1.0, xTab[i], yTab[i]);
}
HarmonicFunction fitted = fitter.fit();
assertEquals(f.getAmplitude(), fitted.getAmplitude(), 7.6e-4);
assertEquals(f.getPulsation(), fitted.getPulsation(), 3.5e-3);
assertEquals(f.getPhase(), MathUtils.normalizeAngle(fitted.getPhase(), f.getPhase()), 1.5e-2);
}
}

134
src/test/org/apache/commons/math/optimization/fitting/PolynomialFitterTest.java Normal file
View File

@ -0,0 +1,134 @@
// 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 static org.junit.Assert.assertTrue;
import java.util.Random;
import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.general.GaussNewtonOptimizer;
import org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizer;
import org.junit.Test;
public class PolynomialFitterTest {
@Test
public void testNoError() throws OptimizationException {
Random randomizer = new Random(64925784252l);
for (int degree = 1; degree < 10; ++degree) {
PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
PolynomialFitter fitter =
new PolynomialFitter(degree, new LevenbergMarquardtOptimizer());
for (int i = 0; i <= degree; ++i) {
fitter.addObservedPoint(1.0, i, p.value(i));
}
PolynomialFunction fitted = fitter.fit();
for (double x = -1.0; x < 1.0; x += 0.01) {
double error = Math.abs(p.value(x) - fitted.value(x)) /
(1.0 + Math.abs(p.value(x)));
assertEquals(0.0, error, 1.0e-6);
}
}
}
@Test
public void testSmallError() throws OptimizationException {
Random randomizer = new Random(53882150042l);
double maxError = 0;
for (int degree = 0; degree < 10; ++degree) {
PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
PolynomialFitter fitter =
new PolynomialFitter(degree, new LevenbergMarquardtOptimizer());
for (double x = -1.0; x < 1.0; x += 0.01) {
fitter.addObservedPoint(1.0, x,
p.value(x) + 0.1 * randomizer.nextGaussian());
}
PolynomialFunction fitted = fitter.fit();
for (double x = -1.0; x < 1.0; x += 0.01) {
double error = Math.abs(p.value(x) - fitted.value(x)) /
(1.0 + Math.abs(p.value(x)));
maxError = Math.max(maxError, error);
assertTrue(Math.abs(error) < 0.1);
}
}
assertTrue(maxError > 0.01);
}
@Test
public void testRedundantSolvable() {
// Levenberg-Marquardt should handle redundant information gracefully
checkUnsolvableProblem(new LevenbergMarquardtOptimizer(), true);
}
@Test
public void testRedundantUnsolvable() {
// Gauss-Newton should not be able to solve redundant information
DifferentiableMultivariateVectorialOptimizer optimizer =
new GaussNewtonOptimizer(true);
checkUnsolvableProblem(optimizer, false);
}
private void checkUnsolvableProblem(DifferentiableMultivariateVectorialOptimizer optimizer,
boolean solvable) {
Random randomizer = new Random(1248788532l);
for (int degree = 0; degree < 10; ++degree) {
PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
PolynomialFitter fitter = new PolynomialFitter(degree, optimizer);
// reusing the same point over and over again does not bring
// information, the problem cannot be solved in this case for
// degrees greater than 1 (but one point is sufficient for
// degree 0)
for (double x = -1.0; x < 1.0; x += 0.01) {
fitter.addObservedPoint(1.0, 0.0, p.value(0.0));
}
try {
fitter.fit();
assertTrue(solvable || (degree == 0));
} catch(OptimizationException e) {
assertTrue((! solvable) && (degree > 0));
}
}
}
private PolynomialFunction buildRandomPolynomial(int degree, Random randomizer) {
final double[] coefficients = new double[degree + 1];
for (int i = 0; i <= degree; ++i) {
coefficients[i] = randomizer.nextGaussian();
}
return new PolynomialFunction(coefficients);
}
}