Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Add wrapper for FixedPointAcceleration
- Loading branch information
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,3 +1,56 @@ | ||
| module NonlinearSolveFixedPointAccelerationExt | ||
|
|
||
| using NonlinearSolve, FixedPointAcceleration, DiffEqBase, SciMLBase | ||
|
|
||
| function SciMLBase.__solve(prob::NonlinearProblem, alg::FixedPointAccelerationJL, args...; | ||
| abstol = nothing, maxiters = 1000, alias_u0::Bool = false, | ||
| show_trace::Val{PrintReports} = Val(false), termination_condition = nothing, | ||
| kwargs...) where {PrintReports} | ||
| @assert (termination_condition === | ||
| nothing)||(termination_condition isa AbsNormTerminationMode) "SpeedMappingJL does not support termination conditions!" | ||
|
|
||
| u0 = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0) | ||
| u_size = size(u0) | ||
| T = eltype(u0) | ||
| iip = isinplace(prob) | ||
| p = prob.p | ||
|
|
||
| if !iip && prob.u0 isa Number | ||
| # FixedPointAcceleration makes the scalar problem into a vector problem | ||
| f = (u) -> [prob.f(u[1], p) .+ u[1]] | ||
| elseif !iip && prob.u0 isa AbstractVector | ||
| f = (u) -> (prob.f(u, p) .+ u) | ||
| elseif !iip && prob.u0 isa AbstractArray | ||
| f = (u) -> vec(prob.f(reshape(u, u_size), p) .+ u) | ||
| elseif iip && prob.u0 isa AbstractVector | ||
| du = similar(u0) | ||
| f = (u) -> (prob.f(du, u, p); du .+ u) | ||
| else | ||
| du = similar(u0) | ||
| f = (u) -> (prob.f(du, reshape(u, u_size), p); vec(du) .+ u) | ||
| end | ||
|
|
||
| tol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol | ||
|
|
||
| sol = fixed_point(f, NonlinearSolve._vec(u0); Algorithm = alg.algorithm, | ||
| ConvergenceMetricThreshold = tol, MaxIter = maxiters, MaxM = alg.m, | ||
| ExtrapolationPeriod = alg.extrapolation_period, Dampening = alg.dampening, | ||
| PrintReports, ReplaceInvalids = alg.replace_invalids, | ||
| ConditionNumberThreshold = alg.condition_number_threshold, quiet_errors = true) | ||
|
|
||
| res = prob.u0 isa Number ? first(sol.FixedPoint_) : sol.FixedPoint_ | ||
| if res === missing | ||
| resid = NonlinearSolve.evaluate_f(prob, u0) | ||
| res = u0 | ||
| converged = false | ||
| else | ||
| resid = NonlinearSolve.evaluate_f(prob, res) | ||
| converged = maximum(abs, resid) ≤ tol | ||
| end | ||
| return SciMLBase.build_solution(prob, alg, res, resid; | ||
| retcode = converged ? ReturnCode.Success : ReturnCode.Failure, | ||
| stats = SciMLBase.NLStats(sol.Iterations_, 0, 0, 0, sol.Iterations_), | ||
| original = sol) | ||
| end | ||
|
|
||
| end | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,68 @@ | ||
| using NonlinearSolve, FixedPointAcceleration, LinearAlgebra, Test | ||
|
|
||
| # Simple Scalar Problem | ||
| @testset "Simple Scalar Problem" begin | ||
| f1(x, p) = cos(x) - x | ||
| prob = NonlinearProblem(f1, 1.1) | ||
|
|
||
| for alg in (:Anderson, :MPE, :RRE, :VEA, :SEA, :Simple, :Aitken, :Newton) | ||
| @test abs(solve(prob, FixedPointAccelerationJL()).resid) ≤ 1e-10 | ||
| end | ||
| end | ||
|
|
||
| # Simple Vector Problem | ||
| @testset "Simple Vector Problem" begin | ||
| f2(x, p) = cos.(x) .- x | ||
| prob = NonlinearProblem(f2, [1.1, 1.1]) | ||
|
|
||
| for alg in (:Anderson, :MPE, :RRE, :VEA, :SEA, :Simple, :Aitken, :Newton) | ||
| @test maximum(abs.(solve(prob, FixedPointAccelerationJL()).resid)) ≤ 1e-10 | ||
| end | ||
| end | ||
|
|
||
| # Fixed Point for Power Method | ||
| # Taken from https://github.com/NicolasL-S/SpeedMapping.jl/blob/95951db8f8a4457093090e18802ad382db1c76da/test/runtests.jl | ||
| @testset "Power Method" begin | ||
| C = [1 2 3; 4 5 6; 7 8 9] | ||
| A = C + C' | ||
| B = Hermitian(ones(10) * ones(10)' .* im + Diagonal(1:10)) | ||
|
|
||
| function power_method!(du, u, A) | ||
| mul!(du, A, u) | ||
| du ./= norm(du, Inf) | ||
| du .-= u # Convert to a root finding problem | ||
| return nothing | ||
| end | ||
|
|
||
| prob = NonlinearProblem(power_method!, ones(3), A) | ||
|
|
||
| for alg in (:Anderson, :MPE, :RRE, :VEA, :SEA, :Simple, :Aitken, :Newton) | ||
| sol = solve(prob, FixedPointAccelerationJL(; algorithm = alg)) | ||
| if SciMLBase.successful_retcode(sol) | ||
| @test sol.u' * A[:, 3] ≈ 32.916472867168096 | ||
| else | ||
| @warn "Power Method failed for FixedPointAccelerationJL(; algorithm = $alg)" | ||
| @test_broken sol.u' * A[:, 3] ≈ 32.916472867168096 | ||
| end | ||
| end | ||
|
|
||
| # Non vector inputs | ||
| function power_method_nonvec!(du, u, A) | ||
| mul!(vec(du), A, vec(u)) | ||
| du ./= norm(du, Inf) | ||
| du .-= u # Convert to a root finding problem | ||
| return nothing | ||
| end | ||
|
|
||
| prob = NonlinearProblem(power_method_nonvec!, ones(1, 3, 1), A) | ||
|
|
||
| for alg in (:Anderson, :MPE, :RRE, :VEA, :SEA, :Simple, :Aitken, :Newton) | ||
| sol = solve(prob, FixedPointAccelerationJL(; algorithm = alg)) | ||
| if SciMLBase.successful_retcode(sol) | ||
| @test sol.u' * A[:, 3] ≈ 32.916472867168096 | ||
| else | ||
| @warn "Power Method failed for FixedPointAccelerationJL(; algorithm = $alg)" | ||
| @test_broken sol.u' * A[:, 3] ≈ 32.916472867168096 | ||
| end | ||
| end | ||
| end |