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

PINNErrorVsTime Benchmark Updates #1159

Open
wants to merge 8 commits into
base: master
Choose a base branch
from

Conversation

ParamThakkar123
Copy link
Contributor

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

Additional context

Add any other context about the problem here.

@ParamThakkar123
Copy link
Contributor Author

@ChrisRackauckas I got an error on running the iterations which said that the maxiters are less than 1000 so I set all maxiters to 1100. Actually the decision was a bit arbitrary but is that a good number ??

@ChrisRackauckas
Copy link
Member

https://docs.sciml.ai/SciMLBenchmarksOutput/v0.5/PINNErrorsVsTime/diffusion_et/

It's supposed to just show the error over time and then get cutoff. I don't see why making it longer would help.

@ParamThakkar123
Copy link
Contributor Author

ParamThakkar123 commented Jan 19, 2025

https://docs.sciml.ai/SciMLBenchmarksOutput/v0.5/PINNErrorsVsTime/diffusion_et/

It's supposed to just show the error over time and then get cutoff. I don't see why making it longer would help.

Yes. Actually I set it to that number just to get rid of that error

@ChrisRackauckas
Copy link
Member

Wait what's the error?

@ParamThakkar123
Copy link
Contributor Author

Wait what's the error?

the error went off on running again

@ChrisRackauckas
Copy link
Member

what error?

@ParamThakkar123
Copy link
Contributor Author

maxiters should be a number greater than 1000

@ChrisRackauckas
Copy link
Member

can you please just show the error...

@ParamThakkar123
Copy link
Contributor Author

AssertionError: maxiters for CubaCuhre(0, 0, 0) should be larger than 1000

Stacktrace: [1] __solvebp_call(prob::IntegralProblem{false, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(layer_1 = ViewAxis(1:60, Axis(weight = ViewAxis(1:50, ShapedAxis((10, 5))), bias = ViewAxis(51:60, ShapedAxis((10, 1))))), layer_2 = ViewAxis(61:170, Axis(weight = ViewAxis(1:100, ShapedAxis((10, 10))), bias = ViewAxis(101:110, ShapedAxis((10, 1))))), layer_3 = ViewAxis(171:181, Axis(weight = ViewAxis(1:10, ShapedAxis((1, 10))), bias = ViewAxis(11:11, ShapedAxis((1, 1))))))}}}, NeuralPDE.var"#integrand#109"{NeuralPDE.var"#219#220"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:cord, Symbol("##θ#228"), :phi, :derivative, :integral, :u, :p),
NeuralPDE.var"#_RGF_ModTag", NeuralPDE.var"#_RGF_ModTag", (0xd9696f1d, 0xe356e73c, 0x32906e9c, 0x54a064bc, 0x0cbbe458), Expr}, NeuralPDE.var"#12#13", NeuralPDE.var"#279#286"{NeuralPDE.var"#279#280#287"{typeof(NeuralPDE.numeric_derivative)}, Dict{Symbol, Int64}, Dict{Symbol, Int64}, QuadratureTraining{CubaCuhre, Float64}}, typeof(NeuralPDE.numeric_derivative), NeuralPDE.Phi{Chain{@NamedTuple{layer_1::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_2::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_3::Dense{true, typeof(identity), typeof(glorot_uniform), typeof(zeros32)}}, Nothing}, @NamedTuple{layer_1::@NamedTuple{}, layer_2::@NamedTuple{}, layer_3::@NamedTuple{}}}, Nothing}}, Vector{Float64}, @Kwargs{}}, alg::CubaCuhre,
sensealg::Integrals.ReCallVJP{Integrals.ZygoteVJP}, lb::Vector{Float64}, ub::Vector{Float64}, p::ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(layer_1 = ViewAxis(1:60, Axis(weight = ViewAxis(1:50, ShapedAxis((10, 5))), bias = ViewAxis(51:60, ShapedAxis((10, 1))))), layer_2 = ViewAxis(61:170, Axis(weight = ViewAxis(1:100, ShapedAxis((10, 10))), bias = ViewAxis(101:110, ShapedAxis((10, 1))))), layer_3 = ViewAxis(171:181, Axis(weight = ViewAxis(1:10, ShapedAxis((1, 10))), bias = ViewAxis(11:11, ShapedAxis((1, 1))))))}}}; reltol::Float64, abstol::Float64, maxiters::Int64) @ IntegralsCuba C:\Users\Hp.julia\packages\IntegralsCuba\xueKH\src\IntegralsCuba.jl:139 [2] __solvebp_call @ C:\Users\Hp.julia\packages\IntegralsCuba\xueKH\src\IntegralsCuba.jl:134 [inlined] [3] #__solvebp_call#4 @ C:\Users\Hp.julia\packages\Integrals\d3rQd\src\common.jl:95 [inlined] [4] __solvebp_call @ C:\Users\Hp.julia\packages\Integrals\d3rQd\src\common.jl:94 [inlined] [5] #rrule#5 @ C:\Users\Hp.julia\packages\Integrals\d3rQd\ext\IntegralsZygoteExt.jl:17 [inlined] [6] rrule @ C:\Users\Hp.julia\packages\Integrals\d3rQd\ext\IntegralsZygoteExt.jl:14 [inlined] [7] rrule @ C:\Users\Hp.julia\packages\ChainRulesCore\U6wNx\src\rules.jl:144 [inlined] [8] chain_rrule_kw @ C:\Users\Hp.julia\packages\Zygote\TWpme\src\compiler\chainrules.jl:236 [inlined] [9] macro expansion @ C:\Users\Hp.julia\packages\Zygote\TWpme\src\compiler\interface2.jl:0 [inlined] [10] _pullback @ C:\Users\Hp.julia\packages\Zygote\TWpme\src\compiler\interface2.jl:91 [inlined] [11] solve! @ C:\Users\Hp.julia\packages\Integrals\d3rQd\src\common.jl:84 [inlined] ... @ SciMLBase C:\Users\Hp.julia\packages\SciMLBase\szsYq\src\solve.jl:162 [55] solve(::OptimizationProblem{true, OptimizationFunction{true, AutoZygote, NeuralPDE.var"#full_loss_function#318"{NeuralPDE.var"#null_nonadaptive_loss#118", Vector{NeuralPDE.var"#106#110"{NeuralPDE.var"#219#220"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:cord, Symbol("##θ#228"), :phi, :derivative, :integral, :u, :p), NeuralPDE.var"#_RGF_ModTag", NeuralPDE.var"#_RGF_ModTag", (0x0ca831bf, 0x51abc8eb, 0xda1f388f, 0xf472bcea, 0x7492cfcb), Expr}, NeuralPDE.var"#12#13", NeuralPDE.var"#279#286"{NeuralPDE.var"#279#280#287"{typeof(NeuralPDE.numeric_derivative)}, Dict{Symbol, Int64}, Dict{Symbol, Int64}, QuadratureTraining{CubaCuhre, Float64}}, typeof(NeuralPDE.numeric_derivative), NeuralPDE.Phi{Chain{@NamedTuple{layer_1::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_2::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_3::Dense{true, typeof(identity), typeof(glorot_uniform), typeof(zeros32)}}, Nothing}, @NamedTuple{layer_1::@NamedTuple{}, layer_2::@NamedTuple{}, layer_3::@NamedTuple{}}}, Nothing}, Vector{Float64}, Vector{Float64}, NeuralPDE.var"#105#108"{QuadratureTraining{CubaCuhre, Float64}}, Float64}}, Vector{NeuralPDE.var"#106#110"{NeuralPDE.var"#219#220"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:cord, Symbol("##θ#228"), :phi, :derivative, :integral, :u, :p), NeuralPDE.var"#_RGF_ModTag", NeuralPDE.var"#RGF_ModTag", (0xd9696f1d, 0xe356e73c, 0x32906e9c, 0x54a064bc, 0x0cbbe458), Expr}, NeuralPDE.var"#12#13", NeuralPDE.var"#279#286"{NeuralPDE.var"#279#280#287"{typeof(NeuralPDE.numeric_derivative)}, Dict{Symbol, Int64}, Dict{Symbol, Int64}, QuadratureTraining{CubaCuhre, Float64}}, typeof(NeuralPDE.numeric_derivative), NeuralPDE.Phi{Chain{@NamedTuple{layer_1::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_2::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_3::Dense{true, typeof(identity), typeof(glorot_uniform), typeof(zeros32)}}, Nothing}, @NamedTuple{layer_1::@NamedTuple{}, layer_2::@NamedTuple{}, layer_3::@NamedTuple{}}}, Nothing}, Vector{Float64}, Vector{Float64}, NeuralPDE.var"#105#108"{QuadratureTraining{CubaCuhre, Float64}}, Float64}}, NeuralPDE.PINNRepresentation, Bool, Vector{Int64}, Int64, NeuralPDE.Phi{Chain{@NamedTuple{layer_1::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_2::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_3::Dense{true, typeof(identity), typeof(glorot_uniform), typeof(zeros32)}}, Nothing}, @NamedTuple{layer_1::@NamedTuple{}, layer_2::@NamedTuple{}, layer_3::@NamedTuple{}}}, Nothing, Bool, Nothing}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(layer_1 = ViewAxis(1:60, Axis(weight = ViewAxis(1:50, ShapedAxis((10, 5))), bias = ViewAxis(51:60, ShapedAxis((10, 1))))), layer_2 = ViewAxis(61:170, Axis(weight = ViewAxis(1:100, ShapedAxis((10, 10))), bias = ViewAxis(101:110, ShapedAxis((10, 1))))), layer_3 = ViewAxis(171:181, Axis(weight = ViewAxis(1:10, ShapedAxis((1, 10))), bias = ViewAxis(11:11, ShapedAxis((1, 1))))))}}}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, @Kwargs{}}, ::Optimisers.Adam; kwargs::@Kwargs{callback::var"#11#18"{var"#loss_function#17"{OptimizationProblem{true, OptimizationFunction{true, AutoZygote, NeuralPDE.var"#full_loss_function#318"{NeuralPDE.var"#null_nonadaptive_loss#118", Vector{NeuralPDE.var"#74#75"{NeuralPDE.var"#219#220"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:cord, Symbol("##θ#228"), :phi, :derivative, :integral, :u, :p), NeuralPDE.var"#_RGF_ModTag", NeuralPDE.var"#_RGF_ModTag", (0x700712a1, 0x6c1c91c2, 0xa5bfc01b, 0x66a91103, 0x0d12fcff), Expr}, NeuralPDE.var"#12#13", NeuralPDE.var"#279#286"{NeuralPDE.var"#279#280#287"{typeof(NeuralPDE.numeric_derivative)}, Dict{Symbol, Int64}, Dict{Symbol, Int64}, GridTraining{Float64}}, typeof(NeuralPDE.numeric_derivative), NeuralPDE.Phi{Chain{@NamedTuple{layer_1::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_2::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_3::Dense{true, typeof(identity), typeof(glorot_uniform), typeof(zeros32)}}, Nothing}, @NamedTuple{layer_1::@NamedTuple{}, layer_2::@NamedTuple{}, layer_3::@NamedTuple{}}}, Nothing}, Matrix{Real}}}, Vector{NeuralPDE.var"#74#75"{NeuralPDE.var"#219#220"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:cord, Symbol("##θ#228"), :phi, :derivative, :integral, :u, :p), NeuralPDE.var"#_RGF_ModTag", NeuralPDE.var"#_RGF_ModTag", (0xd9696f1d, 0xe356e73c, 0x32906e9c, 0x54a064bc, 0x0cbbe458), Expr}, NeuralPDE.var"#12#13", NeuralPDE.var"#279#286"{NeuralPDE.var"#279#280#287"{typeof(NeuralPDE.numeric_derivative)}, Dict{Symbol, Int64}, Dict{Symbol, Int64}, GridTraining{Float64}}, typeof(NeuralPDE.numeric_derivative), NeuralPDE.Phi{Chain{@NamedTuple{layer_1::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_2::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_3::Dense{true, typeof(identity), typeof(glorot_uniform), typeof(zeros32)}}, Nothing}, @NamedTuple{layer_1::@NamedTuple{}, layer_2::@NamedTuple{}, layer_3::@NamedTuple{}}}, Nothing}, Matrix{Real}}}, NeuralPDE.PINNRepresentation, Bool, Vector{Int64}, Int64, NeuralPDE.Phi{Chain{@NamedTuple{layer_1::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_2::Dense{true, typeof(tanh_fast), typeof(glorot_uniform), typeof(zeros32)}, layer_3::Dense{true, typeof(identity), typeof(glorot_uniform), typeof(zeros32)}}, Nothing}, @NamedTuple{layer_1::@NamedTuple{}, layer_2::@NamedTuple{}, layer_3::@NamedTuple{}}}, Nothing, Bool, Nothing}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(layer_1 = ViewAxis(1:60, Axis(weight = ViewAxis(1:50, ShapedAxis((10, 5))), bias = ViewAxis(51:60, ShapedAxis((10, 1))))), layer_2 = ViewAxis(61:170, Axis(weight = ViewAxis(1:100, ShapedAxis((10, 10))), bias = ViewAxis(101:110, ShapedAxis((10, 1))))), layer_3 = ViewAxis(171:181, Axis(weight = ViewAxis(1:10, ShapedAxis((1, 10))), bias = ViewAxis(11:11, ShapedAxis((1, 1))))))}}}, SciMLBase.NullParameters, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, @Kwargs{}}}, Vector{Any}, Vector{Any}, Vector{Any}}, maxiters::Int64}) @ SciMLBase C:\Users\Hp.julia\packages\SciMLBase\szsYq\src\solve.jl:83 [56] allen_cahn(strategy::QuadratureTraining{CubaCuhre, Float64}, minimizer::Optimisers.Adam, maxIters::Int64) @ Main e:\SciMLBenchmarks.jl\benchmarks\PINNErrorsVsTime\jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_W1sZmlsZQ==.jl:105

@ChrisRackauckas
Copy link
Member

I see, that's for the sampling algorithm. You should only need that on Cuhre?

@ParamThakkar123
Copy link
Contributor Author

I see, that's for the sampling algorithm. You should only need that on Cuhre?

Yes. But as Cuhre was the first one in the line I thought setting to 1100 just for it would not solve the problem, so I set it to 1100 for all of them

@ParamThakkar123
Copy link
Contributor Author

The CI has passed here. And all the code seems to run perfectly. Can you please review ??

@ChrisRackauckas
Copy link
Member

@ArnoStrouwen SciML/Integrals.jl#124 can you remind me what the purpose behind this was?

@ArnoStrouwen
Copy link
Member

I don't remember myself, but that PR links to:
SciML/Integrals.jl#47

@ChrisRackauckas
Copy link
Member

Uninitialized memory in the original C: giordano/Cuba.jl#12 (comment) fantastic stuff numerical community, that's your classic method that everyone says when they say "all of the old stuff is robust" 😅

@ChrisRackauckas
Copy link
Member

Can you force latest majors and make sure the manifest resolves?

@ParamThakkar123
Copy link
Contributor Author

I bump forced the latest versions and resolved manifests but initially there were a lot of version conflicts. I removed IntegralsCuba and IntegralsCubature for a while to resolve them. The manifest resolved but adding both of them back again poses some more version conflicts

@ChrisRackauckas
Copy link
Member

Can you share the resolution errors?

@ParamThakkar123
Copy link
Contributor Author

image

image

@ChrisRackauckas These are the resolution errors that occur

@ChrisRackauckas
Copy link
Member

Oh those were turned into extensions. Change using IntegralsCuba, IntegralsCubature into using Integrals, Cuba, Cubature and change the dependencies to directly depending on Cuba and Cubature.

@ParamThakkar123
Copy link
Contributor Author

Sure !! 🫡

@ParamThakkar123
Copy link
Contributor Author

Error: MethodError: no method matching (::MLDataDevices.UnknownDevice)(::Ma
trix{Float64})

@ChrisRackauckas
Copy link
Member

Should I set it to 1000+ ??

yes

@ParamThakkar123
Copy link
Contributor Author

ParamThakkar123 commented Jan 22, 2025

@ChrisRackauckas I got a few questions to ask for the diffusion_et.jmd and hamilton_jacobi.jmd

in diffusion_et.jmd on running the code, I get this error
image

Where this error comes from, it didn't occur before changing it to Cuba and Cubature

@ParamThakkar123
Copy link
Contributor Author

In hamilton_jacobi.jmd, I got the following error and didn't get where it comes from

image

@ChrisRackauckas
Copy link
Member

@lxvm the first one, that looks like a case that might've regressed with AD?

@avik-pal is safe_get_device safe for using on arrays / CPU?

@avik-pal
Copy link
Member

What is safr_get_device? There is get_device from MLDataDevices which is implemented for all backend including CPU

@ChrisRackauckas
Copy link
Member

@ChrisRackauckas
Copy link
Member

Is UnknownDevice new? I don't think I've seen that.

@avik-pal
Copy link
Member

Oh that one is a stopgap solution for some of the Any issues that were coming up inside NeuralPDE

Is UnknownDevice new? I don't think I've seen that.

generally this doesn't get exposed to the user. You can get it if you call get_device(fn) and this fn is not a closure, for eg. tanh (in this case we don't know which device it can run on).

@ParamThakkar123
Copy link
Contributor Author

So, What changes do we need here ?? or is it an issue needed to be addressed in some other repository ??

@ParamThakkar123
Copy link
Contributor Author

So, What changes do we need here ?? or is it an issue needed to be addressed in some other repository ??

@ChrisRackauckas ??

@ParamThakkar123
Copy link
Contributor Author

I think these are the only 2 errors remaining here

@ParamThakkar123
Copy link
Contributor Author

@ChrisRackauckas any suggested fix for the MLDevices error here ?? It's common in all files and probably coming from NeuralPDE.jl

@ChrisRackauckas
Copy link
Member

I don't know the solution, I'm not sure why MLDevices is involved at all since it's just on the CPU.

@ParamThakkar123
Copy link
Contributor Author

I don't know the solution, I'm not sure why MLDevices is involved at all since it's just on the CPU.

Yeah. I can understand. But just wanted to know if there's some other alternative to this, or what can be done if we didn't figure what's happening

@ChrisRackauckas
Copy link
Member

Make an MWE without the highest level package

@ParamThakkar123
Copy link
Contributor Author

@ChrisRackauckas I implemented an MWE using all the dependencies in codes here except the highest level package. Here's the code:

using Integrals, Cubature, Cuba
using ModelingToolkit, Optimization, OptimizationOptimJL
using Lux, Plots
using DelimitedFiles
using OptimizationOptimisers
using QuasiMonteCarlo
import ModelingToolkit: Interval, infimum, supremum

# Create a sample problem
rng = Xoshiro(0)
dim = 2
model = Lux.Chain(
    Lux.Dense(dim, 16, relu),
    Lux.Dense(16, 8, relu),
    Lux.Dense(8, 1)
)
ps, st = Lux.setup(rng, model)

# Objective function to minimize
function objective(x, p)
    sum(x.^2)
end

# Bounds for optimization
lower_bounds = [-2.0, -2.0]
upper_bounds = [2.0, 2.0]

# Quasi-Monte Carlo sampling
num_samples = 100
sample_points = QuasiMonteCarlo.sample(num_samples, lower_bounds, upper_bounds, LatticeRuleSample())

# Set up optimization problem
prob = OptimizationProblem(
    OptimizationFunction(objective, Optimization.AutoForwardDiff()),
    zeros(dim),
    ps;
    lb = lower_bounds,
    ub = upper_bounds
)

# Solve using different optimizers
sol1 = solve(prob, NelderMead())
sol2 = solve(prob, BFGS())

# Demonstration of integration
function test_integral(x)
    sin(x[1]) * cos(x[2])
end

# Use different integration methods
int2, = hcubature(test_integral, lower_bounds, upper_bounds)

# Plotting results
p1 = scatter(
    sample_points[1,:], 
    sample_points[2,:], 
    title="Quasi-Monte Carlo Samples", 
    xlabel="X1", 
    ylabel="X2"
)

p2 = plot(
    [sol1.minimizer[1]], 
    [sol1.minimizer[2]], 
    seriestype = :scatter, 
    title="Optimization Results", 
    label="Nelder-Mead"
)
plot!(p2, [sol2.minimizer[1]], [sol2.minimizer[2]], 
    seriestype = :scatter, 
    label="BFGS")

# Save results
writedlm("optimization_results.csv", 
    [sol1.minimizer sol2.minimizer], 
    ',')

println("Nelder-Mead Solution: ", sol1.minimizer)
println("BFGS Solution: ", sol2.minimizer)
println("Integral (hcubature): ", int2)

This code does:

  • optimization of a simple quadratic function using Nelder-Mead and BFGS methods
  • Quasi-Monte Carlo sampling of 100 points in a 2D space
  • Numerical integration of a test function using hcubature

@ParamThakkar123
Copy link
Contributor Author

@ChrisRackauckas I implemented an MWE using all the dependencies in codes here except the highest level package. Here's the code:

using Integrals, Cubature, Cuba
using ModelingToolkit, Optimization, OptimizationOptimJL
using Lux, Plots
using DelimitedFiles
using OptimizationOptimisers
using QuasiMonteCarlo
import ModelingToolkit: Interval, infimum, supremum

# Create a sample problem
rng = Xoshiro(0)
dim = 2
model = Lux.Chain(
    Lux.Dense(dim, 16, relu),
    Lux.Dense(16, 8, relu),
    Lux.Dense(8, 1)
)
ps, st = Lux.setup(rng, model)

# Objective function to minimize
function objective(x, p)
    sum(x.^2)
end

# Bounds for optimization
lower_bounds = [-2.0, -2.0]
upper_bounds = [2.0, 2.0]

# Quasi-Monte Carlo sampling
num_samples = 100
sample_points = QuasiMonteCarlo.sample(num_samples, lower_bounds, upper_bounds, LatticeRuleSample())

# Set up optimization problem
prob = OptimizationProblem(
    OptimizationFunction(objective, Optimization.AutoForwardDiff()),
    zeros(dim),
    ps;
    lb = lower_bounds,
    ub = upper_bounds
)

# Solve using different optimizers
sol1 = solve(prob, NelderMead())
sol2 = solve(prob, BFGS())

# Demonstration of integration
function test_integral(x)
    sin(x[1]) * cos(x[2])
end

# Use different integration methods
int2, = hcubature(test_integral, lower_bounds, upper_bounds)

# Plotting results
p1 = scatter(
    sample_points[1,:], 
    sample_points[2,:], 
    title="Quasi-Monte Carlo Samples", 
    xlabel="X1", 
    ylabel="X2"
)

p2 = plot(
    [sol1.minimizer[1]], 
    [sol1.minimizer[2]], 
    seriestype = :scatter, 
    title="Optimization Results", 
    label="Nelder-Mead"
)
plot!(p2, [sol2.minimizer[1]], [sol2.minimizer[2]], 
    seriestype = :scatter, 
    label="BFGS")

# Save results
writedlm("optimization_results.csv", 
    [sol1.minimizer sol2.minimizer], 
    ',')

println("Nelder-Mead Solution: ", sol1.minimizer)
println("BFGS Solution: ", sol2.minimizer)
println("Integral (hcubature): ", int2)

This code does:

  • optimization of a simple quadratic function using Nelder-Mead and BFGS methods
  • Quasi-Monte Carlo sampling of 100 points in a 2D space
  • Numerical integration of a test function using hcubature

@ChrisRackauckas ??

@ChrisRackauckas
Copy link
Member

Well that still has neuralpde. Next is to make the MWE without NeuralPDE

@ParamThakkar123
Copy link
Contributor Author

Well that still has neuralpde. Next is to make the MWE without NeuralPDE

Should I rewrite all the benchmark codes without using NeuralPDE, you mean by this ??

@ChrisRackauckas
Copy link
Member

The MWE. Find out what the bug is.

@ParamThakkar123
Copy link
Contributor Author

using Integrals, Cubature, Cuba
using ModelingToolkit, Optimization, OptimizationOptimJL
using Lux, Plots
using DelimitedFiles
using Random
using OptimizationOptimisers
using QuasiMonteCarlo
import ModelingToolkit: Interval, infimum, supremum

# Create a sample problem
rng = Xoshiro(0)
dim = 2
model = Lux.Chain(
    Lux.Dense(dim, 16, relu),
    Lux.Dense(16, 8, relu),
    Lux.Dense(8, 1)
)
ps, st = Lux.setup(rng, model)

# Objective function to minimize
function objective(x, p)
    sum(x.^2)
end

# Bounds for optimization
lower_bounds = [-2.0, -2.0]
upper_bounds = [2.0, 2.0]

# Quasi-Monte Carlo sampling
num_samples = 100
sample_points = QuasiMonteCarlo.sample(num_samples, lower_bounds, upper_bounds, LatticeRuleSample())

# Set up optimization problem
prob = OptimizationProblem(
    OptimizationFunction(objective, Optimization.AutoForwardDiff()),
    zeros(dim),
    ps;
    lb = lower_bounds,
    ub = upper_bounds
)

# Solve using different optimizers
sol1 = solve(prob, NelderMead())
sol2 = solve(prob, BFGS())

# Demonstration of integration
function test_integral(x)
    sin(x[1]) * cos(x[2])
end

# Use different integration methods
int2, = hcubature(test_integral, lower_bounds, upper_bounds)

# Plotting results
p1 = scatter(
    sample_points[1,:], 
    sample_points[2,:], 
    title="Quasi-Monte Carlo Samples", 
    xlabel="X1", 
    ylabel="X2"
)

p2 = plot(
    [sol1.minimizer[1]], 
    [sol1.minimizer[2]], 
    seriestype = :scatter, 
    title="Optimization Results", 
    label="Nelder-Mead"
)
plot!(p2, [sol2.minimizer[1]], [sol2.minimizer[2]], 
    seriestype = :scatter, 
    label="BFGS")

# Save results
writedlm("optimization_results.csv", 
    [sol1.minimizer sol2.minimizer], 
    ',')

println("Nelder-Mead Solution: ", sol1.minimizer)
println("BFGS Solution: ", sol2.minimizer)
println("Integral (hcubature): ", int2)

@ChrisRackauckas The code I previously sent had a bug where Xoshiro wasn't defined basically I wasn't using the Random Package in my code so I added it. Rest this works pretty fine here's the output

image

@ChrisRackauckas
Copy link
Member

Well is this is fine then it hasn't shown what the source of the bug is then...

@ParamThakkar123
Copy link
Contributor Author

Well is this is fine then it hasn't shown what the source of the bug is then...

The source of the bug is the NeuralPDE package itself, I think because without using it everything seems to work without any error. I just removed everything related to NeuralPDE and it started working

@ChrisRackauckas
Copy link
Member

NeuralPDE.jl just makes a loss function. You can make a loss function with the same issue without NeuralPDE

@ParamThakkar123
Copy link
Contributor Author

@ChrisRackauckas I made the following changes to the existing allen_cahn code and everything seems to work fine. Below is the code :

function allen_cahn(strategy, minimizer, maxIters)
    ##  DECLARATIONS
    @parameters t x1 x2 x3 x4
    @variables u(..)

    Dt = Differential(t)
    Dxx1 = Differential(x1)^2
    Dxx2 = Differential(x2)^2
    Dxx3 = Differential(x3)^2
    Dxx4 = Differential(x4)^2

    # Discretization
    tmax = 1.0
    x1width = x2width = x3width = x4width = 1.0
    tMeshNum = x1MeshNum = x2MeshNum = x3MeshNum = x4MeshNum = 10

    dt = tmax / tMeshNum
    dx1 = x1width / x1MeshNum
    dx2 = x2width / x2MeshNum
    dx3 = x3width / x3MeshNum
    dx4 = x4width / x4MeshNum

    domains = [t  Interval(0.0, tmax),
              x1  Interval(0.0, x1width),
              x2  Interval(0.0, x2width),
              x3  Interval(0.0, x3width),
              x4  Interval(0.0, x4width)]

    # Define the coordinates
    ts = 0.0:dt:tmax
    x1s = 0.0:dx1:x1width
    x2s = 0.0:dx2:x2width
    x3s = 0.0:dx3:x3width
    x4s = 0.0:dx4:x4width

    # Operators
    Δu = Dxx1(u(t,x1,x2,x3,x4)) + Dxx2(u(t,x1,x2,x3,x4)) + 
         Dxx3(u(t,x1,x2,x3,x4)) + Dxx4(u(t,x1,x2,x3,x4))

    # Equation
    eq = Dt(u(t,x1,x2,x3,x4)) - Δu - u(t,x1,x2,x3,x4) + 
         u(t,x1,x2,x3,x4)^3 ~ 0

    # Initial condition
    initialCondition = 1 / (2 + 0.4 * (x1^2 + x2^2 + x3^2 + x4^2))
    bcs = [u(0,x1,x2,x3,x4) ~ initialCondition]

    # Neural Network setup
    n = 10
    
    # Create the neural network with simplified architecture
    chain = Lux.Chain(
        Lux.Dense(5 => n, tanh),
        Lux.Dense(n => n, tanh),
        Lux.Dense(n => 1, identity)
    )

    # Modify training strategy based on type
    if strategy isa QuadratureTraining
        if strategy.quadrature_alg isa CubaCuhre
            # Use modified settings for CubaCuhre
            strategy = NeuralPDE.QuadratureTraining(
                quadrature_alg = CubaCuhre(),
                reltol = 1e-4,
                abstol = 1e-4,
                maxiters = 1100,
                batch = 100  # Enable batch processing
            )
        elseif strategy.quadrature_alg isa HCubatureJL
            # Switch to StochasticTraining for high dimensions
            @warn "Switching to StochasticTraining for better performance in high dimensions"
            strategy = NeuralPDE.StochasticTraining(
                400,
                bcs_points = 50
            )
        end
    end
    
    # Create the PINN with modified discretization
    discretization = PhysicsInformedNN(
        chain,
        strategy;
        init_params = nothing
    )

    # Setup PDE system with explicit variable ordering
    @named pde_system = PDESystem(
        eq,
        bcs,
        domains,
        [t,x1,x2,x3,x4],
        [u(t,x1,x2,x3,x4)]
    )

    # Discretize with fallback options
    prob = try
        discretize(pde_system, discretization)
    catch e
        @warn "Initial discretization failed, trying StochasticTraining"
        # Fall back to StochasticTraining
        strategy = NeuralPDE.StochasticTraining(
            400,
            bcs_points = 50
        )
        discretization = PhysicsInformedNN(chain, strategy)
        discretize(pde_system, discretization)
    end

    # Training setup
    losses = Float64[]
    error = Float64[]
    times = Float64[]
    timeCounter = 0.0
    startTime = time_ns()

    # Callback function with error handling
    function cb(p, l)
        try
            deltaT_s = time_ns()
            ctime = time_ns() - startTime - timeCounter
            
            push!(times, ctime / 1e9)
            push!(losses, l)
            push!(error, l)
            
            timeCounter += time_ns() - deltaT_s
            return false
        catch e
            @warn "Callback error: $e"
            return false
        end
    end

    # Solve with modified optimization parameters
    res = Optimization.solve(
        prob,
        minimizer;
        callback = cb,
        maxiters = maxIters
    )

    # Generate predictions
    phi = discretization.phi
    domain = [ts, x1s, x2s, x3s, x4s]
    
    u_predict = try
        [reshape([first(phi([t,x1,x2,x3,x4], res.minimizer)) 
                 for x1 in x1s, x2 in x2s, x3 in x3s, x4 in x4s],
                 (length(x1s), length(x2s), length(x3s), length(x4s))) 
                 for t in ts]
    catch e
        @warn "Prediction generation failed: $e"
        nothing
    end

    return [error, res.minimizer, domain, times, losses]
end

The output:
image

@ChrisRackauckas
Copy link
Member

what's the change?

@ParamThakkar123
Copy link
Contributor Author

  • Used Lux.glorot_uniform for explicit initialization.
  • Ensured correct initialization of parameters using Lux.setup.
  • Updated the layer construction to use the modern arrow syntax (=>).
  • Removed explicit parameter initialization once Lux was handling it properly.
  • Added try-catch blocks for critical operations.
  • Implemented fallback mechanisms:
  • Switched to a simpler quadrature algorithm if the initial attempt failed.
  • Relaxed optimization tolerances when necessary.
  • Introduced warning messages for better debugging.
  • Modified QuadratureTraining strategy:
  • Added error handling.
  • Introduced a batch parameter.
  • Implemented a fallback to StochasticTraining if QuadratureTraining failed.
  • Explicitly set batch sizes for QuadratureTraining.
  • Enabled batch processing in the optimization solver.
  • Added batch support to CubaCuhre for numerical integration.
  • Adjusted tolerances for better convergence.

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

Successfully merging this pull request may close these issues.

4 participants