For now, the package is not an official Julia package. Therefore, the user should manually add it by running:
import Pkg
Pkg.add("https://github.com/d4v1dcub/ApproxMasterEqs.jl.git")
This is a package to perform the numerical integration of some Approximated Master Equations for the dynamics of binary variables in graphs. It contains two different methods: the Cavity Master Equation (CME) [1] [2] [3] and the Conditional Dynamic Approximation (CDA) [4]. Provided a model and some network representing the interactions, the package estimates the local probabilities of observing a specific configuration of the variables at time
The two main functions implemented in the package so far are CME_KSAT
and CDA_KSAT
, with the following arguments
(ratefunc::Function, rargs_cst, rarg_build::Function;
graph::HGraph=build_empty_graph(), N::Int64=0, K::Int64=0,
alpha::Union{Float64, Int64}=0.0, seed_g::Int64=rand(1:typemax(Int64)),
links::Matrix{Int8}=Matrix{Int8}(undef, 0, 0), seed_l::Int64=rand(1:typemax(Int64)),
tspan::Vector{Float64}=[0.0, 1.0], p0::Float64=0.5, method=VCABM(), eth::Float64=1e6,
cbs_save::CallbackSet=CallbackSet(), dt_s::Float64=0.1, abstol::Float64=1e6, reltol::Float64=1e3)
In its first version, the package works for models on hypergraphs where only one configuration in the factor nodes unsatisfies the interaction (KSAT like). This contains as a particular case a model with pairwise interactions (2SAT like)
The package has its own implementation of hypergraph structures.
mutable struct HGraph
# This structure holds the hypergraph information
N::Int64 # number of variable nodes
M::Int64 # number of hyperedges
K::Int64 # number of nodes in each hyperedge
var_2_he::Array{Array{Int64, 1}, 1} # list of hyperedges for each variable
he_2_var::Matrix{Int64} # list of variables per hyperedge
degrees::Vector{Int64} # list of variable nodes' degrees
nodes_in::Array{Dict{Int64, Int64}, 1} # dictionary with key=node_index and val=place_in_he
#.......
end
that can be initialized with
build_HGraph(var_2_he::Array{Array{Int64, 1}, 1} , he_2_var::Matrix{Int64}, degrees::Vector{Int64})
There are a couple of implemented examples:
build_RR_HGraph(N::Int64, c::Int64, K::Int64, idum::Int64=rand(1:typemax(Int64))) # Random Regular Graphs
build_ER_HGraph(N::Int64, c::Union{Float64, Int64}, K::Int64, idum::Int64=rand(1:typemax(Int64))) # ErdösRényi
If the user does not provide a graph to the functions CME_KSAT
or CDA_KSAT
, an empty graph is taken as default. The functions then expect to receive three numbers:
The couplings or links for the KSAT boolean formula can be input via the parameter links::Matrix{Float64}
. The indexes are links[he, i]
with
For a neat example see the end of the Section 5
The information about the interactions goes into the transition rates. These are related to the first three arguments of the functions CME_KSAT
and CDA_KSAT
:
ratefunc::Function # function that computes the transition rate
rargs_cst # constant arguments that 'ratefunc' will receive when evaluated
rarg_build::Function # function that computes the nonconstant arguments to be obtained
# at each integration step and then passed to 'ratefunc'
The structure of ratefunc
must fit the following
ratefunc(Ep::Int64, Em::Int64, rateargs...)
where Ep
is the local energy when the variable is positive Em
is the local energy when the variable is negative
To build the rest of the arguments the user should use a function like
rarg_build(graph::HGraph, st::State_CME, rargs_cst...)
rarg_build(graph::HGraph, st::State_CDA, rargs_cst...)
The first argument is the graph, with all the associated information (see previous section).
The second argument is a struct
with the following information:
mutable struct State_CME
p_cav::Array{Float64, 4} # cavity conditional probabilities p_cav[hyperedge, site_cond, s_cond, chain]
probi::Vector{Float64} # singlesite probabilities probi[node]
pu_cond::Array{Float64, 3} # cavity conditional probabilities of partially unsatisfied
# hyperedges pu_cond[hyperedge, site_cond, s_cond]
p_joint_u::Vector{Float64} # probability of having the hyperedge unsatisfied p_joint_u[hyperedge]
end
mutable struct State_CDA
p_joint::Matrix{Float64} # joint probabilities p_joint[hyperedge, chain]
pu_cond::Array{Float64, 3} # conditional probabilities of partially unsatisfied
# hyperedges pu_cond[hyperedge, site_cond, s_cond]
p_joint_u::Vector{Float64} # probability of having the hyperedge unsatisfied p_joint_u[hyperedge]
end
All the other constant arguments of the transition rates must go into rargs_cst...
The function rarg_build()
must return a list of arguments ready to be inserted into the function ratefunc
As an example, let us see the implementation of the Focused Metropolis Search algorithm for KSAT optimization
# This function computes the energy in a KSAT formula, where there is only one unsatisfied
# configuration of the variables in a clause
function ener(p_joint_u::Vector{Float64})
e = 0.0
for he in eachindex(p_joint_u)
e += p_joint_u[he]
end
return e
end
# This is the rate of FMS algorithm for KSAT
function rate_FMS_KSAT(Ep::Int64, Em::Int64, eta::Float64, avE::Float64, K::Number, N::Int64)
dE = Em  Ep
if dE > 0
return Ep * N / K / avE * eta ^ dE
else
return Ep * N / K / avE
end
end
# Each rate function needs some specific arguments. The general form of these functions
function build_args_rate_FMS(graph::HGraph, st::State_CME, eta::Float64)
avE = ener(st.p_joint_u)
return eta, avE, graph.K, graph.N
end
# Each rate function needs some specific arguments. The general form of these functions
function build_args_rate_FMS(graph::HGraph, st::State_CDA, eta::Float64)
avE = ener(st.p_joint_u)
return eta, avE, graph.K, graph.N
end
A simple example is written in the file "package_dir/test/runtests.jl"
using ApproxMasEq
using Test
N = 10
alpha = 2
K = 3
seed = 1
eta = 1.0
rargs = [eta]
answ_CME = CME_KSAT(rate_FMS_KSAT, rargs, build_args_rate_FMS, N=N, K=K, alpha=alpha, seed_g=seed)
answ_CDA = CDA_KSAT(rate_FMS_KSAT, rargs, build_args_rate_FMS, N=N, K=K, alpha=alpha, seed_g=seed)
Here, a hypergraph with
The model is the one in the example of Section 4: Focused Metropolis Search algorithm in KSAT. The constant argument for the transition rates is eta=1 (
Other keyword arguments that can be specified before the numerical integration:

p0::Float=0.5
is the probability for generating the initial conditions. Every variable is set to 1 or 1 with$p(1) = 1p(1) = p_0$ . 
tspan::Vector{Float64}=[0.0, 1.0]
is the time interval for the integration. 
method=VCABM()
is the numerical method for the integration performed with the package OrdinaryDiffEqs. See the full list here. 
abstol::Float64=1e6
absolute tolerance for the numerical integration with OrdinaryDiffEqs. 
reltol::Float64=1e6
relative tolerance for the numerical integration with OrdinaryDiffEqs.
At this point, the user should be able to run a simple example and collect the output of the functions CME_KSAT
or CDA_KSAT
. These functions return an object of type SciMLBase.ODESolution
(see here) with the solution saved in tspan
sampled at intervals of length dts
. The latter is a parameter of the functions CME_KSAT
and CDA_KSAT
(dt_s::Float64=0.1
).
To get other quantities the user should use callbacks. This is allowed by the package DiffEqCallbacks. A stopping condition is implemented by default that stops the integration when the energy goes below a given value eth
. This is also a parameter of the functions CME_KSAT
and CDA_KSAT
(eth::Float64=1e6
).
The following shows how to save some quantity at each integration step. This is implemented via a SavingCallback
(see here). The general form of a function to be used as a callback is
function save_something(u, t, integrator)
# compute something
return #something
end
Thus, the package implements certain predefined requests for the numerical integrator.
The user could:

get_graph(integrator)
returns thegraph::HGraph
variable 
reshape_u_to_probs_CME(u::Vector{Float64}, integrator)
reshapes the vector$u$ to the probabilities involved in the CME: 'p_cav' and 'probi' (see Section 4). 
reshape_u_to_probs_CDA(u::Vector{Float64}, integrator)
reshapes the vector$u$ to the probability involved in the CDA: 'p_joint' (see Section 4). 
get_pju_CME(u::Vector{Float64}, integrator)
computes the probability of the unsatisfied configuration for each hyperedge in the CME (see Section 4). 
get_pju_CDA(u::Vector{Float64}, integrator)
computes the probability of the unsatisfied configuration for each hyperedge in the CDA (see Section 4).
With this, the user can define a callback that, for example, saves the system's energy at each step of the CDA's integration:
function save_ener_CDA(u, t, integrator)
p_joint_u = get_pju_CDA(u, integrator)
e = ener(p_joint_u)
return e
end
where ener()
is simply
function ener(p_joint_u::Vector{Float64})
e = 0.0
for he in eachindex(p_joint_u)
e += p_joint_u[he]
end
return e
end
The numerical integration that saves the energy at each step can be performed by running:
using ApproxMasEq
using OrdinaryDiffEq, DiffEqCallbacks, Sundials
N = 100
alpha = 3.5
K = 3
eta = 0.7
rargs = [eta]
t0 = 0.0
tlim = 10.0
saved_eners_CDA = SavedValues(Float64, Float64)
cb_ener_CDA = SavingCallback(save_ener_CDA, saved_eners_CDA)
cbs_save_CDA = CallbackSet(cb_ener_CDA)
answ = CDA_KSAT(rate_FMS_KSAT, rargs, build_args_rate_FMS, N=N, K=K, alpha=alpha, tspan=[t0, tlim], cbs_save=cbs_save_CDA)
where the callback is passed as a parameter. The function receives a CallbackSet
, so different callbacks can be passed at once, via the argument
cbs_save::CallbackSet=CallbackSet()
. In the line cbs_save_CDA = CallbackSet(cb_ener_CDA)
the custom callback is cast to an object of type CallbackSet
and then put as an argument of CDA_KSAT
.
Something analogous can be done for the function CME_KSAT
[1] E. Aurell, G. D. Ferraro, E. Domínguez, and R. Mulet. A cavity master equation for the continuous time dynamics of discrete spins models. Physical Review E, 95:052119, 2017.
[2] E. Aurell, E. Domínguez, D. Machado and R. Mulet. Exploring the diluted ferromagnetic pspin model with a cavity master equation. Physical Review E, 97:050103(R), 2018
[3] E. Aurell, E. Domínguez, D. Machado and R. Mulet. A theory nonequilibrium local search on random satisfaction problems. Physical Review Letters, 123:230602, 2019
[4] D. Machado, R. Mulet and F. RicciTersenghi. Improved meanfield dynamical equations are able to detect the twostep relaxation in glassy dynamics at low temperatures. Journal of Statistical Mechanics: Theory and Experiment, 2023:123301