Experimental: continuous/discrete multiplication

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: "mm(normal(-5,1), normal(0, 1), normal(10, 1), normal(11, 1), normal(16, 1))", // , triangular(30, 40, 60)
+        guesstimatorString: "normal(0, 10) * 100", // , triangular(30, 40, 60)
         domainType: "Complete",
         xPoint: "50.0",
         xPoint2: "60.0",

src/distPlus/distribution/AlgebraicShapeCombination.re

@@ -208,3 +208,175 @@ let combineShapesContinuousContinuous =
 
   {xs: outputXs, ys: outputYs};
 };
+
+let toDiscretePointMassesFromDiscrete = (s: DistTypes.xyShape): pointMassesWithMoments => {
+  let n = s |> XYShape.T.length;
+  let {xs, ys}: XYShape.T.t = s;
+  let n = E.A.length(xs);
+
+  let masses: array(float) = Belt.Array.makeUninitializedUnsafe(n); // doesn't include the fake first and last points
+  let means: array(float) = Belt.Array.makeUninitializedUnsafe(n);
+  let variances: array(float) = Belt.Array.makeUninitializedUnsafe(n);
+
+  for (i in 0 to n - 1) {
+    let _ =
+      Belt.Array.set(
+        masses,
+        i,
+        ys[i]
+      );
+
+    let _ =
+      Belt.Array.set(
+        means,
+        i,
+        xs[i]
+      );
+
+    let _ =
+      Belt.Array.set(
+        variances,
+        i,
+        0.0
+      );
+    ();
+  };
+
+  {n, masses, means, variances};
+};
+
+let combineShapesContinuousDiscreteAdd =
+    (op: ExpressionTypes.algebraicOperation, s1: DistTypes.xyShape, s2: DistTypes.xyShape)
+  : DistTypes.xyShape => {
+  let t1n = s1 |> XYShape.T.length;
+  let t2n = s2 |> XYShape.T.length;
+
+  // each x pair is added/subtracted
+  let fn = Operation.Algebraic.toFn(op);
+
+  let outXYShapes: array(array((float, float))) =
+    Belt.Array.makeUninitializedUnsafe(t2n);
+
+  for (j in 0 to t2n - 1) {
+    // for each one of the discrete points
+    // create a new distribution, as long as the original continuous one
+
+    let dxyShape: array((float, float)) =
+      Belt.Array.makeUninitializedUnsafe(t1n);
+    for (i in 0 to t1n - 1) {
+      let _ =
+        Belt.Array.set(
+          dxyShape,
+          i,
+          (fn(s1.xs[i], s2.xs[j]), s1.ys[i] *. s2.ys[j]),
+        );
+      ();
+    };
+
+    let _ = Belt.Array.set(outXYShapes, j, dxyShape);
+    ();
+  };
+
+  outXYShapes
+  |> E.A.fold_left(XYShape.PointwiseCombination.combineLinear((+.)), XYShape.T.empty);
+};
+
+let combineShapesContinuousDiscreteMultiply =
+    (op: ExpressionTypes.algebraicOperation, s1: DistTypes.xyShape, s2: DistTypes.xyShape)
+  : DistTypes.xyShape => {
+  let t1n = s1 |> XYShape.T.length;
+  let t2n = s2 |> XYShape.T.length;
+
+  let t1m = toDiscretePointMassesFromTriangulars(s1);
+  let t2m = toDiscretePointMassesFromDiscrete(s2);
+
+  let combineMeansFn =
+    switch (op) {
+    | `Add => ((m1, m2) => m1 +. m2)
+    | `Subtract => ((m1, m2) => m1 -. m2)
+    | `Multiply => ((m1, m2) => m1 *. m2)
+    | `Divide => ((m1, m2) => m1 /. m2)
+    };
+
+  let combineVariancesFn =
+    switch (op) {
+    | `Add
+    | `Subtract => ((v1, v2, _, _) => v1 +. v2)
+    | `Multiply
+    | `Divide => (
+        (v1, v2, m1, m2) => v1 *. m2 ** 2.
+      )
+    };
+
+  let outputMinX: ref(float) = ref(infinity);
+  let outputMaxX: ref(float) = ref(neg_infinity);
+  let masses: array(float) =
+    Belt.Array.makeUninitializedUnsafe(t1m.n * t2m.n);
+  let means: array(float) =
+    Belt.Array.makeUninitializedUnsafe(t1m.n * t2m.n);
+  let variances: array(float) =
+    Belt.Array.makeUninitializedUnsafe(t1m.n * t2m.n);
+  // then convolve the two sets of pointMassesWithMoments
+  for (i in 0 to t1m.n - 1) {
+    for (j in 0 to t2m.n - 1) {
+      let k = i * t2m.n + j;
+      let _ = Belt.Array.set(masses, k, t1m.masses[i] *. t2m.masses[j]);
+
+      let mean = combineMeansFn(t1m.means[i], t2m.means[j]);
+      let variance =
+        combineVariancesFn(
+          t1m.variances[i],
+          t2m.variances[j],
+          t1m.means[i],
+          t2m.means[j],
+        );
+      let _ = Belt.Array.set(means, k, mean);
+      let _ = Belt.Array.set(variances, k, variance);
+
+      // update bounds
+      let minX = mean -. 2. *. sqrt(variance) *. 1.644854;
+      let maxX = mean +. 2. *. sqrt(variance) *. 1.644854;
+      if (minX < outputMinX^) {
+        outputMinX := minX;
+      };
+      if (maxX > outputMaxX^) {
+        outputMaxX := maxX;
+      };
+    };
+  };
+  
+
+  // we now want to create a set of target points. For now, let's just evenly distribute 200 points between
+  // between the outputMinX and outputMaxX
+  let nOut = 300;
+  let outputXs: array(float) = E.A.Floats.range(outputMinX^, outputMaxX^, nOut);
+  let outputYs: array(float) = Belt.Array.make(nOut, 0.0);
+  // now, for each of the outputYs, accumulate from a Gaussian kernel over each input point.
+  for (j in 0 to E.A.length(masses) - 1) { // go through all of the result points
+    let _ = if (variances[j] > 0. && masses[j] > 0.) {
+      for (i in 0 to E.A.length(outputXs) - 1) { // go through all of the target points
+        let dx = outputXs[i] -. means[j];
+        let contribution = masses[j] *. exp(-. (dx ** 2.) /. (2. *. variances[j])) /. (sqrt(2. *. 3.14159276 *. variances[j]));
+        let _ = Belt.Array.set(outputYs, i, outputYs[i] +. contribution);
+        ();
+      };
+      ();
+    };
+    ();
+  };
+
+  {xs: outputXs, ys: outputYs};
+};
+
+let combineShapesContinuousDiscrete =
+    (op: ExpressionTypes.algebraicOperation, s1: DistTypes.xyShape, s2: DistTypes.xyShape)
+    : DistTypes.xyShape => {
+
+    switch (op) {
+    | `Add
+    | `Subtract => combineShapesContinuousDiscreteAdd(op, s1, s2);
+    | `Multiply
+    | `Divide => combineShapesContinuousDiscreteMultiply(op, s1, s2);
+    };
+
+};

src/distPlus/distribution/Distributions.re

@@ -235,7 +235,7 @@ module Continuous = {
         t
         |> shapeMap(
              XYShape.XsConversion.proportionByProbabilityMass(
-               length,
+          length,
                integral(~cache, t).xyShape,
              ),
            );
@@ -280,60 +280,37 @@ module Continuous = {
      each discrete data point, and then adds them all together. */
   let combineAlgebraicallyWithDiscrete =
       (
-        ~downsample=false,
         op: ExpressionTypes.algebraicOperation,
         t1: t,
         t2: DistTypes.discreteShape,
       ) => {
-    let t1s = t1 |> getShape;
-    let t2s = t2.xyShape; // would like to use Discrete.getShape here, but current file structure doesn't allow for that
-    let t1n = t1s |> XYShape.T.length;
-    let t2n = t2s |> XYShape.T.length;
 
-    let fn = Operation.Algebraic.toFn(op);
+    let s1 = t1 |> getShape;
-
+    let s2 = t2.xyShape;
-    let outXYShapes: array(array((float, float))) =
+    let t1n = s1 |> XYShape.T.length;
-      Belt.Array.makeUninitializedUnsafe(t2n);
+    let t2n = s2 |> XYShape.T.length;
-
+    if (t1n == 0 || t2n == 0) {
-    for (j in 0 to t2n - 1) {
+      empty;
-      // for each one of the discrete points
+    } else {
-      // create a new distribution, as long as the original continuous one
+      let combinedShape =
-
+        AlgebraicShapeCombination.combineShapesContinuousDiscrete(
-      let dxyShape: array((float, float)) =
+          op,
-        Belt.Array.makeUninitializedUnsafe(t1n);
+          s1,
-      for (i in 0 to t1n - 1) {
+          s2,
-        let _ =
+        );
-          Belt.Array.set(
+      let combinedIntegralSum =
-            dxyShape,
+        Common.combineIntegralSums(
-            i,
+          (a, b) => Some(a *. b),
-            (fn(t1s.xs[i], t2s.xs[j]), t1s.ys[i] *. t2s.ys[j]),
+          t1.knownIntegralSum,
-          );
+          t2.knownIntegralSum,
-        ();
+        );
-      };
+      // return a new Continuous distribution
-
+      make(`Linear, combinedShape, combinedIntegralSum);
-      let _ = Belt.Array.set(outXYShapes, j, dxyShape);
-      ();
     };
-
-    let combinedIntegralSum =
-      Common.combineIntegralSums(
-        (a, b) => Some(a *. b),
-        t1.knownIntegralSum,
-        t2.knownIntegralSum,
-      );
-
-    outXYShapes
-    |> E.A.fmap(s => {
-         let xyShape = XYShape.T.fromZippedArray(s);
-         make(`Linear, xyShape, None);
-       })
-    |> reduce((+.))
-    |> updateKnownIntegralSum(combinedIntegralSum);
   };
 
   let combineAlgebraically =
       (
-        ~downsample=false,
         op: ExpressionTypes.algebraicOperation,
         t1: t,
         t2: t,
@@ -475,6 +452,7 @@ module Discrete = {
       type integral = DistTypes.continuousShape;
       let integral = (~cache, t) =>
         if (t |> getShape |> XYShape.T.length > 0) {
+          Js.log2("Integrating discrete shape", XYShape.T.accumulateYs((+.), getShape(t)));
           switch (cache) {
           | Some(c) => c
           | None =>
@@ -849,7 +827,6 @@ module Mixed = {
 
   let combineAlgebraically =
       (
-        ~downsample=false,
         op: ExpressionTypes.algebraicOperation,
         t1: t,
         t2: t,
@@ -861,33 +838,31 @@ module Mixed = {
     // An alternative (to be explored in the future) may be to first perform the full convolution and then to downsample the result;
     // to use non-uniform fast Fourier transforms (for addition only), add web workers or gpu.js, etc. ...
 
-    let downsampleIfTooLarge = (t: t) => {
+    // we have to figure out where to downsample, and how to effectively
-      let sqtl = sqrt(float_of_int(totalLength(t)));
+    //let downsampleIfTooLarge = (t: t) => {
-      sqtl > 10. && downsample ? T.downsample(int_of_float(sqtl), t) : t;
+    //  let sqtl = sqrt(float_of_int(totalLength(t)));
-    };
+    //  sqtl > 10 ? T.downsample(int_of_float(sqtl), t) : t;
+    //};
 
-    let t1d = downsampleIfTooLarge(t1);
+    let t1d = t1; //downsampleIfTooLarge(t1);
-    let t2d = downsampleIfTooLarge(t2);
+    let t2d = t2; //downsampleIfTooLarge(t2);
 
     // continuous (*) continuous => continuous, but also
     // discrete (*) continuous => continuous (and vice versa). We have to take care of all combos and then combine them:
     let ccConvResult =
       Continuous.combineAlgebraically(
-        ~downsample=false,
         op,
         t1d.continuous,
         t2d.continuous,
       );
     let dcConvResult =
       Continuous.combineAlgebraicallyWithDiscrete(
-        ~downsample=false,
         op,
         t2d.continuous,
         t1d.discrete,
       );
     let cdConvResult =
       Continuous.combineAlgebraicallyWithDiscrete(
-        ~downsample=false,
         op,
         t1d.continuous,
         t2d.discrete,
@@ -931,14 +906,13 @@ module Shape = {
     switch (t1, t2) {
     | (Continuous(m1), Continuous(m2)) =>
       DistTypes.Continuous(
-        Continuous.combineAlgebraically(~downsample=true, op, m1, m2),
+        Continuous.combineAlgebraically(op, m1, m2),
       )
     | (Discrete(m1), Discrete(m2)) =>
       DistTypes.Discrete(Discrete.combineAlgebraically(op, m1, m2))
     | (m1, m2) =>
       DistTypes.Mixed(
         Mixed.combineAlgebraically(
-          ~downsample=true,
           op,
           toMixed(m1),
           toMixed(m2),