//! Solve a G-set Max-Cut instance with simulated bifurcation and report the
//! best cut, how often the target is reached, and the steps to reach it with
//! 99% probability.
//!
//! ```text
//! cargo run --release --example qio_gset --features quantum_qio -- FILE VARIANT STEPS TRIALS [TARGET] [DT] [C1] [GAMMA] [THREADS]
//! ```
//! VARIANT is bsb, dsb, hbsb or hdsb. Omitted settings take the published values.
use wai_quantum::quantum_qio::{descend, read_gset, simulated_annealing, simulated_bifurcation, step_to_solution, SaConfig, SbConfig, Variant};

fn main() {
    let a: Vec<String> = std::env::args().collect();
    let g = read_gset(&std::fs::read_to_string(&a[1]).expect("read")).expect("parse");
    let polish = a[2].ends_with("+ls");
    let name = a[2].trim_end_matches("+ls");
    if name == "sa" {
        // STEPS sweeps; DT and C1 slots carry beta0 and beta1.
        let sweeps: u32 = a[3].parse().unwrap();
        let trials: usize = a[4].parse().unwrap();
        let target: Option<f64> = a.get(5).and_then(|x| x.parse().ok());
        let beta0 = a.get(6).and_then(|x| x.parse().ok()).unwrap_or(0.1);
        let beta1 = a.get(7).and_then(|x| x.parse().ok()).unwrap_or(3.0);
        let threads: usize = a.get(9).map_or(8, |x| x.parse().unwrap());
        let t0 = std::time::Instant::now();
        let run = simulated_annealing(&g.ising(), &SaConfig { sweeps, beta0, beta1, trials, seed: 1 }, threads);
        let cuts: Vec<f64> = run.trials.iter().map(|t| g.cut(&t.spins)).collect();
        let best = cuts.iter().cloned().fold(f64::MIN, f64::max);
        print!("SA beta {beta0}->{beta1} sweeps={sweeps} trials={trials}  best {best} mean {:.1}", cuts.iter().sum::<f64>() / cuts.len() as f64);
        if let Some(t) = target {
            let p = cuts.iter().filter(|&&c| c >= t).count() as f64 / cuts.len() as f64;
            print!("  P(≥{t}) {p:.4}  S99 {:?}", step_to_solution(p, f64::from(sweeps)).map(|s| s.round()));
        }
        println!("  {:.1}s", t0.elapsed().as_secs_f64());
        return;
    }
    let v = match name {
        "bsb" => Variant::Ballistic,
        "dsb" => Variant::Discrete,
        "hbsb" => Variant::HeatedBallistic(0.0),
        "hdsb" => Variant::HeatedDiscrete(0.0),
        other => panic!("variant {other}"),
    };
    let steps: u32 = a[3].parse().unwrap();
    let trials: usize = a[4].parse().unwrap();
    let mut cfg = SbConfig::published(v, steps, trials, 1);
    let target: Option<f64> = a.get(5).and_then(|x| x.parse().ok());
    if let Some(dt) = a.get(6).and_then(|x| x.parse().ok()) {
        cfg.dt = dt;
    }
    if let Some(c1) = a.get(7).and_then(|x| x.parse().ok()) {
        cfg.c1 = c1;
    }
    if let Some(gamma) = a.get(8).and_then(|x| x.parse::<f64>().ok()) {
        cfg.variant = match cfg.variant {
            Variant::HeatedBallistic(_) => Variant::HeatedBallistic(gamma),
            Variant::HeatedDiscrete(_) => Variant::HeatedDiscrete(gamma),
            other => other,
        };
    }
    let threads: usize = a.get(9).map_or(8, |x| x.parse().unwrap());
    let t0 = std::time::Instant::now();
    let run = simulated_bifurcation(&g.ising(), &cfg, threads);
    let dt = t0.elapsed().as_secs_f64();
    let ising = g.ising();
    let cuts: Vec<f64> = run.trials.iter().map(|t| if polish { g.cut(&descend(&ising, &t.spins).spins) } else { g.cut(&t.spins) }).collect();
    let best = cuts.iter().cloned().fold(f64::MIN, f64::max);
    let mean = cuts.iter().sum::<f64>() / cuts.len() as f64;
    print!("n={} edges={} {}{:?} dt={} c1={} steps={} trials={}  best {} mean {:.1}", g.n, g.edges.len(), if polish { "+descent " } else { "" }, cfg.variant, cfg.dt, cfg.c1, steps, trials, best, mean);
    if let Some(t) = target {
        let p = cuts.iter().filter(|&&c| c >= t).count() as f64 / cuts.len() as f64;
        print!("  P(≥{t}) {p:.4}  S99 {:?}", step_to_solution(p, f64::from(steps)).map(|s| s.round()));
    }
    println!("  {:.1}s ({:.2} ms/trial-thread)", dt, 1e3 * dt * threads as f64 / trials as f64);
}