cleanup: Move integrals to its own module

This commit is contained in:
NunoSempere 2022-09-06 13:59:03 +02:00
parent 8ddd1a166a
commit edce22050a

View File

@ -57,6 +57,187 @@ module Combinatorics = {
} }
} }
module Integration = {
let Helpers = {
let integrateFunctionBetweenWithNumIntegrationPoints = (
aLambda,
min: float,
max: float,
numIntegrationPoints: float, // cast as int?
environment,
reducer,
) => {
let applyFunctionAtFloatToFloatOption = (point: float) => {
// Defined here so that it has access to environment, reducer
let pointAsInternalExpression = FunctionRegistry_Helpers.Wrappers.evNumber(point)
let resultAsInternalExpression = Reducer_Expression_Lambda.doLambdaCall(
aLambda,
list{pointAsInternalExpression},
environment,
reducer,
)
let result = switch resultAsInternalExpression {
| Ok(IEvNumber(x)) => Ok(x)
| Error(_) =>
Error(
"Error 1 in Danger.integrate. It's possible that your function doesn't return a number, try definining auxiliaryFunction(x) = mean(yourFunction(x)) and integrate auxiliaryFunction instead",
)
| _ => Error("Error 2 in Danger.integrate")
}
result
}
// worked example in comments below, assuming min=0, max = 10
let numTotalPoints = Belt.Float.toInt(numIntegrationPoints) // superflous declaration, but useful to keep track that we are interpreting "numIntegrationPoints" as the total number on which we evaluate the function, not e.g., as the inner integration points.
let numInnerPoints = numTotalPoints - 2
let numOuterPoints = 2
let totalWeight = max -. min
let weightForAnInnerPoint = totalWeight /. E.I.toFloat(numTotalPoints - 1)
let weightForAnOuterPoint = totalWeight /. E.I.toFloat(numTotalPoints - 1) /. 2.0
let innerPointIncrement = (max -. min) /. E.I.toFloat(numTotalPoints - 1)
let innerXs = Belt.Array.makeBy(numInnerPoints, i =>
min +. Belt_Float.fromInt(i + 1) *. innerPointIncrement
)
// Gotcha: makeBy goes from 0 to (n-1): <https://rescript-lang.org/docs/manual/latest/api/belt/array#makeby>
let ysOptions = Belt.Array.map(innerXs, x => applyFunctionAtFloatToFloatOption(x))
let okYs = E.A.R.filterOk(ysOptions)
/* Logging, with a worked example. */
// Useful for understanding what is happening.
// assuming min = 0, max = 10, numTotalPoints=10, results below:
let verbose = false
if verbose {
Js.Console.log2("numTotalPoints", numTotalPoints) // 5
Js.Console.log2("numInnerPoints", numInnerPoints) // 3
Js.Console.log2("numOuterPoints", numOuterPoints) // always 2
Js.Console.log2("totalWeight", totalWeight) // 10 - 0 = 10
Js.Console.log2("weightForAnInnerPoint", weightForAnInnerPoint) // 10/4 = 2.5
Js.Console.log2("weightForAnOuterPoint", weightForAnOuterPoint) // 10/4/2 = 1.25
Js.Console.log2(
"weightForAnInnerPoint * numInnerPoints + weightForAnOuterPoint * numOuterPoints",
weightForAnInnerPoint *. E.I.toFloat(numInnerPoints) +.
weightForAnOuterPoint *. E.I.toFloat(numOuterPoints),
) // should be 10
Js.Console.log2(
"sum of weights == totalWeight",
weightForAnInnerPoint *. E.I.toFloat(numInnerPoints) +.
weightForAnOuterPoint *. E.I.toFloat(numOuterPoints) == totalWeight,
) // true
Js.Console.log2("innerPointIncrement", innerPointIncrement) // (10-0)/4 = 2.5
Js.Console.log2("innerXs", innerXs) // 2.5, 5, 7.5
Js.Console.log2("ysOptions", ysOptions)
Js.Console.log2("okYs", okYs)
}
//This is pretty hacky. It should use a result type instead of checking that length matches.
let result = if E.A.length(ysOptions) == E.A.length(okYs) {
let innerPointsSum = okYs->E.A.reduce(0.0, (a, b) => a +. b)
let resultWithOuterPoints = switch (
applyFunctionAtFloatToFloatOption(min),
applyFunctionAtFloatToFloatOption(max),
) {
| (Ok(yMin), Ok(yMax)) => {
let result =
(yMin +. yMax) *. weightForAnOuterPoint +. innerPointsSum *. weightForAnInnerPoint
let wrappedResult = result->ReducerInterface_InternalExpressionValue.IEvNumber->Ok
wrappedResult
}
| (Error(b), _) => Error(b)
| (_, Error(b)) => Error(b)
}
resultWithOuterPoints
} else {
Error(
"Integration error 3 in Danger.integrate. It's possible that your function doesn't return a number, try definining auxiliaryFunction(x) = mean(yourFunction(x)) and integrate auxiliaryFunction instead",
)
}
result
}
}
module Lib = {
// Integral in terms of function, min, max, num points
// Note that execution time will be more predictable, because it
// will only depend on num points and the complexity of the function
let integrateFunctionBetweenWithNumIntegrationPoints = Function.make(
~name="integrateFunctionBetweenWithNumIntegrationPoints",
~nameSpace,
~output=EvtNumber,
~requiresNamespace=false,
~examples=[`Danger.integrateFunctionBetweenWithNumIntegrationPoints({|x| x+1}, 1, 10, 10)`],
// should be [x^2/2 + x]1_10 = (100/2 + 10) - (1/2 + 1) = 60 - 1.5 = 58.5
// https://www.wolframalpha.com/input?i=integrate+x%2B1+from+1+to+10
~definitions=[
FnDefinition.make(
~name="integrateFunctionBetweenWithNumIntegrationPoints",
~inputs=[FRTypeLambda, FRTypeNumber, FRTypeNumber, FRTypeNumber],
~run=(inputs, _, env, reducer) => {
let result = switch inputs {
| [_, _, _, IEvNumber(0.0)] =>
Error("Integration error 4 in Danger.integrate: Increment can't be 0.")
| [IEvLambda(aLambda), IEvNumber(min), IEvNumber(max), IEvNumber(numIntegrationPoints)] =>
Helpers.integrateFunctionBetweenWithNumIntegrationPoints(
aLambda,
min,
max,
numIntegrationPoints,
env,
reducer,
)
| _ =>
Error(
"Integration error 5 in Danger.integrate. Remember that inputs are (function, number (min), number (max), number(increment))",
)
}
result
},
(),
),
],
(),
)
// Integral in terms of function, min, max, epsilon (distance between points)
// Execution time will be less predictable, because it
// will depend on min, max and epsilon together,
// as well and the complexity of the function
let integrateFunctionBetweenWithEpsilon = Function.make(
~name="integrateFunctionBetweenWithEpsilon",
~nameSpace,
~output=EvtNumber,
~requiresNamespace=false,
~examples=[`Danger.integrateFunctionBetweenWithEpsilon({|x| x+1}, 1, 10, 0.1)`],
~definitions=[
FnDefinition.make(
~name="integrateFunctionBetweenWithEpsilon",
~inputs=[FRTypeLambda, FRTypeNumber, FRTypeNumber, FRTypeNumber],
~run=(inputs, _, env, reducer) => {
let result = switch inputs {
| [_, _, _, IEvNumber(0.0)] =>
Error("Integration error in Danger.integrate: Increment can't be 0.")
| [IEvLambda(aLambda), IEvNumber(min), IEvNumber(max), IEvNumber(epsilon)] =>
Helpers.integrateFunctionBetweenWithNumIntegrationPoints(
aLambda,
min,
max,
(max -. min) /. epsilon,
env,
reducer,
)->E.R2.errMap(_ =>
"Integration error 7 in Danger.integrate. Something went wrong along the way"
)
| _ =>
Error(
"Integration error 8 in Danger.integrate. Remember that inputs are (function, number (min), number (max), number(increment))",
)
}
result
},
(),
),
],
(),
)
}
}
module Internals = { module Internals = {
// Probability functions // Probability functions
@ -65,100 +246,7 @@ module Internals = {
->Belt.Array.map(FunctionRegistry_Helpers.Wrappers.evNumber) ->Belt.Array.map(FunctionRegistry_Helpers.Wrappers.evNumber)
->FunctionRegistry_Helpers.Wrappers.evArray ->FunctionRegistry_Helpers.Wrappers.evArray
let integrateFunctionBetweenWithNumIntegrationPoints = (
aLambda,
min: float,
max: float,
numIntegrationPoints: float, // cast as int?
environment,
reducer,
) => {
let applyFunctionAtFloatToFloatOption = (point: float) => {
// Defined here so that it has access to environment, reducer
let pointAsInternalExpression = FunctionRegistry_Helpers.Wrappers.evNumber(point)
let resultAsInternalExpression = Reducer_Expression_Lambda.doLambdaCall(
aLambda,
list{pointAsInternalExpression},
environment,
reducer,
)
let result = switch resultAsInternalExpression {
| Ok(IEvNumber(x)) => Ok(x)
| Error(_) =>
Error(
"Error 1 in Danger.integrate. It's possible that your function doesn't return a number, try definining auxiliaryFunction(x) = mean(yourFunction(x)) and integrate auxiliaryFunction instead",
)
| _ => Error("Error 2 in Danger.integrate")
}
result
}
// worked example in comments below, assuming min=0, max = 10
let numTotalPoints = Belt.Float.toInt(numIntegrationPoints) // superflous declaration, but useful to keep track that we are interpreting "numIntegrationPoints" as the total number on which we evaluate the function, not e.g., as the inner integration points.
let numInnerPoints = numTotalPoints - 2
let numOuterPoints = 2
let totalWeight = max -. min
let weightForAnInnerPoint = totalWeight /. E.I.toFloat(numTotalPoints - 1)
let weightForAnOuterPoint = totalWeight /. E.I.toFloat(numTotalPoints - 1) /. 2.0
let innerPointIncrement = (max -. min) /. E.I.toFloat(numTotalPoints - 1)
let innerXs = Belt.Array.makeBy(numInnerPoints, i =>
min +. Belt_Float.fromInt(i + 1) *. innerPointIncrement
)
// Gotcha: makeBy goes from 0 to (n-1): <https://rescript-lang.org/docs/manual/latest/api/belt/array#makeby>
let ysOptions = Belt.Array.map(innerXs, x => applyFunctionAtFloatToFloatOption(x))
let okYs = E.A.R.filterOk(ysOptions)
/* Logging, with a worked example. */
// Useful for understanding what is happening.
// assuming min = 0, max = 10, numTotalPoints=10, results below:
let verbose = false
if verbose {
Js.Console.log2("numTotalPoints", numTotalPoints) // 5
Js.Console.log2("numInnerPoints", numInnerPoints) // 3
Js.Console.log2("numOuterPoints", numOuterPoints) // always 2
Js.Console.log2("totalWeight", totalWeight) // 10 - 0 = 10
Js.Console.log2("weightForAnInnerPoint", weightForAnInnerPoint) // 10/4 = 2.5
Js.Console.log2("weightForAnOuterPoint", weightForAnOuterPoint) // 10/4/2 = 1.25
Js.Console.log2(
"weightForAnInnerPoint * numInnerPoints + weightForAnOuterPoint * numOuterPoints",
weightForAnInnerPoint *. E.I.toFloat(numInnerPoints) +.
weightForAnOuterPoint *. E.I.toFloat(numOuterPoints),
) // should be 10
Js.Console.log2(
"sum of weights == totalWeight",
weightForAnInnerPoint *. E.I.toFloat(numInnerPoints) +.
weightForAnOuterPoint *. E.I.toFloat(numOuterPoints) == totalWeight,
) // true
Js.Console.log2("innerPointIncrement", innerPointIncrement) // (10-0)/4 = 2.5
Js.Console.log2("innerXs", innerXs) // 2.5, 5, 7.5
Js.Console.log2("ysOptions", ysOptions)
Js.Console.log2("okYs", okYs)
}
//This is pretty hacky. It should use a result type instead of checking that length matches.
let result = if E.A.length(ysOptions) == E.A.length(okYs) {
let innerPointsSum = okYs->E.A.reduce(0.0, (a, b) => a +. b)
let resultWithOuterPoints = switch (
applyFunctionAtFloatToFloatOption(min),
applyFunctionAtFloatToFloatOption(max),
) {
| (Ok(yMin), Ok(yMax)) => {
let result =
(yMin +. yMax) *. weightForAnOuterPoint +. innerPointsSum *. weightForAnInnerPoint
let wrappedResult = result->ReducerInterface_InternalExpressionValue.IEvNumber->Ok
wrappedResult
}
| (Error(b), _) => Error(b)
| (_, Error(b)) => Error(b)
}
resultWithOuterPoints
} else {
Error(
"Integration error 3 in Danger.integrate. It's possible that your function doesn't return a number, try definining auxiliaryFunction(x) = mean(yourFunction(x)) and integrate auxiliaryFunction instead",
)
}
result
}
// Diminishing returns // Diminishing returns
// Helpers // Helpers
type diminishingReturnsAccumulatorInner = { type diminishingReturnsAccumulatorInner = {
@ -272,92 +360,13 @@ module Internals = {
} }
let library = [ let library = [
Combinatorics.Lib.laplace, Combinatorics.Lib.laplace,
Combinatorics.Lib.factorial, Combinatorics.Lib.factorial,
Combinatorics.Lib.choose, Combinatorics.Lib.choose,
Combinatorics.Lib.binomial, Combinatorics.Lib.binomial,
// Integral in terms of function, min, max, num points
// Note that execution time will be more predictable, because it
// will only depend on num points and the complexity of the function
Function.make(
~name="integrateFunctionBetweenWithNumIntegrationPoints",
~nameSpace,
~output=EvtNumber,
~requiresNamespace=false,
~examples=[`Danger.integrateFunctionBetweenWithNumIntegrationPoints({|x| x+1}, 1, 10, 10)`],
// should be [x^2/2 + x]1_10 = (100/2 + 10) - (1/2 + 1) = 60 - 1.5 = 58.5
// https://www.wolframalpha.com/input?i=integrate+x%2B1+from+1+to+10
~definitions=[
FnDefinition.make(
~name="integrateFunctionBetweenWithNumIntegrationPoints",
~inputs=[FRTypeLambda, FRTypeNumber, FRTypeNumber, FRTypeNumber],
~run=(inputs, _, env, reducer) => {
let result = switch inputs {
| [_, _, _, IEvNumber(0.0)] =>
Error("Integration error 4 in Danger.integrate: Increment can't be 0.")
| [IEvLambda(aLambda), IEvNumber(min), IEvNumber(max), IEvNumber(numIntegrationPoints)] =>
Internals.integrateFunctionBetweenWithNumIntegrationPoints(
aLambda,
min,
max,
numIntegrationPoints,
env,
reducer,
)
| _ =>
Error(
"Integration error 5 in Danger.integrate. Remember that inputs are (function, number (min), number (max), number(increment))",
)
}
result
},
(),
),
],
(),
),
// Integral in terms of function, min, max, epsilon (distance between points)
// Execution time will be less predictable, because it
// will depend on min, max and epsilon together,
// as well and the complexity of the function
Function.make(
~name="integrateFunctionBetweenWithEpsilon",
~nameSpace,
~output=EvtNumber,
~requiresNamespace=false,
~examples=[`Danger.integrateFunctionBetweenWithEpsilon({|x| x+1}, 1, 10, 0.1)`],
~definitions=[
FnDefinition.make(
~name="integrateFunctionBetweenWithEpsilon",
~inputs=[FRTypeLambda, FRTypeNumber, FRTypeNumber, FRTypeNumber],
~run=(inputs, _, env, reducer) => {
let result = switch inputs {
| [_, _, _, IEvNumber(0.0)] =>
Error("Integration error in Danger.integrate: Increment can't be 0.")
| [IEvLambda(aLambda), IEvNumber(min), IEvNumber(max), IEvNumber(epsilon)] =>
Internals.integrateFunctionBetweenWithNumIntegrationPoints(
aLambda,
min,
max,
(max -. min) /. epsilon,
env,
reducer,
)->E.R2.errMap(_ =>
"Integration error 7 in Danger.integrate. Something went wrong along the way"
)
| _ =>
Error(
"Integration error 8 in Danger.integrate. Remember that inputs are (function, number (min), number (max), number(increment))",
)
}
result
},
(),
),
],
(),
),
// Diminishing marginal return functions // Diminishing marginal return functions
// There are functions diminishingMarginalReturnsForFunctions2 through diminishingMarginalReturnsForFunctions7 // There are functions diminishingMarginalReturnsForFunctions2 through diminishingMarginalReturnsForFunctions7
// Because of this bug: <https://github.com/quantified-uncertainty/squiggle/issues/1090> // Because of this bug: <https://github.com/quantified-uncertainty/squiggle/issues/1090>