Merge pull request #374 from killah-t-cell/gb-zm-heterogeneous-input

[WIP] ZDM / GB heterogeneous input

Author: ChrisRackauckas

.DS_Store

@@ -62,7 +62,7 @@ chain = FastChain(FastDense(dim,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(1
 
 discretization = PhysicsInformedNN(chain, QuadratureTraining())
 
-@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u])
+@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
 prob = discretize(pde_system,discretization)
 
 cb = function (p,l)

README.md

@@ -74,7 +74,7 @@ discretization = PhysicsInformedNN(chain,
                                    strategy;
                                    init_params = initθ)
 
-@named pde_system = PDESystem(eq,bcs,domains,[t,x,y],[u])
+@named pde_system = PDESystem(eq,bcs,domains,[t,x,y],[u(t, x, y)])
 prob = discretize(pde_system,discretization)
 symprob = symbolic_discretize(pde_system,discretization)
 
@@ -126,8 +126,8 @@ plot_(res)
 
 ![3pde](https://user-images.githubusercontent.com/12683885/129949743-9471d230-c14f-4105-945f-6bc52677d40e.gif)
 
-
 ## Performance benchmarks
+
 Here are some performance benchmarks for 2d-pde with various number of input points and the number of neurons in the hidden layer, measuring the time for 100 iterations. Сomparing runtime with GPU and CPU.
 
 ```julia

docs/src/pinn/2D.md

@@ -38,7 +38,7 @@ domains = [x ∈ Interval(0.0,1.0)]
 chain = FastChain(FastDense(1,8,Flux.σ),FastDense(8,1))
 
 discretization = PhysicsInformedNN(chain, QuasiRandomTraining(20))
-@named pde_system = PDESystem(eq,bcs,domains,[x],[u])
+@named pde_system = PDESystem(eq,bcs,domains,[x],[u(x)])
 prob = discretize(pde_system,discretization)
 
 cb = function (p,l)
@@ -66,4 +66,5 @@ x_plot = collect(xs)
 plot(x_plot ,u_real,title = "real")
 plot!(x_plot ,u_predict,title = "predict")
 ```
+
 ![hodeplot](https://user-images.githubusercontent.com/12683885/90276340-69bc3e00-de6c-11ea-89a7-7d291123a38b.png)

docs/src/pinn/3rd.md

@@ -49,7 +49,7 @@ isapprox(dphi[1][2], dphi2, atol=1e-8)
 
 
 indvars = [x,t]
-depvars = [u]
+depvars = [u(x, t)]
 dim = length(domains)
 dx = 0.1
 strategy = NeuralPDE.GridTraining(dx)

docs/src/pinn/debugging.md

@@ -69,7 +69,7 @@ discretization = PhysicsInformedNN(chain,
                                    init_params = initθ,
                                    additional_loss=norm_loss_function)
 
-@named pde_system = PDESystem(eq,bcs,domains,[x],[p])
+@named pde_system = PDESystem(eq,bcs,domains,[x],[p(x)])
 prob = discretize(pde_system,discretization)
 
 pde_inner_loss_functions = prob.f.f.loss_function.pde_loss_function.pde_loss_functions.contents

docs/src/pinn/fp.md

@@ -0,0 +1,37 @@
+# Differential Equations with Heterogeneous Inputs
+
+A differential equation is said to have heterogeneous inputs when its dependent variables depend on different independent variables:
+
+```math
+u(x) + w(x, v) = \frac{\partial w(x, v)}{\partial w}
+```
+
+Here, we write an arbitrary heterogeneous system:
+
+```julia
+@parameters x y
+@variables p(..) q(..) r(..) s(..)
+Dx = Differential(x)
+Dy = Differential(y)
+
+# 2D PDE
+eq  = p(x) + q(y) + Dx(r(x, y)) + Dy(s(y, x)) ~ 0
+
+# Initial and boundary conditions
+bcs = [p(1) ~ 0.f0, q(-1) ~ 0.0f0,
+       r(x, -1) ~ 0.f0, r(1, y) ~ 0.0f0,
+       s(y, 1) ~ 0.0f0, s(-1, x) ~ 0.0f0]
+
+# Space and time domains
+domains = [x ∈ Interval(0.0, 1.0),
+           y ∈ Interval(0.0, 1.0)]
+
+numhid = 3
+fastchains = [[FastChain(FastDense(1, numhid, Flux.σ), FastDense(numhid, numhid, Flux.σ), FastDense(numhid, 1)) for i in 1:2];
+                        [FastChain(FastDense(2, numhid, Flux.σ), FastDense(numhid, numhid, Flux.σ), FastDense(numhid, 1)) for i in 1:2]]
+discretization = NeuralPDE.PhysicsInformedNN(fastchains, QuadratureTraining()
+
+@named pde_system = PDESystem(eq, bcs, domains, [x,y], [p(x), q(y), r(x, y), s(y, x)])
+prob = SciMLBase.discretize(pde_system, discretization)
+res = GalacticOptim.solve(prob, BFGS(); cb=cb, maxiters=100)
+```

docs/src/pinn/heterogeneous.md

@@ -1,4 +1,3 @@
-
 # Integro Differential Equations
 
 The integral of function u(x),
@@ -31,14 +30,16 @@ Similarly a rectangular or cuboidal domain can be defined using `ProductDomain`
 Ix = Integral((x,y) in DomainSets.ProductDomain(ClosedInterval(0 ,1), ClosedInterval(0 ,x)))
 ```
 
-
 ## 1-dimensional example
 
 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
 ```
+
 and boundary condition
+
 ```math
 u(0) = 0
 ```
@@ -63,7 +64,7 @@ discretization = PhysicsInformedNN(chain,
                                    init_params = nothing,
                                    phi = nothing,
                                    derivative = nothing)
-pde_system = PDESystem(eq,bcs,domains,[t],[i])
+pde_system = PDESystem(eq,bcs,domains,[t],[i(t)])
 prob = NeuralPDE.discretize(pde_system,discretization)
 cb = function (p,l)
     println("Current loss is: $l")

docs/src/pinn/integro_diff.md

@@ -6,13 +6,13 @@ Let's consider the Kuramoto–Sivashinsky equation, which contains a 4th-order d
 ∂_t u(x, t) + u(x, t) ∂_x u(x, t) + \alpha ∂^2_x u(x, t) + \beta ∂^3_x u(x, t) + \gamma ∂^4_x u(x, t) =  0 \, ,
 ```
 
-where ``\alpha = \gamma = 1`` and ``\beta = 4``. The exact solution is:
+where `\alpha = \gamma = 1` and `\beta = 4`. The exact solution is:
 
 ```math
 u_e(x, t) = 11 + 15 \tanh \theta - 15 \tanh^2 \theta - 15 \tanh^3 \theta \, ,
 ```
 
-where ``\theta = 1 - x/2`` and with initial and boundary conditions:
+where `\theta = 1 - x/2` and with initial and boundary conditions:
 
 ```math
 \begin{align*}
@@ -62,7 +62,7 @@ dx = 0.4; dt = 0.2
 chain = FastChain(FastDense(2,12,Flux.σ),FastDense(12,12,Flux.σ),FastDense(12,1))
 
 discretization = PhysicsInformedNN(chain, GridTraining([dx,dt]))
-@named pde_system = PDESystem(eq,bcs,domains,[x,t],[u])
+@named pde_system = PDESystem(eq,bcs,domains,[x,t],[u(x, t)])
 prob = discretize(pde_system,discretization)
 
 cb = function (p,l)

docs/src/pinn/ks.md

@@ -44,7 +44,7 @@ discretization = NeuralPDE.PhysicsInformedNN(chain1,
                                              quadrature_strategy;
                                              init_params = initθ)
 
-@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u])
+@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
 prob = NeuralPDE.discretize(pde_system,discretization)
 sym_prob = NeuralPDE.symbolic_discretize(pde_system,discretization)
 
@@ -95,11 +95,12 @@ p5 = plot(xs, ys, diff_u_,linetype=:contourf,title = "error 2");
 plot(p1,p2,p3,p4,p5)
 
 ```
+
 ![neural_adapter](https://user-images.githubusercontent.com/12683885/127645487-24226cc0-fff6-4bb9-8509-77cf3e120563.png)
 
 ## Domain decomposition
 
-In this example, we first obtain a prediction of 2D Poisson equation on subdomains. We split up full domain into 10 sub problems by x, and  create separate neural networks for each sub interval. If x domain ∈ [x_0, x_end] so, it is decomposed on 10 part: sub x domains = {[x_0, x_1], ... [x_i,x_i+1], ..., x_9,x_end]}.
+In this example, we first obtain a prediction of 2D Poisson equation on subdomains. We split up full domain into 10 sub problems by x, and create separate neural networks for each sub interval. If x domain ∈ [x_0, x_end] so, it is decomposed on 10 part: sub x domains = {[x_0, x_1], ... [x_i,x_i+1], ..., x_9,x_end]}.
 And then using the method neural_adapter, we retrain the banch of 10 predictions to the one prediction for full domain of task.
 
 ![domain_decomposition](https://user-images.githubusercontent.com/12683885/127149752-a4ecea50-2984-45d8-b0d4-d2eadecf58e7.png)
@@ -171,7 +172,7 @@ for i in 1:count_decomp
     @register phi_bound(x,y)
     Base.Broadcast.broadcasted(::typeof(phi_bound), x,y) = phi_bound(x,y)
     bcs_ = create_bcs(bcs,domains_[1].domain, phi_bound)
-    @named pde_system_ = PDESystem(eq, bcs_, domains_, [x, y], [u])
+    @named pde_system_ = PDESystem(eq, bcs_, domains_, [x, y], [u(x, y)])
     push!(pde_system_map,pde_system_)
     strategy = NeuralPDE.GridTraining([0.1/count_decomp, 0.1])
 
@@ -225,7 +226,7 @@ chain2 = FastChain(FastDense(2,inner_,af),
 
 initθ2 =Float64.(DiffEqFlux.initial_params(chain2))
 
-@named pde_system = PDESystem(eq, bcs, domains, [x, y], [u])
+@named pde_system = PDESystem(eq, bcs, domains, [x, y], [u(x, y)])
 
 losses = map(1:count_decomp) do i
     loss(cord,θ) = chain2(cord,θ) .- phis[i](cord,reses[i].minimizer)

docs/src/pinn/neural_adapter.md

@@ -1,5 +1,6 @@
 # Optimising Parameters of a Lorenz System
- Consider a Lorenz System ,
+
+Consider a Lorenz System ,
 
 ```math
 \begin{align*}
@@ -9,7 +10,7 @@
 \end{align*}
 ```
 
-with Physics-Informed Neural Networks. Now we would consider the case where we want to optimise the parameters ``\sigma``, ``\beta``,  and ``\rho``.
+with Physics-Informed Neural Networks. Now we would consider the case where we want to optimise the parameters `\sigma`, `\beta`, and `\rho`.
 
 We start by defining the the problem,
 
@@ -27,14 +28,17 @@ bcs = [x(0) ~ 1.0, y(0) ~ 0.0, z(0) ~ 0.0]
 domains = [t ∈ Interval(0.0,1.0)]
 dt = 0.01
 ```
+
 And the neural networks as,
+
 ```julia
 input_ = length(domains)
 n = 8
 chain1 = FastChain(FastDense(input_,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,1))
 chain2 = FastChain(FastDense(input_,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,1))
 chain3 = FastChain(FastDense(input_,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,n,Flux.σ),FastDense(n,1))
 ```
+
 We will add an additional loss term based on the data that we have in order to optimise the parameters.
 
 Here we simply calculate the solution of the lorenz system with [OrdinaryDiffEq.jl](https://diffeq.sciml.ai/v1.10/tutorials/ode_example.html#In-Place-Updates-1) based on the adaptivity of the ODE solver. This is used to introduce non-uniformity to the time series.
@@ -64,17 +68,20 @@ sep = [acum[i]+1 : acum[i+1] for i in 1:length(acum)-1]
 (u_ , t_) = data
 len = length(data[2])
 ```
+
 Then we define the additional loss funciton `additional_loss(phi, θ , p)`, the function has three arguments, `phi` the trial solution, `θ` the parameters of neural networks, and the hyperparameters `p` .
 
 ```julia
 function additional_loss(phi, θ , p)
     return sum(sum(abs2, phi[i](t_ , θ[sep[i]]) .- u_[[i], :])/len for i in 1:1:3)
 end
 ```
+
 Then finally defining and optimising using the `PhysicsInformedNN` interface.
+
 ```julia
 discretization = NeuralPDE.PhysicsInformedNN([chain1 , chain2, chain3],NeuralPDE.GridTraining(dt), param_estim=true, additional_loss=additional_loss)
-@named pde_system = PDESystem(eqs,bcs,domains,[t],[x, y, z],[σ_, ρ, β], defaults=Dict([p .=> 1.0 for p in [σ_, ρ, β]]))
+@named pde_system = PDESystem(eqs,bcs,domains,[t],[x(t), y(t), z(t)],[σ_, ρ, β], defaults=Dict([p .=> 1.0 for p in [σ_, ρ, β]]))
 prob = NeuralPDE.discretize(pde_system,discretization)
 cb = function (p,l)
     println("Current loss is: $l")
@@ -83,7 +90,9 @@ end
 res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=5000)
 p_ = res.minimizer[end-2:end] # p_ = [9.93, 28.002, 2.667]
 ```
+
 And then finally some analyisis by plotting.
+
 ```julia
 initθ = discretization.init_params
 acum =  [0;accumulate(+, length.(initθ))]

docs/src/pinn/parm_estim.md

@@ -54,7 +54,7 @@ chain = FastChain(FastDense(dim,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(1
 dx = 0.05
 discretization = PhysicsInformedNN(chain,GridTraining(dx))
 
-@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u])
+@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
 prob = discretize(pde_system,discretization)
 
 #Optimizer
@@ -117,6 +117,7 @@ chain = FastChain(FastDense(dim,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(1
 ```
 
 Here, we build PhysicsInformedNN algorithm where `dx` is the step of discretization and `strategy` stores information for choosing a training strategy.
+
 ```julia
 # Discretization
 dx = 0.05
@@ -126,7 +127,7 @@ discretization = PhysicsInformedNN(chain, GridTraining(dx))
 As described in the API docs, we now need to define the `PDESystem` and create PINNs problem using the `discretize` method.
 
 ```julia
-@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u])
+@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
 prob = discretize(pde_system,discretization)
 ```
 
@@ -161,4 +162,5 @@ p2 = plot(xs, ys, u_predict, linetype=:contourf,title = "predict");
 p3 = plot(xs, ys, diff_u,linetype=:contourf,title = "error");
 plot(p1,p2,p3)
 ```
+
 ![poissonplot](https://user-images.githubusercontent.com/12683885/90962648-2db35980-e4ba-11ea-8e58-f4f07c77bcb9.png)