From 5cf330c68cc637ae8f6bb1308d9a6022059d2223 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Sun, 26 Jul 2026 21:42:40 +0200 Subject: [PATCH 1/4] test: add benchmarks for distances and similarities Built on benchmark.js, the library the other mljs repositories use: - all.js lists the cost of every exported function - beforeAfter.js compares two implementations of the same function, both copied into one process so they see identical data - exponentiation.js isolates the cost of the `**` operator - elementReads.js isolates the cost of repeated `a[i]` reads - equivalence.js checks `x ** 1.5` against `x * Math.sqrt(x)` Baseline cost of the current implementations, measured with beforeAfter.js on node 24 (n = 10000, Float64Array). These are per-element costs of each function, not speedups: topsoe 21.99 ns/element similarity.pearson 19.06 ns/element kumarJohnson 15.50 ns/element minkowski(p=1) 11.51 ns/element squared 4.51 ns/element cosine 4.38 ns/element pearson 3.74 ns/element clark 3.22 ns/element The next commit optimizes these and reports the same measurements again. Assisted-By: Claude Opus 5 (1M context) --- benchmark/all.js | 84 ++++++++ benchmark/beforeAfter.js | 373 ++++++++++++++++++++++++++++++++++++ benchmark/elementReads.js | 153 +++++++++++++++ benchmark/equivalence.js | 137 +++++++++++++ benchmark/exponentiation.js | 184 ++++++++++++++++++ package.json | 1 + 6 files changed, 932 insertions(+) create mode 100644 benchmark/all.js create mode 100644 benchmark/beforeAfter.js create mode 100644 benchmark/elementReads.js create mode 100644 benchmark/equivalence.js create mode 100644 benchmark/exponentiation.js diff --git a/benchmark/all.js b/benchmark/all.js new file mode 100644 index 0000000..2402f43 --- /dev/null +++ b/benchmark/all.js @@ -0,0 +1,84 @@ +/* + * Benchmarks every distance and similarity function so their relative cost can + * be compared. Use beforeAfter.js instead to compare two implementations of the + * same function. + * + * Run with: node benchmark/all.js (or: bun benchmark/all.js) + */ +import Benchmark from 'benchmark'; + +import { distance, similarity } from '../src/index.ts'; + +const LENGTH = 10000; + +function makeVector(seed) { + const vector = new Float64Array(LENGTH); + let state = seed; + for (let i = 0; i < LENGTH; i++) { + state = (state * 1103515245 + 12345) % 2147483648; + vector[i] = 0.1 + (state / 2147483648) * 2; + } + return vector; +} + +const a = makeVector(42); +const b = makeVector(1337); + +const entries = []; +for (const [namespace, functions] of [ + ['distance', distance], + ['similarity', similarity], +]) { + for (const [name, callback] of Object.entries(functions)) { + if (typeof callback !== 'function') continue; + if (name === 'minkowski') { + for (const p of [1, 2, 3]) { + entries.push([ + `${namespace}.${name}(p=${p})`, + (x, y) => callback(x, y, p), + ]); + } + continue; + } + entries.push([`${namespace}.${name}`, callback]); + } +} + +function log(message) { + // eslint-disable-next-line no-console -- benchmark output + console.log(message); +} + +const results = []; +const suite = new Benchmark.Suite(); +for (const [name, callback] of entries) { + suite.add(name, () => callback(a, b), { minSamples: 30 }); +} + +suite + .on('cycle', (event) => { + const { name, hz, stats } = event.target; + results.push({ + name, + nanoseconds: 1e9 / hz, + rme: stats.rme, + samples: stats.sample.length, + }); + }) + .on('complete', () => { + results.sort((first, second) => second.nanoseconds - first.nanoseconds); + log(`\nn = ${LENGTH}, sorted slowest first\n`); + log( + `${'function'.padEnd(32)}${'per call'.padStart(12)}${'per element'.padStart(14)}${'error'.padStart(9)} result`, + ); + const values = new Map(entries.map(([name, callback]) => [name, callback])); + for (const { name, nanoseconds, rme } of results) { + // A wide confidence interval means the engine kept re-tiering this one; + // treat the number as indicative only. + const flag = rme > 10 ? ' (!)' : ''; + log( + `${name.padEnd(32)}${`${(nanoseconds / 1000).toFixed(1)} µs`.padStart(12)}${`${(nanoseconds / LENGTH).toFixed(2)} ns`.padStart(14)}${`±${rme.toFixed(1)}%${flag}`.padStart(13)} ${values.get(name)(a, b)}`, + ); + } + }) + .run({ async: false }); diff --git a/benchmark/beforeAfter.js b/benchmark/beforeAfter.js new file mode 100644 index 0000000..f80dbab --- /dev/null +++ b/benchmark/beforeAfter.js @@ -0,0 +1,373 @@ +/* + * Compares the previous implementation of each optimized function with the + * current one. Both versions are copied in here so that they run in the same + * process on the same data. + * + * Run with: node benchmark/beforeAfter.js + */ +import Benchmark from 'benchmark'; + +const LENGTHS = [1000, 10000, 100000]; +const KINDS = ['Array', 'Float64Array']; + +function makeArrays(kind, length) { + const a = + kind === 'Float64Array' ? new Float64Array(length) : new Array(length); + const b = + kind === 'Float64Array' ? new Float64Array(length) : new Array(length); + let state = 42; + for (let i = 0; i < length; i++) { + state = (state * 1103515245 + 12345) % 2147483648; + a[i] = 0.1 + (state / 2147483648) * 2; + state = (state * 1103515245 + 12345) % 2147483648; + b[i] = 0.1 + (state / 2147483648) * 2; + } + return [a, b]; +} + +/* ------------------------------------------------ distances/pearson */ +function pearsonBefore(a, b) { + let d = 0; + for (let i = 0; i < a.length; i++) { + d += ((a[i] - b[i]) * (a[i] - b[i])) / b[i]; + } + return d; +} +function pearsonAfter(a, b) { + let d = 0; + for (let i = 0; i < a.length; i++) { + const bi = b[i]; + const diff = a[i] - bi; + d += (diff * diff) / bi; + } + return d; +} + +/* ------------------------------------------------ distances/squared */ +function squaredBefore(a, b) { + let d = 0; + for (let i = 0; i < a.length; i++) { + d += ((a[i] - b[i]) * (a[i] - b[i])) / (a[i] + b[i]); + } + return d; +} +function squaredAfter(a, b) { + let d = 0; + for (let i = 0; i < a.length; i++) { + const ai = a[i]; + const bi = b[i]; + const diff = ai - bi; + d += (diff * diff) / (ai + bi); + } + return d; +} + +/* ------------------------------------------------ similarities/cosine */ +function cosineBefore(a, b) { + let p = 0; + let p2 = 0; + let q2 = 0; + for (let i = 0; i < a.length; i++) { + p += a[i] * b[i]; + p2 += a[i] * a[i]; + q2 += b[i] * b[i]; + } + return p / (Math.sqrt(p2) * Math.sqrt(q2)); +} +function cosineAfter(a, b) { + let p = 0; + let p2 = 0; + let q2 = 0; + for (let i = 0; i < a.length; i++) { + const ai = a[i]; + const bi = b[i]; + p += ai * bi; + p2 += ai * ai; + q2 += bi * bi; + } + return p / (Math.sqrt(p2) * Math.sqrt(q2)); +} + +/* ------------------------------------------------ distances/clark */ +function clarkBefore(a, b) { + let d = 0; + for (let i = 0; i < a.length; i++) { + d += (Math.abs(a[i] - b[i]) / (a[i] + b[i])) ** 2; + } + return Math.sqrt(d); +} +function clarkAfter(a, b) { + let d = 0; + for (let i = 0; i < a.length; i++) { + const ai = a[i]; + const bi = b[i]; + const ratio = (ai - bi) / (ai + bi); + d += ratio * ratio; + } + return Math.sqrt(d); +} + +/* ------------------------------------------------ distances/kumarJohnson */ +function kumarJohnsonBefore(a, b) { + let ans = 0; + for (let i = 0; i < a.length; i++) { + ans += (a[i] * a[i] - b[i] * b[i]) ** 2 / (2 * (a[i] * b[i]) ** 1.5); + } + return ans; +} +function kumarJohnsonAfter(a, b) { + let ans = 0; + for (let i = 0; i < a.length; i++) { + const ai = a[i]; + const bi = b[i]; + const numerator = ai * ai - bi * bi; + const prod = ai * bi; + ans += (numerator * numerator) / (2 * prod * Math.sqrt(prod)); + } + return ans; +} + +/* ------------------------------------------------ distances/minkowski + * Each order gets its own copy on purpose. Calling one shared `minkowski` with + * p = 1 and p = 2 makes the exponent polymorphic and lets the two orders share + * inline caches, which hides the difference the branches are meant to measure. + */ +function minkowskiBeforeP1(a, b, p) { + let d = 0; + for (let i = 0; i < a.length; i++) { + d += Math.abs(a[i] - b[i]) ** p; + } + return d ** (1 / p); +} +function minkowskiAfterP1(a, b, p) { + let d = 0; + if (p === 1) { + for (let i = 0; i < a.length; i++) { + d += Math.abs(a[i] - b[i]); + } + return d; + } + if (p === 2) { + for (let i = 0; i < a.length; i++) { + const diff = a[i] - b[i]; + d += diff * diff; + } + return Math.sqrt(d); + } + for (let i = 0; i < a.length; i++) { + d += Math.abs(a[i] - b[i]) ** p; + } + return d ** (1 / p); +} +function minkowskiBeforeP2(a, b, p) { + let d = 0; + for (let i = 0; i < a.length; i++) { + d += Math.abs(a[i] - b[i]) ** p; + } + return d ** (1 / p); +} +function minkowskiAfterP2(a, b, p) { + let d = 0; + if (p === 1) { + for (let i = 0; i < a.length; i++) { + d += Math.abs(a[i] - b[i]); + } + return d; + } + if (p === 2) { + for (let i = 0; i < a.length; i++) { + const diff = a[i] - b[i]; + d += diff * diff; + } + return Math.sqrt(d); + } + for (let i = 0; i < a.length; i++) { + d += Math.abs(a[i] - b[i]) ** p; + } + return d ** (1 / p); +} + +/* ------------------------------------------------ similarities/pearson */ +function meanOf(input) { + let sumValue = 0; + for (const value of input) sumValue += value; + return sumValue / input.length; +} +// A private copy of cosine: sharing cosineBefore with the cosine pair would +// feed it both Float64Array and Array inputs and make its loads polymorphic. +function cosineForPearson(a, b) { + let p = 0; + let p2 = 0; + let q2 = 0; + for (let i = 0; i < a.length; i++) { + p += a[i] * b[i]; + p2 += a[i] * a[i]; + q2 += b[i] * b[i]; + } + return p / (Math.sqrt(p2) * Math.sqrt(q2)); +} +function pearsonSimilarityBefore(a, b) { + const avgA = meanOf(a); + const avgB = meanOf(b); + const newA = new Array(a.length); + const newB = new Array(b.length); + for (let i = 0; i < newA.length; i++) { + newA[i] = a[i] - avgA; + newB[i] = b[i] - avgB; + } + return cosineForPearson(newA, newB); +} +function pearsonSimilarityAfter(a, b) { + const length = a.length; + let sumA = 0; + let sumB = 0; + for (let i = 0; i < length; i++) { + sumA += a[i]; + sumB += b[i]; + } + const avgA = sumA / length; + const avgB = sumB / length; + let p = 0; + let p2 = 0; + let q2 = 0; + for (let i = 0; i < length; i++) { + const centredA = a[i] - avgA; + const centredB = b[i] - avgB; + p += centredA * centredB; + p2 += centredA * centredA; + q2 += centredB * centredB; + } + return p / (Math.sqrt(p2) * Math.sqrt(q2)); +} + +/* ------------------------------------------------ topsoe */ +function topsoeBefore(a, b) { + let ans = 0; + for (let i = 0; i < a.length; i++) { + ans += + a[i] * Math.log((2 * a[i]) / (a[i] + b[i])) + + b[i] * Math.log((2 * b[i]) / (a[i] + b[i])); + } + return ans; +} +function topsoeAfter(a, b) { + let ans = 0; + for (let i = 0; i < a.length; i++) { + const ai = a[i]; + const bi = b[i]; + const sum = ai + bi; + ans += ai * Math.log((2 * ai) / sum) + bi * Math.log((2 * bi) / sum); + } + return ans; +} + +const PAIRS = [ + { name: 'distances/pearson', before: pearsonBefore, after: pearsonAfter }, + { name: 'distances/squared', before: squaredBefore, after: squaredAfter }, + { name: 'similarities/cosine', before: cosineBefore, after: cosineAfter }, + { name: 'distances/clark', before: clarkBefore, after: clarkAfter }, + { name: 'distances/topsoe', before: topsoeBefore, after: topsoeAfter }, + { + name: 'distances/kumarJohnson', + before: kumarJohnsonBefore, + after: kumarJohnsonAfter, + }, + { + name: 'distances/minkowski p=1', + before: (a, b) => minkowskiBeforeP1(a, b, 1), + after: (a, b) => minkowskiAfterP1(a, b, 1), + }, + { + name: 'distances/minkowski p=2', + before: (a, b) => minkowskiBeforeP2(a, b, 2), + after: (a, b) => minkowskiAfterP2(a, b, 2), + }, + { + name: 'similarities/pearson', + before: pearsonSimilarityBefore, + after: pearsonSimilarityAfter, + }, +]; + +function log(message) { + // eslint-disable-next-line no-console -- benchmark output + console.log(message); +} + +function runPair({ name, before, after }, a, b) { + return new Promise((resolve) => { + const beforeValue = before(a, b); + const afterValue = after(a, b); + const identical = Object.is(beforeValue, afterValue); + let beforeHz = 0; + let afterHz = 0; + let beforeRme = 0; + let afterRme = 0; + new Benchmark.Suite(name) + .add('before', () => before(a, b), { maxTime: 2 }) + .add('after', () => after(a, b), { maxTime: 2 }) + .on('cycle', (event) => { + const { name: which, hz, stats } = event.target; + if (which === 'before') { + beforeHz = hz; + beforeRme = stats.rme; + } else { + afterHz = hz; + afterRme = stats.rme; + } + log( + ` ${which.padEnd(7)}${(1e3 / hz).toFixed(4).padStart(10)} ms/op ±${stats.rme.toFixed(2)}% (${stats.sample.length} samples)`, + ); + }) + .on('complete', () => { + const speedup = afterHz / beforeHz; + log( + ` => ${speedup.toFixed(2)}x ${identical ? 'identical result' : `DIFFERS ${beforeValue} -> ${afterValue}`}\n`, + ); + resolve({ speedup, identical, rme: Math.max(beforeRme, afterRme) }); + }) + .run({ async: false }); + }); +} + +const summary = []; +for (const length of LENGTHS) { + for (const kind of KINDS) { + const [a, b] = makeArrays(kind, length); + log(`\n=== ${kind}, ${length} elements ===\n`); + for (const pair of PAIRS) { + log(pair.name); + // eslint-disable-next-line no-await-in-loop -- suites must not overlap + const result = await runPair(pair, a, b); + summary.push({ ...result, name: pair.name, length, kind }); + } + } +} + +log('\n\n=== speedup summary (after / before) ===\n'); +for (const kind of KINDS) { + log(kind); + log( + ` ${'function'.padEnd(24)}${LENGTHS.map((l) => `${l}`.padStart(10)).join('')}`, + ); + for (const pair of PAIRS) { + const cells = LENGTHS.map((length) => { + const found = summary.find( + (entry) => + entry.name === pair.name && + entry.length === length && + entry.kind === kind, + ); + return `${found.speedup.toFixed(2)}x`.padStart(10); + }); + log(` ${pair.name.padEnd(24)}${cells.join('')}`); + } + log(''); +} + +const differing = summary.filter((entry) => !entry.identical); +log( + differing.length === 0 + ? 'all results identical' + : `results differ for: ${[...new Set(differing.map((entry) => entry.name))].join(', ')}`, +); diff --git a/benchmark/elementReads.js b/benchmark/elementReads.js new file mode 100644 index 0000000..f2b5907 --- /dev/null +++ b/benchmark/elementReads.js @@ -0,0 +1,153 @@ +/* + * Why does caching `a[i]` in a local help at all? + * + * Not because of the `NumberArray` TypeScript type: types are erased and V8 + * never sees them. What matters is the element representation the load site + * actually observes at run time. On a site that has only ever seen one + * Float64Array shape, V8 eliminates the repeated loads by itself and the + * caching is worth nothing. On plain arrays, or on a site fed more than one + * array kind, the loads survive and the caching pays. + * + * Run with: node benchmark/elementReads.js (or: bun benchmark/elementReads.js) + */ +import Benchmark from 'benchmark'; + +const LENGTH = 10000; + +function fill(target) { + let state = 42; + for (let i = 0; i < target.length; i++) { + state = (state * 1103515245 + 12345) % 2147483648; + target[i] = 0.1 + (state / 2147483648) * 2; + } + return target; +} + +const typedA = fill(new Float64Array(LENGTH)); +const typedB = fill(new Float64Array(LENGTH)); +const plainA = fill(Array.from({ length: LENGTH }, () => 0)); +const plainB = fill(Array.from({ length: LENGTH }, () => 0)); + +/* + * Six copies of the same two loops. They must not be shared: a single copy + * called with several array kinds would make every measurement polymorphic, + * which is exactly the variable under test. + */ +function typedBefore(a, b) { + let d = 0; + for (let i = 0; i < a.length; i++) { + d += ((a[i] - b[i]) * (a[i] - b[i])) / (a[i] + b[i]); + } + return d; +} +function typedAfter(a, b) { + let d = 0; + for (let i = 0; i < a.length; i++) { + const ai = a[i]; + const bi = b[i]; + const diff = ai - bi; + d += (diff * diff) / (ai + bi); + } + return d; +} +function plainBefore(a, b) { + let d = 0; + for (let i = 0; i < a.length; i++) { + d += ((a[i] - b[i]) * (a[i] - b[i])) / (a[i] + b[i]); + } + return d; +} +function plainAfter(a, b) { + let d = 0; + for (let i = 0; i < a.length; i++) { + const ai = a[i]; + const bi = b[i]; + const diff = ai - bi; + d += (diff * diff) / (ai + bi); + } + return d; +} +function mixedBefore(a, b) { + let d = 0; + for (let i = 0; i < a.length; i++) { + d += ((a[i] - b[i]) * (a[i] - b[i])) / (a[i] + b[i]); + } + return d; +} +function mixedAfter(a, b) { + let d = 0; + for (let i = 0; i < a.length; i++) { + const ai = a[i]; + const bi = b[i]; + const diff = ai - bi; + d += (diff * diff) / (ai + bi); + } + return d; +} + +for (let i = 0; i < 5000; i++) { + typedBefore(typedA, typedB); + typedAfter(typedA, typedB); + plainBefore(plainA, plainB); + plainAfter(plainA, plainB); + mixedBefore(typedA, typedB); + mixedBefore(plainA, plainB); + mixedAfter(typedA, typedB); + mixedAfter(plainA, plainB); +} + +function log(message) { + // eslint-disable-next-line no-console -- benchmark output + console.log(message); +} + +const times = new Map(); +new Benchmark.Suite() + .add('Float64Array only before', () => typedBefore(typedA, typedB), { + minSamples: 30, + }) + .add('Float64Array only after', () => typedAfter(typedA, typedB), { + minSamples: 30, + }) + .add('Array only before', () => plainBefore(plainA, plainB), { + minSamples: 30, + }) + .add('Array only after', () => plainAfter(plainA, plainB), { + minSamples: 30, + }) + .add('both kinds before', () => mixedBefore(typedA, typedB), { + minSamples: 30, + }) + .add('both kinds after', () => mixedAfter(typedA, typedB), { + minSamples: 30, + }) + .on('cycle', (event) => { + const { name, hz, stats } = event.target; + times.set(name, 1e9 / hz / LENGTH); + log( + `${name.padEnd(28)}${(1e9 / hz / LENGTH).toFixed(2).padStart(7)} ns/element ±${stats.rme.toFixed(1)}%`, + ); + }) + .on('complete', () => { + log('\ngain from caching the element reads:'); + for (const [label, before, after] of [ + [ + 'Float64Array only', + 'Float64Array only before', + 'Float64Array only after', + ], + [ + 'Array only ', + 'Array only before', + 'Array only after', + ], + [ + 'both kinds ', + 'both kinds before', + 'both kinds after', + ], + ]) { + log(` ${label} ${(times.get(before) / times.get(after)).toFixed(2)}x`); + } + }) + .run({ async: false }); diff --git a/benchmark/equivalence.js b/benchmark/equivalence.js new file mode 100644 index 0000000..000ec69 --- /dev/null +++ b/benchmark/equivalence.js @@ -0,0 +1,137 @@ +/* + * Differential test: is `x ** 1.5` bit-identical to `x * Math.sqrt(x)`? + * + * Mathematically x^1.5 = x * sqrt(x), but the two expressions round + * differently: the multiplication form rounds twice (sqrt, then *), while `**` + * calls the engine's pow. Neither is guaranteed correctly rounded by the spec + * (ECMA-262 leaves Math.pow implementation-defined), so this measures the + * actual disagreement instead of assuming it away. + * + * Run with: node benchmark/equivalence.js + */ + +function log(message) { + // eslint-disable-next-line no-console -- benchmark output + console.log(message); +} + +const buffer = new DataView(new ArrayBuffer(8)); + +/** + * Distance in representable doubles (ULPs) between two finite values. + * @param x + * @param y + */ +function ulpDistance(x, y) { + if (x === y) return 0n; + if (Number.isNaN(x) || Number.isNaN(y)) return null; + return ordinal(y) - ordinal(x); +} + +/** + * Maps a double onto a monotonic signed integer, so subtraction counts ULPs. + * @param value + */ +function ordinal(value) { + buffer.setFloat64(0, value); + const bits = buffer.getBigUint64(0); + return bits & 0x8000000000000000n + ? -(bits & 0x7fffffffffffffffn) + : BigInt(bits); +} + +log('--- special values ---'); +for (const x of [ + 0, + -0, + 1, + -1, + Infinity, + -Infinity, + Number.NaN, + Number.MIN_VALUE, + Number.MAX_VALUE, + Number.EPSILON, +]) { + const pow = x ** 1.5; + const mul = x * Math.sqrt(x); + const same = Object.is(pow, mul); + log( + ` x = ${String(x).padEnd(24)} x**1.5 = ${String(pow).padEnd(24)} x*sqrt(x) = ${String(mul).padEnd(24)} ${same ? 'same' : 'DIFFERENT'}`, + ); +} + +log('\n--- random sweep across magnitudes ---'); +let state = 88172645463325252n; +const mask = (1n << 64n) - 1n; +function nextBits() { + state ^= (state << 13n) & mask; + state ^= state >> 7n; + state ^= (state << 17n) & mask; + return state; +} + +for (const [label, low, high] of [ + ['[1e-300, 1e-200)', 1e-300, 1e-200], + ['[1e-10, 1e-5)', 1e-10, 1e-5], + ['[0.1, 2)', 0.1, 2], + ['[1, 1000)', 1, 1000], + ['[1e100, 1e200)', 1e100, 1e200], +]) { + const SAMPLES = 2_000_000; + const logLow = Math.log(low); + const logSpan = Math.log(high) - logLow; + let identical = 0; + let maxUlp = 0n; + let maxRelative = 0; + for (let i = 0; i < SAMPLES; i++) { + const unit = Number(nextBits() >> 11n) / 2 ** 53; + const x = Math.exp(logLow + unit * logSpan); + const pow = x ** 1.5; + const mul = x * Math.sqrt(x); + if (pow === mul) { + identical++; + continue; + } + const ulps = ulpDistance(pow, mul); + const absolute = ulps < 0n ? -ulps : ulps; + if (absolute > maxUlp) maxUlp = absolute; + const relative = Math.abs(pow - mul) / Math.abs(pow); + if (relative > maxRelative) maxRelative = relative; + } + const percent = ((identical / SAMPLES) * 100).toFixed(4); + log( + ` ${label.padEnd(18)} identical: ${percent}% max |diff|: ${maxUlp} ulp max relative: ${maxRelative.toExponential(3)}`, + ); +} + +log('\n--- exhaustive mantissa scan in [1, 2) ---'); +{ + // Walk consecutive doubles from 1.0 upward: every value tested is adjacent to + // the previous one, so this is exhaustive over the scanned prefix. + const STEPS = 5_000_000; + let identical = 0; + let maxUlp = 0n; + let x = 1; + for (let i = 0; i < STEPS; i++) { + const pow = x ** 1.5; + const mul = x * Math.sqrt(x); + if (pow === mul) { + identical++; + } else { + const ulps = ulpDistance(pow, mul); + const absolute = ulps < 0n ? -ulps : ulps; + if (absolute > maxUlp) maxUlp = absolute; + } + x = nextUp(x); + } + log( + ` ${STEPS} consecutive doubles from 1.0: identical ${((identical / STEPS) * 100).toFixed(4)}% max |diff|: ${maxUlp} ulp`, + ); +} + +function nextUp(value) { + buffer.setFloat64(0, value); + buffer.setBigUint64(0, buffer.getBigUint64(0) + 1n); + return buffer.getFloat64(0); +} diff --git a/benchmark/exponentiation.js b/benchmark/exponentiation.js new file mode 100644 index 0000000..0f3b542 --- /dev/null +++ b/benchmark/exponentiation.js @@ -0,0 +1,184 @@ +/* + * The exponentiation operator is the reason minkowski and kumarJohnson were + * slow. This isolates which forms of `**` an engine specializes and which fall + * back to a generic pow call. + * + * Run with: node benchmark/exponentiation.js (or: bun benchmark/exponentiation.js) + */ +import Benchmark from 'benchmark'; + +const LENGTH = 10000; +const values = new Float64Array(LENGTH); +let state = 42; +for (let i = 0; i < LENGTH; i++) { + state = (state * 1103515245 + 12345) % 2147483648; + values[i] = 0.1 + (state / 2147483648) * 2; +} + +// A variable exponent, opaque to the compiler. +function exponent(value) { + return values.length > 0 ? value : 0; +} +const variableOne = exponent(1); +const variableTwo = exponent(2); +const variableOneAndAHalf = exponent(1.5); +const variableHalf = exponent(0.5); + +// `power` reaches powerCosine as a parameter with a default, not as a constant. +function poweredSum(input, power) { + let s = 0; + for (let i = 0; i < input.length; i++) s += input[i] ** power; + return s; +} +function sqrtSum(input) { + let s = 0; + for (let i = 0; i < input.length; i++) s += Math.sqrt(input[i]); + return s; +} + +const CASES = [ + [ + 'x ** 2 (literal)', + () => { + let s = 0; + for (let i = 0; i < LENGTH; i++) s += values[i] ** 2; + return s; + }, + ], + [ + 'x * x', + () => { + let s = 0; + for (let i = 0; i < LENGTH; i++) { + const v = values[i]; + s += v * v; + } + return s; + }, + ], + [ + 'x ** p (p = 2)', + () => { + let s = 0; + for (let i = 0; i < LENGTH; i++) s += values[i] ** variableTwo; + return s; + }, + ], + [ + 'x ** 1 (literal)', + () => { + let s = 0; + for (let i = 0; i < LENGTH; i++) s += values[i] ** 1; + return s; + }, + ], + [ + 'x ** p (p = 1)', + () => { + let s = 0; + for (let i = 0; i < LENGTH; i++) s += values[i] ** variableOne; + return s; + }, + ], + [ + 'x (identity)', + () => { + let s = 0; + for (let i = 0; i < LENGTH; i++) s += values[i]; + return s; + }, + ], + [ + 'x ** 0.5 (literal)', + () => { + let s = 0; + for (let i = 0; i < LENGTH; i++) s += values[i] ** 0.5; + return s; + }, + ], + [ + 'Math.sqrt(x)', + () => { + let s = 0; + for (let i = 0; i < LENGTH; i++) s += Math.sqrt(values[i]); + return s; + }, + ], + [ + 'x ** 1.5 (literal)', + () => { + let s = 0; + for (let i = 0; i < LENGTH; i++) s += values[i] ** 1.5; + return s; + }, + ], + [ + 'x ** p (p = 1.5)', + () => { + let s = 0; + for (let i = 0; i < LENGTH; i++) s += values[i] ** variableOneAndAHalf; + return s; + }, + ], + [ + 'x * Math.sqrt(x)', + () => { + let s = 0; + for (let i = 0; i < LENGTH; i++) { + const v = values[i]; + s += v * Math.sqrt(v); + } + return s; + }, + ], + [ + 'x ** p (p = 0.5)', + () => { + let s = 0; + for (let i = 0; i < LENGTH; i++) s += values[i] ** variableHalf; + return s; + }, + ], + ['x ** p (p = 0.5 param)', () => poweredSum(values, 0.5)], + ['Math.sqrt(x) (in fn)', () => sqrtSum(values)], +]; + +function log(message) { + // eslint-disable-next-line no-console -- benchmark output + console.log(message); +} + +const suite = new Benchmark.Suite(); +for (const [name, callback] of CASES) { + suite.add(name, callback, { maxTime: 2 }); +} + +const timings = new Map(); +suite + .on('cycle', (event) => { + const { name, hz, stats } = event.target; + const nanosecondsPerElement = 1e9 / hz / LENGTH; + timings.set(name, nanosecondsPerElement); + log( + `${name.padEnd(22)}${nanosecondsPerElement.toFixed(3).padStart(8)} ns/element ±${stats.rme.toFixed(2)}%`, + ); + }) + .on('complete', () => { + log('\nrelative cost (vs the multiplication-based equivalent):'); + for (const [slow, fast] of [ + ['x ** 2 (literal)', 'x * x'], + ['x ** p (p = 2)', 'x * x'], + ['x ** 1 (literal)', 'x (identity)'], + ['x ** p (p = 1)', 'x (identity)'], + ['x ** 0.5 (literal)', 'Math.sqrt(x)'], + ['x ** p (p = 0.5)', 'Math.sqrt(x)'], + ['x ** p (p = 0.5 param)', 'Math.sqrt(x) (in fn)'], + ['x ** 1.5 (literal)', 'x * Math.sqrt(x)'], + ['x ** p (p = 1.5)', 'x * Math.sqrt(x)'], + ]) { + log( + ` ${slow.padEnd(22)}${(timings.get(slow) / timings.get(fast)).toFixed(2).padStart(7)}x the cost of ${fast}`, + ); + } + }) + .run({ async: false }); diff --git a/package.json b/package.json index b458f0c..72a88bf 100644 --- a/package.json +++ b/package.json @@ -57,6 +57,7 @@ "@types/node": "^26.1.1", "@vitest/coverage-v8": "^4.1.10", "@zakodium/tsconfig": "^1.0.5", + "benchmark": "^2.1.4", "eslint": "^9.39.5", "eslint-config-cheminfo-typescript": "^22.1.0", "prettier": "^3.9.6", From 6f395228e545ce0cce9863229c14f8c5d8958c87 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Sun, 26 Jul 2026 21:43:05 +0200 Subject: [PATCH 2/4] perf: avoid slow exponentiation and repeated element reads Three algorithmic changes, each large and stable across engines: - kumarJohnson: `prod ** 1.5` -> `prod * Math.sqrt(prod)`. No engine specializes a 1.5 exponent, so this is the single biggest win. - minkowski: special-case p = 1 and p = 2. V8 runs `x ** 1` through generic pow at ~9x the cost of a plain read; p = 2 is worth ~1.5x. - similarity.pearson: drop the two arrays it allocated only to mean-centre before calling cosine, and fold the work into two passes. Measured with beforeAfter.js, before -> after (n = 10000, Float64Array, node 24), against the baseline recorded in the previous commit: kumarJohnson 15.50 -> 2.41 ns/element 6.4x faster similarity.pearson 19.06 -> 3.23 ns/element 5.9x faster minkowski(p=1) 11.51 -> 2.02 ns/element 5.7x faster squared 4.51 -> 2.22 ns/element 2.0x faster cosine 4.38 -> 2.19 ns/element 2.0x faster pearson 3.74 -> 2.14 ns/element 1.7x faster clark 3.22 -> 2.17 ns/element 1.5x faster topsoe 21.99 -> 17.21 ns/element 1.3x faster The remaining change caches repeated `a[i]` reads in locals. Its value is much smaller and depends on what the load site actually observes at run time, not on any TypeScript type: V8 often eliminates the repeated loads by itself. elementReads.js measures anywhere from 1.00x when a site only ever sees one Float64Array shape up to ~1.5x once it has seen several array kinds, and the exact figure moves with how the harness is built. The rows above are the higher end of that range because beforeAfter.js exercises every function with both array kinds in one process. It is never a regression, so it stays, but it is not a reliable 2x. JavaScriptCore eliminates those loads outright: on bun the caching is 1.00x everywhere and only kumarJohnson (4x) and similarity.pearson (6x) improve. kumarJohnson moves by at most 1 ulp, and `(-Infinity) ** 1.5` is Infinity where `-Infinity * Math.sqrt(-Infinity)` is NaN; both are outside the non-negative domain the metric is defined on. Its test now uses toBeCloseTo. ml-array-mean is no longer used. Assisted-By: Claude Opus 5 (1M context) --- package.json | 1 - src/distances/__tests__/kumarJohnson.test.ts | 2 +- src/distances/additiveSymmetric.ts | 5 +++- src/distances/canberra.ts | 4 ++- src/distances/clark.ts | 5 +++- src/distances/dice.ts | 9 ++++-- src/distances/divergence.ts | 6 +++- src/distances/harmonicMean.ts | 4 ++- src/distances/jeffreys.ts | 4 ++- src/distances/jensenDifference.ts | 7 +++-- src/distances/jensenShannon.ts | 7 +++-- src/distances/kdivergence.ts | 3 +- src/distances/kulczynski.ts | 6 ++-- src/distances/kullbackLeibler.ts | 3 +- src/distances/kumarJohnson.ts | 8 ++++- src/distances/minkowski.ts | 15 ++++++++++ src/distances/motyka.ts | 6 ++-- src/distances/neyman.ts | 4 ++- src/distances/pearson.ts | 4 ++- src/distances/probabilisticSymmetric.ts | 5 +++- src/distances/ruzicka.ts | 6 ++-- src/distances/soergel.ts | 6 ++-- src/distances/sorensen.ts | 6 ++-- src/distances/squared.ts | 5 +++- src/distances/squaredChord.ts | 3 +- src/distances/taneja.ts | 7 +++-- src/distances/tanimoto.ts | 8 +++-- src/distances/topsoe.ts | 7 +++-- src/distances/waveHedges.ts | 4 ++- src/similarities/cosine.ts | 8 +++-- src/similarities/czekanowski.ts | 6 ++-- src/similarities/kumarHassebrook.ts | 8 +++-- src/similarities/pearson.ts | 31 ++++++++++++-------- src/similarities/tanimoto.ts | 8 +++-- 34 files changed, 154 insertions(+), 67 deletions(-) diff --git a/package.json b/package.json index 72a88bf..4752b5d 100644 --- a/package.json +++ b/package.json @@ -49,7 +49,6 @@ "homepage": "https://github.com/mljs/distance", "dependencies": { "cheminfo-types": "^1.15.0", - "ml-array-mean": "^2.0.0", "ml-distance-euclidean": "^3.0.1", "ml-tree-similarity": "^1.0.0" }, diff --git a/src/distances/__tests__/kumarJohnson.test.ts b/src/distances/__tests__/kumarJohnson.test.ts index c9bb5de..9dd489e 100644 --- a/src/distances/__tests__/kumarJohnson.test.ts +++ b/src/distances/__tests__/kumarJohnson.test.ts @@ -6,5 +6,5 @@ const v1 = [0.2, 0.4, 0.3, 0.1]; const v2 = [0.3, 0.2, 0.3, 0.2]; test('should be correct', () => { - expect(distance.kumarJohnson(v1, v2)).toBe(0.5623488044808911); + expect(distance.kumarJohnson(v1, v2)).toBeCloseTo(0.5623488044808911, 15); }); diff --git a/src/distances/additiveSymmetric.ts b/src/distances/additiveSymmetric.ts index 45d4e30..8b84ac3 100644 --- a/src/distances/additiveSymmetric.ts +++ b/src/distances/additiveSymmetric.ts @@ -8,7 +8,10 @@ import type { NumberArray } from 'cheminfo-types'; export function additiveSymmetric(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - d += ((a[i] - b[i]) * (a[i] - b[i]) * (a[i] + b[i])) / (a[i] * b[i]); + const ai = a[i]; + const bi = b[i]; + const diff = ai - bi; + d += (diff * diff * (ai + bi)) / (ai * bi); } return d; } diff --git a/src/distances/canberra.ts b/src/distances/canberra.ts index aa1ab3e..15237d6 100644 --- a/src/distances/canberra.ts +++ b/src/distances/canberra.ts @@ -8,7 +8,9 @@ import type { NumberArray } from 'cheminfo-types'; export function canberra(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - ans += Math.abs(a[i] - b[i]) / (a[i] + b[i]); + const ai = a[i]; + const bi = b[i]; + ans += Math.abs(ai - bi) / (ai + bi); } return ans; } diff --git a/src/distances/clark.ts b/src/distances/clark.ts index faea766..fe5a466 100644 --- a/src/distances/clark.ts +++ b/src/distances/clark.ts @@ -8,7 +8,10 @@ import type { NumberArray } from 'cheminfo-types'; export function clark(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - d += (Math.abs(a[i] - b[i]) / (a[i] + b[i])) ** 2; + const ai = a[i]; + const bi = b[i]; + const ratio = (ai - bi) / (ai + bi); + d += ratio * ratio; } return Math.sqrt(d); } diff --git a/src/distances/dice.ts b/src/distances/dice.ts index 7696bc1..674dd86 100644 --- a/src/distances/dice.ts +++ b/src/distances/dice.ts @@ -10,9 +10,12 @@ export function dice(a: NumberArray, b: NumberArray): number { let b2 = 0; let prod2 = 0; for (let i = 0; i < a.length; i++) { - a2 += a[i] * a[i]; - b2 += b[i] * b[i]; - prod2 += (a[i] - b[i]) * (a[i] - b[i]); + const ai = a[i]; + const bi = b[i]; + const diff = ai - bi; + a2 += ai * ai; + b2 += bi * bi; + prod2 += diff * diff; } return prod2 / (a2 + b2); } diff --git a/src/distances/divergence.ts b/src/distances/divergence.ts index 1b62dd4..2039014 100644 --- a/src/distances/divergence.ts +++ b/src/distances/divergence.ts @@ -8,7 +8,11 @@ import type { NumberArray } from 'cheminfo-types'; export function divergence(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - d += ((a[i] - b[i]) * (a[i] - b[i])) / ((a[i] + b[i]) * (a[i] + b[i])); + const ai = a[i]; + const bi = b[i]; + const diff = ai - bi; + const sum = ai + bi; + d += (diff * diff) / (sum * sum); } return 2 * d; } diff --git a/src/distances/harmonicMean.ts b/src/distances/harmonicMean.ts index dde08dd..eff5777 100644 --- a/src/distances/harmonicMean.ts +++ b/src/distances/harmonicMean.ts @@ -8,7 +8,9 @@ import type { NumberArray } from 'cheminfo-types'; export function harmonicMean(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - ans += (a[i] * b[i]) / (a[i] + b[i]); + const ai = a[i]; + const bi = b[i]; + ans += (ai * bi) / (ai + bi); } return 2 * ans; } diff --git a/src/distances/jeffreys.ts b/src/distances/jeffreys.ts index cdaec4e..a5e0220 100644 --- a/src/distances/jeffreys.ts +++ b/src/distances/jeffreys.ts @@ -8,7 +8,9 @@ import type { NumberArray } from 'cheminfo-types'; export function jeffreys(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - ans += (a[i] - b[i]) * Math.log(a[i] / b[i]); + const ai = a[i]; + const bi = b[i]; + ans += (ai - bi) * Math.log(ai / bi); } return ans; } diff --git a/src/distances/jensenDifference.ts b/src/distances/jensenDifference.ts index 6beaf95..582ffe0 100644 --- a/src/distances/jensenDifference.ts +++ b/src/distances/jensenDifference.ts @@ -8,9 +8,10 @@ import type { NumberArray } from 'cheminfo-types'; export function jensenDifference(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - ans += - (a[i] * Math.log(a[i]) + b[i] * Math.log(b[i])) / 2 - - ((a[i] + b[i]) / 2) * Math.log((a[i] + b[i]) / 2); + const ai = a[i]; + const bi = b[i]; + const half = (ai + bi) / 2; + ans += (ai * Math.log(ai) + bi * Math.log(bi)) / 2 - half * Math.log(half); } return ans; } diff --git a/src/distances/jensenShannon.ts b/src/distances/jensenShannon.ts index ededd6e..e625ae7 100644 --- a/src/distances/jensenShannon.ts +++ b/src/distances/jensenShannon.ts @@ -9,8 +9,11 @@ export function jensenShannon(a: NumberArray, b: NumberArray): number { let p = 0; let q = 0; for (let i = 0; i < a.length; i++) { - p += a[i] * Math.log((2 * a[i]) / (a[i] + b[i])); - q += b[i] * Math.log((2 * b[i]) / (a[i] + b[i])); + const ai = a[i]; + const bi = b[i]; + const sum = ai + bi; + p += ai * Math.log((2 * ai) / sum); + q += bi * Math.log((2 * bi) / sum); } return (p + q) / 2; } diff --git a/src/distances/kdivergence.ts b/src/distances/kdivergence.ts index 915dcb9..6dbf995 100644 --- a/src/distances/kdivergence.ts +++ b/src/distances/kdivergence.ts @@ -8,7 +8,8 @@ import type { NumberArray } from 'cheminfo-types'; export function kdivergence(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - ans += a[i] * Math.log((2 * a[i]) / (a[i] + b[i])); + const ai = a[i]; + ans += ai * Math.log((2 * ai) / (ai + b[i])); } return ans; } diff --git a/src/distances/kulczynski.ts b/src/distances/kulczynski.ts index 7e2fe67..f8c2d5f 100644 --- a/src/distances/kulczynski.ts +++ b/src/distances/kulczynski.ts @@ -9,8 +9,10 @@ export function kulczynski(a: NumberArray, b: NumberArray): number { let up = 0; let down = 0; for (let i = 0; i < a.length; i++) { - up += Math.abs(a[i] - b[i]); - down += Math.min(a[i], b[i]); + const ai = a[i]; + const bi = b[i]; + up += Math.abs(ai - bi); + down += Math.min(ai, bi); } return up / down; } diff --git a/src/distances/kullbackLeibler.ts b/src/distances/kullbackLeibler.ts index 0df38bb..e5dac0d 100644 --- a/src/distances/kullbackLeibler.ts +++ b/src/distances/kullbackLeibler.ts @@ -8,7 +8,8 @@ import type { NumberArray } from 'cheminfo-types'; export function kullbackLeibler(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - ans += a[i] * Math.log(a[i] / b[i]); + const ai = a[i]; + ans += ai * Math.log(ai / b[i]); } return ans; } diff --git a/src/distances/kumarJohnson.ts b/src/distances/kumarJohnson.ts index 208309f..a61b103 100644 --- a/src/distances/kumarJohnson.ts +++ b/src/distances/kumarJohnson.ts @@ -8,7 +8,13 @@ import type { NumberArray } from 'cheminfo-types'; export function kumarJohnson(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - ans += (a[i] * a[i] - b[i] * b[i]) ** 2 / (2 * (a[i] * b[i]) ** 1.5); + const ai = a[i]; + const bi = b[i]; + const numerator = ai * ai - bi * bi; + // `prod * Math.sqrt(prod)` is ~6x faster than `prod ** 1.5`, which no + // engine specializes; it costs at most 1 ulp of accuracy + const prod = ai * bi; + ans += (numerator * numerator) / (2 * prod * Math.sqrt(prod)); } return ans; } diff --git a/src/distances/minkowski.ts b/src/distances/minkowski.ts index e0639c4..e241504 100644 --- a/src/distances/minkowski.ts +++ b/src/distances/minkowski.ts @@ -8,6 +8,21 @@ import type { NumberArray } from 'cheminfo-types'; */ export function minkowski(a: NumberArray, b: NumberArray, p: number) { let d = 0; + // `x ** p` is far slower than the equivalent multiplication: ~9x for p = 1 + // and ~1.5x for p = 2, the two orders that are used in practice. + if (p === 1) { + for (let i = 0; i < a.length; i++) { + d += Math.abs(a[i] - b[i]); + } + return d; + } + if (p === 2) { + for (let i = 0; i < a.length; i++) { + const diff = a[i] - b[i]; + d += diff * diff; + } + return Math.sqrt(d); + } for (let i = 0; i < a.length; i++) { d += Math.abs(a[i] - b[i]) ** p; } diff --git a/src/distances/motyka.ts b/src/distances/motyka.ts index b7f1b40..1db4325 100644 --- a/src/distances/motyka.ts +++ b/src/distances/motyka.ts @@ -9,8 +9,10 @@ export function motyka(a: NumberArray, b: NumberArray): number { let up = 0; let down = 0; for (let i = 0; i < a.length; i++) { - up += Math.min(a[i], b[i]); - down += a[i] + b[i]; + const ai = a[i]; + const bi = b[i]; + up += Math.min(ai, bi); + down += ai + bi; } return 1 - up / down; } diff --git a/src/distances/neyman.ts b/src/distances/neyman.ts index e406bd1..3e567a0 100644 --- a/src/distances/neyman.ts +++ b/src/distances/neyman.ts @@ -8,7 +8,9 @@ import type { NumberArray } from 'cheminfo-types'; export function neyman(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - d += ((a[i] - b[i]) * (a[i] - b[i])) / a[i]; + const ai = a[i]; + const diff = ai - b[i]; + d += (diff * diff) / ai; } return d; } diff --git a/src/distances/pearson.ts b/src/distances/pearson.ts index 60ce11e..da4dbe4 100644 --- a/src/distances/pearson.ts +++ b/src/distances/pearson.ts @@ -8,7 +8,9 @@ import type { NumberArray } from 'cheminfo-types'; export function pearson(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - d += ((a[i] - b[i]) * (a[i] - b[i])) / b[i]; + const bi = b[i]; + const diff = a[i] - bi; + d += (diff * diff) / bi; } return d; } diff --git a/src/distances/probabilisticSymmetric.ts b/src/distances/probabilisticSymmetric.ts index 5d310f0..e187255 100644 --- a/src/distances/probabilisticSymmetric.ts +++ b/src/distances/probabilisticSymmetric.ts @@ -8,7 +8,10 @@ import type { NumberArray } from 'cheminfo-types'; export function probabilisticSymmetric(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - d += ((a[i] - b[i]) * (a[i] - b[i])) / (a[i] + b[i]); + const ai = a[i]; + const bi = b[i]; + const diff = ai - bi; + d += (diff * diff) / (ai + bi); } return 2 * d; } diff --git a/src/distances/ruzicka.ts b/src/distances/ruzicka.ts index 3554691..e24d5b5 100644 --- a/src/distances/ruzicka.ts +++ b/src/distances/ruzicka.ts @@ -9,8 +9,10 @@ export function ruzicka(a: NumberArray, b: NumberArray): number { let up = 0; let down = 0; for (let i = 0; i < a.length; i++) { - up += Math.min(a[i], b[i]); - down += Math.max(a[i], b[i]); + const ai = a[i]; + const bi = b[i]; + up += Math.min(ai, bi); + down += Math.max(ai, bi); } return up / down; } diff --git a/src/distances/soergel.ts b/src/distances/soergel.ts index ef00338..2995d8e 100644 --- a/src/distances/soergel.ts +++ b/src/distances/soergel.ts @@ -10,8 +10,10 @@ export function soergel(a: NumberArray, b: NumberArray): number { let up = 0; let down = 0; for (let i = 0; i < a.length; i++) { - up += Math.abs(a[i] - b[i]); - down += Math.max(a[i], b[i]); + const ai = a[i]; + const bi = b[i]; + up += Math.abs(ai - bi); + down += Math.max(ai, bi); } return up / down; } diff --git a/src/distances/sorensen.ts b/src/distances/sorensen.ts index c13cb3c..771f9da 100644 --- a/src/distances/sorensen.ts +++ b/src/distances/sorensen.ts @@ -10,8 +10,10 @@ export function sorensen(a: NumberArray, b: NumberArray): number { let up = 0; let down = 0; for (let i = 0; i < a.length; i++) { - up += Math.abs(a[i] - b[i]); - down += a[i] + b[i]; + const ai = a[i]; + const bi = b[i]; + up += Math.abs(ai - bi); + down += ai + bi; } return up / down; } diff --git a/src/distances/squared.ts b/src/distances/squared.ts index 3565dae..c9cba2a 100644 --- a/src/distances/squared.ts +++ b/src/distances/squared.ts @@ -8,7 +8,10 @@ import type { NumberArray } from 'cheminfo-types'; export function squared(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - d += ((a[i] - b[i]) * (a[i] - b[i])) / (a[i] + b[i]); + const ai = a[i]; + const bi = b[i]; + const diff = ai - bi; + d += (diff * diff) / (ai + bi); } return d; } diff --git a/src/distances/squaredChord.ts b/src/distances/squaredChord.ts index b62e45d..d613f4c 100644 --- a/src/distances/squaredChord.ts +++ b/src/distances/squaredChord.ts @@ -8,7 +8,8 @@ import type { NumberArray } from 'cheminfo-types'; export function squaredChord(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - ans += (Math.sqrt(a[i]) - Math.sqrt(b[i])) ** 2; + const diff = Math.sqrt(a[i]) - Math.sqrt(b[i]); + ans += diff * diff; } return ans; } diff --git a/src/distances/taneja.ts b/src/distances/taneja.ts index 12a0df3..8c9dacf 100644 --- a/src/distances/taneja.ts +++ b/src/distances/taneja.ts @@ -8,9 +8,10 @@ import type { NumberArray } from 'cheminfo-types'; export function taneja(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - ans += - ((a[i] + b[i]) / 2) * - Math.log((a[i] + b[i]) / (2 * Math.sqrt(a[i] * b[i]))); + const ai = a[i]; + const bi = b[i]; + const sum = ai + bi; + ans += (sum / 2) * Math.log(sum / (2 * Math.sqrt(ai * bi))); } return ans; } diff --git a/src/distances/tanimoto.ts b/src/distances/tanimoto.ts index a0b1260..d48336d 100644 --- a/src/distances/tanimoto.ts +++ b/src/distances/tanimoto.ts @@ -20,9 +20,11 @@ export function tanimoto( let q = 0; let m = 0; for (let i = 0; i < a.length; i++) { - p += a[i]; - q += b[i]; - m += Math.min(a[i], b[i]); + const ai = a[i]; + const bi = b[i]; + p += ai; + q += bi; + m += Math.min(ai, bi); } return (p + q - 2 * m) / (p + q - m); } diff --git a/src/distances/topsoe.ts b/src/distances/topsoe.ts index 3df7cd4..126e675 100644 --- a/src/distances/topsoe.ts +++ b/src/distances/topsoe.ts @@ -8,9 +8,10 @@ import type { NumberArray } from 'cheminfo-types'; export function topsoe(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - ans += - a[i] * Math.log((2 * a[i]) / (a[i] + b[i])) + - b[i] * Math.log((2 * b[i]) / (a[i] + b[i])); + const ai = a[i]; + const bi = b[i]; + const sum = ai + bi; + ans += ai * Math.log((2 * ai) / sum) + bi * Math.log((2 * bi) / sum); } return ans; } diff --git a/src/distances/waveHedges.ts b/src/distances/waveHedges.ts index 86ceae1..108c0be 100644 --- a/src/distances/waveHedges.ts +++ b/src/distances/waveHedges.ts @@ -8,7 +8,9 @@ import type { NumberArray } from 'cheminfo-types'; export function waveHedges(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - ans += 1 - Math.min(a[i], b[i]) / Math.max(a[i], b[i]); + const ai = a[i]; + const bi = b[i]; + ans += 1 - Math.min(ai, bi) / Math.max(ai, bi); } return ans; } diff --git a/src/similarities/cosine.ts b/src/similarities/cosine.ts index 3857896..0bedbac 100644 --- a/src/similarities/cosine.ts +++ b/src/similarities/cosine.ts @@ -9,9 +9,11 @@ export function cosine(a: NumberArray, b: NumberArray): number { let p2 = 0; let q2 = 0; for (let i = 0; i < a.length; i++) { - p += a[i] * b[i]; - p2 += a[i] * a[i]; - q2 += b[i] * b[i]; + const ai = a[i]; + const bi = b[i]; + p += ai * bi; + p2 += ai * ai; + q2 += bi * bi; } return p / (Math.sqrt(p2) * Math.sqrt(q2)); } diff --git a/src/similarities/czekanowski.ts b/src/similarities/czekanowski.ts index 7ca8039..3dc6fee 100644 --- a/src/similarities/czekanowski.ts +++ b/src/similarities/czekanowski.ts @@ -9,8 +9,10 @@ export function czekanowski(a: NumberArray, b: NumberArray): number { let up = 0; let down = 0; for (let i = 0; i < a.length; i++) { - up += Math.min(a[i], b[i]); - down += a[i] + b[i]; + const ai = a[i]; + const bi = b[i]; + up += Math.min(ai, bi); + down += ai + bi; } return (2 * up) / down; } diff --git a/src/similarities/kumarHassebrook.ts b/src/similarities/kumarHassebrook.ts index c2bba35..fd0a8e6 100644 --- a/src/similarities/kumarHassebrook.ts +++ b/src/similarities/kumarHassebrook.ts @@ -10,9 +10,11 @@ export function kumarHassebrook(a: NumberArray, b: NumberArray): number { let p2 = 0; let q2 = 0; for (let i = 0; i < a.length; i++) { - p += a[i] * b[i]; - p2 += a[i] * a[i]; - q2 += b[i] * b[i]; + const ai = a[i]; + const bi = b[i]; + p += ai * bi; + p2 += ai * ai; + q2 += bi * bi; } return p / (p2 + q2 - p); } diff --git a/src/similarities/pearson.ts b/src/similarities/pearson.ts index 59d7e4c..442d3e0 100644 --- a/src/similarities/pearson.ts +++ b/src/similarities/pearson.ts @@ -1,7 +1,4 @@ import type { NumberArray } from 'cheminfo-types'; -import mean from 'ml-array-mean'; - -import { cosine } from './cosine.ts'; /** * Returns the Pearson correlation between vectors a and b, i.e. the cosine @@ -10,15 +7,25 @@ import { cosine } from './cosine.ts'; * @param b - second vector */ export function pearson(a: NumberArray, b: NumberArray): number { - const avgA = mean(a); - const avgB = mean(b); - - const newA = new Array(a.length); - const newB = new Array(b.length); - for (let i = 0; i < newA.length; i++) { - newA[i] = a[i] - avgA; - newB[i] = b[i] - avgB; + const length = a.length; + let sumA = 0; + let sumB = 0; + for (let i = 0; i < length; i++) { + sumA += a[i]; + sumB += b[i]; } + const avgA = sumA / length; + const avgB = sumB / length; - return cosine(newA, newB); + let p = 0; + let p2 = 0; + let q2 = 0; + for (let i = 0; i < length; i++) { + const centredA = a[i] - avgA; + const centredB = b[i] - avgB; + p += centredA * centredB; + p2 += centredA * centredA; + q2 += centredB * centredB; + } + return p / (Math.sqrt(p2) * Math.sqrt(q2)); } diff --git a/src/similarities/tanimoto.ts b/src/similarities/tanimoto.ts index 371d1c8..11439ff 100644 --- a/src/similarities/tanimoto.ts +++ b/src/similarities/tanimoto.ts @@ -27,9 +27,11 @@ export function tanimoto( let q = 0; let m = 0; for (let i = 0; i < a.length; i++) { - p += a[i]; - q += b[i]; - m += Math.min(a[i], b[i]); + const ai = a[i]; + const bi = b[i]; + p += ai; + q += bi; + m += Math.min(ai, bi); } return 1 - (p + q - 2 * m) / (p + q - m); } From 8401b5f749cea8d25af86c2b4dfff95b95578ba3 Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Wed, 29 Jul 2026 16:04:46 +0200 Subject: [PATCH 3/4] benchmark: simple benchmark without library --- benchmark/arrayKinds.js | 105 ++++++++++++++++++++++++++++++++++++++++ benchmark/arrayKinds.sh | 41 ++++++++++++++++ 2 files changed, 146 insertions(+) create mode 100644 benchmark/arrayKinds.js create mode 100755 benchmark/arrayKinds.sh diff --git a/benchmark/arrayKinds.js b/benchmark/arrayKinds.js new file mode 100644 index 0000000..5936682 --- /dev/null +++ b/benchmark/arrayKinds.js @@ -0,0 +1,105 @@ +/* + * Cost of `(ai * bi) / (ai + bi)` depending on the element representation. + * One kind per process, so the loads stay monomorphic. + * + * `cached` copies `x[i]` into a local, `direct` reads it on every use. + * + * Run with: node benchmark/arrayKinds.js + */ +import { argv } from 'node:process'; + +const LENGTH = 10000; +const WARMUP_MS = 1000; +const TARGET_MS = 2000; // keep at 1000 or more, shorter runs are too noisy + +const kind = argv[2] ?? 'typed'; +const reads = argv[3] ?? 'cached'; + +function fill(target) { + let state = 42; + for (let i = 0; i < target.length; i++) { + state = (state * 1103515245 + 12345) % 2147483648; + target[i] = 0.1 + (state / 2147483648) * 2; + } + return target; +} + +function plainVector() { + return fill(Array.from({ length: LENGTH }, () => 0)); +} + +function typedVector() { + return fill(new Float64Array(LENGTH)); +} + +let a; +let b; +if (kind === 'array') { + a = plainVector(); + b = plainVector(); +} else if (kind === 'typed') { + a = typedVector(); + b = typedVector(); +} else if (kind === 'mixed') { + a = typedVector(); + b = plainVector(); +} else { + throw new Error(`unknown kind: ${kind}`); +} + +function cachedKernel(x, y) { + let sum = 0; + for (let i = 0; i < x.length; i++) { + const xi = x[i]; + const yi = y[i]; + sum += (xi * yi) / (xi + yi); + } + return sum; +} + +function directKernel(x, y) { + let sum = 0; + for (let i = 0; i < x.length; i++) { + sum += (x[i] * y[i]) / (x[i] + y[i]); + } + return sum; +} + +const kernel = reads === 'cached' ? cachedKernel : directKernel; + +let sink = 0; +// the arguments are swapped after each call, so `mixed` really sees both orders +let x = a; +let y = b; + +let start = performance.now(); +while (performance.now() - start < WARMUP_MS) { + sink += kernel(x, y); + [x, y] = [y, x]; +} + +let rounds = 0; +let elapsed = 0; +start = performance.now(); +while (elapsed < TARGET_MS) { + sink += kernel(x, y); + [x, y] = [y, x]; + rounds++; + elapsed = performance.now() - start; +} + +// keeps the loop from being dropped as dead code +if (!Number.isFinite(sink)) throw new Error('kernel diverged'); + +const operationsPerSecond = (rounds * LENGTH * 1000) / elapsed; +// eslint-disable-next-line no-console -- benchmark output +console.log( + [ + kind, + reads, + operationsPerSecond.toFixed(0), + LENGTH, + kernel(a, b).toFixed(10), + kernel(b, a).toFixed(10), + ].join('\t'), +); diff --git a/benchmark/arrayKinds.sh b/benchmark/arrayKinds.sh new file mode 100755 index 0000000..7f15c1f --- /dev/null +++ b/benchmark/arrayKinds.sh @@ -0,0 +1,41 @@ +#!/bin/bash +# Runs every element representation against both read styles, one process each, +# then prints a summary table. +set -euo pipefail + +directory="$(dirname "$0")" +measurements="" + +for kind in array typed mixed; do + for reads in cached direct; do + measurements="${measurements}$(node "$directory/arrayKinds.js" "$kind" "$reads") +" + done +done + +printf '%s' "$measurements" | awk -F'\t' ' + { + speed[$1 "/" $2] = $3 / 1e6 + length_ = $4 + results[$5] = 1 + results[$6] = 1 + } + END { + printf "\n (ai * bi) / (ai + bi) on %d elements, arguments swapped at each call\n\n", length_ + printf " %-8s %14s %14s %10s\n", "kind", "cached", "direct", "ratio" + printf " %-8s %14s %14s %10s\n", "--------", "--------------", "--------------", "----------" + split("array typed mixed", kinds, " ") + for (i = 1; i <= 3; i++) { + cached = speed[kinds[i] "/cached"] + direct = speed[kinds[i] "/direct"] + printf " %-8s %11.0f M/s %11.0f M/s %9.2fx\n", kinds[i], cached, direct, cached / direct + } + distinct = 0 + for (result in results) { distinct++; sample = result } + if (distinct == 1) { + printf "\n all results identical: %s\n\n", sample + } else { + printf "\n WARNING: %d different results\n\n", distinct + } + } +' From c9ef019f5c78743316d44e7da87a46c2c56f813b Mon Sep 17 00:00:00 2001 From: Luc Patiny Date: Thu, 30 Jul 2026 14:14:59 +0200 Subject: [PATCH 4/4] chore: keep only relevant changes --- src/distances/additiveSymmetric.ts | 5 +---- src/distances/canberra.ts | 4 +--- src/distances/clark.ts | 5 +---- src/distances/dice.ts | 9 +++------ src/distances/divergence.ts | 6 +----- src/distances/harmonicMean.ts | 4 +--- src/distances/jeffreys.ts | 4 +--- src/distances/jensenDifference.ts | 7 +++---- src/distances/jensenShannon.ts | 7 ++----- src/distances/kdivergence.ts | 3 +-- src/distances/kulczynski.ts | 6 ++---- src/distances/kullbackLeibler.ts | 3 +-- src/distances/kumarJohnson.ts | 6 ++---- src/distances/motyka.ts | 6 ++---- src/distances/neyman.ts | 4 +--- src/distances/pearson.ts | 4 +--- src/distances/probabilisticSymmetric.ts | 5 +---- src/distances/ruzicka.ts | 6 ++---- src/distances/soergel.ts | 6 ++---- src/distances/sorensen.ts | 6 ++---- src/distances/squared.ts | 5 +---- src/distances/squaredChord.ts | 3 +-- src/distances/taneja.ts | 7 +++---- src/distances/tanimoto.ts | 8 +++----- src/distances/topsoe.ts | 7 +++---- src/distances/waveHedges.ts | 4 +--- src/similarities/cosine.ts | 8 +++----- src/similarities/czekanowski.ts | 6 ++---- src/similarities/kumarHassebrook.ts | 8 +++----- src/similarities/tanimoto.ts | 8 +++----- 30 files changed, 54 insertions(+), 116 deletions(-) diff --git a/src/distances/additiveSymmetric.ts b/src/distances/additiveSymmetric.ts index 8b84ac3..45d4e30 100644 --- a/src/distances/additiveSymmetric.ts +++ b/src/distances/additiveSymmetric.ts @@ -8,10 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function additiveSymmetric(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - const diff = ai - bi; - d += (diff * diff * (ai + bi)) / (ai * bi); + d += ((a[i] - b[i]) * (a[i] - b[i]) * (a[i] + b[i])) / (a[i] * b[i]); } return d; } diff --git a/src/distances/canberra.ts b/src/distances/canberra.ts index 15237d6..aa1ab3e 100644 --- a/src/distances/canberra.ts +++ b/src/distances/canberra.ts @@ -8,9 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function canberra(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - ans += Math.abs(ai - bi) / (ai + bi); + ans += Math.abs(a[i] - b[i]) / (a[i] + b[i]); } return ans; } diff --git a/src/distances/clark.ts b/src/distances/clark.ts index fe5a466..faea766 100644 --- a/src/distances/clark.ts +++ b/src/distances/clark.ts @@ -8,10 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function clark(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - const ratio = (ai - bi) / (ai + bi); - d += ratio * ratio; + d += (Math.abs(a[i] - b[i]) / (a[i] + b[i])) ** 2; } return Math.sqrt(d); } diff --git a/src/distances/dice.ts b/src/distances/dice.ts index 674dd86..7696bc1 100644 --- a/src/distances/dice.ts +++ b/src/distances/dice.ts @@ -10,12 +10,9 @@ export function dice(a: NumberArray, b: NumberArray): number { let b2 = 0; let prod2 = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - const diff = ai - bi; - a2 += ai * ai; - b2 += bi * bi; - prod2 += diff * diff; + a2 += a[i] * a[i]; + b2 += b[i] * b[i]; + prod2 += (a[i] - b[i]) * (a[i] - b[i]); } return prod2 / (a2 + b2); } diff --git a/src/distances/divergence.ts b/src/distances/divergence.ts index 2039014..1b62dd4 100644 --- a/src/distances/divergence.ts +++ b/src/distances/divergence.ts @@ -8,11 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function divergence(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - const diff = ai - bi; - const sum = ai + bi; - d += (diff * diff) / (sum * sum); + d += ((a[i] - b[i]) * (a[i] - b[i])) / ((a[i] + b[i]) * (a[i] + b[i])); } return 2 * d; } diff --git a/src/distances/harmonicMean.ts b/src/distances/harmonicMean.ts index eff5777..dde08dd 100644 --- a/src/distances/harmonicMean.ts +++ b/src/distances/harmonicMean.ts @@ -8,9 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function harmonicMean(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - ans += (ai * bi) / (ai + bi); + ans += (a[i] * b[i]) / (a[i] + b[i]); } return 2 * ans; } diff --git a/src/distances/jeffreys.ts b/src/distances/jeffreys.ts index a5e0220..cdaec4e 100644 --- a/src/distances/jeffreys.ts +++ b/src/distances/jeffreys.ts @@ -8,9 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function jeffreys(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - ans += (ai - bi) * Math.log(ai / bi); + ans += (a[i] - b[i]) * Math.log(a[i] / b[i]); } return ans; } diff --git a/src/distances/jensenDifference.ts b/src/distances/jensenDifference.ts index 582ffe0..6beaf95 100644 --- a/src/distances/jensenDifference.ts +++ b/src/distances/jensenDifference.ts @@ -8,10 +8,9 @@ import type { NumberArray } from 'cheminfo-types'; export function jensenDifference(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - const half = (ai + bi) / 2; - ans += (ai * Math.log(ai) + bi * Math.log(bi)) / 2 - half * Math.log(half); + ans += + (a[i] * Math.log(a[i]) + b[i] * Math.log(b[i])) / 2 - + ((a[i] + b[i]) / 2) * Math.log((a[i] + b[i]) / 2); } return ans; } diff --git a/src/distances/jensenShannon.ts b/src/distances/jensenShannon.ts index e625ae7..ededd6e 100644 --- a/src/distances/jensenShannon.ts +++ b/src/distances/jensenShannon.ts @@ -9,11 +9,8 @@ export function jensenShannon(a: NumberArray, b: NumberArray): number { let p = 0; let q = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - const sum = ai + bi; - p += ai * Math.log((2 * ai) / sum); - q += bi * Math.log((2 * bi) / sum); + p += a[i] * Math.log((2 * a[i]) / (a[i] + b[i])); + q += b[i] * Math.log((2 * b[i]) / (a[i] + b[i])); } return (p + q) / 2; } diff --git a/src/distances/kdivergence.ts b/src/distances/kdivergence.ts index 6dbf995..915dcb9 100644 --- a/src/distances/kdivergence.ts +++ b/src/distances/kdivergence.ts @@ -8,8 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function kdivergence(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - ans += ai * Math.log((2 * ai) / (ai + b[i])); + ans += a[i] * Math.log((2 * a[i]) / (a[i] + b[i])); } return ans; } diff --git a/src/distances/kulczynski.ts b/src/distances/kulczynski.ts index f8c2d5f..7e2fe67 100644 --- a/src/distances/kulczynski.ts +++ b/src/distances/kulczynski.ts @@ -9,10 +9,8 @@ export function kulczynski(a: NumberArray, b: NumberArray): number { let up = 0; let down = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - up += Math.abs(ai - bi); - down += Math.min(ai, bi); + up += Math.abs(a[i] - b[i]); + down += Math.min(a[i], b[i]); } return up / down; } diff --git a/src/distances/kullbackLeibler.ts b/src/distances/kullbackLeibler.ts index e5dac0d..0df38bb 100644 --- a/src/distances/kullbackLeibler.ts +++ b/src/distances/kullbackLeibler.ts @@ -8,8 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function kullbackLeibler(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - ans += ai * Math.log(ai / b[i]); + ans += a[i] * Math.log(a[i] / b[i]); } return ans; } diff --git a/src/distances/kumarJohnson.ts b/src/distances/kumarJohnson.ts index a61b103..8435c1a 100644 --- a/src/distances/kumarJohnson.ts +++ b/src/distances/kumarJohnson.ts @@ -8,12 +8,10 @@ import type { NumberArray } from 'cheminfo-types'; export function kumarJohnson(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - const numerator = ai * ai - bi * bi; + const numerator = a[i] * a[i] - b[i] * b[i]; // `prod * Math.sqrt(prod)` is ~6x faster than `prod ** 1.5`, which no // engine specializes; it costs at most 1 ulp of accuracy - const prod = ai * bi; + const prod = a[i] * b[i]; ans += (numerator * numerator) / (2 * prod * Math.sqrt(prod)); } return ans; diff --git a/src/distances/motyka.ts b/src/distances/motyka.ts index 1db4325..b7f1b40 100644 --- a/src/distances/motyka.ts +++ b/src/distances/motyka.ts @@ -9,10 +9,8 @@ export function motyka(a: NumberArray, b: NumberArray): number { let up = 0; let down = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - up += Math.min(ai, bi); - down += ai + bi; + up += Math.min(a[i], b[i]); + down += a[i] + b[i]; } return 1 - up / down; } diff --git a/src/distances/neyman.ts b/src/distances/neyman.ts index 3e567a0..e406bd1 100644 --- a/src/distances/neyman.ts +++ b/src/distances/neyman.ts @@ -8,9 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function neyman(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const diff = ai - b[i]; - d += (diff * diff) / ai; + d += ((a[i] - b[i]) * (a[i] - b[i])) / a[i]; } return d; } diff --git a/src/distances/pearson.ts b/src/distances/pearson.ts index da4dbe4..60ce11e 100644 --- a/src/distances/pearson.ts +++ b/src/distances/pearson.ts @@ -8,9 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function pearson(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - const bi = b[i]; - const diff = a[i] - bi; - d += (diff * diff) / bi; + d += ((a[i] - b[i]) * (a[i] - b[i])) / b[i]; } return d; } diff --git a/src/distances/probabilisticSymmetric.ts b/src/distances/probabilisticSymmetric.ts index e187255..5d310f0 100644 --- a/src/distances/probabilisticSymmetric.ts +++ b/src/distances/probabilisticSymmetric.ts @@ -8,10 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function probabilisticSymmetric(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - const diff = ai - bi; - d += (diff * diff) / (ai + bi); + d += ((a[i] - b[i]) * (a[i] - b[i])) / (a[i] + b[i]); } return 2 * d; } diff --git a/src/distances/ruzicka.ts b/src/distances/ruzicka.ts index e24d5b5..3554691 100644 --- a/src/distances/ruzicka.ts +++ b/src/distances/ruzicka.ts @@ -9,10 +9,8 @@ export function ruzicka(a: NumberArray, b: NumberArray): number { let up = 0; let down = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - up += Math.min(ai, bi); - down += Math.max(ai, bi); + up += Math.min(a[i], b[i]); + down += Math.max(a[i], b[i]); } return up / down; } diff --git a/src/distances/soergel.ts b/src/distances/soergel.ts index 2995d8e..ef00338 100644 --- a/src/distances/soergel.ts +++ b/src/distances/soergel.ts @@ -10,10 +10,8 @@ export function soergel(a: NumberArray, b: NumberArray): number { let up = 0; let down = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - up += Math.abs(ai - bi); - down += Math.max(ai, bi); + up += Math.abs(a[i] - b[i]); + down += Math.max(a[i], b[i]); } return up / down; } diff --git a/src/distances/sorensen.ts b/src/distances/sorensen.ts index 771f9da..c13cb3c 100644 --- a/src/distances/sorensen.ts +++ b/src/distances/sorensen.ts @@ -10,10 +10,8 @@ export function sorensen(a: NumberArray, b: NumberArray): number { let up = 0; let down = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - up += Math.abs(ai - bi); - down += ai + bi; + up += Math.abs(a[i] - b[i]); + down += a[i] + b[i]; } return up / down; } diff --git a/src/distances/squared.ts b/src/distances/squared.ts index c9cba2a..3565dae 100644 --- a/src/distances/squared.ts +++ b/src/distances/squared.ts @@ -8,10 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function squared(a: NumberArray, b: NumberArray): number { let d = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - const diff = ai - bi; - d += (diff * diff) / (ai + bi); + d += ((a[i] - b[i]) * (a[i] - b[i])) / (a[i] + b[i]); } return d; } diff --git a/src/distances/squaredChord.ts b/src/distances/squaredChord.ts index d613f4c..b62e45d 100644 --- a/src/distances/squaredChord.ts +++ b/src/distances/squaredChord.ts @@ -8,8 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function squaredChord(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - const diff = Math.sqrt(a[i]) - Math.sqrt(b[i]); - ans += diff * diff; + ans += (Math.sqrt(a[i]) - Math.sqrt(b[i])) ** 2; } return ans; } diff --git a/src/distances/taneja.ts b/src/distances/taneja.ts index 8c9dacf..12a0df3 100644 --- a/src/distances/taneja.ts +++ b/src/distances/taneja.ts @@ -8,10 +8,9 @@ import type { NumberArray } from 'cheminfo-types'; export function taneja(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - const sum = ai + bi; - ans += (sum / 2) * Math.log(sum / (2 * Math.sqrt(ai * bi))); + ans += + ((a[i] + b[i]) / 2) * + Math.log((a[i] + b[i]) / (2 * Math.sqrt(a[i] * b[i]))); } return ans; } diff --git a/src/distances/tanimoto.ts b/src/distances/tanimoto.ts index d48336d..a0b1260 100644 --- a/src/distances/tanimoto.ts +++ b/src/distances/tanimoto.ts @@ -20,11 +20,9 @@ export function tanimoto( let q = 0; let m = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - p += ai; - q += bi; - m += Math.min(ai, bi); + p += a[i]; + q += b[i]; + m += Math.min(a[i], b[i]); } return (p + q - 2 * m) / (p + q - m); } diff --git a/src/distances/topsoe.ts b/src/distances/topsoe.ts index 126e675..3df7cd4 100644 --- a/src/distances/topsoe.ts +++ b/src/distances/topsoe.ts @@ -8,10 +8,9 @@ import type { NumberArray } from 'cheminfo-types'; export function topsoe(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - const sum = ai + bi; - ans += ai * Math.log((2 * ai) / sum) + bi * Math.log((2 * bi) / sum); + ans += + a[i] * Math.log((2 * a[i]) / (a[i] + b[i])) + + b[i] * Math.log((2 * b[i]) / (a[i] + b[i])); } return ans; } diff --git a/src/distances/waveHedges.ts b/src/distances/waveHedges.ts index 108c0be..86ceae1 100644 --- a/src/distances/waveHedges.ts +++ b/src/distances/waveHedges.ts @@ -8,9 +8,7 @@ import type { NumberArray } from 'cheminfo-types'; export function waveHedges(a: NumberArray, b: NumberArray): number { let ans = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - ans += 1 - Math.min(ai, bi) / Math.max(ai, bi); + ans += 1 - Math.min(a[i], b[i]) / Math.max(a[i], b[i]); } return ans; } diff --git a/src/similarities/cosine.ts b/src/similarities/cosine.ts index 0bedbac..3857896 100644 --- a/src/similarities/cosine.ts +++ b/src/similarities/cosine.ts @@ -9,11 +9,9 @@ export function cosine(a: NumberArray, b: NumberArray): number { let p2 = 0; let q2 = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - p += ai * bi; - p2 += ai * ai; - q2 += bi * bi; + p += a[i] * b[i]; + p2 += a[i] * a[i]; + q2 += b[i] * b[i]; } return p / (Math.sqrt(p2) * Math.sqrt(q2)); } diff --git a/src/similarities/czekanowski.ts b/src/similarities/czekanowski.ts index 3dc6fee..7ca8039 100644 --- a/src/similarities/czekanowski.ts +++ b/src/similarities/czekanowski.ts @@ -9,10 +9,8 @@ export function czekanowski(a: NumberArray, b: NumberArray): number { let up = 0; let down = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - up += Math.min(ai, bi); - down += ai + bi; + up += Math.min(a[i], b[i]); + down += a[i] + b[i]; } return (2 * up) / down; } diff --git a/src/similarities/kumarHassebrook.ts b/src/similarities/kumarHassebrook.ts index fd0a8e6..c2bba35 100644 --- a/src/similarities/kumarHassebrook.ts +++ b/src/similarities/kumarHassebrook.ts @@ -10,11 +10,9 @@ export function kumarHassebrook(a: NumberArray, b: NumberArray): number { let p2 = 0; let q2 = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - p += ai * bi; - p2 += ai * ai; - q2 += bi * bi; + p += a[i] * b[i]; + p2 += a[i] * a[i]; + q2 += b[i] * b[i]; } return p / (p2 + q2 - p); } diff --git a/src/similarities/tanimoto.ts b/src/similarities/tanimoto.ts index 11439ff..371d1c8 100644 --- a/src/similarities/tanimoto.ts +++ b/src/similarities/tanimoto.ts @@ -27,11 +27,9 @@ export function tanimoto( let q = 0; let m = 0; for (let i = 0; i < a.length; i++) { - const ai = a[i]; - const bi = b[i]; - p += ai; - q += bi; - m += Math.min(ai, bi); + p += a[i]; + q += b[i]; + m += Math.min(a[i], b[i]); } return 1 - (p + q - 2 * m) / (p + q - m); }