//! this "example" positions ferromagnetic buffers adjacently and uses an ASYMMETRIC coil winding
//! to couple them. i parameterize the entire setup over a bunch of different factors in order to
//! search for the conditions which maximize energy transfer from the one core to the other.

use coremem::{Driver, mat, meas};
use coremem::geom::Meters;
use coremem::geom::region::{
    self,
    Cube,
    Dilate,
    Intersection,
    InvertedRegion,
    Memoize,
    Spiral,
    SwapYZ,
    Torus,
    Translate,
    Wrap
};
use coremem::mat::{Ferroxcube3R1MH, IsoConductorOr};
use coremem::real::{R32, Real as _};
use coremem::render::CsvRenderer;
use coremem::sim::spirv::{SpirvSim, WgpuBackend};
use coremem::sim::units::{Seconds, Time as _};
use coremem::stim::{CurlVectorField, Exp, ModulatedVectorField, Sinusoid, TimeVaryingExt as _};
use log::{error, info, warn};
use rayon::prelude::*;
use serde::{Deserialize, Serialize};

mod cache;
use cache::SyncDiskCache;

type Mat = IsoConductorOr<f32, Ferroxcube3R1MH>;

#[allow(unused)]
use coremem::geom::{Coord as _, Region as _};

#[allow(unused)]
#[derive(Debug, Copy, Clone, PartialEq, Eq)]
enum PulseType {
    Square,
    Sine,
    /// "pulse" window represents 2 half-lifes
    ExpDecay2x,
}

#[allow(unused)]
/// Return just the extrema of some collection
fn extrema(mut meas: Vec<f32>) -> Vec<f32> {
    let mut i = 0;
    while i + 2 < meas.len() {
        let (prev, cur, next) = (meas[i], meas[i+1], meas[i+2]);
        if (prev <= cur && cur <= next) || (prev >= cur && cur >= next) {
            meas.remove(i+1);
        } else {
            i += 1;
        }
    }
    meas
}

/// Return the (signed) peak magnitude and stable value
fn significa(meas: Vec<f32>) -> (f32, f32) {
    let peak = meas.iter().max_by(|a, b| a.abs().partial_cmp(&b.abs()).unwrap()).copied().unwrap_or_default();
    let stable = meas.last().copied().unwrap_or_default();
    (peak, stable)
}

fn min_max(meas: Vec<f32>) -> (f32, f32) {
    let min = meas.iter().min_by(|a, b| a.partial_cmp(b).unwrap()).copied().unwrap_or_default();
    let max = meas.iter().max_by(|a, b| a.partial_cmp(b).unwrap()).copied().unwrap_or_default();
    (min, max)
}

#[derive(Copy, Clone, Debug, Default, PartialEq, PartialOrd, Serialize, Deserialize)]
struct GeomParams {
    feat_size: f32,
    buffer_xy: f32,
    buffer_z: f32,
    boundary_xy: f32,
    boundary_z: f32,

    ferro_major: f32,
    ferro_minor: f32,
    ferro_buffer: f32, // horizontal space between ferros,
    wire_minor: f32,
    wire_wrap_minor: f32, // 2x wire_wrap_minor + feat_size must be < ferro_buffer
    wire_set_major: f32,
    wire_wrap_dilation: f32,
    wire_wrap_iters: usize,

    wraps1: f32,
    wrap1_coverage: f32,
    wraps2: f32,
    wrap2_coverage: f32,
}

// impl GeomParams {
//     fn key(&self) -> String {
//         let GeomParams {
//             feat_size,
//             buffer_xy,
//             buffer_z,
//             boundary_xy,
//             boundary_z,
//             ferro_major,
//             ferro_minor,
//             ferro_buffer,
//             wire_wrap_minor,
//             wire_set_major,
//             wire_wrap_dilation,
//             wire_wrap_iters,
//             wraps1,
//             wrap1_coverage,
//             wraps2,
//             wrap2_coverage
//         } = self;
//         format!(
//             "{feat_size}um-{buffer_xy}x{buffer_z}buf-{boundary_xy}x{boundary_z}bound-{ferro_major}ferromaj-{ferro_minor}min-{ferro_buffer}buf-{wire_wrap_major}:{wire_set_major}:{wire_wrap_dilation}:{wire_wrap_iters}wrap-{wraps1}:{wraps2}wraps-{wrap1_coverage}:{wrap2_coverage}cov",
//             feat_size=feat_size*1e6,
//             buffer_xy=buffer_xy,
//             buffer_z=buffer_z,
//             boundary_xy=boundary_xy,
//             boundary_z=boundary_z,
//             ferro_major=ferro_major,
//             ferro_minor=ferro_minor,
//             ferro_buffer=ferro_buffer,
//             wire_wrap_minor=wire_wrap_minor,
//             wire_set_major=wire_set_major,
//             wire_wrap_dilation=wire_wrap_dilation,
//             wire_wrap_iters=wire_wrap_iters,
//             wraps1=wraps1,
//             wrap1_coverage=wrap1_coverage,
//             wraps2=wraps2,
//             wrap2_coverage=wrap2_coverage,
//         )
//     }
// }

#[derive(Copy, Clone, Debug)]
struct Params {
    dry_run: bool,
    geom: GeomParams,

    wire_conductivity: f32,

    peak_set_current: f32,
    peak_clock_current: f32,
    set_duration: f32,
    clock_duration: f32,
    clock_type: PulseType,
    pre_time: f32, // how long between set and clock
    post_time: f32, // how long to wait after the clock

    dump_frames: (u64, Option<u64>),
}

#[derive(Clone, Default, Serialize, Deserialize)]
struct Geometries {
    dim: Meters,
    ferro1_region: Torus,
    ferro2_region: Torus,
    set1_region: Torus,
    set2_region: Torus,
    coupling_region: region::Union3<
        Memoize<Dilate<Wrap<Translate<SwapYZ<Intersection<Spiral, Cube>>>>>>,
        Memoize<Dilate<Wrap<Translate<SwapYZ<Intersection<Spiral, InvertedRegion<Cube>>>>>>>,
        region::Union3<Cube, Cube, region::Union4<Cube, Cube, Cube, Cube>>
    >,
    coupling_wire_top: Cube,
    coupling_wire_bot: Cube,
    wrap1_len: f32,
    wrap2_len: f32,
}

/// computed measurements which get written to disk for later, manual (or grep-based) analysis.
/// because we only write these (except for the Debug impl reading them to write to disk),
/// rustc thinks all the fields are dead.
#[derive(Clone, Debug, Default)]
struct Results {
    #[allow(dead_code)]
    m1_peak: f32,
    #[allow(dead_code)]
    m2_peak: f32,
    #[allow(dead_code)]
    m1_stable: f32,
    #[allow(dead_code)]
    m2_stable: f32,
    #[allow(dead_code)]
    h1_peak: f32,
    #[allow(dead_code)]
    h2_max: f32,
    #[allow(dead_code)]
    h2_min: f32,
    #[allow(dead_code)]
    h1_stable: f32,
    #[allow(dead_code)]
    h2_stable: f32,
    #[allow(dead_code)]
    iset_min: f32,
    #[allow(dead_code)]
    iset_max: f32,
    #[allow(dead_code)]
    icoupling_peak: f32,
    #[allow(dead_code)]
    peak_m_ratio: f32,
    #[allow(dead_code)]
    stable_m_ratio: f32,
    /// m2_stable divided by m1_peak. i.e. "amplification"
    #[allow(dead_code)]
    m2_stable_m1_peak: f32,
    #[allow(dead_code)]
    t: f32,
}

fn derive_geometries(p: GeomParams) -> Option<Geometries> {
    use std::f32::consts::PI;
    let feat_sizes = Meters::new(p.feat_size, p.feat_size, p.feat_size);
    let width = 4.0*p.ferro_major + 2.0*(p.buffer_xy + p.boundary_xy + p.wire_set_major + p.wire_minor) + p.ferro_buffer;
    let height = 2.0*(p.ferro_major + p.ferro_minor + 4.0*p.wire_wrap_minor + 12.0*p.feat_size + p.buffer_xy + p.boundary_xy);
    let depth = 2.0*(p.wire_set_major.max(4.0*p.wire_wrap_minor + p.ferro_minor + p.feat_size) + p.wire_minor + p.buffer_z + p.boundary_z);
    let dim = Meters::new(width, height, depth);

    let ferro1_center = Meters::new(
        p.buffer_xy + p.boundary_xy + p.wire_set_major + p.wire_minor + p.ferro_major,
        p.buffer_xy + p.boundary_xy + p.ferro_major + p.ferro_minor + 4.0*p.wire_wrap_minor + 12.0*p.feat_size,
        0.5*depth,
        // buffer_z + boundary_z + wire_set_major + wire_minor
    );
    let ferro2_center = ferro1_center + Meters::new(2.0*p.ferro_major + p.ferro_buffer, 0.0, 0.0);
    let ferro_center = (ferro1_center + ferro2_center)*0.5;

    // reserve the left/right locations for the SET wires.
    let set1_center = ferro1_center - Meters::new_x(p.ferro_major);
    let set2_center = ferro2_center + Meters::new_x(p.ferro_major);

    let ferro1_region = Torus::new_xy(ferro1_center, p.ferro_major, p.ferro_minor);
    let ferro2_region = Torus::new_xy(ferro2_center, p.ferro_major, p.ferro_minor);

    let set1_region = Torus::new_xz(set1_center, p.wire_set_major, p.wire_minor);
    let set2_region = Torus::new_xz(set2_center, p.wire_set_major, p.wire_minor);

    let wrap1_rate = 2.0*p.wrap1_coverage/p.wraps1;
    let coupling_region1 = Memoize::new(Dilate::new(
        Wrap::new_about(
            Translate::new(
                SwapYZ::new(region::and(
                    Spiral::new(p.ferro_minor + 2.0*p.wire_wrap_minor + p.feat_size, p.wire_wrap_minor, wrap1_rate),
                    Cube::new(Meters::new(-1.0, -1.0, -p.wrap1_coverage), Meters::new(1.0, 1.0, p.wrap1_coverage))
                )),
                ferro1_center + Meters::new(1.0*p.ferro_major, 0.0, 0.0),
            ),
            1.0, // one half-rev => y=1.0
            ferro1_center,
        ),
        p.wire_wrap_dilation,
        p.wire_wrap_dilation / (p.wire_wrap_iters as f32),
    ));

    let wrap2_rate = 2.0*p.wrap2_coverage/p.wraps2;
    let coupling_region2 = Memoize::new(Dilate::new(
        Wrap::new_about(
            Translate::new(
                SwapYZ::new(region::and_not(
                    Spiral::new(p.ferro_minor + 2.0*p.wire_wrap_minor + p.feat_size, p.wire_wrap_minor, wrap2_rate),
                    Cube::new(Meters::new(-1.0, -1.0, -1.0 + p.wrap2_coverage), Meters::new(1.0, 1.0, 1.0 - p.wrap2_coverage))
                )),
                ferro2_center + Meters::new_x(p.ferro_major),
            ),
            1.0, // one half-rev => y=1.0
            ferro2_center,
        ),
        p.wire_wrap_dilation,
        p.wire_wrap_dilation / (p.wire_wrap_iters as f32),
    ));

    let coupling_wire_top = Cube::new_centered(
        ferro_center - Meters::new_y(p.ferro_major + 4.0*p.wire_wrap_minor + 12.0*p.feat_size),
        Meters::new(p.ferro_buffer + 4.0*p.ferro_major + 4.0*p.feat_size, 2.0*p.feat_size, 2.0*p.feat_size)
    );
    let coupling_wire_bot = Cube::new_centered(
        ferro_center + Meters::new_y(p.ferro_major + 4.0*p.wire_wrap_minor + 12.0*p.feat_size),
        Meters::new(p.ferro_buffer + 4.0*p.ferro_major + 4.0*p.feat_size, 2.0*p.feat_size, 2.0*p.feat_size)
    );
    let wrap1_top = ferro1_center + Meters::new_x(p.ferro_major).rotate_z(-p.wrap1_coverage*PI);
    let wrap1_bot = ferro1_center + Meters::new_x(p.ferro_major).rotate_z(p.wrap1_coverage*PI);
    let wrap2_top = ferro2_center + Meters::new_x(p.ferro_major).rotate_z((1.0+p.wrap2_coverage)*PI);
    let wrap2_bot = ferro2_center + Meters::new_x(p.ferro_major).rotate_z((1.0-p.wrap2_coverage)*PI);
    let coupling_stub_top_left = Cube::new_including_negatives(
        wrap1_top + feat_sizes*2.0,
        wrap1_top.with_y(coupling_wire_top.bot()) - feat_sizes*2.0,
    );
    let coupling_stub_bot_left = Cube::new_including_negatives(
        wrap1_bot - feat_sizes*2.0,
        wrap1_bot.with_y(coupling_wire_bot.top()) + feat_sizes*2.0,
    );
    let coupling_stub_top_right = Cube::new_including_negatives(
        wrap2_top + feat_sizes*2.0,
        wrap2_top.with_y(coupling_wire_top.bot()) - feat_sizes*2.0,
    );
    let coupling_stub_bot_right = Cube::new_including_negatives(
        wrap2_bot - feat_sizes*2.0,
        wrap2_bot.with_y(coupling_wire_bot.top()) + feat_sizes*2.0,
    );
    let coupling_stubs = region::Union::new4(
        coupling_stub_top_left.clone(),
        coupling_stub_top_right.clone(),
        coupling_stub_bot_left.clone(),
        coupling_stub_bot_right.clone(),
    );

    let coupling_wires = region::Union::new3(
        coupling_wire_top.clone(),
        coupling_wire_bot.clone(),
        coupling_stubs.clone(),
    );

    let coupling_region = region::Union::new3(
        coupling_region1.clone(),
        coupling_region2.clone(),
        coupling_wires.clone(),
    );

    let wrap1_with_coupling = region::Union::new2(
        coupling_region1.clone(), coupling_wires.clone()
    );
    let wrap2_with_coupling = region::Union::new2(
        coupling_region2.clone(), coupling_wires.clone()
    );

    // show that the coupling top/bot wires are connected through the wrapping
    if !region::is_connected(
        &wrap1_with_coupling, 
        coupling_wire_top.center(),
        coupling_wire_bot.center(),
        p.feat_size,
    ) {
        warn!("wrap1 not connected for params: {:?}", p);
        return None;
    }
    if !region::is_connected(
        &wrap2_with_coupling, 
        coupling_wire_top.center(),
        coupling_wire_bot.center(),
        p.feat_size,
    ) {
        warn!("wrap2 not connected for params: {:?}", p);
        return None;
    }

    let wrap1_len = region::distance_to(
        &wrap1_with_coupling,
        coupling_stub_top_left.center().to_index(p.feat_size),
        coupling_stub_bot_left.center().to_index(p.feat_size),
        p.feat_size,
    ).unwrap();
    let wrap2_len = region::distance_to(
        &wrap2_with_coupling,
        coupling_stub_top_right.center().to_index(p.feat_size),
        coupling_stub_bot_right.center().to_index(p.feat_size),
        p.feat_size,
    ).unwrap();
    info!("wrap lengths: {}, {} (wraps: {}, {})", wrap1_len, wrap2_len, p.wraps1, p.wraps2);

    Some(Geometries {
        dim,
        ferro1_region,
        ferro2_region,
        set1_region,
        set2_region,
        coupling_region,
        coupling_wire_top,
        coupling_wire_bot,
        wrap1_len,
        wrap2_len,
    })
}

fn ensure_out_dir(id: u32) -> String {
    let base = format!("out/applications/buffer_proto5/{}", id);
    let _ = std::fs::create_dir_all(&base);
    base
}

fn run_sim(id: u32, p: Params, g: Geometries) -> Results {
    info!("run_sim {}: {:?}", id, p);

    let m_to_um = |m: f32| (m * 1e6).round() as u32;
    let feat_vol = p.geom.feat_size * p.geom.feat_size * p.geom.feat_size;

    info!("boundary: {}um; {}um", m_to_um(p.geom.boundary_xy), m_to_um(p.geom.boundary_z));
    info!("size: {:?}", g.dim);
    info!("ferro1: {:?}", g.ferro1_region.center());
    info!("ferro2: {:?}", g.ferro2_region.center());

    let base = ensure_out_dir(id);
    let prefix = format!(
        "{}/{}-{}-{}setmA-{}setps-{}clkmA-{}clkps-{}um-{}ferromaj-{}:{}wraps-{}:{}cov-{:?}clk",
        base,
        id,
        *g.dim.to_index(p.geom.feat_size),
        (p.peak_set_current * 1e3).round() as i64,
        (p.set_duration * 1e12).round() as i64,
        (p.peak_clock_current * 1e3).round() as i64,
        (p.clock_duration * 1e12).round() as i64,
        (p.geom.feat_size * 1e6).round() as i64,
        p.geom.ferro_major,
        p.geom.wraps1,
        p.geom.wraps2,
        p.geom.wrap1_coverage,
        p.geom.wrap2_coverage,
        p.clock_type,
    );

    let mut driver = Driver::new(SpirvSim::<f32, Mat, WgpuBackend>::new(
        g.dim.to_index(p.geom.feat_size), p.geom.feat_size
    ));
    if !driver.add_state_file(&*format!("{}/state.bc", prefix), 4000) {
        // mu_r=881.33, starting at H=25 to H=75.
        let ferro_mat = mat::Ferroxcube3R1MH::new();
        // let ferro_mat = mat::db::conductor(wire_conductivity);
        let wire_mat = mat::IsomorphicConductor::new(p.wire_conductivity);

        driver.fill_region(&g.ferro1_region, ferro_mat);
        driver.fill_region(&g.ferro2_region, ferro_mat);
        driver.fill_region(&g.set1_region, wire_mat);
        driver.fill_region(&g.set2_region, wire_mat);
        driver.fill_region(&g.coupling_region, wire_mat);

        driver.add_classical_boundary(Meters::new(p.geom.boundary_xy, p.geom.boundary_xy, p.geom.boundary_z));

        // assert!(driver.test_region_filled(&g.ferro1_region, ferro_mat));
        // assert!(driver.test_region_filled(&g.ferro2_region, ferro_mat));
        // assert!(driver.test_region_filled(&g.set1_region, wire_mat));
        // assert!(driver.test_region_filled(&g.set2_region, wire_mat));
        // assert!(driver.test_region_filled(&g.coupling_region, wire_mat));
    } else {
        info!("loaded state file: skipping geometry calculations");
    }

    let add_drive_sine_pulse = |driver: &mut Driver<f32, _, _>, region: &Torus, start: f32, duration: f32, amp: f32| {
        let wave = Sinusoid::from_wavelength(duration * 2.0)
            .half_cycle()
            .scaled(amp)
            .shifted(start);
        driver.add_stimulus(ModulatedVectorField::new(
            CurlVectorField::new(region.clone()),
            wave,
        ));
    };

    let add_drive_square_pulse = |driver: &mut Driver<f32, _, _>, region: &Torus, start: f32, duration: f32, amp: f32| {
        let wave = amp.gated(start, start+duration);
        driver.add_stimulus(ModulatedVectorField::new(
            CurlVectorField::new(region.clone()),
            wave,
        ));
    };

    let add_drive_exp_pulse = |driver: &mut Driver<f32, _, _>, region: &Torus, start: f32, duration: f32, amp: f32| {
        let wave = Exp::new_at(amp, start, 0.5*duration);
        driver.add_stimulus(ModulatedVectorField::new(
            CurlVectorField::new(region.clone()),
            wave,
        ));
    };

    // step function: "permanently" increase the current by `amp`.
    let _add_drive_step = |driver: &mut Driver<f32, _, _>, region: &Torus, start: f32, amp: f32| {
        add_drive_square_pulse(driver, region, start, 1.0 /* effectively infinite duration */, amp);
    };

    let add_drive_pulse = |ty: PulseType, driver: &mut Driver<f32, _, _>, region: &Torus, start: f32, duration: f32, amp: f32| {
        match ty {
            PulseType::Square => add_drive_square_pulse(driver, region, start, duration, amp),
            PulseType::Sine => add_drive_sine_pulse(driver, region, start, duration, amp),
            PulseType::ExpDecay2x => add_drive_exp_pulse(driver, region, start, duration, amp),
        }
    };

    // J=\sigma E
    // dJ/dt = \sigma dE/dT
    // dE/dt = dJ/dt / \sigma
    // dE/dt = dI/dt / (A*\sigma)
    // if I = k*sin(w t) then dE/dt = k*w cos(w t) / (A*\sigma)
    // i.e. dE/dt is proportional to I/(A*\sigma), multiplied by w (or, divided by wavelength)
    let peak_set = p.peak_set_current / feat_vol / (g.set1_region.cross_section() * p.wire_conductivity);
    let peak_clock = p.peak_clock_current / feat_vol / (g.set1_region.cross_section() * p.wire_conductivity);

    // SET cores
    add_drive_sine_pulse(&mut driver, &g.set1_region, 0.01*p.set_duration, p.set_duration, -peak_set);
    add_drive_sine_pulse(&mut driver, &g.set2_region, 0.01*p.set_duration, p.set_duration, peak_set);

    // CLEAR core1
    add_drive_pulse(p.clock_type, &mut driver, &g.set1_region, p.set_duration + p.pre_time, p.clock_duration, peak_clock);
    // add_drive_step(&mut driver, &set1_region, set_duration + pre_time, peak_clock);

    let duration = Seconds(p.set_duration + p.pre_time + p.clock_duration + p.post_time)
        .to_frame(driver.timestep())
        .round_up(32000);

    driver.add_measurement(meas::Volume::new("mem1", g.ferro1_region.clone()));
    driver.add_measurement(meas::MagneticLoop::new("mem1", g.ferro1_region.clone()));

    driver.add_measurement(meas::Volume::new("mem2", g.ferro2_region.clone()));
    driver.add_measurement(meas::MagneticLoop::new("mem2", g.ferro2_region.clone()));

    driver.add_measurement(meas::CurrentLoop::new("set1", g.set1_region.clone()));
    driver.add_measurement(meas::Power::new("set1", g.set1_region.clone()));
    driver.add_measurement(meas::CurrentLoop::new("set2", g.set2_region.clone()));

    // driver.add_measurement(meas::CurrentLoop::new("coupling1", coupling_region1.clone()));
    // driver.add_measurement(meas::CurrentLoop::new("coupling2", coupling_region2.clone()));
    driver.add_measurement(meas::Current::new("couplingtop", g.coupling_wire_top.clone()));
    driver.add_measurement(meas::Current::new("couplingbot", g.coupling_wire_bot.clone()));

    if p.dry_run {
        info!("bailing (dry run): {}", prefix);
        return Results::default();
    }

    let _ = std::fs::create_dir_all(&prefix);

    let (frame_freq, frame_limit) = p.dump_frames;
    driver.add_serializer_renderer(&*format!("{}/frame-", prefix), frame_freq, frame_limit);
    let meas_csv = format!("{}/meas.csv", prefix);
    let meas_sparse_csv = format!("{}/meas-sparse.csv", prefix);
    driver.add_csv_renderer(&*meas_csv, 400, None);
    driver.add_csv_renderer(&*meas_sparse_csv, 8000, None);

    driver.set_steps_per_stimulus(20);
    driver.step_until(duration);

    let (m1_peak, m1_stable) = significa(CsvRenderer::new(&*meas_sparse_csv).read_column_as_f32("M(mem1)"));
    let (m2_peak, m2_stable) = significa(CsvRenderer::new(&*meas_sparse_csv).read_column_as_f32("M(mem2)"));
    let (h1_peak, h1_stable) = significa(CsvRenderer::new(&*meas_sparse_csv).read_column_as_f32("H(mem1)"));
    let (_h2_peak, h2_stable) = significa(CsvRenderer::new(&*meas_sparse_csv).read_column_as_f32("H(mem2)"));
    let (h2_min, h2_max) = min_max(CsvRenderer::new(&*meas_sparse_csv).read_column_as_f32("H(mem2)"));
    let (iset_min, iset_max) = min_max(CsvRenderer::new(&*meas_sparse_csv).read_column_as_f32("I(set1)"));
    let (icoupling_peak, _icoupling_stable) = significa(CsvRenderer::new(&*meas_sparse_csv).read_column_as_f32("Imag/cell(couplingtop)"));
    let res = Results {
        m1_peak,
        m2_peak,
        m1_stable,
        m2_stable,
        h1_peak,
        h2_min,
        h2_max,
        h1_stable,
        h2_stable,
        iset_min,
        iset_max,
        icoupling_peak,
        peak_m_ratio: m2_peak / m1_peak,
        stable_m_ratio: m2_stable / m1_stable,
        m2_stable_m1_peak: m2_stable / m1_peak,
        t: driver.time(),
    };
    std::fs::write(
        format!("{}/results.txt", prefix),
        format!("{:#?}\n", res),
    ).unwrap();
    info!("completed sim: {}", prefix);
    res
}


fn main() {
    let um = 1e-6;
    let ns = 1e-9;
    coremem::init_logging();
    let i = 63;
    let dry_run = false;
    let mut variants = Vec::new();
    // old variants:
    // for wrap1_density in [1.0, 0.5, 0.2] {
    // for wrap2_density in [1.0, 0.5, 0.2, 0.0] {
    // for (wrap1_cov, wrap2_cov) in [(0.8, 0.8), (0.8, 0.25), (0.25, 0.8)] {
    let clock_domain = [
        (25600.0, 5.0 * ns),
        // (1600.0, 5.0 * ns), // low relevance for large cores
        // (1600.0, 1.0 * ns), // very poor perf (0.05 m2_stable_m1_peak)
        (25600.0, 1.0 * ns),
        // (6400.0, 1.0 * ns), // low relevance for large cores
        // (51200.0, 1.0 * ns),
        (102400.0, 1.0 * ns),

        (51200.0, 5.0 * ns),
        (12800.0, 5.0 * ns),
        // TODO: RE-ENABLE
        // (6400.0, 25.0 * ns), // suspended because costly
        (12800.0, 25.0 * ns), // suspended because costly
        // (25600.0, 25.0 * ns), // suspended because costly
        // (6400.0, 5.0 * ns),  // low relevance for large cores (for >= 9mm rad, <0.5 m2_stable_m1_peak)

        // (409600.0, 1.0 * ns), // I(set1) shows significant underdamping => bad m2_stable_m1_peak
        // (400.0, 25.0 * ns), // poor perf (0.28 m2_stable_m1_peak)
        // (1600.0, 25.0 * ns), // mediocre perf (0.65 m2_stable_m1_peak)

        // (1600.0, 100.0 * ns),
        // (6400.0, 100.0 * ns),
        // (400.0, 100.0 * ns),
        // (400.0, 25.0 * ns),
        // (100.0, 100e-9),
        // (25.0, 25e-9),
        // (25.0, 200e-9), // if we see a good trend here, scale this further
        // (100.0, 25e-9),
        // (1600.0, 25e-9),
    ];
    let ferro_majors = [
        11000.0 * um,
        // 1360.0 * um,
        9000.0 * um,
        // 13500.0 * um, // fails geom test?
        6000.0 * um,
        4080.0 * um,
        // 2320.0 * um,
        // 3200.0 * um,
        // 1680.0 * um
    ];
    let wrap2_densities = [
        -0.3,
        -0.5,
        -0.15,
        -0.05,
        // 0.3,
        // 0.2,
        // 0.15,
        // 0.05,
        // 0.0,
    ];
    let clock_types = [
        // PulseType::Sine,
        PulseType::ExpDecay2x,
        // PulseType::Square
    ];
    let post_times = [
        // 10ns is enough to see m2 peak
        // 10.0 * ns, // completed
        // 20.0 * ns,
        // 50ns is 4 extra half-lifes for a 25ns exp decay
        50.0 * ns,
        // 125ns is a midpoint. we should be able to guess m2 stabilization *trends*
        125.0 * ns,
        // 275ns is enough to see m2 stabilize
        275.0 * ns,
    ];
    let wrap1_density = 1.0;
    let (wrap1_cov, wrap2_cov) = (0.8, 0.8);
    for post_time in post_times {
        for ferro_major in ferro_majors {
            for wrap2_density in wrap2_densities {
                for (peak_clock_current, clock_duration) in clock_domain {
                    for clock_type in clock_types {
                        variants.push((
                            peak_clock_current,
                            clock_duration,
                            post_time,
                            clock_type,
                            ferro_major,
                            wrap1_cov,
                            wrap2_cov,
                            wrap1_density,
                            wrap2_density
                        ));
                    }
                }
            }
        }
    }
    info!(
        "evaluating {} variants ({} time increments of {} primitives)",
        variants.len(),
        post_times.len(),
        variants.len() / post_times.len(),
    );

    let geom_cache = SyncDiskCache::new_with_supplier(
        &format!("{}/.geom_cache", ensure_out_dir(i)),
        |geom: &GeomParams| derive_geometries(geom.clone())
    );

    for (peak_clock_current, clock_duration, post_time, clock_type, ferro_major, wrap1_coverage, wrap2_coverage, wrap1_density, wrap2_density) in variants {
        info!("{}A/{}s {}s {}m {}:{}cov {}:{}density", peak_clock_current, clock_duration, post_time, ferro_major, wrap1_coverage, wrap2_coverage, wrap1_density, wrap2_density);
        let base_params = Params {
            dry_run,
            geom: GeomParams {
                feat_size: 40e-6f32,

                buffer_xy: 160e-6,
                buffer_z: 160e-6,
                boundary_xy: 320e-6,
                boundary_z: 320e-6,

                ferro_major,
                ferro_minor: 60e-6,
                ferro_buffer: 1320e-6,
                wire_minor: 40e-6,
                wire_wrap_minor: 50e-6,
                wire_set_major: 200e-6,
                wire_wrap_dilation: 50e-6,
                wire_wrap_iters: 3,

                wraps1: 8.0,
                wraps2: 8.0,
                wrap1_coverage,
                wrap2_coverage,
            },

            wire_conductivity: 5e6f32,

            peak_set_current: 60.0,
            peak_clock_current,
            set_duration: 0e-9,
            clock_duration,
            clock_type,
            pre_time: 1e-9,
            post_time,

            dump_frames: (256000, Some(257000)),
        };

        let wraps1_choices: Vec<_> = (-120..120)
            .into_par_iter()
            .filter_map(|wraps1| {
                let params = GeomParams {
                    wraps1: (wraps1 * 4) as f32,
                    ..base_params.geom
                };
                let geoms = geom_cache.get_or_insert_from_supplier(params.clone())?;
                Some((params, geoms))
            })
            .collect();
        let wraps1_max = wraps1_choices
            .iter()
            .max_by_key(|(_p, geoms)| R32::from_primitive(geoms.wrap1_len))
            .cloned()
            .unwrap_or_default()
            .0.wraps1;
        let wraps1_target = wraps1_max * wrap1_density;
        let wraps1 = wraps1_choices
            .iter()
            .min_by_key(|(p, _g)| R32::from_primitive((p.wraps1 - wraps1_target).abs()))
            .cloned()
            .unwrap_or_default()
            .0.wraps1;

        let wraps2_choices: Vec<_> = (-120..120)
            .into_par_iter()
            .filter_map(|wraps2| {
                let params = GeomParams {
                    wraps2: (wraps2 * 4) as f32,
                    ..base_params.geom
                };
                let geoms = geom_cache.get_or_insert_from_supplier(params.clone())?;
                Some((params, geoms))
            })
            .collect();
        let wraps2_max = wraps2_choices
            .iter()
            .max_by_key(|(_p, geoms)| R32::from_primitive(geoms.wrap2_len))
            .cloned()
            .unwrap_or_default()
            .0.wraps2;
        let wraps2_target = wraps2_max * wrap2_density;
        let wraps2 = wraps2_choices
            .iter()
            .min_by_key(|(p, _g)| R32::from_primitive((p.wraps2 - wraps2_target).abs()))
            .cloned()
            .unwrap_or_default()
            .0.wraps2;

        let mut params = base_params;
        params.geom.wraps1 = wraps1;
        params.geom.wraps2 = wraps2;
        match geom_cache.get_or_insert_from_supplier(params.geom.clone()) {
            Some(geoms) => {
                run_sim(i, params, geoms);
            },
            None => error!("skipping sim because no valid geometry: {:?}", params),
        }
    }

    info!("done");
}