Integrals.jl icon indicating copy to clipboard operation
Integrals.jl copied to clipboard

1D example

Open frank-ccc opened this issue 4 years ago • 6 comments

Hello,

great repo! I am trying to add one parameter to the 1D example, see self-explanatory code below

`using Quadrature, Cuba, Cubature, Zygote, FiniteDiff, ForwardDiff using Test

One Dimensional

f(y,p) = sum(sin.(p[1].*y) + p[2]) lb = 1.0 ub = 3.0 p = [2.0,1.0] prob = QuadratureProblem(f,lb,ub,p) sol = solve(prob,QuadGKJL(),reltol=1e-3,abstol=1e-3)

function testf(lb,ub,p) prob = QuadratureProblem(f,lb,ub,p) solve(prob,QuadGKJL(),reltol=1e-3,abstol=1e-3) end

dlb1,dub1,dp1 = Zygote.gradient(testf,lb,ub,p)`

However, I find the following error at the Zygote gradient call:

ERROR: Output should be scalar; gradients are not defined for output u: 1.31184143840124581

Also, is there any way to handle additional input parameters?

Thank you in advance.

frank-ccc avatar Dec 21 '20 07:12 frank-ccc

You need to do solve(prob,QuadGKJL(),reltol=1e-3,abstol=1e-3).u. But it's having an issue that the domain transformation is increasing nout? @agerlach is that supposed to happen?

ChrisRackauckas avatar Dec 21 '20 13:12 ChrisRackauckas

I'll take a look at it. May not be until tomorrow though.

agerlach avatar Dec 21 '20 13:12 agerlach

But it's having an issue that the domain transformation is increasing nout? @agerlach is that supposed to happen?

Yes. B/c the gradient is computed as the integral of the vector of partials. Here length(p)=2, so nout gets changed to 2. QuadGKJL() does not support vector valued integrands. So, @frank-ccc if you swap out QuadGKJL() for any of the other methods it will work as expected, e.q.

function testf(lb,ub,p)
    prob = QuadratureProblem(f,lb,ub,p)
    solve(prob,QuadGKJL(),reltol=1e-3,abstol=1e-3)[1]
end

function testf2(lb,ub,p)
    prob = QuadratureProblem(f,lb,ub,p)
    solve(prob,CubatureJLh(),reltol=1e-3,abstol=1e-3)[1]
end

dp1 = Zygote.gradient(p->testf(lb,ub,p),p)[1]
dp2 = Zygote.gradient(p->testf2(lb,ub,p),p)[1]

agerlach avatar Dec 22 '20 20:12 agerlach

Oh hmm, it would be good to find an informative error message here, or break it apart on QuadGK

ChrisRackauckas avatar Dec 23 '20 00:12 ChrisRackauckas

@ChrisRackauckas @agerlach thanks a lot to help resolve this. I am actually trying to solve a Fredholm problem with differential programming. May I contact you separately if have questions?

frank-ccc avatar Dec 23 '20 02:12 frank-ccc

Yeah sure

ChrisRackauckas avatar Dec 23 '20 02:12 ChrisRackauckas

@ChrisRackauckas It appears this issue was fixed by #189 since we no longer use nout as a check on QuadGKJL. The following translation of the mwe is working for me

using Integrals, Cuba, Cubature, Zygote, FiniteDiff, ForwardDiff
using Test

f(y,p) = sum(sin.(p[1].*y) + p[2])
lb = 1.0
ub = 3.0
p = [2.0,1.0]
prob = IntegralProblem(f,lb,ub,p)
sol = solve(prob,QuadGKJL(),reltol=1e-3,abstol=1e-3)

function testf(lb,ub,p)
    prob = IntegralProblem(f,lb,ub,p)
    solve(prob,QuadGKJL(),reltol=1e-3,abstol=1e-3).u
end

function testf2(lb,ub,p)
    prob = IntegralProblem(f,lb,ub,p)
    solve(prob,CubatureJLh(),reltol=1e-3,abstol=1e-3).u
end

dlb1,dub1,dp1 = Zygote.gradient(testf,lb,ub,p)
dlb2,dub2,dp2 = Zygote.gradient(testf2,lb,ub,p)

@test dlb1 ≈ dlb2
@test dub1 ≈ dub2
@test dp1 ≈ dp2

I will follow up on this with a pr to run the existing AD tests on all the compatible algorithms, which would be helpful at catching edge cases.

lxvm avatar Nov 03 '23 18:11 lxvm

Awesome, really nice job on that update PR. Seems like it fixed most issues 😅

ChrisRackauckas avatar Nov 03 '23 19:11 ChrisRackauckas

Can confirm that this is fixed by #189 and will close now that we have #196 for better testing

lxvm avatar Jan 03 '24 18:01 lxvm