runmat_runtime/builtins/control/
impulse.rs

1//! MATLAB-compatible `impulse` response builtin for supported control models.
2
3use nalgebra::DMatrix;
4use runmat_builtins::{Tensor, Value};
5use runmat_macros::runtime_builtin;
6
7use crate::builtins::common::spec::{
8    BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
9    ReductionNaN, ResidencyPolicy, ShapeRequirements,
10};
11use crate::builtins::common::tensor;
12use crate::builtins::control::type_resolvers::impulse_type;
13use crate::{build_runtime_error, BuiltinResult, RuntimeError};
14
15const BUILTIN_NAME: &str = "impulse";
16const TF_CLASS: &str = "tf";
17const EPS: f64 = 1.0e-12;
18const DEFAULT_POINTS: usize = 100;
19const MAX_DISCRETE_SAMPLES: usize = 1_000_000;
20
21#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::control::impulse")]
22pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
23    name: "impulse",
24    op_kind: GpuOpKind::Custom("control-impulse-response"),
25    supported_precisions: &[],
26    broadcast: BroadcastSemantics::None,
27    provider_hooks: &[],
28    constant_strategy: ConstantStrategy::InlineLiteral,
29    residency: ResidencyPolicy::GatherImmediately,
30    nan_mode: ReductionNaN::Include,
31    two_pass_threshold: None,
32    workgroup_size: None,
33    accepts_nan_mode: false,
34    notes: "Control-system response evaluation runs on the host. GPU-resident metadata is gathered before simulation.",
35};
36
37#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::control::impulse")]
38pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
39    name: "impulse",
40    shape: ShapeRequirements::Any,
41    constant_strategy: ConstantStrategy::InlineLiteral,
42    elementwise: None,
43    reduction: None,
44    emits_nan: false,
45    notes: "Impulse-response simulation materialises host-side time and output vectors and terminates fusion chains.",
46};
47
48fn impulse_error(message: impl Into<String>) -> RuntimeError {
49    build_runtime_error(message)
50        .with_builtin(BUILTIN_NAME)
51        .build()
52}
53
54#[runtime_builtin(
55    name = "impulse",
56    category = "control",
57    summary = "Compute or plot the impulse response of a supported dynamic system model.",
58    keywords = "impulse,control system,transfer function,response,tf",
59    type_resolver(impulse_type),
60    builtin_path = "crate::builtins::control::impulse"
61)]
62async fn impulse_builtin(system: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
63    let system = TfSystem::parse(system).await?;
64    let time = TimeSpec::parse(&system, &rest).await?;
65    let response = evaluate_impulse(&system, time)?;
66
67    if let Some(out_count) = crate::output_count::current_output_count() {
68        if out_count == 0 {
69            return Ok(Value::OutputList(Vec::new()));
70        }
71        if out_count == 1 {
72            return Ok(Value::OutputList(vec![response.y_value()?]));
73        }
74        return Ok(crate::output_count::output_list_with_padding(
75            out_count,
76            vec![response.y_value()?, response.t_value()?],
77        ));
78    }
79
80    if crate::output_context::requested_output_count() == Some(0) {
81        return render_impulse_plot(&response).await;
82    }
83
84    response.y_value()
85}
86
87#[derive(Clone, Debug)]
88struct TfSystem {
89    numerator: Vec<f64>,
90    denominator: Vec<f64>,
91    sample_time: f64,
92}
93
94impl TfSystem {
95    async fn parse(value: Value) -> BuiltinResult<Self> {
96        let gathered = crate::dispatcher::gather_if_needed_async(&value).await?;
97        let Value::Object(object) = gathered else {
98            return Err(impulse_error(format!(
99                "impulse: expected a dynamic system model, got {gathered:?}"
100            )));
101        };
102        if object.class_name != TF_CLASS {
103            return Err(impulse_error(format!(
104                "impulse: unsupported model class '{}'; only SISO tf objects are currently supported",
105                object.class_name
106            )));
107        }
108
109        let numerator = real_coefficients(property(&object, "Numerator")?, "Numerator")?;
110        let denominator = real_coefficients(property(&object, "Denominator")?, "Denominator")?;
111        let sample_time = scalar_property(property(&object, "Ts")?, "Ts")?;
112        let input_delay = scalar_property(property(&object, "InputDelay")?, "InputDelay")?;
113        let output_delay = scalar_property(property(&object, "OutputDelay")?, "OutputDelay")?;
114        if !sample_time.is_finite() || sample_time < 0.0 {
115            return Err(impulse_error(format!(
116                "impulse: Ts must be a finite non-negative scalar, got {sample_time}"
117            )));
118        }
119        if !input_delay.is_finite() || input_delay < 0.0 {
120            return Err(impulse_error(format!(
121                "impulse: InputDelay must be a finite non-negative scalar, got {input_delay}"
122            )));
123        }
124        if !output_delay.is_finite() || output_delay < 0.0 {
125            return Err(impulse_error(format!(
126                "impulse: OutputDelay must be a finite non-negative scalar, got {output_delay}"
127            )));
128        }
129        if input_delay.abs() > EPS || output_delay.abs() > EPS {
130            return Err(impulse_error(
131                "impulse: transfer functions with input or output delays are not supported yet",
132            ));
133        }
134
135        let numerator = trim_leading_zeros(numerator);
136        let denominator = trim_leading_zeros(denominator);
137        if denominator.is_empty() {
138            return Err(impulse_error(
139                "impulse: denominator coefficients cannot be empty",
140            ));
141        }
142        if numerator.is_empty() {
143            return Ok(Self {
144                numerator,
145                denominator,
146                sample_time,
147            });
148        }
149        if denominator.len() <= 1 {
150            return Err(impulse_error(
151                "impulse: static-gain transfer functions do not have a finite impulse-response vector",
152            ));
153        }
154        if numerator.len() >= denominator.len() {
155            return Err(impulse_error(
156                "impulse: only strictly proper SISO tf models are currently supported",
157            ));
158        }
159
160        Ok(Self {
161            numerator,
162            denominator,
163            sample_time,
164        })
165    }
166
167    fn is_discrete(&self) -> bool {
168        self.sample_time > 0.0
169    }
170}
171
172fn property<'a>(
173    object: &'a runmat_builtins::ObjectInstance,
174    name: &str,
175) -> BuiltinResult<&'a Value> {
176    object
177        .properties
178        .get(name)
179        .ok_or_else(|| impulse_error(format!("impulse: tf object is missing {name} property")))
180}
181
182fn real_coefficients(value: &Value, label: &str) -> BuiltinResult<Vec<f64>> {
183    match value {
184        Value::Tensor(tensor) => {
185            ensure_vector(label, &tensor.shape)?;
186            finite_values(label, tensor.data.clone())
187        }
188        Value::Num(n) => finite_values(label, vec![*n]),
189        Value::Int(i) => finite_values(label, vec![i.to_f64()]),
190        Value::Bool(b) => finite_values(label, vec![if *b { 1.0 } else { 0.0 }]),
191        Value::LogicalArray(logical) => {
192            ensure_vector(label, &logical.shape)?;
193            finite_values(
194                label,
195                logical
196                    .data
197                    .iter()
198                    .map(|bit| if *bit == 0 { 0.0 } else { 1.0 })
199                    .collect(),
200            )
201        }
202        Value::Complex(_, _) | Value::ComplexTensor(_) => Err(impulse_error(
203            "impulse: complex-coefficient tf models are not supported yet",
204        )),
205        other => Err(impulse_error(format!(
206            "impulse: {label} must be a real numeric coefficient vector, got {other:?}"
207        ))),
208    }
209}
210
211fn finite_values(label: &str, values: Vec<f64>) -> BuiltinResult<Vec<f64>> {
212    if values.iter().any(|value| !value.is_finite()) {
213        return Err(impulse_error(format!(
214            "impulse: {label} coefficients must be finite"
215        )));
216    }
217    Ok(values)
218}
219
220fn ensure_vector(label: &str, shape: &[usize]) -> BuiltinResult<()> {
221    let non_unit = shape.iter().copied().filter(|&dim| dim > 1).count();
222    if non_unit <= 1 {
223        Ok(())
224    } else {
225        Err(impulse_error(format!(
226            "impulse: {label} coefficients must be a vector"
227        )))
228    }
229}
230
231fn scalar_property(value: &Value, label: &str) -> BuiltinResult<f64> {
232    match value {
233        Value::Num(n) => Ok(*n),
234        Value::Int(i) => Ok(i.to_f64()),
235        Value::Bool(b) => Ok(if *b { 1.0 } else { 0.0 }),
236        Value::Tensor(tensor) if tensor.data.len() == 1 => Ok(tensor.data[0]),
237        other => Err(impulse_error(format!(
238            "impulse: {label} must be a real scalar, got {other:?}"
239        ))),
240    }
241}
242
243fn trim_leading_zeros(values: Vec<f64>) -> Vec<f64> {
244    let first = values.iter().position(|value| value.abs() > EPS);
245    match first {
246        Some(idx) => values[idx..].to_vec(),
247        None => Vec::new(),
248    }
249}
250
251#[derive(Clone, Debug)]
252enum TimeSpec {
253    Values(Vec<f64>),
254}
255
256impl TimeSpec {
257    async fn parse(system: &TfSystem, rest: &[Value]) -> BuiltinResult<Self> {
258        match rest {
259            [] => Ok(Self::Values(default_time_vector(system))),
260            [value] => {
261                let gathered = crate::dispatcher::gather_if_needed_async(value).await?;
262                if let Some(final_time) = scalar_time_from_value(&gathered)? {
263                    return Ok(Self::Values(time_vector_from_final_time(
264                        system, final_time,
265                    )?));
266                }
267                let vector = time_vector_from_value(gathered)?;
268                validate_time_vector(system, &vector)?;
269                Ok(Self::Values(vector))
270            }
271            _ => Err(impulse_error(
272                "impulse: expected impulse(sys), impulse(sys, tFinal), or impulse(sys, t)",
273            )),
274        }
275    }
276}
277
278fn scalar_time_from_value(value: &Value) -> BuiltinResult<Option<f64>> {
279    match value {
280        Value::Num(n) => Ok(Some(*n)),
281        Value::Int(i) => Ok(Some(i.to_f64())),
282        Value::Bool(b) => Ok(Some(if *b { 1.0 } else { 0.0 })),
283        Value::Tensor(tensor) if tensor.data.len() == 1 => Ok(Some(tensor.data[0])),
284        Value::LogicalArray(logical) if logical.data.len() == 1 => {
285            Ok(Some(if logical.data[0] == 0 { 0.0 } else { 1.0 }))
286        }
287        Value::Tensor(_) | Value::LogicalArray(_) => Ok(None),
288        _ => Ok(None),
289    }
290}
291
292fn default_time_vector(system: &TfSystem) -> Vec<f64> {
293    if system.is_discrete() {
294        (0..DEFAULT_POINTS)
295            .map(|idx| idx as f64 * system.sample_time)
296            .collect()
297    } else {
298        linspace(0.0, 10.0, DEFAULT_POINTS)
299    }
300}
301
302fn time_vector_from_final_time(system: &TfSystem, final_time: f64) -> BuiltinResult<Vec<f64>> {
303    if !final_time.is_finite() || final_time < 0.0 {
304        return Err(impulse_error(
305            "impulse: final time must be a finite non-negative scalar",
306        ));
307    }
308    if system.is_discrete() {
309        let count = checked_discrete_sample_count(system, final_time)?;
310        Ok((0..count)
311            .map(|idx| idx as f64 * system.sample_time)
312            .collect())
313    } else if final_time == 0.0 {
314        Ok(vec![0.0])
315    } else {
316        Ok(linspace(0.0, final_time, DEFAULT_POINTS))
317    }
318}
319
320fn checked_discrete_sample_count(system: &TfSystem, final_time: f64) -> BuiltinResult<usize> {
321    let samples = final_time / system.sample_time;
322    if !samples.is_finite() {
323        return Err(impulse_error(
324            "impulse: discrete sample count exceeds platform limits",
325        ));
326    }
327
328    let count = samples.floor() + 1.0;
329    if count > usize::MAX as f64 || count > MAX_DISCRETE_SAMPLES as f64 {
330        return Err(impulse_error(format!(
331            "impulse: discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"
332        )));
333    }
334    Ok(count as usize)
335}
336
337fn checked_discrete_sample_index(system: &TfSystem, time: f64) -> BuiltinResult<usize> {
338    let samples = time / system.sample_time;
339    let index = samples.round();
340    if !index.is_finite() || index > usize::MAX as f64 {
341        return Err(impulse_error(
342            "impulse: discrete sample index exceeds platform limits",
343        ));
344    }
345    if index >= MAX_DISCRETE_SAMPLES as f64 {
346        return Err(impulse_error(format!(
347            "impulse: discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"
348        )));
349    }
350    Ok(index as usize)
351}
352
353fn linspace(start: f64, stop: f64, count: usize) -> Vec<f64> {
354    if count <= 1 {
355        return vec![start];
356    }
357    let step = (stop - start) / ((count - 1) as f64);
358    (0..count).map(|idx| start + idx as f64 * step).collect()
359}
360
361fn time_vector_from_value(value: Value) -> BuiltinResult<Vec<f64>> {
362    let tensor = tensor::value_into_tensor_for(BUILTIN_NAME, value)
363        .map_err(|err| impulse_error(format!("impulse: time vector must be numeric: {err}")))?;
364    ensure_vector("time", &tensor.shape)?;
365    if tensor.data.is_empty() {
366        return Err(impulse_error("impulse: time vector cannot be empty"));
367    }
368    Ok(tensor.data)
369}
370
371fn validate_time_vector(system: &TfSystem, values: &[f64]) -> BuiltinResult<()> {
372    if values
373        .iter()
374        .any(|value| !value.is_finite() || *value < 0.0)
375    {
376        return Err(impulse_error(
377            "impulse: time vector values must be finite and non-negative",
378        ));
379    }
380    if values.windows(2).any(|pair| pair[1] <= pair[0]) {
381        return Err(impulse_error(
382            "impulse: time vector values must be strictly increasing",
383        ));
384    }
385    if system.is_discrete() {
386        for &value in values {
387            let samples = value / system.sample_time;
388            if (samples - samples.round()).abs() > 1.0e-8 {
389                return Err(impulse_error(
390                    "impulse: discrete-time vectors must use integer multiples of the sample time",
391                ));
392            }
393        }
394    }
395    Ok(())
396}
397
398#[derive(Clone, Debug)]
399struct ImpulseResponse {
400    t: Vec<f64>,
401    y: Vec<f64>,
402    discrete: bool,
403}
404
405impl ImpulseResponse {
406    fn y_value(&self) -> BuiltinResult<Value> {
407        let tensor = Tensor::new(self.y.clone(), vec![self.y.len(), 1])
408            .map_err(|err| impulse_error(format!("impulse: {err}")))?;
409        Ok(Value::Tensor(tensor))
410    }
411
412    fn t_value(&self) -> BuiltinResult<Value> {
413        let tensor = Tensor::new(self.t.clone(), vec![self.t.len(), 1])
414            .map_err(|err| impulse_error(format!("impulse: {err}")))?;
415        Ok(Value::Tensor(tensor))
416    }
417}
418
419fn evaluate_impulse(system: &TfSystem, time: TimeSpec) -> BuiltinResult<ImpulseResponse> {
420    let TimeSpec::Values(t) = time;
421    let realization = Realization::from_tf(system)?;
422    let y = if system.is_discrete() {
423        discrete_response(system, &realization, &t)?
424    } else {
425        continuous_response(&realization, &t)
426    };
427    Ok(ImpulseResponse {
428        t,
429        y,
430        discrete: system.is_discrete(),
431    })
432}
433
434#[derive(Clone, Debug)]
435struct Realization {
436    a: DMatrix<f64>,
437    c: Vec<f64>,
438}
439
440impl Realization {
441    fn from_tf(system: &TfSystem) -> BuiltinResult<Self> {
442        if system.numerator.is_empty() {
443            let order = system.denominator.len().saturating_sub(1).max(1);
444            return Ok(Self {
445                a: DMatrix::zeros(order, order),
446                c: vec![0.0; order],
447            });
448        }
449        let leading = system.denominator[0];
450        if leading.abs() <= EPS {
451            return Err(impulse_error(
452                "impulse: denominator leading coefficient must be non-zero",
453            ));
454        }
455        let denominator: Vec<f64> = system
456            .denominator
457            .iter()
458            .map(|value| *value / leading)
459            .collect();
460        let mut numerator: Vec<f64> = system
461            .numerator
462            .iter()
463            .map(|value| *value / leading)
464            .collect();
465        let order = denominator.len() - 1;
466        while numerator.len() < order {
467            numerator.insert(0, 0.0);
468        }
469
470        let mut a = DMatrix::<f64>::zeros(order, order);
471        for row in 0..order.saturating_sub(1) {
472            a[(row, row + 1)] = 1.0;
473        }
474        for col in 0..order {
475            a[(order - 1, col)] = -denominator[order - col];
476        }
477        let c = numerator.into_iter().rev().collect();
478        Ok(Self { a, c })
479    }
480}
481
482fn continuous_response(realization: &Realization, t: &[f64]) -> Vec<f64> {
483    t.iter()
484        .map(|&time| {
485            let exp_at = matrix_exp(&(realization.a.clone() * time));
486            dot_c_with_last_column(&realization.c, &exp_at)
487        })
488        .collect()
489}
490
491fn discrete_response(
492    system: &TfSystem,
493    realization: &Realization,
494    t: &[f64],
495) -> BuiltinResult<Vec<f64>> {
496    if t.len() > MAX_DISCRETE_SAMPLES {
497        return Err(impulse_error(format!(
498            "impulse: discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"
499        )));
500    }
501    let sample_indices: Vec<usize> = t
502        .iter()
503        .map(|value| checked_discrete_sample_index(system, *value))
504        .collect::<BuiltinResult<_>>()?;
505    let max_index = sample_indices.iter().copied().max().unwrap_or(0);
506    let order = realization.c.len();
507    let value_count = max_index
508        .checked_add(1)
509        .ok_or_else(|| impulse_error("impulse: discrete sample index exceeds platform limits"))?;
510    let mut values = vec![0.0; value_count];
511    if order == 0 {
512        return Ok(sample_indices.into_iter().map(|idx| values[idx]).collect());
513    }
514
515    let mut state = vec![0.0; order];
516    state[order - 1] = 1.0;
517    let impulse_scale = 1.0 / system.sample_time;
518    for k in 1..=max_index {
519        values[k] = dot(&realization.c, &state) * impulse_scale;
520        state = mat_vec_mul(&realization.a, &state);
521    }
522    Ok(sample_indices.into_iter().map(|idx| values[idx]).collect())
523}
524
525fn dot_c_with_last_column(c: &[f64], matrix: &DMatrix<f64>) -> f64 {
526    if c.is_empty() {
527        return 0.0;
528    }
529    let last_col = matrix.ncols() - 1;
530    c.iter()
531        .enumerate()
532        .map(|(row, coeff)| coeff * matrix[(row, last_col)])
533        .sum()
534}
535
536fn dot(lhs: &[f64], rhs: &[f64]) -> f64 {
537    lhs.iter().zip(rhs).map(|(a, b)| a * b).sum()
538}
539
540fn mat_vec_mul(matrix: &DMatrix<f64>, vector: &[f64]) -> Vec<f64> {
541    let mut out = vec![0.0; matrix.nrows()];
542    for row in 0..matrix.nrows() {
543        let mut acc = 0.0;
544        for col in 0..matrix.ncols() {
545            acc += matrix[(row, col)] * vector[col];
546        }
547        out[row] = acc;
548    }
549    out
550}
551
552fn matrix_exp(matrix: &DMatrix<f64>) -> DMatrix<f64> {
553    let norm = matrix_one_norm(matrix);
554    let scale_power = if norm <= 0.5 {
555        0usize
556    } else {
557        norm.log2().ceil().max(0.0) as usize + 1
558    };
559    let scale = 2.0_f64.powi(scale_power as i32);
560    let scaled = matrix / scale;
561    let n = matrix.nrows();
562    let mut result = DMatrix::<f64>::identity(n, n);
563    let mut term = DMatrix::<f64>::identity(n, n);
564    for k in 1..=48 {
565        term = (&term * &scaled) / (k as f64);
566        result += &term;
567        if matrix_one_norm(&term) <= 1.0e-14 {
568            break;
569        }
570    }
571    for _ in 0..scale_power {
572        result = &result * &result;
573    }
574    result
575}
576
577fn matrix_one_norm(matrix: &DMatrix<f64>) -> f64 {
578    let mut best = 0.0;
579    for col in 0..matrix.ncols() {
580        let mut sum = 0.0;
581        for row in 0..matrix.nrows() {
582            sum += matrix[(row, col)].abs();
583        }
584        if sum > best {
585            best = sum;
586        }
587    }
588    best
589}
590
591#[cfg(feature = "plot-core")]
592async fn render_impulse_plot(response: &ImpulseResponse) -> BuiltinResult<Value> {
593    let t = response.t_value()?;
594    let y = response.y_value()?;
595    let plot_name = if response.discrete { "stem" } else { "plot" };
596    crate::dispatcher::call_builtin_async(plot_name, &[t, y]).await
597}
598
599#[cfg(not(feature = "plot-core"))]
600async fn render_impulse_plot(_response: &ImpulseResponse) -> BuiltinResult<Value> {
601    Ok(Value::Num(f64::NAN))
602}
603
604#[cfg(test)]
605mod tests {
606    use super::*;
607    use futures::executor::block_on;
608    use runmat_builtins::{CharArray, ObjectInstance};
609
610    fn tf_object(num: Vec<f64>, den: Vec<f64>, ts: f64) -> Value {
611        tf_object_with_delays(num, den, ts, 0.0, 0.0)
612    }
613
614    fn tf_object_with_delays(
615        num: Vec<f64>,
616        den: Vec<f64>,
617        ts: f64,
618        input_delay: f64,
619        output_delay: f64,
620    ) -> Value {
621        let mut object = ObjectInstance::new("tf".to_string());
622        object.properties.insert(
623            "Numerator".to_string(),
624            Value::Tensor(Tensor::new(num.clone(), vec![1, num.len()]).unwrap()),
625        );
626        object.properties.insert(
627            "Denominator".to_string(),
628            Value::Tensor(Tensor::new(den.clone(), vec![1, den.len()]).unwrap()),
629        );
630        object.properties.insert(
631            "Variable".to_string(),
632            Value::CharArray(CharArray::new_row(if ts > 0.0 { "z" } else { "s" })),
633        );
634        object.properties.insert("Ts".to_string(), Value::Num(ts));
635        object
636            .properties
637            .insert("InputDelay".to_string(), Value::Num(input_delay));
638        object
639            .properties
640            .insert("OutputDelay".to_string(), Value::Num(output_delay));
641        Value::Object(object)
642    }
643
644    fn run_impulse(system: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
645        block_on(impulse_builtin(system, rest))
646    }
647
648    fn tensor_data(value: Value) -> Vec<f64> {
649        match value {
650            Value::Tensor(tensor) => tensor.data,
651            other => panic!("expected tensor, got {other:?}"),
652        }
653    }
654
655    #[test]
656    fn impulse_first_order_continuous_explicit_time() {
657        let sys = tf_object(vec![20.0], vec![1.0, 5.0], 0.0);
658        let t = Value::Tensor(Tensor::new(vec![0.0, 0.1, 0.2], vec![1, 3]).unwrap());
659        let y = tensor_data(run_impulse(sys, vec![t]).expect("impulse"));
660        let expected = [20.0, 20.0 * (-0.5f64).exp(), 20.0 * (-1.0f64).exp()];
661        for (actual, expected) in y.iter().zip(expected) {
662            assert!((actual - expected).abs() < 1.0e-8);
663        }
664    }
665
666    #[test]
667    fn impulse_second_order_continuous() {
668        let sys = tf_object(vec![1.0], vec![1.0, 3.0, 2.0], 0.0);
669        let t = Value::Tensor(Tensor::new(vec![0.0, 0.5, 1.0], vec![1, 3]).unwrap());
670        let y = tensor_data(run_impulse(sys, vec![t]).expect("impulse"));
671        for (actual, time) in y.iter().zip([0.0_f64, 0.5, 1.0]) {
672            let expected = (-time).exp() - (-2.0 * time).exp();
673            assert!((actual - expected).abs() < 1.0e-8);
674        }
675    }
676
677    #[test]
678    fn impulse_multi_output_returns_y_and_t_columns() {
679        let _guard = crate::output_count::push_output_count(Some(2));
680        let sys = tf_object(vec![20.0], vec![1.0, 5.0], 0.0);
681        let t_arg = Value::Tensor(Tensor::new(vec![0.0, 0.1], vec![1, 2]).unwrap());
682        let result = run_impulse(sys, vec![t_arg]).expect("impulse");
683        let Value::OutputList(outputs) = result else {
684            panic!("expected output list");
685        };
686        assert_eq!(outputs.len(), 2);
687        match &outputs[0] {
688            Value::Tensor(tensor) => assert_eq!(tensor.shape, vec![2, 1]),
689            other => panic!("expected y tensor, got {other:?}"),
690        }
691        match &outputs[1] {
692            Value::Tensor(tensor) => {
693                assert_eq!(tensor.shape, vec![2, 1]);
694                assert_eq!(tensor.data, vec![0.0, 0.1]);
695            }
696            other => panic!("expected t tensor, got {other:?}"),
697        }
698    }
699
700    #[test]
701    fn impulse_discrete_siso_response() {
702        let sys = tf_object(vec![1.0], vec![1.0, -0.5], 0.1);
703        let t = Value::Tensor(Tensor::new(vec![0.0, 0.1, 0.2, 0.3], vec![1, 4]).unwrap());
704        let y = tensor_data(run_impulse(sys, vec![t]).expect("impulse"));
705        assert_eq!(y.len(), 4);
706        assert!((y[0] - 0.0).abs() < 1.0e-12);
707        assert!((y[1] - 10.0).abs() < 1.0e-12);
708        assert!((y[2] - 5.0).abs() < 1.0e-12);
709        assert!((y[3] - 2.5).abs() < 1.0e-12);
710    }
711
712    #[test]
713    fn impulse_discrete_final_time_rejects_excessive_sample_count() {
714        let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0e-6);
715        let err = run_impulse(sys, vec![Value::Num(2.0)]).expect_err("should fail");
716        assert!(err.message().contains("more than 1000000 samples"));
717    }
718
719    #[test]
720    fn impulse_discrete_time_vector_rejects_excessive_sample_index() {
721        let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0);
722        let t =
723            Value::Tensor(Tensor::new(vec![0.0, MAX_DISCRETE_SAMPLES as f64], vec![1, 2]).unwrap());
724        let err = run_impulse(sys, vec![t]).expect_err("should fail");
725        assert!(err.message().contains("more than 1000000 samples"));
726    }
727
728    #[test]
729    fn impulse_rejects_unsupported_model_type() {
730        let object = ObjectInstance::new("ss".to_string());
731        let err = run_impulse(Value::Object(object), Vec::new()).expect_err("should fail");
732        assert!(err.message().contains("unsupported model class"));
733    }
734
735    #[test]
736    fn impulse_rejects_direct_feedthrough_tf() {
737        let sys = tf_object(vec![1.0, 1.0], vec![1.0, 2.0], 0.0);
738        let err = run_impulse(sys, Vec::new()).expect_err("should fail");
739        assert!(err.message().contains("strictly proper"));
740    }
741
742    #[test]
743    fn impulse_rejects_invalid_time_metadata() {
744        let err = run_impulse(tf_object(vec![1.0], vec![1.0, -0.5], -0.1), Vec::new())
745            .expect_err("negative sample time should fail");
746        assert!(err.message().contains("Ts must be"));
747
748        let err = run_impulse(
749            tf_object_with_delays(vec![1.0], vec![1.0, 5.0], 0.0, f64::NAN, 0.0),
750            Vec::new(),
751        )
752        .expect_err("NaN input delay should fail");
753        assert!(err.message().contains("InputDelay must be"));
754
755        let err = run_impulse(
756            tf_object_with_delays(vec![1.0], vec![1.0, 5.0], 0.0, 0.0, -1.0),
757            Vec::new(),
758        )
759        .expect_err("negative output delay should fail");
760        assert!(err.message().contains("OutputDelay must be"));
761    }
762}
runmat_runtime/builtins/control/impulse.rs

runmat_runtime/builtins/control/
impulse.rs
