---
title: "Stratigraphic graphs"
author: "Joe Roe"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
bibliography: "`r system.file('REFERENCES.bib', package='stratigraphr')`"
vignette: >
  %\VignetteIndexEntry{Stratigraphic graphs}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,
  message = FALSE
)
```

`stratigraphr` provides a tidy framework for working with archaeological stratigraphy as a graph data structure in R.
This vignette introduces the core functions for constructing, validating, and visualising stratigraphic graphs (a 'Harris matrix'), and demonstrates how they integrate with the broader tidygraph ecosystem.

## Terminology

Archaeological stratigraphy records the depositional history of a site as a set of **units**^[Also called contexts, loci, layers, or other terms depending on regional tradition.] and their relative temporal ordering as determining by the law of superposition [@Harris1979-qs].
A **stratigraphic graph** is the formal representation of this information: a directed acyclic graph where nodes represent units and directed edges represent stratigraphic **relations** [@Dye2015-xd].

The **Harris matrix** is the conventional visual representation of a stratigraphic graph, where units are arranged in boxes and lines indicate direct stratigraphic relations [@Harris1979-qs].
While the terms are sometimes used interchangeably in archaeological literature, it is useful to distinguish the abstract data structure (the stratigraphic graph) from its visual representation (the Harris matrix).

Other information commonly featured on Harris matrices -- such as unit equivalencies, the distinction between deposits and interfaces, unit types, or phase assignments -- are not part of the pure stratigraphic graph.
They are associated information that can be used for display, analysis, or modelling, but the graph structure itself is determined solely by the above/below relations between units.

## stratigraph objects

The `stratigraph` object is the package's central data structure.
It is a specialised subclass of `tidygraph::tbl_graph`, representing a stratigraphic graph as tidy data with nodes (units) and edges (relations).

The package includes the dataset `harris12`, a classic stratigraphic sequence from Harris [-@Harris1979-qs, figure 12]:

```{r harris12-data}
library("stratigraphr")

data("harris12")
harris12
```

The data frame contains a `context` column with unit labels and list-columns (`above`, `below`, `equal`) describing the stratigraphic relations between units.
To construct the graph, we specify which column contains the unit labels and which column describes the relations:

```{r harris12-graph}
h12_graph <- stratigraph(harris12, "context", "above")
h12_graph
```

The `direction` argument (default `"above"`) indicates whether the relation column lists units *above* or *below* each unit.
Because `stratigraph` inherits from `tbl_graph`, it can be manipulated using the full suite of tidygraph and igraph functions:

```{r harris12-tidygraph}
library("tidygraph")

# Inspect the nodes (units)
h12_graph |>
  activate("nodes") |>
  as_tibble()

# Inspect the edges (relations)
h12_graph |>
  activate("edges") |>
  as_tibble()
```

See the [tidygraph documentation](https://tidygraph.data-imaginist.com/) for more details on tidy graph manipulation and analysis.

## Importing stratigraphic data

Stratigraphic data is fundamentally simple: a data frame of units and relations.
The `stratigraph()` function expects:

- A **label column** with unique identifiers for each unit
- A **relation column** describing stratigraphic relations between units

The relation column can be either a list-column (where each element is a vector of related units) or a regular column in long format (where each row represents a single relation).

### Constructing from scratch

You can construct a stratigraphic data frame directly in R:

```{r from-scratch}
strat_data <- data.frame(
  unit = c("A", "B", "C"),
  below = c("B", "C", NA)
)

stratigraph(strat_data, "unit", "below", direction = "below")
```

### Reading from CSV

More commonly, stratigraphic data is stored in a CSV file.
With long-format input, you can read the CSV and pass it directly to `stratigraph()`:

```{r csv-workflow}
library("readr")

csv_text <- "context,above
A,
B,A
C,A"

strat_data <- read_csv(csv_text, show_col_types = FALSE)

stratigraph(strat_data, "context", "above")
```

Stratigraphic data often includes both 'above' and 'below' columns.
These should describe the same graph from opposite perspectives.
Use `strat_is_mirror()` to verify they are consistent:

```{r mirror-check}
strat_is_mirror(harris12$context, harris12$above, harris12$below)
```

### Reading LST files

The `read_lst()` function reads stratigraphic data from LST format files, used by BASP Harris, Stratify, and ArchEd:


```{r read-lst}
lst_file <- system.file("extdata", "bonn.lst", package = "stratigraphr")
lst_data <- read_lst(lst_file)
stratigraph(lst_data, "name", "above")
```

It supports both the original BASP format and the extended format used by Stratify and ArchEd – see `?read_lst` for details.

## Validating stratigraphic graphs

Following @Dye2015-xd, a valid stratigraphic graph^[Or 'archaeological sequence diagram', in their terminology.] must satisfy two properties:

1. **It is acyclic** -- it contains no cycles, i.e. no circular chains of relations.
A cycle would mean unit A is above unit B, B is above C, and C is above A, violating the law of stratigraphic superposition.
2. **It contains no redundant relations** -- a relation is redundant if it is already implied by other relations through transitivity.
Transitivity means that if A is above B and B is above C, then A is necessarily above C; the direct relation between A and C adds no new information.
Per Harris' 'Law of Stratigraphical Succession', only the uppermost and undermost relations are significant when placing a unit in a stratigraphic sequence.

The `stratigraph()` function warns on construction if the graph is invalid, but allows you to create and inspect invalid graphs.
Use `strg_is_valid()` for explicit checking:

```{r validity-check}
strg_is_valid(h12_graph)
```

### Cycles

A cycle occurs when the stratigraphic relations are contradictory.
This violates the law of superposition and indicates an error in the data:

```{r cycle-example}
cycle_data <- data.frame(
  unit = c("A", "B", "C"),
  below = c("B", "C", "A")
)

cycle_graph <- stratigraph(cycle_data, "unit", "below", direction = "below")
```

### Redundant relations

A redundant relation is one that is already implied by transitivity.
For example, if A is above B and B is above C, then A's position above C is already established through the chain A--B--C.
An explicit relation between A and C is redundant and should be removed to produce a valid stratigraphic graph.

```{r redundant-example}
redundant_data <- data.frame(
  unit = c("A", "B", "B", "C"),
  above = c("B", "C", "C", NA)
)

redundant_graph <- stratigraph(redundant_data, "unit", "above")
strg_is_valid(redundant_graph)
```

The `strg_prune()` function computes the transitive reduction, removing redundant edges:

```{r prune-example}
pruned_graph <- strg_prune(redundant_graph)
pruned_graph
strg_is_valid(pruned_graph)
```

## Plotting with ggraph

[ggraph](https://ggraph.data-imaginist.com/) extends the grammar of graphics to network visualisation, providing a ggplot2-compatible interface for plotting graph data.
Because `stratigraph` objects inherit from `tbl_graph`, they work seamlessly with ggraph: you can map aesthetics, add layers, and use geoms just as you would with any other tidy data.

The Harris matrix is the conventional visual representation of a stratigraphic graph.
We can reproduce it using ggraph with a Sugiyama (layered) layout:

```{r ggraph-harris}
library("ggraph")

ggraph(h12_graph, layout = "sugiyama") +
  geom_edge_elbow() +
  geom_node_label(aes(label = context), label.r = unit(0, "mm")) +
  theme_graph()
```

The Sugiyama layout arranges units in horizontal layers according to their stratigraphic position, with the earliest (lowest) units at the bottom.
The `geom_edge_elbow()` geom produces the right-angled edges characteristic of Harris matrices.

### Incorporating associated information

While the stratigraphic graph itself contains only units and relations, we can use associated information to enhance the visualisation.
The `shub1` dataset [@Richter2017-xy] includes columns for context `type`, `phase`, and `structure`, which can be mapped to visual aesthetics:

```{r ggraph-shub1, fig.height = 8}
shub1_graph <- stratigraph(shub1, "context", "above")

ggraph(shub1_graph, layout = "sugiyama") +
  geom_edge_elbow() +
  geom_node_point(aes(shape = type), fill = "white", size = 6) +
  geom_node_text(aes(label = context), vjust = 0.5, size = 3) +
  scale_shape_manual(values = c(22, 21, 23, 24)) +
  theme_graph()
```

## Further reading

- **Chronological Query Language (CQL)**: The `vignette("cql")` vignette describes how to convert stratigraphic graphs into OxCal chronological models for Bayesian radiocarbon calibration.
- **Tidy radiocarbon analysis**: Radiocarbon analysis functions have moved to the [c14 package](https://c14.joeroe.io/).
See its [introductory vignette](https://c14.joeroe.io/articles/c14.html) for details.
- **Graph manipulation**: The [tidygraph documentation](https://tidygraph.data-imaginist.com/) provides comprehensive coverage of graph manipulation in R.

## References