From 7676df4a7084275be2e403bfb47d855b6f4c9c20 Mon Sep 17 00:00:00 2001
From: Steve Ataucuri
Date: Wed, 16 Sep 2020 16:37:58 +0000
Subject: [PATCH 1/3] Moved config files
---
config_cnn.json => configs/config_cnn.json | 0
config_cnn_parallel.json => configs/config_cnn_parallel.json | 0
.../config_cnn_parallel_internal.json | 0
config_default.json => configs/config_default.json | 0
config_lstm.json => configs/config_lstm.json | 0
config_lstm_gaussian.json => configs/config_lstm_gaussian.json | 0
config_lstm_parallel.json => configs/config_lstm_parallel.json | 0
.../config_lstm_parallel_internal.json | 0
config_mlp.json => configs/config_mlp.json | 0
config_mlp_gaussian.json => configs/config_mlp_gaussian.json | 0
config_rnn.json => configs/config_rnn.json | 0
11 files changed, 0 insertions(+), 0 deletions(-)
rename config_cnn.json => configs/config_cnn.json (100%)
rename config_cnn_parallel.json => configs/config_cnn_parallel.json (100%)
rename config_cnn_parallel_internal.json => configs/config_cnn_parallel_internal.json (100%)
rename config_default.json => configs/config_default.json (100%)
rename config_lstm.json => configs/config_lstm.json (100%)
rename config_lstm_gaussian.json => configs/config_lstm_gaussian.json (100%)
rename config_lstm_parallel.json => configs/config_lstm_parallel.json (100%)
rename config_lstm_parallel_internal.json => configs/config_lstm_parallel_internal.json (100%)
rename config_mlp.json => configs/config_mlp.json (100%)
rename config_mlp_gaussian.json => configs/config_mlp_gaussian.json (100%)
rename config_rnn.json => configs/config_rnn.json (100%)
diff --git a/config_cnn.json b/configs/config_cnn.json
similarity index 100%
rename from config_cnn.json
rename to configs/config_cnn.json
diff --git a/config_cnn_parallel.json b/configs/config_cnn_parallel.json
similarity index 100%
rename from config_cnn_parallel.json
rename to configs/config_cnn_parallel.json
diff --git a/config_cnn_parallel_internal.json b/configs/config_cnn_parallel_internal.json
similarity index 100%
rename from config_cnn_parallel_internal.json
rename to configs/config_cnn_parallel_internal.json
diff --git a/config_default.json b/configs/config_default.json
similarity index 100%
rename from config_default.json
rename to configs/config_default.json
diff --git a/config_lstm.json b/configs/config_lstm.json
similarity index 100%
rename from config_lstm.json
rename to configs/config_lstm.json
diff --git a/config_lstm_gaussian.json b/configs/config_lstm_gaussian.json
similarity index 100%
rename from config_lstm_gaussian.json
rename to configs/config_lstm_gaussian.json
diff --git a/config_lstm_parallel.json b/configs/config_lstm_parallel.json
similarity index 100%
rename from config_lstm_parallel.json
rename to configs/config_lstm_parallel.json
diff --git a/config_lstm_parallel_internal.json b/configs/config_lstm_parallel_internal.json
similarity index 100%
rename from config_lstm_parallel_internal.json
rename to configs/config_lstm_parallel_internal.json
diff --git a/config_mlp.json b/configs/config_mlp.json
similarity index 100%
rename from config_mlp.json
rename to configs/config_mlp.json
diff --git a/config_mlp_gaussian.json b/configs/config_mlp_gaussian.json
similarity index 100%
rename from config_mlp_gaussian.json
rename to configs/config_mlp_gaussian.json
diff --git a/config_rnn.json b/configs/config_rnn.json
similarity index 100%
rename from config_rnn.json
rename to configs/config_rnn.json
From 644473bcd37a1604ac262f25c2418617ad1e7207 Mon Sep 17 00:00:00 2001
From: Steve Ataucuri
Date: Sun, 4 Oct 2020 01:23:47 +0000
Subject: [PATCH 2/3] Last implementations for gaussian model
---
README.md | 96 +-
configs/config_cnn.json | 10 +-
configs/config_cnn_parallel.json | 15 +-
configs/config_lstm.json | 8 +-
configs/config_lstm_gaussian.json | 34 +-
configs/config_lstm_parallel.json | 29 +-
configs/config_mlp.json | 16 +-
core/data/data_loader.py | 19 +-
core/models/base.py | 207 +++-
core/models/gaussian_loss.py | 14 +-
core/models/lstm.py | 14 +-
core/utils/metrics.py | 8 +-
core/utils/transformation.py | 11 +
imgs/plot_eta.png | Bin 0 -> 18806 bytes
imgs/plot_phi.png | Bin 0 -> 20450 bytes
imgs/plot_polar.png | Bin 0 -> 123142 bytes
main_inference.py | 106 +-
main_train.ipynb | 292 ++---
main_train.py | 19 +-
notebooks/main_test-cylind.ipynb | 624 -----------
notebooks/plot_distributions.ipynb | 1681 ++++++++++++++++++++++++++++
notebooks/plot_prediction.ipynb | 2 +-
22 files changed, 2301 insertions(+), 904 deletions(-)
create mode 100644 imgs/plot_eta.png
create mode 100644 imgs/plot_phi.png
create mode 100644 imgs/plot_polar.png
delete mode 100644 notebooks/main_test-cylind.ipynb
create mode 100644 notebooks/plot_distributions.ipynb
diff --git a/README.md b/README.md
index 2b96da7..693aee7 100644
--- a/README.md
+++ b/README.md
@@ -38,45 +38,60 @@ Considering the mentioned regions and the symmetry of the detector, each dataset
-A short datasets are in `dataset` directory.
+Short datasets are in `dataset` directory. We split a region in three parts (this is too big to be here).
# Running
## Training
There are some predefined config files to train diferents models (MLP, CNN, LSTM, CNN-parallel and others). If you need to change the parameters then change the config_*.json file. We used internal barrel as dataset, this dataset is previously transformed and linked in json file:
```sh
-$ python main_train.py --config config_lstm_parallel_internal.json
+$ python main_train.py --config configs/config_lstm_parallel_internal.json
```
There are other configurations for example a CNN model:
```sh
-$ python main_train.py --config config_cnn_parallel_internal.json
+$ python main_train.py --config configs/config_cnn_parallel_internal.json
```
-If you want to see the training process when ajust any parameters of .json file. Run the notebook:
+If you want to use a gaussian model
```sh
-$ main_train.ipynb
+$ python main_train.py --config configs/config_lstm_gaussian.json
```
-For many trainings and testings with scripts. You can run it, with the default configuration:
+
+If you want to see the training process when ajust any parameters of .json file. Run the notebook:
```sh
-$ ./run_trains.sh
-$ ./run_tests.show
+$ main_train.ipynb
```
## Inference
-You can inference data test:
+You can inference the previous compiled model with the same config file:
```sh
-$ python main_inference.py --config config_lstm_parallel_internal.json
+$ python main_inference.py --config configs/config_lstm_parallel_internal.json --load True
```
-This will produce a `results/encrypt_name/results-test.txt` file.
+This will produce an report at `results/encrypt_name/results-test.txt` file. With the number of reconstructed tracks and quantity of hits correctly predicted by layer.
+If you want to inference with a short dataset:
+```sh
+$ python main_inference.py --config configs/config_lstm_parallel_internal.json --load True --samples 500
+```
+or
+```sh
+$ python main_inference.py --config configs/config_lstm_gaussian.json --load True
+```
## Auxiliary Scripts
+You can automatizate trainings and testings with scripts. there are some default configurations:
+
+```sh
+$ ./run_trains.sh
+$ ./run_tests.sh
+```
# Performance
## Accuracy of Algorithm
+
We are using regressions metrics for accuracy of models. We show 2 groups of metrics.
- The principal metrics is a scoring. Scoring counts how many correct hits were found per layer and comparates with original truth hits. Finally we count the quantity of tracks reconstructed.
@@ -89,18 +104,24 @@ Output test file:
```
[Output] Results
---Parameters---
- Model Name : lstm
- Dataset : phi025-025_eta025-025_train1_lasthit_20200219.csv
- Tracks : 528
- Model saved : /compiled/model-lstm-DCtuvkiXn32hugVsTaokcp-coord-xyz-normalise-true-epochs-21-batch-6.h5
- Test date : 10/06/2020 12:09:34
- Coordenates : xyz
+ Model Name : gaussian-lstm
+ Dataset : eta_n0.5-0.5_phi_n0.5-0.5_20200518171238_tracks.csv
+ Tracks : 3000
+ Model saved : compiled/model-gaussian-lstm-MEnGHa5DCmJGE9C9BpAg4i-coord-cylin-normalise-true-epochs-30-batch-30.h5
+ Test date : 02/10/2020 20:23:55
+ Coordenates : cylin
+ Coordenate 3D : True
Model Scaled : True
Model Optimizer : adam
- Prediction Opt : nearest
- Total correct hits per layer [256. 251. 213. 194. 157. 126.] of 528 tracks tolerance=0.0:
- Total porcentage correct hits : ['48.48%', '47.54%', '40.34%', '36.74%', '29.73%', '23.86%']
- Reconstructed tracks: 74 of 528 tracks
+ Prediction Opt : gaussian
+ Distance metric : ['polar', 'euclidean', 'mahalanobis', 'cosine']
+ Remove hit : False
+ Correct hits per layer (polar) ['191.0(6.37%)', '135.0(4.5%)', '182.0(6.07%)', '167.0(5.57%)', '121.0(4.03%)', '127.0(4.23%)'] of 3000 tracks:
+ Correct hits per layer (euclidean) ['61.0(2.03%)', '54.0(1.8%)', '64.0(2.13%)', '62.0(2.07%)', '49.0(1.63%)', '36.0(1.2%)'] of 3000 tracks:
+ Correct hits per layer (mahalanobis) ['191.0(6.37%)', '193.0(6.43%)', '178.0(5.93%)', '193.0(6.43%)', '165.0(5.5%)', '133.0(4.43%)'] of 3000 tracks:
+ Correct hits per layer (cosine) ['2323.0(77.43%)', '1980.0(66.0%)', '1792.0(59.73%)', '1702.0(56.73%)', '1487.0(49.57%)', '1351.0(45.03%)'] of 3000 tracks:
+ Reconstructed tracks: 887 of 3000 tracks (29.566666666666666)
+
```
Above output shows scoring per layer for example 48% with 256 hits were matched at the first layer, results are 74 tracks reconstructed of 528 tracks(it is a short dataset just). We also write other info like what kind of coordinate, if we use the nearest optimization, epochs, batchs, optimazer used, model name etc.
@@ -138,23 +159,30 @@ layer 7
The last output shows one geral metric for all hits and four (R^2, MSE, RMSE, MAE) metrics per layer.
-## Vizualization
-If you want to see the results with plots, go to the plot_prediction.ipynb file at `notebooks` directory.
+## Plots
+If you want to see the results with plots, go to the plot_prediction.ipynb file at `notebooks` directory. This example plots cartesian coordinates.
-This plot is 10 tracks reconstructed.
-
-
-
+
+ |
+
The next plot shows all hits.
-
-
-
+
+ |
+
-The next plot is the prediction of all hits.
-
-
-
+### Plot Distributions
+There is other notebook for plot the gaussian distributions plot_distributions.ipynb
+
+
+ |
+
+
+The pull of distributions for $etha$ and $phi$.
+
+ |
+ |
+
diff --git a/configs/config_cnn.json b/configs/config_cnn.json
index 29c6f18..c993f10 100644
--- a/configs/config_cnn.json
+++ b/configs/config_cnn.json
@@ -5,7 +5,7 @@
"log_dir": "logs"
},
"data": {
- "filename": "/data/track-ml/eramia/dataset/eta-05_05.csv",
+ "filename": "dataset/eta_n0.5-0.5_phi_n0.5-0.5_short.csv",
"train_split": 0.80,
"normalise": true,
"cylindrical": false,
@@ -21,12 +21,12 @@
"use_gpu": true,
"show_metrics": false,
"shuffle": true,
- "earlystopping": true
+ "earlystopping": true
},
"testing": {
- "type_optimization" : "normal",
+ "type_optimization" : "nearest",
"normalise": false,
- "tolerance": 0.02
+ "tolerance": 0.01
},
"model": {
"name": "cnn",
@@ -53,7 +53,7 @@
},
{
"type": "dense",
- "neurons": 18,
+ "neurons": 12,
"activation": "relu"
},
{
diff --git a/configs/config_cnn_parallel.json b/configs/config_cnn_parallel.json
index b3ba31e..d71542e 100644
--- a/configs/config_cnn_parallel.json
+++ b/configs/config_cnn_parallel.json
@@ -5,17 +5,19 @@
"log_dir": "logs"
},
"data": {
- "filename": "./dataset/phi025-025_eta025-025_train1_lasthit_20200219.csv",
+ "filename": "dataset/eta_n0.5-0.5_phi_n0.5-0.5_short.csv",
"train_split": 0.80,
"normalise": true,
- "cylindrical": false,
+ "type_norm": "maxmin",
+ "cylindrical": true,
+ "points_3d": true,
"num_hits": 10,
"features": 3
},
"training": {
- "epochs": 15,
+ "epochs": 18,
"batch_size": 6,
- "validation": 0.2,
+ "validation": 0.3,
"save_model": true,
"load_model": false,
"use_gpu": true,
@@ -25,8 +27,9 @@
},
"testing": {
"type_optimization" : "nearest",
+ "metric": "cosine",
"normalise": false,
- "tolerance": 0.02
+ "tolerance": 0.01
},
"model": {
"name": "cnn-parallel",
@@ -53,7 +56,7 @@
},
{
"type": "dense",
- "neurons": 64,
+ "neurons": 32,
"activation": "relu"
},
{
diff --git a/configs/config_lstm.json b/configs/config_lstm.json
index 7b4edf3..9c27d2c 100644
--- a/configs/config_lstm.json
+++ b/configs/config_lstm.json
@@ -5,7 +5,7 @@
"log_dir": "logs"
},
"data": {
- "filename": "/data/track-ml/eramia/dataset/phi025-025_eta025-025_train1_lasthit_20200219.csv",
+ "filename": "dataset/eta_n0.5-0.5_phi_n0.5-0.5_short.csv",
"train_split": 0.80,
"normalise": true,
"cylindrical": false,
@@ -13,9 +13,9 @@
"features": 3
},
"training": {
- "epochs": 21,
+ "epochs": 20,
"batch_size": 6,
- "validation": 0.2,
+ "validation": 0.33,
"save_model": true,
"load_model": false,
"use_gpu": true,
@@ -46,7 +46,7 @@
},
{
"type": "dropout",
- "rate": 0.25
+ "rate": 0.2
},
{
"type": "lstm",
diff --git a/configs/config_lstm_gaussian.json b/configs/config_lstm_gaussian.json
index 9a30eae..0659fea 100644
--- a/configs/config_lstm_gaussian.json
+++ b/configs/config_lstm_gaussian.json
@@ -5,19 +5,19 @@
"log_dir": "logs"
},
"data": {
- "filename": "./dataset/eta_n0.5-0.5_phi_n0.5-0.5_internal_short1.csv",
+ "filename": "/data/track-ml/bracis/datasets/eta_n0.5-0.5_phi_n0.5-0.5/eta_n0.5-0.5_phi_n0.5-0.5_20200518171238_tracks.csv",
"train_split": 0.80,
"normalise": true,
- "type_norm": "zscore",
+ "type_norm": "maxmin",
"cylindrical": true,
- "points_3d": false,
+ "points_3d": true,
"num_hits": 10,
- "features": 2
+ "features": 3
},
"training": {
"epochs": 30,
- "batch_size": 18,
- "validation": 0.3,
+ "batch_size": 30,
+ "validation": 0.35,
"save_model": true,
"load_model": false,
"use_gpu": true,
@@ -26,27 +26,29 @@
"earlystopping": true
},
"testing": {
- "type_optimization" : "nearest",
- "metric": "cosine",
+ "type_optimization" : "gaussian",
+ "metric": ["polar", "euclidean", "mahalanobis", "cosine"],
"normalise": false,
- "tolerance": 0.02
- },
+ "tolerance": 0.02,
+ "remove_hit": false,
+ "show_metrics": false
+ },
"model": {
"name": "gaussian-lstm",
"overwrite": true,
- "isparallel": false,
+ "isparallel": false,
"loss": "mse",
"metrics": ["acc"],
"optimizer": "adam",
- "learningrate": 0.01,
+ "learningrate": 0.0005,
"layers": [
{
"type": "input",
- "neurons": 500,
- "input_features": 2,
+ "neurons": 700,
+ "input_features": 3,
"input_timesteps": 4,
- "activation": "tanh"
- }
+ "activation": "tanh"
+ }
]
}
}
diff --git a/configs/config_lstm_parallel.json b/configs/config_lstm_parallel.json
index 3cb6991..8b4e603 100644
--- a/configs/config_lstm_parallel.json
+++ b/configs/config_lstm_parallel.json
@@ -5,27 +5,27 @@
"log_dir": "logs"
},
"data": {
- "filename": "./dataset/phi025-025_eta025-025_train1_lasthit_20200219.csv",
+ "filename": "./dataset/eta_2.0-3.0_phi_n0.5-0.5_external_short.csv",
"train_split": 0.80,
"normalise": true,
"cylindrical": false,
"num_hits": 10,
"features": 3
-
},
"training": {
- "epochs": 18,
- "batch_size": 24,
+ "epochs": 10,
+ "batch_size": 30,
"validation": 0.2,
"save_model": true,
"load_model": false,
"use_gpu": true,
"show_metrics": false,
- "shuffle": true,
+ "shuffle": false,
"earlystopping": true
},
"testing": {
"type_optimization" : "nearest",
+ "metric": "euclidean",
"normalise": false,
"tolerance": 0.01
},
@@ -39,21 +39,30 @@
"layers": [
{
"type": "lstm",
- "neurons": 500,
+ "neurons": 800,
"input_timesteps": 4,
- "input_features": 1,
+ "input_features": 1,
"return_seq": false,
"stateful": false
},
{
"type": "dropout",
- "rate": 0.25
+ "rate": 0.5
},
{
"type": "dense",
- "neurons": 64,
+ "neurons": 100,
"activation": "relu"
- },
+ },
+ {
+ "type": "dropout",
+ "rate": 0.2
+ },
+ {
+ "type": "dense",
+ "neurons": 18,
+ "activation": "relu"
+ },
{
"type": "dense",
"neurons": 3,
diff --git a/configs/config_mlp.json b/configs/config_mlp.json
index b7420c1..5b50646 100644
--- a/configs/config_mlp.json
+++ b/configs/config_mlp.json
@@ -21,10 +21,11 @@
"use_gpu": true,
"show_metrics": false,
"shuffle": true,
- "earlystopping": false
+ "earlystopping": true
},
"testing": {
"type_optimization" : "nearest",
+ "metric": "cosine",
"normalise": false,
"tolerance": 0.02
},
@@ -34,23 +35,24 @@
"isparallel": false,
"loss": "mse",
"metrics": ["acc"],
- "optimizer": "rmsprop",
+ "optimizer": "adam",
+ "learningrate": 0.001,
"layers": [
{
"type": "mlp",
- "neurons": 12,
+ "neurons": 90,
"input_features": 3,
"input_timesteps": 4,
- "activation": "relu"
+ "activation": "tanh"
},
{
"type": "dense",
- "neurons": 6,
- "activation": "relu"
+ "neurons": 60,
+ "activation": "tanh"
},
{
"type": "dense",
- "neurons": 3,
+ "neurons": 2,
"activation": "linear"
}
]
diff --git a/core/data/data_loader.py b/core/data/data_loader.py
index 07f2a4f..78ebde8 100644
--- a/core/data/data_loader.py
+++ b/core/data/data_loader.py
@@ -24,7 +24,7 @@ class KindNormalization(Enum):
Nothing = 4
class Dataset():
- def __init__(self, input_path, train_size, cylindrical, hits, kind_normalization, points_3d=True):
+ def __init__(self, input_path, train_size, cylindrical, hits, kind_normalization, points_3d=True, samples=None):
#np.set_printoptions(suppress=True)
@@ -34,8 +34,8 @@ def __init__(self, input_path, train_size, cylindrical, hits, kind_normalization
self.kind = kind_normalization
if self.kind == KindNormalization.Scaling:
- self.x_scaler = MinMaxScaler(feature_range=(-1, 1))
- self.y_scaler = MinMaxScaler(feature_range=(-1, 1))
+ self.x_scaler = MinMaxScaler(feature_range=(-np.pi, np.pi))
+ self.y_scaler = MinMaxScaler(feature_range=(-np.pi, np.pi))
elif self.kind == KindNormalization.Zscore:
self.x_scaler = StandardScaler() # mean and standart desviation
@@ -53,7 +53,12 @@ def __init__(self, input_path, train_size, cylindrical, hits, kind_normalization
self.interval = 11
self.decimals = 4
- self.data = dataframe.iloc[:, self.start_hits:]
+ # select an specific samples of dataset
+ if samples is not None:
+ self.data = dataframe.iloc[:samples, self.start_hits:]
+ else:
+ self.data = dataframe.iloc[:, self.start_hits:]
+
#self.self = 0
if cylindrical:
@@ -261,7 +266,7 @@ def get_testing_data(self, n_hit_in, n_hit_out, n_features, normalise=False, xsc
sequences = self.data_test.values
rows = sequences.shape[0]
-
+
for i in range(0, rows):
end_ix = n_hit_in*n_features
seq_x, seq_y = sequences[i, 0:end_ix], sequences[i, end_ix:]
@@ -282,7 +287,7 @@ def get_testing_data(self, n_hit_in, n_hit_out, n_features, normalise=False, xsc
y_data = pd.DataFrame(Y)
elif normalise and (xscaler is not None or yscaler is not None):
- print('not is none')
+
# we load a previous scaler mean and std
self.x_scaler = xscaler
self.y_scaler = yscaler
@@ -295,7 +300,7 @@ def get_testing_data(self, n_hit_in, n_hit_out, n_features, normalise=False, xsc
y_data = pd.DataFrame(Y)
elif not normalise and (xscaler is not None or yscaler is not None):
- print('normalise false and is none')
+
# we load a previous scaler mean and std
self.x_scaler = xscaler
self.y_scaler = yscaler
diff --git a/core/models/base.py b/core/models/base.py
index f477d7b..77230ea 100644
--- a/core/models/base.py
+++ b/core/models/base.py
@@ -20,6 +20,7 @@
from core.utils.utils import *
from core.models.gaussian_loss import gaussian_loss, gaussian_nll
+from scipy.stats import norm
class BagOfHits(Enum):
All=1,
@@ -58,7 +59,7 @@ def __init__(self, configs):
'model-%s-%s-coord-%s-normalise-%s-epochs-%s-batch-%s.h5' % (
self.name, self.encryp_ds_name, coord,
str(self.normalise).lower(), self.epochs, self.batch_size))
- print(self.save_fnameh5)
+ #print(self.save_fnameh5)
self.save_fname = os.path.join(configs['paths']['save_dir'], 'architecture-%s.png' % self.name)
@@ -635,6 +636,210 @@ def predict_full_sequences_nearest_parallel(self, x_test, y_test, data, bag_of_h
return pred_sequences, count_correct_nearest, count_correct
+ def predict_prob(self, x, bs):
+ """Make predictions given model and 2d data
+ """
+
+ ypred = self.model.predict(x, verbose=0, batch_size=bs)
+ n_outs = int(ypred.shape[1] / 2)
+ mean = ypred[:, 0:n_outs]
+ sigma = np.exp(ypred[:, n_outs:])
+
+ return mean, sigma
+
+ def predict_prob_full_sequences_nearest(self, x_test, y_test, data, bag_of_hits, y_test_xyz=None,
+ t_steps=4, t_features=1, n_features=3, num_hits=6,
+ normalise=False, cylindrical=False, verbose=False, tol=0.01,
+ remove_hit=False, points_3d=True, metrics = 'euclidean'):
+
+ '''
+ This func predict the mean and variance shift the window by 1 new prediction each time,
+ re-run predictions on new window
+
+ cylindrical : it says the training model was fed with cylindrical data if it is true
+ x_test and y_test are original data coordinates and not normalised
+ normalise : say if the model was compiled with normalised data
+ '''
+
+ timer = Timer()
+ timer.start()
+
+ print('[Model] Predicting Sequences with Nearest Started')
+
+ total = len(y_test)
+
+ pred_tracks = []
+ pred_sequences_orig = []
+
+ # if cylindrical is true then we change to cartesian coordinates
+ '''
+ if cylindrical:
+ y_test = y_test_aux
+ else:
+ y_test = y_test
+ '''
+ # bag_of_hits
+ #bag_of_hits_all = np.array(convert_matrix_to_vec(y_test, n_features))
+
+ layers, layers_xyz = [] , []
+ begin = 0
+ for i in range(num_hits):
+ end = begin+n_features
+ layer_orig = np.array(y_test.iloc[:,begin:end]).reshape(total, n_features)
+ layer_xyz = np.array(y_test_xyz.iloc[:,begin:end]).reshape(total, n_features)
+ begin = end
+ layers.append(layer_orig)
+ layers_xyz.append(layer_xyz)
+
+ corrects_by_layer = np.zeros((len(metrics), num_hits))
+ begin_idx, end_idx = 0, 0
+ num_boh = 6
+
+ for j in range(total):
+ # in original coordinates, if was trained with polar then it must be polar coordinate
+ curr_frame = x_test[j]
+ # bag_of_hit by track
+ curr_track = np.array(y_test.iloc[j,0:]).reshape(num_hits, n_features)
+
+ if verbose:
+ print('curr_track %s , %s:' % (j , curr_track))
+
+ predicted = []
+ predicted_orig = []
+ #begin = 0
+ for i in range(num_hits):
+
+ # this current hit is in original coordinates
+ curr_hit = curr_track[i]
+
+ if verbose:
+ # primeira esta em originais
+ print('input:\n', curr_frame)
+
+ # if the model was normalized, we need to transform
+ if normalise:
+ curr_frame = data.x_scaler.transform(np.reshape(curr_frame,(1, t_steps*n_features)))
+ curr_frame_orig = data.inverse_transform_x(pd.DataFrame(curr_frame).values)
+ curr_frame_orig = np.reshape(curr_frame_orig, (t_steps, n_features))
+ curr_frame = np.reshape(curr_frame, (t_steps, n_features))
+ else:
+ curr_frame = curr_frame
+ curr_frame_orig = curr_frame
+
+ if verbose:
+ print('input:\n', curr_frame)
+ #print('input orig:\n', curr_frame_orig)
+
+ #print('input inv :\n', curr_frame_inv.reshape(4,3))
+ #print('newaxis ', curr_frame[np.newaxis,:,:])
+ # input must be scaled curr_frame
+ # the input and output must be in original coordinates that was compiled
+ mean, sigma = self.predict_prob(curr_frame[np.newaxis,:,:], bs=1)
+
+ if points_3d:
+ rho_mean, rho_sigma = mean[:,0], sigma[:,0]
+ eta_mean, eta_sigma = mean[:,1], sigma[:,1]
+ phi_mean, phi_sigma = mean[:,2], sigma[:,2]
+
+ rho_pred = norm.median(rho_mean, rho_sigma)
+ eta_pred = norm.median(eta_mean, eta_sigma)
+ phi_pred = norm.median(phi_mean, phi_sigma)
+
+ y_pred = np.reshape([rho_pred, eta_pred, phi_pred], (1, n_features))
+
+ else:
+ eta_mean, eta_sigma = mean[:,0], sigma[:,0]
+ phi_mean, phi_sigma = mean[:,1], sigma[:,1]
+
+ eta_pred = norm.median(eta_mean, eta_sigma)
+ phi_pred = norm.median(phi_mean, phi_sigma)
+
+ y_pred = np.reshape([eta_pred, phi_pred], (1, n_features))
+
+ if normalise:
+ y_pred_orig = data.inverse_transform_y(y_pred)
+ else:
+ y_pred_orig = y_pred
+
+
+ # convert cylindrical to xyz for calculate the distance in the euclidean space
+ '''
+ if cylindrical:
+ rho, eta, phi = y_pred_orig[0][0], y_pred_orig[0][1], y_pred_orig[0][2]
+ #print(rho,eta,phi)
+ x, y, z = convert_rhoetaphi_to_xyz(rho, eta, phi)
+
+ #print(x,y,z)
+ y_pred_orig[0][0] = x
+ y_pred_orig[0][1] = y
+ y_pred_orig[0][2] = z
+
+ #print(y_pred_orig)
+ '''
+
+ curr_layer = layers[i]
+
+ for m, metric in enumerate(metrics):
+ # calculate the nearest hit with a polar distance function
+ if metric == 'polar':
+ if points_3d:
+ dist = distance.cdist(curr_layer, y_pred_orig, distance_cylindrical_3D)
+ else:
+ dist = distance.cdist(curr_layer, y_pred_orig, distance_cylindrical_2D)
+ else: # other distances
+ dist = distance.cdist(curr_layer, y_pred_orig, metric)
+
+ # get the hit with lowest distance in orignal coordinates
+ idx = np.argmin(dist)
+ most_nearest = curr_layer[idx]
+
+ # if curr_hit is in cartesian coord, near_pred must be in cartesian coord too
+ # very small numbers are differents or equals
+ # we count how many times we get the correct hit by layer
+ if np.isclose(curr_hit, most_nearest, atol=tol).all():
+ corrects_by_layer[m][i]+=1
+
+ # removing the hit from layer, mayority of times is not good
+ if remove_hit:
+ layers[i] = np.delete(layers[i], idx, 0)
+
+ if verbose:
+ print('pred:', y_pred)
+ print('inv pred:', y_pred_orig)
+ print('current:', curr_hit)
+ print('nearest:', most_nearest)
+
+ # i think it is not necessary
+ '''
+ if cylindrical:
+ y, z = near_pred[0], near_pred[1]
+ #print(x,y,z)
+ _, eta, phi = convert_xyz_to_rhoetaphi(1, y, z)
+ #print(rho,eta,phi)
+ #near_pred[0] = rho
+ near_pred[0] = eta
+ near_pred[1] = phi
+ '''
+ # we change to the original input at the same postion
+ #if cylindrical:
+ # near_pred = curr_layer_[idx]
+
+ #near adiciona em valores originaeis
+ predicted.append(most_nearest)
+
+ curr_frame = curr_frame_orig[1:]
+ curr_frame = np.insert(curr_frame, [3], predicted[-1], axis=0)
+
+
+ pred_tracks.append(predicted)
+
+ if verbose:
+ if j == 10: break
+
+ print('[Model] Prediction Finished.')
+ timer.stop()
+
+ return pred_tracks, corrects_by_layer
def nearest_hit(self, hit, hits,
silent = True,
diff --git a/core/models/gaussian_loss.py b/core/models/gaussian_loss.py
index ddd6ad4..9abb1fb 100644
--- a/core/models/gaussian_loss.py
+++ b/core/models/gaussian_loss.py
@@ -67,17 +67,17 @@ def gaussian_nll(ytrue, ypreds):
"""
- #n_dims = int(int(ypreds.shape[1])/2)
+ n_dims = int(int(ypreds.shape[1])/2)
- #mu = ypreds[:, 0:n_dims]
- #logsigma = ypreds[:, n_dims:]
- dims = ypreds.shape[1]
- mu = ypreds[:,0:1]
- logsigma = ypreds[:, 1:]
+ mu = ypreds[:, 0:n_dims]
+ logsigma = ypreds[:, n_dims:]
+ #dims = ypreds.shape[1]
+ #mu = ypreds[:,0:1]
+ #logsigma = ypreds[:, 1:]
mse = -0.5*K.sum(K.square((ytrue-mu)/K.exp(logsigma)),axis=1)
sigma_trace = -K.sum(logsigma, axis=1)
- log2pi = -0.5*dims*np.log(2*np.pi)
+ log2pi = -0.5*n_dims*np.log(2*np.pi)
log_likelihood = mse+sigma_trace+log2pi
diff --git a/core/models/lstm.py b/core/models/lstm.py
index 84db586..3dbd67f 100644
--- a/core/models/lstm.py
+++ b/core/models/lstm.py
@@ -170,20 +170,22 @@ def build_model(self):
inputs = Input(shape=(input_timesteps, input_features))
outputs = LSTM(neurons, return_sequences=False, activation='tanh')(inputs)
#outputs = LSTM(neurons, return_sequences=False, activation='tanh')(outputs)
- outputs = Dense(2, activation='linear')(outputs)
+ outputs = Dense(neurons/2, activation='linear')(outputs)
+ outputs = Dense(input_features*2, activation='linear')(outputs)
- distributions = Lambda(gaussian_layer_2d)(outputs)
+ #distributions = Lambda(gaussian_layer_2d)(outputs)
####################
- self.model = Model(inputs=inputs, outputs=distributions)
+ self.model = Model(inputs=inputs, outputs=outputs)
if configs['model']['optimizer'] == 'adam':
opt = Adam(lr=configs['model']['learningrate'])
elif configs['model']['optimizer'] == 'rmsprop':
opt = RMSprop(lr=configs['model']['learningrate'])
- print(self.model.summary())
- self.model.compile(loss=gaussian_nll, optimizer=opt, metrics=['accuracy'])
+
+ self.model.compile(loss=gaussian_nll, optimizer=opt, metrics=['accuracy'])
+ print(self.model.summary())
print('[Model] Model Compiled with structure:', self.model.inputs)
self.save_architecture(self.save_fname)
- timer.stop()
\ No newline at end of file
+ timer.stop()
diff --git a/core/utils/metrics.py b/core/utils/metrics.py
index 6e4732e..c4c8d49 100644
--- a/core/utils/metrics.py
+++ b/core/utils/metrics.py
@@ -41,15 +41,15 @@ def evaluate_training(history, save_to, keyword):
# acc=history['val_acc'][-1], loss=history['val_loss'][-1]))
report_string = 'accuracy: {acc}, loss: {loss}, val_acc: {val_acc}, val_loss: {val_loss}'.format(
- acc=history['val_acc'][-1],
+ acc=history['accuracy'][-1],
loss=history['val_loss'][-1],
- val_acc=history['val_acc'][-1],
+ val_acc=history['val_accuracy'][-1],
val_loss=history['val_loss'][-1])
print(report_string)
# summarize history for accuracy
- plt.plot(history['acc'])
- plt.plot(history['val_acc'])
+ plt.plot(history['accuracy'])
+ plt.plot(history['val_accuracy'])
plt.title('model accuracy')
plt.ylabel('accuracy')
plt.xlabel('epoch')
diff --git a/core/utils/transformation.py b/core/utils/transformation.py
index 0f225b8..d01d49a 100644
--- a/core/utils/transformation.py
+++ b/core/utils/transformation.py
@@ -142,7 +142,18 @@ def parts_from_tracks(tracks):
hits = tracks[:, 7:]
return indexes, vertices, momenta, hits
+def distance_cylindrical_2D(a, b):
+ r0, theta0 = a[0] , a[1]
+ r1, theta1 = b[0] , b[1]
+ return np.sqrt(r0**2 + r1**2 - 2*r0*r1*np.cos(theta0 - theta1))
+
+def distance_cylindrical_3D(a, b):
+ r0, eta0, phi0 = a[0] , a[1], a[2]
+ r1, eta1, phi1 = b[0] , b[1], b[2]
+
+ return np.sqrt(r0**2 + r1**2 - 2*r0*r1*np.cos(phi0 - phi1) + (eta0-eta1)**2)
+
# if __name__ == "__main__":
# print("Test 1: rotate hit to horizontal")
# hit = np.array([np.random.random(), np.random.random(), np.random.random()])
diff --git a/imgs/plot_eta.png b/imgs/plot_eta.png
new file mode 100644
index 0000000000000000000000000000000000000000..f74c7647bbeadc1abf16e08c4157c734b7de9af2
GIT binary patch
literal 18806
zcmb8X2RxSj`!{~sdu3&Zq6i70tgNVrLb69u$h_>mca%yQnMHQVy3EKZ3EA1JF4-$&
zJ;!;u@9*<_e$ROQ|Gi$_-Pbuj=W!mNV;{%+jMUP+LPgF>jzXcR)YX)4qEL7n@OP7x
z7>=N#C`I7Egf5EeI;8N&hxBd){ANL`=%a5rSfSm`o$sNn?HwN86LPV1zIV^w<-P-Y
znV?<{E)qa4QgXg$j(*@^e?sTM!+R(vSMw8M=T5lXL!UT*?)E|Mt9UYi8ZA*5l*XT9;ob*W38VT{>lN({&w8
zUgzN+jo;s%G}VPc|DASI$4=QTIIX-$M38zK-)r-C==Qe9$y29@4)&&)>rtqhy|wo3
zgROe5vaMH_+=cBfusumgc$|_FbFjZVJsHlWdhgynvovX!pJqL8Gzu^3+4^-4*89u1
z#IVq7mg1p8Lqi><*f(E!&pf_T;ynAMBu=Au4Bn3C?*8A++@d1tJwG4Q}
zuC89WSU!x}R^@tr_D)br3j0YuK1E&Kjo*C~hkpm^@t5b(7oC2ky7Nd&vjmgSw7in_
z#vf)-`S_Uk3Vz&v7K&NK1Z*S$WQcl(b5pCbTsoK~G0JnnSHsbP4
zSg`%Mkdme*_05|%9bv`BIdu7NTa}Y+$h&XO56}qe6*ROx=G8Vfj*W_nYABbm?Vsqr
z!lmNrDgFNa`<31CL%EWNgC*vJwce87(&XdP_Gb7irypGkBF~$dxf{bJU)OC)c;VrX
zkQ?&*4B<3fswyf(CjR@b=-<6$RaI3l#H^{AU@w0Ewka^JuQM?m_doR2(a~wJcqcCP
z6XQDX*vjEf_~VlAU$lqJE34rbFJ45%vPtsD%F3S0d;IE^i9=dbIF03>g&~ihFRdN_
zEZn*NN`|Fk)QxS#-rX!sa?FE(h=}O4uyD&*t@r9$yGX^~dF7Mz@OB~B)8shbzI_`g
zFH*CjSi0N!LX<#q#F5}E6UEQ3FM~&`-ATH*(}vax_!T-s!JQ1k*ImI#13Q*Hh=o
z&dtq9`5*WM|0>iI%*)Gj#CCC$a(5M*%A=RZYKJUkR~iV~+S(L+e0;Ko!$BOl_OVw=
z3<|dw9sDVBdLj!AD~S8XT51md&c8T+kGv(0izdG>iim!?D{_Rq3ZZ7Fw)z+AI^5|h+ei8ZJC2nrKzP`S|R^YJs-_7~(r%xHO
z-@M6u|Na;#%jL414<8sfF1gvQ>^YDQ)0{Yq|J@vSY_=ur>|=5It;xhy7v6+
zM>dC=Ku`$SvlF1kiV=ETSr`l}BBqW#8j=#2{Xi1Vh9cwY68}R#b~&HpH)J<;N*_La
zOT&NT73JRW{AOl$_Q`WTuj1nal9G}d{y5a{m-n@5+_=#iNzacT7fLJTxpvZP1KYhY
zAMt}1HtwnGd_UExQ>Qk+?Ow>q%OhusCNjpL1M~7c_EuXB3)1(EgqhtIe}vBU3KqkNP`z${k&&dcuto5>U><}{^z
zJxQEwz{Fd~>(4@DWo0E(cuk_PaTI6{dSQ^xebkj1F(#ezhwMN#)t_3&N&KbSEtBlm
zuQRMWUtxRS(ledaN1-Y~x1A;$@L=05{%*__F6}N>j6WJM_2+o-;DN_x-|f|YLx+zx
zW%#IrKNa=nU<=hWG@5_+<c3
zOEF*uEa8@tA-h6*Hf`^W2(y;e2yh;W`Wf_ebm{?tn##&2@rW;idm#kzP*y%>mVR`C
z^RlVWs!BmJt?2uW(ETpvKT7_0=c3QT*jkZBH8ND!Aq296Iw>j^<`#%4vqE%yjr)
zDphZuNO;9w$}nk_s~_qUrDovARb2@I8za^c`*hXV5k)>SH?b%58~f89TT7$cO*d?E12@&%+p%Bj1^;Ef`0q+rb^ljXAR3K4a@Dw7Q7b*Mb8>l(t=Rk*4n2?tg8=>
zh(Kqzyf(kBEn%d|zu@FVcWEhJP8~yy<|bqkHj3bX<$E-Uzw-P0fRf`P
zgGlYhF}J0WNHBDW`=-j9<5mir08tGavI}dS{)mWSLPD_9bc^}NAR_$|Gd$G9-p=~*
zOKxZRG?EzJT_|S?ABl{h%yby9)35W9w)&P52^UXBi_|sr+>nn2jkg4sC*`$q8Z`m>
z+(E$eh^
zT&wKXNU;CGdfMf!!u0eLU1J$a;f?kcmg9B4RDZ5rzaBDVM?7I@)jz!|{v!^InE6cG
zW9;|W*Sa2(CHncFI?0nsd=ZTekR11+XOp!5w7)eTSAVcs?6LBZWVMk(#B676R=?ct
zPHC~c^9(7c&&tOJ%i;kx(l^>!s*@SvoHi>-Jvldc!85K@0^KWKwcsigY;0_5eQ$s<
z2XQ3i1WD%@_uGA_HNGZEP{5W2Uf{flCT&j?X}z>u
zLr^{X2n-u$%*m=FDfGvWAH=Q85aXEmZt6ZwOjIHbe8h3-!hHu#P0iqe#o-E?IGHt4
z%@nD?p;BwhuP;QY{&>z8R;)mDB#?WcN)EAQo5G$_71%)(E$K5Pr~^5&TipbP>sQzt
zPrTz4*4O|
z(LEdpz#}N}+xJS5_vfUMSwSH}X5JJ=RK7e~&HU?IYEy^gID7ZjLm@Cyj{`vvPyyiV
zI;<#U9L~Q?&l5m3*j;Fv-y|+K_TzuwuY-kmr=ff+zqita9jm|!mL{v=uZ)<+zozAPjnyhJt}g$qZYsb1k!A&C)3PA|Ps(pk
z>d)fvfZ6GqQ~a7K3gE)6zP<40`AGLFU11)!)Ul3A?9bK$_30j)whg_
zH;RWJY5W4pYzL+ec9&b>8I$OX-Sz&5#Jh~l%vs=d&b_>z`}XZe5Q?Z_`)Z+*(K)n5
zQy3K;7b7V#?62DSU}CC92pL5|VWilpja)h;q`d#wwY7@9ysnQypW*b~==Mr1+}YuO
zi+FPPabHQlt?F^}z3o+48gb%}pFSl-c`Vm%m7wQtLWF{sEa{-?=O<_4HFxUny?cG$
z^$-GDATG(EP>bOdEf%qNh$S6%sOPx#=K5uOd*QXUwR7-P@Ra5wl}?wp7C(}3eQFG)
zAow*29+vg^yY2mD|2RYYDx&F@*ruNy$pob%wI(9=wVQe)POTh)(}+=f`SK;z-|x>E
zy1&05Dl;%QH_u8q8Kd^Ru7vBa#RorjNDY7p3Yj`I2c~+=(Y?~blLN!qnGxQ2dcShb
zTi^-eAj@t#2N8WHW75&BdYobx=jsyid(sv5c`?4HaB^_Y!03YeBg#jJgQ2dH!O!DM
zKi~EGBV=MWVFf)gIvoQd(%v~Yx(`e&5Rb+U9uDL=;e4-zgK;{=m`0J%_HsBU^G9!#
zQDnkYL`2erqm@9_64Z0MtCcL*@X}Cesf{`{)m=q&pMy&6T
zjHu|kdls^+;Uo9xoYOm#C21ghdWep!Xo!Jt8M{-~8=WVQNeYqA+uOFrakcioy|vnR+f-t(ANbQ?K?%{cnwpvd2;!R?X`EiR
zuHWvvKe2gq7Fh_QoA5^2S^f8a*z~`WpTAVKQZY^$?V8YNg{+yWm%$87KahpU@p|e-
z-QVZpYYaxQobaZ86
z&zNMv$j#2srYJTj|vGMp3lJOTqonJOKUca
z7ozby=0mUC($dZ1kq&e+L+=qrMx4Ao{k>3gem(iCkg()pDiCI^*
z2O;XN*()_KfF8cbh>2vxoW&8NTzE|>rGqd*olM4BlmaLqJEFkwRMA14?|&AO5Ed%A
znvxcpFqM7u7-7@a(N4+sRHdc!M~{l7sxp+C5vIc3p6*GM!s$?Ztsm83LoAaQ@blg?
zdo+$A4RiLDxs(%^!l_pb&@4z5=m*~6#iiA9K>>7vNLOkNKW0+&-=nb|&LqdNM-L!~
z9eqJmsgvQ-HoA?YJfI#wuJoN^A|t{aj+IUe$>e8VorWyf?>K@80u>~U>97af?!^DD
z=$~%-_*Y4?M*C9NNPKmN&D2$=e^zC8PDF}m
zxg7T&w57Tmiq*vv_>R?hVtsuxk>&)
zn$Qh-17AfPs{}dZcBlY4Ue$C+pnEHLcrwpD82j0oL4oI?LyUQRh9}1K^CK?#N087@
z5sd@pPGznDjMB)Yaqu*ieN9coYS>{XT9zC*HRIv6*oB@y_G
z{8o2PG==0nubk1D0Y>6S;J^p+pxw|G#{@}Ps
z-8PRL*m1tC&mA4lEY#{sO3v|cMRERTS*cq~L0i+s<{At5*M7-umcPAM-7UDzIn%g~~LQ4XZYh3;Z
z;yW{_Y9Xnv69l%Z^A95;a)}g4JycvuM|>}LDoB>c3vNJ6npILl3w9iegdPj;S|?UL
zP$;PGAwh=|B!eUVhx>Js@AXSuetm2EkQ5PNpwZ{^t2w!Pn6*Jc5(AsoktES!Zwd@d
zNSJ^}EN&CjzU6jjBQPd6I`b$(%~|L+w(gg)QRiqEq83l88L+Vqzv8RealdZW%7l+`
zB2s)ZCy#RD50W{bl=cs`w^{SU*4SwvR*L8IgQVrm_UV6AN1>>H)^xYM!vD`#{
zr0GwlDT7o;+Qxl+$
z$M71IDZt%G>^$x{GIpW2-9Pl3f}{YBAaN=onyLX^
zto<{kB~{!_DS45Qe(wM3lAf_Q!tK=Cr2Fp_C6iGqcrqyr*7~?(Wk26bsA!v&o_1V+
zFVTbqZR(nvD!l|R$id3}bEIOH+r2?j^bsS*2{E}KA%$!}
zHiL(L2HO#hd6F-<-Sx0X)V!e0)GZys8_$+{r6ehgko>6^6Ke7GO+R&2sp8Q~dU9uGBX`2#O@!rfx9&+Aik8{5jeg>GC&2B%*%vtG9Nn7t
zlLFhHrUZui=+3v7sj6IFOYt6EIS9^DSAkS-PW2PQ+ZX&1?fnIM^dTs5rnxYN|3ru?
zhX?y?qrxE#kdy}$6r7!znVD!8si&-Z3~^vmdf3Q(zxS_CQkLYnXfZk{D2T7ZVWg5y
zMMdSbfPf-IT}#!=bz0+55fO?EV2a<4*|zB`Qg3u;T|o+^8uq3Y_8d!?dEwJhL5NYz
z%*4;Rkw(3_7!h%r@N$Q9#d}4`_qnM`B0@qcs;a6J@0unRz=D1ldlV7DO4r+X-c@Od
zlSUBD`(AmzA-3a*&eN9q3@FO~8$;rZWTBOcuB)pAIR%B~XiP->YmtR;NCqA{n7*Y_
z;L-A=&C55R73yVO*j;mB7Z*JJx1an3+2I(5%J@l)?rXL=u?LV6TQv2JcM=-<0=T{{Vp)0_iH)8
zbWVNGu6({#*1tIOGrpQT)F|Ws9C)_1w>JTa@64NsA?xCXGBdPFqV%&OBXV7Xs<<%z
zH&eoyGudvo1=c}rz^XzpD(F|5zg%wN~kXas_DC%f5*x5^e#cTM@0orrNN{vwZ<84G%jb7>P{{Hkh>m6i~V;w
zQ5rCrNPPA9Q&1{{^=y2#IZkbvHl)#lg6OOx8nAbktp3$l{d3pq@xe~q#c7|XI9~Lr
z2yzJ7!v1-5RJs5<|4!g=6Sq6>`%f4=Pz%|@m$O9=@U6cX*ZF%B8F4yHh@wsocBpBZ
z8DTl(xp{ch$+!(9+DGSBr@G!5pcEjsXH8_Tj>;n#4&P9>NM}PK8QaO(O?<@@xj;$+
zen-0t_DD~kURh72pb_ZmAQz7fv{y!mQuq8FB>h66;sv7z(F7x)Bx7$Eu}~wMfji=Fsu?53iEt!aQ9697xJwso&(Z0OZ_XrO
z#^`X@W>8jGl_G@vovU1HFrbHbT%p+jhqtOcb51P@Kq+WrFd~$d9{!5dXz|Z
zu({|Nc4x^GIMP8>SwB&_%a^K8##iIvwtrSzMXPDi&q^%aM*BXk3eNxBj_;0bl~{LjLSj4x6_SwQMe6EqW#!)YTGmyD
z*U3c;#nK__+5NNF_T`JBl@%W}WRW@`6mI&4uV^#q=;~f34TQuRX;&~spIT4_-EUV+
zWoW)B?6Y@%eKeUS8gn!AZyT4<9~^$<=Q+i#V=KT+YtPK?+sd!{Z0U36VLAka+~<
zg?dn@f{VP_Hv68IO5RT8?P2UX5M_o`*05P|jvH8QXee?@j-u(IP^#rnzqjWz^XnG@
z3IfA$^#s1q9|bS+ii%XBP*adhNB2+9A?u;84xdtG?cmldS#D8_gQIXSEBlU(pibxO
zv^hR^9-TN2n<0BtL!AIUb=XSuWyAAdtiL}MwQLFTjfp4tZ+5!4qBBo7f2DeSpE@e@
zgrFb=7}XZ2^Ex(BN%jtlNJ=t8Q#0d&iX^xQ-n9&D;6K%Y`-a87o}{;3qniRgN~^tU
zDfxF6U{F_Q6e|Y=+}-2~RN}WK=)7}7x50af=WHqoK8kolU9m#kJNbjOj_S_!%ri5q
z%TU+MfEHkLbF*XhsD8DZ&>M|pJ%==i@L8cuMvO|4b|E(P-OB*8VL8x`@g8?2*;!?HX_Hm43PDe$rZLVzP_db%db`-tk=Rh=G2Az#aMnca=%nPCM-b@Y
z41GK}wo}>d{XYJ-vExq9&{4?c84+KdApkm!yV|oCb7#ruZ`mj=xNf6A?)W<#ts)uU^IN8G~Qww>wL!IL^cuE7I_GX^_|8ubl6t3dZe*X(;Dd{H2C*uDcdST59(1zoMm?M7!NU&j^oH7vj~}b*SN%h*S=>
zf4h6Hp-DlK5_b!JVqvHg2gd(-b>8`|Sfu!KlEpui^V?M?JuwIsaI5egW5i^hE_uP2
zTVvmG?ti!RqdAbd5|B_UySZJkuiN?Qh{UZ19~lG5Iaa8uy`^UCk3ZKk%{?;pH1d12
z)ZqkBm2ZO1)k^KMuSm%a7`#}jKlDRd*}wCPTcGaF{1`FUJ$u}CG+BxyK21wC)T2_0
z2OlaN#veh1e%H#X(W=zEzW~>osr5DSUX+0NmEg9iuLLUHcA#jo>^aadQ`}WJsfT3>
z{nbTvaz^35#n-jI)Le2gFd%?*ZPCKLQ>y9EFv$@wX?H>3P3u9-?P;-E(pHK%g@)=
z*FQdQr5y6Fsbw?eiX@a97p_VL_|C{96Gb>U=e|&eO
zsgs)W|72+_r>S2Ku;4jHPmgqMphuZ=
zIsMnqpFEP1lD%)Rh_g|aW03t%E--z6@ywv+Q*&syq5{Nbg%!?9rxWc9QS^#x-Ol
zlu?xrO|@&y_n$L!q{h4+${Ly}ncEk8K*1f?O38MC5=nQ2g~NsyTYi+;iaz}D!E1NK
zZ+CI<{XGgDT@?)ts`7`!VSf96Lr3rP@Yo~`^EP7cc63_kfn5#^G0HGJ%OtZk;_^bu
zi3e)!HJPp#-;$47{Jdv-q`5J}
z)3INr(a{cVweRwi)clXMDxXIYDyV_ej}LG$&a@t{PR4+-_mZjSqYYf;f~De0gH`O;
zY^i`idj>?iKUeCA*gbDEDkt9dyL@6Fmut36)0JFFSF8Z%qdduTMx4}wof&LM
zZ&B?$U4GX-F2^Ost)X0i$X5o=K1^+rd-SU5Vkyx_di$wf$5z>@YJek{qWw4layQQs
z#^xf=K*p{hS`yL-oN@nTs^Rjw@((Vkq`RH2NaNpTtjz#5{_qhhoF>!VAQVJ5+(|p?0H7a;YD*NrXq~N9U|U
zi`;PwRLYp>^NMp*PYCJ6emXW!3_Tc+JR6YMVXZGlBsJ(gUkr+hgjK;#>%wP-E$c5ZLJ>ZbER*{?P>s?$!!Q!^F6fRYV?{{@2H{Z*V;iJlqo8Apt*C
zlshvB0y@5skCNh~wd?nJ$tQ`0cy-9T`6JI&&rMq9JNr0=fdqLM#DJ4^u%#+g++hS)
zf5WqXgmhL
zx*qD$3!(`~1-mH7$-62XOd;)_gvzj;b3F`jL63Sm?QqW;7vFb8^h=|wbpzXZCkol{
zPWVFAL_(tph9%xPO{=m?IijFD^f8ctgmrWoMgtl@MbbCfoQBk2#71dJ&s3*US8vXX
z$vGVcDL5dPX#4s#q{FJI;7)0Y{V0ZbSN6ZO;PjNhR(l4L%aVgTUaoxO4sJ%@z#u9n
zCZ@p>jOD}Oa*4rTFvQVvh08Yt=oI%X`2Q4_^Ip`^P;m8}Ec!s8Ml3+1;A%oX&ijm9
z=f-r?tYYEE$HnCrsyI;JUPTkchukRl34H3|8-^K^OE8ON5#w91gbtXZ_}&WR>jQ^S
zWK81S-OtBfbaagW<^R71;SrCH>~UCb0lpTQcj7_GpfFW~)xA&CU7vrbn&v-WRMk7A
z)37HAu>D@s5P0jj)Hdq0Sw_0<7o
zVk#qs8<>MS=ec$lh9RrQTr`MwpwhS!qMJ&|Z|
z-_;5DN0C9;Bep}N8Y&>jalv}wq?Zt?9g;XX&y<+4hiV5mPJ1#`9<+^SpLUhJZ=OES
zg6QwDJf>7`cT~yn^@1)}q^*t96gU!5dcJEDcOde<9@gmsefY>f=lj@FU9bTTZ!QU<
zeQ{=9v4`+$H>-2S%Za<8MYcHoQ?@l{Nw;rwiy?;y=87X(VC*7#J
z6Fd~>PM{=ho8>md+L8rqzmHad-}g$HUi;K^^RUTiD)@yJzrW7%@H5%^8+u
zC{$u=vv$78w
z=>!U830o=yFc}S38s|nHNMQW}6^ier&m%RXwuwBOtoldQY4Tkn=et~yc;T+2u7pkR
zfYa8slH`m4Mhqvyf6Agfs&J$bN#}Z$C)I$CG(YAa7nT(l*$(
zYyZhq2MDPWw?7Foy!XHP1+U$h|JHS^=|4uNE&GglzrC*9JiYbhq#n@-r?7AqKP>U|
ze=Yrd-bV%pVchjhZv`(;rm7M7Mmd}FjhQkmhGA{0S}}sSb7x?)_IoRre+;v!4_%t?
zpF96clQP9TxE}hGFdy+>`;%?Mw=4+aNr#)2w8~lpf8zW{Q`!*f1g^$on=s7g?OQX9
ze}eIuHFZ|-aImi~AG^5Z!Fqimc210}`KTJ-7P#an<+wKUxnUYcOQAF+Ip_e**ven8
zyi35H!l?S&Zp+!9Kd(CZ)&95AeA{l(3T4Tm;Ys5@uIGHu{W`uxv_g;{Yk#ZVjW_==
zg0DjI&)=sWjDx#EDdi?5Xw6Rna_ZRURy4)D=HGj)j|pFplNra~RGQ$m_dM@A|g6;KVm8i}d;H`r|~9^|U~_e_l2*
zF_C$D>&r92cYvU991DV}B5hA?6=$DejfnMYQIalBEUA-Moj^PE;Ba`zaUv+X8{;%v
zzqz@b3yKc|nPOd=NSWoXmLlnRymyQM+IsYTX*LaFDT4D{y(hI*zJ`t9P4IxRKtjQ`
z^u77YR*SaO{N_2N$PB;R>sszUzf`O;z=3Gob-j58iHrmeEi88OWT=}JK?Ih|D&
z{lu6MgiHkq;qi-0?#y%j=0{`vb>BKS%Y<+I!U2%o`No%!8Up;HxCN
zz+qm;r6B${i@4e|_fDd3tZVwNH6~4_eI5G_X_B5l0Y!!OJ4*o$Wm2fqppWdxO9pJ5
zpLt~_$vtMJ`1y|Y_2J9N7AUP1PjgyFFUC26nXOHTRa&};1jy&p^^2$ds`4-X(7M^8
z5GR2!9xt_c?Sw=N6tzb${m{~zyIgoJ+V$2C+_{gW+**5u^s>cKa0s<6fJU9Z_ly(g
z4ua2APh6f=VH`lm#lbn$T5b8vJxEy}1Z5|_EXkt+_D?+Q;qgQ#M
znW&eeo1!!JWiOz?)pBjX0&%r>&t1{%YwR9+B$X`PuR+#^wfF^6g=5c8>{|n)grc}f
zH~A}xQcs^(c_chPjyNh7HF&`%L
zs*1+w`#jUc-iGdG2rd(YdKM%1mBoRnIag%_JLI4zPgoFqn})ytYIPi;O(f(4mnl9+wHK20acGQxNT`QF<+ZOhS5Zh^lL
z@tm@IO+w~Ie{6cARFNGeC2PWd@n^E{xKKx!(5}F{9r;-`jkVa5lyeE$)#^UyEkc3i
zwAO+$nEY?WGo~v))kdlA<}Tv+q|itB?Y1Y{`(Q*w`eRVhgZ>75sa`i6yFfrm6j8=E
zH;gj~6!g4g>hf16@HwXT7X%9d=t
zJ6&&s8LxP9)@mS|i!NHrjG`Zx?mpFQ^@4HmJ-pq9r(~XOnZD#}T=8{oaBf|A&vTOu
zqL1Tl#?A%FJvIkJORbuvR~-ZqHl+&48Ap@cX8)*?6UOCNFUD>>WigtqH~+%!vcVG~
zQ_iBs0e)MH(G<5pkg{y0-8j=vnR{91&X9BNsNXQ^>*k#2BT~_3h+#`1R>6CvzfSea
zH;0fpmEA9nPk2;GIh$9rBCIdDb7Hac+f?2ow>?i2K2)GXyX-rAgD#Bw`)U4ZOvW}^
z=z_g`_(KwriC?#RL3%@3~0$Q`VLC
zaC_#GQDxhFFYlQda^oIxWZIj{Z=0K)nk0;hE$l(hn@R8ST@$l3XbxpUR{^R|gXIs|
zpwlh~7z806DzFI`a}1jv-zpn)tt3xdF8(ArUezT!6vQIwdO7Z!n~(8SvX+~%s!xm|
zE2jp%+O=}C)x#coQzl*W^+ZOGw-e$f`vVon(Hh6j_e?*2#z_Mzn(hAnqs`>dgcf8e
z!xau!^z@jZS%HUwG^!C`NxzByBDk*sWB1Fz_v51YjEsy@zJJ9r?#rBh?ZGEQ6k^+>
z4HB8%@`F68S(wDaeIj^NNb%zMCbWsZ@)5b^A6vavAd|<;oN5AvjFW>JY7t^re6e~x
zgd=&(CU}EquLh^FCEaew@`A~}ToXKv_5H<7e5O+TWR<`^i=cppNHI-zS>5QY92=+X
z4@~@ZE{%rng6lILVygGt8yX7k&V^-()O$vMzT3E+!6CZ+T5AdJGh;C`gcQdX-bDDu86WJ=YZKKpZn1cU
z2B>StGdKi|vQ}+HL$+88kSzcc2DVa%XICVII-cTwjg9T4iM~{lCWlIDn6w(&)0;hqG>L
zl|!pf`9zyY&S{0r(U?yjVntsY%>H8i<@*HF`}S9J>x$p%9SZrH-a5s2XgVq9p<{X&
z+{gVf=GExbIad$;1b%mQL#Zv1YH%n}W@m^fN8oV6ZZv6LPvoVV-_xI})4v`SIh<=Na
zlbu>^O}j@tG>xT@?h3s>xFhKA@}NTw6`z(C2Q>&(K+l}AQHqmrN2L+qW&q2?{$J|E
zE7BAb_h*vmnX>%(tbZTaxBIcF
z0gERmM~{q*XeDQ|A?O%rkNwwKkiWdwOz+NJky3o+i9kOY7#TAvE05b(utB}N8H(}0
z7sVd_5U*MLB1o%fHs4>Urmi0OdB+7%xI-|A@&7E}A_=|SU=8|c{=_F3Fnk9MV|~dpcRS^F&mw8e!h+pdgAb?sJILDK-)m1lP5$M}7mLxS
z@dz9V)SScF*?9#(Q?#Jy@9b}C85m%OPeYX@5+)rpP3!#twc0Wcz%^LXjWBc(vY9D(
z$kfkd2c}DI3#2lKvx%lr0Eprap@l{@Jzz86*;wR@_r>k+Zbbul761&*u#~cu)f+8c
zi7+9mZ`-TW&0zmRXccP?4|b829z6_J+yGoFf+r%&={XZZI|Tpwlp4Vdb8_q%!EZ|3oe+)W)_MgR;TP%Z$9M#6*N`5zwWOA4L@
z1lX}^x5ZmF&YIS-*YA!ln^o-EXAMDJ>3GFy
zGy*R!0qlgaNY>C30Ax29~l`Xwj?4_1)O@B{(i
z0`v^QLPKZDYGLplZ@&ilXPDNCMlf)INr4EE=jZE0lSTj~Ta3o&iTypCd&4g8$9}le
z>3_~=dqrRJ6rYFuh8||K-*6Rbi5EN$I{gpC?CbacSeaG>a#UucOKI@KgIgbx25d4m
zls-6u;aLI%R^X
zIWB&}NQ(kfUTYD+xGkal-XzWH7s29*B-^4$$NV7XtMUfbAr5CH28yefFoNf_e?Gdj
z)Bt#}9=#571DAaqTJ7r#JQZzZP@v)N`HX@BN-hQ;K%6Gss;S&>#ey}MfNgC8+@I6lwiBRx
ze0F9NF_i1ip{r>2`N=5{z$ZP9kEfmgN%&ggnmb!4Ek1(pK~Twv2EVh{v}^Xb`HEr8
zYhDNj_h`&Zhm;43jGIHq=s_6~I3oiY8QBzw49#Bm`q!B^e?h7J;qw8A@>c%&VM`#^
z5^Dj%X0bByk%yoE$_Vh&0?;L3H`I)T69G>C2+Rc&K#-1sY&-%0Mt`aGjiztT^*@tr
zLz*T@_r?!E=8&jv765k4zJE995T(80|MylwGW~nN?USM`U@|4Ja|!|Axt3Z4Ay5i{
zxgs5m{|RmVzhYiB8M1%;Fc^XCy-U2Xdl78fHh|V}USg4gnbE
zH~|P`c~b6>N`Nh@7#XEh0`|N)h=@X2SNB=v>~q6U>6e3*VMMbTFn>K*r~Lf<=1#vq
zVhF%Y)Nw*#6;P0*iS4J(pFgh*>~02png|GzAao=7R3EXQA*NVL4;VEg9!OCGs{ZAw^NQetjitxG_I)gMpbDoXTgBBf2qDWmUmxI4~FLc?j-
zoOA6i%=BiSs5BDpE;N*cz-a}cesa?*^vwec(4r}Te?Oo-&lz_GEH62C2te8jDiN6)
zn1ATF0MK6j8o=!V&%d;{tAPF}ySlmtcNKz~hs?p*(}2J0zUy7=gJTS+^;rOBvW_m%
zY=BFNO18(S>)gIg!5xB3bu|O>hiPybAlp=DGj_o_MKwn=!#V7Def{Y(XUH}(-feBq
ziukig^?NVZ@}tCAu($n$+RJ+)tUd)NWNeWurHKPg~M3^fkps4881sHlFF>L@osPQ+%P}?5*%_A
zSVHki?$0o5v;xV2#nAhEC=`6i09?eI%WMLIf~%&7o2Jm<2+BOY;vsW+)BMMW2N7{`
z4sL1itpjB2EHIh8YScc4UHW`UErdc+uD=;8I6Zg5S1wgXmuqPrcSCX-M9alR
zlhZ@iz}yQX5a;{=TRAZN3%va%CsIL$vS~AXi2y;o<=+;OO4XD(+*|jU3}r^4*e=)w
z(BdN=4+0j$inkL
z5-{lvB1GWMTi0gW@leQ@GeA9707hvEC0O`ChI>{q%`d&I!Uy#O^D40a3ClN5@@hu$
zk}Cq7R9qt%FrCn-Tm61bBm5-7Ymi*mL1K2Q_ctiw*#_3+=OMMyjOr{F*2f
zGS~;7s;~sMz!Y&Sa{EGQr(*A@$ejG)4h>4&s|I{93T5KAGix>MJ!@Os3$MN~R?B8m
z%nY+zp8zI`hXO5Zs6X6y0No%&Y2CUN3@)9-#C?0k94RG&Rj7o)YJ7UaI8NCp@CMVX
ztz1{Y3OGT@l=$IsSueDi>K@GCX`7gEI`3A)mn9Ig;ez1y*+vji5)|@T6xh00f7o4=
zPklNvMKM&&__!OA{2+0Nu6+IsS-&`Ej_48Ohr`vD455LiW2O)JLSspv@VqW-F=M?&c9{^6
z*XFZ~43}yHUn!8%&BqJmW`8EKs|^ahLR|5()NF@OP4?Fu_?oDdPzN6e_~!{
z@c14}(Y&(^r>QDKz3jl+OY`Gr6R%-Lg|oV`8RfhCuNK{+r=`S`HKwEvD)e}kHOF(J
zX3%nJcZvrV+^RFP7%WPMvkni{#~PR=@kFP0qC^u=@JC%mQ@KddEa-m#dOdIJ
literal 0
HcmV?d00001
diff --git a/imgs/plot_phi.png b/imgs/plot_phi.png
new file mode 100644
index 0000000000000000000000000000000000000000..43fa673d2ac93c7382fbb00682b2d3c47f90e9a5
GIT binary patch
literal 20450
zcmb5W1z43`w>G*!3=|9mkq}T8(v75mgu)`EyOr)v2@3>~MGI2W4T6-kB7$@y4T5y1
zg1{LM@ArNG-uplMI@ftGUzBG(YtAvp9C6PPGf?rN6ybUD^C%RG5F;(Aj6&fOp-?!x
z=gz}KR>f-a|o
z#HG4vTIZ{or19jwaB*ERr8H4&VmS9~Qj%Wm)~jdydCyAtx$wWnym*1GNgrG+C|K$9
zIL<1N5gbdy2NZF4QiS*&B14c(LYNu}fL7a&=><(?^|x
z*SdLV8mF;MWBv@PIbo|))Zk!m?XA&uCPv2Yv!cJn#ypqunV6Y{mJdzaVv_mB-5<+g
zt*-lr2+Uq*)vHu(jb@D)j62$EUMZgtndKe#m`f1cYU8FKP28y18}}-;8e|)D>J%Ot
z@mYUPJu_2;(?yF@(Xm&%CtS1fjbiwN(ev>qTk_s}RTJ$C(;|DLPTEC}UJm5zWa`*{
zqNB%pe=z2X-D<)M6~)7^PhVVGa@id+vl=P085-GL9m{QHuhT^qSWvmhB{>ZZnddf9YP9b@ecmAn^A=flov1Uf3p7ZL2NG3tKgyN;dcCK?t0d9
zM1Hb&?%Y}UGnX{@<>iHs4a0S*=sfJsyF<7HeMSL2+Y+TAJid?
z72fWW@SKe`BMy}GC+Owj~NeUQ~jNPOWL-Xzzcfp5)2;T=_0
zgIe7+-{S+9_B(^dSUn(7aDP6PUh_0g|q
zjj_hBX!G71hb|Q3*{XkTPd!VbLAQsNz3*6x;da!SM$5adK`!rA&B<
z^$?ey_f9`hwa}YMz3VO^A+fte#%;9%voxbX_vahePYV*u!dU$+NFVM|4OY&BqI2hmCu0M5Y
zK1xGhV0rq@z~k<|KAw52x{1QYMa!tjNCxsZrr6kw*n|Vl=%)PK+}`mTua6UJp#LIC
z3Z(z9vh?#Nn!&44PtA4e+TxXs0f-<7oLaxhQ7MTedcdQ1aWeRU;)p(ZL*6gUqW@LBz?0$Uik728v
zzjgSNw61F*GpRu9qbVVyUmBGO{
zOsC#b>D5ZAg2viWXJOmQYu~*GxvB8*!Q31D-#?d2aZtQ(rO_p3om6E0$>rmoxn1qw
z@sm2|-dTk$dkZHkN`EwM51&};Y@NP}M%B+EmMCGp?!rJ{DqQDm&f=hWu*tMaL
zQq;%4n|2m+3q$^J?yXG}7C!mKJY3~g`0Vq0MA&NPgxAO%^?iNo5TUUdt5k#GemutZ
zC%BDgZ*FcjP{qmZ?!O^>;N?|SW<8X0mdgEJc!}|MqA`c28;kSvit_UEPj$|sBG`$z
z`z!MK=CD=sKOA*V)~r9BiB6c@%`#5??3ruCnQB<=6-qti9cqKTq^7HM`Q(TAPODP$K|
z^jPznIB*_w8C-yUN~o-=a#_kRUo0Ke?*;drXUvriAi;9H+GD?w-6xq;-(3t*Z8YH@
z5hr=tp=ythe|_q%vFT>j){1z~&{B7Y)5$u+02c-JN_rov+@;dEOEn1W&zzk#ec52-
zQKRosGrsHu_O*3-BdfezLu6u7^yymlR)~6`e!e+R4)}Okwoj6X*QN5_R!Uc9L{^EI
zG!c-yrA*PoY$}hLEAk0~2u~Md3g!$Igm)r&ZHCi;1XkS&MV~qY#Bn-*!KWt>WCs&p
zSg0^IHWr>{&z}c6@^s>P81@De750(`I}0F$wiRK@k+XqCm47!|5r@~jwzri2R=>w(
zP%jf)kuqq{AjvQSenvao?tfi&p8fav9H{%{h%+}n-
z!66ef#mMwSpbx-X7lJHct6KMG`RcNxm^CcNDqTF+Kc960lg@CRpk;2$0negO=pvFyfaOGteKFI5S|7G*aQoMs_5zINkoVGY(^yY
zMYHN`yZP?9g{OV&5R?Uwc^Jf8(1olhU#EQ0pTf3_%q*UBhWZ_XSxQ}3
zv_YW~eD)XiN9qnojFP4iS^hMan47JVDH9I#XMo@^dHpSFqG}!5sfKM2Pa`#Wq7=Rj
z5L640)}Q4!tm#uod|TylX`=YpyuR^%GH{@wW>KO0}~fC<-Nuh|Yq
z@ORSj-#Nra`QP<7(P|KJ{bnvq6pM0N*R9LQxGkR~;@#{`<+aQSX8@jqt)s30k+e4y
ze2d&kf7*<&^0!CM;$yTwn!Z5*3}Wqa^78tIOReZI%<6evq_yj(6P$Z6`{lj|d^Y3N
zTKa_=pi>jSzww>ffdk;}K_#gh>A-r>A=&UmpM488b#)nX9qdp+D;@bnOS2Lu#Z
zD>Nr3daVk8ebMure0uf_A(?c~?#jp=-=im(9`GHx;5%iTD4p6un|=VOFzrmtmu0U}
zb((EEF}fHgc6JR}S=rs{UIkG-(chmh^^K1g3=eA`E}Mx`Z1X=}2&y`H^}G-3Z4&wwsp9;N%*_0IINiHDR?)MaNrh!)WtJs8;slh2gYN#DSM7)aCS
zji14OF(Q_K^WYFb^AL^0)JxLEjV9_rWW1?MzM2{a-lYYfX36y&8c9G|EnPLu?$-YO
z#lG~7!mjVz8wwJFGdpkWo=E81RX&(uEgwxm4v3_>_D>7yjy&_>L@6815{72g{6JHM
z`^z2#x`{j=d!}qEI`DCLRl$A7!5w?n-K_x~I};7I9}c1cf%f1qcIy(`))Yk)Z&BP+
z^7OXC`qZUYb+w}?JFiGkphB%KiT0+fN-KkGJZ6@`Gj!jUMls0bk2&0N``X_h$Kb4j
zlA?u2f*X=D2~*odjpY+>h%4?qZgA$6Q?2!NJ4b{H>Nr|GRdd;Ht;{#{+)Hr)KJ;>i
zqb#k;;J3w-6~_~w3&2YcWi=`bo=&+uz<$ls`iKCnm7$OHI!~Gxes27JFRuR3Y;Fnv
zrMQHp)eMms?!9h{C7C_x>a`!|4<|ZAVi{%q2f15n&NqQIpHIcAKYsjpJbGQe?Ay0*
z)hmaJ=WMGIHw?{B!Pi!TTBpN01dC)2HzHljhkfS^>I)DeA6jYX-XRf5aCffs3%xpv
zLAl-!MJYOKR}sDArF-xH&X1{;K&mYQwH!o)^M1a>^_!Z;jnJm)?q@?Lo?dPS$#Ho6
z*MKf1Ogc-c3%#nwueWGc{v6dRoyUtSd>2N!>SK{vWHjuQF}=!qYM$BS&xX>&)9RLu
z+)u@vi}uB3V9Gyr2c4nzha?28re@nryRcKYZ9d!g*##x@+@DG2PYmA00Jn2~DLsOU
z5-GJz@1@)NbgB;2^rm{L%XcgXN*JZAUvzYwNZ_HJE!_<^EhwK?JV8rKII67zeuOh?)Q?m*@*L7bl+*w$Q6Z-km
zoM_S<;PCzQQs#n6#v`$EHI3ghdmiVA>hUh4QN&RLMxN+uD(qtJDa%!~G-?{*8P6?F
zm72rZ_kq!5Y9=B<$~af=Ay2tR0?cx>saUSBGhI_Cay)gN9!@|h-CP{bN_=uy$*wiE
zmT;omVz}2Yoje~9b__jtx`#TT$lX$T1dsm)vgzh4a?R@*8WxmH^Lx_aDfrLuI1tQ!
zKNZNcP{HM@PQ2Ol3?%fbI5MMl)z5dlZIZH>^MB4hlYAyXRPT-u{idwo>xV5`v?tmx
z@D_;fdDsvyg-1WqlX{*~;m^V<8jdPZy3DHj`OzIQY1=_>BBS&Z!alM3baD7k%!M^c
zVDUKgVcLQyIWo)90cOl_mFNvBLC=Vxk!x*i-xHTM+=;inNK&P6uao-oqKIGGqigQF
zH;Rn}P&3W!ah2c%pD_E}rR7v5B|%3iYmP}H^MCQ;!4$~8
zoRbRAVYuO;DIufqHMlVw9-lpVycfqL+V=MxMy?F+WP2GTQ!-i}Tp+ThyNGXr~5fJJB7rYvhO19cJz`I}mt6*5J2vKZP+WL|Hzjv1KV&Gc
zQv)R$6jp?R!Iy=J>MLNME40zPWt%_95MEqD-ru?3q)sJ7dwrcD`a|R+W(Pw#WGYtf
z63zjtz^;2rIsEst{Wxz*qwGs-&1LjEpKcDO;ulBp|3g7^_RXmc%G$fG4>F{lMN!oM
zIGMtj@(w!f`R)k9jH?IV+Plth+(!67vZ_ED
z#dP7#S}L2~?YL*E^_mD7>Sa=oC?zO}bOKME?4@?iDPf}NSZl`g*UCkIubz-GcwS0n
z^Wg60!EKI{t@>xU5KShEIZ0VY;h}IUd_!x5A*%#*0hl8jGiW4t>u
z{~8fwe-ON=O3$m|JjvgqW0MI%cLfplEE3Su$>&MklK!`p^>xQxBr>@gV@=kwE6o91
zM0T6m#lKgA;;MEjXBX`i2+v%PqJ8$YBSEMXf;|v)(1gfmF9agEe||!O4qoG1RQY;U
zuAAiN>RkmJFdpbAyl=;xX70k96jc7aew_A_p_&klK%$J9Nj}=GZ`lWP)N=iWJ9*^u
z^=q`@+%ttPzw9^Wyn9>q+Gg|8j@Dvp
z0H#F8;n7}QA%fSP<~kB))Nwz=b=WKaqPu$Gin8agGa#7r@a=z}@T9R$pb^Z7GF+W{t2<;3H)@D?#Rx(rY
zf3w<}Z4Y{oDqyR%3dfae|0n2STM-@YDzQa@J6vNm_-d{#Coz4jYL~#PP0qbPfya99
z;|CqDkQZ?P3dVQfSj0e&eSARE-PQFT&bF}`j13&d?=rM?bo#rxLWMK)^X04D)&?Nh
z!-`p4SO_!GT7lW6urS|e=7`H<^1tp+&ksj
zekG^bV0V2FkMFOVj*IIg+HW}M9I}kA4;xo3XruS^I^BPWNcRm4SWbO;i6CwyM1cTs
z_%$ax`+Fky<8K5t`-{0lYil)o`uA^L8FTKHe`(I;hikR?ia9DgL+xQ+fHHZE_4f?D
z`YdZvDe-VP=A;Ht`)Ix7I|Y|pxcu~4VDHAsT~h0U6%ZXkGGPJpVY?T)qmS9OP>>x&*kZ(iu>`H
z@&`_zs!1%Z2-DMN_L5jURB%(34QIdVZ{qWEPyU$I^C{krn1Kgbjw3YiFzJs
zV0Hm%Jw}W*Q*3(D#}M}rNl)|VG&q!0h(b0ShA=AtOI1i1sn%0CyGws(pt`R@CfayHe=6>kusqeug~S
zEoCogSLD9S*4bhdO`q>0w>7|oiT0nhaKg(cJeCO!>pW9dMv$L#4TCe*c3UP&Qqdo7t956JxjW
zyHpDLN%Zfe7glZ@hNl_JSgt(Ga?r8c{hji-6`zQ~6GYxVel{d#k={X5A{^HT`xD>w
z81IC44ACap>D-*z0Q?!}x0n(Q`_I@jAcD8BD9;ovZQ&_zubjj;xyKX@=a>Lky?x{?
zMt?u{3gz|T4=ey{o+AvC*{bj9{GLP`c2p5e4}$s`*@1LLv0MC?(T`x8&L9dGIU7=u
z!~Ep@>DznL%oOOKc=(fs2-A9ct|V=fP=phoPqR{&TnPL3+a15tE?G${bKC=?b0Uz^
z@1ln*d(~L!br^)Co_&x0XPVS}p#oKz9QJ#m&D+xd`FgN5Vc_C-)(3PLp9b2(05yvG
zWP~Q#DsOijPxM_X>@y#BP?MXKMK)syUUK{;1NqVHoE#ae!GiwLPflkOw!p>dBB6_!
z9nUD4bmt+Hw+j(Te~pdE324~g3H{`%?&*{96-yMoz79Z+^2`68Oqh>vd`-I4Bg+m+
z3?)0eLbgvcLEnsGZ=4}0-fU)=^3@ZqrQ-ywWD*=XtT
zXcny$R-Lk})YQja9*{jT>RJ`}cP+rC`aC4Xxo+NI01pSxOt@ozh=~J^SE|
zl~hrQfV8y9b$|4C4gv_}W@l$(0A4-s@PkHYy``te>nU;G6yGET!YQnNrZY$&)8$C^X!G9{q
z_OjT2eixZ-3({U0D)Oy<>ma$jKP>k@VHGD85~E?On8Za&0FXFshrk>qB&qD=_Jk4m
zG+IeJ_~S>+N7CA^f6Zu*@W%2dleJ8U)i-GFl;`ea%@I4k{M02U7qPZk6*Gb~y-2p|JD
zFra>ihbP!Q2@<38X%P>5)_nG=IM3pvc8ULSuDgjZ(2S!=MF@9oi{9H1innk6@C{;p9}V1`o#z28_Cq3oFWcsHsHm_v
z2BU0}t@MG{I7jDOwk^K#Q0>1^8U)D~g(m)lKve@0%f<8rkuYz+WPc(?pWqzsv)`Mv
z2+{fQ&bbH^EhHuA0{K9aieD(lbT1ldIY#EuCUw)ZUFLW1EL34}U4Q*waX#;@c0L8u
z2$WoL;w#3wDaLu~^)@FJSzLZUX~Oy|p8MfBEZpZS^Mrb7#aQttQHVj94~)rSVYQ~W
zgoBIC;1im&0IB7)SsFA$G``osv_{
z=y_nYa;4NWumdNAUwfTN{3cUX=}`pDm31Z~Z7vl8I}|_xzxL9AtKr7E$zVOA@121M
z`QY-cq$5YtrJey8r9r3xN4*gCK{iynkP?nu(%`T@lpYG>!)YpD`3f4rvWQbe&4nx{_aiuUlO*H-gc9<@B||6z$Nbi~j^@ht$B8qcfP@WoZ1f;%X3bM_
z6LS+hU@`K=)-qut9CjD9-w6l$2ci@)^>KQPMx@KA^OS~#ZuBV1Gzb_C%lIb8aY{0T
z>X#a6vEZCOnG>GBlksZDhV}>(naTu#jYk1*Y)W0|c>1G-i59RO(tRxUvVLC>W~3j*
zWZTzelT(}toV*u0PbZB68h85!67=`Ata_6PF+q~C_W0nDB_7m
zkeEgmE{HI_wI8yJjiQE7%8r08NKu
z4*f$VQsOcS6=_nq;DkRU!^zZC=AF991XbW=35tLA!jj^L$moDD5sZdx7vahpW-#=_
z7G(j@!MbIJPLNV@xvwq@Q|
zfXn>xP$4NBfRRXp55Oj_epau8So)_wA#ragBpjB;)Ve%EaWj2RT3Xr7ptznD=RBfC
zuAa9r@eH?-pgq)!A
zw+Z?}vp)Rk@<+=va`%zlmi}h}G71KjYz>jhX4}#nB~Kp#3^rtrrNXZo3=8cM1Mke>
zib!TSDa-uCgihA(EFa6t3#g9BMR8SjAN=_RNEb;SnF5c$a*+jKXI9?_;&erf7D^UouG3v-F^
zoK|3i12WPZI6P^7X?gg_`PrhQXfChf2_X%bUq#$G%smG-n$c?W8K~2Uno$T}WWpvW
z=}oEJwO}MlmE#S95T>8aHbG>QtIQ32M&Xj1MR>%G$lmdC2@|PM4$$6Iiah_%gUSSc
zsg(FbUd7S)E!)$RGDwMmrl)jjw(xCinlulb;k9nU6XfV;DBg2*y#nYn6=kX-;1^JO8&LjoO3NUMLQO0L#qQpKZvoD?vq(Di?rfl>)enok1G?GH;*7-h5)by1bKek&it)lY1%CcrGRa9tP-iw}j$YHIX!%B+!0J{kiFq;5!=_H}o2&R2iz6w#rhTljJZ>P4y{U!$Qhm{v|&
z|CPT;>ef5yY}lOLWh^HAsH8^W&4?#bVxg0N)m~_%=@%p=`v(VGIeMlu{!_#Np&YVZ
z5eh!rJ|v$Am08$79=vwrdFDLtB6#OiAY+pp7s~xl_=r1p@F#OZ&RErOX=VfM_uP)d
zcA0G|evw?yXxFu2?w^O3-0pf?558CxN}R@52`5;w87>TeivO0I*|OBE>FCOukka48
zp|`ySX5A_G1w4?1|AyyE>EL$0oww}Xgzs_1svFb`%w*sH)cfCnmsxRPdOCMT>X%=O
za{hEKR^#0t&*qOZu(Nlc#ITMVsQIvywU0$@wW5Rh_Lwxa0f(t#!ZQ@uHN1D4_Y%(e
zbS}5rKqXPxnuN8vxgwO0$#IX#oBunQ#xXMrMVk`vf9h)0nh1gci{s8uo&%e)DSdul
z_@!-2g~wg0D1U;T`zkMO&vPt5!0>;8T96Yq^M>NTzVf<58>wwJ)O8j>B-htlrXBP{
zLwh2aV@aPrGCrWWiGj+3ZX^nZfb*#U>H(YbvD+A=v=gcJ#PqaacrbQmo`0e^79Gl4-%4H(rBED3v9~)+>YRj3Nt|
zzw;f%EIUHZ7dC=>6a}RvQ{*Zn9UNAxuK8J^!=sfM1ssu2o|f%hK`R=cj6(j91?!y!
zc_ZK8WI|%nxmqiniU-a#()^z0kKZG>je)nrdl1+Dja;bxW*4Mwq9)cAWpVE}osiWn
zx?bo`UpO3kCP=xP_USD+H!IrPULOvH>;hsP9d==AW@J-5<3AY~0oq}&>cQa8e)|!}h
zuFeUEU)#WU-(ZbSjI=H6Js981ko^&R$u0H;A|Occh#I)78-Lt!deR(xD%x!Aqr-V20YZ9dL^_w)Qo%+0*d#uDZ|R+UG&OVD&dY6Jpqdg~mTY
zL?u%+`~3uM!_u-hFPF@j7GrvOSQwmg!7%ebs!T!hAabtmhD3~$i2FO;Ke84Ao(uR_
zbr1voeKOIF9VK;MN={?>V^|oYOg}^&gI;)z*PtlBdrMG&{{%mE3W)iviG-tiOirHp
z^D0t~Zmpi4N;8OGQ~{%}7BX}RiF!+Fc9SoOyeLPt
zvM5hnL3i&6aaD2(5OO1e9B`OZY>H|Ec(8K_qvYI@guKwO0C3%$CseJ9sQjX1BOr>U
z3gnQWuos0)=d*=X*_@u2C;ND_4l)eHr;(?x=m12h4Ge4>9wX^sqKbDQe&k6aC#H~vE(1c36MnDH~obR93Wj5mSJdj^Ce
z$&9;YUUSmF`>rO-@s46aAZfr`+vd^U#Hh4n(LMc%WNI&ud1E|kR)JZm#VljJ;<45h
z2|^Q}QPzc3?AeUzA|Aqt;(T=M83Y7Fpi#DfhZHXhUr|wM85Ix}7vw1E7ITrUphV_(
z-lH44fo|y?9O#Da9|f%Po~H5uj$BTVbHbeeU7f>nr*HiKeB?DY29%UW`8>nZ1>y{Q
zAbn|HR|Vh@@ZVWr3x;%6p_2mm&;he@EO4o@j@^ZGBt+j97v)>hD(Q6LKw{WqlA~Q&
zBlXBK-GK46HjzCbM?Z`Pinqm$x^8U41(koZEUuc!AkbJ{Uaw3doe!r?uMNT|;uSLRjJpYq=rKM3O$GF2%I0^j#edIF9khP
zyvC(Uv<{SLBa~SBx{3GZ1kp=N6Q@6`Iw7O%79-Z1BS`u437g~rg|WZu>ARJEGc1fv
zrg2l@pUDLjgL5i8(UxaD(JZY@&MJQKfGA^l8Y$7?8R?1V`dOz$knuC}h1jy)Aw~l{
zvYmZJQ{*)F5Mn!WNM5Gcq;dmLPdMa94x^W>9Eb3Ds{%30DS+70;d@z$I$}=Zfv?oZ>!3DgI)AH`K%C{{}ZT)-C@o
z`%A88J;gx>#TX|-qFjmN3Wva
zUNZ}2VCU?kpp&lm2GE1HXn%Ld^Y>Ihsr$ww;lJw*l~8YtKB-uP;z>q~=i$zv3&hdp
zpPr@adW;p$fW8PR_Ud)25$lw~
z3019T=7uGx)INWU`+Lb-{I~XGh#Fp5<$LX_DUl+Dy9)Zq!RSew(IOksYsA?$0m^kB%M5n>`8#`>Lo}vlXK0m_B>K0%NzBh08=nC
z1}ezZ8dSx1>v%B;Ye`!4=J@IhpMnj7OOu^tC`Sz$j!4
zX$p(QQX-QHkthDVlSl&2qrd0M+bOkf1knu3Pzw{;!@GXfmaDld)CdK)TA2J?+$>Dw)RbcCJkclkFG>gQrUk#xepF87BNmUeVG#s(v
zgH`FXP%K1nCzUr#uPw~)0z)`f7^7mak$}
zdbb{Ee!|H_h>42^Vg~VMlKsD4=C*tu4$tRKo&%AKxmrb>@$b$jwTSDlm@l$}g{4JA
zgA_?)&=&ZUKl^GA!TpCZIt8L@^UVMG7##882C*cP@Y8oPtH_$&9eI9=AHkRUJZ|4&%mvP@B^cbq
zY|GRaAtR-HqV14u<>b#{%gS-27Q+qs#q4|vRe_By?t_?BEWI44l65PR1j=wbzH4h5
zIR4v8v;9W5HAz>``{W>HkrX4}ztOnOpF$E}glF8xel?Q3%2-oFtJhzg_f^r%PPwzT
zJAE-|`;6J^5LSBW`YX^<5$e>=a!C0Wn0{cpJVwn7v7!xx!Jp(fbM_YvpX&T1Q2Vvu
z&9_@}1wDBeq)lM3hwN<=ZNJ*Vg)w5Ob@#jP?y#l_A3QE*fqoC1Z-C%w<{*WUlfBM7
zIZbPtXDUyaBO0wN#jgc{lIbs1l=nJThNydkwP&S*oAN9cjA+1OASq`7lXuasygOD$
z!QtY3jY>?dPXPLC+v=2b$^SEKrfd0q5(RBRS5u4RxS#w^ixR`V5n^nhlC{acumm>dpVD>|s5eg-S$@YN<^*`DZuZ4H6&nH}iSt*fOOT7&|
zeO^rM5LYsWOmo9(?qnsjAF`Z!`52Mb;Mn34Zv)jPj&8>lL|FWSY|?~qe06UbQnpQT!#_M|J1Qhr@|r$Zec8-7&5steyuMdT@P
zTT+%ud-e;ID}uz=?7Av|oeXY4`XKFVA){Tl`kH}UA-kZZFuLRO>8PLh(Oe+`gG=2Q
zskcj)bQ)yRSv(e|p1uKfrYwTONGIP#^54VC(`BETFJ5cQc>Ab&oV!l0&No-ecj#*O
zv(%RJW0Z!H5It%REGu~N#9!(T$=?Gs=|jOvU488Ku=*^fejMRK3dDl}C00&3rB#3J
z`XXTB>EXnu#%5#%+?8U1E}o@_=GswpUGff5P7^BWdaOz-=8x70e?XM4{)*MYE9S5g
zRSiIdD!RX((`$S;bKhR37CJvhY$%BgV67qR8TGBi(i82Nozj=tdrfuG?RbQ)*Vz5k
zE<}0lzm8i3wyewCA;AmOt(uIqblJVV(^SW`K=jY2t~@(di^QTNtXq?ZO&VD#U%Ni~
zF0+6sK&mq+!^&Ov7TCS^@3;w;-rD@Uwzc?tExU(Vr;G%;ZEN1$WW%jU%eV!pCi3+*$fL4**6+d
z6z}wJbS3S^42}F#pL$MeYD4Fop*Vz_&nDXT1P}X(xI0E2xk87(8XTR;q|2z=Gv#*92-rHBI<^SyeBUZO5t5|Ggc+`+SDuVU
zA^F!}ox5Q@Sl`~f$7a#)uD+TT+7q?{&JIdOv&a3fUZv=pRVv
zQX}R2Q2q@q=%6kWn3n#AeMQ5K+6Z?S&N&=xb}OsBw4Ej9-9;pxYkz~p`7*Oqqh#L3
z3#io`W~?3+x-G=|v4YSC4?&%+bL*jys@iJr@!rM|^m8e_Q=;e#I4)BMZi-eW-*nlp7x(%hXw$NY3F(8qgOuP{~xhSfv<}b75g0+tyX)eFfgXF+3rk&P03>
z#3bveB43K{*gl{|hG1;{xI>*Os=y)i4|}4aMD};rmLd{uk3W<)6Qup66T>;AvJz9>F2s>Kk_s)Ebk0H{gDx?cg`$HH9a-+?R?_#
zaDL&j@mn_jbg4lfFp`hw;H6avr&R5<(g=%I6akK!6wnsy`lRVo9PquV-?nA
ziS8zKiSFrozN6O$!t3(;r8?rjExp7{!HoZc@v?8(`k~tqf$b4C@X%b;VM?5se^<}alky8Yh*z1uyn>9y>3PF@v?Nk4BT`Q5{C+Q1kBCE+h=+LHV}=`@Y^8S)~mvf4nHU_xO(1QNWxC6XQAiSfM`s8uHOVL
z7VUw?!gG}`%D$tjEQ!iIVaiV?T>|N$sz^m9W+BO}rvzxn-e@s%S5=$s-X11OD0{Db
z0m~F7MMW1LUADSF?dbFEuYYJ6-+g-LztXEh!y5}Ka(CRYe+xICh`5PRc|q6BI{_z&
z;@e>;)TP_Vy#TXWdxdO#+xrGz_78Q0WzejDys7%tPOBQMXy~aa;v}6Jr=SAdlCB)B
z98KpmdFcsHqRcVd1C>Q=CbU=)*zM~!EGgBmDFNd#9EV!JxbF^io^B0)gfAPGl?}*z
z%6#96o%|#5o46o4WtkhBiR+$dkBc>|hgEAk?#Axv;eN?Ip&BN1{w?BXMd%@m`wXYE
z-VC$v9pPvP^fFcGx>QD20SSrYA!_uEc)oi&z7BIMOL@ak#%DldcS%mtD}-mI}dU(vh~5W-S_nNcO*w*1nH;D2U0hSlk`48-#Ws^ht0U&gkls(Wg+8
zwjtn(QT2(|q9r&e)G{QQpeqRiJ^$lp+wSS?pznxHq5eC07CPOlY-|6RA~}(v$A^c9
zNT*H?1l(IW2$1?wL1;qs_l`f_FFfvnLQjwgq$mTqUt!bpLL8asPPRh!$B!z?o|ywS
zBV`Lno6fV-7tSQ?LBHxE{NJ&yhX&`vyY6gXNEB6Q&~v
zp_xj!5*kiFP6RX@C4n1l^Dyl)4L~
zjdwSR>~bOw8*cPyK$qDrQ@Ky#=sA9;t6Z&IDgos->MT02;*6J~@{ok4!qry|?5BczVZQ@pWM$8`N3~S^zxVla3E4HV
zA)Uog1jRswT}nFsljL%6JTYDq!1L$O0p)AeQEe2>r%N6f_0pIvL(7O(&!bf~R5Ogm
z-{=vkxV&;&@Q1WW77UAQmUP6|zRgP)64W_}HDaz?OwrBJal2F+d;B??>I@htakj18
zIjE*_ckS5Fz597<$Rnf^gZbQ_tMRr{=fdvH_XsTw+7c66)j-_9`z4`v_SU=<(2d(h
zSdp76Ym4tZyv>-^`}I-L7fkjd$6ygO0^bx7@fmn>_3G8x;f{`rV4a~2@wMm2{SnQ)
z{QTIZ$w|YsjEod$1h#EFr&;;w_n7W%MoJ1U85vn~YwPouFB^kRc#Ld)%CBdhVt>i3
zQSd%LpKxesh`ZpTmbNz6ppe_@i#KlwL8_-_XUk&7{^l-tbPWzRO*I6VjYj6#IQ91S
zK8KFXwDk1$m_k;kRGDpa3EAJjf4_hK9#!AB=_n!b{L88tYJGS28-#hgCMK|ff%wGN
zu6-FCbStV7N`<|JE{0{mngd)lH8nARCpOqmojT=yxL^9kz?G#mF_CJvG9%**U_i&u
zjXd&+uMME%^wym_B%eNg`cPV$P#S-`qrJ{u(Eh=LGvjr>*t|TJ`K6^!&Dp`@=in~q
zmzQ7E)Equ?ad?=f%!nq!M;SUg5_KktHi5$;Zn5YcMX5X$z&6FXycDCDeb7JO#YU*C
zuI^t_!b3nv__XoS=dP~HZfoO>U{h$KWWqv2{lMGKFD~9PnVC
z`ThO-1Sz7dSCx#9kI%c(x{ASIIJmhbtgZ8O_X`hw*C%nPggxG!CnwI?vuDwCbe~&VWSHkn-<6m15na9dtM;V%X`&aG#m^~`3lG!4KVJ1P
zg#Lg2Eh(W}R4%XTdpu8w3C>XgIAh@KSe(LPS_n$s#EhHH%I-Gy`
zMYLAseZ`zrluLM$A1*SzH{s=
zzLw<4i{OQ)78hHdBpyA8OGtPF{`+%p?^QZFI!94%9-bGEn?jy8Owpjxc;MaS5>r!Y
z;YJ%%YuN;sQ}@KgaKLbc#l>B?81OtRivi`AnMv>H=y>_^WfUnXY4DpjjD0%?G7n5l
z=$<`$R{zVVw~WX?<8y3*3pGP$~Xk(QQrwAPzKU0wZhpd|eHP*(O0
zFje>XcnFy9nc3NtiVAmqn=2{VAViXy|nSLa9SQ7Vkk}d44I-2zQ;ZvG+6+d
zMl(LPwhm5CzB%gz*PncOOz0A#0;9ks1lOskzCrcqF{RhTJy5GyTwzeN1AG3e{3EYZwX5t&5A3D0a4t;!lJU#b0
z@L(wp9v&73VADMTc&ba$@)@hD=Kzh}QTke}?%$uzc7FLsUtV4o;9xQCyMO6Y(7M0t
z*6B_C+B{*Yi?&cp27_FGJFxt2{Q_Kzun9Wtw`@yGi-NtqJg_5n_{kF$p!0zv1_m}ZJPZ>~repw9RiCUi
z!4aZH#uuMaqcp=e_x*dznCnwf!VUqD8tXVxsKXB!XgH}Z2bk$pQ+!NnPPfJfTO-3vmw
zxVdL$ce*qIgE;bgnt;=lc3|G`2>SE)ucBDD=q3G_9S3q6VsyTJ`^NSDVSuQwl)B=q
t4!K-lsjcXex#i|05W{P+LlX-9o4-EkSN_Y64XYV|z|+;wWt~$(695}}@-zSd
literal 0
HcmV?d00001
diff --git a/imgs/plot_polar.png b/imgs/plot_polar.png
new file mode 100644
index 0000000000000000000000000000000000000000..97e884c8b8502a2986804a4cee386e20b6ada92b
GIT binary patch
literal 123142
zcmb@ucU;eV{61W1skG42j*v7+OM4RuZS6r@+8Uxwk~kzG$!@2lsWd4nm6o2R6kss3ZcCzzzIqfNE4Yr?|o+jj~oA0oelW-*yE$Ge#$hXZ<=qq8mBX*Ivb~~Zo9D9
zo=mYzj5eDuB(oW5IjOPgs@;ln+Cm5WEVY+
zB>Bp^C2dffe6<(ZR7SpPxJQO_lCQDxhO@%|{x%zfHTfzet5Fq4zSijcUtWZQ;tH=q
z`0T7#a8heV=*m#ct-YtNls5`ScCq0p?wmK=Av5sGQeq{f;^)AA|L?g0XFG?yQ*&s1JVy1;#E~NR-!Ga~H0|xl54803T)Mb%`APT@rQbe(vr
zYN>!MNm*GA5uM9wc6Mn$E0lwmY-b0nW)~N={rzPntjnDLntOQcetamID$k55EG(?-
z$Gcq=6w!9GRdfP^Vqy#zE?fxw`*S?aa0kVzO>@VmtP^b9+!xQEr^qzRd-M6EwKf08
zwzi?EsrdW%1;s62u~K*qgf3H3goK6)ZDKg{Yr)06JD&!(uiSFFj51L#(y8&LU3b1y
z;B@hL(^OyAasQ0Y<5?ynt-1o*nwk;Cr}ER%h91B4pXrZ(_;BmWuV3~}x1?Wrk6!wD
ziVW_++TQ>s9vONS6_s0xeqz3geO5L;8nF)PgW9pz-CSM64~c5}`R%J%{h=?#6p_65
z)Ps@|o_+iF4bJ}@qgY4L_T)&6PI}cN{iK$UAH(JP14kaq4uq|Rl4sOS;MWcelsCz-
zj>nqO+sZB@CpS{fqDmc8UF_zuZQC~abF1?94i1zQG&D4~_PKG$dh}2~GrGiS9Y}e7
z`;lntoS3(cZ)z@XWZS$s%B?dyt)hbei0U613fsn;ySwt87?}lBE~ke=R(hsQhgFDmqZTU=OJsF8fV8s5@k_|ynj7mJLJwqAafpTFl&itGkF
zNR*av)ymw6w9Ch>zP`SEGQ__Yx
zARu7&XUNL$-y;{d-#f1+cs?nK^QMH=rF7i{Dhi%0Tjm-iU)zq<#m~;quS=+}#=3s*
zktZG(#mG5ad10M=$cj{ry)?!7`jND}JhqOGj={M=u6?fUd=Z^3X_`Y56Bi>RsV{Ty
zrJWydrlSa1>@W=+isF=V?~;1;>XmKXWp1y4fO<^d)ug1D=xEwZi$aEb_wLQk&Q>-g
zicnHuyH>u5Vr1jv+mL1$dt;aBaOq4%spq**_`;}?n$X^E^ziZH2Tz|iSO#lr(^#B1
zfyF2tvNF5B;$B)>npv(brKOcs)%g6(N@+l_@{Hs7hr73=?5OJEx6tt?#K;YVY_!PY
zcw%sO&)SEKoQg{r+gL5(L$_rnO73ZD#QHoS%WY#^Egnc}?L>O(U=*jdnvv0_vhwmX
zSpMn~N!5;f(q6vgz&6DGEppeFNE*S2)eP6h4ULRARIFtc^=r7hi}MG471Gkx9Z646
zDRDTAg=P56sH@m#`-gkVjD;2(D4u3#*L~~A7@C}PtLXUj=|Nf=Wnp1q?4?U3qj#QT
zV6?Qglh^n*Z>C9ajU5|1IUib8RW&w4Lv`V-^=H+!pGW?=al>w~@)Yj6XcHG1L#rx;(-XHVx)wrOP&)A_mN@qWYyvFh%I|u9g%h;Ix``i1-!k)1(
zH85b#&(EL53a4CAZWp&I+2mdPIi$^5=9a7nkFc;ql`>jW3uqH8qb~S$X^oIDB{m
zW`=FkCg0A^?Gh5qYk!wtljW`*$1(JFdupUbR)t4GU;J`#P>@D?>#eoE_IAg=e+_cR
z`=aenhZu-AKESBk+1Ukz{?FT7wSHV)&KII0P*G8F^{t~>MM&|_Q$j*QfeWp=(qVr?
z3JMFWY^82Rm_LzrsI#rV!WXzYSAS;dyE(<{Be6M8pHejD&n_)R$Hg&PT3VJEHSeD(
zi@CyL}9p}S0zxpcrtfb&=
zRsJj-w=56ZEOvuI|E73V`?KSvzI8)OQ@uiodX=Nzm}5o?uQO*Zn2FOSHEJ6hckl_Q
z2F}`&7=Q?7xMP!9vA0u_bu$93-lJpl2Pov;y?tBNUmiR+ne+OMt(5KJw+?I@U4MW7
zBx}?ER;ixk{IAIz&5RCVK|zY@7?ygG1SRjj<|2*$1~q*xIn%Adub8L6|I
zAA5gG>FgfN?on&&E68H38fmQN=H|Cl!&I;-+a3iSJ9bPeXu+9m$cc%Gm2Vk2v+s80
zR5mxWAj6HVACZudXn6l#8=Fp8M1+l>|H|F{0aO$x-`|d!XnQhJ@q1cx?eCwVp`iE-ExD>RH8f_UfR;$l2cR7%Q;
zqRO^M8v|n4S7c{PYYAbnxDvvaJ{?uqwIaJb8zQ
zhow9}%VRwzhy5EKf`SyU@GHCKTw5>nu#sRNUY#>%q$_@X+x2C1G!2*{cQ(>RHGjk<^#8p?9fGA)$ch&|HPx$<
z68nX(tt%5xQfKFzRUI(F6z943o*YErV43cDY%It@%C*7)6V@fk3Z$uoKF8nFZ8AVD9}(%TGKt{EeMD5
z{@=N@gxN4QteZAP*461ys5v-@P_0`xWw}i#a)AjzE>?4bS10q?v#JEuzp*ZDkL`YZ
zdD~Uuw~NAd@>5pe_q_V==EC?dFdWGxK_lUpTR${}H=}UeVG2g`zyJTmS#Tr2j4$(_l_DdLgE2wMWO)@yRklvbbP(
zE5%d1F&snQIQpixwsuuHsgWkNImU2DBx0+jl}o1%{h31r0m6}H2!c8e4(UHXwX{$q
zABtXDTAK8{J9Np(JXf5NIc6O`l9H12Q5u0kLE-fy$uj|1Heom3Is)|3R7*)E9~zUr
ztRqrjvD-*&^FIf-9TIKxEBqtDIOia>;6f<#Iy1T09Hq%oN
zdY7uU4IyiD?LB#@h`%~t+*gs$NUXzeCXnICQ
zO3Q_%
zQJMAw^u@bzn};l|lZ!hfBy_NU-ZeHJJ$5V{(NR68Q$$8)(>huvHA_q0;gJz(+e(Va
zsHnxKr`}lveH?MxE#9ahcH;pa4v=n#q@=F7Id`^IDQ&bJ*3{WPw`44H0Re&Fq{H-i
zdr{Xc98>Z2l0DwXvl8n*v!Qj)H%J<$%l6A+Ex8?e*DbP8=6g<
zHmRiiJt>q70jLp(2=p5UenE;-KX{Pg8uZ)6RDHj!jCZ^_RTpn2VfhM+;!0PS8HohT17WHN9(=A1Nvt1zQ1+dnLRym#
zTG&B=F(aphptv{_f>0Hyk#67qCZid({9KMU=fnh)Lw$lk{ZMtxT<_`ELa1iQPu-Gs
z*b&P%Tth(Xc~2$_teP>FJ|WHMh3*$)zYflpur=nBiYrk>*Xv6|;$(Q&pXGrB1SSUi
zNcr`c19=L7_jPp4gBZkDUPDyj5fO8h>sS_k-uFwd;3@NL=dK|SCZc=g{;OA9$ZoIJ
z{{EIew;EVeQ`1$drP1JlO0~w|puT=$N=h8+s!C^>VU#hoBvbi*HDs_TQm>l)1l)ue
zeeTbqDM1GtHf#t)B}%Y3(-v7KWXmOqMAjd#@fP+XI_YeHoEn~<63l#xkyzCLoa;Ak
zH_$RPyjtGq!?3BWL~a|R^XuF?QW
z#1YX$x6UK}VuC*PKVtj1oeU2ZqKvnX&)ihez=Psq?jy>Jo9>1DIH0AY(^Fbt9TlGB
z^mf$yd-n^vfpfoiZ(^W%d^jcg+3_d-lece1UY!O>pvRJYNR(H>o4U8?bUmsIsWZb`
zO25Cq(9_e)4R9&vppl$(@${_waCd*6&$t1p7k@2G+>&u&S5j7dcXTTCJx!`2J
zY4c{g-%C^SfwO}7FJ8=j_GuQ{zkdrabAA8lV$h~C_jGJtloe_ieXLUQ_d^%x`0u~)
zSf59M(ePC#W#ylx6RRsD310_ZtKJJ;GFSTbH9hoXZJ@4tRpU#uo3-69PUSh&GxI73
zZ9-zHMPVLSP#~k&BKR-{6YRZ(HiD82SAok$#NT_t(+Y=@_s}5VE=gpW+L?2RUSE$7
zQ#C?eu6S_c3?l=d2(*un
zq#Df7X_Q>V?8`kxE;O|XjYSy8hf7-^Xip+aMI{3LW!;}-9pfx!$nuFFmo}{
zg3`c8ipi(Xfn0xh*q+(2h`|{5TwGkF->{+FL2YmHAt|2t&cfy)FWjHO)YKHY^Hp(i
zRbOA<>SBVNh5N$v!(%&NdiGtK>Me3W5DuKHFTpJm*+8&@uB?a@vGY
z-Z&Tq?!*4wKRbc*FoDdNnXijM01M8A7WnV{G6G}e%+Q#L`s2NWxRKAwgv+am2}6-;
zsj-$*GCZfZur3CbD{T4^s00Xu531DpWqyAC0`Z15dvgxmE4)=$L?|mF%wL=wN2%iY@*9%r_1?{q5n-{M^@f
znn4If11>B(Ig#7JD|yq&PKqm=@z&jeNnb@ITK)vOG(&oPkZS|ap6y;J^pbi1{{70h
z>e@xCQvV^)CpD`-=TV}1%;eLbk>6CNmuh5bdAV72jrZ)?vt-u(u8ylNKPwjFq+w24A$Ppa
z>gs&+%5a>dAS->YZ55T$%3wshzFxx`s?=XhH=WL`%-V&Zs6h!HIC{f$wndX=@fQpA
zT?4u(4NBx){1>Gr-Uk_KWo6}-CPza3i-2bVqewpfheS8FyYB1h?M*Wh0|#`!!*ytS
zni-g2esX$ZqI$5Q}i{{_w(NXLG-k
z`9B6v_Qqs_LX%=zr>d%o(j6ZxlqVa^i-w!3CQUHgO>1GeKSD0{AL~V}+2!ScZ&yQh
zfqdB_>t1<-
z1SxN2Uyl!>*yqomn}iyhn{zr@MiyQU`iVr_fy8;Qeh5&nt!*PnJ^Q>@oz&8I4Pthm
zcsWRj0VIvtx-F*#<>kBjo!IO*V!7mCyk~1UZ3aPVR2R6k&3@FERQMBlHocWT`4H_z
zt?*BuOn~7@Wrr)Gye&U>t|vYy_}|oLn9*@nZ>zJ5%LSeE3m8#AW;Rw<3eY;iNr6WY
ze*C9DU-X;&B&4j&hd?5Fcm07NWd`e-iNb;c+o_(ymC+;{seozAraKB7LsSO(ofhdE
z^PNwh)&<_fOX31#jyU3mW@ZM)b8d8YhIhs*2G=MCD<2l!2xtYkPMxUtr9l0TlQ~Iz
zsUbhsqukF8oU`A$eS2@IxJG7pIOiD)jCY1#7FHVARz(6P`H#hf6qs`xNJmbIj
zLgF(3q_+c^VD{qVhQ|hXFJYm)e`at%j*H3J-CZj*R2eIfR5XKww%BqDm)5tvRIAFE
zl;PQc2We?4bc&U`+48MO*4W^=5l&Nch|9+egCwj)x77jyfT%6mQ>D7F+FL
z&Vf>}>i5j+gw}Ljhq>WeEa)YP41LVf?UaFXdgPn
z1d{N|zuda7DL~(&M>I$nWPKPKvc7rqW^iiCW9H)>ewv)%Cv6jkHZ}<$t3d4$j0f_?
zATKb(BCD>%Z!*`dQ)0DwZI$KT*)eXO+uIFduqyOW|Cp8d?gSU9vO24|b?ery(ubwT
zc5DJ27HuRJWhSlzejXIF3s-CSjzII9I_c6eQWw;zE^a3KAFnlNvrV5j&M9LXyB!^iME-U{2a&zy>
zr6e5B@gvHkxkFJ-=Ae_wb+SxEl$3Zuf5?p2(Bi>U-#Tibo<;46rAr7CcJ|>Gw}#R2
z@e$Bxx+hPH_J=N8@qF6O{Qk_qfrrfvjg8eHd#(V~4uffC-aluV)0qy68-tHIkG*&?
zElum=BnJ8xizpTo3M4WFpmEbtR4|`4u|o+!uzc8`NgFN&@2!NR!lsMwPpJ_?;hmTA
z!LAnjaTw6Wr|YaGV4tzEF_~77n{TSC`F?$UEtk$O`Py=?fkbmt$F`K!g)ipKOk?KfJP7d5>^}xWu
z(6{eyI|93Hf`j*im1Y1ikBy}VNL2j7ohY-a$9u!ax+_jTtE^1vGS5B+Jd*s`A6vb>
zy?Mo&fFAz&=bu3o4M?C|ioQGc?A?n#{X~4q=
zAARh9Wc9>nMH*BD7$*S&>kL)bouO_Y(RL@H1+2IZh)8H1gd*!HeK-jo&&bH=K9bg-
z-^;p?39?weoPS#1IOluM{g}BqxjM=mjB{^k=46B8p4CywFm)7xAb9jS2SxostD5y8
zE+w$2tW;tUHZ#zmyi9AGVBNTJBzNAVVP;35f-H)qOP3fx7}nI)3Gdj^kk)?$97rJ8
z@wveW7W=o?wlOPc9Sl;%pt@RESj_%tlqC5i_1v$WtE;Qut=P{GB(yv_wiB1G8uDW0
zu(3|YR;d3bk=9cRVhk{M9-+Z^+
za+gVFC05sz74!Llb*2c8s6Bx#a6R3Po_Fz6;}5W>jg6zgIHT|u5)*r87dx4Yj($y6USp4wOqoZGiYu+g-DP6sF
zD;7%{944WxUinY)a7E_@xab`^RD*99j6r4XDcH&nIehtJZ4QA=vA7K-1W$Q2ZR%feOfWue_^Gz73B?&rgCev)o3mF-^
zz#zED=Ov}_=k#l1|2}v?1(-wF1wFmkzUJ^U-wAdC+G7OHabYPw?w7Pb|&hs3w-$!Gi|{
z*)M{ap(vGDynpPO6exPHV#^Rnkif(E?FHtA1KK58VmF5K>iLP*ha%m*v!w$m8a0>t
zlQuh4`(HbA1TbMiiAc#LV)Xa-lZYXL|kvbAz#4ueOpl6tzxA~w`8!F0d@WRkw&R#8z6agN@sn;3l}
zSV3H79{X@tzl}#P0RW(=R*FgG?}`Jjycfkl$ovMp=hQI>gc=|NkdqVt`6m(@YH7Rr
z_k?iMfk!48M=dRxI5|1DZrj#9RjHrRF^mKdqmfpbWOMGoDTC{t>2uuP`iQw)s^^rz
z+kqi)Ec*AYN9q$IurAx59^F(zKaq~QR~j;!mX_8a>S{0`xs9`b~Gn7H^C3_L1hhx2}?mV+zgYfH<~QUj!8U<37%F~R(M)!78VvfEidC@
zAlu-w^(^8PHdc2bXBk#6Y8oYia!_^Nh~__#g!ezmRE^cIQaLwLKlfg40L*ZW#6A*2
z^YRj-H@`i+BQR0**9fJ>!Go8qKAdZxTnzj@Rdn>kiMXs2MQ<7#8^`vsj%#SfQ&-OH
z&%JYpx1gxV0cpm6y%l!Rc7r<<6g(zdBR?CD&pLp3Hn^vB5ilmj=3GNxA0Md~K0iW6>{Pz6^qWuhW+-
z*ZM8vy%Bq$i&X9y1A8cvU_oC@yxc&$F59x03Y9Zd%~BhlNO?L`vb%uoBy1|&i~ht}
z=hVkZR?s7xX4<^&&)zk&v0igTA@nMMTsRG**TUDMyh?#BV`=)>Ghg0bpPidiH!`+($+(+0Yxm7YvS!RJFV`d3qjUkoM0eNV)C&(EdETFL>k2#d14`{jC}&?E%RW5UN?k3y!7rOjqC`IFPXJce%pPql*6XXWCG{qn_eK9pb_Y3CL}
zCVnL@LSBKUrl5fKy$(A7w3jrrwK9~=WA%hyf|x~CV6bH)ku^hq`MM8O@L_Jh|GfT9
z^rq8MWIIAp?%TH+OkjJXHho?(bU4RTr_vukrj3e<0!1R*yL6s%o!8eEq1E}+wVL3W
z@>6+sN3Oq(#M-*GJOP7?d|Oi&CFuoywPy~{Hir%4eQMd=aJH|7Vq{TO4@wUnc?`4MnMpvt5Ik@qUQei9ls*B(cDodo`PM9GwyRK#eI
ze^pXVhjrzO9QGs=1eFeh!v28ig6=|Q?r35SA;<>gP;9MzihB~PMhSsf-QC?%Tuj7%
zu@F>#_+{ge*X*Ygr|s}qL`M$(Jh=p!b>#XHm138?2}1Ei&~A8x#n|%c6D