---
title: "powdR: Full Pattern Summation of X-ray Powder Diffraction Data"
bibliography: bibliography.bib
author: "Benjamin Butler"
date: "`r Sys.Date()`"
output:
      rmarkdown::html_vignette:
        number_sections: true
        fig_caption: yes
        fig_width: 7
        fig_height: 6
vignette: >
  %\VignetteIndexEntry{powdR: Full Pattern Summation of X-ray Powder Diffraction Data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<style type="text/css">

body{ /* Normal  */
      font-size: 16px;
  }
code.r{ /* Code block */
    font-size: 16px;
}
p {line-height: 1.5em;}
</style>

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

# Overview
`powdR` is an R implementation of the full pattern summation approach to quantitative mineralogy from X-ray powder diffraction (XRPD) data. Whilst this is available in Excel spreadsheets such as FULLPAT [@Chipera2002] and RockJock [@Eberl2003], implementation in R allows for faster computation than is currently available, and provides a user-friendly `Shiny` application to help with the often iterative process of mineral quantification. Furthermore, the `afps()` function in `powdR` is designed to automate the full pattern summation procedure, which is particularly advantageous in high-throughput XRPD datasets.


# Full pattern summation
A powerful property of XRPD data is that it can provide quantitative estimates of phase concentrations in mineral mixtures. Though several methods are available for such quantitative analysis, full pattern summation (also referred to as "full pattern fitting of prior measured standards") is particularly suitable for mixtures containing crystalline mineral components in combination with disordered and/or X-ray amorphous phases. Soil is a prime example of such mixtures, where crystalline minerals such as quartz and feldspars can be present in combination with clay minerals (i.e. disordered phases), and organic matter (i.e. amorphous phases).

The full pattern summation implemented in `powdR` is based on the principle that an observed pattern is comprised of the sum of signals from the individual components within it. A key component of this approach is a library containing measured or calculated patterns of the pure phases that may be encountered within the samples. These "reference" patterns would ideally be measured on the same instrument as the sample. To quantify a given sample, suitable phases from the library are selected that together account for all peaks within the data, and their relative contributions to the observed signal optimised until an appropriate fit is achieved [@Chipera2002]. This fit is usually refined using least squares optimisation of an objective parameter. The scaled intensities of the optimised patterns are then converted to weight % using reference intensity ratios, which are a measure of diffracting power relative to a standard mineral (usually corundum).  

---

# Using `powdR`

## The `powdRlib` object
A key component of the functions within `powdR` is the library of reference patterns. These are stored within a `powdRlib` object created using the `powdRlib()` constructor function. `powdRlib()` builds a `powdRlib` object from two components. The first component, specified via the `xrd_table` argument of `powdRlib()`, is a data frame of the count intensities of the reference patterns, with their 2$\theta$ axis as the first column. The column for a given reference pattern must be named using a unique identifier (a phase ID). An example of such a format is provided in the `minerals_xrd` data:

```{r, message = FALSE, warning = FALSE}
library(powdR)

data(minerals_xrd)

head(minerals_xrd[1:9])
```
The second component required to build a `powdRlib` object (specified via the `phases_table` argument of `powdRlib()`) is a data frame containing 3 columns. The first is a string of unique ID's representing each reference pattern in the data provided to the `xrd_table` argument. The second column is the name of the phase group that this reference pattern belongs to (e.g. quartz or dolomite). The third column is the reference intensity ratio (RIR) of that reference pattern (relative to a known standard, usually corundum). An example of the format required for the `phases_table` argument of `powRlib()` is provided in the `minerals_phases` data.

```{r, message = FALSE, warning = FALSE}
data(minerals_phases)

minerals_phases[1:8,]
```

Crucially when building the `powdRlib` object, all phase ID's in the first column of the `phases_table` must match the column names of the `xrd_table` (except the name of the first column which is the 2$\theta$ scale), for example.

```{r, message = FALSE, warning = FALSE}
identical(names(minerals_xrd[-1]),
          minerals_phases$phase_id)
```

Once created, `powdRlib` objects can easily be visualised using `plot()`, which when used for `powdRlib` objects accepts the arguments `wavelength` and `refs` that are used to specify the X-ray wavelength and the reference patterns to plot, respectively. **In all cases where `plot()` is used hereafter, the addition of `interactive = TRUE` to the function call will produce an interactive html graph that can be viewed in RStudio or web browsers.**.

```{r, message = FALSE, warning = FALSE}
my_lib <- powdRlib(minerals_xrd, minerals_phases)

plot(my_lib, wavelength = "Cu",
     refs = c("ALB", "DOL.1",
              "QUA.1", "GOE.2"))
```

### Subsetting a `powdRlib` object
Occasionally it may be useful to subset a reference library to a smaller selection. This can be achieved using `subset()`, which for `powdRlib` objects accepts three arguments; `x`, `refs` and `mode`. The `x` argument specifies the `powdRlib` object to be subset, `refs` specifies the ID's and/or names of phases to select, and `mode` specifies whether these phases are kept (`mode = "keep"`) or removed (`mode = "remove"`).

```{r, message = FALSE, warning = FALSE}
data(rockjock)

#Have a look at the phase ID's in rockjock
rockjock$phases$phase_id[1:10]

#Remove reference patterns from rockjock
rockjock_1 <- subset(rockjock,
                     refs = c("ALUNITE", #phase ID
                              "AMPHIBOLE", #phase ID
                              "ANALCIME", #phase ID
                              "Plagioclase"), #phase name
                     mode = "remove")

#Check number of reference patterns remaining in library
nrow(rockjock_1$phases)

#Keep certain reference patterns of rockjock
rockjock_2 <- subset(rockjock,
                     refs = c("ALUNITE", #phase ID
                              "AMPHIBOLE", #phase ID
                              "ANALCIME", #phase ID
                              "Plagioclase"), #phase name
                     mode = "keep")

#Check number of reference patterns remaining
nrow(rockjock_2$phases)
```

## RockJock
There are two `powdRlib` objects provided as part of the `powdR` package. The first is `minerals` (accessed via `data(minerals)`), which is a simple and low resolution library designed to facilitate fast computation of basic examples. The second is `rockjock` (accessed via `data(rockjock)`), which is a comprehensive library of 169 reference patterns covering most phases that might be encountered in geological and soil samples. The `rockjock` library in `powdR` uses data from the original RockJock program [@Eberl2003] thanks to the permission of Dennis Eberl. In `rockjock`, each reference pattern from the original RockJock program has been scaled to a maximum intensity of 10000 counts, and the RIR's normalised relative to Corundum. All `rockjock` data were analysed using Cu K$\alpha$ radiation.


## RockJock synthetic mixtures
To accompany the `rockjock` reference library, a list of eight synthetic mixtures from the original RockJock program [@Eberl2003] are also included in `powdR` in the `rockjock_mixtures` data (accessed via `data(rockjock_mixtures)`. Their known weights (see `?rockjock_mixtures`) can be compared to full pattern summation outputs (i.e. from `fps()` and `afps()`, detailed below) to assess accuracy.


## Full pattern summation:
Full pattern summation in `powdR` is provided via `fps()`, and an automated version provided in `afps()`. Here the `rockjock` and `rockjock_mixtures` data will be used to demonstrate the use of these functions.

### Full pattern summation **with** internal standard

In some cases samples are prepared for XRPD with an internal standard of known concentration. If this is the case, then the `std_conc` argument of `fps()` and `afps()` can be used to define the concentration of the internal standard (in weight %), which is used in combination with the reference intensity ratios to compute phase concentrations. For example, all samples in the `rockjock_mixtures` data were prepared with 20 % corundum as the internal standard, thus this can be specified using using `std = "CORUNDUM"` and `std_conc = 20` in the call to `fps()` or `afps()`.

```{r, message = FALSE, warning = FALSE}
data("rockjock_mixtures")

fit_1 <- fps(lib = rockjock,
             smpl = rockjock_mixtures$Mix1,
             refs = c("ORDERED_MICROCLINE",
                      "LABRADORITE",
                      "KAOLINITE_DRY_BRANCH",
                      "MONTMORILLONITE_WYO",
                      "ILLITE_1M_RM30",
                      "CORUNDUM"),
             std = "CORUNDUM",
             std_conc = 20,
             align = 0.3)
```

The computed concentrations can be accessed in the `phases` dataframe of the `powdRlib` output:
```{r, message = FALSE, warning = FALSE}
fit_1$phases
```

Further, notice that if the concentration of the internal standard is specified then the phase
concentrations do not necessarily sum to 100 %:
```{r, message = FALSE, warning = FALSE}
sum(fit_1$phases$phase_percent)
```

Unlike other software where only certain phases can be used as an internal standard, any phase can be defined in `powdR`. For example, the `rockjock_mixtures$Mix5` sample contains 20 % quartz (see `?rockjock_mixtures`), thus adding `"QUARTZ"` as the `std` argument results in this reference pattern becoming the internal standard instead.

```{r, message = FALSE, warning = FALSE}
fit_2 <- fps(lib = rockjock,
             smpl = rockjock_mixtures$Mix5,
             refs = c("ORDERED_MICROCLINE",
                      "LABRADORITE",
                      "KAOLINITE_DRY_BRANCH",
                      "MONTMORILLONITE_WYO",
                      "CORUNDUM",
                      "QUARTZ"),
             std = "QUARTZ",
             std_conc = 20,
             align = 0.3)

fit_2$phases

sum(fit_2$phases$phase_percent)
```
### Full pattern summation **without** internal standard

In cases where an internal standard is not added to a sample, phase quantification can be achieved by assuming that all detectable phases can be identified and that they sum to 100 weight %. By setting the `std_conc` argument of `fps()` or `afps()` to `NA`, or leaving it out of the function call, it will be assumed that the sample has been prepared without an internal standard and the phase concentrations computed accordingly.

```{r, message = FALSE, warning = FALSE}
fit_3 <- fps(lib = rockjock,
             smpl = rockjock_mixtures$Mix1,
             refs = c("ORDERED_MICROCLINE",
                      "LABRADORITE",
                      "KAOLINITE_DRY_BRANCH",
                      "MONTMORILLONITE_WYO",
                      "ILLITE_1M_RM30",
                      "CORUNDUM"),
             std_conc = NA,
             std = "CORUNDUM",
             align = 0.3)
```

In this case the phase specified in the `std` argument is only used for sample alignment, and is included in the computed phase concentrations.

```{r, message = FALSE, warning = FALSE}
fit_3$phases
```

Furthermore, the phase concentrations sum to 100 %.
```{r, message = FALSE, warning = FALSE}
sum(fit_3$phases$phase_percent)
```

### Full pattern summation with data harmonisation

It is usually recommended that the reference library used for full pattern summation is measured on the same instrument as the sample using an identical 2$\theta$ range and resolution. In some cases this is not feasible, and the reference library patterns may be from a different instrument to the sample. To allow for seamless use of samples and libraries from different instruments (measured on the same wavelength), `fps()` and `afps()` contain a logical `harmonise` argument (default `= TRUE`). When the sample and library contain non-identical 2$\theta$ axes, `harmonise = TRUE` will convert the data onto the same scale by determining the overlapping 2$\theta$ range and interpolating to the coarsest resolution available.

```{r, message = FALSE, warning = FALSE}
#Create a sample with a shorter 2theta axis than the library
Mix1_short <- subset(rockjock_mixtures$Mix1, tth > 10 & tth < 55)

#Reduce the resolution by selecting only odd rows of the data
Mix1_short <- Mix1_short[seq(1, nrow(Mix1_short), 2),]

fit_4 <- fps(lib = rockjock,
             smpl = Mix1_short,
             refs = c("ORDERED_MICROCLINE",
                      "LABRADORITE",
                      "KAOLINITE_DRY_BRANCH",
                      "MONTMORILLONITE_WYO",
                      "ILLITE_1M_RM30",
                      "CORUNDUM"),
             std = "CORUNDUM",
             align = 0.3)

fit_4$phases
```

## Automated full pattern summation
The selection of suitable reference patterns for full pattern summation can often be challenging and time consuming. An attempt to automate this process is provided in `afps()`, which can select appropriate reference patterns from a reference library and subsequently exclude reference patterns based on limit of detection estimates. Such an approach is considered particularly advantageous when quantifying high-throughput XRPD datasets that display considerable mineralogical variation.

Here the `rockjock` library, containing 169 reference patterns, will be used to quantify one of the samples in the `rockjock_mixtures` data.

```{r, message = FALSE, warning = FALSE}
fit_5 <- afps(lib = rockjock,
              smpl = rockjock_mixtures$Mix1,
              std = "CORUNDUM",
              align = 0.3,
              lod = 1)

fit_5$phases_grouped
```

## Plotting
Plotting results from `fps()` or `afps()` (`powdRfps` and `powdRafps` objects, respectively) is achieved using `plot()`. Static `ggplot()` plots can be created using:
```{r, message = FALSE, warning = FALSE}
plot(fit_5, wavelength = "Cu")
```

Alternatively, interactive `ggplotly()` plots can be created by adding `interactive = TRUE` to the function call, e.g. `plot(fit_5, wavelength = "Cu", interactive = TRUE)`.

In addition to above, plotting for `powdRfps` and `powdRafps` objects can be further adjusted by the `mode` and `xlim` arguments. The `mode` argument can be one of "fit" (the default), "residuals" or "both", for example:
```{r, message = FALSE, warning = FALSE}
plot(fit_5, wavelength = "Cu", mode = "residuals")
```

or alternatively:
```{r, message = FALSE, warning = FALSE}
plot(fit_5, wavelength = "Cu", mode = "both", xlim = c(20,30))
```

## Quantifying multiple samples
The easiest way to quantify multiple samples is with `lapply()`

```{r, message = FALSE, warning = FALSE}
multi_fit <- lapply(rockjock_mixtures[1:3], fps,
                    lib = rockjock,
                    std = "CORUNDUM",
                    refs = c("ORDERED_MICROCLINE",
                             "LABRADORITE",
                             "KAOLINITE_DRY_BRANCH",
                             "MONTMORILLONITE_WYO",
                             "ILLITE_1M_RM30",
                             "CORUNDUM",
                             "QUARTZ"),
                    align = 0.3)

names(multi_fit)
```

## Summarising mineralogy
When multiple samples are quantified, it is often useful to report the phase concentrations of all of the samples in a single table. For a given list of `powdRfps` and/or `powdRafps` objects, the `summarise_mineralogy()` function yields such summary tables, for example:

```{r, message = FALSE, warning = FALSE}
summarise_mineralogy(multi_fit, type = "grouped", order = TRUE)
```

where `type = "grouped"` denotes that phases with the same `phase_name` will be summed together, and `order = TRUE` specifies that the columns will be ordered from most common to least common (assessed by the sum of each phase across the samples).

## The powdR Shiny app

To run `powdR` via the Shiny app, use `run_powdR()`. This loads the application in your default web browser. The application has six tabs:

1. **Reference Library Builder:** Allows you to create and export a `powdRlib` reference library from two .csv files: one for the XRPD measurements, and the other for the ID, name and reference intensity ratio of each pattern.
2. **Reference Library Viewer:** Facilitates quick inspection of the phases within a `powdRlib` reference library.
3. **Reference Library Editor:** Allows the user to easily subset a `powdRlib` reference library .
4. **Full Pattern Summation:** A user friendly interface for iterative full pattern summation of single samples using `fps()` or `afps()`.
5. **Results Viewer/Editor:** Allows for results from previously saved `powdRfps` and `powdRafps` objects to be viewed and edited via addition or removal of reference patterns.
6. **Help** Provides a series of video tutorials (via YouTube) detailing the use of the `powdR` Shiny application.

**Note: that the Shiny application may not work properly in Microsoft Edge web browsers. Please get in touch if you experience any other issues with the Shiny application.**

## References
