text:
% Preamble
\documentclass[12pt,a4paper]{article}
\usepackage{enumerate}
\usepackage{setspace}
\usepackage{authblk}
\usepackage{graphicx}
%\usepackage[nomarkers, nolists, figuresonly]{endfloat}
\usepackage{pdflscape}
\usepackage{mathtools}
\usepackage[osf]{mathpazo}
\usepackage{lineno}
\usepackage{ms}
\usepackage{hyperref}
\usepackage[round]{natbib}
\usepackage{setspace}

\setcounter{secnumdepth}{0}
\raggedright
\pagenumbering{arabic}
\linenumbers

% First order headings upper case bold
\usepackage{titlesec}
\titleformat*{\section}{\small\bfseries\uppercase}

% Second order headings normal case italics
\titleformat*{\subsection}{\small\itshape}

% Third order, italics, paragraph style
\titleformat*{\paragraph}{\small\itshape}

% Title page information

\title{Time for a rethink: time sub-sampling methods in disparity-through-time analyses}

\author{
	Thomas Guillerme$^{1}$ and Natalie Cooper$^{2*}$
}

\date{}

\affiliation{\noindent{\footnotesize
	$^1$School of Biological Sciences, University of Queensland, St. Lucia, Queensland, Australia.\\
	$^2$Department of Life Sciences, Natural History Museum, Cromwell Road, London, SW7 5BD, UK. natalie.cooper@nhm.ac.uk}\\
	$^*$Corresponding author\\}

\vfill

\runninghead{Time sub-samples in disparity-through-time analyses}
\keywords{time bin, time-slice, disparity, sub-sampling}
% up to 6

\begin{document}

\mstitlepage
\parindent = 1.5em
\addtolength{\parskip}{.3em}

\section{Abstract}
% 300 words max

Disparity-through-time analyses can be used to determine how morphological diversity changes in response to mass extinctions, or to investigate the drivers of morphological change.
These analyses are routinely applied to palaeobiological datasets, yet although there is much discussion about how to best calculate disparity, there has been little consideration of how taxa should be sub-sampled through time.
Standard practice is to group taxa into discrete time bins, often based on stratigraphic periods.
However, this can introduce biases when bins are of unequal size, and implicitly assumes a punctuated model of evolution.
In addition, many time bins may have few or no taxa, meaning that disparity cannot be calculated for the bin and making it harder to complete downstream analyses.
Here we describe a different method to complement the disparity-through-time tool-kit: time-slicing.
This method uses a time-calibrated phylogenetic tree to sample disparity-through-time at any fixed point in time rather than binning taxa.
It uses all available data (tips, nodes and branches) to increase the power of the analyses, specifies the implied model of evolution (punctuated or gradual), and is implemented in \texttt{R}.
We test the time-slicing method on four example datasets and compare its performance in common disparity-through-time analyses.
We find that the way you time sub-sample your taxa can change your interpretations of the results of disparity-through-time analyses.
We advise using multiple methods for time sub-sampling taxa, rather than just time binning, to gain a better understanding disparity-through-time.

\newpage
\raggedright
\doublespacing
\setlength{\parindent}{1cm}

\section{Introduction}
Disparity-through-time analyses are common in palaeontology \citep{gould1991disparity,briggs1992morphological,Wills1994,Foote01071994}.
They reveal how the morphological diversity of clades has changed through time, and allow us to make inferences about the breadth of ecological niches extinct taxa occupied \citep{foote1997evolution}.
Results from disparity-through-time studies also provide insights into the ecological impacts of mass extinctions, competitive replacements, and the drivers of morphological evolution \citep{Brusatte12092008,Foote29111996,friedmanexplosive2010}.
Unfortunately, the way we perform these analyses may have profound effects on our conclusions.

Disparity-through-time analyses have two main analysis components: calculating disparity, and creating time sub-subsets of the data.
Here we focus on the latter.
The nature of disparity (i.e. it is a diversity metric), means it cannot be calculated using a single individual, so some way of sub-sampling taxa is required.
Changes in disparity-through-time are generally investigated by calculating the disparity of taxa present during specific time intervals or time bins \citep[e.g.][]{cisneros2010,prentice2011,Hughes20082013,hopkinsdecoupling2013,bentonmodels2014,bensonfaunal2014}.
These time bins are usually defined based on stratigraphy \citep[e.g.][]{cisneros2010,prentice2011,Hughes20082013,bentonmodels2014} but can also be arbitrarily chosen time bins of equal (or approximately equal) duration \citep{Butler2012,hopkinsdecoupling2013,bensonfaunal2014}.
However, this approach has several limitations.

First, time bins defined by stratigraphy are not of equal size, biasing higher disparity towards longer stratigraphic periods.
This can be dealt with using rarefaction methods, i.e. repeating the analysis while resampling the taxa to have the same number of taxa in each bin (e.g. using bootstrapping with limited resampling).
This can, however, lead to large confidence intervals when there are stratigraphic periods with few species.
Other studies split large time bins so they are of roughly equal size, but this is often an \textit{ad hoc} procedure that can introduce more bias depending on where bins are split.
Second, the time binning approaches (whether bins are equally sized or not) favour punctuated equilibrium modes of evolution.
Whether the disparity represents an average across the interval (with no interpretation of if or how it varied within the time bin), or it is \textit{effectively} postulated to be constant, when analysing the changes in disparity-through-time, this method will only allow changes in disparity to occur \textit{between} intervals rather than also allowing for gradual changes within intervals \citep[a pattern that is fairly common in the fossil record;][]{Hunt21042015}.
Third, when investigating changes in disparity due to events at a specific time point (e.g. a mass extinction), time bins may have not have high enough resolution to resolve changes at the event, for example if time bins are every 20 million years it may be hard to capture the effects of an event five million years into the bin.
Finally, time bin analyses are often limited by the number of taxa in each bin.
If there are insufficient taxa in a time bin, disparity cannot be calculated, so further analyses, e.g. correlations of disparity with hypothesised drivers of morphological evolution, are not possible.

To address these issues, we propose a ``time-slicing'' approach that takes advantage of the wealth of palaeontological datasets which now have associated phylogenies.
Time-slicing uses a phylogenetic tree and considers subsets of taxa at specific equidistant points in time, as opposed to considering subsets of taxa between two points in time \citep[a similar approach is outlined in][]{halliday2016eutherian}.
This results in even-sampling across time and permits us to define the underlying model of character evolution (punctuated or gradual).
Time-slicing also includes any element present in the phylogeny (branches, nodes and tips) at the time-slice in question as part of the disparity calculation.
This allows us to measure disparity at time points where there are no sampled terminal taxa, and increases the sample size at each time point, making downstream analyses of the drivers of disparity much more feasible.

Here we present our time-slicing methods using four datasets taken from the literature.
We calculate disparity-through-time for each dataset using a range of time binning and time-slicing methods, and then compare these approaches with respect to the relative disparities calculated, but also investigate how the different approaches influence biological conclusions.
We find that the choice of time sub-sampling method can have profound effects on the conclusions of disparity-through-time analyses.

\section{Materials and Methods}
\subsection{Overview}
\label{overview-section}
To test the different time sub-sampling methods, we followed the protocol below (Fig. \ref{fig:overview}).
All the code needed to reproduce these analyses (along with detailed instructions) is provided on GitHub (\url{https://github.com/nhcooper123/time-slice}).
% NC: Add proper reference to this before publication.

% overview figure
\begin{figure}[!htbp]
    \centering
    \includegraphics[width=1\linewidth, height=1\textheight, keepaspectratio]{figures/Figure1-outlinesvg.pdf}
    \caption[Outline of the disparity-through-time pipeline]
    {Outline of the disparity-through-time pipeline. 1) We use ancestral character estimation to infer nodal character states; 2) we measure the pairwise Gower distance between the tip character states and nodal character states; 3) we ordinate the distance matrix using principal coordinates analysis (PCoA/PCO); 4) we time sub-sample the PCoA matrix using time bins defined by stratigraphic periods, equally sized time bins and time-slices (using six methods to estimate ordination scores for branches); and finally 5) we measure disparity-through-time for each of these methods.}
    \label{fig:overview}
\end{figure}

\subsection{Example datasets}
\label{datasets}
To test the different time binning/slicing methods we selected four datasets: a mammal dataset from \cite{beckancient2014}, two theropod datasets from \cite{brusatte2014gradual} and \cite{bapst2016topology}, and a crinoid dataset from \cite{wright2017bayesian}.
Table \ref{table:datasets} and the Supporting Information Appendix S1 provide more details.
Each dataset consists of first and last occurrence dates for all taxa, a matrix of morphological characters in NEXUS format, and a time-scaled phylogeny.
These datasets are freely available with their accompanying papers (Table \ref{table:datasets}), but for reproducibility purposes we also provide the data we used on GitHub (\url{https://github.com/nhcooper123/time-slice/data}).

% Dataset details table
\begin{landscape}
\input{tables/table_datasets.tex}
  % Label in input or referencing doesn't work
\end{landscape}

\subsection{Preparing the data for disparity-through-time analysis}

\paragraph{Estimating ancestral character states.}
\label{ace}
For each dataset we estimated the ancestral character states at each node using the \texttt{AncStatesEstMatrix} function from the \texttt{Claddis} \texttt{R} package \citep{Claddis,R}.
This function uses the re-rooting method \citep{Yang01121995,Garland2000} to get Maximum Likelihood estimates of the ancestral states for each character at every node in the phylogeny \citep[based on the \texttt{rerootingMethod} function in \texttt{phytools};][]{phytools}.
Inapplicable and missing characters for any taxon were treated as ambiguous characters (i.e. any possible observed state for the character).
To prevent poor ancestral state estimations from biasing our results, especially when a lot of error is associated with the estimations, we only included ancestral state estimations with a scaled Likelihood $\geq$ $0.95$.
Ancestral state estimations with scaled Likelihoods below this threshold were recoded as missing data (``?'').
This allowed our results to be less dependent on the quality (or the absence thereof) of the ancestral state estimations methods, especially in parts of the datasets where data were sparse.
This approach is similar to \cite{brusatte2011phylogenetic} but uses model based estimations (rather than parsimony) allowing us to control for ambiguous (i.e. poorly estimated) nodes.

\paragraph{Building morphospaces.}
To explore disparity-through-time in our datasets, we used a morphospace approach \citep[e.g.][]{Foote01071994,Foote29111996,Wesley-Hunt2005,Brusatte12092008,friedmanexplosive2010,toljagictriassic-jurassic2013,Hughes20082013}.
Morphospaces can be obtained from any multidimensional morphological data set but can differ in the data used (e.g. discrete or continuous), and whether they include phylogenetic data or not.
Although empirical morphospaces from discrete or continuous data have been shown to have similar properties \citep{foth2012different,hetherington2015cladistic}, our morphospaces are based on discrete morphological data (originally collected for phylogenetic analysis; \textit{c.f.} geometric morphometric data) and include some phylogenetic information (see above).
Mathematically, our morphospaces are $n$ dimensional objects that summarise the distances between discrete morphological characters of the taxa present and their ancestors.

\paragraph{Constructing distance matrices.}
To estimate the morphospaces for each of our datasets we first constructed pairwise distance matrices of length $n$, where $n$ is the total number of tips and nodes in the dataset.
We calculated the $n$$\times$$n$ distances using the Gower distance \citep{Gower71}, i.e. the number of mismatched characters over the number of shared characters.
This allows us to correct for distances between two taxa that share many characters and could be closer to each other than to taxa with fewer characters in common (i.e. because some pairs of taxa share more characters in common than others, they are more likely to be similar).
For discrete morphological matrices, using this corrected distance is preferable to the raw Euclidean distance because of its ability to deal with discrete or/and ordinated characters as well as with missing data \citep{anderson2012using}.
However, the Gower distance cannot calculate distances when taxa have no overlapping data.
Therefore, we used the \texttt{TrimMorphDistMatrix} function from the \texttt{Claddis} \texttt{R} package to remove pairs of taxa with no cladistic characters in common.
This led to us removing nine taxa from the \cite{bapst2016topology} dataset, and 19 from the \cite{brusatte2014gradual} dataset, but none from the other two datasets (see Supporting Information Appendix S1 for details of which species).

\paragraph{Ordination.}
After constructing our distance matrices we transformed them using classical multidimensional scaling \citep[MDS;][]{torgerson1965multidimensional,GOWER01121966,cailliez1983analytical}.
This method (also referred to as PCO; e.g. \citealt{Brusatte2015}; or PCoA; e.g. \citealt{paradisape:2004}; but distinguished in \citealt{legendre2012numerical}) is an eigen decomposition of the distance matrix.
Because we used Gower distances instead of raw Euclidean distances, negative eigenvalues can be calculated.
To avoid this problem, we first transformed the distance matrices by applying the Cailliez correction \citep{cailliez1983analytical} which adds a constant $c^*$ to the values in a distance matrix (apart from the diagonal) so that all the Gower distances become Euclidean ($d_{Gower}+c^*=d_{Euclidean}$; \citealt{cailliez1983analytical}).
We were then able to extract $k$ eigenvectors for each matrix (representing the $k$ dimensions of the morphospace) where $k$ is equal to $n-2$, i.e. the number of taxa in the matrix ($n$) minus the last two eigenvectors that are always null after applying the Cailliez correction.
Contrary to previous studies \citep[e.g][]{brusatte50,cisneros2010,prentice2011,anderson2012using,Hughes20082013,bentonmodels2014}, we use all $k$ dimensions of our morphospaces and not a sub-sample representing the majority of the variance in the distance matrix (e.g. selecting only $x$ dimensions that represent up to 90\% of the variance in the distance matrix; \citealt{Brusatte12092008,toljagictriassic-jurassic2013}).
Note that our morphospaces represent an ordination of all possible morphologies coded in each study through time.
It is unlikely that all morphologies will co-occur at each time point, therefore, the disparity of the whole morphospace is expected to be greater than the disparity at any specific point in time.

\subsection{Disparity-through-time analyses}

Disparity-through-time analyses were performed using the \texttt{dispRity} \texttt{R} package \citep{dispRity.v02}.
% NC: Update reference here prior to publication

\paragraph{Calculating disparity.}
\label{disparity_calc}
Disparity can be calculated in many different ways \citep[e.g.][]{Wills1994,Ciampaglio2004,thorneresetting2011,hopkinsdecoupling2013,huang2015origins}, however a majority of studies in palaeobiology estimate disparity using four metrics: the sum and products of ranges and variances, each of which gives a slightly different estimate of how the data fits within the morphospace \citep{Foote01071994,Wills1994,brusatte50,Brusatte12092008,cisneros2010,thorneresetting2011,prentice2011,brusattedinosaur2012,toljagictriassic-jurassic2013,ruta2013,bentonmodels2014,bensonfaunal2014}.
However, these metrics have limitations.
First, the range metrics are affected by the uneven sampling of the fossil record \citep{Butler2012}.
Second, because we include all $k$ dimensions in the analysis (see above), the products of ranges and variances will tend towards zero since the scores of the last dimension are usually really close to zero themselves.
We therefore use the sum of variances metric to estimate disparity here:

\begin{equation}
\text{disparity} = \sum\sigma^{2}{k_i}
\end{equation}

\noindent where $\sigma^{2}{k_i}$ is the variance for the $k_i^{th}$ dimension ranging from $n$ to $n-2$ with $n$ being the number of taxa in the dataset.
Note that there are still statistical issues with this metric (such as the co-variance between dimensions not being measured), but for the purposes of comparison with previous work we decided to use a standard metric for these analyses.

\paragraph{Time sub-sampling}
\label{time_sub-samples}

To estimate disparity-through-time we first need to split the data into time sub-samples.
Here we use three time sub-sampling methods.

\begin{enumerate}
  \item Stratigraphic time bins. This is the traditional method, where all the taxa within each stratigraphic period are included in the disparity calculation. This often leads to bins of unequal duration. Here we use stratigraphic stages and epochs.
  \item Equally sized time bins. This is another commonly used method, where the time frame of interest is split into equally sized time bins, then all the taxa within each time bin are included in the disparity calculation.
  \item Time-slicing. We describe this in more detail below, but in brief, time-slicing uses a phylogeny, and rather than binning the data, it takes slices through a phylogeny and includes all the taxa and nodes in that slice within the disparity calculation.
\end{enumerate}

\paragraph{Time-slicing.}
\label{time_slicing}
The ``time-slicing'' approach considers subsets of taxa in the morphospace at specific equidistant points in time, as opposed to considering subsets of taxa between two points in time.
This results in even-sampling of the morphospace across time and allows us to use different underlying models of character evolution (punctuated or gradual).

In practice, time-slicing considers the disparity of any element present in the phylogeny (branches, nodes and tips) at any point in time.
When the phylogenetic elements are nodes or tips, the ordination scores for the nodes (estimated using ancestral state reconstruction as described above) or tips are directly used for calculating disparity.
When the phylogenetic elements are branches we choose the ordination score for the branch using one of two evolutionary models:

\begin{enumerate}

    \item{\textbf{Punctuated evolution.}}
    This model selects the ordination score from either the ancestral node or the descendant node/tip of the branch regardless of the position of the slice along the branch.
    Similarly to the time bin approach, this reflects a model of punctuated evolution where changes in disparity occur either at the start or at the end of a branch over a relatively short time period, and clades undergo long periods of stasis during their evolution \citep{Gould1977,Hunt20112007}.
    We apply this model in four ways:

    \begin{enumerate}[(i)]

      \item The ``acctran'' model, always selecting the ordination score of the descendant node/tip of the branch.
      \item The ``deltran'' model, always selecting the ordination score of the ancestral node of the branch.
      \item The ``random'' model, randomly selecting the ordination score of either the ancestor or the descendant of the branch.
      \item The ``proximity'' model, selecting the ordination score of the ancestor if the slice occurs in the first half of the branch, and the descendant if the slice occurs in the second half of the branch.

    \end{enumerate}

    The two first models assume that changes always occur early (\textit{acc}elerated \textit{tran}sition) or late along the branches (\textit{del}ayed \textit{tran}sition).
    The third model makes neither assumption and simply selects data from the ancestor or the descendant at random, and the fourth bases the selection of either the ancestor or the descendant on where the slice occurs along the branch.
    These punctuated models only select either the ordination score from the ancestor and the descendant once in the whole disparity analysis.
    For example, if using the ``random'' model, if the data of the ancestor has been randomly chosen, only this data will be used during the bootstrapping (see below) and for the disparity calculation.

    \item{\textbf{Gradual evolution.}}
    Unlike the punctuated models, the following models do not select the ordination score of either the ancestor or the descendant but associate a probability to both.
    This reflects a model of gradual evolution where changes in disparity are gradual and cumulative along the branch.

    \begin{enumerate}[(i)]
    \setcounter{enumii}{4}
      \item The ``equal splits'' model (probabilistic), selects the ordination score from both the ancestor and the descendant with an equal probability:
          \begin{equation}
          p(\text{ancestor}) = p(\text{descendant}) = 0.5
          \end{equation}

    \item The ``gradual splits'' model (probabilistic), selects the ordination score from both the ancestor and the descendant with a probability function of the distance between the nodes/tip at the ends of the branch and the slice:

          \begin{equation}
              p(\text{ancestor}) = \frac{d(\text{ancestor},\text{slice})}{d(\text{ancestor},\text{descendant})}
          \end{equation}
          \begin{equation}
              p(\text{descendant}) = 1 - p(\text{ancestor})
          \end{equation}
    \noindent where $d(x,y)$ is the distance between the two elements $x$, $y$ (ancestor, slice or descendant) measured in units of branch length.
    \end{enumerate}

    In these models, the ordination scores of both the ancestor and descendant contribute to the disparity calculation.
    For example, using the ``gradual splits'' model, if the slice occurs in the third quarter of a branch joining node A to node/tip B (75\% of the total branch length), after bootstrapping, the disparity results will be based on 25\% of the data from A and 75\% of the data from B.
    Because of the probabilistic nature of these models, they are only meaningful when calculating disparity from bootstrapped datasets.
\end{enumerate}

It is important to note that the time-slicing method is not an ancestral states estimation method \textit{per se}.
This method does not estimate values along a branch applying a model (\textit{c.f.} methods described for ancestral character estimation in the ``Preparing the data for disparity-through-time analysis'' section above) but rather chooses between the two available pieces of information (the ordination score of the descendant or the ancestor) using the methods described above.
This allows the method to be used in post-ordination analysis where the data used in each time-slice is data already present in the morphospace.
In other words, this method does not require a re-ordination of the morphospace every time a slice goes through a branch, thus allowing the properties of the morphospace (e.g. distance between species, variance of each axis, etc.) to remain constant.
For example, using the ``equal.splits'' model on an ancestor and a descendant with PCO1 values of respectively $0.04$ and $0.03$, after a sufficient number of bootstrap replicates (e.g. 100) the value along the branch will be close to $0.5\times0.04 + 0.5\times0.03 = 0.035$.
By estimating this value rather than generating it (i.e. creating a new element mid-way along the branch that would be the average of the descendant and ancestor - $0.035$) we obtain the same results without modifying the morphospace properties.

\subsection{Comparing time sub-sampling methods}

To compare the time binning and time-slicing approaches we applied the methods as follows (see Fig. \ref{figure:bins}).

\begin{figure}[!htbp]
    \centering
    \includegraphics[width=1\linewidth, height=1\textheight, keepaspectratio]{figures/Figure2-binning.pdf}
    \caption[Outline of the three time sub-sampling methods]
    {Outline of the three time sub-sampling methods.
    Stratigraphy: time sub-samples are defined by stratigraphic periods.
    Here there are five stratigraphic periods in the 20 My (million years) time frame of interest, i.e. five bins/slices with variable sizes/intervals.
    Duration: time sub-samples are defined based on the mean duration of stratigraphic periods in the time frame of interest.
    Here, the mean duration of stratigraphic periods is 5 My, so there are four bins/slices of 5 My duration (or four slices with 5 My intervals between them) in the 20 My time frame of interest.
    Number: time sub-samples are defined based on the number of stratigraphic periods in the time frame of interest.
    Here, there are five stratigraphic periods, so there are five bins/slices of 4 My duration (or five slices with 4 My intervals between them) in the 20 My time frame of interest.}
    \label{fig:bins}
\end{figure}

\begin{enumerate}
  \item Stratigraphy: time sub-samples defined by stratigraphic periods (Fig. \ref{fig:bins}).
  \begin{enumerate}[(i)]
    \item Time bins (unequal).
    We calculated disparity for the taxa in each stratigraphic period (stage or epoch).
    To reduce the influence of outliers on our disparity estimates, we bootstrapped each disparity measurement for each time bin by randomly resampling with replacement a new sub-sample of taxa from the observed taxa in the bin 100 times.
    We then calculated the median disparity value for each time bin along with the 50\% and 95\% confidence intervals.
    \item Time-slices (non-equidistant).
    We calculated disparity using our time-slicing approach with slices occurring at the midpoint of each stratigraphic period (stage or epoch), and using all six time-slicing methods (acctran, deltran, random, proximity, equal splits and gradual splits).
    To reduce the influence of outliers on our disparity estimates, we bootstrapped each disparity measurement as described above for the stratigraphic time bins.
    \end{enumerate}

  \item Duration: time sub-samples defined by the duration of stratigraphic periods (Fig. \ref{fig:bins}).
  \begin{enumerate}[(i)]
    \item Time bins (equal).
    We calculated disparity for the taxa in each time bin where time bin size was defined by the mean \textit{duration} of the stratigraphic period (stage or epoch), and bootstrapped the disparity values as described above.
    \item Time-slices (equidistant).
    We calculated disparity using our time-slicing approach where the interval between slices, was defined by the mean \textit{duration} of the stratigraphic period (stage or epoch).
    We used the six time-slicing methods and bootstrapped as described above.
  \end{enumerate}

  \item Number: time sub-samples defined by the number of stratigraphic periods (Fig. \ref{fig:bins}).
  \begin{enumerate}[(i)]
    \item Time bins (equal).
    We calculated disparity for the taxa in each time bin where the number of time bins was defined by the \textit{number} of stratigraphic periods (ages or epochs) in the time frame of interest, and bootstrapped the disparity values as described above.
    \item Time-slices (equidistant).
    We calculated disparity using our time-slicing approach where the number of slices, was defined by the \textit{number} of stratigraphic periods (ages or epochs) in the time frame of interest.
    We used the six time-slicing methods and bootstrapped as described above.
  \end{enumerate}

\end{enumerate}

We also recorded the number of taxa (or taxa and nodes for time-slicing methods) in each sub-sample as a proxy for taxonomic diversity.

\subsection{Testing for differences in the time sub-sampling methods}
\label{testing}
Testing for statistical differences among the time sub-sampling methods described above is difficult, as we need to compare similar units, and also to tackle questions important to the interpretation of disparity-through-time analyses.
We therefore present three different, simple ways of comparing the time sub-sampling methods as follows.

\paragraph{Systematic differences in disparity-through-time.}
To test whether using time bins or time-slices resulted in significantly different disparity values at common time points, we used paired Wilcoxon tests to compare the median bootstrapped disparities obtained in the stratigraphy (time sub-samples defined by stratigraphic periods), duration (time sub-samples defined by the duration of stratigraphic periods), and number (time sub-samples defined by the number of stratigraphic periods) analyses described above.

Due to the uneven spread of taxa across phylogenies, some time bins will contain one or no species, meaning that we cannot estimate disparity for that time bin.
We first, therefore, removed the time bins, and corresponding time-slices, without disparity estimates.
We then performed paired Wilcoxon tests with Bonferroni corrected p-values, so that bins and slices for the same time period are being compared.
Significant results suggest that there is a systematic difference in disparity values at each time point, depending on whether bins or slices are used.

\paragraph{Disparity peaks.}
We are perhaps more interested in how the conclusions of disparity-through-time analyses are influenced by the choice of time sub-sampling method, rather than the disparities estimated by each method \textit{per se}, especially as these will be influenced by the number of taxa (and/or nodes) included in each sub-sample.
Therefore, we also investigated where peaks of disparity occurred in each of our datasets for the different time sub-sampling methods.
We calculated the maximum bootstrapped disparities for each dataset and for each time sub-sampling method, along with their associated confidence intervals.
Significant shifts in disparity peaks suggest that the choice of time sub-sampling method will influence our conclusions about relative changes in the disparity of our groups through time.

\paragraph{Effects of mass extinction events.}
Many analyses of disparity-through-time aim to demonstrate differences in disparity before and after mass extinction events.
Two of our four datasets contain taxa before and immediately after a mass extinction (Cretaceous-Paleogene 66 MYA; \citealt{beckancient2014}; Ordovician-Silurian 455-430 MYA; \citealt{wright2017bayesian}), so we used Wilcoxon tests with Bonferroni corrected p-values to compare disparity in the time bin/slice prior to the appropriate mass extinction, to that of the time bin/slice following the extinction event.
Significant results suggest an effect of the mass extinction on disparity in the group.
We then compare these results across the time sub-sampling methods to determine if our conclusions change depending on the method used.
We repeated these analyses using the two time bins/slices after the one immediately following the mass extinction event to account for any lag effects of the mass extinction on disparity.

\section{Results}

\subsection{Disparity-through-time analyses}

Disparity changes through time for each of our four datasets (Fig. \ref{figure:dtt1}, Supporting Information Appendix S2: Figs A1-A2).
Relative disparities tend to be lower with time binning methods, likely because these contain fewer taxa than time-slicing methods.
The six different time-slicing methods (acctran, deltran, random, proximity, equal splits and gradual splits) show similar patterns, so we focus only on the results for one method with a punctuated model of evolution (specifically the `proximity' method), and one method with a gradual model of evolution (specifically the `gradual splits' method).
Results for all six methods can be found in Supporting Information Appendix S2: Figures A1-A2.

\begin{figure}[!htbp]
    \centering
    \includegraphics[width=1\linewidth, height=1\textheight, keepaspectratio]{figures/Figure3-fig-dtt-epoch-revision.pdf}
    \caption[Relative disparity-through-time.]
    {Relative disparity-through-time. Median bootstrapped disparities were calculated using time binning and time-slicing approaches.
    Green points represent time binning methods, purple points are time-slices with a punctuated model of evolution (`proximity' method), and blue points are time-slices with a gradual model of evolution (`gradual splits' method).
    Relative disparities (median bootstrapped disparity divided by the maximum median bootstrapped disparity for a dataset and analysis method) are presented so they can be compared across datasets/methods.
    Stratigraphy uses unequal time bins or non-equidistant time-slices, where the width of the bin, or the interval between slices, is equivalent to stratigraphic epochs.
    Duration uses equal time bins or equidistant time-slices, where the width of the bin, or the interval between slices, is the mean duration of stratigraphic epochs in the time frame of the dataset.
    Number uses equal time bins or equidistant time-slices, where the number of bins, or the number of slices, is the mean number of stratigraphic epochs in the time frame of the dataset.
    In all cases, time bin disparities are plotted at the midpoint of the bin, and error bars represent the 95\% confidence intervals around the bootstrapped median disparity.
    The four dataset names are on the first plot for each dataset (see Table \ref{table:datasets} for details).
    Results for stratigraphic stages, and for other time-slicing methods, are in the Supporting Information Appendix S2: Figures A1-A2.}
    \label{figure:dtt1}
\end{figure}

\subsection{Testing for differences in the time sub-sampling methods}

\paragraph{Systematic differences in disparity-through-time.}
There is no overall significant systematic difference among the disparities calculated using time bins and those calculated using the time-slicing methods (Table \ref{table:wilcox}, Supporting Information Appendix S2: Table A1).
Instead, the differences depend on the dataset and method in question.
For example, the \cite{brusatte2014gradual}, \cite{bapst2016topology} and \cite{wright2017bayesian} datasets, show significant differences when using bins versus time-slices defined by stratigraphy, but the \cite{beckancient2014} dataset appears robust to these different approaches.
Likewise, the \cite{beckancient2014}, \cite{brusatte2014gradual} and \cite{bapst2016topology} datasets have different disparities when the number of bins or slices is the mean number of stratigraphic periods, but this is not seen in the \cite{wright2017bayesian} dataset.
Note that for epochs, we find fewer significant differences simply because the smaller numbers of bins and slices being compared means we have low power to detect a significant difference.

% Wilcoxon results table
\input{tables/table_wilcox.tex}

\paragraph{Disparity peaks.}
In the \cite{beckancient2014} and \cite{bapst2016topology} datasets, disparity peaks occur much at much older ages when time-slicing rather than time binning approaches are used (Fig. \ref{figure:peak1}; Supporting Information Appendix S2: Figs A3-A4).
This is also true for stratigraphic time bins in the \cite{wright2017bayesian} dataset, although when using equal time bins the peaks are later than the time-slicing methods, or very similar (Fig. \ref{figure:peak1}; Supporting Information Appendix S2: Figs A3-A4).
Across the three time binning methods, the \cite{brusatte2014gradual} dataset has similar disparity peaks whichever method is used, the \cite{wright2017bayesian} dataset only had variation in peaks when using unequal time bins (stratigraphy), whereas in the \cite{bapst2016topology} and \cite{beckancient2014} datasets, stratigraphic (unequal) versus equally sized time bins make a large difference to where the disparity peak occurs (Fig. \ref{figure:peak1}; Supporting Information Appendix S2: Figs A3-A4).
Additionally, there seem to be small discrepancies within the time-slicing methods (gradual.split \textit{vs.} proximity) except in the \cite{beckancient2014} dataset where the gradual split model recovered disparity peaks at younger ages than the proximity model (Fig. \ref{figure:peak1}; Supporting Information Appendix S2: Figs A3-A4)

% Plots showing the differences
\begin{figure}[!htbp]
    \centering
    \includegraphics[width=1\linewidth, height=1\textheight, keepaspectratio]{figures/Figure4-fig-peaks-epoch-revision.pdf}
    \caption[Timing of peak disparity for four example datasets.]
    {Timing of peak disparity.
    Median bootstrapped disparities were calculated using time binning and time-slicing approaches.
    Green points represent time binning methods, blue points are time-slices with a punctuated model of evolution (`proximity' method), and purple points are time-slices with a gradual model of evolution (`gradual splits' method).
    Stratigraphy uses unequal time bins or non-equidistant time-slices, where the width of the bin, or the interval between slices, is equivalent to stratigraphic epochs.
    Duration uses equal time bins or equidistant time-slices, where the width of the bin, or the interval between slices, is the mean duration of stratigraphic epochs in the time frame of the dataset.
    Number uses equal time bins or equidistant time-slices, where the number of bins, or the number of slices, is the mean number of stratigraphic epochs in the time frame of the dataset.
    For time bins the points indicate the maximum and minimum ages of the time bin within which peak disparities appeared.
    The four dataset names are on the first plot for each dataset (see Table \ref{table:datasets} for details).
    Results for stratigraphic stages, and for other time-slicing methods, are in the Supporting Information Appendix S2: Figures A3-A4.}
    \label{figure:peak1}
\end{figure}

\paragraph{Effects of mass extinction events.}
Mass extinction events influence disparity in both the \cite{beckancient2014} and \cite{wright2017bayesian} datasets (Fig. \ref{figure:extinction1}).
However, whether this change in disparity is significant or not depends on the method used to create time sub-samples (Fig. \ref{figure:extinction1}), and whether stages or epochs are used.
In general, for the \cite{beckancient2014} dataset, time binning tended to give more significant results than time-slicing methods, but this was not the case for the \cite{wright2017bayesian} dataset.

% Big heat map table thing
\begin{figure}[!htbp]
    \centering
    \includegraphics[width=1\linewidth, height=1\textheight, keepaspectratio]{figures/Figure5-fig-extinction-effect.pdf}
    \caption[Extinction example datasets.]
    {Effects of mass extinction events on disparity.
    Pink cells and blue cells indicate respectively a significant or non significant change in disparity before and after the mass extinction event (Cretaceous-Paleogene 66 MYA; \citealt{beckancient2014}; Ordovician-Silurian 455-430 MYA; \citealt{wright2017bayesian}).
    \textit{e:1}, \textit{e:2}, and \textit{e:3} denote whether the comparison was between the time bin or time-slice immediately after the mass extinction (\textit{e:1}), or the second (\textit{e:2}) or third (\textit{e:3}) bin/slice after the mass extinction to account for any lag effect.
    The top seven rows use stratigraphic stages and the bottom seven rows use stratigraphic epochs.
    Labels on the left hand side indicate whether time bins (``bins'') were used, or which of the six time-slicing methods was used.}
    \label{figure:extinction1}
\end{figure}

\section{Discussion}
Disparity-through-time analyses are influenced by the choice of time sub-sampling method used to divide the taxa.
While differences in the relative disparities calculated among time sub-sampling methods may not be of much biological importance, these changes can have important implications for the conclusions of downstream analyses.
For example, using stratigraphic epochs as our reference time period, there are 21 potential methods for time sub-sampling our data (splitting by stratigraphy, number and duration, see methods, and using time bins or one of six time-slicing methods).
Of these 21 methods, in 16 (76\%) we show that placental mammals \citep{beckancient2014} significantly increased in disparity in the time bin/slice immediately after the K-Pg mass extinction event, and in 20 (95\%) we show that crinoids \citep{wright2017bayesian} significantly decreased in disparity in the time bin/slice immediately after the O-S mass extinction event.
Given the high congruence (76\% and 95\%) of these results, one could argue that time-sub-sampling methods are not important.
However, if we had chosen to investigate crinoid disparity only using time bins and splitting these so the number of time bins was equal to the number of epochs (number), we would have concluded that the O-S extinction had no effect on crinoid disparity.
Likewise, the timing of peak disparity differs among methods.
This is particularly evident when comparing stratigraphic time bins to time slicing methods, where for most of our datasets we see a much later time to peak disparity.
This could have major implications for our understanding of how morphological diversity changes through time, for example in response to climate.
These results highlight the sensitivity of disparity-through-time analyses to the choice of time sub-sampling method.
Fortunately this issue is easy to solve; either disparity-through-time analyses should use, and report results from, multiple time sub-sampling methods (as demonstrated here), or great care should be taken in determining the appropriate time sub-samples to answer the question of interest.

Time-slicing has several advantages over time binning (using either equally or unequally sized bins) approaches.
First, it allows us to use as much of the information available to us, in the form of phylogenetic relationships and ancestral taxa, as possible.
This increases our ability to investigate key biological questions, such as how do various drivers influence morphological diversity through time, and how do mass extinctions influence disparity \citep{Brusatte12092008,Foote29111996,friedmanexplosive2010}, both by increasing the statistical power of analyses and through the availability of data at key time points in the history of our groups.
Second, we are able to be more explicit about the mode of evolution in our clades; in time-slicing we can apply punctuated or gradual models of trait change rather than making an assumption of punctuated evolution.
This may be important, as gradual change is a common pattern of trait evolution in the fossil record \citep{Hunt20112007}.

Of course the method also has limitations.
The main one of these is a practical one; it requires a time-calibrated phylogeny and these are not available for all palaeontological datasets.
Furthermore, like most phylogeny based methods, time-slicing depends on ancestral state estimations.
Care should be taken in interpreting these, as they are highly dependent on the data and models used for the estimations \citep{Slateretal2012,doi:10.1080/10635150801910451}.
The difference between the time-binning and time-slicing results could also simply be due to the nature of the fossil record.
Rates of sedimentation vary in time and space influencing the groups found within the rock record and their temporal distribution.
In this case, different beds could represent different ``packages'' of fauna through time separated by gaps, resulting in natural ``bins'' rather than slices of the data.
Slicing through such strata will yield similar results no matter where in time the slice occurs.
It is important to note however, that the time slicing method also includes ancestral estimations (either through the nodes or the branches) that are by definition not available in the fossil record and thus are not influenced by its nature.
Additionally, this effect is likely to be most obvious in groups where the fossil record is ``patchy'', e.g. vertebrates, but less problematic for groups with a more continuous record like Foraminifera.
Finally, \cite{Hunt21042015} found that time series are best characterized by gradual directional changes (biased random walks).
In fact, homogeneous directional changes are more likely to be supported than heterogeneous ones (e.g. punctuated changes) in longer duration series with few samples in each series.
In our implementation of time-slicing, the models are not selected based on any model fit criterion (e.g. AIC) but merely on researcher assumptions.
We thus suggest that both types of models (punctuated and gradual) are tested during analysis, unless there is strong independent support for one or the other.

\section{Conclusions}
The choice of time sub-sampling method can alter the conclusions we obtain from disparity-through-time studies.
Time-slicing methods, with explicit models of evolution, provide an alternative to traditional time binning approaches.
Note that while we introduce the time-slicing methods here, and describe their advantages, we do not suggest time-slicing is necessarily the ``best'' method for time sub-sampling in all cases.
As with all methods, the choice of methodology should be appropriate for the question and data at hand.
However, we do strongly recommend performing disparity-through-time analyses using a series of appropriate time sub-sampling methods, and comparing these, to ensure that results are not merely a consequence of the time sub-sampling method employed.

\section{Data archiving statement}
Data for this study are available from \cite{dryad_n2g80,dryad_6hb7j,dryad_84t75,beckancient2014} but for ease of reproducibility can also be found on GitHub (\url{https://github.com/nhcooper123/time-slice/data}). % add proper link to zenodo
Supporting information can be found on Dryad at https://doi.org/10.5061/dryad.vp4q518 \citep{data-time-slice}.

\section{Supporting information}
Additional Supporting information can be found on Dryad at: https://doi.org/10.5061/dryad.vp4q518 \citep{data-time-slice}.\\
Appendix S1: Extra details of datasets. \\
Appendix S2: Additional figures. \\
Appendix S3: Additional tables.

\paragraph{Acknowledgments.}
  NC thanks Mark Sutton and Philip Mannion for the invitation to contribute to the `Evolutionary Modelling' symposium at The Palaeontological Association Annual Meeting 2017.
  TG acknowledges support from the Australian Discovery Project Grant number DP170103227 awarded to Vera Weisbecker.
  We thank Dave Bapst, Graeme Lloyd, April Wright and David Wright for assistance in gathering data for the analyses and/or discussions about the approach; and Steve Brusatte, Sally Thomas and one anonymous reviewer for helpful comments on the manuscript.

\section{References}

\bibliographystyle{palaeo}
\bibliography{time-refs}

\end{document}


text:
\documentclass[11pt,a4paper]{report}
\usepackage{fullpage}
\usepackage{hyperref}

\renewcommand{\baselinestretch}{2}
\author{Abdalraheem Abdullah Yousef Ijjeh \\
	Supervised by: \\
	Professor Pawel Kudela}
\title{\Huge Quarterly report \\
	\huge For the period \\ 2021-03-02 to 2021-06-02
	\\ \Large In accordance with Project\\
	Feasibility studies of artificial intelligence-driven diagnostics\\ Institute of Fluid-Flow Machinery, Polish Academy of Sciences
	\\ Gdańsk, Poland}

\begin{document}
	\maketitle
	% \tableofcontents
	\newpage
	\section{Project work flow}
	\begin{itemize}
		\item Several datasets for RNN models were prepared (of full wavefield frames).
		\item Different RNN/LSTM and CNN models (using full wavefield frames) were implemented and tested to estimate the size and location of delamination.
		\item Second chapter (monograph) entitled: "Data-driven based approach for damage
		detection" was submitted and accepted by the TSD of the Polish Academy of Sciences.
		\item Currently, I am working on implementing models that artificially generate Full wavefield frames of Lamb waves resembling its propagation in CFRP and its interaction with delamination and the boundaries of the specimen.
		Accordingly, we will be able to create a large dataset of propagating Lamb waves interacting with delaminations of various locations, sizes, and angles in a short time compared with the current numerical method.
		Furthermore, this artificially created dataset can be used in training other deep learning models.
		Moreover, the idea can be generalized to further types of damage, or in general to study different phenomena.
	\end{itemize}
%	\bibliography{quarter_report.bib}
%	\bibliographystyle{unsrt}
\end{document}

ext: tex

max_stars_repo_name: IFFM-PAS-MISD/aidd

max_issues_repo_name: IFFM-PAS-MISD/aidd

max_forks_repo_name: IFFM-PAS-MISD/aidd


text:
\chapter*{Colophon}

\begin{center}
\parbox{200pt}{\raggedright\lettrine[lines=3,slope=-2pt,nindent=-4pt]{\textcolor{ISIMA-Gold}{T}}{his report was typeset} using \LaTeX,was modified by Yann Feunteun originally developed by Anis Ayari and Leslie Lamport and based on Donald Knuth's for Harvard \TeX. The body text is set in 11 point Arial, designed by Microsoft. A template, which can be used to format an ISIMA report, has been released under the permissive \textsc{mit} (\textsc{x}11) license, and can be found online at \href{https://github.com/yafeunteun/ISIMA-Latex}{github.com/yafeunteun/ISIMA-Latex} or from the author at \href{mailto:yannfeunteun@gmail.com}{yannfeunteun@gmail.com}
}
\end{center}

max_stars_repo_name: yafeunteun/ISIMA-Latex

max_issues_repo_name: yafeunteun/ISIMA-Latex

max_forks_repo_name: yafeunteun/ISIMA-Latex


text:
%!TEX root=report.tex
\section{Appendix}

\subsection{Time series analysis results}
\label{first-R}

The first model is based on interpretations of graphs (R-output):

\begin{minipage}{1.3\textwidth}
\begin{lstlisting}
Series: xdata
ARIMA(3,1,3)(1,1,4)[36]

Coefficients:
          ar1     ar2     ar3     ma1      ma2     ma3    sar1     sma1    sma2
      -0.3860  0.1239  0.0471  0.2225  -0.2783  0.0035  0.0926  -0.9147  0.1808
s.e.   0.1556     NaN     NaN  0.1673      NaN     NaN     NaN      NaN  0.0508
        sma3     sma4
      0.0980  -0.1206
s.e.  0.0466   0.0749

sigma^2 estimated as 0.0005546:  log likelihood=727.89
AIC=-1431.78   AICc=-1430.76   BIC=-1386.6

$AIC
[1] -6.435554

$AICc
[1] -6.427381

$BIC
[1] -7.315823

Warning message:
In sqrt(diag(x$var.coef)) : NaNs produced
\end{lstlisting}
\end{minipage}

max_stars_repo_name: AndreasMadsen/grace

max_issues_repo_name: AndreasMadsen/grace

max_forks_repo_name: AndreasMadsen/grace


text:
\section{Model}
\subsection{Basic Functions}
The netwok model relies on three types of entities (at present, as we model more complex functionality we might get to
other things):

\begin{itemize}
\item $node$: These are network elements including endhosts, switches (when modelled), middleboxes, etc. In the rest of
the document identifiers $e, e_1, e_2, \ldots$ are assumed to be $node$s.
\item $address$: Addresses assigned to nodes. In the rest of the document identifiers $a, b, a_1, a_2\ldots$ are assumed
to be $address$es
\item $packet$: Packets. We currently model them as tuples with the following fields:
\begin{itemize}
\item $src$: $src\in address$.
\item $dest$: $dest\in address$.
\item $origin$: $origin\in node$. This is a pseudofield in the packet tracking its true origin (which is important for
answering questions about isolation).
\item $body$: $body\in \mathbb{N}$. Payload
\item $seq$: $seq\in \mathbb{N}$. Sequence number.
\item $options$: IP Options: represented as integer.
\item $src_port$: Source port
\item $dest_port$: Destination port
\end{itemize}
For the rest of the document identifiers $p, p_1,\ldots$ represent packets.
\item $etype = \left\{R, S\right\}$ The event type, this is used to record times at which packets are received and sent. $R$ represents
the event of a node receiving a packet, while $S$ represents a node sending a packet. Event type variables are represented as $E$ below.
\item $time \subset \mathbb{N}$: Time/ordering information for an event, represented by the variable $\tau$.
\item Events: Events are represented by the tuple $\left( e, p, E \right)$ and represent the act of packet $p$ either being received at or
sent from node $e$.
\end{itemize}

We also create a few basic functions to model relationships between these types and also network functionality. These
are:
\begin{itemize}
\item $hostHasAddr: e, a \rightarrow boolean$: True if $node$ $e$ is assigned address $a$.
\item $addrToHost: a \rightarrow e$: Returns $node$ $e$ associated with address $a$. (Note this implies that the address
to node mapping is onto).
\item $send: e_1, e_2, p \rightarrow boolean$ True if $e_1$ sent $e_2$ packet $p$.
\item $recv: e_1, e_2, p \rightarrow boolean$ True if $e_2$ received packet $p$ from $e_1$.
\item $etime: e, p, E \rightarrow \tau$ The time at which event $\left( e, p, E \right)$ occurs.
\item $failed: e, \tau \rightarrow boolean$ True if $e$ is failed at time $\tau$
\end{itemize}

Given these functions and types we layer them on to produce the test network. We go bottom up starting with the network:

\subsection{Network Layer}
\subsubsection{Basic Conditions}
These are basic conditions we assume for the entire network:
\begin{itemize}
\item $\forall e_1, e_2, p:\ recv(e_1, e_2, p) \iff send(e_1, e_2, p)$: No packet loss
\item Don't consider loopback packets, we have really no control over them.
\begin{align*}
\forall e_1, e_2, p:& send(e_1, e_2, p) \Rightarrow e_1 \neq e_2\\
\forall e_1, e_2, p:& recv(e_1, e_2. p) \Rightarrow e_1 \neq e_2\\
\forall e_1, e_2, p:& send(e_1, e_2, p) \Rightarrow p.src \neq p.dest\\
\end{align*}
\item And a set of rules for time:
\begin{itemize}
\item Record the event that a packet is received and make sure packets are not received before being sent:
\begin{align*}
    \forall e_1, e_2, p:& recv(e_1, e_2, p)\Rightarrow etime(e_2, p, R) > 0 \land etime(e_2, p, R) > etime(e_1, p, S)
\end{align*}
\item Record the event that packet is sent. Make sure packets are received after being sent (yes this is redundant given above)
\begin{align*}
    \forall e_1, e_2, p:& send(e_1, e_2, p)\Rightarrow etime(e_1, p, S) > 0 \land etime(e_2, p, R) > etime(e_1, p, S)
\end{align*}
\item Make sure there are no other causes of receive events
\begin{align*}
    \forall e_1, p:& \neg(\exists e_2:\ recv(e_2, e_1, p)) \implies etime(e_1, p, R) = 0
\end{align*}
\item There are similarly no other causes for send events
\begin{align*}
    \forall e_1, p:& \neg(\exists e_2:\ send(e_1, e_2, p)) \implies etime(e_1, p, S) = 0
\end{align*}
\item No sends or receives when failed
\end{itemize}
\end{itemize}

\subsection{End Hosts}
End hosts are just nodes with a couple of extra constraints. For endhost $e$:
\begin{itemize}
\item $\forall e_i, p:\ send(e, e_i, p) \Rightarrow hostHasAddr(e, p.src)$: An end host does not forward packets on anyone
elses behalf. Note that this prevents host spoofing. One could obviously turn this off, but...
\item $\forall e_i, p:\ send(e, e_i, p) \Rightarrow p.origin = e$ Track origin correctly.
\end{itemize}

\subsection{Stateless Firewall}
Stateless firewalls just implement ACLs based on packet source and destination. ACLs are specified as a table $A \subset
address\times address = \left\{(a, b)\right\}$ which specifies the set of addresses that are denied (it is equally easy
to build one based on allowed addresses). Given a firewall $f$ and a table $A$, the firewall model says:

\begin{itemize}
\item The correct sending condition from equation~\ref{eq:sanesend}.
\item $\forall e_1, p:\ send(f, e_1, p) \Rightarrow \exists e_2 recv(e_2, f, p)$: The firewall does not invent packets.
\item $\forall e, p:\ send(f, e, p) \Rightarrow (p.src, p.dest) \not \in A \land (p.dest, p.src) \not \in A$ Sending the packet is not disallowed by an
ACL.
\end{itemize}

\subsection{Stateful Firewall}
Stateful firewally can cache previous decisions. In particular this means we don't need to have ACLs that are symmetric
(one side can punch a hole through the firewall). Cached rules are modeled as a function: $cached: a_1, a_2 \rightarrow
boolean$ which is true when packets with either source $a_1$ and destination $a_2$ or source $a_2$ and destination $a_1$
are allowed through the firewall. We also track the time at which a rule gets cached using a separate function $ctime:
a_1, a_2 \rightarrow \tau$ which models the time at which rules allowing traffic between $a_1$ and $a_2$ get cached.
(This also allows us to potentially model expiration times).

The logical model for a stateful firewall $f$ with ACL table $A$ is:

\begin{itemize}
\item The correct sending condition from equation~\ref{eq:sanesend}.
\item ACL based constraint for caching
    \begin{align*}
        \forall a, x, b, y:\ cached(a, x, b, y) \iff \\
        \ \exists e, p:\ recv(e, f, p)\\
        \ \ \land p.src = a\\
        \ \ \land p.src\_port = x\\
        \ \ \land p.dest = b\\
        \ \ \land p.dest\_port = y\\
        \ \ \land ctime(a, x, b, y) = etime(f, p, R)\\
        \ \ \land \lnot \left(a, b\right) \in A
    \end{align*}
\item Do not invent packets
     \begin{align*}
         \forall e, p:\ send(f, e, p) \implies\\
         \ \exists e_2:\ recv(e_2, f, p)\\
         \ \ \left(\nexists e_3:\ e_3 \neq e \land send(f, e_3, p)\right)\\
         \ \ etime(f, p, R) < etime(f, p, S)
     \end{align*}
\item Actual sends
        \begin{align*}
            \forall e, p:\ send(f, e, p) \implies\\
            \ (cached(p.src, p.src\_port, p.dest, p.dest\_port) \land\\
             \  \     ctime(p.src, p.src\_port, p.dest, p.dest\_port) \leq\\
             \ \ \    etime(f, p, R))
             \ \lor (cached(p.dest, p.dest\_port, p.src, p.src\_port) \land\\
             \  \     ctime(p.dest, p.dest\_port, p.src, p.src\_port) \leq\\
             \ \ \    etime(f, p, R))
        \end{align*}
\end{itemize}

\subsection{Web Proxy}
Web proxies modify packet headers (and were our original source of violations). Given a web proxy $w$ we model it as
follows:
\begin{itemize}
\item $\forall e, p:\ send(w, e, p) \implies hostHasAddr(w, p.src)$ Send all packets so source address belongs to the
proxy, thus allowing caching.
\item $\forall e_1, p_1:\ send(w, e, p) \implies \exists e_2, p_2:\ recv(e_2, w, p_2) \land p_2.origin = p.origin \land
p_2.dest = p.dest \land hostHasAddr(p_2.origin, p_2.src) \land etime(w, p_1, S) \geq etime(w, p_1, R)$ rules for what
packets a webproxy sends. These are actually somewhat incomplete (we are not accounting for responses here). The
reason is that currently I model all endhosts as endhosts instead of as servers and client, which would make this
somewhat more reasonable.
\end{itemize}

\subsection{Compression Box}
First we define a couple of global functions (actually we can define many of them they are just a kind of algorithm)
$compress$ and $uncompress$. Compression box is represented by $c$.

\begin{itemize}
    \item $\forall x, y \in \mathbb{Z}: compress(x) = y \implies uncompress(y) = x$
    \item
        \begin{align*}
            \forall e_1, p_1:\ send(c, e_1, p_1) \implies\\
            \ \exists e_2, p_2:&\  recv(e_2, c, p_2)\\
                               &\land PacketHeaderEqual(p_1, p_2)\\
                               &\land p_1.body = compress(p_2.body)\\
                               &\land \left(\nexists e_3:\ e_3 \neq e_1 \land send(c, e_3, p_1)\right)\\
                               &\land etime(c, p_1, S) > etime(c, p_2, R)
        \end{align*}
\end{itemize}

max_stars_repo_name: apanda/modeling

max_issues_repo_name: apanda/modeling

max_forks_repo_name: apanda/modeling


text:
\chapter{Future Work and Conclusion}
\label{chapter:conclusion}

Our work in synthesizing term rewriting systems suggests several further directions. We detail some intended future work here and conclude.

\section{Future Work}

\subsection{Replacing the Halide variable solver}
\label{sec:varsolverreplacement}

In our evaluation (section~\ref{sec:varsolvereval}), we showed that we were able to synthesize effective term rewriting systems to address the variable solver's task. However, our ultimate, practical aim is to synthesize a term rewriting system that can replace the existing variable solver in the Halide compiler altogether. There are three main issues that need to be addressed before we can achieve that goal.

\subsubsection{Support for rule predicates}
In our previous evaluation, we did not attempt to synthesize predicates that must be fulfilled for a rule to be applied, as we did in the simplifier work. For that reason, we replaced all constant values in our benchmarks with fresh variables. However, calls to the variable solver frequently involve constants that are known at compiler time, as seen in the \texttt{TrimNoOps} example in section~\ref{sec:thevariablesolver}. Rules that can exploit information about those constants to apply further rewrites could be crucial in solving expressions that our synthesized TRSs cannot solve now.

Our process for synthesizing rule predicates for the simplifier was fairly involved and computationally expensive. However, we believe we could take a much simpler approach for the variable solver. For the simplifier, we first synthesized a RHS for a LHS expression that contained constants, then replaced those constants with fresh variables. If the equivalence between LHS and RHS no longer held, we attempted to synthesize a predicate that implied that equivalence. For the variable solver, rather than attempt to synthesize arbitrary predicates, we believe we could get similar performance with a small set of fixed predicates, such as $c > 0$, $c < 0$, $c \neq 0$, and conjunctions of the above. For each candidate LHS, we might choose one of those fixed predicates and a (symbolic) constant in the LHS to apply the predicate to, and ask the solver to find a RHS that is equivalent to the LHS assuming that predicate is true. This would multiply the number of queries to the solvers up to the size of the predicate set, but for a relatively small set this is still highly tractable.

\subsubsection{Support for division and modulo}
Integer division and modulo operations are difficult for SMT solvers such as Z3 to reason about, as discussed in section~\ref{sec:rhssynthesis}. In earlier experiments, candidate LHSs that contained division and modulo operations almost always resulted in solver queries that timed out, and almost no RHS expressions that contained those operations were ever synthesized. For that reason, in the variable solver experiments in section~\ref{sec:varsolvereval}, we skipped all candidate LHSs that contained division or modulo and removed those operations from the language grammar from which RHSs were synthesized. This saved large amounts of processing time and made virtually no difference to the synthesized results. However, the variable solver is frequently invoked on expressions that contain those operations, so it will be important for the variable solver TRS to be able to handle them. (Indeed, the TRS based on Coq axioms described in section~\ref{sec:varsolvereval} included many rules which used division and modulo, which may be why it was able to solve many benchmarks that our synthesized TRSs could not, even though it solved few benchmarks than the synthesized TRSs overall.) We plan to investigate more performant ways to reason about these operations, such as approximating integers with bitwidth encodings.

\subsubsection{Conclusive benchmarks}
Our evaluation showed that a relatively small, randomly sampled training set is sufficient to learn a ruleset that is capable of solving a sizable fraction of a much larger test set. However, in order to propose a TRS that can serve as a substitute for the existing variable solver, we need a benchmark set that acts as a good proxy for realistic workloads the Halide compiler will encounter in the wild. We propose to work with Halide developers to gather such a benchmark set. We also did not demonstrate that a TRS-based variable solver with greater solving power would produce beneficial downstream effects in compiled pipelines, as we did for the simplifier; such a demonstration would be a useful argument for replacing the existing code with a synthesized TRS.

No matter how confident we are in our choice of benchmarks, it is difficult to draw a bright line and declare a synthesized TRS to be good enough---it is always possible to continue learning new rules, as our problem remains undecidable in the general case.

\subsection{Completion}
\label{sec:completion}

We know by observation that the Halide simplifier TRS cannot prove certain equalities
that are in fact true, or reduce certain expressions that can be further simplified.
Our goal is to learn new rules that would strengthen the TRS and allow it to make
further progress on these ``stuck'' expressions. This goal seems similar to that of
\emph{completion}, which constructs a decision procedure through syntatic rewrites
for a set of identities. We do not use completion directly, although
our synthesis algorithm could be considered analogous to completion in some ways.

In the standard Knuth-Bendix completion algorithm~\citep{knuth1983simple}, we take a
finite set of identities $E$ and a reduction order $>$ on terms as inputs; if successful, the
algorithm will return a finite, convergent set of rules $R$ that is equivalent to
$E$. The algorithm may also fail, or fail to terminate. At each step the algorithm
maintains a set of identities and a set of rules, both of which can be updated; the
algorithm may transform an identity into a new rule, find a new identity as a
consequence of the ruleset, or use the present ruleset to refine either an identity
or a rule. The algorithm runs until the ruleset has converged; specific implementations
may use some conditions under which to terminate with a failure.

No finite set of identities exists for the theory of integers. We could
fix a set of identities to use in a completion procedure, but choosing these axioms
is a non-trivial task. One issue is that the theory contains axioms such as commutativity;
an identity such as $x + y \equiv y + x$ cannot be oriented by any possible reduction
order, so our completion algorithm cannot make use of this fact. Another is that any
sufficiently powerful set of identities would almost certainly result in a non-terminating
completion procedure.  In addition, even if we use a subset of the Halide TRS for our
identities (thus yielding a confluent Halide TRS), our experience shows that many
failures in the compiler's use of the TRS are due to missing semantic facts that are not
derivable from the current Halide ruleset.

%%% SAK: I integrated this into a sentence above.  Not sure we need this much text.
%% \modified{One approach might be to transform the existing Halide TRS into a set of undirected
%% identities $E_R$ and use it together with a reduction order as inputs to a completion procedure.
%% If successful, this would yield a confluent TRS, for which any sequence of rewrites to
%% an expression would ultimately yield the same result. (It is not immediately clear how
%% to handle rule predicates in this case, however.) While this is certainly desirable,
%% our experience is that many proof failures in the current Halide TRS is not due to a
%% lack of confluence, but because of knowledge of the Halide expression langauge semantics
%% not encompassed by the identities in $E_R$.}

%Instead, our synthesis algorithm seeks to learn rules which cannot (necessarily) be
  %deduced from the existing ruleset.
In the absence of a finite set of identities, we
treat an SMT solver as a decision procedure to determine if a suggested identity
holds in our theory (of course, the solver itself is also incomplete; we only make use of the
soundness of the solver and cannot derive any information in the case where the solver cannot
show that an identity holds). If the identity holds and can be oriented using our reduction order,
it is added as a rule. It is possible that the newly-synthesized rule may be a consequence
of the existing ruleset and thus could have been found by running completion,
but we know that many synthesized rules contain information that is previously
unknown to the TRS.

If we consider our solver as standing in for an infinite set of identities that make up
our theory, we clearly could synthesize an infinite number of rules. Here we make use
of the fact that expressions encountered by the compiler have some bias
and are not sampled randomly from the entire expression space. In a preliminary
experiment, we tried generating LHSs at random within a certain expression size and
synthesizing RHSs to serve as new rules. We were able to find an extremely large number of
``missing'' rules not represented in the current TRS, but the new ruleset had
no measurable performance impact on benchmarks. Thus, we only synthesize rules if their LHS could be
applied to at least one expression observed by the compiler under realistic usecases. %Even
%without a formal definition of the bias on encountered expressions, we are able to
%target synthesis to produce new rules that have some impact on compiler output.

%\sak{Should we say something about how this applies to things other than Halide?}

As we observed in the variable solver evaluation, adding a new rule is not guaranteed to strictly increase the number of benchmarks a TRS can solve. Particularly, we saw a case in which the synthesized goal form TRS was able to solve a particular expression when the synthesized gradual TRS was not; the application of a rule that only existed in the gradual TRS blocked the use of a rule that could have fully solved the expression. This occurred because the gradual TRS was not \emph{confluent}. In a confluent ruleset $R$, for any expression $s$, if some sequence of rule applications from $R$ can transform $s$ to $t$, and some other sequence of rule applications can transform $s$ to $t'$, then both $t$ and $t'$ must be eventually rewritten by $R$ to the same form $u$. In other words, $R$ will always normalize its inputs to the same form no matter what order rewrites are applied in.

If a ruleset $R$ is finite and terminating, then a decision procedure exists to determine whether or not $R$ is confluent. Thus, we are able to detect if adding a new rule to a ruleset we are synthesizing would change it from confluent to non-confluent, giving rise to potential issues like the failing benchmark above. If we can detect these issues, can we address them at the time the rule is synthesized?

One strategy is suggested by \emph{completion} algorithms. A completion algorithm takes a set of (undirected) equalities and a reduction order and returns a terminating, confluent term rewriting system that has the same expressive power as the set of equalities. As a term rewriting system is also a set of (directed) equalities, these algorithms can also transform a non-confluent TRS into a confluent one. However, success is not guaranteed; the algorithm may succeed, fail, or fail to terminate. Thus naively running a completion algorithm on a synthesized ruleset will not necessarily prevent these issues and may also be computationally expensive. However, we plan to investigate how a modified completion algorithm could be used to prevent some instances of non-confluence without paying a large penalty in either performance or in ruleset size.

\subsection{Other Future Work}

As discussed earlier, choosing a gradual order that describes a strategy or gradient towards a desired goal form is not a trivial task and requires human intuition. Not only must the order describe a gradient over terms towards the goal, but it also must be a valid reduction order; if it is not, then the resulting term rewriting system (as well as the synthesis pipeline that produces it) is not guaranteed to terminate or indeed to actually make progress towards the goal. Proving an order is a valid reduction order is sometimes tricky, and appealing orders are often found not to be valid reduction orders after closer inspection. Reduction orders over terms often follow a set of patterns, such as orders that reduce the size of expressions or reduce the number of particular operations or variables. It would be a useful contribution to produce a handbook of these reduction order families, as well as strategies for efficiently encoding them into sketches for synthesis.

Finally, we intend to explore other applications of this work to other problems, other types of rewriting, and other domains. For example, synthesizing rewrite rules for a very different rewrite algorithm such as that used in an e-graph might require a very different strategy. Synthesizing a term rewriting system for a decidable theory, such as bitwidths of finite width or linear integer arithmetic, might also look very different.

\section{Conclusion}

Term rewriting systems can serve as a key reasoning component within compilers, allowing them to successfully apply optimizations to produce performant code. In this work, we have demonstrated how formal methods such as verification and synthesis can assist compiler authors in creating and maintaining correct and powerful TRS-based reasoning engines. We have shown how these techniques can be used on existing, handwritten TRSs, by first providing proofs of correctness and termination for an existing system, and then using synthesis to grow the system's reasoning power. Our work has resulted in bug fixes and improvements to the system that have since been merged into the main Halide codebase. We then showed how synthesis can be used to create a new term rewriting system, saving the programming effort of writing rules by hand. We argued that an efficient way to write specifications for a term rewriting system is through a reduction order, and discussed ways of translating the purpose of a TRS into this formalism. We showed that our refined synthesis pipeline is able to produce general and robust TRSs that perform well when compared to TRSs created using other means for the same purpose. Finally, we presented some directions for future work, including the goal of synthesizing a TRS that can act as a drop-in replacement for an existing compiler component written in general-purpose code.


text:
\documentclass{article}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{amssymb}
\usepackage{algorithm}
\usepackage{graphicx}
\usepackage[noend]{algpseudocode}
\usepackage{url}
\usepackage{listings}

\newlength\tindent
\setlength{\tindent}{\parindent}
\setlength{\parindent}{0pt}
\renewcommand{\indent}{\hspace*{\tindent}}

\title{COMP 312 Assignment 9}
\author{Daniel Braithwaite}

\begin{document}
	\pagenumbering{gobble}
	\maketitle
	\newpage
  	\pagenumbering{arabic}

  	\section{Python}
		\subsection{Code}
				\lstinputlisting[language=Python]{Programs/01.py}
		\subsection{Output}
			Lambda:  1.0\\
	Estimate of W: 3.07115102268\\
	Conf int of W: (3.0545014142296387, 3.087800631127156)\\
	Estimate of L: 3.07293392046\\
	Conf int of L: (3.0505709116845887, 3.0952969292354391)\\
\\
Lambda:  1.1\\
	Estimate of W: 3.30291442125\\
	Conf int of W: (3.2837205998350485, 3.322108242669946)\\
	Estimate of L: 3.63566350108\\
	Conf int of L: (3.6075450879612228, 3.6637819142041272)\\
\\
Lambda:  1.2\\
	Estimate of W: 3.62809887072\\
	Conf int of W: (3.6021168652532589, 3.6540808761833237)\\
	Estimate of L: 4.36000156289\\
	Conf int of L: (4.3220178721139915, 4.3979852536574624)\\
\\
Lambda:  1.3\\
	Estimate of W: 4.07158480319\\
	Conf int of W: (4.0340691696071893, 4.1091004367628461)\\
	Estimate of L: 5.30478973571\\
	Conf int of L: (5.2468732130936129, 5.3627062583358907)\\
\\
Lambda:  1.4\\
	Estimate of W: 4.69411430471\\
	Conf int of W: (4.6332972140664737, 4.7549313953467323)\\
	Estimate of L: 6.57650750236\\
	Conf int of L: (6.4790046360596518, 6.6740103686512215)\\
\\
Lambda:  1.5\\
	Estimate of W: 5.65164280855\\
	Conf int of W: (5.5619584714599943, 5.7413271456394588)\\
	Estimate of L: 8.49841455215\\
	Conf int of L: (8.3475595502920541, 8.6492695539997584)\\
\\
Lambda:  1.6\\
	Estimate of W: 6.87174636875\\
	Conf int of W: (6.7320208108733182, 7.0114719266196799)\\
	Estimate of L: 10.9660720227\\
	Conf int of L: (10.724309048752211, 11.207834996743063)\\

Lambda:  1.7\\
	Estimate of W: 9.80069631029\\
	Conf int of W: (9.4732894506847209, 10.12810316989798)\\
	Estimate of L: 16.6594119506\\
	Conf int of L: (16.06997802383087, 17.248845877432245)\\
\\
Lambda:  1.8\\
	Estimate of W: 16.76146322\\
	Conf int of W: (15.835988850257873, 17.686937589686902)\\
	Estimate of L: 30.1149881389\\
	Conf int of L: (28.392912294383457, 31.837063983325834)\\
\\
Lambda:  1.9\\
	Estimate of W: 51.2452902908\\
	Conf int of W: (45.592299283965971, 56.898281297587403)\\
	Estimate of L: 96.2641781498\\
	Conf int of L: (85.643801313328282, 106.88455498621624)\\
\\
Lambda:  2.0\\
	Estimate of W: 151.156888049\\
	Conf int of W: (140.73420055369161, 161.57957554384637)\\
	Estimate of L: 286.424293786\\
	Conf int of L: (266.92387870816964, 305.92470886372922)\\
\\
Lambda:  2.1\\
	Estimate of W: 265.889338596\\
	Conf int of W: (255.27250813040555, 276.50616906077812)\\
	Estimate of L: 503.940054731\\
	Conf int of L: (484.4649347904483, 523.4151746720172)\\
\\
Lambda:  2.2\\
	Estimate of W: 373.765552379\\
	Conf int of W: (363.61554795108702, 383.91555680696695)\\
	Estimate of L: 710.667473673\\
	Conf int of L: (692.00715517220692, 729.32779217390498)\\

\end{document}

max_stars_repo_name: danielbraithwt/University

max_issues_repo_name: danielbraithwt/University

max_forks_repo_name: danielbraithwt/University