---
title: "Chronological Query Language (CQL)"
output: rmarkdown::html_vignette
bibliography: references.bib
vignette: >
  %\VignetteIndexEntry{Chronological Query Language (CQL)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

The Chronological Query Language (CQL) is a tool for formally describing chronological models [@Bronk_Ramsey1998-zw].
It is most commonly used to input data for Bayesian radiocarbon calibration in [OxCal](https://c14.arch.ox.ac.uk/oxcal.html) [@Bronk_Ramsey2009-yk].
`stratigraphr` includes an R interface for CQL2, the version used in OxCal v3+.

This vignette describes how to use this interface to generate CQL models in R.

## Overview

* Individual CQL commands (see <https://c14.arch.ox.ac.uk/oxcalhelp/hlp_commands.html>) are implemented as `cql_*` functions.
* Use `cql()` to group together CQL functions or include arbitrary CQL code.
* Execute CQL scripts in OxCal using `write_oxcal()` or the [oxcAAR package](https://CRAN.R-project.org/package=oxcAAR).


## Writing CQL models

```{r setup}
library("stratigraphr")
```

## Generating CQL models from other types of data

Used in the simple way above, `stratigraphr`'s CQL interface offers little benefit over writing CQL directly.
Its real power is in combining `cql()` with other R tools to build models based on other data.

For the following examples, we will use a simple example dataset of radiocarbon dates:

```{r example-data}
library("purrr")
library("dplyr")

dates <- data.frame(
  lab_id = c("ABC-001", "ABC-002", "ABC-003", "ABC-004", "ABC-005"),
  age = c(5050, 5000, 4950, 4900, 4850),
  error = c(30, 30, 30, 30, 30),
  phase = c("Phase 1", "Phase 1", "Phase 1", "Phase 2", "Phase 2")
)
```

### Tabular data

With radiocarbon data in a tabular format, you can take advantage of `dplyr`'s powerful tools for data manipulation to build CQL models programmatically.

For example, use `dplyr::mutate()` to concisely express a table of dates as CQL `R_Date` commands:

```{r r-date}
dates |>
  mutate(r_date = cql_r_date(lab_id, age, error)) |>
  pluck("r_date") |>
  cql()
```

Or use `dplyr::group_by()` and `dplyr::summarise()` to build phase models.
This is a three stage process:

1. Preprocess the dates, e.g. to order phases.
2. Group dates by phase and summarise them with `cql_phase()`.
3. Group phases into a sequence, using `boundaries` to automatically add boundary constraints between them.

```{r phase-model}
dates |>
  group_by(phase) |>
  summarise(cql = cql_phase(phase, cql_r_date(lab_id, age, error))) |>
  arrange(desc(phase)) |>
  summarise(cql = cql_sequence("Example Sequence", cql, boundaries = TRUE)) |>
  pluck("cql") |>
  cql() ->
  example_cql

example_cql
```

### Stratigraphs

...

## Running CQL models

You can run models generated by `cql()` using the desktop or online versions of OxCal by simply copying the output into the program.
Alternatively, use `write_oxcal()` to create a .oxcal file:

```{r write-oxcal}
oxcal_cql <- cql(
  cql_r_date("ABC-001", 9100, 30),
  cql_r_date("ABC-002", 9200, 30),
  cql_r_date("ABC-003", 9300, 30)
)

write_oxcal(oxcal_cql, "cql.oxcal")
```
```{r write-oxcal-cleanup, include=FALSE}
file.remove("cql.oxcal")
```

You can also run OxCal directly through R using the [oxcAAR](https://CRAN.R-project.org/package=oxcAAR) package.
This depends on a local installation of OxCal.
If you already have one installed, you can set the path to the executable using `oxcAAR::setOxcalExecutablePath()`.
Otherwise, use `oxcAAR::quickSetupOxcal()` to download one, for example to a temporary directory:

```{r oxcaar-setup, results=FALSE, message=FALSE}
library("oxcAAR")
quickSetupOxcal(path = tempdir())
```

You can then use `oxcAAR::executeOxcalScript()` to run the CQL script and `oxcAAR::readOxcalOutput()` to read the output back into R.

```{r oxcaar-execute, results=FALSE, message=FALSE}
executeOxcalScript(oxcal_cql) |>
  readOxcalOutput() ->
  oxcal_output
```

You can parse the output with `oxcAAR::parseOxcalOutput()` and visualise it using `oxcAAR`'s built-in plotting functions:

```{r oxcaar-plot, fig.show="hold"}
oxcal_parsed <- oxcAAR::parseOxcalOutput(oxcal_output)
plot(oxcal_parsed)
calcurve_plot(oxcal_parsed)
```

The current CRAN version of oxcAAR (v. 1.0.0) does not read the posterior probabilities produced by a model with Bayesian calibration, so to work with these you need to install the latest development version (`devtools::install_github("ISAAKiel/oxcAAR")`).
With this, `oxcAAR::parseOxcalOutput()` also contains the modelled results in `$posterior_sigma_ranges` and `$posterior_probabilities`.
Again, you can quickly visualise these with the built-in plotting functions:

```{r oxcaar-bayes, eval=FALSE}
# Not run: slow
example_oxcal <- executeOxcalScript(example_cql)
readOxcalOutput(example_oxcal) |>
  parseOxcalOutput() |>
  plot()
```

## References