//! A quench of the transverse-field Ising chain by TDVP, held at every
//! report to the exact free-fermion solution: the chain starts with every
//! σˣ = +1 and evolves under `−J Σ σᶻσᶻ − h Σ σˣ`.
//!
//! ```text
//! cargo run --release --example tdvp_quench --features quantum_tdvp -- N BOND DT T [h] [cutoff] [threads]
//! ```
use wai_quantum::quantum_dmrg::{Chain, Mps, Op, Scalar, C};
use wai_quantum::quantum_tdvp::{ising_quench, Evolution, TdvpConfig};

fn main() {
    let a: Vec<String> = std::env::args().collect();
    let n: usize = a[1].parse().unwrap();
    let bond: usize = a[2].parse().unwrap();
    let dt: f64 = a[3].parse().unwrap();
    let t_end: f64 = a[4].parse().unwrap();
    let h: f64 = a.get(5).map_or(1.0, |x| x.parse().unwrap());
    let cutoff: f64 = a.get(6).map_or(1e-20, |x| x.parse().unwrap());
    let threads: usize = a.get(7).map_or(TdvpConfig::default().threads, |x| x.parse().unwrap());
    let chain = Chain::ising(n, 1.0, h);
    let plus = Mps::product(&vec![[C { re: 1.0, im: 0.0 }, C { re: 1.0, im: 0.0 }]; n]);
    let cfg = TdvpConfig { dt, max_bond: bond, cutoff, threads, ..TdvpConfig::default() };
    let mut ev = Evolution::new(&chain, &plus, cfg);
    let e0 = ev.energy();
    let steps = (t_end / dt).round() as usize;
    let every = ((0.5 / dt).round() as usize).max(1);
    let t0 = std::time::Instant::now();
    let mut discarded = 0.0;
    println!("ising n={n} h={h} bond={bond} dt={dt} cutoff={cutoff:e}");
    println!("    t   max|Δσˣ|   max|Δσᶻσᶻ|  σˣ(mid)        exact          rate      exact rate  S(mid)  bond  |ΔE|     discarded  seconds");
    for s in 1..=steps {
        discarded += ev.step().discarded;
        if s % every == 0 || s == steps {
            let t = ev.time();
            let ff = ising_quench(n, 1.0, h, t);
            let x = ev.expect(Op::x());
            let zz = ev.bond_expect(Op::z(), Op::z());
            let dx = x.iter().zip(&ff.x).map(|(p, q)| (p - q).abs()).fold(0.0, f64::max);
            let dzz = zz.iter().zip(&ff.zz).map(|(p, q)| (p - q).abs()).fold(0.0, f64::max);
            let echo = ev.overlap(&plus).norm2();
            let rate = -echo.ln() / n as f64;
            let exact_rate = -ff.echo.ln() / n as f64;
            let s_mid = ev.entropies()[n / 2 - 1];
            println!(
                "{t:5.2}  {dx:.2e}   {dzz:.2e}    {:+.10}  {:+.10}  {rate:.8}  {exact_rate:.8}  {s_mid:.4}  {:4}  {:.1e}  {discarded:.1e}  {:.1}",
                x[n / 2],
                ff.x[n / 2],
                ev.max_bond(),
                (ev.energy() - e0).abs(),
                t0.elapsed().as_secs_f64()
            );
        }
    }
}