/*
* 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.math3.optimization.general;
import java.io.Serializable;
import java.util.ArrayList;
import java.util.List;
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
import org.apache.commons.math3.analysis.differentiation.MultivariateDifferentiableVectorFunction;
import org.apache.commons.math3.exception.ConvergenceException;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.TooManyEvaluationsException;
import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
import org.apache.commons.math3.linear.SingularMatrixException;
import org.apache.commons.math3.optimization.PointVectorValuePair;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.Precision;
import org.junit.Assert;
import org.junit.Test;
import org.junit.Ignore;
/**
* <p>Some of the unit tests are re-implementations of the MINPACK <a
* href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a
* href="http://www.netlib.org/minpack/ex/file22">file22</a> test files.
* The redistribution policy for MINPACK is available <a
* href="http://www.netlib.org/minpack/disclaimer">here</a>, for
* convenience, it is reproduced below.</p>
* <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
* <tr><td>
* Minpack Copyright Notice (1999) University of Chicago.
* All rights reserved
* </td></tr>
* <tr><td>
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* <ol>
* <li>Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.</li>
* <li>Redistributions in binary form must reproduce the above
* copyright notice, this list of conditions and the following
* disclaimer in the documentation and/or other materials provided
* with the distribution.</li>
* <li>The end-user documentation included with the redistribution, if any,
* must include the following acknowledgment:
* <code>This product includes software developed by the University of
* Chicago, as Operator of Argonne National Laboratory.</code>
* Alternately, this acknowledgment may appear in the software itself,
* if and wherever such third-party acknowledgments normally appear.</li>
* <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
* WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
* UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
* THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
* OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
* OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
* USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
* THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
* DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
* UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
* BE CORRECTED.</strong></li>
* <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
* HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
* ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
* INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
* ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
* PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
* SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
* (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
* EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
* POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
* <ol></td></tr>
* </table>
* @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests)
* @author Burton S. Garbow (original fortran minpack tests)
* @author Kenneth E. Hillstrom (original fortran minpack tests)
* @author Jorge J. More (original fortran minpack tests)
* @author Luc Maisonobe (non-minpack tests and minpack tests Java translation)
*/
@Deprecated
public class LevenbergMarquardtOptimizerTest extends AbstractLeastSquaresOptimizerAbstractTest {
@Override
public AbstractLeastSquaresOptimizer createOptimizer() {
return new LevenbergMarquardtOptimizer();
}
@Override
@Test(expected=SingularMatrixException.class)
public void testNonInvertible() {
/*
* Overrides the method from parent class, since the default singularity
* threshold (1e-14) does not trigger the expected exception.
*/
LinearProblem problem = new LinearProblem(new double[][] {
{ 1, 2, -3 },
{ 2, 1, 3 },
{ -3, 0, -9 }
}, new double[] { 1, 1, 1 });
AbstractLeastSquaresOptimizer optimizer = createOptimizer();
PointVectorValuePair optimum = optimizer.optimize(100, problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0, 0 });
Assert.assertTrue(FastMath.sqrt(problem.target.length) * optimizer.getRMS() > 0.6);
optimizer.computeCovariances(optimum.getPoint(), 1.5e-14);
}
@Test
public void testControlParameters() {
CircleVectorial circle = new CircleVectorial();
circle.addPoint( 30.0, 68.0);
circle.addPoint( 50.0, -6.0);
circle.addPoint(110.0, -20.0);
circle.addPoint( 35.0, 15.0);
circle.addPoint( 45.0, 97.0);
checkEstimate(circle, 0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false);
checkEstimate(circle, 0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true);
checkEstimate(circle, 0.1, 5, 1.0e-15, 1.0e-16, 1.0e-10, true);
circle.addPoint(300, -300);
checkEstimate(circle, 0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true);
}
private void checkEstimate(MultivariateDifferentiableVectorFunction problem,
double initialStepBoundFactor, int maxCostEval,
double costRelativeTolerance, double parRelativeTolerance,
double orthoTolerance, boolean shouldFail) {
try {
LevenbergMarquardtOptimizer optimizer
= new LevenbergMarquardtOptimizer(initialStepBoundFactor,
costRelativeTolerance,
parRelativeTolerance,
orthoTolerance,
Precision.SAFE_MIN);
optimizer.optimize(maxCostEval, problem, new double[] { 0, 0, 0, 0, 0 },
new double[] { 1, 1, 1, 1, 1 },
new double[] { 98.680, 47.345 });
Assert.assertTrue(!shouldFail);
} catch (DimensionMismatchException ee) {
Assert.assertTrue(shouldFail);
} catch (TooManyEvaluationsException ee) {
Assert.assertTrue(shouldFail);
}
}
// Test is skipped because it fails with the latest code update.
@Ignore@Test
public void testMath199() {
try {
QuadraticProblem problem = new QuadraticProblem();
problem.addPoint (0, -3.182591015485607);
problem.addPoint (1, -2.5581184967730577);
problem.addPoint (2, -2.1488478161387325);
problem.addPoint (3, -1.9122489313410047);
problem.addPoint (4, 1.7785661310051026);
LevenbergMarquardtOptimizer optimizer
= new LevenbergMarquardtOptimizer(100, 1e-10, 1e-10, 1e-10, 0);
optimizer.optimize(100, problem,
new double[] { 0, 0, 0, 0, 0 },
new double[] { 0.0, 4.4e-323, 1.0, 4.4e-323, 0.0 },
new double[] { 0, 0, 0 });
Assert.fail("an exception should have been thrown");
} catch (ConvergenceException ee) {
// expected behavior
}
}
/**
* Non-linear test case: fitting of decay curve (from Chapter 8 of
* Bevington's textbook, "Data reduction and analysis for the physical sciences").
* XXX The expected ("reference") values may not be accurate and the tolerance too
* relaxed for this test to be currently really useful (the issue is under
* investigation).
*/
@Test
public void testBevington() {
final double[][] dataPoints = {
// column 1 = times
{ 15, 30, 45, 60, 75, 90, 105, 120, 135, 150,
165, 180, 195, 210, 225, 240, 255, 270, 285, 300,
315, 330, 345, 360, 375, 390, 405, 420, 435, 450,
465, 480, 495, 510, 525, 540, 555, 570, 585, 600,
615, 630, 645, 660, 675, 690, 705, 720, 735, 750,
765, 780, 795, 810, 825, 840, 855, 870, 885, },
// column 2 = measured counts
{ 775, 479, 380, 302, 185, 157, 137, 119, 110, 89,
74, 61, 66, 68, 48, 54, 51, 46, 55, 29,
28, 37, 49, 26, 35, 29, 31, 24, 25, 35,
24, 30, 26, 28, 21, 18, 20, 27, 17, 17,
14, 17, 24, 11, 22, 17, 12, 10, 13, 16,
9, 9, 14, 21, 17, 13, 12, 18, 10, },
};
final BevingtonProblem problem = new BevingtonProblem();
final int len = dataPoints[0].length;
final double[] weights = new double[len];
for (int i = 0; i < len; i++) {
problem.addPoint(dataPoints[0][i],
dataPoints[1][i]);
weights[i] = 1 / dataPoints[1][i];
}
final LevenbergMarquardtOptimizer optimizer
= new LevenbergMarquardtOptimizer();
final PointVectorValuePair optimum
= optimizer.optimize(100, problem, dataPoints[1], weights,
new double[] { 10, 900, 80, 27, 225 });
final double[] solution = optimum.getPoint();
final double[] expectedSolution = { 10.4, 958.3, 131.4, 33.9, 205.0 };
final double[][] covarMatrix = optimizer.computeCovariances(solution, 1e-14);
final double[][] expectedCovarMatrix = {
{ 3.38, -3.69, 27.98, -2.34, -49.24 },
{ -3.69, 2492.26, 81.89, -69.21, -8.9 },
{ 27.98, 81.89, 468.99, -44.22, -615.44 },
{ -2.34, -69.21, -44.22, 6.39, 53.80 },
{ -49.24, -8.9, -615.44, 53.8, 929.45 }
};
final int numParams = expectedSolution.length;
// Check that the computed solution is within the reference error range.
for (int i = 0; i < numParams; i++) {
final double error = FastMath.sqrt(expectedCovarMatrix[i][i]);
Assert.assertEquals("Parameter " + i, expectedSolution[i], solution[i], error);
}
// Check that each entry of the computed covariance matrix is within 10%
// of the reference matrix entry.
for (int i = 0; i < numParams; i++) {
for (int j = 0; j < numParams; j++) {
Assert.assertEquals("Covariance matrix [" + i + "][" + j + "]",
expectedCovarMatrix[i][j],
covarMatrix[i][j],
FastMath.abs(0.1 * expectedCovarMatrix[i][j]));
}
}
}
@Test
public void testCircleFitting2() {
final double xCenter = 123.456;
final double yCenter = 654.321;
final double xSigma = 10;
final double ySigma = 15;
final double radius = 111.111;
// The test is extremely sensitive to the seed.
final long seed = 59421061L;
final RandomCirclePointGenerator factory
= new RandomCirclePointGenerator(xCenter, yCenter, radius,
xSigma, ySigma,
seed);
final CircleProblem circle = new CircleProblem(xSigma, ySigma);
final int numPoints = 10;
for (Vector2D p : factory.generate(numPoints)) {
circle.addPoint(p);
// System.out.println(p.x + " " + p.y);
}
// First guess for the center's coordinates and radius.
final double[] init = { 90, 659, 115 };
final LevenbergMarquardtOptimizer optimizer
= new LevenbergMarquardtOptimizer();
final PointVectorValuePair optimum = optimizer.optimize(100, circle,
circle.target(), circle.weight(),
init);
final double[] paramFound = optimum.getPoint();
// Retrieve errors estimation.
final double[][] covMatrix = optimizer.computeCovariances(paramFound, 1e-14);
final double[] asymptoticStandardErrorFound = optimizer.guessParametersErrors();
final double[] sigmaFound = new double[covMatrix.length];
for (int i = 0; i < covMatrix.length; i++) {
sigmaFound[i] = FastMath.sqrt(covMatrix[i][i]);
// System.out.println("i=" + i + " value=" + paramFound[i]
// + " sigma=" + sigmaFound[i]
// + " ase=" + asymptoticStandardErrorFound[i]);
}
// System.out.println("chi2=" + optimizer.getChiSquare());
// Check that the parameters are found within the assumed error bars.
Assert.assertEquals(xCenter, paramFound[0], asymptoticStandardErrorFound[0]);
Assert.assertEquals(yCenter, paramFound[1], asymptoticStandardErrorFound[1]);
Assert.assertEquals(radius, paramFound[2], asymptoticStandardErrorFound[2]);
}
private static class QuadraticProblem implements MultivariateDifferentiableVectorFunction, Serializable {
private static final long serialVersionUID = 7072187082052755854L;
private List<Double> x;
private List<Double> y;
public QuadraticProblem() {
x = new ArrayList<Double>();
y = new ArrayList<Double>();
}
public void addPoint(double x, double y) {
this.x.add(x);
this.y.add(y);
}
public double[] value(double[] variables) {
double[] values = new double[x.size()];
for (int i = 0; i < values.length; ++i) {
values[i] = (variables[0] * x.get(i) + variables[1]) * x.get(i) + variables[2];
}
return values;
}
public DerivativeStructure[] value(DerivativeStructure[] variables) {
DerivativeStructure[] values = new DerivativeStructure[x.size()];
for (int i = 0; i < values.length; ++i) {
values[i] = (variables[0].multiply(x.get(i)).add(variables[1])).multiply(x.get(i)).add(variables[2]);
}
return values;
}
}
private static class BevingtonProblem
implements MultivariateDifferentiableVectorFunction {
private List<Double> time;
private List<Double> count;
public BevingtonProblem() {
time = new ArrayList<Double>();
count = new ArrayList<Double>();
}
public void addPoint(double t, double c) {
time.add(t);
count.add(c);
}
public double[] value(double[] params) {
double[] values = new double[time.size()];
for (int i = 0; i < values.length; ++i) {
final double t = time.get(i);
values[i] = params[0]
+ params[1] * FastMath.exp(-t / params[3])
+ params[2] * FastMath.exp(-t / params[4]);
}
return values;
}
public DerivativeStructure[] value(DerivativeStructure[] params) {
DerivativeStructure[] values = new DerivativeStructure[time.size()];
for (int i = 0; i < values.length; ++i) {
final double t = time.get(i);
values[i] = params[0].add(
params[1].multiply(params[3].reciprocal().multiply(-t).exp())).add(
params[2].multiply(params[4].reciprocal().multiply(-t).exp()));
}
return values;
}
}
}