# ApproxMasterEqs

## 1. Installation

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/d4v1d-cub/ApproxMasterEqs.jl.git")

## 2. Description

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 $t$.

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=1e-6,
cbs_save::CallbackSet=CallbackSet(), dt_s::Float64=0.1, abstol::Float64=1e-6, reltol::Float64=1e-3)

In its first version, the package works for models on hypergraphs where only one configuration in the factor nodes unsatisfies the interaction (K-SAT like). This contains as a particular case a model with pairwise interactions (2-SAT like)

## 3. How to provide a graph?

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ös-Ré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: $N$, $K$, and $\alpha$. These, together with the optional seed_g (seed for the random generator), are used to create an Erdös-Rényi hypergraph with size $N$, $K$ nodes per hyperedge and mean node connectivity $c=\alpha K$.

The couplings or links for the K-SAT boolean formula can be input via the parameter `links::Matrix{Float64}`

. The indexes are `links[he, i]`

with $he=1, \ldots, M$ and $i=1, \ldots, K$, where $M$ is the number of hyperedges in the graph. If the user does not input any matrix of links, the default is to randomly choose one.

For a neat example see the end of the Section 5

## 4. How to provide a model?

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 non-constant 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 $(\sigma=1)$ and `Em`

is the local energy when the variable is negative $(\sigma=-1)$.

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} # single-site 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 K-SAT optimization

# This function computes the energy in a K-SAT 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

## 5. How to numerically integrate?

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 $N=10$ nodes, $K=3$ node per hyperedge and mean connectivity $c=\alpha K = 6$ is randomly built with seed_g=1.

The model is the one in the example of Section 4: Focused Metropolis Search algorithm in K-SAT. The constant argument for the transition rates is eta=1 ($\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) = 1-p(-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=1e-6`

absolute tolerance for the numerical integration with OrdinaryDiffEqs.`reltol::Float64=1e-6`

relative tolerance for the numerical integration with OrdinaryDiffEqs.

## 6. Accessing the results

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=1e-6`

).

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 pre-defined requests for the numerical integrator.

The user could:

`get_graph(integrator)`

returns the`graph::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`

## References

[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 p-spin 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 non-equilibrium local search on random satisfaction problems. Physical Review Letters, 123:230602, 2019

[4] D. Machado, R. Mulet and F. Ricci-Tersenghi. Improved mean-field dynamical equations are able to detect the two-step relaxation in glassy dynamics at low temperatures. Journal of Statistical Mechanics: Theory and Experiment, 2023:123301