/*
    * 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.core.jdkmath.JdkMath;
  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 = JdkMath.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),
                          JdkMath.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()
         );
     }
 }