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 discrete: bool,
602}
603
604impl ImpulseResponse {
605 fn y_value(&self) -> BuiltinResult<Value> {
606 let tensor = Tensor::new(self.y.clone(), vec![self.y.len(), 1]).map_err(|err| {
607 impulse_error_with_detail(
608 &IMPULSE_ERROR_INTERNAL,
609 format!("failed to build response tensor: {err}"),
610 )
611 })?;
612 Ok(Value::Tensor(tensor))
613 }
614
615 fn t_value(&self) -> BuiltinResult<Value> {
616 let tensor = Tensor::new(self.t.clone(), vec![self.t.len(), 1]).map_err(|err| {
617 impulse_error_with_detail(
618 &IMPULSE_ERROR_INTERNAL,
619 format!("failed to build time tensor: {err}"),
620 )
621 })?;
622 Ok(Value::Tensor(tensor))
623 }
624}
625
626fn evaluate_impulse(system: &TfSystem, time: TimeSpec) -> BuiltinResult<ImpulseResponse> {
627 let TimeSpec::Values(t) = time;
628 let realization = Realization::from_tf(system)?;
629 let y = if system.is_discrete() {
630 discrete_response(system, &realization, &t)?
631 } else {
632 continuous_response(&realization, &t)
633 };
634 Ok(ImpulseResponse {
635 t,
636 y,
637 discrete: system.is_discrete(),
638 })
639}
640
641#[derive(Clone, Debug)]
642struct Realization {
643 a: DMatrix<f64>,
644 c: Vec<f64>,
645}
646
647impl Realization {
648 fn from_tf(system: &TfSystem) -> BuiltinResult<Self> {
649 if system.numerator.is_empty() {
650 let order = system.denominator.len().saturating_sub(1).max(1);
651 return Ok(Self {
652 a: DMatrix::zeros(order, order),
653 c: vec![0.0; order],
654 });
655 }
656 let leading = system.denominator[0];
657 if leading.abs() <= EPS {
658 return Err(impulse_error_with_detail(
659 &IMPULSE_ERROR_INVALID_MODEL,
660 "denominator leading coefficient must be non-zero",
661 ));
662 }
663 let denominator: Vec<f64> = system
664 .denominator
665 .iter()
666 .map(|value| *value / leading)
667 .collect();
668 let mut numerator: Vec<f64> = system
669 .numerator
670 .iter()
671 .map(|value| *value / leading)
672 .collect();
673 let order = denominator.len() - 1;
674 while numerator.len() < order {
675 numerator.insert(0, 0.0);
676 }
677
678 let mut a = DMatrix::<f64>::zeros(order, order);
679 for row in 0..order.saturating_sub(1) {
680 a[(row, row + 1)] = 1.0;
681 }
682 for col in 0..order {
683 a[(order - 1, col)] = -denominator[order - col];
684 }
685 let c = numerator.into_iter().rev().collect();
686 Ok(Self { a, c })
687 }
688}
689
690fn continuous_response(realization: &Realization, t: &[f64]) -> Vec<f64> {
691 t.iter()
692 .map(|&time| {
693 let exp_at = matrix_exp(&(realization.a.clone() * time));
694 dot_c_with_last_column(&realization.c, &exp_at)
695 })
696 .collect()
697}
698
699fn discrete_response(
700 system: &TfSystem,
701 realization: &Realization,
702 t: &[f64],
703) -> BuiltinResult<Vec<f64>> {
704 if t.len() > MAX_DISCRETE_SAMPLES {
705 return Err(impulse_error_with_detail(
706 &IMPULSE_ERROR_DISCRETE_LIMIT,
707 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
708 ));
709 }
710 let sample_indices: Vec<usize> = t
711 .iter()
712 .map(|value| checked_discrete_sample_index(system, *value))
713 .collect::<BuiltinResult<_>>()?;
714 let max_index = sample_indices.iter().copied().max().unwrap_or(0);
715 let order = realization.c.len();
716 let value_count = max_index.checked_add(1).ok_or_else(|| {
717 impulse_error_with_detail(
718 &IMPULSE_ERROR_DISCRETE_LIMIT,
719 "discrete sample index exceeds platform limits",
720 )
721 })?;
722 let mut values = vec![0.0; value_count];
723 if order == 0 {
724 return Ok(sample_indices.into_iter().map(|idx| values[idx]).collect());
725 }
726
727 let mut state = vec![0.0; order];
728 state[order - 1] = 1.0;
729 let impulse_scale = 1.0 / system.sample_time;
730 for k in 1..=max_index {
731 values[k] = dot(&realization.c, &state) * impulse_scale;
732 state = mat_vec_mul(&realization.a, &state);
733 }
734 Ok(sample_indices.into_iter().map(|idx| values[idx]).collect())
735}
736
737fn dot_c_with_last_column(c: &[f64], matrix: &DMatrix<f64>) -> f64 {
738 if c.is_empty() {
739 return 0.0;
740 }
741 let last_col = matrix.ncols() - 1;
742 c.iter()
743 .enumerate()
744 .map(|(row, coeff)| coeff * matrix[(row, last_col)])
745 .sum()
746}
747
748fn dot(lhs: &[f64], rhs: &[f64]) -> f64 {
749 lhs.iter().zip(rhs).map(|(a, b)| a * b).sum()
750}
751
752fn mat_vec_mul(matrix: &DMatrix<f64>, vector: &[f64]) -> Vec<f64> {
753 let mut out = vec![0.0; matrix.nrows()];
754 for row in 0..matrix.nrows() {
755 let mut acc = 0.0;
756 for col in 0..matrix.ncols() {
757 acc += matrix[(row, col)] * vector[col];
758 }
759 out[row] = acc;
760 }
761 out
762}
763
764fn matrix_exp(matrix: &DMatrix<f64>) -> DMatrix<f64> {
765 let norm = matrix_one_norm(matrix);
766 let scale_power = if norm <= 0.5 {
767 0usize
768 } else {
769 norm.log2().ceil().max(0.0) as usize + 1
770 };
771 let scale = 2.0_f64.powi(scale_power as i32);
772 let scaled = matrix / scale;
773 let n = matrix.nrows();
774 let mut result = DMatrix::<f64>::identity(n, n);
775 let mut term = DMatrix::<f64>::identity(n, n);
776 for k in 1..=48 {
777 term = (&term * &scaled) / (k as f64);
778 result += &term;
779 if matrix_one_norm(&term) <= 1.0e-14 {
780 break;
781 }
782 }
783 for _ in 0..scale_power {
784 result = &result * &result;
785 }
786 result
787}
788
789fn matrix_one_norm(matrix: &DMatrix<f64>) -> f64 {
790 let mut best = 0.0;
791 for col in 0..matrix.ncols() {
792 let mut sum = 0.0;
793 for row in 0..matrix.nrows() {
794 sum += matrix[(row, col)].abs();
795 }
796 if sum > best {
797 best = sum;
798 }
799 }
800 best
801}
802
803#[cfg(feature = "plot-core")]
804async fn render_impulse_plot(response: &ImpulseResponse) -> BuiltinResult<Value> {
805 let t = response.t_value()?;
806 let y = response.y_value()?;
807 let plot_name = if response.discrete { "stem" } else { "plot" };
808 crate::dispatcher::call_builtin_async(plot_name, &[t, y]).await
809}
810
811#[cfg(not(feature = "plot-core"))]
812async fn render_impulse_plot(_response: &ImpulseResponse) -> BuiltinResult<Value> {
813 Ok(Value::Num(f64::NAN))
814}
815
816#[cfg(test)]
817mod tests {
818 use super::*;
819 use futures::executor::block_on;
820 use runmat_builtins::{CharArray, ObjectInstance};
821
822 fn tf_object(num: Vec<f64>, den: Vec<f64>, ts: f64) -> Value {
823 tf_object_with_delays(num, den, ts, 0.0, 0.0)
824 }
825
826 fn tf_object_with_delays(
827 num: Vec<f64>,
828 den: Vec<f64>,
829 ts: f64,
830 input_delay: f64,
831 output_delay: f64,
832 ) -> Value {
833 let mut object = ObjectInstance::new("tf".to_string());
834 object.properties.insert(
835 "Numerator".to_string(),
836 Value::Tensor(Tensor::new(num.clone(), vec![1, num.len()]).unwrap()),
837 );
838 object.properties.insert(
839 "Denominator".to_string(),
840 Value::Tensor(Tensor::new(den.clone(), vec![1, den.len()]).unwrap()),
841 );
842 object.properties.insert(
843 "Variable".to_string(),
844 Value::CharArray(CharArray::new_row(if ts > 0.0 { "z" } else { "s" })),
845 );
846 object.properties.insert("Ts".to_string(), Value::Num(ts));
847 object
848 .properties
849 .insert("InputDelay".to_string(), Value::Num(input_delay));
850 object
851 .properties
852 .insert("OutputDelay".to_string(), Value::Num(output_delay));
853 Value::Object(object)
854 }
855
856 fn run_impulse(system: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
857 block_on(impulse_builtin(system, rest))
858 }
859
860 fn tensor_data(value: Value) -> Vec<f64> {
861 match value {
862 Value::Tensor(tensor) => tensor.data,
863 other => panic!("expected tensor, got {other:?}"),
864 }
865 }
866
867 #[test]
868 fn impulse_descriptor_signatures_cover_core_forms() {
869 let labels: Vec<&str> = IMPULSE_DESCRIPTOR
870 .signatures
871 .iter()
872 .map(|sig| sig.label)
873 .collect();
874 assert!(labels.contains(&"y = impulse(sys)"));
875 assert!(labels.contains(&"y = impulse(sys, tFinal)"));
876 assert!(labels.contains(&"y = impulse(sys, t)"));
877 assert!(labels.contains(&"[y,t] = impulse(sys)"));
878 assert!(labels.contains(&"[y,t] = impulse(sys, tFinal)"));
879 assert!(labels.contains(&"[y,t] = impulse(sys, t)"));
880 }
881
882 #[test]
883 fn impulse_first_order_continuous_explicit_time() {
884 let sys = tf_object(vec![20.0], vec![1.0, 5.0], 0.0);
885 let t = Value::Tensor(Tensor::new(vec![0.0, 0.1, 0.2], vec![1, 3]).unwrap());
886 let y = tensor_data(run_impulse(sys, vec![t]).expect("impulse"));
887 let expected = [20.0, 20.0 * (-0.5f64).exp(), 20.0 * (-1.0f64).exp()];
888 for (actual, expected) in y.iter().zip(expected) {
889 assert!((actual - expected).abs() < 1.0e-8);
890 }
891 }
892
893 #[test]
894 fn impulse_second_order_continuous() {
895 let sys = tf_object(vec![1.0], vec![1.0, 3.0, 2.0], 0.0);
896 let t = Value::Tensor(Tensor::new(vec![0.0, 0.5, 1.0], vec![1, 3]).unwrap());
897 let y = tensor_data(run_impulse(sys, vec![t]).expect("impulse"));
898 for (actual, time) in y.iter().zip([0.0_f64, 0.5, 1.0]) {
899 let expected = (-time).exp() - (-2.0 * time).exp();
900 assert!((actual - expected).abs() < 1.0e-8);
901 }
902 }
903
904 #[test]
905 fn impulse_multi_output_returns_y_and_t_columns() {
906 let _guard = crate::output_count::push_output_count(Some(2));
907 let sys = tf_object(vec![20.0], vec![1.0, 5.0], 0.0);
908 let t_arg = Value::Tensor(Tensor::new(vec![0.0, 0.1], vec![1, 2]).unwrap());
909 let result = run_impulse(sys, vec![t_arg]).expect("impulse");
910 let Value::OutputList(outputs) = result else {
911 panic!("expected output list");
912 };
913 assert_eq!(outputs.len(), 2);
914 match &outputs[0] {
915 Value::Tensor(tensor) => assert_eq!(tensor.shape, vec![2, 1]),
916 other => panic!("expected y tensor, got {other:?}"),
917 }
918 match &outputs[1] {
919 Value::Tensor(tensor) => {
920 assert_eq!(tensor.shape, vec![2, 1]);
921 assert_eq!(tensor.data, vec![0.0, 0.1]);
922 }
923 other => panic!("expected t tensor, got {other:?}"),
924 }
925 }
926
927 #[test]
928 fn impulse_zero_output_count_emits_no_values() {
929 let _guard = crate::output_count::push_output_count(Some(0));
930 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
931 let result = run_impulse(sys, Vec::new()).expect("impulse");
932 let Value::OutputList(outputs) = result else {
933 panic!("expected output list");
934 };
935 assert!(outputs.is_empty());
936 }
937
938 #[test]
939 fn impulse_requested_zero_outputs_emits_no_values() {
940 let _guard = crate::output_context::push_output_count(0);
941 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
942 let result = run_impulse(sys, Vec::new()).expect("impulse");
943 let Value::OutputList(outputs) = result else {
944 panic!("expected output list");
945 };
946 assert!(outputs.is_empty());
947 }
948
949 #[test]
950 fn impulse_discrete_siso_response() {
951 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 0.1);
952 let t = Value::Tensor(Tensor::new(vec![0.0, 0.1, 0.2, 0.3], vec![1, 4]).unwrap());
953 let y = tensor_data(run_impulse(sys, vec![t]).expect("impulse"));
954 assert_eq!(y.len(), 4);
955 assert!((y[0] - 0.0).abs() < 1.0e-12);
956 assert!((y[1] - 10.0).abs() < 1.0e-12);
957 assert!((y[2] - 5.0).abs() < 1.0e-12);
958 assert!((y[3] - 2.5).abs() < 1.0e-12);
959 }
960
961 #[test]
962 fn impulse_discrete_final_time_rejects_excessive_sample_count() {
963 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0e-6);
964 let err = run_impulse(sys, vec![Value::Num(2.0)]).expect_err("should fail");
965 assert!(err.message().contains("more than 1000000 samples"));
966 assert_eq!(err.identifier(), IMPULSE_ERROR_DISCRETE_LIMIT.identifier);
967 }
968
969 #[test]
970 fn impulse_discrete_time_vector_rejects_excessive_sample_index() {
971 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0);
972 let t =
973 Value::Tensor(Tensor::new(vec![0.0, MAX_DISCRETE_SAMPLES as f64], vec![1, 2]).unwrap());
974 let err = run_impulse(sys, vec![t]).expect_err("should fail");
975 assert!(err.message().contains("more than 1000000 samples"));
976 }
977
978 #[test]
979 fn impulse_rejects_unsupported_model_type() {
980 let object = ObjectInstance::new("ss".to_string());
981 let err = run_impulse(Value::Object(object), Vec::new()).expect_err("should fail");
982 assert!(err.message().contains("unsupported model class"));
983 assert_eq!(err.identifier(), IMPULSE_ERROR_UNSUPPORTED_MODEL.identifier);
984 }
985
986 #[test]
987 fn impulse_rejects_direct_feedthrough_tf() {
988 let sys = tf_object(vec![1.0, 1.0], vec![1.0, 2.0], 0.0);
989 let err = run_impulse(sys, Vec::new()).expect_err("should fail");
990 assert!(err.message().contains("strictly proper"));
991 }
992
993 #[test]
994 fn impulse_rejects_invalid_time_metadata() {
995 let err = run_impulse(tf_object(vec![1.0], vec![1.0, -0.5], -0.1), Vec::new())
996 .expect_err("negative sample time should fail");
997 assert!(err.message().contains("Ts must be"));
998
999 let err = run_impulse(
1000 tf_object_with_delays(vec![1.0], vec![1.0, 5.0], 0.0, f64::NAN, 0.0),
1001 Vec::new(),
1002 )
1003 .expect_err("NaN input delay should fail");
1004 assert!(err.message().contains("InputDelay must be"));
1005
1006 let err = run_impulse(
1007 tf_object_with_delays(vec![1.0], vec![1.0, 5.0], 0.0, 0.0, -1.0),
1008 Vec::new(),
1009 )
1010 .expect_err("negative output delay should fail");
1011 assert!(err.message().contains("OutputDelay must be"));
1012 }
1013}