Implement algorithms for finding concordance shells

ref #471

src/__tests__/analysis.spec.ts

@@ -1,7 +1,16 @@
 import { describe, it, expect } from "vitest";
-import { arraysEqual } from "xen-dev-utils";
+import { arraysEqual, valueToCents } from "xen-dev-utils";
 
-import { otonalFundamental, utonalFundamental } from "../analysis";
+import {
+  alignValues,
+  concordanceShell,
+  minimaxConcordanceShells,
+  misalignment,
+  otonalFundamental,
+  utonalFundamental,
+} from "../analysis";
+
+const EPSILON = 1e-4;
 
 describe("Otonal balancer", () => {
   it("can figure out that the major chord in 12edo approximates 4:5:6", () => {
@@ -28,3 +37,73 @@ describe("Utonal balancer", () => {
     ).toBeTruthy();
   });
 });
+
+describe("Equal-division deviation minimizer", () => {
+  it("can figure out the optimal alignment of 4:5:6 on 12edo", () => {
+    const minimumAchievableError = alignValues([4, 5, 6], 100.0).error;
+    expect(minimumAchievableError).closeTo(7.8206, EPSILON);
+
+    // Attempt (and fail) to find a counter-example
+    const pitches = [4, 5, 6].map(valueToCents);
+    for (let i = 0; i < 100; ++i) {
+      const offset = Math.random() * 1200;
+      const candidate = pitches.map((pitch) => pitch + offset);
+      expect(misalignment(candidate, 100.0)).toBeGreaterThanOrEqual(
+        minimumAchievableError
+      );
+    }
+  });
+});
+
+describe("Concordance shell finder", () => {
+  it("finds the 37 : 44 : 66 : 70 : 83 : 93 : 99 : 111 : 157 : 187 : 197 : 209 : 235 : 249 shell in 12edo", () => {
+    const shell = concordanceShell(37, 255, 5.0, 100.0);
+    expect(shell).toEqual([
+      37, 44, 66, 70, 83, 93, 99, 111, 157, 187, 197, 209, 235, 249,
+    ]);
+  });
+});
+
+describe("Minimaxing concordance shell finder", () => {
+  it("finds more variants for the shell rooted at 37 in 12edo", () => {
+    const shells = minimaxConcordanceShells(37, 127, 3, 100.0, 4);
+
+    // Just some way to create a deterministic order.
+    function cmp(a: number[], b: number[]) {
+      if (!a.length && b.length) {
+        return -1;
+      } else if (a.length && !b.length) {
+        return +1;
+      } else if (!a.length && !b.length) {
+        return 0;
+      }
+      if (a[0] < b[0]) {
+        return -1;
+      } else if (a[0] > b[0]) {
+        return +1;
+      }
+      return cmp(a.slice(1), b.slice(1));
+    }
+    shells.sort((a, b) => cmp(a.harmonics, b.harmonics));
+    for (const shell of shells) {
+      shell.error = Math.round(shell.error * 1000) / 1000;
+    }
+    expect(shells).toEqual([
+      {
+        harmonics: [37, 44, 66, 70, 83, 99, 111],
+        error: 2.59,
+        gridPitches: [0, 300, 1000, 1100, 1400, 1700, 1900],
+      },
+      {
+        harmonics: [37, 44, 66, 70, 99, 105, 111], // The harmonic 105 is novel here compared to the rooted error model.
+        error: 2.888,
+        gridPitches: [0, 300, 1000, 1100, 1700, 1800, 1900],
+      },
+      {
+        harmonics: [37, 44, 83, 93],
+        error: 2.177,
+        gridPitches: [0, 300, 1400, 1600],
+      },
+    ]);
+  });
+});

src/analysis.ts

@@ -1,4 +1,5 @@
 import type { Scale } from "scale-workshop-core";
+import { arraysEqual, valueToCents } from "xen-dev-utils";
 
 const EPSILON = 1e-6;
 
@@ -112,3 +113,226 @@ export function intervalMatrix(scale: Scale) {
   }
   return result;
 }
+
+/**
+ * Calculate the maximum deviation of a set of pitches from an equal grid in pitch-space.
+ * @param pitches Array of pitches measured in cents.
+ * @param gridCents The distance in cents between two grid lines of the equal division e.g. `100.0`.
+ * @returns The maximum distance in cents from a grid "line" in the set.
+ */
+export function misalignment(pitches: number[], gridCents: number) {
+  const gridPitches = pitches.map(
+    (pitch) => Math.round(pitch / gridCents) * gridCents
+  );
+  let error = 0;
+  for (let i = 0; i < pitches.length; ++i) {
+    error = Math.max(error, Math.abs(pitches[i] - gridPitches[i]));
+  }
+  return error;
+}
+
+/**
+ * Align a set of pitches on an equal division of pitch-space such that the maximum absolute error is minimized.
+ * @param pitches An array of pitches measured in cents.
+ * @param gridCents The distance in cents between two grid lines of the equal division e.g. `100.0`.
+ * @returns The minimum misalignment achievable measured in cents and the pitches snapped to the grid.
+ */
+export function alignCents(pitches: number[], gridCents: number) {
+  // The error function we're trying to optimize is piecewise linear.
+
+  // Find the segment where the global optimum lies.
+  let optimalPitches: number[];
+  let minError = Infinity;
+  for (let i = 0; i < pitches.length; ++i) {
+    const aligned = pitches.map((pitch) => pitch - pitches[i]);
+    const error = misalignment(aligned, gridCents);
+    if (error < minError) {
+      optimalPitches = aligned;
+      minError = error;
+    }
+  }
+
+  // Calculate the shape of the segment.
+  const gridPitches = optimalPitches!.map(
+    (pitch) => Math.round(pitch / gridCents) * gridCents
+  );
+  let minOffset = Infinity;
+  let maxOffset = -Infinity;
+  let root = Infinity;
+  for (let i = 0; i < gridPitches.length; ++i) {
+    const offset = optimalPitches![i] - gridPitches[i];
+    minOffset = Math.min(minOffset, offset);
+    maxOffset = Math.max(maxOffset, offset);
+    if (gridPitches[i] < root) {
+      root = gridPitches[i];
+    }
+  }
+
+  // Calculate minimum achievable absolute error.
+  const error = 0.5 * Math.abs(minOffset) + 0.5 * Math.abs(maxOffset);
+  // Move root to grid origin.
+  for (let i = 0; i < gridPitches.length; ++i) {
+    gridPitches[i] -= root;
+  }
+  return {
+    error,
+    gridPitches,
+  };
+}
+
+/**
+ * Align a set of ratios in frequency-space on an equal division of pitch-space such that the maximum absolute error is minimized.
+ * @param ratios An array of frequency ratios.
+ * @param gridCents The distance in cents between two grid lines of the equal division e.g. `100.0`.
+ * @returns The minimum misalignment achievable measured in cents and the pitches snapped to the grid.
+ */
+export function alignValues(ratios: number[], gridCents: number) {
+  return alignCents(ratios.map(valueToCents), gridCents);
+}
+
+/**
+ * Find the concordance shell centered around a root harmonic.
+ * @param rootHarmonic Root harmonic.
+ * @param maxHarmonic Maximum harmonic to consider for inclusion.
+ * @param tolerance Maximum deviation in cents from the equal division of pitch-space.
+ * @param gridCents The distance in cents between two grid lines of the equal division e.g. `100.0`.
+ * @param exclusionEquave The interval of equivalence. If non-null equivalent harmonics are filtered out.
+ * @returns An enumerated chord as an array of integers.
+ */
+export function concordanceShell(
+  rootHarmonic: number,
+  maxHarmonic: number,
+  tolerance: number,
+  gridCents: number,
+  exclusionEquave: number | null = 2
+) {
+  const harmonics = [rootHarmonic];
+  search: for (
+    let harmonic = rootHarmonic + 1;
+    harmonic <= maxHarmonic;
+    ++harmonic
+  ) {
+    if (exclusionEquave !== null) {
+      for (const existing of harmonics) {
+        if (harmonic % existing === 0) {
+          let temp = harmonic;
+          while (temp > existing && temp % exclusionEquave === 0) {
+            temp /= exclusionEquave;
+            if (temp === existing) {
+              continue search;
+            }
+          }
+        }
+      }
+    }
+    const pitch = valueToCents(harmonic / rootHarmonic);
+    const gridPitch = Math.round(pitch / gridCents) * gridCents;
+    if (Math.abs(pitch - gridPitch) < tolerance) {
+      harmonics.push(harmonic);
+    }
+  }
+  return harmonics;
+}
+
+function* subshells(shell: number[]) {
+  if (shell.length <= 1) {
+    return;
+  }
+  for (let i = 1; i < shell.length; ++i) {
+    const subshell = [...shell];
+    subshell.splice(i, 1);
+    yield subshell;
+  }
+}
+
+function subsetOf(a: number[], b: number[]) {
+  for (const value of a) {
+    if (!b.includes(value)) {
+      return false;
+    }
+  }
+  return true;
+}
+
+export type Shell = {
+  harmonics: number[];
+  gridPitches: number[];
+  error: number;
+};
+
+/**
+ * Find concordance shells based on the given harmonic.
+ * @param rootHarmonic Root harmonic.
+ * @param maxHarmonic Maximum harmonic to consider for inclusion.
+ * @param tolerance Maximum deviation in cents from the equal division of pitch-space.
+ * @param gridCents The distance in cents between two grid lines of the equal division e.g. `100.0`.
+ * @param minSize Minimum size of the shell to search for.
+ * @param maxCount Maximum number of shells to discover.
+ * @param exclusionEquave The interval of equivalence. If non-null equivalent harmonics are filtered out.
+ * @returns An array of objects with
+ *  harmonics of the shell as an array of integers,
+ *  the harmonics snapped to the grid as an array of cents,
+ *  the minimum achievable misalignment of the harmonics within the grid measured in cents.
+ */
+export function minimaxConcordanceShells(
+  rootHarmonic: number,
+  maxHarmonic: number,
+  tolerance: number,
+  gridCents: number,
+  minSize = 3,
+  maxCount = 100,
+  exclusionEquave: number | null = 2
+): Shell[] {
+  // Calculate the largest possible shell if everything aligns perfectly.
+  const supershell = concordanceShell(
+    rootHarmonic,
+    maxHarmonic,
+    2 * tolerance,
+    gridCents,
+    exclusionEquave
+  );
+
+  if (supershell.length < minSize) {
+    return [];
+  }
+
+  const superAlignment = alignValues(supershell, gridCents);
+  // Not gonna happen, but let's pretend to be hopeful.
+  if (superAlignment.error <= tolerance) {
+    return [{ harmonics: supershell, ...superAlignment }];
+  }
+
+  // Start breaking the super-shell into smaller and smaller sub-shells that fit within the tolerance.
+  const result: Shell[] = [];
+  const badShells = [supershell];
+  while (badShells.length) {
+    const shell = badShells.shift()!; // Break bigger shells first.
+    search: for (const subshell of subshells(shell)) {
+      for (const existing of result) {
+        if (subsetOf(subshell, existing.harmonics)) {
+          continue search;
+        }
+      }
+      for (const existing of badShells) {
+        if (arraysEqual(existing, subshell)) {
+          continue search;
+        }
+      }
+
+      const alignment = alignValues(subshell, gridCents);
+      if (alignment.error <= tolerance) {
+        result.push({
+          harmonics: subshell,
+          ...alignment,
+        });
+        if (result.length >= maxCount) {
+          return result;
+        }
+      } else if (subshell.length > minSize) {
+        badShells.push(subshell);
+      }
+    }
+  }
+
+  return result;
+}

src/views/AboutView.vue

@@ -39,7 +39,8 @@
         Lajos Mészáros - <i>developer</i><br />
         Forrest Cahoon - <i>developer</i> <br />
         Videco - <i>developer</i> <br />
-        Kraig Grady - <i>lattice advisor</i>
+        Kraig Grady - <i>lattice advisor</i> <br />
+        Ben Lampert - <i>developer</i>
       </p>
     </div>
   </div>