diff --git a/Project.toml b/Project.toml index e4edff8..91a3a44 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,16 @@ +authors = ["gszep "] name = "BifurcationInference" uuid = "7fe238d6-d31e-4646-aa16-9d8429fd6da8" -authors = ["gszep "] -version = "0.1.3" +version = "0.1.4" [deps] BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" +Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" InvertedIndices = "41ab1584-1d38-5bbf-9106-f11c6c58b48f" +JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" @@ -15,18 +18,20 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] -BifurcationKit = "0.1.5, 0.1" -Flux = "0.12" +BifurcationKit = "0.3" +Flux = "0.14" ForwardDiff = "0.10" InvertedIndices = "1" LaTeXStrings = "1" Parameters = "0.12" Plots = "1" -Setfield = "0.7, 0.8" -SpecialFunctions = "1.5, 2" -StaticArrays = "1.2" +Setfield = "1" +SpecialFunctions = "2" +StaticArrays = "1" julia = "1" [extras] diff --git a/src/BifurcationInference.jl b/src/BifurcationInference.jl index a516e7f..ef52a92 100644 --- a/src/BifurcationInference.jl +++ b/src/BifurcationInference.jl @@ -1,140 +1,154 @@ module BifurcationInference - using BifurcationKit: ContIterable, newton, ContinuationPar, NewtonPar, DeflationOperator - using BifurcationKit: BorderedArray, AbstractLinearSolver, AbstractEigenSolver, BorderingBLS - using BifurcationKit: ContState, detectBifucation +using BifurcationKit: BifurcationProblem, re_make, PALC, ContIterable, newton, ContinuationPar, NewtonPar, DeflationOperator +using BifurcationKit: BorderedArray, AbstractLinearSolver, AbstractEigenSolver, BorderingBLS +using BifurcationKit: ContState, detect_bifurcation - using ForwardDiff: Dual,tagtype,derivative,gradient,jacobian - using Flux: Momentum,update! +using ForwardDiff: Dual, tagtype, derivative, gradient, jacobian +using Flux: Momentum, update! - using Setfield: @lens,@set,setproperties - using Parameters: @unpack +using Setfield: @lens, @set, setproperties +using Parameters: @unpack - using InvertedIndices: Not - using LinearAlgebra, StaticArrays +using InvertedIndices: Not +using LinearAlgebra, StaticArrays - include("Structures.jl") - include("Utils.jl") +include("Structures.jl") +include("Utils.jl") - include("Objectives.jl") - include("Gradients.jl") - include("Plots.jl") +include("Objectives.jl") +include("Gradients.jl") +include("Plots.jl") - export plot,@unpack,BorderedArray,SizedVector - export StateSpace,deflationContinuation,train! - export getParameters,loss,∇loss,norm +export plot, @unpack, BorderedArray, SizedVector +export StateSpace, deflationContinuation, train! +export getParameters, loss, ∇loss, norm - """ root finding with newton deflation method""" - function findRoots!( f::Function, J::Function, roots::AbstractVector{<:AbstractVector}, - parameters::NamedTuple, hyperparameters::ContinuationPar; - maxRoots::Int = 3, maxIter::Int=500, verbosity=0 ) +""" root finding with newton deflation method""" +function findRoots!(f::Function, J::Function, roots::AbstractVector{<:AbstractVector}, + parameters::NamedTuple, hyperparameters::ContinuationPar; + maxRoots::Int=3, max_iterations::Int=500, verbosity=0) - hyperparameters = @set hyperparameters.newtonOptions = setproperties( - hyperparameters.newtonOptions; maxIter = maxIter, verbose = verbosity ) + hyperparameters = @set hyperparameters.newton_options = setproperties( + hyperparameters.newton_options; max_iterations=max_iterations, verbose=verbosity) - # search for roots across parameter range - pRange = range(hyperparameters.pMin,hyperparameters.pMax,length=length(roots)) - roots .= findRoots.( Ref(f), Ref(J), roots, pRange, Ref(parameters), Ref(hyperparameters); maxRoots=maxRoots ) - end + # search for roots across parameter range + pRange = range(hyperparameters.p_min, hyperparameters.p_max, length=length(roots)) + roots .= findRoots.(Ref(f), Ref(J), roots, pRange, Ref(parameters), Ref(hyperparameters); maxRoots=maxRoots) +end + +function findRoots(f::Function, J::Function, roots::AbstractVector{V}, p::T, + parameters::NamedTuple, hyperparameters::ContinuationPar{T,S,E}; maxRoots::Int=3, converged=false +) where {T<:Number,V<:AbstractVector{T},S<:AbstractLinearSolver,E<:AbstractEigenSolver} + + Zero = zero(first(roots)) + inf = Zero .+ Inf - function findRoots( f::Function, J::Function, roots::AbstractVector{V}, p::T, - parameters::NamedTuple, hyperparameters::ContinuationPar{T, S, E}; maxRoots::Int = 3, converged = false - ) where { T<:Number, V<:AbstractVector{T}, S<:AbstractLinearSolver, E<:AbstractEigenSolver } + # search for roots at specific parameter value + deflation = DeflationOperator(one(T), dot, one(T), [inf]) # dummy deflation at infinity + parameters = @set parameters.p = p - Zero = zero(first(roots)) - inf = Zero .+ Inf + problem = BifurcationProblem(f, roots[begin] .+ hyperparameters.ds, parameters; J=J) + solution = newton(problem, deflation, hyperparameters.newton_options) - # search for roots at specific parameter value - deflation = DeflationOperator(one(T), dot, one(T), [inf] ) # dummy deflation at infinity - parameters = @set parameters.p = p + for u ∈ roots # update existing roots + solution = newton(re_make(problem; u0=u .+ hyperparameters.ds), deflation, hyperparameters.newton_options) - for u ∈ roots # update existing roots - u, residual, converged, niter = newton( f, J, u.+hyperparameters.ds, parameters, - hyperparameters.newtonOptions, deflation) + i = 0 + while any(isnan.(solution.residuals)) & (i < hyperparameters.newton_options.max_iterations) + u .= randn(length(u)) + + solution = newton(re_make(problem; u0=u .+ hyperparameters.ds), deflation, hyperparameters.newton_options) + i += 1 + end - i = 0 - while any(isnan.(residual)) & (i root ≠ inf, deflation.roots) # remove dummy deflation at infinity + @assert(length(deflation.roots) > 0, "No roots f(u,p)=0 found at p = $(parameters.p), θ = $(parameters.θ); try increasing max_iterations") + return deflation.roots +end - @assert( !any(isnan.(residual)), "f(u,p) = $(residual[end]) at u = $u, p = $(parameters.p), θ = $(parameters.θ)") - if converged push!(deflation,u) else break end +""" deflation continuation method """ +function deflationContinuation(f::Function, roots::AbstractVector{<:AbstractVector{V}}, + parameters::NamedTuple, hyperparameters::ContinuationPar{T,S,E}; + maxRoots::Int=3, max_iterations::Int=500, resolution=400, verbosity=0, kwargs... +) where {T<:Number,V<:AbstractVector{T},S<:AbstractLinearSolver,E<:AbstractEigenSolver} + + max_iterationsContinuation, ds = hyperparameters.newton_options.max_iterations, hyperparameters.ds + J(u, p) = jacobian(x -> f(x, p), u) + + findRoots!(f, J, roots, parameters, hyperparameters; maxRoots=maxRoots, max_iterations=max_iterations, verbosity=verbosity) + pRange = range(hyperparameters.p_min, hyperparameters.p_max, length=length(roots)) + intervals = ([zero(T), step(pRange)], [-step(pRange), zero(T)]) + + branches = Vector{Branch{V,T}}() + problem = BifurcationProblem(f, roots[begin][begin], parameters, (@lens _.p); J=J) + + hyperparameters = @set hyperparameters.newton_options.max_iterations = max_iterationsContinuation + linsolver = BorderingBLS(hyperparameters.newton_options.linsolver) + algorithm = PALC() + + for (i, us) ∈ enumerate(roots) + for u ∈ us # perform continuation for each root + + # forwards and backwards branches + for (p_min, p_max) ∈ intervals + + hyperparameters = setproperties(hyperparameters; + p_min=pRange[i] + p_min, p_max=pRange[i] + p_max, + ds=sign(hyperparameters.ds) * ds) + + # main continuation method + branch = Branch{V,T}() + parameters = @set parameters.p = pRange[i] + hyperparameters.ds + + try + iterator = ContIterable(re_make(problem; u0=u, params=parameters), algorithm, hyperparameters; verbosity=verbosity) + for state ∈ iterator + push!(branch, state) + end + + midpoint = sum(s -> s.z.p, branch) / length(branch) + if minimum(pRange) < midpoint < maximum(pRange) + push!(branches, p_min < 0 ? reverse(branch) : branch) + end + + catch error + printstyled(color=:red, "Continuation Error at f(u,p)=$(f(u,parameters))\nu=$u, p=$(parameters.p), θ=$(parameters.θ)\n") + rethrow(error) + end + hyperparameters = @set hyperparameters.ds = -hyperparameters.ds + end end + end - u = Zero - if converged || length(deflation)==1 # search for new roots - while length(deflation)-1 < maxRoots - - u, residual, converged, niter = newton( f, J, u.+hyperparameters.ds, parameters, - hyperparameters.newtonOptions, deflation) - - # make sure new roots are different from existing - if any( isapprox.( Ref(u), deflation.roots, atol=2*hyperparameters.ds ) ) break end - if converged push!(deflation,u) else break end - end - end - - filter!( root->root≠inf, deflation.roots ) # remove dummy deflation at infinity - @assert( length(deflation.roots)>0, "No roots f(u,p)=0 found at p = $(parameters.p), θ = $(parameters.θ); try increasing maxIter") - return deflation.roots - end - - """ deflation continuation method """ - function deflationContinuation( f::Function, roots::AbstractVector{<:AbstractVector{V}}, - parameters::NamedTuple, hyperparameters::ContinuationPar{T, S, E}; - maxRoots::Int = 3, maxIter::Int=500, resolution=400, verbosity=0, kwargs... - ) where {T<:Number, V<:AbstractVector{T}, S<:AbstractLinearSolver, E<:AbstractEigenSolver} - - maxIterContinuation,ds = hyperparameters.newtonOptions.maxIter,hyperparameters.ds - J(u,p) = jacobian(x->f(x,p),u) - - findRoots!( f, J, roots, parameters, hyperparameters; maxRoots=maxRoots, maxIter=maxIter, verbosity=verbosity) - pRange = range(hyperparameters.pMin,hyperparameters.pMax,length=length(roots)) - intervals = ([zero(T),step(pRange)],[-step(pRange),zero(T)]) - - branches = Vector{Branch{V,T}}() - hyperparameters = @set hyperparameters.newtonOptions.maxIter = maxIterContinuation - linsolver = BorderingBLS(hyperparameters.newtonOptions.linsolver) - - for (i,us) ∈ enumerate(roots) - for u ∈ us # perform continuation for each root - - # forwards and backwards branches - for (pMin,pMax) ∈ intervals - - hyperparameters = setproperties(hyperparameters; - pMin=pRange[i]+pMin, pMax=pRange[i]+pMax, - ds=sign(hyperparameters.ds)*ds) - - # main continuation method - branch = Branch{V,T}() - parameters = @set parameters.p = pRange[i]+hyperparameters.ds - - try - iterator = ContIterable( f, J, u, parameters, (@lens _.p), hyperparameters, linsolver, verbosity=verbosity) - for state ∈ iterator - push!(branch,state) - end - - midpoint = sum( s -> s.z.p, branch ) / length(branch) - if minimum(pRange) < midpoint < maximum(pRange) - push!(branches,pMin < 0 ? reverse(branch) : branch) end - - catch error - printstyled(color=:red,"Continuation Error at f(u,p)=$(f(u,parameters))\nu=$u, p=$(parameters.p), θ=$(parameters.θ)\n") - rethrow(error) - end - hyperparameters = @set hyperparameters.ds = -hyperparameters.ds - end - end - end - - hyperparameters = @set hyperparameters.ds = ds - updateParameters!(hyperparameters,branches;resolution=resolution) - return unique(branches; atol=10*hyperparameters.ds) - end + hyperparameters = @set hyperparameters.ds = ds + updateParameters!(hyperparameters, branches; resolution=resolution) + return unique(branches; atol=10 * hyperparameters.ds) end +end # module diff --git a/src/Gradients.jl b/src/Gradients.jl index 78e4efe..54ca02c 100644 --- a/src/Gradients.jl +++ b/src/Gradients.jl @@ -1,71 +1,71 @@ ################################################################################ -function ∇loss( F::Function, branches::AbstractVector{<:Branch}, θ::AbstractVector, targets::StateSpace; kwargs...) +function ∇loss(F::Function, branches::AbstractVector{<:Branch}, θ::AbstractVector, targets::StateSpace; kwargs...) - predictions = unique([ s.z for branch ∈ branches for s ∈ branch if s.bif ]; atol=2*step(targets.parameter)) - λ = length(targets.targets)-length(predictions) + predictions = unique([s.z for branch ∈ branches for s ∈ branch if s.bif]; atol=2 * step(targets.parameter)) + λ = length(targets.targets) - length(predictions) - if λ≠0 - Φ,∇Φ = measure(F,branches,θ,targets), ∇measure(F,branches,θ,targets) - return errors(predictions,targets) - λ*log(Φ), ∇errors(F,predictions,θ,targets) - λ*∇Φ/Φ - else - return errors(predictions,targets), ∇errors(F,predictions,θ,targets) - end + if λ ≠ 0 + Φ, ∇Φ = measure(F, branches, θ, targets), ∇measure(F, branches, θ, targets) + return errors(predictions, targets) - λ * log(Φ), ∇errors(F, predictions, θ, targets) - λ * ∇Φ / Φ + else + return errors(predictions, targets), ∇errors(F, predictions, θ, targets) + end end ################################################################################ -function ∇errors( F::Function, predictions::AbstractVector{<:BorderedArray}, θ::AbstractVector, targets::StateSpace) - return mean( p′-> mean( z->(z.p-p′)^2, predictions; type=:geometric )*mean( z-> 2velocity(F,z,θ)/(z.p-p′), predictions; type=:arithmetic ), targets.targets; type=:arithmetic ) +function ∇errors(F::Function, predictions::AbstractVector{<:BorderedArray}, θ::AbstractVector, targets::StateSpace) + return mean(p′ -> mean(z -> (z.p - p′)^2, predictions; type=:geometric) * mean(z -> 2velocity(F, z, θ) / (z.p - p′), predictions; type=:arithmetic), targets.targets; type=:arithmetic) end -function ∇measure( F::Function, z::BorderedArray, θ::AbstractVector, targets::StateSpace; newtonOptions=NewtonPar(verbose=false,maxIter=800,tol=1e-6) ) - ∂implicit, _, _ = newtonOptions.linsolver( -∂Fz(F,z,θ)', gradient(z->measure(F,z,θ,targets),z) ) - return gradient( θ -> measure(F,z,θ,targets) + F(z,θ)'∂implicit , θ ) + measure(F,z,θ,targets)*∇region(F,z,θ) +function ∇measure(F::Function, z::BorderedArray, θ::AbstractVector, targets::StateSpace; newton_options=NewtonPar(verbose=false, max_iterations=800, tol=1e-6)) + ∂implicit, _, _ = newton_options.linsolver(-∂Fz(F, z, θ)', gradient(z -> measure(F, z, θ, targets), z)) + return gradient(θ -> measure(F, z, θ, targets) + F(z, θ)'∂implicit, θ) + measure(F, z, θ, targets) * ∇region(F, z, θ) end ########################################################################### -function ∇measure( F::Function, branch::Branch, θ::AbstractVector, targets::StateSpace ) - return sum( s -> ∇measure(F,s.z,θ,targets)*s.ds, branch ) +function ∇measure(F::Function, branch::Branch, θ::AbstractVector, targets::StateSpace) + return sum(s -> ∇measure(F, s.z, θ, targets) * s.ds, branch) end -function ∇measure( F::Function, branches::AbstractVector{<:Branch}, θ::AbstractVector, targets::StateSpace ) - return sum( branch -> ∇measure(F,branch,θ,targets), branches ) +function ∇measure(F::Function, branches::AbstractVector{<:Branch}, θ::AbstractVector, targets::StateSpace) + return sum(branch -> ∇measure(F, branch, θ, targets), branches) end ############################################## gradient term due to changing integration region dz -deformation( F::Function, z::BorderedArray, θ::AbstractVector ) = -∂Fz(F,z,θ)\∂Fθ(F,z,θ) # todo@(gszep) matrix inverse is computational bottleneck -function ∇region( F::Function, z::BorderedArray, θ::AbstractVector ) +deformation(F::Function, z::BorderedArray, θ::AbstractVector) = -∂Fz(F, z, θ) \ ∂Fθ(F, z, θ) # todo@(gszep) matrix inverse is computational bottleneck +function ∇region(F::Function, z::BorderedArray, θ::AbstractVector) - ∇deformation = reshape( jacobian(z->deformation(F,z,θ),z), length(z),length(θ),length(z) ) - tangent = tangent_field(F,z,θ) - t = [tangent.u;tangent.p] + ∇deformation = reshape(jacobian(z -> deformation(F, z, θ), z), length(z), length(θ), length(z)) + tangent = tangent_field(F, z, θ) + t = [tangent.u; tangent.p] - ∇region = similar(θ) - for k ∈ 1:length(θ) - ∇region[k] = t'∇deformation[:,k,:]t - end + ∇region = similar(θ) + for k ∈ 1:length(θ) + ∇region[k] = t'∇deformation[:, k, :]t + end - return ∇region + return ∇region end ########################################################################### jacobians # statespace jacobian -∂Fu(F::Function,z::BorderedArray,θ::AbstractVector) = ∂Fz(F,z,θ)[:,Not(end)] +∂Fu(F::Function, z::BorderedArray, θ::AbstractVector) = ∂Fz(F, z, θ)[:, Not(end)] # parameter jacobian -∂Fθ(F::Function,z::BorderedArray,θ::AbstractVector) = jacobian( θ->F(z,θ), θ ) +∂Fθ(F::Function, z::BorderedArray, θ::AbstractVector) = jacobian(θ -> F(z, θ), θ) # augmented jacobian -∂Fz(F::Function,z::BorderedArray,θ::AbstractVector) = jacobian( z->F(z,θ), z ) +∂Fz(F::Function, z::BorderedArray, θ::AbstractVector) = jacobian(z -> F(z, θ), z) ############################################################# autodiff wrappers for BorderedArray -import ForwardDiff: derivative,gradient,jacobian -derivative(f,x::BorderedArray,s::BorderedArray) = derivative( λ -> f(x + λ*s), 0.0 ) +import ForwardDiff: derivative, gradient, jacobian +derivative(f, x::BorderedArray, s::BorderedArray) = derivative(λ -> f(x + λ * s), 0.0) -function gradient( f, z::BorderedArray ) - return gradient( z -> f( BorderedArray(z[Not(end)],z[end]) ), [z.u; z.p] ) +function gradient(f, z::BorderedArray) + return gradient(z -> f(BorderedArray(z[Not(end)], z[end])), [z.u; z.p]) end -function jacobian( f, z::BorderedArray ) - return jacobian( z -> f( BorderedArray(z[Not(end)],z[end]) ), [z.u; z.p] ) +function jacobian(f, z::BorderedArray) + return jacobian(z -> f(BorderedArray(z[Not(end)], z[end])), [z.u; z.p]) end \ No newline at end of file diff --git a/src/Objectives.jl b/src/Objectives.jl index c221385..e32c3c3 100644 --- a/src/Objectives.jl +++ b/src/Objectives.jl @@ -1,26 +1,26 @@ ################################################################################ -function loss( F::Function, branches::AbstractVector{<:Branch}, θ::AbstractVector, targets::StateSpace; kwargs...) +function loss(F::Function, branches::AbstractVector{<:Branch}, θ::AbstractVector, targets::StateSpace; kwargs...) - predictions = unique([ s.z for branch ∈ branches for s ∈ branch if s.bif ]; atol=2*step(targets.parameter)) - λ = length(targets.targets)-length(predictions) + predictions = unique([s.z for branch ∈ branches for s ∈ branch if s.bif]; atol=2 * step(targets.parameter)) + λ = length(targets.targets) - length(predictions) - if λ≠0 - Φ = measure(F,branches,θ,targets) - return errors(predictions,targets) - λ*log(Φ) - else - return errors(predictions,targets) - end + if λ ≠ 0 + Φ = measure(F, branches, θ, targets) + return errors(predictions, targets) - λ * log(Φ) + else + return errors(predictions, targets) + end end ################################################################################ -function errors( predictions::AbstractVector{<:BorderedArray}, targets::StateSpace) - return mean( p′-> mean( z->(z.p-p′)^2, predictions; type=:geometric ), targets.targets; type=:arithmetic ) +function errors(predictions::AbstractVector{<:BorderedArray}, targets::StateSpace) + return mean(p′ -> mean(z -> (z.p - p′)^2, predictions; type=:geometric), targets.targets; type=:arithmetic) end realeigvals(J) = real(eigen(J).values) -function measure( F::Function, z::BorderedArray, θ::AbstractVector, targets::StateSpace ) - λ,dλ = realeigvals(∂Fu(F,z,θ)), derivative( z -> realeigvals(∂Fu(F,z,θ)), z, tangent_field(F,z,θ) ) - return window_function(targets.parameter,z) * mapreduce( (λ,dλ) -> 2 / ( 2 + abs(sinh(2λ)/dλ) ), +, λ, dλ ) +function measure(F::Function, z::BorderedArray, θ::AbstractVector, targets::StateSpace) + λ, dλ = realeigvals(∂Fu(F, z, θ)), derivative(z -> realeigvals(∂Fu(F, z, θ)), z, tangent_field(F, z, θ)) + return window_function(targets.parameter, z) * mapreduce((λ, dλ) -> 2 / (2 + abs(sinh(2λ) / dλ)), +, λ, dλ) end # Giles, M. 2008. An extended collection of matrix derivative results for forward and reverse mode automatic differentiation @@ -39,13 +39,13 @@ end function eigen(A::StridedMatrix{Dual{T,V,N}}) where {T,V,N} Av = map(a -> a.value, A) - λ,U = eigen(Av) + λ, U = eigen(Av) dΛ = U \ A * U dλ = diag(dΛ) F = similar(Av, eltype(λ)) - for i ∈ axes(A,1), j ∈ axes(A,2) + for i ∈ axes(A, 1), j ∈ axes(A, 2) if i == j F[i, j] = 0 else @@ -61,30 +61,30 @@ function eigen(A::StridedMatrix{Dual{T,V,N}}) where {T,V,N} end ################################################################################ -function measure( F::Function, branch::Branch, θ::AbstractVector, targets::StateSpace ) - return sum( s -> measure(F,s.z,θ,targets)*s.ds, branch ) +function measure(F::Function, branch::Branch, θ::AbstractVector, targets::StateSpace) + return sum(s -> measure(F, s.z, θ, targets) * s.ds, branch) end -function measure( F::Function, branches::AbstractVector{<:Branch}, θ::AbstractVector, targets::StateSpace ) - return sum( branch -> measure(F,branch,θ,targets), branches ) +function measure(F::Function, branches::AbstractVector{<:Branch}, θ::AbstractVector, targets::StateSpace) + return sum(branch -> measure(F, branch, θ, targets), branches) end ########################################################## bifurcation distance and velocity -distance(F::Function,z::BorderedArray,θ::AbstractVector) = [ F(z,θ); det(∂Fu(F,z,θ)) ] -function velocity( F::Function, z::BorderedArray, θ::AbstractVector; newtonOptions=NewtonPar(verbose=false,maxIter=800,tol=1e-6) ) - ∂implicit, _, _ = newtonOptions.linsolver( -jacobian( z->distance(F,z,θ), z )', [zero(z.u);one(z.p)] ) - return gradient( θ->distance(F,z,θ)'∂implicit, θ ) +distance(F::Function, z::BorderedArray, θ::AbstractVector) = [F(z, θ); det(∂Fu(F, z, θ))] +function velocity(F::Function, z::BorderedArray, θ::AbstractVector; newton_options=NewtonPar(verbose=false, max_iterations=800, tol=1e-6)) + ∂implicit, _, _ = newton_options.linsolver(-jacobian(z -> distance(F, z, θ), z)', [zero(z.u); one(z.p)]) + return gradient(θ -> distance(F, z, θ)'∂implicit, θ) end ################################################################################ -function tangent_field( F::Function, z::BorderedArray, θ::AbstractVector ) - field = kernel(∂Fz(F,z,θ);nullity=length(z.p)) - return norm(field)^(-1) * BorderedArray( field[Not(end)], field[end] ) # unit tangent field +function tangent_field(F::Function, z::BorderedArray, θ::AbstractVector) + field = kernel(∂Fz(F, z, θ); nullity=length(z.p)) + return norm(field)^(-1) * BorderedArray(field[Not(end)], field[end]) # unit tangent field end ################################################################################ using SpecialFunctions: erf -function window_function( parameter::AbstractVector, z::BorderedArray; β::Real=10 ) - pMin,pMax = extrema(parameter) - return ( 1 + erf((β*(z.p-pMin)-3)/√2) )/2 * ( 1 - ( 1 + erf((β*(z.p-pMax)+3)/√2) )/2 ) +function window_function(parameter::AbstractVector, z::BorderedArray; β::Real=10) + p_min, p_max = extrema(parameter) + return (1 + erf((β * (z.p - p_min) - 3) / √2)) / 2 * (1 - (1 + erf((β * (z.p - p_max) + 3) / √2)) / 2) end \ No newline at end of file diff --git a/src/Structures.jl b/src/Structures.jl index 851528e..275c5c6 100644 --- a/src/Structures.jl +++ b/src/Structures.jl @@ -1,4 +1,4 @@ -import Base: length,show,push! +import Base: length, show, push! struct StateSpace{N,T} roots::AbstractVector{<:AbstractVector{<:AbstractVector}} @@ -26,9 +26,9 @@ Define state space with targets to be used in optimisation - `parameter` one dimensional bifurcation parameter grid `p` - `targets` vector of target locations """ -function StateSpace(dimension::Integer,parameter::AbstractRange,targets::AbstractVector; nRoots::Integer=2, eltype::DataType=Float64) - roots, nTargets = fill( [zero(SizedVector{dimension,eltype})] ,nRoots ), length(targets) - return StateSpace{dimension,eltype}( SizedVector{nRoots}(roots), StepRangeLen{eltype}(parameter), SVector{nTargets,eltype}(targets) ) +function StateSpace(dimension::Integer, parameter::AbstractRange, targets::AbstractVector; nRoots::Integer=2, eltype::DataType=Float64) + roots, nTargets = fill([zero(SizedVector{dimension,eltype})], nRoots), length(targets) + return StateSpace{dimension,eltype}(SizedVector{nRoots}(roots), StepRangeLen{eltype}(parameter), SVector{nTargets,eltype}(targets)) end """ @@ -42,38 +42,38 @@ Initialises vector of named tuples that contain the following fields - `ds` arclength steps sizes between continuation points - `bif` boolean telling us if point is a bifurcation """ -BranchPoint{V,T} = NamedTuple{(:z,:λ,:ds,:bif),Tuple{BorderedArray{V,T},Vector{Complex{T}},T,Bool}} +BranchPoint{V,T} = NamedTuple{(:z, :λ, :ds, :bif),Tuple{BorderedArray{V,T},Vector{Complex{T}},T,Bool}} Branch{V,T} = Vector{BranchPoint{V,T}} -dim(branch::Branch) = length(branch) > 0 ? length(first(branch).z)-1 : 0 +dim(branch::Branch) = length(branch) > 0 ? length(first(branch).z) - 1 : 0 -function push!(branch::Branch,state::ContState) - z,∂z = copy(state.z), [state.τ.u;state.τ.p] - push!(branch, ( z=z, λ=state.eigvals, ds=norm(∂z)*abs(state.ds), bif=detectBifucation(state) )) +function push!(branch::Branch, state::ContState) + z, ∂z = copy(state.z), [state.τ.u; state.τ.p] + !any(isnan.(∂z)) && push!(branch, (z=z, λ=state.eigvals, ds=norm(∂z) * abs(state.ds), bif=detect_bifurcation(state))) end -import Base: zero,isapprox,unique,+ ################################## methods for BorderedArray +import Base: zero, isapprox, unique, + ################################## methods for BorderedArray import LinearAlgebra: norm -+(x::BorderedArray, y::BorderedArray) = BorderedArray(x.u+y.u,x.p+y.p) -norm(x::Branch, y::Branch) = mapreduce( (x,y)->norm(x.z.u-y.z.u)+norm(x.z.p-y.z.p),+,x,y)/min(length(x),length(y)) ++(x::BorderedArray, y::BorderedArray) = BorderedArray(x.u + y.u, x.p + y.p) +norm(x::Branch, y::Branch) = mapreduce((x, y) -> norm(x.z.u - y.z.u) + norm(x.z.p - y.z.p), +, x, y) / min(length(x), length(y)) -zero(::Type{BorderedArray{T,U}}; ϵ::Bool=true) where {T<:Union{Number,StaticArray},U<:Union{Number,StaticArray}} = BorderedArray(zero(T).+ϵ*eps(),zero(U).+ϵ*eps()) -isapprox( x::BorderedArray, y::BorderedArray; kwargs... ) = isapprox( [x.u;x.p], [y.u;y.p] ; kwargs...) -function unique(X::AbstractVector{T}; kwargs...) where T<:BorderedArray +zero(::Type{BorderedArray{T,U}}; ϵ::Bool=true) where {T<:Union{Number,StaticArray},U<:Union{Number,StaticArray}} = BorderedArray(zero(T) .+ ϵ * eps(), zero(U) .+ ϵ * eps()) +isapprox(x::BorderedArray, y::BorderedArray; kwargs...) = isapprox([x.u; x.p], [y.u; y.p]; kwargs...) +function unique(X::AbstractVector{T}; kwargs...) where {T<:BorderedArray} z = T[] - for x ∈ X - if all( zi -> ~isapprox( x, zi; kwargs... ), z ) - push!(z,x) + for x ∈ X + if all(zi -> ~isapprox(x, zi; kwargs...), z) + push!(z, x) end end return z end -function unique(X::AbstractVector{T}; kwargs...) where T<:Branch +function unique(X::AbstractVector{T}; kwargs...) where {T<:Branch} branches = T[] - for branch ∈ X - if all( branchᵢ -> ~isapprox( norm(branch,branchᵢ), 0; kwargs... ), branches ) - push!(branches,branch) + for branch ∈ X + if all(branchᵢ -> ~isapprox(norm(branch, branchᵢ), 0; kwargs...), branches) + push!(branches, branch) end end return branches @@ -84,13 +84,13 @@ function kernel(A::AbstractMatrix; nullity::Int=0) if iszero(nullity) # default to SVD return nullspace(A) - else - return qr(A').Q[:,end+1-nullity:end] + else + return qr(A').Q[:, end+1-nullity:end] end end #################################################### display methods show(io::IO, branches::Vector{Branch{V,T}}) where {V,T} = print(io, "Vector{Branch}[dim=$(dim(first(branches))) bifurcations=$(length(unique([ s.z for branch ∈ branches for s ∈ branch if s.bif ]; atol=0.1))) branches=$(length(branches)), states=$(sum(branch->length(branch),branches))]") -show(io::IO, M::MIME"text/plain", branches::Vector{Branch{V,T}}) where {V,T} = show(io,branches) -show(io::IO, states::StateSpace{N,T}) where {N,T} = print(io,"StateSpace{$N,$T}(parameters=$(states.parameter),targets=$(states.targets))") \ No newline at end of file +show(io::IO, M::MIME"text/plain", branches::Vector{Branch{V,T}}) where {V,T} = show(io, branches) +show(io::IO, states::StateSpace{N,T}) where {N,T} = print(io, "StateSpace{$N,$T}(parameters=$(states.parameter),targets=$(states.targets))") \ No newline at end of file diff --git a/src/Utils.jl b/src/Utils.jl index c4d843a..13bbdff 100644 --- a/src/Utils.jl +++ b/src/Utils.jl @@ -1,100 +1,103 @@ ############################################################################ hyperparameter updates -function getParameters(data::StateSpace{N,T}; maxIter::Int=800, tol=1e-12, kwargs...) where {N,T<:Number} - newtonOptions = NewtonPar(verbose=false,maxIter=maxIter,tol=T(tol)) - - # support for StaticArrays github.com/JuliaArrays/StaticArrays.jl/issues/73 - newtonOptions = @set newtonOptions.linsolver.useFactorization = false - return ContinuationPar( - - pMin=minimum(data.parameter), pMax=maximum(data.parameter), maxSteps=10*length(data.parameter), - ds=step(data.parameter), dsmax=step(data.parameter), dsmin=step(data.parameter), - - newtonOptions=newtonOptions, detectFold=false, detectBifurcation=true, - saveEigenvectors=false, nev=N ) +function getParameters(data::StateSpace{N,T}; max_iterations::Int=800, tol=1e-12, kwargs...) where {N,T<:Number} + newton_options = NewtonPar(verbose=false, max_iterations=max_iterations, tol=T(tol)) + + # support for StaticArrays github.com/JuliaArrays/StaticArrays.jl/issues/73 + newton_options = @set newton_options.linsolver.useFactorization = false + return ContinuationPar(p_min=minimum(data.parameter), p_max=maximum(data.parameter), max_steps=10 * length(data.parameter), + ds=step(data.parameter), dsmax=step(data.parameter), dsmin=step(data.parameter), newton_options=newton_options, detect_fold=false, detect_bifurcation=true, + save_eigenvectors=false, nev=N) end -function updateParameters!(parameters::ContinuationPar{T, S, E}, steady_states::Vector{Branch{V,T}}; - resolution=400 ) where {T<:Number, V<:AbstractVector{T}, S<:AbstractLinearSolver, E<:AbstractEigenSolver} +function updateParameters!(parameters::ContinuationPar{T,S,E}, steady_states::Vector{Branch{V,T}}; + resolution=400) where {T<:Number,V<:AbstractVector{T},S<:AbstractLinearSolver,E<:AbstractEigenSolver} # estimate scale from steady state curves - branch_points = map(length,steady_states) - ds = maximum(branch_points)*parameters.ds/resolution - parameters = setproperties(parameters;ds=ds,dsmin=ds,dsmax=ds) + branch_points = map(length, steady_states) + ds = maximum(branch_points) * parameters.ds / resolution + parameters = setproperties(parameters; ds=ds, dsmin=ds, dsmax=ds) end -function mean(f::Function,x::AbstractArray{T}; type::Symbol=:geometric, kwargs...) where T - return type == :geometric ? mapreduce(y->f(y)^(1/length(x)),*,x;init=1,kwargs...) : mapreduce(y->f(y)/length(x),+,x;init=zero(f(zero(T))),kwargs...) +function mean(f::Function, x::AbstractArray{T}; type::Symbol=:geometric, kwargs...) where {T} + return type == :geometric ? mapreduce(y -> f(y)^(1 / length(x)), *, x; init=1, kwargs...) : mapreduce(y -> f(y) / length(x), +, x; init=zero(f(zero(T))), kwargs...) end ############################################################################# training loop -function train!( F::Function, parameters::NamedTuple, data::StateSpace; - iter::Int=200, optimiser=Momentum(0.001), plot_solution = false, - convergence = 0.02, kwargs...) - - hyperparameters = getParameters(data;kwargs...) - Loss = steady_states = nothing - - ∇Loss = similar(parameters.θ) - trajectory = typeof(parameters.θ)[] - - for i=1:iter - try - steady_states = deflationContinuation(F,data.roots,parameters,hyperparameters;kwargs...) - Loss,∇Loss = ∇loss(F,steady_states,parameters.θ,data;kwargs...) - - catch error - printstyled(color=:red, "Iteration $i\tError = $error\n") end - printstyled(color=:yellow,"Iteration $i\tLoss = $Loss\n") - - printstyled(color=:blue,"$steady_states\n") - println("Parameters\t$(parameters.θ)") - println("Gradients\t$(∇Loss)") - - if isinf(Loss) throw("infinite loss") end - push!(trajectory,copy(parameters.θ)) - - predictions = unique([ s.z for branch ∈ steady_states for s ∈ branch if s.bif ]; atol=2*step(data.parameter)) - if (Loss < convergence) & (length(predictions) == length(data.targets)) - break - end - - if plot_solution>0 if i%plot_solution==0 plot(steady_states,data) end end - update!(optimiser, parameters.θ, ∇Loss ) - end - - return trajectory +function train!(F::Function, parameters::NamedTuple, data::StateSpace; + iter::Int=200, optimiser=Momentum(0.001), plot_solution=false, + convergence=0.02, kwargs...) + + hyperparameters = getParameters(data; kwargs...) + Loss = steady_states = nothing + + ∇Loss = similar(parameters.θ) + trajectory = typeof(parameters.θ)[] + + for i = 1:iter + try + steady_states = deflationContinuation(F, data.roots, parameters, hyperparameters; kwargs...) + Loss, ∇Loss = ∇loss(F, steady_states, parameters.θ, data; kwargs...) + + catch error + printstyled(color=:red, "Iteration $i\tError = $error\n") + end + printstyled(color=:yellow, "Iteration $i\tLoss = $Loss\n") + + printstyled(color=:blue, "$steady_states\n") + println("Parameters\t$(parameters.θ)") + println("Gradients\t$(∇Loss)") + + if isinf(Loss) + throw("infinite loss") + end + push!(trajectory, copy(parameters.θ)) + + predictions = unique([s.z for branch ∈ steady_states for s ∈ branch if s.bif]; atol=2 * step(data.parameter)) + if (Loss < convergence) & (length(predictions) == length(data.targets)) + break + end + + if plot_solution > 0 + if i % plot_solution == 0 + plot(steady_states, data) + end + end + update!(optimiser, parameters.θ, ∇Loss) + end + + return trajectory end ############################################################################## loss evaluation helper function loss(F::Function, θ::AbstractVector, data::StateSpace, hyperparameters::ContinuationPar; kwargs...) - try - parameters = (θ=θ,p=minimum(data.parameter)) - steady_states = deflationContinuation(F,data.roots,parameters,hyperparameters;kwargs...) - return loss(F,steady_states,θ,data;kwargs...) + try + parameters = (θ=θ, p=minimum(data.parameter)) + steady_states = deflationContinuation(F, data.roots, parameters, hyperparameters; kwargs...) + return loss(F, steady_states, θ, data; kwargs...) - catch error - printstyled(color=:red,"$error\n") - return NaN - end + catch error + printstyled(color=:red, "$error\n") + return NaN + end end function loss(F::Function, θ::AbstractVector, data::StateSpace; kwargs...) - return loss(F,θ,data,getParameters(data); kwargs...) + return loss(F, θ, data, getParameters(data); kwargs...) end function ∇loss(F::Function, θ::AbstractVector, data::StateSpace, hyperparameters::ContinuationPar; kwargs...) - try - parameters = (θ=θ,p=minimum(data.parameter)) - steady_states = deflationContinuation(F,data.roots,parameters,hyperparameters;kwargs...) - return ∇loss(F,steady_states,θ,data;kwargs...) - - catch error - printstyled(color=:red,"$error\n") - return NaN,fill(NaN,length(θ)) - end + try + parameters = (θ=θ, p=minimum(data.parameter)) + steady_states = deflationContinuation(F, data.roots, parameters, hyperparameters; kwargs...) + return ∇loss(F, steady_states, θ, data; kwargs...) + + catch error + printstyled(color=:red, "$error\n") + return NaN, fill(NaN, length(θ)) + end end function ∇loss(F::Function, θ::AbstractVector, data::StateSpace; kwargs...) - return ∇loss(F,θ,data,getParameters(data); kwargs...) + return ∇loss(F, θ, data, getParameters(data); kwargs...) end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 0584a0f..30a1cc7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,66 +5,66 @@ using Flux: Optimise error_tolerance = 0.004 @testset "Gradients" begin - @testset "Saddle Node" begin - include("minimal/saddle-node.jl") - - @time for θ ∈ [[5.0,1.0],[5.0,-1.0],[-1.0,-1.0],[2.5,-1.0],[-1.0,2.5],[-4.0,1.0]] - x, y = finite_difference_gradient(θ->loss(F,θ,X),θ), ∇loss(F,θ,X)[2] - @test acos(x'y/(norm(x)*norm(y)))/π < error_tolerance - end - end - - @testset "Pitchfork" begin - include("minimal/pitchfork.jl") - - @time for θ ∈ [[1.0,1.0],[3.0,3.0],[1.0,-1.0],[3.0,-3.0],[-1.0,-1.0],[-3.0,-3.0],[-1.0,1.0],[-3.0,3.0]] - x, y = finite_difference_gradient(θ->loss(F,θ,X),θ), ∇loss(F,θ,X)[2] - @test acos(x'y/(norm(x)*norm(y)))/π < error_tolerance - end - end - - @testset "Two State" begin - include("applied/two-state.jl") - - @time for θ ∈ [[0.5, 0.5, 0.5470, 2.0, 7.5],[0.5, 0.5, 1.5470, 2.0, 7.5],[0.5, 0.88, 1.5470, 2.0, 7.5]] - x, y = finite_difference_gradient(θ->loss(F,θ,X),θ), ∇loss(F,θ,X)[2] - @test acos(x'y/(norm(x)*norm(y)))/π < error_tolerance - end - end + @testset "Saddle Node" begin + include("minimal/saddle-node.jl") + + @time for θ ∈ [[5.0, 1.0], [5.0, -1.0], [-1.0, -1.0], [2.5, -1.0], [-1.0, 2.5], [-4.0, 1.0]] + x, y = finite_difference_gradient(θ -> loss(F, θ, X), θ), ∇loss(F, θ, X)[2] + @test acos(x'y / (norm(x) * norm(y))) / π < error_tolerance + end + end + + @testset "Pitchfork" begin + include("minimal/pitchfork.jl") + + @time for θ ∈ [[1.0, 1.0], [3.0, 3.0], [1.0, -1.0], [3.0, -3.0], [-1.0, -1.0], [-3.0, -3.0], [-1.0, 1.0], [-3.0, 3.0]] + x, y = finite_difference_gradient(θ -> loss(F, θ, X), θ), ∇loss(F, θ, X)[2] + @test acos(x'y / (norm(x) * norm(y))) / π < error_tolerance + end + end + + @testset "Two State" begin + include("applied/two-state.jl") + + @time for θ ∈ [[0.5, 0.5, 0.5470, 2.0, 7.5], [0.5, 0.5, 1.5470, 2.0, 7.5], [0.5, 0.88, 1.5470, 2.0, 7.5]] + x, y = finite_difference_gradient(θ -> loss(F, θ, X), θ), ∇loss(F, θ, X)[2] + @test acos(x'y / (norm(x) * norm(y))) / π < error_tolerance + end + end end error_tolerance = 0.1 @testset "Optimisations" begin -include("applied/two-state.jl") + include("applied/two-state.jl") - @testset "Mutual Inhibition" begin - parameters = ( θ=SizedVector{5}(0.5,0.5,0.5470,1.0,1.5), p=minimum(X.parameter) ) - trajectory = train!( F, parameters, X; iter=300, optimiser=Optimise.ADAM(0.01) ) + @testset "Mutual Inhibition" begin + parameters = (θ=SizedVector{5}(0.5, 0.5, 0.5470, 1.0, 1.5), p=minimum(X.parameter)) + trajectory = train!(F, parameters, X; iter=300, optimiser=Optimise.ADAM(0.01)) - steady_states = deflationContinuation(F,X.roots,(p=minimum(X.parameter),θ=trajectory[end]),getParameters(X)) - bifurcations = [ s.z for branch ∈ steady_states for s ∈ branch if s.bif ] + steady_states = deflationContinuation(F, X.roots, (p=minimum(X.parameter), θ=trajectory[end]), getParameters(X)) + bifurcations = [s.z for branch ∈ steady_states for s ∈ branch if s.bif] - @test all( z->any(isapprox.(X.targets,z.p; atol=error_tolerance)), bifurcations) - end + @test all(z -> any(isapprox.(X.targets, z.p; atol=error_tolerance)), bifurcations) + end - @testset "Mutual Activation" begin - parameters = ( θ=SizedVector{5}(0.246,-1.113,-0.87,-1.36,1.28), p=minimum(X.parameter) ) - trajectory = train!( F, parameters, X; iter=300, optimiser=Optimise.ADAM(0.01) ) + @testset "Mutual Activation" begin + parameters = (θ=SizedVector{5}(0.246, -1.113, -0.87, -1.36, 1.28), p=minimum(X.parameter)) + trajectory = train!(F, parameters, X; iter=300, optimiser=Optimise.ADAM(0.01)) - steady_states = deflationContinuation(F,X.roots,(p=minimum(X.parameter),θ=trajectory[end]),getParameters(X)) - bifurcations = [ s.z for branch ∈ steady_states for s ∈ branch if s.bif ] + steady_states = deflationContinuation(F, X.roots, (p=minimum(X.parameter), θ=trajectory[end]), getParameters(X)) + bifurcations = [s.z for branch ∈ steady_states for s ∈ branch if s.bif] - @test all( z->any(isapprox.(X.targets,z.p; atol=error_tolerance)), bifurcations) - end + @test all(z -> any(isapprox.(X.targets, z.p; atol=error_tolerance)), bifurcations) + end - include("applied/double-exclusive.jl") - @testset "Double Exclusive" begin - parameters = ( θ=SizedVector{21}(-0.5498107083383627, -0.30983502459648393, 0.631939862866151, -1.3407390261868093, 0.6331343612914606, -0.4049028383285403, -0.8524704190293824, -0.11432918583984555, -1.7299987854365662, -2.079032179004907, 0.4916870522978877, -0.010788745553567142, -0.3878513750140796, -0.08334187371052072, -0.2772816729274815, 0.8285363415116678, 0.054887187754647675, -0.49494957503215126, -0.8693294194842291, 0.9289938806228071, 0.22800165770588346), p=minimum(X.parameter) ) - trajectory = train!( F, parameters, X; iter=300, optimiser=Optimise.ADAM(0.01) ) + include("applied/double-exclusive.jl") + @testset "Double Exclusive" begin + parameters = (θ=SizedVector{21}(-0.5498107083383627, -0.30983502459648393, 0.631939862866151, -1.3407390261868093, 0.6331343612914606, -0.4049028383285403, -0.8524704190293824, -0.11432918583984555, -1.7299987854365662, -2.079032179004907, 0.4916870522978877, -0.010788745553567142, -0.3878513750140796, -0.08334187371052072, -0.2772816729274815, 0.8285363415116678, 0.054887187754647675, -0.49494957503215126, -0.8693294194842291, 0.9289938806228071, 0.22800165770588346), p=minimum(X.parameter)) + trajectory = train!(F, parameters, X; iter=300, optimiser=Optimise.ADAM(0.01)) - steady_states = deflationContinuation(F,X.roots,(p=minimum(X.parameter),θ=trajectory[end]),getParameters(X)) - bifurcations = [ s.z for branch ∈ steady_states for s ∈ branch if s.bif ] + steady_states = deflationContinuation(F, X.roots, (p=minimum(X.parameter), θ=trajectory[end]), getParameters(X)) + bifurcations = [s.z for branch ∈ steady_states for s ∈ branch if s.bif] - @test all( z->any(isapprox.(X.targets,z.p; atol=error_tolerance)), bifurcations) - end + @test all(z -> any(isapprox.(X.targets, z.p; atol=error_tolerance)), bifurcations) + end end \ No newline at end of file