1use nalgebra::DMatrix;
4use runmat_builtins::{
5 BuiltinCompletionPolicy, BuiltinDescriptor, BuiltinErrorDescriptor, BuiltinExtensionDescriptor,
6 BuiltinExtensionMode, BuiltinIntegerBackendRule, BuiltinIntegerCapabilityDescriptor,
7 BuiltinIntegerComputationDomain, BuiltinIntegerInputAvailability,
8 BuiltinIntegerInputCapability, BuiltinIntegerOutputClassRule, BuiltinIntegerOverflowRule,
9 BuiltinIntegerOverloadKind, BuiltinIntegerScalarDoubleRule, BuiltinOutputMode,
10 BuiltinParamArity, BuiltinParamDescriptor, BuiltinParamType, BuiltinSignatureDescriptor,
11};
12use runmat_macros::runtime_builtin;
13use runmat_value::{IntValue, NumericDType, NumericScalar};
14use runmat_value::{Tensor, Value};
15
16use crate::builtins::common::spec::{
17 BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
18 ReductionNaN, ResidencyPolicy, ShapeRequirements,
19};
20use crate::builtins::common::tensor;
21use crate::builtins::control::type_resolvers::impulse_type;
22use crate::{build_runtime_error, BuiltinResult, RuntimeError};
23
24const BUILTIN_NAME: &str = "impulse";
25const TF_CLASS: &str = "tf";
26const EPS: f64 = 1.0e-12;
27const DEFAULT_POINTS: usize = 100;
28const MAX_DISCRETE_SAMPLES: usize = 1_000_000;
29
30const IMPULSE_INTEGER_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
31 id: "impulse-integer-numeric-role",
32 mode: BuiltinExtensionMode::RunMatOnly,
33 description: "impulse with typed-integer model metadata or time input is a RunMat extension",
34 error_identifier: Some("RunMat:compatibility:ImpulseIntegerExtension"),
35};
36const IMPULSE_LOGICAL_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
37 id: "impulse-logical-numeric-role",
38 mode: BuiltinExtensionMode::RunMatOnly,
39 description: "impulse with logical model metadata or time input is a RunMat extension",
40 error_identifier: Some("RunMat:compatibility:ImpulseLogicalExtension"),
41};
42const IMPULSE_SINGLE_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
43 id: "impulse-single-numeric-role",
44 mode: BuiltinExtensionMode::RunMatOnly,
45 description: "impulse with single-precision model metadata or time input is a RunMat extension",
46 error_identifier: Some("RunMat:compatibility:ImpulseSingleExtension"),
47};
48const IMPULSE_RESIDENT_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
49 id: "impulse-resident-input",
50 mode: BuiltinExtensionMode::RunMatOnly,
51 description: "impulse with resident model metadata or time input is a RunMat extension",
52 error_identifier: Some("RunMat:compatibility:ImpulseResidentExtension"),
53};
54pub const IMPULSE_EXTENSIONS: [BuiltinExtensionDescriptor; 4] = [
55 IMPULSE_INTEGER_EXTENSION,
56 IMPULSE_LOGICAL_EXTENSION,
57 IMPULSE_SINGLE_EXTENSION,
58 IMPULSE_RESIDENT_EXTENSION,
59];
60
61const IMPULSE_INTEGER_INPUTS: [BuiltinIntegerInputCapability; 2] = [
62 BuiltinIntegerInputCapability {
63 name: "sys numeric metadata",
64 classes: &crate::builtins::common::integer_capability::ALL_INTEGER_CLASSES,
65 availability: BuiltinIntegerInputAvailability::RunMatOnly,
66 scalar_double: BuiltinIntegerScalarDoubleRule::NotApplicable,
67 notes: "Typed-integer coefficients and timing metadata are outside the documented dynamic-system/time surface, are gated before model parsing, and must be exactly representable at the host floating simulation boundary when the RunMat extension is enabled.",
68 },
69 BuiltinIntegerInputCapability {
70 name: "t or tFinal",
71 classes: &crate::builtins::common::integer_capability::ALL_INTEGER_CLASSES,
72 availability: BuiltinIntegerInputAvailability::RunMatOnly,
73 scalar_double: BuiltinIntegerScalarDoubleRule::NotApplicable,
74 notes: "The public impulse reference specifies numeric time values and shapes but no typed-integer class support, so every typed-integer time control is consistently classified as a RunMat-only extension.",
75 },
76];
77pub const IMPULSE_INTEGER_CAPABILITIES: [BuiltinIntegerCapabilityDescriptor; 1] =
78 [BuiltinIntegerCapabilityDescriptor {
79 form: "[y, t] = impulse(sys, integer_time?)",
80 inputs: &IMPULSE_INTEGER_INPUTS,
81 computation_domain: BuiltinIntegerComputationDomain::FloatingPoint,
82 output_class: BuiltinIntegerOutputClassRule::Double,
83 overflow: BuiltinIntegerOverflowRule::Error,
84 backend: BuiltinIntegerBackendRule::GatherFallback,
85 overload: BuiltinIntegerOverloadKind::Multiple,
86 notes: "Non-double and resident numeric roles are independent RunMat extensions. Typed integers cross to binary64 only after exact-representability validation; this bounded SISO tf implementation simulates on the host and returns double response/time arrays.",
87 }];
88
89const IMPULSE_OUTPUT_Y: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
90 name: "y",
91 ty: BuiltinParamType::NumericArray,
92 arity: BuiltinParamArity::Required,
93 default: None,
94 description: "Impulse response samples (column vector).",
95}];
96const IMPULSE_OUTPUT_Y_T: [BuiltinParamDescriptor; 2] = [
97 BuiltinParamDescriptor {
98 name: "y",
99 ty: BuiltinParamType::NumericArray,
100 arity: BuiltinParamArity::Required,
101 default: None,
102 description: "Impulse response samples (column vector).",
103 },
104 BuiltinParamDescriptor {
105 name: "t",
106 ty: BuiltinParamType::NumericArray,
107 arity: BuiltinParamArity::Required,
108 default: None,
109 description: "Time samples (column vector).",
110 },
111];
112const IMPULSE_INPUTS_SYS: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
113 name: "sys",
114 ty: BuiltinParamType::Any,
115 arity: BuiltinParamArity::Required,
116 default: None,
117 description: "SISO tf model.",
118}];
119const IMPULSE_INPUTS_SYS_TIME: [BuiltinParamDescriptor; 2] = [
120 BuiltinParamDescriptor {
121 name: "sys",
122 ty: BuiltinParamType::Any,
123 arity: BuiltinParamArity::Required,
124 default: None,
125 description: "SISO tf model.",
126 },
127 BuiltinParamDescriptor {
128 name: "time",
129 ty: BuiltinParamType::Any,
130 arity: BuiltinParamArity::Optional,
131 default: None,
132 description: "Final time scalar or explicit time vector.",
133 },
134];
135const IMPULSE_SIGNATURES: [BuiltinSignatureDescriptor; 6] = [
136 BuiltinSignatureDescriptor {
137 label: "y = impulse(sys)",
138 inputs: &IMPULSE_INPUTS_SYS,
139 outputs: &IMPULSE_OUTPUT_Y,
140 },
141 BuiltinSignatureDescriptor {
142 label: "y = impulse(sys, tFinal)",
143 inputs: &IMPULSE_INPUTS_SYS_TIME,
144 outputs: &IMPULSE_OUTPUT_Y,
145 },
146 BuiltinSignatureDescriptor {
147 label: "y = impulse(sys, t)",
148 inputs: &IMPULSE_INPUTS_SYS_TIME,
149 outputs: &IMPULSE_OUTPUT_Y,
150 },
151 BuiltinSignatureDescriptor {
152 label: "[y,t] = impulse(sys)",
153 inputs: &IMPULSE_INPUTS_SYS,
154 outputs: &IMPULSE_OUTPUT_Y_T,
155 },
156 BuiltinSignatureDescriptor {
157 label: "[y,t] = impulse(sys, tFinal)",
158 inputs: &IMPULSE_INPUTS_SYS_TIME,
159 outputs: &IMPULSE_OUTPUT_Y_T,
160 },
161 BuiltinSignatureDescriptor {
162 label: "[y,t] = impulse(sys, t)",
163 inputs: &IMPULSE_INPUTS_SYS_TIME,
164 outputs: &IMPULSE_OUTPUT_Y_T,
165 },
166];
167const IMPULSE_ERROR_INVALID_ARGUMENT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
168 code: "RM.IMPULSE.INVALID_ARGUMENT",
169 identifier: Some("RunMat:impulse:InvalidArgument"),
170 when: "Inputs do not match supported impulse invocation forms.",
171 message: "impulse: invalid argument",
172};
173const IMPULSE_ERROR_INVALID_MODEL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
174 code: "RM.IMPULSE.INVALID_MODEL",
175 identifier: Some("RunMat:impulse:InvalidModel"),
176 when: "Input system is not a supported tf object with valid required metadata.",
177 message: "impulse: invalid model",
178};
179const IMPULSE_ERROR_INVALID_TIME: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
180 code: "RM.IMPULSE.INVALID_TIME",
181 identifier: Some("RunMat:impulse:InvalidTime"),
182 when: "Time argument is invalid for the model class or sampling mode.",
183 message: "impulse: invalid time input",
184};
185const IMPULSE_ERROR_UNSUPPORTED_MODEL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
186 code: "RM.IMPULSE.UNSUPPORTED_MODEL",
187 identifier: Some("RunMat:impulse:UnsupportedModel"),
188 when: "Model is well-formed but unsupported by the current impulse implementation.",
189 message: "impulse: unsupported model",
190};
191const IMPULSE_ERROR_DISCRETE_LIMIT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
192 code: "RM.IMPULSE.DISCRETE_LIMIT",
193 identifier: Some("RunMat:impulse:DiscreteLimit"),
194 when: "Discrete simulation would exceed platform or configured sample limits.",
195 message: "impulse: discrete simulation limit exceeded",
196};
197const IMPULSE_ERROR_INTERNAL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
198 code: "RM.IMPULSE.INTERNAL",
199 identifier: Some("RunMat:impulse:Internal"),
200 when: "Internal response assembly failed.",
201 message: "impulse: internal error",
202};
203const IMPULSE_ERRORS: [BuiltinErrorDescriptor; 6] = [
204 IMPULSE_ERROR_INVALID_ARGUMENT,
205 IMPULSE_ERROR_INVALID_MODEL,
206 IMPULSE_ERROR_INVALID_TIME,
207 IMPULSE_ERROR_UNSUPPORTED_MODEL,
208 IMPULSE_ERROR_DISCRETE_LIMIT,
209 IMPULSE_ERROR_INTERNAL,
210];
211pub const IMPULSE_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
212 signatures: &IMPULSE_SIGNATURES,
213 output_mode: BuiltinOutputMode::ByRequestedOutputCount,
214 completion_policy: BuiltinCompletionPolicy::Public,
215 errors: &IMPULSE_ERRORS,
216};
217
218#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::control::impulse")]
219pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
220 name: "impulse",
221 op_kind: GpuOpKind::Custom("control-impulse-response"),
222 supported_precisions: &[],
223 broadcast: BroadcastSemantics::None,
224 provider_hooks: &[],
225 constant_strategy: ConstantStrategy::InlineLiteral,
226 residency: ResidencyPolicy::GatherImmediately,
227 nan_mode: ReductionNaN::Include,
228 two_pass_threshold: None,
229 workgroup_size: None,
230 accepts_nan_mode: false,
231 notes: "Control-system response evaluation runs on the host. GPU-resident metadata is gathered before simulation.",
232};
233
234#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::control::impulse")]
235pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
236 name: "impulse",
237 shape: ShapeRequirements::Any,
238 constant_strategy: ConstantStrategy::InlineLiteral,
239 elementwise: None,
240 reduction: None,
241 emits_nan: false,
242 notes: "Impulse-response simulation materialises host-side time and output vectors and terminates fusion chains.",
243};
244
245fn impulse_error_with_detail(
246 error: &'static BuiltinErrorDescriptor,
247 detail: impl AsRef<str>,
248) -> RuntimeError {
249 impulse_error_with_message(format!("{}: {}", error.message, detail.as_ref()), error)
250}
251
252fn impulse_error_with_message(
253 message: impl Into<String>,
254 error: &'static BuiltinErrorDescriptor,
255) -> RuntimeError {
256 let mut builder = build_runtime_error(message).with_builtin(BUILTIN_NAME);
257 if let Some(identifier) = error.identifier {
258 builder = builder.with_identifier(identifier);
259 }
260 builder.build()
261}
262
263#[runtime_builtin(
264 name = "impulse",
265 category = "control",
266 summary = "Compute or plot impulse responses.",
267 keywords = "impulse,control system,transfer function,response,tf",
268 type_resolver(impulse_type),
269 descriptor(crate::builtins::control::impulse::IMPULSE_DESCRIPTOR),
270 extensions(crate::builtins::control::impulse::IMPULSE_EXTENSIONS),
271 integer_capabilities(crate::builtins::control::impulse::IMPULSE_INTEGER_CAPABILITIES),
272 builtin_path = "crate::builtins::control::impulse"
273)]
274async fn impulse_builtin(system: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
275 ensure_impulse_extensions(&system, &rest)?;
276 let system = TfSystem::parse(system).await?;
277 let time = TimeSpec::parse(&system, &rest).await?;
278 let response = evaluate_impulse(&system, time)?;
279
280 if let Some(out_count) = crate::output_count::current_output_count() {
281 if out_count == 0 {
282 emit_impulse_plot(&response).await?;
283 return Ok(Value::OutputList(Vec::new()));
284 }
285 if out_count == 1 {
286 return Ok(Value::OutputList(vec![response.y_value()?]));
287 }
288 return Ok(crate::output_count::output_list_with_padding(
289 out_count,
290 vec![response.y_value()?, response.t_value()?],
291 ));
292 }
293
294 if crate::output_context::requested_output_count() == Some(0) {
295 emit_impulse_plot(&response).await?;
296 return Ok(Value::OutputList(Vec::new()));
297 }
298
299 response.y_value()
300}
301
302fn ensure_impulse_extensions(system: &Value, rest: &[Value]) -> BuiltinResult<()> {
303 let mut values = std::iter::once(system).chain(rest.iter());
304 let integer = values.clone().any(value_contains_typed_integer);
305 let logical = values.clone().any(value_contains_logical);
306 let single = values.clone().any(value_contains_single);
307 let resident = values.any(crate::dispatcher::value_contains_gpu);
308 if integer {
309 crate::compatibility::ensure_builtin_extension_enabled(
310 &IMPULSE_INTEGER_EXTENSION,
311 BUILTIN_NAME,
312 )?;
313 }
314 if logical {
315 crate::compatibility::ensure_builtin_extension_enabled(
316 &IMPULSE_LOGICAL_EXTENSION,
317 BUILTIN_NAME,
318 )?;
319 }
320 if single {
321 crate::compatibility::ensure_builtin_extension_enabled(
322 &IMPULSE_SINGLE_EXTENSION,
323 BUILTIN_NAME,
324 )?;
325 }
326 if resident {
327 crate::compatibility::ensure_builtin_extension_enabled(
328 &IMPULSE_RESIDENT_EXTENSION,
329 BUILTIN_NAME,
330 )?;
331 }
332 Ok(())
333}
334
335fn value_contains_typed_integer(value: &Value) -> bool {
336 match value {
337 Value::Int(_) => true,
338 Value::Tensor(tensor) => tensor.integer_storage().is_some(),
339 Value::GpuTensor(handle) => runmat_accelerate_api::handle_integer_type(handle).is_some(),
340 Value::Object(object) => object.properties.values().any(value_contains_typed_integer),
341 _ => false,
342 }
343}
344
345fn value_contains_logical(value: &Value) -> bool {
346 match value {
347 Value::Bool(_) | Value::LogicalArray(_) => true,
348 Value::GpuTensor(handle) => runmat_accelerate_api::handle_is_logical(handle),
349 Value::Object(object) => object.properties.values().any(value_contains_logical),
350 _ => false,
351 }
352}
353
354fn value_contains_single(value: &Value) -> bool {
355 match value {
356 Value::Tensor(tensor) => tensor.numeric_dtype() == NumericDType::F32,
357 Value::GpuTensor(handle) => {
358 runmat_accelerate_api::handle_integer_type(handle).is_none()
359 && !runmat_accelerate_api::handle_is_logical(handle)
360 && runmat_accelerate_api::handle_precision(handle)
361 == Some(runmat_accelerate_api::ProviderPrecision::F32)
362 }
363 Value::Object(object) => object.properties.values().any(value_contains_single),
364 _ => false,
365 }
366}
367
368async fn emit_impulse_plot(response: &ImpulseResponse) -> BuiltinResult<()> {
369 if let Err(err) = render_impulse_plot(response).await {
370 if is_nonfatal_plot_setup_error(&err) {
371 return Ok(());
372 }
373 return Err(err);
374 }
375 Ok(())
376}
377
378fn is_nonfatal_plot_setup_error(err: &RuntimeError) -> bool {
379 let lower = err.message().to_ascii_lowercase();
380 lower.contains("plotting is unavailable")
381 || lower.contains("non-main thread")
382 || lower.contains("interactive plotting failed")
383}
384
385#[derive(Clone, Debug)]
386struct TfSystem {
387 numerator: Vec<f64>,
388 denominator: Vec<f64>,
389 sample_time: f64,
390}
391
392impl TfSystem {
393 async fn parse(value: Value) -> BuiltinResult<Self> {
394 let gathered = crate::dispatcher::gather_if_needed_async(&value).await?;
395 let Value::Object(object) = gathered else {
396 return Err(impulse_error_with_detail(
397 &IMPULSE_ERROR_INVALID_MODEL,
398 format!("expected a dynamic system model, got {gathered:?}"),
399 ));
400 };
401 if object.class_name != TF_CLASS {
402 return Err(impulse_error_with_detail(
403 &IMPULSE_ERROR_UNSUPPORTED_MODEL,
404 format!(
405 "unsupported model class '{}'; only SISO tf objects are currently supported",
406 object.class_name
407 ),
408 ));
409 }
410
411 let numerator = real_coefficients(property(&object, "Numerator")?, "Numerator")?;
412 let denominator = real_coefficients(property(&object, "Denominator")?, "Denominator")?;
413 let sample_time = scalar_property(property(&object, "Ts")?, "Ts")?;
414 let input_delay = scalar_property(property(&object, "InputDelay")?, "InputDelay")?;
415 let output_delay = scalar_property(property(&object, "OutputDelay")?, "OutputDelay")?;
416 if !sample_time.is_finite() || sample_time < 0.0 {
417 return Err(impulse_error_with_detail(
418 &IMPULSE_ERROR_INVALID_MODEL,
419 format!("Ts must be a finite non-negative scalar, got {sample_time}"),
420 ));
421 }
422 if !input_delay.is_finite() || input_delay < 0.0 {
423 return Err(impulse_error_with_detail(
424 &IMPULSE_ERROR_INVALID_MODEL,
425 format!("InputDelay must be a finite non-negative scalar, got {input_delay}"),
426 ));
427 }
428 if !output_delay.is_finite() || output_delay < 0.0 {
429 return Err(impulse_error_with_detail(
430 &IMPULSE_ERROR_INVALID_MODEL,
431 format!("OutputDelay must be a finite non-negative scalar, got {output_delay}"),
432 ));
433 }
434 if input_delay.abs() > EPS || output_delay.abs() > EPS {
435 return Err(impulse_error_with_detail(
436 &IMPULSE_ERROR_UNSUPPORTED_MODEL,
437 "transfer functions with input or output delays are not supported yet",
438 ));
439 }
440
441 let numerator = trim_leading_zeros(numerator);
442 let denominator = trim_leading_zeros(denominator);
443 if denominator.is_empty() {
444 return Err(impulse_error_with_detail(
445 &IMPULSE_ERROR_INVALID_MODEL,
446 "denominator coefficients cannot be empty",
447 ));
448 }
449 if numerator.is_empty() {
450 return Ok(Self {
451 numerator,
452 denominator,
453 sample_time,
454 });
455 }
456 if denominator.len() <= 1 {
457 return Err(impulse_error_with_detail(
458 &IMPULSE_ERROR_UNSUPPORTED_MODEL,
459 "static-gain transfer functions do not have a finite impulse-response vector",
460 ));
461 }
462 if numerator.len() >= denominator.len() {
463 return Err(impulse_error_with_detail(
464 &IMPULSE_ERROR_UNSUPPORTED_MODEL,
465 "only strictly proper SISO tf models are currently supported",
466 ));
467 }
468
469 Ok(Self {
470 numerator,
471 denominator,
472 sample_time,
473 })
474 }
475
476 fn is_discrete(&self) -> bool {
477 self.sample_time > 0.0
478 }
479}
480
481fn property<'a>(object: &'a runmat_value::ObjectInstance, name: &str) -> BuiltinResult<&'a Value> {
482 object.properties.get(name).ok_or_else(|| {
483 impulse_error_with_detail(
484 &IMPULSE_ERROR_INVALID_MODEL,
485 format!("tf object is missing {name} property"),
486 )
487 })
488}
489
490fn real_coefficients(value: &Value, label: &str) -> BuiltinResult<Vec<f64>> {
491 match value {
492 Value::Tensor(tensor) => {
493 ensure_vector(label, &tensor.shape)?;
494 ensure_exact_integer_tensor(tensor, label)?;
495 finite_values(label, tensor::tensor_values_f64(tensor))
496 }
497 Value::Num(n) => finite_values(label, vec![*n]),
498 Value::Int(i) => finite_values(label, vec![exact_integer_as_f64(i, label)?]),
499 Value::Bool(b) => finite_values(label, vec![if *b { 1.0 } else { 0.0 }]),
500 Value::LogicalArray(logical) => {
501 ensure_vector(label, &logical.shape)?;
502 finite_values(
503 label,
504 logical
505 .data
506 .iter()
507 .map(|bit| if *bit == 0 { 0.0 } else { 1.0 })
508 .collect(),
509 )
510 }
511 Value::Complex(_, _) | Value::ComplexTensor(_) => Err(impulse_error_with_detail(
512 &IMPULSE_ERROR_UNSUPPORTED_MODEL,
513 "complex-coefficient tf models are not supported yet",
514 )),
515 other => Err(impulse_error_with_detail(
516 &IMPULSE_ERROR_INVALID_MODEL,
517 format!("{label} must be a real numeric coefficient vector, got {other:?}"),
518 )),
519 }
520}
521
522fn finite_values(label: &str, values: Vec<f64>) -> BuiltinResult<Vec<f64>> {
523 if values.iter().any(|value| !value.is_finite()) {
524 return Err(impulse_error_with_detail(
525 &IMPULSE_ERROR_INVALID_MODEL,
526 format!("{label} coefficients must be finite"),
527 ));
528 }
529 Ok(values)
530}
531
532fn ensure_vector(label: &str, shape: &[usize]) -> BuiltinResult<()> {
533 let non_unit = shape.iter().copied().filter(|&dim| dim > 1).count();
534 if non_unit <= 1 {
535 Ok(())
536 } else {
537 Err(impulse_error_with_detail(
538 &IMPULSE_ERROR_INVALID_MODEL,
539 format!("{label} coefficients must be a vector"),
540 ))
541 }
542}
543
544fn scalar_property(value: &Value, label: &str) -> BuiltinResult<f64> {
545 match value {
546 Value::Num(n) => Ok(*n),
547 Value::Int(i) => exact_integer_as_f64(i, label),
548 Value::Bool(b) => Ok(if *b { 1.0 } else { 0.0 }),
549 Value::Tensor(tensor) if tensor::is_scalar_tensor(tensor) => {
550 ensure_exact_integer_tensor(tensor, label)?;
551 Ok(tensor::tensor_value_f64(tensor, 0))
552 }
553 other => Err(impulse_error_with_detail(
554 &IMPULSE_ERROR_INVALID_MODEL,
555 format!("{label} must be a real scalar, got {other:?}"),
556 )),
557 }
558}
559
560fn exact_integer_as_f64(value: &IntValue, label: &str) -> BuiltinResult<f64> {
561 const MAX_EXACT: i128 = 1_i128 << 53;
562 let exact = match value {
563 IntValue::I8(value) => i128::from(*value),
564 IntValue::I16(value) => i128::from(*value),
565 IntValue::I32(value) => i128::from(*value),
566 IntValue::I64(value) => i128::from(*value),
567 IntValue::U8(value) => i128::from(*value),
568 IntValue::U16(value) => i128::from(*value),
569 IntValue::U32(value) => i128::from(*value),
570 IntValue::U64(value) => i128::from(*value),
571 };
572 if !(-MAX_EXACT..=MAX_EXACT).contains(&exact) {
573 return Err(impulse_error_with_detail(
574 &IMPULSE_ERROR_INVALID_ARGUMENT,
575 format!("{label} integer values must be exactly representable as double"),
576 ));
577 }
578 Ok(exact as f64)
579}
580
581fn ensure_exact_integer_tensor(tensor: &Tensor, label: &str) -> BuiltinResult<()> {
582 if tensor.integer_storage().is_none() {
583 return Ok(());
584 }
585 for index in 0..tensor.len() {
586 let value = match tensor.numeric_value_at(index) {
587 Some(NumericScalar::I8(value)) => IntValue::I8(value),
588 Some(NumericScalar::I16(value)) => IntValue::I16(value),
589 Some(NumericScalar::I32(value)) => IntValue::I32(value),
590 Some(NumericScalar::I64(value)) => IntValue::I64(value),
591 Some(NumericScalar::U8(value)) => IntValue::U8(value),
592 Some(NumericScalar::U16(value)) => IntValue::U16(value),
593 Some(NumericScalar::U32(value)) => IntValue::U32(value),
594 Some(NumericScalar::U64(value)) => IntValue::U64(value),
595 _ => continue,
596 };
597 exact_integer_as_f64(&value, label)?;
598 }
599 Ok(())
600}
601
602fn trim_leading_zeros(values: Vec<f64>) -> Vec<f64> {
603 let first = values.iter().position(|value| value.abs() > EPS);
604 match first {
605 Some(idx) => values[idx..].to_vec(),
606 None => Vec::new(),
607 }
608}
609
610#[derive(Clone, Debug)]
611enum TimeSpec {
612 Values(Vec<f64>),
613}
614
615impl TimeSpec {
616 async fn parse(system: &TfSystem, rest: &[Value]) -> BuiltinResult<Self> {
617 match rest {
618 [] => Ok(Self::Values(default_time_vector(system))),
619 [value] => {
620 let gathered = crate::dispatcher::gather_if_needed_async(value).await?;
621 if let Some(final_time) = scalar_time_from_value(&gathered)? {
622 return Ok(Self::Values(time_vector_from_final_time(
623 system, final_time,
624 )?));
625 }
626 let vector = time_vector_from_value(gathered)?;
627 validate_time_vector(system, &vector)?;
628 Ok(Self::Values(vector))
629 }
630 _ => Err(impulse_error_with_detail(
631 &IMPULSE_ERROR_INVALID_ARGUMENT,
632 "expected impulse(sys), impulse(sys, tFinal), or impulse(sys, t)",
633 )),
634 }
635 }
636}
637
638fn scalar_time_from_value(value: &Value) -> BuiltinResult<Option<f64>> {
639 match value {
640 Value::Num(n) => Ok(Some(*n)),
641 Value::Int(i) => Ok(Some(exact_integer_as_f64(i, "time")?)),
642 Value::Bool(b) => Ok(Some(if *b { 1.0 } else { 0.0 })),
643 Value::Tensor(tensor) if tensor::is_scalar_tensor(tensor) => {
644 ensure_exact_integer_tensor(tensor, "time")?;
645 Ok(Some(tensor::tensor_value_f64(tensor, 0)))
646 }
647 Value::LogicalArray(logical) if logical.data.len() == 1 => {
648 Ok(Some(if logical.data[0] == 0 { 0.0 } else { 1.0 }))
649 }
650 Value::Tensor(_) | Value::LogicalArray(_) => Ok(None),
651 _ => Ok(None),
652 }
653}
654
655fn default_time_vector(system: &TfSystem) -> Vec<f64> {
656 if system.is_discrete() {
657 (0..DEFAULT_POINTS)
658 .map(|idx| idx as f64 * system.sample_time)
659 .collect()
660 } else {
661 linspace(0.0, 10.0, DEFAULT_POINTS)
662 }
663}
664
665fn time_vector_from_final_time(system: &TfSystem, final_time: f64) -> BuiltinResult<Vec<f64>> {
666 if !final_time.is_finite() || final_time < 0.0 {
667 return Err(impulse_error_with_detail(
668 &IMPULSE_ERROR_INVALID_TIME,
669 "final time must be a finite non-negative scalar",
670 ));
671 }
672 if system.is_discrete() {
673 let count = checked_discrete_sample_count(system, final_time)?;
674 Ok((0..count)
675 .map(|idx| idx as f64 * system.sample_time)
676 .collect())
677 } else if final_time == 0.0 {
678 Ok(vec![0.0])
679 } else {
680 Ok(linspace(0.0, final_time, DEFAULT_POINTS))
681 }
682}
683
684fn checked_discrete_sample_count(system: &TfSystem, final_time: f64) -> BuiltinResult<usize> {
685 let samples = final_time / system.sample_time;
686 if !samples.is_finite() {
687 return Err(impulse_error_with_detail(
688 &IMPULSE_ERROR_DISCRETE_LIMIT,
689 "discrete sample count exceeds platform limits",
690 ));
691 }
692
693 let count = samples.floor() + 1.0;
694 if count > usize::MAX as f64 || count > MAX_DISCRETE_SAMPLES as f64 {
695 return Err(impulse_error_with_detail(
696 &IMPULSE_ERROR_DISCRETE_LIMIT,
697 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
698 ));
699 }
700 Ok(count as usize)
701}
702
703fn checked_discrete_sample_index(system: &TfSystem, time: f64) -> BuiltinResult<usize> {
704 let samples = time / system.sample_time;
705 let index = samples.round();
706 if !index.is_finite() || index > usize::MAX as f64 {
707 return Err(impulse_error_with_detail(
708 &IMPULSE_ERROR_DISCRETE_LIMIT,
709 "discrete sample index exceeds platform limits",
710 ));
711 }
712 if index >= MAX_DISCRETE_SAMPLES as f64 {
713 return Err(impulse_error_with_detail(
714 &IMPULSE_ERROR_DISCRETE_LIMIT,
715 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
716 ));
717 }
718 Ok(index as usize)
719}
720
721fn linspace(start: f64, stop: f64, count: usize) -> Vec<f64> {
722 if count <= 1 {
723 return vec![start];
724 }
725 let step = (stop - start) / ((count - 1) as f64);
726 (0..count).map(|idx| start + idx as f64 * step).collect()
727}
728
729fn time_vector_from_value(value: Value) -> BuiltinResult<Vec<f64>> {
730 let tensor = tensor::value_into_tensor_for(BUILTIN_NAME, value).map_err(|err| {
731 impulse_error_with_detail(
732 &IMPULSE_ERROR_INVALID_TIME,
733 format!("time vector must be numeric: {err}"),
734 )
735 })?;
736 ensure_vector("time", &tensor.shape)?;
737 ensure_exact_integer_tensor(&tensor, "time")?;
738 let values = tensor::tensor_into_values_f64(tensor);
739 if values.is_empty() {
740 return Err(impulse_error_with_detail(
741 &IMPULSE_ERROR_INVALID_TIME,
742 "time vector cannot be empty",
743 ));
744 }
745 Ok(values)
746}
747
748fn validate_time_vector(system: &TfSystem, values: &[f64]) -> BuiltinResult<()> {
749 if values
750 .iter()
751 .any(|value| !value.is_finite() || *value < 0.0)
752 {
753 return Err(impulse_error_with_detail(
754 &IMPULSE_ERROR_INVALID_TIME,
755 "time vector values must be finite and non-negative",
756 ));
757 }
758 if values.windows(2).any(|pair| pair[1] <= pair[0]) {
759 return Err(impulse_error_with_detail(
760 &IMPULSE_ERROR_INVALID_TIME,
761 "time vector values must be strictly increasing",
762 ));
763 }
764 if system.is_discrete() {
765 for &value in values {
766 let samples = value / system.sample_time;
767 if (samples - samples.round()).abs() > 1.0e-8 {
768 return Err(impulse_error_with_detail(
769 &IMPULSE_ERROR_INVALID_TIME,
770 "discrete-time vectors must use integer multiples of the sample time",
771 ));
772 }
773 }
774 }
775 Ok(())
776}
777
778#[derive(Clone, Debug)]
779struct ImpulseResponse {
780 t: Vec<f64>,
781 y: Vec<f64>,
782 #[cfg(feature = "plot-core")]
783 discrete: bool,
784}
785
786impl ImpulseResponse {
787 fn y_value(&self) -> BuiltinResult<Value> {
788 let tensor = Tensor::new(self.y.clone(), vec![self.y.len(), 1]).map_err(|err| {
789 impulse_error_with_detail(
790 &IMPULSE_ERROR_INTERNAL,
791 format!("failed to build response tensor: {err}"),
792 )
793 })?;
794 Ok(Value::Tensor(tensor))
795 }
796
797 fn t_value(&self) -> BuiltinResult<Value> {
798 let tensor = Tensor::new(self.t.clone(), vec![self.t.len(), 1]).map_err(|err| {
799 impulse_error_with_detail(
800 &IMPULSE_ERROR_INTERNAL,
801 format!("failed to build time tensor: {err}"),
802 )
803 })?;
804 Ok(Value::Tensor(tensor))
805 }
806}
807
808fn evaluate_impulse(system: &TfSystem, time: TimeSpec) -> BuiltinResult<ImpulseResponse> {
809 let TimeSpec::Values(t) = time;
810 let realization = Realization::from_tf(system)?;
811 let y = if system.is_discrete() {
812 discrete_response(system, &realization, &t)?
813 } else {
814 continuous_response(&realization, &t)
815 };
816 Ok(ImpulseResponse {
817 t,
818 y,
819 #[cfg(feature = "plot-core")]
820 discrete: system.is_discrete(),
821 })
822}
823
824#[derive(Clone, Debug)]
825struct Realization {
826 a: DMatrix<f64>,
827 c: Vec<f64>,
828}
829
830impl Realization {
831 fn from_tf(system: &TfSystem) -> BuiltinResult<Self> {
832 if system.numerator.is_empty() {
833 let order = system.denominator.len().saturating_sub(1).max(1);
834 return Ok(Self {
835 a: DMatrix::zeros(order, order),
836 c: vec![0.0; order],
837 });
838 }
839 let leading = system.denominator[0];
840 if leading.abs() <= EPS {
841 return Err(impulse_error_with_detail(
842 &IMPULSE_ERROR_INVALID_MODEL,
843 "denominator leading coefficient must be non-zero",
844 ));
845 }
846 let denominator: Vec<f64> = system
847 .denominator
848 .iter()
849 .map(|value| *value / leading)
850 .collect();
851 let mut numerator: Vec<f64> = system
852 .numerator
853 .iter()
854 .map(|value| *value / leading)
855 .collect();
856 let order = denominator.len() - 1;
857 while numerator.len() < order {
858 numerator.insert(0, 0.0);
859 }
860
861 let mut a = DMatrix::<f64>::zeros(order, order);
862 for row in 0..order.saturating_sub(1) {
863 a[(row, row + 1)] = 1.0;
864 }
865 for col in 0..order {
866 a[(order - 1, col)] = -denominator[order - col];
867 }
868 let c = numerator.into_iter().rev().collect();
869 Ok(Self { a, c })
870 }
871}
872
873fn continuous_response(realization: &Realization, t: &[f64]) -> Vec<f64> {
874 t.iter()
875 .map(|&time| {
876 let exp_at = matrix_exp(&(realization.a.clone() * time));
877 dot_c_with_last_column(&realization.c, &exp_at)
878 })
879 .collect()
880}
881
882fn discrete_response(
883 system: &TfSystem,
884 realization: &Realization,
885 t: &[f64],
886) -> BuiltinResult<Vec<f64>> {
887 if t.len() > MAX_DISCRETE_SAMPLES {
888 return Err(impulse_error_with_detail(
889 &IMPULSE_ERROR_DISCRETE_LIMIT,
890 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
891 ));
892 }
893 let sample_indices: Vec<usize> = t
894 .iter()
895 .map(|value| checked_discrete_sample_index(system, *value))
896 .collect::<BuiltinResult<_>>()?;
897 let max_index = sample_indices.iter().copied().max().unwrap_or(0);
898 let order = realization.c.len();
899 let value_count = max_index.checked_add(1).ok_or_else(|| {
900 impulse_error_with_detail(
901 &IMPULSE_ERROR_DISCRETE_LIMIT,
902 "discrete sample index exceeds platform limits",
903 )
904 })?;
905 let mut values = vec![0.0; value_count];
906 if order == 0 {
907 return Ok(sample_indices.into_iter().map(|idx| values[idx]).collect());
908 }
909
910 let mut state = vec![0.0; order];
911 state[order - 1] = 1.0;
912 let impulse_scale = 1.0 / system.sample_time;
913 for k in 1..=max_index {
914 values[k] = dot(&realization.c, &state) * impulse_scale;
915 state = mat_vec_mul(&realization.a, &state);
916 }
917 Ok(sample_indices.into_iter().map(|idx| values[idx]).collect())
918}
919
920fn dot_c_with_last_column(c: &[f64], matrix: &DMatrix<f64>) -> f64 {
921 if c.is_empty() {
922 return 0.0;
923 }
924 let last_col = matrix.ncols() - 1;
925 c.iter()
926 .enumerate()
927 .map(|(row, coeff)| coeff * matrix[(row, last_col)])
928 .sum()
929}
930
931fn dot(lhs: &[f64], rhs: &[f64]) -> f64 {
932 lhs.iter().zip(rhs).map(|(a, b)| a * b).sum()
933}
934
935fn mat_vec_mul(matrix: &DMatrix<f64>, vector: &[f64]) -> Vec<f64> {
936 let mut out = vec![0.0; matrix.nrows()];
937 for row in 0..matrix.nrows() {
938 let mut acc = 0.0;
939 for col in 0..matrix.ncols() {
940 acc += matrix[(row, col)] * vector[col];
941 }
942 out[row] = acc;
943 }
944 out
945}
946
947fn matrix_exp(matrix: &DMatrix<f64>) -> DMatrix<f64> {
948 let norm = matrix_one_norm(matrix);
949 let scale_power = if norm <= 0.5 {
950 0usize
951 } else {
952 norm.log2().ceil().max(0.0) as usize + 1
953 };
954 let scale = 2.0_f64.powi(scale_power as i32);
955 let scaled = matrix / scale;
956 let n = matrix.nrows();
957 let mut result = DMatrix::<f64>::identity(n, n);
958 let mut term = DMatrix::<f64>::identity(n, n);
959 for k in 1..=48 {
960 term = (&term * &scaled) / (k as f64);
961 result += &term;
962 if matrix_one_norm(&term) <= 1.0e-14 {
963 break;
964 }
965 }
966 for _ in 0..scale_power {
967 result = &result * &result;
968 }
969 result
970}
971
972fn matrix_one_norm(matrix: &DMatrix<f64>) -> f64 {
973 let mut best = 0.0;
974 for col in 0..matrix.ncols() {
975 let mut sum = 0.0;
976 for row in 0..matrix.nrows() {
977 sum += matrix[(row, col)].abs();
978 }
979 if sum > best {
980 best = sum;
981 }
982 }
983 best
984}
985
986#[cfg(feature = "plot-core")]
987async fn render_impulse_plot(response: &ImpulseResponse) -> BuiltinResult<Value> {
988 let t = response.t_value()?;
989 let y = response.y_value()?;
990 let plot_name = if response.discrete { "stem" } else { "plot" };
991 crate::dispatcher::call_builtin_async(plot_name, &[t, y]).await
992}
993
994#[cfg(not(feature = "plot-core"))]
995async fn render_impulse_plot(_response: &ImpulseResponse) -> BuiltinResult<Value> {
996 Ok(Value::Num(f64::NAN))
997}
998
999#[cfg(test)]
1000mod tests {
1001 use super::*;
1002 use futures::executor::block_on;
1003 use runmat_value::{CharArray, IntegerStorage, ObjectInstance};
1004
1005 fn tf_object(num: Vec<f64>, den: Vec<f64>, ts: f64) -> Value {
1006 tf_object_with_delays(num, den, ts, 0.0, 0.0)
1007 }
1008
1009 fn tf_object_with_delays(
1010 num: Vec<f64>,
1011 den: Vec<f64>,
1012 ts: f64,
1013 input_delay: f64,
1014 output_delay: f64,
1015 ) -> Value {
1016 let mut object = ObjectInstance::new("tf".to_string());
1017 object.properties.insert(
1018 "Numerator".to_string(),
1019 Value::Tensor(Tensor::new(num.clone(), vec![1, num.len()]).unwrap()),
1020 );
1021 object.properties.insert(
1022 "Denominator".to_string(),
1023 Value::Tensor(Tensor::new(den.clone(), vec![1, den.len()]).unwrap()),
1024 );
1025 object.properties.insert(
1026 "Variable".to_string(),
1027 Value::CharArray(CharArray::new_row(if ts > 0.0 { "z" } else { "s" })),
1028 );
1029 object.properties.insert("Ts".to_string(), Value::Num(ts));
1030 object
1031 .properties
1032 .insert("InputDelay".to_string(), Value::Num(input_delay));
1033 object
1034 .properties
1035 .insert("OutputDelay".to_string(), Value::Num(output_delay));
1036 Value::Object(object)
1037 }
1038
1039 fn run_impulse(system: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
1040 block_on(impulse_builtin(system, rest))
1041 }
1042
1043 fn tensor_data(value: Value) -> Vec<f64> {
1044 match value {
1045 Value::Tensor(tensor) => tensor.materialize_f64(),
1046 other => panic!("expected tensor, got {other:?}"),
1047 }
1048 }
1049
1050 fn integer_tensor(storage: IntegerStorage, shape: Vec<usize>) -> Value {
1051 let tensor = Tensor::new_integer(storage, shape).expect("integer tensor");
1052 Value::Tensor(tensor)
1053 }
1054
1055 #[test]
1056 fn impulse_descriptor_signatures_cover_core_forms() {
1057 let labels: Vec<&str> = IMPULSE_DESCRIPTOR
1058 .signatures
1059 .iter()
1060 .map(|sig| sig.label)
1061 .collect();
1062 assert!(labels.contains(&"y = impulse(sys)"));
1063 assert!(labels.contains(&"y = impulse(sys, tFinal)"));
1064 assert!(labels.contains(&"y = impulse(sys, t)"));
1065 assert!(labels.contains(&"[y,t] = impulse(sys)"));
1066 assert!(labels.contains(&"[y,t] = impulse(sys, tFinal)"));
1067 assert!(labels.contains(&"[y,t] = impulse(sys, t)"));
1068 }
1069
1070 #[test]
1071 fn impulse_first_order_continuous_explicit_time() {
1072 let sys = tf_object(vec![20.0], vec![1.0, 5.0], 0.0);
1073 let t = Value::Tensor(Tensor::new(vec![0.0, 0.1, 0.2], vec![1, 3]).unwrap());
1074 let y = tensor_data(run_impulse(sys, vec![t]).expect("impulse"));
1075 let expected = [20.0, 20.0 * (-0.5f64).exp(), 20.0 * (-1.0f64).exp()];
1076 for (actual, expected) in y.iter().zip(expected) {
1077 assert!((actual - expected).abs() < 1.0e-8);
1078 }
1079 }
1080
1081 #[test]
1082 fn impulse_second_order_continuous() {
1083 let sys = tf_object(vec![1.0], vec![1.0, 3.0, 2.0], 0.0);
1084 let t = Value::Tensor(Tensor::new(vec![0.0, 0.5, 1.0], vec![1, 3]).unwrap());
1085 let y = tensor_data(run_impulse(sys, vec![t]).expect("impulse"));
1086 for (actual, time) in y.iter().zip([0.0_f64, 0.5, 1.0]) {
1087 let expected = (-time).exp() - (-2.0 * time).exp();
1088 assert!((actual - expected).abs() < 1.0e-8);
1089 }
1090 }
1091
1092 #[test]
1093 fn impulse_multi_output_returns_y_and_t_columns() {
1094 let _guard = crate::output_count::push_output_count(Some(2));
1095 let sys = tf_object(vec![20.0], vec![1.0, 5.0], 0.0);
1096 let t_arg = Value::Tensor(Tensor::new(vec![0.0, 0.1], vec![1, 2]).unwrap());
1097 let result = run_impulse(sys, vec![t_arg]).expect("impulse");
1098 let Value::OutputList(outputs) = result else {
1099 panic!("expected output list");
1100 };
1101 assert_eq!(outputs.len(), 2);
1102 match &outputs[0] {
1103 Value::Tensor(tensor) => assert_eq!(tensor.shape, vec![2, 1]),
1104 other => panic!("expected y tensor, got {other:?}"),
1105 }
1106 match &outputs[1] {
1107 Value::Tensor(tensor) => {
1108 assert_eq!(tensor.shape, vec![2, 1]);
1109 assert_eq!(tensor.materialize_f64(), vec![0.0, 0.1]);
1110 }
1111 other => panic!("expected t tensor, got {other:?}"),
1112 }
1113 }
1114
1115 #[test]
1116 fn impulse_zero_output_count_emits_no_values() {
1117 let _guard = crate::output_count::push_output_count(Some(0));
1118 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
1119 let result = run_impulse(sys, Vec::new()).expect("impulse");
1120 let Value::OutputList(outputs) = result else {
1121 panic!("expected output list");
1122 };
1123 assert!(outputs.is_empty());
1124 }
1125
1126 #[test]
1127 fn impulse_requested_zero_outputs_emits_no_values() {
1128 let _guard = crate::output_context::push_output_count(0);
1129 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
1130 let result = run_impulse(sys, Vec::new()).expect("impulse");
1131 let Value::OutputList(outputs) = result else {
1132 panic!("expected output list");
1133 };
1134 assert!(outputs.is_empty());
1135 }
1136
1137 #[test]
1138 fn impulse_discrete_siso_response() {
1139 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 0.1);
1140 let t = Value::Tensor(Tensor::new(vec![0.0, 0.1, 0.2, 0.3], vec![1, 4]).unwrap());
1141 let y = tensor_data(run_impulse(sys, vec![t]).expect("impulse"));
1142 assert_eq!(y.len(), 4);
1143 assert!((y[0] - 0.0).abs() < 1.0e-12);
1144 assert!((y[1] - 10.0).abs() < 1.0e-12);
1145 assert!((y[2] - 5.0).abs() < 1.0e-12);
1146 assert!((y[3] - 2.5).abs() < 1.0e-12);
1147 }
1148
1149 #[test]
1150 fn impulse_typed_integer_coefficients_and_time_cross_double_boundary_exactly() {
1151 let _compat = crate::compatibility::push_runmat_extensions_enabled(true);
1152 let mut object = ObjectInstance::new("tf".to_string());
1153 object.properties.insert(
1154 "Numerator".to_string(),
1155 integer_tensor(IntegerStorage::U16(vec![20]), vec![1, 1]),
1156 );
1157 object.properties.insert(
1158 "Denominator".to_string(),
1159 integer_tensor(IntegerStorage::I16(vec![1, 5]), vec![1, 2]),
1160 );
1161 object.properties.insert(
1162 "Variable".to_string(),
1163 Value::CharArray(CharArray::new_row("s")),
1164 );
1165 object.properties.insert("Ts".to_string(), Value::Num(0.0));
1166 object
1167 .properties
1168 .insert("InputDelay".to_string(), Value::Num(0.0));
1169 object
1170 .properties
1171 .insert("OutputDelay".to_string(), Value::Num(0.0));
1172
1173 let t = integer_tensor(IntegerStorage::U64(vec![0, 1]), vec![1, 2]);
1174 let _guard = crate::output_count::push_output_count(Some(2));
1175 let result = run_impulse(Value::Object(object), vec![t]).expect("impulse");
1176 let Value::OutputList(outputs) = result else {
1177 panic!("expected output list");
1178 };
1179 assert_eq!(tensor_data(outputs[1].clone()), vec![0.0, 1.0]);
1180 }
1181
1182 #[test]
1183 fn impulse_time_vector_parser_ignores_poisoned_integer_mirrors_for_all_classes() {
1184 let storages = [
1185 IntegerStorage::I8(vec![0, 1]),
1186 IntegerStorage::I16(vec![0, 1]),
1187 IntegerStorage::I32(vec![0, 1]),
1188 IntegerStorage::I64(vec![0, 1]),
1189 IntegerStorage::U8(vec![0, 1]),
1190 IntegerStorage::U16(vec![0, 1]),
1191 IntegerStorage::U32(vec![0, 1]),
1192 IntegerStorage::U64(vec![0, 1]),
1193 ];
1194
1195 for storage in storages {
1196 assert_eq!(
1197 time_vector_from_value(integer_tensor(storage, vec![1, 2]))
1198 .expect("typed integer time vector"),
1199 vec![0.0, 1.0]
1200 );
1201 }
1202 }
1203
1204 #[test]
1205 fn impulse_scalar_final_time_reads_typed_integer_storage_length_exactly() {
1206 let _compat = crate::compatibility::push_runmat_extensions_enabled(true);
1207 let mut object = ObjectInstance::new("tf".to_string());
1208 object.properties.insert(
1209 "Numerator".to_string(),
1210 Value::Tensor(Tensor::new(vec![1.0], vec![1, 1]).unwrap()),
1211 );
1212 object.properties.insert(
1213 "Denominator".to_string(),
1214 Value::Tensor(Tensor::new(vec![1.0, 1.0], vec![1, 2]).unwrap()),
1215 );
1216 object.properties.insert(
1217 "Variable".to_string(),
1218 Value::CharArray(CharArray::new_row("s")),
1219 );
1220 object.properties.insert("Ts".to_string(), Value::Num(0.0));
1221 object
1222 .properties
1223 .insert("InputDelay".to_string(), Value::Num(0.0));
1224 object
1225 .properties
1226 .insert("OutputDelay".to_string(), Value::Num(0.0));
1227
1228 let final_time =
1229 Tensor::new_integer(IntegerStorage::U16(vec![2]), vec![1, 1]).expect("final time");
1230 let _guard = crate::output_count::push_output_count(Some(2));
1231 let result =
1232 run_impulse(Value::Object(object), vec![Value::Tensor(final_time)]).expect("impulse");
1233 let Value::OutputList(outputs) = result else {
1234 panic!("expected output list");
1235 };
1236 let time = tensor_data(outputs[1].clone());
1237 assert_eq!(time.first().copied(), Some(0.0));
1238 assert_eq!(time.last().copied(), Some(2.0));
1239 }
1240
1241 #[test]
1242 fn impulse_discrete_final_time_rejects_excessive_sample_count() {
1243 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0e-6);
1244 let err = run_impulse(sys, vec![Value::Num(2.0)]).expect_err("should fail");
1245 assert!(err.message().contains("more than 1000000 samples"));
1246 assert_eq!(err.identifier(), IMPULSE_ERROR_DISCRETE_LIMIT.identifier);
1247 }
1248
1249 #[test]
1250 fn impulse_discrete_time_vector_rejects_excessive_sample_index() {
1251 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0);
1252 let t =
1253 Value::Tensor(Tensor::new(vec![0.0, MAX_DISCRETE_SAMPLES as f64], vec![1, 2]).unwrap());
1254 let err = run_impulse(sys, vec![t]).expect_err("should fail");
1255 assert!(err.message().contains("more than 1000000 samples"));
1256 }
1257
1258 #[test]
1259 fn impulse_rejects_unsupported_model_type() {
1260 let object = ObjectInstance::new("ss".to_string());
1261 let err = run_impulse(Value::Object(object), Vec::new()).expect_err("should fail");
1262 assert!(err.message().contains("unsupported model class"));
1263 assert_eq!(err.identifier(), IMPULSE_ERROR_UNSUPPORTED_MODEL.identifier);
1264 }
1265
1266 #[test]
1267 fn impulse_rejects_direct_feedthrough_tf() {
1268 let sys = tf_object(vec![1.0, 1.0], vec![1.0, 2.0], 0.0);
1269 let err = run_impulse(sys, Vec::new()).expect_err("should fail");
1270 assert!(err.message().contains("strictly proper"));
1271 }
1272
1273 #[test]
1274 fn impulse_rejects_invalid_time_metadata() {
1275 let err = run_impulse(tf_object(vec![1.0], vec![1.0, -0.5], -0.1), Vec::new())
1276 .expect_err("negative sample time should fail");
1277 assert!(err.message().contains("Ts must be"));
1278
1279 let err = run_impulse(
1280 tf_object_with_delays(vec![1.0], vec![1.0, 5.0], 0.0, f64::NAN, 0.0),
1281 Vec::new(),
1282 )
1283 .expect_err("NaN input delay should fail");
1284 assert!(err.message().contains("InputDelay must be"));
1285
1286 let err = run_impulse(
1287 tf_object_with_delays(vec![1.0], vec![1.0, 5.0], 0.0, 0.0, -1.0),
1288 Vec::new(),
1289 )
1290 .expect_err("negative output delay should fail");
1291 assert!(err.message().contains("OutputDelay must be"));
1292 }
1293
1294 #[test]
1295 fn impulse_runmat_numeric_extensions_are_gated_and_exact() {
1296 {
1297 let _compat = crate::compatibility::push_runmat_extensions_enabled(false);
1298 let integer = run_impulse(
1299 tf_object(vec![1.0], vec![1.0, 1.0], 0.0),
1300 vec![Value::Int(IntValue::U16(2))],
1301 )
1302 .unwrap_err();
1303 assert_eq!(
1304 integer.identifier(),
1305 Some("RunMat:compatibility:ImpulseIntegerExtension")
1306 );
1307 let single = run_impulse(
1308 tf_object(vec![1.0], vec![1.0, 1.0], 0.0),
1309 vec![Value::Tensor(
1310 Tensor::from_f32(vec![0.0, 1.0], vec![1, 2]).unwrap(),
1311 )],
1312 )
1313 .unwrap_err();
1314 assert_eq!(
1315 single.identifier(),
1316 Some("RunMat:compatibility:ImpulseSingleExtension")
1317 );
1318 }
1319
1320 let _compat = crate::compatibility::push_runmat_extensions_enabled(true);
1321 let wide = run_impulse(
1322 tf_object(vec![1.0], vec![1.0, 1.0], 0.0),
1323 vec![Value::Int(IntValue::U64((1_u64 << 53) + 1))],
1324 )
1325 .unwrap_err();
1326 assert!(wide.message().contains("exactly representable as double"));
1327 }
1328
1329 #[test]
1330 fn impulse_resident_time_is_an_explicit_host_output_extension() {
1331 crate::builtins::common::test_support::with_test_provider(|provider| {
1332 let tensor = Tensor::new(vec![0.0, 1.0], vec![1, 2]).unwrap();
1333 let handle = crate::builtins::common::gpu_helpers::upload_tensor(provider, &tensor)
1334 .expect("resident time");
1335 runmat_accelerate::fusion_residency::mark(&handle);
1336 {
1337 let _compat = crate::compatibility::push_runmat_extensions_enabled(false);
1338 let error = run_impulse(
1339 tf_object(vec![1.0], vec![1.0, 1.0], 0.0),
1340 vec![Value::GpuTensor(handle.clone())],
1341 )
1342 .unwrap_err();
1343 assert_eq!(
1344 error.identifier(),
1345 Some("RunMat:compatibility:ImpulseResidentExtension")
1346 );
1347 }
1348 let _compat = crate::compatibility::push_runmat_extensions_enabled(true);
1349 let output = run_impulse(
1350 tf_object(vec![1.0], vec![1.0, 1.0], 0.0),
1351 vec![Value::GpuTensor(handle.clone())],
1352 )
1353 .expect("resident host fallback");
1354 assert!(matches!(output, Value::Tensor(_)));
1355 assert!(runmat_accelerate::fusion_residency::is_resident(&handle));
1356 let _ = provider.free(&handle);
1357 });
1358 }
1359}