Huge News!Announcing our $40M Series B led by Abstract Ventures.Learn More →
Socket
Sign inDemoInstall
Socket
Sign inDemoInstall

nmr-simulation

Package Overview
Dependencies
Maintainers
8
Versions
27
Alerts
File Explorer

Advanced tools

License
Socket logo

Install Socket

Detect and block malicious and high-risk dependencies

Install

nmr-simulation - npm Package Compare versions

Comparing version 1.0.13 to 1.0.14

48

lib/pauli.js
'use strict';
Object.defineProperty(exports, "__esModule", {
value: true
value: true
});

@@ -15,28 +15,28 @@ exports.default = getPauli;

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 };
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 _mlSparseMatrix2.default(n, m, { initialCapacity: 20 });
for (var i = 0; i < A.length; i++) {
if (i - d >= 0 && i - d < n && i < m) {
diag.set(i - d, i, A[i]);
}
const diag = new _mlSparseMatrix2.default(n, m, { initialCapacity: 20 });
for (var i = 0; i < A.length; i++) {
if (i - d >= 0 && i - d < n && i < m) {
diag.set(i - d, i, A[i]);
}
return diag;
}
return diag;
}

@@ -47,3 +47,3 @@

function getPauli(mult) {
if (mult === 2) return pauli2;else return createPauli(mult);
if (mult === 2) return pauli2;else return createPauli(mult);
}

456

lib/simulate1D.js
'use strict';
Object.defineProperty(exports, "__esModule", {
value: true
value: true
});

@@ -31,278 +31,278 @@ exports.default = simulate1d;

function simulate1d(spinSystem, options) {
var i, j;
var {
lineWidth = 1,
nbPoints = 1024,
maxClusterSize = 10,
output = 'y',
frequency: frequencyMHz = 400,
noiseFactor = 1
} = options;
var i, j;
var {
lineWidth = 1,
nbPoints = 1024,
maxClusterSize = 10,
output = 'y',
frequency: frequencyMHz = 400,
noiseFactor = 1
} = options;
nbPoints = Number(nbPoints);
nbPoints = Number(nbPoints);
const from = options.from * frequencyMHz || 0;
const to = (options.to || 10) * frequencyMHz;
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;
}
const chemicalShifts = spinSystem.chemicalShifts.slice();
for (i = 0; i < chemicalShifts.length; i++) {
chemicalShifts[i] = chemicalShifts[i] * frequencyMHz;
}
let lineWidthPoints = nbPoints * lineWidth / Math.abs(to - from) / 2.355;
let lnPoints = lineWidthPoints * 20;
let lineWidthPoints = nbPoints * lineWidth / Math.abs(to - from) / 2.355;
let lnPoints = lineWidthPoints * 20;
const gaussianLength = lnPoints | 0;
const gaussian = new Array(gaussianLength);
const b = lnPoints / 2;
const c = lineWidthPoints * lineWidthPoints * 2;
for (i = 0; i < gaussianLength; i++) {
gaussian[i] = 1e9 * Math.exp(-((i - b) * (i - b)) / c);
}
const gaussianLength = lnPoints | 0;
const gaussian = new Array(gaussianLength);
const b = lnPoints / 2;
const c = lineWidthPoints * lineWidthPoints * 2;
for (i = 0; i < gaussianLength; i++) {
gaussian[i] = 1e9 * Math.exp(-((i - b) * (i - b)) / c);
}
var result = options.withNoise ? [...new Array(nbPoints)].map(() => Math.random() * noiseFactor) : new Array(nbPoints).fill(0);
var 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++) {
const cluster = spinSystem.clusters[h];
const multiplicity = spinSystem.multiplicity;
for (var h = 0; h < spinSystem.clusters.length; h++) {
const cluster = spinSystem.clusters[h];
var clusterFake = new Array(cluster.length);
for (i = 0; i < cluster.length; i++) {
clusterFake[i] = cluster[i] < 0 ? -cluster[i] - 1 : cluster[i];
var clusterFake = new Array(cluster.length);
for (i = 0; i < cluster.length; i++) {
clusterFake[i] = cluster[i] < 0 ? -cluster[i] - 1 : cluster[i];
}
var 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.
var 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[index][clusterFake[i]] / 2;
currentSize = frequencies.length;
for (j = 0; j < currentSize; j++) {
frequencies.push(frequencies[j] + jc);
frequencies[j] -= jc;
}
}
}
var 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.
var 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[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;
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);
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 _mlMatrix2.default.DC.EVD(hamiltonian);
const V = evd.eigenvectorMatrix;
const diagB = evd.realEigenvalues;
const assignmentMatrix = new _mlSparseMatrix2.default(hamSize, hamSize);
const multLen = cluster.length;
weight = 0;
for (var n = 0; n < multLen; n++) {
const L = (0, _pauli2.default)(multiplicity[clusterFake[n]]);
const hamSize = hamiltonian.rows;
const evd = new _mlMatrix2.default.DC.EVD(hamiltonian);
const V = evd.eigenvectorMatrix;
const diagB = evd.realEigenvalues;
const assignmentMatrix = new _mlSparseMatrix2.default(hamSize, hamSize);
const multLen = cluster.length;
weight = 0;
for (var n = 0; n < multLen; n++) {
const L = (0, _pauli2.default)(multiplicity[clusterFake[n]]);
let temp = 1;
for (j = 0; j < n; j++) {
temp *= multiplicity[clusterFake[j]];
}
const A = _mlSparseMatrix2.default.eye(temp);
let temp = 1;
for (j = 0; j < n; j++) {
temp *= multiplicity[clusterFake[j]];
}
const A = _mlSparseMatrix2.default.eye(temp);
temp = 1;
for (j = n + 1; j < multLen; j++) {
temp *= multiplicity[clusterFake[j]];
}
const B = _mlSparseMatrix2.default.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]));
}
}
temp = 1;
for (j = n + 1; j < multLen; j++) {
temp *= multiplicity[clusterFake[j]];
}
const B = _mlSparseMatrix2.default.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 = _mlMatrix2.default.zeros(hamSize, hamSize);
assignmentMatrix.forEachNonZero((i, j, v) => {
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]);
}
}
}
return v;
});
let rhoip = _mlMatrix2.default.zeros(hamSize, hamSize);
assignmentMatrix.forEachNonZero((i, j, v) => {
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]);
}
}
}
return v;
});
let rhoip2 = rhoip.clone();
assignmentMatrix.forEachNonZero((i, j, v) => {
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]);
}
}
}
return v;
});
let rhoip2 = rhoip.clone();
assignmentMatrix.forEachNonZero((i, j, v) => {
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]);
}
}
}
return v;
});
const tV = V.transpose();
rhoip = tV.mmul(rhoip);
rhoip = new _mlSparseMatrix2.default(rhoip, { threshold: smallValue });
triuTimesAbs(rhoip, smallValue);
rhoip2 = tV.mmul(rhoip2);
rhoip2 = new _mlSparseMatrix2.default(rhoip2, { threshold: smallValue });
triuTimesAbs(rhoip2, smallValue);
// eslint-disable-next-line no-loop-func
rhoip2.forEachNonZero((i, j, v) => {
var val = rhoip.get(i, j);
val = Math.min(Math.abs(val), Math.abs(v));
val *= val;
const tV = V.transpose();
rhoip = tV.mmul(rhoip);
rhoip = new _mlSparseMatrix2.default(rhoip, { threshold: smallValue });
triuTimesAbs(rhoip, smallValue);
rhoip2 = tV.mmul(rhoip2);
rhoip2 = new _mlSparseMatrix2.default(rhoip2, { threshold: smallValue });
triuTimesAbs(rhoip2, smallValue);
rhoip2.forEachNonZero((i, j, v) => {
var val = rhoip.get(i, j);
val = Math.min(Math.abs(val), Math.abs(v));
val *= val;
sumI += val;
var valFreq = diagB[i] - diagB[j];
var insertIn = (0, _binarySearch2.default)(frequencies, valFreq, _numSort.asc);
if (insertIn < 0) {
frequencies.splice(-1 - insertIn, 0, valFreq);
intensities.splice(-1 - insertIn, 0, val);
} else {
intensities[insertIn] += val;
}
});
sumI += val;
var valFreq = diagB[i] - diagB[j];
var insertIn = (0, _binarySearch2.default)(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 / 32;
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);
});
}
const numFreq = frequencies.length;
if (numFreq > 0) {
weight = weight / sumI;
const diff = lineWidth / 32;
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');
}
if (output === 'xy') {
return { x: _getX(options.from, options.to, nbPoints), y: result };
}
if (output === 'y') {
return result;
}
throw new RangeError('wrong output option');
}
function addPeak(result, freq, height, from, to, nbPoints, gaussian) {
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++) {
index = i | 0;
if (i >= 0 && i < nbPoints) {
result[index] = result[index] + gaussian[indexLorentz] * height;
}
indexLorentz++;
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++) {
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;
});
A.forEachNonZero((i, j, v) => {
if (i > j) return 0;
if (Math.abs(v) <= val) return 0;
return v;
});
}
function getHamiltonian(chemicalShifts, couplingConstants, multiplicity, conMatrix, cluster) {
let hamSize = 1;
for (var i = 0; i < cluster.length; i++) {
hamSize *= multiplicity[cluster[i]];
let hamSize = 1;
for (var i = 0; i < cluster.length; i++) {
hamSize *= multiplicity[cluster[i]];
}
const clusterHam = new _mlSparseMatrix2.default(hamSize, hamSize);
for (var pos = 0; pos < cluster.length; pos++) {
var n = cluster[pos];
const L = (0, _pauli2.default)(multiplicity[n]);
let A1, B1;
let temp = 1;
for (let i = 0; i < pos; i++) {
temp *= multiplicity[cluster[i]];
}
A1 = _mlSparseMatrix2.default.eye(temp);
const clusterHam = new _mlSparseMatrix2.default(hamSize, hamSize);
temp = 1;
for (let i = pos + 1; i < cluster.length; i++) {
temp *= multiplicity[cluster[i]];
}
B1 = _mlSparseMatrix2.default.eye(temp);
for (var pos = 0; pos < cluster.length; pos++) {
var n = cluster[pos];
const alpha = chemicalShifts[n];
const kronProd = A1.kroneckerProduct(L.z).kroneckerProduct(B1);
clusterHam.add(kronProd.mul(alpha));
const L = (0, _pauli2.default)(multiplicity[n]);
for (var pos2 = 0; pos2 < cluster.length; pos2++) {
const k = cluster[pos2];
if (conMatrix[n][k] === 1) {
const S = (0, _pauli2.default)(multiplicity[k]);
let A1, B1;
let A2, B2;
let temp = 1;
for (let i = 0; i < pos; i++) {
temp *= multiplicity[cluster[i]];
for (let i = 0; i < pos2; i++) {
temp *= multiplicity[cluster[i]];
}
A1 = _mlSparseMatrix2.default.eye(temp);
A2 = _mlSparseMatrix2.default.eye(temp);
temp = 1;
for (let i = pos + 1; i < cluster.length; i++) {
temp *= multiplicity[cluster[i]];
for (let i = pos2 + 1; i < cluster.length; i++) {
temp *= multiplicity[cluster[i]];
}
B1 = _mlSparseMatrix2.default.eye(temp);
B2 = _mlSparseMatrix2.default.eye(temp);
const alpha = chemicalShifts[n];
const kronProd = A1.kroneckerProduct(L.z).kroneckerProduct(B1);
clusterHam.add(kronProd.mul(alpha));
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)));
for (var pos2 = 0; pos2 < cluster.length; pos2++) {
const k = cluster[pos2];
if (conMatrix[n][k] === 1) {
const S = (0, _pauli2.default)(multiplicity[k]);
let A2, B2;
let temp = 1;
for (let i = 0; i < pos2; i++) {
temp *= multiplicity[cluster[i]];
}
A2 = _mlSparseMatrix2.default.eye(temp);
temp = 1;
for (let i = pos2 + 1; i < cluster.length; i++) {
temp *= multiplicity[cluster[i]];
}
B2 = _mlSparseMatrix2.default.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[n][k] / 2));
}
}
clusterHam.add(kron1.mul(couplingConstants[n][k] / 2));
}
}
}
return clusterHam;
return clusterHam;
}
function _getX(from, to, nbPoints) {
const x = new Array(nbPoints);
const dx = (to - from) / (nbPoints - 1);
for (var i = 0; i < nbPoints; i++) {
x[i] = from + i * dx;
}
return x;
const x = new Array(nbPoints);
const dx = (to - from) / (nbPoints - 1);
for (var i = 0; i < nbPoints; i++) {
x[i] = from + i * dx;
}
return x;
}

136

lib/simulate2D.js
'use strict';
Object.defineProperty(exports, "__esModule", {
value: true
value: true
});

@@ -15,85 +15,85 @@ exports.default = simule2DNmrSpectrum;

let defOptions = {
H: { frequency: 400, lineWidth: 10 },
C: { frequency: 100, lineWidth: 10 }
H: { frequency: 400, lineWidth: 10 },
C: { frequency: 100, lineWidth: 10 }
};
function simule2DNmrSpectrum(table, options) {
var 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;
var lineWidthX = options.lineWidthX || defOptions[fromLabel].lineWidth;
var lineWidthY = options.lineWidthY || defOptions[toLabel].lineWidth;
var 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;
var lineWidthX = options.lineWidthX || defOptions[fromLabel].lineWidth;
var lineWidthY = options.lineWidthY || defOptions[toLabel].lineWidth;
var sigmaX = lineWidthX / frequencyX;
var sigmaY = lineWidthY / frequencyY;
var sigmaX = lineWidthX / frequencyX;
var sigmaY = lineWidthY / frequencyY;
var minX = table[0].fromChemicalShift;
var maxX = table[0].fromChemicalShift;
var minY = table[0].toChemicalShift;
var 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++;
}
var minX = table[0].fromChemicalShift;
var maxX = table[0].fromChemicalShift;
var minY = table[0].toChemicalShift;
var 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;
}
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;
}
var nbPointsX = options.nbPointsX || 512;
var nbPointsY = options.nbPointsY || 512;
var nbPointsX = options.nbPointsX || 512;
var nbPointsY = options.nbPointsY || 512;
var spectraMatrix = new _mlMatrix2.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(spectraMatrix, peak);
i++;
}
return spectraMatrix;
var spectraMatrix = new _mlMatrix2.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(spectraMatrix, peak);
i++;
}
return spectraMatrix;
}
function unitsToArrayPoints(x, from, to, nbPoints) {
return ((x - from) * nbPoints - 1) / (to - from);
return ((x - from) * nbPoints - 1) / (to - from);
}
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));
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));
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;
}
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;
}
}
}

462

lib/SpinSystem.js
'use strict';
Object.defineProperty(exports, "__esModule", {
value: true
value: true
});

@@ -26,160 +26,160 @@

class SpinSystem {
constructor(chemicalShifts, couplingConstants, multiplicity) {
this.chemicalShifts = chemicalShifts;
this.couplingConstants = couplingConstants;
this.multiplicity = multiplicity;
this.nSpins = chemicalShifts.length;
this._initConnectivity();
this._initClusters();
constructor(chemicalShifts, couplingConstants, multiplicity) {
this.chemicalShifts = chemicalShifts;
this.couplingConstants = couplingConstants;
this.multiplicity = multiplicity;
this.nSpins = chemicalShifts.length;
this._initConnectivity();
this._initClusters();
}
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 = _mlMatrix2.default.zeros(nspins, nspins);
for (let i = 0; i < nspins; i++) {
var tokens = lines[i].split('\t');
cs[i] = +tokens[2];
ids[tokens[0] - 1] = i;
integrals[i] = +tokens[5]; // Is it always 1??
}
for (let i = 0; i < nspins; i++) {
tokens = lines[i].split('\t');
var 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];
}
}
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 = _mlMatrix2.default.zeros(nspins, nspins);
for (let i = 0; i < nspins; i++) {
var tokens = lines[i].split('\t');
cs[i] = +tokens[2];
ids[tokens[0] - 1] = i;
integrals[i] = +tokens[5]; //Is it always 1??
}
for (let i = 0; i < nspins; i++) {
tokens = lines[i].split('\t');
var 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];
}
}
for (var j = 0; j < nspins; j++) {
for (var i = j; i < nspins; i++) {
jc[j][i] = jc[i][j];
}
}
return new SpinSystem(cs, jc, (0, _newArray2.default)(nspins, 2));
}
for (var j = 0; j < nspins; j++) {
for (var i = j; i < nspins; i++) {
jc[j][i] = jc[i][j];
}
}
return new SpinSystem(cs, jc, (0, _newArray2.default)(nspins, 2));
static fromPrediction(input) {
let predictions = SpinSystem.ungroupAtoms(input);
const nSpins = predictions.length;
const cs = new Array(nSpins);
const jc = _mlMatrix2.default.zeros(nSpins, nSpins);
const multiplicity = new Array(nSpins);
const ids = {};
var 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;
}
multiplicity[i] = predictions[i].integral + 1;
}
static fromPrediction(input) {
let predictions = SpinSystem.ungroupAtoms(input);
const nSpins = predictions.length;
const cs = new Array(nSpins);
const jc = _mlMatrix2.default.zeros(nSpins, nSpins);
const multiplicity = new Array(nSpins);
const ids = {};
var 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;
}
multiplicity[i] = predictions[i].integral + 1;
}
return new SpinSystem(cs, jc, multiplicity);
}
return new SpinSystem(cs, jc, multiplicity);
}
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);
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;
}
return result;
}
_initClusters() {
this.clusters = (0, _mlSimpleClustering2.default)(this.connectivity, { out: 'indexes' });
}
_initClusters() {
this.clusters = (0, _mlSimpleClustering2.default)(this.connectivity, { out: 'indexes' });
}
_initConnectivity() {
const couplings = this.couplingConstants;
const connectivity = _mlMatrix2.default.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;
}
}
_initConnectivity() {
const couplings = this.couplingConstants;
const connectivity = _mlMatrix2.default.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;
}
this.connectivity = connectivity;
}
}
this.connectivity = connectivity;
}
_calculateBetas(J, frequency) {
var betas = _mlMatrix2.default.zeros(J.length, J.length);
//Before clustering, we must add hidden J, we could use molecular information if available
var 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;
}
}
_calculateBetas(J, frequency) {
var betas = _mlMatrix2.default.zeros(J.length, J.length);
// Before clustering, we must add hidden J, we could use molecular information if available
var 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;
}
return betas;
}
}
return betas;
}
ensureClusterSize(options) {
var betas = this._calculateBetas(this.couplingConstants, options.frequency || 400);
var cluster = _mlHclust2.default.agnes(betas, { isDistanceMatrix: true });
var list = [];
this._splitCluster(cluster, list, options.maxClusterSize || 8, false);
var 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++) {
this.clusters[j] = [];
for (var i = 0; i < this.nSpins; i++) {
if (clusters[j][i] !== 0) {
if (clusters[j][i] < 0) {
this.clusters[j].push(-(i + 1));
} else {
this.clusters[j].push(i);
}
}
}
ensureClusterSize(options) {
var betas = this._calculateBetas(this.couplingConstants, options.frequency || 400);
var cluster = _mlHclust2.default.agnes(betas, { isDistanceMatrix: true });
var list = [];
this._splitCluster(cluster, list, options.maxClusterSize || 8, false);
var 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++) {
this.clusters[j] = [];
for (var i = 0; i < this.nSpins; i++) {
if (clusters[j][i] !== 0) {
if (clusters[j][i] < 0) {
this.clusters[j].push(-(i + 1));
} else {
this.clusters[j].push(i);
}
}
}
}
}
/**
/**
* Recursively split the clusters until the maxClusterSize criteria has been ensured.

@@ -191,41 +191,41 @@ * @param {Array} cluster

*/
_splitCluster(cluster, clusterList, maxClusterSize, force) {
if (!force && cluster.index.length <= maxClusterSize) {
clusterList.push(this._getMembers(cluster));
} else {
for (var child of cluster.children) {
if (!isNaN(child.index) || child.index.length <= maxClusterSize) {
var 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++) {
if (members[i] === 1) {
count++;
for (var j = 0; j < this.nSpins; j++) {
if (this.connectivity[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);
_splitCluster(cluster, clusterList, maxClusterSize, force) {
if (!force && cluster.index.length <= maxClusterSize) {
clusterList.push(this._getMembers(cluster));
} else {
for (var child of cluster.children) {
if (!isNaN(child.index) || child.index.length <= maxClusterSize) {
var 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++) {
if (members[i] === 1) {
count++;
for (var j = 0; j < this.nSpins; j++) {
if (this.connectivity[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

@@ -236,69 +236,69 @@ * @param cluster

_getMembers(cluster) {
var members = new Array(this.nSpins);
for (var i = 0; i < this.nSpins; i++) {
members[i] = 0;
}
if (!isNaN(cluster.index)) {
members[cluster.index * 1] = 1;
} else {
for (var index of cluster.index) {
members[index.index * 1] = 1;
}
}
return members;
_getMembers(cluster) {
var members = new Array(this.nSpins);
for (var i = 0; i < this.nSpins; i++) {
members[i] = 0;
}
if (!isNaN(cluster.index)) {
members[cluster.index * 1] = 1;
} else {
for (var index of cluster.index) {
members[index.index * 1] = 1;
}
}
return members;
}
_mergeClusters(list) {
var nElements = 0;
var clusterA, clusterB, i, j, index, common, count;
for (i = list.length - 1; i >= 0; i--) {
clusterA = list[i];
nElements = clusterA.length;
index = 0;
_mergeClusters(list) {
var nElements = 0;
var 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);
// 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 (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++;
}
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;
}
return list;
}
}
exports.default = SpinSystem;

8

package.json
{
"name": "nmr-simulation",
"version": "1.0.13",
"version": "1.0.14",
"description": "Simulate NMR spectra from spin systems",

@@ -35,7 +35,9 @@ "main": "lib/index.js",

"ml-sparse-matrix": "^0.2.1",
"num-sort": "^1.0.0"
"new-array": "^1.0.0",
"num-sort": "^1.0.0",
"should": "^13.2.3"
},
"devDependencies": {
"nmr-predictor": "^1.1.6"
"nmr-predictor": "^1.1.7"
}
}

66

README.md
# nmr-simulation
[![NPM version][npm-image]][npm-url]
[![build status][travis-image]][travis-url]
[![David deps][david-image]][david-url]
[![npm download][download-image]][download-url]
Simulate NMR spectra from spin systems
The computational cost for the simulation of NMR spectra grows exponentially with the number of nuclei. Today, the memory available to store the Hamiltonian limits the size of the system that can be studied. Modern computers enable to tackle systems containing up to 13 spins [1], which obviously does not allow to study most molecules of interest in research. This issue can be addressed by identifying groups of spins or fragments that are not or only weakly interacting together, i.e., that only share weakly coupled spin pairs. Such a fragmentation is only permitted in the weak coupling regime, i.e., when the coupling interaction is weak compared to the difference in chemical shift of the coupled spins. Here, we propose a procedure that removes weak coupling interactions in order to split the spin system efficiently and to correct a posteriori for the effect of the neglected couplings. This approach yields accurate spectra when the adequate interactions are removed, i.e., between spins only involved in weak coupling interactions, but fails otherwise. As a result, the computational time for the simulation of 1D spectra grows linearly with the size of the spin system.
This library is a JavaScript implementation of such algorithm. The spin system is specified as a set of magnetically equivalent spins. The description of the input format is ilustrated by means of the bellow example. There we have a spin system composed of 4 spins(11, 12, 5 and 6). Spins 11 and 12 are chemically equivalent, so they must be specified on the same group. nbAtoms stands for the number of atoms in the group. It should math the atomIDs.length. AtomLabel stands for the kind of atom of the group. It could be H, C, N, P refering to nuclei 1H, 13C, 15N and 31P respectively. Coupling between atoms should be specified in the j array.
## Installation
`$ npm i nmr-simulation`
## [API Documentation](https://mljs.github.io/spectra/packages/nmr-simulation)
## Example
```js
const simulation = require('nmr-simulation');
const prediction = [
{
atomIDs: ['11', '12'],
nbAtoms: 2,
delta: 1,
atomLabel: 'H'
}, {
atomIDs: ['5', '6'],
nbAtoms: 2,
delta: 1,
atomLabel: 'H'
}
];
var options1h = {
frequency: 400.082470657773,
from: 0,
to: 3,
lineWidth: 3,
nbPoints: 16384,
maxClusterSize: 8,
output: 'xy'
};
const spinSystem = nmr.SpinSystem.fromPrediction(prediction1h);
spinSystem.ensureClusterSize(options1h);
var spectrum = nmr.simulate1D(spinSystem, options1h);
console.log(spectrum.x);// x in PPM
console.log(spectrum.y);// y in arbitrary units
```
## License
[MIT](./LICENSE)
[npm-image]: https://img.shields.io/npm/v/nmr-simulation.svg?style=flat-square
[npm-url]: https://www.npmjs.com/package/nmr-simulation.
[travis-image]: https://img.shields.io/travis/mljs/nmr-simulation./master.svg?style=flat-square
[travis-url]: https://travis-ci.org/mljs/nmr-simulation
[david-image]: https://img.shields.io/david/mljs/nmr-simulation.svg?style=flat-square
[david-url]: https://david-dm.org/mljs/nmr-simulation
[download-image]: https://img.shields.io/npm/dm/nmr-simulation.svg?style=flat-square
[download-url]: https://www.npmjs.com/package/nmr-simulation

6

src/index.js

@@ -1,3 +0,3 @@

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

48

src/pauli.js
import SparseMatrix from 'ml-sparse-matrix';
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};
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 (var i = 0; i < A.length; i++) {
if ((i - d) >= 0 && (i - d) < n && i < m) {
diag.set(i - d, i, A[i]);
}
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) {
diag.set(i - d, i, A[i]);
}
return diag;
}
return diag;
}

@@ -35,4 +35,4 @@

export default function getPauli(mult) {
if (mult === 2) return pauli2;
else return createPauli(mult);
if (mult === 2) return pauli2;
else return createPauli(mult);
}

469

src/simulate1D.js
import Matrix from 'ml-matrix';
import SparseMatrix from 'ml-sparse-matrix';
import binarySearch from 'binary-search';
import {asc as sortAsc} from 'num-sort';
import { asc as sortAsc } from 'num-sort';

@@ -11,285 +11,284 @@ import getPauli from './pauli';

export default function simulate1d(spinSystem, options) {
var i, j;
var {
lineWidth = 1,
nbPoints = 1024,
maxClusterSize = 10,
output = 'y',
frequency: frequencyMHz = 400,
noiseFactor = 1
} = options;
var i, j;
var {
lineWidth = 1,
nbPoints = 1024,
maxClusterSize = 10,
output = 'y',
frequency: frequencyMHz = 400,
noiseFactor = 1
} = options;
nbPoints = Number(nbPoints);
nbPoints = Number(nbPoints);
const from = options.from * frequencyMHz || 0;
const to = (options.to || 10) * frequencyMHz;
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;
}
const chemicalShifts = spinSystem.chemicalShifts.slice();
for (i = 0; i < chemicalShifts.length; i++) {
chemicalShifts[i] = chemicalShifts[i] * frequencyMHz;
}
let lineWidthPoints = (nbPoints * lineWidth / Math.abs(to - from)) / 2.355;
let lnPoints = lineWidthPoints * 20;
let lineWidthPoints = (nbPoints * lineWidth / Math.abs(to - from)) / 2.355;
let lnPoints = lineWidthPoints * 20;
const gaussianLength = lnPoints | 0;
const gaussian = new Array(gaussianLength);
const b = lnPoints / 2;
const c = lineWidthPoints * lineWidthPoints * 2;
for (i = 0; i < gaussianLength; i++) {
gaussian[i] = 1e9 * Math.exp(-((i - b) * (i - b)) / c);
}
const gaussianLength = lnPoints | 0;
const gaussian = new Array(gaussianLength);
const b = lnPoints / 2;
const c = lineWidthPoints * lineWidthPoints * 2;
for (i = 0; i < gaussianLength; i++) {
gaussian[i] = 1e9 * Math.exp(-((i - b) * (i - b)) / c);
}
var result = options.withNoise ? [...new Array(nbPoints)].map(() => Math.random() * noiseFactor) : new Array(nbPoints).fill(0);
var 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++) {
const cluster = spinSystem.clusters[h];
const multiplicity = spinSystem.multiplicity;
for (var h = 0; h < spinSystem.clusters.length; h++) {
const cluster = spinSystem.clusters[h];
var clusterFake = new Array(cluster.length);
for (i = 0; i < cluster.length; i++) {
clusterFake[i] = cluster[i] < 0 ? -cluster[i] - 1 : cluster[i];
var clusterFake = new Array(cluster.length);
for (i = 0; i < cluster.length; i++) {
clusterFake[i] = cluster[i] < 0 ? -cluster[i] - 1 : cluster[i];
}
var 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.
var 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[index][clusterFake[i]] / 2;
currentSize = frequencies.length;
for (j = 0; j < currentSize; j++) {
frequencies.push(frequencies[j] + jc);
frequencies[j] -= jc;
}
}
}
var 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.
var 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[index][clusterFake[i]] / 2;
currentSize = frequencies.length;
for (j = 0; j < currentSize; j++) {
frequencies.push(frequencies[j] + jc);
frequencies[j] -= jc;
}
}
}
frequencies.sort(sortAsc);
sumI = frequencies.length;
weight = 1;
frequencies.sort(sortAsc);
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
);
for (i = 0; i < sumI; i++) {
intensities.push(1);
}
const hamSize = hamiltonian.rows;
const evd = new Matrix.DC.EVD(hamiltonian);
const V = evd.eigenvectorMatrix;
const diagB = evd.realEigenvalues;
const assignmentMatrix = new SparseMatrix(hamSize, hamSize);
const multLen = cluster.length;
weight = 0;
for (var 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 {
const hamiltonian = getHamiltonian(
chemicalShifts,
spinSystem.couplingConstants,
multiplicity,
spinSystem.connectivity,
clusterFake
);
assignmentMatrix.add(tempMat.mul(cluster[n]));
}
}
const hamSize = hamiltonian.rows;
const evd = new Matrix.DC.EVD(hamiltonian);
const V = evd.eigenvectorMatrix;
const diagB = evd.realEigenvalues;
const assignmentMatrix = new SparseMatrix(hamSize, hamSize);
const multLen = cluster.length;
weight = 0;
for (var n = 0; n < multLen; n++) {
const L = getPauli(multiplicity[clusterFake[n]]);
let rhoip = Matrix.zeros(hamSize, hamSize);
assignmentMatrix.forEachNonZero((i, j, v) => {
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]);
}
}
}
return v;
});
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 rhoip2 = rhoip.clone();
assignmentMatrix.forEachNonZero((i, j, v) => {
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]);
}
}
}
return v;
});
let rhoip = Matrix.zeros(hamSize, hamSize);
assignmentMatrix.forEachNonZero((i, j, v) => {
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]);
}
}
}
return v;
});
const tV = V.transpose();
rhoip = tV.mmul(rhoip);
rhoip = new SparseMatrix(rhoip, { threshold: smallValue });
triuTimesAbs(rhoip, smallValue);
rhoip2 = tV.mmul(rhoip2);
rhoip2 = new SparseMatrix(rhoip2, { threshold: smallValue });
triuTimesAbs(rhoip2, smallValue);
// eslint-disable-next-line no-loop-func
rhoip2.forEachNonZero((i, j, v) => {
var val = rhoip.get(i, j);
val = Math.min(Math.abs(val), Math.abs(v));
val *= val;
let rhoip2 = rhoip.clone();
assignmentMatrix.forEachNonZero((i, j, v) => {
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]);
}
}
}
return v;
});
const tV = V.transpose();
rhoip = tV.mmul(rhoip);
rhoip = new SparseMatrix(rhoip, {threshold: smallValue});
triuTimesAbs(rhoip, smallValue);
rhoip2 = tV.mmul(rhoip2);
rhoip2 = new SparseMatrix(rhoip2, {threshold: smallValue});
triuTimesAbs(rhoip2, smallValue);
rhoip2.forEachNonZero((i, j, v) => {
var val = rhoip.get(i, j);
val = Math.min(Math.abs(val), Math.abs(v));
val *= val;
sumI += val;
var valFreq = diagB[i] - diagB[j];
var insertIn = binarySearch(frequencies, valFreq, sortAsc);
if (insertIn < 0) {
frequencies.splice(-1 - insertIn, 0, valFreq);
intensities.splice(-1 - insertIn, 0, val);
} else {
intensities[insertIn] += val;
}
});
sumI += val;
var valFreq = diagB[i] - diagB[j];
var insertIn = binarySearch(frequencies, valFreq, sortAsc);
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 / 32;
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);
});
}
const numFreq = frequencies.length;
if (numFreq > 0) {
weight = weight / sumI;
const diff = lineWidth / 32;
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');
}
if (output === 'xy') {
return { x: _getX(options.from, options.to, nbPoints), y: result };
}
if (output === 'y') {
return result;
}
throw new RangeError('wrong output option');
}
function addPeak(result, freq, height, from, to, nbPoints, gaussian) {
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++) {
index = i | 0;
if (i >= 0 && i < nbPoints) {
result[index] = result[index] + gaussian[indexLorentz] * height;
}
indexLorentz++;
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++) {
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;
});
A.forEachNonZero((i, j, v) => {
if (i > j) return 0;
if (Math.abs(v) <= val) return 0;
return v;
});
}
function getHamiltonian(chemicalShifts, couplingConstants, multiplicity, conMatrix, cluster) {
let hamSize = 1;
for (var i = 0; i < cluster.length; i++) {
hamSize *= multiplicity[cluster[i]];
let hamSize = 1;
for (var i = 0; i < cluster.length; i++) {
hamSize *= multiplicity[cluster[i]];
}
const clusterHam = new SparseMatrix(hamSize, hamSize);
for (var pos = 0; pos < cluster.length; pos++) {
var 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);
const clusterHam = new SparseMatrix(hamSize, hamSize);
temp = 1;
for (let i = pos + 1; i < cluster.length; i++) {
temp *= multiplicity[cluster[i]];
}
B1 = SparseMatrix.eye(temp);
for (var pos = 0; pos < cluster.length; pos++) {
var n = cluster[pos];
const alpha = chemicalShifts[n];
const kronProd = A1.kroneckerProduct(L.z).kroneckerProduct(B1);
clusterHam.add(kronProd.mul(alpha));
const L = getPauli(multiplicity[n]);
for (var pos2 = 0; pos2 < cluster.length; pos2++) {
const k = cluster[pos2];
if (conMatrix[n][k] === 1) {
const S = getPauli(multiplicity[k]);
let A1, B1;
let A2, B2;
let temp = 1;
for (let i = 0; i < pos; i++) {
temp *= multiplicity[cluster[i]];
for (let i = 0; i < pos2; i++) {
temp *= multiplicity[cluster[i]];
}
A1 = SparseMatrix.eye(temp);
A2 = SparseMatrix.eye(temp);
temp = 1;
for (let i = pos + 1; i < cluster.length; i++) {
temp *= multiplicity[cluster[i]];
for (let i = pos2 + 1; i < cluster.length; i++) {
temp *= multiplicity[cluster[i]];
}
B1 = SparseMatrix.eye(temp);
B2 = SparseMatrix.eye(temp);
const alpha = chemicalShifts[n];
const kronProd = A1.kroneckerProduct(L.z).kroneckerProduct(B1);
clusterHam.add(kronProd.mul(alpha));
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)));
for (var pos2 = 0; pos2 < cluster.length; pos2++) {
const k = cluster[pos2];
if (conMatrix[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[n][k] / 2));
}
}
clusterHam.add(kron1.mul(couplingConstants[n][k] / 2));
}
}
}
return clusterHam;
return clusterHam;
}
function _getX(from, to, nbPoints) {
const x = new Array(nbPoints);
const dx = (to - from) / (nbPoints - 1);
for (var i = 0; i < nbPoints; i++) {
x[i] = from + i * dx;
}
return x;
const x = new Array(nbPoints);
const dx = (to - from) / (nbPoints - 1);
for (var i = 0; i < nbPoints; i++) {
x[i] = from + i * dx;
}
return x;
}

134

src/simulate2D.js
import Matrix from 'ml-matrix';
let defOptions = {
H: {frequency: 400, lineWidth: 10},
C: {frequency: 100, lineWidth: 10}
H: { frequency: 400, lineWidth: 10 },
C: { frequency: 100, lineWidth: 10 }
};
export default function simule2DNmrSpectrum(table, options) {
var 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;
var lineWidthX = options.lineWidthX || defOptions[fromLabel].lineWidth;
var lineWidthY = options.lineWidthY || defOptions[toLabel].lineWidth;
var 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;
var lineWidthX = options.lineWidthX || defOptions[fromLabel].lineWidth;
var lineWidthY = options.lineWidthY || defOptions[toLabel].lineWidth;
var sigmaX = lineWidthX / frequencyX;
var sigmaY = lineWidthY / frequencyY;
var sigmaX = lineWidthX / frequencyX;
var sigmaY = lineWidthY / frequencyY;
var minX = table[0].fromChemicalShift;
var maxX = table[0].fromChemicalShift;
var minY = table[0].toChemicalShift;
var 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++;
}
var minX = table[0].fromChemicalShift;
var maxX = table[0].fromChemicalShift;
var minY = table[0].toChemicalShift;
var 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;
}
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;
}
var nbPointsX = options.nbPointsX || 512;
var nbPointsY = options.nbPointsY || 512;
var nbPointsX = options.nbPointsX || 512;
var nbPointsY = options.nbPointsY || 512;
var spectraMatrix = new Matrix(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(spectraMatrix, peak);
i++;
}
return spectraMatrix;
var spectraMatrix = new Matrix(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(spectraMatrix, peak);
i++;
}
return spectraMatrix;
}
function unitsToArrayPoints(x, from, to, nbPoints) {
return ((x - from) * nbPoints - 1) / (to - from);
return ((x - from) * nbPoints - 1) / (to - from);
}
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));
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));
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 +
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;
}
var result = 10000 * peak.z * Math.exp(-exponent);
matrix[j][i] += result;
}
}
}

460

src/SpinSystem.js

@@ -7,163 +7,163 @@ import Matrix from 'ml-matrix';

export default class SpinSystem {
constructor(chemicalShifts, couplingConstants, multiplicity) {
this.chemicalShifts = chemicalShifts;
this.couplingConstants = couplingConstants;
this.multiplicity = multiplicity;
this.nSpins = chemicalShifts.length;
this._initConnectivity();
this._initClusters();
constructor(chemicalShifts, couplingConstants, multiplicity) {
this.chemicalShifts = chemicalShifts;
this.couplingConstants = couplingConstants;
this.multiplicity = multiplicity;
this.nSpins = chemicalShifts.length;
this._initConnectivity();
this._initClusters();
}
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);
for (let i = 0; i < nspins; i++) {
var tokens = lines[i].split('\t');
cs[i] = +tokens[2];
ids[tokens[0] - 1] = i;
integrals[i] = +tokens[5];// Is it always 1??
}
for (let i = 0; i < nspins; i++) {
tokens = lines[i].split('\t');
var 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];
}
}
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);
for (let i = 0; i < nspins; i++) {
var tokens = lines[i].split('\t');
cs[i] = +tokens[2];
ids[tokens[0] - 1] = i;
integrals[i] = +tokens[5];//Is it always 1??
}
for (let i = 0; i < nspins; i++) {
tokens = lines[i].split('\t');
var 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];
}
}
for (var j = 0; j < nspins; j++) {
for (var i = j; i < nspins; i++) {
jc[j][i] = jc[i][j];
}
}
return new SpinSystem(cs, jc, newArray(nspins, 2));
}
for (var j = 0; j < nspins; j++) {
for (var i = j; i < nspins; i++) {
jc[j][i] = jc[i][j];
}
}
return new SpinSystem(cs, jc, newArray(nspins, 2));
static fromPrediction(input) {
let predictions = SpinSystem.ungroupAtoms(input);
const nSpins = predictions.length;
const cs = new Array(nSpins);
const jc = Matrix.zeros(nSpins, nSpins);
const multiplicity = new Array(nSpins);
const ids = {};
var 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;
}
multiplicity[i] = predictions[i].integral + 1;
}
static fromPrediction(input) {
let predictions = SpinSystem.ungroupAtoms(input);
const nSpins = predictions.length;
const cs = new Array(nSpins);
const jc = Matrix.zeros(nSpins, nSpins);
const multiplicity = new Array(nSpins);
const ids = {};
var 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;
}
multiplicity[i] = predictions[i].integral + 1;
}
return new SpinSystem(cs, jc, multiplicity);
}
return new SpinSystem(cs, jc, multiplicity);
}
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);
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;
}
return result;
}
_initClusters() {
this.clusters = simpleClustering(this.connectivity, {out: 'indexes'});
}
_initClusters() {
this.clusters = simpleClustering(this.connectivity, { out: 'indexes' });
}
_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;
}
}
_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;
}
this.connectivity = connectivity;
}
}
this.connectivity = connectivity;
}
_calculateBetas(J, frequency) {
var betas = Matrix.zeros(J.length, J.length);
//Before clustering, we must add hidden J, we could use molecular information if available
var 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;
}
}
_calculateBetas(J, frequency) {
var betas = Matrix.zeros(J.length, J.length);
// Before clustering, we must add hidden J, we could use molecular information if available
var 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;
}
return betas;
}
}
return betas;
}
ensureClusterSize(options) {
var betas = this._calculateBetas(this.couplingConstants, options.frequency || 400);
var cluster = hlClust.agnes(betas, {isDistanceMatrix: true});
var list = [];
this._splitCluster(cluster, list, options.maxClusterSize || 8, false);
var 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++) {
this.clusters[j] = [];
for (var i = 0; i < this.nSpins; i++) {
if (clusters[j][i] !== 0) {
if (clusters[j][i] < 0) {
this.clusters[j].push(-(i + 1));
} else {
this.clusters[j].push(i);
}
}
}
ensureClusterSize(options) {
var betas = this._calculateBetas(this.couplingConstants, options.frequency || 400);
var cluster = hlClust.agnes(betas, { isDistanceMatrix: true });
var list = [];
this._splitCluster(cluster, list, options.maxClusterSize || 8, false);
var 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++) {
this.clusters[j] = [];
for (var i = 0; i < this.nSpins; i++) {
if (clusters[j][i] !== 0) {
if (clusters[j][i] < 0) {
this.clusters[j].push(-(i + 1));
} else {
this.clusters[j].push(i);
}
}
}
}
}
/**
/**
* Recursively split the clusters until the maxClusterSize criteria has been ensured.

@@ -175,41 +175,41 @@ * @param {Array} cluster

*/
_splitCluster(cluster, clusterList, maxClusterSize, force) {
if (!force && cluster.index.length <= maxClusterSize) {
clusterList.push(this._getMembers(cluster));
} else {
for (var child of cluster.children) {
if (!isNaN(child.index) || child.index.length <= maxClusterSize) {
var 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++) {
if (members[i] === 1) {
count++;
for (var j = 0; j < this.nSpins; j++) {
if (this.connectivity[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);
_splitCluster(cluster, clusterList, maxClusterSize, force) {
if (!force && cluster.index.length <= maxClusterSize) {
clusterList.push(this._getMembers(cluster));
} else {
for (var child of cluster.children) {
if (!isNaN(child.index) || child.index.length <= maxClusterSize) {
var 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++) {
if (members[i] === 1) {
count++;
for (var j = 0; j < this.nSpins; j++) {
if (this.connectivity[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

@@ -220,68 +220,68 @@ * @param cluster

_getMembers(cluster) {
var members = new Array(this.nSpins);
for (var i = 0; i < this.nSpins; i++) {
members[i] = 0;
}
if (!isNaN(cluster.index)) {
members[cluster.index * 1] = 1;
} else {
for (var index of cluster.index) {
members[index.index * 1] = 1;
}
}
return members;
_getMembers(cluster) {
var members = new Array(this.nSpins);
for (var i = 0; i < this.nSpins; i++) {
members[i] = 0;
}
if (!isNaN(cluster.index)) {
members[cluster.index * 1] = 1;
} else {
for (var index of cluster.index) {
members[index.index * 1] = 1;
}
}
return members;
}
_mergeClusters(list) {
var nElements = 0;
var clusterA, clusterB, i, j, index, common, count;
for (i = list.length - 1; i >= 0; i--) {
clusterA = list[i];
nElements = clusterA.length;
index = 0;
_mergeClusters(list) {
var nElements = 0;
var 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);
// 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 (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++;
}
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;
}
return list;
}
}
SocketSocket SOC 2 Logo

Product

  • Package Alerts
  • Integrations
  • Docs
  • Pricing
  • FAQ
  • Roadmap
  • Changelog

About

Packages

npm

Go

Maven

PyPI

Rubygems

Stay in touch

Get open source security insights delivered straight into your inbox.

  • Terms
  • Privacy
  • Security

Made with ⚡️ by Socket Inc