From dcf12c4cec1fa53d60f963ff20288283eb99e194 Mon Sep 17 00:00:00 2001 From: David Contreras Date: Wed, 20 Sep 2023 01:48:37 -0600 Subject: [PATCH] Quantile-refactor (#3003) * Included math to syntax when missing * Included solveODE * renamed initialStep as firstStep * Included tests for solveODE * Test the full state instead of the final state * Fixed issue with tolerance * Added unit signature for y0 * Included units test also for y0 * Included embedded docs and more tests * Included error for tspan * It works with bignumbers * reduced calling bignumber * extended the search for bignumbers * The jsdocs is less ambiguous * included tests for step options * Allowed for 0 minStep * Optimization to avoid checking the sign every step * refactor * Typo * removed unnecesary error * Fixes conflict with develop * Merge logic on _quantileSeqProbNumber * Reduced _quantileSeqProbCollection * Merged logic of _quantileSeq * Fixed issue with transform and browser --------- Co-authored-by: David Contreras Co-authored-by: Jos de Jong --- .../transform/quantileSeq.transform.js | 16 +- src/function/statistics/quantileSeq.js | 255 ++++++------------ 2 files changed, 88 insertions(+), 183 deletions(-) diff --git a/src/expression/transform/quantileSeq.transform.js b/src/expression/transform/quantileSeq.transform.js index aa83d759fa..97b967e255 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', 'subtract', 'divide', 'multiply', 'partitionSelect', 'compare', 'isInteger', 'smaller', 'smallerEq', 'larger'] /** * Attach a transform function to math.quantileSeq @@ -12,12 +12,18 @@ 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, 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)), - '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 6d772c8fbc..f647309ddd 100644 --- a/src/function/statistics/quantileSeq.js +++ b/src/function/statistics/quantileSeq.js @@ -1,12 +1,12 @@ -import { isBigNumber, isCollection, 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' const name = 'quantileSeq' -const dependencies = ['typed', 'add', 'multiply', 'partitionSelect', 'compare', 'isInteger'] +const dependencies = ['typed', '?bignumber', 'add', 'subtract', '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, 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. @@ -43,133 +43,77 @@ export const createQuantileSeq = /* #__PURE__ */ factory(name, dependencies, ({ */ const apply = createApply({ typed, isInteger }) - /** - * 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 typed(name, { - 'Array|Matrix, number|BigNumber|Array': (data, prob) => quantileSeq(data, prob, false), - 'Array|Matrix, number|BigNumber|Array, boolean': quantileSeq, - 'Array|Matrix, number|BigNumber|Array, number': (data, prob, dim) => _quantileSeqDim(data, prob, false, dim), - 'Array|Matrix, number|BigNumber|Array, boolean, number': (data, prob, sorted, dim) => _quantileSeqDim(data, prob, sorted, dim) + '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, + 'Array | Matrix, Array | Matrix, boolean, number': (data, prob, sorted, dim) => _quantileSeqDim(data, prob, sorted, dim, _quantileSeqProbCollection) }) - function _quantileSeqDim (data, prob, sorted, dim) { - // return [1.3, 1.2] - return apply(data, dim, x => quantileSeq(x, prob, sorted)) + function _quantileSeqDim (data, prob, sorted, dim, fn) { + return apply(data, dim, x => fn(x, prob, sorted)) } - 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') + function _quantileSeqProbNumber (data, probOrN, sorted) { + let probArr + const dataArr = data.valueOf() + if (smaller(probOrN, 0)) { + 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') + } - 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') - } - - 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 (isBigNumber(probOrN)) { - const BigNumber = probOrN.constructor - - if (probOrN.isNegative()) { - throw new Error('N/prob must be non-negative') - } - - one = new BigNumber(1) - - 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 (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 - } + // 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') + } - probArr[i] = _quantileSeq(dataArr, currProb, sorted) - } - return probArr - } + const nPlusOne = add(probOrN, 1) + probArr = [] - throw new TypeError('Unexpected type of argument in function quantileSeq') // FIXME: becomes redundant when converted to typed-function + for (let i = 0; smaller(i, probOrN); i++) { + const prob = divide(i + 1, nPlusOne) + probArr.push(_quantileSeq(dataArr, prob, sorted)) } - throw new TypeError('Unexpected type of argument in function quantileSeq') // FIXME: becomes redundant when converted to typed-function + return isNumber(probOrN) ? probArr : bignumber(probArr) } + } - throw new TypeError('Unexpected type of argument in function quantileSeq') // FIXME: becomes redundant when converted to typed-function + /** + * 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 _quantileSeqProbCollection (data, probOrN, sorted) { + const dataArr = data.valueOf() + // quantileSeq([a, b, c, d, ...], [prob1, prob2, ...][,sorted]) + const probOrNArr = probOrN.valueOf() + const probArr = [] + for (let i = 0; i < probOrNArr.length; ++i) { + probArr.push(_quantileSeq(dataArr, probOrNArr[i], sorted)) + } + return probArr } /** @@ -188,80 +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) - - validate(value) - - 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] - } - } - } - - 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)) - } - - // 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) - - validate(value) - - return value + 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() + ) } - - 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] } } } - - 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)) + return add(multiply(left, subtract(1, fracPart)), multiply(right, fracPart)) } })