Merge branch 'master' of github.com:SciML/SciMLBase.jl into ap/bvp

.github/workflows/Downstream.yml

@@ -43,6 +43,11 @@ jobs:
           - {user: SciML, repo: StochasticDelayDiffEq.jl, group: All}
           - {user: SciML, repo: SimpleNonlinearSolve.jl, group: All}
           - {user: SciML, repo: SimpleDiffEq.jl, group: All}
+          - {user: SciML, repo: SciMLSensitivity.jl, group: Core1}
+          - {user: SciML, repo: SciMLSensitivity.jl, group: Core2}
+          - {user: SciML, repo: SciMLSensitivity.jl, group: Core3}
+          - {user: SciML, repo: SciMLSensitivity.jl, group: Core4}
+          - {user: SciML, repo: SciMLSensitivity.jl, group: Core5}
 
     steps:
       - uses: actions/checkout@v4

Project.toml

@@ -1,7 +1,7 @@
 name = "SciMLBase"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 authors = ["Chris Rackauckas <[email protected]> and contributors"]
-version = "1.98.0"  # Remember to make it 2.0.0
+version = "1.98.2"  # Remember to make it 2.0.0
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -12,6 +12,7 @@ ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
 IteratorInterfaceExtensions = "82899510-4779-5014-852e-03e436cf321d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -45,6 +46,7 @@ CommonSolve = "0.2.4"
 ConstructionBase = "1"
 DocStringExtensions = "0.8, 0.9"
 EnumX = "1"
+FillArrays = "1.6"
 FunctionWrappersWrappers = "0.1.3"
 IteratorInterfaceExtensions = "^0.1, ^1"
 PrecompileTools = "1"

src/SciMLBase.jl

@@ -23,6 +23,7 @@ import TruncatedStacktraces
 import ADTypes: AbstractADType
 import ChainRulesCore
 import ZygoteRules: @adjoint
+import FillArrays
 
 using Reexport
 using SciMLOperators
@@ -760,9 +761,7 @@ export solve, solve!, init, discretize, symbolic_discretize
 
 export LinearProblem,
     NonlinearProblem, IntervalNonlinearProblem,
-    IntegralProblem, OptimizationProblem
-
-export IntegralProblem
+    IntegralProblem, SampledIntegralProblem, OptimizationProblem
 
 export DiscreteProblem, ImplicitDiscreteProblem
 export SteadyStateProblem, SteadyStateSolution

src/problems/basic_problems.jl

@@ -334,10 +334,12 @@ which are `Number`s or `AbstractVector`s with the same geometry as `u`.
 
 ### Constructors
 
+```
 IntegralProblem{iip}(f,lb,ub,p=NullParameters();
                   nout=1, batch = 0, kwargs...)
+```
 
-- f: the integrand, `y = f(u,p)` for out-of-place or `f(y,u,p)` for in-place.
+- f: the integrand, callable function `y = f(u,p)` for out-of-place or `f(y,u,p)` for in-place.
 - lb: Either a number or vector of lower bounds.
 - ub: Either a number or vector of upper bounds.
 - p: The parameters associated with the problem.
@@ -375,14 +377,16 @@ struct IntegralProblem{isinplace, P, F, B, K} <: AbstractIntegralProblem{isinpla
         batch = 0, kwargs...) where {iip}
         @assert typeof(lb)==typeof(ub) "Type of lower and upper bound must match"
         warn_paramtype(p)
-        new{iip, typeof(p), typeof(f), typeof(lb), typeof(kwargs)}(f, lb, ub, nout, p,
+        new{iip, typeof(p), typeof(f), typeof(lb), typeof(kwargs)}(f,
+            lb, ub, nout, p,
             batch, kwargs)
     end
 end
 
 TruncatedStacktraces.@truncate_stacktrace IntegralProblem 1 4
 
-function IntegralProblem(f, lb, ub, args...; kwargs...)
+function IntegralProblem(f, lb, ub, args...;
+    kwargs...)
     IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...)
 end
 
@@ -391,6 +395,56 @@ struct QuadratureProblem end
 
 @doc doc"""
 
+Defines a integral problem over pre-sampled data.
+Documentation Page: https://docs.sciml.ai/Integrals/stable/
+
+## Mathematical Specification of a data Integral Problem
+
+Sampled integral problems are defined as:
+
+```math
+\sum_i w_i y_i
+```
+where `y_i` are sampled values of the integrand, and `w_i` are weights 
+assigned by a quadrature rule, which depend on sampling points `x`. 
+
+## Problem Type
+
+### Constructors
+
+```
+SampledIntegralProblem(y::AbstractArray, x::AbstractVector; dim=ndims(y), kwargs...)
+```
+- y: The sampled integrand, must be a subtype of `AbstractArray`. 
+  It is assumed that the values of `y` along dimension `dim` 
+  correspond to the integrand evaluated at sampling points `x`
+- x: Sampling points, must be a subtype of `AbstractVector`.   
+- dim: Dimension along which to integrate. Defaults to the last dimension of `y`.
+- kwargs: Keyword arguments copied to the solvers.
+
+### Fields
+
+The fields match the names of the constructor arguments.
+"""
+struct SampledIntegralProblem{Y, X, K} <: AbstractIntegralProblem{false}
+    y::Y
+    x::X
+    dim::Int
+    kwargs::K
+    @add_kwonly function SampledIntegralProblem(y::AbstractArray, x::AbstractVector;
+        dim = ndims(y),
+        kwargs...)
+        @assert dim <= ndims(y) "The integration dimension `dim` is larger than the number of dimensions of the integrand `y`"
+        @assert length(x)==size(y, dim) "The integrand `y` must have the same length as the sampling points `x` along the integrated dimension."
+        @assert axes(x, 1)==axes(y, dim) "The integrand `y` must obey the same indexing as the sampling points `x` along the integrated dimension."
+        new{typeof(y), typeof(x), typeof(kwargs)}(y, x, dim, kwargs)
+    end
+end
+
+TruncatedStacktraces.@truncate_stacktrace SampledIntegralProblem 1 4
+
+@doc doc"""
+
 Defines an optimization problem.
 Documentation Page: https://docs.sciml.ai/Optimization/stable/API/optimization_problem/
 

src/scimlfunctions.jl

@@ -440,28 +440,28 @@ and exponential integrators.
 
 ```julia
 SplitFunction{iip,specialize}(f1,f2;
-                             mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I,
-                             analytic = __has_analytic(f) ? f.analytic : nothing,
-                             tgrad= __has_tgrad(f) ? f.tgrad : nothing,
-                             jac = __has_jac(f) ? f.jac : nothing,
-                             jvp = __has_jvp(f) ? f.jvp : nothing,
-                             vjp = __has_vjp(f) ? f.vjp : nothing,
-                             jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing,
-                             sparsity = __has_sparsity(f) ? f.sparsity : jac_prototype,
-                             paramjac = __has_paramjac(f) ? f.paramjac : nothing,
-                             syms = __has_syms(f) ? f.syms : nothing,
-                             indepsym= __has_indepsym(f) ? f.indepsym : nothing,
-                             paramsyms = __has_paramsyms(f) ? f.paramsyms : nothing,
-                             colorvec = __has_colorvec(f) ? f.colorvec : nothing,
-                             sys = __has_sys(f) ? f.sys : nothing)
+                             mass_matrix = __has_mass_matrix(f1) ? f1.mass_matrix : I,
+                             analytic = __has_analytic(f1) ? f1.analytic : nothing,
+                             tgrad= __has_tgrad(f1) ? f1.tgrad : nothing,
+                             jac = __has_jac(f1) ? f1.jac : nothing,
+                             jvp = __has_jvp(f1) ? f1.jvp : nothing,
+                             vjp = __has_vjp(f1) ? f1.vjp : nothing,
+                             jac_prototype = __has_jac_prototype(f1) ? f1.jac_prototype : nothing,
+                             sparsity = __has_sparsity(f1) ? f1.sparsity : jac_prototype,
+                             paramjac = __has_paramjac(f1) ? f1.paramjac : nothing,
+                             syms = __has_syms(f1) ? f1.syms : nothing,
+                             indepsym= __has_indepsym(f1) ? f1.indepsym : nothing,
+                             paramsyms = __has_paramsyms(f1) ? f1.paramsyms : nothing,
+                             colorvec = __has_colorvec(f1) ? f1.colorvec : nothing,
+                             sys = __has_sys(f1) ? f1.sys : nothing)
 ```
 
 Note that only the functions `f_i` themselves are required. These functions should
 be given as `f_i!(du,u,p,t)` or `du = f_i(u,p,t)`. See the section on `iip`
 for more details on in-place vs out-of-place handling.
 
 All of the remaining functions are optional for improving or accelerating
-the usage of `f`. These include:
+the usage of the `SplitFunction`. These include:
 
 - `mass_matrix`: the mass matrix `M` represented in the ODE function. Can be used
   to determine that the equation is actually a differential-algebraic equation (DAE)

src/solutions/zygote.jl

@@ -1,6 +1,6 @@
 @adjoint function getindex(VA::ODESolution, i::Int)
     function ODESolution_getindex_pullback(Δ)
-        Δ′ = [[i == k ? Δ[j] : zero(x[1]) for k in 1:length(x)]
+        Δ′ = [(i == j ? Δ : FillArrays.Fill(zero(eltype(x)), size(x)))
               for (x, j) in zip(VA.u, 1:length(VA))]
         (Δ′, nothing)
     end

test/function_building_error_messages.jl

@@ -453,7 +453,17 @@ NonlinearFunction(nfoop, vjp = nvjp)
 intf(u) = 1.0
 @test_throws SciMLBase.TooFewArgumentsError IntegralProblem(intf, 0.0, 1.0)
 intf(u, p) = 1.0
+p = 2.0
+
 IntegralProblem(intf, 0.0, 1.0)
+IntegralProblem(intf, 0.0, 1.0, p)
+IntegralProblem(intf, [0.0], [1.0])
+IntegralProblem(intf, [0.0], [1.0], p)
+
+x = [1.0, 2.0]
+y = rand(2, 2)
+SampledIntegralProblem(y, x)
+SampledIntegralProblem(y, x; dim=2)
 
 # Optimization