# Strategy: LPJmL benchmarking

## History

* (-2017) there used to be an [R benchmarking script](https://gitlab.pik-potsdam.de/lpjml/LPJmL_internal/-/wikis/Benchmarking) that could compare the output of two different LPJmL versions regarding:
1. global sums (with literature data)
2. global sums time series
3. difference maps
4. cell time series
* with eddy flux data (NEP, ET)
* with runoff/dischage
6. regional validation with FAO data

* (2018-today) [LPJmLBenchmark function](https://gitlab.pik-potsdam.de/lpjml/LPJmL_internal/-/wikis/LPJmLBenchmark) was established to serve same functionality but be more flexible
* flexibility: in principle every LPJmL output could be used
* new versions of LPJmL (nitrogen, 5arcmin, ...) could be benchmarked
* functionality as before (with exception of FAO data)

## Status quo
`LPJmLBenchmark()` is *outdated*, *unflexible*, *smelly*, *slow* and not *maintainable* ...
* *outdated* it is based on the `LandMark` package, an internal extension of the `raster` and the meanwhile deprecated `sp` package
* *unflexible* it uses a procedural approach that handles the processing of data in long functions without any abstraction to improve understandability. Also function can handle only the rendering of static pdf reports
* [*code smells*](https://en.wikipedia.org/wiki/Code_smell) the code is not DRY, functions are way to long and complicated, too much parameters and object handling, abuse of global environment and more
* *slow* `raster` package lacks performance when it comes to time series data, 30 minutes for a benchmark is too long
* *maintainable* missing strategy in the beginning of `LPJmLBenchmark()` development lead to an acretion of new features rather than a structure that allowed new features to be added, adjusted and maintained with little effort

## Requirements

#### Revised PDF report
* better overview: e.g. table of contents, more accessible and intiuitive content (usage of latest rmarkdown tools)
* reduced complexity or better description of metrics
* additional metrics
* country maps: e.g. FAO data
* lat plots: e.g. NPP data
* ...
* assessment column of sums data (e.g. <5% of soil carbon deviation OK , > 5% not OK)
* revised literature data with additional bibliography at end of report
* ...

####

#### OOP approach based on lpjmlkit
* Why object oriented programming (OOP)?
* human thinking is bound to things/objects and its *natural* composition, e.g. *a room object usually has furniture objects inside of it.*
* modularity facilitates troubleshooting
* reuse of code through inheritance
* flexibility through polymorphism
* in general effective problem solving
* ...
* inherit lpmlkit's `LPJmLData` class to write a new and more generic `LandData` class as well as `LandaDataSet` to execute methods for multiple outputs at once
* to read and process other data types, e.g. region/country-wise FAO data
* inherit `LPJmLMetaData` class to write `LandMetaData` class with less and more generic attributes
* add arithmetics (+, -, *, /) to `LandData` class
* `aggregate()` & `disaggregate()` methods to switch between (to be introduced) new space formats: global, country, region, ... , cell
* `Metric` class to include required data as `LandDataSet`s for each metric with generic (polymorph) method `create_metric()` to e.g. create a sum table or a sum time series
* `MetricSet` class to (analogously to LandDataSet) execute methods for multiple Metrics at once
* `LPJmLBenchmark` as the main benchmark object to store all metrics and the required meta data for a benchmark as well as methods to `print()`, `report()` and more

#### Fast Metrics based on preprocessed data
Each metric is based on a required data format. The underlying conversion and processing should be easily doable also outside of the benchmark reporting to also calculate them easily independently:

Regarding the mentioned metrics this could look like the following:
* **global sums** *following lpjmlkit syntax* e.g.:

    # test dataset (the output to be benchmarked)
    # area contains the area for every cell of x   
   test_sums <- aggregate(x = test_dataset * area,
                           cell = sum,
                           year = mean)

    # reference dataset (the benchmark)
    ref_sums <- aggregate(x = ref_dataset * area,
                          cell = sum,
                          year = mean)

* **global sums timeseries** *following lpjmlkit syntax* e.g.:

    # plot ts of test dataset
    plot(x = test_dataset * area,
         aggregate = list(cell = sum))

    # add ts of reference dataset (the benchmark) to plot
    plot(x = ref_dataset * area,
         aggregate = list(cell = sum),
         add = TRUE)

* **raster (difference) maps** *following lpjmlkit syntax* e.g.:

    # absolute maps -------------------------------------- #
    # plot map of test dataset
    plot(x = test_dataset * area,
         aggregate = list(year = mean))

    # add ts of reference dataset (the benchmark) to plot
    plot(x = ref_dataset * area,
         aggregate = list(year = mean),
         add = TRUE)

    # difference maps ------------------------------------ #
    # plot difference map of test and reference dataset
    plot(x = test_dataset * area - ref_dataset * area,
         aggregate = list(year = mean))

* ...

## Implementation

#### Generic `LandData` model

- use the [R6 class system of package R6](https://r6.r-lib.org) as already done for [`lpjmlkit::LPJmLData`](https://github.com/PIK-LPJmL/lpjmlkit/blob/master/R/LPJmLData.R)
- `LandData`: child class of the `LPJmLData` class that is to be extended thus inheriting its functionality that can then be adjusted if required. An example is the [`LPJmLGridData` class](https://github.com/PIK-LPJmL/lpjmlkit/blob/master/R/LPJmLGridData.R)
- Arithmetics: add arithmetics to `LandData` objects e.g. (`harvest <- yield * cftfrac * cellarea`)
- `aggregate()` function to calculate metrics e.g. for
- global sums `aggregate(harvest, cell = sum, time = mean)`
- global sums time series `aggregate(harvest, cell = sum)`
- regional validation, e.g. FAO data: `aggregate(harvest, cell = sum, use_regions = TRUE)`
- ...
- `disaggregate()` function (regions to cells)

- `LandMetaData`: child class of the [`LPJmLMetaData` class](https://github.com/PIK-LPJmL/lpjmlkit/blob/master/R/LPJmLMetaData.R) (e.g. [`LPJmLGridData` class](https://github.com/PIK-LPJmL/lpjmlkit/blob/master/R/LPJmLGridData.R))
- you may delete attributes such as `sim_name`, `order` or `map`
- you may add new attributes such as ...
- `LandDataSet` is a set of multiple `LandData` objects

- `units` package implementation for data and meta data

%%{ init: { 'theme': 'base', 'themeVariables': {'primaryColor': '#ffffff', 'primaryBorderColor': '#000000'} } }%%

flowchart TD
    classDef new fill:#FFFF00
    subgraph " "

    lpjmlkit::LPJmLMetaData-."inherit".-LandMetaData:::new
    LandMetaData-."compose \n (1:1)".-LandData:::new
    lpjmlkit::LPJmLData-."inherit".-LandData:::new
    grid_file[grid file]-.-> lpjmlkit::LPJmLGridData

    lpjmlkit::LPJmLGridData-."compose \n (1:1) \n optional".-LandData
    cow_file[COW file]-.-> lpjmlkit::LPJmLRegionsData:::new
    lpjmlkit::LPJmLRegionsData:::new-."compose \n (1:1) \n optional".-LandData:::new

    LandData---Arithmetics("Arithmetics (+, -, *, /)"):::new
    LandData---aggregate("aggregate()"):::new
    LandData---disaggregate("disaggregate()"):::new
    LandData---create_table("create_table()"):::new

    LandData-."compose \n (n:1)".-LandDataSet:::new

    NEW{{TODO}}:::new 
    end

#### Generic `read_data()` function

%%{ init: { 'theme': 'base', 'themeVariables': {'primaryColor': '#ffffff', 'primaryBorderColor': '#000000'} } }%%

flowchart
    classDef new fill:#FFFF00

    subgraph " "


    %% data type handling
    read_data:::new-.->type{type}

    type{type}-. "output (clm/meta/cdf)" .->read_io("lpjmlkit::read_io()")
    type{type}-. "tabular (csv)" .->read_table("read_table()"):::new
    cow_file[lpjmlkit::LPJmLRegionsData]:::new-.-> read_table
    grid_file[lpjmlkit::LPJmLGridData]-.-> read_table


    read_io-.->LandData$new("LandData$new()"):::new-.->LandData:::new
    read_table:::new -.-> LandData:::new

    NEW{{TODO}}:::new 
    EXISTING{{EXISTING}}


    end

#### `benchmark()` initialization

The `benchmark()` function could look like the following:

# Benchmarking function definition, returns an object of class LPJmLBenchmark, that includes required LandDataSets for every metric
benchmark <- function(x,
                      ref,
                      subset = NULL, # list(year = as.character(1990:2019)),
                      # define either output set flags ("default" benchmarking output set, "all" intersection of x and ref outputs, "landuse" with additional/only landuse outputs) or define vector of output ids to be benchmarked
                      output = "default" #, "all", "landuse", ..., c("vegc", "soilc", "litc", "runoff", ...)
                      metric = c("sum_table",
                                 "sum_ts",
                                 "diff_map",
                                 "abs_map",
                                 "site_ts",
                                 "region_map",
                                 "lat_dist") # ...
                      dataset = c("fao", 
                                  "flux",
                                  "runoff"), # ...
                      ...) {
    ... 
}

The `benchmark()` initialization will initialize the benchmarking by creating a `LPJmLBenchmark` object, that includes the processed data and the methods to e.g. generate a benchmark report based on the data
1. `benchmark()` will use the information available in the configurations (json files) to generate meta information as attributes (`author`, `data`, `description`, etc.) of the returned `LPJmLBenchmark` object
2. `benchmark()` will intialize a `MetricSet`, a set of multiple (or a single) `Metric`s that are set via the metric argument, e.g. sum table or sum time series
* to save memory benchmark will iterate over the outputs, read each and process it directly be store it in the corresponding `Metric`'s `LandDataSet` (x, ref, ...) - this step may be parallelized.
* some objects may be hold in memory that are used to calculate intermediate products/sums, like `cftfrac` and `grid`/cellarea for `harvest`
3. If all metrics in `MetricSet` are filled with processed `LPJmLDataSet`s the LPJmLBenchmark data object is returned

Idealized structure of `benchmark(x, ref, subset, output, metric, dataset, ...)` as a flow chart:

%%{ init: { 'theme': 'base', 'themeVariables': {'primaryColor': '#ffffff', 'primaryBorderColor': '#000000'} } }%%

flowchart TD
    classDef new fill:#FFFF00

    subgraph " "
    %% output reading


    test_outputs["to test \n outputs"]-->read_data2("read_data()"):::new
    read_data2("read_data()"):::new-->LandData2["LandData"]:::new
    LandData2:::new-."compose \n (n:1)".-LandDataSet:::new

    ref_outputs["reference \n outputs"]-->read_data3("read_data()"):::new
    read_data3("read_data()"):::new-->LandData3["LandData"]:::new
    LandData3:::new-."compose \n (n:1)".-LandDataSet2["LandDataSet"]:::new


    %% LandData:::new<-."to be harmonized \n for each metric".->LandData2["LandData"]:::new

    fao_data[FAO data]-->read_data("read_data()"):::new
    discharge_data[Discharge data]-->read_data("read_data()"):::new
    flux_data[Flux data]-->read_data("read_data()"):::new-->LandData:::new
    dots["..."]-->read_data("read_data()"):::new


    Metric:::new-."inherit".-SumTable

    LandDataSet:::new--"x"-->SumTable["SumTable(x,ref, ...)"]
    LandDataSet2:::new--"ref"-->SumTable["SumTable(x,ref, ...)"]
    LandData:::new--"..."-->SumTable["SumTable(x,ref, ...)"]

    SumTable:::new--"compose \n (n:1)"-->MetricSet:::new
    MetricSet:::new--"compose \n (1:1)"-->LPJmLBenchmark:::new


    SumTable:::new-."OR".-SumTS:::new
    SumTable:::new-."OR".-DiffMap:::new
    SumTable:::new-."OR".-AbsMap:::new
    SumTable:::new-."OR".-SiteTS:::new
    SumTable:::new-."OR".-RegionMap:::new
    SumTable:::new-."OR".-LatDist:::new


    EXISTING{{EXISTING}}
    NEW{{TODO}}:::new 
    end

Notes:
* `benchmark()` knows about ideal order of data processing to save additional processing (information could be stored in a yaml file), e.g.:
1. `SumTS` object first (space aggregation)
2. `SumTable` based on `SumTS$x`, `SumTS$ref` and `SumTS$...` (+ time aggregation)
* `Metric` (and child classes) initialization (`$new()`) uses harmonization function to ensure LPJmLData formats are equal/compatible between `x`, `ref` and `...`
* `Metric` (and child classes) initialization (`$new()`) also ensures that required intermediate products/sums of variables (e.g. "harvest" or "net ecosystem production") are calculated

* ...

#### Benchmark datasets storage

* Examples are:
* FAO data
* flux data
* dischage/runoff data
* lat distribution data (@sybills)
* face data (?)
* remote sensing data (?)
* ...

* the data would either have to be readable (within `read_data()`) via internal `read_table()` or lpjmlkit::`read_io()` function

* it will be tested to store these datasets as RData/RDS or CSV (or bin/clm/NetCDF) files in the the package
* using lazy loading this should not affect the package loading and memory: https://rstudio.github.io/r-manuals/r-exts/Creating-R-packages.html#data-in-packages
* otherwise easy, non error-prone alternatives have to be found (cluster, database, ...)!

* Also additional functionality should be provided to generate an overview over the already available data.

* for time series that would have to steadily be updated (e.g. FAO data) dataset deprecation warnings could be printed to request maintainers to update the data sets from time to time

#### Create benchmark report (and more)

An `LPJmLBenchmark` object is the main benchmarking interface that should provide multiple methods to enable different types of output generation
* `print()` method
* `report()` generation
* and more (e.g. [shiny dashboard](https://shiny.rstudio.com)):

%%{ init: { 'theme': 'base', 'themeVariables': {'primaryColor': '#ffffff', 'primaryBorderColor': '#000000', 'lineColor': '#000000'} } }%%

flowchart TD
    classDef new fill:#FFFF00
    classDef future fill:#FF6800

    subgraph " "
    benchmark("benchmark(x, ref)"):::new --> LPJmLBenchmark["LPJmLBenchmark"]

    LPJmLBenchmark["LPJmLBenchmark"]:::new --> print("print()"):::new

    LPJmLBenchmark["LPJmLBenchmark"]:::new --> report("report(format = 'pdf')"):::new

    LPJmLBenchmark["LPJmLBenchmark"]:::new --> dashboard("dashboard(...)"):::future

    LPJmLBenchmark["LPJmLBenchmark"]:::new --> gitlab_runner("gitlab_runner(...)"):::future

    LPJmLBenchmark["LPJmLBenchmark"]:::new --> other("..."):::future

    NEW{{TODO}}:::new 
    FUTURE{{FUTURE \n WORK}}:::future 
    end
     

*Not all of functions are to be implemented now but might be at a later stage (FUTURE WORK)*

A typical benchmark workflow could look like this:

# Usually only x and ref argument for configs are required
my_benchmark <- benchmark(x = "./config_x.json", ref = "./config_ref.json")

# Check auto generated information
my_benchmark$author
# jannesbr
my_benchmark$date
my_benchmark$description

# And adjust if required
my_benchmark$description <- "Compare LPJmL 4 with LPJmL 5"

# Save LPJmLBenchmark object to later reproduce report, use in gitlab runner or for a shiny dashboard
saveRDS(my_benchmark, "my_benchmark.rds")

# Generate pdf report
report(my_benchmark, "./my_benchmark_report.pdf")

> TEST COMMENT [color=#e2449e]
> [name=Jannes Breier]
> [time=Tue, Apr 4, 2023 1:18 PM]