1use num_complex::Complex64;
4use runmat_builtins::{
5 BuiltinCompletionPolicy, BuiltinDescriptor, BuiltinErrorDescriptor, BuiltinOutputMode,
6 BuiltinParamArity, BuiltinParamDescriptor, BuiltinParamType, BuiltinSignatureDescriptor,
7 Tensor, Value,
8};
9use runmat_macros::runtime_builtin;
10
11use crate::builtins::common::spec::{
12 BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
13 ReductionNaN, ResidencyPolicy, ShapeRequirements,
14};
15use crate::builtins::control::tf_model::{
16 control_error, output_complex_column, polynomial_roots, ss_poles_from_object, TfModel, EPS,
17 SS_CLASS, TF_CLASS,
18};
19use crate::builtins::control::type_resolvers::damp_type;
20use crate::{dispatcher, BuiltinResult};
21
22const BUILTIN_NAME: &str = "damp";
23const DAMP_PARAM_WN: BuiltinParamDescriptor = BuiltinParamDescriptor {
24 name: "wn",
25 ty: BuiltinParamType::NumericArray,
26 arity: BuiltinParamArity::Required,
27 default: None,
28 description: "Natural frequencies as an N-by-1 column vector.",
29};
30const DAMP_PARAM_ZETA: BuiltinParamDescriptor = BuiltinParamDescriptor {
31 name: "zeta",
32 ty: BuiltinParamType::NumericArray,
33 arity: BuiltinParamArity::Required,
34 default: None,
35 description: "Damping ratios as an N-by-1 column vector.",
36};
37const DAMP_PARAM_P: BuiltinParamDescriptor = BuiltinParamDescriptor {
38 name: "p",
39 ty: BuiltinParamType::Any,
40 arity: BuiltinParamArity::Required,
41 default: None,
42 description: "Model poles as an N-by-1 real or complex column vector.",
43};
44const DAMP_OUTPUT_WN: [BuiltinParamDescriptor; 1] = [DAMP_PARAM_WN];
45const DAMP_OUTPUT_WNZETA: [BuiltinParamDescriptor; 2] = [DAMP_PARAM_WN, DAMP_PARAM_ZETA];
46const DAMP_OUTPUT_WNZETAP: [BuiltinParamDescriptor; 3] =
47 [DAMP_PARAM_WN, DAMP_PARAM_ZETA, DAMP_PARAM_P];
48const DAMP_INPUTS: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
49 name: "sys",
50 ty: BuiltinParamType::Any,
51 arity: BuiltinParamArity::Required,
52 default: None,
53 description: "SISO tf model or ss state-space model.",
54}];
55const DAMP_SIGNATURES: [BuiltinSignatureDescriptor; 3] = [
56 BuiltinSignatureDescriptor {
57 label: "wn = damp(sys)",
58 inputs: &DAMP_INPUTS,
59 outputs: &DAMP_OUTPUT_WN,
60 },
61 BuiltinSignatureDescriptor {
62 label: "[wn, zeta] = damp(sys)",
63 inputs: &DAMP_INPUTS,
64 outputs: &DAMP_OUTPUT_WNZETA,
65 },
66 BuiltinSignatureDescriptor {
67 label: "[wn, zeta, p] = damp(sys)",
68 inputs: &DAMP_INPUTS,
69 outputs: &DAMP_OUTPUT_WNZETAP,
70 },
71];
72const DAMP_ERRORS: [BuiltinErrorDescriptor; 4] = [
73 BuiltinErrorDescriptor {
74 code: "RM.DAMP.INVALID_MODEL",
75 identifier: Some("RunMat:damp:InvalidModel"),
76 when: "Input system is malformed or missing model metadata.",
77 message: "damp: invalid model",
78 },
79 BuiltinErrorDescriptor {
80 code: "RM.DAMP.UNSUPPORTED_MODEL",
81 identifier: Some("RunMat:damp:UnsupportedModel"),
82 when: "Model class is unsupported.",
83 message: "damp: unsupported model",
84 },
85 BuiltinErrorDescriptor {
86 code: "RM.DAMP.OUTPUT_COUNT",
87 identifier: Some("RunMat:damp:OutputCount"),
88 when: "More than three outputs are requested.",
89 message: "damp: too many output arguments",
90 },
91 BuiltinErrorDescriptor {
92 code: "RM.DAMP.INTERNAL",
93 identifier: Some("RunMat:damp:Internal"),
94 when: "Pole extraction or output construction failed.",
95 message: "damp: internal error",
96 },
97];
98pub const DAMP_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
99 signatures: &DAMP_SIGNATURES,
100 output_mode: BuiltinOutputMode::ByRequestedOutputCount,
101 completion_policy: BuiltinCompletionPolicy::Public,
102 errors: &DAMP_ERRORS,
103};
104
105#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::control::damp")]
106pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
107 name: "damp",
108 op_kind: GpuOpKind::Custom("control-damping"),
109 supported_precisions: &[],
110 broadcast: BroadcastSemantics::None,
111 provider_hooks: &[],
112 constant_strategy: ConstantStrategy::InlineLiteral,
113 residency: ResidencyPolicy::GatherImmediately,
114 nan_mode: ReductionNaN::Include,
115 two_pass_threshold: None,
116 workgroup_size: None,
117 accepts_nan_mode: false,
118 notes: "damp computes pole-derived control metadata on the host.",
119};
120
121#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::control::damp")]
122pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
123 name: "damp",
124 shape: ShapeRequirements::Any,
125 constant_strategy: ConstantStrategy::InlineLiteral,
126 elementwise: None,
127 reduction: None,
128 emits_nan: true,
129 notes: "damp is model analysis and is not fused.",
130};
131
132#[runtime_builtin(
133 name = "damp",
134 category = "control",
135 summary = "Return natural frequencies, damping ratios, and poles for control models.",
136 keywords = "damp,damping ratio,natural frequency,poles,control system,tf,ss",
137 type_resolver(damp_type),
138 descriptor(crate::builtins::control::damp::DAMP_DESCRIPTOR),
139 builtin_path = "crate::builtins::control::damp"
140)]
141async fn damp_builtin(sys: Value) -> BuiltinResult<Value> {
142 let eval = DampingEval::from_value(sys).await?;
143 match crate::output_count::current_output_count() {
144 Some(0) => Ok(Value::OutputList(Vec::new())),
145 Some(1) => Ok(Value::OutputList(vec![eval.wn_value()?])),
146 Some(2) => Ok(Value::OutputList(vec![
147 eval.wn_value()?,
148 eval.zeta_value()?,
149 ])),
150 Some(3) => Ok(Value::OutputList(vec![
151 eval.wn_value()?,
152 eval.zeta_value()?,
153 eval.poles_value()?,
154 ])),
155 Some(_) => Err(damp_error(
156 "RunMat:damp:OutputCount",
157 "damp: at most three outputs are supported",
158 )),
159 None => eval.wn_value(),
160 }
161}
162
163#[derive(Clone, Debug)]
164struct DampingEval {
165 modes: Vec<DampingMode>,
166}
167
168#[derive(Clone, Copy, Debug)]
169struct DampingMode {
170 pole: Complex64,
171 wn: f64,
172 zeta: f64,
173}
174
175impl DampingEval {
176 async fn from_value(value: Value) -> BuiltinResult<Self> {
177 let gathered = dispatcher::gather_if_needed_async(&value).await?;
178 match gathered {
179 Value::Object(object) if object.is_class(TF_CLASS) => {
180 let model = TfModel::from_value(Value::Object(object), BUILTIN_NAME)?;
181 let poles = polynomial_roots(&model.denominator, BUILTIN_NAME)?;
182 Ok(Self::from_poles(poles, model.sample_time))
183 }
184 Value::Object(object) if object.is_class(SS_CLASS) => {
185 let (poles, sample_time) = ss_poles_from_object(&object, BUILTIN_NAME)?;
186 Ok(Self::from_poles(poles, sample_time))
187 }
188 Value::Object(object) => Err(damp_error(
189 "RunMat:damp:UnsupportedModel",
190 format!(
191 "damp: unsupported model class '{}'; supported classes are tf and ss",
192 object.class_name
193 ),
194 )),
195 other => Err(damp_error(
196 "RunMat:damp:InvalidModel",
197 format!("damp: expected a tf or ss object, got {other:?}"),
198 )),
199 }
200 }
201
202 fn from_poles(poles: Vec<Complex64>, sample_time: f64) -> Self {
203 let mut modes = poles
204 .into_iter()
205 .map(|pole| {
206 if sample_time > 0.0 && pole.norm() <= EPS {
207 return DampingMode {
208 pole,
209 wn: f64::INFINITY,
210 zeta: 1.0,
211 };
212 }
213 let equivalent = if sample_time > 0.0 {
214 pole.ln() / sample_time
215 } else {
216 pole
217 };
218 let wn = equivalent.norm();
219 let zeta = if wn <= EPS {
220 f64::NAN
221 } else {
222 -equivalent.re / wn
223 };
224 DampingMode { pole, wn, zeta }
225 })
226 .collect::<Vec<_>>();
227 modes.sort_by(compare_modes);
228 Self { modes }
229 }
230
231 fn wn_value(&self) -> BuiltinResult<Value> {
232 real_column(self.modes.iter().map(|mode| mode.wn).collect())
233 }
234
235 fn zeta_value(&self) -> BuiltinResult<Value> {
236 real_column(self.modes.iter().map(|mode| mode.zeta).collect())
237 }
238
239 fn poles_value(&self) -> BuiltinResult<Value> {
240 output_complex_column(
241 self.modes.iter().map(|mode| mode.pole).collect(),
242 BUILTIN_NAME,
243 )
244 }
245}
246
247fn compare_modes(lhs: &DampingMode, rhs: &DampingMode) -> std::cmp::Ordering {
248 lhs.wn
249 .total_cmp(&rhs.wn)
250 .then_with(|| lhs.zeta.total_cmp(&rhs.zeta))
251 .then_with(|| lhs.pole.re.total_cmp(&rhs.pole.re))
252 .then_with(|| lhs.pole.im.total_cmp(&rhs.pole.im))
253}
254
255fn real_column(data: Vec<f64>) -> BuiltinResult<Value> {
256 let rows = data.len();
257 Tensor::new(data, vec![rows, 1])
258 .map(Value::Tensor)
259 .map_err(|err| {
260 damp_error(
261 "RunMat:damp:Internal",
262 format!("damp: failed to build output tensor: {err}"),
263 )
264 })
265}
266
267fn damp_error(identifier: &'static str, message: impl Into<String>) -> crate::RuntimeError {
268 control_error(BUILTIN_NAME, identifier, message)
269}
270
271#[cfg(test)]
272mod tests {
273 use super::*;
274 use futures::executor::block_on;
275 use runmat_builtins::ComplexTensor;
276
277 fn tf(numerator: Value, denominator: Value, rest: &[Value]) -> Value {
278 let mut args = vec![numerator, denominator];
279 args.extend_from_slice(rest);
280 block_on(crate::call_builtin_async("tf", &args)).expect("tf")
281 }
282
283 fn ss(a: Value, b: Value, c: Value, d: Value, rest: &[Value]) -> Value {
284 let mut args = vec![a, b, c, d];
285 args.extend_from_slice(rest);
286 block_on(crate::call_builtin_async("ss", &args)).expect("ss")
287 }
288
289 fn damp_outputs(sys: Value, count: usize) -> Vec<Value> {
290 let _guard = crate::output_count::push_output_count(Some(count));
291 let Value::OutputList(outputs) = block_on(damp_builtin(sys)).expect("damp") else {
292 panic!("expected output list");
293 };
294 outputs
295 }
296
297 fn tensor_data(value: &Value) -> &[f64] {
298 let Value::Tensor(tensor) = value else {
299 panic!("expected tensor, got {value:?}");
300 };
301 &tensor.data
302 }
303
304 #[test]
305 fn damp_descriptor_signatures_cover_multi_output_forms() {
306 let labels: Vec<&str> = DAMP_DESCRIPTOR
307 .signatures
308 .iter()
309 .map(|sig| sig.label)
310 .collect();
311 assert!(labels.contains(&"wn = damp(sys)"));
312 assert!(labels.contains(&"[wn, zeta] = damp(sys)"));
313 assert!(labels.contains(&"[wn, zeta, p] = damp(sys)"));
314 }
315
316 #[test]
317 fn damp_tf_second_order_returns_frequency_ratio_and_poles() {
318 let sys = tf(
319 Value::Num(1.0),
320 Value::Tensor(Tensor::new(vec![100.0, 50.0, 400.0], vec![1, 3]).unwrap()),
321 &[],
322 );
323
324 let outputs = damp_outputs(sys, 3);
325 let wn = tensor_data(&outputs[0]);
326 let zeta = tensor_data(&outputs[1]);
327 assert_eq!(wn.len(), 2);
328 assert_eq!(zeta.len(), 2);
329 for value in wn {
330 assert!((*value - 2.0).abs() < 1.0e-8);
331 }
332 for value in zeta {
333 assert!((*value - 0.125).abs() < 1.0e-8);
334 }
335 match &outputs[2] {
336 Value::ComplexTensor(poles) => {
337 assert_eq!(poles.shape, vec![2, 1]);
338 assert!(poles.data.iter().all(|(re, _)| (*re + 0.25).abs() < 1.0e-8));
339 assert!(poles
340 .data
341 .iter()
342 .any(|(_, im)| (*im - 1.984313483298443).abs() < 1.0e-8));
343 assert!(poles
344 .data
345 .iter()
346 .any(|(_, im)| (*im + 1.984313483298443).abs() < 1.0e-8));
347 }
348 other => panic!("expected complex poles, got {other:?}"),
349 }
350 }
351
352 #[test]
353 fn damp_tf_first_order_returns_real_pole() {
354 let sys = tf(
355 Value::Num(1.0),
356 Value::Tensor(Tensor::new(vec![1.0, 2.0], vec![1, 2]).unwrap()),
357 &[],
358 );
359
360 let outputs = damp_outputs(sys, 3);
361 assert!((tensor_data(&outputs[0])[0] - 2.0).abs() < 1.0e-12);
362 assert!((tensor_data(&outputs[1])[0] - 1.0).abs() < 1.0e-12);
363 assert_eq!(tensor_data(&outputs[2]), &[-2.0]);
364 }
365
366 #[test]
367 fn damp_discrete_tf_uses_log_equivalent_poles() {
368 let sys = tf(
369 Value::Num(1.0),
370 Value::Tensor(Tensor::new(vec![1.0, -0.5], vec![1, 2]).unwrap()),
371 &[Value::Num(0.1)],
372 );
373
374 let outputs = damp_outputs(sys, 3);
375 assert!((tensor_data(&outputs[0])[0] - 6.931471805599453).abs() < 1.0e-10);
376 assert!((tensor_data(&outputs[1])[0] - 1.0).abs() < 1.0e-12);
377 assert!((tensor_data(&outputs[2])[0] - 0.5).abs() < 1.0e-12);
378 }
379
380 #[test]
381 fn damp_discrete_zero_pole_reports_deadbeat_damping() {
382 let report = DampingEval::from_poles(vec![Complex64::new(0.0, 0.0)], 0.1);
383 assert_eq!(report.modes.len(), 1);
384 assert!(report.modes[0].wn.is_infinite());
385 assert_eq!(report.modes[0].zeta, 1.0);
386 }
387
388 #[test]
389 fn damp_ss_uses_state_matrix_eigenvalues() {
390 let sys = ss(
391 Value::Tensor(Tensor::new(vec![0.0, -4.0, 1.0, -0.5], vec![2, 2]).unwrap()),
392 Value::Tensor(Tensor::new(vec![0.0, 1.0], vec![2, 1]).unwrap()),
393 Value::Tensor(Tensor::new(vec![1.0, 0.0], vec![1, 2]).unwrap()),
394 Value::Num(0.0),
395 &[],
396 );
397
398 let outputs = damp_outputs(sys, 3);
399 for value in tensor_data(&outputs[0]) {
400 assert!((*value - 2.0).abs() < 1.0e-8);
401 }
402 for value in tensor_data(&outputs[1]) {
403 assert!((*value - 0.125).abs() < 1.0e-8);
404 }
405 assert!(
406 matches!(&outputs[2], Value::ComplexTensor(ComplexTensor { shape, .. }) if shape == &vec![2, 1])
407 );
408 }
409
410 #[test]
411 fn damp_rejects_more_than_three_outputs() {
412 let sys = tf(
413 Value::Num(1.0),
414 Value::Tensor(Tensor::new(vec![1.0, 2.0], vec![1, 2]).unwrap()),
415 &[],
416 );
417
418 let _guard = crate::output_count::push_output_count(Some(4));
419 let err = block_on(damp_builtin(sys)).expect_err("too many outputs should fail");
420 assert_eq!(err.identifier(), Some("RunMat:damp:OutputCount"));
421 }
422}