diff --git a/.ergo/plans.jsonl b/.ergo/plans.jsonl new file mode 100644 index 0000000..6e8a719 --- /dev/null +++ b/.ergo/plans.jsonl @@ -0,0 +1,90 @@ +{"type":"new_task","ts":"2026-07-18T20:09:01.023735Z","data":{"id":"FDNFWB","uuid":"b29e262f-3f32-4561-9dbc-bd7ee8e2b00d","epic_id":"","state":"todo","title":"Architectural and mathematical correctness","body":"","created_at":"2026-07-18T20:09:01.023735Z"}} +{"type":"new_task","ts":"2026-07-18T20:09:01.02797Z","data":{"id":"2H6QTC","uuid":"02f5b1b3-2fed-42d4-9793-79a169303da0","epic_id":"FDNFWB","state":"todo","title":"Correct ModelArray stacking semantics","body":"## Goal\n- Make stacking create a leading axis while preserving every source array shape and model.\n- Keep concatenation behavior explicit instead of hiding it behind stack terminology.\n\n## Acceptance Criteria\n- StackModelArrays requires equal non-empty shapes and returns [count]+shape.\n- Model ordering and metadata copying are preserved.\n- Shape mismatch, nil arrays, and empty input have deterministic errors.\n- README examples and behavioral tests describe true stacking.\n\n## Validation Gates\n- go test -run 'ModelArray|Stack'\n- go vet ./...","created_at":"2026-07-18T20:09:01.02797Z"}} +{"type":"new_task","ts":"2026-07-18T20:09:01.027977Z","data":{"id":"XIMOJQ","uuid":"d8ca4cc0-a62e-4010-9046-93c452307a8c","epic_id":"FDNFWB","state":"todo","title":"Make passivity results mathematically honest","body":"## Goal\n- Distinguish finite-grid passivity evidence from analytic certification.\n- Preserve compatibility while preventing sampled checks from presenting themselves as proofs.\n\n## Acceptance Criteria\n- Results disclose whether they are sampled or certified and retain the checked frequency range/grid.\n- State-space and FRD checks use correct Hermitian-part eigenvalue criteria for MIMO systems.\n- Tests cover narrow-band or out-of-range counterexamples and invalid grids.\n- Documentation states the exact guarantee.\n\n## Validation Gates\n- go test -run 'Passiv|Spectral'\n- go vet ./...","created_at":"2026-07-18T20:09:01.027977Z"}} +{"type":"new_task","ts":"2026-07-18T20:09:01.02798Z","data":{"id":"KFFJUD","uuid":"9c568d6f-cbb7-425b-9b86-f79988d40d7d","epic_id":"FDNFWB","state":"todo","title":"Implement basis-invariant modal truncation","body":"## Goal\n- Replace state-index deletion with a real ordered modal/Schur reduction workflow.\n- Preserve complex-conjugate blocks and descriptor-system correctness or reject unsupported cases explicitly.\n\n## Acceptance Criteria\n- Mode selection controls retained modal blocks rather than source coordinates.\n- Similarity-equivalent realizations retain the same selected poles and equivalent reduced responses.\n- Complex conjugate pairs are never split; unstable-mode policy is explicit.\n- Tests use coupled non-symmetric systems and compare pole/response invariants.\n\n## Validation Gates\n- go test -run 'Modal|Modred'\n- go test -race -run 'Modal|Modred'\n- go vet ./...","created_at":"2026-07-18T20:09:01.02798Z"}} +{"type":"new_task","ts":"2026-07-18T20:09:01.027982Z","data":{"id":"V7TYI7","uuid":"a43a226c-82ba-4bb1-bbde-18543108fc56","epic_id":"FDNFWB","state":"todo","title":"Assemble physical port connections into dynamics","body":"## Goal\n- Make physical connections and grounding alter the assembled equations rather than only validate metadata.\n- Give assembly one owner for topology, constraints, elimination, and external-port construction.\n\n## Acceptance Criteria\n- Connected and disconnected component assemblies have demonstrably different transfer behavior.\n- Compatible through/across port equations and grounding are represented mathematically.\n- Invalid, duplicate, underdetermined, and overconstrained topologies return actionable errors.\n- Descriptor constraints are preserved when algebraic elimination is not valid.\n- Tests cover coupling, grounding, multi-component topology, and non-symmetric dynamics.\n\n## Validation Gates\n- go test -run 'Physical|Assembly|Port'\n- go test -race -run 'Physical|Assembly|Port'\n- go vet ./...","created_at":"2026-07-18T20:09:01.027982Z"}} +{"type":"new_task","ts":"2026-07-18T20:09:01.027985Z","data":{"id":"JA5SHN","uuid":"0417cefc-9a27-4c1d-8b9f-63a952377f8c","epic_id":"FDNFWB","state":"todo","title":"Deepen generalized loops and tuning goals","body":"## Goal\n- Bind analysis points to actual interconnection channels and extract the requested loop.\n- Make MIMO sensitivity and tuning-goal evaluation mathematically correct.\n- Give grid tuning an honest, explicit policy while preserving useful extension seams.\n\n## Acceptance Criteria\n- Different analysis points can identify different loops in one topology.\n- Sensitivity reuses dimension-aware loop-sensitivity behavior.\n- MIMO gain goals use maximum singular value and honor frequency bands/weights.\n- Systune and Looptune either have distinct documented policies or delegate to an explicitly named shared grid tuner.\n- Tests cover multiple points, MIMO loops, custom TunableBlock implementations, and deterministic tuning limits.\n\n## Validation Gates\n- go test -run 'Generalized|AnalysisPoint|Tuning|Systune|Looptune|Loopsens'\n- go test -race -run 'Generalized|Tuning'\n- go vet ./...","created_at":"2026-07-18T20:09:01.027985Z"}} +{"type":"new_task","ts":"2026-07-18T20:09:01.027987Z","data":{"id":"CMLZFF","uuid":"d1dad137-0c41-411b-8c6f-82c3233ac95c","epic_id":"FDNFWB","state":"todo","title":"Integrate architectural correctness changes","body":"## Goal\n- Align package documentation and examples with the corrected mathematical contracts.\n- Run repository-wide static, race, numerical, and benchmark checks before publication.\n\n## Acceptance Criteria\n- README and package docs no longer overstate tracer capabilities.\n- All exported API docs match implemented guarantees and compatibility notes.\n- Full tests, race tests, go fix, go vet, and relevant benchmarks pass.\n- Commits remain grouped by feature and the branch is ready for a draft PR.\n\n## Validation Gates\n- go fix ./...\n- go vet ./...\n- go test -v -count=1 ./...\n- go test -race ./...\n- go test -bench=. -benchmem ./...","created_at":"2026-07-18T20:09:01.027987Z"}} +{"type":"link","ts":"2026-07-18T20:09:11.615311Z","data":{"from_id":"CMLZFF","to_id":"2H6QTC","type":"depends"}} +{"type":"link","ts":"2026-07-18T20:09:11.712916Z","data":{"from_id":"CMLZFF","to_id":"XIMOJQ","type":"depends"}} +{"type":"link","ts":"2026-07-18T20:09:11.862506Z","data":{"from_id":"CMLZFF","to_id":"KFFJUD","type":"depends"}} +{"type":"link","ts":"2026-07-18T20:09:11.955566Z","data":{"from_id":"CMLZFF","to_id":"V7TYI7","type":"depends"}} +{"type":"link","ts":"2026-07-18T20:09:12.049981Z","data":{"from_id":"CMLZFF","to_id":"JA5SHN","type":"depends"}} +{"type":"claim","ts":"2026-07-18T20:09:21.963819Z","data":{"id":"2H6QTC","agent_id":"gpt-5@codex","ts":"2026-07-18T20:09:21.963819Z"}} +{"type":"state","ts":"2026-07-18T20:09:21.963819Z","data":{"id":"2H6QTC","state":"doing","ts":"2026-07-18T20:09:21.963819Z"}} +{"type":"result","ts":"2026-07-18T20:11:27.026516Z","data":{"task_id":"2H6QTC","summary":"model_array.go","path":"model_array.go","sha256_at_attach":"e305638757980cb032cfa577899f61559b360137c69a46db0c92acf575e33796","mtime_at_attach":"2026-07-18T20:10:13.854733865Z","git_commit_at_attach":"ed73d670c345ec12d9e845cfaceb2128c111cdef","ts":"2026-07-18T20:11:27.026516Z"}} +{"type":"body","ts":"2026-07-18T20:11:27.026516Z","data":{"id":"2H6QTC","body":"## Goal\n- Make stacking create a leading axis while preserving every source array shape and model.\n- Keep concatenation behavior explicit instead of hiding it behind stack terminology.\n\n## Acceptance Criteria\n- StackModelArrays requires equal non-empty shapes and returns [count]+shape.\n- Model ordering and metadata copying are preserved.\n- Shape mismatch, nil arrays, and empty input have deterministic errors.\n- README examples and behavioral tests describe true stacking.\n\n## Validation Gates\n- go test -run 'ModelArray|Stack'\n- go vet ./...\n\n## Completion\n- StackModelArrays now creates a leading axis and rejects unequal shapes.\n- ConcatModelArrays preserves the prior flatten-and-append operation under an honest name.\n- Added shape, ordering, void-entry, incompatibility, and concatenation coverage.\n- Validated with gopls diagnostics, focused tests, the full suite, and go vet.\n","ts":"2026-07-18T20:11:27.026516Z"}} +{"type":"state","ts":"2026-07-18T20:11:27.026516Z","data":{"id":"2H6QTC","state":"done","ts":"2026-07-18T20:11:27.026516Z"}} +{"type":"claim","ts":"2026-07-18T20:11:27.119203Z","data":{"id":"XIMOJQ","agent_id":"gpt-5@codex","ts":"2026-07-18T20:11:27.119203Z"}} +{"type":"state","ts":"2026-07-18T20:11:27.119203Z","data":{"id":"XIMOJQ","state":"doing","ts":"2026-07-18T20:11:27.119203Z"}} +{"type":"result","ts":"2026-07-18T20:15:56.288913Z","data":{"task_id":"XIMOJQ","summary":"passivity.go","path":"passivity.go","sha256_at_attach":"fdbcd7cf267a25e8857ba5145f40483f809088cf224c03d818d4c7ae6fe07d68","mtime_at_attach":"2026-07-18T20:14:46.772346488Z","git_commit_at_attach":"b105aa0cdeecd6c3d5c76b6df28ebe3ae6435f32","ts":"2026-07-18T20:15:56.288913Z"}} +{"type":"body","ts":"2026-07-18T20:15:56.288913Z","data":{"id":"XIMOJQ","body":"## Goal\n- Distinguish finite-grid passivity evidence from analytic certification.\n- Preserve compatibility while preventing sampled checks from presenting themselves as proofs.\n\n## Acceptance Criteria\n- Results disclose whether they are sampled or certified and retain the checked frequency range/grid.\n- State-space and FRD checks use correct Hermitian-part eigenvalue criteria for MIMO systems.\n- Tests cover narrow-band or out-of-range counterexamples and invalid grids.\n- Documentation states the exact guarantee.\n\n## Validation Gates\n- go test -run 'Passiv|Spectral'\n- go vet ./...\n\n## Completion\n- Added explicit violated, sampled-pass, and certified result statuses and retained the exact grid and tolerance.\n- Added SampledPassive while preserving Passive as a compatibility alias.\n- Included DC, bounded discrete-time checks at Nyquist, and rejected malformed grids, tolerances, FRD shapes, and non-finite data.\n- Added a persistent restricted-grid counterexample and completed full, race, diagnostic, and vet validation.\n","ts":"2026-07-18T20:15:56.288913Z"}} +{"type":"state","ts":"2026-07-18T20:15:56.288913Z","data":{"id":"XIMOJQ","state":"done","ts":"2026-07-18T20:15:56.288913Z"}} +{"type":"claim","ts":"2026-07-18T20:15:56.38886Z","data":{"id":"KFFJUD","agent_id":"gpt-5@codex","ts":"2026-07-18T20:15:56.38886Z"}} +{"type":"state","ts":"2026-07-18T20:15:56.38886Z","data":{"id":"KFFJUD","state":"doing","ts":"2026-07-18T20:15:56.38886Z"}} +{"type":"result","ts":"2026-07-18T20:24:42.638055Z","data":{"task_id":"KFFJUD","summary":"modal_reduction.go","path":"modal_reduction.go","sha256_at_attach":"a8eb1f89dabe851fd369880a046b852679a1ec3299414e46770d521f22f35a5f","mtime_at_attach":"2026-07-18T20:20:16.446837181Z","git_commit_at_attach":"e029f71a069f98e965f0f9778881ac0e184097ad","ts":"2026-07-18T20:24:42.638055Z"}} +{"type":"body","ts":"2026-07-18T20:24:42.638055Z","data":{"id":"KFFJUD","body":"## Goal\n- Replace state-index deletion with a real ordered modal/Schur reduction workflow.\n- Preserve complex-conjugate blocks and descriptor-system correctness or reject unsupported cases explicitly.\n\n## Acceptance Criteria\n- Mode selection controls retained modal blocks rather than source coordinates.\n- Similarity-equivalent realizations retain the same selected poles and equivalent reduced responses.\n- Complex conjugate pairs are never split; unstable-mode policy is explicit.\n- Tests use coupled non-symmetric systems and compare pole/response invariants.\n\n## Validation Gates\n- go test -run 'Modal|Modred'\n- go test -race -run 'Modal|Modred'\n- go vet ./...\n\n## Completion\n- Replaced source-coordinate deletion with ordered real Schur decomposition and Sylvester block separation.\n- Preserved unstable modes and conjugate pairs and exposed retained poles plus left/right projection bases.\n- Added persistent similarity-invariance, projection, pair-preservation, unstable-mode, and threshold-selection tests.\n- Full tests, focused race tests, diagnostics, and vet pass.\n- Benchstat records 492.8 us/op and 146.8 KiB/op for the corrected N=50/order=10 algorithm versus 5.75 us/op and 25.5 KiB/op for the previous non-modal state slice.\n","ts":"2026-07-18T20:24:42.638055Z"}} +{"type":"state","ts":"2026-07-18T20:24:42.638055Z","data":{"id":"KFFJUD","state":"done","ts":"2026-07-18T20:24:42.638055Z"}} +{"type":"new_task","ts":"2026-07-18T20:24:42.731138Z","data":{"id":"4NPDOO","uuid":"8ed87628-27b4-4194-800d-8de87d884765","epic_id":"FDNFWB","state":"todo","title":"Update Gonum fork dependency","body":"","created_at":"2026-07-18T20:24:42.731138Z"}} +{"type":"claim","ts":"2026-07-18T20:24:42.731138Z","data":{"id":"4NPDOO","agent_id":"gpt-5@codex","ts":"2026-07-18T20:24:42.731138Z"}} +{"type":"state","ts":"2026-07-18T20:24:42.731138Z","data":{"id":"4NPDOO","state":"doing","ts":"2026-07-18T20:24:42.731138Z"}} +{"type":"body","ts":"2026-07-18T20:25:03.048658Z","data":{"id":"4NPDOO","body":"## Goal\n- Update the Gonum replacement to the latest live release on github.com/jamestjsp/gonum.\n- Verify that numerical and architectural changes remain correct against the updated fork.\n\n## Context\n- Live tag inspection on 2026-07-18 found v0.17.7-fork newer than the current v0.17.6-fork.\n\n## Acceptance Criteria\n- go.mod resolves gonum.org/v1/gonum through github.com/jamestjsp/gonum v0.17.7-fork.\n- Module files are tidy and no unrelated dependency changes occur.\n- Workspace vulnerability, static, full test, race, and modal numerical checks pass.\n\n## Validation Gates\n- go mod tidy\n- go vulncheck ./...\n- go fix ./...\n- go vet ./...\n- go test -v -count=1 ./...\n- go test -race ./...\n","ts":"2026-07-18T20:25:03.048658Z"}} +{"type":"link","ts":"2026-07-18T20:25:03.144556Z","data":{"from_id":"CMLZFF","to_id":"4NPDOO","type":"depends"}} +{"type":"result","ts":"2026-07-18T20:27:45.387877Z","data":{"task_id":"4NPDOO","summary":"go.mod","path":"go.mod","sha256_at_attach":"4757c07b5ab37cbaea17c69fd4dc06762d2ed0e6c9216b65e76f4de2dfee3811","mtime_at_attach":"2026-07-18T20:25:07.260545573Z","git_commit_at_attach":"e7552cf47698163a12eb138009db7dd71596b18e","ts":"2026-07-18T20:27:45.387877Z"}} +{"type":"body","ts":"2026-07-18T20:27:45.387877Z","data":{"id":"4NPDOO","body":"## Goal\n- Update the Gonum replacement to the latest live release on github.com/jamestjsp/gonum.\n- Verify that numerical and architectural changes remain correct against the updated fork.\n\n## Context\n- Live tag inspection on 2026-07-18 found v0.17.7-fork newer than the previous v0.17.6-fork.\n\n## Acceptance Criteria\n- go.mod resolves gonum.org/v1/gonum through github.com/jamestjsp/gonum v0.17.7-fork.\n- Module files are tidy and no unrelated dependency changes occur.\n- Workspace vulnerability, static, full test, race, and modal numerical checks pass.\n\n## Validation Gates\n- go mod tidy\n- go vulncheck ./...\n- go fix ./...\n- go vet ./...\n- go test -v -count=1 ./...\n- go test -race ./...\n\n## Completion\n- Updated the replace target and checksums to v0.17.7-fork, the latest live fork tag.\n- go mod tidy changed no unrelated modules.\n- Full tests, race tests, go fix, go vet, and diagnostics pass.\n- Vulnerability analysis found only two Go 1.26.4 standard-library advisories and no Gonum/dependency findings.\n- Canon/modal benchmark comparison showed no significant timing or allocation change.\n","ts":"2026-07-18T20:27:45.387877Z"}} +{"type":"state","ts":"2026-07-18T20:27:45.387877Z","data":{"id":"4NPDOO","state":"done","ts":"2026-07-18T20:27:45.387877Z"}} +{"type":"claim","ts":"2026-07-18T20:27:45.479569Z","data":{"id":"V7TYI7","agent_id":"gpt-5@codex","ts":"2026-07-18T20:27:45.479569Z"}} +{"type":"state","ts":"2026-07-18T20:27:45.479569Z","data":{"id":"V7TYI7","state":"doing","ts":"2026-07-18T20:27:45.479569Z"}} +{"type":"result","ts":"2026-07-18T20:40:29.598727Z","data":{"task_id":"V7TYI7","summary":"physical_assembly.go","path":"physical_assembly.go","sha256_at_attach":"0cc166244f3a0d904d0df1cf625ff075c91473983efc13aa639f20a71ad7a876","mtime_at_attach":"2026-07-18T20:38:52.841498285Z","git_commit_at_attach":"244a2068b0cc2fb04af3247770fb255504d6e64d","ts":"2026-07-18T20:40:29.598727Z"}} +{"type":"body","ts":"2026-07-18T20:40:29.598727Z","data":{"id":"V7TYI7","body":"## Goal\n- Make physical connections and grounding alter the assembled equations rather than only validate metadata.\n- Give assembly one owner for topology, constraints, elimination, and external-port construction.\n\n## Acceptance Criteria\n- Connected and disconnected component assemblies have demonstrably different transfer behavior.\n- Compatible through/across port equations and grounding are represented mathematically.\n- Invalid, duplicate, underdetermined, and overconstrained topologies return actionable errors.\n- Descriptor constraints are preserved when algebraic elimination is not valid.\n- Tests cover coupling, grounding, multi-component topology, and non-symmetric dynamics.\n\n## Validation Gates\n- go test -run 'Physical|Assembly|Port'\n- go test -race -run 'Physical|Assembly|Port'\n- go vet ./...\n\n## Completion\n- Added explicit physical-port input/output channel ownership with compatibility mapping for legacy positional ports.\n- Added an assembly-plan owner for component validation, node unioning, external-channel selection, descriptor constraints, and metadata.\n- Connected nodes enforce across-variable equality and through-variable conservation; grounded nodes preserve reaction variables.\n- Added descriptor frequency response through the regular pencil sE-A and direct scalar-oracle coverage.\n- Added pair, grounded, three-component, descriptor-preservation, and invalid-topology behavior tests.\n- Full and focused race tests, diagnostics, and vet pass.\n- Correct assembly benchmark is 24.9 us/op and 92.2 KiB/op versus 14.9 us/op and 55.1 KiB/op for the prior validation-only append.\n","ts":"2026-07-18T20:40:29.598727Z"}} +{"type":"state","ts":"2026-07-18T20:40:29.598727Z","data":{"id":"V7TYI7","state":"done","ts":"2026-07-18T20:40:29.598727Z"}} +{"type":"claim","ts":"2026-07-18T20:40:29.696252Z","data":{"id":"JA5SHN","agent_id":"gpt-5@codex","ts":"2026-07-18T20:40:29.696252Z"}} +{"type":"state","ts":"2026-07-18T20:40:29.696252Z","data":{"id":"JA5SHN","state":"doing","ts":"2026-07-18T20:40:29.696252Z"}} +{"type":"result","ts":"2026-07-18T20:57:32.859488Z","data":{"task_id":"JA5SHN","summary":"generalized.go","path":"generalized.go","sha256_at_attach":"04348604c2ce16487fcc6c5c920df195c2246f523037c689ead6fe0efd91157d","mtime_at_attach":"2026-07-18T20:56:28.108894409Z","git_commit_at_attach":"244a2068b0cc2fb04af3247770fb255504d6e64d","ts":"2026-07-18T20:57:32.859488Z"}} +{"type":"state","ts":"2026-07-18T20:57:32.859488Z","data":{"id":"JA5SHN","state":"done","ts":"2026-07-18T20:57:32.859488Z"}} +{"type":"claim","ts":"2026-07-18T20:57:52.859874Z","data":{"id":"CMLZFF","agent_id":"codex@gpt-5","ts":"2026-07-18T20:57:52.859874Z"}} +{"type":"state","ts":"2026-07-18T20:57:52.859874Z","data":{"id":"CMLZFF","state":"doing","ts":"2026-07-18T20:57:52.859874Z"}} +{"type":"result","ts":"2026-07-18T21:09:30.134851Z","data":{"task_id":"CMLZFF","summary":"docs/codebase-interface-diagram.md","path":"docs/codebase-interface-diagram.md","sha256_at_attach":"926f3544aa5bf37f9e532cb81a2595c257a2d5f0691118000c4cf9348c4fb00b","mtime_at_attach":"2026-07-18T20:59:46.593472404Z","git_commit_at_attach":"9efb609bbe5d07ccc3eab19a74d919ae6caffd4e","ts":"2026-07-18T21:09:30.134851Z"}} +{"type":"state","ts":"2026-07-18T21:09:30.134851Z","data":{"id":"CMLZFF","state":"done","ts":"2026-07-18T21:09:30.134851Z"}} +{"type":"new_task","ts":"2026-07-19T04:19:00.343832Z","data":{"id":"JU3PKB","uuid":"7d0dfdc6-4de5-431a-a104-9b2f1180a113","epic_id":"","state":"todo","title":"Deep numerical architecture and performance","body":"","created_at":"2026-07-19T04:19:00.343832Z"}} +{"type":"new_task","ts":"2026-07-19T04:19:00.343868Z","data":{"id":"2CSJE4","uuid":"285c5ee9-7a14-423a-932e-f529471df5ae","epic_id":"JU3PKB","state":"todo","title":"Replace inverse-based DARE with a generalized Schur solve","body":"\n## Goal\nMake the discrete algebraic Riccati solver robust when A is singular or ill-conditioned by owning the full generalized symplectic-pencil workflow inside the Riccati module.\n\n## Context\nThe current DARE path explicitly forms inverse-transpose(A), which rejects solvable singular-A problems and amplifies conditioning errors. SLICOT SB02OD uses an ordered generalized-Schur/deflating-subspace method, and the pinned Gonum fork provides the required Dgges/Dtgsen-class machinery.\n\n## Acceptance criteria\n- The implementation no longer computes an explicit inverse of A.\n- Persistent tests cover singular-A and ill-conditioned, non-symmetric DARE cases with fixed values produced by an independent oracle.\n- The returned X is symmetric, the scaled Riccati residual meets the existing tolerance policy, and the closed-loop eigenvalues are strictly inside the unit circle for stabilizing cases.\n- Existing CARE/DARE APIs, workspace reuse, error semantics, and allocation expectations remain compatible.\n- Focused benchmarks show no material regression on existing representative DARE sizes.\n\n## Validation gates\n- Run focused Riccati tests and oracle fixtures.\n- Run focused before/after DARE benchmarks through benchstat.\n- Run go fix ./..., go vet ./..., and go test -v -count=1.","created_at":"2026-07-19T04:19:00.343868Z"}} +{"type":"new_task","ts":"2026-07-19T04:19:00.343872Z","data":{"id":"O7HHVY","uuid":"83d3df33-7f0f-44d6-9669-3e22c4dcc60e","epic_id":"JU3PKB","state":"todo","title":"Deepen the frequency-sweep evaluator","body":"\n## Goal\nMake frequencyEvaluator own backend selection and reusable sweep workspaces so repeated frequencies do not rebuild transfer-function realizations unnecessarily.\n\n## Context\nProfiles show newTFWorkspace dominates allocation in FreqResponse and tuning. The evaluator currently converts a state-space model to transfer functions for multi-point sweeps, while single-point evaluation uses a reusable state-space solve. The LFT kernel also performs one duplicate diagonal multiply.\n\n## Acceptance criteria\n- frequencyEvaluator selects an evidence-backed state-space or transfer-function kernel based on model/grid shape, with the decision hidden from callers.\n- Sweep workspaces are allocated once per response operation and reused across frequencies.\n- The duplicate LFT diagonal multiply is removed.\n- Persistent parity tests cover both selected kernels, non-symmetric A, continuous/discrete models, descriptor models, delays, and LFT responses.\n- Numerical error is no worse than the current implementation against the existing reference behavior.\n- Focused FreqResponse and LFT benchmarks show lower allocations and no material time regression across representative small and larger systems.\n\n## Validation gates\n- Run focused frequency tests.\n- Run focused before/after frequency and LFT benchmarks through benchstat.\n- Run go fix ./..., go vet ./..., and go test -v -count=1.","created_at":"2026-07-19T04:19:00.343872Z"}} +{"type":"new_task","ts":"2026-07-19T04:19:00.343877Z","data":{"id":"7AJL7Y","uuid":"72cc9afe-9b31-4894-a654-7083f5d9d0e4","epic_id":"JU3PKB","state":"todo","title":"Unify complex singular-value evaluation","body":"\n## Goal\nGive the frequency-analysis numerical module one reusable singular-value implementation and make tuning goals consume flat response storage without per-frequency matrix cloning.\n\n## Context\ntuning_goal.go duplicates complex singular-value logic using mat.SVD and Values(nil), while frequency_svd.go already owns reusable buffers for Sigma and FRD operations. The duplicate path allocates heavily and can drift numerically.\n\n## Acceptance criteria\n- One internal complex-SVD workspace owns singular-value evaluation for Sigma, FRD, and tuning goals.\n- Tuning consumes stride-correct flat frequency-response slices without constructing a new matrix per frequency.\n- Persistent tests compare maximum singular values against Gonum SVD over rectangular, rank-deficient, scaled, and non-symmetric complex matrices.\n- Existing tuning-goal values and tolerances are preserved or improved.\n- Focused tuning benchmarks show lower allocations and no material time regression.\n\n## Validation gates\n- Run focused frequency-SVD and tuning tests.\n- Run focused before/after tuning benchmarks through benchstat.\n- Run go fix ./..., go vet ./..., and go test -v -count=1.","created_at":"2026-07-19T04:19:00.343877Z"}} +{"type":"new_task","ts":"2026-07-19T04:19:00.34388Z","data":{"id":"6VEROU","uuid":"4e6c134e-f057-4af1-98c7-8b909afe3805","epic_id":"JU3PKB","state":"todo","title":"Batch all-input step and impulse responses","body":"\n## Goal\nMake timeResponsePlan simulate all input channels in one state propagation pass for step and impulse responses.\n\n## Context\nStep and Impulse currently allocate an input history and call Simulate once per input. This repeats A-times-state work and intermediate allocations. The planner is the cohesive owner of discretization, delays, output layout, and simulation policy.\n\n## Acceptance criteria\n- Step and Impulse use a batched all-input kernel where applicable, with delay or exceptional paths safely routed through equivalent behavior.\n- The public row ordering j*p+i, time grids, feedthrough behavior, continuous/discrete semantics, and delay semantics remain unchanged.\n- The kernel uses stride-aware raw storage and preallocated buffers.\n- Persistent parity tests cover non-symmetric MIMO systems, continuous and discrete models, gain-only systems, feedthrough, and delays.\n- Focused MIMO benchmarks show fewer allocations and lower runtime; SISO has no material regression.\n\n## Validation gates\n- Run focused response tests.\n- Run focused before/after SISO and MIMO step/impulse benchmarks through benchstat.\n- Run go fix ./..., go vet ./..., and go test -v -count=1.","created_at":"2026-07-19T04:19:00.34388Z"}} +{"type":"new_task","ts":"2026-07-19T04:19:00.343882Z","data":{"id":"KD2R4T","uuid":"1e419cc6-bd35-4a2a-bf29-b2c11f5ef072","epic_id":"JU3PKB","state":"todo","title":"Make physical assembly own matrices once","body":"\n## Goal\nMove physical-model copy and validation ownership into physicalAssemblyPlan so assembly crosses one ownership boundary instead of cloning components and aggregate matrices repeatedly.\n\n## Context\nThe current path copies every component, copies into aggregate matrices, then the Descriptor constructor copies the aggregate again. The synchronous assembly plan can validate inputs while building an owned aggregate result directly.\n\n## Acceptance criteria\n- The plan produces an independently owned Descriptor model with no mutable aliasing to component matrices or port metadata.\n- Redundant component and final aggregate clones are removed through an internal validated ownership-transfer seam.\n- Constraint assembly uses stride-aware raw storage rather than repeated At/Set calls in hot loops.\n- Persistent mutation tests prove that changing source components after assembly cannot change the result, and changing the result cannot change sources.\n- Existing constraint, port, and validation behavior is unchanged.\n- Focused physical-assembly benchmarks show fewer allocations and no material time regression.\n\n## Validation gates\n- Run focused physical-assembly tests.\n- Run focused before/after assembly benchmarks through benchstat.\n- Run go fix ./..., go vet ./..., and go test -v -count=1.","created_at":"2026-07-19T04:19:00.343882Z"}} +{"type":"new_task","ts":"2026-07-19T04:19:00.343885Z","data":{"id":"7OSGT6","uuid":"0cb72a9c-7d18-415b-83ec-4da7646a9132","epic_id":"JU3PKB","state":"todo","title":"Integrate the numerical architecture and performance audit","body":"\n## Goal\nClose the audit with repository-wide correctness, static-analysis, and benchmark evidence, and record why other routines remain unchanged.\n\n## Context\nThe codebase-wide inventory found these five evidence-backed opportunities. Most remaining solver time is inside LAPACK/BLAS calls with existing reusable workspaces; changing those without a demonstrated bottleneck would increase risk without an architectural payoff.\n\n## Acceptance criteria\n- Every production routine family is classified as changed, covered by an existing deep owner, or intentionally unchanged with profile evidence.\n- All task-specific accuracy and benchmark gates pass.\n- Benchmark comparisons are recorded with benchstat for every changed hot path.\n- No go.mod change is made unless required; if it changes, go vulncheck passes.\n- The final worktree contains coherent commits for each completed task and one final commit for the completed ergo graph.\n\n## Validation gates\n- Run go fix ./....\n- Run go vet ./....\n- Run go test -v -count=1.\n- Run the representative full benchmark sweep and review benchstat summaries.","created_at":"2026-07-19T04:19:00.343885Z"}} +{"type":"link","ts":"2026-07-19T04:19:08.941268Z","data":{"from_id":"7AJL7Y","to_id":"O7HHVY","type":"depends"}} +{"type":"link","ts":"2026-07-19T04:19:09.053186Z","data":{"from_id":"7OSGT6","to_id":"2CSJE4","type":"depends"}} +{"type":"link","ts":"2026-07-19T04:19:09.16955Z","data":{"from_id":"7OSGT6","to_id":"7AJL7Y","type":"depends"}} +{"type":"link","ts":"2026-07-19T04:19:09.280233Z","data":{"from_id":"7OSGT6","to_id":"6VEROU","type":"depends"}} +{"type":"link","ts":"2026-07-19T04:19:09.458482Z","data":{"from_id":"7OSGT6","to_id":"KD2R4T","type":"depends"}} +{"type":"claim","ts":"2026-07-19T04:20:00.951578Z","data":{"id":"2CSJE4","agent_id":"gpt-5@codex","ts":"2026-07-19T04:20:00.951578Z"}} +{"type":"state","ts":"2026-07-19T04:20:00.951578Z","data":{"id":"2CSJE4","state":"doing","ts":"2026-07-19T04:20:00.951578Z"}} +{"type":"title","ts":"2026-07-19T04:37:21.810949Z","data":{"id":"2CSJE4","title":"Add generalized-Schur DARE for singular and ill-conditioned models","ts":"2026-07-19T04:37:21.810949Z"}} +{"type":"body","ts":"2026-07-19T04:37:21.810949Z","data":{"id":"2CSJE4","body":"## Goal\nMake the discrete algebraic Riccati solver robust when A is singular or ill-conditioned while preserving the fast path for well-conditioned models.\n\n## Context\nThe former DARE path rejected solvable singular-A problems. An always-generalized implementation fixed that but caused a material pure-Go QZ allocation and runtime regression. The Riccati module now owns a condition-based kernel choice: the existing symplectic Schur reduction for well-conditioned A and a SLICOT/SciPy-style extended symplectic pencil for singular or poorly conditioned A.\n\n## Acceptance criteria\n- Singular and poorly conditioned A use an ordered generalized-Schur deflating-subspace solve.\n- Well-conditioned models retain the established fast Schur kernel behind the same private owner.\n- Persistent tests cover singular-A cross-term and ill-conditioned, non-symmetric DARE cases with fixed values produced by an independent SciPy oracle.\n- The returned X is symmetric, the scaled Riccati residual meets the existing tolerance policy, and the closed-loop eigenvalues are strictly inside the unit circle for stabilizing cases.\n- Existing CARE/DARE APIs, workspace reuse, and error semantics remain compatible.\n- Focused benchmarks show no material regression on existing representative DARE sizes.\n\n## Validation gates\n- Run focused Riccati tests and oracle fixtures.\n- Run focused before/after DARE benchmarks through benchstat.\n- Run go fix ./..., go vet ./..., and go test -v -count=1.\n\n## Completion note\n- Added a private discrete-stable-subspace owner that chooses the regular or generalized kernel from the reciprocal condition estimate of A.\n- Implemented the extended symplectic pencil, QR input-column deflation, ordered generalized Schur solve, and generalized eigenvalue recovery.\n- Added singular cross-term and ill-conditioned non-symmetric SciPy oracle fixtures and corrected the cross-term residual helper.\n- Ten-sample benchstat: n=10 timing unchanged; n=50 +2.37%; allocation counts unchanged and bytes decreased slightly.\n","ts":"2026-07-19T04:37:21.810949Z"}} +{"type":"result","ts":"2026-07-19T04:38:21.148783Z","data":{"task_id":"2CSJE4","summary":"riccati.go","path":"riccati.go","sha256_at_attach":"013af0553718fdd1e4a1ecec7995e831fe42267a5260537de3a054dccbfeff3d","mtime_at_attach":"2026-07-19T04:35:36.157088804Z","git_commit_at_attach":"9c75eddbb50326459fad811dcd1fa7d90f03f86b","ts":"2026-07-19T04:38:21.148783Z"}} +{"type":"state","ts":"2026-07-19T04:38:21.148783Z","data":{"id":"2CSJE4","state":"done","ts":"2026-07-19T04:38:21.148783Z"}} +{"type":"claim","ts":"2026-07-19T04:38:26.256556Z","data":{"id":"O7HHVY","agent_id":"gpt-5@codex","ts":"2026-07-19T04:38:26.256556Z"}} +{"type":"state","ts":"2026-07-19T04:38:26.256556Z","data":{"id":"O7HHVY","state":"doing","ts":"2026-07-19T04:38:26.256556Z"}} +{"type":"body","ts":"2026-07-19T04:49:17.845033Z","data":{"id":"O7HHVY","body":"## Goal\nMake frequencyEvaluator own backend selection and reusable sweep workspaces so short sweeps avoid rebuilding transfer-function realizations while long sweeps keep polynomial-evaluation throughput.\n\n## Context\nProfiles showed newTFWorkspace dominated FreqResponse allocation, but always using direct state-space solves materially slowed longer grids. Persistent kernel benchmarks established a crossover governed primarily by state count times frequency count. The evaluator now owns that choice and falls back to transfer-function evaluation if a direct sweep encounters a pole.\n\n## Acceptance criteria\n- frequencyEvaluator selects the state-space or transfer-function kernel from an evidence-backed work limit hidden from callers.\n- State-space sweep workspaces are allocated once per response operation and reused across frequencies.\n- Persistent parity tests cover both selected kernels, a non-symmetric A, delay phases, descriptor evaluation, continuous/discrete behavior, and pole fallback.\n- Numerical error is no worse than EvalFr within the established tolerance.\n- Focused FreqResponse benchmarks show lower allocations and runtime for short sweeps with no material regression for longer, delayed, or LFT sweeps.\n\n## Validation gates\n- Run focused frequency tests.\n- Run focused before/after frequency and LFT benchmarks through benchstat.\n- Run go fix ./..., go vet ./..., and go test -v -count=1.\n\n## Completion note\n- Added private state-space and transfer-function sweep kernels under frequencyEvaluator and a conservative state-count by frequency-count work policy.\n- Added transfer-function fallback when a short direct sweep reaches a pole, preserving multi-frequency behavior.\n- Added persistent kernel-crossover and public short-sweep benchmarks plus parity coverage for both kernels and all delay forms.\n- Ten-sample benchstat for the 10-state, 8-frequency case: runtime -42.19%, bytes -85.74%, allocations -94.44%. Long, delayed, B747, and LFT sweeps were unchanged within noise.\n- The duplicate LFT multiply noted during initial planning was already absent after synchronizing with origin/main, so no LFT code change was needed.\n","ts":"2026-07-19T04:49:17.845033Z"}} +{"type":"result","ts":"2026-07-19T04:50:53.10985Z","data":{"task_id":"O7HHVY","summary":"frequency.go","path":"frequency.go","sha256_at_attach":"ee0bbd87cb24867a8101ac06cdd4177e2d6c100777871e9b64d7b2fca15e06a1","mtime_at_attach":"2026-07-19T04:48:48.329286428Z","git_commit_at_attach":"c717d4e322418e8a4d15997ba550c1a799edd26e","ts":"2026-07-19T04:50:53.10985Z"}} +{"type":"state","ts":"2026-07-19T04:50:53.10985Z","data":{"id":"O7HHVY","state":"done","ts":"2026-07-19T04:50:53.10985Z"}} +{"type":"claim","ts":"2026-07-19T04:51:03.117216Z","data":{"id":"7AJL7Y","agent_id":"gpt-5@codex","ts":"2026-07-19T04:51:03.117216Z"}} +{"type":"state","ts":"2026-07-19T04:51:03.117216Z","data":{"id":"7AJL7Y","state":"doing","ts":"2026-07-19T04:51:03.117216Z"}} +{"type":"body","ts":"2026-07-19T04:59:25.975012Z","data":{"id":"7AJL7Y","body":"# Unify complex singular-value evaluation\n\n## Goal\nGive the frequency-analysis numerical module one reusable singular-value implementation and make tuning goals consume flat response storage without per-frequency matrix cloning.\n\n## Context\ntuning_goal.go duplicated complex singular-value logic using mat.SVD and Values(nil), while frequency_svd.go already owned reusable buffers for Sigma and FRD operations. The duplicate path allocated heavily and could drift numerically.\n\n## Acceptance criteria\n- One internal complex-SVD workspace owns singular-value evaluation for Sigma, FRD, and tuning goals.\n- Tuning consumes stride-correct flat frequency-response slices without constructing a new matrix per frequency.\n- Persistent tests compare maximum singular values against Gonum SVD over rectangular, rank-deficient, scaled, and non-symmetric complex matrices.\n- Existing tuning-goal values and tolerances are preserved or improved.\n- Focused tuning benchmarks show lower allocations and no material time regression.\n\n## Completion note\n- `complexSVDWorkspace` now owns full and maximum singular-value evaluation for Sigma, FRD, and tuning goals.\n- Tuning reads frequency-major response blocks directly and allocates only weighted matrix products when weights are present.\n- Scale normalization prevents overflow and underflow in Gram formation; persistent Gonum-oracle tests cover rectangular, rank-deficient, non-symmetric, `1e-150`, `1e150`, and vector `1e200` inputs.\n- Ten-sample benchstat results: MIMO weighted gain time -52.11%, bytes -30.91%, allocations -80.91%; dynamic weighted MIMO time -53.25%, bytes -14.61%, allocations -66.07%; SISO time -7.93% with unchanged allocations.\n\n## Validation gates\n- Focused frequency-SVD and tuning tests passed.\n- Focused before/after tuning benchmarks compared through benchstat.\n- `go fix ./...`, `go vet ./...`, and `go test -v -count=1` passed.\n","ts":"2026-07-19T04:59:25.975012Z"}} +{"type":"result","ts":"2026-07-19T04:59:42.279704Z","data":{"task_id":"7AJL7Y","summary":"frequency_svd.go","path":"frequency_svd.go","sha256_at_attach":"76a01cd070e5839c9cebbd5ca7eee85b1ac3bde44db7af17ca17f7bc0656ebb3","mtime_at_attach":"2026-07-19T04:53:41.447520265Z","git_commit_at_attach":"208dbe3880230d7ee7d1ac50e3a72d116c5ebdca","ts":"2026-07-19T04:59:42.279704Z"}} +{"type":"state","ts":"2026-07-19T04:59:42.279704Z","data":{"id":"7AJL7Y","state":"done","ts":"2026-07-19T04:59:42.279704Z"}} +{"type":"claim","ts":"2026-07-19T04:59:57.858388Z","data":{"id":"6VEROU","agent_id":"gpt-5@codex","ts":"2026-07-19T04:59:57.858388Z"}} +{"type":"state","ts":"2026-07-19T04:59:57.858388Z","data":{"id":"6VEROU","state":"doing","ts":"2026-07-19T04:59:57.858388Z"}} +{"type":"body","ts":"2026-07-19T05:08:57.717542Z","data":{"id":"6VEROU","body":"# Batch all-input step and impulse responses\n\n## Goal\nMake `timeResponsePlan` simulate all input channels in one state propagation pass for step and impulse responses.\n\n## Context\nStep and Impulse allocated an input history and called Simulate once per input. That repeated dispatch, intermediate allocation, and output assembly. The planner now owns discretization, delays, output layout, and the all-input simulation policy.\n\n## Acceptance criteria\n- Step and Impulse use a batched all-input kernel where applicable, with delay or exceptional paths safely routed through equivalent behavior.\n- The public row ordering `j*p+i`, time grids, feedthrough behavior, continuous/discrete semantics, and delay semantics remain unchanged.\n- The kernel uses stride-aware raw storage and preallocated buffers.\n- Persistent parity tests cover non-symmetric MIMO systems, continuous and discrete models, gain-only systems, feedthrough, and delays.\n- Focused MIMO benchmarks show fewer allocations and lower runtime; SISO has no material regression.\n\n## Completion note\n- `timeResponsePlan` now owns standard-input response policy and batches delay-free, standard multi-input systems with contiguous input-major state buffers.\n- SISO, descriptor, zero-output, and delayed/internal-delay systems retain the established `Simulate` path.\n- Independent-simulation parity tests cover non-symmetric continuous and discrete MIMO models, feedthrough, gain-only models, delays, Step, Impulse, and output row ordering.\n- Ten-sample benchstat results: MIMO Step time -50.29%, bytes -70.51%, allocations -66.49%; MIMO Impulse time -50.27%, bytes -70.51%, allocations -66.49%. SISO Step and Impulse were about 2% faster with unchanged allocation counts.\n\n## Validation gates\n- Focused response tests passed.\n- Focused before/after SISO and MIMO Step/Impulse benchmarks compared through benchstat.\n- `go fix ./...`, `go vet ./...`, and `go test -v -count=1` passed.\n","ts":"2026-07-19T05:08:57.717542Z"}} +{"type":"result","ts":"2026-07-19T05:09:13.256925Z","data":{"task_id":"6VEROU","summary":"response.go","path":"response.go","sha256_at_attach":"b6b5962e76a9927cc81e1683d7c2e12746947d139706f0d326afb2c35877321a","mtime_at_attach":"2026-07-19T05:08:37.954259528Z","git_commit_at_attach":"c9da0c155d15ac0b04a4396d68ae3a8d3045f44a","ts":"2026-07-19T05:09:13.256925Z"}} +{"type":"state","ts":"2026-07-19T05:09:13.256925Z","data":{"id":"6VEROU","state":"done","ts":"2026-07-19T05:09:13.256925Z"}} +{"type":"claim","ts":"2026-07-19T05:09:13.262558Z","data":{"id":"KD2R4T","agent_id":"gpt-5@codex","ts":"2026-07-19T05:09:13.262558Z"}} +{"type":"state","ts":"2026-07-19T05:09:13.262558Z","data":{"id":"KD2R4T","state":"doing","ts":"2026-07-19T05:09:13.262558Z"}} +{"type":"body","ts":"2026-07-19T05:13:21.217114Z","data":{"id":"KD2R4T","body":"# Make physical assembly own matrices once\n\n## Goal\nMove physical-model copy and validation ownership into `physicalAssemblyPlan` so assembly crosses one ownership boundary instead of cloning components and aggregate matrices repeatedly.\n\n## Context\nThe connected path copied every component, copied into aggregate matrices, then copied the aggregate again through the public Descriptor constructor. The synchronous assembly plan can validate borrowed inputs while building an owned aggregate result directly.\n\n## Acceptance criteria\n- The plan produces an independently owned Descriptor model with no mutable aliasing to component matrices or port metadata.\n- Redundant component and final aggregate clones are removed through an internal validated ownership-transfer seam.\n- Constraint assembly uses stride-aware raw storage rather than repeated At/Set calls in hot loops.\n- Persistent mutation tests prove that changing source components after assembly cannot change the result, and changing the result cannot change sources.\n- Existing constraint, port, and validation behavior is unchanged.\n- Focused physical-assembly benchmarks show fewer allocations and no material time regression.\n\n## Completion note\n- The synchronous plan validates and borrows component inputs, then creates the result at the aggregate-matrix ownership boundary.\n- Internal `newDescriptorOwned` validates dimensions and descriptor shape while adopting plan-owned matrices; public `NewDescriptor` retains defensive copies.\n- Augmented matrices and constraint rows now use stride-aware raw storage.\n- Mutation tests cover matrices, signal metadata, and port metadata in both source-to-result and result-to-source directions.\n- Ten-sample benchstat results for eight connected components: time -30.96%, bytes -38.00%, allocations -29.00%.\n\n## Validation gates\n- Focused physical-assembly tests passed.\n- Focused before/after assembly benchmarks compared through benchstat.\n- `go fix ./...`, `go vet ./...`, and `go test -v -count=1` passed.\n","ts":"2026-07-19T05:13:21.217114Z"}} +{"type":"result","ts":"2026-07-19T05:13:39.412919Z","data":{"task_id":"KD2R4T","summary":"physical_assembly.go","path":"physical_assembly.go","sha256_at_attach":"22bf46ae9e42b2ae2fcfc66cf7daf44be2530b46587351a819c511ed69a0da0a","mtime_at_attach":"2026-07-19T05:13:13.667629616Z","git_commit_at_attach":"2d490c23f45b99b0cb5d77de1d063c37d7013ca5","ts":"2026-07-19T05:13:39.412919Z"}} +{"type":"state","ts":"2026-07-19T05:13:39.412919Z","data":{"id":"KD2R4T","state":"done","ts":"2026-07-19T05:13:39.412919Z"}} +{"type":"claim","ts":"2026-07-19T05:13:39.425093Z","data":{"id":"7OSGT6","agent_id":"gpt-5@codex","ts":"2026-07-19T05:13:39.425093Z"}} +{"type":"state","ts":"2026-07-19T05:13:39.425093Z","data":{"id":"7OSGT6","state":"doing","ts":"2026-07-19T05:13:39.425093Z"}} +{"type":"body","ts":"2026-07-19T05:23:22.131093Z","data":{"id":"7OSGT6","body":"# Integrate the numerical architecture and performance audit\n\n## Goal\nClose the audit with repository-wide correctness, static-analysis, and benchmark evidence, and record why other routines remain unchanged.\n\n## Context\nThe codebase-wide inventory found five evidence-backed opportunities. Most remaining solver time is inside LAPACK/BLAS calls with existing reusable workspaces; changing those without a demonstrated bottleneck would increase risk without an architectural payoff.\n\n## Production routine-family classification\n\n### Changed after measurement or oracle failure\n- Riccati and optimal-control kernel: `riccati.go` now routes singular or poorly conditioned DARE problems through an extended symplectic pencil and ordered generalized Schur decomposition while retaining the fast regular path.\n- Frequency sweep: `frequency.go` now owns the state-space versus transfer-function kernel decision and reusable short-sweep workspace.\n- Complex singular values and tuning: `frequency_svd.go` and `tuning_goal.go` share one scaled, reusable complex-SVD owner and flat frequency-response storage.\n- Standard time responses: `response.go` owns Step/Impulse policy and batches delay-free MIMO propagation while retaining established exceptional paths.\n- Physical assembly: `physical_assembly.go` and the internal seam in `state_space_utils.go` establish one validated matrix-ownership boundary.\n\n### Existing deep owners; intentionally unchanged\n- Core model, validation, domain, descriptor, metadata, and sampled layouts: `ss.go`, `descriptor.go`, `errors.go`, `names.go`, `time_domain.go`, `signal_metadata.go`, `sampled_data.go`, `sampled_response.go`, `sampled_signal.go`, `doc.go`.\n- Delay representation, conversion, topology, and safe feedback: `delay.go`, `delay_bank.go`, `delay_conversion_policy.go`, `delay_topology.go`, `delay_visibility.go`, `pade.go`, `thiran.go`, `feedback_lft.go`, `safe_feedback.go`, `wellposed.go`.\n- Interconnection and generalized-model workflows: `connect.go`, `interconnection_topology.go`, `lft.go`, `generalized.go`, `loopsens.go`, `sumblk.go`.\n- Conversion and realization workflows: `augstate.go`, `canon.go`, `conversion_plan.go`, `convert.go`, `inv.go`, `poly.go`, `polynomial_channel.go`, `pzmap.go`, `realization_policy.go`, `transfer.go`, `transform.go`, `zpk.go`.\n- Matrix equations, decompositions, balancing, reduction, and zeros: `balance.go`, `ctrbobsv.go`, `eigdecomp.go`, `gramian.go`, `lyapunov.go`, `matlog.go`, `matrix_equation_problem.go`, `minimal.go`, `modal_reduction.go`, `modsep.go`, `prescale.go`, `qr_colpivot.go`, `rq_rowpivot.go`, `sminreal.go`, `ssbal.go`, `stability.go`, `stabsep.go`, `staircase.go`, `zeros.go`. Profiles place the remaining cost in LAPACK/BLAS kernels; these routines already centralize workspace, policy, or result construction.\n- Frequency and robust analysis outside the changed sweep/SVD paths: `energy_policy.go`, `frd.go`, `freqrestest.go`, `loop_analysis.go`, `margin.go`, `norms.go`, `nyquist.go`, `passivity.go`, `rlocus.go`.\n- Simulation and time-response analysis outside standard-input batching: `gensig.go`, `simulate.go`, `response_info.go`; delay dispatch and general input histories remain cohesive and benchmark-stable.\n- Estimation, identification, local approximation, and random models: `covar.go`, `covariance_policy.go`, `ekf.go`, `linearize.go`, `local_approx.go`, `observer.go`, `random.go`, `sysid.go`.\n- Controllers and synthesis outside DARE: `controller.go`, `controller_observer_policy.go`, `h2syn.go`, `hinfsyn.go`, `lqg.go`, `pid.go`, `pidtune.go`, `smith.go`, `synthesis_partition.go`, `tunable.go`, `tuning.go`.\n- Collections and matrix utilities: `model_array.go`, `matutil.go`; ModelArray inherits the Step improvement without a second implementation.\n\nThis classification accounts for all 97 production Go files.\n\n## Accuracy evidence\n- SciPy/SLICOT-style generalized-pencil DARE tests cover singular A with cross term and ill-conditioned non-symmetric A; regular cases retain the existing ordered-Schur path.\n- Frequency kernels agree for non-symmetric continuous/discrete systems and all delay forms, with pole fallback preserving prior behavior.\n- Complex singular values agree with Gonum SVD for rectangular, rank-deficient, non-symmetric, `1e-150`, `1e150`, and vector `1e200` inputs.\n- Batched Step/Impulse agree with independent `Simulate` calls for continuous/discrete non-symmetric MIMO, feedthrough, gains, and delays.\n- Physical assembly mutation tests prove source/result matrix and metadata independence.\n\n## Benchmark evidence\n- DARE regular path, 10 samples: N10 statistically unchanged; N50 +2.37% time, -0.26% bytes, unchanged allocations. The generalized path is condition-selected for previously unsupported inputs.\n- Short frequency sweep, 10 samples: -42.19% time, -85.74% bytes, -94.44% allocations; long/delayed/LFT/B747 paths unchanged.\n- MIMO tuning, 10 samples: -52.11% to -53.25% time and -66.07% to -80.91% allocations.\n- MIMO Step/Impulse, 10 samples: about -50.3% time, -70.51% bytes, -66.49% allocations; SISO allocation counts unchanged.\n- Eight-component physical assembly, 10 samples: -30.96% time, -38.00% bytes, -29.00% allocations.\n- Representative full sweep: 234 benchmarks passed in 340.820s. Against the 217-benchmark `origin/main` baseline, no shared benchmark had a runtime regression above 30%; changes above that threshold were improvements in MIMO tuning, Step/Impulse, and ModelArray Step.\n\n## Validation gates\n- `go fix ./...` passed.\n- `go vet ./...` passed.\n- `go test -v -count=1` passed.\n- Task-specific 10-sample before/after benchmarks passed through benchstat.\n- Representative full benchmark sweep passed and was reviewed through benchstat.\n- `go.mod` and `go.sum` are unchanged, so `go vulncheck` was not required.\n- Five coherent implementation commits are present after synchronized `origin/main`.\n","ts":"2026-07-19T05:23:22.131093Z"}} +{"type":"state","ts":"2026-07-19T05:23:28.783951Z","data":{"id":"7OSGT6","state":"done","ts":"2026-07-19T05:23:28.783951Z"}} diff --git a/bench_test.go b/bench_test.go index 93b3ed4..9749528 100644 --- a/bench_test.go +++ b/bench_test.go @@ -273,6 +273,53 @@ func BenchmarkFreqResponse(b *testing.B) { } } +func BenchmarkFreqResponse_ShortSweep(b *testing.B) { + sys := benchSys(10, 2, 3) + omega := logspace(-2, 2, 8) + b.ResetTimer() + for b.Loop() { + sys.FreqResponse(omega) + } +} + +func BenchmarkFrequencySweepKernels(b *testing.B) { + tests := []struct { + name string + n, m, p, nw int + }{ + {name: "N4_W2", n: 4, m: 2, p: 2, nw: 2}, + {name: "N4_W8", n: 4, m: 2, p: 2, nw: 8}, + {name: "N4_W32", n: 4, m: 2, p: 2, nw: 32}, + {name: "N4_W100", n: 4, m: 2, p: 2, nw: 100}, + {name: "N10_W2", n: 10, m: 2, p: 3, nw: 2}, + {name: "N10_W8", n: 10, m: 2, p: 3, nw: 8}, + {name: "N10_W32", n: 10, m: 2, p: 3, nw: 32}, + {name: "N10_W100", n: 10, m: 2, p: 3, nw: 100}, + } + for _, test := range tests { + sys := benchSys(test.n, test.m, test.p) + evaluator := newFrequencyEvaluator(sys) + omega := logspace(-2, 2, test.nw) + size := test.nw * test.p * test.m + b.Run(test.name+"/StateSpace", func(b *testing.B) { + for b.Loop() { + data := make([]complex128, size) + if err := evaluator.evalStateSpaceSweepInto(omega, data); err != nil { + b.Fatal(err) + } + } + }) + b.Run(test.name+"/TransferFunction", func(b *testing.B) { + for b.Loop() { + data := make([]complex128, size) + if err := evaluator.evalTransferFunctionSweepInto(omega, data); err != nil { + b.Fatal(err) + } + } + }) + } +} + func BenchmarkBode(b *testing.B) { sys := benchSys(10, 2, 3) b.ResetTimer() diff --git a/frequency.go b/frequency.go index f1f0e70..f821cf9 100644 --- a/frequency.go +++ b/frequency.go @@ -144,6 +144,8 @@ type frequencyEvaluator struct { p int } +const directFrequencySweepWorkLimit = 80 + func newFrequencyEvaluator(sys *System) frequencyEvaluator { n, m, p := sys.Dims() return frequencyEvaluator{sys: sys, n: n, m: m, p: p} @@ -171,33 +173,21 @@ func (e frequencyEvaluator) response(omega []float64) (*FreqResponseMatrix, erro pm := e.p * e.m data := make([]complex128, nw*pm) if e.sys.IsDescriptor() { - ws := newSSEvalWorkspace(e.n, e.p, e.m) - for k, w := range omega { - s := e.sAt(w) - if err := evalFrSSInto(ws, e.sys, s, e.n, e.p, e.m); err != nil { - return nil, err - } - copy(data[k*pm:(k+1)*pm], ws.g[:pm]) + if err := e.evalStateSpaceSweepInto(omega, data); err != nil { + return nil, err } applyIODelayPhase(e.sys, omega, data, e.p, e.m, true) return e.matrix(data, omega), nil } - if nw == 1 { - s := e.sAt(omega[0]) - if err := e.evalStateSpaceInto(s, data); err != nil { - return nil, err + if e.useStateSpaceSweep(nw) { + if err := e.evalStateSpaceSweepInto(omega, data); err == nil { + applyIODelayPhase(e.sys, omega, data, e.p, e.m, true) + return e.matrix(data, omega), nil } - applyIODelayAtS(e.sys, s, data, e.p, e.m, true) - return e.matrix(data, omega), nil } - - res, err := e.sys.TransferFunction(nil) - if err != nil { + if err := e.evalTransferFunctionSweepInto(omega, data); err != nil { return nil, err } - for k, w := range omega { - res.TF.evalInto(e.sAt(w), data[k*pm:(k+1)*pm]) - } applyIODelayPhase(e.sys, omega, data, e.p, e.m, false) return e.matrix(data, omega), nil } @@ -238,6 +228,37 @@ func (e frequencyEvaluator) evalStateSpaceInto(s complex128, dst []complex128) e return nil } +func (e frequencyEvaluator) useStateSpaceSweep(nw int) bool { + if nw <= 1 || e.n == 0 { + return true + } + return nw <= directFrequencySweepWorkLimit/e.n +} + +func (e frequencyEvaluator) evalStateSpaceSweepInto(omega []float64, dst []complex128) error { + pm := e.p * e.m + ws := newSSEvalWorkspace(e.n, e.p, e.m) + for k, w := range omega { + if err := evalFrSSInto(ws, e.sys, e.sAt(w), e.n, e.p, e.m); err != nil { + return err + } + copy(dst[k*pm:(k+1)*pm], ws.g[:pm]) + } + return nil +} + +func (e frequencyEvaluator) evalTransferFunctionSweepInto(omega []float64, dst []complex128) error { + res, err := e.sys.TransferFunction(nil) + if err != nil { + return err + } + pm := e.p * e.m + for k, w := range omega { + res.TF.evalInto(e.sAt(w), dst[k*pm:(k+1)*pm]) + } + return nil +} + func (e frequencyEvaluator) matrix(data []complex128, omega []float64) *FreqResponseMatrix { return newFreqResponseMatrix(data, omega, e.p, e.m, e.sys.InputName, e.sys.OutputName) } diff --git a/frequency_svd.go b/frequency_svd.go index 3227655..b7745e3 100644 --- a/frequency_svd.go +++ b/frequency_svd.go @@ -12,6 +12,7 @@ type complexSVDWorkspace struct { block []float64 eig []float64 work []float64 + scale float64 nSV int gramN int blockN int @@ -54,20 +55,35 @@ func (ws *complexSVDWorkspace) singularValuesFromNested(dst []float64, data [][] func (ws *complexSVDWorkspace) fillBlock(at func(i, j int) complex128, p, m int) { n := ws.gramN stride := ws.blockN - for i := range ws.block { - ws.block[i] = 0 + clear(ws.block) + ws.scale = 0 + for i := range p { + for j := range m { + magnitude := cmplx.Abs(at(i, j)) + if math.IsNaN(magnitude) || magnitude > ws.scale { + ws.scale = magnitude + } + } + } + if ws.scale == 0 { + return } + invScale := complex(1/ws.scale, 0) for a := range n { for b := range n { var g complex128 if ws.useCol { for row := range p { - g += cmplx.Conj(at(row, a)) * at(row, b) + left := at(row, a) * invScale + right := at(row, b) * invScale + g += cmplx.Conj(left) * right } } else { for col := range m { - g += at(a, col) * cmplx.Conj(at(b, col)) + left := at(a, col) * invScale + right := at(b, col) * invScale + g += left * cmplx.Conj(right) } } @@ -81,8 +97,7 @@ func (ws *complexSVDWorkspace) fillBlock(at func(i, j int) complex128, p, m int) } func (ws *complexSVDWorkspace) singularValues(dst []float64) { - ok := impl.Dsyev(lapack.EVNone, blas.Upper, ws.blockN, ws.block, ws.blockN, ws.eig, ws.work, len(ws.work)) - if !ok { + if !ws.factorize() { for i := range dst { dst[i] = math.NaN() } @@ -94,8 +109,63 @@ func (ws *complexSVDWorkspace) singularValues(dst []float64) { } for i := range dst { lambda := nonnegativeGramEigenvalue(ws.eig[ws.blockN-1-2*i], scale) - dst[i] = math.Sqrt(lambda) + dst[i] = ws.scale * math.Sqrt(lambda) + } +} + +func (ws *complexSVDWorkspace) maximumFromFlat(data []complex128, base, p, m int) (float64, bool) { + if p == 0 || m == 0 { + return 0, true + } + values := data[base : base+p*m] + if p == 1 || m == 1 { + var norm float64 + for _, value := range values { + norm = math.Hypot(norm, cmplx.Abs(value)) + } + return norm, true + } + if p == 2 && m == 2 { + return maximumComplex2x2SingularValue(values), true + } + ws.fillBlock(func(i, j int) complex128 { + return data[base+i*m+j] + }, p, m) + if !ws.factorize() { + return 0, false } + scale := 1.0 + for _, lambda := range ws.eig { + scale = max(scale, math.Abs(lambda)) + } + lambda := nonnegativeGramEigenvalue(ws.eig[ws.blockN-1], scale) + return ws.scale * math.Sqrt(lambda), true +} + +func (ws *complexSVDWorkspace) factorize() bool { + return impl.Dsyev(lapack.EVNone, blas.Upper, ws.blockN, ws.block, ws.blockN, ws.eig, ws.work, len(ws.work)) +} + +func maximumComplex2x2SingularValue(data []complex128) float64 { + scale := 0.0 + for _, value := range data[:4] { + scale = math.Max(scale, cmplx.Abs(value)) + } + if scale == 0 { + return 0 + } + a00 := data[0] / complex(scale, 0) + a01 := data[1] / complex(scale, 0) + a10 := data[2] / complex(scale, 0) + a11 := data[3] / complex(scale, 0) + frobeniusSquared := complexMagnitudeSquared(a00) + complexMagnitudeSquared(a01) + complexMagnitudeSquared(a10) + complexMagnitudeSquared(a11) + determinantSquared := complexMagnitudeSquared(a00*a11 - a01*a10) + discriminant := math.Max(0, frobeniusSquared*frobeniusSquared-4*determinantSquared) + return scale * math.Sqrt((frobeniusSquared+math.Sqrt(discriminant))/2) +} + +func complexMagnitudeSquared(value complex128) float64 { + return real(value)*real(value) + imag(value)*imag(value) } func nonnegativeGramEigenvalue(lambda, scale float64) float64 { diff --git a/frequency_svd_test.go b/frequency_svd_test.go new file mode 100644 index 0000000..3a75a38 --- /dev/null +++ b/frequency_svd_test.go @@ -0,0 +1,112 @@ +package controlsys + +import ( + "math" + "testing" + + "gonum.org/v1/gonum/mat" +) + +func TestComplexSVDWorkspaceMaximumMatchesGonum(t *testing.T) { + tests := []struct { + name string + p, m int + data []complex128 + }{ + { + name: "Rectangular", + p: 3, m: 2, + data: []complex128{1 + 2i, -3 + 0.5i, 0.25 - 4i, 2 + 1i, -1.5 + 0.75i, 0.2 - 0.1i}, + }, + { + name: "RankDeficient", + p: 3, m: 3, + data: []complex128{ + 1 + 2i, 2 - 1i, -0.5 + 0.25i, + 2 + 4i, 4 - 2i, -1 + 0.5i, + -1 - 2i, -2 + 1i, 0.5 - 0.25i, + }, + }, + { + name: "NonSymmetric2x2", + p: 2, m: 2, + data: []complex128{1e-150 + 2e-150i, -3e-150 + 0.5e-150i, 4e-150 - 2e-150i, 0.25e-150 + 1e-150i}, + }, + { + name: "LargeScale", + p: 3, m: 2, + data: []complex128{1e150 + 2e150i, -3e150, 0.5e150 - 0.25e150i, 2e150 + 1e150i, -1e150, 0.75e150i}, + }, + { + name: "Vector", + p: 1, m: 4, + data: []complex128{1e200, 2e200i, -3e200 + 0.5e200i, 0.25e200 - 0.75e200i}, + }, + } + for _, test := range tests { + t.Run(test.name, func(t *testing.T) { + var workspace *complexSVDWorkspace + if test.p > 1 && test.m > 1 && (test.p != 2 || test.m != 2) { + workspace = newComplexSVDWorkspace(test.p, test.m) + } + got, ok := workspace.maximumFromFlat(test.data, 0, test.p, test.m) + if !ok { + t.Fatal("maximum singular-value decomposition failed") + } + want := gonumComplexSingularValues(t, test.data, test.p, test.m)[0] + if !sameRelative(got, want, 2e-12) { + t.Fatalf("maximum singular value = %.17g, want %.17g", got, want) + } + }) + } +} + +func TestComplexSVDWorkspaceScaledSingularValues(t *testing.T) { + p, m := 4, 3 + data := []complex128{ + 1e150 + 2e150i, -3e150, 0.5e150 - 0.25e150i, + 2e150 + 1e150i, -1e150, 0.75e150i, + -0.5e150, 0.2e150 + 0.1e150i, 4e150 - 2e150i, + 1.5e150 - 0.5e150i, -2e150 + 3e150i, 0.25e150, + } + want := gonumComplexSingularValues(t, data, p, m) + got := make([]float64, min(p, m)) + workspace := newComplexSVDWorkspace(p, m) + workspace.singularValuesFromFlat(got, data, 0, p, m) + for i := range got { + if !sameRelative(got[i], want[i], 3e-12) { + t.Fatalf("singular value %d = %.17g, want %.17g", i, got[i], want[i]) + } + } +} + +func gonumComplexSingularValues(t *testing.T, data []complex128, p, m int) []float64 { + t.Helper() + realForm := mat.NewDense(2*p, 2*m, nil) + for i := range p { + for j := range m { + value := data[i*m+j] + realForm.Set(i, j, real(value)) + realForm.Set(i, j+m, -imag(value)) + realForm.Set(i+p, j, imag(value)) + realForm.Set(i+p, j+m, real(value)) + } + } + var decomposition mat.SVD + if ok := decomposition.Factorize(realForm, mat.SVDNone); !ok { + t.Fatal("Gonum SVD factorization failed") + } + values := decomposition.Values(nil) + result := make([]float64, min(p, m)) + for i := range result { + result[i] = values[2*i] + } + return result +} + +func sameRelative(got, want, tolerance float64) bool { + if got == want { + return true + } + return math.Abs(got-want) <= tolerance*math.Max(math.Abs(got), math.Abs(want)) +} diff --git a/frequency_test.go b/frequency_test.go index ce32c95..7277bd5 100644 --- a/frequency_test.go +++ b/frequency_test.go @@ -270,6 +270,88 @@ func TestEvalFr_MatchesFreqResponse(t *testing.T) { } } +func TestFrequencyEvaluatorSweepKernelParity(t *testing.T) { + sys, err := New( + mat.NewDense(4, 4, []float64{ + -0.8, 0.4, 0, 0.1, + -0.2, -1.1, 0.3, 0, + 0.1, 0, -1.5, 0.2, + 0, -0.1, 0.25, -2, + }), + mat.NewDense(4, 2, []float64{1, 0.2, 0, 1, 0.4, -0.3, 0.1, 0.5}), + mat.NewDense(2, 4, []float64{1, 0, 0.3, -0.2, 0.1, 0.7, 0, 0.4}), + mat.NewDense(2, 2, []float64{0.1, 0, -0.05, 0.2}), + 0, + ) + if err != nil { + t.Fatal(err) + } + sys.Delay = mat.NewDense(2, 2, []float64{0.02, 0.04, 0.01, 0.03}) + sys.InputDelay = []float64{0.05, 0.02} + sys.OutputDelay = []float64{0.01, 0.03} + discrete := sys.Copy() + discrete.Dt = 0.1 + discrete.Delay = mat.NewDense(2, 2, []float64{1, 2, 3, 1}) + discrete.InputDelay = []float64{2, 1} + discrete.OutputDelay = []float64{1, 3} + + tests := []struct { + name string + system *System + omega []float64 + wantStateSpace bool + }{ + {name: "StateSpace", system: sys, omega: logspace(-2, 2, 8), wantStateSpace: true}, + {name: "TransferFunction", system: sys, omega: logspace(-2, 2, 100), wantStateSpace: false}, + {name: "DiscreteTransferFunction", system: discrete, omega: logspace(-2, 1, 100), wantStateSpace: false}, + } + for _, test := range tests { + t.Run(test.name, func(t *testing.T) { + evaluator := newFrequencyEvaluator(test.system) + if got := evaluator.useStateSpaceSweep(len(test.omega)); got != test.wantStateSpace { + t.Fatalf("useStateSpaceSweep = %t, want %t", got, test.wantStateSpace) + } + response, err := test.system.FreqResponse(test.omega) + if err != nil { + t.Fatal(err) + } + for k, w := range test.omega { + want, err := test.system.EvalFr(evaluator.sAt(w)) + if err != nil { + t.Fatal(err) + } + for i := range response.P { + for j := range response.M { + if diff := cmplx.Abs(response.At(k, i, j) - want[i][j]); diff > 2e-9 { + t.Fatalf("w=%g output=%d input=%d diff=%g", w, i, j, diff) + } + } + } + } + }) + } +} + +func TestFreqResponse_ShortSweepFallsBackAtPole(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, nil), + 0, + ) + if err != nil { + t.Fatal(err) + } + response, err := sys.FreqResponse([]float64{0, 1}) + if err != nil { + t.Fatal(err) + } + if got := response.At(1, 0, 0); cmplx.Abs(got+1i) > 1e-12 { + t.Fatalf("response at w=1 = %v, want -1i", got) + } +} + func TestFreqResponse_Empty(t *testing.T) { sys, err := New( mat.NewDense(1, 1, []float64{-1}), diff --git a/hardening_test.go b/hardening_test.go index 8d8a399..f1cbbb8 100644 --- a/hardening_test.go +++ b/hardening_test.go @@ -12,7 +12,7 @@ import ( func TestStateSpace_RejectsImproperTF(t *testing.T) { tf := &TransferFunc{ Num: [][][]float64{{{1, 0}}}, // s (degree 1) - Den: [][]float64{{1}}, // 1 (degree 0) + Den: [][]float64{{1}}, // 1 (degree 0) } _, err := tf.StateSpace(nil) if !errors.Is(err, ErrImproperTF) { @@ -224,7 +224,7 @@ func TestDare_GeneralizedWithCrossTerm(t *testing.T) { t.Fatal(err) } - if r := dareResidual(A, B, Q, R, res.X); r > 1e-9 { + if r := dareResidualWithCrossTerm(A, B, Q, R, S, res.X); r > 1e-9 { t.Errorf("DARE residual = %e", r) } checkSymmetric(t, res.X, 1e-10) @@ -258,6 +258,71 @@ func TestDare_IllConditionedA(t *testing.T) { } } +func TestDare_SingularA_SciPyOracle(t *testing.T) { + A := mat.NewDense(2, 2, []float64{0, 1, 0, 0.2}) + B := mat.NewDense(2, 1, []float64{0.2, 1}) + Q := mat.NewDense(2, 2, []float64{2, 0.1, 0.1, 1}) + R := mat.NewDense(1, 1, []float64{1.5}) + S := mat.NewDense(2, 1, []float64{0.1, 0.2}) + wantX := mat.NewDense(2, 2, []float64{ + 1.9977280744658075, 0.07198755827337247, + 0.07198755827337247, 2.7928485672615015, + }) + wantK := mat.NewDense(1, 2, []float64{0.022719255341908287, 0.2801244172662588}) + + res, err := Dare(A, B, Q, R, &RiccatiOpts{S: S}) + if err != nil { + t.Fatal(err) + } + if !mat.EqualApprox(res.X, wantX, 2e-12) { + t.Errorf("X differs from SciPy oracle:\ngot:\n%v\nwant:\n%v", mat.Formatted(res.X), mat.Formatted(wantX)) + } + if !mat.EqualApprox(res.K, wantK, 2e-12) { + t.Errorf("K differs from SciPy oracle:\ngot:\n%v\nwant:\n%v", mat.Formatted(res.K), mat.Formatted(wantK)) + } + if residual := dareResidualWithCrossTerm(A, B, Q, R, S, res.X); residual > 2e-14 { + t.Errorf("residual = %e", residual) + } + checkSymmetric(t, res.X, 1e-13) + for _, eig := range res.Eig { + if cmplx.Abs(eig) >= 1 { + t.Errorf("non-stable closed-loop eigenvalue: %v", eig) + } + } +} + +func TestDare_IllConditionedNonSymmetricA_SciPyOracle(t *testing.T) { + A := mat.NewDense(2, 2, []float64{1.05, 0.4, 0, 1e-12}) + B := mat.NewDense(2, 1, []float64{1, 0.3}) + Q := mat.NewDense(2, 2, []float64{2, 0.25, 0.25, 1}) + R := mat.NewDense(1, 1, []float64{0.75}) + wantX := mat.NewDense(2, 2, []float64{ + 2.64626274734184, 0.49619533232058405, + 0.49619533232058405, 1.0937886980268443, + }) + wantK := mat.NewDense(1, 2, []float64{0.7738796586414393, 0.2948112985302894}) + + res, err := Dare(A, B, Q, R, nil) + if err != nil { + t.Fatal(err) + } + if !mat.EqualApprox(res.X, wantX, 2e-12) { + t.Errorf("X differs from SciPy oracle:\ngot:\n%v\nwant:\n%v", mat.Formatted(res.X), mat.Formatted(wantX)) + } + if !mat.EqualApprox(res.K, wantK, 2e-12) { + t.Errorf("K differs from SciPy oracle:\ngot:\n%v\nwant:\n%v", mat.Formatted(res.K), mat.Formatted(wantK)) + } + if residual := dareResidual(A, B, Q, R, res.X); residual > 1e-13 { + t.Errorf("residual = %e", residual) + } + checkSymmetric(t, res.X, 1e-13) + for _, eig := range res.Eig { + if cmplx.Abs(eig) >= 1 { + t.Errorf("non-stable closed-loop eigenvalue: %v", eig) + } + } +} + func TestCare_MATLABValidated_3x3(t *testing.T) { A := mat.NewDense(3, 3, []float64{ -1.017, -0.224, 0.043, diff --git a/physical_assembly.go b/physical_assembly.go index bf7d590..c50dba9 100644 --- a/physical_assembly.go +++ b/physical_assembly.go @@ -92,7 +92,7 @@ func newPhysicalAssemblyPlan(name string, components []PhysicalComponent, connec } plan := &physicalAssemblyPlan{ name: name, - components: make([]PhysicalComponent, len(components)), + components: components, stateOffsets: make([]int, len(components)), inputOffsets: make([]int, len(components)), outputOffsets: make([]int, len(components)), @@ -118,7 +118,6 @@ func newPhysicalAssemblyPlan(name string, components []PhysicalComponent, connec plan.totalStates += n plan.totalInputs += m plan.totalOutputs += p - plan.components[i] = NewPhysicalComponent(component.Name, component.System, component.Ports) byComponent[component.Name] = i } portIndex, err := plan.bindPorts() @@ -288,12 +287,15 @@ func (p *physicalAssemblyPlan) assemble() (*System, error) { bAugmented := newDense(nAugmented, len(p.externalInputs)) setBlock(aAugmented, 0, 0, a) setBlock(eAugmented, 0, 0, e) + aAugmentedRaw := aAugmented.RawMatrix() + bAugmentedRaw := bAugmented.RawMatrix() + bRaw := b.RawMatrix() for state := range n { for position, channel := range p.internalInputs { - aAugmented.Set(state, n+position, b.At(state, channel)) + aAugmentedRaw.Data[state*aAugmentedRaw.Stride+n+position] = bRaw.Data[state*bRaw.Stride+channel] } for position, channel := range p.externalInputs { - bAugmented.Set(state, position, b.At(state, channel)) + bAugmentedRaw.Data[state*bAugmentedRaw.Stride+position] = bRaw.Data[state*bRaw.Stride+channel] } } constraint := 0 @@ -321,7 +323,8 @@ func (p *physicalAssemblyPlan) assemble() (*System, error) { row := n + constraint for _, portIndex := range group.ports { channel := p.ports[portIndex].inputs[coordinate] - aAugmented.Set(row, n+p.internalInputPos[channel], aAugmented.At(row, n+p.internalInputPos[channel])+1) + column := n + p.internalInputPos[channel] + aAugmentedRaw.Data[row*aAugmentedRaw.Stride+column]++ } constraint++ } @@ -331,15 +334,19 @@ func (p *physicalAssemblyPlan) assemble() (*System, error) { } cAugmented := newDense(len(p.externalOutputs), nAugmented) dAugmented := newDense(len(p.externalOutputs), len(p.externalInputs)) + cRaw := c.RawMatrix() + dRaw := d.RawMatrix() + cAugmentedRaw := cAugmented.RawMatrix() + dAugmentedRaw := dAugmented.RawMatrix() for row, output := range p.externalOutputs { for state := range n { - cAugmented.Set(row, state, c.At(output, state)) + cAugmentedRaw.Data[row*cAugmentedRaw.Stride+state] = cRaw.Data[output*cRaw.Stride+state] } for position, input := range p.internalInputs { - cAugmented.Set(row, n+position, d.At(output, input)) + cAugmentedRaw.Data[row*cAugmentedRaw.Stride+n+position] = dRaw.Data[output*dRaw.Stride+input] } for column, input := range p.externalInputs { - dAugmented.Set(row, column, d.At(output, input)) + dAugmentedRaw.Data[row*dAugmentedRaw.Stride+column] = dRaw.Data[output*dRaw.Stride+input] } } var constructorB, constructorC, constructorD *mat.Dense @@ -352,7 +359,7 @@ func (p *physicalAssemblyPlan) assemble() (*System, error) { if len(p.externalInputs) > 0 && len(p.externalOutputs) > 0 { constructorD = dAugmented } - result, err := NewDescriptor(aAugmented, constructorB, constructorC, constructorD, eAugmented, p.components[0].System.Dt) + result, err := newDescriptorOwned(aAugmented, constructorB, constructorC, constructorD, eAugmented, p.components[0].System.Dt) if err != nil { return nil, err } @@ -365,14 +372,19 @@ func (p *physicalAssemblyPlan) assemble() (*System, error) { func (p *physicalAssemblyPlan) addAcrossConstraint(a, b, c, d *mat.Dense, row, portIndex, coordinate int, factor float64) { output := p.ports[portIndex].outputs[coordinate] + aRaw := a.RawMatrix() + bRaw := b.RawMatrix() + cRaw := c.RawMatrix() + dRaw := d.RawMatrix() for state := range p.totalStates { - a.Set(row, state, a.At(row, state)+factor*c.At(output, state)) + aRaw.Data[row*aRaw.Stride+state] += factor * cRaw.Data[output*cRaw.Stride+state] } for position, input := range p.internalInputs { - a.Set(row, p.totalStates+position, a.At(row, p.totalStates+position)+factor*d.At(output, input)) + column := p.totalStates + position + aRaw.Data[row*aRaw.Stride+column] += factor * dRaw.Data[output*dRaw.Stride+input] } for position, input := range p.externalInputs { - b.Set(row, position, b.At(row, position)+factor*d.At(output, input)) + bRaw.Data[row*bRaw.Stride+position] += factor * dRaw.Data[output*dRaw.Stride+input] } } @@ -382,6 +394,7 @@ func (p *physicalAssemblyPlan) aggregateMatrices() (a, b, c, d, e *mat.Dense) { c = newDense(p.totalOutputs, p.totalStates) d = newDense(p.totalOutputs, p.totalInputs) e = newDense(p.totalStates, p.totalStates) + eRaw := e.RawMatrix() for i, component := range p.components { n, _, _ := component.System.Dims() stateOffset := p.stateOffsets[i] @@ -393,7 +406,7 @@ func (p *physicalAssemblyPlan) aggregateMatrices() (a, b, c, d, e *mat.Dense) { setBlock(d, outputOffset, inputOffset, component.System.D) if component.System.E == nil { for state := range n { - e.Set(stateOffset+state, stateOffset+state, 1) + eRaw.Data[(stateOffset+state)*eRaw.Stride+stateOffset+state] = 1 } } else { setBlock(e, stateOffset, stateOffset, component.System.E) diff --git a/physical_assembly_test.go b/physical_assembly_test.go index 1548805..d16cf90 100644 --- a/physical_assembly_test.go +++ b/physical_assembly_test.go @@ -107,6 +107,35 @@ func TestPhysicalAssemblyPreservesUnconnectedDescriptorComponent(t *testing.T) { } } +func TestPhysicalAssemblyOwnsConnectedResult(t *testing.T) { + left := physicalCoupledComponent(t, "left", 1) + right := physicalCoupledComponent(t, "right", 2) + assembled, err := AssemblePhysical("pair", []PhysicalComponent{left, right}, []PhysicalConnection{ + {FromComponent: "left", FromPort: "mount", ToComponent: "right", ToPort: "mount"}, + }) + if err != nil { + t.Fatal(err) + } + + assembledA := assembled.A.At(0, 0) + assembledInputName := assembled.InputName[0] + left.System.A.Set(0, 0, 99) + left.System.InputName[0] = "changed" + left.Ports[0].Name = "changed" + left.Ports[0].Input[0] = 0 + if assembled.A.At(0, 0) != assembledA || assembled.InputName[0] != assembledInputName { + t.Fatal("source mutation changed assembled result") + } + + sourceA := left.System.A.At(0, 0) + sourceInputName := left.System.InputName[0] + assembled.A.Set(0, 0, -123) + assembled.InputName[0] = "result.changed" + if left.System.A.At(0, 0) != sourceA || left.System.InputName[0] != sourceInputName { + t.Fatal("assembled result mutation changed source component") + } +} + func TestPhysicalAssemblyImplicitPortsSkipExplicitBindings(t *testing.T) { explicit := PhysicalPort{Name: "external", Kind: PhysicalPortDisplacement, Dimension: 1, Input: []int{0}, Output: []int{0}} implicit := PhysicalPort{Name: "mount", Kind: PhysicalPortDisplacement, Dimension: 2} diff --git a/response.go b/response.go index c91e691..41b03df 100644 --- a/response.go +++ b/response.go @@ -27,6 +27,13 @@ type timeResponsePlanner struct { sys *System } +type standardInputResponse uint8 + +const ( + stepResponse standardInputResponse = iota + impulseResponse +) + func newTimeResponsePlanner(sys *System) timeResponsePlanner { return timeResponsePlanner{sys: sys} } @@ -173,6 +180,119 @@ func (p timeResponsePlan) response(Y *mat.Dense) *TimeResponse { return &TimeResponse{T: p.t, Y: Y, OutputName: copyStringSlice(p.original.OutputName)} } +func (p timeResponsePlan) allInputResponse(kind standardInputResponse) (*TimeResponse, error) { + _, m, outputs := p.sim.Dims() + if m <= 1 || outputs == 0 || p.sim.HasDelay() || p.sim.IsDescriptor() { + return p.simulatedInputResponse(kind) + } + if err := p.sim.Validate(); err != nil { + return nil, fmt.Errorf("input 0: %w", err) + } + return p.batchedInputResponse(kind), nil +} + +func (p timeResponsePlan) simulatedInputResponse(kind standardInputResponse) (*TimeResponse, error) { + _, m, outputs := p.sim.Dims() + Y := mat.NewDense(outputs*m, p.steps, nil) + amplitude := kind.amplitude(p) + + for input := range m { + u := mat.NewDense(m, p.steps, nil) + if kind == stepResponse { + uRaw := u.RawMatrix() + for sample := range p.steps { + uRaw.Data[input*uRaw.Stride+sample] = amplitude + } + } else { + u.Set(input, 0, amplitude) + } + resp, err := p.sim.Simulate(u, nil, nil) + if err != nil { + return nil, fmt.Errorf("input %d: %w", input, err) + } + if resp.Y == nil { + continue + } + yRaw := Y.RawMatrix() + responseRaw := resp.Y.RawMatrix() + for output := range outputs { + dstBase := (input*outputs + output) * yRaw.Stride + srcBase := output * responseRaw.Stride + copy(yRaw.Data[dstBase:dstBase+p.steps], responseRaw.Data[srcBase:srcBase+p.steps]) + } + } + + return p.response(Y), nil +} + +func (p timeResponsePlan) batchedInputResponse(kind standardInputResponse) *TimeResponse { + states, inputs, outputs := p.sim.Dims() + Y := mat.NewDense(outputs*inputs, p.steps, nil) + yRaw := Y.RawMatrix() + dRaw := p.sim.D.RawMatrix() + amplitude := kind.amplitude(p) + + if states == 0 { + for sample := range p.steps { + if kind == impulseResponse && sample > 0 { + break + } + for input := range inputs { + for output := range outputs { + yRaw.Data[(input*outputs+output)*yRaw.Stride+sample] = amplitude * dRaw.Data[output*dRaw.Stride+input] + } + } + } + return p.response(Y) + } + + x := make([]float64, states*inputs) + nextX := make([]float64, states*inputs) + aRaw := p.sim.A.RawMatrix() + bRaw := p.sim.B.RawMatrix() + cRaw := p.sim.C.RawMatrix() + + for sample := range p.steps { + inputActive := kind == stepResponse || sample == 0 + for input := range inputs { + xInput := x[input*states : (input+1)*states] + nextInput := nextX[input*states : (input+1)*states] + for output := range outputs { + value := 0.0 + if inputActive { + value = amplitude * dRaw.Data[output*dRaw.Stride+input] + } + cRow := cRaw.Data[output*cRaw.Stride : output*cRaw.Stride+states] + for state, coefficient := range cRow { + value += coefficient * xInput[state] + } + yRaw.Data[(input*outputs+output)*yRaw.Stride+sample] = value + } + for state := range states { + value := 0.0 + if inputActive { + value = amplitude * bRaw.Data[state*bRaw.Stride+input] + } + aRow := aRaw.Data[state*aRaw.Stride : state*aRaw.Stride+states] + for column, coefficient := range aRow { + value += coefficient * xInput[column] + } + nextInput[state] = value + } + } + x, nextX = nextX, x + } + + return p.response(Y) +} + +func (kind standardInputResponse) amplitude(p timeResponsePlan) float64 { + if kind == impulseResponse && p.wasContinuous { + return 1 / p.dt + } + return 1 +} + func prepareLsimResponse(sys *System, u *mat.Dense, t []float64) (timeResponsePlan, *mat.Dense, error) { return newTimeResponsePlanner(sys).lsim(u, t) } @@ -544,30 +664,11 @@ func Step(sys *System, tFinal float64) (*TimeResponse, error) { if err != nil { return nil, err } - - _, m, p := plan.sim.Dims() - rows := p * m - Y := mat.NewDense(rows, plan.steps, nil) - - for j := range m { - u := mat.NewDense(m, plan.steps, nil) - for k := 0; k < plan.steps; k++ { - u.Set(j, k, 1) - } - resp, err := plan.sim.Simulate(u, nil, nil) - if err != nil { - return nil, fmt.Errorf("Step: input %d: %w", j, err) - } - if resp.Y != nil { - for i := range p { - for k := 0; k < plan.steps; k++ { - Y.Set(j*p+i, k, resp.Y.At(i, k)) - } - } - } + resp, err := plan.allInputResponse(stepResponse) + if err != nil { + return nil, fmt.Errorf("Step: %w", err) } - - return plan.response(Y), nil + return resp, nil } func Impulse(sys *System, tFinal float64) (*TimeResponse, error) { @@ -576,32 +677,11 @@ func Impulse(sys *System, tFinal float64) (*TimeResponse, error) { return nil, err } - _, m, p := plan.sim.Dims() - rows := p * m - Y := mat.NewDense(rows, plan.steps, nil) - - amp := 1.0 - if plan.wasContinuous { - amp = 1.0 / plan.dt - } - - for j := range m { - u := mat.NewDense(m, plan.steps, nil) - u.Set(j, 0, amp) - resp, err := plan.sim.Simulate(u, nil, nil) - if err != nil { - return nil, fmt.Errorf("Impulse: input %d: %w", j, err) - } - if resp.Y != nil { - for i := range p { - for k := 0; k < plan.steps; k++ { - Y.Set(j*p+i, k, resp.Y.At(i, k)) - } - } - } + resp, err := plan.allInputResponse(impulseResponse) + if err != nil { + return nil, fmt.Errorf("Impulse: %w", err) } - - return plan.response(Y), nil + return resp, nil } func Initial(sys *System, x0 *mat.VecDense, tFinal float64) (*TimeResponse, error) { diff --git a/response_batch_test.go b/response_batch_test.go new file mode 100644 index 0000000..a49c5af --- /dev/null +++ b/response_batch_test.go @@ -0,0 +1,123 @@ +package controlsys + +import ( + "math" + "testing" + + "gonum.org/v1/gonum/mat" +) + +func TestAllInputResponseMatchesIndependentSimulations(t *testing.T) { + continuous, err := New( + mat.NewDense(2, 2, []float64{-1.2, 0.7, -0.4, -2.1}), + mat.NewDense(2, 3, []float64{1, -0.5, 0.25, 0.3, 1.2, -0.8}), + mat.NewDense(2, 2, []float64{1.1, -0.2, 0.4, 0.9}), + mat.NewDense(2, 3, []float64{0.1, -0.3, 0.2, 0.5, 0.25, -0.4}), + 0, + ) + if err != nil { + t.Fatal(err) + } + discrete, err := New( + mat.NewDense(2, 2, []float64{0.7, 0.2, -0.1, 0.5}), + mat.NewDense(2, 3, []float64{1, -0.5, 0.25, 0.3, 1.2, -0.8}), + mat.NewDense(2, 2, []float64{1.1, -0.2, 0.4, 0.9}), + mat.NewDense(2, 3, []float64{0.1, -0.3, 0.2, 0.5, 0.25, -0.4}), + 0.2, + ) + if err != nil { + t.Fatal(err) + } + gain, err := NewGain(mat.NewDense(2, 3, []float64{1, -2, 0.5, 0.25, 3, -0.75}), 0) + if err != nil { + t.Fatal(err) + } + delayed := discrete.Copy() + delayed.Delay = mat.NewDense(2, 3, []float64{0, 1, 2, 2, 0, 1}) + + tests := []struct { + name string + sys *System + final float64 + kind standardInputResponse + public func(*System, float64) (*TimeResponse, error) + }{ + {name: "ContinuousStep", sys: continuous, final: 1.2, kind: stepResponse, public: Step}, + {name: "ContinuousImpulse", sys: continuous, final: 1.2, kind: impulseResponse, public: Impulse}, + {name: "DiscreteStep", sys: discrete, final: 2, kind: stepResponse, public: Step}, + {name: "DiscreteImpulse", sys: discrete, final: 2, kind: impulseResponse, public: Impulse}, + {name: "GainStep", sys: gain, final: 0.2, kind: stepResponse, public: Step}, + {name: "GainImpulse", sys: gain, final: 0.2, kind: impulseResponse, public: Impulse}, + {name: "DelayedStep", sys: delayed, final: 2, kind: stepResponse, public: Step}, + {name: "DelayedImpulse", sys: delayed, final: 2, kind: impulseResponse, public: Impulse}, + } + + for _, test := range tests { + t.Run(test.name, func(t *testing.T) { + got, err := test.public(test.sys, test.final) + if err != nil { + t.Fatal(err) + } + want := independentInputResponse(t, test.sys, test.final, test.kind) + compareTimeResponses(t, got, want, 2e-13) + }) + } +} + +func independentInputResponse(t *testing.T, sys *System, tFinal float64, kind standardInputResponse) *TimeResponse { + t.Helper() + plan, err := prepareAutoTimeResponse(sys, tFinal, 0) + if err != nil { + t.Fatal(err) + } + _, inputs, outputs := plan.sim.Dims() + Y := mat.NewDense(outputs*inputs, plan.steps, nil) + amplitude := kind.amplitude(plan) + for input := range inputs { + u := mat.NewDense(inputs, plan.steps, nil) + if kind == stepResponse { + for sample := range plan.steps { + u.Set(input, sample, amplitude) + } + } else { + u.Set(input, 0, amplitude) + } + response, err := plan.sim.Simulate(u, nil, nil) + if err != nil { + t.Fatal(err) + } + for output := range outputs { + for sample := range plan.steps { + Y.Set(input*outputs+output, sample, response.Y.At(output, sample)) + } + } + } + return plan.response(Y) +} + +func compareTimeResponses(t *testing.T, got, want *TimeResponse, tolerance float64) { + t.Helper() + if len(got.T) != len(want.T) { + t.Fatalf("time-grid lengths differ: got %d, want %d", len(got.T), len(want.T)) + } + for i := range got.T { + if got.T[i] != want.T[i] { + t.Fatalf("time %d = %.17g, want %.17g", i, got.T[i], want.T[i]) + } + } + gotRows, gotCols := got.Y.Dims() + wantRows, wantCols := want.Y.Dims() + if gotRows != wantRows || gotCols != wantCols { + t.Fatalf("response dimensions = %dx%d, want %dx%d", gotRows, gotCols, wantRows, wantCols) + } + for row := range gotRows { + for col := range gotCols { + gotValue := got.Y.At(row, col) + wantValue := want.Y.At(row, col) + scale := math.Max(1, math.Max(math.Abs(gotValue), math.Abs(wantValue))) + if math.Abs(gotValue-wantValue) > tolerance*scale { + t.Fatalf("response(%d,%d) = %.17g, want %.17g", row, col, gotValue, wantValue) + } + } + } +} diff --git a/riccati.go b/riccati.go index 09bde82..c3e428f 100644 --- a/riccati.go +++ b/riccati.go @@ -1,6 +1,8 @@ package controlsys import ( + "math" + "gonum.org/v1/gonum/blas" "gonum.org/v1/gonum/blas/blas64" "gonum.org/v1/gonum/lapack" @@ -8,67 +10,65 @@ import ( ) type RiccatiWorkspace struct { - rChol []float64 - aWork []float64 - qWork []float64 - rinvBt []float64 - rinvSt []float64 - g []float64 - h []float64 - wr []float64 - wi []float64 - vs []float64 - bwork []bool - work []float64 - u11 []float64 - u21 []float64 - ipiv []int - xData []float64 - eig []complex128 - kData []float64 - ait []float64 - aipiv []int - aiWork []float64 - aitq []float64 - gait []float64 - gaitq []float64 - z []float64 - btx []float64 - rbar []float64 - iwork []int + rChol []float64 + aWork []float64 + qWork []float64 + rinvBt []float64 + rinvSt []float64 + g []float64 + h []float64 + wr []float64 + wi []float64 + vs []float64 + bwork []bool + work []float64 + u11 []float64 + u21 []float64 + ipiv []int + xData []float64 + eig []complex128 + kData []float64 + z []float64 + beta []float64 + pencilH []float64 + pencilJ []float64 + pencilInput []float64 + tau []float64 + btx []float64 + rbar []float64 + iwork []int } func NewRiccatiWorkspace(n, m int) *RiccatiWorkspace { nn := 2 * n return &RiccatiWorkspace{ - rChol: make([]float64, m*m), - aWork: make([]float64, n*n), - qWork: make([]float64, n*n), - rinvBt: make([]float64, m*n), - rinvSt: make([]float64, m*n), - g: make([]float64, n*n), - h: make([]float64, nn*nn), - wr: make([]float64, nn), - wi: make([]float64, nn), - vs: make([]float64, nn*nn), - bwork: make([]bool, nn), - work: make([]float64, nn*50), - u11: make([]float64, n*n), - u21: make([]float64, n*n), - ipiv: make([]int, n), - xData: make([]float64, n*n), - eig: make([]complex128, n), - kData: make([]float64, m*n), - ait: make([]float64, n*n), - aipiv: make([]int, n), - aiWork: make([]float64, n*50), - aitq: make([]float64, n*n), - gait: make([]float64, n*n), - gaitq: make([]float64, n*n), - z: make([]float64, nn*nn), - btx: make([]float64, m*n), - rbar: make([]float64, m*m), - iwork: make([]int, n), + rChol: make([]float64, m*m), + aWork: make([]float64, n*n), + qWork: make([]float64, n*n), + rinvBt: make([]float64, m*n), + rinvSt: make([]float64, m*n), + g: make([]float64, n*n), + h: make([]float64, nn*nn), + wr: make([]float64, nn), + wi: make([]float64, nn), + vs: make([]float64, nn*nn), + bwork: make([]bool, nn), + work: make([]float64, nn*50), + u11: make([]float64, n*n), + u21: make([]float64, n*n), + ipiv: make([]int, n), + xData: make([]float64, n*n), + eig: make([]complex128, n), + kData: make([]float64, m*n), + z: make([]float64, nn*nn), + beta: make([]float64, nn), + pencilH: make([]float64, (nn+m)*nn), + pencilJ: make([]float64, (nn+m)*nn), + pencilInput: make([]float64, (nn+m)*m), + tau: make([]float64, m), + btx: make([]float64, m*n), + rbar: make([]float64, m*m), + iwork: make([]int, n), } } @@ -291,39 +291,135 @@ func Dare(A, B, Q, R *mat.Dense, opts *RiccatiOpts) (*RiccatiResult, error) { return nil, ErrSingularR } - aWork := ws.aWork[:n*n] aRaw := A.RawMatrix() + bRaw := B.RawMatrix() + + subspace, err := problem.discreteStableSubspace() + if err != nil { + return nil, err + } + vs := subspace.vectors + + // Extract U11 = vs[0:n, 0:n], U21 = vs[n:2n, 0:n] + u11 := ws.u11[:n*n] + u21 := ws.u21[:n*n] + copyStrided(u11, n, vs, 2*n, n, n) + copyBlock(u21, n, 0, 0, vs, 2*n, n, 0, n, n) + + work := ws.work + anorm := impl.Dlange(lapack.MaxColumnSum, n, n, u11, n, work[:n]) + ipiv := ws.ipiv[:n] + if !impl.Dgetrf(n, n, u11, n, ipiv) { + return nil, ErrNoStabilizing + } + + iwork2 := ws.iwork[:n] + rcnd := impl.Dgecon(lapack.MaxColumnSum, n, u11, n, anorm, work[:4*n], iwork2) + + xData := ws.xData[:n*n] + for i := range n { + for j := range n { + xData[i*n+j] = u21[j*n+i] + } + } + impl.Dgetrs(blas.Trans, n, n, u11, n, ipiv, xData, n) + + symmetrize(xData, n, n) + X := mat.NewDense(n, n, xData) + + // Closed-loop eigenvalues + eig := ws.eig[:n] + for i := range n { + eig[i] = complex(subspace.alphaR[i]/subspace.beta[i], subspace.alphaI[i]/subspace.beta[i]) + } + + // Gain K = (R + B'XB)⁻¹ * (B'XA + S') + btx := ws.btx[:m*n] + blas64.Gemm(blas.Trans, blas.NoTrans, + 1, blas64.General{Rows: n, Cols: m, Data: bRaw.Data, Stride: bRaw.Stride}, + blas64.General{Rows: n, Cols: n, Data: xData, Stride: n}, + 0, blas64.General{Rows: m, Cols: n, Data: btx, Stride: n}) + + rbar := ws.rbar[:m*m] + copyStrided(rbar, m, rRaw.Data, rRaw.Stride, m, m) + blas64.Gemm(blas.NoTrans, blas.NoTrans, + 1, blas64.General{Rows: m, Cols: n, Data: btx, Stride: n}, + blas64.General{Rows: n, Cols: m, Data: bRaw.Data, Stride: bRaw.Stride}, + 1, blas64.General{Rows: m, Cols: m, Data: rbar, Stride: m}) + + if !impl.Dpotrf(blas.Upper, m, rbar, m) { + return nil, ErrSingularR + } + + // BtXA = BtX * A (m×n) + kData := ws.kData[:m*n] + blas64.Gemm(blas.NoTrans, blas.NoTrans, + 1, blas64.General{Rows: m, Cols: n, Data: btx, Stride: n}, + blas64.General{Rows: n, Cols: n, Data: aRaw.Data, Stride: aRaw.Stride}, + 0, blas64.General{Rows: m, Cols: n, Data: kData, Stride: n}) + if S != nil { + sRaw := S.RawMatrix() + for j := range m { + row := kData[j*n:] + for i := range n { + row[i] += sRaw.Data[i*sRaw.Stride+j] + } + } + } + impl.Dpotrs(blas.Upper, m, n, rbar, m, kData, n) + K := mat.NewDense(m, n, kData) + + return &RiccatiResult{X: X, K: K, Eig: eig, Rcnd: rcnd}, nil +} + +type discreteRiccatiSubspace struct { + vectors []float64 + alphaR []float64 + alphaI []float64 + beta []float64 +} + +func (problem riccatiProblem) discreteStableSubspace() (discreteRiccatiSubspace, error) { + if subspace, suitable, err := problem.regularDiscreteStableSubspace(); suitable || err != nil { + return subspace, err + } + return problem.generalizedDiscreteStableSubspace() +} + +func (problem riccatiProblem) regularDiscreteStableSubspace() (subspace discreteRiccatiSubspace, suitable bool, err error) { + n, m, ws := problem.n, problem.m, problem.ws + nn := 2 * n + aRaw := problem.A.RawMatrix() + bRaw := problem.B.RawMatrix() + + aWork := ws.aWork[:n*n] copyStrided(aWork, n, aRaw.Data, aRaw.Stride, n, n) qWork := ws.qWork[:n*n] - qRaw := Q.RawMatrix() + qRaw := problem.Q.RawMatrix() copyStrided(qWork, n, qRaw.Data, qRaw.Stride, n, n) - bRaw := B.RawMatrix() - - // R⁻¹*B' (m×n) rinvBt := ws.rinvBt[:m*n] for i := range n { for j := range m { rinvBt[j*n+i] = bRaw.Data[i*bRaw.Stride+j] } } - impl.Dpotrs(blas.Upper, m, n, rChol, m, rinvBt, n) + impl.Dpotrs(blas.Upper, m, n, ws.rChol[:m*m], m, rinvBt, n) - if S != nil { + if problem.S != nil { rinvSt := ws.rinvSt[:m*n] - sRaw := S.RawMatrix() + sRaw := problem.S.RawMatrix() for i := range n { for j := range m { rinvSt[j*n+i] = sRaw.Data[i*sRaw.Stride+j] } } - impl.Dpotrs(blas.Upper, m, n, rChol, m, rinvSt, n) + impl.Dpotrs(blas.Upper, m, n, ws.rChol[:m*m], m, rinvSt, n) blas64.Gemm(blas.NoTrans, blas.NoTrans, -1, blas64.General{Rows: n, Cols: m, Data: bRaw.Data, Stride: bRaw.Stride}, blas64.General{Rows: m, Cols: n, Data: rinvSt, Stride: n}, 1, blas64.General{Rows: n, Cols: n, Data: aWork, Stride: n}) - blas64.Gemm(blas.NoTrans, blas.NoTrans, -1, blas64.General{Rows: n, Cols: m, Data: sRaw.Data, Stride: sRaw.Stride}, blas64.General{Rows: m, Cols: n, Data: rinvSt, Stride: n}, @@ -331,7 +427,6 @@ func Dare(A, B, Q, R *mat.Dense, opts *RiccatiOpts) (*RiccatiResult, error) { symmetrize(qWork, n, n) } - // G = B * R⁻¹ * B' g := ws.g[:n*n] blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, blas64.General{Rows: n, Cols: m, Data: bRaw.Data, Stride: bRaw.Stride}, @@ -339,49 +434,47 @@ func Dare(A, B, Q, R *mat.Dense, opts *RiccatiOpts) (*RiccatiResult, error) { 0, blas64.General{Rows: n, Cols: n, Data: g, Stride: n}) symmetrize(g, n, n) - // Form symplectic matrix Z = M⁻¹*L (requires A invertible): - // Z = [[A + G*Ait*Q, -G*Ait], [-Ait*Q, Ait]] where Ait = (A')⁻¹ - // Compute Ait = (A')⁻¹ via LU factorization of A' - ait := ws.ait[:n*n] + scratch := ws.pencilH + ait := scratch[:n*n] + aitq := scratch[n*n : 2*n*n] + gait := scratch[2*n*n : 3*n*n] + gaitq := scratch[3*n*n : 4*n*n] for i := range n { for j := range n { ait[i*n+j] = aWork[j*n+i] } } - aipiv := ws.aipiv[:n] - if !impl.Dgetrf(n, n, ait, n, aipiv) { - return nil, ErrNoStabilizing + work := ws.work + anorm := impl.Dlange(lapack.MaxColumnSum, n, n, ait, n, work[:n]) + ipiv := ws.ipiv[:n] + if !impl.Dgetrf(n, n, ait, n, ipiv) { + return discreteRiccatiSubspace{}, false, nil + } + rcnd := impl.Dgecon(lapack.MaxColumnSum, n, ait, n, anorm, work[:4*n], ws.iwork[:n]) + if rcnd < math.Sqrt(eps()) { + return discreteRiccatiSubspace{}, false, nil } - var aiWorkQ [1]float64 - impl.Dgetri(n, ait, n, aipiv, aiWorkQ[:], -1) - aiWork := ws.aiWork - if len(aiWork) < int(aiWorkQ[0]) { - aiWork = make([]float64, int(aiWorkQ[0])) - ws.aiWork = aiWork + var inverseQuery [1]float64 + impl.Dgetri(n, ait, n, ipiv, inverseQuery[:], -1) + lwork := int(inverseQuery[0]) + if len(ws.work) < lwork { + ws.work = make([]float64, lwork) } - impl.Dgetri(n, ait, n, aipiv, aiWork, len(aiWork)) + impl.Dgetri(n, ait, n, ipiv, ws.work, lwork) - // Ait*Q (n×n) - aitq := ws.aitq[:n*n] blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, blas64.General{Rows: n, Cols: n, Data: ait, Stride: n}, blas64.General{Rows: n, Cols: n, Data: qWork, Stride: n}, 0, blas64.General{Rows: n, Cols: n, Data: aitq, Stride: n}) - - // G*Ait (n×n) - gait := ws.gait[:n*n] blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, blas64.General{Rows: n, Cols: n, Data: g, Stride: n}, blas64.General{Rows: n, Cols: n, Data: ait, Stride: n}, 0, blas64.General{Rows: n, Cols: n, Data: gait, Stride: n}) - - gaitq := ws.gaitq[:n*n] blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, blas64.General{Rows: n, Cols: n, Data: gait, Stride: n}, blas64.General{Rows: n, Cols: n, Data: qWork, Stride: n}, 0, blas64.General{Rows: n, Cols: n, Data: gaitq, Stride: n}) - nn := 2 * n z := ws.z[:nn*nn] for i := range n { for j := range n { @@ -392,102 +485,140 @@ func Dare(A, B, Q, R *mat.Dense, opts *RiccatiOpts) (*RiccatiResult, error) { } } - // Schur decomposition with sorting: |λ| < 1 to top-left - wr := ws.wr[:nn] - wi := ws.wi[:nn] - vs := ws.vs[:nn*nn] + alphaR := ws.wr[:nn] + alphaI := ws.wi[:nn] + beta := ws.beta[:nn] + vectors := ws.vs[:nn*nn] bwork := ws.bwork[:nn] - - selctg := func(wr, wi float64) bool { - return wr*wr+wi*wi < 1 - } - - var workQuery2 [1]float64 - impl.Dgees(lapack.SchurHess, lapack.SortSelected, selctg, - nn, z, nn, wr, wi, vs, nn, workQuery2[:], -1, bwork) - lwork := int(workQuery2[0]) - work := ws.work - if len(work) < lwork { - work = make([]float64, lwork) - ws.work = work + insideUnitCircle := func(real, imag float64) bool { + return math.Hypot(real, imag) < 1 } - sdim, ok := impl.Dgees(lapack.SchurHess, lapack.SortSelected, selctg, - nn, z, nn, wr, wi, vs, nn, work, lwork, bwork) + var workQuery [1]float64 + impl.Dgees(lapack.SchurHess, lapack.SortSelected, insideUnitCircle, + nn, z, nn, alphaR, alphaI, vectors, nn, workQuery[:], -1, bwork) + lwork = int(workQuery[0]) + if len(ws.work) < lwork { + ws.work = make([]float64, lwork) + } + sdim, ok := impl.Dgees(lapack.SchurHess, lapack.SortSelected, insideUnitCircle, + nn, z, nn, alphaR, alphaI, vectors, nn, ws.work, lwork, bwork) if !ok { - return nil, ErrSchurFailed + return discreteRiccatiSubspace{}, true, ErrSchurFailed } if sdim != n { - return nil, ErrNoStabilizing + return discreteRiccatiSubspace{}, true, ErrNoStabilizing + } + for i := range nn { + beta[i] = 1 } + return discreteRiccatiSubspace{ + vectors: vectors, + alphaR: alphaR, + alphaI: alphaI, + beta: beta, + }, true, nil +} - // Extract U11 = vs[0:n, 0:n], U21 = vs[n:2n, 0:n] - u11 := ws.u11[:n*n] - u21 := ws.u21[:n*n] - copyStrided(u11, n, vs, nn, n, n) - copyBlock(u21, n, 0, 0, vs, nn, n, 0, n, n) +func (problem riccatiProblem) generalizedDiscreteStableSubspace() (discreteRiccatiSubspace, error) { + n, m, ws := problem.n, problem.m, problem.ws + nn := 2 * n + rows := nn + m - ipiv := ws.ipiv[:n] - if !impl.Dgetrf(n, n, u11, n, ipiv) { - return nil, ErrNoStabilizing - } + hLeft := ws.pencilH[:rows*nn] + jLeft := ws.pencilJ[:rows*nn] + input := ws.pencilInput[:rows*m] + clear(hLeft) + clear(jLeft) + clear(input) - anorm := impl.Dlange(lapack.MaxColumnSum, n, n, u11, n, work[:n]) - iwork2 := ws.iwork[:n] - rcnd := impl.Dgecon(lapack.MaxColumnSum, n, u11, n, anorm, work[:4*n], iwork2) + aRaw := problem.A.RawMatrix() + bRaw := problem.B.RawMatrix() + qRaw := problem.Q.RawMatrix() + rRaw := problem.R.RawMatrix() + var sRaw blas64.General + if problem.S != nil { + sRaw = problem.S.RawMatrix() + } - xData := ws.xData[:n*n] for i := range n { for j := range n { - xData[i*n+j] = u21[j*n+i] + hLeft[i*nn+j] = aRaw.Data[i*aRaw.Stride+j] + hLeft[(n+i)*nn+j] = -qRaw.Data[i*qRaw.Stride+j] + jLeft[(n+i)*nn+n+j] = aRaw.Data[j*aRaw.Stride+i] } - } - impl.Dgetrs(blas.Trans, n, n, u11, n, ipiv, xData, n) - - symmetrize(xData, n, n) - X := mat.NewDense(n, n, xData) + hLeft[(n+i)*nn+n+i] = 1 + jLeft[i*nn+i] = 1 - // Closed-loop eigenvalues - eig := ws.eig[:n] - for i := range n { - eig[i] = complex(wr[i], wi[i]) + for j := range m { + input[i*m+j] = bRaw.Data[i*bRaw.Stride+j] + if problem.S != nil { + input[(n+i)*m+j] = -sRaw.Data[i*sRaw.Stride+j] + hLeft[(nn+j)*nn+i] = sRaw.Data[i*sRaw.Stride+j] + } + jLeft[(nn+j)*nn+n+i] = -bRaw.Data[i*bRaw.Stride+j] + } + } + for i := range m { + for j := range m { + input[(nn+i)*m+j] = rRaw.Data[i*rRaw.Stride+j] + } } - // Gain K = (R + B'XB)⁻¹ * (B'XA + S') - btx := ws.btx[:m*n] - blas64.Gemm(blas.Trans, blas.NoTrans, - 1, blas64.General{Rows: n, Cols: m, Data: bRaw.Data, Stride: bRaw.Stride}, - blas64.General{Rows: n, Cols: n, Data: xData, Stride: n}, - 0, blas64.General{Rows: m, Cols: n, Data: btx, Stride: n}) + if m != 0 { + tau := ws.tau[:m] + var qrQuery, applyQuery [1]float64 + impl.Dgeqrf(rows, m, nil, m, nil, qrQuery[:], -1) + impl.Dormqr(blas.Left, blas.Trans, rows, nn, m, nil, m, nil, nil, nn, applyQuery[:], -1) + lwork := max(int(qrQuery[0]), int(applyQuery[0])) + if len(ws.work) < lwork { + ws.work = make([]float64, lwork) + } + impl.Dgeqrf(rows, m, input, m, tau, ws.work, lwork) + impl.Dormqr(blas.Left, blas.Trans, rows, nn, m, input, m, tau, hLeft, nn, ws.work, lwork) + impl.Dormqr(blas.Left, blas.Trans, rows, nn, m, input, m, tau, jLeft, nn, ws.work, lwork) + } - rbar := ws.rbar[:m*m] - copyStrided(rbar, m, rRaw.Data, rRaw.Stride, m, m) - blas64.Gemm(blas.NoTrans, blas.NoTrans, - 1, blas64.General{Rows: m, Cols: n, Data: btx, Stride: n}, - blas64.General{Rows: n, Cols: m, Data: bRaw.Data, Stride: bRaw.Stride}, - 1, blas64.General{Rows: m, Cols: m, Data: rbar, Stride: m}) + h := ws.h[:nn*nn] + j := ws.z[:nn*nn] + copy(h, hLeft[m*nn:]) + copy(j, jLeft[m*nn:]) + + alphaR := ws.wr[:nn] + alphaI := ws.wi[:nn] + beta := ws.beta[:nn] + vectors := ws.vs[:nn*nn] + bwork := ws.bwork[:nn] + insideUnitCircle := func(alphaR, alphaI, beta float64) bool { + return math.Hypot(alphaR, alphaI) < math.Abs(beta) + } - if !impl.Dpotrf(blas.Upper, m, rbar, m) { - return nil, ErrSingularR + var workQuery [1]float64 + impl.Dgges(lapack.SchurNone, lapack.SchurHess, lapack.SortSelected, insideUnitCircle, + nn, h, nn, j, nn, alphaR, alphaI, beta, nil, 1, vectors, nn, workQuery[:], -1, bwork) + lwork := int(workQuery[0]) + if len(ws.work) < lwork { + ws.work = make([]float64, lwork) } - // BtXA = BtX * A (m×n) - kData := ws.kData[:m*n] - blas64.Gemm(blas.NoTrans, blas.NoTrans, - 1, blas64.General{Rows: m, Cols: n, Data: btx, Stride: n}, - blas64.General{Rows: n, Cols: n, Data: aRaw.Data, Stride: aRaw.Stride}, - 0, blas64.General{Rows: m, Cols: n, Data: kData, Stride: n}) - if S != nil { - sRaw := S.RawMatrix() - for j := range m { - row := kData[j*n:] - for i := range n { - row[i] += sRaw.Data[i*sRaw.Stride+j] - } + sdim, ok := impl.Dgges(lapack.SchurNone, lapack.SchurHess, lapack.SortSelected, insideUnitCircle, + nn, h, nn, j, nn, alphaR, alphaI, beta, nil, 1, vectors, nn, ws.work, lwork, bwork) + if !ok { + return discreteRiccatiSubspace{}, ErrSchurFailed + } + if sdim != n { + return discreteRiccatiSubspace{}, ErrNoStabilizing + } + for i := range n { + if beta[i] == 0 { + return discreteRiccatiSubspace{}, ErrNoStabilizing } } - impl.Dpotrs(blas.Upper, m, n, rbar, m, kData, n) - K := mat.NewDense(m, n, kData) - return &RiccatiResult{X: X, K: K, Eig: eig, Rcnd: rcnd}, nil + return discreteRiccatiSubspace{ + vectors: vectors, + alphaR: alphaR, + alphaI: alphaI, + beta: beta, + }, nil } diff --git a/riccati_test.go b/riccati_test.go index 5cb100c..d103fc4 100644 --- a/riccati_test.go +++ b/riccati_test.go @@ -32,6 +32,10 @@ func careResidual(A, B, Q, R, X *mat.Dense) float64 { } func dareResidual(A, B, Q, R, X *mat.Dense) float64 { + return dareResidualWithCrossTerm(A, B, Q, R, nil, X) +} + +func dareResidualWithCrossTerm(A, B, Q, R, S, X *mat.Dense) float64 { n, _ := A.Dims() var atx, atxa, atxb, btx, btxb, rbar, btxa, mid, res mat.Dense @@ -42,6 +46,10 @@ func dareResidual(A, B, Q, R, X *mat.Dense) float64 { btx.Mul(B.T(), X) btxb.Mul(&btx, B) btxa.Mul(&btx, A) + if S != nil { + atxb.Add(&atxb, S) + btxa.Add(&btxa, S.T()) + } rbar.Add(R, &btxb) var luRbar mat.LU diff --git a/state_space_utils.go b/state_space_utils.go index 45f405b..e84c146 100644 --- a/state_space_utils.go +++ b/state_space_utils.go @@ -20,6 +20,19 @@ func NewDescriptor(A, B, C, D, E *mat.Dense, dt float64) (*System, error) { return sys, nil } +func newDescriptorOwned(A, B, C, D, E *mat.Dense, dt float64) (*System, error) { + sys, err := newNoCopy(A, B, C, D, dt) + if err != nil { + return nil, err + } + n, _, _ := sys.Dims() + if err := (descriptorPolicy{E: E}).validate(n); err != nil { + return nil, err + } + sys.E = E + return sys, nil +} + func (sys *System) DescriptorE() *mat.Dense { if sys == nil { return nil diff --git a/tuning_goal.go b/tuning_goal.go index 5938074..7b83a4b 100644 --- a/tuning_goal.go +++ b/tuning_goal.go @@ -3,8 +3,6 @@ package controlsys import ( "fmt" "math" - - "gonum.org/v1/gonum/mat" ) type TuningGoalType int @@ -400,36 +398,38 @@ func frequencyGainRange(sys *System, omega []float64, outputWeight, inputWeight } maxGain := 0.0 minGain := math.Inf(1) - responseData := make([]complex128, resp.P*resp.M) - var outputData, outputProduct, inputData, inputProduct []complex128 + var outputProduct, inputProduct []complex128 weightedRows := resp.P + weightedCols := resp.M if outputResponse != nil { - outputData = make([]complex128, outputResponse.P*outputResponse.M) outputProduct = make([]complex128, outputResponse.P*resp.M) weightedRows = outputResponse.P } if inputResponse != nil { - inputData = make([]complex128, inputResponse.P*inputResponse.M) inputProduct = make([]complex128, weightedRows*inputResponse.M) + weightedCols = inputResponse.M + } + var singularValues *complexSVDWorkspace + if weightedRows > 1 && weightedCols > 1 && (weightedRows != 2 || weightedCols != 2) { + singularValues = newComplexSVDWorkspace(weightedRows, weightedCols) } - var singularValues complexSingularValueWorkspace for k := range omega { - gain := complexResponseAt(resp, k, responseData) + gain := complexResponseAt(resp, k) if outputResponse != nil { - gain, err = multiplyComplexMatricesInto(outputProduct, complexResponseAt(outputResponse, k, outputData), gain) + gain, err = multiplyComplexMatricesInto(outputProduct, complexResponseAt(outputResponse, k), gain) if err != nil { return 0, 0, fmt.Errorf("output weight: %w", err) } } if inputResponse != nil { - gain, err = multiplyComplexMatricesInto(inputProduct, gain, complexResponseAt(inputResponse, k, inputData)) + gain, err = multiplyComplexMatricesInto(inputProduct, gain, complexResponseAt(inputResponse, k)) if err != nil { return 0, 0, fmt.Errorf("input weight: %w", err) } } - sigma, err := singularValues.maximum(gain) - if err != nil { - return 0, 0, err + sigma, ok := singularValues.maximumFromFlat(gain.data, 0, gain.rows, gain.cols) + if !ok { + return 0, 0, fmt.Errorf("maximum singular value: decomposition failed: %w", ErrSingularTransform) } if sigma > maxGain { maxGain = sigma @@ -447,16 +447,10 @@ type complexMatrix struct { data []complex128 } -func complexResponseAt(response *FreqResponseMatrix, frequency int, data []complex128) complexMatrix { - if len(data) != response.P*response.M { - data = make([]complex128, response.P*response.M) - } - for i := range response.P { - for j := range response.M { - data[i*response.M+j] = response.At(frequency, i, j) - } - } - return complexMatrix{rows: response.P, cols: response.M, data: data} +func complexResponseAt(response *FreqResponseMatrix, frequency int) complexMatrix { + blockSize := response.P * response.M + base := frequency * blockSize + return complexMatrix{rows: response.P, cols: response.M, data: response.Data[base : base+blockSize]} } func multiplyComplexMatricesInto(dst []complex128, a, b complexMatrix) (complexMatrix, error) { @@ -480,67 +474,3 @@ func multiplyComplexMatricesInto(dst []complex128, a, b complexMatrix) (complexM } return result, nil } - -type complexSingularValueWorkspace struct { - realForm *mat.Dense - svd mat.SVD -} - -func (w *complexSingularValueWorkspace) maximum(a complexMatrix) (float64, error) { - if a.rows == 0 || a.cols == 0 { - return 0, nil - } - if a.rows == 1 || a.cols == 1 { - norm := 0.0 - for _, value := range a.data { - norm = math.Hypot(norm, math.Hypot(real(value), imag(value))) - } - return norm, nil - } - if a.rows == 2 && a.cols == 2 { - return maximumComplex2x2SingularValue(a), nil - } - rows, cols := 2*a.rows, 2*a.cols - if w.realForm == nil { - w.realForm = mat.NewDense(rows, cols, nil) - } else if r, c := w.realForm.Dims(); r != rows || c != cols { - w.realForm.Reset() - w.realForm.ReuseAs(rows, cols) - } - for i := range a.rows { - for j := range a.cols { - value := a.data[i*a.cols+j] - w.realForm.Set(i, j, real(value)) - w.realForm.Set(i, j+a.cols, -imag(value)) - w.realForm.Set(i+a.rows, j, imag(value)) - w.realForm.Set(i+a.rows, j+a.cols, real(value)) - } - } - if ok := w.svd.Factorize(w.realForm, mat.SVDNone); !ok { - return 0, fmt.Errorf("maximum singular value: decomposition failed: %w", ErrSingularTransform) - } - values := w.svd.Values(nil) - return values[0], nil -} - -func maximumComplex2x2SingularValue(a complexMatrix) float64 { - scale := 0.0 - for _, value := range a.data { - scale = math.Max(scale, math.Hypot(real(value), imag(value))) - } - if scale == 0 { - return 0 - } - a00 := a.data[0] / complex(scale, 0) - a01 := a.data[1] / complex(scale, 0) - a10 := a.data[2] / complex(scale, 0) - a11 := a.data[3] / complex(scale, 0) - frobeniusSquared := complexMagnitudeSquared(a00) + complexMagnitudeSquared(a01) + complexMagnitudeSquared(a10) + complexMagnitudeSquared(a11) - determinantSquared := complexMagnitudeSquared(a00*a11 - a01*a10) - discriminant := math.Max(0, frobeniusSquared*frobeniusSquared-4*determinantSquared) - return scale * math.Sqrt((frobeniusSquared+math.Sqrt(discriminant))/2) -} - -func complexMagnitudeSquared(value complex128) float64 { - return real(value)*real(value) + imag(value)*imag(value) -}