commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/fitting/leastsquares/DifferentiatorVectorMultivariateJacobianFunctionTest.java

/*
 * 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.math4.legacy.fitting.leastsquares;

import org.apache.commons.math4.legacy.analysis.differentiation.FiniteDifferencesDifferentiator;
import org.apache.commons.math4.legacy.analysis.differentiation.UnivariateVectorFunctionDifferentiator;
import org.apache.commons.math4.legacy.linear.DiagonalMatrix;
import org.apache.commons.math4.legacy.linear.RealMatrix;
import org.apache.commons.math4.legacy.linear.RealVector;
import org.apache.commons.math4.legacy.optim.SimpleVectorValueChecker;
import org.apache.commons.math4.legacy.core.jdkmath.AccurateMath;
import org.junit.Assert;
import org.junit.Test;

public class DifferentiatorVectorMultivariateJacobianFunctionTest {

    private static final int POINTS = 20;
    private static final double STEP_SIZE = 0.2;

    private final UnivariateVectorFunctionDifferentiator differentiator = new FiniteDifferencesDifferentiator(POINTS, STEP_SIZE);
    private final LeastSquaresOptimizer optimizer = this.getOptimizer();

    public LeastSquaresBuilder base() {
        return new LeastSquaresBuilder()
                .checkerPair(new SimpleVectorValueChecker(1e-6, 1e-6))
                .maxEvaluations(100)
                .maxIterations(getMaxIterations());
    }

    public LeastSquaresBuilder builder(BevingtonProblem problem, boolean automatic){
        if(automatic) {
            return base()
                    .model(new DifferentiatorVectorMultivariateJacobianFunction(problem.getModelFunction(), differentiator));
        } else {
            return base()
                    .model(problem.getModelFunction(), problem.getModelFunctionJacobian());
        }
    }

    public int getMaxIterations() {
        return 25;
    }

    public LeastSquaresOptimizer getOptimizer() {
        return new LevenbergMarquardtOptimizer();
    }

    /**
     * Non-linear test case: fitting of decay curve (from Chapter 8 of
     * Bevington's textbook, "Data reduction and analysis for the physical sciences").
     */
    @Test
    public void testBevington() {

        // the analytical optimum to compare to
        final LeastSquaresOptimizer.Optimum analyticalOptimum = findBevington(false);
        final RealVector analyticalSolution = analyticalOptimum.getPoint();
        final RealMatrix analyticalCovarianceMatrix = analyticalOptimum.getCovariances(1e-14);

        // the automatic DifferentiatorVectorMultivariateJacobianFunction optimum
        final LeastSquaresOptimizer.Optimum automaticOptimum = findBevington(true);
        final RealVector automaticSolution = automaticOptimum.getPoint();
        final RealMatrix automaticCovarianceMatrix = automaticOptimum.getCovariances(1e-14);

        final int numParams = analyticalOptimum.getPoint().getDimension();

        // Check that the automatic solution is within the reference error range.
        for (int i = 0; i < numParams; i++) {
            final double error = AccurateMath.sqrt(analyticalCovarianceMatrix.getEntry(i, i));
            Assert.assertEquals("Parameter " + i, analyticalSolution.getEntry(i), automaticSolution.getEntry(i), error);
        }

        // Check that each entry of the computed covariance matrix is within 1%
        // of the reference analytical matrix entry.
        for (int i = 0; i < numParams; i++) {
            for (int j = 0; j < numParams; j++) {
                Assert.assertEquals("Covariance matrix [" + i + "][" + j + "]",
                        analyticalCovarianceMatrix.getEntry(i, j),
                        automaticCovarianceMatrix.getEntry(i, j),
                        AccurateMath.abs(0.01 * analyticalCovarianceMatrix.getEntry(i, j)));
            }
        }

        // Check various measures of goodness-of-fit.
        final double tol = 1e-40;
        Assert.assertEquals(analyticalOptimum.getChiSquare(), automaticOptimum.getChiSquare(), tol);
        Assert.assertEquals(analyticalOptimum.getCost(), automaticOptimum.getCost(), tol);
        Assert.assertEquals(analyticalOptimum.getRMS(), automaticOptimum.getRMS(), tol);
        Assert.assertEquals(analyticalOptimum.getReducedChiSquare(automaticOptimum.getPoint().getDimension()), automaticOptimum.getReducedChiSquare(automaticOptimum.getPoint().getDimension()), tol);
    }

    /**
     * Build the problem and return the optimum, doesn't actually test the results.
     *
     * Pass in if you want to test using analytical derivatives,
     * or the automatic {@link DifferentiatorVectorMultivariateJacobianFunction}
     *
     * @param automatic automatic {@link DifferentiatorVectorMultivariateJacobianFunction}, as opposed to analytical
     * @return the optimum for this test
     */
    private LeastSquaresOptimizer.Optimum findBevington(boolean automatic) {
        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 double[] start = {10, 900, 80, 27, 225};

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

        return optimizer.optimize(
                builder(problem, automatic)
                        .target(dataPoints[1])
                        .weight(new DiagonalMatrix(weights))
                        .start(start)
                        .build()
        );
    }
}