1use runmat_builtins::{
4 BuiltinCompletionPolicy, BuiltinDescriptor, BuiltinErrorDescriptor, BuiltinOutputMode,
5 BuiltinParamArity, BuiltinParamDescriptor, BuiltinParamType, BuiltinSignatureDescriptor,
6 LogicalArray, 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::math::optim::common::call_function;
15use crate::builtins::math::optim::type_resolvers::numerical_integral_type;
16use crate::{build_runtime_error, BuiltinResult, RuntimeError};
17
18const NAME: &str = "integral";
19const DEFAULT_ABS_TOL: f64 = 1.0e-10;
20const DEFAULT_REL_TOL: f64 = 1.0e-6;
21const DEFAULT_MAX_FUN_EVALS: usize = 10_000;
22const MAX_DEPTH: usize = 30;
23
24const INTEGRAL_OUTPUT_Q: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
25 name: "q",
26 ty: BuiltinParamType::NumericScalar,
27 arity: BuiltinParamArity::Required,
28 default: None,
29 description: "Numerical integral estimate.",
30}];
31
32const INTEGRAL_INPUTS_CORE: [BuiltinParamDescriptor; 3] = [
33 BuiltinParamDescriptor {
34 name: "fun",
35 ty: BuiltinParamType::Any,
36 arity: BuiltinParamArity::Required,
37 default: None,
38 description: "Scalar integrand callback.",
39 },
40 BuiltinParamDescriptor {
41 name: "xmin",
42 ty: BuiltinParamType::Any,
43 arity: BuiltinParamArity::Required,
44 default: None,
45 description: "Lower integration bound.",
46 },
47 BuiltinParamDescriptor {
48 name: "xmax",
49 ty: BuiltinParamType::Any,
50 arity: BuiltinParamArity::Required,
51 default: None,
52 description: "Upper integration bound.",
53 },
54];
55
56const INTEGRAL_INPUTS_OPTIONS_STRUCT: [BuiltinParamDescriptor; 4] = [
57 BuiltinParamDescriptor {
58 name: "fun",
59 ty: BuiltinParamType::Any,
60 arity: BuiltinParamArity::Required,
61 default: None,
62 description: "Scalar integrand callback.",
63 },
64 BuiltinParamDescriptor {
65 name: "xmin",
66 ty: BuiltinParamType::Any,
67 arity: BuiltinParamArity::Required,
68 default: None,
69 description: "Lower integration bound.",
70 },
71 BuiltinParamDescriptor {
72 name: "xmax",
73 ty: BuiltinParamType::Any,
74 arity: BuiltinParamArity::Required,
75 default: None,
76 description: "Upper integration bound.",
77 },
78 BuiltinParamDescriptor {
79 name: "options",
80 ty: BuiltinParamType::Any,
81 arity: BuiltinParamArity::Optional,
82 default: None,
83 description: "Options struct for AbsTol/RelTol/MaxFunEvals.",
84 },
85];
86
87const INTEGRAL_INPUTS_NAME_VALUE: [BuiltinParamDescriptor; 5] = [
88 BuiltinParamDescriptor {
89 name: "fun",
90 ty: BuiltinParamType::Any,
91 arity: BuiltinParamArity::Required,
92 default: None,
93 description: "Scalar integrand callback.",
94 },
95 BuiltinParamDescriptor {
96 name: "xmin",
97 ty: BuiltinParamType::Any,
98 arity: BuiltinParamArity::Required,
99 default: None,
100 description: "Lower integration bound.",
101 },
102 BuiltinParamDescriptor {
103 name: "xmax",
104 ty: BuiltinParamType::Any,
105 arity: BuiltinParamArity::Required,
106 default: None,
107 description: "Upper integration bound.",
108 },
109 BuiltinParamDescriptor {
110 name: "name",
111 ty: BuiltinParamType::PropertyName,
112 arity: BuiltinParamArity::Optional,
113 default: None,
114 description: "Option name.",
115 },
116 BuiltinParamDescriptor {
117 name: "value",
118 ty: BuiltinParamType::PropertyValue,
119 arity: BuiltinParamArity::Variadic,
120 default: None,
121 description: "Option value and additional name/value pairs.",
122 },
123];
124
125const INTEGRAL_SIGNATURES: [BuiltinSignatureDescriptor; 3] = [
126 BuiltinSignatureDescriptor {
127 label: "q = integral(fun, xmin, xmax)",
128 inputs: &INTEGRAL_INPUTS_CORE,
129 outputs: &INTEGRAL_OUTPUT_Q,
130 },
131 BuiltinSignatureDescriptor {
132 label: "q = integral(fun, xmin, xmax, options)",
133 inputs: &INTEGRAL_INPUTS_OPTIONS_STRUCT,
134 outputs: &INTEGRAL_OUTPUT_Q,
135 },
136 BuiltinSignatureDescriptor {
137 label: "q = integral(fun, xmin, xmax, name, value, ...)",
138 inputs: &INTEGRAL_INPUTS_NAME_VALUE,
139 outputs: &INTEGRAL_OUTPUT_Q,
140 },
141];
142
143const INTEGRAL_ERROR_INVALID_ARGUMENT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
144 code: "RM.INTEGRAL.INVALID_ARGUMENT",
145 identifier: Some("RunMat:integral:InvalidArgument"),
146 when: "Option grammar/name-value parsing is invalid.",
147 message: "integral: invalid argument",
148};
149
150const INTEGRAL_ERROR_INVALID_INPUT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
151 code: "RM.INTEGRAL.INVALID_INPUT",
152 identifier: Some("RunMat:integral:InvalidInput"),
153 when: "Bounds/integrand/adaptive solver semantics are invalid.",
154 message: "integral: invalid input",
155};
156
157const INTEGRAL_ERRORS: [BuiltinErrorDescriptor; 2] = [
158 INTEGRAL_ERROR_INVALID_ARGUMENT,
159 INTEGRAL_ERROR_INVALID_INPUT,
160];
161
162pub const INTEGRAL_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
163 signatures: &INTEGRAL_SIGNATURES,
164 output_mode: BuiltinOutputMode::Fixed,
165 completion_policy: BuiltinCompletionPolicy::Public,
166 errors: &INTEGRAL_ERRORS,
167};
168
169fn integral_error_with_detail(
170 error: &'static BuiltinErrorDescriptor,
171 detail: impl AsRef<str>,
172) -> RuntimeError {
173 let detail = detail.as_ref();
174 let message = if detail.starts_with("integral:") {
175 detail.to_string()
176 } else {
177 format!("{}: {detail}", error.message)
178 };
179 let mut builder = build_runtime_error(message).with_builtin(NAME);
180 if let Some(identifier) = error.identifier {
181 builder = builder.with_identifier(identifier);
182 }
183 builder.build()
184}
185
186fn integral_map_error(
187 err: RuntimeError,
188 fallback: &'static BuiltinErrorDescriptor,
189) -> RuntimeError {
190 if err.identifier().is_some() {
191 err
192 } else {
193 integral_error_with_detail(fallback, err.message())
194 }
195}
196
197#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::math::optim::integral")]
198pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
199 name: "integral",
200 op_kind: GpuOpKind::Custom("adaptive-quadrature"),
201 supported_precisions: &[],
202 broadcast: BroadcastSemantics::None,
203 provider_hooks: &[],
204 constant_strategy: ConstantStrategy::InlineLiteral,
205 residency: ResidencyPolicy::GatherImmediately,
206 nan_mode: ReductionNaN::Include,
207 two_pass_threshold: None,
208 workgroup_size: None,
209 accepts_nan_mode: false,
210 notes: "Host adaptive quadrature solver. Callback computations may use GPU-aware builtins, but the adaptive integration loop runs on the CPU.",
211};
212
213#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::math::optim::integral")]
214pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
215 name: "integral",
216 shape: ShapeRequirements::Any,
217 constant_strategy: ConstantStrategy::InlineLiteral,
218 elementwise: None,
219 reduction: None,
220 emits_nan: false,
221 notes: "Adaptive integration repeatedly invokes user code and terminates fusion planning.",
222};
223
224#[runtime_builtin(
225 name = "integral",
226 category = "math/optim",
227 summary = "Approximate finite scalar definite integrals using adaptive quadrature.",
228 keywords = "integral,numerical integration,adaptive quadrature,quadrature,function handle",
229 accel = "sink",
230 type_resolver(numerical_integral_type),
231 descriptor(crate::builtins::math::optim::integral::INTEGRAL_DESCRIPTOR),
232 builtin_path = "crate::builtins::math::optim::integral"
233)]
234async fn integral_builtin(
235 function: Value,
236 a: Value,
237 b: Value,
238 rest: Vec<Value>,
239) -> BuiltinResult<Value> {
240 let options = IntegralOptions::parse(rest)
241 .map_err(|err| integral_map_error(err, &INTEGRAL_ERROR_INVALID_ARGUMENT))?;
242 let a = scalar_bound("lower bound", a)
243 .await
244 .map_err(|err| integral_map_error(err, &INTEGRAL_ERROR_INVALID_INPUT))?;
245 let b = scalar_bound("upper bound", b)
246 .await
247 .map_err(|err| integral_map_error(err, &INTEGRAL_ERROR_INVALID_INPUT))?;
248 if a == b {
249 return Ok(Value::Num(0.0));
250 }
251
252 let sign = if b < a { -1.0 } else { 1.0 };
253 let lo = a.min(b);
254 let hi = a.max(b);
255 let result = integrate_finite_scalar(&function, lo, hi, &options)
256 .await
257 .map_err(|err| integral_map_error(err, &INTEGRAL_ERROR_INVALID_INPUT))?;
258 Ok(Value::Num(sign * result))
259}
260
261#[derive(Clone, Copy)]
262struct IntegralOptions {
263 abs_tol: f64,
264 rel_tol: f64,
265 max_fun_evals: usize,
266}
267
268impl IntegralOptions {
269 fn parse(rest: Vec<Value>) -> BuiltinResult<Self> {
270 let mut options = Self {
271 abs_tol: DEFAULT_ABS_TOL,
272 rel_tol: DEFAULT_REL_TOL,
273 max_fun_evals: DEFAULT_MAX_FUN_EVALS,
274 };
275 if rest.is_empty() {
276 return Ok(options);
277 }
278 if rest.len() == 1 {
279 return match &rest[0] {
280 Value::Struct(fields) => {
281 options.apply_struct(fields)?;
282 Ok(options)
283 }
284 other => Err(integral_error_with_detail(
285 &INTEGRAL_ERROR_INVALID_ARGUMENT,
286 format!("expected option name/value pairs, got {other:?}"),
287 )),
288 };
289 }
290 if !rest.len().is_multiple_of(2) {
291 return Err(integral_error_with_detail(
292 &INTEGRAL_ERROR_INVALID_ARGUMENT,
293 "expected option name/value pairs",
294 ));
295 }
296 for pair in rest.chunks(2) {
297 let name = option_name(&pair[0])?;
298 options.apply_option(&name, &pair[1])?;
299 }
300 options.validate()?;
301 Ok(options)
302 }
303
304 fn apply_struct(&mut self, fields: &StructValue) -> BuiltinResult<()> {
305 for (name, value) in &fields.fields {
306 self.apply_option(name, value)?;
307 }
308 self.validate()
309 }
310
311 fn apply_option(&mut self, name: &str, value: &Value) -> BuiltinResult<()> {
312 match name.to_ascii_lowercase().as_str() {
313 "abstol" => self.abs_tol = numeric_option("AbsTol", value)?,
314 "reltol" => self.rel_tol = numeric_option("RelTol", value)?,
315 "maxfunevals" | "maxintervalcount" => {
316 let parsed = integer_option(name, value)?;
317 if parsed < 5 {
318 return Err(integral_error_with_detail(
319 &INTEGRAL_ERROR_INVALID_ARGUMENT,
320 "MaxFunEvals must be an integer scalar >= 5",
321 ));
322 }
323 self.max_fun_evals = parsed;
324 }
325 "arrayvalued" => {
326 if bool_option("ArrayValued", value)? {
327 return Err(integral_error_with_detail(
328 &INTEGRAL_ERROR_INVALID_ARGUMENT,
329 "ArrayValued true is not supported yet",
330 ));
331 }
332 }
333 other => {
334 return Err(integral_error_with_detail(
335 &INTEGRAL_ERROR_INVALID_ARGUMENT,
336 format!("unsupported option {other}"),
337 ))
338 }
339 }
340 Ok(())
341 }
342
343 fn validate(&self) -> BuiltinResult<()> {
344 if self.abs_tol < 0.0 {
345 return Err(integral_error_with_detail(
346 &INTEGRAL_ERROR_INVALID_ARGUMENT,
347 "AbsTol must be nonnegative",
348 ));
349 }
350 if self.rel_tol < 0.0 {
351 return Err(integral_error_with_detail(
352 &INTEGRAL_ERROR_INVALID_ARGUMENT,
353 "RelTol must be nonnegative",
354 ));
355 }
356 if self.abs_tol == 0.0 && self.rel_tol == 0.0 {
357 return Err(integral_error_with_detail(
358 &INTEGRAL_ERROR_INVALID_ARGUMENT,
359 "AbsTol and RelTol cannot both be zero",
360 ));
361 }
362 Ok(())
363 }
364}
365
366fn option_name(value: &Value) -> BuiltinResult<String> {
367 match value {
368 Value::String(s) => Ok(s.clone()),
369 Value::StringArray(sa) if sa.data.len() == 1 => Ok(sa.data[0].clone()),
370 Value::CharArray(chars) if chars.rows == 1 => Ok(chars.data.iter().collect()),
371 other => Err(integral_error_with_detail(
372 &INTEGRAL_ERROR_INVALID_ARGUMENT,
373 format!("option names must be strings, got {other:?}"),
374 )),
375 }
376}
377
378async fn scalar_bound(label: &str, value: Value) -> BuiltinResult<f64> {
379 let value = crate::dispatcher::gather_if_needed_async(&value).await?;
380 let parsed = match value {
381 Value::Num(n) => n,
382 Value::Int(i) => i.to_f64(),
383 Value::Bool(b) => {
384 if b {
385 1.0
386 } else {
387 0.0
388 }
389 }
390 Value::Tensor(tensor) if tensor.data.len() == 1 => tensor.data[0],
391 Value::LogicalArray(LogicalArray { data, .. }) if data.len() == 1 => {
392 if data[0] != 0 {
393 1.0
394 } else {
395 0.0
396 }
397 }
398 other => {
399 return Err(integral_error_with_detail(
400 &INTEGRAL_ERROR_INVALID_INPUT,
401 format!("{label} must be a finite real scalar, got {other:?}"),
402 ))
403 }
404 };
405 if parsed.is_finite() {
406 Ok(parsed)
407 } else {
408 Err(integral_error_with_detail(
409 &INTEGRAL_ERROR_INVALID_INPUT,
410 format!("{label} must be finite"),
411 ))
412 }
413}
414
415fn numeric_option(name: &str, value: &Value) -> BuiltinResult<f64> {
416 let parsed = match value {
417 Value::Num(n) => *n,
418 Value::Int(i) => i.to_f64(),
419 Value::Bool(b) => {
420 if *b {
421 1.0
422 } else {
423 0.0
424 }
425 }
426 Value::Tensor(Tensor { data, .. }) if data.len() == 1 => data[0],
427 Value::LogicalArray(LogicalArray { data, .. }) if data.len() == 1 => {
428 if data[0] != 0 {
429 1.0
430 } else {
431 0.0
432 }
433 }
434 other => {
435 return Err(integral_error_with_detail(
436 &INTEGRAL_ERROR_INVALID_ARGUMENT,
437 format!("option {name} must be numeric, got {other:?}"),
438 ))
439 }
440 };
441 if parsed.is_finite() {
442 Ok(parsed)
443 } else {
444 Err(integral_error_with_detail(
445 &INTEGRAL_ERROR_INVALID_ARGUMENT,
446 format!("option {name} must be finite"),
447 ))
448 }
449}
450
451fn integer_option(name: &str, value: &Value) -> BuiltinResult<usize> {
452 let parsed = numeric_option(name, value)?;
453 if parsed < 0.0 {
454 return Err(integral_error_with_detail(
455 &INTEGRAL_ERROR_INVALID_ARGUMENT,
456 format!("option {name} must be nonnegative"),
457 ));
458 }
459 if parsed.fract() != 0.0 {
460 return Err(integral_error_with_detail(
461 &INTEGRAL_ERROR_INVALID_ARGUMENT,
462 format!("option {name} must be an integer scalar"),
463 ));
464 }
465 Ok(parsed as usize)
466}
467
468fn bool_option(name: &str, value: &Value) -> BuiltinResult<bool> {
469 match value {
470 Value::Bool(flag) => Ok(*flag),
471 Value::Num(n) if *n == 0.0 || *n == 1.0 => Ok(*n != 0.0),
472 Value::Int(i) => {
473 let raw = i.to_i64();
474 if raw == 0 || raw == 1 {
475 Ok(raw != 0)
476 } else {
477 Err(integral_error_with_detail(
478 &INTEGRAL_ERROR_INVALID_ARGUMENT,
479 format!("option {name} must be logical scalar"),
480 ))
481 }
482 }
483 other => Err(integral_error_with_detail(
484 &INTEGRAL_ERROR_INVALID_ARGUMENT,
485 format!("option {name} must be logical scalar, got {other:?}"),
486 )),
487 }
488}
489
490async fn integrate_finite_scalar(
491 function: &Value,
492 a: f64,
493 b: f64,
494 options: &IntegralOptions,
495) -> BuiltinResult<f64> {
496 let fa = call_integrand(function, a).await?;
497 let m = 0.5 * (a + b);
498 let fm = call_integrand(function, m).await?;
499 let fb = call_integrand(function, b).await?;
500 let mut evals = 3usize;
501 let whole = simpson(a, b, fa, fm, fb);
502 let tol = options.abs_tol.max(options.rel_tol * whole.abs());
503 adaptive_simpson(
504 function,
505 SimpsonState {
506 a,
507 b,
508 fa,
509 fm,
510 fb,
511 whole,
512 tol,
513 depth: MAX_DEPTH,
514 },
515 &mut evals,
516 options.max_fun_evals,
517 )
518 .await
519}
520
521#[derive(Clone, Copy)]
522struct SimpsonState {
523 a: f64,
524 b: f64,
525 fa: f64,
526 fm: f64,
527 fb: f64,
528 whole: f64,
529 tol: f64,
530 depth: usize,
531}
532
533#[async_recursion::async_recursion(?Send)]
534async fn adaptive_simpson(
535 function: &Value,
536 state: SimpsonState,
537 evals: &mut usize,
538 max_fun_evals: usize,
539) -> BuiltinResult<f64> {
540 if *evals + 2 > max_fun_evals {
541 return Err(integral_error_with_detail(
542 &INTEGRAL_ERROR_INVALID_INPUT,
543 "exceeded maximum function evaluations",
544 ));
545 }
546
547 let c = 0.5 * (state.a + state.b);
548 let d = 0.5 * (state.a + c);
549 let e = 0.5 * (c + state.b);
550 let fd = call_integrand(function, d).await?;
551 let fe = call_integrand(function, e).await?;
552 *evals += 2;
553
554 let left = simpson(state.a, c, state.fa, fd, state.fm);
555 let right = simpson(c, state.b, state.fm, fe, state.fb);
556 let refined = left + right;
557 let error = refined - state.whole;
558 if error.abs() <= 15.0 * state.tol {
559 return Ok(refined + error / 15.0);
560 }
561 if state.depth == 0 {
562 return Err(integral_error_with_detail(
563 &INTEGRAL_ERROR_INVALID_INPUT,
564 "adaptive quadrature did not converge",
565 ));
566 }
567
568 let left_value = adaptive_simpson(
569 function,
570 SimpsonState {
571 a: state.a,
572 b: c,
573 fa: state.fa,
574 fm: fd,
575 fb: state.fm,
576 whole: left,
577 tol: state.tol * 0.5,
578 depth: state.depth - 1,
579 },
580 evals,
581 max_fun_evals,
582 )
583 .await?;
584 let right_value = adaptive_simpson(
585 function,
586 SimpsonState {
587 a: c,
588 b: state.b,
589 fa: state.fm,
590 fm: fe,
591 fb: state.fb,
592 whole: right,
593 tol: state.tol * 0.5,
594 depth: state.depth - 1,
595 },
596 evals,
597 max_fun_evals,
598 )
599 .await?;
600 Ok(left_value + right_value)
601}
602
603fn simpson(a: f64, b: f64, fa: f64, fm: f64, fb: f64) -> f64 {
604 (b - a) * (fa + 4.0 * fm + fb) / 6.0
605}
606
607async fn call_integrand(function: &Value, x: f64) -> BuiltinResult<f64> {
608 let value = call_function(function, vec![Value::Num(x)]).await?;
609 let value = crate::dispatcher::gather_if_needed_async(&value).await?;
610 match value {
611 Value::Num(n) if n.is_finite() => Ok(n),
612 Value::Int(i) => Ok(i.to_f64()),
613 Value::Bool(b) => Ok(if b { 1.0 } else { 0.0 }),
614 Value::Tensor(tensor) if tensor.data.len() == 1 && tensor.data[0].is_finite() => {
615 Ok(tensor.data[0])
616 }
617 Value::LogicalArray(logical) if logical.data.len() == 1 => {
618 Ok(if logical.data[0] != 0 { 1.0 } else { 0.0 })
619 }
620 Value::Num(_) | Value::Tensor(_) => Err(integral_error_with_detail(
621 &INTEGRAL_ERROR_INVALID_INPUT,
622 "function value must be a finite real scalar",
623 )),
624 other => Err(integral_error_with_detail(
625 &INTEGRAL_ERROR_INVALID_INPUT,
626 format!("function value must be real numeric scalar, got {other:?}"),
627 )),
628 }
629}
630
631#[cfg(test)]
632mod tests {
633 use super::*;
634 use futures::executor::block_on;
635
636 const INTEGRAL_HELPER_OUTPUT: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
637 name: "fx",
638 ty: BuiltinParamType::NumericScalar,
639 arity: BuiltinParamArity::Required,
640 default: None,
641 description: "Integrand scalar value.",
642 }];
643
644 const INTEGRAL_HELPER_INPUTS: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
645 name: "x",
646 ty: BuiltinParamType::NumericScalar,
647 arity: BuiltinParamArity::Required,
648 default: None,
649 description: "Integrand sample location.",
650 }];
651
652 const INTEGRAL_HELPER_SIGNATURES: [BuiltinSignatureDescriptor; 1] =
653 [BuiltinSignatureDescriptor {
654 label: "fx = __integral_helper(x)",
655 inputs: &INTEGRAL_HELPER_INPUTS,
656 outputs: &INTEGRAL_HELPER_OUTPUT,
657 }];
658
659 const INTEGRAL_HELPER_ERRORS: [BuiltinErrorDescriptor; 0] = [];
660
661 pub const INTEGRAL_TEST_HELPER_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
662 signatures: &INTEGRAL_HELPER_SIGNATURES,
663 output_mode: BuiltinOutputMode::Fixed,
664 completion_policy: BuiltinCompletionPolicy::HiddenInternal,
665 errors: &INTEGRAL_HELPER_ERRORS,
666 };
667
668 #[runtime_builtin(
669 name = "__integral_square",
670 type_resolver(crate::builtins::math::optim::type_resolvers::numerical_integral_type),
671 descriptor(crate::builtins::math::optim::integral::tests::INTEGRAL_TEST_HELPER_DESCRIPTOR),
672 builtin_path = "crate::builtins::math::optim::integral::tests"
673 )]
674 async fn square_helper(x: Value) -> crate::BuiltinResult<Value> {
675 let x = scalar_bound("x", x).await?;
676 Ok(Value::Num(x * x))
677 }
678
679 #[runtime_builtin(
680 name = "__integral_vector",
681 type_resolver(crate::builtins::math::optim::type_resolvers::numerical_integral_type),
682 descriptor(crate::builtins::math::optim::integral::tests::INTEGRAL_TEST_HELPER_DESCRIPTOR),
683 builtin_path = "crate::builtins::math::optim::integral::tests"
684 )]
685 async fn vector_helper(_x: Value) -> crate::BuiltinResult<Value> {
686 Ok(Value::Tensor(
687 Tensor::new(vec![1.0, 2.0], vec![1, 2]).unwrap(),
688 ))
689 }
690
691 #[runtime_builtin(
692 name = "__integral_nan",
693 type_resolver(crate::builtins::math::optim::type_resolvers::numerical_integral_type),
694 descriptor(crate::builtins::math::optim::integral::tests::INTEGRAL_TEST_HELPER_DESCRIPTOR),
695 builtin_path = "crate::builtins::math::optim::integral::tests"
696 )]
697 async fn nan_helper(_x: Value) -> crate::BuiltinResult<Value> {
698 Ok(Value::Num(f64::NAN))
699 }
700
701 fn run(function: Value, a: f64, b: f64) -> crate::BuiltinResult<Value> {
702 block_on(integral_builtin(
703 function,
704 Value::Num(a),
705 Value::Num(b),
706 Vec::new(),
707 ))
708 }
709
710 #[test]
711 fn integral_test_helper_descriptor_is_attached_shape() {
712 assert_eq!(
713 INTEGRAL_TEST_HELPER_DESCRIPTOR.signatures[0].label,
714 "fx = __integral_helper(x)"
715 );
716 }
717
718 #[test]
719 fn integrates_named_sine_function() {
720 let result = run(
721 Value::FunctionHandle("sin".into()),
722 0.0,
723 std::f64::consts::PI,
724 )
725 .expect("integral");
726 match result {
727 Value::Num(value) => assert!((value - 2.0).abs() < 1.0e-7),
728 other => panic!("unexpected value {other:?}"),
729 }
730 }
731
732 #[test]
733 fn integrates_polynomial_helper() {
734 let result =
735 run(Value::FunctionHandle("__integral_square".into()), 0.0, 1.0).expect("integral");
736 match result {
737 Value::Num(value) => assert!((value - (1.0 / 3.0)).abs() < 1.0e-9),
738 other => panic!("unexpected value {other:?}"),
739 }
740 }
741
742 #[test]
743 fn reversed_bounds_negate_result() {
744 let result = run(
745 Value::FunctionHandle("sin".into()),
746 std::f64::consts::PI,
747 0.0,
748 )
749 .expect("integral");
750 match result {
751 Value::Num(value) => assert!((value + 2.0).abs() < 1.0e-7),
752 other => panic!("unexpected value {other:?}"),
753 }
754 }
755
756 #[test]
757 fn zero_width_interval_returns_zero_without_callback() {
758 let result =
759 run(Value::FunctionHandle("__integral_nan".into()), 1.0, 1.0).expect("integral");
760 assert!(matches!(result, Value::Num(0.0)));
761 }
762
763 #[test]
764 fn rejects_vector_valued_integrand_for_initial_scope() {
765 let err = run(Value::FunctionHandle("__integral_vector".into()), 0.0, 1.0).unwrap_err();
766 assert!(err.message().contains("finite real scalar"));
767 }
768
769 #[test]
770 fn rejects_nonfinite_integrand_values() {
771 let err = run(Value::FunctionHandle("__integral_nan".into()), 0.0, 1.0).unwrap_err();
772 assert!(err.message().contains("finite real scalar"));
773 }
774
775 #[test]
776 fn accepts_tolerance_name_value_options() {
777 let result = block_on(integral_builtin(
778 Value::FunctionHandle("sin".into()),
779 Value::Num(0.0),
780 Value::Num(std::f64::consts::PI),
781 vec![
782 Value::from("AbsTol"),
783 Value::Num(1.0e-12),
784 Value::from("RelTol"),
785 Value::Num(1.0e-8),
786 ],
787 ))
788 .expect("integral");
789 match result {
790 Value::Num(value) => assert!((value - 2.0).abs() < 1.0e-8),
791 other => panic!("unexpected value {other:?}"),
792 }
793 }
794
795 #[test]
796 fn rejects_too_small_max_fun_evals() {
797 let err = block_on(integral_builtin(
798 Value::FunctionHandle("sin".into()),
799 Value::Num(0.0),
800 Value::Num(1.0),
801 vec![Value::from("MaxFunEvals"), Value::Num(4.0)],
802 ))
803 .unwrap_err();
804 assert!(err.message().contains("integer scalar >= 5"));
805 }
806
807 #[test]
808 fn rejects_fractional_max_fun_evals() {
809 let err = block_on(integral_builtin(
810 Value::FunctionHandle("sin".into()),
811 Value::Num(0.0),
812 Value::Num(1.0),
813 vec![Value::from("MaxFunEvals"), Value::Num(5.5)],
814 ))
815 .unwrap_err();
816 assert!(err.message().contains("integer scalar"));
817 }
818
819 #[test]
820 fn integral_descriptor_signatures_cover_core_forms() {
821 let labels: Vec<&str> = INTEGRAL_DESCRIPTOR
822 .signatures
823 .iter()
824 .map(|signature| signature.label)
825 .collect();
826 assert_eq!(
827 labels,
828 vec![
829 "q = integral(fun, xmin, xmax)",
830 "q = integral(fun, xmin, xmax, options)",
831 "q = integral(fun, xmin, xmax, name, value, ...)",
832 ]
833 );
834
835 let codes: Vec<&str> = INTEGRAL_DESCRIPTOR
836 .errors
837 .iter()
838 .map(|error| error.code)
839 .collect();
840 assert_eq!(
841 codes,
842 vec!["RM.INTEGRAL.INVALID_ARGUMENT", "RM.INTEGRAL.INVALID_INPUT"]
843 );
844 }
845
846 #[test]
847 fn integral_bad_name_value_pairs_use_stable_identifier() {
848 let err = block_on(integral_builtin(
849 Value::FunctionHandle("sin".into()),
850 Value::Num(0.0),
851 Value::Num(1.0),
852 vec![Value::from("AbsTol")],
853 ))
854 .unwrap_err();
855 assert_eq!(err.identifier(), Some("RunMat:integral:InvalidArgument"));
856 }
857}