From b97574db4f752f01636225ea86fbcbacacf7789f Mon Sep 17 00:00:00 2001 From: vmakarichev Date: Tue, 21 Jan 2025 18:29:03 +0200 Subject: [PATCH] Added methods to the lib --- libraries/diff-studio-utils/index.ts | 1 - libraries/diff-studio-utils/src/constants.ts | 52 + libraries/diff-studio-utils/src/math-tools.js | 219 --- libraries/diff-studio-utils/src/math-tools.ts | 251 ---- .../diff-studio-utils/src/model-error.ts | 19 + .../diff-studio-utils/src/scripting-tools.ts | 1310 +++++++++++++++++ .../solver-tools/callbacks/callback-base.js | 15 + .../solver-tools/callbacks/callback-base.ts | 8 + .../solver-tools/callbacks/callback-tools.ts | 20 + .../callbacks/iter-checker-callback.ts | 21 + .../callbacks/time-checker-callback.ts | 19 + .../src/solver-tools/lin-alg-tools.js | 61 + .../src/solver-tools/lin-alg-tools.ts | 75 + .../src/solver-tools/mrt-method.js | 201 +++ .../src/solver-tools/mrt-method.ts | 254 ++++ .../src/solver-tools/ros34prw-method.js | 256 ++++ .../src/solver-tools/ros34prw-method.ts | 329 +++++ .../src/solver-tools/ros3prw-method.js | 224 +++ .../src/solver-tools/ros3prw-method.ts | 289 ++++ .../src/solver-tools/solver-defs.js | 83 ++ .../src/solver-tools/solver-defs.ts | 95 ++ .../src/solver-tools/tests.js | 392 +++++ .../src/solver-tools/tests.ts | 437 ++++++ libraries/diff-studio-utils/src/solver.ts | 39 + libraries/diff-studio-utils/src/templates.ts | 159 ++ .../diff-studio-utils/src/ui-constants.ts | 326 ++++ libraries/diff-studio-utils/src/use-cases.ts | 608 ++++++++ libraries/diff-studio-utils/tsconfig.json | 2 +- .../diff-studio-utils/wasm/get-inv-matr.ts | 5 - .../diff-studio-utils/wasm/matr-api.d.ts | 1 - libraries/diff-studio-utils/wasm/matr-api.js | 8 - .../wasm/matrix-operations.js | 19 - .../wasm/matrix-operations.wasm | Bin 59309 -> 0 bytes 33 files changed, 5293 insertions(+), 505 deletions(-) create mode 100644 libraries/diff-studio-utils/src/constants.ts delete mode 100644 libraries/diff-studio-utils/src/math-tools.js delete mode 100644 libraries/diff-studio-utils/src/math-tools.ts create mode 100644 libraries/diff-studio-utils/src/model-error.ts create mode 100644 libraries/diff-studio-utils/src/scripting-tools.ts create mode 100644 libraries/diff-studio-utils/src/solver-tools/callbacks/callback-base.js create mode 100644 libraries/diff-studio-utils/src/solver-tools/callbacks/callback-base.ts create mode 100644 libraries/diff-studio-utils/src/solver-tools/callbacks/callback-tools.ts create mode 100644 libraries/diff-studio-utils/src/solver-tools/callbacks/iter-checker-callback.ts create mode 100644 libraries/diff-studio-utils/src/solver-tools/callbacks/time-checker-callback.ts create mode 100644 libraries/diff-studio-utils/src/solver-tools/lin-alg-tools.js create mode 100644 libraries/diff-studio-utils/src/solver-tools/lin-alg-tools.ts create mode 100644 libraries/diff-studio-utils/src/solver-tools/mrt-method.js create mode 100644 libraries/diff-studio-utils/src/solver-tools/mrt-method.ts create mode 100644 libraries/diff-studio-utils/src/solver-tools/ros34prw-method.js create mode 100644 libraries/diff-studio-utils/src/solver-tools/ros34prw-method.ts create mode 100644 libraries/diff-studio-utils/src/solver-tools/ros3prw-method.js create mode 100644 libraries/diff-studio-utils/src/solver-tools/ros3prw-method.ts create mode 100644 libraries/diff-studio-utils/src/solver-tools/solver-defs.js create mode 100644 libraries/diff-studio-utils/src/solver-tools/solver-defs.ts create mode 100644 libraries/diff-studio-utils/src/solver-tools/tests.js create mode 100644 libraries/diff-studio-utils/src/solver-tools/tests.ts create mode 100644 libraries/diff-studio-utils/src/solver.ts create mode 100644 libraries/diff-studio-utils/src/templates.ts create mode 100644 libraries/diff-studio-utils/src/ui-constants.ts create mode 100644 libraries/diff-studio-utils/src/use-cases.ts delete mode 100644 libraries/diff-studio-utils/wasm/get-inv-matr.ts delete mode 100644 libraries/diff-studio-utils/wasm/matr-api.d.ts delete mode 100644 libraries/diff-studio-utils/wasm/matr-api.js delete mode 100644 libraries/diff-studio-utils/wasm/matrix-operations.js delete mode 100644 libraries/diff-studio-utils/wasm/matrix-operations.wasm diff --git a/libraries/diff-studio-utils/index.ts b/libraries/diff-studio-utils/index.ts index a13c30a3c6..e69de29bb2 100644 --- a/libraries/diff-studio-utils/index.ts +++ b/libraries/diff-studio-utils/index.ts @@ -1 +0,0 @@ -export {invMatrix} from './wasm/get-inv-matr'; \ No newline at end of file diff --git a/libraries/diff-studio-utils/src/constants.ts b/libraries/diff-studio-utils/src/constants.ts new file mode 100644 index 0000000000..63493e1a7c --- /dev/null +++ b/libraries/diff-studio-utils/src/constants.ts @@ -0,0 +1,52 @@ +// Control constants +export const CONTROL_TAG = '#'; +export const CONTROL_TAG_LEN = CONTROL_TAG.length; +export const DF_NAME = 'df'; +const META = `${CONTROL_TAG}meta`; +export const MAX_LINE_CHART = 4; + +/** Control expressions for the problem specifying */ +export enum CONTROL_EXPR { + NAME = `${CONTROL_TAG}name`, + TAGS = `${CONTROL_TAG}tags`, + DESCR = `${CONTROL_TAG}description`, + DIF_EQ = `${CONTROL_TAG}equations`, + EXPR = `${CONTROL_TAG}expressions`, + ARG = `${CONTROL_TAG}argument`, + INITS = `${CONTROL_TAG}inits`, + CONSTS = `${CONTROL_TAG}constants`, + PARAMS = `${CONTROL_TAG}parameters`, + TOL = `${CONTROL_TAG}tolerance`, + LOOP = `${CONTROL_TAG}loop`, + UPDATE = `${CONTROL_TAG}update`, + RUN_ON_OPEN = `${META}.runOnOpen`, + RUN_ON_INPUT = `${META}.runOnInput`, + OUTPUT = `${CONTROL_TAG}output`, + COMMENT = `${CONTROL_TAG}comment`, + SOLVER = `${META}.solver`, + INPUTS = `${META}.inputs`, +}; + +/** Loop consts */ +export enum LOOP { + MIN_LINES_COUNT = 1, + COUNT_IDX = 0, + COUNT_NAME = '_count', + MIN_COUNT = 1, +}; + +/** UPDATE consts */ +export enum UPDATE { + MIN_LINES_COUNT = 1, + DURATION_IDX = 0, + DURATION = '_duration', +}; + +/** Ranges of the solver options */ +export const SOLVER_OPTIONS_RANGES = new Map([ + ['maxTime', {min: 1, max: 10000}], + ['scale', {min: 0.5, max: 1}], +]); + +export const TINY = 0.0001; +export const STEP_RATIO = 0.5; diff --git a/libraries/diff-studio-utils/src/math-tools.js b/libraries/diff-studio-utils/src/math-tools.js deleted file mode 100644 index dc578d13ad..0000000000 --- a/libraries/diff-studio-utils/src/math-tools.js +++ /dev/null @@ -1,219 +0,0 @@ -/*function luDecomp(A: Float64Array, n: number) { - const L = new Float64Array(n * n).fill(0); - const U = new Float64Array(n * n).fill(0); - - for (let i = 0; i < n; ++i) - L[i * n + i] = 1; - - let sumU = 0; - let sumL = 0; - let k = 0; - let j = 0; - let p = 0; - let i =0; - - for (k = 0; k < n; ++k) { - for (j = k; j < n; ++j) { - sumU = 0; - - for (p = 0; p < k; ++p) - sumU += L[k * n + p] * U[p * n + j]; - - U[k * n + j] = A[k * n +j] - sumU; - } - - for (i = k + 1; i < n; ++i) { - sumL = 0; - - for (p = 0; p < k; ++p) - sumL += L[i * n + p] * U[p * n + k]; - - L[i * n + k] = (A[i * n + k] - sumL) / U[k * n + k]; - } - } - - return {L: L, U: U}; -} - -function forwSubs(L: Float64Array, b: Float64Array, n: number) { - const y = new Float64Array(n); - let sumLy = 0; - - for (let i = 0; i < n; ++i) { - sumLy = 0; - - for (let j = 0; j < i; ++j) - sumLy += L[i * n + j] * y[j]; - - y[i] = (b[i] - sumLy) / L[i * n + i]; - } - - return y; -} - -function backSubs(U: Float64Array, y: Float64Array, n: number) { - const x = new Float64Array(n); - let sumUx = 0; - - for (let i = n - 1; i > -1; --i) { - sumUx = 0; - - for (let j = i + 1; j < n; ++j) - sumUx += U[i * n + j] * x[j]; - - x[i] = (y[i] - sumUx) / U[i * n + i]; - } - - return x; -}*/ -function luDecomp(A, L, U, n) { - L.fill(0); - U.fill(0); - for (var i_1 = 0; i_1 < n; ++i_1) - L[i_1 * n + i_1] = 1; - var sumU = 0; - var sumL = 0; - var k = 0; - var j = 0; - var p = 0; - var i = 0; - for (k = 0; k < n; ++k) { - for (j = k; j < n; ++j) { - sumU = 0; - for (p = 0; p < k; ++p) - sumU += L[k * n + p] * U[p * n + j]; - U[k * n + j] = A[k * n + j] - sumU; - } - for (i = k + 1; i < n; ++i) { - sumL = 0; - for (p = 0; p < k; ++p) - sumL += L[i * n + p] * U[p * n + k]; - L[i * n + k] = (A[i * n + k] - sumL) / U[k * n + k]; - } - } -} -function forwSubs(L, b, y, n) { - var sumLy = 0; - for (var i = 0; i < n; ++i) { - sumLy = 0; - for (var j = 0; j < i; ++j) - sumLy += L[i * n + j] * y[j]; - y[i] = (b[i] - sumLy) / L[i * n + i]; - } -} -function backSubs(U, y, x, n) { - var sumUx = 0; - for (var i = n - 1; i > -1; --i) { - sumUx = 0; - for (var j = i + 1; j < n; ++j) - sumUx += U[i * n + j] * x[j]; - x[i] = (y[i] - sumUx) / U[i * n + i]; - } -} -function solve(L, U, b, y, x, n) { - var sumLy = 0; - for (var i = 0; i < n; ++i) { - sumLy = 0; - for (var j = 0; j < i; ++j) - sumLy += L[i * n + j] * y[j]; - y[i] = (b[i] - sumLy) / L[i * n + i]; - } - var sumUx = 0; - for (var i = n - 1; i > -1; --i) { - sumUx = 0; - for (var j = i + 1; j < n; ++j) - sumUx += U[i * n + j] * x[j]; - x[i] = (y[i] - sumUx) / U[i * n + i]; - } -} -var A = new Float64Array([ - 2, 3, 5, - 4, -5, 2, - 1, 2, 3, -]); -var b = new Float64Array([ - 10, - 1, - 6, -]); -var L = new Float64Array(3 * 3); -var U = new Float64Array(3 * 3); -var x = new Float64Array(3); -var y = new Float64Array(3); -luDecomp(A, L, U, 3); -solve(L, U, b, y, x, 3); -console.log("Hello!"); -console.log(x); -/* - - -def backward_substitution(U, y): - """ - Solves Ux = y for x using backward substitution. - """ - n = U.shape[0] - x = np.zeros(n) - - for i in range(n-1, -1, -1): - sum_ux = sum(U[i][j] * x[j] for j in range(i+1, n)) - x[i] = (y[i] - sum_ux) / U[i][i] - - return x - -def matrix_inverse_lu(A): - """ - Computes the inverse of matrix A using LU decomposition. - - Parameters: - A (numpy.ndarray): Input square matrix - - Returns: - A_inv (numpy.ndarray): Inverse of matrix A - - Raises: - ValueError: If matrix is not square or is singular - """ - n = A.shape[0] - if A.shape[1] != n: - raise ValueError("Matrix must be square") - - # Get LU decomposition - L, U = lu_decomposition(A) - - # Initialize inverse matrix - A_inv = np.zeros((n, n)) - - # Compute inverse column by column - for j in range(n): - # Create unit vector ej - ej = np.zeros(n) - ej[j] = 1.0 - - # Solve Ly = ej for y - y = forward_substitution(L, ej) - - # Solve Ux = y for x (jth column of inverse) - x = backward_substitution(U, y) - - # Store the column in the inverse matrix - A_inv[:, j] = x - - return A_inv - -def verify_inverse(A, A_inv, tolerance=1e-10): - """ - Verifies if computed inverse is correct by checking if A * A_inv ≈ I. - - Parameters: - A (numpy.ndarray): Original matrix - A_inv (numpy.ndarray): Computed inverse matrix - tolerance (float): Maximum allowed deviation from identity matrix - - Returns: - bool: True if inverse is correct within tolerance - """ - n = A.shape[0] - I = np.eye(n) - product = np.dot(A, A_inv) - - return np.allclose(product, I, atol=tolerance)*/ diff --git a/libraries/diff-studio-utils/src/math-tools.ts b/libraries/diff-studio-utils/src/math-tools.ts deleted file mode 100644 index 113c8159b1..0000000000 --- a/libraries/diff-studio-utils/src/math-tools.ts +++ /dev/null @@ -1,251 +0,0 @@ -/*function luDecomp(A: Float64Array, n: number) { - const L = new Float64Array(n * n).fill(0); - const U = new Float64Array(n * n).fill(0); - - for (let i = 0; i < n; ++i) - L[i * n + i] = 1; - - let sumU = 0; - let sumL = 0; - let k = 0; - let j = 0; - let p = 0; - let i =0; - - for (k = 0; k < n; ++k) { - for (j = k; j < n; ++j) { - sumU = 0; - - for (p = 0; p < k; ++p) - sumU += L[k * n + p] * U[p * n + j]; - - U[k * n + j] = A[k * n +j] - sumU; - } - - for (i = k + 1; i < n; ++i) { - sumL = 0; - - for (p = 0; p < k; ++p) - sumL += L[i * n + p] * U[p * n + k]; - - L[i * n + k] = (A[i * n + k] - sumL) / U[k * n + k]; - } - } - - return {L: L, U: U}; -} - -function forwSubs(L: Float64Array, b: Float64Array, n: number) { - const y = new Float64Array(n); - let sumLy = 0; - - for (let i = 0; i < n; ++i) { - sumLy = 0; - - for (let j = 0; j < i; ++j) - sumLy += L[i * n + j] * y[j]; - - y[i] = (b[i] - sumLy) / L[i * n + i]; - } - - return y; -} - -function backSubs(U: Float64Array, y: Float64Array, n: number) { - const x = new Float64Array(n); - let sumUx = 0; - - for (let i = n - 1; i > -1; --i) { - sumUx = 0; - - for (let j = i + 1; j < n; ++j) - sumUx += U[i * n + j] * x[j]; - - x[i] = (y[i] - sumUx) / U[i * n + i]; - } - - return x; -}*/ - -function luDecomp(A: Float64Array, L: Float64Array, U: Float64Array, n: number) { - L.fill(0); - U.fill(0); - - for (let i = 0; i < n; ++i) - L[i * n + i] = 1; - - let sumU = 0; - let sumL = 0; - let k = 0; - let j = 0; - let p = 0; - let i =0; - - for (k = 0; k < n; ++k) { - for (j = k; j < n; ++j) { - sumU = 0; - - for (p = 0; p < k; ++p) - sumU += L[k * n + p] * U[p * n + j]; - - U[k * n + j] = A[k * n +j] - sumU; - } - - for (i = k + 1; i < n; ++i) { - sumL = 0; - - for (p = 0; p < k; ++p) - sumL += L[i * n + p] * U[p * n + k]; - - L[i * n + k] = (A[i * n + k] - sumL) / U[k * n + k]; - } - } -} - -function forwSubs(L: Float64Array, b: Float64Array, y: Float64Array, n: number) { - let sumLy = 0; - - for (let i = 0; i < n; ++i) { - sumLy = 0; - - for (let j = 0; j < i; ++j) - sumLy += L[i * n + j] * y[j]; - - y[i] = (b[i] - sumLy) / L[i * n + i]; - } -} - -function backSubs(U: Float64Array, y: Float64Array, x: Float64Array, n: number) { - let sumUx = 0; - - for (let i = n - 1; i > -1; --i) { - sumUx = 0; - - for (let j = i + 1; j < n; ++j) - sumUx += U[i * n + j] * x[j]; - - x[i] = (y[i] - sumUx) / U[i * n + i]; - } -} - -function solve(L: Float64Array, U: Float64Array, b: Float64Array, y: Float64Array, x: Float64Array, n: number) { - let sumLy = 0; - - for (let i = 0; i < n; ++i) { - sumLy = 0; - - for (let j = 0; j < i; ++j) - sumLy += L[i * n + j] * y[j]; - - y[i] = (b[i] - sumLy) / L[i * n + i]; - } - - let sumUx = 0; - - for (let i = n - 1; i > -1; --i) { - sumUx = 0; - - for (let j = i + 1; j < n; ++j) - sumUx += U[i * n + j] * x[j]; - - x[i] = (y[i] - sumUx) / U[i * n + i]; - } -} - -let A = new Float64Array([ - 2, 3, 5, - 4, -5, 2, - 1, 2, 3, -]); - -let b = new Float64Array([ - 10, - 1, - 6, -]); - -let L = new Float64Array(3 * 3); -let U = new Float64Array(3 * 3); - -let x = new Float64Array(3); -let y = new Float64Array(3); - -luDecomp(A, L, U, 3); -solve(L, U, b, y, x, 3); - -console.log("Hello!"); -console.log(x); - -/* - - -def backward_substitution(U, y): - """ - Solves Ux = y for x using backward substitution. - """ - n = U.shape[0] - x = np.zeros(n) - - for i in range(n-1, -1, -1): - sum_ux = sum(U[i][j] * x[j] for j in range(i+1, n)) - x[i] = (y[i] - sum_ux) / U[i][i] - - return x - -def matrix_inverse_lu(A): - """ - Computes the inverse of matrix A using LU decomposition. - - Parameters: - A (numpy.ndarray): Input square matrix - - Returns: - A_inv (numpy.ndarray): Inverse of matrix A - - Raises: - ValueError: If matrix is not square or is singular - """ - n = A.shape[0] - if A.shape[1] != n: - raise ValueError("Matrix must be square") - - # Get LU decomposition - L, U = lu_decomposition(A) - - # Initialize inverse matrix - A_inv = np.zeros((n, n)) - - # Compute inverse column by column - for j in range(n): - # Create unit vector ej - ej = np.zeros(n) - ej[j] = 1.0 - - # Solve Ly = ej for y - y = forward_substitution(L, ej) - - # Solve Ux = y for x (jth column of inverse) - x = backward_substitution(U, y) - - # Store the column in the inverse matrix - A_inv[:, j] = x - - return A_inv - -def verify_inverse(A, A_inv, tolerance=1e-10): - """ - Verifies if computed inverse is correct by checking if A * A_inv ≈ I. - - Parameters: - A (numpy.ndarray): Original matrix - A_inv (numpy.ndarray): Computed inverse matrix - tolerance (float): Maximum allowed deviation from identity matrix - - Returns: - bool: True if inverse is correct within tolerance - """ - n = A.shape[0] - I = np.eye(n) - product = np.dot(A, A_inv) - - return np.allclose(product, I, atol=tolerance)*/ diff --git a/libraries/diff-studio-utils/src/model-error.ts b/libraries/diff-studio-utils/src/model-error.ts new file mode 100644 index 0000000000..534b0b9a41 --- /dev/null +++ b/libraries/diff-studio-utils/src/model-error.ts @@ -0,0 +1,19 @@ +/** Diff Studio model error */ +export class ModelError extends Error { + private helpUrl: string; + private toHighlight: string | undefined = undefined; + + constructor(message: string, helpUrl: string, toHighlight?: string) { + super(message); + this.helpUrl = helpUrl; + this.toHighlight = toHighlight; + } + + public getHelpUrl() { + return this.helpUrl; + } + + public getToHighlight() { + return this.toHighlight; + } +}; // ModelError \ No newline at end of file diff --git a/libraries/diff-studio-utils/src/scripting-tools.ts b/libraries/diff-studio-utils/src/scripting-tools.ts new file mode 100644 index 0000000000..8177b147e0 --- /dev/null +++ b/libraries/diff-studio-utils/src/scripting-tools.ts @@ -0,0 +1,1310 @@ +/* eslint-disable max-len */ +/* Scripting tools for the Initial Value Problem (IVP) solver: + - parser of formulas defining IVP; + - JS-script generator. + + The parser processes IVP formulas given in the special form (see "Project structure" in README.md). + + JS-script generator creates DATAGROK JavaScript script: annotation & code. +*/ + +import {CONTROL_TAG, CONTROL_TAG_LEN, DF_NAME, CONTROL_EXPR, LOOP, UPDATE, MAX_LINE_CHART, + SOLVER_OPTIONS_RANGES, TINY, STEP_RATIO} from './constants'; + +import {ModelError} from './model-error'; + +// Scripting specific constants +export const CONTROL_SEP = ':'; +const COMMA = ','; +const EQUAL_SIGN = '='; +const DIV_SIGN = '/'; +const SERVICE = '_'; +export const BRACE_OPEN = '{'; +export const BRACE_CLOSE = '}'; +export const BRACKET_OPEN = '['; +export const BRACKET_CLOSE = ']'; +export const ANNOT_SEPAR = ';'; +const DEFAULT_TOL = '0.00005'; +export const DEFAULT_SOLVER_SETTINGS: string = '{}'; +const COLUMNS = `${SERVICE}columns`; +const COMMENT_SEQ = '//'; +export const STAGE_COL_NAME = `${SERVICE}Stage`; +const INCEPTION = 'Inception'; + +/** Solver package name */ +const PACKAGE_NAME = 'DiffStudio'; + +/** Numerical solver function */ +const SOLVER_FUNC = 'solveEquations'; + +/** Elementary math tools */ +const MATH_FUNCS = ['pow', 'sin', 'cos', 'tan', 'asin', 'acos', 'atan', 'sqrt', 'exp', 'log', 'sinh', 'cosh', 'tanh']; +const POW_IDX = MATH_FUNCS.indexOf('pow'); +const MATH_CONSTS = ['PI', 'E']; + +/** Default meta */ +const defaultMetas = `//meta.runOnOpen: true +//meta.runOnInput: true +//meta.features: {"sens-analysis": true, "fitting": true}`; + +/** Numerical input specification */ +export type Input = { + value: number, + annot: string | null, +}; + +/** Argument of IVP specification */ +type Arg = { + name: string, + initial: Input, + final: Input, + step: Input, + updateName: string | undefined, +}; + +/** Input keys of Arg */ +export const ARG_INPUT_KEYS = ['initial', 'final', 'step']; + +/** Scripting specific constants */ +export enum SCRIPTING { + ARG_NAME = 'name', + COUNT = 'count', + DURATION = 'duration', +}; + +/** Differential equations specification */ +type DifEqs = { + equations: Map, + solutionNames: string[] +}; + +/** Loop specification */ +type Loop = { + count: Input, + updates: string[], +}; + +/** Update specification */ +type Update = { + name: string, + durationFormula: string, + updates: string[], +}; + +/** Output specification */ +type Output = { + caption: string, + formula: string | null, +}; + +/** Initial Value Problem (IVP) specification type */ +export type IVP = { + name: string, + tags: string | null, + descr: string | null, + deqs: DifEqs, + exprs: Map | null, + arg: Arg, + inits: Map, + consts: Map | null, + params: Map | null, + tolerance: string, + usedMathFuncs: number[], + usedMathConsts: number[], + loop: Loop | null, + updates: Update[] | null, + metas: string[], + outputs: Map | null, + solverSettings: string, + inputsLookup: string | null, +}; + +/** Help links for model errors */ +enum ERROR_LINK { + MAIN_DOCS = '/help/compute/diff-studio', + CORE_BLOCKS = `${MAIN_DOCS}#core-blocks`, + COMPS_SYNTAX = `${MAIN_DOCS}#model-components-and-syntax`, + LOOP = `${MAIN_DOCS}#cyclic-processes`, + UPDATE = `${MAIN_DOCS}#multistage-processes`, + UI_OPTS = `${MAIN_DOCS}#user-interface-options`, + LOOP_VS_UPDATE = `${MAIN_DOCS}#advanced-features`, + SOLVER_CONFIG = `${MAIN_DOCS}#solver-configuration`, + MODEL_PARAMS = `${MAIN_DOCS}#model-parameters`, +}; + +/** Specific error messages */ +enum ERROR_MSG { + CTRL_EXPR = `Unsupported control expression with the tag **"${CONTROL_TAG}"**`, + ARG = `'The **${CONTROL_EXPR.ARG}** block must consist of 3 lines specifying initial and final time, and solution grid step.`, + LOOP = `The **${CONTROL_EXPR.LOOP}** block must contain at least one line.`, + COUNT = 'Incorrect loop count', + LOOP_VS_UPDATE = `The **${CONTROL_EXPR.LOOP}** and **${CONTROL_EXPR.UPDATE}** blocks cannot be used simultaneously.`, + UPDATE_LINES_COUNT = `The **${CONTROL_EXPR.UPDATE}** block must contain at least one line.`, + DURATION = 'Incorrect update duration', + BRACES = ' Missing one of the braces (**{**, **}**).', + COLON = `Incorrect position of **"${CONTROL_SEP}"**.`, + CASE_INSENS = 'Non-unique name (case-insensitive): use different caption for ', + MISSING_INIT = `Correct the **${CONTROL_EXPR.INITS}** block.`, + UNDEF_NAME = `Model name missing. Specify the model name in the **${CONTROL_EXPR.NAME}** block.`, + UNDEF_DEQS = `Differential equation(s) are required for this model. Add equation(s) under the **${CONTROL_EXPR.DIF_EQ}** block.`, + UNDEF_INITS = `Initial conditions are required for this model. Add initial conditions under the **${CONTROL_EXPR.INITS}** block.`, + UNDEF_ARG = `Argument specification is required for this model. Specify an argument, its range, and a grid step in the **${CONTROL_EXPR.ARG}** block.`, + CORRECT_ARG_LIM = `Correct limits in the **${CONTROL_EXPR.ARG}** block.`, + INTERVAL = `Incorrect range for`, + NEGATIVE_STEP = `Solution grid step must be positive. Correct the **${CONTROL_EXPR.ARG}** block.`, + INCOR_STEP = `Grid step must less than the length of solution interval. Correct the **${CONTROL_EXPR.ARG}** block.`, + MISS_COLON = `Missing **"${CONTROL_SEP}"**`, + NAN = `is not a valid number. Correct the line`, + SERVICE_START = `Variable names must not begin with **"${SERVICE}"**.`, + REUSE_NAME = 'Variable reuse (case-insensitive): rename ', + SOLVER = `Incorrect solver options. Correct the **${CONTROL_EXPR.SOLVER}** line.`, +} // ERROR_MSG + +/** Datagrok annotations */ +enum ANNOT { + NAME = '//name:', + DESCR = '//description:', + TAGS = '//tags:', + LANG = '//language: javascript', + DOUBLE_INPUT = '//input: double', + INT_INPUT = '//input: int', + OUTPUT = `//output: dataframe ${DF_NAME}`, + EDITOR = '//editor: Compute:RichFunctionViewEditor', + SIDEBAR = '//sidebar: @compute', + CAPTION = 'caption:', + ARG_INIT = '{caption: Initial; category: Argument}', + ARG_FIN = '{caption: Final; category: Argument}', + ARG_STEP = '{caption: Step; category: Argument}', + INITS = 'category: Initial values', + PARAMS = 'category: Parameters', +} + +/** JS-scripting components */ +enum SCRIPT { + CONSTS = '// constants', + ODE_COM = '// the problem definition', + ODE = 'let odes = {', + SOLVER_COM = '// solve the problem', + SOLVER = `const solver = await grok.functions.eval('${PACKAGE_NAME}:${SOLVER_FUNC}');`, + PREPARE = 'let call = solver.prepare({problem: odes, options: opts});', + CALL = 'await call.call();', + OUTPUT = `let ${DF_NAME} = call.getParamValue('${DF_NAME}');`, + SPACE2 = ' ', + SPACE4 = ' ', + SPACE6 = ' ', + SPACE8 = ' ', + FUNC_VALS = '// extract function values', + EVAL_EXPR = '// evaluate expressions', + COMP_OUT = '// compute output', + MATH_FUNC_COM = '// used Math-functions', + MATH_CONST_COM = '// used Math-constants', + ONE_STAGE_COM = '\n// one stage solution', + ONE_STAGE_BEGIN = 'let _oneStage = async (', + ONE_STAGE_END = ') => {', + ASYNC_OUTPUT = `let ${DF_NAME} = await _oneStage(`, + RETURN_OUTPUT = `return call.getParamValue('${DF_NAME}');`, + EMPTY_OUTPUT = `let ${DF_NAME} = DG.DataFrame.create();`, + APPEND = `${DF_NAME}.append(`, + SOLUTION_DF_COM = '// solution dataframe', + LOOP_INTERVAL_COM = '// loop interval', + LOOP_INTERVAL = `${SERVICE}interval`, + LAST_IDX = `${SERVICE}lastIdx`, + UPDATE_COM = '// update ', + CUSTOM_OUTPUT_COM = '// create custom output', + CUSTOM_COLUMNS = `let ${COLUMNS} = [`, + ONE_STAGE = 'await _oneStage(', + SEGMENT_COM = '// add segment category', +} + +/** Limits of the problem specification */ +type Block = { + begin: number, + end: number +} + +/** Get start of the problem skipping note-lines */ +function getStartOfProblemDef(lines: string[]): number { + const linesCount = lines.length; + let idx = 0; + + while (!lines[idx].startsWith(CONTROL_TAG)) { + ++idx; + + if (idx === linesCount) + throw new ModelError(ERROR_MSG.UNDEF_NAME, ERROR_LINK.CORE_BLOCKS); + } + + return idx; +} + +/** Get limits of blocks specifying IVP */ +function getBlocks(lines: string[], start: number): Block[] { + const linesCount = lines.length; + let beg = start; + let idx = start + 1; + const blocks = [] as Block[]; + + while (true) { + if (idx === linesCount) { + blocks.push({begin: beg, end: idx}); + break; + } + + if (lines[idx].startsWith(CONTROL_TAG)) { + blocks.push({begin: beg, end: idx}); + beg = idx; + } + + ++idx; + } + + return blocks; +} + +/** Concatenate multi-line formulas & expressions */ +function concatMultilineFormulas(source: string[]): string[] { + const res: string[] = [source[0]]; + + let idxRes = 0; + let idxSource = 1; + const size = source.length; + + while (idxSource < size) { + const str = source[idxSource]; + + if (!str.includes(EQUAL_SIGN)) + res[idxRes] += str; + else { + res.push(str); + ++idxRes; + } + + ++idxSource; + } + + return res; +} + +/** Get indeces of math functions used in the current text */ +function getUsedMathIds(text: string, mathIds: string[]) { + const res = [] as number[]; + const size = mathIds.length; + + for (let i = 0; i < size; ++i) { + if (text.includes(`${mathIds[i]}`)) + res.push(i); + } + + return res; +} + +/** Get differential equations */ +function getDifEquations(lines: string[]): DifEqs { + const deqs = new Map(); + const names = [] as string[]; + + let divIdx = 0; + let eqIdx = 0; + + for (const line of lines) { + if (line === undefined) + throw new ModelError(ERROR_MSG.UNDEF_DEQS, ERROR_LINK.CORE_BLOCKS, CONTROL_EXPR.DIF_EQ); + + divIdx = line.indexOf(DIV_SIGN); + eqIdx = line.indexOf(EQUAL_SIGN); + const name = line.slice(line.indexOf('d') + 1, divIdx).replace(/[() ]/g, ''); + names.push(name); + deqs.set(name, line.slice(eqIdx + 1).trim()); + } + + return { + equations: deqs, + solutionNames: names, + }; +} + +/** Get expressions of IVP */ +function getExpressions(lines: string[]): Map { + const exprs = new Map(); + let eqIdx = 0; + + for (const line of lines) { + if (line === undefined) + break; + + eqIdx = line.indexOf(EQUAL_SIGN); + exprs.set(line.slice(0, eqIdx).replaceAll(' ', ''), line.slice(eqIdx + 1).trim()); + } + + return exprs; +} + +/** Get input specification */ +function getInput(line: string) : Input { + const str = line.slice(line.indexOf(EQUAL_SIGN) + 1).trim(); + const braceIdx = str.indexOf(BRACE_OPEN); + let val: number; + + if (braceIdx === -1) { // no annotation + val = Number(str); + + // Check right-hand side + if (isNaN(val)) + throw new ModelError(`**${str}** ${ERROR_MSG.NAN}\n **${line}**.`, ERROR_LINK.COMPS_SYNTAX, line); + + return { + value: val, + annot: null, + }; + } + + // There is an annotation + val = Number(str.slice(0, braceIdx)); + + // Check right-hand side + if (isNaN(val)) + throw new ModelError(`**${str.slice(0, braceIdx)}** ${ERROR_MSG.NAN}: **${line}**`, ERROR_LINK.COMPS_SYNTAX, line); + + return { + value: val, + annot: str.slice(braceIdx), + }; +} + +/** Get argument (independent variable) of IVP */ +function getArg(lines: string[]): Arg { + if (lines.length !== 4) + throw new ModelError(ERROR_MSG.ARG, ERROR_LINK.CORE_BLOCKS, CONTROL_EXPR.ARG); + + const sepIdx = lines[0].indexOf(CONTROL_SEP); + if (sepIdx === -1) + throw new ModelError(`${ERROR_MSG.MISS_COLON}. Correct the line **${lines[0]}**`, ERROR_LINK.CORE_BLOCKS, lines[0]); + + const commaIdx = lines[0].indexOf(COMMA); + + return { + name: ((commaIdx > -1) ? lines[0].slice(sepIdx + 1, commaIdx) : lines[0].slice(sepIdx + 1)).trim(), + initial: getInput(lines[1]), + final: getInput(lines[2]), + step: getInput(lines[3]), + updateName: (commaIdx > -1) ? lines[0].slice(commaIdx + 1).trim() : undefined, + }; +} + +/** Get equalities specification */ +function getEqualities(lines: string[], begin: number, end: number): Map { + const source = concatMultilineFormulas(lines.slice(begin, end)); + const eqs = new Map(); + let eqIdx = 0; + + for (const line of source) { + if (line === undefined) + break; + + eqIdx = line.indexOf(EQUAL_SIGN); + eqs.set(line.slice(0, eqIdx).replaceAll(' ', ''), getInput(line)); + } + + return eqs; +} + +/** Get loop specification */ +function getLoop(lines: string[], begin: number, end: number): Loop { + if (begin >= end) + throw new ModelError(ERROR_MSG.LOOP, ERROR_LINK.LOOP, CONTROL_EXPR.LOOP); + + const source = concatMultilineFormulas(lines.slice(begin, end)); + const size = source.length; + + if (size < LOOP.MIN_LINES_COUNT) + throw new ModelError(ERROR_MSG.LOOP, ERROR_LINK.LOOP, CONTROL_EXPR.LOOP); + + const count = getInput(source[LOOP.COUNT_IDX]); + if (count.value < LOOP.MIN_COUNT) + throw new ModelError(`${ERROR_MSG.COUNT}: **${count.value}**.`, ERROR_LINK.LOOP, source[LOOP.COUNT_IDX]); + + return {count: count, updates: source.slice(LOOP.COUNT_IDX + 1)}; +} + +/** Get update specification */ +function getUpdate(lines: string[], begin: number, end: number): Update { + const colonIdx = lines[begin].indexOf(CONTROL_SEP); + if (colonIdx === -1) { + throw new ModelError( + `${ERROR_MSG.MISS_COLON}. Correct the line: **${lines[begin]}**`, + ERROR_LINK.UPDATE, + lines[begin], + ); + } + + const source = concatMultilineFormulas(lines.slice(begin + 1, end)); + const size = source.length; + + if (size < UPDATE.MIN_LINES_COUNT) + throw new ModelError(ERROR_MSG.UPDATE_LINES_COUNT, ERROR_LINK.UPDATE, CONTROL_EXPR.UPDATE); + + if (source[UPDATE.DURATION_IDX] == undefined) + throw new ModelError(ERROR_MSG.UPDATE_LINES_COUNT, ERROR_LINK.UPDATE, CONTROL_EXPR.UPDATE); + + const eqIdx = source[UPDATE.DURATION_IDX].indexOf(EQUAL_SIGN); + + if (eqIdx === -1) { + throw new ModelError( + `Missing **${EQUAL_SIGN}**. Correct the line: **${source[UPDATE.DURATION_IDX]}**`, + ERROR_LINK.UPDATE, + source[UPDATE.DURATION_IDX], + ); + } + + return { + name: lines[begin].slice(colonIdx + 1), + durationFormula: source[UPDATE.DURATION_IDX].slice(eqIdx + 1).trim(), + updates: source.slice(UPDATE.DURATION_IDX + 1), + }; +} // getUpdate + +/** Get output specification */ +function getOutput(lines: string[], begin: number, end: number): Map { + const res = new Map(); + + let line: string; + let token: string; + let eqIdx: number; + let openIdx: number; + let closeIdx: number; + let colonIdx: number; + + for (let i = begin; i < end; ++i) { + line = lines[i].trim(); + + if (line.length < 1) + continue; + + eqIdx = line.indexOf(EQUAL_SIGN); + openIdx = line.indexOf(BRACE_OPEN); + closeIdx = line.indexOf(BRACE_CLOSE); + colonIdx = line.indexOf(CONTROL_SEP); + + // check braces + if (openIdx * closeIdx <= 0) { + throw new ModelError( + `${ERROR_MSG.BRACES} Correct the line **${line}** in the **${CONTROL_EXPR.OUTPUT}** block.`, + ERROR_LINK.UI_OPTS, + line, + ); + } + + // check ":" + if (openIdx * colonIdx <= 0) { + throw new ModelError( + `${ERROR_MSG.COLON} Correct the line **${line}** in the **#output** block.`, + ERROR_LINK.UI_OPTS, + line, + ); + } + + if (eqIdx === -1) {// no formula + if (openIdx === -1) { // no annotation + token = line.trim(); + res.set(token, { + caption: token, + formula: null, + }); + } else { // there is an annotation + token = line.slice(0, openIdx).trim(); + res.set(token, { + caption: line.slice(colonIdx + 1, closeIdx).trim(), + formula: null, + }); + } + } else // there is a formula + if (openIdx === -1) { // no annotation + token = line.slice(0, eqIdx).trim(); + res.set(token, { + caption: token, + formula: line.slice(eqIdx + 1), + }); + } else { // there is an annotation + token = line.slice(0, eqIdx).trim(); + res.set(token, { + caption: line.slice(colonIdx + 1, closeIdx), + formula: line.slice(eqIdx + 1, openIdx), + }); + } + } // for i + + return res; +} // getOutput + +/** Get initial value problem specification given in the text */ +export function getIVP(text: string): IVP { + // The current Initial Value Problem (IVP) specification + let name: string; + let tags: string | null = null; + let descr: string | null = null; + let deqs: DifEqs; + let exprs: Map | null = null; + let arg: Arg; + let inits: Map; + let consts: Map | null = null; + let params: Map | null = null; + let tolerance = DEFAULT_TOL; + let loop: Loop | null = null; + const updates = [] as Update[]; + const metas = [] as string[]; + let outputs: Map | null = null; + let solverSettings = DEFAULT_SOLVER_SETTINGS; + let inputsLookup: string | null = null; + + // 0. Split text into lines & remove comments + const lines = text.replaceAll('\t', ' ').split('\n') + .map((s) => { + const idx = s.indexOf(COMMENT_SEQ); + return s.slice(0, (idx >= 0) ? idx : undefined).trimStart(); + }) + .filter((s) => s !== ''); + + // 1. Skip first lines without the control tag + const start = getStartOfProblemDef(lines); + + // 2. Get blocks limits + const blocks = getBlocks(lines, start); + + // 3. Process blocks + for (const block of blocks) { + const firstLine = lines[block.begin]; + + if (firstLine.startsWith(CONTROL_EXPR.NAME)) { // the 'name' block + name = firstLine.slice(firstLine.indexOf(CONTROL_SEP) + 1).trim(); + } else if (firstLine.startsWith(CONTROL_EXPR.TAGS)) { // the 'tags' block + tags = firstLine.slice(firstLine.indexOf(CONTROL_SEP) + 1).trim(); + } else if (firstLine.startsWith(CONTROL_EXPR.DESCR)) { // the 'description' block + descr = firstLine.slice(firstLine.indexOf(CONTROL_SEP) + 1).trim(); + } else if (firstLine.startsWith(CONTROL_EXPR.DIF_EQ)) { // the 'differential equations' block + deqs = getDifEquations(concatMultilineFormulas(lines.slice(block.begin + 1, block.end))); + } else if (firstLine.startsWith(CONTROL_EXPR.EXPR)) { // the 'expressions' block + exprs = getExpressions(concatMultilineFormulas(lines.slice(block.begin + 1, block.end))); + } else if (firstLine.startsWith(CONTROL_EXPR.ARG)) { // the 'argument' block + arg = getArg(concatMultilineFormulas(lines.slice(block.begin, block.end))); + } else if (firstLine.startsWith(CONTROL_EXPR.INITS)) { // the 'initial values' block + inits = getEqualities(lines, block.begin + 1, block.end); + } else if (firstLine.startsWith(CONTROL_EXPR.CONSTS)) { // the 'constants' block + consts = getEqualities(lines, block.begin + 1, block.end); + } else if (firstLine.startsWith(CONTROL_EXPR.PARAMS)) { // the 'parameters' block + params = getEqualities(lines, block.begin + 1, block.end); + } else if (firstLine.startsWith(CONTROL_EXPR.TOL)) { // the 'tolerance' block + tolerance = firstLine.slice( firstLine.indexOf(CONTROL_SEP) + 1).trim(); + } else if (firstLine.startsWith(CONTROL_EXPR.SOLVER)) { // the 'solver settings' block + solverSettings = firstLine.slice(firstLine.indexOf(CONTROL_SEP) + 1).trim(); + } else if (firstLine.startsWith(CONTROL_EXPR.LOOP)) { // the 'loop' block + loop = getLoop(lines, block.begin + 1, block.end); + } else if (firstLine.startsWith(CONTROL_EXPR.UPDATE)) { // the 'update' block + updates.push(getUpdate(lines, block.begin, block.end)); + } else if (firstLine.startsWith(CONTROL_EXPR.OUTPUT)) { // the 'output' block + outputs = getOutput(lines, block.begin + 1, block.end); + } else if (firstLine.startsWith(CONTROL_EXPR.INPUTS)) { // the 'inputs' block + inputsLookup = firstLine.slice(firstLine.indexOf(CONTROL_SEP) + 1).trim(); + } else if (firstLine.startsWith(CONTROL_EXPR.COMMENT)) { // the 'comment' block + // just skip it + } else + metas.push(firstLine.slice(CONTROL_TAG_LEN)); + } // for + + // check loop- & update-features: just one is supported + if ( (loop !== null) && (updates.length > 0)) + throw new ModelError(ERROR_MSG.LOOP_VS_UPDATE, ERROR_LINK.LOOP_VS_UPDATE, CONTROL_EXPR.LOOP); + + // obtained model + const ivp = { + name: name!, + tags: tags, + descr: descr, + deqs: deqs!, + exprs: exprs, + arg: arg!, + inits: inits!, + consts: consts, + params: params, + tolerance: tolerance, + usedMathFuncs: getUsedMathIds(text, MATH_FUNCS), + usedMathConsts: getUsedMathIds(text, MATH_CONSTS), + loop: loop, + updates: (updates.length === 0) ? null : updates, + metas: metas, + outputs: outputs, + solverSettings: solverSettings, + inputsLookup: inputsLookup, + }; + + checkCorrectness(ivp); + + return ivp; +} // getIVP + +/** Return input specification, required for annotation generating */ +function getInputSpec(inp: Input): string { + if (inp.annot) + return `${inp.value} ${inp.annot}`; + + return `${inp.value}`; +} + +/** Return annotation line specifying viewers */ +function getViewersLine(ivp: IVP): string { + const outputColsCount = (ivp.outputs) ? ivp.outputs.size : ivp.inits.size; + const multiAxis = (outputColsCount > MAX_LINE_CHART - 1) ? 'true' : 'false'; + + const segments = (ivp.updates) ? ` segmentColumnName: "${STAGE_COL_NAME}",` : ''; + + // eslint-disable-next-line max-len + return `viewer: Line chart(block: 100, multiAxis: "${multiAxis}",${segments} multiAxisLegendPosition: "RightCenter", autoLayout: "false", showAggrSelectors: "false") | Grid(block: 100)`; +} + +/** Generate annotation lines */ +function getAnnot(ivp: IVP, toAddViewers = true, toAddEditor = false): string[] { + const res = [] as string[]; + + // the 'name' line + res.push(`${ANNOT.NAME} ${ivp.name}`); + + // the 'tags' line + if (ivp.tags) + res.push(`${ANNOT.TAGS} ${ivp.tags}`); + + // the 'description' line + if (ivp.descr) + res.push(`${ANNOT.DESCR} ${ivp.descr}`); + + // the 'language' line + res.push(ANNOT.LANG); + + // the 'loop' lines + if (ivp.loop) + res.push(`${ANNOT.INT_INPUT} ${LOOP.COUNT_NAME} = ${getInputSpec(ivp.loop.count)}`); + + // argument lines + const arg = ivp.arg; + const t0 = `${SERVICE}${arg.name}0`; + const t1 = `${SERVICE}${arg.name}1`; + const h = `${SERVICE}h`; + res.push(`${ANNOT.DOUBLE_INPUT} ${t0} = ${getInputSpec(arg.initial)}`); + res.push(`${ANNOT.DOUBLE_INPUT} ${t1} = ${getInputSpec(arg.final)}`); + res.push(`${ANNOT.DOUBLE_INPUT} ${h} = ${getInputSpec(arg.step)}`); + + // initial values lines + ivp.inits.forEach((val, key) => res.push(`${ANNOT.DOUBLE_INPUT} ${key} = ${getInputSpec(val)}`)); + + // parameters lines + if (ivp.params !== null) + ivp.params.forEach((val, key) => res.push(`${ANNOT.DOUBLE_INPUT} ${key} = ${getInputSpec(val)}`)); + + // the 'output' line + if (toAddViewers) + res.push(`${ANNOT.OUTPUT} {${ANNOT.CAPTION} ${ivp.name}; ${getViewersLine(ivp)}}`); + else + res.push(`${ANNOT.OUTPUT} {${ANNOT.CAPTION} ${ivp.name}`); + + // the 'editor' line + if (toAddEditor) { + res.push(ANNOT.EDITOR); + res.push(ANNOT.SIDEBAR); + } + + if (ivp.metas.length >0) + ivp.metas.forEach((line) => res.push(`//${line}`)); + else + res.push(defaultMetas); + + return res; +} // getAnnot + +/** Returns math functions arguments string */ +function getMathArg(funcIdx: number): string { + return (funcIdx > POW_IDX) ? '(x)' : '(x, y)'; +} + +/** Return custom output lines: no expressions */ +function getCustomOutputLinesNoExpressions(name: string, + outputs: Map, toAddUpdateCol: boolean): string[] { + const res = ['']; + + res.push(SCRIPT.CUSTOM_OUTPUT_COM); + res.push(SCRIPT.CUSTOM_COLUMNS); + + outputs.forEach((val, key) => { + if (!val.formula) + res.push(`${SCRIPT.SPACE2}${DF_NAME}.col('${key}'),`); + }); + + if (toAddUpdateCol) + res.push(`${SCRIPT.SPACE2}${DF_NAME}.col('${STAGE_COL_NAME}'),`); + + res.push('];'); + + outputs.forEach((val, key) => { + if (!val.formula) + res.push(`${DF_NAME}.col('${key}').name = '${val.caption}';`); + }); + + res.push(`${DF_NAME} = DG.DataFrame.fromColumns(${COLUMNS});`); + res.push(`${DF_NAME}.name = '${name}';`); + return res; +} // getCustomOutputLinesNoExpressions + +/** Return custom output lines: with expressions */ +function getCustomOutputLinesWithExpressions(ivp: IVP): string[] { + const res = ['']; + + res.push(SCRIPT.CUSTOM_OUTPUT_COM); + + // 1. Solution raw data + res.push(`const ${ivp.arg.name}RawData = ${DF_NAME}.col('${ivp.arg.name}').getRawData();`); + res.push(`let ${ivp.arg.name};`); + res.push(`const len = ${DF_NAME}.rowCount;\n`); + + ivp.inits.forEach((val, key) => { + res.push(`const ${key}RawData = ${DF_NAME}.col('${key}').getRawData();`); + }); + + res.push(''); + + // 2. Expressions raw data & variables + ivp.exprs!.forEach((val, key) => { + if (ivp.outputs!.has(key)) + res.push(`const ${key}RawData = new Float64Array(len);`); + + res.push(`let ${key};`); + }); + + res.push(''); + + // 3. Computations + res.push('for (let i = 0; i < len; ++i) {'); + res.push(`${SCRIPT.SPACE2}${ivp.arg.name} = ${ivp.arg.name}RawData[i];`); + + ivp.inits.forEach((val, key) => { + res.push(`${SCRIPT.SPACE2}${key} = ${key}RawData[i];`); + }); + + ivp.exprs!.forEach((val, key) => { + res.push(`${SCRIPT.SPACE2}${key} = ${val};`); + + if (ivp.outputs!.has(key)) + res.push(`${SCRIPT.SPACE2}${key}RawData[i] = ${key};`); + }); + + res.push('}\n'); + + // 4. Form output + res.push(`${DF_NAME} = DG.DataFrame.fromColumns([`); + ivp.outputs!.forEach((val, key) => { + if (!val.formula) + res.push(`${SCRIPT.SPACE2}DG.Column.fromFloat64Array('${val.caption}', ${key}RawData.slice(0, len)),`); + }); + + if (ivp.updates !== null) + res.push(`${SCRIPT.SPACE2}${DF_NAME}.col('${STAGE_COL_NAME}'),`); + + res.push(']);'); + res.push(`${DF_NAME}.name = '${ivp.name}';`); + + return res; +} // getCustomOutputLinesWithExpressions + +/** Return custom output lines */ +function getCustomOutputLines(ivp: IVP): string[] { + if (ivp.exprs !== null) { + for (const key of ivp.exprs.keys()) { + if (ivp.outputs?.has(key)) + return getCustomOutputLinesWithExpressions(ivp); + } + } + + return getCustomOutputLinesNoExpressions(ivp.name, ivp.outputs!, (ivp.updates !== null)); +} // getCustomOutputLinesWithExpressions + +/** Return main body of JS-script: basic variant */ +function getScriptMainBodyBasic(ivp: IVP): string[] { + const res = [] as string[]; + + // 1. Constants lines + if (ivp.consts !== null) { + res.push(''); + res.push(SCRIPT.CONSTS); + ivp.consts.forEach((val, key) => res.push(`const ${key} = ${val.value};`)); + } + + // 2. The problem definition lines + res.push(''); + res.push(SCRIPT.ODE_COM); + res.push(SCRIPT.ODE); + res.push(`${SCRIPT.SPACE4}name: '${ivp.name}',`); + + // 2.1) argument + const t = ivp.arg.name; + const t0 = `${SERVICE}${t}0`; + const t1 = `${SERVICE}${t}1`; + const h = `${SERVICE}h`; + res.push(`${SCRIPT.SPACE4}arg: {name: '${t}', start: ${t0}, finish: ${t1}, step: ${h}},`); + + const names = ivp.deqs.solutionNames; + + // 2.2) initial values + res.push(`${SCRIPT.SPACE4}initial: [${names.join(', ')}],`); + + // 2.3) the right-hand side of the problem + res.push(`${SCRIPT.SPACE4}func: (${t}, ${SERVICE}y, ${SERVICE}output) => {`); + + res.push(`${SCRIPT.SPACE6}${SCRIPT.FUNC_VALS}`); + names.forEach((name, idx) => res.push(`${SCRIPT.SPACE6}const ${name} = ${SERVICE}y[${idx}];`)); + + if (ivp.exprs !== null) { + res.push(`\n${SCRIPT.SPACE6}${SCRIPT.EVAL_EXPR}`); + ivp.exprs.forEach((val, key) => res.push(`${SCRIPT.SPACE6}const ${key} = ${val};`)); + } + + res.push(`\n${SCRIPT.SPACE6}${SCRIPT.COMP_OUT}`); + names.forEach((name, idx) => res.push(`${SCRIPT.SPACE6}${SERVICE}output[${idx}] = ${ivp.deqs.equations.get(name)};`)); + + res.push(`${SCRIPT.SPACE4}},`); + + // 2.4) final lines of the problem specification + res.push(`${SCRIPT.SPACE4}tolerance: ${ivp.tolerance},`); + //res.push(`${SCRIPT.SPACE6}solverOptions: ${ivp.solverSettings.replaceAll(ANNOT_SEPAR, COMMA)},`); + res.push(`${SCRIPT.SPACE4}solutionColNames: [${names.map((key) => `'${key}'`).join(', ')}]`); + res.push('};'); + + // 2.5) solver options + res.push(''); + res.push(`let opts = ${ivp.solverSettings.replaceAll(';', ',')};`); + + // 3. Math functions + if (ivp.usedMathFuncs.length > 0) { + res.push(''); + res.push(SCRIPT.MATH_FUNC_COM); + ivp.usedMathFuncs.forEach((i) => + res.push(`const ${MATH_FUNCS[i]} = ${getMathArg(i)} => Math.${MATH_FUNCS[i]}${getMathArg(i)};`, + )); + } + + // 4. Math constants + if (ivp.usedMathConsts.length > 0) { + res.push(''); + res.push(SCRIPT.MATH_CONST_COM); + ivp.usedMathConsts.forEach((i) => res.push(`const ${MATH_CONSTS[i]} = Math.${MATH_CONSTS[i]};`)); + } + + // 5. The 'call solver' lines + res.push(''); + res.push(SCRIPT.SOLVER_COM); + res.push(SCRIPT.SOLVER); + res.push(SCRIPT.PREPARE); + res.push(SCRIPT.CALL); + res.push(SCRIPT.OUTPUT); + + return res; +} // getScriptMainBodyBasic + +/** Return function for JS-script */ +function getScriptFunc(ivp: IVP, funcParamsNames: string): string[] { + const res = [] as string[]; + + // 0. Function declaration + res.push(SCRIPT.ONE_STAGE_COM); + res.push(`${SCRIPT.ONE_STAGE_BEGIN}${funcParamsNames}${SCRIPT.ONE_STAGE_END}`); + + // 1. Constants lines + if (ivp.consts !== null) { + res.push(''); + res.push(SCRIPT.SPACE2 + SCRIPT.CONSTS); + ivp.consts.forEach((val, key) => res.push(`${SCRIPT.SPACE2}const ${key} = ${val.value};`)); + } + + // 2. The problem definition lines + res.push(SCRIPT.SPACE2 + SCRIPT.ODE_COM); + res.push(SCRIPT.SPACE2 + SCRIPT.ODE); + res.push(`${SCRIPT.SPACE6}name: '${ivp.name}',`); + + // 2.1) argument + const t = ivp.arg.name; + const t0 = `${SERVICE}${t}0`; + const t1 = `${SERVICE}${t}1`; + const h = `${SERVICE}h`; + res.push(`${SCRIPT.SPACE6}arg: {name: '${t}', start: ${t0}, finish: ${t1}, step: ${h}},`); + + const names = ivp.deqs.solutionNames; + + // 2.2) initial values + res.push(`${SCRIPT.SPACE6}initial: [${names.join(', ')}],`); + + // 2.3) the right-hand side of the problem + res.push(`${SCRIPT.SPACE6}func: (${t}, ${SERVICE}y, ${SERVICE}output) => {`); + + res.push(`${SCRIPT.SPACE8}${SCRIPT.FUNC_VALS}`); + names.forEach((name, idx) => res.push(`${SCRIPT.SPACE8}const ${name} = ${SERVICE}y[${idx}];`)); + + if (ivp.exprs !== null) { + res.push(`\n${SCRIPT.SPACE8}${SCRIPT.EVAL_EXPR}`); + ivp.exprs.forEach((val, key) => res.push(`${SCRIPT.SPACE8}const ${key} = ${val};`)); + } + + res.push(`\n${SCRIPT.SPACE8}${SCRIPT.COMP_OUT}`); + names.forEach((name, idx) => res.push(`${SCRIPT.SPACE8}${SERVICE}output[${idx}] = ${ivp.deqs.equations.get(name)};`)); + + res.push(`${SCRIPT.SPACE6}},`); + + // 2.4) final lines of the problem specification + res.push(`${SCRIPT.SPACE6}tolerance: ${ivp.tolerance},`); + //res.push(`${SCRIPT.SPACE6}solverOptions: ${ivp.solverSettings.replaceAll(ANNOT_SEPAR, COMMA)},`); + res.push(`${SCRIPT.SPACE6}solutionColNames: [${names.map((key) => `'${key}'`).join(', ')}]`); + res.push(`${SCRIPT.SPACE2}};`); + + // 2.5) solver options + res.push(''); + res.push(`${SCRIPT.SPACE2}let opts = ${ivp.solverSettings.replaceAll(';', ',')};`); + + // 3. Math functions + if (ivp.usedMathFuncs.length > 0) { + res.push(''); + res.push(SCRIPT.SPACE2 + SCRIPT.MATH_FUNC_COM); + // eslint-disable-next-line max-len + ivp.usedMathFuncs.forEach((i) => res.push(`${SCRIPT.SPACE2}const ${MATH_FUNCS[i]} = ${getMathArg(i)} => Math.${MATH_FUNCS[i]}${getMathArg(i)};`)); + } + + // 4. Math constants + if (ivp.usedMathConsts.length > 0) { + res.push(''); + res.push(SCRIPT.SPACE2 + SCRIPT.MATH_CONST_COM); + ivp.usedMathConsts.forEach((i) => res.push(`${SCRIPT.SPACE2}const ${MATH_CONSTS[i]} = Math.${MATH_CONSTS[i]};`)); + } + + // 5. The 'call solver' lines + res.push(''); + res.push(SCRIPT.SPACE2 + SCRIPT.SOLVER_COM); + res.push(SCRIPT.SPACE2 + SCRIPT.SOLVER); + res.push(SCRIPT.SPACE2 + SCRIPT.PREPARE); + res.push(SCRIPT.SPACE2 + SCRIPT.CALL); + res.push(SCRIPT.SPACE2 + SCRIPT.RETURN_OUTPUT); + + // 6. Close the function + res.push('};'); + + return res; +} // getScriptFunc + +/** Return main body of JS-script: loop case */ +function getScriptMainBodyLoopCase(ivp: IVP): string[] { + const funcParamsNames = getFuncParamsNames(ivp); + const res = getScriptFunc(ivp, funcParamsNames); + res.push(''); + //res.push(`${SCRIPT.ASYNC_OUTPUT}${funcParamsNames});`); + + res.push(SCRIPT.SOLUTION_DF_COM); + const dfNames = getSolutionDfColsNames(ivp); + + res.push(`let ${DF_NAME} = DG.DataFrame.fromColumns([`); + dfNames.forEach((name) => res.push(`${SCRIPT.SPACE2}DG.Column.fromFloat64Array('${name}', []),`)); + res.push(`]);`); + res.push(`${DF_NAME}.name = '${ivp.name}';`); + res.push(''); + + res.push(SCRIPT.LOOP_INTERVAL_COM); + res.push(`const ${SCRIPT.LOOP_INTERVAL} = ${SERVICE}${ivp.arg.name}1 - ${SERVICE}${ivp.arg.name}0;`); + res.push(''); + res.push(`let ${SCRIPT.LAST_IDX} = 0;\n`); + + res.push(SCRIPT.SOLVER_COM); + res.push(`for (let ${SERVICE}idx = 0; ${SERVICE}idx < ${LOOP.COUNT_NAME}; ++${SERVICE}idx) {`); + ivp.loop!.updates.forEach((upd) => res.push(`${SCRIPT.SPACE2}${upd};`)); + res.push(`${SCRIPT.SPACE2}${SCRIPT.APPEND}${SCRIPT.ONE_STAGE}${funcParamsNames}), true);`); + res.push(`${SCRIPT.SPACE2}${SERVICE}${ivp.arg.name}0 = ${SERVICE}${ivp.arg.name}1;`); + res.push(`${SCRIPT.SPACE2}${SERVICE}${ivp.arg.name}1 += ${SCRIPT.LOOP_INTERVAL};`); + res.push(`${SCRIPT.SPACE2}${SCRIPT.LAST_IDX} = ${DF_NAME}.rowCount - 1;`); + + dfNames.forEach((name, idx) => { + if (idx !== 0) + res.push(`${SCRIPT.SPACE2}${name} = ${DF_NAME}.get('${name}', ${SCRIPT.LAST_IDX});`); + }); + + // eslint-disable-next-line max-len + res.push(`${SCRIPT.SPACE2}${DF_NAME}.set('${ivp.arg.name}', ${SCRIPT.LAST_IDX}, ${SERVICE}${ivp.arg.name}0 - Math.min(${SERVICE}h * ${STEP_RATIO}, ${TINY}));`); + + res.push('};'); + + return res; +} // getScriptMainBodyLoopCase + +/** Return main body of JS-script: update case */ +function getScriptMainBodyUpdateCase(ivp: IVP): string[] { + const funcParamsNames = getFuncParamsNames(ivp); + const res = getScriptFunc(ivp, funcParamsNames); + + res.push(''); + + res.push(SCRIPT.SOLUTION_DF_COM); + const dfNames = getSolutionDfColsNames(ivp); + + res.push(`${SCRIPT.ASYNC_OUTPUT}${funcParamsNames});`); + + res.push(''); + + res.push(SCRIPT.SEGMENT_COM); + // eslint-disable-next-line max-len + res.push(`${DF_NAME}.columns.add(DG.Column.fromList('string', '${STAGE_COL_NAME}', new Array(${DF_NAME}.rowCount).fill('${ivp.arg.updateName ?? INCEPTION}')));`); + + res.push(''); + res.push(`let ${SCRIPT.LAST_IDX} = 0;`); + + ivp.updates!.forEach((upd, idx) => { + res.push(''); + res.push(`${SCRIPT.UPDATE_COM} ${idx + 1}`); + res.push(`const ${UPDATE.DURATION}${idx + 1} = ${upd.durationFormula};`); + res.push(`${SCRIPT.LAST_IDX} = ${DF_NAME}.rowCount - 1;`); + + // eslint-disable-next-line max-len + res.push(`${DF_NAME}.set('${ivp.arg.name}', ${SCRIPT.LAST_IDX}, ${SERVICE}${ivp.arg.name}1 - Math.min(${SERVICE}h * ${STEP_RATIO}, ${TINY}));`); + + dfNames.forEach((name, idx) => { + if (idx !== 0) + res.push(`${name} = ${DF_NAME}.get('${name}', ${SCRIPT.LAST_IDX});`); + }); + + upd.updates.forEach((upd) => res.push(`${upd};`)); + + res.push(`${SERVICE}${ivp.arg.name}0 = ${SERVICE}${ivp.arg.name}1;`); + res.push(`${SERVICE}${ivp.arg.name}1 += ${UPDATE.DURATION}${idx + 1};`); + + res.push(`const ${SERVICE}DF${idx + 1} = ${SCRIPT.ONE_STAGE}${funcParamsNames});`); + // eslint-disable-next-line max-len + res.push(`${SERVICE}DF${idx + 1}.columns.add(DG.Column.fromList('string', '${STAGE_COL_NAME}', new Array(${SERVICE}DF${idx + 1}.rowCount).fill('${upd.name}')));`); + res.push(`${SCRIPT.APPEND}${SERVICE}DF${idx + 1}, true);`); + }); + + return res; +} // getScriptMainBodyUpdateCase + +/** Return main body of JS-script */ +function getScriptMainBody(ivp: IVP): string[] { + if (ivp.loop) + return getScriptMainBodyLoopCase(ivp); + + if (ivp.updates) + return getScriptMainBodyUpdateCase(ivp); + + return getScriptMainBodyBasic(ivp); +} + +/** Return JS-script lines */ +export function getScriptLines(ivp: IVP, toAddViewers = true, toAddEditor = false): string[] { + const res = getAnnot(ivp, toAddViewers, toAddEditor).concat(getScriptMainBody(ivp)); + + if (ivp.outputs) + return res.concat(getCustomOutputLines(ivp)); + + return res; +} + +/** Return parameters of JS-script */ +export function getScriptParams(ivp: IVP): Record { + const res = {} as Record; + + if (ivp.loop) + res[LOOP.COUNT_NAME] = ivp.loop.count.value; + + const arg = ivp.arg; + + res[`${SERVICE}${arg.name}0`] = arg.initial.value; + res[`${SERVICE}${arg.name}1`] = arg.final.value; + res[`${SERVICE}h`] = arg.step.value; + + ivp.inits.forEach((val, key) => res[key] = val.value); + + if (ivp.params) + ivp.params.forEach((val, key) => res[key] = val.value); + + return res; +} + +/** Return func parameters names string */ +function getFuncParamsNames(ivp: IVP): string { + const names = [] as string []; + + const arg = ivp.arg.name; + + names.push(`${SERVICE}${arg}0`); + names.push(`${SERVICE}${arg}1`); + names.push(`${SERVICE}h`); + + ivp.inits.forEach((val, key) => names.push(key)); + + if (ivp.params) + ivp.params.forEach((val, key) => names.push(key)); + + return names.join(', '); +} + +/** Return solution dataframe columns names*/ +function getSolutionDfColsNames(ivp: IVP): string[] { + const res = [] as string[]; + + res.push(ivp.arg.name); + + ivp.inits.forEach((val, key) => res.push(key)); + + return res; +} + +/** Check solver settings */ +function checkSolverSettings(line: string): void { + const settings = new Map(); + + const openBraceIdx = line.indexOf(BRACE_OPEN); + const closeBraceIdx = line.indexOf(BRACE_CLOSE); + let sepIdx: number; + + if ((openBraceIdx < 0 )|| (closeBraceIdx < 0)) + throw new ModelError(`${ERROR_MSG.BRACES}. Correct the line **${line}**.`, ERROR_LINK.SOLVER_CONFIG, line); + + for (const item of line.slice(openBraceIdx + 1, closeBraceIdx).split(ANNOT_SEPAR)) { + sepIdx = item.indexOf(CONTROL_SEP); + + if (sepIdx > 1) + settings.set(item.slice(0, sepIdx).trim(), item.slice(sepIdx + 1).trim()); + } + + SOLVER_OPTIONS_RANGES.forEach((range, opt) => { + if (settings.has(opt)) { + const val = Number(settings.get(opt)); + + if ((val < range.min) || (val > range.max)) { + throw new ModelError( + `${ERROR_MSG.SOLVER}: **${opt}** must be in the range **${range.min}..${range.max}**.`, + ERROR_LINK.SOLVER_CONFIG, + line, + ); + } + } + }); +} // checkSolverSettings + +/** Check IVP correctness */ +function checkCorrectness(ivp: IVP): void { + // 0. Check basic elements + if (ivp.name === undefined) + throw new ModelError(ERROR_MSG.UNDEF_NAME, ERROR_LINK.CORE_BLOCKS); + + if ((ivp.deqs === undefined) || (ivp.deqs.equations.size === 0)) + throw new ModelError(ERROR_MSG.UNDEF_DEQS, ERROR_LINK.CORE_BLOCKS); + + if ((ivp.inits === undefined) || (ivp.inits.size === 0)) + throw new ModelError(ERROR_MSG.UNDEF_INITS, ERROR_LINK.CORE_BLOCKS); + + if (ivp.arg === undefined) + throw new ModelError(ERROR_MSG.UNDEF_ARG, ERROR_LINK.CORE_BLOCKS); + + // 1. Check initial values + ivp.deqs.equations.forEach((ignore, name) => { + if (!ivp.inits.has(name)) { + throw new ModelError( + `Initial value for **${name}** is missing. ${ERROR_MSG.MISSING_INIT}`, + ERROR_LINK.CORE_BLOCKS, + CONTROL_EXPR.INITS, + ); + } + }); + + // 2. Check names of output columns + const usedNames = [] as string[]; + let lowCase: string; + + if (ivp.outputs !== null) { + ivp.outputs.forEach((val) => { + lowCase = val.caption.toLowerCase(); + + if (usedNames.includes(lowCase)) + throw new ModelError(`${ERROR_MSG.CASE_INSENS}**${val.caption}**.`, ERROR_LINK.UI_OPTS, CONTROL_EXPR.OUTPUT); + else + usedNames.push(lowCase); + }); + } else { + const usedNames = [ivp.arg.name]; + + ivp.deqs.solutionNames.forEach((name) => { + lowCase = name.toLowerCase(); + + if (usedNames.includes(lowCase)) + throw new ModelError(`${ERROR_MSG.CASE_INSENS}**${name}**.`, ERROR_LINK.UI_OPTS); + else + usedNames.push(lowCase); + }); + } + + // 3. Check argument + if (ivp.arg.initial.value >= ivp.arg.final.value) { + throw new ModelError( + `${ERROR_MSG.INTERVAL} **${ivp.arg.name}**. ${ERROR_MSG.CORRECT_ARG_LIM}`, + ERROR_LINK.CORE_BLOCKS, + CONTROL_EXPR.ARG, + ); + } + + if (ivp.arg.step.value <= 0) { + throw new ModelError( + `${ERROR_MSG.NEGATIVE_STEP}`, + ERROR_LINK.CORE_BLOCKS); + } + + if (ivp.arg.step.value > ivp.arg.final.value - ivp.arg.initial.value) + throw new ModelError(ERROR_MSG.INCOR_STEP, ERROR_LINK.CORE_BLOCKS); + + // 4. Check script inputs, due to https://reddata.atlassian.net/browse/GROK-15152 + const scriptInputs = [] as string[]; + let current: string; + + // 4.1) initial values + ivp.inits.forEach((_, key) => { + if (key[0] === SERVICE) { + throw new ModelError( + `${ERROR_MSG.SERVICE_START} Correct **"${key}"** in the **${CONTROL_EXPR.INITS}** block.`, + ERROR_LINK.CORE_BLOCKS, + ); + } + + current = key.toLocaleLowerCase(); + + if (scriptInputs.includes(current)) { + throw new ModelError( + `${ERROR_MSG.REUSE_NAME} **"${key}"** in the **${CONTROL_EXPR.INITS}** block.`, + ERROR_LINK.CORE_BLOCKS, + ); + } + + scriptInputs.push(current); + }); + + // 4.2) parameters + if (ivp.params !== null) { + ivp.params.forEach((_, key) => { + if (key[0] === SERVICE) { + throw new ModelError( + `${ERROR_MSG.SERVICE_START}: correct **"${key}"** in the **${CONTROL_EXPR.PARAMS}** block.`, + ERROR_LINK.MODEL_PARAMS, + ); + } + + current = key.toLocaleLowerCase(); + + if (scriptInputs.includes(current)) { + throw new ModelError( + `${ERROR_MSG.REUSE_NAME} **"${key}"** in the **${CONTROL_EXPR.PARAMS}** block.`, + ERROR_LINK.MODEL_PARAMS, + ); + } + + scriptInputs.push(current); + }); + } + + // 5. Check solver settings + checkSolverSettings(ivp.solverSettings); +} // checkCorrectness diff --git a/libraries/diff-studio-utils/src/solver-tools/callbacks/callback-base.js b/libraries/diff-studio-utils/src/solver-tools/callbacks/callback-base.js new file mode 100644 index 0000000000..cdda3ed4a5 --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/callbacks/callback-base.js @@ -0,0 +1,15 @@ +"use strict"; +// Solver's callback base +Object.defineProperty(exports, "__esModule", { value: true }); +exports.Callback = void 0; +/** Solver callback */ +var Callback = /** @class */ (function () { + function Callback() { + } + ; + Callback.prototype.onIterationStart = function () { }; + Callback.prototype.onComputationsCompleted = function () { }; + return Callback; +}()); +exports.Callback = Callback; +; diff --git a/libraries/diff-studio-utils/src/solver-tools/callbacks/callback-base.ts b/libraries/diff-studio-utils/src/solver-tools/callbacks/callback-base.ts new file mode 100644 index 0000000000..11b7f58f54 --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/callbacks/callback-base.ts @@ -0,0 +1,8 @@ +// Solver's callback base + +/** Solver callback */ +export class Callback { + constructor() {}; + public onIterationStart(): void {} + public onComputationsCompleted(): void {} +}; diff --git a/libraries/diff-studio-utils/src/solver-tools/callbacks/callback-tools.ts b/libraries/diff-studio-utils/src/solver-tools/callbacks/callback-tools.ts new file mode 100644 index 0000000000..68644a0ce9 --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/callbacks/callback-tools.ts @@ -0,0 +1,20 @@ +// Manager of callbacks + +import {SolverOptions} from '../solver-defs'; +import {Callback} from './callback-base'; +import {TimeCheckerCallback} from './time-checker-callback'; +import {IterCheckerCallback} from './iter-checker-callback'; + +/** Get callback corresponding to the options */ +export function getCallback(options?: Partial): Callback | undefined { + if (options === undefined) + return undefined; + + if (options.maxIterations !== undefined) + return new IterCheckerCallback(options.maxIterations); + + if (options.maxTimeMs !== undefined) + return new TimeCheckerCallback(options.maxTimeMs); + + return undefined; +} diff --git a/libraries/diff-studio-utils/src/solver-tools/callbacks/iter-checker-callback.ts b/libraries/diff-studio-utils/src/solver-tools/callbacks/iter-checker-callback.ts new file mode 100644 index 0000000000..c6ef489ef7 --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/callbacks/iter-checker-callback.ts @@ -0,0 +1,21 @@ +import {Callback} from './callback-base'; +import {CallbackAction} from '../solver-defs'; + +/** This callback terminates computations if the maximum iterations is exceeded */ +export class IterCheckerCallback extends Callback { + private maxIter: number; + private currentIter: number; + + constructor(maxIter: number) { + super(); + this.maxIter = maxIter; + this.currentIter = 0; + } + + public onIterationStart(): void { + ++this.currentIter; + + if (this.currentIter > this.maxIter) + throw new CallbackAction(`Max iterations count exceeded (${this.maxIter})`); + } +} diff --git a/libraries/diff-studio-utils/src/solver-tools/callbacks/time-checker-callback.ts b/libraries/diff-studio-utils/src/solver-tools/callbacks/time-checker-callback.ts new file mode 100644 index 0000000000..27f91883b4 --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/callbacks/time-checker-callback.ts @@ -0,0 +1,19 @@ +import {Callback} from './callback-base'; +import {CallbackAction} from '../solver-defs'; + +/** This callback terminates computations if the maximum calculation time is exceeded */ +export class TimeCheckerCallback extends Callback { + private maxTime = 0; + private startingTime = 0; + + constructor(maxTime: number) { + super(); + this.maxTime = maxTime; + this.startingTime = performance.now(); + } + + public onIterationStart(): void { + if (performance.now() - this.startingTime > this.maxTime) + throw new CallbackAction(`Max time exceeded (${this.maxTime} ms)`); + } +} diff --git a/libraries/diff-studio-utils/src/solver-tools/lin-alg-tools.js b/libraries/diff-studio-utils/src/solver-tools/lin-alg-tools.js new file mode 100644 index 0000000000..4d3f5d1ead --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/lin-alg-tools.js @@ -0,0 +1,61 @@ +"use strict"; +// LU-decomposition tools +Object.defineProperty(exports, "__esModule", { value: true }); +exports.solve1d2d = exports.luSolve = exports.luDecomp = void 0; +/** Compute LU-decomposition of matrix */ +function luDecomp(A, L, U, n) { + L.fill(0); + U.fill(0); + for (var i_1 = 0; i_1 < n; ++i_1) + L[i_1 * n + i_1] = 1; + var sumU = 0; + var sumL = 0; + var k = 0; + var j = 0; + var p = 0; + var i = 0; + for (k = 0; k < n; ++k) { + for (j = k; j < n; ++j) { + sumU = 0; + for (p = 0; p < k; ++p) + sumU += L[k * n + p] * U[p * n + j]; + U[k * n + j] = A[k * n + j] - sumU; + } + for (i = k + 1; i < n; ++i) { + sumL = 0; + for (p = 0; p < k; ++p) + sumL += L[i * n + p] * U[p * n + k]; + L[i * n + k] = (A[i * n + k] - sumL) / U[k * n + k]; + } + } +} // luDecomp +exports.luDecomp = luDecomp; +/** Solve the system Ax = b using pre-computed LU-decomposition */ +function luSolve(L, U, b, y, x, n) { + var sumLy = 0; + for (var i = 0; i < n; ++i) { + sumLy = 0; + for (var j = 0; j < i; ++j) + sumLy += L[i * n + j] * y[j]; + y[i] = (b[i] - sumLy) / L[i * n + i]; + } + var sumUx = 0; + for (var i = n - 1; i > -1; --i) { + sumUx = 0; + for (var j = i + 1; j < n; ++j) + sumUx += U[i * n + j] * x[j]; + x[i] = (y[i] - sumUx) / U[i * n + i]; + } +} // luSolve +exports.luSolve = luSolve; +/** Solve Ax = b for 1D & 2D cases */ +function solve1d2d(A, b, x) { + if (x.length === 1) { + x[0] = b[0] / A[0]; + return; + } + var delta = A[0] * A[3] - A[1] * A[2]; + x[0] = (b[0] * A[3] - b[1] * A[1]) / delta; + x[1] = (A[0] * b[1] - A[2] * b[0]) / delta; +} +exports.solve1d2d = solve1d2d; diff --git a/libraries/diff-studio-utils/src/solver-tools/lin-alg-tools.ts b/libraries/diff-studio-utils/src/solver-tools/lin-alg-tools.ts new file mode 100644 index 0000000000..c2801513be --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/lin-alg-tools.ts @@ -0,0 +1,75 @@ +// LU-decomposition tools + +/** Compute LU-decomposition of matrix */ +export function luDecomp(A: Float64Array, L: Float64Array, U: Float64Array, n: number) { + L.fill(0); + U.fill(0); + + for (let i = 0; i < n; ++i) + L[i * n + i] = 1; + + let sumU = 0; + let sumL = 0; + let k = 0; + let j = 0; + let p = 0; + let i =0; + + for (k = 0; k < n; ++k) { + for (j = k; j < n; ++j) { + sumU = 0; + + for (p = 0; p < k; ++p) + sumU += L[k * n + p] * U[p * n + j]; + + U[k * n + j] = A[k * n +j] - sumU; + } + + for (i = k + 1; i < n; ++i) { + sumL = 0; + + for (p = 0; p < k; ++p) + sumL += L[i * n + p] * U[p * n + k]; + + L[i * n + k] = (A[i * n + k] - sumL) / U[k * n + k]; + } + } +} // luDecomp + +/** Solve the system Ax = b using pre-computed LU-decomposition */ +export function luSolve(L: Float64Array, U: Float64Array, b: Float64Array, y: Float64Array, + x: Float64Array, n: number) { + let sumLy = 0; + + for (let i = 0; i < n; ++i) { + sumLy = 0; + + for (let j = 0; j < i; ++j) + sumLy += L[i * n + j] * y[j]; + + y[i] = (b[i] - sumLy) / L[i * n + i]; + } + + let sumUx = 0; + + for (let i = n - 1; i > -1; --i) { + sumUx = 0; + + for (let j = i + 1; j < n; ++j) + sumUx += U[i * n + j] * x[j]; + + x[i] = (y[i] - sumUx) / U[i * n + i]; + } +} // luSolve + +/** Solve Ax = b for 1D & 2D cases */ +export function solve1d2d(A: Float64Array, b: Float64Array, x: Float64Array) { + if (x.length === 1) { + x[0] = b[0] / A[0]; + return; + } + + const delta = A[0] * A[3] - A[1] * A[2]; + x[0] = (b[0] * A[3] - b[1] * A[1]) / delta; + x[1] = (A[0] * b[1] - A[2] * b[0]) / delta; +} diff --git a/libraries/diff-studio-utils/src/solver-tools/mrt-method.js b/libraries/diff-studio-utils/src/solver-tools/mrt-method.js new file mode 100644 index 0000000000..8cb3b2a4cb --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/mrt-method.js @@ -0,0 +1,201 @@ +"use strict"; +/* The modified Rosenbrock triple method (MTR). + + References: + [1] https://doi.org/10.1137/S1064827594276424 + [2] https://doi.org/10.1016/S0898-1221(00)00175-9 +*/ +Object.defineProperty(exports, "__esModule", { value: true }); +exports.mrt = void 0; +var solver_defs_1 = require("./solver-defs"); +var lin_alg_tools_1 = require("./lin-alg-tools"); +// Quantities used in Rosenbrock method (see [1], [2] for more details) +var D = 1.0 - Math.sqrt(2.0) / 2.0; +var E32 = 6.0 + Math.sqrt(2.0); +/** Solve initial value problem the modified Rosenbrock triple (MRT) method [1, 2] */ +function mrt(odes, callback) { + /** right-hand side of the IVP solved */ + var f = odes.func; + // operating variables + var t0 = odes.arg.start; + var t1 = odes.arg.finish; + var h = odes.arg.step; + var hDataframe = h; + var tolerance = odes.tolerance; + /** number of solution dataframe rows */ + var rowCount = Math.trunc((t1 - t0) / h) + 1; + /** dimension of the problem */ + var dim = odes.initial.length; + var dimSquared = dim * dim; + /** independent variable values */ + var tArr = new Float64Array(rowCount); + /** arrays of solution values */ + var yArrs = Array(dim); + for (var i = 0; i < dim; ++i) + yArrs[i] = new Float64Array(rowCount); + // method routine + var timeDataframe = t0 + hDataframe; + var t = t0; + var tPrev = t0; + var hNext = 0.0; + var flag = true; + var index = 1; + var errmax = 0; + var hTemp = 0; + var tNew = 0; + // 0 BUFFERS & TEMP STRUCTURES + /** identity matrix */ + var I = new Float64Array(dim * dim); + // compute identity matrix + for (var i = 0; i < dim; ++i) { + for (var j = 0; j < dim; ++j) + I[j + i * dim] = (i === j) ? 1 : 0; + } + var y = new Float64Array(odes.initial); + var yPrev = new Float64Array(odes.initial); + var dydt = new Float64Array(dim); + var yScale = new Float64Array(dim); + var yTemp = new Float64Array(dim); + var yErr = new Float64Array(dim); + var W = new Float64Array(dimSquared); + var f0 = new Float64Array(dim); + var k1 = new Float64Array(dim); + var f1 = new Float64Array(dim); + var k2 = new Float64Array(dim); + var f2 = new Float64Array(dim); + var k3 = new Float64Array(dim); + var yDer = new Float64Array(dim); + var hdT = new Float64Array(dim); + var f0Buf = new Float64Array(dim); + var f1Buf = new Float64Array(dim); + var hd = 0; + var hDivNum = 0; + var L = new Float64Array(dimSquared); + var U = new Float64Array(dimSquared); + var luBuf = new Float64Array(dim); + var toUseLU = dim > 2; + // 1. SOLUTION AT THE POINT t0 + tArr[0] = t0; + for (var i = 0; i < dim; ++i) + yArrs[i][0] = y[i]; + // 2. COMPUTE NUMERICAL SOLUTION FOR THE POINTS FROM THE INTERVAL (t0, t1) + while (flag) { + // compute derivative + f(t, y, dydt); + // check whether to go on computations + if (callback) + callback.onIterationStart(); + // compute scale vector + for (var i = 0; i < dim; ++i) + yScale[i] = (0, solver_defs_1.abs)(y[i]) + h * (0, solver_defs_1.abs)(dydt[i]) + solver_defs_1.TINY; + // check end point + if (t + h > t1) { + h = t1 - t; + flag = false; + } + // call of adaptive step modified Rosenbrok triple method + // computation of solution (y), time (t) and next step (hNext) + while (true) { + // one stage of the modified Rosenbrok triple approach + // hdT = h * d * T(t, y, EPS); + (0, solver_defs_1.tDerivative)(t, y, f, solver_defs_1.EPS, f0Buf, f1Buf, hdT); + hd = h * D; + for (var i = 0; i < dim; ++i) + hdT[i] *= hd; + // The main computations + // f0 = f(t, y); + f(t, y, f0); + // W = I - h * d * J(t, y, EPS); + (0, solver_defs_1.jacobian)(t, y, f, solver_defs_1.EPS, f0Buf, f1Buf, W); + for (var i = 0; i < dimSquared; ++i) + W[i] = I[i] - hd * W[i]; + // compute LU-decomposition + if (toUseLU) + (0, lin_alg_tools_1.luDecomp)(W, L, U, dim); + // compute k1: solve the system W * k1 = f0 + hdT + for (var i = 0; i < dim; ++i) + f0Buf[i] = f0[i] + hdT[i]; + if (toUseLU) + (0, lin_alg_tools_1.luSolve)(L, U, f0Buf, luBuf, k1, dim); + else + (0, lin_alg_tools_1.solve1d2d)(W, f0Buf, k1); + hDivNum = 0.5 * h; + // yDer = y + 0.5 * h * k1; + for (var i = 0; i < dim; ++i) + yDer[i] = y[i] + hDivNum * k1[i]; + // f1 = f(t + 0.5 * h, yDer); + f(t + hDivNum, yDer, f1); + // compute k2: solve the system W * (k2 - k1) = f1 - k1 + for (var i = 0; i < dim; ++i) + f1Buf[i] = f1[i] - k1[i]; + if (toUseLU) + (0, lin_alg_tools_1.luSolve)(L, U, f1Buf, luBuf, k2, dim); + else + (0, lin_alg_tools_1.solve1d2d)(W, f1Buf, k2); + for (var i = 0; i < dim; ++i) + k2[i] = k2[i] + k1[i]; + // yOut = y + k2 * h; <--> yTemp + for (var i = 0; i < dim; ++i) + yTemp[i] = y[i] + h * k2[i]; + // f2 = f(t + h, yOut); + f(t + h, yTemp, f2); + // compute k3: solve the system W * k3 = f2 - e32 * (k2 - f1) - 2.0 * (k1 - f0) + hdT + for (var i = 0; i < dim; ++i) + f1Buf[i] = f2[i] - E32 * (k2[i] - f1[i]) - 2.0 * (k1[i] - f0[i]) + hdT[i]; + if (toUseLU) + (0, lin_alg_tools_1.luSolve)(L, U, f1Buf, luBuf, k3, dim); + else + (0, lin_alg_tools_1.solve1d2d)(W, f1Buf, k3); + // yErr = (k1 - 2.0 * k2 + k3) * h / 6; + hDivNum = h / 6; + for (var i = 0; i < dim; ++i) + yErr[i] = (k1[i] - 2.0 * k2[i] + k3[i]) * hDivNum; + // estimating error + errmax = 0; + for (var i = 0; i < dim; ++i) + errmax = (0, solver_defs_1.max)(errmax, (0, solver_defs_1.abs)(yErr[i] / yScale[i])); + errmax /= tolerance; + // processing the error obtained + if (errmax > 1) { + hTemp = solver_defs_1.SAFETY * h * Math.pow(errmax, solver_defs_1.PSHRNK); + h = (0, solver_defs_1.max)(hTemp, solver_defs_1.REDUCE_COEF * h); + tNew = t + h; + if (tNew == t) + throw new Error(solver_defs_1.ERROR_MSG.MRT_FAILS); + } + else { + if (errmax > solver_defs_1.ERR_CONTR) + hNext = solver_defs_1.SAFETY * h * Math.pow(errmax, solver_defs_1.PSGROW); + else + hNext = solver_defs_1.GROW_COEF * h; + t = t + h; + for (var i = 0; i < dim; ++i) + y[i] = yTemp[i]; + break; + } + } // while (true) + // compute lineraly interpolated results and store them in dataframe + while (timeDataframe < t) { + var cLeft = (t - timeDataframe) / (t - tPrev); + var cRight = 1.0 - cLeft; + tArr[index] = timeDataframe; + for (var j = 0; j < dim; ++j) + yArrs[j][index] = cRight * y[j] + cLeft * yPrev[j]; + timeDataframe += hDataframe; + ++index; + } + h = hNext; + tPrev = t; + for (var i = 0; i < dim; ++i) + yPrev[i] = y[i]; + } // while (flag) + // perform final callback actions + if (callback) + callback.onComputationsCompleted(); + // 3. solution at the point t1 + tArr[rowCount - 1] = t1; + for (var i = 0; i < dim; ++i) + yArrs[i][rowCount - 1] = y[i]; + return [tArr].concat(yArrs); +} // MTR +exports.mrt = mrt; diff --git a/libraries/diff-studio-utils/src/solver-tools/mrt-method.ts b/libraries/diff-studio-utils/src/solver-tools/mrt-method.ts new file mode 100644 index 0000000000..3b5096791f --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/mrt-method.ts @@ -0,0 +1,254 @@ +/* The modified Rosenbrock triple method (MTR). + + References: + [1] https://doi.org/10.1137/S1064827594276424 + [2] https://doi.org/10.1016/S0898-1221(00)00175-9 +*/ + +import {ODEs, max, abs, SAFETY, PSHRNK, PSGROW, REDUCE_COEF, GROW_COEF, + ERR_CONTR, TINY, EPS, tDerivative, jacobian, ERROR_MSG} from './solver-defs'; +import {Callback} from './callbacks/callback-base'; +import {luDecomp, luSolve, solve1d2d} from './lin-alg-tools'; + +// Quantities used in Rosenbrock method (see [1], [2] for more details) +const D = 1.0 - Math.sqrt(2.0) / 2.0; +const E32 = 6.0 + Math.sqrt(2.0); + +/** Solve initial value problem the modified Rosenbrock triple (MRT) method [1, 2] */ +export function mrt(odes: ODEs, callback?: Callback): Float64Array[] { + /** right-hand side of the IVP solved */ + const f = odes.func; + + // operating variables + const t0 = odes.arg.start; + const t1 = odes.arg.finish; + let h = odes.arg.step; + const hDataframe = h; + const tolerance = odes.tolerance; + + /** number of solution dataframe rows */ + const rowCount = Math.trunc((t1 - t0) / h) + 1; + + /** dimension of the problem */ + const dim = odes.initial.length; + const dimSquared = dim * dim; + + /** independent variable values */ + const tArr = new Float64Array(rowCount); + + /** arrays of solution values */ + const yArrs = Array(dim); + + for (let i = 0; i < dim; ++i) + yArrs[i] = new Float64Array(rowCount); + + // method routine + let timeDataframe = t0 + hDataframe; + let t = t0; + let tPrev = t0; + let hNext = 0.0; + let flag = true; + let index = 1; + let errmax = 0; + let hTemp = 0; + let tNew = 0; + + // 0 BUFFERS & TEMP STRUCTURES + + /** identity matrix */ + const I = new Float64Array(dim * dim); + + // compute identity matrix + for (let i = 0; i < dim; ++i) { + for (let j = 0; j < dim; ++j) + I[j + i * dim] = (i === j) ? 1 : 0; + } + + const y = new Float64Array(odes.initial); + const yPrev = new Float64Array(odes.initial); + const dydt = new Float64Array(dim); + const yScale = new Float64Array(dim); + const yTemp = new Float64Array(dim); + const yErr = new Float64Array(dim); + + const W = new Float64Array(dimSquared); + + const f0 = new Float64Array(dim); + const k1 = new Float64Array(dim); + const f1 = new Float64Array(dim); + const k2 = new Float64Array(dim); + const f2 = new Float64Array(dim); + const k3 = new Float64Array(dim); + const yDer = new Float64Array(dim); + const hdT = new Float64Array(dim); + + const f0Buf = new Float64Array(dim); + const f1Buf = new Float64Array(dim); + let hd = 0; + let hDivNum = 0; + + const L = new Float64Array(dimSquared); + const U = new Float64Array(dimSquared); + const luBuf = new Float64Array(dim); + const toUseLU = dim > 2; + + // 1. SOLUTION AT THE POINT t0 + tArr[0] = t0; + for (let i = 0; i < dim; ++i) + yArrs[i][0] = y[i]; + + // 2. COMPUTE NUMERICAL SOLUTION FOR THE POINTS FROM THE INTERVAL (t0, t1) + while (flag) { + // compute derivative + f(t, y, dydt); + + // check whether to go on computations + if (callback) + callback.onIterationStart(); + + // compute scale vector + for (let i = 0; i < dim; ++i) + yScale[i] = abs(y[i]) + h * abs(dydt[i]) + TINY; + + // check end point + if (t + h > t1) { + h = t1 - t; + flag = false; + } + + // call of adaptive step modified Rosenbrok triple method + // computation of solution (y), time (t) and next step (hNext) + while (true) { + // one stage of the modified Rosenbrok triple approach + // hdT = h * d * T(t, y, EPS); + tDerivative(t, y, f, EPS, f0Buf, f1Buf, hdT); + hd = h * D; + for (let i = 0; i < dim; ++i) + hdT[i] *= hd; + + // The main computations + + // f0 = f(t, y); + f(t, y, f0); + + // W = I - h * d * J(t, y, EPS); + jacobian(t, y, f, EPS, f0Buf, f1Buf, W); + for (let i = 0; i < dimSquared; ++i) + W[i] = I[i] - hd * W[i]; + + // compute LU-decomposition + if (toUseLU) + luDecomp(W, L, U, dim); + + // compute k1: solve the system W * k1 = f0 + hdT + for (let i = 0; i < dim; ++i) + f0Buf[i] = f0[i] + hdT[i]; + + if (toUseLU) + luSolve(L, U, f0Buf, luBuf, k1, dim); + else + solve1d2d(W, f0Buf, k1); + + hDivNum = 0.5 * h; + + // yDer = y + 0.5 * h * k1; + for (let i = 0; i < dim; ++i) + yDer[i] = y[i] + hDivNum * k1[i]; + + // f1 = f(t + 0.5 * h, yDer); + f(t + hDivNum, yDer, f1); + + // compute k2: solve the system W * (k2 - k1) = f1 - k1 + for (let i = 0; i < dim; ++i) + f1Buf[i] = f1[i] - k1[i]; + + if (toUseLU) + luSolve(L, U, f1Buf, luBuf, k2, dim); + else + solve1d2d(W, f1Buf, k2); + + for (let i = 0; i < dim; ++i) + k2[i] = k2[i] + k1[i]; + + // yOut = y + k2 * h; <--> yTemp + for (let i = 0; i < dim; ++i) + yTemp[i] = y[i] + h * k2[i]; + + // f2 = f(t + h, yOut); + f(t + h, yTemp, f2); + + // compute k3: solve the system W * k3 = f2 - e32 * (k2 - f1) - 2.0 * (k1 - f0) + hdT + for (let i = 0; i < dim; ++i) + f1Buf[i] = f2[i] - E32 * (k2[i] - f1[i]) - 2.0 * (k1[i] - f0[i]) + hdT[i]; + + if (toUseLU) + luSolve(L, U, f1Buf, luBuf, k3, dim); + else + solve1d2d(W, f1Buf, k3); + + // yErr = (k1 - 2.0 * k2 + k3) * h / 6; + hDivNum = h / 6; + + for (let i = 0; i < dim; ++i) + yErr[i] = (k1[i] - 2.0 * k2[i] + k3[i]) * hDivNum; + + // estimating error + errmax = 0; + for (let i = 0; i < dim; ++i) + errmax = max(errmax, abs(yErr[i] / yScale[i])); + errmax /= tolerance; + + // processing the error obtained + if (errmax > 1) { + hTemp = SAFETY * h * errmax**PSHRNK; + h = max(hTemp, REDUCE_COEF * h); + tNew = t + h; + if (tNew == t) + throw new Error(ERROR_MSG.MRT_FAILS); + } else { + if (errmax > ERR_CONTR) + hNext = SAFETY * h * errmax**PSGROW; + else + hNext = GROW_COEF * h; + t = t + h; + + for (let i = 0; i < dim; ++i) + y[i] = yTemp[i]; + + break; + } + } // while (true) + + // compute lineraly interpolated results and store them in dataframe + while (timeDataframe < t) { + const cLeft = (t - timeDataframe) / (t - tPrev); + const cRight = 1.0 - cLeft; + + tArr[index] = timeDataframe; + + for (let j = 0; j < dim; ++j) + yArrs[j][index] = cRight * y[j] + cLeft * yPrev[j]; + + timeDataframe += hDataframe; + ++index; + } + + h = hNext; + tPrev = t; + + for (let i = 0; i < dim; ++i) + yPrev[i] = y[i]; + } // while (flag) + + // perform final callback actions + if (callback) + callback.onComputationsCompleted(); + + // 3. solution at the point t1 + tArr[rowCount - 1] = t1; + + for (let i = 0; i < dim; ++i) + yArrs[i][rowCount - 1] = y[i]; + + return [tArr].concat(yArrs); +} // MTR diff --git a/libraries/diff-studio-utils/src/solver-tools/ros34prw-method.js b/libraries/diff-studio-utils/src/solver-tools/ros34prw-method.js new file mode 100644 index 0000000000..a43a44095d --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/ros34prw-method.js @@ -0,0 +1,256 @@ +"use strict"; +/* The ROS34PRw method implementation + + References: + [1] https://doi.org/10.1016/j.cam.2015.03.010 */ +Object.defineProperty(exports, "__esModule", { value: true }); +exports.ros34prw = void 0; +var solver_defs_1 = require("./solver-defs"); +var lin_alg_tools_1 = require("./lin-alg-tools"); +// The method specific constants (see Table 3 [1]): +var GAMMA = 0.435866521508459; +var GAMMA_21 = -1.3075995645253771; +var GAMMA_21_SCALED = GAMMA_21 / GAMMA; +var GAMMA_2 = GAMMA_21 + GAMMA; +var GAMMA_31 = -0.7098857586097217; +var GAMMA_31_SCALED = GAMMA_31 / GAMMA; +var GAMMA_32 = -0.55996735960277766; +var GAMMA_32_SCALED = GAMMA_32 / GAMMA; +var GAMMA_3 = GAMMA_31 + GAMMA_32 + GAMMA; +var GAMMA_41 = -0.15550856807552085; +var GAMMA_41_SCALED = GAMMA_41 / GAMMA; +var GAMMA_42 = -0.95388516575112225; +var GAMMA_42_SCALED = GAMMA_42 / GAMMA; +var GAMMA_43 = 0.67352721231818413; +var GAMMA_43_SCALED = GAMMA_43 / GAMMA; +var GAMMA_4 = GAMMA_41 + GAMMA_42 + GAMMA_43 + GAMMA; +var ALPHA_21 = 1.3075995645253771; +var ALPHA_2 = ALPHA_21; +/* REMARK. Expressions with ALPHAs are simplified */ +/*const ALPHA_31 = 0.5; +const ALPHA_32 = 0.5; +const ALPHA_3 = ALPHA_31 + ALPHA_32; +const ALPHA_41 = 0.5; +const ALPHA_42 = 0.5; +const ALPHA_43 = 0; +const ALPHA_4 = ALPHA_41 + ALPHA_42;*/ +var B_1 = 0.34449143192447917; +var B_2 = -0.45388516575112231; +var B_3 = 0.67352721231818413; +var B_4 = 0.435866521508459; +var B_HAT_1 = 0.5; +var B_HAT_2 = -0.25738812086522078; +var B_HAT_3 = 0.43542008724775044; +var B_HAT_4 = 0.32196803361747034; +var R_1 = B_1 - B_HAT_1; +var R_2 = B_2 - B_HAT_2; +var R_3 = B_3 - B_HAT_3; +var R_4 = B_4 - B_HAT_4; +/** Solve initial value problem using the ROS34PRw method [5]. */ +function ros34prw(odes, callback) { + /** right-hand side of the IVP solved */ + var f = odes.func; + // operating variables + var t0 = odes.arg.start; + var t1 = odes.arg.finish; + var h = odes.arg.step; + var hDataframe = h; + var tolerance = odes.tolerance; + /** number of solution dataframe rows */ + var rowCount = Math.trunc((t1 - t0) / h) + 1; + /** dimension of the problem */ + var dim = odes.initial.length; + var dimSquared = dim * dim; + /** independent variable values */ + var tArr = new Float64Array(rowCount); + /** arrays of solution values */ + var yArrs = Array(dim); + for (var i = 0; i < dim; ++i) + yArrs[i] = new Float64Array(rowCount); + // method routine + var timeDataframe = t0 + hDataframe; + var t = t0; + var tPrev = t0; + var flag = true; + var index = 1; + var errmax = 0; + var hTemp = 0; + var tNew = 0; + var hNext = 0.0; + // 0 BUFFERS & TEMP STRUCTURES + /** identity matrix */ + var I = new Float64Array(dim * dim); + // compute identity matrix + for (var i = 0; i < dim; ++i) { + for (var j = 0; j < dim; ++j) + I[j + i * dim] = (i === j) ? 1 : 0; + } + var y = new Float64Array(odes.initial); + var yPrev = new Float64Array(odes.initial); + var dydt = new Float64Array(dim); + var yScale = new Float64Array(dim); + var yTemp = new Float64Array(dim); + var yErr = new Float64Array(dim); + var W = new Float64Array(dimSquared); + var k1 = new Float64Array(dim); + var k2 = new Float64Array(dim); + var k3 = new Float64Array(dim); + var k4 = new Float64Array(dim); + var HT = new Float64Array(dim); + var f0Buf = new Float64Array(dim); + var f1Buf = new Float64Array(dim); + var hByGamma = 0; + var fBuf = new Float64Array(dim); + var kBuf = new Float64Array(dim); + var L = new Float64Array(dimSquared); + var U = new Float64Array(dimSquared); + var luBuf = new Float64Array(dim); + var bBuf = new Float64Array(dim); + var toUseLU = dim > 2; + // 1. SOLUTION AT THE POINT t0 + tArr[0] = t0; + for (var i = 0; i < dim; ++i) + yArrs[i][0] = y[i]; + // 2. COMPUTE NUMERICAL SOLUTION FOR THE POINTS FROM THE INTERVAL (t0, t1) + while (flag) { + // compute derivative + f(t, y, dydt); + // check whether to go on computations + if (callback) + callback.onIterationStart(); + // compute scale vector + for (var i = 0; i < dim; ++i) + yScale[i] = (0, solver_defs_1.abs)(y[i]) + h * (0, solver_defs_1.abs)(dydt[i]) + solver_defs_1.TINY; + // check end point + if (t + h > t1) { + h = t1 - t; + flag = false; + } + // call of adaptive step the ROS3Pw [1] method + // computation of solution (y), time (t) and next step (hNext) + while (true) { + hByGamma = h * GAMMA; + // one stage of the ROS3Pw method + // 1) Jacobian & dF/dt matrices + (0, solver_defs_1.jacobian)(t, y, f, solver_defs_1.EPS, f0Buf, f1Buf, W); + (0, solver_defs_1.tDerivative)(t, y, f, solver_defs_1.EPS, f0Buf, f1Buf, HT); + // 2) W & LU-decomposition + for (var i = 0; i < dimSquared; ++i) + W[i] = I[i] - hByGamma * W[i]; + if (toUseLU) + (0, lin_alg_tools_1.luDecomp)(W, L, U, dim); + // 3) Scale dF/dt: HT = j * T + for (var i = 0; i < dim; ++i) + HT[i] *= h; + // 4) F1 = F(t, y) <-- Fbuf + f(t, y, fBuf); + // 5) k1 = W_inv * (F1 + gamma * HT) + for (var i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + GAMMA * HT[i]; + if (toUseLU) + (0, lin_alg_tools_1.luSolve)(L, U, bBuf, luBuf, k1, dim); + else + (0, lin_alg_tools_1.solve1d2d)(W, bBuf, k1); + // 6) F2 = F(t + alpha2 * h, y + alpha21 * k1) <-- Fbuf + for (var i = 0; i < dim; ++i) + kBuf[i] = y[i] + ALPHA_21 * h * k1[i]; + f(t + ALPHA_2 * h, kBuf, fBuf); + // 7) kBuf = gamma21 / gamma * k1 + for (var i = 0; i < dim; ++i) + kBuf[i] = GAMMA_21_SCALED * k1[i]; + // 8) k2 = W_inv * [Fbuf + kBuf + gamma2 * HT] - kBuf + for (var i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + kBuf[i] + GAMMA_2 * HT[i]; + if (toUseLU) + (0, lin_alg_tools_1.luSolve)(L, U, bBuf, luBuf, k2, dim); + else + (0, lin_alg_tools_1.solve1d2d)(W, bBuf, k2); + for (var i = 0; i < dim; ++i) + k2[i] -= kBuf[i]; + // 9) F3 = F(t + alpha3 * h, y + h * (alpha31 * k1 + alpha32 * k2)) <-- Fbuf + for (var i = 0; i < dim; ++i) + kBuf[i] = y[i] + h * 0.5 * (k1[i] + k2[i]); + f(t + h, kBuf, fBuf); + // 10) kBuf = gamma31 / gamma * k1 + gamma32 / gamma * k2 + for (var i = 0; i < dim; ++i) + kBuf[i] = GAMMA_31_SCALED * k1[i] + GAMMA_32_SCALED * k2[i]; + // 11) k3 = W_inv * (F3 + kBuf + gamma3 * HT) - kBuf + for (var i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + kBuf[i] + GAMMA_3 * HT[i]; + if (toUseLU) + (0, lin_alg_tools_1.luSolve)(L, U, bBuf, luBuf, k3, dim); + else + (0, lin_alg_tools_1.solve1d2d)(W, bBuf, k3); + for (var i = 0; i < dim; ++i) + k3[i] -= kBuf[i]; + // 12) F4 = F(t + alpha4 * h, y + h * (alpha41 * k1 + alpha42 * k2 + alpha43 * k3)) <-- Fbuf + for (var i = 0; i < dim; ++i) + kBuf[i] = y[i] + h * 0.5 * (k1[i] + k2[i]); + f(t + h, kBuf, fBuf); + // 13) kBuf = gamma41 / gamma * k1 + gamma42 / gamma * k2 + gamma43 / gamma * k3 + for (var i = 0; i < dim; ++i) + kBuf[i] = GAMMA_41_SCALED * k1[i] + GAMMA_42_SCALED * k2[i] + GAMMA_43_SCALED * k3[i]; + // 14) k4 = W_inv * (F4 + kBuf + gamma4 * HT) - kBuf + for (var i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + kBuf[i] + GAMMA_4 * HT[i]; + if (toUseLU) + (0, lin_alg_tools_1.luSolve)(L, U, bBuf, luBuf, k4, dim); + else + (0, lin_alg_tools_1.solve1d2d)(W, bBuf, k4); + for (var i = 0; i < dim; ++i) + k4[i] -= kBuf[i]; + // 15) yNext = y + h * (b1 * k1 + b2 * k2 + b3 * k3 + b4 * k4) <-- yTemp + for (var i = 0; i < dim; ++i) + yTemp[i] = y[i] + h * (B_1 * k1[i] + B_2 * k2[i] + B_3 * k3[i] + B_4 * k4[i]); + // 16) yErr = h * (r1 * k1 + r2 * k2 + r3 * k3 + r4 * k4) + for (var i = 0; i < dim; ++i) + yErr[i] = h * (R_1 * k1[i] + R_2 * k2[i] + R_3 * k3[i] + R_4 * k4[i]); + // estimating error + errmax = 0; + for (var i = 0; i < dim; ++i) + errmax = (0, solver_defs_1.max)(errmax, (0, solver_defs_1.abs)(yErr[i] / yScale[i])); + errmax /= tolerance; + // processing the error obtained + if (errmax > 1) { + hTemp = solver_defs_1.SAFETY * h * Math.pow(errmax, solver_defs_1.PSHRNK); + h = (0, solver_defs_1.max)(hTemp, solver_defs_1.REDUCE_COEF * h); + tNew = t + h; + if (tNew == t) + throw new Error(solver_defs_1.ERROR_MSG.ROS34PRW_FAILS); + } + else { + if (errmax > solver_defs_1.ERR_CONTR) + hNext = solver_defs_1.SAFETY * h * Math.pow(errmax, solver_defs_1.PSGROW); + else + hNext = solver_defs_1.GROW_COEF * h; + t = t + h; + for (var i = 0; i < dim; ++i) + y[i] = yTemp[i]; + break; + } + } // while (true) + // compute lineraly interpolated results and store them in dataframe + while (timeDataframe < t) { + var cLeft = (t - timeDataframe) / (t - tPrev); + var cRight = 1.0 - cLeft; + tArr[index] = timeDataframe; + for (var j = 0; j < dim; ++j) + yArrs[j][index] = cRight * y[j] + cLeft * yPrev[j]; + timeDataframe += hDataframe; + ++index; + } + h = hNext; + tPrev = t; + for (var i = 0; i < dim; ++i) + yPrev[i] = y[i]; + } // while (flag) + // perform final callback actions + if (callback) + callback.onComputationsCompleted(); + // 3. solution at the point t1 + tArr[rowCount - 1] = t1; + for (var i = 0; i < dim; ++i) + yArrs[i][rowCount - 1] = y[i]; + return [tArr].concat(yArrs); +} // ros34prw +exports.ros34prw = ros34prw; diff --git a/libraries/diff-studio-utils/src/solver-tools/ros34prw-method.ts b/libraries/diff-studio-utils/src/solver-tools/ros34prw-method.ts new file mode 100644 index 0000000000..09c6b5b745 --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/ros34prw-method.ts @@ -0,0 +1,329 @@ +/* The ROS34PRw method implementation + + References: + [1] https://doi.org/10.1016/j.cam.2015.03.010 */ + +import {ODEs, max, abs, SAFETY, PSHRNK, PSGROW, REDUCE_COEF, GROW_COEF, + ERR_CONTR, TINY, EPS, tDerivative, jacobian, ERROR_MSG} from './solver-defs'; +import {Callback} from './callbacks/callback-base'; +import {luDecomp, luSolve, solve1d2d} from './lin-alg-tools'; + +// The method specific constants (see Table 3 [1]): +const GAMMA = 0.435866521508459; + +const GAMMA_21 = -1.3075995645253771; +const GAMMA_21_SCALED = GAMMA_21 / GAMMA; +const GAMMA_2 = GAMMA_21 + GAMMA; + +const GAMMA_31 = -0.7098857586097217; +const GAMMA_31_SCALED = GAMMA_31 / GAMMA; +const GAMMA_32 = -0.55996735960277766; +const GAMMA_32_SCALED = GAMMA_32 / GAMMA; +const GAMMA_3 = GAMMA_31 + GAMMA_32 + GAMMA; + +const GAMMA_41 = -0.15550856807552085; +const GAMMA_41_SCALED = GAMMA_41 / GAMMA; +const GAMMA_42 = -0.95388516575112225; +const GAMMA_42_SCALED = GAMMA_42 / GAMMA; +const GAMMA_43 = 0.67352721231818413; +const GAMMA_43_SCALED = GAMMA_43 / GAMMA; +const GAMMA_4 = GAMMA_41 + GAMMA_42 + GAMMA_43 + GAMMA; + +const ALPHA_21 = 1.3075995645253771; +const ALPHA_2 = ALPHA_21; + +/* REMARK. Expressions with ALPHAs are simplified */ +/*const ALPHA_31 = 0.5; +const ALPHA_32 = 0.5; +const ALPHA_3 = ALPHA_31 + ALPHA_32; +const ALPHA_41 = 0.5; +const ALPHA_42 = 0.5; +const ALPHA_43 = 0; +const ALPHA_4 = ALPHA_41 + ALPHA_42;*/ + +const B_1 = 0.34449143192447917; +const B_2 = -0.45388516575112231; +const B_3 = 0.67352721231818413; +const B_4 = 0.435866521508459; + +const B_HAT_1 = 0.5; +const B_HAT_2 = -0.25738812086522078; +const B_HAT_3 = 0.43542008724775044; +const B_HAT_4 = 0.32196803361747034; + +const R_1 = B_1 - B_HAT_1; +const R_2 = B_2 - B_HAT_2; +const R_3 = B_3 - B_HAT_3; +const R_4 = B_4 - B_HAT_4; + + +/** Solve initial value problem using the ROS34PRw method [5]. */ +export function ros34prw(odes: ODEs, callback?: Callback): Float64Array[] { + /** right-hand side of the IVP solved */ + const f = odes.func; + + // operating variables + const t0 = odes.arg.start; + const t1 = odes.arg.finish; + let h = odes.arg.step; + const hDataframe = h; + const tolerance = odes.tolerance; + + /** number of solution dataframe rows */ + const rowCount = Math.trunc((t1 - t0) / h) + 1; + + /** dimension of the problem */ + const dim = odes.initial.length; + const dimSquared = dim * dim; + + /** independent variable values */ + const tArr = new Float64Array(rowCount); + + /** arrays of solution values */ + const yArrs = Array(dim); + + for (let i = 0; i < dim; ++i) + yArrs[i] = new Float64Array(rowCount); + + // method routine + let timeDataframe = t0 + hDataframe; + let t = t0; + let tPrev = t0; + let flag = true; + let index = 1; + let errmax = 0; + let hTemp = 0; + let tNew = 0; + let hNext = 0.0; + + // 0 BUFFERS & TEMP STRUCTURES + + /** identity matrix */ + const I = new Float64Array(dim * dim); + + // compute identity matrix + for (let i = 0; i < dim; ++i) { + for (let j = 0; j < dim; ++j) + I[j + i * dim] = (i === j) ? 1 : 0; + } + + const y = new Float64Array(odes.initial); + const yPrev = new Float64Array(odes.initial); + const dydt = new Float64Array(dim); + const yScale = new Float64Array(dim); + const yTemp = new Float64Array(dim); + const yErr = new Float64Array(dim); + + const W = new Float64Array(dimSquared); + + const k1 = new Float64Array(dim); + const k2 = new Float64Array(dim); + const k3 = new Float64Array(dim); + const k4 = new Float64Array(dim); + const HT = new Float64Array(dim); + + const f0Buf = new Float64Array(dim); + const f1Buf = new Float64Array(dim); + + let hByGamma = 0; + const fBuf = new Float64Array(dim); + const kBuf = new Float64Array(dim); + + const L = new Float64Array(dimSquared); + const U = new Float64Array(dimSquared); + const luBuf = new Float64Array(dim); + const bBuf = new Float64Array(dim); + const toUseLU = dim > 2; + + // 1. SOLUTION AT THE POINT t0 + tArr[0] = t0; + for (let i = 0; i < dim; ++i) + yArrs[i][0] = y[i]; + + // 2. COMPUTE NUMERICAL SOLUTION FOR THE POINTS FROM THE INTERVAL (t0, t1) + while (flag) { + // compute derivative + f(t, y, dydt); + + // check whether to go on computations + if (callback) + callback.onIterationStart(); + + // compute scale vector + for (let i = 0; i < dim; ++i) + yScale[i] = abs(y[i]) + h * abs(dydt[i]) + TINY; + + // check end point + if (t + h > t1) { + h = t1 - t; + flag = false; + } + + // call of adaptive step the ROS3Pw [1] method + // computation of solution (y), time (t) and next step (hNext) + while (true) { + hByGamma = h * GAMMA; + + // one stage of the ROS3Pw method + + // 1) Jacobian & dF/dt matrices + jacobian(t, y, f, EPS, f0Buf, f1Buf, W); + tDerivative(t, y, f, EPS, f0Buf, f1Buf, HT); + + // 2) W & LU-decomposition + for (let i = 0; i < dimSquared; ++i) + W[i] = I[i] - hByGamma * W[i]; + + if (toUseLU) + luDecomp(W, L, U, dim); + + // 3) Scale dF/dt: HT = j * T + for (let i = 0; i < dim; ++i) + HT[i] *= h; + + // 4) F1 = F(t, y) <-- Fbuf + f(t, y, fBuf); + + // 5) k1 = W_inv * (F1 + gamma * HT) + for (let i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + GAMMA * HT[i]; + + if (toUseLU) + luSolve(L, U, bBuf, luBuf, k1, dim); + else + solve1d2d(W, bBuf, k1); + + // 6) F2 = F(t + alpha2 * h, y + alpha21 * k1) <-- Fbuf + for (let i = 0; i < dim; ++i) + kBuf[i] = y[i] + ALPHA_21 * h * k1[i]; + + f(t + ALPHA_2 * h, kBuf, fBuf); + + // 7) kBuf = gamma21 / gamma * k1 + for (let i = 0; i < dim; ++i) + kBuf[i] = GAMMA_21_SCALED * k1[i]; + + // 8) k2 = W_inv * [Fbuf + kBuf + gamma2 * HT] - kBuf + for (let i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + kBuf[i] + GAMMA_2 * HT[i]; + + if (toUseLU) + luSolve(L, U, bBuf, luBuf, k2, dim); + else + solve1d2d(W, bBuf, k2); + + for (let i = 0; i < dim; ++i) + k2[i] -= kBuf[i]; + + // 9) F3 = F(t + alpha3 * h, y + h * (alpha31 * k1 + alpha32 * k2)) <-- Fbuf + for (let i = 0; i < dim; ++i) + kBuf[i] = y[i] + h * 0.5 * (k1[i] + k2[i]); + + f(t + h, kBuf, fBuf); + + // 10) kBuf = gamma31 / gamma * k1 + gamma32 / gamma * k2 + for (let i = 0; i < dim; ++i) + kBuf[i] = GAMMA_31_SCALED * k1[i] + GAMMA_32_SCALED * k2[i]; + + // 11) k3 = W_inv * (F3 + kBuf + gamma3 * HT) - kBuf + for (let i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + kBuf[i] + GAMMA_3 * HT[i]; + + if (toUseLU) + luSolve(L, U, bBuf, luBuf, k3, dim); + else + solve1d2d(W, bBuf, k3); + + for (let i = 0; i < dim; ++i) + k3[i] -= kBuf[i]; + + // 12) F4 = F(t + alpha4 * h, y + h * (alpha41 * k1 + alpha42 * k2 + alpha43 * k3)) <-- Fbuf + for (let i = 0; i < dim; ++i) + kBuf[i] = y[i] + h * 0.5 * (k1[i] + k2[i]); + + f(t + h, kBuf, fBuf); + + // 13) kBuf = gamma41 / gamma * k1 + gamma42 / gamma * k2 + gamma43 / gamma * k3 + for (let i = 0; i < dim; ++i) + kBuf[i] = GAMMA_41_SCALED * k1[i] + GAMMA_42_SCALED * k2[i] + GAMMA_43_SCALED * k3[i]; + + // 14) k4 = W_inv * (F4 + kBuf + gamma4 * HT) - kBuf + for (let i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + kBuf[i] + GAMMA_4 * HT[i]; + + if (toUseLU) + luSolve(L, U, bBuf, luBuf, k4, dim); + else + solve1d2d(W, bBuf, k4); + + for (let i = 0; i < dim; ++i) + k4[i] -= kBuf[i]; + + // 15) yNext = y + h * (b1 * k1 + b2 * k2 + b3 * k3 + b4 * k4) <-- yTemp + for (let i = 0; i < dim; ++i) + yTemp[i] = y[i] + h * (B_1 * k1[i] + B_2 * k2[i] + B_3 * k3[i] + B_4 * k4[i]); + + // 16) yErr = h * (r1 * k1 + r2 * k2 + r3 * k3 + r4 * k4) + for (let i = 0; i < dim; ++i) + yErr[i] = h * (R_1 * k1[i] + R_2 * k2[i] + R_3 * k3[i] + R_4 * k4[i]); + + // estimating error + errmax = 0; + for (let i = 0; i < dim; ++i) + errmax = max(errmax, abs(yErr[i] / yScale[i])); + errmax /= tolerance; + + // processing the error obtained + if (errmax > 1) { + hTemp = SAFETY * h * errmax**PSHRNK; + h = max(hTemp, REDUCE_COEF * h); + tNew = t + h; + if (tNew == t) + throw new Error(ERROR_MSG.ROS34PRW_FAILS); + } else { + if (errmax > ERR_CONTR) + hNext = SAFETY * h * errmax**PSGROW; + else + hNext = GROW_COEF * h; + t = t + h; + + for (let i = 0; i < dim; ++i) + y[i] = yTemp[i]; + + break; + } + } // while (true) + + // compute lineraly interpolated results and store them in dataframe + while (timeDataframe < t) { + const cLeft = (t - timeDataframe) / (t - tPrev); + const cRight = 1.0 - cLeft; + + tArr[index] = timeDataframe; + + for (let j = 0; j < dim; ++j) + yArrs[j][index] = cRight * y[j] + cLeft * yPrev[j]; + + timeDataframe += hDataframe; + ++index; + } + + h = hNext; + tPrev = t; + + for (let i = 0; i < dim; ++i) + yPrev[i] = y[i]; + } // while (flag) + + // perform final callback actions + if (callback) + callback.onComputationsCompleted(); + + // 3. solution at the point t1 + tArr[rowCount - 1] = t1; + + for (let i = 0; i < dim; ++i) + yArrs[i][rowCount - 1] = y[i]; + + return [tArr].concat(yArrs); +} // ros34prw diff --git a/libraries/diff-studio-utils/src/solver-tools/ros3prw-method.js b/libraries/diff-studio-utils/src/solver-tools/ros3prw-method.js new file mode 100644 index 0000000000..d47c6cae8e --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/ros3prw-method.js @@ -0,0 +1,224 @@ +"use strict"; +/* The ROS3PRw method implementation + + References: + [1] https://doi.org/10.1016/j.cam.2015.03.010 */ +Object.defineProperty(exports, "__esModule", { value: true }); +exports.ros3prw = void 0; +var solver_defs_1 = require("./solver-defs"); +var lin_alg_tools_1 = require("./lin-alg-tools"); +// The method specific constants (see Table 2 [1]) +var GAMMA = 0.78867513459481287; +var GAMMA_21 = -2.3660254037844388; +var GAMMA_21_SCALED = GAMMA_21 / GAMMA; +var GAMMA_2 = GAMMA_21 + GAMMA; +var GAMMA_31 = -0.86791218280355165; +var GAMMA_31_SCALED = GAMMA_31 / GAMMA; +var GAMMA_32 = -0.87306695894642317; +var GAMMA_32_SCALED = GAMMA_32 / GAMMA; +var GAMMA_3 = GAMMA_31 + GAMMA_32 + GAMMA; +var ALPHA_21 = 2.3660254037844388; +var ALPHA_2 = ALPHA_21; +var ALPHA_31 = 0.5; +var ALPHA_32 = 0.76794919243112270; +var ALPHA_3 = ALPHA_31 + ALPHA_32; +var B_1 = 0.50544867840851759; +var B_2 = -0.11571687603637559; +var B_3 = 0.610268197627858; +var B_HAT_1 = 0.28973180237214197; +var B_HAT_2 = 0.10000000000000001; +var B_HAT_3 = 0.610268197627858; +var R_1 = B_1 - B_HAT_1; +var R_2 = B_2 - B_HAT_2; +var R_3 = B_3 - B_HAT_3; +/** Solve initial value problem using the ROS3Pw method [5]. */ +function ros3prw(odes, callback) { + /** right-hand side of the IVP solved */ + var f = odes.func; + // operating variables + var t0 = odes.arg.start; + var t1 = odes.arg.finish; + var h = odes.arg.step; + var hDataframe = h; + var tolerance = odes.tolerance; + /** number of solution dataframe rows */ + var rowCount = Math.trunc((t1 - t0) / h) + 1; + /** dimension of the problem */ + var dim = odes.initial.length; + var dimSquared = dim * dim; + /** independent variable values */ + var tArr = new Float64Array(rowCount); + /** arrays of solution values */ + var yArrs = Array(dim); + for (var i = 0; i < dim; ++i) + yArrs[i] = new Float64Array(rowCount); + // method routine + var timeDataframe = t0 + hDataframe; + var t = t0; + var tPrev = t0; + var flag = true; + var index = 1; + var errmax = 0; + var hTemp = 0; + var tNew = 0; + var hNext = 0.0; + // 0 BUFFERS & TEMP STRUCTURES + /** identity matrix */ + var I = new Float64Array(dim * dim); + // compute identity matrix + for (var i = 0; i < dim; ++i) { + for (var j = 0; j < dim; ++j) + I[j + i * dim] = (i === j) ? 1 : 0; + } + var y = new Float64Array(odes.initial); + var yPrev = new Float64Array(odes.initial); + var dydt = new Float64Array(dim); + var yScale = new Float64Array(dim); + var yTemp = new Float64Array(dim); + var yErr = new Float64Array(dim); + var W = new Float64Array(dimSquared); + var k1 = new Float64Array(dim); + var k2 = new Float64Array(dim); + var k3 = new Float64Array(dim); + var HT = new Float64Array(dim); + var f0Buf = new Float64Array(dim); + var f1Buf = new Float64Array(dim); + var hByGamma = 0; + var fBuf = new Float64Array(dim); + var kBuf = new Float64Array(dim); + var L = new Float64Array(dimSquared); + var U = new Float64Array(dimSquared); + var luBuf = new Float64Array(dim); + var bBuf = new Float64Array(dim); + var toUseLU = dim > 2; + // 1. SOLUTION AT THE POINT t0 + tArr[0] = t0; + for (var i = 0; i < dim; ++i) + yArrs[i][0] = y[i]; + // 2. COMPUTE NUMERICAL SOLUTION FOR THE POINTS FROM THE INTERVAL (t0, t1) + while (flag) { + // compute derivative + f(t, y, dydt); + // check whether to go on computations + if (callback) + callback.onIterationStart(); + // compute scale vector + for (var i = 0; i < dim; ++i) + yScale[i] = (0, solver_defs_1.abs)(y[i]) + h * (0, solver_defs_1.abs)(dydt[i]) + solver_defs_1.TINY; + // check end point + if (t + h > t1) { + h = t1 - t; + flag = false; + } + // call of adaptive step the ROS3Pw [1] method + // computation of solution (y), time (t) and next step (hNext) + while (true) { + hByGamma = h * GAMMA; + // one stage of the ROS3Pw method + // 1) Jacobian & dF/dt matrices + (0, solver_defs_1.jacobian)(t, y, f, solver_defs_1.EPS, f0Buf, f1Buf, W); + (0, solver_defs_1.tDerivative)(t, y, f, solver_defs_1.EPS, f0Buf, f1Buf, HT); + // 2) W & LU-decomposition + for (var i = 0; i < dimSquared; ++i) + W[i] = I[i] - hByGamma * W[i]; + if (toUseLU) + (0, lin_alg_tools_1.luDecomp)(W, L, U, dim); + // 3) Scale dF/dt: HT = j * T + for (var i = 0; i < dim; ++i) + HT[i] *= h; + // 4) F1 = F(t, y) <-- Fbuf + f(t, y, fBuf); + // 5) k1 = W_inv * (F1 + gamma * HT) + for (var i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + GAMMA * HT[i]; + if (toUseLU) + (0, lin_alg_tools_1.luSolve)(L, U, bBuf, luBuf, k1, dim); + else + (0, lin_alg_tools_1.solve1d2d)(W, bBuf, k1); + // 6) F2 = F(t + alpha2 * h, y + alpha21 * k1) <-- Fbuf + for (var i = 0; i < dim; ++i) + kBuf[i] = y[i] + ALPHA_21 * h * k1[i]; + f(t + ALPHA_2 * h, kBuf, fBuf); + // 7) kBuf = gamma21 / gamma * k1 + for (var i = 0; i < dim; ++i) + kBuf[i] = GAMMA_21_SCALED * k1[i]; + // 8) k2 = W_inv * [Fbuf + kBuf + gamma2 * HT] - kBuf + for (var i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + kBuf[i] + GAMMA_2 * HT[i]; + if (toUseLU) + (0, lin_alg_tools_1.luSolve)(L, U, bBuf, luBuf, k2, dim); + else + (0, lin_alg_tools_1.solve1d2d)(W, bBuf, k2); + for (var i = 0; i < dim; ++i) + k2[i] -= kBuf[i]; + // 9) F3 = F(t + alpha3 * h, y + h * (alpha31 * k1 + alpha32 * k2)) <-- Fbuf + for (var i = 0; i < dim; ++i) + kBuf[i] = y[i] + h * (ALPHA_31 * k1[i] + ALPHA_32 * k2[i]); + f(t + ALPHA_3 * h, kBuf, fBuf); + // 10) kBuf = gamma31 / gamma * k1 + gamma32 / gamma * k2 + for (var i = 0; i < dim; ++i) + kBuf[i] = GAMMA_31_SCALED * k1[i] + GAMMA_32_SCALED * k2[i]; + // 11) k3 = W_inv * (F3 + kBuf + gamma3 * HT) - kBuf + for (var i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + kBuf[i] + GAMMA_3 * HT[i]; + if (toUseLU) + (0, lin_alg_tools_1.luSolve)(L, U, bBuf, luBuf, k3, dim); + else + (0, lin_alg_tools_1.solve1d2d)(W, bBuf, k3); + for (var i = 0; i < dim; ++i) + k3[i] -= kBuf[i]; + // 12) yNext = y + h * (b1 * k1 + b2 * k2 + b3 * k3) <-- yTemp + for (var i = 0; i < dim; ++i) + yTemp[i] = y[i] + h * (B_1 * k1[i] + B_2 * k2[i] + B_3 * k3[i]); + // 13) yErr = h * (r1 * k1 + r2 * k2 + r3 * k3) + for (var i = 0; i < dim; ++i) + yErr[i] = h * (R_1 * k1[i] + R_2 * k2[i] + R_3 * k3[i]); + // estimating error + errmax = 0; + for (var i = 0; i < dim; ++i) + errmax = (0, solver_defs_1.max)(errmax, (0, solver_defs_1.abs)(yErr[i] / yScale[i])); + errmax /= tolerance; + // processing the error obtained + if (errmax > 1) { + hTemp = solver_defs_1.SAFETY * h * Math.pow(errmax, solver_defs_1.PSHRNK); + h = (0, solver_defs_1.max)(hTemp, solver_defs_1.REDUCE_COEF * h); + tNew = t + h; + if (tNew == t) + throw new Error(solver_defs_1.ERROR_MSG.ROS3PRW_FAILS); + } + else { + if (errmax > solver_defs_1.ERR_CONTR) + hNext = solver_defs_1.SAFETY * h * Math.pow(errmax, solver_defs_1.PSGROW); + else + hNext = solver_defs_1.GROW_COEF * h; + t = t + h; + for (var i = 0; i < dim; ++i) + y[i] = yTemp[i]; + break; + } + } // while (true) + // compute lineraly interpolated results and store them in dataframe + while (timeDataframe < t) { + var cLeft = (t - timeDataframe) / (t - tPrev); + var cRight = 1.0 - cLeft; + tArr[index] = timeDataframe; + for (var j = 0; j < dim; ++j) + yArrs[j][index] = cRight * y[j] + cLeft * yPrev[j]; + timeDataframe += hDataframe; + ++index; + } + h = hNext; + tPrev = t; + for (var i = 0; i < dim; ++i) + yPrev[i] = y[i]; + } // while (flag) + // perform final callback actions + if (callback) + callback.onComputationsCompleted(); + // 3. solution at the point t1 + tArr[rowCount - 1] = t1; + for (var i = 0; i < dim; ++i) + yArrs[i][rowCount - 1] = y[i]; + return [tArr].concat(yArrs); +} // ros3pw +exports.ros3prw = ros3prw; diff --git a/libraries/diff-studio-utils/src/solver-tools/ros3prw-method.ts b/libraries/diff-studio-utils/src/solver-tools/ros3prw-method.ts new file mode 100644 index 0000000000..d51d3bdbfb --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/ros3prw-method.ts @@ -0,0 +1,289 @@ +/* The ROS3PRw method implementation + + References: + [1] https://doi.org/10.1016/j.cam.2015.03.010 */ + +import {ODEs, max, abs, SAFETY, PSHRNK, PSGROW, REDUCE_COEF, GROW_COEF, + ERR_CONTR, TINY, EPS, tDerivative, jacobian, ERROR_MSG} from './solver-defs'; +import {Callback} from './callbacks/callback-base'; +import {luDecomp, luSolve, solve1d2d} from './lin-alg-tools'; + +// The method specific constants (see Table 2 [1]) +const GAMMA = 0.78867513459481287; + +const GAMMA_21 = -2.3660254037844388; +const GAMMA_21_SCALED = GAMMA_21 / GAMMA; +const GAMMA_2 = GAMMA_21 + GAMMA; + +const GAMMA_31 = -0.86791218280355165; +const GAMMA_31_SCALED = GAMMA_31 / GAMMA; +const GAMMA_32 = -0.87306695894642317; +const GAMMA_32_SCALED = GAMMA_32 / GAMMA; +const GAMMA_3 = GAMMA_31 + GAMMA_32 + GAMMA; + +const ALPHA_21 = 2.3660254037844388; +const ALPHA_2 = ALPHA_21; +const ALPHA_31 = 0.5; +const ALPHA_32 = 0.76794919243112270; +const ALPHA_3 = ALPHA_31 + ALPHA_32; + +const B_1 = 0.50544867840851759; +const B_2 = -0.11571687603637559; +const B_3 = 0.610268197627858; + +const B_HAT_1 = 0.28973180237214197; +const B_HAT_2 = 0.10000000000000001; +const B_HAT_3 = 0.610268197627858; + +const R_1 = B_1 - B_HAT_1; +const R_2 = B_2 - B_HAT_2; +const R_3 = B_3 - B_HAT_3; + + +/** Solve initial value problem using the ROS3Pw method [5]. */ +export function ros3prw(odes: ODEs, callback?: Callback): Float64Array[] { + /** right-hand side of the IVP solved */ + const f = odes.func; + + // operating variables + const t0 = odes.arg.start; + const t1 = odes.arg.finish; + let h = odes.arg.step; + const hDataframe = h; + const tolerance = odes.tolerance; + + /** number of solution dataframe rows */ + const rowCount = Math.trunc((t1 - t0) / h) + 1; + + /** dimension of the problem */ + const dim = odes.initial.length; + const dimSquared = dim * dim; + + /** independent variable values */ + const tArr = new Float64Array(rowCount); + + /** arrays of solution values */ + const yArrs = Array(dim); + + for (let i = 0; i < dim; ++i) + yArrs[i] = new Float64Array(rowCount); + + // method routine + let timeDataframe = t0 + hDataframe; + let t = t0; + let tPrev = t0; + let flag = true; + let index = 1; + let errmax = 0; + let hTemp = 0; + let tNew = 0; + let hNext = 0.0; + + // 0 BUFFERS & TEMP STRUCTURES + + /** identity matrix */ + const I = new Float64Array(dim * dim); + + // compute identity matrix + for (let i = 0; i < dim; ++i) { + for (let j = 0; j < dim; ++j) + I[j + i * dim] = (i === j) ? 1 : 0; + } + + const y = new Float64Array(odes.initial); + const yPrev = new Float64Array(odes.initial); + const dydt = new Float64Array(dim); + const yScale = new Float64Array(dim); + const yTemp = new Float64Array(dim); + const yErr = new Float64Array(dim); + + const W = new Float64Array(dimSquared); + + const k1 = new Float64Array(dim); + const k2 = new Float64Array(dim); + const k3 = new Float64Array(dim); + const HT = new Float64Array(dim); + + const f0Buf = new Float64Array(dim); + const f1Buf = new Float64Array(dim); + + let hByGamma = 0; + const fBuf = new Float64Array(dim); + const kBuf = new Float64Array(dim); + + const L = new Float64Array(dimSquared); + const U = new Float64Array(dimSquared); + const luBuf = new Float64Array(dim); + const bBuf = new Float64Array(dim); + const toUseLU = dim > 2; + + // 1. SOLUTION AT THE POINT t0 + tArr[0] = t0; + for (let i = 0; i < dim; ++i) + yArrs[i][0] = y[i]; + + // 2. COMPUTE NUMERICAL SOLUTION FOR THE POINTS FROM THE INTERVAL (t0, t1) + while (flag) { + // compute derivative + f(t, y, dydt); + + // check whether to go on computations + if (callback) + callback.onIterationStart(); + + // compute scale vector + for (let i = 0; i < dim; ++i) + yScale[i] = abs(y[i]) + h * abs(dydt[i]) + TINY; + + // check end point + if (t + h > t1) { + h = t1 - t; + flag = false; + } + + // call of adaptive step the ROS3Pw [1] method + // computation of solution (y), time (t) and next step (hNext) + while (true) { + hByGamma = h * GAMMA; + + // one stage of the ROS3Pw method + + // 1) Jacobian & dF/dt matrices + jacobian(t, y, f, EPS, f0Buf, f1Buf, W); + tDerivative(t, y, f, EPS, f0Buf, f1Buf, HT); + + // 2) W & LU-decomposition + for (let i = 0; i < dimSquared; ++i) + W[i] = I[i] - hByGamma * W[i]; + + if (toUseLU) + luDecomp(W, L, U, dim); + + // 3) Scale dF/dt: HT = j * T + for (let i = 0; i < dim; ++i) + HT[i] *= h; + + // 4) F1 = F(t, y) <-- Fbuf + f(t, y, fBuf); + + // 5) k1 = W_inv * (F1 + gamma * HT) + for (let i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + GAMMA * HT[i]; + + if (toUseLU) + luSolve(L, U, bBuf, luBuf, k1, dim); + else + solve1d2d(W, bBuf, k1); + + // 6) F2 = F(t + alpha2 * h, y + alpha21 * k1) <-- Fbuf + for (let i = 0; i < dim; ++i) + kBuf[i] = y[i] + ALPHA_21 * h * k1[i]; + + f(t + ALPHA_2 * h, kBuf, fBuf); + + // 7) kBuf = gamma21 / gamma * k1 + for (let i = 0; i < dim; ++i) + kBuf[i] = GAMMA_21_SCALED * k1[i]; + + // 8) k2 = W_inv * [Fbuf + kBuf + gamma2 * HT] - kBuf + for (let i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + kBuf[i] + GAMMA_2 * HT[i]; + + if (toUseLU) + luSolve(L, U, bBuf, luBuf, k2, dim); + else + solve1d2d(W, bBuf, k2); + + for (let i = 0; i < dim; ++i) + k2[i] -= kBuf[i]; + + // 9) F3 = F(t + alpha3 * h, y + h * (alpha31 * k1 + alpha32 * k2)) <-- Fbuf + for (let i = 0; i < dim; ++i) + kBuf[i] = y[i] + h * (ALPHA_31 * k1[i] + ALPHA_32 * k2[i]); + + f(t + ALPHA_3 * h, kBuf, fBuf); + + // 10) kBuf = gamma31 / gamma * k1 + gamma32 / gamma * k2 + for (let i = 0; i < dim; ++i) + kBuf[i] = GAMMA_31_SCALED * k1[i] + GAMMA_32_SCALED * k2[i]; + + // 11) k3 = W_inv * (F3 + kBuf + gamma3 * HT) - kBuf + for (let i = 0; i < dim; ++i) + bBuf[i] = fBuf[i] + kBuf[i] + GAMMA_3 * HT[i]; + + if (toUseLU) + luSolve(L, U, bBuf, luBuf, k3, dim); + else + solve1d2d(W, bBuf, k3); + + for (let i = 0; i < dim; ++i) + k3[i] -= kBuf[i]; + + // 12) yNext = y + h * (b1 * k1 + b2 * k2 + b3 * k3) <-- yTemp + for (let i = 0; i < dim; ++i) + yTemp[i] = y[i] + h * (B_1 * k1[i] + B_2 * k2[i] + B_3 * k3[i]); + + // 13) yErr = h * (r1 * k1 + r2 * k2 + r3 * k3) + for (let i = 0; i < dim; ++i) + yErr[i] = h * (R_1 * k1[i] + R_2 * k2[i] + R_3 * k3[i]); + + // estimating error + errmax = 0; + for (let i = 0; i < dim; ++i) + errmax = max(errmax, abs(yErr[i] / yScale[i])); + errmax /= tolerance; + + // processing the error obtained + if (errmax > 1) { + hTemp = SAFETY * h * errmax**PSHRNK; + h = max(hTemp, REDUCE_COEF * h); + tNew = t + h; + if (tNew == t) + throw new Error(ERROR_MSG.ROS3PRW_FAILS); + } else { + if (errmax > ERR_CONTR) + hNext = SAFETY * h * errmax**PSGROW; + else + hNext = GROW_COEF * h; + t = t + h; + + for (let i = 0; i < dim; ++i) + y[i] = yTemp[i]; + + break; + } + } // while (true) + + // compute lineraly interpolated results and store them in dataframe + while (timeDataframe < t) { + const cLeft = (t - timeDataframe) / (t - tPrev); + const cRight = 1.0 - cLeft; + + tArr[index] = timeDataframe; + + for (let j = 0; j < dim; ++j) + yArrs[j][index] = cRight * y[j] + cLeft * yPrev[j]; + + timeDataframe += hDataframe; + ++index; + } + + h = hNext; + tPrev = t; + + for (let i = 0; i < dim; ++i) + yPrev[i] = y[i]; + } // while (flag) + + // perform final callback actions + if (callback) + callback.onComputationsCompleted(); + + // 3. solution at the point t1 + tArr[rowCount - 1] = t1; + + for (let i = 0; i < dim; ++i) + yArrs[i][rowCount - 1] = y[i]; + + return [tArr].concat(yArrs); +} // ros3pw diff --git a/libraries/diff-studio-utils/src/solver-tools/solver-defs.js b/libraries/diff-studio-utils/src/solver-tools/solver-defs.js new file mode 100644 index 0000000000..4824adfb52 --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/solver-defs.js @@ -0,0 +1,83 @@ +"use strict"; +// Solver definitions +var __extends = (this && this.__extends) || (function () { + var extendStatics = function (d, b) { + extendStatics = Object.setPrototypeOf || + ({ __proto__: [] } instanceof Array && function (d, b) { d.__proto__ = b; }) || + function (d, b) { for (var p in b) if (Object.prototype.hasOwnProperty.call(b, p)) d[p] = b[p]; }; + return extendStatics(d, b); + }; + return function (d, b) { + if (typeof b !== "function" && b !== null) + throw new TypeError("Class extends value " + String(b) + " is not a constructor or null"); + extendStatics(d, b); + function __() { this.constructor = d; } + d.prototype = b === null ? Object.create(b) : (__.prototype = b.prototype, new __()); + }; +})(); +Object.defineProperty(exports, "__esModule", { value: true }); +exports.DEFAULT_OPTIONS = exports.CallbackAction = exports.ERROR_MSG = exports.jacobian = exports.tDerivative = exports.EPS = exports.TINY = exports.ERR_CONTR = exports.GROW_COEF = exports.REDUCE_COEF = exports.PSGROW = exports.PSHRNK = exports.SAFETY = exports.min = exports.max = exports.abs = void 0; +/** The abs function */ +var abs = function (x) { return (x > 0) ? x : -x; }; +exports.abs = abs; +/** The max function */ +var max = function (x, y) { return (x > y) ? x : y; }; +exports.max = max; +/** The max function */ +var min = function (x, y) { return (x < y) ? x : y; }; +exports.min = min; +/** Routine constants of adaptive step method */ +exports.SAFETY = 0.9; +exports.PSHRNK = -0.25; +exports.PSGROW = -0.2; +exports.REDUCE_COEF = 0.25; +exports.GROW_COEF = 4.0; +exports.ERR_CONTR = 1.89e-4; +/** Misc */ +exports.TINY = 1e-20; +exports.EPS = 1.0e-10; +/** Returns derivative with respect to t. */ +function tDerivative(t, y, f, eps, f0Buf, f1Buf, output) { + var size = y.length; + f(t, y, f0Buf); + f(t + eps, y, f1Buf); + for (var i = 0; i < size; ++i) + output[i] = (f1Buf[i] - f0Buf[i]) / eps; +} +exports.tDerivative = tDerivative; +/** Returns Jacobian. */ +function jacobian(t, y, f, eps, f0Buf, f1Buf, output) { + var size = y.length; + f(t, y, f0Buf); + for (var j = 0; j < size; ++j) { + y[j] += eps; + f(t, y, f1Buf); + for (var i = 0; i < size; ++i) + output[j + i * size] = (f1Buf[i] - f0Buf[i]) / eps; + y[j] -= eps; + } +} +exports.jacobian = jacobian; +/** Error messeges */ +var ERROR_MSG; +(function (ERROR_MSG) { + ERROR_MSG["MRT_FAILS"] = "The modified Rosenbrock triple method fails"; + ERROR_MSG["ROS3PRW_FAILS"] = "The ROS3PRw method fails"; + ERROR_MSG["ROS34PRW_FAILS"] = "The ROS34PRw method fails"; +})(ERROR_MSG = exports.ERROR_MSG || (exports.ERROR_MSG = {})); +; +/** Callback action */ +var CallbackAction = /** @class */ (function (_super) { + __extends(CallbackAction, _super); + function CallbackAction(msg) { + return _super.call(this, msg) || this; + } + return CallbackAction; +}(Error)); +exports.CallbackAction = CallbackAction; +/** Default options of the solver */ +var DEFAULT_OPTIONS; +(function (DEFAULT_OPTIONS) { + DEFAULT_OPTIONS["SCRIPTING"] = "{maxIterations: 1}"; + DEFAULT_OPTIONS["NO_CHECKS"] = "{ }"; +})(DEFAULT_OPTIONS = exports.DEFAULT_OPTIONS || (exports.DEFAULT_OPTIONS = {})); diff --git a/libraries/diff-studio-utils/src/solver-tools/solver-defs.ts b/libraries/diff-studio-utils/src/solver-tools/solver-defs.ts new file mode 100644 index 0000000000..aec30fe19e --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/solver-defs.ts @@ -0,0 +1,95 @@ +// Solver definitions + +/** Right-hand side of IVP */ +type Func = (t: number, y: Float64Array, output: Float64Array) => void; + +/** Solver optional settings */ +export type SolverOptions = { + maxIterations: number, + maxTimeMs: number, + method: string, +}; + +/** Initial Value Problem (IVP) type */ +export type ODEs = { + name: string, + arg: { + name: string, + start: number, + finish: number, + step: number + }, + initial: number[] | Float32Array | Float64Array | Int32Array, + func: Func, + tolerance: number, + solutionColNames: string[], +}; + +/** The abs function */ +export const abs = (x: number) => (x > 0) ? x : -x; + +/** The max function */ +export const max = (x: number, y: number) => (x > y) ? x : y; + +/** The max function */ +export const min = (x: number, y: number) => (x < y) ? x : y; + +/** Routine constants of adaptive step method */ +export const SAFETY = 0.9; +export const PSHRNK = -0.25; +export const PSGROW = -0.2; +export const REDUCE_COEF = 0.25; +export const GROW_COEF = 4.0; +export const ERR_CONTR = 1.89e-4; + +/** Misc */ +export const TINY = 1e-20; +export const EPS = 1.0e-10; + +/** Returns derivative with respect to t. */ +export function tDerivative(t: number, y: Float64Array, f: Func, eps: number, + f0Buf: Float64Array, f1Buf: Float64Array, output: Float64Array): void { + const size = y.length; + f(t, y, f0Buf); + f(t + eps, y, f1Buf); + + for (let i = 0; i < size; ++i) + output[i] = (f1Buf[i] - f0Buf[i]) / eps; +} + +/** Returns Jacobian. */ +export function jacobian(t: number, y: Float64Array, f: Func, eps: number, + f0Buf: Float64Array, f1Buf: Float64Array, output: Float64Array): void { + const size = y.length; + f(t, y, f0Buf); + + for (let j = 0; j < size; ++j) { + y[j] += eps; + f(t, y, f1Buf); + + for (let i = 0; i < size; ++i) + output[j + i * size] = (f1Buf[i] - f0Buf[i]) / eps; + + y[j] -= eps; + } +} + +/** Error messeges */ +export enum ERROR_MSG { + MRT_FAILS = 'The modified Rosenbrock triple method fails', + ROS3PRW_FAILS = 'The ROS3PRw method fails', + ROS34PRW_FAILS = 'The ROS34PRw method fails', +}; + +/** Callback action */ +export class CallbackAction extends Error { + constructor(msg: string) { + super(msg); + } +} + +/** Default options of the solver */ +export enum DEFAULT_OPTIONS { + SCRIPTING = '{maxIterations: 1}', + NO_CHECKS = '{ }', +} diff --git a/libraries/diff-studio-utils/src/solver-tools/tests.js b/libraries/diff-studio-utils/src/solver-tools/tests.js new file mode 100644 index 0000000000..da0e4b1d91 --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/tests.js @@ -0,0 +1,392 @@ +"use strict"; +Object.defineProperty(exports, "__esModule", { value: true }); +var mrt_method_1 = require("./mrt-method"); +var ros34prw_method_1 = require("./ros34prw-method"); +var ros3prw_method_1 = require("./ros3prw-method"); +var methods = new Map([ + ['MRT', mrt_method_1.mrt], + ['ROS3PRw', ros3prw_method_1.ros3prw], + ['ROS34PRw', ros34prw_method_1.ros34prw], +]); +/** Robertson chemical reaction, updated version (see https://archimede.uniba.it/~testset/problems/rober.php) */ +var robertson = { + name: 'Robertson', + arg: { name: 't', start: 0, finish: 10e11, step: 2.5e7 }, + initial: [1, 0, 0], + func: function (t, y, output) { + output[0] = -0.04 * y[0] + 1e4 * y[1] * y[2]; + output[1] = 0.04 * y[0] - 1e4 * y[1] * y[2] - 3e7 * Math.pow(y[1], 2); + output[2] = 3e7 * Math.pow(y[1], 2); + }, + tolerance: 1e-7, + solutionColNames: ['A', 'B', 'C'], +}; +/** High Irradiance Responses of photomorphogenesis (see https://archimede.uniba.it/~testset/problems/hires.php) */ +var hires = { + name: 'HIRES', + arg: { name: 't', start: 0, finish: 321.8122, step: 0.01 }, + initial: [1, 0, 0, 0, 0, 0, 0, 0.0057], + func: function (t, y, output) { + // extract function values + var y1 = y[0]; + var y2 = y[1]; + var y3 = y[2]; + var y4 = y[3]; + var y5 = y[4]; + var y6 = y[5]; + var y7 = y[6]; + var y8 = y[7]; + // compute output + output[0] = -1.71 * y1 + 0.43 * y2 + 8.32 * y3 + 0.0007; + output[1] = 1.71 * y1 - 8.75 * y2; + output[2] = -10.03 * y3 + 0.43 * y4 + 0.035 * y5; + output[3] = 8.32 * y2 + 1.71 * y3 - 1.12 * y4; + output[4] = -1.745 * y5 + 0.43 * y6 + 0.43 * y7; + output[5] = -280 * y6 * y8 + 0.69 * y4 + 1.71 * y5 - 0.43 * y6 + 0.69 * y7; + output[6] = 280 * y6 * y8 - 1.81 * y7; + output[7] = -280 * y6 * y8 + 1.81 * y7; + }, + tolerance: 1e-10, + solutionColNames: ['y1', 'y2', 'y3', 'y4', 'y5', 'y6', 'y7', 'y8'], +}; // hires +/** Van der Pol oscillator (see https://archimede.uniba.it/~testset/problems/vdpol.php) */ +var vanDerPol = { + name: 'van der Pol', + arg: { name: 't', start: 0, finish: 2000, step: 0.1 }, + initial: [-1, 1], + func: function (t, y, output) { + output[0] = y[1]; + output[1] = -y[0] + 1000 * (1 - y[0] * y[0]) * y[1]; + }, + tolerance: 1e-12, + solutionColNames: ['x1', 'x2'], +}; +/** The OREGO model (see https://archimede.uniba.it/~testset/report/orego.pdf) */ +var orego = { + name: 'OREGO', + arg: { name: 't', start: 0, finish: 360, step: 0.01 }, + initial: [1, 2, 3], + func: function (t, y, output) { + // extract function values + var y1 = y[0]; + var y2 = y[1]; + var y3 = y[2]; + // compute output + output[0] = 77.27 * (y2 - y1 * y2 + y1 - 0.000008375 * y1 * y1); + output[1] = 1 / 77.27 * (-y2 - y1 * y2 + y3); + output[2] = 0.161 * (y1 - y3); + }, + tolerance: 1e-8, + solutionColNames: ['y1', 'y2', 'y3'], +}; +/** Kintetic constants for the E5 model */ +var E5; +(function (E5) { + E5[E5["K1"] = 7.89e-10] = "K1"; + E5[E5["K2"] = 1130000000] = "K2"; + E5[E5["K3"] = 11000000] = "K3"; + E5[E5["K4"] = 1130] = "K4"; +})(E5 || (E5 = {})); +; +/** The E5 model (chemical pyrolysis: https://archimede.uniba.it/~testset/report/e5.pdf) */ +var e5 = { + name: 'E5', + arg: { name: 't', start: 0, finish: 1e13, step: 2.5e8 }, + initial: [0.00176, 0, 0, 0], + func: function (t, y, output) { + // extract function values + var y1 = y[0]; + var y2 = y[1]; + var y3 = y[2]; + var y4 = y[3]; + // compute output + output[0] = -E5.K1 * y1 - E5.K3 * y1 * y3; + output[1] = E5.K1 * y1 - E5.K2 * y2 * y3; + output[2] = E5.K1 * y1 - E5.K2 * y2 * y3 - E5.K3 * y1 * y3 + E5.K4 * y4; + output[3] = E5.K3 * y1 * y3 - E5.K4 * y4; + }, + tolerance: 1e-6, + solutionColNames: ['y1', 'y2', 'y3', 'y4'], +}; +/** Kintetic constants for the Pollution model */ +var POL; +(function (POL) { + POL[POL["K1"] = 0.35] = "K1"; + POL[POL["K2"] = 26.6] = "K2"; + POL[POL["K3"] = 12300] = "K3"; + POL[POL["K4"] = 0.00086] = "K4"; + POL[POL["K5"] = 0.00082] = "K5"; + POL[POL["K6"] = 15000] = "K6"; + POL[POL["K7"] = 0.00013] = "K7"; + POL[POL["K8"] = 24000] = "K8"; + POL[POL["K9"] = 16500] = "K9"; + POL[POL["K10"] = 9000] = "K10"; + POL[POL["K11"] = 0.022] = "K11"; + POL[POL["K12"] = 12000] = "K12"; + POL[POL["K13"] = 1.88] = "K13"; + POL[POL["K14"] = 16300] = "K14"; + POL[POL["K15"] = 4800000] = "K15"; + POL[POL["K16"] = 0.00035] = "K16"; + POL[POL["K17"] = 0.0175] = "K17"; + POL[POL["K18"] = 100000000] = "K18"; + POL[POL["K19"] = 444000000000] = "K19"; + POL[POL["K20"] = 1240] = "K20"; + POL[POL["K21"] = 2.1] = "K21"; + POL[POL["K22"] = 5.78] = "K22"; + POL[POL["K23"] = 0.0474] = "K23"; + POL[POL["K24"] = 1780] = "K24"; + POL[POL["K25"] = 3.12] = "K25"; +})(POL || (POL = {})); +; // POL +/** The chemical reaction part of the air pollution model (https://archimede.uniba.it/~testset/report/pollu.pdf) */ +var pollution = { + name: 'Pollution', + arg: { name: 't', start: 0, finish: 60, step: 0.002 }, + initial: [0, 0.2, 0, 0.04, 0, 0, 0.1, 0.3, 0.01, 0, 0, 0, 0, 0, 0, 0, 0.007, 0, 0, 0], + func: function (t, y, output) { + // extract function values + var y1 = y[0]; + var y2 = y[1]; + var y3 = y[2]; + var y4 = y[3]; + var y5 = y[4]; + var y6 = y[5]; + var y7 = y[6]; + var y8 = y[7]; + var y9 = y[8]; + var y10 = y[9]; + var y11 = y[10]; + var y12 = y[11]; + var y13 = y[12]; + var y14 = y[13]; + var y15 = y[14]; + var y16 = y[15]; + var y17 = y[16]; + var y18 = y[17]; + var y19 = y[18]; + var y20 = y[19]; + // evaluate expressions + var r1 = POL.K1 * y1; + var r2 = POL.K2 * y2 * y4; + var r3 = POL.K3 * y5 * y2; + var r4 = POL.K4 * y7; + var r5 = POL.K5 * y7; + var r6 = POL.K6 * y7 * y6; + var r7 = POL.K7 * y9; + var r8 = POL.K8 * y9 * y6; + var r9 = POL.K9 * y11 * y2; + var r10 = POL.K10 * y11 * y1; + var r11 = POL.K11 * y13; + var r12 = POL.K12 * y10 * y2; + var r13 = POL.K13 * y14; + var r14 = POL.K14 * y1 * y6; + var r15 = POL.K15 * y3; + var r16 = POL.K16 * y4; + var r17 = POL.K17 * y4; + var r18 = POL.K18 * y16; + var r19 = POL.K19 * y16; + var r20 = POL.K20 * y17 * y6; + var r21 = POL.K21 * y19; + var r22 = POL.K22 * y19; + var r23 = POL.K23 * y1 * y4; + var r24 = POL.K24 * y19 * y1; + var r25 = POL.K25 * y20; + // compute output + output[0] = -(r1 + r10 + r14 + r23 + r24) + (r2 + r3 + r9 + r11 + r12 + r22 + r25); + output[1] = -r2 - r3 - r9 - r12 + r1 + r21; + output[2] = -r15 + r1 + r17 + r19 + r22; + output[3] = -r2 - r16 - r17 - r23 + r15; + output[4] = -r3 + 2 * r4 + r6 + r7 + r13 + r20; + output[5] = -r6 - r8 - r14 - r20 + r3 + 2 * r18; + output[6] = -r4 - r5 - r6 + r13; + output[7] = r4 + r5 + r6 + r7; + output[8] = -r7 - r8; + output[9] = -r12 + r7 + r9; + output[10] = -r9 - r10 + r8 + r11; + output[11] = r9; + output[12] = -r11 + r10; + output[13] = -r13 + r12; + output[14] = r14; + output[15] = -r18 - r19 + r16; + output[16] = -r20; + output[17] = r20; + output[18] = -r21 - r22 - r24 + r23 + r25; + output[19] = -r25 + r24; + }, + tolerance: 1e-6, + solutionColNames: ['y1', 'y2', 'y3', 'y4', 'y5', 'y6', 'y7', 'y8', 'y9', 'y10', 'y11', + 'y12', 'y13', 'y14', 'y15', 'y16', 'y17', 'y18', 'y19', 'y20'], +}; // pollution +/** Problems for testing solvers' performance */ +var performanceProblems = new Map([ + [robertson.name, robertson], + [hires.name, hires], + [vanDerPol.name, vanDerPol], + [orego.name, orego], + [e5.name, e5], + [pollution.name, pollution], +]); +/** Non-stiff 1D test problem: equation (see [1], p. 736) */ +var nonStiff1D = { + name: 'non-stiff 1D', + arg: { name: 't', start: 0, finish: 4, step: 0.01 }, + initial: [2], + func: function (t, y, output) { + output[0] = 4 * Math.exp(0.8 * t) - 0.5 * y[0]; + }, + tolerance: 0.00001, + solutionColNames: ['y'], +}; +/** Non-stiff 1D test problem: exact solution (see [1], p. 736) */ +var exactNonStiff1D = function (t) { + return new Float64Array([ + (Math.exp(0.8 * t) - Math.exp(-0.5 * t)) * 4 / 1.3 + 2 * Math.exp(-0.5 * t), + ]); +}; +/** Non-stiff 2D test problem: equation */ +var nonStiff2D = { + name: 'non-stiff 2D', + arg: { name: 't', start: 0, finish: 4, step: 0.01 }, + initial: [1, 1], + func: function (t, y, output) { + output[0] = y[0] + y[1]; + output[1] = y[1] - y[0]; + }, + tolerance: 0.000000001, + solutionColNames: ['x', 'y'], +}; +/** Non-stiff 2D test problem: exact solution */ +var exactNonStiff2D = function (t) { + return new Float64Array([ + Math.exp(t) * (Math.cos(t) + Math.sin(t)), + Math.exp(t) * (Math.cos(t) - Math.sin(t)), + ]); +}; +/** Non-stiff 3D test problem: equations */ +var nonStiff3D = { + name: 'non-stiff 3D', + arg: { name: 't', start: 0, finish: 2, step: 0.001 }, + initial: [0.3, -0.8, 0], + func: function (t, y, output) { + output[0] = 5 * y[0] + 2 * y[1] + Math.sin(t); + output[1] = -4 * y[0] - y[1] + Math.exp(2 * t); + output[2] = 5 * Math.pow(t, 4) - 3 * Math.pow(t, 2) + 2 * t; + }, + tolerance: 0.00000001, + solutionColNames: ['x', 'y', 'z'], +}; +/** Non-stiff 3D test problem: exact solution */ +var exactNonStiff3D = function (t) { + var e1 = Math.exp(t); + var e2 = Math.exp(2 * t); + var e3 = Math.exp(3 * t); + var c = Math.cos(t); + var s = Math.sin(t); + return new Float64Array([ + e1 + e3 + 0.1 * (3 * c - s) - 2 * e2, + -2 * e1 - e3 - 0.2 * (2 * s + 4 * c) + 3 * e2, + Math.pow(t, 5) - Math.pow(t, 3) + Math.pow(t, 2), + ]); +}; +/** Stiff 1D test problem: equation (see [1], p. 767) */ +var stiff1D = { + name: 'stiff 1D', + arg: { name: 't', start: 0, finish: 4, step: 0.01 }, + initial: [0], + func: function (t, y, output) { + output[0] = -1000 * y[0] + 3000 - 2000 * Math.exp(-t); + }, + tolerance: 0.0000005, + solutionColNames: ['y'], +}; +/** Stiff 1D test problem: exact solution (see [1], p. 767) */ +var exactStiff1D = function (t) { + return new Float64Array([ + 3 - 0.998 * Math.exp(-1000 * t) - 2.002 * Math.exp(-t), + ]); +}; +/** Stiff 2D test problem: equations (see [1], p. 770) */ +var stiff2D = { + name: 'stiff 2D', + arg: { name: 't', start: 0, finish: 4, step: 0.01 }, + initial: [52.29, 83.82], + func: function (t, y, output) { + output[0] = -5 * y[0] + 3 * y[1]; + output[1] = 100 * y[0] - 301 * y[1]; + }, + tolerance: 0.0000005, + solutionColNames: ['x', 'y'], +}; +/** Stiff 2D test problem: exact solution (see [1], p. 770) */ +var exactStiff2D = function (t) { + return new Float64Array([ + 52.96 * Math.exp(-3.9899 * t) - 0.67 * Math.exp(-302.0101 * t), + 17.83 * Math.exp(-3.9899 * t) + 65.99 * Math.exp(-302.0101 * t), + ]); +}; +/** Stiff 3D test problem: equations */ +var stiff3D = { + name: 'stiff 3D', + arg: { name: 't', start: 0, finish: 4, step: 0.01 }, + initial: [2, 1, 0], + func: function (t, y, output) { + output[0] = -100 * y[0] + 100 * y[1] + y[2]; + output[1] = y[2]; + output[2] = -y[1]; + }, + tolerance: 0.00000001, + solutionColNames: ['x', 'y', 'z'], +}; +/** Stiff 3D test problem: exact solution */ +var exactStiff3D = function (t) { + return new Float64Array([ + Math.exp(-100 * t) + Math.cos(t), + Math.cos(t), + -Math.sin(t), + ]); +}; +var correctnessProblems = [ + { odes: nonStiff1D, exact: exactNonStiff1D }, + { odes: nonStiff2D, exact: exactNonStiff2D }, + { odes: nonStiff3D, exact: exactNonStiff3D }, + { odes: stiff1D, exact: exactStiff1D }, + { odes: stiff2D, exact: exactStiff2D }, + { odes: stiff3D, exact: exactStiff3D }, +]; +/** Return numerical solution error: maximum absolute deviation between approximate & exact solutions */ +function getError(method, corProb) { + var error = 0; + // Get numerical solution + var approxSolution = method(corProb.odes); + var exact = corProb.exact; + var arg = approxSolution[0]; + var pointsCount = arg.length; + var funcsCount = approxSolution.length - 1; + // Compute error + for (var i = 0; i < pointsCount; ++i) { + var exactSolution = exact(arg[i]); + for (var j = 0; j < funcsCount; ++j) + error = Math.max(error, Math.abs(exactSolution[j] - approxSolution[j + 1][i])); + } + return error; +} +console.log('Performance:\n'); +methods.forEach(function (method, name) { + console.log(' ', name); + performanceProblems.forEach(function (odes, name) { + var start = Date.now(); + method(odes); + var finish = Date.now(); + console.log(" ".concat(name, ": ").concat(finish - start, " ms.")); + }); + console.log(); +}); +console.log('Correctness:\n'); +methods.forEach(function (method, name) { + console.log(' ', name); + correctnessProblems.forEach(function (problem) { + var error = getError(method, problem); + console.log(" ".concat(problem.odes.name, ": ").concat(error)); + }); + console.log(); +}); diff --git a/libraries/diff-studio-utils/src/solver-tools/tests.ts b/libraries/diff-studio-utils/src/solver-tools/tests.ts new file mode 100644 index 0000000000..0de6c1f046 --- /dev/null +++ b/libraries/diff-studio-utils/src/solver-tools/tests.ts @@ -0,0 +1,437 @@ +import {mrt} from './mrt-method'; +import { ros34prw } from './ros34prw-method'; +import {ros3prw} from './ros3prw-method'; +import {ODEs} from './solver-defs'; + +const methods = new Map([ + ['MRT', mrt], + ['ROS3PRw', ros3prw], + ['ROS34PRw', ros34prw], +]); + +/** Robertson chemical reaction, updated version (see https://archimede.uniba.it/~testset/problems/rober.php) */ +const robertson = { + name: 'Robertson', + arg: {name: 't', start: 0, finish: 10e11, step: 2.5e7}, + initial: [1, 0, 0], + func: (t: number, y: Float64Array, output: Float64Array) => { + output[0] = -0.04 * y[0] + 1e4 * y[1] * y[2]; + output[1] = 0.04 * y[0] - 1e4 * y[1] * y[2] - 3e7 * y[1]**2; + output[2] = 3e7 * y[1]**2; + }, + tolerance: 1e-7, + solutionColNames: ['A', 'B', 'C'], + }; + + /** High Irradiance Responses of photomorphogenesis (see https://archimede.uniba.it/~testset/problems/hires.php) */ + const hires = { + name: 'HIRES', + arg: {name: 't', start: 0, finish: 321.8122, step: 0.01}, + initial: [1, 0, 0, 0, 0, 0, 0, 0.0057], + func: (t: number, y: Float64Array, output: Float64Array) => { + // extract function values + const y1 = y[0]; + const y2 = y[1]; + const y3 = y[2]; + const y4 = y[3]; + const y5 = y[4]; + const y6 = y[5]; + const y7 = y[6]; + const y8 = y[7]; + + // compute output + output[0] = -1.71 * y1 + 0.43 * y2 + 8.32 * y3 + 0.0007; + output[1] = 1.71 * y1 - 8.75 * y2; + output[2] = -10.03 * y3 + 0.43 * y4 + 0.035 * y5; + output[3] = 8.32 * y2 + 1.71 * y3 - 1.12 * y4; + output[4] = -1.745 * y5 + 0.43 * y6 + 0.43 * y7; + output[5] = -280 * y6 * y8 + 0.69 * y4 + 1.71 * y5 - 0.43 * y6 + 0.69 * y7; + output[6] = 280 * y6 * y8 - 1.81 * y7; + output[7] = -280 * y6 * y8 + 1.81 * y7; + }, + tolerance: 1e-10, + solutionColNames: ['y1', 'y2', 'y3', 'y4', 'y5', 'y6', 'y7', 'y8'], + }; // hires + + /** Van der Pol oscillator (see https://archimede.uniba.it/~testset/problems/vdpol.php) */ + const vanDerPol = { + name: 'van der Pol', + arg: {name: 't', start: 0, finish: 2000, step: 0.1}, + initial: [-1, 1], + func: (t: number, y: Float64Array, output: Float64Array) => { + output[0] = y[1]; + output[1] = -y[0] + 1000 * (1 - y[0] * y[0]) * y[1]; + }, + tolerance: 1e-12, + solutionColNames: ['x1', 'x2'], + }; + + /** The OREGO model (see https://archimede.uniba.it/~testset/report/orego.pdf) */ + const orego = { + name: 'OREGO', + arg: {name: 't', start: 0, finish: 360, step: 0.01}, + initial: [1, 2, 3], + func: (t: number, y: Float64Array, output: Float64Array) => { + // extract function values + const y1 = y[0]; + const y2 = y[1]; + const y3 = y[2]; + + // compute output + output[0] = 77.27 * (y2 - y1 * y2 + y1 - 0.000008375 * y1 * y1 ); + output[1] = 1 / 77.27 * (-y2 - y1 * y2 + y3); + output[2] = 0.161 * (y1 - y3); + }, + tolerance: 1e-8, + solutionColNames: ['y1', 'y2', 'y3'], + }; + + /** Kintetic constants for the E5 model */ + enum E5 { + K1 = 7.89e-10, + K2 = 1.13e9, + K3 = 1.1e7, + K4 = 1.13e3, + }; + + /** The E5 model (chemical pyrolysis: https://archimede.uniba.it/~testset/report/e5.pdf) */ + const e5 = { + name: 'E5', + arg: {name: 't', start: 0, finish: 1e13, step: 2.5e8}, + initial: [0.00176, 0, 0, 0], + func: (t: number, y: Float64Array, output: Float64Array) => { + // extract function values + const y1 = y[0]; + const y2 = y[1]; + const y3 = y[2]; + const y4 = y[3]; + + // compute output + output[0] = -E5.K1 * y1 - E5.K3 * y1 * y3; + output[1] = E5.K1 * y1 - E5.K2 * y2 * y3; + output[2] = E5.K1 * y1 - E5.K2 * y2 * y3 - E5.K3 * y1 * y3 + E5.K4 * y4; + output[3] = E5.K3 * y1 * y3 - E5.K4 * y4; + }, + tolerance: 1e-6, + solutionColNames: ['y1', 'y2', 'y3', 'y4'], + }; + + /** Kintetic constants for the Pollution model */ + enum POL { + K1 = 0.35, + K2 = 26.6, + K3 = 1.23e4, + K4 = 8.6e-4, + K5 = 8.2e-4, + K6 = 1.5e4, + K7 = 1.3e-4, + K8 = 2.4e4, + K9 = 1.65e4, + K10 = 9e3, + K11 = 0.022, + K12 = 1.2e4, + K13 = 1.88, + K14 = 1.63e4, + K15 = 4.8e6, + K16 = 3.5e-4, + K17 = 0.0175, + K18 = 1e8, + K19 = 4.44e11, + K20 = 1240, + K21 = 2.1, + K22 = 5.78, + K23 = 0.0474, + K24 = 1780, + K25 = 3.12, + }; // POL + + /** The chemical reaction part of the air pollution model (https://archimede.uniba.it/~testset/report/pollu.pdf) */ + const pollution = { + name: 'Pollution', + arg: {name: 't', start: 0, finish: 60, step: 0.002}, + initial: [0, 0.2, 0, 0.04, 0, 0, 0.1, 0.3, 0.01, 0, 0, 0, 0, 0, 0, 0, 0.007, 0, 0, 0], + func: (t: number, y: Float64Array, output: Float64Array) => { + // extract function values + const y1 = y[0]; + const y2 = y[1]; + const y3 = y[2]; + const y4 = y[3]; + const y5 = y[4]; + const y6 = y[5]; + const y7 = y[6]; + const y8 = y[7]; + const y9 = y[8]; + const y10 = y[9]; + const y11 = y[10]; + const y12 = y[11]; + const y13 = y[12]; + const y14 = y[13]; + const y15 = y[14]; + const y16 = y[15]; + const y17 = y[16]; + const y18 = y[17]; + const y19 = y[18]; + const y20 = y[19]; + + // evaluate expressions + const r1 = POL.K1 * y1; + const r2 = POL.K2 * y2 * y4; + const r3 = POL.K3 * y5 * y2; + const r4 = POL.K4 * y7; + const r5 = POL.K5 * y7; + const r6 = POL.K6 * y7 * y6; + const r7 = POL.K7 * y9; + const r8 = POL.K8 * y9 * y6; + const r9 = POL.K9 * y11 * y2; + const r10 = POL.K10 * y11 * y1; + const r11 = POL.K11 * y13; + const r12 = POL.K12 * y10 * y2; + const r13 = POL.K13 * y14; + const r14 = POL.K14 * y1 * y6; + const r15 = POL.K15 * y3; + const r16 = POL.K16 * y4; + const r17 = POL.K17 * y4; + const r18 = POL.K18 * y16; + const r19 = POL.K19 * y16; + const r20 = POL.K20 * y17 * y6; + const r21 = POL.K21 * y19; + const r22 = POL.K22 * y19; + const r23 = POL.K23 * y1 * y4; + const r24 = POL.K24 * y19 * y1; + const r25 = POL.K25 * y20; + + // compute output + output[0] = -(r1 + r10 + r14 + r23 + r24) + (r2 + r3 + r9 + r11 + r12 + r22 + r25); + output[1] = -r2 - r3 - r9 - r12 + r1 + r21; + output[2] = -r15 + r1 + r17 + r19 + r22; + output[3] = -r2 - r16 - r17 - r23 + r15; + output[4] = -r3 +2 * r4 + r6 +r7 +r13 + r20; + output[5] = -r6 - r8 - r14 - r20 + r3 + 2 * r18; + output[6] = -r4 - r5 - r6 + r13; + output[7] = r4 + r5 + r6 + r7; + output[8] = -r7 - r8; + output[9] = -r12 +r7 + r9; + output[10] = -r9 - r10 + r8 + r11; + output[11] = r9; + output[12] = -r11 + r10; + output[13] = -r13 + r12; + output[14] = r14; + output[15] = -r18 - r19 + r16; + output[16] = -r20; + output[17] = r20; + output[18] = -r21 - r22 - r24 + r23 + r25; + output[19] = -r25 + r24; + }, + tolerance: 1e-6, + solutionColNames: ['y1', 'y2', 'y3', 'y4', 'y5', 'y6', 'y7', 'y8', 'y9', 'y10', 'y11', + 'y12', 'y13', 'y14', 'y15', 'y16', 'y17', 'y18', 'y19', 'y20'], + }; // pollution + + /** Problems for testing solvers' performance */ + const performanceProblems = new Map([ + [robertson.name, robertson], + [hires.name, hires], + [vanDerPol.name, vanDerPol], + [orego.name, orego], + [e5.name, e5], + [pollution.name, pollution], + ]); + + /** Non-stiff 1D test problem: equation (see [1], p. 736) */ +const nonStiff1D = { + name: 'non-stiff 1D', + arg: {name: 't', start: 0, finish: 4, step: 0.01}, + initial: [2], + func: (t: number, y: Float64Array, output: Float64Array) => { + output[0] = 4 * Math.exp(0.8 * t) - 0.5 * y[0]; + }, + tolerance: 0.00001, + solutionColNames: ['y'], + }; + + /** Non-stiff 1D test problem: exact solution (see [1], p. 736) */ + const exactNonStiff1D = (t: number) => { + return new Float64Array([ + (Math.exp(0.8 * t) - Math.exp(-0.5 * t)) * 4 / 1.3 + 2 * Math.exp(-0.5 * t), + ]); + }; + + /** Non-stiff 2D test problem: equation */ + const nonStiff2D = { + name: 'non-stiff 2D', + arg: {name: 't', start: 0, finish: 4, step: 0.01}, + initial: [1, 1], + func: (t: number, y: Float64Array, output: Float64Array) => { + output[0] = y[0] + y[1]; + output[1] = y[1] - y[0]; + }, + tolerance: 0.000000001, + solutionColNames: ['x', 'y'], + }; + + /** Non-stiff 2D test problem: exact solution */ + const exactNonStiff2D = (t: number) => { + return new Float64Array([ + Math.exp(t) * (Math.cos(t) + Math.sin(t)), + Math.exp(t) * (Math.cos(t) - Math.sin(t)), + ]); + }; + + /** Non-stiff 3D test problem: equations */ + const nonStiff3D = { + name: 'non-stiff 3D', + arg: {name: 't', start: 0, finish: 2, step: 0.001}, + initial: [0.3, -0.8, 0], + func: (t: number, y: Float64Array, output: Float64Array) => { + output[0] = 5 * y[0] + 2 * y[1] + Math.sin(t); + output[1] = -4 * y[0] - y[1] + Math.exp(2 * t); + output[2] = 5 * t**4 - 3 * t**2 + 2 * t; + }, + tolerance: 0.00000001, + solutionColNames: ['x', 'y', 'z'], + }; + + /** Non-stiff 3D test problem: exact solution */ + const exactNonStiff3D = (t: number) => { + const e1 = Math.exp(t); + const e2 = Math.exp(2 * t); + const e3 = Math.exp(3 * t); + const c = Math.cos(t); + const s = Math.sin(t); + + return new Float64Array([ + e1 + e3 + 0.1 * (3 * c - s) - 2 * e2, + -2 * e1 - e3 - 0.2 * (2 * s + 4 * c) + 3 * e2, + t**5 - t**3 + t**2, + ]); + }; + + /** Stiff 1D test problem: equation (see [1], p. 767) */ + const stiff1D = { + name: 'stiff 1D', + arg: {name: 't', start: 0, finish: 4, step: 0.01}, + initial: [0], + func: (t: number, y: Float64Array, output: Float64Array) => { + output[0] = -1000 * y[0] + 3000 - 2000 * Math.exp(-t); + }, + tolerance: 0.0000005, + solutionColNames: ['y'], + }; + + /** Stiff 1D test problem: exact solution (see [1], p. 767) */ + const exactStiff1D = (t: number) => { + return new Float64Array([ + 3 - 0.998 * Math.exp(-1000 * t) - 2.002 * Math.exp(-t), + ]); + }; + + /** Stiff 2D test problem: equations (see [1], p. 770) */ + const stiff2D = { + name: 'stiff 2D', + arg: {name: 't', start: 0, finish: 4, step: 0.01}, + initial: [52.29, 83.82], + func: (t: number, y: Float64Array, output: Float64Array) => { + output[0] = -5 * y[0] + 3 * y[1]; + output[1] = 100 * y[0] - 301 * y[1]; + }, + tolerance: 0.0000005, + solutionColNames: ['x', 'y'], + }; + + /** Stiff 2D test problem: exact solution (see [1], p. 770) */ + const exactStiff2D = (t: number) => { + return new Float64Array([ + 52.96 * Math.exp(-3.9899 * t) - 0.67 * Math.exp(-302.0101 * t), + 17.83 * Math.exp(-3.9899 * t) + 65.99 * Math.exp(-302.0101 * t), + ]); + }; + + /** Stiff 3D test problem: equations */ + const stiff3D = { + name: 'stiff 3D', + arg: {name: 't', start: 0, finish: 4, step: 0.01}, + initial: [2, 1, 0], + func: (t: number, y: Float64Array, output: Float64Array) => { + output[0] = -100 * y[0] + 100 * y[1] + y[2]; + output[1] = y[2]; + output[2] = -y[1]; + }, + tolerance: 0.00000001, + solutionColNames: ['x', 'y', 'z'], + }; + + /** Stiff 3D test problem: exact solution */ +const exactStiff3D = (t: number) => { + return new Float64Array([ + Math.exp(-100 * t) + Math.cos(t), + Math.cos(t), + -Math.sin(t), + ]); + }; + + type CorrectnessProblem = { + odes: ODEs, + exact: (t: number) => Float64Array, + }; + + const correctnessProblems = [ + {odes: nonStiff1D, exact: exactNonStiff1D}, + {odes: nonStiff2D, exact: exactNonStiff2D}, + {odes: nonStiff3D, exact: exactNonStiff3D}, + {odes: stiff1D, exact: exactStiff1D}, + {odes: stiff2D, exact: exactStiff2D}, + {odes: stiff3D, exact: exactStiff3D}, + ]; + + /** Return numerical solution error: maximum absolute deviation between approximate & exact solutions */ + function getError(method: (odes: ODEs) => Float64Array[], corProb: CorrectnessProblem): number { + let error = 0; + + // Get numerical solution + const approxSolution = method(corProb.odes); + + const exact = corProb.exact; + + const arg = approxSolution[0]; + + const pointsCount = arg.length; + const funcsCount = approxSolution.length - 1; + + // Compute error + for (let i = 0; i < pointsCount; ++i) { + const exactSolution = exact(arg[i]); + + for (let j = 0; j < funcsCount; ++j) + error = Math.max(error, Math.abs(exactSolution[j] - approxSolution[j + 1][i])); + } + + return error; + } + + console.log('Performance:\n'); + + methods.forEach((method, name) => { + console.log(' ', name); + + performanceProblems.forEach((odes, name) => { + const start = Date.now(); + method(odes); + const finish = Date.now(); + console.log(` ${name}: ${finish - start} ms.`); + }); + + console.log(); +}); + +console.log('Correctness:\n'); + +methods.forEach((method, name) => { + console.log(' ', name); + + correctnessProblems.forEach((problem) => { + const error = getError(method, problem) + console.log(` ${problem.odes.name}: ${error}`); + }); + + console.log(); +}); \ No newline at end of file diff --git a/libraries/diff-studio-utils/src/solver.ts b/libraries/diff-studio-utils/src/solver.ts new file mode 100644 index 0000000000..a6e756609c --- /dev/null +++ b/libraries/diff-studio-utils/src/solver.ts @@ -0,0 +1,39 @@ +// Solver of initial value problem + +import {ODEs, SolverOptions} from './solver-tools/solver-defs'; +import {mrt} from './solver-tools/mrt-method'; +import {ros3prw} from './solver-tools/ros3prw-method'; +import {ros34prw} from './solver-tools/ros34prw-method'; +import {getCallback} from './solver-tools/callbacks/callback-tools'; +import {METHOD} from './ui-constants'; + +/** Default solver of initial value problem. */ +export const solveDefault = (odes: ODEs): DG.DataFrame => ros34prw(odes); + +/** Return method specified by options */ +const getMethod = (options?: Partial) => { + if (options === undefined) + return ros34prw; + + switch (options.method) { + case METHOD.MRT: + return mrt; + + case METHOD.ROS3PRw: + return ros3prw; + + case METHOD.ROS34PRw: + return ros34prw; + + default: + return ros34prw; + } +}; + +/** Customizable solver of initial value problem. */ +export function solveIVP(odes: ODEs, options?: Partial): DG.DataFrame { + const callback = getCallback(options); + const method = getMethod(options); + + return method(odes, callback); +} diff --git a/libraries/diff-studio-utils/src/templates.ts b/libraries/diff-studio-utils/src/templates.ts new file mode 100644 index 0000000000..5890caf3c0 --- /dev/null +++ b/libraries/diff-studio-utils/src/templates.ts @@ -0,0 +1,159 @@ +import {CONTROL_EXPR} from './constants'; + +/** Basic template illustrating the simplest features */ +const TEMPLATE_BASIC = `${CONTROL_EXPR.NAME}: Template +${CONTROL_EXPR.DIF_EQ}: + dy/dt = -y + sin(t) / t + +${CONTROL_EXPR.ARG}: t + initial = 0.01 {min: 0.01; max: 10} + final = 15 {min: 15; max: 150} + step = 0.01 {min: 0.001; max: 0.1} + +${CONTROL_EXPR.INITS}: + y = 0 {min: 0; max: 9}`; + +/** Advanced template illustrating extanded features */ +const TEMPLATE_ADVANCED = `${CONTROL_EXPR.NAME}: Advanced +${CONTROL_EXPR.COMMENT}: + This is an advanced template. Modify it. Use multi-line formulas if needed. + Add new equations, expressions, constants & parameters. Edit these comment lines if required. +${CONTROL_EXPR.DIF_EQ}: + dx/dt = E1 * y + sin(t) + + dy/dt = E2 * x - pow(t, 5) + +${CONTROL_EXPR.EXPR}: + E1 = C1 * exp(-t) + P1 + E2 = C2 * cos(2 * t) + P2 + +${CONTROL_EXPR.ARG}: t + start = 0 + finish = 10 {min: 10; max: 100} + step = 0.01 + +${CONTROL_EXPR.INITS}: + x = 2 + y = 0 + +${CONTROL_EXPR.CONSTS}: + C1 = 1 + C2 = 3 + +${CONTROL_EXPR.PARAMS}: + P1 = 1 + P2 = -1 + +${CONTROL_EXPR.TOL}: 5e-5`; + +/** Extended template illustrating extanded features */ +const TEMPLATE_EXTENDED = `${CONTROL_EXPR.NAME}: Extended +${CONTROL_EXPR.TAGS}: model +${CONTROL_EXPR.DESCR}: 2D ordinary differential equations system sample +${CONTROL_EXPR.COMMENT}: + This is an extended template. It has additional scripting annotations. + +${CONTROL_EXPR.DIF_EQ}: + dx/dt = E1 * y + sin(t) + dy/dt = E2 * x - pow(t, 5) + +${CONTROL_EXPR.EXPR}: + E1 = C1 * exp(-t) + P1 + E2 = C2 * cos(2 * t) + P2 + +${CONTROL_EXPR.CONSTS}: + C1 = 1 + C2 = 3 + +${CONTROL_EXPR.PARAMS}: + P1 = 1 {category: Parameters; min: 1; max: 10} [P1 parameter] + P2 = -1 {category: Parameters; min: -10; max: -1} [P2 parameter] + +${CONTROL_EXPR.INITS}: + x = 2 {category: Initial values; min: 0; max: 5} [Initial value of x] + y = 0 {category: Initial values; min: -2; max: 2} [Initial value of y] + +${CONTROL_EXPR.ARG}: t + start = 0 {caption: Initial; category: Time; min: 0; max: 10} [Initial time of simulation] + finish = 10 {caption: Final; category: Time; min: 10; max: 20} [Final time of simulation] + step = 0.01 {caption: Step; category: Time; min: 0.01; max: 0.1; step: 0.001} [Time step of simlulation] + +${CONTROL_EXPR.TOL}: 5e-5`; + +/** Initial value problem templates */ +export enum TEMPLATES { + EMPTY = '', + BASIC = TEMPLATE_BASIC, + ADVANCED = TEMPLATE_ADVANCED, + EXTENDED = TEMPLATE_EXTENDED, +}; + +/** Template for the Demo app */ +export const DEMO_TEMPLATE = `${CONTROL_EXPR.NAME}: My model +${CONTROL_EXPR.DIF_EQ}: + dx/dt = y * cos(p * t) + sin(t) + dy/dt = x * sin(q * t) - exp(-t) + +${CONTROL_EXPR.ARG}: t + initial = 0 {min: 0; max: 3} + final = 8 {min: 4; max: 20} + step = 0.01 {min: 0.01; max: 0.1} + +${CONTROL_EXPR.INITS}: + x = 0 {min: 0; max: 5} + y = 1 {min: 1; max: 7} + +${CONTROL_EXPR.PARAMS}: + p = -2 {min: -2; max: 2} + q = -1 {min: -1; max: 1}`; + +/** Model for testing JS code use & expressions in the output */ +export const ENERGY_N_CONTROL = `${CONTROL_EXPR.NAME}: Energy-n-Control +${CONTROL_EXPR.TAGS}#tags: model +${CONTROL_EXPR.DESCR}#description: The energy-n-control demo model +${CONTROL_EXPR.COMMENT}#comment: + This is an artificial example. It illustrates the use of the Math module + and creating custom functions using JavaScript. + +${CONTROL_EXPR.DIF_EQ}#equations: + dx/dt = E1 * y * weight + dy/dt = E2 * x * weight + +${CONTROL_EXPR.EXPR}#expressions: + E1 = P2 * exp(-t) + E2 = P2 * cos(t) + weight = sin(PI * P1) + + // this makes further code shorter + ceil = Math.ceil + + // simple if-then-else custom function + func1 = (p, t) => (p < 8) ? ceil(t) : ceil(t**2 / 20) + + // another custom function using the previouse one + func2 = (p, t) => (p < 5) ? ceil(exp(-t) * 10) : func1(p, t); + + // usage of custom function + control = (P1 < 3) ? ceil(cos(10 * t)) : func2(P1, t) + + energy = x**2 + y**2 + +${CONTROL_EXPR.PARAMS}#parameters: + P1 = 1.1 {caption: frequency; category: Parameters; min: 1; max: 10} [The control parameter] + P2 = -1.1 {caption: shift; category: Parameters; min: -10; max: -1} [The shift parameter] + +${CONTROL_EXPR.INITS}#inits: + x = 2 {category: Initial state; min: 0; max: 5} [Initial value of the x coordinate] + y = 0 {category: Initial state; min: -2; max: 2} [Initial value of the y coordinate] + +${CONTROL_EXPR.OUTPUT}#output: + t {caption: time} + energy + control + +${CONTROL_EXPR.ARG}#argument: t + start = 0 {caption: Initial; category: Time; min: 0; max: 10} [Initial time of simulation] + finish = 10 {caption: Final; category: Time; min: 10; max: 20} [Final time of simulation] + step = 0.01 {caption: Step; category: Time; min: 0.01; max: 0.1; step: 0.001} [Time step of simlulation] + +${CONTROL_EXPR.TOL}#tolerance: 5e-5`; diff --git a/libraries/diff-studio-utils/src/ui-constants.ts b/libraries/diff-studio-utils/src/ui-constants.ts new file mode 100644 index 0000000000..e45c43052c --- /dev/null +++ b/libraries/diff-studio-utils/src/ui-constants.ts @@ -0,0 +1,326 @@ +/** Diff Studio application UI constants */ + +/** Hot keys */ +export enum HOT_KEY { + RUN = 'F5', +}; + +const COMPUTATION_TIME_UNITS = 'sec'; + +/** Tooltips messages */ +export enum HINT { + HELP = 'Open help in a new tab', + OPEN = 'Open model', + SAVE_MODEL = 'Save model', + SAVE_LOC = 'Save model to local file', + SAVE_MY = 'Save model to My files', + LOAD = 'Load model from local file', + BASIC = 'Basic template', + ADV = 'Advanced template', + EXT = 'Extended template', + CHEM = 'Mass-action kinetics illustration', + ROB = `Robertson's chemical reaction model`, + FERM = 'Fermentation process simulation', + PK = 'Pharmacokinetic model', + PKPD = 'Pharmacokinetic-pharmacodynamic model', + ACID = 'Gluconic acid production model', + NIM = 'Nimotuzumab disposition model', + BIO = 'Bioreactor simulation', + POLL = 'The chemical reaction part of the air pollution model', + CLEAR = 'Clear model', + TO_JS = 'Open in script editor', + APP = 'Export model to platform application with user interface', + OPEN_DS = 'Open model in Diff Studio', + SAVE = 'Save changes', + SENS_AN = 'Run sensitivity analysis', + FITTING = 'Run fitting inputs', + CONTINUE = 'Continue computations', + ABORT = 'Abort computations', + MAX_TIME = `Max computation time, ${COMPUTATION_TIME_UNITS}.`, + CLICK_RUN = `Click to run`, + SOLVE = `Solve equations (${HOT_KEY.RUN})`, + NO_MODELS = 'No models found', + EDIT = 'Edit', +}; // HINT + +/** UI titles */ +export enum TITLE { + LOAD = 'Load...', + IMPORT = 'Import...', + SAVE_TO_MY_FILES = 'Save to My files', + LOCAL_FILE = 'Local file...', + MY_FILES = 'My files...', + TO_MY_FILES = 'Save to My Files...', + AS_LOCAL = 'Save as Local File...', + TEMPL = 'Templates', + BASIC = 'Basic', + ADV = 'Advanced', + EXT = 'Extended', + LIBRARY = 'Library', + RECENT = 'Recent', + CHEM = 'Chem reactions', + ROB = `Robertson's model`, + FERM = 'Fermentation', + PK = 'PK', + PKPD = 'PK-PD', + ACID = 'Acid production', + NIM = 'Nimotuzumab', + BIO = 'Bioreactor', + POLL = 'Pollution', + CLEAR = 'Clear', + TO_JS = 'js', + MISC = 'Misc', + VARY = 'Vary inputs', + EDIT = 'Edit', + SOLVE = 'Solve', + FIT = 'Fit', + SOLUTION = 'Solution', + OPEN = 'Open', + BROWSE = 'Browse', + SAVE = 'Save', + APPS = 'Apps', + COMP = 'Compute', + DIF_ST = 'Diff Studio', + NAME = 'Name', + TYPE = 'Type', + INFO = 'Info', + IS_CUST = 'Custom model', + MY_MODELS = 'My Models', + NO_MODELS = 'None', +}; // TITLE + +/** Titles of template models */ +export const TEMPLATE_TITLES = [TITLE.BASIC, TITLE.ADV, TITLE.EXT]; + +/** Titles of example models */ +export const EXAMPLE_TITLES = [TITLE.CHEM, TITLE.ROB, TITLE.FERM, TITLE.PK, + TITLE.PKPD, TITLE.ACID, TITLE.NIM, TITLE.BIO, TITLE.POLL]; + +/** Models' tooltips map */ +export const MODEL_HINT = new Map([ + [TITLE.BASIC, HINT.BASIC], + [TITLE.ADV, HINT.ADV], + [TITLE.EXT, HINT.EXT], + [TITLE.CHEM, HINT.CHEM], + [TITLE.ROB, HINT.ROB], + [TITLE.FERM, HINT.FERM], + [TITLE.PK, HINT.PK], + [TITLE.PKPD, HINT.PKPD], + [TITLE.ACID, HINT.ACID], + [TITLE.NIM, HINT.NIM], + [TITLE.BIO, HINT.BIO], + [TITLE.POLL, HINT.POLL], +]); + +/** Help links */ +export enum LINK { + DIF_STUDIO_REL = '/help/compute/diff-studio', + MODELS = '/help/compute/models', + DIF_STUDIO = 'https://datagrok.ai/help/compute/diff-studio', + SENS_AN = '/help/compute/function-analysis#sensitivity-analysis', + FITTING = '/help/compute/function-analysis#parameter-optimization', + CHEM_REACT = `${MODELS}#chem-reactions`, + FERMENTATION = `${MODELS}#fermentation`, + GA_PRODUCTION = `${MODELS}#acid-production`, + NIMOTUZUMAB = `${MODELS}#nimotuzumab`, + PK = `${MODELS}#pk`, + PKPD = `${MODELS}#pk-pd`, + ROBERTSON = `${MODELS}#robertson-model`, + BIOREACTOR = `${MODELS}#bioreactor`, + POLLUTION = `${MODELS}#pollution`, + COMPUTE = 'https://datagrok.ai/help/compute', + LOAD_SAVE = `${DIF_STUDIO_REL}#working-with-models`, + MODEL_COMPONENTS = `${DIF_STUDIO_REL}#model-components-and-syntax`, + INTERFACE = `${DIF_STUDIO_REL}#user-interface-options`, +}; + +/** Error messages */ +export enum ERROR_MSG { + SOLVING_FAILS = 'Solving fails', + APP_CREATING_FAILS = 'Application creating fails', + EXPORT_TO_SCRIPT_FAILS = 'Export to JavaScript script fails', + SCRIPTING_ISSUE = 'Platform scripting issue', + UI_ISSUE = 'UI creating issue', + MISSING_CLOSING_BRACKET = 'Annotation error: **"]"** is missing', + INCORRECT_BRACES_USE = 'Annotation error: incorrect use of **"{}"**', + MISSING_COLON = 'Annotation error: **":"** is missing', + CHECK_ARGUMENTS = ' (check the **#argument** block)', + INCORRECT_ARG_SEGM = 'Incorrect limits for the argument', + OPEN_IN_DIF_STUD = 'To change the model, open it in Diff Studio.', + FAILED_TO_SAVE = 'Failed to save model to the file', + INCORRECT_MODEL = 'Incorrect model', + SENS_AN_FAILS = 'Sensitivity Analysis fails', + FITTING_FAILS = 'Fitting fails', + PLATFORM_ISSUE = 'Platform issue', + INCORTECT_INPUT = 'Incorrect input caption', + NANS_OBTAINED = 'Failed computations: obtained NaN-s. Check inputs and formulas.', +}; + +/** Lookup table fails */ +export enum LOOKUP_DF_FAIL { + LOAD = 'Failed to load lookup table: ', + PLATFORM = 'the platform issue', + FUNCTION = 'incorrect function', + NO_DF = 'no dataframe', + INCORRECT = 'incorrect dataframe, ', + ROWS = 'at least one row is needed.', + NULLS = 'missing values are not allowed.', + NUMS = 'no numerical columns.', + CHOICES = 'first column must contain strings.', +}; + +/** Lookup table specification fails */ +export enum LOOKUP_EXPR_FAIL { + MISSING = 'Incorrect input: missing ', +}; + +/** Other UI constants */ +export enum MISC { + VIEW_DEFAULT_NAME = 'Template', + MODEL_FILE_EXT = 'ivp', + FILE_DEFAULT_NAME = `equations.${MODEL_FILE_EXT}`, + DEFAULT = 'Default', + CHOICES = 'choices', + NAME = 'name', + CATEGORY = 'category', + CAPTION = 'caption', + TOOLTIP = 'tooltip', + IS_NOT_DEF = 'is not defined', + UNEXPECTED = 'Unexpected identifier', + PROP_OF_NULL = 'Cannot set properties of null', +}; + +/** Warning dialog lines */ +export enum WARNING { + TITLE = 'WARNING', + CHECK = 'Show this warning', + OVERWRITE_MODEL = 'Overwrite the current model?', + OVERWRITE_FILE = 'Overwrite existing file?', + CONTINUE = 'Continue?', + CHECK_PERF = 'Check time', + TIME_LIM = 'Time limit', + UNITS = COMPUTATION_TIME_UNITS, + PREVIEW = `Model preview is unavailble for this type. Use "${MISC.MODEL_FILE_EXT}" instead`, +}; + +/** Code completion infos */ +export enum INFO { + NAME = 'name of the model', + TAGS = 'scripting tags', + DESCR = 'descritpion of the model', + DIF_EQ = 'block of differential equation(s) specification', + EXPR = 'block of auxiliary expressions & computations', + ARG = 'independent variable specification', + INITS = 'initial values of the model', + PARAMS = 'parameters of the model', + CONSTS = 'constants definition', + TOL = 'tolerance of numerical solution', + EPS_SCALE = 'scale used when computing Jacobian', + LOOP = 'loop feature', + UPDATE = 'update model feature', + OUTPUT = 'output specification', + COMMENT = 'block with comments', + SOLVER = 'solver options', + INPUS = 'path to table with inputs', +}; + +/** Demo app help info */ +export const demoInfo = `# Diff Studio +In-browser solver of ordinary differential equations +([ODEs](https://en.wikipedia.org/wiki/Ordinary_differential_equation)) + +# Interactivity +Play with the model inputs. Move sliders to explore its behavior. + +# Model +Turn on the **${TITLE.EDIT}** toggle on the top panel to modify formulas. + +# No-code +Define equations in a declarative form. + +# Examples +Click icon and +explore **Library**. + +# Scripting +Click **** icon to export model to JavaScript script. + +# Analysis +Turn off the **${TITLE.EDIT}** toggle, and perform analysis: +* Click the **Fit** icon on the top panel to [optimize inputs](${LINK.FITTING}). +* Click the **Sensitivity** icon to run [sensitivity analysis](${LINK.SENS_AN}). + +# Learn more +* [Diff Studio](${LINK.DIF_STUDIO}) +* [Compute](${LINK.COMPUTE})`; + +/** Inputs types */ +export enum INPUT_TYPE { + FLOAT = 'Float', + INT = 'Int', +}; + +/** Path related consts */ +export enum PATH { + APP_DATA_DS = '/files/system.appdata/diffstudio', + APPS_DS = '/apps/DiffStudio', + MODEL = `?model=`, + CUSTOM = `${MODEL}custom`, + EMPTY = `${MODEL}empty`, + EQ = '=', + AND = '&', + PARAM = `?params:`, + BROWSE = 'browse', + RECENT = 'diff-studio-recent.d42', + MY_FILES = 'Myfiles', + HOME = 'Home', + SYSTEM = 'System', +}; + +/** UI time constants */ +export enum UI_TIME { + DOCK_EDITOR_TIMEOUT = 100, + PREVIEW_RUN_SOLVING = 105, + APP_RUN_SOLVING = 100, + SOLV_DEFAULT_TIME_SEC = 5, + SOLV_TIME_MIN_SEC = 1, + BROWSING = APP_RUN_SOLVING + 500, + SWITCH_TO_FOLDER = 1, + WGT_CLICK = 10, +}; + +/** Numerical methods names */ +export enum METHOD { + MRT = 'mrt', + ROS3PRw = 'ros3prw', + ROS34PRw = 'ros34prw', +}; + +export const DOCK_RATIO = 0.3; + +export const MAX_RECENT_COUNT = 10; + +export const CUSTOM_MODEL_IMAGE_LINK = 'images/custom.png'; + +/** Model image link */ +export const modelImageLink = new Map([ + [TITLE.BASIC, 'images/basic.png'], + [TITLE.ADV, 'images/advanced.png'], + [TITLE.EXT, 'images/extended.png'], + [TITLE.CHEM, 'images/chem-react.png'], + [TITLE.ROB, 'images/robertson.png'], + [TITLE.FERM, 'images/fermentation.png'], + [TITLE.PK, 'images/pk.png'], + [TITLE.PKPD, 'images/pk-pd.png'], + [TITLE.ACID, 'images/ga-production.png'], + [TITLE.NIM, 'images/nimotuzumab.png'], + [TITLE.BIO, 'images/bioreactor.png'], + [TITLE.POLL, 'images/pollution.png'], +]); + +/** Inputs table constants */ +export enum INPUTS_DF { + MIN_ROWS_COUNT = 1, + INP_NAMES_IDX = 0, + INPUT_SETS_COL_IDX = 0, +}; diff --git a/libraries/diff-studio-utils/src/use-cases.ts b/libraries/diff-studio-utils/src/use-cases.ts new file mode 100644 index 0000000000..589b00e1e0 --- /dev/null +++ b/libraries/diff-studio-utils/src/use-cases.ts @@ -0,0 +1,608 @@ +/* eslint-disable max-len */ +import {CONTROL_EXPR} from './constants'; + +/** Chemical reactions, mass-action kinetics */ +const CHEM_REACT_MODEL = `${CONTROL_EXPR.NAME}: Chem react +${CONTROL_EXPR.TAGS}: model +${CONTROL_EXPR.DESCR}: Mass-action kinetics illustration +${CONTROL_EXPR.COMMENT}: + Source: https://doi.org/10.1002/ijch.201800003. +${CONTROL_EXPR.DIF_EQ}: + dx1/dt = -k1 * x1 + k2 * (x2)**2 + k3 * x1 * x3 + - k4 * (x1)**2 - 2 * k5 * (x1)**2 + k6 * x2 * x4 + + dx2/dt = 2 * k1 * x1 - 2 * k2 * (x2)**2 + + k5 * (x1)**2 - k6 * x2 * x4 + + dx3/dt = -k3 * x1 * x3 + k4 * (x1)**2 + k6 * x2 * x4 + + dx4/dt = k5 * (x1)**2 - k6 * x2 * x4 + +${CONTROL_EXPR.PARAMS}: + k1 = 0.7 {category: Reaction parameters; min: 0.1; max: 5} + k2 = 0.9 {category: Reaction parameters; min: 0.1; max: 5} + k3 = 1.2 {category: Reaction parameters; min: 0.1; max: 5} + k4 = 3.4 {category: Reaction parameters; min: 0.1; max: 5} + k5 = 2.3 {category: Reaction parameters; min: 0.1; max: 5} + k6 = 4.5 {category: Reaction parameters; min: 0.1; max: 5} + +${CONTROL_EXPR.INITS}: + x1 = 1 {units: mol/L; category: Initial concentrations; min: 0; max: 2} + x2 = 0 {units: mol/L; category: Initial concentrations; min: 0; max: 2} + x3 = 0 {units: mol/L; category: Initial concentrations; min: 0; max: 2} + x4 = 0 {units: mol/L; category: Initial concentrations; min: 0; max: 2} + +${CONTROL_EXPR.ARG}: t + initial = 0 {units: min; caption: Initial; category: Time; min: 0; max: 5} [Initial time of simulation] + final = 6 {units: min; caption: Final; category: Time; min: 6; max: 10} [Final time of simulation] + step = 0.01 {units: min; caption: Step; category: Time; min: 0.01; max: 0.1; step: 0.001} [Time step of simlulation] + +${CONTROL_EXPR.TOL}: 5e-5`; + +/** Robertson's chemical reaction model - stiff ODEs */ +const ROBERTSON_MODEL = `${CONTROL_EXPR.NAME}: Robertson +${CONTROL_EXPR.TAGS}: model +${CONTROL_EXPR.DESCR}: Robertson chemical reaction model +${CONTROL_EXPR.COMMENT}: This is classic example of stiff ODEs. +${CONTROL_EXPR.DIF_EQ}: + dA/dt = -0.04 * A + 1e4 * B * C + dB/dt = 0.04 * A - 1e4 * B * C - 3e7 * B**2 + dC/dt = 3e7 * B**2 + +${CONTROL_EXPR.INITS}: + A = 1 {units: mol/L; category: Initial concentrations; min: 0; max: 5} + B = 0 {units: mol/L; category: Initial concentrations; min: 0; max: 5} + C = 0 {units: mol/L; category: Initial concentrations; min: 0; max: 5} + +${CONTROL_EXPR.ARG}: t + start = 0 {units: sec; caption: Initial; category: Time; min: 0; max: 1} [Initial time of simulation] + finish = 40 {units: sec; caption: Final; category: Time; min: 2; max: 50} [Final time of simulation] + step = 0.01 {units: sec; caption: Step; category: Time; min: 0.01; max: 0.1} [Time step of simlulation] + +${CONTROL_EXPR.TOL}: 1e-7`; + +/** Fermentation process simulation */ +const FERMENTATION_MODEL = `${CONTROL_EXPR.NAME}: Fermentation +${CONTROL_EXPR.TAGS}: model +${CONTROL_EXPR.DESCR}: Simulation of fermentation process in the ethanol production +${CONTROL_EXPR.COMMENT}: + Source: https://core.ac.uk/download/pdf/11737483.pdf. +${CONTROL_EXPR.DIF_EQ}: + dP/dt = r * X + dS/dt = -q * X + dX/dt = V * S / (K + S) * X + +${CONTROL_EXPR.ARG}: t + initial = 0 {units: d; caption: Initial; category: Time; min: 0; max: 10} [Initial time of simulation] + final = 70 {units: d; caption: Final; category: Time; min: 15; max: 100} [Final time of simulation] + step = 0.01 {units: d; caption: Step; category: Time; min: 0.01; max: 1} [Time step of simlulation] + +${CONTROL_EXPR.INITS}: + P = 4.276 {units: ml/ml; caption: ethanol, P; category: Initials; min: 3; max: 6} [Concentration of ethanol] + S = 0.3185 {units: mg/ml; caption: glucose, S; category: Initials; min: 0.1; max: 0.5} [Concentration of glucose] + X = 9.2e-2 {units: mg/ml; caption: saccharomyces, X; category: Initials; min: 0.01; max: 0.2} [Saccharomyces wet weight] + +${CONTROL_EXPR.PARAMS}: + r = 1.3455 {units: mg/ml; category: Parameters; min: 1; max: 2} [The growth rate of ethanol] + q = 1.1129e-2 {units: mg/ml; category: Parameters; min: 0.01; max: 0.2} [The rate of glucose consumption] + V = 8.7e-2 {units: mg/ml; category: Parameters; min: 0.01; max: 0.2} [The maximal growth rate of Saccharomyces] + K = 6.28e-2 {category: Parameters} [Michaelis-Menten constant] + +${CONTROL_EXPR.TOL}: 1e-7`; + +/** PK simulation */ +const PK_MODEL = `${CONTROL_EXPR.NAME}: PK +${CONTROL_EXPR.TAGS}: model +${CONTROL_EXPR.DESCR}: Pharmacokinetic (PK) simulation: one-compartment model +${CONTROL_EXPR.DIF_EQ}: + d(depot)/dt = -KA * depot + d(centr)/dt = KA * depot - CL * C2 + +${CONTROL_EXPR.EXPR}: + C2 = centr / V2 + +${CONTROL_EXPR.LOOP}: + count = 1 {category: Dosing; min: 1; max: 10} [Number of doses] + depot += dose + +${CONTROL_EXPR.ARG}: t + start = 0 {units: h; caption: begin; category: Dosing; min: 0; max: 1} [Begin of dosing interval] + final = 12 {units: h; caption: end; category: Dosing; min: 5; max: 15} [End of dosing interval] + step = 0.01 {units: h; caption: step; category: Dosing; min: 0.01; max: 0.1} [Time step of simlulation] + +${CONTROL_EXPR.INITS}: + depot = 0 {category: Initial values} + centr = 0 {category: Initial values} [Central] + +${CONTROL_EXPR.PARAMS}: + dose = 1e4 {category: Dosing; min: 1e3; max: 2e4; step: 1e3} [Dosage] + KA = 0.3 {caption: rate constant; category: Parameters; min: 0.1; max: 1} + CL = 2 {caption: clearance; category: Parameters; min: 1; max: 5} + V2 = 4 {caption: central volume; category: Parameters; min: 1; max: 10} [Central compartment volume] + +${CONTROL_EXPR.SOLVER}: {method: 'mrt'; maxTimeMs: 50} + +${CONTROL_EXPR.TOL}: 1e-9`; + +/** PK-PD simulation */ +const PK_PD_MODEL = `${CONTROL_EXPR.NAME}: PK-PD +${CONTROL_EXPR.TAGS}: model +${CONTROL_EXPR.DESCR}: Pharmacokinetic-pharmacodynamic (PK-PD) simulation: two-compartment model +${CONTROL_EXPR.DIF_EQ}: + d(depot)/dt = -KA * depot + d(centr)/dt = KA * depot - CL * C2 - Q * C2 + Q * C3 + d(peri)/dt = Q * C2 - Q * C3 + d(eff)/dt = Kin - Kout * (1 - C2/(EC50 + C2)) * eff + +${CONTROL_EXPR.EXPR}: + C2 = centr / V2 + C3 = peri / V3 + +${CONTROL_EXPR.LOOP}: + count = 10 {caption: count; category: Dosing; min: 1; max: 20} [Number of doses] + depot += dose + +${CONTROL_EXPR.ARG}: t + start = 0 {units: h; caption: begin; category: Dosing; min: 0; max: 1} [Begin of dosing interval] + final = 12 {units: h; caption: end; category: Dosing; min: 5; max: 15} [End of dosing interval] + step = 0.1 {units: h; caption: step; category: Dosing; min: 0.01; max: 0.1} [Time step of simlulation] + +${CONTROL_EXPR.INITS}: + depot = 0 {category: Initial values} + centr = 0 {category: Initial values} [Central] + peri = 0 {category: Initial values} [Peripheral] + eff = 0.2 {category: Initial values} [Effective compartment rate] + +${CONTROL_EXPR.PARAMS}: + dose = 1e4 {category: Dosing; min: 1e3; max: 2e4; step: 1e3} [Dosage] + KA = 0.3 {caption: rate constant; category: Parameters; min: 0.1; max: 1} + CL = 2 {caption: clearance; category: Parameters; min: 1; max: 5} + V2 = 4 {caption: central volume; category: Parameters; min: 1; max: 10} [Central compartment volume] + Q = 1 {caption: inter rate; category: Parameters; min: 0.1; max: 1} [Intercompartmental rate] + V3 = 30 {caption: peri volume; category: Parameters; min: 20; max: 40} [Peripheral compartment volume] + EC50 = 8 {caption: effect; category: Parameters; min: 1; max: 10} + Kin = 0.2 {caption: Kin; category: Parameters; min: 0.1; max: 0.5} [The first-order production constant] + Kout = 0.2 {caption: Kout; category: Parameters; min: 0.1; max: 0.5} [The first-order dissipation rate constant] + +${CONTROL_EXPR.TOL}: 1e-9`; + +/** Gluconic acid production */ +const ACID_PROD_MODEL = `${CONTROL_EXPR.NAME}: GA-production +${CONTROL_EXPR.TAGS}: model +${CONTROL_EXPR.DESCR}: Gluconic acid (GA) production by Aspergillus niger modeling +${CONTROL_EXPR.DIF_EQ}: + dX/dt = rX + dS/dt = -gamma * rX - lambda * X + dO/dt = Kla * (Cod - O) - delta * rX - phi * X + dP/dt = alpha * rX + beta * X + +${CONTROL_EXPR.EXPR}: + mu = muM * S / (Ks + S) * O / (Ko + O) + rX = mu * X + +${CONTROL_EXPR.ARG}: t, 1-st stage + _t0 = 0 {units: h; caption: initial; category: Misc} [Start of the process] + _t1 = 60 {units: h; caption: 1-st stage; category: Durations; min: 20; max: 80} [Duration of the 1-st stage] + step = 0.1 {units: h; caption: step; category: Misc; min: 0.01; max: 1} [Time step of simlulation] + +${CONTROL_EXPR.UPDATE}: 2-nd stage + duration = overall - _t1 + S += 70 + +${CONTROL_EXPR.INITS}: + X = 5 {units: kg/m³; caption: biomass; category: Initial concentrations; min: 1; max: 10} [Aspergillus niger biomass] + S = 150 {units: kg/m³; caption: glucose; category: Initial concentrations; min: 50; max: 200} [Glucose] + O = 7 {units: kg/m³; caption: oxygen; category: Initial concentrations; min: 1; max: 10} [Dissolved oxygen] + P = 0 {units: kg/m³; caption: acid; category: Initial concentrations; min: 0; max: 0.1} [Gluconic acid] + +${CONTROL_EXPR.OUTPUT}: + t {caption: time} + X {caption: biomass} + S {caption: glucose} + O {caption: oxygen} + P {caption: acid} + +${CONTROL_EXPR.PARAMS}: + overall = 100 {units: h; category: Durations; min: 100; max: 140} [Overall duration] + muM = 0.668 {units: 1/h; category: Parameters} [Monod type model parameter] + alpha = 2.92 {category: Parameters} [Monod type model parameter] + beta = 0.131 {units: 1/h; category: Parameters} [Monod type model parameter] + gamma = 2.12 {category: Parameters} [Monod type model parameter] + lambda = 0.232 {units: 1/h; category: Parameters} [Monod type model parameter] + delta = 0.278 {category: Parameters} [Monod type model parameter] + phi = 4.87e-3 {units: 1/h; category: Parameters} [Monod type model parameter] + Ks = 1.309e2 {units: g/L; category: Parameters} [Monod type model parameter] + Ko = 3.63e-4 {units: g/L; category: Parameters} [Monod type model parameter] + Kla = 1.7e-2 {units: 1/s; category: Parameters} [Volumetric mass transfer coefficient] + Cod = 15 {units: kg/m³; category: Parameters} [Liquid phase dissolved oxygen saturation concentration] + +${CONTROL_EXPR.TOL}: 1e-9`; + +/** Nimotuzumab disposition model */ +const NIMOTUZUMAB_MODEL = `${CONTROL_EXPR.NAME}: Nimotuzumab +${CONTROL_EXPR.TAGS}: model +${CONTROL_EXPR.DESCR}: Nimotuzumab disposition model +${CONTROL_EXPR.COMMENT}: + Source: https://www.mdpi.com/1999-4923/12/12/1147 +${CONTROL_EXPR.DIF_EQ}: + dA1/dt = (-(CL * A3 / V1 + Q / V1) * A1 + Q / V2 * A2 - kint * Rtot * A1 / (Kss + A1 / V1)) + / (1 + Rtot * Kss / (Kss + A1 / V1)**2) + + dA2/dt = (Q / V1 * A1 - Q / V2 * A2) / (1 + Rtotp * Kss / (Kss + A2 / V2)**2) + + dA3/dt = Kin * (1 + Smax * A1**gamma / (S50**gamma + A1**gamma)) - Kout * A3 + +${CONTROL_EXPR.ARG}: t + start = 0 {caption: Initial; category: Time; min: 0; max: 10} [Initial time of simulation] + finish = 70 {caption: Final; category: Time; min: 50; max: 100} [Final time of simulation] + step = 0.01 {caption: Step; category: Time; min: 0.01; max: 1} [Time step of simlulation] + +${CONTROL_EXPR.INITS}: + A1 = 0.43 {category: Initials; min: 0.2; max: 0.6} [Central] + A2 = 0.55 {category: Initials; min: 0.3; max: 0.8} [Periferal] + A3 = 10.32 {category: Initials; min: 5; max: 15} [Mediator] + +${CONTROL_EXPR.OUTPUT}: + t {caption: Time, h} + A1 {caption: Central} + A2 {caption: Periferal} + +${CONTROL_EXPR.PARAMS}: + CL = 9.64e-3 {category: Parameters; min: 0.005; max: 0.015} [Non-specific clearance] + V1 = 2.63 {category: Parameters; min: 0.7; max: 3} [Apparent volume of distribution of the central compartment] + Q = 2.88e-2 {category: Parameters; min: 0.015; max: 0.035} [Distribution clearance of free nimotuzumab between the central and peripheral compartment] + V2 = 9.92e-3 {category: Parameters; min: 0.0067; max: 0.0265} [Apparent volume of distribution of the peripheral compartment] + Kss = 15.5 {category: Parameters; min: 10; max: 45} [Quasi steady state constant for the EGFR] + kint = 1e-2 {category: Parameters; min: 0.01; max: 0.05} [Internalization rate for nimotuzumab-EGFR] + Rtot = 1.05e-2 {category: Parameters; min: 0.007; max: 0.032} [EGFR in the central compartment] + Rtotp = 956 {category: Parameters; min: 120; max: 3500} [EGFR in the peripheral compartment] + Kin = 1.33e-2 {category: Parameters; min: 0.01; max: 0.05} [Zero-order synthesis rate constant] + Kout = 1.33e-2 {category: Parameters; min: 0.001; max: 0.02} [First-order elimination rate constant] + S50 = 8.57 {category: Parameters; min: 5; max: 12} [Concentration of free nimotuzumab in the central compartment that achieves the half of Smax] + Smax = 3.18 {category: Parameters; min: 1; max: 5} [Maximal effect of the stimulation] + gamma = 0.5 {category: Parameters; min: 0.1; max: 1} [Hill coefficient of the sigmoid function] + +${CONTROL_EXPR.TOL}: 1e-7`; + +/** Bioreactor simulation */ +const BIOREACTOR_MODEL = `${CONTROL_EXPR.NAME}: Bioreactor +${CONTROL_EXPR.TAGS}: model +${CONTROL_EXPR.DESCR}: Bioreactor simulation +${CONTROL_EXPR.COMMENT}: + Source: https://doi.org/10.1074/jbc.RA117.000303. +${CONTROL_EXPR.DIF_EQ}: + + d(FFox)/dt = -E11 + E12 + + d(KKox)/dt = -E21 + E22 + + d(FFred)/dt = E11 - E12 - E31 + E32 + + d(KKred)/dt = E21 - E22 - E41 + E42 + + d(Ffree)/dt = 2 * (E31 - E32) - E51 + E52 + + d(Kfree)/dt = 2 * (E41 - E42) - E51 + E52 + + d(FKred)/dt = E51 - E52 - E61 + E62 + + d(FKox)/dt = E61 - E62 + + d(MEAthiol)/dt = 2 * (-E11 + E12 - E21 + E22 + E31 + E41 - E32 - E42 - E62 - ktox * E71 * E72) + - (MEAthiol + MA) * (Fin + Fper) / VL + + d(CO2)/dt = (Fin * pO2sat * 0.0022 - 2 * Fper * CO2) / VL + OTR - 0.5 * ktox * E71 * E72 + + d(yO2P)/dt = -OTR * (VL / (VTV - VL)) * R * T * P + yO2in * qin - yO2P * qout + + d(CYST)/dt = ktox * E71 * E72 - krcyst * CYST * E0 - (Fin + Fper) * CYST / VL + + d(VL)/dt = Fin - Fper + +${CONTROL_EXPR.EXPR}: + + KF = pow(VL, -0.65) * 0.065 * pow(speed**3 * diam**5 * power / 2.16e12, 0.361) + + MA = MEAthiol * pow(10, pH - pKa2MEA) + + qout = qin - KF * (yO2P * H - CO2) * VL * R * T / (P * 1000) + + OTR = KF * (yO2P * H - CO2) + + Fin = t < switchTime ? 0 : 0.025 + + Fper = t < switchTime ? 0.025 : Fin + + E0 = VLinit / VL + + E1 = (MA * E0)**2 + + E11 = k1red * FFox * E0 * E1 + + E12 = k1ox * FFred * E0 + + E31 = k2Fd * FFred * E0 + + E32 = k2Fa * (Ffree * E0)**2 * E1 + + E21 = k1red * KKox * E0 * E1 + + E22 = k1ox * KKred * E0 + + E41 = k2Kd * KKred * E0 + + E42 = k2Ka * (Kfree * E0)**2 * E1 + + E51 = k3FKa * Ffree * E0 * Kfree * E0 + + E52 = k3FKd * FKred * E0 + + E61 = k4ox * FKred * E0 * (CYST * E0)**2 + + E62 = k4red * FKox * E0 * E1 + + E70 = E0 * CO2 + + E71 = (MEAthiol * E0)**2 + + E72 = (E70 >= 0) ? sqrt(E70) : 0 + +${CONTROL_EXPR.ARG}: t + t0 = 0.0 {units: min; caption: Initial; category: Time} [Initial time of simulation] + t1 = 1000 {units: min; caption: Final; category: Time; min: 500; max: 1000} [Final time of simulation] + h = 1 {units: min; caption: Step; category: Time; min: 0.1; max: 2} [Time step of simlulation] + +${CONTROL_EXPR.INITS}: + FFox = 0.2 {units: mmol/L; category: Initial values; min: 0.15; max: 0.25; step: 0.01} [FF oxidized] + KKox = 0.2 {units: mmol/L; category: Initial values; min: 0.15; max: 0.25; step: 0.01} [KK oxidized] + FFred = 0.1 {units: mmol/L; category: Initial values; min: 0.08; max: 0.12; step: 0.01} [FF reduced] + KKred = 0.1 {units: mmol/L; category: Initial values; min: 0.08; max: 0.12; step: 0.01} [KK reduced] + Ffree = 0 {units: mmol/L; category: Initial values} [F free] + Kfree = 0 {units: mmol/L; category: Initial values} [K free] + FKred = 0 {units: mmol/L; category: Initial values} [FK reduced] + FKox = 0 {units: mmol/L; category: Initial values} [FK oxidized] + MEAthiol = 15 {units: mmol/L; category: Initial values; min: 10; max: 16} [MEAthiol] + CO2 = 0.12 {units: mmol/L; category: Initial values; min: 0.09; max: 0.15} [Dissolved oxygen] + yO2P = 0.209 {units: atm; category: Initial values} [Atm headspace] + CYST = 0 {units: mmol/L; category: Initial values} [Cystamine] + VL = 7.2 {units: L; category: Initial values} [Liquid volume] + +${CONTROL_EXPR.OUTPUT}: + t + FFox {caption: FFox(t)} + KKox {caption: KKox(t)} + FFred {caption: FFred(t)} + KKred {caption: KKred(t)} + Ffree {caption: Ffree(t)} + Kfree {caption: Kfree(t)} + FKred {caption: FKred(t)} + FKox {caption: FKox(t)} + MEAthiol {caption: MEAthiol(t)} + CO2 {caption: CO2(t)} + yO2P {caption: yO2P(t)} + CYST {caption: CYST(t)} + VL {caption: VL(t)} + +${CONTROL_EXPR.CONSTS}: + VLinit = 7.2 + VTV = 10 + speed = 400 + diam = 6 + power = 2.1 + pH = 7.4 + k1red = 5.604e-2 + k1ox = 1.08e-2 + k2Fd = 1.35 + k2Fa = 1.104e8 + k2Kd = 4.038e-2 + k2Ka = 1.2e8 + k3FKa = 1.812e8 + k3FKd = 1.188e-2 + k4ox = 1.08e-2 + k4red = 5.604e-2 + ktox = 5e-3 + krcyst = 0 + pO2sat = 100 + pKa2MEA = 8.18 + H = 1.072069378 + R = 8.2e-2 + +${CONTROL_EXPR.PARAMS}: + qin = 1 {units: L/min; caption: Gas; category: Parameters; min: 0.5; max: 1.5} [Gas to headspace] + yO2in = 0.21 { caption: O2 fraction; category: Parameters; min: 0.1; max: 0.9} [Oxygen mole fraction] + T = 300 {units: K; caption: temperature; category: Parameters; min: 250; max: 350} [System temperature] + P = 1 {units: atm; caption: pressure; category: Parameters; min: 1; max: 2} [Headspace pressure] + switchTime = 135 {units: min; caption: switch at; category: Time; min: 70; max: 180; step: 10} [Switch mode time] + +${CONTROL_EXPR.INPUTS}: mode {caption: Process mode; category: Process parameters; choices: OpenFile("System:AppData/DiffStudio/library/bioreactor-inputs.csv")} [Reactions flow mode]`; + +/** Pollution model */ +const POLLUTION_MODEL = `${CONTROL_EXPR.NAME}: Pollution +${CONTROL_EXPR.TAGS}: model +${CONTROL_EXPR.DESCR}: The chemical reaction part of the air pollution model developed at The Dutch National Institute of Public Health and Environmental Protection +${CONTROL_EXPR.COMMENT}: + Source: https://archimede.uniba.it/~testset/report/pollu.pdf +${CONTROL_EXPR.DIF_EQ}: + dy1/dt = -(r1 + r10 + r14 + r23 + r24) + (r2 + r3 + r9 + r11 + r12 + r22 + r25) + + dy2/dt = -r2 - r3 - r9 - r12 + r1 + r21 + + dy3/dt = -r15 + r1 + r17 + r19 + r22 + + dy4/dt = -r2 - r16 - r17 - r23 + r15 + + dy5/dt = -r3 +2 * r4 + r6 +r7 +r13 + r20 + + dy6/dt = -r6 - r8 - r14 - r20 + r3 + 2 * r18 + + dy7/dt = -r4 - r5 - r6 + r13 + + dy8/dt = r4 + r5 + r6 + r7 + + dy9/dt = -r7 - r8 + + dy10/dt = -r12 +r7 + r9 + + dy11/dt = -r9 - r10 + r8 + r11 + + dy12/dt = r9 + + dy13/dt = -r11 + r10 + + dy14/dt = -r13 + r12 + + dy15/dt = r14 + + dy16/dt = -r18 - r19 + r16 + + dy17/dt = -r20 + + dy18/dt = r20 + + dy19/dt = -r21 - r22 - r24 + r23 + r25 + + dy20/dt = -r25 + r24 + +${CONTROL_EXPR.EXPR}: + r1 = k1 * y1 + + r2 = k2 * y2 * y4 + + r3 = k3 * y5 * y2 + + r4 = k4 * y7 + + r5 = k5 * y7 + + r6 = k6 * y7 * y6 + + r7 = k7 * y9 + + r8 = k8 * y9 * y6 + + r9 = k9 * y11 * y2 + + r10 = k10 * y11 * y1 + + r11 = k11 * y13 + + r12 = k12 * y10 * y2 + + r13 = k13 * y14 + + r14 = k14 * y1 * y6 + + r15 = k15 * y3 + + r16 = k16 * y4 + + r17 = k17 * y4 + + r18 = k18 * y16 + + r19 = k19 * y16 + + r20 = k20 * y17 * y6 + + r21 = k21 * y19 + + r22 = k22 * y19 + + r23 = k23 * y1 * y4 + + r24 = k24 * y19 * y1 + + r25 = k25 * y20 + +${CONTROL_EXPR.ARG}: t + t0 = 0 {units: min; caption: Initial; category: Time; min: 0; max: 0.9} [Initial time of simulation] + t1 = 60 {units: min; caption: Final; category: Time; min: 1; max: 100; step: 1} [Final time of simulation] + h = 0.01 {units: min; caption: Step; category: Time; min: 0.001; max: 0.1; step: 0.001} [Time step of simlulation] + + +${CONTROL_EXPR.INITS}: + y1 = 0 {caption: NO2; category: Initial concentrations; min: 0; max: 0.1} [Initial concentration of NO2] + y2 = 0.2 {caption: NO; category: Initial concentrations; min: 0; max: 0.4} [Initial concentration of NO] + y3 = 0 {caption: O3P; category: Initial concentrations; min: 0; max: 0.1} [Initial concentration of O3P] + y4 = 0.04 {caption: O3; category: Initial concentrations; min: 0; max: 0.1} [Initial concentration of O3] + y5 = 0 {caption: HO2; category: Initial concentrations} [Initial concentration of HO2] + y6 = 0 {caption: OH; category: Initial concentrations} [Initial concentration of OH] + y7 = 0.1 {caption: HCHO; category: Initial concentrations} [Initial concentration of HCHO] + y8 = 0.3 {caption: CO; category: Initial concentrations} [Initial concentration of CO] + y9 = 0.01 {caption: ALD; category: Initial concentrations} [Initial concentration of ALD] + y10 = 0 {caption: MEO2; category: Initial concentrations} [Initial concentration of MEO2] + y11 = 0 {caption: C2O3; category: Initial concentrations} [Initial concentration of C2O3] + y12 = 0 {caption: CO2; category: Initial concentrations} [Initial concentration of CO2] + y13 = 0 {caption: PAN; category: Initial concentrations} [Initial concentration of PAN] + y14 = 0 {caption: CH3O; category: Initial concentrations} [Initial concentration of CH3O] + y15 = 0 {caption: HNO3; category: Initial concentrations} [Initial concentration of HNO3] + y16 = 0 {caption: O1D; category: Initial concentrations} [Initial concentration of O1D] + y17 = 0.007 {caption: SO2; category: Initial concentrations} [Initial concentration of SO2] + y18 = 0 {caption: SO4; category: Initial concentrations} [Initial concentration of SO4] + y19 = 0 {caption: NO3; category: Initial concentrations} [Initial concentration of NO3] + y20 = 0 {caption: N2O5; category: Initial concentrations} [Initial concentration of N2O5] + +${CONTROL_EXPR.OUTPUT}: + t {caption: t, min} + y1 {caption: NO2} + y2 {caption: NO} + y3 {caption: O3P} + y4 {caption: O3} + y5 {caption: HO2} + y6 {caption: OH} + y7 {caption: HCHO} + y8 {caption: CO} + y9 {caption: ALD} + y10 {caption: MEO2} + y11 {caption: C2O3} + y12 {caption: CO2} + y13 {caption: PAN} + y14 {caption: CH3O} + y15 {caption: HNO3} + y16 {caption: O1D} + y17 {caption: SO2} + y18 {caption: SO4} + y19 {caption: NO3} + y20 {caption: N2O5} + +${CONTROL_EXPR.PARAMS}: + k1 = 0.35 {category: Reaction constants} [NO2 -> NO + O3P] + k2 = 26.6 {category: Reaction constants} [NO + O3 -> NO2] + k3 = 1.23e4 {category: Reaction constants} [HO2 + NO -> NO2 + OH] + k4 = 8.6e-4 {category: Reaction constants} [HCHO -> 2 HO2 + CO] + k5 = 8.2e-4 {category: Reaction constants} [HCHO -> CO] + k6 = 1.5e4 {category: Reaction constants} [HCHO + OH -> HO2 + CO] + k7 = 1.3e-4 {category: Reaction constants} [ALD -> MEO2 + HO2 + CO] + k8 = 2.4e4 {category: Reaction constants} [ALD + OH -> C2O3] + k9 = 1.65e4 {category: Reaction constants} [C2O3 + NO-> NO2 + MEO2 + CO2] + k10 = 9e3 {category: Reaction constants} [C2O3 + NO2-> PAN] + k11 = 0.022 {category: Reaction constants} [PAN-> CH3O + NO2] + k12 = 1.2e4 {category: Reaction constants} [MEO2 + NO-> CH3O + NO2] + k13 = 1.88 {category: Reaction constants} [CH3O-> HCHO + HO2] + k14 = 1.63e4 {category: Reaction constants} [NO2 + OH -> HNO3] + k15 = 4.8e6 {category: Reaction constants} [O3P -> O3] + k16 = 3.5e-4 {category: Reaction constants} [O3 -> O1D] + k17 = 0.0175 {category: Reaction constants} [O3 -> O3P] + k18 = 1e8 {category: Reaction constants} [O1D -> 2 OH] + k19 = 4.44e11 {category: Reaction constants} [O1D -> O3P] + k20 = 1240 {category: Reaction constants} [SO2 + OH -> SO4 + HO2] + k21 = 2.1 {category: Reaction constants} [NO3 -> NO] + k22 = 5.78 {category: Reaction constants} [NO3 -> NO2 + O3P] + k23 = 0.0474 {category: Reaction constants} [NO2 + O3 -> NO3] + k24 = 1780 {category: Reaction constants} [NO3 + NO2 -> N2O5] + k25 = 3.12 {category: Reaction constants} [N2O5 -> NO3 + NO2] + +${CONTROL_EXPR.TOL}: 1e-6`; + +/** Initial value problem use cases */ +export enum USE_CASES { + CHEM_REACT = CHEM_REACT_MODEL, + ROBERTSON = ROBERTSON_MODEL, + FERMENTATION = FERMENTATION_MODEL, + PK = PK_MODEL, + PK_PD = PK_PD_MODEL, + ACID_PROD = ACID_PROD_MODEL, + NIMOTUZUMAB = NIMOTUZUMAB_MODEL, + BIOREACTOR = BIOREACTOR_MODEL, + POLLUTION = POLLUTION_MODEL, +} diff --git a/libraries/diff-studio-utils/tsconfig.json b/libraries/diff-studio-utils/tsconfig.json index af1869b185..fa1c08a95a 100644 --- a/libraries/diff-studio-utils/tsconfig.json +++ b/libraries/diff-studio-utils/tsconfig.json @@ -7,7 +7,7 @@ "target": "es6", /* Specify ECMAScript target version: 'ES3' (default), 'ES5', 'ES2015', 'ES2016', 'ES2017', 'ES2018', 'ES2019', 'ES2020', or 'ESNEXT'. */ "module": "es2020", /* Specify module code generation: 'none', 'commonjs', 'amd', 'system', 'umd', 'es2015', 'es2020', or 'ESNext'. */ "lib": ["ES2022", "dom"], /* Specify library files to be included in the compilation. */ - "allowJs": true, /* Allow javascript files to be compiled. */ + // "allowJs": true, /* Allow javascript files to be compiled. */ // "checkJs": true, /* Report errors in .js files. */ // "jsx": "preserve", /* Specify JSX code generation: 'preserve', 'react-native', 'react', 'react-jsx' or 'react-jsxdev'. */ "declaration": true, /* Generates corresponding '.d.ts' file. */ diff --git a/libraries/diff-studio-utils/wasm/get-inv-matr.ts b/libraries/diff-studio-utils/wasm/get-inv-matr.ts deleted file mode 100644 index 2daafc51f4..0000000000 --- a/libraries/diff-studio-utils/wasm/get-inv-matr.ts +++ /dev/null @@ -1,5 +0,0 @@ -import {invMatr} from './matr-api.js'; - -export async function invMatrix(matr: Float64Array) { - await invMatr(matr); -} diff --git a/libraries/diff-studio-utils/wasm/matr-api.d.ts b/libraries/diff-studio-utils/wasm/matr-api.d.ts deleted file mode 100644 index 2e7a306f0a..0000000000 --- a/libraries/diff-studio-utils/wasm/matr-api.d.ts +++ /dev/null @@ -1 +0,0 @@ -export declare async function invMatr(matr: Float64Array): void; \ No newline at end of file diff --git a/libraries/diff-studio-utils/wasm/matr-api.js b/libraries/diff-studio-utils/wasm/matr-api.js deleted file mode 100644 index 04396a8572..0000000000 --- a/libraries/diff-studio-utils/wasm/matr-api.js +++ /dev/null @@ -1,8 +0,0 @@ -import {exportMatrOper} from './matrix-operations.js'; - -export async function invMatr(matr) { - let wasmInstance = await exportMatrOper(); - const mem = wasmInstance._malloc(matr.length * 4); - wasmInstance._free(mem); - console.log('Done!'); -} diff --git a/libraries/diff-studio-utils/wasm/matrix-operations.js b/libraries/diff-studio-utils/wasm/matrix-operations.js deleted file mode 100644 index 18d5f21347..0000000000 --- a/libraries/diff-studio-utils/wasm/matrix-operations.js +++ /dev/null @@ -1,19 +0,0 @@ - -export var exportMatrOper = (() => { - var _scriptDir = typeof document !== 'undefined' && document.currentScript ? document.currentScript.src : undefined; - - return ( -function(exportMatrOper = {}) { -var Module=typeof exportMatrOper!="undefined"?exportMatrOper:{};var readyPromiseResolve,readyPromiseReject;Module["ready"]=new Promise(function(resolve,reject){readyPromiseResolve=resolve;readyPromiseReject=reject});var moduleOverrides=Object.assign({},Module);var arguments_=[];var thisProgram="./this.program";var quit_=(status,toThrow)=>{throw toThrow};var ENVIRONMENT_IS_WEB=typeof window=="object";var ENVIRONMENT_IS_WORKER=typeof importScripts=="function";var ENVIRONMENT_IS_NODE=typeof process=="object"&&typeof process.versions=="object"&&typeof process.versions.node=="string";var scriptDirectory="";function locateFile(path){if(Module["locateFile"]){return Module["locateFile"](path,scriptDirectory)}return scriptDirectory+path}var read_,readAsync,readBinary,setWindowTitle;if(ENVIRONMENT_IS_WEB||ENVIRONMENT_IS_WORKER){if(ENVIRONMENT_IS_WORKER){scriptDirectory=self.location.href}else if(typeof document!="undefined"&&document.currentScript){scriptDirectory=document.currentScript.src}if(_scriptDir){scriptDirectory=_scriptDir}if(scriptDirectory.indexOf("blob:")!==0){scriptDirectory=scriptDirectory.substr(0,scriptDirectory.replace(/[?#].*/,"").lastIndexOf("/")+1)}else{scriptDirectory=""}{read_=url=>{var xhr=new XMLHttpRequest;xhr.open("GET",url,false);xhr.send(null);return xhr.responseText};if(ENVIRONMENT_IS_WORKER){readBinary=url=>{var xhr=new XMLHttpRequest;xhr.open("GET",url,false);xhr.responseType="arraybuffer";xhr.send(null);return new Uint8Array(xhr.response)}}readAsync=(url,onload,onerror)=>{var xhr=new XMLHttpRequest;xhr.open("GET",url,true);xhr.responseType="arraybuffer";xhr.onload=()=>{if(xhr.status==200||xhr.status==0&&xhr.response){onload(xhr.response);return}onerror()};xhr.onerror=onerror;xhr.send(null)}}setWindowTitle=title=>document.title=title}else{}var out=Module["print"]||console.log.bind(console);var err=Module["printErr"]||console.warn.bind(console);Object.assign(Module,moduleOverrides);moduleOverrides=null;if(Module["arguments"])arguments_=Module["arguments"];if(Module["thisProgram"])thisProgram=Module["thisProgram"];if(Module["quit"])quit_=Module["quit"];var wasmBinary;if(Module["wasmBinary"])wasmBinary=Module["wasmBinary"];var noExitRuntime=Module["noExitRuntime"]||true;if(typeof WebAssembly!="object"){abort("no native wasm support detected")}var wasmMemory;var ABORT=false;var EXITSTATUS;var UTF8Decoder=typeof TextDecoder!="undefined"?new TextDecoder("utf8"):undefined;function UTF8ArrayToString(heapOrArray,idx,maxBytesToRead){var endIdx=idx+maxBytesToRead;var endPtr=idx;while(heapOrArray[endPtr]&&!(endPtr>=endIdx))++endPtr;if(endPtr-idx>16&&heapOrArray.buffer&&UTF8Decoder){return UTF8Decoder.decode(heapOrArray.subarray(idx,endPtr))}var str="";while(idx>10,56320|ch&1023)}}return str}function UTF8ToString(ptr,maxBytesToRead){return ptr?UTF8ArrayToString(HEAPU8,ptr,maxBytesToRead):""}function stringToUTF8Array(str,heap,outIdx,maxBytesToWrite){if(!(maxBytesToWrite>0))return 0;var startIdx=outIdx;var endIdx=outIdx+maxBytesToWrite-1;for(var i=0;i=55296&&u<=57343){var u1=str.charCodeAt(++i);u=65536+((u&1023)<<10)|u1&1023}if(u<=127){if(outIdx>=endIdx)break;heap[outIdx++]=u}else if(u<=2047){if(outIdx+1>=endIdx)break;heap[outIdx++]=192|u>>6;heap[outIdx++]=128|u&63}else if(u<=65535){if(outIdx+2>=endIdx)break;heap[outIdx++]=224|u>>12;heap[outIdx++]=128|u>>6&63;heap[outIdx++]=128|u&63}else{if(outIdx+3>=endIdx)break;heap[outIdx++]=240|u>>18;heap[outIdx++]=128|u>>12&63;heap[outIdx++]=128|u>>6&63;heap[outIdx++]=128|u&63}}heap[outIdx]=0;return outIdx-startIdx}function stringToUTF8(str,outPtr,maxBytesToWrite){return stringToUTF8Array(str,HEAPU8,outPtr,maxBytesToWrite)}var buffer,HEAP8,HEAPU8,HEAP16,HEAPU16,HEAP32,HEAPU32,HEAPF32,HEAPF64;function updateGlobalBufferAndViews(buf){buffer=buf;Module["HEAP8"]=HEAP8=new Int8Array(buf);Module["HEAP16"]=HEAP16=new Int16Array(buf);Module["HEAP32"]=HEAP32=new Int32Array(buf);Module["HEAPU8"]=HEAPU8=new Uint8Array(buf);Module["HEAPU16"]=HEAPU16=new Uint16Array(buf);Module["HEAPU32"]=HEAPU32=new Uint32Array(buf);Module["HEAPF32"]=HEAPF32=new Float32Array(buf);Module["HEAPF64"]=HEAPF64=new Float64Array(buf)}var INITIAL_MEMORY=Module["INITIAL_MEMORY"]||2147483648;var wasmTable;var __ATPRERUN__=[];var __ATINIT__=[];var __ATPOSTRUN__=[];var runtimeInitialized=false;function preRun(){if(Module["preRun"]){if(typeof Module["preRun"]=="function")Module["preRun"]=[Module["preRun"]];while(Module["preRun"].length){addOnPreRun(Module["preRun"].shift())}}callRuntimeCallbacks(__ATPRERUN__)}function initRuntime(){runtimeInitialized=true;callRuntimeCallbacks(__ATINIT__)}function postRun(){if(Module["postRun"]){if(typeof Module["postRun"]=="function")Module["postRun"]=[Module["postRun"]];while(Module["postRun"].length){addOnPostRun(Module["postRun"].shift())}}callRuntimeCallbacks(__ATPOSTRUN__)}function addOnPreRun(cb){__ATPRERUN__.unshift(cb)}function addOnInit(cb){__ATINIT__.unshift(cb)}function addOnPostRun(cb){__ATPOSTRUN__.unshift(cb)}var runDependencies=0;var runDependencyWatcher=null;var dependenciesFulfilled=null;function addRunDependency(id){runDependencies++;if(Module["monitorRunDependencies"]){Module["monitorRunDependencies"](runDependencies)}}function removeRunDependency(id){runDependencies--;if(Module["monitorRunDependencies"]){Module["monitorRunDependencies"](runDependencies)}if(runDependencies==0){if(runDependencyWatcher!==null){clearInterval(runDependencyWatcher);runDependencyWatcher=null}if(dependenciesFulfilled){var callback=dependenciesFulfilled;dependenciesFulfilled=null;callback()}}}function abort(what){if(Module["onAbort"]){Module["onAbort"](what)}what="Aborted("+what+")";err(what);ABORT=true;EXITSTATUS=1;what+=". Build with -sASSERTIONS for more info.";var e=new WebAssembly.RuntimeError(what);readyPromiseReject(e);throw e}var dataURIPrefix="data:application/octet-stream;base64,";function isDataURI(filename){return filename.startsWith(dataURIPrefix)}var wasmBinaryFile;wasmBinaryFile="matrix-operations.wasm";if(!isDataURI(wasmBinaryFile)){wasmBinaryFile=locateFile(wasmBinaryFile)}function getBinary(file){try{if(file==wasmBinaryFile&&wasmBinary){return new Uint8Array(wasmBinary)}if(readBinary){return readBinary(file)}throw"both async and sync fetching of the wasm failed"}catch(err){abort(err)}}function getBinaryPromise(){if(!wasmBinary&&(ENVIRONMENT_IS_WEB||ENVIRONMENT_IS_WORKER)){if(typeof fetch=="function"){return fetch(wasmBinaryFile,{credentials:"same-origin"}).then(function(response){if(!response["ok"]){throw"failed to load wasm binary file at '"+wasmBinaryFile+"'"}return response["arrayBuffer"]()}).catch(function(){return getBinary(wasmBinaryFile)})}}return Promise.resolve().then(function(){return getBinary(wasmBinaryFile)})}function createWasm(){var info={"a":asmLibraryArg};function receiveInstance(instance,module){var exports=instance.exports;Module["asm"]=exports;wasmMemory=Module["asm"]["e"];updateGlobalBufferAndViews(wasmMemory.buffer);wasmTable=Module["asm"]["k"];addOnInit(Module["asm"]["f"]);removeRunDependency("wasm-instantiate")}addRunDependency("wasm-instantiate");function receiveInstantiationResult(result){receiveInstance(result["instance"])}function instantiateArrayBuffer(receiver){return getBinaryPromise().then(function(binary){return WebAssembly.instantiate(binary,info)}).then(function(instance){return instance}).then(receiver,function(reason){err("failed to asynchronously prepare wasm: "+reason);abort(reason)})}function instantiateAsync(){if(!wasmBinary&&typeof WebAssembly.instantiateStreaming=="function"&&!isDataURI(wasmBinaryFile)&&typeof fetch=="function"){return fetch(wasmBinaryFile,{credentials:"same-origin"}).then(function(response){var result=WebAssembly.instantiateStreaming(response,info);return result.then(receiveInstantiationResult,function(reason){err("wasm streaming compile failed: "+reason);err("falling back to ArrayBuffer instantiation");return instantiateArrayBuffer(receiveInstantiationResult)})})}else{return instantiateArrayBuffer(receiveInstantiationResult)}}if(Module["instantiateWasm"]){try{var exports=Module["instantiateWasm"](info,receiveInstance);return exports}catch(e){err("Module.instantiateWasm callback failed with error: "+e);readyPromiseReject(e)}}instantiateAsync().catch(readyPromiseReject);return{}}function callRuntimeCallbacks(callbacks){while(callbacks.length>0){callbacks.shift()(Module)}}function ___assert_fail(condition,filename,line,func){abort("Assertion failed: "+UTF8ToString(condition)+", at: "+[filename?UTF8ToString(filename):"unknown filename",line,func?UTF8ToString(func):"unknown function"])}function ExceptionInfo(excPtr){this.excPtr=excPtr;this.ptr=excPtr-24;this.set_type=function(type){HEAPU32[this.ptr+4>>2]=type};this.get_type=function(){return HEAPU32[this.ptr+4>>2]};this.set_destructor=function(destructor){HEAPU32[this.ptr+8>>2]=destructor};this.get_destructor=function(){return HEAPU32[this.ptr+8>>2]};this.set_refcount=function(refcount){HEAP32[this.ptr>>2]=refcount};this.set_caught=function(caught){caught=caught?1:0;HEAP8[this.ptr+12>>0]=caught};this.get_caught=function(){return HEAP8[this.ptr+12>>0]!=0};this.set_rethrown=function(rethrown){rethrown=rethrown?1:0;HEAP8[this.ptr+13>>0]=rethrown};this.get_rethrown=function(){return HEAP8[this.ptr+13>>0]!=0};this.init=function(type,destructor){this.set_adjusted_ptr(0);this.set_type(type);this.set_destructor(destructor);this.set_refcount(0);this.set_caught(false);this.set_rethrown(false)};this.add_ref=function(){var value=HEAP32[this.ptr>>2];HEAP32[this.ptr>>2]=value+1};this.release_ref=function(){var prev=HEAP32[this.ptr>>2];HEAP32[this.ptr>>2]=prev-1;return prev===1};this.set_adjusted_ptr=function(adjustedPtr){HEAPU32[this.ptr+16>>2]=adjustedPtr};this.get_adjusted_ptr=function(){return HEAPU32[this.ptr+16>>2]};this.get_exception_ptr=function(){var isPointer=___cxa_is_pointer_type(this.get_type());if(isPointer){return HEAPU32[this.excPtr>>2]}var adjusted=this.get_adjusted_ptr();if(adjusted!==0)return adjusted;return this.excPtr}}var exceptionLast=0;var uncaughtExceptionCount=0;function ___cxa_throw(ptr,type,destructor){var info=new ExceptionInfo(ptr);info.init(type,destructor);exceptionLast=ptr;uncaughtExceptionCount++;throw ptr}function _emscripten_memcpy_big(dest,src,num){HEAPU8.copyWithin(dest,src,src+num)}function getHeapMax(){return 2147483648}function emscripten_realloc_buffer(size){try{wasmMemory.grow(size-buffer.byteLength+65535>>>16);updateGlobalBufferAndViews(wasmMemory.buffer);return 1}catch(e){}}function _emscripten_resize_heap(requestedSize){var oldSize=HEAPU8.length;requestedSize=requestedSize>>>0;var maxHeapSize=getHeapMax();if(requestedSize>maxHeapSize){return false}let alignUp=(x,multiple)=>x+(multiple-x%multiple)%multiple;for(var cutDown=1;cutDown<=4;cutDown*=2){var overGrownHeapSize=oldSize*(1+.2/cutDown);overGrownHeapSize=Math.min(overGrownHeapSize,requestedSize+100663296);var newSize=Math.min(maxHeapSize,alignUp(Math.max(requestedSize,overGrownHeapSize),65536));var replacement=emscripten_realloc_buffer(newSize);if(replacement){return true}}return false}function getCFunc(ident){var func=Module["_"+ident];return func}function writeArrayToMemory(array,buffer){HEAP8.set(array,buffer)}function ccall(ident,returnType,argTypes,args,opts){var toC={"string":str=>{var ret=0;if(str!==null&&str!==undefined&&str!==0){var len=(str.length<<2)+1;ret=stackAlloc(len);stringToUTF8(str,ret,len)}return ret},"array":arr=>{var ret=stackAlloc(arr.length);writeArrayToMemory(arr,ret);return ret}};function convertReturnValue(ret){if(returnType==="string"){return UTF8ToString(ret)}if(returnType==="boolean")return Boolean(ret);return ret}var func=getCFunc(ident);var cArgs=[];var stack=0;if(args){for(var i=0;itype==="number"||type==="boolean");var numericRet=returnType!=="string";if(numericRet&&numericArgs&&!opts){return getCFunc(ident)}return function(){return ccall(ident,returnType,argTypes,arguments,opts)}}var asmLibraryArg={"a":___assert_fail,"b":___cxa_throw,"d":_emscripten_memcpy_big,"c":_emscripten_resize_heap};var asm=createWasm();var ___wasm_call_ctors=Module["___wasm_call_ctors"]=function(){return(___wasm_call_ctors=Module["___wasm_call_ctors"]=Module["asm"]["f"]).apply(null,arguments)};var _invMatrixD=Module["_invMatrixD"]=function(){return(_invMatrixD=Module["_invMatrixD"]=Module["asm"]["g"]).apply(null,arguments)};var _invMatrixF=Module["_invMatrixF"]=function(){return(_invMatrixF=Module["_invMatrixF"]=Module["asm"]["h"]).apply(null,arguments)};var _free=Module["_free"]=function(){return(_free=Module["_free"]=Module["asm"]["i"]).apply(null,arguments)};var _malloc=Module["_malloc"]=function(){return(_malloc=Module["_malloc"]=Module["asm"]["j"]).apply(null,arguments)};var stackSave=Module["stackSave"]=function(){return(stackSave=Module["stackSave"]=Module["asm"]["l"]).apply(null,arguments)};var stackRestore=Module["stackRestore"]=function(){return(stackRestore=Module["stackRestore"]=Module["asm"]["m"]).apply(null,arguments)};var stackAlloc=Module["stackAlloc"]=function(){return(stackAlloc=Module["stackAlloc"]=Module["asm"]["n"]).apply(null,arguments)};var ___cxa_is_pointer_type=Module["___cxa_is_pointer_type"]=function(){return(___cxa_is_pointer_type=Module["___cxa_is_pointer_type"]=Module["asm"]["o"]).apply(null,arguments)};Module["ccall"]=ccall;Module["cwrap"]=cwrap;var calledRun;dependenciesFulfilled=function runCaller(){if(!calledRun)run();if(!calledRun)dependenciesFulfilled=runCaller};function run(args){args=args||arguments_;if(runDependencies>0){return}preRun();if(runDependencies>0){return}function doRun(){if(calledRun)return;calledRun=true;Module["calledRun"]=true;if(ABORT)return;initRuntime();readyPromiseResolve(Module);if(Module["onRuntimeInitialized"])Module["onRuntimeInitialized"]();postRun()}if(Module["setStatus"]){Module["setStatus"]("Running...");setTimeout(function(){setTimeout(function(){Module["setStatus"]("")},1);doRun()},1)}else{doRun()}}if(Module["preInit"]){if(typeof Module["preInit"]=="function")Module["preInit"]=[Module["preInit"]];while(Module["preInit"].length>0){Module["preInit"].pop()()}}run(); - - - return exportMatrOper.ready -} -); -})(); -if (typeof exports === 'object' && typeof module === 'object') - module.exports = exportMatrOper; -else if (typeof define === 'function' && define['amd']) - define([], function() { return exportMatrOper; }); -else if (typeof exports === 'object') - exports["exportMatrOper"] = exportMatrOper; diff --git a/libraries/diff-studio-utils/wasm/matrix-operations.wasm b/libraries/diff-studio-utils/wasm/matrix-operations.wasm deleted file mode 100644 index dc23d97c2257dee45a08987b352df1dc747d3382..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 59309 zcmeIbdyHMledl-1xwpG--$(QIi)^wvWSv8Dra9!yP#nHQ4J96F9vo?AG?3*$!AcnQ zkUg4y@YUTj8ac8^BVI#T`H#(5&c?9=%*t8YF=RwAHrP!TS!oex;{*mQ>}mrASp-Vr zEc}N92#JF@8*QA=_g8h!J@?UMb2Q^fiwuW6r|Q(>SHD;N>R0uv8!c~LkE1AxH}`zL zeuw{}&yVY&m-tS6=kvAj8rd5eW8rI}A`uz9@b83W_q|4+Z~4rh*P~2**{7l?*f zKU(7Tc{2IDmPYyA|zaCLoXXSHYRm-9bc z;(sdUUvcplZ=Xw|cCAH$+PS3dYOa1RX|$7;$}gnjz+(Wkvsk@Y>opgvXsjBIHuxkz zpIgi7g0lwbN#54d8RRJCYK(G+l(lCG+@+(aG$QOk`pb~U_|9L(x@ zvYr3NiM6a|!(|M>2@GJN6X!o%;(tG!dO4i&yqnLtahi3>vx{q?^fMg>7Q69t$%0y@ zOJPpq==v5HY}|8MO_;Y2bKwI+aAQkJT}Zv1*PiP{`JXR2!Rp4q_E*3oIEJ9{rrtop zntk0vi*D?bU?tR!p=!j);yk(i4Al5Vqv4nc^NYt_LWB8VEb+gF8m(p4L$BS$)BL$k zjqVLqWY}|0E3kt#bvHCcp+DcDsVa!#KYB|&V4jnJcQhyQVMk5^zx|x#(dId6OgIRf zWUbc%M<4J=?&qnoK1LaPZ+J=rr0YSR63*_;Qw*%aQ)5-0Vu*pKXn7=0k#9T&!;0Ao zo)T*eF$~gUV*MPf9XJdgx+jmdXvNbkf{v2-G*0e#AI!C#x7YIz zBgz^RD&(L`J{p}ztP0spCz8wh<7L&Ls&8iq1APC(C3^QDO_&7G#so4d@S=P_ZNzV` z=Mm4t&n<2P$ppfY#;BB7Dxq%7`)r4f^M~F{W0hAabVOdiy_QYDI#kGF(=X}(r64*; z+Kb_8@am8D8d*R>O&=Q&1SFs^R7IE(N?Q3ZALoCmop?Qq@;feC1f{Lk?=QuzJ4d3t zyPma-oiOOqx3~<$b#j!WhJzE#Vktkf6{}KoA%Z zMeDa;x2cAdAQMjq$f7~Wm~S$)ZzN2kE}e%NSbVIDvnbm9lVwWnV;(R)iH|$uSBI`^GX|DzmtxElyBep1jWRyQ2MRg=~zHV zg9!~}Xw-qz?ML!XbC9>4K#jsk&okx^pQww zfPIFp3iVA70}9MdpgePN`lf5Z^unc2SIM7(GNW$oa+SWQjptK36`4JSF?SZFPObm? z9T$UKplcb{kqg;}n%+SKNY$92$pXexT`%p->^!E2q4ihO~5-~>&4DQ{vWFCtWGb#VQ&-I^Vy_B z**3t}~}dSel>a~^eD z(CE&y&>yNZ5_NA=7t~>zkyNaM!TK7IqR`h;5!{V`5~*&QJS>tEst`lpt9EmNwR)B zOST9w;N{nkTwpatnv_;R0 z(k4Ih@?knH|JF$@L?YVrLW4nZN9n(L+_RwG|Gwv=8W{3ZKNBmQtZCL{E6(2`!(O4L z^s-}R(vmog=Swn7-eS1<#BCuHcIRyoo6Ra~xdA;8FnZ#xH_Ik#`4Wqm{Hw?F?+F=v zXFx7`ENbVC+l-U>>JvO5)QcK5JUG6l3ba_G$6DyIR`iJJVuD#NY```3fX{ag>?ZFt zIHK#C@KuMgc!>5Mvo85W1*5}VLRRzybFs)0Y2pmziCNy3qtSF+kQA*>$t)H} zTCF;0)=+~mF_%OwX?Y_U4~mDOgWX|FV|o^P4x+t<3h#MQJOc%O?9{T+eD+#m2Vk(NsQ}T!wOf z4Z-&YY)y$bu@FiR0^(|ukGC^R4?qEJGRA6!bWk#hS(K!(8u! zn}p<`A9y3)r3$fHg7aVmxuE!t>$)U|MF-b)$)Kex!G=4uuH(_xI?_MDQA4tHSg`D3 zVF%-kX#sT5jroOLEub`wlVVXPkg(Zp!qTFS`NYFx3@V#ntGQaR<7(cHV@84|X@aCC zvRaGPT#bUP=1>W;%P0k=%pdJtpY^|rW{KY=C#!!=lu6G|K!5RjeVySOdKZCrtM-GJF_8_tBrNamA;6m@7U2Q7zvi z-;w+>EDcaykX05WlWb%&X_(zm6ZLZ11WHKoEr_YhE*AZ+MSaKT~V=s`vmP z!ziXkUyG^3WCXmy6`<+Mn3`UK%c}-cW6z9cY8z|m@+#0R&B(Y>DSZq*?izuiHKj-^ z)^=p*F&KKRWM~9lY5<1D3@#XYyvNXEo}pmVmV45@+!y2LaBWlJx(Kl?mXF=Fqk!9LMfvxl?pXfaXsxv! z$9G(!CA%4yP^zL8=V3l4dS%PV zBuVCy+Jjqs!>zv_(<_YuIgJNXnY|pKsQDeLuIztirs#z_&K-H9Z6W+gEe51rx@%qd zl>HTww{B;Jq9ea$V>bOz-_z=jBl$miggiL-Ps8(DOL^Rlqt+j{lKDIBpZDI=gP;br z`Eg)~(^;3KbFje#s_D+Ni|IV?4J1qNdm(WBIzY%nf;Pz78tFqg=*T)m!8S@2b2r@w zZ1L^%uxsN{h4|@$4L3!5Y4Vr`YwSL$!&%pH^BY)bGy&4ks(bi2=k~b08}5+LZG9}_ zzcdB1hyt-s9`%o7{&B=Vmi*(me?00RCjdj8(eMcuJqk@9cCAN~gI|AJsdn~=QcwBR zT$U;IS)V$XIi6B}G;~1WtGPQu@1_T>AK(L|B5}?g+V+n9$2(SbQ0O*HW zICmFN25{C9Hzg`%E{3H*2B|ufAO)XUllF-w!8^IIk+p5PwC!_Y8tcSz>BQ&4BJ?&T zWi|yaL2ZmgqJi-Z$qCVNGxVK0w;uwHf8FgT_4Wk>e}?);$+5%%UuF@epv+=X<^W}m z6lL~c`DT0Ega{>be`17;x{0(=lZ-Ybqs<_rJ;5|N?TH-BKvxnYjfeNT-Zf$|vfhOg z-rIT?>Ui(!onh;JPVXQ_?+bbdb$UOjcOjP^3U&iE{oKatvX7&Bi~+!K1DM=rngN9+ zwDv(u%PlPhP{@8oFh!wk%I)(=U&w2$4?sm)&omp?~=&Nuk$9Ds<(xo{qwQAR;LBUe?jx>#^SJ zk>2Z4@AY`^_0isIzrsIJYe;#Xs9pHRJ1#m=TafZRQG3g$mZUsS)c$_#lNY2tPt^W` zWTEm;)K2lu>fGYRihoPRt>h0?GAR8MwWXru>qW`mDNFvTN(QxmqUMT{pD#+jQGvOKQ;Jfu)oZU!TvYPe!i=cLHUEF zqU4pLDw)FGm*r_=0!Ri$QWJo>Wdcx}OaQ9D1OU{% z37|GS0r;%G3BYIdO#nWtC+9m%0H#Y@UY9Zf2+>AQ5M4|FhCIWMeE~nk1Yjrt=z$49 zji9+bwWh>(sCwhOVgm5=uw0nR+O`V9z2z>2PApe=wp>T9CIEDWCIC=7Gy%|eF#&wT z?I-ol2W#a7@C{2G@MUZQ_=a07%7h8Pm$}abU@|>40hml{0_ch0@FYNYk}Qo0Y@_%+ z$BC)+EeBPw}vKxZ@3v6HR|nR68NN6Qs2{VCUXfcEuD49kiRHl zm6As;a`B6cTDp9M*=HN4zKgnvx#hoyOGGOLRz0PI<4)_B;<#pu6rjV|&RwAAF2i|5@pX+JUmTxKheOhv8>HU@45HZohg!O>;U{R8-3=;fh z>p$1gyBcUg$XJgs%W9KwVuKwIrJ0XqRVh~sAX%K159}FZlgV$cN9-nLp-sA$pQqcP zFEl%{QQ03X_F5QJG{NBtZCUs^9-91T^Z7s6lmEkR{$DS4qv)>~6Nq;oCP>4SQZECT zC_aF+U0{OJ6-*R$g`)?+#Q$E4n|I*)KUrd74ry7r|3CtePUipIp3VIGddkgo_pT%R z2TOW^+4=C%yLW$~+-Clhr98e!BP!P9-?*idRJLM;6yhQn91lvAjlmRSALjs(!$)yt zF+T+#auQ-Jmlwlv94UV%zAo95XFYylF{{apj?%V#pw!~u>sftcJ!7jESwKu!CbPMy zEuhBQbsp;}-c9xpOeQrPo%T}K4W={m+HI&*7O^59+kQ!{%lWc>!;kYFc^}5MZ?v)+ z1Nniz)wgd@8bSQ{V*X_y6GZH;@CR7%UKliYJT&XZ{+Og%>#wsJA5a^rCx|*A7lk{B!c-THc%-xAwMR{yCL9XlMHL2A4Yj_2{|e zsDE2JmuS0(uP4tXC;c1#=u`gf%(>)@f4gulx!~WPI+r}f$joR+c;MvaM1j^00Lz)P zKIpX&F$^#jF+z=8NUP;_cIgWg!bSM17?xUU&)e z1}8%*T3)$2ydW)K`5=cmT!sqpDlV+UC%u+$3M2Z2ZX4AolZKt#6hlkutXAuv)#5gy z7kq2YM4VQaPD%U_`^E->bUtZ(L!F>f1SOhd)2BNjjf2F#PlYt2u=F0GfF`68H1L7?tr)1o8_T&4R2yd2^GF z%HsQe|K}#A4Q3n7I>p?S^vzAoMwPjVd5h9$mG^LNf;E1OxvBlG=ce~mJ7E;{X{Vl% zG@zY|4>*Z-p`9p=S}0UrQCE8FJ=9J#+ozqBu4*TxtJ+EFs&*nBv=e4VtFSn|b~5wR zE2Js^r5m)9ePPG=ebr6^bw}+axV?4)?(qWdQX2)(eY(3@K=0EFA-%9m?L=Ruop=U| zeD~U^7R9aKtocAGsb&nG$O|hobytUpv_pkaYl|z}x2dCvwaM3z*)@iNn~tL<0bS}* zD?M0oKoY$uoX zhg@p_Lp;s9FJLv_VcX?TiBj56ukX-n#1`S6rQ}vRt;YuQF?_*rtoAgrH|eIWHoKOW zImB*8BMablwgM+XUV2EG$>kcog+6KLbvv%B2D2{c^sv(l2F1b)of)=uCfx$|>N*Rl z#1Jl-U(gZNLwrDekYeeS;r5pqghL+TQ_K%-649U)%33y$Xr!kZC8%e4pU3zm@ZsEM ztMl})rKiQyeB4L7pzd_~r0m@(ISXc#;&TTkpLaFQ?R%YXlqP>IJ3}t`v|cRZtb39_ zXRl`M#Y}H!u4c`}l&3rEPP_J^BX!#C1->Mtv`*8S#dUbF$e!^mI}0QOZJ(h9YjIv} zs0ChFOVt^*NuOJQ`NUzTOE26qh+<_WcLK+Uh2hSE(GSxe1cN9su@)dTQXqjcGT3x< zf5C5J{f=p&^{1I4<#E4HvX+-Wnk)^ye>_O!<`R$j-=dMr1vJu{tMguVgz4(6MO3eeey z%sQc+PSFlJqvqzt3-hiK6pZ5BkSQHC)9KxUBrNHw)Ob~gTk&t=PSj?4i!)E{`Bwh- z4(5NOlmBl?{vVk1o0y3@@+Cwsp_Mt^u=q@87HOD}AWWDbuoix%Lyp$XI`$YSG%s=i zNg$sM#Rvv_>#%|GwT$Y-c)F+uPBc(=`Zq-ZHVJqx1#X1JPhNPFbR$x$MP*H z+tZ*AU6ERnXxB70b{OG2m~wov7#J?KFfcJ|&6su8Gi%3JnRZUA@?JfJU>z!W%`oQB zfaj85j{Wm4Pid$xlY|a$@z{^P<2cyM;>q*BnvJS8-TkJZheOTH=)=}bqtDAndM zj-233Gx~Z3`a_YWJ5@U?q&%(Vt zmi$qSIpt3&<5~(-I)KM&Sub6wX7mN{^XyDB3_Ht9u2@rHxEvFK)o!P*#^`CtCDPeV zL}|+FjtN*6MOWxYQ?=A*cBGznw4$|z?D@59*3BBn2>&u>zVWxkiW3k6WwGP;acNI{ zQR2sjoAdJM;I+^&r{6iqQZ#JVm^EyEV6m_*_|W+687zeqU5a2YD7d25YPH)ojcEqd z=^-VE_em{S(ocWlfh*|OlrwPo4>`aV-_9^T+{B{H5B^Mid@*(I&@&mlz@)pGKE^z_ zAcVPju{+347CUW5NjqZ6c^>5yNmXE)R}jUxbIh*|aFb24h;va{TN6H3u?Y{WYu8$L zlUF)ARb4bcjg!xMBNt9yS0^f2yzLQpe+C1A$ESHPwLlbMZr&Z*fQ`(Ml>ktLnHP)c zS$A@iSxzFtHb)5ye6q*kqmv2JqNGYZZ7oh(i~I4HOG>P(r-X#{4Cpk#MzxV~HyOr* zY=rTkGZ^O~agzzBqZN*Q9rcYZ=|PV;&JDH+dM4^)5H>2bDG7)iEHMyOP{)nhZa<(k zoR9@ufM`JsH`c+1CRpu^<=yW7^l?^BQutC%b?V^tbn_q2W*mJj2U@1 z#r%xi-|4`NjH^bHJH>+&6~XTDM==7}Q+;zWDRxv1>?sen_DE=Gpx z0`=`~%TSD-Fosq*vrvrQ6^cPB4m6=yS{Q`mp;rU&pZTynJ!3bk>5ZV6Ag@zMrE^bg zpd;nCcO=JqLOQVyrN#Y_yeS>PC((VMsBs$$8oDQ!0gE2xg%1n>WTta1%m=_^l+vZ_ zF(NWu)04eHV!TM8;45UG%WJI5XqV%WUPv;jROmVDZ5gFRcI;j~8l2HVx6W4VHB2Kk z@o+K7zMEo@9ec08({^aBT@bBG_HdA9+S_0~xJ7B;W@}b9*iaKyx=(wZPeh+l{^LP) zqrC~OMnXRWF?oENdmKO1A9GH)!Gxl!%qM&D3m{ruj<2odv)W%pnTYui8AkD3rwQk=g>h2cC|iaHEUTRd zd*jmu-e``Uho}78)6)L%tE?r4PwUThv}Ndl>-a)oc)`D6IXT`qsXWI6^QNl@(SxFa z-WH^x;dvxhi%GZNB8>&~l!sNAmjz=&qk%D^NF!|{eUxq}Ca96l8p2V;e2KEXIT~Ew zG-EDyx|)6oY0|Vx)9Ivb9gHJasFjz%uVJa7PfW5jKy|5x#SRMsnvqiB{Fm6xtCFTT zG_a?eQbN%VbSB}~Txvo);bu4RkL{(GeO}Z}q0_CPkZu|q>N%S(``%mjJq9$KOs=#1 z8rMlm&E#zJAyKZYf2M(pJYKztwKJ43A>c`SVeBnxIJeUSWF^n02Z5_4h#vytRu$q- z0dY%c@Ej#vVGjw}f>;Vy*b*&VN5v+zl&BT($>3ParUfBNfKu&r&Zw?1>Vpu5E>2@x z7UR&6!c3?-@7O1__~9^4zR$VChE)5mM!nt0%&XaA(d&={(Boe1rl(s>sSgRITCgJT zH0Ms6p&Ab5d0ImN`E5&Di-Kax*y96Zxt_+dnZ-M?VPnRuZ(j(@41n6xt^B_~kpJze z{LiEO|Cy90C>W&s*|xU-{^Hwr&exDcL?j(ASN-W}#GReVslrp_WgM%NerZ_Au0-Bc zLTHr zsa1VSRqMVg&8wb`H&k8lRl5{GfwXB)Cpf)oZU}w@+GNDAjWb1BCu5)4LtXIFlyaF5 z9WUsxg2NmeGh=E+L!#^pu-&AnrsXc&g)<1{PWybOChV$p5*h1wq+s_Z?{*15S27@p zl*x7)i?3uptd-sx2jJ^~}dj7x2D zT<$86O$4YEDR0M?E;ooKF9x)JS zr)L}{#kOvy0)umG>RY(3=DqG%#;sBED{<-v1c-DUolct;G?Gk??%%qNanEcuXAvmU zl;}4e9?8Ig(1K{QT6mL%1a5>h^Yox0<~j6y>#Q0HR=ZAWXw>y@E&CWZlFKH^+DzI6Vm_2r=GGC7!s?678 z2g09&;r0&m;DfQlS{1cJ?9lghh#mUAR@osOnwa%#<$anpL7U840uygrU>I}gVNa=9 zBax_8rBX!$*~W{h5|P{Ees2%=r~B1`NWq3m;YidXw?z0^lk2@bblFEHj{^ls&QoHq ztCB=$_Vs2e^;B-AubdEXKF3Abj*IRgKYbte)AzlANJn;&K6!sx?jt~zc<4`Uk{|R0 zm~{Ij7lwO7jm&i~C8Kd)3cw(m!R>4nw>@Sa!mK`l*@FSIhSw@nm$vq-V4W>(ZI7>= z8}`O_pN{JvIm(a$d+e9S$#if1+{GR_8p#I~-8&Xis~jz>g$=(PlLb3hSg?(ZXul51FMdqzu+2$F`2yIE4 zE%u}=Jlf06PGub%EL_>zW)1W5R90BS*tFj6#Yzhc`FS0bXQ8Wk&g@>MK5zF@^Jy$I zQtKEM2(0WkGkV*SV)+)wQE8{OgfFO!yTx@1`wOEttZ$pW@6GCWSYJymEds>_SazOn zUC#wOH!Ln-;BJ)Bs}T!sMnwB6rDI~UeawY(f!@s3fER)@Zr+6hePwmANQW>COoQ7i z%1CKp!RCCEafQWY#?|wSYlzk(|fI0Z12|efnIZ7IZSMt zVjka2A3?;&Hq(qpeKTF+QQJ(9^GG(+6FlP0^ielq#6RP3c)0@w#hKEk-9zN;W9Po7 zE?F_ibeYDgeWa2r#p-jjO72-RlvH-9l6&5aD&-!pWRE(3MT`Dj6TDE~XgL#gJ~WxPnjU+3TGyX(R*qe0qU2H`6c*pLV1j zX=ZV=T$V^LkuDZ%r%0b7%|^a`zd-r|SjH9CbjLhYO?MO`^$34|9_Z1vwgwwPyEW`@ zEl8kpK|8c8zM3GLM!N&#>`1$+`tH%LYLgVS>&Nl=9dTUX$VamDa8%&9z)^wY0!JRl z#-4a{o9@VNu!FN_A{$4bk-W3a-qo?6*D{ikD{NKdmq!e<|14BFFZUBsNXEpLASkP6 zXASX9LFZyj#Tbe)6=U#Yg4@zN(#U%*zt zTEOOEHKqZ!IZUEOHZ|5^BuR@H$+h~yWig2+d5KvC=5hy^E>N5^4J3nUsCY5zElJhU&n3i0h0OiU@Ze5##zQrYA(M5_Mbn{XH@G&TIkIlkrK); z)m$EH+3~)SQtBahyjRLtbo96}I4g@NjT==H;*s`>4mx_%G1Vi}$LuAHpkhXq*&?SN zrpn2}{xMIV8x?TV9UXAfeV~Aw4p?OhKPB%erxY@CY@I|)vz2mX(Vnd$;HHHU6f-@u zi!+@s*_md}p;^3*twC5QQi0)I&l-f2w4spBsp(1z0*?@XRIiV?G@ueOV$f(WtL-MQ1JWeW zGzbS1H0;4+5^dvp3b8hmn2aZ7iUK-Dl4~GlgaV1)bkG-UgFEs$2K);aY{SK>9R=Iy zt_};fVeVQK!?6i~C8A&(Yg1U1)dg8)27wzw2{=&zJi;=k2e6SW*k)V?uP|c-VXGG0 zZ0K+cxH1Hhusy8Hg}2KtGHUoUf)|*7%(4tFuqAVVVKI@}z+&Ptm+tkR-U<0!3vadO zfF?Cs5Iy#|CMeIOjo4y3c@A3Z?d zN7Q$NN7;9y=zG?BKH%3F)-x_?9r2=PvrQF-LC^Tb*~Uu`Fx$ynJ?7Xj;Q`b@q1QSr zIrE$_^Ahe7>@`$q!#%A1hNjf@$FmcQa8-JcD)JR_v+^^#UzcKSw@=1z8~tG#nyV3+ z1HB8bb##Mk9>chRs$D?U z0V@7p$xaEX>;qltzWa~p1AXv!fd5nhzXk#yL6{S~aUGfduiLS{M1wM_uyGGEnpd(( zql$Ly$NL9t=132jBEz`I(1Hy7y^_rt8arK9aCASX?zw@(O{#kcd_<0U;Rh?;OnkRV)7|QT|JIvbn|96*n;*WkmDtaSko&G%wdgFSy4!`nZ@bxFk;hvXkTl>b&c|DU4#w+)i~+xztWUq;3EzcZJAZz}(9qWqVAdD{JZ zoYEZ2{}``dBlnx}^L_Uaotf2ln`T`BAc zl?nuMwguRC*;C>#dt!z%P0_pTsqU}6Qh>EQ3p&ol@ju$vdBXeM`SkC{kq)6NY)Hlw}P3n|i?a;0=A%V=!@*C)CR~I|Gl#ILuPG^ zX)W;xRVux-{?mjwcBUFOMHBdjPf?}_MxH4#>-<2B-Z)OGZA|AM1X=n`O+BLXKE(}c zI#kMu>HIs#gcR_hNrMP6PX>&9%`i1Ug;wZ19qM54nkf~0z-bZs&-tQz{QzWc+Ni=X zm~cQ6_m*ppD9Tn!E=-R{$z>X*6Jb|tDN3Xc$gO`AE9RECh=vZB+8;}D38WC}uR0VCmEs5;7(1 zqA29BvC9QGgi=Er5+G?dnw`<&m>JfC=3;zhA?8wgT z60v(rXITigjF=L1GRp!GHO~1m^YwIGn@nv|U@e6p z-tTmhNAASBQt#AyHi0V=W7Cefu;z;S)8zavPBiJb*^TTWGb1gi+RnusNT-c2Wkx;( zUci}(B;JH=7G~shFk#J*Z2IXT{6+?b_^VXkC z=0mE7(tPpD^!`p7hGpra+I8O@v9JvQoD2^%j$j74S?zhPU0!GKrl&~!T49q2**w9E{ z((!a0F)u8pClKIDtLYKadvPLWC*4@lReIDiVfM!=na2Qsf9*0O-LLf|)i^;-yS-{m zBV>w_**>Er3M|F4v-5jGbxKniS8qs)ZA28 zN80PmBpmfZ;_yPq{Unp8(n3eh{~-CH`jX9wQI}d?_+K>VVWpA@f6070;msO}o>AKp zkx~z6cPl6XTpK0~Qv(vygT-3W{5@shj{|S-Nuq0XYZ1;)PGmdI0QvRZA#?E9`fgsUH-GtTxqLwW9_p2t)(UW-1Xwu5 za%oSC|L$oo{qeu4KvT_N?`@#z;ZcF6hx!9e7f=f|dVV*-rgqtKHQ4kpK{gE?vG>Lv zD8k#cW3VarH59?7y3u_df(Y%G5LiS9z zHui^`?uX4#`?K!!*VuXW0a`V;PYapTm2gv{G0Je$8EGdSFCZF>MrrbEgqzMV&+PAI zgmBZ($v0);rk~3W^}jKjObX&(A3TCE`66$~?aYV?IO62n@nJ_~fq5Ti3N5gh3Nn0Fiz#ArnfMi8s!R2E3^0W( z@;(j+n6j6%DQb;X0!$Uy8Azk;G?~VMZZ-y^6b!A!r3IKmUJEeo@sy%SDeF%Eszztr zZ2_hdn7&VeKUs$;Ixh^o2r#X*hfASsPlML3npo}@=wfUm0!#_73c|@EVhS== z0!-O!d0C!;bIA*op~WjQr4>M74Qq860!#@ZJwuAA=_C))#;>G&fN67w0MjYXvdL)* z4o%23JrYhPlzx_AC!`U{Edxv$n-u4vn(xn9fT?LSuO}v$mJ_b` z4Jy)eRyf_yI+(I@IhQ)6j17U1k?afQv(9#k+ey8?YY82Dn*}NCTLxR&cjH4D) zx@g#z<<_uVD9A41ud0H569P}p z>o89UFvTnv-2!bbz!bFj0MnUjfaxKZ5~mhcXlYoXDZsQj6kysk^2?*sTptx+>e)l9 zO}4b?Bm|h+n}vIRh9FCJ5`BQF5e3s6PX>`096N$W3JvukriamBK8&*)LQJi97>BVa z#1sZa+Oe28TUZG(l{Vo>2JNsfrOp&$%IJ0funq2&MV?ehl?>2@1~faiTgVpn2EKZ@2xVzjSs!}Qotsnv@!=P zzBIr?e5t?=t{>Lw4;Wu6o#0Wt+xSvKB+$BE6~p3F2xZ!fFO`gj_)^Iz2hQe=uOPX$ zbW!U4@ukpxm+_^+4wHK41MlzmxTfkUyR__2rbY$s67 zysP-q1;2+s9ADbpReb4!k1vHAWy z6#~O}mQ>N=OBZ~6DOpm`JP4i>4fIxxFJ=Ac7N}#tBG(sRinCgdJwAr`(vFNDE?9vq zKE4#n$oI-`n+WfUFXg3?O`5^uk19+`o0EIb_|m4GL!!D=Lw|f}vlm}#2YTerlPS`R zIn2FFmC4HdZ0yME+NQm53Gt;8b(3r?1wkpOZBNx5-+2mT{ z$W|Ee+iF%8qc&Vv!HftkE$j8Qt8>mxMIy0+V2z5X3*yt}0j-|fl3}C`jFQYEp0}7S z)E2|Kp2o7Tqa1VI*rq?UR2P-*99p`dZGs>DTPA-9Mmlcm``ty9BIGOtXdcNy7ON4Z z7Tm}5Rr?}JJ%`>aqO^%;iDP0z5v3qGeJ>(PMJ%mv zVJqV$I7+@-j`9pel$u{vE9DUprGYm>=wrcCn54#Y(*@6kh|=Qw5YL5(QVZifZM8=? zY)mH&Wrvi%Zanuklc<&-cCrc;xx0u`oo83n?qEded}Y^3xaP(*1`#4k=jFDPcjtj4 zO6PZuDCK6Z5K(G(trQWZAOJJeB1%0jWGhHh5xPQPN{A@s3~T>aB1-S^m59#HG1QbaVd`muW1ium5QqCjlr}|JMU?VVN4^wM%B%4s zp{j~rm{V|2UqmSf5l2LnTI6ZOZj@%eUPP%^4NnFQ;U8!HPa>bSrp=GWj(Rx#>k0Rri7``LWL}K&=8rEl`-PE-)#n;mV zUz5JCX9vF4`o5kY_zIOOEngV;NjwwE5)oSYdU4<@5uuf@AF6zvWkkBf#P8Tn zaSZrEGga-B)$r{ph7IAo=2+;en7;JQ2_D~>MFXFQ#R*&2^$(wofcG@F1 z(G_Lc9=VD5EOQ@}+-p_T4UxO=>kzs7zE;T{B0j?j`)b5z%SU{oxNfJ;Q(~@m8B;*i zs9)|we=pKgWW39X`U71`D%nu7}_my(2q!Zn!_P zv!~g`wemkGvQy}J?;|_CkyDNAgdFQAt)`FcBveb0oz3CMPC^0v3`5l3mC=!%T2Hn{ zMs|Arb|y?E{;?Yq1-YBoaacY zQW?4F={FF$S&$jayb`%t=Julhh}J6=hG{(g#WNERS>)!uXxyV?O_7=mLIHNrsv4<@ z*3?Z_9wh^jnkAlkw5r+?Jz7<5iGo)BI4vVJhj3cp#Q*SQ4<`jq3!D@ z2#eGli4FWEmFUd=@s`n->=L=xU+v4U2A~I*mEBT0puEqO!Aca$_`Znm3ZdZ`ld@mZaZ-nOmBF?f`^=(jq zaA-~_H(n`a(U_G|MjFflyj+7~+dpg*0EVg1CA6FYbl+%9L_yJ*BQe0SLf~mB64OhJ zTk0FKMPXLPD-U0Fghdo)A#Or}Mr%=+1xnEeK# zFwIuEe-vh`|7Xv38HK6eNaJVEI`9;>$SJmgd5nz0BrebY>{&BJVdDCk&gX`sFq=+M zn1u4Z2T_>%RT)vk|Ei4Ms=q43+N3WEQ@<+n9z;XR1LETTq36eefMUC*lr zdyI_2bRddj7%Dk&`6$d;MPasx!faNeFi~StWmg_33KJ#5YCu{8=8S97{k+gO7=@YmD9pLhQJAiGO7|x%3Ui@|!kkeQrf%9*6z2XAg}J|o z!sMioOF56&ABDO9U)Lzix{tzSx%gg2Vg94|uOg0);7x57)b#kNjH4qky)TYV>1rIE z($zRRrK@ptq%$8!$Gr(1UBjnaLKQ7H>k~e7ZpJ4h?uw(cVl4Ii*`cu4`cn``r_m?! zJn)xIUzA}bj!vlXade;I`y-2^)6ZS)K90`cSt~bMe;ggxG59#Ta23RR6-QU}HyB4} z{V9%au)p^+j!wkCZyX(kbi}92e%pk2n&(MUQJOQ+1CSUP3iJC;r*`eW(HLA=sYv2^z9 zC>f`jD1cTwqJr*mf7QEU={RcVW9>$K9~})Q-y_2|<)K(QRa7ipui?R1y4{0!tkMd- z41w48{XxL1Rve6_yEk|~R!)OcEM0GqfX6kH4b|v2^^RjaLBoilw9K=vX>FASM7Hsb43OFT~P~kQ#ipAw15` zIacH8^s7Bpov}kaT^K+R)DTZs_Eg5yA)mXBsS9_e+#{wgl=!h?>blX9$VKbA=10Gw z7`4XvhYL^Y^jp2){z>0J(~&=XYK0DhcAObY8r> z*lL~M>604JMiLc$X##w8Ak1Gj(faMUtXNl@teO?M?=J)4?VvxqL}V6GApATf-wgyp ze7m|vpENT(0USES#FNV~?COY`P6ZO^t`HUn40NtDs@9In_BiAwQ=nXr)-oGP6p~cOegcY&TYkc^r4s-n;sSk=HQfm z4~s5-EIO<4=Ttta)EUks>32N>^u6Kt#w`?tD>rxx7nt}_57W%|rAt3rtq&B4$o+*P zM%V7a!Aic$!D(*VSr}cT)%xRFJQl7(sdHNQ4&$=%b|*T~ft*(;Y|s3M$GQ*M)jkKf z9D<`3x^#ixK#eX&s&IxV79F(;NghHe59rPQna{7Wp;jru&!nQR$_z(!PDBw(NJZ_Y z(kg+Lx-o4`t*a>6P-ZNjWNEGTr3CHauWj-5YI)lirQrUupngM_@)lpDqrR3R(u+fC zMYXx|ZJ~TSSTse2`fUv0Y>Z!B%$>E%Eg)%(7be!#Fm#{?E$abXdI(B>!a`W!z+o*v z%=P?ux9j3Iu7Jt!2qHL15V1#P$G({k5Q+#==-4n42nBULjLXQ=kx4Eqv7dLq=b)Tc z%4dH;OVs*5mTIT(BtIX2rEHpaXwA zbnBUbWHAxb@(&^+$IxGd5osJV8tN%*$1WjGzzje$1Q~OI%P~5d;f%KaIkQpxdN7Go zEJpcNX8+5Hbd!kP?UCqGoJ6MzbSjT4WW;H8W!9^i%)P^zgaq_ghiZPkAv;v3 ze}(=rG1?fc&2W-l$CK(~vgiXzRLz7(#(ghV69I|RItakM-L~ z+(zKVI2QU7FU8L(gN~X^0LfLj_a?sXY6M{^RElC7!b10`Z!PvM=*Q_fu}3N=!V}fm zg<9$@d=D{m`aT!Fi%uf=hhnFP0{DzDcmM)bU@QQnl?^=Dc7k2{g1}q+=Az1^Vu!{- zjszJu_n0D+#tP)p9Os&ko+em_Rx zG8`iStsO>p7H=_3Q*7>j?QNwQ)le4vJXosMmzAESqH#G&5>4rU#1hleBEo8;`m!p< z-ZdezR@?##T3_IIC+pAP1jFtFf)kSU7rAXRPV4m5 z>EKZ|q-flBNaSnjIpNvx<NgM5bA#{lD0o#H>&0Tx7{oEhsFJM_wKFc`~i}`706pO-g0NTolU$fU7!&{-X3U{EsN;Cmv!jc!sc@CjtsZ^0uhCS+r9w3ja37KR={ z>*>UfHOTL#qA~nCT1N%v%w3S*O}a-8V5P@y^05fF+U-T(HvaBhJ^wv_>E3^JNaYh&1_WWyt`cQ# z^;e0aT3SN4tZCG*)fHBcX!LMHG84%TZhhU7!6Vw9RnWNk|V;1nhrOVbh zO?U#6HnjX^c6sKu=BRLDqE<0^F)vSYK9)OFx(x0|t#AFb4r}{IqI& zip^o8kMWqls81;7ExF8kw!Q^>(9|3BGYa&ezSrjBCmGHl^DfI(Lp8M6A@7P7%OobxE7fCAxUU|sIK%S3+ z>>HUN7U&s_j8-g1#v;m(Y`-5FSTkH?+|nQmAJJ&@;KwXk#rlz1=0WBcMy8b&|Jcat zvhusaJD9SOQBny%GL1?XvkQ#`$SdW8^3o}94-v>F*FtjH8i4A$DPh0epKJ|h8`&dX zZ|N7_$cnNf>7uEEBUq@9zzfhvMQqy~aum1Daix%LoF6e&58Paw+|zcTLxW^twU#bg zyXg^~D5AfmzW(aj(e$A8ca;8?e19r3*xy0x?P#SpSPAtc-3oWAw?w^xL~#@Mh_Wd! z3+^zCn|=V+KB)fnm*F?VLHP7BGMelBwuhCSaNJ9z&GrNM^^kN?=rL`Ha2DI`(=~HQ z3F4%}wAgekfnXeP%D>#8hKW}jY6+ytOvQJ&#ytZ3YqpM>)`c4W#Ix^1Q)?5ij%+MO zhNwq=UzPhYG!&rD>o?I`Fu@{MGMJX?0?2;HS|q}g&tI;Z>tL8@p}{>C=)u!G(1VZj zAg1y-502wMs>gbI!c8K0Vn3QX~viUA_-lc3_VYbWX1T*91pnHjIyrH{+dCyvQ$Q|-{w);tY zc*gC#neF&}G@eED5s0p~+!FhJ1nB162Z}}3uNoWdDOmyG_(YW(rWLnv)-Hzmx@LB+ ziUYV2ncOcrT(V8(G}SGjSk)!Q#v@2G_Yrs=dA!o0vd=nQSzz74{3^dke6hzp?#sh2 zYy9TKlRia6_YOkH73k7c$B;b)j*IEE7(~D|xzQ*WuX=q_s=KnGZlUgi96x4F>@)E# zUFlp;*VVQ1{06f>z=<+l&!lH%SRjZ?z4I&w-O-_2=P2<=uf(D|GE@Ti&Csvk_hLFf zS@oOhruqlWcV<%EWZuHC)IvZ@h<|a^mJD(JOS5(Z`QW|Z1NVn>Vy&n-cmw!1kLLel zH~)9z{P+Bg;JAzK^-6>6|128@@B2V-9U^#5}xeDGuyHZie;nuk ziC+&j3>~?y{sR~Hy!&tH0`&jlNd6C}^Is}1TFwV9IRCxlayx$ceRezkubrJ3 z(s6J*hURwXcI4OD@#o#c?O2oH_n>aa1AT7AAIt5y}8-qR!918`op)!Ca3^ABvYGOvh`EAI|VjLieDq$Kwy|dIU@Er{DF6tnTc3 zO@%w1^S#VR*uy<3cLAfb= z4e#WpEC4q*nopH{le?p(aW1R zuf1_~JG!>D{rKk1FK#Uzcb6_XQlBE>9)4J_k6*pH(fjZ*64BQ7@}`Qo4_|VpRP5D_ zn^)JZPLWd(*w``eT)%)U#4F+24@*D09jpp`idxXvn~Mbz^z!m21n}%dcE# zv@4s@$G*6_wer!`8_Szt`s6L2dg+ag4WHb)y1cQH-?)}vyY}gu(dHXBqF1iIzH)W_ zm93i_Z?3$uytTEmDKs(4YZouBymfWumL)%`F-DI+@p$F`XI8(ka^ngA@7CtkCq8y_ zbLELo3HOh`zMG=AJpQ+ycy8s!%I5OMi_6=at8e*#ch|zpE1TEf*j^Td4^r*vo7=l= z-$VF7g`atIdE<@c?VFohyQ}tZtiqT#H+NS!-`ZOJ!i}#R@+^jH?~(dp7V+}`RYn! zI^y{o;Ddy({>tiWuPkqD-2CFoHLo&WUA|@+S0i*$A)GJM`qj-#mp%|JE&U9g=Gz~m z*VT=cPp@8Iky(wgb2HZ~~4uEtO z`WRl9+|rT>s?ZPn_)N9ex6y#)NF*>maA@+%v!Zxv>W zz$+~uecR-lZ#-Ilp;rz@TK-foNj3ZGON!0xYpB@k&7^Vj>d4-TZ=upnbcm)>pcPMO z_-o<0`jiCOV|XpRM-02yI}~_jyK4GvJ?^VSp(YS|l!xr?*H^b5`>-e6hpjW8;^42XTHMrQVymf11 z_3G-jd+m)Iac5?w4AU{jCPdBv@6xtljOzT{q`v#r~22t!NDHvkoz@Fo+B zv*~4f_06sYLxZqBh6?&V1Va!} z6%DW1^MNc9wJrC>*H>gs32c7uNyG_^40A(uu{I{%v5pPH+|#wYCJ%58mBK-4n$dGv4vzWMB z&8+smQgc-KL75n?_0k7Uz4FS{x87QQb@j~;eBj(GuWY@(di|AcCNa#c*KR(eP1j$Z zjiTR}iJ}+gqNu<8dFo;;Z}nGD`9GrkclbM0{?@D5iNoc8G_w2$UwP%$O(u$!&1wy` zf0DnBzxwZo)J=HS_*>`i7=Oq4JHg)ue<%5?|Gr0gm22|1TCPR9&EF1x``|qdZnjoO z4bgD&P1;l6`tKCDiFr=(caJjoZz1Gy3@jqRGWwbA5B65Dq_0uGK^gLUW-9p_pCe!O X`BoBxPgdutxhVhMQmgj)z48AGCJ20E