From 1e4f2ec1a579d275400c59bda61977e48a53db4c Mon Sep 17 00:00:00 2001
From: Patrick Kofod Mogensen <patrick.mogensen@gmail.com>
Date: Sat, 20 Jan 2024 16:24:37 +0100
Subject: [PATCH 1/3] Adjust Trust region solver to allow for nonallocating
 dogleg

---
 Project.toml                                  |  4 ++
 src/NLSolvers.jl                              | 16 ++++++--
 .../trs_solvers/solvers/Dogleg.jl             |  9 ++--
 src/globalization/trs_solvers/solvers/NTR.jl  |  1 +
 src/globalization/trs_solvers/solvers/NWI.jl  |  2 +-
 src/globalization/trs_solvers/solvers/TCG.jl  |  1 +
 src/nlsolve/trustregions/newton.jl            | 31 +++++---------
 src/objectives.jl                             | 15 -------
 .../trustregions/optimize/inplace_loop.jl     | 16 ++++----
 src/quasinewton/approximations/newton.jl      |  2 +-
 .../truncatedconjugategradient.jl             | 10 ++---
 test/nlsolve/interface.jl                     | 41 ++++++++++++++++---
 test/nlsolve/problems.jl                      | 28 +++++++++++++
 test/runtests.jl                              |  5 ++-
 14 files changed, 115 insertions(+), 66 deletions(-)

diff --git a/Project.toml b/Project.toml
index 3a038f1..04d17bb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,10 +4,14 @@ authors = ["Patrick Kofod Mogensen <patrick.mogensen@gmail.com>"]
 version = "0.4.0"
 
 [deps]
+Cthulhu = "f68482b8-f384-11e8-15f7-abe071a5a75f"
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 PositiveFactorizations = "85a6dd25-e78a-55b7-8502-1745935b8125"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
diff --git a/src/NLSolvers.jl b/src/NLSolvers.jl
index 2aad83a..dda7ff7 100644
--- a/src/NLSolvers.jl
+++ b/src/NLSolvers.jl
@@ -185,9 +185,6 @@ include("lsqsolve/problem_types.jl")
 include("lsqsolve/optimization.jl")
 export LeastSquaresProblem, LeastSquaresOptions
 
-function negate(problem::AbstractProblem, A)
-
-end
 function negate(::InPlace, A::AbstractArray)
     A .= .-A
     return A
@@ -195,7 +192,18 @@ end
 function negate(::OutOfPlace, A::AbstractArray)
     -A
 end
-
+function scale(::InPlace, B, A, m)
+    B .= A .* m
+end
+function scale(::OutOfPlace, B, A, m)
+    A .* m
+end
+function _copyto(::InPlace, z, x)
+    copyto!(z, x)
+end
+function _copyto(::OutOfPlace, z, x)
+    copy(x)
+end
 mstyle(problem::AbstractProblem) = problem.mstyle
 
 """
diff --git a/src/globalization/trs_solvers/solvers/Dogleg.jl b/src/globalization/trs_solvers/solvers/Dogleg.jl
index 4554013..134d001 100644
--- a/src/globalization/trs_solvers/solvers/Dogleg.jl
+++ b/src/globalization/trs_solvers/solvers/Dogleg.jl
@@ -20,7 +20,7 @@ struct Dogleg{T} <: TRSPSolver
 end
 Dogleg() = Dogleg(nothing)
 
-function (dogleg::Dogleg)(∇f, H, Δ, p, scheme; abstol = 1e-10, maxiter = 50)
+function (dogleg::Dogleg)(∇f, H, Δ, p, scheme, mstyle; abstol = 1e-10, maxiter = 50)
     T = eltype(p)
     n = length(∇f)
 
@@ -32,7 +32,8 @@ function (dogleg::Dogleg)(∇f, H, Δ, p, scheme; abstol = 1e-10, maxiter = 50)
     norm_d_cauchy = norm(d_cauchy)
     if norm_d_cauchy ≥ Δ
         shrink = Δ / norm_d_cauchy # inv(Δ/norm_d_cauchy) puts it on the border
-        p .= d_cauchy .* shrink
+
+        p = scale(mstyle, p, d_cauchy, shrink)
         interior = false
     else
         # Else, calculate (Quasi-)Newton step. If this is interior, then take the
@@ -42,7 +43,6 @@ function (dogleg::Dogleg)(∇f, H, Δ, p, scheme; abstol = 1e-10, maxiter = 50)
         p = find_direction!(p, H, nothing, ∇f, scheme)
         norm_p = norm(p)
         if norm_p ≤ Δ # fixme really need to add the 20% slack here (see TR book and NTR)
-            p .= p
             if norm_p < Δ
                 interior = true
             else
@@ -75,7 +75,8 @@ function (dogleg::Dogleg)(∇f, H, Δ, p, scheme; abstol = 1e-10, maxiter = 50)
             else
                 t = T(0)
             end
-            p .= d_cauchy .+ t .* p
+
+            p, _ = move(mstyle, p, d_cauchy, p, p, t)
             interior = false
         end
     end
diff --git a/src/globalization/trs_solvers/solvers/NTR.jl b/src/globalization/trs_solvers/solvers/NTR.jl
index 367f4c3..9dc42cc 100644
--- a/src/globalization/trs_solvers/solvers/NTR.jl
+++ b/src/globalization/trs_solvers/solvers/NTR.jl
@@ -51,6 +51,7 @@ function (ms::NTR)(
     Δ::T,
     s,
     scheme,
+    mstyle,
     λ0 = 0;
     abstol = 1e-10,
     maxiter = 50,
diff --git a/src/globalization/trs_solvers/solvers/NWI.jl b/src/globalization/trs_solvers/solvers/NWI.jl
index 94f86a2..c7ec478 100644
--- a/src/globalization/trs_solvers/solvers/NWI.jl
+++ b/src/globalization/trs_solvers/solvers/NWI.jl
@@ -148,7 +148,7 @@ Returns:
     solved - Whether or not a solution was reached (as opposed to
       terminating early due to maxiter)
 """
-function (ms::NWI)(∇f, H, Δ, p, scheme; abstol = 1e-10, maxiter = 50)
+function (ms::NWI)(∇f, H, Δ, p, scheme, mstyle; abstol = 1e-10, maxiter = 50)
     T = eltype(p)
     n = length(∇f)
     H = H isa UniformScaling ? Diagonal(copy(∇f) .* 0 .+ 1) : H
diff --git a/src/globalization/trs_solvers/solvers/TCG.jl b/src/globalization/trs_solvers/solvers/TCG.jl
index a65d040..94d1405 100644
--- a/src/globalization/trs_solvers/solvers/TCG.jl
+++ b/src/globalization/trs_solvers/solvers/TCG.jl
@@ -22,6 +22,7 @@ function (ms::TCG)(
     Δ::T,
     s,
     scheme,
+    mstyle,
     λ0 = 0;
     abstol = 1e-10,
     maxiter = min(5, length(s)),
diff --git a/src/nlsolve/trustregions/newton.jl b/src/nlsolve/trustregions/newton.jl
index bc983da..af25da8 100644
--- a/src/nlsolve/trustregions/newton.jl
+++ b/src/nlsolve/trustregions/newton.jl
@@ -12,10 +12,6 @@ has_param(::NormedResiduals) = false
 function value(nr::NormedResiduals, x)
     Fx = nr.F.F(nr.Fx, x)
 
-    # this is just to grab them outside, but this hsould come from the convergence info perhaps?
-    nr.x .= x
-    nr.Fx .= Fx
-
     f = (norm(Fx)^2) / 2
     return f
 end
@@ -24,23 +20,17 @@ function upto_gradient(nr::NormedResiduals, Fx, x)
     Fx, Jx = nr.F.FJ(nr.Fx, nr.Fx * x', x)
 
     # this is just to grab them outside, but this hsould come from the convergence info perhaps?
-    nr.x .= x
-    nr.Fx .= Fx
-
     f = (norm(Fx)^2) / 2
-    Fx .= Jx' * Fx
-    return f, Fx
+    JxtFx = Jx' * Fx
+    return f, JxtFx
 end
 function upto_hessian(nr::NormedResiduals, Fx, Jx, x)  #Fx is the gradient and Jx is the Hessian
     Fx, Jx = nr.F.FJ(nr.Fx, Jx, x)
 
     # this is just to grab them outside, but this hsould come from the convergence info perhaps?
-    nr.x .= x
-    nr.Fx .= Fx
-
     f = (norm(Fx)^2) / 2
     # this is the gradient
-    Fx .= Jx' * Fx
+    JxtFx = Jx' * Fx
     # As you may notice, this can be expensive... Because the gradient
     # is going to be very simple. May want to create a
     # special type or way to hook into trust regions here. We can exploit
@@ -48,8 +38,8 @@ function upto_hessian(nr::NormedResiduals, Fx, Jx, x)  #Fx is the gradient and J
     # systems. There is no need to get the full hessian. because these two
     # steps are don't need these multiplies
     # This is the Hessian
-    Jx .= Jx' * Jx
-    return f, Fx, Jx
+    Jx2 = Jx' * Jx
+    return f, JxtFx, Jx2
 end
 
 function solve(
@@ -58,7 +48,7 @@ function solve(
     approach::TrustRegion{<:Union{SR1,DBFGS,BFGS,Newton},<:Any,<:Any},
     options::NEqOptions,
 )
-    if !(mstyle(prob) === InPlace())
+  if !(mstyle(prob) === InPlace()) && !(approach.spsolve isa Dogleg)
         throw(
             ErrorException(
                 "solve() not defined for OutOfPlace() with Trustregion for NEqProblem",
@@ -75,17 +65,16 @@ function solve(
     normed_residual = NormedResiduals(x_outer, Fx_outer, F)
     ρ2F0 = sqrt(value(normed_residual, x_outer) * 2)
     ρF0 = norm(normed_residual.Fx, Inf)
-    td = OptimizationProblem(normed_residual)
+    td = OptimizationProblem(normed_residual, inplace = mstyle(prob) == InPlace())
+    options.maxiter
     res = solve(td, x, approach, OptimizationOptions(maxiter = options.maxiter))
-    # normed_residual to get Fx
-    # f0*2 is wrong because it's 2norm
     newinfo = (
         solution = solution(res),
-        best_residual = Fx_outer,
+        best_residual = value(F, Fx_outer, solution(res)),
         ρF0 = ρF0,
         ρ2F0 = ρ2F0,
         time = res.info.time,
         iter = res.info.iter,
     )
-    return ConvergenceInfo(approach, newinfo, options)
+    return  ConvergenceInfo(approach, newinfo, options)
 end
diff --git a/src/objectives.jl b/src/objectives.jl
index 4b4ed8e..41f2037 100644
--- a/src/objectives.jl
+++ b/src/objectives.jl
@@ -177,21 +177,6 @@ function _value(mo, R::VectorObjective, Fx, x)
     (ϕ = (norm(Fx)^2) / 2, Fx = Fx)
 end
 
-struct LsqWrapper{Tobj,TF,TJ} <: ObjWrapper
-    R::Tobj
-    F::TF
-    J::TJ
-end
-function (lw::LsqWrapper)(x)
-    F = lw.R(lw.F, x)
-    sum(abs2, F) / 2
-end
-function (lw::LsqWrapper)(∇f, x)
-    _F, _J = lw.R(lw.F, lw.J, x)
-    copyto!(∇f, sum(_J; dims = 1))
-    sum(abs2, _F), ∇f
-end
-
 struct LeastSquaresObjective{TFx,TJx,Tf,Tfj,Td}
     Fx::TFx
     Jx::TJx
diff --git a/src/optimize/trustregions/optimize/inplace_loop.jl b/src/optimize/trustregions/optimize/inplace_loop.jl
index c42f02a..2b195ad 100644
--- a/src/optimize/trustregions/optimize/inplace_loop.jl
+++ b/src/optimize/trustregions/optimize/inplace_loop.jl
@@ -16,7 +16,7 @@ function solve(
     objvars;
     initial_Δ = 20.0,
 )
-    if !(mstyle(problem) === InPlace())
+    if !(mstyle(problem) === InPlace()) && !(approach.spsolve == Dogleg())
         throw(
             ErrorException("solve() not defined for OutOfPlace() with TrustRegion solvers"),
         )
@@ -113,9 +113,11 @@ function iterate!(
     scheme, subproblemsolver = modelscheme(approach), algorithm(approach)
     y, d, s = qnvars.y, qnvars.d, qnvars.s
     fx = fz
-    copyto!(x, z)
-    copyto!(∇fx, ∇fz)
-    spr = subproblemsolver(∇fx, B, Δk, p, scheme; abstol = 1e-10)
+
+    x = _copyto(mstyle(problem), x, z)
+    ∇fx = _copyto(mstyle(problem), ∇fx, ∇fz)
+
+    spr = subproblemsolver(∇fx, B, Δk, p, scheme, problem.mstyle; abstol = 1e-10)
     Δm = -spr.mz
 
     # Grab the model value, m. If m is zero, the solution, z, does not improve
@@ -128,7 +130,7 @@ function iterate!(
         # set flag to check for problems
     end
 
-    z = retract(problem, z, x, p)
+    z = retract(problem, z, x, spr.p)
     # Update before acceptance, to keep adding information about the hessian
     # even when the step is not "good" enough.
 
@@ -144,9 +146,9 @@ function iterate!(
     Δkp1, reject_step = update_trust_region(spr, R, p)
 
     if reject_step
-        z .= x
+        z = _copyto(mstyle(problem), z, x)
+        ∇fz = _copyto(mstyle(problem), ∇fz, ∇fx)
         fz = fx
-        ∇fz .= ∇fx
     end
     return (x = x, fx = fx, ∇fx = ∇fx, z = z, fz = fz, ∇fz = ∇fz, B = B, Pg = nothing),
     Δkp1,
diff --git a/src/quasinewton/approximations/newton.jl b/src/quasinewton/approximations/newton.jl
index d421d3f..90eb0f9 100644
--- a/src/quasinewton/approximations/newton.jl
+++ b/src/quasinewton/approximations/newton.jl
@@ -12,7 +12,7 @@ function default_newton_linsolve(B, g)
     B \ g
 end
 function default_newton_linsolve(d, B, g)
-    d .= (B \ g)
+    d = (B \ g)
 end
 
 function init_f∇fB(prob, scheme::Newton, ∇fz, B, x)
diff --git a/test/globalization/truncatedconjugategradient.jl b/test/globalization/truncatedconjugategradient.jl
index 3f8641c..1f9a8d3 100644
--- a/test/globalization/truncatedconjugategradient.jl
+++ b/test/globalization/truncatedconjugategradient.jl
@@ -5,29 +5,29 @@ using Test, NLSolvers, LinearAlgebra
     H = [0.3 0.0; 0.0 0.9]
     g = [0.2, 0.4]
 
-    m(g, H, 0.7, rand(2), 1)
+    m(g, H, 0.7, rand(2), 1, NLSolvers.InPlace())
 
     # Since -H is negative definite the solution is guaranteed
     # to be at the boundary unless g = 0
     for Δ in range(0, 100; step = 0.1)
-        @test norm(m(g, -H, Δ, rand(2), 1).p, 2) ≈ Δ
+        @test norm(m(g, -H, Δ, rand(2), 1, NLSolvers.InPlace()).p, 2) ≈ Δ
     end
 
     # Small gradient entries
     g .= [1e-12, 1e-9]
     for Δ in range(0, 100; step = 0.1)
-        @test norm(m(g, -H, Δ, rand(2), 1).p, 2) ≈ Δ
+        @test norm(m(g, -H, Δ, rand(2), 1, NLSolvers.InPlace()).p, 2) ≈ Δ
     end
 
     # Mixed gradient entries
     g .= [1e12, 1e-9]
     for Δ in range(0, 100; step = 0.1)
-        @test norm(m(g, -H, Δ, rand(2), 1).p, 2) ≈ Δ
+        @test norm(m(g, -H, Δ, rand(2), 1, NLSolvers.InPlace()).p, 2) ≈ Δ
     end
 
     # Zero case
     g = [0.0, 0.0]
     for Δ in range(0, 100; step = 0.1)
-        @test norm(m(g, -H, Δ, rand(2), 1).p, 2) == 0
+        @test norm(m(g, -H, Δ, rand(2), 1, NLSolvers.InPlace()).p, 2) == 0
     end
 end
diff --git a/test/nlsolve/interface.jl b/test/nlsolve/interface.jl
index 6996400..55b9a2c 100644
--- a/test/nlsolve/interface.jl
+++ b/test/nlsolve/interface.jl
@@ -341,9 +341,6 @@ end
     # Fc = zeros(6)
     # minimize!(lsqwrap, [100.0, 1.0], NelderMead(), OptimizationOptions())
 
-    # lsqwrap = NLSolvers.LsqWrapper(F, zeros(6), zeros(6,2))
-    # minimize!(MinProblem(;obj=lsqwrap,bounds=([0.0,0.0],[250.0,2.0])), [100.0, 1.0], ParticleSwarm(), OptimizationOptions())
-
     # function F(F, b)
     #   @. F = b[1]*(1 - exp(-b[2]*xdata)) - ydata
 
@@ -363,8 +360,6 @@ end
     # Fc = zeros(6)
 
 
-    # lsqwrap = NLSolvers.LsqWrapper(OnceDiffed(F), zeros(6), zeros(6,2))
-    # minimize!(lsqwrap, [100.0, 1.0], LineSearch(LBFGS()), OptimizationOptions())
 
 
     # #using Plots
@@ -467,4 +462,38 @@ end
 
     x0 = 0.3
     res = solve(prob, x0, LineSearch(Newton(), Backtracking()))
-end
\ No newline at end of file
+end
+
+function solve_static()
+    function F_rosenbrock_static(Fx, x)
+        Fx1 = 1 - x[1]
+        Fx2 = 10(x[2] - x[1]^2)
+        return @SVector([Fx1,Fx2])
+    end
+    function J_rosenbrock_static(Jx, x)
+        Jx11 = -1
+        Jx12 = 0
+        Jx21 = -20 * x[1]
+        Jx22 = 10
+        return @SArray([Jx11 Jx12; Jx21 Jx22])
+    end
+    function FJ_rosenbrock_static(Fx, Jx, x)
+        Fx = F_rosenbrock_static(Fx, x)
+        Jx = J_rosenbrock_static(Jx, x)
+        Fx, Jx
+    end
+    obj = NLSolvers.VectorObjective(
+        F_rosenbrock_static,
+        J_rosenbrock_static,
+        FJ_rosenbrock_static,
+        nothing,
+    )
+
+    prob_static =  NEqProblem(obj; inplace=false)
+    x0_static = @SVector([-1.2, 1.0])
+    res = solve(prob_static, x0_static, TrustRegion(Newton(), Dogleg()), NEqOptions())
+end
+
+solve_static()
+alloced = @allocated solve_static()
+@test alloced == 0
\ No newline at end of file
diff --git a/test/nlsolve/problems.jl b/test/nlsolve/problems.jl
index 4239fa2..988f31a 100644
--- a/test/nlsolve/problems.jl
+++ b/test/nlsolve/problems.jl
@@ -41,6 +41,34 @@ NLE_PROBS["rosenbrock"]["array"]["mutating"] = NLSolvers.VectorObjective(
 )
 
 
+NLE_PROBS["rosenbrock"]["static"] = Dict()
+function F_rosenbrock_static(Fx, x)
+    Fx1 = 1 - x[1]
+    Fx2 = 10(x[2] - x[1]^2)
+    return @SVector([Fx1,Fx2])
+end
+function J_rosenbrock_static(Jx, x)
+    Jx11 = -1
+    Jx12 = 0
+    Jx21 = -20 * x[1]
+    Jx22 = 10
+    return @SArray([Jx11 Jx12; Jx21 Jx22])
+end
+function FJ_rosenbrock_static(Fx, Jx, x)
+    F_rosenbrock_static(Fx, x)
+    J_rosenbrock_static(Jx, x)
+    Fx, Jx
+end
+
+NLE_PROBS["rosenbrock"]["static"]["x0"] = @SVector([-1.2, 1.0])
+NLE_PROBS["rosenbrock"]["static"]["mutating"] = NLSolvers.VectorObjective(
+    F_rosenbrock_static,
+    J_rosenbrock_static,
+    FJ_rosenbrock_static,
+    nothing,
+)
+
+
 function F_powell_singular!(x::Vector, Fx::Vector, Jx::Union{Nothing,Matrix} = nothing)
     if !(Fx isa Nothing)
         Fx[1] = x[1] + 10x[2]
diff --git a/test/runtests.jl b/test/runtests.jl
index cc12bb8..2eeb21d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,11 +7,12 @@ using Printf
 using LinearAlgebra: norm, I
 import Random
 Random.seed!(41234)
+
+include("nlsolve/problems.jl")
+include("nlsolve/interface.jl")
 include("optimize/problems.jl")
 include("optimize/interface.jl")
 include("optimize/preconditioning.jl")
 include("optimize/complex.jl")
-include("nlsolve/problems.jl")
-include("nlsolve/interface.jl")
 include("lsqfit/interface.jl")
 include("globalization/runtests.jl")

From 8160bfa45df0b93a1b5beae02c3b17a65e53f4cc Mon Sep 17 00:00:00 2001
From: Patrick Kofod Mogensen <patrick.mogensen@gmail.com>
Date: Sat, 20 Jan 2024 17:35:24 +0100
Subject: [PATCH 2/3] Small fixes

---
 src/NLSolvers.jl                               |  4 ++--
 src/globalization/linesearches/backtracking.jl |  2 +-
 .../trs_solvers/solvers/Dogleg.jl              |  2 +-
 src/objectives.jl                              | 12 ++++++++++--
 .../trustregions/optimize/inplace_loop.jl      | 18 ++++++++++++++----
 5 files changed, 28 insertions(+), 10 deletions(-)

diff --git a/src/NLSolvers.jl b/src/NLSolvers.jl
index dda7ff7..72617c0 100644
--- a/src/NLSolvers.jl
+++ b/src/NLSolvers.jl
@@ -192,10 +192,10 @@ end
 function negate(::OutOfPlace, A::AbstractArray)
     -A
 end
-function scale(::InPlace, B, A, m)
+function _scale(::InPlace, B, A, m)
     B .= A .* m
 end
-function scale(::OutOfPlace, B, A, m)
+function _scale(::OutOfPlace, B, A, m)
     A .* m
 end
 function _copyto(::InPlace, z, x)
diff --git a/src/globalization/linesearches/backtracking.jl b/src/globalization/linesearches/backtracking.jl
index 2e64cb2..5de80e3 100644
--- a/src/globalization/linesearches/backtracking.jl
+++ b/src/globalization/linesearches/backtracking.jl
@@ -158,7 +158,7 @@ function find_steplength(mstyle, ls::Backtracking, φ::T, λ) where {T}
         Tf(ls.ratio), Tf(ls.decrease), ls.maxiter, ls.verbose
 
     #== factor in Armijo condition ==#
-    t = -decrease * dφ0
+    t = decrease * dφ0
 
     iter, α, β = 0, λ, λ # iteration variables
     f_α = φ(α) # initial function value
diff --git a/src/globalization/trs_solvers/solvers/Dogleg.jl b/src/globalization/trs_solvers/solvers/Dogleg.jl
index 134d001..d61e916 100644
--- a/src/globalization/trs_solvers/solvers/Dogleg.jl
+++ b/src/globalization/trs_solvers/solvers/Dogleg.jl
@@ -33,7 +33,7 @@ function (dogleg::Dogleg)(∇f, H, Δ, p, scheme, mstyle; abstol = 1e-10, maxite
     if norm_d_cauchy ≥ Δ
         shrink = Δ / norm_d_cauchy # inv(Δ/norm_d_cauchy) puts it on the border
 
-        p = scale(mstyle, p, d_cauchy, shrink)
+        p = _scale(mstyle, p, d_cauchy, shrink)
         interior = false
     else
         # Else, calculate (Quasi-)Newton step. If this is interior, then take the
diff --git a/src/objectives.jl b/src/objectives.jl
index 41f2037..8390e64 100644
--- a/src/objectives.jl
+++ b/src/objectives.jl
@@ -19,8 +19,16 @@ struct ScalarObjective{Tf,Tg,Tfg,Tfgh,Th,Thv,Tbf,P}
     batched_f::Tbf
     param::P
 end
-ScalarObjective(; f = nothing, g = nothing, fg = nothing, fgh = nothing, h = nothing) =
-    ScalarObjective(f, g, fg, fgh, h, nothing, nothing, nothing)
+ScalarObjective(;
+     f = nothing,
+     g = nothing,
+     fg = nothing,
+     fgh = nothing,
+     h = nothing,
+     hv = nothing,
+     batched_f = nothing,
+     param = nothing,
+ ) = ScalarObjective(f, g, fg, fgh, h, hv, batched_f, param)
 has_param(so::ScalarObjective) = so.param === nothing ? false : true
 function value(so::ScalarObjective, x)
     if has_param(so)
diff --git a/src/optimize/trustregions/optimize/inplace_loop.jl b/src/optimize/trustregions/optimize/inplace_loop.jl
index 2b195ad..f1b7eb8 100644
--- a/src/optimize/trustregions/optimize/inplace_loop.jl
+++ b/src/optimize/trustregions/optimize/inplace_loop.jl
@@ -49,7 +49,7 @@ function solve(
     qnvars = QNVars(objvars.z, objvars.z)
     p = copy(objvars.x)
 
-    objvars, Δkp1, reject =
+    objvars, Δkp1, reject, qnvars =
         iterate!(p, objvars, Δk, approach, problem, options, qnvars, false)
 
     iter = 1
@@ -57,7 +57,7 @@ function solve(
     is_converged = converged(approach, objvars, ∇f0, options, reject, Δkp1)
     while iter <= options.maxiter && !any(is_converged)
         iter += 1
-        objvars, Δkp1, reject =
+        objvars, Δkp1, reject, qnvars =
             iterate!(p, objvars, Δkp1, approach, problem, options, qnvars, false)
 
         # Check for convergence
@@ -126,7 +126,7 @@ function iterate!(
     # tive value, we may conclude that "something" is wrong. We might be at a
     # ridge (positive-indefinite case) for example, or the scaling of the model
     # is such that we cannot satisfy ||∇f|| < tol.
-    if abs(spr.mz) < eps(T)
+    if abs(spr.mz) < eps(spr.mz)
         # set flag to check for problems
     end
 
@@ -134,6 +134,10 @@ function iterate!(
     # Update before acceptance, to keep adding information about the hessian
     # even when the step is not "good" enough.
 
+    y, d, s = qnvars.y, qnvars.d, qnvars.s
+    fx = fz
+    # Should build a good code for picking update model.
+
     fz, ∇fz, B, s, y = update_obj!(problem, spr.p, y, ∇fx, z, ∇fz, B, scheme, scale)
 
     # Δf is often called ared or Ared for actual reduction. I prefer "change in"
@@ -149,10 +153,16 @@ function iterate!(
         z = _copyto(mstyle(problem), z, x)
         ∇fz = _copyto(mstyle(problem), ∇fz, ∇fx)
         fz = fx
+        # This is correct because 
+        s = _scale(mstyle(problem), s, s, 0) # z - x
+        y = _scale(mstyle(problem), y, y, 0) # ∇fz - ∇fx
+        # and will cause quasinewton updates to not update
+        # this seems incorrect as it's already updated, should hold off here
     end
     return (x = x, fx = fx, ∇fx = ∇fx, z = z, fz = fz, ∇fz = ∇fz, B = B, Pg = nothing),
     Δkp1,
-    reject_step
+    reject_step,
+    QNVars(d, s, y)
 end
 
 function update_trust_region(spr, R, p)

From 5c5f598cbb27852129b37013dfaa757694e2c6cb Mon Sep 17 00:00:00 2001
From: Patrick Kofod Mogensen <patrick.mogensen@gmail.com>
Date: Sat, 20 Jan 2024 17:49:35 +0100
Subject: [PATCH 3/3] Update Project.toml

---
 Project.toml | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index 04d17bb..3a038f1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,14 +4,10 @@ authors = ["Patrick Kofod Mogensen <patrick.mogensen@gmail.com>"]
 version = "0.4.0"
 
 [deps]
-Cthulhu = "f68482b8-f384-11e8-15f7-abe071a5a75f"
-DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 PositiveFactorizations = "85a6dd25-e78a-55b7-8502-1745935b8125"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]