Metadata-Version: 2.1
Name: trufl
Version: 0.0.4
Summary: Optimize adaptive sampling
Home-page: https://github.com/franckalbinet/trufl
Author: Floris Abrams, Franck Albinet
Author-email: franckalbinet@gmail.com
License: Apache Software License 2.0
Keywords: nbdev jupyter notebook python
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Developers
Classifier: Natural Language :: English
Classifier: Programming Language :: Python :: 3.7
Classifier: Programming Language :: Python :: 3.8
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: License :: OSI Approved :: Apache Software License
Requires-Python: >=3.7
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: fastcore
Requires-Dist: geopandas
Requires-Dist: rasterio
Requires-Dist: h3==4.0.0b5
Requires-Dist: contextily
Requires-Dist: pysal==24.1
Requires-Dist: gstools
Requires-Dist: ipywidgets
Provides-Extra: dev

# Trufl


<!-- WARNING: THIS FILE WAS AUTOGENERATED! DO NOT EDIT! -->

**Trufl** was initiated in the context of the [IAEA (International
Atomic Energy Agency)](https://www.iaea.org) Coordinated Research
Project titled [“Monitoring and Predicting Radionuclide Uptake and
Dynamics for Optimizing Remediation of Radioactive Contamination in
Agriculture”](https://www.iaea.org/newscenter/news/new-crp-monitoring-and-predicting-radionuclide-uptake-and-dynamics-for-optimizing-remediation-of-radioactive-contamination-in-agriculture-crp-d15019).

While **Trufl** was originally developed to address the remediation of
farmland affected by nuclear accidents, its approach and algorithms are
**applicable to a wide range of application domains**. This includes
managing **legacy contaminants or monitoring any phenomenon that
requires consideration of multiple decision criteria**, potentially
involving a large set of data.

This package leverages the work done by [Floris
Abrams](https://www.linkedin.com/in/floris-abrams-59080a15a) in the
context of his PhD at [KU Leuven](https://www.kuleuven.be) and [Franck
Albinet](https://www.linkedin.com/in/franckalbinet), International
Consultant in Geospatial Data Science and currently PhD researcher in AI
applied to nuclear remedation at KU Leuven.

## Install

`pip install trufl`

## Getting started

In highly sensitive and high-stakes situations, it is **essential that
decision making is informed, transparent, and accountable**, with
decisions being based on a thorough and objective analysis of the
available data and the needs and concerns of affected communities being
taken into account.

Given the time constraints and limited budgets that are often associated
with data surveys (in particular ones supposed to informed highligh
sensitive situation), it is **crucial to make informed decisions about
how to allocate resources**. This is even more important when
considering the many variables that can be taken into account, such as
prior knowledge of the area, health and economic impacts, land use,
whether remediation has already taken place, population density, and
more. Our approach leverages Linear Programming techniques such as
**Multiple-criteria decision-making** to optimize the data survey
workflow:

In this demo, we will walk you through a **typical workflow** using the
`Trufl` package. To help illustrate the process, we will use a “toy”
dataset that represents a typical spatial pattern of soil contaminants.

1.  We **assume that we have access to the ground truth**, which is a
    raster file that shows the spatial distribution of a soil
    contaminant;
2.  We will make decisions about how to sample and how much at
    **administrative units** level, which in this case are **simulated
    as a grid** (using the
    [`gridder`](https://franckalbinet.github.io/trufl/utils.html#gridder)
    utilities function);
3.  Based on prior knowledge of the phenomenon of interest, such as
    prior airborne surveys or other data, an
    [`Optimizer`](https://franckalbinet.github.io/trufl/optimizer.html#optimizer)
    will **rank each grid cell according to its priority for sampling**;
4.  We will then **perform random sampling on the designated grid
    cells** (using a
    [`Sampler`](https://franckalbinet.github.io/trufl/sampler.html#sampler)).
    To simulate the measurement process, we will use the ground truth to
    emulate measurements at each location (using a
    [`DataCollector`](https://franckalbinet.github.io/trufl/collector.html#datacollector));
5.  We will **evaluate the new state of each cell based on the
    measurements** and **pass it to a new round of optimization**. This
    process will be repeated iteratively to refine the sampling
    strategy.

### Imports

``` python
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import rasterio
import geopandas as gpd

from trufl.utils import gridder
from trufl.sampler import Sampler, rank_to_sample
from trufl.collector import DataCollector
from trufl.callbacks import (State, MaxCB, MinCB, StdCB, CountCB, MoranICB, PriorCB)
from trufl.optimizer import Optimizer

red, black = '#BF360C', '#263238'
```

### Our simulated ground truth

The assumed ground truth reveals a typical spatial pattern of
contaminant such as `Cs137` after a nuclear accident for instance.

``` python
fname_raster = './files/ground-truth-01-4326-simulated.tif'
with rasterio.open(fname_raster) as src:
    plt.axis('off')
    plt.imshow(src.read(1))
```

![](index_files/figure-commonmark/cell-3-output-1.png)

### Simulate administrative units

The sampling strategy will be determined on a per-grid-cell basis within
the administrative unit. We define below a 10 x 10 grid over the area of
interest:

``` python
gdf_grid = gridder(fname_raster, nrows=10, ncols=10)
gdf_grid.head()
```

<div>

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }
&#10;    .dataframe tbody tr th {
        vertical-align: top;
    }
&#10;    .dataframe thead th {
        text-align: right;
    }
</style>

|        | geometry                                          |
|--------|---------------------------------------------------|
| loc_id |                                                   |
| 0      | POLYGON ((-1.20830 43.26950, -1.20830 43.26042... |
| 1      | POLYGON ((-1.20830 43.27858, -1.20830 43.26950... |
| 2      | POLYGON ((-1.20830 43.28766, -1.20830 43.27858... |
| 3      | POLYGON ((-1.20830 43.29673, -1.20830 43.28766... |
| 4      | POLYGON ((-1.20830 43.30581, -1.20830 43.29673... |

</div>

</div>

> [!TIP]
>
> Note how each administrative unit is uniquely identified by its
> `loc_id`.

``` python
gdf_grid.boundary.plot(color=black, lw=0.5)
plt.axis('off');
```

![](index_files/figure-commonmark/cell-5-output-1.png)

### What prior knowledge do we have?

At the initial time $t_0$, data sampling has not yet begun, but we can
often **leverage existing prior knowledge of our phenomenon** of
interest to inform our sampling strategy/policy. In the context of
nuclear remediation, this prior knowledge can often be obtained through
mobile surveys, such as airborne or carborne surveys, which can provide
a **coarse estimation** of soil contamination levels.

In the example below, we **simulate prior information about the soil
property of interest by calculating the average value of the property
over each grid cell**.

At this stage, we have no measurements, so we simply create an empty
[Geopandas
GeoDataFrame](https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoDataFrame.html).

``` python
samples_t0 = gpd.GeoDataFrame(index=pd.Index([], name='loc_id'), 
                              geometry=None, data={'value': None})
```

> [!TIP]
>
> We need to set an index `loc_id` and have a `geometry` and `value`
> columns.

Now we get/“sense” the state of our grid cells based on the simulated
prior (Mean over each grid cell
[`PriorCB`](https://franckalbinet.github.io/trufl/callbacks.html#priorcb)):

``` python
state = State(samples_t0, gdf_grid, cbs=[PriorCB(fname_raster)])

# You have to call the instance
state_t0 = state(); state_t0.head()
```

<div>

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }
&#10;    .dataframe tbody tr th {
        vertical-align: top;
    }
&#10;    .dataframe thead th {
        text-align: right;
    }
</style>

|        | Prior    |
|--------|----------|
| loc_id |          |
| 0      | 0.102492 |
| 1      | 0.125727 |
| 2      | 0.161802 |
| 3      | 0.184432 |
| 4      | 0.201405 |

</div>

</div>

> [!TIP]
>
> We get the `Prior` for each individual `loc_id` (here only the first 5
> shown). The current
> [`State`](https://franckalbinet.github.io/trufl/callbacks.html#state)
> is only composed of a single
> [`PriorCB`](https://franckalbinet.github.io/trufl/callbacks.html#priorcb)
> variable but can include many more variables as we will see below.

### Optimize sampling based on prior at $t_0$

``` python
# Bugs?
#   - fails when ran twice

benefit_criteria = [True]
optimizer = Optimizer(state=state_t0)
df_rank = optimizer.get_rank(is_benefit_x=benefit_criteria, w_vector = [1],  
                             n_method=None, c_method = None, 
                             w_method=None, s_method="CP")

df_rank.head()
```

<div>

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }
&#10;    .dataframe tbody tr th {
        vertical-align: top;
    }
&#10;    .dataframe thead th {
        text-align: right;
    }
</style>

|        | rank |
|--------|------|
| loc_id |      |
| 92     | 1    |
| 93     | 2    |
| 91     | 3    |
| 94     | 4    |
| 82     | 5    |

</div>

</div>

``` python
gdf_grid.join(df_rank, how='left').plot(column='rank',
                                        cmap='viridis_r', 
                                        legend_kwds={'label': 'Rank'}, 
                                        legend=True)
plt.axis('off');
```

![](index_files/figure-commonmark/cell-9-output-1.png)

### Start sampling

> [!TIP]
>
> It’s worth noting that in the absence of any prior knowledge, a
> uniform sampling strategy over the area of interest may be used.
> However, this approach may not be the most efficient use of the
> available data collection and analysis budget.

Based on the **ranks (sampling priority)** calculated by the
[`Optimizer`](https://franckalbinet.github.io/trufl/optimizer.html#optimizer)
and given sampling **budget**, let’s calculate the number of samples to
be collected for each administrative unit (`loc_id`):

``` python
n = rank_to_sample(df_rank['rank'].sort_index().values, 
                   budget=300, min=1); n
```

    array([ 1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
            1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
            1,  1,  2,  3,  2,  1,  1,  1,  1,  1,  1,  3,  4,  7,  3,  1,  1,
            1,  1,  1,  2,  3,  4,  4,  2,  1,  1,  1,  2,  1,  2,  2,  2,  1,
            1,  1,  1,  1,  3,  2,  2,  2,  2,  2,  1,  1,  2,  5, 12,  8, 10,
            5,  4,  3,  1,  1,  6, 19, 58, 29, 14,  6,  3,  2,  1,  1])

We can now decide where to sample based on this sampling schema:

``` python
sampler = Sampler(gdf_grid)
sample_locs_t0 = sampler.sample(n, method='uniform')

print(sample_locs_t0.head())
ax = sample_locs_t0.plot(markersize=2, color=red)

gdf_grid.boundary.plot(color=black, lw=0.5, ax=ax)
plt.axis('off');
```

                             geometry
    loc_id                           
    0       POINT (-1.22083 43.26342)
    1       POINT (-1.22006 43.27588)
    2       POINT (-1.21642 43.28474)
    3       POINT (-1.21536 43.29169)
    4       POINT (-1.20880 43.30279)

![](index_files/figure-commonmark/cell-11-output-2.png)

### Emulating data collection and analysis

The data collector collects measurements at the random sampling
locations in the field. In our case, we emulate this process by
extracting measurements from the provided raster file.

“Measuring” variable of interest from a given raster:

``` python
dc_emulator = DataCollector(fname_raster)
samples_t0 = dc_emulator.collect(sample_locs_t0)

print(samples_t0.head())
ax = samples_t0.plot(column='value', s=2, legend=True)
gdf_grid.boundary.plot(color=black, lw=0.5, ax=ax);
plt.axis('off');
```

                             geometry     value
    loc_id                                     
    0       POINT (-1.22083 43.26342)  0.144753
    1       POINT (-1.22006 43.27588)  0.149566
    2       POINT (-1.21642 43.28474)  0.174893
    3       POINT (-1.21536 43.29169)  0.192808
    4       POINT (-1.20880 43.30279)  0.208858

![](index_files/figure-commonmark/cell-12-output-2.png)

### Getting current state

``` python
state = State(samples_t0, gdf_grid, cbs=[
    MaxCB(), MinCB(), StdCB(), CountCB(), MoranICB(k=5), PriorCB(fname_raster)
])

# You have to call the instance
state_t0 = state(); state_t0
```

<div>

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }
&#10;    .dataframe tbody tr th {
        vertical-align: top;
    }
&#10;    .dataframe thead th {
        text-align: right;
    }
</style>

|        | Max      | Min      | Standard Deviation | Count | Moran.I  | Prior    |
|--------|----------|----------|--------------------|-------|----------|----------|
| loc_id |          |          |                    |       |          |          |
| 0      | 0.145231 | 0.074017 | 0.025631           | 8     | 0.915078 | 0.102492 |
| 1      | 0.152964 | 0.111763 | 0.015033           | 7     | 0.679907 | 0.125727 |
| 2      | 0.160229 | 0.160229 | 0.000000           | 1     | NaN      | 0.161802 |
| 3      | 0.188177 | 0.164264 | 0.008248           | 6     | 0.940746 | 0.184432 |
| 4      | 0.261005 | 0.241690 | 0.009657           | 2     | NaN      | 0.201405 |
| ...    | ...      | ...      | ...                | ...   | ...      | ...      |
| 95     | 0.877876 | 0.800729 | 0.027397           | 5     | NaN      | 0.803670 |
| 96     | 0.804376 | 0.795942 | 0.004217           | 2     | NaN      | 0.763408 |
| 97     | 0.799064 | 0.672111 | 0.047581           | 8     | 0.948441 | 0.727797 |
| 98     | 0.708847 | 0.661171 | 0.013101           | 9     | 0.731343 | 0.646002 |
| 99     | 0.704167 | 0.673213 | 0.015477           | 2     | NaN      | 0.655185 |

<p>100 rows × 6 columns</p>
</div>

</div>

## Build the ranking of polygons based on several criteria

### Criteria

- MaxCB()
- MinCB()
- StdCB()
- CountCB()
- MoranICB(k=5) – Gives 2 values (value , p-value)
- PriorCB

### Criteria type

- Benefit (high values –\> high score –\> rank high –\> prioritized
  sampling needed)

- Cost (high values –\> low score –\> low high –\> Less sampling needed)

- MaxCB() – Benefit

- MinCB() – ???

- StdCB() – Benefit

- CountCB() – Cost (Low count – higher priority because more samples
  need)

- MoranICB(k=5) – Cost (high value – highly correlated – less need for
  sampling ?? )

- PriorCB – Benefit

### MCDM techniques

- CP – low values – good alternative
- TOPSIS – High Value – good alternative

! Everything is converted to rank to account for these differences !

``` python
benefit_criteria = [True, True, True]
state = State(samples_t0, gdf_grid, cbs=[MaxCB(), MinCB(), StdCB()])
```

``` python
optimizer = Optimizer(state=state())
df = optimizer.rank(is_benefit_x=benefit_criteria, w_vector = [0.3, 0.3, 0.4],  
                    n_method=None, c_method = None, w_method=None, s_method="CP")

df.head()
```

<div>

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }
&#10;    .dataframe tbody tr th {
        vertical-align: top;
    }
&#10;    .dataframe thead th {
        text-align: right;
    }
</style>

|        | rank |
|--------|------|
| loc_id |      |
| 93     | 1    |
| 92     | 2    |
| 84     | 3    |
| 91     | 4    |
| 83     | 5    |

</div>

</div>

``` python
# https://kapernikov.com/ipywidgets-with-matplotlib/
```
