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