1use nalgebra::DMatrix;
4use runmat_builtins::{
5 BuiltinCompletionPolicy, BuiltinDescriptor, BuiltinErrorDescriptor, BuiltinOutputMode,
6 BuiltinParamArity, BuiltinParamDescriptor, BuiltinParamType, BuiltinSignatureDescriptor,
7 Tensor, Value,
8};
9use runmat_macros::runtime_builtin;
10
11use crate::builtins::common::spec::{
12 BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
13 ReductionNaN, ResidencyPolicy, ShapeRequirements,
14};
15use crate::builtins::common::tensor;
16use crate::builtins::control::type_resolvers::impulse_type;
17use crate::{build_runtime_error, BuiltinResult, RuntimeError};
18
19const BUILTIN_NAME: &str = "impulse";
20const TF_CLASS: &str = "tf";
21const EPS: f64 = 1.0e-12;
22const DEFAULT_POINTS: usize = 100;
23const MAX_DISCRETE_SAMPLES: usize = 1_000_000;
24
25const IMPULSE_OUTPUT_Y: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
26 name: "y",
27 ty: BuiltinParamType::NumericArray,
28 arity: BuiltinParamArity::Required,
29 default: None,
30 description: "Impulse response samples (column vector).",
31}];
32const IMPULSE_OUTPUT_Y_T: [BuiltinParamDescriptor; 2] = [
33 BuiltinParamDescriptor {
34 name: "y",
35 ty: BuiltinParamType::NumericArray,
36 arity: BuiltinParamArity::Required,
37 default: None,
38 description: "Impulse response samples (column vector).",
39 },
40 BuiltinParamDescriptor {
41 name: "t",
42 ty: BuiltinParamType::NumericArray,
43 arity: BuiltinParamArity::Required,
44 default: None,
45 description: "Time samples (column vector).",
46 },
47];
48const IMPULSE_INPUTS_SYS: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
49 name: "sys",
50 ty: BuiltinParamType::Any,
51 arity: BuiltinParamArity::Required,
52 default: None,
53 description: "SISO tf model.",
54}];
55const IMPULSE_INPUTS_SYS_TIME: [BuiltinParamDescriptor; 2] = [
56 BuiltinParamDescriptor {
57 name: "sys",
58 ty: BuiltinParamType::Any,
59 arity: BuiltinParamArity::Required,
60 default: None,
61 description: "SISO tf model.",
62 },
63 BuiltinParamDescriptor {
64 name: "time",
65 ty: BuiltinParamType::Any,
66 arity: BuiltinParamArity::Optional,
67 default: None,
68 description: "Final time scalar or explicit time vector.",
69 },
70];
71const IMPULSE_SIGNATURES: [BuiltinSignatureDescriptor; 6] = [
72 BuiltinSignatureDescriptor {
73 label: "y = impulse(sys)",
74 inputs: &IMPULSE_INPUTS_SYS,
75 outputs: &IMPULSE_OUTPUT_Y,
76 },
77 BuiltinSignatureDescriptor {
78 label: "y = impulse(sys, tFinal)",
79 inputs: &IMPULSE_INPUTS_SYS_TIME,
80 outputs: &IMPULSE_OUTPUT_Y,
81 },
82 BuiltinSignatureDescriptor {
83 label: "y = impulse(sys, t)",
84 inputs: &IMPULSE_INPUTS_SYS_TIME,
85 outputs: &IMPULSE_OUTPUT_Y,
86 },
87 BuiltinSignatureDescriptor {
88 label: "[y,t] = impulse(sys)",
89 inputs: &IMPULSE_INPUTS_SYS,
90 outputs: &IMPULSE_OUTPUT_Y_T,
91 },
92 BuiltinSignatureDescriptor {
93 label: "[y,t] = impulse(sys, tFinal)",
94 inputs: &IMPULSE_INPUTS_SYS_TIME,
95 outputs: &IMPULSE_OUTPUT_Y_T,
96 },
97 BuiltinSignatureDescriptor {
98 label: "[y,t] = impulse(sys, t)",
99 inputs: &IMPULSE_INPUTS_SYS_TIME,
100 outputs: &IMPULSE_OUTPUT_Y_T,
101 },
102];
103const IMPULSE_ERROR_INVALID_ARGUMENT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
104 code: "RM.IMPULSE.INVALID_ARGUMENT",
105 identifier: Some("RunMat:impulse:InvalidArgument"),
106 when: "Inputs do not match supported impulse invocation forms.",
107 message: "impulse: invalid argument",
108};
109const IMPULSE_ERROR_INVALID_MODEL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
110 code: "RM.IMPULSE.INVALID_MODEL",
111 identifier: Some("RunMat:impulse:InvalidModel"),
112 when: "Input system is not a supported tf object with valid required metadata.",
113 message: "impulse: invalid model",
114};
115const IMPULSE_ERROR_INVALID_TIME: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
116 code: "RM.IMPULSE.INVALID_TIME",
117 identifier: Some("RunMat:impulse:InvalidTime"),
118 when: "Time argument is invalid for the model class or sampling mode.",
119 message: "impulse: invalid time input",
120};
121const IMPULSE_ERROR_UNSUPPORTED_MODEL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
122 code: "RM.IMPULSE.UNSUPPORTED_MODEL",
123 identifier: Some("RunMat:impulse:UnsupportedModel"),
124 when: "Model is well-formed but unsupported by the current impulse implementation.",
125 message: "impulse: unsupported model",
126};
127const IMPULSE_ERROR_DISCRETE_LIMIT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
128 code: "RM.IMPULSE.DISCRETE_LIMIT",
129 identifier: Some("RunMat:impulse:DiscreteLimit"),
130 when: "Discrete simulation would exceed platform or configured sample limits.",
131 message: "impulse: discrete simulation limit exceeded",
132};
133const IMPULSE_ERROR_INTERNAL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
134 code: "RM.IMPULSE.INTERNAL",
135 identifier: Some("RunMat:impulse:Internal"),
136 when: "Internal response assembly failed.",
137 message: "impulse: internal error",
138};
139const IMPULSE_ERRORS: [BuiltinErrorDescriptor; 6] = [
140 IMPULSE_ERROR_INVALID_ARGUMENT,
141 IMPULSE_ERROR_INVALID_MODEL,
142 IMPULSE_ERROR_INVALID_TIME,
143 IMPULSE_ERROR_UNSUPPORTED_MODEL,
144 IMPULSE_ERROR_DISCRETE_LIMIT,
145 IMPULSE_ERROR_INTERNAL,
146];
147pub const IMPULSE_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
148 signatures: &IMPULSE_SIGNATURES,
149 output_mode: BuiltinOutputMode::ByRequestedOutputCount,
150 completion_policy: BuiltinCompletionPolicy::Public,
151 errors: &IMPULSE_ERRORS,
152};
153
154#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::control::impulse")]
155pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
156 name: "impulse",
157 op_kind: GpuOpKind::Custom("control-impulse-response"),
158 supported_precisions: &[],
159 broadcast: BroadcastSemantics::None,
160 provider_hooks: &[],
161 constant_strategy: ConstantStrategy::InlineLiteral,
162 residency: ResidencyPolicy::GatherImmediately,
163 nan_mode: ReductionNaN::Include,
164 two_pass_threshold: None,
165 workgroup_size: None,
166 accepts_nan_mode: false,
167 notes: "Control-system response evaluation runs on the host. GPU-resident metadata is gathered before simulation.",
168};
169
170#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::control::impulse")]
171pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
172 name: "impulse",
173 shape: ShapeRequirements::Any,
174 constant_strategy: ConstantStrategy::InlineLiteral,
175 elementwise: None,
176 reduction: None,
177 emits_nan: false,
178 notes: "Impulse-response simulation materialises host-side time and output vectors and terminates fusion chains.",
179};
180
181fn impulse_error_with_detail(
182 error: &'static BuiltinErrorDescriptor,
183 detail: impl AsRef<str>,
184) -> RuntimeError {
185 impulse_error_with_message(format!("{}: {}", error.message, detail.as_ref()), error)
186}
187
188fn impulse_error_with_message(
189 message: impl Into<String>,
190 error: &'static BuiltinErrorDescriptor,
191) -> RuntimeError {
192 let mut builder = build_runtime_error(message).with_builtin(BUILTIN_NAME);
193 if let Some(identifier) = error.identifier {
194 builder = builder.with_identifier(identifier);
195 }
196 builder.build()
197}
198
199#[runtime_builtin(
200 name = "impulse",
201 category = "control",
202 summary = "Compute or plot impulse responses.",
203 keywords = "impulse,control system,transfer function,response,tf",
204 type_resolver(impulse_type),
205 descriptor(crate::builtins::control::impulse::IMPULSE_DESCRIPTOR),
206 builtin_path = "crate::builtins::control::impulse"
207)]
208async fn impulse_builtin(system: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
209 let system = TfSystem::parse(system).await?;
210 let time = TimeSpec::parse(&system, &rest).await?;
211 let response = evaluate_impulse(&system, time)?;
212
213 if let Some(out_count) = crate::output_count::current_output_count() {
214 if out_count == 0 {
215 emit_impulse_plot(&response).await?;
216 return Ok(Value::OutputList(Vec::new()));
217 }
218 if out_count == 1 {
219 return Ok(Value::OutputList(vec![response.y_value()?]));
220 }
221 return Ok(crate::output_count::output_list_with_padding(
222 out_count,
223 vec![response.y_value()?, response.t_value()?],
224 ));
225 }
226
227 if crate::output_context::requested_output_count() == Some(0) {
228 emit_impulse_plot(&response).await?;
229 return Ok(Value::OutputList(Vec::new()));
230 }
231
232 response.y_value()
233}
234
235async fn emit_impulse_plot(response: &ImpulseResponse) -> BuiltinResult<()> {
236 if let Err(err) = render_impulse_plot(response).await {
237 if is_nonfatal_plot_setup_error(&err) {
238 return Ok(());
239 }
240 return Err(err);
241 }
242 Ok(())
243}
244
245fn is_nonfatal_plot_setup_error(err: &RuntimeError) -> bool {
246 let lower = err.message().to_ascii_lowercase();
247 lower.contains("plotting is unavailable")
248 || lower.contains("non-main thread")
249 || lower.contains("interactive plotting failed")
250}
251
252#[derive(Clone, Debug)]
253struct TfSystem {
254 numerator: Vec<f64>,
255 denominator: Vec<f64>,
256 sample_time: f64,
257}
258
259impl TfSystem {
260 async fn parse(value: Value) -> BuiltinResult<Self> {
261 let gathered = crate::dispatcher::gather_if_needed_async(&value).await?;
262 let Value::Object(object) = gathered else {
263 return Err(impulse_error_with_detail(
264 &IMPULSE_ERROR_INVALID_MODEL,
265 format!("expected a dynamic system model, got {gathered:?}"),
266 ));
267 };
268 if object.class_name != TF_CLASS {
269 return Err(impulse_error_with_detail(
270 &IMPULSE_ERROR_UNSUPPORTED_MODEL,
271 format!(
272 "unsupported model class '{}'; only SISO tf objects are currently supported",
273 object.class_name
274 ),
275 ));
276 }
277
278 let numerator = real_coefficients(property(&object, "Numerator")?, "Numerator")?;
279 let denominator = real_coefficients(property(&object, "Denominator")?, "Denominator")?;
280 let sample_time = scalar_property(property(&object, "Ts")?, "Ts")?;
281 let input_delay = scalar_property(property(&object, "InputDelay")?, "InputDelay")?;
282 let output_delay = scalar_property(property(&object, "OutputDelay")?, "OutputDelay")?;
283 if !sample_time.is_finite() || sample_time < 0.0 {
284 return Err(impulse_error_with_detail(
285 &IMPULSE_ERROR_INVALID_MODEL,
286 format!("Ts must be a finite non-negative scalar, got {sample_time}"),
287 ));
288 }
289 if !input_delay.is_finite() || input_delay < 0.0 {
290 return Err(impulse_error_with_detail(
291 &IMPULSE_ERROR_INVALID_MODEL,
292 format!("InputDelay must be a finite non-negative scalar, got {input_delay}"),
293 ));
294 }
295 if !output_delay.is_finite() || output_delay < 0.0 {
296 return Err(impulse_error_with_detail(
297 &IMPULSE_ERROR_INVALID_MODEL,
298 format!("OutputDelay must be a finite non-negative scalar, got {output_delay}"),
299 ));
300 }
301 if input_delay.abs() > EPS || output_delay.abs() > EPS {
302 return Err(impulse_error_with_detail(
303 &IMPULSE_ERROR_UNSUPPORTED_MODEL,
304 "transfer functions with input or output delays are not supported yet",
305 ));
306 }
307
308 let numerator = trim_leading_zeros(numerator);
309 let denominator = trim_leading_zeros(denominator);
310 if denominator.is_empty() {
311 return Err(impulse_error_with_detail(
312 &IMPULSE_ERROR_INVALID_MODEL,
313 "denominator coefficients cannot be empty",
314 ));
315 }
316 if numerator.is_empty() {
317 return Ok(Self {
318 numerator,
319 denominator,
320 sample_time,
321 });
322 }
323 if denominator.len() <= 1 {
324 return Err(impulse_error_with_detail(
325 &IMPULSE_ERROR_UNSUPPORTED_MODEL,
326 "static-gain transfer functions do not have a finite impulse-response vector",
327 ));
328 }
329 if numerator.len() >= denominator.len() {
330 return Err(impulse_error_with_detail(
331 &IMPULSE_ERROR_UNSUPPORTED_MODEL,
332 "only strictly proper SISO tf models are currently supported",
333 ));
334 }
335
336 Ok(Self {
337 numerator,
338 denominator,
339 sample_time,
340 })
341 }
342
343 fn is_discrete(&self) -> bool {
344 self.sample_time > 0.0
345 }
346}
347
348fn property<'a>(
349 object: &'a runmat_builtins::ObjectInstance,
350 name: &str,
351) -> BuiltinResult<&'a Value> {
352 object.properties.get(name).ok_or_else(|| {
353 impulse_error_with_detail(
354 &IMPULSE_ERROR_INVALID_MODEL,
355 format!("tf object is missing {name} property"),
356 )
357 })
358}
359
360fn real_coefficients(value: &Value, label: &str) -> BuiltinResult<Vec<f64>> {
361 match value {
362 Value::Tensor(tensor) => {
363 ensure_vector(label, &tensor.shape)?;
364 finite_values(label, tensor.data.clone())
365 }
366 Value::Num(n) => finite_values(label, vec![*n]),
367 Value::Int(i) => finite_values(label, vec![i.to_f64()]),
368 Value::Bool(b) => finite_values(label, vec![if *b { 1.0 } else { 0.0 }]),
369 Value::LogicalArray(logical) => {
370 ensure_vector(label, &logical.shape)?;
371 finite_values(
372 label,
373 logical
374 .data
375 .iter()
376 .map(|bit| if *bit == 0 { 0.0 } else { 1.0 })
377 .collect(),
378 )
379 }
380 Value::Complex(_, _) | Value::ComplexTensor(_) => Err(impulse_error_with_detail(
381 &IMPULSE_ERROR_UNSUPPORTED_MODEL,
382 "complex-coefficient tf models are not supported yet",
383 )),
384 other => Err(impulse_error_with_detail(
385 &IMPULSE_ERROR_INVALID_MODEL,
386 format!("{label} must be a real numeric coefficient vector, got {other:?}"),
387 )),
388 }
389}
390
391fn finite_values(label: &str, values: Vec<f64>) -> BuiltinResult<Vec<f64>> {
392 if values.iter().any(|value| !value.is_finite()) {
393 return Err(impulse_error_with_detail(
394 &IMPULSE_ERROR_INVALID_MODEL,
395 format!("{label} coefficients must be finite"),
396 ));
397 }
398 Ok(values)
399}
400
401fn ensure_vector(label: &str, shape: &[usize]) -> BuiltinResult<()> {
402 let non_unit = shape.iter().copied().filter(|&dim| dim > 1).count();
403 if non_unit <= 1 {
404 Ok(())
405 } else {
406 Err(impulse_error_with_detail(
407 &IMPULSE_ERROR_INVALID_MODEL,
408 format!("{label} coefficients must be a vector"),
409 ))
410 }
411}
412
413fn scalar_property(value: &Value, label: &str) -> BuiltinResult<f64> {
414 match value {
415 Value::Num(n) => Ok(*n),
416 Value::Int(i) => Ok(i.to_f64()),
417 Value::Bool(b) => Ok(if *b { 1.0 } else { 0.0 }),
418 Value::Tensor(tensor) if tensor.data.len() == 1 => Ok(tensor.data[0]),
419 other => Err(impulse_error_with_detail(
420 &IMPULSE_ERROR_INVALID_MODEL,
421 format!("{label} must be a real scalar, got {other:?}"),
422 )),
423 }
424}
425
426fn trim_leading_zeros(values: Vec<f64>) -> Vec<f64> {
427 let first = values.iter().position(|value| value.abs() > EPS);
428 match first {
429 Some(idx) => values[idx..].to_vec(),
430 None => Vec::new(),
431 }
432}
433
434#[derive(Clone, Debug)]
435enum TimeSpec {
436 Values(Vec<f64>),
437}
438
439impl TimeSpec {
440 async fn parse(system: &TfSystem, rest: &[Value]) -> BuiltinResult<Self> {
441 match rest {
442 [] => Ok(Self::Values(default_time_vector(system))),
443 [value] => {
444 let gathered = crate::dispatcher::gather_if_needed_async(value).await?;
445 if let Some(final_time) = scalar_time_from_value(&gathered)? {
446 return Ok(Self::Values(time_vector_from_final_time(
447 system, final_time,
448 )?));
449 }
450 let vector = time_vector_from_value(gathered)?;
451 validate_time_vector(system, &vector)?;
452 Ok(Self::Values(vector))
453 }
454 _ => Err(impulse_error_with_detail(
455 &IMPULSE_ERROR_INVALID_ARGUMENT,
456 "expected impulse(sys), impulse(sys, tFinal), or impulse(sys, t)",
457 )),
458 }
459 }
460}
461
462fn scalar_time_from_value(value: &Value) -> BuiltinResult<Option<f64>> {
463 match value {
464 Value::Num(n) => Ok(Some(*n)),
465 Value::Int(i) => Ok(Some(i.to_f64())),
466 Value::Bool(b) => Ok(Some(if *b { 1.0 } else { 0.0 })),
467 Value::Tensor(tensor) if tensor.data.len() == 1 => Ok(Some(tensor.data[0])),
468 Value::LogicalArray(logical) if logical.data.len() == 1 => {
469 Ok(Some(if logical.data[0] == 0 { 0.0 } else { 1.0 }))
470 }
471 Value::Tensor(_) | Value::LogicalArray(_) => Ok(None),
472 _ => Ok(None),
473 }
474}
475
476fn default_time_vector(system: &TfSystem) -> Vec<f64> {
477 if system.is_discrete() {
478 (0..DEFAULT_POINTS)
479 .map(|idx| idx as f64 * system.sample_time)
480 .collect()
481 } else {
482 linspace(0.0, 10.0, DEFAULT_POINTS)
483 }
484}
485
486fn time_vector_from_final_time(system: &TfSystem, final_time: f64) -> BuiltinResult<Vec<f64>> {
487 if !final_time.is_finite() || final_time < 0.0 {
488 return Err(impulse_error_with_detail(
489 &IMPULSE_ERROR_INVALID_TIME,
490 "final time must be a finite non-negative scalar",
491 ));
492 }
493 if system.is_discrete() {
494 let count = checked_discrete_sample_count(system, final_time)?;
495 Ok((0..count)
496 .map(|idx| idx as f64 * system.sample_time)
497 .collect())
498 } else if final_time == 0.0 {
499 Ok(vec![0.0])
500 } else {
501 Ok(linspace(0.0, final_time, DEFAULT_POINTS))
502 }
503}
504
505fn checked_discrete_sample_count(system: &TfSystem, final_time: f64) -> BuiltinResult<usize> {
506 let samples = final_time / system.sample_time;
507 if !samples.is_finite() {
508 return Err(impulse_error_with_detail(
509 &IMPULSE_ERROR_DISCRETE_LIMIT,
510 "discrete sample count exceeds platform limits",
511 ));
512 }
513
514 let count = samples.floor() + 1.0;
515 if count > usize::MAX as f64 || count > MAX_DISCRETE_SAMPLES as f64 {
516 return Err(impulse_error_with_detail(
517 &IMPULSE_ERROR_DISCRETE_LIMIT,
518 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
519 ));
520 }
521 Ok(count as usize)
522}
523
524fn checked_discrete_sample_index(system: &TfSystem, time: f64) -> BuiltinResult<usize> {
525 let samples = time / system.sample_time;
526 let index = samples.round();
527 if !index.is_finite() || index > usize::MAX as f64 {
528 return Err(impulse_error_with_detail(
529 &IMPULSE_ERROR_DISCRETE_LIMIT,
530 "discrete sample index exceeds platform limits",
531 ));
532 }
533 if index >= MAX_DISCRETE_SAMPLES as f64 {
534 return Err(impulse_error_with_detail(
535 &IMPULSE_ERROR_DISCRETE_LIMIT,
536 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
537 ));
538 }
539 Ok(index as usize)
540}
541
542fn linspace(start: f64, stop: f64, count: usize) -> Vec<f64> {
543 if count <= 1 {
544 return vec![start];
545 }
546 let step = (stop - start) / ((count - 1) as f64);
547 (0..count).map(|idx| start + idx as f64 * step).collect()
548}
549
550fn time_vector_from_value(value: Value) -> BuiltinResult<Vec<f64>> {
551 let tensor = tensor::value_into_tensor_for(BUILTIN_NAME, value).map_err(|err| {
552 impulse_error_with_detail(
553 &IMPULSE_ERROR_INVALID_TIME,
554 format!("time vector must be numeric: {err}"),
555 )
556 })?;
557 ensure_vector("time", &tensor.shape)?;
558 if tensor.data.is_empty() {
559 return Err(impulse_error_with_detail(
560 &IMPULSE_ERROR_INVALID_TIME,
561 "time vector cannot be empty",
562 ));
563 }
564 Ok(tensor.data)
565}
566
567fn validate_time_vector(system: &TfSystem, values: &[f64]) -> BuiltinResult<()> {
568 if values
569 .iter()
570 .any(|value| !value.is_finite() || *value < 0.0)
571 {
572 return Err(impulse_error_with_detail(
573 &IMPULSE_ERROR_INVALID_TIME,
574 "time vector values must be finite and non-negative",
575 ));
576 }
577 if values.windows(2).any(|pair| pair[1] <= pair[0]) {
578 return Err(impulse_error_with_detail(
579 &IMPULSE_ERROR_INVALID_TIME,
580 "time vector values must be strictly increasing",
581 ));
582 }
583 if system.is_discrete() {
584 for &value in values {
585 let samples = value / system.sample_time;
586 if (samples - samples.round()).abs() > 1.0e-8 {
587 return Err(impulse_error_with_detail(
588 &IMPULSE_ERROR_INVALID_TIME,
589 "discrete-time vectors must use integer multiples of the sample time",
590 ));
591 }
592 }
593 }
594 Ok(())
595}
596
597#[derive(Clone, Debug)]
598struct ImpulseResponse {
599 t: Vec<f64>,
600 y: Vec<f64>,
601 #[cfg(feature = "plot-core")]
602 discrete: bool,
603}
604
605impl ImpulseResponse {
606 fn y_value(&self) -> BuiltinResult<Value> {
607 let tensor = Tensor::new(self.y.clone(), vec![self.y.len(), 1]).map_err(|err| {
608 impulse_error_with_detail(
609 &IMPULSE_ERROR_INTERNAL,
610 format!("failed to build response tensor: {err}"),
611 )
612 })?;
613 Ok(Value::Tensor(tensor))
614 }
615
616 fn t_value(&self) -> BuiltinResult<Value> {
617 let tensor = Tensor::new(self.t.clone(), vec![self.t.len(), 1]).map_err(|err| {
618 impulse_error_with_detail(
619 &IMPULSE_ERROR_INTERNAL,
620 format!("failed to build time tensor: {err}"),
621 )
622 })?;
623 Ok(Value::Tensor(tensor))
624 }
625}
626
627fn evaluate_impulse(system: &TfSystem, time: TimeSpec) -> BuiltinResult<ImpulseResponse> {
628 let TimeSpec::Values(t) = time;
629 let realization = Realization::from_tf(system)?;
630 let y = if system.is_discrete() {
631 discrete_response(system, &realization, &t)?
632 } else {
633 continuous_response(&realization, &t)
634 };
635 Ok(ImpulseResponse {
636 t,
637 y,
638 #[cfg(feature = "plot-core")]
639 discrete: system.is_discrete(),
640 })
641}
642
643#[derive(Clone, Debug)]
644struct Realization {
645 a: DMatrix<f64>,
646 c: Vec<f64>,
647}
648
649impl Realization {
650 fn from_tf(system: &TfSystem) -> BuiltinResult<Self> {
651 if system.numerator.is_empty() {
652 let order = system.denominator.len().saturating_sub(1).max(1);
653 return Ok(Self {
654 a: DMatrix::zeros(order, order),
655 c: vec![0.0; order],
656 });
657 }
658 let leading = system.denominator[0];
659 if leading.abs() <= EPS {
660 return Err(impulse_error_with_detail(
661 &IMPULSE_ERROR_INVALID_MODEL,
662 "denominator leading coefficient must be non-zero",
663 ));
664 }
665 let denominator: Vec<f64> = system
666 .denominator
667 .iter()
668 .map(|value| *value / leading)
669 .collect();
670 let mut numerator: Vec<f64> = system
671 .numerator
672 .iter()
673 .map(|value| *value / leading)
674 .collect();
675 let order = denominator.len() - 1;
676 while numerator.len() < order {
677 numerator.insert(0, 0.0);
678 }
679
680 let mut a = DMatrix::<f64>::zeros(order, order);
681 for row in 0..order.saturating_sub(1) {
682 a[(row, row + 1)] = 1.0;
683 }
684 for col in 0..order {
685 a[(order - 1, col)] = -denominator[order - col];
686 }
687 let c = numerator.into_iter().rev().collect();
688 Ok(Self { a, c })
689 }
690}
691
692fn continuous_response(realization: &Realization, t: &[f64]) -> Vec<f64> {
693 t.iter()
694 .map(|&time| {
695 let exp_at = matrix_exp(&(realization.a.clone() * time));
696 dot_c_with_last_column(&realization.c, &exp_at)
697 })
698 .collect()
699}
700
701fn discrete_response(
702 system: &TfSystem,
703 realization: &Realization,
704 t: &[f64],
705) -> BuiltinResult<Vec<f64>> {
706 if t.len() > MAX_DISCRETE_SAMPLES {
707 return Err(impulse_error_with_detail(
708 &IMPULSE_ERROR_DISCRETE_LIMIT,
709 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
710 ));
711 }
712 let sample_indices: Vec<usize> = t
713 .iter()
714 .map(|value| checked_discrete_sample_index(system, *value))
715 .collect::<BuiltinResult<_>>()?;
716 let max_index = sample_indices.iter().copied().max().unwrap_or(0);
717 let order = realization.c.len();
718 let value_count = max_index.checked_add(1).ok_or_else(|| {
719 impulse_error_with_detail(
720 &IMPULSE_ERROR_DISCRETE_LIMIT,
721 "discrete sample index exceeds platform limits",
722 )
723 })?;
724 let mut values = vec![0.0; value_count];
725 if order == 0 {
726 return Ok(sample_indices.into_iter().map(|idx| values[idx]).collect());
727 }
728
729 let mut state = vec![0.0; order];
730 state[order - 1] = 1.0;
731 let impulse_scale = 1.0 / system.sample_time;
732 for k in 1..=max_index {
733 values[k] = dot(&realization.c, &state) * impulse_scale;
734 state = mat_vec_mul(&realization.a, &state);
735 }
736 Ok(sample_indices.into_iter().map(|idx| values[idx]).collect())
737}
738
739fn dot_c_with_last_column(c: &[f64], matrix: &DMatrix<f64>) -> f64 {
740 if c.is_empty() {
741 return 0.0;
742 }
743 let last_col = matrix.ncols() - 1;
744 c.iter()
745 .enumerate()
746 .map(|(row, coeff)| coeff * matrix[(row, last_col)])
747 .sum()
748}
749
750fn dot(lhs: &[f64], rhs: &[f64]) -> f64 {
751 lhs.iter().zip(rhs).map(|(a, b)| a * b).sum()
752}
753
754fn mat_vec_mul(matrix: &DMatrix<f64>, vector: &[f64]) -> Vec<f64> {
755 let mut out = vec![0.0; matrix.nrows()];
756 for row in 0..matrix.nrows() {
757 let mut acc = 0.0;
758 for col in 0..matrix.ncols() {
759 acc += matrix[(row, col)] * vector[col];
760 }
761 out[row] = acc;
762 }
763 out
764}
765
766fn matrix_exp(matrix: &DMatrix<f64>) -> DMatrix<f64> {
767 let norm = matrix_one_norm(matrix);
768 let scale_power = if norm <= 0.5 {
769 0usize
770 } else {
771 norm.log2().ceil().max(0.0) as usize + 1
772 };
773 let scale = 2.0_f64.powi(scale_power as i32);
774 let scaled = matrix / scale;
775 let n = matrix.nrows();
776 let mut result = DMatrix::<f64>::identity(n, n);
777 let mut term = DMatrix::<f64>::identity(n, n);
778 for k in 1..=48 {
779 term = (&term * &scaled) / (k as f64);
780 result += &term;
781 if matrix_one_norm(&term) <= 1.0e-14 {
782 break;
783 }
784 }
785 for _ in 0..scale_power {
786 result = &result * &result;
787 }
788 result
789}
790
791fn matrix_one_norm(matrix: &DMatrix<f64>) -> f64 {
792 let mut best = 0.0;
793 for col in 0..matrix.ncols() {
794 let mut sum = 0.0;
795 for row in 0..matrix.nrows() {
796 sum += matrix[(row, col)].abs();
797 }
798 if sum > best {
799 best = sum;
800 }
801 }
802 best
803}
804
805#[cfg(feature = "plot-core")]
806async fn render_impulse_plot(response: &ImpulseResponse) -> BuiltinResult<Value> {
807 let t = response.t_value()?;
808 let y = response.y_value()?;
809 let plot_name = if response.discrete { "stem" } else { "plot" };
810 crate::dispatcher::call_builtin_async(plot_name, &[t, y]).await
811}
812
813#[cfg(not(feature = "plot-core"))]
814async fn render_impulse_plot(_response: &ImpulseResponse) -> BuiltinResult<Value> {
815 Ok(Value::Num(f64::NAN))
816}
817
818#[cfg(test)]
819mod tests {
820 use super::*;
821 use futures::executor::block_on;
822 use runmat_builtins::{CharArray, ObjectInstance};
823
824 fn tf_object(num: Vec<f64>, den: Vec<f64>, ts: f64) -> Value {
825 tf_object_with_delays(num, den, ts, 0.0, 0.0)
826 }
827
828 fn tf_object_with_delays(
829 num: Vec<f64>,
830 den: Vec<f64>,
831 ts: f64,
832 input_delay: f64,
833 output_delay: f64,
834 ) -> Value {
835 let mut object = ObjectInstance::new("tf".to_string());
836 object.properties.insert(
837 "Numerator".to_string(),
838 Value::Tensor(Tensor::new(num.clone(), vec![1, num.len()]).unwrap()),
839 );
840 object.properties.insert(
841 "Denominator".to_string(),
842 Value::Tensor(Tensor::new(den.clone(), vec![1, den.len()]).unwrap()),
843 );
844 object.properties.insert(
845 "Variable".to_string(),
846 Value::CharArray(CharArray::new_row(if ts > 0.0 { "z" } else { "s" })),
847 );
848 object.properties.insert("Ts".to_string(), Value::Num(ts));
849 object
850 .properties
851 .insert("InputDelay".to_string(), Value::Num(input_delay));
852 object
853 .properties
854 .insert("OutputDelay".to_string(), Value::Num(output_delay));
855 Value::Object(object)
856 }
857
858 fn run_impulse(system: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
859 block_on(impulse_builtin(system, rest))
860 }
861
862 fn tensor_data(value: Value) -> Vec<f64> {
863 match value {
864 Value::Tensor(tensor) => tensor.data,
865 other => panic!("expected tensor, got {other:?}"),
866 }
867 }
868
869 #[test]
870 fn impulse_descriptor_signatures_cover_core_forms() {
871 let labels: Vec<&str> = IMPULSE_DESCRIPTOR
872 .signatures
873 .iter()
874 .map(|sig| sig.label)
875 .collect();
876 assert!(labels.contains(&"y = impulse(sys)"));
877 assert!(labels.contains(&"y = impulse(sys, tFinal)"));
878 assert!(labels.contains(&"y = impulse(sys, t)"));
879 assert!(labels.contains(&"[y,t] = impulse(sys)"));
880 assert!(labels.contains(&"[y,t] = impulse(sys, tFinal)"));
881 assert!(labels.contains(&"[y,t] = impulse(sys, t)"));
882 }
883
884 #[test]
885 fn impulse_first_order_continuous_explicit_time() {
886 let sys = tf_object(vec![20.0], vec![1.0, 5.0], 0.0);
887 let t = Value::Tensor(Tensor::new(vec![0.0, 0.1, 0.2], vec![1, 3]).unwrap());
888 let y = tensor_data(run_impulse(sys, vec![t]).expect("impulse"));
889 let expected = [20.0, 20.0 * (-0.5f64).exp(), 20.0 * (-1.0f64).exp()];
890 for (actual, expected) in y.iter().zip(expected) {
891 assert!((actual - expected).abs() < 1.0e-8);
892 }
893 }
894
895 #[test]
896 fn impulse_second_order_continuous() {
897 let sys = tf_object(vec![1.0], vec![1.0, 3.0, 2.0], 0.0);
898 let t = Value::Tensor(Tensor::new(vec![0.0, 0.5, 1.0], vec![1, 3]).unwrap());
899 let y = tensor_data(run_impulse(sys, vec![t]).expect("impulse"));
900 for (actual, time) in y.iter().zip([0.0_f64, 0.5, 1.0]) {
901 let expected = (-time).exp() - (-2.0 * time).exp();
902 assert!((actual - expected).abs() < 1.0e-8);
903 }
904 }
905
906 #[test]
907 fn impulse_multi_output_returns_y_and_t_columns() {
908 let _guard = crate::output_count::push_output_count(Some(2));
909 let sys = tf_object(vec![20.0], vec![1.0, 5.0], 0.0);
910 let t_arg = Value::Tensor(Tensor::new(vec![0.0, 0.1], vec![1, 2]).unwrap());
911 let result = run_impulse(sys, vec![t_arg]).expect("impulse");
912 let Value::OutputList(outputs) = result else {
913 panic!("expected output list");
914 };
915 assert_eq!(outputs.len(), 2);
916 match &outputs[0] {
917 Value::Tensor(tensor) => assert_eq!(tensor.shape, vec![2, 1]),
918 other => panic!("expected y tensor, got {other:?}"),
919 }
920 match &outputs[1] {
921 Value::Tensor(tensor) => {
922 assert_eq!(tensor.shape, vec![2, 1]);
923 assert_eq!(tensor.data, vec![0.0, 0.1]);
924 }
925 other => panic!("expected t tensor, got {other:?}"),
926 }
927 }
928
929 #[test]
930 fn impulse_zero_output_count_emits_no_values() {
931 let _guard = crate::output_count::push_output_count(Some(0));
932 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
933 let result = run_impulse(sys, Vec::new()).expect("impulse");
934 let Value::OutputList(outputs) = result else {
935 panic!("expected output list");
936 };
937 assert!(outputs.is_empty());
938 }
939
940 #[test]
941 fn impulse_requested_zero_outputs_emits_no_values() {
942 let _guard = crate::output_context::push_output_count(0);
943 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
944 let result = run_impulse(sys, Vec::new()).expect("impulse");
945 let Value::OutputList(outputs) = result else {
946 panic!("expected output list");
947 };
948 assert!(outputs.is_empty());
949 }
950
951 #[test]
952 fn impulse_discrete_siso_response() {
953 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 0.1);
954 let t = Value::Tensor(Tensor::new(vec![0.0, 0.1, 0.2, 0.3], vec![1, 4]).unwrap());
955 let y = tensor_data(run_impulse(sys, vec![t]).expect("impulse"));
956 assert_eq!(y.len(), 4);
957 assert!((y[0] - 0.0).abs() < 1.0e-12);
958 assert!((y[1] - 10.0).abs() < 1.0e-12);
959 assert!((y[2] - 5.0).abs() < 1.0e-12);
960 assert!((y[3] - 2.5).abs() < 1.0e-12);
961 }
962
963 #[test]
964 fn impulse_discrete_final_time_rejects_excessive_sample_count() {
965 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0e-6);
966 let err = run_impulse(sys, vec![Value::Num(2.0)]).expect_err("should fail");
967 assert!(err.message().contains("more than 1000000 samples"));
968 assert_eq!(err.identifier(), IMPULSE_ERROR_DISCRETE_LIMIT.identifier);
969 }
970
971 #[test]
972 fn impulse_discrete_time_vector_rejects_excessive_sample_index() {
973 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0);
974 let t =
975 Value::Tensor(Tensor::new(vec![0.0, MAX_DISCRETE_SAMPLES as f64], vec![1, 2]).unwrap());
976 let err = run_impulse(sys, vec![t]).expect_err("should fail");
977 assert!(err.message().contains("more than 1000000 samples"));
978 }
979
980 #[test]
981 fn impulse_rejects_unsupported_model_type() {
982 let object = ObjectInstance::new("ss".to_string());
983 let err = run_impulse(Value::Object(object), Vec::new()).expect_err("should fail");
984 assert!(err.message().contains("unsupported model class"));
985 assert_eq!(err.identifier(), IMPULSE_ERROR_UNSUPPORTED_MODEL.identifier);
986 }
987
988 #[test]
989 fn impulse_rejects_direct_feedthrough_tf() {
990 let sys = tf_object(vec![1.0, 1.0], vec![1.0, 2.0], 0.0);
991 let err = run_impulse(sys, Vec::new()).expect_err("should fail");
992 assert!(err.message().contains("strictly proper"));
993 }
994
995 #[test]
996 fn impulse_rejects_invalid_time_metadata() {
997 let err = run_impulse(tf_object(vec![1.0], vec![1.0, -0.5], -0.1), Vec::new())
998 .expect_err("negative sample time should fail");
999 assert!(err.message().contains("Ts must be"));
1000
1001 let err = run_impulse(
1002 tf_object_with_delays(vec![1.0], vec![1.0, 5.0], 0.0, f64::NAN, 0.0),
1003 Vec::new(),
1004 )
1005 .expect_err("NaN input delay should fail");
1006 assert!(err.message().contains("InputDelay must be"));
1007
1008 let err = run_impulse(
1009 tf_object_with_delays(vec![1.0], vec![1.0, 5.0], 0.0, 0.0, -1.0),
1010 Vec::new(),
1011 )
1012 .expect_err("negative output delay should fail");
1013 assert!(err.message().contains("OutputDelay must be"));
1014 }
1015}