Add Integro Differential Equations Support #330

ashutosh-b-b · 2021-06-23T18:39:10Z

How do I add support for closure inside the get_numeric_integral? get_numeric_integral cannot be a fixed function since the number of dependent variables may vary.

KirillZubov · 2021-06-23T20:49:24Z

why can't? , probably better name parse_integral(or something like this) then get_numeric_integral . You need to call this function for every symbol form of integral in the model and pass dependent variables for every integral.

KirillZubov · 2021-06-23T20:54:24Z

my pseudo code for this :)

eqs = integral1(depvars1,...) + integral2(depvars2,...)+ somthing ->
inner_representation_integral1(depvars1) + inner_representation_integral2(depvars2)+ somthing

parser(integral1(depvars1, ...),...) -> inner_representation_integral1

KirillZubov · 2021-07-06T17:27:56Z

src/pinns_pde_solve.jl

                break
+            elseif e isa Symbolics.Integral
+                integrating_variable = toexpr(_args[1].x)
+                integrand = transform_expression(_args[2],dict_indvars,dict_depvars,chain,eltypeθ,strategy)


_transform_expression

KirillZubov · 2021-07-06T18:54:00Z

src/pinns_pde_solve.jl

+    integral =
+        (u, cord, phi, θ, integrating_variable, integrand, lb, ub, strategy=strategy, indvars=indvars, depvars=depvars)->
+            begin
+                integrand_ = NeuralPDE.build_symbolic_loss_function(nothing, indvars, depvars, phi, derivative, nothing, chain, θ, strategy, integrand = integrand)


all actions with the parser must be performed before starting training. To call the build_symbolic_loss in integral function means that it will be called every time during training, so it will take a very long time.

KirillZubov · 2021-07-12T09:04:18Z

src/pinns_pde_solve.jl

    phi = discretization.phi
    derivative = discretization.derivative
    strategy = discretization.strategy
+    global integral = get_numeric_integral(strategy, pde_system.indvars, pde_system.depvars, dict_indvars, chain, derivative)


global is bad. it will always lead to trouble. Avoid it

in this case I needed to use integral outside the symbolic_discretize and discretize. And hence global declaration. What are other ways to do it?

KirillZubov · 2021-07-19T11:51:51Z

test/NNPDE_tests.jl

+eq = Di(i(t)) + 2*i(t) + 5*Ii(i(t)) ~ 1
+bcs = [i(0.) ~ 0.0]
+domains = [t ∈ Interval(0.0,2.0)]
+chain = FastChain(FastDense(1,15,Flux.σ),FastDense(15,1))


it's some trouble with FastChain after the new release of Flux, just use Chain for a test while FastChain won't have been fixed.

initial_params change behavior after the update in Flux and now it works incorrectly.

chain = FastChain(FastDense(1,12,Flux.σ),FastDense(12,1)) julia> initθ = DiffEqFlux.initial_params(chain) 26-element Vector{Union{Nothing, Float32}}: 0.44163412f0 0.3313627f0 -0.24543865f0 0.16519924f0 0.18937507f0 -0.047894925f0 ⋮ -0.50703824f0 -0.13543157f0 0.036089323f0 0.0010503983f0 nothing

Can you try and narrow that down?

ok, I will open the issue in DiffEqFlux as soon as I figure out what exactly is the problem with DiffEqFlux.initial_params.

KirillZubov · 2021-07-20T18:16:02Z

test/NNPDE_tests.jl

+                                             derivative = nothing,
+                                             )
+pde_system = PDESystem(eq,bcs,domains,[t],[i])
+prob = NeuralPDE.discretize(pde_system,discretization)


it work in forward but

prob.f.f.loss_function(initθ) 1.1225828958232524

donesn't work in backward

julia> Zygote.gradient(prob.f.f.loss_function,initθ) ERROR: MethodError: no method matching get(::Tuple{Int64, Int64}, ::Int64, ::Int64)

I suppose it is related to integral function.
So you need that work in both way:
integral = get_numeric_integral(strategy, _indvars, _depvars, dict_indvars, chain, derivative)

integral(initθ,... ) integral_ = (θ) -> integral(θ, ...) Zygote.gradient(integral_,initθ)

All symbolic things should do before calculation because Zygote.gradient probably can't handle all of it.

ChrisRackauckas · 2021-07-26T12:33:44Z

Is it still WIP?

ChrisRackauckas · 2021-07-26T12:43:15Z

test/NNPDE_tests.jl

+## Simple Integral Test
+@parameters x
+@variables u(..)
+Ix = Integral(x, DomainSets.ClosedInterval(0, x))


needs a multidimensional test

Yes, regarding that I wanted to ask this, for a multidimensional integral, I wanted to change definition of the symbolic integral to have a vector of symbols. For eg. I wanted it to be like this : I = Integral([x , y, z], UnitCube()). Should I proceed?

It should be done in the same way we're describing sets, i.e. (x , y, z) \in UnitCube()

I am talking about the Integral.

yes. I = Integral(x^2 + y^2 + z^2,(x,y,z) in UnitCube())

We first define the Integral as Integral(integrating_variable, domain) and then pass the function. What I meant was using vector of integrating_variables instead. Is it fine? Or do I have to change it?

KirillZubov · 2021-07-27T13:41:26Z

I think it will need to split tests files and move tests with integral to integro_differential_test.jl

KirillZubov · 2021-08-09T12:50:06Z

Project.toml

 Distributions = "0.23, 0.24, 0.25"
 DocStringExtensions = "0.8"
-Flux = "0.10.1, 0.11, 0.12"
+Flux = "0.10.1, 0.11, < 0.12.5"


"Resolve conflicts" to Flux = "0.12.6"

just press the button “resolve conflict” and merge it

KirillZubov · 2021-08-09T13:46:06Z

src/pinns_pde_solve.jl

-
-    loss_function = parse_equation(eqs,dict_indvars,dict_depvars,chain,eltypeθ,strategy)
-    vars = :(cord, $θ, phi, derivative,u,p)
+    if integrand isa Nothing


why is it here?

So we can re use the build_symbolic_loss_function for getting the integral function in our inner representation.

KirillZubov · 2021-08-10T16:45:23Z

src/pinns_pde_solve.jl

    phi = discretization.phi
    derivative = discretization.derivative
    strategy = discretization.strategy
+    global integral = get_numeric_integral(strategy, pde_system.indvars, pde_system.depvars, chain, derivative)


why isn't it removed global yet

KirillZubov · 2021-08-10T16:55:01Z

src/pinns_pde_solve.jl

+function get_numeric_integral(strategy, _indvars, _depvars, chain, derivative)
+    depvars,indvars,dict_indvars,dict_depvars = get_vars(_indvars, _depvars)
+    integral =
+        (u, cord, phi, integrating_var_id, integrand_func, lb, ub, θ ;strategy=strategy, indvars=indvars, depvars=depvars, dict_indvars=dict_indvars, dict_depvars=dict_depvars)->


these parameters can be remove as no use indvars=indvars, depvars=depvars, dict_indvars=dict_indvars, dict_depvars=dict_depvars

These are used inside the integral function, so shouldn't we pass them?

we don't use these for initializing function inside get_numeric_integral so we don't need to pass them.

KirillZubov · 2021-08-10T17:57:47Z

src/pinns_pde_solve.jl

+    integral =
+        (u, cord, phi, integrating_var_id, integrand_func, lb, ub, θ ;strategy=strategy, indvars=indvars, depvars=depvars, dict_indvars=dict_indvars, dict_depvars=dict_depvars)->
+            begin
+                flat_θ = if (typeof(chain) <: AbstractVector) reduce(vcat,θ) else θ end


why do you need it? θ should already be flat

I checked it, it wasnt. I ll recheck

KirillZubov · 2021-08-10T18:27:54Z

src/pinns_pde_solve.jl

+                        ub_ = repeat(ub, 1, size(cord)[2])
+                    end
+                end
+                for i in 1:size(cord)[2]


if cord = ones(2,300)
so then
julia> size(cord)[2]
300
why we calculate it for each index? it is very slow

in order to give the vectorise output. We have each index is a cord point in the data, so we need to calculate integral seperately for each point

I did try mapslices but it threw a Mutating arrays not supported error

it is a vector we do not need an index, we need a broadcast here for fast performance (see how it works in get_numeric_derivative).
In any case, it works fully correctly, but slow

using BenchmarkTools f(x) =x^2 function f_index(x) for i in 1:size(x)[2] x[:,i].^2 end end julia> @benchmark f.(data) setup=(data=rand(2,1000)) BenchmarkTools.Trial: 10000 samples with 9 evaluations. Range (min … max): 2.525 μs … 3.818 ms ┊ GC (min … max): 0.00% … 99.82% Time (median): 3.732 μs ┊ GC (median): 0.00% Time (mean ± σ): 7.412 μs ± 85.571 μs ┊ GC (mean ± σ): 28.22% ± 2.44% ▁ ▂ █▅▄▇▂ ▁ ▃▄▅▄▄▄▃▃▃▃▂▂▂▁▁ ▂ █▄█████████▆▅▄▅▄▅▆▅▅▅▃▅▆▅▅▅▇█████████████████▇▇▇▇▆▄▅▅▆▅▅▅▅ █ 2.52 μs Histogram: log(frequency) by time 12.2 μs < Memory estimate: 15.75 KiB, allocs estimate: 1. julia> @benchmark f_index(data) setup=(data=rand(2,1000)) BenchmarkTools.Trial: 10000 samples with 1 evaluation. Range (min … max): 56.758 μs … 36.600 ms ┊ GC (min … max): 0.00% … 99.53% Time (median): 78.680 μs ┊ GC (median): 0.00% Time (mean ± σ): 118.139 μs ± 1.100 ms ┊ GC (mean ± σ): 29.35% ± 3.15% ▄▃ ▆█▆▆▆▄▅█▇▃▃▂▂▃▂▃▃▃▂▂▂▂▁▁▁▁▁▁▁▁▁ ▁▁ ▂ ████████████████████████████████████████▇▇█▇▇▇▇▇▇▆▆▆▇▆▆▄▅▄▅▅ █ 56.8 μs Histogram: log(frequency) by time 188 μs < Memory estimate: 187.50 KiB, allocs estimate: 2000.

ChrisRackauckas · 2021-08-15T22:32:06Z

Project.toml

 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"


why is this added?

Oh that got pushed by mistake. I ll remove it.

KirillZubov · 2021-08-16T20:12:43Z

test/NNPDE_tests.jl

+                                             phi = nothing,
+                                             derivative = nothing,
+                                             )
+pde_system = PDESystem(eq,bcs,domains,[x,y],[u])


@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u])

KirillZubov · 2021-08-16T20:13:33Z

test/NNPDE_tests.jl

+                                             phi = nothing,
+                                             derivative = nothing,
+                                             )
+pde_system = PDESystem(eq,bcs,domains,[x,y],[u])


KirillZubov · 2021-08-16T20:13:45Z

test/NNPDE_tests.jl

+                                             phi = nothing,
+                                             derivative = nothing,
+                                             )
+pde_system = PDESystem(eq,bcs,domains,[x],[u])


KirillZubov · 2021-08-16T20:14:37Z

docs/src/pinn/integro_diff.md

+                                             phi = nothing,
+                                             derivative = nothing,
+                                             )
+pde_system = PDESystem(eq,bcs,domains,[t],[i])


KirillZubov · 2021-08-16T20:16:42Z

test/NNPDE_tests.jl

+                                             phi = nothing,
+                                             derivative = nothing,
+                                             )
+pde_system = PDESystem(eq,bcs,domains,[t],[i])


and here @named

KirillZubov · 2021-08-16T20:22:20Z

docs/src/pinn/integro_diff.md

+We use the grid training strategy and then use the `PhysicsInformedNN` interface to define the optimization problem.
+```julia
+strategy_ = NeuralPDE.GridTraining(0.05)
+discretization = NeuralPDE.PhysicsInformedNN(chain,


remove all NeuralPDE. from doc

KirillZubov · 2021-08-16T20:51:18Z

docs/src/pinn/integro_diff.md

+# Integro Differential Equations
+Lets take an example of an integro differential equation :
+```math
+\frac{∂}{∂x} u(x)  + 2u(x) + 5 \int_{0}{x}u(t)dt = 1 for x \geq 0


\int_{0}^{x}u(t)dt

KirillZubov · 2021-08-16T22:19:33Z

docs/src/pinn/integro_diff.md

+                                             derivative = nothing,
+                                             )
+pde_system = PDESystem(eq,bcs,domains,[t],[i])
+prob = NeuralPDE.discretize(pde_system,discretization)


add

cb = function (p,l) println("Current loss is: $l") return false end

ChrisRackauckas · 2021-08-22T13:30:49Z

test/NNPDE_tests.jl

+@variables u(..)
+Dx = Differential(x)
+Dy = Differential(y)
+Ix = Integral([x,y], (x,y) in DomainSets.ProductDomain(UnitInterval(),ClosedInterval(0 ,x)))


why is this two outputs? The equations below don't have the right size.

ChrisRackauckas · 2021-08-22T13:31:04Z

test/NNPDE_tests.jl

+Ix = Integral([x,y], (x,y) in DomainSets.UnitSquare())
+eq = Ix(u(x,y)) ~ 1/3
+bcs = [u(0., 0.) ~ 1, Dx(u(x,y)) ~ -2*x , Dy(u(x ,y)) ~ -2*y ]


same sizing issue.

ChrisRackauckas · 2021-08-23T11:12:50Z

@sdwfrost

ashutosh-b-b force-pushed the Integro_diff branch from 6988d52 to 0416373 Compare July 2, 2021 19:04

ashutosh-b-b mentioned this pull request Jul 5, 2021

Volterra integral equations SciML/DifferentialEquations.jl#760

Open

KirillZubov reviewed Jul 6, 2021

View reviewed changes

ashutosh-b-b force-pushed the Integro_diff branch from 0416373 to e481ca6 Compare July 11, 2021 13:18

KirillZubov reviewed Jul 12, 2021

View reviewed changes

ashutosh-b-b requested a review from KirillZubov July 18, 2021 17:20

KirillZubov reviewed Jul 19, 2021

View reviewed changes

KirillZubov reviewed Jul 20, 2021

View reviewed changes

ashutosh-b-b requested review from ChrisRackauckas and KirillZubov July 26, 2021 12:30

ChrisRackauckas reviewed Jul 26, 2021

View reviewed changes

KirillZubov reviewed Aug 9, 2021

View reviewed changes

ashutosh-b-b requested a review from KirillZubov August 10, 2021 15:00

KirillZubov reviewed Aug 10, 2021

View reviewed changes

ChrisRackauckas reviewed Aug 15, 2021

View reviewed changes

ashutosh-b-b requested a review from KirillZubov August 16, 2021 17:14

KirillZubov reviewed Aug 16, 2021

View reviewed changes

docs/src/pinn/integro_diff.md

phi = nothing,

derivative = nothing,

)

pde_system = PDESystem(eq,bcs,domains,[t],[i])

Copy link

Member

KirillZubov Aug 16, 2021

Choose a reason for hiding this comment

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

@named

KirillZubov reviewed Aug 16, 2021

View reviewed changes

ashutosh-b-b added 13 commits August 17, 2021 16:41

update Project.toml

74dc2c4

Remove dict_indvars from parameters

501f6ba

Fix test

705c85b

add multidim integral

cd4f969

Add simple 2d test for integral

eb1901b

update Project.toml

91b03ca

Remove global integral

1cb3d8d

Change for symbolic multi integral support

7f4dbd0

Update Tests

6bec680

Add docs

baeec1a

Update Project.toml

fa86808

Fix after rebase

ed96585

fix tests

d4f71af

ashutosh-b-b force-pushed the Integro_diff branch from 73190de to d4f71af Compare August 17, 2021 11:15

ashutosh-b-b and others added 6 commits August 18, 2021 18:53

Edit Docs

abc9659

Merge branch 'master' into Integro_diff

68b591c

Update pinns_pde_solve.jl

a3ab03b

fix

61310cb

fix test

e34b670

small update tests

532e559

ashutosh-b-b changed the title ~~[WIP] Add Integro Differential Equations Support~~ Add Integro Differential Equations Support Aug 21, 2021

ChrisRackauckas reviewed Aug 22, 2021

View reviewed changes

Shift integrating_var from x to VarDomainPairing

a5ab956

ChrisRackauckas approved these changes Aug 23, 2021

View reviewed changes

ChrisRackauckas closed this Aug 23, 2021

ChrisRackauckas reopened this Aug 23, 2021

ChrisRackauckas merged commit 5114808 into SciML:master Aug 23, 2021

KirillZubov mentioned this pull request Sep 8, 2021

Broadcast integral calculation and increase performance #390

Closed

Uh oh!

Add Integro Differential Equations Support #330

Add Integro Differential Equations Support #330

Uh oh!

Conversation

ashutosh-b-b commented Jun 23, 2021

Uh oh!

KirillZubov commented Jun 23, 2021

Uh oh!

KirillZubov commented Jun 23, 2021

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ChrisRackauckas commented Jul 26, 2021

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KirillZubov commented Jul 27, 2021

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KirillZubov Aug 10, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

KirillZubov Aug 10, 2021 •

edited

Loading

KirillZubov Aug 16, 2021 •

edited

Loading

KirillZubov Aug 16, 2021 •

edited

Loading

KirillZubov Aug 16, 2021 •

edited

Loading