Add hessian_sparsity method for NonlinearSystem (#1602)

Author: Vaibhavdixit02

Project.toml

@@ -58,7 +58,7 @@ DocStringExtensions = "0.7, 0.8, 0.9"
 DomainSets = "0.5"
 Graphs = "1.5.2"
 IfElse = "0.1"
-JuliaFormatter = "0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 1"
+JuliaFormatter = "1"
 LabelledArrays = "1.3"
 Latexify = "0.11, 0.12, 0.13, 0.14, 0.15"
 MacroTools = "0.5"

src/systems/nonlinear/nonlinearsystem.jl

@@ -172,6 +172,11 @@ function jacobian_sparsity(sys::NonlinearSystem)
                       states(sys))
 end
 
+function hessian_sparsity(sys::NonlinearSystem)
+    [hessian_sparsity(eq.rhs,
+                      states(sys)) for eq in equations(sys)]
+end
+
 function DiffEqBase.NonlinearFunction(sys::NonlinearSystem, args...; kwargs...)
     NonlinearFunction{true}(sys, args...; kwargs...)
 end

src/systems/optimization/optimizationsystem.jl

@@ -184,15 +184,10 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
         _hess = nothing
     end
 
-    _f = DiffEqBase.OptimizationFunction{iip, SciMLBase.NoAD, typeof(f), typeof(_grad),
-                                         typeof(_hess), Nothing, Nothing, Nothing, Nothing}(f,
-                                                                                            SciMLBase.NoAD(),
-                                                                                            _grad,
-                                                                                            _hess,
-                                                                                            nothing,
-                                                                                            nothing,
-                                                                                            nothing,
-                                                                                            nothing)
+    _f = DiffEqBase.OptimizationFunction{iip}(f,
+                                              SciMLBase.NoAD();
+                                              grad = _grad,
+                                              hess = _hess)
 
     defs = defaults(sys)
     defs = mergedefaults(defs, parammap, ps)
@@ -272,13 +267,9 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0,
         hess = $_hess
         lb = $lb
         ub = $ub
-        _f = OptimizationFunction{$iip, typeof(f), typeof(grad), typeof(hess),
-                                  SciMLBase.NoAD, Nothing, Nothing, Nothing}(f, grad, hess,
-                                                                             nothing,
-                                                                             SciMLBase.NoAD(),
-                                                                             nothing,
-                                                                             nothing,
-                                                                             nothing, 0)
+        _f = OptimizationFunction{iip}(f, SciMLBase.NoAD();
+                                       grad = grad,
+                                       hess = hess)
         OptimizationProblem{$iip}(_f, u0, p; lb = lb, ub = ub, kwargs...)
     end
 end

test/nonlinearsystem.jl

@@ -72,6 +72,12 @@ jac = calculate_jacobian(ns)
 
 jac = generate_jacobian(ns)
 
+sH = calculate_hessian(ns)
+@test getfield.(ModelingToolkit.hessian_sparsity(ns), :colptr) ==
+      getfield.(sparse.(sH), :colptr)
+@test getfield.(ModelingToolkit.hessian_sparsity(ns), :rowval) ==
+      getfield.(sparse.(sH), :rowval)
+
 prob = NonlinearProblem(ns, ones(3), ones(3))
 sol = solve(prob, NewtonRaphson())
 @test sol.u[1] ≈ sol.u[2]