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

added the estimation package from Mantissa 

git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@512061 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2007-02-26 22:59:45 +00:00
parent 7f8c5e2562
commit f3b02ccea3
16 changed files with 87 additions and 446 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/mantissa/src/org/spaceroots/mantissa/estimation/EstimatedParameter.java → src/java/org/apache/commons/math/estimation/EstimatedParameter.java
View File

@ -15,7 +15,7 @@
// specific language governing permissions and limitations
// under the License.
package org.spaceroots.mantissa.estimation;
package org.apache.commons.math.estimation;
import java.io.Serializable;

17
src/mantissa/src/org/spaceroots/mantissa/estimation/EstimationException.java → src/java/org/apache/commons/math/estimation/EstimationException.java
View File

@ -15,9 +15,9 @@
// specific language governing permissions and limitations
// under the License.
package org.spaceroots.mantissa.estimation;
package org.apache.commons.math.estimation;
import org.spaceroots.mantissa.MantissaException;
import org.apache.commons.math.MathException;
/** This class represents exceptions thrown by the estimation solvers.
@ -27,15 +27,10 @@ import org.spaceroots.mantissa.MantissaException;
*/
public class EstimationException
extends MantissaException {
extends MathException {
/** Simple constructor.
* Build an exception by translating the specified message
* @param message message to translate
*/
public EstimationException(String message) {
super(message);
}
/** Serializable version identifier. */
private static final long serialVersionUID = -7414806622114810487L;
/** Simple constructor.
* Build an exception by translating and formating a message
@ -54,6 +49,4 @@ public class EstimationException
super(cause);
}
private static final long serialVersionUID = 1613719630569355278L;
}

2
src/mantissa/src/org/spaceroots/mantissa/estimation/EstimationProblem.java → src/java/org/apache/commons/math/estimation/EstimationProblem.java
View File

@ -15,7 +15,7 @@
// specific language governing permissions and limitations
// under the License.
package org.spaceroots.mantissa.estimation;
package org.apache.commons.math.estimation;
/** This interface represents an estimation problem.

2
src/mantissa/src/org/spaceroots/mantissa/estimation/Estimator.java → src/java/org/apache/commons/math/estimation/Estimator.java
View File

@ -15,7 +15,7 @@
// specific language governing permissions and limitations
// under the License.
package org.spaceroots.mantissa.estimation;
package org.apache.commons.math.estimation;
/** This interface represents solvers for estimation problems.

53
src/mantissa/src/org/spaceroots/mantissa/estimation/GaussNewtonEstimator.java → src/java/org/apache/commons/math/estimation/GaussNewtonEstimator.java
View File

@ -15,14 +15,13 @@
// specific language governing permissions and limitations
// under the License.
package org.spaceroots.mantissa.estimation;
package org.apache.commons.math.estimation;
import java.io.Serializable;
import org.spaceroots.mantissa.linalg.Matrix;
import org.spaceroots.mantissa.linalg.GeneralMatrix;
import org.spaceroots.mantissa.linalg.SymetricalMatrix;
import org.spaceroots.mantissa.linalg.SingularMatrixException;
import org.apache.commons.math.linear.InvalidMatrixException;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.RealMatrixImpl;
/** This class implements a solver for estimation problems.
@ -67,18 +66,13 @@ public class GaussNewtonEstimator
* <code>Jn</code> and <code>Jn-1</code> are the current and
* preceding criterion value (square sum of the weighted residuals
* of considered measurements).
* @param epsilon threshold under which the matrix of the linearized
* problem is considered singular (see {@link
* org.spaceroots.mantissa.linalg.SquareMatrix#solve(Matrix,double)
* SquareMatrix.solve}). */
*/
public GaussNewtonEstimator(int maxIterations,
double convergence,
double steadyStateThreshold,
double epsilon) {
double steadyStateThreshold) {
this.maxIterations = maxIterations;
this.steadyStateThreshold = steadyStateThreshold;
this.convergence = convergence;
this.epsilon = epsilon;
}
/** Solve an estimation problem using a least squares criterion.
@ -153,24 +147,37 @@ public class GaussNewtonEstimator
WeightedMeasurement[] measurements = problem.getMeasurements();
// build the linear problem
GeneralMatrix b = new GeneralMatrix(parameters.length, 1);
SymetricalMatrix a = new SymetricalMatrix(parameters.length);
RealMatrix b = new RealMatrixImpl(parameters.length, 1);
RealMatrix a = new RealMatrixImpl(parameters.length, parameters.length);
double[] grad = new double[parameters.length];
RealMatrixImpl bDecrement = new RealMatrixImpl(parameters.length, 1);
double[][] bDecrementData = bDecrement.getDataRef();
RealMatrixImpl wGradGradT = new RealMatrixImpl(parameters.length, parameters.length);
double[][] wggData = wGradGradT.getDataRef();
for (int i = 0; i < measurements.length; ++i) {
if (! measurements [i].isIgnored()) {
double weight = measurements[i].getWeight();
double residual = measurements[i].getResidual();
// compute the normal equation
double[] grad = new double[parameters.length];
Matrix bDecrement = new GeneralMatrix(parameters.length, 1);
for (int j = 0; j < parameters.length; ++j) {
grad[j] = measurements[i].getPartial(parameters[j]);
bDecrement.setElement(j, 0, weight * residual * grad[j]);
bDecrementData[j][0] = weight * residual * grad[j];
}
// build the contribution matrix for measurement i
for (int k = 0; k < parameters.length; ++k) {
double[] wggRow = wggData[k];
double gk = grad[k];
for (int l = 0; l < parameters.length; ++l) {
wggRow[l] = weight * gk * grad[l];
}
}
// update the matrices
a.selfAddWAAt(weight, grad);
b.selfAdd(bDecrement);
a = a.add(wGradGradT);
b = b.add(bDecrement);
}
}
@ -178,15 +185,14 @@ public class GaussNewtonEstimator
try {
// solve the linearized least squares problem
Matrix dX = a.solve(b, epsilon);
RealMatrix dX = a.solve(b);
// update the estimated parameters
for (int i = 0; i < parameters.length; ++i) {
parameters[i].setEstimate(parameters[i].getEstimate()
+ dX.getElement(i, 0));
parameters[i].setEstimate(parameters[i].getEstimate() + dX.getEntry(i, 0));
}
} catch(SingularMatrixException e) {
} catch(InvalidMatrixException e) {
throw new EstimationException(e);
}
@ -223,7 +229,6 @@ public class GaussNewtonEstimator
private int maxIterations;
private double steadyStateThreshold;
private double convergence;
private double epsilon;
private static final long serialVersionUID = -7606628156644194170L;

2
src/mantissa/src/org/spaceroots/mantissa/estimation/LevenbergMarquardtEstimator.java → src/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java
View File

@ -14,7 +14,7 @@
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.spaceroots.mantissa.estimation;
package org.apache.commons.math.estimation;
import java.io.Serializable;
import java.util.Arrays;

2
src/mantissa/src/org/spaceroots/mantissa/estimation/WeightedMeasurement.java → src/java/org/apache/commons/math/estimation/WeightedMeasurement.java
View File

@ -15,7 +15,7 @@
// specific language governing permissions and limitations
// under the License.
package org.spaceroots.mantissa.estimation;
package org.apache.commons.math.estimation;
import java.io.Serializable;

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

73
src/mantissa/src/org/spaceroots/mantissa/estimation/LeastSquaresEstimator.java
View File

@ -1,73 +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.estimation;
import java.io.Serializable;
/** This class implements a solver for estimation problems.
* @deprecated this class has been replaced by the {@link
* org.spaceroots.mantissa.estimation.GaussNewtonEstimator GaussNewtonEstimator}
* class. It is now a simple wrapper delegating everything to {@link
* org.spaceroots.mantissa.estimation.GaussNewtonEstimator GaussNewtonEstimator}
* @version $Id: LeastSquaresEstimator.java 1705 2006-09-17 19:57:39Z luc $
* @author L. Maisonobe
*/
public class LeastSquaresEstimator implements Estimator, Serializable {
/** Simple constructor.
* @see org.spaceroots.mantissa.estimation.GaussNewtonEstimator#GaussNewtonEstimator(int,
* double, double, double)
*/
public LeastSquaresEstimator(int maxIterations,
double convergence,
double steadyStateThreshold,
double epsilon) {
estimator = new GaussNewtonEstimator(maxIterations,
convergence,
steadyStateThreshold,
epsilon);
}
/** Solve an estimation problem using a least squares criterion.
* @see org.spaceroots.mantissa.estimation.GaussNewtonEstimator#estimate
*/
public void estimate(EstimationProblem problem)
throws EstimationException {
estimator.estimate(problem);
}
/** Estimate the solution of a linear least square problem.
* @see org.spaceroots.mantissa.estimation.GaussNewtonEstimator#linearEstimate
*/
public void linearEstimate(EstimationProblem problem)
throws EstimationException {
estimator.linearEstimate(problem);
}
/** Get the Root Mean Square value.
* @see org.spaceroots.mantissa.estimation.GaussNewtonEstimator#getRMS
*/
public double getRMS(EstimationProblem problem) {
return estimator.getRMS(problem);
}
private GaussNewtonEstimator estimator;
private static final long serialVersionUID = -7542643494637247770L;
}

BIN
src/mantissa/src/org/spaceroots/mantissa/estimation/doc-files/org_spaceroots_mantissa_estimation_classes.png
View File

Binary file not shown.

39
src/mantissa/tests-src/org/spaceroots/mantissa/estimation/AllTests.java
View File

@ -1,39 +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.estimation;
import junit.framework.Test;
import junit.framework.TestSuite;
public class AllTests {
public static Test suite() {
TestSuite suite = new TestSuite("org.spaceroots.mantissa.estimation");
suite.addTest(EstimatedParameterTest.suite());
suite.addTest(WeightedMeasurementTest.suite());
suite.addTest(GaussNewtonEstimatorTest.suite());
suite.addTest(LevenbergMarquardtEstimatorTest.suite());
suite.addTest(MinpackTest.suite());
return suite;
}
}

262
src/mantissa/tests-src/org/spaceroots/mantissa/estimation/GaussNewtonEstimatorTest.java
View File

@ -1,262 +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.estimation;
import java.util.Random;
import junit.framework.*;
public class GaussNewtonEstimatorTest
extends TestCase
implements EstimationProblem {
public GaussNewtonEstimatorTest(String name) {
super(name);
}
public void testNoMeasurementError()
throws EstimationException {
initRandomizedGrid(2.3);
initProblem(0.0);
GaussNewtonEstimator estimator =
new GaussNewtonEstimator(100, 1.0e-7, 1.0e-10, 1.0e-10);
estimator.estimate(this);
checkGrid(0.01);
}
public void testSmallMeasurementError()
throws EstimationException {
initRandomizedGrid(2.3);
initProblem(0.02);
GaussNewtonEstimator estimator =
new GaussNewtonEstimator(100, 1.0e-7, 1.0e-10, 1.0e-10);
estimator.estimate(this);
checkGrid(0.1);
}
public void testNoError()
throws EstimationException {
initRandomizedGrid(0.0);
initProblem(0.0);
GaussNewtonEstimator estimator =
new GaussNewtonEstimator(100, 1.0e-7, 1.0e-10, 1.0e-10);
estimator.estimate(this);
checkGrid(1.0e-10);
}
public void testUnsolvableProblem() {
initRandomizedGrid(2.3);
initProblem(0.0);
// reduce the number of measurements below the limit threshold
int unknowns = unboundPars.length;
WeightedMeasurement[] reducedSet = new WeightedMeasurement[unknowns - 1];
for (int i = 0; i < reducedSet.length; ++i) {
reducedSet[i] = measurements[i];
}
measurements = reducedSet;
boolean gotIt = false;
try {
GaussNewtonEstimator estimator =
new GaussNewtonEstimator(100, 1.0e-7, 1.0e-10, 1.0e-10);
estimator.estimate(this);
} catch(EstimationException e) {
gotIt = true;
}
assertTrue(gotIt);
}
public static Test suite() {
return new TestSuite(GaussNewtonEstimatorTest.class);
}
public void setUp() {
initPerfectGrid(5);
}
public void tearDown() {
perfectPars = null;
randomizedPars = null;
unboundPars = null;
measurements = null;
}
private void initPerfectGrid(int gridSize) {
perfectPars = new EstimatedParameter[gridSize * gridSize * 2];
int k = 0;
for (int i = 0; i < gridSize; ++i) {
for (int j = 0; j < gridSize; ++j) {
String name = Integer.toString(k);
perfectPars[2 * k] = new EstimatedParameter("x" + name, i);
perfectPars[2 * k + 1] = new EstimatedParameter("y" + name, j);
++k;
}
}
}
private void initRandomizedGrid(double initialGuessError) {
Random randomizer = new Random(2353995334l);
randomizedPars = new EstimatedParameter[perfectPars.length];
// add an error to every point coordinate
for (int k = 0; k < randomizedPars.length; ++k) {
String name = perfectPars[k].getName();
double value = perfectPars[k].getEstimate();
double error = randomizer.nextGaussian() * initialGuessError;
randomizedPars[k] = new EstimatedParameter(name, value + error);
}
}
private void initProblem(double measurementError) {
int pointsNumber = randomizedPars.length / 2;
int measurementsNumber = pointsNumber * (pointsNumber - 1) / 2;
measurements = new WeightedMeasurement[measurementsNumber];
Random randomizer = new Random(5785631926l);
// for the test, we consider that the perfect grid is the reality
// and that the randomized grid is the first (wrong) estimate.
int i = 0;
for (int l = 0; l < (pointsNumber - 1); ++l) {
for (int m = l + 1; m < pointsNumber; ++m) {
// perfect measurements on the real data
double dx = perfectPars[2 * l].getEstimate()
- perfectPars[2 * m].getEstimate();
double dy = perfectPars[2 * l + 1].getEstimate()
- perfectPars[2 * m + 1].getEstimate();
double d = Math.sqrt(dx * dx + dy * dy);
// adding a noise to the measurements
d += randomizer.nextGaussian() * measurementError;
// add the measurement to the current problem
measurements[i++] = new Distance(1.0, d,
randomizedPars[2 * l],
randomizedPars[2 * l + 1],
randomizedPars[2 * m],
randomizedPars[2 * m + 1]);
}
}
// fix three values in the randomized grid and bind them (there
// are two abscissas and one ordinate, so if there were no error
// at all, the estimated grid should be correctly centered on the
// perfect grid)
int oddNumber = 2 * (randomizedPars.length / 4) - 1;
for (int k = 0; k < 2 * oddNumber + 1; k += oddNumber) {
randomizedPars[k].setEstimate(perfectPars[k].getEstimate());
randomizedPars[k].setBound(true);
}
// store the unbound parameters in a specific table
unboundPars = new EstimatedParameter[randomizedPars.length - 3];
for (int src = 0, dst = 0; src < randomizedPars.length; ++src) {
if (! randomizedPars[src].isBound()) {
unboundPars[dst++] = randomizedPars[src];
}
}
}
private void checkGrid(double threshold) {
double rms = 0;
for (int i = 0; i < perfectPars.length; ++i) {
rms += perfectPars[i].getEstimate() - randomizedPars[i].getEstimate();
}
rms = Math.sqrt(rms / perfectPars.length);
assertTrue(rms <= threshold);
}
private static class Distance extends WeightedMeasurement {
public Distance(double weight, double measuredValue,
EstimatedParameter x1, EstimatedParameter y1,
EstimatedParameter x2, EstimatedParameter y2) {
super(weight, measuredValue);
this.x1 = x1;
this.y1 = y1;
this.x2 = x2;
this.y2 = y2;
}
public double getTheoreticalValue() {
double dx = x2.getEstimate() - x1.getEstimate();
double dy = y2.getEstimate() - y1.getEstimate();
return Math.sqrt(dx * dx + dy * dy);
}
public double getPartial(EstimatedParameter p) {
// first quick answer for most parameters
if ((p != x1) && (p != y1) && (p != x2) && (p != y2)) {
return 0.0;
}
// compute the value now as we know we depend on the specified parameter
double distance = getTheoreticalValue();
if (p == x1) {
return (x1.getEstimate() - x2.getEstimate()) / distance;
} else if (p == x2) {
return (x2.getEstimate() - x1.getEstimate()) / distance;
} else if (p == y1) {
return (y1.getEstimate() - y2.getEstimate()) / distance;
} else {
return (y2.getEstimate() - y1.getEstimate()) / distance;
}
}
private EstimatedParameter x1;
private EstimatedParameter y1;
private EstimatedParameter x2;
private EstimatedParameter y2;
private static final long serialVersionUID = 4090004243280980746L;
}
public WeightedMeasurement[] getMeasurements() {
return (WeightedMeasurement[]) measurements.clone();
}
public EstimatedParameter[] getUnboundParameters() {
return (EstimatedParameter[]) unboundPars.clone();
}
public EstimatedParameter[] getAllParameters() {
return (EstimatedParameter[]) randomizedPars.clone();
}
private EstimatedParameter[] perfectPars;
private EstimatedParameter[] randomizedPars;
private EstimatedParameter[] unboundPars;
private WeightedMeasurement[] measurements;
}

4
src/mantissa/tests-src/org/spaceroots/mantissa/estimation/EstimatedParameterTest.java → src/test/org/apache/commons/math/estimation/EstimatedParameterTest.java
View File

@ -15,7 +15,9 @@
// specific language governing permissions and limitations
// under the License.
package org.spaceroots.mantissa.estimation;
package org.apache.commons.math.estimation;
import org.apache.commons.math.estimation.EstimatedParameter;
import junit.framework.*;

12
src/mantissa/tests-src/org/spaceroots/mantissa/estimation/LevenbergMarquardtEstimatorTest.java → src/test/org/apache/commons/math/estimation/LevenbergMarquardtEstimatorTest.java
View File

@ -15,13 +15,19 @@
* limitations under the License.
*/
package org.spaceroots.mantissa.estimation;
package org.apache.commons.math.estimation;
import java.util.ArrayList;
import java.util.IdentityHashMap;
import java.util.HashMap;
import java.util.Iterator;
import java.util.Set;
import org.apache.commons.math.estimation.EstimatedParameter;
import org.apache.commons.math.estimation.EstimationException;
import org.apache.commons.math.estimation.EstimationProblem;
import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
import org.apache.commons.math.estimation.WeightedMeasurement;
import junit.framework.*;
/**
@ -519,7 +525,7 @@ public class LevenbergMarquardtEstimatorTest
}
public EstimatedParameter[] getAllParameters() {
IdentityHashMap map = new IdentityHashMap();
HashMap map = new HashMap();
for (int i = 0; i < measurements.length; ++i) {
EstimatedParameter[] parameters = measurements[i].getParameters();
for (int j = 0; j < parameters.length; ++j) {

8
src/mantissa/tests-src/org/spaceroots/mantissa/estimation/MinpackTest.java → src/test/org/apache/commons/math/estimation/MinpackTest.java
View File

@ -1,7 +1,13 @@
package org.spaceroots.mantissa.estimation;
package org.apache.commons.math.estimation;
import java.util.Arrays;
import org.apache.commons.math.estimation.EstimatedParameter;
import org.apache.commons.math.estimation.EstimationException;
import org.apache.commons.math.estimation.EstimationProblem;
import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
import org.apache.commons.math.estimation.WeightedMeasurement;
import junit.framework.*;
/**

5
src/mantissa/tests-src/org/spaceroots/mantissa/estimation/WeightedMeasurementTest.java → src/test/org/apache/commons/math/estimation/WeightedMeasurementTest.java
View File

@ -15,7 +15,10 @@
// specific language governing permissions and limitations
// under the License.
package org.spaceroots.mantissa.estimation;
package org.apache.commons.math.estimation;
import org.apache.commons.math.estimation.EstimatedParameter;
import org.apache.commons.math.estimation.WeightedMeasurement;
import junit.framework.*;