From dca5b21f3ab38eb439f604c82c5b066259f43bbf Mon Sep 17 00:00:00 2001 From: Sebastian Kosch Date: Thu, 16 Jul 2020 18:14:42 -0700 Subject: [PATCH] Move integralCache into shapes; clean up interpolator functions for PointwiseCombinations --- src/components/DistBuilder.re | 5 +- src/components/DistBuilder2.re | 2 +- src/components/Drawer.re | 8 +- src/components/charts/DistPlusPlot.re | 64 ++---- .../distribution/AlgebraicShapeCombination.re | 9 +- src/distPlus/distribution/Continuous.re | 170 +++++++++------ src/distPlus/distribution/Discrete.re | 119 ++++++----- src/distPlus/distribution/DistPlus.re | 26 +-- src/distPlus/distribution/DistPlusTime.re | 4 +- src/distPlus/distribution/DistTypes.re | 22 +- src/distPlus/distribution/Distributions.re | 31 ++- src/distPlus/distribution/Mixed.re | 193 ++++++++++-------- .../distribution/MixedShapeBuilder.re | 8 +- src/distPlus/distribution/Shape.re | 77 ++++--- src/distPlus/distribution/XYShape.re | 114 ++++++----- .../expressionTree/ExpressionTreeEvaluator.re | 16 +- src/distPlus/expressionTree/Operation.re | 8 +- .../renderers/samplesRenderer/Guesstimator.re | 2 +- .../renderers/samplesRenderer/Samples.re | 4 +- src/distPlus/symbolic/SymbolicDist.re | 4 +- src/distPlus/utility/E.re | 2 +- 21 files changed, 513 insertions(+), 375 deletions(-) diff --git a/src/components/DistBuilder.re b/src/components/DistBuilder.re index 556dc728..87746905 100644 --- a/src/components/DistBuilder.re +++ b/src/components/DistBuilder.re @@ -147,7 +147,10 @@ module DemoDist = { ); let response = DistPlusRenderer.run(inputs); switch (RenderTypes.DistPlusRenderer.Outputs.distplus(response)) { - | Some(distPlus) => + | Some(distPlus) => { + let normalizedDistPlus = DistPlus.T.normalize(distPlus); + ; + } | _ => "Correct Guesstimator string input to show a distribution." |> R.ste diff --git a/src/components/DistBuilder2.re b/src/components/DistBuilder2.re index 2627b65b..fd4ef8fc 100644 --- a/src/components/DistBuilder2.re +++ b/src/components/DistBuilder2.re @@ -44,7 +44,7 @@ module DemoDist = { DistPlus.make( ~shape= Continuous( - Continuous.make(`Linear, {xs, ys}, None), + Continuous.make(`Linear, {xs, ys}, None, None), ), ~domain=Complete, ~unit=UnspecifiedDistribution, diff --git a/src/components/Drawer.re b/src/components/Drawer.re index 58ec97b3..ffd4d550 100644 --- a/src/components/Drawer.re +++ b/src/components/Drawer.re @@ -177,7 +177,8 @@ module Convert = { let continuousShape: Types.continuousShape = { xyShape, interpolation: `Linear, - knownIntegralSum: None, + integralSumCache: None, + integralCache: None, }; let integral = XYShape.Analysis.integrateContinuousShape(continuousShape); @@ -189,7 +190,8 @@ module Convert = { ys, }, interpolation: `Linear, - knownIntegralSum: Some(1.0), + integralSumCache: Some(1.0), + integralCache: None, }; continuousShape; }; @@ -673,7 +675,7 @@ module State = { pdf, ); - let cdf = Continuous.T.integral(~cache=None, _pdf); + let cdf = Continuous.T.integral(_pdf); let xs = [||]; let ys = [||]; for (i in 1 to 999) { diff --git a/src/components/charts/DistPlusPlot.re b/src/components/charts/DistPlusPlot.re index 44b7d172..efc5689c 100644 --- a/src/components/charts/DistPlusPlot.re +++ b/src/components/charts/DistPlusPlot.re @@ -51,13 +51,13 @@ let table = (distPlus, x) => { {distPlus - |> DistPlus.T.Integral.xToY(~cache=None, x) + |> DistPlus.T.Integral.xToY(x) |> E.Float.with2DigitsPrecision |> ReasonReact.string} {distPlus - |> DistPlus.T.Integral.sum(~cache=None) + |> DistPlus.T.Integral.sum |> E.Float.with2DigitsPrecision |> ReasonReact.string} @@ -70,15 +70,9 @@ let table = (distPlus, x) => { {"Continuous Total" |> ReasonReact.string} - - {"Scaled Continuous Total" |> ReasonReact.string} - {"Discrete Total" |> ReasonReact.string} - - {"Scaled Discrete Total" |> ReasonReact.string} - @@ -87,17 +81,7 @@ let table = (distPlus, x) => { {distPlus |> DistPlus.T.toContinuous |> E.O.fmap( - Continuous.T.Integral.sum(~cache=None), - ) - |> E.O.fmap(E.Float.with2DigitsPrecision) - |> E.O.default("") - |> ReasonReact.string} - - - {distPlus - |> DistPlus.T.normalizedToContinuous - |> E.O.fmap( - Continuous.T.Integral.sum(~cache=None), + Continuous.T.Integral.sum ) |> E.O.fmap(E.Float.with2DigitsPrecision) |> E.O.default("") @@ -106,15 +90,7 @@ let table = (distPlus, x) => { {distPlus |> DistPlus.T.toDiscrete - |> E.O.fmap(Discrete.T.Integral.sum(~cache=None)) - |> E.O.fmap(E.Float.with2DigitsPrecision) - |> E.O.default("") - |> ReasonReact.string} - - - {distPlus - |> DistPlus.T.normalizedToDiscrete - |> E.O.fmap(Discrete.T.Integral.sum(~cache=None)) + |> E.O.fmap(Discrete.T.Integral.sum) |> E.O.fmap(E.Float.with2DigitsPrecision) |> E.O.default("") |> ReasonReact.string} @@ -143,42 +119,42 @@ let percentiles = distPlus => { {distPlus - |> DistPlus.T.Integral.yToX(~cache=None, 0.01) + |> DistPlus.T.Integral.yToX(0.01) |> showFloat} {distPlus - |> DistPlus.T.Integral.yToX(~cache=None, 0.05) + |> DistPlus.T.Integral.yToX(0.05) |> showFloat} {distPlus - |> DistPlus.T.Integral.yToX(~cache=None, 0.25) + |> DistPlus.T.Integral.yToX(0.25) |> showFloat} {distPlus - |> DistPlus.T.Integral.yToX(~cache=None, 0.5) + |> DistPlus.T.Integral.yToX(0.5) |> showFloat} {distPlus - |> DistPlus.T.Integral.yToX(~cache=None, 0.75) + |> DistPlus.T.Integral.yToX(0.75) |> showFloat} {distPlus - |> DistPlus.T.Integral.yToX(~cache=None, 0.95) + |> DistPlus.T.Integral.yToX(0.95) |> showFloat} {distPlus - |> DistPlus.T.Integral.yToX(~cache=None, 0.99) + |> DistPlus.T.Integral.yToX(0.99) |> showFloat} {distPlus - |> DistPlus.T.Integral.yToX(~cache=None, 0.99999) + |> DistPlus.T.Integral.yToX(0.99999) |> showFloat} @@ -225,10 +201,11 @@ module DistPlusChart = { [@react.component] let make = (~distPlus: DistTypes.distPlus, ~config: chartConfig, ~onHover) => { open DistPlus; - let discrete = distPlus |> T.normalizedToDiscrete |> E.O.fmap(Discrete.getShape); + + let discrete = distPlus |> T.toDiscrete |> E.O.fmap(Discrete.getShape); let continuous = distPlus - |> T.normalizedToContinuous + |> T.toContinuous |> E.O.fmap(Continuous.getShape); let range = T.xTotalRange(distPlus); @@ -236,7 +213,7 @@ module DistPlusChart = { // let minX = // switch ( // distPlus - // |> DistPlus.T.Integral.yToX(~cache=None, 0.0001), + // |> DistPlus.T.Integral.yToX(0.0001), // range, // ) { // | (min, Some(range)) => Some(min -. range *. 0.001) @@ -244,11 +221,11 @@ module DistPlusChart = { // }; let minX = { - distPlus |> DistPlus.T.Integral.yToX(~cache=None, 0.00001); + distPlus |> DistPlus.T.Integral.yToX(0.00001); }; let maxX = { - distPlus |> DistPlus.T.Integral.yToX(~cache=None, 0.99999); + distPlus |> DistPlus.T.Integral.yToX(0.99999); }; let timeScale = distPlus.unit |> DistTypes.DistributionUnit.toJson; @@ -283,11 +260,11 @@ module IntegralChart = { |> Continuous.toLinear |> E.O.fmap(Continuous.getShape); let minX = { - distPlus |> DistPlus.T.Integral.yToX(~cache=None, 0.00001); + distPlus |> DistPlus.T.Integral.yToX(0.00001); }; let maxX = { - distPlus |> DistPlus.T.Integral.yToX(~cache=None, 0.99999); + distPlus |> DistPlus.T.Integral.yToX(0.99999); }; let timeScale = distPlus.unit |> DistTypes.DistributionUnit.toJson; { let (x, setX) = React.useState(() => 0.); let (state, dispatch) = React.useReducer(DistPlusPlotReducer.reducer, DistPlusPlotReducer.init); +
{state.distributions |> E.L.fmapi((index, config) => diff --git a/src/distPlus/distribution/AlgebraicShapeCombination.re b/src/distPlus/distribution/AlgebraicShapeCombination.re index 8301d5e4..1a2a5945 100644 --- a/src/distPlus/distribution/AlgebraicShapeCombination.re +++ b/src/distPlus/distribution/AlgebraicShapeCombination.re @@ -247,7 +247,7 @@ let toDiscretePointMassesFromDiscrete = (s: DistTypes.xyShape): pointMassesWithM let combineShapesContinuousDiscrete = (op: ExpressionTypes.algebraicOperation, s1: DistTypes.xyShape, s2: DistTypes.xyShape) - : DistTypes.xyShape => { + : array(DistTypes.xyShape) => { let t1n = s1 |> XYShape.T.length; let t2n = s2 |> XYShape.T.length; @@ -303,10 +303,5 @@ let combineShapesContinuousDiscrete = }; outXYShapes - |> E.A.fmap(XYShape.T.fromZippedArray) - |> E.A.fold_left( - XYShape.PointwiseCombination.combine((+.), - XYShape.XtoY.linearBetweenPointsExtrapolateZero, - XYShape.XtoY.linearBetweenPointsExtrapolateZero), - XYShape.T.empty); + |> E.A.fmap(XYShape.T.fromZippedArray); }; diff --git a/src/distPlus/distribution/Continuous.re b/src/distPlus/distribution/Continuous.re index b484371c..930cc2af 100644 --- a/src/distPlus/distribution/Continuous.re +++ b/src/distPlus/distribution/Continuous.re @@ -3,30 +3,45 @@ open Distributions; type t = DistTypes.continuousShape; let getShape = (t: t) => t.xyShape; let interpolation = (t: t) => t.interpolation; -let make = (interpolation, xyShape, knownIntegralSum): t => { +let make = (interpolation, xyShape, integralSumCache, integralCache): t => { xyShape, interpolation, - knownIntegralSum, + integralSumCache, + integralCache, }; -let shapeMap = (fn, {xyShape, interpolation, knownIntegralSum}: t): t => { +let shapeMap = (fn, {xyShape, interpolation, integralSumCache, integralCache}: t): t => { xyShape: fn(xyShape), interpolation, - knownIntegralSum, + integralSumCache, + integralCache, }; let lastY = (t: t) => t |> getShape |> XYShape.T.lastY; let oShapeMap = - (fn, {xyShape, interpolation, knownIntegralSum}: t) + (fn, {xyShape, interpolation, integralSumCache, integralCache}: t) : option(DistTypes.continuousShape) => - fn(xyShape) |> E.O.fmap(make(interpolation, _, knownIntegralSum)); + fn(xyShape) |> E.O.fmap(make(interpolation, _, integralSumCache, integralCache)); +let emptyIntegral: DistTypes.continuousShape = { + xyShape: {xs: [|neg_infinity|], ys: [|0.0|]}, + interpolation: `Linear, + integralSumCache: Some(0.0), + integralCache: None, +}; let empty: DistTypes.continuousShape = { xyShape: XYShape.T.empty, interpolation: `Linear, - knownIntegralSum: Some(0.0), + integralSumCache: Some(0.0), + integralCache: Some(emptyIntegral), }; + +let stepwiseToLinear = (t: t): t => + make(`Linear, XYShape.Range.stepwiseToLinear(t.xyShape), t.integralSumCache, t.integralCache); + let combinePointwise = ( - ~knownIntegralSumsFn, + ~integralSumCachesFn=(_, _) => None, + ~integralCachesFn: (t, t) => option(t) =(_, _) => None, + ~extrapolation=`UseZero, fn: (float, float) => float, t1: DistTypes.continuousShape, t2: DistTypes.continuousShape, @@ -36,61 +51,89 @@ let combinePointwise = // can just sum them up. Otherwise, all bets are off. let combinedIntegralSum = Common.combineIntegralSums( - knownIntegralSumsFn, - t1.knownIntegralSum, - t2.knownIntegralSum, + integralSumCachesFn, + t1.integralSumCache, + t2.integralSumCache, ); + // TODO: does it ever make sense to pointwise combine the integrals here? + // It could be done for pointwise additions, but is that ever needed? + + // If combining stepwise and linear, we must convert the stepwise to linear first, + // i.e. add a point at the bottom of each step + let (t1, t2) = switch (t1.interpolation, t2.interpolation) { + | (`Linear, `Linear) => (t1, t2); + | (`Stepwise, `Stepwise) => (t1, t2); + | (`Linear, `Stepwise) => (t1, stepwiseToLinear(t2)); + | (`Stepwise, `Linear) => (stepwiseToLinear(t1), t2); + }; + + let interpolator = XYShape.XtoY.continuousInterpolator(t1.interpolation, extrapolation); + make( `Linear, XYShape.PointwiseCombination.combine( (+.), - XYShape.XtoY.linearBetweenPointsExtrapolateZero, - XYShape.XtoY.linearBetweenPointsExtrapolateZero, + interpolator, + interpolator, t1.xyShape, t2.xyShape, ), combinedIntegralSum, + None, ); }; let toLinear = (t: t): option(t) => { switch (t) { - | {interpolation: `Stepwise, xyShape, knownIntegralSum} => + | {interpolation: `Stepwise, xyShape, integralSumCache, integralCache} => xyShape |> XYShape.Range.stepsToContinuous - |> E.O.fmap(make(`Linear, _, knownIntegralSum)) + |> E.O.fmap(make(`Linear, _, integralSumCache, integralCache)) | {interpolation: `Linear} => Some(t) }; }; let shapeFn = (fn, t: t) => t |> getShape |> fn; -let updateKnownIntegralSum = (knownIntegralSum, t: t): t => { + +let updateIntegralSumCache = (integralSumCache, t: t): t => { ...t, - knownIntegralSum, + integralSumCache, +}; + +let updateIntegralCache = (integralCache, t: t): t => { + ...t, + integralCache, }; let reduce = ( - ~knownIntegralSumsFn: (float, float) => option(float)=(_, _) => None, + ~integralSumCachesFn: (float, float) => option(float)=(_, _) => None, + ~integralCachesFn: (t, t) => option(t)=(_, _) => None, fn, continuousShapes, ) => continuousShapes - |> E.A.fold_left(combinePointwise(~knownIntegralSumsFn, fn), empty); + |> E.A.fold_left(combinePointwise(~integralSumCachesFn, ~integralCachesFn, fn), empty); -let mapY = (~knownIntegralSumFn=_ => None, fn, t: t) => { - let u = E.O.bind(_, knownIntegralSumFn); +let mapY = (~integralSumCacheFn=_ => None, + ~integralCacheFn=_ => None, + fn, t: t) => { let yMapFn = shapeMap(XYShape.T.mapY(fn)); - t |> yMapFn |> updateKnownIntegralSum(u(t.knownIntegralSum)); + t + |> yMapFn + |> updateIntegralSumCache(E.O.bind(t.integralSumCache, integralSumCacheFn)) + |> updateIntegralCache(E.O.bind(t.integralCache, integralCacheFn)); }; -let scaleBy = (~scale=1.0, t: t): t => { +let rec scaleBy = (~scale=1.0, t: t): t => { + let scaledIntegralSumCache = E.O.bind(t.integralSumCache, v => Some(scale *. v)); + let scaledIntegralCache = E.O.bind(t.integralCache, v => Some(scaleBy(~scale, v))); + t |> mapY((r: float) => r *. scale) - |> updateKnownIntegralSum( - E.O.bind(t.knownIntegralSum, v => Some(scale *. v)), - ); + |> updateIntegralSumCache(scaledIntegralSumCache) + |> updateIntegralCache(scaledIntegralCache) }; module T = @@ -135,46 +178,48 @@ module T = let truncatedShape = XYShape.T.fromZippedArray(truncatedZippedPairsWithNewPoints); - make(`Linear, truncatedShape, None); + make(`Linear, truncatedShape, None, None); }; // TODO: This should work with stepwise plots. - let integral = (~cache, t) => - if (t |> getShape |> XYShape.T.length > 0) { - switch (cache) { + let integral = (t) => { + if (t |> getShape |> XYShape.T.isEmpty) { + make(`Linear, {xs: [|neg_infinity|], ys: [|0.0|]}, None, None); + } else { + switch (t.integralCache) { | Some(cache) => cache | None => t |> getShape |> XYShape.Range.integrateWithTriangles |> E.O.toExt("This should not have happened") - |> make(`Linear, _, None) + |> make(`Linear, _, None, None) }; - } else { - make(`Linear, {xs: [|neg_infinity|], ys: [|0.0|]}, None); + }; }; - let downsample = (~cache=None, length, t): t => + let downsample = (length, t): t => t |> shapeMap( XYShape.XsConversion.proportionByProbabilityMass( length, - integral(~cache, t).xyShape, + integral(t).xyShape, ), ); - let integralEndY = (~cache, t: t) => - t.knownIntegralSum |> E.O.default(t |> integral(~cache) |> lastY); - let integralXtoY = (~cache, f, t: t) => - t |> integral(~cache) |> shapeFn(XYShape.XtoY.linear(f)); - let integralYtoX = (~cache, f, t: t) => - t |> integral(~cache) |> shapeFn(XYShape.YtoX.linear(f)); + let integralEndY = (t: t) => + t.integralSumCache |> E.O.default(t |> integral |> lastY); + let integralXtoY = (f, t: t) => + t |> integral |> shapeFn(XYShape.XtoY.linear(f)); + let integralYtoX = (f, t: t) => + t |> integral |> shapeFn(XYShape.YtoX.linear(f)); let toContinuous = t => Some(t); let toDiscrete = _ => None; let normalize = (t: t): t => { t - |> scaleBy(~scale=1. /. integralEndY(~cache=None, t)) - |> updateKnownIntegralSum(Some(1.0)); + |> updateIntegralCache(Some(integral(t))) + |> scaleBy(~scale=1. /. integralEndY(t)) + |> updateIntegralSumCache(Some(1.0)); }; let normalizedToContinuous = t => Some(t |> normalize); @@ -203,7 +248,6 @@ module T = each discrete data point, and then adds them all together. */ let combineAlgebraicallyWithDiscrete = ( - ~downsample=false, op: ExpressionTypes.algebraicOperation, t1: t, t2: DistTypes.discreteShape, @@ -211,27 +255,35 @@ let combineAlgebraicallyWithDiscrete = 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; + if (XYShape.T.isEmpty(t1s) || XYShape.T.isEmpty(t2s)) { + empty; + } else { + let shapeArray = AlgebraicShapeCombination.combineShapesContinuousDiscrete(op, t1s, t2s); - if (t1n > 0 && t2n > 0) { - let combinedShape = AlgebraicShapeCombination.combineShapesContinuousDiscrete(op, t1s, t2s); + let t1Interpolator = XYShape.XtoY.continuousInterpolator(t1.interpolation, `UseZero); + let t2Interpolator = XYShape.XtoY.discreteInterpolator; + + let combinedShape = + shapeArray + |> E.A.fold_left( + XYShape.PointwiseCombination.combine((+.), + t1Interpolator, + t2Interpolator), + XYShape.T.empty); let combinedIntegralSum = Common.combineIntegralSums( (a, b) => Some(a *. b), - t1.knownIntegralSum, - t2.knownIntegralSum, + t1.integralSumCache, + t2.integralSumCache, ); - make(`Linear, combinedShape, combinedIntegralSum); - } else { - empty; - } + // TODO: It could make sense to automatically transform the integrals here (shift or scale) + make(t1.interpolation, combinedShape, combinedIntegralSum, None) + }; }; -let combineAlgebraically = - (~downsample=false, op: ExpressionTypes.algebraicOperation, t1: t, t2: t) => { +let combineAlgebraically = (op: ExpressionTypes.algebraicOperation, t1: t, t2: t) => { let s1 = t1 |> getShape; let s2 = t2 |> getShape; let t1n = s1 |> XYShape.T.length; @@ -244,10 +296,10 @@ let combineAlgebraically = let combinedIntegralSum = Common.combineIntegralSums( (a, b) => Some(a *. b), - t1.knownIntegralSum, - t2.knownIntegralSum, + t1.integralSumCache, + t2.integralSumCache, ); // return a new Continuous distribution - make(`Linear, combinedShape, combinedIntegralSum); + make(`Linear, combinedShape, combinedIntegralSum, None); }; }; diff --git a/src/distPlus/distribution/Discrete.re b/src/distPlus/distribution/Discrete.re index 1ef32ad5..ae204479 100644 --- a/src/distPlus/distribution/Discrete.re +++ b/src/distPlus/distribution/Discrete.re @@ -2,23 +2,25 @@ open Distributions; type t = DistTypes.discreteShape; -let make = (xyShape, knownIntegralSum): t => {xyShape, knownIntegralSum}; -let shapeMap = (fn, {xyShape, knownIntegralSum}: t): t => { +let make = (xyShape, integralSumCache, integralCache): t => {xyShape, integralSumCache, integralCache}; +let shapeMap = (fn, {xyShape, integralSumCache, integralCache}: t): t => { xyShape: fn(xyShape), - knownIntegralSum, + integralSumCache, + integralCache }; let getShape = (t: t) => t.xyShape; -let oShapeMap = (fn, {xyShape, knownIntegralSum}: t): option(t) => - fn(xyShape) |> E.O.fmap(make(_, knownIntegralSum)); +let oShapeMap = (fn, {xyShape, integralSumCache, integralCache}: t): option(t) => + fn(xyShape) |> E.O.fmap(make(_, integralSumCache, integralCache)); -let empty: t = {xyShape: XYShape.T.empty, knownIntegralSum: Some(0.0)}; +let empty: t = {xyShape: XYShape.T.empty, integralSumCache: Some(0.0), integralCache: None}; let shapeFn = (fn, t: t) => t |> getShape |> fn; let lastY = (t: t) => t |> getShape |> XYShape.T.lastY; let combinePointwise = ( - ~knownIntegralSumsFn, + ~integralSumCachesFn = (_, _) => None, + ~integralCachesFn: (DistTypes.continuousShape, DistTypes.continuousShape) => option(DistTypes.continuousShape) = (_, _) => None, fn, t1: DistTypes.discreteShape, t2: DistTypes.discreteShape, @@ -26,38 +28,49 @@ let combinePointwise = : DistTypes.discreteShape => { let combinedIntegralSum = Common.combineIntegralSums( - knownIntegralSumsFn, - t1.knownIntegralSum, - t2.knownIntegralSum, + integralSumCachesFn, + t1.integralSumCache, + t2.integralSumCache, ); + // TODO: does it ever make sense to pointwise combine the integrals here? + // It could be done for pointwise additions, but is that ever needed? + make( XYShape.PointwiseCombination.combine( (+.), - XYShape.XtoY.assumeZeroBetweenPoints, - XYShape.XtoY.assumeZeroBetweenPoints, + XYShape.XtoY.discreteInterpolator, + XYShape.XtoY.discreteInterpolator, t1.xyShape, t2.xyShape, ), combinedIntegralSum, + None, ); }; let reduce = - (~knownIntegralSumsFn=(_, _) => None, fn, discreteShapes) - : DistTypes.discreteShape => + (~integralSumCachesFn=(_, _) => None, + ~integralCachesFn=(_, _) => None, + fn, discreteShapes) + : DistTypes.discreteShape => discreteShapes - |> E.A.fold_left(combinePointwise(~knownIntegralSumsFn, fn), empty); + |> E.A.fold_left(combinePointwise(~integralSumCachesFn, ~integralCachesFn, fn), empty); -let updateKnownIntegralSum = (knownIntegralSum, t: t): t => { +let updateIntegralSumCache = (integralSumCache, t: t): t => { ...t, - knownIntegralSum, + integralSumCache, +}; + +let updateIntegralCache = (integralCache, t: t): t => { + ...t, + integralCache, }; /* This multiples all of the data points together and creates a new discrete distribution from the results. Data points at the same xs get added together. It may be a good idea to downsample t1 and t2 before and/or the result after. */ let combineAlgebraically = - (op: ExpressionTypes.algebraicOperation, t1: t, t2: t) => { + (op: ExpressionTypes.algebraicOperation, t1: t, t2: t): t => { let t1s = t1 |> getShape; let t2s = t2 |> getShape; let t1n = t1s |> XYShape.T.length; @@ -66,8 +79,8 @@ let combineAlgebraically = let combinedIntegralSum = Common.combineIntegralSums( (s1, s2) => Some(s1 *. s2), - t1.knownIntegralSum, - t2.knownIntegralSum, + t1.integralSumCache, + t2.integralSumCache, ); let fn = Operation.Algebraic.toFn(op); @@ -87,31 +100,45 @@ let combineAlgebraically = let combinedShape = XYShape.T.fromZippedArray(rxys); - make(combinedShape, combinedIntegralSum); + make(combinedShape, combinedIntegralSum, None); }; -let mapY = (~knownIntegralSumFn=previousKnownIntegralSum => None, fn, t: t) => { - let u = E.O.bind(_, knownIntegralSumFn); +let mapY = (~integralSumCacheFn=_ => None, + ~integralCacheFn=_ => None, + fn, t: t) => { let yMapFn = shapeMap(XYShape.T.mapY(fn)); - t |> yMapFn |> updateKnownIntegralSum(u(t.knownIntegralSum)); + t + |> yMapFn + |> updateIntegralSumCache(E.O.bind(t.integralSumCache, integralSumCacheFn)) + |> updateIntegralCache(E.O.bind(t.integralCache, integralCacheFn)); }; + let scaleBy = (~scale=1.0, t: t): t => { + let scaledIntegralSumCache = E.O.bind(t.integralSumCache, v => Some(scale *. v)); + let scaledIntegralCache = E.O.bind(t.integralCache, v => Some(Continuous.scaleBy(~scale, v))); + t |> mapY((r: float) => r *. scale) - |> updateKnownIntegralSum( - E.O.bind(t.knownIntegralSum, v => Some(scale *. v)), - ); + |> updateIntegralSumCache(scaledIntegralSumCache) + |> updateIntegralCache(scaledIntegralCache) }; module T = Dist({ type t = DistTypes.discreteShape; type integral = DistTypes.continuousShape; - let integral = (~cache, t) => - if (t |> getShape |> XYShape.T.length > 0) { - switch (cache) { + let integral = (t) => + if (t |> getShape |> XYShape.T.isEmpty) { + Continuous.make( + `Stepwise, + {xs: [|neg_infinity|], ys: [|0.0|]}, + None, + None, + ); + } else { + switch (t.integralCache) { | Some(c) => c | None => { let ts = getShape(t); @@ -123,20 +150,14 @@ module T = |> XYShape.T.concat(prependedZeroPoint) |> XYShape.T.accumulateYs((+.)); - Continuous.make(`Stepwise, integralShape, None); + Continuous.make(`Stepwise, integralShape, None, None); } }; - } else { - Continuous.make( - `Stepwise, - {xs: [|neg_infinity|], ys: [|0.0|]}, - None, - ); }; - let integralEndY = (~cache, t: t) => - t.knownIntegralSum - |> E.O.default(t |> integral(~cache) |> Continuous.lastY); + let integralEndY = (t: t) => + t.integralSumCache + |> E.O.default(t |> integral |> Continuous.lastY); let minX = shapeFn(XYShape.T.minX); let maxX = shapeFn(XYShape.T.maxX); let toDiscreteProbabilityMassFraction = _ => 1.0; @@ -147,14 +168,14 @@ module T = let normalize = (t: t): t => { t - |> scaleBy(~scale=1. /. integralEndY(~cache=None, t)) - |> updateKnownIntegralSum(Some(1.0)); + |> scaleBy(~scale=1. /. integralEndY(t)) + |> updateIntegralSumCache(Some(1.0)); }; let normalizedToContinuous = _ => None; let normalizedToDiscrete = t => Some(t); // TODO: this should be normalized! - let downsample = (~cache=None, i, t: t): t => { + let downsample = (i, t: t): t => { // It's not clear how to downsample a set of discrete points in a meaningful way. // The best we can do is to clip off the smallest values. let currentLength = t |> getShape |> XYShape.T.length; @@ -170,7 +191,7 @@ module T = |> XYShape.Zipped.sortByX |> XYShape.T.fromZippedArray; - make(clippedShape, None); // if someone needs the sum, they'll have to recompute it + make(clippedShape, None, None); // if someone needs the sum, they'll have to recompute it } else { t; }; @@ -188,7 +209,7 @@ module T = ) |> XYShape.T.fromZippedArray; - make(truncatedShape, None); + make(truncatedShape, None, None); }; let xToY = (f, t) => @@ -198,11 +219,11 @@ module T = |> E.O.default(0.0) |> DistTypes.MixedPoint.makeDiscrete; - let integralXtoY = (~cache, f, t) => - t |> integral(~cache) |> Continuous.getShape |> XYShape.XtoY.linear(f); + let integralXtoY = (f, t) => + t |> integral |> Continuous.getShape |> XYShape.XtoY.linear(f); - let integralYtoX = (~cache, f, t) => - t |> integral(~cache) |> Continuous.getShape |> XYShape.YtoX.linear(f); + let integralYtoX = (f, t) => + t |> integral |> Continuous.getShape |> XYShape.YtoX.linear(f); let mean = (t: t): float => { let s = getShape(t); diff --git a/src/distPlus/distribution/DistPlus.re b/src/distPlus/distribution/DistPlus.re index ff3d459b..e0ecbdce 100644 --- a/src/distPlus/distribution/DistPlus.re +++ b/src/distPlus/distribution/DistPlus.re @@ -2,8 +2,7 @@ open DistTypes; type t = DistTypes.distPlus; -let shapeIntegral = shape => - Shape.T.Integral.get(~cache=None, shape); +let shapeIntegral = shape => Shape.T.Integral.get(shape); let make = ( ~shape, @@ -34,6 +33,7 @@ let update = }; let updateShape = (shape, t) => { + Js.log("Updating the shape, recalculating the integral"); let integralCache = shapeIntegral(shape); update(~shape, ~integralCache, t); }; @@ -59,7 +59,6 @@ module T = let normalize = (t: t): t => { let normalizedShape = t |> toShape |> Shape.T.normalize; t |> updateShape(normalizedShape); - // TODO: also adjust for domainIncludedProbabilityMass here. }; let truncate = (leftCutoff, rightCutoff, t: t): t => { @@ -71,6 +70,7 @@ module T = t |> updateShape(truncatedShape); }; + // TODO: is this still needed? let normalizedToContinuous = (t: t) => { t |> toShape @@ -82,6 +82,7 @@ module T = ); }; + // TODO: is this still needed? let normalizedToDiscrete = (t: t) => { t |> toShape @@ -105,34 +106,33 @@ module T = shapeFn(Shape.T.toDiscreteProbabilityMassFraction); // This bit is kind of awkward, could probably use rethinking. - let integral = (~cache, t: t) => + let integral = (t: t) => updateShape(Continuous(t.integralCache), t); - let downsample = (~cache=None, i, t): t => + let downsample = (i, t): t => updateShape(t |> toShape |> Shape.T.downsample(i), t); // todo: adjust for limit, maybe? let mapY = ( - ~knownIntegralSumFn=previousIntegralSum => None, + ~integralSumCacheFn=previousIntegralSum => None, + ~integralCacheFn=previousIntegralCache => None, fn, {shape, _} as t: t, ) : t => - Shape.T.mapY(~knownIntegralSumFn, fn, shape) + Shape.T.mapY(~integralSumCacheFn, fn, shape) |> updateShape(_, t); // get the total of everything - let integralEndY = (~cache as _, t: t) => { + let integralEndY = (t: t) => { Shape.T.Integral.sum( - ~cache=Some(t.integralCache), toShape(t), ); }; // TODO: Fix this below, obviously. Adjust for limits - let integralXtoY = (~cache as _, f, t: t) => { + let integralXtoY = (f, t: t) => { Shape.T.Integral.xToY( - ~cache=Some(t.integralCache), f, toShape(t), ) @@ -140,8 +140,8 @@ module T = }; // TODO: This part is broken when there is a limit, if this is supposed to be taken into account. - let integralYtoX = (~cache as _, f, t: t) => { - Shape.T.Integral.yToX(~cache=None, f, toShape(t)); + let integralYtoX = (f, t: t) => { + Shape.T.Integral.yToX(f, toShape(t)); }; let mean = (t: t) => { diff --git a/src/distPlus/distribution/DistPlusTime.re b/src/distPlus/distribution/DistPlusTime.re index 064ec32c..d5b36f8e 100644 --- a/src/distPlus/distribution/DistPlusTime.re +++ b/src/distPlus/distribution/DistPlusTime.re @@ -23,6 +23,6 @@ include DistPlus.T.Integral; let xToY = (f: TimeTypes.timeInVector, t: t) => { timeInVectorToX(f, t) - |> E.O.fmap(x => DistPlus.T.Integral.xToY(~cache=None, x, t)); + |> E.O.fmap(x => DistPlus.T.Integral.xToY(x, t)); }; - }; \ No newline at end of file + }; diff --git a/src/distPlus/distribution/DistTypes.re b/src/distPlus/distribution/DistTypes.re index 8e1b7b10..10d91667 100644 --- a/src/distPlus/distribution/DistTypes.re +++ b/src/distPlus/distribution/DistTypes.re @@ -14,20 +14,36 @@ type xyShape = { ys: array(float), }; +type interpolation = [ + | `Stepwise + | `Linear +]; +type extrapolation = [ + | `UseZero + | `UseOutermostPoints +]; + +type interpolator = (xyShape, int, float) => float; + type continuousShape = { xyShape, - interpolation: [ | `Stepwise | `Linear], - knownIntegralSum: option(float), + interpolation: interpolation, + integralSumCache: option(float), + integralCache: option(continuousShape), }; type discreteShape = { xyShape, - knownIntegralSum: option(float), + /* interpolation is always `Discrete */ + integralSumCache: option(float), + integralCache: option(continuousShape), }; type mixedShape = { continuous: continuousShape, discrete: discreteShape, + integralSumCache: option(float), + integralCache: option(continuousShape), }; type shapeMonad('a, 'b, 'c) = diff --git a/src/distPlus/distribution/Distributions.re b/src/distPlus/distribution/Distributions.re index 00d366e0..d060ded5 100644 --- a/src/distPlus/distribution/Distributions.re +++ b/src/distPlus/distribution/Distributions.re @@ -4,7 +4,7 @@ module type dist = { let minX: t => float; let maxX: t => float; let mapY: - (~knownIntegralSumFn: float => option(float)=?, float => float, t) => t; + (~integralSumCacheFn: float => option(float)=?, ~integralCacheFn: DistTypes.continuousShape => option(DistTypes.continuousShape)=?, float => float, t) => t; let xToY: (float, t) => DistTypes.mixedPoint; let toShape: t => DistTypes.shape; let toContinuous: t => option(DistTypes.continuousShape); @@ -13,13 +13,13 @@ module type dist = { let normalizedToContinuous: t => option(DistTypes.continuousShape); let normalizedToDiscrete: t => option(DistTypes.discreteShape); let toDiscreteProbabilityMassFraction: t => float; - let downsample: (~cache: option(integral)=?, int, t) => t; + let downsample: (int, t) => t; let truncate: (option(float), option(float), t) => t; - let integral: (~cache: option(integral), t) => integral; - let integralEndY: (~cache: option(integral), t) => float; - let integralXtoY: (~cache: option(integral), float, t) => float; - let integralYtoX: (~cache: option(integral), float, t) => float; + let integral: (t) => integral; + let integralEndY: (t) => float; + let integralXtoY: (float, t) => float; + let integralYtoX: (float, t) => float; let mean: t => float; let variance: t => float; @@ -59,10 +59,23 @@ module Common = { let combineIntegralSums = ( combineFn: (float, float) => option(float), - t1KnownIntegralSum: option(float), - t2KnownIntegralSum: option(float), + t1IntegralSumCache: option(float), + t2IntegralSumCache: option(float), ) => { - switch (t1KnownIntegralSum, t2KnownIntegralSum) { + switch (t1IntegralSumCache, t2IntegralSumCache) { + | (None, _) + | (_, None) => None + | (Some(s1), Some(s2)) => combineFn(s1, s2) + }; + }; + + let combineIntegrals = + ( + combineFn: (DistTypes.continuousShape, DistTypes.continuousShape) => option(DistTypes.continuousShape), + t1IntegralCache: option(DistTypes.continuousShape), + t2IntegralCache: option(DistTypes.continuousShape), + ) => { + switch (t1IntegralCache, t2IntegralCache) { | (None, _) | (_, None) => None | (Some(s1), Some(s2)) => combineFn(s1, s2) diff --git a/src/distPlus/distribution/Mixed.re b/src/distPlus/distribution/Mixed.re index 9e271c39..82d697ee 100644 --- a/src/distPlus/distribution/Mixed.re +++ b/src/distPlus/distribution/Mixed.re @@ -1,7 +1,7 @@ open Distributions; type t = DistTypes.mixedShape; -let make = (~continuous, ~discrete): t => {continuous, discrete}; +let make = (~continuous, ~discrete, integralSumCache, integralCache): t => {continuous, discrete, integralSumCache, integralCache}; let totalLength = (t: t): int => { let continuousLength = @@ -11,31 +11,17 @@ let totalLength = (t: t): int => { continuousLength + discreteLength; }; -let scaleBy = (~scale=1.0, {discrete, continuous}: t): t => { - let scaledDiscrete = Discrete.scaleBy(~scale, discrete); - let scaledContinuous = Continuous.scaleBy(~scale, continuous); - make(~discrete=scaledDiscrete, ~continuous=scaledContinuous); +let scaleBy = (~scale=1.0, t: t): t => { + let scaledDiscrete = Discrete.scaleBy(~scale, t.discrete); + let scaledContinuous = Continuous.scaleBy(~scale, t.continuous); + let scaledIntegralCache = E.O.bind(t.integralCache, v => Some(Continuous.scaleBy(~scale, v))); + let scaledIntegralSumCache = E.O.bind(t.integralSumCache, s => Some(s *. scale)); + make(~discrete=scaledDiscrete, ~continuous=scaledContinuous, scaledIntegralSumCache, scaledIntegralCache); }; let toContinuous = ({continuous}: t) => Some(continuous); let toDiscrete = ({discrete}: t) => Some(discrete); -let combinePointwise = (~knownIntegralSumsFn, fn, t1: t, t2: t) => { - let reducedDiscrete = - [|t1, t2|] - |> E.A.fmap(toDiscrete) - |> E.A.O.concatSomes - |> Discrete.reduce(~knownIntegralSumsFn, fn); - - let reducedContinuous = - [|t1, t2|] - |> E.A.fmap(toContinuous) - |> E.A.O.concatSomes - |> Continuous.reduce(~knownIntegralSumsFn, fn); - - make(~discrete=reducedDiscrete, ~continuous=reducedContinuous); -}; - module T = Dist({ type t = DistTypes.mixedShape; @@ -61,30 +47,35 @@ module T = let truncatedDiscrete = Discrete.T.truncate(leftCutoff, rightCutoff, discrete); - make(~discrete=truncatedDiscrete, ~continuous=truncatedContinuous); + make(~discrete=truncatedDiscrete, ~continuous=truncatedContinuous, None, None); }; let normalize = (t: t): t => { + let continuousIntegral = Continuous.T.Integral.get(t.continuous); + let discreteIntegral = Discrete.T.Integral.get(t.discrete); + + let continuous = t.continuous |> Continuous.updateIntegralCache(Some(continuousIntegral)); + let discrete = t.discrete |> Discrete.updateIntegralCache(Some(discreteIntegral)); let continuousIntegralSum = - Continuous.T.Integral.sum(~cache=None, t.continuous); + Continuous.T.Integral.sum(continuous); let discreteIntegralSum = - Discrete.T.Integral.sum(~cache=None, t.discrete); + Discrete.T.Integral.sum(discrete); let totalIntegralSum = continuousIntegralSum +. discreteIntegralSum; let newContinuousSum = continuousIntegralSum /. totalIntegralSum; let newDiscreteSum = discreteIntegralSum /. totalIntegralSum; let normalizedContinuous = - t.continuous + continuous |> Continuous.scaleBy(~scale=newContinuousSum /. continuousIntegralSum) - |> Continuous.updateKnownIntegralSum(Some(newContinuousSum)); + |> Continuous.updateIntegralSumCache(Some(newContinuousSum)); let normalizedDiscrete = - t.discrete + discrete |> Discrete.scaleBy(~scale=newDiscreteSum /. discreteIntegralSum) - |> Discrete.updateKnownIntegralSum(Some(newDiscreteSum)); + |> Discrete.updateIntegralSumCache(Some(newDiscreteSum)); - make(~continuous=normalizedContinuous, ~discrete=normalizedDiscrete); + make(~continuous=normalizedContinuous, ~discrete=normalizedDiscrete, Some(1.0), None); }; let xToY = (x, t: t) => { @@ -98,23 +89,22 @@ module T = let toDiscreteProbabilityMassFraction = ({discrete, continuous}: t) => { let discreteIntegralSum = - Discrete.T.Integral.sum(~cache=None, discrete); + Discrete.T.Integral.sum(discrete); let continuousIntegralSum = - Continuous.T.Integral.sum(~cache=None, continuous); + Continuous.T.Integral.sum(continuous); let totalIntegralSum = discreteIntegralSum +. continuousIntegralSum; discreteIntegralSum /. totalIntegralSum; }; - let downsample = (~cache=None, count, {discrete, continuous}: t): t => { + let downsample = (count, t: t): t => { // We will need to distribute the new xs fairly between the discrete and continuous shapes. // The easiest way to do this is to simply go by the previous probability masses. - // The cache really isn't helpful here, because we would need two separate caches let discreteIntegralSum = - Discrete.T.Integral.sum(~cache=None, discrete); + Discrete.T.Integral.sum(t.discrete); let continuousIntegralSum = - Continuous.T.Integral.sum(~cache=None, continuous); + Continuous.T.Integral.sum(t.continuous); let totalIntegralSum = discreteIntegralSum +. continuousIntegralSum; // TODO: figure out what to do when the totalIntegralSum is zero. @@ -124,7 +114,7 @@ module T = int_of_float( float_of_int(count) *. (discreteIntegralSum /. totalIntegralSum), ), - discrete, + t.discrete, ); let downsampledContinuous = @@ -132,10 +122,10 @@ module T = int_of_float( float_of_int(count) *. (continuousIntegralSum /. totalIntegralSum), ), - continuous, + t.continuous, ); - {discrete: downsampledDiscrete, continuous: downsampledContinuous}; + {...t, discrete: downsampledDiscrete, continuous: downsampledContinuous}; }; let normalizedToContinuous = (t: t) => Some(normalize(t).continuous); @@ -143,66 +133,69 @@ module T = let normalizedToDiscrete = ({discrete} as t: t) => Some(normalize(t).discrete); - let integral = (~cache, {continuous, discrete}: t) => { - switch (cache) { + let integral = (t: t) => { + switch (t.integralCache) { | Some(cache) => cache | None => // note: if the underlying shapes aren't normalized, then these integrals won't be either -- but that's the way it should be. - let continuousIntegral = - Continuous.T.Integral.get(~cache=None, continuous); - let discreteIntegral = Discrete.T.Integral.get(~cache=None, discrete); + let continuousIntegral = Continuous.T.Integral.get(t.continuous); + let discreteIntegral = Continuous.stepwiseToLinear(Discrete.T.Integral.get(t.discrete)); Continuous.make( `Linear, XYShape.PointwiseCombination.combine( (+.), - XYShape.XtoY.linearBetweenPointsExtrapolateFlat, - XYShape.XtoY.stepwiseBetweenPointsExtrapolateFlat, + XYShape.XtoY.continuousInterpolator(`Linear, `UseOutermostPoints), + XYShape.XtoY.continuousInterpolator(`Linear, `UseOutermostPoints), Continuous.getShape(continuousIntegral), Continuous.getShape(discreteIntegral), ), None, + None, ); }; }; - let integralEndY = (~cache, t: t) => { - integral(~cache, t) |> Continuous.lastY; + let integralEndY = (t: t) => { + t |> integral |> Continuous.lastY; }; - let integralXtoY = (~cache, f, t) => { - t |> integral(~cache) |> Continuous.getShape |> XYShape.XtoY.linear(f); + let integralXtoY = (f, t) => { + t |> integral |> Continuous.getShape |> XYShape.XtoY.linear(f); }; - let integralYtoX = (~cache, f, t) => { - t |> integral(~cache) |> Continuous.getShape |> XYShape.YtoX.linear(f); + let integralYtoX = (f, t) => { + t |> integral |> Continuous.getShape |> XYShape.YtoX.linear(f); }; // This pipes all ys (continuous and discrete) through fn. - // If mapY is a linear operation, we might be able to update the knownIntegralSums as well; + // If mapY is a linear operation, we might be able to update the integralSumCaches as well; // if not, they'll be set to None. let mapY = ( - ~knownIntegralSumFn=previousIntegralSum => None, + ~integralSumCacheFn=previousIntegralSum => None, + ~integralCacheFn=previousIntegral => None, fn, - {discrete, continuous}: t, + t: t, ) : t => { - let u = E.O.bind(_, knownIntegralSumFn); - - let yMappedDiscrete = - discrete + let yMappedDiscrete: DistTypes.discreteShape = + t.discrete |> Discrete.T.mapY(fn) - |> Discrete.updateKnownIntegralSum(u(discrete.knownIntegralSum)); + |> Discrete.updateIntegralSumCache(E.O.bind(t.discrete.integralSumCache, integralSumCacheFn)) + |> Discrete.updateIntegralCache(E.O.bind(t.discrete.integralCache, integralCacheFn)); - let yMappedContinuous = - continuous + let yMappedContinuous: DistTypes.continuousShape = + t.continuous |> Continuous.T.mapY(fn) - |> Continuous.updateKnownIntegralSum(u(continuous.knownIntegralSum)); + |> Continuous.updateIntegralSumCache(E.O.bind(t.continuous.integralSumCache, integralSumCacheFn)) + |> Continuous.updateIntegralCache(E.O.bind(t.continuous.integralCache, integralCacheFn)); { discrete: yMappedDiscrete, - continuous: Continuous.T.mapY(fn, continuous), + continuous: yMappedContinuous, + integralSumCache: E.O.bind(t.integralSumCache, integralSumCacheFn), + integralCache: E.O.bind(t.integralCache, integralCacheFn), }; }; @@ -211,10 +204,8 @@ module T = let continuousMean = Continuous.T.mean(continuous); // the combined mean is the weighted sum of the two: - let discreteIntegralSum = - Discrete.T.Integral.sum(~cache=None, discrete); - let continuousIntegralSum = - Continuous.T.Integral.sum(~cache=None, continuous); + let discreteIntegralSum = Discrete.T.Integral.sum(discrete); + let continuousIntegralSum = Continuous.T.Integral.sum(continuous); let totalIntegralSum = discreteIntegralSum +. continuousIntegralSum; ( @@ -228,10 +219,8 @@ module T = let variance = ({discrete, continuous} as t: t): float => { // the combined mean is the weighted sum of the two: - let discreteIntegralSum = - Discrete.T.Integral.sum(~cache=None, discrete); - let continuousIntegralSum = - Continuous.T.Integral.sum(~cache=None, continuous); + let discreteIntegralSum = Discrete.T.Integral.sum(discrete); + let continuousIntegralSum = Continuous.T.Integral.sum(continuous); let totalIntegralSum = discreteIntegralSum +. continuousIntegralSum; let getMeanOfSquares = ({discrete, continuous}: t) => { @@ -260,7 +249,7 @@ module T = }); let combineAlgebraically = - (~downsample=false, op: ExpressionTypes.algebraicOperation, t1: t, t2: t) + (op: ExpressionTypes.algebraicOperation, t1: t, t2: t) : t => { // Discrete convolution can cause a huge increase in the number of samples, // so we'll first downsample. @@ -268,43 +257,79 @@ let combineAlgebraically = // 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) => { + // TODO: figure out when to downsample strategically. Could be an evaluationParam? + /*let downsampleIfTooLarge = (t: t) => { let sqtl = sqrt(float_of_int(totalLength(t))); sqtl > 10. && downsample ? T.downsample(int_of_float(sqtl), t) : t; }; let t1d = downsampleIfTooLarge(t1); let t2d = 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, + t1.continuous, + t2.continuous, ); let dcConvResult = Continuous.combineAlgebraicallyWithDiscrete( - ~downsample=false, op, - t2d.continuous, - t1d.discrete, + t2.continuous, + t1.discrete, ); let cdConvResult = Continuous.combineAlgebraicallyWithDiscrete( - ~downsample=false, op, - t1d.continuous, - t2d.discrete, + t1.continuous, + t2.discrete, ); let continuousConvResult = Continuous.reduce((+.), [|ccConvResult, dcConvResult, cdConvResult|]); // ... finally, discrete (*) discrete => discrete, obviously: let discreteConvResult = - Discrete.combineAlgebraically(op, t1d.discrete, t2d.discrete); + Discrete.combineAlgebraically(op, t1.discrete, t2.discrete); - {discrete: discreteConvResult, continuous: continuousConvResult}; + let combinedIntegralSum = + Common.combineIntegralSums( + (a, b) => Some(a *. b), + t1.integralSumCache, + t2.integralSumCache, + ); + + {discrete: discreteConvResult, continuous: continuousConvResult, integralSumCache: combinedIntegralSum, integralCache: None}; +}; + +let combinePointwise = (~integralSumCachesFn = (_, _) => None, ~integralCachesFn = (_, _) => None, fn, t1: t, t2: t): t => { + let reducedDiscrete = + [|t1, t2|] + |> E.A.fmap(toDiscrete) + |> E.A.O.concatSomes + |> Discrete.reduce(~integralSumCachesFn, ~integralCachesFn, fn); + + let reducedContinuous = + [|t1, t2|] + |> E.A.fmap(toContinuous) + |> E.A.O.concatSomes + |> Continuous.reduce(~integralSumCachesFn, ~integralCachesFn, fn); + + let combinedIntegralSum = + Common.combineIntegralSums( + integralSumCachesFn, + t1.integralSumCache, + t2.integralSumCache, + ); + + let combinedIntegral = + Common.combineIntegrals( + integralCachesFn, + t1.integralCache, + t2.integralCache, + ); + + make(~discrete=reducedDiscrete, ~continuous=reducedContinuous, combinedIntegralSum, combinedIntegral); }; diff --git a/src/distPlus/distribution/MixedShapeBuilder.re b/src/distPlus/distribution/MixedShapeBuilder.re index 0a930f25..69a725df 100644 --- a/src/distPlus/distribution/MixedShapeBuilder.re +++ b/src/distPlus/distribution/MixedShapeBuilder.re @@ -9,8 +9,8 @@ type assumptions = { }; let buildSimple = (~continuous: option(DistTypes.continuousShape), ~discrete: option(DistTypes.discreteShape)): option(DistTypes.shape) => { - let continuous = continuous |> E.O.default(Continuous.make(`Linear, {xs: [||], ys: [||]}, Some(0.0))); - let discrete = discrete |> E.O.default(Discrete.make({xs: [||], ys: [||]}, Some(0.0))); + let continuous = continuous |> E.O.default(Continuous.make(`Linear, {xs: [||], ys: [||]}, Some(0.0), None)); + let discrete = discrete |> E.O.default(Discrete.make({xs: [||], ys: [||]}, Some(0.0), None)); let cLength = continuous |> Continuous.getShape @@ -25,7 +25,9 @@ let buildSimple = (~continuous: option(DistTypes.continuousShape), ~discrete: op let mixedDist = Mixed.make( ~continuous, - ~discrete + ~discrete, + None, + None, ); Some(Mixed(mixedDist)); }; diff --git a/src/distPlus/distribution/Shape.re b/src/distPlus/distribution/Shape.re index 6b722ec2..140fa29b 100644 --- a/src/distPlus/distribution/Shape.re +++ b/src/distPlus/distribution/Shape.re @@ -15,11 +15,12 @@ let fmap = ((fn1, fn2, fn3), t: t): t => | Continuous(m) => Continuous(fn3(m)) }; + let toMixed = mapToAll(( m => m, - d => Mixed.make(~discrete=d, ~continuous=Continuous.empty), - c => Mixed.make(~discrete=Discrete.empty, ~continuous=c), + d => Mixed.make(~discrete=d, ~continuous=Continuous.empty, d.integralSumCache, d.integralCache), + c => Mixed.make(~discrete=Discrete.empty, ~continuous=c, c.integralSumCache, c.integralCache), )); let combineAlgebraically = @@ -27,19 +28,18 @@ let combineAlgebraically = switch (t1, t2) { | (Continuous(m1), Continuous(m2)) => DistTypes.Continuous( - Continuous.combineAlgebraically(~downsample=true, op, m1, m2), + Continuous.combineAlgebraically(op, m1, m2), ) | (Continuous(m1), Discrete(m2)) | (Discrete(m2), Continuous(m1)) => DistTypes.Continuous( - Continuous.combineAlgebraicallyWithDiscrete(~downsample=false, op, m1, m2), + Continuous.combineAlgebraicallyWithDiscrete(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), @@ -49,20 +49,25 @@ let combineAlgebraically = }; let combinePointwise = - (~knownIntegralSumsFn=(_, _) => None, fn, t1: t, t2: t) => + (~integralSumCachesFn: (float, float) => option(float) = (_, _) => None, + ~integralCachesFn: (DistTypes.continuousShape, DistTypes.continuousShape) => option(DistTypes.continuousShape) = (_, _) => None, + fn, + t1: t, + t2: t) => switch (t1, t2) { | (Continuous(m1), Continuous(m2)) => DistTypes.Continuous( - Continuous.combinePointwise(~knownIntegralSumsFn, fn, m1, m2), + Continuous.combinePointwise(~integralSumCachesFn, ~integralCachesFn, fn, m1, m2), ) | (Discrete(m1), Discrete(m2)) => DistTypes.Discrete( - Discrete.combinePointwise(~knownIntegralSumsFn, fn, m1, m2), + Discrete.combinePointwise(~integralSumCachesFn, ~integralCachesFn, fn, m1, m2), ) | (m1, m2) => DistTypes.Mixed( Mixed.combinePointwise( - ~knownIntegralSumsFn, + ~integralSumCachesFn, + ~integralCachesFn, fn, toMixed(m1), toMixed(m2), @@ -87,7 +92,7 @@ module T = let toContinuous = t => None; let toDiscrete = t => None; - let downsample = (~cache=None, i, t) => + let downsample = (i, t) => fmap( ( Mixed.T.downsample(i), @@ -108,8 +113,14 @@ module T = ); let toDiscreteProbabilityMassFraction = t => 0.0; + let normalize = - fmap((Mixed.T.normalize, Discrete.T.normalize, Continuous.T.normalize)); + fmap(( + Mixed.T.normalize, + Discrete.T.normalize, + Continuous.T.normalize + )); + let toContinuous = mapToAll(( Mixed.T.toContinuous, @@ -143,38 +154,38 @@ module T = Continuous.T.normalizedToContinuous, )); let minX = mapToAll((Mixed.T.minX, Discrete.T.minX, Continuous.T.minX)); - let integral = (~cache) => + let integral = mapToAll(( - Mixed.T.Integral.get(~cache=None), - Discrete.T.Integral.get(~cache=None), - Continuous.T.Integral.get(~cache=None), + Mixed.T.Integral.get, + Discrete.T.Integral.get, + Continuous.T.Integral.get, )); - let integralEndY = (~cache) => + let integralEndY = mapToAll(( - Mixed.T.Integral.sum(~cache=None), - Discrete.T.Integral.sum(~cache), - Continuous.T.Integral.sum(~cache=None), + Mixed.T.Integral.sum, + Discrete.T.Integral.sum, + Continuous.T.Integral.sum, )); - let integralXtoY = (~cache, f) => { + let integralXtoY = (f) => { mapToAll(( - Mixed.T.Integral.xToY(~cache, f), - Discrete.T.Integral.xToY(~cache, f), - Continuous.T.Integral.xToY(~cache, f), + Mixed.T.Integral.xToY(f), + Discrete.T.Integral.xToY(f), + Continuous.T.Integral.xToY(f), )); }; - let integralYtoX = (~cache, f) => { + let integralYtoX = (f) => { mapToAll(( - Mixed.T.Integral.yToX(~cache, f), - Discrete.T.Integral.yToX(~cache, f), - Continuous.T.Integral.yToX(~cache, f), + Mixed.T.Integral.yToX(f), + Discrete.T.Integral.yToX(f), + Continuous.T.Integral.yToX(f), )); }; let maxX = mapToAll((Mixed.T.maxX, Discrete.T.maxX, Continuous.T.maxX)); - let mapY = (~knownIntegralSumFn=previousIntegralSum => None, fn) => + let mapY = (~integralSumCacheFn=previousIntegralSum => None, ~integralCacheFn=previousIntegral=>None, fn) => fmap(( - Mixed.T.mapY(~knownIntegralSumFn, fn), - Discrete.T.mapY(~knownIntegralSumFn, fn), - Continuous.T.mapY(~knownIntegralSumFn, fn), + Mixed.T.mapY(~integralSumCacheFn, ~integralCacheFn, fn), + Discrete.T.mapY(~integralSumCacheFn, ~integralCacheFn, fn), + Continuous.T.mapY(~integralSumCacheFn, ~integralCacheFn, fn), )); let mean = (t: t): float => @@ -197,8 +208,8 @@ let pdf = (f: float, t: t) => { mixedPoint.continuous +. mixedPoint.discrete; }; -let inv = T.Integral.yToX(~cache=None); -let cdf = T.Integral.xToY(~cache=None); +let inv = T.Integral.yToX; +let cdf = T.Integral.xToY; let sample = (t: t): float => { // this can go, already taken care of in Ozzie's sampling branch diff --git a/src/distPlus/distribution/XYShape.re b/src/distPlus/distribution/XYShape.re index fd2e59ca..fcfd92e5 100644 --- a/src/distPlus/distribution/XYShape.re +++ b/src/distPlus/distribution/XYShape.re @@ -19,6 +19,7 @@ module T = { let ys = (t: t) => t.ys; let length = (t: t) => E.A.length(t.xs); let empty = {xs: [||], ys: [||]}; + let isEmpty = (t: t) => length(t) == 0; let minX = (t: t) => t |> xs |> E.A.Sorted.min |> extImp; let maxX = (t: t) => t |> xs |> E.A.Sorted.max |> extImp; let firstY = (t: t) => t |> ys |> E.A.first |> extImp; @@ -143,60 +144,66 @@ module XtoY = { n; }; - /* The extrapolators below can be used e.g. with PointwiseCombination.combine */ - let linearBetweenPointsExtrapolateFlat = (t: T.t, leftIndex: int, x: float) => { - if (leftIndex < 0) { - t.ys[0]; - } else if (leftIndex >= T.length(t) - 1) { - T.lastY(t); - } else { - let x1 = t.xs[leftIndex]; - let x2 = t.xs[leftIndex + 1]; - let y1 = t.ys[leftIndex]; - let y2 = t.ys[leftIndex + 1]; - let fraction = (x -. x1) /. (x2 -. x1); - y1 *. (1. -. fraction) +. y2 *. fraction; + /* Returns a between-points-interpolating function that can be used with PointwiseCombination.combine. + Interpolation can either be stepwise (using the value on the left) or linear. Extrapolation can be `UseZero or `UseOutermostPoints. */ + let continuousInterpolator = (interpolation: DistTypes.interpolation, extrapolation: DistTypes.extrapolation): interpolator => { + switch (interpolation) { + | `Linear => switch (extrapolation) { + | `UseZero => (t: T.t, leftIndex: int, x: float) => { + if (leftIndex < 0) { + 0.0 + } else if (leftIndex >= T.length(t) - 1) { + 0.0 + } else { + let x1 = t.xs[leftIndex]; + let x2 = t.xs[leftIndex + 1]; + let y1 = t.ys[leftIndex]; + let y2 = t.ys[leftIndex + 1]; + let fraction = (x -. x1) /. (x2 -. x1); + y1 *. (1. -. fraction) +. y2 *. fraction; + }; + } + | `UseOutermostPoints => (t: T.t, leftIndex: int, x: float) => { + if (leftIndex < 0) { + t.ys[0]; + } else if (leftIndex >= T.length(t) - 1) { + T.lastY(t); + } else { + let x1 = t.xs[leftIndex]; + let x2 = t.xs[leftIndex + 1]; + let y1 = t.ys[leftIndex]; + let y2 = t.ys[leftIndex + 1]; + let fraction = (x -. x1) /. (x2 -. x1); + y1 *. (1. -. fraction) +. y2 *. fraction; + }; + } + } + | `Stepwise => switch (extrapolation) { + | `UseZero => (t: T.t, leftIndex: int, x: float) => { + if (leftIndex < 0) { + 0.0 + } else if (leftIndex >= T.length(t) - 1) { + 0.0 + } else { + t.ys[leftIndex]; + } + } + | `UseOutermostPoints => (t: T.t, leftIndex: int, x: float) => { + if (leftIndex < 0) { + t.ys[0]; + } else if (leftIndex >= T.length(t) - 1) { + T.lastY(t); + } else { + t.ys[leftIndex]; + } + } + } } }; - let linearBetweenPointsExtrapolateZero = (t: T.t, leftIndex: int, x: float) => { - if (leftIndex < 0) { - 0.0 - } else if (leftIndex >= T.length(t) - 1) { - 0.0 - } else { - let x1 = t.xs[leftIndex]; - let x2 = t.xs[leftIndex + 1]; - let y1 = t.ys[leftIndex]; - let y2 = t.ys[leftIndex + 1]; - let fraction = (x -. x1) /. (x2 -. x1); - y1 *. (1. -. fraction) +. y2 *. fraction; - } - }; - - let stepwiseBetweenPointsExtrapolateFlat = (t: T.t, leftIndex: int, x: float) => { - if (leftIndex < 0) { - t.ys[0]; - } else if (leftIndex >= T.length(t) - 1) { - T.lastY(t); - } else { - t.ys[leftIndex]; - } - }; - - let stepwiseBetweenPointsExtrapolateZero = (t: T.t, leftIndex: int, x: float) => { - if (leftIndex < 0) { - 0.0 - } else if (leftIndex >= T.length(t) - 1) { - 0.0 - } else { - t.ys[leftIndex]; - } - }; - - let assumeZeroBetweenPoints = (t: T.t, leftIndex: int, x: float) => { - 0.0; - }; + /* Returns a between-points-interpolating function that can be used with PointwiseCombination.combine. + For discrete distributions, the probability density between points is zero, so we just return zero here. */ + let discreteInterpolator: interpolator = (t: T.t, leftIndex: int, x: float) => 0.0; }; module XsConversion = { @@ -371,6 +378,11 @@ module Range = { let derivative = mapYsBasedOnRanges(delta_y_over_delta_x); + let stepwiseToLinear = t => { + // adds points at the bottom of each step. + t; + }; + // TODO: It would be nicer if this the diff didn't change the first element, and also maybe if there were a more elegant way of doing this. let stepsToContinuous = t => { let diff = T.xTotalRange(t) |> (r => r *. 0.00001); diff --git a/src/distPlus/expressionTree/ExpressionTreeEvaluator.re b/src/distPlus/expressionTree/ExpressionTreeEvaluator.re index 6c6e6c6f..b9c9c1a5 100644 --- a/src/distPlus/expressionTree/ExpressionTreeEvaluator.re +++ b/src/distPlus/expressionTree/ExpressionTreeEvaluator.re @@ -61,7 +61,8 @@ module VerticalScaling = { (evaluationParams: evaluationParams, scaleOp, t, scaleBy) => { // scaleBy has to be a single float, otherwise we'll return an error. let fn = Operation.Scale.toFn(scaleOp); - let knownIntegralSumFn = Operation.Scale.toKnownIntegralSumFn(scaleOp); + let integralSumCacheFn = Operation.Scale.toIntegralSumCacheFn(scaleOp); + let integralCacheFn = Operation.Scale.toIntegralCacheFn(scaleOp); let renderedShape = render(evaluationParams, t); switch (renderedShape, scaleBy) { @@ -69,7 +70,8 @@ module VerticalScaling = { Ok( `RenderedDist( Shape.T.mapY( - ~knownIntegralSumFn=knownIntegralSumFn(sm), + ~integralSumCacheFn=integralSumCacheFn(sm), + ~integralCacheFn=integralCacheFn(sm), fn(sm), rs, ), @@ -82,13 +84,14 @@ module VerticalScaling = { }; module PointwiseCombination = { - let pointwiseAdd = (evaluationParams: evaluationParams, t1, t2) => { + let pointwiseAdd = (evaluationParams: evaluationParams, t1: t, t2: t) => { switch (render(evaluationParams, t1), render(evaluationParams, t2)) { | (Ok(`RenderedDist(rs1)), Ok(`RenderedDist(rs2))) => Ok( `RenderedDist( Shape.combinePointwise( - ~knownIntegralSumsFn=(a, b) => Some(a +. b), + ~integralSumCachesFn=(a, b) => Some(a +. b), + ~integralCachesFn=(a, b) => Some(Continuous.combinePointwise(~extrapolation=`UseOutermostPoints, (+.), a, b)), (+.), rs1, rs2, @@ -101,7 +104,7 @@ module PointwiseCombination = { }; }; - let pointwiseMultiply = (evaluationParams: evaluationParams, t1, t2) => { + let pointwiseMultiply = (evaluationParams: evaluationParams, t1: t, t2: t) => { // TODO: construct a function that we can easily sample from, to construct // a RenderedDist. Use the xMin and xMax of the rendered shapes to tell the sampling function where to look. Error( @@ -110,7 +113,7 @@ module PointwiseCombination = { }; let operationToLeaf = - (evaluationParams: evaluationParams, pointwiseOp, t1, t2) => { + (evaluationParams: evaluationParams, pointwiseOp: pointwiseOperation, t1: t, t2: t) => { switch (pointwiseOp) { | `Add => pointwiseAdd(evaluationParams, t1, t2) | `Multiply => pointwiseMultiply(evaluationParams, t1, t2) @@ -227,7 +230,6 @@ let toLeaf = node: t, ) : result(t, string) => { - Js.log2("EVALUATION PARAMS", evaluationParams); switch (node) { // Leaf nodes just stay leaf nodes | `SymbolicDist(_) diff --git a/src/distPlus/expressionTree/Operation.re b/src/distPlus/expressionTree/Operation.re index 06194894..26826f3f 100644 --- a/src/distPlus/expressionTree/Operation.re +++ b/src/distPlus/expressionTree/Operation.re @@ -64,11 +64,17 @@ module Scale = { | `Log => {j|verticalLog($value, $scaleBy) |j} }; - let toKnownIntegralSumFn = + let toIntegralSumCacheFn = fun | `Multiply => ((a, b) => Some(a *. b)) | `Exponentiate => ((_, _) => None) | `Log => ((_, _) => None); + + let toIntegralCacheFn = + fun + | `Multiply => ((a, b) => None) // TODO: this could probably just be multiplied out (using Continuous.scaleBy) + | `Exponentiate => ((_, _) => None) + | `Log => ((_, _) => None); }; module T = { diff --git a/src/distPlus/renderers/samplesRenderer/Guesstimator.re b/src/distPlus/renderers/samplesRenderer/Guesstimator.re index a08fb591..6cdc4400 100644 --- a/src/distPlus/renderers/samplesRenderer/Guesstimator.re +++ b/src/distPlus/renderers/samplesRenderer/Guesstimator.re @@ -7,7 +7,7 @@ type discrete = { let jsToDistDiscrete = (d: discrete): DistTypes.discreteShape => {xyShape: { xs: xsGet(d), ys: ysGet(d), -}, knownIntegralSum: None}; +}, integralSumCache: None, integralCache: None}; [@bs.module "./GuesstimatorLibrary.js"] external stringToSamples: (string, int) => array(float) = "stringToSamples"; diff --git a/src/distPlus/renderers/samplesRenderer/Samples.re b/src/distPlus/renderers/samplesRenderer/Samples.re index e96bbdce..f4c024c8 100644 --- a/src/distPlus/renderers/samplesRenderer/Samples.re +++ b/src/distPlus/renderers/samplesRenderer/Samples.re @@ -120,7 +120,7 @@ module T = { |> E.FloatFloatMap.fmap(r => r /. length) |> E.FloatFloatMap.toArray |> XYShape.T.fromZippedArray - |> Discrete.make(_, None); + |> Discrete.make(_, None, None); let pdf = continuousPart |> E.A.length > 5 @@ -150,7 +150,7 @@ module T = { ~outputXYPoints=samplingInputs.outputXYPoints, formatUnitWidth(usedUnitWidth), ) - |> Continuous.make(`Linear, _, None) + |> Continuous.make(`Linear, _, None, None) |> (r => Some((r, foo))); } : None; diff --git a/src/distPlus/symbolic/SymbolicDist.re b/src/distPlus/symbolic/SymbolicDist.re index cdb0f209..0ed65de9 100644 --- a/src/distPlus/symbolic/SymbolicDist.re +++ b/src/distPlus/symbolic/SymbolicDist.re @@ -299,13 +299,13 @@ module T = { switch (d) { | `Float(v) => Discrete( - Discrete.make({xs: [|v|], ys: [|1.0|]}, Some(1.0)), + Discrete.make({xs: [|v|], ys: [|1.0|]}, Some(1.0), None), ) | _ => let xs = interpolateXs(~xSelection=`ByWeight, d, sampleCount); let ys = xs |> E.A.fmap(x => pdf(x, d)); Continuous( - Continuous.make(`Linear, {xs, ys}, Some(1.0)), + Continuous.make(`Linear, {xs, ys}, Some(1.0), None), ); }; }; diff --git a/src/distPlus/utility/E.re b/src/distPlus/utility/E.re index a0e30e79..f5ff541c 100644 --- a/src/distPlus/utility/E.re +++ b/src/distPlus/utility/E.re @@ -464,4 +464,4 @@ module JsArray = { Rationale.Option.toExn("Warning: This should not have happened"), ); let filter = Js.Array.filter; -}; \ No newline at end of file +};