use nalgebra::DMatrix;
use num_complex::Complex64;
use runmat_builtins::{ComplexTensor, ObjectInstance, Tensor, Value};
use runmat_macros::runtime_builtin;
use crate::builtins::common::spec::{
BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
ReductionNaN, ResidencyPolicy, ShapeRequirements,
};
use crate::builtins::control::type_resolvers::step_type;
use crate::{build_runtime_error, BuiltinResult, RuntimeError};
const BUILTIN_NAME: &str = "step";
const EPS: f64 = 1.0e-12;
const DEFAULT_SAMPLES: usize = 101;
const MAX_DISCRETE_SAMPLES: usize = 1_000_000;
#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::control::step")]
pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
name: "step",
op_kind: GpuOpKind::Custom("control-step-response"),
supported_precisions: &[],
broadcast: BroadcastSemantics::None,
provider_hooks: &[],
constant_strategy: ConstantStrategy::InlineLiteral,
residency: ResidencyPolicy::GatherImmediately,
nan_mode: ReductionNaN::Include,
two_pass_threshold: None,
workgroup_size: None,
accepts_nan_mode: false,
notes: "Step-response simulation runs on the host from transfer-function metadata.",
};
#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::control::step")]
pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
name: "step",
shape: ShapeRequirements::Any,
constant_strategy: ConstantStrategy::InlineLiteral,
elementwise: None,
reduction: None,
emits_nan: false,
notes: "step simulates a dynamic system and terminates numeric fusion chains.",
};
fn step_error(message: impl Into<String>) -> RuntimeError {
build_runtime_error(message)
.with_builtin(BUILTIN_NAME)
.build()
}
#[runtime_builtin(
name = "step",
category = "control",
summary = "Compute or plot the step response of a SISO transfer-function model.",
keywords = "step,response,control system,transfer function,tf",
sink = true,
suppress_auto_output = true,
type_resolver(step_type),
builtin_path = "crate::builtins::control::step"
)]
async fn step_builtin(sys: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
if rest.len() > 1 {
return Err(step_error(
"step: expected step(sys), step(sys, tFinal), or step(sys, t)",
));
}
let sys = crate::gather_if_needed_async(&sys).await?;
let mut rest_host = Vec::with_capacity(rest.len());
for arg in &rest {
rest_host.push(crate::gather_if_needed_async(arg).await?);
}
let model = TransferFunction::from_value(sys)?;
let time = TimeSpec::parse(rest_host.first(), model.sample_time)?;
let eval = evaluate_step(&model, time)?;
if crate::output_context::requested_output_count() == Some(0)
&& crate::output_count::current_output_count().is_none()
{
plot_response(&eval).await?;
return Ok(Value::OutputList(Vec::new()));
}
if let Some(out_count) = crate::output_count::current_output_count() {
if out_count == 0 {
plot_response(&eval).await?;
return Ok(Value::OutputList(Vec::new()));
}
if out_count == 1 {
return Ok(Value::OutputList(vec![eval.y_value()?]));
}
return Ok(crate::output_count::output_list_with_padding(
out_count,
eval.outputs()?,
));
}
eval.y_value()
}
#[derive(Clone, Debug)]
struct TransferFunction {
numerator: Vec<f64>,
denominator: Vec<f64>,
sample_time: f64,
}
impl TransferFunction {
fn from_value(value: Value) -> BuiltinResult<Self> {
let Value::Object(object) = value else {
return Err(step_error("step: expected a tf object"));
};
if !object.is_class("tf") {
return Err(step_error(format!(
"step: expected a tf object, got {}",
object.class_name
)));
}
let numerator = property_coefficients(&object, "Numerator")?;
let denominator = property_coefficients(&object, "Denominator")?;
let sample_time = property_scalar(&object, "Ts")?;
if !sample_time.is_finite() || sample_time < 0.0 {
return Err(step_error(
"step: tf sample time must be finite and non-negative",
));
}
let numerator = trim_leading_zeros(numerator);
let denominator = trim_leading_zeros(denominator);
if denominator.is_empty() {
return Err(step_error("step: denominator coefficients cannot be empty"));
}
if denominator[0].abs() <= EPS {
return Err(step_error(
"step: leading denominator coefficient must be non-zero",
));
}
if numerator.len().saturating_sub(1) > denominator.len().saturating_sub(1) {
return Err(step_error(
"step: improper transfer functions are not supported yet",
));
}
Ok(Self {
numerator,
denominator,
sample_time,
})
}
fn normalized(&self) -> (Vec<f64>, Vec<f64>) {
let leading = self.denominator[0];
let den = self
.denominator
.iter()
.map(|value| value / leading)
.collect::<Vec<_>>();
let num = self
.numerator
.iter()
.map(|value| value / leading)
.collect::<Vec<_>>();
(num, den)
}
}
fn property_coefficients(object: &ObjectInstance, name: &str) -> BuiltinResult<Vec<f64>> {
let value = object
.properties
.get(name)
.ok_or_else(|| step_error(format!("step: tf object is missing {name}")))?;
match value {
Value::Tensor(tensor) => Ok(tensor.data.clone()),
Value::ComplexTensor(tensor) => real_complex_coefficients(tensor, name),
Value::Num(n) => Ok(vec![*n]),
Value::Int(i) => Ok(vec![i.to_f64()]),
other => Err(step_error(format!(
"step: tf {name} coefficients must be numeric, got {other:?}"
))),
}
}
fn real_complex_coefficients(tensor: &ComplexTensor, name: &str) -> BuiltinResult<Vec<f64>> {
let mut out = Vec::with_capacity(tensor.data.len());
for &(re, im) in &tensor.data {
if im.abs() > EPS {
return Err(step_error(format!(
"step: complex tf {name} coefficients are not supported yet"
)));
}
out.push(re);
}
Ok(out)
}
fn property_scalar(object: &ObjectInstance, name: &str) -> BuiltinResult<f64> {
let value = object
.properties
.get(name)
.ok_or_else(|| step_error(format!("step: tf object is missing {name}")))?;
match value {
Value::Num(n) => Ok(*n),
Value::Int(i) => Ok(i.to_f64()),
other => Err(step_error(format!(
"step: tf {name} property must be a scalar, got {other:?}"
))),
}
}
#[derive(Clone, Debug)]
enum TimeSpec {
Auto,
FinalTime(f64),
Vector(Vec<f64>),
}
impl TimeSpec {
fn parse(value: Option<&Value>, sample_time: f64) -> BuiltinResult<Self> {
let Some(value) = value else {
return Ok(Self::Auto);
};
match value {
Value::Num(n) => Self::final_time(*n),
Value::Int(i) => Self::final_time(i.to_f64()),
Value::Tensor(tensor) => {
ensure_time_vector_shape(&tensor.shape)?;
if tensor.data.len() == 1 {
return Self::final_time(tensor.data[0]);
}
Self::vector(tensor.data.clone(), sample_time)
}
other => Err(step_error(format!(
"step: time input must be a scalar final time or numeric vector, got {other:?}"
))),
}
}
fn final_time(value: f64) -> BuiltinResult<Self> {
if !value.is_finite() || value <= 0.0 {
return Err(step_error(
"step: final time must be a positive finite scalar",
));
}
Ok(Self::FinalTime(value))
}
fn vector(values: Vec<f64>, sample_time: f64) -> BuiltinResult<Self> {
validate_time_vector(&values)?;
if sample_time > 0.0 {
for &t in &values {
let k = (t / sample_time).round();
if (t - k * sample_time).abs() > 1.0e-8 * sample_time.max(1.0) {
return Err(step_error(
"step: discrete-time sample vector must align with the model sample time",
));
}
}
}
Ok(Self::Vector(values))
}
}
fn ensure_time_vector_shape(shape: &[usize]) -> BuiltinResult<()> {
let non_unit = shape.iter().copied().filter(|&dim| dim > 1).count();
if non_unit <= 1 {
Ok(())
} else {
Err(step_error("step: time input must be a vector"))
}
}
fn validate_time_vector(values: &[f64]) -> BuiltinResult<()> {
if values.is_empty() {
return Err(step_error("step: time vector must not be empty"));
}
let mut previous = None;
for &value in values {
if !value.is_finite() || value < 0.0 {
return Err(step_error(
"step: time vector values must be finite and non-negative",
));
}
if let Some(prev) = previous {
if value < prev {
return Err(step_error("step: time vector must be nondecreasing"));
}
}
previous = Some(value);
}
Ok(())
}
#[derive(Clone, Debug)]
struct StepEval {
y: Vec<f64>,
t: Vec<f64>,
}
impl StepEval {
fn y_value(&self) -> BuiltinResult<Value> {
column_tensor(self.y.clone())
}
fn t_value(&self) -> BuiltinResult<Value> {
column_tensor(self.t.clone())
}
fn outputs(&self) -> BuiltinResult<Vec<Value>> {
Ok(vec![self.y_value()?, self.t_value()?])
}
}
fn evaluate_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
if model.sample_time > 0.0 {
evaluate_discrete_step(model, time)
} else {
evaluate_continuous_step(model, time)
}
}
fn evaluate_continuous_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
let t = continuous_time_vector(model, time)?;
let (num, den) = model.normalized();
let response = continuous_response(&num, &den, &t)?;
Ok(StepEval { y: response, t })
}
fn continuous_time_vector(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<Vec<f64>> {
match time {
TimeSpec::Auto => Ok(linspace(0.0, automatic_final_time(model), DEFAULT_SAMPLES)),
TimeSpec::FinalTime(final_time) => Ok(linspace(0.0, final_time, DEFAULT_SAMPLES)),
TimeSpec::Vector(values) => Ok(values),
}
}
fn continuous_response(num: &[f64], den: &[f64], t: &[f64]) -> BuiltinResult<Vec<f64>> {
let order = den.len() - 1;
if order == 0 {
let gain = num.last().copied().unwrap_or(0.0) / den[0];
return Ok(vec![gain; t.len()]);
}
let mut padded_num = vec![0.0; order + 1 - num.len()];
padded_num.extend_from_slice(num);
let direct = padded_num[0];
let a = &den[1..];
let mut c = Vec::with_capacity(order);
for state_idx in 0..order {
let coeff_idx = order - state_idx;
c.push(padded_num[coeff_idx] - direct * den[coeff_idx]);
}
let mut state = vec![0.0; order];
let mut current_t = 0.0;
let mut response = Vec::with_capacity(t.len());
for &target_t in t {
if target_t > current_t {
integrate_to(&mut state, a, current_t, target_t);
current_t = target_t;
}
response.push(dot(&c, &state) + direct);
}
Ok(response)
}
fn integrate_to(state: &mut [f64], a: &[f64], start: f64, end: f64) {
let duration = end - start;
if duration <= 0.0 {
return;
}
let steps = ((duration / 0.01).ceil() as usize).clamp(1, 10_000);
let h = duration / steps as f64;
for _ in 0..steps {
rk4_step(state, a, h);
}
}
fn rk4_step(state: &mut [f64], a: &[f64], h: f64) {
let k1 = derivative(state, a);
let s2 = add_scaled(state, &k1, h * 0.5);
let k2 = derivative(&s2, a);
let s3 = add_scaled(state, &k2, h * 0.5);
let k3 = derivative(&s3, a);
let s4 = add_scaled(state, &k3, h);
let k4 = derivative(&s4, a);
for idx in 0..state.len() {
state[idx] += h * (k1[idx] + 2.0 * k2[idx] + 2.0 * k3[idx] + k4[idx]) / 6.0;
}
}
fn derivative(state: &[f64], a: &[f64]) -> Vec<f64> {
let order = state.len();
let mut dx = vec![0.0; order];
if order > 1 {
dx[..(order - 1)].copy_from_slice(&state[1..order]);
}
let mut last = 1.0;
for state_idx in 0..order {
let coeff = a[order - 1 - state_idx];
last -= coeff * state[state_idx];
}
dx[order - 1] = last;
dx
}
fn add_scaled(state: &[f64], delta: &[f64], scale: f64) -> Vec<f64> {
state
.iter()
.zip(delta)
.map(|(value, delta)| value + scale * delta)
.collect()
}
fn evaluate_discrete_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
let t = discrete_time_vector(model.sample_time, time)?;
if t.len() > MAX_DISCRETE_SAMPLES {
return Err(step_error(format!(
"step: discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"
)));
}
let sample_indices = t
.iter()
.map(|&value| checked_discrete_sample_index(model.sample_time, value))
.collect::<BuiltinResult<Vec<_>>>()?;
let max_k = sample_indices.iter().copied().max().unwrap_or(0);
let count = max_k
.checked_add(1)
.ok_or_else(|| step_error("step: discrete sample index exceeds platform limits"))?;
let (num, den) = model.normalized();
let all_y = discrete_response(&num, &den, count)?;
let y = sample_indices
.into_iter()
.map(|idx| {
all_y
.get(idx)
.copied()
.ok_or_else(|| step_error("step: discrete sample index exceeds response length"))
})
.collect::<BuiltinResult<Vec<_>>>()?;
Ok(StepEval { y, t })
}
fn discrete_time_vector(sample_time: f64, time: TimeSpec) -> BuiltinResult<Vec<f64>> {
match time {
TimeSpec::Auto => Ok((0..DEFAULT_SAMPLES)
.map(|idx| idx as f64 * sample_time)
.collect()),
TimeSpec::FinalTime(final_time) => {
let steps = checked_discrete_sample_steps(sample_time, final_time)?;
Ok((0..=steps).map(|idx| idx as f64 * sample_time).collect())
}
TimeSpec::Vector(values) => Ok(values),
}
}
fn checked_discrete_sample_steps(sample_time: f64, final_time: f64) -> BuiltinResult<usize> {
let steps = (final_time / sample_time).floor();
if !steps.is_finite() || steps < 0.0 || steps > usize::MAX as f64 {
return Err(step_error(
"step: discrete sample count exceeds platform limits",
));
}
if steps >= MAX_DISCRETE_SAMPLES as f64 {
return Err(step_error(format!(
"step: discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"
)));
}
Ok(steps as usize)
}
fn checked_discrete_sample_index(sample_time: f64, time: f64) -> BuiltinResult<usize> {
let index = (time / sample_time).round();
if !index.is_finite() || index < 0.0 || index > usize::MAX as f64 {
return Err(step_error(
"step: discrete sample index exceeds platform limits",
));
}
if index >= MAX_DISCRETE_SAMPLES as f64 {
return Err(step_error(format!(
"step: discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"
)));
}
Ok(index as usize)
}
fn discrete_response(num: &[f64], den: &[f64], count: usize) -> BuiltinResult<Vec<f64>> {
let order = den.len() - 1;
let mut padded_num = vec![0.0; order + 1 - num.len()];
padded_num.extend_from_slice(num);
let mut y = vec![0.0; count];
for k in 0..count {
let mut value = 0.0;
for (idx, &coeff) in padded_num.iter().enumerate() {
if k >= idx {
value += coeff;
}
}
for idx in 1..den.len() {
if k >= idx {
value -= den[idx] * y[k - idx];
}
}
y[k] = value;
}
Ok(y)
}
async fn plot_response(eval: &StepEval) -> BuiltinResult<()> {
let t = eval.t_value()?;
let y = eval.y_value()?;
if let Err(err) = crate::call_builtin_async("plot", &[t, y]).await {
if is_nonfatal_plot_setup_error(&err) {
return Ok(());
}
return Err(err);
}
let _ = crate::call_builtin_async("title", &[Value::from("Step Response")]).await;
let _ = crate::call_builtin_async("xlabel", &[Value::from("Time")]).await;
let _ = crate::call_builtin_async("ylabel", &[Value::from("Amplitude")]).await;
Ok(())
}
fn is_nonfatal_plot_setup_error(err: &RuntimeError) -> bool {
let lower = err.to_string().to_ascii_lowercase();
lower.contains("plotting is unavailable")
|| lower.contains("non-main thread")
|| lower.contains("interactive plotting failed")
|| lower.contains("eventloop can't be recreated")
}
fn automatic_final_time(model: &TransferFunction) -> f64 {
let (_, den) = model.normalized();
let poles = polynomial_roots(&den).unwrap_or_default();
let slowest_decay = poles
.iter()
.filter_map(|pole| if pole.re < -EPS { Some(-pole.re) } else { None })
.fold(f64::INFINITY, f64::min);
if slowest_decay.is_finite() && slowest_decay > EPS {
(5.0 / slowest_decay).clamp(1.0, 100.0)
} else {
10.0
}
}
fn polynomial_roots(coeffs: &[f64]) -> BuiltinResult<Vec<Complex64>> {
let trimmed = trim_leading_zeros(coeffs.to_vec());
if trimmed.len() <= 1 {
return Ok(Vec::new());
}
if trimmed.len() == 2 {
return Ok(vec![Complex64::new(-trimmed[1] / trimmed[0], 0.0)]);
}
let degree = trimmed.len() - 1;
let leading = trimmed[0];
let mut companion = DMatrix::<Complex64>::zeros(degree, degree);
for row in 1..degree {
companion[(row, row - 1)] = Complex64::new(1.0, 0.0);
}
for (idx, coeff) in trimmed.iter().enumerate().skip(1) {
companion[(0, idx - 1)] = Complex64::new(-coeff / leading, 0.0);
}
let eigenvalues = companion
.eigenvalues()
.ok_or_else(|| step_error("step: failed to compute transfer-function poles"))?;
Ok(eigenvalues.iter().copied().collect())
}
fn trim_leading_zeros(coeffs: Vec<f64>) -> Vec<f64> {
let first_nonzero = coeffs
.iter()
.position(|value| value.abs() > EPS)
.unwrap_or(coeffs.len());
coeffs[first_nonzero..].to_vec()
}
fn linspace(start: f64, end: f64, count: usize) -> Vec<f64> {
if count <= 1 {
return vec![end];
}
let step = (end - start) / (count - 1) as f64;
(0..count).map(|idx| start + idx as f64 * step).collect()
}
fn dot(lhs: &[f64], rhs: &[f64]) -> f64 {
lhs.iter().zip(rhs).map(|(a, b)| a * b).sum()
}
fn column_tensor(data: Vec<f64>) -> BuiltinResult<Value> {
let rows = data.len();
let tensor =
Tensor::new(data, vec![rows, 1]).map_err(|err| step_error(format!("step: {err}")))?;
Ok(Value::Tensor(tensor))
}
#[cfg(test)]
mod tests {
use super::*;
use futures::executor::block_on;
use runmat_builtins::{CharArray, ObjectInstance};
fn tf_object(num: Vec<f64>, den: Vec<f64>, sample_time: f64) -> Value {
let mut object = ObjectInstance::new("tf".to_string());
object.properties.insert(
"Numerator".to_string(),
Value::Tensor(Tensor::new(num.clone(), vec![1, num.len()]).unwrap()),
);
object.properties.insert(
"Denominator".to_string(),
Value::Tensor(Tensor::new(den.clone(), vec![1, den.len()]).unwrap()),
);
object.properties.insert(
"Variable".to_string(),
Value::CharArray(CharArray::new_row(if sample_time > 0.0 {
"z"
} else {
"s"
})),
);
object
.properties
.insert("Ts".to_string(), Value::Num(sample_time));
Value::Object(object)
}
fn run_step(sys: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
block_on(step_builtin(sys, rest))
}
fn tensor_data(value: Value) -> Vec<f64> {
match value {
Value::Tensor(tensor) => tensor.data,
other => panic!("expected tensor, got {other:?}"),
}
}
#[test]
fn first_order_continuous_response_matches_closed_form_for_explicit_time() {
let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
let time = Value::Tensor(Tensor::new(vec![0.0, 0.5, 1.0, 2.0], vec![1, 4]).unwrap());
let y = tensor_data(run_step(sys, vec![time]).expect("step"));
for (actual, t) in y.iter().zip([0.0_f64, 0.5, 1.0, 2.0]) {
let expected = 1.0 - (-t).exp();
assert!(
(actual - expected).abs() < 1.0e-5,
"t={t} actual={actual} expected={expected}"
);
}
}
#[test]
fn multi_output_returns_y_then_time() {
let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
let _guard = crate::output_count::push_output_count(Some(2));
let time = Value::Tensor(Tensor::new(vec![0.0, 1.0], vec![2, 1]).unwrap());
let result = run_step(sys, vec![time]).expect("step");
let Value::OutputList(outputs) = result else {
panic!("expected output list");
};
assert_eq!(outputs.len(), 2);
assert_eq!(tensor_data(outputs[1].clone()), vec![0.0, 1.0]);
let y = tensor_data(outputs[0].clone());
assert!((y[1] - (1.0 - (-1.0_f64).exp())).abs() < 1.0e-5);
}
#[test]
fn single_requested_output_returns_only_response() {
let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
let _guard = crate::output_count::push_output_count(Some(1));
let time = Value::Tensor(Tensor::new(vec![0.0, 1.0], vec![1, 2]).unwrap());
let result = run_step(sys, vec![time]).expect("step");
let Value::OutputList(outputs) = result else {
panic!("expected output list");
};
assert_eq!(outputs.len(), 1);
let y = tensor_data(outputs[0].clone());
assert_eq!(y.len(), 2);
assert!((y[1] - (1.0 - (-1.0_f64).exp())).abs() < 1.0e-5);
}
#[test]
fn scalar_final_time_generates_column_time_vector_ending_at_final_time() {
let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
let _guard = crate::output_count::push_output_count(Some(2));
let result = run_step(sys, vec![Value::Num(5.0)]).expect("step");
let Value::OutputList(outputs) = result else {
panic!("expected output list");
};
let t = tensor_data(outputs[1].clone());
assert_eq!(t.len(), DEFAULT_SAMPLES);
assert_eq!(t[0], 0.0);
assert!((t[t.len() - 1] - 5.0).abs() < 1.0e-12);
}
#[test]
fn discrete_response_uses_sample_time_grid() {
let sys = tf_object(vec![1.0], vec![1.0, -0.5], 0.1);
let time = Value::Tensor(Tensor::new(vec![0.0, 0.1, 0.2], vec![1, 3]).unwrap());
let y = tensor_data(run_step(sys, vec![time]).expect("step"));
assert_eq!(y, vec![0.0, 1.0, 1.5]);
}
#[test]
fn discrete_final_time_rejects_excessive_sample_count() {
let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0e-6);
let err = run_step(sys, vec![Value::Num(2.0)]).expect_err("should fail");
assert!(err.message().contains("more than 1000000 samples"));
}
#[test]
fn discrete_time_vector_rejects_excessive_sample_index() {
let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0);
let time =
Value::Tensor(Tensor::new(vec![0.0, MAX_DISCRETE_SAMPLES as f64], vec![1, 2]).unwrap());
let err = run_step(sys, vec![time]).expect_err("should fail");
assert!(err.message().contains("more than 1000000 samples"));
}
#[test]
fn rejects_non_tf_input() {
let err = run_step(Value::Num(1.0), Vec::new()).expect_err("expected error");
assert!(err.message().contains("expected a tf object"));
}
}