feat: FreqResponsePointwise — batch sweep with guaranteed single-point arithmetic

README.md

@@ -138,6 +138,7 @@ func main() {
 | Method | Description |
 |--------|-------------|
 | `FreqResponse` | H(jw) at given frequencies |
+| `FreqResponsePointwise` | H(jw) with each frequency bit-identical to a single-point `FreqResponse` call, regardless of sweep length |
 | `Bode` | Magnitude (dB) and phase (deg) vs frequency |
 | `Nyquist` | Nyquist plot with encirclement counting |
 | `Nichols` | Nichols chart (magnitude vs phase) |

frequency.go

@@ -113,6 +113,19 @@ func (sys *System) FreqResponse(omega []float64) (*FreqResponseMatrix, error) {
 	return newFrequencyEvaluator(sys).response(omega)
 }
 
+// FreqResponsePointwise evaluates the frequency response with guaranteed
+// per-frequency single-point arithmetic: the value at each omega[k] is
+// bit-identical to FreqResponse([]float64{omega[k]}), regardless of
+// len(omega). FreqResponse may switch long sweeps of delay-free
+// state-space models to a transfer-function conversion whose values are
+// close but not bit-identical to the single-point path;
+// FreqResponsePointwise never does, at the cost of one dense solve per
+// frequency. Use it when downstream comparisons require sweep results to
+// reproduce single-point evaluations exactly.
+func (sys *System) FreqResponsePointwise(omega []float64) (*FreqResponseMatrix, error) {
+	return newFrequencyEvaluator(sys).responsePointwise(omega)
+}
+
 func (sys *System) Bode(omega []float64, nPoints int) (*BodeResult, error) {
 	if omega == nil {
 		var err2 error
@@ -192,6 +205,45 @@ func (e frequencyEvaluator) response(omega []float64) (*FreqResponseMatrix, erro
 	return e.matrix(data, omega), nil
 }
 
+// responsePointwise evaluates each frequency exactly as response would for
+// a one-element sweep: direct state-space solve first, per-point
+// transfer-function fallback on solve failure, with the delay phase applied
+// per point using the flag of whichever path produced the value.
+func (e frequencyEvaluator) responsePointwise(omega []float64) (*FreqResponseMatrix, error) {
+	if len(omega) == 0 {
+		return nil, nil
+	}
+	if e.sys.HasInternalDelay() || e.sys.IsDescriptor() {
+		// These paths already evaluate one frequency at a time with
+		// batch-size-independent arithmetic.
+		return e.response(omega)
+	}
+
+	pm := e.p * e.m
+	data := make([]complex128, len(omega)*pm)
+	ws := newSSEvalWorkspace(e.n, e.p, e.m)
+	var tf *TransferFunc
+	for k, w := range omega {
+		s := e.sAt(w)
+		dst := data[k*pm : (k+1)*pm]
+		if err := evalFrSSInto(ws, e.sys, s, e.n, e.p, e.m); err == nil {
+			copy(dst, ws.g[:pm])
+			applyIODelayAtS(e.sys, s, dst, e.p, e.m, true)
+			continue
+		}
+		if tf == nil {
+			res, err := e.sys.TransferFunction(nil)
+			if err != nil {
+				return nil, err
+			}
+			tf = res.TF
+		}
+		tf.evalInto(s, dst)
+		applyIODelayAtS(e.sys, s, dst, e.p, e.m, false)
+	}
+	return e.matrix(data, omega), nil
+}
+
 func (e frequencyEvaluator) eval(s complex128) ([][]complex128, error) {
 	pm := e.p * e.m
 

frequency_pointwise_test.go

@@ -0,0 +1,196 @@
+package controlsys
+
+import (
+	"fmt"
+	"math"
+	"testing"
+
+	"gonum.org/v1/gonum/mat"
+)
+
+// FreqResponsePointwise guarantees that the value at each omega[k] is
+// bit-identical to FreqResponse([]float64{omega[k]}), regardless of sweep
+// length. These tests pin that guarantee across the sweep-length regimes
+// where FreqResponse itself changes evaluation strategy (the direct
+// state-space work limit) and across the structurally distinct paths:
+// delay-free, I/O delay, discrete, MIMO, internal delay (LFT), and the
+// per-point transfer-function fallback when the state-space solve fails.
+
+func pointwiseTestGrid(wMin, wMax float64, n int) []float64 {
+	return logspace(math.Log10(wMin), math.Log10(wMax), n)
+}
+
+func assertBitIdenticalToSinglePoint(t *testing.T, sys *System, omega []float64) {
+	t.Helper()
+
+	got, err := sys.FreqResponsePointwise(omega)
+	if err != nil {
+		t.Fatalf("FreqResponsePointwise: %v", err)
+	}
+	if got.NFreq != len(omega) {
+		t.Fatalf("NFreq = %d, want %d", got.NFreq, len(omega))
+	}
+
+	for k, w := range omega {
+		want, err := sys.FreqResponse([]float64{w})
+		if err != nil {
+			t.Fatalf("FreqResponse([%v]): %v", w, err)
+		}
+		for i := range got.P {
+			for j := range got.M {
+				g := got.At(k, i, j)
+				e := want.At(0, i, j)
+				if math.Float64bits(real(g)) != math.Float64bits(real(e)) ||
+					math.Float64bits(imag(g)) != math.Float64bits(imag(e)) {
+					t.Errorf("omega[%d]=%v out=%d in=%d: pointwise %v, single-point %v", k, w, i, j, g, e)
+				}
+			}
+		}
+	}
+}
+
+func TestFreqResponsePointwise_BitIdenticalAcrossSweepLengths(t *testing.T) {
+	firstOrder, err := New(
+		mat.NewDense(1, 1, []float64{-0.1}),
+		mat.NewDense(1, 1, []float64{0.2}),
+		mat.NewDense(1, 1, []float64{1}),
+		mat.NewDense(1, 1, []float64{0}),
+		0,
+	)
+	if err != nil {
+		t.Fatal(err)
+	}
+
+	secondOrderDelay, err := New(
+		mat.NewDense(2, 2, []float64{-0.125, 0, 0.125, -1.0 / 3.0}),
+		mat.NewDense(2, 1, []float64{0.1875, 0}),
+		mat.NewDense(1, 2, []float64{0, 1}),
+		mat.NewDense(1, 1, []float64{0}),
+		0,
+	)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if err := secondOrderDelay.SetInputDelay([]float64{1.5}); err != nil {
+		t.Fatal(err)
+	}
+
+	gain, err := NewGain(mat.NewDense(1, 1, []float64{3.5}), 0)
+	if err != nil {
+		t.Fatal(err)
+	}
+
+	systems := []struct {
+		name string
+		sys  *System
+	}{
+		{"first_order", firstOrder},
+		{"second_order_input_delay", secondOrderDelay},
+		{"static_gain", gain},
+	}
+
+	// 1 and 40 stay within the direct state-space work limit for these
+	// orders; 137 and 500 force FreqResponse onto its transfer-function
+	// sweep for any n >= 1.
+	for _, tc := range systems {
+		for _, nPoints := range []int{1, 40, 137, 500} {
+			t.Run(fmt.Sprintf("%s/%d_points", tc.name, nPoints), func(t *testing.T) {
+				assertBitIdenticalToSinglePoint(t, tc.sys, pointwiseTestGrid(0.001, 10.0, nPoints))
+			})
+		}
+	}
+}
+
+func TestFreqResponsePointwise_Discrete(t *testing.T) {
+	sysc, err := New(
+		mat.NewDense(2, 2, []float64{-0.125, 0, 0.125, -1.0 / 3.0}),
+		mat.NewDense(2, 1, []float64{0.1875, 0}),
+		mat.NewDense(1, 2, []float64{0, 1}),
+		mat.NewDense(1, 1, []float64{0}),
+		0,
+	)
+	if err != nil {
+		t.Fatal(err)
+	}
+	sysd, err := sysc.Discretize(0.5)
+	if err != nil {
+		t.Fatal(err)
+	}
+	assertBitIdenticalToSinglePoint(t, sysd, pointwiseTestGrid(0.001, 2.0, 137))
+}
+
+func TestFreqResponsePointwise_MIMO(t *testing.T) {
+	sys, err := New(
+		mat.NewDense(3, 3, []float64{-0.5, 0.1, 0, 0, -0.25, 0.2, 0, 0, -1}),
+		mat.NewDense(3, 2, []float64{1, 0, 0, 1, 0.5, 0.5}),
+		mat.NewDense(2, 3, []float64{1, 0, 0, 0, 1, 1}),
+		mat.NewDense(2, 2, []float64{0, 0.1, 0, 0}),
+		0,
+	)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if err := sys.SetOutputDelay([]float64{0, 2.0}); err != nil {
+		t.Fatal(err)
+	}
+	assertBitIdenticalToSinglePoint(t, sys, pointwiseTestGrid(0.001, 10.0, 137))
+}
+
+func TestFreqResponsePointwise_InternalDelay(t *testing.T) {
+	sys, err := New(
+		mat.NewDense(2, 2, []float64{-0.5, 0.1, 0, -0.25}),
+		mat.NewDense(2, 1, []float64{1, 0.5}),
+		mat.NewDense(1, 2, []float64{1, 0}),
+		mat.NewDense(1, 1, []float64{0}),
+		0,
+	)
+	if err != nil {
+		t.Fatal(err)
+	}
+	B2 := mat.NewDense(2, 1, []float64{0.5, 0.3})
+	C2 := mat.NewDense(1, 2, []float64{0.2, 0.4})
+	zero := mat.NewDense(1, 1, []float64{0})
+	if err := sys.SetInternalDelay([]float64{3}, B2, C2, zero, zero, zero); err != nil {
+		t.Fatal(err)
+	}
+	assertBitIdenticalToSinglePoint(t, sys, pointwiseTestGrid(0.001, 10.0, 137))
+}
+
+// A pure integrator makes the state-space pencil singular at omega=0, so
+// FreqResponse([]float64{0}) falls back to transfer-function evaluation
+// for that point while every other point solves directly. The pointwise
+// sweep must reproduce that mixed per-point behavior bit-for-bit.
+func TestFreqResponsePointwise_PerPointFallback(t *testing.T) {
+	sys, err := New(
+		mat.NewDense(1, 1, []float64{0}),
+		mat.NewDense(1, 1, []float64{1}),
+		mat.NewDense(1, 1, []float64{1}),
+		mat.NewDense(1, 1, []float64{0}),
+		0,
+	)
+	if err != nil {
+		t.Fatal(err)
+	}
+	omega := append([]float64{0}, pointwiseTestGrid(0.001, 10.0, 40)...)
+	assertBitIdenticalToSinglePoint(t, sys, omega)
+}
+
+func TestFreqResponsePointwise_EmptyOmega(t *testing.T) {
+	sys, err := New(
+		mat.NewDense(1, 1, []float64{-1}),
+		mat.NewDense(1, 1, []float64{1}),
+		mat.NewDense(1, 1, []float64{1}),
+		mat.NewDense(1, 1, []float64{0}),
+		0,
+	)
+	if err != nil {
+		t.Fatal(err)
+	}
+	resp, err := sys.FreqResponsePointwise(nil)
+	if err != nil {
+		t.Fatalf("FreqResponsePointwise(nil): %v", err)
+	}
+	if resp != nil {
+		t.Fatalf("FreqResponsePointwise(nil) = %v, want nil (matching FreqResponse)", resp)
+	}
+}