Merge pull request #24 from DJ4Earth/testoptimization

Testoptimization

Author: enigne

Project.toml

@@ -9,11 +9,12 @@ ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
 Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+JSOSolvers = "10dff2fc-5484-5881-a0e0-c90441020f8a"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
 MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+ManualNLPModels = "30dfa513-9b2f-4fb3-9796-781eabac1617"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
-OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
@@ -23,9 +24,11 @@ Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
 ColorSchemes = "3.25"
 Enzyme = "0.13"
 Flux = "0.16"
+JSOSolvers = "0.14.1"
 MAT = "0.10"
 MadNLP = "0.8"
 Makie = "0.22"
+ManualNLPModels = "0.2.0"
 SparseArrays = "1"
 StatsBase = "0.34"
 Triangulate = "2.3"

src/usr/classes.jl

@@ -357,10 +357,10 @@ function model(md::model; mesh::AbstractMesh=md.mesh, friction::AbstractFriction
 					 md.masstransport, md.transient, md.inversion, md.calving, 
 					 md.levelset, md.frontalforcings)
 end#}}}
-function model(matmd::Dict,verbose::Bool=true) #{{{
+function model(matmd::Dict; verbose::Bool=true, friction::AbstractFriction=BuddFriction()) #{{{
 
 	#initialize output
-	md = model()
+	md = model(model(), friction=friction)
 
 	#Loop over all possible fields
 	for name1 in keys(matmd)

test/testoptimization.jl

@@ -3,7 +3,10 @@ using DJUICE
 using MAT
 using Test
 
-using Optimization, OptimizationOptimJL
+#using Optimization, OptimizationNLopt
+using ManualNLPModels, JSOSolvers
+using MadNLP
+
 
 #Load model from MATLAB file
 file = matopen(joinpath(@__DIR__, "..", "data","temp.mat")) #SMALL model (35 elements)
@@ -26,13 +29,36 @@ md.inversion.dependent_string = ["SurfaceAbsVelMisfit"]
 α = md.inversion.independent
 ∂J_∂α = zero(α)
 
-femmodel=DJUICE.ModelProcessor(md, :StressbalanceSolution)
-n = length(α)
 
-DJUICE.costfunction(α, femmodel)
-# test Enzyme autodiff only
-dfemmodel = Enzyme.Compiler.make_zero(Base.Core.Typeof(femmodel), IdDict(), femmodel)
-autodiff(set_runtime_activity(Enzyme.Reverse), DJUICE.costfunction, Active, Duplicated(α, ∂J_∂α), Duplicated(femmodel,dfemmodel))
+# cost function f
+f(x) = begin
+	femmodel=DJUICE.ModelProcessor(md, :StressbalanceSolution)
+	DJUICE.costfunction(x, femmodel)
+end
+
+# Enzyme gradient g
+g!(gx, x) = begin
+	femmodel=DJUICE.ModelProcessor(md, :StressbalanceSolution)
+	dfemmodel = Enzyme.Compiler.make_zero(Base.Core.Typeof(femmodel), IdDict(), femmodel)
+	Enzyme.autodiff(set_runtime_activity(Enzyme.Reverse), DJUICE.costfunction, Active, Duplicated(x, gx), Duplicated(femmodel,dfemmodel))
+	gx
+end
+
+nlp = NLPModel(
+  α,
+  f,
+  grad = g!,
+)
+
+#output = lbfgs(nlp)
+
+results_qn = madnlp(
+           nlp;
+           linear_solver=LapackCPUSolver,
+           hessian_approximation=MadNLP.CompactLBFGS,
+       )
+
+
 
 # use user defined grad, errors!
 #optprob = OptimizationFunction(DJUICE.costfunction, Optimization.AutoEnzyme(; mode=set_runtime_activity(Enzyme.Reverse)))
@@ -44,3 +70,9 @@ autodiff(set_runtime_activity(Enzyme.Reverse), DJUICE.costfunction, Active, Dupl
 #sol = Optimization.solve(prob, Optim.NelderMead())
 #sol = Optimization.solve(prob, Optimization.LBFGS(), maxiters = 100)
 
+#using JuMP, Optim
+#m = Model(Optim.Optimizer);
+#set_optimizer_attribute(m, "method", LBFGS())
+#@variable(m, α)
+#@variable(m, femmodel)
+#@objective(m, Min, DJUICE.costfunction(α, femmodel))