From a38b3d4fd39b3f7a2743380778506781c19cea12 Mon Sep 17 00:00:00 2001 From: Lumi Pakkanen Date: Fri, 28 Jun 2024 20:59:50 +0300 Subject: [PATCH] Replace tune2, tune3 and tune4 with a single search over integer multiples --- documentation/BUILTIN.md | 12 +-- documentation/tempering.md | 7 +- examples/keen12.sw | 2 +- src/__tests__/temper.spec.ts | 25 ++++- src/parser/__tests__/expression.spec.ts | 37 +++++++ src/parser/__tests__/stdlib.spec.ts | 41 +------- src/stdlib/builtin.ts | 57 ++++++++++- src/stdlib/prelude.ts | 127 ------------------------ src/temper.ts | 65 +++++++++++- 9 files changed, 188 insertions(+), 185 deletions(-) diff --git a/documentation/BUILTIN.md b/documentation/BUILTIN.md index 49672887..cb0fcc77 100644 --- a/documentation/BUILTIN.md +++ b/documentation/BUILTIN.md @@ -411,6 +411,9 @@ Transpose a matrix. For modal transposition see rotate(). ### trunc(*interval*) Truncate value towards zero to the nearest integer. +### tune(*vals*, *searchRadius*, *weights*) +Attempt to combine the given vals into a more Tenney-Euclid optimal val. Weights are applied multiplicatively on top of Tenney weights of the subgroup basis. + ### unshift(*interval*, *scale = $$*) Prepend an interval at the beginning of the current/given scale. @@ -728,15 +731,6 @@ Obtain a copy of the current/given scale quantized to subharmonics of the given ### trap(*message*) Produce a function that fails with the given message when called. -### tune2(*a*, *b*, *numIter = 1*, *weights = niente*) -Find a combination of two vals that is closer to just intonation. - -### tune3(*a*, *b*, *c*, *numIter = 1*, *weights = niente*) -Find a combination of three vals that is closer to just intonation. - -### tune4(*a*, *b*, *c*, *d*, *numIter = 1*, *weights = niente*) -Find a combination of four vals that is closer to just intonation. - ### u(*scale = ££*) Obtain a undertonal reflection of the popped/given overtonal scale. diff --git a/documentation/tempering.md b/documentation/tempering.md index 323efc34..66dc2609 100644 --- a/documentation/tempering.md +++ b/documentation/tempering.md @@ -267,11 +267,10 @@ const S = @2.7.11 ### errorTE In order to compare the quality of vals w.r.t. just intonation SonicWeave provides the `errorTE` helper that measures [RMS TE error](https://en.xen.wiki/w/Tenney-Euclidean_temperament_measures#TE_error) in cents. Providing explicit subgroups is highly recommended so that irrelevant higher primes do not interfere with the measure. -### tune2 -The helper `tune2` takes two vals and tries to find a combination that's closer to just intonation, effectively performing constrained Tenney-Euclidean optimization (CTE). E.g. `tune2(12@.5, 19@.5)` finds 31p. Given more iterations `tune2(5@.5, 7@.5, 7)` finds an equal temperament in the thousands that's virtually indistinguishable from the true meantone CTE tuning. +### tune +The helper `tune` takes an array of vals and tries to find a combination that's closer to just intonation, effectively performing constrained Tenney-Euclidean optimization (CTE). E.g. `tune([12@.5, 19@.5])` finds 31p. Given a large search radius `tune2(5@.5, 7@.5, 200)` finds an equal temperament in the thousands that's virtually indistinguishable from the true meantone CTE tuning. -### tune3 and tune4 -The helpers `tune3` and `tune4` do the same but try combinations of 3 or 4 vals instead. Assuming the vals are linearly independent the process corresponds to rank-3 and rank-4 CTE optimization. +With more than two linearly independent vals the process corresponds to higher rank CTE optimization. ### Discovering vals Every just intonation subgroup has a generalized patent val sequence that dances around approximations to the just intonation point i.e. the perfectly pure tuning. diff --git a/examples/keen12.sw b/examples/keen12.sw index 4aa021ba..e942ce29 100644 --- a/examples/keen12.sw +++ b/examples/keen12.sw @@ -18,4 +18,4 @@ labelAbsoluteFJS (* See https://en.xen.wiki/w/Val#Sparse_Offset_Val_notation for the alternative notation for 12d. *) (* Doesn't actually matter here because we're stacking fifths, but communicates the intended interpretation. *) (* Same as 56@. *) -tune(12[v7]@, 22@, 2) +tune([12[v7]@, 22@], 2) diff --git a/src/__tests__/temper.spec.ts b/src/__tests__/temper.spec.ts index f57f1fe2..089d45fb 100644 --- a/src/__tests__/temper.spec.ts +++ b/src/__tests__/temper.spec.ts @@ -1,8 +1,13 @@ import {describe, expect, it} from 'vitest'; -import {TuningMap, combineTuningMaps, vanishCommas} from '../temper'; +import { + TuningMap, + combineTuningMaps, + intCombineTuningMaps, + vanishCommas, +} from '../temper'; import {LOG_PRIMES, applyWeights, dot, unapplyWeights} from 'xen-dev-utils'; -describe('Val combination optimizer', () => { +describe('Mapping combination optimizer', () => { it('computes TE meantone', () => { const p12: TuningMap = [12, 19, 28].map(c => (c / 12) * LOG_PRIMES[0]); const p19: TuningMap = [19, 30, 44].map(c => (c / 19) * LOG_PRIMES[0]); @@ -47,6 +52,22 @@ describe('Val combination optimizer', () => { }); }); +describe('Val combination search', () => { + it('combines 12p and 19p into 31p', () => { + const p12: TuningMap = [12, 19 / Math.log2(3), 28 / Math.log2(5)]; + const p19: TuningMap = [19, 30 / Math.log2(3), 44 / Math.log2(5)]; + const coeffs = intCombineTuningMaps([1, 1, 1], [p12, p19], 1); + expect(coeffs).toEqual([1, 1]); + }); + + it('combines 5p and 7p into 31p', () => { + const p5: TuningMap = [5, 8 / Math.log2(3), 12 / Math.log2(5)]; + const p7: TuningMap = [7, 11 / Math.log2(3), 16 / Math.log2(5)]; + const coeffs = intCombineTuningMaps([1, 1, 1], [p5, p7], 3); + expect(coeffs).toEqual([2, 3]); + }); +}); + describe('Comma vanisher', () => { it('computes TE meantone', () => { const syntonic = [-4, 4, -1]; diff --git a/src/parser/__tests__/expression.spec.ts b/src/parser/__tests__/expression.spec.ts index 80be5ba7..a1735294 100644 --- a/src/parser/__tests__/expression.spec.ts +++ b/src/parser/__tests__/expression.spec.ts @@ -2703,4 +2703,41 @@ describe('SonicWeave expression evaluator', () => { ); expect(fraction).toBe('11/9'); }); + + it('can combine two vals to approach the JIP', () => { + const thirtyOne = evaluate('tune([12@.5, 19@.5])') as Val; + expect(thirtyOne.value.toIntegerMonzo()).toEqual([31, 49, 72]); + }); + + it('can combine two vals to approach the JIP (Wilson metric)', () => { + const eighty = evaluate( + 'tune([12@.5, 22@.5], 5, [1 % 2, 3 /_ 2 % 3, 5 /_ 2 % 5])' + ) as Val; + expect(eighty.value.toIntegerMonzo()).toEqual([126, 200, 293]); + }); + + it('can combine three vals to approach the JIP', () => { + const fourtyOne = evaluate('tune([5@.7, 17@.7, 19@.7])') as Val; + expect(fourtyOne.value.toIntegerMonzo()).toEqual([41, 65, 95, 115]); + }); + + it('can combine four vals to approach the JIP', () => { + const val = evaluate('tune([5@.11, 17@.11, 19@.11, 31@.11])') as Val; + expect(val.value.toIntegerMonzo()).toEqual([72, 114, 167, 202, 249]); + }); + + it("doesn't move from 31p when tuned with 5p", () => { + const p31 = evaluate('str(tune([31@.5, 5@.5]))'); + expect(p31).toBe('<31 49 72]'); + }); + + it('moves from 31p when tuned with 5p if given a large enough radius', () => { + const p31 = evaluateExpression('str(tune([31@.5, 5@.5], 6))'); + expect(p31).toBe('<191 302 444]'); + }); + + it('tunes close to CTE meantone if given a large enough radius', () => { + const cents = evaluate('str(cents(3/2 tmpr tune([5@.5, 7@.5], 200), 4))'); + expect(cents).toEqual('697.2143'); + }); }); diff --git a/src/parser/__tests__/stdlib.spec.ts b/src/parser/__tests__/stdlib.spec.ts index 4e10a5f0..9e1636ec 100644 --- a/src/parser/__tests__/stdlib.spec.ts +++ b/src/parser/__tests__/stdlib.spec.ts @@ -6,7 +6,7 @@ import { getSourceVisitor, parseAST, } from '../../parser'; -import {Interval, Val} from '../../interval'; +import {Interval} from '../../interval'; import {builtinNode, track} from '../../stdlib'; import {Fraction} from 'xen-dev-utils'; @@ -366,30 +366,6 @@ describe('SonicWeave standard library', () => { ); }); - it('can combine two vals to approach the JIP', () => { - const thirtyOne = evaluateExpression('tune2(12@.5, 19@.5)') as Val; - expect(thirtyOne.value.toIntegerMonzo()).toEqual([31, 49, 72]); - }); - - it('can combine two vals to approach the JIP (Wilson metric)', () => { - const eighty = evaluateExpression( - 'tune2(12@.5, 22@.5, 2, [log(2)/2, log(3)/3, log(5)/5])' - ) as Val; - expect(eighty.value.toIntegerMonzo()).toEqual([126, 200, 293]); - }); - - it('can combine three vals to approach the JIP', () => { - const fourtyOne = evaluateExpression('tune3(5@.7, 17@.7, 19@.7)') as Val; - expect(fourtyOne.value.toIntegerMonzo()).toEqual([41, 65, 95, 115]); - }); - - it('can combine four vals to approach the JIP', () => { - const val = evaluateExpression( - 'tune4(5@.11, 17@.11, 19@.11, 31@.11)' - ) as Val; - expect(val.value.toIntegerMonzo()).toEqual([94, 149, 218, 264, 325]); - }); - // Remember that unison (0.0 c) is implicit in SonicWeave it('has harmonic segments sounding downwards', () => { @@ -1705,19 +1681,4 @@ describe('SonicWeave standard library', () => { '[6b@2.3.5.7, 10@2.3.5.7, 16@2.3.5.7, 20@2.3.5.7, 26@2.3.5.7]' ); }); - - it("doesn't move from 31p when tuned with 5p", () => { - const p31 = evaluateExpression('str(tune2(31@.5, 5@.5))'); - expect(p31).toBe('<31 49 72]'); - }); - - it('moves from 31p when tuned with 5p if given enough time', () => { - const p31 = evaluateExpression('str(tune2(31@.5, 5@.5, 3))'); - expect(p31).toBe('<191 302 444]'); - }); - - it('tunes close to CTE meantone if given enough time', () => { - const scale = expand('cents(3/2 tmpr tune2(5@.5, 7@.5, 7), 4)'); - expect(scale).toEqual(['697.2143']); - }); }); diff --git a/src/stdlib/builtin.ts b/src/stdlib/builtin.ts index 4d10ab4a..da62fc76 100644 --- a/src/stdlib/builtin.ts +++ b/src/stdlib/builtin.ts @@ -80,7 +80,12 @@ import { } from './public'; import {scaleMonzos} from '../diamond-mos'; import {valToSparseOffset, valToWarts} from '../warts'; -import {TuningMap, combineTuningMaps, vanishCommas} from '../temper'; +import { + TuningMap, + combineTuningMaps, + intCombineTuningMaps, + vanishCommas, +} from '../temper'; const {version: VERSION} = require('../../package.json'); // === Library === @@ -1770,6 +1775,55 @@ function nextGPV( nextGPV.__doc__ = 'Obtain the next generalized patent val in the sequence.'; nextGPV.__node__ = builtinNode(nextGPV); +function tune( + this: ExpressionVisitor, + vals: SonicWeaveValue, + searchRadius: SonicWeaveValue, + weights: SonicWeaveValue +) { + if (!Array.isArray(vals)) { + throw new Error('An array of vals is required.'); + } + if (!vals.length) { + throw new Error('At least one val is required.'); + } + for (const val of vals) { + if (!(val instanceof Val)) { + throw new Error('An array of vals is required.'); + } + } + const vs = vals as Val[]; + const basis = vs[0].basis; + for (const val of vs.slice(1)) { + if (!val.basis.equals(basis)) { + throw new Error('Vals bases must agree in tuning.'); + } + } + const radius = + searchRadius === undefined ? 1 : upcastBool(searchRadius).toInteger(); + this.spendGas((2 * radius + 1) ** vals.length); + const jip = basis.value.map(m => m.totalCents()); + const ws = valWeights(weights, basis.size); + const wvals = vs.map(val => + applyWeights( + unapplyWeights( + basis.value.map(m => m.dot(val.value).valueOf()), + jip + ), + ws + ) + ); + const coeffs = intCombineTuningMaps(ws, wvals, radius); + let result = vs[0].mul(fromInteger(coeffs[0])); + for (let i = 1; i < coeffs.length; ++i) { + result = result.add(vs[i].mul(fromInteger(coeffs[i]))); + } + return result.abs(); +} +tune.__doc__ = + 'Attempt to combine the given vals into a more Tenney-Euclid optimal val. Weights are applied multiplicatively on top of Tenney weights of the subgroup basis.'; +tune.__node__ = builtinNode(tune); + function tenneyHeight( this: ExpressionVisitor, interval: SonicWeaveValue @@ -3080,6 +3134,7 @@ export const BUILTIN_CONTEXT: Record = { TE, errorTE, nextGPV, + tune, tenneyHeight, wilsonHeight, respell, diff --git a/src/stdlib/prelude.ts b/src/stdlib/prelude.ts index 49ceb6b7..7540d449 100644 --- a/src/stdlib/prelude.ts +++ b/src/stdlib/prelude.ts @@ -333,133 +333,6 @@ riff enumerate(array = $$) { return [[i, array[i]] for i in array]; } -riff tune2(a, b, numIter = 1, weights = niente) { - "Find a combination of two vals that is closer to just intonation."; - - let error = (v => errorTE(v, weights, true)); - if (isFunction(weights)) { - error = weights; - } - - numIter = real(numIter) + 1r; - while (--numIter) { - const combos = [ - a, - b, - a + b, - 2*a + b, - 2*b + a, - - abs (a - b), - abs (2*a - b), - abs (2*b - a), - ]; - - [a, b] = sort(combos, (u, v) => error(u) ~- error(v)); - } - return abs a; -} - -riff tune3(a, b, c, numIter = 1, weights = niente) { - "Find a combination of three vals that is closer to just intonation."; - - let error = (v => errorTE(v, weights, true)); - if (isFunction(weights)) { - error = weights; - } - - numIter = real(numIter) + 1r; - while (--numIter) { - const combos = [ - (* Corners *) - a, - b, - c, - - (* Edge midpoints *) - a + b, - a + c, - b + c, - - (* Midpoint *) - a + b + c, - - (* Corner extensions *) - 2 * a + b + c, - 2 * b + a + c, - 2 * c + a + b, - - (* Edge extensions *) - 2 * (a + b) + c, - 2 * (a + c) + b, - 2 * (b + c) + a, - ]; - - [a, b, c] = sort(combos, (u, v) => error(u) ~- error(v)); - } - return abs a; -} - -riff tune4(a, b, c, d, numIter = 1, weights = niente) { - "Find a combination of four vals that is closer to just intonation."; - - let error = (v => errorTE(v, weights, true)); - if (isFunction(weights)) { - error = weights; - } - - numIter = real(numIter) + 1r; - while (--numIter) { - const combos = [ - (* Corners *) - a, - b, - c, - d, - - (* Edge midpoints *) - a + b, - a + c, - a + d, - b + c, - b + d, - c + d, - - (* Face midpoints *) - a + b + c, - a + b + d, - a + c + d, - b + c + d, - - (* Midpoint *) - a + b + c + d, - - (* Corner extensions *) - 2 * a + b + c + d, - 2 * b + a + c + d, - 2 * c + a + b + d, - 2 * d + a + b + c, - - (* Edge extensions *) - 2 * (a + b) + c + d, - 2 * (a + c) + b + d, - 2 * (a + d) + b + c, - 2 * (b + c) + a + b, - 2 * (b + d) + a + c, - 2 * (c + d) + a + b, - - (* Face extensions *) - 2 * (a + b + c) + d, - 2 * (a + b + d) + c, - 2 * (a + c + d) + b, - 2 * (b + c + d) + a, - ]; - - [a, b, c, d] = sort(combos, (u, v) => error(u) ~- error(v)); - } - return abs a; -} - riff supportingGPVs(initialVal, commas, count = 5, weights = niente, maxIter = 1000) { "Obtain generalized patent vals in the same sequence as the initial val that make the given commas vanish."; if (not isArray(commas)) { diff --git a/src/temper.ts b/src/temper.ts index 0dd849bc..6a3a5b09 100644 --- a/src/temper.ts +++ b/src/temper.ts @@ -1,4 +1,14 @@ -import {Monzo, inv, matmul, sub, transpose} from 'xen-dev-utils'; +import { + Monzo, + add, + gcd, + inv, + matmul, + norm, + scale, + sub, + transpose, +} from 'xen-dev-utils'; export type TuningMap = number[]; @@ -28,3 +38,56 @@ export function vanishCommas(jip: TuningMap, commas: Monzo[]): TuningMap { const nullProjector = matmul(pinv, commas); return sub(jip, matmul(nullProjector, jip)); } + +/** + * Combine tuning maps to minimize the distance to the JIP. + * @param jip Just Intonation Point in (co-)weighted coordinates. + * @param vals Vals to combine in (co-)weighted coordinates. + * @param searchRadius Width of the search space. + * @returns Integer coefficients of the linear combination closest to the JIP. + */ +export function intCombineTuningMaps( + jip: TuningMap, + vals: Monzo[], + searchRadius: number +): number[] { + if (!vals.length) { + return []; + } + if (vals.length === 1) { + return [1]; + } + let leastError = Infinity; + let result: number[] = []; + function combine(coeffs: number[]) { + if (coeffs.length < vals.length) { + for (let i = -searchRadius; i <= searchRadius; ++i) { + const newCoeffs = [...coeffs]; + newCoeffs.push(i); + combine(newCoeffs); + } + } else { + if (Math.abs(coeffs.reduce(gcd)) !== 1) { + return; + } + let val = scale(vals[0], coeffs[0]); + for (let i = 1; i < coeffs.length; ++i) { + val = add(val, scale(vals[i], coeffs[i])); + } + if (!val[0]) { + return; + } + const map = val.map(v => (v / val[0]) * jip[0]); + const error = norm(sub(jip, map), 'L2'); + if (error < leastError) { + leastError = error; + result = coeffs; + } + } + } + // The search space is projective so we only need to search half of it compared to euclidean space. + for (let i = 0; i <= searchRadius; ++i) { + combine([i]); + } + return result; +}