Apache/commons-math
1
0
Fork 0
mirror of https://github.com/apache/commons-math.git synced 2025-03-08 17:39:14 +00:00
Code Issues Packages Projects Releases Wiki Activity

merged curve fitting from mantissa into commons-math

Browse Source 
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk/src/mantissa@786479 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-06-19 12:36:16 +00:00
parent c6361017d4
commit 9c152a81cf
13 changed files with 0 additions and 1493 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

222
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/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/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/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/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/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/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/org/spaceroots/mantissa/fitting/doc-files/org_spaceroots_mantissa_fitting_classes.png
View File

Binary file not shown.

26
src/org/spaceroots/mantissa/fitting/package.html
View File

@ -1,26 +0,0 @@
<html>
<body>
This package provides classes to perform curve fitting.
<p>Curve fitting is a special case of an {@link
org.spaceroots.mantissa.estimation.EstimationProblem estimation problem}
were the parameters are the coefficients of a function <code>f</code>
whose graph <code>y=f(x)</code> should pass through sample points, and
were the measurements are the ordinates of the curve
<code>yi=f(xi)</code>.</p>
<p>The organisation of this package is explained in the class diagram
below. The {@link org.spaceroots.mantissa.fitting.AbstractCurveFitter
AbstractCurveFitter} class is the base class for all curve fitting
classes, it handles all common features like sample points
handling. Each specific curve fitting class should sub-class this
abstract class and implement the {@link
org.spaceroots.mantissa.fitting.AbstractCurveFitter#valueAt valueAt} and
{@link org.spaceroots.mantissa.fitting.AbstractCurveFitter#partial
partial} methods.</p>
<img src="doc-files/org_spaceroots_mantissa_fitting_classes.png" />
@author L. Maisonobe
</body>
</html>

108
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
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
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
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;
}
}