1use nalgebra::DMatrix;
4use num_complex::Complex64;
5use runmat_builtins::{ObjectInstance, Tensor, Value};
6use runmat_macros::runtime_builtin;
7
8use crate::builtins::common::spec::{
9 BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
10 ReductionNaN, ResidencyPolicy, ShapeRequirements,
11};
12use crate::builtins::control::type_resolvers::nyquist_type;
13use crate::{build_runtime_error, BuiltinResult, RuntimeError};
14
15const BUILTIN_NAME: &str = "nyquist";
16const TF_CLASS: &str = "tf";
17const EPS: f64 = 1.0e-12;
18const DEFAULT_FREQUENCY_POINTS: usize = 200;
19
20#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::control::nyquist")]
21pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
22 name: "nyquist",
23 op_kind: GpuOpKind::Custom("control-nyquist-frequency-response"),
24 supported_precisions: &[],
25 broadcast: BroadcastSemantics::None,
26 provider_hooks: &[],
27 constant_strategy: ConstantStrategy::InlineLiteral,
28 residency: ResidencyPolicy::GatherImmediately,
29 nan_mode: ReductionNaN::Include,
30 two_pass_threshold: None,
31 workgroup_size: None,
32 accepts_nan_mode: false,
33 notes: "Nyquist frequency-response evaluation runs on the host from transfer-function metadata. GPU-resident coefficient inputs are gathered by tf before model construction.",
34};
35
36#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::control::nyquist")]
37pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
38 name: "nyquist",
39 shape: ShapeRequirements::Any,
40 constant_strategy: ConstantStrategy::InlineLiteral,
41 elementwise: None,
42 reduction: None,
43 emits_nan: false,
44 notes: "nyquist materialises host response vectors and terminates numeric fusion chains.",
45};
46
47fn nyquist_error(message: impl Into<String>) -> RuntimeError {
48 build_runtime_error(message)
49 .with_builtin(BUILTIN_NAME)
50 .build()
51}
52
53#[runtime_builtin(
54 name = "nyquist",
55 category = "control",
56 summary = "Compute or plot Nyquist frequency responses of SISO transfer-function models.",
57 keywords = "nyquist,frequency response,control system,transfer function,tf",
58 sink = true,
59 suppress_auto_output = true,
60 type_resolver(nyquist_type),
61 builtin_path = "crate::builtins::control::nyquist"
62)]
63async fn nyquist_builtin(system: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
64 if rest.len() > 1 {
65 return Err(nyquist_error(
66 "nyquist: expected nyquist(sys) or nyquist(sys, w)",
67 ));
68 }
69
70 let system = TfSystem::parse(system).await?;
71 let frequencies = FrequencySpec::parse(&system, rest.first()).await?;
72 let response = evaluate_nyquist(&system, frequencies)?;
73
74 if let Some(out_count) = crate::output_count::current_output_count() {
75 if out_count == 0 {
76 render_nyquist_plot(&response).await?;
77 return Ok(Value::OutputList(Vec::new()));
78 }
79 if out_count == 1 {
80 return Ok(Value::OutputList(vec![response.re_value()?]));
81 }
82 if out_count == 2 {
83 return Ok(Value::OutputList(vec![
84 response.re_value()?,
85 response.im_value()?,
86 ]));
87 }
88 return Ok(crate::output_count::output_list_with_padding(
89 out_count,
90 response.outputs()?,
91 ));
92 }
93
94 if crate::output_context::requested_output_count() == Some(0) {
95 render_nyquist_plot(&response).await?;
96 return Ok(Value::OutputList(Vec::new()));
97 }
98
99 response.re_value()
100}
101
102#[derive(Clone, Debug)]
103struct TfSystem {
104 numerator: Vec<Complex64>,
105 denominator: Vec<Complex64>,
106 sample_time: f64,
107 is_real: bool,
108}
109
110impl TfSystem {
111 async fn parse(value: Value) -> BuiltinResult<Self> {
112 let gathered = crate::dispatcher::gather_if_needed_async(&value).await?;
113 let Value::Object(object) = gathered else {
114 return Err(nyquist_error(format!(
115 "nyquist: expected a dynamic system model, got {gathered:?}"
116 )));
117 };
118 if object.class_name != TF_CLASS {
119 return Err(nyquist_error(format!(
120 "nyquist: unsupported model class '{}'; only SISO tf objects are currently supported",
121 object.class_name
122 )));
123 }
124
125 let numerator = coefficients(property(&object, "Numerator")?, "Numerator")?;
126 let denominator = coefficients(property(&object, "Denominator")?, "Denominator")?;
127 let sample_time = scalar_property(property(&object, "Ts")?, "Ts")?;
128 let input_delay = scalar_property(property(&object, "InputDelay")?, "InputDelay")?;
129 let output_delay = scalar_property(property(&object, "OutputDelay")?, "OutputDelay")?;
130
131 if !sample_time.is_finite() || sample_time < 0.0 {
132 return Err(nyquist_error(format!(
133 "nyquist: Ts must be a finite non-negative scalar, got {sample_time}"
134 )));
135 }
136 if !input_delay.is_finite() || input_delay < 0.0 {
137 return Err(nyquist_error(format!(
138 "nyquist: InputDelay must be a finite non-negative scalar, got {input_delay}"
139 )));
140 }
141 if !output_delay.is_finite() || output_delay < 0.0 {
142 return Err(nyquist_error(format!(
143 "nyquist: OutputDelay must be a finite non-negative scalar, got {output_delay}"
144 )));
145 }
146 if input_delay.abs() > EPS || output_delay.abs() > EPS {
147 return Err(nyquist_error(
148 "nyquist: transfer functions with input or output delays are not supported yet",
149 ));
150 }
151
152 let numerator = trim_leading_zeros(numerator);
153 let denominator = trim_leading_zeros(denominator);
154 if denominator.is_empty() {
155 return Err(nyquist_error(
156 "nyquist: denominator coefficients cannot be empty",
157 ));
158 }
159 if denominator[0].norm() <= EPS {
160 return Err(nyquist_error(
161 "nyquist: leading denominator coefficient must be non-zero",
162 ));
163 }
164
165 let is_real = numerator
166 .iter()
167 .chain(&denominator)
168 .all(|value| value.im.abs() <= EPS);
169 Ok(Self {
170 numerator,
171 denominator,
172 sample_time,
173 is_real,
174 })
175 }
176
177 fn is_discrete(&self) -> bool {
178 self.sample_time > 0.0
179 }
180}
181
182fn property<'a>(object: &'a ObjectInstance, name: &str) -> BuiltinResult<&'a Value> {
183 object
184 .properties
185 .get(name)
186 .ok_or_else(|| nyquist_error(format!("nyquist: tf object is missing {name} property")))
187}
188
189fn coefficients(value: &Value, label: &str) -> BuiltinResult<Vec<Complex64>> {
190 match value {
191 Value::Tensor(tensor) => {
192 ensure_vector(label, &tensor.shape)?;
193 finite_complex_values(
194 label,
195 tensor
196 .data
197 .iter()
198 .map(|&value| Complex64::new(value, 0.0))
199 .collect(),
200 )
201 }
202 Value::ComplexTensor(tensor) => {
203 ensure_vector(label, &tensor.shape)?;
204 finite_complex_values(
205 label,
206 tensor
207 .data
208 .iter()
209 .map(|&(re, im)| Complex64::new(re, im))
210 .collect(),
211 )
212 }
213 Value::Num(n) => finite_complex_values(label, vec![Complex64::new(*n, 0.0)]),
214 Value::Int(i) => finite_complex_values(label, vec![Complex64::new(i.to_f64(), 0.0)]),
215 Value::Bool(b) => {
216 finite_complex_values(label, vec![Complex64::new(if *b { 1.0 } else { 0.0 }, 0.0)])
217 }
218 Value::Complex(re, im) => finite_complex_values(label, vec![Complex64::new(*re, *im)]),
219 other => Err(nyquist_error(format!(
220 "nyquist: {label} must be a numeric coefficient vector, got {other:?}"
221 ))),
222 }
223}
224
225fn finite_complex_values(label: &str, values: Vec<Complex64>) -> BuiltinResult<Vec<Complex64>> {
226 if values
227 .iter()
228 .any(|value| !value.re.is_finite() || !value.im.is_finite())
229 {
230 return Err(nyquist_error(format!(
231 "nyquist: {label} coefficients must be finite"
232 )));
233 }
234 Ok(values)
235}
236
237fn ensure_vector(label: &str, shape: &[usize]) -> BuiltinResult<()> {
238 let non_unit = shape.iter().copied().filter(|&dim| dim > 1).count();
239 if non_unit <= 1 {
240 Ok(())
241 } else {
242 Err(nyquist_error(format!(
243 "nyquist: {label} coefficients must be a vector"
244 )))
245 }
246}
247
248fn scalar_property(value: &Value, label: &str) -> BuiltinResult<f64> {
249 match value {
250 Value::Num(n) => Ok(*n),
251 Value::Int(i) => Ok(i.to_f64()),
252 Value::Bool(b) => Ok(if *b { 1.0 } else { 0.0 }),
253 Value::Tensor(tensor) if tensor.data.len() == 1 => Ok(tensor.data[0]),
254 other => Err(nyquist_error(format!(
255 "nyquist: {label} must be a real scalar, got {other:?}"
256 ))),
257 }
258}
259
260#[derive(Clone, Debug)]
261enum FrequencySpec {
262 Values(Vec<f64>),
263}
264
265impl FrequencySpec {
266 async fn parse(system: &TfSystem, value: Option<&Value>) -> BuiltinResult<Self> {
267 let Some(value) = value else {
268 return Ok(Self::Values(default_frequency_vector(system)));
269 };
270 let gathered = crate::dispatcher::gather_if_needed_async(value).await?;
271 let tensor = frequency_tensor_from_value(gathered)?;
272 ensure_vector("frequency", &tensor.shape)?;
273 validate_frequency_vector(&tensor.data)?;
274 Ok(Self::Values(tensor.data))
275 }
276}
277
278fn frequency_tensor_from_value(value: Value) -> BuiltinResult<Tensor> {
279 match value {
280 Value::Tensor(tensor) => Ok(tensor),
281 Value::Num(n) => {
282 Tensor::new(vec![n], vec![1, 1]).map_err(|err| nyquist_error(format!("nyquist: {err}")))
283 }
284 Value::Int(i) => Tensor::new(vec![i.to_f64()], vec![1, 1])
285 .map_err(|err| nyquist_error(format!("nyquist: {err}"))),
286 Value::Bool(b) => Tensor::new(vec![if b { 1.0 } else { 0.0 }], vec![1, 1])
287 .map_err(|err| nyquist_error(format!("nyquist: {err}"))),
288 other => Err(nyquist_error(format!(
289 "nyquist: frequency input must be a real numeric vector, got {other:?}"
290 ))),
291 }
292}
293
294fn validate_frequency_vector(values: &[f64]) -> BuiltinResult<()> {
295 if values.is_empty() {
296 return Err(nyquist_error("nyquist: frequency vector cannot be empty"));
297 }
298 if values
299 .iter()
300 .any(|value| !value.is_finite() || *value < 0.0)
301 {
302 return Err(nyquist_error(
303 "nyquist: frequency values must be finite and non-negative",
304 ));
305 }
306 Ok(())
307}
308
309fn default_frequency_vector(system: &TfSystem) -> Vec<f64> {
310 if system.is_discrete() {
311 return open_linspace(
312 0.0,
313 std::f64::consts::PI / system.sample_time,
314 DEFAULT_FREQUENCY_POINTS,
315 );
316 }
317
318 let mut breakpoints = Vec::new();
319 for coeffs in [&system.numerator, &system.denominator] {
320 if let Ok(roots) = polynomial_roots(coeffs) {
321 breakpoints.extend(
322 roots
323 .into_iter()
324 .map(|root| root.norm())
325 .filter(|value| value.is_finite() && *value > EPS),
326 );
327 }
328 }
329
330 if breakpoints.is_empty() {
331 return logspace(-2.0, 2.0, DEFAULT_FREQUENCY_POINTS);
332 }
333
334 let min_w = breakpoints.iter().copied().fold(f64::INFINITY, f64::min);
335 let max_w = breakpoints.iter().copied().fold(0.0, f64::max);
336 let start = (min_w / 100.0).max(1.0e-4);
337 let stop = (max_w * 100.0).max(start * 10.0);
338 logspace(start.log10(), stop.log10(), DEFAULT_FREQUENCY_POINTS)
339}
340
341#[derive(Clone, Debug)]
342struct NyquistResponse {
343 re: Vec<f64>,
344 im: Vec<f64>,
345 w: Vec<f64>,
346 mirror_negative_frequency: bool,
347}
348
349impl NyquistResponse {
350 fn re_value(&self) -> BuiltinResult<Value> {
351 column_tensor(self.re.clone())
352 }
353
354 fn im_value(&self) -> BuiltinResult<Value> {
355 column_tensor(self.im.clone())
356 }
357
358 fn w_value(&self) -> BuiltinResult<Value> {
359 column_tensor(self.w.clone())
360 }
361
362 fn outputs(&self) -> BuiltinResult<Vec<Value>> {
363 Ok(vec![self.re_value()?, self.im_value()?, self.w_value()?])
364 }
365}
366
367fn evaluate_nyquist(
368 system: &TfSystem,
369 frequencies: FrequencySpec,
370) -> BuiltinResult<NyquistResponse> {
371 let FrequencySpec::Values(w) = frequencies;
372 let mut re = Vec::with_capacity(w.len());
373 let mut im = Vec::with_capacity(w.len());
374 for &frequency in &w {
375 let point = if system.is_discrete() {
376 let phase = frequency * system.sample_time;
377 Complex64::new(phase.cos(), phase.sin())
378 } else {
379 Complex64::new(0.0, frequency)
380 };
381 let value = transfer_response(system, point)?;
382 re.push(zero_small(value.re));
383 im.push(zero_small(value.im));
384 }
385 Ok(NyquistResponse {
386 re,
387 im,
388 w,
389 mirror_negative_frequency: system.is_real,
390 })
391}
392
393fn transfer_response(system: &TfSystem, point: Complex64) -> BuiltinResult<Complex64> {
394 let numerator = polynomial_eval(&system.numerator, point);
395 let denominator = polynomial_eval(&system.denominator, point);
396 if denominator.norm() <= EPS {
397 return Err(nyquist_error(
398 "nyquist: frequency response is singular at one or more requested frequencies",
399 ));
400 }
401 Ok(numerator / denominator)
402}
403
404fn polynomial_eval(coeffs: &[Complex64], point: Complex64) -> Complex64 {
405 let mut acc = Complex64::new(0.0, 0.0);
406 for &coeff in coeffs {
407 acc = acc * point + coeff;
408 }
409 acc
410}
411
412fn polynomial_roots(coeffs: &[Complex64]) -> BuiltinResult<Vec<Complex64>> {
413 let trimmed = trim_leading_zeros(coeffs.to_vec());
414 if trimmed.len() <= 1 {
415 return Ok(Vec::new());
416 }
417 if trimmed.len() == 2 {
418 return Ok(vec![-trimmed[1] / trimmed[0]]);
419 }
420
421 let degree = trimmed.len() - 1;
422 let leading = trimmed[0];
423 let mut companion = DMatrix::<Complex64>::zeros(degree, degree);
424 for row in 1..degree {
425 companion[(row, row - 1)] = Complex64::new(1.0, 0.0);
426 }
427 for (idx, coeff) in trimmed.iter().enumerate().skip(1) {
428 companion[(0, idx - 1)] = -*coeff / leading;
429 }
430 let eigenvalues = companion
431 .eigenvalues()
432 .ok_or_else(|| nyquist_error("nyquist: failed to compute transfer-function roots"))?;
433 Ok(eigenvalues.iter().copied().collect())
434}
435
436fn trim_leading_zeros(values: Vec<Complex64>) -> Vec<Complex64> {
437 let first = values.iter().position(|value| value.norm() > EPS);
438 match first {
439 Some(idx) => values[idx..].to_vec(),
440 None => Vec::new(),
441 }
442}
443
444async fn render_nyquist_plot(response: &NyquistResponse) -> BuiltinResult<()> {
445 let mut args = vec![response.re_value()?, response.im_value()?];
446 if response.mirror_negative_frequency {
447 let negative_im = response.im.iter().map(|value| -*value).collect::<Vec<_>>();
448 args.push(response.re_value()?);
449 args.push(column_tensor(negative_im)?);
450 }
451
452 if let Some((&re0, &im0)) = response.re.first().zip(response.im.first()) {
453 args.push(column_tensor(vec![re0])?);
454 args.push(column_tensor(vec![im0])?);
455 args.push(Value::from("x"));
456 if response.mirror_negative_frequency {
457 args.push(column_tensor(vec![re0])?);
458 args.push(column_tensor(vec![-im0])?);
459 args.push(Value::from("o"));
460 }
461 }
462
463 if let Err(err) = crate::call_builtin_async("plot", &args).await {
464 if super::is_nonfatal_plot_setup_error(&err) {
465 return Ok(());
466 }
467 return Err(err);
468 }
469 let _ = crate::call_builtin_async("title", &[Value::from("Nyquist Diagram")]).await;
470 let _ = crate::call_builtin_async("xlabel", &[Value::from("Real Axis")]).await;
471 let _ = crate::call_builtin_async("ylabel", &[Value::from("Imaginary Axis")]).await;
472 let _ = crate::call_builtin_async("grid", &[Value::from("on")]).await;
473 Ok(())
474}
475
476fn column_tensor(data: Vec<f64>) -> BuiltinResult<Value> {
477 let rows = data.len();
478 let tensor =
479 Tensor::new(data, vec![rows, 1]).map_err(|err| nyquist_error(format!("nyquist: {err}")))?;
480 Ok(Value::Tensor(tensor))
481}
482
483fn linspace(start: f64, stop: f64, count: usize) -> Vec<f64> {
484 if count <= 1 {
485 return vec![start];
486 }
487 let step = (stop - start) / ((count - 1) as f64);
488 (0..count).map(|idx| start + idx as f64 * step).collect()
489}
490
491fn open_linspace(start: f64, stop: f64, count: usize) -> Vec<f64> {
492 if count == 0 {
493 return Vec::new();
494 }
495 let step = (stop - start) / ((count + 1) as f64);
496 (0..count)
497 .map(|idx| start + (idx + 1) as f64 * step)
498 .collect()
499}
500
501fn logspace(start_exp: f64, stop_exp: f64, count: usize) -> Vec<f64> {
502 linspace(start_exp, stop_exp, count)
503 .into_iter()
504 .map(|value| 10.0_f64.powf(value))
505 .collect()
506}
507
508fn zero_small(value: f64) -> f64 {
509 if value.abs() <= EPS {
510 0.0
511 } else {
512 value
513 }
514}
515
516#[cfg(test)]
517mod tests {
518 use super::*;
519 use futures::executor::block_on;
520 use runmat_builtins::{CharArray, ComplexTensor};
521
522 fn tf_object(num: Vec<f64>, den: Vec<f64>, ts: f64) -> Value {
523 let mut object = ObjectInstance::new("tf".to_string());
524 object.properties.insert(
525 "Numerator".to_string(),
526 Value::Tensor(Tensor::new(num.clone(), vec![1, num.len()]).unwrap()),
527 );
528 object.properties.insert(
529 "Denominator".to_string(),
530 Value::Tensor(Tensor::new(den.clone(), vec![1, den.len()]).unwrap()),
531 );
532 object.properties.insert(
533 "Variable".to_string(),
534 Value::CharArray(CharArray::new_row(if ts > 0.0 { "z" } else { "s" })),
535 );
536 object.properties.insert("Ts".to_string(), Value::Num(ts));
537 object
538 .properties
539 .insert("InputDelay".to_string(), Value::Num(0.0));
540 object
541 .properties
542 .insert("OutputDelay".to_string(), Value::Num(0.0));
543 Value::Object(object)
544 }
545
546 fn run_nyquist(system: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
547 block_on(nyquist_builtin(system, rest))
548 }
549
550 fn tensor_data(value: Value) -> Vec<f64> {
551 match value {
552 Value::Tensor(tensor) => tensor.data,
553 other => panic!("expected tensor, got {other:?}"),
554 }
555 }
556
557 #[test]
558 fn nyquist_first_order_continuous_explicit_frequency() {
559 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
560 let w = Value::Tensor(Tensor::new(vec![0.0, 1.0, 2.0], vec![1, 3]).unwrap());
561 let _guard = crate::output_count::push_output_count(Some(3));
562 let result = run_nyquist(sys, vec![w]).expect("nyquist");
563 let Value::OutputList(outputs) = result else {
564 panic!("expected output list");
565 };
566 let re = tensor_data(outputs[0].clone());
567 let im = tensor_data(outputs[1].clone());
568 let w_out = tensor_data(outputs[2].clone());
569 let expected_re = [1.0, 0.5, 0.2];
570 let expected_im = [0.0, -0.5, -0.4];
571 for ((actual_re, actual_im), (expected_re, expected_im)) in re
572 .iter()
573 .zip(&im)
574 .zip(expected_re.into_iter().zip(expected_im))
575 {
576 assert!((actual_re - expected_re).abs() < 1.0e-12);
577 assert!((actual_im - expected_im).abs() < 1.0e-12);
578 }
579 assert_eq!(w_out, vec![0.0, 1.0, 2.0]);
580 }
581
582 #[test]
583 fn nyquist_two_outputs_returns_real_and_imaginary_columns() {
584 let sys = tf_object(vec![1.0], vec![1.0, 2.0, 1.0], 0.0);
585 let w = Value::Tensor(Tensor::new(vec![0.0, 1.0], vec![2, 1]).unwrap());
586 let _guard = crate::output_count::push_output_count(Some(2));
587 let result = run_nyquist(sys, vec![w]).expect("nyquist");
588 let Value::OutputList(outputs) = result else {
589 panic!("expected output list");
590 };
591 assert_eq!(outputs.len(), 2);
592 match &outputs[0] {
593 Value::Tensor(tensor) => assert_eq!(tensor.shape, vec![2, 1]),
594 other => panic!("expected real tensor, got {other:?}"),
595 }
596 let re = tensor_data(outputs[0].clone());
597 let im = tensor_data(outputs[1].clone());
598 assert!((re[0] - 1.0).abs() < 1.0e-12);
599 assert!(im[0].abs() < 1.0e-12);
600 assert!(re[1].abs() < 1.0e-12);
601 assert!((im[1] + 0.5).abs() < 1.0e-12);
602 }
603
604 #[test]
605 fn nyquist_discrete_uses_unit_circle_frequency_mapping() {
606 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 0.1);
607 let w = Value::Tensor(Tensor::new(vec![0.0], vec![1, 1]).unwrap());
608 let _guard = crate::output_count::push_output_count(Some(2));
609 let result = run_nyquist(sys, vec![w]).expect("nyquist");
610 let Value::OutputList(outputs) = result else {
611 panic!("expected output list");
612 };
613 let re = tensor_data(outputs[0].clone());
614 let im = tensor_data(outputs[1].clone());
615 assert!((re[0] - 2.0).abs() < 1.0e-12);
616 assert!(im[0].abs() < 1.0e-12);
617 }
618
619 #[test]
620 fn nyquist_discrete_default_grid_excludes_singular_unit_circle_endpoints() {
621 let system = TfSystem {
622 numerator: vec![Complex64::new(1.0, 0.0)],
623 denominator: vec![Complex64::new(1.0, 0.0), Complex64::new(-1.0, 0.0)],
624 sample_time: 0.1,
625 is_real: true,
626 };
627 let w = default_frequency_vector(&system);
628 assert_eq!(w.len(), DEFAULT_FREQUENCY_POINTS);
629 assert!(w[0] > 0.0);
630 assert!(w[w.len() - 1] < std::f64::consts::PI / system.sample_time);
631
632 let _guard = crate::output_count::push_output_count(Some(3));
633 run_nyquist(tf_object(vec![1.0], vec![1.0, -1.0], 0.1), Vec::new())
634 .expect("pole at z=1 should not be evaluated at w=0");
635 run_nyquist(tf_object(vec![1.0], vec![1.0, 1.0], 0.1), Vec::new())
636 .expect("pole at z=-1 should not be evaluated at the Nyquist frequency");
637 }
638
639 #[test]
640 fn nyquist_statement_form_plots_without_error() {
641 let sys = tf_object(vec![1.0], vec![1.0, 2.0, 1.0], 0.0);
642 let _guard = crate::output_count::push_output_count(Some(0));
643 let result = run_nyquist(sys, Vec::new()).expect("nyquist");
644 assert!(matches!(result, Value::OutputList(outputs) if outputs.is_empty()));
645 }
646
647 #[test]
648 fn nyquist_rejects_invalid_frequency_vector() {
649 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
650 let w = Value::Tensor(Tensor::new(vec![0.0, f64::INFINITY], vec![1, 2]).unwrap());
651 let err = run_nyquist(sys, vec![w]).expect_err("should fail");
652 assert!(err.message().contains("frequency values must be finite"));
653 }
654
655 #[test]
656 fn nyquist_rejects_unsupported_model_type() {
657 let object = ObjectInstance::new("ss".to_string());
658 let err = run_nyquist(Value::Object(object), Vec::new()).expect_err("should fail");
659 assert!(err.message().contains("unsupported model class"));
660 }
661
662 #[test]
663 fn nyquist_complex_coefficients_are_supported() {
664 let mut object = ObjectInstance::new("tf".to_string());
665 object.properties.insert(
666 "Numerator".to_string(),
667 Value::ComplexTensor(
668 ComplexTensor::new(vec![(1.0, 1.0)], vec![1, 1]).expect("numerator"),
669 ),
670 );
671 object.properties.insert(
672 "Denominator".to_string(),
673 Value::Tensor(Tensor::new(vec![1.0, 1.0], vec![1, 2]).unwrap()),
674 );
675 object.properties.insert(
676 "Variable".to_string(),
677 Value::CharArray(CharArray::new_row("s")),
678 );
679 object.properties.insert("Ts".to_string(), Value::Num(0.0));
680 object
681 .properties
682 .insert("InputDelay".to_string(), Value::Num(0.0));
683 object
684 .properties
685 .insert("OutputDelay".to_string(), Value::Num(0.0));
686
687 let w = Value::Tensor(Tensor::new(vec![0.0], vec![1, 1]).unwrap());
688 let _guard = crate::output_count::push_output_count(Some(2));
689 let result = run_nyquist(Value::Object(object), vec![w]).expect("nyquist");
690 let Value::OutputList(outputs) = result else {
691 panic!("expected output list");
692 };
693 let re = tensor_data(outputs[0].clone());
694 let im = tensor_data(outputs[1].clone());
695 assert!((re[0] - 1.0).abs() < 1.0e-12);
696 assert!((im[0] - 1.0).abs() < 1.0e-12);
697 }
698}