cl_machinelearning.qmd

---
title: "Machine learning algorithms for image time series"
format: html
--- 

<a href="https://www.kaggle.com/code/esensing/machine-learning-algorithms-for-image-time-series" target="_blank">
<img src="https://kaggle.com/static/images/open-in-kaggle.svg"/>
</a>

### Configurations to run this chapter{-}

:::{.panel-tabset}
## R
```{r}
#| echo: true
#| eval: true
#| output: false
# load package "tibble"
library(tibble)
# load packages "sits" and "sitsdata"
library(sits)
library(sitsdata)
# set tempdir if it does not exist 
tempdir_r <- "~/sitsbook/tempdir/R/cl_machinelearning"
dir.create(tempdir_r, showWarnings = FALSE)
```
## Python
```{python}
#| echo: true
#| eval: true
#| output: false
# load "pysits" library
from pysits import *
from pathlib import Path
# set tempdir if it does not exist 
tempdir_py = Path.home() / "sitsbook/tempdir/Python/cl_machinelearning"
tempdir_py.mkdir(parents=True, exist_ok=True)
```
:::

## Data used in this chapter

The following examples show how to train machine learning methods and apply them to classify a single time series. We use the set `samples_matogrosso_mod13q1`, containing time series samples from the Brazilian Mato Grosso state obtained from the MODIS MOD13Q1 product. It has 1,892 samples and nine classes (`Cerrado`, `Forest`, `Pasture`, `Soy_Corn`, `Soy_Cotton`, `Soy_Fallow`, `Soy_Millet`). Each time series covers 12 months (23 data points) with six bands (NDVI, EVI, BLUE, RED, NIR, MIR). The samples are arranged along an agricultural year, starting in September and ending in August. The dataset was used in the paper "Big Earth observation time series analysis for monitoring Brazilian agriculture" [@Picoli2018], and is available in the R package `sitsdata`.

## Common interface to machine learning and deep learning models

The `sits_train()` function provides a standard interface to all machine learning models. This function takes two mandatory parameters: the training data (`samples`) and the ML algorithm (`ml_method`). After the model is estimated, it can classify individual time series or data cubes with `sits_classify()`. In what follows, we show how to apply each method to classify a single time series. Then, in Chapter [Image classification in data cubes](https://e-sensing.github.io/sitsbook/cl_rasterclassification.html), we discuss how to classify data cubes.

Since `sits` is aimed at remote sensing users who are not machine learning experts, it provides a set of default values for all classification models. These settings have been chosen based on testing by the authors. Nevertheless, users can control all parameters for each model. Novice users can rely on the default values, while experienced ones can fine-tune model parameters to meet their needs. Model tuning is discussed at the end of this Chapter. 

When a set of time series organized as tibble is taken as input to the classifier, the result is the same tibble with one additional column (`predicted`), which contains the information on the labels assigned for each interval. The results can be shown in text format using the function `sits_show_prediction()` or graphically using `plot()`.

## Random Forest

Random Forest is a machine learning algorithm that uses an ensemble learning method for classification tasks. The algorithm consists of multiple decision trees, each trained on a different subset of the training data and with a different subset of features. To make a prediction, each decision tree in the forest independently classifies the input data. The final prediction is made based on the majority vote of all the decision trees. The randomness in the algorithm comes from the random subsets of data and features used to train each decision tree, which helps to reduce overfitting and improve the accuracy of the model. This classifier measures the importance of each feature in the classification task, which can be helpful in feature selection and data visualization. For an in-depth discussion of the robustness of Random Forest for satellite image time series classification, please see Pelletier et al [@Pelletier2016]. 

```{r}
#| label: fig-ml-rf
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: "Random forest algorithm (source: Venkata Jagannath in Wikipedia)."
knitr::include_graphics("./images/random_forest.png") 
```

`sits` provides `sits_rfor()`, which uses the R `randomForest` package [@Wright2017]; its main parameter is `num_trees`, which is the number of trees to grow with a default value of 100. The model can be visualized using `plot()`.

:::{.panel-tabset}
## R
```{r}
#| eval: false
set.seed(290356)

# Train the Mato Grosso samples with Random Forest algorithm
rfor_model <- sits_train(
    samples = samples_matogrosso_mod13q1, 
    ml_method = sits_rfor(num_trees = 100))

# Plot the most important variables of the model
plot(rfor_model)
```
## Python
```{python}
#| eval: false
r_set_seed(290356)

# Load samples
samples_matogrosso_mod13q1 = load_samples(
    name = "samples_matogrosso_mod13q1", 
    package = "sitsdata"
)

# Train the Mato Grosso samples with Random Forest algorithm
rfor_model = sits_train(
    samples = samples_matogrosso_mod13q1, 
    ml_method = sits_rfor(num_trees = 100))

# Plot the most important variables of the model
plot(rfor_model)
```
:::

```{r}
#| label: fig-ml-rfimp
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: "Most important variables in Random Forest model."
knitr::include_graphics("./images/mlrformodel.png") 
```

The model plot shows the most important explanatory variables, which are the NIR (near infrared) band on date 17 (2007-05-25) and the MIR (middle infrared) band on date 22 (2007-08-13). The NIR value at the end of May captures the growth of the second crop for double cropping classes.  Values of the MIR band at the end of the period (late July to late August) capture bare soil signatures to distinguish between agricultural and natural classes. This corresponds to summertime when the ground is drier after harvesting crops.

:::{.panel-tabset}
## R
```{r}
#| eval: false
# Classify using Random Forest model and plot the result
point_class <- sits_classify(
    data = point_mt_mod13q1, 
    ml_model  = rfor_model)

plot(point_class, bands = c("NDVI", "EVI"))
```
## Python
```{python}
#| eval: false
# Load samples
point_mt_mod13q1 = load_samples(
    name = "point_mt_mod13q1", 
    package = "sitsdata"
)

# Classify using Random Forest model and plot the result
point_class = sits_classify(
    data = point_mt_mod13q1, 
    ml_model  = rfor_model)

plot(point_class, bands = ("NDVI", "EVI"))
```
:::

```{r}
#| label: fig-ml-rfplot
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: "Classification of time series using Random Forest."
knitr::include_graphics("./images/mlrforplot.png") 
```


The result shows that the area started as a forest in 2000, was deforested from 2004 to 2005, used as pasture from 2006 to 2007, and for double-cropping agriculture from 2009 onwards. This behavior is consistent with expert evaluation of land change process in this region of Amazonia.

Random Forest is robust to outliers and can deal with irrelevant inputs [@Hastie2009]. The method tends to overemphasize some variables because its performance tends to stabilize after part of the trees is grown [@Hastie2009]. In cases where abrupt change occurs, such as deforestation mapping, Random Forest (if properly trained) will emphasize the temporal instances and bands that capture such quick change. Before using Random Forest, it is recommended that users balance their training samples as explained in ["Reducing imbalances in training samples"](https://e-sensing.github.io/sitsbook/cl_rasterclassification.html)

## Support vector machine

The support vector machine (SVM) classifier is a generalization of a linear classifier that finds an optimal separation hyperplane which minimizes misclassification [@Cortes1995]. Since a set of samples with $n$ features defines an n-dimensional feature space, hyperplanes are linear ${(n-1)}$-dimensional boundaries that define linear partitions in that space. If the classes are linearly separable on the feature space, there will be an optimal solution defined by the maximal margin hyperplane, which is the separating hyperplane that is farthest from the training observations [@James2013]. The maximal margin is computed as the smallest distance from the observations to the hyperplane. The solution for the hyperplane coefficients depends only on the samples that define the maximum margin criteria, the so-called support vectors.

```{r}
#| label: fig-ml-svm
#| echo: false
#| out.width: "50%"
#| fig.align: "center"
#| fig.cap: "Maximum-margin hyperplane and margins for an SVM trained with samples from two classes. Samples on the margin are called the support vectors. (source: Larhmam in Wikipedia)."
knitr::include_graphics("./images/svm_margin.png") 
```

For data that is not linearly separable, SVM provides kernel functions that map the original feature space into a higher dimensional space. Despite having a linear boundary on the enlarged feature space, the new classification model generally translates its hyperplane to a nonlinear boundary in the original attribute space. Kernels are an efficient computational strategy to produce nonlinear boundaries in the input attribute space; thus, they improve training-class separation. SVM is one of the most widely used algorithms in machine learning applications and has been applied to classify remote sensing data [@Mountrakis2011].

In `sits`, SVM is implemented as a wrapper of `e1071` R package that uses the `LIBSVM` implementation [@Chang2011]. The `sits` package adopts the *one-against-one* method for multiclass classification. For a $q$ class problem, this method creates ${q(q-1)/2}$ SVM binary models, one for each class pair combination, testing any unknown input vectors throughout all those models. A voting scheme computes the overall result.

The example below shows how to apply SVM to classify time series using default values. The main parameters are `kernel`, which controls whether to use a nonlinear transformation (default is `radial`), `cost`, which measures the punishment for wrongly-classified samples (default is 10), and `cross`, which sets the value of the k-fold cross validation (default is 10).

:::{.panel-tabset}
## R
```{r}
#| eval: false
set.seed(290356)

# Train an SVM model
svm_model <- sits_train(
    samples = samples_matogrosso_mod13q1, 
    ml_method = sits_svm())

# Classify using the SVM model and plot the result
point_class <- sits_classify(
    data = point_mt_mod13q1, 
    ml_model = svm_model)

# Plot the result
plot(point_class, bands = c("NDVI", "EVI"))
```
## Python
```{python}
#| eval: false
r_set_seed(290356)

# Train an SVM model
svm_model = sits_train(
    samples = samples_matogrosso_mod13q1, 
    ml_method = sits_svm())

# Classify using the SVM model and plot the result
point_class = sits_classify(
    data = point_mt_mod13q1, 
    ml_model = svm_model)

# Plot the result
plot(point_class, bands = ("NDVI", "EVI"))
```
:::

```{r}
#| label: fig-ml-svmplot
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: "Classification of time series using SVM."
knitr::include_graphics("./images/mlsvmplot.png") 
```

The SVM classifier is less stable and less robust to outliers than the Random Forest method. In this example, it tends to misclassify some of the data. In 2008, it is likely that the correct land class was still `Pasture` rather than `Soy_Millet` as produced by the algorithm, while the `Soy_Cotton` class in 2012 is also inconsistent with the previous and latter classification of `Soy_Corn`.

## Extreme gradient boosting

XGBoost (eXtreme Gradient Boosting) [@Chen2016] is an implementation of gradient boosted decision trees designed for speed and performance.  It is an ensemble learning method, meaning it combines the predictions from multiple models to produce a final prediction. XGBoost builds trees one at a time, where each new tree helps to correct errors made by previously trained tree. Each tree builds a new model to correct the errors made by previous models. Using gradient descent, the algorithm iteratively adjusts the predictions of each tree by focusing on instances where previous trees made errors. Models are added sequentially until no further improvements can be made. 

Although Random Forest and boosting use trees for classification, there are significant differences. While Random Forest builds multiple decision trees in parallel and merges them together later, XGBoost builds trees one at a time. In XGBoost, each new tree helps to correct errors made by previously trained tree. XGBoost is often preferred for its speed and performance, particularly on large datasets and is well-suited for problems where precision is paramount. Random Forest, on the other hand, is simpler to implement, more interpretable, and can be more robust to overfitting, making it a good choice for general-purpose applications.

```{r}
#| label: fig-ml-xgb
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: "Flow chart of XGBoost algorithm (source: Guo et al., Applied Sciences, 2020)." 
knitr::include_graphics("./images/flow_chart_xgboost.png") 
```

The boosting method starts from a weak predictor and then improves performance sequentially by fitting a better model at each iteration. It fits a simple classifier to the training data and uses the residuals of the fit to build a predictor. Typically, the base classifier is a regression tree. Although random forest and boosting use trees for classification, there are significant differences. The performance of Random Forest generally increases with the number of trees until it becomes stable. Boosting trees apply finer divisions over previous results to improve performance [@Hastie2009]. Some recent papers show that it outperforms Random Forest for remote sensing image classification [@Jafarzadeh2021]. However, this result is not generalizable since the quality of the training dataset controls actual performance.

In `sits`, the XGBoost method is implemented by the `sits_xbgoost()` function, based on `XGBoost` R package, and has five hyperparameters that require tuning. The `sits_xbgoost()` function takes the user choices as input to a cross-validation to determine suitable values for the predictor.

The learning rate `eta` varies from 0.0 to 1.0 and should be kept small (default is 0.3) to avoid overfitting. The minimum loss value `gamma` specifies the minimum reduction required to make a split. Its default is 0; increasing it makes the algorithm more conservative. The `max_depth` value controls the maximum depth of the trees. Increasing this value will make the model more complex and likely to overfit (default is 6). The `subsample` parameter controls the percentage of samples supplied to a tree. Its default is 1 (maximum). Setting it to lower values means that xgboost randomly collects only part of the data instances to grow trees, thus preventing overfitting. The `nrounds` parameter controls the maximum number of boosting interactions; its default is 100, which has proven to be enough in most cases. To follow the convergence of the algorithm, users can turn the `verbose` parameter on. In general, the results using the extreme gradient boosting algorithm are similar to the Random Forest method.

:::{.panel-tabset}
## R
```{r}
#| eval: false
set.seed(290356)

# Train using  XGBoost
xgb_model <- sits_train(
    samples = samples_matogrosso_mod13q1, 
    ml_method = sits_xgboost(verbose = FALSE))

# Classify using SVM model and plot the result
point_class_xgb <- sits_classify(
    data = point_mt_mod13q1, 
    ml_model = xgb_model)

# View classification
plot(point_class_xgb, bands = c("NDVI", "EVI"))
```
## Python
```{python}
#| eval: false
r_set_seed(290356)

# Train using  XGBoost
xgb_model = sits_train(
    samples = samples_matogrosso_mod13q1, 
    ml_method = sits_xgboost(verbose = False))

# Classify using SVM model and plot the result
point_class_xgb = sits_classify(
    data = point_mt_mod13q1, 
    ml_model = xgb_model)

# View classification
plot(point_class_xgb, bands = ("NDVI", "EVI"))
```
:::

```{r}
#| label: fig-ml-xgbplot
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: "Classification of time series using XGBoost."

knitr::include_graphics("./images/mlxgbplot.png") 
```

## Deep learning using multilayer perceptron

To support deep learning methods, `sits` uses the `torch` R package, which takes the Facebook `torch` C++ library as a back-end. Machine learning algorithms that use the R `torch` package are similar to those developed using `PyTorch`. The simplest deep learning method is multilayer perceptron (MLP), which are feedforward artificial neural networks. An MLP consists of three kinds of nodes: an input layer, a set of hidden layers, and an output layer. The input layer has the same dimension as the number of features in the dataset. The hidden layers attempt to approximate the best classification function. The output layer decides which class should be assigned to the input. In an MLP, all inputs are treated equally at first; based on iterative matching of training and test data, the backpropagation technique feeds information back to the initial layers to identify the most suitable combination of inputs that produces the best output.

In `sits`, MLP models can be built using `sits_mlp()`. Since there is no established model for generic classification of satellite image time series, designing MLP models requires parameter customization. The most important decisions are the number of layers in the model and the number of neurons per layer. These values are set by the `layers` parameter, which is a list of integer values. The size of the list is the number of layers, and each element indicates the number of nodes per layer.

The choice of the number of layers depends on the inherent separability of the dataset to be classified. For datasets where the classes have different signatures, a shallow model (with three layers) may provide appropriate responses. More complex situations require models of deeper hierarchy. Models with many hidden layers may take a long time to train and may not converge. We suggest to start with three layers and test different options for the number of neurons per layer before increasing the number of layers. In our experience, using three to five layers is a reasonable compromise if the training data has a good quality. Further increase in the number of layers will not improve the model. 

MLP models also need to include the activation function. The activation function of a node defines the output of that node given an input or set of inputs. Following standard practices [@Goodfellow2016], we use the `relu` activation function.

The optimization method (`optimizer`) represents the gradient descent algorithm to be used. These methods aim to maximize an objective function by updating the parameters in the opposite direction of the gradient of the objective function [@Ruder2016]. Since gradient descent plays a key role in deep learning model fitting, developing optimizers is an important topic of research [@Bottou2018]. Many optimizers have been proposed in the literature, and recent results are reviewed by Schmidt et al. [@Schmidt2021]. The Adamw optimizer  provides a good baseline and reliable performance for general deep learning applications [@Kingma2017]. By default, all deep learning algorithms in `sits` use Adamw. 

Another relevant parameter is the list of dropout rates (`dropout`). Dropout is a technique for randomly dropping units from the neural network during training [@Srivastava2014]. By randomly discarding some neurons, dropout reduces overfitting. Since a cascade of neural nets aims to improve learning as more data is acquired, discarding some neurons may seem like a waste of resources. In practice, dropout prevents an early convergence to a local minimum [@Goodfellow2016]. We suggest users experiment with different dropout rates, starting from small values (10-30%) and increasing as required.

The following example shows how to use `sits_mlp()`. The default parameters have been chosen based on a modified version of @Wang2017, which proposes using multilayer perceptron as a baseline for time series classification. These parameters are: (a) Three layers with 512 neurons each, specified by the parameter `layers`; (b) Using the "relu" activation function; (c) dropout rates of 40%, 30%, and 20% for the layers; (d) the "optimizer_adamw" as optimizer (default value); (e) a number of training steps (`epochs`) of 100; (f) a `batch_size` of 64, which indicates how many time series are used for input at a given step; and (g) a validation percentage of 20%, which means 20% of the samples will be randomly set aside for validation.

To simplify the output, the `verbose` option has been turned off. After the model has been generated, we plot its training history.

:::{.panel-tabset}

## R
```{r}
#| eval: false
set.seed(290356)

# Train using an MLP model
# This is an example of how to set parameters
# First-time users should test default options first
mlp_model <- sits_train(
    samples = samples_matogrosso_mod13q1, 
    ml_method = sits_mlp(
        optimizer        = torch::optim_adamw, 
        layers           = c(512, 512, 512),
        dropout_rates    = c(0.40, 0.30, 0.20),
        epochs           = 80,
        batch_size       = 64,
        verbose          = FALSE,
        validation_split = 0.2))

# Show training evolution
plot(mlp_model)
```

## Python
```{python}
#| eval: false
r_set_seed(290356)

# Train using an MLP model
# This is an example of how to set parameters
# First-time users should test default options first
mlp_model = sits_train(
    samples = samples_matogrosso_mod13q1,
    ml_method = sits_mlp(
        optimizer        = "torch::optim_adamw",
        layers           = (512, 512, 512),
        dropout_rates    = (0.40, 0.30, 0.20),
        epochs           = 80,
        batch_size       = 64,
        verbose          = False,
        validation_split = 0.2))

# Show training evolution
plot(mlp_model)
```
:::

```{r}
#| label: fig-ml-mlp
#| echo: false
#| out.width: "80%"
#| fig.align: "center"
#| fig.cap: "Evolution of training accuracy of MLP model."
knitr::include_graphics("./images/mlmlpmodel.png") 
```

Then, we classify a 16-year time series using the multilayer perceptron model.

:::{.panel-tabset}
## R
```{r}
#| eval: false
# Classify using MLP model and plot the result
point_mt_mod13q1 |>  
    sits_classify(mlp_model) |>  
    plot(bands = c("NDVI", "EVI"))
```
## Python
```{python}
#| eval: false
# Classify using MLP model
point_class = sits_classify(point_mt_mod13q1, mlp_model)

# Plot
plot(
    point_class, bands = ("NDVI", "EVI")
)
```
:::

```{r}
#| label: fig-ml-mlplot
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: "Classification of time series using MLP."
knitr::include_graphics("./images/mlmlpplot.png") 
```

In theory, multilayer perceptron model can capture more subtle changes than Random Forest and XGBoost. In this specific case, the result is similar to theirs. Although the model mixes the `Soy_Corn` and `Soy_Millet` classes, the distinction between their temporal signatures is quite subtle. Also it suggests the need to improve the number of samples. In this example, the MLP model shows an increase in sensitivity compared to previous models. We recommend to compare different configurations since the MLP model is sensitive to changes in its parameters.

## Temporal Convolutional Neural Network (TempCNN)

Convolutional neural networks (CNN) are deep learning methods that apply convolution filters (sliding windows) to the input data sequentially. The Temporal Convolutional Neural Network (TempCNN) is a neural network architecture specifically designed to process sequential data such as time series. In the case of time series, a 1D CNN applies a moving temporal window to the time series to produce another time series as the result of the convolution. 

The TempCNN architecture for satellite image time series classification is proposed by Pelletier et al. [@Pelletier2019].  It has three 1D convolutional layers and a final softmax layer for classification. The authors combine different methods to avoid overfitting and reduce the vanishing gradient effect, including dropout, regularization, and batch normalization. In the TempCNN reference paper [@Pelletier2019], the authors favourably compare their model with the Recurrent Neural Network proposed by Russwurm and Körner [@Russwurm2018]. @fig-ml-tcnn shows the architecture of the TempCNN model. TempCNN applies one-dimensional convolutions on the input sequence to capture temporal dependencies, allowing the network to learn long-term dependencies in the input sequence. Each layer of the model captures temporal dependencies at a different scale. Due to its multi-scale approach, TempCNN can capture complex temporal patterns in the data and produce accurate predictions.

```{r}
#| label: fig-ml-tcnn
#| echo: false
#| out.width: "80%"
#| fig.align: "center"
#| fig.cap: |
#|   Structure of tempCNN architecture (source: [@Pelletier2019]).
knitr::include_graphics("./images/tempcnn.png")
```

The function `sits_tempcnn()` implements the model. The first parameter is the `optimizer` used in the backpropagation phase for gradient descent. The default is `adamw` which is considered as a stable and reliable optimization function. The parameter `cnn_layers` controls the number of 1D-CNN layers and the size of the filters applied at each layer; the default values are three CNNs with 128 units. The parameter `cnn_kernels` indicates the size of the convolution kernels; the default is kernels of size 7. Activation for all 1D-CNN layers uses the "relu" function. The dropout rates for each 1D-CNN layer are controlled individually by the parameter `cnn_dropout_rates`. The `validation_split` controls the size of the test set relative to the full dataset. We recommend setting aside at least 20% of the samples for validation.

:::{.panel-tabset}
## R
```{r}
#| eval: false
set.seed(290356)
library(torch)

# Train using tempCNN
tempcnn_model <- sits_train(
    samples_matogrosso_mod13q1, 
    sits_tempcnn(
        optimizer            = torch::optim_adamw,
        cnn_layers           = c(256, 256, 256),
        cnn_kernels          = c(7, 7, 7),
        cnn_dropout_rates    = c(0.2, 0.2, 0.2),
        epochs               = 80,
        batch_size           = 64,
        validation_split     = 0.2,
        verbose              = TRUE))

# Show training evolution
plot(tempcnn_model)
```
## Python
```{python}
#| eval: false
r_set_seed(290356)

# Train using tempCNN
tempcnn_model = sits_train(
    samples_matogrosso_mod13q1, 
    sits_tempcnn(
        optimizer            = "torch::optim_adamw",
        cnn_layers           = (256, 256, 256),
        cnn_kernels          = (7, 7, 7),
        cnn_dropout_rates    = (0.2, 0.2, 0.2),
        epochs               = 80,
        batch_size           = 64,
        validation_split     = 0.2,
        verbose              = False))

# Show training evolution
plot(tempcnn_model)
```
:::

```{r}
#| label: fig-ml-tcnnmodel
#| echo: false
#| out.width: "80%"
#| fig.align: "center"
#| fig.cap: |
#|    Training evolution of TempCNN model.
knitr::include_graphics("./images/mltcnnmodel.png")
```

Using the TempCNN model, we classify a 16-year time series.

:::{.panel-tabset}
## R
```{r}
#| eval: false
# Classify using TempCNN model and plot the result
point_mt_mod13q1 |>  
    sits_classify(tempcnn_model) |>  
    plot(bands = c("NDVI", "EVI"))
```
## Python
```{python}
#| eval: false
# Classify using TempCNN model
point_class = sits_classify(point_mt_mod13q1, tempcnn_model)

# Plot
plot(
    point_class, bands = ("NDVI", "EVI")
)
```
:::

```{r}
#| label: fig-ml-tcnnplot
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: |
#|   Classification of time series using TempCNN
knitr::include_graphics("./images/mltcnnplot.png")
```


The result has important differences from the previous ones. The TempCNN model indicates the `Soy_Cotton` class as the most likely one in 2004. While this result is possibly wrong, it shows that the time series for 2004 is different from those of Forest and Pasture classes. One possible explanation is that there was forest degradation in 2004, leading to a signature that is a mix of forest and bare soil. In this case, including forest degradation samples could improve the training data. In our experience, TempCNN models are a reliable way of classifying image time series [@Simoes2021]. Recent work which compares different models also provides evidence that TempCNN models have satisfactory behavior, especially in the case of crop classes [@Russwurm2020].

## Residual 1D CNN networks (ResNet)

A residual 1D CNN network, also known as ResNet, is an extension of the standard 1D CNN architecture,  adding residual connections between the layers. Residual connections allow the network to learn residual mappings, which are the difference between the input and output of a layer. By adding these residual connections, the network can learn to bypass specific layers and still capture essential features in the data.

The Residual Network (ResNet) for time series classification was proposed by Wang et al. [@Wang2017], based on the idea of deep residual networks for 2D image recognition [@He2016]. The ResNet architecture comprises 11 layers, with three blocks of three 1D CNN layers each (see @fig-ml-rnet). Each block corresponds to a 1D CNN architecture. The output of each block is combined with a shortcut that links its output to its input, called a skip connection. The purpose of combining the input layer of each block with its output layer (after the convolutions) is to avoid the so-called "vanishing gradient problem". This issue occurs in deep networks as the neural network's weights are updated based on the partial derivative of the error function. If the gradient is too small, the weights will not be updated, stopping the training [@Hochreiter1998]. Skip connections aim to avoid vanishing gradients from occurring, allowing deep networks to be trained.

```{r}
#| label: fig-ml-rnet
#| echo: false
#| out.width: "100%"
#| fig.align: "center"
#| fig.cap: |
#|   Structure of ResNet architecture (source: [@Wang2017]).
knitr::include_graphics("./images/resnet.png")
```


In `sits`, the Residual Network is implemented using `sits_resnet()`. The default parameters are those proposed by Wang et al. [@Wang2017], as implemented by Fawaz et al. [@Fawaz2019]. The first parameter is `blocks`, which controls the number of blocks and the size of filters in each block. By default, the model implements three blocks, the first with 64 filters and the others with 128. The parameter `kernels` controls the size of the kernels of the three layers inside each block. It is useful to experiment a bit with these kernel sizes in the case of satellite image time series. The default activation is "relu", which is recommended in the literature to reduce the problem of vanishing gradients. The default optimizer is `optim_adamw`, available in package `torch`.

:::{.panel-tabset}
## R
```{r}
#| eval: false
# Train using ResNet
resnet_model <- sits_train(samples_matogrosso_mod13q1, 
                       sits_resnet(
                          blocks               = c(64, 128, 128),
                          kernels              = c(7, 5, 3),
                          epochs               = 100,
                          batch_size           = 64,
                          validation_split     = 0.2,
                          verbose              = FALSE))

# Show training evolution
plot(resnet_model)
```
## Python
```{python}
#| eval: false
# Train using ResNet
resnet_model = sits_train(samples_matogrosso_mod13q1, 
                       sits_resnet(
                          blocks               = (64, 128, 128),
                          kernels              = (7, 5, 3),
                          epochs               = 100,
                          batch_size           = 64,
                          validation_split     = 0.2,
                          verbose              = False))

# Show training evolution
plot(resnet_model)
```
:::

```{r}
#| label: fig-ml-resnetmodel
#| echo: false
#| out.width: "80%"
#| fig.align: "center"
#| fig.cap: |
#|    Training evolution of ResNet model.
knitr::include_graphics("./images/mlresnetmodel.png")
```

Using the ResNet model, we classify a 16-year time series.

:::{.panel-tabset}
## R
```{r}
#| eval: false
# Classify using Resnet model and plot the result
point_mt_mod13q1 |>  
    sits_classify(resnet_model) |>  
    plot(bands = c("NDVI", "EVI"))
```
## Python
```{python}
#| eval: false
# Classify using Resnet model and plot the result
point_class = sits_classify(point_mt_mod13q1, resnet_model)

# Plot
plot(
    point_class, bands = ("NDVI", "EVI")
)
```
:::

```{r}
#| label: fig-ml-resnetplot
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: |
#|   Classification of time series using TempCNN
knitr::include_graphics("./images/mlresnetplot.png")
```

In this case, the result of the ResNet model is quite similar to that of the TempCNN model. In the same way as TempCNN, the ResNet algorithm tends to be better at detecting agricultural or natural forest classes. When it comes to situation where transitions happen during the classification period, as in the case of the transition from forest to pasture in 2004, random forest models tend to be more efficient than TempCNN or ResNet. 

## Attention-based models

Attention-based deep learning models are a class of models that use a mechanism inspired by human attention to focus on specific parts of input during processing. These models have been shown to be effective for various tasks such as machine translation, image captioning, and speech recognition.

The basic idea behind attention-based models is to allow the model to selectively focus on different input parts at different times. This is done by introducing a mechanism that assigns weights to each element of the input, indicating the relative importance of that element to the current processing step. The model uses them to compute a weighted sum of the input. The results capture the model's attention on specific parts of the input.

Attention-based models have become one of the most used deep learning architectures for problems that involve sequential data inputs, e.g., text recognition and automatic translation. The general idea is that not all inputs are alike in applications such as language translation. Consider the English sentence "Look at all the lonely people". A sound translation system needs to relate the words "look" and "people" as the key parts of this sentence to ensure such link is captured in the translation. A specific type of attention models, called transformers, enables the recognition of such complex relationships between input and output sequences [@Vaswani2017]. 

The basic structure of transformers is the same as other neural network algorithms. They have an encoder that transforms textual input values into numerical vectors and a decoder that processes these vectors to provide suitable answers. The difference is how the values are handled internally. The two main differences between transformer models and other algorithms are positional encoding and self-attention. Positional encoding assigns an index to each input value, ensuring that the relative locations of the inputs are maintained throughout the learning and processing phases. Self-attention compares every word in a sentence to every other word in the same sentence, including itself. In this way, it learns contextual information about the relation between the words. This conception has been validated in large language models such as BERT [@Devlin2019] and GPT-4 [@Bubeck2023].

The application of attention-based models for satellite image time series analysis is proposed by Garnot et al. [@Garnot2020a] and Russwurm and Körner [@Russwurm2020]. A self-attention network can learn to focus on specific time steps and image features most relevant for distinguishing between different classes. The algorithm tries to identify which combination of individual temporal observations is most relevant to identify each class. For example, crop identification will use observations that capture the onset of the growing season, the date of maximum growth, and the end of the growing season. In the case of deforestation, the algorithm tries to identify the dates when the forest is being cut. Attention-based models are a means to identify events that characterize each land class.

The first model proposed by Garnot et al. is a full transformer-based model [@Garnot2020a]. Considering that image time series classification is easier than natural language processing, Garnot et al. also propose a simplified version of the full transformer model [@Garnot2020]. This simpler model uses a reduced way to compute the attention matrix, reducing time for training and classification without loss of quality of the result. 

In `sits`, the full transformer-based model proposed by Garnot et al. [@Garnot2020a]   is implemented using `sits_tae()`. The default parameters are those proposed by the authors. The default optimizer is `optim_adamw`, as also used in the `sits_tempcnn()` function.

:::{.panel-tabset}
## R
```{r}
#| eval: false
# Train a machine learning model using TAE
tae_model <- sits_train(samples_matogrosso_mod13q1, 
                       sits_tae(
                          epochs               = 80,
                          batch_size           = 64,
                          optimizer            = torch::optim_adamw,
                          validation_split     = 0.2,
                          verbose              = FALSE))

# Show training evolution
plot(tae_model)
```
## Python
```{python}
#| eval: false
# Train a machine learning model using TAE
tae_model = sits_train(samples_matogrosso_mod13q1, 
                       sits_tae(
                          epochs               = 80,
                          batch_size           = 64,
                          optimizer            = "torch::optim_adamw",
                          validation_split     = 0.2,
                          verbose              = False))

# Show training evolution
plot(tae_model)
```
:::

```{r}
#| label: fig-ml-taemodel
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: |
#|   Training evolution of Temporal Self-Attention model.
knitr::include_graphics("./images/mltaemodel.png")
```

Then, we classify a 16-year time series using the TAE model. 

:::{.panel-tabset}
## R
```{r}
#| eval: false
# Classify using DL model and plot the result
point_mt_mod13q1 |>  
    sits_classify(tae_model) |>  
    plot(bands = c("NDVI", "EVI"))
```
## Python
```{python}
#| eval: false
# Classify using DL model
point_class = sits_classify(point_mt_mod13q1, tae_model)

# Plot
plot(
    point_class, bands = ("NDVI", "EVI")
)
```
:::

```{r}
#| label: fig-ml-taeplot
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: |
#|   Classification of time series using TAE.
knitr::include_graphics("./images/mltaeplot.png")
```


Garnot and co-authors also proposed the Lightweight Temporal Self-Attention Encoder (LTAE) [@Garnot2020], which the authors claim can achieve high classification accuracy with fewer parameters compared to other neural network models. It is a good choice for applications where computational resources are limited. The `sits_lighttae()` function implements this algorithm. The most important parameter to be set is the learning rate `lr`. Values ranging from 0.001 to 0.005 should produce good results. See also the section below on model tuning. 

:::{.panel-tabset}
## R
```{r}
#| eval: false
# Train a machine learning model using TAE
ltae_model <- sits_train(samples_matogrosso_mod13q1, 
                       sits_lighttae(
                          epochs               = 80,
                          batch_size           = 64,
                          optimizer            = torch::optim_adamw,
                          opt_hparams = list(lr = 0.001),
                          validation_split     = 0.2))

# Show training evolution
plot(ltae_model)
```
## Python
```{python}
#| eval: false
# Train a machine learning model using TAE
ltae_model = sits_train(samples_matogrosso_mod13q1, 
                       sits_lighttae(
                          epochs               = 80,
                          batch_size           = 64,
                          optimizer            = "torch::optim_adamw",
                          opt_hparams = dict(lr = 0.001),
                          validation_split      = 0.2))

# Show training evolution
plot(ltae_model)
```
:::

```{r}
#| label: fig-ml-ltaemodel
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: |
#|   Training evolution of Lightweight Temporal Self-Attention model.
knitr::include_graphics("./images/mlltaemodel.png")
```

Then, we classify a 16-year time series using the LightTAE model. 

:::{.panel-tabset}
## R
```{r}
#| eval: false
# Classify using DL model and plot the result
point_mt_mod13q1 |>  
    sits_classify(ltae_model) |>  
    plot(bands = c("NDVI", "EVI"))
```
## Python
```{python}
#| eval: false
# Classify using DL model
point_class = sits_classify(point_mt_mod13q1, ltae_model)

# Plot
plot(
    point_class, bands = ("NDVI", "EVI")
)
```
:::

```{r}
#| label: fig-ml-ltaeplot
#| echo: false
#| out.width: "90%"
#| fig.align: "center"
#| fig.cap: |
#|   Classification of time series using LightTAE.
knitr::include_graphics("./images/mlltaeplot.png")
```

The behaviour of both `sits_tae()` and `sits_lighttae()` is similar to that of `sits_tempcnn()`. It points out the possible need for more classes and training data to better represent the transition period between 2004 and 2010. One possibility is that the training data associated with the Pasture class is only consistent with the time series between the years 2005 to 2008. However, the transition from Forest to Pasture in 2004 and from Pasture to Agriculture in 2009-2010 is subject to uncertainty since the classifiers do not agree on the resulting classes. In general, deep learning temporal-aware models are more sensitive to class variability than Random Forest and extreme gradient boosters. 

## Summary

In this chapter, we present a basic description of the machine learning algorithms available in `sits`. The basic distinction is between time-sensitive algorithms such as LightTAE and TempCNN and those who do not consider temporal order of values, such as Random Forest. In practice, we suggest that users take Random Forest as their baseline and then use LightTAE or TempCNN to try to improve classification accuracy. 

## References{-}