package com.imsl.test.example.math;

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

/**
 * <p>
 * Approximates the Hessian of a function using
 * numerical differentiation. </p>
 *
 * This example uses the same data as in {@link NumericalDerivativesEx4}. In
 * this example numerical differentiation is used to approximate the Hessian
 * matrix of \(f(y_1 ,y_2 )\).
 *
 * This symmetric Hessian matrix is $$Hf = \left[ {\begin{array}{*{20}c}
 * {\frac{{\partial ^2 f}} {{\partial y_1^2 }}} & {\frac{{\partial ^2 f}}
 * {{\partial y_1 \partial y_2 }}} \\ {\frac{{\partial ^2 f}} {{\partial y_1
 * \partial y_2 }}} & {\frac{{\partial ^2 f}} {{\partial y_2^2 }}} \\
 * \end{array} } \right]$$
 *
 * Our method is based on casting the matrix \(Hf\) as the 2 by 2 Jacobian
 * matrix of the gradient function. Each inner evaluation of the gradient
 * function is itself computed using <code>NumericalDerivatives</code>. Central
 * differences are used for both the inner and outer numerical differentiation.
 * Because of the inherent error in both processes, the expected accuracy is
 * about the \(4/9 = \left( {2/3} \right)^2\) power of machine precision. Note
 * that the approximation obtained is not symmetric. However, the difference
 * between the off-diagonal elements provides an error estimate of that term.
 *
 * @see <a href="NumericalDerivativesEx5.java">Code</a>
 * @see <a href="doc-files/NumericalDerivativesEx5.out">Output</a>
 */
public class NumericalDerivativesEx5 extends NumericalDerivatives {

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

    class InnerNumericalDerivatives extends NumericalDerivatives {

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

        // Override evaluateF.
        @Override
        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 NumericalDerivativesEx5(NumericalDerivatives.Function fcn) {
        super(fcn);
    }

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

        // This is the analytic gradient.  for comparison only.
        //        stateh(1)=a*b*exp(b*y(1))+c*y(2)**2
        //        stateh(2)=2*c*y(1)*y(2)
        //
        // Each request for a gradient evaluation uses the
        // functionality of numerical evaluation, but with
        // the same numerical code.
        iopt[0] = NumericalDerivatives.CENTRAL;
        iopt[1] = NumericalDerivatives.CENTRAL;

        // Set the increment used at the default value.
        // Set defaults for increments and scaling:
        fac[0] = 1.4901161193847656E-8;
        fac[1] = 1.4901161193847656E-8;
        // Change scale to account for large value of a.
        switch (varIndex) {
            case 1:
                scale[0] = 1.e0;
                scale[1] = 8.e8;
                break;
            case 2:
                scale[0] = 1.e4;
                scale[1] = 8.e8;
                break;
        }

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

        InnerNumericalDerivatives derv = new InnerNumericalDerivatives(fcn);
        derv.setPercentageFactor(fac);
        derv.setDifferencingMethods(iopt);
        derv.setScalingFactors(scale);
        derv.setInitialF(valueF);
        double[][] fjac = derv.evaluateJ(y);

        // Since the function is never evaluated at the
        // initial point, hold back until the request is made.
        // Copy gradient value into array expected by
        // outer loop computing the Hessian matrix.
        double[] tmp = new double[n];
        for (int i = 0; i < n; i++) {
            tmp[i] = fjac[0][i];
        }
        return tmp;
    }

    public static void main(String args[]) {
        int[] iopth = new int[n];
        double u;
        double[] fach = new double[n];
        double[] stateh = new double[n];
        double[] scaleh = new double[n];
        double[][] actual = new double[n][n];
        double[] y = new double[n];

        // 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;

        // Compute expected relative error using two applications
        // of central differences.
        v = Math.pow(3.e0 * u, 2.e0 / 3.e0);
        v = Math.pow(3 * v, 2. / 3.);

        // Set increments and scaling:
        fach[0] = 1.4901161193847656E-8;
        fach[1] = 1.4901161193847656E-8;
        iopth[0] = NumericalDerivatives.CENTRAL;
        iopth[1] = NumericalDerivatives.CENTRAL;

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

        // Set the increment used at the default value.
        scaleh[0] = 1;
        scaleh[1] = 8.e5;

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

        NumericalDerivativesEx5 derv2 = new NumericalDerivativesEx5(fcn);
        derv2.setPercentageFactor(fach);
        derv2.setDifferencingMethods(iopth);
        derv2.setScalingFactors(scaleh);
        derv2.setInitialF(stateh);
        double[][] h = derv2.evaluateJ(y);

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

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

        // Since the function is never evaluated at the
        // initial point, hold back until the request is made.
        // Subtract the actual hessian matrix values and check.
        h[0][0] = (h[0][0] - actual[0][0]) / h[0][0] / v;
        h[1][0] = (h[1][0] - actual[1][0]) / h[1][0] / v;
        h[0][1] = (h[0][1] - actual[0][1]) / h[0][1] / v;
        h[1][1] = (h[1][1] - actual[1][1]) / h[1][1] / v;

        new PrintMatrix("Hessian Matrix, Expected Normalized "
                + "Relative Error, |all entries|").print(h);
    }
}