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

1D example #57

Closed
frank-ccc opened this issue Dec 21, 2020 · 9 comments
Closed

1D example #57

frank-ccc opened this issue Dec 21, 2020 · 9 comments

Comments

@frank-ccc
Copy link

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.

@ChrisRackauckas
Copy link
Member

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?

@agerlach
Copy link
Collaborator

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

@agerlach
Copy link
Collaborator

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]

@ChrisRackauckas
Copy link
Member

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

@frank-ccc
Copy link
Author

@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?

@ChrisRackauckas
Copy link
Member

Yeah sure

@lxvm
Copy link
Collaborator

lxvm commented Nov 3, 2023

@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.

@ChrisRackauckas
Copy link
Member

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

@lxvm
Copy link
Collaborator

lxvm commented Jan 3, 2024

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

@lxvm lxvm closed this as completed Jan 3, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants