<commit_before>#![feature (i128_type)]

extern crate serde;
#[macro_use]
extern crate serde_derive;

#[macro_use]
extern crate time_steward;

// Imports for the UI
#[macro_use]
extern crate glium;
#[path = "../dev-shared/emscripten_compatibility.rs"] mod emscripten_compatibility;
pub use emscripten_compatibility::canvas_click;

macro_rules! printlnerr(
    ($($arg:tt)*) => { {use std::io::Write;
        let r = writeln!(&mut ::std::io::stderr(), $($arg)*);
        r.expect("failed printing to stderr");
    } }
);


use std::cmp::{min, max};
use time_steward::{DeterministicRandomId};
use time_steward::rowless::api::{self, PersistentTypeId, PersistentlyIdentifiedType, ListedType, StewardData, QueryOffset, DataHandleTrait, DataTimelineCellTrait, ExtendedTime, Basics as BasicsTrait};
use time_steward::rowless::stewards::{simple_full as steward_module};
use steward_module::{TimeSteward, ConstructibleTimeSteward, IncrementalTimeSteward, Event, DataHandle, DataTimelineCell, EventHandle, Accessor, EventAccessor, FutureCleanupAccessor, SnapshotAccessor, simple_timeline};
use simple_timeline::{SimpleTimeline, GetVarying, IterateUniquelyOwnedPredictions, tracking_query, modify_simple_timeline, unmodify_simple_timeline};


type Time = i64;
type Distance = i64;
type Amount = i128;
const SECOND: Time = (1 as Time) << 20;
const METER: Distance = (1 as Distance) << 20;
const GENERIC_DENSITY: Amount = (SECOND as Amount) << 20;
const GENERIC_INK_AMOUNT: Amount = (METER as Amount)*(METER as Amount)*GENERIC_DENSITY;

// Velocity is "ink units per space unit per time unit"
// and slope is "ink density per space unit", i.e. "ink units per space unit^3"
// 
// Let's say we want flow across a meter^2 cell, from a density of GENERIC_INK_AMOUNT/m^2 to 0, to drain GENERIC_INK_AMOUNT per second. Then
//   slope = GENERIC_INK_AMOUNT/m^3
//   velocity = GENERIC_INK_AMOUNT/m*s
// so
const VELOCITY_PER_SLOPE: Amount = (METER as Amount)*(METER as Amount) / (SECOND as Amount);

type Steward = steward_module::Steward <Basics>;

macro_rules! serialization_cheat {
  (
    [$($bounds:tt)*]
    [$($concrete:tt)*]
  ) => {

/*thread_local! {
  static SERIALIZATION_CONTEXT: RefCell<Option <SerializationContext>> = RefCell::new (None);
  static DESERIALIZATION_CONTEXT: RefCell<Option <DeserializationContext>> = RefCell::new (None);
}*/


impl <$($bounds)*> $crate::serde::Serialize for $($concrete)* {
  fn serialize <S: $crate::serde::Serializer> (&self, _: S)->Result <S::Ok, S::Error> {
    unimplemented!()
  }
}

impl <'a, $($bounds)*> $crate::serde::Deserialize <'a> for $($concrete)* {
  fn deserialize <D: $crate::serde::Deserializer<'a>> (_: D)->Result <Self, D::Error> {
    unimplemented!()
  }
}

  };
}

#[derive (Clone, PartialEq, Eq, Serialize, Deserialize, Debug)]
struct NodeData {
  width: Distance,
  center: [Distance ; 2],
  parent: Option <NodeHandle>,
  varying: DataTimelineCell <SimpleTimeline <NodeVarying, Steward>>,
}
#[derive (Clone, PartialEq, Eq, Debug)]
struct NodeVarying {
  last_change: Time,
  ink_at_last_change: Amount,
  children: Vec<NodeHandle>,
  boundaries: [[Vec<BoundaryHandle>;2];2],
  
  accumulation_rate: Amount,
  slope: [Amount; 2],
}
type NodeHandle = DataHandle <NodeData>;
impl StewardData for NodeData {}
impl StewardData for NodeVarying {}
impl PersistentlyIdentifiedType for NodeData {
  const ID: PersistentTypeId = PersistentTypeId(0x734528f6aefdc1b9);
}
impl IterateUniquelyOwnedPredictions <Steward> for NodeVarying {}
serialization_cheat!([][NodeVarying]);

#[derive (Clone, PartialEq, Eq, Serialize, Deserialize, Debug)]
struct BoundaryData {
  length: Distance,
  center: [Distance ; 2],
  nodes: [NodeHandle; 2],
  varying: DataTimelineCell <SimpleTimeline <BoundaryVarying, Steward>>,
}
#[derive (Clone, PartialEq, Eq, Debug)]
struct BoundaryVarying {
  transfer_velocity: Amount,
  next_change: Option <EventHandle <Basics>>,
}
type BoundaryHandle = DataHandle <BoundaryData>;
impl StewardData for BoundaryData {}
impl StewardData for BoundaryVarying {}
impl PersistentlyIdentifiedType for BoundaryData {
  const ID: PersistentTypeId = PersistentTypeId(0x9d643214b58b24dc);
}
impl IterateUniquelyOwnedPredictions <Steward> for BoundaryVarying {
  fn iterate_predictions <F: FnMut (& <Steward as TimeSteward>::EventHandle)> (&self, callback: &mut F) {
    if let Some (prediction) = self.next_change.as_ref() {
      callback (prediction);
    }
  }
}
serialization_cheat!([][BoundaryVarying]);

macro_rules! get {
  ($accessor: expr, $cell: expr) => {
    $accessor.query ($cell, &GetVarying, QueryOffset::After).unwrap().1
  }
}
macro_rules! set {
  ($accessor: expr, $cell: expr, $field: ident, $value: expr) => {
    {
      let mut value = $accessor.query ($cell, &GetVarying, QueryOffset::After).unwrap().1;
      value.$field = $value;
      modify_simple_timeline ($accessor, $cell, Some (value));
    }
  }
}
macro_rules! set_with {
  ($accessor: expr, $cell: expr, | $varying: ident | $actions: expr) => {
    {
      let mut $varying = $accessor.query ($cell, &GetVarying, QueryOffset::After).unwrap().1;
      $actions;
      modify_simple_timeline ($accessor, $cell, Some ($varying));
    }
  }
}

/*struct EventContext <'a, A: EventAccessor <Steward = Steward>> {
  accessor: &'a A,
  node_clones: HashMap <DataHandle <NodeData>, Rc<RefCell<Node>>,
  boundary_clones: HashMap <DataHandle <NodeData>, Rc<RefCell<Node>>,
}
impl<A: EventAccessor <Steward = Steward>> EventContext <A> {
  fn new (accessor: & A)->Self {
    EventContext {
      accessor: accessor,
      node_clones: HashMap::new(),
      boundary_clones: HashMap::new(),
    }
  }
  fn get_node (&self, handle: & DataHandle <NodeData>)->Rc {
    self.node_clones.insert_or (self.accessor.query (
  }
}
impl<A: EventAccessor <Steward = Steward>> Drop for EventContext <A> {
  fn drop (&mut self) {
    for (handle, node) in self.node_clones.iter() {
      if node.modified {
        modify_simple_timeline (self.accessor, handle, Some(node.varying));
      }
    }
    for (handle, boundary) in self.node_clones.iter() {
      if boundary.modified {
        modify_simple_timeline (self.accessor, handle, Some(boundary.varying));
      }
    }
  }
}*/

fn update_inferred_node_properties <A: EventAccessor <Steward = Steward >> (accessor: &A, node: &NodeHandle) {
  let mut varying = accessor.query (&node.varying, &GetVarying, QueryOffset::After).unwrap().1;
  varying.slope = [0, 0];
  varying.accumulation_rate = 0;
  for dimension in 0..2 {
    for direction in 0..2 {
      let direction_signum = (direction as Amount*2)-1;
      for boundary in varying.boundaries [dimension] [direction].iter() {
        let boundary_varying = get!(accessor, & boundary.varying);
        let transfer_rate = boundary_varying.transfer_velocity*boundary.length as Amount;
        varying.accumulation_rate -= transfer_rate*direction_signum;
        varying.slope [dimension] -= transfer_rate; // divided by
        // (2*node.width*VELOCITY_PER_SLOPE), but we do that later as an optimization
      }
    }
    varying.slope [dimension] /= (2*node.width as Amount)*VELOCITY_PER_SLOPE;
  }
  modify_simple_timeline (accessor, &node.varying, Some (varying));
}

fn update_node <A: EventAccessor <Steward = Steward >> (accessor: &A, node: &NodeHandle) {
  set_with!(accessor, &node.varying, | varying | {
    varying.ink_at_last_change = ink_at (accessor, node, *accessor.now());
    varying.last_change = *accessor.now();
  });
}

fn ink_at <A: Accessor <Steward = Steward >> (accessor: &A, node: &NodeHandle, time: Time)->Amount {
  let varying = get!(accessor, &node.varying);
  varying.ink_at_last_change + varying.accumulation_rate*((time - varying.last_change) as Amount)
}
fn density_at <A: Accessor <Steward = Steward >> (accessor: &A, node: &NodeHandle, time: Time)->Amount {
  ink_at (accessor, node, time)/((node.width*node.width) as Amount)
}

/*fn inferred_density_trajectory_at <A: Accessor <Steward = Steward >> (accessor: &A, node: &NodeHandle, location: [Amount; 2], time: Time)->[Amount; 2] {
  let area = node.width*node.width;
  [
    ink_at(node, time)/area
      + node.slope [0]*(location [0] - node.center [0])
      + node.slope [1]*(location [1] - node.center [1]),
    node.accumulation_rate/area
  ]
}*/


fn density_difference_to_ideal_transfer_velocity (nodes: [& NodeHandle; 2], density_difference: Amount)->Amount {
  //let ideal_transfer_slope = density_difference*2/((nodes [0].width + nodes [1].width) as Amount);
  //ideal_transfer_slope*VELOCITY_PER_SLOPE
  density_difference*(2*VELOCITY_PER_SLOPE)/((nodes [0].width + nodes [1].width) as Amount)
}
fn transfer_velocity_to_ideal_density_difference (nodes: [&NodeHandle;2], transfer_velocity: Amount)->Amount {
  //let ideal_transfer_slope = transfer_velocity/VELOCITY_PER_SLOPE;
  //ideal_transfer_slope*((nodes [0].width + nodes [1].width) as Amount)/2
  transfer_velocity*((nodes [0].width + nodes [1].width) as Amount)/(2*VELOCITY_PER_SLOPE)
}


fn safer_scaling (value: Amount, numerator: Amount, denominator: Amount)->Amount {
  if numerator.abs() > denominator.abs()*50 {
    value*(numerator/denominator)
  }
  else if value.abs() > denominator.abs()*50 {
    numerator*(value/denominator)
  }
  else if denominator.abs() > numerator.abs()*50 {
    value/(denominator/numerator)
  }
  else if denominator.abs() > value.abs()*50 {
    numerator/(denominator/value)
  }
  else {
    (value*numerator)/denominator
  }
}


fn update_transfer_change_prediction <A: EventAccessor <Steward = Steward >> (accessor: &A, boundary: &BoundaryHandle) {
  set_with!(accessor, &boundary.varying, | boundary_varying | {
  let nodes = &boundary.nodes;
  
  let areas = [
    (nodes [0].width*nodes [0].width) as Amount,
    (nodes [1].width*nodes [1].width) as Amount,
  ];
  let densities = [
    density_at (accessor, &nodes[0], *accessor.now()),
    density_at (accessor, &nodes[1], *accessor.now()),
  ];
  let accumulation_rates = [
    get!(accessor, &nodes[0].varying).accumulation_rate,
    get!(accessor, &nodes[1].varying).accumulation_rate,
  ];
  let scale_factor = areas [0]*areas [1];
  let density_accumulations_scaled = [
    accumulation_rates [0]*areas [1],
    accumulation_rates [1]*areas [0],
  ];
  let difference_now = densities[0] - densities[1];
  let density_accumulation_difference_scaled = density_accumulations_scaled [0] - density_accumulations_scaled [1];
  let numerator = accumulation_rates [0] * areas[1] - accumulation_rates [1] * areas[0];
  let denominator = 4*(nodes[0].width as Amount * areas[1] + nodes[1].width as Amount * areas[0]);
  let max_reasonable_velocity_skew = numerator/denominator;
  
  // the more unstable the boundary is, the more error we permit
  let max_rounding_error = max(VELOCITY_PER_SLOPE, (nodes [0].width + nodes [1].width) as Amount);
  let permissible_error = 2*max_rounding_error + max_reasonable_velocity_skew.abs();// + GENERIC_INK_AMOUNT/(SECOND as Amount)/(METER as Amount)/100000000;// + density_accumulation_difference.abs()*boundary.length as Amount/8;
  
  // Note to self: it may look weird for the boundary conditions to be based on
  // the CURRENT transfer_velocity, when it can change, but remember
  // that "when the current gets too far from the ideal" is the same as
  // "when the ideal gets too far from the current"
  let (min_difference, max_difference) = (
    transfer_velocity_to_ideal_density_difference([&nodes[0], &nodes[1]], boundary_varying.transfer_velocity - permissible_error),
    transfer_velocity_to_ideal_density_difference([&nodes[0], &nodes[1]], boundary_varying.transfer_velocity + permissible_error),
  );
  let time = if difference_now < min_difference || difference_now > max_difference {
    Some(*accessor.now())
  } else if density_accumulation_difference_scaled > 0 {
    Some(*accessor.now() + safer_scaling (max_difference-difference_now, scale_factor, density_accumulation_difference_scaled) as Time)
  }
  else if density_accumulation_difference_scaled < 0 {
    Some(*accessor.now() + safer_scaling (min_difference-difference_now, scale_factor, density_accumulation_difference_scaled) as Time)
  }
  else {
    None
  };
  //printlnerr!("{:?}", (density_accumulation_difference.signum(), max_reasonable_velocity_skew.signum()));
  //printlnerr!("{:?}", (*accessor.now()*1000/SECOND, time.map(|time|time*1000/SECOND), difference_now, density_accumulation_difference_scaled*SECOND as Amount/scale_factor, max_reasonable_velocity_skew*boundary.length as Amount / areas[0], - max_reasonable_velocity_skew*boundary.length as Amount / areas[1], (min_difference, difference_now, max_difference), boundary_varying.transfer_velocity, density_difference_to_ideal_transfer_velocity ([&nodes[0], &nodes[1]], difference_now), max_reasonable_velocity_skew, permissible_error));
  
  if let Some (discarded) = boundary_varying.next_change.take() {accessor.destroy_prediction (&discarded);}
  boundary_varying.next_change = time.map (|time| {
    accessor.create_prediction (
      time,
      DeterministicRandomId::new (&(accessor.id(), boundary.center)),
      TransferChange {boundary: boundary.clone()}
    )
  });
  
  });
}

fn maybe_split <A: EventAccessor <Steward = Steward >> (accessor: &A, node: &NodeHandle) {
  if node.width <METER {return;}
  let mut varying = get!(accessor, &node.varying);
  let mut averaged_ideal_velocities = [0; 2];
  for dimension in 0..2 {
    for direction in 0..2 {
      let direction_signum = (direction as Amount*2)-1;
      for boundary in varying.boundaries [dimension] [direction].iter() {
        let density_difference =
          density_at (accessor, &boundary.nodes[0], *accessor.now()) - 
          density_at (accessor, &boundary.nodes[1], *accessor.now());
        let ideal_velocity = density_difference_to_ideal_transfer_velocity ([&boundary.nodes[0], &boundary.nodes[0]], density_difference);
        let ideal_rate = ideal_velocity*boundary.length as Amount;
        averaged_ideal_velocities [dimension] -= ideal_rate;
      }
    }
    averaged_ideal_velocities [dimension] /= 2*node.width as Amount;
  }
  let mut close_enough = true;
  for dimension in 0..2 {
    for direction in 0..2 {
      for boundary in varying.boundaries [dimension] [direction].iter() {
        let density_difference =
          density_at (accessor, &boundary.nodes[0], *accessor.now()) - 
          density_at (accessor, &boundary.nodes[1], *accessor.now());
        let ideal_velocity = density_difference_to_ideal_transfer_velocity ([&boundary.nodes[0], &boundary.nodes[0]], density_difference);
        let averaged_velocity = averaged_ideal_velocities [dimension];
        let discrepancy = (ideal_velocity - averaged_velocity).abs();
        // TODO: should this be proportional in some way to the size of the node(s)?
        if discrepancy*(boundary.length as Amount) > GENERIC_INK_AMOUNT/(SECOND as Amount) {
          close_enough = false;
        }
      }
    }
  }
  
  if !close_enough {
    split (accessor, node);
  }
}

fn split <A: EventAccessor <Steward = Steward >> (accessor: &A, node: &NodeHandle) {
  let mut varying = get!(accessor, &node.varying);
  //printlnerr!("{:?}", (node.width, node.center, varying.children.len()));
  //if !varying.children.is_empty() {
    //printlnerr!("{:?}", (node.width, node.center, varying.children.len()));
  //}
  assert!(varying.children.is_empty());
  let mut velocities = [[[0;2];2];2];
  let mut old_boundaries = [[[None,None],[None,None]],[[None,None],[None,None]]];
  let mut middle_velocities = [[0;2];2];
  for dimension in 0..2 {
    for direction in 0..2 {
      //let direction_signum = (direction*2)-1;
      let other_direction = (direction + 1) & 1;
      
      match varying.boundaries [dimension] [direction].len() {
        1 => {
          loop {
            {
              let boundary = &varying.boundaries [dimension] [direction] [0];
              let neighbor = &boundary.nodes [direction];
              if neighbor.width <= node.width {
                break;
              }
              split (accessor, &neighbor);
            }
            varying = get!(accessor, &node.varying);
          }
          let boundary = &varying.boundaries [dimension] [direction] [0];
          assert!(&boundary.nodes[other_direction] == node);
          velocities [dimension] [direction] = [
            get!(accessor, &boundary.varying).transfer_velocity,
            get!(accessor, &boundary.varying).transfer_velocity
          ];
          old_boundaries [dimension] [direction] = [
            Some(boundary.clone()),
            Some(boundary.clone()),
          ];
        },
        2 => {
          let boundaries = &varying.boundaries [dimension] [direction];
          assert!(&boundaries[0].nodes[other_direction] == node);
          assert!(&boundaries[1].nodes[other_direction] == node);
          velocities [dimension] [direction] = [
            get!(accessor, &boundaries [0].varying).transfer_velocity,
            get!(accessor, &boundaries [1].varying).transfer_velocity,
          ];
          old_boundaries [dimension] [direction] = [
            Some(boundaries [0].clone()),
            Some(boundaries [1].clone()),
          ];
        },
        0 => (),
        _ => unreachable!(),
      };
    }
    middle_velocities [dimension] = [
      (velocities [dimension] [0] [0] + velocities [dimension] [1] [0]) / 2,
      (velocities [dimension] [0] [1] + velocities [dimension] [1] [1]) / 2,
    ];
  }
  
  // TODO: more accurate computation of ideal middle velocities and ink amounts,
  // and conserve ink correctly
  let mut new_children = [[[None,None],[None,None]],[[None,None],[None,None]]];
  for x in 0..2 {
    let x_signum = (x as Distance*2)-1;
    for y in 0..2 {
      let y_signum = (y as Distance*2)-1;
      let center = [
        node.center[0] + (node.width >> 2)*x_signum,
        node.center[1] + (node.width >> 2)*y_signum,
      ];
      let new_child = DataHandle::new (NodeData {
        width: node.width >> 1,
        center: center,
        parent: Some (node.clone()),
        varying: DataTimelineCell::new (SimpleTimeline::new ()),
      });
      modify_simple_timeline (accessor, &new_child.varying, Some (NodeVarying {
        last_change: *accessor.now(),
        ink_at_last_change: varying.ink_at_last_change >> 2,
        children: Vec::new(),
        boundaries: [[Vec::new(), Vec::new()], [Vec::new(), Vec::new()]],
        accumulation_rate: 0,
        slope: [0; 2],
      }));
      varying.children.push (new_child.clone());
      new_children [0][x][y] = Some(new_child.clone());
      new_children [1][y][x] = Some(new_child.clone());
    }
  }
  for dimension in 0..2 {
    for direction in 0..2 {
      let other_direction = (direction + 1) & 1;
      for which in 0..2 {
        if let Some(boundary) = old_boundaries [dimension] [direction] [which].as_ref() {
          let neighbor = &boundary.nodes [direction];
          assert!(neighbor.width <= node.width);
          set_with!(accessor, &neighbor.varying, | neighbor_varying | {
            for b in neighbor_varying.boundaries [dimension] [other_direction].iter() {
              assert!(&b.nodes[other_direction] == node);
            }
            neighbor_varying.boundaries [dimension] [other_direction].clear();
          });
        }
      }
    }
  }
  let mut new_boundaries = Vec::new();
  for dimension in 0..2 {
    let other_dimension = (dimension + 1) & 1;
    for direction in 0..2 {
      let other_direction = (direction + 1) & 1;
      let direction_signum = (direction as Distance*2)-1;
      for which in 0..2 {
        let which_signum = (which as Distance*2)-1;
        if let Some(boundary) = old_boundaries [dimension] [direction] [which].as_ref() {
          let child = new_children [dimension] [direction] [which].as_ref().unwrap();
          let neighbor = &boundary.nodes [direction];
          let mut center = node.center;
          center [dimension] += (node.width >> 2)*direction_signum;
          center [other_dimension] += (node.width >> 4)*which_signum;
          let new_boundary = DataHandle::new (BoundaryData {
            length: node.width >> 1,
            center: center,
            nodes: if direction == 0 { [neighbor.clone(), child.clone()] } else { [child.clone(), neighbor.clone()] },
            varying: DataTimelineCell::new (SimpleTimeline::new ()),
          });
          modify_simple_timeline (accessor, &new_boundary.varying, Some (BoundaryVarying {
            transfer_velocity: velocities [dimension] [direction] [which],
            next_change: None,
          }));
          set_with!(accessor, &child.varying, | child_varying | {
            child_varying.boundaries [dimension] [direction].push(new_boundary.clone());
          });
          set_with!(accessor, &neighbor.varying, | neighbor_varying | {
            neighbor_varying.boundaries [dimension] [other_direction].push(new_boundary.clone());
          });
          new_boundaries.push (new_boundary);
        }
      }
    }
  }
  for dimension in 0..2 {
    let other_dimension = (dimension + 1) & 1;
    for which in 0..2 {
      let which_signum = (which as Distance*2)-1;
      let child0 = new_children [dimension] [0] [which].as_ref().unwrap();
      let child1 = new_children [dimension] [1] [which].as_ref().unwrap();
      let mut center = node.center;
      center [other_dimension] += (node.width >> 4)*which_signum;
      let new_boundary = DataHandle::new (BoundaryData {
        length: node.width >> 1,
        center: center,
        nodes: [child0.clone(), child1.clone()],
        varying: DataTimelineCell::new (SimpleTimeline::new ()),
      });
      modify_simple_timeline (accessor, &new_boundary.varying, Some (BoundaryVarying {
        transfer_velocity: middle_velocities [dimension] [which],
        next_change: None,
      }));
      set_with!(accessor, &child0.varying, | child0_varying | {
        child0_varying.boundaries [dimension] [1].push(new_boundary.clone());
      });
      set_with!(accessor, &child1.varying, | child1_varying | {
        child1_varying.boundaries [dimension] [0].push(new_boundary.clone());
      });
      new_boundaries.push (new_boundary);
    }
  }
  for dimension in 0..2 {
    for direction in 0..2 {
      for which in 0..2 {
        if let Some(boundary) = old_boundaries [dimension] [direction] [which].as_ref() {
          if let Some ((_, boundary_varying)) = accessor.query (&boundary.varying, &GetVarying, QueryOffset::After) {
            if let Some (discarded) = boundary_varying.next_change {accessor.destroy_prediction (&discarded);}
            modify_simple_timeline (accessor, &boundary.varying, None);
          }
        }
      }
    }
  }
  
  varying.boundaries = [[Vec::new(), Vec::new()], [Vec::new(), Vec::new()]];
  
  modify_simple_timeline (accessor, &node.varying, Some(varying));
  
  for boundary in new_boundaries {
    update_transfer_change_prediction (accessor, &boundary);
  }
}


#[derive (Clone, PartialEq, Eq, Serialize, Deserialize, Debug)]
struct Globals {
  size: [Distance; 2],

  
  root: NodeHandle,
}
impl StewardData for Globals {}



/// The TransferChange event, as used above.
#[derive (Clone, PartialEq, Eq, Serialize, Deserialize, Debug)]
struct TransferChange {boundary: BoundaryHandle}
impl StewardData for TransferChange {}
impl PersistentlyIdentifiedType for TransferChange {
  const ID: PersistentTypeId = PersistentTypeId(0xd6621e9cfad1c765);
}
impl Event for TransferChange {
  type Steward = Steward;
  type ExecutionData = ();
  fn execute <Accessor: EventAccessor <Steward = Self::Steward>> (&self, accessor: &mut Accessor) {
    let nodes = &self.boundary.nodes;
    for node in nodes.iter() {
      update_node (accessor, node);
    }
    let areas = [
      (nodes [0].width*nodes [0].width) as Amount,
      (nodes [1].width*nodes [1].width) as Amount,
    ];
    let densities = [
      density_at (accessor, &nodes[0], *accessor.now()),
      density_at (accessor, &nodes[1], *accessor.now()),
    ];
    let mut accumulation_rates = [
      get!(accessor, &nodes[0].varying).accumulation_rate,
      get!(accessor, &nodes[1].varying).accumulation_rate,
    ];
    let difference_now = densities[0] - densities[1];
    let physics_ideal_velocity = density_difference_to_ideal_transfer_velocity ([&nodes[0], &nodes[1]], difference_now);
    set_with!(accessor, &self.boundary.varying, | boundary_varying | {
      // for the sake of this math, assume that the accumulation rates have already been updated based on the physics-wise ideal value
      accumulation_rates [0] -= (physics_ideal_velocity - boundary_varying.transfer_velocity)*self.boundary.length as Amount;
      accumulation_rates [1] += (physics_ideal_velocity - boundary_varying.transfer_velocity)*self.boundary.length as Amount;
      // Select a velocity that would neutralize the relative density-accumulation if the rest of each node's perimeter was shifted in a similar way.
      // (accumulation0 - shift*perimeter0)/w0^2 == (accumulation1 + shift*perimeter1)/w1^2.
      // After a little algebra, 
      // (accumulation0*w1^2-accumulation1*w0^2)/(perimeter0*w1^2+perimeter1*w0^2) == shift.
      let numerator = accumulation_rates [0] * areas[1] - accumulation_rates [1] * areas[0];
      let denominator = 4*(nodes[0].width as Amount * areas[1] + nodes[1].width as Amount * areas[0]);
      let max_reasonable_velocity_skew = numerator/denominator;
      assert!(
        (accumulation_rates [0]*areas[1]-accumulation_rates [1]*areas[0]).abs()
        >=
        (
         (accumulation_rates [0] - max_reasonable_velocity_skew*self.boundary.length as Amount)*areas[1]
        -(accumulation_rates [1] + max_reasonable_velocity_skew*self.boundary.length as Amount)*areas[0]
        ).abs()
      );
      //printlnerr!("{:?}", max_reasonable_velocity_skew);
      boundary_varying.transfer_velocity = physics_ideal_velocity + numerator.signum()+max_reasonable_velocity_skew/2;
    });
    for node in nodes.iter() {
      update_inferred_node_properties (accessor, node);
    }
    for node in nodes.iter() {
      for whatever in get!(accessor, &node.varying).boundaries.iter() {for something in whatever.iter() {for other_boundary in something.iter() {
        update_transfer_change_prediction (accessor, other_boundary);
      }}}
    }
    
    for node in nodes.iter() {
      if get!(accessor, &node.varying).children.is_empty() {
        maybe_split (accessor, node);
      }
    }
    
    // TODO: invalidation
  }

  fn undo <Accessor: FutureCleanupAccessor <Steward = Self::Steward>> (&self, _accessor: &mut Accessor, _: ()) {
    unimplemented!()
  }
}


#[derive (Clone, PartialEq, Eq, Debug)]
struct AddInk {coordinates: [Distance; 2], amount: Amount, accumulation: Amount}
serialization_cheat!([][AddInk]);
impl StewardData for AddInk {}
impl PersistentlyIdentifiedType for AddInk {
  const ID: PersistentTypeId = PersistentTypeId(0x3e6d029c3da8b9a2);
}
impl Event for AddInk {
  type Steward = Steward;
  type ExecutionData = ();
  fn execute <Accessor: EventAccessor <Steward = Self::Steward>> (&self, accessor: &mut Accessor) {
    let mut node = accessor.globals().root.clone();<commit_msg>Implemented automatic merging nodes<commit_after>#![feature (i128_type)]

extern crate serde;
#[macro_use]
extern crate serde_derive;

#[macro_use]
extern crate time_steward;

// Imports for the UI
#[macro_use]
extern crate glium;
#[path = "../dev-shared/emscripten_compatibility.rs"] mod emscripten_compatibility;
pub use emscripten_compatibility::canvas_click;

macro_rules! printlnerr(
    ($($arg:tt)*) => { {use std::io::Write;
        let r = writeln!(&mut ::std::io::stderr(), $($arg)*);
        r.expect("failed printing to stderr");
    } }
);


use std::cmp::{min, max};
use std::collections::HashSet;

use time_steward::{DeterministicRandomId};
use time_steward::rowless::api::{self, PersistentTypeId, PersistentlyIdentifiedType, ListedType, StewardData, QueryOffset, DataHandleTrait, DataTimelineCellTrait, ExtendedTime, Basics as BasicsTrait};
use time_steward::rowless::stewards::{simple_full as steward_module};
use steward_module::{TimeSteward, ConstructibleTimeSteward, IncrementalTimeSteward, Event, DataHandle, DataTimelineCell, EventHandle, Accessor, EventAccessor, FutureCleanupAccessor, SnapshotAccessor, simple_timeline};
use simple_timeline::{SimpleTimeline, GetVarying, IterateUniquelyOwnedPredictions, tracking_query, modify_simple_timeline, unmodify_simple_timeline};


type Time = i64;
type Distance = i64;
type Amount = i128;
const SECOND: Time = (1 as Time) << 20;
const METER: Distance = (1 as Distance) << 20;
const GENERIC_DENSITY: Amount = (SECOND as Amount) << 20;
const GENERIC_INK_AMOUNT: Amount = (METER as Amount)*(METER as Amount)*GENERIC_DENSITY;

// Velocity is "ink units per space unit per time unit"
// and slope is "ink density per space unit", i.e. "ink units per space unit^3"
// 
// Let's say we want flow across a meter^2 cell, from a density of GENERIC_INK_AMOUNT/m^2 to 0, to drain GENERIC_INK_AMOUNT per second. Then
//   slope = GENERIC_INK_AMOUNT/m^3
//   velocity = GENERIC_INK_AMOUNT/m*s
// so
const VELOCITY_PER_SLOPE: Amount = (METER as Amount)*(METER as Amount) / (SECOND as Amount);

type Steward = steward_module::Steward <Basics>;

macro_rules! serialization_cheat {
  (
    [$($bounds:tt)*]
    [$($concrete:tt)*]
  ) => {

/*thread_local! {
  static SERIALIZATION_CONTEXT: RefCell<Option <SerializationContext>> = RefCell::new (None);
  static DESERIALIZATION_CONTEXT: RefCell<Option <DeserializationContext>> = RefCell::new (None);
}*/


impl <$($bounds)*> $crate::serde::Serialize for $($concrete)* {
  fn serialize <S: $crate::serde::Serializer> (&self, _: S)->Result <S::Ok, S::Error> {
    unimplemented!()
  }
}

impl <'a, $($bounds)*> $crate::serde::Deserialize <'a> for $($concrete)* {
  fn deserialize <D: $crate::serde::Deserializer<'a>> (_: D)->Result <Self, D::Error> {
    unimplemented!()
  }
}

  };
}

#[derive (Clone, PartialEq, Eq, Serialize, Deserialize, Debug)]
struct NodeData {
  width: Distance,
  center: [Distance ; 2],
  parent: Option <NodeHandle>,
  varying: DataTimelineCell <SimpleTimeline <NodeVarying, Steward>>,
}
#[derive (Clone, PartialEq, Eq, Debug)]
struct NodeVarying {
  last_change: Time,
  ink_at_last_change: Amount,
  children: Vec<NodeHandle>,
  boundaries: [[Vec<BoundaryHandle>;2];2],
  
  accumulation_rate: Amount,
  slope: [Amount; 2],
}
type NodeHandle = DataHandle <NodeData>;
impl StewardData for NodeData {}
impl StewardData for NodeVarying {}
impl PersistentlyIdentifiedType for NodeData {
  const ID: PersistentTypeId = PersistentTypeId(0x734528f6aefdc1b9);
}
impl IterateUniquelyOwnedPredictions <Steward> for NodeVarying {}
serialization_cheat!([][NodeVarying]);

#[derive (Clone, PartialEq, Eq, Serialize, Deserialize, Debug)]
struct BoundaryData {
  length: Distance,
  center: [Distance ; 2],
  nodes: [NodeHandle; 2],
  varying: DataTimelineCell <SimpleTimeline <BoundaryVarying, Steward>>,
}
#[derive (Clone, PartialEq, Eq, Debug)]
struct BoundaryVarying {
  transfer_velocity: Amount,
  next_change: Option <EventHandle <Basics>>,
}
type BoundaryHandle = DataHandle <BoundaryData>;
impl StewardData for BoundaryData {}
impl StewardData for BoundaryVarying {}
impl PersistentlyIdentifiedType for BoundaryData {
  const ID: PersistentTypeId = PersistentTypeId(0x9d643214b58b24dc);
}
impl IterateUniquelyOwnedPredictions <Steward> for BoundaryVarying {
  fn iterate_predictions <F: FnMut (& <Steward as TimeSteward>::EventHandle)> (&self, callback: &mut F) {
    if let Some (prediction) = self.next_change.as_ref() {
      callback (prediction);
    }
  }
}
serialization_cheat!([][BoundaryVarying]);

macro_rules! get {
  ($accessor: expr, $cell: expr) => {
    $accessor.query ($cell, &GetVarying, QueryOffset::After).unwrap().1
  }
}
macro_rules! set {
  ($accessor: expr, $cell: expr, $field: ident, $value: expr) => {
    {
      let mut value = $accessor.query ($cell, &GetVarying, QueryOffset::After).unwrap().1;
      value.$field = $value;
      modify_simple_timeline ($accessor, $cell, Some (value));
    }
  }
}
macro_rules! set_with {
  ($accessor: expr, $cell: expr, | $varying: ident | $actions: expr) => {
    {
      let mut $varying = $accessor.query ($cell, &GetVarying, QueryOffset::After).unwrap().1;
      $actions;
      modify_simple_timeline ($accessor, $cell, Some ($varying));
    }
  }
}
macro_rules! exists {
  ($accessor: expr, $cell: expr) => {
    $accessor.query ($cell, &GetVarying, QueryOffset::After).is_some()
  }
}


/*struct EventContext <'a, A: EventAccessor <Steward = Steward>> {
  accessor: &'a A,
  node_clones: HashMap <DataHandle <NodeData>, Rc<RefCell<Node>>,
  boundary_clones: HashMap <DataHandle <NodeData>, Rc<RefCell<Node>>,
}
impl<A: EventAccessor <Steward = Steward>> EventContext <A> {
  fn new (accessor: & A)->Self {
    EventContext {
      accessor: accessor,
      node_clones: HashMap::new(),
      boundary_clones: HashMap::new(),
    }
  }
  fn get_node (&self, handle: & DataHandle <NodeData>)->Rc {
    self.node_clones.insert_or (self.accessor.query (
  }
}
impl<A: EventAccessor <Steward = Steward>> Drop for EventContext <A> {
  fn drop (&mut self) {
    for (handle, node) in self.node_clones.iter() {
      if node.modified {
        modify_simple_timeline (self.accessor, handle, Some(node.varying));
      }
    }
    for (handle, boundary) in self.node_clones.iter() {
      if boundary.modified {
        modify_simple_timeline (self.accessor, handle, Some(boundary.varying));
      }
    }
  }
}*/

fn update_inferred_node_properties <A: EventAccessor <Steward = Steward >> (accessor: &A, node: &NodeHandle) {
  let mut varying = accessor.query (&node.varying, &GetVarying, QueryOffset::After).unwrap().1;
  varying.slope = [0, 0];
  varying.accumulation_rate = 0;
  for dimension in 0..2 {
    for direction in 0..2 {
      let direction_signum = (direction as Amount*2)-1;
      for boundary in varying.boundaries [dimension] [direction].iter() {
        let boundary_varying = get!(accessor, & boundary.varying);
        let transfer_rate = boundary_varying.transfer_velocity*boundary.length as Amount;
        varying.accumulation_rate -= transfer_rate*direction_signum;
        varying.slope [dimension] -= transfer_rate; // divided by
        // (2*node.width*VELOCITY_PER_SLOPE), but we do that later as an optimization
      }
    }
    varying.slope [dimension] /= (2*node.width as Amount)*VELOCITY_PER_SLOPE;
  }
  modify_simple_timeline (accessor, &node.varying, Some (varying));
}

fn update_node <A: EventAccessor <Steward = Steward >> (accessor: &A, node: &NodeHandle) {
  set_with!(accessor, &node.varying, | varying | {
    varying.ink_at_last_change = ink_at (accessor, node, *accessor.now());
    varying.last_change = *accessor.now();
  });
}

fn ink_at <A: Accessor <Steward = Steward >> (accessor: &A, node: &NodeHandle, time: Time)->Amount {
  let varying = get!(accessor, &node.varying);
  varying.ink_at_last_change + varying.accumulation_rate*((time - varying.last_change) as Amount)
}
fn density_at <A: Accessor <Steward = Steward >> (accessor: &A, node: &NodeHandle, time: Time)->Amount {
  ink_at (accessor, node, time)/((node.width*node.width) as Amount)
}

/*fn inferred_density_trajectory_at <A: Accessor <Steward = Steward >> (accessor: &A, node: &NodeHandle, location: [Amount; 2], time: Time)->[Amount; 2] {
  let area = node.width*node.width;
  [
    ink_at(node, time)/area
      + node.slope [0]*(location [0] - node.center [0])
      + node.slope [1]*(location [1] - node.center [1]),
    node.accumulation_rate/area
  ]
}*/


fn density_difference_to_ideal_transfer_velocity (nodes: [& NodeHandle; 2], density_difference: Amount)->Amount {
  //let ideal_transfer_slope = density_difference*2/((nodes [0].width + nodes [1].width) as Amount);
  //ideal_transfer_slope*VELOCITY_PER_SLOPE
  density_difference*(2*VELOCITY_PER_SLOPE)/((nodes [0].width + nodes [1].width) as Amount)
}
fn transfer_velocity_to_ideal_density_difference (nodes: [&NodeHandle;2], transfer_velocity: Amount)->Amount {
  //let ideal_transfer_slope = transfer_velocity/VELOCITY_PER_SLOPE;
  //ideal_transfer_slope*((nodes [0].width + nodes [1].width) as Amount)/2
  transfer_velocity*((nodes [0].width + nodes [1].width) as Amount)/(2*VELOCITY_PER_SLOPE)
}


fn safer_scaling (value: Amount, numerator: Amount, denominator: Amount)->Amount {
  if numerator.abs() > denominator.abs()*50 {
    value*(numerator/denominator)
  }
  else if value.abs() > denominator.abs()*50 {
    numerator*(value/denominator)
  }
  else if denominator.abs() > numerator.abs()*50 {
    value/(denominator/numerator)
  }
  else if denominator.abs() > value.abs()*50 {
    numerator/(denominator/value)
  }
  else {
    (value*numerator)/denominator
  }
}


fn update_transfer_change_prediction <A: EventAccessor <Steward = Steward >> (accessor: &A, boundary: &BoundaryHandle) {
  set_with!(accessor, &boundary.varying, | boundary_varying | {
  let nodes = &boundary.nodes;
  
  let areas = [
    (nodes [0].width*nodes [0].width) as Amount,
    (nodes [1].width*nodes [1].width) as Amount,
  ];
  let densities = [
    density_at (accessor, &nodes[0], *accessor.now()),
    density_at (accessor, &nodes[1], *accessor.now()),
  ];
  let accumulation_rates = [
    get!(accessor, &nodes[0].varying).accumulation_rate,
    get!(accessor, &nodes[1].varying).accumulation_rate,
  ];
  let scale_factor = areas [0]*areas [1];
  let density_accumulations_scaled = [
    accumulation_rates [0]*areas [1],
    accumulation_rates [1]*areas [0],
  ];
  let difference_now = densities[0] - densities[1];
  let density_accumulation_difference_scaled = density_accumulations_scaled [0] - density_accumulations_scaled [1];
  let numerator = accumulation_rates [0] * areas[1] - accumulation_rates [1] * areas[0];
  let denominator = 4*(nodes[0].width as Amount * areas[1] + nodes[1].width as Amount * areas[0]);
  let max_reasonable_velocity_skew = numerator/denominator;
  
  // the more unstable the boundary is, the more error we permit
  let max_rounding_error = max(VELOCITY_PER_SLOPE, (nodes [0].width + nodes [1].width) as Amount);
  let permissible_error = 2*max_rounding_error + max_reasonable_velocity_skew.abs() + GENERIC_INK_AMOUNT/(SECOND as Amount)/(METER as Amount)/1000000;// + density_accumulation_difference.abs()*boundary.length as Amount/8;
  
  // Note to self: it may look weird for the boundary conditions to be based on
  // the CURRENT transfer_velocity, when it can change, but remember
  // that "when the current gets too far from the ideal" is the same as
  // "when the ideal gets too far from the current"
  let (min_difference, max_difference) = (
    transfer_velocity_to_ideal_density_difference([&nodes[0], &nodes[1]], boundary_varying.transfer_velocity - permissible_error),
    transfer_velocity_to_ideal_density_difference([&nodes[0], &nodes[1]], boundary_varying.transfer_velocity + permissible_error),
  );
  let time = if difference_now < min_difference || difference_now > max_difference {
    Some(*accessor.now())
  } else if density_accumulation_difference_scaled > 0 {
    Some(*accessor.now() + safer_scaling (max_difference-difference_now, scale_factor, density_accumulation_difference_scaled) as Time)
  }
  else if density_accumulation_difference_scaled < 0 {
    Some(*accessor.now() + safer_scaling (min_difference-difference_now, scale_factor, density_accumulation_difference_scaled) as Time)
  }
  else {
    None
  };
  //printlnerr!("{:?}", (density_accumulation_difference.signum(), max_reasonable_velocity_skew.signum()));
  //printlnerr!("{:?}", (*accessor.now()*1000/SECOND, time.map(|time|time*1000/SECOND), difference_now, density_accumulation_difference_scaled*SECOND as Amount/scale_factor, max_reasonable_velocity_skew*boundary.length as Amount / areas[0], - max_reasonable_velocity_skew*boundary.length as Amount / areas[1], (min_difference, difference_now, max_difference), boundary_varying.transfer_velocity, density_difference_to_ideal_transfer_velocity ([&nodes[0], &nodes[1]], difference_now), max_reasonable_velocity_skew, permissible_error));
  
  if let Some (discarded) = boundary_varying.next_change.take() {accessor.destroy_prediction (&discarded);}
  boundary_varying.next_change = time.map (|time| {
    accessor.create_prediction (
      time,
      DeterministicRandomId::new (&(accessor.id(), boundary.center)),
      TransferChange {boundary: boundary.clone()}
    )
  });
  
  });
}


fn audit <A: EventAccessor <Steward = Steward >> (accessor: &A) {
  audit_node (accessor, &accessor.globals().root);
}
fn audit_node <A: EventAccessor <Steward = Steward >> (accessor: &A, node: &NodeHandle) {
  let varying = get!(accessor, &node.varying);
  for dimension in 0..2 {
    for direction in 0..2 {
      for boundary in varying.boundaries [dimension] [direction].iter() {
        assert!(&boundary.nodes[(direction+1)&1] == node);
        audit_boundary (accessor, boundary);
      }
    }
  }
  for child in varying.children.iter() {
    audit_node (accessor, child);
  }
}
fn audit_boundary <A: EventAccessor <Steward = Steward >> (accessor: &A, boundary: &BoundaryHandle) {
  get!(accessor, &boundary.varying);
  for direction in 0..2 {
    assert!(get!(accessor, &boundary.nodes [direction].varying).boundaries.iter().any(|whatever| whatever[(direction+1)&1].iter().any(| other_boundary| other_boundary == boundary)));
  }
}

fn maybe_split <A: EventAccessor <Steward = Steward >> (accessor: &A, node: &NodeHandle)->bool {
  if node.width <METER {return false;}
  let varying = get!(accessor, &node.varying);
  let mut averaged_ideal_velocities = [0; 2];
  for dimension in 0..2 {
    for direction in 0..2 {
      for boundary in varying.boundaries [dimension] [direction].iter() {
        let density_difference =
          density_at (accessor, &boundary.nodes[0], *accessor.now()) - 
          density_at (accessor, &boundary.nodes[1], *accessor.now());
        let ideal_velocity = density_difference_to_ideal_transfer_velocity ([&boundary.nodes[0], &boundary.nodes[0]], density_difference);
        let ideal_rate = ideal_velocity*boundary.length as Amount;
        averaged_ideal_velocities [dimension] -= ideal_rate;
      }
    }
    averaged_ideal_velocities [dimension] /= 2*node.width as Amount;
  }
  let mut close_enough = true;
  for dimension in 0..2 {
    for direction in 0..2 {
      for boundary in varying.boundaries [dimension] [direction].iter() {
        let density_difference =
          density_at (accessor, &boundary.nodes[0], *accessor.now()) - 
          density_at (accessor, &boundary.nodes[1], *accessor.now());
        let ideal_velocity = density_difference_to_ideal_transfer_velocity ([&boundary.nodes[0], &boundary.nodes[0]], density_difference);
        let averaged_velocity = averaged_ideal_velocities [dimension];
        let discrepancy = (ideal_velocity - averaged_velocity).abs();
        // TODO: should this be proportional in some way to the size of the node(s)?
        if discrepancy*(boundary.length as Amount) > GENERIC_INK_AMOUNT/(SECOND as Amount) {
          close_enough = false;
        }
      }
    }
  }
  
  if !close_enough {
    split (accessor, node);
  }
  !close_enough
}

fn split <A: EventAccessor <Steward = Steward >> (accessor: &A, node: &NodeHandle) {
  let mut varying = get!(accessor, &node.varying);
  //printlnerr!("{:?}", (node.width, node.center, varying.children.len()));
  //if !varying.children.is_empty() {
    //printlnerr!("{:?}", (node.width, node.center, varying.children.len()));
  //}
  assert!(varying.children.is_empty());
  let mut velocities = [[[0;2];2];2];
  let mut old_boundaries = [[[None,None],[None,None]],[[None,None],[None,None]]];
  let mut middle_velocities = [[0;2];2];
  for dimension in 0..2 {
    for direction in 0..2 {
      //let direction_signum = (direction*2)-1;
      let other_direction = (direction + 1) & 1;
      
      match varying.boundaries [dimension] [direction].len() {
        1 => {
          loop {
            {
              let boundary = &varying.boundaries [dimension] [direction] [0];
              let neighbor = &boundary.nodes [direction];
              if neighbor.width <= node.width {
                break;
              }
              split (accessor, &neighbor);
            }
            varying = get!(accessor, &node.varying);
          }
          let boundary = &varying.boundaries [dimension] [direction] [0];
          assert!(&boundary.nodes[other_direction] == node);
          velocities [dimension] [direction] = [
            get!(accessor, &boundary.varying).transfer_velocity,
            get!(accessor, &boundary.varying).transfer_velocity
          ];
          old_boundaries [dimension] [direction] = [
            Some(boundary.clone()),
            Some(boundary.clone()),
          ];
        },
        2 => {
          let boundaries = &varying.boundaries [dimension] [direction];
          assert!(&boundaries[0].nodes[other_direction] == node);
          assert!(&boundaries[1].nodes[other_direction] == node);
          velocities [dimension] [direction] = [
            get!(accessor, &boundaries [0].varying).transfer_velocity,
            get!(accessor, &boundaries [1].varying).transfer_velocity,
          ];
          old_boundaries [dimension] [direction] = [
            Some(boundaries [0].clone()),
            Some(boundaries [1].clone()),
          ];
        },
        0 => (),
        _ => unreachable!(),
      };
    }
    middle_velocities [dimension] = [
      (velocities [dimension] [0] [0] + velocities [dimension] [1] [0]) / 2,
      (velocities [dimension] [0] [1] + velocities [dimension] [1] [1]) / 2,
    ];
  }
  
  // TODO: more accurate computation of ideal middle velocities and ink amounts,
  // and conserve ink correctly
  let mut new_children = [[[None,None],[None,None]],[[None,None],[None,None]]];
  for x in 0..2 {
    let x_signum = (x as Distance*2)-1;
    for y in 0..2 {
      let y_signum = (y as Distance*2)-1;
      let center = [
        node.center[0] + (node.width >> 2)*x_signum,
        node.center[1] + (node.width >> 2)*y_signum,
      ];
      let new_child = DataHandle::new (NodeData {
        width: node.width >> 1,
        center: center,
        parent: Some (node.clone()),
        varying: DataTimelineCell::new (SimpleTimeline::new ()),
      });
      modify_simple_timeline (accessor, &new_child.varying, Some (NodeVarying {
        last_change: *accessor.now(),
        ink_at_last_change: varying.ink_at_last_change >> 2,
        children: Vec::new(),
        boundaries: [[Vec::new(), Vec::new()], [Vec::new(), Vec::new()]],
        accumulation_rate: 0,
        slope: [0; 2],
      }));
      varying.children.push (new_child.clone());
      new_children [0][x][y] = Some(new_child.clone());
      new_children [1][y][x] = Some(new_child.clone());
    }
  }
  for dimension in 0..2 {
    for direction in 0..2 {
      let other_direction = (direction + 1) & 1;
      for which in 0..2 {
        if let Some(boundary) = old_boundaries [dimension] [direction] [which].as_ref() {
          let neighbor = &boundary.nodes [direction];
          assert!(neighbor.width <= node.width);
          set_with!(accessor, &neighbor.varying, | neighbor_varying | {
            for b in neighbor_varying.boundaries [dimension] [other_direction].iter() {
              assert!(&b.nodes[other_direction] == node);
            }
            neighbor_varying.boundaries [dimension] [other_direction].clear();
          });
        }
      }
    }
  }
  let mut new_boundaries = Vec::new();
  for dimension in 0..2 {
    let other_dimension = (dimension + 1) & 1;
    for direction in 0..2 {
      let other_direction = (direction + 1) & 1;
      let direction_signum = (direction as Distance*2)-1;
      for which in 0..2 {
        let which_signum = (which as Distance*2)-1;
        if let Some(boundary) = old_boundaries [dimension] [direction] [which].as_ref() {
          let child = new_children [dimension] [direction] [which].as_ref().unwrap();
          let neighbor = &boundary.nodes [direction];
          let mut center = node.center;
          center [dimension] += (node.width >> 2)*direction_signum;
          center [other_dimension] += (node.width >> 4)*which_signum;
          let new_boundary = DataHandle::new (BoundaryData {
            length: node.width >> 1,
            center: center,
            nodes: if direction == 0 { [neighbor.clone(), child.clone()] } else { [child.clone(), neighbor.clone()] },
            varying: DataTimelineCell::new (SimpleTimeline::new ()),
          });
          modify_simple_timeline (accessor, &new_boundary.varying, Some (BoundaryVarying {
            transfer_velocity: velocities [dimension] [direction] [which],
            next_change: None,
          }));
          set_with!(accessor, &child.varying, | child_varying | {
            child_varying.boundaries [dimension] [direction].push(new_boundary.clone());
          });
          set_with!(accessor, &neighbor.varying, | neighbor_varying | {
            neighbor_varying.boundaries [dimension] [other_direction].push(new_boundary.clone());
          });
          new_boundaries.push (new_boundary);
        }
      }
    }
  }
  for dimension in 0..2 {
    let other_dimension = (dimension + 1) & 1;
    for which in 0..2 {
      let which_signum = (which as Distance*2)-1;
      let child0 = new_children [dimension] [0] [which].as_ref().unwrap();
      let child1 = new_children [dimension] [1] [which].as_ref().unwrap();
      let mut center = node.center;
      center [other_dimension] += (node.width >> 4)*which_signum;
      let new_boundary = DataHandle::new (BoundaryData {
        length: node.width >> 1,
        center: center,
        nodes: [child0.clone(), child1.clone()],
        varying: DataTimelineCell::new (SimpleTimeline::new ()),
      });
      modify_simple_timeline (accessor, &new_boundary.varying, Some (BoundaryVarying {
        transfer_velocity: middle_velocities [dimension] [which],
        next_change: None,
      }));
      set_with!(accessor, &child0.varying, | child0_varying | {
        child0_varying.boundaries [dimension] [1].push(new_boundary.clone());
      });
      set_with!(accessor, &child1.varying, | child1_varying | {
        child1_varying.boundaries [dimension] [0].push(new_boundary.clone());
      });
      new_boundaries.push (new_boundary);
    }
  }
  for dimension in 0..2 {
    for direction in 0..2 {
      for which in 0..2 {
        if let Some(boundary) = old_boundaries [dimension] [direction] [which].as_ref() {
          if let Some ((_, boundary_varying)) = accessor.query (&boundary.varying, &GetVarying, QueryOffset::After) {
            if let Some (discarded) = boundary_varying.next_change {accessor.destroy_prediction (&discarded);}
            modify_simple_timeline (accessor, &boundary.varying, None);
          }
        }
      }
    }
  }
  
  varying.boundaries = [[Vec::new(), Vec::new()], [Vec::new(), Vec::new()]];
  
  modify_simple_timeline (accessor, &node.varying, Some(varying));
  
  for boundary in new_boundaries {
    update_transfer_change_prediction (accessor, &boundary);
  }
  
}


fn maybe_merge <A: EventAccessor <Steward = Steward >> (accessor: &A, node: &NodeHandle)->bool {
  if !exists!(accessor, &node.varying) { return false; }
  let varying = get!(accessor, &node.varying);
  if varying.children.is_empty() {
    for child in varying.children.iter() {
      if maybe_merge(accessor, child) { return true; }
    }
    return false;
  }
  let mut averaged_ideal_velocities = [0; 2];
  for dimension in 0..2 {
    for direction in 0..2 {
      for child in varying.children.iter() {
        let child_varying = get!(accessor, &child.varying);
        if child_varying.children.len() >0 {return false;}
        if child_varying.boundaries [dimension] [direction].len() >1 {return false;}
        for boundary in child_varying.boundaries [dimension] [direction].iter() {
          let density_difference =
            density_at (accessor, &boundary.nodes[0], *accessor.now()) - 
            density_at (accessor, &boundary.nodes[1], *accessor.now());
          let ideal_velocity = density_difference_to_ideal_transfer_velocity ([&boundary.nodes[0], &boundary.nodes[0]], density_difference);
          let ideal_rate = ideal_velocity*boundary.length as Amount;
          averaged_ideal_velocities [dimension] -= ideal_rate;
        }
      }
    }
    averaged_ideal_velocities [dimension] /= 4*node.width as Amount;
  }
  let mut close_enough = true;
  for dimension in 0..2 {
    for direction in 0..2 {
      for child in varying.children.iter() {
        let child_varying = get!(accessor, &child.varying);
        for boundary in child_varying.boundaries [dimension] [direction].iter() {
          let density_difference =
            density_at (accessor, &boundary.nodes[0], *accessor.now()) - 
            density_at (accessor, &boundary.nodes[1], *accessor.now());
          let ideal_velocity = density_difference_to_ideal_transfer_velocity ([&boundary.nodes[0], &boundary.nodes[0]], density_difference);
          let averaged_velocity = averaged_ideal_velocities [dimension];
          let discrepancy = (ideal_velocity - averaged_velocity).abs();
          // TODO: should this be proportional in some way to the size of the node(s)?
          if discrepancy*(boundary.length as Amount) > GENERIC_INK_AMOUNT/(SECOND as Amount)/4 {
            close_enough = false;
          }
        }
      }
    }
  }
  
  if close_enough {
    merge (accessor, node);
  }
  close_enough
}

fn merge <A: EventAccessor <Steward = Steward >> (accessor: &A, node: &NodeHandle) {
  let mut varying = get!(accessor, &node.varying);
  assert!(!varying.children.is_empty());
  
  varying.last_change =*accessor.now();
  varying.ink_at_last_change = 0;
  let mut prior_boundaries = Vec::new();
  let mut discovered_prior_boundaries = HashSet::new();
  let mut prior_children = HashSet::new();
  let mut neighbors = Vec::new();
  let mut discovered_neighbors = HashSet::new();
  for child in varying.children.iter() {
    update_node(accessor, child);
    let child_varying = get!(accessor, &child.varying);
    assert!(child_varying.children.is_empty());
    assert!(child_varying.last_change == *accessor.now());
    varying.ink_at_last_change += child_varying.ink_at_last_change;
    for (dimension, whatever) in child_varying.boundaries.iter().enumerate() {for (direction, something) in whatever.iter().enumerate() {for other_boundary in something.iter() {
      if discovered_prior_boundaries.insert (other_boundary.clone()) {
        prior_boundaries.push ((dimension, direction, other_boundary.clone()));
      }
    }}}
    prior_children.insert (child.clone()) ;

    modify_simple_timeline (accessor, &child.varying, None);
  }
  
  let mut new_boundaries = Vec::new();
  for (dimension, direction, boundary) in prior_boundaries {
    let boundary_varying = get!(accessor, &boundary.varying);
    if let Some (discarded) = boundary_varying.next_change {accessor.destroy_prediction (&discarded);}
    let neighbor = if !prior_children.contains (&boundary.nodes [1]) {
      Some (boundary.nodes [1].clone())
    } else if !prior_children.contains (&boundary.nodes [0]) {
      Some (boundary.nodes [0].clone())
    } else {None};
    if let Some(neighbor) = neighbor {
      if discovered_neighbors.insert (neighbor.clone()) {
        neighbors.push (neighbor.clone());
        let other_direction = (direction+1)&1;
        set_with!(accessor, &neighbor.varying, | neighbor_varying | {
          for b in neighbor_varying.boundaries [dimension] [other_direction].iter() {
            assert!(prior_children.contains(&b.nodes[other_direction]));
          }
          neighbor_varying.boundaries [dimension] [other_direction].clear();
        });
        let mut center = boundary.center;
        if neighbor.width == node.width {
          center = [
            (node.center[0] + neighbor.center [0]) >> 1,
            (node.center[1] + neighbor.center [1]) >> 1,
          ];
        }
        let nodes = if direction == 1 {[node.clone(), neighbor.clone()]}else{[neighbor.clone(), node.clone()]};
        let new_boundary = DataHandle::new (BoundaryData {
          length: neighbor.width,
          center: center,
          nodes: nodes.clone(),
          varying: DataTimelineCell::new (SimpleTimeline::new ()),
        });
        modify_simple_timeline (accessor, &new_boundary.varying, Some (BoundaryVarying {
          transfer_velocity: boundary_varying.transfer_velocity,
          next_change: None,
        }));
        set_with!(accessor, &neighbor.varying, | neighbor_varying | {
          neighbor_varying.boundaries [dimension] [other_direction].push(new_boundary.clone());
        });
        varying.boundaries [dimension] [direction].push(new_boundary.clone());
        new_boundaries.push (new_boundary);
      }
    }
    modify_simple_timeline (accessor, &boundary.varying, None);
  }
  
  varying.children.clear();
  modify_simple_timeline (accessor, &node.varying, Some(varying));
  update_inferred_node_properties (accessor, node);
  
  for boundary in new_boundaries {
    update_transfer_change_prediction (accessor, &boundary);
  }

  if let Some(parent) = node.parent.as_ref() {
    //audit(accessor);
    maybe_merge (accessor, parent);
    //audit(accessor);
  }
  
  for neighbor in neighbors {
    maybe_merge (accessor, &neighbor);
  }
}


#[derive (Clone, PartialEq, Eq, Serialize, Deserialize, Debug)]
struct Globals {
  size: [Distance; 2],

  
  root: NodeHandle,
}
impl StewardData for Globals {}



/// The TransferChange event, as used above.
#[derive (Clone, PartialEq, Eq, Serialize, Deserialize, Debug)]
struct TransferChange {boundary: BoundaryHandle}
impl StewardData for TransferChange {}
impl PersistentlyIdentifiedType for TransferChange {
  const ID: PersistentTypeId = PersistentTypeId(0xd6621e9cfad1c765);
}
impl Event for TransferChange {
  type Steward = Steward;
  type ExecutionData = ();
  fn execute <Accessor: EventAccessor <Steward = Self::Steward>> (&self, accessor: &mut Accessor) {
    let nodes = &self.boundary.nodes;
    for node in nodes.iter() {
      update_node (accessor, node);
    }
    let areas = [
      (nodes [0].width*nodes [0].width) as Amount,
      (nodes [1].width*nodes [1].width) as Amount,
    ];
    let densities = [
      density_at (accessor, &nodes[0], *accessor.now()),
      density_at (accessor, &nodes[1], *accessor.now()),
    ];
    let mut accumulation_rates = [
      get!(accessor, &nodes[0].varying).accumulation_rate,
      get!(accessor, &nodes[1].varying).accumulation_rate,
    ];
    let difference_now = densities[0] - densities[1];
    let physics_ideal_velocity = density_difference_to_ideal_transfer_velocity ([&nodes[0], &nodes[1]], difference_now);
    set_with!(accessor, &self.boundary.varying, | boundary_varying | {
      // for the sake of this math, assume that the accumulation rates have already been updated based on the physics-wise ideal value
      accumulation_rates [0] -= (physics_ideal_velocity - boundary_varying.transfer_velocity)*self.boundary.length as Amount;
      accumulation_rates [1] += (physics_ideal_velocity - boundary_varying.transfer_velocity)*self.boundary.length as Amount;
      // Select a velocity that would neutralize the relative density-accumulation if the rest of each node's perimeter was shifted in a similar way.
      // (accumulation0 - shift*perimeter0)/w0^2 == (accumulation1 + shift*perimeter1)/w1^2.
      // After a little algebra, 
      // (accumulation0*w1^2-accumulation1*w0^2)/(perimeter0*w1^2+perimeter1*w0^2) == shift.
      let numerator = accumulation_rates [0] * areas[1] - accumulation_rates [1] * areas[0];
      let denominator = 4*(nodes[0].width as Amount * areas[1] + nodes[1].width as Amount * areas[0]);
      let max_reasonable_velocity_skew = numerator/denominator;
      assert!(
        (accumulation_rates [0]*areas[1]-accumulation_rates [1]*areas[0]).abs()
        >=
        (
         (accumulation_rates [0] - max_reasonable_velocity_skew*self.boundary.length as Amount)*areas[1]
        -(accumulation_rates [1] + max_reasonable_velocity_skew*self.boundary.length as Amount)*areas[0]
        ).abs()
      );
      //printlnerr!("{:?}", max_reasonable_velocity_skew);
      boundary_varying.transfer_velocity = physics_ideal_velocity + numerator.signum()+max_reasonable_velocity_skew/2;
    });
    for node in nodes.iter() {
      update_inferred_node_properties (accessor, node);
    }
    for node in nodes.iter() {
      for whatever in get!(accessor, &node.varying).boundaries.iter() {for something in whatever.iter() {for other_boundary in something.iter() {
        update_transfer_change_prediction (accessor, other_boundary);
      }}}
    }
    
    let mut split = false;
    for node in nodes.iter() {
      let varying = get!(accessor, &node.varying);
      if varying.children.is_empty() {
        //audit(accessor);
        if maybe_split (accessor, node) {split=true;}
        //audit(accessor);
      }
    }
    if !split {
      for node in nodes.iter() {
        if let Some(parent) = node.parent.as_ref() {
          //audit(accessor);
          maybe_merge (accessor, parent);
          //audit(accessor);
        }
      }
    }
    
    // TODO: invalidation
  }

  fn undo <Accessor: FutureCleanupAccessor <Steward = Self::Steward>> (&self, _accessor: &mut Accessor, _: ()) {
    unimplemented!()
  }
}


#[derive (Clone, PartialEq, Eq, Debug)]
struct AddInk {coordinates: [Distance; 2], amount: Amount, accumulation: Amount}
serialization_cheat!([][AddInk]);
impl StewardData for AddInk {}
impl PersistentlyIdentifiedType for AddInk {
  const ID: PersistentTypeId = PersistentTypeId(0x3e6d029c3da8b9a2);
}
impl Event for AddInk {
  type Steward = Steward;
  type ExecutionData = ();
  fn execute <Accessor: EventAccessor <Steward = Self::Steward>> (&self, accessor: &mut Accessor) {
    let mut node = accessor.globals().root.clone();
{'fraction_non_alphanumeric': 0.13725048656125846, 'fraction_numerical': 0.012256518312396455, 'mean_word_length': 3.637970906378217, 'pattern_counts': {'":': 0, '<': 169, '<?xml version=': 0, '>': 243, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2950519', 'n_tokens_mistral': 21000, 'n_tokens_neox': 19919, 'n_words': 4536}