text:
\documentclass{jors}

%% Set the header information
\pagestyle{fancy}
\definecolor{mygray}{gray}{0.6}
\renewcommand\headrule{}
\rhead{\footnotesize 3}
\rhead{\textcolor{gray}{UP JORS software Latex paper template version 0.1}}

\usepackage{graphics,color,graphicx,amsmath,hyperref}
\usepackage{xr,natbib}
\bibliographystyle{vancouver}
\usepackage{amsbsy}

\newcommand{\norm}[1]{\left\Vert#1\right\Vert}
\newcommand{\abs}[1]{\left\vert#1\right\vert}
\newcommand{\set}[1]{\left\{#1\right\}}
\newcommand{\eps}{\varepsilon}
%\bibliographystyle{abbrvnat}
\graphicspath{{images}}

\begin{document}

% COMMANDS -------------------------------------------------------
\newcommand{\mr}[1]{\mathrm{#1}}
\newcommand{\mb}[1]{\mathbf{#1}}
\newcommand{\br}[1]{\left<#1\right>}
\newcommand{\bl}[1]{\left|#1\right|}
\newcommand{\mc}[1]{\mathcal{#1}}
\newcommand{\tb}[1]{\textcolor{blue}{#1}}
\newcommand{\tr}[1]{\textcolor{red}{#1}}
\newcommand{\tg}[1]{\textcolor{green}{#1}}
\newcommand{\si}[0]{{\rm s}_{\rm i}}
\newcommand{\sj}[0]{{\rm s}_{\rm j}}
\newcommand{\bs}[1]{\boldsymbol{#1}}
\newcommand{\rs}[0]{{\rm s}}
\newcommand{\rk}[0]{{\rm k}}

\rule{\textwidth}{1pt}

\section*{(1) Overview}
\vspace{0.5cm}

\section*{Convenient interface to inverse Ising (ConIII): A Python 3 package for solving Ising-type maximum entropy models}

\section*{Paper Authors}
1. Lee, Edward D.; \\
2. Daniels, Bryan C.

\section*{Paper Author Roles and Affiliations}
1. Department of Physics, 142 Sciences Dr, Cornell University, Ithaca NY 14853 \\
2. ASU--SFI Center for Biosocial Complex Systems, Arizona State University, Tempe, AZ 85287

\section*{Abstract}

ConIII (pronounced CON-ee) is an open-source Python project providing a simple interface to solving the pairwise and higher order Ising model and a base for extension to other maximum entropy models. We describe the maximum entropy problem and give an overview of the algorithms that are implemented as part of ConIII (\url{https://github.com/eltrompetero/coniii}) including Monte Carlo histogram, pseudolikelihood,  minimum probability flow, a regularized mean field method, and a cluster expansion method. Our goal is to make a variety of maximum entropy techniques accessible to those unfamiliar with the techniques and accelerate workflow for users.

\section*{Keywords}
inverse Ising, maxent, maximum entropy, inference

\section*{Introduction}
Many biological and social systems are characterized by collective behavior: the correlated pattern of neural firing \cite{Schneidman:2006he}, protein diversity in the immune system \cite{Mora:2010jxb}, conflict participation in monkeys \cite{Daniels:2017cq}, flocking in birds \cite{Bialek:2012cs}, statistics of letters in words \cite{Stephens:2010hi}, or consensus voting in the US Supreme Court \cite{Lee:2015ev,Lee:ky}. Statistical physics is a natural approach to probing such systems precisely because they are collective \cite{Castellano:2009ce}.
Recently, the development of numerical, analytic, and computational tools have made it feasible in these large collective systems to solve for the maximum entropy (maxent) model that reproduces system behavior, corresponding to solving an ``inverse problem.''
This approach contrasts with the typical problem in statistical physics where one postulates the microscopic model (the Hamiltonian) and works out the physical behavior of the system. In the inverse problem, we  find the parameters that correspond to observed behavior of a known system. In many cases, this is a very difficult problem to solve and does not have an analytical solution, and we must rely on analytic approximation and numerical techniques to estimate the parameters.

The pairwise maxent model, the Ising model, has been of particular interest because of its simplicity and generality. A variety of algorithms have been proposed to solve the inverse Ising problem, but different approaches are disparately available on separate code bases in different coding languages, which makes comparison difficult and pedagogy more complicated. With ConIII, it is possible to solve the inverse Ising problem with a variety of algorithms in just a few lines of code.

ConIII is intended to provide a centralized resource for the inverse Ising problem and easy extension to other maxent problems. Although some of the implemented algorithms are specific to the pairwise Ising model, maxent models with arbitrary combinations of higher order constraints can be solved as well by specifying the particular constraints of the maxent model of interest.
%The modular structure of ConIII was inspired by the machine learning package scikit-learn, where each algorithm is implemented in a class.
%There are a variety of techniques for approaching the inverse problem|often in context of a particular model|but they are often presented as different approaches without clear relations to one another. Instead of focusing on the details of each technique separately, we relate these techniques to one another and point out when one would fail or be difficult to apply.

%We hope to provide a (soft) introduction to the analytic, numerical, and computational techniques used to solve these problems to make them accessible to students especially at the graduate level in physics and beyond (and even advanced undergraduates). With the goal of furthering general understanding, we aim not to be rigorous but to be accurate and comprehensible, pointing to the relevant literature whenever possible. Although we strive to be as clear and explicit as possible across this guide, a good background in mathematics or physics will be extremely helpful. In some places, we have left out steps in the derivations, and we encourage the reader to work out the missing steps.

In the first few sections of this paper, we give a brief overview of maxent and describe at a high level the algorithms implemented in this package. For those unfamiliar with maxent, we also provide some useful references like Ref~\cite{Bialek:2012ueb} and the appendix of \cite{Lee:2015ev}. For those seeking more detail about the implemented algorithms, we provide references specific to each algorithm section. Then, we describe the architecture of the package and how to contribute.

\section*{What is maximum entropy?}
Shannon introduced the concept of information entropy in his seminal paper about communication over a noisy channel \cite{Shannon:1948wk}. Information entropy is the unique measure of uncertainty that follows from insisting on elementary principles of consistency. According to Shannon, the entropy over the probability distribution $p(\rs)$ of possible discrete configurations $\mathcal S$ of a system is
\begin{align}
	S[p] &= -\sum_{{\rm s}\in \mathcal{S}} p({\rm s}) \log p({\rm s}).
\end{align}
These configurations could be on-off patterns of firing in neurons, the arrangement of letters in a word, or the orientation of spins in a material.

When there is no structure in the distribution, meaning that the probability is uniform, entropy is at a maximum. In the context of communication theory as Shannon first discussed, this means that there is no structure to exploit to make a prediction about the next part of an incoming message; thus, maximum entropy means that each new part of the message is maximally ``surprising.'' At the other extreme, when the message consists of the same bit over and over again, we can always guess at the following part of the message and the signal has zero entropy. In the context of modeling, we use entropy not to refer to the difficulty of the message, but to our state of knowledge about it. Entropy precisely measures our uncertainty about the configuration in which we expect to find the system.

Maximum entropy, or maxent, is the formal framework for building models that are consistent with statistics from the data but otherwise as structureless as possible \cite{Bretthorst:2003ua,Jaynes:1957fy}.
We begin by choosing a set of $K$ useful or important features from the data $f_{\rm k}(\rs)$ that should be true for the model that we are trying to build. These could be whether or not a set of neurons fire together in a temporal bin or the pairwise coincidence for primates in a conflicts. The average of this feature across the data set $\mathcal{D}$ with $R$ samples is
\begin{align}
	\br{f_{\rm k}}_{\rm data} &= \frac{1}{R}\sum_{{\rm s}\in \mathcal{D}} f_{\rm k}({\rm s}).
\end{align}
According to the model in which each observation s occurs with some probability $p(\rm s)$, the same average is calculated over all possible states
\begin{align}
	\br{f_{\rm k}} &= \sum_{\rm s\in\mathcal{S}} p(\rs)f_{\rm k}(\rs).
\end{align}

We assert that the model should fit the $K$ features while maximizing entropy. The standard procedure is to solve this by the method of Langrangian multipliers. We construct the Langrangian functional $\mathcal{L}$ by introducing the multipliers $\lambda_{\rm k}$.
\begin{align}
	\mathcal{L}[p] &= -\sum_{\rm s\in\mathcal{S}} p(\rs)\log p(\rs) - \sum_{\rm k=1}^{K} \lambda_{\rm k} \big(\br{f_{\rm k}}-\br{f_{\rm k}}_{\rm data}\big)\label{eq:Lagrang}
\end{align}

%In the notation of statistical physics, the Langrangian is the Helmholtz free energy describing the competition between entropy and the structure described in the Hamiltonian $E$ in equilibrium (See Appendix) \footnote{If we fix the average energy of the system, we get the microcanonical ensemble (See Appendix). which is equivalent to the axiom that all microstates have equal probability (that is the maxent distribution).}.
%\begin{align}
%	F &= S - \br{E}\\
%	E &= -\sum_{\rm k}^K\lambda_{\rm k}f_{\rm k}(s)
%\end{align}
%This formulation makes clear the fundamental connection that statistical mechanics is an inference procedure using the maximum entropy principle \cite{Jaynes:1957fy}.

Then, we solve for the fixed point by taking the derivative with respect to $\lambda_{\rm k}$.
The resulting maxent model is a Boltzmann distribution over states:
\begin{align}
    \label{eq:energy}
	p(\rs) &= \left.e^{-E(\rs)}\right/Z,\\
\intertext{with relative negative log-likelihood (also known as the energy or Hamiltonian)}
E(\rs) &= -\sum_{\rm k=1}^K\lambda_{\rm k}f_{\rm k}(\rs),
\intertext{and normalization factor (also known as the partition function)}
	Z &= \sum_{\rm s\in\mathcal{S}} e^{-E(\rs)}.
\end{align}
Entropy is a convex function of $p$ and the constraints are linear with respect to $p$, so the problem is convex and the maxent distribution unique.
Readers familiar with statistical physics will recognize this as an alternative derivation of the microcanonical ensemble, demonstrating that statistical mechanics can be viewed as an inference procedure using the maxent principle \cite{Jaynes:1957fy}.

Finding the parameters $\lambda_{\rm k}$ that match the constraints $\br{f_{\rm k}}_{\rm data}$ is equivalent to minimizing the Kullback-Leibler divergence between the model and the data \cite{Cover:2006tl}
\begin{align}
	D_{\rm KL}(p_{\rm data}||p) &= \sum_\rs p_{\rm data} \log\left(\frac{p_{\rm data}(\rs)}{p(\rs)}\right)\\
	\frac{\partial D_{\rm KL}}{\partial \lambda_{\rm k}} &= \sum_{\rs} p_{\rm data}(\rs) \frac{\partial (-E(\rs)-\log Z)}{\partial \lambda_{\rm k}} = 0 \nonumber \\
	\implies  \br{f_{\rm k}}_{\rm data} &= \br{f_{\rm k}}\label{eq:fk-fk}.
\end{align}
 In other words, the parameters of the maxent model are the ones that minimize the information theoretic ``distance'' to the distribution of the data given the constraints. Note that these parameters are given by the data: once the constraints have been chosen, there is a single maxent solution, with no free parameters.

\subsection*{The Ising model}
The Ising model is a statistical physics model of magnetism \cite{Ising:1924vf}. It consists of a set of spins $\{\rs_{\rm i}\}$ with 2 possible orientations (up and down), each responds to its own external magnetic field $h_{\rm i}$ and each pair is coupled to each other with pairwise coupling $J_{\rm ij}$. The strength of the magnetic field determines the tendency of each of the spins to orient in a particular direction and the couplings determine whether the spins tend to point together ($J_{\rm ij}>0$) or against each other ($J_{\rm ij}<0$). Typically, neighbors are defined as spins that interact with one another given by some underlying network structure. Figure \ref{gr:ising} shows a fully-connected example.

\begin{figure}[tb]\centering
	\includegraphics[width=.65\linewidth,clip,trim={100 70 70 60}]{ising_example.pdf}
\caption{Example of a fully connected pairwise Ising model with positive and negative couplings. Each spin $\rs_{\rm i}$ (circle) can take one of two states (black or white, corresponding to $-1$ and $1$) and is connected to every other spin in the system with a positive (red) or negative (blue) coupling. These states could describe the on-off patterns of firing in neurons, the orientation of spins in a material, or if each spin is no longer binary the arrangement of letters in a word (a Potts model).}
\label{gr:ising}
\end{figure}

The energy of each configuration determines its probability via Eq \eqref{eq:energy},
\begin{align}
	E(\rs) &= -\sum_{\rm i<j}^N J_{\rm ij}\si\sj -\sum_{\rm i=1}^Nh_{\rm i}\si,
\end{align}
such that lower energy states are more probable.
%States with lower energy are more likely so spin tend to orient along the direction of the local field $h_{\rm i}^{\rm local} = \sum_{\rm j}J_{\rm ij}\sj +h_{\rm i}$.

We can derive the Ising model from the perspective of maxent. Fixing the the means and pairwise correlations to those observed in the data
\begin{align}
	\br{\si} &= \br{\si}_{\rm data} \\
	\br{\si\sj} &= \br{\si\sj}_{\rm data}
\end{align}
we go through the procedure of constructing the Langrangian from Eq~\ref{eq:Lagrang}
\begin{align}
	\mathcal{L}[p] &= -\sum_\rs p(\rs)\log p(\rs) +\sum_{\rm i<j}^N J_{\rm ij}\br{\si\sj} +\sum_{\rm i}^Nh_{\rm i}\br{\si}\\
	\frac{\partial\mathcal{L}[p]}{\partial p(\rs)} &= -\log p(\rs)-1 +\sum_{\rm i<j}^N J_{\rm ij}\si\sj +\sum_{\rm i}^Nh_{\rm i}\si\\
	\log p(\rs) &= -1 +\sum_{\rm i<j}^N J_{\rm ij}\si\sj +\sum_{\rm i}^Nh_{\rm i}\si\label{eq:E-1}\\
	p(\rs) &= \left.e^{-E(\rs)}\right/Z\\
\intertext{where the $-1$ in Eq \ref{eq:E-1} has been absorbed into the normalization factor}
	Z &= \sum_{\rs} e^{-E(\rs)}.
\end{align}
such that the probability distribution is normalized $\sum_{\rs} p(s)=1$.
Thus, the resulting model is exactly the Ising model mentioned earlier.

Despite the simplicity of the Ising model, the structure imposed by the discrete nature of the spins means that finding the parameters is challenging analytically and computationally. In the last few years, numerous techniques have been suggested for solving the inverse Ising problem exactly or approximately \cite{Nguyen:2017ww}. We have implemented some of them in ConIII and designed a package structure to make it easily extensible to include more methods.  Here, we briefly describe the algorithms that are part of the first official version of the package. The goal is to give the user a sense of how they work without getting bogged down in heavy detail. For more detail, we suggest perusing the papers referenced in each section or the review \cite{Nguyen:2017ww}. For a complete beginner, it may be useful to first get familiar with a slower introduction like in the Appendices of Ref \cite{Lee:2015ev}, Ref \cite{Bialek:2012ueb}, or Ref \cite{Bretthorst:2003ua}.

\section*{Inverse Ising methods implemented in ConIII}
\subsection*{Enumeration}
The na\"{i}ve approach that only works for small systems is to write out the equations from Eq~\ref{eq:fk-fk} and solve them numerically. After writing out all $K$ equations,
\begin{equation}
	\br{f_{\rm k}} = -\frac{\partial \ln Z}{\partial \lambda_{\rm k}} = \br{f_{\rm k}}_{\rm data} \\,
\end{equation}
we can use any standard root-finding algorithm to find the parameters $\lambda_{\rm k}$.
This approach, however, involves enumerating all states of the system, whose number grows exponentially with system size.
%[Simple example with logical gates?]

For the Ising model, writing down the equations has a number of steps $\mathcal{O}(K^2 2^{N})$, where $K$ is the number of constraints and $N$ the number of spins.  Each evaluation of the objective in the root-finding algorithm will be of the same order. For relatively small systems, around $N\leq15$, this approach is feasible on a typical desktop computer and is a good way to test the results of a more complicated algorithm.

This approach is implemented by the {\tt Enumerate} class.

\subsection*{Monte Carlo Histogram (MCH)}
Perhaps the most straightforward and most expensive computational approach is Monte Carlo Markov Chain (MCMC) sampling.  A series of states sampled from a proposed $p(\rs)$ is produced by MCMC to approximate $\br{f_{\rm k}}$ and determine how close we are to matching $\br{f_{\rm k}}_{\rm data}$.  The parameters are then adjusted using a learning rule, and both sampling and learning are repeated until a stopping criterion is met.
%The magic behind MCMC approaches is that only relative probabilities of states need to be known to return a sample of the distribution. Thus, we do not need to calculate the partition function, but can instead compute the difference in energy between two states.
This can be combined with a variety of approximate gradient descent methods to reduce the number of sampling steps by predicting how the distribution will change if we modify the parameters slightly.
The particular technique implemented in ConIII is the Monte Carlo Histogram (MCH) method \cite{Broderick:2007wq}.

Since the sampling step is expensive, the idea behind MCH is to reuse a sample for more than one gradient descent step \cite{Broderick:2007wq}. Given that we have a sample with probability distribution $p(\rs)$ generated with parameters $\lambda_{\rm k}$, we would like to estimate the proposed distribution $p'(\rs)$ from adjusting our parameters $\lambda_{\rm k}' = \lambda_{\rm k}+\Delta\lambda_{\rm k}$. We can leverage our current sample to make this extrapolation.
\begin{align}
	p' &= \frac{p'}{p}p\\
	p'(\rs)	&= \frac{Z}{Z'}e^{\sum_\rk \Delta\lambda_k f_\rk(\rs)} p(\rs)\\
\intertext{To estimate the average,}
	\sum_\rs p'(\rs) f_{\rm k}(\rs) &= \frac{Z}{Z'} \sum_\rs p(\rs) e^{\sum_{\rm k} \Delta\lambda_{\rm k} f_{\rm k}(\rs)} f_{\rm k}(\rs)\\
\intertext{To be explicit about the fact that we only have a sampled approximation to $p$, we replace $p$ with the sample distribution.}
	\br{f_{\rm k}}' &= \frac{Z}{Z'} \br{e^{\sum_{\rm k} \Delta\lambda_{\rm k} f_{\rm k}(\rs)} f_{\rm k}(\rs)}_{\rm sample}
\end{align}
Likewise, the ratio of the partition function can be estimated
\begin{align}
	\frac{Z}{Z'} \approx 1\left/\br{e^{\sum_{\rm k} \Delta\lambda_{\rm k} f_{\rm k}(\rs)}}_{\rm sample}\right.
\end{align}

At each step, we update the Lagrangian multipliers $\{\lambda_\rk\}$ while being careful to stay within the bounds of a reasonable extrapolation. One suggestion is to update the parameters with some inertia \cite{Tkacik:2006vq}
\begin{align}
	\Delta\lambda_\rk(t+1) &= \Delta \lambda_\rk(t) + \epsilon \Delta\lambda_\rk(t-1)\label{eq:mch learn1}\\
	\Delta \lambda_\rk(t) &= \eta\left(\br{f_\rk}'-\br{f_\rk}\right)\label{eq:mch learn2}
\end{align}
This has a fixed point at the correct parameters.
%Another heuristic is to begin with a small sample size and increase the size of the sample set as we refine our estimates of the parameters.

In practice, MCH can be difficult to tune properly and one must check in on the progress of the algorithm often. One issue is choosing how to set the learning rule parameters $\eta$ and $\epsilon$. One suggestion for $\eta$ is to shrink it as the inverse of the number of iterations \cite{Tkacik:2006vq}. Another issue is that parameters cannot be changed by too much when using the MCH approximation step or the extrapolation to $\lambda_{\rm k}'$ will be inaccurate and the algorithm will fail to converge. In ConIII, this can be controlled by setting a bound on the maximum possible change in each parameter  $\Delta\lambda_{\rm max}$ and restricting the norm of the vector of change in parameters $\sum_{\rm k} \sqrt{\Delta\lambda_{\rm k}^2}$. Another issue is setting the parameters of the MCMC sampling routine. Both the burn time (the number of iterations before starting to sample) and sampling iterations (number of iterations between samples) must be large enough that we are sampling from the equilibrium distribution.  Typically, these are found by measuring how long the energy or individual parameter values remain correlated as MCMC progresses. The parameters may need to be updated during the course of MCH because the sampling parameters may need to change with the estimated parameters of the model. For some regimes of parameter space, samples are correlated over long times and alternative sampling methods like Wolff or Swendsen-Wang would vastly reduce time to reach the equilibrium distribution although these are not included in the current release of ConIII. We do not discuss these sampling details here, but see Refs \cite{MacKay:2005wc,Newman:1999wu} for examples.

The main computational cost for MCH lies in the sampling step. For each iteration of MCH, the runtime is proportional to the number of samples $n$, number of MCMC iterations $T$, and the number of constraints for the Ising model $N^2$,
$\mathcal{O}(T n N^2)$, whereas the MCH estimate is relatively quick $\mathcal{O}(t n N^2)$ because the number of MCH approximation steps needed to converge is much smaller than the number of MCMC sampling iterations $t<<T$.
%What is the runtime for the learning rules Eqs \ref{eq:mch learn1} and \ref{eq:mch learn2}? This a much harder question that gets into optimization theory.

MCH is implemented in the {\tt MCH} class.

\subsection*{Pseudolikelihood}
The pseudolikelihood approach is an analytic approximation to the likelihood that drastically reduces the computational complexity of the problem and is exact as $N \rightarrow \infty$ \cite{Aurell:2012hi}. We calculate the conditional probability of each spin $\rs_{\rm i}$ given the rest of the system $\{\rs_{\rm j\neq i}\}$
\begin{align}
	p\left(\rs_{\rm i}|\{\rs_{\rm j\neq i}\}\right) &= \left( 1+e^{-2\rs_{\rm i} \left(h_{\rm i}+\sum_{j\neq i}J_{\rm ij}\rs_{\rm j}\right)} \right)^{-1}
\end{align}
Taking the logarithm, we define the approximate log-likelihood by summing over data points indexed by r:
\begin{align}
	f(h_{\rm i},\{J_{\rm ij}\}) &= \sum_{\rm r=1}^R \ln p\left(\left.\rs_{\rm i}^{(\rm r)}\right|\{\rs_{\rm j\neq i}\}^{(\rm r)}\right).
\end{align}
In the limit where the ensemble is well sampled, the average over the data can be replaced by an average over the ensemble:
\begin{align}
	f(h_{\rm i},\{J_{\rm ij}\}) &= \sum_{\rm s} \ln p\left(\left.\rs_{\rm i}\right|\{\rs_{\rm j\neq i}\}\right)p\left(\rs;h_{\rm i},\{J_{\rm ij}\}\right).
\end{align}

To find the point of maximum likelihood for a single spin $\rs_{\rm i}$, we calculate the analytical gradient and Hessian, $\partial f/\partial J_{\rm ij}$ and $\partial^2 f/\partial J_{\rm ij}\partial J_{\rm i'j'}$ for a Newton conjugate-gradient descent method. After maximizing likelihood for all spins, the maximum likelihood parameters may not satisfy the symmetry $J_{\rm ij}=J_{\rm ji}$. We impose the symmetry by insisting that
\begin{align}
	J_{\rm ij}' &= (J_{\rm ij}+J_{\rm ji})/2.
\end{align}

Pseudolikelihood is extremely fast and often surprisingly accurate. Each calculation of the gradient is order $\mathcal{O}(RN^2)$ and Hessian $\mathcal{O}(RN^3)$, which must be done for all $N$. With analytic forms for the gradient and Hessian, the conjugate-gradient descent method tends to converge quickly.

Pseudolikelihood for the Ising model is implemented in {\tt Pseudo}.

\subsection*{Minimum Probability Flow (MPF)}
Minimum probability flow involves analytically approximating how the probability distribution \textit{changes} as we modify the \textit{configurations} \cite{Sohl-Dickstein:2009tt,SohlDickstein:2011im}. In the methods so far mentioned, the approach has been to maximize the objective (the likelihood function) by immediately taking the derivative with respect to the parameters. With MPF, we first posit a set of dynamics that will lead the data distribution to equilibrate to that of the model. When these distributions are equivalent, then there is no ``probability flow'' between them. This technique is closely related to score matching, where we instead have a continuous state space and can directly take the derivative with respect to the states without specifying dynamics \cite{Hyvarinen:2007ed}.

First note that Monte Carlo dynamics (satisfying ergodicity and detailed balance) would lead to equilibration to the stationary distribution. One such transition matrix suggested in Ref \cite{SohlDickstein:2011im} is
\begin{align}
	\dot{p}_{\rm s} &= \sum_{\rm s'\neq s} \Gamma_{\rm ss'} p_{\rm s'} -\sum_{\rm s'\neq s} \Gamma_{\rm s's} p_{\rm s}\\
	\Gamma_{\rm ss'} &= g_{\rm ss'}\exp\left[ \frac{1}{2}\left( E_{\rm s'}-E_{\rm s} \right) \right]
\end{align}
with transition probabilities $\Gamma_{\rm ss'}$ from state ${\rm s'}$ to state s. The connectivity matrix $g_{\rm ss'}$ specifies whether there is edge between states s and ${\rm s'}$ such that probability can flow between them. By choosing a sparse $g_{\rm ss'}$ while not breaking ergodicity, we can drastically reduce the computational cost of computing this matrix.

Imagine that we start with the distribution over the states as given by the data and run the Monte Carlo dynamics. When data and model distributions are different, probability will flow between them and indicate that the parameters must be changed.  By minimizing a derivative of the Kullback-Leibler divergence, we measure how the difference between the model and the states in the data $\mathcal{D}$ changes when the dynamics are run for an infinitesimal amount of time.
\begin{align}
	L(\{\lambda_{\rm k}\}) \equiv \partial_t D_{\rm KL}(p^{(0)}||p^{(t)}\left(\{\lambda_{\rm k}\}\right)) &= \sum_{\rm s \not\in \mathcal{D}} \dot{p}_{\rm s}(\lambda_{\rm k})
\end{align}
The idea is that this derivative is also minimized with optimal parameters: the MPF algorithm looks for a minimum of the objective function $L$.

%MPF satisfies a number of useful properties:

For the Ising model, each evaluation of the objective function where $\Gamma_{\rm s s'}$ connects each data state with $G$ neighbors has runtime $\mathcal{O}(RGN^2)$. In a large fully connected system, $G\sim 2^N$ would be prohibitively large so a sparse choice is necessary.

MPF is implemented in the {\tt MPF} class.

\subsection*{Regularized mean-field method}
One attractively simple and efficient approach uses a regularized version of
mean-field theory.  In the inverse Ising problem, mean-field theory is equivalent
to treating each binary individual as instead having a continuously varying state
(corresponding to its mean value).  The inverse problem then turns into simply inverting
the correlation matrix $C$ \cite{Monasson:2011fo}:
\begin{equation}
\label{meanFieldSolution}
J^{\mathrm{mean-field}}_{\rm ij} =
    - \frac{ (C^{-1})_{\rm ij} }{ \sqrt{p_{\rm i}(1-p_{\rm i})p_{\rm j}(1-p_{\rm j})} },
\end{equation}
where
\begin{equation}
C_{\rm ij} = \frac{ p_{\rm ij} - p_{\rm i} p_{\rm j} }{ \sqrt{p_{\rm i}(1-p_{\rm i})p_{\rm j}(1-p_{\rm j})} },
\end{equation}
and where $p_{\rm i}$ corresponds to the frequency of individual i being
in the active ($+1$) state and $p_{\rm ij}$ is the frequency of the pair
i and j being simultanously in the active state.

A simple regularization scheme in this case is to discourage large values in the interaction
matrix $J_{\rm ij}$.  This corresponds to putting more weight on solutions that are closer to
the case with no interactions (independent individuals).  A particularly convenient form
adds the following term, quadratic in $J_{\rm ij}$, to the negative log-likelihood:
\begin{equation}
\gamma \sum_{\rm i} \sum_{\rm i<j} J_{\rm ij}^2 p_{\rm i} (1-p_{\rm i}) p_{\rm j} (1-p_{\rm j}).
\end{equation}
In this case, the regularized version of the mean-field solution in \eqref{meanFieldSolution}
can be solved analytically, with the slowest computational step coming from the inversion
of the correlation matrix.  For details, see Refs.~\cite{Daniels:2017cq,Barton:2013fja}.

The idea is then to vary the regularization strength $\gamma$ to move between the
non-interacting case ($\gamma \rightarrow \infty$) and the naively calculated
mean-field solution \eqref{meanFieldSolution} ($\gamma \rightarrow 0$).
While there is no guarantee that varying this one parameter will produce solutions that are
good enough to ``fit within error bars,'' this approach has been successful in at least
one case of fitting social interactions \cite{Daniels:2017cq}.

The inversion of the correlation matrix is relatively fast, bounded by $\mathcal{O}(N^3)$. Finding the optimal $\gamma$, involves Monte Carlo sampling from the model distribution, which has computational cost similar to MCH. It is, however, much more efficient because we are only optimizing a single parameter.

This is implemented in {\tt RegularizedMeanField}.

\subsection*{Cluster expansion}

Adaptive cluster expansion \cite{Monasson:2011fo,Monasson:2011fo,Barton:2013fja}
iteratively calculates terms in the
cluster expansion of the entropy $S$:
\begin{equation}
S - S_0 = \sum_\Gamma \Delta S_\Gamma,
\end{equation}
where the sum is over clusters $\Gamma$ and in the exact case
includes all $2^N - 1$ possible nonempty subsets of individuals in the system.  In the simplest version of the expansion,
one expands around $S_0 = 0$.  In some cases it can be more advantageous to
expand around
the independent individual solution or one of the mean-field solutions
described in the previous section \cite{Barton:2013fja}.

The inverse Ising problem is solved independently
on each of the clusters, which can be done exactly when the
clusters are small.  These results are used to construct a full
interaction matrix $J_{\rm ij}$.
The expansion starts with small clusters and expands to use larger
clusters, neglecting any clusters whose
contribution $\Delta S_\Gamma$ to the entropy falls below a threshold.
To find the best solution that does not overfit,
the threshold is initially set at a large value and then lowered,
gradually including more clusters in the expansion, until samples from
the resulting $J_{\rm ij}$ fit the desired statistics of the data sufficiently well.

The runtime will depend on the size of clusters included in the expansion.
If the expansion is truncated at clusters of size $n$, the worst-case runtime
would be $\mathcal{O}\left(\binom{N}{n} 2^n\right)$.  The point is that $S$ can often
be accurately estimated even when $n \ll N$.

The adaptive cluster expansion method is implemented in the {\tt ClusterExpansion} class.

%\subsection*{Bethe/Kikuchi free energy/cavity methods}
%Another approach involves a cluster-expansion of the free energy, also known as the cavity method. This has not yet been implemented in ConIII.
%The cavity method involves use of the marginalized distribution.

%\section*{Samplers}
%For estimating statistics of a given probability distribution $p(\rs)$, we have implemented in ConIII two versions of the well-known Metropolis algorithm. One is specific to the Ising model ({\tt FastMCIsing}) and the other ({\tt Metropolis}) can sample a system given an arbitrary function that calculates the energy of each state.
%
%In some parameter regimes where spins are tightly correlated, the Metropolis algorithm is very inefficient because the time to flip large numbers of spins simultaneously is extremely long. Cluster sampling like Wolff or Swendsen-Wang are much more efficient and will be included in future releases of ConIII.

%\section*{The finite sample problem}
%A fundamental problem in model inference is that uncertainty coming
%from the finiteness of data translates into uncertainties in parameters.
%In the exposition of the maxent formulation derived from Eq~\ref{eq:Lagrang}, there is no ambiguity in the model parameters because they are fully specified by the constraints calculated from data. In reality, these constraints are typically estimated from a finite sample, and they are noisy.
%The straightforward answer to this problem is to take more data---in a pairwise
%maxent problem, we might insist that we have enough samples to well-constrain
%the correlation between every pair of individuals.  But it is not always possible
%to take enough data.  For instance, in a social system in which we are trying to
%measure stable social structure that lasts on the order of months, there are only
%a finite number of social interactions that occur over those months, which may
%not be enough to tightly constrain parameters. When we fit exactly to constraints measured from a small amount of data, we run the danger of overfitting and poor generalization.
%
%When we find the parameters that minimize the Kullback-Leilber divergence between the model and the data distributions, we are maximizing the likelihood of the data.
%Sample size influences curvature of the likelihood function around the maximum-likelihood peak.  Specifically, small sample size leads to a smaller curvature, implying a broad peak and larger uncertainty in parameters.
%We can estimate this uncertainty from the data: Assuming that the data is independently and identically distributed, the errors are given by the standard error of the mean of a binomial distribution \mbox{$p_{\rm ij...k} = p(s_{\rm i}=s_{\rm j}=...=s_{\rm k}=1)$} with $R$ data points.
%\begin{align}
%	\delta_{\rm ij...k} &= \sqrt{(1-p_{\rm ij...k})p_{\rm ij...k}/R}.
%\end{align}
%In algorithms that iteratively bring the model closer to the data, we can stop once we have gotten close enough, where close enough is given by the expected fluctuations in the data distribution.
%
%An alternative approach is to regularize the problem as in the regularized mean field or cluster expansion algorithms. Here, we restrict the search space in some principled way so that more complicated
%solutions are disallowed.  We then check that the regularized solutions fit the
%data within expected statistical fluctuations.  If not, a more lax regularization
%can be used to allow more complicated solutions that are able to fit the
%remaining signal in the data.
%In a Bayesian formulation, this approach is equivalent to including a prior distribution over  parameters that more heavily weights simpler models.
%%\begin{align}
%%    \log p({\rm model}|{\rm data}) &\propto \log p({\rm data}|{\rm model}) + \log p({\rm model})\label{eq:prior}
%%\end{align}
%
%\section*{Model Validation}
%Besides just fitting the model while accounting for the noisiness of the data, how do we know if we need a more complicated model? To answer this question, we need to quantify how well the model fits and the corresponding tradeoff between the complexity and fit of the model.
%
%One way to formalize this is to quantify the tradeoff between a good description of the data with the cost of describing the model. In the Bayesian formulation, we would write this as maximizing the posterior probability of the data given the parameters $\theta$,
%\begin{align}
%    \log p(\theta|\mathcal{D}) &= \log p(\mathcal{D}|\theta) +\log p(\theta) + {\rm const}
%\end{align}
%and the tradeoff manifests as a competition between the likelihood $p(\mathcal{D}|\theta)$ and prior $p(\theta)$. For example, in the regularized mean-field method, the quadratic prior on the couplings favors sparse solutions over more complex solutions that might capture the data better but fail to generalize.
%%in Eq~\ref{eq:prior}. Roughly speaking, we can imagine that specifying the $L$ parameters is localizing a region in a high-dimensional space where each dimension shrinks with the number of data points $K$ as $\sim K^{-1/2}$. The volume of the ``ball'' grows like $K^{-L/2}$, so that the information cost goes like $-\frac{L}{2}\log(K)$ \cite{Lee:2015ev}.
%This picture is formalized by quantities like the Akaike information criterion (AIC) or the Bayesian information criterion (BIC) \cite{Anonymous:mVL3xTtr}.
%
%Although likelihood tells us how well the models are doing relative to one another, it does not tell us in detail how the model is fitting the data.
%Basic checks would be to compare against other correlations in the data.  This could include higher-order correlations than those used to fit the model or other features that are relevant to the question at hand. For instance, an important coarse-grained feature in  macaque conflict is the distribution of fight sizes, which is well reproduced by fitting only pairwise correlations \cite{Daniels:2017cq}.
%% In Figure~??, we show an example of a pairwise maxent model tested against higher order correlations.
%The most strict check would be to compare the entire model probability distribution with that of the data \cite{Lee:2015ev}, but this is only feasible when the data set is reasonably large.
%
%A summary of how well the model captures the distribution across the entire probability distribution is the multi-information \cite{Nemenman:2004wz,Schneidman:2006he}. If we take a maxent model and add further constraints, the models can be ordered in terms of entropy $S_{\rm 1}>S_{2}>...S_{\rm m}>...>S_{\rm data}$, where the minimum entropy the most constrained model could have is equal to the data. To measure how much correlation in the data our model has captured, we can calculate the multi-information $I_{\rm m} = S_{\rm 1}-S_{\rm m}$, the amount of correlation captured by the model relative to the independent model. The fraction of multi-information captured is $F=I_{\rm m}/I_{\rm data}$. This is a measure of how much of the correlation in the data is captured by the model.\footnote{See Refs \cite{Bialek:2012ueb} and \cite{Lee:2015ev} for details and further references on how to estimate information quantities.}
%
%These tools to quantify model fit then inform the important choice of which constraints to impose. Typically, the approach is to constrain the lowest order interactions that are sufficient to produce collective behavior. The intuition from physics is from the observation that many physical systems are very well described by pairwise interactions \cite{Ranganathan:2007wz}.
%From the model fitting perspective, however, we might choose to constrain other parts of the probability distribution. Ref \cite{Ganmor:2011ct} explores choosing correlations to constrain depending on whether they are signicantly large. Ref \cite{Nemenman:2016kl} discusses how a pairwise model becomes an increasingly effective description of a system with higher order interactions as the system gets larger.

\section*{Implementation and architecture}
The package is divided into three principal modules containing the algorithms for solving the inverse maxent problem ({\tt solvers.py}), the Monte Carlo Markov Chain (MCMC) sampling algorithms ({\tt samplers.py}), and supporting ``utility'' functions for the remaining modules ({\tt utils.py}) as shown in Fig.~\ref{gr:architecture}.
Besides the {\tt utils.py} module, the package is organized around classes that correspond to different algorithms. This class-based structure ensures that the state of the solver or sampler, including the data it was fit to and the current guess for the parameters, are all contained within the instance of the algorithm class. As a result, the current state of work can be saved and moved between workstations using the Python package {\tt dill}.

For the solvers, the different algorithms available are accessible from the {\tt coniii.solvers} module as listed in Fig.~\ref{gr:architecture}. These algorithm classes are all derived from a base {\tt Solver} class as shown in Fig.~\ref{gr:architecture}. In the accompanying release version of the package, v1.0.3, the method {\tt Solver.solve} serves as the interface for solving the inverse maxent problem. To keep the solution algorithms generic enough to solve a variety of different maxent problems, they all require that the user define the maxent model upon instantiation through the definition of keyword arguments like {\tt calc\_observables}. The particular methods required to specify the maxent problem differ by algorithm, but for the pairwise maxent problem we have made it easy by defining those functions as part of the package. These helper functions are available as part of the {\tt utils.py} module and their use is demonstrated in the Jupyter notebook usage guide.

The MCMC sampling algorithms are likewise based on a class architecture derived from {\tt Sampler} as shown in Fig.~\ref{gr:architecture}. Each instance of {\tt Solver} automatically instantiates this class under {\tt Solver.sampler} and wraps calls to it. For the Ising model, this is an instance of {\tt Metropolis}. Other sampling algorithms listed in the {\tt samplers.py} box in Fig.~\ref{gr:architecture} will be released with later versions of this package.

\begin{figure}[tb]
	\includegraphics[width=\linewidth,clip=True,trim=70 450 60 0]{architecture}
	\caption{ConIII architecture. The principal modules are {\tt solvers.py}, {\tt samplers.py}, and {\tt utils.py}. The module {\tt solvers.py} contains classes based on {\tt Solver} that each implement a different algorithm for solving the relevant inverse maxent problem accessible through the method {\tt Solver.solve()}. The {\tt samplers.py} module contains the Metropolis algorithm for Monte Carlo Markov Chain sampling and will support other samplers in future versions (gray font) including Wolff sampling, Swendsen-Wang sampling, and parallel tempering. The {\tt utils.py} module contains supporting functions for the other modules such as the few examples listed. ConIII's modularized structure ensures that contributed algorithms can be appended independently of existing code.}
	\label{gr:architecture}
\end{figure}

\section*{Quality Control}
For checks of basic functionality, the package is released with unit tests that can be run with the Python package \href{https://pytest.org}{pytest}.

The most direct test of the algorithms is to generate a system where the parameters are known, sample from the system to generate a data set, and run the inverse solution to make sure that the correlations and parameters match the known values. With a finite sample, exact correspondence to the correct parameters is not expected although differences should decrease with a larger sample. Furthermore, most of the algorithms only return an approximate solution such that the fidelity of the found parameters to the original ones will depend on the sample size and whether or not the approximation is valid. The Jupyter notebook released with the software provides examples for using the algorithms included in ConIII for a random system of five spins. We recommend that the user run this notebook to check how well different algorithms converge to the solutions depending on the algorithm and sample size.

More importantly, the user can check if the algorithms match the expected correlations closely or not. How one checks the validity of a particular maxent model for data is beyond the scope of this paper, but we point the reader to the appendix of Ref~\cite{Lee:2015ev} where the methodology is explained in detail for a broad audience.

If there are any issues or bugs in the software, we organize improvements and patches through  the \href{https://github.com/eltrompetero/coniii}{GitHub repository} where both issues can be filed and pull requests made.

\section*{(2) Availability}
\vspace{0.5cm}
\section*{Operating system}
Linux, MacOS, Windows

\section*{Programming language}
Python 3.6, 3.7

\section*{Dependencies}
Python packages multiprocess $\geq$v0.70.5 and $<$v1, numpy, scipy, joblib, matplotlib, numba $\geq$ v0.39.0, dill.

\section*{Software location:}
\begin{description}[noitemsep,topsep=0pt]
	\item[Name:] PyPI
	\item[Persistent identifier:] \url{https://pypi.org/project/coniii/}
	\item[Licence:] MIT License
	\item[Publisher:]  Edward D. Lee
	\item[Version published:] v1.1.1
	\item[Date published:] 12/12/18
\end{description}

{\bf Code repository}

\begin{description}[noitemsep,topsep=0pt]
	\item[Name:] GitHub
	\item[Persistent identifier:] \url{https://github.com/eltrompetero/coniii}
	\item[Licence:] MIT License
	\item[Version published:] v1.1.1
\end{description}

\begin{description}[noitemsep,topsep=0pt]
	\item[Name:] Zenodo
	\item[Persistent identifier:] \url{https://doi.org/10.5281/zenodo.2236632}
	\item[Licence:] MIT License
	\item[Version published:] v1.1.1
\end{description}

\section*{Language}
English

\section*{(3) Reuse potential}
To contribute either an algorithm for the inverse maxent problem or a sampling technique, we suggest following the template for the classes described in the base {\tt Solver} and {\tt Sampler} classes. New algorithms should be filed as a pull request to the \href{https://github.com/eltrompetero/coniii}{GitHub repository} along with an example solution that can be included in the usage guide Jupyter notebook and unit tests.

Documentation for the package is included as part of the GitHub repository and also hosted online at \url{https://eddielee.co/coniii/index.html}.

\section*{Acknowledgements}
We thank the anonymous reviewers for helpful feedback on both the manuscript and the accessibility of the software.

\section*{Funding statement}
EDL was supported by an NSF Graduate Fellowship under grant no. DGE-1650441. This research was supported in part by a congressional research grant provided by the Dirksen Congressional Center.

\section*{Competing interests}
The authors declare that they have no competing interests.

%\section*{Conclusion}
%Maxent provides a coherent framework for modeling a diverse body of biological and social systems.  The approach is structured to consider easily issues of model selection,  embodying a canonical approach for the scientific process: we start with the simplest models and the fewest constraints possible and then increase the complexity of the model until it makes good predictions. In principle, we could add as many constraints as would allow us to better fit the data, but complex models are more computationally expensive and often do not generalize well.  Maxent therefore stops before adding all possible constraints, leaving the model otherwise as structureless as possible.
%In this sense, the maxent approach is a principled method of statistical inference that is relevant beyond its roots in statistical physics.  We have built ConIII as an accessible software package in the hope that it will inspire those unfamiliar with maxent approaches to experiment and apply this technique to their own research questions.

\bibliography{refs}

\vspace{2cm}

\rule{\textwidth}{1pt}

%\section*{How to use}
%Each algorithm class instance must be instantiated with a few parameters that describe the particular maxent problem at hand. For example, all of the classes require a {\tt calc\_observables} keyword argument to be defined. This function will take a 2-dimensional sample (each row is a single state vector and each column is an entry in the state vector) from the ensemble and calculate the observables that correspond to the maxent model of interest. Some algorithms require other functions that are specific both to the problem and the algorithm. For example, the pseudolikelihood algorithm needs the likelihood of a single entry in the state vector ({\it i.e.} a single spin or a single vote) conditional on the rest of the state which depends on the maxent model. We have written helper methods for defining necessary inputs for all the available algorithms for the pairwise maxent problem. Examples of how to use these are given in the accompanying Jupyter notebook.
%
%Once an algorithm class is instantiated, it can be run by calling the {\tt solve} method and by providing the data in the corresponding format for the algorithm. Some algorithms like {\tt MCH} only require the correlations that are constrained in the maxent model. Others like {\tt MPF} require the full data set.
%
%To check if an algorithm has converged to reasonably close agreement with the data,
%
%Examples of this workflow for all the provided algorithms are given in the usage guide accompanying the Github package.

%\section*{The microcanonical ensemble and maximum entropy}
%The conventional textbook in statistical mechanics first introduces the concept of entropy as a way of counting the phase volume available to the system at a given energy
%\begin{align}
%	S(E) &= k_B \log\Omega(E)
%\end{align}
%where $k_B$ is Boltzmann's constant and $\Omega$ the number of states between $E$ and $E+\delta E$. Temperature is defined as
%\begin{align}
%	\frac{1}{T} &= \frac{\partial S}{\partial E}
%\end{align}
%
% begins with the concept of a small system coupled to a heat bath. In this limit, we can linearly expand thermodynamic quantities about the energy of the bath $E_{\rm bath} = E-E_s$
%\begin{align}
%	S_{\rm bath}(E-E_s) &\approx S_{\rm bath}(E) -E_s\left.\frac{\partial S}{\partial E_s}\right|_{E-E_s}\\
%		&= S_{\rm bath}(E) -\frac{E_s}{T}\\
%	e^{S_{\rm bath}(E-E_s)/k_B} &= e^{S_{\rm bath}(E)/k_B} e^{-E_s/k_BT}
%\intertext{Since the entropy is proportional to the density of states at a particular energy $E_s$, Eq ?? corresponds to}
%	p(E_s) &= e^{-E_s/k_BT}/Z
%\end{align}
%With partition function $Z$. In other words, the Gibbs measure takes exponential form. Eq ?? can be rewritten in terms of free energy
%\begin{align}
%	\log p(E_s) &= -E_s/k_BT -\log(Z)
%\intertext{Averaging both sides over all energy configurations and rearranging}
%	-k_B T\log Z &= \sum_s p(E_s) E_s +k_BT\sum_s p(E_s)\log p(E_s)
%\intertext{Remembering the fundamental postulate of statistical mechanics that all states are equally likely}
%	-k_BT\log Z &= \br{E_s} - TS\\
%	F &= \br{E_s} -TS
%\end{align}
%This equation tells us how the Helmholtz free energy is related to the internal energy and the entropy of the system.
%
%When we say that free energy is minimized for a system with Hamiltonian $E_s$ at equilibrium, we are equivalently saying that entropy is maximized. Entropy maximization is the crux of statistical physics models.

\end{document}

max_stars_repo_name: mjbommar/coniii

max_issues_repo_name: mjbommar/coniii

max_forks_repo_name: mjbommar/coniii


text: \chapter{Planning}

ext: tex

max_stars_repo_name: adamdboult/nodeHomePage

max_issues_repo_name: adamdboult/nodeHomePage

max_forks_repo_name: adamdboult/nodeHomePage


text:
\documentclass[class=NCU_thesis, crop=false]{standalone}
\begin{document}

\chapter{}

\section{}

\subsection{}

\end{document}


text:
\onecolumn
%\usepackage{MnSymbol}
%\usepackage{catchfilebetweentags}
%\input{macros}
%
%
%%commands with single parameter
%\newcommand{\p}[1]{\mathtt{process}_{#1}}
%\newcommand{\ta}[1]{\mathtt{take}_{#1}}
%\newcommand{\e}[1]{\mathtt{emit}_{#1}}
%\newcommand{\add}[1]{\mathtt{add}_{#1}}
%\newcommand{\f}[1]{\mathtt{fail}_{#1}}
%\newcommand{\startf}[1]{\mathtt{startFailure}_{#1}}
%\newcommand{\startid}[1]{\mathtt{startIdle}_{#1}}
%\newcommand{\id}[1]{\mathtt{idle}_{#1}}
%\newcommand{\cl}[1]{\mathtt{clock}_{#1}}
%\newcommand{\cltf}[1]{ \cl{to\f{#1}}}
%\newcommand{\ph}[1]{\mathtt{phase}_{#1}}
%
%%commands with two parameters (index, rate)
%\newcommand{\pr}[2]{\p{#1}(#2)}
%\newcommand{\tar}[2]{\ta{#1}(#2)}
%\newcommand{\er}[2]{\e{#1}(#2)}
%\newcommand{\addr}[2]{\add{#1}(#2)}
%
%\newcommand{\ra}[1]{r_{\add{#1}}}
%\newcommand{\rp}[1]{r_{\p{#1}}}
%\newcommand{\re}[1]{r_{\e{#1}}}
%\newcommand{\rt}[1]{r_{\ta{#1}}}
%\newcommand{\rf}[1]{r_{\mathtt{failure}_{#1}}}
%\newcommand{\rff}[2]{r_{\mathtt{failure}_{#1#2}}}
%\newcommand{\rr}[1]{r_{\mathtt{replay}_{#1}}}
%\newcommand{\reb}[1]{\bar{r}_{\e{#1}}}
%\newcommand{\rth}[1]{\hat{r}_{\ta{#1}}}
%
%\newcommand{\tph}[2]{t_{\ph{#1}}^{#2} }

\section{Appendix-extended results}\label{app:extended-exp-results}
\noindent
This appendix shows the some extended result of the verification tool.
We exercised on applying OSTIA-BASED formal verification to support continuous architecting of the ``focused-crawler topology''.
All the experiments were carried out with a random initial configuration for the queue occupancy which allowed us to detect possible executions leading to unbounded queues.
The three experiments were parameterized with different values of $ \alpha $ (time to process  one tuple).

%\begin{figure}[h]
%	\centering
%	\subfigure[{\footnotesize Starting configuration of the model: All the three bolt queues show an unbounded behavior. }]{\includegraphics[width=0.45\linewidth]{images/appendix/3bolts-1}}\label{vb1}
%	\hspace{0.5cm}
%	\subfigure[{\footnotesize WpDeserializer seems to better manage the incoming tuples flow, but it still unbounded, as the other bolts.}]{\includegraphics[width=0.45\linewidth]{images/appendix/3bolts-2}}\label{vb2}
%	%\hspace{0.5cm}
%	\subfigure[{\footnotesize This configuration allows WpDeserializer bolt to process the incoming tuple rate from the spout. Problems are still present for expander and articleExtraction }]{\includegraphics[width=0.45\linewidth]{images/appendix/3bolts-3}}\label{vb3}
%\end{figure}
\begin{figure}[h]
	\centering
	\includegraphics[width=0.7\linewidth]{images/appendix/3bolts-1}
	\caption{Starting configuration of the model: All the three bolt queues show an unbounded behavior.}
	\label{fig:3bolts-1}
\end{figure}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.7\linewidth]{images/appendix/3bolts-2}
	\caption{WpDeserializer seems to better manage the incoming tuples flow, but it still unbounded, as the other bolts.}
	\label{fig:3bolts-2}
\end{figure}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.7\linewidth]{images/appendix/3bolts-3}
	\caption{This configuration allows WpDeserializer bolt to process the incoming tuple rate from the spout. Problems are still present for expander and articleExtraction}
	\label{fig:3bolts-3}
\end{figure}

\clearpage
\section{Appendix-Formal model}\label{app:extended-formal-model}
\noindent
This appendix is an extract of the full CLTLoc model to provide the user with a clear overview of the formal background behind the verification in OSTIA.
Therefore, some part are missing, in particular, those related to the modeling of node failures.

A Storm topology formal definition is already given in Section~\ref{rs}.
%A Storm Topology is a directed graph $\mathbf{G} = \{ \mathbf{N}, \mathbf{E}\}$ such that the set of nodes $\mathbf{N} = \mathbf{S}\bigcup \mathbf{B}$ includes in the sets of spouts (\textbf{S}) and bolts (\textbf{B}), while the set of edges $\mathbf{E} = \{ Sub_{i,j} | i \in \mathbf{B}, j \in \{\mathbf{S}\bigcup \mathbf{B}\} \}$ defines how the nodes are connected each other expressing the subscription ($Sub_{i,j}$) of the bolts $i$ to the streams emitted by the spouts/bolts $j$.\\
%The main execution steps of each bolt (respectively spout) are: \texttt{take()} ( respectively \texttt{nextTuple()}), \texttt{execute()} (respectively \texttt{transform()}) and \texttt{emit()}. These steps are represented in a more generic way by the propositions:
The possible states and ``actions'' performed by the nodes are represented by the propositions:
\begin{itemize}
%\item $\ta{i}$ indicates that the \texttt{take()}/\texttt{nextTuple()} step is performed by the $i^{th}$ bolt/spout;
\item $\ta{i}$ indicates that the $i^{th}$ bolt is reading tuples from its queue;
%\item $\p{i}$ corresponds to the \texttt{execute()}/\texttt{transform()} step performed by the $i^{th}$ bolt/spout;
\item $\e{i}$ represents the $ i^{th} $ bolt/spout  emitting tuples;
\item $\p{i}$ corresponds to processing state of the $i^{th}$ bolt, that is, from the \texttt{take()} until the \texttt{emit()};
\end{itemize}

%\newcommand{\ori}{\mathtt{orig}}

\subsection*{Single node behaviour}
Let $\ori$ be a shorthand for $\neg\Y{\top}$.
The formula is true only in the origin.
The sequence of states of bolts is defined by the following properties, according to Fig.~\ref{figure-fsa} of Section~\ref{ver}:
\begin{align*}
 \bigwedge_{
i \in \mathbf{B} }
&(\ta{i} \Rightarrow
\p{i} \wedge \mathbf{X}(\lnot \ta{i} \U (\e{i} \lor \f{i})) \land \lnot \e{i}
\land \mathbf{Y}(\lnot \p{i} \, \Snc \, (\e{i}\lor \ori \lor \f{i}))
) \\
%
 \bigwedge_{
i \in \mathbf{B} }
&(\p{i} \Rightarrow
\p{i} \, \Snc \, ( \ta{i} \lor (\ori \land \p{i})) \land \p{i} \, \U \, (\e{i} \lor \f{i}) \land \lnot \f{i}) \\
%
 \bigwedge_{
i \in \mathbf{B} }
&(\e{i} \Rightarrow \p{i}
\wedge (\ori \lor \mathbf{Y}(\lnot \e{i} \, \Snc \, \ta{i})
\land \mathbf{X}(\lnot \p{i} \, \U \, \ta{i}) \\
%
\bigwedge_{
	\begin{subarray}{c}
	\,i \in \mathbf{B}
	\end{subarray}
}
&(\startid{i} \iff ( (\e{i} \land \mathbf{X}(q_i>0)) \lor (\lnot \f{i} \land \lnot \p{i} \land q_i>0 \land \lnot \mathbf{Y}(q>0 \lor \e{i}) ) ))\\
%
%
\bigwedge_{
	i \in \mathbf{B} }
&\left(
\id{i} \Rightarrow
\id{i} \, \Snc \, \startid{j} \land
\id{i} \, \U \, (\ta{i} \lor \f{i}) \land \lnot \f{i}
\right)  \\
%
%
\bigwedge_{
i \in \mathbf{S}} \,
&(\G{\F{\e{i}}}) \\
%
%
 \bigwedge_{
\begin{subarray}{c}
i \in \mathbf{B}
\end{subarray}} \,
&(\G{\F{\ta{i}}})
%
\end{align*}

%%%\subsubsection*{Node Failure}
%%%
%%%
%%%\begin{align*}
%%%%
%%% &\bigwedge_{
%%%\begin{subarray}{c}
%%%\,j \in B
%%%\end{subarray}
%%%}
%%%(\startf{j} \iff (\lnot \mathbf{Y}(\f{j}) \land \f{j}))\\
%%%%
%%%%
%%%%
%%%&\bigwedge_{
%%%\begin{subarray}{c}
%%%\,j \in B
%%%\end{subarray}
%%%}
%%%(\startf{j} \Rightarrow
%%%(
%%%\bigwedge_{
%%%	\begin{subarray}{c}
%%%	\,i \in B \\
%%%	i \neq j
%%%	\end{subarray}
%%%}
%%%(\lnot \startf{i})
%%%))
%%%%
%%%\end{align*}

\subsection*{Single queue behaviour}

Each bolt has a receive queue where the incoming tuples are collected before being read and processed by the node.
The queues have infinite size and the level of occupation of each $j^{th}$ queue is described by the variable $q_j$.\\
In order to express the connections among nodes in the topology, we defined the set of propositions $Sub(j,i)$ where $i \in \{ \mathbf{S} \bigcup \mathbf{B}\}$ and $j \in \mathbf{B}$. $Sub(j,i)$ means that "bolt $j$ submits to the stream emitted by spout/bolt $i$". This allows to establish a direct connection between the emission of a tuple by a node and the arrival of the tuple in the receive queues of the submitting nodes.\\
We also define $In_j=\{i_0,i_1, \dots, i_n\ | Sub(j,i) \}$ where $j\in \mathbf{B}$ and $i \in \{ \mathbf{S}\bigcup\mathbf{B}\}$, as the set of all the spouts/bolts whose streams are subscribed by the bolt $j$ (i.e. the nodes whose tuple are sent to the input queue of the bolt $j$).\\
The proposition $\add{j}$ is used to indicate that some of the nodes belonging to $In(j)$ are emitting tuples, that is, those tuples are added to the receive queue of bolt $j$.
\begin{align*}
 \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
} ( \add{j} \iff
%\lnot \f{i} \land FM -> not fail eliminato (se il bolt è fail la coda può riempirsi comunque)
\bigvee_{
\begin{subarray}{c}
\,i \in \{S \bigcup B \}\\
Sub(j,i)
\end{subarray}
 } \e{i})
\end{align*}
\subsubsection*{Quantites and Rates}
To represent the quantities of tuples that are taken/processed/emitted by each node in the topology, we introduce the concept of rate defined as "the quantity of exchanged tuples  in the current time unit".
Therefore, to express rates we use variables representing quantities and to measure time elapsing we exploit clocks of CLTLoc, as later explained in this section.

Following quantities define relevant metrics related to the message exchange process implemented by message brokers:
\begin{itemize}
\item $\re{j}$: number of tuples emitted in the current time unit by bolt/spout $j$.
\item $\rt{j}$: number of tuples taken in the current time unit by spout $j$.
\item $\ra{j}$: number of tuples that are added to the receive queue of spout $j$ in the current time unit, consisting in the sum of all the tuples currently emitted by the spouts/bolts belonging to $In(j)$.
\end{itemize}
%
Therefore, the evolution of each queue occupation, given the occurrence of adding and/or receiving event, is determined by the following formualae:
%
\begin{align*}
\bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
} & q_j \geq  0\\ % FM eliminato \wedge r \geq 0
%
 \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
} &( \add{j}  \wedge \lnot \ta{j}  \land \lnot \startf{j}
\Rightarrow (\aX q_j = q_j + \ra{j} ))\\ % FM eliminato \wedge r \geq 0
%
%
 \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
} &( \ta{j}{}  \Rightarrow (\aX q_j = q_j + \ra{j} - \rp{j} )) \\% \land q_j + \ra{j} \geq \rp{j})\\
%
%
%*** \bigwedge_{
%\begin{subarray}{c}
%\,j \in B
%\end{subarray}
%} &( \lnot \add{j}  \wedge \lnot \ta{j}{} \Rightarrow (\aX q_j = q_j))\\
%
%
% \bigwedge_{
%\begin{subarray}{c}
%\,j \in B
%\end{subarray}
%} &(\lnot \add{j}  \wedge \ta{j}{} \Rightarrow (\aX q_j = q_j - \rt{j} ) \land q_j \geq \rt{j})\\
%
%
 \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
} &( \startf{j}
\Rightarrow \mathbf{X}( q_j = 0 ))
\end{align*}
%
Following formulae define the necessary conditions ruling the occurrence of a variation of queue occupancy, i.e., if the queue occupancy increases then an $\add{}$ event occurs; if the queue occupancy decreases then a $\ta{j}$ action or a failure  $\startf{j}$ occurs.
%
\begin{align*}
 \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
} &\left( (\aX q_j > q_j)
\Rightarrow
\add{j}\right)\\
%
%
 \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
} &\left( (\aX q_j < q_j)
\Rightarrow
\ta{j} \lor \startf{j} \right)
\end{align*}
%
The following properties define the constraints applying to the tuple quantities of all bolts, distinguishing, when needed, between final ("leaf" node of the topology graph) and non-final bolt. $\mathit{final}(i)$ is equivalent to $\lnot \exists j: Sub(j,i)$:
%
%
\begin{align*}
& \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
} ( \ra{j} \geq 0 \land \rp{j} \geq 0 \land \re{j} \geq 0 ) \\
%
%
& \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
} ( \ra{j} = \sum_{
\begin{subarray}{c}
\,i \in \{S \bigcup B \}\\
Sub(j,i) \\
\e{i}
\end{subarray}
 } \re{i})
%
%
\end{align*}
%
%
\begin{align*}
& \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
} ( \ra{j} > 0 \iff \add{j}) \\
%
%
%& \bigwedge_{
%\begin{subarray}{c}
%\,j \in B
%\end{subarray}
%} ( \ra{j} =0 \iff \lnot \add{j}) \\
%
%
%
%& \bigwedge_{
%\begin{subarray}{c}
%\,j \in B
%\end{subarray}
%} ( \rt{j} > 0 \iff \ta{j}) \\
%%
%%
& \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
} ( \p{j} \Rightarrow \rp{j} > 0 ) \\
%
%
%
%& \bigwedge_{
%\begin{subarray}{c}
%\,j \in B
%\end{subarray}
%} ( \ta{j} \Rightarrow \rp{j}=\rt{j} ) \\
%%
%%
%%
& \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
} ( \rp{j} > 0 \Rightarrow (\p{j} \land \aX \rp{j}=\rp{j}) \, \U \, (\e{j}\lor \f{j})) \\
%
%
%
& \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
} ( \rp{j} = 0 \iff (\ori \lor (\neg\p{j} \land  \mathbf{Y}(\neg\p{j} \lor \e{j}) )) \\
%
%
%ELIMINATO DOPO INSERIMENTO IFF
%%& \bigwedge_{
%%\begin{subarray}{c}
%%\,j \in B
%%\end{subarray}
%%} ( ( (\startf{j} \lor \e{j}) \land \neg\X\ta{j}) \Rightarrow \X(\rp{j}=0) ) \\
%%%
%
%
& \bigwedge_{
\begin{subarray}{c}
\,j \in B\\
\lnot \mathit{final}(j)
\end{subarray}
} ( \e{j} \Rightarrow \re{j} = \sigma_j \rp{j} + d_j ) \land (\lnot \e{j} \Rightarrow \re{j} = 0 )\\
%
%
& \bigwedge_{
\begin{subarray}{c}
\,j \in B\\
\mathit{final}(j)
\end{subarray}
} ( \re{j} = 0)
%
%
%& \bigwedge_{
%\begin{subarray}{c}
%\,j \in B
%\end{subarray}
%} ( \re{j} =0 \iff \lnot \e{j})
%
\end{align*}
Each bolt has a maximum take capability $\rth{j}$ that limits the number of tuples that can be taken from the receive queue at any moment. That is, if the elements in the queue are greater than or equal to $\rth{j}$ and a $take$ is performed by the bolt, then exactly $\rth{j}$ tuples are taken from the queue. Otherwise, all the tuples are taken.
\begin{align*}
%
%
& \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
}((\ta{j} \land \rth{j} \geq q_j + \ra{j}) \Rightarrow (\rp{j} = q_j + \ra{j})) \\
%
%
%
%
& \bigwedge_{
\begin{subarray}{c}
\,j \in B
\end{subarray}
}((\ta{j} \land \rth{j} < q_j + \ra{j}) \Rightarrow (\rp{j} = \rth{j}))
%
%
\end{align*}

Since the spouts are modeled as node which only emit tuples into the topology, we define the amount of emitted tuples $\re{j}$ as the "nominal" emitting quantity $\reb{j}$ (i.e. the quantity of tuples that the spout would emit if no failure happens) plus the additional quantity $\rf{j}$ expressing the number of re-emitted tuples due to failures in the topology.
\begin{align*}
%
& \bigwedge_{
\begin{subarray}{c}
\,j \in S\\
\end{subarray}
} ( \re{j} = \reb{j} + \rr{j}) \\
%
%
& \bigwedge_{
\begin{subarray}{c}
\,j \in S\\
\end{subarray}
} ( \re{j} > 0 \iff \e{j}) \\
%
%
\end{align*}

%For final bolts the last formula is replaced with: $\re{j} = 0$

%\newpage

%%%%%%\subsubsection*{Failure Propagation}
%%%%%%In our model, whenever a node fails, the tuples being processed by the the node, together with the tuples in its receive queue, are considered as failed (not fully processed by the topology). According to the reliable implementation of Storm, the spout tuples that generated them must be resubmitted to the topology.\\
%%%%%%Since we do not keep track of the single tuples, but we only consider quantities of tuples throughout the topology, given an arbitrary amount of failed tuples, we can estimate the amount of spout tuples that have to be re-emitted by the connected spouts.
%%%%%%
%%%%%%In order to express this relationship between the failing tuples in a specific (failing) node and the new tuples having to be re-emitted, we introduce the concept of \textit{impact} of the node failure with respect to another (connected) node.
%%%%%%
%%%%%%Such impact can be precomputed given the topology and we define it as follows:
%%%%%%
%%%%%%$Imp(j,i)$ ("impact of node $j$ failure on node $i$") is the coefficient expressing the ratio $\frac{tuples\_to\_be\_replayed(i)}{failed\_tuples(j)}$ where $j \in B$ is the failing bolt and $i \in \{S \bigcup B\}$ is another node in the topology.\\
%%%%%%If exist a path $path(j,i)=\{p_0, \dots, p_n | n>0, p_0=j, p_n=i\}$ connecting the two nodes such that $\forall k \in [0,n-1] Sub(p_k,p_{k+1})$, then a failure of node $j$ has an impact on node $i$ and $Imp(j,i)>0$. If such a path does not exist, $Imp(j,i)=0$.\\
%%%%%%In order to define how to calculate $Imp(j,i)$ over a generic path we first show how to obtain its value for two basic cases:
%%%%%%%
%%%%%%\begin{itemize}
%%%%%%\item if two nodes $j$ and $i$ are directly connected:
%%%%%%\begin{align*}
%%%%%%&Sub(j,i) \Rightarrow
%%%%%%(Imp(j,i)= \frac{r_{out_i}}
%%%%%%{\sum_{
%%%%%%\begin{subarray}
%%%%%%\,k \in \{S \bigcup B\} \\
%%%%%%Sub(j,k)
%%%%%%\end{subarray}
%%%%%%} r_{out_k}})
%%%%%%%
%%%%%%\end{align*}
%%%%%%\ExecuteMetaData[storm-topology-example.tex]{failure-1}\\
%%%%%%
%%%%%%\item if the nodes $j$ and $i$ are connected by path passing through another node $h$:
%%%%%%\begin{align*}
%%%%%%%%
%%%%%%%& (Sub(j,h) \land Sub(h,i)) \Rightarrow
%%%%%%%\left( Imp(j,i)= Imp(j,h) \cdot \frac{1}{\sigma_h} \cdot
%%%%%%%\frac{r_{out_i}}
%%%%%%%{\sum_{
%%%%%%%\begin{subarray}
%%%%%%%\,k \in \{S \bigcup B\} \\
%%%%%%%Sub(j,k)
%%%%%%%\end{subarray}
%%%%%%%} r_{out_k}} \right)\\
%%%%%%%%
%%%%%%%
%%%%%%& (Sub(j,h) \land Sub(h,i)) \Rightarrow
%%%%%%\left( Imp(j,i)= Imp(j,h) \cdot \frac{1}{\sigma_h} \cdot
%%%%%%Imp(h,i)
%%%%%% \right)
%%%%%%\end{align*}
%%%%%%\ExecuteMetaData[storm-topology-example.tex]{failure-2}\\
%%%%%%\end{itemize}
%%%%%%
%%%%%%In general, if there is a path $path(j,i)=\{p_0, \dots, p_n | n>1, p_0=j, p_n=i\}$ defined as above:
%%%%%%\begin{align*}
%%%%%%Imp(j,i)=Imp(p_0,p_1)\cdot \prod_{k=1}^{n-1}\frac{1}{\sigma_k}\cdot Imp(k,k+1)
%%%%%%\end{align*}
%%%%%%
%%%%%%Once this coefficient is calculated for all the couples of $(bolt,spout)$ in the topology, it allows to determine the number of tuples to be re-emitted by each spout after a bolt failure by simply multiplying the number of failed tuples by the appropriate coefficient.\\
%%%%%%
%%%%%%%\subsubsection*{Failure Rates Behaviour - OLD}
%%%%%%%The impact of failure is therefore employed to calculate $\rf{j}$.
%%%%%%%The behaviour of $\rf{j}$ is defined as follows:
%%%%%%%\begin{align*}
%%%%%%%& \bigwedge_{
%%%%%%%\begin{subarray}{c}
%%%%%%%i \in S\\
%%%%%%%j \in B\\
%%%%%%%Imp(j,i)>0
%%%%%%%\end{subarray}
%%%%%%%} ((\f{j} \land \lnot \e{i}) \Rightarrow (\aX \rf{i} = \rf{i} + (q_j + \ra{j})\cdot Imp(j,i))))\\
%%%%%%%%
%%%%%%%& \bigwedge_{
%%%%%%%\begin{subarray}{c}
%%%%%%%i \in S\\
%%%%%%%j \in B\\
%%%%%%%Imp(j,i)>0
%%%%%%%\end{subarray}
%%%%%%%} ((\f{j} \land \e{i}) \Rightarrow (\aX \rf{i} = (q_j + \ra{j})\cdot Imp(j,i)))\\
%%%%%%%%
%%%%%%%%
%%%%%%%& \bigwedge_{
%%%%%%%\begin{subarray}{c}
%%%%%%%i \in S\\
%%%%%%%j \in B\\
%%%%%%%Imp(j,i)>0
%%%%%%%\end{subarray}
%%%%%%%} ((\lnot \f{j} \land \e{i}) \Rightarrow (\aX \rf{i} = 0))\\
%%%%%%%%
%%%%%%%& \bigwedge_{
%%%%%%%\begin{subarray}{c}
%%%%%%%i \in S\\
%%%%%%%j \in B\\
%%%%%%%Imp(j,i)>0
%%%%%%%\end{subarray}
%%%%%%%} \left( \rf{i}>0 \Rightarrow
%%%%%%%(\rf{i}>0 \, \Snc
%%%%%%%\bigvee_{
%%%%%%%\begin{subarray}{c}
%%%%%%%j \in B\\
%%%%%%%Imp(j,i)>0
%%%%%%%\end{subarray}
%%%%%%%} \f{j})
%%%%%%%\land (\rf{i}>0 \, \U \, \e{i})\right)\\
%%%%%%%%
%%%%%%%\end{align*}
%%%%%%
%%%%%%\subsubsection*{Failure Rates Behaviour - A (Single bolt formulae)}
%%%%%%%For each bolt we have:
%%%%%%%\begin{align*}
%%%%%%%%
%%%%%%%& \bigwedge_{
%%%%%%%\begin{subarray}{c}
%%%%%%%\,j \in S\\
%%%%%%%\end{subarray}
%%%%%%%} ( \re{j} = \reb{j} + \rr{j}) \\
%%%%%%%%
%%%%%%%\end{align*}
%%%%%%$\rr{j}$ is the quantity of tuples that need to be replayed due to a failure:
%%%%%%\begin{align*}
%%%%%%%
%%%%%%& \bigwedge_{
%%%%%%\begin{subarray}{c}
%%%%%%i \in S\\
%%%%%%\end{subarray}
%%%%%%}
%%%%%%(\rr{i} = \sum_{\begin{subarray}{c}
%%%%%%j \in B\\
%%%%%%Imp(j,i)>0
%%%%%%\end{subarray}} \rff{j}{i} \cdot Imp(j,i) )
%%%%%%%
%%%%%%\end{align*}
%%%%%%The impact of failure is therefore employed to calculate $\rff{j}{i}$ for each failing bolt.
%%%%%%The behaviour of $\rff{j}{i}$ is defined as follows:
%%%%%%\begin{align*}
%%%%%%& \bigwedge_{
%%%%%%\begin{subarray}{c}
%%%%%%i \in S\\
%%%%%%j \in B\\
%%%%%%Imp(j,i)>0
%%%%%%\end{subarray}
%%%%%%} ((\startf{j} \land \lnot \e{i}) \Rightarrow (\aX \rff{j}{i} = \rff{j}{i} + (q_j + \rp{j} + \ra{j}))))\\ %ELIMINATO \cdot Imp(j,i) (ANCHE SOTTO)
%%%%%%%
%%%%%%& \bigwedge_{
%%%%%%\begin{subarray}{c}
%%%%%%i \in S\\
%%%%%%j \in B\\
%%%%%%Imp(j,i)>0
%%%%%%\end{subarray}
%%%%%%} ((\startf{j} \land \e{i}) \Rightarrow (\aX \rff{j}{i} = q_j + \rp{j} + \ra{j}))\\
%%%%%%%
%%%%%%%
%%%%%%& \bigwedge_{
%%%%%%\begin{subarray}{c}
%%%%%%i \in S\\
%%%%%%j \in B\\
%%%%%%Imp(j,i)>0
%%%%%%\end{subarray}
%%%%%%} ((\lnot \startf{j} \land \e{i}) \Rightarrow (\aX \rff{j}{i} = 0))\\
%%%%%%%
%%%%%%%
%%%%%%%
%%%%%%& \bigwedge_{
%%%%%%\begin{subarray}{c}
%%%%%%i \in S\\
%%%%%%j \in B\\
%%%%%%Imp(j,i)>0
%%%%%%\end{subarray}
%%%%%%} ((\lnot \startf{j} \land \lnot \e{i}) \Rightarrow (\aX \rff{j}{i} = \rff{j}{i}))\\
%%%%%%%
%%%%%%%
%%%%%%%& \bigwedge_{
%%%%%%%\begin{subarray}{c}
%%%%%%%i \in S\\
%%%%%%%j \in B\\
%%%%%%%Imp(j,i)>0
%%%%%%%\end{subarray}
%%%%%%%} \left( \rf{i}>0 \Rightarrow
%%%%%%%(\rf{i}>0 \, \Snc
%%%%%%%\bigvee_{
%%%%%%%\begin{subarray}{c}
%%%%%%%j \in B\\
%%%%%%%Imp(j,i)>0
%%%%%%%\end{subarray}
%%%%%%%} \f{j})
%%%%%%%\land (\rf{i}>0 \, \U \, \e{i})\right)\\
%%%%%%%%
%%%%%%\end{align*}

%\subsubsection*{Failure Rates Behaviour - B (Spout perspective)}
%
%
%The impact of failure is therefore employed to calculate $\rf{j}$.
%The behaviour of $\rf{j}$ is defined as follows:
%\begin{align*}
%&
%\bigwedge_{
%\begin{subarray}{c}
%\mathtt{FAILED} \subseteq B
%\end{subarray}
%}\left(
% \bigwedge_{
%\begin{subarray}{c}
%i \in S\\
%\end{subarray}
%}
% \left(
%(
%\bigwedge_{
%\begin{subarray}{c}
%j \in \mathtt{FAILED}
%\end{subarray}
%}\startf{j}
%\land \lnot \e{i})
%\Rightarrow
%(\aX \rf{i} = \rf{i} + \sum_{j \in \mathtt{FAILED}}(x_j\cdot Imp(j,i))))
%\right)
%\right)\\
%%
%%
%%
%%
%&
%\bigwedge_{
%\begin{subarray}{c}
%\mathtt{FAILED} \subseteq B
%\end{subarray}
%}\left(
% \bigwedge_{
%\begin{subarray}{c}
%i \in S\\
%\end{subarray}
%}
% \left(
%(
%\bigwedge_{
%\begin{subarray}{c}
%j \in \mathtt{FAILED}
%\end{subarray}
%}\startf{j}
%\land \e{i})
%\Rightarrow
%(\aX \rf{i} = \sum_{j \in \mathtt{FAILED}}(x_j\cdot Imp(j,i))))
%\right)
%\right)\\
%%
%%
%%
%& \bigwedge_{
%\begin{subarray}{c}
%i \in S
%\end{subarray}
%} \left((
%\bigwedge_{
%\begin{subarray}{c}
%j \in B\\
%Imp(j,i)>0
%\end{subarray}
%}
%(\lnot \startf{j}) \land \e{i}) \Rightarrow (\aX \rf{i} = 0)
%\right)\\
%%
%%
%%
%& \bigwedge_{
%\begin{subarray}{c}
%i \in S
%\end{subarray}
%} \left((
%\bigwedge_{
%\begin{subarray}{c}
%j \in B\\
%Imp(j,i)>0
%\end{subarray}
%}
%(\lnot \startf{j}) \land \lnot \e{i}) \Rightarrow (\aX \rf{i} =  \rf{i})
%\right)\\
%%
%%
%\end{align*}
%
%
%\subsection*{Average rates}
%Given a an average rate \texttt{M} and a period $\tau$ , we want to ensure that every $\tau$ time intervals, the average value of the rate $r_i$ registered so far is exactly \texttt{M}.
%
%\begin{align*}
% \bigwedge_{
%\begin{subarray}{c}
%i > 0
%\end{subarray}
%}
%(r_{i \tau} = \mathtt{M}\tau - \sum_{j=(i-1)\tau+1}^{i\tau - 1}r_j)
%\end{align*}

%OLD QUEUE FORMULAE
%\newpage
%\input{queue_behaviour_old}

\subsection*{Clocks Formulae}
In order to represent the duration of the various processing phases of each bolt we introduce different clocks:
\begin{itemize}
	\item \textbf{$ \tph{j}{0} $} and \textbf{$ \tph{j}{1} $} measure the duration of the $ \p{j} $, $ \id{j} $ and $ \f{j} $ phases for each bolt $ j $ and the time elapsed between one $ \e{j} $ and the next one for each spout $ i $.
	\item \textbf{$ \cltf{j} $} measures the \textit{time to failure}, i.e. the time elapsing between the end of a failure and the beginning of the next one for each bolt $ j $.
\end{itemize}
%
Following formulae impose that $\tph{j}{0}$ and $\tph{j}{1}$ are reset alternatively and never at the same time.
\begin{align*}
%
& \bigwedge_{
	\begin{subarray}{c}
	j \in \{ S \bigcup  B \}\\
	\end{subarray}
	}
	(\tph{j}{0} =  0 \Rightarrow \X{(\tph{j}{1} = 0) \R (\tph{j}{0}>0)} \land (\tph{j}{1} > 0) \\
	%
	& \bigwedge_{
	\begin{subarray}{c}
	j \in \{ S \bigcup  B \}\\
	\end{subarray}
	}
	(\tph{j}{1} =  0 \Rightarrow \X{(\tph{j}{0} = 0) \R (\tph{j}{1}>0)})
		%
		\end{align*}

Each clock $ \tph{j}{0} $ will be initially set to 0.
In the following formula, we will use the shorthand $ \tph{j}{} \sim c$ to indicate the formula:

\begin{align*}
%
&  (\tph{j}{0} >  0 \land  (\tph{j}{1} \lor (\tph{j}{1} > 0) \Snc (\tph{j}{0})) \Rightarrow \tph{j}{0} \sim c )\\
%
& \land \\
%
& (\tph{j}{1} >  0 \land  (\tph{j}{0} \lor (\tph{j}{0} > 0) \Snc (\tph{j}{1})) \Rightarrow \tph{j}{1} \sim c )\\
%
\end{align*}
%
which, intuitively, has the ability of choosing the last clock which was reset before the current position.
The last clock reset keeps the current time delay since the last event associated with it and, therefore, it is the clock to employ for measuring the next deadline.

Reset conditions for $\tph{j}{0}$ and $\tph{j}{1}$ are defined by next formulae which bind the occurrence of specific combination of events with the reset itself.

\begin{itemize}
 \item \textbf{ $ \tph{j}{} $ for bolts} - start of processing, failure or idle phase.
\begin{align*}
%
& \bigwedge_{
	\begin{subarray}{c}
	j \in B\\
	\end{subarray}
}
\left(
%\begin{array}{l}
(\tph{j}{} =  0) \iff %\\
\left(\ori \lor \ta{j} \lor %\\
(\f{j} \land \lnot \Y (\f{j}) )  \lor %\\
(\id{j} \land \lnot \Y (\id{j}) ) \right)
%\end{array}
\right)\\
%
\end{align*}
 \item \textbf{$ \tph{j}{} $ for spouts } - clock resets every time the corresponding spout emits
\begin{align*}
%
& \bigwedge_{
	\begin{subarray}{c}
	i \in S\\
	\end{subarray}
}
((\tph{j}{} =  0) \iff
\e{i})\\
%
\end{align*}
 \item \textbf{$ \cltf{j} $ (bolts)} - time-to-failure clock resets every time a failure ends.
\begin{align*}
%
& \bigwedge_{
	\begin{subarray}{c}
	j \in B\\
	\end{subarray}
}
((\cltf{j} =  0) \iff
 (\lnot \f{j} \land \lnot \Y (\lnot \f{j})))\\
%
\end{align*}
\end{itemize}

Following formulae determine the duration of permanence in one state.

Processing duration: $\alpha_j$ is the time required by bolt $j$ to process one tuple and $\epsilon$ is a fixed defined tolerance.
\begin{align*}
%
%
&\bigwedge_{
	\begin{subarray}{c}
	j \in B\\
	\end{subarray}
}
\left( \begin{array}{l}
(\p{j} \Rightarrow \\
(\p{j} \land \lnot \e{j}) \U (( \tph{j}{} \geq \alpha_{j} - \epsilon) \land ( \tph{j}{} \leq \alpha_{j} + \epsilon) \land (\f{j} \lor \e{j}) )
\end{array} \right)\\
\end{align*}

Failure duration:

\begin{align*}
%
\bigwedge_{
	\begin{subarray}{c}
	j \in B\\
	\end{subarray}
}
\left( \begin{array}{l}
\f{j} \Rightarrow \\
(\f{j}  \U (( \tph{j}{} \geq \mathtt{MIN\_REBOOT\_TIME_j}) \land ( \tph{j}{} < \mathtt{MAX\_REBOOT\_TIME_j}) \land \lnot \f{j})
\end{array} \right)\\
%
%
\end{align*}

Idle duration:

\begin{align*}
%
\bigwedge_{
	\begin{subarray}{c}
	j \in B\\
	\end{subarray}
}
\left( \begin{array}{l}
\id{j}  \Rightarrow \\
(\id{j}  \U (( \tph{j}{} \geq \mathtt{MIN\_IDLE\_TIME_j}) \land ( \tph{j}{} < \mathtt{MAX\_IDLE\_TIME_j}) \land \lnot \id{j})
\end{array} \right)\\
%
%
\end{align*}

Minimum time to failure:

\begin{align*}
%
\bigwedge_{
	\begin{subarray}{c}
	j \in B\\
	\end{subarray}
}
\left( \begin{array}{l}
\lnot \f{j} \Rightarrow \\
(\lnot \f{j}  \U ( \cltf{j} > \mathtt{MIN\_TTF_j}))
\end{array} \right)\\
%
%
\end{align*}

Spout emitting intervals: to define emitting rates we here define the ratio between the number of emitted tuples in an time interval and the duration of  time interval.
Clock $\tph{j}{}$ measure the duration of the interval where a certain amount of tuples is emitted and the number of emitted tuples is expressed by $\mathtt{proc\_times[r][\cdot]}$.
The amount of tuples $r$ that can be emitted is discretized and each possible amount in $(0, \reh{j}]$ is associated with a proper delay.
This allows us to impose a specific emitting rate for each spout.
To set up a proper tolerance in measuring time duration, we employ two time bounds $\mathtt{proc\_times[r][0]}$ and $\mathtt{proc\_times[r][1]}$.
Hence, the ratio imposed by the formula is between $\frac{r}{\mathtt{proc\_times[r][1]}}$ and $\frac{r}{\mathtt{proc\_times[r][0]}}$, for each $r \in (0, \reh{j}]$.

\begin{align*}
%
&\bigwedge_{
	\begin{subarray}{c}
	j \in \mathbf{S}\\
	\end{subarray}
}
\left( \begin{array}{l}
(\lnot \e{j} \mathbf{S} \ori) \Rightarrow \\
(\lnot \e{j} \U
(( \tph{j}{} \leq \mathtt{proc\_times[1][0]}) \land  \e{j}) )
\end{array} \right)\\
%
%
&\bigwedge_{
	\begin{subarray}{c}
	r \in (0, \reh{j}]\\
	\end{subarray}
}
\left(
\bigwedge_{
	\begin{subarray}{c}
	j \in \mathbf{S}\\
	\end{subarray}
}
\left( \begin{array}{l}
(\e{j} \land \reb{j} = r) \Rightarrow \\
(\mathbf{X}(\lnot \e{j} \U
(( \tph{j}{} \geq \mathtt{proc\_times[r][0]}) \land ( \tph{j}{} < \mathtt{proc\_times[r][1]}) \land  \e{j}) ))
\end{array} \right) \right)\\
\end{align*}

max_stars_repo_path: [after-JSS] OSTIA/appendix.tex

max_stars_repo_name: maelstromdat/OSTIA

max_issues_repo_path: [after-JSS] OSTIA/appendix.tex

max_issues_repo_name: maelstromdat/OSTIA

max_forks_repo_path: [after-JSS] OSTIA/appendix.tex

max_forks_repo_name: maelstromdat/OSTIA


text:
\subsection*{Pythagorean identities}

\begin{align*}
  \sin^2 x + \cos^2 x &= 1\\
  1 + \tan^2 x &= \sec^2 x\\
  1 + \cot^2 x &= \csc^2 x
\end{align*}


text:
\section{Data Reduction}
% description of how reduced and calibrated  your data

The CCD observations are registered as astronomical images (FITS, Flexibble
Image Transport System). All the
analysis was done in IDL and ATV \cite{atv} - the source codes are
included in the appendix.

\subsection{CCD Dark Frame Calibration}

The first step of the data reduction is to remove the various instrumental
artifacts of the CCD, by calibrating  the science images. We create {\it dark
master frames} for every exposure times, \ie median combined frames of the CCD
closed, the high signal-to-noise dark. We subtract these frames from the
respective science frames with same exposure times and median combine these
last to form a deep exposure frame.

\subsection{Spectra Trace Calibration}

Spectra are arranged on the CCD in a manner that most efficiently makes
use of the detector size/area. Before these spectra can be extracted, their
exact positions on the CCD must be mapped. This is the process of {\it
aperture tracing}. To this, we calculated the  spectra's {\it trace} from two
methods:
\begin{enumerate}
 \item Method 1: From a standard point source;
\item Method 2: By integrating the 1 + 1/2 strips.
\end{enumerate}

The trace of the point source, \eg a bright
(standard) star in the sky, is taken by positioning this star in each of the
three slits.
With ATV \cite{atv}, it is possible to extract the {\it spectral flux
distribution} \footnote{Loading the image with \texttt{decomposed, device=0}
option and pressing \texttt{x} on it.}, $F_{\lambda}$,  as shown in the Fig.
\ref{f}.

\begin{figure}[htb]
\begin{center}
\includegraphics[scale=0.25]{plots/trace/point2.png}
\includegraphics[scale=0.25]{plots/trace/point3.png}
\includegraphics[scale=0.25]{plots/trace/point1.png}
\caption{Trace extraction from a point source  star in each of the three
slits
(from ATV).}
\label{f}
\end{center}
\end{figure}

 For the Moonshine, the flux is very high and using only the the point
source trace (method 1) was enough.
(see Figs. \ref{star-cal}).

\begin{figure}[htb]
\begin{center}
\includegraphics[scale=0.33]{plots/trace/moon.png}
\includegraphics[scale=0.33]{plots/trace/moonsky.png}
\includegraphics[scale=0.33]{plots/trace/earth.png}
\includegraphics[scale=0.33]{plots/trace/earthsky.png}
\caption{Extraction of the point source spectra for the bright (above) and dark
(below) side of the
waning crescent moon and its adjacent sky (from ATV).}
\label{star-cal}
\end{center}
\end{figure}

Only for the Earthshine, due
to their much lower signal (in the same order of the adjacent sky), we determine
a second set of data (method 2), aiming to extract more statistics from the
data. We proceeded integrating over the whole
Earthshine strip and adding half of the middle strip (which has the adjacent
sky on the other half).

\subsection{Adjacent Sky Subtraction}
Both the dark and bright side of the Moon are corrected for scattered
light in the telescope by subtracting the adjacent sky spectrum. The adjacent
skies were always extracted in the same way the their respective set of data. The
sky-subtracted spectra can been in the see Figs. \ref{comparing}.

\begin{figure}[htb]
\begin{center}
\includegraphics[scale=0.31]{plots/earth_point.png}
\includegraphics[scale=0.31]{plots/spectra1.png}
\includegraphics[scale=0.31]{plots/earth_strip.png}
\includegraphics[scale=0.31]{plots/spectra2.png}
\caption{Earthshine spectra from method 1(top) and method 2 (down), comparing
to the scaled Moonshine spectra, in function of CCD pixel counts, and
wavelength calibrated.}
\label{comparing}
\end{center}
\end{figure}

After these
reductions we were able to calculate the signal
to background for
the Earthshine (methods 1 and 2) and for the Moonshine spectra, obtaining
$S/N^{point}_{earth} = 1.8$,  $S/N^{strip}_{earth} = 3.4$ and
$S/N_{moon} = 148$, respectively.

\subsection{Obtaining the Reflectivity of Earth}

The reflectivity of Earth is the ratio of the Earthshine to the Moonlight
spectra, \ie by reducing the contribution of the extra
passage of the Sun through the Earth's atmosphere in the spectrum of the
Earthshine. Dividing  the dark side by the bright side (the earthshine by the
 moonshine) for methods 1 and 2, we obtain the Fig. \ref{comparing2}. We
confirm that they show a very good agreement, with  a very small enhancement
from the method 2. We shall use only this spectra on the following analysis and
refer to it as the Earthshine reflectivity.

\begin{figure}[htb]
\begin{center}
\includegraphics[scale=0.6]{plots/reflectivity0.png}
\caption{Earthshine reflectivity from method 1 and 2, scaled to each
other.}
\label{comparing2}
\end{center}
\end{figure}

\subsection{Absolute Spectra Calibration}
To obtain the absolute wavelengths of the Earthshine spectrum, we measure the
spectrum of a neon calibration lamp and compare it to the well-known wavelengths
of neon
transitions, see Fig. \ref{fig-neon}.

\begin{figure}[htb]
\begin{center}
\includegraphics[scale=0.2]{figs/neon.jpg}
\includegraphics[scale=0.45]{plots/neon1.png}
\caption{Wavelength transitions of neon: (left) from the literature
\cite{neon-wave}, (right)  experimentally obtained.}
\label{fig-neon}
\end{center}
\end{figure}

We find the {\it $\lambda$/pixels scale dispersion} and we reduce the
Earthshine spectra from it.This results on the Earthshine spectra
correlation to light wavelengths, \ie the reflectivity of the Earth.

\begin{figure}[htb]
\begin{center}
\includegraphics[scale=0.6]{plots/reflectivity2.png}
\caption{The Reflectivity of the Earth.}
\label{fig-es}
\end{center}
\end{figure}


text:
\section{Introduction}

% Motivation
Main-sequence dwarfs are the most common stars in the Milky Way, and their
ages could reveal the evolution of Galactic stellar populations and planetary
systems.
However, the ages of low-mass stars, particularly K and M dwarfs are difficult
to measure because their luminosities and temperatures evolve slowly on the
main sequence \citep[see][for a review of stellar ages]{soderblom2010}.
Fortunately, rotation-dating, or `gyrochronology’ provides a promising means
to measure precise ages for these cool dwarfs
\citep[\eg][]{schatzman1962, weber1967, kraft1967, skumanich1972, kawaler1988,
pinsonneault1989, barnes2003, barnes2007, mamajek2008, barnes2010, meibom2011,
epstein2014, meibom2015, vansaders2016, vansaders2018, claytor2020}.
% Gyro 101
Stars with significant convective envelopes ($\lesssim$ 1.3 M$_\odot$) have
strong magnetic fields, generated by a Solar-type $\alpha-\Omega$ dynamo.
These stars lose angular momentum over long timescales via magnetic braking
\citep[\eg][]{schatzman1962, weber1967, kraft1967, skumanich1972, kawaler1988,
pinsonneault1989}.
Stars arrive on the main sequence with a random distribution of rotation
periods, ranging from around 1 to 10 days \citep[\eg][]{rebull2018}, however
observations of open clusters show that stellar rotation periods
converge onto a unique sequence by around 500-700 million years
\citep[\eg][]{irwin2009, gallet2013}.
This sequence is often called the `slow rotator sequence'.
After stars reach the slow rotator sequence, the rotation period of a star can
be determined, to first order, by its photometric color and age
\citep[\eg][]{barnes2003, barnes2007, barnes2010, meibom2011, meibom2015}.

A well-calibrated gyrochronology model that captures the time and
mass-dependence of stellar angular momentum evolution could provide ages that
are precise to within 20\% for millions of Milky Way stars in the time-domain
era \citep{epstein2014, najita2016, angus2019, claytor2020}.
However, gyrochronology models are not yet fully calibrated, especially for K
and M dwarfs, as they are limited by a lack of low-mass stars with precisely
measured ages and rotation periods.
For example, the oldest K dwarfs in an open cluster with measured rotation
periods are those in Ruprecht 147 (2.7 Gyr) \citep{curtis2020}.
The oldest cluster M dwarfs with measured rotation periods are those in
Praesepe (600-700 Myr) \citep{douglas2017, rebull2017}.
Gyrochronology relations can therefore be considered untested for K dwarfs
older than 2.7 Gyr and M dwarfs older than 600-700 Myr.

Once stars have converged to the slow rotator sequence, their rotation periods
broadly follow a power law relation in time with an index of around 0.5
\citep{skumanich1972}.
However, new Kepler/K2 rotation periods of stars in open clusters and
asteroseismic stars, show that the exact value of this power law index varies
as a function of mass and age \citep{angus2015, vansaders2016, curtis2019,
spada2019}.
Two major phenomena in particular, `weakened magnetic braking' and `stalled
braking' are both responsible for slowing the rate of magnetic
braking.
These two phenomena affect stars of different masses and ages and are thought
to be generated by different physical processes.
`Weakened magnetic braking' was first observed in asteroseismic stars with
rotation rates that appeared too rapid for their age
\citep{angus2015}.
It affects the rotational evolution of G dwarfs older
than around Solar age (4-5 Gyrs) and is thought to be caused by a shutting down
of the magnetic dynamo when stars reach a Rossby number of around 2
\citep{vansaders2016, vansaders2018}.
\footnote{Rossby number is the ratio of rotation period to convective overturn
time.}
`Stalled braking' is a different behavior that has been observed in
the K dwarfs of at least three open clusters: NGC 752, NGC 6811, and Ruprecht
147 \citep{agueros2018, curtis2019, curtis2020}.
The K dwarfs in these clusters again appear to rotate too rapidly for their
age, however in this case the underlying physical cause is thought to be the
redistribution of angular momentum in the stellar interior \citep{spada2019}.

With the discovery of weakened and stalled braking, it is clear that stars of
different masses and ages spin-down at different rates.
Because of the non-linear and time-variable relation between rotation period
and age, the gyrochronology relations should not, if possible, be extrapolated
to old ages or low masses where little calibration data exist.
Instead, they should be actively calibrated for stars of all masses and ages.
Our goal in this paper is to provide a new fully-empirical gyrochronology
relation that captures the stalled braking behavior of old K and early M
dwarfs.
% However, our model will not capture the behavior of stars undergoing weakened
% braking because our calibration is not yet suitable for weakened braking
% stars.
% Similarly, we do not have any mid or late M dwarfs in our calibration sample,
% so our model is also unsuitable for those stars.
% We hope to extend our model to include weakened magnetic braking and mid/late
% M dwarfs in future.
% We also note that a fully-empirical model has certain limitations, and we
% describe the full list of caveats later in this paper.

% The rate at which a star loses angular momentum is thought to chiefly depend
% on its magnetic field strength.
% Lower-mass stars with deeper convection zones have stronger magnetic fields,
% thus larger Alfvén radii, and therefore experience greater angular momentum
% loss rates \citep[\eg][]{schatzman1962, kraft1967, parker1970, kawaler1988,
% charbonneau2010}.
% This effect dominates the evolutionary sequence of young stars and, for
% example in the $\sim$ 700 Myr Praesepe cluster, rotation rate increases with
% increasing mass (rotation {\it period decreases} with increasing mass).
% However, at ages greater than around 1 Gyr \citep{spada2019, curtis2019,
% angus2020}, internal angular momentum transport becomes important and starts
% to influence the evolution of surface rotation periods.

% The rotation periods of these calibration stars can often be measured with
% precise time-series photometry, and the Kepler spacecraft has played a leading
% role in providing light curves for stellar rotation studies
% \citep[\eg][]{borucki2010, meibom2011, mcquillan2013, mcquillan2014,
% howell2014, reinhold2014, garcia2014, aigrain2015, meibom2015, douglas2016,
% douglas2017, rebull2016, rebull2017, santos2019, reinhold2019, rebull2020,
% gordon2020, curtis2020, breton2021}.
% Magnetically active regions on the surfaces of stars create dark and bright
% surface features which produce periodic variability in the overall brightness
% of a rotating star.
% This variability is, at most, on the order of 1\% of the star's overall flux.
% The rotation periods of stars can often be measured from their light curves
% using signal-processing techniques.
% TESS, the Transiting Exoplanet Survey Satellite, is continuing the legacy of
% Kepler and will provide rotation periods for many more stars
% \citep{ricker2014}.

% For the purposes of calibrating gyrochronology, oen clusters provide good
% mass coverage for young stars: rotation periods have been measured for F to
% mid M dwarfs up to ages of around 700 Myr.
% In contrast, asteroseismic stars provide reasonable age coverage for hot
% stars: ages and photometric surface rotation periods have been measured for F,
% G and early K dwarfs up to ages of 10 Gyr.
% However, neither asteroseismology nor cluster analysis can provide rotation
% periods and ages for old, late K and M dwarfs.
% In addition, cluster and asteroseismic stars generally provide sparse coverage
% of the rotation period-effective temperature plane, and cannot reveal the
% detailed evolution of stellar rotation rates.
% As a result, most empirical gyrochronology relations are only reliable for G
% dwarfs up to Solar age, K dwarfs up to 2-3 Gyr, and early M dwarfs up to $<$ 1
% Gyr.

Historically, gyrochronology calibration samples have mainly been comprised of
stars in open clusters and stars with detectable acoustic pulsations, both of
which can be precisely dated with stellar evolution models.
However, to extend gyrochronology relations to older K and M dwarfs, it may be
necessary to use alternative dating methods.
Old open clusters are generally too distant for photometric rotation period
measurements of faint M dwarfs, and asteroseismic analyses are impacted by the
large magnetic activity signals and low-amplitudes of oscillations for
low-mass stars.
There is great promise, however, in using binaries as an alternative method
for calibrating gyrochronology.
For example, low-mass K and M dwarfs with a precisely dateable co-moving
companion (\eg\ a white dwarf, a subgiant, or a star with detectable acoustic
oscillations) could be good calibrators.
In addition, Galactic kinematics is a dating method with good potential for
gyrochronology calibration.
In \citet{angus2020} we explored the kinematic properties of Kepler field
stars and found that the velocity dispersion of stars increase with stellar
rotation period, as expected (because both quantities increase with age).
In \citet{lu2021} we used the velocity dispersions of Kepler field stars to
estimate their ages using an Age-Velocity dispersion Relation (AVR).
In this paper, we use those kinematic ages for around 30,000 Kepler field
stars to extend the gyrochronology relations into the old K and early M dwarf
regime.

The star forming molecular gas clouds observed in the Milky Way have a low
out-of-plane, or vertical, velocity \citep[\eg][]{stark1989, stark2005,
aumer2009, martig2014, aumer2016}.
In contrast, the vertical velocities of older stars are observed to be larger
in magnitude on average \citep{stromberg1946, wielen1977, nordstrom2004,
holmberg2007, holmberg2009, aumer2009, casagrande2011, ting2019, yu2018}.
There are two possible explanations for this observed increase in velocity
dispersion with age: either stars are born kinematically `cool' and their
orbits are heated over time via interactions with giant molecular clouds
\citep[see][for a review of secular evolution in the MW]{sellwood2014}, or
stars formed kinematically `hotter' in the past \citep[\eg][]{bird2013}.
Either way, the vertical velocity dispersions of thin disk stars are observed
to increase with stellar age.
This behavior is codified by Age-Velocity dispersion Relations (AVRs), which
typically express the relationship between age and velocity dispersion as a
power law: $\sigma_v \propto t^\beta$, with free parameter, $\beta$
\citep[\eg][]{holmberg2009, yu2018}.
These expressions can be used to infer the ages of groups of stars from their
velocity dispersions.

Kinematic ages have been used to explore the evolution of cool dwarfs for over
a decade.
For example, \citet{west2004, west2006} found that the fraction of
magnetically active M dwarfs decreases over time, by using the vertical
distances of stars from the Galactic mid-plane as an age proxy, and
\citet{west2008} used kinematic ages to calculate the expected activity
lifetime for M dwarfs of different spectral types.
\citet{faherty2009} used tangential velocities to infer the ages of M, L and T
dwarfs, and showed that dwarfs with lower surface gravities tended to be
kinematically younger, and \citet{kiman2019} used velocity dispersion as an
age proxy to explore the evolution of H$\alpha$ equivalent width (a magnetic
activity indicator), in M dwarfs.

AVRs are usually calibrated in Galactocentric velocity coordinates (\vx, \vy,
\vz\ or $UVW$), and these velocities can only be calculated with full 6D
positional and velocity information, however most \kepler\ rotators do not
have RV measurements\footnote{Although RVs for most will be released in \gaia\
DR3}.
In \citep{angus2020} we used velocity in the direction of Galactic latitude
(\vb) as a stand-in for \vz\ because, in the {\it Galactic} coordinate system,
velocities can be calculated from 3D positions and {\it 2D} proper motions.
The \kepler\ field lies at low Galactic latitude, so \vb\ is a close
approximation to \vz\ for Kepler stars.
Although \vb\ velocity dispersion does not equal \vz\ velocity dispersion, it
still increases monotonically over time and provides accurate age rankings for
\kepler\ stars.
In this paper however, we calculate or infer \vz\ for each Kepler star in our
sample.
For stars without RV measurements we infer {\it vertical} velocity, \vz, by
marginalizing over missing RVs using a hierarchical Bayesian model.
These vertical velocities were used by \citet{lu2021} to calculate velocity
dispersions of groups of stars with similar properties, which they then
converted into ages using an AVR \citep{yu2018}.
In this paper, we use the kinematic ages calculated by \citet{lu2021} to
calibrate a new gyrochronology relation.

This paper is laid out as follows: in section \ref{sec:data} we describe the
data used to calibrate our new gyrochronology model.
In section \ref{sec:velocities} we describe how we calculated the kinematic
ages of Kepler stars from their positions, proper motions and parallaxes (and
RVs where available).
In section \ref{sec:gyro_calibration} we describe the new model and how we fit
it to the data.
In section \ref{sec:results} we test our new model using cross validation and
test its precision and accuracy, and in section \ref{sec:discussion} we
discuss the successes and caveats of our new model.

max_stars_repo_name: RuthAngus/aviary

max_issues_repo_name: RuthAngus/aviary

max_forks_repo_name: RuthAngus/aviary