/*
 * Copyright (c) NM LTD.
 * https://nm.dev/
 * 
 * THIS SOFTWARE IS LICENSED, NOT SOLD.
 * 
 * YOU MAY USE THIS SOFTWARE ONLY AS DESCRIBED IN THE LICENSE.
 * IF YOU ARE NOT AWARE OF AND/OR DO NOT AGREE TO THE TERMS OF THE LICENSE,
 * DO NOT USE THIS SOFTWARE.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITH NO WARRANTY WHATSOEVER,
 * EITHER EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION,
 * ANY WARRANTIES OF ACCURACY, ACCESSIBILITY, COMPLETENESS,
 * FITNESS FOR A PARTICULAR PURPOSE, MERCHANTABILITY, NON-INFRINGEMENT, 
 * TITLE AND USEFULNESS.
 * 
 * IN NO EVENT AND UNDER NO LEGAL THEORY,
 * WHETHER IN ACTION, CONTRACT, NEGLIGENCE, TORT, OR OTHERWISE,
 * SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
 * ANY CLAIMS, DAMAGES OR OTHER LIABILITIES,
 * ARISING AS A RESULT OF USING OR OTHER DEALINGS IN THE SOFTWARE.
 */
package dev.nm.nmj;

import dev.nm.algebra.linear.vector.doubles.Vector;
import dev.nm.algebra.linear.vector.doubles.dense.DenseVector;
import dev.nm.analysis.differentialequation.ode.ivp.problem.DerivativeFunction;
import dev.nm.analysis.differentialequation.ode.ivp.problem.ODE1stOrder;
import dev.nm.analysis.differentialequation.ode.ivp.solver.EulerMethod;
import dev.nm.analysis.differentialequation.ode.ivp.solver.ODESolution;
import dev.nm.analysis.differentialequation.ode.ivp.solver.ODESolver;
import dev.nm.analysis.differentialequation.ode.ivp.solver.multistep.adamsbashforthmoulton.ABMPredictorCorrector1;
import dev.nm.analysis.differentialequation.ode.ivp.solver.multistep.adamsbashforthmoulton.ABMPredictorCorrector2;
import dev.nm.analysis.differentialequation.ode.ivp.solver.multistep.adamsbashforthmoulton.ABMPredictorCorrector3;
import dev.nm.analysis.differentialequation.ode.ivp.solver.multistep.adamsbashforthmoulton.ABMPredictorCorrector4;
import dev.nm.analysis.differentialequation.ode.ivp.solver.multistep.adamsbashforthmoulton.ABMPredictorCorrector5;
import dev.nm.analysis.differentialequation.ode.ivp.solver.multistep.adamsbashforthmoulton.AdamsBashforthMoulton;
import dev.nm.analysis.differentialequation.ode.ivp.solver.rungekutta.RungeKutta;
import dev.nm.analysis.differentialequation.ode.ivp.solver.rungekutta.RungeKutta1;
import dev.nm.analysis.differentialequation.ode.ivp.solver.rungekutta.RungeKutta2;
import dev.nm.analysis.differentialequation.ode.ivp.solver.rungekutta.RungeKutta3;
import dev.nm.analysis.differentialequation.ode.ivp.solver.rungekutta.RungeKuttaStepper;
import dev.nm.analysis.function.rn2r1.univariate.AbstractUnivariateRealFunction;
import dev.nm.analysis.function.rn2r1.univariate.UnivariateRealFunction;
import dev.nm.analysis.function.rn2rm.AbstractRealVectorFunction;
import dev.nm.analysis.function.rn2rm.RealVectorFunction;
import java.io.IOException;
import static java.lang.Math.exp;

/**
 * Numerical Methods Using Java: For Data Science, Analysis, and Engineering
 *
 * @author haksunli
 * @see
 * https://www.amazon.com/Numerical-Methods-Using-Java-Engineering/dp/1484267966
 * https://nm.dev/
 */
public class Chapter7 {

    public static void main(String[] args) throws IOException {
        System.out.println("Chapter 7 demos");

        Chapter7 chapter7 = new Chapter7();
        chapter7.EulerMethod();
        chapter7.RungeKutta();
        chapter7.ABM();
        chapter7.system_of_ODEs();
        chapter7.higher_order_ODEs();
    }

    private void higher_order_ODEs() {
        System.out.println("solve a higher-order ODE using Euler's method");

        // define the equivalent system of ODEs to solve
        DerivativeFunction dY = new DerivativeFunction() {

            @Override
            public Vector evaluate(double t, Vector v) {
                double y_1 = v.get(1);
                double y_2 = v.get(2);

                double dy_1 = y_2;
                double dy_2 = y_2 + 6. * y_1;
                return new DenseVector(new double[]{dy_1, dy_2});
            }

            @Override
            public int dimension() {
                return 2;
            }
        };
        // initial condition, y1(0) = 1, y2(0) = 2
        Vector Y0 = new DenseVector(1., 2.);

        double t0 = 0, t1 = 1.; // solution domain
        double h = 0.1; // step size

        // the analytical solution
        UnivariateRealFunction f
                = new AbstractUnivariateRealFunction() {
            @Override
            public double evaluate(double t) {
                double f = 0.8 * exp(3. * t);
                f += 0.2 * exp(-2. * t);
                return f;
            }
        };

        // define an IVP
        ODE1stOrder ivp = new ODE1stOrder(dY, Y0, t0, t1);
        // construt an ODE solver using Euler's method
//        ODESolver solver = new EulerMethod(h);
        // construt an ODE solver using the third order Runge-Kutta formula
        RungeKuttaStepper stepper3 = new RungeKutta3();
        ODESolver solver = new RungeKutta(stepper3, h);

        // solve the ODE
        ODESolution soln = solver.solve(ivp);
        // print out the solution function, y, at discrete points
        double[] t = soln.x();
        Vector[] v = soln.y();
        for (int i = 0; i < t.length; ++i) {
            double y1 = v[i].get(1); // the numerical solution
            System.out.println(String.format(
                    "y(%f) = %f vs %f",
                    t[i],
                    y1,
                    f.evaluate(t[i])
            ));
        }
    }

    private void system_of_ODEs() {
        System.out.println("solve a system of ODEs using Euler's method");

        // define the system of ODEs to solve
        DerivativeFunction dY = new DerivativeFunction() {

            @Override
            public Vector evaluate(double t, Vector v) {
                double x = v.get(1);
                double y = v.get(2);

                double dx = 3. * x - 4. * y;
                double dy = 4. * x - 7. * y;
                return new DenseVector(new double[]{dx, dy});
            }

            @Override
            public int dimension() {
                return 2;
            }
        };
        // initial condition, x0=y0=1
        Vector Y0 = new DenseVector(1., 1.);

        double x0 = 0, x1 = 1.; // solution domain
        double h = 0.1; // step size

        // the analytical solution
        RealVectorFunction F
                = new AbstractRealVectorFunction(1, 2) {
            @Override
            public Vector evaluate(Vector v) {
                double t = v.get(1);
                double x = 2. / 3 * exp(t) + 1. / 3 * exp(-5. * t);
                double y = 1. / 3 * exp(t) + 2. / 3 * exp(-5. * t);
                return new DenseVector(x, y);
            }
        };

        // define an IVP
        ODE1stOrder ivp = new ODE1stOrder(dY, Y0, x0, x1);
        // construt an ODE solver using Euler's method
        ODESolver solver = new EulerMethod(h);
        // solve the ODE
        ODESolution soln = solver.solve(ivp);
        // print out the solution function, y, at discrete points
        double[] t = soln.x();
        Vector[] v = soln.y();
        for (int i = 0; i < t.length; ++i) {
            System.out.println(String.format(
                    "y(%f) = %s vs %s",
                    t[i],
                    v[i],
                    F.evaluate(new DenseVector(t[i]))
            ));
        }
    }

    private void ABM() {
        System.out.println("solve an ODE using Adams-Bashforth methods");

        // define the ODE to solve
        DerivativeFunction dy = new DerivativeFunction() {

            @Override
            public Vector evaluate(double x, Vector v) {
                double y = v.get(1);
                double dy = y - x + 1;
                return new DenseVector(dy);
            }

            @Override
            public int dimension() {
                return 1;
            }
        };
        // initial condition, y0=1
        Vector y0 = new DenseVector(1.);

        double x0 = 0, x1 = 1.; // solution domain
        double h = 0.1; // step size

        // the analytical solution
        UnivariateRealFunction y = new AbstractUnivariateRealFunction() {
            @Override
            public double evaluate(double x) {
                double y = exp(x) + x;
                return y;
            }
        };

        // define an IVP
        ODE1stOrder ivp = new ODE1stOrder(dy, y0, x0, x1);

        // using first order Adams-Bashforth formula
        ODESolver solver1 = new AdamsBashforthMoulton(new ABMPredictorCorrector1(), h);
        ODESolution soln1 = solver1.solve(ivp);

        // using second order Adams-Bashforth formula
        ODESolver solver2 = new AdamsBashforthMoulton(new ABMPredictorCorrector2(), h);
        ODESolution soln2 = solver2.solve(ivp);

        // using third order Adams-Bashforth formula
        ODESolver solver3 = new AdamsBashforthMoulton(new ABMPredictorCorrector3(), h);
        ODESolution soln3 = solver3.solve(ivp);

        // using forth order Adams-Bashforth formula
        ODESolver solver4 = new AdamsBashforthMoulton(new ABMPredictorCorrector4(), h);
        ODESolution soln4 = solver4.solve(ivp);

        // using fifth order Adams-Bashforth formula
        ODESolver solver5 = new AdamsBashforthMoulton(new ABMPredictorCorrector5(), h);
        ODESolution soln5 = solver5.solve(ivp);

        double[] x = soln1.x();
        Vector[] y1 = soln1.y();
        Vector[] y2 = soln2.y();
        Vector[] y3 = soln3.y();
        Vector[] y4 = soln4.y();
        Vector[] y5 = soln5.y();
        for (int i = 0; i < x.length; ++i) {
            double yx = y.evaluate(x[i]); // the analytical solution
            double diff1 = yx - y1[i].get(1); // the first order error
            double diff2 = yx - y2[i].get(1); // the second order error
            double diff3 = yx - y3[i].get(1); // the third order error
            double diff4 = yx - y4[i].get(1); // the forth order error
            double diff5 = yx - y5[i].get(1); // the fifth order error
            System.out.println(
                    String.format("y(%f) = %s (%.16f); = %s (%.16f); = %s (%.16f); = %s (%.16f); = %s (%.16f)",
                            x[i], y1[i], diff1,
                            y2[i], diff2,
                            y3[i], diff3,
                            y4[i], diff4,
                            y5[i], diff5
                    ));
        }
    }

    private void RungeKutta() {
        System.out.println("solve an ODE using Runge-Kutta methods");

        // define the ODE to solve
        DerivativeFunction dy = new DerivativeFunction() {

            @Override
            public Vector evaluate(double x, Vector v) {
                double y = v.get(1);
                double dy = y - x + 1;
                return new DenseVector(dy);
            }

            @Override
            public int dimension() {
                return 1;
            }
        };
        // initial condition, y0=1
        Vector y0 = new DenseVector(1.);

        double x0 = 0, x1 = 1.; // solution domain
        double h = 0.1; // step size

        // the analytical solution
        UnivariateRealFunction y = new AbstractUnivariateRealFunction() {
            @Override
            public double evaluate(double x) {
                double y = exp(x) + x;
                return y;
            }
        };

        // define an IVP
        ODE1stOrder ivp = new ODE1stOrder(dy, y0, x0, x1);

        // using first order Runge-Kutta formula
        RungeKuttaStepper stepper1 = new RungeKutta1();
        ODESolver solver1 = new RungeKutta(stepper1, h);
        ODESolution soln1 = solver1.solve(ivp);

        // using second order Runge-Kutta formula
        RungeKuttaStepper stepper2 = new RungeKutta2();
        ODESolver solver2 = new RungeKutta(stepper2, h);
        ODESolution soln2 = solver2.solve(ivp);

        // using third order Runge-Kutta formula
        RungeKuttaStepper stepper3 = new RungeKutta3();
        ODESolver solver3 = new RungeKutta(stepper3, h);
        ODESolution soln3 = solver3.solve(ivp);

        double[] x = soln1.x();
        Vector[] y1 = soln1.y();
        Vector[] y2 = soln2.y();
        Vector[] y3 = soln3.y();
        for (int i = 0; i < x.length; ++i) {
            double yx = y.evaluate(x[i]); // the analytical solution
            double diff1 = yx - y1[i].get(1); // the first order error
            double diff2 = yx - y2[i].get(1); // the second order error
            double diff3 = yx - y3[i].get(1); // the third order error
            System.out.println(
                    String.format("y(%f) = %s (%.16f); = %s (%.16f); = %s (%.16f)",
                            x[i], y1[i], diff1,
                            y2[i], diff2,
                            y3[i], diff3
                    ));
        }
    }

    private void EulerMethod() {
        System.out.println("solve an ODE using Euler's method");

        // define the ODE to solve
        DerivativeFunction dy = new DerivativeFunction() {

            @Override
            public Vector evaluate(double x, Vector y) {

                Vector dy = y.scaled(-2. * x);
                return dy.add(1); // y' = 1 - 2xy
            }

            @Override
            public int dimension() {
                return 1;
            }
        };
        // initial condition, y0=0
        Vector y0 = new DenseVector(0.);

        double x0 = 0, x1 = 1.; // solution domain
        double h = 0.1; // step size

        // define an IVP
        ODE1stOrder ivp = new ODE1stOrder(dy, y0, x0, x1);
        // construt an ODE solver using Euler's method
        ODESolver solver = new EulerMethod(h);
        // solve the ODE
        ODESolution soln = solver.solve(ivp);
        // print out the solution function, y, at discrete points
        double[] x = soln.x();
        Vector[] y = soln.y();
        for (int i = 0; i < x.length; ++i) {
            System.out.println(String.format("y(%f) = %s", x[i], y[i]));
        }
    }

}