From 104310fce75ea0e740ba51eb52469aaca9b73713 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Mon, 24 Apr 2023 16:40:18 -0600 Subject: [PATCH 01/25] Included math to syntax when missing --- src/function/algebra/derivative.js | 4 ++-- src/function/algebra/leafCount.js | 2 +- src/function/algebra/polynomialRoot.js | 2 +- src/function/algebra/rationalize.js | 8 ++++---- src/function/algebra/resolve.js | 2 +- src/function/algebra/simplify.js | 14 +++++++------- src/function/algebra/simplifyConstant.js | 4 ++-- src/function/algebra/simplifyCore.js | 4 ++-- src/function/algebra/symbolicEqual.js | 16 ++++++++-------- src/function/matrix/sqrtm.js | 2 +- src/type/unit/function/splitUnit.js | 2 +- 11 files changed, 30 insertions(+), 30 deletions(-) diff --git a/src/function/algebra/derivative.js b/src/function/algebra/derivative.js index 7c56a2d8df..503a8cf6d1 100644 --- a/src/function/algebra/derivative.js +++ b/src/function/algebra/derivative.js @@ -43,8 +43,8 @@ export const createDerivative = /* #__PURE__ */ factory(name, dependencies, ({ * * Syntax: * - * derivative(expr, variable) - * derivative(expr, variable, options) + * math.derivative(expr, variable) + * math.derivative(expr, variable, options) * * Examples: * diff --git a/src/function/algebra/leafCount.js b/src/function/algebra/leafCount.js index 0aefa277de..d37cab665a 100644 --- a/src/function/algebra/leafCount.js +++ b/src/function/algebra/leafCount.js @@ -31,7 +31,7 @@ export const createLeafCount = /* #__PURE__ */ factory(name, dependencies, ({ * * Syntax: * - * leafCount(expr) + * math.leafCount(expr) * * Examples: * diff --git a/src/function/algebra/polynomialRoot.js b/src/function/algebra/polynomialRoot.js index 71f0d550bd..2d331f5c18 100644 --- a/src/function/algebra/polynomialRoot.js +++ b/src/function/algebra/polynomialRoot.js @@ -39,7 +39,7 @@ export const createPolynomialRoot = /* #__PURE__ */ factory(name, dependencies, * * Syntax: * - * polynomialRoot(constant, linearCoeff, quadraticCoeff, cubicCoeff) + * math.polynomialRoot(constant, linearCoeff, quadraticCoeff, cubicCoeff) * * Examples: * // linear diff --git a/src/function/algebra/rationalize.js b/src/function/algebra/rationalize.js index 5c560e4411..be7415db37 100644 --- a/src/function/algebra/rationalize.js +++ b/src/function/algebra/rationalize.js @@ -67,10 +67,10 @@ export const createRationalize = /* #__PURE__ */ factory(name, dependencies, ({ * * Syntax: * - * rationalize(expr) - * rationalize(expr, detailed) - * rationalize(expr, scope) - * rationalize(expr, scope, detailed) + * math.rationalize(expr) + * math.rationalize(expr, detailed) + * math.rationalize(expr, scope) + * math.rationalize(expr, scope, detailed) * * Examples: * diff --git a/src/function/algebra/resolve.js b/src/function/algebra/resolve.js index 5e8f80f155..6a1249ed93 100644 --- a/src/function/algebra/resolve.js +++ b/src/function/algebra/resolve.js @@ -25,7 +25,7 @@ export const createResolve = /* #__PURE__ */ factory(name, dependencies, ({ * * Syntax: * - * resolve(expr, scope) + * math.resolve(expr, scope) * * Examples: * diff --git a/src/function/algebra/simplify.js b/src/function/algebra/simplify.js index 9c1f0f983b..82432069d4 100644 --- a/src/function/algebra/simplify.js +++ b/src/function/algebra/simplify.js @@ -153,13 +153,13 @@ export const createSimplify = /* #__PURE__ */ factory(name, dependencies, ( * * Syntax: * - * simplify(expr) - * simplify(expr, rules) - * simplify(expr, rules) - * simplify(expr, rules, scope) - * simplify(expr, rules, scope, options) - * simplify(expr, scope) - * simplify(expr, scope, options) + * math.simplify(expr) + * math.simplify(expr, rules) + * math.simplify(expr, rules) + * math.simplify(expr, rules, scope) + * math.simplify(expr, rules, scope, options) + * math.simplify(expr, scope) + * math.simplify(expr, scope, options) * * Examples: * diff --git a/src/function/algebra/simplifyConstant.js b/src/function/algebra/simplifyConstant.js index 1196eb2e1f..d5f3888b57 100644 --- a/src/function/algebra/simplifyConstant.js +++ b/src/function/algebra/simplifyConstant.js @@ -48,8 +48,8 @@ export const createSimplifyConstant = /* #__PURE__ */ factory(name, dependencies * * Syntax: * - * simplifyConstant(expr) - * simplifyConstant(expr, options) + * math.simplifyConstant(expr) + * math.simplifyConstant(expr, options) * * Examples: * diff --git a/src/function/algebra/simplifyCore.js b/src/function/algebra/simplifyCore.js index 9a7c1b87e5..86ebe850da 100644 --- a/src/function/algebra/simplifyCore.js +++ b/src/function/algebra/simplifyCore.js @@ -82,8 +82,8 @@ export const createSimplifyCore = /* #__PURE__ */ factory(name, dependencies, ({ * * Syntax: * - * simplifyCore(expr) - * simplifyCore(expr, options) + * math.simplifyCore(expr) + * math.simplifyCore(expr, options) * * Examples: * diff --git a/src/function/algebra/symbolicEqual.js b/src/function/algebra/symbolicEqual.js index 4c846bda34..8494ae93ad 100644 --- a/src/function/algebra/symbolicEqual.js +++ b/src/function/algebra/symbolicEqual.js @@ -31,20 +31,20 @@ export const createSymbolicEqual = /* #__PURE__ */ factory(name, dependencies, ( * * Syntax: * - * symbolicEqual(expr1, expr2) - * symbolicEqual(expr1, expr2, options) + * math.symbolicEqual(expr1, expr2) + * math.symbolicEqual(expr1, expr2, options) * * Examples: * - * symbolicEqual('x*y', 'y*x') // Returns true - * symbolicEqual('x*y', 'y*x', {context: {multiply: {commutative: false}}}) // Returns false - * symbolicEqual('x/y', '(y*x^(-1))^(-1)') // Returns true - * symbolicEqual('abs(x)','x') // Returns false - * symbolicEqual('abs(x)','x', simplify.positiveContext) // Returns true + * math.symbolicEqual('x*y', 'y*x') // Returns true + * math.symbolicEqual('x*y', 'y*x', {context: {multiply: {commutative: false}}}) // Returns false + * math.symbolicEqual('x/y', '(y*x^(-1))^(-1)') // Returns true + * math.symbolicEqual('abs(x)','x') // Returns false + * math.symbolicEqual('abs(x)','x', simplify.positiveContext) // Returns true * * See also: * - * simplify, evaluate + * simplify, evaluate * * @param {Node|string} expr1 The first expression to compare * @param {Node|string} expr2 The second expression to compare diff --git a/src/function/matrix/sqrtm.js b/src/function/matrix/sqrtm.js index e59b255d66..ac1b22c925 100644 --- a/src/function/matrix/sqrtm.js +++ b/src/function/matrix/sqrtm.js @@ -49,7 +49,7 @@ export const createSqrtm = /* #__PURE__ */ factory(name, dependencies, ({ typed, * * Syntax: * - * X = math.sqrtm(A) + * math.sqrtm(A) * * Examples: * diff --git a/src/type/unit/function/splitUnit.js b/src/type/unit/function/splitUnit.js index d6a787a5cf..8baae61303 100644 --- a/src/type/unit/function/splitUnit.js +++ b/src/type/unit/function/splitUnit.js @@ -9,7 +9,7 @@ export const createSplitUnit = /* #__PURE__ */ factory(name, dependencies, ({ ty * * Syntax: * - * splitUnit(unit: Unit, parts: Array.) + * math.splitUnit(unit: Unit, parts: Array.) * * Example: * From 431ee46b1f3b4b9914ba7e98c0c4330242cd62b1 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Fri, 26 May 2023 17:29:10 -0600 Subject: [PATCH 02/25] Included solveODE --- src/factoriesAny.js | 1 + src/function/numeric/solveODE.js | 254 +++++++++++++++++++++++++++++++ 2 files changed, 255 insertions(+) create mode 100644 src/function/numeric/solveODE.js diff --git a/src/factoriesAny.js b/src/factoriesAny.js index 5915cb57c1..f349e39462 100644 --- a/src/factoriesAny.js +++ b/src/factoriesAny.js @@ -94,6 +94,7 @@ export { createCtranspose } from './function/matrix/ctranspose.js' export { createZeros } from './function/matrix/zeros.js' export { createFft } from './function/matrix/fft.js' export { createIfft } from './function/matrix/ifft.js' +export { createSolveODE } from './function/numeric/solveODE.js' export { createErf } from './function/special/erf.js' export { createMode } from './function/statistics/mode.js' export { createProd } from './function/statistics/prod.js' diff --git a/src/function/numeric/solveODE.js b/src/function/numeric/solveODE.js new file mode 100644 index 0000000000..aa09813967 --- /dev/null +++ b/src/function/numeric/solveODE.js @@ -0,0 +1,254 @@ +import { isUnit } from '../../utils/is.js' +import { factory } from '../../utils/factory.js' + +const name = 'solveODE' +const dependencies = [ + 'typed', + 'add', + 'subtract', + 'multiply', + 'divide', + 'max', + 'map', + 'abs', + 'isPositive', + 'isNegative', + 'larger', + 'smaller' +] + +export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( + { + typed, + add, + subtract, + multiply, + divide, + max, + map, + abs, + isPositive, + isNegative, + larger, + smaller + } +) => { + /** + * Numerical Integration of Ordinary Differential Equations + * + * Two variable step methods are privded: + * - "RK23": Bogacki–Shampine method + * - "RK45": Dormand-Prince method RK5(4)7M (default) + * + * The arguments are expected as follows. + * + * - `func` should be the forcing function `f(t,y)` + * - `tspan` should be a vector of two numbers or units `[tStart, tEnd]` + * - `y0` the initial values, should be a scalar or a flat array + * - `options` should be an object with the following information: + * - `method` ('RK45'): ['RK23', 'RK45'] + * - `tol` (1e-3): A numeric value + * - `initialStep`: A numeric value or unit + * - `minStep`: minimum step size of the method + * - `maxStep`: maximum step size of the method + * - `minDelta` (0.2): minimum ratio of change for the step + * - `maxDelta` (5): maximum ratio of change for the step + * - `maxIter` (1e4): maximum number of iterations + * + * Syntax: + * + * math.solveODE(func, tspan, y0) + * math.solveODE(func, tspan, y0, options) + * + * Examples: + * + * function func(t, y) {return y} + * const tspan = [0, 4] + * const y0 = 1 + * math.solveODE(func, tspan, y0) + * math.solveODE(func, tspan, [1, 1.1]) + * math.solveODE(func, tspan, y0, { method:"RK23", maxStep:0.1 }) + * + * See also: + * + * simplifyCore, derivative, evaluate, parse, rationalize, resolve + * + * @param {function} func The forcing function f(t,y) + * @param {Array | Matrix} tspan The time span + * @param {number | BigNumber | Unit | Array | Matrix} y0 The initial value + * @param {Object} [options] Optional configuration options + * @return {Object} Return an object with t and y values as arrays + */ + + function _rk (butcherTableau) { + // generates an adaptive runge kutta method from it's butcher tableau + + return function (f, T, y0, options) { + // adaptive runge kutta methods + + const t0 = T[0] // initial time + const tf = T[1] // final time + const steps = 1 // divide the time in this number of steps + const tol = options.tol ? options.tol : 1e-4 // define a tolerance (must be an option) + const maxStep = options.maxStep + const minStep = options.minStep + const minDelta = options.minDelta ? options.minDelta : 0.2 + const maxDelta = options.maxDelta ? options.maxDelta : 5 + const maxIter = options.maxIter ? options.maxIter : 10_000 // stop inifite evaluation if something goes wrong + const [a, c, b, bp] = [butcherTableau.a, butcherTableau.c, butcherTableau.b, butcherTableau.bp] + + let h = options.initialStep + ? options.initialStep + : divide(subtract(tf, t0), steps) // define the initial step size + const t = [t0] // start the time array + const y = [y0] // start the solution array + + const dletaB = subtract(b, bp) // b - bp + + let n = 0 + let iter = 0 + // iterate unitil it reaches either the final time or maximum iterations + while (ongoing(t[n], tf, h) && iter < maxIter) { + const k = [] + + // trim the time step so that it doesn't overshoot + h = trimStep(t[n], tf, h) + + // calculate the first value of k + k.push(f(t[n], y[n])) + + // calculate the rest of the values of k + for (let i = 1; i < c.length; ++i) { + k.push( + f( + add(t[n], multiply(c[i], h)), + add(y[n], multiply(h, a[i], k)) + ) + ) + } + + // estimate the error by comparing solutions of different orders + const TE = max( + abs( + map(multiply(dletaB, k), (X) => + isUnit(X) ? X.value : X + ) + ) + ) + + if (TE < tol && tol / TE > 1 / 4) { + // if within tol then push everything + t.push(add(t[n], h)) + y.push(add(y[n], multiply(h, b, k))) + n++ + } + + // estimate the delta value that will affect the step size + const delta = 0.84 * (tol / TE) ** (1 / 5) + + if (delta < minDelta) { + h = multiply(h, minDelta) + } else if (delta > maxDelta) { + h = multiply(h, maxDelta) + } else { + h = multiply(h, delta) + } + + if (maxStep && larger(abs(h), abs(maxStep))) { + h = isPositive(h) ? abs(maxStep) : multiply(-1, abs(maxStep)) + } else if (minStep && smaller(abs(h), abs(minStep))) { + h = isPositive(h) ? abs(minStep) : multiply(-1, abs(minStep)) + } + iter++ + } + return { t, y } + } + } + + function _rk23 (f, T, y0, options) { + // Bogacki–Shampine method + + // Define the butcher table + const a = [ + [], + [1 / 2], + [0, 3 / 4], + [2 / 9, 1 / 3, 4 / 9] + ] + + const c = [null, 1 / 2, 3 / 4, 1] + const b = [2 / 9, 1 / 3, 4 / 9, 0] + const bp = [7 / 24, 1 / 4, 1 / 3, 1 / 8] + + const butcherTableau = { a, c, b, bp } + + // Solve an adaptive step size rk method + return _rk(butcherTableau)(f, T, y0, options) + } + + function _rk45 (f, T, y0, options) { + // Dormand Prince method + + // Define the butcher tableau + const a = [ + [], + [1 / 5], + [3 / 40, 9 / 40], + [44 / 45, -56 / 15, 32 / 9], + [19372 / 6561, -25360 / 2187, 64448 / 6561, -212 / 729], + [9017 / 3168, -355 / 33, 46732 / 5247, 49 / 176, -5103 / 18656], + [35 / 384, 0, 500 / 1113, 125 / 192, -2187 / 6784, 11 / 84] + ] + + const c = [null, 1 / 5, 3 / 10, 4 / 5, 8 / 9, 1, 1] + const b = [35 / 384, 0, 500 / 1113, 125 / 192, -2187 / 6784, 11 / 84, 0] + const bp = [5179 / 57600, 0, 7571 / 16695, 393 / 640, -92097 / 339200, 187 / 2100, 1 / 40] + + const butcherTableau = { a, c, b, bp } + + // Solve an adaptive step size rk method + return _rk(butcherTableau)(f, T, y0, options) + } + + function _solveODE (f, T, y0, opt) { + const method = opt.method ? opt.method : 'RK45' + const methods = { + RK23: _rk23, + RK45: _rk45 + } + const methodOptions = { ...opt } // clone the options object + delete methodOptions.method // delete the method as it won't be needed + return methods[method.toUpperCase()](f, T, y0, methodOptions) + } + + function ongoing (t, tf, h) { + // returns true if the time has not reached tf for both postitive an negative step (h) + return isPositive(h) + ? smaller(t, tf) + : larger(t, tf) + } + + function trimStep (t, tf, h) { + // Trims the time step so that the next step doesn't overshoot + const next = add(t, h) + return ( + (isPositive(h) && larger(next, tf)) || + (isNegative(h) && smaller(next, tf)) + ) + ? subtract(tf, t) + : h + } + + return typed('solveODE', { + 'function, Array, Array, Object': _solveODE, + 'function, Array, Array': (f, T, y0) => _solveODE(f, T, y0, {}), + 'function, Array, number | BigNumber': (f, T, y0) => { + const sol = _solveODE(f, T, [y0], {}) + return { t: sol.t, y: sol.y.map((Y) => Y[0]) } + }, + 'function, Array, number | BigNumber, Object': (f, T, y0, options) => { + const sol = _solveODE(f, T, [y0], options) + return { t: sol.t, y: sol.y.map((Y) => Y[0]) } + } + }) +}) From 1475b41d014b8764937f2a28fca84b47b7f9ad83 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Mon, 29 May 2023 11:39:48 -0600 Subject: [PATCH 03/25] renamed initialStep as firstStep --- src/function/numeric/solveODE.js | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/function/numeric/solveODE.js b/src/function/numeric/solveODE.js index aa09813967..1a96f350c4 100644 --- a/src/function/numeric/solveODE.js +++ b/src/function/numeric/solveODE.js @@ -36,7 +36,7 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( /** * Numerical Integration of Ordinary Differential Equations * - * Two variable step methods are privded: + * Two variable step methods are provided: * - "RK23": Bogacki–Shampine method * - "RK45": Dormand-Prince method RK5(4)7M (default) * @@ -48,7 +48,7 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( * - `options` should be an object with the following information: * - `method` ('RK45'): ['RK23', 'RK45'] * - `tol` (1e-3): A numeric value - * - `initialStep`: A numeric value or unit + * - `firstStep`: A numeric value or unit * - `minStep`: minimum step size of the method * - `maxStep`: maximum step size of the method * - `minDelta` (0.2): minimum ratio of change for the step @@ -66,7 +66,7 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( * const tspan = [0, 4] * const y0 = 1 * math.solveODE(func, tspan, y0) - * math.solveODE(func, tspan, [1, 1.1]) + * math.solveODE(func, tspan, [1, 2]) * math.solveODE(func, tspan, y0, { method:"RK23", maxStep:0.1 }) * * See also: @@ -97,9 +97,9 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( const maxIter = options.maxIter ? options.maxIter : 10_000 // stop inifite evaluation if something goes wrong const [a, c, b, bp] = [butcherTableau.a, butcherTableau.c, butcherTableau.b, butcherTableau.bp] - let h = options.initialStep - ? options.initialStep - : divide(subtract(tf, t0), steps) // define the initial step size + let h = options.firstStep + ? options.firstStep + : divide(subtract(tf, t0), steps) // define the first step size const t = [t0] // start the time array const y = [y0] // start the solution array From b4cac474e0a651fcc4cf498e756e09aaa51234f6 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Mon, 29 May 2023 11:40:00 -0600 Subject: [PATCH 04/25] Included tests for solveODE --- .../function/numeric/solveODE.test.js | 184 ++++++++++++++++++ 1 file changed, 184 insertions(+) create mode 100644 test/unit-tests/function/numeric/solveODE.test.js diff --git a/test/unit-tests/function/numeric/solveODE.test.js b/test/unit-tests/function/numeric/solveODE.test.js new file mode 100644 index 0000000000..5d8f861b6f --- /dev/null +++ b/test/unit-tests/function/numeric/solveODE.test.js @@ -0,0 +1,184 @@ +import assert from 'assert' +import math from '../../../../src/defaultInstance.js' + +const solveODE = math.solveODE +const matrix = math.matrix +const subtract = math.subtract +const exp = math.exp +const min = math.min +const max = math.max +const abs = math.abs +const tol = 1e-3 +const smaller = math.smaller +const smallerEq = math.smallerEq +const largerEq = math.largerEq +const multiply = math.multiply +const divide = math.divide +const diff = math.diff +const unit = math.unit + +function withinTolerance (A, B, tol) { + const maxAbsError = abs(subtract(A, B)) + return smaller(maxAbsError, tol) +} + +describe('solveODE', function () { + function f (t, y) { return y } + const tspan = [0, 4] + const y0 = [1, 2] + const exactSol = multiply(y0, exp(tspan[1])) // this is only valid for y' = y + + it('should throw an error if the first argument is not a function', function () { + assert.throws(function () { + const sol = solveODE('notAFunction', tspan, y0) + assert.deepStrictEqual( + withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), + [true, true] + ) + }, /TypeError: Unexpected type of argument in function solveODE \(expected: function.*/) + }) + + it('should throw an error if the second argument is not an array like', function () { + assert.throws(function () { + const sol = solveODE(f, 'NotAnArrayLike', y0) + assert.deepStrictEqual( + withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), + [true, true] + ) + }, /TypeError: Unexpected type of argument in function solveODE \(expected: Array or Matrix,.*/) + }) + + it('should throw an error if the second argument does not have a second member', function () { + assert.throws(function () { + const wrongTSpan = [tspan[0]] + const sol = solveODE(f, wrongTSpan, y0) + assert.deepStrictEqual( + withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), + [true, true] + ) + }, /TypeError: Unexpected type of argument in function.*/) + }) + + it('should throw an error if the name of the method is wrong', function () { + assert.throws(function () { + const method = 'wrongMethodName' + const sol = solveODE(f, tspan, y0, { method }) + assert.deepStrictEqual( + withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), + [true, true] + ) + }, /TypeError: .* is not a function/) + }) + + it('should solve close to the analytical solution', function () { + const sol = solveODE(f, tspan, y0) + assert.deepStrictEqual( + withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), + [true, true] + ) + }) + + it('should solve backwards', function () { + const sol = solveODE(f, [tspan[1], tspan[0]], exactSol) + assert.deepStrictEqual( + withinTolerance(sol.y[sol.y.length - 1], y0, tol), + [true, true] + ) + }) + + it('should solve if the arguments are matrices', function () { + const sol = solveODE(f, matrix(tspan), matrix(y0)) + assert.deepStrictEqual( + withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), + [true, true] + ) + }) + + it('should solve if with options even if they are empty', function () { + const options = {} + const sol = solveODE(f, matrix(tspan), matrix(y0), options) + assert.deepStrictEqual( + withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), + [true, true] + ) + }) + + it('should solve when y0 is a scalar', function () { + const sol = solveODE(f, tspan, y0[0]) + assert.deepStrictEqual( + withinTolerance(sol.y[sol.y.length - 1], exactSol[0], tol), + true + ) + }) + + it('should solve with units', function () { + const seconds = unit('s') + function fWithUnits (t, y) { return divide(y, seconds) } + const tspanWithUnits = multiply(tspan, seconds) + const sol = solveODE(fWithUnits, tspanWithUnits, y0) + + assert.deepStrictEqual( + withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), + [true, true] + ) + }) + + it('should solve close to the analytical solution with RK23 method', function () { + const sol = solveODE(f, tspan, y0, { method: 'RK23' }) + assert.deepStrictEqual( + withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), + [true, true] + ) + }) + + it('should solve close to the analytical solution with RK45 method', function () { + const sol = solveODE(f, tspan, y0, { method: 'RK45' }) + assert.deepStrictEqual( + withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), + [true, true] + ) + }) + + it('should solve with method name in lower case', function () { + const sol = solveODE(f, tspan, y0, { method: 'rk45' }) + assert.deepStrictEqual( + withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), + [true, true] + ) + }) + + it('should solve with less steps if a higher tolerance is specified', function () { + const sol = solveODE(f, tspan, y0, { method: 'RK45', tol: 1e-2 }) + assert.deepStrictEqual( + smaller(sol.y.length - 1, 6), + true + ) + }) + + it('should respect maxStep', function () { + const maxStep = 0.4 + const sol = solveODE(f, tspan, y0, { method: 'RK45', maxStep }) + assert.deepStrictEqual( + smallerEq(max(diff(sol.t)), maxStep), + true + ) + }) + + it('should respect minStep', function () { + const minStep = 0.2 + const sol = solveODE(f, tspan, y0, { method: 'RK45', minStep }) + assert.deepStrictEqual( + largerEq(min(diff(sol.t)), minStep), + true + ) + }) + + it('should respect first step', function () { + const firstStep = 0.1 + const sol = solveODE(f, tspan, y0, { firstStep }) + assert.deepStrictEqual( + subtract(sol.t[1], sol.t[0]), + firstStep + ) + }) +}) From ca78866380cc4c8eb6693f3830f5c12365203450 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Mon, 29 May 2023 14:40:29 -0600 Subject: [PATCH 05/25] Test the full state instead of the final state --- .../function/numeric/solveODE.test.js | 63 ++++++++++--------- 1 file changed, 33 insertions(+), 30 deletions(-) diff --git a/test/unit-tests/function/numeric/solveODE.test.js b/test/unit-tests/function/numeric/solveODE.test.js index 5d8f861b6f..b2a7807313 100644 --- a/test/unit-tests/function/numeric/solveODE.test.js +++ b/test/unit-tests/function/numeric/solveODE.test.js @@ -13,27 +13,30 @@ const smaller = math.smaller const smallerEq = math.smallerEq const largerEq = math.largerEq const multiply = math.multiply +const dotMultiply = math.dotMultiply const divide = math.divide const diff = math.diff const unit = math.unit +const map = math.map +const transpose = math.transpose function withinTolerance (A, B, tol) { - const maxAbsError = abs(subtract(A, B)) + const maxAbsError = max(abs(subtract(A, B))) return smaller(maxAbsError, tol) } describe('solveODE', function () { function f (t, y) { return y } + function exactSol (T, y0) { return dotMultiply(y0, map(transpose([T]), exp)) } // this is only valid for y' = y const tspan = [0, 4] const y0 = [1, 2] - const exactSol = multiply(y0, exp(tspan[1])) // this is only valid for y' = y it('should throw an error if the first argument is not a function', function () { assert.throws(function () { const sol = solveODE('notAFunction', tspan, y0) assert.deepStrictEqual( - withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), - [true, true] + withinTolerance(sol.y, exactSol(sol.t, y0), tol), + true ) }, /TypeError: Unexpected type of argument in function solveODE \(expected: function.*/) }) @@ -42,8 +45,8 @@ describe('solveODE', function () { assert.throws(function () { const sol = solveODE(f, 'NotAnArrayLike', y0) assert.deepStrictEqual( - withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), - [true, true] + withinTolerance(sol.y, exactSol(sol.t, y0), tol), + true ) }, /TypeError: Unexpected type of argument in function solveODE \(expected: Array or Matrix,.*/) }) @@ -53,8 +56,8 @@ describe('solveODE', function () { const wrongTSpan = [tspan[0]] const sol = solveODE(f, wrongTSpan, y0) assert.deepStrictEqual( - withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), - [true, true] + withinTolerance(sol.y, exactSol(sol.t, y0), tol), + true ) }, /TypeError: Unexpected type of argument in function.*/) }) @@ -64,8 +67,8 @@ describe('solveODE', function () { const method = 'wrongMethodName' const sol = solveODE(f, tspan, y0, { method }) assert.deepStrictEqual( - withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), - [true, true] + withinTolerance(sol.y, exactSol(sol.t, y0), tol), + true ) }, /TypeError: .* is not a function/) }) @@ -73,24 +76,25 @@ describe('solveODE', function () { it('should solve close to the analytical solution', function () { const sol = solveODE(f, tspan, y0) assert.deepStrictEqual( - withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), - [true, true] + withinTolerance(sol.y, exactSol(sol.t, y0), tol), + true ) }) it('should solve backwards', function () { - const sol = solveODE(f, [tspan[1], tspan[0]], exactSol) + const exactSolEnd = exactSol([4], y0)[0] + const sol = solveODE(f, [tspan[1], tspan[0]], exactSolEnd) assert.deepStrictEqual( - withinTolerance(sol.y[sol.y.length - 1], y0, tol), - [true, true] + withinTolerance(sol.y, exactSol(sol.t, sol.y[sol.y.length - 1]), tol), + true ) }) it('should solve if the arguments are matrices', function () { const sol = solveODE(f, matrix(tspan), matrix(y0)) assert.deepStrictEqual( - withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), - [true, true] + withinTolerance(sol.y, exactSol(sol.t, y0), tol), + true ) }) @@ -98,15 +102,15 @@ describe('solveODE', function () { const options = {} const sol = solveODE(f, matrix(tspan), matrix(y0), options) assert.deepStrictEqual( - withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), - [true, true] + withinTolerance(sol.y, exactSol(sol.t, y0), tol), + true ) }) it('should solve when y0 is a scalar', function () { const sol = solveODE(f, tspan, y0[0]) assert.deepStrictEqual( - withinTolerance(sol.y[sol.y.length - 1], exactSol[0], tol), + withinTolerance(sol.y, exactSol(sol.t, [y0[0]]).map(x => x[0]), 2 * tol), true ) }) @@ -116,41 +120,40 @@ describe('solveODE', function () { function fWithUnits (t, y) { return divide(y, seconds) } const tspanWithUnits = multiply(tspan, seconds) const sol = solveODE(fWithUnits, tspanWithUnits, y0) - assert.deepStrictEqual( - withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), - [true, true] + withinTolerance(sol.y, exactSol(divide(sol.t, seconds), y0), tol), + true ) }) it('should solve close to the analytical solution with RK23 method', function () { const sol = solveODE(f, tspan, y0, { method: 'RK23' }) assert.deepStrictEqual( - withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), - [true, true] + withinTolerance(sol.y, exactSol(sol.t, y0), tol), + true ) }) it('should solve close to the analytical solution with RK45 method', function () { const sol = solveODE(f, tspan, y0, { method: 'RK45' }) assert.deepStrictEqual( - withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), - [true, true] + withinTolerance(sol.y, exactSol(sol.t, y0), tol), + true ) }) it('should solve with method name in lower case', function () { const sol = solveODE(f, tspan, y0, { method: 'rk45' }) assert.deepStrictEqual( - withinTolerance(sol.y[sol.y.length - 1], exactSol, tol), - [true, true] + withinTolerance(sol.y, exactSol(sol.t, y0), tol), + true ) }) it('should solve with less steps if a higher tolerance is specified', function () { const sol = solveODE(f, tspan, y0, { method: 'RK45', tol: 1e-2 }) assert.deepStrictEqual( - smaller(sol.y.length - 1, 6), + smallerEq(sol.y.length, 6), true ) }) From bc0dd3f5dcaddb3acece8ff07a80cd49d7a3b73e Mon Sep 17 00:00:00 2001 From: David Contreras Date: Mon, 29 May 2023 15:21:18 -0600 Subject: [PATCH 06/25] Fixed issue with tolerance --- test/unit-tests/function/numeric/solveODE.test.js | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit-tests/function/numeric/solveODE.test.js b/test/unit-tests/function/numeric/solveODE.test.js index b2a7807313..b850667fc5 100644 --- a/test/unit-tests/function/numeric/solveODE.test.js +++ b/test/unit-tests/function/numeric/solveODE.test.js @@ -110,7 +110,7 @@ describe('solveODE', function () { it('should solve when y0 is a scalar', function () { const sol = solveODE(f, tspan, y0[0]) assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(sol.t, [y0[0]]).map(x => x[0]), 2 * tol), + withinTolerance(sol.y, exactSol(sol.t, [y0[0]]).map(x => x[0]), tol), true ) }) From 35a875e39204f2b9a449a2ad75bf0ef83329ae7f Mon Sep 17 00:00:00 2001 From: David Contreras Date: Tue, 30 May 2023 12:56:02 -0600 Subject: [PATCH 07/25] Added unit signature for y0 --- src/function/numeric/solveODE.js | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/function/numeric/solveODE.js b/src/function/numeric/solveODE.js index 1a96f350c4..da206141da 100644 --- a/src/function/numeric/solveODE.js +++ b/src/function/numeric/solveODE.js @@ -242,11 +242,11 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( return typed('solveODE', { 'function, Array, Array, Object': _solveODE, 'function, Array, Array': (f, T, y0) => _solveODE(f, T, y0, {}), - 'function, Array, number | BigNumber': (f, T, y0) => { + 'function, Array, number | BigNumber | Unit': (f, T, y0) => { const sol = _solveODE(f, T, [y0], {}) return { t: sol.t, y: sol.y.map((Y) => Y[0]) } }, - 'function, Array, number | BigNumber, Object': (f, T, y0, options) => { + 'function, Array, number | BigNumber | Unit, Object': (f, T, y0, options) => { const sol = _solveODE(f, T, [y0], options) return { t: sol.t, y: sol.y.map((Y) => Y[0]) } } From 9eb77ac2f35d6b352399f21b049a2c2408a883d5 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Fri, 2 Jun 2023 14:19:13 -0600 Subject: [PATCH 08/25] Included units test also for y0 --- test/unit-tests/function/numeric/solveODE.test.js | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/unit-tests/function/numeric/solveODE.test.js b/test/unit-tests/function/numeric/solveODE.test.js index b850667fc5..a872723f01 100644 --- a/test/unit-tests/function/numeric/solveODE.test.js +++ b/test/unit-tests/function/numeric/solveODE.test.js @@ -117,11 +117,13 @@ describe('solveODE', function () { it('should solve with units', function () { const seconds = unit('s') + const meters = unit('m') function fWithUnits (t, y) { return divide(y, seconds) } const tspanWithUnits = multiply(tspan, seconds) - const sol = solveODE(fWithUnits, tspanWithUnits, y0) + const y0withUnits = multiply(y0, meters) + const sol = solveODE(fWithUnits, tspanWithUnits, y0withUnits) assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(divide(sol.t, seconds), y0), tol), + withinTolerance(divide(sol.y, meters), exactSol(divide(sol.t, seconds), y0), tol), true ) }) From 727c34e4f5408257d17724d5452bc1a7ea1263e9 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Sat, 1 Jul 2023 17:24:45 -0600 Subject: [PATCH 09/25] Included embedded docs and more tests --- src/expression/embeddedDocs/embeddedDocs.js | 4 + .../embeddedDocs/function/numeric/solveODE.js | 17 +++ src/function/numeric/solveODE.js | 84 +++++++++---- .../function/numeric/solveODE.test.js | 112 +++++++----------- 4 files changed, 127 insertions(+), 90 deletions(-) create mode 100644 src/expression/embeddedDocs/function/numeric/solveODE.js diff --git a/src/expression/embeddedDocs/embeddedDocs.js b/src/expression/embeddedDocs/embeddedDocs.js index e1e74b528d..feebdf2a7b 100644 --- a/src/expression/embeddedDocs/embeddedDocs.js +++ b/src/expression/embeddedDocs/embeddedDocs.js @@ -236,6 +236,7 @@ import { numericDocs } from './function/utils/numeric.js' import { octDocs } from './function/utils/oct.js' import { printDocs } from './function/utils/print.js' import { typeOfDocs } from './function/utils/typeOf.js' +import { solveODEDocs } from './function/numeric/solveODE.js' export const embeddedDocs = { @@ -476,6 +477,9 @@ export const embeddedDocs = { schur: schurDocs, lyap: lyapDocs, + // functions - numeric + solveODE: solveODEDocs, + // functions - probability combinations: combinationsDocs, combinationsWithRep: combinationsWithRepDocs, diff --git a/src/expression/embeddedDocs/function/numeric/solveODE.js b/src/expression/embeddedDocs/function/numeric/solveODE.js new file mode 100644 index 0000000000..f1e566b94e --- /dev/null +++ b/src/expression/embeddedDocs/function/numeric/solveODE.js @@ -0,0 +1,17 @@ +export const solveODEDocs = { + name: 'solveODE', + category: 'Numeric', + syntax: [ + 'solveODE(func, tspan, y0)', + 'solveODE(func, tspan, y0, options)' + ], + description: 'Numerical Integration of Ordinary Differential Equations.', + examples: [ + 'f(t,y) = y', + 'tspan = [0, 4]', + 'solveODE(f, tspan, 1)', + 'solveODE(f, tspan, [1, 2])', + 'solveODE(f, tspan, 1, { method:"RK23", maxStep:0.1 })' + ], + seealso: ['derivative', 'simplifyCore'] +} diff --git a/src/function/numeric/solveODE.js b/src/function/numeric/solveODE.js index da206141da..af0587d757 100644 --- a/src/function/numeric/solveODE.js +++ b/src/function/numeric/solveODE.js @@ -1,4 +1,4 @@ -import { isUnit } from '../../utils/is.js' +import { isUnit, isBigNumber } from '../../utils/is.js' import { factory } from '../../utils/factory.js' const name = 'solveODE' @@ -14,7 +14,9 @@ const dependencies = [ 'isPositive', 'isNegative', 'larger', - 'smaller' + 'smaller', + 'matrix', + 'bignumber' ] export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( @@ -30,7 +32,9 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( isPositive, isNegative, larger, - smaller + smaller, + matrix, + bignumber } ) => { /** @@ -42,9 +46,9 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( * * The arguments are expected as follows. * - * - `func` should be the forcing function `f(t,y)` + * - `func` should be the forcing function `f(t, y)` * - `tspan` should be a vector of two numbers or units `[tStart, tEnd]` - * - `y0` the initial values, should be a scalar or a flat array + * - `y0` the initial state values, should be a scalar or a flat array * - `options` should be an object with the following information: * - `method` ('RK45'): ['RK23', 'RK45'] * - `tol` (1e-3): A numeric value @@ -55,6 +59,12 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( * - `maxDelta` (5): maximum ratio of change for the step * - `maxIter` (1e4): maximum number of iterations * + * The returned value is an object with `{t, y}` please note that even though `t` means time, it can represent any other independant variable like `x`: + * - `t` an array of size `[n]` + * - `y` the states array can be in two ways + * - **if `y0` is a scalar:** returns an array-like of size `[n]` + * - **if `y0` is a flat array-like of size [m]:** returns an array like of size `[n, m]` + * * Syntax: * * math.solveODE(func, tspan, y0) @@ -71,7 +81,7 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( * * See also: * - * simplifyCore, derivative, evaluate, parse, rationalize, resolve + * derivative, simplifyCore * * @param {function} func The forcing function f(t,y) * @param {Array | Matrix} tspan The time span @@ -85,30 +95,37 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( return function (f, T, y0, options) { // adaptive runge kutta methods + const hasBigNumbers = T.some(isBigNumber) || y0.some(isBigNumber) - const t0 = T[0] // initial time - const tf = T[1] // final time - const steps = 1 // divide the time in this number of steps + const t0 = hasBigNumbers ? bignumber(T[0]) : T[0] // initial time + const tf = hasBigNumbers ? bignumber(T[1]) : T[1] // final time + const steps = 1 // divide time in this number of steps const tol = options.tol ? options.tol : 1e-4 // define a tolerance (must be an option) - const maxStep = options.maxStep - const minStep = options.minStep + const maxStep = hasBigNumbers ? bignumber(options.maxStep) : options.maxStep + const minStep = hasBigNumbers ? bignumber(options.minStep) : options.minStep const minDelta = options.minDelta ? options.minDelta : 0.2 const maxDelta = options.maxDelta ? options.maxDelta : 5 const maxIter = options.maxIter ? options.maxIter : 10_000 // stop inifite evaluation if something goes wrong - const [a, c, b, bp] = [butcherTableau.a, butcherTableau.c, butcherTableau.b, butcherTableau.bp] - + const [a, c, b, bp] = hasBigNumbers + ? [ + bignumber(butcherTableau.a), + bignumber(butcherTableau.c), + bignumber(butcherTableau.b), + bignumber(butcherTableau.bp) + ] + : [butcherTableau.a, butcherTableau.c, butcherTableau.b, butcherTableau.bp] let h = options.firstStep - ? options.firstStep + ? hasBigNumbers ? bignumber(options.firstStep) : options.firstStep : divide(subtract(tf, t0), steps) // define the first step size - const t = [t0] // start the time array - const y = [y0] // start the solution array + const t = hasBigNumbers ? [bignumber(t0)] : [t0] // start the time array + const y = hasBigNumbers ? [bignumber(y0)] : [y0] // start the solution array const dletaB = subtract(b, bp) // b - bp let n = 0 let iter = 0 // iterate unitil it reaches either the final time or maximum iterations - while (ongoing(t[n], tf, h) && iter < maxIter) { + while (ongoing(t[n], tf, h)) { const k = [] // trim the time step so that it doesn't overshoot @@ -160,6 +177,9 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( h = isPositive(h) ? abs(minStep) : multiply(-1, abs(minStep)) } iter++ + if (iter > maxIter) { + throw new Error('Maximum number of iterations reached, try changing options') + } } return { t, y } } @@ -216,9 +236,17 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( RK23: _rk23, RK45: _rk45 } - const methodOptions = { ...opt } // clone the options object - delete methodOptions.method // delete the method as it won't be needed - return methods[method.toUpperCase()](f, T, y0, methodOptions) + if (method.toUpperCase() in methods) { + const methodOptions = { ...opt } // clone the options object + delete methodOptions.method // delete the method as it won't be needed + return methods[method.toUpperCase()](f, T, y0, methodOptions) + } else { + // throw an error indicating there is no such method + const methodsWithQuotes = Object.keys(methods).map(x => `"${x}"`) + // generates a string of methods like: "BDF", "RK23" and "RK45" + const availableMethodsString = `${methodsWithQuotes.slice(0, -1).join(', ')} and ${methodsWithQuotes.slice(-1)}` + throw new Error(`Unavailable method "${method}". Available methods are ${availableMethodsString}`) + } } function ongoing (t, tf, h) { @@ -239,16 +267,32 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( : h } + function _matrixSolveODE (f, T, y0, options) { + // receives matrices and returns matrices + const sol = _solveODE(f, T.toArray(), y0.toArray(), options) + return { t: matrix(sol.t), y: matrix(sol.y) } + } + return typed('solveODE', { 'function, Array, Array, Object': _solveODE, + 'function, Matrix, Matrix, Object': _matrixSolveODE, 'function, Array, Array': (f, T, y0) => _solveODE(f, T, y0, {}), + 'function, Matrix, Matrix': (f, T, y0) => _matrixSolveODE(f, T, y0, {}), 'function, Array, number | BigNumber | Unit': (f, T, y0) => { const sol = _solveODE(f, T, [y0], {}) return { t: sol.t, y: sol.y.map((Y) => Y[0]) } }, + 'function, Matrix, number | BigNumber | Unit': (f, T, y0) => { + const sol = _solveODE(f, T.toArray(), [y0], {}) + return { t: matrix(sol.t), y: matrix(sol.y.map((Y) => Y[0])) } + }, 'function, Array, number | BigNumber | Unit, Object': (f, T, y0, options) => { const sol = _solveODE(f, T, [y0], options) return { t: sol.t, y: sol.y.map((Y) => Y[0]) } + }, + 'function, Matrix, number | BigNumber | Unit, Object': (f, T, y0, options) => { + const sol = _solveODE(f, T.toArray(), [y0], options) + return { t: matrix(sol.t), y: matrix(sol.y.map((Y) => Y[0])) } } }) }) diff --git a/test/unit-tests/function/numeric/solveODE.test.js b/test/unit-tests/function/numeric/solveODE.test.js index a872723f01..e8a70dd03a 100644 --- a/test/unit-tests/function/numeric/solveODE.test.js +++ b/test/unit-tests/function/numeric/solveODE.test.js @@ -1,15 +1,15 @@ import assert from 'assert' import math from '../../../../src/defaultInstance.js' +import approx from '../../../../tools/approx.js' const solveODE = math.solveODE const matrix = math.matrix +const add = math.add const subtract = math.subtract const exp = math.exp const min = math.min const max = math.max -const abs = math.abs const tol = 1e-3 -const smaller = math.smaller const smallerEq = math.smallerEq const largerEq = math.largerEq const multiply = math.multiply @@ -19,137 +19,108 @@ const diff = math.diff const unit = math.unit const map = math.map const transpose = math.transpose - -function withinTolerance (A, B, tol) { - const maxAbsError = max(abs(subtract(A, B))) - return smaller(maxAbsError, tol) -} +const isMatrix = math.isMatrix +const bignumber = math.bignumber describe('solveODE', function () { - function f (t, y) { return y } - function exactSol (T, y0) { return dotMultiply(y0, map(transpose([T]), exp)) } // this is only valid for y' = y + function f (t, y) { return subtract(y, t) } + function exactSol (T, y0) { return add(dotMultiply(subtract(y0, 1), map(transpose([T]), exp)), transpose([T]), 1) } // this is only valid for y' = y - t const tspan = [0, 4] const y0 = [1, 2] it('should throw an error if the first argument is not a function', function () { assert.throws(function () { - const sol = solveODE('notAFunction', tspan, y0) - assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(sol.t, y0), tol), - true - ) + solveODE('notAFunction', tspan, y0) }, /TypeError: Unexpected type of argument in function solveODE \(expected: function.*/) }) it('should throw an error if the second argument is not an array like', function () { assert.throws(function () { - const sol = solveODE(f, 'NotAnArrayLike', y0) - assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(sol.t, y0), tol), - true - ) + solveODE(f, 'NotAnArrayLike', y0) }, /TypeError: Unexpected type of argument in function solveODE \(expected: Array or Matrix,.*/) }) it('should throw an error if the second argument does not have a second member', function () { assert.throws(function () { const wrongTSpan = [tspan[0]] - const sol = solveODE(f, wrongTSpan, y0) - assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(sol.t, y0), tol), - true - ) + solveODE(f, wrongTSpan, y0) }, /TypeError: Unexpected type of argument in function.*/) }) it('should throw an error if the name of the method is wrong', function () { assert.throws(function () { const method = 'wrongMethodName' - const sol = solveODE(f, tspan, y0, { method }) - assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(sol.t, y0), tol), - true - ) - }, /TypeError: .* is not a function/) + solveODE(f, tspan, y0, { method }) + }, /Unavailable method.*/) + }) + + it('should throw an error if the maximum number of iterations is reached', function () { + assert.throws(function () { + solveODE(f, tspan, y0, { maxIter: 1 }) + }, /Maximum number of iterations reached.*/) }) it('should solve close to the analytical solution', function () { const sol = solveODE(f, tspan, y0) - assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(sol.t, y0), tol), - true - ) + approx.deepEqual(sol.y, exactSol(sol.t, y0), tol) }) it('should solve backwards', function () { const exactSolEnd = exactSol([4], y0)[0] const sol = solveODE(f, [tspan[1], tspan[0]], exactSolEnd) - assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(sol.t, sol.y[sol.y.length - 1]), tol), - true - ) + approx.deepEqual(sol.y, exactSol(sol.t, sol.y[sol.y.length - 1]), tol) }) it('should solve if the arguments are matrices', function () { const sol = solveODE(f, matrix(tspan), matrix(y0)) + approx.deepEqual(sol.y, matrix(exactSol(sol.t, y0)), tol) assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(sol.t, y0), tol), + isMatrix(sol.t) && isMatrix(sol.y), true ) }) - it('should solve if with options even if they are empty', function () { + /* + it('should solve if the arguments have bignumbers', function () { + const sol = solveODE(f, bignumber(tspan), bignumber(y0)) + approx.deepEqual(sol.y, bignumber(exactSol(sol.t, y0)), tol) + }) + */ + + it('should solve with options even if they are empty', function () { const options = {} - const sol = solveODE(f, matrix(tspan), matrix(y0), options) - assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(sol.t, y0), tol), - true - ) + const sol = solveODE(f, tspan, y0, options) + approx.deepEqual(sol.y, exactSol(sol.t, y0), tol) }) it('should solve when y0 is a scalar', function () { const sol = solveODE(f, tspan, y0[0]) - assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(sol.t, [y0[0]]).map(x => x[0]), tol), - true - ) + approx.deepEqual(sol.y, exactSol(sol.t, [y0[0]]).map(x => x[0]), tol) }) it('should solve with units', function () { const seconds = unit('s') const meters = unit('m') - function fWithUnits (t, y) { return divide(y, seconds) } + function fWithUnits (t, y) { return subtract(divide(y, seconds), multiply(t, divide(meters, multiply(seconds, seconds)))) } const tspanWithUnits = multiply(tspan, seconds) const y0withUnits = multiply(y0, meters) const sol = solveODE(fWithUnits, tspanWithUnits, y0withUnits) - assert.deepStrictEqual( - withinTolerance(divide(sol.y, meters), exactSol(divide(sol.t, seconds), y0), tol), - true - ) + approx.deepEqual(divide(sol.y, meters), exactSol(divide(sol.t, seconds), y0), tol) }) it('should solve close to the analytical solution with RK23 method', function () { const sol = solveODE(f, tspan, y0, { method: 'RK23' }) - assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(sol.t, y0), tol), - true - ) + approx.deepEqual(sol.y, exactSol(sol.t, y0), tol) }) it('should solve close to the analytical solution with RK45 method', function () { const sol = solveODE(f, tspan, y0, { method: 'RK45' }) - assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(sol.t, y0), tol), - true - ) + approx.deepEqual(sol.y, exactSol(sol.t, y0), tol) }) it('should solve with method name in lower case', function () { const sol = solveODE(f, tspan, y0, { method: 'rk45' }) - assert.deepStrictEqual( - withinTolerance(sol.y, exactSol(sol.t, y0), tol), - true - ) + approx.deepEqual(sol.y, exactSol(sol.t, y0), tol) }) it('should solve with less steps if a higher tolerance is specified', function () { @@ -162,7 +133,7 @@ describe('solveODE', function () { it('should respect maxStep', function () { const maxStep = 0.4 - const sol = solveODE(f, tspan, y0, { method: 'RK45', maxStep }) + const sol = solveODE(f, tspan, y0, { maxStep }) assert.deepStrictEqual( smallerEq(max(diff(sol.t)), maxStep), true @@ -170,10 +141,11 @@ describe('solveODE', function () { }) it('should respect minStep', function () { - const minStep = 0.2 - const sol = solveODE(f, tspan, y0, { method: 'RK45', minStep }) + const minStep = 0.1 + const sol = solveODE(f, tspan, y0, { minStep }) assert.deepStrictEqual( - largerEq(min(diff(sol.t)), minStep), + // slicing to exclude the last step + largerEq(min(diff(sol.t.slice(0, sol.t.length - 1))), minStep), true ) }) From 9873270bed06eed0ca3f9282c8bcc2254fa79ea0 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Sat, 1 Jul 2023 18:01:29 -0600 Subject: [PATCH 10/25] Included error for tspan --- AUTHORS | 1 + src/function/numeric/solveODE.js | 7 ++++++- test/unit-tests/function/numeric/solveODE.test.js | 2 +- 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/AUTHORS b/AUTHORS index 94bccd9c50..a8cf64d0dd 100644 --- a/AUTHORS +++ b/AUTHORS @@ -223,5 +223,6 @@ David Contreras Jakub Riegel Angus Comrie TMTron +Maxim Mazurok # Generated by tools/update-authors.js diff --git a/src/function/numeric/solveODE.js b/src/function/numeric/solveODE.js index af0587d757..6b4ecf560f 100644 --- a/src/function/numeric/solveODE.js +++ b/src/function/numeric/solveODE.js @@ -13,6 +13,7 @@ const dependencies = [ 'abs', 'isPositive', 'isNegative', + 'isNumeric', 'larger', 'smaller', 'matrix', @@ -31,6 +32,7 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( abs, isPositive, isNegative, + isNumeric, larger, smaller, matrix, @@ -96,7 +98,10 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( return function (f, T, y0, options) { // adaptive runge kutta methods const hasBigNumbers = T.some(isBigNumber) || y0.some(isBigNumber) - + const wrongTSpan = !((isNumeric(T[0]) && isNumeric(T[1])) || (isUnit(T[0]) && isUnit(T[1]))) + if(wrongTSpan){ + throw new Error('"tspan" must be an Array of two numeric values or two units [tStart, tEnd]') + } const t0 = hasBigNumbers ? bignumber(T[0]) : T[0] // initial time const tf = hasBigNumbers ? bignumber(T[1]) : T[1] // final time const steps = 1 // divide time in this number of steps diff --git a/test/unit-tests/function/numeric/solveODE.test.js b/test/unit-tests/function/numeric/solveODE.test.js index e8a70dd03a..e7ac4a266b 100644 --- a/test/unit-tests/function/numeric/solveODE.test.js +++ b/test/unit-tests/function/numeric/solveODE.test.js @@ -44,7 +44,7 @@ describe('solveODE', function () { assert.throws(function () { const wrongTSpan = [tspan[0]] solveODE(f, wrongTSpan, y0) - }, /TypeError: Unexpected type of argument in function.*/) + }, /Error: "tspan" must be an Array of two numeric values.*/) }) it('should throw an error if the name of the method is wrong', function () { From 8071e9109092ba42a8fe7da050e761bda745e0c8 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Sun, 2 Jul 2023 08:50:27 -0600 Subject: [PATCH 11/25] It works with bignumbers --- src/function/numeric/solveODE.js | 27 ++++++++++++------- .../function/numeric/solveODE.test.js | 9 +++---- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/src/function/numeric/solveODE.js b/src/function/numeric/solveODE.js index 6b4ecf560f..3384457c38 100644 --- a/src/function/numeric/solveODE.js +++ b/src/function/numeric/solveODE.js @@ -99,15 +99,15 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( // adaptive runge kutta methods const hasBigNumbers = T.some(isBigNumber) || y0.some(isBigNumber) const wrongTSpan = !((isNumeric(T[0]) && isNumeric(T[1])) || (isUnit(T[0]) && isUnit(T[1]))) - if(wrongTSpan){ + if (wrongTSpan) { throw new Error('"tspan" must be an Array of two numeric values or two units [tStart, tEnd]') } const t0 = hasBigNumbers ? bignumber(T[0]) : T[0] // initial time const tf = hasBigNumbers ? bignumber(T[1]) : T[1] // final time const steps = 1 // divide time in this number of steps const tol = options.tol ? options.tol : 1e-4 // define a tolerance (must be an option) - const maxStep = hasBigNumbers ? bignumber(options.maxStep) : options.maxStep - const minStep = hasBigNumbers ? bignumber(options.minStep) : options.minStep + const maxStep = options.maxStep ? hasBigNumbers ? bignumber(options.maxStep) : options.maxStep : null + const minStep = options.minStep ? hasBigNumbers ? bignumber(options.minStep) : options.minStep : null const minDelta = options.minDelta ? options.minDelta : 0.2 const maxDelta = options.maxDelta ? options.maxDelta : 5 const maxIter = options.maxIter ? options.maxIter : 10_000 // stop inifite evaluation if something goes wrong @@ -167,13 +167,22 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( // estimate the delta value that will affect the step size const delta = 0.84 * (tol / TE) ** (1 / 5) - - if (delta < minDelta) { - h = multiply(h, minDelta) - } else if (delta > maxDelta) { - h = multiply(h, maxDelta) + if (hasBigNumbers) { + if (smaller(delta, minDelta)) { + h = multiply(h, bignumber(minDelta)) + } else if (larger(delta, maxDelta)) { + h = multiply(h, bignumber(maxDelta)) + } else { + h = multiply(h, bignumber(delta)) + } } else { - h = multiply(h, delta) + if (delta < minDelta) { + h = multiply(h, minDelta) + } else if (delta > maxDelta) { + h = multiply(h, maxDelta) + } else { + h = multiply(h, delta) + } } if (maxStep && larger(abs(h), abs(maxStep))) { diff --git a/test/unit-tests/function/numeric/solveODE.test.js b/test/unit-tests/function/numeric/solveODE.test.js index e7ac4a266b..977c8ef178 100644 --- a/test/unit-tests/function/numeric/solveODE.test.js +++ b/test/unit-tests/function/numeric/solveODE.test.js @@ -21,12 +21,13 @@ const map = math.map const transpose = math.transpose const isMatrix = math.isMatrix const bignumber = math.bignumber +const number = math.number describe('solveODE', function () { function f (t, y) { return subtract(y, t) } function exactSol (T, y0) { return add(dotMultiply(subtract(y0, 1), map(transpose([T]), exp)), transpose([T]), 1) } // this is only valid for y' = y - t const tspan = [0, 4] - const y0 = [1, 2] + const y0 = [1.5, 2] it('should throw an error if the first argument is not a function', function () { assert.throws(function () { @@ -80,12 +81,10 @@ describe('solveODE', function () { ) }) - /* it('should solve if the arguments have bignumbers', function () { const sol = solveODE(f, bignumber(tspan), bignumber(y0)) - approx.deepEqual(sol.y, bignumber(exactSol(sol.t, y0)), tol) + approx.deepEqual(number(sol.y), exactSol(number(sol.t), y0), tol) }) - */ it('should solve with options even if they are empty', function () { const options = {} @@ -123,7 +122,7 @@ describe('solveODE', function () { approx.deepEqual(sol.y, exactSol(sol.t, y0), tol) }) - it('should solve with less steps if a higher tolerance is specified', function () { + it('should solve with few steps if a higher tolerance is specified', function () { const sol = solveODE(f, tspan, y0, { method: 'RK45', tol: 1e-2 }) assert.deepStrictEqual( smallerEq(sol.y.length, 6), From b364fc558acf0325a118593f50ee4972e0d0f9ff Mon Sep 17 00:00:00 2001 From: David Contreras Date: Sun, 2 Jul 2023 11:22:17 -0600 Subject: [PATCH 12/25] reduced calling bignumber --- src/function/numeric/solveODE.js | 58 ++++++++++++++------------------ 1 file changed, 25 insertions(+), 33 deletions(-) diff --git a/src/function/numeric/solveODE.js b/src/function/numeric/solveODE.js index 3384457c38..1ff719fae9 100644 --- a/src/function/numeric/solveODE.js +++ b/src/function/numeric/solveODE.js @@ -95,19 +95,19 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( function _rk (butcherTableau) { // generates an adaptive runge kutta method from it's butcher tableau - return function (f, T, y0, options) { + return function (f, tspan, y0, options) { // adaptive runge kutta methods - const hasBigNumbers = T.some(isBigNumber) || y0.some(isBigNumber) - const wrongTSpan = !((isNumeric(T[0]) && isNumeric(T[1])) || (isUnit(T[0]) && isUnit(T[1]))) + const hasBigNumbers = tspan.some(isBigNumber) || y0.some(isBigNumber) + const wrongTSpan = !((tspan.length === 2) && ((isNumeric(tspan[0]) && isNumeric(tspan[1])) || (isUnit(tspan[0]) && isUnit(tspan[1])))) if (wrongTSpan) { throw new Error('"tspan" must be an Array of two numeric values or two units [tStart, tEnd]') } - const t0 = hasBigNumbers ? bignumber(T[0]) : T[0] // initial time - const tf = hasBigNumbers ? bignumber(T[1]) : T[1] // final time + const t0 = tspan[0] // initial time + const tf = tspan[1] // final time const steps = 1 // divide time in this number of steps const tol = options.tol ? options.tol : 1e-4 // define a tolerance (must be an option) - const maxStep = options.maxStep ? hasBigNumbers ? bignumber(options.maxStep) : options.maxStep : null - const minStep = options.minStep ? hasBigNumbers ? bignumber(options.minStep) : options.minStep : null + const maxStep = options.maxStep + const minStep = options.minStep const minDelta = options.minDelta ? options.minDelta : 0.2 const maxDelta = options.maxDelta ? options.maxDelta : 5 const maxIter = options.maxIter ? options.maxIter : 10_000 // stop inifite evaluation if something goes wrong @@ -120,10 +120,10 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( ] : [butcherTableau.a, butcherTableau.c, butcherTableau.b, butcherTableau.bp] let h = options.firstStep - ? hasBigNumbers ? bignumber(options.firstStep) : options.firstStep + ? options.firstStep : divide(subtract(tf, t0), steps) // define the first step size - const t = hasBigNumbers ? [bignumber(t0)] : [t0] // start the time array - const y = hasBigNumbers ? [bignumber(y0)] : [y0] // start the solution array + const t = [t0] // start the time array + const y = [y0] // start the solution array const dletaB = subtract(b, bp) // b - bp @@ -166,25 +166,17 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( } // estimate the delta value that will affect the step size - const delta = 0.84 * (tol / TE) ** (1 / 5) - if (hasBigNumbers) { - if (smaller(delta, minDelta)) { - h = multiply(h, bignumber(minDelta)) - } else if (larger(delta, maxDelta)) { - h = multiply(h, bignumber(maxDelta)) - } else { - h = multiply(h, bignumber(delta)) - } - } else { - if (delta < minDelta) { - h = multiply(h, minDelta) - } else if (delta > maxDelta) { - h = multiply(h, maxDelta) - } else { - h = multiply(h, delta) - } + let delta = 0.84 * (tol / TE) ** (1 / 5) + + if (smaller(delta, minDelta)) { + delta = minDelta + } else if (larger(delta, maxDelta)) { + delta = maxDelta } + delta = hasBigNumbers ? bignumber(delta) : delta + h = multiply(h, delta) + if (maxStep && larger(abs(h), abs(maxStep))) { h = isPositive(h) ? abs(maxStep) : multiply(-1, abs(maxStep)) } else if (minStep && smaller(abs(h), abs(minStep))) { @@ -199,7 +191,7 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( } } - function _rk23 (f, T, y0, options) { + function _rk23 (f, tspan, y0, options) { // Bogacki–Shampine method // Define the butcher table @@ -217,10 +209,10 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( const butcherTableau = { a, c, b, bp } // Solve an adaptive step size rk method - return _rk(butcherTableau)(f, T, y0, options) + return _rk(butcherTableau)(f, tspan, y0, options) } - function _rk45 (f, T, y0, options) { + function _rk45 (f, tspan, y0, options) { // Dormand Prince method // Define the butcher tableau @@ -241,10 +233,10 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( const butcherTableau = { a, c, b, bp } // Solve an adaptive step size rk method - return _rk(butcherTableau)(f, T, y0, options) + return _rk(butcherTableau)(f, tspan, y0, options) } - function _solveODE (f, T, y0, opt) { + function _solveODE (f, tspan, y0, opt) { const method = opt.method ? opt.method : 'RK45' const methods = { RK23: _rk23, @@ -253,7 +245,7 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( if (method.toUpperCase() in methods) { const methodOptions = { ...opt } // clone the options object delete methodOptions.method // delete the method as it won't be needed - return methods[method.toUpperCase()](f, T, y0, methodOptions) + return methods[method.toUpperCase()](f, tspan, y0, methodOptions) } else { // throw an error indicating there is no such method const methodsWithQuotes = Object.keys(methods).map(x => `"${x}"`) From d2a9cb813cb0c1774f3215782c85c7d2243e3693 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Sun, 2 Jul 2023 11:36:49 -0600 Subject: [PATCH 13/25] extended the search for bignumbers --- src/function/numeric/solveODE.js | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/function/numeric/solveODE.js b/src/function/numeric/solveODE.js index 1ff719fae9..08585bbdd3 100644 --- a/src/function/numeric/solveODE.js +++ b/src/function/numeric/solveODE.js @@ -97,7 +97,7 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( return function (f, tspan, y0, options) { // adaptive runge kutta methods - const hasBigNumbers = tspan.some(isBigNumber) || y0.some(isBigNumber) + const hasBigNumbers = [...tspan, ...y0, options.maxStep, options.minStep].some(isBigNumber) const wrongTSpan = !((tspan.length === 2) && ((isNumeric(tspan[0]) && isNumeric(tspan[1])) || (isUnit(tspan[0]) && isUnit(tspan[1])))) if (wrongTSpan) { throw new Error('"tspan" must be an Array of two numeric values or two units [tStart, tEnd]') From 6f56580a9d9e45cfcf0a6f0518cd4891f408c878 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Wed, 5 Jul 2023 12:53:20 -0600 Subject: [PATCH 14/25] The jsdocs is less ambiguous --- src/function/numeric/solveODE.js | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/function/numeric/solveODE.js b/src/function/numeric/solveODE.js index 08585bbdd3..f4b47b08bb 100644 --- a/src/function/numeric/solveODE.js +++ b/src/function/numeric/solveODE.js @@ -53,8 +53,8 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( * - `y0` the initial state values, should be a scalar or a flat array * - `options` should be an object with the following information: * - `method` ('RK45'): ['RK23', 'RK45'] - * - `tol` (1e-3): A numeric value - * - `firstStep`: A numeric value or unit + * - `tol` (1e-3): Numeric tolerance of the method, the solver keeps the error estimates less than this value + * - `firstStep`: Initial step size * - `minStep`: minimum step size of the method * - `maxStep`: maximum step size of the method * - `minDelta` (0.2): minimum ratio of change for the step From 9448e8883b5c2c8d9907254f23b3344bdeb31c72 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Wed, 5 Jul 2023 22:34:44 -0600 Subject: [PATCH 15/25] included tests for step options --- src/function/numeric/solveODE.js | 67 ++++++++++------- .../function/numeric/solveODE.test.js | 72 +++++++++++++++++++ 2 files changed, 115 insertions(+), 24 deletions(-) diff --git a/src/function/numeric/solveODE.js b/src/function/numeric/solveODE.js index f4b47b08bb..47a039624c 100644 --- a/src/function/numeric/solveODE.js +++ b/src/function/numeric/solveODE.js @@ -1,4 +1,4 @@ -import { isUnit, isBigNumber } from '../../utils/is.js' +import { isUnit, isNumber, isBigNumber } from '../../utils/is.js' import { factory } from '../../utils/factory.js' const name = 'solveODE' @@ -12,12 +12,11 @@ const dependencies = [ 'map', 'abs', 'isPositive', - 'isNegative', - 'isNumeric', 'larger', 'smaller', 'matrix', - 'bignumber' + 'bignumber', + 'unaryMinus' ] export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( @@ -31,12 +30,11 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( map, abs, isPositive, - isNegative, - isNumeric, larger, smaller, matrix, - bignumber + bignumber, + unaryMinus } ) => { /** @@ -97,20 +95,35 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( return function (f, tspan, y0, options) { // adaptive runge kutta methods - const hasBigNumbers = [...tspan, ...y0, options.maxStep, options.minStep].some(isBigNumber) - const wrongTSpan = !((tspan.length === 2) && ((isNumeric(tspan[0]) && isNumeric(tspan[1])) || (isUnit(tspan[0]) && isUnit(tspan[1])))) + const wrongTSpan = !((tspan.length === 2) && (tspan.every(isNumOrBig) || tspan.every(isUnit))) if (wrongTSpan) { throw new Error('"tspan" must be an Array of two numeric values or two units [tStart, tEnd]') } const t0 = tspan[0] // initial time const tf = tspan[1] // final time - const steps = 1 // divide time in this number of steps - const tol = options.tol ? options.tol : 1e-4 // define a tolerance (must be an option) + const isForwards = larger(tf, t0) + const firstStep = options.firstStep + if (firstStep !== undefined && !isPositive(firstStep)) { + throw new Error('"firstStep" must be positive') + } const maxStep = options.maxStep + if (maxStep !== undefined && !isPositive(maxStep)) { + throw new Error('"maxStep" must be positive') + } const minStep = options.minStep + if (minStep !== undefined && !isPositive(minStep)) { + throw new Error('"minStep" must be positive') + } + const timeVars = [t0, tf, firstStep, minStep, maxStep].filter(x => x !== undefined) + if (!(timeVars.every(isNumOrBig) || timeVars.every(isUnit))) { + throw new Error('Inconsistent type of "t" dependant variables') + } + const steps = 1 // divide time in this number of steps + const tol = options.tol ? options.tol : 1e-4 // define a tolerance (must be an option) const minDelta = options.minDelta ? options.minDelta : 0.2 const maxDelta = options.maxDelta ? options.maxDelta : 5 const maxIter = options.maxIter ? options.maxIter : 10_000 // stop inifite evaluation if something goes wrong + const hasBigNumbers = [t0, tf, ...y0, maxStep, minStep].some(isBigNumber) const [a, c, b, bp] = hasBigNumbers ? [ bignumber(butcherTableau.a), @@ -119,13 +132,14 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( bignumber(butcherTableau.bp) ] : [butcherTableau.a, butcherTableau.c, butcherTableau.b, butcherTableau.bp] - let h = options.firstStep - ? options.firstStep + + let h = firstStep + ? isForwards ? firstStep : unaryMinus(firstStep) : divide(subtract(tf, t0), steps) // define the first step size const t = [t0] // start the time array const y = [y0] // start the solution array - const dletaB = subtract(b, bp) // b - bp + const deltaB = subtract(b, bp) // b - bp let n = 0 let iter = 0 @@ -134,7 +148,7 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( const k = [] // trim the time step so that it doesn't overshoot - h = trimStep(t[n], tf, h) + h = trimStep(t[n], tf, h, isForwards) // calculate the first value of k k.push(f(t[n], y[n])) @@ -152,14 +166,14 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( // estimate the error by comparing solutions of different orders const TE = max( abs( - map(multiply(dletaB, k), (X) => + map(multiply(deltaB, k), (X) => isUnit(X) ? X.value : X ) ) ) if (TE < tol && tol / TE > 1 / 4) { - // if within tol then push everything + // push solution if within tol t.push(add(t[n], h)) y.push(add(y[n], multiply(h, b, k))) n++ @@ -177,10 +191,10 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( delta = hasBigNumbers ? bignumber(delta) : delta h = multiply(h, delta) - if (maxStep && larger(abs(h), abs(maxStep))) { - h = isPositive(h) ? abs(maxStep) : multiply(-1, abs(maxStep)) - } else if (minStep && smaller(abs(h), abs(minStep))) { - h = isPositive(h) ? abs(minStep) : multiply(-1, abs(minStep)) + if (maxStep && larger(abs(h), maxStep)) { + h = isForwards ? maxStep : unaryMinus(maxStep) + } else if (minStep && smaller(abs(h), minStep)) { + h = isForwards ? minStep : unaryMinus(minStep) } iter++ if (iter > maxIter) { @@ -262,17 +276,22 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( : larger(t, tf) } - function trimStep (t, tf, h) { + function trimStep (t, tf, h, isForwards) { // Trims the time step so that the next step doesn't overshoot const next = add(t, h) return ( - (isPositive(h) && larger(next, tf)) || - (isNegative(h) && smaller(next, tf)) + (isForwards && larger(next, tf)) || + (!isForwards && smaller(next, tf)) ) ? subtract(tf, t) : h } + function isNumOrBig (x) { + // checks if it's a number or bignumber + return isBigNumber(x) || isNumber(x) + } + function _matrixSolveODE (f, T, y0, options) { // receives matrices and returns matrices const sol = _solveODE(f, T.toArray(), y0.toArray(), options) diff --git a/test/unit-tests/function/numeric/solveODE.test.js b/test/unit-tests/function/numeric/solveODE.test.js index 977c8ef178..16fc8852eb 100644 --- a/test/unit-tests/function/numeric/solveODE.test.js +++ b/test/unit-tests/function/numeric/solveODE.test.js @@ -61,6 +61,78 @@ describe('solveODE', function () { }, /Maximum number of iterations reached.*/) }) + it('should throw an error if the firstStep is not positive', function () { + assert.throws(function () { + solveODE(f, tspan, y0, { firstStep: -1 }) + }, /"firstStep" must be positive/) + assert.throws(function () { + solveODE(f, tspan, y0, { firstStep: 0 }) + }, /"firstStep" must be positive/) + assert.throws(function () { + solveODE(f, tspan, y0, { firstStep: unit(-1, 's') }) + }, /"firstStep" must be positive/) + assert.throws(function () { + solveODE(f, tspan, y0, { firstStep: unit(0, 's') }) + }, /"firstStep" must be positive/) + }) + + it('should throw an error if the minStep is not positive', function () { + assert.throws(function () { + solveODE(f, tspan, y0, { minStep: -1 }) + }, /"minStep" must be positive/) + assert.throws(function () { + solveODE(f, tspan, y0, { minStep: 0 }) + }, /"minStep" must be positive/) + assert.throws(function () { + solveODE(f, tspan, y0, { minStep: unit(-1, 's') }) + }, /"minStep" must be positive/) + assert.throws(function () { + solveODE(f, tspan, y0, { minStep: unit(0, 's') }) + }, /"minStep" must be positive/) + }) + + it('should throw an error if the maxStep is not positive', function () { + assert.throws(function () { + solveODE(f, tspan, y0, { maxStep: -1 }) + }, /"maxStep" must be positive/) + assert.throws(function () { + solveODE(f, tspan, y0, { maxStep: 0 }) + }, /"maxStep" must be positive/) + assert.throws(function () { + solveODE(f, tspan, y0, { maxStep: unit(-1, 's') }) + }, /"maxStep" must be positive/) + assert.throws(function () { + solveODE(f, tspan, y0, { maxStep: unit(0, 's') }) + }, /"maxStep" must be positive/) + }) + + it('should throw an error if the time dependant variables do not match', function () { + assert.throws(function () { + solveODE(f, tspan, y0, { firstStep: unit(1, 's') }) + }, /Inconsistent type of "t" dependant variables/) + assert.throws(function () { + solveODE(f, [unit(tspan[0], 's'), unit(tspan[1], 's')], y0, { maxStep: unit(2, 's'), firstStep: 1 }) + }, /Inconsistent type of "t" dependant variables/) + assert.throws(function () { + solveODE(f, tspan, y0, { firstStep: unit(1, 's') }) + }, /Inconsistent type of "t" dependant variables/) + assert.throws(function () { + solveODE(f, [unit(tspan[0], 's'), unit(tspan[1], 's')], y0, { firstStep: 1 }) + }, /Inconsistent type of "t" dependant variables/) + assert.throws(function () { + solveODE(f, tspan, y0, { maxStep: unit(1, 's') }) + }, /Inconsistent type of "t" dependant variables/) + assert.throws(function () { + solveODE(f, [unit(tspan[0], 's'), unit(tspan[1], 's')], y0, { maxStep: 1 }) + }, /Inconsistent type of "t" dependant variables/) + assert.throws(function () { + solveODE(f, tspan, y0, { minStep: unit(1, 's') }) + }, /Inconsistent type of "t" dependant variables/) + assert.throws(function () { + solveODE(f, [unit(tspan[0], 's'), unit(tspan[1], 's')], y0, { minStep: 1 }) + }, /Inconsistent type of "t" dependant variables/) + }) + it('should solve close to the analytical solution', function () { const sol = solveODE(f, tspan, y0) approx.deepEqual(sol.y, exactSol(sol.t, y0), tol) From a827c7b49a99b51a443c794e239da793862e25ac Mon Sep 17 00:00:00 2001 From: David Contreras Date: Wed, 5 Jul 2023 23:02:59 -0600 Subject: [PATCH 16/25] Allowed for 0 minStep --- src/function/numeric/solveODE.js | 6 +++-- .../function/numeric/solveODE.test.js | 24 ++++++++++++------- 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/src/function/numeric/solveODE.js b/src/function/numeric/solveODE.js index 47a039624c..7d54985034 100644 --- a/src/function/numeric/solveODE.js +++ b/src/function/numeric/solveODE.js @@ -12,6 +12,7 @@ const dependencies = [ 'map', 'abs', 'isPositive', + 'isNegative', 'larger', 'smaller', 'matrix', @@ -30,6 +31,7 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( map, abs, isPositive, + isNegative, larger, smaller, matrix, @@ -111,8 +113,8 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( throw new Error('"maxStep" must be positive') } const minStep = options.minStep - if (minStep !== undefined && !isPositive(minStep)) { - throw new Error('"minStep" must be positive') + if (minStep && isNegative(minStep)) { + throw new Error('"minStep" must be positive or zero') } const timeVars = [t0, tf, firstStep, minStep, maxStep].filter(x => x !== undefined) if (!(timeVars.every(isNumOrBig) || timeVars.every(isUnit))) { diff --git a/test/unit-tests/function/numeric/solveODE.test.js b/test/unit-tests/function/numeric/solveODE.test.js index 16fc8852eb..bf02922449 100644 --- a/test/unit-tests/function/numeric/solveODE.test.js +++ b/test/unit-tests/function/numeric/solveODE.test.js @@ -76,19 +76,25 @@ describe('solveODE', function () { }, /"firstStep" must be positive/) }) - it('should throw an error if the minStep is not positive', function () { + it('should throw an error if the minStep is not positive or zero', function () { assert.throws(function () { solveODE(f, tspan, y0, { minStep: -1 }) - }, /"minStep" must be positive/) - assert.throws(function () { - solveODE(f, tspan, y0, { minStep: 0 }) - }, /"minStep" must be positive/) + }, /"minStep" must be positive or zero/) assert.throws(function () { solveODE(f, tspan, y0, { minStep: unit(-1, 's') }) - }, /"minStep" must be positive/) - assert.throws(function () { - solveODE(f, tspan, y0, { minStep: unit(0, 's') }) - }, /"minStep" must be positive/) + }, /"minStep" must be positive or zero/) + }) + + it('should solve when minStep is zero', function () { + const sol = solveODE(f, tspan, y0, { minStep: 0 }) + approx.deepEqual(sol.y, exactSol(sol.t, y0), tol) + const seconds = unit('s') + const meters = unit('m') + function fWithUnits (t, y) { return subtract(divide(y, seconds), multiply(t, divide(meters, multiply(seconds, seconds)))) } + const tspanWithUnits = multiply(tspan, seconds) + const y0withUnits = multiply(y0, meters) + const solU = solveODE(fWithUnits, tspanWithUnits, y0withUnits, { minStep: unit(0, 's') }) + approx.deepEqual(divide(solU.y, meters), exactSol(divide(solU.t, seconds), y0), tol) }) it('should throw an error if the maxStep is not positive', function () { From 61f063fc440352e002e8feb962b7a5727cd58fe9 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Mon, 10 Jul 2023 15:55:33 -0600 Subject: [PATCH 17/25] Optimization to avoid checking the sign every step --- src/function/numeric/solveODE.js | 29 +++++++++++++---------------- 1 file changed, 13 insertions(+), 16 deletions(-) diff --git a/src/function/numeric/solveODE.js b/src/function/numeric/solveODE.js index 7d54985034..76b61f4c66 100644 --- a/src/function/numeric/solveODE.js +++ b/src/function/numeric/solveODE.js @@ -145,12 +145,14 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( let n = 0 let iter = 0 + const ongoing = _createOngoing(isForwards) + const trimStep = _createTrimStep(isForwards) // iterate unitil it reaches either the final time or maximum iterations - while (ongoing(t[n], tf, h)) { + while (ongoing(t[n], tf)) { const k = [] // trim the time step so that it doesn't overshoot - h = trimStep(t[n], tf, h, isForwards) + h = trimStep(t[n], tf, h) // calculate the first value of k k.push(f(t[n], y[n])) @@ -271,22 +273,17 @@ export const createSolveODE = /* #__PURE__ */ factory(name, dependencies, ( } } - function ongoing (t, tf, h) { - // returns true if the time has not reached tf for both postitive an negative step (h) - return isPositive(h) - ? smaller(t, tf) - : larger(t, tf) + function _createOngoing (isForwards) { + // returns the correct function to test if it's still iterating + return isForwards ? smaller : larger } - function trimStep (t, tf, h, isForwards) { - // Trims the time step so that the next step doesn't overshoot - const next = add(t, h) - return ( - (isForwards && larger(next, tf)) || - (!isForwards && smaller(next, tf)) - ) - ? subtract(tf, t) - : h + function _createTrimStep (isForwards) { + const outOfBounds = isForwards ? larger : smaller + return function (t, tf, h) { + const next = add(t, h) + return outOfBounds(next, tf) ? subtract(tf, t) : h + } } function isNumOrBig (x) { From 41bd5c8dd9a6ae83fea6620a84cfb5fbbf7aac69 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Sat, 29 Jul 2023 12:15:15 -0600 Subject: [PATCH 18/25] refactor --- src/function/statistics/quantileSeq.js | 203 ++++++++---------- .../function/statistics/quantileSeq.test.js | 6 +- 2 files changed, 94 insertions(+), 115 deletions(-) diff --git a/src/function/statistics/quantileSeq.js b/src/function/statistics/quantileSeq.js index dbf8a351dd..6e68563932 100644 --- a/src/function/statistics/quantileSeq.js +++ b/src/function/statistics/quantileSeq.js @@ -1,4 +1,4 @@ -import { isBigNumber, isCollection, isNumber } from '../../utils/is.js' +import { isBigNumber, isNumber } from '../../utils/is.js' import { isInteger } from '../../utils/number.js' import { flatten } from '../../utils/array.js' import { factory } from '../../utils/factory.js' @@ -41,109 +41,114 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ * @param {Boolean} sorted=false is data sorted in ascending order * @return {Number, BigNumber, Unit, Array} Quantile(s) */ - function quantileSeq (data, probOrN, sorted) { - let probArr, dataArr, one - if (arguments.length < 2 || arguments.length > 3) { - throw new SyntaxError('Function quantileSeq requires two or three parameters') - } + return typed(name, { + 'Array | Matrix, number | BigNumber | Unit': (data, p) => _quantileSeqProbNumber(data, p, false), + 'Array | Matrix, number | BigNumber | Unit, boolean': _quantileSeqProbNumber, + 'Array | Matrix, Array | Matrix': (data, p) => _quantileSeqProbCoolection(data, p, false), + 'Array | Matrix, Array | Matrix, boolean': _quantileSeqProbCoolection + }) - if (isCollection(data)) { - sorted = sorted || false - if (typeof sorted === 'boolean') { - dataArr = data.valueOf() - if (isNumber(probOrN)) { - if (probOrN < 0) { - throw new Error('N/prob must be non-negative') - } + function _quantileSeqProbNumber (data, probOrN, sorted) { + let probArr, one + const dataArr = data.valueOf() + if (isNumber(probOrN)) { + if (probOrN < 0) { + throw new Error('N/prob must be non-negative') + } - if (probOrN <= 1) { - // quantileSeq([a, b, c, d, ...], prob[,sorted]) - return _quantileSeq(dataArr, probOrN, sorted) - } + if (probOrN <= 1) { + // quantileSeq([a, b, c, d, ...], prob[,sorted]) + return _quantileSeq(dataArr, probOrN, sorted) + } - if (probOrN > 1) { - // quantileSeq([a, b, c, d, ...], N[,sorted]) - if (!isInteger(probOrN)) { - throw new Error('N must be a positive integer') - } - - const nPlusOne = probOrN + 1 - probArr = new Array(probOrN) - for (let i = 0; i < probOrN;) { - probArr[i] = _quantileSeq(dataArr, (++i) / nPlusOne, sorted) - } - return probArr - } + if (probOrN > 1) { + // quantileSeq([a, b, c, d, ...], N[,sorted]) + if (!isInteger(probOrN)) { + throw new Error('N must be a positive integer') + } + + const nPlusOne = probOrN + 1 + probArr = new Array(probOrN) + for (let i = 0; i < probOrN; i++) { + probArr[i] = _quantileSeq(dataArr, (i + 1) / nPlusOne, sorted) } + return probArr + } + } - if (isBigNumber(probOrN)) { - const BigNumber = probOrN.constructor + if (isBigNumber(probOrN)) { + const BigNumber = probOrN.constructor - if (probOrN.isNegative()) { - throw new Error('N/prob must be non-negative') - } + if (probOrN.isNegative()) { + throw new Error('N/prob must be non-negative') + } - one = new BigNumber(1) + one = new BigNumber(1) - if (probOrN.lte(one)) { - // quantileSeq([a, b, c, d, ...], prob[,sorted]) - return new BigNumber(_quantileSeq(dataArr, probOrN, sorted)) - } + if (probOrN.lte(one)) { + // quantileSeq([a, b, c, d, ...], prob[,sorted]) + return new BigNumber(_quantileSeq(dataArr, probOrN, sorted)) + } - if (probOrN.gt(one)) { - // quantileSeq([a, b, c, d, ...], N[,sorted]) - if (!probOrN.isInteger()) { - throw new Error('N must be a positive integer') - } - - // largest possible Array length is 2^32-1 - // 2^32 < 10^15, thus safe conversion guaranteed - const intN = probOrN.toNumber() - if (intN > 4294967295) { - throw new Error('N must be less than or equal to 2^32-1, as that is the maximum length of an Array') - } - - const nPlusOne = new BigNumber(intN + 1) - probArr = new Array(intN) - for (let i = 0; i < intN;) { - probArr[i] = new BigNumber(_quantileSeq(dataArr, new BigNumber(++i).div(nPlusOne), sorted)) - } - return probArr - } + if (probOrN.gt(one)) { + // quantileSeq([a, b, c, d, ...], N[,sorted]) + if (!probOrN.isInteger()) { + throw new Error('N must be a positive integer') } - if (isCollection(probOrN)) { - // quantileSeq([a, b, c, d, ...], [prob1, prob2, ...][,sorted]) - const probOrNArr = probOrN.valueOf() - probArr = new Array(probOrNArr.length) - for (let i = 0; i < probArr.length; ++i) { - const currProb = probOrNArr[i] - if (isNumber(currProb)) { - if (currProb < 0 || currProb > 1) { - throw new Error('Probability must be between 0 and 1, inclusive') - } - } else if (isBigNumber(currProb)) { - one = new currProb.constructor(1) - if (currProb.isNegative() || currProb.gt(one)) { - throw new Error('Probability must be between 0 and 1, inclusive') - } - } else { - throw new TypeError('Unexpected type of argument in function quantileSeq') // FIXME: becomes redundant when converted to typed-function - } - - probArr[i] = _quantileSeq(dataArr, currProb, sorted) - } - return probArr + // largest possible Array length is 2^32-1 + // 2^32 < 10^15, thus safe conversion guaranteed + const intN = probOrN.toNumber() + if (intN > 4294967295) { + throw new Error('N must be less than or equal to 2^32-1, as that is the maximum length of an Array') + } + + const nPlusOne = new BigNumber(intN + 1) + probArr = new Array(intN) + for (let i = 0; i < intN; i++) { + probArr[i] = new BigNumber(_quantileSeq(dataArr, new BigNumber(i + 1).div(nPlusOne), sorted)) } + return probArr + } + } + } + /** + * Calculate the prob order quantile of an n-dimensional array. + * + * @param {Array, Matrix} array + * @param {Array, Matrix} prob + * @param {Boolean} sorted + * @return {Number, BigNumber, Unit} prob order quantile + * @private + */ + + function _quantileSeqProbCoolection (data, probOrN, sorted) { + const dataArr = data.valueOf() + let one + + // quantileSeq([a, b, c, d, ...], [prob1, prob2, ...][,sorted]) + const probOrNArr = probOrN.valueOf() + const probArr = new Array(probOrNArr.length) + for (let i = 0; i < probArr.length; ++i) { + const currProb = probOrNArr[i] + if (isNumber(currProb)) { + if (currProb < 0 || currProb > 1) { + throw new Error('Probability must be between 0 and 1, inclusive') + } + } else if (isBigNumber(currProb)) { + one = new currProb.constructor(1) + if (currProb.isNegative() || currProb.gt(one)) { + throw new Error('Probability must be between 0 and 1, inclusive') + } + } else { throw new TypeError('Unexpected type of argument in function quantileSeq') // FIXME: becomes redundant when converted to typed-function } - throw new TypeError('Unexpected type of argument in function quantileSeq') // FIXME: becomes redundant when converted to typed-function + probArr[i] = _quantileSeq(dataArr, currProb, sorted) } - - throw new TypeError('Unexpected type of argument in function quantileSeq') // FIXME: becomes redundant when converted to typed-function + return probArr } /** @@ -167,9 +172,6 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ const fracPart = index % 1 if (fracPart === 0) { const value = sorted ? flat[index] : partitionSelect(flat, index) - - validate(value) - return value } @@ -192,9 +194,6 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ } } - validate(left) - validate(right) - // Q(prob) = (1-f)*A[floor(index)] + f*A[floor(index)+1] return add(multiply(left, 1 - fracPart), multiply(right, fracPart)) } @@ -204,9 +203,6 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ if (index.isInteger()) { index = index.toNumber() const value = sorted ? flat[index] : partitionSelect(flat, index) - - validate(value) - return value } @@ -231,25 +227,8 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ } } - validate(left) - validate(right) - // Q(prob) = (1-f)*A[floor(index)] + f*A[floor(index)+1] const one = new fracPart.constructor(1) return add(multiply(left, one.minus(fracPart)), multiply(right, fracPart)) } - - /** - * Check if array value types are valid, throw error otherwise. - * @param {number | BigNumber | Unit} x - * @param {number | BigNumber | Unit} x - * @private - */ - const validate = typed({ - 'number | BigNumber | Unit': function (x) { - return x - } - }) - - return quantileSeq }) diff --git a/test/unit-tests/function/statistics/quantileSeq.test.js b/test/unit-tests/function/statistics/quantileSeq.test.js index 65872cbd5e..7d89d0bed0 100644 --- a/test/unit-tests/function/statistics/quantileSeq.test.js +++ b/test/unit-tests/function/statistics/quantileSeq.test.js @@ -151,9 +151,9 @@ describe('quantileSeq', function () { }) it('should throw an error if called with invalid number of arguments', function () { - assert.throws(function () { quantileSeq() }, SyntaxError) - assert.throws(function () { quantileSeq(2) }, SyntaxError) - assert.throws(function () { quantileSeq([], 2, 3, 1) }, SyntaxError) + assert.throws(function () { quantileSeq() }, TypeError) + assert.throws(function () { quantileSeq(2) }, TypeError) + assert.throws(function () { quantileSeq([], 2, 3, 1) }, TypeError) }) it('should throw an error if called with unsupported type of arguments', function () { From c9517b04bae2ddd6fc0edc026635aeda01c4318e Mon Sep 17 00:00:00 2001 From: David Contreras Date: Sat, 29 Jul 2023 12:25:17 -0600 Subject: [PATCH 19/25] Typo --- src/function/statistics/quantileSeq.js | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/function/statistics/quantileSeq.js b/src/function/statistics/quantileSeq.js index 6e68563932..4c29ba4368 100644 --- a/src/function/statistics/quantileSeq.js +++ b/src/function/statistics/quantileSeq.js @@ -45,8 +45,8 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ return typed(name, { 'Array | Matrix, number | BigNumber | Unit': (data, p) => _quantileSeqProbNumber(data, p, false), 'Array | Matrix, number | BigNumber | Unit, boolean': _quantileSeqProbNumber, - 'Array | Matrix, Array | Matrix': (data, p) => _quantileSeqProbCoolection(data, p, false), - 'Array | Matrix, Array | Matrix, boolean': _quantileSeqProbCoolection + 'Array | Matrix, Array | Matrix': (data, p) => _quantileSeqProbCollection(data, p, false), + 'Array | Matrix, Array | Matrix, boolean': _quantileSeqProbCollection }) function _quantileSeqProbNumber (data, probOrN, sorted) { @@ -124,7 +124,7 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ * @private */ - function _quantileSeqProbCoolection (data, probOrN, sorted) { + function _quantileSeqProbCollection (data, probOrN, sorted) { const dataArr = data.valueOf() let one From 8e4d8997a2fa11001f61eff72b9ae546736d49ab Mon Sep 17 00:00:00 2001 From: David Contreras Date: Sat, 29 Jul 2023 12:33:20 -0600 Subject: [PATCH 20/25] removed unnecesary error --- src/function/statistics/quantileSeq.js | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/function/statistics/quantileSeq.js b/src/function/statistics/quantileSeq.js index 4c29ba4368..6f2bff923e 100644 --- a/src/function/statistics/quantileSeq.js +++ b/src/function/statistics/quantileSeq.js @@ -142,8 +142,6 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ if (currProb.isNegative() || currProb.gt(one)) { throw new Error('Probability must be between 0 and 1, inclusive') } - } else { - throw new TypeError('Unexpected type of argument in function quantileSeq') // FIXME: becomes redundant when converted to typed-function } probArr[i] = _quantileSeq(dataArr, currProb, sorted) From 6134f1880148946c1ce51a1904a4db08a08b5d50 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Thu, 24 Aug 2023 02:04:04 +0000 Subject: [PATCH 21/25] Fixes conflict with develop --- .../transform/quantileSeq.transform.js | 27 ++++++++++++++ src/factoriesAny.js | 1 + src/function/statistics/quantileSeq.js | 19 +++++++--- .../function/statistics/quantileSeq.test.js | 36 ++++++++++++++++++- 4 files changed, 78 insertions(+), 5 deletions(-) create mode 100644 src/expression/transform/quantileSeq.transform.js diff --git a/src/expression/transform/quantileSeq.transform.js b/src/expression/transform/quantileSeq.transform.js new file mode 100644 index 0000000000..aa83d759fa --- /dev/null +++ b/src/expression/transform/quantileSeq.transform.js @@ -0,0 +1,27 @@ +import { factory } from '../../utils/factory.js' +import { createQuantileSeq } from '../../function/statistics/quantileSeq.js' +import { lastDimToZeroBase } from './utils/lastDimToZeroBase.js' + +const name = 'quantileSeq' +const dependencies = ['typed', 'add', 'multiply', 'partitionSelect', 'compare', 'isInteger'] + +/** + * Attach a transform function to math.quantileSeq + * Adds a property transform containing the transform function. + * + * This transform changed the `dim` parameter of function std + * from one-based to zero based + */ +export const createQuantileSeqTransform = /* #__PURE__ */ factory(name, dependencies, ({ typed, add, multiply, partitionSelect, compare, isInteger }) => { + const quantileSeq = createQuantileSeq({ typed, add, multiply, partitionSelect, compare, isInteger }) + + return typed('quantileSeq', { + 'Array|Matrix, number|BigNumber|Array, number': (arr, prob, dim) => quantileSeq(arr, prob, dimToZeroBase(dim)), + 'Array|Matrix, number|BigNumber|Array, boolean, number': (arr, prob, sorted, dim) => quantileSeq(arr, prob, sorted, dimToZeroBase(dim)) + }) + + function dimToZeroBase (dim) { + // TODO: find a better way, maybe lastDimToZeroBase could apply to more cases. + return lastDimToZeroBase([[], dim])[1] + } +}, { isTransformFunction: true }) diff --git a/src/factoriesAny.js b/src/factoriesAny.js index 704c25c904..00b5fe195c 100644 --- a/src/factoriesAny.js +++ b/src/factoriesAny.js @@ -352,5 +352,6 @@ export { createConcatTransform } from './expression/transform/concat.transform.j export { createDiffTransform } from './expression/transform/diff.transform.js' export { createStdTransform } from './expression/transform/std.transform.js' export { createSumTransform } from './expression/transform/sum.transform.js' +export { createQuantileSeqTransform } from './expression/transform/quantileSeq.transform.js' export { createCumSumTransform } from './expression/transform/cumsum.transform.js' export { createVarianceTransform } from './expression/transform/variance.transform.js' diff --git a/src/function/statistics/quantileSeq.js b/src/function/statistics/quantileSeq.js index 6f2bff923e..3117073cc6 100644 --- a/src/function/statistics/quantileSeq.js +++ b/src/function/statistics/quantileSeq.js @@ -1,12 +1,12 @@ import { isBigNumber, isNumber } from '../../utils/is.js' -import { isInteger } from '../../utils/number.js' import { flatten } from '../../utils/array.js' import { factory } from '../../utils/factory.js' +import { createApply } from '../matrix/apply.js' const name = 'quantileSeq' -const dependencies = ['typed', 'add', 'multiply', 'partitionSelect', 'compare'] +const dependencies = ['typed', 'add', 'multiply', 'partitionSelect', 'compare', 'isInteger'] -export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ typed, add, multiply, partitionSelect, compare }) => { +export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ typed, add, multiply, partitionSelect, compare, isInteger }) => { /** * Compute the prob order quantile of a matrix or a list with values. * The sequence is sorted and the middle value is returned. @@ -42,13 +42,24 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ * @return {Number, BigNumber, Unit, Array} Quantile(s) */ + const apply = createApply({ typed, isInteger }) + return typed(name, { 'Array | Matrix, number | BigNumber | Unit': (data, p) => _quantileSeqProbNumber(data, p, false), + 'Array | Matrix, number | BigNumber | Unit, number': (data, prob, dim) => _quantileSeqDim(data, prob, false, dim, _quantileSeqProbNumber), 'Array | Matrix, number | BigNumber | Unit, boolean': _quantileSeqProbNumber, + 'Array | Matrix, number | BigNumber | Unit, boolean, number': (data, prob, sorted, dim) => _quantileSeqDim(data, prob, sorted, dim, _quantileSeqProbNumber), 'Array | Matrix, Array | Matrix': (data, p) => _quantileSeqProbCollection(data, p, false), - 'Array | Matrix, Array | Matrix, boolean': _quantileSeqProbCollection + 'Array | Matrix, Array | Matrix, number': (data, prob, dim) => _quantileSeqDim(data, prob, false, dim, _quantileSeqProbCollection), + 'Array | Matrix, Array | Matrix, boolean': _quantileSeqProbCollection, + 'Array | Matrix, Array | Matrix, boolean, number': (data, prob, sorted, dim) => _quantileSeqDim(data, prob, sorted, dim, _quantileSeqProbCollection) }) + function _quantileSeqDim (data, prob, sorted, dim, fn) { + // return [1.3, 1.2] + return apply(data, dim, x => fn(x, prob, sorted)) + } + function _quantileSeqProbNumber (data, probOrN, sorted) { let probArr, one const dataArr = data.valueOf() diff --git a/test/unit-tests/function/statistics/quantileSeq.test.js b/test/unit-tests/function/statistics/quantileSeq.test.js index 7d89d0bed0..f64680e6e0 100644 --- a/test/unit-tests/function/statistics/quantileSeq.test.js +++ b/test/unit-tests/function/statistics/quantileSeq.test.js @@ -23,6 +23,40 @@ describe('quantileSeq', function () { assert.strictEqual(quantileSeq(lst, 1), 3.7) }) + it('should return the quantileSeq from a multidimensional array in the specified dimension', function () { + const arr = [ + [3.7, 2.7, 3.3, 1.3, 2.2, 3.1], + [3.8, 2.5, 3.2, 1.2, 3.2, 4.1] + ] + assert.deepStrictEqual(quantileSeq(arr, 0, 1), [1.3, 1.2]) + assert.deepStrictEqual(quantileSeq(arr, 0.25, 1), [2.325, 2.675]) + assert.deepStrictEqual(quantileSeq(arr, 0.5, 1), [2.9000000000000004, 3.2]) + assert.deepStrictEqual(quantileSeq(arr, 0.75, 1), [3.2499999999999996, 3.6499999999999995]) + assert.deepStrictEqual(quantileSeq(arr, 1, 1), [3.7, 4.1]) + assert.deepStrictEqual(quantileSeq(arr, 0, false, 1), [1.3, 1.2]) + assert.deepStrictEqual(quantileSeq(arr, 0.25, false, 1), [2.325, 2.675]) + assert.deepStrictEqual(quantileSeq(arr, 0.5, false, 1), [2.9000000000000004, 3.2]) + assert.deepStrictEqual(quantileSeq(arr, 0.75, false, 1), [3.2499999999999996, 3.6499999999999995]) + assert.deepStrictEqual(quantileSeq(arr, 1, false, 1), [3.7, 4.1]) + }) + + it('should return the quantileSeq from a multidimensional array in the specified dimension in the parser', function () { + const arr = [ + [3.7, 2.7, 3.3, 1.3, 2.2, 3.1], + [3.8, 2.5, 3.2, 1.2, 3.2, 4.1] + ] + assert.deepStrictEqual(math.evaluate('quantileSeq(arr, 0, 2)', { arr }), [1.3, 1.2]) + assert.deepStrictEqual(quantileSeq(arr, 0.25, 1), [2.325, 2.675]) + assert.deepStrictEqual(quantileSeq(arr, 0.5, 1), [2.9000000000000004, 3.2]) + assert.deepStrictEqual(quantileSeq(arr, 0.75, 1), [3.2499999999999996, 3.6499999999999995]) + assert.deepStrictEqual(quantileSeq(arr, 1, 1), [3.7, 4.1]) + assert.deepStrictEqual(quantileSeq(arr, 0, false, 1), [1.3, 1.2]) + assert.deepStrictEqual(quantileSeq(arr, 0.25, false, 1), [2.325, 2.675]) + assert.deepStrictEqual(quantileSeq(arr, 0.5, false, 1), [2.9000000000000004, 3.2]) + assert.deepStrictEqual(quantileSeq(arr, 0.75, false, 1), [3.2499999999999996, 3.6499999999999995]) + assert.deepStrictEqual(quantileSeq(arr, 1, false, 1), [3.7, 4.1]) + }) + it('should return the quantileSeq from an ascending array with number probability', function () { const lst = [1.3, 2.2, 2.7, 3.1, 3.3, 3.7] assert.strictEqual(quantileSeq(lst, 0, true), 1.3) @@ -78,7 +112,7 @@ describe('quantileSeq', function () { assert.deepStrictEqual(quantileSeq([math.unit('5mm'), math.unit('15mm'), math.unit('10mm')], 0.5), math.unit('10mm')) }) - it('should return the quantileSeq from an 1d matrix', function () { + it('should return the quantileSeq from a 1d matrix', function () { assert.strictEqual(quantileSeq(math.matrix([2, 4, 6, 8, 10, 12, 14]), 0.25), 5) }) From 8b4d23bc6513d3ea1e8f05c2881d0f136b1b9b1d Mon Sep 17 00:00:00 2001 From: David Contreras Date: Sat, 16 Sep 2023 17:05:59 +0000 Subject: [PATCH 22/25] Merge logic on _quantileSeqProbNumber --- .../transform/quantileSeq.transform.js | 6 +- src/function/statistics/quantileSeq.js | 96 +++++++------------ 2 files changed, 36 insertions(+), 66 deletions(-) diff --git a/src/expression/transform/quantileSeq.transform.js b/src/expression/transform/quantileSeq.transform.js index aa83d759fa..7d71f76f4d 100644 --- a/src/expression/transform/quantileSeq.transform.js +++ b/src/expression/transform/quantileSeq.transform.js @@ -3,7 +3,7 @@ import { createQuantileSeq } from '../../function/statistics/quantileSeq.js' import { lastDimToZeroBase } from './utils/lastDimToZeroBase.js' const name = 'quantileSeq' -const dependencies = ['typed', 'add', 'multiply', 'partitionSelect', 'compare', 'isInteger'] +const dependencies = ['typed', 'bignumber', 'add', 'divide', 'multiply', 'partitionSelect', 'compare', 'isInteger', 'smaller', 'smallerEq', 'larger'] /** * Attach a transform function to math.quantileSeq @@ -12,8 +12,8 @@ const dependencies = ['typed', 'add', 'multiply', 'partitionSelect', 'compare', * This transform changed the `dim` parameter of function std * from one-based to zero based */ -export const createQuantileSeqTransform = /* #__PURE__ */ factory(name, dependencies, ({ typed, add, multiply, partitionSelect, compare, isInteger }) => { - const quantileSeq = createQuantileSeq({ typed, add, multiply, partitionSelect, compare, isInteger }) +export const createQuantileSeqTransform = /* #__PURE__ */ factory(name, dependencies, ({ typed, bignumber, add, divide, multiply, partitionSelect, compare, isInteger, smaller, smallerEq, larger }) => { + const quantileSeq = createQuantileSeq({ typed, bignumber, add, divide, multiply, partitionSelect, compare, isInteger, smaller, smallerEq, larger }) return typed('quantileSeq', { 'Array|Matrix, number|BigNumber|Array, number': (arr, prob, dim) => quantileSeq(arr, prob, dimToZeroBase(dim)), diff --git a/src/function/statistics/quantileSeq.js b/src/function/statistics/quantileSeq.js index 3117073cc6..d4e6a9cc90 100644 --- a/src/function/statistics/quantileSeq.js +++ b/src/function/statistics/quantileSeq.js @@ -4,9 +4,9 @@ import { factory } from '../../utils/factory.js' import { createApply } from '../matrix/apply.js' const name = 'quantileSeq' -const dependencies = ['typed', 'add', 'multiply', 'partitionSelect', 'compare', 'isInteger'] +const dependencies = ['typed', 'bignumber', 'add', 'divide', 'multiply', 'partitionSelect', 'compare', 'isInteger', 'smaller', 'smallerEq', 'larger'] -export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ typed, add, multiply, partitionSelect, compare, isInteger }) => { +export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ typed, bignumber, add, divide, multiply, partitionSelect, compare, isInteger, smaller, smallerEq, larger }) => { /** * Compute the prob order quantile of a matrix or a list with values. * The sequence is sorted and the middle value is returned. @@ -45,10 +45,10 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ const apply = createApply({ typed, isInteger }) return typed(name, { - 'Array | Matrix, number | BigNumber | Unit': (data, p) => _quantileSeqProbNumber(data, p, false), - 'Array | Matrix, number | BigNumber | Unit, number': (data, prob, dim) => _quantileSeqDim(data, prob, false, dim, _quantileSeqProbNumber), - 'Array | Matrix, number | BigNumber | Unit, boolean': _quantileSeqProbNumber, - 'Array | Matrix, number | BigNumber | Unit, boolean, number': (data, prob, sorted, dim) => _quantileSeqDim(data, prob, sorted, dim, _quantileSeqProbNumber), + 'Array | Matrix, number | BigNumber': (data, p) => _quantileSeqProbNumber(data, p, false), + 'Array | Matrix, number | BigNumber, number': (data, prob, dim) => _quantileSeqDim(data, prob, false, dim, _quantileSeqProbNumber), + 'Array | Matrix, number | BigNumber, boolean': _quantileSeqProbNumber, + 'Array | Matrix, number | BigNumber, boolean, number': (data, prob, sorted, dim) => _quantileSeqDim(data, prob, sorted, dim, _quantileSeqProbNumber), 'Array | Matrix, Array | Matrix': (data, p) => _quantileSeqProbCollection(data, p, false), 'Array | Matrix, Array | Matrix, number': (data, prob, dim) => _quantileSeqDim(data, prob, false, dim, _quantileSeqProbCollection), 'Array | Matrix, Array | Matrix, boolean': _quantileSeqProbCollection, @@ -56,72 +56,42 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ }) function _quantileSeqDim (data, prob, sorted, dim, fn) { - // return [1.3, 1.2] return apply(data, dim, x => fn(x, prob, sorted)) } function _quantileSeqProbNumber (data, probOrN, sorted) { - let probArr, one + let probArr const dataArr = data.valueOf() - if (isNumber(probOrN)) { - if (probOrN < 0) { - throw new Error('N/prob must be non-negative') - } - - if (probOrN <= 1) { - // quantileSeq([a, b, c, d, ...], prob[,sorted]) - return _quantileSeq(dataArr, probOrN, sorted) - } - - if (probOrN > 1) { - // quantileSeq([a, b, c, d, ...], N[,sorted]) - if (!isInteger(probOrN)) { - throw new Error('N must be a positive integer') - } - - const nPlusOne = probOrN + 1 - probArr = new Array(probOrN) - for (let i = 0; i < probOrN; i++) { - probArr[i] = _quantileSeq(dataArr, (i + 1) / nPlusOne, sorted) - } - return probArr - } + if (smaller(probOrN, 0)) { + throw new Error('N/prob must be non-negative') } - - if (isBigNumber(probOrN)) { - const BigNumber = probOrN.constructor - - if (probOrN.isNegative()) { - throw new Error('N/prob must be non-negative') + if (smallerEq(probOrN, 1)) { + // quantileSeq([a, b, c, d, ...], prob[,sorted]) + return isNumber(probOrN) + ? _quantileSeq(dataArr, probOrN, sorted) + : bignumber(_quantileSeq(dataArr, probOrN, sorted)) + } + if (larger(probOrN, 1)) { + // quantileSeq([a, b, c, d, ...], N[,sorted]) + if (!isInteger(probOrN)) { + throw new Error('N must be a positive integer') } - one = new BigNumber(1) - - if (probOrN.lte(one)) { - // quantileSeq([a, b, c, d, ...], prob[,sorted]) - return new BigNumber(_quantileSeq(dataArr, probOrN, sorted)) + // largest possible Array length is 2^32-1 + // 2^32 < 10^15, thus safe conversion guaranteed + if (larger(probOrN, 4294967295)) { + throw new Error('N must be less than or equal to 2^32-1, as that is the maximum length of an Array') } - if (probOrN.gt(one)) { - // quantileSeq([a, b, c, d, ...], N[,sorted]) - if (!probOrN.isInteger()) { - throw new Error('N must be a positive integer') - } - - // largest possible Array length is 2^32-1 - // 2^32 < 10^15, thus safe conversion guaranteed - const intN = probOrN.toNumber() - if (intN > 4294967295) { - throw new Error('N must be less than or equal to 2^32-1, as that is the maximum length of an Array') - } + const nPlusOne = add(probOrN, 1) + probArr = [] - const nPlusOne = new BigNumber(intN + 1) - probArr = new Array(intN) - for (let i = 0; i < intN; i++) { - probArr[i] = new BigNumber(_quantileSeq(dataArr, new BigNumber(i + 1).div(nPlusOne), sorted)) - } - return probArr + for (let i = 0; smaller(i, probOrN); i++) { + const prob = divide(i + 1, nPlusOne) + probArr.push(_quantileSeq(dataArr, prob, sorted)) } + + return isNumber(probOrN) ? probArr : bignumber(probArr) } } @@ -141,8 +111,8 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ // quantileSeq([a, b, c, d, ...], [prob1, prob2, ...][,sorted]) const probOrNArr = probOrN.valueOf() - const probArr = new Array(probOrNArr.length) - for (let i = 0; i < probArr.length; ++i) { + const probArr = [] + for (let i = 0; i < probOrNArr.length; ++i) { const currProb = probOrNArr[i] if (isNumber(currProb)) { if (currProb < 0 || currProb > 1) { @@ -155,7 +125,7 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ } } - probArr[i] = _quantileSeq(dataArr, currProb, sorted) + probArr.push(_quantileSeq(dataArr, currProb, sorted)) } return probArr } From 119d35337133dc6e1fb572dcb0e7338c119c352a Mon Sep 17 00:00:00 2001 From: David Contreras Date: Sat, 16 Sep 2023 17:17:18 +0000 Subject: [PATCH 23/25] Reduced _quantileSeqProbCollection --- src/function/statistics/quantileSeq.js | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/src/function/statistics/quantileSeq.js b/src/function/statistics/quantileSeq.js index d4e6a9cc90..2352ed1e3b 100644 --- a/src/function/statistics/quantileSeq.js +++ b/src/function/statistics/quantileSeq.js @@ -1,4 +1,4 @@ -import { isBigNumber, isNumber } from '../../utils/is.js' +import { isNumber } from '../../utils/is.js' import { flatten } from '../../utils/array.js' import { factory } from '../../utils/factory.js' import { createApply } from '../matrix/apply.js' @@ -107,25 +107,11 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ function _quantileSeqProbCollection (data, probOrN, sorted) { const dataArr = data.valueOf() - let one - // quantileSeq([a, b, c, d, ...], [prob1, prob2, ...][,sorted]) const probOrNArr = probOrN.valueOf() const probArr = [] for (let i = 0; i < probOrNArr.length; ++i) { - const currProb = probOrNArr[i] - if (isNumber(currProb)) { - if (currProb < 0 || currProb > 1) { - throw new Error('Probability must be between 0 and 1, inclusive') - } - } else if (isBigNumber(currProb)) { - one = new currProb.constructor(1) - if (currProb.isNegative() || currProb.gt(one)) { - throw new Error('Probability must be between 0 and 1, inclusive') - } - } - - probArr.push(_quantileSeq(dataArr, currProb, sorted)) + probArr.push(_quantileSeq(dataArr, probOrNArr[i], sorted)) } return probArr } From 8c2d4d45f45a71985ad1a31c2aaa9b38898a0af9 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Sat, 16 Sep 2023 18:23:12 +0000 Subject: [PATCH 24/25] Merged logic of _quantileSeq --- .../transform/quantileSeq.transform.js | 6 +- src/function/statistics/quantileSeq.js | 71 +++++-------------- 2 files changed, 22 insertions(+), 55 deletions(-) diff --git a/src/expression/transform/quantileSeq.transform.js b/src/expression/transform/quantileSeq.transform.js index 7d71f76f4d..db93f67639 100644 --- a/src/expression/transform/quantileSeq.transform.js +++ b/src/expression/transform/quantileSeq.transform.js @@ -3,7 +3,7 @@ import { createQuantileSeq } from '../../function/statistics/quantileSeq.js' import { lastDimToZeroBase } from './utils/lastDimToZeroBase.js' const name = 'quantileSeq' -const dependencies = ['typed', 'bignumber', 'add', 'divide', 'multiply', 'partitionSelect', 'compare', 'isInteger', 'smaller', 'smallerEq', 'larger'] +const dependencies = ['typed', 'bignumber', 'add', 'subtract', 'divide', 'multiply', 'partitionSelect', 'compare', 'isInteger', 'smaller', 'smallerEq', 'larger'] /** * Attach a transform function to math.quantileSeq @@ -12,8 +12,8 @@ const dependencies = ['typed', 'bignumber', 'add', 'divide', 'multiply', 'partit * This transform changed the `dim` parameter of function std * from one-based to zero based */ -export const createQuantileSeqTransform = /* #__PURE__ */ factory(name, dependencies, ({ typed, bignumber, add, divide, multiply, partitionSelect, compare, isInteger, smaller, smallerEq, larger }) => { - const quantileSeq = createQuantileSeq({ typed, bignumber, add, divide, multiply, partitionSelect, compare, isInteger, smaller, smallerEq, larger }) +export const createQuantileSeqTransform = /* #__PURE__ */ factory(name, dependencies, ({ typed, bignumber, add, subtract, divide, multiply, partitionSelect, compare, isInteger, smaller, smallerEq, larger }) => { + const quantileSeq = createQuantileSeq({ typed, bignumber, add, subtract, divide, multiply, partitionSelect, compare, isInteger, smaller, smallerEq, larger }) return typed('quantileSeq', { 'Array|Matrix, number|BigNumber|Array, number': (arr, prob, dim) => quantileSeq(arr, prob, dimToZeroBase(dim)), diff --git a/src/function/statistics/quantileSeq.js b/src/function/statistics/quantileSeq.js index 2352ed1e3b..3e4cb9e13c 100644 --- a/src/function/statistics/quantileSeq.js +++ b/src/function/statistics/quantileSeq.js @@ -4,9 +4,9 @@ import { factory } from '../../utils/factory.js' import { createApply } from '../matrix/apply.js' const name = 'quantileSeq' -const dependencies = ['typed', 'bignumber', 'add', 'divide', 'multiply', 'partitionSelect', 'compare', 'isInteger', 'smaller', 'smallerEq', 'larger'] +const dependencies = ['typed', 'bignumber', 'add', 'subtract', 'divide', 'multiply', 'partitionSelect', 'compare', 'isInteger', 'smaller', 'smallerEq', 'larger'] -export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ typed, bignumber, add, divide, multiply, partitionSelect, compare, isInteger, smaller, smallerEq, larger }) => { +export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ typed, bignumber, add, subtract, divide, multiply, partitionSelect, compare, isInteger, smaller, smallerEq, larger }) => { /** * Compute the prob order quantile of a matrix or a list with values. * The sequence is sorted and the middle value is returned. @@ -132,68 +132,35 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ throw new Error('Cannot calculate quantile of an empty sequence') } - if (isNumber(prob)) { - const index = prob * (len - 1) - const fracPart = index % 1 - if (fracPart === 0) { - const value = sorted ? flat[index] : partitionSelect(flat, index) - return value - } - - const integerPart = Math.floor(index) - - let left - let right - if (sorted) { - left = flat[integerPart] - right = flat[integerPart + 1] - } else { - right = partitionSelect(flat, integerPart + 1) - - // max of partition is kth largest - left = flat[integerPart] - for (let i = 0; i < integerPart; ++i) { - if (compare(flat[i], left) > 0) { - left = flat[i] - } - } - } - - // Q(prob) = (1-f)*A[floor(index)] + f*A[floor(index)+1] - return add(multiply(left, 1 - fracPart), multiply(right, fracPart)) + const index = isNumber(prob) ? prob * (len - 1) : prob.times(len - 1) + const integerPart = isNumber(prob) ? Math.floor(index) : index.floor().toNumber() + const fracPart = isNumber(prob) ? index % 1 : index.minus(integerPart) + + if (isInteger(index)) { + return sorted + ? flat[index] + : partitionSelect( + flat, + isNumber(prob) ? index : index.valueOf() + ) } - - // If prob is a BigNumber - let index = prob.times(len - 1) - if (index.isInteger()) { - index = index.toNumber() - const value = sorted ? flat[index] : partitionSelect(flat, index) - return value - } - - const integerPart = index.floor() - const fracPart = index.minus(integerPart) - const integerPartNumber = integerPart.toNumber() - let left let right if (sorted) { - left = flat[integerPartNumber] - right = flat[integerPartNumber + 1] + left = flat[integerPart] + right = flat[integerPart + 1] } else { - right = partitionSelect(flat, integerPartNumber + 1) + right = partitionSelect(flat, integerPart + 1) // max of partition is kth largest - left = flat[integerPartNumber] - for (let i = 0; i < integerPartNumber; ++i) { + left = flat[integerPart] + for (let i = 0; i < integerPart; ++i) { if (compare(flat[i], left) > 0) { left = flat[i] } } } - // Q(prob) = (1-f)*A[floor(index)] + f*A[floor(index)+1] - const one = new fracPart.constructor(1) - return add(multiply(left, one.minus(fracPart)), multiply(right, fracPart)) + return add(multiply(left, subtract(1, fracPart)), multiply(right, fracPart)) } }) From 5d9f685f8d3e52be2b83ee6dd4678a8a9af0ec2a Mon Sep 17 00:00:00 2001 From: David Contreras Date: Sun, 17 Sep 2023 14:42:40 +0000 Subject: [PATCH 25/25] Fixed issue with transform and browser --- src/expression/transform/quantileSeq.transform.js | 10 ++++++++-- src/function/statistics/quantileSeq.js | 2 +- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/expression/transform/quantileSeq.transform.js b/src/expression/transform/quantileSeq.transform.js index db93f67639..97b967e255 100644 --- a/src/expression/transform/quantileSeq.transform.js +++ b/src/expression/transform/quantileSeq.transform.js @@ -16,8 +16,14 @@ export const createQuantileSeqTransform = /* #__PURE__ */ factory(name, dependen const quantileSeq = createQuantileSeq({ typed, bignumber, add, subtract, divide, multiply, partitionSelect, compare, isInteger, smaller, smallerEq, larger }) return typed('quantileSeq', { - 'Array|Matrix, number|BigNumber|Array, number': (arr, prob, dim) => quantileSeq(arr, prob, dimToZeroBase(dim)), - 'Array|Matrix, number|BigNumber|Array, boolean, number': (arr, prob, sorted, dim) => quantileSeq(arr, prob, sorted, dimToZeroBase(dim)) + 'Array | Matrix, number | BigNumber': quantileSeq, + 'Array | Matrix, number | BigNumber, number': (arr, prob, dim) => quantileSeq(arr, prob, dimToZeroBase(dim)), + 'Array | Matrix, number | BigNumber, boolean': quantileSeq, + 'Array | Matrix, number | BigNumber, boolean, number': (arr, prob, sorted, dim) => quantileSeq(arr, prob, sorted, dimToZeroBase(dim)), + 'Array | Matrix, Array | Matrix': quantileSeq, + 'Array | Matrix, Array | Matrix, number': (data, prob, dim) => quantileSeq(data, prob, dimToZeroBase(dim)), + 'Array | Matrix, Array | Matrix, boolean': quantileSeq, + 'Array | Matrix, Array | Matrix, boolean, number': (data, prob, sorted, dim) => quantileSeq(data, prob, sorted, dimToZeroBase(dim)) }) function dimToZeroBase (dim) { diff --git a/src/function/statistics/quantileSeq.js b/src/function/statistics/quantileSeq.js index 3e4cb9e13c..f647309ddd 100644 --- a/src/function/statistics/quantileSeq.js +++ b/src/function/statistics/quantileSeq.js @@ -4,7 +4,7 @@ import { factory } from '../../utils/factory.js' import { createApply } from '../matrix/apply.js' const name = 'quantileSeq' -const dependencies = ['typed', 'bignumber', 'add', 'subtract', 'divide', 'multiply', 'partitionSelect', 'compare', 'isInteger', 'smaller', 'smallerEq', 'larger'] +const dependencies = ['typed', '?bignumber', 'add', 'subtract', 'divide', 'multiply', 'partitionSelect', 'compare', 'isInteger', 'smaller', 'smallerEq', 'larger'] export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ typed, bignumber, add, subtract, divide, multiply, partitionSelect, compare, isInteger, smaller, smallerEq, larger }) => { /**