From 8f88211c1e8be4ae53325e02533b6613d69da873 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Wed, 10 Apr 2024 01:46:48 +0000 Subject: [PATCH 01/22] Update Project with new dependencies --- Project.toml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a6a29c2..f4bdf00 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,8 @@ version = "0.1.1" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" +FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Lux = "b2108857-7c20-44ae-9111-449ecde12c47" @@ -19,6 +21,7 @@ PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] @@ -26,7 +29,6 @@ BenchmarkTools = "1" ComponentArrays = "0.15" Distributions = "0.25" Infiltrator = "1.2" -Lux = "0.5" Optimization = "3.12" OptimizationOptimJL = "0.1.5" OptimizationOptimisers = "0.1.2" From ad54394e8b9751de95f20524715618170176876a Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Wed, 17 Apr 2024 02:26:47 +0000 Subject: [PATCH 02/22] remove FiniteDifferences from dependencies --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index f4bdf00..6557314 100644 --- a/Project.toml +++ b/Project.toml @@ -8,7 +8,6 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" -FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Lux = "b2108857-7c20-44ae-9111-449ecde12c47" From a4790d3d4dc52d52da1ea53a06f5a8a0a0779dbf Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Wed, 17 Apr 2024 19:03:09 +0000 Subject: [PATCH 03/22] Example of APWP fit based on Jupp1987 --- Project.toml | 2 + .../APWP_Jupp1987/Jupp-etal-1987_dataset.csv | 27 ++++++++++++ examples/APWP_Jupp1987/apwp.jl | 41 +++++++++++++++++++ src/plot.jl | 4 +- src/train.jl | 24 +++++------ src/utils.jl | 33 ++++++++++++++- 6 files changed, 114 insertions(+), 17 deletions(-) create mode 100644 examples/APWP_Jupp1987/Jupp-etal-1987_dataset.csv create mode 100644 examples/APWP_Jupp1987/apwp.jl diff --git a/Project.toml b/Project.toml index 6557314..968e8a7 100644 --- a/Project.toml +++ b/Project.toml @@ -5,7 +5,9 @@ version = "0.1.1" [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" diff --git a/examples/APWP_Jupp1987/Jupp-etal-1987_dataset.csv b/examples/APWP_Jupp1987/Jupp-etal-1987_dataset.csv new file mode 100644 index 0000000..322e2d4 --- /dev/null +++ b/examples/APWP_Jupp1987/Jupp-etal-1987_dataset.csv @@ -0,0 +1,27 @@ +Obs,Time,Lat,Lon +A,0.3,80.6,349.2 +B,1,86.2,123.7 +C,1,86.6,205.5 +D,14.9,83.3,294 +E,26.6,84.8,67.6 +F,32.5,81,41 +G,55.5,86,178 +H,68.9,76,147 +I,98,86,298 +J,106,70.8,15.4 +K,107.8,61.5,244.2 +L,109.6,87,49.5 +M,118,36,296.3 +N,144,86,344 +O,171.9,56.5,12 +P,172,58,38 +Q,172,45,39 +R,172,59,41 +S,189,54.2,40.2 +T,189,44.1,51.5 +U,189,53.8,42.6 +V,209.5,82,290 +W,249.5,51,206 +X,481.5,9.3,206.7 +Y,481.5,28,190 +Z,519.7,1.5,208.5 \ No newline at end of file diff --git a/examples/APWP_Jupp1987/apwp.jl b/examples/APWP_Jupp1987/apwp.jl new file mode 100644 index 0000000..df03292 --- /dev/null +++ b/examples/APWP_Jupp1987/apwp.jl @@ -0,0 +1,41 @@ +using Pkg; Pkg.activate(".") +using Revise + +using SphereUDE +using DataFrames, CSV +using Random +rng = Random.default_rng() +Random.seed!(rng, 666) + + +df = CSV.read("examples/APWP_Jupp1987/Jupp-etal-1987_dataset.csv", DataFrame; header=true) + +times = df.Time +times .+= rand(sampler(Normal(0,0.1)), length(times)) # Needs to fix this! +X = sph2cart(Matrix(df[:,["Lat","Lon"]])'; radians=false) + +data = SphereData(times=times, directions=X, kappas=nothing, L=nothing) + +# Expected maximum angular deviation in one unit of time (degrees) +Δω₀ = 1.0 +# Angular velocity +ω₀ = Δω₀ * π / 180.0 + +# regs = [Regularization(order=1, power=1.0, λ=0.001, diff_mode="Finite Differences"), +# Regularization(order=0, power=2.0, λ=0.1, diff_mode="Finite Differences")] +# regs = [Regularization(order=1, power=2.0, λ=0.1, diff_mode="Finite Differences")] +regs = nothing + +params = SphereParameters(tmin=0., tmax=520., + reg=regs, + u0=[0.0, 0.0, 1.0], ωmax=ω₀, reltol=1e-7, abstol=1e-7, + niter_ADAM=1000, niter_LBFGS=400) + +results = train(data, params, rng, nothing) + +############################################################## +###################### PyCall Plots ######################### +############################################################## + +plot_sphere(data, results, mean(df.Lat), mean(df.Lon), saveas="examples/APWP_Jupp1987/plot_sphere.pdf", title="Double rotation") +plot_L(data, results, saveas="examples/APWP_Jupp1987/plot_L.pdf", title="Double rotation") \ No newline at end of file diff --git a/src/plot.jl b/src/plot.jl index e37d03b..bbf3697 100644 --- a/src/plot.jl +++ b/src/plot.jl @@ -30,7 +30,7 @@ function plot_sphere(# ax::PyCall.PyObject, plt.figure(figsize=(10,10)) ax = plt.axes(projection=ccrs.Orthographic(central_latitude=central_latitude, central_longitude=central_longitude)) - ax.coastlines() + # ax.coastlines() ax.gridlines() ax.set_global() @@ -39,7 +39,7 @@ function plot_sphere(# ax::PyCall.PyObject, X_fit_path = cart2sph(results.fit_directions, radians=false) sns.scatterplot(ax=ax, x = X_true_points[1,:], y=X_true_points[2, :], - hue = data.times, s=50, + hue = data.times, s=150, palette="viridis", transform = ccrs.PlateCarree()); diff --git a/src/train.jl b/src/train.jl index 641b72e..1b846df 100644 --- a/src/train.jl +++ b/src/train.jl @@ -1,11 +1,15 @@ export train +# For L1 regularization relu_cap works better, but for L2 I think is better to include sigmoid function get_NN(params, rng, θ_trained) # Define neural network U = Lux.Chain( Lux.Dense(1, 5, relu_cap), # explore discontinuity function for activation Lux.Dense(5, 10, relu_cap), - Lux.Dense(10, 5, relu_cap), + Lux.Dense(10, 5, relu_cap), + # Lux.Dense(1, 5, sigmoid), + # Lux.Dense(5, 10, sigmoid), + # Lux.Dense(10, 5, sigmoid), Lux.Dense(5, 3, x->sigmoid_cap(x; ω₀=params.ωmax)) ) θ, st = Lux.setup(rng, U) @@ -19,9 +23,6 @@ function train(data::AD, U, θ, st = get_NN(params, rng, θ_trained) - # one option is to restrict where the NN is evaluated to discrete t to - # generate piece-wise dynamics. - function ude_rotation!(du, u, p, t) # Angular momentum given by network prediction L = U([t], p, st)[1] @@ -66,17 +67,14 @@ function train(data::AD, # Create (uniform) spacing time # Δt = (params.tmax - params.tmin) / n_nodes # times_reg = collect(params.tmin:Δt:params.tmax) - # LinRange does not propagate thought the backprop step! # times_reg = collect(LinRange(params.tmin, params.tmax, n_nodes)) l_ = 0.0 if reg.order==0 l_ += quadrature(t -> norm(U([t], θ, st)[1])^reg.power, params.tmin, params.tmax, n_nodes) - # for t in times_reg - # l_ += norm(U([t], θ, st)[1])^reg.power - # end elseif reg.order==1 if reg.diff_mode=="AD" + throw("Method not working well.") # Compute gradient using automatic differentiaion in the NN # This currently doesn't run... too slow. for t in times_reg @@ -85,10 +83,10 @@ function train(data::AD, l_ += norm(grad)^reg.power end elseif reg.diff_mode=="Finite Differences" - # Compute finite differences - L_estimated = map(t -> (first ∘ U)([t], θ, st), times_reg) - dLdt = diff(L_estimated) ./ diff(times_reg) - l_ += sum(norm.(dLdt).^reg.power) + # Using FiniteDifferences break precompilation becuase of name collission + # l_ += quadrature(t -> norm(FiniteDifferences.jacobian(FiniteDifferences.central_fdm(2,1), τ -> (first ∘ U)([τ], θ, st), t)[1])^reg.power, params.tmin, params.tmax, n_nodes) + ϵ = 0.1 * (params.tmax - params.tmin) / n_nodes + l_ += quadrature(t -> norm(central_fdm(τ -> (first ∘ U)([τ], θ, st), t, ϵ=ϵ))^reg.power, params.tmin, params.tmax, n_nodes) else throw("Method not implemented.") end @@ -122,7 +120,7 @@ function train(data::AD, θ_trained = res2.u # Final Fit - fit_times = collect(params.tmin:0.1:params.tmax) + fit_times = collect(range(params.tmin,params.tmax, length=200)) fit_directions = predict(θ_trained, T=fit_times) return Results(θ_trained=θ_trained, U=U, st=st, diff --git a/src/utils.jl b/src/utils.jl index d441b79..9dfc1fb 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,7 +1,7 @@ export sigmoid_cap, relu_cap, step_cap -export cart2sph +export cart2sph, sph2cart export AbstractNoise, FisherNoise -export quadrature +export quadrature, central_fdm # Normalization of the NN. Ideally we want to do this with L2 norm . @@ -37,6 +37,23 @@ function cart2sph(X::AbstractArray{<:Number}; radians::Bool=true) return Y end + +""" + sph2cart(X::AbstractArray{<:Number}; radians::Bool=true) + +Convert spherical coordinates to cartesian +""" +function sph2cart(X::AbstractArray{<:Number}; radians::Bool=true) + @assert size(X)[1] == 2 "Input array must have two rows corresponding to Latitude and Longitude." + if !radians + X *= π / 180. + end + Y = mapslices(x -> [cos(x[1])*cos(x[2]), + cos(x[1])*sin(x[2]), + sin(x[1])] , X, dims=1) + return Y +end + """ Add Fisher noise to matrix of three dimensional unit vectors @@ -78,4 +95,16 @@ function quadrature(f::Function, t₀, t₁, n_nodes::Int) nodes = (t₀+t₁)/2 .+ nodes * (t₁-t₀)/2 weigths = (t₁-t₀) / 2 * weigths return dot(weigths, f.(nodes)) +end + +""" + central_fdm(f::Function, x::Float64; ϵ=0.01) + +Simple central differences implementation. + +FiniteDifferences.jl does not work with AD so I implemented this manually. +Still remains to test this with FiniteDiff.jl +""" +function central_fdm(f::Function, x::Float64; ϵ=0.01) + return (f(x+ϵ)-f(x-ϵ)) / (2ϵ) end \ No newline at end of file From 29b8c78450fa93d3a586db334734f75237d6d8c4 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Thu, 18 Apr 2024 19:46:37 -0700 Subject: [PATCH 04/22] add complex-step method --- src/train.jl | 14 +++++++++++--- src/utils.jl | 11 ++++++++++- 2 files changed, 21 insertions(+), 4 deletions(-) diff --git a/src/train.jl b/src/train.jl index 1b846df..647c33d 100644 --- a/src/train.jl +++ b/src/train.jl @@ -29,6 +29,12 @@ function train(data::AD, du .= cross(L, u) end + # function ude_rotation!(du::Array{Complex{Float64}}, u::Array{Complex{Float64}}, p, t) + # # Angular momentum given by network prediction + # L = U([t], p, st)[1] + # du .= cross(L, u) + # end + prob_nn = ODEProblem(ude_rotation!, params.u0, [params.tmin, params.tmax], θ) function predict(θ::ComponentVector; u0=params.u0, T=data.times) @@ -82,11 +88,13 @@ function train(data::AD, grad = Zygote.jacobian(first ∘ U, [t], θ, st)[1] l_ += norm(grad)^reg.power end - elseif reg.diff_mode=="Finite Differences" - # Using FiniteDifferences break precompilation becuase of name collission - # l_ += quadrature(t -> norm(FiniteDifferences.jacobian(FiniteDifferences.central_fdm(2,1), τ -> (first ∘ U)([τ], θ, st), t)[1])^reg.power, params.tmin, params.tmax, n_nodes) + elseif reg.diff_mode=="FD" + # Finite differences ϵ = 0.1 * (params.tmax - params.tmin) / n_nodes l_ += quadrature(t -> norm(central_fdm(τ -> (first ∘ U)([τ], θ, st), t, ϵ=ϵ))^reg.power, params.tmin, params.tmax, n_nodes) + elseif reg.diff_mode=="CS" + # Complex step differentiation + l_ += quadrature(t -> norm(complex_step_differentiation(τ -> (first ∘ U)([τ], θ, st), t))^reg.power, params.tmin, params.tmax, n_nodes) else throw("Method not implemented.") end diff --git a/src/utils.jl b/src/utils.jl index 9dfc1fb..9ad00c9 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,7 +1,7 @@ export sigmoid_cap, relu_cap, step_cap export cart2sph, sph2cart export AbstractNoise, FisherNoise -export quadrature, central_fdm +export quadrature, central_fdm, complex_step_differentiation # Normalization of the NN. Ideally we want to do this with L2 norm . @@ -107,4 +107,13 @@ Still remains to test this with FiniteDiff.jl """ function central_fdm(f::Function, x::Float64; ϵ=0.01) return (f(x+ϵ)-f(x-ϵ)) / (2ϵ) +end + +""" + complex_step_differentiation(f::Function, x::Float64; ϵ=1e-10) + +Manual implementation of comple-step differentiation +""" +function complex_step_differentiation(f::Function, x::Float64; ϵ=1e-10) + return imag(f(x + ϵ * im)) / ϵ end \ No newline at end of file From a3852d16258370884346cab348700b6b104262a5 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Thu, 18 Apr 2024 23:10:00 -0700 Subject: [PATCH 05/22] Double differentiation working with complex-step method --- examples/double_rotation/double_rotation.jl | 6 ++-- src/SphereUDE.jl | 5 +++- src/train.jl | 11 +++---- src/utils.jl | 32 ++++++++++++++++++++- 4 files changed, 42 insertions(+), 12 deletions(-) diff --git a/examples/double_rotation/double_rotation.jl b/examples/double_rotation/double_rotation.jl index 1e0e5b8..d23e5b0 100644 --- a/examples/double_rotation/double_rotation.jl +++ b/examples/double_rotation/double_rotation.jl @@ -64,9 +64,9 @@ X_true = X_noiseless + FisherNoise(kappa=200.) data = SphereData(times=times_samples, directions=X_true, kappas=nothing, L=L_true) -# regs = [Regularization(order=1, power=1.0, λ=0.001, diff_mode="Finite Differences"), - # Regularization(order=0, power=2.0, λ=0.1, diff_mode="Finite Differences")] -regs = [Regularization(order=0, power=2.0, λ=0.1, diff_mode="AD")] +regs = [Regularization(order=1, power=1.0, λ=0.001, diff_mode="CS"), + Regularization(order=0, power=2.0, λ=0.1, diff_mode="AD")] +# regs = [Regularization(order=1, power=1.0, λ=0.1, diff_mode="CS")] params = SphereParameters(tmin=tspan[1], tmax=tspan[2], reg=regs, diff --git a/src/SphereUDE.jl b/src/SphereUDE.jl index ac216c0..ebce6ae 100644 --- a/src/SphereUDE.jl +++ b/src/SphereUDE.jl @@ -14,11 +14,14 @@ using Optimization, OptimizationOptimisers, OptimizationOptimJL using ComponentArrays: ComponentVector using PyPlot, PyCall +# Testing double-differentiation +# using BatchedRoutines + # Debugging using Infiltrator -include("utils.jl") include("types.jl") +include("utils.jl") include("train.jl") include("plot.jl") diff --git a/src/train.jl b/src/train.jl index 647c33d..d0f561a 100644 --- a/src/train.jl +++ b/src/train.jl @@ -21,6 +21,9 @@ function train(data::AD, rng, θ_trained=[]) where{AD <: AbstractData, AP <: AbstractParameters} + # Raise warnings + raise_warnings(data::AD, params::AP) + U, θ, st = get_NN(params, rng, θ_trained) function ude_rotation!(du, u, p, t) @@ -29,12 +32,6 @@ function train(data::AD, du .= cross(L, u) end - # function ude_rotation!(du::Array{Complex{Float64}}, u::Array{Complex{Float64}}, p, t) - # # Angular momentum given by network prediction - # L = U([t], p, st)[1] - # du .= cross(L, u) - # end - prob_nn = ODEProblem(ude_rotation!, params.u0, [params.tmin, params.tmax], θ) function predict(θ::ComponentVector; u0=params.u0, T=data.times) @@ -107,7 +104,7 @@ function train(data::AD, losses = Float64[] callback = function (p, l) push!(losses, l) - if length(losses) % 200 == 0 + if length(losses) % 100 == 0 println("Current loss after $(length(losses)) iterations: $(losses[end])") end return false diff --git a/src/utils.jl b/src/utils.jl index 9ad00c9..c8f2eed 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -2,6 +2,7 @@ export sigmoid_cap, relu_cap, step_cap export cart2sph, sph2cart export AbstractNoise, FisherNoise export quadrature, central_fdm, complex_step_differentiation +export raise_warnings # Normalization of the NN. Ideally we want to do this with L2 norm . @@ -14,6 +15,10 @@ function sigmoid_cap(x; ω₀=1.0) return min_value + (max_value - min_value) / ( 1.0 + exp(-x) ) end +function sigmoid(x::Complex) + return 1 / ( 1.0 + exp(-x) ) +end + """ relu_cap(x; ω₀=1.0) """ @@ -23,6 +28,12 @@ function relu_cap(x; ω₀=1.0) return min_value + (max_value - min_value) * max(0.0, min(x, 1.0)) end +function relu_cap(x::Complex; ω₀=1.0) + min_value = - ω₀ + max_value = + ω₀ + return min_value + (max_value - min_value) * max(0.0, min(real(x), 1.0)) + (max_value - min_value) * max(0.0, min(imag(x), 1.0)) * im +end + """ cart2sph(X::AbstractArray{<:Number}; radians::Bool=true) @@ -92,7 +103,7 @@ function quadrature(f::Function, t₀, t₁, n_nodes::Int) # Ignore AD here since FastGaussQuadrature is using mutating arrays nodes, weigths = gausslegendre(n_nodes) end - nodes = (t₀+t₁)/2 .+ nodes * (t₁-t₀)/2 + nodes = (t₀+t₁)/2 .+ nodes * (t₁-t₀)/2 weigths = (t₁-t₀) / 2 * weigths return dot(weigths, f.(nodes)) end @@ -116,4 +127,23 @@ Manual implementation of comple-step differentiation """ function complex_step_differentiation(f::Function, x::Float64; ϵ=1e-10) return imag(f(x + ϵ * im)) / ϵ +end + +""" + raise_warnings(data::AD, params::AP) + +Raise warnings. +""" +function raise_warnings(data::SphereData, params::SphereParameters) + if length(unique(data.times)) < length(data.times) + @warn "[SphereUDE] Timeseries includes duplicated times. \n This can produce unexpected errors." + end + if !isnothing(params.reg) + for reg in params.reg + if reg.diff_mode=="CS" + @warn "[SphereUDE] Complex-step differentiation inside the loss function \n This just work for cases where the activation function of the neural network admits complex numbers \n Change predefined activation functions to be complex numbers." + end + end + end + nothing end \ No newline at end of file From cc5ce31b3d82d1109e44e7d397fb4727066c494e Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Fri, 19 Apr 2024 12:53:51 -0700 Subject: [PATCH 06/22] Initial condition u0 fitting implemented --- examples/double_rotation/double_rotation.jl | 13 +++-- src/plot.jl | 2 +- src/train.jl | 56 +++++++++++++-------- src/types.jl | 3 +- 4 files changed, 47 insertions(+), 27 deletions(-) diff --git a/examples/double_rotation/double_rotation.jl b/examples/double_rotation/double_rotation.jl index d23e5b0..04ddc3f 100644 --- a/examples/double_rotation/double_rotation.jl +++ b/examples/double_rotation/double_rotation.jl @@ -64,16 +64,19 @@ X_true = X_noiseless + FisherNoise(kappa=200.) data = SphereData(times=times_samples, directions=X_true, kappas=nothing, L=L_true) -regs = [Regularization(order=1, power=1.0, λ=0.001, diff_mode="CS"), - Regularization(order=0, power=2.0, λ=0.1, diff_mode="AD")] -# regs = [Regularization(order=1, power=1.0, λ=0.1, diff_mode="CS")] +# regs = [Regularization(order=1, power=1.0, λ=0.001, diff_mode="CS"), +# Regularization(order=0, power=2.0, λ=0.1, diff_mode="AD")] +regs = [Regularization(order=0, power=2.0, λ=0.1, diff_mode=nothing), + Regularization(order=1, power=1.1, λ=0.01, diff_mode="CS")] +# regs = nothing params = SphereParameters(tmin=tspan[1], tmax=tspan[2], reg=regs, u0=[0.0, 0.0, -1.0], ωmax=ω₀, reltol=reltol, abstol=abstol, - niter_ADAM=1000, niter_LBFGS=600) + niter_ADAM=1000, niter_LBFGS=500) -results = train(data, params, rng, nothing) +# results = train(data, params, rng, nothing) +results = train(data, params, rng, nothing; train_initial_condition=true) ############################################################## ###################### PyCall Plots ######################### diff --git a/src/plot.jl b/src/plot.jl index bbf3697..8fafd15 100644 --- a/src/plot.jl +++ b/src/plot.jl @@ -69,7 +69,7 @@ function plot_L(data::AbstractData, fig, ax = plt.subplots(figsize=(10,5)) times_smooth = collect(LinRange(results.fit_times[begin], results.fit_times[end], 1000)) - Ls = reduce(hcat, (t -> results.U([t], results.θ_trained, results.st)[1]).(times_smooth)) + Ls = reduce(hcat, (t -> results.U([t], results.θ, results.st)[1]).(times_smooth)) angular_velocity = mapslices(x -> norm(x), Ls, dims=1)[1,:] diff --git a/src/train.jl b/src/train.jl index d0f561a..60a837a 100644 --- a/src/train.jl +++ b/src/train.jl @@ -19,32 +19,46 @@ end function train(data::AD, params::AP, rng, - θ_trained=[]) where{AD <: AbstractData, AP <: AbstractParameters} + θ_trained=[]; + train_initial_condition::Bool=false) where {AD <: AbstractData, AP <: AbstractParameters} # Raise warnings raise_warnings(data::AD, params::AP) U, θ, st = get_NN(params, rng, θ_trained) + # Set component vector for Optimization + if train_initial_condition + β = ComponentVector{Float64}(θ=θ, u0=params.u0) + else + β = ComponentVector{Float64}(θ=θ) + end + function ude_rotation!(du, u, p, t) # Angular momentum given by network prediction L = U([t], p, st)[1] du .= cross(L, u) end - prob_nn = ODEProblem(ude_rotation!, params.u0, [params.tmin, params.tmax], θ) + prob_nn = ODEProblem(ude_rotation!, params.u0, [params.tmin, params.tmax], β.θ) - function predict(θ::ComponentVector; u0=params.u0, T=data.times) - _prob = remake(prob_nn, u0=u0, - tspan=(min(T[1], params.tmin), max(T[end], params.tmax)), - p = θ) + function predict(β::ComponentVector; T=data.times) + if train_initial_condition + _prob = remake(prob_nn, u0=β.u0 / norm(β.u0), + tspan=(min(T[1], params.tmin), max(T[end], params.tmax)), + p = β.θ) + else + _prob = remake(prob_nn, u0=params.u0, + tspan=(min(T[1], params.tmin), max(T[end], params.tmax)), + p = β.θ) + end Array(solve(_prob, params.solver, saveat=T, abstol=params.abstol, reltol=params.reltol, sensealg=params.sensealg)) end - function loss(θ::ComponentVector) - u_ = predict(θ) + function loss(β::ComponentVector) + u_ = predict(β) # Empirical error # l_emp = mean(abs2, u_ .- data.directions) if isnothing(data.kappas) @@ -58,8 +72,7 @@ function train(data::AD, if !isnothing(params.reg) # for (order, power, λ) in params.reg for reg in params.reg - # l_reg += reg.λ * regularization(θ; order=reg.order, power=reg.power) - l_reg += regularization(θ, reg) + l_reg += regularization(β.θ, reg) end end return l_emp + l_reg @@ -67,11 +80,6 @@ function train(data::AD, function regularization(θ::ComponentVector, reg::AbstractRegularization; n_nodes=100) - # Create (uniform) spacing time - # Δt = (params.tmax - params.tmin) / n_nodes - # times_reg = collect(params.tmin:Δt:params.tmax) - # LinRange does not propagate thought the backprop step! - # times_reg = collect(LinRange(params.tmin, params.tmax, n_nodes)) l_ = 0.0 if reg.order==0 l_ += quadrature(t -> norm(U([t], θ, st)[1])^reg.power, params.tmin, params.tmax, n_nodes) @@ -111,8 +119,8 @@ function train(data::AD, end adtype = Optimization.AutoZygote() - optf = Optimization.OptimizationFunction((x, θ) -> loss(x), adtype) - optprob = Optimization.OptimizationProblem(optf, ComponentVector{Float64}(θ)) + optf = Optimization.OptimizationFunction((x, β) -> loss(x), adtype) + optprob = Optimization.OptimizationProblem(optf, β) res1 = Optimization.solve(optprob, ADAM(0.002), callback=callback, maxiters=params.niter_ADAM) println("Training loss after $(length(losses)) iterations: $(losses[end])") @@ -122,12 +130,20 @@ function train(data::AD, println("Final training loss after $(length(losses)) iterations: $(losses[end])") # Optimized NN parameters - θ_trained = res2.u + β_trained = res2.u + θ_trained = β_trained.θ + + # Optimized initial condition + if train_initial_condition + u0_trained = Array(β_trained.u0) # β.u0 is a view type + else + u0_trained = params.u0 + end # Final Fit fit_times = collect(range(params.tmin,params.tmax, length=200)) - fit_directions = predict(θ_trained, T=fit_times) + fit_directions = predict(β_trained, T=fit_times) - return Results(θ_trained=θ_trained, U=U, st=st, + return Results(θ=θ_trained, u0=u0_trained, U=U, st=st, fit_times=fit_times, fit_directions=fit_directions) end diff --git a/src/types.jl b/src/types.jl index b1c3001..a8c6237 100644 --- a/src/types.jl +++ b/src/types.jl @@ -39,7 +39,8 @@ end Final results """ @kwdef struct Results{F <: AbstractFloat} <: AbstractResult - θ_trained::ComponentVector + θ::ComponentVector + u0::Vector{F} U::Lux.Chain st::NamedTuple fit_times::Vector{F} From 92f00d6cc70ccd737f43d9cc21ab48a6176211a1 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Fri, 19 Apr 2024 17:48:50 -0700 Subject: [PATCH 07/22] Projected gradient descent working for u0. Some more tests --- src/plot.jl | 8 +++++--- src/train.jl | 14 +++++++++----- src/utils.jl | 6 +++--- test/runtests.jl | 5 ++++- test/utils.jl | 22 +++++++++++++++++++--- 5 files changed, 40 insertions(+), 15 deletions(-) diff --git a/src/plot.jl b/src/plot.jl index 8fafd15..7deb758 100644 --- a/src/plot.jl +++ b/src/plot.jl @@ -38,14 +38,16 @@ function plot_sphere(# ax::PyCall.PyObject, # X_true_path = cart2sph(X_path, radians=false) X_fit_path = cart2sph(results.fit_directions, radians=false) - sns.scatterplot(ax=ax, x = X_true_points[1,:], y=X_true_points[2, :], + # Plots in Python follow the long, lat ordering + + sns.scatterplot(ax=ax, x = X_true_points[2,:], y=X_true_points[1, :], hue = data.times, s=150, palette="viridis", transform = ccrs.PlateCarree()); for i in 1:(length(results.fit_times)-1) - plt.plot([X_fit_path[1,i], X_fit_path[1,i+1]], - [X_fit_path[2,i], X_fit_path[2,i+1]], + plt.plot([X_fit_path[2,i], X_fit_path[2,i+1]], + [X_fit_path[1,i], X_fit_path[1,i+1]], linewidth=2, color="black",#cmap(norm(results.fit_times[i])), transform = ccrs.Geodetic()) end diff --git a/src/train.jl b/src/train.jl index 60a837a..2587e1d 100644 --- a/src/train.jl +++ b/src/train.jl @@ -44,7 +44,7 @@ function train(data::AD, function predict(β::ComponentVector; T=data.times) if train_initial_condition - _prob = remake(prob_nn, u0=β.u0 / norm(β.u0), + _prob = remake(prob_nn, u0=β.u0 / norm(β.u0), # We enforced the norm=1 condition again here tspan=(min(T[1], params.tmin), max(T[end], params.tmax)), p = β.θ) else @@ -60,12 +60,13 @@ function train(data::AD, function loss(β::ComponentVector) u_ = predict(β) # Empirical error - # l_emp = mean(abs2, u_ .- data.directions) if isnothing(data.kappas) + l_emp = 3.0 * mean(abs2.(u_ .- data.directions)) # The 3 is needed since the mean is computen on a 3xN matrix - l_emp = 1 - 3.0 * mean(u_ .* data.directions) + # l_emp = 1 - 3.0 * mean(u_ .* data.directions) else - l_emp = norm(data.kappas)^2 - 3.0 * mean(data.kappas .* u_ .* data.directions) + # l_emp = norm(data.kappas)^2 - 3.0 * mean(data.kappas .* u_ .* data.directions) + l_emp = mean(data.kappas .* abs2.(u_ .- data.directions), dims=1) end # Regularization l_reg = 0.0 @@ -112,9 +113,12 @@ function train(data::AD, losses = Float64[] callback = function (p, l) push!(losses, l) - if length(losses) % 100 == 0 + if length(losses) % 10 == 0 println("Current loss after $(length(losses)) iterations: $(losses[end])") end + if train_initial_condition + p.u0 ./= norm(p.u0) + end return false end diff --git a/src/utils.jl b/src/utils.jl index c8f2eed..4ca6894 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,4 +1,4 @@ -export sigmoid_cap, relu_cap, step_cap +export sigmoid_cap, sigmoid, relu_cap, step_cap export cart2sph, sph2cart export AbstractNoise, FisherNoise export quadrature, central_fdm, complex_step_differentiation @@ -41,7 +41,7 @@ Convert cartesian coordinates to spherical """ function cart2sph(X::AbstractArray{<:Number}; radians::Bool=true) @assert size(X)[1] == 3 "Input array must have three rows." - Y = mapslices(x -> [angle(x[1] + x[2]*im), asin(x[3])] , X, dims=1) + Y = mapslices(x -> [asin(x[3]), angle(x[1] + x[2]*im)] , X, dims=1) if !radians Y *= 180. / π end @@ -61,7 +61,7 @@ function sph2cart(X::AbstractArray{<:Number}; radians::Bool=true) end Y = mapslices(x -> [cos(x[1])*cos(x[2]), cos(x[1])*sin(x[2]), - sin(x[1])] , X, dims=1) + sin(x[1])], X, dims=1) return Y end diff --git a/test/runtests.jl b/test/runtests.jl index f0b11a1..f33bf74 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,6 @@ using SphereUDE using Test +using Lux include("constructors.jl") include("utils.jl") @@ -11,5 +12,7 @@ include("utils.jl") end @testset "Utils" begin - test_cart2sph() + test_coordinate() + test_complex_activation() + test_integration() end diff --git a/test/utils.jl b/test/utils.jl index e12e0e9..172ba3b 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,6 +1,22 @@ -function test_cart2sph() - X₀ = [1 0 0 √2; 0 1 0 √2; 0 0 1 0] - Y₀ = [0.0 90.0 0.0 45.0; 0.0 0.0 90.0 0.0] +function test_coordinate() + X₀ = [1 0 0 1/√2; 0 0 1 1/√2; 0 1 0 0] + Y₀ = [0.0 90.0 0.0 0.0; 0.0 0.0 90.0 45.0] @test all(isapprox.(Y₀, cart2sph(X₀, radians=false), atol=1e-6)) + @test all(isapprox.(Y₀ * π / 180., cart2sph(X₀, radians=true), atol=1e-6)) + @test all(isapprox.(X₀, sph2cart(Y₀, radians=false), atol=1e-6)) + @test all(isapprox.(X₀, sph2cart(Y₀ * π / 180., radians=true), atol=1e-6)) +end + +function test_complex_activation() + pure_real = 1.0 + pure_complex = 1.0 + 1.0 * im + @test isapprox(Lux.sigmoid(pure_real), 0.7310585786300049, atol=1e-6) + @test isapprox(imag(SphereUDE.sigmoid(pure_complex)), 0.2019482276580129, atol=1e-6) +end + +function test_integration() + @test isapprox(quadrature(x->1, 0.0, 1.0, 100), 1.0, rtol=1e-6) + @test isapprox(quadrature(x->x, -1.0, 1.0, 100), 0.0, atol=1e-6) + @test isapprox(quadrature(x->x^2, -1.0, 1.0, 100), 2/3., rtol=1e-6) end \ No newline at end of file From c69a0596f25ae3a3089b4d0ee42910d06e303c88 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Sun, 28 Apr 2024 19:34:32 -0700 Subject: [PATCH 08/22] Testing activation functions with complex-step --- Project.toml | 1 + examples/double_rotation/double_rotation.jl | 15 ++++---- src/train.jl | 30 +++++++++++---- src/utils.jl | 42 +++++++++++++++++---- 4 files changed, 66 insertions(+), 22 deletions(-) diff --git a/Project.toml b/Project.toml index 968e8a7..e59eec6 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +DiffEqFlux = "aae7a2af-3d4f-5e19-a356-7da93b79d9d0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" diff --git a/examples/double_rotation/double_rotation.jl b/examples/double_rotation/double_rotation.jl index 04ddc3f..f927deb 100644 --- a/examples/double_rotation/double_rotation.jl +++ b/examples/double_rotation/double_rotation.jl @@ -56,7 +56,7 @@ true_sol = solve(prob, Tsit5(), reltol=reltol, abstol=abstol, saveat=times_samp # Add Fisher noise to true solution X_noiseless = Array(true_sol) -X_true = X_noiseless + FisherNoise(kappa=200.) +X_true = X_noiseless + FisherNoise(kappa=50.) ############################################################## ####################### Training ########################### @@ -64,19 +64,18 @@ X_true = X_noiseless + FisherNoise(kappa=200.) data = SphereData(times=times_samples, directions=X_true, kappas=nothing, L=L_true) -# regs = [Regularization(order=1, power=1.0, λ=0.001, diff_mode="CS"), -# Regularization(order=0, power=2.0, λ=0.1, diff_mode="AD")] -regs = [Regularization(order=0, power=2.0, λ=0.1, diff_mode=nothing), - Regularization(order=1, power=1.1, λ=0.01, diff_mode="CS")] +regs = [Regularization(order=1, power=1.0, λ=1.0, diff_mode="CS"), + Regularization(order=0, power=2.0, λ=0.1, diff_mode=nothing)] +# regs = [Regularization(order=0, power=2.0, λ=0.1, diff_mode=nothing)] + # Regularization(order=1, power=1.1, λ=0.01, diff_mode="CS")] # regs = nothing params = SphereParameters(tmin=tspan[1], tmax=tspan[2], reg=regs, u0=[0.0, 0.0, -1.0], ωmax=ω₀, reltol=reltol, abstol=abstol, - niter_ADAM=1000, niter_LBFGS=500) + niter_ADAM=1000, niter_LBFGS=800) -# results = train(data, params, rng, nothing) -results = train(data, params, rng, nothing; train_initial_condition=true) +results = train(data, params, rng, nothing; train_initial_condition=false) ############################################################## ###################### PyCall Plots ######################### diff --git a/src/train.jl b/src/train.jl index 2587e1d..aec7ee7 100644 --- a/src/train.jl +++ b/src/train.jl @@ -4,18 +4,29 @@ export train function get_NN(params, rng, θ_trained) # Define neural network U = Lux.Chain( - Lux.Dense(1, 5, relu_cap), # explore discontinuity function for activation - Lux.Dense(5, 10, relu_cap), - Lux.Dense(10, 5, relu_cap), - # Lux.Dense(1, 5, sigmoid), - # Lux.Dense(5, 10, sigmoid), - # Lux.Dense(10, 5, sigmoid), + Lux.Dense(1, 5, sigmoid), # explore discontinuity function for activation + Lux.Dense(5, 10, sigmoid), + Lux.Dense(10, 5, sigmoid), + # Lux.Dense(1, 5, relu_cap), + # Lux.Dense(5, 10, relu_cap), + # Lux.Dense(10, 5, relu_cap), Lux.Dense(5, 3, x->sigmoid_cap(x; ω₀=params.ωmax)) ) θ, st = Lux.setup(rng, U) return U, θ, st end +""" + predict_L(t, NN, θ, st) + +Predict value of rotation given by L given by the neural network. + +% To do: replace all the calls in U by predict_L +""" +function predict_L(t, NN, θ, st) + return NN([t], θ, st)[1] +end + function train(data::AD, params::AP, rng, @@ -37,6 +48,7 @@ function train(data::AD, function ude_rotation!(du, u, p, t) # Angular momentum given by network prediction L = U([t], p, st)[1] + # L = predict_L(t, U, p, st) du .= cross(L, u) end @@ -65,8 +77,8 @@ function train(data::AD, # The 3 is needed since the mean is computen on a 3xN matrix # l_emp = 1 - 3.0 * mean(u_ .* data.directions) else - # l_emp = norm(data.kappas)^2 - 3.0 * mean(data.kappas .* u_ .* data.directions) l_emp = mean(data.kappas .* abs2.(u_ .- data.directions), dims=1) + # l_emp = norm(data.kappas)^2 - 3.0 * mean(data.kappas .* u_ .* data.directions) end # Regularization l_reg = 0.0 @@ -89,6 +101,10 @@ function train(data::AD, throw("Method not working well.") # Compute gradient using automatic differentiaion in the NN # This currently doesn't run... too slow. + + # Test this with the new implementation in Lux.jl: + # https://lux.csail.mit.edu/stable/manual/nested_autodiff + for t in times_reg # Try ReverseDiff grad = Zygote.jacobian(first ∘ U, [t], θ, st)[1] diff --git a/src/utils.jl b/src/utils.jl index 4ca6894..b8cb4f5 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,22 +1,38 @@ -export sigmoid_cap, sigmoid, relu_cap, step_cap +export sigmoid, sigmoid_cap +export relu, relu_cap export cart2sph, sph2cart export AbstractNoise, FisherNoise export quadrature, central_fdm, complex_step_differentiation export raise_warnings -# Normalization of the NN. Ideally we want to do this with L2 norm . """ sigmoid_cap(x; ω₀=1.0) + +Normalization of the neural network last layer """ function sigmoid_cap(x; ω₀=1.0) min_value = - ω₀ max_value = + ω₀ - return min_value + (max_value - min_value) / ( 1.0 + exp(-x) ) + return min_value + (max_value - min_value) * sigmoid(x) +end + +function sigmoid(x) + return 1.0 / (1.0 + exp(-x)) +# if x > 0 +# return 1 / ( 1.0 + exp(-x) ) +# else +# return exp(x) / (1.0 + exp(x)) +# end end -function sigmoid(x::Complex) - return 1 / ( 1.0 + exp(-x) ) +function sigmoid(z::Complex) + return 1.0 / ( 1.0 + exp(-z) ) + # if real(z) > 0 + # return 1 / ( 1.0 + exp(-z) ) + # else + # return exp(z) / (1.0 + exp(z)) + # end end """ @@ -28,10 +44,22 @@ function relu_cap(x; ω₀=1.0) return min_value + (max_value - min_value) * max(0.0, min(x, 1.0)) end -function relu_cap(x::Complex; ω₀=1.0) + +""" + relu(x::Complex) + +Extension of ReLU function to complex numbers based on the complex cardioid introduced in +Virtue et al. (2017), "Better than Real: Complex-valued Neural Nets for MRI Fingerprinting". +This function is equivalent to relu when x is real (and hence angle(x)=0 or angle(x)=π). +""" +function relu(z::Complex) + return 0.5 * (1 + cos(angle(z))) * z +end + +function relu_cap(z::Complex; ω₀=1.0) min_value = - ω₀ max_value = + ω₀ - return min_value + (max_value - min_value) * max(0.0, min(real(x), 1.0)) + (max_value - min_value) * max(0.0, min(imag(x), 1.0)) * im + return min_value + (max_value - min_value) * relu(z - relu(z-1)) end """ From 3e3016070ba30248e6061b45aaf0d835788d390a Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Mon, 29 Apr 2024 19:39:19 -0700 Subject: [PATCH 09/22] predict function, return multiple losses --- Project.toml | 1 + examples/double_rotation/double_rotation.jl | 16 +++-- src/SphereUDE.jl | 1 + src/train.jl | 67 +++++++++++++-------- src/utils.jl | 14 ++--- 5 files changed, 63 insertions(+), 36 deletions(-) diff --git a/Project.toml b/Project.toml index e59eec6..2edae2f 100644 --- a/Project.toml +++ b/Project.toml @@ -18,6 +18,7 @@ Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" diff --git a/examples/double_rotation/double_rotation.jl b/examples/double_rotation/double_rotation.jl index f927deb..00bca5a 100644 --- a/examples/double_rotation/double_rotation.jl +++ b/examples/double_rotation/double_rotation.jl @@ -17,10 +17,12 @@ using Random rng = Random.default_rng() Random.seed!(rng, 666) +function run() + # Total time simulation tspan = [0, 160.0] # Number of sample points -N_samples = 50 +N_samples = 100 # Times where we sample points times_samples = sort(rand(sampler(Uniform(tspan[1], tspan[2])), N_samples)) @@ -35,8 +37,8 @@ L0 = ω₀ .* [1.0, 0.0, 0.0] L1 = 0.5ω₀ .* [0.0, 1/sqrt(2), 1/sqrt(2)] # Solver tolerances -reltol = 1e-7 -abstol = 1e-7 +reltol = 1e-12 +abstol = 1e-12 function L_true(t::Float64; τ₀=τ₀, p=[L0, L1]) if t < τ₀ @@ -73,7 +75,8 @@ regs = [Regularization(order=1, power=1.0, λ=1.0, diff_mode="CS"), params = SphereParameters(tmin=tspan[1], tmax=tspan[2], reg=regs, u0=[0.0, 0.0, -1.0], ωmax=ω₀, reltol=reltol, abstol=abstol, - niter_ADAM=1000, niter_LBFGS=800) + niter_ADAM=1000, niter_LBFGS=600, + sensealg=GaussAdjoint(autojacvec=ReverseDiffVJP(true))) results = train(data, params, rng, nothing; train_initial_condition=false) @@ -82,4 +85,7 @@ results = train(data, params, rng, nothing; train_initial_condition=false) ############################################################## plot_sphere(data, results, -20., 150., saveas="examples/double_rotation/plot_sphere.pdf", title="Double rotation") -plot_L(data, results, saveas="examples/double_rotation/plot_L.pdf", title="Double rotation") \ No newline at end of file +plot_L(data, results, saveas="examples/double_rotation/plot_L.pdf", title="Double rotation") + +end +run() \ No newline at end of file diff --git a/src/SphereUDE.jl b/src/SphereUDE.jl index ebce6ae..746d615 100644 --- a/src/SphereUDE.jl +++ b/src/SphereUDE.jl @@ -13,6 +13,7 @@ using SciMLSensitivity using Optimization, OptimizationOptimisers, OptimizationOptimJL using ComponentArrays: ComponentVector using PyPlot, PyCall +using PrettyTables # Testing double-differentiation # using BatchedRoutines diff --git a/src/train.jl b/src/train.jl index aec7ee7..ba946c2 100644 --- a/src/train.jl +++ b/src/train.jl @@ -20,13 +20,16 @@ end predict_L(t, NN, θ, st) Predict value of rotation given by L given by the neural network. - -% To do: replace all the calls in U by predict_L """ function predict_L(t, NN, θ, st) return NN([t], θ, st)[1] end +""" + train() + +Training function. +""" function train(data::AD, params::AP, rng, @@ -47,8 +50,7 @@ function train(data::AD, function ude_rotation!(du, u, p, t) # Angular momentum given by network prediction - L = U([t], p, st)[1] - # L = predict_L(t, U, p, st) + L = predict_L(t, U, p, st) du .= cross(L, u) end @@ -64,13 +66,24 @@ function train(data::AD, tspan=(min(T[1], params.tmin), max(T[end], params.tmax)), p = β.θ) end - Array(solve(_prob, params.solver, saveat=T, + sol = solve(_prob, params.solver, saveat=T, abstol=params.abstol, reltol=params.reltol, - sensealg=params.sensealg)) + sensealg=params.sensealg) + return Array(sol), sol.retcode end function loss(β::ComponentVector) - u_ = predict(β) + u_, retcode = predict(β) + + # If numerical integration fails or bad choice of parameter, return infinity + if retcode != :Success + @warn "[SphereUDE] Numerical solver not converging. This can be causes by numerical innestabilities around a bad choice of parameter." + return Inf + end + + # Record the value of each individual loss to the total loss function for hyperparameter selection. + loss_dict = Dict() + # Empirical error if isnothing(data.kappas) l_emp = 3.0 * mean(abs2.(u_ .- data.directions)) @@ -80,22 +93,26 @@ function train(data::AD, l_emp = mean(data.kappas .* abs2.(u_ .- data.directions), dims=1) # l_emp = norm(data.kappas)^2 - 3.0 * mean(data.kappas .* u_ .* data.directions) end + loss_dict["Empirical"] = l_emp + # Regularization l_reg = 0.0 if !isnothing(params.reg) # for (order, power, λ) in params.reg for reg in params.reg - l_reg += regularization(β.θ, reg) + reg₀ = regularization(β.θ, reg) + l_reg += reg₀ + loss_dict["Regularization (order=$(reg.order), power=$(reg.power)"] = reg₀ end end - return l_emp + l_reg + return l_emp + l_reg, loss_dict end function regularization(θ::ComponentVector, reg::AbstractRegularization; n_nodes=100) l_ = 0.0 if reg.order==0 - l_ += quadrature(t -> norm(U([t], θ, st)[1])^reg.power, params.tmin, params.tmax, n_nodes) + l_ += quadrature(t -> norm(predict_L(t, U, θ, st))^reg.power, params.tmin, params.tmax, n_nodes) elseif reg.order==1 if reg.diff_mode=="AD" throw("Method not working well.") @@ -104,19 +121,13 @@ function train(data::AD, # Test this with the new implementation in Lux.jl: # https://lux.csail.mit.edu/stable/manual/nested_autodiff - - for t in times_reg - # Try ReverseDiff - grad = Zygote.jacobian(first ∘ U, [t], θ, st)[1] - l_ += norm(grad)^reg.power - end elseif reg.diff_mode=="FD" # Finite differences ϵ = 0.1 * (params.tmax - params.tmin) / n_nodes - l_ += quadrature(t -> norm(central_fdm(τ -> (first ∘ U)([τ], θ, st), t, ϵ=ϵ))^reg.power, params.tmin, params.tmax, n_nodes) + l_ += quadrature(t -> norm(central_fdm(τ -> predict_L(τ, U, θ, st), t, ϵ=ϵ))^reg.power, params.tmin, params.tmax, n_nodes) elseif reg.diff_mode=="CS" # Complex step differentiation - l_ += quadrature(t -> norm(complex_step_differentiation(τ -> (first ∘ U)([τ], θ, st), t))^reg.power, params.tmin, params.tmax, n_nodes) + l_ += quadrature(t -> norm(complex_step_differentiation(τ -> predict_L(τ, U, θ, st), t))^reg.power, params.tmin, params.tmax, n_nodes) else throw("Method not implemented.") end @@ -129,7 +140,7 @@ function train(data::AD, losses = Float64[] callback = function (p, l) push!(losses, l) - if length(losses) % 10 == 0 + if length(losses) % 50 == 0 println("Current loss after $(length(losses)) iterations: $(losses[end])") end if train_initial_condition @@ -139,15 +150,19 @@ function train(data::AD, end adtype = Optimization.AutoZygote() - optf = Optimization.OptimizationFunction((x, β) -> loss(x), adtype) + optf = Optimization.OptimizationFunction((x, β) -> (first ∘ loss)(x), adtype) optprob = Optimization.OptimizationProblem(optf, β) res1 = Optimization.solve(optprob, ADAM(0.002), callback=callback, maxiters=params.niter_ADAM) println("Training loss after $(length(losses)) iterations: $(losses[end])") - optprob2 = Optimization.OptimizationProblem(optf, res1.u) - res2 = Optimization.solve(optprob2, Optim.LBFGS(), callback=callback, maxiters=params.niter_LBFGS) - println("Final training loss after $(length(losses)) iterations: $(losses[end])") + if params.niter_LBFGS > 0 + optprob2 = Optimization.OptimizationProblem(optf, res1.u) + res2 = Optimization.solve(optprob2, Optim.LBFGS(), callback=callback, maxiters=params.niter_LBFGS) + println("Final training loss after $(length(losses)) iterations: $(losses[end])") + else + res2 = res1 + end # Optimized NN parameters β_trained = res2.u @@ -162,7 +177,11 @@ function train(data::AD, # Final Fit fit_times = collect(range(params.tmin,params.tmax, length=200)) - fit_directions = predict(β_trained, T=fit_times) + fit_directions, _ = predict(β_trained, T=fit_times) + + # Recover final balance between different terms involved in the loss function to assess hyperparameter selection. + _, loss_dict = loss(β_trained) + pretty_table(loss_dict, sortkeys=true, header=["Loss term", "Value"]) return Results(θ=θ_trained, u0=u0_trained, U=U, st=st, fit_times=fit_times, fit_directions=fit_directions) diff --git a/src/utils.jl b/src/utils.jl index b8cb4f5..b7dfb54 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -166,12 +166,12 @@ function raise_warnings(data::SphereData, params::SphereParameters) if length(unique(data.times)) < length(data.times) @warn "[SphereUDE] Timeseries includes duplicated times. \n This can produce unexpected errors." end - if !isnothing(params.reg) - for reg in params.reg - if reg.diff_mode=="CS" - @warn "[SphereUDE] Complex-step differentiation inside the loss function \n This just work for cases where the activation function of the neural network admits complex numbers \n Change predefined activation functions to be complex numbers." - end - end - end + # if !isnothing(params.reg) + # for reg in params.reg + # if reg.diff_mode=="CS" + # @warn "[SphereUDE] Complex-step differentiation inside the loss function \n This just work for cases where the activation function of the neural network admits complex numbers \n Change predefined activation functions to be complex numbers." + # end + # end + # end nothing end \ No newline at end of file From 10ae8cd995badca8889a2d20685cffbec2d90069 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Mon, 29 Apr 2024 19:43:13 -0700 Subject: [PATCH 10/22] Co-authored-by: Jordi Bolibar --- examples/double_rotation/double_rotation.jl | 1 - src/utils.jl | 5 ----- 2 files changed, 6 deletions(-) diff --git a/examples/double_rotation/double_rotation.jl b/examples/double_rotation/double_rotation.jl index 6510fa9..00bca5a 100644 --- a/examples/double_rotation/double_rotation.jl +++ b/examples/double_rotation/double_rotation.jl @@ -80,7 +80,6 @@ params = SphereParameters(tmin=tspan[1], tmax=tspan[2], results = train(data, params, rng, nothing; train_initial_condition=false) - ############################################################## ###################### PyCall Plots ######################### ############################################################## diff --git a/src/utils.jl b/src/utils.jl index 4c17b43..b7dfb54 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -35,10 +35,6 @@ function sigmoid(z::Complex) # end end -function sigmoid(x::Complex) - return 1 / ( 1.0 + exp(-x) ) -end - """ relu_cap(x; ω₀=1.0) """ @@ -64,7 +60,6 @@ function relu_cap(z::Complex; ω₀=1.0) min_value = - ω₀ max_value = + ω₀ return min_value + (max_value - min_value) * relu(z - relu(z-1)) - end """ From 3d5269180a9388efb2d01a342e548a758ceba7f8 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Thu, 2 May 2024 19:44:57 -0700 Subject: [PATCH 11/22] Multiple shooting working once sensealg specified --- src/train.jl | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/src/train.jl b/src/train.jl index e819b28..863d869 100644 --- a/src/train.jl +++ b/src/train.jl @@ -85,8 +85,8 @@ function train(data::AD, # Empirical error if isnothing(data.kappas) - l_emp = 3.0 * mean(abs2.(u_ .- data.directions)) # The 3 is needed since the mean is computen on a 3xN matrix + l_emp = 3.0 * mean(abs2.(u_ .- data.directions)) # l_emp = 1 - 3.0 * mean(u_ .* data.directions) else l_emp = mean(data.kappas .* abs2.(u_ .- data.directions), dims=1) @@ -112,13 +112,10 @@ function train(data::AD, # Empirical error l_emp = 3.0 * mean(abs2.(pred .- data)) - # The 3 is needed since the mean is computen on a 3xN matrix - # l_emp = 1 - 3.0 * mean(u_ .* data.directions) # Regularization l_reg = 0.0 if !isnothing(params.reg) - # for (order, power, λ) in params.reg for reg in params.reg reg₀ = regularization(β.θ, reg) l_reg += reg₀ @@ -129,11 +126,19 @@ function train(data::AD, end # Define parameters for Multiple Shooting - group_size = 50 + group_size = 10 continuity_term = 100 - ps = ComponentArray(θ) - pd, pax = getdata(ps), getaxes(ps) + ps = ComponentArray(θ) # are these necesary? + pd, pax = getdata(ps), getaxes(ps) + + function continuity_loss(u_pred, u_initial) + if !isapprox(norm(u_initial), 1.0, atol=1e-6) || !isapprox(norm(u_pred), 1.0, atol=1e-6) + @warn "Directions during multiple shooting are not in the sphere. Small deviations from unit norm observed:" + @show norm(u_initial), norm(u_pred) + end + return sum(abs2, u_pred - u_initial) + end function loss_multiple_shooting(β::ComponentVector) @@ -149,11 +154,11 @@ function train(data::AD, p = β.θ) end - return multiple_shoot(β.θ, data.directions, data.times, _prob, loss_function, Tsit5(), - group_size; continuity_term) + return multiple_shoot(β.θ, data.directions, data.times, _prob, loss_function, continuity_loss, params.solver, + group_size; continuity_term, sensealg=params.sensealg) end - function regularization(θ::ComponentVector, reg::AG; n_nodes=100) where{AG <: AbstractRegularization} + function regularization(θ::ComponentVector, reg::AG; n_nodes=100) where {AG <: AbstractRegularization} l_ = 0.0 if reg.order==0 From 9fa33541fb8703f74ce186294ac922380ab6cc90 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Fri, 19 Jul 2024 21:49:28 -0300 Subject: [PATCH 12/22] Double rotation example with small changes in src --- examples/double_rotation/double_rotation.jl | 73 +++++++++++++-------- src/SphereUDE.jl | 1 + src/plot.jl | 14 +++- src/setup/config.jl | 5 +- src/train.jl | 39 +++++++---- src/types.jl | 6 +- src/utils.jl | 25 ++++++- 7 files changed, 117 insertions(+), 46 deletions(-) diff --git a/examples/double_rotation/double_rotation.jl b/examples/double_rotation/double_rotation.jl index 27682f1..5587565 100644 --- a/examples/double_rotation/double_rotation.jl +++ b/examples/double_rotation/double_rotation.jl @@ -17,7 +17,7 @@ using Random rng = Random.default_rng() Random.seed!(rng, 666) -function run() +function run(;kappa=50., regs=regs, title="plot") # Total time simulation tspan = [0, 160.0] @@ -31,10 +31,10 @@ times_samples = sort(rand(sampler(Uniform(tspan[1], tspan[2])), N_samples)) # Angular velocity ω₀ = Δω₀ * π / 180.0 # Change point -τ₀ = 65.0 +τ₀ = 70.0 # Angular momentum L0 = ω₀ .* [1.0, 0.0, 0.0] -L1 = 0.5ω₀ .* [0.0, 1/sqrt(2), 1/sqrt(2)] +L1 = 0.6ω₀ .* [0.0, 1/sqrt(2), 1/sqrt(2)] # Solver tolerances reltol = 1e-12 @@ -49,37 +49,30 @@ function L_true(t::Float64; τ₀=τ₀, p=[L0, L1]) end function true_rotation!(du, u, p, t) - L = L_true(t; τ₀=τ₀, p=p) + L = L_true(t; τ₀ = τ₀, p = p) du .= cross(L, u) end prob = ODEProblem(true_rotation!, [0.0, 0.0, -1.0], tspan, [L0, L1]) -true_sol = solve(prob, Tsit5(), reltol=reltol, abstol=abstol, saveat=times_samples) +true_sol = solve(prob, Tsit5(), reltol = reltol, abstol = abstol, saveat = times_samples) # Add Fisher noise to true solution X_noiseless = Array(true_sol) -X_true = X_noiseless + FisherNoise(kappa=50.) +X_true = X_noiseless + FisherNoise(kappa=kappa) ############################################################## ####################### Training ########################### ############################################################## -data = SphereData(times=times_samples, directions=X_true, kappas=nothing, L=L_true) +data = SphereData(times=times_samples, directions=X_true, kappas=nothing, L=L_true) -regs = [Regularization(order=0, power=2.0, λ=0.1, diff_mode=nothing)] -# regs = [Regularization(order=1, power=1.0, λ=1.0, diff_mode="CS"), -# Regularization(order=0, power=2.0, λ=0.1, diff_mode=nothing)] -# regs = [Regularization(order=0, power=2.0, λ=0.1, diff_mode=nothing)] - # Regularization(order=1, power=1.1, λ=0.01, diff_mode="CS")] -# regs = nothing - -params = SphereParameters(tmin=tspan[1], tmax=tspan[2], - reg=regs, - train_initial_condition=false, - multiple_shooting=true, - u0=[0.0, 0.0, -1.0], ωmax=ω₀, reltol=reltol, abstol=abstol, - niter_ADAM=1000, niter_LBFGS=600, - sensealg=GaussAdjoint(autojacvec=ReverseDiffVJP(true))) +params = SphereParameters(tmin = tspan[1], tmax = tspan[2], + reg = regs, + train_initial_condition = false, + multiple_shooting = false, + u0 = [0.0, 0.0, -1.0], ωmax = ω₀, reltol = reltol, abstol = abstol, + niter_ADAM = 2000, niter_LBFGS = 1000, + sensealg = GaussAdjoint(autojacvec = ReverseDiffVJP(true))) results = train(data, params, rng, nothing) @@ -87,8 +80,36 @@ results = train(data, params, rng, nothing) ###################### PyCall Plots ######################### ############################################################## -plot_sphere(data, results, -20., 150., saveas="examples/double_rotation/plot_sphere.pdf", title="Double rotation") -plot_L(data, results, saveas="examples/double_rotation/plot_L.pdf", title="Double rotation") - -end -run() \ No newline at end of file +plot_sphere(data, results, -20., 125., saveas="examples/double_rotation/" * title * "_sphere.pdf", title="Double rotation") # , matplotlib_rcParams=Dict("font.size"=> 50)) +plot_L(data, results, saveas="examples/double_rotation/" * title * "_L.pdf", title="Double rotation") + +end # run + +# Run different experiments + +λ₀ = 0.1 +λ₁ = 0.1 +# λ₀ = 0.1 +# λ₁ = 0.01 +run(; kappa = 50., + regs = [Regularization(order=1, power=1.0, λ=λ₁, diff_mode="CS"), + Regularization(order=0, power=2.0, λ=λ₀, diff_mode=nothing)], + title = "plot_50_lambda$(λ₁)") + +λ₀ = 0.1 +λ₁ = 0.1 +# λ₀ = 0.1 +# λ₁ = 0.01 +run(; kappa = 200., + regs = [Regularization(order=1, power=1.0, λ=λ₁, diff_mode="CS"), + Regularization(order=0, power=2.0, λ=λ₀)], + title = "plot_200_lambda$(λ₁)") + +λ₀ = 0.1 +λ₁ = 0.1 +# λ₀ = 0.1 +# λ₁ = 0.01 +run(; kappa = 1000., + regs = [Regularization(order=1, power=1.0, λ=λ₁, diff_mode="CS"), + Regularization(order=0, power=2.0, λ=λ₀)], + title = "plot_1000_lambda$(λ₁)") \ No newline at end of file diff --git a/src/SphereUDE.jl b/src/SphereUDE.jl index be87a6b..7a49453 100644 --- a/src/SphereUDE.jl +++ b/src/SphereUDE.jl @@ -28,6 +28,7 @@ include("train.jl") include("plot.jl") # Python libraries +const mpl_base::PyObject = isdefined(SphereUDE, :mpl_base) ? SphereUDE.mpl_base : PyNULL() const mpl_colors::PyObject = isdefined(SphereUDE, :mpl_colors) ? SphereUDE.mpl_colors : PyNULL() const mpl_colormap::PyObject = isdefined(SphereUDE, :mpl_colormap) ? SphereUDE.mpl_colormap : PyNULL() const sns::PyObject = isdefined(SphereUDE, :sns) ? SphereUDE.sns : PyNULL() diff --git a/src/plot.jl b/src/plot.jl index 7deb758..c33df4e 100644 --- a/src/plot.jl +++ b/src/plot.jl @@ -23,12 +23,22 @@ function plot_sphere(# ax::PyCall.PyObject, central_latitude::Float64, central_longitude::Float64; saveas::Union{String, Nothing}, - title::String) + title::String, + matplotlib_rcParams::Union{Dict, Nothing} = nothing) # cmap = mpl_colormap.get_cmap("viridis") plt.figure(figsize=(10,10)) ax = plt.axes(projection=ccrs.Orthographic(central_latitude=central_latitude, central_longitude=central_longitude)) + + # Set default plot parameters. + # See https://matplotlib.org/stable/users/explain/customizing.html for customizable optionsz + if !isnothing(matplotlib_rcParams) + for (key, value) in matplotlib_rcParams + @warn "Setting Matplotlib parameters with rcParams currently not working. See following GitHub issue: https://github.com/JuliaPy/PyPlot.jl/issues/525" + mpl_base.rcParams[key] = value + end + end # ax.coastlines() ax.gridlines() @@ -51,7 +61,7 @@ function plot_sphere(# ax::PyCall.PyObject, linewidth=2, color="black",#cmap(norm(results.fit_times[i])), transform = ccrs.Geodetic()) end - plt.title(title) + plt.title(title, fontsize=20) if !isnothing(saveas) plt.savefig(saveas, format="pdf") end diff --git a/src/setup/config.jl b/src/setup/config.jl index 7472e96..824b52d 100644 --- a/src/setup/config.jl +++ b/src/setup/config.jl @@ -1,8 +1,9 @@ -export mpl_colormap, mpl_colormap, sns, ccrs, feature +export mpl_base, mpl_colormap, mpl_colormap, sns, ccrs, feature function __init__() try + copy!(mpl_base, pyimport("matplotlib")) copy!(mpl_colors, pyimport("matplotlib.colors")) copy!(mpl_colormap, pyimport("matplotlib.cm")) copy!(sns, pyimport("seaborn")) @@ -10,7 +11,7 @@ function __init__() copy!(feature, pyimport("cartopy.feature")) catch e @warn "It looks like you have not installed and/or activated the virtual Python environment. \n - Please follow the guidelines in: https://github.com/facusapienza21/SphereUDE.jl#readme" + Please follow the guidelines in: https://github.com/ODINN-SciML/SphereUDE.jl" @warn exception=(e, catch_backtrace()) end diff --git a/src/train.jl b/src/train.jl index 863d869..3d48eac 100644 --- a/src/train.jl +++ b/src/train.jl @@ -1,17 +1,29 @@ export train -# For L1 regularization relu_cap works better, but for L2 I think is better to include sigmoid function get_NN(params, rng, θ_trained) # Define neural network - U = Lux.Chain( - Lux.Dense(1, 5, sigmoid), # explore discontinuity function for activation - Lux.Dense(5, 10, sigmoid), - Lux.Dense(10, 5, sigmoid), - # Lux.Dense(1, 5, relu_cap), - # Lux.Dense(5, 10, relu_cap), - # Lux.Dense(10, 5, relu_cap), - Lux.Dense(5, 3, x->sigmoid_cap(x; ω₀=params.ωmax)) - ) + + # For L1 regularization relu_cap works better, but for L2 I think is better to include sigmoid + if isL1reg(params.reg) + @warn "[SphereUDE] Using ReLU activation functions for neural network due to L1 regularization." + U = Lux.Chain( + Lux.Dense(1, 5, sigmoid), + Lux.Dense(5, 10, sigmoid), + Lux.Dense(10, 5, sigmoid), + # Lux.Dense(1, 5, relu_cap), + # Lux.Dense(5, 10, relu_cap), + # Lux.Dense(10, 5, relu_cap), + Lux.Dense(5, 3, x -> sigmoid_cap(x; ω₀=params.ωmax)) + # Lux.Dense(5, 3, x -> relu_cap(x; ω₀=params.ωmax)) + ) + else + U = Lux.Chain( + Lux.Dense(1, 5, sigmoid), + Lux.Dense(5, 10, sigmoid), + Lux.Dense(10, 5, sigmoid), + Lux.Dense(5, 3, x -> sigmoid_cap(x; ω₀=params.ωmax)) + ) + end θ, st = Lux.setup(rng, U) return U, θ, st end @@ -67,7 +79,8 @@ function train(data::AD, end sol = solve(_prob, params.solver, saveat=T, abstol=params.abstol, reltol=params.reltol, - sensealg=params.sensealg) + sensealg=params.sensealg, + dtmin=1e-4 * (params.tmax - params.tmin), force_dtmin=true) # Force minimum step in case L(t) changes drastically due to bad behaviour of neural network return Array(sol), sol.retcode end @@ -101,7 +114,7 @@ function train(data::AD, for reg in params.reg reg₀ = regularization(β.θ, reg) l_reg += reg₀ - loss_dict["Regularization (order=$(reg.order), power=$(reg.power)"] = reg₀ + loss_dict["Regularization (order=$(reg.order), power=$(reg.power))"] = reg₀ end end return l_emp + l_reg, loss_dict @@ -211,7 +224,7 @@ function train(data::AD, if params.niter_LBFGS > 0 optprob2 = Optimization.OptimizationProblem(optf, res1.u) - res2 = Optimization.solve(optprob2, Optim.LBFGS(), callback=callback, maxiters=params.niter_LBFGS) + res2 = Optimization.solve(optprob2, Optim.LBFGS(), callback=callback, maxiters=params.niter_LBFGS) #, reltol=1e-6) println("Final training loss after $(length(losses)) iterations: $(losses[end])") else res2 = res1 diff --git a/src/types.jl b/src/types.jl index 8e3783a..315720e 100644 --- a/src/types.jl +++ b/src/types.jl @@ -56,5 +56,9 @@ Regularization information order::I # Order of derivative power::F # Power of the Euclidean norm λ::F # Regularization hyperparameter - diff_mode::Union{Nothing, String} # AD differentiation mode used in regulatization + # AD differentiation mode used in regulatization + diff_mode::Union{Nothing, String} = nothing + + # Include this in the constructor + # @assert (order == 0) || (!isnothing(diff_mode)) "Diffentiation methods needs to be provided for regularization with order larger than zero." end \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl index ac049a2..c632e19 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -4,7 +4,7 @@ export cart2sph, sph2cart export AbstractNoise, FisherNoise export quadrature, central_fdm, complex_step_differentiation export raise_warnings - +export isL1reg """ sigmoid_cap(x; ω₀=1.0) @@ -41,9 +41,12 @@ end function relu_cap(x; ω₀=1.0) min_value = - ω₀ max_value = + ω₀ - return min_value + (max_value - min_value) * max(0.0, min(x, 1.0)) + return relu_cap(x, min_value, max_value) end +function relu_cap(x, min_value::Float64, max_value::Float64) + return min_value + (max_value - min_value) * max(0.0, min(x, 1.0)) +end """ relu(x::Complex) @@ -59,6 +62,11 @@ end function relu_cap(z::Complex; ω₀=1.0) min_value = - ω₀ max_value = + ω₀ + return relu_cap(z, min_value, max_value) + # return min_value + (max_value - min_value) * relu(z - relu(z-1)) +end + +function relu_cap(z::Complex, min_value::Float64, max_value::Float64) return min_value + (max_value - min_value) * relu(z - relu(z-1)) end @@ -174,4 +182,17 @@ function raise_warnings(data::SphereData, params::SphereParameters) # end # end nothing +end + +""" + +Function to check for the presence of L1 regularization in the loss function. +""" +function isL1reg(regs::Vector{R}) where {R <: AbstractRegularization} + for reg in regs + if reg.power == 1 + return true + end + end + return false end \ No newline at end of file From e2116c39128915ffe3a4c31755a7b98ebcbc2989 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Sat, 21 Sep 2024 15:01:51 -0400 Subject: [PATCH 13/22] feat: update to support Lux 1.0 --- Project.toml | 13 ++++++++----- src/SphereUDE.jl | 4 +--- src/types.jl | 4 ++-- 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index 3db8f61..6241b8a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SphereUDE" uuid = "d7416ba7-148a-4110-b27d-9087fcebab2d" authors = ["Facundo Sapienza ", "Jordi Bolibar "] -version = "0.1.1" +version = "0.1.2" [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" @@ -19,32 +19,35 @@ Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1" OptimizationPolyalgorithms = "500b13db-7e66-49ce-bda4-eed966be6282" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] BenchmarkTools = "1" ComponentArrays = "0.15" +DiffEqFlux = "4" Distributions = "0.25" Infiltrator = "1.2" +Lux = "1" Optimization = "3.12" OptimizationOptimJL = "0.1.5" OptimizationOptimisers = "0.1.2" -OrdinaryDiffEq = "5, 6" +OrdinaryDiffEqCore = "1.6.0" +OrdinaryDiffEqTsit5 = "1.1.0" PyCall = "1.9" PyPlot = "2.11" Revise = "3.1" SciMLSensitivity = "7.20" Statistics = "1" Zygote = "0.6" -julia = "1.7" +julia = "1.10" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/SphereUDE.jl b/src/SphereUDE.jl index be87a6b..d341fc7 100644 --- a/src/SphereUDE.jl +++ b/src/SphereUDE.jl @@ -1,15 +1,13 @@ __precompile__() module SphereUDE -# types -using Base: @kwdef # utils # training using LinearAlgebra, Statistics, Distributions using FastGaussQuadrature using Lux, Zygote, DiffEqFlux using ChainRules: @ignore_derivatives -using OrdinaryDiffEq +using OrdinaryDiffEqCore, OrdinaryDiffEqTsit5 using SciMLSensitivity using Optimization, OptimizationOptimisers, OptimizationOptimJL, OptimizationPolyalgorithms using ComponentArrays diff --git a/src/types.jl b/src/types.jl index 8e3783a..fe43787 100644 --- a/src/types.jl +++ b/src/types.jl @@ -23,7 +23,7 @@ Training parameters niter_LBFGS::I reltol::F abstol::F - solver::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = Tsit5() + solver::OrdinaryDiffEqCore.OrdinaryDiffEqAlgorithm = Tsit5() sensealg::SciMLBase.AbstractAdjointSensitivityAlgorithm = QuadratureAdjoint(autojacvec=ReverseDiffVJP(true)) end @@ -43,7 +43,7 @@ Final results @kwdef struct Results{F <: AbstractFloat} <: AbstractResult θ::ComponentVector u0::Vector{F} - U::Lux.Chain + U::Lux.AbstractLuxLayer st::NamedTuple fit_times::Vector{F} fit_directions::Matrix{F} From f4eee7a53bc6495e17e5e48718531968a1d35c9b Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Thu, 26 Sep 2024 15:05:05 -0700 Subject: [PATCH 14/22] Example with double rotation working with non-updated Lux --- Project.toml | 3 ++- README.md | 8 ++++++-- examples/double_rotation/double_rotation.jl | 19 ++++++++----------- 3 files changed, 16 insertions(+), 14 deletions(-) diff --git a/Project.toml b/Project.toml index 3db8f61..48bc457 100644 --- a/Project.toml +++ b/Project.toml @@ -34,6 +34,7 @@ BenchmarkTools = "1" ComponentArrays = "0.15" Distributions = "0.25" Infiltrator = "1.2" +Lux = "<0.5.49" Optimization = "3.12" OptimizationOptimJL = "0.1.5" OptimizationOptimisers = "0.1.2" @@ -50,4 +51,4 @@ julia = "1.7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test"] +test = ["Test"] \ No newline at end of file diff --git a/README.md b/README.md index 9ee6936..481ae68 100644 --- a/README.md +++ b/README.md @@ -64,7 +64,11 @@ To make plots using Matplotlib, Cartopy, and PMagPy, we install both [PyCall.jl] - Create a Python conda environment, based on [this conda environment file](https://raw.githubusercontent.com/facusapienza21/SphereUDE.jl/main/environment.yml), with all the required packages using `conda env create -f environment.yml`. - Inside the Julia REPL, install both `PyCall.jl` and `PyPlot.jl` with `] add PyCall, Pyplot`. -- Specify the Python path of the new environment with `ENV["PYTHON"] = ...`, where you should complete the path of the Python installation that shows when you do `conda activate SphereUDE`, `which python`. -- Inside the Julia REPL, execute `Pkg.build("PyCall")` to re-build PyCall with the new Python path. +- Specify the Python path of the new environment with `ENV["PYTHON"] = ...`, where you should complete the path of the Python installation that shows when you do `conda activate SphereUDE`, `which python`. Inside the Julia REPL, execute `Pkg.build("PyCall")` to re-build PyCall with the new Python path: +``` +julia> ENV["PYTHON"] = read(`which python`, String)[1:end-1] # trim backspace +julia> import Pkg; Pkg.build("PyCall") +julia> exit() +``` You are ready to use Python from your Julia session! diff --git a/examples/double_rotation/double_rotation.jl b/examples/double_rotation/double_rotation.jl index 5587565..597d53a 100644 --- a/examples/double_rotation/double_rotation.jl +++ b/examples/double_rotation/double_rotation.jl @@ -88,28 +88,25 @@ end # run # Run different experiments λ₀ = 0.1 -λ₁ = 0.1 -# λ₀ = 0.1 -# λ₁ = 0.01 +λ₁ = 0.001 + run(; kappa = 50., - regs = [Regularization(order=1, power=1.0, λ=λ₁, diff_mode="CS"), + regs = [Regularization(order=1, power=1.0, λ=λ₁, diff_mode="FD"), Regularization(order=0, power=2.0, λ=λ₀, diff_mode=nothing)], - title = "plot_50_lambda$(λ₁)") + title = "plots/plot_50_lambda$(λ₁)") λ₀ = 0.1 λ₁ = 0.1 -# λ₀ = 0.1 -# λ₁ = 0.01 + run(; kappa = 200., regs = [Regularization(order=1, power=1.0, λ=λ₁, diff_mode="CS"), Regularization(order=0, power=2.0, λ=λ₀)], - title = "plot_200_lambda$(λ₁)") + title = "plots/plot_200_lambda$(λ₁)") λ₀ = 0.1 λ₁ = 0.1 -# λ₀ = 0.1 -# λ₁ = 0.01 + run(; kappa = 1000., regs = [Regularization(order=1, power=1.0, λ=λ₁, diff_mode="CS"), Regularization(order=0, power=2.0, λ=λ₀)], - title = "plot_1000_lambda$(λ₁)") \ No newline at end of file + title = "plots/plot_1000_lambda$(λ₁)") \ No newline at end of file From 8a01ca5a16da67448bf173718d971766b0099075 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Thu, 26 Sep 2024 16:18:13 -0700 Subject: [PATCH 15/22] Integration test of inversion --- test/rotation.jl | 63 ++++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 6 ++++- 2 files changed, 68 insertions(+), 1 deletion(-) create mode 100644 test/rotation.jl diff --git a/test/rotation.jl b/test/rotation.jl new file mode 100644 index 0000000..5e78b1c --- /dev/null +++ b/test/rotation.jl @@ -0,0 +1,63 @@ +using LinearAlgebra, Statistics, Distributions +using OrdinaryDiffEq +using SciMLSensitivity +using Optimization, OptimizationOptimisers, OptimizationOptimJL + +using Random +rng = Random.default_rng() +Random.seed!(rng, 666) + +############################################################## +############### Simulation of Simple Rotation ############### +############################################################## + +function test_single_rotation() + + # Total time simulation + tspan = [0, 160.0] + # Number of sample points + N_samples = 10 + # Times where we sample points + times_samples = sort(rand(sampler(Uniform(tspan[1], tspan[2])), N_samples)) + + # Expected maximum angular deviation in one unit of time (degrees) + Δω₀ = 1.0 + # Angular velocity + ω₀ = Δω₀ * π / 180.0 + + # Create simple example + X = zeros(3, N_samples) + X[3, :] .= 1 + X[1, :] = LinRange(0,1,N_samples) + X = X ./ norm.(eachcol(X))' + + ############################################################## + ####################### Training ########################### + ############################################################## + + data = SphereData(times=times_samples, directions=X, kappas=nothing, L=nothing) + + regs = [Regularization(order=1, power=1.0, λ=0.1, diff_mode="FD"), + Regularization(order=0, power=2.0, λ=0.001, diff_mode=nothing)] + + params = SphereParameters(tmin = tspan[1], tmax = tspan[2], + reg = regs, + train_initial_condition = false, + multiple_shooting = false, + u0 = [0.0, 0.0, -1.0], ωmax = ω₀, reltol = 1e-12, abstol = 1e-12, + niter_ADAM = 20, niter_LBFGS = 10, + sensealg = GaussAdjoint(autojacvec = ReverseDiffVJP(true))) + + results = train(data, params, rng, nothing) + + @test true + +end + +############################################################## +###################### PyCall Plots ######################### +############################################################## + +# plot_sphere(data, results, -20., 125., saveas="examples/double_rotation/" * title * "_sphere.pdf", title="Double rotation") # , matplotlib_rcParams=Dict("font.size"=> 50)) +# plot_L(data, results, saveas="examples/double_rotation/" * title * "_L.pdf", title="Double rotation") + diff --git a/test/runtests.jl b/test/runtests.jl index f33bf74..db621eb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,7 +4,7 @@ using Lux include("constructors.jl") include("utils.jl") - +include("rotation.jl") @testset "Constructors" begin test_reg_constructor() @@ -16,3 +16,7 @@ end test_complex_activation() test_integration() end + +@testset "Inversion" begin + test_single_rotation() +end \ No newline at end of file From 8ee17ed3454c8139d240a7b6640b6ea32a682dd2 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Thu, 26 Sep 2024 16:37:07 -0700 Subject: [PATCH 16/22] Added Random as test dependency --- Project.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 48bc457..fe0c86b 100644 --- a/Project.toml +++ b/Project.toml @@ -49,6 +49,7 @@ julia = "1.7" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [targets] -test = ["Test"] \ No newline at end of file +test = ["Test", "Random"] From aad95b41ebe07bf01cda2e001de59aab8d3e0804 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Fri, 27 Sep 2024 14:15:20 -0700 Subject: [PATCH 17/22] Fix Lux version Co-authored-by: Avik Pal --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ddc6c1f..4b3f55d 100644 --- a/Project.toml +++ b/Project.toml @@ -35,7 +35,7 @@ ComponentArrays = "0.15" DiffEqFlux = "4" Distributions = "0.25" Infiltrator = "1.2" -Lux = "1" +Lux = "1.0" Optimization = "3.12" OptimizationOptimJL = "0.1.5" OptimizationOptimisers = "0.1.2" From f7c85f4035f18d1e166c8512bedf12089bb83929 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Mon, 14 Oct 2024 16:00:46 -0700 Subject: [PATCH 18/22] Reorganization of loss function in different module --- .../Torsvik-etal-2012_dataset.csv | 595 ++++++++++++++++++ examples/Torsvik_2012/plots/plot_results.jl | 140 +++++ src/SphereUDE.jl | 1 + src/losses.jl | 214 +++++++ src/train.jl | 262 +------- src/types.jl | 2 + src/utils.jl | 43 +- test/rotation.jl | 2 +- 8 files changed, 1018 insertions(+), 241 deletions(-) create mode 100644 examples/Torsvik_2012/Torsvik-etal-2012_dataset.csv create mode 100644 examples/Torsvik_2012/plots/plot_results.jl create mode 100644 src/losses.jl diff --git a/examples/Torsvik_2012/Torsvik-etal-2012_dataset.csv b/examples/Torsvik_2012/Torsvik-etal-2012_dataset.csv new file mode 100644 index 0000000..18ff110 --- /dev/null +++ b/examples/Torsvik_2012/Torsvik-etal-2012_dataset.csv @@ -0,0 +1,595 @@ +Q,a95,Plate,Plate_code,Lat,Lon,CLat,CLon,RLat,RLon,Eplat,Eplong,Epang,Age +5,4.8,north_america,101,-86.3,5.7,,,,,,,,0.5 +5,9.1,north_america,101,-86.4,8.4,,,-86.4,9.2,79.2,23,0.2,0.8 +5,7.1,north_america,101,-85.3,265.9,,,-85.3,-94.1,79.9,22.7,0.3,1 +5,8.7,north_america,101,-88.4,225.5,,,-88.4,-135.7,79.9,22.7,0.3,1 +5,4.3,north_america,101,-88.9,285,,,-88.9,-75.1,79.9,22.7,0.3,1 +3,4.9,north_america,101,-77.9,303.7,-80,3.5,-80,4.2,80.8,22.8,0.4,1.5 +3,7,north_america,101,-85.7,122.7,-81.2,73.9,-81.1,74.5,80.8,22.8,0.4,1.5 +5,5.7,north_america,101,-84.5,241,,,-84.5,-119.1,80.8,22.8,0.4,1.5 +5,6.2,north_america,101,-89.5,214.8,,,-89.5,-157.9,80.8,22.9,0.8,3 +5,4.1,north_america,101,-84.6,332.3,,,-84.8,-25.6,80.9,22.8,1,4 +5,6.7,north_america,101,-87.9,275.9,,,-88,-84.4,80.9,22.8,1,4 +5,5,north_america,101,-85.5,197.4,,,-85.5,-163.9,80.9,22.8,1.3,5 +5,5,north_america,101,-85.3,262.1,,,-85.5,-97.9,80.9,22.8,1.8,7 +5,6.9,north_america,101,-85.9,357.7,,,-86,3.9,80.9,22.8,2,8 +5,7.8,north_america,101,-82.9,309.1,,,-83.3,-47.3,80.9,22.9,2.6,10 +3,12.9,north_america,101,-81.1,225.3,-83.1,93.2,-82.4,98.8,80.9,23.2,4.2,15.5 +5,3.7,north_america,101,-82.9,346.2,,,-83.5,-1.3,80.3,25.4,5.8,21 +5,9.8,north_america,101,-82.7,301.7,,,-83.7,-51,80.2,26,6,21.7 +4,6.7,north_america,101,-87.1,189.5,,,-86.7,177.4,80,26.6,6.2,22.5 +3,8.4,north_america,101,-76.4,30.3,,,-76.2,41,79.9,26.9,6.3,23 +4,5.2,north_america,101,-80.9,331.2,,,-81.8,-17.8,79.8,27.3,6.4,23.5 +4,6.8,north_america,101,-79.7,342.8,,,-80.5,-4.3,79.3,27.5,7.1,26 +4,5.4,north_america,101,-80.2,315.4,,,-81.4,-32.9,78.8,22.7,7.5,27 +3,5,north_america,101,-82.8,316.2,,,-84.2,-24,77.4,12.5,8.6,30 +4,4.4,north_america,101,-81.9,323.6,,,-83.1,-16.9,77.4,12.5,8.6,30 +6,6.2,north_america,101,-80.6,22.1,,,-79.5,43.5,76,6.1,9.7,33 +4,4.3,north_america,101,-79.8,5.9,,,-78.9,31.9,74.9,1.4,11.4,37 +5,7,north_america,101,-78.2,297.5,,,-80.7,-40.3,74.5,0.1,12,38.5 +3,2.4,north_america,101,-85.2,296.6,,,-87.2,-9.1,74.5,-1.1,12.6,40 +6,3,north_america,101,-79.2,326.4,,,-80.1,-1.7,74.4,-2.5,13.4,42 +5,7.7,north_america,101,-83.7,323.7,-76.7,50.4,-73.2,70.3,74.3,-3.7,14.2,44 +5,8,north_america,101,-82.7,333.6,,,-82.6,19.8,74.2,-4.9,15,46 +4,9.6,north_america,101,-79.3,324.3,,,-80.2,0.8,74.2,-4.9,15,46 +4,10.1,north_america,101,-85.2,62,,,-81.2,84,74.5,-4.8,15.3,47 +6,5.6,north_america,101,-72,343.4,,,-71.9,10.7,74.5,-4.8,15.3,47 +5,4,north_america,101,-77.1,325.8,,,-78.1,-0.9,75.9,-3.5,16.2,50 +6,6.2,north_america,101,-83.1,326.3,,,-83.5,15.9,75.9,-3.5,16.2,50 +6,4.9,north_america,101,-77.6,309.1,-78.3,22.7,-75.8,51.7,76.5,-2.8,16.5,51 +6,2.6,north_america,101,-82.7,347.2,,,-81.8,31.2,76.5,-2.8,16.5,51 +6,4.4,north_america,101,-81.4,347.7,-73.9,48.8,-71.3,71.7,79.8,4.1,17.6,55 +5,14,north_america,101,-68.1,9.4,,,-67.3,34.2,80.9,6.2,18.2,57 +6,3,north_america,101,-75.9,326.4,-74.9,24.7,-73.4,52,81.8,4.8,19.4,61 +4,1.1,north_america,101,-77.1,21.3,,,-75.6,50.6,82.2,4,20,63 +4,6.6,north_america,101,-72,5.3,-68,38.1,-66.1,62.3,82.2,4,20,63 +6,3.7,north_america,101,-81.9,350.6,,,-81.7,29.2,82.2,4,20,63 +5,3.9,north_america,101,-80.5,359.4,,,-79.9,34.6,82.4,3.6,20.4,64 +5,5.8,north_america,101,-71.8,28.3,,,-69.2,56.8,81.4,-8.2,22.7,71 +5,6.2,north_america,101,-69.8,354.9,,,-68.6,26.4,81.3,-9.2,23.1,72 +6,4.6,north_america,101,-83.7,15.4,,,-80.7,58.1,80.7,-12.3,24.3,74.5 +6,9,north_america,101,-80.8,358.1,,,-78,43.4,79.5,-15.9,25.7,77 +5,9.6,north_america,101,-59.8,17.7,-57.8,41.8,-52.8,70.3,78.7,-18.1,26.9,79 +6,6.2,north_america,101,-70.5,27.6,,,-65.5,60.8,78.2,-18.8,27.5,80 +7,6.6,north_america,101,-81.8,7.8,,,-77.2,55.4,77.8,-19.4,28.1,81 +5,4.4,north_america,101,-73.6,11.1,,,-61.2,65.2,69.4,-23.5,40.5,100 +5,8.3,north_america,101,-72.4,20.7,,,-59.3,70,69.4,-23.5,40.5,100 +5,4,north_america,101,-76.4,319.4,,,-69.4,45.2,68.8,-23.1,42.6,103 +5,4.4,north_america,101,-71.5,2.4,,,-56,67.6,67.3,-22,48.2,111 +5,6.5,north_america,101,-77.4,5.9,,,-60.7,74.5,67.2,-21.8,48.9,112 +5,5.6,north_america,101,-75.7,29.1,,,-57.1,84.1,67,-21.7,49.6,113 +5,3.6,north_america,101,-74.2,30.3,,,-53.1,88.2,66,-20.6,54.2,120 +5,5.3,north_america,101,-74.8,354.8,,,-56.5,73.9,65.8,-20.2,54.8,121 +5,4.6,north_america,101,-65.9,27.8,,,-44.4,86.2,65.7,-19.8,55.3,122 +6,2.4,north_america,101,-72,11.2,,,-51.7,79.7,65.7,-19.8,55.3,122 +5,3.3,north_america,101,-71.3,7.5,,,-51.5,77.6,65.7,-19.8,55.3,122 +5,7.5,north_america,101,-71,16.9,,,-49.4,83.7,65.4,-18.9,56.8,125 +5,3.6,north_america,101,-67.2,30.8,,,-43.6,91.6,64.9,-17.7,58.7,129 +7,2.6,north_america,101,-58,23.1,,,-33.4,89.2,64.9,-16.7,63.2,144 +6,4.1,north_america,101,-64.2,338.8,-65.5,12.5,-41.3,86.4,64.8,-16.8,64.4,148 +5,3.6,north_america,101,-63.4,344.8,-64,17.5,-39.3,89.3,65.1,-16.1,64.9,150 +6,7,north_america,101,-58.7,315.1,,,-49.8,50.8,65.3,-15.7,65.2,151 +5,5.3,north_america,101,-57.4,327.5,-62.5,351,-42.6,76.4,65.8,-14.2,66.8,156 +6,6,north_america,101,-56.8,328.1,-58.1,350.6,-37.8,76.3,65.8,-13.5,69.6,163 +6,7.4,north_america,101,-52.6,318.2,-58.6,334.6,-42,69.4,65.9,-13,71.5,168 +5,4.3,north_america,101,-64,299,-71.4,321.9,-54.9,77,65.9,-13,71.5,168 +4,7.4,north_america,101,-78.7,270.3,,,-68.5,87.8,65.9,-12.9,71.9,169 +5,7.8,north_america,101,-58,303.8,,,-51.9,53,65.6,-12.9,72.9,172 +4,1.4,north_america,101,-63.2,283.2,,,-61.7,54.7,64.8,-13.4,75,180 +3,1,north_america,101,-75.7,264.7,,,-67.6,85,64.4,-13.6,75.7,183 +6,3.3,north_america,101,-58.8,256.4,-62,257.2,-72,47.8,64.1,-13.8,76.5,186 +5,7,north_america,101,-61.3,264.8,-67.8,268.3,-67.4,64.9,64.1,-13.8,76.5,186 +4,9.8,north_america,101,-58.7,277.8,-64.4,283,-61,59.4,64.1,-13.8,76.5,186 +3,7.2,north_america,101,-73.1,275.4,,,-64.1,78.6,63.8,-14,77.3,189 +5,3.1,north_america,101,-72.8,268.1,,,-66.1,79.2,63.7,-14,77.6,190 +5,8.9,north_america,101,-66,266,,,-66.7,63.8,63.3,-14.2,78.6,194 +4,2.3,north_america,101,-63,263.1,,,-67.6,56.7,63.2,-14.2,79,197 +4,11.1,north_america,101,-65.5,267.5,,,-65.9,63.2,63.2,-14.2,79,197 +6,2.8,north_america,101,-62.5,251,,,-73.2,55.4,63.2,-14.2,79,197 +5,4,north_america,101,-68,268.5,,,-65.4,69.2,63.2,-14.2,79,197 +5,6.2,north_america,101,-63.6,268.7,,,-65.2,58.6,63.2,-14.2,79,197 +5,6,north_america,101,-55.6,274.6,-59.8,273.3,-62.2,51.5,63.2,-14.2,79.1,198 +5,7.9,north_america,101,-61.5,234,,,-81,61.6,63.2,-14.2,79.2,199 +6,10.7,north_america,101,-66.4,252,,,-71.9,68.2,63.2,-14.1,79.2,200 +7,3.2,north_america,101,-66.6,268.2,,,-65.5,66.3,63.2,-14.1,79.3,201 +6,5,north_america,101,-67.8,275.8,,,-62.5,69.5,63.2,-14,79.5,204 +5,8,north_america,101,-58.5,256.9,-61.2,257.1,-70.3,51.8,63.2,-14,79.5,204 +6,10.7,north_america,101,-58.7,250.9,-62.8,251.2,-73,57.3,63.2,-14,79.5,204 +6,4.2,north_america,101,-57.8,259.3,-59.6,259.5,-68.9,47.8,63.2,-14,79.5,204 +6,6.5,north_america,101,-64.9,276.6,,,-61.9,63.3,63.2,-14,79.5,204 +6,2.5,north_america,101,-58.1,271.8,,,-62.3,48.4,63.2,-13.9,79.7,207 +6,5,north_america,101,-66.9,267.2,,,-65.8,67.6,63.2,-13.9,79.7,207 +6,4.7,north_america,101,-65.5,255.1,,,-70.7,65.7,63.2,-13.9,79.9,211 +7,5.6,north_america,101,-55.6,274.6,-59.9,273.3,-62.1,52.8,63.2,-13.9,79.9,211 +6,3.4,north_america,101,-56.6,255.9,-58.6,256.2,-70.3,45,63.2,-13.9,79.9,211 +6,4,north_america,101,-61.7,274.7,,,-61.9,57,63.2,-13.9,79.9,211 +5,3,north_america,101,-57.6,269.6,-63.3,266.9,-65.7,59,63.2,-13.9,79.9,211 +4,5.2,north_america,101,-61.4,282.2,,,-58.4,58.5,63.2,-13.9,79.9,212 +6,3.1,north_america,101,-60.1,277.1,,,-60.3,54.5,63.2,-13.9,79.9,214 +4,7,north_america,101,-60.1,271.8,,,-62.9,52.7,63.2,-13.9,79.9,215 +3,10,north_america,101,-59,267.6,,,-64.6,49,63.2,-13.9,79.9,215 +3,14,north_america,101,-56.1,276,-62.4,280.5,-59.5,59.9,63.2,-13.9,79.9,215 +7,7.8,north_america,101,-50.5,267.6,-53.4,268.7,-62.1,37.9,63.2,-13.9,79.9,215 +6,5.6,north_america,101,-57.4,267.7,-59.3,268.4,-64.3,50,63.2,-13.9,79.9,216 +6,3.2,north_america,101,-59.9,279.4,,,-59.2,54.9,63.2,-13.9,79.9,217 +6,5.1,north_america,101,-52.9,282,-53.5,282.5,-54.9,45.6,63.2,-13.9,79.9,218 +6,7.7,north_america,101,-56.4,276.8,-58,277.6,-59.4,50.7,63.2,-13.9,79.9,218 +7,5,north_america,101,-59.6,277.5,-63.5,280.6,-59.7,61.9,63.2,-13.9,79.9,220 +6,5,north_america,101,-53.4,281.7,-55.9,281.5,-56.6,48.9,63.2,-13.9,79.9,220 +6,2,north_america,101,-58.5,279.5,,,-58.6,52.4,63.2,-13.9,79.9,221 +5,5,north_america,101,-60.5,281.6,,,-58.4,56.7,63.2,-13.9,79.9,221 +4,3.9,north_america,101,-48.3,272.3,,,-57.6,31.5,63.2,-13.9,79.9,221 +6,5,north_america,101,-49.1,285.1,-52.7,287.8,-51.7,47.3,63.2,-13.9,79.9,225 +6,4,north_america,101,-54.1,285.2,,,-53.8,47.9,63.2,-13.9,79.9,227 +6,3.2,north_america,101,-48.4,278.5,,,-54.2,35.8,63.2,-13.9,79.9,228 +5,2.5,north_america,101,-54.1,288.3,-58.6,293.1,-52.3,58,63.2,-13.9,79.9,230 +6,4,north_america,101,-45.2,295.4,-49.1,299,-44,49.4,63.2,-13.9,79.9,230 +5,3.1,north_america,101,-52.5,290.7,-56.6,295,-50.4,56.1,63.2,-13.9,79.9,230 +6,4.5,north_america,101,-56.5,283.2,-59.8,286.6,-55.8,57.3,63.2,-13.9,79.9,230 +6,3.3,north_america,101,-46.1,293.6,-49.8,296.9,-45.5,48.9,63.2,-13.9,79.9,230 +7,3.4,north_america,101,-55.5,287.9,-58.2,290.7,-53.2,56.5,63.2,-13.9,79.9,234 +7,5,north_america,101,-54.6,284.5,-57.8,287.8,-54.4,54.5,63.2,-13.9,79.9,234 +6,4.9,north_america,101,-44.7,301.4,-44.9,301.7,-39.8,46.9,63.2,-13.9,79.9,234 +5,4.9,north_america,101,-55.6,285.8,-60.5,290.9,-54.2,59.9,63.2,-13.9,79.9,234 +4,5.3,north_america,101,-41.1,305.6,-40.8,305.2,-35,45.7,63.2,-13.9,79.9,234 +5,7,north_america,101,-46.1,301,-50.6,306.1,-41.5,55.1,63.2,-13.9,79.9,235 +6,7.2,north_america,101,-44.3,271.6,-45.5,271.1,-56.6,26.5,63.2,-13.9,79.9,246 +6,5,north_america,101,-51,306.5,-53.8,310.3,-41.8,60.5,63.2,-13.9,79.9,250 +4,8,north_america,101,-49.9,298.1,-51.6,300,-45.1,52.6,63.2,-13.9,79.9,255 +4,15,north_america,101,-54.8,299.3,-57.3,302.2,-47.6,60.4,63.2,-13.9,79.9,258 +5,5,north_america,101,-51.5,306.7,-54.8,311.6,-41.9,62.2,63.2,-13.9,79.9,263 +5,3.8,north_america,101,-56.3,302.9,,,-46.6,59.5,63.2,-13.9,79.9,270 +4,3.8,north_america,101,-53,308.7,,,-42,58.9,63.2,-13.9,79.9,272 +4,8.6,north_america,101,-54.8,292.1,,,-50.8,52.3,63.2,-13.9,79.9,272 +3,10,north_america,101,-51.9,303,-56.5,311,-43.2,63.6,63.2,-13.9,79.9,277 +3,5,north_america,101,-51.7,302.1,-53.7,304,-44.5,57.1,63.2,-13.9,79.9,277 +7,3.6,north_america,101,-42.1,306.5,-41.4,306.4,-34.9,46.9,63.2,-13.9,79.9,280 +5,16.3,north_america,101,-33.5,306.3,,,-29,40.8,63.2,-13.9,79.9,282 +4,13.1,north_america,101,-44.6,305.3,-47.8,308.8,-38.3,54,63.2,-13.9,79.9,283 +5,2.1,north_america,101,-46.8,304,-48.3,305.9,-40.1,52.7,63.2,-13.9,79.9,285 +3,10.2,north_america,101,-38.9,300.8,,,-35.9,41,63.2,-13.9,79.9,289 +5,1.5,north_america,101,-50.5,303,-56.2,310.3,-43.4,63,63.2,-13.9,79.9,291 +4,5,north_america,101,-37.5,296.6,-35.8,295.2,-36.7,34.4,63.2,-13.9,79.9,292 +5,7.1,north_america,101,-40.1,307.7,-42.1,310.1,-33.5,50,63.2,-13.9,79.9,292 +5,2,north_america,101,-43.1,307.9,-46.5,311.7,-36,54.6,63.2,-13.9,79.9,292 +4,2,north_america,101,-41.6,300.4,-42.6,301.5,-38.3,44.7,63.2,-13.9,79.9,292 +4,2.8,north_america,101,-40.1,300.5,-40.3,300.8,-37,42.2,63.2,-13.9,79.9,298 +4,12.8,north_america,101,-55.3,279.8,-60.6,285.1,-56.8,58,63.2,-13.9,79.9,299 +5,3.9,north_america,101,-44.1,301.5,-41.5,300.4,-38,43,63.2,-13.9,79.9,300 +5,2.1,north_america,101,-42.1,312.1,-43,313.4,-32.7,52.9,63.2,-13.9,79.9,301 +5,3.4,north_america,101,-44.1,303.9,-46.3,306.8,-38.2,51.4,63.2,-13.9,79.9,301 +6,3.1,north_america,101,-28.6,299.9,,,-28.6,32.4,63.2,-13.9,79.9,303 +5,1.8,north_america,101,-45.7,308.6,-50.5,314.6,-37.6,59.8,63.2,-13.9,79.9,303 +5,6,north_america,101,-36,302,-30.2,301.5,-29,34.8,63.2,-13.9,79.9,310 +7,7.7,north_america,101,-27.2,298.3,,,-28.4,30.2,63.2,-13.9,79.9,317 +6,8.3,north_america,101,-22.6,294.4,,,-26.9,23.8,63.2,-13.9,79.9,320 +7,15.3,north_america,101,-27.9,297.2,,,-29.5,29.8,63.2,-13.9,79.9,322 +4,6.5,north_america,101,-19.5,315.8,,,-12.6,39.2,63.2,-13.9,79.9,330 +6,8,north_america,101,-27,311,-17.8,309.8,,,,,,333 +7,9,north_america,101,-18.6,304.2,,,,,,,,335 +3,16,north_america,101,-27.4,303,-16.6,299.6,,,,,,370 +4,9,north_america,101,-13,285,1.5,284.8,,,,,,415 +7,5.3,north_america,101,-17,305,,,,,,,,420 +6,5.8,north_america,101,-19.1,308.3,,,,,,,,425 +4,7.3,north_america,101,-24,326.6,-16.9,321.7,,,,,,438 +4,3.9,north_america,101,-13.4,329.3,,,,,,,,470 +4,4.3,north_america,101,-17.5,332.4,,,,,,,,480 +6,11.9,north_america,101,-10.3,346.5,,,,,,,,490 +5,8.5,north_america,101,0.6,343,3.1,338.9,,,,,,495 +5,9.7,north_america,101,-10.6,338,-8.4,334.6,,,,,,495 +5,9,north_america,101,-5.2,345.8,-4.7,345,,,,,,495 +5,7.1,north_america,101,3.4,355.1,,,,,,,,500 +6,4.3,north_america,101,-12.6,337.3,,,,,,,,500 +6,10,north_america,101,5.4,348.7,,,,,,,,503 +5,3.3,north_america,101,-0.6,341.1,-1.7,342.6,,,,,,508 +5,6.2,north_america,101,11.9,4.5,,,,,,,,532 +6,7.4,greenland,102,-76.3,21.5,,,-75,44.1,72.8,9.1,11.2,39 +5,6.3,greenland,102,-74.6,339.4,,,-77.6,11.6,71.1,30.9,16,54 +4,6,greenland,102,-63,0,,,-64.6,25.6,71.4,30.2,16.4,54.5 +4,15,greenland,102,-63,354,,,-65,19.8,71.7,29.6,16.8,55 +3,4.2,greenland,102,-61,345,,,-63.8,9.9,72,29,17.2,55.5 +4,8.9,greenland,102,-63.4,5.1,,,-64.4,32,72,29,17.2,55.5 +6,3,greenland,102,-68,358,,,-69.6,28,72,30.1,17.5,59 +5,9,greenland,102,-56,3,,,-57.4,27.7,71.9,30.3,17.5,59.5 +6,3.2,greenland,102,-67.5,15,,,-67.6,44.9,71.8,30.8,17.5,60.5 +5,6,greenland,102,-76.2,37.9,,,-73.9,73.2,71.8,30.8,17.5,60.5 +5,6.9,greenland,102,-64.8,321.5,,,-69.4,-14.9,71.8,30.8,17.5,60.5 +6,5,greenland,102,-52.7,278.7,,,-67.9,36.6,60.5,1.4,69.3,209 +3,11,greenland,102,1,302,,,,,,,,400 +4,3.6,europe,301,-80.6,267.5,,,-80.6,-92.3,18.3,-47,0.1,0.5 +4,4.4,europe,301,-86.4,296.1,,,-86.4,-63.1,18.3,-47,0.1,0.5 +3,12.9,europe,301,-84.3,357.7,,,-83.9,5.8,18,-26.7,1,8 +5,1.8,europe,301,-78.9,328.3,,,-78.9,-25.8,17.9,-26.4,1.1,9.5 +5,3.5,europe,301,-77.4,314.2,,,-77.7,-40.4,18.5,-26.3,1.2,10 +5,6.9,europe,301,-84.1,351.2,,,-83.6,3,19.7,-25.7,1.4,11.5 +5,1.5,europe,301,-72.4,352,,,-71.9,-3.5,19.6,-25.5,1.5,12 +4,4.4,europe,301,-77.8,310.8,,,-78.9,-36.2,20,-19.7,2.9,24 +3,3.4,europe,301,-80.8,2,,,-78,26.5,26.8,-19.5,5.5,34 +5,6.8,europe,301,-63.7,358.6,-76.7,357.2,-72.2,30.6,30.5,-15.7,10.3,52 +5,1.5,europe,301,-83,335,,,-78.2,41.3,34.5,-15.7,12.4,59 +5,3.5,europe,301,-76,340,,,-72.9,22.7,34.5,-15.7,12.4,59 +5,1.2,europe,301,-81.7,359.8,,,-74.9,46.3,34.6,-15.7,12.5,59.4 +5,2.8,europe,301,-80.2,339.6,,,-76,33.2,34.6,-15.7,12.5,59.5 +4,1.5,europe,301,-82.5,338,,,-77.6,40.6,34.7,-15.7,12.5,59.8 +5,5.7,europe,301,-77.7,325.4,,,-76.2,18,34.7,-15.7,12.5,59.8 +4,2.7,europe,301,-77,355,,,-71.5,34.4,34.8,-15.7,12.6,60 +6,4.5,europe,301,-71.4,334.7,,,-69.8,11.3,34.8,-15.7,12.6,60.1 +6,2.4,europe,301,-81,359,,,-74.3,44.8,34.9,-15.7,12.7,60.7 +6,4.7,europe,301,-73.3,346.2,,,-69.5,23.1,35,-15.7,12.8,60.9 +5,4.7,europe,301,-78.9,347,,,-74,34.1,35,-15.7,12.8,61 +5,2.7,europe,301,-74,351,,,-69.4,27.5,35.1,-15.7,12.8,61.2 +6,3.4,europe,301,-79.1,344.5,-84.9,222.9,-82.4,95.9,35.8,-16,14.6,68 +7,8,europe,301,-73,336,,,-69.8,21.1,35.4,-15.7,15.2,74 +5,3,europe,301,-74,341,,,-67.7,35.4,37.5,-14,19.5,86 +5,3,europe,301,-74,328,,,-70.4,28.7,37.5,-14,19.5,86 +4,5,europe,301,-68,329,,,-65.1,22.3,39.2,-13.6,21.8,89.5 +5,4,europe,301,-76,1,,,-62.9,54.7,39.4,-13.7,24.2,93 +6,2.5,europe,301,-80.8,338.4,,,-63.2,70.6,41.6,-10.1,34.1,108 +6,2.9,europe,301,-74,3,,,-48,82.4,45.9,-3.4,48.4,140 +5,6,europe,301,-78,328,,,-56,87.9,48.4,-1.1,52.8,156.5 +3,3.9,europe,301,-70,327,,,-54.5,74.3,48.5,-1,53,157 +5,7,europe,301,-78,310,,,-59.5,88.7,48.6,-0.9,53.4,158 +4,4,europe,301,-72,312,,,-59.3,77.2,48.7,-0.8,53.8,159 +6,7.3,europe,301,-72,330,,,-53.9,78.9,48.7,-0.8,53.8,160 +7,6.3,europe,301,-63,300,,,-62.2,62.5,49.8,0.3,57,168 +6,6.8,europe,301,-69,283,,,-65.8,82.6,49.4,0.2,60.3,179 +6,12,europe,301,-71,276,,,-65.6,90.9,49.2,0,61.8,184 +5,5.1,europe,301,-66,295.2,,,-60.6,75.2,49.1,-0.1,62.3,186 +4,3,europe,301,-77,315,,,-52.1,94.8,48.7,0,63.6,192 +5,7.5,europe,301,-61,259,,,-76,88.7,48.1,0.7,63.1,198 +6,9,europe,301,-55,280,,,-69.1,54,47.8,1.1,62.5,201 +7,3,europe,301,-51,285,,,-65.8,44.1,47.8,1.1,62.5,201 +6,4.5,europe,301,-50,286.4,,,-65.3,41.7,47.6,1.6,61.9,204.2 +6,8,europe,301,-50,292,-58,272.9,-73.6,65.6,47.3,2.2,61.1,208 +5,5.1,europe,301,-50,308,-58.2,298.2,-61.7,60.7,46.6,3.2,59.5,215 +5,4.6,europe,301,-50,305,,,-56.9,45.8,46,3.9,58.2,221 +6,6,europe,301,-49,311,,,-52.9,46.9,46,3.9,58.2,226 +5,2.9,europe,301,-47.1,301.6,,,-57.8,39.2,46,3.9,58.2,228 +6,3,europe,301,-54,321,,,-49.4,58.3,46,3.9,58.2,234 +6,12,europe,301,-53,303,,,-59.1,50.4,46,3.9,58.2,234 +5,15,europe,301,-49,326,-56.5,318.6,-51.6,61.2,46,3.9,58.2,239 +6,5,europe,301,-43,326,-47.8,322.2,-45.8,50.6,46,3.9,58.2,243 +6,3.8,europe,301,-49,348.2,-57,344.7,-39.2,71.2,46,3.9,58.2,246 +6,7.8,europe,301,-59.3,325.8,,,-49,67.4,46,3.9,58.2,248 +7,3.3,europe,301,-50.6,345.6,-58.8,341,-41.8,71.7,46,3.9,58.2,249 +3,10,europe,301,-59,330,,,-46.9,68.2,46,3.9,58.2,250 +6,3.3,europe,301,-56.2,326,,,-47.7,63.2,46,3.9,58.2,251 +3,13.9,europe,301,-52.7,328.4,,,-44.8,59.6,46,3.9,58.2,251 +5,2.1,europe,301,-56.4,321.7,,,-49.9,62,46,3.9,58.2,251 +6,5,europe,301,-50,343,-59.3,336.5,-44.1,70.7,46,3.9,58.2,251 +3,5.3,europe,301,-53.3,330.2,,,-44.2,61.1,46,3.9,58.2,251 +3,5,europe,301,-54.3,323,,,-48.4,59.5,46,3.9,58.2,251 +3,2.7,europe,301,-58.5,314.5,,,-54.3,63.2,46,3.9,58.2,251 +6,9.7,europe,301,-52.8,334.4,,,-41.8,62.3,46,3.9,58.2,251 +5,2.7,europe,301,-46,327,-50,324.7,-45.5,54.7,46,3.9,58.2,255 +6,4,europe,301,-51,341,-55.4,338.1,-41.4,66.8,46,3.9,58.2,255 +5,3.5,europe,301,-45.6,350.2,-54.9,338.2,-41,66.2,46,3.9,58.2,260 +5,5,europe,301,-47,331,-50.6,327.7,-44.2,56.7,46,3.9,58.2,261 +6,4,europe,301,-49,343,-52.6,341.7,-38,65.3,46,3.9,58.2,261 +4,0,europe,301,-53,331,-58.6,326.4,-48.5,66.5,46,3.9,58.2,264 +5,1.5,europe,301,-49,334,-52.5,331.5,-43.1,60.7,46,3.9,58.2,264 +3,4.6,europe,301,-47,336,-49.3,334.7,-39.6,58.5,46,3.9,58.2,264 +5,4,europe,301,-51,324,-56.4,318.8,-51.4,61,46,3.9,58.2,264 +5,6.1,europe,301,-51.5,322,,,-47.7,55.3,46,3.9,58.2,264 +4,4.1,europe,301,-50,344,,,-35.3,63.8,46,3.9,58.2,269 +5,2.5,europe,301,-51,343,,,-36.4,64.3,46,3.9,58.2,271 +5,5.9,europe,301,-53,344,,,-37.1,66.7,46,3.9,58.2,271 +4,8.6,europe,301,-51,345,,,-35.4,65.3,46,3.9,58.2,275 +4,11,europe,301,-54,352,,,-34.1,71.4,46,3.9,58.2,279 +4,7,europe,301,-37,341,,,-28.3,50.6,46,3.9,58.2,280 +3,14,europe,301,-47,337,,,-37.1,57.3,46,3.9,58.2,280 +5,10,europe,301,-42,346,,,-29,57.8,46,3.9,58.2,280 +5,1,europe,301,-47,337,,,-37.1,57.3,46,3.9,58.2,281 +4,13.4,europe,301,-44.6,337.4,,,-35.4,55.2,46,3.9,58.2,281 +4,6.9,europe,301,-38,346,,,-26.2,54.6,46,3.9,58.2,281 +5,6.5,europe,301,-49.4,359.7,,,-27.4,71.4,46,3.9,58.2,282.6 +4,6.7,europe,301,-41,352,,,-25.1,60.6,46,3.9,58.2,285 +5,3.2,europe,301,-43,352,,,-26.5,62.1,46,3.9,58.2,285 +5,5.1,europe,301,-44,4,-48.6,4.1,-24.9,73.2,46,3.9,58.2,285 +5,2,europe,301,-40,346,-42.4,345.2,-29.6,57.7,46,3.9,58.2,285 +3,7.7,europe,301,-44,350,-39.8,-9.4,-25,58.8,46,3.9,58.2,285 +5,6.3,europe,301,-38,346,,,-26.2,54.6,46,3.9,58.2,285 +4,8.1,europe,301,-42,354,,,-24.8,62.5,46,3.9,58.2,285 +5,2,europe,301,-42,349,-40.5,-9.8,-25.7,59,46,3.9,58.2,285 +4,17,europe,301,-49,342,-52.5,340.1,-38.7,64.5,46,3.9,58.2,285 +4,6.8,europe,301,-37,340,-38.5,339.4,-30.2,50.9,46,3.9,58.2,285 +4,7.9,europe,301,-43,345,,,-30.2,58.1,46,3.9,58.2,285 +5,4,europe,301,-41,345,-44,343.5,-31.6,58.1,46,3.9,58.2,285 +3,13.2,europe,301,-40,352,,,-24.3,59.8,46,3.9,58.2,285 +5,4,europe,301,-50,330,,,-42.6,57.1,46,3.9,58.2,286 +4,5.9,europe,301,-49,356,,,-28.8,69.1,46,3.9,58.2,286 +4,10,europe,301,-48,343,,,-34.5,61.4,46,3.9,58.2,286 +3,1,europe,301,-42,353,,,-25.3,61.9,46,3.9,58.2,286 +4,5.8,europe,301,-41.5,340,-45.3,338.1,-35.5,56.3,46,3.9,58.2,287 +4,2.4,europe,301,-32,354,,,-17.3,55.4,46,3.9,58.2,291 +3,15.9,europe,301,-41,349,,,-26.7,58.8,46,3.9,58.2,291 +3,13,europe,301,-46,347,,,-31.2,61.8,46,3.9,58.2,291 +6,13,europe,301,-42,346,,,-29,57.8,46,3.9,58.2,293 +4,4.8,europe,301,-44,339,,,-34.1,55.5,46,3.9,58.2,294 +5,3.5,europe,301,-32.9,347.1,,,-21.9,51.4,46,3.9,58.2,294 +5,6.3,europe,301,-35.4,346.8,,,-23.9,53.1,46,3.9,58.2,294 +4,4,europe,301,-47,348,,,-31.3,63.2,46,3.9,58.2,294 +3,19,europe,301,-42,348,,,-27.9,59,46,3.9,58.2,294 +4,4.8,europe,301,-44,355,,,-25.7,64.6,46,3.9,58.2,294 +5,8.1,europe,301,-47.1,337.1,,,-37.1,57.5,46,3.9,58.2,294 +5,6.5,europe,301,-38,348,,,-25.1,55.8,46,3.9,58.2,294 +5,11,europe,301,-37,354,,,-21.1,58.9,46,3.9,58.2,294 +4,7.1,europe,301,-37.1,350,,,-23.3,56.4,46,3.9,58.2,295 +4,13.6,europe,301,-43,354,,,-25.5,63.3,46,3.9,58.2,296 +3,7.1,europe,301,-42.5,339.6,,,-32.8,54.5,46,3.9,58.2,297 +4,2.9,europe,301,-39,341,,,-29.7,52.3,46,3.9,58.2,297 +5,1.3,europe,301,-41,342,,,-30.5,54.6,46,3.9,58.2,297 +6,3,europe,301,-43,345,-49.9,337.3,-38.6,60.4,46,3.9,58.2,297 +5,2.4,europe,301,-48.4,349.8,-56.1,341.2,-40.3,68.8,46,3.9,58.2,299 +5,4,europe,301,-31,354,,,-16.6,54.7,46,3.9,58.2,299 +5,2.2,europe,301,-48.2,348.3,-55.9,339.4,-41,67.9,46,3.9,58.2,299 +6,4,europe,301,-42,359,-45.7,356.5,-26.3,66.8,46,3.9,58.2,301 +5,2,europe,301,-48.2,342.3,,,-35,61.3,46,3.9,58.2,303 +3,3,europe,301,-49,349,,,-32.2,65.4,46,3.9,58.2,303 +5,5.2,europe,301,-38.3,354,,,-22,59.8,46,3.9,58.2,305 +5,9,europe,301,-38,343,-40.3,341.5,-30.3,53.7,46,3.9,58.2,305 +5,2.9,europe,301,-38.4,339.5,,,-30.1,50.9,46,3.9,58.2,312 +6,7,europe,301,-14,332,,,,,,,,332 +6,8.2,europe,301,-14.3,335.9,,,,,,,,335 +6,11,europe,301,4,323,,,,,,,,396 +6,7,europe,301,-5,320,,,,,,,,410 +7,3,europe,301,2,318,,,,,,,,410 +6,7.3,europe,301,3.7,325.5,,,,,,,,411 +5,6,europe,301,2,321,,,,,,,,412 +4,2.5,europe,301,-8,335,,,,,,,,415 +4,4.3,europe,301,-1,325,,,,,,,,415 +4,4.9,europe,301,-17,350,,,,,,,,419 +7,9.1,europe,301,-19,344,-13.6,-16,,,,,,421 +3,5.7,europe,301,-27,346,,,,,,,,421 +3,8,europe,301,-23,351,,,,,,,,422 +4,4.6,europe,301,-15,350,,,,,,,,424 +5,7.9,europe,301,-15,347,,,,,,,,425 +3,6,europe,301,-21,344,,,,,,,,425 +4,2,europe,301,-19,349,,,,,,,,426 +5,5.1,europe,301,-19,352,,,,,,,,427 +3,5,europe,301,-21,344,,,,,,,,430 +4,5.2,europe,301,-21,358,,,,,,,,432 +5,13.4,europe,301,3,35,,,,,,,,458 +6,4.8,europe,301,5,34,,,,,,,,459 +5,4.9,europe,301,12,41.9,,,,,,,,463 +6,4.4,europe,301,14,49,,,,,,,,466 +4,6.7,europe,301,11.3,39.1,,,,,,,,467 +6,9.2,europe,301,19,51,,,,,,,,471 +6,6.8,europe,301,18.7,54,,,,,,,,472 +4,7.1,europe,301,25,50.9,,,,,,,,472 +5,9,europe,301,30,55,,,,,,,,475 +6,5.1,europe,301,18,46,,,,,,,,475 +6,4,europe,301,18,55,,,,,,,,477 +5,3.6,europe,301,34.7,59.1,,,,,,,,478 +5,5,europe,301,22,87,34.2,79.8,,,,,,500 +3,6.8,europe,301,52,111,,,,,,,,500.1 +4,6.9,europe,301,58.4,122.5,68.7,102.2,,,,,,535 +3,10,Amazonia,201,-85.4,303.8,,,-85.5,-37.4,62.1,-40.5,2.8,8.5 +3,10,Amazonia,201,-85.7,80.5,,,-84.6,75.6,62.1,-40.5,2.8,8.5 +5,11.4,Amazonia,201,-84.5,316,,,-84.4,-26.3,62,-40.6,3.1,9.5 +5,5.9,Amazonia,201,-79.5,0,,,-69.6,46.8,63.5,-33.4,26.2,70.5 +6,2.6,Amazonia,201,-83.2,320.1,,,-71.8,48.8,61.6,-34.3,33.8,84 +5,4.2,Amazonia,201,-79.4,331.9,,,-67.5,44.3,61.2,-34.3,34.4,85 +3,4.8,Amazonia,201,-87.6,315.1,,,-69.5,65.2,58.2,-34.6,38.7,92 +4,3,Amazonia,201,-87.9,335.9,,,-63.5,69.8,54.9,-34.8,44.9,102 +4,2.8,Amazonia,201,-83.6,81,,,-53.7,84.1,51.8,-35,52.4,118 +6,2.6,Amazonia,201,-89.1,3.3,,,-56.3,76.3,51,-34.3,53.6,123.5 +6,2,Amazonia,201,-82.4,30.3,,,-48.1,77.9,50.1,-32.8,54.8,130 +6,2.4,Amazonia,201,-83,71.4,,,-49.6,85.6,50,-32.5,55.1,132 +4,14.1,Amazonia,201,-80.6,275,,,-59.3,63.3,50,-32.5,55.1,146 +3,9.3,Amazonia,201,-85.3,262.5,,,-58.5,72.7,50,-32.5,55.1,175 +5,3.8,Amazonia,201,-65.5,250.3,,,-70.6,34.3,50,-32.5,55.1,196.6 +5,4,Amazonia,201,-81.2,235.1,,,-63.6,72.7,50,-32.5,55.1,198 +4,4.9,Amazonia,201,-66.9,245.6,,,-71.9,40.3,50,-32.5,55.1,202.5 +4,10,Amazonia,201,-82,320,,,-52.7,66.4,50,-32.5,55.1,232 +6,6,Amazonia,201,-71.4,303.6,-60,294.1,-50.4,29.1,50,-32.5,55.1,248.5 +6,6.6,Amazonia,201,-80.7,7,-70.7,325.5,-45.5,52.9,50,-32.5,55.1,260 +4,4.5,Amazonia,201,-68.2,321.3,-56.1,305.2,-43.4,29.2,50,-32.5,55.1,280 +4,4.1,Amazonia,201,-65.7,330.9,-53.2,324,-33.7,36,50,-32.5,55.1,300 +3,11.2,Amazonia,201,-57,357,-48.3,338.2,-24.3,41.3,50,-32.5,55.1,310 +3,6,Parana,202,-77,18,,,-66.8,52.3,63.9,-33.6,24.7,65.5 +4,3.7,Parana,202,-84.6,115.4,,,-57,86.1,51.7,-35,52.6,119 +3,10.4,Parana,202,-81,14,,,-49.2,71.8,51.5,-34.9,53,121 +5,5.9,Parana,202,-72,25,,,-39.1,73.6,50.9,-34.2,53.7,124 +3,11,Parana,202,-79,8,,,-43,71.3,47.5,-33.3,56.1,139.5 +3,18.1,Parana,202,5.9,338.1,,,24.1,12.1,47.6,-33.3,56.2,510 +4,8,Colorado,290,-85,222,,,-73,71.4,56.9,-34.7,40.9,95.5 +5,8.7,Colorado,290,-83,138,,,-53.3,90.1,47.5,-33.3,57.3,183 +6,4.5,Colorado,290,-74,67,,,-38.1,89.4,47.5,-33.3,57.3,183 +5,6.8,Colorado,290,-75.5,129.4,,,-51,102,47.5,-33.3,57.3,187 +6,4.5,Colorado,290,-51,223,,,-84.7,-21.3,47.5,-33.3,57.3,195 +4,7.6,Colorado,290,-81.8,298.3,,,-52.4,65.2,47.5,-33.3,57.3,216 +4,7,Colorado,290,-83,317,-69.1,298.5,-49.7,45.5,47.5,-33.3,57.3,240 +6,6.4,Colorado,290,-76,313.4,,,-48.1,57.7,47.5,-33.3,57.3,240 +6,4.9,Colorado,290,-89.2,346.1,-75.1,293.5,-52.6,54.2,47.5,-33.3,57.3,240 +6,3.3,Colorado,290,-80.1,348.6,,,-44.9,68.6,47.5,-33.3,57.3,263 +6,4.1,Colorado,290,-75.7,326,,,-45.3,59.5,47.5,-33.3,57.3,264 +4,5.2,Colorado,290,-80.6,308.3,,,-50.7,63.7,47.5,-33.3,57.3,267 +4,2.8,Colorado,290,-80.6,268.8,-66.7,285.5,-53.7,39.6,47.5,-33.3,57.3,283 +3,2.5,Colorado,290,-59.5,357.5,-55.2,332.3,-29.2,43.2,47.5,-33.3,57.3,283 +3,3.1,Colorado,290,-74,313,-60.9,301.3,-45.1,35.2,47.5,-33.3,57.3,283 +4,4.9,Colorado,290,-76.8,293.7,-62.4,292.2,-49.6,34.4,47.5,-33.3,57.3,290 +4,7,Colorado,290,-73.1,272.4,,,-58.4,50.8,47.5,-33.3,57.3,290 +3,5,Colorado,290,-66,348,,,-33.4,57.8,47.5,-33.3,57.3,300 +4,5.7,Colorado,290,-51,347,,,-20.8,48.8,47.5,-33.3,57.3,310 +6,5,Colorado,290,-49,343,-45.4,323.6,-25.1,31.5,47.5,-33.3,57.3,310 +6,9.6,Colorado,290,-57,350,,,-25.2,53.6,47.5,-33.3,57.3,310 +5,5.7,Patagonia,291,-81,337.4,,,-74.4,33.2,57.2,-31.2,19.4,47 +6,4.3,Patagonia,291,-78.7,358.4,,,-68.8,45.5,63.4,-33.4,26.5,71.5 +5,2,Patagonia,291,-86.8,35.2,,,-66.7,71.5,58.2,-34.6,38.7,92 +4,3.8,Patagonia,291,-87,159,-77.5,283.2,-60.9,51.4,52,-35,52.1,116 +5,5.5,Patagonia,291,-81,172,,,-57.9,90.4,47.5,-33.3,58,156.5 +4,4.9,Patagonia,291,-85,197,,,-55.8,82,47.5,-33.3,59.2,167 +5,5.2,Southern_Africa,701,-64.1,46.1,,,,,,,,90.5 +6,9.7,Southern_Africa,701,-47.6,89.9,,,,,,,,129 +5,3.1,Southern_Africa,701,-48.3,86.6,,,,,,,,132 +4,13.3,Southern_Africa,701,-64,80.6,,,,,,,,180 +3,15.8,Southern_Africa,701,-61.9,71.9,,,,,,,,183 +5,3.2,Southern_Africa,701,-71.6,93.5,,,,,,,,183 +4,11,Southern_Africa,701,-70.5,88.7,,,,,,,,183 +3,7,Southern_Africa,701,-57,84,,,,,,,,183 +5,9.5,Southern_Africa,701,-65.4,75.1,,,,,,,,183 +5,8.7,Southern_Africa,701,-70.7,106.7,,,,,,,,186 +3,4.6,Southern_Africa,701,-68,50.5,-54.7,39.5,,,,,,221.5 +3,6,Southern_Africa,701,-54,80,-49,62.6,,,,,,248.5 +5,7.6,Southern_Africa,701,-50.9,86.3,-48,63.8,,,,,,251 +6,5.6,Southern_Africa,701,-26.5,26.5,-21.6,27.7,,,,,,281.5 +6,12,Southern_Africa,701,-25,67,-25.2,53.6,,,,,,315 +4,7,Southern_Africa,701,10,15,-3.3,16,,,,,,398.5 +5,18,Southern_Africa,701,25,343,17,-11,,,,,,446.5 +6,9,Southern_Africa,701,28,14,13.9,13.8,,,,,,482.5 +4,14,Meseta,707,-46,78,,,,,,,,120 +6,9,Meseta,707,-44,71,,,-42.6,73.6,33.6,26,2.3,173.5 +4,11,Meseta,707,-45,68,,,-43.7,70.7,33.6,26,2.3,173.5 +5,6,Meseta,707,-69.2,55.5,,,-68.2,61,33.6,26,2.3,201 +3,7,Meseta,707,-71,36,,,-70.5,42.7,33.6,26,2.3,201 +5,19,Meseta,707,-73,61.3,,,-71.8,67.4,33.6,26,2.3,201 +3,4.6,Meseta,707,-38.7,56.8,,,-37.7,59.4,33.6,26,2.3,273 +5,4.7,Meseta,707,-32.2,64.1,-33.4,66,-32.1,68.3,33.6,26,2.3,273 +3,7.8,Meseta,707,-24,63.8,-23.3,62.3,-22.1,64.3,33.6,26,2.3,273 +3,20.9,Meseta,707,-36,58,,,-34.9,60.5,33.6,26,2.3,280.5 +3,4.1,Somalia,709,-87.5,359.3,,,-87.5,-2.3,52.4,6.3,-0.1,1 +3,5.7,Somalia,709,-87.2,37.1,,,-87.3,34.3,50.2,6.5,-0.2,2 +3,4.1,Somalia,709,-79.7,350.2,,,-79.6,-11,50.3,6.4,-0.3,2.5 +5,3.8,Somalia,709,-85.7,75.8,-84.4,65.5,-85.1,59.6,50.4,6.4,-1.3,11.5 +4,8.8,Somalia,709,-86.5,6.6,,,-86.4,-7.5,50.4,6.4,-1.3,12 +3,10,Somalia,709,-80.1,214.2,,,-79.7,-142.7,50.3,6.4,-1.3,13.5 +5,3.1,Somalia,709,-84.6,343.3,,,-84.2,-25.5,50.4,6.4,-1.3,17 +4,8.4,Somalia,709,-75.1,350.3,,,-74.8,-14.2,50.3,6.4,-1.5,34 +4,3,Somalia,709,-61.8,79.5,,,-62.7,77.8,50.3,6.4,-1.5,111 +3,9.3,Somalia,709,-60,82,,,-60.9,80.4,50.3,6.4,-1.5,124.5 +3,5,Somalia,709,-46,40,-33.3,36.9,-33.7,35.2,50.3,6.4,-1.5,257 +5,1.9,Somalia,709,27.5,344.8,,,27.8,-15.9,50.3,6.4,-1.5,522 +4,5,Somalia,709,-28.4,319.1,,,-27.7,-42.4,50.3,6.4,-1.5,547 +4,6.7,Northwest_Africa,714,-87.5,358.2,,,,,,,,7.5 +3,5.2,Northwest_Africa,714,-77.8,326.2,,,,,,,,8 +4,4.1,Northwest_Africa,714,-81.9,294.4,,,,,,,,13 +3,2.3,Northwest_Africa,714,-86.8,202.9,,,,,,,,13 +5,8,Northwest_Africa,714,-72,71.2,,,,,,,,81 +6,6.3,Northwest_Africa,714,-65.2,20.3,,,-65.3,25.8,33.6,26,2.3,152.5 +3,19.2,Northwest_Africa,714,-62.5,61.6,,,-61.3,65.8,33.6,26,2.3,160 +5,7.4,Northwest_Africa,714,-68.5,62.4,,,-67.3,67.4,33.6,26,2.3,185.5 +4,4.1,Northwest_Africa,714,-69.4,52,,,-68.5,57.7,33.6,26,2.3,187 +4,6.1,Northwest_Africa,714,-71.4,60.2,,,-70.2,66,33.6,26,2.3,187 +4,6.2,Northwest_Africa,714,-82.9,32.7,,,-82.4,48.7,33.6,26,2.3,193 +5,4.1,Northwest_Africa,714,-73,64.7,,,-71.7,70.6,33.6,26,2.3,200 +6,2.6,Northwest_Africa,714,-70.9,55.1,-76.2,78.9,-74.6,84.4,33.6,26,2.3,206.5 +5,6,Northwest_Africa,714,-29,60,-26.8,56.5,-25.7,58.6,33.6,26,2.3,273 +5,3.6,Northwest_Africa,714,-29.1,57.8,-26.3,54.2,-25.4,56.3,33.6,26,2.3,275 +5,2.8,Northwest_Africa,714,-38.5,57.5,-33.7,52.4,-32.8,54.8,33.6,26,2.3,286.5 +4,4.1,Northwest_Africa,714,-33.8,61.4,-29,55.5,-28,57.7,33.6,26,2.3,290 +5,3.5,Northwest_Africa,714,-28.7,55.8,-20.9,48.3,-20.1,50.3,33.6,26,2.3,307 +6,4.6,Northwest_Africa,714,-28.3,58.9,-21.1,51.5,-20.2,53.4,33.6,26,2.3,309 +5,2.6,Northwest_Africa,714,-32.8,55.7,-27.8,49.9,-27,52.1,33.6,26,2.3,310 +4,4.5,Northwest_Africa,714,-28.2,55.5,-20.3,47.8,-19.6,49.8,33.6,26,2.3,317 +7,5.3,Northwest_Africa,714,-26.6,44.7,-17.5,36.2,-17.1,38.1,33.6,26,2.3,320 +5,5.9,Northwest_Africa,714,-18.8,31.2,,,-18.6,33.2,33.6,26,2.3,348 +5,3.7,Northwest_Africa,714,-21,19.9,,,-21.2,21.9,33.6,26,2.3,365 +6,3.7,Northwest_Africa,714,-19.2,19.8,,,-19.4,21.8,33.6,26,2.3,365 +5,6.6,Northwest_Africa,714,-43.4,8.6,,,-43.9,11.7,33.6,26,2.3,409 +5,2.5,Northeast_Africa,715,-87.6,346.9,,,,,,,,1.5 +5,3.1,Northeast_Africa,715,-84.9,127.9,,,,,,,,3 +5,7.4,Northeast_Africa,715,-86.1,336.5,,,,,,,,4 +5,11.2,Northeast_Africa,715,-69,4,,,,,,,,11.5 +5,8.3,Northeast_Africa,715,-78.4,16.1,,,,,,,,11.5 +5,2.6,Northeast_Africa,715,-73,0.5,-83.2,341.7,,,,,,13.5 +4,12.7,Northeast_Africa,715,-83,13.3,,,,,,,,26.5 +6,6,Northeast_Africa,715,-79.6,332.2,-79.5,258.6,,,,,,29 +5,5.4,Northeast_Africa,715,-77.9,32.8,,,,,,,,30 +5,6.4,Northeast_Africa,715,-83.5,318.6,,,,,,,,37 +6,4,Northeast_Africa,715,-71,340,-76.5,308.4,,,,,,37.5 +6,4.2,Northeast_Africa,715,-78.1,342.8,,,,,,,,42.5 +6,6,Northeast_Africa,715,-68,338,-73,314.8,,,,,,44.5 +3,5.8,Northeast_Africa,715,-69.4,9.4,,,,,,,,44.5 +5,8.5,Northeast_Africa,715,-69.3,78.1,,,-69.4,78.3,39.9,-61.4,-0.2,93 +3,18.1,Northeast_Africa,715,-75.7,48.3,,,-75.9,48.4,40.1,-61.4,-0.2,94.5 +3,11.5,Northeast_Africa,715,-54.9,43.3,-59.7,27,-60.2,26.5,40.5,-61.4,-0.7,221.5 +4,5.5,Northeast_Africa,715,-54.5,45.8,,,-55,45.6,40.5,-61.4,-0.7,231 +5,3.8,Northeast_Africa,715,-59.3,34.1,,,-59.8,33.7,40.5,-61.4,-0.7,231 +4,6,Northeast_Africa,715,-40.8,71.3,,,-41.2,71.2,40.4,-61.4,-0.7,280 +3,10.8,Northeast_Africa,715,25.9,11.6,,,25.4,11.2,40.4,-61.4,-0.7,377 +6,9.2,Northeast_Africa,715,39.6,329.5,,,39.3,-30.6,40.4,-61.4,-0.7,463 +4,5.4,East_Gondwana,501,-39.2,105.6,,,-73.4,66.6,18.8,22.6,-39.2,64 +6,5.7,East_Gondwana,501,-39,100.8,,,-72,54.5,19,21.9,-40.2,65 +4,6.7,East_Gondwana,501,-40,96,,,-70.6,42.6,19.2,21.5,-40.7,65.5 +6,5.9,East_Gondwana,501,-38.4,102.4,,,-72.7,57.9,19.2,21.5,-40.7,65.5 +3,3.8,East_Gondwana,501,-34.5,103.6,,,-69.8,66.8,19.2,21.5,-40.7,65.5 +5,9.4,East_Gondwana,501,-37.2,100.5,,,-70.8,56.2,19.2,21.5,-40.7,65.5 +3,3.8,East_Gondwana,501,-39,99,,,-71.5,49.9,19.2,21.5,-40.7,65.5 +4,10.1,East_Gondwana,501,-34.6,94,,,-67.6,42.9,20.2,19.3,-43.8,69 +4,12,East_Gondwana,501,-21.6,119.4,,,-74.9,73.7,19.8,27.2,-59.2,88 +5,7.5,East_Gondwana,501,-14.2,117.8,,,-66.7,79.4,20.2,27.6,-58.5,91.2 +3,4,East_Gondwana,501,-3,118,,,-49.8,86.8,23.4,31.3,-53.8,116 +6,3.5,East_Gondwana,501,-7,117,,,-53.3,83.2,23.4,31.3,-53.8,116 +6,8.3,East_Gondwana,501,-9.3,124.8,,,-57.6,95.2,23.4,31.3,-53.8,116 +4,5.5,East_Gondwana,501,-6.5,120.2,,,-53.8,88.5,23.4,31.3,-53.8,116 +3,7,East_Gondwana,501,-16,121,,,-63,84.4,23.5,31.4,-53.7,116.5 +5,2.4,East_Gondwana,501,-9.4,116.6,,,-55.1,81.1,23.5,31.5,-53.6,117 +5,4.6,East_Gondwana,501,-10.1,130.1,-2.4,118.5,-45.1,71.3,29.8,42.1,-60.5,206.5 +4,4.6,East_Gondwana,501,7.3,124.3,12.2,110.9,-28.7,73.5,29.8,42.1,-60.5,243 +5,6,East_Gondwana,501,7.5,120.5,13.4,109.1,-26.7,72.5,29.8,42.1,-60.5,248.5 +6,4.3,East_Gondwana,501,2.2,125.8,,,-44.6,83.5,29.8,42.1,-60.5,250.5 +3,1.8,East_Gondwana,501,4.1,102.8,10.2,94.5,-20.3,58.1,29.8,42.1,-60.5,250.5 +6,6.5,East_Gondwana,501,4,129,9.5,115.1,-33.2,75.8,29.8,42.1,-60.5,250.5 +3,12.1,East_Gondwana,501,18.1,111,,,-23.8,77,29.8,42.1,-60.5,289.5 +5,4.1,Arabia,503,-82.4,62.2,,,-82.9,57.8,36.5,18,-0.8,2.5 +5,3.4,Arabia,503,-76,13.5,,,-75.2,-3.1,36.5,18,-4.5,19 +5,3.6,Arabia,503,-74.2,69.1,,,-77.7,46.7,35.7,20.4,-7.1,29 +3,7.2,Arabia,503,-25.6,64,-13.4,58.4,-17.7,51.5,37.1,17.1,-8.9,306.5 +5,4.4,Madagascar,702,-74,43.7,,,,,,,,86.5 +4,7.6,Madagascar,702,-64,63,,,,,,,,87 +4,4.9,Madagascar,702,-66.1,49.7,,,,,,,,87 +4,4.4,Madagascar,702,-65.8,35.6,,,,,,,,87 +4,8.9,Madagascar,702,-73.7,73.1,,,,,,,,87 +5,4.3,Madagascar,702,-65.5,38,,,,,,,,87 +4,6.9,Madagascar,702,-70.3,63.1,,,,,,,,87 +5,2.4,Madagascar,702,-76.8,68.2,,,,,,,,87 +5,10.7,Madagascar,702,-66.7,43.5,,,,,,,,91 +3,5.9,Madagascar,702,-74,97.1,-65.2,70.8,-50.9,59.5,14.8,137.5,-15.4,206.5 +3,7.6,Madagascar,702,-76.7,110.8,-68.4,73.5,-54.3,60.4,14.8,137.5,-15.4,250.5 +3,9.5,Madagascar,702,-51.3,72.6,-42.5,60.2,-27.8,54.4,14.8,137.5,-15.4,305 +5,11,Madagascar,702,-6.8,1,,,3.1,-2.6,14.8,137.5,-15.4,509 +4,14,Madagascar,702,-6.8,352.7,,,1.6,-10.7,14.8,137.5,-15.4,521 +5,1.8,Australia,801,-78.4,103.1,-70.3,148.6,-80.2,155.3,12,48.8,-10.4,17.5 +5,3.6,Australia,801,-70.5,125.6,,,-80.2,19,14.1,57,-24.4,53 +5,1.8,Australia,801,-48.9,148.7,-45.1,146.8,-56.4,93.9,17.2,103.2,-36,112 +5,3.7,Australia,801,-53,179.6,,,-59.2,64.5,19.5,117.8,-56.2,168 +4,6,Australia,801,-46.1,175.2,,,-57.3,78.3,19.5,117.8,-56.2,180 +5,2.9,Australia,801,-50.7,174.5,,,-56.6,69.9,19.5,117.8,-56.2,183 +6,5.1,Australia,801,-63.8,124.5,,,-31.5,48.8,19.5,117.8,-56.2,321 +6,6,Australia,801,-47.1,41,-49.3,54.7,2.2,28.5,19.5,117.8,-56.2,367 +6,7.8,Australia,801,-49.1,38,,,3.6,17.7,19.5,117.8,-56.2,370 +6,15.2,Australia,801,-62,23.2,,,-10,10.7,19.5,117.8,-56.2,370 +5,5,Australia,801,-63,38.6,,,-10.3,18,19.5,117.8,-56.2,376 +5,3,Australia,801,-26.7,33.7,-30.4,46,21.9,25.2,19.5,117.8,-56.2,465 +4,13,Australia,801,-13,25,,,38.3,1.5,19.5,117.8,-56.2,493 +6,7.4,Australia,801,3.1,54.1,,,52.7,45,19.5,117.8,-56.2,495 +3,3.8,Australia,801,-37.5,34.4,,,15.2,14.8,19.5,117.8,-56.2,500 +6,10,Australia,801,-19.3,39.1,,,33.4,19,19.5,117.8,-56.2,510 +3,10.1,Australia,801,-31.4,26.9,-31.3,26.7,20.6,7.4,19.5,117.8,-56.2,510 +3,10.4,Australia,801,-38.3,24.5,,,13.4,6.9,19.5,117.8,-56.2,510 +3,12.3,Australia,801,-33.8,15.1,-32.8,13.9,16.4,-3.2,19.5,117.8,-56.2,515 +7,6.7,Australia,801,-43.2,339.9,,,-6.2,-20.8,19.5,117.8,-56.2,522 +6,14.4,Australia,801,-37.4,20.1,,,13.5,3.2,19.5,117.8,-56.2,523 +6,7.3,Australia,801,-32.7,11.5,-28.7,5.7,17.4,-11.7,19.5,117.8,-56.2,525 +6,4.1,Australia,801,-46.6,337.3,-37.4,331.7,-6.5,-29.4,19.5,117.8,-56.2,534 +5,11.4,Australia,801,-21.3,14.9,,,27.5,-6.6,19.5,117.8,-56.2,535 +6,16,Australia,801,-33,328,-20.2,326.5,1.6,-45.3,19.5,117.8,-56.2,550 +5,4,East_Antarctica,802,-85.5,143.6,,,-85.5,141.7,8.3,-49.4,0.2,1 +5,6.3,East_Antarctica,802,-87.3,137.3,,,-87.3,130.8,8.2,-49.4,0.3,2 +5,7,East_Antarctica,802,-86.4,28.5,,,-85.7,30.7,8.2,-49.4,0.8,5 +5,2.3,East_Antarctica,802,-85.5,9.3,,,-81.7,25.3,11.6,-48.2,4.2,27 +5,4.4,East_Antarctica,802,-51.4,203.4,,,-55,98.6,10.5,148.8,-58.2,164 +4,3.4,East_Antarctica,802,-47.8,225.5,,,-69.3,91.3,10.5,148.8,-58.2,183 +3,3.3,East_Antarctica,802,-57.8,224.3,,,-61.5,75.7,10.5,148.8,-58.2,183 +5,2.4,East_Antarctica,802,-45.3,208,,,-59.4,108.3,10.5,148.8,-58.2,183 +5,10.2,East_Antarctica,802,-50.5,211.4,,,-60,97.1,10.5,148.8,-58.2,183 +4,6.9,East_Antarctica,802,-44.1,231.5,,,-74.9,91.3,10.5,148.8,-58.2,193 +5,3.8,East_Antarctica,802,-41.8,226.5,,,-73.2,105.8,10.5,148.8,-58.2,195 +5,5.2,East_Antarctica,802,-2.5,23.8,,,37.8,-2.5,10.5,148.8,-58.2,471.5 +5,7.6,East_Antarctica,802,-3.5,22.7,,,36.3,-2.7,10.5,148.8,-58.2,475 +5,7.2,East_Antarctica,802,-11,21,,,29.6,2,10.5,148.8,-58.2,479 +4,10.9,East_Antarctica,802,-9.3,26.7,,,34.5,5.7,10.5,148.8,-58.2,484 +4,12,East_Antarctica,802,-7,21.4,,,32.8,-0.9,10.5,148.8,-58.2,499 +4,8.1,East_Antarctica,802,-5.4,18.5,,,32,-4.6,10.5,148.8,-58.2,500 +3,4.5,East_Antarctica,802,-28.5,9.5,,,9.5,6.1,10.5,148.8,-58.2,515 diff --git a/examples/Torsvik_2012/plots/plot_results.jl b/examples/Torsvik_2012/plots/plot_results.jl new file mode 100644 index 0000000..eaa314a --- /dev/null +++ b/examples/Torsvik_2012/plots/plot_results.jl @@ -0,0 +1,140 @@ +using Pkg; Pkg.activate(".") + +using SphereUDE +using Serialization +using Plots +using Plots.PlotMeasures +using JLD2 +using LinearAlgebra + +JLD2.@load "examples/Torsvik_2012/results/results_dict.jld2" + +data_directions_sph = cart2sph(results_dict["directions"], radians=false) +fit_directions_sph = cart2sph(results_dict["fit_directions"], radians=false) + +data_lat = data_directions_sph[1,:] +fit_lat = fit_directions_sph[1,:] + +data_lon = data_directions_sph[2,:] +fit_lon = fit_directions_sph[2,:] + +α_error = 140 ./ sqrt.(results_dict["kappas"]) + +# Latitude Plot + +Plots.scatter(results_dict["times"], data_lat, label="Paleopole Latitudes", yerr=α_error, c=:lightsteelblue2, markersize=5) +plot_latitude = Plots.plot!(results_dict["fit_times"], fit_lat, label="Estimated APWP using SphereUDE", + xlabel="Age (Myr)", yticks=[-90, -60, -30, 0, 30, 60], + ylabel="Latitude (degrees)", ylims=(-90,60), lw = 4, c=:brown, + legend=:topleft) +plot!(fontfamily="Computer Modern", + #title="PIL51", + titlefontsize=18, + tickfontsize=15, + legendfontsize=15, + guidefontsize=18, + #ylimits=(0.1,10), + #xlimits=(10^(-4),10^(-1)), + legend=true, + margin= 7mm, + size=(1200,500), + dpi=600) + +Plots.savefig(plot_latitude, "examples/Torsvik_2012/plots/latitude.pdf") + +### Longitude Plot + +α_error_lon = α_error ./ cos.(π ./ 180. .* data_lat) + +Plots.scatter(results_dict["times"], data_lon, label="Paleopole Longitudes", yerr=α_error_lon, c=:lightsteelblue2, markersize=5) +plot_longitude = Plots.plot!(results_dict["fit_times"], fit_lon, label="Estimated APWP using SphereUDE", + xlabel="Age (Myr)", yticks=[-180, -120, -60 , 0, 60, 120, 180], + ylabel="Longitude (degrees)", ylims=(-180,180), lw = 4, c=:brown, + legend=:bottomright) +plot!(fontfamily="Computer Modern", + #title="PIL51", + titlefontsize=18, + tickfontsize=15, + legendfontsize=15, + guidefontsize=18, + #ylimits=(0.1,10), + #xlimits=(10^(-4),10^(-1)), + margin= 7mm, + size=(1200,500), + dpi=600) + +Plots.savefig(plot_longitude, "examples/Torsvik_2012/plots/longitude.pdf") + +### Lat and long combined + +combo_plot = plot(plot_latitude, plot_longitude, layout = (2, 1)) +plot!(fontfamily="Computer Modern", + #title="PIL51", + titlefontsize=18, + tickfontsize=15, + legendfontsize=15, + guidefontsize=18, + #ylimits=(0.1,10), + #xlimits=(10^(-4),10^(-1)), + margin= 7mm, + size=(1200,800), + dpi=600) +Plots.savefig(combo_plot, "examples/Torsvik_2012/plots/latitude_longitude.pdf") + + +### Angular velocity Plot + +angular_velocity = mapslices(x -> norm(x), results_dict["fit_rotations"], dims=1)[:] +angular_velocity_path = [norm(cross(results_dict["fit_directions"][:,i], results_dict["fit_rotations"][:,i] )) for i in axes(results_dict["fit_rotations"],2)] + +plot_angular = Plots.plot(results_dict["fit_times"], angular_velocity, label="Maximum total angular velocity", + xlabel="Age (Myr)", + ylabel="Angular velocity (degrees/My)", lw = 5, c=:darkgreen, + legend=:topleft) +plot!(results_dict["fit_times"], angular_velocity_path, label="Pole angular velocity", lw = 4, c=:lightseagreen, ls=:dot) +plot!(fontfamily="Computer Modern", + #title="PIL51", + titlefontsize=18, + tickfontsize=15, + legendfontsize=15, + guidefontsize=18, + #ylimits=(0.1,10), + #xlimits=(10^(-4),10^(-1)), + margin= 10mm, + size=(1200,500), + dpi=300) + + +Plots.savefig(plot_angular, "examples/Torsvik_2012/plots/angular.pdf") + + +### Loss function + +losses = results_dict["losses"] + +plot_loss = Plots.plot(1:length(losses), losses, label="Loss Function", + xlabel="Epoch", + ylabel="Loss", lw = 5, c=:indigo, + yaxis=:log, + # yticks=[1,10,100], + xlimits=(0,10000), + ylimits=(10, 1000), + xticks=[0,2500,5000,7500,10000], + legend=:topright) + +vspan!(plot_loss, [0,5000], color = :navajowhite3, alpha = 0.2, labels = "ADAM"); +vspan!(plot_loss, [5000,10000], color = :navajowhite4, alpha = 0.2, labels = "BFGS"); + +plot!(fontfamily="Computer Modern", + #title="PIL51", + titlefontsize=18, + tickfontsize=15, + legendfontsize=15, + guidefontsize=18, + #ylimits=(0.1,10), + #xlimits=(10^(-4),10^(-1)), + margin= 10mm, + size=(700,500), + dpi=300) + +Plots.savefig(plot_loss, "examples/Torsvik_2012/plots/loss.pdf") diff --git a/src/SphereUDE.jl b/src/SphereUDE.jl index 6003420..23194ee 100644 --- a/src/SphereUDE.jl +++ b/src/SphereUDE.jl @@ -20,6 +20,7 @@ using Infiltrator include("types.jl") include("utils.jl") +include("losses.jl") include("train.jl") include("plot.jl") diff --git a/src/losses.jl b/src/losses.jl new file mode 100644 index 0000000..c48083e --- /dev/null +++ b/src/losses.jl @@ -0,0 +1,214 @@ +export loss +export loss_empirical +export regularization +export loss_multiple_shooting + +""" +General Loss Function +""" +function loss(β::ComponentVector, + data::AD, + params::AP, + # predict::Function, + # predict_L::Function, + U::Chain, + st::NamedTuple) where {AD <: AbstractData, AP <: AbstractParameters} + + # Record the value of each individual loss to the total loss function for hyperparameter selection. + loss_dict = Dict() + + l_emp = loss_empirical(β, data, params) + + loss_dict["Empirical"] = l_emp + + # Regularization + l_reg = 0.0 + if !isnothing(params.reg) + for reg in params.reg + reg₀ = regularization(β.θ, U, st, reg, params) + l_reg += reg₀ + loss_dict["Regularization (order=$(reg.order), power=$(reg.power))"] = reg₀ + end + end + return l_emp + l_reg, loss_dict +end + +""" +Empirical loss function +""" +function loss_empirical(β::ComponentVector, + data::AD, + params::AP) where {AD <: AbstractData, AP <: AbstractParameters} + + # Predict trajectory on times associated to dataset + u_, retcode = predict(β, params, data.times) + + # If numerical integration fails or bad choice of parameter, return infinity + if retcode != :Success + @warn "[SphereUDE] Numerical solver not converging. This can be causes by numerical innestabilities around a bad choice of parameter. This can be due to just a bad initial condition of the neural network, so it is worth changing the randon number used for initialization. " + return Inf + end + + # Empirical error + if isnothing(data.kappas) + # The 3 is needed since the mean is computen on a 3xN matrix + return 3.0 * mean(abs2.(u_ .- data.directions)) + # l_emp = 1 - 3.0 * mean(u_ .* data.directions) + else + return mean(data.kappas .* sum(abs2.(u_ .- data.directions), dims=1)) + # l_emp = norm(data.kappas)^2 - 3.0 * mean(data.kappas .* u_ .* data.directions) + end + +end + +""" +Empirical Prediction function +""" +function predict(β::ComponentVector, + # U::Chain, + # st::NamedTuple, + params::AP, + T::Vector) where {AP <: AbstractParameters} + + if params.train_initial_condition + _prob = remake(prob_nn, u0=β.u0 / norm(β.u0), # We enforce the norm=1 condition again here + tspan=(min(T[1], params.tmin), max(T[end], params.tmax)), + p = β.θ) + else + _prob = remake(prob_nn, u0=params.u0, + tspan=(min(T[1], params.tmin), max(T[end], params.tmax)), + p = β.θ) + end + sol = solve(_prob, params.solver, saveat=T, + abstol=params.abstol, reltol=params.reltol, + sensealg=params.sensealg, + dtmin=1e-4 * (params.tmax - params.tmin), force_dtmin=true) # Force minimum step in case L(t) changes drastically due to bad behaviour of neural network + return Array(sol), sol.retcode +end + +""" +Regularization +""" +function regularization(θ::ComponentVector, + U::Chain, + st::NamedTuple, + reg::AG, + params::AP) where {AG <: AbstractRegularization, AP <: AbstractParameters} + + # Define Statefull Lux NN + smodel = StatefulLuxLayer{true}(U, θ, st) + + l_ = 0.0 + if reg.order==0 + + # l_ += quadrature(t -> norm(predict_L(t, U, θ, st))^reg.power, params.tmin, params.tmax, params.n_quadrature) + l_ += quadrature(t -> norm(smodel([t]))^reg.power, params.tmin, params.tmax, params.n_quadrature) + + elseif reg.order==1 + + if typeof(reg.diff_mode) <: LuxNestedAD + # Automatic Differentiation + nodes, weights = quadrature(params.tmin, params.tmax, params.n_quadrature) + + if reg.diff_mode.method == "ForwardDiff" + # Jac = ForwardDiff.jacobian(smodel, reshape(nodes, 1, params.n_quadrature)) + Jac = batched_jacobian(smodel, AutoForwardDiff(), reshape(nodes, 1, params.n_quadrature)) + elseif reg.diff_mode.method == "Zygote" + # This can also be done with Zygote in reverse mode + # Jac = Zygote.jacobian(smodel, reshape(nodes, 1, params.n_quadrature))[1] + Jac = batched_jacobian(smodel, AutoZygote(), reshape(nodes, 1, params.n_quadrature)) + else + throw("Method for AD backend no implemented.") + end + + # Compute the final agregation to the loss + l_ += sum([weights[j] * norm(Jac[:,1,j])^reg.power for j in 1:params.n_quadrature]) + + # Test every a few iterations that AD is working properly + if rand(Bernoulli(0.001)) + l_AD = sum([weights[j] * norm(Jac[:,1,j])^reg.power for j in 1:params.n_quadrature]) + l_FD = quadrature(t -> norm(central_fdm(τ -> smodel([τ]), t, 1e-5))^reg.power, params.tmin, params.tmax, params.n_quadrature) + if abs(l_AD - l_FD) < 1e-2 * l_FD + @warn "[SphereUDE] Nested AD is giving significant different results than Finite Differences." + @printf "[SphereUDE] Regularization with AD: %.9f vs %.9f using Finite Differences" l_AD l_FD + end + end + + elseif typeof(reg.diff_mode) <: FiniteDifferences + # Finite differences + # l_ += quadrature(t -> norm(central_fdm(τ -> predict_L(τ, U, θ, st), t, reg.diff_mode.ϵ))^reg.power, params.tmin, params.tmax, params.n_quadrature) + l_ += quadrature(t -> norm(central_fdm(τ -> smodel([τ]), t, reg.diff_mode.ϵ))^reg.power, params.tmin, params.tmax, params.n_quadrature) + + elseif typeof(reg.diff_mode) <: ComplexStepDifferentiation + # Complex step differentiation + # l_ += quadrature(t -> norm(complex_step_differentiation(τ -> predict_L(τ, U, θ, st), t, reg.diff_mode.ϵ))^reg.power, params.tmin, params.tmax, params.n_quadrature) + l_ += quadrature(t -> norm(complex_step_differentiation(τ -> smodel([τ]), t, reg.diff_mode.ϵ))^reg.power, params.tmin, params.tmax, params.n_quadrature) + + else + throw("Method not implemented.") + end + + else + throw("Method not implemented.") + end + + return reg.λ * l_ +end + + +""" +Loss function to be called for multiple shooting +This seems duplicated from before, so be careful with this +""" +function loss_multiple_shooting(data, pred) + + throw("[SphereUDE] Method to be implemented.") + # # Empirical error + # l_emp = 3.0 * mean(abs2.(pred .- data)) + + # # Regularization + # l_reg = 0.0 + # if !isnothing(params.reg) + # for reg in params.reg + # reg₀ = regularization(β.θ, reg) + # l_reg += reg₀ + # end + # end + + # return l_emp + l_reg + # end + + # # Define parameters for Multiple Shooting + # group_size = 10 + # continuity_term = 100 + + # ps = ComponentArray(θ) # are these necesary? + # pd, pax = getdata(ps), getaxes(ps) + + # function continuity_loss(u_pred, u_initial) + # if !isapprox(norm(u_initial), 1.0, atol=1e-6) || !isapprox(norm(u_pred), 1.0, atol=1e-6) + # @warn "Directions during multiple shooting are not in the sphere. Small deviations from unit norm observed:" + # @show norm(u_initial), norm(u_pred) + # end + # return sum(abs2, u_pred - u_initial) + # end + + # function loss_multiple_shooting(β::ComponentVector) + + # ps = ComponentArray(β.θ, pax) + + # if params.train_initial_condition + # _prob = remake(prob_nn, u0=β.u0 / norm(β.u0), # We enforce the norm=1 condition again here + # tspan=(min(data.times[1], params.tmin), max(data.times[end], params.tmax)), + # p = β.θ) + # else + # _prob = remake(prob_nn, u0=params.u0, + # tspan=(min(data.times[1], params.tmin), max(data.times[end], params.tmax)), + # p = β.θ) + # end + + # return multiple_shoot(β.θ, data.directions, data.times, _prob, _loss_multiple_shooting, continuity_loss, params.solver, + # group_size; continuity_term, sensealg=params.sensealg) + + +end \ No newline at end of file diff --git a/src/train.jl b/src/train.jl index 23f120e..5d4bc55 100644 --- a/src/train.jl +++ b/src/train.jl @@ -1,42 +1,5 @@ export train -function get_NN(params, rng, θ_trained) - # Define default neural network - - # For L1 regularization relu_cap works better, but for L2 I think is better to include sigmoid - if isL1reg(params.reg) - @warn "[SphereUDE] Using ReLU activation functions for neural network due to L1 regularization." - U = Lux.Chain( - Lux.Dense(1, 5, sigmoid), - Lux.Dense(5, 10, sigmoid), - Lux.Dense(10, 10, sigmoid), - Lux.Dense(10, 10, sigmoid), - Lux.Dense(10, 5, sigmoid), - Lux.Dense(5, 3, Base.Fix2(sigmoid_cap, params.ωmax)) - ) - else - U = Lux.Chain( - Lux.Dense(1, 5, gelu), - Lux.Dense(5, 10, gelu), - Lux.Dense(10, 10, gelu), - Lux.Dense(10, 10, gelu), - Lux.Dense(10, 5, gelu), - Lux.Dense(5, 3, Base.Fix2(sigmoid_cap, params.ωmax)) - ) - end - return U -end - -""" - predict_L(t, NN, θ, st) - -Predict value of rotation given by L given by the neural network. -""" -function predict_L(t, NN, θ, st) - return NN([t], θ, st)[1] - # return 0.2 * ( NN([t-1.0], θ, st)[1] + NN([t-0.5], θ, st)[1] + NN([t], θ, st)[1] + NN([t+0.5], θ, st)[1] + NN([t+1.0], θ, st)[1] ) -end - """ train() @@ -51,14 +14,13 @@ function train(data::AD, # Raise warnings raise_warnings(data::AD, params::AP) + # Set Neural Network if isnothing(model) - U = get_NN(params, rng, θ_trained) + U = get_default_NN(params, rng, θ_trained) else U = model end θ, st = Lux.setup(rng, U) - # Make it a stateful layer (this I don't know where is best to add it, it repeats) - smodel = StatefulLuxLayer{true}(U, θ, st) # Set component vector for Optimization if params.train_initial_condition @@ -67,202 +29,28 @@ function train(data::AD, β = ComponentVector{Float64}(θ=θ) end + # Function to predict angular momentum. Used just for ODE part + # of the loss function + function predict_L(t, NN, θ, st) + smodel = StatefulLuxLayer{true}(NN, θ, st) + return smodel([t]) + end + + # Sphere-constrained ODE function ude_rotation!(du, u, p, t) # Angular momentum given by network prediction L = predict_L(t, U, p, st) du .= cross(L, u) end - prob_nn = ODEProblem(ude_rotation!, params.u0, [params.tmin, params.tmax], β.θ) - - function predict(β::ComponentVector; T=data.times) - if params.train_initial_condition - _prob = remake(prob_nn, u0=β.u0 / norm(β.u0), # We enforce the norm=1 condition again here - tspan=(min(T[1], params.tmin), max(T[end], params.tmax)), - p = β.θ) - else - _prob = remake(prob_nn, u0=params.u0, - tspan=(min(T[1], params.tmin), max(T[end], params.tmax)), - p = β.θ) - end - sol = solve(_prob, params.solver, saveat=T, - abstol=params.abstol, reltol=params.reltol, - sensealg=params.sensealg, - dtmin=1e-4 * (params.tmax - params.tmin), force_dtmin=true) # Force minimum step in case L(t) changes drastically due to bad behaviour of neural network - return Array(sol), sol.retcode - end - - ##### Definition of loss functions to be used ##### - - """ - General Loss Function - """ - function loss(β::ComponentVector) - - # Record the value of each individual loss to the total loss function for hyperparameter selection. - loss_dict = Dict() - - l_emp = loss_empirical(β) - - loss_dict["Empirical"] = l_emp - - # Regularization - l_reg = 0.0 - if !isnothing(params.reg) - for reg in params.reg - reg₀ = regularization(β.θ, reg) - l_reg += reg₀ - loss_dict["Regularization (order=$(reg.order), power=$(reg.power))"] = reg₀ - end - end - return l_emp + l_reg, loss_dict - end - - """ - Empirical loss function - """ - function loss_empirical(β::ComponentVector) - - u_, retcode = predict(β) - - # If numerical integration fails or bad choice of parameter, return infinity - if retcode != :Success - @warn "[SphereUDE] Numerical solver not converging. This can be causes by numerical innestabilities around a bad choice of parameter. This can be due to just a bad initial condition of the neural network, so it is worth changing the randon number used for initialization. " - return Inf - end - - # Empirical error - if isnothing(data.kappas) - # The 3 is needed since the mean is computen on a 3xN matrix - return 3.0 * mean(abs2.(u_ .- data.directions)) - # l_emp = 1 - 3.0 * mean(u_ .* data.directions) - else - return mean(data.kappas .* sum(abs2.(u_ .- data.directions), dims=1)) - # l_emp = norm(data.kappas)^2 - 3.0 * mean(data.kappas .* u_ .* data.directions) - end + global prob_nn = ODEProblem(ude_rotation!, params.u0, [params.tmin, params.tmax], β.θ) - end - - """ - Regularization - """ - function regularization(θ::ComponentVector, reg::AG; n_nodes=100) where {AG <: AbstractRegularization} - - l_ = 0.0 - if reg.order==0 - - l_ += quadrature(t -> norm(predict_L(t, U, θ, st))^reg.power, params.tmin, params.tmax, n_nodes) - - elseif reg.order==1 - - if typeof(reg.diff_mode) <: LuxNestedAD - # Automatic Differentiation - nodes, weights = quadrature(params.tmin, params.tmax, n_nodes) - - if reg.diff_mode.method == "ForwardDiff" - # Jac = ForwardDiff.jacobian(smodel, reshape(nodes, 1, n_nodes)) - Jac = batched_jacobian(smodel, AutoForwardDiff(), reshape(nodes, 1, n_nodes)) - elseif reg.diff_mode.method == "Zygote" - # This can also be done with Zygote in reverse mode - # Jac = Zygote.jacobian(smodel, reshape(nodes, 1, n_nodes))[1] - Jac = batched_jacobian(smodel, AutoZygote(), reshape(nodes, 1, n_nodes)) - else - throw("Method for AD backend no implemented.") - end - - # Compute the final agregation to the loss - l_ += sum([weights[j] * norm(Jac[:,1,j])^reg.power for j in 1:n_nodes]) - - # Test every a few iterations that AD is working properly - if rand(Bernoulli(0.001)) - l_AD = sum([weights[j] * norm(Jac[:,1,j])^reg.power for j in 1:n_nodes]) - l_FD = quadrature(t -> norm(central_fdm(τ -> predict_L(τ, U, θ, st), t, 1e-5))^reg.power, params.tmin, params.tmax, n_nodes) - if abs(l_AD - l_FD) < 1e-2 * l_FD - @warn "[SphereUDE] Nested AD is giving significant different results than Finite Differences." - @printf "[SphereUDE] Regularization with AD: %.9f vs %.9f using Finite Differences" l_AD l_FD - end - end - - elseif typeof(reg.diff_mode) <: FiniteDifferences - # Finite differences - l_ += quadrature(t -> norm(central_fdm(τ -> predict_L(τ, U, θ, st), t, reg.diff_mode.ϵ))^reg.power, params.tmin, params.tmax, n_nodes) - - elseif typeof(reg.diff_mode) <: ComplexStepDifferentiation - # Complex step differentiation - l_ += quadrature(t -> norm(complex_step_differentiation(τ -> predict_L(τ, U, θ, st), t, reg.diff_mode.ϵ))^reg.power, params.tmin, params.tmax, n_nodes) - - else - throw("Method not implemented.") - end - - else - throw("Method not implemented.") - end - return reg.λ * l_ - end - - """ - Loss function to be called for multiple shooting - This seems duplicated from before, so be careful with this - """ - function _loss_multiple_shooting(data, pred) - - # Empirical error - l_emp = 3.0 * mean(abs2.(pred .- data)) - - # Regularization - l_reg = 0.0 - if !isnothing(params.reg) - for reg in params.reg - reg₀ = regularization(β.θ, reg) - l_reg += reg₀ - end - end - - return l_emp + l_reg - end - - # Define parameters for Multiple Shooting - group_size = 10 - continuity_term = 100 - - ps = ComponentArray(θ) # are these necesary? - pd, pax = getdata(ps), getaxes(ps) - - function continuity_loss(u_pred, u_initial) - if !isapprox(norm(u_initial), 1.0, atol=1e-6) || !isapprox(norm(u_pred), 1.0, atol=1e-6) - @warn "Directions during multiple shooting are not in the sphere. Small deviations from unit norm observed:" - @show norm(u_initial), norm(u_pred) - end - return sum(abs2, u_pred - u_initial) - end - - function loss_multiple_shooting(β::ComponentVector) - - ps = ComponentArray(β.θ, pax) - - if params.train_initial_condition - _prob = remake(prob_nn, u0=β.u0 / norm(β.u0), # We enforce the norm=1 condition again here - tspan=(min(data.times[1], params.tmin), max(data.times[end], params.tmax)), - p = β.θ) - else - _prob = remake(prob_nn, u0=params.u0, - tspan=(min(data.times[1], params.tmin), max(data.times[end], params.tmax)), - p = β.θ) - end - - return multiple_shoot(β.θ, data.directions, data.times, _prob, _loss_multiple_shooting, continuity_loss, params.solver, - group_size; continuity_term, sensealg=params.sensealg) - end - - - ### Callback function + ### Callback losses = Float64[] callback = function (p, l) push!(losses, l) if length(losses) % 50 == 0 @printf "Iteration: [%5d / %5d] \t Loss: %.9f \n" length(losses) (params.niter_ADAM+params.niter_LBFGS) losses[end] - # println("Current loss after $(length(losses)) iterations: $(losses[end])") end if params.train_initial_condition p.u0 ./= norm(p.u0) @@ -270,13 +58,15 @@ function train(data::AD, return false end - println("Loss after initalization: ", loss(β)[1]) - # Dispatch the right loss function - f_loss = params.multiple_shooting ? loss_multiple_shooting : loss - - # Optimization setting with AD - adtype = Optimization.AutoZygote() + if params.multiple_shooting + throw("[SphereUDE] Method not implemented.") + # f_loss(β) = loss_multiple_shooting + else + f_loss(β) = loss(β, data, params, U, st) + end + + println("Loss after initalization: ", f_loss(β)[1]) """ Pretraining to find parameters without impossing regularization @@ -291,7 +81,7 @@ function train(data::AD, end return false end - optf₀ = Optimization.OptimizationFunction((x, β) -> loss_empirical(x), adtype) + optf₀ = Optimization.OptimizationFunction((x, β) -> loss_empirical(x), params.adtype) optprob₀ = Optimization.OptimizationProblem(optf₀, β) res₀ = Optimization.solve(optprob₀, ADAM(), callback=callback_pretrain, maxiters=params.niter_ADAM, verbose=false) optprob₁ = Optimization.OptimizationProblem(optf₀, res₀.u) @@ -302,7 +92,7 @@ function train(data::AD, # To do: implement this with polyoptimizaion to put ADAM and BFGS in one step. # Maybe better to keep like this for the line search. - optf = Optimization.OptimizationFunction((x, β) -> (first ∘ f_loss)(x), adtype) + optf = Optimization.OptimizationFunction((x, β) -> (first ∘ f_loss)(x), params.adtype) optprob = Optimization.OptimizationProblem(optf, β) res1 = Optimization.solve(optprob, ADAM(), callback=callback, maxiters=params.niter_ADAM, verbose=true) @@ -330,14 +120,14 @@ function train(data::AD, # Final Fit fit_times = collect(range(params.tmin,params.tmax, length=1000)) - fit_directions, _ = predict(β_trained, T=fit_times) - fit_rotations = reduce(hcat, (t -> U([t], θ_trained, st)[1]).(fit_times)) + fit_directions, _ = predict(β_trained, params, fit_times) + fit_rotations = reduce(hcat, (t -> predict_L(t, U, β_trained.θ, st)).(fit_times)) # Recover final balance between different terms involved in the loss function to assess hyperparameter selection. - _, loss_dict = loss(β_trained) + _, loss_dict = f_loss(β_trained) pretty_table(loss_dict, sortkeys=true, header=["Loss term", "Value"]) return Results(θ=θ_trained, u0=u0_trained, U=U, st=st, fit_times=fit_times, fit_directions=fit_directions, fit_rotations=fit_rotations, losses=losses) -end +end \ No newline at end of file diff --git a/src/types.jl b/src/types.jl index 5d5ed8c..3278ee7 100644 --- a/src/types.jl +++ b/src/types.jl @@ -25,7 +25,9 @@ Training parameters niter_LBFGS::I reltol::F abstol::F + n_quadrature::I = 100 solver::OrdinaryDiffEqCore.OrdinaryDiffEqAlgorithm = Tsit5() + adtype::Optimization.AbstractADType = AutoZygote() sensealg::SciMLBase.AbstractAdjointSensitivityAlgorithm = QuadratureAdjoint(autojacvec=ReverseDiffVJP(true)) pretrain::Bool = false end diff --git a/src/utils.jl b/src/utils.jl index 47b86bc..7258a36 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -7,9 +7,11 @@ export quadrature, central_fdm, complex_step_differentiation export raise_warnings export isL1reg export convert2dict +export get_default_NN # Import activation function for complex extension import Lux: relu, gelu +using Lux: Chain # import Lux: sigmoid, relu, gelu ### Custom Activation Funtions @@ -109,6 +111,38 @@ function gelu(z::Complex) return 0.5 * z * (1 + tanh((sqrt(2/π))*(z + 0.044715 * z^3))) end +### Machine Learning Utils + +""" + +Return a default neural network for those cases where NN has not being provided by user. +""" +function get_default_NN(params::AP, rng, θ_trained) where {AP <: AbstractParameters} + # Define default neural network + + # For L1 regularization relu_cap works better, but for L2 I think is better to include sigmoid + if isL1reg(params.reg) + @warn "[SphereUDE] Using ReLU activation functions for neural network due to L1 regularization." + U = Lux.Chain( + Lux.Dense(1, 5, sigmoid), + Lux.Dense(5, 10, sigmoid), + Lux.Dense(10, 10, sigmoid), + Lux.Dense(10, 10, sigmoid), + Lux.Dense(10, 5, sigmoid), + Lux.Dense(5, 3, Base.Fix2(sigmoid_cap, params.ωmax)) + ) + else + U = Lux.Chain( + Lux.Dense(1, 5, gelu), + Lux.Dense(5, 10, gelu), + Lux.Dense(10, 10, gelu), + Lux.Dense(10, 10, gelu), + Lux.Dense(10, 5, gelu), + Lux.Dense(5, 3, Base.Fix2(sigmoid_cap, params.ωmax)) + ) + end + return U +end ### Spherical Utils @@ -178,15 +212,16 @@ end Numerical integral using Gaussian quadrature """ -function quadrature(f::Function, t₀, t₁, n_nodes::Int) - nodes, weigths = quadrature(t₀, t₁, n_nodes) +function quadrature(f::Function, t₀, t₁, n_quadrature::Int) + nodes, weigths = quadrature(t₀, t₁, n_quadrature) return dot(weigths, f.(nodes)) end -function quadrature(t₀, t₁, n_nodes::Int) +function quadrature(t₀, t₁, n_quadrature::Int) ignore() do # Ignore AD here since FastGaussQuadrature is using mutating arrays - nodes, weigths = gausslegendre(n_nodes) + nodes, weigths = gausslegendre(n_quadrature) + # nodes .+ rand(sampler(Uniform(-0.1,0.1)), n_quadrature) end nodes = (t₀+t₁)/2 .+ nodes * (t₁-t₀)/2 weigths = (t₁-t₀) / 2 * weigths diff --git a/test/rotation.jl b/test/rotation.jl index 43b8b55..bb3cdd1 100644 --- a/test/rotation.jl +++ b/test/rotation.jl @@ -47,7 +47,7 @@ function test_single_rotation() niter_ADAM = 20, niter_LBFGS = 10, sensealg = GaussAdjoint(autojacvec = ReverseDiffVJP(true))) - results = train(data, params, rng, nothing) + results = train(data, params, rng, nothing, nothing) @test true From f22f8096af0aa5a6b9d237c5c4b705f650845190 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Thu, 17 Oct 2024 17:21:14 -0700 Subject: [PATCH 19/22] Implementation of cubic splines as Jupp 1987 --- src/losses.jl | 102 +++++++++++++++++++++++++++++++++++++------------- src/train.jl | 2 +- src/types.jl | 6 +++ 3 files changed, 84 insertions(+), 26 deletions(-) diff --git a/src/losses.jl b/src/losses.jl index c48083e..127164f 100644 --- a/src/losses.jl +++ b/src/losses.jl @@ -1,6 +1,7 @@ export loss export loss_empirical export regularization +export cubic_regularization export loss_multiple_shooting """ @@ -9,8 +10,6 @@ General Loss Function function loss(β::ComponentVector, data::AD, params::AP, - # predict::Function, - # predict_L::Function, U::Chain, st::NamedTuple) where {AD <: AbstractData, AP <: AbstractParameters} @@ -24,10 +23,18 @@ function loss(β::ComponentVector, # Regularization l_reg = 0.0 if !isnothing(params.reg) - for reg in params.reg - reg₀ = regularization(β.θ, U, st, reg, params) - l_reg += reg₀ - loss_dict["Regularization (order=$(reg.order), power=$(reg.power))"] = reg₀ + for reg in params.reg + if typeof(reg) <: Regularization + reg₀ = regularization(β.θ, U, st, reg, params) + loss_dict["Regularization (order=$(reg.order), power=$(reg.power))"] = reg₀ + elseif typeof(reg) <: CubicSplinesRegularization + reg₀ = cubic_regularization(β, U, st, reg, params) + loss_dict["Regularization with cubic splines"] = reg₀ + else + throw("Regularization not implemented.") + end + # Add contribution to regularization + l_reg += reg₀ end end return l_emp + l_reg, loss_dict @@ -41,13 +48,7 @@ function loss_empirical(β::ComponentVector, params::AP) where {AD <: AbstractData, AP <: AbstractParameters} # Predict trajectory on times associated to dataset - u_, retcode = predict(β, params, data.times) - - # If numerical integration fails or bad choice of parameter, return infinity - if retcode != :Success - @warn "[SphereUDE] Numerical solver not converging. This can be causes by numerical innestabilities around a bad choice of parameter. This can be due to just a bad initial condition of the neural network, so it is worth changing the randon number used for initialization. " - return Inf - end + u_ = predict(β, params, data.times) # Empirical error if isnothing(data.kappas) @@ -65,8 +66,6 @@ end Empirical Prediction function """ function predict(β::ComponentVector, - # U::Chain, - # st::NamedTuple, params::AP, T::Vector) where {AP <: AbstractParameters} @@ -79,11 +78,20 @@ function predict(β::ComponentVector, tspan=(min(T[1], params.tmin), max(T[end], params.tmax)), p = β.θ) end + + # Force minimum step in case L(t) changes drastically due to bad behaviour of neural network sol = solve(_prob, params.solver, saveat=T, abstol=params.abstol, reltol=params.reltol, sensealg=params.sensealg, - dtmin=1e-4 * (params.tmax - params.tmin), force_dtmin=true) # Force minimum step in case L(t) changes drastically due to bad behaviour of neural network - return Array(sol), sol.retcode + dtmin=1e-4 * (params.tmax - params.tmin), force_dtmin=true) + + # If numerical integration fails or bad choice of parameter, return infinity + if sol.retcode != :Success + @warn "[SphereUDE] Numerical solver not converging. This can be causes by numerical innestabilities around a bad choice of parameter. This can be due to just a bad initial condition of the neural network, so it is worth changing the randon number used for initialization. " + return Inf + end + + return Array(sol) end """ @@ -92,8 +100,8 @@ Regularization function regularization(θ::ComponentVector, U::Chain, st::NamedTuple, - reg::AG, - params::AP) where {AG <: AbstractRegularization, AP <: AbstractParameters} + reg::Regularization, + params::AP) where {AP <: AbstractParameters} # Define Statefull Lux NN smodel = StatefulLuxLayer{true}(U, θ, st) @@ -125,12 +133,14 @@ function regularization(θ::ComponentVector, l_ += sum([weights[j] * norm(Jac[:,1,j])^reg.power for j in 1:params.n_quadrature]) # Test every a few iterations that AD is working properly - if rand(Bernoulli(0.001)) - l_AD = sum([weights[j] * norm(Jac[:,1,j])^reg.power for j in 1:params.n_quadrature]) - l_FD = quadrature(t -> norm(central_fdm(τ -> smodel([τ]), t, 1e-5))^reg.power, params.tmin, params.tmax, params.n_quadrature) - if abs(l_AD - l_FD) < 1e-2 * l_FD - @warn "[SphereUDE] Nested AD is giving significant different results than Finite Differences." - @printf "[SphereUDE] Regularization with AD: %.9f vs %.9f using Finite Differences" l_AD l_FD + ignore() do + if rand(Bernoulli(0.001)) + l_AD = sum([weights[j] * norm(Jac[:,1,j])^reg.power for j in 1:params.n_quadrature]) + l_FD = quadrature(t -> norm(central_fdm(τ -> smodel([τ]), t, 1e-5))^reg.power, params.tmin, params.tmax, params.n_quadrature) + if abs(l_AD - l_FD) > 1e-2 * abs(l_FD) + @warn "[SphereUDE] Nested AD is giving significant different results than Finite Differences." + @printf "[SphereUDE] Regularization with AD: %.9f vs %.9f using Finite Differences" l_AD l_FD + end end end @@ -155,6 +165,48 @@ function regularization(θ::ComponentVector, return reg.λ * l_ end +""" +Cubic regularization from Jupp (1987) +""" +function cubic_regularization(β::ComponentVector, + U::Chain, + st::NamedTuple, + reg::CubicSplinesRegularization, + params::AP) where {AP <: AbstractParameters} + + smodel = StatefulLuxLayer{true}(U, β.θ, st) + + # Create prediction of solution time series in integration points + nodes, weights = quadrature(params.tmin, params.tmax, params.n_quadrature) + u_ = predict(β, params, nodes) + + if typeof(reg.diff_mode) <: LuxNestedAD + # Automatic Differentiation + + if reg.diff_mode.method == "ForwardDiff" + Jac = batched_jacobian(smodel, AutoForwardDiff(), reshape(nodes, 1, params.n_quadrature)) + elseif reg.diff_mode.method == "Zygote" + Jac = batched_jacobian(smodel, AutoZygote(), reshape(nodes, 1, params.n_quadrature)) + else + throw("Method for AD backend no implemented.") + end + + L_cross_u = [cross(Jac[:,1,j], u_[:,j]) for j in 1:params.n_quadrature] + + elseif typeof(reg.diff_mode) <: FiniteDifferences + L_ = [central_fdm(τ -> smodel([τ]), t, reg.diff_mode.ϵ) for t in nodes] + L_cross_u = [cross(L_[j], u_[:,j]) for j in 1:params.n_quadrature] + + elseif typeof(reg.diff_mode) <: ComplexStepDifferentiation + L_ = [complex_step_differentiation(τ -> smodel([τ]), t, reg.diff_mode.ϵ) for t in nodes] + L_cross_u = [cross(L_[j], u_[:,j]) for j in 1:params.n_quadrature] + + else + throw("Method not implemented.") + end + + return reg.λ * sum([weights[j] * norm(L_cross_u[j])^2.0 for j in 1:params.n_quadrature]) +end """ Loss function to be called for multiple shooting diff --git a/src/train.jl b/src/train.jl index 5d4bc55..21fd947 100644 --- a/src/train.jl +++ b/src/train.jl @@ -120,7 +120,7 @@ function train(data::AD, # Final Fit fit_times = collect(range(params.tmin,params.tmax, length=1000)) - fit_directions, _ = predict(β_trained, params, fit_times) + fit_directions = predict(β_trained, params, fit_times) fit_rotations = reduce(hcat, (t -> predict_L(t, U, β_trained.θ, st)).(fit_times)) # Recover final balance between different terms involved in the loss function to assess hyperparameter selection. diff --git a/src/types.jl b/src/types.jl index 3278ee7..ea261ea 100644 --- a/src/types.jl +++ b/src/types.jl @@ -3,6 +3,7 @@ export SphereData, AbstractData export Results, AbstractResult export Regularization, AbstractRegularization export FiniteDifferences, ComplexStepDifferentiation, LuxNestedAD, AbstractDifferentiation +export CubicSplinesRegularization abstract type AbstractParameters end abstract type AbstractData end @@ -70,6 +71,11 @@ Regularization information # @assert (order == 0) || (!isnothing(diff_mode)) "Diffentiation methods needs to be provided for regularization with order larger than zero." end +@kwdef struct CubicSplinesRegularization{F <: AbstractFloat} <: AbstractRegularization + λ::F + diff_mode::Union{Nothing, AbstractDifferentiation} = nothing +end + """ Differentiation methods """ From dd1dc655ae88423b600a8e8104850d31567cc813 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Tue, 29 Oct 2024 12:47:04 -0700 Subject: [PATCH 20/22] Curl examples and reweighted loss experiment --- examples/Torsvik_2012/Gondwana.jl | 90 +++++++++++++++++++ .../{APWP-Torsvik.jl => Laurentia.jl} | 16 ++-- examples/curl/curl.jl | 60 ++++++++----- examples/double_rotation/double_rotation.jl | 16 ++-- src/losses.jl | 9 ++ src/train.jl | 14 +-- src/types.jl | 9 +- src/utils.jl | 22 +++-- 8 files changed, 180 insertions(+), 56 deletions(-) create mode 100644 examples/Torsvik_2012/Gondwana.jl rename examples/Torsvik_2012/{APWP-Torsvik.jl => Laurentia.jl} (77%) diff --git a/examples/Torsvik_2012/Gondwana.jl b/examples/Torsvik_2012/Gondwana.jl new file mode 100644 index 0000000..2eab675 --- /dev/null +++ b/examples/Torsvik_2012/Gondwana.jl @@ -0,0 +1,90 @@ +using Pkg; Pkg.activate(".") +using Revise +using Lux + +using LinearAlgebra, Statistics, Distributions +using SciMLSensitivity +# using OrdinaryDiffEqCore, OrdinaryDiffEqTsit5 +using Optimization, OptimizationOptimisers, OptimizationOptimJL + +using SphereUDE + +# Random seed +using Random +rng = Random.default_rng() +Random.seed!(rng, 613) + +using DataFrames, CSV +using Serialization, JLD2 + +df = CSV.read("./examples/Torsvik_2012/Torsvik-etal-2012_dataset.csv", DataFrame, delim=",") + +# Filter the plates that were once part of the supercontinent Gondwana + +Gondwana = ["Amazonia", "Parana", "Colorado", "Southern_Africa", + "East_Antarctica", "Madagascar", "Patagonia", "Northeast_Africa", + "Northwest_Africa", "Somalia", "Arabia", "East_Gondwana"] + +df = filter(row -> row.Plate ∈ Gondwana, df) +df.Times = df.Age .+= rand(sampler(Normal(0,0.1)), nrow(df)) # Needs to fix this! + +df = sort(df, :Times) +times = df.Times + +# Fill missing values +df.RLat .= coalesce.(df.RLat, df.Lat) +df.RLon .= coalesce.(df.RLon, df.Lon) + +X = sph2cart(Matrix(df[:,["RLat","RLon"]])'; radians=false) + +# Retrieve uncertanties from poles and convert α95 into κ +kappas = (140.0 ./ df.a95).^2 + +data = SphereData(times=times, directions=X, kappas=kappas, L=nothing) + +# Training + +# Expected maximum angular deviation in one unit of time (degrees) +Δω₀ = 1.5 +# Angular velocity +ω₀ = Δω₀ * π / 180.0 + +tspan = [times[begin], times[end]] + +# params = SphereParameters(tmin = tspan[1], tmax = tspan[2], +# reg = [Regularization(order=1, power=2.0, λ=1e5, diff_mode=FiniteDifferences(1e-4))], +# # reg = nothing, +# pretrain = false, +# u0 = [0.0, 0.0, -1.0], ωmax = ω₀, +# reltol = 1e-6, abstol = 1e-6, +# niter_ADAM = 5000, niter_LBFGS = 5000, +# sensealg = InterpolatingAdjoint(autojacvec = ReverseDiffVJP(true))) +params = SphereParameters(tmin = tspan[1], tmax = tspan[2], + reg = [Regularization(order=1, power=2.0, λ=1.0, diff_mode=FiniteDifferences(1e-4))], + # reg = nothing, + pretrain = false, + u0 = [0.0, 0.0, -1.0], ωmax = ω₀, + reltol = 1e-6, abstol = 1e-6, + niter_ADAM = 2000, niter_LBFGS = 2000, + sensealg = InterpolatingAdjoint(autojacvec = ReverseDiffVJP(true)), + hyperparameter_balance = true) + + +init_bias(rng, in_dims) = LinRange(tspan[1], tspan[2], in_dims) +init_weight(rng, out_dims, in_dims) = 0.1 * ones(out_dims, in_dims) + +# Customized neural network to similate weighted moving window in L +U = Lux.Chain( + Lux.Dense(1, 200, rbf, init_bias=init_bias, init_weight=init_weight, use_bias=true), + Lux.Dense(200,10, gelu), + Lux.Dense(10, 3, Base.Fix2(sigmoid_cap, params.ωmax), use_bias=false) +) + +results = train(data, params, rng, nothing, U) +results_dict = convert2dict(data, results) + + +# JLD2.@save "examples/Torsvik_2012/results/results_dict.jld2" results_dict + +plot_sphere(data, results, -30., 0., saveas="examples/Torsvik_2012/plots/plot_sphere.pdf", title="Double rotation") +plot_L(data, results, saveas="examples/Torsvik_2012/plots/plot_L.pdf", title="Double rotation") diff --git a/examples/Torsvik_2012/APWP-Torsvik.jl b/examples/Torsvik_2012/Laurentia.jl similarity index 77% rename from examples/Torsvik_2012/APWP-Torsvik.jl rename to examples/Torsvik_2012/Laurentia.jl index dbb30fb..da77b77 100644 --- a/examples/Torsvik_2012/APWP-Torsvik.jl +++ b/examples/Torsvik_2012/Laurentia.jl @@ -21,11 +21,9 @@ df = CSV.read("./examples/Torsvik_2012/Torsvik-etal-2012_dataset.csv", DataFrame # Filter the plates that were once part of the supercontinent Gondwana -Gondwana = ["Amazonia", "Parana", "Colorado", "Southern_Africa", - "East_Antarctica", "Madagascar", "Patagonia", "Northeast_Africa", - "Northwest_Africa", "Somalia", "Arabia", "East_Gondwana"] +Laurentia = ["north_america", "greenland"] -df = filter(row -> row.Plate ∈ Gondwana, df) +df = filter(row -> row.Plate ∈ Laurentia, df) df.Times = df.Age .+= rand(sampler(Normal(0,0.1)), nrow(df)) # Needs to fix this! df = sort(df, :Times) @@ -52,7 +50,7 @@ data = SphereData(times=times, directions=X, kappas=kappas, L=nothing) tspan = [times[begin], times[end]] params = SphereParameters(tmin = tspan[1], tmax = tspan[2], - reg = [Regularization(order=1, power=2.0, λ=1e5, diff_mode=FiniteDifferences(1e-4))], + reg = [Regularization(order=1, power=2.0, λ=1e6, diff_mode=FiniteDifferences(1e-4))], # reg = nothing, pretrain = false, u0 = [0.0, 0.0, -1.0], ωmax = ω₀, @@ -74,10 +72,8 @@ U = Lux.Chain( results = train(data, params, rng, nothing, U) results_dict = convert2dict(data, results) -# JLD2.@save "examples/Torsvik_2012/results/data.jld2" data -# JLD2.@save "examples/Torsvik_2012/results/results.jld2" results -JLD2.@save "examples/Torsvik_2012/results/results_dict.jld2" results_dict +# JLD2.@save "examples/Torsvik_2012/results/results_dict.jld2" results_dict -plot_sphere(data, results, -30., 0., saveas="examples/Torsvik_2012/plots/plot_sphere.pdf", title="Double rotation") -plot_L(data, results, saveas="examples/Torsvik_2012/plots/plot_L.pdf", title="Double rotation") +plot_sphere(data, results, -30., 0., saveas="examples/Torsvik_2012/plots/Laurentia_plot_sphere.pdf", title="Double rotation") +plot_L(data, results, saveas="examples/Torsvik_2012/plots/Laurentia_plot_L.pdf", title="Double rotation") diff --git a/examples/curl/curl.jl b/examples/curl/curl.jl index 76d81f0..cab8fd9 100644 --- a/examples/curl/curl.jl +++ b/examples/curl/curl.jl @@ -2,9 +2,10 @@ using Pkg; Pkg.activate(".") using Revise using LinearAlgebra, Statistics, Distributions -using OrdinaryDiffEq using SciMLSensitivity +using OrdinaryDiffEqCore, OrdinaryDiffEqTsit5 using Optimization, OptimizationOptimisers, OptimizationOptimJL +using Lux using SphereUDE @@ -15,27 +16,29 @@ using SphereUDE # Random seed using Random rng = Random.default_rng() -Random.seed!(rng, 616) +Random.seed!(rng, 333) # Total time simulation tspan = [0, 100.0] # Number of sample points N_samples = 300 # Times where we sample points -times_samples = sort(rand(sampler(Uniform(tspan[1], tspan[2])), N_samples)) +times_samples = collect(LinRange(tspan[1], tspan[2], N_samples)) # Expected maximum angular deviation in one unit of time (degrees) -Δω₀ = 1.0 -# Angular velocity s +Δω₀ = 10.0 +# Angular velocity ω₀ = Δω₀ * π / 180.0 +# Angular momentum # Solver tolerances -reltol = 1e-7 -abstol = 1e-7 +reltol = 1e-6 +abstol = 1e-6 function L_true(t::Float64) - τ = π * (t / tspan[2]) - ω₀ * [sin(τ), cos(τ), 0] + lon = 0.0 + lat = -40.0 * (t - tspan[2]) / (tspan[2] - tspan[1]) - 40.0 * (t - tspan[1]) / (tspan[2] - tspan[1]) + return ω₀ * sph2cart([lat, lon], radians=false) end function true_rotation!(du, u, p, t) @@ -43,32 +46,45 @@ function true_rotation!(du, u, p, t) du .= cross(L, u) end -prob = ODEProblem(true_rotation!, [0.0, 0.0, 1.0], tspan) -true_sol = solve(prob, Tsit5(), reltol=reltol, abstol=abstol, saveat=times_samples) +x0 = [0.5, 0.0, 0.7] +x0 /= norm(x0) + +prob = ODEProblem(true_rotation!, x0, tspan) +true_sol = solve(prob, Tsit5(), reltol = reltol, abstol = abstol, saveat = times_samples) # Add Fisher noise to true solution X_noiseless = Array(true_sol) -X_true = X_noiseless #+ FisherNoise(kappa=10000.) +X_true = X_noiseless #+ FisherNoise(kappa=1000.0) ############################################################## ####################### Training ########################### ############################################################## -data = SphereData(times=times_samples, directions=X_true, kappas=nothing, L=L_true) +data = SphereData(times=times_samples, directions=X_true, kappas=nothing, L=L_true) + +params = SphereParameters(tmin = tspan[1], tmax = tspan[2], + reg = [Regularization(order=1, power=2.0, λ=1e1, diff_mode=LuxNestedAD())], + u0 = [0.0, 0.0, -1.0], + train_initial_condition = true, + ωmax = ω₀, reltol = reltol, abstol = abstol, + niter_ADAM = 2000, niter_LBFGS = 1000, + pretrain = true, + sensealg = QuadratureAdjoint(autojacvec = ReverseDiffVJP(true))) -regs = [Regularization(order=1, power=1.0, λ=0.001, diff_mode="Finite Differences"), - Regularization(order=0, power=2.0, λ=0.1, diff_mode="Finite Differences")] +init_bias(rng, in_dims) = LinRange(tspan[1], tspan[2], in_dims) +init_weight(rng, out_dims, in_dims) = 0.1 * ones(out_dims, in_dims) -params = SphereParameters(tmin=tspan[1], tmax=tspan[2], - reg=regs, - u0=[0.0, 0.0, 1.0], ωmax=10*ω₀, reltol=reltol, abstol=abstol, - niter_ADAM=100, niter_LBFGS=60) +U = Lux.Chain( + Lux.Dense(1, 200, rbf, init_bias=init_bias, init_weight=init_weight, use_bias=true), + Lux.Dense(200,10, relu), + Lux.Dense(10, 3, Base.Fix2(sigmoid_cap, params.ωmax), use_bias=false) +) -results = train(data, params, rng, nothing) +results = train(data, params, rng, nothing, U) ############################################################## ###################### PyCall Plots ######################### ############################################################## -plot_sphere(data, results, 0., 0., saveas="examples/curl/plot_sphere.pdf", title="Double rotation") -plot_L(data, results, saveas="examples/curl/plot_L.pdf", title="Double rotation") \ No newline at end of file +plot_sphere(data, results, 0.0, 0.0, saveas="examples/curl/_sphere.pdf", title="Curl") # , matplotlib_rcParams=Dict("font.size"=> 50)) +plot_L(data, results, saveas="examples/curl/_L.pdf", title="Curl") diff --git a/examples/double_rotation/double_rotation.jl b/examples/double_rotation/double_rotation.jl index 0e4c5cb..77c256a 100644 --- a/examples/double_rotation/double_rotation.jl +++ b/examples/double_rotation/double_rotation.jl @@ -72,15 +72,15 @@ params = SphereParameters(tmin = tspan[1], tmax = tspan[2], train_initial_condition = false, multiple_shooting = false, u0 = [0.0, 0.0, -1.0], ωmax = ω₀, reltol = reltol, abstol = abstol, - niter_ADAM = 5000, niter_LBFGS = 5000, - sensealg = GaussAdjoint(autojacvec = ReverseDiffVJP(true))) + niter_ADAM = 2000, niter_LBFGS = 1000, + sensealg = QuadratureAdjoint(autojacvec = ReverseDiffVJP(true))) init_bias(rng, in_dims) = LinRange(tspan[1], tspan[2], in_dims) init_weight(rng, out_dims, in_dims) = 0.1 * ones(out_dims, in_dims) U = Lux.Chain( Lux.Dense(1, 200, rbf, init_bias=init_bias, init_weight=init_weight, use_bias=true), - Lux.Dense(200,10, gelu), + Lux.Dense(200,10, relu), Lux.Dense(10, 3, Base.Fix2(sigmoid_cap, params.ωmax), use_bias=false) ) @@ -140,18 +140,18 @@ end # run ### AD run(; kappa = 50., - regs = [Regularization(order=1, power=1.0, λ=0.01, diff_mode=LuxNestedAD())], %, + regs = [Regularization(order=1, power=1.0, λ=0.1, diff_mode=LuxNestedAD())], # Regularization(order=0, power=2.0, λ=0.1)], title = "plots/AD_plot_50") run(; kappa = 200., - regs = [Regularization(order=1, power=1.0, λ=0.1, diff_mode=LuxNestedAD()), - Regularization(order=0, power=2.0, λ=0.1)], + regs = [Regularization(order=1, power=1.0, λ=0.1, diff_mode=LuxNestedAD())], + # Regularization(order=0, power=2.0, λ=0.1)], title = "plots/AD_plot_200") run(; kappa = 1000., - regs = [Regularization(order=1, power=1.0, λ=0.1, diff_mode=LuxNestedAD()), - Regularization(order=0, power=2.0, λ=0.1)], + regs = [Regularization(order=1, power=1.0, λ=0.1, diff_mode=LuxNestedAD())], + # Regularization(order=0, power=2.0, λ=0.1)], title = "plots/AD_plot_1000") diff --git a/src/losses.jl b/src/losses.jl index 127164f..9e9f4c3 100644 --- a/src/losses.jl +++ b/src/losses.jl @@ -19,6 +19,11 @@ function loss(β::ComponentVector, l_emp = loss_empirical(β, data, params) loss_dict["Empirical"] = l_emp + + if params.hyperparameter_balance + # @ignore_derivatives l_emp /= l_emp + l_emp = log(l_emp) + end # Regularization l_reg = 0.0 @@ -34,6 +39,10 @@ function loss(β::ComponentVector, throw("Regularization not implemented.") end # Add contribution to regularization + if params.hyperparameter_balance + reg₀ = log(reg₀) + # @ignore_derivatives reg₀ /= reg₀ + end l_reg += reg₀ end end diff --git a/src/train.jl b/src/train.jl index 21fd947..04ed441 100644 --- a/src/train.jl +++ b/src/train.jl @@ -49,8 +49,9 @@ function train(data::AD, losses = Float64[] callback = function (p, l) push!(losses, l) - if length(losses) % 50 == 0 - @printf "Iteration: [%5d / %5d] \t Loss: %.9f \n" length(losses) (params.niter_ADAM+params.niter_LBFGS) losses[end] + if length(losses) % 100 == 0 + _, l_dict = f_loss(p) + @printf "Iteration: [%5d / %5d] \t Loss: %.9f = Empirical: %.9f + Regularization: %.9f \n" length(losses) (params.niter_ADAM+params.niter_LBFGS) sum(values(l_dict)) l_dict["Empirical"] (sum(values(l_dict))-l_dict["Empirical"]) end if params.train_initial_condition p.u0 ./= norm(p.u0) @@ -66,8 +67,6 @@ function train(data::AD, f_loss(β) = loss(β, data, params, U, st) end - println("Loss after initalization: ", f_loss(β)[1]) - """ Pretraining to find parameters without impossing regularization """ @@ -81,7 +80,9 @@ function train(data::AD, end return false end - optf₀ = Optimization.OptimizationFunction((x, β) -> loss_empirical(x), params.adtype) + # Define the loss function with just empirical component + f_loss_empirical(β) = loss_empirical(β, data, params) + optf₀ = Optimization.OptimizationFunction((x, β) -> f_loss_empirical(x), params.adtype) optprob₀ = Optimization.OptimizationProblem(optf₀, β) res₀ = Optimization.solve(optprob₀, ADAM(), callback=callback_pretrain, maxiters=params.niter_ADAM, verbose=false) optprob₁ = Optimization.OptimizationProblem(optf₀, res₀.u) @@ -89,8 +90,7 @@ function train(data::AD, β = res₁.u end - # To do: implement this with polyoptimizaion to put ADAM and BFGS in one step. - # Maybe better to keep like this for the line search. + println("Loss after initalization: ", f_loss(β)[1]) optf = Optimization.OptimizationFunction((x, β) -> (first ∘ f_loss)(x), params.adtype) optprob = Optimization.OptimizationProblem(optf, β) diff --git a/src/types.jl b/src/types.jl index ea261ea..418f455 100644 --- a/src/types.jl +++ b/src/types.jl @@ -22,15 +22,16 @@ Training parameters reg::Union{Nothing, Array} train_initial_condition::Bool = false multiple_shooting::Bool = false - niter_ADAM::I - niter_LBFGS::I - reltol::F - abstol::F + niter_ADAM::I = 2000 + niter_LBFGS::I = 2000 + reltol::F = 1e-6 + abstol::F = 1e-6 n_quadrature::I = 100 solver::OrdinaryDiffEqCore.OrdinaryDiffEqAlgorithm = Tsit5() adtype::Optimization.AbstractADType = AutoZygote() sensealg::SciMLBase.AbstractAdjointSensitivityAlgorithm = QuadratureAdjoint(autojacvec=ReverseDiffVJP(true)) pretrain::Bool = false + hyperparameter_balance::Bool = false end """ diff --git a/src/utils.jl b/src/utils.jl index 7258a36..1ae2164 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -8,6 +8,7 @@ export raise_warnings export isL1reg export convert2dict export get_default_NN +export fisher_mean # Import activation function for complex extension import Lux: relu, gelu @@ -177,6 +178,16 @@ function sph2cart(X::AbstractArray{<:Number}; radians::Bool=true) return Y end +""" +Return Fisher mean on the sphere +""" +function fisher_mean(latitudes, longitudes; radians::Bool=true) + Y = sph2cart(hcat(latitudes, longitudes)', radians=radians) + Ŷ = mean(Y, dims=2) + Ŷ ./= norm(Ŷ) + return cart2sph(Ŷ, radians=radians) +end + """ Add Fisher noise to matrix of three dimensional unit vectors @@ -218,11 +229,12 @@ function quadrature(f::Function, t₀, t₁, n_quadrature::Int) end function quadrature(t₀, t₁, n_quadrature::Int) - ignore() do - # Ignore AD here since FastGaussQuadrature is using mutating arrays - nodes, weigths = gausslegendre(n_quadrature) - # nodes .+ rand(sampler(Uniform(-0.1,0.1)), n_quadrature) - end + # Ignore AD here since FastGaussQuadrature is using mutating arrays + @ignore_derivatives nodes, weigths = gausslegendre(n_quadrature) + # ignore() do + # nodes, weigths = gausslegendre(n_quadrature) + # # nodes .+ rand(sampler(Uniform(-0.1,0.1)), n_quadrature) + # end nodes = (t₀+t₁)/2 .+ nodes * (t₁-t₀)/2 weigths = (t₁-t₀) / 2 * weigths return nodes, weigths From ee48555e56fd62fe2f643c7978c869e222a73bb1 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Wed, 6 Nov 2024 13:21:51 -0800 Subject: [PATCH 21/22] minimal changes --- src/plot.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/plot.jl b/src/plot.jl index ddfa06d..2e7ff30 100644 --- a/src/plot.jl +++ b/src/plot.jl @@ -61,7 +61,7 @@ function plot_sphere(# ax::PyCall.PyObject, linewidth=2, color="black",#cmap(norm(results.fit_times[i])), transform = ccrs.Geodetic()) end - plt.title(title, fontsize=20) + # plt.title(title, fontsize=20) if !isnothing(saveas) plt.savefig(saveas, format="pdf") end @@ -97,7 +97,7 @@ function plot_L(data::AbstractData, plt.xlabel("Time") plt.ylabel("Angular Velocity") plt.legend() - plt.title(title) + # plt.title(title) if !isnothing(saveas) plt.savefig(saveas, format="pdf") From 4610185f218061a2a50d41643e928b5ced01f1f5 Mon Sep 17 00:00:00 2001 From: Facundo Sapienza Date: Wed, 6 Nov 2024 13:22:33 -0800 Subject: [PATCH 22/22] Bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 98d68a1..d2ed2f5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SphereUDE" uuid = "d7416ba7-148a-4110-b27d-9087fcebab2d" authors = ["Facundo Sapienza ", "Jordi Bolibar "] -version = "0.1.2" +version = "0.1.3" [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"