/*

* 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.matrix.doubles.Matrix;

import dev.nm.algebra.linear.matrix.doubles.matrixtype.dense.DenseMatrix;

import dev.nm.algebra.linear.vector.doubles.Vector;

import dev.nm.algebra.linear.vector.doubles.dense.DenseVector;

import dev.nm.analysis.function.matrix.RntoMatrix;

import dev.nm.analysis.function.polynomial.Polynomial;

import dev.nm.analysis.function.rn2r1.AbstractBivariateRealFunction;

import dev.nm.analysis.function.rn2r1.RealScalarFunction;

import dev.nm.analysis.function.rn2rm.RealVectorFunction;

import dev.nm.analysis.function.special.gamma.LogGamma;

import dev.nm.solver.IterativeSolution;

import dev.nm.solver.multivariate.unconstrained.BruteForceMinimizer;

import dev.nm.solver.multivariate.unconstrained.DoubleBruteForceMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.conjugatedirection.ConjugateGradientMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.conjugatedirection.FletcherReevesMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.conjugatedirection.PowellMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.conjugatedirection.ZangwillMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.quasinewton.BFGSMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.quasinewton.DFPMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.quasinewton.HuangMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.quasinewton.PearsonMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.quasinewton.QuasiNewtonMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.quasinewton.RankOneMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.steepestdescent.FirstOrderMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.steepestdescent.GaussNewtonMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.steepestdescent.NewtonRaphsonMinimizer;

import dev.nm.solver.multivariate.unconstrained.c2.steepestdescent.SteepestDescentMinimizer;

import dev.nm.solver.problem.C2OptimProblem;

import dev.nm.solver.problem.C2OptimProblemImpl;

import dev.nm.solver.problem.OptimProblem;

import dev.nm.root.univariate.UnivariateMinimizer;

import dev.nm.root.univariate.bracketsearch.BracketSearchMinimizer;

import dev.nm.root.univariate.bracketsearch.BrentMinimizer;

import dev.nm.root.univariate.bracketsearch.FibonaccMinimizer;

import dev.nm.root.univariate.bracketsearch.GoldenMinimizer;

import static java.lang.Math.pow;

import static java.lang.Math.sqrt;

import java.util.ArrayList;

import java.util.List;

/**

* 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 Chapter9 {

public static void main(String[] args) throws Exception {

System.out.println("Chapter 9 demos");

Chapter9 chapter9 = new Chapter9();

chapter9.solve_by_brute_force_search_1();

chapter9.solve_by_brute_force_search_2();

chapter9.solve_by_brute_force_search_3();

chapter9.solve_by_brute_force_search_4();

chapter9.solve_loggamma_by_bracketing();

chapter9.solve_by_steepest_descent();

chapter9.solve_by_Newton_Raphson();

chapter9.solve_by_Gauss_Newton();

chapter9.solve_by_conjugate_direction_methods();

chapter9.solve_by_quasi_Newton();

}

public void solve_by_brute_force_search_1() throws Exception {

System.out.println("solve uniivariate function by brute force search");

// define the optimization problem using an objective function

OptimProblem problem = new OptimProblem() {

@Override

public int dimension() {

return 1;

}

@Override

public RealScalarFunction f() {

return new RealScalarFunction() {

// the objective function

@Override

public Double evaluate(Vector v) {

double x = v.get(1);

Polynomial polynomial = new Polynomial(1, 0, -4); // f(x) = x^2 - 4

double fx = polynomial.evaluate(x);

return fx;

}

@Override

public int dimensionOfDomain() {

return 1;

}

@Override

public int dimensionOfRange() {

return 1;

}

};

}

};

// set up the solver to use and the solution

DoubleBruteForceMinimizer solver = new DoubleBruteForceMinimizer(false);

BruteForceMinimizer.Solution soln = solver.solve(problem);

// for brute force search, we need to explicitly enumerate the values in the domain

List<Vector> domain = new ArrayList<>();

domain.add(new DenseVector(-2.));

domain.add(new DenseVector(-1.));

domain.add(new DenseVector(0.)); // the minimizer

domain.add(new DenseVector(1.));

domain.add(new DenseVector(2.));

soln.setDomain(domain);

System.out.println(String.format("f(%s) = %f", soln.minimizer(), soln.min()));

}

public void solve_by_brute_force_search_2() throws Exception {

System.out.println("solve multivariate function by brute force search");

OptimProblem problem = new OptimProblem() {

@Override

public int dimension() {

return 2;

}

@Override

public RealScalarFunction f() {

return new RealScalarFunction() {

@Override

public Double evaluate(Vector v) {

double x = v.get(1);

double y = v.get(2);

double fx = x * x + y * y;

return fx;

}

@Override

public int dimensionOfDomain() {

return 2;

}

@Override

public int dimensionOfRange() {

return 1;

}

};

}

};

DoubleBruteForceMinimizer bf = new DoubleBruteForceMinimizer(true);

BruteForceMinimizer.Solution soln = bf.solve(problem);

List<Vector> domain = new ArrayList<>();

domain.add(new DenseVector(-2., -2.));

domain.add(new DenseVector(-1., -1.));

domain.add(new DenseVector(0., 0.)); // the minimizer

domain.add(new DenseVector(1., 1.));

domain.add(new DenseVector(2., 2.));

soln.setDomain(domain);

System.out.println(String.format("f(%s) = %f", soln.minimizer(), soln.min()));

}

public void solve_by_brute_force_search_3() throws Exception {

System.out.println("solve uniivariate function by brute force search");

// set up the solver to use and the solution

DoubleBruteForceMinimizer solver = new DoubleBruteForceMinimizer(false);

BruteForceMinimizer.Solution soln

= solver.solve(new C2OptimProblemImpl(new Polynomial(1, 0, -4))); // f(x) = x^2 - 4

// for brute force search, we need to explicitly enumerate the values in the domain

List<Vector> domain = new ArrayList<>();

domain.add(new DenseVector(-2.));

domain.add(new DenseVector(-1.));

domain.add(new DenseVector(0.)); // the minimizer

domain.add(new DenseVector(1.));

domain.add(new DenseVector(2.));

soln.setDomain(domain);

System.out.println(String.format("f(%s) = %f", soln.minimizer(), soln.min()));

}

public void solve_by_brute_force_search_4() throws Exception {

System.out.println("solve multivariate function by brute force search");

DoubleBruteForceMinimizer bf = new DoubleBruteForceMinimizer(true);

BruteForceMinimizer.Solution soln = bf.solve(

new C2OptimProblemImpl(

new AbstractBivariateRealFunction() {

@Override

public double evaluate(double x, double y) {

double fx = x * x + y * y;

return fx;

}

}));

List<Vector> domain = new ArrayList<>();

domain.add(new DenseVector(-2., -2.));

domain.add(new DenseVector(-1., -1.));

domain.add(new DenseVector(0., 0.)); // the minimizer

domain.add(new DenseVector(1., 1.));

domain.add(new DenseVector(2., 2.));

soln.setDomain(domain);

System.out.println(String.format("f(%s) = %f", soln.minimizer(), soln.min()));

}

public void solve_loggamma_by_bracketing() throws Exception {

System.out.println("solve loggamma function by bracketing");

LogGamma logGamma = new LogGamma(); // the log-gamma function

BracketSearchMinimizer solver1 = new FibonaccMinimizer(1e-8, 15);

UnivariateMinimizer.Solution soln1 = solver1.solve(logGamma);

double x_min_1 = soln1.search(0, 5);

System.out.println(String.format("f(%f) = %f", x_min_1, logGamma.evaluate(x_min_1)));

BracketSearchMinimizer solver2 = new GoldenMinimizer(1e-8, 15);

UnivariateMinimizer.Solution soln2 = solver2.solve(logGamma);

double x_min_2 = soln2.search(0, 5);

System.out.println(String.format("f(%f) = %f", x_min_2, logGamma.evaluate(x_min_2)));

BracketSearchMinimizer solver3 = new BrentMinimizer(1e-8, 10);

UnivariateMinimizer.Solution soln3 = solver3.solve(logGamma);

double x_min_3 = soln3.search(0, 5);

System.out.println(String.format("f(%f) = %f", x_min_3, logGamma.evaluate(x_min_3)));

}

public void solve_by_steepest_descent() throws Exception {

System.out.println("solve multivariate function by steepest-descent");

// the objective function

// the global minimizer is at x = [0,0,0,0]

RealScalarFunction f = new RealScalarFunction() {

@Override

public Double evaluate(Vector x) {

double x1 = x.get(1);

double x2 = x.get(2);

double x3 = x.get(3);

double x4 = x.get(4);

double result = pow(x1 - 4 * x2, 4);

result += 12 * pow(x3 - x4, 4);

result += 3 * pow(x2 - 10 * x3, 2);

result += 55 * pow(x1 - 2 * x4, 2);

return result;

}

@Override

public int dimensionOfDomain() {

return 4;

}

@Override

public int dimensionOfRange() {

return 1;

}

};

// the gradient function

RealVectorFunction g = new RealVectorFunction() {

@Override

public Vector evaluate(Vector x) {

double x1 = x.get(1);

double x2 = x.get(2);

double x3 = x.get(3);

double x4 = x.get(4);

double[] result = new double[4];

result[0] = 4 * pow(x1 - 4 * x2, 3) + 110 * (x1 - 2 * x4);

result[1] = -16 * pow(x1 - 4 * x2, 3) + 6 * (x2 - 10 * x3);

result[2] = 48 * pow(x3 - x4, 3) - 60 * (x2 - 10 * x3);

result[3] = -48 * pow(x3 - x4, 3) - 220 * (x1 - 2 * x4);

return new DenseVector(result);

}

@Override

public int dimensionOfDomain() {

return 4;

}

@Override

public int dimensionOfRange() {

return 4;

}

};

C2OptimProblem problem = new C2OptimProblemImpl(f, g); // only gradient information

SteepestDescentMinimizer firstOrderMinimizer

= new FirstOrderMinimizer(

FirstOrderMinimizer.Method.IN_EXACT_LINE_SEARCH, // FirstOrderMinimizer.Method.ANALYTIC

1e-8,

40000

);

IterativeSolution<Vector> soln = firstOrderMinimizer.solve(problem);

Vector xmin = soln.search(new DenseVector(new double[]{1, -1, -1, 1}));

double f_xmin = f.evaluate(xmin);

System.out.println(String.format("f(%s) = %f", xmin.toString(), f_xmin));

}

public void solve_by_Newton_Raphson() throws Exception {

System.out.println("solve multivariate function by Newton-Raphson");

// the objective function

// the global minimizer is at x = [0,0,0,0]

RealScalarFunction f = new RealScalarFunction() {

@Override

public Double evaluate(Vector x) {

double x1 = x.get(1);

double x2 = x.get(2);

double x3 = x.get(3);

double x4 = x.get(4);

double result = pow(x1 - 4 * x2, 4);

result += 12 * pow(x3 - x4, 4);

result += 3 * pow(x2 - 10 * x3, 2);

result += 55 * pow(x1 - 2 * x4, 2);

return result;

}

@Override

public int dimensionOfDomain() {

return 4;

}

@Override

public int dimensionOfRange() {

return 1;

}

};

// the gradient function

RealVectorFunction g = new RealVectorFunction() {

@Override

public Vector evaluate(Vector x) {

double x1 = x.get(1);

double x2 = x.get(2);

double x3 = x.get(3);

double x4 = x.get(4);

double[] result = new double[4];

result[0] = 4 * pow(x1 - 4 * x2, 3) + 110 * (x1 - 2 * x4);

result[1] = -16 * pow(x1 - 4 * x2, 3) + 6 * (x2 - 10 * x3);

result[2] = 48 * pow(x3 - x4, 3) - 60 * (x2 - 10 * x3);

result[3] = -48 * pow(x3 - x4, 3) - 220 * (x1 - 2 * x4);

return new DenseVector(result);

}

@Override

public int dimensionOfDomain() {

return 4;

}

@Override

public int dimensionOfRange() {

return 4;

}

};

C2OptimProblem problem = new C2OptimProblemImpl(f, g); // use numerical Hessian

SteepestDescentMinimizer newtonRaphsonMinimizer

= new NewtonRaphsonMinimizer(

1e-8,

20

);

IterativeSolution<Vector> soln = newtonRaphsonMinimizer.solve(problem);

Vector xmin = soln.search(new DenseVector(new double[]{1, -1, -1, 1}));

double f_xmin = f.evaluate(xmin);

System.out.println(String.format("f(%s) = %f", xmin.toString(), f_xmin));

}

public void solve_by_Gauss_Newton() throws Exception {

System.out.println("solve multivariate function by Gauss-Newton");

// the objective function

// the global minimizer is at x = [0,0,0,0]

RealVectorFunction f = new RealVectorFunction() {

@Override

public Vector evaluate(Vector x) {

double x1 = x.get(1);

double x2 = x.get(2);

double x3 = x.get(3);

double x4 = x.get(4);

double[] fx = new double[4];

fx[0] = pow(x1 - 4 * x2, 2);

fx[1] = sqrt(12) * pow(x3 - x4, 2);

fx[2] = sqrt(3) * (x2 - 10 * x3);

fx[3] = sqrt(55) * (x1 - 2 * x4);

return new DenseVector(fx);

}

@Override

public int dimensionOfDomain() {

return 4;

}

@Override

public int dimensionOfRange() {

return 4;

}

};

// the Jacobian

RntoMatrix J = new RntoMatrix() {

@Override

public Matrix evaluate(Vector x) {

double x1 = x.get(1);

double x2 = x.get(2);

double x3 = x.get(3);

double x4 = x.get(4);

Matrix Jx = new DenseMatrix(4, 4);

double value = 2 * (x1 - 4 * x2);

Jx.set(1, 1, value);

value = -8 * (x1 - 4 * x2);

Jx.set(1, 2, value);

value = 2 * sqrt(12) * (x3 - x4);

Jx.set(2, 3, value);

Jx.set(2, 4, -value);

Jx.set(3, 2, sqrt(3));

Jx.set(3, 3, -10 * sqrt(3));

Jx.set(4, 1, sqrt(55));

Jx.set(4, 4, -2 * sqrt(55));

return Jx;

}

@Override

public int dimensionOfDomain() {

return 4;

}

@Override

public int dimensionOfRange() {

return 1;

}

};

GaussNewtonMinimizer optim1 = new GaussNewtonMinimizer(1e-8, 10);

IterativeSolution<Vector> soln = optim1.solve(f, J);//analytical gradient

Vector xmin = soln.search(new DenseVector(new double[]{1, -1, -1, 1}));

System.out.println(String.format("f(%s) = %s", xmin.toString(), f.evaluate(xmin).toString()));

}

public void solve_by_conjugate_direction_methods() throws Exception {

System.out.println("solve multivariate function by conjugate-direction methods");

/**

* The Himmelblau function: f(x) = (x1^2 + x2 - 11)^2 + (x1 + x2^2 -

* 7)^2

*/

RealScalarFunction f = new RealScalarFunction() {

@Override

public Double evaluate(Vector x) {

double x1 = x.get(1);

double x2 = x.get(2);

double result = pow(x1 * x1 + x2 - 11, 2);

result += pow(x1 + x2 * x2 - 7, 2);

return result;

}

@Override

public int dimensionOfDomain() {

return 2;

}

@Override

public int dimensionOfRange() {

return 1;

}

};

RealVectorFunction g = new RealVectorFunction() {

@Override

public Vector evaluate(Vector x) {

double x1 = x.get(1);

double x2 = x.get(2);

double w1 = x1 * x1 + x2 - 11;

double w2 = x1 + x2 * x2 - 7;

double[] result = new double[2];

result[0] = 4 * w1 * x1 + 2 * w2;

result[1] = 2 * w1 + 4 * w2 * x2;

return new DenseVector(result);

}

@Override

public int dimensionOfDomain() {

return 2;

}

@Override

public int dimensionOfRange() {

return 2;

}

};

C2OptimProblemImpl problem = new C2OptimProblemImpl(f, g);

ConjugateGradientMinimizer ConjugateGradientMinimizer

= new ConjugateGradientMinimizer(1e-16, 40);

IterativeSolution<Vector> soln1 = ConjugateGradientMinimizer.solve(problem);

Vector xmin1 = soln1.search(new DenseVector(new double[]{6, 6}));

double f_xmin1 = f.evaluate(xmin1);

System.out.println(String.format("f(%s) = %.16f", xmin1.toString(), f_xmin1));

ConjugateGradientMinimizer fletcherReevesMinimizer

= new FletcherReevesMinimizer(1e-16, 20);

IterativeSolution<Vector> soln2 = fletcherReevesMinimizer.solve(problem);

Vector xmin2 = soln2.search(new DenseVector(new double[]{6, 6}));

double f_xmin2 = f.evaluate(xmin2);

System.out.println(String.format("f(%s) = %.16f", xmin2.toString(), f_xmin2));

SteepestDescentMinimizer powellMinimizer

= new PowellMinimizer(1e-16, 20);

IterativeSolution<Vector> soln3 = powellMinimizer.solve(problem);

Vector xmin3 = soln3.search(new DenseVector(new double[]{6, 6}));

double f_xmin3 = f.evaluate(xmin3);

System.out.println(String.format("f(%s) = %.16f", xmin3.toString(), f_xmin3));

SteepestDescentMinimizer zangwillMinimizer

= new ZangwillMinimizer(1e-16, 1e-16, 20);

IterativeSolution<Vector> soln4 = zangwillMinimizer.solve(problem);

Vector xmin4 = soln4.search(new DenseVector(new double[]{6, 6}));

double f_xmin4 = f.evaluate(xmin4);

System.out.println(String.format("f(%s) = %.16f", xmin4.toString(), f_xmin4));

}

public void solve_by_quasi_Newton() throws Exception {

System.out.println("solve multivariate function by quasi-Newton");

/**

* The Himmelblau function: f(x) = (x1^2 + x2 - 11)^2 + (x1 + x2^2 -

* 7)^2

*/

RealScalarFunction f = new RealScalarFunction() {

@Override

public Double evaluate(Vector x) {

double x1 = x.get(1);

double x2 = x.get(2);

double result = pow(x1 * x1 + x2 - 11, 2);

result += pow(x1 + x2 * x2 - 7, 2);

return result;

}

@Override

public int dimensionOfDomain() {

return 2;

}

@Override

public int dimensionOfRange() {

return 1;

}

};

RealVectorFunction g = new RealVectorFunction() {

@Override

public Vector evaluate(Vector x) {

double x1 = x.get(1);

double x2 = x.get(2);

double w1 = x1 * x1 + x2 - 11;

double w2 = x1 + x2 * x2 - 7;

double[] result = new double[2];

result[0] = 4 * w1 * x1 + 2 * w2;

result[1] = 2 * w1 + 4 * w2 * x2;

return new DenseVector(result);

}

@Override

public int dimensionOfDomain() {

return 2;

}

@Override

public int dimensionOfRange() {

return 2;

}

};

C2OptimProblemImpl problem = new C2OptimProblemImpl(f, g);

QuasiNewtonMinimizer rankOneMinimizer = new RankOneMinimizer(1e-16, 15);

IterativeSolution<Vector> soln1 = rankOneMinimizer.solve(problem);

Vector xmin = soln1.search(new DenseVector(new double[]{6, 6}));

double f_xmin = f.evaluate(xmin);

System.out.println(String.format("f(%s) = %.16f", xmin.toString(), f_xmin));

QuasiNewtonMinimizer dfpMinimizer = new DFPMinimizer(1e-16, 15);

IterativeSolution<Vector> soln2 = dfpMinimizer.solve(problem);

xmin = soln2.search(new DenseVector(new double[]{6, 6}));

f_xmin = f.evaluate(xmin);

System.out.println(String.format("f(%s) = %.16f", xmin.toString(), f_xmin));

QuasiNewtonMinimizer bfgsMinimizer = new BFGSMinimizer(false, 1e-16, 15);

IterativeSolution<Vector> soln3 = bfgsMinimizer.solve(problem);

xmin = soln3.search(new DenseVector(new double[]{6, 6}));

f_xmin = f.evaluate(xmin);

System.out.println(String.format("f(%s) = %.16f", xmin.toString(), f_xmin));

QuasiNewtonMinimizer huangMinimizer = new HuangMinimizer(0, 1, 0, 1, 1e-16, 15);

IterativeSolution<Vector> soln4 = huangMinimizer.solve(problem);

xmin = soln4.search(new DenseVector(new double[]{6, 6}));

f_xmin = f.evaluate(xmin);

System.out.println(String.format("f(%s) = %.16f", xmin.toString(), f_xmin));

QuasiNewtonMinimizer pearsonMinimizer = new PearsonMinimizer(1e-16, 15);

IterativeSolution<Vector> soln5 = pearsonMinimizer.solve(problem);

xmin = soln5.search(new DenseVector(new double[]{6, 6}));

f_xmin = f.evaluate(xmin);

System.out.println(String.format("f(%s) = %.16f", xmin.toString(), f_xmin));

}

}