From 2057b63baa40dd9ce4a2432aa9132d6a8077a33a Mon Sep 17 00:00:00 2001 From: sora-ryu Date: Thu, 26 Mar 2026 15:40:43 -0600 Subject: [PATCH 01/10] add copyright holder --- LICENSE | 2 ++ 1 file changed, 2 insertions(+) diff --git a/LICENSE b/LICENSE index 261eeb9..46a3717 100644 --- a/LICENSE +++ b/LICENSE @@ -2,6 +2,8 @@ Version 2.0, January 2004 http://www.apache.org/licenses/ + Copyright (c) 2026 Alliance for Energy Innovation, LLC + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION 1. Definitions. From 732992984fc070a476bc848e3bd63f1c3757b415 Mon Sep 17 00:00:00 2001 From: sora-ryu Date: Thu, 26 Mar 2026 15:43:49 -0600 Subject: [PATCH 02/10] remove placeholder --- LICENSE | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/LICENSE b/LICENSE index 46a3717..7d8166e 100644 --- a/LICENSE +++ b/LICENSE @@ -3,7 +3,7 @@ http://www.apache.org/licenses/ Copyright (c) 2026 Alliance for Energy Innovation, LLC - + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION 1. Definitions. @@ -188,8 +188,6 @@ same "printed page" as the copyright notice for easier identification within third-party archives. - Copyright [yyyy] [name of copyright owner] - Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at From 98de40f1d89ed11236be520597d37c703bbff0fc Mon Sep 17 00:00:00 2001 From: sora-ryu Date: Thu, 26 Mar 2026 15:51:01 -0600 Subject: [PATCH 03/10] add software record number --- README.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/README.md b/README.md index f8f09d8..11c3dc4 100644 --- a/README.md +++ b/README.md @@ -20,5 +20,8 @@ Executables can be downloaded for Windows, MacOS, and Linux from the Releases pa Please open issues on Github if you encounter problems with the software. If possible, provide a minimal case and instructions to reproduce the failure. +## Acknowledgements + +Released under software record NLR/SWR-26-042 From 49d03ac6ebfe70ec9149757f02769a70422f5a7e Mon Sep 17 00:00:00 2001 From: Sora Ryu <64963755+sora-ryu@users.noreply.github.com> Date: Thu, 2 Apr 2026 15:47:47 -0600 Subject: [PATCH 04/10] Add NOTICE file --- NOTICE | 13 +++++++++++++ 1 file changed, 13 insertions(+) create mode 100644 NOTICE diff --git a/NOTICE b/NOTICE new file mode 100644 index 0000000..fa39c04 --- /dev/null +++ b/NOTICE @@ -0,0 +1,13 @@ +Copyright (c) 2026 Alliance for Energy Innovation, LLC + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at +       +      http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. From f5047fa02ff32d93da86cd62576ec69d16246b0d Mon Sep 17 00:00:00 2001 From: Sora Ryu <64963755+sora-ryu@users.noreply.github.com> Date: Thu, 2 Apr 2026 15:49:48 -0600 Subject: [PATCH 05/10] Add copyright notice separately --- LICENSE | 2 -- 1 file changed, 2 deletions(-) diff --git a/LICENSE b/LICENSE index 7d8166e..cbd5c75 100644 --- a/LICENSE +++ b/LICENSE @@ -2,8 +2,6 @@ Version 2.0, January 2004 http://www.apache.org/licenses/ - Copyright (c) 2026 Alliance for Energy Innovation, LLC - TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION 1. Definitions. From e3b50bcafeb26ae2ebb5a0a1eb32d1393121b564 Mon Sep 17 00:00:00 2001 From: Sora Ryu <64963755+sora-ryu@users.noreply.github.com> Date: Thu, 2 Apr 2026 15:55:43 -0600 Subject: [PATCH 06/10] Keep verbatim Apache license --- LICENSE | 2 ++ 1 file changed, 2 insertions(+) diff --git a/LICENSE b/LICENSE index cbd5c75..261eeb9 100644 --- a/LICENSE +++ b/LICENSE @@ -186,6 +186,8 @@ same "printed page" as the copyright notice for easier identification within third-party archives. + Copyright [yyyy] [name of copyright owner] + Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at From 11f7b66aedacecb307c7e9e3c4ece903e299c09a Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Wed, 26 Aug 2026 11:40:03 -0600 Subject: [PATCH 07/10] lin: fix MACX denominator to use the first mode in both terms The second term of the first denominator factor summed md2.EigenVector[i] * md2.EigenVector[i] instead of the md1 product, so the factor evaluated to (phi1^H phi1 + |phi2^T phi2|) rather than (phi1^H phi1 + |phi1^T phi1|). That made the criterion asymmetric in its arguments and allowed it to exceed 1, which matters because connectModesMAC and spectralClustering compare modes in whichever order they happen to iterate. Checked against the formula from the referenced ISMA 2010 paper over 5000 random complex mode shapes. The previous expression disagreed with itself under argument swap on every pair, with a largest discrepancy of 1.70, and fell outside [0, 1] on 217 of them. With the fix the criterion matches the reference, is symmetric, stays within [0, 1], and returns exactly 1 for a mode compared against itself. Assisted-by: Kiro:claude-opus-5 --- lin/mode.go | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/lin/mode.go b/lin/mode.go index d8b8be0..2af6403 100644 --- a/lin/mode.go +++ b/lin/mode.go @@ -169,7 +169,11 @@ func (md1 Mode) MACX(md2 *Mode) (float64, error) { numer2 += md1.EigenVector[i] * md2.EigenVector[i] denom11 += md1.EigenVector[i] * cmplx.Conj(md1.EigenVector[i]) - denom12 += md2.EigenVector[i] * md2.EigenVector[i] + // NOTE: this term belongs to md1. Using md2 here makes the first + // denominator factor (|phi1^H phi1| + |phi2^T phi2|) instead of + // (|phi1^H phi1| + |phi1^T phi1|), so the criterion is no longer + // symmetric in its arguments and does not match the reference. + denom12 += md1.EigenVector[i] * md1.EigenVector[i] denom21 += md2.EigenVector[i] * cmplx.Conj(md2.EigenVector[i]) denom22 += md2.EigenVector[i] * md2.EigenVector[i] From 6ce28a9e40e5220e5caad81772f31d13464891fc Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Wed, 26 Aug 2026 11:42:03 -0600 Subject: [PATCH 08/10] diagram: fix row reduction and augmenting path buffer in assignment Three defects in the minimum cost assignment implementation. Step1 held a pointer into the row it was about to modify rather than a copy of the row minimum. Once the subtraction loop reached the position of the minimum, that element became 0, so every later element in the row had 0 subtracted instead of the row minimum. That is not a uniform row offset, so it changes which assignment is optimal rather than only shifting the dual variables. Against an exhaustive search over 3000 random square matrices the old code returned a suboptimal assignment 422 times. Step5 wrote into a path buffer of N entries, but the alternating series of primed and starred zeros reaches 2N-1 entries (N primes and N-1 stars). This panicked with an index out of range on a 5x5 cost matrix reached through MinCostAssignment, and on 170 of 6000 random matrices. The reference implementation allocates 2N. MinCostAssignment indexed cost[0] without checking for an empty matrix, so a caller with nothing to assign panicked instead of getting an empty result. Also replaced ':=' with '=' on Step4's findZero call. The declaration shadowed the outer row and col, which made the later 'col = star_col' dead and restarted every search from (0, 0). This one is not a correctness fix: over 3000 random matrices the shadowed version still returned an optimal assignment every time, so it removes redundant rescanning and brings the step back in line with the reference. Verified by comparing against exhaustive search over 4000 random cost matrices up to 6x6, including rectangular shapes. Assisted-by: Kiro:claude-opus-5 --- diagram/mcassign.go | 37 ++++++++++++++++++++++++++++--------- 1 file changed, 28 insertions(+), 9 deletions(-) diff --git a/diagram/mcassign.go b/diagram/mcassign.go index 4dfef97..c1f68c8 100644 --- a/diagram/mcassign.go +++ b/diagram/mcassign.go @@ -88,6 +88,11 @@ func NewIntMatrix(m, n, v int) IntMatrix { // through the matrix func MinCostAssignment(cost IntMatrix) (results [][2]int, err error) { + // Nothing to assign; return an empty result instead of indexing cost[0] + if len(cost) == 0 || len(cost[0]) == 0 { + return nil, nil + } + // Pad cost matrix so it is square, get size costSq := padMatrix(cost, 0) N := len(costSq) @@ -102,7 +107,10 @@ func MinCostAssignment(cost IntMatrix) (results [][2]int, err error) { Z0_r: 0, Z0_c: 0, Marked: NewIntMatrix(N, N, 0), - path: make([][2]int, N), + // Step5 builds an alternating series of primed and starred zeros that + // can reach 2N-1 entries (N primes and N-1 stars), so N slots is not + // enough and overflows for larger augmenting paths. + path: make([][2]int, 2*N), } done := false @@ -160,23 +168,30 @@ func (m *MCA) Step1() (int, error) { // Loop through rows in C for i := range m.C { - // Find minimum value in row, ignore invalid values - var minVal *int - for j, v := range m.C[i] { - if (v != INVALID) && (minVal == nil || v < *minVal) { - minVal = &m.C[i][j] + // Find minimum value in row, ignore invalid values. + // NOTE: this must be a copy of the minimum, not a pointer into the + // row. The subtraction loop below writes to the same row, so a + // pointer would be read back as 0 once the loop passes the position + // of the minimum, and every element after it would have 0 subtracted + // instead of the row minimum. That is not a uniform row offset, so it + // changes which assignment is optimal rather than only shifting the + // dual variables. + minVal := INVALID + for _, v := range m.C[i] { + if v != INVALID && v < minVal { + minVal = v } } // If no min value found, return error - if minVal == nil { + if minVal == INVALID { return 0, fmt.Errorf("all values in row %d are INVALID", i+1) } // Subtract minimum value from all values in row for j, v := range m.C[i] { if v != INVALID { - m.C[i][j] -= *minVal + m.C[i][j] -= minVal } } } @@ -265,7 +280,11 @@ func (m *MCA) Step4() (int, error) { col := 0 for { - row, col := m.findZero(row, col) + // NOTE: assign with '=' rather than ':='. Declaring new variables here + // shadows the outer row/col, which makes the 'col = star_col' + // assignment below dead and restarts every search from (0, 0) instead + // of resuming from the starred column. + row, col = m.findZero(row, col) if row < 0 { return 6, nil } From 2f8e1ecd024e540ba0745dfc0fda60a1d01b97f6 Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Wed, 26 Aug 2026 11:44:27 -0600 Subject: [PATCH 09/10] diagram: make mode set labels deterministic and handle empty operating points Unpaired modes were tracked in a map keyed by column index, and the loop that turned leftovers into new mode sets ranged over that map. Go randomizes map iteration order, so identical input produced different ID and Label assignments from run to run. Confirmed on a two operating point case with two unpaired modes, where labels 2 and 3 swapped between consecutive runs. The candidate modes are now held in a slice and walked in index order, with a parallel bool slice marking the ones that paired. An operating point whose modes were all filtered out left modeSets empty, and the following mat.NewDense call panicked with "zero length in matrix dimension". Candidate modes are now gathered up front, an operating point with no candidates is skipped, and if no mode sets exist yet each candidate seeds one. Sizing the weight matrix from len(op.Modes) rather than the number of candidates also left trailing all-zero columns that mat.Max could report as the maximum, and dividing by a zero maximum produced int(NaN) in the cost matrix, which the language leaves undefined. Both are now avoided. These two are hardening; I did not find an input where they change the resulting mode sets. Assisted-by: Kiro:claude-opus-5 --- diagram/connect.go | 94 +++++++++++++++++++++++++++++----------------- 1 file changed, 60 insertions(+), 34 deletions(-) diff --git a/diagram/connect.go b/diagram/connect.go index 521c2a1..30f1673 100644 --- a/diagram/connect.go +++ b/diagram/connect.go @@ -42,29 +42,49 @@ func connectModesMAC(OPs []lin.LinOP, freqRangeHz [2]float64, structMax bool) ([ continue } - // Create empty weighting matrix - w := mat.NewDense(len(modeSets), len(op.Modes), nil) + // Collect the modes in this operating point that pass the filter. + // NOTE: these are gathered up front so the weight matrix has exactly + // one column per candidate. Sizing it from len(op.Modes) leaves + // trailing all-zero columns that are never written but are still + // visible to mat.Max below. + filteredModes := []*lin.Mode{} + for l := range op.Modes { + mn := &op.Modes[l] + if mn.Filter(freqRangeHz, structMax) { + filteredModes = append(filteredModes, mn) + } + } + + // No candidate modes in this operating point, nothing to connect + if len(filteredModes) == 0 { + continue + } + + // No mode sets to connect to yet, because no mode in any earlier + // operating point passed the filter. Seed one set per candidate mode; + // building a zero-row weight matrix below would panic. + if len(modeSets) == 0 { + for _, mn := range filteredModes { + modeSets = append(modeSets, &ModeSet{ + ID: len(modeSets), + Label: fmt.Sprintf("%d", len(modeSets)), + Modes: []*lin.Mode{mn}, + }) + } + continue + } - // Create map mapping mode index to mode - modeIndexMap := map[int]*lin.Mode{} + // Create empty weighting matrix + w := mat.NewDense(len(modeSets), len(filteredModes), nil) - // Loop through modes in mode set map + // Loop through mode sets for j, modeSet := range modeSets { // Get last mode in mode set mp := modeSet.Modes[len(modeSet.Modes)-1] - // Loop through modes in current operating point - k := 0 - for l := range op.Modes { - - // Get mode - mn := &op.Modes[l] - - // If mode should not be filtered, continue - if !mn.Filter(freqRangeHz, structMax) { - continue - } + // Loop through candidate modes in current operating point + for k, mn := range filteredModes { // Calculate MAC between modes mac, err := mp.MAC(mn) @@ -77,23 +97,22 @@ func connectModesMAC(OPs []lin.LinOP, freqRangeHz [2]float64, structMax bool) ([ // Add MAC to weight matrix w.Set(j, k, mac) - - // Add mode to index map - modeIndexMap[k] = mn - - k++ } } // Get max weight value wMax := mat.Max(w) - // Create cost matrix (ints) from weights (rescale to maximize precision) - cost := NewIntMatrix(len(modeSets), len(modeIndexMap), 0) - for j := range cost { - for k := range cost[j] { - v := w.At(j, k) - cost[j][k] = int(1e7 * (1 - v/wMax)) + // Create cost matrix (ints) from weights (rescale to maximize + // precision). If nothing correlates at all then every cost is equal, + // and the guard is required because dividing by a zero wMax yields + // NaN, whose conversion to int is not defined by the language spec. + cost := NewIntMatrix(len(modeSets), len(filteredModes), 0) + if wMax > 0 { + for j := range cost { + for k := range cost[j] { + cost[j][k] = int(1e7 * (1 - w.At(j, k)/wMax)) + } } } @@ -103,25 +122,32 @@ func connectModesMAC(OPs []lin.LinOP, freqRangeHz [2]float64, structMax bool) ([ return nil, err } - // Add connected modes to sets + // Add connected modes to sets, tracking which candidates were paired + paired := make([]bool, len(filteredModes)) for _, pair := range pairs { // Look up mode set from previous mode index modeSet := modeSets[pair[0]] // Add paired mode to slice of modes - modeSet.Modes = append(modeSet.Modes, modeIndexMap[pair[1]]) + modeSet.Modes = append(modeSet.Modes, filteredModes[pair[1]]) - // Remove paired mode from map - delete(modeIndexMap, pair[1]) + // Mark paired candidate mode + paired[pair[1]] = true } - // Loop through unpaired modes and create new mode sets - for _, m := range modeIndexMap { + // Loop through unpaired modes and create new mode sets. + // NOTE: walk the slice in index order rather than ranging over a map. + // Go randomizes map iteration order, so the IDs and labels assigned to + // these new mode sets varied between runs on identical input. + for k, mn := range filteredModes { + if paired[k] { + continue + } modeSets = append(modeSets, &ModeSet{ ID: len(modeSets), Label: fmt.Sprintf("%d", len(modeSets)), - Modes: []*lin.Mode{m}, + Modes: []*lin.Mode{mn}, }) } } From 70a808544215be1047ac0f7c52c0d258a7161bfe Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Wed, 26 Aug 2026 11:46:00 -0600 Subject: [PATCH 10/10] diagram: guard against zero degree modes in spectral clustering A mode with no MAC correlation to any other mode in its group has a zero degree, so 1/sqrt(di) stored +Inf in D^-1/2 and NaN propagated through Lsym and the eigen solve. The damage was silent rather than a crash: on a group of two mode sets containing one such mode, both sets came back with zero modes instead of their assigned modes. The scaling is now left at zero for those rows. Rows of the feature matrix are only normalized when their norm is non-zero. Scaling by 1/0 fills the observation with NaN, and because every comparison against NaN is false, the point lands in whichever cluster is checked first rather than the nearest one. Also corrected the comment above the eigenvalue sort, which said largest to smallest while the comparator sorts smallest first. The smallest eigenvalues of the symmetric Laplacian carry the cluster structure, so those are the ones that belong in the feature matrix. No behavior change. Assisted-by: Kiro:claude-opus-5 --- diagram/cluster.go | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/diagram/cluster.go b/diagram/cluster.go index a4724df..4b83b66 100644 --- a/diagram/cluster.go +++ b/diagram/cluster.go @@ -91,7 +91,12 @@ func spectralClustering(modeSets []*ModeSet) error { } di := mat.Sum(W.RowView(i)) D.Set(i, i, di) - D_isr.Set(i, i, 1/math.Sqrt(di)) + // A mode with no MAC correlation to any other mode in the group has a + // zero degree. Leave its scaling at zero instead of storing +Inf, + // which would propagate NaN through Lsym and the eigen solve. + if di > 0 { + D_isr.Set(i, i, 1/math.Sqrt(di)) + } } // Calculate Laplacian matrix (D - W) @@ -112,7 +117,10 @@ func spectralClustering(modeSets []*ModeSet) error { eigenVectors := &mat.CDense{} eig.VectorsTo(eigenVectors) - // Get indices that would sort from largest to smallest eigenvalues + // Get indices that would sort from smallest to largest eigenvalues. + // The comparator below is '<', and the smallest eigenvalues of the + // symmetric Laplacian are the ones that carry the cluster structure, so + // those are what get selected for the feature matrix. indices := argsort.SortSlice(eigenValues, func(i, j int) bool { return real(eigenValues[i]) < real(eigenValues[j]) }) @@ -126,7 +134,13 @@ func spectralClustering(modeSets []*ModeSet) error { for j, ind := range indices[:numDims] { row[j] = real(eigenVectors.At(i, ind)) } - floats.Scale(1/floats.Norm(row, 2), row) + // Only normalize a row with a non-zero norm. Scaling by 1/0 fills the + // observation with NaN, and every distance comparison against NaN is + // false, so the point silently lands in whichever cluster is checked + // first instead of the nearest one. + if norm := floats.Norm(row, 2); norm > 0 { + floats.Scale(1/norm, row) + } d[i] = Observation(row) }