From 8ee5ca4db3a11430f73e098f33fe03919436c664 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 10 Aug 2023 16:19:14 -0400 Subject: [PATCH 01/70] TwoPointBVPFunction restructure --- Project.toml | 2 +- src/problems/bvp_problems.jl | 10 +++++++--- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index a5efdc86d..ebf74000c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SciMLBase" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" authors = ["Chris Rackauckas and contributors"] -version = "1.94.0" +version = "2.0.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" diff --git a/src/problems/bvp_problems.jl b/src/problems/bvp_problems.jl index aa5d887b4..84c9715a4 100644 --- a/src/problems/bvp_problems.jl +++ b/src/problems/bvp_problems.jl @@ -133,9 +133,13 @@ $(TYPEDEF) struct TwoPointBVPFunction{bF} bc::bF end -TwoPointBVPFunction(; bc = error("No argument bc")) = TwoPointBVPFunction(bc) -(f::TwoPointBVPFunction)(residual, ua, ub, p) = f.bc(residual, ua, ub, p) -(f::TwoPointBVPFunction)(residual, u, p) = f.bc(residual, u[1], u[end], p) +TwoPointBVPFunction(; bc = error("No argument `bc`")) = TwoPointBVPFunction(bc) +function (f::TwoPointBVPFunction)(residuala, residualb, ua, ub, p) + return f.bc(residuala, residualb, ua, ub, p) +end +function (f::TwoPointBVPFunction)(residual::Tuple, u, p) + return f(residual[1], residual[2], u[1], u[end], p) +end """ $(TYPEDEF) From 75f258e34e2bc3a484ed772b37b18460f88e9b3a Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 5 Sep 2023 16:53:31 -0400 Subject: [PATCH 02/70] This is breaking --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a14d73463..951f6b7ca 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SciMLBase" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" authors = ["Chris Rackauckas and contributors"] -version = "1.96.1" +version = "2.0.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 7bfe639c3ab701b3a816fa9cbc958d2aad6502f0 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 5 Sep 2023 17:02:49 -0400 Subject: [PATCH 03/70] Extend nonconforming to `bc` --- src/scimlfunctions.jl | 36 ++++++++++-------------------------- 1 file changed, 10 insertions(+), 26 deletions(-) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index 102e0c041..d31677fef 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -2233,11 +2233,7 @@ For more details on this argument, see the ODEFunction documentation. The fields of the BVPFunction type directly match the names of the inputs. """ struct BVPFunction{iip, specialize, F, BF, TMM, Ta, Tt, TJ, BCTJ, JVP, VJP, JP, - BCJP, SP, TW, TWt, - TPJ, - S, S2, S3, O, TCV, BCTCV, - SYS} <: - AbstractBVPFunction{iip} + BCJP, SP, TW, TWt, TPJ, S, S2, S3, O, TCV, BCTCV, SYS} <: AbstractBVPFunction{iip} f::F bc::BF mass_matrix::TMM @@ -3815,31 +3811,23 @@ function OptimizationFunction{iip}(f, adtype::AbstractADType = NoAD(); end function BVPFunction{iip, specialize}(f, bc; - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : - I, + 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, bcjac = __has_jac(bc) ? bc.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, - bcjac_prototype = __has_jac_prototype(bc) ? - bc.jac_prototype : - nothing, - sparsity = __has_sparsity(f) ? f.sparsity : - jac_prototype, + jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing, + bcjac_prototype = __has_jac_prototype(bc) ? bc.jac_prototype : nothing, + sparsity = __has_sparsity(f) ? f.sparsity : jac_prototype, Wfact = __has_Wfact(f) ? f.Wfact : nothing, Wfact_t = __has_Wfact_t(f) ? f.Wfact_t : nothing, 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, - observed = __has_observed(f) ? f.observed : - DEFAULT_OBSERVED, + paramsyms = __has_paramsyms(f) ? f.paramsyms : nothing, + observed = __has_observed(f) ? f.observed : DEFAULT_OBSERVED, colorvec = __has_colorvec(f) ? f.colorvec : nothing, bccolorvec = __has_colorvec(bc) ? bc.colorvec : nothing, sys = __has_sys(f) ? f.sys : nothing) where {iip, specialize} @@ -3892,17 +3880,13 @@ function BVPFunction{iip, specialize}(f, bc; Wfact_tiip = Wfact_t !== nothing ? isinplace(Wfact_t, 5, "Wfact_t", iip) : iip paramjaciip = paramjac !== nothing ? isinplace(paramjac, 4, "paramjac", iip) : iip - nonconforming = (jaciip, - tgradiip, - jvpiip, - vjpiip, - Wfactiip, - Wfact_tiip, + nonconforming = (bciip, jaciip, tgradiip, jvpiip, vjpiip, Wfactiip, Wfact_tiip, paramjaciip) .!= iip bc_nonconforming = bcjaciip .!= bciip if any(nonconforming) nonconforming = findall(nonconforming) - functions = ["jac", "bcjac", "tgrad", "jvp", "vjp", "Wfact", "Wfact_t", "paramjac"][nonconforming] + functions = ["bc", "jac", "bcjac", "tgrad", "jvp", "vjp", "Wfact", "Wfact_t", + "paramjac"][nonconforming] throw(NonconformingFunctionsError(functions)) end From e47627367861475f3625d40b91a3f38b106bb2d3 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 5 Sep 2023 18:11:55 -0400 Subject: [PATCH 04/70] Don't dump BVPFunction Information --- src/problems/bvp_problems.jl | 42 +++++++++++++++--------------------- src/scimlfunctions.jl | 18 ++++++++++------ 2 files changed, 28 insertions(+), 32 deletions(-) diff --git a/src/problems/bvp_problems.jl b/src/problems/bvp_problems.jl index a25ea8f74..d9787f9bc 100644 --- a/src/problems/bvp_problems.jl +++ b/src/problems/bvp_problems.jl @@ -89,15 +89,12 @@ struct BVProblem{uType, tType, isinplace, P, F, BF, PT, K} <: kwargs::K @add_kwonly function BVProblem{iip}(f::AbstractBVPFunction{iip}, bc, u0, tspan, - p = NullParameters(), - problem_type = StandardBVProblem(); - kwargs...) where {iip} + p = NullParameters(), problem_type = StandardBVProblem(); kwargs...) where {iip} _tspan = promote_tspan(tspan) warn_paramtype(p) - new{typeof(u0), typeof(_tspan), iip, typeof(p), - typeof(f.f), typeof(bc), - typeof(problem_type), typeof(kwargs)}(f.f, bc, u0, _tspan, p, - problem_type, kwargs) + return new{typeof(u0), typeof(_tspan), iip, typeof(p), typeof(f), typeof(bc), + typeof(problem_type), typeof(kwargs)}(f, bc, u0, _tspan, p, problem_type, + kwargs) end function BVProblem{iip}(f, bc, u0, tspan, p = NullParameters(); kwargs...) where {iip} @@ -108,11 +105,12 @@ end TruncatedStacktraces.@truncate_stacktrace BVProblem 3 1 2 function BVProblem(f, bc, u0, tspan, p = NullParameters(); kwargs...) - BVProblem(BVPFunction(f, bc), u0, tspan, p; kwargs...) + iip = isinplace(f, 4) + return BVProblem{iip}(BVPFunction{iip}(f, bc), u0, tspan, p; kwargs...) end function BVProblem(f::AbstractBVPFunction, u0, tspan, p = NullParameters(); kwargs...) - BVProblem{isinplace(f)}(f.f, f.bc, u0, tspan, p; kwargs...) + return BVProblem{isinplace(f)}(f, f.bc, u0, tspan, p; kwargs...) end """ @@ -132,31 +130,25 @@ end """ $(TYPEDEF) """ -struct TwoPointBVProblem{iip} end +struct TwoPointBVProblem end + function TwoPointBVProblem(f, bc, u0, tspan, p = NullParameters(); kwargs...) - iip = isinplace(f, 4) - TwoPointBVProblem{iip}(f, bc, u0, tspan, p; kwargs...) + return TwoPointBVProblem(BVPFunction(f, bc), u0, tspan, p; kwargs...) end -function TwoPointBVProblem{iip}(f, bc, u0, tspan, p = NullParameters(); - kwargs...) where {iip} - BVProblem{iip}(f, TwoPointBVPFunction(bc), u0, tspan, p; kwargs...) +function TwoPointBVProblem(f::AbstractBVPFunction, u0, tspan, p = NullParameters(); + kwargs...) + iip = isinplace(f) + return BVProblem{iip}(f, f.bc, u0, tspan, p, TwoPointBVProblem(); kwargs...) end # Allow previous timeseries solution -function TwoPointBVProblem(f::AbstractODEFunction, - bc, - sol::T, - tspan::Tuple, +function TwoPointBVProblem(f::AbstractODEFunction, bc, sol::T, tspan::Tuple, p = NullParameters()) where {T <: AbstractTimeseriesSolution} TwoPointBVProblem(f, bc, sol.u, tspan, p) end # Allow initial guess function for the initial guess -function TwoPointBVProblem(f::AbstractODEFunction, - bc, - initialGuess, - tspan::AbstractVector, - p = NullParameters(); - kwargs...) +function TwoPointBVProblem(f::AbstractODEFunction, bc, initialGuess, tspan::AbstractVector, + p = NullParameters(); kwargs...) u0 = [initialGuess(i) for i in tspan] TwoPointBVProblem(f, bc, u0, (tspan[1], tspan[end]), p) end diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index d31677fef..6cd7a6b4b 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -2233,7 +2233,7 @@ For more details on this argument, see the ODEFunction documentation. The fields of the BVPFunction type directly match the names of the inputs. """ struct BVPFunction{iip, specialize, F, BF, TMM, Ta, Tt, TJ, BCTJ, JVP, VJP, JP, - BCJP, SP, TW, TWt, TPJ, S, S2, S3, O, TCV, BCTCV, SYS} <: AbstractBVPFunction{iip} + BCJP, BCRP, SP, TW, TWt, TPJ, S, S2, S3, O, TCV, BCTCV, SYS} <: AbstractBVPFunction{iip} f::F bc::BF mass_matrix::TMM @@ -2245,6 +2245,7 @@ struct BVPFunction{iip, specialize, F, BF, TMM, Ta, Tt, TJ, BCTJ, JVP, VJP, JP, vjp::VJP jac_prototype::JP bcjac_prototype::BCJP + bcresid_prototype::BCRP sparsity::SP Wfact::TW Wfact_t::TWt @@ -3820,7 +3821,8 @@ function BVPFunction{iip, specialize}(f, bc; vjp = __has_vjp(f) ? f.vjp : nothing, jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing, bcjac_prototype = __has_jac_prototype(bc) ? bc.jac_prototype : nothing, - sparsity = __has_sparsity(f) ? f.sparsity : jac_prototype, + bcresid_prototype = nothing, + sparsity = __has_sparsity(f) ? f.sparsity : jac_prototype, Wfact = __has_Wfact(f) ? f.Wfact : nothing, Wfact_t = __has_Wfact_t(f) ? f.Wfact_t : nothing, paramjac = __has_paramjac(f) ? f.paramjac : nothing, @@ -3831,6 +3833,7 @@ function BVPFunction{iip, specialize}(f, bc; colorvec = __has_colorvec(f) ? f.colorvec : nothing, bccolorvec = __has_colorvec(bc) ? bc.colorvec : nothing, sys = __has_sys(f) ? f.sys : nothing) where {iip, specialize} + # TODO: Add checks for bcresid_prototype for TwoPointBVPFunction if mass_matrix === I && typeof(f) <: Tuple mass_matrix = ((I for i in 1:length(f))...,) end @@ -3898,14 +3901,14 @@ function BVPFunction{iip, specialize}(f, bc; if specialize === NoSpecialize BVPFunction{iip, specialize, Any, Any, Any, Any, Any, - Any, Any, Any, Any, Any, Any, Any, Any, Any, + Any, Any, Any, Any, Any, Any, Any, Any, Any, Any, Any, typeof(syms), typeof(indepsym), typeof(paramsyms), Any, typeof(_colorvec), typeof(_bccolorvec), Any}(f, bc, mass_matrix, analytic, tgrad, jac, bcjac, jvp, vjp, jac_prototype, - bcjac_prototype, + bcjac_prototype, bcresid_prototype, sparsity, Wfact, Wfact_t, paramjac, syms, @@ -3916,11 +3919,12 @@ function BVPFunction{iip, specialize}(f, bc; BVPFunction{iip, FunctionWrapperSpecialize, typeof(f), typeof(bc), typeof(mass_matrix), typeof(analytic), typeof(tgrad), typeof(jac), typeof(jvp), typeof(vjp), typeof(jac_prototype), + typeof(bcjac_prototype), typeof(bcresid_prototype), typeof(sparsity), typeof(Wfact), typeof(Wfact_t), typeof(paramjac), typeof(syms), typeof(indepsym), typeof(paramsyms), typeof(observed), typeof(_colorvec), typeof(_bccolorvec), typeof(sys)}(f, bc, mass_matrix, analytic, tgrad, jac, bcjac, - jvp, vjp, jac_prototype, bcjac_prototype, sparsity, Wfact, + jvp, vjp, jac_prototype, bcjac_prototype, bcresid_prototype, sparsity, Wfact, Wfact_t, paramjac, syms, indepsym, paramsyms, observed, _colorvec, _bccolorvec, sys) else @@ -3928,13 +3932,13 @@ function BVPFunction{iip, specialize}(f, bc; typeof(analytic), typeof(tgrad), typeof(jac), typeof(bcjac), typeof(jvp), typeof(vjp), typeof(jac_prototype), - typeof(bcjac_prototype), + typeof(bcjac_prototype), typeof(bcresid_prototype), typeof(sparsity), typeof(Wfact), typeof(Wfact_t), typeof(paramjac), typeof(syms), typeof(indepsym), typeof(paramsyms), typeof(observed), typeof(_colorvec), typeof(_bccolorvec), typeof(sys)}(f, bc, mass_matrix, analytic, tgrad, jac, bcjac, jvp, vjp, - jac_prototype, bcjac_prototype, sparsity, + jac_prototype, bcjac_prototype, bcresid_prototype, sparsity, Wfact, Wfact_t, paramjac, syms, indepsym, paramsyms, observed, _colorvec, _bccolorvec, sys) From 2b73aefd267af7f213d8fc2f5ef72e3083d213ae Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Wed, 6 Sep 2023 14:16:21 -0400 Subject: [PATCH 05/70] Fix TwoPointBVP Dispatches --- src/problems/bvp_problems.jl | 81 +++++++++++++++++++++--------------- src/scimlfunctions.jl | 71 ++++++++++++------------------- 2 files changed, 73 insertions(+), 79 deletions(-) diff --git a/src/problems/bvp_problems.jl b/src/problems/bvp_problems.jl index d9787f9bc..b46b0479c 100644 --- a/src/problems/bvp_problems.jl +++ b/src/problems/bvp_problems.jl @@ -3,6 +3,11 @@ $(TYPEDEF) """ struct StandardBVProblem end +""" +$(TYPEDEF) +""" +struct TwoPointBVProblem end + @doc doc""" Defines an BVP problem. @@ -17,7 +22,7 @@ condition ``u_0`` which define an ODE: \frac{du}{dt} = f(u,p,t) ``` -along with an implicit function `bc!` which defines the residual equation, where +along with an implicit function `bc` which defines the residual equation, where ```math bc(u,p,t) = 0 @@ -36,22 +41,27 @@ u(t_f) = b ### Constructors ```julia -TwoPointBVProblem{isinplace}(f,bc!,u0,tspan,p=NullParameters();kwargs...) -BVProblem{isinplace}(f,bc!,u0,tspan,p=NullParameters();kwargs...) +TwoPointBVProblem{isinplace}(f,bc,u0,tspan,p=NullParameters();kwargs...) +BVProblem{isinplace}(f,bc,u0,tspan,p=NullParameters();kwargs...) ``` or if we have an initial guess function `initialGuess(t)` for the given BVP, we can pass the initial guess to the problem constructors: ```julia -TwoPointBVProblem{isinplace}(f,bc!,initialGuess,tspan,p=NullParameters();kwargs...) -BVProblem{isinplace}(f,bc!,initialGuess,tspan,p=NullParameters();kwargs...) +TwoPointBVProblem{isinplace}(f,bc,initialGuess,tspan,p=NullParameters();kwargs...) +BVProblem{isinplace}(f,bc,initialGuess,tspan,p=NullParameters();kwargs...) ``` -For any BVP problem type, `bc!` is the inplace function: +For any BVP problem type, `bc` must be inplace if `f` is inplace. Otherwise it must be +out-of-place. + +If the bvp is a StandardBVProblem (also known as a Multi-Point BV Problem) it must define +either of the following functions ```julia bc!(residual, u, p, t) +residual = bc(u, p, t) ``` where `residual` computed from the current `u`. `u` is an array of solution values @@ -61,6 +71,16 @@ time points, and for shooting type methods `u=sol` the ODE solution. Note that all features of the `ODESolution` are present in this form. In both cases, the size of the residual matches the size of the initial condition. +If the bvp is a TwoPointBVProblem it must define either of the following functions + +```julia +bc!((resid_a, resid_b), (u_a, u_b), p) +resid_a, resid_b = bc((u_a, u_b), p) +``` + +where `resid_a` and `resid_b` are the residuals at the two endpoints, `u_a` and `u_b` are +the solution values at the two endpoints, and `p` are the parameters. + Parameters are optional, and if not given, then a `NullParameters()` singleton will be used which will throw nice errors if you try to index non-existent parameters. Any extra keyword arguments are passed on to the solvers. For example, @@ -88,10 +108,11 @@ struct BVProblem{uType, tType, isinplace, P, F, BF, PT, K} <: problem_type::PT kwargs::K - @add_kwonly function BVProblem{iip}(f::AbstractBVPFunction{iip}, bc, u0, tspan, - p = NullParameters(), problem_type = StandardBVProblem(); kwargs...) where {iip} + @add_kwonly function BVProblem{iip}(f::AbstractBVPFunction{iip, TP}, bc, u0, tspan, + p = NullParameters(); kwargs...) where {iip, TP} _tspan = promote_tspan(tspan) warn_paramtype(p) + problem_type = TP ? TwoPointBVProblem() : StandardBVProblem() return new{typeof(u0), typeof(_tspan), iip, typeof(p), typeof(f), typeof(bc), typeof(problem_type), typeof(kwargs)}(f, bc, u0, _tspan, p, problem_type, kwargs) @@ -113,42 +134,34 @@ function BVProblem(f::AbstractBVPFunction, u0, tspan, p = NullParameters(); kwar return BVProblem{isinplace(f)}(f, f.bc, u0, tspan, p; kwargs...) end -""" -$(TYPEDEF) -""" -struct TwoPointBVPFunction{bF} - bc::bF -end -TwoPointBVPFunction(; bc = error("No argument `bc`")) = TwoPointBVPFunction(bc) -function (f::TwoPointBVPFunction)(residuala, residualb, ua, ub, p) - return f.bc(residuala, residualb, ua, ub, p) -end -function (f::TwoPointBVPFunction)(residual::Tuple, u, p) - return f(residual[1], residual[2], u[1], u[end], p) -end +# This is mostly a fake stuct and isn't used anywhere +# But we need it for function calls like TwoPointBVProblem{iip}(...) = ... +struct TwoPointBVPFunction{iip} end -""" -$(TYPEDEF) -""" -struct TwoPointBVProblem end +@inline TwoPointBVPFunction(args...; kwargs...) = BVPFunction(args...; kwargs..., twopoint=true) +@inline function TwoPointBVPFunction{iip}(args...; kwargs...) where {iip} + return BVPFunction{iip}(args...; kwargs..., twopoint=true) +end -function TwoPointBVProblem(f, bc, u0, tspan, p = NullParameters(); kwargs...) - return TwoPointBVProblem(BVPFunction(f, bc), u0, tspan, p; kwargs...) +function TwoPointBVProblem(f, bc, u0, tspan, p = NullParameters(); + bcresid_prototype=nothing, kwargs...) + return TwoPointBVProblem(TwoPointBVPFunction(f, bc; bcresid_prototype), u0, tspan, p; + kwargs...) end -function TwoPointBVProblem(f::AbstractBVPFunction, u0, tspan, p = NullParameters(); - kwargs...) - iip = isinplace(f) - return BVProblem{iip}(f, f.bc, u0, tspan, p, TwoPointBVProblem(); kwargs...) +function TwoPointBVProblem(f::AbstractBVPFunction{iip, twopoint}, u0, tspan, + p = NullParameters(); kwargs...) where {iip, twopoint} + @assert twopoint "`TwoPointBVProblem` can only be used with a `TwoPointBVPFunction`. Instead of using `BVPFunction`, use `TwoPointBVPFunction` or pass a kwarg `twopoint=true` during the construction of the `BVPFunction`." + return BVProblem{iip}(f, f.bc, u0, tspan, p; kwargs...) end # Allow previous timeseries solution function TwoPointBVProblem(f::AbstractODEFunction, bc, sol::T, tspan::Tuple, - p = NullParameters()) where {T <: AbstractTimeseriesSolution} - TwoPointBVProblem(f, bc, sol.u, tspan, p) + p = NullParameters(); kwargs...) where {T <: AbstractTimeseriesSolution} + return TwoPointBVProblem(f, bc, sol.u, tspan, p; kwargs...) end # Allow initial guess function for the initial guess function TwoPointBVProblem(f::AbstractODEFunction, bc, initialGuess, tspan::AbstractVector, p = NullParameters(); kwargs...) u0 = [initialGuess(i) for i in tspan] - TwoPointBVProblem(f, bc, u0, (tspan[1], tspan[end]), p) + return TwoPointBVProblem(f, bc, u0, (tspan[1], tspan[end]), p; kwargs...) end diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index 38edeb094..4c7c0caba 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -2126,8 +2126,7 @@ TruncatedStacktraces.@truncate_stacktrace OptimizationFunction 1 2 """ $(TYPEDEF) """ -abstract type AbstractBVPFunction{iip} <: - AbstractDiffEqFunction{iip} end +abstract type AbstractBVPFunction{iip, twopoint} <: AbstractDiffEqFunction{iip} end @doc doc""" BVPFunction{iip,F,BF,TMM,Ta,Tt,TJ,BCTJ,JVP,VJP,JP,BCJP,SP,TW,TWt,TPJ,S,S2,S3,O,TCV,BCTCV} <: AbstractBVPFunction{iip,specialize} @@ -2232,12 +2231,9 @@ For more details on this argument, see the ODEFunction documentation. The fields of the BVPFunction type directly match the names of the inputs. """ -struct BVPFunction{iip, specialize, F, BF, TMM, Ta, Tt, TJ, BCTJ, JVP, VJP, JP, - BCJP, SP, TW, TWt, - TPJ, - S, S2, S3, O, TCV, BCTCV, - SYS} <: - AbstractBVPFunction{iip} +struct BVPFunction{iip, specialize, twopoint, F, BF, TMM, Ta, Tt, TJ, BCTJ, JVP, VJP, + JP, BCJP, BCRP, SP, TW, TWt, TPJ, S, S2, S3, O, TCV, BCTCV, + SYS} <: AbstractBVPFunction{iip, twopoint} f::F bc::BF mass_matrix::TMM @@ -3817,7 +3813,7 @@ function OptimizationFunction{iip}(f, adtype::AbstractADType = NoAD(); cons_expr, sys) end -function BVPFunction{iip, specialize}(f, bc; +function BVPFunction{iip, specialize, twopoint}(f, bc; mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, analytic = __has_analytic(f) ? f.analytic : nothing, tgrad = __has_tgrad(f) ? f.tgrad : nothing, @@ -3838,8 +3834,7 @@ function BVPFunction{iip, specialize}(f, bc; observed = __has_observed(f) ? f.observed : DEFAULT_OBSERVED, colorvec = __has_colorvec(f) ? f.colorvec : nothing, bccolorvec = __has_colorvec(bc) ? bc.colorvec : nothing, - sys = __has_sys(f) ? f.sys : nothing) where {iip, specialize} - # TODO: Add checks for bcresid_prototype for TwoPointBVPFunction + sys = __has_sys(f) ? f.sys : nothing) where {iip, specialize, twopoint} if mass_matrix === I && typeof(f) <: Tuple mass_matrix = ((I for i in 1:length(f))...,) end @@ -3879,7 +3874,7 @@ function BVPFunction{iip, specialize}(f, bc; _bccolorvec = bccolorvec end - bciip = isinplace(bc, 4, "bc", iip) + bciip = !twopoint ? isinplace(bc, 4, "bc", iip) : isinplace(bc, 3, "bc", iip) jaciip = jac !== nothing ? isinplace(jac, 4, "jac", iip) : iip bcjaciip = bcjac !== nothing ? isinplace(bcjac, 4, "bcjac", bciip) : bciip tgradiip = tgrad !== nothing ? isinplace(tgrad, 4, "tgrad", iip) : iip @@ -3899,6 +3894,12 @@ function BVPFunction{iip, specialize}(f, bc; throw(NonconformingFunctionsError(functions)) end + if iip && twopoint + if bcresid_prototype === nothing || length(bcresid_prototype) != 2 + error("bcresid_prototype must be a tuple / indexable collection of length 2 for a TwoPointBVPFunction") + end + end + if any(bc_nonconforming) bc_nonconforming = findall(bc_nonconforming) functions = ["bcjac"][bc_nonconforming] @@ -3906,41 +3907,21 @@ function BVPFunction{iip, specialize}(f, bc; end if specialize === NoSpecialize - BVPFunction{iip, specialize, Any, Any, Any, Any, Any, + BVPFunction{iip, specialize, twopoint, Any, Any, Any, Any, Any, Any, Any, Any, Any, Any, Any, Any, Any, Any, Any, Any, typeof(syms), typeof(indepsym), typeof(paramsyms), Any, typeof(_colorvec), typeof(_bccolorvec), Any}(f, bc, mass_matrix, - analytic, - tgrad, - jac, bcjac, jvp, vjp, - jac_prototype, + analytic, tgrad, jac, bcjac, jvp, vjp, jac_prototype, bcjac_prototype, bcresid_prototype, - sparsity, Wfact, - Wfact_t, - paramjac, syms, - indepsym, paramsyms, - observed, + sparsity, Wfact, Wfact_t, paramjac, syms, indepsym, paramsyms, observed, _colorvec, _bccolorvec, sys) - elseif specialize === false - BVPFunction{iip, FunctionWrapperSpecialize, - typeof(f), typeof(bc), typeof(mass_matrix), typeof(analytic), typeof(tgrad), - typeof(jac), typeof(jvp), typeof(vjp), typeof(jac_prototype), - typeof(sparsity), typeof(Wfact), typeof(Wfact_t), typeof(paramjac), - typeof(syms), typeof(indepsym), typeof(paramsyms), typeof(observed), - typeof(_colorvec), typeof(_bccolorvec), - typeof(sys)}(f, bc, mass_matrix, analytic, tgrad, jac, bcjac, - jvp, vjp, jac_prototype, bcjac_prototype, sparsity, Wfact, - Wfact_t, paramjac, syms, indepsym, paramsyms, - observed, _colorvec, _bccolorvec, sys) else - BVPFunction{iip, specialize, typeof(f), typeof(bc), typeof(mass_matrix), - typeof(analytic), - typeof(tgrad), - typeof(jac), typeof(bcjac), typeof(jvp), typeof(vjp), typeof(jac_prototype), - typeof(bcjac_prototype), typeof(bcresid_prototype), - typeof(sparsity), typeof(Wfact), typeof(Wfact_t), - typeof(paramjac), typeof(syms), typeof(indepsym), typeof(paramsyms), - typeof(observed), + BVPFunction{iip, specialize, twopoint, typeof(f), typeof(bc), typeof(mass_matrix), + typeof(analytic), typeof(tgrad), typeof(jac), typeof(bcjac), typeof(jvp), + typeof(vjp), typeof(jac_prototype), + typeof(bcjac_prototype), typeof(bcresid_prototype), typeof(sparsity), + typeof(Wfact), typeof(Wfact_t), typeof(paramjac), typeof(syms), + typeof(indepsym), typeof(paramsyms), typeof(observed), typeof(_colorvec), typeof(_bccolorvec), typeof(sys)}(f, bc, mass_matrix, analytic, tgrad, jac, bcjac, jvp, vjp, jac_prototype, bcjac_prototype, bcresid_prototype, sparsity, @@ -3950,12 +3931,12 @@ function BVPFunction{iip, specialize}(f, bc; end end -function BVPFunction{iip}(f, bc; kwargs...) where {iip} - BVPFunction{iip, FullSpecialize}(f, bc; kwargs...) +function BVPFunction{iip}(f, bc; twopoint::Bool=false, kwargs...) where {iip} + BVPFunction{iip, FullSpecialize, twopoint}(f, bc; kwargs...) end BVPFunction{iip}(f::BVPFunction, bc; kwargs...) where {iip} = f -function BVPFunction(f, bc; kwargs...) - BVPFunction{isinplace(f, 4), FullSpecialize}(f, bc; kwargs...) +function BVPFunction(f, bc; twopoint::Bool=false, kwargs...) + BVPFunction{isinplace(f, 4), FullSpecialize, twopoint}(f, bc; kwargs...) end BVPFunction(f::BVPFunction; kwargs...) = f From 07e569877035b1d2490684c92ad9794c865dca76 Mon Sep 17 00:00:00 2001 From: Ilian Pihlajamaa Date: Mon, 11 Sep 2023 17:18:20 +0200 Subject: [PATCH 06/70] add IntegralDataProblem --- src/SciMLBase.jl | 2 +- src/problems/basic_problems.jl | 53 ++++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+), 1 deletion(-) diff --git a/src/SciMLBase.jl b/src/SciMLBase.jl index e4d1e371f..b583be5ab 100644 --- a/src/SciMLBase.jl +++ b/src/SciMLBase.jl @@ -762,7 +762,7 @@ export LinearProblem, NonlinearProblem, IntervalNonlinearProblem, IntegralProblem, OptimizationProblem -export IntegralProblem +export IntegralDataProblem export DiscreteProblem, ImplicitDiscreteProblem export SteadyStateProblem, SteadyStateSolution diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 6cb329290..90f5bfefa 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -386,6 +386,59 @@ function IntegralProblem(f, lb, ub, args...; kwargs...) IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...) end +@doc doc""" + +Defines a sampled integral problem. +Documentation Page: https://docs.sciml.ai/Integrals/stable/ + +## Mathematical Specification of a Sampled Integral Problem + +Integral problems define integrals over a set of data `y` +with corresponding sampling points `x`. + +`y` is an `AbstractArray` and `x` is a `AbstractVector`, +whose elements define the space being integrated along the dimension of `y` that is integrated. +This dimension is given by `dim`. + +## Problem Type + +### Constructors + +IntegralProblem(y,x; + dim=1, kwargs...) + +- `y`: The integrand +- `x`: The sampled integration domain +- `dim`: The dimension of `y`being integrated +- kwargs: Keyword arguments copied to the solvers. + +### Fields + +The fields match the names of the constructor arguments. +""" +struct IntegralDataProblem{Y<:AbstractArray, X<:AbstractVector, DIMVAL} <: AbstractIntegralProblem{false} + x::X + y::Y + dim::DIMVAL + kwargs::K + @add_kwonly function IntegralProblem(x, y; + dim = 1, + kwargs...) + @assert length(x) > 1 + @assert isfinite(first(x)) "Infinite limits are not supported" + @assert isfinite(last(x)) "Infinite limits are not supported" + @assert size(y, dim) == length(x) "Integrand and grid must be of equal length along the integrated dimension" + @assert axes(y, dim) == axes(x,1) "Grid and integrand array must use the same indexing along integrated dimension" + new{typeof(x), typeof(y), Val{dim}, typeof(kwargs)}(x, y, Val(dim), kwargs) + end +end + +TruncatedStacktraces.@truncate_stacktrace IntegralDataProblem 1 4 + +function IntegralProblem(x::AbstractArray, y::AbstractVector; dim=1, kwargs...) + IntegralDataProblem(x, y; dim=dim) +end + struct QuadratureProblem end @deprecate QuadratureProblem(args...; kwargs...) IntegralProblem(args...; kwargs...) From 181847d4c6a7d639ca4e51b3fe7291630914c026 Mon Sep 17 00:00:00 2001 From: Ilian Pihlajamaa Date: Mon, 11 Sep 2023 17:55:06 +0200 Subject: [PATCH 07/70] fix mistake --- src/problems/basic_problems.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 90f5bfefa..dd94f21a0 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -416,7 +416,7 @@ IntegralProblem(y,x; The fields match the names of the constructor arguments. """ -struct IntegralDataProblem{Y<:AbstractArray, X<:AbstractVector, DIMVAL} <: AbstractIntegralProblem{false} +struct IntegralDataProblem{Y<:AbstractArray, X<:AbstractVector, DIMVAL, K} <: AbstractIntegralProblem{false} x::X y::Y dim::DIMVAL From 2fcfbb549097ecd446304b127b4e0cb032b135b7 Mon Sep 17 00:00:00 2001 From: Ilian Pihlajamaa Date: Mon, 11 Sep 2023 18:05:13 +0200 Subject: [PATCH 08/70] another small mistake --- src/problems/basic_problems.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index dd94f21a0..24b955965 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -421,7 +421,7 @@ struct IntegralDataProblem{Y<:AbstractArray, X<:AbstractVector, DIMVAL, K} <: Ab y::Y dim::DIMVAL kwargs::K - @add_kwonly function IntegralProblem(x, y; + @add_kwonly function IntegralDataProblem(x, y; dim = 1, kwargs...) @assert length(x) > 1 From cb2ad9fc1fe124983eb9b98f002557ba3d25051c Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 12 Sep 2023 13:52:02 -0400 Subject: [PATCH 09/70] Construct ArrayPartition --- Project.toml | 2 +- src/problems/bvp_problems.jl | 2 +- src/scimlfunctions.jl | 15 +++++---------- 3 files changed, 7 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index 951f6b7ca..7b8a2b09a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SciMLBase" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" authors = ["Chris Rackauckas and contributors"] -version = "2.0.0" +version = "1.96.0" # Remember to make it 2.0.0 [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" diff --git a/src/problems/bvp_problems.jl b/src/problems/bvp_problems.jl index b46b0479c..b6fe87866 100644 --- a/src/problems/bvp_problems.jl +++ b/src/problems/bvp_problems.jl @@ -127,7 +127,7 @@ TruncatedStacktraces.@truncate_stacktrace BVProblem 3 1 2 function BVProblem(f, bc, u0, tspan, p = NullParameters(); kwargs...) iip = isinplace(f, 4) - return BVProblem{iip}(BVPFunction{iip}(f, bc), u0, tspan, p; kwargs...) + return BVProblem{iip}(BVPFunction{iip}(f, bc), bc, u0, tspan, p; kwargs...) end function BVProblem(f::AbstractBVPFunction, u0, tspan, p = NullParameters(); kwargs...) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index 4c7c0caba..f7f2178bb 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -1884,10 +1884,7 @@ For more details on this argument, see the ODEFunction documentation. The fields of the NonlinearFunction type directly match the names of the inputs. """ struct NonlinearFunction{iip, specialize, F, TMM, Ta, Tt, TJ, JVP, VJP, JP, SP, TW, TWt, - TPJ, - S, S2, O, TCV, - SYS, -} <: AbstractNonlinearFunction{iip} + TPJ, S, S2, O, TCV, SYS} <: AbstractNonlinearFunction{iip} f::F mass_matrix::TMM analytic::Ta @@ -3646,10 +3643,7 @@ function NonlinearFunction{iip, specialize}(f; DEFAULT_OBSERVED_NO_TIME, colorvec = __has_colorvec(f) ? f.colorvec : nothing, - sys = __has_sys(f) ? f.sys : nothing) where { - iip, - specialize, -} + sys = __has_sys(f) ? f.sys : nothing) where {iip, specialize} if mass_matrix === I && typeof(f) <: Tuple mass_matrix = ((I for i in 1:length(f))...,) end @@ -3701,8 +3695,8 @@ function NonlinearFunction{iip, specialize}(f; typeof(Wfact_t), typeof(paramjac), typeof(syms), typeof(paramsyms), typeof(observed), - typeof(_colorvec), typeof(sys)}(f, mass_matrix, analytic, tgrad, - jac, + typeof(_colorvec), typeof(sys)}(f, mass_matrix, + analytic, tgrad, jac, jvp, vjp, jac_prototype, sparsity, Wfact, Wfact_t, paramjac, syms, @@ -3898,6 +3892,7 @@ function BVPFunction{iip, specialize, twopoint}(f, bc; if bcresid_prototype === nothing || length(bcresid_prototype) != 2 error("bcresid_prototype must be a tuple / indexable collection of length 2 for a TwoPointBVPFunction") end + bcresid_prototype = ArrayPartition(bcresid_prototype[1], bcresid_prototype[2]) end if any(bc_nonconforming) From 8b006588b62c420262a5ff212a88f9dcb36abec7 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 12 Sep 2023 18:07:17 -0400 Subject: [PATCH 10/70] Proper resid prototype handling --- src/scimlfunctions.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index 644545948..9da24fb14 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -3892,11 +3892,13 @@ function BVPFunction{iip, specialize, twopoint}(f, bc; throw(NonconformingFunctionsError(functions)) end - if iip && twopoint - if bcresid_prototype === nothing || length(bcresid_prototype) != 2 - error("bcresid_prototype must be a tuple / indexable collection of length 2 for a TwoPointBVPFunction") + if twopoint + if iip && (bcresid_prototype === nothing || length(bcresid_prototype) != 2) + error("bcresid_prototype must be a tuple / indexable collection of length 2 for a inplace TwoPointBVPFunction") + end + if bcresid_prototype !== nothing && length(bcresid_prototype) == 2 + bcresid_prototype = ArrayPartition(bcresid_prototype[1], bcresid_prototype[2]) end - bcresid_prototype = ArrayPartition(bcresid_prototype[1], bcresid_prototype[2]) end if any(bc_nonconforming) From 705a524b07251978fb125b20e6f003110f801492 Mon Sep 17 00:00:00 2001 From: "William S. Moses" Date: Tue, 12 Sep 2023 21:11:27 -0400 Subject: [PATCH 11/70] Make ODEProblem into a mutable struct --- src/problems/ode_problems.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/problems/ode_problems.jl b/src/problems/ode_problems.jl index 20949d644..a4dcbc544 100644 --- a/src/problems/ode_problems.jl +++ b/src/problems/ode_problems.jl @@ -94,7 +94,7 @@ prob = ODEProblemLibrary.prob_ode_linear sol = solve(prob) ``` """ -struct ODEProblem{uType, tType, isinplace, P, F, K, PT} <: +mutable struct ODEProblem{uType, tType, isinplace, P, F, K, PT} <: AbstractODEProblem{uType, tType, isinplace} """The ODE is `du = f(u,p,t)` for out-of-place and f(du,u,p,t) for in-place.""" f::F From dfd917a38c9dd3d4acec6c804932e6601752ebee Mon Sep 17 00:00:00 2001 From: Ilian Pihlajamaa Date: Wed, 13 Sep 2023 09:44:32 +0200 Subject: [PATCH 12/70] removed `IntegralDataProblem`, added x and dim field to `IntegralProblem` --- src/problems/basic_problems.jl | 85 ++++++++++++---------------------- 1 file changed, 29 insertions(+), 56 deletions(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 24b955965..877320ae5 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -334,10 +334,16 @@ 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. +IntegralProblem{iip}(f, x::AbstractVector{<:Number}, p = NullParameters(); + nout = 1, + batch = 0, kwargs...) +``` + +- f: the integrand, callable function `y = f(u,p)` for out-of-place or `f(y,u,p)` for in-place, or `AbstractArray` for sampled problems. - 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. @@ -353,6 +359,8 @@ IntegralProblem{iip}(f,lb,ub,p=NullParameters(); if `f` is a scalar valued function `y[i]` is the evaluation of `f` at a different point. Note that batch is a suggestion for the number of points, and it is not necessarily true that batch is the same as batchsize in all algorithms. +- x: Vector of gridpoints on which to evaluate the integrand for integrators that support this. Defaults to `nothing` +- dim: For integrals over sampled data (when `f isa AbstractArray`), this specifies the dimension along which is integrated. - kwargs: Keyword arguments copied to the solvers. Additionally, we can supply iip like IntegralProblem{iip}(...) as true or false to declare at @@ -362,21 +370,33 @@ compile time whether the integrator function is in-place. The fields match the names of the constructor arguments. """ -struct IntegralProblem{isinplace, P, F, B, K} <: AbstractIntegralProblem{isinplace} +struct IntegralProblem{isinplace, P, F, B, K, X, D} <: AbstractIntegralProblem{isinplace} f::F lb::B ub::B nout::Int p::P batch::Int + x::X + dim::D kwargs::K @add_kwonly function IntegralProblem{iip}(f, lb, ub, p = NullParameters(); nout = 1, - batch = 0, kwargs...) where {iip} + batch = 0, x=nothing, dim=Val(1), 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, - batch, kwargs) + new{iip, typeof(p), typeof(f), typeof(lb), typeof(x), typeof(dim) typeof(kwargs)}(f, lb, ub, nout, p, + batch, x, dim, kwargs) + end + @add_kwonly function IntegralProblem{iip}(f, x::AbstractVector{<:Number}, p = NullParameters(); + dim=Val(1), + nout = 1, + batch = 0, kwargs...) where {iip} + lb = x[begin] + ub = x[end] + warn_paramtype(p) + new{iip, typeof(p), typeof(f), typeof(lb), typeof(x), typeof(dim), typeof(kwargs)}(f, lb, ub, nout, p, + batch, x, dim, kwargs) end end @@ -385,58 +405,11 @@ TruncatedStacktraces.@truncate_stacktrace IntegralProblem 1 4 function IntegralProblem(f, lb, ub, args...; kwargs...) IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...) end - -@doc doc""" - -Defines a sampled integral problem. -Documentation Page: https://docs.sciml.ai/Integrals/stable/ - -## Mathematical Specification of a Sampled Integral Problem - -Integral problems define integrals over a set of data `y` -with corresponding sampling points `x`. - -`y` is an `AbstractArray` and `x` is a `AbstractVector`, -whose elements define the space being integrated along the dimension of `y` that is integrated. -This dimension is given by `dim`. - -## Problem Type - -### Constructors - -IntegralProblem(y,x; - dim=1, kwargs...) - -- `y`: The integrand -- `x`: The sampled integration domain -- `dim`: The dimension of `y`being integrated -- kwargs: Keyword arguments copied to the solvers. - -### Fields - -The fields match the names of the constructor arguments. -""" -struct IntegralDataProblem{Y<:AbstractArray, X<:AbstractVector, DIMVAL, K} <: AbstractIntegralProblem{false} - x::X - y::Y - dim::DIMVAL - kwargs::K - @add_kwonly function IntegralDataProblem(x, y; - dim = 1, - kwargs...) - @assert length(x) > 1 - @assert isfinite(first(x)) "Infinite limits are not supported" - @assert isfinite(last(x)) "Infinite limits are not supported" - @assert size(y, dim) == length(x) "Integrand and grid must be of equal length along the integrated dimension" - @assert axes(y, dim) == axes(x,1) "Grid and integrand array must use the same indexing along integrated dimension" - new{typeof(x), typeof(y), Val{dim}, typeof(kwargs)}(x, y, Val(dim), kwargs) - end +function IntegralProblem(f, x::AbstractVector{<:Number}, args...; kwargs...) + IntegralProblem{isinplace(f, 3)}(f, x, args...; kwargs...) end - -TruncatedStacktraces.@truncate_stacktrace IntegralDataProblem 1 4 - -function IntegralProblem(x::AbstractArray, y::AbstractVector; dim=1, kwargs...) - IntegralDataProblem(x, y; dim=dim) +function IntegralProblem(y::AbstractArray, x::AbstractVector{<:Number}, args...; dim::Int=1, kwargs...) + IntegralProblem{false}(y, x, args...; dim=Val(dim), kwargs...) end struct QuadratureProblem end From e60e921f5e7f2e538bbb6cc70275a41b492096ed Mon Sep 17 00:00:00 2001 From: Ilian Pihlajamaa Date: Wed, 13 Sep 2023 09:45:37 +0200 Subject: [PATCH 13/70] Stop exporting IntegralDataProblem --- src/SciMLBase.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/SciMLBase.jl b/src/SciMLBase.jl index b583be5ab..f52d73b2c 100644 --- a/src/SciMLBase.jl +++ b/src/SciMLBase.jl @@ -762,8 +762,6 @@ export LinearProblem, NonlinearProblem, IntervalNonlinearProblem, IntegralProblem, OptimizationProblem -export IntegralDataProblem - export DiscreteProblem, ImplicitDiscreteProblem export SteadyStateProblem, SteadyStateSolution export NoiseProblem From 0c9051442c14832e3980cfaba326dfca75bf77ca Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Thu, 14 Sep 2023 10:42:39 +1000 Subject: [PATCH 14/70] Update scimlfunctions.jl --- src/scimlfunctions.jl | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index 1867e1d58..45d727074 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -440,20 +440,20 @@ 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 @@ -461,7 +461,7 @@ 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) From 98f81240a6d54073d05ea0b7bce0388dacc08a6b Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Thu, 14 Sep 2023 07:27:32 -0400 Subject: [PATCH 15/70] Fix adjoints for ODESolution indexing --- Project.toml | 1 + src/SciMLBase.jl | 1 + src/solutions/zygote.jl | 2 +- 3 files changed, 3 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 2aba2b87b..0d085b4e1 100644 --- a/Project.toml +++ b/Project.toml @@ -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" diff --git a/src/SciMLBase.jl b/src/SciMLBase.jl index e4d1e371f..a0a9bc900 100644 --- a/src/SciMLBase.jl +++ b/src/SciMLBase.jl @@ -23,6 +23,7 @@ import TruncatedStacktraces import ADTypes: AbstractADType import ChainRulesCore import ZygoteRules: @adjoint +import FillArrays using Reexport using SciMLOperators diff --git a/src/solutions/zygote.jl b/src/solutions/zygote.jl index 08090bdbf..d41d07e0f 100644 --- a/src/solutions/zygote.jl +++ b/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 From f1328246781f86319c1e0fc378b686043acdb795 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 14 Sep 2023 08:13:37 -0400 Subject: [PATCH 16/70] Update Downstream.yml --- .github/workflows/Downstream.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml index 0854c0e04..4afc74b32 100644 --- a/.github/workflows/Downstream.yml +++ b/.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 From a92e7f701ebd2814f2bcb1a4efdafb6d202abd5e Mon Sep 17 00:00:00 2001 From: Ilian Pihlajamaa Date: Thu, 14 Sep 2023 14:18:50 +0200 Subject: [PATCH 17/70] fix typo in IntegralProblem --- src/problems/basic_problems.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 877320ae5..cc106f3ec 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -385,7 +385,7 @@ struct IntegralProblem{isinplace, P, F, B, K, X, D} <: AbstractIntegralProblem{i batch = 0, x=nothing, dim=Val(1), 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(x), typeof(dim) typeof(kwargs)}(f, lb, ub, nout, p, + new{iip, typeof(p), typeof(f), typeof(lb), typeof(x), typeof(dim), typeof(kwargs)}(f, lb, ub, nout, p, batch, x, dim, kwargs) end @add_kwonly function IntegralProblem{iip}(f, x::AbstractVector{<:Number}, p = NullParameters(); From af4cfb8e328f6b60420d26fb24ba50697901c47a Mon Sep 17 00:00:00 2001 From: Ilian Pihlajamaa Date: Thu, 14 Sep 2023 14:51:51 +0200 Subject: [PATCH 18/70] removed one unnecessary inner constructor --- src/problems/basic_problems.jl | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index cc106f3ec..7fe5a58fc 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -370,7 +370,7 @@ compile time whether the integrator function is in-place. The fields match the names of the constructor arguments. """ -struct IntegralProblem{isinplace, P, F, B, K, X, D} <: AbstractIntegralProblem{isinplace} +struct IntegralProblem{isinplace, P, F, B, X, D, K} <: AbstractIntegralProblem{isinplace} f::F lb::B ub::B @@ -382,22 +382,12 @@ struct IntegralProblem{isinplace, P, F, B, K, X, D} <: AbstractIntegralProblem{i kwargs::K @add_kwonly function IntegralProblem{iip}(f, lb, ub, p = NullParameters(); nout = 1, - batch = 0, x=nothing, dim=Val(1), kwargs...) where {iip} + batch = 0, x=nothing, dim=nothing, 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(x), typeof(dim), typeof(kwargs)}(f, lb, ub, nout, p, batch, x, dim, kwargs) end - @add_kwonly function IntegralProblem{iip}(f, x::AbstractVector{<:Number}, p = NullParameters(); - dim=Val(1), - nout = 1, - batch = 0, kwargs...) where {iip} - lb = x[begin] - ub = x[end] - warn_paramtype(p) - new{iip, typeof(p), typeof(f), typeof(lb), typeof(x), typeof(dim), typeof(kwargs)}(f, lb, ub, nout, p, - batch, x, dim, kwargs) - end end TruncatedStacktraces.@truncate_stacktrace IntegralProblem 1 4 @@ -406,10 +396,10 @@ function IntegralProblem(f, lb, ub, args...; kwargs...) IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...) end function IntegralProblem(f, x::AbstractVector{<:Number}, args...; kwargs...) - IntegralProblem{isinplace(f, 3)}(f, x, args...; kwargs...) + IntegralProblem(f, x[begin], x[end], args...; x=x, kwargs...) end function IntegralProblem(y::AbstractArray, x::AbstractVector{<:Number}, args...; dim::Int=1, kwargs...) - IntegralProblem{false}(y, x, args...; dim=Val(dim), kwargs...) + IntegralProblem{false}(y, x[begin], x[end], args...; dim=Val(dim), kwargs...) end struct QuadratureProblem end From d1585f0799c4607fd1493919241027d647c4af7a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 14 Sep 2023 16:08:00 -0400 Subject: [PATCH 19/70] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 0d085b4e1..e80a5a8ce 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SciMLBase" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" authors = ["Chris Rackauckas and contributors"] -version = "1.97.0" +version = "1.97.1" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 1827d300340d536b601ebf9f47a9c2336b779ca8 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 14 Sep 2023 16:20:05 -0400 Subject: [PATCH 20/70] Update Project.toml --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index e80a5a8ce..f5ab30163 100644 --- a/Project.toml +++ b/Project.toml @@ -46,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" From caedc7c5ad772a510ae23842420ec0262ad4176d Mon Sep 17 00:00:00 2001 From: Ilian Pihlajamaa Date: Fri, 15 Sep 2023 09:43:35 +0200 Subject: [PATCH 21/70] remove method ambiguity --- src/problems/basic_problems.jl | 5 ++++- test/function_building_error_messages.jl | 1 + 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 7fe5a58fc..1fde2ced6 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -392,7 +392,10 @@ end TruncatedStacktraces.@truncate_stacktrace IntegralProblem 1 4 -function IntegralProblem(f, lb, ub, args...; kwargs...) +function IntegralProblem(f, lb::AbstractVector{<:Number}, ub::AbstractVector{<:Number}, args...; kwargs...) + IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...) +end +function IntegralProblem(f, lb::Number, ub::Number, args...; kwargs...) IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...) end function IntegralProblem(f, x::AbstractVector{<:Number}, args...; kwargs...) diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index a150371d8..8a1dbdf47 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -455,6 +455,7 @@ intf(u) = 1.0 intf(u, p) = 1.0 IntegralProblem(intf, 0.0, 1.0) + # Optimization optf(u) = 1.0 From 7cee282d8fc3a5c21638f2e1613c5db53a976e48 Mon Sep 17 00:00:00 2001 From: Ilian Pihlajamaa Date: Fri, 15 Sep 2023 09:48:19 +0200 Subject: [PATCH 22/70] add problem building tests --- test/function_building_error_messages.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index 8a1dbdf47..ab8ddcbc4 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -453,8 +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 +x = [1.0,2.0] +y = rand(2,2) 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) +IntegralProblem(intf, x) +IntegralProblem(intf, x, p) +IntegralProblem(y, x) +IntegralProblem(y, x, p) # Optimization From 6b21a0d8b29e1013319ad9f2b3fa7457a4b7a4cd Mon Sep 17 00:00:00 2001 From: Ilian Pihlajamaa Date: Fri, 15 Sep 2023 09:58:47 +0200 Subject: [PATCH 23/70] ran the format_document.jl thing in VScode --- src/problems/basic_problems.jl | 21 +++++++++++++++------ test/function_building_error_messages.jl | 6 +++--- 2 files changed, 18 insertions(+), 9 deletions(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 1fde2ced6..3cf8bb819 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -382,27 +382,36 @@ struct IntegralProblem{isinplace, P, F, B, X, D, K} <: AbstractIntegralProblem{i kwargs::K @add_kwonly function IntegralProblem{iip}(f, lb, ub, p = NullParameters(); nout = 1, - batch = 0, x=nothing, dim=nothing, kwargs...) where {iip} + batch = 0, x = nothing, dim = nothing, 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(x), typeof(dim), typeof(kwargs)}(f, lb, ub, nout, p, + new{iip, typeof(p), typeof(f), typeof(lb), typeof(x), typeof(dim), typeof(kwargs)}(f, + lb, ub, nout, p, batch, x, dim, kwargs) end end TruncatedStacktraces.@truncate_stacktrace IntegralProblem 1 4 -function IntegralProblem(f, lb::AbstractVector{<:Number}, ub::AbstractVector{<:Number}, args...; kwargs...) +function IntegralProblem(f, + lb::AbstractVector{<:Number}, + ub::AbstractVector{<:Number}, + args...; + kwargs...) IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...) end function IntegralProblem(f, lb::Number, ub::Number, args...; kwargs...) IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...) end function IntegralProblem(f, x::AbstractVector{<:Number}, args...; kwargs...) - IntegralProblem(f, x[begin], x[end], args...; x=x, kwargs...) + IntegralProblem(f, x[begin], x[end], args...; x = x, kwargs...) end -function IntegralProblem(y::AbstractArray, x::AbstractVector{<:Number}, args...; dim::Int=1, kwargs...) - IntegralProblem{false}(y, x[begin], x[end], args...; dim=Val(dim), kwargs...) +function IntegralProblem(y::AbstractArray, + x::AbstractVector{<:Number}, + args...; + dim::Int = 1, + kwargs...) + IntegralProblem{false}(y, x[begin], x[end], args...; dim = Val(dim), kwargs...) end struct QuadratureProblem end diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index ab8ddcbc4..3a7b339d7 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -454,10 +454,10 @@ intf(u) = 1.0 @test_throws SciMLBase.TooFewArgumentsError IntegralProblem(intf, 0.0, 1.0) intf(u, p) = 1.0 p = 2.0 -x = [1.0,2.0] -y = rand(2,2) +x = [1.0, 2.0] +y = rand(2, 2) IntegralProblem(intf, 0.0, 1.0) -IntegralProblem(intf, 0.0, 1.0, p) +IntegralProblem(intf, 0.0, 1.0, p) IntegralProblem(intf, [0.0], [1.0]) IntegralProblem(intf, [0.0], [1.0], p) IntegralProblem(intf, x) From 09561b62e56c39a89b5dea18a5b86d404d6fb5e5 Mon Sep 17 00:00:00 2001 From: Ilian Pihlajamaa Date: Fri, 15 Sep 2023 11:33:42 +0200 Subject: [PATCH 24/70] Update basic_problems.jl --- src/problems/basic_problems.jl | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 3cf8bb819..63ad114e0 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -393,17 +393,24 @@ end TruncatedStacktraces.@truncate_stacktrace IntegralProblem 1 4 -function IntegralProblem(f, +function IntegralProblem(f::Function, lb::AbstractVector{<:Number}, ub::AbstractVector{<:Number}, args...; kwargs...) IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...) end -function IntegralProblem(f, lb::Number, ub::Number, args...; kwargs...) +function IntegralProblem(f::Function, + lb::Number, + ub::Number, + args...; + kwargs...) IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...) end -function IntegralProblem(f, x::AbstractVector{<:Number}, args...; kwargs...) +function IntegralProblem(f::Function, + x::AbstractVector{<:Number}, + args...; + kwargs...) IntegralProblem(f, x[begin], x[end], args...; x = x, kwargs...) end function IntegralProblem(y::AbstractArray, @@ -411,7 +418,7 @@ function IntegralProblem(y::AbstractArray, args...; dim::Int = 1, kwargs...) - IntegralProblem{false}(y, x[begin], x[end], args...; dim = Val(dim), kwargs...) + IntegralProblem{false}(y, x[begin], x[end], p = nothing; dim = Val(dim), kwargs...) end struct QuadratureProblem end From cdc36635707177d6118aa560c4823ba6faa658b8 Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Sat, 16 Sep 2023 21:52:03 +0200 Subject: [PATCH 25/70] revert to DataIntegralProblem --- src/SciMLBase.jl | 2 +- src/problems/basic_problems.jl | 92 ++++++++++++++---------- test/function_building_error_messages.jl | 12 ++-- 3 files changed, 62 insertions(+), 44 deletions(-) diff --git a/src/SciMLBase.jl b/src/SciMLBase.jl index f52d73b2c..1de4703d6 100644 --- a/src/SciMLBase.jl +++ b/src/SciMLBase.jl @@ -760,7 +760,7 @@ export solve, solve!, init, discretize, symbolic_discretize export LinearProblem, NonlinearProblem, IntervalNonlinearProblem, - IntegralProblem, OptimizationProblem + IntegralProblem, DataIntegralProblem, OptimizationProblem export DiscreteProblem, ImplicitDiscreteProblem export SteadyStateProblem, SteadyStateSolution diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 63ad114e0..f794b0e16 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -337,13 +337,9 @@ which are `Number`s or `AbstractVector`s with the same geometry as `u`. ``` IntegralProblem{iip}(f,lb,ub,p=NullParameters(); nout=1, batch = 0, kwargs...) - -IntegralProblem{iip}(f, x::AbstractVector{<:Number}, p = NullParameters(); - nout = 1, - batch = 0, kwargs...) ``` -- f: the integrand, callable function `y = f(u,p)` for out-of-place or `f(y,u,p)` for in-place, or `AbstractArray` for sampled problems. +- 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. @@ -359,8 +355,6 @@ IntegralProblem{iip}(f, x::AbstractVector{<:Number}, p = NullParameters(); if `f` is a scalar valued function `y[i]` is the evaluation of `f` at a different point. Note that batch is a suggestion for the number of points, and it is not necessarily true that batch is the same as batchsize in all algorithms. -- x: Vector of gridpoints on which to evaluate the integrand for integrators that support this. Defaults to `nothing` -- dim: For integrals over sampled data (when `f isa AbstractArray`), this specifies the dimension along which is integrated. - kwargs: Keyword arguments copied to the solvers. Additionally, we can supply iip like IntegralProblem{iip}(...) as true or false to declare at @@ -370,62 +364,86 @@ compile time whether the integrator function is in-place. The fields match the names of the constructor arguments. """ -struct IntegralProblem{isinplace, P, F, B, X, D, K} <: AbstractIntegralProblem{isinplace} +struct IntegralProblem{isinplace, P, F, B, K} <: AbstractIntegralProblem{isinplace} f::F lb::B ub::B nout::Int p::P batch::Int - x::X - dim::D kwargs::K @add_kwonly function IntegralProblem{iip}(f, lb, ub, p = NullParameters(); nout = 1, - batch = 0, x = nothing, dim = nothing, kwargs...) where {iip} + 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(x), typeof(dim), typeof(kwargs)}(f, + new{iip, typeof(p), typeof(f), typeof(lb), typeof(kwargs)}(f, lb, ub, nout, p, - batch, x, dim, kwargs) + batch, kwargs) end end TruncatedStacktraces.@truncate_stacktrace IntegralProblem 1 4 -function IntegralProblem(f::Function, - lb::AbstractVector{<:Number}, - ub::AbstractVector{<:Number}, - args...; - kwargs...) - IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...) -end -function IntegralProblem(f::Function, - lb::Number, - ub::Number, - args...; +function IntegralProblem(f, lb, ub, args...; kwargs...) IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...) end -function IntegralProblem(f::Function, - x::AbstractVector{<:Number}, - args...; - kwargs...) - IntegralProblem(f, x[begin], x[end], args...; x = x, kwargs...) -end -function IntegralProblem(y::AbstractArray, - x::AbstractVector{<:Number}, - args...; - dim::Int = 1, - kwargs...) - IntegralProblem{false}(y, x[begin], x[end], p = nothing; dim = Val(dim), kwargs...) -end struct QuadratureProblem end @deprecate QuadratureProblem(args...; kwargs...) IntegralProblem(args...; kwargs...) @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 + +Data 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 + +``` +DataIntegralProblem(x::AbstractVector, y::AbstractArray; dim=1, 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. +- kwargs: Keyword arguments copied to the solvers. + +### Fields + +The fields match the names of the constructor arguments. +""" +struct DataIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false} + y::Y + x::X + dim::D + kwargs::K + @add_kwonly function DataIntegralProblem(y::AbstractArray, x::AbstractVector; + dim = 1, + kwargs...) + @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), Val{dim}, typeof(kwargs)}(y, x, Val(dim), kwargs) + end +end + +TruncatedStacktraces.@truncate_stacktrace DataIntegralProblem 1 4 + +@doc doc""" + Defines an optimization problem. Documentation Page: https://docs.sciml.ai/Optimization/stable/API/optimization_problem/ diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index 3a7b339d7..8fffc806b 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -454,16 +454,16 @@ intf(u) = 1.0 @test_throws SciMLBase.TooFewArgumentsError IntegralProblem(intf, 0.0, 1.0) intf(u, p) = 1.0 p = 2.0 -x = [1.0, 2.0] -y = rand(2, 2) + 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) -IntegralProblem(intf, x) -IntegralProblem(intf, x, p) -IntegralProblem(y, x) -IntegralProblem(y, x, p) + +x = [1.0, 2.0] +y = rand(2, 2) +DataIntegralProblem(y, x) +DataIntegralProblem(y, x; dim=2) # Optimization From d3619dba982275026aa1b0e91ecc3b627cd0599b Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Sat, 16 Sep 2023 21:57:18 +0200 Subject: [PATCH 26/70] add additional assert --- src/problems/basic_problems.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index f794b0e16..b2caeabe9 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -434,6 +434,7 @@ struct DataIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false} @add_kwonly function DataIntegralProblem(y::AbstractArray, x::AbstractVector; dim = 1, 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), Val{dim}, typeof(kwargs)}(y, x, Val(dim), kwargs) From 124d3bec83466215bab456c3106e94be25e07481 Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Sat, 16 Sep 2023 22:45:01 +0200 Subject: [PATCH 27/70] rename `DataIntegralProblem` to `SampledIntegralProblem` --- src/SciMLBase.jl | 2 +- src/problems/basic_problems.jl | 10 +++++----- test/function_building_error_messages.jl | 4 ++-- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/SciMLBase.jl b/src/SciMLBase.jl index 1de4703d6..f0811160f 100644 --- a/src/SciMLBase.jl +++ b/src/SciMLBase.jl @@ -760,7 +760,7 @@ export solve, solve!, init, discretize, symbolic_discretize export LinearProblem, NonlinearProblem, IntervalNonlinearProblem, - IntegralProblem, DataIntegralProblem, OptimizationProblem + IntegralProblem, SampledIntegralProblem, OptimizationProblem export DiscreteProblem, ImplicitDiscreteProblem export SteadyStateProblem, SteadyStateSolution diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index b2caeabe9..0217c3f37 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -400,7 +400,7 @@ Documentation Page: https://docs.sciml.ai/Integrals/stable/ ## Mathematical Specification of a data Integral Problem -Data integral problems are defined as: +Sampled integral problems are defined as: ```math \sum_i w_i y_i @@ -413,7 +413,7 @@ assigned by a quadrature rule, which depend on sampling points `x`. ### Constructors ``` -DataIntegralProblem(x::AbstractVector, y::AbstractArray; dim=1, kwargs...) +SampledIntegralProblem(x::AbstractVector, y::AbstractArray; dim=1, kwargs...) ``` - y: The sampled integrand, must be a subtype of `AbstractArray`. It is assumed that the values of `y` along dimension `dim` @@ -426,12 +426,12 @@ DataIntegralProblem(x::AbstractVector, y::AbstractArray; dim=1, kwargs...) The fields match the names of the constructor arguments. """ -struct DataIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false} +struct SampledIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false} y::Y x::X dim::D kwargs::K - @add_kwonly function DataIntegralProblem(y::AbstractArray, x::AbstractVector; + @add_kwonly function SampledIntegralProblem(y::AbstractArray, x::AbstractVector; dim = 1, kwargs...) @assert dim <= ndims(y) "The integration dimension `dim` is larger than the number of dimensions of the integrand `y`" @@ -441,7 +441,7 @@ struct DataIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false} end end -TruncatedStacktraces.@truncate_stacktrace DataIntegralProblem 1 4 +TruncatedStacktraces.@truncate_stacktrace SampledIntegralProblem 1 4 @doc doc""" diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index 8fffc806b..c424a6f84 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -462,8 +462,8 @@ IntegralProblem(intf, [0.0], [1.0], p) x = [1.0, 2.0] y = rand(2, 2) -DataIntegralProblem(y, x) -DataIntegralProblem(y, x; dim=2) +SampledIntegralProblem(y, x) +SampledIntegralProblem(y, x; dim=2) # Optimization From 1e960d185092f00a4fb460b0573f7e65d6cdcb31 Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Sun, 17 Sep 2023 11:28:26 +0200 Subject: [PATCH 28/70] change default dimension to `ndims(y)` --- src/problems/basic_problems.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 0217c3f37..e0385140d 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -413,13 +413,13 @@ assigned by a quadrature rule, which depend on sampling points `x`. ### Constructors ``` -SampledIntegralProblem(x::AbstractVector, y::AbstractArray; dim=1, kwargs...) +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. +- dim: Dimension along which to integrate. Defaults to the last dimension of `y`. - kwargs: Keyword arguments copied to the solvers. ### Fields @@ -432,7 +432,7 @@ struct SampledIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false} dim::D kwargs::K @add_kwonly function SampledIntegralProblem(y::AbstractArray, x::AbstractVector; - dim = 1, + 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." From 8d293bfc0b859b1ae870772aaf023d2a93f6eb45 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 17 Sep 2023 09:53:57 -0400 Subject: [PATCH 29/70] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f5ab30163..fa982a687 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SciMLBase" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" authors = ["Chris Rackauckas and contributors"] -version = "1.97.1" +version = "1.98.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From ac73376f443bc832071ee54ee257ea978ac5db22 Mon Sep 17 00:00:00 2001 From: IlianPihlajamaa <73794090+IlianPihlajamaa@users.noreply.github.com> Date: Mon, 18 Sep 2023 21:14:40 +0200 Subject: [PATCH 30/70] remove Value type requirement from SampledIntegralProblem --- src/problems/basic_problems.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index e0385140d..1dc56ec9f 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -437,7 +437,7 @@ struct SampledIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false} @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), Val{dim}, typeof(kwargs)}(y, x, Val(dim), kwargs) + new{typeof(y), typeof(x), typeof(dim)}, typeof(kwargs)}(y, x, dim, kwargs) end end From 57ba9547dd148c2378c3acfaa24ac23ba99de956 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 16 Sep 2023 19:29:03 -0400 Subject: [PATCH 31/70] add integrand interface --- src/integrand_interface.jl | 55 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 src/integrand_interface.jl diff --git a/src/integrand_interface.jl b/src/integrand_interface.jl new file mode 100644 index 000000000..31e67249c --- /dev/null +++ b/src/integrand_interface.jl @@ -0,0 +1,55 @@ +""" + InplaceIntegrand(f!, result::AbstractArray) + +Constructor for a `InplaceIntegrand` accepting an integrand of the form `f!(y,x,p)`. The +caller also provides an output array needed to store the result of the quadrature. +Intermediate `y` arrays are allocated during the calculation, and the final result is is +written to `result`, so use the IntegralSolution immediately after the calculation to read +the result, and don't expect it to persist if the same integrand is used for another +calculation. +""" +struct InplaceIntegrand{F,T<:AbstractArray} + # in-place function f!(y, x, p) that takes one x value and outputs an array of results in-place + f!::F + I::T +end + + +""" + BatchIntegrand(f!, y::AbstractArray, x::AbstractVector, max_batch=typemax(Int)) + +Constructor for a `BatchIntegrand` accepting an integrand of the form `f!(y,x,p) = y .= f!.(x, Ref(p))` +that can evaluate the integrand at multiple quadrature nodes using, for example, threads, +the GPU, or distributed-memory. The `max_batch` keyword is a soft limit on the number of +nodes passed to the integrand. The buffers `y,x` must both be `resize!`-able since the +number of evaluation points may vary between calls to `f!`. +""" +struct BatchIntegrand{F,Y,X} + # in-place function f!(y, x, p) that takes an array of x values and outputs an array of results in-place + f!::F + y::Y + x::X + max_batch::Int # maximum number of x to supply in parallel + function BatchIntegrand(f!, y::AbstractArray, x::AbstractVector, max_batch::Integer=typemax(Int)) + max_batch > 0 || throw(ArgumentError("maximum batch size must be positive")) + return new{typeof(f!),typeof(y),typeof(x)}(f!, y, x, max_batch) + end +end + + +""" + BatchIntegrand(f!, y, x; max_batch=typemax(Int)) + +Constructor for a `BatchIntegrand` with pre-allocated buffers. +""" +BatchIntegrand(f!, y, x; max_batch::Integer=typemax(Int)) = + BatchIntegrand(f!, y, x, max_batch) + +""" + BatchIntegrand(f!, y::Type, x::Type=Nothing; max_batch=typemax(Int)) + +Constructor for a `BatchIntegrand` whose range type is known. The domain type is optional. +Array buffers for those types are allocated internally. +""" +BatchIntegrand(f!, Y::Type, X::Type=Nothing; max_batch::Integer=typemax(Int)) = + BatchIntegrand(f!, Y[], X[], max_batch) From 21b789516179ef2670652d545a1b98c786419768 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 16 Sep 2023 19:51:31 -0400 Subject: [PATCH 32/70] add InplaceBatchIntegrand --- src/integrand_interface.jl | 35 ++++++++++++++++++++++++++++++++++- 1 file changed, 34 insertions(+), 1 deletion(-) diff --git a/src/integrand_interface.jl b/src/integrand_interface.jl index 31e67249c..c091f2366 100644 --- a/src/integrand_interface.jl +++ b/src/integrand_interface.jl @@ -30,7 +30,7 @@ struct BatchIntegrand{F,Y,X} y::Y x::X max_batch::Int # maximum number of x to supply in parallel - function BatchIntegrand(f!, y::AbstractArray, x::AbstractVector, max_batch::Integer=typemax(Int)) + function BatchIntegrand(f!, y::AbstractVector, x::AbstractVector, max_batch::Integer=typemax(Int)) max_batch > 0 || throw(ArgumentError("maximum batch size must be positive")) return new{typeof(f!),typeof(y),typeof(x)}(f!, y, x, max_batch) end @@ -53,3 +53,36 @@ Array buffers for those types are allocated internally. """ BatchIntegrand(f!, Y::Type, X::Type=Nothing; max_batch::Integer=typemax(Int)) = BatchIntegrand(f!, Y[], X[], max_batch) + +""" + InplaceBatchIntegrand(f!, result::AbstractArray, y::AbstractArray, x::AbstractVector, max_batch=typemax(Int)) + +Constructor for a `InplaceBatchIntegrand` accepting an integrand of the form `f!(y,x,p) = y +.= f!.(x, Ref(p))` that can evaluate an inplace, array-valued integrand at multiple +quadrature nodes simultaneously using, for example, threads, the GPU, or distributed-memory. +The `max_batch` keyword is a soft limit on the number of nodes passed to the integrand. The +buffers `y,x` must both be `resize!`-able since the number of evaluation points may vary +between calls to `f!`. In particular, for a resizeable `y` buffer see ElasticArrays.jl . The +solution is written inplace to `result`. +""" +struct InplaceBatchIntegrand{F,T,Y,X} + # in-place function f!(y, x, p) that takes an array of x values and outputs an array of results in-place + f!::F + I::T + y::Y + x::X + max_batch::Int # maximum number of x to supply in parallel + function InplaceBatchIntegrand(f!, I::AbstractArray, y::AbstractArray, x::AbstractVector, max_batch::Integer=typemax(Int)) + max_batch > 0 || throw(ArgumentError("maximum batch size must be positive")) + return new{typeof(f!),typeof(I),typeof(y),typeof(x)}(f!, I, y, x, max_batch) + end +end + + +""" + InplaceBatchIntegrand(f!, result, y, x; max_batch=typemax(Int)) + +Constructor for a `InplaceBatchIntegrand` with pre-allocated buffers. +""" +InplaceBatchIntegrand(f!, result, y, x; max_batch::Integer=typemax(Int)) = + InplaceBatchIntegrand(f!, result, y, x, max_batch) From 6a2038ae520a0c3644c3b505b28be1f0cbddce35 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 16 Sep 2023 20:06:50 -0400 Subject: [PATCH 33/70] format and include --- src/SciMLBase.jl | 1 + src/integrand_interface.jl | 33 ++++++++++++++++++++------------- 2 files changed, 21 insertions(+), 13 deletions(-) diff --git a/src/SciMLBase.jl b/src/SciMLBase.jl index 011dfd690..b28aa5422 100644 --- a/src/SciMLBase.jl +++ b/src/SciMLBase.jl @@ -714,6 +714,7 @@ include("ensemble/ensemble_analysis.jl") include("solve.jl") include("interpolation.jl") +include("integrand_interface.jl") include("integrator_interface.jl") include("tabletraits.jl") include("remake.jl") diff --git a/src/integrand_interface.jl b/src/integrand_interface.jl index c091f2366..7b7ebc18e 100644 --- a/src/integrand_interface.jl +++ b/src/integrand_interface.jl @@ -8,13 +8,12 @@ written to `result`, so use the IntegralSolution immediately after the calculati the result, and don't expect it to persist if the same integrand is used for another calculation. """ -struct InplaceIntegrand{F,T<:AbstractArray} +struct InplaceIntegrand{F, T <: AbstractArray} # in-place function f!(y, x, p) that takes one x value and outputs an array of results in-place f!::F I::T end - """ BatchIntegrand(f!, y::AbstractArray, x::AbstractVector, max_batch=typemax(Int)) @@ -24,26 +23,29 @@ the GPU, or distributed-memory. The `max_batch` keyword is a soft limit on the n nodes passed to the integrand. The buffers `y,x` must both be `resize!`-able since the number of evaluation points may vary between calls to `f!`. """ -struct BatchIntegrand{F,Y,X} +struct BatchIntegrand{F, Y, X} # in-place function f!(y, x, p) that takes an array of x values and outputs an array of results in-place f!::F y::Y x::X max_batch::Int # maximum number of x to supply in parallel - function BatchIntegrand(f!, y::AbstractVector, x::AbstractVector, max_batch::Integer=typemax(Int)) + function BatchIntegrand(f!, + y::AbstractVector, + x::AbstractVector, + max_batch::Integer = typemax(Int)) max_batch > 0 || throw(ArgumentError("maximum batch size must be positive")) - return new{typeof(f!),typeof(y),typeof(x)}(f!, y, x, max_batch) + return new{typeof(f!), typeof(y), typeof(x)}(f!, y, x, max_batch) end end - """ BatchIntegrand(f!, y, x; max_batch=typemax(Int)) Constructor for a `BatchIntegrand` with pre-allocated buffers. """ -BatchIntegrand(f!, y, x; max_batch::Integer=typemax(Int)) = +function BatchIntegrand(f!, y, x; max_batch::Integer = typemax(Int)) BatchIntegrand(f!, y, x, max_batch) +end """ BatchIntegrand(f!, y::Type, x::Type=Nothing; max_batch=typemax(Int)) @@ -51,8 +53,9 @@ BatchIntegrand(f!, y, x; max_batch::Integer=typemax(Int)) = Constructor for a `BatchIntegrand` whose range type is known. The domain type is optional. Array buffers for those types are allocated internally. """ -BatchIntegrand(f!, Y::Type, X::Type=Nothing; max_batch::Integer=typemax(Int)) = +function BatchIntegrand(f!, Y::Type, X::Type = Nothing; max_batch::Integer = typemax(Int)) BatchIntegrand(f!, Y[], X[], max_batch) +end """ InplaceBatchIntegrand(f!, result::AbstractArray, y::AbstractArray, x::AbstractVector, max_batch=typemax(Int)) @@ -65,24 +68,28 @@ buffers `y,x` must both be `resize!`-able since the number of evaluation points between calls to `f!`. In particular, for a resizeable `y` buffer see ElasticArrays.jl . The solution is written inplace to `result`. """ -struct InplaceBatchIntegrand{F,T,Y,X} +struct InplaceBatchIntegrand{F, T, Y, X} # in-place function f!(y, x, p) that takes an array of x values and outputs an array of results in-place f!::F I::T y::Y x::X max_batch::Int # maximum number of x to supply in parallel - function InplaceBatchIntegrand(f!, I::AbstractArray, y::AbstractArray, x::AbstractVector, max_batch::Integer=typemax(Int)) + function InplaceBatchIntegrand(f!, + I::AbstractArray, + y::AbstractArray, + x::AbstractVector, + max_batch::Integer = typemax(Int)) max_batch > 0 || throw(ArgumentError("maximum batch size must be positive")) - return new{typeof(f!),typeof(I),typeof(y),typeof(x)}(f!, I, y, x, max_batch) + return new{typeof(f!), typeof(I), typeof(y), typeof(x)}(f!, I, y, x, max_batch) end end - """ InplaceBatchIntegrand(f!, result, y, x; max_batch=typemax(Int)) Constructor for a `InplaceBatchIntegrand` with pre-allocated buffers. """ -InplaceBatchIntegrand(f!, result, y, x; max_batch::Integer=typemax(Int)) = +function InplaceBatchIntegrand(f!, result, y, x; max_batch::Integer = typemax(Int)) InplaceBatchIntegrand(f!, result, y, x, max_batch) +end From 3f777598d2bd16793b5db18b122e38ce72c26cf5 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 18 Sep 2023 22:39:11 -0400 Subject: [PATCH 34/70] make the IntegralFunctions --- src/SciMLBase.jl | 16 +- src/integrand_interface.jl | 95 --------- src/scimlfunctions.jl | 256 ++++++++++++++++++++++- test/function_building_error_messages.jl | 17 ++ 4 files changed, 285 insertions(+), 99 deletions(-) delete mode 100644 src/integrand_interface.jl diff --git a/src/SciMLBase.jl b/src/SciMLBase.jl index b28aa5422..d206ca07d 100644 --- a/src/SciMLBase.jl +++ b/src/SciMLBase.jl @@ -589,6 +589,14 @@ abstract type AbstractDiffEqFunction{iip} <: """ $(TYPEDEF) +Base for types defining integrand functions. +""" +abstract type AbstractIntegralFunction{iip} <: + AbstractSciMLFunction{iip} end + +""" +$(TYPEDEF) + Base for types defining optimization functions. """ abstract type AbstractOptimizationFunction{iip} <: AbstractSciMLFunction{iip} end @@ -659,7 +667,9 @@ function specialization(::Union{ODEFunction{iip, specialize}, RODEFunction{iip, specialize}, NonlinearFunction{iip, specialize}, OptimizationFunction{iip, specialize}, - BVPFunction{iip, specialize}}) where {iip, + BVPFunction{iip, specialize}, + IntegralFunction{iip, specialize}, + BatchIntegralFunction{iip, specialize}}) where {iip, specialize} specialize end @@ -714,7 +724,6 @@ include("ensemble/ensemble_analysis.jl") include("solve.jl") include("interpolation.jl") -include("integrand_interface.jl") include("integrator_interface.jl") include("tabletraits.jl") include("remake.jl") @@ -788,7 +797,8 @@ export remake export ODEFunction, DiscreteFunction, ImplicitDiscreteFunction, SplitFunction, DAEFunction, DDEFunction, SDEFunction, SplitSDEFunction, RODEFunction, SDDEFunction, - IncrementingODEFunction, NonlinearFunction, IntervalNonlinearFunction, BVPFunction + IncrementingODEFunction, NonlinearFunction, IntervalNonlinearFunction, BVPFunction, + IntegralFunction, BatchIntegralFunction export OptimizationFunction diff --git a/src/integrand_interface.jl b/src/integrand_interface.jl deleted file mode 100644 index 7b7ebc18e..000000000 --- a/src/integrand_interface.jl +++ /dev/null @@ -1,95 +0,0 @@ -""" - InplaceIntegrand(f!, result::AbstractArray) - -Constructor for a `InplaceIntegrand` accepting an integrand of the form `f!(y,x,p)`. The -caller also provides an output array needed to store the result of the quadrature. -Intermediate `y` arrays are allocated during the calculation, and the final result is is -written to `result`, so use the IntegralSolution immediately after the calculation to read -the result, and don't expect it to persist if the same integrand is used for another -calculation. -""" -struct InplaceIntegrand{F, T <: AbstractArray} - # in-place function f!(y, x, p) that takes one x value and outputs an array of results in-place - f!::F - I::T -end - -""" - BatchIntegrand(f!, y::AbstractArray, x::AbstractVector, max_batch=typemax(Int)) - -Constructor for a `BatchIntegrand` accepting an integrand of the form `f!(y,x,p) = y .= f!.(x, Ref(p))` -that can evaluate the integrand at multiple quadrature nodes using, for example, threads, -the GPU, or distributed-memory. The `max_batch` keyword is a soft limit on the number of -nodes passed to the integrand. The buffers `y,x` must both be `resize!`-able since the -number of evaluation points may vary between calls to `f!`. -""" -struct BatchIntegrand{F, Y, X} - # in-place function f!(y, x, p) that takes an array of x values and outputs an array of results in-place - f!::F - y::Y - x::X - max_batch::Int # maximum number of x to supply in parallel - function BatchIntegrand(f!, - y::AbstractVector, - x::AbstractVector, - max_batch::Integer = typemax(Int)) - max_batch > 0 || throw(ArgumentError("maximum batch size must be positive")) - return new{typeof(f!), typeof(y), typeof(x)}(f!, y, x, max_batch) - end -end - -""" - BatchIntegrand(f!, y, x; max_batch=typemax(Int)) - -Constructor for a `BatchIntegrand` with pre-allocated buffers. -""" -function BatchIntegrand(f!, y, x; max_batch::Integer = typemax(Int)) - BatchIntegrand(f!, y, x, max_batch) -end - -""" - BatchIntegrand(f!, y::Type, x::Type=Nothing; max_batch=typemax(Int)) - -Constructor for a `BatchIntegrand` whose range type is known. The domain type is optional. -Array buffers for those types are allocated internally. -""" -function BatchIntegrand(f!, Y::Type, X::Type = Nothing; max_batch::Integer = typemax(Int)) - BatchIntegrand(f!, Y[], X[], max_batch) -end - -""" - InplaceBatchIntegrand(f!, result::AbstractArray, y::AbstractArray, x::AbstractVector, max_batch=typemax(Int)) - -Constructor for a `InplaceBatchIntegrand` accepting an integrand of the form `f!(y,x,p) = y -.= f!.(x, Ref(p))` that can evaluate an inplace, array-valued integrand at multiple -quadrature nodes simultaneously using, for example, threads, the GPU, or distributed-memory. -The `max_batch` keyword is a soft limit on the number of nodes passed to the integrand. The -buffers `y,x` must both be `resize!`-able since the number of evaluation points may vary -between calls to `f!`. In particular, for a resizeable `y` buffer see ElasticArrays.jl . The -solution is written inplace to `result`. -""" -struct InplaceBatchIntegrand{F, T, Y, X} - # in-place function f!(y, x, p) that takes an array of x values and outputs an array of results in-place - f!::F - I::T - y::Y - x::X - max_batch::Int # maximum number of x to supply in parallel - function InplaceBatchIntegrand(f!, - I::AbstractArray, - y::AbstractArray, - x::AbstractVector, - max_batch::Integer = typemax(Int)) - max_batch > 0 || throw(ArgumentError("maximum batch size must be positive")) - return new{typeof(f!), typeof(I), typeof(y), typeof(x)}(f!, I, y, x, max_batch) - end -end - -""" - InplaceBatchIntegrand(f!, result, y, x; max_batch=typemax(Int)) - -Constructor for a `InplaceBatchIntegrand` with pre-allocated buffers. -""" -function InplaceBatchIntegrand(f!, result, y, x; max_batch::Integer = typemax(Int)) - InplaceBatchIntegrand(f!, result, y, x, max_batch) -end diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index 45d727074..a24cd5238 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -2261,6 +2261,184 @@ end TruncatedStacktraces.@truncate_stacktrace BVPFunction 1 2 +@doc doc""" + IntegralFunction{iip,specialize,F,T,J,TJ,TPJ,Ta,S,JP,SP,TCV,O} <: AbstractIntegralFunction{iip} + +A representation of an integrand `f` defined by: + +```math +f(u, p) +``` + +and its related functions, such as its Jacobian and gradient with respect to parameters. For +an in-place form of `f` see the `iip` section below for details on in-place or out-of-place +handling. + +```julia +IntegralFunction{iip,specialize}(f, [I]; + jac = __has_jac(f) ? f.jac : nothing, + paramjac = __has_paramjac(f) ? f.paramjac : nothing, + analytic = __has_analytic(f) ? f.analytic : nothing, + syms = __has_syms(f) ? f.syms : nothing, + jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing, + sparsity = __has_sparsity(f) ? f.sparsity : nothing, + colorvec = __has_colorvec(f) ? f.colorvec : nothing, + observed = __has_observed(f) ? f.observed : nothing) +``` + +Note that only `f` is required, and in the case of inplace integrands a mutable container +`I` to store the result of the integral. + +The remaining functions are optional and mainly used for accelerating the usage of `f`: +- `jac`: unused +- `paramjac`: unused +- `analytic`: unused +- `syms`: unused +- `jac_prototype`: unused +- `sparsity`: unused +- `colorvec`: unused +- `observed`: unused + +Since most arguments are unused, the following constructor provides the essential behavior: + +```julia +IntegralFunction(f, [I]; kws..) +``` + +If `I` is present, `f` is interpreted as in-place, and otherwise `f` is assumed to be +out-of-place. + +## iip: In-Place vs Out-Of-Place + +Out-of-place functions must be of the form ``f(u, p)`` and in-place functions of the form +``f(y, u, p)``. Since `f` is allowed to return any type (e.g. real or complex numbers or +arrays), in-place functions must provide a container `I` that is of the right type for the +final result of the integral, and the result is written to this container in-place. When +in-place forms are used, in-place array operations may be used by algorithms to reduce +allocations. If `I` is not provided, `f` is assumed to be out-of-place and quadrature is +performed assuming immutable return types. + +## specialize + +This field is currently unused + +## Fields + +The fields of the IntegralFunction type directly match the names of the inputs. +""" +struct IntegralFunction{iip, specialize, F, T, TJ, TPJ, Ta, S, JP, SP, TCV, O} <: + AbstractIntegralFunction{iip} + f::F + I::T + jac::TJ + paramjac::TPJ + analytic::Ta + syms::S + jac_prototype::JP + sparsity::SP + colorvec::TCV + observed::O +end + +TruncatedStacktraces.@truncate_stacktrace IntegralFunction 1 2 + +@doc doc""" +BatchIntegralFunction{iip,specialize,F,T,Y,J,TJ,TPJ,Ta,S,JP,SP,TCV,O} <: AbstractIntegralFunction{iip} + +A representation of an integrand `f` that can be evaluated at multiple points simultaneously +using threads, the gpu, or distributed memory defined by: + +```math +f(y, u, p) +``` + +and its related functions, such as its Jacobian and gradient with respect to parameters. For +an in-place form of `f` see the `iip` section below for details on in-place or out-of-place +handling. + +``u`` is a vector whose elements correspond to distinct evaluation points to `f`, whose +output must be returned in the corresponding entries of ``y``. In general, the integration +algorithm is allowed to vary the number of evaluation points between subsequent calls to `f` + +```julia +BatchIntegralFunction{iip,specialize}(f, y, [I]; + max_batch = typemax(Int), + jac = __has_jac(f) ? f.jac : nothing, + paramjac = __has_paramjac(f) ? f.paramjac : nothing, + analytic = __has_analytic(f) ? f.analytic : nothing, + syms = __has_syms(f) ? f.syms : nothing, + jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing, + sparsity = __has_sparsity(f) ? f.sparsity : nothing, + colorvec = __has_colorvec(f) ? f.colorvec : nothing, + observed = __has_observed(f) ? f.observed : nothing) +``` + +Note that `f` is required and a `resize`-able buffer `y` to store the output, or range of +`f`, and in the case of inplace integrands a mutable container `I` to store the result of +the integral. These buffers can be reused across multiple compatible integrals to reduce +allocations. + +The keyword `max_batch` is used to set a soft limit on the number of points to batch at the +same time so that memory usage is controlled. + +The remaining functions are optional and mainly used for accelerating the usage of `f`: +- `jac`: unused +- `paramjac`: unused +- `analytic`: unused +- `syms`: unused +- `jac_prototype`: unused +- `sparsity`: unused +- `colorvec`: unused +- `observed`: unused + +Since most arguments are unused, the following constructor provides the essential behavior: + +```julia +BatchIntegralFunction(f, y, [I]; max_batch=typemax(Int), kws..) +``` + +If `I` is present, `f` is interpreted as in-place, and otherwise `f` is assumed to be +out-of-place. + +## iip: In-Place vs Out-Of-Place + +Out-of-place and in-place functions are both of the form ``f(y, u, p)``, but differ in the +element type of ``y``. Since `f` is allowed to return any type (e.g. real or complex numbers +or arrays), in-place functions must provide a container `I` that is of the right type for +the final result of the integral, and the result is written to this container in-place. When +`f` is in-place, the output buffer ``y`` is assumed to have a mutable element type, and the +last dimension of ``y`` should correspond to the batch index. For example, ``y`` would have +to be an `ElasticArray` or a `VectorOfSimilarArrays` of an `ElasticArray`. When in-place +forms are used, in-place array operations may be used by algorithms to reduce allocations. +If `I` is not provided, `f` is assumed to be out-of-place and quadrature is performed +assuming ``y`` is an `AbstractVector` with an immutable element type. + +## specialize + +This field is currently unused + +## Fields + +The fields of the BatchIntegralFunction type directly match the names of the inputs. +""" +struct BatchIntegralFunction{iip, specialize, F, Y, T, TJ, TPJ, Ta, S, JP, SP, TCV, O} <: + AbstractIntegralFunction{iip} + f::F + y::Y + I::T + max_batch::Int + jac::TJ + paramjac::TPJ + analytic::Ta + syms::S + jac_prototype::JP + sparsity::SP + colorvec::TCV + observed::O +end + +TruncatedStacktraces.@truncate_stacktrace BatchIntegralFunction 1 2 + ######### Backwards Compatibility Overloads (f::ODEFunction)(args...) = f.f(args...) @@ -3955,6 +4133,80 @@ function BVPFunction(f, bc; kwargs...) end BVPFunction(f::BVPFunction; kwargs...) = f +function IntegralFunction{iip, specialize}(f, I; + jac = __has_jac(f) ? f.jac : nothing, + paramjac = __has_paramjac(f) ? f.paramjac : nothing, + analytic = __has_analytic(f) ? f.analytic : nothing, + syms = __has_syms(f) ? f.syms : nothing, + jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing, + sparsity = __has_sparsity(f) ? f.sparsity : nothing, + colorvec = __has_colorvec(f) ? f.colorvec : nothing, + observed = __has_observed(f) ? f.observed : nothing) where {iip, specialize} + IntegralFunction{ + iip, + specialize, + typeof(f), + typeof(I), + typeof(jac), + typeof(paramjac), + typeof(analytic), + typeof(syms), + typeof(jac_prototype), + typeof(sparsity), + typeof(colorvec), + typeof(observed), + }(f, + I, + jac, + paramjac, + analytic, + syms, + jac_prototype, + sparsity, + colorvec, + observed) +end + +function IntegralFunction{iip}(f, I; kws...) where {iip} + return IntegralFunction{iip, FullSpecialize}(f, I; kws...) +end +IntegralFunction(f; kws...) = IntegralFunction{false}(f, nothing; kws...) +IntegralFunction(f, I; kws...) = IntegralFunction{true}(f, I; kws...) + +function BatchIntegralFunction{iip, specialize}(f, y, I; + max_batch::Integer = typemax(Int), + jac = __has_jac(f) ? f.jac : nothing, + paramjac = __has_paramjac(f) ? f.paramjac : nothing, + analytic = __has_analytic(f) ? f.analytic : nothing, + syms = __has_syms(f) ? f.syms : nothing, + jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing, + sparsity = __has_sparsity(f) ? f.sparsity : nothing, + colorvec = __has_colorvec(f) ? f.colorvec : nothing, + observed = __has_observed(f) ? f.observed : nothing) where {iip, specialize} + BatchIntegralFunction{ + iip, + specialize, + typeof(f), + typeof(y), + typeof(I), + typeof(jac), + typeof(paramjac), + typeof(analytic), + typeof(syms), + typeof(jac_prototype), + typeof(sparsity), + typeof(colorvec), + typeof(observed), + }(f, y, I, max_batch, jac, paramjac, analytic, syms, jac_prototype, sparsity, colorvec, + observed) +end + +function BatchIntegralFunction{iip}(f, y, I; kws...) where {iip} + return BatchIntegralFunction{iip, FullSpecialize}(f, y, I; kws...) +end +BatchIntegralFunction(f, y; kws...) = BatchIntegralFunction{false}(f, y, nothing; kws...) +BatchIntegralFunction(f, y, I; kws...) = BatchIntegralFunction{true}(f, y, I; kws...) + ########## Existence Functions # Check that field/property exists (may be nothing) @@ -4064,7 +4316,9 @@ for S in [:ODEFunction :NonlinearFunction :IntervalNonlinearFunction :IncrementingODEFunction - :BVPFunction] + :BVPFunction + :IntegralFunction + :BatchIntegralFunction] @eval begin function ConstructionBase.constructorof(::Type{<:$S{iip}}) where { iip, diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index c424a6f84..ccbac5bb8 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -601,3 +601,20 @@ BVPFunction(bfoop, bciip, vjp = bvjp) bvjp(du, u, v, p, t) = [1.0] BVPFunction(bfiip, bciip, vjp = bvjp) BVPFunction(bfoop, bciip, vjp = bvjp) + +# IntegralFunction + +ioop(u, p) = p * u +iiip(y, u, p) = y .= u * p + +IntegralFunction(ioop) +IntegralFunction(iiip, Float64[]) + +# BatchIntegralFunction + +boop(y, u, p) = y .= p .* u +biip(y, u, p) = y .= p .* u # this example is not realistic + +BatchIntegralFunction(boop, Float64[]) +BatchIntegralFunction(boop, Float64[], max_batch = 20) +BatchIntegralFunction(biip, Float64[], Float64[]) # the 2nd argument should be an ElasticArray From f56d6545b88f47aa55925f82fc6f31000e8f0718 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Tue, 19 Sep 2023 02:51:33 -0400 Subject: [PATCH 35/70] canonicalize --- src/scimlfunctions.jl | 263 +++++++++++++++++------------------------- 1 file changed, 106 insertions(+), 157 deletions(-) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index a24cd5238..e24533809 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -203,6 +203,48 @@ function Base.showerror(io::IO, e::NonconformingFunctionsError) printstyled(io, e.nonconforming; bold = true, color = :red) end +const INTEGRAND_MISMATCH_FUNCTIONS_ERROR_MESSAGE = """ + Nonconforming functions detected. If an integrand function `f` is defined + as out-of-place (`f(u,p)`), then no integrand_prototype can be passed into the + function constructor. Likewise if `f` is defined as in-place (`f(out,u,p)`), then + an integrand_prototype is required. Either change the use of the function + constructor or define the appropriate dispatch for `f`. + """ + +struct IntegrandMismatchFunctionError <: Exception + iip::Bool + integrand_passed::Bool +end + +function Base.showerror(io::IO, e::IntegrandMismatchFunctionError) + println(io, INTEGRAND_MISMATCH_FUNCTIONS_ERROR_MESSAGE) + print(io, "Mismatch: IIP=") + printstyled(io, e.iip; bold = true, color = :red) + print(io, ", Integrand passed=") + printstyled(io, e.integrand_passed; bold = true, color = :red) +end + +const BATCH_INTEGRAND_MISMATCH_FUNCTIONS_ERROR_MESSAGE = """ + Nonconforming functions detected. If an integrand function `f` is defined + as out-of-place (`f(u,p)`), then no integrand_prototype can be passed into the + function constructor. Likewise if `f` is defined as in-place (`f(out,u,p)`), then + an integrand_prototype is required. Either change the use of the function + constructor or define the appropriate dispatch for `f`. + """ + +struct BatchedIntegrandMismatchFunctionError <: Exception + iip::Bool + integrand_passed::Bool +end + +function Base.showerror(io::IO, e::BatchedIntegrandMismatchFunctionError) + println(io, BATCH_INTEGRAND_MISMATCH_FUNCTIONS_ERROR_MESSAGE) + print(io, "Mismatch: IIP=") + printstyled(io, e.iip; bold = true, color = :red) + print(io, ", Integrand passed=") + printstyled(io, e.integrand_passed; bold = true, color = :red) +end + """ $(TYPEDEF) """ @@ -2262,7 +2304,7 @@ end TruncatedStacktraces.@truncate_stacktrace BVPFunction 1 2 @doc doc""" - IntegralFunction{iip,specialize,F,T,J,TJ,TPJ,Ta,S,JP,SP,TCV,O} <: AbstractIntegralFunction{iip} + IntegralFunction{iip,specialize,F} <: AbstractIntegralFunction{iip} A representation of an integrand `f` defined by: @@ -2270,53 +2312,26 @@ A representation of an integrand `f` defined by: f(u, p) ``` -and its related functions, such as its Jacobian and gradient with respect to parameters. For -an in-place form of `f` see the `iip` section below for details on in-place or out-of-place +For an in-place form of `f` see the `iip` section below for details on in-place or out-of-place handling. ```julia -IntegralFunction{iip,specialize}(f, [I]; - jac = __has_jac(f) ? f.jac : nothing, - paramjac = __has_paramjac(f) ? f.paramjac : nothing, - analytic = __has_analytic(f) ? f.analytic : nothing, - syms = __has_syms(f) ? f.syms : nothing, - jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing, - sparsity = __has_sparsity(f) ? f.sparsity : nothing, - colorvec = __has_colorvec(f) ? f.colorvec : nothing, - observed = __has_observed(f) ? f.observed : nothing) +IntegralFunction{iip,specialize}(f, [integrand_prototype]) ``` Note that only `f` is required, and in the case of inplace integrands a mutable container -`I` to store the result of the integral. - -The remaining functions are optional and mainly used for accelerating the usage of `f`: -- `jac`: unused -- `paramjac`: unused -- `analytic`: unused -- `syms`: unused -- `jac_prototype`: unused -- `sparsity`: unused -- `colorvec`: unused -- `observed`: unused - -Since most arguments are unused, the following constructor provides the essential behavior: - -```julia -IntegralFunction(f, [I]; kws..) -``` - -If `I` is present, `f` is interpreted as in-place, and otherwise `f` is assumed to be -out-of-place. +`integrand_prototype` to store the result of the integral. If `integrand_prototype` is present, +`f` is interpreted as in-place, and otherwise `f` is assumed to be out-of-place. ## iip: In-Place vs Out-Of-Place Out-of-place functions must be of the form ``f(u, p)`` and in-place functions of the form ``f(y, u, p)``. Since `f` is allowed to return any type (e.g. real or complex numbers or -arrays), in-place functions must provide a container `I` that is of the right type for the -final result of the integral, and the result is written to this container in-place. When -in-place forms are used, in-place array operations may be used by algorithms to reduce -allocations. If `I` is not provided, `f` is assumed to be out-of-place and quadrature is -performed assuming immutable return types. +arrays), in-place functions must provide a container `integrand_prototype` that is of the +right type for the final result of the integral, and the result is written to this container +in-place. When in-place forms are used, in-place array operations may be used by algorithms +to reduce allocations. If `integrand_prototype` is not provided, `f` is assumed to be +out-of-place and quadrature is performed assuming immutable return types. ## specialize @@ -2326,18 +2341,10 @@ This field is currently unused The fields of the IntegralFunction type directly match the names of the inputs. """ -struct IntegralFunction{iip, specialize, F, T, TJ, TPJ, Ta, S, JP, SP, TCV, O} <: +struct IntegralFunction{iip, specialize, F, T} <: AbstractIntegralFunction{iip} f::F - I::T - jac::TJ - paramjac::TPJ - analytic::Ta - syms::S - jac_prototype::JP - sparsity::SP - colorvec::TCV - observed::O + integrand_prototype::T end TruncatedStacktraces.@truncate_stacktrace IntegralFunction 1 2 @@ -2352,25 +2359,13 @@ using threads, the gpu, or distributed memory defined by: f(y, u, p) ``` -and its related functions, such as its Jacobian and gradient with respect to parameters. For -an in-place form of `f` see the `iip` section below for details on in-place or out-of-place -handling. - ``u`` is a vector whose elements correspond to distinct evaluation points to `f`, whose output must be returned in the corresponding entries of ``y``. In general, the integration algorithm is allowed to vary the number of evaluation points between subsequent calls to `f` ```julia -BatchIntegralFunction{iip,specialize}(f, y, [I]; - max_batch = typemax(Int), - jac = __has_jac(f) ? f.jac : nothing, - paramjac = __has_paramjac(f) ? f.paramjac : nothing, - analytic = __has_analytic(f) ? f.analytic : nothing, - syms = __has_syms(f) ? f.syms : nothing, - jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing, - sparsity = __has_sparsity(f) ? f.sparsity : nothing, - colorvec = __has_colorvec(f) ? f.colorvec : nothing, - observed = __has_observed(f) ? f.observed : nothing) +BatchIntegralFunction{iip,specialize}(f, y, [integrand_prototype]; + max_batch=typemax(Int)) ``` Note that `f` is required and a `resize`-able buffer `y` to store the output, or range of @@ -2381,37 +2376,27 @@ allocations. The keyword `max_batch` is used to set a soft limit on the number of points to batch at the same time so that memory usage is controlled. -The remaining functions are optional and mainly used for accelerating the usage of `f`: -- `jac`: unused -- `paramjac`: unused -- `analytic`: unused -- `syms`: unused -- `jac_prototype`: unused -- `sparsity`: unused -- `colorvec`: unused -- `observed`: unused - Since most arguments are unused, the following constructor provides the essential behavior: ```julia -BatchIntegralFunction(f, y, [I]; max_batch=typemax(Int), kws..) +BatchIntegralFunction(f, y, [integrand_prototype]; max_batch=typemax(Int), kws..) ``` -If `I` is present, `f` is interpreted as in-place, and otherwise `f` is assumed to be +If `integrand_prototype` is present, `f` is interpreted as in-place, and otherwise `f` is assumed to be out-of-place. ## iip: In-Place vs Out-Of-Place Out-of-place and in-place functions are both of the form ``f(y, u, p)``, but differ in the element type of ``y``. Since `f` is allowed to return any type (e.g. real or complex numbers -or arrays), in-place functions must provide a container `I` that is of the right type for -the final result of the integral, and the result is written to this container in-place. When -`f` is in-place, the output buffer ``y`` is assumed to have a mutable element type, and the -last dimension of ``y`` should correspond to the batch index. For example, ``y`` would have -to be an `ElasticArray` or a `VectorOfSimilarArrays` of an `ElasticArray`. When in-place -forms are used, in-place array operations may be used by algorithms to reduce allocations. -If `I` is not provided, `f` is assumed to be out-of-place and quadrature is performed -assuming ``y`` is an `AbstractVector` with an immutable element type. +or arrays), in-place functions must provide a container `integrand_prototype` that is of the +right type for the final result of the integral, and the result is written to this container +in-place. When `f` is in-place, the output buffer ``y`` is assumed to have a mutable element +type, and the last dimension of ``y`` should correspond to the batch index. For example, +``y`` would have to be an `ElasticArray` or a `VectorOfSimilarArrays` of an `ElasticArray`. +When in-place forms are used, in-place array operations may be used by algorithms to reduce +allocations. If `integrand_prototype` is not provided, `f` is assumed to be out-of-place and +quadrature is performed assuming ``y`` is an `AbstractVector` with an immutable element type. ## specialize @@ -2421,20 +2406,12 @@ This field is currently unused The fields of the BatchIntegralFunction type directly match the names of the inputs. """ -struct BatchIntegralFunction{iip, specialize, F, Y, T, TJ, TPJ, Ta, S, JP, SP, TCV, O} <: +struct BatchIntegralFunction{iip, specialize, F, Y, T} <: AbstractIntegralFunction{iip} f::F y::Y - I::T + integrand_prototype::T max_batch::Int - jac::TJ - paramjac::TPJ - analytic::Ta - syms::S - jac_prototype::JP - sparsity::SP - colorvec::TCV - observed::O end TruncatedStacktraces.@truncate_stacktrace BatchIntegralFunction 1 2 @@ -4133,79 +4110,51 @@ function BVPFunction(f, bc; kwargs...) end BVPFunction(f::BVPFunction; kwargs...) = f -function IntegralFunction{iip, specialize}(f, I; - jac = __has_jac(f) ? f.jac : nothing, - paramjac = __has_paramjac(f) ? f.paramjac : nothing, - analytic = __has_analytic(f) ? f.analytic : nothing, - syms = __has_syms(f) ? f.syms : nothing, - jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing, - sparsity = __has_sparsity(f) ? f.sparsity : nothing, - colorvec = __has_colorvec(f) ? f.colorvec : nothing, - observed = __has_observed(f) ? f.observed : nothing) where {iip, specialize} - IntegralFunction{ - iip, - specialize, - typeof(f), - typeof(I), - typeof(jac), - typeof(paramjac), - typeof(analytic), - typeof(syms), - typeof(jac_prototype), - typeof(sparsity), - typeof(colorvec), - typeof(observed), - }(f, - I, - jac, - paramjac, - analytic, - syms, - jac_prototype, - sparsity, - colorvec, - observed) +function IntegralFunction{iip, specialize}(f, integrand_prototype) where {iip, specialize} + IntegralFunction{iip,specialize,typeof(f),typeof(I)}(f,integrand_prototype) end -function IntegralFunction{iip}(f, I; kws...) where {iip} - return IntegralFunction{iip, FullSpecialize}(f, I; kws...) +function IntegralFunction{iip}(f, integrand_prototype) where {iip} + return IntegralFunction{iip, FullSpecialize}(f, integrand_prototype) +end +function IntegralFunction(f) + calcuated_iip = isinplace(f, 3, "integral", iip) + if !calcuated_iip + throw(IntegrandMismatchFunctionError(calculated_iip, false)) + end + IntegralFunction{false}(f, nothing) +end +function IntegralFunction(f, integrand_prototype) + calcuated_iip = isinplace(f, 3, "integral", iip) + if !calcuated_iip + throw(IntegrandMismatchFunctionError(calculated_iip, true)) + end + IntegralFunction{true}(f, integrand_prototype) end -IntegralFunction(f; kws...) = IntegralFunction{false}(f, nothing; kws...) -IntegralFunction(f, I; kws...) = IntegralFunction{true}(f, I; kws...) -function BatchIntegralFunction{iip, specialize}(f, y, I; - max_batch::Integer = typemax(Int), - jac = __has_jac(f) ? f.jac : nothing, - paramjac = __has_paramjac(f) ? f.paramjac : nothing, - analytic = __has_analytic(f) ? f.analytic : nothing, - syms = __has_syms(f) ? f.syms : nothing, - jac_prototype = __has_jac_prototype(f) ? f.jac_prototype : nothing, - sparsity = __has_sparsity(f) ? f.sparsity : nothing, - colorvec = __has_colorvec(f) ? f.colorvec : nothing, - observed = __has_observed(f) ? f.observed : nothing) where {iip, specialize} - BatchIntegralFunction{ - iip, - specialize, - typeof(f), - typeof(y), - typeof(I), - typeof(jac), - typeof(paramjac), - typeof(analytic), - typeof(syms), - typeof(jac_prototype), - typeof(sparsity), - typeof(colorvec), - typeof(observed), - }(f, y, I, max_batch, jac, paramjac, analytic, syms, jac_prototype, sparsity, colorvec, - observed) +function BatchIntegralFunction{iip, specialize}(f, output_prototype, integrand_prototype; + max_batch::Integer = typemax(Int)) where {iip, specialize} + BatchIntegralFunction{iip,specialize,typeof(f),typeof(output_prototype),typeof(integrand_prototype) + }(f, output_prototype, integrand_prototype, max_batch) end -function BatchIntegralFunction{iip}(f, y, I; kws...) where {iip} - return BatchIntegralFunction{iip, FullSpecialize}(f, y, I; kws...) +function BatchIntegralFunction{iip}(f, output_prototype, integrand_prototype; kwargs...) where {iip} + return BatchIntegralFunction{iip, FullSpecialize}(f, output_prototype, integrand_prototype; kwargs...) +end +function BatchIntegralFunction(f, output_prototype; kwargs...) + calcuated_iip = isinplace(f, 4, "batchintegral", iip) + if !calcuated_iip + throw(BatchedIntegrandMismatchFunctionError(calculated_iip, false)) + end + BatchIntegralFunction{false}(f, output_prototype, nothing; kwargs...) +end +function BatchIntegralFunction(f, output_prototype, integrand_prototype; kwargs...) + calcuated_iip = isinplace(f, 4, "batchintegral", iip) + if !calcuated_iip + throw(BatchIntegrandMismatchFunctionError(calculated_iip, true)) + end + BatchIntegralFunction{true}(f, output_prototype, integrand_prototype; kwargs...) end -BatchIntegralFunction(f, y; kws...) = BatchIntegralFunction{false}(f, y, nothing; kws...) -BatchIntegralFunction(f, y, I; kws...) = BatchIntegralFunction{true}(f, y, I; kws...) ########## Existence Functions From e3ee4539316a8fac403053dbcad40a2a840dd6cf Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Tue, 19 Sep 2023 02:52:31 -0400 Subject: [PATCH 36/70] Remove error checking on function definition of batch integral --- src/scimlfunctions.jl | 31 ++----------------------------- 1 file changed, 2 insertions(+), 29 deletions(-) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index e24533809..d874c04e2 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -224,27 +224,6 @@ function Base.showerror(io::IO, e::IntegrandMismatchFunctionError) printstyled(io, e.integrand_passed; bold = true, color = :red) end -const BATCH_INTEGRAND_MISMATCH_FUNCTIONS_ERROR_MESSAGE = """ - Nonconforming functions detected. If an integrand function `f` is defined - as out-of-place (`f(u,p)`), then no integrand_prototype can be passed into the - function constructor. Likewise if `f` is defined as in-place (`f(out,u,p)`), then - an integrand_prototype is required. Either change the use of the function - constructor or define the appropriate dispatch for `f`. - """ - -struct BatchedIntegrandMismatchFunctionError <: Exception - iip::Bool - integrand_passed::Bool -end - -function Base.showerror(io::IO, e::BatchedIntegrandMismatchFunctionError) - println(io, BATCH_INTEGRAND_MISMATCH_FUNCTIONS_ERROR_MESSAGE) - print(io, "Mismatch: IIP=") - printstyled(io, e.iip; bold = true, color = :red) - print(io, ", Integrand passed=") - printstyled(io, e.integrand_passed; bold = true, color = :red) -end - """ $(TYPEDEF) """ @@ -4142,17 +4121,11 @@ function BatchIntegralFunction{iip}(f, output_prototype, integrand_prototype; kw return BatchIntegralFunction{iip, FullSpecialize}(f, output_prototype, integrand_prototype; kwargs...) end function BatchIntegralFunction(f, output_prototype; kwargs...) - calcuated_iip = isinplace(f, 4, "batchintegral", iip) - if !calcuated_iip - throw(BatchedIntegrandMismatchFunctionError(calculated_iip, false)) - end + calcuated_iip = isinplace(f, 3, "batchintegral", iip) BatchIntegralFunction{false}(f, output_prototype, nothing; kwargs...) end function BatchIntegralFunction(f, output_prototype, integrand_prototype; kwargs...) - calcuated_iip = isinplace(f, 4, "batchintegral", iip) - if !calcuated_iip - throw(BatchIntegrandMismatchFunctionError(calculated_iip, true)) - end + calcuated_iip = isinplace(f, 3, "batchintegral", iip) BatchIntegralFunction{true}(f, output_prototype, integrand_prototype; kwargs...) end From 318e79f8106b38602fe290a0fd3d93f8539eadd0 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Tue, 19 Sep 2023 02:55:09 -0400 Subject: [PATCH 37/70] add error test on incorrect integral function dispatches --- test/function_building_error_messages.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index ccbac5bb8..bf68b0fca 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -610,6 +610,9 @@ iiip(y, u, p) = y .= u * p IntegralFunction(ioop) IntegralFunction(iiip, Float64[]) +@test_throws SciMLBase.IntegrandMismatchFunctionError IntegralFunction(ioop, Float64[]) +@test_throws SciMLBase.IntegrandMismatchFunctionError IntegralFunction(iiip) + # BatchIntegralFunction boop(y, u, p) = y .= p .* u From 783b88e768058387793859ed13032fb63945b7ea Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Tue, 19 Sep 2023 02:59:40 -0400 Subject: [PATCH 38/70] argument amounts testing --- test/function_building_error_messages.jl | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index bf68b0fca..41b31bcff 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -606,18 +606,28 @@ BVPFunction(bfoop, bciip, vjp = bvjp) ioop(u, p) = p * u iiip(y, u, p) = y .= u * p +i1(u) = u +itoo(y, u, p, a) = y .= u * p IntegralFunction(ioop) IntegralFunction(iiip, Float64[]) @test_throws SciMLBase.IntegrandMismatchFunctionError IntegralFunction(ioop, Float64[]) -@test_throws SciMLBase.IntegrandMismatchFunctionError IntegralFunction(iiip) +@test_throws SciMLBase.IntegrandMismatchFunctionError IntegralFunction(i1) +@test_throws SciMLBase.TooFewArgumentsError IntegralFunction(i1) +@test_throws SciMLBase.TooManyArgumentsError IntegralFunction(itoo) +@test_throws SciMLBase.TooManyArgumentsError IntegralFunction(itoo, Float64[]) # BatchIntegralFunction boop(y, u, p) = y .= p .* u biip(y, u, p) = y .= p .* u # this example is not realistic +bi1(y, u) = y .= p .* u +bitoo(y, u, p, a) = y .= p .* u BatchIntegralFunction(boop, Float64[]) BatchIntegralFunction(boop, Float64[], max_batch = 20) BatchIntegralFunction(biip, Float64[], Float64[]) # the 2nd argument should be an ElasticArray +@test_throws SciMLBase.TooFewArgumentsError IntegralFunction(bi1) +@test_throws SciMLBase.TooManyArgumentsError IntegralFunction(bitoo) +@test_throws SciMLBase.TooManyArgumentsError IntegralFunction(bitoo, Float64[]) \ No newline at end of file From 819ec920ac894040c0687de94bbcc97b419a98d4 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 19 Sep 2023 03:08:44 -0400 Subject: [PATCH 39/70] Update src/problems/basic_problems.jl --- src/problems/basic_problems.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 1dc56ec9f..efa1bf3a0 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -437,7 +437,7 @@ struct SampledIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false} @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(dim)}, typeof(kwargs)}(y, x, dim, kwargs) + new{typeof(y), typeof(x), typeof(dim), typeof(kwargs)}(y, x, dim, kwargs) end end From fcc7edb1bd1d7418fd23f992878dcd5180ac670a Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Tue, 19 Sep 2023 05:41:52 -0400 Subject: [PATCH 40/70] some better utils checks --- src/scimlfunctions.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index d874c04e2..126c263c1 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -4090,21 +4090,21 @@ end BVPFunction(f::BVPFunction; kwargs...) = f function IntegralFunction{iip, specialize}(f, integrand_prototype) where {iip, specialize} - IntegralFunction{iip,specialize,typeof(f),typeof(I)}(f,integrand_prototype) + IntegralFunction{iip,specialize,typeof(f),typeof(integrand_prototype)}(f,integrand_prototype) end function IntegralFunction{iip}(f, integrand_prototype) where {iip} return IntegralFunction{iip, FullSpecialize}(f, integrand_prototype) end function IntegralFunction(f) - calcuated_iip = isinplace(f, 3, "integral", iip) + calcuated_iip = isinplace(f, 3, "integral", true) if !calcuated_iip throw(IntegrandMismatchFunctionError(calculated_iip, false)) end IntegralFunction{false}(f, nothing) end function IntegralFunction(f, integrand_prototype) - calcuated_iip = isinplace(f, 3, "integral", iip) + calcuated_iip = isinplace(f, 3, "integral", true) if !calcuated_iip throw(IntegrandMismatchFunctionError(calculated_iip, true)) end @@ -4121,11 +4121,11 @@ function BatchIntegralFunction{iip}(f, output_prototype, integrand_prototype; kw return BatchIntegralFunction{iip, FullSpecialize}(f, output_prototype, integrand_prototype; kwargs...) end function BatchIntegralFunction(f, output_prototype; kwargs...) - calcuated_iip = isinplace(f, 3, "batchintegral", iip) + calcuated_iip = isinplace(f, 3, "batchintegral", true; has_two_dispatches = false) BatchIntegralFunction{false}(f, output_prototype, nothing; kwargs...) end function BatchIntegralFunction(f, output_prototype, integrand_prototype; kwargs...) - calcuated_iip = isinplace(f, 3, "batchintegral", iip) + calcuated_iip = isinplace(f, 3, "batchintegral", true; has_two_dispatches = false) BatchIntegralFunction{true}(f, output_prototype, integrand_prototype; kwargs...) end From 5a37040f87b27080967840b1f992ebb6aab636a2 Mon Sep 17 00:00:00 2001 From: lxvm Date: Tue, 19 Sep 2023 06:59:23 -0400 Subject: [PATCH 41/70] apply format --- src/problems/basic_problems.jl | 2 +- src/scimlfunctions.jl | 57 +++++++++++++++--------- test/function_building_error_messages.jl | 4 +- 3 files changed, 39 insertions(+), 24 deletions(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index e0385140d..e3efddcd9 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -434,7 +434,7 @@ struct SampledIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false} @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 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), Val{dim}, typeof(kwargs)}(y, x, Val(dim), kwargs) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index 126c263c1..aaa2f10c6 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -206,7 +206,7 @@ end const INTEGRAND_MISMATCH_FUNCTIONS_ERROR_MESSAGE = """ Nonconforming functions detected. If an integrand function `f` is defined as out-of-place (`f(u,p)`), then no integrand_prototype can be passed into the - function constructor. Likewise if `f` is defined as in-place (`f(out,u,p)`), then + function constructor. Likewise if `f` is defined as in-place (`f(out,u,p)`), then an integrand_prototype is required. Either change the use of the function constructor or define the appropriate dispatch for `f`. """ @@ -2299,17 +2299,17 @@ IntegralFunction{iip,specialize}(f, [integrand_prototype]) ``` Note that only `f` is required, and in the case of inplace integrands a mutable container -`integrand_prototype` to store the result of the integral. If `integrand_prototype` is present, +`integrand_prototype` to store the result of the integral. If `integrand_prototype` is present, `f` is interpreted as in-place, and otherwise `f` is assumed to be out-of-place. ## iip: In-Place vs Out-Of-Place Out-of-place functions must be of the form ``f(u, p)`` and in-place functions of the form ``f(y, u, p)``. Since `f` is allowed to return any type (e.g. real or complex numbers or -arrays), in-place functions must provide a container `integrand_prototype` that is of the -right type for the final result of the integral, and the result is written to this container -in-place. When in-place forms are used, in-place array operations may be used by algorithms -to reduce allocations. If `integrand_prototype` is not provided, `f` is assumed to be +arrays), in-place functions must provide a container `integrand_prototype` that is of the +right type for the final result of the integral, and the result is written to this container +in-place. When in-place forms are used, in-place array operations may be used by algorithms +to reduce allocations. If `integrand_prototype` is not provided, `f` is assumed to be out-of-place and quadrature is performed assuming immutable return types. ## specialize @@ -2368,13 +2368,13 @@ out-of-place. Out-of-place and in-place functions are both of the form ``f(y, u, p)``, but differ in the element type of ``y``. Since `f` is allowed to return any type (e.g. real or complex numbers -or arrays), in-place functions must provide a container `integrand_prototype` that is of the -right type for the final result of the integral, and the result is written to this container -in-place. When `f` is in-place, the output buffer ``y`` is assumed to have a mutable element -type, and the last dimension of ``y`` should correspond to the batch index. For example, -``y`` would have to be an `ElasticArray` or a `VectorOfSimilarArrays` of an `ElasticArray`. -When in-place forms are used, in-place array operations may be used by algorithms to reduce -allocations. If `integrand_prototype` is not provided, `f` is assumed to be out-of-place and +or arrays), in-place functions must provide a container `integrand_prototype` that is of the +right type for the final result of the integral, and the result is written to this container +in-place. When `f` is in-place, the output buffer ``y`` is assumed to have a mutable element +type, and the last dimension of ``y`` should correspond to the batch index. For example, +``y`` would have to be an `ElasticArray` or a `VectorOfSimilarArrays` of an `ElasticArray`. +When in-place forms are used, in-place array operations may be used by algorithms to reduce +allocations. If `integrand_prototype` is not provided, `f` is assumed to be out-of-place and quadrature is performed assuming ``y`` is an `AbstractVector` with an immutable element type. ## specialize @@ -4090,13 +4090,14 @@ end BVPFunction(f::BVPFunction; kwargs...) = f function IntegralFunction{iip, specialize}(f, integrand_prototype) where {iip, specialize} - IntegralFunction{iip,specialize,typeof(f),typeof(integrand_prototype)}(f,integrand_prototype) + IntegralFunction{iip, specialize, typeof(f), typeof(integrand_prototype)}(f, + integrand_prototype) end function IntegralFunction{iip}(f, integrand_prototype) where {iip} return IntegralFunction{iip, FullSpecialize}(f, integrand_prototype) end -function IntegralFunction(f) +function IntegralFunction(f) calcuated_iip = isinplace(f, 3, "integral", true) if !calcuated_iip throw(IntegrandMismatchFunctionError(calculated_iip, false)) @@ -4113,18 +4114,32 @@ end function BatchIntegralFunction{iip, specialize}(f, output_prototype, integrand_prototype; max_batch::Integer = typemax(Int)) where {iip, specialize} - BatchIntegralFunction{iip,specialize,typeof(f),typeof(output_prototype),typeof(integrand_prototype) - }(f, output_prototype, integrand_prototype, max_batch) + BatchIntegralFunction{ + iip, + specialize, + typeof(f), + typeof(output_prototype), + typeof(integrand_prototype), + }(f, + output_prototype, + integrand_prototype, + max_batch) end -function BatchIntegralFunction{iip}(f, output_prototype, integrand_prototype; kwargs...) where {iip} - return BatchIntegralFunction{iip, FullSpecialize}(f, output_prototype, integrand_prototype; kwargs...) +function BatchIntegralFunction{iip}(f, + output_prototype, + integrand_prototype; + kwargs...) where {iip} + return BatchIntegralFunction{iip, FullSpecialize}(f, + output_prototype, + integrand_prototype; + kwargs...) end -function BatchIntegralFunction(f, output_prototype; kwargs...) +function BatchIntegralFunction(f, output_prototype; kwargs...) calcuated_iip = isinplace(f, 3, "batchintegral", true; has_two_dispatches = false) BatchIntegralFunction{false}(f, output_prototype, nothing; kwargs...) end -function BatchIntegralFunction(f, output_prototype, integrand_prototype; kwargs...) +function BatchIntegralFunction(f, output_prototype, integrand_prototype; kwargs...) calcuated_iip = isinplace(f, 3, "batchintegral", true; has_two_dispatches = false) BatchIntegralFunction{true}(f, output_prototype, integrand_prototype; kwargs...) end diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index 41b31bcff..755617583 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -463,7 +463,7 @@ IntegralProblem(intf, [0.0], [1.0], p) x = [1.0, 2.0] y = rand(2, 2) SampledIntegralProblem(y, x) -SampledIntegralProblem(y, x; dim=2) +SampledIntegralProblem(y, x; dim = 2) # Optimization @@ -630,4 +630,4 @@ BatchIntegralFunction(boop, Float64[], max_batch = 20) BatchIntegralFunction(biip, Float64[], Float64[]) # the 2nd argument should be an ElasticArray @test_throws SciMLBase.TooFewArgumentsError IntegralFunction(bi1) @test_throws SciMLBase.TooManyArgumentsError IntegralFunction(bitoo) -@test_throws SciMLBase.TooManyArgumentsError IntegralFunction(bitoo, Float64[]) \ No newline at end of file +@test_throws SciMLBase.TooManyArgumentsError IntegralFunction(bitoo, Float64[]) From b02e47058905e83e0c728c59638828c690b8c093 Mon Sep 17 00:00:00 2001 From: lxvm Date: Tue, 19 Sep 2023 07:03:32 -0400 Subject: [PATCH 42/70] fix integralfunction iip --- src/scimlfunctions.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index aaa2f10c6..da514ac29 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -4098,8 +4098,8 @@ function IntegralFunction{iip}(f, integrand_prototype) where {iip} return IntegralFunction{iip, FullSpecialize}(f, integrand_prototype) end function IntegralFunction(f) - calcuated_iip = isinplace(f, 3, "integral", true) - if !calcuated_iip + calculated_iip = isinplace(f, 3, "integral", true) + if calculated_iip throw(IntegrandMismatchFunctionError(calculated_iip, false)) end IntegralFunction{false}(f, nothing) From 95bdb1d2cc26415154f742a8f01337275d36bc20 Mon Sep 17 00:00:00 2001 From: lxvm Date: Tue, 19 Sep 2023 07:12:27 -0400 Subject: [PATCH 43/70] rename integrand_prototype to integral_prototype --- src/scimlfunctions.jl | 84 +++++++++++++++++++++---------------------- 1 file changed, 40 insertions(+), 44 deletions(-) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index da514ac29..43f9c8b6f 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -205,9 +205,9 @@ end const INTEGRAND_MISMATCH_FUNCTIONS_ERROR_MESSAGE = """ Nonconforming functions detected. If an integrand function `f` is defined - as out-of-place (`f(u,p)`), then no integrand_prototype can be passed into the + as out-of-place (`f(u,p)`), then no integral_prototype can be passed into the function constructor. Likewise if `f` is defined as in-place (`f(out,u,p)`), then - an integrand_prototype is required. Either change the use of the function + an integral_prototype is required. Either change the use of the function constructor or define the appropriate dispatch for `f`. """ @@ -2295,21 +2295,21 @@ For an in-place form of `f` see the `iip` section below for details on in-place handling. ```julia -IntegralFunction{iip,specialize}(f, [integrand_prototype]) +IntegralFunction{iip,specialize}(f, [integral_prototype]) ``` Note that only `f` is required, and in the case of inplace integrands a mutable container -`integrand_prototype` to store the result of the integral. If `integrand_prototype` is present, +`integral_prototype` to store the result of the integral. If `integral_prototype` is present, `f` is interpreted as in-place, and otherwise `f` is assumed to be out-of-place. ## iip: In-Place vs Out-Of-Place Out-of-place functions must be of the form ``f(u, p)`` and in-place functions of the form ``f(y, u, p)``. Since `f` is allowed to return any type (e.g. real or complex numbers or -arrays), in-place functions must provide a container `integrand_prototype` that is of the +arrays), in-place functions must provide a container `integral_prototype` that is of the right type for the final result of the integral, and the result is written to this container in-place. When in-place forms are used, in-place array operations may be used by algorithms -to reduce allocations. If `integrand_prototype` is not provided, `f` is assumed to be +to reduce allocations. If `integral_prototype` is not provided, `f` is assumed to be out-of-place and quadrature is performed assuming immutable return types. ## specialize @@ -2323,13 +2323,14 @@ The fields of the IntegralFunction type directly match the names of the inputs. struct IntegralFunction{iip, specialize, F, T} <: AbstractIntegralFunction{iip} f::F - integrand_prototype::T + integral_prototype::T end TruncatedStacktraces.@truncate_stacktrace IntegralFunction 1 2 @doc doc""" -BatchIntegralFunction{iip,specialize,F,T,Y,J,TJ,TPJ,Ta,S,JP,SP,TCV,O} <: AbstractIntegralFunction{iip} +BatchIntegralFunction{iip,specialize,F,T,Y,J,TJ,TPJ,Ta,S,JP,SP,TCV,O} <: +AbstractIntegralFunction{iip} A representation of an integrand `f` that can be evaluated at multiple points simultaneously using threads, the gpu, or distributed memory defined by: @@ -2343,39 +2344,34 @@ output must be returned in the corresponding entries of ``y``. In general, the i algorithm is allowed to vary the number of evaluation points between subsequent calls to `f` ```julia -BatchIntegralFunction{iip,specialize}(f, y, [integrand_prototype]; +BatchIntegralFunction{iip,specialize}(f, output_prototype, [integral_prototype]; max_batch=typemax(Int)) ``` -Note that `f` is required and a `resize`-able buffer `y` to store the output, or range of -`f`, and in the case of inplace integrands a mutable container `I` to store the result of -the integral. These buffers can be reused across multiple compatible integrals to reduce -allocations. +Note that `f` is required and a `resize`-able buffer `output_prototype` to store the output, +or range of `f`, and in the case of inplace integrands a mutable container +`integral_prototype` to store the result of the integral. These buffers can be reused across +multiple compatible integrals to reduce allocations. The keyword `max_batch` is used to set a soft limit on the number of points to batch at the same time so that memory usage is controlled. -Since most arguments are unused, the following constructor provides the essential behavior: - -```julia -BatchIntegralFunction(f, y, [integrand_prototype]; max_batch=typemax(Int), kws..) -``` - -If `integrand_prototype` is present, `f` is interpreted as in-place, and otherwise `f` is assumed to be -out-of-place. +If `integral_prototype` is present, `f` is interpreted as in-place, and otherwise `f` is +assumed to be out-of-place. ## iip: In-Place vs Out-Of-Place Out-of-place and in-place functions are both of the form ``f(y, u, p)``, but differ in the element type of ``y``. Since `f` is allowed to return any type (e.g. real or complex numbers -or arrays), in-place functions must provide a container `integrand_prototype` that is of the +or arrays), in-place functions must provide a container `integral_prototype` that is of the right type for the final result of the integral, and the result is written to this container -in-place. When `f` is in-place, the output buffer ``y`` is assumed to have a mutable element -type, and the last dimension of ``y`` should correspond to the batch index. For example, -``y`` would have to be an `ElasticArray` or a `VectorOfSimilarArrays` of an `ElasticArray`. -When in-place forms are used, in-place array operations may be used by algorithms to reduce -allocations. If `integrand_prototype` is not provided, `f` is assumed to be out-of-place and -quadrature is performed assuming ``y`` is an `AbstractVector` with an immutable element type. +in-place. When `f` is in-place, the buffer `output_prototype` is assumed to have a mutable +element type, and the last dimension of `output_prototype` should correspond to the batch +index. For example, `output_prototype` would have to be an `ElasticArray` or a +`VectorOfSimilarArrays` of an `ElasticArray`. When in-place forms are used, in-place array +operations may be used by algorithms to reduce allocations. If `integral_prototype` is not +provided, `f` is assumed to be out-of-place and quadrature is performed assuming +`output_prototype` is an `AbstractVector` with an immutable element type. ## specialize @@ -2388,8 +2384,8 @@ The fields of the BatchIntegralFunction type directly match the names of the inp struct BatchIntegralFunction{iip, specialize, F, Y, T} <: AbstractIntegralFunction{iip} f::F - y::Y - integrand_prototype::T + output_prototype::Y + integral_prototype::T max_batch::Int end @@ -4089,13 +4085,13 @@ function BVPFunction(f, bc; kwargs...) end BVPFunction(f::BVPFunction; kwargs...) = f -function IntegralFunction{iip, specialize}(f, integrand_prototype) where {iip, specialize} - IntegralFunction{iip, specialize, typeof(f), typeof(integrand_prototype)}(f, - integrand_prototype) +function IntegralFunction{iip, specialize}(f, integral_prototype) where {iip, specialize} + IntegralFunction{iip, specialize, typeof(f), typeof(integral_prototype)}(f, + integral_prototype) end -function IntegralFunction{iip}(f, integrand_prototype) where {iip} - return IntegralFunction{iip, FullSpecialize}(f, integrand_prototype) +function IntegralFunction{iip}(f, integral_prototype) where {iip} + return IntegralFunction{iip, FullSpecialize}(f, integral_prototype) end function IntegralFunction(f) calculated_iip = isinplace(f, 3, "integral", true) @@ -4104,44 +4100,44 @@ function IntegralFunction(f) end IntegralFunction{false}(f, nothing) end -function IntegralFunction(f, integrand_prototype) +function IntegralFunction(f, integral_prototype) calcuated_iip = isinplace(f, 3, "integral", true) if !calcuated_iip throw(IntegrandMismatchFunctionError(calculated_iip, true)) end - IntegralFunction{true}(f, integrand_prototype) + IntegralFunction{true}(f, integral_prototype) end -function BatchIntegralFunction{iip, specialize}(f, output_prototype, integrand_prototype; +function BatchIntegralFunction{iip, specialize}(f, output_prototype, integral_prototype; max_batch::Integer = typemax(Int)) where {iip, specialize} BatchIntegralFunction{ iip, specialize, typeof(f), typeof(output_prototype), - typeof(integrand_prototype), + typeof(integral_prototype), }(f, output_prototype, - integrand_prototype, + integral_prototype, max_batch) end function BatchIntegralFunction{iip}(f, output_prototype, - integrand_prototype; + integral_prototype; kwargs...) where {iip} return BatchIntegralFunction{iip, FullSpecialize}(f, output_prototype, - integrand_prototype; + integral_prototype; kwargs...) end function BatchIntegralFunction(f, output_prototype; kwargs...) calcuated_iip = isinplace(f, 3, "batchintegral", true; has_two_dispatches = false) BatchIntegralFunction{false}(f, output_prototype, nothing; kwargs...) end -function BatchIntegralFunction(f, output_prototype, integrand_prototype; kwargs...) +function BatchIntegralFunction(f, output_prototype, integral_prototype; kwargs...) calcuated_iip = isinplace(f, 3, "batchintegral", true; has_two_dispatches = false) - BatchIntegralFunction{true}(f, output_prototype, integrand_prototype; kwargs...) + BatchIntegralFunction{true}(f, output_prototype, integral_prototype; kwargs...) end ########## Existence Functions From d70980ed563d985cc145923dc9249c11fc792190 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 20 Sep 2023 16:59:10 -0400 Subject: [PATCH 44/70] Update basic_problems.jl --- src/problems/basic_problems.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index efa1bf3a0..cd057eba5 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -426,10 +426,10 @@ SampledIntegralProblem(y::AbstractArray, x::AbstractVector; dim=ndims(y), kwargs The fields match the names of the constructor arguments. """ -struct SampledIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false} +struct SampledIntegralProblem{Y, X, K} <: AbstractIntegralProblem{false} y::Y x::X - dim::D + dim::Int kwargs::K @add_kwonly function SampledIntegralProblem(y::AbstractArray, x::AbstractVector; dim = ndims(y), @@ -437,7 +437,7 @@ struct SampledIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false} @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(dim), typeof(kwargs)}(y, x, dim, kwargs) + new{typeof(y), typeof(x), typeof(kwargs)}(y, x, dim, kwargs) end end From eb9299515d9c8f50d08dc6406f2a929eb23d8a06 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 20 Sep 2023 16:59:47 -0400 Subject: [PATCH 45/70] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index fa982a687..95a458cb1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SciMLBase" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" authors = ["Chris Rackauckas and contributors"] -version = "1.98.0" +version = "1.98.1" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 5675e6ffd76169b59355e77e4d4d634ef43ec68b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 20 Sep 2023 20:15:57 -0400 Subject: [PATCH 46/70] Update test/function_building_error_messages.jl --- test/function_building_error_messages.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index 755617583..bdf0308c1 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -613,7 +613,7 @@ IntegralFunction(ioop) IntegralFunction(iiip, Float64[]) @test_throws SciMLBase.IntegrandMismatchFunctionError IntegralFunction(ioop, Float64[]) -@test_throws SciMLBase.IntegrandMismatchFunctionError IntegralFunction(i1) +@test_throws SciMLBase.IntegrandMismatchFunctionError IntegralFunction(iiip) @test_throws SciMLBase.TooFewArgumentsError IntegralFunction(i1) @test_throws SciMLBase.TooManyArgumentsError IntegralFunction(itoo) @test_throws SciMLBase.TooManyArgumentsError IntegralFunction(itoo, Float64[]) From 774e4bed8c0a5c6a43f7ddfb1c0f12fde1583cb9 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Wed, 20 Sep 2023 23:15:41 -0400 Subject: [PATCH 47/70] fix typos --- src/scimlfunctions.jl | 2 +- test/function_building_error_messages.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index 43f9c8b6f..30ab6d4aa 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -4103,7 +4103,7 @@ end function IntegralFunction(f, integral_prototype) calcuated_iip = isinplace(f, 3, "integral", true) if !calcuated_iip - throw(IntegrandMismatchFunctionError(calculated_iip, true)) + throw(IntegrandMismatchFunctionError(calcuated_iip, true)) end IntegralFunction{true}(f, integral_prototype) end diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index bdf0308c1..cd47a6a33 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -628,6 +628,6 @@ bitoo(y, u, p, a) = y .= p .* u BatchIntegralFunction(boop, Float64[]) BatchIntegralFunction(boop, Float64[], max_batch = 20) BatchIntegralFunction(biip, Float64[], Float64[]) # the 2nd argument should be an ElasticArray -@test_throws SciMLBase.TooFewArgumentsError IntegralFunction(bi1) -@test_throws SciMLBase.TooManyArgumentsError IntegralFunction(bitoo) -@test_throws SciMLBase.TooManyArgumentsError IntegralFunction(bitoo, Float64[]) +@test_throws SciMLBase.TooFewArgumentsError BatchIntegralFunction(bi1, Float64[]) +@test_throws SciMLBase.TooManyArgumentsError BatchIntegralFunction(bitoo, Float64[], Float64[]) +@test_throws SciMLBase.TooManyArgumentsError BatchIntegralFunction(bitoo, Float64[]) From bbe691bb569c3f077cbc89ce3e846730d6b12045 Mon Sep 17 00:00:00 2001 From: lxvm Date: Wed, 20 Sep 2023 23:27:21 -0400 Subject: [PATCH 48/70] revert naming to integrand_prototype --- src/scimlfunctions.jl | 87 ++++++++++++++++++++++--------------------- 1 file changed, 44 insertions(+), 43 deletions(-) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index 30ab6d4aa..997069573 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -205,9 +205,9 @@ end const INTEGRAND_MISMATCH_FUNCTIONS_ERROR_MESSAGE = """ Nonconforming functions detected. If an integrand function `f` is defined - as out-of-place (`f(u,p)`), then no integral_prototype can be passed into the + as out-of-place (`f(u,p)`), then no integrand_prototype can be passed into the function constructor. Likewise if `f` is defined as in-place (`f(out,u,p)`), then - an integral_prototype is required. Either change the use of the function + an integrand_prototype is required. Either change the use of the function constructor or define the appropriate dispatch for `f`. """ @@ -2291,26 +2291,26 @@ A representation of an integrand `f` defined by: f(u, p) ``` -For an in-place form of `f` see the `iip` section below for details on in-place or out-of-place -handling. +For an in-place form of `f` see the `iip` section below for details on in-place or +out-of-place handling. ```julia -IntegralFunction{iip,specialize}(f, [integral_prototype]) +IntegralFunction{iip,specialize}(f, [integrand_prototype]) ``` Note that only `f` is required, and in the case of inplace integrands a mutable container -`integral_prototype` to store the result of the integral. If `integral_prototype` is present, -`f` is interpreted as in-place, and otherwise `f` is assumed to be out-of-place. +`integrand_prototype` to store the result of the integrand. If `integrand_prototype` is +present, `f` is interpreted as in-place, and otherwise `f` is assumed to be out-of-place. ## iip: In-Place vs Out-Of-Place Out-of-place functions must be of the form ``f(u, p)`` and in-place functions of the form ``f(y, u, p)``. Since `f` is allowed to return any type (e.g. real or complex numbers or -arrays), in-place functions must provide a container `integral_prototype` that is of the -right type for the final result of the integral, and the result is written to this container -in-place. When in-place forms are used, in-place array operations may be used by algorithms -to reduce allocations. If `integral_prototype` is not provided, `f` is assumed to be -out-of-place and quadrature is performed assuming immutable return types. +arrays), in-place functions must provide a container `integrand_prototype` that is of the +right type for the variable ``y``, and the result is written to this container in-place. +When in-place forms are used, in-place array operations, i.e. broadcasting, may be used by +algorithms to reduce allocations. If `integrand_prototype` is not provided, `f` is assumed +to be out-of-place and quadrature is performed assuming immutable return types. ## specialize @@ -2323,7 +2323,7 @@ The fields of the IntegralFunction type directly match the names of the inputs. struct IntegralFunction{iip, specialize, F, T} <: AbstractIntegralFunction{iip} f::F - integral_prototype::T + integrand_prototype::T end TruncatedStacktraces.@truncate_stacktrace IntegralFunction 1 2 @@ -2344,34 +2344,35 @@ output must be returned in the corresponding entries of ``y``. In general, the i algorithm is allowed to vary the number of evaluation points between subsequent calls to `f` ```julia -BatchIntegralFunction{iip,specialize}(f, output_prototype, [integral_prototype]; +BatchIntegralFunction{iip,specialize}(f, output_prototype, [integrand_prototype]; max_batch=typemax(Int)) ``` Note that `f` is required and a `resize`-able buffer `output_prototype` to store the output, -or range of `f`, and in the case of inplace integrands a mutable container -`integral_prototype` to store the result of the integral. These buffers can be reused across -multiple compatible integrals to reduce allocations. +or range of `f`, consisting of multiple integrand evaluations, and in the case of inplace +integrands a mutable container `integrand_prototype` to store the result of one integrand +evaluation. These buffers can be reused across multiple compatible integrals to reduce +allocations. The keyword `max_batch` is used to set a soft limit on the number of points to batch at the same time so that memory usage is controlled. -If `integral_prototype` is present, `f` is interpreted as in-place, and otherwise `f` is +If `integrand_prototype` is present, `f` is interpreted as in-place, and otherwise `f` is assumed to be out-of-place. ## iip: In-Place vs Out-Of-Place Out-of-place and in-place functions are both of the form ``f(y, u, p)``, but differ in the -element type of ``y``. Since `f` is allowed to return any type (e.g. real or complex numbers -or arrays), in-place functions must provide a container `integral_prototype` that is of the -right type for the final result of the integral, and the result is written to this container -in-place. When `f` is in-place, the buffer `output_prototype` is assumed to have a mutable -element type, and the last dimension of `output_prototype` should correspond to the batch -index. For example, `output_prototype` would have to be an `ElasticArray` or a -`VectorOfSimilarArrays` of an `ElasticArray`. When in-place forms are used, in-place array -operations may be used by algorithms to reduce allocations. If `integral_prototype` is not -provided, `f` is assumed to be out-of-place and quadrature is performed assuming -`output_prototype` is an `AbstractVector` with an immutable element type. +type of ``y``. Since `f` is allowed to return any type (e.g. real or complex numbers or +arrays), in-place functions must provide a container `output_prototype` of the right type +for ``y`` that stores multiple integrand evaluations. Typically, this means +`output_prototype` is a vector in the out-of-place case, or a matrix/array whose last +dimension is a batch index in the in-place case. For the in-place case, an +`integrand_prototype` is required and must be a container of the right type for a single +integrand evaluation. When `f` is in-place, the buffer `output_prototype` should be of the +type used by the integrand. When in-place forms are used, in-place array operations may be +used by algorithms to reduce allocations. If `integrand_prototype` is not provided, `f` is +assumed to be out-of-place. ## specialize @@ -2385,7 +2386,7 @@ struct BatchIntegralFunction{iip, specialize, F, Y, T} <: AbstractIntegralFunction{iip} f::F output_prototype::Y - integral_prototype::T + integrand_prototype::T max_batch::Int end @@ -4085,13 +4086,13 @@ function BVPFunction(f, bc; kwargs...) end BVPFunction(f::BVPFunction; kwargs...) = f -function IntegralFunction{iip, specialize}(f, integral_prototype) where {iip, specialize} - IntegralFunction{iip, specialize, typeof(f), typeof(integral_prototype)}(f, - integral_prototype) +function IntegralFunction{iip, specialize}(f, integrand_prototype) where {iip, specialize} + IntegralFunction{iip, specialize, typeof(f), typeof(integrand_prototype)}(f, + integrand_prototype) end -function IntegralFunction{iip}(f, integral_prototype) where {iip} - return IntegralFunction{iip, FullSpecialize}(f, integral_prototype) +function IntegralFunction{iip}(f, integrand_prototype) where {iip} + return IntegralFunction{iip, FullSpecialize}(f, integrand_prototype) end function IntegralFunction(f) calculated_iip = isinplace(f, 3, "integral", true) @@ -4100,44 +4101,44 @@ function IntegralFunction(f) end IntegralFunction{false}(f, nothing) end -function IntegralFunction(f, integral_prototype) +function IntegralFunction(f, integrand_prototype) calcuated_iip = isinplace(f, 3, "integral", true) if !calcuated_iip throw(IntegrandMismatchFunctionError(calcuated_iip, true)) end - IntegralFunction{true}(f, integral_prototype) + IntegralFunction{true}(f, integrand_prototype) end -function BatchIntegralFunction{iip, specialize}(f, output_prototype, integral_prototype; +function BatchIntegralFunction{iip, specialize}(f, output_prototype, integrand_prototype; max_batch::Integer = typemax(Int)) where {iip, specialize} BatchIntegralFunction{ iip, specialize, typeof(f), typeof(output_prototype), - typeof(integral_prototype), + typeof(integrand_prototype), }(f, output_prototype, - integral_prototype, + integrand_prototype, max_batch) end function BatchIntegralFunction{iip}(f, output_prototype, - integral_prototype; + integrand_prototype; kwargs...) where {iip} return BatchIntegralFunction{iip, FullSpecialize}(f, output_prototype, - integral_prototype; + integrand_prototype; kwargs...) end function BatchIntegralFunction(f, output_prototype; kwargs...) calcuated_iip = isinplace(f, 3, "batchintegral", true; has_two_dispatches = false) BatchIntegralFunction{false}(f, output_prototype, nothing; kwargs...) end -function BatchIntegralFunction(f, output_prototype, integral_prototype; kwargs...) +function BatchIntegralFunction(f, output_prototype, integrand_prototype; kwargs...) calcuated_iip = isinplace(f, 3, "batchintegral", true; has_two_dispatches = false) - BatchIntegralFunction{true}(f, output_prototype, integral_prototype; kwargs...) + BatchIntegralFunction{true}(f, output_prototype, integrand_prototype; kwargs...) end ########## Existence Functions From c0f40629f101b9383bbc3c657e1414f5c755a247 Mon Sep 17 00:00:00 2001 From: lxvm Date: Thu, 21 Sep 2023 00:45:54 -0400 Subject: [PATCH 49/70] wrap integrand with IntegralFunction in IntegralProblem --- src/problems/basic_problems.jl | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index e3efddcd9..41a27a137 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -335,11 +335,12 @@ which are `Number`s or `AbstractVector`s with the same geometry as `u`. ### Constructors ``` -IntegralProblem{iip}(f,lb,ub,p=NullParameters(); +IntegralProblem(f,lb,ub,p=NullParameters(); nout=1, batch = 0, kwargs...) ``` -- f: the integrand, callable function `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 (default) or an + `IntegralFunction` or `BatchIntegralFunction` for inplace and batching alternatives. - 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. @@ -372,7 +373,7 @@ struct IntegralProblem{isinplace, P, F, B, K} <: AbstractIntegralProblem{isinpla p::P batch::Int kwargs::K - @add_kwonly function IntegralProblem{iip}(f, lb, ub, p = NullParameters(); + @add_kwonly function IntegralProblem{iip}(f::AbstractIntegralFunction{iip}, lb, ub, p = NullParameters(); nout = 1, batch = 0, kwargs...) where {iip} @assert typeof(lb)==typeof(ub) "Type of lower and upper bound must match" @@ -385,9 +386,12 @@ end TruncatedStacktraces.@truncate_stacktrace IntegralProblem 1 4 -function IntegralProblem(f, lb, ub, args...; - kwargs...) - IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...) +function IntegralProblem(f::AbstractIntegralFunction, lb, ub, p = NullParameters(); kwargs...) + IntegralProblem{isinplace(f)}(f, lb, ub, p; kwargs...) +end + +function IntegralProblem(f, lb, ub, p = NullParameters(); kwargs...) + IntegralProblem(IntegralFunction(f), lb, ub, p; kwargs...) end struct QuadratureProblem end @@ -405,8 +409,8 @@ 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`. +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 @@ -415,10 +419,10 @@ assigned by a quadrature rule, which depend on sampling points `x`. ``` 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` +- 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`. +- 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. From 740576a76b0acf320632d64e12335c952b833d3e Mon Sep 17 00:00:00 2001 From: lxvm Date: Thu, 21 Sep 2023 01:05:56 -0400 Subject: [PATCH 50/70] make integral functions callable --- src/scimlfunctions.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index 997069573..c9aef6d0a 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -2397,6 +2397,8 @@ TruncatedStacktraces.@truncate_stacktrace BatchIntegralFunction 1 2 (f::ODEFunction)(args...) = f.f(args...) (f::NonlinearFunction)(args...) = f.f(args...) (f::IntervalNonlinearFunction)(args...) = f.f(args...) +(f::IntegralFunction)(args...) = f.f(args...) +(f::BatchIntegralFunction)(args...) = f.f(args...) function (f::DynamicalODEFunction)(u, p, t) ArrayPartition(f.f1(u.x[1], u.x[2], p, t), f.f2(u.x[1], u.x[2], p, t)) From 3f7d1fb825957f4e514ad83362bb19bbdcd26b94 Mon Sep 17 00:00:00 2001 From: lxvm Date: Thu, 21 Sep 2023 01:18:50 -0400 Subject: [PATCH 51/70] simplify IntegralProblem definition --- src/problems/basic_problems.jl | 46 +++++++++++++--------------------- 1 file changed, 17 insertions(+), 29 deletions(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 41a27a137..0401dc228 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -335,27 +335,18 @@ which are `Number`s or `AbstractVector`s with the same geometry as `u`. ### Constructors ``` +IntegralProblem(f,domain,p=NullParameters(); + nout=1, batch = 0, kwargs...) IntegralProblem(f,lb,ub,p=NullParameters(); nout=1, batch = 0, kwargs...) ``` - f: the integrand, callable function `y = f(u,p)` for out-of-place (default) or an - `IntegralFunction` or `BatchIntegralFunction` for inplace and batching alternatives. + `IntegralFunction` or `BatchIntegralFunction` for inplace and batching optimizations. +- domain: an object representing an integration domain, i.e. the tuple `(lb, ub)`. - 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. -- nout: The output size of the function f. Defaults to 1, i.e., a scalar valued function. - If `nout > 1` f is a vector valued function . -- batch: The preferred number of points to batch. This allows user-side parallelization - of the integrand. If `batch == 0` no batching is performed. - If `batch > 0` both `u` and `y` get an additional dimension added to it. - This means that: - if `f` is a multi variable function each `u[:,i]` is a different point to evaluate `f` at, - if `f` is a single variable function each `u[i]` is a different point to evaluate `f` at, - if `f` is a vector valued function each `y[:,i]` is the evaluation of `f` at a different point, - if `f` is a scalar valued function `y[i]` is the evaluation of `f` at a different point. - Note that batch is a suggestion for the number of points, - and it is not necessarily true that batch is the same as batchsize in all algorithms. - kwargs: Keyword arguments copied to the solvers. Additionally, we can supply iip like IntegralProblem{iip}(...) as true or false to declare at @@ -365,33 +356,30 @@ compile time whether the integrator function is in-place. The fields match the names of the constructor arguments. """ -struct IntegralProblem{isinplace, P, F, B, K} <: AbstractIntegralProblem{isinplace} +struct IntegralProblem{isinplace, P, F, T, K} <: AbstractIntegralProblem{isinplace} f::F - lb::B - ub::B - nout::Int + domain::T p::P - batch::Int kwargs::K - @add_kwonly function IntegralProblem{iip}(f::AbstractIntegralFunction{iip}, lb, ub, p = NullParameters(); - nout = 1, - batch = 0, kwargs...) where {iip} - @assert typeof(lb)==typeof(ub) "Type of lower and upper bound must match" + @add_kwonly function IntegralProblem{iip}(f::AbstractIntegralFunction{iip}, domain, p = NullParameters(); + kwargs...) where {iip} warn_paramtype(p) - new{iip, typeof(p), typeof(f), typeof(lb), typeof(kwargs)}(f, - lb, ub, nout, p, - batch, kwargs) + new{iip, typeof(p), typeof(f), typeof(domain), typeof(kwargs)}(f, + domain, p, kwargs) end end TruncatedStacktraces.@truncate_stacktrace IntegralProblem 1 4 -function IntegralProblem(f::AbstractIntegralFunction, lb, ub, p = NullParameters(); kwargs...) - IntegralProblem{isinplace(f)}(f, lb, ub, p; kwargs...) +function IntegralProblem(f::AbstractIntegralFunction, domain, p = NullParameters(); kwargs...) + IntegralProblem{isinplace(f)}(f, domain, p; kwargs...) end -function IntegralProblem(f, lb, ub, p = NullParameters(); kwargs...) - IntegralProblem(IntegralFunction(f), lb, ub, p; kwargs...) +function IntegralProblem(f::AbstractIntegralFunction, lb::B, ub::B, p = NullParameters(); kwargs...) where {B} + IntegralProblem(f, (lb, ub), p; kwargs...) +end +function IntegralProblem(f, args...; kwargs...) + IntegralProblem(IntegralFunction(f), args...; kwargs...) end struct QuadratureProblem end From 0deeefb3281bb16a52c50c4884d1f384f9a37cc8 Mon Sep 17 00:00:00 2001 From: lxvm Date: Thu, 21 Sep 2023 01:21:22 -0400 Subject: [PATCH 52/70] update docstrings --- src/problems/basic_problems.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 0401dc228..6e6b6269f 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -335,10 +335,8 @@ which are `Number`s or `AbstractVector`s with the same geometry as `u`. ### Constructors ``` -IntegralProblem(f,domain,p=NullParameters(); - nout=1, batch = 0, kwargs...) -IntegralProblem(f,lb,ub,p=NullParameters(); - nout=1, batch = 0, kwargs...) +IntegralProblem(f,domain,p=NullParameters(); kwargs...) +IntegralProblem(f,lb,ub,p=NullParameters(); kwargs...) ``` - f: the integrand, callable function `y = f(u,p)` for out-of-place (default) or an From e27965d52a8dfd11ff874379fa8ad1ca2ca1e757 Mon Sep 17 00:00:00 2001 From: lxvm Date: Thu, 21 Sep 2023 01:21:53 -0400 Subject: [PATCH 53/70] apply format --- src/problems/basic_problems.jl | 14 +++++++++++--- test/function_building_error_messages.jl | 4 +++- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 6e6b6269f..50621ccb9 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -359,7 +359,8 @@ struct IntegralProblem{isinplace, P, F, T, K} <: AbstractIntegralProblem{isinpla domain::T p::P kwargs::K - @add_kwonly function IntegralProblem{iip}(f::AbstractIntegralFunction{iip}, domain, p = NullParameters(); + @add_kwonly function IntegralProblem{iip}(f::AbstractIntegralFunction{iip}, domain, + p = NullParameters(); kwargs...) where {iip} warn_paramtype(p) new{iip, typeof(p), typeof(f), typeof(domain), typeof(kwargs)}(f, @@ -369,11 +370,18 @@ end TruncatedStacktraces.@truncate_stacktrace IntegralProblem 1 4 -function IntegralProblem(f::AbstractIntegralFunction, domain, p = NullParameters(); kwargs...) +function IntegralProblem(f::AbstractIntegralFunction, + domain, + p = NullParameters(); + kwargs...) IntegralProblem{isinplace(f)}(f, domain, p; kwargs...) end -function IntegralProblem(f::AbstractIntegralFunction, lb::B, ub::B, p = NullParameters(); kwargs...) where {B} +function IntegralProblem(f::AbstractIntegralFunction, + lb::B, + ub::B, + p = NullParameters(); + kwargs...) where {B} IntegralProblem(f, (lb, ub), p; kwargs...) end function IntegralProblem(f, args...; kwargs...) diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index cd47a6a33..008758ddb 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -629,5 +629,7 @@ BatchIntegralFunction(boop, Float64[]) BatchIntegralFunction(boop, Float64[], max_batch = 20) BatchIntegralFunction(biip, Float64[], Float64[]) # the 2nd argument should be an ElasticArray @test_throws SciMLBase.TooFewArgumentsError BatchIntegralFunction(bi1, Float64[]) -@test_throws SciMLBase.TooManyArgumentsError BatchIntegralFunction(bitoo, Float64[], Float64[]) +@test_throws SciMLBase.TooManyArgumentsError BatchIntegralFunction(bitoo, + Float64[], + Float64[]) @test_throws SciMLBase.TooManyArgumentsError BatchIntegralFunction(bitoo, Float64[]) From ee4931ceb04730185b0a8c84d6c94bb2f1f4a7dc Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 21 Sep 2023 13:28:45 -0400 Subject: [PATCH 54/70] Fix the remake issue with problem_type --- src/problems/bvp_problems.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/problems/bvp_problems.jl b/src/problems/bvp_problems.jl index b6fe87866..2e9f1feab 100644 --- a/src/problems/bvp_problems.jl +++ b/src/problems/bvp_problems.jl @@ -109,10 +109,16 @@ struct BVProblem{uType, tType, isinplace, P, F, BF, PT, K} <: kwargs::K @add_kwonly function BVProblem{iip}(f::AbstractBVPFunction{iip, TP}, bc, u0, tspan, - p = NullParameters(); kwargs...) where {iip, TP} + p = NullParameters(); problem_type=nothing, kwargs...) where {iip, TP} _tspan = promote_tspan(tspan) warn_paramtype(p) - problem_type = TP ? TwoPointBVProblem() : StandardBVProblem() + prob_type = TP ? TwoPointBVProblem() : StandardBVProblem() + # Needed to ensure that `problem_type` doesn't get passed in kwargs + if problem_type === nothing + problem_type = prob_type + else + @assert prob_type === problem_type "This indicates incorrect problem type specification! Users should never pass in `problem_type` kwarg, this exists exclusively for internal use." + end return new{typeof(u0), typeof(_tspan), iip, typeof(p), typeof(f), typeof(bc), typeof(problem_type), typeof(kwargs)}(f, bc, u0, _tspan, p, problem_type, kwargs) From 7b2d7f74aa5d2b77d91b00aa997f4c497b321299 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 21 Sep 2023 13:44:02 -0400 Subject: [PATCH 55/70] Fix BVP tests --- Project.toml | 2 +- test/downstream/ensemble_bvp.jl | 2 +- test/function_building_error_messages.jl | 49 +++++++++++++++--------- 3 files changed, 32 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index c63cead23..492af0df1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SciMLBase" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" authors = ["Chris Rackauckas and contributors"] -version = "1.98.2" # Remember to make it 2.0.0 +version = "2.0.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" diff --git a/test/downstream/ensemble_bvp.jl b/test/downstream/ensemble_bvp.jl index 4f3b2bc2a..ad08236c7 100644 --- a/test/downstream/ensemble_bvp.jl +++ b/test/downstream/ensemble_bvp.jl @@ -19,4 +19,4 @@ tspan = (0.0, pi / 2) p = [rand()] bvp = BVProblem(ode!, bc!, initial_guess, tspan, p) ensemble_prob = EnsembleProblem(bvp, prob_func = prob_func) -sim = solve(ensemble_prob, GeneralMIRK4(), trajectories = 10, dt = 0.1) +sim = solve(ensemble_prob, MIRK4(), trajectories = 10, dt = 0.1) diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index c424a6f84..279c80e73 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -463,7 +463,7 @@ IntegralProblem(intf, [0.0], [1.0], p) x = [1.0, 2.0] y = rand(2, 2) SampledIntegralProblem(y, x) -SampledIntegralProblem(y, x; dim=2) +SampledIntegralProblem(y, x; dim = 2) # Optimization @@ -516,10 +516,7 @@ bcjac(u, p, t) = [1.0] bcoop, jac = bjac, bcjac = bcjac) -@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip, - bciip, - jac = bjac, - bcjac = bcjac) +BVPFunction(bfiip, bciip, jac = bjac, bcjac = bcjac) @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, jac = bjac, @@ -528,8 +525,14 @@ BVPFunction(bfoop, bcoop, jac = bjac) bjac(du, u, p, t) = [1.0] bcjac(du, u, p, t) = [1.0] BVPFunction(bfiip, bciip, jac = bjac, bcjac = bcjac) -BVPFunction(bfoop, bciip, jac = bjac, bcjac = bcjac) -BVPFunction(bfiip, bcoop, jac = bjac, bcjac = bcjac) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, + bciip, + jac = bjac, + bcjac = bcjac) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip, + bcoop, + jac = bjac, + bcjac = bcjac) BVPFunction(bfoop, bcoop, jac = bjac, bcjac = bcjac) bWfact(u, t) = [1.0] @@ -540,10 +543,10 @@ bWfact(u, p, t) = [1.0] @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfoop, bciip, Wfact = bWfact) bWfact(u, p, gamma, t) = [1.0] @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip, bciip, Wfact = bWfact) -BVPFunction(bfoop, bciip, Wfact = bWfact) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, Wfact = bWfact) bWfact(du, u, p, gamma, t) = [1.0] BVPFunction(bfiip, bciip, Wfact = bWfact) -BVPFunction(bfoop, bciip, Wfact = bWfact) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, Wfact = bWfact) bWfact_t(u, t) = [1.0] @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfiip, bciip, Wfact_t = bWfact_t) @@ -555,20 +558,24 @@ bWfact_t(u, p, gamma, t) = [1.0] @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip, bciip, Wfact_t = bWfact_t) -BVPFunction(bfoop, bciip, Wfact_t = bWfact_t) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, + bciip, + Wfact_t = bWfact_t) bWfact_t(du, u, p, gamma, t) = [1.0] BVPFunction(bfiip, bciip, Wfact_t = bWfact_t) -BVPFunction(bfoop, bciip, Wfact_t = bWfact_t) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, + bciip, + Wfact_t = bWfact_t) btgrad(u, t) = [1.0] @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfiip, bciip, tgrad = btgrad) @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfoop, bciip, tgrad = btgrad) btgrad(u, p, t) = [1.0] @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip, bciip, tgrad = btgrad) -BVPFunction(bfoop, bciip, tgrad = btgrad) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, tgrad = btgrad) btgrad(du, u, p, t) = [1.0] BVPFunction(bfiip, bciip, tgrad = btgrad) -BVPFunction(bfoop, bciip, tgrad = btgrad) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, tgrad = btgrad) bparamjac(u, t) = [1.0] @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfiip, bciip, paramjac = bparamjac) @@ -577,27 +584,31 @@ bparamjac(u, p, t) = [1.0] @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip, bciip, paramjac = bparamjac) -BVPFunction(bfoop, bciip, paramjac = bparamjac) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, + bciip, + paramjac = bparamjac) bparamjac(du, u, p, t) = [1.0] BVPFunction(bfiip, bciip, paramjac = bparamjac) -BVPFunction(bfoop, bciip, paramjac = bparamjac) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, + bciip, + paramjac = bparamjac) bjvp(u, p, t) = [1.0] @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfiip, bciip, jvp = bjvp) @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfoop, bciip, jvp = bjvp) bjvp(u, v, p, t) = [1.0] @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip, bciip, jvp = bjvp) -BVPFunction(bfoop, bciip, jvp = bjvp) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, jvp = bjvp) bjvp(du, u, v, p, t) = [1.0] BVPFunction(bfiip, bciip, jvp = bjvp) -BVPFunction(bfoop, bciip, jvp = bjvp) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, jvp = bjvp) bvjp(u, p, t) = [1.0] @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfiip, bciip, vjp = bvjp) @test_throws SciMLBase.TooFewArgumentsError BVPFunction(bfoop, bciip, vjp = bvjp) bvjp(u, v, p, t) = [1.0] @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip, bciip, vjp = bvjp) -BVPFunction(bfoop, bciip, vjp = bvjp) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, vjp = bvjp) bvjp(du, u, v, p, t) = [1.0] BVPFunction(bfiip, bciip, vjp = bvjp) -BVPFunction(bfoop, bciip, vjp = bvjp) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, vjp = bvjp) From 8b173b96e74bab2b775981775b20d8c0b94bb138 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 21 Sep 2023 13:58:08 -0400 Subject: [PATCH 56/70] Fix 1 last test --- test/function_building_error_messages.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index 279c80e73..1db0fb468 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -516,7 +516,10 @@ bcjac(u, p, t) = [1.0] bcoop, jac = bjac, bcjac = bcjac) -BVPFunction(bfiip, bciip, jac = bjac, bcjac = bcjac) +@test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfiip, + bciip, + jac = bjac, + bcjac = bcjac) @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, jac = bjac, From a64db934fbefe1e581d981f22dca9700c139b201 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 21 Sep 2023 13:58:46 -0400 Subject: [PATCH 57/70] DONT MERGE TEST COMMIT --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 492af0df1..b15c4569b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SciMLBase" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" authors = ["Chris Rackauckas and contributors"] -version = "2.0.0" +version = "1.99.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 593e7f05ee51aea1800e86aed0b25f54de5eb17f Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 21 Sep 2023 14:33:20 -0400 Subject: [PATCH 58/70] Bump Project version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index b15c4569b..492af0df1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SciMLBase" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" authors = ["Chris Rackauckas and contributors"] -version = "1.99.0" +version = "2.0.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 5be7d7afbca32c2b84a688389c311d9efc83fed1 Mon Sep 17 00:00:00 2001 From: lxvm Date: Thu, 21 Sep 2023 14:35:12 -0400 Subject: [PATCH 59/70] remove output_prototype --- src/scimlfunctions.jl | 78 ++++++++++++------------ test/function_building_error_messages.jl | 22 ++++--- 2 files changed, 52 insertions(+), 48 deletions(-) diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index c9aef6d0a..71da8be5a 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -2283,7 +2283,7 @@ end TruncatedStacktraces.@truncate_stacktrace BVPFunction 1 2 @doc doc""" - IntegralFunction{iip,specialize,F} <: AbstractIntegralFunction{iip} + IntegralFunction{iip,specialize,F,T} <: AbstractIntegralFunction{iip} A representation of an integrand `f` defined by: @@ -2304,7 +2304,7 @@ present, `f` is interpreted as in-place, and otherwise `f` is assumed to be out- ## iip: In-Place vs Out-Of-Place -Out-of-place functions must be of the form ``f(u, p)`` and in-place functions of the form +Out-of-place functions must be of the form ``y = f(u, p)`` and in-place functions of the form ``f(y, u, p)``. Since `f` is allowed to return any type (e.g. real or complex numbers or arrays), in-place functions must provide a container `integrand_prototype` that is of the right type for the variable ``y``, and the result is written to this container in-place. @@ -2329,30 +2329,30 @@ end TruncatedStacktraces.@truncate_stacktrace IntegralFunction 1 2 @doc doc""" -BatchIntegralFunction{iip,specialize,F,T,Y,J,TJ,TPJ,Ta,S,JP,SP,TCV,O} <: -AbstractIntegralFunction{iip} +BatchIntegralFunction{iip,specialize,F,T} <: AbstractIntegralFunction{iip} A representation of an integrand `f` that can be evaluated at multiple points simultaneously using threads, the gpu, or distributed memory defined by: ```math -f(y, u, p) +y = f(u, p) ``` ``u`` is a vector whose elements correspond to distinct evaluation points to `f`, whose -output must be returned in the corresponding entries of ``y``. In general, the integration -algorithm is allowed to vary the number of evaluation points between subsequent calls to `f` +output must be returned as an array whose last "batching" dimension corresponds to integrand +evaluations at the different points in ``u``. In general, the integration algorithm is +allowed to vary the number of evaluation points between subsequent calls to `f`. + +For an in-place form of `f` see the `iip` section below for details on in-place or +out-of-place handling. ```julia -BatchIntegralFunction{iip,specialize}(f, output_prototype, [integrand_prototype]; +BatchIntegralFunction{iip,specialize}(f, [integrand_prototype]; max_batch=typemax(Int)) ``` - -Note that `f` is required and a `resize`-able buffer `output_prototype` to store the output, -or range of `f`, consisting of multiple integrand evaluations, and in the case of inplace -integrands a mutable container `integrand_prototype` to store the result of one integrand -evaluation. These buffers can be reused across multiple compatible integrals to reduce -allocations. +Note that only `f` is required, and in the case of inplace integrands a mutable container +`integrand_prototype` to store the result of the integrand of one integrand, without a last +"batching" dimension. The keyword `max_batch` is used to set a soft limit on the number of points to batch at the same time so that memory usage is controlled. @@ -2362,17 +2362,17 @@ assumed to be out-of-place. ## iip: In-Place vs Out-Of-Place -Out-of-place and in-place functions are both of the form ``f(y, u, p)``, but differ in the -type of ``y``. Since `f` is allowed to return any type (e.g. real or complex numbers or -arrays), in-place functions must provide a container `output_prototype` of the right type -for ``y`` that stores multiple integrand evaluations. Typically, this means -`output_prototype` is a vector in the out-of-place case, or a matrix/array whose last -dimension is a batch index in the in-place case. For the in-place case, an -`integrand_prototype` is required and must be a container of the right type for a single -integrand evaluation. When `f` is in-place, the buffer `output_prototype` should be of the -type used by the integrand. When in-place forms are used, in-place array operations may be -used by algorithms to reduce allocations. If `integrand_prototype` is not provided, `f` is -assumed to be out-of-place. +Out-of-place functions must be of the form ``y = f(u,p)`` and in-place functions of the form +``f(y, u, p)``. Since `f` is allowed to return any type (e.g. real or complex numbers or +arrays), in-place functions must provide a container `integrand_prototype` of the right type +for a single integrand evaluation. The integration algorithm will then allocate a ``y`` +array with the same element type as `integrand_prototype` and an additional last "batching" +dimension to store multiple integrand evaluations. In the out-of-place case, the algorithm +may infer the type of ``y`` by passing `f` an empty array of input points. This means ``y`` +is a vector in the out-of-place case, or a matrix/array in the in-place case. The number of +batched points may vary between subsequent calls to `f`. When in-place forms are used, +in-place array operations may be used by algorithms to reduce allocations. If +`integrand_prototype` is not provided, `f` is assumed to be out-of-place. ## specialize @@ -2382,10 +2382,9 @@ This field is currently unused The fields of the BatchIntegralFunction type directly match the names of the inputs. """ -struct BatchIntegralFunction{iip, specialize, F, Y, T} <: +struct BatchIntegralFunction{iip, specialize, F, T} <: AbstractIntegralFunction{iip} f::F - output_prototype::Y integrand_prototype::T max_batch::Int end @@ -4111,36 +4110,39 @@ function IntegralFunction(f, integrand_prototype) IntegralFunction{true}(f, integrand_prototype) end -function BatchIntegralFunction{iip, specialize}(f, output_prototype, integrand_prototype; +function BatchIntegralFunction{iip, specialize}(f, integrand_prototype; max_batch::Integer = typemax(Int)) where {iip, specialize} BatchIntegralFunction{ iip, specialize, typeof(f), - typeof(output_prototype), typeof(integrand_prototype), }(f, - output_prototype, integrand_prototype, max_batch) end function BatchIntegralFunction{iip}(f, - output_prototype, integrand_prototype; kwargs...) where {iip} return BatchIntegralFunction{iip, FullSpecialize}(f, - output_prototype, integrand_prototype; kwargs...) end -function BatchIntegralFunction(f, output_prototype; kwargs...) - calcuated_iip = isinplace(f, 3, "batchintegral", true; has_two_dispatches = false) - BatchIntegralFunction{false}(f, output_prototype, nothing; kwargs...) + +function BatchIntegralFunction(f; kwargs...) + calculated_iip = isinplace(f, 3, "batchintegral", true) + if calculated_iip + throw(IntegrandMismatchFunctionError(calculated_iip, false)) + end + BatchIntegralFunction{false}(f, nothing; kwargs...) end -function BatchIntegralFunction(f, output_prototype, integrand_prototype; kwargs...) - calcuated_iip = isinplace(f, 3, "batchintegral", true; has_two_dispatches = false) - BatchIntegralFunction{true}(f, output_prototype, integrand_prototype; kwargs...) +function BatchIntegralFunction(f, integrand_prototype; kwargs...) + calculated_iip = isinplace(f, 3, "batchintegral", true) + if !calculated_iip + throw(IntegrandMismatchFunctionError(calculated_iip, true)) + end + BatchIntegralFunction{true}(f, integrand_prototype; kwargs...) end ########## Existence Functions diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index 008758ddb..52b2054cc 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -620,16 +620,18 @@ IntegralFunction(iiip, Float64[]) # BatchIntegralFunction -boop(y, u, p) = y .= p .* u -biip(y, u, p) = y .= p .* u # this example is not realistic -bi1(y, u) = y .= p .* u +boop(u, p) = p .* u +biip(y, u, p) = y .= p .* u +bi1(u) = u bitoo(y, u, p, a) = y .= p .* u -BatchIntegralFunction(boop, Float64[]) -BatchIntegralFunction(boop, Float64[], max_batch = 20) -BatchIntegralFunction(biip, Float64[], Float64[]) # the 2nd argument should be an ElasticArray -@test_throws SciMLBase.TooFewArgumentsError BatchIntegralFunction(bi1, Float64[]) -@test_throws SciMLBase.TooManyArgumentsError BatchIntegralFunction(bitoo, - Float64[], - Float64[]) +BatchIntegralFunction(boop) +BatchIntegralFunction(boop, max_batch = 20) +BatchIntegralFunction(biip, Float64[]) +BatchIntegralFunction(biip, Float64[], max_batch = 20) + +@test_throws SciMLBase.IntegrandMismatchFunctionError BatchIntegralFunction(boop, Float64[]) +@test_throws SciMLBase.IntegrandMismatchFunctionError BatchIntegralFunction(biip) +@test_throws SciMLBase.TooFewArgumentsError BatchIntegralFunction(bi1) +@test_throws SciMLBase.TooManyArgumentsError BatchIntegralFunction(bitoo) @test_throws SciMLBase.TooManyArgumentsError BatchIntegralFunction(bitoo, Float64[]) From 8ebfe42cb07cbaee338c692585ef42b2504389ee Mon Sep 17 00:00:00 2001 From: lxvm Date: Thu, 21 Sep 2023 14:35:36 -0400 Subject: [PATCH 60/70] add deprecation method --- src/problems/basic_problems.jl | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 50621ccb9..5fcd36d9b 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -384,8 +384,28 @@ function IntegralProblem(f::AbstractIntegralFunction, kwargs...) where {B} IntegralProblem(f, (lb, ub), p; kwargs...) end + +# deprecation methods, which assume integrands return Float64 values (same as C libraries) +function IntegralProblem{iip}(f, args...; nout = 1, batch = 0, kwargs...) where {iip} + @warn "`nout` and `batch` keywords are deprecated in favor of inplace `IntegralFunction`s or `BatchIntegralFunction`s" + g = if iip + output_prototype = Vector{Float64}(undef, nout) + if batch == 0 + IntegralFunction(f, output_prototype) + else + BatchIntegralFunction(f, output_prototype, max_batch=batch) + end + else + if batch == 0 + IntegralFunction(f) + else + BatchIntegralFunction(f, max_batch=batch) + end + end + IntegralProblem(g, args...; kwargs...) +end function IntegralProblem(f, args...; kwargs...) - IntegralProblem(IntegralFunction(f), args...; kwargs...) + IntegralProblem{isinplace(f, 3)}(f, args...; kwargs...) end struct QuadratureProblem end From e6a0547a80437118f8223d73aba3901db6411de7 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Sep 2023 15:27:48 -0400 Subject: [PATCH 61/70] Update src/problems/basic_problems.jl --- src/problems/basic_problems.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 5fcd36d9b..ab4b05077 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -385,10 +385,14 @@ function IntegralProblem(f::AbstractIntegralFunction, IntegralProblem(f, (lb, ub), p; kwargs...) end -# deprecation methods, which assume integrands return Float64 values (same as C libraries) -function IntegralProblem{iip}(f, args...; nout = 1, batch = 0, kwargs...) where {iip} - @warn "`nout` and `batch` keywords are deprecated in favor of inplace `IntegralFunction`s or `BatchIntegralFunction`s" +function IntegralProblem(f, args...; nout = nothing, batch = nothing, kwargs...) + if nout !== nothing || batch !== nothing + @warn "`nout` and `batch` keywords are deprecated in favor of inplace `IntegralFunction`s or `BatchIntegralFunction`s. See the updated Integrals.jl documentation for details." + end + + iip = isinplace(f, 3) g = if iip + nout = nout === nothing ? 1 : nout output_prototype = Vector{Float64}(undef, nout) if batch == 0 IntegralFunction(f, output_prototype) @@ -402,10 +406,7 @@ function IntegralProblem{iip}(f, args...; nout = 1, batch = 0, kwargs...) where BatchIntegralFunction(f, max_batch=batch) end end - IntegralProblem(g, args...; kwargs...) -end -function IntegralProblem(f, args...; kwargs...) - IntegralProblem{isinplace(f, 3)}(f, args...; kwargs...) + IntegralProblem{iip}(g, args...; kwargs...) end struct QuadratureProblem end From 619ac07298b3195f64c792b965bc26f57732dfce Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Sep 2023 15:43:15 -0400 Subject: [PATCH 62/70] Update test/function_building_error_messages.jl --- test/function_building_error_messages.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index cf1dddfe6..54a81a724 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -614,7 +614,6 @@ bvjp(u, v, p, t) = [1.0] @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, vjp = bvjp) bvjp(du, u, v, p, t) = [1.0] BVPFunction(bfiip, bciip, vjp = bvjp) -BVPFunction(bfoop, bciip, vjp = bvjp) @test_throws SciMLBase.NonconformingFunctionsError BVPFunction(bfoop, bciip, vjp = bvjp) From 83f933dd3fc34ac25d743a0f94356f5fbf220d44 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Sep 2023 15:43:36 -0400 Subject: [PATCH 63/70] Update function_building_error_messages.jl --- test/function_building_error_messages.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/function_building_error_messages.jl b/test/function_building_error_messages.jl index 54a81a724..269985a39 100644 --- a/test/function_building_error_messages.jl +++ b/test/function_building_error_messages.jl @@ -649,4 +649,4 @@ BatchIntegralFunction(biip, Float64[], max_batch = 20) @test_throws SciMLBase.IntegrandMismatchFunctionError BatchIntegralFunction(biip) @test_throws SciMLBase.TooFewArgumentsError BatchIntegralFunction(bi1) @test_throws SciMLBase.TooManyArgumentsError BatchIntegralFunction(bitoo) -@test_throws SciMLBase.TooManyArgumentsError BatchIntegralFunction(bitoo, Float64[]) \ No newline at end of file +@test_throws SciMLBase.TooManyArgumentsError BatchIntegralFunction(bitoo, Float64[]) From 7740dd4844a0ec9e551eb97345bbf9344c300bbd Mon Sep 17 00:00:00 2001 From: lxvm Date: Thu, 21 Sep 2023 18:51:28 -0400 Subject: [PATCH 64/70] fix default batch and dispatch --- src/problems/basic_problems.jl | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/src/problems/basic_problems.jl b/src/problems/basic_problems.jl index 3a67cd4d3..aece5fa7b 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/basic_problems.jl @@ -389,24 +389,23 @@ function IntegralProblem(f, args...; nout = nothing, batch = nothing, kwargs...) if nout !== nothing || batch !== nothing @warn "`nout` and `batch` keywords are deprecated in favor of inplace `IntegralFunction`s or `BatchIntegralFunction`s. See the updated Integrals.jl documentation for details." end - - iip = isinplace(f, 3) - g = if iip - nout = nout === nothing ? 1 : nout - output_prototype = Vector{Float64}(undef, nout) - if batch == 0 + + max_batch = batch === nothing ? 0 : batch + g = if isinplace(f, 3) + output_prototype = Vector{Float64}(undef, nout === nothing ? 1 : nout) + if max_batch == 0 IntegralFunction(f, output_prototype) else - BatchIntegralFunction(f, output_prototype, max_batch=batch) + BatchIntegralFunction(f, output_prototype, max_batch=max_batch) end else - if batch == 0 + if max_batch == 0 IntegralFunction(f) else - BatchIntegralFunction(f, max_batch=batch) + BatchIntegralFunction(f, max_batch=max_batch) end end - IntegralProblem{iip}(g, args...; kwargs...) + IntegralProblem(g, args...; kwargs...) end struct QuadratureProblem end From a6fd63a5ec0f8dda7fb7fc73b66a28ceb7c788cd Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Sep 2023 19:13:41 -0400 Subject: [PATCH 65/70] Change version just to run tests --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 492af0df1..b15c4569b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SciMLBase" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" authors = ["Chris Rackauckas and contributors"] -version = "2.0.0" +version = "1.99.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From bb384b3826177a598368124cbf0d1300c83cf1b7 Mon Sep 17 00:00:00 2001 From: "William S. Moses" Date: Tue, 12 Sep 2023 21:11:27 -0400 Subject: [PATCH 66/70] Make ODEProblem into a mutable struct --- src/problems/ode_problems.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/problems/ode_problems.jl b/src/problems/ode_problems.jl index 20949d644..a4dcbc544 100644 --- a/src/problems/ode_problems.jl +++ b/src/problems/ode_problems.jl @@ -94,7 +94,7 @@ prob = ODEProblemLibrary.prob_ode_linear sol = solve(prob) ``` """ -struct ODEProblem{uType, tType, isinplace, P, F, K, PT} <: +mutable struct ODEProblem{uType, tType, isinplace, P, F, K, PT} <: AbstractODEProblem{uType, tType, isinplace} """The ODE is `du = f(u,p,t)` for out-of-place and f(du,u,p,t) for in-place.""" f::F From 97517ddb96ea6a48f13e3e15d5d3056a564d3799 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Thu, 21 Sep 2023 20:27:04 -0400 Subject: [PATCH 67/70] add a big warning on mutation of ODEProblem --- src/problems/ode_problems.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/problems/ode_problems.jl b/src/problems/ode_problems.jl index a4dcbc544..ce4425db3 100644 --- a/src/problems/ode_problems.jl +++ b/src/problems/ode_problems.jl @@ -162,6 +162,16 @@ mutable struct ODEProblem{uType, tType, isinplace, P, F, K, PT} <: end TruncatedStacktraces.@truncate_stacktrace ODEProblem 3 1 2 +function Base.setproperty!(prob::ODEProblem, s::Symbol, v) + @warn "Mutation of ODEProblem detected. SciMLBase v2.0 has made ODEProblem temporarily mutable in order to allow for interfacing with EnzymeRules due to a current limitation in the rule system. This change is only intended to be temporary and ODEProblem will return to being a struct in a later non-breaking release. Do not rely on this behavior, use with caution." + Base.setfield!(prob, s, v) +end + +function Base.setproperty!(prob::ODEProblem, s::Symbol, v, order::Symbol) + @warn "Mutation of ODEProblem detected. SciMLBase v2.0 has made ODEProblem temporarily mutable in order to allow for interfacing with EnzymeRules due to a current limitation in the rule system. This change is only intended to be temporary and ODEProblem will return to being a struct in a later non-breaking release. Do not rely on this behavior, use with caution." + Base.setfield!(prob, s, v, order) +end + """ ODEProblem(f::ODEFunction,u0,tspan,p=NullParameters(),callback=CallbackSet()) From bb9b42355e8cd3c30fb79bdfe070c0e771a02fe9 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Sep 2023 21:37:46 -0400 Subject: [PATCH 68/70] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index b15c4569b..492af0df1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SciMLBase" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" authors = ["Chris Rackauckas and contributors"] -version = "1.99.0" +version = "2.0.0" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" From 4ce6764cc6c5655f90543326e1cde420200c282c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 21 Sep 2023 21:41:29 -0400 Subject: [PATCH 69/70] Update README.md --- README.md | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/README.md b/README.md index 87749abf7..eed2c69f7 100644 --- a/README.md +++ b/README.md @@ -12,3 +12,13 @@ SciMLBase.jl is the core interface definition of the SciML ecosystem. It is a low dependency library made to be depended on by the downstream libraries to supply the common interface and allow for interexchange of mathematical problems. + +## v2.0 Breaking Changes + +The breaking changes in v2.0 are: + +* `IntegralProblem` has moved to an interface with `IntegralFunction` and `BatchedIntegralFunction` which requires specifying `prototype`s for the values to be modified + instead of `nout` and `batch`. +* `ODEProblem` was made temporarily into a `mutable struct` to allow for EnzymeRules support. Using the mutation throws a warning that this is only experimental and should not be relied on. +* `BVProblem` now has a new interface for `TwoPointBVProblem` which splits the bc terms for the two sides, forcing a true two-point BVProblem to allow for further specializations and to allow + for wrapping Fortran solvers in the interface. From 600fd4c778fed7dae928978dd6a9af2dc4c57ced Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 22 Sep 2023 09:22:45 -0400 Subject: [PATCH 70/70] Update README.md --- README.md | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index eed2c69f7..b294f98c8 100644 --- a/README.md +++ b/README.md @@ -18,7 +18,10 @@ supply the common interface and allow for interexchange of mathematical problems The breaking changes in v2.0 are: * `IntegralProblem` has moved to an interface with `IntegralFunction` and `BatchedIntegralFunction` which requires specifying `prototype`s for the values to be modified - instead of `nout` and `batch`. + instead of `nout` and `batch`. https://github.com/SciML/SciMLBase.jl/pull/497 * `ODEProblem` was made temporarily into a `mutable struct` to allow for EnzymeRules support. Using the mutation throws a warning that this is only experimental and should not be relied on. + https://github.com/SciML/SciMLBase.jl/pull/501 * `BVProblem` now has a new interface for `TwoPointBVProblem` which splits the bc terms for the two sides, forcing a true two-point BVProblem to allow for further specializations and to allow - for wrapping Fortran solvers in the interface. + for wrapping Fortran solvers in the interface. https://github.com/SciML/SciMLBase.jl/pull/477 +* `SDEProblem` constructor was changed to remove an anti-pattern which required passing the diffusion function `g` twice, i.e. `SDEProblem(SDEFunction(f,g),g, ...)`. + Now this is simply `SDEProblem(SDEFunction(f,g),...)`. https://github.com/SciML/SciMLBase.jl/pull/489