From eacd314a7a055b84fe3baae1222d2a65e23c4003 Mon Sep 17 00:00:00 2001
From: Lumi Pakkanen <lumi.pakkanen@gmail.com>
Date: Sun, 3 Mar 2024 20:24:18 +0200
Subject: [PATCH] WIP: Delta rationalize chord

ref #133
---
 AUTHORS                     |  5 ++
 src/__tests__/utils.spec.ts | 10 ++++
 src/utils.ts                | 95 +++++++++++++++++++++++++++++++++++++
 3 files changed, 110 insertions(+)
 create mode 100644 AUTHORS
 create mode 100644 src/__tests__/utils.spec.ts

diff --git a/AUTHORS b/AUTHORS
new file mode 100644
index 00000000..5fff1b81
--- /dev/null
+++ b/AUTHORS
@@ -0,0 +1,5 @@
+# Lead developer
+Lumi Pakkanen <lumi.pakkanen@gmail.com> (https://github.com/frostburn/)
+
+# Contributors
+Akselai <underscore.taroumaru@gmail.com> (https://github.com/akselai/)
diff --git a/src/__tests__/utils.spec.ts b/src/__tests__/utils.spec.ts
new file mode 100644
index 00000000..3827b6db
--- /dev/null
+++ b/src/__tests__/utils.spec.ts
@@ -0,0 +1,10 @@
+import {describe, it, expect} from 'vitest';
+import {deltaRationalize} from '../utils';
+
+describe('Chord delta rationalizer', () => {
+  it('delta rationalizes with max root = 4', () => {
+    const result = deltaRationalize([1, 1.26, 1.498], 4);
+    expect(result.error).toBeLessThan(0.008);
+    expect(result.signature).toHaveLength(2); // It's broken, plz fix.
+  });
+});
diff --git a/src/utils.ts b/src/utils.ts
index 33b649a3..d55069f4 100644
--- a/src/utils.ts
+++ b/src/utils.ts
@@ -1,5 +1,7 @@
 import {Fraction, PRIMES} from 'xen-dev-utils';
 
+const EPSILON = 1e-5;
+
 export const ZERO = new Fraction(0);
 export const ONE = new Fraction(1);
 export const NEGATIVE_ONE = new Fraction(-1);
@@ -214,3 +216,96 @@ export function validateBigInt(n: bigint) {
     throw new Error('Integer overflow.');
   }
 }
+
+/**
+ * Return the index of the minimum value of the array.
+ * @param array Array of values to compare.
+ * @returns The index of the minimum value.
+ */
+export function argMin(array: number[]) {
+  if (!array.length) {
+    return NaN;
+  }
+  let indexOfMinimum = 0;
+  for (let i = 1; i < array.length; i++) {
+    if (array[i] < array[indexOfMinimum]) {
+      indexOfMinimum = i;
+    }
+  }
+  return indexOfMinimum;
+}
+
+export type MinimizationResult = {
+  x: number;
+  value: number;
+};
+
+export function minimizeFunction(
+  fn: (x: number) => number,
+  minX: number,
+  maxX: number,
+  numPartitions = 40
+): MinimizationResult {
+  if (maxX - minX < EPSILON) {
+    return {
+      x: (minX + maxX) / 2,
+      value: fn((minX + maxX) / 2),
+    };
+  }
+  const values = [];
+  for (let i = 0; i <= numPartitions; i++) {
+    const h = fn(minX + ((maxX - minX) / numPartitions) * i);
+    values.push(h);
+  }
+  const m = argMin(values);
+  return minimizeFunction(
+    fn,
+    minX + ((maxX - minX) / numPartitions) * (m - 1),
+    minX + ((maxX - minX) / numPartitions) * (m + 1),
+    5
+  );
+}
+// the objective function to be minimized
+function deltaRationalObjective(
+  chord: number[],
+  maxRoot: number,
+  offset: number
+) {
+  let signature: number[] = [];
+  let minError = Infinity;
+  for (let l = 1; l <= maxRoot; l++) {
+    const jiChord = [];
+    for (const x of chord) {
+      const b = chord[0] / x;
+      // utonal version: x / chord[chord.length - 1];
+      jiChord.push(Math.round(l / b));
+    }
+    let maxDiff = 0;
+    const root = jiChord[0] + offset;
+    for (let i = 0; i < chord.length; i++) {
+      const d = Math.abs((jiChord[i] + offset) / root - chord[i]);
+      if (d > maxDiff) {
+        maxDiff = d;
+      }
+    }
+    if (maxDiff < minError) {
+      signature = jiChord;
+      minError = maxDiff;
+    }
+  }
+  return {
+    signature,
+    error: minError,
+  };
+}
+
+export function deltaRationalize(chord: number[], maxRoot: number) {
+  const f = (x: number) => deltaRationalObjective(chord, maxRoot, x).error;
+  const minimum = minimizeFunction(f, -0.5, 0.5);
+  const offset = minimum.x;
+  return {
+    ...deltaRationalObjective(chord, maxRoot, offset),
+    offset,
+    chord: chord.map(c => c + offset),
+  };
+}