Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add references to other Julia libraries #395

Merged
merged 5 commits into from
Jan 23, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions paper/paper.tex
Original file line number Diff line number Diff line change
Expand Up @@ -366,13 +366,13 @@ \subsection{Queuing methods for multivariate processes}

\section{Implementation} \label{sec:implementation}

\texttt{JumpProcesses.jl} is a Julia library for simulating jump --- or point --- processes which is part of Julia's SciML organization. Jumps are handled via callbacks that are checked at the end of each time-step of some time evolution algorithm, \eg an ODE solver from \texttt{OrdinaryDiffEq.jl}, a stochastic differential equation solver from \texttt{StochasticDiffEq.jl}, or the pure-jump process \texttt{SSAStepper} provided by \texttt{JumpProcesses.jl}. In simple terms, callbacks are functions that can be arbitrarily called at each step of the main loop of a time-stepping method.
\texttt{JumpProcesses.jl} is a Julia library for simulating jump --- or point --- processes which is part of Julia's SciML organization. In the Julia ecosystem, there are other libraries that can sample certain TPPs including \path{Hawkes.jl} \footnote{\url{https://github.com/em1234321/Hawkes.jl}}, \path{HawkesProcesses.jl}~\footnote{\url{https://github.com/dm13450/HawkesProcesses.jl}}, \path{NetworkHawkesProcesses.jl} \footnote{\url{https://github.com/cswaney/NetworkHawkesProcesses.jl}}, \path{PointProcessInference.jl} \cite{schauer2020} \footnote{\url{https://github.com/mschauer/PointProcessInference.jl}}, \path{GeoStats.jl} \cite{hoffimann2020} \footnote{\url{https://github.com/JuliaEarth/GeoStats.jl}}, \path{PiecewiseDeterministicMarkovProcesses.jl}~\cite{veltz2015}~\footnote{\url{https://github.com/rveltz/PiecewiseDeterministicMarkovProcesses.jl}}, and \path{PointProcesses.jl}~\cite{dalle2024} \footnote{\url{https://github.com/gdalle/PointProcesses.jl}}. Apart from \texttt{PiecewiseDeterministicMarkovProcesses.jl}, these other libraries can only sample the Poisson and/or the Hawkes processes. \texttt{PointProcesses.jl} also offers a formalized interface that other packages can implement to leverage its TPP modelling functionality. While \texttt{JumpProcesses.jl} can be used to directly simulate TPPs, in its documentation we also show how it can be wrapped to conform to this interface \footnote{\url{https://docs.sciml.ai/JumpProcesses/stable/applications/advanced_point_process}}.

Our discussion in Section~\ref{sec:act} identified three exact methods for simulating point processes. In all the cases, we identified two mathematical constructs required for simulation: the intensity rate and the mark distribution. In \texttt{JumpProcesses.jl}, these can be mapped to user defined functions \texttt{rate(u, p, t)} and \texttt{affect!(integrator)}. The former takes the current state of the system, \texttt{u}, user provided parameters, \texttt{p}, and the current time, \texttt{t}, and returns the value of the intensity function at this time. The latter takes the solver \texttt{integrator} object, which stores all solution information, and updates it, including the state \texttt{integrator.u}, for whatever changes should occur when the jump it encodes fires at the time \texttt{integrator.t}. The library provides APIs for defining processes based on the nature of the intensity rate and the intended simulation algorithm. Processes simulated using exact sampling methods can choose between \texttt{ConstantRateJump} and \texttt{VariableRateJump}. While the former expects the rate between jumps to be constant, the latter allows for time-dependent rates. The library also provides the \texttt{MassActionJump} API to define large systems of point processes that can be expressed as mass action type reaction equations. Finally, \texttt{RegularJump} is intended for tau-leaping methods.

The \textit{inverse} method as described around Equation~\ref{eqn:inverse} uses root finding to calculate the next jump time. Jumps to be simulated via the \textit{inverse} method must be initialized as a \texttt{VariableRateJump}. \texttt{JumpProcesses.jl} builds a continuous callback following the algorithm in~\cite{salis2005} and passes the problem to an \texttt{OrdinaryDiffEq.jl} integrator, which easily interoperates with \texttt{JumpProcesses.jl} (both libraries are part of the \textit{SciML} organization, and by design built to easily compose). \texttt{JumpProcesses.jl} does not yet support the CHV ODE based approach.

Alternatively, \textit{thinning} methods can be simulated via discrete steps. In the context of the library, any method that uses thinning via a discrete callback is called an \textit{aggregator}. There are twelve different aggregators which we discuss below and are summarized in Table~\ref{tab:aggregators} in the \hyperref[sec:annex]{Annex}. Aggregator's handle sampling the next jump time and type, which is then read via the callback by the user-selected time-stepper.
Alternatively, \textit{thinning} methods can be simulated via discrete steps. In \texttt{JumpProcesses.jl}, simulation approaches that take discrete steps are handled via discrete callbacks that are checked at the end of each time-step of some time evolution algorithm, \eg an ODE solver from \texttt{OrdinaryDiffEq.jl}, a stochastic differential equation solver from \texttt{StochasticDiffEq.jl}, or the pure-jump process \texttt{SSAStepper} provided by \texttt{JumpProcesses.jl}. In simple terms, discrete callbacks involve two functions. Condition functions are checked at each step of the main loop of a time-stepping method to see if the callback should be executed, and if it should, an associated affect function is called. In the context of the library, any method that uses thinning via a discrete callback is called an \textit{aggregator}. There are twelve different aggregators which we discuss below and are summarized in Table~\ref{tab:aggregators} in the \hyperref[sec:annex]{Annex}.

We start with constant rate \textit{thinning} aggregators for marked TPPs. Algorithm~\ref{algo:sim-thinning} assumes that there is a single process. In reality, all the implementations first assume a finite multivariate point process with \( M \) interdependent sub-processes. This can be easily conciliated, as we do now, using Definition 6.4.1~\cite{daley2003} which states the equivalence of such process with a point process with a finite space of marks.

Expand Down Expand Up @@ -526,7 +526,7 @@ \subsection{Benchmarks} \label{subsec:benchmark}

We fix the Hawkes parameters at \( \lambda = 0.5 , \alpha = 0.1 , \beta = 5.0 \) ensuring the process does not explode and simulate models in the range from \( 1 \) to \( 95 \) nodes for \( 25 \) units of time. We simulate \( 50 \) trajectories with a limit of ten seconds to complete execution. For this benchmark, we save the state of the system exactly after each jump.

We assess the benchmark in eight different settings. First, we run the \textit{inverse} method, \texttt{Coevolve} and \textit{CHV simple} using the brute force formula of the intensity rate which loops through the whole history of past events --- Equation~\ref{eqn:hawkes-brute}. Second, we simulate the same three methods with the recursive formula --- Equation~\ref{eqn:hawkes-recursive}. Next, we run the benchmark against \textit{CHV full}. All \textit{CHV} specifications are implemented with \texttt{PiecewiseDeterministicMarkovProcesses.jl}~\footnote{\url{https://github.com/rveltz/PiecewiseDeterministicMarkovProcesses.jl}} which is developed by Veltz, the author of the \textit{CHV} algorithm discussed in Subsection~\ref{subsec:sim-inverse}. Finally, we run the benchmark using the Python library Tick\footnote{\url{https://github.com/X-DataInitiative/tick}}. This library implements a version of the thinning method for simulating the Hawkes process and implements a recursive algorithm for computing the intensity rate.
We assess the benchmark in eight different settings. First, we run the \textit{inverse} method, \texttt{Coevolve} and \textit{CHV simple} using the brute force formula of the intensity rate which loops through the whole history of past events --- Equation~\ref{eqn:hawkes-brute}. Second, we simulate the same three methods with the recursive formula --- Equation~\ref{eqn:hawkes-recursive}. Next, we run the benchmark against \textit{CHV full}. All \textit{CHV} specifications are implemented with \texttt{PiecewiseDeterministicMarkovProcesses.jl} which is developed by Veltz, the author of the \textit{CHV} algorithm discussed in Subsection~\ref{subsec:sim-inverse}. Finally, we run the benchmark using the Python library Tick\footnote{\url{https://github.com/X-DataInitiative/tick}}. This library implements a version of the thinning method for simulating the Hawkes process and implements a recursive algorithm for computing the intensity rate.
gzagatti marked this conversation as resolved.
Show resolved Hide resolved

Table~\ref{tab:benchmark-hawkes} shows that the \textit{Inverse} method which relies on root finding is the most inefficient of all methods for any system size. For large system size this method is unable to complete all \( 50 \) simulation runs because it needs to find an ever larger number of roots of an ever larger system of differential equations.

Expand Down
Loading
Loading