runmat_runtime/builtins/control/
step.rs

1//! MATLAB-compatible `step` response builtin for RunMat.
2
3use nalgebra::DMatrix;
4use num_complex::Complex64;
5use runmat_builtins::{ComplexTensor, ObjectInstance, Tensor, Value};
6use runmat_macros::runtime_builtin;
7
8use crate::builtins::common::spec::{
9    BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
10    ReductionNaN, ResidencyPolicy, ShapeRequirements,
11};
12use crate::builtins::control::type_resolvers::step_type;
13use crate::{build_runtime_error, BuiltinResult, RuntimeError};
14
15const BUILTIN_NAME: &str = "step";
16const EPS: f64 = 1.0e-12;
17const DEFAULT_SAMPLES: usize = 101;
18const MAX_DISCRETE_SAMPLES: usize = 1_000_000;
19
20#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::control::step")]
21pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
22    name: "step",
23    op_kind: GpuOpKind::Custom("control-step-response"),
24    supported_precisions: &[],
25    broadcast: BroadcastSemantics::None,
26    provider_hooks: &[],
27    constant_strategy: ConstantStrategy::InlineLiteral,
28    residency: ResidencyPolicy::GatherImmediately,
29    nan_mode: ReductionNaN::Include,
30    two_pass_threshold: None,
31    workgroup_size: None,
32    accepts_nan_mode: false,
33    notes: "Step-response simulation runs on the host from transfer-function metadata.",
34};
35
36#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::control::step")]
37pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
38    name: "step",
39    shape: ShapeRequirements::Any,
40    constant_strategy: ConstantStrategy::InlineLiteral,
41    elementwise: None,
42    reduction: None,
43    emits_nan: false,
44    notes: "step simulates a dynamic system and terminates numeric fusion chains.",
45};
46
47fn step_error(message: impl Into<String>) -> RuntimeError {
48    build_runtime_error(message)
49        .with_builtin(BUILTIN_NAME)
50        .build()
51}
52
53#[runtime_builtin(
54    name = "step",
55    category = "control",
56    summary = "Compute or plot the step response of a SISO transfer-function model.",
57    keywords = "step,response,control system,transfer function,tf",
58    sink = true,
59    suppress_auto_output = true,
60    type_resolver(step_type),
61    builtin_path = "crate::builtins::control::step"
62)]
63async fn step_builtin(sys: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
64    if rest.len() > 1 {
65        return Err(step_error(
66            "step: expected step(sys), step(sys, tFinal), or step(sys, t)",
67        ));
68    }
69
70    let sys = crate::gather_if_needed_async(&sys).await?;
71    let mut rest_host = Vec::with_capacity(rest.len());
72    for arg in &rest {
73        rest_host.push(crate::gather_if_needed_async(arg).await?);
74    }
75    let model = TransferFunction::from_value(sys)?;
76    let time = TimeSpec::parse(rest_host.first(), model.sample_time)?;
77    let eval = evaluate_step(&model, time)?;
78
79    if crate::output_context::requested_output_count() == Some(0)
80        && crate::output_count::current_output_count().is_none()
81    {
82        plot_response(&eval).await?;
83        return Ok(Value::OutputList(Vec::new()));
84    }
85
86    if let Some(out_count) = crate::output_count::current_output_count() {
87        if out_count == 0 {
88            plot_response(&eval).await?;
89            return Ok(Value::OutputList(Vec::new()));
90        }
91        if out_count == 1 {
92            return Ok(Value::OutputList(vec![eval.y_value()?]));
93        }
94        return Ok(crate::output_count::output_list_with_padding(
95            out_count,
96            eval.outputs()?,
97        ));
98    }
99
100    eval.y_value()
101}
102
103#[derive(Clone, Debug)]
104struct TransferFunction {
105    numerator: Vec<f64>,
106    denominator: Vec<f64>,
107    sample_time: f64,
108}
109
110impl TransferFunction {
111    fn from_value(value: Value) -> BuiltinResult<Self> {
112        let Value::Object(object) = value else {
113            return Err(step_error("step: expected a tf object"));
114        };
115        if !object.is_class("tf") {
116            return Err(step_error(format!(
117                "step: expected a tf object, got {}",
118                object.class_name
119            )));
120        }
121
122        let numerator = property_coefficients(&object, "Numerator")?;
123        let denominator = property_coefficients(&object, "Denominator")?;
124        let sample_time = property_scalar(&object, "Ts")?;
125        if !sample_time.is_finite() || sample_time < 0.0 {
126            return Err(step_error(
127                "step: tf sample time must be finite and non-negative",
128            ));
129        }
130
131        let numerator = trim_leading_zeros(numerator);
132        let denominator = trim_leading_zeros(denominator);
133        if denominator.is_empty() {
134            return Err(step_error("step: denominator coefficients cannot be empty"));
135        }
136        if denominator[0].abs() <= EPS {
137            return Err(step_error(
138                "step: leading denominator coefficient must be non-zero",
139            ));
140        }
141        if numerator.len().saturating_sub(1) > denominator.len().saturating_sub(1) {
142            return Err(step_error(
143                "step: improper transfer functions are not supported yet",
144            ));
145        }
146
147        Ok(Self {
148            numerator,
149            denominator,
150            sample_time,
151        })
152    }
153
154    fn normalized(&self) -> (Vec<f64>, Vec<f64>) {
155        let leading = self.denominator[0];
156        let den = self
157            .denominator
158            .iter()
159            .map(|value| value / leading)
160            .collect::<Vec<_>>();
161        let num = self
162            .numerator
163            .iter()
164            .map(|value| value / leading)
165            .collect::<Vec<_>>();
166        (num, den)
167    }
168}
169
170fn property_coefficients(object: &ObjectInstance, name: &str) -> BuiltinResult<Vec<f64>> {
171    let value = object
172        .properties
173        .get(name)
174        .ok_or_else(|| step_error(format!("step: tf object is missing {name}")))?;
175    match value {
176        Value::Tensor(tensor) => Ok(tensor.data.clone()),
177        Value::ComplexTensor(tensor) => real_complex_coefficients(tensor, name),
178        Value::Num(n) => Ok(vec![*n]),
179        Value::Int(i) => Ok(vec![i.to_f64()]),
180        other => Err(step_error(format!(
181            "step: tf {name} coefficients must be numeric, got {other:?}"
182        ))),
183    }
184}
185
186fn real_complex_coefficients(tensor: &ComplexTensor, name: &str) -> BuiltinResult<Vec<f64>> {
187    let mut out = Vec::with_capacity(tensor.data.len());
188    for &(re, im) in &tensor.data {
189        if im.abs() > EPS {
190            return Err(step_error(format!(
191                "step: complex tf {name} coefficients are not supported yet"
192            )));
193        }
194        out.push(re);
195    }
196    Ok(out)
197}
198
199fn property_scalar(object: &ObjectInstance, name: &str) -> BuiltinResult<f64> {
200    let value = object
201        .properties
202        .get(name)
203        .ok_or_else(|| step_error(format!("step: tf object is missing {name}")))?;
204    match value {
205        Value::Num(n) => Ok(*n),
206        Value::Int(i) => Ok(i.to_f64()),
207        other => Err(step_error(format!(
208            "step: tf {name} property must be a scalar, got {other:?}"
209        ))),
210    }
211}
212
213#[derive(Clone, Debug)]
214enum TimeSpec {
215    Auto,
216    FinalTime(f64),
217    Vector(Vec<f64>),
218}
219
220impl TimeSpec {
221    fn parse(value: Option<&Value>, sample_time: f64) -> BuiltinResult<Self> {
222        let Some(value) = value else {
223            return Ok(Self::Auto);
224        };
225        match value {
226            Value::Num(n) => Self::final_time(*n),
227            Value::Int(i) => Self::final_time(i.to_f64()),
228            Value::Tensor(tensor) => {
229                ensure_time_vector_shape(&tensor.shape)?;
230                if tensor.data.len() == 1 {
231                    return Self::final_time(tensor.data[0]);
232                }
233                Self::vector(tensor.data.clone(), sample_time)
234            }
235            other => Err(step_error(format!(
236                "step: time input must be a scalar final time or numeric vector, got {other:?}"
237            ))),
238        }
239    }
240
241    fn final_time(value: f64) -> BuiltinResult<Self> {
242        if !value.is_finite() || value <= 0.0 {
243            return Err(step_error(
244                "step: final time must be a positive finite scalar",
245            ));
246        }
247        Ok(Self::FinalTime(value))
248    }
249
250    fn vector(values: Vec<f64>, sample_time: f64) -> BuiltinResult<Self> {
251        validate_time_vector(&values)?;
252        if sample_time > 0.0 {
253            for &t in &values {
254                let k = (t / sample_time).round();
255                if (t - k * sample_time).abs() > 1.0e-8 * sample_time.max(1.0) {
256                    return Err(step_error(
257                        "step: discrete-time sample vector must align with the model sample time",
258                    ));
259                }
260            }
261        }
262        Ok(Self::Vector(values))
263    }
264}
265
266fn ensure_time_vector_shape(shape: &[usize]) -> BuiltinResult<()> {
267    let non_unit = shape.iter().copied().filter(|&dim| dim > 1).count();
268    if non_unit <= 1 {
269        Ok(())
270    } else {
271        Err(step_error("step: time input must be a vector"))
272    }
273}
274
275fn validate_time_vector(values: &[f64]) -> BuiltinResult<()> {
276    if values.is_empty() {
277        return Err(step_error("step: time vector must not be empty"));
278    }
279    let mut previous = None;
280    for &value in values {
281        if !value.is_finite() || value < 0.0 {
282            return Err(step_error(
283                "step: time vector values must be finite and non-negative",
284            ));
285        }
286        if let Some(prev) = previous {
287            if value < prev {
288                return Err(step_error("step: time vector must be nondecreasing"));
289            }
290        }
291        previous = Some(value);
292    }
293    Ok(())
294}
295
296#[derive(Clone, Debug)]
297struct StepEval {
298    y: Vec<f64>,
299    t: Vec<f64>,
300}
301
302impl StepEval {
303    fn y_value(&self) -> BuiltinResult<Value> {
304        column_tensor(self.y.clone())
305    }
306
307    fn t_value(&self) -> BuiltinResult<Value> {
308        column_tensor(self.t.clone())
309    }
310
311    fn outputs(&self) -> BuiltinResult<Vec<Value>> {
312        Ok(vec![self.y_value()?, self.t_value()?])
313    }
314}
315
316fn evaluate_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
317    if model.sample_time > 0.0 {
318        evaluate_discrete_step(model, time)
319    } else {
320        evaluate_continuous_step(model, time)
321    }
322}
323
324fn evaluate_continuous_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
325    let t = continuous_time_vector(model, time)?;
326    let (num, den) = model.normalized();
327    let response = continuous_response(&num, &den, &t)?;
328    Ok(StepEval { y: response, t })
329}
330
331fn continuous_time_vector(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<Vec<f64>> {
332    match time {
333        TimeSpec::Auto => Ok(linspace(0.0, automatic_final_time(model), DEFAULT_SAMPLES)),
334        TimeSpec::FinalTime(final_time) => Ok(linspace(0.0, final_time, DEFAULT_SAMPLES)),
335        TimeSpec::Vector(values) => Ok(values),
336    }
337}
338
339fn continuous_response(num: &[f64], den: &[f64], t: &[f64]) -> BuiltinResult<Vec<f64>> {
340    let order = den.len() - 1;
341    if order == 0 {
342        let gain = num.last().copied().unwrap_or(0.0) / den[0];
343        return Ok(vec![gain; t.len()]);
344    }
345
346    let mut padded_num = vec![0.0; order + 1 - num.len()];
347    padded_num.extend_from_slice(num);
348    let direct = padded_num[0];
349    let a = &den[1..];
350    let mut c = Vec::with_capacity(order);
351    for state_idx in 0..order {
352        let coeff_idx = order - state_idx;
353        c.push(padded_num[coeff_idx] - direct * den[coeff_idx]);
354    }
355
356    let mut state = vec![0.0; order];
357    let mut current_t = 0.0;
358    let mut response = Vec::with_capacity(t.len());
359    for &target_t in t {
360        if target_t > current_t {
361            integrate_to(&mut state, a, current_t, target_t);
362            current_t = target_t;
363        }
364        response.push(dot(&c, &state) + direct);
365    }
366    Ok(response)
367}
368
369fn integrate_to(state: &mut [f64], a: &[f64], start: f64, end: f64) {
370    let duration = end - start;
371    if duration <= 0.0 {
372        return;
373    }
374    let steps = ((duration / 0.01).ceil() as usize).clamp(1, 10_000);
375    let h = duration / steps as f64;
376    for _ in 0..steps {
377        rk4_step(state, a, h);
378    }
379}
380
381fn rk4_step(state: &mut [f64], a: &[f64], h: f64) {
382    let k1 = derivative(state, a);
383    let s2 = add_scaled(state, &k1, h * 0.5);
384    let k2 = derivative(&s2, a);
385    let s3 = add_scaled(state, &k2, h * 0.5);
386    let k3 = derivative(&s3, a);
387    let s4 = add_scaled(state, &k3, h);
388    let k4 = derivative(&s4, a);
389    for idx in 0..state.len() {
390        state[idx] += h * (k1[idx] + 2.0 * k2[idx] + 2.0 * k3[idx] + k4[idx]) / 6.0;
391    }
392}
393
394fn derivative(state: &[f64], a: &[f64]) -> Vec<f64> {
395    let order = state.len();
396    let mut dx = vec![0.0; order];
397    if order > 1 {
398        dx[..(order - 1)].copy_from_slice(&state[1..order]);
399    }
400    let mut last = 1.0;
401    for state_idx in 0..order {
402        let coeff = a[order - 1 - state_idx];
403        last -= coeff * state[state_idx];
404    }
405    dx[order - 1] = last;
406    dx
407}
408
409fn add_scaled(state: &[f64], delta: &[f64], scale: f64) -> Vec<f64> {
410    state
411        .iter()
412        .zip(delta)
413        .map(|(value, delta)| value + scale * delta)
414        .collect()
415}
416
417fn evaluate_discrete_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
418    let t = discrete_time_vector(model.sample_time, time)?;
419    if t.len() > MAX_DISCRETE_SAMPLES {
420        return Err(step_error(format!(
421            "step: discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"
422        )));
423    }
424    let sample_indices = t
425        .iter()
426        .map(|&value| checked_discrete_sample_index(model.sample_time, value))
427        .collect::<BuiltinResult<Vec<_>>>()?;
428    let max_k = sample_indices.iter().copied().max().unwrap_or(0);
429    let count = max_k
430        .checked_add(1)
431        .ok_or_else(|| step_error("step: discrete sample index exceeds platform limits"))?;
432    let (num, den) = model.normalized();
433    let all_y = discrete_response(&num, &den, count)?;
434    let y = sample_indices
435        .into_iter()
436        .map(|idx| {
437            all_y
438                .get(idx)
439                .copied()
440                .ok_or_else(|| step_error("step: discrete sample index exceeds response length"))
441        })
442        .collect::<BuiltinResult<Vec<_>>>()?;
443    Ok(StepEval { y, t })
444}
445
446fn discrete_time_vector(sample_time: f64, time: TimeSpec) -> BuiltinResult<Vec<f64>> {
447    match time {
448        TimeSpec::Auto => Ok((0..DEFAULT_SAMPLES)
449            .map(|idx| idx as f64 * sample_time)
450            .collect()),
451        TimeSpec::FinalTime(final_time) => {
452            let steps = checked_discrete_sample_steps(sample_time, final_time)?;
453            Ok((0..=steps).map(|idx| idx as f64 * sample_time).collect())
454        }
455        TimeSpec::Vector(values) => Ok(values),
456    }
457}
458
459fn checked_discrete_sample_steps(sample_time: f64, final_time: f64) -> BuiltinResult<usize> {
460    let steps = (final_time / sample_time).floor();
461    if !steps.is_finite() || steps < 0.0 || steps > usize::MAX as f64 {
462        return Err(step_error(
463            "step: discrete sample count exceeds platform limits",
464        ));
465    }
466    if steps >= MAX_DISCRETE_SAMPLES as f64 {
467        return Err(step_error(format!(
468            "step: discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"
469        )));
470    }
471    Ok(steps as usize)
472}
473
474fn checked_discrete_sample_index(sample_time: f64, time: f64) -> BuiltinResult<usize> {
475    let index = (time / sample_time).round();
476    if !index.is_finite() || index < 0.0 || index > usize::MAX as f64 {
477        return Err(step_error(
478            "step: discrete sample index exceeds platform limits",
479        ));
480    }
481    if index >= MAX_DISCRETE_SAMPLES as f64 {
482        return Err(step_error(format!(
483            "step: discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"
484        )));
485    }
486    Ok(index as usize)
487}
488
489fn discrete_response(num: &[f64], den: &[f64], count: usize) -> BuiltinResult<Vec<f64>> {
490    let order = den.len() - 1;
491    let mut padded_num = vec![0.0; order + 1 - num.len()];
492    padded_num.extend_from_slice(num);
493    let mut y = vec![0.0; count];
494    for k in 0..count {
495        let mut value = 0.0;
496        for (idx, &coeff) in padded_num.iter().enumerate() {
497            if k >= idx {
498                value += coeff;
499            }
500        }
501        for idx in 1..den.len() {
502            if k >= idx {
503                value -= den[idx] * y[k - idx];
504            }
505        }
506        y[k] = value;
507    }
508    Ok(y)
509}
510
511async fn plot_response(eval: &StepEval) -> BuiltinResult<()> {
512    let t = eval.t_value()?;
513    let y = eval.y_value()?;
514    if let Err(err) = crate::call_builtin_async("plot", &[t, y]).await {
515        if is_nonfatal_plot_setup_error(&err) {
516            return Ok(());
517        }
518        return Err(err);
519    }
520    let _ = crate::call_builtin_async("title", &[Value::from("Step Response")]).await;
521    let _ = crate::call_builtin_async("xlabel", &[Value::from("Time")]).await;
522    let _ = crate::call_builtin_async("ylabel", &[Value::from("Amplitude")]).await;
523    Ok(())
524}
525
526fn is_nonfatal_plot_setup_error(err: &RuntimeError) -> bool {
527    let lower = err.to_string().to_ascii_lowercase();
528    lower.contains("plotting is unavailable")
529        || lower.contains("non-main thread")
530        || lower.contains("interactive plotting failed")
531        || lower.contains("eventloop can't be recreated")
532}
533
534fn automatic_final_time(model: &TransferFunction) -> f64 {
535    let (_, den) = model.normalized();
536    let poles = polynomial_roots(&den).unwrap_or_default();
537    let slowest_decay = poles
538        .iter()
539        .filter_map(|pole| if pole.re < -EPS { Some(-pole.re) } else { None })
540        .fold(f64::INFINITY, f64::min);
541    if slowest_decay.is_finite() && slowest_decay > EPS {
542        (5.0 / slowest_decay).clamp(1.0, 100.0)
543    } else {
544        10.0
545    }
546}
547
548fn polynomial_roots(coeffs: &[f64]) -> BuiltinResult<Vec<Complex64>> {
549    let trimmed = trim_leading_zeros(coeffs.to_vec());
550    if trimmed.len() <= 1 {
551        return Ok(Vec::new());
552    }
553    if trimmed.len() == 2 {
554        return Ok(vec![Complex64::new(-trimmed[1] / trimmed[0], 0.0)]);
555    }
556    let degree = trimmed.len() - 1;
557    let leading = trimmed[0];
558    let mut companion = DMatrix::<Complex64>::zeros(degree, degree);
559    for row in 1..degree {
560        companion[(row, row - 1)] = Complex64::new(1.0, 0.0);
561    }
562    for (idx, coeff) in trimmed.iter().enumerate().skip(1) {
563        companion[(0, idx - 1)] = Complex64::new(-coeff / leading, 0.0);
564    }
565    let eigenvalues = companion
566        .eigenvalues()
567        .ok_or_else(|| step_error("step: failed to compute transfer-function poles"))?;
568    Ok(eigenvalues.iter().copied().collect())
569}
570
571fn trim_leading_zeros(coeffs: Vec<f64>) -> Vec<f64> {
572    let first_nonzero = coeffs
573        .iter()
574        .position(|value| value.abs() > EPS)
575        .unwrap_or(coeffs.len());
576    coeffs[first_nonzero..].to_vec()
577}
578
579fn linspace(start: f64, end: f64, count: usize) -> Vec<f64> {
580    if count <= 1 {
581        return vec![end];
582    }
583    let step = (end - start) / (count - 1) as f64;
584    (0..count).map(|idx| start + idx as f64 * step).collect()
585}
586
587fn dot(lhs: &[f64], rhs: &[f64]) -> f64 {
588    lhs.iter().zip(rhs).map(|(a, b)| a * b).sum()
589}
590
591fn column_tensor(data: Vec<f64>) -> BuiltinResult<Value> {
592    let rows = data.len();
593    let tensor =
594        Tensor::new(data, vec![rows, 1]).map_err(|err| step_error(format!("step: {err}")))?;
595    Ok(Value::Tensor(tensor))
596}
597
598#[cfg(test)]
599mod tests {
600    use super::*;
601    use futures::executor::block_on;
602    use runmat_builtins::{CharArray, ObjectInstance};
603
604    fn tf_object(num: Vec<f64>, den: Vec<f64>, sample_time: f64) -> Value {
605        let mut object = ObjectInstance::new("tf".to_string());
606        object.properties.insert(
607            "Numerator".to_string(),
608            Value::Tensor(Tensor::new(num.clone(), vec![1, num.len()]).unwrap()),
609        );
610        object.properties.insert(
611            "Denominator".to_string(),
612            Value::Tensor(Tensor::new(den.clone(), vec![1, den.len()]).unwrap()),
613        );
614        object.properties.insert(
615            "Variable".to_string(),
616            Value::CharArray(CharArray::new_row(if sample_time > 0.0 {
617                "z"
618            } else {
619                "s"
620            })),
621        );
622        object
623            .properties
624            .insert("Ts".to_string(), Value::Num(sample_time));
625        Value::Object(object)
626    }
627
628    fn run_step(sys: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
629        block_on(step_builtin(sys, rest))
630    }
631
632    fn tensor_data(value: Value) -> Vec<f64> {
633        match value {
634            Value::Tensor(tensor) => tensor.data,
635            other => panic!("expected tensor, got {other:?}"),
636        }
637    }
638
639    #[test]
640    fn first_order_continuous_response_matches_closed_form_for_explicit_time() {
641        let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
642        let time = Value::Tensor(Tensor::new(vec![0.0, 0.5, 1.0, 2.0], vec![1, 4]).unwrap());
643        let y = tensor_data(run_step(sys, vec![time]).expect("step"));
644        for (actual, t) in y.iter().zip([0.0_f64, 0.5, 1.0, 2.0]) {
645            let expected = 1.0 - (-t).exp();
646            assert!(
647                (actual - expected).abs() < 1.0e-5,
648                "t={t} actual={actual} expected={expected}"
649            );
650        }
651    }
652
653    #[test]
654    fn multi_output_returns_y_then_time() {
655        let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
656        let _guard = crate::output_count::push_output_count(Some(2));
657        let time = Value::Tensor(Tensor::new(vec![0.0, 1.0], vec![2, 1]).unwrap());
658        let result = run_step(sys, vec![time]).expect("step");
659        let Value::OutputList(outputs) = result else {
660            panic!("expected output list");
661        };
662        assert_eq!(outputs.len(), 2);
663        assert_eq!(tensor_data(outputs[1].clone()), vec![0.0, 1.0]);
664        let y = tensor_data(outputs[0].clone());
665        assert!((y[1] - (1.0 - (-1.0_f64).exp())).abs() < 1.0e-5);
666    }
667
668    #[test]
669    fn single_requested_output_returns_only_response() {
670        let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
671        let _guard = crate::output_count::push_output_count(Some(1));
672        let time = Value::Tensor(Tensor::new(vec![0.0, 1.0], vec![1, 2]).unwrap());
673        let result = run_step(sys, vec![time]).expect("step");
674        let Value::OutputList(outputs) = result else {
675            panic!("expected output list");
676        };
677        assert_eq!(outputs.len(), 1);
678        let y = tensor_data(outputs[0].clone());
679        assert_eq!(y.len(), 2);
680        assert!((y[1] - (1.0 - (-1.0_f64).exp())).abs() < 1.0e-5);
681    }
682
683    #[test]
684    fn scalar_final_time_generates_column_time_vector_ending_at_final_time() {
685        let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
686        let _guard = crate::output_count::push_output_count(Some(2));
687        let result = run_step(sys, vec![Value::Num(5.0)]).expect("step");
688        let Value::OutputList(outputs) = result else {
689            panic!("expected output list");
690        };
691        let t = tensor_data(outputs[1].clone());
692        assert_eq!(t.len(), DEFAULT_SAMPLES);
693        assert_eq!(t[0], 0.0);
694        assert!((t[t.len() - 1] - 5.0).abs() < 1.0e-12);
695    }
696
697    #[test]
698    fn discrete_response_uses_sample_time_grid() {
699        let sys = tf_object(vec![1.0], vec![1.0, -0.5], 0.1);
700        let time = Value::Tensor(Tensor::new(vec![0.0, 0.1, 0.2], vec![1, 3]).unwrap());
701        let y = tensor_data(run_step(sys, vec![time]).expect("step"));
702        assert_eq!(y, vec![0.0, 1.0, 1.5]);
703    }
704
705    #[test]
706    fn discrete_final_time_rejects_excessive_sample_count() {
707        let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0e-6);
708        let err = run_step(sys, vec![Value::Num(2.0)]).expect_err("should fail");
709        assert!(err.message().contains("more than 1000000 samples"));
710    }
711
712    #[test]
713    fn discrete_time_vector_rejects_excessive_sample_index() {
714        let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0);
715        let time =
716            Value::Tensor(Tensor::new(vec![0.0, MAX_DISCRETE_SAMPLES as f64], vec![1, 2]).unwrap());
717        let err = run_step(sys, vec![time]).expect_err("should fail");
718        assert!(err.message().contains("more than 1000000 samples"));
719    }
720
721    #[test]
722    fn rejects_non_tf_input() {
723        let err = run_step(Value::Num(1.0), Vec::new()).expect_err("expected error");
724        assert!(err.message().contains("expected a tf object"));
725    }
726}
runmat_runtime/builtins/control/step.rs

runmat_runtime/builtins/control/
step.rs
