nmr-simulation - npm Package Compare versions

Comparing version 1.0.19 to 1.0.21

lib/index.js

 'use strict';
 
-Object.defineProperty(exports, "__esModule", {
-  value: true
-});
+Object.defineProperty(exports, '__esModule', { value: true });
 
-var _SpinSystem = require('./SpinSystem');
+function _interopDefault (ex) { return (ex && (typeof ex === 'object') && 'default' in ex) ? ex['default'] : ex; }
 
-Object.defineProperty(exports, 'SpinSystem', {
-  enumerable: true,
-  get: function () {
-    return _interopRequireDefault(_SpinSystem).default;
+var Matrix = require('ml-matrix');
+var Matrix__default = _interopDefault(Matrix);
+var newArray = _interopDefault(require('new-array'));
+var simpleClustering = _interopDefault(require('ml-simple-clustering'));
+var hlClust = _interopDefault(require('ml-hclust'));
+var SparseMatrix = _interopDefault(require('ml-sparse-matrix'));
+var binarySearch = _interopDefault(require('binary-search'));
+var numSort = require('num-sort');
+
+class SpinSystem {
+  constructor(chemicalShifts, couplingConstants, multiplicity) {
+    this.chemicalShifts = chemicalShifts;
+    this.couplingConstants = couplingConstants;
+    this.multiplicity = multiplicity;
+    this.nSpins = chemicalShifts.length;
+    this._initConnectivity();
+    this._initClusters();
   }
-});
 
-var _simulate1D = require('./simulate1D');
+  static fromSpinusPrediction(result) {
+    let lines = result.split('\n');
+    let nspins = lines.length - 1;
+    let cs = new Array(nspins);
+    let integrals = new Array(nspins);
+    let ids = {};
+    let jc = Matrix__default.zeros(nspins, nspins);
+    for (let i = 0; i < nspins; i++) {
+      var tokens = lines[i].split('\t');
+      cs[i] = Number(tokens[2]);
+      ids[tokens[0] - 1] = i;
+      integrals[i] = Number(tokens[5]); // Is it always 1??
+    }
+    for (let i = 0; i < nspins; i++) {
+      tokens = lines[i].split('\t');
+      let nCoup = (tokens.length - 4) / 3;
+      for (j = 0; j < nCoup; j++) {
+        let withID = tokens[4 + 3 * j] - 1;
+        let idx = ids[withID];
+        // jc[i][idx] = +tokens[6 + 3 * j];
+        jc.set(i, idx, Number(tokens[6 + 3 * j]));
+      }
+    }
 
-Object.defineProperty(exports, 'simulate1D', {
-  enumerable: true,
-  get: function () {
-    return _interopRequireDefault(_simulate1D).default;
+    for (var j = 0; j < nspins; j++) {
+      for (let i = j; i < nspins; i++) {
+        jc.set(j, i, jc.get(i, j));
+      }
+    }
+    return new SpinSystem(cs, jc, newArray(nspins, 2));
   }
-});
 
-var _simulate2D = require('./simulate2D');
+  static fromPrediction(input) {
+    let predictions = SpinSystem.ungroupAtoms(input);
+    const nSpins = predictions.length;
+    const cs = new Array(nSpins);
+    const jc = Matrix__default.zeros(nSpins, nSpins);
+    const multiplicity = new Array(nSpins);
+    const ids = {};
+    let i, k, j;
+    for (i = 0; i < nSpins; i++) {
+      cs[i] = predictions[i].delta;
+      ids[predictions[i].atomIDs[0]] = i;
+    }
+    for (i = 0; i < nSpins; i++) {
+      cs[i] = predictions[i].delta;
+      j = predictions[i].j;
+      for (k = 0; k < j.length; k++) {
+        // jc[ids[predictions[i].atomIDs[0]]][ids[j[k].assignment]] = j[k].coupling;
+        // jc[ids[j[k].assignment]][ids[predictions[i].atomIDs[0]]] = j[k].coupling;
+        jc.set(
+          ids[predictions[i].atomIDs[0]],
+          ids[j[k].assignment],
+          j[k].coupling,
+        );
+        jc.set(
+          ids[j[k].assignment],
+          ids[predictions[i].atomIDs[0]],
+          j[k].coupling,
+        );
+      }
+      multiplicity[i] = predictions[i].integral + 1;
+    }
+    return new SpinSystem(cs, jc, multiplicity);
+  }
 
-Object.defineProperty(exports, 'simulate2D', {
-  enumerable: true,
-  get: function () {
-    return _interopRequireDefault(_simulate2D).default;
+  static ungroupAtoms(prediction) {
+    let result = [];
+    prediction.forEach((pred) => {
+      let atomIDs = pred.atomIDs;
+      if (atomIDs instanceof Array) {
+        for (let i = 0; i < atomIDs.length; i++) {
+          let tempPred = JSON.parse(JSON.stringify(pred));
+          let nmrJ = [];
+          tempPred.atomIDs = [atomIDs[i]];
+          tempPred.integral = 1;
+          if (tempPred.j instanceof Array) {
+            for (let j = 0; j < tempPred.j.length; j++) {
+              let assignment = tempPred.j[j].assignment;
+              if (assignment instanceof Array) {
+                for (let k = 0; k < assignment.length; k++) {
+                  let tempJ = JSON.parse(JSON.stringify(tempPred.j[j]));
+                  tempJ.assignment = assignment[k];
+                  nmrJ.push(tempJ);
+                }
+              }
+            }
+          }
+          tempPred.j = nmrJ;
+          delete tempPred.nbAtoms;
+          result.push(tempPred);
+        }
+      }
+    });
+
+    return result;
   }
-});
 
-function _interopRequireDefault(obj) { return obj && obj.__esModule ? obj : { default: obj }; }
+  _initClusters() {
+    this.clusters = simpleClustering(this.connectivity.to2DArray(), {
+      out: 'indexes',
+    });
+  }
+
+  //I am asumming that couplingConstants is a square matrix
+  _initConnectivity() {
+    const couplings = this.couplingConstants;
+    const connectivity = Matrix__default.ones(couplings.rows, couplings.rows);
+    for (let i = 0; i < couplings.rows; i++) {
+      for (let j = i; j < couplings.columns; j++) {
+        if (couplings.get(i, j) === 0) {
+          connectivity.set(i, j, 0);
+          connectivity.set(j, i, 0);
+        }
+      }
+    }
+    this.connectivity = connectivity;
+  }
+
+  _calculateBetas(J, frequency) {
+    let betas = Matrix__default.zeros(J.rows, J.rows);
+    // Before clustering, we must add hidden J, we could use molecular information if available
+    let i, j;
+    for (i = 0; i < J.rows; i++) {
+      for (j = i; j < J.columns; j++) {
+        let element = J.get(i, j);
+        if (this.chemicalShifts[i] - this.chemicalShifts[j] !== 0) {
+          let value =
+            1 -
+            Math.abs(
+              element /
+                ((this.chemicalShifts[i] - this.chemicalShifts[j]) * frequency),
+            );
+          betas.set(i, j, value);
+          betas.set(j, i, value);
+        } else if (!(i === j || element !== 0)) {
+          betas.set(i, j, 1);
+          betas.set(j, i, 1);
+        }
+      }
+    }
+    return betas;
+  }
+
+  ensureClusterSize(options) {
+    let betas = this._calculateBetas(
+      this.couplingConstants,
+      options.frequency || 400,
+    );
+    let cluster = hlClust.agnes(betas.to2DArray(), { isDistanceMatrix: true });
+    let list = [];
+    this._splitCluster(cluster, list, options.maxClusterSize || 8, false);
+    let clusters = this._mergeClusters(list);
+    this.nClusters = clusters.length;
+    this.clusters = new Array(clusters.length);
+
+    for (let j = 0; j < this.nClusters; j++) {
+      this.clusters[j] = [];
+      for (let i = 0; i < this.nSpins; i++) {
+        let element = clusters[j][i];
+        if (element !== 0) {
+          if (element < 0) {
+            this.clusters[j].push(-(i + 1));
+          } else {
+            this.clusters[j].push(i);
+          }
+        }
+      }
+    }
+  }
+
+  /**
+   * Recursively split the clusters until the maxClusterSize criteria has been ensured.
+   * @param {Array} cluster
+   * @param {Array} clusterList
+   * @param {number} maxClusterSize
+   * @param  {boolean} force
+   */
+  _splitCluster(cluster, clusterList, maxClusterSize, force) {
+    if (!force && cluster.index.length <= maxClusterSize) {
+      clusterList.push(this._getMembers(cluster));
+    } else {
+      for (let child of cluster.children) {
+        if (!isNaN(child.index) || child.index.length <= maxClusterSize) {
+          let members = this._getMembers(child);
+          // Add the neighbors that shares at least 1 coupling with the given cluster
+          let count = 0;
+          for (let i = 0; i < this.nSpins; i++) {
+            if (members[i] === 1) {
+              count++;
+              for (let j = 0; j < this.nSpins; j++) {
+                if (this.connectivity.get(i, j) === 1 && members[j] === 0) {
+                  members[j] = -1;
+                  count++;
+                }
+              }
+            }
+          }
+
+          if (count <= maxClusterSize) {
+            clusterList.push(members);
+          } else {
+            if (isNaN(child.index)) {
+              this._splitCluster(child, clusterList, maxClusterSize, true);
+            } else {
+              // We have to threat this spin alone and use the resurrection algorithm instead of the simulation
+              members[child.index] = 2;
+              clusterList.push(members);
+            }
+          }
+        } else {
+          this._splitCluster(child, clusterList, maxClusterSize, false);
+        }
+      }
+    }
+  }
+  /**
+   * Recursively gets the cluster members
+   * @param cluster
+   * @param members
+   */
+
+  _getMembers(cluster) {
+    let members = new Array(this.nSpins);
+    for (let i = 0; i < this.nSpins; i++) {
+      members[i] = 0;
+    }
+    if (!isNaN(cluster.index)) {
+      members[cluster.index * 1] = 1;
+    } else {
+      for (let index of cluster.index) {
+        members[index.index * 1] = 1;
+      }
+    }
+    return members;
+  }
+
+  _mergeClusters(list) {
+    let nElements = 0;
+    let clusterA, clusterB, i, j, index, common, count;
+    for (i = list.length - 1; i >= 0; i--) {
+      clusterA = list[i];
+      nElements = clusterA.length;
+      index = 0;
+
+      // Is it a candidate to be merged?
+      while (index < nElements && clusterA[index++] !== -1);
+
+      if (index < nElements) {
+        for (j = list.length - 1; j >= i + 1; j--) {
+          clusterB = list[j];
+          // Do they have common elements?
+          index = 0;
+          common = 0;
+          count = 0;
+          while (index < nElements) {
+            if (clusterA[index] * clusterB[index] === -1) {
+              common++;
+            }
+            if (clusterA[index] !== 0 || clusterB[index] !== 0) {
+              count++;
+            }
+            index++;
+          }
+
+          if (common > 0 && count <= this.maxClusterSize) {
+            // Then we can merge those 2 clusters
+            index = 0;
+            while (index < nElements) {
+              if (clusterB[index] === 1) {
+                clusterA[index] = 1;
+              } else {
+                if (clusterB[index] === -1 && clusterA[index] !== 1) {
+                  clusterA[index] = -1;
+                }
+              }
+              index++;
+            }
+            // list.remove(clusterB);
+            list.splice(j, 1);
+            j++;
+          }
+        }
+      }
+    }
+
+    return list;
+  }
+}
+
+function createPauli(mult) {
+  const spin = (mult - 1) / 2;
+  const prjs = new Array(mult);
+  const temp = new Array(mult);
+  for (var i = 0; i < mult; i++) {
+    prjs[i] = mult - 1 - i - spin;
+    temp[i] = Math.sqrt(spin * (spin + 1) - prjs[i] * (prjs[i] + 1));
+  }
+  const p = diag(temp, 1, mult, mult);
+  for (i = 0; i < mult; i++) {
+    temp[i] = Math.sqrt(spin * (spin + 1) - prjs[i] * (prjs[i] - 1));
+  }
+  const m = diag(temp, -1, mult, mult);
+  const x = p
+    .clone()
+    .add(m)
+    .mul(0.5);
+  const y = m
+    .clone()
+    .mul(-1)
+    .add(p)
+    .mul(-0.5);
+  const z = diag(prjs, 0, mult, mult);
+  return { x, y, z, m, p };
+}
+
+function diag(A, d, n, m) {
+  const diag = new SparseMatrix(n, m, { initialCapacity: 20 });
+  for (let i = 0; i < A.length; i++) {
+    if (i - d >= 0 && i - d < n && i < m) {
+      diag.set(i - d, i, A[i]);
+    }
+  }
+  return diag;
+}
+
+const pauli2 = createPauli(2);
+
+function getPauli(mult) {
+  if (mult === 2) return pauli2;
+  else return createPauli(mult);
+}
+
+const smallValue = 1e-2;
+
+/**
+ * This function simulates a one dimensional nmr spectrum. This function returns an array containing the relative intensities of the spectrum in the specified simulation window (from-to).
+ * @param {object} spinSystem - The SpinSystem object to be simulated
+ * @param {object} options
+ * @param {number} options.frequency - The frequency in Mhz of the fake spectrometer that records the spectrum. 400 by default
+ * @param {number} options.from - The low limit of the ordinate variable. 0 by default
+ * @param {number} options.to - The upper limit of the ordinate variable. 10 by default|
+ * @param {number} options.lineWidth - The linewidth of the output spectrum, expresed in Hz. 1Hz by default
+ * @param {number} options.nbPoints - Number of points of the output spectrum. 1024 by default
+ * @param {number} options.maxClusterSize - Maximum number of atoms on each cluster that can be considered to be simulated together. It affects the the quality and speed of the simulation. 10 by default
+ * @param {number} options.output - ['y' or 'xy'] it specify the output format. if 'y' is specified, the output of the simulation will be a single vector containing the y data of the spectrum. if 'xy' is specified, the output of the simulation will be an object containing {x,[], y:[]}, the x, y of the spectrum. 'y' by default
+ * @param {number} options.lortogauRatio - How much of the shape is Lorentzian relative to Gaussian - default : 0.5
+ * @return {object}
+ */
+function simulate1d(spinSystem, options) {
+  let i, j;
+  let {
+    lineWidth = 1,
+    nbPoints = 1024,
+    maxClusterSize = 10,
+    output = 'y',
+    frequency: frequencyMHz = 400,
+    noiseFactor = 1,
+    lortogauRatio = 0.5
+  } = options;
+
+  nbPoints = Number(nbPoints);
+
+  const from = options.from * frequencyMHz || 0;
+  const to = (options.to || 10) * frequencyMHz;
+
+  const chemicalShifts = spinSystem.chemicalShifts.slice();
+  for (i = 0; i < chemicalShifts.length; i++) {
+    chemicalShifts[i] = chemicalShifts[i] * frequencyMHz;
+  }
+
+  // Prepare pseudo voigt
+  let lineWidthPointsG = lortogauRatio * (nbPoints * lineWidth) / Math.abs(to - from) / 2.355;
+  let lineWidthPointsL = (1 - lortogauRatio) * (nbPoints * lineWidth) / Math.abs(to - from) / 2;
+  let lnPoints = lineWidthPointsL * 40;
+
+  const gaussianLength = lnPoints | 0;
+  const gaussian = new Array(gaussianLength);
+  const b = lnPoints / 2;
+  const c = lineWidthPointsG * lineWidthPointsG * 2;
+  const l2 = lineWidthPointsL * lineWidthPointsL;
+  const g2pi = lineWidthPointsG * Math.sqrt(2 * Math.PI);
+  for (i = 0; i < gaussianLength; i++) {
+    let x2 = (i - b) * (i - b);
+    gaussian[i] =
+      10e9 *
+      (Math.exp(-x2 / c) / g2pi + lineWidthPointsL / ((x2 + l2) * Math.PI));
+  }
+
+  let result = options.withNoise
+    ? [...new Array(nbPoints)].map(() => Math.random() * noiseFactor)
+    : new Array(nbPoints).fill(0);
+
+  const multiplicity = spinSystem.multiplicity;
+  for (let h = 0; h < spinSystem.clusters.length; h++) {
+    const cluster = spinSystem.clusters[h];
+
+    let clusterFake = new Array(cluster.length);
+    for (i = 0; i < cluster.length; i++) {
+      clusterFake[i] = cluster[i] < 0 ? -cluster[i] - 1 : cluster[i];
+    }
+
+    let weight = 1;
+    var sumI = 0;
+    const frequencies = [];
+    const intensities = [];
+    if (cluster.length > maxClusterSize) {
+      // This is a single spin, but the cluster exceeds the maxClusterSize criteria
+      // we use the simple multiplicity algorithm
+      // Add the central peak. It will be split with every single J coupling.
+      let index = 0;
+      while (cluster[index++] < 0);
+      index = cluster[index - 1];
+      var currentSize, jc;
+      frequencies.push(-chemicalShifts[index]);
+      for (i = 0; i < cluster.length; i++) {
+        if (cluster[i] < 0) {
+          jc = spinSystem.couplingConstants.get(index, clusterFake[i]) / 2;
+          currentSize = frequencies.length;
+          for (j = 0; j < currentSize; j++) {
+            frequencies.push(frequencies[j] + jc);
+            frequencies[j] -= jc;
+          }
+        }
+      }
+
+      frequencies.sort(numSort.asc);
+      sumI = frequencies.length;
+      weight = 1;
+
+      for (i = 0; i < sumI; i++) {
+        intensities.push(1);
+      }
+    } else {
+      const hamiltonian = getHamiltonian(
+        chemicalShifts,
+        spinSystem.couplingConstants,
+        multiplicity,
+        spinSystem.connectivity,
+        clusterFake,
+      );
+      const hamSize = hamiltonian.rows;
+      const evd = new Matrix.EVD(hamiltonian);
+      const V = evd.eigenvectorMatrix;
+      const diagB = evd.realEigenvalues;
+      const assignmentMatrix = new SparseMatrix(hamSize, hamSize);
+      const multLen = cluster.length;
+      weight = 0;
+      for (let n = 0; n < multLen; n++) {
+        const L = getPauli(multiplicity[clusterFake[n]]);
+
+        let temp = 1;
+        for (j = 0; j < n; j++) {
+          temp *= multiplicity[clusterFake[j]];
+        }
+        const A = SparseMatrix.eye(temp);
+
+        temp = 1;
+        for (j = n + 1; j < multLen; j++) {
+          temp *= multiplicity[clusterFake[j]];
+        }
+        const B = SparseMatrix.eye(temp);
+        const tempMat = A.kroneckerProduct(L.m).kroneckerProduct(B);
+        if (cluster[n] >= 0) {
+          assignmentMatrix.add(tempMat.mul(cluster[n] + 1));
+          weight++;
+        } else {
+          assignmentMatrix.add(tempMat.mul(cluster[n]));
+        }
+      }
+
+      let rhoip = Matrix.Matrix.zeros(hamSize, hamSize);
+      assignmentMatrix.forEachNonZero((i, j, v) => {
+        if (v > 0) {
+          for (let k = 0; k < V.columns; k++) {
+            let element = V.get(j, k);
+            if (element !== 0) {
+              rhoip.set(i, k, rhoip.get(i, k) + element);
+            }
+          }
+        }
+        return v;
+      });
+
+      let rhoip2 = rhoip.clone();
+      assignmentMatrix.forEachNonZero((i, j, v) => {
+        if (v < 0) {
+          for (let k = 0; k < V.columns; k++) {
+            let element = V.get(j, k);
+            if (element !== 0) {
+              rhoip2.set(i, k, rhoip2.get(i, k) + element);
+            }
+          }
+        }
+        return v;
+      });
+      const tV = V.transpose();
+
+      rhoip = tV.mmul(rhoip);
+      rhoip = new SparseMatrix(rhoip.to2DArray(), { threshold: smallValue });
+      triuTimesAbs(rhoip, smallValue);
+      rhoip2 = tV.mmul(rhoip2);
+
+      rhoip2 = new SparseMatrix(rhoip2.to2DArray(), { threshold: smallValue });
+      rhoip2.forEachNonZero((i, j, v) => {
+        return v;
+      });
+      triuTimesAbs(rhoip2, smallValue);
+      // eslint-disable-next-line no-loop-func
+      rhoip2.forEachNonZero((i, j, v) => {
+        let val = rhoip.get(i, j);
+        val = Math.min(Math.abs(val), Math.abs(v));
+        val *= val;
+
+        sumI += val;
+        let valFreq = diagB[i] - diagB[j];
+        let insertIn = binarySearch(frequencies, valFreq, numSort.asc);
+        if (insertIn < 0) {
+          frequencies.splice(-1 - insertIn, 0, valFreq);
+          intensities.splice(-1 - insertIn, 0, val);
+        } else {
+          intensities[insertIn] += val;
+        }
+      });
+    }
+    const numFreq = frequencies.length;
+    if (numFreq > 0) {
+      weight = weight / sumI;
+      const diff = lineWidth / 64;
+      let valFreq = frequencies[0];
+      let inte = intensities[0];
+      let count = 1;
+      for (i = 1; i < numFreq; i++) {
+        if (Math.abs(frequencies[i] - valFreq / count) < diff) {
+          inte += intensities[i];
+          valFreq += frequencies[i];
+          count++;
+        } else {
+          addPeak(
+            result,
+            valFreq / count,
+            inte * weight,
+            from,
+            to,
+            nbPoints,
+            gaussian,
+          );
+          valFreq = frequencies[i];
+          inte = intensities[i];
+          count = 1;
+        }
+      }
+      addPeak(
+        result,
+        valFreq / count,
+        inte * weight,
+        from,
+        to,
+        nbPoints,
+        gaussian,
+      );
+    }
+  }
+  if (output === 'xy') {
+    return { x: _getX(options.from, options.to, nbPoints), y: result };
+  }
+  if (output === 'y') {
+    return result;
+  }
+  throw new RangeError('wrong output option');
+}
+/**
+ * Add a new peak to the current array
+ * @param {Array} result - Array of numbers
+ * @param {number} freq - center of the peak
+ * @param {*} height - peak height
+ * @param {*} from - start point of the peak
+ * @param {*} to - end point of the peak
+ * @param {*} nbPoints - number of points to add
+ * @param {*} gaussian - Shape to fill with
+ * @return {undefined}
+ */
+function addPeak(result, freq, height, from, to, nbPoints, gaussian) {
+  const center = ((nbPoints * (-freq - from)) / (to - from)) | 0;
+  const lnPoints = gaussian.length;
+  let index = 0;
+  let indexLorentz = 0;
+  for (let i = center - lnPoints / 2; i < center + lnPoints / 2; i++) {
+    index = i | 0;
+    if (i >= 0 && i < nbPoints) {
+      result[index] = result[index] + gaussian[indexLorentz] * height;
+    }
+    indexLorentz++;
+  }
+}
+
+function triuTimesAbs(A, val) {
+  A.forEachNonZero((i, j, v) => {
+    if (i > j) return 0;
+    if (Math.abs(v) <= val) return 0;
+    return v;
+  });
+}
+/**
+ * Create a hamiltonian matrix for the given spinsystem
+ * @param {Array} chemicalShifts - An array containing the chemical shift in Hz
+ * @param {Array} couplingConstants - An array containing the coupling constants in Hz
+ * @param {Array} multiplicity - An array specifiying the multiplicities of each scalar coupling
+ * @param {Array} conMatrix - A one step connectivity matrix for the given spin system
+ * @param {Array} cluster - An binary array specifiying the spins to be considered for this hamiltonial
+ * @return {object}
+ */
+function getHamiltonian(
+  chemicalShifts,
+  couplingConstants,
+  multiplicity,
+  conMatrix,
+  cluster,
+) {
+  let hamSize = 1;
+  for (let i = 0; i < cluster.length; i++) {
+    hamSize *= multiplicity[cluster[i]];
+  }
+
+  const clusterHam = new SparseMatrix(hamSize, hamSize);
+
+  for (let pos = 0; pos < cluster.length; pos++) {
+    let n = cluster[pos];
+
+    const L = getPauli(multiplicity[n]);
+
+    let A1, B1;
+    let temp = 1;
+    for (let i = 0; i < pos; i++) {
+      temp *= multiplicity[cluster[i]];
+    }
+    A1 = SparseMatrix.eye(temp);
+
+    temp = 1;
+    for (let i = pos + 1; i < cluster.length; i++) {
+      temp *= multiplicity[cluster[i]];
+    }
+    B1 = SparseMatrix.eye(temp);
+
+    const alpha = chemicalShifts[n];
+    const kronProd = A1.kroneckerProduct(L.z).kroneckerProduct(B1);
+    clusterHam.add(kronProd.mul(alpha));
+    for (let pos2 = 0; pos2 < cluster.length; pos2++) {
+      const k = cluster[pos2];
+      if (conMatrix.get(n, k) === 1) {
+        const S = getPauli(multiplicity[k]);
+
+        let A2, B2;
+        let temp = 1;
+        for (let i = 0; i < pos2; i++) {
+          temp *= multiplicity[cluster[i]];
+        }
+        A2 = SparseMatrix.eye(temp);
+
+        temp = 1;
+        for (let i = pos2 + 1; i < cluster.length; i++) {
+          temp *= multiplicity[cluster[i]];
+        }
+        B2 = SparseMatrix.eye(temp);
+
+        const kron1 = A1.kroneckerProduct(L.x)
+          .kroneckerProduct(B1)
+          .mmul(A2.kroneckerProduct(S.x).kroneckerProduct(B2));
+        kron1.add(
+          A1.kroneckerProduct(L.y)
+            .kroneckerProduct(B1)
+            .mul(-1)
+            .mmul(A2.kroneckerProduct(S.y).kroneckerProduct(B2)),
+        );
+        kron1.add(
+          A1.kroneckerProduct(L.z)
+            .kroneckerProduct(B1)
+            .mmul(A2.kroneckerProduct(S.z).kroneckerProduct(B2)),
+        );
+
+        clusterHam.add(kron1.mul(couplingConstants.get(n, k) / 2));
+      }
+    }
+  }
+  return clusterHam;
+}
+
+function _getX(from, to, nbPoints) {
+  const x = new Array(nbPoints);
+  const dx = (to - from) / (nbPoints - 1);
+  for (let i = 0; i < nbPoints; i++) {
+    x[i] = from + i * dx;
+  }
+  return x;
+}
+
+let defOptions = {
+  H: { frequency: 400, lineWidth: 10 },
+  C: { frequency: 100, lineWidth: 10 },
+};
+
+function simule2DNmrSpectrum(table, options) {
+  let i;
+  const fromLabel = table[0].fromAtomLabel;
+  const toLabel = table[0].toLabel;
+  const frequencyX = options.frequencyX || defOptions[fromLabel].frequency;
+  const frequencyY = options.frequencyY || defOptions[toLabel].frequency;
+  let lineWidthX = options.lineWidthX || defOptions[fromLabel].lineWidth;
+  let lineWidthY = options.lineWidthY || defOptions[toLabel].lineWidth;
+  let symmetrize = options.symmetrize || false;
+
+  let sigmaX = lineWidthX / frequencyX;
+  let sigmaY = lineWidthY / frequencyY;
+
+  let minX = table[0].fromChemicalShift;
+  let maxX = table[0].fromChemicalShift;
+  let minY = table[0].toChemicalShift;
+  let maxY = table[0].toChemicalShift;
+  i = 1;
+  while (i < table.length) {
+    minX = Math.min(minX, table[i].fromChemicalShift);
+    maxX = Math.max(maxX, table[i].fromChemicalShift);
+    minY = Math.min(minY, table[i].toChemicalShift);
+    maxY = Math.max(maxY, table[i].toChemicalShift);
+    i++;
+  }
+
+  if (options.firstX !== null && !isNaN(options.firstX)) {
+    minX = options.firstX;
+  }
+  if (options.firstY !== null && !isNaN(options.firstY)) {
+    minY = options.firstY;
+  }
+  if (options.lastX !== null && !isNaN(options.lastX)) {
+    maxX = options.lastX;
+  }
+  if (options.lastY !== null && !isNaN(options.lastY)) {
+    maxY = options.lastY;
+  }
+
+  let nbPointsX = options.nbPointsX || 512;
+  let nbPointsY = options.nbPointsY || 512;
+
+  let spectraMatrix = new Matrix__default(nbPointsY, nbPointsX).fill(0);
+  i = 0;
+  while (i < table.length) {
+    // parameters.couplingConstant = table[i].j;
+    // parameters.pathLength = table[i].pathLength;
+    let peak = {
+      x: unitsToArrayPoints(table[i].fromChemicalShift, minX, maxX, nbPointsX),
+      y: unitsToArrayPoints(table[i].toChemicalShift, minY, maxY, nbPointsY),
+      z: table[i].fromAtoms.length + table[i].toAtoms.length,
+      widthX: unitsToArrayPoints(sigmaX + minX, minX, maxX, nbPointsX),
+      widthY: unitsToArrayPoints(sigmaY + minY, minY, maxY, nbPointsY),
+    };
+    addPeak$1(spectraMatrix, peak);
+    if (symmetrize) {
+      addPeak$1(spectraMatrix, {
+        x: peak.y,
+        y: peak.x,
+        z: peak.z,
+        widthX: peak.widthY,
+        widthY: peak.widthX,
+      });
+    }
+    i++;
+  }
+  return spectraMatrix;
+}
+
+function unitsToArrayPoints(x, from, to, nbPoints) {
+  return ((x - from) * (nbPoints - 1)) / (to - from);
+}
+
+function addPeak$1(matrix, peak) {
+  let nSigma = 4;
+  let fromX = Math.max(0, Math.round(peak.x - peak.widthX * nSigma));
+  // var toX = Math.min(matrix[0].length - 1, Math.round(peak.x + peak.widthX * nSigma));
+  let toX = Math.min(
+    matrix.columns - 1,
+    Math.round(peak.x + peak.widthX * nSigma),
+  );
+  let fromY = Math.max(0, Math.round(peak.y - peak.widthY * nSigma));
+  // var toY = Math.min(matrix.length - 1, Math.round(peak.y + peak.widthY * nSigma));
+  let toY = Math.min(
+    matrix.rows - 1,
+    Math.round(peak.y + peak.widthY * nSigma),
+  );
+
+  let squareSigmaX = peak.widthX * peak.widthX;
+  let squareSigmaY = peak.widthY * peak.widthY;
+  for (let j = fromY; j < toY; j++) {
+    for (let i = fromX; i < toX; i++) {
+      let exponent =
+        Math.pow(peak.x - i, 2) / squareSigmaX +
+        Math.pow(peak.y - j, 2) / squareSigmaY;
+      let result = 10000 * peak.z * Math.exp(-exponent);
+      // matrix[j][i] += result;
+      matrix.set(j, i, matrix.get(j, i) + result);
+    }
+  }
+}
+
+exports.SpinSystem = SpinSystem;
+exports.simulate1D = simulate1d;
+exports.simulate2D = simule2DNmrSpectrum;

package.json

 {
   "name": "nmr-simulation",
-  "version": "1.0.19",
+  "version": "1.0.21",
   "description": "Simulate NMR spectra from spin systems",
@@ -32,13 +32,17 @@ "main": "lib/index.js",
     "ml-hclust": "^1.3.0",
-    "ml-matrix": "^2.3.0",
+    "ml-matrix": "^6.2.0",
     "ml-simple-clustering": "0.1.0",
     "ml-sparse-matrix": "^0.2.1",
+    "ml-stat": "^1.3.3",
     "new-array": "^1.0.0",
-    "num-sort": "^1.0.0",
-    "ml-stat": "^1.3.3"
+    "num-sort": "^1.0.0"
   },
   "devDependencies": {
-    "nmr-predictor": "^1.2.0",
-    "should": "^13.2.3"
-  }
+    "jest": "^24.8.0",
+    "jest-matcher-deep-close-to": "^1.3.0",
+    "ml-array-xy-filter-x": "^1.0.1",
+    "ml-array-xy-max-y": "^1.0.1",
+    "nmr-predictor": "^1.2.1"
+  },
+  "gitHead": "c1a4b6ef4aa3d1ca9927461a0b85d9ebaee0d1fb"
 }

README.md

@@ -0,0 +0,0 @@ # nmr-simulation

src/index.js

 export { default as SpinSystem } from './SpinSystem';
 export { default as simulate1D } from './simulate1D';
 export { default as simulate2D } from './simulate2D';

src/pauli.js

@@ -8,3 +8,3 @@ import SparseMatrix from 'ml-sparse-matrix';
   for (var i = 0; i < mult; i++) {
-    prjs[i] = (mult - 1) - i - spin;
+    prjs[i] = mult - 1 - i - spin;
     temp[i] = Math.sqrt(spin * (spin + 1) - prjs[i] * (prjs[i] + 1));
@@ -17,4 +17,11 @@ }
   const m = diag(temp, -1, mult, mult);
-  const x = p.clone().add(m).mul(0.5);
-  const y = m.clone().mul(-1).add(p).mul(-0.5);
+  const x = p
+    .clone()
+    .add(m)
+    .mul(0.5);
+  const y = m
+    .clone()
+    .mul(-1)
+    .add(p)
+    .mul(-0.5);
   const z = diag(prjs, 0, mult, mult);
@@ -26,4 +33,4 @@ return { x, y, z, m, p };
   const diag = new SparseMatrix(n, m, { initialCapacity: 20 });
-  for (var i = 0; i < A.length; i++) {
-    if ((i - d) >= 0 && (i - d) < n && i < m) {
+  for (let i = 0; i < A.length; i++) {
+    if (i - d >= 0 && i - d < n && i < m) {
       diag.set(i - d, i, A[i]);
@@ -30,0 +37,0 @@ }

src/simulate1D.js

@@ -1,2 +0,2 @@
-import Matrix from 'ml-matrix';
+import { Matrix, EVD } from 'ml-matrix';
 import SparseMatrix from 'ml-sparse-matrix';
@@ -10,3 +10,2 @@ import binarySearch from 'binary-search';
 
-
 /**
@@ -23,7 +22,8 @@ * This function simulates a one dimensional nmr spectrum. This function returns an array containing the relative intensities of the spectrum in the specified simulation window (from-to).
  * @param {number} options.output - ['y' or 'xy'] it specify the output format. if 'y' is specified, the output of the simulation will be a single vector containing the y data of the spectrum. if 'xy' is specified, the output of the simulation will be an object containing {x,[], y:[]}, the x, y of the spectrum. 'y' by default
+ * @param {number} options.lortogauRatio - How much of the shape is Lorentzian relative to Gaussian - default : 0.5
  * @return {object}
  */
 export default function simulate1d(spinSystem, options) {
-  var i, j;
-  var {
+  let i, j;
+  let {
     lineWidth = 1,
@@ -34,3 +34,4 @@ nbPoints = 1024,
     frequency: frequencyMHz = 400,
-    noiseFactor = 1
+    noiseFactor = 1,
+    lortogauRatio = 0.5
   } = options;
@@ -48,5 +49,5 @@
 
-  //Prepare pseudo voigt
-  let lineWidthPointsG = (nbPoints * lineWidth / Math.abs(to - from)) / 2.355;
-  let lineWidthPointsL = (nbPoints * lineWidth / Math.abs(to - from)) / 2;
+  // Prepare pseudo voigt
+  let lineWidthPointsG = lortogauRatio * (nbPoints * lineWidth) / Math.abs(to - from) / 2.355;
+  let lineWidthPointsL = (1 - lortogauRatio) * (nbPoints * lineWidth) / Math.abs(to - from) / 2;
   let lnPoints = lineWidthPointsL * 40;
@@ -62,12 +63,16 @@
     let x2 = (i - b) * (i - b);
-    gaussian[i] = 10e9 * (Math.exp(-x2 / c) / g2pi + lineWidthPointsL / ((x2 + l2) * Math.PI));
+    gaussian[i] =
+      10e9 *
+      (Math.exp(-x2 / c) / g2pi + lineWidthPointsL / ((x2 + l2) * Math.PI));
   }
 
-  var result = options.withNoise ? [...new Array(nbPoints)].map(() => Math.random() * noiseFactor) : new Array(nbPoints).fill(0);
+  let result = options.withNoise
+    ? [...new Array(nbPoints)].map(() => Math.random() * noiseFactor)
+    : new Array(nbPoints).fill(0);
 
   const multiplicity = spinSystem.multiplicity;
-  for (var h = 0; h < spinSystem.clusters.length; h++) {
+  for (let h = 0; h < spinSystem.clusters.length; h++) {
     const cluster = spinSystem.clusters[h];
 
-    var clusterFake = new Array(cluster.length);
+    let clusterFake = new Array(cluster.length);
     for (i = 0; i < cluster.length; i++) {
@@ -77,3 +82,3 @@ clusterFake[i] = cluster[i] < 0 ? -cluster[i] - 1 : cluster[i];
 
-    var weight = 1;
+    let weight = 1;
     var sumI = 0;
@@ -86,3 +91,3 @@ const frequencies = [];
       // Add the central peak. It will be split with every single J coupling.
-      var index = 0;
+      let index = 0;
       while (cluster[index++] < 0);
@@ -94,3 +99,3 @@ index = cluster[index - 1];
         if (cluster[i] < 0) {
-          jc = spinSystem.couplingConstants[index][clusterFake[i]] / 2;
+          jc = spinSystem.couplingConstants.get(index, clusterFake[i]) / 2;
           currentSize = frequencies.length;
@@ -117,7 +122,6 @@ for (j = 0; j < currentSize; j++) {
         spinSystem.connectivity,
-        clusterFake
+        clusterFake,
       );
-
       const hamSize = hamiltonian.rows;
-      const evd = new Matrix.DC.EVD(hamiltonian);
+      const evd = new EVD(hamiltonian);
       const V = evd.eigenvectorMatrix;
@@ -128,3 +132,3 @@ const diagB = evd.realEigenvalues;
       weight = 0;
-      for (var n = 0; n < multLen; n++) {
+      for (let n = 0; n < multLen; n++) {
         const L = getPauli(multiplicity[clusterFake[n]]);
@@ -155,6 +159,6 @@
         if (v > 0) {
-          const row = V[j];
-          for (var k = 0; k < row.length; k++) {
-            if (row[k] !== 0) {
-              rhoip.set(i, k, rhoip.get(i, k) + row[k]);
+          for (let k = 0; k < V.columns; k++) {
+            let element = V.get(j, k);
+            if (element !== 0) {
+              rhoip.set(i, k, rhoip.get(i, k) + element);
             }
@@ -169,6 +173,6 @@ }
         if (v < 0) {
-          const row = V[j];
-          for (var k = 0; k < row.length; k++) {
-            if (row[k] !== 0) {
-              rhoip2.set(i, k, rhoip2.get(i, k) + row[k]);
+          for (let k = 0; k < V.columns; k++) {
+            let element = V.get(j, k);
+            if (element !== 0) {
+              rhoip2.set(i, k, rhoip2.get(i, k) + element);
             }
@@ -179,13 +183,17 @@ }
       });
+      const tV = V.transpose();
 
-      const tV = V.transpose();
       rhoip = tV.mmul(rhoip);
-      rhoip = new SparseMatrix(rhoip, { threshold: smallValue });
+      rhoip = new SparseMatrix(rhoip.to2DArray(), { threshold: smallValue });
       triuTimesAbs(rhoip, smallValue);
       rhoip2 = tV.mmul(rhoip2);
-      rhoip2 = new SparseMatrix(rhoip2, { threshold: smallValue });
+
+      rhoip2 = new SparseMatrix(rhoip2.to2DArray(), { threshold: smallValue });
+      rhoip2.forEachNonZero((i, j, v) => {
+        return v;
+      });
       triuTimesAbs(rhoip2, smallValue);
       // eslint-disable-next-line no-loop-func
       rhoip2.forEachNonZero((i, j, v) => {
-        var val = rhoip.get(i, j);
+        let val = rhoip.get(i, j);
         val = Math.min(Math.abs(val), Math.abs(v));
@@ -195,4 +203,4 @@ val *= val;
         sumI += val;
-        var valFreq = diagB[i] - diagB[j];
-        var insertIn = binarySearch(frequencies, valFreq, sortAsc);
+        let valFreq = diagB[i] - diagB[j];
+        let insertIn = binarySearch(frequencies, valFreq, sortAsc);
         if (insertIn < 0) {
@@ -219,3 +227,11 @@ frequencies.splice(-1 - insertIn, 0, valFreq);
         } else {
-          addPeak(result, valFreq / count, inte * weight, from, to, nbPoints, gaussian);
+          addPeak(
+            result,
+            valFreq / count,
+            inte * weight,
+            from,
+            to,
+            nbPoints,
+            gaussian,
+          );
           valFreq = frequencies[i];
@@ -226,3 +242,11 @@ inte = intensities[i];
       }
-      addPeak(result, valFreq / count, inte * weight, from, to, nbPoints, gaussian);
+      addPeak(
+        result,
+        valFreq / count,
+        inte * weight,
+        from,
+        to,
+        nbPoints,
+        gaussian,
+      );
     }
@@ -247,9 +271,10 @@ }
  * @param {*} gaussian - Shape to fill with
+ * @return {undefined}
  */
 function addPeak(result, freq, height, from, to, nbPoints, gaussian) {
-  const center = (nbPoints * (-freq - from) / (to - from)) | 0;
+  const center = ((nbPoints * (-freq - from)) / (to - from)) | 0;
   const lnPoints = gaussian.length;
-  var index = 0;
-  var indexLorentz = 0;
-  for (var i = center - lnPoints / 2; i < center + lnPoints / 2; i++) {
+  let index = 0;
+  let indexLorentz = 0;
+  for (let i = center - lnPoints / 2; i < center + lnPoints / 2; i++) {
     index = i | 0;
@@ -279,5 +304,11 @@ if (i >= 0 && i < nbPoints) {
  */
-function getHamiltonian(chemicalShifts, couplingConstants, multiplicity, conMatrix, cluster) {
+function getHamiltonian(
+  chemicalShifts,
+  couplingConstants,
+  multiplicity,
+  conMatrix,
+  cluster,
+) {
   let hamSize = 1;
-  for (var i = 0; i < cluster.length; i++) {
+  for (let i = 0; i < cluster.length; i++) {
     hamSize *= multiplicity[cluster[i]];
@@ -288,4 +319,4 @@ }
 
-  for (var pos = 0; pos < cluster.length; pos++) {
-    var n = cluster[pos];
+  for (let pos = 0; pos < cluster.length; pos++) {
+    let n = cluster[pos];
 
@@ -310,6 +341,5 @@ const L = getPauli(multiplicity[n]);
     clusterHam.add(kronProd.mul(alpha));
-
-    for (var pos2 = 0; pos2 < cluster.length; pos2++) {
+    for (let pos2 = 0; pos2 < cluster.length; pos2++) {
       const k = cluster[pos2];
-      if (conMatrix[n][k] === 1) {
+      if (conMatrix.get(n, k) === 1) {
         const S = getPauli(multiplicity[k]);
@@ -330,11 +360,21 @@
 
-        const kron1 = A1.kroneckerProduct(L.x).kroneckerProduct(B1).mmul(A2.kroneckerProduct(S.x).kroneckerProduct(B2));
-        kron1.add(A1.kroneckerProduct(L.y).kroneckerProduct(B1).mul(-1).mmul(A2.kroneckerProduct(S.y).kroneckerProduct(B2)));
-        kron1.add(A1.kroneckerProduct(L.z).kroneckerProduct(B1).mmul(A2.kroneckerProduct(S.z).kroneckerProduct(B2)));
+        const kron1 = A1.kroneckerProduct(L.x)
+          .kroneckerProduct(B1)
+          .mmul(A2.kroneckerProduct(S.x).kroneckerProduct(B2));
+        kron1.add(
+          A1.kroneckerProduct(L.y)
+            .kroneckerProduct(B1)
+            .mul(-1)
+            .mmul(A2.kroneckerProduct(S.y).kroneckerProduct(B2)),
+        );
+        kron1.add(
+          A1.kroneckerProduct(L.z)
+            .kroneckerProduct(B1)
+            .mmul(A2.kroneckerProduct(S.z).kroneckerProduct(B2)),
+        );
 
-        clusterHam.add(kron1.mul(couplingConstants[n][k] / 2));
+        clusterHam.add(kron1.mul(couplingConstants.get(n, k) / 2));
       }
     }
   }
-
   return clusterHam;
@@ -346,3 +386,3 @@ }
   const dx = (to - from) / (nbPoints - 1);
-  for (var i = 0; i < nbPoints; i++) {
+  for (let i = 0; i < nbPoints; i++) {
     x[i] = from + i * dx;
@@ -349,0 +389,0 @@ }

src/simulate2D.js

@@ -5,7 +5,7 @@ import Matrix from 'ml-matrix';
   H: { frequency: 400, lineWidth: 10 },
-  C: { frequency: 100, lineWidth: 10 }
+  C: { frequency: 100, lineWidth: 10 },
 };
 
 export default function simule2DNmrSpectrum(table, options) {
-  var i;
+  let i;
   const fromLabel = table[0].fromAtomLabel;
@@ -15,13 +15,13 @@ const toLabel = table[0].toLabel;
   const frequencyY = options.frequencyY || defOptions[toLabel].frequency;
-  var lineWidthX = options.lineWidthX || defOptions[fromLabel].lineWidth;
-  var lineWidthY = options.lineWidthY || defOptions[toLabel].lineWidth;
+  let lineWidthX = options.lineWidthX || defOptions[fromLabel].lineWidth;
+  let lineWidthY = options.lineWidthY || defOptions[toLabel].lineWidth;
   let symmetrize = options.symmetrize || false;
 
-  var sigmaX = lineWidthX / frequencyX;
-  var sigmaY = lineWidthY / frequencyY;
+  let sigmaX = lineWidthX / frequencyX;
+  let sigmaY = lineWidthY / frequencyY;
 
-  var minX = table[0].fromChemicalShift;
-  var maxX = table[0].fromChemicalShift;
-  var minY = table[0].toChemicalShift;
-  var maxY = table[0].toChemicalShift;
+  let minX = table[0].fromChemicalShift;
+  let maxX = table[0].fromChemicalShift;
+  let minY = table[0].toChemicalShift;
+  let maxY = table[0].toChemicalShift;
   i = 1;
@@ -49,6 +49,6 @@ while (i < table.length) {
 
-  var nbPointsX = options.nbPointsX || 512;
-  var nbPointsY = options.nbPointsY || 512;
+  let nbPointsX = options.nbPointsX || 512;
+  let nbPointsY = options.nbPointsY || 512;
 
-  var spectraMatrix = new Matrix(nbPointsY, nbPointsX).fill(0);
+  let spectraMatrix = new Matrix(nbPointsY, nbPointsX).fill(0);
   i = 0;
@@ -63,7 +63,13 @@ while (i < table.length) {
       widthX: unitsToArrayPoints(sigmaX + minX, minX, maxX, nbPointsX),
-      widthY: unitsToArrayPoints(sigmaY + minY, minY, maxY, nbPointsY)
+      widthY: unitsToArrayPoints(sigmaY + minY, minY, maxY, nbPointsY),
     };
     addPeak(spectraMatrix, peak);
     if (symmetrize) {
-      addPeak(spectraMatrix, { x: peak.y, y: peak.x, z: peak.z, widthX: peak.widthY, widthY: peak.widthX });
+      addPeak(spectraMatrix, {
+        x: peak.y,
+        y: peak.x,
+        z: peak.z,
+        widthX: peak.widthY,
+        widthY: peak.widthX,
+      });
     }
@@ -80,18 +86,28 @@ i++;
 function addPeak(matrix, peak) {
-  var nSigma = 4;
-  var fromX = Math.max(0, Math.round(peak.x - peak.widthX * nSigma));
-  var toX = Math.min(matrix[0].length - 1, Math.round(peak.x + peak.widthX * nSigma));
-  var fromY = Math.max(0, Math.round(peak.y - peak.widthY * nSigma));
-  var toY = Math.min(matrix.length - 1, Math.round(peak.y + peak.widthY * nSigma));
+  let nSigma = 4;
+  let fromX = Math.max(0, Math.round(peak.x - peak.widthX * nSigma));
+  // var toX = Math.min(matrix[0].length - 1, Math.round(peak.x + peak.widthX * nSigma));
+  let toX = Math.min(
+    matrix.columns - 1,
+    Math.round(peak.x + peak.widthX * nSigma),
+  );
+  let fromY = Math.max(0, Math.round(peak.y - peak.widthY * nSigma));
+  // var toY = Math.min(matrix.length - 1, Math.round(peak.y + peak.widthY * nSigma));
+  let toY = Math.min(
+    matrix.rows - 1,
+    Math.round(peak.y + peak.widthY * nSigma),
+  );
 
-  var squareSigmaX = peak.widthX * peak.widthX;
-  var squareSigmaY = peak.widthY * peak.widthY;
-  for (var j = fromY; j < toY; j++) {
-    for (var i = fromX; i < toX; i++) {
-      var exponent = Math.pow(peak.x - i, 2) / squareSigmaX +
-                Math.pow(peak.y - j, 2) / squareSigmaY;
-      var result = 10000 * peak.z * Math.exp(-exponent);
-      matrix[j][i] += result;
+  let squareSigmaX = peak.widthX * peak.widthX;
+  let squareSigmaY = peak.widthY * peak.widthY;
+  for (let j = fromY; j < toY; j++) {
+    for (let i = fromX; i < toX; i++) {
+      let exponent =
+        Math.pow(peak.x - i, 2) / squareSigmaX +
+        Math.pow(peak.y - j, 2) / squareSigmaY;
+      let result = 10000 * peak.z * Math.exp(-exponent);
+      // matrix[j][i] += result;
+      matrix.set(j, i, matrix.get(j, i) + result);
     }
   }
 }

src/SpinSystem.js

@@ -17,21 +17,22 @@ import Matrix from 'ml-matrix';
   static fromSpinusPrediction(result) {
-    var lines = result.split('\n');
-    var nspins = lines.length - 1;
-    var cs = new Array(nspins);
-    var integrals = new Array(nspins);
-    var ids = {};
-    var jc = Matrix.zeros(nspins, nspins);
+    let lines = result.split('\n');
+    let nspins = lines.length - 1;
+    let cs = new Array(nspins);
+    let integrals = new Array(nspins);
+    let ids = {};
+    let jc = Matrix.zeros(nspins, nspins);
     for (let i = 0; i < nspins; i++) {
       var tokens = lines[i].split('\t');
-      cs[i] = +tokens[2];
+      cs[i] = Number(tokens[2]);
       ids[tokens[0] - 1] = i;
-      integrals[i] = +tokens[5];// Is it always 1??
+      integrals[i] = Number(tokens[5]); // Is it always 1??
     }
     for (let i = 0; i < nspins; i++) {
       tokens = lines[i].split('\t');
-      var nCoup = (tokens.length - 4) / 3;
+      let nCoup = (tokens.length - 4) / 3;
       for (j = 0; j < nCoup; j++) {
-        var withID = tokens[4 + 3 * j] - 1;
-        var idx = ids[withID];
-        jc[i][idx] = +tokens[6 + 3 * j];
+        let withID = tokens[4 + 3 * j] - 1;
+        let idx = ids[withID];
+        // jc[i][idx] = +tokens[6 + 3 * j];
+        jc.set(i, idx, Number(tokens[6 + 3 * j]));
       }
@@ -41,4 +42,4 @@ }
     for (var j = 0; j < nspins; j++) {
-      for (var i = j; i < nspins; i++) {
-        jc[j][i] = jc[i][j];
+      for (let i = j; i < nspins; i++) {
+        jc.set(j, i, jc.get(i, j));
       }
@@ -56,3 +57,3 @@ }
     const ids = {};
-    var i, k, j;
+    let i, k, j;
     for (i = 0; i < nSpins; i++) {
@@ -66,12 +67,20 @@ cs[i] = predictions[i].delta;
       for (k = 0; k < j.length; k++) {
-        jc[ids[predictions[i].atomIDs[0]]][ids[j[k].assignment]] = j[k].coupling;
-        jc[ids[j[k].assignment]][ids[predictions[i].atomIDs[0]]] = j[k].coupling;
+        // jc[ids[predictions[i].atomIDs[0]]][ids[j[k].assignment]] = j[k].coupling;
+        // jc[ids[j[k].assignment]][ids[predictions[i].atomIDs[0]]] = j[k].coupling;
+        jc.set(
+          ids[predictions[i].atomIDs[0]],
+          ids[j[k].assignment],
+          j[k].coupling,
+        );
+        jc.set(
+          ids[j[k].assignment],
+          ids[predictions[i].atomIDs[0]],
+          j[k].coupling,
+        );
       }
       multiplicity[i] = predictions[i].integral + 1;
     }
-
     return new SpinSystem(cs, jc, multiplicity);
   }
 
-
   static ungroupAtoms(prediction) {
@@ -109,15 +118,17 @@ let result = [];
 
-
   _initClusters() {
-    this.clusters = simpleClustering(this.connectivity, { out: 'indexes' });
+    this.clusters = simpleClustering(this.connectivity.to2DArray(), {
+      out: 'indexes',
+    });
   }
 
+  //I am asumming that couplingConstants is a square matrix
   _initConnectivity() {
     const couplings = this.couplingConstants;
-    const connectivity = Matrix.ones(couplings.length, couplings.length);
-    for (var i = 0; i < couplings.length; i++) {
-      for (var j = i; j < couplings[i].length; j++) {
-        if (couplings[i][j] === 0) {
-          connectivity[i][j] = 0;
-          connectivity[j][i] = 0;
+    const connectivity = Matrix.ones(couplings.rows, couplings.rows);
+    for (let i = 0; i < couplings.rows; i++) {
+      for (let j = i; j < couplings.columns; j++) {
+        if (couplings.get(i, j) === 0) {
+          connectivity.set(i, j, 0);
+          connectivity.set(j, i, 0);
         }
@@ -129,15 +140,21 @@ }
 
-
   _calculateBetas(J, frequency) {
-    var betas = Matrix.zeros(J.length, J.length);
+    let betas = Matrix.zeros(J.rows, J.rows);
     // Before clustering, we must add hidden J, we could use molecular information if available
-    var i, j;
+    let i, j;
     for (i = 0; i < J.rows; i++) {
       for (j = i; j < J.columns; j++) {
-        if ((this.chemicalShifts[i] - this.chemicalShifts[j]) !== 0) {
-          betas[i][j] = 1 - Math.abs(J[i][j] / ((this.chemicalShifts[i] - this.chemicalShifts[j]) * frequency));
-          betas[j][i] = betas[i][j];
-        } else if (!(i === j || J[i][j] !== 0)) {
-          betas[i][j] = 1;
-          betas[j][i] = 1;
+        let element = J.get(i, j);
+        if (this.chemicalShifts[i] - this.chemicalShifts[j] !== 0) {
+          let value =
+            1 -
+            Math.abs(
+              element /
+                ((this.chemicalShifts[i] - this.chemicalShifts[j]) * frequency),
+            );
+          betas.set(i, j, value);
+          betas.set(j, i, value);
+        } else if (!(i === j || element !== 0)) {
+          betas.set(i, j, 1);
+          betas.set(j, i, 1);
         }
@@ -150,16 +167,19 @@ }
   ensureClusterSize(options) {
-    var betas = this._calculateBetas(this.couplingConstants, options.frequency || 400);
-    var cluster = hlClust.agnes(betas, { isDistanceMatrix: true });
-    var list = [];
+    let betas = this._calculateBetas(
+      this.couplingConstants,
+      options.frequency || 400,
+    );
+    let cluster = hlClust.agnes(betas.to2DArray(), { isDistanceMatrix: true });
+    let list = [];
     this._splitCluster(cluster, list, options.maxClusterSize || 8, false);
-    var clusters = this._mergeClusters(list);
+    let clusters = this._mergeClusters(list);
     this.nClusters = clusters.length;
-    // console.log(clusters);
     this.clusters = new Array(clusters.length);
-    // System.out.println(this.conmatrix);
-    for (var j = 0; j < this.nClusters; j++) {
+
+    for (let j = 0; j < this.nClusters; j++) {
       this.clusters[j] = [];
-      for (var i = 0; i < this.nSpins; i++) {
-        if (clusters[j][i] !== 0) {
-          if (clusters[j][i] < 0) {
+      for (let i = 0; i < this.nSpins; i++) {
+        let element = clusters[j][i];
+        if (element !== 0) {
+          if (element < 0) {
             this.clusters[j].push(-(i + 1));
@@ -175,8 +195,8 @@ } else {
   /**
-     * Recursively split the clusters until the maxClusterSize criteria has been ensured.
-     * @param {Array} cluster
-     * @param {Array} clusterList
-     * @param {number} maxClusterSize
-     * @param  {boolean} force
-     */
+   * Recursively split the clusters until the maxClusterSize criteria has been ensured.
+   * @param {Array} cluster
+   * @param {Array} clusterList
+   * @param {number} maxClusterSize
+   * @param  {boolean} force
+   */
   _splitCluster(cluster, clusterList, maxClusterSize, force) {
@@ -186,12 +206,12 @@ if (!force && cluster.index.length <= maxClusterSize) {
     } else {
-      for (var child of cluster.children) {
+      for (let child of cluster.children) {
         if (!isNaN(child.index) || child.index.length <= maxClusterSize) {
-          var members = this._getMembers(child);
+          let members = this._getMembers(child);
           // Add the neighbors that shares at least 1 coupling with the given cluster
-          var count = 0;
-          for (var i = 0; i < this.nSpins; i++) {
+          let count = 0;
+          for (let i = 0; i < this.nSpins; i++) {
             if (members[i] === 1) {
               count++;
-              for (var j = 0; j < this.nSpins; j++) {
-                if (this.connectivity[i][j] === 1 && members[j] === 0) {
+              for (let j = 0; j < this.nSpins; j++) {
+                if (this.connectivity.get(i, j) === 1 && members[j] === 0) {
                   members[j] = -1;
@@ -222,10 +242,10 @@ count++;
   /**
-     * Recursively gets the cluster members
-     * @param cluster
-     * @param members
-     */
+   * Recursively gets the cluster members
+   * @param cluster
+   * @param members
+   */
 
   _getMembers(cluster) {
-    var members = new Array(this.nSpins);
-    for (var i = 0; i < this.nSpins; i++) {
+    let members = new Array(this.nSpins);
+    for (let i = 0; i < this.nSpins; i++) {
       members[i] = 0;
@@ -236,3 +256,3 @@ }
     } else {
-      for (var index of cluster.index) {
+      for (let index of cluster.index) {
         members[index.index * 1] = 1;
@@ -245,4 +265,4 @@ }
   _mergeClusters(list) {
-    var nElements = 0;
-    var clusterA, clusterB, i, j, index, common, count;
+    let nElements = 0;
+    let clusterA, clusterB, i, j, index, common, count;
     for (i = list.length - 1; i >= 0; i--) {
@@ -249,0 +269,0 @@ clusterA = list[i];

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

1147061

increased by1579.03%

Lines of code

1819

increased by30.68%

Worsened metrics

Dev dependency count

5

increased by150%

Number of medium supply chain risk alerts

1

increased byInfinity%

Dependency changes

• + Addedis-any-array@2.0.1(transitive)

• + Addedml-array-max@1.2.4(transitive)

• + Addedml-array-min@1.2.3(transitive)

• + Addedml-array-rescale@1.3.7(transitive)

• + Addedml-matrix@6.12.0(transitive)

• - Removedml-array-utils@0.3.0(transitive)

• - Removedml-matrix@2.3.0(transitive)

• Updatedml-matrix@^6.2.0