Fix bugs in ChordSpace code for porting HarmonyIFS and ChordSpace to CsoundAC in C++ and WebAssembly

gogins commented 4 years ago

I have discovered that my concepts and code for the OPTI and OPTTI symmetries of chord space were wrong, and therefore wrongly implemented. This is to fix that in both JavaScript and C++ and thus make the ChordSpaceGroup class usable as intended.

gogins commented 4 years ago

I might conceivably have fixed the bugs in eX and iseX. Now to fix the ChordSpaceGroup.

gogins commented 4 years ago

I have simplified the ChordSpaceTest.html unit tests.

Tests for OPTTI for tetrachords are now way off, apparently because iseOPTTI produces false positives.

gogins commented 4 years ago

How does eOPTT major seventh chord compare to eOPTT minor seventh chord? The latter should also be isOPTTI. The data implies that eOPT and eOPTT are correct for trichords and incorrect for higher arities.

CM7: [  60.0000000   64.0000000   67.0000000   71.0000000 ]
        Chord:   [  60.0000000   64.0000000   67.0000000   71.0000000 ]  name:     CM7 
        pcs:     [   0.0000000    4.0000000    7.0000000   11.0000000 ]
        pitches: C4 E4 G4 B4
        type:    [   0.0000000    1.0000000    5.0000000    8.0000000 ]
        eO:      [   0.0000000    4.0000000    7.0000000   -1.0000000 ]  iseO:     false
        eP:      [  60.0000000   64.0000000   67.0000000   71.0000000 ]  iseP:     true
        eT:      [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseT:     false
        eTT:     [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseTT:    false
        eI:      [  60.0000000   64.0000000   67.0000000   71.0000000 ]  iseI:     true
        eV:      [  60.0000000   64.0000000   67.0000000   71.0000000 ]  iseV:     true
        eOP:     [  -1.0000000    0.0000000    4.0000000    7.0000000 ]  iseOP:    false
        eOPT:    [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseOPT:   false
        eOPTI:   [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseOPTI:  false
        eOPTT:   [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseOPTT:  false
        eOPTTI:  [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseOPTTI: false
        sum:     262

CM7.eOP():
        Chord:   [  -1.0000000    0.0000000    4.0000000    7.0000000 ]  name:     CM7 
        pcs:     [   0.0000000    4.0000000    7.0000000   11.0000000 ]
        pitches: B-2 C-1 E-1 G-1
        type:    [   0.0000000    1.0000000    5.0000000    8.0000000 ]
        eO:      [  -1.0000000    0.0000000    4.0000000    7.0000000 ]  iseO:     true
        eP:      [  -1.0000000    0.0000000    4.0000000    7.0000000 ]  iseP:     true
        eT:      [  -3.5000000   -2.5000000    1.5000000    4.5000000 ]  iseT:     false
        eTT:     [  -3.0000000   -2.0000000    2.0000000    5.0000000 ]  iseTT:    false
        eI:      [  -1.0000000    0.0000000    4.0000000    7.0000000 ]  iseI:     true
        eV:      [   0.0000000    4.0000000    7.0000000   11.0000000 ]  iseV:     false
        eOP:     [  -1.0000000    0.0000000    4.0000000    7.0000000 ]  iseOP:    true
        eOPT:    [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseOPT:   false
        eOPTI:   [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseOPTI:  false
        eOPTT:   [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseOPTT:  false
        eOPTTI:  [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseOPTTI: false
        sum:     10

CM7.eOPT():
        Chord:   [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  name:     B-2 D-1 F#-1 G-1 
        pcs:     [   1.5000000    5.5000000    6.5000000   10.5000000 ]
        pitches: G-2 B-2 D-1 F#-1
        type:    [   0.0000000    1.0000000    5.0000000    8.0000000 ]
        eO:      [   6.5000000   -1.5000000    1.5000000    5.5000000 ]  iseO:     true
        eP:      [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseP:     true
        eT:      [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseT:     true
        eTT:     [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseTT:    true
        eI:      [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseI:     true
        eV:      [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseV:     true
        eOP:     [  -1.5000000    1.5000000    5.5000000    6.5000000 ]  iseOP:    true
        eOPT:    [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseOPT:   true
        eOPTI:   [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseOPTI:  true
        eOPTT:   [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseOPTT:  true
        eOPTTI:  [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseOPTTI: true
        sum:     0

CM7.eOPTT():
        Chord:   [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  name:     GM7 
        pcs:     [   2.0000000    6.0000000    7.0000000   11.0000000 ]
        pitches: G-2 B-2 D-1 F#-1
        type:    [   0.0000000    1.0000000    5.0000000    8.0000000 ]
        eO:      [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseO:     true
        eP:      [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseP:     true
        eT:      [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseT:     false
        eTT:     [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseTT:    true
        eI:      [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseI:     true
        eV:      [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseV:     true
        eOP:     [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseOP:    true
        eOPT:    [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseOPT:   false
        eOPTI:   [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseOPTI:  false
        eOPTT:   [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseOPTT:  true
        eOPTTI:  [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseOPTTI: true
        sum:     2

CM7.eOPTTI():
        Chord:   [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  name:     GM7 
        pcs:     [   2.0000000    6.0000000    7.0000000   11.0000000 ]
        pitches: G-2 B-2 D-1 F#-1
        type:    [   0.0000000    1.0000000    5.0000000    8.0000000 ]
        eO:      [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseO:     true
        eP:      [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseP:     true
        eT:      [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseT:     false
        eTT:     [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseTT:    true
        eI:      [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseI:     true
        eV:      [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseV:     true
        eOP:     [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseOP:    true
        eOPT:    [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseOPT:   false
        eOPTI:   [  -5.5000000   -1.5000000    1.5000000    5.5000000 ]  iseOPTI:  false
        eOPTT:   [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseOPTT:  true
        eOPTTI:  [  -5.0000000   -1.0000000    2.0000000    6.0000000 ]  iseOPTTI: true
        sum:     2

Cm7: [  60.0000000   63.0000000   67.0000000   70.0000000 ]
        Chord:   [  60.0000000   63.0000000   67.0000000   70.0000000 ]  name:     Eb6 
        pcs:     [   0.0000000    3.0000000    7.0000000   10.0000000 ]
        pitches: C4 D#4 G4 A#4
        type:    [   0.0000000    2.0000000    5.0000000    9.0000000 ]
        eO:      [   0.0000000    3.0000000    7.0000000   -2.0000000 ]  iseO:     false
        eP:      [  60.0000000   63.0000000   67.0000000   70.0000000 ]  iseP:     true
        eT:      [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseT:     false
        eTT:     [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseTT:    false
        eI:      [  60.0000000   63.0000000   67.0000000   70.0000000 ]  iseI:     true
        eV:      [  60.0000000   63.0000000   67.0000000   70.0000000 ]  iseV:     true
        eOP:     [  -2.0000000    0.0000000    3.0000000    7.0000000 ]  iseOP:    false
        eOPT:    [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPT:   false
        eOPTI:   [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTI:  false
        eOPTT:   [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTT:  false
        eOPTTI:  [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTTI: false
        sum:     260

Cm7.eOP():
        Chord:   [  -2.0000000    0.0000000    3.0000000    7.0000000 ]  name:     Eb6 
        pcs:     [   0.0000000    3.0000000    7.0000000   10.0000000 ]
        pitches: A#-2 C-1 D#-1 G-1
        type:    [   0.0000000    2.0000000    5.0000000    9.0000000 ]
        eO:      [  -2.0000000    0.0000000    3.0000000    7.0000000 ]  iseO:     true
        eP:      [  -2.0000000    0.0000000    3.0000000    7.0000000 ]  iseP:     true
        eT:      [  -4.0000000   -2.0000000    1.0000000    5.0000000 ]  iseT:     false
        eTT:     [  -4.0000000   -2.0000000    1.0000000    5.0000000 ]  iseTT:    false
        eI:      [  -2.0000000    0.0000000    3.0000000    7.0000000 ]  iseI:     true
        eV:      [   0.0000000    3.0000000    7.0000000   10.0000000 ]  iseV:     false
        eOP:     [  -2.0000000    0.0000000    3.0000000    7.0000000 ]  iseOP:    true
        eOPT:    [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPT:   false
        eOPTI:   [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTI:  false
        eOPTT:   [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTT:  false
        eOPTTI:  [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTTI: false
        sum:     8

Cm7.eOPT():
        Chord:   [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  name:     Bb6 
        pcs:     [   2.0000000    5.0000000    7.0000000   10.0000000 ]
        pitches: G-2 A#-2 D-1 F-1
        type:    [   0.0000000    2.0000000    5.0000000    9.0000000 ]
        eO:      [   7.0000000   -2.0000000    2.0000000    5.0000000 ]  iseO:     true
        eP:      [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseP:     true
        eT:      [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseT:     true
        eTT:     [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseTT:    true
        eI:      [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseI:     true
        eV:      [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseV:     true
        eOP:     [  -2.0000000    2.0000000    5.0000000    7.0000000 ]  iseOP:    true
        eOPT:    [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPT:   true
        eOPTI:   [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTI:  true
        eOPTT:   [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTT:  true
        eOPTTI:  [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTTI: true
        sum:     0

Cm7.eOPTT():
        Chord:   [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  name:     Bb6 
        pcs:     [   2.0000000    5.0000000    7.0000000   10.0000000 ]
        pitches: G-2 A#-2 D-1 F-1
        type:    [   0.0000000    2.0000000    5.0000000    9.0000000 ]
        eO:      [   7.0000000   -2.0000000    2.0000000    5.0000000 ]  iseO:     true
        eP:      [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseP:     true
        eT:      [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseT:     true
        eTT:     [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseTT:    true
        eI:      [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseI:     true
        eV:      [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseV:     true
        eOP:     [  -2.0000000    2.0000000    5.0000000    7.0000000 ]  iseOP:    true
        eOPT:    [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPT:   true
        eOPTI:   [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTI:  true
        eOPTT:   [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTT:  true
        eOPTTI:  [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTTI: true
        sum:     0

Cm7.eOPTTI():
        Chord:   [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  name:     Bb6 
        pcs:     [   2.0000000    5.0000000    7.0000000   10.0000000 ]
        pitches: G-2 A#-2 D-1 F-1
        type:    [   0.0000000    2.0000000    5.0000000    9.0000000 ]
        eO:      [   7.0000000   -2.0000000    2.0000000    5.0000000 ]  iseO:     true
        eP:      [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseP:     true
        eT:      [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseT:     true
        eTT:     [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseTT:    true
        eI:      [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseI:     true
        eV:      [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseV:     true
        eOP:     [  -2.0000000    2.0000000    5.0000000    7.0000000 ]  iseOP:    true
        eOPT:    [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPT:   true
        eOPTI:   [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTI:  true
        eOPTT:   [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTT:  true
        eOPTTI:  [  -5.0000000   -2.0000000    2.0000000    5.0000000 ]  iseOPTTI: true
        sum:     0

gogins commented 4 years ago

"Central inversion" see https://en.wikipedia.org/wiki/Point_reflection#Reflection_through_the_origin

gogins commented 4 years ago

I need (it is obvious by now) to be crystal clear at each and every step. Any blank or fuzzy step or guess (it is obvious by now) turns into an error with probability approaching 1. I believe the other symmetries are now clear enough to me. But not inversional equivalence (my iseI function, the I in OPTI).

Tymoczko's predicate for OPTI is reasonably simple and I have correctly coded it. The code for sending a chord to its equivalent in an orbifold is also simple for most of the orbifolds, but not for I and especially not for OPTI. The code I had written for doing this was simply wrong, even if I was able to use it to generate interesting music.

_It is clear that iseI must reflect in a hyperplane that is an inversion flat of OPT._ Everything I thought about inversion was wrong. OPT is a hyperplane in OP of codimension 1. The inversion flat is also a hyperplane; to be precise, it is in at least on case an affine hyperplane, consisting of those chords that are invariant under reflection. The "cyclical region" is tiled by n fundamental domains for OPT and by n x 2 fundamental domains for OPTI. There can thus be more than one inversion flat in some cyclical regions, and also the inversion flat must be defined in the cyclical region that is picked out by Tymoczko's "smallest wraparound interval" predicate.

The math for doing the reflection is now clear enough. If the inversion flat intersects the origin, let the reflection of a non-eI chord c be eI(c) = H c where H is the Householder reflector. H can be written either as I - 2 u x u where I is the identity matrix, u is a unit vector normal to the inversion flat, and x is outer product, or as I - 2 u u^T where u^T is the transpose of u, or as I - 2 [( u u^T)/(u^T u)]. If the inversion flat does not intersect the origin, H can still be used by moving the flat to the origin, moving the chord by the same translation, reflecting the chord with H, and moving the chord back by the inverse translation.

So far, so good. The remaining questions are:

[x] How to identify the correct fundamental domain in the cyclical region.
[x] How to define the correct inversion flat in that funamental domain.

But this doesn't actually seem so easy. The alternative solution is to use iseX to build maps of equivalents in OP, OPT, OPTT, and OPTTI for all chords in some equal temperament within R (ops_for_chords, optts_for_chords, opttis_for_chords). Then one simply looks up eOP, eOPTT, eOPTTI from the chord or, more simply, from eOP which I do know how to derive. The only way I can see to do that is to start with the "representative" iseX chords and then invert and/or translate them. But the inversion will need to be in the inversion flat, and if I have that, I don't need this approach.

But Tymoczko's code must be doing the business with the inversion flat, or some trick that simplifies that. Maybe it's the following. I'm about to cave and look at Tymoczko's Max code -- if I can. (I did, I could, and it doesn't do the eX transformations.)

The cyclic region for n-note chords is the (n-1)-simplicial region of R^n / T 
with n vertices at A_i = (0^{n − i}, 12^n)_T, 0 ≤ i < n. Replacing A_{n − 1} 
with a vertex at the maximally-even set yields a fundamental domain of OPT in R^n. 
Starting with this fundamental domain, replacing A_{n − 2} with a vertex at 
(0, 0, 6, ..., 6)_T yields a fundamental domain of OPTI in R_n . This construction is 
due to Noam Elkies.

Applying Elkies' construction to tetrachords we get 4 vertices (all these are transformed to T): (0, 0, 0, 0)_1, (0, 0, 0, 12)_2, (0, 0, 12, 12)_3, (0, 12, 12, 12)_4. Replace A_3 with the maximally even chord (0, 3, 6, 9). Replace A_2 with (0, 0, 6, 6). This gives vertices (0, 0, 0, 0)_1, (0, 0, 6, 6)_2, (0, 3, 6, 9)_3, (0, 12, 12, 12)_4. Does this yield the fundamental domain for OPTI that is used in the Science article? Not that I can see.

gogins commented 4 years ago

Maybe the flat is orthogonal to the line from the origin to the evenly spaced chord, that would simplify things. No, that's too simple.

I installed Max 8 on Heidi's Windows 7 computer and tried to run Chord Geometries. Then I installed Max on my Ubuntu computer with Wine. Max runs, and Tymoczko's patches load, and I can edit them, but they do not run.

To the extent I could read the Max patches, Tymoczko's logic seems much the same as mine, and does not seem to go as far as producing equivalents of chords (my eX functions), but only has the predicates (my iseX functions). But I need to look more closely, of course. There is a best_so_far object that might be used in a search for a point in the inversion flat, or for the eOPTTI chord.

I had another thought, I can simply put the inversion flats (in the form of normal vectors) identified in the draft paper for the Science article into a table, and that should enable eOPTTI for chords of up to 7 voices in 12TET (i.e., diatonic scales) (assuming the draft is correct!). Possibly doing this would lead to better insight and a more general solution, but this would get my ChordSpaceGroup going in 12TET for all immediate purposes, and perhaps I'd never need more.

[x] Examine Tymoczko's Max code further, especially to see what the best_so_far thing is used for. It turns out to be a variable that is lifted into the namespace of the patcher. I do not see anything resembling my eX code for transforming arbitrary chords into their equivalents in some orbifold. Indeed, Tymoczko's code seems to be in an alpha state.
[x] Try alternatives to the tensor form of the reflection. No need to do that, u x u is the same as u u^T.
[x] Try Noam Elkies' construction for OPTI fundamental domains, and compare results to those obtained using a table of normal vectors to the inversion flat from the Science draft. But what I really need is points to define inversion flats, not points to define fundamental domains.

And then I will call it a day for the JavaScript and move on to the C++.

gogins commented 4 years ago

There are some GitHub repositories turned up by https://github.com/search?q=tymoczko that I should look at. I will also search the following resources for others' implementations of something like eOPTI or eOPTTI using search terms like "tymoczko", "music geometry", "music orbifold", and "music symmetry".

[x] GitHub, mostly seems to be a dozen or so student notebooks.
[x] Assembla, can't browse or search repositories.
[x] Bitbucket, no chord space matches.
[x] Gitlab, nothing obvious.
[x] SourceForge, no matches.
[x] https://common-lisp.net/phub, see Gitlab.
[x] https://hackage.haskell.org/ There's https://hackage.haskell.org/package/Kulitta. The names on this code include both Donya Quick and Paul Hudak. The last edit was in 2016. Although the file is named Kulitta.ChordSpaces.OPTIC the code does not in fact define OPTI, only OP and OPT.

gogins commented 4 years ago

Some lessons...

This project is the hardest thing I have ever done, and appears to be at the very limit of what I am capable of. Although I did indeed have the idea of chords as points in chord space well before Tymoczko (I wrote it down at a woof group meeting and have it in a notebook from I think 1989), I never made anything of this idea, and I am certain that I never would have been capable of making anything of the idea, because I just do not have the requisite mathematical training and, possibly, ability.
When working in this mode of just trying to implement recipes, the slightest degree of unclarity is fatal and will produce incorrect code.
Oh yes, as always the code will be incorrect without thorough unit tests.
At times, the authors of the recipes will themselves have made typos, errors, hasty assumptions, or failed to consider the general case. When that happens, it is imperative to talk to the authors. I would not have been able to implement the other orbifolds (OP, OPT) without going back and forth quite a few times with Dmitri Tymoczko to clear up not only my, but also his unclarities.
The difficulties associated with the inversion flats in OPTI for different cardinalities probably explain why this approach was not fully implemented in other projects that use chord space, such as the Haskell School of Music.
In spite of all that, I can actually do this.
I think the musical payoff might be huge, because my distinct impression is that a big part of the emotional charge in Western music comes from the song reflecting through those damned inversion flats. But we will see....

gogins commented 4 years ago

[x] Find the n - 2 points that define the inversion flats that bound OPTI, either from the Science draft or starting from Noam Elkies' construction. I give up, I'm using the Science draft.
[x] Find vectors normal to the inversion flats from their vertices. I can't use the cross product which works only in 3 dimensions. but I can set up and solve a system of linear equations the solution to which is the plane equation, and if the plane equation is c_1 x_1 + c_2 x_2 ... + c_n x_n = a then (c_1, c_2, ... cn) is a normal vector. Note that a is the distance from an affine hyperplane to the origin. These two steps can be combined into one step, see https://math.stackexchange.com/questions/2723294/how-to-determine-the-equation-of-the-hyperplane-that-contains-several-points, https://en.wikipedia.org/wiki/Kernel(linear_algebra), https://math.stackexchange.com/questions/211392/how-do-you-compute-the-normal-vector-to-a-hyperplane-in-mathbbrn-given-n, and https://madoshakalaka.github.io/2019/03/02/generalized-cross-product-for-high-dimensions.html.

You make an (n − 1) by n matrix with row i given by x_i − x_n. Then the normal vector is the null space of this matrix. The fast way to do this is Gauss elimination.

Let V . n = z, and solve for a normal vector n. V is a matrix formed from n linearly independent points of a hyperplane; here, n non-collinear vertices of the fundamental domain. There is one row in the matrix for each vertex, with a final element of 1. z is zeros, i.e. the origin. Then n = numeric.solve(V, z) is a normal vector.

Actually, only the singular value decomposition from vectors derived from points works in all cases.

gogins commented 4 years ago

Does this define the inversion flat for n voices?

Create n - 2 chords in the inversion flat by adding tritones to the origin, e.g. (0, ..., 6), (0, ..., 6, 6).
Remove equivalent chords, e.g. (0, 0, 0, 6, 6) and (0, 0, 6, 6, 6) are inverses of each other in OPT.
Add the maximally even chord.

Doesn't work for 7 voices....

gogins commented 4 years ago

I have not been able to find a general procedure for identifying inversion flats for chords with any number of voices. Come to think of it, neither were Callender, Quinn, or Tymoczko. I will however check the mathematics arXiv.

I checked the home pages of Callender, Quinn, and Tymoczko and searched the mathematics arXiv. Interesting stuff, but nothing to help me out here.

Assuming I don't find something published, I will just put in the flats as a table. I now have the math to do the rest.

gogins commented 4 years ago

[x] The derivation of a vector normal to the inversion flat from points in it is still not 100% clear to me, although I am quite certain I can clear it up. I will take care of this first. Using https://en.wikipedia.org/wiki/Kernel_(linear_algebra), create a matrix A with one homogeneous column vector of n + 1 elements for each of n - 1 linearly independent points in the inversion flat. Then A . n = 0, where 0 is the origin and n is a homogeneous normal vector of the inversion flat. Matrix A must be square. No, use the singular value decomposition with vectors derived from points.
[x] I would prefer to do the reflection directly taking into account that the inversion flat is an affine hyperplane.
[x] The inversion flats in the Science material are given in pitch-class coordinates, but they need to be translated to pitch coordinates. Or do they? If so, this can presumably be done by normalizing the chords in OP or OPT (but which? the material is not clear on this). What does the subscript T on these chords mean, anyway? "...we can use the subscript 'T' to indicate equivalence-classes under transposition." So this is just like representing OPTI pitches as set-class pitch-set classes, first voice always 0. Note, the inversion flat is in OPT. Just going to use eT.
[x] If these flats do not reflect to form the fundamental domain that has the least wraparound interval, I am at a loss. They don't, but I will try finding points in a list of OPT chords.
[ ] Once this is working, I will go over the documentation comments again, thoroughly.

These are the chords that define the inversion flats.

gogins commented 4 years ago

Either the inversion flats listed in the Science material are not defined by linearly independent points, or I do not understand how to implement the math. I am backing up and trying again.

[x] Code the examples given online in Python until I get the correct answers. I get the correct answers using vectors in the hyperplane with both the "generalized cross product" and the matrix null space. A system of homogeneous linear equations would work, but seems to take one more defining point, so I did not implement that.
[ ] Code the Science flats in Python. The flats are given as n - 1 points, but n points are needed to obtain the n - 1 linearly independent vectors for the calculations. For 4 voices, a 4th point can be found by taking the midpoint of the inversion line that evenly divides the inversion flat, which is one face of the OPTI fundamental domain. The question is now whether the inversion flat thus defined works for the correct sector of the cyclical region.
[x] Translate the Python to JavaScript.
[ ] Translate the Python and JavaScript to C++. Eigen 3 has a null space function.

gogins commented 4 years ago

Now trying with inversion flats from Science...

subtrahend: [-1.5, 0, 1.5, 0] vector[ 0 ]: [ 0. -1.5 -3. 4.5] vector[ 1 ]: [-1.5 -3. 1.5 3. ] vector[ 2 ]: [1.5 1.5 1.5 4.5] points: [[-1.5, -1.5, -1.5, 4.5], [-3, -3, 3, 3], [0, 1.5, 3, 4.5], [-1.5, 0, 1.5, 0]] vectors: [array([ 0. , -1.5, -3. , 4.5]), array([-1.5, -3. , 1.5, 3. ]), array([1.5, 1.5, 1.5, 4.5])] nullspace: [[-0.84266484] [ 0.51856298] [-0.06482037] [ 0.12964074]] basis_vector: [1, 0, 0, 0] Matrix: [[ 0. -1.5 -3. 4.5] [-1.5 -3. 1.5 3. ] [ 1.5 1.5 1.5 4.5] [ 1. 0. 0. 0. ]] basis_vector: [0, 1, 0, 0] Matrix: [[ 0. -1.5 -3. 4.5] [-1.5 -3. 1.5 3. ] [ 1.5 1.5 1.5 4.5] [ 0. 1. 0. 0. ]] basis_vector: [0, 0, 1, 0] Matrix: [[ 0. -1.5 -3. 4.5] [-1.5 -3. 1.5 3. ] [ 1.5 1.5 1.5 4.5] [ 0. 0. 1. 0. ]] basis_vector: [0, 0, 0, 1] Matrix: [[ 0. -1.5 -3. 4.5] [-1.5 -3. 1.5 3. ] [ 1.5 1.5 1.5 4.5] [ 0. 0. 0. 1. ]] generalized_cross_product: [87.75, -53.99999999999999, 6.750000000000005, -13.499999999999991] unit normal vector: [ 0.84266484 -0.51856298 0.06482037 -0.12964074] Array(4)0: (4) [0, -1.5, -3, 4.5]1: (4) [-1.5, -3, 1.5, 3]2: (4) [-4.5, 0, -1.5, 6]3: (4) [1, 0, 0, 0]length: 4proto: Array(0) ChordSpace.js:3111 determinant: 54 ChordSpace.js:3108 matrix: ChordSpace.js:3109 Array(4)0: (4) [0, -1.5, -3, 4.5]1: (4) [-1.5, -3, 1.5, 3]2: (4) [-4.5, 0, -1.5, 6]3: (4) [0, 1, 0, 0]length: 4proto: Array(0) ChordSpace.js:3111 determinant: 54 ChordSpace.js:3108 matrix: ChordSpace.js:3109 Array(4)0: (4) [0, -1.5, -3, 4.5]1: (4) [-1.5, -3, 1.5, 3]2: (4) [-4.5, 0, -1.5, 6]3: (4) [0, 0, 1, 0]length: 4proto: Array(0) ChordSpace.js:3111 determinant: 54 ChordSpace.js:3108 matrix: ChordSpace.js:3109 Array(4)0: (4) [0, -1.5, -3, 4.5]1: (4) [-1.5, -3, 1.5, 3]2: (4) [-4.5, 0, -1.5, 6]3: (4) [0, 0, 0, 1]length: 4proto: Array(0)

gogins commented 4 years ago

Now it's clear that just deriving a Householder reflector from the Science inversion flats is not enough, because they are affine hyperplanes, and therefore I need the distance of the flats from the origin.

gogins commented 4 years ago

From https://courses.csail.mit.edu/6.034s/handouts/spring12/recitation10-LinearAlgebra.pdf:

A hyperplane is a higher-dimensional generalization of lines and planes. The equation of a hyperplane is w · x + b = 0, where w is a vector normal to the hyperplane and b is an offset. Note that we can multiply by any constant and preserve the equality; if we multiply by 1/kwk, we get a new equation ˆw · x + b 0 = 0, where ˆw = w/kwk is the unit normal vector and b 0 = b/kwk is the distance from the hyperplane to the origin.

gogins commented 4 years ago

The points defining an inversion flat must be linearly independent, so a matrix formed of these points must have a determinant of 0. This is not the case here. If I can't get round this I am stopped.

I do now understand why this happening. The vectors do not increment the span. The graphics and chords in the Science material imply linearly independent points that are not listed.

gogins commented 4 years ago

https://math.stackexchange.com/questions/99299/best-fitting-plane-given-a-set-of-points

gogins commented 4 years ago

OK, I'm starting to dimly light up. There are many ways to find the equation for the inversion flat. I have implemented functions to do this with the generalized cross product, with the matrix nullspace both by vectors and by homogenized points, with singular value decomposition, and with least squares. I can't get the least squares method to work, but the others all work and all agree with each other, so I must have solved this. I will use the matrix nullspace or the singular value decomposition as they work with overdetermined systems. For ChordSpace.js I must use singular value decomposition.

gogins commented 4 years ago

Distance to origin: https://math.stackexchange.com/questions/1210545/distance-from-a-point-to-a-hyperplane

gogins commented 4 years ago

Back at work on this after a break for taking care of Donna's apartment.

In general, it's a problem that the Science material gives chords in pc-set form, but I need them in chord space coordinate form. It would definitely also be a problem if the pc-sets are rounded off to semitones from OT equivalence.

[x] Review and further test alternative code for deriving the reflection. For n points, the generalized cross product works. For more than n (overdetermined) points, nullspace by vectors, nullspace by points, and SVD by vectors work. The vector math for reflection in a hyperplane from the Wikipedia seems to work and is an involution. I think I can proceed on this basis.
[x] In particular, get the least squares procedure for fitting a hyperplane to work. Can't get it to work, forget about it unless my current approach proves unworkable.
[x] See if Elkies' procedure produces the flats identified in the Science material. Oops, not flats, fundamental domains!

gogins commented 4 years ago

I am trying again to define an inversion flat for 4 voices from the Science material. From Figure S5 it appears that [0,0,0,0], [0, 0, 0, 6], [0, 3, 6, 0], and [0, 1, 2, 9] should define this flat. It does seem to.

gogins commented 4 years ago

For 4 voices, my reflection code to divide OPT into OPTI seems to be working. At least, of the chords I tested, those that I know are on the inversion flat remain invariant, and the others are involutions. I think this proves that I do have the right equation for the inversion flat. No, there were two vanishing singular values.

[x] Complete code for 4, 5, 6, and 7 voices based on the Science inversion flats.
[ ] Complete tests to ensure that reflections are actually in OPT (i.e., are in the correct segment of the cyclical region, the one with the least "wraparound" interval). Chords in OPT must upon reflection always remain in OPT. They should upon reflection either remain invariant and in OPTI (chords on inversion flat) or alternate between OPTI and ~OPTI (and of course between I and ~I). To perform these tests adequately, I may do them only in JavaScript and/or C++ after updating the code, or I must ensure that my Python code uses mpmath which has its own sufficient linear algebra routines. I am going to do this first in C++ and open a separate issue for that, then probably also do it in JavaScript.
[x] Only the SVD-based code for determining the inversion flat from chords seems to work with the Science points. Find out why. It seems the SVD is a deep method that can be used to compute eigenvalue decompositions, total least squares, and matrix kernels among other things. Algorithms for computing the SVD are fairly recent (1954, current method 1970). There is always an SVD for any matrix. The singular values and vectors should be in descending order. If the least singular value is near 0 then the method probably does identify the inversion flat. This can be tested by reflecting chords on and off the inversion flat.
[x] See if Elkies' construction produces the same inversion flats as the body of the Science material. I can only do this, of course, if I can find a procedure for determining the inversion flat from the fundamental domain. If so, use both Elkies' construction and the Science flats to determine the reflections. If they agree, implement the final code using Elkies' construction. I am not doing this now, but I may come back to it.
[ ] Translate the Python code to JavaScript, along with enough instrumentation to prove it does the same things as the Python.
[ ] Correct the ChordSpaceGroup JavaScript code to use the new reflection.
[x] Translate the Python code to C++, along with enough instrumentation to prove it does the same things as the Python. I am creating a new issue for this.
[x] Correct the ChordSpaceGroup C++ code to use the new reflection. I am creating a new issue for this.
[ ] Review and correct the comments and documentation.
[ ] Go back to the original motivation -- compare the C++ and JavaScript implementations of the HarmonyIfs code.

gogins commented 4 years ago

I have continued to work on these issues in #145, which I just closed, so I an closing this also. Any problems that crop up after this will be dealt with in new issues.

gogins / csound-extended

Fix bugs in ChordSpace code for porting HarmonyIFS and ChordSpace to CsoundAC in C++ and WebAssembly #142