Region: remove Clone from the trait, and also parameterize everything

i didn't necessarily *want* to parameterize it all,
but it turned out to be easier to do that than to force all users to
workaround the lack of Clone.

.gitignore

@@ -1,3 +1,2 @@
 out/
 target/
-*.swp

Cargo.lock

@@ -317,7 +317,6 @@ dependencies = [
  "crossterm",
  "csv",
  "dashmap",
- "dyn-clone",
  "env_logger",
  "float_eq",
  "font8x8",
@@ -339,7 +338,6 @@ dependencies = [
  "spirv_backend",
  "spirv_backend_runner",
  "threadpool",
- "typetag",
  "wgpu",
  "y4m",
 ]
@@ -502,16 +500,6 @@ dependencies = [
  "memchr",
 ]
 
-[[package]]
-name = "ctor"
-version = "0.1.22"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f877be4f7c9f246b183111634f75baa039715e3f46ce860677d3b19a69fb229c"
-dependencies = [
- "quote",
- "syn",
-]
-
 [[package]]
 name = "cty"
 version = "0.2.2"
@@ -565,12 +553,6 @@ dependencies = [
  "libc",
 ]
 
-[[package]]
-name = "dyn-clone"
-version = "1.0.6"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "140206b78fb2bc3edbcfc9b5ccbd0b30699cfe8d348b8b31b330e47df5291a5a"
-
 [[package]]
 name = "either"
 version = "1.7.0"
@@ -590,15 +572,6 @@ dependencies = [
  "termcolor",
 ]
 
-[[package]]
-name = "erased-serde"
-version = "0.3.21"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "81d013529d5574a60caeda29e179e695125448e5de52e3874f7b4c1d7360e18e"
-dependencies = [
- "serde",
-]
-
 [[package]]
 name = "exr"
 version = "1.4.2"
@@ -802,17 +775,6 @@ dependencies = [
  "wasm-bindgen",
 ]
 
-[[package]]
-name = "ghost"
-version = "0.1.5"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "b93490550b1782c589a350f2211fff2e34682e25fed17ef53fc4fa8fe184975e"
-dependencies = [
- "proc-macro2",
- "quote",
- "syn",
-]
-
 [[package]]
 name = "gif"
 version = "0.11.4"
@@ -1021,16 +983,6 @@ dependencies = [
  "web-sys",
 ]
 
-[[package]]
-name = "inventory"
-version = "0.3.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e0257e268c91daba296499206db5dd5a058875936bec3d8429b5f3e20ec9dc9a"
-dependencies = [
- "ctor",
- "ghost",
-]
-
 [[package]]
 name = "itertools"
 version = "0.10.3"
@@ -2236,30 +2188,6 @@ version = "1.15.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "dcf81ac59edc17cc8697ff311e8f5ef2d99fcbd9817b34cec66f90b6c3dfd987"
 
-[[package]]
-name = "typetag"
-version = "0.2.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "36afde536989c95b0c35d2b6720618ed268b271a002fa1ebd97e821cf9ff28c6"
-dependencies = [
- "erased-serde",
- "inventory",
- "once_cell",
- "serde",
- "typetag-impl",
-]
-
-[[package]]
-name = "typetag-impl"
-version = "0.2.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "2c935ee7eb11a684e863c6a3332ff17caa2a46fc855719a90c6493aee213852b"
-dependencies = [
- "proc-macro2",
- "quote",
- "syn",
-]
-
 [[package]]
 name = "unicode-ident"
 version = "1.0.1"

crates/applications/buffer_proto5/src/main.rs

@@ -3,7 +3,20 @@
 //! search for the conditions which maximize energy transfer from the one core to the other.
 
 use coremem::{Driver, mat, meas};
-use coremem::geom::{region, Cube, Dilate, Memoize, Meters, Spiral, SwapYZ, Torus, Translate, Wrap};
+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;
@@ -14,7 +27,6 @@ use log::{error, info, warn};
 use serde::{Deserialize, Serialize};
 
 mod cache;
-
 use cache::DiskCache;
 
 type Mat = IsoConductorOr<f32, Ferroxcube3R1MH>;
@@ -149,7 +161,11 @@ struct Geometries {
     ferro2_region: Torus,
     set1_region: Torus,
     set2_region: Torus,
-    coupling_region: region::Union,
+    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,
@@ -285,29 +301,29 @@ fn derive_geometries(p: GeomParams) -> Option<Geometries> {
         wrap2_bot - feat_sizes*2.0,
         wrap2_bot.with_y(coupling_wire_bot.top()) + feat_sizes*2.0,
     );
-    let coupling_stubs = region::Union::new()
+    let coupling_stubs = region::Union::new4(
-        .with(coupling_stub_top_left.clone())
+        coupling_stub_top_left.clone(),
-        .with(coupling_stub_top_right.clone())
+        coupling_stub_top_right.clone(),
-        .with(coupling_stub_bot_left.clone())
+        coupling_stub_bot_left.clone(),
-        .with(coupling_stub_bot_right.clone())
+        coupling_stub_bot_right.clone(),
-    ;
+    );
 
-    let coupling_wires = region::Union::new()
+    let coupling_wires = region::Union::new3(
-        .with(coupling_wire_top.clone())
+        coupling_wire_top.clone(),
-        .with(coupling_wire_bot.clone())
+        coupling_wire_bot.clone(),
-        .with(coupling_stubs.clone())
+        coupling_stubs.clone(),
-    ;
+    );
 
-    let coupling_region = region::Union::new()
+    let coupling_region = region::Union::new3(
-        .with(coupling_region1.clone())
+        coupling_region1.clone(),
-        .with(coupling_region2.clone())
+        coupling_region2.clone(),
-        .with(coupling_wires.clone())
+        coupling_wires.clone(),
-    ;
+    );
 
-    let wrap1_with_coupling = region::union(
+    let wrap1_with_coupling = region::Union::new2(
         coupling_region1.clone(), coupling_wires.clone()
     );
-    let wrap2_with_coupling = region::union(
+    let wrap2_with_coupling = region::Union::new2(
         coupling_region2.clone(), coupling_wires.clone()
     );
 
@@ -429,8 +445,6 @@ fn run_sim(id: u32, p: Params, g: Geometries) -> Results {
         driver.add_stimulus(CurlStimulus::new(
             region.clone(),
             wave.clone(),
-            region.center(),
-            region.axis()
         ));
     };
 
@@ -439,8 +453,6 @@ fn run_sim(id: u32, p: Params, g: Geometries) -> Results {
         driver.add_stimulus(CurlStimulus::new(
             region.clone(),
             wave.clone(),
-            region.center(),
-            region.axis()
         ));
     };
 
@@ -449,8 +461,6 @@ fn run_sim(id: u32, p: Params, g: Geometries) -> Results {
         driver.add_stimulus(CurlStimulus::new(
             region.clone(),
             wave.clone(),
-            region.center(),
-            region.axis()
         ));
     };
 

crates/applications/multi_core_inverter/src/main.rs

@@ -42,134 +42,227 @@ use coremem::meas;
 use coremem::real::{self, Real as _};
 use coremem::sim::spirv::{self, SpirvSim};
 use coremem::sim::units::Seconds;
-use coremem::stim::{CurlStimulus, Gated, Sinusoid1, TimeVarying as _};
+use coremem::stim::{CurlStimulus, Exp, Gated, Sinusoid1, TimeVarying as _};
 use coremem::Driver;
 
 type R = real::R32;
 type Mat = IsoConductorOr<R, Ferroxcube3R1MH>;
 // type Backend = spirv::CpuBackend;
 type Backend = spirv::WgpuBackend;
+type Sim = SpirvSim::<R, Mat, Backend>;
+
+#[derive(Clone, Copy)]
+enum DriveType {
+    HoldHigh,
+    ReleaseHigh,
+    HoldLow,
+    ReleaseLow,
+    Float,
+}
+
+impl DriveType {
+    fn direction(&self) -> i32 {
+        use DriveType::*;
+        match self {
+            HoldHigh | ReleaseHigh => 1,
+            HoldLow | ReleaseLow => -1,
+            Float => 0,
+        }
+    }
+    fn is_hold(&self) -> bool {
+        match self {
+            Self::HoldHigh | Self::HoldLow => true,
+            _ => false,
+        }
+    }
+    fn is_release(&self) -> bool {
+        match self {
+            Self::ReleaseHigh | Self::ReleaseLow => true,
+            _ => false,
+        }
+    }
+}
 
 fn main() {
     coremem::init_logging();
-    coremem::init_debug();
+    // coremem::init_debug();
     let um = |n| n as f32 * 1e-6;
     // let ns = |n| n as f32 * 1e-9;
     let ps = |n| n as f32 * 1e-12;
     let feat_size = um(10);
 
     let input_magnitude = 1.0e9;
-    let clock_phase_duration = ps(1000);
+    let clock_phase_duration = ps(2000);
 
+    // 's' = core (ferromagnetic part)
     let s_major = um(160);
     let s_minor = um(30);
+    // 'io' = drive/control wire
     let io_major = um(80);
     let io_minor = um(30);
     let coupling_major = um(130);
     let coupling_minor = um(30);
 
+    // coords for core 'n'
     let sx = |n| um((n+1) * 400);
     let sy = um(400);
     let sz = um(280);
+    // x-coord for the loop coupling core 'n' to core 'n+1'
     let couplingx = |n| sx(n) + um(200);
 
-    let sim_bounds = Meters::new(sx(4), sy * 2.0, sz * 2.0);
+    let sim_bounds = |num_cores| Meters::new(sx(num_cores), sy * 2.0, sz * 2.0);
     let sim_padding = Meters::new(um(80), um(80), um(80));
 
-    let drive0 = Torus::new_xz(Meters::new(sx(0) - s_major, sy, sz), io_major, io_minor);
+    // core '0' gets an external input in place of its coupling loop
-    let sense3 = Torus::new_xz(Meters::new(sx(3) + s_major, sy, sz), io_major, io_minor);
+    let input0 = Torus::new_xz(Meters::new(sx(0) - s_major, sy, sz), io_major, io_minor);
+    // the last core gets an external output in place of its coupling loop
+    let sense = |n| Torus::new_xz(Meters::new(sx(n) + s_major, sy, sz), io_major, io_minor);
+    // control loop for core n (alternately called "drive" loop)
     let ctl = |n| Torus::new_yz(Meters::new(sx(n), sy + s_major, sz), io_major, io_minor);
     let s = |n| Torus::new_xy(Meters::new(sx(n), sy, sz), s_major, s_minor);
     // coupling(n) is the wire which couples core n into core n+1
     let coupling = |n| Torus::new_xz(Meters::new(couplingx(n), sy, sz), coupling_major, coupling_minor);
 
-    let input = |region: &Torus, cycle: u32, direction: i32| {
+    let control_signal = |driver: &mut Driver<Sim>, region: &Torus, cycle: u32, ty: DriveType| {
         let area = region.cross_section();
-        let amp = direction as f32 * input_magnitude / area;
+        let amp = ty.direction() as f32 * input_magnitude / area;
         let start = clock_phase_duration * cycle as f32;
-        let wave = Gated::new(amp, start, start + clock_phase_duration);
+        if ty.is_hold() {
+            // simple square wave:
+            let wave = Gated::new(amp, start, start + clock_phase_duration);
+            let stim = CurlStimulus::new(
+                region.clone(),
+                wave.clone(),
+            );
+            driver.add_stimulus(stim);
+        } else if ty.is_release() {
+            // decaying exponential wave:
+            let wave = Exp::new_at(amp, start, 0.125 * clock_phase_duration /* half-life */);
+            let stim = CurlStimulus::new(
+                region.clone(),
+                wave.clone(),
+            );
+            driver.add_stimulus(stim);
+        } // else: Float
+        // sine wave (untested):
         // let wave = Sinusoid1::from_wavelength(direction as f32 * input_magnitude / area, clock_phase_duration * 2.0)
         //     .half_cycle()
         //     .shifted(clock_phase_duration * cycle as f32);
-        CurlStimulus::new(
-            region.clone(),
-            wave.clone(),
-            region.center(),
-            region.axis(),
-        )
     };
 
+    //////// define the control signals/transitions
+    // each row N denotes the drive currents at clock cycle N.
+    // each col M denotes the drive current at core M.
+    use DriveType::*;
+    // let s0_init_state = HoldHigh; // logic high
+    let s0_init_state = HoldLow; // logic low
+    let s0_rel_state = ReleaseLow;
+    let drive_map = [
+        // charge each device to '1'
+        [s0_init_state,HoldHigh,     HoldHigh,     HoldHigh,     HoldHigh],
+        // this is when we'd ordinarily open S0 for write
+        [s0_rel_state, HoldHigh,     HoldHigh,     HoldHigh,     HoldHigh],
+
+        // this is when S0 would ordinarily receive a write
+        [Float,        HoldHigh,     HoldHigh,     HoldHigh,     HoldHigh],
+        // open S1 for write
+        [Float,        ReleaseHigh,  HoldHigh,     HoldHigh,     HoldHigh],
+
+        // write S0' => S1
+        [HoldLow,      Float,        HoldHigh,     HoldHigh,     HoldHigh],
+        // open S2 for write
+        [HoldLow,      Float,        ReleaseHigh,  HoldHigh,     HoldHigh],
+
+        // write S1' => S2
+        [HoldLow,      HoldLow,      Float,        HoldHigh,     HoldHigh],
+        // open S3 for write
+        [HoldLow,      HoldLow,      Float,        ReleaseHigh,  HoldHigh],
+
+        // write S2' => S3; RESET S0
+        [HoldHigh,     HoldLow,      HoldLow,      Float,        HoldHigh],
+        // open S4 for write
+        [HoldHigh,     HoldLow,      HoldLow,      Float,        ReleaseHigh],
+
+        // write S3' => S4; RESET S1
+        [HoldHigh,     HoldHigh,     HoldLow,      HoldLow,      Float],
+        // open S0 for write
+        [ReleaseHigh,  HoldHigh,     HoldLow,      HoldLow,      Float],
+
+        // write S4' => output; RESET S2
+        [Float,        HoldHigh,     HoldHigh,     HoldLow,      HoldLow],
+        // open S1 for write
+        [Float,        ReleaseHigh,  HoldHigh,     HoldLow,      HoldLow],
+    ];
+
     let wire_mat = IsomorphicConductor::new(1e6f32.cast::<R>());
     // let ferro_mat = wire_mat;
     let ferro_mat = Ferroxcube3R1MH::new();
 
-    let mut driver = Driver::new(SpirvSim::<R, Mat, Backend>::new(
+    let num_cores = drive_map[0].len();
-        sim_bounds.to_index(feat_size), feat_size,
+    let last_core = num_cores - 1;
+
+    let mut driver = Driver::new(Sim::new(
+        sim_bounds(num_cores).to_index(feat_size), feat_size,
     ));
     driver.add_classical_boundary_explicit::<R, _>(sim_padding);
 
     //////// create the wires and toroids
-    driver.fill_region(&drive0, wire_mat);
+    driver.fill_region(&input0, wire_mat);
-    driver.fill_region(&sense3, wire_mat);
+    driver.fill_region(&sense(last_core), wire_mat);
-    for core in 0..4 {
+    for core in 0..num_cores {
         driver.fill_region(&s(core), ferro_mat);
         driver.fill_region(&ctl(core), wire_mat);
-        if core != 3 {
+        if core != last_core {
             driver.fill_region(&coupling(core), wire_mat);
         }
     }
 
     //////// monitor some measurements
-    // driver.add_measurement(meas::CurrentLoop::new("drv0", drive0.clone()));
+    // driver.add_measurement(meas::CurrentLoop::new("input0", input0.clone()));
-    for core in 0..4 {
+    for core in 0..num_cores {
         driver.add_measurement(meas::CurrentLoop::new(
             &format!("drive{}", core),
             ctl(core),
         ));
     }
-    for core in 0..3 {
+    for core in 0..num_cores {
-        driver.add_measurement(meas::CurrentLoop::new(
+        let name = format!("sense{}", core);
-            &format!("sense{}", core),
+        if core != last_core {
-            coupling(core),
+            driver.add_measurement(meas::CurrentLoop::new(
-        ));
+                &name, coupling(core),
+            ));
+        } else {
+            driver.add_measurement(meas::CurrentLoop::new(
+                &name, sense(core),
+            ));
+        }
     }
-    for core in 0..4 {
+    for core in 0..num_cores {
         driver.add_measurement(meas::MagneticLoop::new(
             &format!("state{}", core),
             s(core),
         ));
     }
-    driver.add_measurement(meas::CurrentLoop::new("sense3", sense3.clone()));
 
-    //////// create the stimuli
+
-    // CTL{n} effectively leads CTL{n-1}
+    let cycles = drive_map.len() as u32;
-    // or: at time t, CTL{n} is at cycle[t+n]
+    let duration = Seconds(clock_phase_duration * (cycles + 3) as f32);
-    // where cycle[t] is defined by CTL[0](t):
+
-    //   0, +Vdd, +Vdd, +Vdd
-    // TODO: this is wrong (as is the diagram in the blog)! CTL0, being an inverter,
-    //   needs -Vdd to recharge to +polarization
-    let cycles = 1;
-    let duration = Seconds(clock_phase_duration * (cycles + 2) as f32);
     for cycle in 0..cycles {
-        for core in 0..1 { // TODO: core 0
+        for core in 0..num_cores {
-            let dir = 1;
+            // current polarity of the drive wires is inverted in this model v.s. the blog article,
-            //let dir = if (cycle+core) % 4 == 0 {
+            // so invert our currents in order to achieve the magnetic states of the article.
-            //    0
+            let sig = drive_map[cycle as usize][core as usize];
-            //} else {
+            control_signal(&mut driver, &ctl(core), cycle, sig);
-            //    1
-            //};
-            if dir != 0 {
-                // micro opt/safety: don't place zero-magnitude stimuli
-                driver.add_stimulus(input(&ctl(core), cycle, dir));
-            }
         }
     }
 
-    let prefix = "out/applications/multi_core_inverter/2/";
+    let prefix = "out/applications/multi_core_inverter/11/";
     let _ = std::fs::create_dir_all(&prefix);
-    // driver.add_state_file(&*format!("{}state.bc", prefix), 9600);
+    driver.add_state_file(&*format!("{}state.bc", prefix), 6400);
-    driver.add_serializer_renderer(&*format!("{}frame-", prefix), 400, None);
+    driver.add_serializer_renderer(&*format!("{}frame-", prefix), 3200, None);
-    driver.add_csv_renderer(&*format!("{}meas.csv", prefix), 100, None);
+    driver.add_csv_renderer(&*format!("{}meas-detailed.csv", prefix), 100, None);
+    driver.add_csv_renderer(&*format!("{}meas.csv", prefix), 800, None);
     driver.set_steps_per_stim(100);
 
     driver.step_until(duration);

crates/applications/sr_latch/src/main.rs

@@ -84,8 +84,6 @@ fn main() {
         driver.add_stimulus(CurlStimulus::new(
             region.clone(),
             wave.clone(),
-            region.center(),
-            region.axis()
         ));
     };
 

crates/coremem/Cargo.toml

@@ -15,7 +15,6 @@ common_macros = "0.1"  # MIT or Apache 2.0
 crossterm = "0.24"  # MIT
 csv = "1.1"  # MIT or Unlicense
 dashmap = "5.3"  # MIT
-dyn-clone = "1.0"  # MIT or Apache 2.0
 env_logger = "0.9"  # MIT or Apache 2.0
 float_eq = "1.0"  # MIT or Apache 2.0
 font8x8 = "0.3"  # MIT
@@ -33,7 +32,6 @@ rand = "0.8"  # MIT or Apache 2.0
 rayon = "1.5"  # MIT or Apache 2.0
 serde = "1.0"  # MIT or Apache 2.0
 threadpool = "1.8"  # MIT or Apache 2.0
-typetag = "0.2"  # MIT or Apache 2.0
 y4m = "0.7"  # MIT
 
 wgpu = "0.12"

crates/coremem/src/diagnostics.rs

@@ -0,0 +1,111 @@
+use std::sync::{Arc, Mutex};
+use std::time::{Duration, Instant};
+
+/// this is just a big dumb bag of perf-related metrics,
+/// gathered with the intention of identifying areas for optimization
+pub struct Diagnostics {
+    frames_completed: u64,
+    time_spent_stepping: Duration,
+    time_spent_on_stimuli: Duration,
+    time_spent_prepping_render: Duration,
+    time_spent_blocked_on_render: Duration,
+    time_spent_rendering: Duration,
+    start_time: Instant,
+}
+
+#[derive(Clone, Default)]
+pub struct SyncDiagnostics(Arc<Mutex<Diagnostics>>);
+
+impl Default for Diagnostics {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+impl Diagnostics {
+    pub fn new() -> Self {
+        Self {
+            frames_completed: 0,
+            time_spent_stepping: Default::default(),
+            time_spent_on_stimuli: Default::default(),
+            time_spent_prepping_render: Default::default(),
+            time_spent_blocked_on_render: Default::default(),
+            time_spent_rendering: Default::default(),
+            start_time: Instant::now(),
+        }
+    }
+
+    pub fn format(&self) -> String {
+        let step_time = self.time_spent_stepping.as_secs_f64();
+        let stim_time = self.time_spent_on_stimuli.as_secs_f64();
+        let render_time = self.time_spent_rendering.as_secs_f64();
+        let render_prep_time = self.time_spent_prepping_render.as_secs_f64();
+        let block_time = self.time_spent_blocked_on_render.as_secs_f64();
+        let overall_time = self.start_time.elapsed().as_secs_f64();
+        let fps = (self.frames_completed as f64) / overall_time;
+        format!("fps: {:6.2} (sim: {:.1}s, stim: {:.1}s, [render: {:.1}s], blocked: {:.1}s, render_prep: {:.1}s, other: {:.1}s)",
+            fps,
+            step_time,
+            stim_time,
+            render_time,
+            block_time,
+            render_prep_time,
+            overall_time - step_time /*- stim_time*/ - block_time - render_prep_time,
+        )
+    }
+}
+
+impl SyncDiagnostics {
+    pub fn new() -> Self {
+        Self(Arc::new(Mutex::new(Diagnostics::new())))
+    }
+    pub fn format(&self) -> String {
+        self.0.lock().unwrap().format()
+    }
+    /// measure the duration of some arbitrary chunk of code.
+    /// used internally.
+    pub fn measure<R, F: FnOnce() -> R>(f: F) -> (Duration, R) {
+        let start = Instant::now();
+        let r = f();
+        (start.elapsed(), r)
+    }
+
+    /// record the duration of the sim step operation.
+    pub fn instrument_step<R, F: FnOnce() -> R>(&self, frames: u64, f: F) -> R {
+        let (elapsed, ret) = Self::measure(f);
+        let mut me = self.0.lock().unwrap();
+        me.time_spent_stepping += elapsed;
+        me.frames_completed += frames as u64;
+        ret
+    }
+
+    /// record the duration spent preparing for render (i.e. cloning stuff and moving it into a
+    /// render pool).
+    pub fn instrument_render_prep<R, F: FnOnce() -> R>(&self, f: F) -> R {
+        let (elapsed, ret) = Self::measure(f);
+        self.0.lock().unwrap().time_spent_prepping_render += elapsed;
+        ret
+    }
+
+    /// record the duration actually spent doing CPU render work
+    pub fn instrument_render_cpu_side<R, F: FnOnce() -> R>(&self, f: F) -> R {
+        let (elapsed, ret) = Self::measure(f);
+        self.0.lock().unwrap().time_spent_rendering += elapsed;
+        ret
+    }
+
+    /// record the duration spent blocking the simulation because the render queue is full.
+    pub fn instrument_render_blocked<R, F: FnOnce() -> R>(&self, f: F) -> R{
+        let (elapsed, ret) = Self::measure(f);
+        self.0.lock().unwrap().time_spent_blocked_on_render += elapsed;
+        ret
+    }
+
+    pub fn instrument_stimuli<R, F: FnOnce() -> R>(&self, f: F) -> R {
+        let (elapsed, ret) = Self::measure(f);
+        self.0.lock().unwrap().time_spent_on_stimuli += elapsed;
+        ret
+    }
+
+}
+

crates/coremem/src/driver.rs

@@ -1,3 +1,4 @@
+use crate::diagnostics::SyncDiagnostics;
 use crate::geom::{Coord, Index, Meters, Region};
 use crate::mat;
 use crate::meas::{self, AbstractMeasurement};
@@ -10,9 +11,9 @@ use crate::stim::{self, AbstractStimulus};
 use log::{info, trace};
 use serde::{Deserialize, Serialize};
 use std::path::PathBuf;
-use std::sync::{Arc, Mutex};
+use std::sync::Arc;
 use std::sync::mpsc::{sync_channel, SyncSender, Receiver};
-use std::time::{Duration, Instant};
+use std::time::Instant;
 use threadpool::ThreadPool;
 
 pub struct Driver<S> {
@@ -21,32 +22,23 @@ pub struct Driver<S> {
     // TODO: use Rayon's thread pool?
     render_pool: ThreadPool,
     render_channel: (SyncSender<()>, Receiver<()>),
-    // TODO: split out into a Diagnostics struct
-    time_spent_stepping: Duration,
-    time_spent_on_stimuli: Duration,
-    time_spent_prepping_render: Duration,
-    time_spent_blocked_on_render: Duration,
-    time_spent_rendering: Arc<Mutex<Duration>>,
     measurements: Vec<Arc<dyn AbstractMeasurement<S>>>,
     stimuli: StimuliAdapter,
-    start_time: Instant,
-    last_diag_time: Instant,
     /// simulation end time
     sim_end_time: Option<Frame>,
+    diag: SyncDiagnostics,
+    last_diag_time: Instant,
 }
 
 impl<S: AbstractSim> Driver<S> {
-    pub fn new(state: S) -> Self {
+    pub fn new(mut state: S) -> Self {
+        let diag = SyncDiagnostics::new();
+        state.use_diagnostics(diag.clone());
         Self {
             state,
             renderer: Arc::new(MultiRenderer::new()),
             render_pool: ThreadPool::new(3),
             render_channel: sync_channel(0),
-            time_spent_stepping: Default::default(),
-            time_spent_on_stimuli: Default::default(),
-            time_spent_prepping_render: Default::default(),
-            time_spent_blocked_on_render: Default::default(),
-            time_spent_rendering: Default::default(),
             measurements: vec![
                 Arc::new(meas::Time),
                 Arc::new(meas::Meta),
@@ -54,9 +46,9 @@ impl<S: AbstractSim> Driver<S> {
                 Arc::new(meas::Power::world()),
             ],
             stimuli: StimuliAdapter::new(),
-            start_time: Instant::now(),
-            last_diag_time: Instant::now(),
             sim_end_time: None,
+            diag,
+            last_diag_time: Instant::now(),
         }
     }
 
@@ -135,6 +127,7 @@ impl<S: AbstractSim + Send + Sync + Serialize + for<'a> Deserialize<'a> + 'stati
         let ser = render::SerializerRenderer::new(state_file);
         let loaded = ser.try_load().map(|s| {
             self.state = s.state;
+            self.state.use_diagnostics(self.diag.clone());
         }).is_some();
         self.add_renderer(ser, state_file, snapshot_frequency, None);
         loaded
@@ -143,26 +136,26 @@ impl<S: AbstractSim + Send + Sync + Serialize + for<'a> Deserialize<'a> + 'stati
 
 impl<S: AbstractSim + Clone + Default + Send + 'static> Driver<S> {
     fn render(&mut self) {
-        let prep_start = Instant::now();
+        self.diag.instrument_render_prep(|| {
-        let their_state = self.state.clone();
+            let diag_handle = self.diag.clone();
-        let their_measurements = self.measurements.clone();
+            let their_state = self.state.clone();
-        let renderer = self.renderer.clone();
+            let their_measurements = self.measurements.clone();
-        let time_spent_rendering = self.time_spent_rendering.clone();
+            let renderer = self.renderer.clone();
-        let sender = self.render_channel.0.clone();
+            let sender = self.render_channel.0.clone();
-        self.render_pool.execute(move || {
+            self.render_pool.execute(move || {
-            // unblock the main thread (this limits the number of renders in-flight at any time
+                // unblock the main thread (this limits the number of renders in-flight at any time
-            sender.send(()).unwrap();
+                sender.send(()).unwrap();
-            trace!("render begin");
+                trace!("render begin");
-            let start_time = Instant::now();
+                diag_handle.instrument_render_cpu_side(|| {
-            let meas: Vec<&dyn AbstractMeasurement<S>> = their_measurements.iter().map(|m| &**m).collect();
+                    let meas: Vec<&dyn AbstractMeasurement<S>> = their_measurements.iter().map(|m| &**m).collect();
-            renderer.render(&their_state, &*meas, Default::default());
+                    renderer.render(&their_state, &*meas, Default::default());
-            *time_spent_rendering.lock().unwrap() += start_time.elapsed();
+                });
-            trace!("render end");
+                trace!("render end");
+            });
+        });
+        self.diag.instrument_render_blocked(|| {
+            self.render_channel.1.recv().unwrap();
         });
-        self.time_spent_prepping_render += prep_start.elapsed();
-        let block_start = Instant::now();
-        self.render_channel.1.recv().unwrap();
-        self.time_spent_blocked_on_render += block_start.elapsed();
     }
     /// Return the number of steps actually stepped
     fn step_at_most(&mut self, at_most: u32) -> u32 {
@@ -183,37 +176,22 @@ impl<S: AbstractSim + Clone + Default + Send + 'static> Driver<S> {
             can_step += 1;
         }
         trace!("step begin");
-        let start_time = Instant::now();
+        self.diag.instrument_step(can_step as u64, || {
-        self.state.step_multiple(can_step, &self.stimuli);
+            self.state.step_multiple(can_step, &self.stimuli);
-        self.time_spent_stepping += start_time.elapsed();
+        });
         trace!("step end");
         if self.last_diag_time.elapsed().as_secs_f64() >= 5.0 {
             self.last_diag_time = Instant::now();
             let step = self.state.step_no();
-            let step_time = self.time_spent_stepping.as_secs_f64();
+            let diagstr = self.diag.format();
-            let stim_time = self.time_spent_on_stimuli.as_secs_f64();
-            let render_time = self.time_spent_rendering.lock().unwrap().as_secs_f64();
-            let render_prep_time = self.time_spent_prepping_render.as_secs_f64();
-            let block_time = self.time_spent_blocked_on_render.as_secs_f64();
-            let overall_time = self.start_time.elapsed().as_secs_f64();
-            let fps = (self.state.step_no() as f64) / overall_time;
             let sim_time = self.state.time() as f64;
             let percent_complete = match self.sim_end_time {
                 Some(t) => format!("[{:.1}%] ", 100.0 * self.state.time() / *t.to_seconds(self.timestep())),
                 None => "".to_owned(),
             };
             info!(
-                "{}t={:.2e} frame {:06} fps: {:6.2} (sim: {:.1}s, stim: {:.1}s, [render: {:.1}s], blocked: {:.1}s, render_prep: {:.1}s, other: {:.1}s)",
+                "{}t={:.2e} frame {:06} {}",
-                percent_complete,
+                percent_complete, sim_time, step, diagstr
-                sim_time,
-                step,
-                fps,
-                step_time,
-                stim_time,
-                render_time,
-                block_time,
-                render_prep_time,
-                overall_time - step_time - stim_time - block_time - render_prep_time,
             );
         }
         can_step as u32

crates/coremem/src/geom/mod.rs

@@ -6,7 +6,22 @@ mod units;
 pub use line::Line2d;
 pub use polygon::Polygon2d;
 pub use region::{
-    Cube, CylinderZ, Dilate, InvertedRegion, Memoize, Region, Sphere, Spiral, SwapXZ, SwapYZ, Torus, Translate, Union, WorldRegion, Wrap
+    Cube,
+    CylinderZ,
+    Dilate,
+    HasCrossSection,
+    InvertedRegion,
+    Memoize,
+    Region,
+    Sphere,
+    Spiral,
+    SwapXZ,
+    SwapYZ,
+    Torus,
+    Translate,
+    Union,
+    WorldRegion,
+    Wrap,
 };
 pub use units::{Coord, Meters, OrdMeters, Index};
 

crates/coremem/src/geom/region/mod.rs

@@ -1,6 +1,7 @@
 use crate::geom::{Coord, Meters, OrdMeters};
 
-use dyn_clone::{self, DynClone};
+use coremem_cross::vec::Vec3;
+
 use rayon::prelude::*;
 use serde::{Serialize, Deserialize};
 use std::collections::BTreeSet;
@@ -10,22 +11,27 @@ mod primitives;
 pub use primitives::*;
 
 
-#[typetag::serde(tag = "type")]
+pub trait Region: Send + Sync {
-pub trait Region: Send + Sync + DynClone {
     fn contains(&self, p: Meters) -> bool;
 }
-dyn_clone::clone_trait_object!(Region);
 
-pub fn and<T1: Region + 'static, T2: Region + 'static>(r1: T1, r2: T2) -> Intersection {
+/// some (volume) which has a tangent vector everywhere inside/on it.
-    Intersection::new().and(r1).and(r2)
+/// for example, a cylinder has tangents everywhere except its axis.
+/// the returned vector should represent the area of the cross section.
+pub trait HasCrossSection {
+    fn cross_section_normal(&self, p: Meters) -> Vec3<f32>;
 }
 
-pub fn and_not<T1: Region + 'static, T2: Region + 'static>(r1: T1, r2: T2) -> Intersection {
+pub fn and<T1: Region + 'static, T2: Region + 'static>(r1: T1, r2: T2) -> Intersection<T1, T2> {
+    Intersection::new2(r1, r2)
+}
+
+pub fn and_not<T1: Region + 'static, T2: Region + 'static>(r1: T1, r2: T2) -> Intersection<T1, InvertedRegion<T2>> {
     and(r1, InvertedRegion::new(r2))
 }
 
-pub fn union<T1: Region + 'static, T2: Region + 'static>(r1: T1, r2: T2) -> Union {
+pub fn union<T1: Region + 'static, T2: Region + 'static>(r1: T1, r2: T2) -> Union<T1, T2> {
-    Union::new().with(r1).with(r2)
+    Union::new2(r1, r2)
 }
 
 /// returns true if there's a path (via the cardinal directions) from p0 to p1 within this region.
@@ -67,137 +73,124 @@ pub fn distance_to<R: Region, C: Coord>(r: &R, p0: C, p1: C, feat_size: f32) ->
 }
 
 /// Region describing the entire simulation space
-#[derive(Copy, Clone, Serialize, Deserialize)]
+#[derive(Copy, Clone, Default, Serialize, Deserialize)]
 pub struct WorldRegion;
 
-#[typetag::serde]
 impl Region for WorldRegion {
     fn contains(&self, _: Meters) -> bool {
         true
     }
 }
 
-#[derive(Clone, Serialize, Deserialize)]
+#[derive(Clone, Default, Serialize, Deserialize)]
-pub struct InvertedRegion(Box<dyn Region>);
+pub struct InvertedRegion<R>(R);
 
-impl InvertedRegion {
+impl<R> InvertedRegion<R> {
-    pub fn new<R: Region + 'static>(r: R) -> Self {
+    pub fn new(r: R) -> Self {
-        Self(Box::new(r))
+        Self(r)
     }
 }
 
-#[typetag::serde]
+impl<R: Region> Region for InvertedRegion<R> {
-impl Region for InvertedRegion {
     fn contains(&self, p: Meters) -> bool {
         !self.0.contains(p)
     }
 }
 
 #[derive(Clone, Default, Serialize, Deserialize)]
-pub struct Union(Vec<Box<dyn Region>>);
+pub struct Union<R1, R2>(R1, R2);
 
-impl Union {
+pub type Union3<R1, R2, R3> = Union<Union<R1, R2>, R3>;
-    pub fn new() -> Self {
+pub type Union4<R1, R2, R3, R4> = Union<Union3<R1, R2, R3>, R4>;
-        Self(Vec::new())
+
+impl<R1, R2> Union<R1, R2> {
+    pub fn with<R: Region>(self, r: R) -> Union<Self, R> {
+        Union::new2(self, r)
     }
-    pub fn new_with<R: Region + 'static>(r: R) -> Self {
+    pub fn new2(r1: R1, r2: R2) -> Self {
-        Self::new().with(r)
+        Self(r1, r2)
     }
-    pub fn with<R: Region + 'static>(self, r: R) -> Self {
+    pub fn new3<R3: Region>(r1: R1, r2: R2, r3: R3) -> Union<Self, R3> {
-        self.with_box(Box::new(r))
+        Union::new2(r1, r2).with(r3)
     }
-    pub fn with_box(mut self, r: Box<dyn Region>) -> Self {
+    pub fn new4<R3: Region, R4: Region>(r1: R1, r2: R2, r3: R3, r4: R4) -> Union<Union<Self, R3>, R4> {
-        self.0.push(r);
+        Union::new2(r1, r2).with(r3).with(r4)
-        self
     }
 }
 
-#[typetag::serde]
+impl<R1: Region, R2: Region> Region for Union<R1, R2> {
-impl Region for Union {
     fn contains(&self, p: Meters) -> bool {
-        self.0.iter().any(|r| r.contains(p))
+        self.0.contains(p) || self.1.contains(p)
     }
 }
 
-#[derive(Clone, Serialize, Deserialize)]
+#[derive(Clone, Default, Serialize, Deserialize)]
-pub struct Intersection(Vec<Box<dyn Region>>);
+pub struct Intersection<R1, R2>(R1, R2);
 
 
-impl Intersection {
+impl<R1, R2> Intersection<R1, R2> {
-    pub fn new() -> Self {
+    pub fn and<R3: Region>(self, r: R3) -> Intersection<Self, R3> {
-        Self(Vec::new())
+        Intersection::new2(self, r)
     }
-    pub fn new_with<R: Region + 'static>(r: R) -> Self {
+    pub fn new2(r1: R1, r2: R2) -> Self {
-        Self::new().and(r)
+        Self(r1, r2)
-    }
-    pub fn and<R: Region + 'static>(self, r: R) -> Self {
-        self.and_box(Box::new(r))
-    }
-    pub fn and_box(mut self, r: Box<dyn Region>) -> Self {
-        self.0.push(r);
-        self
     }
 }
 
-#[typetag::serde]
+impl<R1: Region, R2: Region> Region for Intersection<R1, R2> {
-impl Region for Intersection {
     fn contains(&self, p: Meters) -> bool {
-        self.0.iter().all(|r| r.contains(p))
+        self.0.contains(p) && self.1.contains(p)
     }
 }
 
-#[derive(Clone, Serialize, Deserialize)]
+#[derive(Clone, Default, Serialize, Deserialize)]
-pub struct Translate {
+pub struct Translate<R> {
-    inner: Box<dyn Region>,
+    inner: R,
     shift: Meters,
 }
 
-impl Translate {
+impl<R> Translate<R> {
-    pub fn new<T: Region + 'static>(inner: T, shift: Meters) -> Self {
+    pub fn new(inner: R, shift: Meters) -> Self {
-        Self { inner: Box::new(inner), shift }
+        Self { inner, shift }
     }
 }
 
-#[typetag::serde]
+impl<R: Region> Region for Translate<R> {
-impl Region for Translate {
     fn contains(&self, p: Meters) -> bool {
         self.inner.contains(p - self.shift)
     }
 }
 
-#[derive(Clone, Serialize, Deserialize)]
+#[derive(Clone, Default, Serialize, Deserialize)]
-pub struct SwapXZ {
+pub struct SwapXZ<R> {
-    inner: Box<dyn Region>,
+    inner: R,
 }
 
-impl SwapXZ {
+impl<R> SwapXZ<R> {
-    pub fn new<T: Region + 'static>(inner: T) -> Self {
+    pub fn new(inner: R) -> Self {
-        Self { inner: Box::new(inner) }
+        Self { inner }
     }
 }
 
-#[typetag::serde]
+impl<R: Region> Region for SwapXZ<R> {
-impl Region for SwapXZ {
     fn contains(&self, p: Meters) -> bool {
         let p = Meters::new(p.z(), p.y(), p.z());
         self.inner.contains(p)
     }
 }
 
-
+#[derive(Clone, Default, Serialize, Deserialize)]
-#[derive(Clone, Serialize, Deserialize)]
+pub struct SwapYZ<R> {
-pub struct SwapYZ {
+    inner: R,
-    inner: Box<dyn Region>,
 }
 
-impl SwapYZ {
+impl<R> SwapYZ<R> {
-    pub fn new<T: Region + 'static>(inner: T) -> Self {
+    pub fn new(inner: R) -> Self {
-        Self { inner: Box::new(inner) }
+        Self { inner }
     }
 }
 
-#[typetag::serde]
+impl<R: Region> Region for SwapYZ<R> {
-impl Region for SwapYZ {
     fn contains(&self, p: Meters) -> bool {
         let mapped = Meters::new(p.x(), p.z(), p.y());
         self.inner.contains(mapped)
@@ -210,20 +203,20 @@ impl Region for SwapYZ {
 /// the resulting region is mapped onto the original region y=[0, y_max]. x is just the radius
 /// so that (0, 0) is mapped to (0, 0), and (1, 0) is mapped to (1, 0) and (0, 1) is mapped to
 /// (1, 0.5*y_max) and (-5, 0) is mapped to (5, 0.5*y_max).
-#[derive(Clone, Serialize, Deserialize)]
+#[derive(Clone, Default, Serialize, Deserialize)]
-pub struct Wrap {
+pub struct Wrap<R> {
-    inner: Box<dyn Region>,
+    inner: R,
     y_max: f32,
     about: Meters,
 }
 
-impl Wrap {
+impl<R> Wrap<R> {
-    pub fn new<T: Region + 'static>(inner: T, y_max: f32) -> Self {
+    pub fn new(inner: R, y_max: f32) -> Self {
         Self::new_about(inner, y_max, Meters::new(0.0, 0.0, 0.0))
     }
 
-    pub fn new_about<T: Region + 'static>(inner: T, y_max: f32, about: Meters) -> Self {
+    pub fn new_about(inner: R, y_max: f32, about: Meters) -> Self {
-        Self { inner: Box::new(inner), y_max, about }
+        Self { inner, y_max, about }
     }
 
     fn map(&self, p: Meters) -> Meters {
@@ -235,28 +228,26 @@ impl Wrap {
     }
 }
 
-#[typetag::serde]
+impl<R: Region> Region for Wrap<R> {
-impl Region for Wrap {
     fn contains(&self, p: Meters) -> bool {
         self.inner.contains(self.map(p))
     }
 }
 
-#[derive(Clone, Serialize, Deserialize)]
+#[derive(Clone, Default, Serialize, Deserialize)]
-pub struct Dilate {
+pub struct Dilate<R> {
-    inner: Box<dyn Region>,
+    inner: R,
     rad: f32,
     res: f32,
 }
 
-impl Dilate {
+impl<R> Dilate<R> {
-    pub fn new<T: Region + 'static>(inner: T, rad: f32, res: f32) -> Self {
+    pub fn new(inner: R, rad: f32, res: f32) -> Self {
-        Self { inner: Box::new(inner), rad, res }
+        Self { inner, rad, res }
     }
 }
 
-#[typetag::serde]
+impl<R: Region> Region for Dilate<R> {
-impl Region for Dilate {
     fn contains(&self, p: Meters) -> bool {
         let rad_iters = (self.rad / self.res).ceil() as i32;
         let rad_range = -rad_iters..=rad_iters;
@@ -279,24 +270,23 @@ impl Region for Dilate {
     }
 }
 
-#[derive(Clone, Serialize, Deserialize)]
+#[derive(Clone, Default, Serialize, Deserialize)]
-pub struct Memoize {
+pub struct Memoize<R> {
     #[serde(skip)]
     lut: Arc<dashmap::DashMap<OrdMeters, bool>>,
-    inner: Box<dyn Region>,
+    inner: R,
 }
 
-impl Memoize {
+impl<R> Memoize<R> {
-    pub fn new<R: Region + 'static>(inner: R) -> Self {
+    pub fn new(inner: R) -> Self {
         Self {
             lut: Arc::new(dashmap::DashMap::new()),
-            inner: Box::new(inner),
+            inner,
         }
     }
 }
 
-#[typetag::serde]
+impl<R: Region> Region for Memoize<R> {
-impl Region for Memoize {
     fn contains(&self, p: Meters) -> bool {
         *self.lut.entry(OrdMeters(p)).or_insert_with(|| self.inner.contains(p))
     }
@@ -307,7 +297,7 @@ mod test {
     use super::*;
     use float_eq::assert_float_eq;
 
-    fn assert_map(w: &Wrap, from: Meters, to: Meters) {
+    fn assert_map<R>(w: &Wrap<R>, from: Meters, to: Meters) {
         let mapped = w.map(from);
         assert_float_eq!(mapped.x(), to.x(), abs <= 0.01);
         assert_float_eq!(mapped.y(), to.y(), abs <= 0.01);

crates/coremem/src/geom/region/primitives.rs

@@ -6,9 +6,9 @@ use serde::{Serialize, Deserialize};
 use std::fmt::{self, Display};
 use std::ops::Range;
 
-use super::Region;
+use super::{HasCrossSection, Region};
 
-#[derive(Copy, Clone, Serialize, Deserialize)]
+#[derive(Copy, Clone, Default, Serialize, Deserialize)]
 pub struct CylinderZ {
     center: Vec2<f32>,
     radius: f32,
@@ -24,7 +24,6 @@ impl CylinderZ {
     }
 }
 
-#[typetag::serde]
 impl Region for CylinderZ {
     fn contains(&self, p: Meters) -> bool {
         p.xy().distance_sq(self.center) <= self.radius * self.radius
@@ -77,7 +76,6 @@ impl Torus {
     }
 }
 
-#[typetag::serde]
 impl Region for Torus {
     fn contains(&self, p: Meters) -> bool {
         // a torus is the set of all points < distance `r` from the circle of radius `R`,
@@ -104,7 +102,16 @@ impl Region for Torus {
     }
 }
 
-#[derive(Copy, Clone, Serialize, Deserialize)]
+impl HasCrossSection for Torus {
+    fn cross_section_normal(&self, coord: Meters) -> Vec3<f32> {
+        let axis = self.axis();
+        let to_coord = *coord - *self.center();
+        // this creates a normal which always points "counter-clockwise" along the shape
+        axis.cross(to_coord)
+    }
+}
+
+#[derive(Copy, Clone, Default, Serialize, Deserialize)]
 pub struct Sphere {
     center: Meters,
     rad: f32,
@@ -119,7 +126,6 @@ impl Sphere {
     }
 }
 
-#[typetag::serde]
 impl Region for Sphere {
     fn contains(&self, p: Meters) -> bool {
         p.distance_sq(*self.center) < self.rad * self.rad
@@ -228,7 +234,6 @@ impl Cube {
     }
 }
 
-#[typetag::serde]
 impl Region for Cube {
     fn contains(&self, p: Meters) -> bool {
         self.x_range().contains(&p.x()) &&
@@ -238,7 +243,7 @@ impl Region for Cube {
 }
 
 /// a Spiral traces out a circle on the xy plane as z increases.
-#[derive(Copy, Clone, Serialize, Deserialize)]
+#[derive(Copy, Clone, Default, Serialize, Deserialize)]
 pub struct Spiral {
     /// radius of the spiral
     major: f32,
@@ -257,7 +262,6 @@ impl Spiral {
     }
 }
 
-#[typetag::serde]
 impl Region for Spiral {
     fn contains(&self, p: Meters) -> bool {
         let revs = p.z() / self.period;

crates/coremem/src/lib.rs

@@ -7,6 +7,7 @@
 
 use log::info;
 
+mod diagnostics;
 pub mod driver;
 pub mod geom;
 pub mod meas;

crates/coremem/src/meas.rs

@@ -1,13 +1,13 @@
-use crate::geom::{Meters, Region, Torus, WorldRegion};
+use crate::geom::{HasCrossSection, Meters, Region, Torus, WorldRegion};
 use crate::real::{Real as _, ToFloat as _};
-use crate::cross::vec::Vec3;
+use crate::cross::vec::{Vec3, Vec3u};
 use crate::sim::AbstractSim;
-use indexmap::IndexMap;
 use serde::{Serialize, Deserialize};
+use std::ops::AddAssign;
 
+// TODO: do we really need both Send and Sync?
 pub trait AbstractMeasurement<S>: Send + Sync {
-    fn eval(&self, state: &S) -> String;
+    fn key_value(&self, state: &S) -> Vec<Measurement>;
-    fn key_value(&self, state: &S) -> IndexMap<String, String>;
 }
 
 pub fn as_dyn_measurements<S, M: AbstractMeasurement<S>>(meas: &[M]) -> Vec<&dyn AbstractMeasurement<S>> {
@@ -15,20 +15,104 @@ pub fn as_dyn_measurements<S, M: AbstractMeasurement<S>>(meas: &[M]) -> Vec<&dyn
 }
 
 
-/// create one IndexMap out of several measurements
+/// combine several measurements
-pub fn eval_multiple_kv<S>(state: &S, meas: &[&dyn AbstractMeasurement<S>]) -> IndexMap<String, String> {
+pub fn eval_multiple<S>(state: &S, meas: &[&dyn AbstractMeasurement<S>]) -> Vec<Measurement> {
-    let mut r = IndexMap::new();
+    meas.into_iter().flat_map(|m| m.key_value(state).into_iter()).collect()
-    for m in meas {
-        let other = m.key_value(state);
-        r.extend(other.into_iter());
-    }
-    r
 }
 
-pub fn eval_to_vec<S>(state: &S, meas: &[&dyn AbstractMeasurement<S>]) -> Vec<Evaluated> {
+#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
-    eval_multiple_kv(state, meas).into_iter().map(|(k, v)| {
+pub enum MeasurementValue {
-        Evaluated::new(k, v)
+    Field(Vec3<f32>),
-    }).collect()
+    Float(f32),
+    Int(u64),
+    Dim(Vec3u),
+}
+
+impl From<Vec3<f32>> for MeasurementValue {
+    fn from(v: Vec3<f32>) -> Self {
+        Self::Field(v)
+    }
+}
+impl From<f32> for MeasurementValue {
+    fn from(v: f32) -> Self {
+        Self::Float(v)
+    }
+}
+impl From<u64> for MeasurementValue {
+    fn from(v: u64) -> Self {
+        Self::Int(v)
+    }
+}
+impl From<Vec3u> for MeasurementValue {
+    fn from(v: Vec3u) -> Self {
+        Self::Dim(v)
+    }
+}
+
+#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub struct Measurement {
+    name: String,
+    value: MeasurementValue,
+    /// e.g. "A" for Amps
+    unit: String,
+}
+
+impl Measurement {
+    fn new<T: Into<MeasurementValue>>(name: &str, value: T, unit: &str) -> Self {
+        Self {
+            name: name.to_owned(),
+            value: value.into(),
+            unit: unit.to_owned(),
+        }
+    }
+    fn new_unitless<T: Into<MeasurementValue>>(name: &str, value: T) -> Self {
+        Self::new(name, value, "")
+    }
+
+    pub fn name(&self) -> &str {
+        &self.name
+    }
+
+    pub fn pretty_print(&self) -> String {
+        use MeasurementValue::*;
+        match self.value {
+            Field(v) => format!("{}{}", v, self.unit),
+            Float(f) => if self.unit != "" {
+                SiScale::format_short(f, &self.unit)
+            } else {
+                f.to_string()
+            },
+            Int(u)   => format!("{}{}", u, self.unit),
+            Dim(v)   => format!("{}x{}x{}{}", v.x(), v.y(), v.z(), self.unit),
+        }
+    }
+
+    /// format the Measurement in a way that could be parseable later.
+    /// one major use case for this is in dumping the type to a CSV.
+    pub fn machine_readable(&self) -> String {
+        use MeasurementValue::*;
+        match self.value {
+            Field(v) => format!("{},{},{}", v.x(), v.y(), v.z()),
+            Float(f) => f.to_string(),
+            Int(u)   => u.to_string(),
+            Dim(v)   => format!("{},{},{}", v.x(), v.y(), v.z()),
+        }
+    }
+
+    /// retrieve the float value of this measurement -- if it's of float type.
+    /// useful for tests
+    pub fn get_float(&self) -> Option<f32> {
+        match self.value {
+            MeasurementValue::Float(f) => Some(f),
+            _ => None,
+        }
+    }
+}
+
+impl<S> AbstractMeasurement<S> for Measurement {
+    fn key_value(&self, _state: &S) -> Vec<Measurement> {
+        vec![self.clone()]
+    }
 }
 
 enum SiScale {
@@ -46,7 +130,7 @@ enum SiScale {
 impl SiScale {
     fn for_value(v: f32) -> Self {
         use SiScale::*;
-        match v {
+        match v.abs() {
             v if v < 1e-12 => Unit,
             v if v < 1e-9  => Pico,
             v if v < 1e-6  => Nano,
@@ -99,7 +183,7 @@ impl SiScale {
     /// e.g. `format_short(1234, "A") -> "1.23 kA"
     fn format_short(v: f32, unit: &str) -> String {
         let si = SiScale::for_value(v);
-        let scaled = si.scale() * v;
+        let scaled = v / si.scale();
         format!("{:.2} {}{}", scaled, si.shortcode(), unit)
     }
 }
@@ -108,14 +192,11 @@ impl SiScale {
 pub struct Time;
 
 impl<S: AbstractSim> AbstractMeasurement<S> for Time {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
-        format!("{} (step {})", SiScale::format_short(state.time(), "s"), state.step_no())
+        vec![
-    }
+            Measurement::new_unitless("step", state.step_no()),
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
+            Measurement::new("time", state.time(), "s"),
-        [
+        ]
-            ("step".to_string(), state.step_no().to_string()),
-            ("time".to_string(), state.time().to_string()),
-        ].into_iter().collect()
     }
 }
 
@@ -123,49 +204,14 @@ impl<S: AbstractSim> AbstractMeasurement<S> for Time {
 pub struct Meta;
 
 impl<S: AbstractSim> AbstractMeasurement<S> for Meta {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
-        format!("{}x{}x{} feat: {:.1e}m", state.width(), state.height(), state.depth(), state.feature_size())
+        vec![
-    }
+            Measurement::new_unitless("dim", state.size().0),
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
+            Measurement::new("feature_size", state.feature_size(), "m"),
-        [
+        ]
-            ("width".to_string(), state.width().to_string()),
-            ("height".to_string(), state.height().to_string()),
-            ("depth".to_string(), state.depth().to_string()),
-            ("feature_size".to_string(), state.feature_size().to_string()),
-        ].into_iter().collect()
     }
 }
 
-/// some measurement which has already been evaluated.
-/// this is used particularly if we need to monomorphize a measurement (e.g. for serialization)
-/// and know it won't be applied to a new/different state.
-#[derive(Clone, Serialize, Deserialize)]
-pub struct Evaluated(String, String);
-
-impl Evaluated {
-    pub fn new<S1: Into<String>, S2: Into<String>>(key: S1, value: S2) -> Self {
-        Self(key.into(), value.into())
-    }
-    pub fn key(&self) -> &str {
-        &*self.0
-    }
-    pub fn value(&self) -> &str {
-        &*self.1
-    }
-}
-
-impl<S> AbstractMeasurement<S> for Evaluated {
-    fn eval(&self, _state: &S) -> String {
-        format!("{}: {}", self.0, self.1)
-    }
-    fn key_value(&self, _state: &S) -> IndexMap<String, String> {
-        [
-            (self.0.clone(), self.1.clone()),
-        ].into_iter().collect()
-    }
-}
-
-#[derive(Clone, Serialize, Deserialize)]
 pub struct Volume {
     name: String,
     region: Box<dyn Region>,
@@ -187,20 +233,13 @@ impl Volume {
 }
 
 impl<S: AbstractSim> AbstractMeasurement<S> for Volume {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
-        format!("Vol({}): {:.2e} um^3",
+        vec![
-            self.name,
+            Measurement::new(&format!("Vol({})", self.name), self.data(state), "um^3"),
-            self.data(state),
+        ]
-        )
-    }
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
-        [
-            (format!("Vol({})", self.name), self.data(state).to_string()),
-        ].into_iter().collect()
     }
 }
 
-#[derive(Clone, Serialize, Deserialize)]
 pub struct Current {
     name: String,
     region: Box<dyn Region>,
@@ -224,6 +263,24 @@ impl Current {
     }
 }
 
+// TODO: clean up these FieldSample types
+#[derive(Default)]
+struct TupleSum<T>(T);
+
+impl<T0: Default + AddAssign, T1: Default + AddAssign> std::iter::Sum for TupleSum<(T0, T1)> {
+    fn sum<I>(iter: I) -> Self
+        where I: Iterator<Item = Self>
+    {
+        let mut s = Self::default();
+        for TupleSum((a0, a1)) in iter {
+            s.0.0 += a0;
+            s.0.1 += a1;
+        }
+        s
+    }
+}
+
+
 #[derive(Default)]
 struct FieldSample(u32, f64, Vec3<f64>);
 
@@ -270,62 +327,78 @@ impl std::iter::Sum for FieldSamples<[FieldSample; 3]> {
 }
 
 impl<S: AbstractSim> AbstractMeasurement<S> for Current {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let (mean_current_mag, mean_current_vec) = self.data(state);
-        format!("I/cell({}): {:.2e} {:.2e}",
+        vec![
-            self.name,
+            Measurement::new(
-            mean_current_mag,
+                &format!("Imag/cell({})", self.name),
-            mean_current_vec)
+                mean_current_mag,
-    }
+                "A",
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
+            ),
-        let (mean_current_mag, mean_current_vec) = self.data(state);
+            Measurement::new(
-        [
+                &format!("I/cell({})", self.name),
-            (format!("Imag/cell({})", self.name), mean_current_mag.to_string()),
+                mean_current_vec,
-            (format!("I/cell({})", self.name), mean_current_vec.to_string()),
+                "A",
-        ].into_iter().collect()
+            ),
+        ]
     }
 }
 
 /// Measures the current directed around a closed loop
 #[derive(Clone, Serialize, Deserialize)]
-pub struct CurrentLoop {
+pub struct CurrentLoop<R> {
     name: String,
-    region: Torus
+    region: R,
 }
 
-impl CurrentLoop {
+impl<R> CurrentLoop<R> {
-    pub fn new(name: &str, r: Torus) -> Self {
+    pub fn new(name: &str, r: R) -> Self {
         Self {
             name: name.into(),
             region: r,
         }
     }
+}
+impl<R: Region + HasCrossSection> CurrentLoop<R> {
     fn data<S: AbstractSim>(&self, state: &S) -> f32 {
-        let FieldSample(volume, directed_current, _current_vec) = state.map_sum_over_enumerated(&self.region, |coord: Meters, _cell| {
+        // - current exists as a property of a 2d surface.
-            let normal = self.region.axis();
+        // - the user has provided us a 3d volume which behaves as though it's an extruded surface:
-            let to_coord = *coord - *self.region.center();
+        //   for any point in the volume we can query the normal vector of the cross section
-            let tangent = normal.cross(to_coord).norm();
+        //   containing that point.
-            let current = state.current(coord);
+        // - we choose that measuring the "current" on such a volume means to measure the average
-            let directed_current = current.dot(tangent.cast());
+        //   current through all its cross sections (for most boring materials, each
-            FieldSample(1, directed_current.cast(), current.cast())
+        //   cross section has nearly identical current).
+        // - therefore, enumerate the entire volume and compute the "net" current (the sum over
+        //   each cell of whatever current in that cell is along the cross-section normal).
+        //   then divide by the number of complete cross sections we measured, to average.
+        let feature_area = state.feature_size() * state.feature_size();
+        let TupleSum((net_current, cross_sections)) = state.map_sum_over_enumerated(&self.region, move |coord: Meters, _cell| {
+            // `normal` represents both the size of the cross section (m^2) this cell belongs to,
+            // and the normal direction of the cross section.
+            let normal = self.region.cross_section_normal(coord);  // [m^2]
+            // calculate the amount of normal current through this specific cell
+            let current_density = state.current_density(coord);    // [A/m^2]
+            let cross_sectional_current = feature_area * current_density.dot(normal.norm()); // [A]
+            // keep track of how many cross sections we enumerate, since each additional cross
+            // sections represents a double-count of the current.
+            let num_cross_sections_filled = feature_area / normal.mag();
+            TupleSum((cross_sectional_current, num_cross_sections_filled))
         });
-        let mean_directed_current = directed_current.cast::<f32>() / f32::from_primitive(volume);
+        let mean_cross_sectional_current = net_current.cast::<f32>() / cross_sections;
-        let cross_section = self.region.cross_section() / (state.feature_size() * state.feature_size());
+        mean_cross_sectional_current
-        let cross_sectional_current = mean_directed_current * cross_section;
-        cross_sectional_current
     }
 }
 
-impl<S: AbstractSim> AbstractMeasurement<S> for CurrentLoop {
+impl<R: Region + HasCrossSection, S: AbstractSim> AbstractMeasurement<S> for CurrentLoop<R> {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let cross_sectional_current = self.data(state);
-        format!("I({}): {:.2e}", self.name, cross_sectional_current)
+        vec![
-    }
+            Measurement::new(
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
+                &format!("I({})", self.name),
-        let cross_sectional_current = self.data(state);
+                cross_sectional_current,
-        [
+                "A"
-            (format!("I({})", self.name), cross_sectional_current.to_string()),
+            ),
-        ].into_iter().collect()
+        ]
     }
 }
 
@@ -385,27 +458,17 @@ impl MagneticLoop {
 }
 
 impl<S: AbstractSim> AbstractMeasurement<S> for MagneticLoop {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let (mean_directed_m, mean_directed_b, mean_directed_h) = self.data(state);
-        format!(
+        vec![
-            "M({}): {:.2e}; B({}): {:.2e}; H({}): {:.2e}",
+            Measurement::new_unitless(&format!("M({})", self.name), mean_directed_m),
-            self.name, mean_directed_m,
+            Measurement::new_unitless(&format!("B({})", self.name), mean_directed_b),
-            self.name, mean_directed_b,
+            Measurement::new_unitless(&format!("H({})", self.name), mean_directed_h),
-            self.name, mean_directed_h,
+        ]
-        )
-    }
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
-        let (mean_directed_m, mean_directed_b, mean_directed_h) = self.data(state);
-        [
-            (format!("M({})", self.name), mean_directed_m.to_string()),
-            (format!("B({})", self.name), mean_directed_b.to_string()),
-            (format!("H({})", self.name), mean_directed_h.to_string()),
-        ].into_iter().collect()
     }
 }
 
 /// mean M over a region
-#[derive(Clone, Serialize, Deserialize)]
 pub struct MagneticFlux {
     name: String,
     region: Box<dyn Region>,
@@ -430,20 +493,18 @@ impl MagneticFlux {
 }
 
 impl<S: AbstractSim> AbstractMeasurement<S> for MagneticFlux {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let mean_mag = self.data(state);
-        format!("Bavg({}): {:.2e}", self.name, mean_mag)
+        vec![
-    }
+            Measurement::new_unitless(
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
+                &format!("Bavg({})", self.name),
-        let mean_mag = self.data(state);
+                mean_mag,
-        [
+            )
-            (format!("Bavg({})", self.name), mean_mag.to_string()),
+        ]
-        ].into_iter().collect()
     }
 }
 
 /// mean B over a region
-#[derive(Clone, Serialize, Deserialize)]
 pub struct Magnetization {
     name: String,
     region: Box<dyn Region>,
@@ -468,15 +529,13 @@ impl Magnetization {
 }
 
 impl<S: AbstractSim> AbstractMeasurement<S> for Magnetization {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let mean_mag = self.data(state);
-        format!("Mavg({}): {:.2e}", self.name, mean_mag)
+        vec![
-    }
+            Measurement::new_unitless(
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
+                &format!("Mavg({})", self.name), mean_mag
-        let mean_mag = self.data(state);
+            ),
-        [
+        ]
-            (format!("Mavg({})", self.name), mean_mag.to_string()),
-        ].into_iter().collect()
     }
 }
 
@@ -489,15 +548,11 @@ fn loc(v: Meters) -> String {
 pub struct MagnetizationAt(pub Meters);
 
 impl<S: AbstractSim> AbstractMeasurement<S> for MagnetizationAt {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let m = state.sample(self.0).m();
-        format!("M{}: {:.2e}", loc(self.0), m)
+        vec![
-    }
+            Measurement::new_unitless(&format!("M{}", loc(self.0)), m.cast())
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
+        ]
-        let m = state.sample(self.0).m();
-        [
-            (format!("M{}", loc(self.0)), m.to_string()),
-        ].into_iter().collect()
     }
 }
 
@@ -506,15 +561,13 @@ impl<S: AbstractSim> AbstractMeasurement<S> for MagnetizationAt {
 pub struct MagneticFluxAt(pub Meters);
 
 impl<S: AbstractSim> AbstractMeasurement<S> for MagneticFluxAt {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let b = state.sample(self.0).b();
-        format!("B{}: {:.2e}", loc(self.0), b)
+        vec![
-    }
+            Measurement::new_unitless(
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
+                &format!("B{}", loc(self.0)), b.cast()
-        let b = state.sample(self.0).b();
+            )
-        [
+        ]
-            (format!("B{}", loc(self.0)), b.to_string()),
-        ].into_iter().collect()
     }
 }
 
@@ -523,15 +576,13 @@ impl<S: AbstractSim> AbstractMeasurement<S> for MagneticFluxAt {
 pub struct MagneticStrengthAt(pub Meters);
 
 impl<S: AbstractSim> AbstractMeasurement<S> for MagneticStrengthAt {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let h = state.sample(self.0).h();
-        format!("H{}: {:.2e}", loc(self.0), h)
+        vec![
-    }
+            Measurement::new_unitless(
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
+                &format!("H{}", loc(self.0)), h.cast()
-        let h = state.sample(self.0).h();
+            )
-        [
+        ]
-            (format!("H{}", loc(self.0)), h.to_string()),
-        ].into_iter().collect()
     }
 }
 
@@ -539,19 +590,16 @@ impl<S: AbstractSim> AbstractMeasurement<S> for MagneticStrengthAt {
 pub struct ElectricField(pub Meters);
 
 impl<S: AbstractSim> AbstractMeasurement<S> for ElectricField {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let e = state.sample(self.0).e();
-        format!("E{}: {}", loc(self.0), e)
+        vec![
-    }
+            Measurement::new_unitless(
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
+                &format!("E{}", loc(self.0)), e.cast()
-        let e = state.sample(self.0).e();
+            )
-        [
+        ]
-            (format!("E{}", loc(self.0)), e.to_string()),
-        ].into_iter().collect()
     }
 }
 
-#[derive(Clone, Serialize, Deserialize)]
 pub struct Energy {
     name: String,
     region: Box<dyn Region>,
@@ -585,19 +633,16 @@ impl Energy {
 }
 
 impl<S: AbstractSim> AbstractMeasurement<S> for Energy {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let e = self.data(state);
-        format!("U({}): {}", self.name, SiScale::format_short(e, "J"))
+        vec![
-    }
+            Measurement::new(
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
+                &format!("U({})", self.name), e, "J"
-        let e = self.data(state);
+            )
-        [
+        ]
-            (format!("U({})", self.name), e.to_string()),
-        ].into_iter().collect()
     }
 }
 
-#[derive(Clone, Serialize, Deserialize)]
 pub struct Power {
     name: String,
     region: Box<dyn Region>
@@ -626,14 +671,188 @@ impl Power {
 }
 
 impl<S: AbstractSim> AbstractMeasurement<S> for Power {
-    fn eval(&self, state: &S) -> String {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let power = self.data(state);
-        format!("P({}): {}", self.name, SiScale::format_short(power, "W"))
+        vec![
-    }
+            Measurement::new(
-    fn key_value(&self, state: &S) -> IndexMap<String, String> {
+                &format!("P({})", self.name), power, "W"
-        let power = self.data(state);
+            )
-        [
+        ]
-            (format!("P({})", self.name), power.to_string()),
+    }
-        ].into_iter().collect()
+}
+
+#[cfg(test)]
+pub mod test {
+    use super::*;
+    use crate::cross::mat::AnisomorphicConductor;
+    use crate::cross::step::SimMeta;
+    use crate::geom::Index;
+    use crate::sim::{Fields, GenericSim, StaticSim};
+    use crate::stim::AbstractStimulus;
+
+    struct MockSim {
+        e_field: Vec3<f32>,
+        dim: Vec3u,
+        feature_size: f32,
+        mat: AnisomorphicConductor<f32>,
+    }
+    impl AbstractSim for MockSim {
+        type Real = f32;
+        type Material = AnisomorphicConductor<f32>;
+
+        fn meta(&self) -> SimMeta<f32> {
+            SimMeta::new(self.dim, self.feature_size, 1e-9)
+        }
+        fn step_no(&self) -> u64 {
+            unimplemented!()
+        }
+        fn fields_at_index(&self, _pos: Index) -> Fields<Self::Real> {
+            Fields::new(self.e_field, Vec3::zero(), Vec3::zero())
+        }
+        fn get_material_index(&self, _at: Index) -> &Self::Material {
+            &self.mat
+        }
+        fn put_material_index(&mut self, _at: Index, _m: Self::Material) {
+            unimplemented!()
+        }
+        fn step_multiple<S: AbstractStimulus>(&mut self, _num_steps: u32, _s: &S) {
+            unimplemented!()
+        }
+        fn to_static(&self) -> StaticSim {
+            unimplemented!()
+        }
+        fn to_generic(&self) -> GenericSim<Self::Real> {
+            unimplemented!()
+        }
+    }
+
+    struct MockRegion {
+        normal: Vec3<f32>,
+    }
+    impl HasCrossSection for MockRegion {
+        fn cross_section_normal(&self, _p: Meters) -> Vec3<f32> {
+            self.normal
+        }
+    }
+    impl Region for MockRegion {
+        fn contains(&self, _p: Meters) -> bool {
+            true
+        }
+    }
+
+    #[test]
+    fn current_loop_trivial() {
+        let sim = MockSim {
+            e_field: Vec3::new(1.0, 0.0, 0.0),
+            dim: Vec3u::new(1, 1, 1),
+            feature_size: 1.0,
+            mat: AnisomorphicConductor::new(Vec3::new(1.0, 1.0, 1.0)),
+        };
+        let region = MockRegion {
+            normal: Vec3::new(1.0, 0.0, 0.0),
+        };
+
+        let kv = CurrentLoop::new("test", region).key_value(&sim);
+        assert_eq!(kv.len(), 1);
+        // measured area is 1 m^2
+        // region cross-section is 1 m^2
+        // conductivity is 1 S/m
+        assert_eq!(kv[0].get_float().unwrap(), 1.0);
+    }
+    #[test]
+    fn current_loop_multi_cell() {
+        let sim = MockSim {
+            e_field: Vec3::new(1.0, 0.0, 0.0),
+            dim: Vec3u::new(4, 4, 4),
+            feature_size: 0.25,
+            mat: AnisomorphicConductor::new(Vec3::new(1.0, 1.0, 1.0)),
+        };
+        let region = MockRegion {
+            normal: Vec3::new(1.0, 0.0, 0.0),
+        };
+
+        let kv = CurrentLoop::new("test", region).key_value(&sim);
+        assert_eq!(kv.len(), 1);
+        // measured area is 1 m^2
+        // region cross-section is 1 m^2
+        // conductivity is 1 S/m
+        assert_eq!(kv[0].get_float().unwrap(), 1.0);
+    }
+    #[test]
+    fn current_loop_off_conductor() {
+        let sim = MockSim {
+            e_field: Vec3::new(1.0, 1.0, 1.0),
+            dim: Vec3u::new(4, 4, 4),
+            feature_size: 0.25,
+            mat: AnisomorphicConductor::new(Vec3::new(0.0, 1.0, 1.0)),
+        };
+        let region = MockRegion {
+            normal: Vec3::new(1.0, 0.0, 0.0),
+        };
+
+        let kv = CurrentLoop::new("test", region).key_value(&sim);
+        assert_eq!(kv.len(), 1);
+        // material is not conductive in the direction being queried
+        assert_eq!(kv[0].get_float().unwrap(), 0.0);
+    }
+    #[test]
+    fn current_loop_e_field() {
+        let sim = MockSim {
+            e_field: Vec3::new(4.0, 2.0, 1.0),
+            dim: Vec3u::new(4, 4, 4),
+            feature_size: 0.25,
+            mat: AnisomorphicConductor::new(Vec3::new(1.0, 1.0, 1.0)),
+        };
+        let region = MockRegion {
+            normal: Vec3::new(1.0, 0.0, 0.0),
+        };
+
+        let kv = CurrentLoop::new("test", region).key_value(&sim);
+        assert_eq!(kv.len(), 1);
+        // measured area is 1 m^2
+        // region cross-section is 1 m^2
+        // conductivity is 1 S/m
+        // e field is 4 V/m
+        assert_eq!(kv[0].get_float().unwrap(), 4.0);
+    }
+    #[test]
+    fn current_loop_conductivity() {
+        let sim = MockSim {
+            e_field: Vec3::new(4.0, 2.0, 1.0),
+            dim: Vec3u::new(4, 4, 4),
+            feature_size: 0.25,
+            mat: AnisomorphicConductor::new(Vec3::new(3.0, 1.0, 1.0)),
+        };
+        let region = MockRegion {
+            normal: Vec3::new(1.0, 0.0, 0.0),
+        };
+
+        let kv = CurrentLoop::new("test", region).key_value(&sim);
+        assert_eq!(kv.len(), 1);
+        // measured area is 1 m^2
+        // region cross-section is 1 m^2
+        // conductivity is 3 S/m
+        // e field is 4 V/m
+        assert_eq!(kv[0].get_float().unwrap(), 3.0*4.0);
+    }
+    #[test]
+    fn current_loop_cross_section() {
+        let sim = MockSim {
+            e_field: Vec3::new(4.0, 2.0, 1.0),
+            dim: Vec3u::new(4, 4, 4),
+            feature_size: 0.5,
+            mat: AnisomorphicConductor::new(Vec3::new(3.0, 1.0, 1.0)),
+        };
+        let region = MockRegion {
+            normal: Vec3::new(16.0, 0.0, 0.0),
+        };
+
+        let kv = CurrentLoop::new("test", region).key_value(&sim);
+        assert_eq!(kv.len(), 1);
+        // measured area is 2 m^2
+        // region cross-section is 16 m^2
+        // conductivity is 3 S/m
+        // e field is 4 V/m
+        assert_eq!(kv[0].get_float().unwrap(), 3.0*4.0*16.0);
     }
 }

crates/coremem/src/render.rs

@@ -2,7 +2,7 @@ use crate::geom::Index;
 use crate::real::ToFloat as _;
 use crate::cross::vec::{Vec2, Vec3};
 use crate::sim::{AbstractSim, GenericSim, Sample};
-use crate::meas::{self, AbstractMeasurement};
+use crate::meas::{self, AbstractMeasurement, Measurement};
 use crossterm::{cursor, QueueableCommand as _};
 use crossterm::style::{style, Color, PrintStyledContent, Stylize as _};
 use font8x8::{BASIC_FONTS, GREEK_FONTS, UnicodeFonts as _};
@@ -294,8 +294,8 @@ impl<'a, S: AbstractSim> RenderSteps<'a, S> {
         self.im = im;
     }
     fn render_measurements(&mut self) {
-        for (meas_no, m) in self.meas.iter().enumerate() {
+        for (meas_no, m) in meas::eval_multiple(self.sim, &self.meas).into_iter().enumerate() {
-            let meas_string = m.eval(self.sim);
+            let meas_string = m.pretty_print();
             for (i, c) in meas_string.chars().enumerate() {
                 let glyph = BASIC_FONTS.get(c)
                     .or_else(|| GREEK_FONTS.get(c))
@@ -425,6 +425,7 @@ impl<S: AbstractSim> Renderer<S> for ColorTermRenderer {
         measurements: &[&dyn AbstractMeasurement<S>],
         config: RenderConfig,
     ) {
+        let measurements = meas::eval_multiple(state, measurements);
         let (max_w, mut max_h) = crossterm::terminal::size().unwrap();
         max_h = max_h.saturating_sub(2 + measurements.len() as u16);
         let im = RenderSteps::render_configured(state, &[], z, (max_w as _, max_h as _), config);
@@ -453,7 +454,7 @@ impl<S: AbstractSim> Renderer<S> for ColorTermRenderer {
         for m in measurements {
             // Measurements can be slow to compute
             stdout.flush().unwrap();
-            let meas_string = m.eval(state);
+            let meas_string = m.pretty_print();
             stdout.queue(cursor::MoveDown(1)).unwrap();
             stdout.queue(cursor::MoveToColumn(1)).unwrap();
             stdout.queue(PrintStyledContent(style(meas_string))).unwrap();
@@ -592,7 +593,7 @@ pub struct SerializedFrame<S> {
     pub state: S,
     /// although not generally necessary to load the sim, saving the measurements is beneficial for
     /// post-processing.
-    pub measurements: Vec<meas::Evaluated>,
+    pub measurements: Vec<Measurement>,
 }
 
 impl<S: AbstractSim> SerializedFrame<S> {
@@ -633,15 +634,18 @@ impl SerializerRenderer {
 }
 
 impl SerializerRenderer {
-    fn serialize<S: AbstractSim + Serialize>(&self, state: &S, measurements: Vec<meas::Evaluated>) {
+    fn serialize<S: AbstractSim + Serialize>(&self, state: &S, measurements: Vec<Measurement>) {
         let frame = SerializedFrame {
             state,
             measurements,
         };
         let name = self.fmt_str.replace("{step_no}", &*frame.state.step_no().to_string());
-        let out = BufWriter::new(File::create(name).unwrap());
+        // serialize to a temporary file -- in case we run out of disk space, etc.
-        //serde_cbor::to_writer(out, &snap).unwrap();
+        let temp_name = format!("{}.incomplete", name);
+        let out = BufWriter::new(File::create(&temp_name).unwrap());
         bincode::serialize_into(out, &frame).unwrap();
+        // atomically complete the write.
+        std::fs::rename(temp_name, name).unwrap();
     }
 
     pub fn try_load<S: AbstractSim + for <'a> Deserialize<'a>>(&self) -> Option<SerializedFrame<S>> {
@@ -656,9 +660,9 @@ impl<S: AbstractSim + Serialize> Renderer<S> for SerializerRenderer {
     }
     fn render(&self, state: &S, measurements: &[&dyn AbstractMeasurement<S>], _config: RenderConfig) {
         if self.prefer_generic {
-            self.serialize(&state.to_generic(), meas::eval_to_vec(state, measurements));
+            self.serialize(&state.to_generic(), meas::eval_multiple(state, measurements));
         } else {
-            self.serialize(state, meas::eval_to_vec(state, measurements));
+            self.serialize(state, meas::eval_multiple(state, measurements));
         }
     }
 }
@@ -698,7 +702,7 @@ impl<S: AbstractSim> Renderer<S> for CsvRenderer {
         default_render_z_slice(self, state, z, measurements, config)
     }
     fn render(&self, state: &S, measurements: &[&dyn AbstractMeasurement<S>], _config: RenderConfig) {
-        let row = meas::eval_multiple_kv(state, measurements);
+        let row = meas::eval_multiple(state, measurements);
         let step = state.step_no();
         let mut lock = self.state.lock().unwrap();
         let mut writer = match lock.take().unwrap() {
@@ -732,13 +736,13 @@ impl<S: AbstractSim> Renderer<S> for CsvRenderer {
                     file.set_len(0).unwrap();
                     let mut writer = csv::Writer::from_writer(BufWriter::new(file));
                     // write the header
-                    writer.write_record(row.keys()).unwrap();
+                    writer.write_record(row.iter().map(|m| m.name())).unwrap();
                     writer
                 }
             },
             CsvState::Writing(writer) => writer,
         };
-        writer.write_record(row.values()).unwrap();
+        writer.write_record(row.iter().map(|m| m.machine_readable())).unwrap();
         writer.flush().unwrap();
         *lock = Some(CsvState::Writing(writer));
     }

crates/coremem/src/sim/mod.rs

@@ -1,3 +1,4 @@
+use crate::diagnostics::SyncDiagnostics;
 use crate::geom::{Coord, Cube, Index, InvertedRegion, Region};
 use crate::cross::mat::{FullyGenericMaterial, Material, Vacuum};
 use crate::cross::real::Real;
@@ -209,6 +210,9 @@ pub trait AbstractSim: Sync {
     /// Take a "snapshot" of the simulation, dropping all material-specific information.
     fn to_static(&self) -> StaticSim;
     fn to_generic(&self) -> GenericSim<Self::Real>;
+    fn use_diagnostics(&mut self, _diag: SyncDiagnostics) {
+        // optional
+    }
 
 
     //--- HELPER METHODS below (derived) ---//
@@ -324,8 +328,15 @@ pub trait AbstractSim: Sync {
         }))).flatten().flatten().sum()
     }
 
+    /// returns the directed current at `c`, in `A / m^2`
+    fn current_density<C: Coord>(&self, c: C) -> Vec3<f32> {
+        self.sample(c).current_density().cast::<f32>()
+    }
+
+    /// returns the directed current at `c` in absolute units, `A`, or rather, `A` per cell, since
+    /// this looks at just a single cell. you probably want to use `current_density`.
     fn current<C: Coord>(&self, c: C) -> Vec3<f32> {
-        self.sample(c).current_density().cast::<f32>() * self.feature_size() * self.feature_size()
+        self.current_density(c) * self.feature_size() * self.feature_size()
     }
 
     fn fill_region_using<C, Reg, F, M>(&mut self, region: &Reg, f: F)

crates/coremem/src/sim/spirv/mod.rs

@@ -2,6 +2,7 @@ use ndarray::{self, Array3};
 use serde::{Deserialize, Serialize};
 use log::{info, trace, warn};
 
+use crate::diagnostics::SyncDiagnostics;
 use crate::geom::{Coord, Index};
 use crate::real::Real;
 use crate::sim::{AbstractSim, Fields, GenericSim, StaticSim};
@@ -49,6 +50,8 @@ where M: 'static
     step_no: u64,
     #[serde(skip)]
     backend: B,
+    #[serde(skip)]
+    diag: SyncDiagnostics,
 }
 
 // B isn't always clonable (e.g. gpu backend) so rust can't auto-derive this.
@@ -65,6 +68,7 @@ impl<R: Clone, M: Clone, B: Default> Clone for SpirvSim<R, M, B> {
             // we require that the caller explicitly init the backend if they need this.
             // TODO: this probably shouldn't be a `Clone` method, but like "read_only_clone()".
             backend: Default::default(),
+            diag: self.diag.clone(),
         }
     }
 }
@@ -131,6 +135,7 @@ where
             mat,
             step_no: self.step_no,
             backend: Default::default(),
+            diag: self.diag.clone(),
         }
     }
 
@@ -144,8 +149,13 @@ where
             mat,
             step_no: self.step_no,
             backend: Default::default(),
+            diag: self.diag.clone(),
         }
     }
+
+    fn use_diagnostics(&mut self, diag: SyncDiagnostics) {
+        self.diag = diag;
+    }
 }
 
 impl<R: Real, M: Default, B> SpirvSim<R, M, B>
@@ -174,6 +184,7 @@ impl<R: Real, M: Default, B> SpirvSim<R, M, B>
             mat,
             step_no: 0,
             backend,
+            diag: SyncDiagnostics::new(),
         }
     }
 }
@@ -217,27 +228,30 @@ where
         -> (Array3<Vec3<R>>, Array3<Vec3<R>>)
     {
         trace!("eval_stimulus begin");
-        let dim = self.size();
+        let (e, h) = self.diag.instrument_stimuli(|| {
-        let feature_size = self.feature_size();
+            let dim = self.size();
-        let t_sec = self.time();
+            let feature_size = self.feature_size();
-        let timestep = self.meta.time_step;
+            let t_sec = self.time();
+            let timestep = self.meta.time_step;
 
-        // TODO(perf): do this in one loop!
+            // TODO(perf): do this in one loop!
-        let e = ndarray::Zip::from(ndarray::indices(
+            let e = ndarray::Zip::from(ndarray::indices(
-            [dim.z() as usize, dim.y() as usize, dim.x() as usize]
+                [dim.z() as usize, dim.y() as usize, dim.x() as usize]
-        )).par_map_collect(|(z, y, x)| {
+            )).par_map_collect(|(z, y, x)| {
-            let pos_idx = Index::new(x as _, y as _, z as _);
+                let pos_idx = Index::new(x as _, y as _, z as _);
-            let pos_meters = pos_idx.to_meters(feature_size);
+                let pos_meters = pos_idx.to_meters(feature_size);
-            let densities = stim.at(t_sec, pos_meters);
+                let densities = stim.at(t_sec, pos_meters);
-            densities.e.cast::<R>() * timestep
+                densities.e.cast::<R>() * timestep
-        });
+            });
-        let h = ndarray::Zip::from(ndarray::indices(
+            let h = ndarray::Zip::from(ndarray::indices(
-            [dim.z() as usize, dim.y() as usize, dim.x() as usize]
+                [dim.z() as usize, dim.y() as usize, dim.x() as usize]
-        )).par_map_collect(|(z, y, x)| {
+            )).par_map_collect(|(z, y, x)| {
-            let pos_idx = Index::new(x as _, y as _, z as _);
+                let pos_idx = Index::new(x as _, y as _, z as _);
-            let pos_meters = pos_idx.to_meters(feature_size);
+                let pos_meters = pos_idx.to_meters(feature_size);
-            let densities = stim.at(t_sec, pos_meters);
+                let densities = stim.at(t_sec, pos_meters);
-            densities.h.cast::<R>() * timestep
+                densities.h.cast::<R>() * timestep
+            });
+            (e, h)
         });
         trace!("eval_stimulus end");
         (e, h)

crates/coremem/src/stim.rs

@@ -1,6 +1,6 @@
-use crate::real::*;
+use crate::real::Real as _;
 use crate::cross::vec::Vec3;
-use crate::geom::{Meters, Region};
+use crate::geom::{HasCrossSection, Meters, Region};
 use rand;
 
 /// field densities
@@ -150,23 +150,19 @@ impl<R: Region + Sync, T: TimeVarying3 + Sync> AbstractStimulus for Stimulus<R,
 pub struct CurlStimulus<R, T> {
     region: R,
     stim: T,
-    center: Meters,
-    axis: Meters,
 }
 
 impl<R, T> CurlStimulus<R, T> {
-    pub fn new(region: R, stim: T, center: Meters, axis: Meters) -> Self {
+    pub fn new(region: R, stim: T) -> Self {
-        Self { region, stim, center, axis }
+        Self { region, stim }
     }
 }
 
-impl<R: Region + Sync, T: TimeVarying1 + Sync> AbstractStimulus for CurlStimulus<R, T> {
+impl<R: Region + HasCrossSection, T: TimeVarying1 + Sync> AbstractStimulus for CurlStimulus<R, T> {
     fn at(&self, t_sec: f32, pos: Meters) -> Fields {
         if self.region.contains(pos) {
             let FieldMags { e, h } = self.stim.at(t_sec);
-            let from_center_to_point = *pos - *self.center;
+            let rotational = self.region.cross_section_normal(pos).norm();
-            // TODO: is this inverted?
-            let rotational = from_center_to_point.cross(*self.axis);
             let impulse_e = rotational.with_mag(e.cast()).unwrap_or_default();
             let impulse_h = rotational.with_mag(h.cast()).unwrap_or_default();
             Fields {
@@ -398,4 +394,60 @@ mod test {
         assert_approx_eq!(s.at(0.00050), Vec3::zero(), Vec3::zero());
         assert_approx_eq!(s.at(0.00075), Vec3::new(-10.0, -1.0, 100.0), Vec3::zero());
     }
+
+    struct MockRegion {
+        normal: Vec3<f32>,
+    }
+    impl HasCrossSection for MockRegion {
+        fn cross_section_normal(&self, _p: Meters) -> Vec3<f32> {
+            self.normal
+        }
+    }
+    impl Region for MockRegion {
+        fn contains(&self, _p: Meters) -> bool {
+            true
+        }
+    }
+
+    #[test]
+    fn curl_stimulus_trivial() {
+        let region = MockRegion {
+            normal: Vec3::new(1.0, 0.0, 0.0)
+        };
+        // stim acts as e: 1.0, h: 0.0
+        let stim = CurlStimulus::new(region, 1.0);
+        assert_eq!(stim.at(0.0, Meters::new(0.0, 0.0, 0.0)), Fields {
+            e: Vec3::new(1.0, 0.0, 0.0),
+            h: Vec3::new(0.0, 0.0, 0.0),
+        });
+    }
+
+    #[test]
+    fn curl_stimulus_scales_right() {
+        let region = MockRegion {
+            // the magnitude of this normal shouldn't influence the curl.
+            // though, maybe it would be more useful if it *did*?
+            normal: Vec3::new(5.0, 0.0, 0.0)
+        };
+        // stim acts as e: -3.0, h: 0.0
+        let stim = CurlStimulus::new(region, -3.0);
+        assert_eq!(stim.at(0.0, Meters::new(0.0, 0.0, 0.0)), Fields {
+            e: Vec3::new(-3.0, 0.0, 0.0),
+            h: Vec3::new(0.0, 0.0, 0.0),
+        });
+    }
+
+    #[test]
+    fn curl_stimulus_multi_axis() {
+        let region = MockRegion {
+            normal: Vec3::new(0.0, -1.0, 1.0)
+        };
+        // stim acts as e: 2.0, h: 0.0
+        let stim = CurlStimulus::new(region, 2.0);
+        let Fields { e, h: _ } = stim.at(0.0, Meters::new(0.0, 0.0, 0.0));
+        assert!(e.distance(
+                Vec3::new(0.0, -1.0, 1.0).with_mag(2.0).unwrap()
+            ) < 1e-6
+        );
+    }
 }

crates/cross/src/real.rs

@@ -107,6 +107,10 @@ pub trait Real:
         self == Self::zero()
     }
 
+    fn inv(self) -> Self {
+        Self::one() / self
+    }
+
     fn zero() -> Self;
     fn one() -> Self;
     fn two() -> Self;

crates/cross/src/step.rs

@@ -20,6 +20,17 @@ pub struct SimMeta<R> {
     pub feature_size: R,
 }
 
+impl<R: Real> SimMeta<R> {
+    pub fn new(dim: Vec3u, feature_size: R, time_step: R) -> Self {
+        Self {
+            dim,
+            inv_feature_size: feature_size.inv(),
+            time_step,
+            feature_size
+        }
+    }
+}
+
 impl<R: Copy> SimMeta<R> {
     pub fn dim(&self) -> Vec3u {
         self.dim

crates/cross/src/vec.rs

@@ -333,6 +333,8 @@ impl<R: Real> Vec3<R> {
         Self::new(self.x().exp(), self.y().exp(), self.z().exp())
     }
 
+    /// the only condition upon which this returns `None` is if the current magnitude is zero
+    /// and the new magnitude and NON-zero.
     pub fn with_mag(&self, new_mag: R) -> Option<Self> {
         if new_mag.is_zero() {
             // avoid div-by-zero if self.mag() == 0 and new_mag == 0

crates/post/src/bin/csv.rs

@@ -1,4 +1,5 @@
-// use crate::{Loader, LoaderCache};
+//! extracts Measurements from rendered .bc files and dumps them into a CSV
+
 use coremem_post::{Loader, LoaderCache};
 use std::path::PathBuf;
 use structopt::StructOpt;
@@ -15,13 +16,13 @@ fn main() {
     
     let mut frame = cache.load_first();
     for meas in frame.measurements() {
-        print!("\"{}\",", meas.key());
+        print!("\"{}\",", meas.name());
     }
     println!("");
 
     loop {
         for meas in frame.measurements() {
-            print!("\"{}\",", meas.value());
+            print!("\"{}\",", meas.machine_readable().replace(",", "\\,"));
         }
         println!("");
 

crates/post/src/bin/decimate.rs

@@ -1,4 +1,5 @@
-// use crate::Loader;
+//! "decimates" the binary rendered form of a simulation
+//! by deleting all but every N rendered files
 use coremem_post::Loader;
 use std::fs;
 use std::path::PathBuf;

crates/post/src/bin/viewer.rs

@@ -1,3 +1,7 @@
+//! interactive CLI viewer for .bc files.
+//! navigate through the simulation over time and space,
+//! and toggle views to see more detail over material, electric, or magnetic changes.
+
 use coremem_post::{Loader, Viewer};
 use std::path::PathBuf;
 use std::io::Write as _;

crates/post/src/lib.rs

@@ -36,7 +36,7 @@ pub struct Frame {
 }
 
 impl Frame {
-    pub fn measurements(&self) -> &[meas::Evaluated] {
+    pub fn measurements(&self) -> &[meas::Measurement] {
         &*self.data.measurements
     }
     pub fn sim(&self) -> &GenericSim<f32> {