From ca6182329c69798fe202962f700ff88200047950 Mon Sep 17 00:00:00 2001 From: James Joseph <171958596+jamestjat@users.noreply.github.com> Date: Mon, 24 Aug 2026 23:48:36 -0700 Subject: [PATCH] =?UTF-8?q?feat:=20FreqResponsePointwise=20=E2=80=94=20swe?= =?UTF-8?q?ep=20with=20guaranteed=20single-point=20arithmetic?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit FreqResponse switches long sweeps of delay-free state-space models to a transfer-function conversion (directFrequencySweepWorkLimit) whose values are close but not bit-identical to the single-point path. Downstream consumers that need sweep results to reproduce single-point evaluations exactly had to mirror the internal work limit and chunk their sweeps. FreqResponsePointwise contracts per-frequency single-point arithmetic: the value at each omega[k] is bit-identical to FreqResponse([]float64{omega[k]}) regardless of sweep length, including the per-point transfer-function fallback when the state-space solve fails at an isolated frequency. Internal-delay (LFT) and descriptor paths already evaluate one frequency at a time and delegate unchanged. Co-Authored-By: Claude Mythos 5 --- README.md | 1 + frequency.go | 52 ++++++++++ frequency_pointwise_test.go | 196 ++++++++++++++++++++++++++++++++++++ 3 files changed, 249 insertions(+) create mode 100644 frequency_pointwise_test.go diff --git a/README.md b/README.md index 0dc294a..dffc6ad 100644 --- a/README.md +++ b/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) | diff --git a/frequency.go b/frequency.go index d20383e..8858ef9 100644 --- a/frequency.go +++ b/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 diff --git a/frequency_pointwise_test.go b/frequency_pointwise_test.go new file mode 100644 index 0000000..4156e13 --- /dev/null +++ b/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) + } +}