optimization-solvers 0.1.1

use nalgebra::DVector;
use optimization_solvers::{
    FuncEvalMultivariate, LineSearchSolver, MoreThuente, ProjectedGradientDescent, Tracer,
};

fn main() {
    // Setting up logging
    std::env::set_var("RUST_LOG", "info");
    let _ = Tracer::default().with_normal_stdout_layer().build();

    // Convex quadratic function: f(x,y) = (x-2)^2 + (y-3)^2
    // This function has a minimum at (2, 3), but we'll constrain it to a box
    let f_and_g = |x: &DVector<f64>| -> FuncEvalMultivariate {
        let x1 = x[0];
        let x2 = x[1];

        // Function value
        let f = (x1 - 2.0).powi(2) + (x2 - 3.0).powi(2);

        // Gradient
        let g1 = 2.0 * (x1 - 2.0);
        let g2 = 2.0 * (x2 - 3.0);
        let g = DVector::from_vec(vec![g1, g2]);

        FuncEvalMultivariate::new(f, g)
    };

    // Setting up the line search (backtracking)
    let mut ls = MoreThuente::default();

    // Setting up the solver with box constraints
    let tol = 1e-6;
    let x0 = DVector::from_vec(vec![0.0, 0.0]); // Starting point
    let lower_bound = DVector::from_vec(vec![0.0, 0.0]); // x >= 0, y >= 0
    let upper_bound = DVector::from_vec(vec![1.0, 1.0]); // x <= 1, y <= 1
    let mut solver =
        ProjectedGradientDescent::new(tol, x0.clone(), lower_bound.clone(), upper_bound.clone());

    // Running the solver
    let max_iter_solver = 100;
    let max_iter_line_search = 20;

    println!("=== Projected Gradient Descent Example ===");
    println!("Objective: f(x,y) = (x-2)^2 + (y-3)^2 (convex quadratic)");
    println!("Unconstrained minimum: (2, 3) with f(2,3) = 0");
    println!("Constraints: 0 <= x <= 1, 0 <= y <= 1");
    println!("Constrained minimum: (1, 1) with f(1,1) = 5");
    println!("Starting point: {:?}", x0);
    println!("Lower bounds: {:?}", lower_bound);
    println!("Upper bounds: {:?}", upper_bound);
    println!("Tolerance: {}", tol);
    println!();

    match solver.minimize(
        &mut ls,
        f_and_g,
        max_iter_solver,
        max_iter_line_search,
        None,
    ) {
        Ok(()) => {
            let x = solver.x();
            let eval = f_and_g(x);
            println!("✅ Optimization completed successfully!");
            println!("Final iterate: {:?}", x);
            println!("Function value: {:.6}", eval.f());
            println!("Gradient norm: {:.6}", eval.g().norm());
            println!("Iterations: {}", solver.k());

            // Check constraint satisfaction
            println!("Constraint satisfaction:");
            for i in 0..x.len() {
                println!(
                    "  x[{}] = {:.6} (bounds: [{:.1}, {:.1}])",
                    i, x[i], lower_bound[i], upper_bound[i]
                );
            }

            // The constrained minimum should be at (1, 1) since the unconstrained minimum (2, 3) is outside the box
            let expected_min = DVector::from_vec(vec![1.0, 1.0]);
            let distance_to_expected = (x - expected_min).norm();
            println!(
                "Distance to expected constrained minimum (1,1): {:.6}",
                distance_to_expected
            );
            println!("Expected function value at (1,1): 5.0");
        }
        Err(e) => {
            println!("❌ Optimization failed: {:?}", e);
        }
    }
}