1use nalgebra::DMatrix;
4use num_complex::Complex64;
5use runmat_builtins::{
6 BuiltinCompletionPolicy, BuiltinDescriptor, BuiltinErrorDescriptor, BuiltinExtensionDescriptor,
7 BuiltinExtensionMode, BuiltinIntegerBackendRule, BuiltinIntegerCapabilityDescriptor,
8 BuiltinIntegerComputationDomain, BuiltinIntegerInputAvailability,
9 BuiltinIntegerInputCapability, BuiltinIntegerOutputClassRule, BuiltinIntegerOverflowRule,
10 BuiltinIntegerOverloadKind, BuiltinIntegerScalarDoubleRule, BuiltinOutputMode,
11 BuiltinParamArity, BuiltinParamDescriptor, BuiltinParamType, BuiltinSignatureDescriptor,
12};
13use runmat_macros::runtime_builtin;
14use runmat_value::{ComplexTensor, ObjectInstance, Tensor, Value};
15
16use crate::builtins::common::spec::{
17 BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
18 ReductionNaN, ResidencyPolicy, ShapeRequirements,
19};
20use crate::builtins::common::tensor;
21use crate::builtins::control::type_resolvers::step_type;
22use crate::{build_runtime_error, BuiltinResult, RuntimeError};
23
24const BUILTIN_NAME: &str = "step";
25const EPS: f64 = 1.0e-12;
26const DEFAULT_SAMPLES: usize = 101;
27const MAX_DISCRETE_SAMPLES: usize = 1_000_000;
28
29const STEP_OUTPUT_Y: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
30 name: "y",
31 ty: BuiltinParamType::NumericArray,
32 arity: BuiltinParamArity::Required,
33 default: None,
34 description: "Step response samples (column vector).",
35}];
36const STEP_OUTPUT_Y_T: [BuiltinParamDescriptor; 2] = [
37 BuiltinParamDescriptor {
38 name: "y",
39 ty: BuiltinParamType::NumericArray,
40 arity: BuiltinParamArity::Required,
41 default: None,
42 description: "Step response samples (column vector).",
43 },
44 BuiltinParamDescriptor {
45 name: "t",
46 ty: BuiltinParamType::NumericArray,
47 arity: BuiltinParamArity::Required,
48 default: None,
49 description: "Time samples (column vector).",
50 },
51];
52const STEP_INPUTS_SYS: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
53 name: "sys",
54 ty: BuiltinParamType::Any,
55 arity: BuiltinParamArity::Required,
56 default: None,
57 description: "SISO tf model.",
58}];
59const STEP_INPUTS_SYS_TIME: [BuiltinParamDescriptor; 2] = [
60 BuiltinParamDescriptor {
61 name: "sys",
62 ty: BuiltinParamType::Any,
63 arity: BuiltinParamArity::Required,
64 default: None,
65 description: "SISO tf model.",
66 },
67 BuiltinParamDescriptor {
68 name: "time",
69 ty: BuiltinParamType::Any,
70 arity: BuiltinParamArity::Optional,
71 default: None,
72 description: "Final time scalar or explicit time vector.",
73 },
74];
75const STEP_SIGNATURES: [BuiltinSignatureDescriptor; 6] = [
76 BuiltinSignatureDescriptor {
77 label: "y = step(sys)",
78 inputs: &STEP_INPUTS_SYS,
79 outputs: &STEP_OUTPUT_Y,
80 },
81 BuiltinSignatureDescriptor {
82 label: "y = step(sys, tFinal)",
83 inputs: &STEP_INPUTS_SYS_TIME,
84 outputs: &STEP_OUTPUT_Y,
85 },
86 BuiltinSignatureDescriptor {
87 label: "y = step(sys, t)",
88 inputs: &STEP_INPUTS_SYS_TIME,
89 outputs: &STEP_OUTPUT_Y,
90 },
91 BuiltinSignatureDescriptor {
92 label: "[y,t] = step(sys)",
93 inputs: &STEP_INPUTS_SYS,
94 outputs: &STEP_OUTPUT_Y_T,
95 },
96 BuiltinSignatureDescriptor {
97 label: "[y,t] = step(sys, tFinal)",
98 inputs: &STEP_INPUTS_SYS_TIME,
99 outputs: &STEP_OUTPUT_Y_T,
100 },
101 BuiltinSignatureDescriptor {
102 label: "[y,t] = step(sys, t)",
103 inputs: &STEP_INPUTS_SYS_TIME,
104 outputs: &STEP_OUTPUT_Y_T,
105 },
106];
107const STEP_ERROR_INVALID_ARGUMENT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
108 code: "RM.STEP.INVALID_ARGUMENT",
109 identifier: Some("RunMat:step:InvalidArgument"),
110 when: "Inputs do not match supported step invocation forms.",
111 message: "step: invalid argument",
112};
113const STEP_ERROR_INVALID_MODEL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
114 code: "RM.STEP.INVALID_MODEL",
115 identifier: Some("RunMat:step:InvalidModel"),
116 when: "Input system is not a supported tf object with valid required properties.",
117 message: "step: invalid model",
118};
119const STEP_ERROR_INVALID_TIME: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
120 code: "RM.STEP.INVALID_TIME",
121 identifier: Some("RunMat:step:InvalidTime"),
122 when: "Time argument is invalid for the model class or sampling mode.",
123 message: "step: invalid time input",
124};
125const STEP_ERROR_UNSUPPORTED_MODEL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
126 code: "RM.STEP.UNSUPPORTED_MODEL",
127 identifier: Some("RunMat:step:UnsupportedModel"),
128 when: "Model is well-formed but unsupported by the current step implementation.",
129 message: "step: unsupported model",
130};
131const STEP_ERROR_DISCRETE_LIMIT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
132 code: "RM.STEP.DISCRETE_LIMIT",
133 identifier: Some("RunMat:step:DiscreteLimit"),
134 when: "Discrete simulation would exceed platform or configured sample limits.",
135 message: "step: discrete simulation limit exceeded",
136};
137const STEP_ERROR_PLOT_FAILED: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
138 code: "RM.STEP.PLOT_FAILED",
139 identifier: Some("RunMat:step:PlotFailed"),
140 when: "Statement-form plotting failed for reasons other than known nonfatal setup conditions.",
141 message: "step: plotting failed",
142};
143const STEP_ERROR_INTERNAL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
144 code: "RM.STEP.INTERNAL",
145 identifier: Some("RunMat:step:Internal"),
146 when: "Internal response assembly failed.",
147 message: "step: internal error",
148};
149const STEP_ERRORS: [BuiltinErrorDescriptor; 7] = [
150 STEP_ERROR_INVALID_ARGUMENT,
151 STEP_ERROR_INVALID_MODEL,
152 STEP_ERROR_INVALID_TIME,
153 STEP_ERROR_UNSUPPORTED_MODEL,
154 STEP_ERROR_DISCRETE_LIMIT,
155 STEP_ERROR_PLOT_FAILED,
156 STEP_ERROR_INTERNAL,
157];
158pub const STEP_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
159 signatures: &STEP_SIGNATURES,
160 output_mode: BuiltinOutputMode::ByRequestedOutputCount,
161 completion_policy: BuiltinCompletionPolicy::Public,
162 errors: &STEP_ERRORS,
163};
164
165const STEP_INTEGER_TIME_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
166 id: "step-integer-time",
167 mode: BuiltinExtensionMode::RunMatOnly,
168 description: "step with a native typed-integer time input is a RunMat extension",
169 error_identifier: Some("RunMat:compatibility:StepIntegerTimeExtension"),
170};
171pub const STEP_EXTENSIONS: [BuiltinExtensionDescriptor; 1] = [STEP_INTEGER_TIME_EXTENSION];
172const STEP_INTEGER_TIME_INPUTS: [BuiltinIntegerInputCapability; 1] =
173 [BuiltinIntegerInputCapability {
174 name: "t or tFinal",
175 classes: &crate::builtins::common::integer_capability::ALL_INTEGER_CLASSES,
176 availability: BuiltinIntegerInputAvailability::RunMatOnly,
177 scalar_double: BuiltinIntegerScalarDoubleRule::Allowed,
178 notes: "The compatibility target documents positive scalar, two-element, and vector time forms without publishing native integer storage classes. RunMat conservatively gates typed integers and requires exact binary64 representation.",
179 }];
180pub const STEP_INTEGER_CAPABILITIES: [BuiltinIntegerCapabilityDescriptor; 1] =
181 [BuiltinIntegerCapabilityDescriptor {
182 form: "[y,tOut] = step(sys, integer_t_or_tFinal)",
183 inputs: &STEP_INTEGER_TIME_INPUTS,
184 computation_domain: BuiltinIntegerComputationDomain::FloatingPoint,
185 output_class: BuiltinIntegerOutputClassRule::Double,
186 overflow: BuiltinIntegerOverflowRule::Error,
187 backend: BuiltinIntegerBackendRule::GatherFallback,
188 overload: BuiltinIntegerOverloadKind::Multiple,
189 notes: "Compatibility admission and exactness checks occur before provider access. Automatic residency may gather through the exact owner; the host simulator and time/output arrays use the model's binary64 computation domain.",
190 }];
191
192#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::control::step")]
193pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
194 name: "step",
195 op_kind: GpuOpKind::Custom("control-step-response"),
196 supported_precisions: &[],
197 broadcast: BroadcastSemantics::None,
198 provider_hooks: &[],
199 constant_strategy: ConstantStrategy::InlineLiteral,
200 residency: ResidencyPolicy::GatherImmediately,
201 nan_mode: ReductionNaN::Include,
202 two_pass_threshold: None,
203 workgroup_size: None,
204 accepts_nan_mode: false,
205 notes: "Step-response simulation runs on the host from transfer-function metadata.",
206};
207
208#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::control::step")]
209pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
210 name: "step",
211 shape: ShapeRequirements::Any,
212 constant_strategy: ConstantStrategy::InlineLiteral,
213 elementwise: None,
214 reduction: None,
215 emits_nan: false,
216 notes: "step simulates a dynamic system and terminates numeric fusion chains.",
217};
218
219fn step_error_with_detail(
220 error: &'static BuiltinErrorDescriptor,
221 detail: impl AsRef<str>,
222) -> RuntimeError {
223 step_error_with_message(format!("{}: {}", error.message, detail.as_ref()), error)
224}
225
226fn step_error_with_message(
227 message: impl Into<String>,
228 error: &'static BuiltinErrorDescriptor,
229) -> RuntimeError {
230 let mut builder = build_runtime_error(message).with_builtin(BUILTIN_NAME);
231 if let Some(identifier) = error.identifier {
232 builder = builder.with_identifier(identifier);
233 }
234 builder.build()
235}
236
237#[runtime_builtin(
238 name = "step",
239 category = "control",
240 summary = "Compute or plot step responses of SISO transfer-function models.",
241 keywords = "step,response,control system,transfer function,tf",
242 sink = true,
243 suppress_auto_output = true,
244 type_resolver(step_type),
245 descriptor(crate::builtins::control::step::STEP_DESCRIPTOR),
246 extensions(crate::builtins::control::step::STEP_EXTENSIONS),
247 integer_capabilities(crate::builtins::control::step::STEP_INTEGER_CAPABILITIES),
248 builtin_path = "crate::builtins::control::step"
249)]
250async fn step_builtin(sys: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
251 for value in &rest {
252 crate::builtins::common::validation::ensure_runmat_integer_f64_boundary(
253 value,
254 &STEP_INTEGER_TIME_EXTENSION,
255 BUILTIN_NAME,
256 "time",
257 )
258 .await?;
259 }
260 if is_statement_form_call() {
261 plot_multiple_step_responses(sys, rest).await?;
262 return Ok(Value::OutputList(Vec::new()));
263 }
264
265 if rest.len() > 1 {
266 return Err(step_error_with_detail(
267 &STEP_ERROR_INVALID_ARGUMENT,
268 "expected step(sys), step(sys, tFinal), or step(sys, t)",
269 ));
270 }
271
272 let sys = crate::gather_if_needed_async(&sys).await?;
273 let mut rest_host = Vec::with_capacity(rest.len());
274 for arg in &rest {
275 rest_host.push(crate::gather_if_needed_async(arg).await?);
276 }
277 let model = TransferFunction::from_value(sys)?;
278 let time = TimeSpec::parse(rest_host.first(), model.sample_time)?;
279 let eval = evaluate_step(&model, time)?;
280
281 if crate::output_context::requested_output_count() == Some(0)
282 && crate::output_count::current_output_count().is_none()
283 {
284 plot_response(&eval).await?;
285 return Ok(Value::OutputList(Vec::new()));
286 }
287
288 if let Some(out_count) = crate::output_count::current_output_count() {
289 if out_count == 0 {
290 plot_response(&eval).await?;
291 return Ok(Value::OutputList(Vec::new()));
292 }
293 if out_count == 1 {
294 return Ok(Value::OutputList(vec![eval.y_value()?]));
295 }
296 return Ok(crate::output_count::output_list_with_padding(
297 out_count,
298 eval.outputs()?,
299 ));
300 }
301
302 eval.y_value()
303}
304
305fn is_statement_form_call() -> bool {
306 matches!(crate::output_count::current_output_count(), Some(0))
307 || (crate::output_context::requested_output_count() == Some(0)
308 && crate::output_count::current_output_count().is_none())
309}
310
311async fn plot_multiple_step_responses(first_sys: Value, rest: Vec<Value>) -> BuiltinResult<()> {
312 let first_sys = crate::gather_if_needed_async(&first_sys).await?;
313 let mut systems = vec![(first_sys, None)];
314 let mut time_arg = None;
315 for arg in rest {
316 let gathered = crate::gather_if_needed_async(&arg).await?;
317 if is_plot_style_arg(&gathered) {
318 if let Some((_, style)) = systems.last_mut() {
319 if style.is_some() {
320 return Err(step_error_with_detail(
321 &STEP_ERROR_INVALID_ARGUMENT,
322 "only one style argument is supported per system",
323 ));
324 }
325 *style = Some(gathered);
326 continue;
327 }
328 continue;
329 }
330 if is_tf_object(&gathered) {
331 if time_arg.is_some() {
332 return Err(step_error_with_detail(
333 &STEP_ERROR_INVALID_ARGUMENT,
334 "time argument must follow all systems in statement-form step plots",
335 ));
336 }
337 systems.push((gathered, None));
338 continue;
339 }
340 if time_arg.is_none() {
341 time_arg = Some(gathered);
342 continue;
343 }
344 return Err(step_error_with_detail(
345 &STEP_ERROR_INVALID_ARGUMENT,
346 "unsupported statement-form step plot argument",
347 ));
348 }
349
350 if systems.is_empty() {
351 return Err(step_error_with_detail(
352 &STEP_ERROR_INVALID_ARGUMENT,
353 "at least one system is required",
354 ));
355 }
356
357 let mut first = true;
358 let mut hold_enabled = false;
359 let result: BuiltinResult<()> = async {
360 for (system, style) in systems {
361 let model = TransferFunction::from_value(system)?;
362 let time = TimeSpec::parse(time_arg.as_ref(), model.sample_time)?;
363 let eval = evaluate_step(&model, time)?;
364 plot_response_with_style(&eval, style.as_ref()).await?;
365 if first {
366 first = false;
367 let _ = crate::call_builtin_async("hold", &[Value::from("on")]).await;
368 hold_enabled = true;
369 }
370 }
371 Ok(())
372 }
373 .await;
374 if hold_enabled {
375 let _ = crate::call_builtin_async("hold", &[Value::from("off")]).await;
376 }
377 result
378}
379
380fn is_plot_style_arg(value: &Value) -> bool {
381 matches!(
382 value,
383 Value::String(_) | Value::StringArray(_) | Value::CharArray(_)
384 )
385}
386
387fn is_tf_object(value: &Value) -> bool {
388 matches!(value, Value::Object(object) if object.is_class("tf"))
389}
390
391#[derive(Clone, Debug)]
392struct TransferFunction {
393 numerator: Vec<f64>,
394 denominator: Vec<f64>,
395 sample_time: f64,
396}
397
398impl TransferFunction {
399 fn from_value(value: Value) -> BuiltinResult<Self> {
400 let Value::Object(object) = value else {
401 return Err(step_error_with_detail(
402 &STEP_ERROR_INVALID_MODEL,
403 "expected a tf object",
404 ));
405 };
406 if !object.is_class("tf") {
407 return Err(step_error_with_detail(
408 &STEP_ERROR_INVALID_MODEL,
409 format!("expected a tf object, got {}", object.class_name),
410 ));
411 }
412
413 let numerator = property_coefficients(&object, "Numerator")?;
414 let denominator = property_coefficients(&object, "Denominator")?;
415 let sample_time = property_scalar(&object, "Ts")?;
416 if !sample_time.is_finite() || sample_time < 0.0 {
417 return Err(step_error_with_detail(
418 &STEP_ERROR_INVALID_MODEL,
419 "tf sample time must be finite and non-negative",
420 ));
421 }
422
423 let numerator = trim_leading_zeros(numerator);
424 let denominator = trim_leading_zeros(denominator);
425 if denominator.is_empty() {
426 return Err(step_error_with_detail(
427 &STEP_ERROR_INVALID_MODEL,
428 "denominator coefficients cannot be empty",
429 ));
430 }
431 if denominator[0].abs() <= EPS {
432 return Err(step_error_with_detail(
433 &STEP_ERROR_INVALID_MODEL,
434 "leading denominator coefficient must be non-zero",
435 ));
436 }
437 if numerator.len().saturating_sub(1) > denominator.len().saturating_sub(1) {
438 return Err(step_error_with_detail(
439 &STEP_ERROR_UNSUPPORTED_MODEL,
440 "improper transfer functions are not supported yet",
441 ));
442 }
443
444 Ok(Self {
445 numerator,
446 denominator,
447 sample_time,
448 })
449 }
450
451 fn normalized(&self) -> (Vec<f64>, Vec<f64>) {
452 let leading = self.denominator[0];
453 let den = self
454 .denominator
455 .iter()
456 .map(|value| value / leading)
457 .collect::<Vec<_>>();
458 let num = self
459 .numerator
460 .iter()
461 .map(|value| value / leading)
462 .collect::<Vec<_>>();
463 (num, den)
464 }
465}
466
467fn property_coefficients(object: &ObjectInstance, name: &str) -> BuiltinResult<Vec<f64>> {
468 let value = object.properties.get(name).ok_or_else(|| {
469 step_error_with_detail(
470 &STEP_ERROR_INVALID_MODEL,
471 format!("tf object is missing {name}"),
472 )
473 })?;
474 match value {
475 Value::Tensor(tensor) => Ok(tensor::tensor_values_f64(tensor)),
476 Value::ComplexTensor(tensor) => real_complex_coefficients(tensor, name),
477 Value::Num(n) => Ok(vec![*n]),
478 Value::Int(i) => Ok(vec![i.to_f64()]),
479 other => Err(step_error_with_detail(
480 &STEP_ERROR_INVALID_MODEL,
481 format!("tf {name} coefficients must be numeric, got {other:?}"),
482 )),
483 }
484}
485
486fn real_complex_coefficients(tensor: &ComplexTensor, name: &str) -> BuiltinResult<Vec<f64>> {
487 tensor::complex_tensor_values_complex64(tensor)
488 .into_iter()
489 .map(|value| {
490 if value.im.abs() > EPS {
491 return Err(step_error_with_detail(
492 &STEP_ERROR_UNSUPPORTED_MODEL,
493 format!("complex tf {name} coefficients are not supported yet"),
494 ));
495 }
496 Ok(value.re)
497 })
498 .collect()
499}
500
501fn property_scalar(object: &ObjectInstance, name: &str) -> BuiltinResult<f64> {
502 let value = object.properties.get(name).ok_or_else(|| {
503 step_error_with_detail(
504 &STEP_ERROR_INVALID_MODEL,
505 format!("tf object is missing {name}"),
506 )
507 })?;
508 match value {
509 Value::Num(n) => Ok(*n),
510 Value::Int(i) => Ok(i.to_f64()),
511 Value::Tensor(tensor) if tensor::is_scalar_tensor(tensor) => {
512 Ok(tensor::tensor_value_f64(tensor, 0))
513 }
514 other => Err(step_error_with_detail(
515 &STEP_ERROR_INVALID_MODEL,
516 format!("tf {name} property must be a scalar, got {other:?}"),
517 )),
518 }
519}
520
521#[derive(Clone, Debug)]
522enum TimeSpec {
523 Auto,
524 FinalTime(f64),
525 Vector(Vec<f64>),
526}
527
528impl TimeSpec {
529 fn parse(value: Option<&Value>, sample_time: f64) -> BuiltinResult<Self> {
530 let Some(value) = value else {
531 return Ok(Self::Auto);
532 };
533 match value {
534 Value::Num(n) => Self::final_time(*n),
535 Value::Int(i) => Self::final_time(i.to_f64()),
536 Value::Tensor(tensor) => {
537 ensure_time_vector_shape(&tensor.shape)?;
538 let data = tensor::tensor_values_f64(tensor);
539 if tensor::is_scalar_tensor(tensor) {
540 return Self::final_time(tensor::tensor_value_f64(tensor, 0));
541 }
542 Self::vector(data, sample_time)
543 }
544 other => Err(step_error_with_detail(
545 &STEP_ERROR_INVALID_TIME,
546 format!("time input must be a scalar final time or numeric vector, got {other:?}"),
547 )),
548 }
549 }
550
551 fn final_time(value: f64) -> BuiltinResult<Self> {
552 if !value.is_finite() || value <= 0.0 {
553 return Err(step_error_with_detail(
554 &STEP_ERROR_INVALID_TIME,
555 "final time must be a positive finite scalar",
556 ));
557 }
558 Ok(Self::FinalTime(value))
559 }
560
561 fn vector(values: Vec<f64>, sample_time: f64) -> BuiltinResult<Self> {
562 validate_time_vector(&values)?;
563 if sample_time > 0.0 {
564 for &t in &values {
565 let k = (t / sample_time).round();
566 if (t - k * sample_time).abs() > 1.0e-8 * sample_time.max(1.0) {
567 return Err(step_error_with_detail(
568 &STEP_ERROR_INVALID_TIME,
569 "discrete-time sample vector must align with the model sample time",
570 ));
571 }
572 }
573 }
574 Ok(Self::Vector(values))
575 }
576}
577
578fn ensure_time_vector_shape(shape: &[usize]) -> BuiltinResult<()> {
579 let non_unit = shape.iter().copied().filter(|&dim| dim > 1).count();
580 if non_unit <= 1 {
581 Ok(())
582 } else {
583 Err(step_error_with_detail(
584 &STEP_ERROR_INVALID_TIME,
585 "time input must be a vector",
586 ))
587 }
588}
589
590fn validate_time_vector(values: &[f64]) -> BuiltinResult<()> {
591 if values.is_empty() {
592 return Err(step_error_with_detail(
593 &STEP_ERROR_INVALID_TIME,
594 "time vector must not be empty",
595 ));
596 }
597 let mut previous = None;
598 for &value in values {
599 if !value.is_finite() || value < 0.0 {
600 return Err(step_error_with_detail(
601 &STEP_ERROR_INVALID_TIME,
602 "time vector values must be finite and non-negative",
603 ));
604 }
605 if let Some(prev) = previous {
606 if value < prev {
607 return Err(step_error_with_detail(
608 &STEP_ERROR_INVALID_TIME,
609 "time vector must be nondecreasing",
610 ));
611 }
612 }
613 previous = Some(value);
614 }
615 Ok(())
616}
617
618#[derive(Clone, Debug)]
619struct StepEval {
620 y: Vec<f64>,
621 t: Vec<f64>,
622}
623
624impl StepEval {
625 fn y_value(&self) -> BuiltinResult<Value> {
626 column_tensor(self.y.clone())
627 }
628
629 fn t_value(&self) -> BuiltinResult<Value> {
630 column_tensor(self.t.clone())
631 }
632
633 fn outputs(&self) -> BuiltinResult<Vec<Value>> {
634 Ok(vec![self.y_value()?, self.t_value()?])
635 }
636}
637
638fn evaluate_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
639 if model.sample_time > 0.0 {
640 evaluate_discrete_step(model, time)
641 } else {
642 evaluate_continuous_step(model, time)
643 }
644}
645
646fn evaluate_continuous_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
647 let t = continuous_time_vector(model, time)?;
648 let (num, den) = model.normalized();
649 let response = continuous_response(&num, &den, &t)?;
650 Ok(StepEval { y: response, t })
651}
652
653fn continuous_time_vector(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<Vec<f64>> {
654 match time {
655 TimeSpec::Auto => Ok(linspace(0.0, automatic_final_time(model), DEFAULT_SAMPLES)),
656 TimeSpec::FinalTime(final_time) => Ok(linspace(0.0, final_time, DEFAULT_SAMPLES)),
657 TimeSpec::Vector(values) => Ok(values),
658 }
659}
660
661fn continuous_response(num: &[f64], den: &[f64], t: &[f64]) -> BuiltinResult<Vec<f64>> {
662 let order = den.len() - 1;
663 if order == 0 {
664 let gain = num.last().copied().unwrap_or(0.0) / den[0];
665 return Ok(vec![gain; t.len()]);
666 }
667
668 let mut padded_num = vec![0.0; order + 1 - num.len()];
669 padded_num.extend_from_slice(num);
670 let direct = padded_num[0];
671 let a = &den[1..];
672 let mut c = Vec::with_capacity(order);
673 for state_idx in 0..order {
674 let coeff_idx = order - state_idx;
675 c.push(padded_num[coeff_idx] - direct * den[coeff_idx]);
676 }
677
678 let mut state = vec![0.0; order];
679 let mut current_t = 0.0;
680 let mut response = Vec::with_capacity(t.len());
681 for &target_t in t {
682 if target_t > current_t {
683 integrate_to(&mut state, a, current_t, target_t);
684 current_t = target_t;
685 }
686 response.push(dot(&c, &state) + direct);
687 }
688 Ok(response)
689}
690
691fn integrate_to(state: &mut [f64], a: &[f64], start: f64, end: f64) {
692 let duration = end - start;
693 if duration <= 0.0 {
694 return;
695 }
696 let steps = ((duration / 0.01).ceil() as usize).clamp(1, 10_000);
697 let h = duration / steps as f64;
698 for _ in 0..steps {
699 rk4_step(state, a, h);
700 }
701}
702
703fn rk4_step(state: &mut [f64], a: &[f64], h: f64) {
704 let k1 = derivative(state, a);
705 let s2 = add_scaled(state, &k1, h * 0.5);
706 let k2 = derivative(&s2, a);
707 let s3 = add_scaled(state, &k2, h * 0.5);
708 let k3 = derivative(&s3, a);
709 let s4 = add_scaled(state, &k3, h);
710 let k4 = derivative(&s4, a);
711 for idx in 0..state.len() {
712 state[idx] += h * (k1[idx] + 2.0 * k2[idx] + 2.0 * k3[idx] + k4[idx]) / 6.0;
713 }
714}
715
716fn derivative(state: &[f64], a: &[f64]) -> Vec<f64> {
717 let order = state.len();
718 let mut dx = vec![0.0; order];
719 if order > 1 {
720 dx[..(order - 1)].copy_from_slice(&state[1..order]);
721 }
722 let mut last = 1.0;
723 for state_idx in 0..order {
724 let coeff = a[order - 1 - state_idx];
725 last -= coeff * state[state_idx];
726 }
727 dx[order - 1] = last;
728 dx
729}
730
731fn add_scaled(state: &[f64], delta: &[f64], scale: f64) -> Vec<f64> {
732 state
733 .iter()
734 .zip(delta)
735 .map(|(value, delta)| value + scale * delta)
736 .collect()
737}
738
739fn evaluate_discrete_step(model: &TransferFunction, time: TimeSpec) -> BuiltinResult<StepEval> {
740 let t = discrete_time_vector(model.sample_time, time)?;
741 if t.len() > MAX_DISCRETE_SAMPLES {
742 return Err(step_error_with_detail(
743 &STEP_ERROR_DISCRETE_LIMIT,
744 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
745 ));
746 }
747 let sample_indices = t
748 .iter()
749 .map(|&value| checked_discrete_sample_index(model.sample_time, value))
750 .collect::<BuiltinResult<Vec<_>>>()?;
751 let max_k = sample_indices.iter().copied().max().unwrap_or(0);
752 let count = max_k.checked_add(1).ok_or_else(|| {
753 step_error_with_detail(
754 &STEP_ERROR_DISCRETE_LIMIT,
755 "discrete sample index exceeds platform limits",
756 )
757 })?;
758 let (num, den) = model.normalized();
759 let all_y = discrete_response(&num, &den, count)?;
760 let y = sample_indices
761 .into_iter()
762 .map(|idx| {
763 all_y.get(idx).copied().ok_or_else(|| {
764 step_error_with_detail(
765 &STEP_ERROR_DISCRETE_LIMIT,
766 "discrete sample index exceeds response length",
767 )
768 })
769 })
770 .collect::<BuiltinResult<Vec<_>>>()?;
771 Ok(StepEval { y, t })
772}
773
774fn discrete_time_vector(sample_time: f64, time: TimeSpec) -> BuiltinResult<Vec<f64>> {
775 match time {
776 TimeSpec::Auto => Ok((0..DEFAULT_SAMPLES)
777 .map(|idx| idx as f64 * sample_time)
778 .collect()),
779 TimeSpec::FinalTime(final_time) => {
780 let steps = checked_discrete_sample_steps(sample_time, final_time)?;
781 Ok((0..=steps).map(|idx| idx as f64 * sample_time).collect())
782 }
783 TimeSpec::Vector(values) => Ok(values),
784 }
785}
786
787fn checked_discrete_sample_steps(sample_time: f64, final_time: f64) -> BuiltinResult<usize> {
788 let steps = (final_time / sample_time).floor();
789 if !steps.is_finite() || steps < 0.0 || steps > usize::MAX as f64 {
790 return Err(step_error_with_detail(
791 &STEP_ERROR_DISCRETE_LIMIT,
792 "discrete sample count exceeds platform limits",
793 ));
794 }
795 if steps >= MAX_DISCRETE_SAMPLES as f64 {
796 return Err(step_error_with_detail(
797 &STEP_ERROR_DISCRETE_LIMIT,
798 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
799 ));
800 }
801 Ok(steps as usize)
802}
803
804fn checked_discrete_sample_index(sample_time: f64, time: f64) -> BuiltinResult<usize> {
805 let index = (time / sample_time).round();
806 if !index.is_finite() || index < 0.0 || index > usize::MAX as f64 {
807 return Err(step_error_with_detail(
808 &STEP_ERROR_DISCRETE_LIMIT,
809 "discrete sample index exceeds platform limits",
810 ));
811 }
812 if index >= MAX_DISCRETE_SAMPLES as f64 {
813 return Err(step_error_with_detail(
814 &STEP_ERROR_DISCRETE_LIMIT,
815 format!("discrete response would require more than {MAX_DISCRETE_SAMPLES} samples"),
816 ));
817 }
818 Ok(index as usize)
819}
820
821fn discrete_response(num: &[f64], den: &[f64], count: usize) -> BuiltinResult<Vec<f64>> {
822 let order = den.len() - 1;
823 let mut padded_num = vec![0.0; order + 1 - num.len()];
824 padded_num.extend_from_slice(num);
825 let mut y = vec![0.0; count];
826 for k in 0..count {
827 let mut value = 0.0;
828 for (idx, &coeff) in padded_num.iter().enumerate() {
829 if k >= idx {
830 value += coeff;
831 }
832 }
833 for idx in 1..den.len() {
834 if k >= idx {
835 value -= den[idx] * y[k - idx];
836 }
837 }
838 y[k] = value;
839 }
840 Ok(y)
841}
842
843async fn plot_response(eval: &StepEval) -> BuiltinResult<()> {
844 plot_response_with_style(eval, None).await
845}
846
847async fn plot_response_with_style(eval: &StepEval, style: Option<&Value>) -> BuiltinResult<()> {
848 let t = eval.t_value()?;
849 let y = eval.y_value()?;
850 let args = if let Some(style) = style {
851 vec![t, y, style.clone()]
852 } else {
853 vec![t, y]
854 };
855 if let Err(err) = crate::call_builtin_async("plot", &args).await {
856 if super::is_nonfatal_plot_setup_error(&err) {
857 return Ok(());
858 }
859 return Err(step_error_with_detail(
860 &STEP_ERROR_PLOT_FAILED,
861 err.message(),
862 ));
863 }
864 let _ = crate::call_builtin_async("title", &[Value::from("Step Response")]).await;
865 let _ = crate::call_builtin_async("xlabel", &[Value::from("Time")]).await;
866 let _ = crate::call_builtin_async("ylabel", &[Value::from("Amplitude")]).await;
867 Ok(())
868}
869
870fn automatic_final_time(model: &TransferFunction) -> f64 {
871 let (_, den) = model.normalized();
872 let poles = polynomial_roots(&den).unwrap_or_default();
873 let slowest_decay = poles
874 .iter()
875 .filter_map(|pole| if pole.re < -EPS { Some(-pole.re) } else { None })
876 .fold(f64::INFINITY, f64::min);
877 if slowest_decay.is_finite() && slowest_decay > EPS {
878 (5.0 / slowest_decay).clamp(1.0, 100.0)
879 } else {
880 10.0
881 }
882}
883
884fn polynomial_roots(coeffs: &[f64]) -> BuiltinResult<Vec<Complex64>> {
885 let trimmed = trim_leading_zeros(coeffs.to_vec());
886 if trimmed.len() <= 1 {
887 return Ok(Vec::new());
888 }
889 if trimmed.len() == 2 {
890 return Ok(vec![Complex64::new(-trimmed[1] / trimmed[0], 0.0)]);
891 }
892 let degree = trimmed.len() - 1;
893 let leading = trimmed[0];
894 let mut companion = DMatrix::<Complex64>::zeros(degree, degree);
895 for row in 1..degree {
896 companion[(row, row - 1)] = Complex64::new(1.0, 0.0);
897 }
898 for (idx, coeff) in trimmed.iter().enumerate().skip(1) {
899 companion[(0, idx - 1)] = Complex64::new(-coeff / leading, 0.0);
900 }
901 let eigenvalues = companion.eigenvalues().ok_or_else(|| {
902 step_error_with_detail(
903 &STEP_ERROR_INTERNAL,
904 "failed to compute transfer-function poles",
905 )
906 })?;
907 Ok(eigenvalues.iter().copied().collect())
908}
909
910fn trim_leading_zeros(coeffs: Vec<f64>) -> Vec<f64> {
911 let first_nonzero = coeffs
912 .iter()
913 .position(|value| value.abs() > EPS)
914 .unwrap_or(coeffs.len());
915 coeffs[first_nonzero..].to_vec()
916}
917
918fn linspace(start: f64, end: f64, count: usize) -> Vec<f64> {
919 if count <= 1 {
920 return vec![end];
921 }
922 let step = (end - start) / (count - 1) as f64;
923 (0..count).map(|idx| start + idx as f64 * step).collect()
924}
925
926fn dot(lhs: &[f64], rhs: &[f64]) -> f64 {
927 lhs.iter().zip(rhs).map(|(a, b)| a * b).sum()
928}
929
930fn column_tensor(data: Vec<f64>) -> BuiltinResult<Value> {
931 let rows = data.len();
932 let tensor = Tensor::new(data, vec![rows, 1]).map_err(|err| {
933 step_error_with_detail(
934 &STEP_ERROR_INTERNAL,
935 format!("failed to build response tensor: {err}"),
936 )
937 })?;
938 Ok(Value::Tensor(tensor))
939}
940
941#[cfg(test)]
942mod tests {
943 use super::*;
944 use futures::executor::block_on;
945 use runmat_value::{CharArray, IntegerComplexStorage, IntegerStorage, ObjectInstance};
946
947 fn tf_object(num: Vec<f64>, den: Vec<f64>, sample_time: f64) -> Value {
948 let mut object = ObjectInstance::new("tf".to_string());
949 object.properties.insert(
950 "Numerator".to_string(),
951 Value::Tensor(Tensor::new(num.clone(), vec![1, num.len()]).unwrap()),
952 );
953 object.properties.insert(
954 "Denominator".to_string(),
955 Value::Tensor(Tensor::new(den.clone(), vec![1, den.len()]).unwrap()),
956 );
957 object.properties.insert(
958 "Variable".to_string(),
959 Value::CharArray(CharArray::new_row(if sample_time > 0.0 {
960 "z"
961 } else {
962 "s"
963 })),
964 );
965 object
966 .properties
967 .insert("Ts".to_string(), Value::Num(sample_time));
968 Value::Object(object)
969 }
970
971 fn run_step(sys: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
972 block_on(step_builtin(sys, rest))
973 }
974
975 fn tensor_data(value: Value) -> Vec<f64> {
976 match value {
977 Value::Tensor(tensor) => tensor.materialize_f64(),
978 other => panic!("expected tensor, got {other:?}"),
979 }
980 }
981
982 fn integer_tensor(storage: IntegerStorage, shape: Vec<usize>) -> Value {
983 Value::Tensor(Tensor::new_integer(storage, shape).expect("integer tensor"))
984 }
985
986 fn poisoned_complex_integer_tensor(
987 real: IntegerStorage,
988 imag: IntegerStorage,
989 shape: Vec<usize>,
990 ) -> ComplexTensor {
991 let storage = IntegerComplexStorage::new(real, imag).expect("complex integer storage");
992
993 ComplexTensor::new_integer(storage, shape).expect("complex integer tensor")
994 }
995
996 #[test]
997 fn step_descriptor_signatures_cover_core_forms() {
998 let labels: Vec<&str> = STEP_DESCRIPTOR
999 .signatures
1000 .iter()
1001 .map(|sig| sig.label)
1002 .collect();
1003 assert!(labels.contains(&"y = step(sys)"));
1004 assert!(labels.contains(&"y = step(sys, tFinal)"));
1005 assert!(labels.contains(&"y = step(sys, t)"));
1006 assert!(labels.contains(&"[y,t] = step(sys)"));
1007 assert!(labels.contains(&"[y,t] = step(sys, tFinal)"));
1008 assert!(labels.contains(&"[y,t] = step(sys, t)"));
1009 }
1010
1011 #[test]
1012 fn first_order_continuous_response_matches_closed_form_for_explicit_time() {
1013 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
1014 let time = Value::Tensor(Tensor::new(vec![0.0, 0.5, 1.0, 2.0], vec![1, 4]).unwrap());
1015 let y = tensor_data(run_step(sys, vec![time]).expect("step"));
1016 for (actual, t) in y.iter().zip([0.0_f64, 0.5, 1.0, 2.0]) {
1017 let expected = 1.0 - (-t).exp();
1018 assert!(
1019 (actual - expected).abs() < 1.0e-5,
1020 "t={t} actual={actual} expected={expected}"
1021 );
1022 }
1023 }
1024
1025 #[test]
1026 fn multi_output_returns_y_then_time() {
1027 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
1028 let _guard = crate::output_count::push_output_count(Some(2));
1029 let time = Value::Tensor(Tensor::new(vec![0.0, 1.0], vec![2, 1]).unwrap());
1030 let result = run_step(sys, vec![time]).expect("step");
1031 let Value::OutputList(outputs) = result else {
1032 panic!("expected output list");
1033 };
1034 assert_eq!(outputs.len(), 2);
1035 assert_eq!(tensor_data(outputs[1].clone()), vec![0.0, 1.0]);
1036 let y = tensor_data(outputs[0].clone());
1037 assert!((y[1] - (1.0 - (-1.0_f64).exp())).abs() < 1.0e-5);
1038 }
1039
1040 #[test]
1041 fn single_requested_output_returns_only_response() {
1042 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
1043 let _guard = crate::output_count::push_output_count(Some(1));
1044 let time = Value::Tensor(Tensor::new(vec![0.0, 1.0], vec![1, 2]).unwrap());
1045 let result = run_step(sys, vec![time]).expect("step");
1046 let Value::OutputList(outputs) = result else {
1047 panic!("expected output list");
1048 };
1049 assert_eq!(outputs.len(), 1);
1050 let y = tensor_data(outputs[0].clone());
1051 assert_eq!(y.len(), 2);
1052 assert!((y[1] - (1.0 - (-1.0_f64).exp())).abs() < 1.0e-5);
1053 }
1054
1055 #[test]
1056 fn scalar_final_time_generates_column_time_vector_ending_at_final_time() {
1057 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
1058 let _guard = crate::output_count::push_output_count(Some(2));
1059 let result = run_step(sys, vec![Value::Num(5.0)]).expect("step");
1060 let Value::OutputList(outputs) = result else {
1061 panic!("expected output list");
1062 };
1063 let t = tensor_data(outputs[1].clone());
1064 assert_eq!(t.len(), DEFAULT_SAMPLES);
1065 assert_eq!(t[0], 0.0);
1066 assert!((t[t.len() - 1] - 5.0).abs() < 1.0e-12);
1067 }
1068
1069 #[test]
1070 fn discrete_response_uses_sample_time_grid() {
1071 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 0.1);
1072 let time = Value::Tensor(Tensor::new(vec![0.0, 0.1, 0.2], vec![1, 3]).unwrap());
1073 let y = tensor_data(run_step(sys, vec![time]).expect("step"));
1074 assert_eq!(y, vec![0.0, 1.0, 1.5]);
1075 }
1076
1077 #[test]
1078 fn step_typed_integer_time_inputs_cross_double_boundary_exactly() {
1079 let _compat = crate::compatibility::push_runmat_extensions_enabled(true);
1080 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
1081 let _guard = crate::output_count::push_output_count(Some(2));
1082 let result = run_step(
1083 sys,
1084 vec![integer_tensor(
1085 IntegerStorage::U64(vec![0, 1, 2]),
1086 vec![1, 3],
1087 )],
1088 )
1089 .expect("step");
1090 let Value::OutputList(outputs) = result else {
1091 panic!("expected output list");
1092 };
1093 assert_eq!(tensor_data(outputs[1].clone()), vec![0.0, 1.0, 2.0]);
1094 }
1095
1096 #[test]
1097 fn step_scalar_final_time_reads_typed_integer_storage_length_exactly() {
1098 let _compat = crate::compatibility::push_runmat_extensions_enabled(true);
1099 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
1100 let final_time =
1101 Tensor::new_integer(IntegerStorage::U16(vec![2]), vec![1, 1]).expect("final time");
1102 let _guard = crate::output_count::push_output_count(Some(2));
1103 let result = run_step(sys, vec![Value::Tensor(final_time)]).expect("step");
1104 let Value::OutputList(outputs) = result else {
1105 panic!("expected output list");
1106 };
1107 let time = tensor_data(outputs[1].clone());
1108 assert_eq!(time.first().copied(), Some(0.0));
1109 assert_eq!(time.last().copied(), Some(2.0));
1110 }
1111
1112 #[test]
1113 fn step_typed_integer_time_is_mode_gated() {
1114 let _compat = crate::compatibility::push_runmat_extensions_enabled(false);
1115 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
1116 let error = run_step(
1117 sys,
1118 vec![integer_tensor(IntegerStorage::U16(vec![2]), vec![1, 1])],
1119 )
1120 .expect_err("MATLAB-compatible mode must reject typed integer time");
1121 assert_eq!(
1122 error.identifier(),
1123 Some("RunMat:compatibility:StepIntegerTimeExtension")
1124 );
1125 }
1126
1127 #[test]
1128 fn step_rejects_wide_integer_time_before_binary64_rounding() {
1129 let _compat = crate::compatibility::push_runmat_extensions_enabled(true);
1130 let sys = tf_object(vec![1.0], vec![1.0, 1.0], 0.0);
1131 let error = run_step(
1132 sys,
1133 vec![integer_tensor(
1134 IntegerStorage::U64(vec![9_007_199_254_740_993]),
1135 vec![1, 1],
1136 )],
1137 )
1138 .expect_err("inexact binary64 conversion must fail");
1139 assert!(error.message().contains("exactly representable"));
1140 }
1141
1142 #[test]
1143 fn step_sample_time_parser_ignores_poisoned_integer_mirrors_for_all_classes() {
1144 let storages = [
1145 IntegerStorage::I8(vec![1]),
1146 IntegerStorage::I16(vec![1]),
1147 IntegerStorage::I32(vec![1]),
1148 IntegerStorage::I64(vec![1]),
1149 IntegerStorage::U8(vec![1]),
1150 IntegerStorage::U16(vec![1]),
1151 IntegerStorage::U32(vec![1]),
1152 IntegerStorage::U64(vec![1]),
1153 ];
1154
1155 for storage in storages {
1156 let sample_time = Tensor::new_integer(storage, vec![1, 1]).expect("sample time");
1157 let mut object = ObjectInstance::new("tf".to_string());
1158 object
1159 .properties
1160 .insert("Ts".to_string(), Value::Tensor(sample_time));
1161
1162 assert_eq!(property_scalar(&object, "Ts").expect("sample time"), 1.0);
1163 }
1164 }
1165
1166 #[test]
1167 fn step_tf_coefficients_read_complex_typed_integer_storage_exactly() {
1168 let tensor = poisoned_complex_integer_tensor(
1169 IntegerStorage::I16(vec![2, 4]),
1170 IntegerStorage::I16(vec![0, 0]),
1171 vec![1, 2],
1172 );
1173
1174 assert_eq!(
1175 real_complex_coefficients(&tensor, "Numerator").expect("coefficients"),
1176 vec![2.0, 4.0]
1177 );
1178 }
1179
1180 #[test]
1181 fn step_rejects_exact_nonreal_complex_typed_integer_coefficients() {
1182 let tensor = poisoned_complex_integer_tensor(
1183 IntegerStorage::I16(vec![2, 4]),
1184 IntegerStorage::I16(vec![0, 1]),
1185 vec![1, 2],
1186 );
1187
1188 let err = real_complex_coefficients(&tensor, "Numerator").expect_err("nonreal coeff");
1189 assert_eq!(err.identifier(), STEP_ERROR_UNSUPPORTED_MODEL.identifier);
1190 assert!(err.message().contains("complex tf Numerator coefficients"));
1191 }
1192
1193 #[test]
1194 fn discrete_final_time_rejects_excessive_sample_count() {
1195 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0e-6);
1196 let err = run_step(sys, vec![Value::Num(2.0)]).expect_err("should fail");
1197 assert!(err.message().contains("more than 1000000 samples"));
1198 assert_eq!(err.identifier(), STEP_ERROR_DISCRETE_LIMIT.identifier);
1199 }
1200
1201 #[test]
1202 fn discrete_time_vector_rejects_excessive_sample_index() {
1203 let sys = tf_object(vec![1.0], vec![1.0, -0.5], 1.0);
1204 let time =
1205 Value::Tensor(Tensor::new(vec![0.0, MAX_DISCRETE_SAMPLES as f64], vec![1, 2]).unwrap());
1206 let err = run_step(sys, vec![time]).expect_err("should fail");
1207 assert!(err.message().contains("more than 1000000 samples"));
1208 }
1209
1210 #[test]
1211 fn rejects_non_tf_input() {
1212 let err = run_step(Value::Num(1.0), Vec::new()).expect_err("expected error");
1213 assert!(err.message().contains("expected a tf object"));
1214 assert_eq!(err.identifier(), STEP_ERROR_INVALID_MODEL.identifier);
1215 }
1216}