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" 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/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/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/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..9e9f4c3 --- /dev/null +++ b/src/losses.jl @@ -0,0 +1,275 @@ +export loss +export loss_empirical +export regularization +export cubic_regularization +export loss_multiple_shooting + +""" +General Loss Function +""" +function loss(β::ComponentVector, + data::AD, + params::AP, + 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 + + if params.hyperparameter_balance + # @ignore_derivatives l_emp /= l_emp + l_emp = log(l_emp) + end + + # Regularization + l_reg = 0.0 + if !isnothing(params.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 + if params.hyperparameter_balance + reg₀ = log(reg₀) + # @ignore_derivatives reg₀ /= reg₀ + end + l_reg += 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_ = predict(β, params, data.times) + + # 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, + 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 + + # 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) + + # 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 + +""" +Regularization +""" +function regularization(θ::ComponentVector, + U::Chain, + st::NamedTuple, + reg::Regularization, + params::AP) where {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 + 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 + + 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 + +""" +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 +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/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") diff --git a/src/train.jl b/src/train.jl index 23f120e..04ed441 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,29 @@ 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 - - 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 - + global prob_nn = ODEProblem(ude_rotation!, params.u0, [params.tmin, params.tmax], β.θ) - ### 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])") + 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) @@ -270,14 +59,14 @@ 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 + """ Pretraining to find parameters without impossing regularization """ @@ -291,7 +80,9 @@ function train(data::AD, end return false end - optf₀ = Optimization.OptimizationFunction((x, β) -> loss_empirical(x), 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) @@ -299,10 +90,9 @@ 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), 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..418f455 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 @@ -21,13 +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 """ @@ -68,6 +72,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 """ diff --git a/src/utils.jl b/src/utils.jl index 47b86bc..1ae2164 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -7,9 +7,12 @@ export quadrature, central_fdm, complex_step_differentiation export raise_warnings export isL1reg export convert2dict +export get_default_NN +export fisher_mean # Import activation function for complex extension import Lux: relu, gelu +using Lux: Chain # import Lux: sigmoid, relu, gelu ### Custom Activation Funtions @@ -109,6 +112,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 @@ -143,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 @@ -178,16 +223,18 @@ 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) - ignore() do - # Ignore AD here since FastGaussQuadrature is using mutating arrays - nodes, weigths = gausslegendre(n_nodes) - end +function quadrature(t₀, t₁, n_quadrature::Int) + # 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 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