1use runmat_builtins::{
4 BuiltinCompletionPolicy, BuiltinDescriptor, BuiltinErrorDescriptor, BuiltinExtensionDescriptor,
5 BuiltinExtensionMode, BuiltinIntegerBackendRule, BuiltinIntegerCapabilityDescriptor,
6 BuiltinIntegerComputationDomain, BuiltinIntegerInputAvailability,
7 BuiltinIntegerInputCapability, BuiltinIntegerOutputClassRule, BuiltinIntegerOverflowRule,
8 BuiltinIntegerOverloadKind, BuiltinIntegerScalarDoubleRule, BuiltinOutputMode,
9 BuiltinParamArity, BuiltinParamDescriptor, BuiltinParamType, BuiltinSignatureDescriptor,
10};
11use runmat_macros::runtime_builtin;
12use runmat_value::{LogicalArray, StructValue, Tensor, Value};
13
14use crate::builtins::common::spec::{
15 BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
16 ReductionNaN, ResidencyPolicy, ShapeRequirements,
17};
18use crate::builtins::common::tensor;
19use crate::builtins::math::optim::common::call_function;
20use crate::builtins::math::optim::type_resolvers::numerical_integral_type;
21use crate::{build_runtime_error, BuiltinResult, RuntimeError};
22
23const NAME: &str = "integral";
24const DEFAULT_ABS_TOL: f64 = 1.0e-10;
25const DEFAULT_REL_TOL: f64 = 1.0e-6;
26const DEFAULT_MAX_FUN_EVALS: usize = 10_000;
27const MAX_DEPTH: usize = 30;
28
29const INTEGRAL_OUTPUT_Q: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
30 name: "q",
31 ty: BuiltinParamType::NumericScalar,
32 arity: BuiltinParamArity::Required,
33 default: None,
34 description: "Numerical integral estimate.",
35}];
36
37const INTEGRAL_INPUTS_CORE: [BuiltinParamDescriptor; 3] = [
38 BuiltinParamDescriptor {
39 name: "fun",
40 ty: BuiltinParamType::Any,
41 arity: BuiltinParamArity::Required,
42 default: None,
43 description: "Scalar integrand callback.",
44 },
45 BuiltinParamDescriptor {
46 name: "xmin",
47 ty: BuiltinParamType::Any,
48 arity: BuiltinParamArity::Required,
49 default: None,
50 description: "Lower integration bound.",
51 },
52 BuiltinParamDescriptor {
53 name: "xmax",
54 ty: BuiltinParamType::Any,
55 arity: BuiltinParamArity::Required,
56 default: None,
57 description: "Upper integration bound.",
58 },
59];
60
61const INTEGRAL_INPUTS_OPTIONS_STRUCT: [BuiltinParamDescriptor; 4] = [
62 BuiltinParamDescriptor {
63 name: "fun",
64 ty: BuiltinParamType::Any,
65 arity: BuiltinParamArity::Required,
66 default: None,
67 description: "Scalar integrand callback.",
68 },
69 BuiltinParamDescriptor {
70 name: "xmin",
71 ty: BuiltinParamType::Any,
72 arity: BuiltinParamArity::Required,
73 default: None,
74 description: "Lower integration bound.",
75 },
76 BuiltinParamDescriptor {
77 name: "xmax",
78 ty: BuiltinParamType::Any,
79 arity: BuiltinParamArity::Required,
80 default: None,
81 description: "Upper integration bound.",
82 },
83 BuiltinParamDescriptor {
84 name: "options",
85 ty: BuiltinParamType::Any,
86 arity: BuiltinParamArity::Optional,
87 default: None,
88 description: "Options struct for AbsTol/RelTol/MaxFunEvals.",
89 },
90];
91
92const INTEGRAL_INPUTS_NAME_VALUE: [BuiltinParamDescriptor; 5] = [
93 BuiltinParamDescriptor {
94 name: "fun",
95 ty: BuiltinParamType::Any,
96 arity: BuiltinParamArity::Required,
97 default: None,
98 description: "Scalar integrand callback.",
99 },
100 BuiltinParamDescriptor {
101 name: "xmin",
102 ty: BuiltinParamType::Any,
103 arity: BuiltinParamArity::Required,
104 default: None,
105 description: "Lower integration bound.",
106 },
107 BuiltinParamDescriptor {
108 name: "xmax",
109 ty: BuiltinParamType::Any,
110 arity: BuiltinParamArity::Required,
111 default: None,
112 description: "Upper integration bound.",
113 },
114 BuiltinParamDescriptor {
115 name: "name",
116 ty: BuiltinParamType::PropertyName,
117 arity: BuiltinParamArity::Optional,
118 default: None,
119 description: "Option name.",
120 },
121 BuiltinParamDescriptor {
122 name: "value",
123 ty: BuiltinParamType::PropertyValue,
124 arity: BuiltinParamArity::Variadic,
125 default: None,
126 description: "Option value and additional name/value pairs.",
127 },
128];
129
130const INTEGRAL_SIGNATURES: [BuiltinSignatureDescriptor; 3] = [
131 BuiltinSignatureDescriptor {
132 label: "q = integral(fun, xmin, xmax)",
133 inputs: &INTEGRAL_INPUTS_CORE,
134 outputs: &INTEGRAL_OUTPUT_Q,
135 },
136 BuiltinSignatureDescriptor {
137 label: "q = integral(fun, xmin, xmax, options)",
138 inputs: &INTEGRAL_INPUTS_OPTIONS_STRUCT,
139 outputs: &INTEGRAL_OUTPUT_Q,
140 },
141 BuiltinSignatureDescriptor {
142 label: "q = integral(fun, xmin, xmax, name, value, ...)",
143 inputs: &INTEGRAL_INPUTS_NAME_VALUE,
144 outputs: &INTEGRAL_OUTPUT_Q,
145 },
146];
147
148const INTEGRAL_ERROR_INVALID_ARGUMENT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
149 code: "RM.INTEGRAL.INVALID_ARGUMENT",
150 identifier: Some("RunMat:integral:InvalidArgument"),
151 when: "Option grammar/name-value parsing is invalid.",
152 message: "integral: invalid argument",
153};
154
155const INTEGRAL_ERROR_INVALID_INPUT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
156 code: "RM.INTEGRAL.INVALID_INPUT",
157 identifier: Some("RunMat:integral:InvalidInput"),
158 when: "Bounds/integrand/adaptive solver semantics are invalid.",
159 message: "integral: invalid input",
160};
161
162const INTEGRAL_ERRORS: [BuiltinErrorDescriptor; 2] = [
163 INTEGRAL_ERROR_INVALID_ARGUMENT,
164 INTEGRAL_ERROR_INVALID_INPUT,
165];
166
167pub const INTEGRAL_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
168 signatures: &INTEGRAL_SIGNATURES,
169 output_mode: BuiltinOutputMode::Fixed,
170 completion_policy: BuiltinCompletionPolicy::Public,
171 errors: &INTEGRAL_ERRORS,
172};
173
174const INTEGRAL_INTEGER_BOUND_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
175 id: "integral-integer-bound",
176 mode: BuiltinExtensionMode::RunMatOnly,
177 description: "integral with typed-integer bounds is a RunMat extension",
178 error_identifier: Some("RunMat:compatibility:IntegralIntegerBoundExtension"),
179};
180const INTEGRAL_LOGICAL_BOUND_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
181 id: "integral-logical-bound",
182 mode: BuiltinExtensionMode::RunMatOnly,
183 description: "integral with logical bounds is a RunMat extension",
184 error_identifier: Some("RunMat:compatibility:IntegralLogicalBoundExtension"),
185};
186const INTEGRAL_INTEGER_OPTION_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
187 id: "integral-integer-option",
188 mode: BuiltinExtensionMode::RunMatOnly,
189 description: "integral with typed-integer numeric options is a RunMat extension",
190 error_identifier: Some("RunMat:compatibility:IntegralIntegerOptionExtension"),
191};
192const INTEGRAL_LOGICAL_NUMERIC_OPTION_EXTENSION: BuiltinExtensionDescriptor =
193 BuiltinExtensionDescriptor {
194 id: "integral-logical-numeric-option",
195 mode: BuiltinExtensionMode::RunMatOnly,
196 description:
197 "integral with logical tolerance or evaluation-count options is a RunMat extension",
198 error_identifier: Some("RunMat:compatibility:IntegralLogicalNumericOptionExtension"),
199 };
200pub const INTEGRAL_EXTENSIONS: [BuiltinExtensionDescriptor; 4] = [
201 INTEGRAL_INTEGER_BOUND_EXTENSION,
202 INTEGRAL_LOGICAL_BOUND_EXTENSION,
203 INTEGRAL_INTEGER_OPTION_EXTENSION,
204 INTEGRAL_LOGICAL_NUMERIC_OPTION_EXTENSION,
205];
206
207const INTEGRAL_INTEGER_BOUND_INPUTS: [BuiltinIntegerInputCapability; 1] =
208 [BuiltinIntegerInputCapability {
209 name: "xmin or xmax",
210 classes: &crate::builtins::common::integer_capability::ALL_INTEGER_CLASSES,
211 availability: BuiltinIntegerInputAvailability::RunMatOnly,
212 scalar_double: BuiltinIntegerScalarDoubleRule::NotApplicable,
213 notes: "A typed-integer bound is accepted only when exactly representable at the binary64 quadrature boundary.",
214 }];
215const INTEGRAL_INTEGER_OPTION_INPUTS: [BuiltinIntegerInputCapability; 1] =
216 [BuiltinIntegerInputCapability {
217 name: "AbsTol, RelTol, MaxFunEvals, or MaxIntervalCount value",
218 classes: &crate::builtins::common::integer_capability::ALL_INTEGER_CLASSES,
219 availability: BuiltinIntegerInputAvailability::RunMatOnly,
220 scalar_double: BuiltinIntegerScalarDoubleRule::NotApplicable,
221 notes: "Typed-integer tolerance values cross a checked binary64 boundary; count controls remain exact through usize validation.",
222 }];
223pub const INTEGRAL_INTEGER_CAPABILITIES: [BuiltinIntegerCapabilityDescriptor; 2] = [
224 BuiltinIntegerCapabilityDescriptor {
225 form: "q = integral(fun, integer_xmin, integer_xmax, ...)",
226 inputs: &INTEGRAL_INTEGER_BOUND_INPUTS,
227 computation_domain: BuiltinIntegerComputationDomain::FloatingPoint,
228 output_class: BuiltinIntegerOutputClassRule::Double,
229 overflow: BuiltinIntegerOverflowRule::Error,
230 backend: BuiltinIntegerBackendRule::HostOnly,
231 overload: BuiltinIntegerOverloadKind::ScalarOnly,
232 notes: "RunMat-only integer bounds are classified before resident access and converted only after exactness validation; MATLAB-compatible modes retain documented single/double bounds.",
233 },
234 BuiltinIntegerCapabilityDescriptor {
235 form: "q = integral(fun, xmin, xmax, integer_option)",
236 inputs: &INTEGRAL_INTEGER_OPTION_INPUTS,
237 computation_domain: BuiltinIntegerComputationDomain::FunctionSpecific,
238 output_class: BuiltinIntegerOutputClassRule::Double,
239 overflow: BuiltinIntegerOverflowRule::Error,
240 backend: BuiltinIntegerBackendRule::HostOnly,
241 overload: BuiltinIntegerOverloadKind::StructuralParameter,
242 notes: "RunMat-only typed-integer controls are independently gated from bounds and retain exact count parsing.",
243 },
244];
245
246fn integral_error_with_detail(
247 error: &'static BuiltinErrorDescriptor,
248 detail: impl AsRef<str>,
249) -> RuntimeError {
250 let detail = detail.as_ref();
251 let message = if detail.starts_with("integral:") {
252 detail.to_string()
253 } else {
254 format!("{}: {detail}", error.message)
255 };
256 let mut builder = build_runtime_error(message).with_builtin(NAME);
257 if let Some(identifier) = error.identifier {
258 builder = builder.with_identifier(identifier);
259 }
260 builder.build()
261}
262
263fn integral_map_error(
264 err: RuntimeError,
265 fallback: &'static BuiltinErrorDescriptor,
266) -> RuntimeError {
267 if err.identifier().is_some() {
268 err
269 } else {
270 integral_error_with_detail(fallback, err.message())
271 }
272}
273
274#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::math::optim::integral")]
275pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
276 name: "integral",
277 op_kind: GpuOpKind::Custom("adaptive-quadrature"),
278 supported_precisions: &[],
279 broadcast: BroadcastSemantics::None,
280 provider_hooks: &[],
281 constant_strategy: ConstantStrategy::InlineLiteral,
282 residency: ResidencyPolicy::GatherImmediately,
283 nan_mode: ReductionNaN::Include,
284 two_pass_threshold: None,
285 workgroup_size: None,
286 accepts_nan_mode: false,
287 notes: "Host adaptive quadrature solver. Callback computations may use GPU-aware builtins, but the adaptive integration loop runs on the CPU.",
288};
289
290#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::math::optim::integral")]
291pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
292 name: "integral",
293 shape: ShapeRequirements::Any,
294 constant_strategy: ConstantStrategy::InlineLiteral,
295 elementwise: None,
296 reduction: None,
297 emits_nan: false,
298 notes: "Adaptive integration repeatedly invokes user code and terminates fusion planning.",
299};
300
301#[runtime_builtin(
302 name = "integral",
303 category = "math/optim",
304 summary = "Approximate finite scalar definite integrals using adaptive quadrature.",
305 keywords = "integral,numerical integration,adaptive quadrature,quadrature,function handle",
306 accel = "sink",
307 type_resolver(numerical_integral_type),
308 descriptor(crate::builtins::math::optim::integral::INTEGRAL_DESCRIPTOR),
309 extensions(crate::builtins::math::optim::integral::INTEGRAL_EXTENSIONS),
310 integer_capabilities(crate::builtins::math::optim::integral::INTEGRAL_INTEGER_CAPABILITIES),
311 builtin_path = "crate::builtins::math::optim::integral"
312)]
313async fn integral_builtin(
314 function: Value,
315 a: Value,
316 b: Value,
317 rest: Vec<Value>,
318) -> BuiltinResult<Value> {
319 preflight_integral_inputs(&a, &b, &rest)?;
320 let mut gathered_rest = Vec::with_capacity(rest.len());
321 for value in rest {
322 gathered_rest.push(crate::dispatcher::gather_if_needed_async(&value).await?);
323 }
324 let options = IntegralOptions::parse(gathered_rest)
325 .map_err(|err| integral_map_error(err, &INTEGRAL_ERROR_INVALID_ARGUMENT))?;
326 let a = scalar_bound("lower bound", a)
327 .await
328 .map_err(|err| integral_map_error(err, &INTEGRAL_ERROR_INVALID_INPUT))?;
329 let b = scalar_bound("upper bound", b)
330 .await
331 .map_err(|err| integral_map_error(err, &INTEGRAL_ERROR_INVALID_INPUT))?;
332 if a == b {
333 return Ok(Value::Num(0.0));
334 }
335
336 let sign = if b < a { -1.0 } else { 1.0 };
337 let lo = a.min(b);
338 let hi = a.max(b);
339 let result = integrate_finite_scalar(&function, lo, hi, &options)
340 .await
341 .map_err(|err| integral_map_error(err, &INTEGRAL_ERROR_INVALID_INPUT))?;
342 Ok(Value::Num(sign * result))
343}
344
345fn preflight_integral_inputs(a: &Value, b: &Value, rest: &[Value]) -> BuiltinResult<()> {
346 use crate::builtins::common::validation::{
347 value_has_logical_class, value_has_native_integer_class,
348 };
349 for bound in [a, b] {
350 if value_has_native_integer_class(bound) {
351 crate::compatibility::ensure_builtin_extension_enabled(
352 &INTEGRAL_INTEGER_BOUND_EXTENSION,
353 NAME,
354 )?;
355 }
356 if value_has_logical_class(bound) {
357 crate::compatibility::ensure_builtin_extension_enabled(
358 &INTEGRAL_LOGICAL_BOUND_EXTENSION,
359 NAME,
360 )?;
361 }
362 }
363 for (name, value) in integral_option_values(rest) {
364 if value_has_native_integer_class(value) {
365 crate::compatibility::ensure_builtin_extension_enabled(
366 &INTEGRAL_INTEGER_OPTION_EXTENSION,
367 NAME,
368 )?;
369 }
370 if !name.eq_ignore_ascii_case("ArrayValued") && value_has_logical_class(value) {
371 crate::compatibility::ensure_builtin_extension_enabled(
372 &INTEGRAL_LOGICAL_NUMERIC_OPTION_EXTENSION,
373 NAME,
374 )?;
375 }
376 }
377 Ok(())
378}
379
380fn integral_option_values(rest: &[Value]) -> Vec<(String, &Value)> {
381 if let [Value::Struct(fields)] = rest {
382 return fields
383 .fields
384 .iter()
385 .map(|(name, value)| (name.clone(), value))
386 .collect();
387 }
388 rest.chunks_exact(2)
389 .filter_map(|pair| option_name(&pair[0]).ok().map(|name| (name, &pair[1])))
390 .collect()
391}
392
393#[derive(Clone, Copy)]
394struct IntegralOptions {
395 abs_tol: f64,
396 rel_tol: f64,
397 max_fun_evals: usize,
398}
399
400impl IntegralOptions {
401 fn parse(rest: Vec<Value>) -> BuiltinResult<Self> {
402 let mut options = Self {
403 abs_tol: DEFAULT_ABS_TOL,
404 rel_tol: DEFAULT_REL_TOL,
405 max_fun_evals: DEFAULT_MAX_FUN_EVALS,
406 };
407 if rest.is_empty() {
408 return Ok(options);
409 }
410 if rest.len() == 1 {
411 return match &rest[0] {
412 Value::Struct(fields) => {
413 options.apply_struct(fields)?;
414 Ok(options)
415 }
416 other => Err(integral_error_with_detail(
417 &INTEGRAL_ERROR_INVALID_ARGUMENT,
418 format!("expected option name/value pairs, got {other:?}"),
419 )),
420 };
421 }
422 if !rest.len().is_multiple_of(2) {
423 return Err(integral_error_with_detail(
424 &INTEGRAL_ERROR_INVALID_ARGUMENT,
425 "expected option name/value pairs",
426 ));
427 }
428 for pair in rest.chunks(2) {
429 let name = option_name(&pair[0])?;
430 options.apply_option(&name, &pair[1])?;
431 }
432 options.validate()?;
433 Ok(options)
434 }
435
436 fn apply_struct(&mut self, fields: &StructValue) -> BuiltinResult<()> {
437 for (name, value) in &fields.fields {
438 self.apply_option(name, value)?;
439 }
440 self.validate()
441 }
442
443 fn apply_option(&mut self, name: &str, value: &Value) -> BuiltinResult<()> {
444 match name.to_ascii_lowercase().as_str() {
445 "abstol" => self.abs_tol = numeric_option("AbsTol", value)?,
446 "reltol" => self.rel_tol = numeric_option("RelTol", value)?,
447 "maxfunevals" | "maxintervalcount" => {
448 let parsed = integer_option(name, value)?;
449 if parsed < 5 {
450 return Err(integral_error_with_detail(
451 &INTEGRAL_ERROR_INVALID_ARGUMENT,
452 "MaxFunEvals must be an integer scalar >= 5",
453 ));
454 }
455 self.max_fun_evals = parsed;
456 }
457 "arrayvalued" => {
458 if bool_option("ArrayValued", value)? {
459 return Err(integral_error_with_detail(
460 &INTEGRAL_ERROR_INVALID_ARGUMENT,
461 "ArrayValued true is not supported yet",
462 ));
463 }
464 }
465 other => {
466 return Err(integral_error_with_detail(
467 &INTEGRAL_ERROR_INVALID_ARGUMENT,
468 format!("unsupported option {other}"),
469 ))
470 }
471 }
472 Ok(())
473 }
474
475 fn validate(&self) -> BuiltinResult<()> {
476 if self.abs_tol < 0.0 {
477 return Err(integral_error_with_detail(
478 &INTEGRAL_ERROR_INVALID_ARGUMENT,
479 "AbsTol must be nonnegative",
480 ));
481 }
482 if self.rel_tol < 0.0 {
483 return Err(integral_error_with_detail(
484 &INTEGRAL_ERROR_INVALID_ARGUMENT,
485 "RelTol must be nonnegative",
486 ));
487 }
488 if self.abs_tol == 0.0 && self.rel_tol == 0.0 {
489 return Err(integral_error_with_detail(
490 &INTEGRAL_ERROR_INVALID_ARGUMENT,
491 "AbsTol and RelTol cannot both be zero",
492 ));
493 }
494 Ok(())
495 }
496}
497
498fn option_name(value: &Value) -> BuiltinResult<String> {
499 match value {
500 Value::String(s) => Ok(s.clone()),
501 Value::StringArray(sa) if sa.data.len() == 1 => Ok(sa.data[0].clone()),
502 Value::CharArray(chars) if chars.rows == 1 => Ok(chars.data.iter().collect()),
503 other => Err(integral_error_with_detail(
504 &INTEGRAL_ERROR_INVALID_ARGUMENT,
505 format!("option names must be strings, got {other:?}"),
506 )),
507 }
508}
509
510async fn scalar_bound(label: &str, value: Value) -> BuiltinResult<f64> {
511 let value = crate::dispatcher::gather_if_needed_async(&value).await?;
512 if !crate::builtins::common::validation::native_integer_value_is_exact_f64(&value) {
513 return Err(integral_error_with_detail(
514 &INTEGRAL_ERROR_INVALID_INPUT,
515 format!("{label} must be exactly representable as double"),
516 ));
517 }
518 let parsed = match value {
519 Value::Num(n) => n,
520 Value::Int(i) => i.to_f64(),
521 Value::Bool(b) => {
522 if b {
523 1.0
524 } else {
525 0.0
526 }
527 }
528 Value::Tensor(tensor) if tensor::is_scalar_tensor(&tensor) => scalar_tensor_value(&tensor)
529 .ok_or_else(|| {
530 integral_error_with_detail(
531 &INTEGRAL_ERROR_INVALID_INPUT,
532 format!("{label} must be a finite real scalar"),
533 )
534 })?,
535 Value::LogicalArray(LogicalArray { data, .. }) if data.len() == 1 => {
536 if data[0] != 0 {
537 1.0
538 } else {
539 0.0
540 }
541 }
542 other => {
543 return Err(integral_error_with_detail(
544 &INTEGRAL_ERROR_INVALID_INPUT,
545 format!("{label} must be a finite real scalar, got {other:?}"),
546 ))
547 }
548 };
549 if parsed.is_finite() {
550 Ok(parsed)
551 } else {
552 Err(integral_error_with_detail(
553 &INTEGRAL_ERROR_INVALID_INPUT,
554 format!("{label} must be finite"),
555 ))
556 }
557}
558
559fn numeric_option(name: &str, value: &Value) -> BuiltinResult<f64> {
560 if !crate::builtins::common::validation::native_integer_value_is_exact_f64(value) {
561 return Err(integral_error_with_detail(
562 &INTEGRAL_ERROR_INVALID_ARGUMENT,
563 format!("option {name} must be exactly representable as double"),
564 ));
565 }
566 let parsed = match value {
567 Value::Num(n) => *n,
568 Value::Int(i) => i.to_f64(),
569 Value::Bool(b) => {
570 if *b {
571 1.0
572 } else {
573 0.0
574 }
575 }
576 Value::Tensor(tensor) if tensor::is_scalar_tensor(tensor) => {
577 tensor::tensor_value_f64(tensor, 0)
578 }
579 Value::LogicalArray(LogicalArray { data, .. }) if data.len() == 1 => {
580 if data[0] != 0 {
581 1.0
582 } else {
583 0.0
584 }
585 }
586 other => {
587 return Err(integral_error_with_detail(
588 &INTEGRAL_ERROR_INVALID_ARGUMENT,
589 format!("option {name} must be numeric, got {other:?}"),
590 ))
591 }
592 };
593 if parsed.is_finite() {
594 Ok(parsed)
595 } else {
596 Err(integral_error_with_detail(
597 &INTEGRAL_ERROR_INVALID_ARGUMENT,
598 format!("option {name} must be finite"),
599 ))
600 }
601}
602
603fn integer_option(name: &str, value: &Value) -> BuiltinResult<usize> {
604 if let Some(integer) = tensor::scalar_integer_value(value) {
605 return integer.try_to_usize().ok_or_else(|| {
606 integral_error_with_detail(
607 &INTEGRAL_ERROR_INVALID_ARGUMENT,
608 format!("option {name} must be nonnegative"),
609 )
610 });
611 }
612 let parsed = numeric_option(name, value)?;
613 if parsed < 0.0 {
614 return Err(integral_error_with_detail(
615 &INTEGRAL_ERROR_INVALID_ARGUMENT,
616 format!("option {name} must be nonnegative"),
617 ));
618 }
619 if parsed.fract() != 0.0 {
620 return Err(integral_error_with_detail(
621 &INTEGRAL_ERROR_INVALID_ARGUMENT,
622 format!("option {name} must be an integer scalar"),
623 ));
624 }
625 if parsed > usize::MAX as f64 || (usize::BITS == 64 && parsed == usize::MAX as f64) {
626 return Err(integral_error_with_detail(
627 &INTEGRAL_ERROR_INVALID_ARGUMENT,
628 format!("option {name} exceeds maximum supported size"),
629 ));
630 }
631 Ok(parsed as usize)
632}
633
634fn bool_option(name: &str, value: &Value) -> BuiltinResult<bool> {
635 match value {
636 Value::Bool(flag) => Ok(*flag),
637 Value::Num(n) if *n == 0.0 || *n == 1.0 => Ok(*n != 0.0),
638 Value::Int(i) => {
639 let raw = i.to_i64();
640 if raw == 0 || raw == 1 {
641 Ok(raw != 0)
642 } else {
643 Err(integral_error_with_detail(
644 &INTEGRAL_ERROR_INVALID_ARGUMENT,
645 format!("option {name} must be logical scalar"),
646 ))
647 }
648 }
649 other => Err(integral_error_with_detail(
650 &INTEGRAL_ERROR_INVALID_ARGUMENT,
651 format!("option {name} must be logical scalar, got {other:?}"),
652 )),
653 }
654}
655
656async fn integrate_finite_scalar(
657 function: &Value,
658 a: f64,
659 b: f64,
660 options: &IntegralOptions,
661) -> BuiltinResult<f64> {
662 let fa = call_integrand(function, a).await?;
663 let m = 0.5 * (a + b);
664 let fm = call_integrand(function, m).await?;
665 let fb = call_integrand(function, b).await?;
666 let mut evals = 3usize;
667 let whole = simpson(a, b, fa, fm, fb);
668 let tol = options.abs_tol.max(options.rel_tol * whole.abs());
669 adaptive_simpson(
670 function,
671 SimpsonState {
672 a,
673 b,
674 fa,
675 fm,
676 fb,
677 whole,
678 tol,
679 depth: MAX_DEPTH,
680 },
681 &mut evals,
682 options.max_fun_evals,
683 )
684 .await
685}
686
687#[derive(Clone, Copy)]
688struct SimpsonState {
689 a: f64,
690 b: f64,
691 fa: f64,
692 fm: f64,
693 fb: f64,
694 whole: f64,
695 tol: f64,
696 depth: usize,
697}
698
699#[async_recursion::async_recursion(?Send)]
700async fn adaptive_simpson(
701 function: &Value,
702 state: SimpsonState,
703 evals: &mut usize,
704 max_fun_evals: usize,
705) -> BuiltinResult<f64> {
706 if *evals + 2 > max_fun_evals {
707 return Err(integral_error_with_detail(
708 &INTEGRAL_ERROR_INVALID_INPUT,
709 "exceeded maximum function evaluations",
710 ));
711 }
712
713 let c = 0.5 * (state.a + state.b);
714 let d = 0.5 * (state.a + c);
715 let e = 0.5 * (c + state.b);
716 let fd = call_integrand(function, d).await?;
717 let fe = call_integrand(function, e).await?;
718 *evals += 2;
719
720 let left = simpson(state.a, c, state.fa, fd, state.fm);
721 let right = simpson(c, state.b, state.fm, fe, state.fb);
722 let refined = left + right;
723 let error = refined - state.whole;
724 if error.abs() <= 15.0 * state.tol {
725 return Ok(refined + error / 15.0);
726 }
727 if state.depth == 0 {
728 return Err(integral_error_with_detail(
729 &INTEGRAL_ERROR_INVALID_INPUT,
730 "adaptive quadrature did not converge",
731 ));
732 }
733
734 let left_value = adaptive_simpson(
735 function,
736 SimpsonState {
737 a: state.a,
738 b: c,
739 fa: state.fa,
740 fm: fd,
741 fb: state.fm,
742 whole: left,
743 tol: state.tol * 0.5,
744 depth: state.depth - 1,
745 },
746 evals,
747 max_fun_evals,
748 )
749 .await?;
750 let right_value = adaptive_simpson(
751 function,
752 SimpsonState {
753 a: c,
754 b: state.b,
755 fa: state.fm,
756 fm: fe,
757 fb: state.fb,
758 whole: right,
759 tol: state.tol * 0.5,
760 depth: state.depth - 1,
761 },
762 evals,
763 max_fun_evals,
764 )
765 .await?;
766 Ok(left_value + right_value)
767}
768
769fn scalar_tensor_value(tensor: &Tensor) -> Option<f64> {
770 tensor::tensor_values_f64(tensor).into_iter().next()
771}
772
773fn simpson(a: f64, b: f64, fa: f64, fm: f64, fb: f64) -> f64 {
774 (b - a) * (fa + 4.0 * fm + fb) / 6.0
775}
776
777async fn call_integrand(function: &Value, x: f64) -> BuiltinResult<f64> {
778 let value = call_function(function, vec![Value::Num(x)]).await?;
779 let value = crate::dispatcher::gather_if_needed_async(&value).await?;
780 match value {
781 Value::Num(n) if n.is_finite() => Ok(n),
782 Value::Int(i) => Ok(i.to_f64()),
783 Value::Bool(b) => Ok(if b { 1.0 } else { 0.0 }),
784 Value::Tensor(tensor)
785 if tensor::is_scalar_tensor(&tensor)
786 && scalar_tensor_value(&tensor)
787 .map(|value| value.is_finite())
788 .unwrap_or(false) =>
789 {
790 Ok(scalar_tensor_value(&tensor).expect("finite scalar tensor value"))
791 }
792 Value::LogicalArray(logical) if logical.data.len() == 1 => {
793 Ok(if logical.data[0] != 0 { 1.0 } else { 0.0 })
794 }
795 Value::Num(_) | Value::Tensor(_) => Err(integral_error_with_detail(
796 &INTEGRAL_ERROR_INVALID_INPUT,
797 "function value must be a finite real scalar",
798 )),
799 other => Err(integral_error_with_detail(
800 &INTEGRAL_ERROR_INVALID_INPUT,
801 format!("function value must be real numeric scalar, got {other:?}"),
802 )),
803 }
804}
805
806#[cfg(test)]
807mod tests {
808 use super::*;
809 use futures::executor::block_on;
810 use runmat_value::IntegerStorage;
811
812 const INTEGRAL_HELPER_OUTPUT: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
813 name: "fx",
814 ty: BuiltinParamType::NumericScalar,
815 arity: BuiltinParamArity::Required,
816 default: None,
817 description: "Integrand scalar value.",
818 }];
819
820 const INTEGRAL_HELPER_INPUTS: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
821 name: "x",
822 ty: BuiltinParamType::NumericScalar,
823 arity: BuiltinParamArity::Required,
824 default: None,
825 description: "Integrand sample location.",
826 }];
827
828 const INTEGRAL_HELPER_SIGNATURES: [BuiltinSignatureDescriptor; 1] =
829 [BuiltinSignatureDescriptor {
830 label: "fx = __integral_helper(x)",
831 inputs: &INTEGRAL_HELPER_INPUTS,
832 outputs: &INTEGRAL_HELPER_OUTPUT,
833 }];
834
835 const INTEGRAL_HELPER_ERRORS: [BuiltinErrorDescriptor; 0] = [];
836
837 pub const INTEGRAL_TEST_HELPER_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
838 signatures: &INTEGRAL_HELPER_SIGNATURES,
839 output_mode: BuiltinOutputMode::Fixed,
840 completion_policy: BuiltinCompletionPolicy::HiddenInternal,
841 errors: &INTEGRAL_HELPER_ERRORS,
842 };
843
844 #[runtime_builtin(
845 name = "__integral_square",
846 type_resolver(crate::builtins::math::optim::type_resolvers::numerical_integral_type),
847 descriptor(crate::builtins::math::optim::integral::tests::INTEGRAL_TEST_HELPER_DESCRIPTOR),
848 builtin_path = "crate::builtins::math::optim::integral::tests"
849 )]
850 async fn square_helper(x: Value) -> crate::BuiltinResult<Value> {
851 let x = scalar_bound("x", x).await?;
852 Ok(Value::Num(x * x))
853 }
854
855 #[runtime_builtin(
856 name = "__integral_vector",
857 type_resolver(crate::builtins::math::optim::type_resolvers::numerical_integral_type),
858 descriptor(crate::builtins::math::optim::integral::tests::INTEGRAL_TEST_HELPER_DESCRIPTOR),
859 builtin_path = "crate::builtins::math::optim::integral::tests"
860 )]
861 async fn vector_helper(_x: Value) -> crate::BuiltinResult<Value> {
862 Ok(Value::Tensor(
863 Tensor::new(vec![1.0, 2.0], vec![1, 2]).unwrap(),
864 ))
865 }
866
867 #[runtime_builtin(
868 name = "__integral_nan",
869 type_resolver(crate::builtins::math::optim::type_resolvers::numerical_integral_type),
870 descriptor(crate::builtins::math::optim::integral::tests::INTEGRAL_TEST_HELPER_DESCRIPTOR),
871 builtin_path = "crate::builtins::math::optim::integral::tests"
872 )]
873 async fn nan_helper(_x: Value) -> crate::BuiltinResult<Value> {
874 Ok(Value::Num(f64::NAN))
875 }
876
877 #[runtime_builtin(
878 name = "__integral_integer_tensor",
879 type_resolver(crate::builtins::math::optim::type_resolvers::numerical_integral_type),
880 descriptor(crate::builtins::math::optim::integral::tests::INTEGRAL_TEST_HELPER_DESCRIPTOR),
881 builtin_path = "crate::builtins::math::optim::integral::tests"
882 )]
883 async fn integer_tensor_helper(_x: Value) -> crate::BuiltinResult<Value> {
884 let tensor = Tensor::new_integer(IntegerStorage::I16(vec![3]), vec![1, 1]).expect("tensor");
885 Ok(Value::Tensor(tensor))
886 }
887
888 fn run(function: Value, a: f64, b: f64) -> crate::BuiltinResult<Value> {
889 block_on(integral_builtin(
890 function,
891 Value::Num(a),
892 Value::Num(b),
893 Vec::new(),
894 ))
895 }
896
897 #[test]
898 fn integral_test_helper_descriptor_is_attached_shape() {
899 assert_eq!(
900 INTEGRAL_TEST_HELPER_DESCRIPTOR.signatures[0].label,
901 "fx = __integral_helper(x)"
902 );
903 }
904
905 #[test]
906 fn integral_bounds_read_typed_integer_storage_exactly() {
907 let bounds = || {
908 (
909 Tensor::new_integer(IntegerStorage::I16(vec![1]), vec![1, 1]).expect("lower"),
910 Tensor::new_integer(IntegerStorage::U16(vec![3]), vec![1, 1]).expect("upper"),
911 )
912 };
913 let (lower, upper) = bounds();
914 let error = block_on(integral_builtin(
915 Value::FunctionHandle("__integral_square".to_string()),
916 Value::Tensor(lower),
917 Value::Tensor(upper),
918 Vec::new(),
919 ))
920 .expect_err("compatible mode rejects integer bounds");
921 assert_eq!(
922 error.identifier(),
923 INTEGRAL_INTEGER_BOUND_EXTENSION.error_identifier
924 );
925
926 let _guard = crate::compatibility::push_runmat_extensions_enabled(true);
927 let (lower, upper) = bounds();
928 let result = block_on(integral_builtin(
929 Value::FunctionHandle("__integral_square".to_string()),
930 Value::Tensor(lower),
931 Value::Tensor(upper),
932 Vec::new(),
933 ))
934 .expect("integral");
935
936 match result {
937 Value::Num(value) => assert!((value - (26.0 / 3.0)).abs() < 1.0e-6),
938 other => panic!("expected numeric result, got {other:?}"),
939 }
940 }
941
942 #[test]
943 fn integral_integrand_reads_typed_integer_scalar_storage_exactly() {
944 let result = run(
945 Value::FunctionHandle("__integral_integer_tensor".to_string()),
946 0.0,
947 2.0,
948 )
949 .expect("integral");
950
951 match result {
952 Value::Num(value) => assert!((value - 6.0).abs() < 1.0e-9),
953 other => panic!("expected numeric result, got {other:?}"),
954 }
955 }
956
957 #[test]
958 fn integrates_named_sine_function() {
959 let result = run(
960 Value::FunctionHandle("sin".into()),
961 0.0,
962 std::f64::consts::PI,
963 )
964 .expect("integral");
965 match result {
966 Value::Num(value) => assert!((value - 2.0).abs() < 1.0e-7),
967 other => panic!("unexpected value {other:?}"),
968 }
969 }
970
971 #[test]
972 fn integrates_polynomial_helper() {
973 let result =
974 run(Value::FunctionHandle("__integral_square".into()), 0.0, 1.0).expect("integral");
975 match result {
976 Value::Num(value) => assert!((value - (1.0 / 3.0)).abs() < 1.0e-9),
977 other => panic!("unexpected value {other:?}"),
978 }
979 }
980
981 #[test]
982 fn reversed_bounds_negate_result() {
983 let result = run(
984 Value::FunctionHandle("sin".into()),
985 std::f64::consts::PI,
986 0.0,
987 )
988 .expect("integral");
989 match result {
990 Value::Num(value) => assert!((value + 2.0).abs() < 1.0e-7),
991 other => panic!("unexpected value {other:?}"),
992 }
993 }
994
995 #[test]
996 fn zero_width_interval_returns_zero_without_callback() {
997 let result =
998 run(Value::FunctionHandle("__integral_nan".into()), 1.0, 1.0).expect("integral");
999 assert!(matches!(result, Value::Num(0.0)));
1000 }
1001
1002 #[test]
1003 fn rejects_vector_valued_integrand_for_initial_scope() {
1004 let err = run(Value::FunctionHandle("__integral_vector".into()), 0.0, 1.0).unwrap_err();
1005 assert!(err.message().contains("finite real scalar"));
1006 }
1007
1008 #[test]
1009 fn rejects_nonfinite_integrand_values() {
1010 let err = run(Value::FunctionHandle("__integral_nan".into()), 0.0, 1.0).unwrap_err();
1011 assert!(err.message().contains("finite real scalar"));
1012 }
1013
1014 #[test]
1015 fn accepts_tolerance_name_value_options() {
1016 let result = block_on(integral_builtin(
1017 Value::FunctionHandle("sin".into()),
1018 Value::Num(0.0),
1019 Value::Num(std::f64::consts::PI),
1020 vec![
1021 Value::from("AbsTol"),
1022 Value::Num(1.0e-12),
1023 Value::from("RelTol"),
1024 Value::Num(1.0e-8),
1025 ],
1026 ))
1027 .expect("integral");
1028 match result {
1029 Value::Num(value) => assert!((value - 2.0).abs() < 1.0e-8),
1030 other => panic!("unexpected value {other:?}"),
1031 }
1032 }
1033
1034 #[test]
1035 fn rejects_too_small_max_fun_evals() {
1036 let err = block_on(integral_builtin(
1037 Value::FunctionHandle("sin".into()),
1038 Value::Num(0.0),
1039 Value::Num(1.0),
1040 vec![Value::from("MaxFunEvals"), Value::Num(4.0)],
1041 ))
1042 .unwrap_err();
1043 assert!(err.message().contains("integer scalar >= 5"));
1044 }
1045
1046 #[test]
1047 fn rejects_fractional_max_fun_evals() {
1048 let err = block_on(integral_builtin(
1049 Value::FunctionHandle("sin".into()),
1050 Value::Num(0.0),
1051 Value::Num(1.0),
1052 vec![Value::from("MaxFunEvals"), Value::Num(5.5)],
1053 ))
1054 .unwrap_err();
1055 assert!(err.message().contains("integer scalar"));
1056 }
1057
1058 #[test]
1059 fn max_fun_evals_option_reads_typed_integer_storage_exactly() {
1060 let option = || {
1061 Value::Tensor(
1062 Tensor::new_integer(IntegerStorage::U16(vec![50]), vec![1, 1])
1063 .expect("MaxFunEvals"),
1064 )
1065 };
1066
1067 let error = block_on(integral_builtin(
1068 Value::FunctionHandle("sin".into()),
1069 Value::Num(0.0),
1070 Value::Num(1.0),
1071 vec![Value::from("MaxFunEvals"), option()],
1072 ))
1073 .expect_err("compatible mode rejects typed integer options");
1074 assert_eq!(
1075 error.identifier(),
1076 INTEGRAL_INTEGER_OPTION_EXTENSION.error_identifier
1077 );
1078 let _guard = crate::compatibility::push_runmat_extensions_enabled(true);
1079 let result = block_on(integral_builtin(
1080 Value::FunctionHandle("sin".into()),
1081 Value::Num(0.0),
1082 Value::Num(1.0),
1083 vec![Value::from("MaxFunEvals"), option()],
1084 ))
1085 .expect("RunMat integer option");
1086 assert!(matches!(result, Value::Num(_)));
1087 }
1088
1089 #[test]
1090 fn resident_integer_tolerance_option_gathers_after_compatibility_admission() {
1091 use crate::builtins::common::test_support;
1092 use runmat_accelerate_api::{HostIntegerDataView, HostIntegerTensorView};
1093
1094 let _guard = crate::compatibility::push_runmat_extensions_enabled(true);
1095 test_support::with_test_provider(|provider| {
1096 let value = [1_u64];
1097 let handle = provider
1098 .upload_integer(&HostIntegerTensorView {
1099 data: HostIntegerDataView::U64(&value),
1100 shape: &[1, 1],
1101 })
1102 .expect("upload resident tolerance");
1103 let result = block_on(integral_builtin(
1104 Value::FunctionHandle("sin".into()),
1105 Value::Num(0.0),
1106 Value::Num(1.0),
1107 vec![Value::from("AbsTol"), Value::GpuTensor(handle.clone())],
1108 ))
1109 .expect("resident integer tolerance");
1110 assert!(matches!(result, Value::Num(_)));
1111 provider.free(&handle).expect("free resident tolerance");
1112 runmat_accelerate_api::clear_handle_metadata(&handle);
1113 });
1114 }
1115
1116 #[test]
1117 fn max_fun_evals_option_rejects_negative_typed_integer_storage_exactly() {
1118 let _guard = crate::compatibility::push_runmat_extensions_enabled(true);
1119 let max_fun_evals =
1120 Tensor::new_integer(IntegerStorage::I16(vec![-1]), vec![1, 1]).expect("MaxFunEvals");
1121
1122 let err = block_on(integral_builtin(
1123 Value::FunctionHandle("sin".into()),
1124 Value::Num(0.0),
1125 Value::Num(1.0),
1126 vec![Value::from("MaxFunEvals"), Value::Tensor(max_fun_evals)],
1127 ))
1128 .unwrap_err();
1129 assert!(err.message().contains("MaxFunEvals"));
1130 }
1131
1132 #[test]
1133 fn max_fun_evals_option_rejects_unrepresentable_double_boundary() {
1134 let boundary = if usize::BITS == 64 {
1135 usize::MAX as f64
1136 } else {
1137 (usize::MAX as f64) + 1.0
1138 };
1139 let err = block_on(integral_builtin(
1140 Value::FunctionHandle("sin".into()),
1141 Value::Num(0.0),
1142 Value::Num(1.0),
1143 vec![Value::from("MaxFunEvals"), Value::Num(boundary)],
1144 ))
1145 .unwrap_err();
1146 assert!(err.message().contains("MaxFunEvals"));
1147 }
1148
1149 #[test]
1150 fn integral_descriptor_signatures_cover_core_forms() {
1151 let labels: Vec<&str> = INTEGRAL_DESCRIPTOR
1152 .signatures
1153 .iter()
1154 .map(|signature| signature.label)
1155 .collect();
1156 assert_eq!(
1157 labels,
1158 vec![
1159 "q = integral(fun, xmin, xmax)",
1160 "q = integral(fun, xmin, xmax, options)",
1161 "q = integral(fun, xmin, xmax, name, value, ...)",
1162 ]
1163 );
1164
1165 let codes: Vec<&str> = INTEGRAL_DESCRIPTOR
1166 .errors
1167 .iter()
1168 .map(|error| error.code)
1169 .collect();
1170 assert_eq!(
1171 codes,
1172 vec!["RM.INTEGRAL.INVALID_ARGUMENT", "RM.INTEGRAL.INVALID_INPUT"]
1173 );
1174 }
1175
1176 #[test]
1177 fn integral_bad_name_value_pairs_use_stable_identifier() {
1178 let err = block_on(integral_builtin(
1179 Value::FunctionHandle("sin".into()),
1180 Value::Num(0.0),
1181 Value::Num(1.0),
1182 vec![Value::from("AbsTol")],
1183 ))
1184 .unwrap_err();
1185 assert_eq!(err.identifier(), Some("RunMat:integral:InvalidArgument"));
1186 }
1187}