ml-airpls - npm Package Compare versions

Comparing version 0.0.2 to 0.2.0

History.md

@@ -0,1 +1,5 @@
+# [0.2.0](https://github.com/mljs/airpls/compare/v0.0.2...v0.2.0) (2020-02-05)
+
+
+
 <a name="0.0.2"></a>
@@ -2,0 +6,0 @@ ## 0.0.2 (2018-11-03)

lib/index.js

@@ -5,6 +5,344 @@ 'use strict';
 
-var Cholesky = _interopDefault(require('cholesky-solve'));
 var cuthillMckee = _interopDefault(require('cuthill-mckee'));
 
-function airPLS(yData, options = {}) {
+function ldlSymbolic(
+  n /* A and L are n-by-n, where n >= 0 */,
+  Ap /* input of size n + 1, not modified */,
+  Ai /* input of size nz=Ap[n], not modified */,
+  Lp /* output of size n + 1, not defined on input */,
+  Parent /* output of size n, not defined on input */,
+  Lnz /* output of size n, not defined on input */,
+  Flag /* workspace of size n, not defn. on input or output */,
+) {
+  let i, k, p, kk, p2;
+
+  for (k = 0; k < n; k++) {
+    /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
+    Parent[k] = -1; /* parent of k is not yet known */
+    Flag[k] = k; /* mark node k as visited */
+    Lnz[k] = 0; /* count of nonzeros in column k of L */
+    kk = k; /* kth original, or permuted, column */
+    p2 = Ap[kk + 1];
+    for (p = Ap[kk]; p < p2; p++) {
+      /* A (i,k) is nonzero (original or permuted A) */
+      i = Ai[p];
+
+      if (i < k) {
+        /* follow path from i to root of etree, stop at flagged node */
+        for (; Flag[i] !== k; i = Parent[i]) {
+          /* find parent of i if not yet determined */
+          if (Parent[i] === -1) Parent[i] = k;
+          Lnz[i]++; /* L (k,i) is nonzero */
+          Flag[i] = k; /* mark i as visited */
+        }
+      }
+    }
+  }
+  /* construct Lp index array from Lnz column counts */
+  Lp[0] = 0;
+  for (k = 0; k < n; k++) {
+    Lp[k + 1] = Lp[k] + Lnz[k];
+  }
+}
+
+function ldlNumeric(
+  n /* A and L are n-by-n, where n >= 0 */,
+  Ap /* input of size n+1, not modified */,
+  Ai /* input of size nz=Ap[n], not modified */,
+  Ax /* input of size nz=Ap[n], not modified */,
+  Lp /* input of size n+1, not modified */,
+  Parent /* input of size n, not modified */,
+  Lnz /* output of size n, not defn. on input */,
+  Li /* output of size lnz=Lp[n], not defined on input */,
+  Lx /* output of size lnz=Lp[n], not defined on input */,
+  D /* output of size n, not defined on input */,
+  Y /* workspace of size n, not defn. on input or output */,
+  Pattern /* workspace of size n, not defn. on input or output */,
+  Flag /* workspace of size n, not defn. on input or output */,
+) {
+  let yi, lKi;
+  let i, k, p, kk, p2, len, top;
+  for (k = 0; k < n; k++) {
+    /* compute nonzero Pattern of kth row of L, in topological order */
+    Y[k] = 0.0; /* Y(0:k) is now all zero */
+    top = n; /* stack for pattern is empty */
+    Flag[k] = k; /* mark node k as visited */
+    Lnz[k] = 0; /* count of nonzeros in column k of L */
+    kk = k; /* kth original, or permuted, column */
+    p2 = Ap[kk + 1];
+    for (p = Ap[kk]; p < p2; p++) {
+      i = Ai[p]; /* get A(i,k) */
+      if (i <= k) {
+        Y[i] += Ax[p]; /* scatter A(i,k) into Y (sum duplicates) */
+        for (len = 0; Flag[i] !== k; i = Parent[i]) {
+          Pattern[len++] = i; /* L(k,i) is nonzero */
+          Flag[i] = k; /* mark i as visited */
+        }
+        while (len > 0) Pattern[--top] = Pattern[--len];
+      }
+    }
+    /* compute numerical values kth row of L (a sparse triangular solve) */
+    D[k] = Y[k]; /* get D(k,k) and clear Y(k) */
+    Y[k] = 0.0;
+    for (; top < n; top++) {
+      i = Pattern[top]; /* Pattern[top:n-1] is pattern of L(:,k) */
+      yi = Y[i]; /* get and clear Y(i) */
+      Y[i] = 0.0;
+      p2 = Lp[i] + Lnz[i];
+      for (p = Lp[i]; p < p2; p++) {
+        Y[Li[p]] -= Lx[p] * yi;
+      }
+      lKi = yi / D[i]; /* the nonzero entry L(k,i) */
+      D[k] -= lKi * yi;
+      Li[p] = k; /* store L(k,i) in column form of L */
+      Lx[p] = lKi;
+      Lnz[i]++; /* increment count of nonzeros in col i */
+    }
+
+    if (D[k] === 0.0) return k; /* failure, D(k,k) is zero */
+  }
+
+  return n; /* success, diagonal of D is all nonzero */
+}
+
+function ldlLsolve(
+  n /* L is n-by-n, where n >= 0 */,
+  X /* size n. right-hand-side on input, soln. on output */,
+  Lp /* input of size n+1, not modified */,
+  Li /* input of size lnz=Lp[n], not modified */,
+  Lx /* input of size lnz=Lp[n], not modified */,
+) {
+  let j, p, p2;
+  for (j = 0; j < n; j++) {
+    p2 = Lp[j + 1];
+    for (p = Lp[j]; p < p2; p++) {
+      X[Li[p]] -= Lx[p] * X[j];
+    }
+  }
+}
+function ldlDsolve(
+  n /* D is n-by-n, where n >= 0 */,
+  X /* size n. right-hand-side on input, soln. on output */,
+  D /* input of size n, not modified */,
+) {
+  let j;
+  for (j = 0; j < n; j++) {
+    X[j] /= D[j];
+  }
+}
+function ldlLTsolve(
+  n /* L is n-by-n, where n >= 0 */,
+  X /* size n. right-hand-side on input, soln. on output */,
+  Lp /* input of size n+1, not modified */,
+  Li /* input of size lnz=Lp[n], not modified */,
+  Lx /* input of size lnz=Lp[n], not modified */,
+) {
+  let j, p, p2;
+  for (j = n - 1; j >= 0; j--) {
+    p2 = Lp[j + 1];
+    for (p = Lp[j]; p < p2; p++) {
+      X[j] -= Lx[p] * X[Li[p]];
+    }
+  }
+}
+
+function ldlPerm(
+  n /* size of X, B, and P */,
+  X /* output of size n. */,
+  B /* input of size n. */,
+  P /* input permutation array of size n. */,
+) {
+  let j;
+  for (j = 0; j < n; j++) {
+    X[j] = B[P[j]];
+  }
+}
+
+function ldlPermt(
+  n /* size of X, B, and P */,
+  X /* output of size n. */,
+  B /* input of size n. */,
+  P /* input permutation array of size n. */,
+) {
+  let j;
+  for (j = 0; j < n; j++) {
+    X[P[j]] = B[j];
+  }
+}
+
+function prepare(M, n, P) {
+  // if a permutation was specified, apply it.
+  if (P) {
+    let Pinv = new Array(n);
+
+    for (let k = 0; k < n; k++) {
+      Pinv[P[k]] = k;
+    }
+
+    let Mt = []; // scratch memory
+    // Apply permutation. We make M into P*M*P^T
+    for (let a = 0; a < M.length; ++a) {
+      let ar = Pinv[M[a][0]];
+      let ac = Pinv[M[a][1]];
+
+      // we only store the upper-diagonal elements(since we assume matrix is symmetric, we only need to store these)
+      // if permuted element is below diagonal, we simply transpose it.
+      if (ac < ar) {
+        let t = ac;
+        ac = ar;
+        ar = t;
+      }
+
+      Mt[a] = [];
+      Mt[a][0] = ar;
+      Mt[a][1] = ac;
+      Mt[a][2] = M[a][2];
+    }
+
+    M = Mt; // copy scratch memory.
+  } else {
+    // if P argument is null, we just use an identity permutation.
+    P = [];
+    for (let i = 0; i < n; ++i) {
+      P[i] = i;
+    }
+  }
+
+  // The sparse matrix we are decomposing is A.
+  // Now we shall create A from M.
+  let Ap = new Array(n + 1);
+  let Ai = new Array(M.length);
+  let Ax = new Array(M.length);
+
+  // count number of non-zero elements in columns.
+  let LNZ = [];
+  for (let i = 0; i < n; ++i) {
+    LNZ[i] = 0;
+  }
+  for (let a = 0; a < M.length; ++a) {
+    LNZ[M[a][1]]++;
+  }
+
+  Ap[0] = 0;
+  for (let i = 0; i < n; ++i) {
+    Ap[i + 1] = Ap[i] + LNZ[i];
+  }
+
+  let coloffset = [];
+  for (let a = 0; a < n; ++a) {
+    coloffset[a] = 0;
+  }
+
+  // go through all elements in M, and add them to sparse matrix A.
+  for (let i = 0; i < M.length; ++i) {
+    let e = M[i];
+    let col = e[1];
+
+    let adr = Ap[col] + coloffset[col];
+    Ai[adr] = e[0];
+    Ax[adr] = e[2];
+
+    coloffset[col]++;
+  }
+
+  let D = new Array(n);
+  let Y = new Array(n);
+  let Lp = new Array(n + 1);
+  let Parent = new Array(n);
+  let Lnz = new Array(n);
+  let Flag = new Array(n);
+  let Pattern = new Array(n);
+  let bp1 = new Array(n);
+  let x = new Array(n);
+  let d;
+
+  ldlSymbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag);
+
+  let Lx = new Array(Lp[n]);
+  let Li = new Array(Lp[n]);
+
+  d = ldlNumeric(n, Ap, Ai, Ax, Lp, Parent, Lnz, Li, Lx, D, Y, Pattern, Flag);
+
+  if (d === n) {
+    return function(b) {
+      ldlPerm(n, bp1, b, P);
+      ldlLsolve(n, bp1, Lp, Li, Lx);
+      ldlDsolve(n, bp1, D);
+      ldlLTsolve(n, bp1, Lp, Li, Lx);
+      ldlPermt(n, x, bp1, P);
+
+      return x;
+    };
+  } else {
+    return null;
+  }
+}
+
+const getClosestNumber = (array = [], goal = 0) => {
+  const closest = array.reduce((prev, curr) => {
+    return Math.abs(curr - goal) < Math.abs(prev - goal) ? curr : prev;
+  });
+  return closest;
+};
+
+const getCloseIndex = (array = [], goal = 0) => {
+  const closest = getClosestNumber(array, goal);
+  return array.indexOf(closest);
+};
+
+const updateSystem = (matrix, y, weights) => {
+  let nbPoints = y.length;
+  let l = nbPoints - 1;
+  let newMatrix = new Array(matrix.length);
+  let newVector = new Float64Array(nbPoints);
+  for (let i = 0; i < l; i++) {
+    let w = weights[i];
+    let diag = i * 2;
+    let next = diag + 1;
+    newMatrix[diag] = matrix[diag].slice();
+    newMatrix[next] = matrix[next].slice();
+    if (w === 0) {
+      newVector[i] = 0;
+    } else {
+      newVector[i] = y[i] * w;
+      newMatrix[diag][2] += w;
+    }
+  }
+  newVector[l] = y[l] * weights[l];
+  newMatrix[l * 2] = matrix[l * 2].slice();
+  newMatrix[l * 2][2] += weights[l];
+
+  return [newMatrix, newVector];
+};
+
+const getDeltaMatrix = (nbPoints, lambda) => {
+  let matrix = [];
+  let last = nbPoints - 1;
+  for (let i = 0; i < last; i++) {
+    matrix.push([i, i, lambda * 2]);
+    matrix.push([i + 1, i, -1 * lambda]);
+  }
+  matrix[0][2] = lambda;
+  matrix.push([last, last, lambda]);
+  return {
+    lowerTriangularNonZeros: matrix,
+    permutationEncodedArray: cuthillMckee(matrix, nbPoints),
+  };
+};
+
+/**
+ * Fit the baseline drift by iteratively changing weights of sum square error between the fitted baseline and original signals,
+ * for further information about the parameters you can get the [paper of airPLS](https://github.com/zmzhang/airPLS/blob/master/airPLS_manuscript.pdf)
+ * @param {Array<number>} x - x axis data useful when control points or zones are submitted
+ * @param {Array<number>} y - Original data
+ * @param {object} [options={}] - Options object
+ * @param {number} [options.maxIterations = 100] - Maximal number of iterations if the method does not reach the stop criterion
+ * @param {number} [options.factorCriterion = 0.001] - Factor of the sum of absolute value of original data, to compute stop criterion
+ * @param {Array<number>} [options.weights = [1,1,...]] - Initial weights vector, default each point has the same weight
+ * @param {number} [options.lambda = 100] - Factor of weights matrix in -> [I + lambda D'D]z = x
+ * @param {Array<number>} [options.controlPoints = []] - Array of x axis values to force that baseline cross those points.
+ * @param {Array<number>} [options.baseLineZones = []] - Array of x axis values (as from - to), to force that baseline cross those zones.
+ * @returns {{corrected: Array<number>, error: number, iteration: number, baseline: Array<number>}}
+ */
+function airPLS(x, y, options = {}) {
   let {
@@ -14,28 +352,54 @@ maxIterations = 100,
     factorCriterion = 0.001,
-    weights
+    weights = new Array(y.length).fill(1),
+    controlPoints = [],
+    baseLineZones = [],
   } = options;
 
-  var nbPoints = yData.length;
-  var stopCriterion = factorCriterion * yData.reduce((sum, e) => Math.abs(e) + sum, 0);
-
-  if (!weights) {
-    weights = new Array(nbPoints).fill(1);
+  if (controlPoints.length > 0) {
+    controlPoints.forEach((e, i, arr) => (arr[i] = getCloseIndex(x, e)));
   }
+  if (baseLineZones.length > 0) {
+    baseLineZones.forEach((range) => {
+      let indexFrom = getCloseIndex(x, range.from);
+      let indexTo = getCloseIndex(x, range.to);
+      if (indexFrom > indexTo) [indexFrom, indexTo] = [indexTo, indexFrom];
+      for (let i = indexFrom; i < indexTo; i++) {
+        controlPoints.push(i);
+      }
+    });
+  }
 
-  var { lowerTriangularNonZeros, permutationEncodedArray } = getDeltaMatrix(nbPoints, lambda);
+  let baseline, iteration;
+  let nbPoints = y.length;
+  let l = nbPoints - 1;
+  let sumNegDifferences = Number.MAX_SAFE_INTEGER;
+  let stopCriterion =
+    factorCriterion * y.reduce((sum, e) => Math.abs(e) + sum, 0);
 
-  var sumNegDifferences = Number.MAX_SAFE_INTEGER;
-  for (var iteration = 0; (iteration < maxIterations && Math.abs(sumNegDifferences) > stopCriterion); iteration++) {
-    let [leftHandSide, rightHandSide] = updateSystem(lowerTriangularNonZeros, yData, weights);
+  let { lowerTriangularNonZeros, permutationEncodedArray } = getDeltaMatrix(
+    nbPoints,
+    lambda,
+  );
 
-    let cho = Cholesky.prepare(leftHandSide, nbPoints, permutationEncodedArray);
+  for (
+    iteration = 0;
+    iteration < maxIterations && Math.abs(sumNegDifferences) > stopCriterion;
+    iteration++
+  ) {
+    let [leftHandSide, rightHandSide] = updateSystem(
+      lowerTriangularNonZeros,
+      y,
+      weights,
+    );
 
-    var baseline = cho(rightHandSide);
+    let cho = prepare(leftHandSide, nbPoints, permutationEncodedArray);
 
+    baseline = cho(rightHandSide);
+
     sumNegDifferences = 0;
 
-    let difference = yData.map(calculateError);
+    let difference = y.map(calculateError);
 
     let maxNegativeDiff = -1 * Number.MAX_SAFE_INTEGER;
-    for (var i = 1, l = nbPoints - 1; i < l; i++) {
+    for (let i = 1; i < l; i++) {
       let diff = difference[i];
@@ -45,3 +409,3 @@ if (diff >= 0) {
       } else {
-        weights[i] = Math.exp(iteration * diff / sumNegDifferences);
+        weights[i] = Math.exp((iteration * diff) / sumNegDifferences);
         if (maxNegativeDiff < diff) maxNegativeDiff = diff;
@@ -51,12 +415,13 @@ }
 
-    let value = Math.exp(iteration * maxNegativeDiff / sumNegDifferences);
+    let value = Math.exp((iteration * maxNegativeDiff) / sumNegDifferences);
     weights[0] = value;
     weights[l] = value;
+    controlPoints.forEach((i) => (weights[i] = value));
   }
 
   return {
-    corrected: yData.map((e, i) => e - baseline[i]),
+    corrected: y.map((e, i) => e - baseline[i]),
     baseline,
     iteration,
-    error: sumNegDifferences
+    error: sumNegDifferences,
   };
@@ -71,38 +436,2 @@
 
-
-function getDeltaMatrix(nbPoints, lambda) {
-  var matrix = [];
-  for (var i = 0, last = nbPoints - 1; i < last; i++) {
-    matrix.push([i, i, lambda * 2]);
-    matrix.push([i + 1, i, -1 * lambda]);
-  }
-  matrix[0][2] = lambda;
-  matrix.push([last, last, lambda]);
-  return { lowerTriangularNonZeros: matrix, permutationEncodedArray: cuthillMckee(matrix, nbPoints) };
-}
-
-function updateSystem(matrix, yData, weights) {
-  let nbPoints = yData.length;
-  var newMatrix = new Array(matrix.length);
-  var newVector = new Float64Array(nbPoints);
-  for (var i = 0, l = nbPoints - 1; i < l; i++) {
-    let w = weights[i];
-    let diag = i * 2;
-    let next = diag + 1;
-    newMatrix[diag] = matrix[diag].slice();
-    newMatrix[next] = matrix[next].slice();
-    if (w === 0) {
-      newVector[i] = 0;
-    } else {
-      newVector[i] = yData[i] * w;
-      newMatrix[diag][2] += w;
-    }
-  }
-  newVector[l] = yData[l] * weights[l];
-  newMatrix[l * 2] = matrix[l * 2].slice();
-  newMatrix[l * 2][2] += weights[l];
-
-  return [newMatrix, newVector];
-}
-
 module.exports = airPLS;

package.json

 {
   "name": "ml-airpls",
-  "version": "0.0.2",
+  "version": "0.2.0",
   "description": "Baseline correction using adaptive iteratively reweighted penalized least",
-  "main": "src/index.js",
+  "main": "lib/index.js",
   "module": "src/index.js",
   "files": [
     "lib",
-    "runkit.js",
     "src"
   ],
   "scripts": {
+    "compile": "rollup -c",
+    "build": "cheminfo-build",
+    "prepublishOnly": "npm run compile",
     "eslint": "eslint src",
     "eslint-fix": "npm run eslint -- --fix",
-    "prepublish": "rollup -c",
-    "test": "run-s testonly eslint",
+    "test": "run-s test-only eslint",
     "test-travis": "eslint src && jest --coverage && codecov",
-    "testonly": "jest"
+    "test-only": "jest"
   },
   "repository": {
     "type": "git",
-    "url": "https://github.com/mljs/airPLS.git"
+    "url": "https://github.com/mljs/airpls.git"
   },
@@ -28,9 +29,8 @@ "keywords": [
   ],
-  "author": "Jose Alejandro Bolanos Arroyave",
+  "author": "Jose Alejandro Bolanos Arroyave (jobo322)",
   "license": "MIT",
   "bugs": {
-    "url": "https://github.com/mljs/airPLS/issues"
+    "url": "https://github.com/mljs/airpls/issues"
   },
-  "runkitExampleFilename": "./runkit.js",
-  "homepage": "https://github.com/mljs/airPLS#readme",
+  "homepage": "https://github.com/mljs/airpls#readme",
   "jest": {
@@ -40,18 +40,20 @@ "testEnvironment": "node"
   "devDependencies": {
-    "babel-plugin-transform-es2015-modules-commonjs": "^6.26.0",
-    "codecov": "^2.3.0",
-    "eslint": "^4.11.0",
-    "eslint-config-cheminfo": "^1.18.0",
-    "eslint-plugin-import": "^2.14.0",
-    "eslint-plugin-jest": "^21.26.2",
-    "eslint-plugin-no-only-tests": "^2.0.0",
-    "jest": "^21.2.1",
-    "jest-matcher-deep-close-to": "^1.0.2",
-    "npm-run-all": "^4.1.1",
-    "rollup": "^0.42.0"
+    "@babel/plugin-transform-modules-commonjs": "^7.8.3",
+    "cheminfo-build": "^1.0.5",
+    "codecov": "^3.6.4",
+    "eslint": "^6.8.0",
+    "eslint-config-cheminfo": "^2.0.4",
+    "eslint-plugin-import": "^2.20.1",
+    "eslint-plugin-jest": "^23.6.0",
+    "eslint-plugin-prettier": "^3.1.2",
+    "esm": "^3.2.25",
+    "jest": "^25.1.0",
+    "jest-matcher-deep-close-to": "^1.3.0",
+    "npm-run-all": "^4.1.5",
+    "prettier": "^1.19.1",
+    "rollup": "^1.31.0"
   },
   "dependencies": {
-    "cholesky-solve": "^0.2.1",
     "cuthill-mckee": "^1.0.0"
   }
 }

README.md

@@ -1,2 +0,2 @@
-# airPLS
+# airpls
 
@@ -11,7 +11,9 @@ [![NPM version][npm-image]][npm-url]
 
+It is an javascript implementation of [airpls](https://github.com/zmzhang/airPLS/blob/master/airPLS_manuscript.pdf) using cholesky decomposition and reverse Cuthill-Mckee method for reducing the bandwidth of sparse linear systems, obtaining a fast baseline fitter. 
+
 ## Installation
 
-`$ npm install ml-airPLS`
+`$ npm install ml-airpls`
 
-## [API Documentation](https://mljs.github.io/airPLS/)
+## [API Documentation](https://mljs.github.io/airpls/)
 
@@ -21,7 +23,10 @@ ## Example
 ```js
-const airPls = require('ml-airPLS');
+const airpls = require('ml-airpls');
+
+let y = [1, 1, 1, 1, 3, 6, 3, 1, 1, 1];
+let x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
+
+var {baseline, corrected, iteration, error} = airpls(x, y);
 ```
 
-Or test it in [Runkit](https://runkit.com/npm/ml-airPLS)
-
 ## License
@@ -31,11 +36,11 @@
 
-[npm-image]: https://img.shields.io/npm/v/ml-airPLS.svg?style=flat-square
-[npm-url]: https://www.npmjs.com/package/ml-airPLS
-[travis-image]: https://img.shields.io/travis/mljs/airPLS/master.svg?style=flat-square
-[travis-url]: https://travis-ci.org/mljs/airPLS
-[codecov-image]: https://img.shields.io/codecov/c/github/mljs/airPLS.svg?style=flat-square
-[codecov-url]: https://codecov.io/gh/mljs/airPLS
-[david-image]: https://img.shields.io/david/mljs/airPLS.svg?style=flat-square
-[david-url]: https://david-dm.org/mljs/airPLS
-[download-image]: https://img.shields.io/npm/dm/ml-airPLS.svg?style=flat-square
-[download-url]: https://www.npmjs.com/package/ml-airPLS
+[npm-image]: https://img.shields.io/npm/v/ml-airpls.svg?style=flat-square
+[npm-url]: https://www.npmjs.com/package/ml-airpls
+[travis-image]: https://img.shields.io/travis/mljs/airpls/master.svg?style=flat-square
+[travis-url]: https://travis-ci.org/mljs/airpls
+[codecov-image]: https://img.shields.io/codecov/c/github/mljs/airpls.svg?style=flat-square
+[codecov-url]: https://codecov.io/gh/mljs/airpls
+[david-image]: https://img.shields.io/david/mljs/airpls.svg?style=flat-square
+[david-url]: https://david-dm.org/mljs/airpls
+[download-image]: https://img.shields.io/npm/dm/ml-airpls.svg?style=flat-square
+[download-url]: https://www.npmjs.com/package/ml-airpls

src/__tests__/test.js

@@ -7,7 +7,8 @@ import { toBeDeepCloseTo } from 'jest-matcher-deep-close-to';
 
-var vector = [1, 1, 1, 1, 3, 6, 3, 1, 1, 1];
+let y = [1, 1, 1, 1, 3, 6, 3, 1, 1, 1];
+let x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
 
 describe('airPLS test', () => {
   it('Small vector to find baseline', () => {
-    var result = airPLS(vector);
+    let result = airPLS(x, y);
     expect(result.baseline).toBeDeepCloseTo([1, 1, 1, 1, 1, 1, 1, 1, 1, 1], 2);
@@ -14,0 +15,0 @@ expect(result.corrected).toBeDeepCloseTo([0, 0, 0, 0, 2, 5, 2, 0, 0, 0], 2);

src/index.js

@@ -1,5 +0,19 @@
-import Cholesky from 'cholesky-solve';
-import cuthillMckee from 'cuthill-mckee';
+import Cholesky from './choleskySolver';
+import { updateSystem, getDeltaMatrix, getCloseIndex } from './utils';
 
-function airPLS(yData, options = {}) {
+/**
+ * Fit the baseline drift by iteratively changing weights of sum square error between the fitted baseline and original signals,
+ * for further information about the parameters you can get the [paper of airPLS](https://github.com/zmzhang/airPLS/blob/master/airPLS_manuscript.pdf)
+ * @param {Array<number>} x - x axis data useful when control points or zones are submitted
+ * @param {Array<number>} y - Original data
+ * @param {object} [options={}] - Options object
+ * @param {number} [options.maxIterations = 100] - Maximal number of iterations if the method does not reach the stop criterion
+ * @param {number} [options.factorCriterion = 0.001] - Factor of the sum of absolute value of original data, to compute stop criterion
+ * @param {Array<number>} [options.weights = [1,1,...]] - Initial weights vector, default each point has the same weight
+ * @param {number} [options.lambda = 100] - Factor of weights matrix in -> [I + lambda D'D]z = x
+ * @param {Array<number>} [options.controlPoints = []] - Array of x axis values to force that baseline cross those points.
+ * @param {Array<number>} [options.baseLineZones = []] - Array of x axis values (as from - to), to force that baseline cross those zones.
+ * @returns {{corrected: Array<number>, error: number, iteration: number, baseline: Array<number>}}
+ */
+function airPLS(x, y, options = {}) {
   let {
@@ -9,28 +23,54 @@ maxIterations = 100,
     factorCriterion = 0.001,
-    weights
+    weights = new Array(y.length).fill(1),
+    controlPoints = [],
+    baseLineZones = [],
   } = options;
 
-  var nbPoints = yData.length;
-  var stopCriterion = factorCriterion * yData.reduce((sum, e) => Math.abs(e) + sum, 0);
-
-  if (!weights) {
-    weights = new Array(nbPoints).fill(1);
+  if (controlPoints.length > 0) {
+    controlPoints.forEach((e, i, arr) => (arr[i] = getCloseIndex(x, e)));
   }
+  if (baseLineZones.length > 0) {
+    baseLineZones.forEach((range) => {
+      let indexFrom = getCloseIndex(x, range.from);
+      let indexTo = getCloseIndex(x, range.to);
+      if (indexFrom > indexTo) [indexFrom, indexTo] = [indexTo, indexFrom];
+      for (let i = indexFrom; i < indexTo; i++) {
+        controlPoints.push(i);
+      }
+    });
+  }
 
-  var { lowerTriangularNonZeros, permutationEncodedArray } = getDeltaMatrix(nbPoints, lambda);
+  let baseline, iteration;
+  let nbPoints = y.length;
+  let l = nbPoints - 1;
+  let sumNegDifferences = Number.MAX_SAFE_INTEGER;
+  let stopCriterion =
+    factorCriterion * y.reduce((sum, e) => Math.abs(e) + sum, 0);
 
-  var sumNegDifferences = Number.MAX_SAFE_INTEGER;
-  for (var iteration = 0; (iteration < maxIterations && Math.abs(sumNegDifferences) > stopCriterion); iteration++) {
-    let [leftHandSide, rightHandSide] = updateSystem(lowerTriangularNonZeros, yData, weights);
+  let { lowerTriangularNonZeros, permutationEncodedArray } = getDeltaMatrix(
+    nbPoints,
+    lambda,
+  );
 
-    let cho = Cholesky.prepare(leftHandSide, nbPoints, permutationEncodedArray);
+  for (
+    iteration = 0;
+    iteration < maxIterations && Math.abs(sumNegDifferences) > stopCriterion;
+    iteration++
+  ) {
+    let [leftHandSide, rightHandSide] = updateSystem(
+      lowerTriangularNonZeros,
+      y,
+      weights,
+    );
 
-    var baseline = cho(rightHandSide);
+    let cho = Cholesky(leftHandSide, nbPoints, permutationEncodedArray);
 
+    baseline = cho(rightHandSide);
+
     sumNegDifferences = 0;
 
-    let difference = yData.map(calculateError);
+    let difference = y.map(calculateError);
 
     let maxNegativeDiff = -1 * Number.MAX_SAFE_INTEGER;
-    for (var i = 1, l = nbPoints - 1; i < l; i++) {
+    for (let i = 1; i < l; i++) {
       let diff = difference[i];
@@ -40,3 +80,3 @@ if (diff >= 0) {
       } else {
-        weights[i] = Math.exp(iteration * diff / sumNegDifferences);
+        weights[i] = Math.exp((iteration * diff) / sumNegDifferences);
         if (maxNegativeDiff < diff) maxNegativeDiff = diff;
@@ -46,12 +86,13 @@ }
 
-    let value = Math.exp(iteration * maxNegativeDiff / sumNegDifferences);
+    let value = Math.exp((iteration * maxNegativeDiff) / sumNegDifferences);
     weights[0] = value;
     weights[l] = value;
+    controlPoints.forEach((i) => (weights[i] = value));
   }
 
   return {
-    corrected: yData.map((e, i) => e - baseline[i]),
+    corrected: y.map((e, i) => e - baseline[i]),
     baseline,
     iteration,
-    error: sumNegDifferences
+    error: sumNegDifferences,
   };
@@ -66,40 +107,2 @@
 
-
-function getDeltaMatrix(nbPoints, lambda) {
-  var matrix = [];
-  for (var i = 0, last = nbPoints - 1; i < last; i++) {
-    matrix.push([i, i, lambda * 2]);
-    matrix.push([i + 1, i, -1 * lambda]);
-  }
-  matrix[0][2] = lambda;
-  matrix.push([last, last, lambda]);
-  return { lowerTriangularNonZeros: matrix, permutationEncodedArray: cuthillMckee(matrix, nbPoints) };
-}
-
-function updateSystem(matrix, yData, weights) {
-  let nbPoints = yData.length;
-  var newMatrix = new Array(matrix.length);
-  var newVector = new Float64Array(nbPoints);
-  for (var i = 0, l = nbPoints - 1; i < l; i++) {
-    let w = weights[i];
-    let diag = i * 2;
-    let next = diag + 1;
-    newMatrix[diag] = matrix[diag].slice();
-    newMatrix[next] = matrix[next].slice();
-    if (w === 0) {
-      newVector[i] = 0;
-    } else {
-      newVector[i] = yData[i] * w;
-      newMatrix[diag][2] += w;
-    }
-  }
-  newVector[l] = yData[l] * weights[l];
-  newMatrix[l * 2] = matrix[l * 2].slice();
-  newMatrix[l * 2][2] += weights[l];
-
-  return [newMatrix, newVector];
-}
-
-
 export { airPLS as default };
-

New alerts

License Policy Violation

License

This package is not allowed per your license policy. Review the package's license to ensure compliance.

Found 1 instance in 1 package

Major refactor

Supply chain risk

Package has recently undergone a major refactor. It may be unstable or indicate significant internal changes. Use caution when updating to versions that include significant changes.

Found 1 instance in 1 package

Fixed alerts

License Policy Violation

License

This package is not allowed per your license policy. Review the package's license to ensure compliance.

Found 1 instance in 1 package

Improved metrics

Total package byte prevSize

30136

increased by195.71%

Dependency count

1

decreased by-50%

Number of package files

9

increased by12.5%

Lines of code

780

increased by326.23%

Number of lines in readme file

44

increased by12.82%

Worsened metrics

Dev dependency count

14

increased by27.27%

Number of medium supply chain risk alerts

1

increased byInfinity%

Dependency changes

• - Removedcholesky-solve@^0.2.1

• - Removedcholesky-solve@0.2.1(transitive)