1use runmat_builtins::{
4 BuiltinCompletionPolicy, BuiltinDescriptor, BuiltinErrorDescriptor, BuiltinOutputMode,
5 BuiltinParamArity, BuiltinParamDescriptor, BuiltinParamType, BuiltinSignatureDescriptor,
6 StructValue, Tensor, Value,
7};
8use runmat_macros::runtime_builtin;
9
10use crate::builtins::common::spec::{
11 BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
12 ReductionNaN, ResidencyPolicy, ShapeRequirements,
13};
14use crate::builtins::control::tf_model::{control_error, scalar_f64, scalar_text, EPS};
15use crate::builtins::control::type_resolvers::stepinfo_type;
16use crate::{BuiltinResult, RuntimeError};
17
18const BUILTIN_NAME: &str = "stepinfo";
19const DEFAULT_SETTLING_THRESHOLD: f64 = 0.02;
20
21const STEPINFO_OUTPUT: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
22 name: "info",
23 ty: BuiltinParamType::Any,
24 arity: BuiltinParamArity::Required,
25 default: None,
26 description: "Step-response metrics struct.",
27}];
28const STEPINFO_INPUT_SYS: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
29 name: "sys",
30 ty: BuiltinParamType::Any,
31 arity: BuiltinParamArity::Required,
32 default: None,
33 description: "SISO tf model.",
34}];
35const STEPINFO_PARAM_Y: BuiltinParamDescriptor = BuiltinParamDescriptor {
36 name: "y",
37 ty: BuiltinParamType::NumericArray,
38 arity: BuiltinParamArity::Required,
39 default: None,
40 description: "Response samples.",
41};
42const STEPINFO_PARAM_T: BuiltinParamDescriptor = BuiltinParamDescriptor {
43 name: "t",
44 ty: BuiltinParamType::NumericArray,
45 arity: BuiltinParamArity::Optional,
46 default: None,
47 description: "Time samples.",
48};
49const STEPINFO_PARAM_YFINAL: BuiltinParamDescriptor = BuiltinParamDescriptor {
50 name: "yfinal",
51 ty: BuiltinParamType::NumericScalar,
52 arity: BuiltinParamArity::Optional,
53 default: None,
54 description: "Steady-state value.",
55};
56const STEPINFO_INPUT_Y_T: [BuiltinParamDescriptor; 2] = [STEPINFO_PARAM_Y, STEPINFO_PARAM_T];
57const STEPINFO_INPUT_Y_T_FINAL: [BuiltinParamDescriptor; 3] =
58 [STEPINFO_PARAM_Y, STEPINFO_PARAM_T, STEPINFO_PARAM_YFINAL];
59const STEPINFO_SIGNATURES: [BuiltinSignatureDescriptor; 3] = [
60 BuiltinSignatureDescriptor {
61 label: "info = stepinfo(sys)",
62 inputs: &STEPINFO_INPUT_SYS,
63 outputs: &STEPINFO_OUTPUT,
64 },
65 BuiltinSignatureDescriptor {
66 label: "info = stepinfo(y, t)",
67 inputs: &STEPINFO_INPUT_Y_T,
68 outputs: &STEPINFO_OUTPUT,
69 },
70 BuiltinSignatureDescriptor {
71 label: "info = stepinfo(y, t, yfinal)",
72 inputs: &STEPINFO_INPUT_Y_T_FINAL,
73 outputs: &STEPINFO_OUTPUT,
74 },
75];
76const STEPINFO_ERRORS: [BuiltinErrorDescriptor; 5] = [
77 BuiltinErrorDescriptor {
78 code: "RM.STEPINFO.INVALID_ARGUMENT",
79 identifier: Some("RunMat:stepinfo:InvalidArgument"),
80 when: "Inputs do not match supported stepinfo invocation forms.",
81 message: "stepinfo: invalid argument",
82 },
83 BuiltinErrorDescriptor {
84 code: "RM.STEPINFO.INVALID_DATA",
85 identifier: Some("RunMat:stepinfo:InvalidData"),
86 when: "Response, time, or final-value data is malformed.",
87 message: "stepinfo: invalid response data",
88 },
89 BuiltinErrorDescriptor {
90 code: "RM.STEPINFO.INVALID_SYSTEM",
91 identifier: Some("RunMat:stepinfo:InvalidSystem"),
92 when: "System input is not a supported SISO tf object.",
93 message: "stepinfo: invalid system",
94 },
95 BuiltinErrorDescriptor {
96 code: "RM.STEPINFO.UNSUPPORTED_MODEL",
97 identifier: Some("RunMat:stepinfo:UnsupportedModel"),
98 when: "System model form is unsupported.",
99 message: "stepinfo: unsupported model",
100 },
101 BuiltinErrorDescriptor {
102 code: "RM.STEPINFO.INTERNAL",
103 identifier: Some("RunMat:stepinfo:Internal"),
104 when: "Response simulation or metric assembly failed.",
105 message: "stepinfo: internal error",
106 },
107];
108pub const STEPINFO_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
109 signatures: &STEPINFO_SIGNATURES,
110 output_mode: BuiltinOutputMode::Fixed,
111 completion_policy: BuiltinCompletionPolicy::Public,
112 errors: &STEPINFO_ERRORS,
113};
114
115#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::control::stepinfo")]
116pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
117 name: "stepinfo",
118 op_kind: GpuOpKind::Custom("control-stepinfo"),
119 supported_precisions: &[],
120 broadcast: BroadcastSemantics::None,
121 provider_hooks: &[],
122 constant_strategy: ConstantStrategy::InlineLiteral,
123 residency: ResidencyPolicy::GatherImmediately,
124 nan_mode: ReductionNaN::Include,
125 two_pass_threshold: None,
126 workgroup_size: None,
127 accepts_nan_mode: false,
128 notes: "stepinfo gathers sampled response data and computes scalar metrics on the host.",
129};
130
131#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::control::stepinfo")]
132pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
133 name: "stepinfo",
134 shape: ShapeRequirements::Any,
135 constant_strategy: ConstantStrategy::InlineLiteral,
136 elementwise: None,
137 reduction: None,
138 emits_nan: false,
139 notes: "stepinfo computes response metrics and is not fused.",
140};
141
142#[runtime_builtin(
143 name = "stepinfo",
144 category = "control",
145 summary = "Compute step-response metrics from SISO models or sampled responses.",
146 keywords = "stepinfo,step response,rise time,settling time,overshoot,control system",
147 type_resolver(stepinfo_type),
148 descriptor(crate::builtins::control::stepinfo::STEPINFO_DESCRIPTOR),
149 builtin_path = "crate::builtins::control::stepinfo"
150)]
151async fn stepinfo_builtin(input: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
152 let gathered_input = crate::dispatcher::gather_if_needed_async(&input).await?;
153 if matches!(gathered_input, Value::Object(_)) {
154 if !rest.is_empty() && !is_name_value_start(&rest[0]) {
155 return Err(stepinfo_error(
156 "stepinfo(sys) does not accept positional response data",
157 ));
158 }
159 let threshold = parse_options(&rest, 0)?;
160 let response = crate::call_builtin_async_with_outputs(
161 "step",
162 std::slice::from_ref(&gathered_input),
163 2,
164 )
165 .await?;
166 let (y, t) = output_list_two(response)?;
167 let y = numeric_vector(y, "y").await?;
168 let t = numeric_vector(t, "t").await?;
169 let gain = crate::call_builtin_async("dcgain", &[gathered_input]).await?;
170 let y_final = match output_complex_scalar_from_value(gain)? {
171 Some(value) => value,
172 None => *y
173 .last()
174 .ok_or_else(|| stepinfo_error("response vector cannot be empty"))?,
175 };
176 return metrics_to_struct(compute_metrics(&y, &t, y_final, threshold)?);
177 }
178
179 let y = numeric_vector(gathered_input, "y").await?;
180 let (t, y_final, option_start) = parse_sampled_response_tail(&y, &rest).await?;
181 let threshold = parse_options(&rest, option_start)?;
182 metrics_to_struct(compute_metrics(&y, &t, y_final, threshold)?)
183}
184
185async fn parse_sampled_response_tail(
186 y: &[f64],
187 rest: &[Value],
188) -> BuiltinResult<(Vec<f64>, f64, usize)> {
189 if rest.is_empty() || is_name_value_start(&rest[0]) {
190 let t = (0..y.len()).map(|idx| idx as f64).collect::<Vec<_>>();
191 let y_final = *y
192 .last()
193 .ok_or_else(|| stepinfo_error("response vector cannot be empty"))?;
194 return Ok((t, y_final, 0));
195 }
196 let t = numeric_vector(rest[0].clone(), "t").await?;
197 let mut option_start = 1;
198 let y_final = if rest.len() > 1 && !is_name_value_start(&rest[1]) {
199 let gathered = crate::dispatcher::gather_if_needed_async(&rest[1]).await?;
200 option_start = 2;
201 scalar_f64(&gathered, "yfinal", BUILTIN_NAME)?
202 } else {
203 *y.last()
204 .ok_or_else(|| stepinfo_error("response vector cannot be empty"))?
205 };
206 Ok((t, y_final, option_start))
207}
208
209fn parse_options(rest: &[Value], start: usize) -> BuiltinResult<f64> {
210 if rest.len() == start {
211 return Ok(DEFAULT_SETTLING_THRESHOLD);
212 }
213 if !(rest.len() - start).is_multiple_of(2) {
214 return Err(stepinfo_error("name-value options must come in pairs"));
215 }
216 let mut threshold = DEFAULT_SETTLING_THRESHOLD;
217 let mut idx = start;
218 while idx < rest.len() {
219 let name = scalar_text(&rest[idx], "option name", BUILTIN_NAME)?;
220 match name.trim().to_ascii_lowercase().as_str() {
221 "settlingtimethreshold" => {
222 let value = scalar_f64(&rest[idx + 1], "SettlingTimeThreshold", BUILTIN_NAME)?;
223 if !value.is_finite() || value <= 0.0 || value >= 1.0 {
224 return Err(stepinfo_error(
225 "SettlingTimeThreshold must be a finite scalar between 0 and 1",
226 ));
227 }
228 threshold = value;
229 }
230 other => {
231 return Err(stepinfo_error(format!(
232 "unsupported stepinfo option '{other}'"
233 )))
234 }
235 }
236 idx += 2;
237 }
238 Ok(threshold)
239}
240
241fn is_name_value_start(value: &Value) -> bool {
242 matches!(
243 value,
244 Value::String(_) | Value::StringArray(_) | Value::CharArray(_)
245 )
246}
247
248async fn numeric_vector(value: Value, label: &str) -> BuiltinResult<Vec<f64>> {
249 let gathered = crate::dispatcher::gather_if_needed_async(&value).await?;
250 match gathered {
251 Value::Tensor(tensor) => {
252 ensure_vector_shape(&tensor, label)?;
253 finite_vector(tensor.data, label)
254 }
255 Value::Num(n) => finite_vector(vec![n], label),
256 Value::Int(i) => finite_vector(vec![i.to_f64()], label),
257 Value::Bool(b) => finite_vector(vec![if b { 1.0 } else { 0.0 }], label),
258 Value::LogicalArray(logical) => {
259 let data = logical
260 .data
261 .iter()
262 .map(|bit| if *bit == 0 { 0.0 } else { 1.0 })
263 .collect::<Vec<_>>();
264 finite_vector(data, label)
265 }
266 other => Err(stepinfo_error(format!(
267 "{label} must be a real numeric vector, got {other:?}"
268 ))),
269 }
270}
271
272fn ensure_vector_shape(tensor: &Tensor, label: &str) -> BuiltinResult<()> {
273 let non_unit = tensor.shape.iter().copied().filter(|&dim| dim > 1).count();
274 if non_unit <= 1 {
275 Ok(())
276 } else {
277 Err(stepinfo_error(format!("{label} must be a vector")))
278 }
279}
280
281fn finite_vector(values: Vec<f64>, label: &str) -> BuiltinResult<Vec<f64>> {
282 if values.is_empty() {
283 return Err(stepinfo_error(format!("{label} vector cannot be empty")));
284 }
285 if values.iter().any(|value| !value.is_finite()) {
286 return Err(stepinfo_error(format!("{label} values must be finite")));
287 }
288 Ok(values)
289}
290
291fn output_list_two(value: Value) -> BuiltinResult<(Value, Value)> {
292 let Value::OutputList(outputs) = value else {
293 return Err(control_error(
294 BUILTIN_NAME,
295 "RunMat:stepinfo:Internal",
296 "stepinfo: step did not return an output list",
297 ));
298 };
299 if outputs.len() < 2 {
300 return Err(control_error(
301 BUILTIN_NAME,
302 "RunMat:stepinfo:Internal",
303 "stepinfo: step returned too few outputs",
304 ));
305 }
306 Ok((outputs[0].clone(), outputs[1].clone()))
307}
308
309fn output_complex_scalar_from_value(value: Value) -> BuiltinResult<Option<f64>> {
310 match value {
311 Value::Num(n) if n.is_finite() => Ok(Some(n)),
312 Value::Int(i) => Ok(Some(i.to_f64())),
313 Value::Complex(re, im) if im.abs() <= EPS && re.is_finite() => Ok(Some(re)),
314 Value::Complex(_, _) => Ok(None),
315 _ => Ok(None),
316 }
317}
318
319#[derive(Clone, Debug)]
320struct StepMetrics {
321 rise_time: f64,
322 transient_time: f64,
323 settling_time: f64,
324 settling_min: f64,
325 settling_max: f64,
326 overshoot: f64,
327 undershoot: f64,
328 peak: f64,
329 peak_time: f64,
330 steady_state_value: f64,
331}
332
333fn compute_metrics(
334 y: &[f64],
335 t: &[f64],
336 y_final: f64,
337 settling_threshold: f64,
338) -> BuiltinResult<StepMetrics> {
339 if y.len() != t.len() {
340 return Err(stepinfo_error(
341 "y and t must have the same number of samples",
342 ));
343 }
344 if y.is_empty() {
345 return Err(stepinfo_error("response vector cannot be empty"));
346 }
347 if !y_final.is_finite() {
348 return Err(stepinfo_error("yfinal must be finite"));
349 }
350 validate_time(t)?;
351
352 let y0 = y[0];
353 let amplitude = y_final - y0;
354 let amplitude_abs = amplitude.abs();
355 let scale = amplitude_abs.max(y_final.abs()).max(1.0);
356 let direction = if amplitude >= 0.0 { 1.0 } else { -1.0 };
357
358 let (rise_time, rise_idx) = if amplitude_abs <= EPS {
359 (f64::NAN, 0)
360 } else {
361 let t10 = crossing_time(y, t, y0 + 0.1 * amplitude, direction);
362 let (t90, idx90) = crossing_time_with_index(y, t, y0 + 0.9 * amplitude, direction);
363 match (t10, t90) {
364 (Some(start), Some(end)) => (end - start, idx90.unwrap_or(0)),
365 _ => (f64::NAN, 0),
366 }
367 };
368
369 let settling_band = settling_threshold * scale;
370 let mut last_outside = None;
371 for (idx, value) in y.iter().enumerate() {
372 if (*value - y_final).abs() > settling_band {
373 last_outside = Some(idx);
374 }
375 }
376 let settling_time = match last_outside {
377 None => t[0],
378 Some(idx) if idx + 1 < t.len() => t[idx + 1],
379 Some(_) => f64::NAN,
380 };
381
382 let settle_window = &y[rise_idx.min(y.len() - 1)..];
383 let settling_min = settle_window.iter().copied().fold(f64::INFINITY, f64::min);
384 let settling_max = settle_window
385 .iter()
386 .copied()
387 .fold(f64::NEG_INFINITY, f64::max);
388
389 let max_y = y.iter().copied().fold(f64::NEG_INFINITY, f64::max);
390 let min_y = y.iter().copied().fold(f64::INFINITY, f64::min);
391 let overshoot = if amplitude_abs <= EPS {
392 0.0
393 } else if amplitude >= 0.0 {
394 ((max_y - y_final) / amplitude_abs * 100.0).max(0.0)
395 } else {
396 ((y_final - min_y) / amplitude_abs * 100.0).max(0.0)
397 };
398 let undershoot = if amplitude_abs <= EPS {
399 0.0
400 } else if amplitude >= 0.0 {
401 ((y0 - min_y) / amplitude_abs * 100.0).max(0.0)
402 } else {
403 ((max_y - y0) / amplitude_abs * 100.0).max(0.0)
404 };
405
406 let (peak_idx, peak) = y
407 .iter()
408 .enumerate()
409 .map(|(idx, value)| (idx, value.abs()))
410 .max_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal))
411 .unwrap_or((0, y[0].abs()));
412
413 Ok(StepMetrics {
414 rise_time,
415 transient_time: settling_time,
416 settling_time,
417 settling_min,
418 settling_max,
419 overshoot,
420 undershoot,
421 peak,
422 peak_time: t[peak_idx],
423 steady_state_value: y_final,
424 })
425}
426
427fn validate_time(t: &[f64]) -> BuiltinResult<()> {
428 if t.iter().any(|value| !value.is_finite()) {
429 return Err(stepinfo_error("t values must be finite"));
430 }
431 if t.windows(2).any(|pair| pair[1] < pair[0]) {
432 return Err(stepinfo_error("t must be nondecreasing"));
433 }
434 Ok(())
435}
436
437fn crossing_time(y: &[f64], t: &[f64], threshold: f64, direction: f64) -> Option<f64> {
438 crossing_time_with_index(y, t, threshold, direction).0
439}
440
441fn crossing_time_with_index(
442 y: &[f64],
443 t: &[f64],
444 threshold: f64,
445 direction: f64,
446) -> (Option<f64>, Option<usize>) {
447 for idx in 0..y.len() {
448 let current = direction * (y[idx] - threshold);
449 if current >= -EPS {
450 if idx == 0 {
451 return (Some(t[0]), Some(0));
452 }
453 let prev = direction * (y[idx - 1] - threshold);
454 let denom = current - prev;
455 if denom.abs() <= EPS {
456 return (Some(t[idx]), Some(idx));
457 }
458 let alpha = (-prev / denom).clamp(0.0, 1.0);
459 return (Some(t[idx - 1] + alpha * (t[idx] - t[idx - 1])), Some(idx));
460 }
461 }
462 (None, None)
463}
464
465fn metrics_to_struct(metrics: StepMetrics) -> BuiltinResult<Value> {
466 let mut out = StructValue::new();
467 out.insert("RiseTime", Value::Num(metrics.rise_time));
468 out.insert("TransientTime", Value::Num(metrics.transient_time));
469 out.insert("SettlingTime", Value::Num(metrics.settling_time));
470 out.insert("SettlingMin", Value::Num(metrics.settling_min));
471 out.insert("SettlingMax", Value::Num(metrics.settling_max));
472 out.insert("Overshoot", Value::Num(metrics.overshoot));
473 out.insert("Undershoot", Value::Num(metrics.undershoot));
474 out.insert("Peak", Value::Num(metrics.peak));
475 out.insert("PeakTime", Value::Num(metrics.peak_time));
476 out.insert("SteadyStateValue", Value::Num(metrics.steady_state_value));
477 Ok(Value::Struct(out))
478}
479
480fn stepinfo_error(message: impl Into<String>) -> RuntimeError {
481 control_error(BUILTIN_NAME, "RunMat:stepinfo:InvalidData", message)
482}
483
484#[cfg(test)]
485mod tests {
486 use super::*;
487 use futures::executor::block_on;
488
489 #[test]
490 fn sampled_response_reports_basic_metrics() {
491 let y = Value::Tensor(Tensor::new(vec![0.0, 0.5, 0.9, 1.0], vec![1, 4]).unwrap());
492 let t = Value::Tensor(Tensor::new(vec![0.0, 1.0, 2.0, 3.0], vec![1, 4]).unwrap());
493 let Value::Struct(info) =
494 block_on(stepinfo_builtin(y, vec![t, Value::Num(1.0)])).expect("stepinfo")
495 else {
496 panic!("expected struct");
497 };
498 assert!(
499 matches!(info.fields.get("RiseTime"), Some(Value::Num(v)) if (*v - 1.8).abs() < 1.0e-12)
500 );
501 assert!(matches!(info.fields.get("Overshoot"), Some(Value::Num(v)) if v.abs() < 1.0e-12));
502 }
503
504 #[test]
505 fn system_response_form_runs_step_and_dcgain() {
506 let sys = block_on(crate::call_builtin_async(
507 "tf",
508 &[
509 Value::Num(1.0),
510 Value::Tensor(Tensor::new(vec![1.0, 1.0], vec![1, 2]).unwrap()),
511 ],
512 ))
513 .expect("tf");
514 let Value::Struct(info) = block_on(stepinfo_builtin(sys, Vec::new())).expect("stepinfo")
515 else {
516 panic!("expected struct");
517 };
518 assert!(
519 matches!(info.fields.get("SteadyStateValue"), Some(Value::Num(v)) if (*v - 1.0).abs() < 1.0e-12)
520 );
521 }
522}