Skip to content

Commit

Permalink
Refactor basic problems into sets of problems
Browse files Browse the repository at this point in the history
  • Loading branch information
ChrisRackauckas committed Jul 21, 2024
1 parent 8b20758 commit 58592a8
Show file tree
Hide file tree
Showing 4 changed files with 250 additions and 250 deletions.
4 changes: 3 additions & 1 deletion src/SciMLBase.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
169 changes: 169 additions & 0 deletions src/problems/integral_problems.jl
Original file line number Diff line number Diff line change
@@ -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
78 changes: 78 additions & 0 deletions src/problems/linear_problems.jl
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit 58592a8

Please sign in to comment.