Table of Contents
- [Overview](#overview)
- [Preparation](#preparation)
- [Set cell density](#set-cell-density)
- [Interaction function](#interaction-function)
- [Simulation](#simulation)
- [Creating a random mosaic](#creating-a-random-mosaic)
- [Check usable features](#check-usable-features)
- [Optimization Target](#optimization-target)
- [Optimization Routine](#optimization-routine)
- [Annealing schedule](#annealing-schedule)
- [Visualization on optimization](#visualization-on-optimization)
- [Extention](#extention)
- [Features and Entropy](#features-and-entropy)
- [Annealing Schedule](#annealing-schedule-1)
# Overview
Before generating a new artificial mosaic, we show information on `pattern` through the previous tutorial, as
```python
print(pattern)
```
The outputs are
```texture
Spatial pattern of Mouse Horizontal Cell,
- Density: Unknown,
- Natural mosaics: 1 case(s),
- Simulated mosaics: total 30 case(s)
0 case(s) in tag 'default',
10 case(s) in tag 'O-PIPP',
20 case(s) in tag 'PIPP',
- Features: 4
Label | Has target probabilities
NN | True
VD | True
NNRI | False
VDRI | False .
```
# Preparation
The `Pattern` require more information on natural mosaics for simulation and optimization, including the density of cells and the interaction function between two cells.
## Set cell density
`Pattern` has an `estimate_density` to analyze cell density in natural mosaics, as
```python
density = pattern.estimate_density()
# 87/90000.
```
The cell density tells the number of cells in a mosaic, which is essential in mosaic simulation. We use the `set_density` method to add the density into the `pattern`, as
```python
pattern.set_density(density=density)
pattern.set_density(87/90000.) # direct input the density
print(pattern)
```
Therefore, the "Density" in the outputs has its value.
```
Spatial pattern of Mouse Horizontal Cell,
- Density: 0.0009667 (cells/unit^2),
- Natural mosaics: 1 case(s),
- Simulated mosaics: total 30 case(s)
0 case(s) in tag 'default',
10 case(s) in tag 'O-PIPP',
20 case(s) in tag 'PIPP',
- Features: 4
Label | Has target probabilities
NN | True
VD | True
NNRI | False
VDRI | False .
```
## Interaction function
The simulation of mosaic requires an **interaction function**, denoted as `h(u)` to estimate the probability of a distance between any two points in a spatial pattern. With theoretical works from the field of [spatial point pattern analysis](), the formation of the interaction function is flexible. It only has one constraint as `0 <= h(u) <=1` for `u >= 0`.
Here we recommend a well-known formation of `h(u)` for retinal mosaics, as
$$h(u)= \begin{cases} 0,\quad u\leq δ\\ 1-exp(-((u-δ)/φ)^{α}), \quad u>δ \end{cases},$$
where δ, φ, α are parameters estimated by the Poisson point process. Besides, we recommend a [`R` script](estimate_inter_ps.md) to get parameters in $h(u)$.
With estimated parameters, we can get the callable function through the `get_interaction_func` method, as
```python
parameters = [7.5, 32.1206741, 2.64876305] # [δ, φ, α]
h_func = pattern.get_interaction_func(parameters)
```
# Simulation
The routine of mosaic simulation is the **insert-update-optimize** framework. It creates a random pattern and updates a cell's position following the probability yield by the interaction function. After an update iteration with several cells, it calculates the performance of the simulated mosaic and uses an adaptive simulated annealing algorithm to ensure the road of the mosaic is towards to given features and their probability distributions.
In this section, we use NN and VD features as the optimization target and show how to generate a new mosaic in `OPIPP`.
## Creating a random mosaic
The 1st step of the simulation is to create a random mosaic through the `new_mosaic` method, as
```python
mosaic = pattern.new_mosaic(scope=scope)
mosaic.draw_points()
```
Random points at initialization.
You need input a `Scope` to tell the `pattern` how large is the plane.
## Check usable features
Next, you need to decide on features in optimization. A usable feature should have the optimization target through `Distribution.set_target` or `Pattern.set_feature_target`. The `pattern` has four features and only the "NN" and the "VD" feature have target probabilities. You can use the `get_usable_features` to check all usable features in the `pattern`.
```python
usable_features = pattern.get_usable_features()
# It is ["NN", "VD"]
```
## Optimization Target
Once we have decided on features for optimization, we can evaluate the `Entropy` of a mosaic by calculating its KL divergence to target features, as
$$Entropy(Simulated)=\sum_{feature}KL(Feature(Simulated), Target_{feature})),$$
where $feature$ and $Feature$ denote a spatial feature and the related method that calculates the probabilistic distribution of features in a given mosaic. The $Entropy$ is the sum of KL divergence between distributions from the simulated mosaic and the $target$.
In `OPIPP`, you can use the `evaluate` method to calculate the entropy of a mosaic or a series of mosaics. For example, we calculate the entropy of the random mosaic as
```python
# input the mosaic and target features
print(pattern.evaluate([mosaic], features=usable_features))
# return is 5.233049009275092
# If the input is a list of mosaics, it will calculate distributions of these mosaics and yield a single value.
```
## Optimization Routine
The `simulate` method in `pattern` play the role of update and optimization.
```python
from OPIPP import AdaptiveSchedule
mosaic, losses = pattern.simulate(mosaic=mosaic, interaction_func=h_func, features=None, schedule=AdaptiveSchedule(), max_step=None, update_ratio=None, save_prefix="examples/simulated/HC/Sample", save_step=500, verbose=True)
# the entropy of the final mosaic
print(pattern.evaluate([mosaic], features=usable_features))
# 0.0288679366606423
mosaic.draw_points()
```
Simulated mosaic after optimization.
Arguments in `simulate` are
- `mosaic`: The input mosaic.
- `interaction_func`: The interaction function, default=`None`. If it is `None`, the method will use `h(u)=1.0` to accept all updates in cells.
- `features`: The target features for optimization, default=`None`. If it is `None`, the method will use all usable features.
- `schedule`: Annealing schedule for simulated annealing algorithm in optimization, default=`AdaptiveSchedule()`. [Here](3.simulation.md#annealing-schedule) is a more complete description. If it is `None`, the method will accept all updates after an iteration.
- `max_step`: The maximum number of iteration steps, default=`None`. If it is `None` and `schedule` is `None` as well, the method will use `20` as the default value. Otherwise, the method will loop until the schedule is terminated or reaches the maximum step of iterations.
- `update_ratio`: The ratio of cells in an iteration step, default=`None`. If it is `None` and `schedule` is `None` as well, the method will use `1.0` as the default value. If it is `None` but `schedule` is not `None`, the method will use `0.01` in an iteration.
- `save_prefix`: Save the mosaics into local files if given. The output file is `save-prefix_index-of-iteraction.points`.
- `save_step`: The step of saving into local files, default=`1`.
- `verbose`: Whether print the change of entropy during optimization, default=`True`.
`loss` is the trace of entropies alongside the optimization. You can plot it as
```python
plt.plot(losses)
plt.show()
```
The trace of Entropy during optimization
Besides, we recommend using `MPI` to simulate multiple mosaics in parallel. [Here](parallel_processing.md) is an example.
## Annealing schedule
The `OPIPP` use the [simulated annealing algorithm](https://en.wikipedia.org/wiki/Simulated_annealing) to optimize the performance of the simulated mosaic. The algorithm uses a `temperature` to estimate the probability of accepting worse cases and requires an annealing schedule to control the process of iteration. We recommend using the `AdaptiveSchedule` in practice. The creation of an `AdaptiveSchedule` object is
```python
from OPIPP import AdaptiveSchedule
schedule = AdaptiveSchedule(alpha=0.95, init_t=0.5, min_t=1e-4)
```
Arguments are
- `alpha`: The descent parameter in the adaptive schedule, defalut=`0.95`.
- `init_t`: The value of temperature at initialization, default=`0.5`.
- `min_t`: The value of temperature for termination, default=`0.0001`.
Please check the [adaptive simulated annealing algorithm](https://optimization-online.org/wp-content/uploads/2001/03/291.pdf) for more information.
## Visualization on optimization
After optimization, you load files and view simulation results in the `pattern`. Here, we simulate 10 mosaics and 20 mosaics with different simulation parameters. Output files obtaining points in mosaics are stored in [examples/simulated/HC/](examples/simulated/HC). We use [glob](https://docs.python.org/3/library/glob.html) and `Pattern.load_from_files` to load multiple mosaics, as
```python
from glob import glob
# load simulated mosaics by the O-PIPP method
points_files = glob("examples/simulated/HC/W1_*.points")
pattern.load_from_files(points_files, scope=scope, is_natural=False, simulated_tag="O-PIPP")
# load simulated mosaics by the PIPP method
points_files = glob("examples/simulated/HC/PIPP_*.points")
pattern.load_from_files(points_files, scope=scope, is_natural=False, simulated_tag="PIPP")
```
Then, you can use visualization methods in `pattern` to show values of features and entropy. For example, we compare features in two simulated groups as
```python
# boxplot with feature values
pattern.draw_feature_boxes(feature_label="NN", draw_natural=True, simulated_tags=["O-PIPP", "PIPP"])
pattern.draw_feature_boxes(feature_label="VD", draw_natural=True, simulated_tags=["O-PIPP", "PIPP"])
```
Values of nearest Neighbor distances among mosaics
Values of Voronoi domain areas among mosaics
Furthermore, you can let `draw_loss=True` in `draw_value_bars` to draw values of entropy, as
```python
# bars indicate the mean entropy of features in two groups of mosaics
pattern.draw_value_bars(value_method=np.mean, feature_colors={"VD": "skyblue", "NN": "blue"}, draw_loss=True, draw_natural=False, simulated_tags=["O-PIPP", "PIPP"])
```
Compare entropy of two groupds of simulated mosaics
# Extention
The **insert-update-optimize** framework in `OPIPP` is welcome for spatial features and annealing schedules proposed by users. Here we summarize how to implement customized features and annealing schedules in `OPIPP`.
## Features and Entropy
The core of a feature is the `Distribution` class, deciding how to extract features, target values for optimization, and how to calculate entropy during optimization.
A feature that each cell in the mosaic yields a value and the statistics of a population is significant can follow our [previous examples](2.analysis.md#define-a-distribution). However, there are several features not fit the diagram. For instance, the regularity index of a mosaic is a single value. If you want to use these features in mosaic optimization, you should define a new class that inherits the `Distribution` and several methods are overridden, including
- `set_target`: Let it know the target value (or values) for optimization. Arguments are flexible. No return value.
- `has_target`: Let it judge itself if it has the target for optimization. No argument for this method. The return is `True` or `False`.
- `extract_mosaic`: Let it know how to calculate features from a single mosaic. It has an argument, the mosaic(`OPIPP.Mosaic`) for processing. The return is a numpy.darray.
- `extract_mosaics`: Let it know how to calculate features from a list of mosaics. It has an argument, the list of mosaics(`OPIPP.Mosaic`) for processing. The return is a numpy.darray.
- `entropy`: Let it know how to calculate the entropy with given values. It has an argument, the list of values return by `extract_mosaic` or `extract_mosaics`. The return is a single value.
## Annealing Schedule
`Schedule` class in `OPIPP` is a general definition of the annealing schedule. You can import it and create a new schedule. For example, we define a schedule for log annealing, as
```python
from OPIPP.cooling import Schedule
import numpy as np
class LogSchedule(Schedule): # inherits the original class
def __init__(self, base: float=2, min_t: float=1e-4, max_update: int=None):
"""
Override the initialization method.
The following attributes must be decided,
- The `min_t` is the threshold. The simulation will be terminated if the temperature is below it.
- The `max_update` is the max number of iterations. The simulation will be terminated if it reaches the given value. If it is `None`, there is no limitation on the number of iterations.
"""
# parameters for the log schedule
self.base = base
# use the initialization method in `Schedule`
Schedule.__init__(self, min_t=min_t, max_update=max_update)
def init(self):
"""
This method is called before simulation.
You need to set an initial temperature (self.t) and finish the other preparation.
"""
self.t = self.c/(np.log(2)/np.log(self.base))
Schedule.init(self)
def update(self, loss):
"""
The schedule needs to process a new loss (entropy) and update the temperature (self.t) inside.
Besides, the Schedule has an `i_update` attribute to indicate the index of the new loss since the latest `init`.
"""
self.t = self.c/(np.log(1+self.i_update)/np.log(self.base))
```