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 is_statement_form_call() {
218 plot_multiple_step_responses(sys, rest).await?;
219 return Ok(Value::OutputList(Vec::new()));
220 }
221
222 if rest.len() > 1 {
223 return Err(step_error_with_detail(
224 &STEP_ERROR_INVALID_ARGUMENT,
225 "expected step(sys), step(sys, tFinal), or step(sys, t)",
226 ));
227 }
228
229 let sys = crate::gather_if_needed_async(&sys).await?;
230 let mut rest_host = Vec::with_capacity(rest.len());
231 for arg in &rest {
232 rest_host.push(crate::gather_if_needed_async(arg).await?);
233 }
234 let model = TransferFunction::from_value(sys)?;
235 let time = TimeSpec::parse(rest_host.first(), model.sample_time)?;
236 let eval = evaluate_step(&model, time)?;
237
238 if crate::output_context::requested_output_count() == Some(0)
239 && crate::output_count::current_output_count().is_none()
240 {
241 plot_response(&eval).await?;
242 return Ok(Value::OutputList(Vec::new()));
243 }
244
245 if let Some(out_count) = crate::output_count::current_output_count() {
246 if out_count == 0 {
247 plot_response(&eval).await?;
248 return Ok(Value::OutputList(Vec::new()));
249 }
250 if out_count == 1 {
251 return Ok(Value::OutputList(vec![eval.y_value()?]));
252 }
253 return Ok(crate::output_count::output_list_with_padding(
254 out_count,
255 eval.outputs()?,
256 ));
257 }
258
259 eval.y_value()
260}
261
262fn is_statement_form_call() -> bool {
263 matches!(crate::output_count::current_output_count(), Some(0))
264 || (crate::output_context::requested_output_count() == Some(0)
265 && crate::output_count::current_output_count().is_none())
266}
267
268async fn plot_multiple_step_responses(first_sys: Value, rest: Vec<Value>) -> BuiltinResult<()> {
269 let first_sys = crate::gather_if_needed_async(&first_sys).await?;
270 let mut systems = vec![(first_sys, None)];
271 let mut time_arg = None;
272 for arg in rest {
273 let gathered = crate::gather_if_needed_async(&arg).await?;
274 if is_plot_style_arg(&gathered) {
275 if let Some((_, style)) = systems.last_mut() {
276 if style.is_some() {
277 return Err(step_error_with_detail(
278 &STEP_ERROR_INVALID_ARGUMENT,
279 "only one style argument is supported per system",
280 ));
281 }
282 *style = Some(gathered);
283 continue;
284 }
285 continue;
286 }
287 if is_tf_object(&gathered) {
288 if time_arg.is_some() {
289 return Err(step_error_with_detail(
290 &STEP_ERROR_INVALID_ARGUMENT,
291 "time argument must follow all systems in statement-form step plots",
292 ));
293 }
294 systems.push((gathered, None));
295 continue;
296 }
297 if time_arg.is_none() {
298 time_arg = Some(gathered);
299 continue;
300 }
301 return Err(step_error_with_detail(
302 &STEP_ERROR_INVALID_ARGUMENT,
303 "unsupported statement-form step plot argument",
304 ));
305 }
306
307 if systems.is_empty() {
308 return Err(step_error_with_detail(
309 &STEP_ERROR_INVALID_ARGUMENT,
310 "at least one system is required",
311 ));
312 }
313
314 let mut first = true;
315 let mut hold_enabled = false;
316 let result: BuiltinResult<()> = async {
317 for (system, style) in systems {
318 let model = TransferFunction::from_value(system)?;
319 let time = TimeSpec::parse(time_arg.as_ref(), model.sample_time)?;
320 let eval = evaluate_step(&model, time)?;
321 plot_response_with_style(&eval, style.as_ref()).await?;
322 if first {
323 first = false;
324 let _ = crate::call_builtin_async("hold", &[Value::from("on")]).await;
325 hold_enabled = true;
326 }
327 }
328 Ok(())
329 }
330 .await;
331 if hold_enabled {
332 let _ = crate::call_builtin_async("hold", &[Value::from("off")]).await;
333 }
334 result
335}
336
337fn is_plot_style_arg(value: &Value) -> bool {
338 matches!(
339 value,
340 Value::String(_) | Value::StringArray(_) | Value::CharArray(_)
341 )
342}
343
344fn is_tf_object(value: &Value) -> bool {
345 matches!(value, Value::Object(object) if object.is_class("tf"))
346}
347
348#[derive(Clone, Debug)]
349struct TransferFunction {
350 numerator: Vec<f64>,
351 denominator: Vec<f64>,
352 sample_time: f64,
353}
354
355impl TransferFunction {
356 fn from_value(value: Value) -> BuiltinResult<Self> {
357 let Value::Object(object) = value else {
358 return Err(step_error_with_detail(
359 &STEP_ERROR_INVALID_MODEL,
360 "expected a tf object",
361 ));
362 };
363 if !object.is_class("tf") {
364 return Err(step_error_with_detail(
365 &STEP_ERROR_INVALID_MODEL,
366 format!("expected a tf object, got {}", object.class_name),
367 ));
368 }
369
370 let numerator = property_coefficients(&object, "Numerator")?;
371 let denominator = property_coefficients(&object, "Denominator")?;
372 let sample_time = property_scalar(&object, "Ts")?;
373 if !sample_time.is_finite() || sample_time < 0.0 {
374 return Err(step_error_with_detail(
375 &STEP_ERROR_INVALID_MODEL,
376 "tf sample time must be finite and non-negative",
377 ));
378 }
379
380 let numerator = trim_leading_zeros(numerator);
381 let denominator = trim_leading_zeros(denominator);
382 if denominator.is_empty() {
383 return Err(step_error_with_detail(
384 &STEP_ERROR_INVALID_MODEL,
385 "denominator coefficients cannot be empty",
386 ));
387 }
388 if denominator[0].abs() <= EPS {
389 return Err(step_error_with_detail(
390 &STEP_ERROR_INVALID_MODEL,
391 "leading denominator coefficient must be non-zero",
392 ));
393 }
394 if numerator.len().saturating_sub(1) > denominator.len().saturating_sub(1) {
395 return Err(step_error_with_detail(
396 &STEP_ERROR_UNSUPPORTED_MODEL,
397 "improper transfer functions are not supported yet",
398 ));
399 }
400
401 Ok(Self {
402 numerator,
403 denominator,
404 sample_time,
405 })
406 }
407
408 fn normalized(&self) -> (Vec<f64>, Vec<f64>) {
409 let leading = self.denominator[0];
410 let den = self
411 .denominator
412 .iter()
413 .map(|value| value / leading)
414 .collect::<Vec<_>>();
415 let num = self
416 .numerator
417 .iter()
418 .map(|value| value / leading)
419 .collect::<Vec<_>>();
420 (num, den)
421 }
422}
423
424fn property_coefficients(object: &ObjectInstance, name: &str) -> BuiltinResult<Vec<f64>> {
425 let value = object.properties.get(name).ok_or_else(|| {
426 step_error_with_detail(
427 &STEP_ERROR_INVALID_MODEL,
428 format!("tf object is missing {name}"),
429 )
430 })?;
431 match value {
432 Value::Tensor(tensor) => Ok(tensor.data.clone()),
433 Value::ComplexTensor(tensor) => real_complex_coefficients(tensor, name),
434 Value::Num(n) => Ok(vec![*n]),
435 Value::Int(i) => Ok(vec![i.to_f64()]),
436 other => Err(step_error_with_detail(
437 &STEP_ERROR_INVALID_MODEL,
438 format!("tf {name} coefficients must be numeric, got {other:?}"),
439 )),
440 }
441}
442
443fn real_complex_coefficients(tensor: &ComplexTensor, name: &str) -> BuiltinResult<Vec<f64>> {
444 let mut out = Vec::with_capacity(tensor.data.len());
445 for &(re, im) in &tensor.data {
446 if im.abs() > EPS {
447 return Err(step_error_with_detail(
448 &STEP_ERROR_UNSUPPORTED_MODEL,
449 format!("complex tf {name} coefficients are not supported yet"),
450 ));
451 }
452 out.push(re);
453 }
454 Ok(out)
455}
456
457fn property_scalar(object: &ObjectInstance, name: &str) -> BuiltinResult<f64> {
458 let value = object.properties.get(name).ok_or_else(|| {
459 step_error_with_detail(
460 &STEP_ERROR_INVALID_MODEL,
461 format!("tf object is missing {name}"),
462 )
463 })?;
464 match value {
465 Value::Num(n) => Ok(*n),
466 Value::Int(i) => Ok(i.to_f64()),
467 other => Err(step_error_with_detail(
468 &STEP_ERROR_INVALID_MODEL,
469 format!("tf {name} property must be a scalar, got {other:?}"),
470 )),
471 }
472}
473
474#[derive(Clone, Debug)]
475enum TimeSpec {
476 Auto,
477 FinalTime(f64),
478 Vector(Vec<f64>),
479}
480
481impl TimeSpec {
482 fn parse(value: Option<&Value>, sample_time: f64) -> BuiltinResult<Self> {
483 let Some(value) = value else {
484 return Ok(Self::Auto);
485 };
486 match value {
487 Value::Num(n) => Self::final_time(*n),
488 Value::Int(i) => Self::final_time(i.to_f64()),
489 Value::Tensor(tensor) => {
490 ensure_time_vector_shape(&tensor.shape)?;
491 if tensor.data.len() == 1 {
492 return Self::final_time(tensor.data[0]);
493 }
494 Self::vector(tensor.data.clone(), sample_time)
495 }
496 other => Err(step_error_with_detail(
497 &STEP_ERROR_INVALID_TIME,
498 format!("time input must be a scalar final time or numeric vector, got {other:?}"),
499 )),
500 }
501 }
502
503 fn final_time(value: f64) -> BuiltinResult<Self> {
504 if !value.is_finite() || value <= 0.0 {
505 return Err(step_error_with_detail(
506 &STEP_ERROR_INVALID_TIME,
507 "final time must be a positive finite scalar",
508 ));
509 }
510 Ok(Self::FinalTime(value))
511 }
512
513 fn vector(values: Vec<f64>, sample_time: f64) -> BuiltinResult<Self> {
514 validate_time_vector(&values)?;
515 if sample_time > 0.0 {
516 for &t in &values {
517 let k = (t / sample_time).round();
518 if (t - k * sample_time).abs() > 1.0e-8 * sample_time.max(1.0) {
519 return Err(step_error_with_detail(
520 &STEP_ERROR_INVALID_TIME,
521 "discrete-time sample vector must align with the model sample time",
522 ));
523 }
524 }
525 }
526 Ok(Self::Vector(values))
527 }
528}
529
530fn ensure_time_vector_shape(shape: &[usize]) -> BuiltinResult<()> {
531 let non_unit = shape.iter().copied().filter(|&dim| dim > 1).count();
532 if non_unit <= 1 {
533 Ok(())
534 } else {
535 Err(step_error_with_detail(
536 &STEP_ERROR_INVALID_TIME,
537 "time input must be a vector",
538 ))
539 }
540}
541
542fn validate_time_vector(values: &[f64]) -> BuiltinResult<()> {
543 if values.is_empty() {
544 return Err(step_error_with_detail(
545 &STEP_ERROR_INVALID_TIME,
546 "time vector must not be empty",
547 ));
548 }
549 let mut previous = None;
550 for &value in values {
551 if !value.is_finite() || value < 0.0 {
552 return Err(step_error_with_detail(
553 &STEP_ERROR_INVALID_TIME,
554 "time vector values must be finite and non-negative",
555 ));
556 }
557 if let Some(prev) = previous {
558 if value < prev {
559 return Err(step_error_with_detail(
560 &STEP_ERROR_INVALID_TIME,
561 "time vector must be nondecreasing",
562 ));
563 }
564 }
565 previous = Some(value);
566 }
567 Ok(())
568}
569
570#[derive(Clone, Debug)]
571struct StepEval {
572 y: Vec<f64>,
573 t: Vec<f64>,
574}
575
576impl StepEval {
577 fn y_value(&self) -> BuiltinResult<Value> {
578 column_tensor(self.y.clone())
579 }
580
581 fn t_value(&self) -> BuiltinResult<Value> {
582 column_tensor(self.t.clone())
583 }
584
585 fn outputs(&self) -> BuiltinResult<Vec<Value>> {
586 Ok(vec![self.y_value()?, self.t_value()?])
587 }
588}
589
590fn evaluate_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
591 if model.sample_time > 0.0 {
592 evaluate_discrete_step(model, time)
593 } else {
594 evaluate_continuous_step(model, time)
595 }
596}
597
598fn evaluate_continuous_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
599 let t = continuous_time_vector(model, time)?;
600 let (num, den) = model.normalized();
601 let response = continuous_response(&num, &den, &t)?;
602 Ok(StepEval { y: response, t })
603}
604
605fn continuous_time_vector(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<Vec<f64>> {
606 match time {
607 TimeSpec::Auto => Ok(linspace(0.0, automatic_final_time(model), DEFAULT_SAMPLES)),
608 TimeSpec::FinalTime(final_time) => Ok(linspace(0.0, final_time, DEFAULT_SAMPLES)),
609 TimeSpec::Vector(values) => Ok(values),
610 }
611}
612
613fn continuous_response(num: &[f64], den: &[f64], t: &[f64]) -> BuiltinResult<Vec<f64>> {
614 let order = den.len() - 1;
615 if order == 0 {
616 let gain = num.last().copied().unwrap_or(0.0) / den[0];
617 return Ok(vec![gain; t.len()]);
618 }
619
620 let mut padded_num = vec![0.0; order + 1 - num.len()];
621 padded_num.extend_from_slice(num);
622 let direct = padded_num[0];
623 let a = &den[1..];
624 let mut c = Vec::with_capacity(order);
625 for state_idx in 0..order {
626 let coeff_idx = order - state_idx;
627 c.push(padded_num[coeff_idx] - direct * den[coeff_idx]);
628 }
629
630 let mut state = vec![0.0; order];
631 let mut current_t = 0.0;
632 let mut response = Vec::with_capacity(t.len());
633 for &target_t in t {
634 if target_t > current_t {
635 integrate_to(&mut state, a, current_t, target_t);
636 current_t = target_t;
637 }
638 response.push(dot(&c, &state) + direct);
639 }
640 Ok(response)
641}
642
643fn integrate_to(state: &mut [f64], a: &[f64], start: f64, end: f64) {
644 let duration = end - start;
645 if duration <= 0.0 {
646 return;
647 }
648 let steps = ((duration / 0.01).ceil() as usize).clamp(1, 10_000);
649 let h = duration / steps as f64;
650 for _ in 0..steps {
651 rk4_step(state, a, h);
652 }
653}
654
655fn rk4_step(state: &mut [f64], a: &[f64], h: f64) {
656 let k1 = derivative(state, a);
657 let s2 = add_scaled(state, &k1, h * 0.5);
658 let k2 = derivative(&s2, a);
659 let s3 = add_scaled(state, &k2, h * 0.5);
660 let k3 = derivative(&s3, a);
661 let s4 = add_scaled(state, &k3, h);
662 let k4 = derivative(&s4, a);
663 for idx in 0..state.len() {
664 state[idx] += h * (k1[idx] + 2.0 * k2[idx] + 2.0 * k3[idx] + k4[idx]) / 6.0;
665 }
666}
667
668fn derivative(state: &[f64], a: &[f64]) -> Vec<f64> {
669 let order = state.len();
670 let mut dx = vec![0.0; order];
671 if order > 1 {
672 dx[..(order - 1)].copy_from_slice(&state[1..order]);
673 }
674 let mut last = 1.0;
675 for state_idx in 0..order {
676 let coeff = a[order - 1 - state_idx];
677 last -= coeff * state[state_idx];
678 }
679 dx[order - 1] = last;
680 dx
681}
682
683fn add_scaled(state: &[f64], delta: &[f64], scale: f64) -> Vec<f64> {
684 state
685 .iter()
686 .zip(delta)
687 .map(|(value, delta)| value + scale * delta)
688 .collect()
689}
690
691fn evaluate_discrete_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
692 let t = discrete_time_vector(model.sample_time, time)?;
693 if t.len() > MAX_DISCRETE_SAMPLES {
694 return Err(step_error_with_detail(
695 &STEP_ERROR_DISCRETE_LIMIT,
696 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
697 ));
698 }
699 let sample_indices = t
700 .iter()
701 .map(|&value| checked_discrete_sample_index(model.sample_time, value))
702 .collect::<BuiltinResult<Vec<_>>>()?;
703 let max_k = sample_indices.iter().copied().max().unwrap_or(0);
704 let count = max_k.checked_add(1).ok_or_else(|| {
705 step_error_with_detail(
706 &STEP_ERROR_DISCRETE_LIMIT,
707 "discrete sample index exceeds platform limits",
708 )
709 })?;
710 let (num, den) = model.normalized();
711 let all_y = discrete_response(&num, &den, count)?;
712 let y = sample_indices
713 .into_iter()
714 .map(|idx| {
715 all_y.get(idx).copied().ok_or_else(|| {
716 step_error_with_detail(
717 &STEP_ERROR_DISCRETE_LIMIT,
718 "discrete sample index exceeds response length",
719 )
720 })
721 })
722 .collect::<BuiltinResult<Vec<_>>>()?;
723 Ok(StepEval { y, t })
724}
725
726fn discrete_time_vector(sample_time: f64, time: TimeSpec) -> BuiltinResult<Vec<f64>> {
727 match time {
728 TimeSpec::Auto => Ok((0..DEFAULT_SAMPLES)
729 .map(|idx| idx as f64 * sample_time)
730 .collect()),
731 TimeSpec::FinalTime(final_time) => {
732 let steps = checked_discrete_sample_steps(sample_time, final_time)?;
733 Ok((0..=steps).map(|idx| idx as f64 * sample_time).collect())
734 }
735 TimeSpec::Vector(values) => Ok(values),
736 }
737}
738
739fn checked_discrete_sample_steps(sample_time: f64, final_time: f64) -> BuiltinResult<usize> {
740 let steps = (final_time / sample_time).floor();
741 if !steps.is_finite() || steps < 0.0 || steps > usize::MAX as f64 {
742 return Err(step_error_with_detail(
743 &STEP_ERROR_DISCRETE_LIMIT,
744 "discrete sample count exceeds platform limits",
745 ));
746 }
747 if steps >= MAX_DISCRETE_SAMPLES as f64 {
748 return Err(step_error_with_detail(
749 &STEP_ERROR_DISCRETE_LIMIT,
750 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
751 ));
752 }
753 Ok(steps as usize)
754}
755
756fn checked_discrete_sample_index(sample_time: f64, time: f64) -> BuiltinResult<usize> {
757 let index = (time / sample_time).round();
758 if !index.is_finite() || index < 0.0 || index > usize::MAX as f64 {
759 return Err(step_error_with_detail(
760 &STEP_ERROR_DISCRETE_LIMIT,
761 "discrete sample index exceeds platform limits",
762 ));
763 }
764 if index >= MAX_DISCRETE_SAMPLES as f64 {
765 return Err(step_error_with_detail(
766 &STEP_ERROR_DISCRETE_LIMIT,
767 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
768 ));
769 }
770 Ok(index as usize)
771}
772
773fn discrete_response(num: &[f64], den: &[f64], count: usize) -> BuiltinResult<Vec<f64>> {
774 let order = den.len() - 1;
775 let mut padded_num = vec![0.0; order + 1 - num.len()];
776 padded_num.extend_from_slice(num);
777 let mut y = vec![0.0; count];
778 for k in 0..count {
779 let mut value = 0.0;
780 for (idx, &coeff) in padded_num.iter().enumerate() {
781 if k >= idx {
782 value += coeff;
783 }
784 }
785 for idx in 1..den.len() {
786 if k >= idx {
787 value -= den[idx] * y[k - idx];
788 }
789 }
790 y[k] = value;
791 }
792 Ok(y)
793}
794
795async fn plot_response(eval: &StepEval) -> BuiltinResult<()> {
796 plot_response_with_style(eval, None).await
797}
798
799async fn plot_response_with_style(eval: &StepEval, style: Option<&Value>) -> BuiltinResult<()> {
800 let t = eval.t_value()?;
801 let y = eval.y_value()?;
802 let args = if let Some(style) = style {
803 vec![t, y, style.clone()]
804 } else {
805 vec![t, y]
806 };
807 if let Err(err) = crate::call_builtin_async("plot", &args).await {
808 if super::is_nonfatal_plot_setup_error(&err) {
809 return Ok(());
810 }
811 return Err(step_error_with_detail(
812 &STEP_ERROR_PLOT_FAILED,
813 err.message(),
814 ));
815 }
816 let _ = crate::call_builtin_async("title", &[Value::from("Step Response")]).await;
817 let _ = crate::call_builtin_async("xlabel", &[Value::from("Time")]).await;
818 let _ = crate::call_builtin_async("ylabel", &[Value::from("Amplitude")]).await;
819 Ok(())
820}
821
822fn automatic_final_time(model: &TransferFunction) -> f64 {
823 let (_, den) = model.normalized();
824 let poles = polynomial_roots(&den).unwrap_or_default();
825 let slowest_decay = poles
826 .iter()
827 .filter_map(|pole| if pole.re < -EPS { Some(-pole.re) } else { None })
828 .fold(f64::INFINITY, f64::min);
829 if slowest_decay.is_finite() && slowest_decay > EPS {
830 (5.0 / slowest_decay).clamp(1.0, 100.0)
831 } else {
832 10.0
833 }
834}
835
836fn polynomial_roots(coeffs: &[f64]) -> BuiltinResult<Vec<Complex64>> {
837 let trimmed = trim_leading_zeros(coeffs.to_vec());
838 if trimmed.len() <= 1 {
839 return Ok(Vec::new());
840 }
841 if trimmed.len() == 2 {
842 return Ok(vec![Complex64::new(-trimmed[1] / trimmed[0], 0.0)]);
843 }
844 let degree = trimmed.len() - 1;
845 let leading = trimmed[0];
846 let mut companion = DMatrix::<Complex64>::zeros(degree, degree);
847 for row in 1..degree {
848 companion[(row, row - 1)] = Complex64::new(1.0, 0.0);
849 }
850 for (idx, coeff) in trimmed.iter().enumerate().skip(1) {
851 companion[(0, idx - 1)] = Complex64::new(-coeff / leading, 0.0);
852 }
853 let eigenvalues = companion.eigenvalues().ok_or_else(|| {
854 step_error_with_detail(
855 &STEP_ERROR_INTERNAL,
856 "failed to compute transfer-function poles",
857 )
858 })?;
859 Ok(eigenvalues.iter().copied().collect())
860}
861
862fn trim_leading_zeros(coeffs: Vec<f64>) -> Vec<f64> {
863 let first_nonzero = coeffs
864 .iter()
865 .position(|value| value.abs() > EPS)
866 .unwrap_or(coeffs.len());
867 coeffs[first_nonzero..].to_vec()
868}
869
870fn linspace(start: f64, end: f64, count: usize) -> Vec<f64> {
871 if count <= 1 {
872 return vec![end];
873 }
874 let step = (end - start) / (count - 1) as f64;
875 (0..count).map(|idx| start + idx as f64 * step).collect()
876}
877
878fn dot(lhs: &[f64], rhs: &[f64]) -> f64 {
879 lhs.iter().zip(rhs).map(|(a, b)| a * b).sum()
880}
881
882fn column_tensor(data: Vec<f64>) -> BuiltinResult<Value> {
883 let rows = data.len();
884 let tensor = Tensor::new(data, vec![rows, 1]).map_err(|err| {
885 step_error_with_detail(
886 &STEP_ERROR_INTERNAL,
887 format!("failed to build response tensor: {err}"),
888 )
889 })?;
890 Ok(Value::Tensor(tensor))
891}
892
893#[cfg(test)]
894mod tests {
895 use super::*;
896 use futures::executor::block_on;
897 use runmat_builtins::{CharArray, ObjectInstance};
898
899 fn tf_object(num: Vec<f64>, den: Vec<f64>, sample_time: f64) -> Value {
900 let mut object = ObjectInstance::new("tf".to_string());
901 object.properties.insert(
902 "Numerator".to_string(),
903 Value::Tensor(Tensor::new(num.clone(), vec![1, num.len()]).unwrap()),
904 );
905 object.properties.insert(
906 "Denominator".to_string(),
907 Value::Tensor(Tensor::new(den.clone(), vec![1, den.len()]).unwrap()),
908 );
909 object.properties.insert(
910 "Variable".to_string(),
911 Value::CharArray(CharArray::new_row(if sample_time > 0.0 {
912 "z"
913 } else {
914 "s"
915 })),
916 );
917 object
918 .properties
919 .insert("Ts".to_string(), Value::Num(sample_time));
920 Value::Object(object)
921 }
922
923 fn run_step(sys: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
924 block_on(step_builtin(sys, rest))
925 }
926
927 fn tensor_data(value: Value) -> Vec<f64> {
928 match value {
929 Value::Tensor(tensor) => tensor.data,
930 other => panic!("expected tensor, got {other:?}"),
931 }
932 }
933
934 #[test]
935 fn step_descriptor_signatures_cover_core_forms() {
936 let labels: Vec<&str> = STEP_DESCRIPTOR
937 .signatures
938 .iter()
939 .map(|sig| sig.label)
940 .collect();
941 assert!(labels.contains(&"y = step(sys)"));
942 assert!(labels.contains(&"y = step(sys, tFinal)"));
943 assert!(labels.contains(&"y = step(sys, t)"));
944 assert!(labels.contains(&"[y,t] = step(sys)"));
945 assert!(labels.contains(&"[y,t] = step(sys, tFinal)"));
946 assert!(labels.contains(&"[y,t] = step(sys, t)"));
947 }
948
949 #[test]
950 fn first_order_continuous_response_matches_closed_form_for_explicit_time() {
951 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
952 let time = Value::Tensor(Tensor::new(vec![0.0, 0.5, 1.0, 2.0], vec![1, 4]).unwrap());
953 let y = tensor_data(run_step(sys, vec![time]).expect("step"));
954 for (actual, t) in y.iter().zip([0.0_f64, 0.5, 1.0, 2.0]) {
955 let expected = 1.0 - (-t).exp();
956 assert!(
957 (actual - expected).abs() < 1.0e-5,
958 "t={t} actual={actual} expected={expected}"
959 );
960 }
961 }
962
963 #[test]
964 fn multi_output_returns_y_then_time() {
965 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
966 let _guard = crate::output_count::push_output_count(Some(2));
967 let time = Value::Tensor(Tensor::new(vec![0.0, 1.0], vec![2, 1]).unwrap());
968 let result = run_step(sys, vec![time]).expect("step");
969 let Value::OutputList(outputs) = result else {
970 panic!("expected output list");
971 };
972 assert_eq!(outputs.len(), 2);
973 assert_eq!(tensor_data(outputs[1].clone()), vec![0.0, 1.0]);
974 let y = tensor_data(outputs[0].clone());
975 assert!((y[1] - (1.0 - (-1.0_f64).exp())).abs() < 1.0e-5);
976 }
977
978 #[test]
979 fn single_requested_output_returns_only_response() {
980 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
981 let _guard = crate::output_count::push_output_count(Some(1));
982 let time = Value::Tensor(Tensor::new(vec![0.0, 1.0], vec![1, 2]).unwrap());
983 let result = run_step(sys, vec![time]).expect("step");
984 let Value::OutputList(outputs) = result else {
985 panic!("expected output list");
986 };
987 assert_eq!(outputs.len(), 1);
988 let y = tensor_data(outputs[0].clone());
989 assert_eq!(y.len(), 2);
990 assert!((y[1] - (1.0 - (-1.0_f64).exp())).abs() < 1.0e-5);
991 }
992
993 #[test]
994 fn scalar_final_time_generates_column_time_vector_ending_at_final_time() {
995 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
996 let _guard = crate::output_count::push_output_count(Some(2));
997 let result = run_step(sys, vec![Value::Num(5.0)]).expect("step");
998 let Value::OutputList(outputs) = result else {
999 panic!("expected output list");
1000 };
1001 let t = tensor_data(outputs[1].clone());
1002 assert_eq!(t.len(), DEFAULT_SAMPLES);
1003 assert_eq!(t[0], 0.0);
1004 assert!((t[t.len() - 1] - 5.0).abs() < 1.0e-12);
1005 }
1006
1007 #[test]
1008 fn discrete_response_uses_sample_time_grid() {
1009 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 0.1);
1010 let time = Value::Tensor(Tensor::new(vec![0.0, 0.1, 0.2], vec![1, 3]).unwrap());
1011 let y = tensor_data(run_step(sys, vec![time]).expect("step"));
1012 assert_eq!(y, vec![0.0, 1.0, 1.5]);
1013 }
1014
1015 #[test]
1016 fn discrete_final_time_rejects_excessive_sample_count() {
1017 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0e-6);
1018 let err = run_step(sys, vec![Value::Num(2.0)]).expect_err("should fail");
1019 assert!(err.message().contains("more than 1000000 samples"));
1020 assert_eq!(err.identifier(), STEP_ERROR_DISCRETE_LIMIT.identifier);
1021 }
1022
1023 #[test]
1024 fn discrete_time_vector_rejects_excessive_sample_index() {
1025 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0);
1026 let time =
1027 Value::Tensor(Tensor::new(vec![0.0, MAX_DISCRETE_SAMPLES as f64], vec![1, 2]).unwrap());
1028 let err = run_step(sys, vec![time]).expect_err("should fail");
1029 assert!(err.message().contains("more than 1000000 samples"));
1030 }
1031
1032 #[test]
1033 fn rejects_non_tf_input() {
1034 let err = run_step(Value::Num(1.0), Vec::new()).expect_err("expected error");
1035 assert!(err.message().contains("expected a tf object"));
1036 assert_eq!(err.identifier(), STEP_ERROR_INVALID_MODEL.identifier);
1037 }
1038}