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. 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 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) } 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}, }) } } 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 } 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]