diff --git a/src/SciMLBase.jl b/src/SciMLBase.jl index 9a021fceb..fc1623d25 100644 --- a/src/SciMLBase.jl +++ b/src/SciMLBase.jl @@ -700,7 +700,9 @@ include("problems/discrete_problems.jl") include("problems/implicit_discrete_problems.jl") include("problems/steady_state_problems.jl") include("problems/analytical_problems.jl") -include("problems/basic_problems.jl") +include("problems/linear_problems.jl") +include("problems/nonlinear_problems.jl") +include("problems/integral_problems.jl") include("problems/ode_problems.jl") include("problems/rode_problems.jl") include("problems/sde_problems.jl") diff --git a/src/problems/integral_problems.jl b/src/problems/integral_problems.jl new file mode 100644 index 000000000..e68b764ea --- /dev/null +++ b/src/problems/integral_problems.jl @@ -0,0 +1,169 @@ +@doc doc""" + +Defines an integral problem. +Documentation Page: https://docs.sciml.ai/Integrals/stable/ + +## Mathematical Specification of an Integral Problem + +Integral problems are multi-dimensional integrals defined as: + +```math +\int_{lb}^{ub} f(u,p) du +``` + +where `p` are parameters. `u` is a `Number` or `AbstractVector` +whose geometry matches the space being integrated. +This space is bounded by the lowerbound `lb` and upperbound `ub`, +which are `Number`s or `AbstractVector`s with the same geometry as `u`. + +## Problem Type + +### Constructors + +``` +IntegralProblem(f::AbstractIntegralFunction,domain,p=NullParameters(); kwargs...) +IntegralProblem(f::AbstractIntegralFunction,lb,ub,p=NullParameters(); kwargs...) +IntegralProblem(f,domain,p=NullParameters(); nout=nothing, batch=nothing, kwargs...) +IntegralProblem(f,lb,ub,p=NullParameters(); nout=nothing, batch=nothing, kwargs...) +``` + +- f: the integrand, callable function `y = f(u,p)` for out-of-place (default) or an + `IntegralFunction` or `BatchIntegralFunction` for inplace and batching optimizations. +- domain: an object representing an integration domain, i.e. the tuple `(lb, ub)`. +- lb: DEPRECATED: Either a number or vector of lower bounds. +- ub: DEPRECATED: Either a number or vector of upper bounds. +- p: The parameters associated with the problem. +- nout: DEPRECATED (see `IntegralFunction`): length of the vector output of the integrand + (by default the integrand is assumed to be scalar) +- batch: DEPRECATED (see `BatchIntegralFunction`): number of points the integrand can + evaluate simultaneously (by default there is no batching) +- kwargs: Keyword arguments copied to the solvers. + +Additionally, we can supply iip like IntegralProblem{iip}(...) as true or false to declare at +compile time whether the integrator function is in-place. + +### Fields + +The fields match the names of the constructor arguments. +""" +struct IntegralProblem{isinplace, P, F, T, K} <: AbstractIntegralProblem{isinplace} + f::F + domain::T + p::P + kwargs::K + @add_kwonly function IntegralProblem{iip}(f::AbstractIntegralFunction{iip}, domain, + p = NullParameters(); nout = nothing, batch = nothing, + kwargs...) where {iip} + warn_paramtype(p) + new{iip, typeof(p), typeof(f), typeof(domain), typeof(kwargs)}(f, + domain, p, kwargs) + end +end + +function IntegralProblem(f::AbstractIntegralFunction, + domain, + p = NullParameters(); + kwargs...) + IntegralProblem{isinplace(f)}(f, domain, p; kwargs...) +end + +@deprecate IntegralProblem{iip}(f::AbstractIntegralFunction, + lb::Union{Number, AbstractVector{<:Number}}, + ub::Union{Number, AbstractVector{<:Number}}, + p = NullParameters(); kwargs...) where {iip} IntegralProblem{iip}( + f, (lb, ub), p; kwargs...) + +function IntegralProblem(f, args...; kwargs...) + IntegralProblem{isinplace(f, 3)}(f, args...; kwargs...) +end +function IntegralProblem{iip}( + f, args...; nout = nothing, batch = nothing, kwargs...) where {iip} + 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 + + g = if iip + if batch === nothing + output_prototype = nout === nothing ? Array{Float64, 0}(undef) : + Vector{Float64}(undef, nout) + IntegralFunction(f, output_prototype) + else + output_prototype = nout === nothing ? Float64[] : + Matrix{Float64}(undef, nout, 0) + BatchIntegralFunction(f, output_prototype, max_batch = batch) + end + else + if batch === nothing + IntegralFunction(f) + else + BatchIntegralFunction(f, max_batch = batch) + end + end + IntegralProblem(g, args...; kwargs...) +end + +function Base.getproperty(prob::IntegralProblem, name::Symbol) + if name === :lb + domain = getfield(prob, :domain) + lb, ub = domain + return lb + elseif name === :ub + domain = getfield(prob, :domain) + lb, ub = domain + return ub + elseif name === :ps + return ParameterIndexingProxy(prob) + end + return Base.getfield(prob, name) +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 + +Sampled integral problems are defined as: + +```math +\sum_i w_i y_i +``` +where `y_i` are sampled values of the integrand, and `w_i` are weights +assigned by a quadrature rule, which depend on sampling points `x`. + +## Problem Type + +### Constructors + +``` +SampledIntegralProblem(y::AbstractArray, x::AbstractVector; dim=ndims(y), kwargs...) +``` +- y: The sampled integrand, must be a subtype of `AbstractArray`. + It is assumed that the values of `y` along dimension `dim` + correspond to the integrand evaluated at sampling points `x` +- x: Sampling points, must be a subtype of `AbstractVector`. +- dim: Dimension along which to integrate. Defaults to the last dimension of `y`. +- kwargs: Keyword arguments copied to the solvers. + +### Fields + +The fields match the names of the constructor arguments. +""" +struct SampledIntegralProblem{Y, X, K} <: AbstractIntegralProblem{false} + y::Y + x::X + dim::Int + kwargs::K + @add_kwonly function SampledIntegralProblem(y::AbstractArray, x::AbstractVector; + dim = ndims(y), + kwargs...) + @assert dim<=ndims(y) "The integration dimension `dim` is larger than the number of dimensions of the integrand `y`" + @assert length(x)==size(y, dim) "The integrand `y` must have the same length as the sampling points `x` along the integrated dimension." + @assert axes(x, 1)==axes(y, dim) "The integrand `y` must obey the same indexing as the sampling points `x` along the integrated dimension." + new{typeof(y), typeof(x), typeof(kwargs)}(y, x, dim, kwargs) + end +end diff --git a/src/problems/linear_problems.jl b/src/problems/linear_problems.jl new file mode 100644 index 000000000..8a1cc5d35 --- /dev/null +++ b/src/problems/linear_problems.jl @@ -0,0 +1,78 @@ +@doc doc""" + +Defines a linear system problem. +Documentation Page: https://docs.sciml.ai/LinearSolve/stable/basics/LinearProblem/ + +## Mathematical Specification of a Linear Problem + +### Concrete LinearProblem + +To define a `LinearProblem`, you simply need to give the `AbstractMatrix` ``A`` +and an `AbstractVector` ``b`` which defines the linear system: + +```math +Au = b +``` + +### Matrix-Free LinearProblem + +For matrix-free versions, the specification of the problem is given by an +operator `A(u,p,t)` which computes `A*u`, or in-place as `A(du,u,p,t)`. These +are specified via the `AbstractSciMLOperator` interface. For more details, see +the [SciMLBase Documentation](https://docs.sciml.ai/SciMLBase/stable/). + +Note that matrix-free versions of LinearProblem definitions are not compatible +with all solvers. To check a solver for compatibility, use the function xxxxx. + +## Problem Type + +### Constructors + +Optionally, an initial guess ``u₀`` can be supplied which is used for iterative +methods. + +```julia +LinearProblem{isinplace}(A,x,p=NullParameters();u0=nothing,kwargs...) +LinearProblem(f::AbstractSciMLOperator,u0,p=NullParameters();u0=nothing,kwargs...) +``` + +`isinplace` optionally sets whether the function is in-place or not, i.e. whether +the solvers are allowed to mutate. By default this is true for `AbstractMatrix`, +and for `AbstractSciMLOperator`s it matches the choice of the operator definition. + +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. + +### Fields + +* `A`: The representation of the linear operator. +* `b`: The right-hand side of the linear system. +* `p`: The parameters for the problem. Defaults to `NullParameters`. Currently unused. +* `u0`: The initial condition used by iterative solvers. +* `kwargs`: The keyword arguments passed on to the solvers. +""" +struct LinearProblem{uType, isinplace, F, bType, P, K} <: + AbstractLinearProblem{bType, isinplace} + A::F + b::bType + u0::uType + p::P + kwargs::K + @add_kwonly function LinearProblem{iip}(A, b, p = NullParameters(); u0 = nothing, + kwargs...) where {iip} + warn_paramtype(p) + new{typeof(u0), iip, typeof(A), typeof(b), typeof(p), typeof(kwargs)}(A, b, u0, p, + kwargs) + end +end + +function LinearProblem(A, b, args...; kwargs...) + if A isa AbstractArray + LinearProblem{true}(A, b, args...; kwargs...) + elseif A isa Number + LinearProblem{false}(A, b, args...; kwargs...) + else + LinearProblem{isinplace(A, 4)}(A, b, args...; kwargs...) + end +end \ No newline at end of file diff --git a/src/problems/basic_problems.jl b/src/problems/nonlinear_problems.jl similarity index 56% rename from src/problems/basic_problems.jl rename to src/problems/nonlinear_problems.jl index dbf8ff5e3..78a6e05ef 100644 --- a/src/problems/basic_problems.jl +++ b/src/problems/nonlinear_problems.jl @@ -1,82 +1,3 @@ -@doc doc""" - -Defines a linear system problem. -Documentation Page: https://docs.sciml.ai/LinearSolve/stable/basics/LinearProblem/ - -## Mathematical Specification of a Linear Problem - -### Concrete LinearProblem - -To define a `LinearProblem`, you simply need to give the `AbstractMatrix` ``A`` -and an `AbstractVector` ``b`` which defines the linear system: - -```math -Au = b -``` - -### Matrix-Free LinearProblem - -For matrix-free versions, the specification of the problem is given by an -operator `A(u,p,t)` which computes `A*u`, or in-place as `A(du,u,p,t)`. These -are specified via the `AbstractSciMLOperator` interface. For more details, see -the [SciMLBase Documentation](https://docs.sciml.ai/SciMLBase/stable/). - -Note that matrix-free versions of LinearProblem definitions are not compatible -with all solvers. To check a solver for compatibility, use the function xxxxx. - -## Problem Type - -### Constructors - -Optionally, an initial guess ``u₀`` can be supplied which is used for iterative -methods. - -```julia -LinearProblem{isinplace}(A,x,p=NullParameters();u0=nothing,kwargs...) -LinearProblem(f::AbstractSciMLOperator,u0,p=NullParameters();u0=nothing,kwargs...) -``` - -`isinplace` optionally sets whether the function is in-place or not, i.e. whether -the solvers are allowed to mutate. By default this is true for `AbstractMatrix`, -and for `AbstractSciMLOperator`s it matches the choice of the operator definition. - -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. - -### Fields - -* `A`: The representation of the linear operator. -* `b`: The right-hand side of the linear system. -* `p`: The parameters for the problem. Defaults to `NullParameters`. Currently unused. -* `u0`: The initial condition used by iterative solvers. -* `kwargs`: The keyword arguments passed on to the solvers. -""" -struct LinearProblem{uType, isinplace, F, bType, P, K} <: - AbstractLinearProblem{bType, isinplace} - A::F - b::bType - u0::uType - p::P - kwargs::K - @add_kwonly function LinearProblem{iip}(A, b, p = NullParameters(); u0 = nothing, - kwargs...) where {iip} - warn_paramtype(p) - new{typeof(u0), iip, typeof(A), typeof(b), typeof(p), typeof(kwargs)}(A, b, u0, p, - kwargs) - end -end - -function LinearProblem(A, b, args...; kwargs...) - if A isa AbstractArray - LinearProblem{true}(A, b, args...; kwargs...) - elseif A isa Number - LinearProblem{false}(A, b, args...; kwargs...) - else - LinearProblem{isinplace(A, 4)}(A, b, args...; kwargs...) - end -end - """ $(TYPEDEF) """ @@ -400,173 +321,3 @@ end function NonlinearLeastSquaresProblem(f, u0, p = NullParameters(); kwargs...) return NonlinearLeastSquaresProblem(NonlinearFunction(f), u0, p; kwargs...) end - -@doc doc""" - -Defines an integral problem. -Documentation Page: https://docs.sciml.ai/Integrals/stable/ - -## Mathematical Specification of an Integral Problem - -Integral problems are multi-dimensional integrals defined as: - -```math -\int_{lb}^{ub} f(u,p) du -``` - -where `p` are parameters. `u` is a `Number` or `AbstractVector` -whose geometry matches the space being integrated. -This space is bounded by the lowerbound `lb` and upperbound `ub`, -which are `Number`s or `AbstractVector`s with the same geometry as `u`. - -## Problem Type - -### Constructors - -``` -IntegralProblem(f::AbstractIntegralFunction,domain,p=NullParameters(); kwargs...) -IntegralProblem(f::AbstractIntegralFunction,lb,ub,p=NullParameters(); kwargs...) -IntegralProblem(f,domain,p=NullParameters(); nout=nothing, batch=nothing, kwargs...) -IntegralProblem(f,lb,ub,p=NullParameters(); nout=nothing, batch=nothing, kwargs...) -``` - -- f: the integrand, callable function `y = f(u,p)` for out-of-place (default) or an - `IntegralFunction` or `BatchIntegralFunction` for inplace and batching optimizations. -- domain: an object representing an integration domain, i.e. the tuple `(lb, ub)`. -- lb: DEPRECATED: Either a number or vector of lower bounds. -- ub: DEPRECATED: Either a number or vector of upper bounds. -- p: The parameters associated with the problem. -- nout: DEPRECATED (see `IntegralFunction`): length of the vector output of the integrand - (by default the integrand is assumed to be scalar) -- batch: DEPRECATED (see `BatchIntegralFunction`): number of points the integrand can - evaluate simultaneously (by default there is no batching) -- kwargs: Keyword arguments copied to the solvers. - -Additionally, we can supply iip like IntegralProblem{iip}(...) as true or false to declare at -compile time whether the integrator function is in-place. - -### Fields - -The fields match the names of the constructor arguments. -""" -struct IntegralProblem{isinplace, P, F, T, K} <: AbstractIntegralProblem{isinplace} - f::F - domain::T - p::P - kwargs::K - @add_kwonly function IntegralProblem{iip}(f::AbstractIntegralFunction{iip}, domain, - p = NullParameters(); nout = nothing, batch = nothing, - kwargs...) where {iip} - warn_paramtype(p) - new{iip, typeof(p), typeof(f), typeof(domain), typeof(kwargs)}(f, - domain, p, kwargs) - end -end - -function IntegralProblem(f::AbstractIntegralFunction, - domain, - p = NullParameters(); - kwargs...) - IntegralProblem{isinplace(f)}(f, domain, p; kwargs...) -end - -@deprecate IntegralProblem{iip}(f::AbstractIntegralFunction, - lb::Union{Number, AbstractVector{<:Number}}, - ub::Union{Number, AbstractVector{<:Number}}, - p = NullParameters(); kwargs...) where {iip} IntegralProblem{iip}( - f, (lb, ub), p; kwargs...) - -function IntegralProblem(f, args...; kwargs...) - IntegralProblem{isinplace(f, 3)}(f, args...; kwargs...) -end -function IntegralProblem{iip}( - f, args...; nout = nothing, batch = nothing, kwargs...) where {iip} - 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 - - g = if iip - if batch === nothing - output_prototype = nout === nothing ? Array{Float64, 0}(undef) : - Vector{Float64}(undef, nout) - IntegralFunction(f, output_prototype) - else - output_prototype = nout === nothing ? Float64[] : - Matrix{Float64}(undef, nout, 0) - BatchIntegralFunction(f, output_prototype, max_batch = batch) - end - else - if batch === nothing - IntegralFunction(f) - else - BatchIntegralFunction(f, max_batch = batch) - end - end - IntegralProblem(g, args...; kwargs...) -end - -function Base.getproperty(prob::IntegralProblem, name::Symbol) - if name === :lb - domain = getfield(prob, :domain) - lb, ub = domain - return lb - elseif name === :ub - domain = getfield(prob, :domain) - lb, ub = domain - return ub - elseif name === :ps - return ParameterIndexingProxy(prob) - end - return Base.getfield(prob, name) -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 - -Sampled integral problems are defined as: - -```math -\sum_i w_i y_i -``` -where `y_i` are sampled values of the integrand, and `w_i` are weights -assigned by a quadrature rule, which depend on sampling points `x`. - -## Problem Type - -### Constructors - -``` -SampledIntegralProblem(y::AbstractArray, x::AbstractVector; dim=ndims(y), kwargs...) -``` -- y: The sampled integrand, must be a subtype of `AbstractArray`. - It is assumed that the values of `y` along dimension `dim` - correspond to the integrand evaluated at sampling points `x` -- x: Sampling points, must be a subtype of `AbstractVector`. -- dim: Dimension along which to integrate. Defaults to the last dimension of `y`. -- kwargs: Keyword arguments copied to the solvers. - -### Fields - -The fields match the names of the constructor arguments. -""" -struct SampledIntegralProblem{Y, X, K} <: AbstractIntegralProblem{false} - y::Y - x::X - dim::Int - kwargs::K - @add_kwonly function SampledIntegralProblem(y::AbstractArray, x::AbstractVector; - dim = ndims(y), - kwargs...) - @assert dim<=ndims(y) "The integration dimension `dim` is larger than the number of dimensions of the integrand `y`" - @assert length(x)==size(y, dim) "The integrand `y` must have the same length as the sampling points `x` along the integrated dimension." - @assert axes(x, 1)==axes(y, dim) "The integrand `y` must obey the same indexing as the sampling points `x` along the integrated dimension." - new{typeof(y), typeof(x), typeof(kwargs)}(y, x, dim, kwargs) - end -end