Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add DataIntegralProblem #491

Merged
merged 15 commits into from
Sep 17, 2023
Merged
Show file tree
Hide file tree
Changes from 14 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 1 addition & 3 deletions src/SciMLBase.jl
Original file line number Diff line number Diff line change
Expand Up @@ -760,9 +760,7 @@ export solve, solve!, init, discretize, symbolic_discretize

export LinearProblem,
NonlinearProblem, IntervalNonlinearProblem,
IntegralProblem, OptimizationProblem

export IntegralProblem
IntegralProblem, SampledIntegralProblem, OptimizationProblem

export DiscreteProblem, ImplicitDiscreteProblem
export SteadyStateProblem, SteadyStateSolution
Expand Down
60 changes: 57 additions & 3 deletions src/problems/basic_problems.jl
Original file line number Diff line number Diff line change
Expand Up @@ -334,10 +334,12 @@ which are `Number`s or `AbstractVector`s with the same geometry as `u`.

### Constructors

```
IntegralProblem{iip}(f,lb,ub,p=NullParameters();
nout=1, batch = 0, kwargs...)
```

- f: the integrand, `y = f(u,p)` for out-of-place or `f(y,u,p)` for in-place.
- f: the integrand, callable function `y = f(u,p)` for out-of-place or `f(y,u,p)` for in-place.
- lb: Either a number or vector of lower bounds.
- ub: Either a number or vector of upper bounds.
- p: The parameters associated with the problem.
Expand Down Expand Up @@ -375,14 +377,16 @@ struct IntegralProblem{isinplace, P, F, B, K} <: AbstractIntegralProblem{isinpla
batch = 0, kwargs...) where {iip}
@assert typeof(lb)==typeof(ub) "Type of lower and upper bound must match"
warn_paramtype(p)
new{iip, typeof(p), typeof(f), typeof(lb), typeof(kwargs)}(f, lb, ub, nout, p,
new{iip, typeof(p), typeof(f), typeof(lb), typeof(kwargs)}(f,
lb, ub, nout, p,
batch, kwargs)
end
end

TruncatedStacktraces.@truncate_stacktrace IntegralProblem 1 4

function IntegralProblem(f, lb, ub, args...; kwargs...)
function IntegralProblem(f, lb, ub, args...;
kwargs...)
IntegralProblem{isinplace(f, 3)}(f, lb, ub, args...; kwargs...)
end

Expand All @@ -391,6 +395,56 @@ struct QuadratureProblem end

@doc doc"""

Defines a integral problem over pre-sampled data.
Documentation Page: https://docs.sciml.ai/Integrals/stable/

## Mathematical Specification of a data Integral Problem

Sampled integral problems are defined as:

```math
\sum_i w_i y_i
```
where `y_i` are sampled values of the integrand, and `w_i` are weights
assigned by a quadrature rule, which depend on sampling points `x`.

## Problem Type

### Constructors

```
SampledIntegralProblem(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.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason to leave it up to the user? I think that just adds more work in the implementation without a tangible benefit. dimension should always be the last one as other wise you'll have a performance hit from column major, so we should just stick to that.

Copy link
Contributor Author

@IlianPihlajamaa IlianPihlajamaa Sep 17, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I disagree, with both (1) that the last dim is always faster, and (2) even if it were that we should force all users to use that.

  1. Here is a simple example:
function trapz(y, dim)
          solve(SampledIntegralProblem(y, axes(y,dim); dim=dim), 
           TrapezoidalRule())
end
y = rand(5, 100000)
julia> @btime trapz(trapz($y,2), 1)
 1.165 ms (8 allocations: 400 bytes)
u: 200010.59503616963

julia> @btime trapz(trapz($y,1), 1)
 395.000 μs (9 allocations: 781.59 KiB)
u: 200010.59503616992

while there are better ways to do a two-dimensional trapezoidal rule of course.

  1. Doing numeric integration over some array isn't always the performance bottleneck of the algorithm. Typically, generating the data in the first place (by some experiment perhaps) takes more time. It would be annoying to force users to stick to a certain data layout to be able to do integration over it if it is just a small part of their workflow. This is also why other packages implement this feature, see SciPy or Trapz.jl for example.

Having said that, you are of course right that the last dimension is typically faster and should be the default, so I'll change that.

- kwargs: Keyword arguments copied to the solvers.

### Fields

The fields match the names of the constructor arguments.
"""
struct SampledIntegralProblem{Y, X, D, K} <: AbstractIntegralProblem{false}
y::Y
x::X
dim::D
kwargs::K
@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`"
@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 SampledIntegralProblem 1 4

@doc doc"""

Defines an optimization problem.
Documentation Page: https://docs.sciml.ai/Optimization/stable/API/optimization_problem/

Expand Down
10 changes: 10 additions & 0 deletions test/function_building_error_messages.jl
Original file line number Diff line number Diff line change
Expand Up @@ -453,7 +453,17 @@ NonlinearFunction(nfoop, vjp = nvjp)
intf(u) = 1.0
@test_throws SciMLBase.TooFewArgumentsError IntegralProblem(intf, 0.0, 1.0)
intf(u, p) = 1.0
p = 2.0

IntegralProblem(intf, 0.0, 1.0)
IntegralProblem(intf, 0.0, 1.0, p)
IntegralProblem(intf, [0.0], [1.0])
IntegralProblem(intf, [0.0], [1.0], p)

x = [1.0, 2.0]
y = rand(2, 2)
SampledIntegralProblem(y, x)
SampledIntegralProblem(y, x; dim=2)

# Optimization

Expand Down