From 2ddf0c02cbdaf13fb5657a85ae2870e069ea7179 Mon Sep 17 00:00:00 2001
From: Sebastian Kosch <sebastian@aldusleaf.org>
Date: Tue, 7 Jul 2020 15:24:30 -0700
Subject: [PATCH] Add linear-time continuous XYShape combine function

---
 src/components/DistBuilder.re                 |  2 +-
 src/distPlus/distribution/Distributions.re    | 11 ++--
 src/distPlus/distribution/XYShape.re          | 58 ++++++++++++++++++-
 .../expressionTree/ExpressionTreeEvaluator.re |  5 +-
 4 files changed, 63 insertions(+), 13 deletions(-)

diff --git a/src/components/DistBuilder.re b/src/components/DistBuilder.re
index 6ab65759..14872f5f 100644
--- a/src/components/DistBuilder.re
+++ b/src/components/DistBuilder.re
@@ -171,7 +171,7 @@ let make = () => {
       ~onSubmit=({state}) => {None},
       ~initialState={
         //guesstimatorString: "mm(normal(-10, 2), uniform(18, 25), lognormal({mean: 10, stdev: 8}), triangular(31,40,50))",
-        guesstimatorString: "uniform(1,2) * uniform(2, 3)", // , triangular(30, 40, 60)
+        guesstimatorString: "mm(normal(-5,1), normal(0, 1), normal(10, 1), normal(11, 1), normal(16, 1))", // , triangular(30, 40, 60)
         domainType: "Complete",
         xPoint: "50.0",
         xPoint2: "60.0",
diff --git a/src/distPlus/distribution/Distributions.re b/src/distPlus/distribution/Distributions.re
index 963508ab..a081d519 100644
--- a/src/distPlus/distribution/Distributions.re
+++ b/src/distPlus/distribution/Distributions.re
@@ -98,7 +98,7 @@ module Continuous = {
   let combinePointwise =
       (
         ~knownIntegralSumsFn,
-        fn,
+        fn: (float => float => float),
         t1: DistTypes.continuousShape,
         t2: DistTypes.continuousShape,
       )
@@ -112,18 +112,15 @@ module Continuous = {
         t2.knownIntegralSum,
       );
 
-    let res = make(
+    make(
       `Linear,
-      XYShape.PointwiseCombination.combine(
-        ~xsSelection=ALL_XS,
-        ~xToYSelection=XYShape.XtoY.linear,
-        ~fn,
+      XYShape.PointwiseCombination.combineLinear(
+        ~fn=(+.),
         t1.xyShape,
         t2.xyShape,
       ),
       combinedIntegralSum,
     );
-    res
   };
 
   let toLinear = (t: t): option(t) => {
diff --git a/src/distPlus/distribution/XYShape.re b/src/distPlus/distribution/XYShape.re
index 7bea8b06..0fe1ebe0 100644
--- a/src/distPlus/distribution/XYShape.re
+++ b/src/distPlus/distribution/XYShape.re
@@ -175,6 +175,62 @@ module PointwiseCombination = {
     | ALL_XS
     | XS_EVENLY_DIVIDED(int);
 
+  let combineLinear = [%raw {| // : (float => float => float, T.t, T.t) => T.t
+      // This function combines two xyShapes by looping through both of them simultaneously.
+      // It always moves on to the next smallest x, whether that's in the first or second input's xs,
+      // and interpolates the value on the other side, thus accumulating xs and ys.
+      // In this implementation (unlike in XtoY.linear above), values outside of a shape's xs are considered 0.0 (instead of firstY or lastY).
+      // This is written in raw JS because this can still be a bottleneck, and using refs for the i and j indices is quite painful.
+
+      function(fn, t1, t2) {
+        let t1n = t1.xs.length;
+        let t2n = t2.xs.length;
+        let outX = [];
+        let outY = [];
+        let i = -1;
+        let j = -1;
+        while (i <= t1n - 1 && j <= t2n - 1) {
+          let x, ya, yb;
+          if (j == t2n - 1 && i < t1n - 1 ||
+              t1.xs[i+1] < t2.xs[j+1]) { // if a has to catch up to b, or if b is already done
+            i++;
+
+            x = t1.xs[i];
+            ya = t1.ys[i];
+
+            let bFraction = (x - (t2.xs[j] || x)) / ((t2.xs[j+1] || x) - (t2.xs[j] || 0));
+            yb = (t2.ys[j] || 0) * (1-bFraction) + (t2.ys[j+1] || 0) * bFraction;
+          } else if (i == t1n - 1 && j < t2n - 1 ||
+                    t1.xs[i+1] > t2.xs[j+1]) { // if b has to catch up to a, or if a is already done
+            j++;
+
+            x = t2.xs[j];
+            yb = t2.ys[j];
+
+            let aFraction = (x - (t1.xs[i] || x)) / ((t1.xs[i+1] || x) - (t1.xs[i] || 0));
+            ya = (t1.ys[i] || 0) * (1-aFraction) + (t1.ys[i+1] || 0) * aFraction;
+          } else if (i < t1n - 1 && j < t2n && t1.xs[i+1] === t2.xs[j+1]) { // if they happen to be equal, move both ahead
+            i++;
+            j++;
+            x = t1.xs[i];
+            ya = t1.ys[i];
+            yb = t2.ys[j];
+          } else if (i === t1n - 1 && j === t2n - 1) {
+            // finished!
+            i = t1n;
+            j = t2n;
+            continue;
+          } else {
+            console.log("Error!", i, j);
+          }
+
+          outX.push(x);
+          outY.push(fn(ya, yb));
+        }
+        return {xs: outX, ys: outY};
+      }
+    |}];
+
   let combine =
       (
         ~xToYSelection: (float, T.t) => 'a,
@@ -204,7 +260,7 @@ module PointwiseCombination = {
     }
   };
 
-  let combineLinear = combine(~xToYSelection=XtoY.linear);
+  //let combineLinear = combine(~xToYSelection=XtoY.linear);
   let combineStepwise = combine(~xToYSelection=XtoY.stepwiseIncremental);
   let combineIfAtX = combine(~xToYSelection=XtoY.stepwiseIfAtX);
 
diff --git a/src/distPlus/expressionTree/ExpressionTreeEvaluator.re b/src/distPlus/expressionTree/ExpressionTreeEvaluator.re
index d5cd807f..8ef231d7 100644
--- a/src/distPlus/expressionTree/ExpressionTreeEvaluator.re
+++ b/src/distPlus/expressionTree/ExpressionTreeEvaluator.re
@@ -168,11 +168,8 @@ module Normalize = {
   let rec operationToLeaf =
           (toLeaf, renderParams, t: node): result(node, string) => {
     switch (t) {
-    | `RenderedDist(s) => {
-      Js.log2("normed", Distributions.Shape.T.normalize(s));
-      //Js.log2("normed integral", Distributions.Shape.T.normalize(s)));
+    | `RenderedDist(s) =>
         Ok(`RenderedDist(Distributions.Shape.T.normalize(s)));
-        }
     | `SymbolicDist(_) => Ok(t)
     | _ =>
       t