From e96881576e44ea4b194c735b82a053e4d7997fa4 Mon Sep 17 00:00:00 2001
From: James Joseph <171958596+jamestjat@users.noreply.github.com>
Date: Tue, 25 Aug 2026 14:21:08 -0700
Subject: [PATCH] perf: reuse delay matrices across FreqResponsePointwise omega
 points
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The pointwise sweep rebuilt the effective I/O delay matrix inside
applyIODelayAtS on every omega point (a []float64 plus a mat.Dense
header), making sweep allocations scale linearly with len(omega):
2 allocations per point on delayed systems. The matrices are pure
functions of the System, which is not mutated across the sweep, so
they are now computed once before the loop (the transfer-function
fallback variant lazily, alongside the fallback TF itself) and passed
to the unchanged applyIODelayMatrixAtS per-point arithmetic.

Bit-identity holds: the same solves run in the same order and every
per-point float operation is unchanged — only the allocation source
of the delay matrix moved outside the loop. The existing
TestFreqResponsePointwise_* pins stay byte-exact, and a new
regression test asserts allocations per sweep no longer grow with
sweep length (5 per sweep at 10 and at 1000 points, down from
23 and 2003).

Co-Authored-By: Claude Mythos 5 <noreply@anthropic.com>
---
 frequency.go                | 15 +++++++++++++--
 frequency_pointwise_test.go | 37 +++++++++++++++++++++++++++++++++++++
 2 files changed, 50 insertions(+), 2 deletions(-)

diff --git a/frequency.go b/frequency.go
index 8858ef9..814bf3d 100644
--- a/frequency.go
+++ b/frequency.go
@@ -222,13 +222,21 @@ func (e frequencyEvaluator) responsePointwise(omega []float64) (*FreqResponseMat
 	pm := e.p * e.m
 	data := make([]complex128, len(omega)*pm)
 	ws := newSSEvalWorkspace(e.n, e.p, e.m)
+	// The effective delay matrices are pure functions of sys, which is
+	// not mutated across the sweep, so hoisting them out of the loop
+	// keeps every per-point value bit-identical while dropping the
+	// per-omega allocations applyIODelayAtS would repeat.
+	delaySS := effectiveIODelayMatrix(e.sys, e.p, e.m, true)
+	var delayTF *mat.Dense
 	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)
+			if delaySS != nil {
+				applyIODelayMatrixAtS(e.sys, s, dst, e.p, e.m, delaySS)
+			}
 			continue
 		}
 		if tf == nil {
@@ -237,9 +245,12 @@ func (e frequencyEvaluator) responsePointwise(omega []float64) (*FreqResponseMat
 				return nil, err
 			}
 			tf = res.TF
+			delayTF = effectiveIODelayMatrix(e.sys, e.p, e.m, false)
 		}
 		tf.evalInto(s, dst)
-		applyIODelayAtS(e.sys, s, dst, e.p, e.m, false)
+		if delayTF != nil {
+			applyIODelayMatrixAtS(e.sys, s, dst, e.p, e.m, delayTF)
+		}
 	}
 	return e.matrix(data, omega), nil
 }
diff --git a/frequency_pointwise_test.go b/frequency_pointwise_test.go
index 4156e13..ab18fa0 100644
--- a/frequency_pointwise_test.go
+++ b/frequency_pointwise_test.go
@@ -175,6 +175,43 @@ func TestFreqResponsePointwise_PerPointFallback(t *testing.T) {
 	assertBitIdenticalToSinglePoint(t, sys, omega)
 }
 
+// The pointwise sweep reuses one workspace (solve buffers and effective
+// delay matrices) across omega points, so allocations per sweep must stay
+// constant as the sweep grows: only the result buffer scales with
+// len(omega), and that is a single allocation. A per-point allocation
+// regression (like rebuilding the delay matrix inside the loop) makes the
+// long-sweep count scale with n and fails this test.
+func TestFreqResponsePointwise_AllocsAmortizeAcrossSweepLengths(t *testing.T) {
+	sys, 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 := sys.SetInputDelay([]float64{1.5}); err != nil {
+		t.Fatal(err)
+	}
+
+	allocsAt := func(n int) float64 {
+		omega := pointwiseTestGrid(0.001, 10.0, n)
+		return testing.AllocsPerRun(20, func() {
+			if _, err := sys.FreqResponsePointwise(omega); err != nil {
+				t.Fatal(err)
+			}
+		})
+	}
+
+	short := allocsAt(10)
+	long := allocsAt(1000)
+	if long > short {
+		t.Fatalf("allocations scale with sweep length: %v allocs at 10 points, %v at 1000 (want constant)", short, long)
+	}
+}
+
 func TestFreqResponsePointwise_EmptyOmega(t *testing.T) {
 	sys, err := New(
 		mat.NewDense(1, 1, []float64{-1}),