mixingcut 0.1.1

#![allow(non_snake_case)] // reasoning: The code is based on linear algebra notation (X is a matrix, x is a vector)

use clap::Parser;
use mixingcut::initialize::make_random_matrix;
use mixingcut::io_operations::write_solution_matrix;
use mixingcut::maxcut_oracle::{compute_rounded_sol, get_Q_norm, obj};
use mixingcut::step_rules::generate_step_rule;
use mixingcut::{io_operations, maxcut_oracle, sdp_local_search, step_rules};
use std::time::{SystemTime, UNIX_EPOCH};

#[derive(Parser, Debug)]
#[command(version, about, long_about = None)]
struct Args {
    // Name of the input file
    #[clap(short, long)]
    input_path: String,

    // Name of the output file
    #[clap(short, long, default_value = "output.txt")]
    output_path: String,

    #[clap(short, long, default_value = "0")]
    rank: usize,

    // The stopping tolerance
    #[clap(short, long, default_value = "1e-2")]
    tolerance: f64,

    // Number of iterations
    #[clap(short, long, default_value = "1000")]
    max_iters: usize,

    // Step Rule
    #[clap(short, long, default_value = "coord_no_step")]
    step_rule: String,

    // index correction
    #[clap(long, default_value = "1")]
    index_correction: usize,

    // compute dual bound
    #[clap(short, long, default_value = "0")]
    dual_bound: usize,

    // verbosity
    #[clap(short, long, default_value = "1")]
    verbose: usize,

    // rounding iters
    #[clap(long, default_value = "100")]
    rounding_iters: usize,

    // beam search width
    #[clap(long, default_value = "128")]
    beam_width: usize,
}

fn current_time() -> f64 {
    SystemTime::now()
        .duration_since(UNIX_EPOCH)
        .unwrap()
        .as_secs_f64()
}

fn main() {
    let args: Args = Args::parse();

    let index_correction = args.index_correction;

    // read in the graph
    let Q = io_operations::read_graph_matrix(&args.input_path, index_correction);

    let Q_norm = get_Q_norm(&Q);

    // compute the safe step size
    let alpha_safe = 1.0 / Q_norm;

    let n = Q.shape().0;

    let step_rule = generate_step_rule(&args.step_rule, alpha_safe);

    let max_iters = args.max_iters;

    let verbose = args.verbose;

    let max_rounding_iters = args.rounding_iters;

    // print the mixing cut vanity header if verbose
    if verbose == 1 {
        println!("------------------------------------------------------------------");
        println!("               MixingCut v0.0.1 - MAX CUT SDP Solver              ");
        println!("         (c) Dustin Kenefake, Texas A&M University, 2024          ");
        println!("------------------------------------------------------------------");
    }

    // set up the rank size of the problem
    let k = match args.rank {
        0 => 2 * (n as f64).log2() as usize,
        1 => (2.0 * n as f64).sqrt() as usize,
        _ => args.rank,
    };

    // generate random initial point
    let mut V = make_random_matrix(n, k, None);

    // print problem statistics if verbose
    if verbose == 1 {
        println!("Problem Statistics:");
        println!("Size of Q {} {}", n, n);
        println!("NNZ(Q) {}", Q.nnz());
        println!("Q norm {}", Q_norm);
        println!("Rank {}", k);
        println!("------------------------------------------------------------------");
        println!(
            "{0: <20} | {1: <20} | {2: <20}",
            "Iteration", "Primal Value", "Time(sec)"
        );
    }

    // get current time
    let start = current_time();

    // get the objective value
    let mut obj_val = obj(&Q, &V);

    // iterate over the number of iterations
    for i in 0..max_iters {
        // apply the step rule
        V = step_rules::apply_step(&Q, V, step_rule);

        // compute the objective value
        let new_obj_val = obj(&Q, &V);

        // if the objective value is not changing, break
        if (new_obj_val - obj_val).abs() < args.tolerance {
            if verbose == 1 {
                println!(
                    "{0: <20} | {1: <20} | {2: <20.6}",
                    i,
                    new_obj_val,
                    current_time() - start
                );
            }
            break;
        }

        // if the objective value is increasing, break
        if new_obj_val > obj_val {
            if verbose == 1 {
                println!("Objective value is increasing");
            }
            break;
        }

        obj_val = new_obj_val;

        // every 10 iterations, print the objective value if verbose
        if verbose == 1 && i % 10 == 0 {
            println!(
                "{0: <20} | {1: <20} | {2: <20.6}",
                i,
                obj_val,
                current_time() - start
            );
        }
    }

    if verbose == 1 {
        println!("------------------------------------------------------------------")
    }

    // compute the rounded solution
    let (x_0, obj_rounded) = compute_rounded_sol(&Q, &V, max_rounding_iters);

    if verbose == 1 {
        // print the rounded solution
        println!("Rounded solution: {:?} {:?}", obj_rounded, x_0.clone());

        // use beam search to generate better solutions
        let rounded_sols = vec![x_0.clone()];

        let (best_obj, best_sol) = sdp_local_search::beam_search(&Q, args.beam_width, rounded_sols);

        println!(
            "Rounded solution with local search: {:?} {:?}",
            best_obj, best_sol
        );
    }

    // print the dual bound
    if args.dual_bound == 1 {
        let dual_bound = maxcut_oracle::dual_bound(&Q, &V);

        if verbose == 1 {
            println!("Dual bound: {:?}", dual_bound);
        }
    }

    // write the solution to a file
    write_solution_matrix(&args.output_path, x_0, obj_rounded, obj(&Q, &V));
}