Add evolutionary algorithm

Project.toml

@@ -4,9 +4,11 @@ authors = ["QC-Devs"]
 version = "0.1.0"
 
 [deps]
+Evolutionary = "86b6b26d-c046-49b6-aa0b-5f0f74682bd6"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 

docs/Project.toml

@@ -1,3 +1,5 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Evolutionary = "86b6b26d-c046-49b6-aa0b-5f0f74682bd6"
 NICE = "170674ed-ea2e-4c97-a263-9f057371600a"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"

docs/make.jl

@@ -9,6 +9,7 @@ makedocs(;
     repo="https://github.com/theochem/NICE.jl/blob/{commit}{path}#{line}",
     sitename="NICE.jl",
     format=Documenter.HTML(;
+        repolink="https://github.com/theochem/NICE.jl/",
         prettyurls=get(ENV, "CI", "false") == "true",
         canonical="https://nice.qcdevs.org/",
         edit_link="main",

docs/src/index.md

@@ -4,7 +4,7 @@ CurrentModule = NICE
 
 # NICE.jl
 
-Documentation for [NICE.jl](nice.qcdevs.org).
+Documentation for NICE.jl.
 
 ```@contents
 Pages = ["api.md"]

src/Exact.jl

@@ -1,11 +1,25 @@
 using LinearAlgebra
 
+using Reexport
+
 using SciMLBase
 using ForwardDiff
 using NonlinearSolve
 
+using Evolutionary
+@reexport using Evolutionary: ES, CMAES, GA, DE, NSGA2, TreeGP
+
+export solve!
+export evolutionary_solve!
+
+struct ConvergenceError <: Exception
+    msg::String
+end
+
+Base.showerror(io::IO, e::ConvergenceError) = print(io, e.msg)
+
 """
-    solve(rxn_system, K_eqs; maxiters=1000, abstol=1.0e-9, reltol=0.0)
+    solve!(rxn_system, K_eqs; maxiters=1000, abstol=1.0e-9, reltol=0.0)
 
 Solve the system deterministically.
 
@@ -21,58 +35,42 @@ K_\\text{eq}^{\\left(m\\right)} =
     \\prod_{n=1}^N {\\left( a_n^{\\left(0\\right)} +
     \\sum_{m=1}^M c_{mn} \\xi_m \\right)}^{c_{mn}}
 ```
-using user-provided ``K_\text{eq}`` values and the Newton's method + Trust
+using user-provided ``K_\\text{eq}`` values and the Newton's method + Trust
 Region method.
 """
-function solve(
+function solve!(
     rxn_system::ReactionSystem,
-    K_eqs::AbstractVector{Float64};
-    maxiters::Integer=1000,
-    abstol::Real=1.0e-9,
-    reltol::Real=0.0,
+    K_eqs::Vector{Float64};
+    maxiter::Int=1000,
+    abstol::Float64=1.0e-9,
+    reltol::Float64=0.0,
 )
-    # Define the objective function
-    function f_K_eqs(f, ξ, _)
-        for i in 1 : rxn_system.n_reaction
-            f_i1 = K_eqs[i]
-            f_i2 = 1.0
-            for j in 1 : rxn_system.n_species
-                # Reactants
-                if rxn_system.stoich[j, i] < 0
-                    f_i1 *= (
-                        rxn_system.concs_init[j] + rxn_system.stoich[j, :] · ξ
-                    ) ^ -rxn_system.stoich[j, i]
-                # Products
-                elseif rxn_system.stoich[j, i] > 0
-                    f_i2 *= (
-                        rxn_system.concs_init[j] + rxn_system.stoich[j, :] · ξ
-                    ) ^ rxn_system.stoich[j, i]
-                end
-            end
-            f[i] = f_i1 - f_i2
-        end
-        nothing
-    end
-    # Compute reaction extents
-    ξ = rxn_system.stoich \ (rxn_system.concs - rxn_system.concs_init)
-    # Create objective function and Jacobian (automatic, with autodiff)
+    # Define nonlinear problem (derives Jacobian via autodiff)
     problem = NonlinearSolve.NonlinearProblem(
-        f_K_eqs, ξ;
-        maxiters=maxiters, abstol=abstol, reltol=reltol,
+        # Objective function
+        (f, ξ, _) -> f_K_eqs(f, ξ, rxn_system, K_eqs),
+        # Initial guess for reaction extents ξ
+        rxn_system.stoich \ (rxn_system.concs - rxn_system.concs_init);
+        # Solver options
+        maxiters=maxiter,
+        abstol=abstol,
+        reltol=reltol,
     )
     # Run the nonlinear optimization
     solution = NonlinearSolve.solve(problem, NonlinearSolve.TrustRegion())
-    # If nonlinear optimization was successful, update `rxn_system.concs`
-    if SciMLBase.successful_retcode(solution)
-        rxn_system.concs .= rxn_system.concs_init
-        rxn_system.concs .+= rxn_system.stoich * solution.u
+    # If nonlinear optimization was successful, update `rxn_system.concs`;
+    # otherwise throw ConvergenceError exception
+    if !SciMLBase.successful_retcode(solution)
+        throw(ConvergenceError("Did not converge"))
     end
+    rxn_system.concs .= rxn_system.concs_init
+    rxn_system.concs .+= rxn_system.stoich * solution.u
     # Return the nonlinear solution
     solution
 end
 
 """
-    solve(rxn_system, K_eqs; maxiters=1000, abstol=1.0e-9, reltol=0.0)
+    solve!(rxn_system, K_eqs; maxiters=1000, abstol=1.0e-9, reltol=0.0)
 
 Solve the system deterministically.
 
@@ -94,13 +92,100 @@ Instead of specifying the ``K_\\text{eq}`` values directly, they are
 approximated here using the forward or reverse rate constants and a parameter
 ``\\phi``.
 """
-function solve(
+function solve!(
     rxn_system::ReactionSystem;
-    φ::Real=1.0,
+    φ::Float64=1.0,
     rate_consts::Symbol=:forward,
-    maxiters::Integer=1000,
-    abstol::Real=1.0e-9,
-    reltol::Real=0.0,
+    maxiter::Int=1000,
+    abstol::Float64=1.0e-9,
+    reltol::Float64=0.0,
+)
+    solve!(
+        rxn_system,
+        approximate_K_eqs(rxn_system, φ, rate_consts);
+        maxiter=maxiter,
+        abstol=abstol,
+        reltol=reltol,
+    )
+end
+
+"""
+    evolutionary_solve!(rxn_system, K_eqs, algorithm; abstol=Inf, reltol=Inf, maxiter=1000)
+
+Solve the system using an evolutionary algorithm.
+
+Solves the nonlinear problem described by [`solve!`](@ref) by using one of the
+evolutionary algorithms from
+[Evolutionary.jl](https://wildart.github.io/Evolutionary.jl/stable/).
+These algorithm types are reexported by NICE.jl.
+"""
+function evolutionary_solve!(
+    rxn_system::ReactionSystem,
+    K_eqs::Vector{Float64},
+    algorithm::TOpt;
+    abstol::Float64=Inf,
+    reltol::Float64=Inf,
+    maxiter::Int=1000,
+) where {TOpt <: Evolutionary.AbstractOptimizer}
+    f = similar(K_eqs)
+    solution = Evolutionary.optimize(
+        # Objective function
+        ξ -> f_K_eqs_cma(f, ξ, rxn_system, K_eqs),
+        # Initial guess for reaction extents ξ
+        rxn_system.stoich \ (rxn_system.concs - rxn_system.concs_init),
+        # Evolutionary algorithm
+        algorithm,
+        # Solver options
+        Evolutionary.Options(iterations=maxiter, abstol=abstol, reltol=reltol),
+    )
+    # If nonlinear optimization was successful, update `rxn_system.concs`
+    # otherwise throw ConvergenceError exception
+    if !solution.converged
+        throw(ConvergenceError("Did not converge"))
+    end
+    rxn_system.concs .= rxn_system.concs_init
+    rxn_system.concs .+= rxn_system.stoich * solution.minimizer
+    # Return the nonlinear solution
+    solution
+end
+
+"""
+    evolutionary_solve!(rxn_system, K_eqs, algorithm; abstol=Inf, reltol=Inf, maxiter=1000)
+
+Solve the system using an evolutionary algorithm.
+
+Solves the nonlinear problem described by [`solve!`](@ref) by using one of the
+evolutionary algorithms from
+[Evolutionary.jl](https://wildart.github.io/Evolutionary.jl/stable/).
+These algorithm types are reexported by NICE.jl.
+
+Instead of specifying the ``K_\\text{eq}`` values directly, they are
+approximated here using the forward or reverse rate constants and a parameter
+``\\phi``.
+"""
+function evolutionary_solve!(
+    rxn_system::ReactionSystem,
+    algorithm::TOpt;
+    φ::Float64=1.0,
+    rate_consts::Symbol=:forward,
+    abstol::Float64=Inf,
+    reltol::Float64=Inf,
+    maxiter::Int=1000,
+) where {TOpt <: Evolutionary.AbstractOptimizer}
+    evolutionary_solve!(
+        rxn_system,
+        approximate_K_eqs(rxn_system, φ, rate_consts),
+        algorithm;
+        abstol=abstol,
+        reltol=reltol,
+        maxiter=maxiter,
+    )
+end
+
+function approximate_K_eqs(
+    rxn_system::ReactionSystem,
+    φ::Float64,
+    rate_consts::Symbol,
 )
     if rate_consts === :forward
         # Evaluate K_eqs from forward rate constants and φ_fwd
@@ -115,5 +200,32 @@ function solve(
         throw(ArgumentError("invalid rate_consts value $rate_consts \
                              (must be :forward or :reverse)"))
     end
-    solve(rxn_system, K_eqs; maxiters=maxiters, abstol=abstol, reltol=reltol)
+    K_eqs
+end
+
+function f_K_eqs(f, ξ, rxn_system, K_eqs)::Nothing
+    for i in 1 : rxn_system.n_reaction
+        f_i1 = K_eqs[i]
+        f_i2 = 1.0
+        for j in 1 : rxn_system.n_species
+            # Reactants
+            if rxn_system.stoich[j, i] < 0
+                f_i1 *= (
+                    rxn_system.concs_init[j] + rxn_system.stoich[j, :] · ξ
+                ) ^ -rxn_system.stoich[j, i]
+            # Products
+            elseif rxn_system.stoich[j, i] > 0
+                f_i2 *= (
+                    rxn_system.concs_init[j] + rxn_system.stoich[j, :] · ξ
+                ) ^ rxn_system.stoich[j, i]
+            end
+        end
+        f[i] = f_i1 - f_i2
+    end
+    nothing
+end
+
+function f_K_eqs_cma(f::Vector{Float64}, ξ::Vector{Float64}, rxn_system::ReactionSystem, K_eqs::Vector{Float64})::Float64
+    f_K_eqs(f, ξ, rxn_system, K_eqs)
+    √(f · f)
 end

src/NEKMC.jl

@@ -1,24 +1,26 @@
 using LinearAlgebra
 
+export simulate!
+
 """
-    simulate(rxn_system;
+    simulate!(rxn_system;
              n_iter=Int(1e+8), chunk_iter=Int(1e+4),
              ε=1.0e-4, ε_scale=1.0, ε_concs=0.0,
              tol_t=Inf, tol_ε=0.0, tol_concs=0.0)
 
 Run a *N*et-*E*vent *K*inetic *M*onte *C*arlo (NEKMC) simulation to find the
 equilibrium concentrations of the reaction system.
 """
-function simulate(
+function simulate!(
     rxn_system::ReactionSystem;
-    n_iter::Integer=Int(1e+8),
-    chunk_iter::Integer=Int(1e+4),
-    ε::Real=1.0e-4,
-    ε_scale::Real=1.0,
-    ε_concs::Real=0.0,
-    tol_t::Real=Inf,
-    tol_ε::Real=0.0,
-    tol_concs::Real=0.0,
+    n_iter::Int=Int(1e+8),
+    chunk_iter::Int=Int(1e+4),
+    ε::Float64=1.0e-4,
+    ε_scale::Float64=1.0,
+    ε_concs::Float64=0.0,
+    tol_t::Float64=Inf,
+    tol_ε::Float64=0.0,
+    tol_concs::Float64=0.0,
 )
     # Allocate probability vector, Δconcs vector, Δt
     pvec = zeros(Float64, rxn_system.n_reaction)
@@ -82,7 +84,7 @@ end
 
 function select_reaction(
     rxn_system::ReactionSystem,
-    pvec::AbstractVector{Float64},
+    pvec::Vector{Float64},
 )
     # Update probability vector
     p = 0.
@@ -103,8 +105,8 @@ end
 
 function do_reaction(
     rxn_system::ReactionSystem,
-    i_rxn::Integer,
-    ε::Real,
+    i_rxn::Int,
+    ε::Float64,
 )
     rate = rxn_system.net_rates[i_rxn]
     if rate >= 0

src/NICE.jl

@@ -1,16 +1,13 @@
 """
 NICE (N-species Ice Table) module -- provides simultaneous equilibria solvers:
-- [`simulate`](@ref) Net-event kinetic Monte Carlo solver
-- [`solve`](@ref) Nonlinear equations solver (Newton trust region method)
+- [`simulate!`](@ref) Net-event kinetic Monte Carlo solver
+- [`solve!`](@ref) Exact nonlinear equations solver (Newton trust region method)
+- [`evolutionary_solve!`](@ref) Evolutionary nonlinear equations solver (Newton trust region method)
 
 The [`ReactionSystem`](@ref) type is provided as an interface for these solvers.
 """
 module NICE
 
-export ReactionSystem
-export simulate
-export solve
-
 include("ReactionSystem.jl")
 include("NEKMC.jl")
 include("Exact.jl")