Some bug fixes and hardening in the Campbell diagrams

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.

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
 
 

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)
 	}
 

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

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
 		}

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]