Defining an array of symbolic variables causes error with NeuralPDE.discretize()

[hpieper14]

I am trying to solve a PDE with variables [t, z1, …, zn]. Here’s my code:

using Test, Flux, Optim, DiffEqFlux, Optimization
using Random, NeuralPDE, DifferentialEquations
using Statistics, Distributions, LinearAlgebra
import ModelingToolkit: Interval 
import DomainSets: UnitInterval
Random.seed!(100)

order = 3

@parameters t z[1:order]
@variables u(..)
Dt = Differential(t)
z = collect(z)

α = 1.2
β = 1.1

function dW(t, z...)
    val = 0 
    for i in 1:length(z)
        val += z[i]*cos((i-1/2)π*t)
    end
    val = √2*val
end

eq  = Dt(u(t,z...))  ~ α*u(t,z...) + β*u(t,z)*dW(t,z...)
bcs = [u(0, z...) ~ 1.0]

# Space and time domains
domains = Vector{Symbolics.VarDomainPairing}(undef, order + 1)
domains[1] = t ∈ Interval(0.0, 1.0)
for i in 1:order
    domains[i+1] = z[i] ∈ Interval(0.0, 1.0)
end

# number of dimensions
dim = order + 1
chain = Flux.Chain(Dense(dim,16,Flux.σ),Dense(16,16,Flux.σ),Dense(16,1))
# Initial parameters of Neural network
initθ = Float64.(DiffEqFlux.initial_params(chain))

# Discretization
dx = 0.05
discretization = PhysicsInformedNN(chain,GridTraining(dx),init_params =initθ)

@named pde_system = PDESystem(eq,bcs,domains,[t,z...],[u(t,z...)])
prob = discretize(pde_system,discretization)

This gives the error

ERROR: MethodError: no method matching nameof(::Term{Real, Base.ImmutableDict{DataType, Any}})
Closest candidates are:
  nameof(::Sym) at ~/.julia/packages/SymbolicUtils/vnuIf/src/types.jl:144
  nameof(::ModelingToolkit.AbstractSystem) at ~/.julia/packages/ModelingToolkit/iHLWM/src/systems/abstractsystem.jl:139
  nameof(::DataType) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/base/reflection.jl:223
  ...
Stacktrace:
 [1] (::NeuralPDE.var"#131#132")(argument::Term{Real, Base.ImmutableDict{DataType, Any}})
   @ NeuralPDE ./none:0
 [2] iterate
   @ ./generator.jl:47 [inlined]
 [3] collect_to!(dest::Vector{Symbol}, itr::Base.Generator{Vector{SymbolicUtils.Symbolic{Real}}, NeuralPDE.var"#131#132"}, offs::Int64, st::Int64)
   @ Base ./array.jl:782
 [4] collect_to_with_first!
   @ ./array.jl:760 [inlined]
 [5] collect(itr::Base.Generator{Vector{SymbolicUtils.Symbolic{Real}}, NeuralPDE.var"#131#132"})
   @ Base ./array.jl:734
 [6] get_vars(indvars_::Vector{Num}, depvars_::Vector{Num})
   @ NeuralPDE ~/Documents/Julia-fork/NeuralPDE.jl/src/pinns_pde_solve.jl:757
 [7] discretize_inner_functions(pde_system::PDESystem, discretization::PhysicsInformedNN{GridTraining{Float64}, Vector{Float64}, NeuralPDE.Phi{Optimisers.Restructure{Chain{Tuple{Dense{typeof(σ), Matrix{Float32}, Vector{Float32}}, Dense{typeof(σ), Matrix{Float32}, Vector{Float32}}, Dense{typeof(identity), Matrix{Float32}, Vector{Float32}}}}, NamedTuple{(:layers,), Tuple{Tuple{NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}, NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}, NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}}}}}}, typeof(NeuralPDE.numeric_derivative), Bool, Nothing, NonAdaptiveLoss{Float64}, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}})
   @ NeuralPDE ~/Documents/Julia-fork/NeuralPDE.jl/src/pinns_pde_solve.jl:1257
 [8] discretize(pde_system::PDESystem, discretization::PhysicsInformedNN{GridTraining{Float64}, Vector{Float64}, NeuralPDE.Phi{Optimisers.Restructure{Chain{Tuple{Dense{typeof(σ), Matrix{Float32}, Vector{Float32}}, Dense{typeof(σ), Matrix{Float32}, Vector{Float32}}, Dense{typeof(identity), Matrix{Float32}, Vector{Float32}}}}, NamedTuple{(:layers,), Tuple{Tuple{NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}, NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}, NamedTuple{(:weight, :bias, :σ), Tuple{Int64, Int64, Tuple{}}}}}}}}, typeof(NeuralPDE.numeric_derivative), Bool, Nothing, NonAdaptiveLoss{Float64}, Nothing, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}})
   @ NeuralPDE ~/Documents/Julia-fork/NeuralPDE.jl/src/pinns_pde_solve.jl:1512
 [9] top-level scope
   @ ~/Documents/Julia-fork/NeuralPDE.jl/test/NNSDE_wp_tests.jl:47

@ChrisRackauckas says that going through the parser and expanding z is a possible fix.

[hpieper14]

@ChrisRackauckas Can I assign this to myself? Unless it makes more sense for someone else to do this fix, I would like to do it.

[ChrisRackauckas]

Yup assigned. Let me know if you run into any issues: I think this might be a hard one but I think you can do it.

[hpieper14]

I've been reading through the code to map out the necessary changes and I want to make sure I understand the big picture before I go in and start changing things.

We want to extend the discretization framework so that it can work with PDESystems that have splatted vector variables, i.e. pde_system = PDESystem(eq, bcs, domains, [x,y...], [u(x,y...)]). It looks like to me that the changes for this need to occur in the files symbolic_utilities.jl and discretize.jl. Symbolic_utilities.jl is where the current error is being thrown and as suggested by name, has symbolic parsing helper functions and if discretize.jl needs to be modified, it will be in the process of constructing the PINNRepresentation.

Does this seem right to you?

[hpieper14]

I'm having a hard time understanding exactly how a PDESystem with independent array variables differs from a PDESystem initialized with other independent variables. For example, suppose I have the following two PDESystems

@parameters t z1 z2 z[1:2]
pde_system1 =  PDESystem(eq, bcs, domains, [t, z1, z2], [u(t, z1, z2)])
pde_system2 = PDESystem(eqs, bcs, domains, [t, z...], [u(t, z...)]) 

I can look at the independent variables and I would get

In: pde_system1.indvars
Out:  Num[t, z1, z2] 

In: pde_system2.indvars 
Out: Num[t, z[1], z[2]]

and both of these have type Vector{Num}. However, when I call ModelingToolkit.getname, I get

In: ModelingToolkit.getname.(pde_system1.indvars) 
Out: [:t, :z1, :z2] 

In: ModelingToolkit.getname.(pde_system2.indvars) 
Out: [:t, :z, :z] 

It seems to me that getname returns the name of the entire array instead of the name of the particular entry. Is there a way to unpack the independent variables of pde_system2 so that ModelingToolkit.getname returns a vector of symbols where each symbol corresponds to the entry of the array variables?

I thought an alternative approach would be to splat z but it errors

In: ModelingToolkit.getname.(pde_system2.indvars)
Out: ERROR: MethodError: no method matching getname(::Num, ::Num, ::Num)
Closest candidates are:
  getname(::Any, ::Any) at ~/.julia/packages/Symbolics/sDAUx/src/variable.jl:364
  getname(::Any) at ~/.julia/packages/Symbolics/sDAUx/src/variable.jl:364

[hpieper14]

@ChrisRackauckas let me know when you get a chance to look at this!

[ChrisRackauckas]

We want to extend the discretization framework so that it can work with PDESystems that have splatted vector variables, i.e. pde_system = PDESystem(eq, bcs, domains, [x,y...], [u(x,y...)]). It looks like to me that the changes for this need to occur in the files symbolic_utilities.jl and discretize.jl. Symbolic_utilities.jl is where the current error is being thrown and as suggested by name, has symbolic parsing helper functions and if discretize.jl needs to be modified, it will be in the process of constructing the PINNRepresentation.

Does this seem right to you?

That seems right. And I think to do that, the canonical form needs to switch to something that is shared between the two. Since it seems like there's issues with the names on the array variables, instead of checking by name (which is ill-advised anyways), it should probably just keep the symbolic value and do direct symbolic isequals checks. So basically, getname just shouldn't be used at all, and anything that's like getname(x) == getname(y) should just be isequals(x,y).