feat!: new output and a few features (#21)

refactor!: main functions return the same type

refactor!: use named exports in functions

feat: intercepting at other than zero

feat: log some execution information

refactor: destructure fn parameters

refactor: increase test coverage

refactor: add function overloads

Author: newresu

README.md

@@ -1,4 +1,4 @@
-cnnls
+# fcnnls
 
 [![NPM version][npm-image]][npm-url]
 [![build status][ci-image]][ci-url]
@@ -8,73 +8,84 @@ cnnls
 
 Fast Combinatorial Non-negative Least Squares.
 
-Fast algorithm for the solution of large‐scale non‐negativity‐constrained least squares problems from Van Benthem and Keenan ([10.1002/cem.889](http://doi.org/10.1002/cem.889)), based on the active-set method algorithm published by Lawson and Hanson.
+As described in the publication by Van Benthem and Keenan ([10.1002/cem.889](http://doi.org/10.1002/cem.889)), which is in turn based on the active-set method algorithm previously published by Lawson and Hanson. The basic active-set method is implemented in the [nnls repository.](https://github.com/mljs/nnls)
 
-It solves the following optimisation problem.
-Given $\mathbf{X}$ an $n \times l$ matrix and $\mathbf{Y}$ an $n\times p$, find $$\mathrm{argmin}_K ||\mathbf{XK} -\mathbf{Y}||^2_F$$ subject to $\mathbf{K}\geq 0$, where $\mathbf{K}$ is an $l \times p$ matrix and $||\ldots||_F$ is the Frobenius norm. In fact, $\mathbf{K}$ is the best solution to the equation: $\mathbf{XK}=\mathbf{Y}$, where $\mathbf{K} \geq 0$, it performs the regular Non-negative Least Squares algorithm and finds a vector as a solution to the problem. Also, performing this algorithm when $\mathbf{Y}$ is a matrix is like running the algorithm on each column of $\mathbf{Y}$, it will give the same result but in a much more efficient way.
+Given the matrices $\mathbf{X}$ and $\mathbf{Y}$, the code finds the matrix $\mathbf{K}$ that minimises the squared Frobenius norm $$\mathrm{argmin}_K ||\mathbf{XK} -\mathbf{Y}||^2_F$$ subject to $\mathbf{K}\geq 0$.
 
 https://en.wikipedia.org/wiki/Non-negative_least_squares
 
 ## Installation
 
-`$ npm i ml-fcnnls`
+```bash
+npm i ml-fcnnls
+```
 
-## [API Documentation](https://mljs.github.io/fcnnls/)
+## Usage Example
 
-## Usage
+1. Single $y$, using arrays as inputs.
+
+```js
+import { fcnnlsVector } from 'ml-fcnnls';
+
+const X = [
+  [1, 1, 2],
+  [10, 11, -9],
+  [-1, 0, 0],
+  [-5, 6, -7],
+];
+const y = [-1, 11, 0, 1];
+
+const k = fcnnlsVector(X, y).K.to1DArray();
+/* k = [0.4610, 0.5611, 0] */
+```
+
+2. Multiple RHS, using `Matrix` instances as inputs.
 
 ```js
-import { Matrix } from 'ml-matrix';
 import { fcnnls } from 'ml-fcnnls';
+import { Matrix } from 'ml-matrix'; //npm i ml-matrix
 
 // Example with multiple RHS
 
-let X = new Matrix([
+const X = new Matrix([
   [1, 1, 2],
   [10, 11, -9],
   [-1, 0, 0],
   [-5, 6, -7],
 ]);
 
-// Y can either be a Matrix of an array of array
-let Y = new Matrix([
+// Y can either be a Matrix or an array of arrays
+const Y = new Matrix([
   [-1, 0, 0, 9],
   [11, -20, 103, 5],
   [0, 0, 0, 0],
   [1, 2, 3, 4],
 ]);
 
-let K = fcnnls(X, Y);
-
+const K = fcnnls(X, Y).K;
+// `K.to2DArray()` converts the matrix to array.
 /*
 K = Matrix([
-  [0.461, 0, 4.9714, 0],
+  [0.4610, 0, 4.9714, 0],
   [0.5611, 0, 4.7362, 2.2404],
   [0, 1.2388, 0, 1.9136],
-    ])
+])
 */
+```
 
-import { fcnnlsVector } from 'ml-fcnnls';
-
-// Example with single RHS and same X
-// Should be giving a vector with the element of the first column of K in the previous example, since y is the first column of Y
-
-let X = new Matrix([
-  [1, 1, 2],
-  [10, 11, -9],
-  [-1, 0, 0],
-  [-5, 6, -7],
-]);
-
-let y = [-1, 11, 0, 1];
-
-let k = fcnnlsVector(X, y);
+3. Using the options
 
-/*
-k = [0.461, 0.5611, 0]
-*/
+```js
+const K = fcnnls(X, Y, {
+  info: true, // returns the error/iteration.
+  maxIterations: 5,
+  gradientTolerance: 0,
+});
+/* same result than 2*/
 ```
 
+## [API Documentation](https://mljs.github.io/fcnnls/)
+
 ## License
 
 [MIT](./LICENSE)

package.json

@@ -42,16 +42,15 @@
   },
   "homepage": "https://github.com/mljs/fcnnls#readme",
   "devDependencies": {
-    "@types/jest": "^29.5.3",
-    "@vitest/coverage-v8": "^0.33.0",
-    "eslint": "^8.45.0",
+    "@vitest/coverage-v8": "^0.34.5",
+    "eslint": "^8.50.0",
     "eslint-config-cheminfo-typescript": "^12.0.4",
-    "prettier": "^3.0.0",
-    "rimraf": "^5.0.1",
-    "typescript": "^5.1.6",
-    "vitest": "^0.33.0"
+    "prettier": "^3.0.3",
+    "rimraf": "^5.0.5",
+    "typescript": "^5.2.2",
+    "vitest": "^0.34.5"
   },
   "dependencies": {
-    "ml-matrix": "^6.10.4"
+    "ml-matrix": "^6.10.5"
   }
 }

src/__tests__/assertResult.ts

@@ -1,8 +1,21 @@
-import { type Matrix } from 'ml-matrix';
+import { Matrix } from 'ml-matrix';
 import { expect } from 'vitest';
 
-// used for most tests here and in fcnnlsVector.test.ts
-export function assertResult(result: Matrix, solution: Matrix, precision = 4) {
+import { type FcnnlsOutput } from '../fcnnls';
+
+/**
+ *  We only use the value K from the output.
+ * @param output - The output of the fcnnls function
+ * @param solution - The expected solution
+ * @param precision - Number of digits to match
+ */
+export function assertResult(
+  output: FcnnlsOutput,
+  solution: Matrix,
+  precision = 4,
+) {
+  const result = output.K;
+  solution = Matrix.checkMatrix(solution);
   for (let i = 0; i < result.rows; i++) {
     for (let j = 0; j < result.columns; j++) {
       const sol = solution.get(i, j);

src/__tests__/cssls.test.ts

@@ -5,14 +5,25 @@ import { cssls } from '../cssls';
 import { initialisation } from '../initialisation';
 
 import { assertResult } from './assertResult';
+import { prepareInput } from './prepareInitInput';
 
+function csslsResult(X: Matrix, Y: Matrix, useNull?: boolean) {
+  const input = prepareInput(X, Y);
+  const init = initialisation(input);
+  return cssls({
+    XtX: input.XtX,
+    XtY: input.XtY,
+    Pset: useNull ? null : init.Pset,
+    nColsX: X.columns,
+    nColsY: Y.columns,
+  });
+}
 describe('cssls test', () => {
   it('identity X, Y 4x1', () => {
     const X = Matrix.eye(4);
     const Y = new Matrix([[0], [1], [2], [3]]);
-    const { l, p, XtX, XtY, Pset } = initialisation(X, Y);
     const solution = new Matrix([[0], [1], [2], [3]]);
-    const result = cssls(XtX, XtY, Pset, l, p);
+    const result = { K: csslsResult(X, Y) };
     assertResult(result, solution);
   });
   it('identity X, Y 5x3', () => {
@@ -24,15 +35,14 @@ describe('cssls test', () => {
       [3, 8, 13],
       [4, 9, 14],
     ]);
-    const init = initialisation(X, Y);
     const solution = new Matrix([
       [0, 5, 10],
       [1, 6, 11],
       [2, 7, 12],
       [3, 8, 13],
       [4, 9, 14],
     ]);
-    const result = cssls(init.XtX, init.XtY, init.Pset, init.l, init.p);
+    const result = { K: csslsResult(X, Y) };
     assertResult(result, solution);
   });
   it('non-singular square X, Y 3x1', () => {
@@ -42,22 +52,20 @@ describe('cssls test', () => {
       [1, 1, 0],
     ]);
     const Y = new Matrix([[-1], [2], [-3]]);
-    const { l, p, XtX, XtY, Pset } = initialisation(X, Y);
     const solution = new Matrix([[-1], [0], [1]]);
-    const result = cssls(XtX, XtY, Pset, l, p);
+    const result = { K: csslsResult(X, Y) };
     assertResult(result, solution);
   });
-  it('ill-conditionned square X rank 2, Y 3x1', () => {
+  it('ill-conditioned square X rank 2, Y 3x1', () => {
     const X = new Matrix([
       [1, 2, 3],
       [4, 5, 6],
       [7, 8, 9],
     ]);
 
     const Y = new Matrix([[-1], [0], [10]]);
-    const { l, p, XtX, XtY, Pset } = initialisation(X, Y);
     const solution = new Matrix([[1.0455], [0], [0]]);
-    const result = cssls(XtX, XtY, Pset, l, p);
+    const result = { K: csslsResult(X, Y) };
     assertResult(result, solution);
   });
   it('6x3 X full-rank, Y 6x7', () => {
@@ -77,15 +85,13 @@ describe('cssls test', () => {
       [1000, 2, 56, 40, 1, 1, 3],
       [7, 6, 5, 4, 3, 2, 1],
     ]);
-    const { l, p, XtX, XtY, Pset } = initialisation(X, Y);
+
+    const result = { K: csslsResult(X, Y) };
     const solution = new Matrix([
       [203.7567, 0, 0, 0, 0, 0, 0],
       [-149.1338, 3.3309, 2.0243, 1.5134, 0, 0, 0],
       [0, 0, 0, 0, 1.0827, 0.3911, 0.4738],
     ]);
-    const result = Matrix.round(cssls(XtX, XtY, Pset, l, p).mul(10000)).mul(
-      0.0001,
-    );
     assertResult(result, solution);
   });
   it('Van Benthem - Keenan example', () => {
@@ -102,23 +108,19 @@ describe('cssls test', () => {
       [41, 61, 39],
     ]);
 
-    const { l, p, XtX, XtY, Pset } = initialisation(X, Y);
+    const result = { K: csslsResult(X, Y) };
     const solution = new Matrix([
       [0, 0.6873, 0.2836],
       [0.6272, 0, 0.2862],
       [0.3517, 0.2873, 0.335],
     ]);
-    const result = Matrix.round(cssls(XtX, XtY, Pset, l, p).mul(10000)).mul(
-      0.0001,
-    );
     assertResult(result, solution);
   });
   it('negative identity X, positive Y', () => {
     const X = Matrix.eye(3).mul(-1);
     const Y = new Matrix([[1], [2], [3]]);
-    const init = initialisation(X, Y);
     const solution = new Matrix([[-1], [-2], [-3]]);
-    const result = cssls(init.XtX, init.XtY, null, init.l, init.p);
+    const result = { K: csslsResult(X, Y, true) };
     assertResult(result, solution);
   });
 
@@ -129,9 +131,8 @@ describe('cssls test', () => {
       [1, 2, 2, 1],
     ]);
     const Y = new Matrix([[-2], [2], [0]]);
-    const init = initialisation(X, Y);
     const solution = new Matrix([[-2], [0], [0], [2]]);
-    const result = cssls(init.XtX, init.XtY, null, init.l, init.p);
+    const result = { K: csslsResult(X, Y, true) };
     assertResult(result, solution);
   });
 
@@ -142,9 +143,8 @@ describe('cssls test', () => {
       [1, 2, 2, 1],
     ]);
     const Y = new Matrix([[-2], [2], [0]]);
-    const init = initialisation(X, Y);
     const solution = new Matrix([[0], [0], [0], [1]]);
-    const result = cssls(init.XtX, init.XtY, init.Pset, init.l, init.p);
+    const result = { K: csslsResult(X, Y) };
     assertResult(result, solution);
   });
 
@@ -162,7 +162,6 @@ describe('cssls test', () => {
       [0, 23, 23, 23, 0, 0, 0, 23, 23],
     ]);
     const Y = new Matrix([[1], [2], [3], [4], [5], [6], [7], [8], [9], [10]]);
-    const init = initialisation(X, Y);
     const solution = new Matrix([
       [-3.62e-1],
       [1.16712],
@@ -174,7 +173,7 @@ describe('cssls test', () => {
       [0.806154],
       [-4.54969],
     ]);
-    const result = cssls(init.XtX, init.XtY, null, init.l, init.p);
+    const result = { K: csslsResult(X, Y, true) };
     assertResult(result, solution, 3);
   });
 });

src/__tests__/deconv-examples/data_example1.js src/__tests__/deconv-examples/data_example1.ts

@@ -0,0 +1,51 @@
+//plane shifted 5 units
+
+const X = [
+  [0, 0],
+  [1, 3],
+  [2, 4],
+  [3, 6],
+  [4, 8],
+  [5, 10],
+  [6, 12],
+  [7, 14],
+  [8, 16],
+  [9, 18],
+  [10, 20],
+  [11, 22],
+  [12, 24],
+  [13, 26],
+  [14, 28],
+  [15, 30],
+  [16, 32],
+  [17, 34],
+  [18, 36],
+  [19, 38],
+  [20, 40],
+];
+const Y = [
+  0 + 5,
+  2.01 + 5,
+  3 + 5,
+  4.5 + 5,
+  5.8 + 5,
+  7.1 + 5,
+  9.05 + 5,
+  10.5 + 5,
+  12 + 5,
+  13.5 + 5,
+  15 + 5,
+  16.5 + 5,
+  18 + 5,
+  19.5 + 5,
+  21.4 + 5,
+  22.5 + 5,
+  24 + 5,
+  25.5 + 5,
+  27 + 5,
+  28.5 + 5,
+  30 + 5,
+];
+
+const X5 = X.map((x) => [1, ...x]);
+export const data = { X, X5, Y };