package com.imsl.test.example.math;

import com.imsl.math.*;
import java.text.*;

/**
 * <p>
 * Approximates the gradient using central
 * divided differences.
 * </p>
 *
 * This example uses the same data as in example {@link NumericalDerivativesEx3}.
 * Instead of the one-sided difference, the central difference method is used.
 *
 * Agreement should be approximately the two-thirds power of machine precision.
 * That agreement is achieved here. Generally this is the <i>most</i> accuracy
 * one can expect using central divided differences. Note that using central
 * differences requires essentially twice the number of evaluations of the
 * function compared with obtaining one-sided differences. This can be a
 * significant issue for functions that are expensive to evaluate. This example
 * shows how to override <code>evaluateF</code>.
 *
 * @see <a href="NumericalDerivativesEx4.java">Code</a>
 * @see <a href="doc-files/NumericalDerivativesEx4.out">Output</a>
 */
public class NumericalDerivativesEx4 extends NumericalDerivatives {

    static private int m = 1, n = 2;
    static private double a, b, c, v = 0.0;

    public NumericalDerivativesEx4(NumericalDerivatives.Function fcn) {
        super(fcn);
    }

    // Override evaluateF.
    public double[] evaluateF(int varIndex, double[] y) {
        double[] valueF = new double[m];

        valueF[0] = a * Math.exp(b * y[0]) + c * y[0] * y[1] * y[1];
        return valueF;
    }

    public static void main(String args[]) {
        int[] options = new int[n];
        double u;
        double[] y = new double[n];
        double[] scale = new double[n];
        double[][] actual = new double[m][n];
        double[] re = new double[2];

        // Define data and point of evaluation:
        a = 2.5e6;
        b = 3.4e0;
        c = 4.5e0;
        y[0] = 2.1e0;
        y[1] = 3.2e0;

        // Machine precision, for measuring errors
        u = 2.220446049250313e-016;
        v = Math.pow(3.e0 * u, 2.e0 / 3.e0);

        // Set scaling:
        scale[0] = 1.e0;
        // Increase scale to account for large value of a.
        scale[1] = 8.e3;

        // Compute true values of partials.
        actual[0][0] = a * b * Math.exp(b * y[0]) + c * y[1] * y[1];
        actual[0][1] = 2 * c * y[0] * y[1];

        options[0] = NumericalDerivatives.CENTRAL;
        options[1] = NumericalDerivatives.CENTRAL;

        // Set the increment used at the default value.
        scale[1] = 8.e3;

        NumericalDerivatives.Function fcn
                = new NumericalDerivatives.Function() {
            public double[] f(int varIndex, double[] y) {
                return new double[m];
            }
        };

        NumericalDerivativesEx4 derv = new NumericalDerivativesEx4(fcn);
        derv.setDifferencingMethods(options);
        derv.setScalingFactors(scale);
        double[][] jacobian = derv.evaluateJ(y);

        NumberFormat nf = NumberFormat.getInstance();
        nf.setMaximumFractionDigits(2);
        nf.setMinimumFractionDigits(2);

        PrintMatrixFormat pmf = new PrintMatrixFormat();
        pmf.setNumberFormat(nf);
        new PrintMatrix("Numerical gradient:").print(pmf, jacobian);
        new PrintMatrix("Analytic gradient:").print(pmf, actual);

        // Since the function is never evaluated at the
        // initial point, hold back until the request is made.
        // Check the relative accuracy of central differences.
        // They should be good to about two thirds-precision.
        jacobian[0][0] = (jacobian[0][0] - actual[0][0]) / actual[0][0];
        jacobian[0][1] = (jacobian[0][1] - actual[0][1]) / actual[0][1];
        re[0] = jacobian[0][0];
        re[1] = jacobian[0][1];

        System.out.println("Relative accuracy:");
        System.out.println("df/dy_1       df/dy_2");
        System.out.printf(" %.2fv        %.2fv\n", re[0] / v, re[1] / v);
        System.out.printf("(%.3e)  (%.3e)\n", re[0], re[1]);
    }
}