Skip to content

Commit

Permalink
fix: eltype conversions in IntegroDiff
Browse files Browse the repository at this point in the history
  • Loading branch information
avik-pal committed Oct 15, 2024
1 parent c297480 commit ff778e4
Show file tree
Hide file tree
Showing 9 changed files with 87 additions and 109 deletions.
8 changes: 4 additions & 4 deletions lib/NeuralPDELogging/test/adaptive_loss_log_tests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -71,12 +71,12 @@ function test_2d_poisson_equation_adaptive_loss(adaptive_loss, run, outdir, hasl
if iteration[1] % 30 == 0
u_predict = reshape([first(phi([x, y], p.u)) for x in xs for y in ys],
(length(xs), length(ys)))
log_value(logger, "outer_error/total_diff",
sum(abs2, u_predict .- u_real), step = iteration[1])
total_diff = sum(abs, u_predict .- u_real)
log_value(logger, "outer_error/total_diff", total_diff, step = iteration[1])
log_value(logger, "outer_error/total_diff_rel",
total_diff / sum(abs2, u_real), step = iteration[1])
log_value(logger, "outer_error/total_diff_sq", sum(diff_u .^ 2),
step = iteration[1])
log_value(logger, "outer_error/total_diff_sq",
sum(abs2, u_predict .- u_real), step = iteration[1])
end
end
return false
Expand Down
3 changes: 1 addition & 2 deletions src/BPINN_ode.jl
Original file line number Diff line number Diff line change
Expand Up @@ -117,8 +117,7 @@ function BNNODE(chain, Kernel = HMC; strategy = nothing, draw_samples = 2000,
numensemble = floor(Int, draw_samples / 3),
estim_collocate = false,
autodiff = false, progress = false, verbose = false)
!(chain isa AbstractLuxLayer) &&
(chain = adapt(FromFluxAdaptor(false, false), chain))
chain isa AbstractLuxLayer || (chain = FromFluxAdaptor()(chain))
BNNODE(chain, Kernel, strategy,
draw_samples, priorsNNw, param, l2std,
phystd, dataset, physdt, MCMCkwargs,
Expand Down
32 changes: 18 additions & 14 deletions src/discretize.jl
Original file line number Diff line number Diff line change
Expand Up @@ -189,17 +189,17 @@ function generate_training_sets(domains, dx, eqs, bcs, eltypeθ, dict_indvars::D
dict_var_span_ = Dict([Symbol(d.variables) => bc for (d, bc) in zip(domains, bc_data)])

bcs_train_sets = map(bound_args) do bt
span = map(b -> get(dict_var_span, b, b), bt)
return adapt(eltypeθ,
hcat(vec(map(points -> collect(points), Iterators.product(span...)))...))
span = get.((dict_var_span,), bt, bt)
return reduce(hcat, vec(map(collect, Iterators.product(span...)))) |>
EltypeAdaptor{eltypeθ}()
end

pde_args = get_argument(eqs, dict_indvars, dict_depvars)

pde_train_sets = map(pde_args) do bt
span = map(b -> get(dict_var_span_, b, b), bt)
return adapt(eltypeθ,
hcat(vec(map(points -> collect(points), Iterators.product(span...)))...))
span = get.((dict_var_span_,), bt, bt)
return reduce(hcat, vec(map(collect, Iterators.product(span...)))) |>
EltypeAdaptor{eltypeθ}()
end

return [pde_train_sets, bcs_train_sets]
Expand Down Expand Up @@ -231,13 +231,16 @@ function get_bounds(domains, eqs, bcs, eltypeθ, dict_indvars, dict_depvars,

pde_args = get_argument(eqs, dict_indvars, dict_depvars)

ϵ = cbrt(eps(eltypeθ))
eltype_adaptor = EltypeAdaptor{eltypeθ}()

pde_lower_bounds = map(pde_args) do pd
span = map(p -> get(dict_lower_bound, p, p), pd)
map(s -> adapt(eltypeθ, s) + cbrt(eps(eltypeθ)), span)
span = get.((dict_lower_bound,), pd, pd) |> eltype_adaptor
return span .+ ϵ
end
pde_upper_bounds = map(pde_args) do pd
span = map(p -> get(dict_upper_bound, p, p), pd)
map(s -> adapt(eltypeθ, s) - cbrt(eps(eltypeθ)), span)
span = get.((dict_upper_bound,), pd, pd) |> eltype_adaptor
return span .+ ϵ
end
pde_bounds = [pde_lower_bounds, pde_upper_bounds]

Expand Down Expand Up @@ -283,7 +286,7 @@ function get_numeric_integral(pinnrep::PINNRepresentation)
function integration_(cord, lb, ub, θ)
cord_ = cord
function integrand_(x, p)
@ignore_derivatives @views(cord_[integrating_var_id]) .= x
@ignore_derivatives cord_[integrating_var_id] .= x
return integrand_func(cord_, p, phi, derivative, nothing, u, nothing)
end
prob_ = IntegralProblem(integrand_, (lb, ub), θ)
Expand All @@ -296,15 +299,15 @@ function get_numeric_integral(pinnrep::PINNRepresentation)
ub_ = zeros(size(ub)[1], size(cord)[2])
for (i, l) in enumerate(lb)
if l isa Number
@ignore_derivatives lb_[i, :] = fill(l, 1, size(cord)[2])
@ignore_derivatives lb_[i, :] .= l
else
@ignore_derivatives lb_[i, :] = l(
cord, θ, phi, derivative, nothing, u, nothing)
end
end
for (i, u_) in enumerate(ub)
if u_ isa Number
@ignore_derivatives ub_[i, :] = fill(u_, 1, size(cord)[2])
@ignore_derivatives ub_[i, :] .= u_
else
@ignore_derivatives ub_[i, :] = u_(cord, θ, phi, derivative,
nothing, u, nothing)
Expand Down Expand Up @@ -333,8 +336,9 @@ upon to give Solution or the Solution Distribution of the PDE.
For more information, see `discretize` and `PINNRepresentation`.
"""
function SciMLBase.symbolic_discretize(pde_system::PDESystem, discretization::AbstractPINN)
(; eqs, bcs, domain, defaults) = pde_system
(; eqs, bcs, domain) = pde_system
eq_params = pde_system.ps
defaults = pde_system.defaults
(; chain, param_estim, additional_loss, multioutput, init_params, phi, derivative, strategy, logger, iteration, self_increment) = discretization
(; log_frequency) = discretization.log_options
adaloss = discretization.adaptive_loss
Expand Down
4 changes: 2 additions & 2 deletions src/eltype_matching.jl
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
struct EltypeAdaptor{T} end

(l::EltypeAdaptor)(x) = fmap(adapt(l), x)
(l::EltypeAdaptor)(x) = fmap(Adapt.adapt(l), x)
function (l::EltypeAdaptor)(x::AbstractArray{T}) where {T}
return isbitstype(T) ? adapt(l, x) : map(l, x)
return (isbitstype(T) || T <: Number) ? Adapt.adapt(l, x) : map(l, x)
end

function Adapt.adapt_storage(::EltypeAdaptor{T}, x::AbstractArray) where {T}
Expand Down
2 changes: 1 addition & 1 deletion src/pinn_types.jl
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,7 @@ function PhysicsInformedNN(
if multioutput
chain = map(chain) do cᵢ
cᵢ isa AbstractLuxLayer && return cᵢ
return adapt(FromFluxAdaptor(), cᵢ)
return FromFluxAdaptor()(cᵢ)
end
else
chain isa AbstractLuxLayer || (chain = FromFluxAdaptor()(chain))
Expand Down
5 changes: 2 additions & 3 deletions src/symbolic_utilities.jl
Original file line number Diff line number Diff line change
Expand Up @@ -424,9 +424,8 @@ function get_argument end

# Get arguments from boundary condition functions
function get_argument(eqs, _indvars::Array, _depvars::Array)
depvars, indvars, dict_indvars, dict_depvars, dict_depvar_input = get_vars(_indvars,
_depvars)
get_argument(eqs, dict_indvars, dict_depvars)
_, _, dict_indvars, dict_depvars, _ = get_vars(_indvars, _depvars)
return get_argument(eqs, dict_indvars, dict_depvars)
end
function get_argument(eqs, dict_indvars, dict_depvars)
exprs = toexpr.(eqs)
Expand Down
35 changes: 16 additions & 19 deletions test/IDE_tests.jl
Original file line number Diff line number Diff line change
@@ -1,10 +1,7 @@
using Test, NeuralPDE
using Optimization, OptimizationOptimJL
using Test, NeuralPDE, Optimization, OptimizationOptimJL, DomainSets, Lux, Random,
Statistics
import ModelingToolkit: Interval
using DomainSets, Flux
import Lux

using Random
Random.seed!(110)

callback = function (p, l)
Expand All @@ -20,7 +17,7 @@ end
eq = Di(i(t)) + 2 * i(t) + 5 * Ii(i(t)) ~ 1
bcs = [i(0.0) ~ 0.0]
domains = [t Interval(0.0, 2.0)]
chain = Lux.Chain(Lux.Dense(1, 15, Lux.σ), Lux.Dense(15, 1))
chain = Chain(Dense(1, 15, σ), Dense(15, 1))
strategy_ = GridTraining(0.1)
discretization = PhysicsInformedNN(chain, strategy_)
@named pde_system = PDESystem(eq, bcs, domains, [t], [i(t)])
Expand All @@ -31,7 +28,7 @@ end
analytic_sol_func(t) = 1 / 2 * (exp(-t)) * (sin(2 * t))
u_real = [analytic_sol_func(t) for t in ts]
u_predict = [first(phi([t], res.u)) for t in ts]
@test Flux.mse(u_real, u_predict) < 0.01
@test mean(abs2, u_real .- u_predict) < 0.01
end

@testset "Example 2 - 1D" begin
Expand All @@ -45,7 +42,7 @@ end

bcs = [u(0.0) ~ 0.0]
domains = [x Interval(0.0, 1.00)]
chain = Lux.Chain(Lux.Dense(1, 15, Lux.σ), Lux.Dense(15, 1))
chain = Chain(Dense(1, 15, σ), Dense(15, 1))
strategy_ = GridTraining(0.1)
discretization = PhysicsInformedNN(chain, strategy_)
@named pde_system = PDESystem(eq, bcs, domains, [x], [u(x)])
Expand All @@ -56,7 +53,7 @@ end
phi = discretization.phi
u_predict = [first(phi([x], res.u)) for x in xs]
u_real = [x^2 / cos(x) for x in xs]
@test Flux.mse(u_real, u_predict) < 0.001
@test mean(abs2, u_real .- u_predict) < 0.01
end

@testset "Example 3 - 2 Inputs, 1 Output" begin
Expand All @@ -68,7 +65,7 @@ end
eq = Ix(u(x, y)) ~ 1 / 3
bcs = [u(0.0, 0.0) ~ 1, Dx(u(x, y)) ~ -2 * x, Dy(u(x, y)) ~ -2 * y]
domains = [x Interval(0.0, 1.00), y Interval(0.0, 1.00)]
chain = Lux.Chain(Lux.Dense(2, 15, Lux.σ), Lux.Dense(15, 1))
chain = Chain(Dense(2, 15, σ), Dense(15, 1))
strategy_ = GridTraining(0.1)
discretization = PhysicsInformedNN(chain, strategy_)
@named pde_system = PDESystem(eq, bcs, domains, [x, y], [u(x, y)])
Expand All @@ -79,7 +76,7 @@ end
phi = discretization.phi
u_real = collect(1 - x^2 - y^2 for y in ys, x in xs)
u_predict = collect(Array(phi([x, y], res.u))[1] for y in ys, x in xs)
@test Flux.mse(u_real, u_predict) < 0.001
@test mean(abs2, u_real .- u_predict) < 0.001
end

@testset "Example 4 - 2 Inputs, 1 Output" begin
Expand All @@ -91,7 +88,7 @@ end
eq = Ix(u(x, y)) ~ 5 / 12
bcs = [u(0.0, 0.0) ~ 0, Dy(u(x, y)) ~ 2 * y, u(x, 0) ~ x]
domains = [x Interval(0.0, 1.00), y Interval(0.0, 1.00)]
chain = Lux.Chain(Lux.Dense(2, 15, Lux.σ), Lux.Dense(15, 1))
chain = Chain(Dense(2, 15, σ), Dense(15, 1))
strategy_ = GridTraining(0.1)
discretization = PhysicsInformedNN(chain, strategy_)
@named pde_system = PDESystem(eq, bcs, domains, [x, y], [u(x, y)])
Expand All @@ -102,7 +99,7 @@ end
phi = discretization.phi
u_real = collect(x + y^2 for y in ys, x in xs)
u_predict = collect(Array(phi([x, y], res.u))[1] for y in ys, x in xs)
@test Flux.mse(u_real, u_predict) < 0.01
@test mean(abs2, u_real .- u_predict) < 0.01
end

@testset "Example 5 - 1 Input, 2 Outputs" begin
Expand All @@ -113,7 +110,7 @@ end
eqs = [Ix(u(x) * w(x)) ~ log(abs(x)), Dx(w(x)) ~ -2 / (x^3), u(x) ~ x]
bcs = [u(1.0) ~ 1.0, w(1.0) ~ 1.0]
domains = [x Interval(1.0, 2.0)]
chains = [Lux.Chain(Lux.Dense(1, 15, Lux.σ), Lux.Dense(15, 1)) for _ in 1:2]
chains = [Chain(Dense(1, 15, σ), Dense(15, 1)) for _ in 1:2]
strategy_ = GridTraining(0.1)
discretization = PhysicsInformedNN(chains, strategy_)
@named pde_system = PDESystem(eqs, bcs, domains, [x], [u(x), w(x)])
Expand All @@ -125,8 +122,8 @@ end
w_predict = [(phi[2]([x], res.u.depvar.w))[1] for x in xs]
u_real = [x for x in xs]
w_real = [1 / x^2 for x in xs]
@test Flux.mse(u_real, u_predict) < 0.001
@test Flux.mse(w_real, w_predict) < 0.001
@test mean(abs2, u_real .- u_predict) < 0.001
@test mean(abs2, w_real .- w_predict) < 0.001
end

@testset "Example 6: Infinity" begin
Expand All @@ -137,7 +134,7 @@ end
eqs = [I(u(x)) ~ Iinf(u(x)) - 1 / x]
bcs = [u(1) ~ 1]
domains = [x Interval(1.0, 2.0)]
chain = Lux.Chain(Lux.Dense(1, 10, Lux.σ), Lux.Dense(10, 1))
chain = Chain(Dense(1, 10, σ), Dense(10, 1))
discretization = PhysicsInformedNN(chain, NeuralPDE.GridTraining(0.1))
@named pde_system = PDESystem(eqs, bcs, domains, [x], [u(x)])
prob = discretize(pde_system, discretization)
Expand All @@ -146,7 +143,7 @@ end
phi = discretization.phi
u_predict = [first(phi([x], res.u)) for x in xs]
u_real = [1 / x^2 for x in xs]
@test u_realu_predict rtol=10^-2
@test u_realu_predict rtol=10^-1
end

@testset "Example 7: Infinity" begin
Expand All @@ -156,7 +153,7 @@ end
eq = I(u(x)) ~ 1 / x
domains = [x Interval(1.0, 2.0)]
bcs = [u(1) ~ 1]
chain = Lux.Chain(Lux.Dense(1, 12, Lux.tanh), Lux.Dense(12, 1))
chain = Chain(Dense(1, 12, tanh), Dense(12, 1))
discretization = PhysicsInformedNN(chain, GridTraining(0.1))
@named pde_system = PDESystem(eq, bcs, domains, [x], [u(x)])
prob = discretize(pde_system, discretization)
Expand Down
Loading

0 comments on commit ff778e4

Please sign in to comment.