1use nalgebra::DMatrix;
4use num_complex::Complex64;
5use runmat_builtins::{
6 BuiltinCompletionPolicy, BuiltinDescriptor, BuiltinErrorDescriptor, BuiltinOutputMode,
7 BuiltinParamArity, BuiltinParamDescriptor, BuiltinParamType, BuiltinSignatureDescriptor,
8 ComplexTensor, ObjectInstance, Tensor, Value,
9};
10use runmat_macros::runtime_builtin;
11
12use crate::builtins::common::spec::{
13 BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
14 ReductionNaN, ResidencyPolicy, ShapeRequirements,
15};
16use crate::builtins::control::type_resolvers::step_type;
17use crate::{build_runtime_error, BuiltinResult, RuntimeError};
18
19const BUILTIN_NAME: &str = "step";
20const EPS: f64 = 1.0e-12;
21const DEFAULT_SAMPLES: usize = 101;
22const MAX_DISCRETE_SAMPLES: usize = 1_000_000;
23
24const STEP_OUTPUT_Y: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
25 name: "y",
26 ty: BuiltinParamType::NumericArray,
27 arity: BuiltinParamArity::Required,
28 default: None,
29 description: "Step response samples (column vector).",
30}];
31const STEP_OUTPUT_Y_T: [BuiltinParamDescriptor; 2] = [
32 BuiltinParamDescriptor {
33 name: "y",
34 ty: BuiltinParamType::NumericArray,
35 arity: BuiltinParamArity::Required,
36 default: None,
37 description: "Step response samples (column vector).",
38 },
39 BuiltinParamDescriptor {
40 name: "t",
41 ty: BuiltinParamType::NumericArray,
42 arity: BuiltinParamArity::Required,
43 default: None,
44 description: "Time samples (column vector).",
45 },
46];
47const STEP_INPUTS_SYS: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
48 name: "sys",
49 ty: BuiltinParamType::Any,
50 arity: BuiltinParamArity::Required,
51 default: None,
52 description: "SISO tf model.",
53}];
54const STEP_INPUTS_SYS_TIME: [BuiltinParamDescriptor; 2] = [
55 BuiltinParamDescriptor {
56 name: "sys",
57 ty: BuiltinParamType::Any,
58 arity: BuiltinParamArity::Required,
59 default: None,
60 description: "SISO tf model.",
61 },
62 BuiltinParamDescriptor {
63 name: "time",
64 ty: BuiltinParamType::Any,
65 arity: BuiltinParamArity::Optional,
66 default: None,
67 description: "Final time scalar or explicit time vector.",
68 },
69];
70const STEP_SIGNATURES: [BuiltinSignatureDescriptor; 6] = [
71 BuiltinSignatureDescriptor {
72 label: "y = step(sys)",
73 inputs: &STEP_INPUTS_SYS,
74 outputs: &STEP_OUTPUT_Y,
75 },
76 BuiltinSignatureDescriptor {
77 label: "y = step(sys, tFinal)",
78 inputs: &STEP_INPUTS_SYS_TIME,
79 outputs: &STEP_OUTPUT_Y,
80 },
81 BuiltinSignatureDescriptor {
82 label: "y = step(sys, t)",
83 inputs: &STEP_INPUTS_SYS_TIME,
84 outputs: &STEP_OUTPUT_Y,
85 },
86 BuiltinSignatureDescriptor {
87 label: "[y,t] = step(sys)",
88 inputs: &STEP_INPUTS_SYS,
89 outputs: &STEP_OUTPUT_Y_T,
90 },
91 BuiltinSignatureDescriptor {
92 label: "[y,t] = step(sys, tFinal)",
93 inputs: &STEP_INPUTS_SYS_TIME,
94 outputs: &STEP_OUTPUT_Y_T,
95 },
96 BuiltinSignatureDescriptor {
97 label: "[y,t] = step(sys, t)",
98 inputs: &STEP_INPUTS_SYS_TIME,
99 outputs: &STEP_OUTPUT_Y_T,
100 },
101];
102const STEP_ERROR_INVALID_ARGUMENT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
103 code: "RM.STEP.INVALID_ARGUMENT",
104 identifier: Some("RunMat:step:InvalidArgument"),
105 when: "Inputs do not match supported step invocation forms.",
106 message: "step: invalid argument",
107};
108const STEP_ERROR_INVALID_MODEL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
109 code: "RM.STEP.INVALID_MODEL",
110 identifier: Some("RunMat:step:InvalidModel"),
111 when: "Input system is not a supported tf object with valid required properties.",
112 message: "step: invalid model",
113};
114const STEP_ERROR_INVALID_TIME: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
115 code: "RM.STEP.INVALID_TIME",
116 identifier: Some("RunMat:step:InvalidTime"),
117 when: "Time argument is invalid for the model class or sampling mode.",
118 message: "step: invalid time input",
119};
120const STEP_ERROR_UNSUPPORTED_MODEL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
121 code: "RM.STEP.UNSUPPORTED_MODEL",
122 identifier: Some("RunMat:step:UnsupportedModel"),
123 when: "Model is well-formed but unsupported by the current step implementation.",
124 message: "step: unsupported model",
125};
126const STEP_ERROR_DISCRETE_LIMIT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
127 code: "RM.STEP.DISCRETE_LIMIT",
128 identifier: Some("RunMat:step:DiscreteLimit"),
129 when: "Discrete simulation would exceed platform or configured sample limits.",
130 message: "step: discrete simulation limit exceeded",
131};
132const STEP_ERROR_PLOT_FAILED: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
133 code: "RM.STEP.PLOT_FAILED",
134 identifier: Some("RunMat:step:PlotFailed"),
135 when: "Statement-form plotting failed for reasons other than known nonfatal setup conditions.",
136 message: "step: plotting failed",
137};
138const STEP_ERROR_INTERNAL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
139 code: "RM.STEP.INTERNAL",
140 identifier: Some("RunMat:step:Internal"),
141 when: "Internal response assembly failed.",
142 message: "step: internal error",
143};
144const STEP_ERRORS: [BuiltinErrorDescriptor; 7] = [
145 STEP_ERROR_INVALID_ARGUMENT,
146 STEP_ERROR_INVALID_MODEL,
147 STEP_ERROR_INVALID_TIME,
148 STEP_ERROR_UNSUPPORTED_MODEL,
149 STEP_ERROR_DISCRETE_LIMIT,
150 STEP_ERROR_PLOT_FAILED,
151 STEP_ERROR_INTERNAL,
152];
153pub const STEP_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
154 signatures: &STEP_SIGNATURES,
155 output_mode: BuiltinOutputMode::ByRequestedOutputCount,
156 completion_policy: BuiltinCompletionPolicy::Public,
157 errors: &STEP_ERRORS,
158};
159
160#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::control::step")]
161pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
162 name: "step",
163 op_kind: GpuOpKind::Custom("control-step-response"),
164 supported_precisions: &[],
165 broadcast: BroadcastSemantics::None,
166 provider_hooks: &[],
167 constant_strategy: ConstantStrategy::InlineLiteral,
168 residency: ResidencyPolicy::GatherImmediately,
169 nan_mode: ReductionNaN::Include,
170 two_pass_threshold: None,
171 workgroup_size: None,
172 accepts_nan_mode: false,
173 notes: "Step-response simulation runs on the host from transfer-function metadata.",
174};
175
176#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::control::step")]
177pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
178 name: "step",
179 shape: ShapeRequirements::Any,
180 constant_strategy: ConstantStrategy::InlineLiteral,
181 elementwise: None,
182 reduction: None,
183 emits_nan: false,
184 notes: "step simulates a dynamic system and terminates numeric fusion chains.",
185};
186
187fn step_error_with_detail(
188 error: &'static BuiltinErrorDescriptor,
189 detail: impl AsRef<str>,
190) -> RuntimeError {
191 step_error_with_message(format!("{}: {}", error.message, detail.as_ref()), error)
192}
193
194fn step_error_with_message(
195 message: impl Into<String>,
196 error: &'static BuiltinErrorDescriptor,
197) -> RuntimeError {
198 let mut builder = build_runtime_error(message).with_builtin(BUILTIN_NAME);
199 if let Some(identifier) = error.identifier {
200 builder = builder.with_identifier(identifier);
201 }
202 builder.build()
203}
204
205#[runtime_builtin(
206 name = "step",
207 category = "control",
208 summary = "Compute or plot step responses of SISO transfer-function models.",
209 keywords = "step,response,control system,transfer function,tf",
210 sink = true,
211 suppress_auto_output = true,
212 type_resolver(step_type),
213 descriptor(crate::builtins::control::step::STEP_DESCRIPTOR),
214 builtin_path = "crate::builtins::control::step"
215)]
216async fn step_builtin(sys: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
217 if rest.len() > 1 {
218 return Err(step_error_with_detail(
219 &STEP_ERROR_INVALID_ARGUMENT,
220 "expected step(sys), step(sys, tFinal), or step(sys, t)",
221 ));
222 }
223
224 let sys = crate::gather_if_needed_async(&sys).await?;
225 let mut rest_host = Vec::with_capacity(rest.len());
226 for arg in &rest {
227 rest_host.push(crate::gather_if_needed_async(arg).await?);
228 }
229 let model = TransferFunction::from_value(sys)?;
230 let time = TimeSpec::parse(rest_host.first(), model.sample_time)?;
231 let eval = evaluate_step(&model, time)?;
232
233 if crate::output_context::requested_output_count() == Some(0)
234 && crate::output_count::current_output_count().is_none()
235 {
236 plot_response(&eval).await?;
237 return Ok(Value::OutputList(Vec::new()));
238 }
239
240 if let Some(out_count) = crate::output_count::current_output_count() {
241 if out_count == 0 {
242 plot_response(&eval).await?;
243 return Ok(Value::OutputList(Vec::new()));
244 }
245 if out_count == 1 {
246 return Ok(Value::OutputList(vec![eval.y_value()?]));
247 }
248 return Ok(crate::output_count::output_list_with_padding(
249 out_count,
250 eval.outputs()?,
251 ));
252 }
253
254 eval.y_value()
255}
256
257#[derive(Clone, Debug)]
258struct TransferFunction {
259 numerator: Vec<f64>,
260 denominator: Vec<f64>,
261 sample_time: f64,
262}
263
264impl TransferFunction {
265 fn from_value(value: Value) -> BuiltinResult<Self> {
266 let Value::Object(object) = value else {
267 return Err(step_error_with_detail(
268 &STEP_ERROR_INVALID_MODEL,
269 "expected a tf object",
270 ));
271 };
272 if !object.is_class("tf") {
273 return Err(step_error_with_detail(
274 &STEP_ERROR_INVALID_MODEL,
275 format!("expected a tf object, got {}", object.class_name),
276 ));
277 }
278
279 let numerator = property_coefficients(&object, "Numerator")?;
280 let denominator = property_coefficients(&object, "Denominator")?;
281 let sample_time = property_scalar(&object, "Ts")?;
282 if !sample_time.is_finite() || sample_time < 0.0 {
283 return Err(step_error_with_detail(
284 &STEP_ERROR_INVALID_MODEL,
285 "tf sample time must be finite and non-negative",
286 ));
287 }
288
289 let numerator = trim_leading_zeros(numerator);
290 let denominator = trim_leading_zeros(denominator);
291 if denominator.is_empty() {
292 return Err(step_error_with_detail(
293 &STEP_ERROR_INVALID_MODEL,
294 "denominator coefficients cannot be empty",
295 ));
296 }
297 if denominator[0].abs() <= EPS {
298 return Err(step_error_with_detail(
299 &STEP_ERROR_INVALID_MODEL,
300 "leading denominator coefficient must be non-zero",
301 ));
302 }
303 if numerator.len().saturating_sub(1) > denominator.len().saturating_sub(1) {
304 return Err(step_error_with_detail(
305 &STEP_ERROR_UNSUPPORTED_MODEL,
306 "improper transfer functions are not supported yet",
307 ));
308 }
309
310 Ok(Self {
311 numerator,
312 denominator,
313 sample_time,
314 })
315 }
316
317 fn normalized(&self) -> (Vec<f64>, Vec<f64>) {
318 let leading = self.denominator[0];
319 let den = self
320 .denominator
321 .iter()
322 .map(|value| value / leading)
323 .collect::<Vec<_>>();
324 let num = self
325 .numerator
326 .iter()
327 .map(|value| value / leading)
328 .collect::<Vec<_>>();
329 (num, den)
330 }
331}
332
333fn property_coefficients(object: &ObjectInstance, name: &str) -> BuiltinResult<Vec<f64>> {
334 let value = object.properties.get(name).ok_or_else(|| {
335 step_error_with_detail(
336 &STEP_ERROR_INVALID_MODEL,
337 format!("tf object is missing {name}"),
338 )
339 })?;
340 match value {
341 Value::Tensor(tensor) => Ok(tensor.data.clone()),
342 Value::ComplexTensor(tensor) => real_complex_coefficients(tensor, name),
343 Value::Num(n) => Ok(vec![*n]),
344 Value::Int(i) => Ok(vec![i.to_f64()]),
345 other => Err(step_error_with_detail(
346 &STEP_ERROR_INVALID_MODEL,
347 format!("tf {name} coefficients must be numeric, got {other:?}"),
348 )),
349 }
350}
351
352fn real_complex_coefficients(tensor: &ComplexTensor, name: &str) -> BuiltinResult<Vec<f64>> {
353 let mut out = Vec::with_capacity(tensor.data.len());
354 for &(re, im) in &tensor.data {
355 if im.abs() > EPS {
356 return Err(step_error_with_detail(
357 &STEP_ERROR_UNSUPPORTED_MODEL,
358 format!("complex tf {name} coefficients are not supported yet"),
359 ));
360 }
361 out.push(re);
362 }
363 Ok(out)
364}
365
366fn property_scalar(object: &ObjectInstance, name: &str) -> BuiltinResult<f64> {
367 let value = object.properties.get(name).ok_or_else(|| {
368 step_error_with_detail(
369 &STEP_ERROR_INVALID_MODEL,
370 format!("tf object is missing {name}"),
371 )
372 })?;
373 match value {
374 Value::Num(n) => Ok(*n),
375 Value::Int(i) => Ok(i.to_f64()),
376 other => Err(step_error_with_detail(
377 &STEP_ERROR_INVALID_MODEL,
378 format!("tf {name} property must be a scalar, got {other:?}"),
379 )),
380 }
381}
382
383#[derive(Clone, Debug)]
384enum TimeSpec {
385 Auto,
386 FinalTime(f64),
387 Vector(Vec<f64>),
388}
389
390impl TimeSpec {
391 fn parse(value: Option<&Value>, sample_time: f64) -> BuiltinResult<Self> {
392 let Some(value) = value else {
393 return Ok(Self::Auto);
394 };
395 match value {
396 Value::Num(n) => Self::final_time(*n),
397 Value::Int(i) => Self::final_time(i.to_f64()),
398 Value::Tensor(tensor) => {
399 ensure_time_vector_shape(&tensor.shape)?;
400 if tensor.data.len() == 1 {
401 return Self::final_time(tensor.data[0]);
402 }
403 Self::vector(tensor.data.clone(), sample_time)
404 }
405 other => Err(step_error_with_detail(
406 &STEP_ERROR_INVALID_TIME,
407 format!("time input must be a scalar final time or numeric vector, got {other:?}"),
408 )),
409 }
410 }
411
412 fn final_time(value: f64) -> BuiltinResult<Self> {
413 if !value.is_finite() || value <= 0.0 {
414 return Err(step_error_with_detail(
415 &STEP_ERROR_INVALID_TIME,
416 "final time must be a positive finite scalar",
417 ));
418 }
419 Ok(Self::FinalTime(value))
420 }
421
422 fn vector(values: Vec<f64>, sample_time: f64) -> BuiltinResult<Self> {
423 validate_time_vector(&values)?;
424 if sample_time > 0.0 {
425 for &t in &values {
426 let k = (t / sample_time).round();
427 if (t - k * sample_time).abs() > 1.0e-8 * sample_time.max(1.0) {
428 return Err(step_error_with_detail(
429 &STEP_ERROR_INVALID_TIME,
430 "discrete-time sample vector must align with the model sample time",
431 ));
432 }
433 }
434 }
435 Ok(Self::Vector(values))
436 }
437}
438
439fn ensure_time_vector_shape(shape: &[usize]) -> BuiltinResult<()> {
440 let non_unit = shape.iter().copied().filter(|&dim| dim > 1).count();
441 if non_unit <= 1 {
442 Ok(())
443 } else {
444 Err(step_error_with_detail(
445 &STEP_ERROR_INVALID_TIME,
446 "time input must be a vector",
447 ))
448 }
449}
450
451fn validate_time_vector(values: &[f64]) -> BuiltinResult<()> {
452 if values.is_empty() {
453 return Err(step_error_with_detail(
454 &STEP_ERROR_INVALID_TIME,
455 "time vector must not be empty",
456 ));
457 }
458 let mut previous = None;
459 for &value in values {
460 if !value.is_finite() || value < 0.0 {
461 return Err(step_error_with_detail(
462 &STEP_ERROR_INVALID_TIME,
463 "time vector values must be finite and non-negative",
464 ));
465 }
466 if let Some(prev) = previous {
467 if value < prev {
468 return Err(step_error_with_detail(
469 &STEP_ERROR_INVALID_TIME,
470 "time vector must be nondecreasing",
471 ));
472 }
473 }
474 previous = Some(value);
475 }
476 Ok(())
477}
478
479#[derive(Clone, Debug)]
480struct StepEval {
481 y: Vec<f64>,
482 t: Vec<f64>,
483}
484
485impl StepEval {
486 fn y_value(&self) -> BuiltinResult<Value> {
487 column_tensor(self.y.clone())
488 }
489
490 fn t_value(&self) -> BuiltinResult<Value> {
491 column_tensor(self.t.clone())
492 }
493
494 fn outputs(&self) -> BuiltinResult<Vec<Value>> {
495 Ok(vec![self.y_value()?, self.t_value()?])
496 }
497}
498
499fn evaluate_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
500 if model.sample_time > 0.0 {
501 evaluate_discrete_step(model, time)
502 } else {
503 evaluate_continuous_step(model, time)
504 }
505}
506
507fn evaluate_continuous_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
508 let t = continuous_time_vector(model, time)?;
509 let (num, den) = model.normalized();
510 let response = continuous_response(&num, &den, &t)?;
511 Ok(StepEval { y: response, t })
512}
513
514fn continuous_time_vector(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<Vec<f64>> {
515 match time {
516 TimeSpec::Auto => Ok(linspace(0.0, automatic_final_time(model), DEFAULT_SAMPLES)),
517 TimeSpec::FinalTime(final_time) => Ok(linspace(0.0, final_time, DEFAULT_SAMPLES)),
518 TimeSpec::Vector(values) => Ok(values),
519 }
520}
521
522fn continuous_response(num: &[f64], den: &[f64], t: &[f64]) -> BuiltinResult<Vec<f64>> {
523 let order = den.len() - 1;
524 if order == 0 {
525 let gain = num.last().copied().unwrap_or(0.0) / den[0];
526 return Ok(vec![gain; t.len()]);
527 }
528
529 let mut padded_num = vec![0.0; order + 1 - num.len()];
530 padded_num.extend_from_slice(num);
531 let direct = padded_num[0];
532 let a = &den[1..];
533 let mut c = Vec::with_capacity(order);
534 for state_idx in 0..order {
535 let coeff_idx = order - state_idx;
536 c.push(padded_num[coeff_idx] - direct * den[coeff_idx]);
537 }
538
539 let mut state = vec![0.0; order];
540 let mut current_t = 0.0;
541 let mut response = Vec::with_capacity(t.len());
542 for &target_t in t {
543 if target_t > current_t {
544 integrate_to(&mut state, a, current_t, target_t);
545 current_t = target_t;
546 }
547 response.push(dot(&c, &state) + direct);
548 }
549 Ok(response)
550}
551
552fn integrate_to(state: &mut [f64], a: &[f64], start: f64, end: f64) {
553 let duration = end - start;
554 if duration <= 0.0 {
555 return;
556 }
557 let steps = ((duration / 0.01).ceil() as usize).clamp(1, 10_000);
558 let h = duration / steps as f64;
559 for _ in 0..steps {
560 rk4_step(state, a, h);
561 }
562}
563
564fn rk4_step(state: &mut [f64], a: &[f64], h: f64) {
565 let k1 = derivative(state, a);
566 let s2 = add_scaled(state, &k1, h * 0.5);
567 let k2 = derivative(&s2, a);
568 let s3 = add_scaled(state, &k2, h * 0.5);
569 let k3 = derivative(&s3, a);
570 let s4 = add_scaled(state, &k3, h);
571 let k4 = derivative(&s4, a);
572 for idx in 0..state.len() {
573 state[idx] += h * (k1[idx] + 2.0 * k2[idx] + 2.0 * k3[idx] + k4[idx]) / 6.0;
574 }
575}
576
577fn derivative(state: &[f64], a: &[f64]) -> Vec<f64> {
578 let order = state.len();
579 let mut dx = vec![0.0; order];
580 if order > 1 {
581 dx[..(order - 1)].copy_from_slice(&state[1..order]);
582 }
583 let mut last = 1.0;
584 for state_idx in 0..order {
585 let coeff = a[order - 1 - state_idx];
586 last -= coeff * state[state_idx];
587 }
588 dx[order - 1] = last;
589 dx
590}
591
592fn add_scaled(state: &[f64], delta: &[f64], scale: f64) -> Vec<f64> {
593 state
594 .iter()
595 .zip(delta)
596 .map(|(value, delta)| value + scale * delta)
597 .collect()
598}
599
600fn evaluate_discrete_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
601 let t = discrete_time_vector(model.sample_time, time)?;
602 if t.len() > MAX_DISCRETE_SAMPLES {
603 return Err(step_error_with_detail(
604 &STEP_ERROR_DISCRETE_LIMIT,
605 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
606 ));
607 }
608 let sample_indices = t
609 .iter()
610 .map(|&value| checked_discrete_sample_index(model.sample_time, value))
611 .collect::<BuiltinResult<Vec<_>>>()?;
612 let max_k = sample_indices.iter().copied().max().unwrap_or(0);
613 let count = max_k.checked_add(1).ok_or_else(|| {
614 step_error_with_detail(
615 &STEP_ERROR_DISCRETE_LIMIT,
616 "discrete sample index exceeds platform limits",
617 )
618 })?;
619 let (num, den) = model.normalized();
620 let all_y = discrete_response(&num, &den, count)?;
621 let y = sample_indices
622 .into_iter()
623 .map(|idx| {
624 all_y.get(idx).copied().ok_or_else(|| {
625 step_error_with_detail(
626 &STEP_ERROR_DISCRETE_LIMIT,
627 "discrete sample index exceeds response length",
628 )
629 })
630 })
631 .collect::<BuiltinResult<Vec<_>>>()?;
632 Ok(StepEval { y, t })
633}
634
635fn discrete_time_vector(sample_time: f64, time: TimeSpec) -> BuiltinResult<Vec<f64>> {
636 match time {
637 TimeSpec::Auto => Ok((0..DEFAULT_SAMPLES)
638 .map(|idx| idx as f64 * sample_time)
639 .collect()),
640 TimeSpec::FinalTime(final_time) => {
641 let steps = checked_discrete_sample_steps(sample_time, final_time)?;
642 Ok((0..=steps).map(|idx| idx as f64 * sample_time).collect())
643 }
644 TimeSpec::Vector(values) => Ok(values),
645 }
646}
647
648fn checked_discrete_sample_steps(sample_time: f64, final_time: f64) -> BuiltinResult<usize> {
649 let steps = (final_time / sample_time).floor();
650 if !steps.is_finite() || steps < 0.0 || steps > usize::MAX as f64 {
651 return Err(step_error_with_detail(
652 &STEP_ERROR_DISCRETE_LIMIT,
653 "discrete sample count exceeds platform limits",
654 ));
655 }
656 if steps >= MAX_DISCRETE_SAMPLES as f64 {
657 return Err(step_error_with_detail(
658 &STEP_ERROR_DISCRETE_LIMIT,
659 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
660 ));
661 }
662 Ok(steps as usize)
663}
664
665fn checked_discrete_sample_index(sample_time: f64, time: f64) -> BuiltinResult<usize> {
666 let index = (time / sample_time).round();
667 if !index.is_finite() || index < 0.0 || index > usize::MAX as f64 {
668 return Err(step_error_with_detail(
669 &STEP_ERROR_DISCRETE_LIMIT,
670 "discrete sample index exceeds platform limits",
671 ));
672 }
673 if index >= MAX_DISCRETE_SAMPLES as f64 {
674 return Err(step_error_with_detail(
675 &STEP_ERROR_DISCRETE_LIMIT,
676 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
677 ));
678 }
679 Ok(index as usize)
680}
681
682fn discrete_response(num: &[f64], den: &[f64], count: usize) -> BuiltinResult<Vec<f64>> {
683 let order = den.len() - 1;
684 let mut padded_num = vec![0.0; order + 1 - num.len()];
685 padded_num.extend_from_slice(num);
686 let mut y = vec![0.0; count];
687 for k in 0..count {
688 let mut value = 0.0;
689 for (idx, &coeff) in padded_num.iter().enumerate() {
690 if k >= idx {
691 value += coeff;
692 }
693 }
694 for idx in 1..den.len() {
695 if k >= idx {
696 value -= den[idx] * y[k - idx];
697 }
698 }
699 y[k] = value;
700 }
701 Ok(y)
702}
703
704async fn plot_response(eval: &StepEval) -> BuiltinResult<()> {
705 let t = eval.t_value()?;
706 let y = eval.y_value()?;
707 if let Err(err) = crate::call_builtin_async("plot", &[t, y]).await {
708 if super::is_nonfatal_plot_setup_error(&err) {
709 return Ok(());
710 }
711 return Err(step_error_with_detail(
712 &STEP_ERROR_PLOT_FAILED,
713 err.message(),
714 ));
715 }
716 let _ = crate::call_builtin_async("title", &[Value::from("Step Response")]).await;
717 let _ = crate::call_builtin_async("xlabel", &[Value::from("Time")]).await;
718 let _ = crate::call_builtin_async("ylabel", &[Value::from("Amplitude")]).await;
719 Ok(())
720}
721
722fn automatic_final_time(model: &TransferFunction) -> f64 {
723 let (_, den) = model.normalized();
724 let poles = polynomial_roots(&den).unwrap_or_default();
725 let slowest_decay = poles
726 .iter()
727 .filter_map(|pole| if pole.re < -EPS { Some(-pole.re) } else { None })
728 .fold(f64::INFINITY, f64::min);
729 if slowest_decay.is_finite() && slowest_decay > EPS {
730 (5.0 / slowest_decay).clamp(1.0, 100.0)
731 } else {
732 10.0
733 }
734}
735
736fn polynomial_roots(coeffs: &[f64]) -> BuiltinResult<Vec<Complex64>> {
737 let trimmed = trim_leading_zeros(coeffs.to_vec());
738 if trimmed.len() <= 1 {
739 return Ok(Vec::new());
740 }
741 if trimmed.len() == 2 {
742 return Ok(vec![Complex64::new(-trimmed[1] / trimmed[0], 0.0)]);
743 }
744 let degree = trimmed.len() - 1;
745 let leading = trimmed[0];
746 let mut companion = DMatrix::<Complex64>::zeros(degree, degree);
747 for row in 1..degree {
748 companion[(row, row - 1)] = Complex64::new(1.0, 0.0);
749 }
750 for (idx, coeff) in trimmed.iter().enumerate().skip(1) {
751 companion[(0, idx - 1)] = Complex64::new(-coeff / leading, 0.0);
752 }
753 let eigenvalues = companion.eigenvalues().ok_or_else(|| {
754 step_error_with_detail(
755 &STEP_ERROR_INTERNAL,
756 "failed to compute transfer-function poles",
757 )
758 })?;
759 Ok(eigenvalues.iter().copied().collect())
760}
761
762fn trim_leading_zeros(coeffs: Vec<f64>) -> Vec<f64> {
763 let first_nonzero = coeffs
764 .iter()
765 .position(|value| value.abs() > EPS)
766 .unwrap_or(coeffs.len());
767 coeffs[first_nonzero..].to_vec()
768}
769
770fn linspace(start: f64, end: f64, count: usize) -> Vec<f64> {
771 if count <= 1 {
772 return vec![end];
773 }
774 let step = (end - start) / (count - 1) as f64;
775 (0..count).map(|idx| start + idx as f64 * step).collect()
776}
777
778fn dot(lhs: &[f64], rhs: &[f64]) -> f64 {
779 lhs.iter().zip(rhs).map(|(a, b)| a * b).sum()
780}
781
782fn column_tensor(data: Vec<f64>) -> BuiltinResult<Value> {
783 let rows = data.len();
784 let tensor = Tensor::new(data, vec![rows, 1]).map_err(|err| {
785 step_error_with_detail(
786 &STEP_ERROR_INTERNAL,
787 format!("failed to build response tensor: {err}"),
788 )
789 })?;
790 Ok(Value::Tensor(tensor))
791}
792
793#[cfg(test)]
794mod tests {
795 use super::*;
796 use futures::executor::block_on;
797 use runmat_builtins::{CharArray, ObjectInstance};
798
799 fn tf_object(num: Vec<f64>, den: Vec<f64>, sample_time: f64) -> Value {
800 let mut object = ObjectInstance::new("tf".to_string());
801 object.properties.insert(
802 "Numerator".to_string(),
803 Value::Tensor(Tensor::new(num.clone(), vec![1, num.len()]).unwrap()),
804 );
805 object.properties.insert(
806 "Denominator".to_string(),
807 Value::Tensor(Tensor::new(den.clone(), vec![1, den.len()]).unwrap()),
808 );
809 object.properties.insert(
810 "Variable".to_string(),
811 Value::CharArray(CharArray::new_row(if sample_time > 0.0 {
812 "z"
813 } else {
814 "s"
815 })),
816 );
817 object
818 .properties
819 .insert("Ts".to_string(), Value::Num(sample_time));
820 Value::Object(object)
821 }
822
823 fn run_step(sys: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
824 block_on(step_builtin(sys, rest))
825 }
826
827 fn tensor_data(value: Value) -> Vec<f64> {
828 match value {
829 Value::Tensor(tensor) => tensor.data,
830 other => panic!("expected tensor, got {other:?}"),
831 }
832 }
833
834 #[test]
835 fn step_descriptor_signatures_cover_core_forms() {
836 let labels: Vec<&str> = STEP_DESCRIPTOR
837 .signatures
838 .iter()
839 .map(|sig| sig.label)
840 .collect();
841 assert!(labels.contains(&"y = step(sys)"));
842 assert!(labels.contains(&"y = step(sys, tFinal)"));
843 assert!(labels.contains(&"y = step(sys, t)"));
844 assert!(labels.contains(&"[y,t] = step(sys)"));
845 assert!(labels.contains(&"[y,t] = step(sys, tFinal)"));
846 assert!(labels.contains(&"[y,t] = step(sys, t)"));
847 }
848
849 #[test]
850 fn first_order_continuous_response_matches_closed_form_for_explicit_time() {
851 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
852 let time = Value::Tensor(Tensor::new(vec![0.0, 0.5, 1.0, 2.0], vec![1, 4]).unwrap());
853 let y = tensor_data(run_step(sys, vec![time]).expect("step"));
854 for (actual, t) in y.iter().zip([0.0_f64, 0.5, 1.0, 2.0]) {
855 let expected = 1.0 - (-t).exp();
856 assert!(
857 (actual - expected).abs() < 1.0e-5,
858 "t={t} actual={actual} expected={expected}"
859 );
860 }
861 }
862
863 #[test]
864 fn multi_output_returns_y_then_time() {
865 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
866 let _guard = crate::output_count::push_output_count(Some(2));
867 let time = Value::Tensor(Tensor::new(vec![0.0, 1.0], vec![2, 1]).unwrap());
868 let result = run_step(sys, vec![time]).expect("step");
869 let Value::OutputList(outputs) = result else {
870 panic!("expected output list");
871 };
872 assert_eq!(outputs.len(), 2);
873 assert_eq!(tensor_data(outputs[1].clone()), vec![0.0, 1.0]);
874 let y = tensor_data(outputs[0].clone());
875 assert!((y[1] - (1.0 - (-1.0_f64).exp())).abs() < 1.0e-5);
876 }
877
878 #[test]
879 fn single_requested_output_returns_only_response() {
880 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
881 let _guard = crate::output_count::push_output_count(Some(1));
882 let time = Value::Tensor(Tensor::new(vec![0.0, 1.0], vec![1, 2]).unwrap());
883 let result = run_step(sys, vec![time]).expect("step");
884 let Value::OutputList(outputs) = result else {
885 panic!("expected output list");
886 };
887 assert_eq!(outputs.len(), 1);
888 let y = tensor_data(outputs[0].clone());
889 assert_eq!(y.len(), 2);
890 assert!((y[1] - (1.0 - (-1.0_f64).exp())).abs() < 1.0e-5);
891 }
892
893 #[test]
894 fn scalar_final_time_generates_column_time_vector_ending_at_final_time() {
895 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
896 let _guard = crate::output_count::push_output_count(Some(2));
897 let result = run_step(sys, vec![Value::Num(5.0)]).expect("step");
898 let Value::OutputList(outputs) = result else {
899 panic!("expected output list");
900 };
901 let t = tensor_data(outputs[1].clone());
902 assert_eq!(t.len(), DEFAULT_SAMPLES);
903 assert_eq!(t[0], 0.0);
904 assert!((t[t.len() - 1] - 5.0).abs() < 1.0e-12);
905 }
906
907 #[test]
908 fn discrete_response_uses_sample_time_grid() {
909 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 0.1);
910 let time = Value::Tensor(Tensor::new(vec![0.0, 0.1, 0.2], vec![1, 3]).unwrap());
911 let y = tensor_data(run_step(sys, vec![time]).expect("step"));
912 assert_eq!(y, vec![0.0, 1.0, 1.5]);
913 }
914
915 #[test]
916 fn discrete_final_time_rejects_excessive_sample_count() {
917 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0e-6);
918 let err = run_step(sys, vec![Value::Num(2.0)]).expect_err("should fail");
919 assert!(err.message().contains("more than 1000000 samples"));
920 assert_eq!(err.identifier(), STEP_ERROR_DISCRETE_LIMIT.identifier);
921 }
922
923 #[test]
924 fn discrete_time_vector_rejects_excessive_sample_index() {
925 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0);
926 let time =
927 Value::Tensor(Tensor::new(vec![0.0, MAX_DISCRETE_SAMPLES as f64], vec![1, 2]).unwrap());
928 let err = run_step(sys, vec![time]).expect_err("should fail");
929 assert!(err.message().contains("more than 1000000 samples"));
930 }
931
932 #[test]
933 fn rejects_non_tf_input() {
934 let err = run_step(Value::Num(1.0), Vec::new()).expect_err("expected error");
935 assert!(err.message().contains("expected a tf object"));
936 assert_eq!(err.identifier(), STEP_ERROR_INVALID_MODEL.identifier);
937 }
938}