From 149a779b5c15796db629adf2959376807e7a0b32 Mon Sep 17 00:00:00 2001 From: HungryAmoeba Date: Thu, 1 Aug 2024 23:03:33 -0400 Subject: [PATCH 1/4] dymag solver version --- DYMAG_solver_version/DYMAG.py | 68 +++++++++++++++++++++ DYMAG_solver_version/PDE_layer.py | 95 +++++++++++++++++++++++++++++ DYMAG_solver_version/__init__.py | 0 DYMAG_solver_version/aggregators.py | 66 ++++++++++++++++++++ DYMAG_solver_version/data_module.py | 75 +++++++++++++++++++++++ DYMAG_solver_version/main.py | 10 +++ 6 files changed, 314 insertions(+) create mode 100644 DYMAG_solver_version/DYMAG.py create mode 100644 DYMAG_solver_version/PDE_layer.py create mode 100644 DYMAG_solver_version/__init__.py create mode 100644 DYMAG_solver_version/aggregators.py create mode 100644 DYMAG_solver_version/data_module.py create mode 100644 DYMAG_solver_version/main.py diff --git a/DYMAG_solver_version/DYMAG.py b/DYMAG_solver_version/DYMAG.py new file mode 100644 index 0000000..eaf2871 --- /dev/null +++ b/DYMAG_solver_version/DYMAG.py @@ -0,0 +1,68 @@ +# make a torch class that applys a PDE_layer, then a KHopSumAggregator, then a GraphMomentAggregator, then flattens the output +# before passing it to a classifier +import os +import sys +import torch +from .aggregators import KHopSumAggregator, GraphMomentAggregator +from .PDE_layer import PDE_layer, heat_derivative_func + + +class DYMAG(torch.nn.Module): + def __init__(self, num_nodes, num_features, num_outputs = 26, num_layers=2, num_lin_layers_after_pde=2, device='cpu'): + super(DYMAG, self).__init__() + self.num_nodes = num_nodes + self.num_features = num_features + self.num_outputs = num_outputs + self.num_layers = num_layers + self.num_lin_layers_after_pde = num_lin_layers_after_pde + self.device = device + + self.pde_layer = PDE_layer(num_nodes, heat_derivative_func) + self.k_hop_sum_aggregator = KHopSumAggregator() + self.graph_moment_aggregator = GraphMomentAggregator() + + self.lin_layers = torch.nn.ModuleList() + + input_size = 26 * 3 * 4 * 4 * 6 # TODO: FIX PARAMETER SPECIFICATION AND INPUT + print('input size is ' , input_size) + for i in range(num_lin_layers_after_pde): + self.lin_layers.append(torch.nn.Linear(input_size, input_size)) + + self.classifier = torch.nn.Linear(input_size, num_outputs) + + def forward(self, x, edge_index): + print('input data has shape ' , x.shape) + x = self.pde_layer(x, edge_index) + print('pde output has shape ' , x.shape) + x = self.k_hop_sum_aggregator(x, edge_index) + print('k hop sum aggregator has shape ' , x.shape) + x = self.graph_moment_aggregator(x, edge_index) + print('graph moment sum aggregator has shape ' , x.shape) + + x = x.view(-1) + print('flattened output has shape ' , x.shape) + + for lin_layer in self.lin_layers: + x = lin_layer(x) + x = torch.relu(x) + + x = self.classifier(x) + return x + + def reset_parameters(self): + for layer in self.children(): + if hasattr(layer, 'reset_parameters'): + layer.reset_parameters() + +if __name__ == '__main__': + # test the model + num_nodes = 10 + num_features = 6 + x = torch.randn(num_nodes, num_features) + edge_index = torch.tensor([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9], + [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]], dtype=torch.long) + model = DYMAG(num_nodes, num_features, 26, 3, 2, 'cpu') + print(model(x, edge_index).shape) + print(model) + model.reset_parameters() + print(model) \ No newline at end of file diff --git a/DYMAG_solver_version/PDE_layer.py b/DYMAG_solver_version/PDE_layer.py new file mode 100644 index 0000000..7158377 --- /dev/null +++ b/DYMAG_solver_version/PDE_layer.py @@ -0,0 +1,95 @@ +import torch +from torch_geometric.nn import MessagePassing +from torch_geometric.utils import add_self_loops, degree + +import torch.nn as nn + +def heat_derivative_func(lx): + return -lx + +class PDE_layer(MessagePassing): + def __init__(self, num_nodes, derivative_func): + super(PDE_layer, self).__init__(aggr='add') + self.num_nodes = num_nodes + self.step_size = .01 + self.solver = 'rk4' + self.sampling_interval = .2 + self.final_t = 5 + self.derivative_func = derivative_func + + self.lin = nn.Linear(num_nodes, num_nodes) + + def get_laplacian(self, x, edge_index): + edge_index, _ = add_self_loops(edge_index, num_nodes=self.num_nodes) + row, col = edge_index + + deg = degree(row, dtype=torch.float) + deg_inv_sqrt = deg.pow(-0.5) + norm = deg_inv_sqrt[row] * deg_inv_sqrt[col] + + return self.propagate(edge_index, x=x, norm=norm) + + def forward(self, x, edge_index): + + if self.solver == 'euler': + num_steps = int(self.final_t / self.step_size) + sampling_interval_steps = self.sampling_interval // self.step_size + num_outputs = (num_steps // sampling_interval_steps) + 1 + + outputs = torch.zeros((int(num_outputs), *x.shape), device = x.device, requires_grad=False) + outputs[0] = x + + output_idx = 1 + for t_step in range(1, num_steps + 1): + lx = self.get_laplacian(x, edge_index) + dt = self.derivative_func(lx) + x = x + self.step_size * dt + if t_step % sampling_interval_steps == 0: + outputs[output_idx] = x + output_idx += 1 + + return outputs + + elif self.solver == 'rk4': + num_steps = int(self.final_t / self.step_size) + sampling_interval_steps = self.sampling_interval // self.step_size + num_outputs = (num_steps // sampling_interval_steps) + 1 + + outputs = torch.zeros((int(num_outputs), *x.shape), device=x.device, requires_grad=False) + outputs[0] = x + + output_idx = 1 + for t_step in range(1, num_steps + 1): + # Compute an RK4 step + # ref: https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods + k1 = self.step_size * self.derivative_func(self.get_laplacian(x, edge_index)) + k2 = self.step_size * self.derivative_func(self.get_laplacian(x + 0.5 * k1, edge_index)) + k3 = self.step_size * self.derivative_func(self.get_laplacian(x + 0.5 * k2, edge_index)) + k4 = self.step_size * self.derivative_func(self.get_laplacian(x + k3, edge_index)) + + x = x + (1/6) * (k1 + 2 * k2 + 2 * k3 + k4) + + if t_step % sampling_interval_steps == 0: + outputs[output_idx] = x + output_idx += 1 + + return outputs + + def message(self, x_j, norm): + return norm.view(-1, 1) * x_j + + def update(self, aggr_out): + return aggr_out + + +if __name__ == '__main__': + num_nodes = 10 + x = torch.randn(num_nodes,5) + edge_index = torch.tensor([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9], + [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]], dtype=torch.long) + + model = PDE_layer(num_nodes, heat_derivative_func) + output = model(x, edge_index) + import pdb; pdb.set_trace() + print(output) + \ No newline at end of file diff --git a/DYMAG_solver_version/__init__.py b/DYMAG_solver_version/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/DYMAG_solver_version/aggregators.py b/DYMAG_solver_version/aggregators.py new file mode 100644 index 0000000..b1725c9 --- /dev/null +++ b/DYMAG_solver_version/aggregators.py @@ -0,0 +1,66 @@ +# the output of PDE_layer is a tensor of shape (num_outputs, num_nodes, num_features) +# make an aggregator that sums the values in the 1, 2, ..., K hop neighborhoods of each node +# the result should be a tensor of shape (num_outputs, num_nodes, K, num_features) + +import torch +from torch_geometric.utils import k_hop_subgraph +#from torch_geometric.nn import GCNConv + +class KHopSumAggregator(torch.nn.Module): + def __init__(self, K=3, M=4): # M and K referenced from appendix B.2 of the paper + super(KHopSumAggregator, self).__init__() + self.K = K # max number of hops + self.M = M # largest moment for aggregation + + def forward(self, x, edge_index): + # x has shape [num_outputs (t), num_nodes, num_features] + # want output to have shape [num_outputs, num_nodes, K, M, num_features] + + # Compute the moments for each k-hop sum for each node + k_hop_sums = torch.zeros(x.size(0), x.size(1), self.K, self.M, x.size(2)) + for node_idx in range(x.size(1)): + # get the k-hop subgraph for each node + for k in range(1, self.K+1): + subset, _, _, _ = k_hop_subgraph(node_idx, k, edge_index, relabel_nodes=True) + for m in range(1, self.M+1): + #print(k,m) + k_hop_sums[:, node_idx, k-1, m-1,:] = self._get_moment_sum_aggr(x[:,subset,:], m) + return k_hop_sums + + def _get_moment_sum_aggr(self, x, m): + # sum absolute value of x to the mth power + # x has shape [num_outputs, |N(node)|, num_features] + return x.abs().pow(m).sum(dim=1) + +class GraphMomentAggregator(torch.nn.Module): + def __init__(self, S=4): + super(GraphMomentAggregator, self).__init__() + self.S = S + + def forward(self, x, edge_index): + # x has shape [num_outputs, num_nodes, K, M, num_features] + # want output to have shape [num_outputs, S, K, M, num_features] + + # Compute the moments for each node + graph_moments = torch.zeros(x.size(0), self.S, x.size(2), x.size(3), x.size(4)) + for s in range(1, self.S+1): + graph_moments[:, s-1, :, :, :] = self._get_graph_moment(x, s) + return graph_moments + + def _get_graph_moment(self, x, s): + # x has shape [num_outputs, num_nodes, K, M, num_features] + # want output to have shape [num_outputs, K, M, num_features] + return x.abs().pow(s).sum(dim=1) + +if __name__ == '__main__': + # test the aggregator + x = torch.randn(5, 10, 3) + edge_index = torch.tensor([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9], + [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]], dtype=torch.long) + + aggregator = KHopSumAggregator() + output = aggregator(x, edge_index) + print(output.size()) + graph_aggregator = GraphMomentAggregator() + graph_output = graph_aggregator(output, edge_index) + print(graph_output.size()) diff --git a/DYMAG_solver_version/data_module.py b/DYMAG_solver_version/data_module.py new file mode 100644 index 0000000..8745654 --- /dev/null +++ b/DYMAG_solver_version/data_module.py @@ -0,0 +1,75 @@ +import torch +from torch_geometric.data import Data, Dataset +import networkx as nx +import numpy as np + +class RandomDataset(Dataset): + def __init__(self, random_graph_model, num_graphs, graph_params_bounds, node_features=None): + super(RandomDataset, self).__init__() + + self.random_graph_model = random_graph_model + self.num_graphs = num_graphs + self.graph_params_bounds = graph_params_bounds + + self.graphs, self.generating_parameters = self.generate_graphs() + self.node_features = node_features # options are degree or none + + def generate_graphs(self): + graphs = [] + generating_parameters = [] + for _ in range(self.num_graphs): + graph_params = self.sample_graph_params() + if self.random_graph_model == 'er': + graph = nx.erdos_renyi_graph(n=int(graph_params['n']), p=graph_params['p']) + elif self.random_graph_model == 'sbm': + sizes = [int(graph_params['n'] / graph_params['k'])] * int(graph_params['k']) + p_matrix = np.full((len(sizes), len(sizes)), graph_params['p_out']) + np.fill_diagonal(p_matrix, graph_params['p_in']) + graph = nx.stochastic_block_model(sizes, p_matrix) + else: + raise ValueError("Invalid random graph model") + graphs.append(graph) + generating_parameters.append(graph_params) + return graphs, generating_parameters + + def sample_graph_params(self): + graph_params = {} + for param, bounds in self.graph_params_bounds.items(): + lower_bound, upper_bound = bounds + if param == 'n' or param == 'k': + value = np.random.randint(lower_bound, upper_bound + 1) + else: + value = np.random.uniform(lower_bound, upper_bound) + graph_params[param] = value + return graph_params + + def len(self): + return self.num_graphs + + def get(self, idx): + graph = self.graphs[idx] + params = self.generating_parameters[idx] + if self.node_features == 'degree': + x = np.array([graph.degree[node] for node in graph.nodes()]).reshape(-1, 1) + else: + x = np.eye(graph.number_of_nodes()) + y = np.array([params[param] for param in self.graph_params_bounds.keys() if param in params]).reshape(-1, 1) + # store the graph structure so a GCN can use it + data = Data(x=torch.tensor(x, dtype=torch.float), y=torch.tensor(y, dtype=torch.float)) + data.edge_index = torch.tensor(list(graph.edges)).t().contiguous() + return data + #return Data(x=torch.tensor(x, dtype=torch.float), y=torch.tensor(y, dtype=torch.float)) + +if __name__ == '__main__': + random_graph_model = 'er' # 'er' or 'sbm' + num_graphs = 10 + graph_params_bounds = { + 'n': (50, 50), + 'p': (0.1, 0.5), + 'k': (1, 5), + 'p_in': (0.5, 0.9), + 'p_out': (0.05, 0.2) + } + + dataset = RandomDataset(random_graph_model, num_graphs, graph_params_bounds) + print(dataset[0]) diff --git a/DYMAG_solver_version/main.py b/DYMAG_solver_version/main.py new file mode 100644 index 0000000..55cde7b --- /dev/null +++ b/DYMAG_solver_version/main.py @@ -0,0 +1,10 @@ +from .data_module import RandomDataset +from torch_geometric.data import Data +import numpy as np +import torch +import networkx as nx +from .PDE_layer import DYMAG, heat_derivative_func + +# run DYMAG on a RandomDataset of Erdos-Renyi graphs + +# TODO: set up training loop on this dataset From 404c6dce8482e63ade40cb0efa09542feb85330d Mon Sep 17 00:00:00 2001 From: Charles Xu Date: Mon, 5 Aug 2024 10:10:41 -0400 Subject: [PATCH 2/4] sprott dynamics, characterization of states and chaos --- .gitignore | 8 +- DYMAG_solver_version/DYMAG.py | 91 +++++--- DYMAG_solver_version/DYMAG_pl.py | 100 +++++++++ DYMAG_solver_version/PDE_layer.py | 169 +++++++++++---- DYMAG_solver_version/__init__.py | 2 + DYMAG_solver_version/aggregators.py | 70 ++++-- DYMAG_solver_version/data_module.py | 88 +++++++- DYMAG_solver_version/lyapunov.py | 101 +++++++++ DYMAG_solver_version/main.py | 79 ++++++- DYMAG_solver_version/outputs_plot_heat.png | Bin 0 -> 40728 bytes DYMAG_solver_version/phate_plot_heat.png | Bin 0 -> 39867 bytes DYMAG_solver_version/see_edge_index.py | 27 +++ DYMAG_solver_version/test_PDE_layer.py | 66 ++++++ DYMAG_solver_version/test_data_module.py | 50 +++++ DYMAG_solver_version/train_pl.py | 47 ++++ .../visualize_trajectories.py | 201 ++++++++++++++++++ 16 files changed, 989 insertions(+), 110 deletions(-) create mode 100644 DYMAG_solver_version/DYMAG_pl.py create mode 100644 DYMAG_solver_version/lyapunov.py create mode 100644 DYMAG_solver_version/outputs_plot_heat.png create mode 100644 DYMAG_solver_version/phate_plot_heat.png create mode 100644 DYMAG_solver_version/see_edge_index.py create mode 100644 DYMAG_solver_version/test_PDE_layer.py create mode 100644 DYMAG_solver_version/test_data_module.py create mode 100644 DYMAG_solver_version/train_pl.py create mode 100644 DYMAG_solver_version/visualize_trajectories.py diff --git a/.gitignore b/.gitignore index d662796..d1fa0e6 100644 --- a/.gitignore +++ b/.gitignore @@ -159,4 +159,10 @@ cython_debug/ # option (not recommended) you can uncomment the following to ignore the entire idea folder. #.idea/ GDeNet_old/ -old \ No newline at end of file +old +DYMAG_solver_version/lyapunov_results +DYMAG_solver_version/wandb +DYMAG_solver_version/figs +DYMAG_solver_version/lightning_logs +DYMAG_solver_version/data +*~ diff --git a/DYMAG_solver_version/DYMAG.py b/DYMAG_solver_version/DYMAG.py index eaf2871..dd5493c 100644 --- a/DYMAG_solver_version/DYMAG.py +++ b/DYMAG_solver_version/DYMAG.py @@ -1,68 +1,89 @@ -# make a torch class that applys a PDE_layer, then a KHopSumAggregator, then a GraphMomentAggregator, then flattens the output -# before passing it to a classifier import os import sys import torch -from .aggregators import KHopSumAggregator, GraphMomentAggregator -from .PDE_layer import PDE_layer, heat_derivative_func +from aggregators import KHopSumAggregator, GraphMomentAggregator +from PDE_layer import PDE_layer + +# make a torch class that applies a PDE_layer, then a KHopSumAggregator, then a GraphMomentAggregator, then flattens the output +# before passing it to a classifier class DYMAG(torch.nn.Module): - def __init__(self, num_nodes, num_features, num_outputs = 26, num_layers=2, num_lin_layers_after_pde=2, device='cpu'): + def __init__(self, + input_feature_dim, + output_dim, + dynamics='sprott', + n_largest_graph=10, + K = 3, + M = 4, + S = 4, + num_layers=2, + num_lin_layers_after_pde=2, + device='cpu', + ): super(DYMAG, self).__init__() - self.num_nodes = num_nodes - self.num_features = num_features - self.num_outputs = num_outputs + + self.input_feature_dim = input_feature_dim + self.output_dim = output_dim + self.K = K + self.M = M + self.S = S + self.dynamics = dynamics + self.num_layers = num_layers self.num_lin_layers_after_pde = num_lin_layers_after_pde self.device = device - self.pde_layer = PDE_layer(num_nodes, heat_derivative_func) - self.k_hop_sum_aggregator = KHopSumAggregator() - self.graph_moment_aggregator = GraphMomentAggregator() + self.pde_layer = PDE_layer(dynamics=dynamics, n_largest_graph=n_largest_graph) + self.k_hop_sum_aggregator = KHopSumAggregator(self.K, self.M) + self.graph_moment_aggregator = GraphMomentAggregator(self.S) - self.lin_layers = torch.nn.ModuleList() + self.time_points = self.pde_layer.output_times + self.aggregated_size = self.S * self.K * self.M * self.input_feature_dim * len(self.time_points) - input_size = 26 * 3 * 4 * 4 * 6 # TODO: FIX PARAMETER SPECIFICATION AND INPUT - print('input size is ' , input_size) - for i in range(num_lin_layers_after_pde): - self.lin_layers.append(torch.nn.Linear(input_size, input_size)) + self.lin_layers = torch.nn.ModuleList() + + print('input size is', self.aggregated_size) + layer_size_list = [self.aggregated_size, 64, 48, 32] + for i in range(len(layer_size_list) - 1): + self.lin_layers.append(torch.nn.Linear(layer_size_list[i], layer_size_list[i+1])) + self.classifier = torch.nn.Linear(layer_size_list[-1], output_dim) - self.classifier = torch.nn.Linear(input_size, num_outputs) + self.nonlin = torch.nn.LeakyReLU() + self.outnonlin = torch.nn.Sigmoid() - def forward(self, x, edge_index): - print('input data has shape ' , x.shape) - x = self.pde_layer(x, edge_index) - print('pde output has shape ' , x.shape) + def forward(self, x, edge_index, batch_index): + x = self.pde_layer(x, edge_index, batch_index) x = self.k_hop_sum_aggregator(x, edge_index) - print('k hop sum aggregator has shape ' , x.shape) - x = self.graph_moment_aggregator(x, edge_index) - print('graph moment sum aggregator has shape ' , x.shape) + x = self.graph_moment_aggregator(x, batch_index) + + # keep first axis but flatten all rest + x = x.view(x.size(0), -1) - x = x.view(-1) - print('flattened output has shape ' , x.shape) - for lin_layer in self.lin_layers: x = lin_layer(x) - x = torch.relu(x) + x = self.nonlin(x) x = self.classifier(x) + # try without sigmoid return x - + #return self.outnonlin(x) + def reset_parameters(self): for layer in self.children(): if hasattr(layer, 'reset_parameters'): layer.reset_parameters() + if __name__ == '__main__': # test the model num_nodes = 10 - num_features = 6 + num_features = 5 x = torch.randn(num_nodes, num_features) edge_index = torch.tensor([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9], - [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]], dtype=torch.long) - model = DYMAG(num_nodes, num_features, 26, 3, 2, 'cpu') - print(model(x, edge_index).shape) + [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]], dtype=torch.long) + batch_index = torch.tensor([0, 0, 0, 0, 0, 1, 1, 1, 1, 1], dtype=torch.long) + model = DYMAG(num_features, 5, heat_derivative_func) + print(model(x, edge_index, batch_index).shape) print(model) - model.reset_parameters() - print(model) \ No newline at end of file + diff --git a/DYMAG_solver_version/DYMAG_pl.py b/DYMAG_solver_version/DYMAG_pl.py new file mode 100644 index 0000000..57da3ae --- /dev/null +++ b/DYMAG_solver_version/DYMAG_pl.py @@ -0,0 +1,100 @@ +import os +import sys +import torch +import pytorch_lightning as pl +from aggregators import KHopSumAggregator, GraphMomentAggregator +from PDE_layer import PDE_layer + +class DYMAG_pl(pl.LightningModule): + def __init__(self, + input_feature_dim, + output_dim, + dynamics='sprott', + K=3, + M=4, + S=4, + num_layers=2, + num_lin_layers_after_pde=2, + learning_rate=1e-3, + custom_device='cpu'): + super(DYMAG_pl, self).__init__() + + self.save_hyperparameters() + self.input_feature_dim = input_feature_dim + self.output_dim = output_dim + self.K = K + self.M = M + self.S = S + self.num_layers = num_layers + self.num_lin_layers_after_pde = num_lin_layers_after_pde + self.custom_device = custom_device + self.learning_rate = learning_rate + self.dynamics = dynamics + + self.pde_layer = PDE_layer(dynamics=dynamics) + self.k_hop_sum_aggregator = KHopSumAggregator(self.K, self.M) + self.graph_moment_aggregator = GraphMomentAggregator(self.S) + + self.time_points = self.pde_layer.output_times + self.aggregated_size = self.S * self.K * self.M * self.input_feature_dim * len(self.time_points) + + self.lin_layers = torch.nn.ModuleList() + print('input size is', self.aggregated_size) + layer_size_list = [self.aggregated_size, 64, 48, 32] + for i in range(len(layer_size_list) - 1): + self.lin_layers.append(torch.nn.Linear(layer_size_list[i], layer_size_list[i+1])) + self.classifier = torch.nn.Linear(layer_size_list[-1], output_dim) + + self.nonlin = torch.nn.LeakyReLU() + + def forward(self, x, edge_index, batch_index): + x = self.pde_layer(x, edge_index, batch_index) + x = self.k_hop_sum_aggregator(x, edge_index) + x = self.graph_moment_aggregator(x, batch_index) + + # keep first axis but flatten all rest + x = x.view(x.size(0), -1) + + for lin_layer in self.lin_layers: + x = lin_layer(x) + x = self.nonlin(x) + + x = self.classifier(x) + return x + + def reset_parameters(self): + for layer in self.children(): + if hasattr(layer, 'reset_parameters'): + layer.reset_parameters() + + def training_step(self, batch, batch_idx): + out = self.forward(batch.x, batch.edge_index, batch.batch) + loss = torch.nn.functional.mse_loss(out, batch.y) + self.log('train_loss', loss) + return loss + + def validation_step(self, batch, batch_idx): + out = self.forward(batch.x, batch.edge_index, batch.batch) + val_loss = torch.nn.functional.mse_loss(out, batch.y) + self.log('val_loss', val_loss) + + def test_step(self, batch, batch_idx): + out = self.forward(batch.x, batch.edge_index, batch.batch) + test_loss = torch.nn.functional.mse_loss(out, batch.y) + self.log('test_loss', test_loss) + + def configure_optimizers(self): + optimizer = torch.optim.Adam(self.parameters(), lr=self.learning_rate) + return optimizer + +if __name__ == '__main__': + # Test the model + num_nodes = 10 + num_features = 5 + x = torch.randn(num_nodes, num_features) + edge_index = torch.tensor([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9], + [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]], dtype=torch.long) + batch_index = torch.tensor([0, 0, 0, 0, 0, 1, 1, 1, 1, 1], dtype=torch.long) + model = DYMAG_pl(num_features, 5, heat_derivative_func) + print(model(x, edge_index, batch_index).shape) + print(model) diff --git a/DYMAG_solver_version/PDE_layer.py b/DYMAG_solver_version/PDE_layer.py index 7158377..f2d13c1 100644 --- a/DYMAG_solver_version/PDE_layer.py +++ b/DYMAG_solver_version/PDE_layer.py @@ -1,48 +1,132 @@ import torch from torch_geometric.nn import MessagePassing from torch_geometric.utils import add_self_loops, degree - +from torch_geometric.data import Batch import torch.nn as nn -def heat_derivative_func(lx): - return -lx - class PDE_layer(MessagePassing): - def __init__(self, num_nodes, derivative_func): + """ + PDE_layer class represents a layer for solving partial differential equations (PDEs) using message passing. + + Args: + derivative_func (callable): A function that computes the derivative of the PDE. + + Attributes: + step_size (float): The step size for numerical integration. + solver (str): The solver method for solving the PDE. Can be 'euler' or 'rk4'. + sampling_interval (float): The interval at which to sample the solution. + final_t (float): The final time for the integration. + dynamics (str): A description of the time derivative of the PDE. + + Methods: + get_laplacian: Computes the Laplacian of the input data. + forward: Performs the forward pass of the PDE solver. + + """ + + def __init__(self, dynamics='sprott', n_largest_graph=100, b=0.25, **kwargs): super(PDE_layer, self).__init__(aggr='add') - self.num_nodes = num_nodes self.step_size = .01 - self.solver = 'rk4' - self.sampling_interval = .2 - self.final_t = 5 - self.derivative_func = derivative_func + self.solver = kwargs.get('solver', 'rk4') + # set up sampling_interval and final_t from kwargs if provided + self.sampling_interval = kwargs.get('sampling_interval', .2) + self.final_t = kwargs.get('final_t', 5) + self.b = b + #self.random_weights = torch.rand((n_largest_graph, n_largest_graph)) - 0.5 + # set random weights to be +1 or -1 + self.random_weights = (torch.randint(0, 2, (n_largest_graph, n_largest_graph)) * 2 - 1) / (n_largest_graph -1)**(1/2) + if dynamics == 'heat': + self.derivative_func = self.heat_dynamic + elif dynamics == 'sprott': + self.derivative_func = self.sprott_dynamic + print(f'Initialized with {dynamics} dynamics') + + self.output_times = torch.arange(0, self.final_t + self.sampling_interval, self.sampling_interval) + + def heat_dynamic(self, x, edge_index, batch): + return -self.get_laplacian(x, edge_index, batch) + + def sprott_dynamic(self, x, edge_index, batch): + """ + Sprott dynamics: + du/dt = -b * u + tanh(sum_ij a_ij * u_j) + """ + # row, col represent source and target nodes of directed edges + row, col = edge_index - self.lin = nn.Linear(num_nodes, num_nodes) + # Create a map from batched node indices to original indices + batch_node_count = (batch == 0).sum().item() # assuming uniform size across graphs in the batch + batch_offset = batch * batch_node_count - def get_laplacian(self, x, edge_index): - edge_index, _ = add_self_loops(edge_index, num_nodes=self.num_nodes) - row, col = edge_index + # Adjust row and col indices by subtracting the batch offset + row_adjusted = row - batch_offset[row] + col_adjusted = col - batch_offset[col] + + # Propagate messages using random weights + #weighted_x = self.random_weights[row_adjusted, col_adjusted][:, None] * x[col] - deg = degree(row, dtype=torch.float) - deg_inv_sqrt = deg.pow(-0.5) - norm = deg_inv_sqrt[row] * deg_inv_sqrt[col] + # Use self.propagate to aggregate the messages + aggregated_message = self.propagate(edge_index, x=x, norm = self.random_weights[row_adjusted, col_adjusted]) - return self.propagate(edge_index, x=x, norm=norm) + # Apply Sprott dynamics + dt = -self.b * x + torch.tanh(aggregated_message) + return dt - def forward(self, x, edge_index): + + def get_laplacian(self, x, edge_index, batch, normalized=True): + """ + Computes the Laplacian of the input data. + + Args: + x (Tensor): The input data. + edge_index (LongTensor): The edge indices of the graph. + batch (LongTensor): The batch indices of the graph. + normalized (bool): Whether to normalize the Laplacian. + + Returns: + Tensor: The Laplacian of the input data. + + """ + #edge_index, _ = add_self_loops(edge_index, num_nodes=x.size(0)) + row, col = edge_index + if normalized: + deg = degree(row, num_nodes=x.size(0), dtype=torch.float) + deg_inv_sqrt = deg.pow(-0.5) + deg_inv_sqrt[deg_inv_sqrt == float('inf')] = 0 + + norm = deg_inv_sqrt[row] * deg_inv_sqrt[col] + else: + norm = None + + adj = self.propagate(edge_index, x=x, norm=norm, size=(x.size(0), x.size(0))) + return x - adj + + def forward(self, x, edge_index, batch): + """ + Performs the forward pass of the PDE solver. + + Args: + x (Tensor): The input data. + edge_index (LongTensor): The edge indices of the graph. + batch (LongTensor): The batch indices of the graph. + + Returns: + Tensor: The solution of the PDE at different time steps. Output has shape [time_steps, num_nodes, num_features] + """ + num_nodes = x.size(0) + batch_size = batch.max().item() + 1 if self.solver == 'euler': num_steps = int(self.final_t / self.step_size) - sampling_interval_steps = self.sampling_interval // self.step_size + sampling_interval_steps = int(self.sampling_interval // self.step_size) num_outputs = (num_steps // sampling_interval_steps) + 1 - outputs = torch.zeros((int(num_outputs), *x.shape), device = x.device, requires_grad=False) + outputs = torch.zeros((int(num_outputs), num_nodes, x.size(1)), device=x.device, requires_grad=False) outputs[0] = x output_idx = 1 for t_step in range(1, num_steps + 1): - lx = self.get_laplacian(x, edge_index) - dt = self.derivative_func(lx) + dt = self.derivative_func(x, edge_index, batch) x = x + self.step_size * dt if t_step % sampling_interval_steps == 0: outputs[output_idx] = x @@ -52,20 +136,19 @@ def forward(self, x, edge_index): elif self.solver == 'rk4': num_steps = int(self.final_t / self.step_size) - sampling_interval_steps = self.sampling_interval // self.step_size + sampling_interval_steps = int(self.sampling_interval // self.step_size) num_outputs = (num_steps // sampling_interval_steps) + 1 - outputs = torch.zeros((int(num_outputs), *x.shape), device=x.device, requires_grad=False) + outputs = torch.zeros((int(num_outputs), num_nodes, x.size(1)), device=x.device, requires_grad=False) outputs[0] = x output_idx = 1 for t_step in range(1, num_steps + 1): # Compute an RK4 step - # ref: https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods - k1 = self.step_size * self.derivative_func(self.get_laplacian(x, edge_index)) - k2 = self.step_size * self.derivative_func(self.get_laplacian(x + 0.5 * k1, edge_index)) - k3 = self.step_size * self.derivative_func(self.get_laplacian(x + 0.5 * k2, edge_index)) - k4 = self.step_size * self.derivative_func(self.get_laplacian(x + k3, edge_index)) + k1 = self.step_size * self.derivative_func(x, edge_index, batch) + k2 = self.step_size * self.derivative_func(x + 0.5 * k1, edge_index, batch) + k3 = self.step_size * self.derivative_func(x + 0.5 * k2, edge_index, batch) + k4 = self.step_size * self.derivative_func(x + k3, edge_index, batch) x = x + (1/6) * (k1 + 2 * k2 + 2 * k3 + k4) @@ -76,20 +159,28 @@ def forward(self, x, edge_index): return outputs def message(self, x_j, norm): + if norm is None: + return x_j return norm.view(-1, 1) * x_j def update(self, aggr_out): return aggr_out - if __name__ == '__main__': - num_nodes = 10 - x = torch.randn(num_nodes,5) + # Create dummy input data + x = torch.randn(10, 3) # Shape: [num_nodes, num_features] edge_index = torch.tensor([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9], - [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]], dtype=torch.long) + [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]], dtype=torch.long) # Shape: [2, num_edges] + # make edge_index undirected + edge_index = torch.cat([edge_index, edge_index[[1, 0]]], dim=1) + + batch = torch.tensor([0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=torch.long) # Shape: [num_nodes] + + # Create an instance of PDE_layer + pde_layer = PDE_layer() + + # Perform the forward pass + solution = pde_layer.forward(x, edge_index, batch) - model = PDE_layer(num_nodes, heat_derivative_func) - output = model(x, edge_index) - import pdb; pdb.set_trace() - print(output) - \ No newline at end of file + # Print the shape of the solution + print(solution.shape) diff --git a/DYMAG_solver_version/__init__.py b/DYMAG_solver_version/__init__.py index e69de29..4316473 100644 --- a/DYMAG_solver_version/__init__.py +++ b/DYMAG_solver_version/__init__.py @@ -0,0 +1,2 @@ +from .PDE_layer import PDE_layer, heat_derivative_func +from .aggregators import KHopSumAggregator, GraphMomentAggregator diff --git a/DYMAG_solver_version/aggregators.py b/DYMAG_solver_version/aggregators.py index b1725c9..5aa176d 100644 --- a/DYMAG_solver_version/aggregators.py +++ b/DYMAG_solver_version/aggregators.py @@ -15,13 +15,19 @@ def __init__(self, K=3, M=4): # M and K referenced from appendix B.2 of the pape def forward(self, x, edge_index): # x has shape [num_outputs (t), num_nodes, num_features] # want output to have shape [num_outputs, num_nodes, K, M, num_features] - + num_nodes = edge_index.max().item() + 1 # Compute the moments for each k-hop sum for each node k_hop_sums = torch.zeros(x.size(0), x.size(1), self.K, self.M, x.size(2)) for node_idx in range(x.size(1)): + if node_idx >= num_nodes: #, "Node index is out of bounds." + # this is an isolated node + # add in a self loop + self_loop = torch.tensor([[node_idx], [node_idx]], dtype=edge_index.dtype, device=edge_index.device) + edge_index = torch.cat([edge_index, self_loop], dim=1) # get the k-hop subgraph for each node for k in range(1, self.K+1): - subset, _, _, _ = k_hop_subgraph(node_idx, k, edge_index, relabel_nodes=True) + subset, _, _, _ = k_hop_subgraph(node_idx, k, edge_index, relabel_nodes=False) + #print(f"Subset for node {node_idx}, k={k}: {subset}") for m in range(1, self.M+1): #print(k,m) k_hop_sums[:, node_idx, k-1, m-1,:] = self._get_moment_sum_aggr(x[:,subset,:], m) @@ -37,30 +43,56 @@ def __init__(self, S=4): super(GraphMomentAggregator, self).__init__() self.S = S - def forward(self, x, edge_index): + def forward(self, x, batch_index): # x has shape [num_outputs, num_nodes, K, M, num_features] - # want output to have shape [num_outputs, S, K, M, num_features] + # batch_index has shape [num_nodes] + # We want output to have shape [num_graphs, num_outputs, S, K, M, num_features] + + # Get the number of graphs from the batch index + num_graphs = batch_index.max().item() + 1 + + # Initialize the output tensor + num_outputs, num_nodes_total, K, M, num_features = x.size() + graph_moments = torch.zeros(num_graphs, num_outputs, self.S, K, M, num_features, device=x.device) - # Compute the moments for each node - graph_moments = torch.zeros(x.size(0), self.S, x.size(2), x.size(3), x.size(4)) - for s in range(1, self.S+1): - graph_moments[:, s-1, :, :, :] = self._get_graph_moment(x, s) + # Compute moments for each graph in the batch + for g in range(num_graphs): + # Get node indices for the current graph + mask = (batch_index == g) + + # Extract relevant data for the current graph + x_graph = x[:, mask] # Shape [num_outputs, num_nodes_in_graph, K, M, num_features] + + for s in range(1, self.S+1): + # Compute the s-th moment for the current graph + graph_moments[g, :, s-1, :, :, :] = self._get_graph_moment(x_graph, s) + return graph_moments def _get_graph_moment(self, x, s): # x has shape [num_outputs, num_nodes, K, M, num_features] - # want output to have shape [num_outputs, K, M, num_features] + # Output should have shape [num_outputs, K, M, num_features] return x.abs().pow(s).sum(dim=1) if __name__ == '__main__': - # test the aggregator - x = torch.randn(5, 10, 3) - edge_index = torch.tensor([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9], - [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]], dtype=torch.long) + # Example usage + num_outputs = 2 + num_nodes = 10 + K = 3 + M = 5 + num_features = 4 + batch_size = 5 + S = 4 + + # Create example input + x = torch.rand((num_outputs, num_nodes, K, M, num_features)) + edge_index = torch.randint(0, num_nodes, (2, num_nodes)) + batch_index = torch.randint(0, batch_size, (num_nodes,)) + + # Initialize the aggregator + aggregator = GraphMomentAggregator(S=S) + + # Compute the graph moments + graph_moments = aggregator(x, batch_index) - aggregator = KHopSumAggregator() - output = aggregator(x, edge_index) - print(output.size()) - graph_aggregator = GraphMomentAggregator() - graph_output = graph_aggregator(output, edge_index) - print(graph_output.size()) + print("Graph moments shape:", graph_moments.shape) diff --git a/DYMAG_solver_version/data_module.py b/DYMAG_solver_version/data_module.py index 8745654..bcaaf83 100644 --- a/DYMAG_solver_version/data_module.py +++ b/DYMAG_solver_version/data_module.py @@ -1,20 +1,47 @@ import torch from torch_geometric.data import Data, Dataset +from torch_geometric.loader import DataLoader import networkx as nx import numpy as np +from torch_geometric.utils import to_undirected class RandomDataset(Dataset): + """ + A dataset class for generating random graphs. + + Args: + random_graph_model (str): The type of random graph model to use. Options are 'er' (Erdos-Renyi) or 'sbm' (Stochastic Block Model). + num_graphs (int): The number of graphs to generate. + graph_params_bounds (dict): A dictionary specifying the bounds for each graph parameter. The keys are the parameter names and the values are tuples (lower_bound, upper_bound). + node_features (str or None): The type of node features to include. Options are 'degree' or None. + + Attributes: + random_graph_model (str): The type of random graph model. + num_graphs (int): The number of graphs to generate. + graph_params_bounds (dict): The bounds for each graph parameter. + node_features (str or None): The type of node features. + graphs (list): The generated graphs. + generating_parameters (list): The generating parameters for each graph. + """ + def __init__(self, random_graph_model, num_graphs, graph_params_bounds, node_features=None): super(RandomDataset, self).__init__() self.random_graph_model = random_graph_model self.num_graphs = num_graphs self.graph_params_bounds = graph_params_bounds + self.node_features = node_features self.graphs, self.generating_parameters = self.generate_graphs() - self.node_features = node_features # options are degree or none def generate_graphs(self): + """ + Generates random graphs based on the specified parameters. + + Returns: + graphs (list): The generated graphs. + generating_parameters (list): The generating parameters for each graph. + """ graphs = [] generating_parameters = [] for _ in range(self.num_graphs): @@ -23,6 +50,8 @@ def generate_graphs(self): graph = nx.erdos_renyi_graph(n=int(graph_params['n']), p=graph_params['p']) elif self.random_graph_model == 'sbm': sizes = [int(graph_params['n'] / graph_params['k'])] * int(graph_params['k']) + # modify the size of the last block to make sure the total number of nodes is n + sizes[-1] += graph_params['n'] - sum(sizes) p_matrix = np.full((len(sizes), len(sizes)), graph_params['p_out']) np.fill_diagonal(p_matrix, graph_params['p_in']) graph = nx.stochastic_block_model(sizes, p_matrix) @@ -33,35 +62,59 @@ def generate_graphs(self): return graphs, generating_parameters def sample_graph_params(self): + """ + Samples random graph parameters based on the specified bounds. + + Returns: + graph_params (dict): The sampled graph parameters. + """ graph_params = {} for param, bounds in self.graph_params_bounds.items(): lower_bound, upper_bound = bounds - if param == 'n' or param == 'k': + if param in ['n', 'k']: value = np.random.randint(lower_bound, upper_bound + 1) else: value = np.random.uniform(lower_bound, upper_bound) graph_params[param] = value return graph_params - def len(self): + def __len__(self): + """ + Returns the number of graphs in the dataset. + + Returns: + int: The number of graphs. + """ return self.num_graphs - def get(self, idx): + def __getitem__(self, idx): + """ + Returns the graph and its corresponding parameters at the given index. + + Args: + idx (int): The index of the graph. + + Returns: + Data: The graph data object containing the node features, target values, and edge indices. + """ graph = self.graphs[idx] params = self.generating_parameters[idx] if self.node_features == 'degree': x = np.array([graph.degree[node] for node in graph.nodes()]).reshape(-1, 1) + if self.node_features == 'random': + n_random_signal = self.graph_params_bounds.get('n_random_signal', 5) + x = np.random.rand(graph.number_of_nodes(), n_random_signal) else: x = np.eye(graph.number_of_nodes()) - y = np.array([params[param] for param in self.graph_params_bounds.keys() if param in params]).reshape(-1, 1) - # store the graph structure so a GCN can use it - data = Data(x=torch.tensor(x, dtype=torch.float), y=torch.tensor(y, dtype=torch.float)) - data.edge_index = torch.tensor(list(graph.edges)).t().contiguous() - return data - #return Data(x=torch.tensor(x, dtype=torch.float), y=torch.tensor(y, dtype=torch.float)) + y = np.array([params[param] for param in self.graph_params_bounds.keys() if param in params and param != 'n']).reshape(-1, 1) + + edge_index = torch.tensor(list(graph.edges)).t().contiguous() + edge_index = to_undirected(edge_index) + + return Data(x=torch.tensor(x, dtype=torch.float), y=torch.tensor(y, dtype=torch.float), edge_index=edge_index) if __name__ == '__main__': - random_graph_model = 'er' # 'er' or 'sbm' + random_graph_model = 'sbm' # 'er' or 'sbm' num_graphs = 10 graph_params_bounds = { 'n': (50, 50), @@ -71,5 +124,18 @@ def get(self, idx): 'p_out': (0.05, 0.2) } + graph_params_bounds = { + 'n': (50, 50), + 'k': (1, 5), + 'p_in': (0.5, 0.9), + 'p_out': (0.05, 0.2) + } + dataset = RandomDataset(random_graph_model, num_graphs, graph_params_bounds) + dataloader = DataLoader(dataset, batch_size = 5, shuffle=True) print(dataset[0]) + print(dataset[0].y) + assert dataset[0].is_undirected() + + for batch in dataloader: + print(batch) diff --git a/DYMAG_solver_version/lyapunov.py b/DYMAG_solver_version/lyapunov.py new file mode 100644 index 0000000..6a23249 --- /dev/null +++ b/DYMAG_solver_version/lyapunov.py @@ -0,0 +1,101 @@ +import nolds +import numpy as np +import matplotlib.pyplot as plt +import networkx as nx +import torch +from PDE_layer import PDE_layer +import phate +import time +import argparse +from data_module import RandomDataset +from torch_geometric.loader import DataLoader +import pandas as pd +import os +import multiprocessing as mp +import numpy as np + +def compute_lyapunov_for_node(traj): + # Calculate the maximum Lyapunov exponent for a given trajectory + return max(nolds.lyap_e(traj)) + +def compute_lyapunov_for_feature(x, num_nodes, traj_ind): + lyap_max = 0 + # Calculate Lyapunov exponents for each node in the feature + for node_ind in range(num_nodes): + traj = x[:, node_ind, traj_ind] + l = compute_lyapunov_for_node(traj) + if l > lyap_max: + lyap_max = l + return lyap_max + +def measure_lyapunov(x): + # Assume that x is outputs + x = x.detach().cpu().numpy() + num_steps, n_nodes, num_features = x.shape + + # Prepare arguments for parallel processing + args = [(x, n_nodes, traj_ind) for traj_ind in range(num_features)] + + # Use multiprocessing to parallelize across features + with mp.Pool() as pool: + lyap_max_list = pool.starmap(compute_lyapunov_for_feature, args) + + return lyap_max_list + + + +if __name__ == '__main__': + argparser = argparse.ArgumentParser() + argparser.add_argument('--dynamics', type=str, default='sprott', help='Dynamics to simulate') + argparser.add_argument('--n_node_list', type=int, nargs='+', default=[25, 50, 100, 317], help='List of number of nodes') + argparser.add_argument('--num_reps_per_n_node', type=int, default=5, help='Number of repetitions per number of nodes') + argparser.add_argument('--num_graphs', type=int, default=10, help='Number of graphs') + argparser.add_argument('--signal_type', type=str, default='dirac', help='Type of signal to use') + argparser.add_argument('--sampling_interval', type=float, default=0.2, help='Sampling interval') + argparser.add_argument('--final_t', type=float, default=20, help='Final time') + argparser.add_argument('--step_size', type=float, default=0.01, help='Step size') + args = argparser.parse_args() + + results = pd.DataFrame(columns=['n_nodes','rep','graph_ind','lyap_mean', 'lyap_min', 'lyap_max', 'lyap_std', 'time']) + for n_nodes in args.n_node_list: + print('Number of nodes:', n_nodes) + for rep in range(args.num_reps_per_n_node): + print('Repetition:', rep) + # get a random integer seed + seed = rep + torch.manual_seed(seed) + np.random.seed(seed) + + graph_params_bounds = {'n': (n_nodes, n_nodes), 'p': (.1, .5)} + # get dataset + dataset = RandomDataset(random_graph_model='er', + num_graphs=args.num_graphs, + graph_params_bounds=graph_params_bounds, + node_features=args.signal_type) + + dataloader = DataLoader(dataset, batch_size=1, shuffle=True) + pde_layer = PDE_layer(dynamics=args.dynamics, n_largest_graph = n_nodes,sampling_interval = args.sampling_interval, final_t = args.final_t, step_size = args.step_size) + + for graph_ind, data in enumerate(dataloader): + x = data.x + edge_index = data.edge_index + batch = data.batch + start_time = time.time() + outputs = pde_layer(x, edge_index, batch) + time_elapsed = time.time() - start_time + lyap_max_list = measure_lyapunov(outputs) + lyap_mean = np.mean(lyap_max_list) + lyap_min = np.min(lyap_max_list) + lyap_max = np.max(lyap_max_list) + lyap_std = np.std(lyap_max_list) + new_row = {'n_nodes': n_nodes, 'rep': rep, 'graph_ind': graph_ind, 'lyap_mean': lyap_mean, 'lyap_min': lyap_min, 'lyap_max': lyap_max, 'lyap_std': lyap_std, 'time': time_elapsed} + results = pd.concat([results, pd.DataFrame([new_row])], ignore_index=True) + print(f'Results for graph {graph_ind}:') + print(f'Lyapunov mean: {lyap_mean}, Lyapunov min: {lyap_min}, Lyapunov max: {lyap_max}, Lyapunov std: {lyap_std}, Time: {time_elapsed}') + + save_string = f'lyapunov_results_{args.dynamics}_{args.signal_type}_{args.sampling_interval}_{args.final_t}_{args.step_size}.csv' + results_dir = 'lyapunov_results' + if not os.path.exists(results_dir): + os.makedirs(results_dir) + save_string = os.path.join(results_dir, save_string) + results.to_csv(save_string, index=False) diff --git a/DYMAG_solver_version/main.py b/DYMAG_solver_version/main.py index 55cde7b..47de781 100644 --- a/DYMAG_solver_version/main.py +++ b/DYMAG_solver_version/main.py @@ -1,10 +1,79 @@ -from .data_module import RandomDataset -from torch_geometric.data import Data +from data_module import RandomDataset +from torch_geometric.loader import DataLoader import numpy as np import torch import networkx as nx -from .PDE_layer import DYMAG, heat_derivative_func +from PDE_layer import PDE_layer +from DYMAG import DYMAG +import argparse -# run DYMAG on a RandomDataset of Erdos-Renyi graphs +if __name__ == '__main__': + args = argparse.ArgumentParser(description='Run DYMAG on a RandomDataset of Erdos-Renyi graphs') + args.add_argument('--num_graphs', type=int, default=25, help='Number of training graphs to generate') + args.add_argument('--num_graphs_test', type=int, default=5, help='Number of test graphs to generate') + args.add_argument('--n_nodes', type=int, default=100, help='Number of nodes in each graph') + args.add_argument('--p_min', type=float, default=0.1, help='Minimum edge probability') + args.add_argument('--p_max', type=float, default=0.5, help='Maximum edge probability') + args.add_argument('--device', type=str, default='cpu', help='Device to run the model on (cpu or cuda)') + args.add_argument('--dynamic', type=str, default='sprott', help='Dynamics to model') + args = args.parse_args() -# TODO: set up training loop on this dataset + # run DYMAG on a RandomDataset of Erdos-Renyi graphs + # set a seed for reproducability + torch.manual_seed(0) + # set up the dataset + num_graphs = args.num_graphs + num_graphs_test = args.num_graphs_test + n_nodes = args.n_nodes + graph_params_bounds = {'n': (n_nodes, n_nodes), 'p': (args.p_min, args.p_max)} + + train_dataset = RandomDataset(random_graph_model='er', num_graphs=num_graphs, graph_params_bounds=graph_params_bounds) + dataloader = DataLoader(train_dataset, batch_size=5, shuffle=True) + + test_dataset = RandomDataset(random_graph_model='er', num_graphs=num_graphs_test, graph_params_bounds=graph_params_bounds) + + model = DYMAG(input_feature_dim=n_nodes, output_dim=1, dynamics=args.dynamic, n_largest_graph=n_nodes, device=args.device) + + # set up the optimizer and loss function + optimizer = torch.optim.Adam(model.parameters(), lr=0.0001) + criterion = torch.nn.MSELoss() + + num_epochs = 100 + # training loop + for epoch in range(num_epochs): + total_loss = 0 + + for data in dataloader: + optimizer.zero_grad() + + # forward pass + output = model(data.x, data.edge_index, data.batch) + + # compute loss + loss = criterion(output, data.y) + + # backward pass and optimization + loss.backward() + optimizer.step() + + total_loss += loss.item() + + # print average loss for the epoch + avg_loss = total_loss / len(dataloader) + print(f"Epoch {epoch+1}/{num_epochs}, Average Loss: {avg_loss}") + + # run DYMAG on the test dataset and compute the average loss + total_loss = 0 + test_dataloader = DataLoader(test_dataset, batch_size=5, shuffle=False) + + for data in test_dataloader: + x = data.x + edge_index = data.edge_index + batch_index = data.batch + # forward pass + output = model(x, edge_index, batch_index) + loss = criterion(output, data.y) + total_loss += loss.item() + + avg_loss = total_loss / len(test_dataloader) + print(f"Average Loss on Test Dataset: {avg_loss}") diff --git a/DYMAG_solver_version/outputs_plot_heat.png b/DYMAG_solver_version/outputs_plot_heat.png new file mode 100644 index 0000000000000000000000000000000000000000..b39e0475ff3a120418d4744a706553233b735908 GIT binary patch literal 40728 zcmeFZbx@UW*EUQjT`Jv(0ty1sT_T8*A|V|T0@B?Qf~1NP(xD(J-Hjk3Eo_i35$SID z*5&;?_c!yzy!Sib%s0HpsX$VxwQ zPgTeUE}fR=U+8ASWlAjWHnjPVZpSY;p*OX3usT#-C4@ z9x0#RxTi)UEsKl8_An@j8s}9I3C_}}NTZ~{;}Q~2fy6DJMgu3_ZN=lY$yx91?_#}~ zIk~x4gGgj>nBEga7GWYUf=FB*T$X}AWI|cw{_l(bUiA|Ko zhR5K0kghUE$!AHDEaI~Dlt3o*qx9dE{GK8{wut+gBxCKjH=RF;%)w7GEj?)u4`I*x zNmTm5_m$Dor^ESFbM2AaE^IQPB8<3cBpIFRq!(j4=a;#INHTi9G8-|qc@*f^#v0Un z8uKJrPSo(P*KPMu9qT!D^E^!N8iFi57%bwqouXzfbOu1(ZdAFtKo5)rXny&B+m zww;=j`z`%WNY&EEL>xT4t0E%r{mzf4=sQwZ>6#)U2y7cqU9Vof%AMGK?fUi2Xup%z z`ueJ>czFh&x3RZ$P57u^MnudXEf-mA|LR`z`S$JG*RNk0B_+>K44pqDi@G_kmZN4D z7v=TzXn8Gq#7r5&zqIi7irq;Rc*ac~+1Amq-F8LAs90NU?}3uiLd|j^GY`+@%>H{J zM0D}I7Ev7%w{$C>bVSkRKQ6_(tENgCEFFlYZ(k=9`a<3cEtq7n%;(tEwrU0oZqXff zs(5TwZ2LQ}?A^N_D{j5w``iw*vcBi1ZqY9Y;1f;u*CxbHcLwQ?ds8J-XJ<|5-KWt> z?;8L5o-pp%b!%m+aq5(9uXeq`^!)5}+-tpFR$ji;w1a5UXJ6lSWu&=R^4wdF`K{<$ zRp8$0*!I8!uFYMP;r6&w?XNshYyAeQ^- zz_lc0RwN|Q-C^6@4`xDead7NetFq6}*@G;fJYotG3VF>epUV=XAnv}AR%-*Fm#dkl zk$28IKi?UX+xg^Ss`r*c^c;i7{H39(DZ^VgZm4N!Xz=p#YDf_h6Kkm?3C&E}^=I6@ zna26Jq)lH_Q`1#W48raFk<{sC2g8IJJYxi2&i<5NW9@g_%DnDGK`pYoYKERm*WEm% zrJq{~Bn}mucGPs`n%q}Xa)~3P6;dT6AlNw=&&|zM&K|6%rKQb3JM8t7*nW^_zWZnB zW|{=V!{WG#hQ{uD>q0ejm&rPhHJAAVb2e4-muxzvI>+8kmbZ+C)mXm^x zTiatgQP+sAEyte~&%D|zva+&xzqLnD-OAj$xZ+qryU!Q8U`~J#kYZRSmXlUq0{N*cG z(7$$dIjtwg#bzOlmrXj*dSasI8Eq`Y`B zd0VE?=@XI?u?Q%eT7vN(7n@u>KSDLSA1=f_{%rL8*rz_+g=teHSy{ww?R#EPgU(#T z=H@2a%2;J|b09Xd>;^u&8nbhAw=IfM)Jstdb+vU$AHxH%a z9~GJfI_Mk7iFccdkEdi&df%3%5ChAT>9xNmc7C!QmxN{zI$WSz>awhfMNIGA)2|)N z@i^d%iSg#_7y9zdA%;s@PdPGcaF{qF#4yNZlHO_NT^ANkC~7zg%g^U<^YHW(nesh; z7TJMpG$|?RMiVBn>cfY#n-XJR3}|}tKY2mo;9T;iV20xI920x18Joe~d+pDr5N#AM zeY%{S++`Y?<$kiR56+lQ5|vZFqViF+GmYn`>r`Xky1H)Mx^>a&42eSpB_;7WAwNGq zA9kB!^ucV+nSzZCI~yR2qZ>URS9Y4YfGD(|61_P!ry|T<+MAbHCXTEhUDbS zmoE$1Poo+7o$~8e*rY>^`_-MuYxRSz#&z{R93+VFCrEugB`~lN6Wlxw+;A_GF`~7A=9#l3u8_%E_Og{XhZzslXF0b4}4UjG4D<3zC$cNKQzKeO^C$WIb_m! zK|6YU3eL9P-e?G~&A7&sFB6p4Zrpgu)o>7yo$CAk?OW1}%uJ8p=^=0;jZO|W=@*kk zTyA0$(?5YHB&DDL?1H+7!6wW6F{gZ!rHn-85wGOl$;hCjjv1|E`g8GMICHQL}K8cu#|%jF)x*&Z1L1eDLh9$&k5 zjpXv>7I>!dxdpgjj^8zi7eKt0CRq&sJ>u0tH^j;wRB&|_c>er3Pt->M$^B9}x%V9$ zZWk99TkUx!&tmN%oSJaAt=(g)Z&c8{ZH0`@kN-*J0o8b@;T+owDOG{&9 zkY`4)HNM)TN3!t#V~I&p1q^qpT>4|Iix(9z*x(=exRnlf3qBrB>mqPG12^V^>faYH zU;p!k8$6T@-_6`ylgcp%A=(C_y9c3~d3Uo1`Ip>IdXRJ+L=yR8$m%*PYe-26k73AF z@9(fZaO4TF0*mf^iu?QM{cJ+X$~?Th{Z{2|?d>Mqlv!9gxDPeR@`k=a3H|#0dln!? zrHax;r@x;pn`u<8;~Sqkgl4cUCZ+i7jXn@r54eaQf7?{9 z;LBr%%irMRZyPnC0Z4@u`8`Q^Fng2Fc>qTimSf{}QWEX{mPM*h`;SlRqvh6dwH%C$ zP4)S?H;FX?0*NKkOgo}Z zmQ?&^U(3_a3~3qG?{GD~jG=Dwv^<`P@h%moD3EL2C@gLulE-19zNVpZqw>(4kEGN< z*V>0fIe{lw%fKTTj`SsJYNP?66tJO0z$o04$M&z=kg{A?f8K?oW+m1E;BK_qi8a;t z=mwk$4!x>Z@Czld$bbf5y>uzX+_?daSPbVs{_QNHS!58|I%5QXS(&U)R8AJT?Q`Ue z1_)iqaUQ4j1)(Xd$M;0RXy`6jl#@a`pKQ51i@UAq(R==88WbVgVP)Ok>f%yYRh6=|w0!10q!f0|e^>QfpxgcQ zxISyOyJP85zhxG)I#W=U_JRO5-ACW$80+<`{V z5mda|RrXBb6kO)85f-$m@7~?Dvby>B@#D7knww+~A3j71&!NZ;NAzWjiPV;!} z+3Y(2G;-#mA;3C0ISGx5a*gXw5l=igaG&(vRwBKdu3zi=#KK~GK1PKHDb#`-5clG9 zKJ{5y*CdX9GpD%E;-WoNSFc{J7!#xS(27e;d>|teApiK$qqkm(*Ze6|=+RMkCUk;k zE&3cYe%$x8eO=TFr*7sB6qp@tt{XR)u3rxZOn2q-Wf`bbOG`^n%*=+|PvA$8Jg%{@ zTqYrzPJ!KMgz78g&3*PcI|_tvLU-@Lg@#2~T4x!j-GLDiug&;9(zk4I=x z^y1%2%(^HgeNXP{8|u-2LeCE&5Fykm$6Gv(5F1|W%r>?8Ty*EMOej`IKB~mYv!u6B zU-#tnRD4+d;ll&yOWhaVX;}__Bs4EgP|K7vF4wUZYskvFP@vgW)OZ%AxlE;{5t%q@ zQ&{e<;^ujRX4)QM{O#?HeAr-@<>A+_UI9SH?D35H{o4}m<&Tn{wB=8E zV|-56{k^t5rr!_$BfM$vs#&exw()rbFtEyDuIlhtBDz+9)(~4J%S4S!>x})YxHz}r zYk5#u{|z+_RRJoVws&t(fs*3YaI*E_k%kKj<`_n71Ry!0P5JE8DJUpx9qeE3WS;i4k{UBdv-vJ29UN_7LN0@& zc3$>uEEa0)XR|Ky3fqZq`fd}!gG*w|XD?sAjGz~Pf7`VEvr#h!fv{%9uLIVLB z6r0^I|6D`Q&;BWkMMCzZ&BSf1>p6yD2voF;6n8*MxP*is3>$q37$lM)J=M6a`}UYZ z;8O9L6GTg%3Bf8v&#~NS@E7+AVZ-hd`KUitUE% zmHYYv9#Fh_^Cn#58WR&b!o1jf=R)H~Kra@*nBvugN>9>IP{4I?aLD@j5jcU$A4iA1 zyIvc6r?pJ4AhoXishB~_pNLtXyNA^B3*_O!hQ5cvkOADC=N_g)#ogN4f)aG$~C@YVo0B!?nN)MjR8 z5d2aQfQK8c#9Wq^qH4$Ux%Q3@aO6L*Vm+EL+fXz09jww z`TX+Q!l_EBttlZ_rN_ia;6=qf4DDWh{ZYE_KkE;-C}I%&VigRDQZF^go@x@^2bH4wJ!)xbX?n98H*T!?XlQ(~tu*ZOXl-c$D7T|I68GJz++e@~ zpaxuuIa;g*E@vg2`|(phi~olYAApH_hfoLYtH5c^0|m4QG=MoA%7f$Mn4~1dyLa!> z$$J`OJrcNZcx(CE^;rC4Vju)Vk8X*G{1m^4jh)FB<|WR|9N5$(1=Uk7nqC~pFodn~ zIo@*^R1&;IMs^ocIYJ*Te0xW~C|=**V^aL;%^MjwE@x+FJ4=Hh&=WS##g~=wVG&XX zrKU2Zr>8^3gmTyf`A@sloUqiq+YzNmFJzDA&u=}9^#8yT-oAZXF@fg`vjQEk1VA_a z>X&Eu0Kf!(mDl&g-L`I1F4}s`M&K$lv&Z@AzV%58r1S({3nFMcuZf5p4K&)%3|oc- z90SeN#BjC|LS6YO_nUQ5eLJux-AaNJ7nFIvyR^ z(RJW4kp1hR~Qfm zbtHlpW?Dl^tx=r7wW!Srv#6we?@1Ot+!<7P85wCd(-KTTC-Qc6ta7m8861Yumg(j| zDoI~4SOuBv!D7!{D_SCK_zS`TLWlTrX)qg;PAZ8Hc13_jA%n&3{gfdiQiRs8( zu6A17>bpbyIyM%oX>+avAq|9uh1=TN&^|jaJp`(=YA#CDYC-hMBFgeglT#<;6hT2$ zCbh7TP`rfC(P!V&BhLghj2*kdh`^3g?OH)Y1i1?_=u*{9XWssCoKQ+D<^usrPmEOL&5q#}-wG1*C-c@9zWGQe0fP&J5*#p(`%t{rhkD zTy;&~6L^=Ukae=S`&)%%lPT+=b(KYv`i(i9YQ z0lIT>-xb&|n~~4Y0r#cs3IW5es`{$q5h|A+GrO~@33C6n^8W&;cp121{;}U8!g&7b zNqztf6SOtr>t~{FYfqp$CJ8wNw6w?|CmK>E=gNaC7M5kv{txcoAE|XK@)GCeje$=h zAt#pufPvH)kG1MW%dv`E5Pb;oaZNy=V$zM%!)468uG*99sn@oHhz1|^`suqHP#(`b z0M{UC;`C@Y6VjmA;e5=Z534vo|2sefQ5_6^NzC$5eff`h3v2$|baUeX)RwPZ9Dp!Z zJ8Kg(rK5_e513x&i3}6*=Bj-6Og|ZsWlHfixmG;1>rKm}`u4qe~ zY}|RIv5;VqGBe>ewE#_-p3m+O$WDonlkgkl)YL9DH#Y-qJN)=r+R)nkw=3fB$}HcsMgWy{>Kb zhv8no!l^*6wW(dDOrHy4mVJzmG{jAhAi~g58=}xau?h;tg@xg17V3GZb*2_kI#Xj+ zxJ2S@P4r|gjsm?Qt8P-PeXtN`I8tmHA}`3f`*LvU8W&eIAV!3%3BAM<Ro8NS^I0 zFQN+9tpBnIY7@|RRISx4Q-U;+DHplrBE%2L_Ug@>VTap`+rYP_4f6l2?Ek}mI82ka z{PS+c5ZjY^@)(!82=bY7)C5v#UswY&a+#E0{bDs*@}aI$i})U)5hm%ZMVeNZDw^hs zm{-RTLR^sWzrOJKc}+uL-a5;_ft~Q{0AgJ3ENC|+5+$XjmHnEW0yYLv(^dF68N)C6 z+TT)1H@#C|BNPb&n*8r!jQDE|VSMxQM=n`K`36nENMjxStV%Q$FL6%H3jty>Rmi(s z!tk%MjL24ZGnpdr%$^2{aJ0lx4)995@z?wb<1^UwrlduJ$>9m^tG#tnM&zFB$!>7& zAPr>>&+fR$93nA+V1z4I0vj78A(dSR(dA|D>{Eg)=HFOrM?`Ahb^TtwTatGW(kP?; zf?~pjo05XcwYGpL5s>^#po+RZmXh-KdaR*AOVn{OG6@7M&D!U~m-B|gnvJS{&jCp+ zAIGWhzSYV4GAztwY2X92$1f0N0}uk2?9T%QSWgryQ<79>)X8ObT_UoBzfg~Y>3;~s zxmY`~b_pRaJRm0^x+~y55JJx)vbxPnFHuo-0!?l>x4#1-O)|t|Uy-4tPNi)(qE#Tu z8W2^Sr(JflKpQ>h{mrL0#`)%pZt3IVZ8q`3oK(METqaq4b&Yq9Uezjv&xRA!-v0Es z%4_zG)WBnbOsCxHuB9a#)Bux_&&CMv4rCG}>W+<>R-m#dZt0ZeY|;YwG6((4`)DZ} z=(!U685!hrA-z9^bghbLah;VMf^=Uo5pc(oHk%ieqN0CW0abB=&k=oQdvCMG6m+BMEM zpwSq2=%AxGsfUy?0CnNPi6(RgkxjdF8QOLHiR-TP3!fyvN@6hup|MrQ4O;jzr0RLK z>eFltqI!f@|Mv)#b_@UoEn5p+K~E3f(ve#m>%U zwKo4>fJ$mG|#0KF8waYJC1%=B3RUhnt`*LbNBx* z4&|c08;AsmrhEbsw|np5moKlWo8%b#vKDhSXHLpkQ(?J3D9s0-Uc3(EZo*>=b69C_ zpy2+>hAFUk&=z;F7uCri(s^Qc1T{ZiUimL$X{dyKR^=zaWm%T%L_gEuw>5WK{o++K zIEPIMow@o)MPn##`MQ-+00ls{vY&J+05TH5G*riv5a1u&na5e~TyBz^eG%3F4%dj2 zn9yn!$!#r*q<24+cYNMe08$eG@*Dt3EInKUjm(Jm&dbAuAa@uTOg@gu=hUzHY&CQf zaG2(tYhbf~>x}TU6rjCB;OEpdhV=2{V`JrjS_R5na?t{ZnAifcGxXOb0A$wc_a~j! zhzkq=g@8fI51M< z!j9-;%Gulb$3Qvoce8AWC1NrByUGW`bV1fha^!b9H8FYV}u* z>uYcB!Z?J63>DG%S24v9`^)iaU*g@>vsL_um?NR%ef>jXeyfw-_iycOuiOGOR$Q#j zaRM>30{LLF;+2a?^8CUAtOkLMj0}K*prSigZ-70Jy_|_*dFfXOK%8vx4sW_BDEJ}( z5P1q~*F@wZD4TQCGQH`PBcbbnM%3i&VDm>4` z>>Z@1<;e?NcfD_Igw*&X1vo`=YNWO`JVgG%jdYt1d$vCA8(bx_Dwh&l$5oHu-|yo zy@92+a&%O#&priZ9ub})LoIOPLqo@vmmAhPQkYY!eg3yr-vVGBp?3%y9!Ax|Fe-q3XNS| z|J+Z|l|j>;^gBPLZp;QX#Q{j-0nif)t7h?#SOinV;Zk<$zk`WO_}qUe2PDEij5!|s zskEu#Bi=-mOf`!~8h3AT?RS=wl~xa@)F2(_{|WL8CF-{Zw5UiH0aw&qKtl|eHcGq? z9c`-@V!IB^69ph{s}pO#6g%q@g;$zy?;Ef7J^e7Ch11WMv#_x8XR?(QrTp`RJuKC* zAjukJG(Er5-L>7%pSdA>{V_-dCKbfzCty!dOc6Vbm@6$Qxu>Lbot-^=L|aR;7_uLS zVM7uK76$G!7r{HC(3>KjAnsLhTm=r1nm9r%OkZDLuyNIPfw&D4u9#iJAxl_B9~d!8 zN=rQe-X?%nO&z)Hi||Hpfu&0zUtDJJ*}0@5ws`>!sQ+=-v1jzrb=Cs_5e4ifBfgl- zoE@*bv1ZNADhTv_B~*-H{Gp=a`!X1K+2enJ;7?%Np#s@(R9I3Qp@UO+%Fv-f!9;~!8jDFv9c6tnV(gJmQ!c5#V zDIrMZKPxjHMTc+UZn*XHu0Rh8#(ibh=C~j_q zjRxGBtoq&jjllUI;L^xJx{Tfh0l(pRO%#EyKg*s5j@p?2$Wa4bl{y;4CF#7Ot>vP| zj3F!W3XP$v1nb+AgUR%gO~7%-^O)f+qlzq`X9e}0+o;KZJB3|Bg5DpfyJ>r&s`ImB zkCV-g-UGea6VMf+zp1pisn5J;Wvn_hOSPSU=iq(fL_T$v^5F0uwt3*m9k4(6hGqQt zX*zOr?ba>C!{`7G8O*{9MY(EnN;NLaFIs09xwzKz6*QqW(Jr?l13yNdLoM)KGxqLp z%k8F8z$9~7Lt~=cANZ9ESXcvfht=!)Gcz-YFabd_Tx8e_c*Cpb)OyxgqsJX0o$Zk{9-Ch%4gr#(y8|OM0zec*8WZ?^gkQAip#|#| zy5(reb?|I~0|FDAUvthP34Bj4L(ecf6T-j(+%QP5o8NhRH%3gOtEbLI<8ZIAaDzxX zFgQpB!j16$DF})O0CN>oR4l-)6E_C^W2sZ4o|}$39O$JILb8M!k6;R;#_J3UO#t^AZ%HCbm-~T zb4ZLd{?pnV5qKi10mJf)df5mRZvh6w@Dk0VaShb3K-gkH}DAj=rA&v z{NnP|zMAEF64c#8kc4bz*}n=}P!?x1+>Qd?-uO00TxIp>m+h+N2P3O^Vku%jYh+Vr`vuxxO*LrOiug!EG?V2qIbRLtefM|1q|j~*XmB}ACR~r{4_Jg~E=OX? zH4y&d{yKp*CFZNJjwy%r{@xlXVN-|X15*+JVHux*(MiV*d#<(~nDKl{;)UU4W7gm! zeyqo33piiMDt<;F3RW~BHD4$&-ak;rp`Kvyd6Bd0!hJx6+A3k4nLY}a_@YtJVtIv1 zDq*NuW{DC2ySL1LMaSf;;4Vzod*KogKmu$-rSn+=GJJG;jzMRKo=qpgfEi2?*9M1f z!bdOA-|u8wdg-k?T#J#FwL`pq@2B;Jix)lKe#pwI+Skw^tDW$bUL?fFXXfVa(=|5a z$s6)~+yPR|THUrX^lhpSF7fLlmk+}btjB##m(E8l(TKeg8!qzYmiT+qVe3(!vLGdF zygJex0BNN!!HA#&bRb*UtU)hWvVkiGLb4xhXP^QB9cczdXxw>F39QnWNJylL4C=M3 z9Y1tkNO z=a0cAqD3?6e$S5(qd+i&a0|`n=Dmm5gw$>SiHP}+@M*k=UmQf=GHX=Sm|f%PY@^ML zH0;awX=1;l3yMtk^I8W-v(&*C5_ZO&%{UnKTcC+z(9kG;Oe>LKwfzaSBX#$s4{`TY zRI&k(0H4ooGr|pc>3{%z zz!E4NZQYwJN)gyB1wAu-w7?Y>MWe&5h57lUKjE%GqrH9i4l=Y`(k`zz=xZla=cf!{ z-WxXa&BS=ZG5vyz16|GsEj1!wN+Mu{JslNdd zW5|jo>-w2lq)MUBdA@ZYO?kOvA*fnT(CyZuK+636vpmns2eJ|3*ZC4o!RxU#&l}YX zOe*x4;N~a1E3LkC|YKcrZ!=Bk~0Jm$f-CRsJ>*pd>grJiG+@1ZSeS zmjD>G>cCpT7A8~cy85-J=T|3d>XSb|K3Vb=ihP9CeNVV)?;Psf9}PF8EZygJU_>AP(C4g$oy6rKI={Q#XO+H)k&d(RK|CV&Vw30i-s1 zwiOze9B2vf`mP$O+heGZfE8z8JJ(9+0Zuur1mRwRd%HLd0G&g(A_4@y)Lk7=Ip?8^ zL`?e>WSK&upU`bMjwJDfb&We@v@%~|4@;@8T%_}KOC`4vcYbq7Q{rx9WL-wRaDolM zN<5p_vdV9mU!vCa+2;QKzP5^j!VF9k%)QexNC!bDzS~?blGde)mMYdZ(%2!8>=tlj7<4k&|}iD@Q@P_*Iv%nLC<0=&bhqF)zHR(2Ml z10`!sI7VKL7xE(h-RFsuj)a%uDlvae5zRHA4<6=y3vr8Fs4f~5)>G|yKjEQ=;>;xi zHzm&cL(&h!^sU3J|E~x|Wq;jWBa+8zpg>m;bW#rw4-2J|8t0{!T#dYVfoD3AV}PWI zLRvHPVL`p2UgzZgDZHro-~lA7|D-`ZhFeDpxbQx>Md~|tf$(fYh9UsLfw~~NSe#)F zHig@ttVm}W(5uI0O8Y5(@Qh@kcydFQ8S)|sU}kAo_rCyJ$vvfOsqyh&VaE{93HVkS zLH&uO5!{*$=h_184s+v&y2^hgM=fIZm6$ycuo(fch5fa+w*xHQFztXy%)q8Af@v4T zqgg~)DejIy^Q}I%IC^0ovO1jar7RkCxC6-u#o@F+R0^SGMkTFz;&rGp#p}d*?dJbk>RTbh0;V zRlFwUXu z6b>5;>KX-NypAi|vh2KX{YHo>YH)#pfZ9MbmjrC#jk#VFCQw+U1PMg?YiL6=f0#Em zxGW}iKeC&?#J?S5Xfku|EABmf`WN8X0}dj!Tw6RiBf9 zoiI$!QxWTT1{3c1JYCbrHcB)%F8^q|=gf^YZ<=n0Gw$7wp7@S+Yk$wKlbzj_cQM2# z;x{;FUW8nc`Qd~Wh#o0-2TaeGE?ol8aI}Nf-NCfno0@s=KzmK^w*;}df#{RL8(%c8 zJAte6IT<7p4P-=i80y^P+hlDXu2*%mZ8OgwOo6_)JsF=uu*BHfip%szGH>XRZ?C*&RG{cp%MVEc@q^BCK6=IapfXN!=#cEK#e;-UL za{*Cy8lriEAqwqIz?A3$c2bU%QBa)i1-ZU&e>8%e$2 ztx^bPR?~g&jZtOMOAprkr$f19;4o%kUtYPB_LzelOlx8v<(j97-F%1&u%4G+kf@>Z zHAIv0j%R4CZA0U>i2MPu4WfWFJp3^jjRB5*1tkS>8|5Kwaa!>Oa*_k*4j#`;`U zph&yW7E@52O!ugyPee|CZlTz3sS%Xc@UT>+L5n?2>1r^&Em>>;g%d__fZfzik?zgX#PqZcd%0Oc*(lWe)Zb=F3cc2|Gv#_+o)MRXiRAY$o59~C4siVs^ z=iIUL7w7~+zU57=?zsdZO-mJlk_ft@f`TysPs9@m-cf`#$Q-h&236jNPUp05NfX7! z^wCJD>^Z9IBS{gh*^Hh1!D(n24#KQb=dvWD(Mx#8!}wkdSaOl!9Q2Y+W~Z4Z>-FLV6$@9#$QVN??+{kcAnK^yN!GjH2bXyx84 z@m|{WBM{Ur8uFC`>f=W+N=W=mZylF(fbUr$j7iviI&y=YNMqKpN z)Cy={f_V>`Z=HWwsaY(}oY7gPzCT$*uDr7E^5pNc;=F+?WU?I)CVd}ZpPGrMRMHZ|9H8(PdQXx@g&#QTY!Z14yrHxMZyZ&_r*O{==b zQXhM}^R%i}Q_|6x^wuSw{B7 zzpvp_86Jy>!JA}~KrNho^OPS?{j873zr)S#>SELa$pqCmm8NXB=2`hBoj-F(<->6h zfhRMoX1N`3_N=8?B&qI%o!gjf$LoTs(YBtu7aHy15Xu+{@eE`M%+yGOZ`>Kf+!%_r z@s-p;7yRyD9L~f>q`VTiZ|Bcoh@Rcr(T}_v18%ZB&1oB^?Rd!jmmavzI_JLPCZWF{ zNW|RyG-TbPDa0@6q`rYr6M1Y{x}TeEeXa*oWj12=0lj|0^EdfITAneqd=-YuoxjVS z&2VzS^;y@4x(JEIw?mIHw1FGfin}a9LSgc&c=-yt`%I#96q&9b4e7;mk- zR*Wf0>~#ao5Wh|1HrOTE!G)m45j5S$Twt=3k+4`W;Svs@nO~L5NXZBtv?7Po7D_c_qrQ;;gZ`!!iOxEDJ z`Yp%A!(N1b{F;Y54Hi9jZ~!B=9-V?={Yl`fK@h8hJ7Wg?U0vyifu`cOD{7{BoFIA_ z(t2yXX#$^d9*<3Ah5G(lX>Izn)C3md4CRwHP(q8Q6-G449yT|TWH3eB8R#7JnW9)P z!N`vt;zatsAA$RF3IrK|-3sgCBjRFC7zx2nVSakQrED z6SseT2qr$HarsU;ddYtuWnRKuyRjllb29T{ak{zum3?I1>#skxW+GEQn@I>IyO5S6 z`{n$*&SMMA4LsnKiJ}z_fDu@*O|1MZ4~KDvVleC>H#Xw3nsH#dNZBx%(yESMQELCt z577GxwyaCg?QkQM;Kr(C8sn7)PaJ_WiWdzRZXNy5-$3@1d2xukDbcp(CCtfwh6xbF z6bZaBcn?$?w@i`MoiC`SkBA)AHdN#B0lC1ABrK)sh zpzVHBQ=_zf5(+Vy_g8hE>dIkUGH}&eC*J|;@?l=!r*PlJjI&9ch=9b&pkqIOfjfG4 zzMa0uD{>25@CzC6F8TtBn+%^hY#={1p-#bGevIV_A9XN=-BN@>-e*YdT!IYM2Qh6R zTge_xn)PL*=B0n_wQJ5VEL6MuO=Z}11X-BOUqkXu0G_wBD@SckVDttto5M^ZVs-LZ z9&nM#eZ@7NIJG|Wj;wZ**Yfev#!y)pPF2~JPd8U|FXov`^^j{N+_-RbM&+sE94Uze zRZ`%1qt6b}^vU7&;Yw+5&pi?v>O*vWp%8Vm9WTl4?XxRD#F;c_`BJ*Xsp0467W!gp zNfLM56eSb9Al18B_=Fp@0}KjSmfv4b$K==13|xjZ!1Vu(B92@Fyh3%Z@`5s_J1bu~ zO@3|Z$GOSeeNMzYo$g}4d`vQWDgXNCH5c~+d=3v%FQ zL7~Ly4|U(YWDnFym}Bm?Axz_Na<;>i5Yh@?VwQ9Nx3%gn4R_m^slg07rqSMuLzl%$ zu-eDKMm4LPzo=(dd3T*4NP<{qR#j_&4KhRLo8!v3D@`zX)W?ke!;%jy^?lhY3^I3U z&bp4aFl+e81h55r4dgwQFKJaiZLiy`?@zLR`|82FQnKMQY*$WsksugLq3eYN}o ze!cqH2;*;|o+@io1DAHpUNiK47)#V@*$^ClXzR2w-{&gEms7sa>AFQ?<4G*g8@z>~ zl`u(MGl}a%?;mR|gwZ@B0$l|I8e$pb!5Wq0Hh;Dnd_gs{KaK61hrG8c z>BW#%wST`_Ztm8`xv%L`1?_C>ckk1$XaR~xHC5cS2c7S5%)JR8ooRNb_bsUOnR3Wh z?EBKf9(v2HO&Q!nBc>N~Z$36~_(E3hu_VGOnbxKC$k=owKUD~N@AMNBqU)F` z-JKio(dbs&$6|bx$5|T%c>}x+nI=!AA-)My#v7N}(7~I-3Ua-6 zfo8Ph5|-Sf!!1;+*qNU$y8G^4;rKv(^?JY!2|Qhi=iS>d%3v|@{C1%PcI`O2TC!fU z8Z3~dJuJ||1$*S13czM)Gu1FLSKpLdXgp6Rr8e<>f$IC2F#}e?i^D6&eWhmcc9;(} zgv!t>e6Sk1{KI0hEAG1bb5`ET*45YVLQ)Yoa<$h|XIpWNFrh7FU_cML8R+nVFed#P ze7#%HC%bPCWT4-#W_p=UmkCvT{Jg&i>2AmTNiXW?T@*aSe(2gtms&-SZB1CH*~Xq4 zj1p|48gQ=Dv06{i0pjkMl*tud%KN77qOt=Wh(s;d_rgOpFtmONaeY~|#U znxS12Kk$Aclo{|SU;|ACw$TGwqHstKlko=eJJiD~675cg*~cPTrT7s623L$$F1vdR zjvr8$QC|onkFpKn=W4D+HRaJ6Eg2g3x*6qdxcq>pcrR2WfpG@O6y~9kR@xXwXvW>= zB9U>rjSYM71o9$Q2Rr>5<$Gi$Znw_GL$D<$37R5aCs!@zV*Rjtvg9yQPwVP~d8`>& zu=7KpUb3u#{Bnp)VlxTpa;;D>jGaR2AU;ADbC|sVg9B@^O@s+2{&njGLtcld%Gy3F z(#e`UnmeLj4n(rTE%NLm#Y&cGv={=v2&arOi&kB39uAxd^G)MDVgZ%K^6KAE+?b1& z`~?=F_NhjPs4PZF#&s~ru{MskFE-N36=}rt)?nWekv$a2e~_J9xPQUzR$x8*%-14! zzN`Jm(p`D_okn@mB(JZX5<(+?E9P%zgE4+E(;{C{0p=$`;G8YXBQ367JVYNN3{hVw zN_}R$Icx({EbCEA&1qTp3_~o5zYRv-Xi(aK)$Vv!InfT((h{mDwUSfUVjlaE+X@jY@%g$Sv`pdB z*)xH~^;I>`c;QoOACX0^*$)w{84^C03_LBvVrIIf*J;#+-j`8it#5@BOfsAC^|Y2~ z3Y~+U-4tdMzQgy8ykplIsPk|Fho2Tj+#Zy)p<9WqomnaZ+u)alV3DY+6om1*eHZOX z)5@k9rB{c{D~6czzUbK#PN``^D6haW5?-$>03wVy&gneIHTjiPxD;^8(J47MO# zVp38#?=-&*H!l`--^et4T)#}|i*;+eeSh!tm4nNdhqJk4-}NxDQPMrJyjNtL;@PPp z)-QH@-GO=bgTrfjkadgZ&-AfzjWoKaU`_xb|3agq7kfn45rq?|`JTkMJ&{B7L%HDY z<~EPUjlmRej0h2ll)e!5xjWu;niHk5HC8ypg(3Ma_Lf&hy2SVUe>3J@;ak7BS+Bp{ z7rQm0jnkjic{l7;^D~LuNMQv{t8QB43pf zCROn66Iu~DOz7bnD^i56;h}Ho<1taQR{I^NG5>Ct^_6PS(leer^JpybJ$G>a3z@z@ zgNaT$%~P{H|IEhL*WIt`NBMiFBION(G%!pVKeEfZeNx1$y_aVSh>3CE4#<%36aPQh z;SF-@Rhnu&nsAFYa0zzF2gSa#3_^jDme*$xJ{h!x5_bxj_>*bhyt{!Y1bJcyj zaP5G`rGY5V;j^xOba?cdy||o2O)G!0gr@@8U3}vL#*$=qrU*;y?#oPdk%_3X4@n5|pU$P^nKZx84tUYh8`6ys|Ygqx~YUqVJ*mo@sz?l@(6n z^_2G^mCu_jJvvlzW?UT4_Z%_e6BE0@o%gldWrzff9AGsTcJ;~kB3A&RIEwQrhDmtHk?Xas7zo~(2%oif$EW2Yjp zUWLzby>%R{;K&w~Lk1u85JOu>=Nl%mSQ7|CVHs>e!PRnBLXaNeYZVrg-L~KmOfRB7 zGhlWS=lF`gGW<63Ba&9k;DvOy@z2nxlTTW8d_(o;SW+8msRIlD3R2wVCXn@!#JSZ@W5@mED~TI}KR81)@-#QW zX3uLsEp-N8jx)UXVfih@OJj_gCjGHD`;^5zy8VMnp_JO{E52ABT%RY4b!X$*l`8{!LiPE*CzLE|7x~1=QZvN1VEyy0pte2a>dWei~fc;QwZV!gY z5L-7M9-h_lfXjVfDZ;?PPg_0F#~x;1IBAvNRlG(AV3b^;_8ce#JvLooQre&ZMWNSA z!dyb`d&8X?&gNPHl$VK{H{WfO7ka2B7j`upCzD;vV7c)CH+)2;$8pG7wd#kJ}>a+`ApbxTs@^-rF+dMZ>fC(cbb=gBFjv!0g zx`&uT;0OiF1k!jf4HqPKRKPcLM83SjIsgVB2^b>$+2EY~rje04#qk0@%&adurSTPTv0*}B{!6ZM?i_=|KnzrDzyW$VGpK5gmA z=GSBi>W4$jh=;+xh~7^^vh?#8O3{imc;iyFpI0?aJj@J(&bi)Cs0dQ;Q&i!L4k5x< zC;rqGse80@P&=6#+rl!M4$g41#$N@EGL}D|%*wqfl$=9^EkwI^_*B%RkR8vGFU?Z5 zjJnTKm9J2kI&xpvRqOi8yYi43H5|1af9!@<3!(VA8sqohlWCy@&R&y{s0l}3I1g2a zXe+EX)^h}C zqTXCi3Jw+u3wIzTpW)zaQF|0oT1IN|n5t)s=f{>~(g!o4t}chzxZ2?Y%^$g)CJtPp z&jdb?swVgF#?=2z8_KXOXA`0O?F|15)5Q0r3o|ejrcxz1QeVcx3^RysH)k%kLAM5wQ9U*4%Jq*&L_vyn2mbU%oaS-P#t zwIyy@?^hngJzl;-@#{X(#PZBL^?2)X^GJeVHd&?lvs>^rAg>uzPc58&f>++uMY~wo zMmt$EpNL&E+j`oQL@a$htE6u?0vgNE8&m7s_c1RM65;)Q zs+`pHSo7DiJq>3vQDkROyN&DhP>q;WsxK>URueWd{p^;`M0{a{k?mnG9QmxMsO4`{ zTU??(Rv`64*RhreGpRm6;gJP>#w!pKd}` zE?eVjE?%q~U%!ioCs%MOvnNzd zT^L1u`>9g=uKD*GURKFu2J>Vpzk10tPrs8L--D$SvEkF*v-mO=H^(Km8_k{@j1~DD z`5{E6z1_{RG2yoaDZG@4qGDLNL|h{qCEuqpmPt!nVq*4RxR{Q+A+P@K>%pN0i-0+w z?^$71yOHxyKGlvT0hX51`{wEQ&LVx5YR)<+)3I+lqGdX$x$#_rY|oGC<$Aj#BTGB+ zTH^OLSBdOL)K*;|)nX3WI*A)pjHx|y)l$mp{`EdLwm|cdt{LHnF9BAca?8uAhS_h% zv2)~*krkF~ZGC%2Cd7VIZ+pv8i^Je!ygq+4btX}yba;DOT7RJ2vQ%byrW%Lk+y95J zw~ngn`QpDXB48k(l(a=iBi$kKK?Fp&bP3Yk(kMy_NGsBaz@?FHlu!`qkZzEcZk~Po z{?_xy^IPj#%a8cEoO{lknLRUmX7+x)ukctxx4{Yf6)anr-jid#qcl|OsFka#@@wz2 z|F4Umi8uo*qrPJ1mz`AY=N*d{t0$Bx8-s8qVX?Z5D;pF#t(HX1CfXhub$2di;hkQT zqLM|pr(gw7{1v$t6N9!?u1d@KQe0=sBu~BEp7hjI6&@z07%duV_UuMWw^JJu^WGbz z))Iv-Ez3`6vN0j!7Hw%lx@DnCA3aMntXNaipH63I#~$!1#l9yd@_kIdcd?m7Vod&; z3S;8MxMXoGhXi~gig1Q>it0`Kx1D#mctfQof0h+oC6uG04B->FsVL#c@+{YE+1L_5 zXDn%~9B%Y)Vk4T9=H~@E+LcAwIdYi`%k|X{FV+3rHP_2$&gDbX$6LDNgl}E>b&eWW z6HE$}m9%=o9NpZ0@2O5Kk9um4Sw;^{&SS~}r+YC4HX*nq`M-FLLhn~PV#=-Yag9@R zf8JCRBy^6giN~f7CP@$PZ#e%ND%+7G^y7Z#v&kRctHMU;xA6uR(i>-0@ zwPBD;Vs0a;%+R=sIA%c-v+m+qP0+g6{bQjzFNy71@eBONxmSEx) zln}(F-WKq07*Fvd#=01H@kR`>wBKz%iHEmd9!Xrd_VOadt4ll=t2_MhAz9r+v2R^M zONE5Q9}6G1w5d#Zq7;M4ChN#?=5lzI+jE<=3zzk<_MdzJ^<-RG-+m90N|K&qkm#rM`Dw2}EGl(EFTqhZ;+ zrIUX{?ii7#ji5xpC6XD|+OtasFnf8yArk^KYdBrAgN#Nb^c0wCL7`T|Nt=&Iw;<|w zp=5Vhy_{yWwcL1g7+tiI7o4n=ako`GK=WMF=z{LH`MJ3fwo&kI1bwggjf}mq;CTM$ z3-Mc9CHQe0_GQoa`OV>&%4)r)8l7K8<{B5PQf!PjKU%I*Z^8mUo{4*uwtAC(p)elh zGuP37`Nf^5eAu2j{o#E$25w&Z)T}KPobQOyWiKc5WY@6G4bzrKefhN-$)v7`Udf8b zqr7mrH|5bS7s^B7QuZ8I&(7c(i^~}qr6gr}_;@$j6c8nFH@HMC?wF%~u0ly5=fXc01M;I|Ch zO;1x)%@m?uEO5`e@1AFIj!en)0v)v{N|{eT+_bb(x=CtdVw8iA9bCQde97I@5_{6;8Q{c z0`>#a^bu?uP`Rwft7Bnnr5LiWAt{udfLdTGUtuWt^BJ~~wnt3qD$gLTR!7%waaC)U zz2O3M>Ae7kUi$Q_dUvaPxoYpSz0&3xrliYKHTh9)ZOT5`5*b~P@bgo)1cTmag1BqA zWLCZ)ms6sM{GWiZn6ujCl3j7*{Asts(2%cVWh@H?dYo=+Avq86czkwiCi!dnlydU> zlbFb~LYa!=UgT-ehi2Kal#{8N7A6h03WbGB)Uj~kcjHw?w_@&?676dyXnvW&?YHE53YaJv8Bgu> z9CVjeAS%MHn*iJkoXye2_%9Y9v zW+mQ<9P9qsN75bR`Hd^gRe5_BWAf?gEq5){N|QUX+OyMEOV<*<#8~C5@_%_i(N7_1 zC19iP{qaogfmT`GK+WC&rAWL>Ktcy^wQ5I`LqQt3m5OSO?@)6tnhOwFFz9i&yKPJ; zN=v`y=kS(9_T&J$!!-ePL<9x~L=0I+#|IE4=SXNU)DMh|@)H}@-xaUz2)#PwJLLEI zaNhj7o>V2YXIT(llRz*{)6Mr%jucL1@d_%7#Uk0{Xo>DV9Qm^hF4Ji-bbh4%E(q;aZxm(gm8BZajmYFh4sN(ygS zE|;nKWh@iY;tZs{$W-2SdWyO6P&lMDJe|!YUxPlJQi(o-Qe!LId_XlgGeGrNy&%P* zcHyVZ$O^6Gahuu7!8bl1{d}xM&Sw|eBIp(YuE@yB`W9$eX)qW=01hz0=E{cjw+Ayh z{)3&@rGZ>>B%1@JlobuMEn1(0XFyVOoS3EBHXsTNAt-h}e;T)L@JQVuu0{l={qx~7 zT^cvuQ=?WJykKV9r8 zT7zFSv;674XV}0Z`cz(hmqW`a*S+eey1G@ErBYUEVkC+B21kiyod+asw&wBI9DU3y zw;lYzWu@m7{p;R_G9M?V)0Nr6d0)=Lfj)wiOicm(Kr5fH!_oxmtfoO>-%$of$;zZ5 zo-A~+q~*fiRch?f(1frh#pa(?Wgfp;z0T{OK(!!wqn(i12B^qSwY7;>1yuoKW4cFs zYaLKg8=VR%@j!)x6p;CxtYpV{3m!^o1>;FOGdxGVBu zXXH7LRM_eqIR@g z)}m~Z3T@cX1(`y*+#!E*ORg?i7lE-KhP%&T7a~h1_l8`qoY`^#!Jx8^8zHI1_D$wq zUft&=?z_|XSQssKi(*0@PqJYFF58#z&c<>Do>vxvkJw7=Ry1&o|2G*0g=53 zDLO_@MkX>cQe~ElEEFenUX4vzF`|Mhy|zSTA**X~dpxCv?xh9EHH{YCmau%==k}pZ zVOx6jp(?-&gsoi?cCC=v=F4tX#WqwI0`PC!eXa9HDV<8Fl3=3O1SRXLmYm z=FQd=wyag()O{6|3J+Ok-{;FucL@to6?_mbZ}l~mw~|=ee{hPDP0rJ@nHv7ZC7$uv z*x7A;oqgzLx7F^tiT_|(fvWSM!@yy#p8V*?;rY4Bzk6A5_Jzq<1Bqi z%pvNFmh072@ zq96Mz&<9#NgrT1mGt6}(z-clp7U$8 z>*itWka*<(J55Mw;JNm%>xuD9OzVHEav8%s6w8WhSSmEa{^*s}G3VqZXsj8ySnqN* z=9a@yCZ{#ow&o$1rb2|DO-njSW3vjf1b5JmV zs;0R+*W|s(odP*Ebc=cxVbqmq!m97riqp%Mx~5U7T2N5Dv7Vqv*^@7_iGFnojz z4`j>Z{f$CsjZ%OnA($$}7ik>0^a;*oAQz}gh3~>|FK6z;ldUVG+ZEUalmSg_xZLHl zn$@z-yAICvnHc6wt3ClPH}*`cq3Vt>AR#_jR(}<6?IAg(Hp|e_>3_G+N8h*bk2*0i z2@iNnJHHDScMX2_Vo0-^k1%e;@9LB@MzdoL{tUAcWr=*u>(6}hkxhCI(;M#>Up80yHFvl|ym2fLc`>!c+e7^5grDNG<2^?m=7^-~%c^_9T*th2YVEbB5pw>LCD=hFPrHs-eC$#GeJlX`%Sm znmdY%rj>|}WMg*8t3yjgayp7W@x3tl%7k};l+10zO_$dZG@(^yUE0dcLM=GkC!I;| zuF`UYYr?-EM9J&elf6oY6p68g*q&aI{-0B>k@%M6@{{V*?{~s#ro+Cie@gF{l?k?D z<|t0m<|uo1AO2Jr-_K2Kx6y0|0)d_i&lT;YSL#`k8GWhwj|567Bjz$Nk9($+WHFX~ zN(1kb^mAP;mC$R{mBN$gv2!QWU-qpAn+5GMZftiUa$VH9Q;+v;mgQ|_a5yJF3$er} z6^zJdZTFXG=c@yq(PmtwtYI00g_g``c1jd#{iD{Uosok{9YH6WsV}JVo9z+$hQfRaMt~(+o0MmOSVu!@rC&f~h!jtWI#xcP zIB=#hm?m(>UU67iF@&bZ(>C|9uT5+;O^xNZg5jbZJ~#IPNA%NGHR^S%TmeakIy+0Z z>R;Z)xN(GeEp~cpe=YerA?0Ktl%qZsVAr(s4QEFz2qot~GVggv{C9S{=nddagl-K- zB`H^h%InW@RIJ$+8b;4L+WycPJB@TRQQF!%aYxrE)~}L9Tub_r{_54?`Rsw`MirHd z`N_#8S{IgQZj)hYGND3YZANaJv+G`Q)}0GANpklqk{(&V6}slF6X;~~^xhDIOz)lA zXSHwhF$E)4Q5UY#GqS~yYBRBvRA+5jcEpu*tXN9r_8D?o^-z^9DGK2;u;E4%)a8`? zOqDU~Pvu|%9FcCFnmmD`~mO%kAv4Yu6%&fRm!7Qh(UC%Y54;Sc=!lg;Xs z^2PD|^UqTczqVZLfoNS_-qM176&qE52QQD?yc z35LIdcE;ka_$Y5{*0lYxhu^bx>U!l@TbRF$$d7s^4%BQN+`b#JD%Ev1uPt-HovO&~ z!T7q5<#x<>iD=~+=iwxhldB{!J5KXDhg^{L6WmoT2&aB@)Trl+3w3Aj?ayd9q`VS1F z^G1&Enu7d>1P`Xjm$dng{yg~ZUUg{EzJE85s*=6&1+FH3P6qp2zjyI{FXnE z1q>rdiOi)0Q*4O6grqnCi!C&>Q(SLeeeYVybL&wq<& zVKC5n1bsfeCKk+Gf4=eDmLl`bpi^r`pVcL6Q-7Xj601AlHIbOZoiw-c;8xuZ zn%u+xo>hez5GETw{LFHD%%^hmEv1xSCD9i+G`{jF^OwXsnHKln_Tbg0;b!rM_eGC| zUo1$7-d6E>Eq!#I+Jb>;?6F~GaH^^4&zDCGdS6(ZGaUuc%*@0^S3YJqI?^j-b+MVT zRMbQjgjgOYCU+-|xLVfiZyd!#zzL^0?@i|dd4npEayy0AJu9u;Aq>a9odG8LA`!Dr zqGU$3*1-|poLZ7s>GoYyo8&j07I&L(e&zN~dG?*Op3b5ro|;iPYCmO20D{3&G3Fm^r$j zeQWvhN8zU2G1s&*g9+1a3r=+gZJ7LxO}0$>EjX7;*G;;r29rLaV<}3L+tgyK;&azj z`&{cj+fr#C9-ZmYy*O>jf|%F8g}$eukj;mB2U$duXM6kr#-p_}s8;D<|FQ6=H0yUv z&1FC6vF|5*J%BuzE0jdSuKV8v%XzEHl;z0{uf9un&iDA8J1O_~(}mDVb8%5$D&F;; zzBw=>jgC=$lzt`4O;@GBKO@9nBrMu1znv>AdRCh>P@5p|YL5v^^P2>V>kD$Q9aL~> zKUa5;k-h)%J9$}FBA16~G^;D){Ge*EmV4P;o`?Ic^xDq7^_^49UJi4vhw~qn^*Q;% zmQUibK}mM{E^j>x1`PW(CJS1!#be)B?FpAzzZ+arOwb^1L=vhHxfEB>KPMXqbw zv|gLY((Z7@;C#8WgJ0WV<2+%}$0z)T&Slf24-LCVR2d`Tmc9J*Av?~9&3bi^HS)^Z z?eymyBn@7uE4Cj_oHHSf{<8ncYD7cRScCbj3HMIHr8`?85_bwO-AR|ZTiP1pJR$rX}pg?1P`v$SR=GI=b{fww}_{)a1#!x zE)#N2nRWC1nS5@Hj*Zd6Flz6q2CGBszbVu1)GX8fveg&Ov+qOD&SKTpuB@=&mB#$z zfri5O#{+Sl1A}2gT{a9!i)XTRWYNF2MP&EI3TDfqUSFo$T)sO6bSD5XwE~wxq#d|o z9YA*0;W)yorFPMKv#OT#b3`PR(B~ZGS8!d=eR`SP{JL6&_uv5IuimFy`X6PQ7*qWE zf?K;HOrsfzUAxA~MRj$&7o0J3G^{iBN0$jQ$?bW*TnZfDUdZmV3bG`Xe!N0-M+f)% z?>j{2@3@6aMJHXt$u#UYUr=Z{ zGg*wod5quZ5^5O5H9j<`gf>V9HmC$jR=WstxZQCb={Q(-YqTa{UNx&z59X5wK2Q!Ju`z%*&m!+S?d=5%uE!XaB7KStbNefz^GP@6mTmQqaBFZumoQ*-&X z+n-%7b!uNs9njaG`P9SIeDP!OWj-f#?k0D+#xV69nfg0bO-WmoDl2_K8aB+WhUWIa zPZdnIrslr?a12jmZCDc6(Z=!3An^4k5cek&_a8hJ$o))7Qm5`OnVZ^Itxr8<7295g zt2aMft{wKO+){@=eB8!hjV^ZX8Co-(H;{b2^V&@LLXM{l~O3{V&(7kPOF* z_;T~lovMW>%LUO8=9v34&B(1(-V9%2D^c1o^g}2@;5}uwF#o2#1a640f`Vo zw?R7(Sy>b|ql325>#eMXtPz0?zVDOT+|q}ZJ`CDXe$yd@ur{%bk+r#oZGyuui+xB>T8RLSQprA9RMq@VJRfnT9rX0_xJ{Dq4DW(;Y-kVIeZ&)Gc zvA)^U^;m|AdUdpsJ+gD=mer!egNivPL(h6SwiMY6q5iou#29Uw7_GG!?bR4hHU{cM z$F45+8WrCrLW=_-peud|f)?N`+^u+IX4*+*35>oTQI4L1mzyrv@*c3wU_&HaoHA&r zOn)sDI2v<;qu}m~7l}@Ps=czUEp8e=U^U^}J9T~gldVl4Mt|xrd&+$YyxDZut@7+QPtnt2FSeIRY}^#9=_?mLjM$v{!rmWdww=QnVqdm>O^`@ zyR*SCTKuIh)o;u!OXXy(0`=Z=@77*sJQq*Iw!?GsU(cKzH}_{(N8We5Tsst4%R2jw z!^<6MZ*noeT9#;BNK;<|as@b|WpBB36Q)+a;O6g~o85`Ps$~1im%O~K7ESXj|IS_i zZ_#IO*@UT_h3?S$yx+XMpJu*tzG;L}rT5&BceN8~ukUoHjd1&BbZXIeH%eVYmWXU+ za>|kmVupFF_2nVQqSM{ZJS5$)BE|b4=X0KhGly8`K}^Hm0n5H~W!4V{JeiIL`s19e zCmkI<20gDTHO|5yLfjzGMZblSdAN%2D>g-ZMBuU~9+Aom(cQ&7sG+@S1Xr(4})j`@+2@hN8YzVg!&M^eeWmopX+c6#QYQ!leuKDp@V!le#NfQF zc1|}NTx}-m9E8-qNf#|iUGmbOzROkap9;TWdW!Rme2&>ey+3Zl0JvzX>IaJl*$&0WL?CF`;3X;V7F4;#OfQT|HmB`LRLtI}{s55+*gO zAE_%?iF2Lpf;V*2=$qp&DN#nW26vY+#Ge&B3Y;kFzNflTiPHL$AZ^J-2ZcR76Or_bgw`- zSPrf|AlQ0a-Ut+)xSe+tYZ&0%+DY~d^02%T(z^dJNsva#)#kg&Y}8{nes)=3nj}2l ztd7)$_t@9O!Vz9dr3l!tIPiz<64iNbKuF}Oc%j3|xEt0b-zCXAwC>G|zC6(57vj8p z^h+aNjXh`_a5Jb_(dn-DeH={{2YH9mv-EB}Dp_EG8C?w$;BKIZz(S^p84d)P1xj-L zb`92M%}uXQVg9agmC`5SwaD|cgXK}u$fIPEAVSA6`Hb<6fj3Tw)DE}m<O%&mMsvG)n!K;ZWqi_DLrF;tQqPMg|^{Q`6tP z$&_V+!*;LkmvFP;p~7U-X8klh@U4Pqdn;+jL1y&^MYp95Vh~I%@whuOMb)T z2KQ3K123qdV?`UNP1Y7YnK;U^y1{5q<;q&br6;zDp@hz$%2b;^)u$}Tx8j0IKRa`K z6*#GxHB`t{{P)|5-s)%~o`iWR9rjtwe_l03LSYY~^$~y}Nm+sUp(A+Ac)Cd7$jhp=xuu5Ovn!}T?ZIZ{QZL^=4%j5YF@&t9W^1mY)`hOI9$FcORIrnvCc z>g-KK=N;KI${doSoN1`L(EMdmX75x)pb~=m;Ds#T&~ryW+!y!LRuAxPsBb>AVh9;L z{t^Fe{rR_<(~~_tvTPzE)90JR1U$-Dc8|)COpbs@kSY6{nx39I`GgZ9%m8L6J+@+( z%5}5WXS_F>EUw{Lxa#;SmiZ8WLav>1%p!ee!!T_ciy|9+FeYpCtY^ z$3Cg^8NXX|^%EIE|7V)~8X1q%YOrEyW2lZw*0xFS_^8<2Cp#Vq4PsaI zq+e?kW?`6hG-!;*UwLOM#co!j0>rw5-~A_EsIP%!6CVN`FB&^7uNRX z4t-G{KJ?KKsIhBk7E{5e94Dyb*#RrM%?|oUGf|p^)VFFdHSyx(l1nxf?zamJK;=v~W_IH&euk7qqL6WYS8M5RUEh|v3noveD4val5Ubn5 zkQ#g4JXDtDTNTm?<3Xs{q3My&Nu+UM9j*jm(`SeBKWY@r771$cx-okCIPmT*c{*8D zk9SK=?1O3FwhV>}cr_fpkYu^zaJ=?M9^j z(hg0sh~8{CSE+Vo3+ef9wZYWGoh zqP}o7SwlUJ&Kv*APek%(NHq!FyPPWE|0VD3U1hP$1bCAC6aRG=JuCC@iJ0TydM&QB z%={w-N}v)5HDe|qDRF9?@!^7rPxtrj0ZmtBGs+jaQWpjvDw3f4u*Lf^O2wqi`^t zvgxim*XLS=h`h=|#EshqIs9qfQWW-SxZY_lkI1>}!b`$4K{+ zJLgB{`H6jh>`P|qXsA1wTS-zETLWv12=Y$3Uc>Nn$>A@+)@H1$s;E6z*cy2U3;FKN zm2~A`n61t#KKeJ%$z=^~wu!slzSp}^(#3+b!w9JB5WDTP#6Rh1=XBdp+U&e*HJ{kP zF8uvh*>9#Ry1oW->BN?WEEG^0DP1vn)W@Qzl)QVoaaW?fohW)a!9?{?DI2b%-sd?GSs3b71?~4ppRXr!K^M zL#`LhI>?=~XE*8=Dc9K^M41vQ*Pm>g#22k5tw9=!ji+qy{=#4xQ6_TZT0nR)d{2mN z`Rg`;bIY%jn%RyF%K^Bd6Oz|Wbuxx2C zFt~8O8Q38;HaBhDXP75zUGG5 zi_7}sDIg?mxrfu;f3_=E5Gm$CU?mT$E=v~ae}zg@D7<2M@w0NuFQD^|uRz?L*xOd& zyg)+L&EIAyX|c0m`i5LHkdruD-x$tUnh}28v(nmOg=DnC`op64vLOJzK~J}y;9G8Z z-I>s^^}+s1AoN(hE%!l=v9=!>CI@&_l|Ej0JP&#U;X6HsyualdQ5QmA~l=P>mj3mIFU^p88mueNU}ffzUt>qU{R+yrGS9VNST$J>>(bV7wFxJ zNC2TkHlw)Jn2qczS!he0{y(u#kvDJNT*MhDRN-v*IgRzt*uc$?Q&))B-O8q$%v?9z z56HMEmrylJSJy+8*+Gspm(KqLJ9f9u^fV~Z*QQ|8;TLE({i;%iKv7H9jlL{G$Wnx; zr*|-+ASuRS8y+!$W=kQ*st5ype7gCSAS%PdybwA$z zXV8I^=nF5>5K10oUR4695D16Qfg*phTj)Cv(tZp)JpDt{dSnBC=7w|6v2I!a6kWI}J?*Gd8xbluE<`(uYvp=;nD`l%0nq6&qKXQD8r=*OrUw~b%0>$D2d zDoNlnBR%eG+hi!;!3^Cf87+(1+W8t*YL{aErHE_U12>V8cG~l4Szl00o|ZqXAbHNt zUQf6T$$WAYMeh(ESST8|P{6NK+e#ehY<-UOeTDiTIFJc;i#ZmD&zO(9*Aj={|C zwns)>lksSWIUQnncR?iNu1fhKf4lD?sZ1fG_3Ch_SUn660WHT^OKw(FuM$333IDHW z5d+id2D;2wl0Qs>&gHM(?@!)Uc<@|gQCS$CzqX0=Ae9KfU|!PPThE+3-jGzeMCUsx zL+5M0czN$EDlol~*Nq|Vdxim2Lb?hJW7h@FCXP`T2`F9_547z5UIdqDJ16(s;xGRM z;+{K2gR{tdCw+=2mY8+?NVEU*6!n!+->bp5!Q1WK*ILJuSNJ=$D%+gy$Sa%c!tgbv z^U)DQ9{wM;w-KM2K@N1^enJpX)44a6eNToD zXp<|athGLKx)D$EXf&Bz{COh0%xy{rsKugJhubPJ^px90oDZCGED!(8Kcqf66`;)j z3a^vS)VqnSEzFtW>8XkH^CVs;mr>Ki>owUIf0;UL_L4=@W;?jB9(sWNf2x)-qMW_P9D`)%nQBC|4 zL3{!SuOdcwod_696u6H5EG2ursrdS@L9uT^S&h27vX*K#Q0$n!E?bu(2YCjR=vnuT zjFbtMzlbA}{tRXeEWyM13q`*3QUL|GoM{4HQuvJ6oOl?YJ3nsx_lrItJ3nh6%QrF2C}9B z4iDcYI2L9)Jl|UvWo{Qb$3E~RdNKC^`_6YXJFfnB37nQ-O&*S6%iOZVKBbx_tq$}v zQm&T0fKm3Qk)Hm+jL~Y#kF+lN{vR_7o=eFyg+vATy<+<8C*R&e6@PegW}ifPo#JW& z#p0>ZqF|PL^2WAE$|1g)*S`0{haa6ds7j?2(;Xii;jYwtG~VP{!KIsHH8t&^SAOT} zH}@U#5~wBK zQ}>Mhj!6})OPk+0318q^ROgDoy-FnT{Y(q~&PMd8YvAb+5@x4R+w_Z(E9o^=5xYG3 zoz3aEb$49%HVn~Pl(-VyZYC>yD8H5d#nh6OldFj`7Z0^6e7Dwm7`>*J{|x8v!gZ1J z&rm$Z70k@`b*ZSL;+HG^IZ!?=XNt1eYei^}TaHGC`(*zC0)#9WtNZQvhf^n7QA*c| ziKKAy>(mo2$BdI-Q+E8LLwo;sHxnZik81VGkV=Cw9%|P2FavWyr7JdhQ702G;O*ay z)@69p?@?T+gq#4qnp7+6*GUcHiRmzG)&`~8?`xqNr*6{iIIo{`yW{ zTjE_iPc8X#VuE3B2w4vC7#_xOKBqQ2hnr;^Lyau826H%ksl3T7EL;mQK;?r^ben0WQ#-m1a*baHp4u(O#3T(m=Eh~jy;=Y5>s^$JHeEWc*>HMfib zHnZXIYd_n9w@@Nl+}exQ)9XLaL?1L`+fFDk6a>JhtXSBoIJxIXLP~T6IE}x>p)nuO z!bif@6<6-~ooCoC5vOvw2qJHN`zDG_tuI!Gsk0ClCH?pkj|~^{u?q6BiYuve{)G)g z@`Klv3DwkyN!-HsBU_R+<_;m$MS({}YCd9wRPqcbN;{Q9OfX{5lt$!aV2NtvkX(is zTHKo)JN1o4{qGQB!3wkUJIPo49%VjokU<_{0gw2Z^N`oPzUK-#A>Pi_ea`CN7lryQ?M|_L)bN*0HoV4=O8iI&LF2--0ri4Bt4vg5u2WeBg8 zT2^(f1&{p#>Z`2NzlQX=33sqxwr8`|vwAa<&v@b6n4%s5V(tQ9?u>T?BAt#GBASZp zi>5UiK!>r1uo-2(V9C7faOg$sG5q@$&q!PUJo1Z&epuZY zqWG~a^1ukNY{C8I;OvxLZWjvpWar%sy&>HG4$C(b(r=%e+(ZF#E+`Kl`xe|MGWHqj z9{jch9s<9WR$oZN_XFV63qa~mFD+ft&=~Rz&%t26g1-*`9|XaZh7gWyZew%maFrgu z_~2#WmM7}t$B&3i2jFHAH9YHzy8QfWK)HP#A5SiJe&VF2p)uP+ZtyxOiALY|)P-2> zF;cQwmEfTF)j{is(jVk2&n^$-;;ulpLl~&%OoJ@t=FZMV6ogfifjq(iM4CObvm5Iw ze<-j z`{>9)>@yk&*GqK#JUpBi;ntBY#qUu|)6-algc0n5aq77;fDQ&#qt{TWC=>`q027ml zl5%7k$N&;5B&V~&`(2w)IgvrRa(iR%p%PMl}C`S}&!eLj`0<6bgdS z?Ok14fqV=?>mi^o7LSN;!OLz!Wta?YeGwkJ`SW=m03QO>Q*c1t z4A@_*1IQ=DIPFKlw|{g_WC| zm^t9i*eNp%mLr@!Ky%#b%o0Fh^UYRa#T))OCC;5~MW zo?b)X7jQ$;-@`|bu7Gjc?<}c+F6Ui#b^?GMJ#evLa9JDSg-?t@2yB0>9^hC2CA$JB z_7>2ZLRu2~#TbzJyEH66&o8H=&*|v=+%6tL)ZbUV^zeX6?SqU z%F3=^8w6{r!Pn@I0T9U!GU6(lnn4iEJ_TfMgu1SxtBddFz5E??ME=4d6%B~Y;;#JW zcnGqCOqkr<-3}lk?~D%T+UZp};=ua7=*yRH>*F;JTXXV&#m0wd$li!)PcMF=2P}7G z-@b#|Je@5sM?m0_02m-vH_=D^r>-7%R;-jo%)!CoJ{!irvr?d6#~hOW9hRPtpQO20 zIsb`5y@Os8}a!otF!w_f7|$j)HYJ5VrJjkVzO*j=WFwYzbOQ1|dK6X3n0l9Jep zi;EwCqTTiB$$@dG&wC;wC&{D3EDvGKffx;C=gOpwzvKSgex{$QswiObeN)&q&nGD? z>>6ShaARV1(+Or478O-h<8FTFQIP6UY=>H2S;0ct+S;Q0KvlG~wDjhV6~MkX_xCH7 zeagzp5M_@_*ELEoiR1p4@U{>o-IL0JIxkqC2gN<;Es_vV0?I`ksTj-^l&4Dtq1XWm zvPa$n4iFN*{GlypKnyIups$&z_D7gaS>CM3CDl+Y)R8v>K0}9yCu~!H5;DFK}q@Qm=usy_R_;Pvf zSp%Mr97OR5Q1xwYd!#uOeEaH2&9~Sl#HMQ4AM;=J{g z9tcAKc^eq#yQHKrcocXvU-kjE{pVbml@%W-nSi-$gg_8G0nk1!dawwz;+SU`glULu^2YKace$z7F)V z+bk?sL4ot~t5>h^%fUTugrO$Y0LCAf1xng|#uc;rra?|~wNDSw?*7oKsoB~8HNb|< z4Co+spYtCN6o3BA0K&?RqoW^vb=*L?!4~u&05T>88sJDnQK?S8kqQ9iETmI`x(p)y zijY7ZTwG!k5-3&UO`k%8JFX6M3A(PnN4b0bhbUz~)zZ2T zGcN+QxuudBJXc<|Ic6hzb87F-v57{z0lICvfN0;NT#;;Gq99W>G@A9}P}tUQjv{V%eIYVvU9R z-r3m#(!8lmS8D6(kOuc}!MKn?iKi?QQ>bCQw-pYTAYvSXMP=Qjn(z7#`LH)luS)umuYK z8OTXYEiZ?Ff{QU|JXRjZ3xws^KxWL+WP><7i4#VrMEG*OH_>&7cew}_2^dN-F|qf; zZtQ-eg?~UGo)f0y)3=*6kjoF7A1Pu4G%LvYn2Y8o077`CBgcYqE#1BrY97d+yExo~E6 zl@+AWVtRX@0v0e9rX30Bx1Y+&(9+V7_j(5!CMqETBhHQz^gb#L3Ia-l@`U%k9>`uq z!u{K!k2Lr7sZvr>?kr;60$G;pv^BRtZEYo`MmwMHg`!YyJ_hrML?#k6%aB?CRgK`RtSr!)ZR?*d?5Wl6+2J&D8Syw-!wO(sQmTY9Z(BCayx=sylpq ze2xd3&wBni&nzr3P*70ZASH$68W5Ck7#K_vSXo+G5tnJR3A?TL7h>N6M-+D!=AWM* zxZ4c%)YM{Buc+iR`A-g7-n_C0gIejU*>r)OgNFX7|I4R>8S zIy#UC7D0T7r@^UiJbb9~?22(RL7s6!yk*N}mrl973@Su0f0Zw8d z6CZ!jqD(s@dD|@lt*;*rl61ynl@74z8ss=X8fw^ye-RoILPSiAEi7E)+Gc5P-UwP$rbhWv;zYJxLr+ygL^Nnj)Ch5{+#Vtw*P~X9q_|DSP`(!jW@ZNX(Lg;hRK8 z(%PDJY;24ZMAU&&kE*Jw(kwE3>0=FALeQTOq%9H(vI*?}CuW122+lEaqK&JYJdhHF5EV<8>-3s!%v(G6YN*vcU8k zO|ysaI|1kuKI2KoFy)bppZ-z)3yMBcDkQq~Uca-C5cj8OeX7a7RlMo^^0jMk(wUGo zuf2UyQ!f+@B>`(Iu$C{L6IHoA=2Nq?{XkZzdyVoncu;uj{;;x3_BpaZ7DwPjAZjhZYeS*HK+LVIxNnMsu2Z6jdx;(mtgMSF zDmFI9@Ca#}K?hYRKQQQV1kAp#VcIur2_Ze`wDa-*^ywbd(``1kYcLWP6<3{~GoH;c zic1+q(P;CbQAQ846ys7VpxD6cyv$+4-N>x1^J}zY>EoVsDRJ0!0!$U z4{rrIZ&+Eq1?v#)-d4*UcmX+g$Xo!^kI3I<;_SjgIbi6*+7dio7_43cZsRJfRoik4 z>#x9+NGM zRn-*FPo+8vQ)m(B#zC=&ki%2Z#4Cnq+~Z_7fayXFJycU8fmsBR!2xId)~#Ep#Bo=DjLONOFYjNKhio;*1L>~qck;8|bff-Q6 z+&tVKv=ef|$H&LJ7cX869Qgnas7&nY;u73h`}ou(RqrF$uZIIQ06UK{;4c36bG{0| z^3&KN-`l(U%xrUi6RF-EyLJJKsV>l5G28l1!>w=c>=|t-eK!vZa#eN=<>`}wdlHdAh6X~s$^gw04(c0{QS&H xHC?{gqBkNx2mwc%fMbFWd9aQRlGfJzS%25tRj9(KjgbKeJYD@<);T3K0RWaEc-a5| literal 0 HcmV?d00001 diff --git a/DYMAG_solver_version/phate_plot_heat.png b/DYMAG_solver_version/phate_plot_heat.png new file mode 100644 index 0000000000000000000000000000000000000000..13b18197709387d09762e462726dcb23cb2fb648 GIT binary patch literal 39867 zcmeFZ^;gwV`z^XDX-N@~1_|j#X$0w3Y3Y*gE&&M%=>`D-QMx;%r9&DvxC!abP29D8 z-}4=3+&jiS|G@cSD}%9JpU+zBS@W55KF=~jO+^+LiyR9Af#Ax^NoznL$bk?D(h9}{ z@EiW|#WnEPb5|K{S4{^CR}W)nbBMCBtD~)htF5&Ojk~$Ci?xG2FUNBZZZ;Y#S64@u z7o41S|MLY72WLypNwR%fFbJljoQ?|wf@h5QAbk@rvW7r1nB}DA;Ji8WB38kzX)?<-{@wQ<*H7y^^fui+mb|y0Z;HrTm-0(t z(1zff{Ny@bf`A{SFu3a&AjF7YrW*aB{~iz7#P7D-;*sM5N*FdMAzFhNSp) zWWe1?4j7N}X3rtO^y+wBnl8>b%jazScgxN3;85oC1fTuJuu;E#7szfs)c=1cJa=k0yaEK^r)!DxTH#P#7mrHQ;*kwU*XNKJ{%lI<@a-Y|a)(?aem$ zesl4UAFwd-8`{k*52~H%Qqzs0w5geylL%GcJg7_i{q2FKrY5h^{mA?a^7M=hB=@;9 z)?TKCx|1bTTIrzA97ePcoc=^GA&-rXaqYS+fyK!5SoX^F+DJ=GNpS%$_z<4yvva=S zV7dA;6g60~P$})_dWxZnrKRNojMs8xZQiCF?Mq#qdwpM@-EvEs%gFOED?SA~J9fkZ zM{IiCX2alO;#UKz+6|7hxoY9p*WO|$X`^G~9qU%IXu2pGlfoU@fqIyd0bptT2NDySUoD1p`0Ob<})%;ps22J zl+u1<-f?uX?8B+Rk{xj8W#m(TzpoD;7K9Fs%Tg=5x$&<$roO}zbsbwpRbVM=_=C0I zd?0oX4G8$T$|XzbM(ogfMgnQS?BPrm_TT_-Nu(ka4gA{?#t&a2?!d93)I5cOpT`nPZ2Xvm{^0`7e8cG~Z3Qp8nhrHvDCzTgLIH@mY# ze*gY`dwy|xIyOK5V&!Vf)qCuuSYd4lYY>Gfu|?Oec+_elc%;yai$CaddVaE$JFQ59u-Ka-{MVi1J-@~Vm)p^D$E;-&FgsyKKEc36s?knK0^aBwuP z5tcP=KO-k6Z#bLQ^SV8o@iC;XGNsHPXVxCjLW0g035D0vqVI6F8W<8-3+$hV9Qeb@u`ZQ5rchw)xMeO_!CnpLm(P)NMe20 z=x?c-H?w_Gc`YA;WuKxlgoA&pIU)wPy*o3!Pf}v+qP{ys9}&OzrSkrpi)81yk!E?h zV|f1y#N5MzVfbC6B$@nfDU)H8s1(nYuQ&1FHTSDA!<|gZlfGV7%!t zDsQWU_8hMbGTgslL5rX#$!t12Mgeqjq=zqZu7)gK83@OVUw`+2l38ZvS`MXhn>_(9 zI#*>*K9*wOuG4;df6RHGhvWwybG{D;k1PxB{Ib8hySp^AxLCqIxP8&G)Oy}o#*#4b zSLPGR7Llq*CG>9bpei>v*8|Gwe}s+uS3ZGpYkz;=e-HIq^dceW#;iSOKvPw;0Lkz7 z8!g<+w6#0NfEyU~!TEt6hPudCY;BB{wFK2daCEFE3y(m^WLok4LIiA~O4j^=*Nnak zE;vb-fqRLXDrAl-Gj3Nf>id4Q$NFQ!2^U}ozLu56fXP~xFL*E5)#G~Jf}l1_=@U8X zCNkq8Ltv*7=>qip{_fyjMMsAWPyFWaT|N3rbHo(6xN;Cu}$Z-IO`yadkc+Qz@ zu5M^JdUJIGM^ku?PT>e)%JDr8y!C+kT?G2Yik_3Qv$In;brS#pDfGDi)AQ--_8?$1 zT5wtDuMN(6IDO14OrBBP0*k>8-iZj#YKO?=3cY_AZmRNyFT0D#t8-~^K_vdHi}}@`2mf> z|8jNFWmHt(d#gaA4t@b1ejlha-t4jm!4o~jhhB|%XxG{5w67vJ(v*~MO%`3iJ-oc0 zLnaHApSa$~l>eTYueC|zF~@rI=8dr52-xQi#MT{esRpE%`|a6ex=zSvg8fFMA)m7g zmJntfLI_spGTNJ$1Ljv^=asRz*O$=yI}Z}hcOrS^X?TU!53eYo7anO>Lul$*Sy>bT zH%Ey{NlpOEKz^c-%Ed2lJtsH3-PY=3d+9bV^Z4j$G$0#fk2sYguhr1U2;f5(aquX+ zO|#5i1;Dz|0E>YkGT!mpz--+xfP&|LL%w8Xxmdz$r(f0?Yk{0`f7vIl-Q>ywqK!*2 znQb*T;D(+l+mjCL&#b*cK7QVz9E7=>IlK=ha=8+$@7Rfkm=6k7`}X%A4FCXp*&vzR zrVk)L-j!gp#5KTn!q6x#{~~t^b=>9q*_}0{oyopr|@bc7A{|ce%femex#Q+Bf<=Q))PA_NG%8t{{%*klWI-24ADCCRD(Jhs?hH_8h4% zDKT-^J>?M0_+bZo{o;hf{m5Iw;^N}f(@E7$kg_+jefC#tCROqa3oYv3*~k`J=Ou1y zky|V^xrrQhp!m48i5`u|`V*J(&j7|4*yk(`J3l{1iluZNlk6C=tjDIs$}xk030y3? z;DN+U;deS7EarbcXL*9)gkjYFsUWmHVn=TK>0@zTkp20aAG@EjwDSX?7lP91dn~zS zy3bo@f@D#mTN}CH(4N(WC+?dPS|t1PsE>LSkyzUO9%dfR)U)VkbY6$U&@xXUB7=Qe&Yy|xn-8y z^q2C0F-deM5UI`}<$+_{1NLw@K~-$d!FD4CpiWu{QVKXvLyxrd zJI`heM*+g&jR?VTK!i!UxbOmAD-7SQ#zU~3SG9gp;Wu#*wQz!IFh0l-ouQB4HB269 zBBCZvE(Gv%pIDR*z!@yYZb8s*s-Q^pPa~C-mab3hIjny&yrBa}|E4eC2L1jPdN0bz zDtB|UDZV-O5^w!+`z6{vAOfN0ItXU1>idU?!exkCqHg~WU~k8++FiM60 z$v<_J#mA5=|FKZJe>#g_`Fi&rSDSDA;yZ9_|)ym z@bAm&3>wA6)fUp$@6SK_ALBDBr%S(l8H^evqX3ji2!5U%$iC_6qU7(gBz3lyV~@A@ zkJ8-!dh?`3p$7i1@AM~1S1&9df0fW5Pmrhir=E$TzyU-&7?yuD2j^RMjCm9)`JWN{ zFE1~5m$eW^`y{FNwp_n`5sSlA>)Bg>Ad_Rgo5Q-E zY%qLYrC`MhMhyfD@M?HoaN224*1u!gjeU+1OXaTHN&j|m)C7H3J>d5I+#N93pC0Y- zCI3!C6&aZxl~}*KtIfvi?H}+x;j<}CuazM5R}D#0QKJ8rH!9MVXLv@XF=+pM_B@L` zyJeH?JCoM&)FvvR%qqQPX99mQU!+O)LLzfP21gtwK7Kw^-Qa*B z22&Q3to-)2Ll7=%3YE_5zeL5{C%(!fj0k93sa`#g^~Be>_u2Sf{|;NYI=5eYvFH9u z8-$nox8VCeLw9Y~(DZ7H|6hRe;1$UMtg>i7+0CU9%mFjvi)B7f4 z%XZ&lBhNbVZH)hnLtE6^DDVe8JP_FWiItIneF;vXTGF@KAFAx?=BEDkt-uF%VpwQm z5CMqX?ZV-lUo}^nz%s(57ff~#>ihJE(V*l1g)%mw zEsW4`u6~@CD)+d8G!jIr5?9*+Oz7c&8DS~^M`Pem{X?9O#&QzO(m^L&hab6r4L<3t z3qO(vSj8BwKRxPZ;0ry#J+7Jls=)$D)M!P@HP=ZkWI@UWE)p#W)>KH%v{20C!d?dU zwsVGcx>c!{J*#uqN5;6fO2t2sYDHtT1zRy6%mb+r>T0SJ)R&nCPf^lFOfW8 zb|_LVNb@8bVq0OP$`ZQi1qS|eBqS;-+HJxeKVaq{SW_a#@@fGXc0E{uBVT42{-8xx+Lk^qhfa z=r`?<8`-Bn@y8hEZwgxBZ`;qQTpp%Bw`kuPkzAO@}h;nSy2Tj7s6^!tl| zosdQ40W2sN@Ig7OE%&pVT6i9BMi@>Sxo92AY}6R|Y-T+pBaxfxY>#II%~yd%f97eN z+}PygOAimhO6{oE&q85Bmu)!*7N3k~tLh>+vmDtCZ&z?#>n|XOi-v|9-}CQ2<@jy{ z^sZejEH975KUHi1jv^c_;5qO{0b2O>M;0H6j5lUYsp$DS6mdorynMVJyJ3YSq- zt{DdQrv!J$2@acg4l6&XTY_Ts&uIN`sKxu-LExKyo0;wI$+lj(|1HI z99a&%^->z7I1^Bsl9!#6r|DYwmyg1c<9CjgizaWegdgX3GDdZjq5O=vHWtqw))$J7 z^xI%fYG zXsv~03=9o?{H8KJwHd9FZ1~RL&Bna=`ND-^Q9!+i(u&oO@8ACbt7jbmAU#zqo@YyS zwTsd`_p}+#7t%eR$)PQt3A3i3b4Od-L{JxV|Unfe@cN#QH}3Y z5)=Jjzt0DiT}cN%3^Xd~*i#r+?mWTYuy~>;C@1<`|P(ryPSL{iFV#KgSOP{d&l6FxBM%h?!8ZoX-1)<~c~D zJ`qk`aL2)Lfnl=`#zVuHNgL-{d;(frACd?I;3Q?b=GbHm>Yd$QXl*l%t*QYkGGIXULu5x$dQODqc`!lEVPhD`KQINvAVb(C_RgVR+F-dt!%a4 z>m|sLe98eyA*!;i_O@o4keE1#lJ|RKAaHJ%>6U|&vf`kv_;iMdN#R_z#qKom4A%Fr zBmdNo?m>*n0b9IRm}wzp>Ec8+xH*$#dV!b~ik}U9JraKq{kFTZI&jR*fTC~9&$o{| zMcRS1d44`?T3k}{706VeGt&SMv`*CH z?4JD$9dy@!JK+bWIsCr^NBHydHyjn5({AvKNMPlQ1q<=jk*quYRmOs zNx@VXK)$;m1nlX!ED=^{!%^>(=j(}@2-iuoNQDyiW0rRUjqM{HfNgrEWEUcOJ4>Wo zn2;Gy+Hzl?sYHXd8~s4`THx*K1Q&cLLxj<;0mwZyg(E<3)itX8DE|PH^|)Bo?|o6< z6439rcem~L@b<-k`&*H#wFHEs%Aa12%rcMped;ssoyAgW%uYF5%C&(p!FM^dEh_=C_SKE-1;@JAYYJ%^{b7%oWTMo6nM zrWl_({}I+b*Sx;{hO3Q?D}=m7$aWxRg+7!L(FHJcKQ`f9voN#ZN&A3v*WF!Poj0j# z1N7CAXJ&RbRa0}bJOg1XkqTC=t=gM6j}V0g1Wq1#`gXv8;npm2hc14JgwUr5!dcPv z@JYej7!kuKUCSUvKr8Tlm<7l49iBhLypztIeHK*|BmBygq<90%F!9jDY~_f8w0IJCq$`0eN+)U^d56ny;r<5N@8Y$Q**rVl9b z_;>$9V`r89s{OM!jY)eOn5!%68LS=!&vDc2vI_1Rd!YPuwDwFX=bwVnn*+BE;nK3b zYe_|bESqjkI4BU5+axJg!0*$rwtbhWv9}a4~{{2lMMC zz{PGF#-GmiAR5Pcv?XQ*nRnq7fFkze0VJB0{Xp@Qpq_Ho%2)D)9?C7IQ$v>a9 zQrSMae5Xmm$>KJI6wh6dW|12|@dCC#$$VoRg*2t{pb!DNs0&&$hkn$|F*{9o`o@&CfEk~WcXUKoCl1oBVOO#~~yv*h|OD_4?h zEHSi~iH*23rj@2ejI;`*t`c+Oynoq_$|VOh!h0$2QcgnXb6AjI)bDp;h_^j{LPGSa zY_L8RTNs;m1)+37an!0^qP(&)jeH9+lvuVVo4STvh{1Ad>A=Wb;m0#0?>e6POjkZQeo0~}M|V1$gvK=u5a8`iKE zulNm=X#r)SoV39=`3jaWQWrY!AAh!^>i^KBMC5)(t6|kt{O5(QK-v6HyzB-GrQibA`425+*X;(hl*Q=KEh<)czri1AbyNJ_%l{+;dN>+z=dy?%CrTo z306mjRdDou=mV`MTK21~OopmjHW#Ta{ysj!v0~TEIq=2J>;KTs z+p2)Z3+XcH#7!^HY*|p)KvZafN&N?S;!hIl@WP6IR|x)FqevAn0o@%;?iIc8qi+s} zDT8HK+0ZAbF*1yf!|D81r8*Txo%Ym|7yjeGvOL|I1Y#w%aA2<Jf4Tk zgc(QynWuu|fKJIBGD%alULO8H`Ou!(CTUc{EB1YFb`CXm>5Z=LuOH=Y|Efm+YJX?H zwN5}ziB*Dy^wSSOBJ*f#LA;7*H(lghX17oF`HfoMMk<#W)Wuw6RX9$~b2}`HBx)c8 zLwSrKPA~`!lXOOeg`!tyyXpglQ!#-$Wxs%Hh_&4co+&>bJ8KZJ=oRj#4=EXc;|9za zE{RMs8~#iyrn#-S;<#*5^X6P1V9WAfqJR|e2h^=r5tu6dM) zFS4js0TTW%jFOQ|nD>{8QY}bQ_do3BAJT3Z{UW(ibO)CbJ733(mDlq>ug+QZqCym; z5nTs5t={|xEhi7f;d7(Fc@?XegtY-*|Kt)Dd$R0vPzl&6!gV=8G*uw1Qqu{0+E)p# zNS{8e!^P5+UcHB7a;vO==qS(h72VOY-QFXl8eG9MLLeXz8omcUP#AE%d4P+Cut@u9 z!yG{=a{si!Q)DDb>)ilwzc5gdhjPSofTf#(&Cy^qG~ZV5E7}r4ad-GpB@+1yur>R! zxO<%k8|9}RJe`qhjOAnlN%a{Ma*K;2K&=^^zw>JmQA+MM=i@WJmPX^D;P$vTl9xBp0grWF573BHP%&4F&78;c|AzuVh+CAPNM9-k^{$t2X( ze4%}A1|p*ne~H0O zb?IN{t1O8Yh~vU;M24b$3K%0Sa_XE&DoQX^;0{j$&xJ5N8B&43)Lt}z`yr&|zR;AT zfV{sz$semU`G(@dbZ9LgexgGhMrmmJhylEk2CwB)%JEHFjz!8>_EmW>rEz5ZyXZKb zb=s8}=>hVF6M$N{UnF2c_0YD$LPjxAIpPH;fG!RL(Xjza>|DUVS^KF&>6OXE74_dfA+GXG>;X&hgTi?Qmd_oim$_u zOn~E!s7Jki{n`;0Em&nWP6v2jWgivn`;Rgc-~ha&3Yb(DbOPuYqy8yOwJGqQmtk56 zR(xVz*#EL^)WN|&ar3?|%1nze&@xZ*oydEtoaXbc5QPt3_PNXt?DXcnyasqwTFs79 z+Yib?*G8g-$oBq?9hJk`^;OOqXnD$au4g;%KzQ5 zyai}J0E&$R7od3effXfT2bl_?{!s%>rVv(1iV~yd=jy+V$ff-3I~hl|*zGmMocP_i zp=0LizQhK@23PiZ%*YrloILZUeorzHuftZKZM83U-`D4EN^>5B zLEkZ6Q2fN54aq|)94dzS_3fn`Y(yX@58io|LGO|H>r9aX{{In|>n}HG1_DNAC5Y~g zk*vWx_<^d73|c`!0fIP682P5i#Y5C^WI%9-14e|9q4ki^U-YOu&%WWI4|Gp9v@z-v znvD1-q>1D+B9_&*>(_9w_!nc?clOJ+UOw4-F>ar8+jss~=7_rC(K=^9T5s>G_ck`q z*FL{Q)TYBk+VE38OyZO$4$dl15?mqhC;6W))tG7wWCqn6{doD5?Q5*peqDf1if{%IZ=q6mO z(ApUkO`|1tLP)7WT>IGdMNkpmJj#|9(Mxf<%agaAo02JEk`c1yb$g5J-a(f(gpc#f zPHMMRVr;~N=wm59Z(!9wcWkcb-GKS+*A_+NYPR6l^TwK{nM%DqZ;V(DgN(*{>#%;D zr{U*!p0s=@MG?w`$vCl=l=IOYe62Hd;${yX$tE5i7L}%x{Qg71kuo;vK!5;4CsEp* zOo$eX8-=-B4Iq?z#?(2hm|lhthAgsD&O3Pz4=m^m=`R{`TmqcO{xezlT}e+i^{tS; znEl5-ZGKvxnUUzh!jTsB=@HkrlCln*-a`@PH|D5lLfPVLM^ePi`bU<4ruMu-mZ((J zwGkQJ$y#UaG<2#&VbmjUb`{rAf7BVAe}Cs{UxOJmu|cO<9!O5F~C^w8cU$rFEa;C z<`dINMm9E4#ct-W69=Br&?&Nnnkx^*egcZ}g|IMYs1(|x*#QQ8ZsMzu$0+aNECmdg z8Ibs5oJzBgjzSP}+pZ1)ocR|p=E*Nrbnyp_SaBXej><7;sC81?QiI4*Ca#aTN>Vum zQkh6gPnKfzvUc3zx(nU>LeB#e$R8mm7nbqK9pyd;I8a3;m#YkYRjnX0}%w0 z_zEf`SXhAq;t1D|F`n_+(hmw~2-QMH#f}OFJ|{4%u_gD0R769d*FyX7=4zo;$Skr zxTq{I57A^XIl}pf#E6B%%pS)d7#g@s+95Rr*+!^$i1rNwB(vjd_h9$*@}r_Wwy0;H z#^xWPygSB}VtVR>C%F8}HwpthX)r750;a0}<7#_Oby#mt!Fv(w6!bOVV^Rlb6%;ZC$v zv6{NNo`ZB7Z;{FVYd5C`nhfw|iOQ1;U669=6P?}hMe1GfMmZ33~{=x*381_=A4E=>3GGp>D{0LMGMO@G~yQkQRscws4 zTdyX21$66N&F@#c48(mI+Hs$m>@8vbq|KO{{Eqp6Jw8U!%Q_rw(&lp|RLab{6x*FJ ze7XyYu(lnMpE#I8+P${2TJKAP(D%JL0-{-*3v*|umWR_f1iTk<%5o!oILJm5&H_UJ z(#OdWpOQ-GU`p?T`Zv*{CnmDXT)VyvabPmrij{mL$I<|CagCxIpVS-cRBfJ@FBO?I zzhiSQ&~{s~Oi7J_Q=?A&0&!Nc5;? z)d3QcF-WJiQn>V=n5uz0>58Q42+PG5T26yV{E!q_xe7~Jul0isER?qhmh`#2(6pZS zCJaF#Nh+Y-X3b-u2rJwtO?s+<{0gH#<{kyB23wmkDBJ*4%ln(jPMC-y+z1y)F!rq` z@=$skWMEp1t2jGOfg_&r4X_pR^jfa_yy@hCv3Rx;94aKS+L(8>35dKdO0^_Mr3g zQ*)z`Sw@mi1QF7`GO7dFzbuYERulGP2)v3?1e@Swyw~*bwHFJP=t}R&1onZGmw!6{ z{$)0p9g_D183mNdR_}F+^kMBQ7OUM|+V!T_Ng=PkH%bxLc#6EMAd)a+5#lDsgfKN6 z9S@(=-#as)Q*(cU-d&R|d-|bRMIQw7tkj%59QHVgROmhX#*8^J`^Z$$g}06eLzj_u zGOw@%wd-?kr&Px|As1|2jHO#HUbvO?vF!o_ex27zbhHsO7fwy?81vovjG$am7}bIR z0cm05U>|*RTx+#*aExnu)GrGpKR$VA7cbw*_74o*W`A5-#z?7m4^;lu{6R)+IO=OQ z>^P}fit}p@5*Gd=a@X%XBC@%1!*5S9Th$nx^dnjYhgxqT zU7V}I2c%E&R)JrPd1_k{d&`7^Xq!?)Q{}MS$b^h8ZyP4i{|PAAuX}jFpd15l>c=41 zl~u*7Ys{J( znLWj?v;htY?4$8B6(M`EB~kfgcBr(L7Qka&oA=wlzI|sK0!Hk_6p^V7+lcXG;pcOa z7ck74X;gM>EqF~Y(Ui>kCL5|*A7pw4PqYLy4uR$#QbFXp9`hXGI6*DjjWjKCu9qa_ zfNLf1qb{S~E0E*NHT&wb$O!SA#=n_!xJVDQR?{P=HA-cmNEvyp&~}DBE=wDWl&AyR zTnHg5j<8~E(yE{V(TmpbSEi&70LOYDGP;R((LSUl7Vy3;BUSw4Pwai)sa|uf%|JVm^g_0WXU#@5Ew>m*!k@9b9wmb;;@RhjO5L`{7N$L5CH8i<2)QghZ4K2-2il zTz0kqOKdAg18P@yYhCEf4eaW7OuyI!~Mmk$G#7 z{UHnUSeUnt>Jze^H&fC&B$GbK#&sFVc96 zVyyUp93-hRkB&C;jU}=e{rWVUCR>GwPgx^daWU@q+mOGk{}N^@QT|8axR)P;iN&v) z#5<$wIaaDB@bBF_;l*W~PlFlN2qpBeMqiiM*%k>W|MDO@VI*BVBkjngLZ6+OeA!>c zj2)V)Y0!g|u#=`~AU|7YEzj`kd$AO8sMI?wmFW*eh)PwI@IYeMz39yg?*+>KvnrX` zNwO3;mA`0?@*t;X?oSen?w^ql>Kq8L=}px2gk1~s=!ofHUbp9fz_kbLSloF5Q?K6B>o^=ml|PghWnBJ3fv)%wDJp`s(hO7$(5~iSc)Zq4Kbt!0Yf2n_g6) zGv0Bq132smZ)XSjfv1Bh>uB3@-~*nRPY0Im@r;uX1H5;mZv;5-bQesm(A;$bjl3#6 zsuk>p2tPUaUcFRGdtEu`a{6`h#=l38&()o)pd?Yo%`L<(qi;s=@U6|@tbRB?8Xd`C za^|xs?!f=nv`ACl`B@S%b0yNjY^hI*_-{_mE{Af=cC;Mq4svhE~>ZsEnNXQ(?@_975eL9Z&4Jes>JNsNB<+W`hQ|1Bns}{o z3LX&q#ApkY%7Z*FZ3u|98gq#GqwT)0CR$yu+z&hf0E2^P{BJbU^V37fz^k znwbb*=}PGO6Yc6^eyB2qrb4m7a%^x=Zr}Zzu*Szu-AUNjG+kYRT>&h0C6h|h zI5{)@W$I8;Ohg3k6l+aa6cjuGDs@1~!Th=g=6sM}zFzNSszeqkC-ZsYyJI`MJ^_Q; zdrdu2SHa(R{Ix^HC1%P#PcnR(wiF7qXMW4@;ji@%X!fo=3t^@}YkkOt=-|=F-Kgt9 zc@5F-p( zs7#hEt&l#U_v97{HW#VxVwq{yOLp_Igr?=mQac2%#9&||8dr&x)%2|lnn4;L3$)<* zBkH*apIRDMT|8h{={5m(5kGPDdm}!8v>d*qc9WBBA_uz)&TD1r3?6NbHYO#Npv&6B z@-MHt314^f-7*!I-8oA5>9v>U;kDu5aJF!9*veUD|3))gm#^k0 zCvxe>_J}-2kQrv~&L^JzDaL7=c4eCMTs}Vw3TSuL@qZ%t>iv&A;G-z8O1*w9<6S=y ziJvBTkW5a7k7q{^r)`;0$5h33_qIrTLQd-PK4{c=s6 z+ahYvV75VzfgtQcHv<=;xlUnmmbqGF^o>it=2 zxFltMFGUbvrZ`p={=hv;k{C9Fc$HG+C6x2c4;|EZP|K6#bq~sv-<&0#8tK_FsIiN1 zWL5D)@+Zr7KdY3tixM@}dBw2VANztH2U?{lRkcb*SR=M4)bIV}lyw-0wWzNLgwYS2 zlJFQ?YMof4_h7s}TpHvHO^tW^O!8kToBpuO`^2~EX^F_tW#ty|3aoi;|0IZ?UtLcg zW^;i4ij1+yw$>d@FMIp2G($;;M(BYS{i}I%OB|}3bz9wqU%#p0SEO~{+C5Snda=-T zTsw$A)a6yzhFvT#e>ku?PkYx>)b3@mJTH(aH46C!v;*LsJq;;NqWib9?pK0+Cr5Q> z%3mpC4MmTf#WDo83*QVqjtF0}SiKT=>UtQcJ&1aQLGCawz9Q3A-{#HNy z$v4dHFth*c2BHS~CJjL8F3N|qd|BE>3kN#Wt+BQ2P)P2%d&bW&DI*pniJM(8yodSV zBh|o>!imMtp>e#tu6QSo8y z-a5>zo`;Ex7uvZ{#6EJcFI;UU*=9dN4`8_OBJ zBTo}b0eL2RW5;Y~XO{_@=MdLhegD3s&a12>2vp>3O+*x_^L;~28W4>59UsNoY-qK1 zWdYfOz&N9)LcVgpRgUoeyTuAI$S z?U%j%L^HGm2rItS&~+17QYzgJy&$3F3RN3hnVoz*!0>sT=K;{j;NB!?&we6c`&Wzi z=mR!~(LnrDdQWZ>J|Jn;@!#oI$7R1Jr&i`k+B<&pOqh6oVBO~GIxLpSC>nYs?en96 zRfR2$>Xf_S3j^Qj8gAc(_7exDQUhOEB0lyPlk)Q#PAZuqVYkn0UNy~`5ZBh79?n-q z86=CVBwp{O zb-c#PQgG_uLO(n(eBq`9UC7G9T=q1j-`jI~zoSrG5;-=WB>g7I=`R`jNhj}eU*=Up zc%tm=W(OIC^D&a?bX{wl90M~uWzO8dk%{Gx0&t;+2uHy@=tyP%CXXh>_LIyciU--+Mba%+pOkM7f@4)oe4nxe%x=H#~8&{wx>&KrAzX~8j> z(Ty2hkB^@ML3OU-XgY@G3o1v)pl!MC0;<5`3*AICnSefIs?JXY`I1e6XrKLHF&Ali zoNe`S`)1^&x%Al#FmL9nw2JoD+H}rLXTEbG4%W(grOA#erCl+DMctqx{@O<`W7WT7 zg#2|%Xql~hDCx-CnWW+#;62#BL5eY9mDTVH3*1Bf0;)`wIezY>oE*JBFYmAzrKx2{ zvqs!DGrK)TMANp;l2E%yem{B&{tWso-^)Itlw8DLcChNXAUJd3y+ix3pig;xozZHP zVQgXhzBHn?;`v->=Z_Ue82@}Wph&hjzC{pIF+U&Y+gy}4F>{BL#% zEY&uby>6C7ICG{G9Hsm_0yWhett2%cpI6%)`6=u9;FT5WRNNCk6mq2A7#f^Tr#R4Q zZ^QD(Ss~S~u_6HXDLD`qf-Vok#Fa#-Ho>Jp#9aeRPz3^JH|Ym6h)gF2x@}CFVAkdn zfeF%`+!hAu4jyY&CLnb#cl3TV99A5OmRYIRo_BA_7orc3JS(u;i{SOq_^lFqQ8PMG zSgRD781HJ(#KYoSib+O7uwbImfZD}1rVI00K(PY-!h&%D!^2A=C%fCr?I?CnE-y*# z&!GzDgLG+6gSfHpAVA$EPiHT?&m~JL`easpCqza9_Veq5)DAfY3Hchw-q0SH(*R`u zc(n4m5UsHwWcyaqci~qvec1+uL(6u~+sced zTklzWmz$I+t@I?l;RmgT_*Tq$wFGUKp$UB@%DOhdIX~H$s`R$wJt$jN(9@$RFq*G< z#yFEg5DT!ZFwM(l=iB#ZnYQWwy5o>JQFyPb1ChW*9mEZx;o)HvZ+vi>dbc&ek9!*= z4738#Y%>ENCfWtSn!0)B2w4;)_2@cwTLV+q zYumy59d1o@l3*)WV{D5}^Jr|HftjB^5lsrx8f0K^U8YZyy1*cuDn8o2*r8l667xe+ z=eytz`W(#|FaOg5$u^%Z*;FyR&AFAhLss9RyIXh7W(;*4PHjv`b*G_;p)`P@b+g;T-25VcTD7&HUNlpHk9q z6*&3fDUv0U3e;6ujhq_kMc)2GH$`^8$m2D;OWDl45*49R;g8OS7%h+naHWTsn% ze1#p(Zr5rj%L^Ih*B0`ptp&XEOd{GCVvC&Xh|Ut#zhBX4B4lYw+`kakkkIk`i}!oo zAeymklg}G zRd6G@QRD8HohI<^AIq{CW^Vox0Fhb&NB3yUJ!C_ZOi<)}kzTFxBx67|BcpkCb}GV} zd{?lDYV>4oo}}wY#A5XO>bIAcl?RJ25;eZ^+Dk^UQYmpj*uHZ1ULpGUju6*Jo5@T# z%(lw$0Aqib~i_? zO05T7VKzAx1PoRd(EYOWEoZ{lQ4$4rS?1_)m1hI*!ju^1h|;14>-{eUr|qq}E|5^z z79@oIK3KoxcClv%OkyRqh9z09I=3?a0vh*4*+n)}j<3Sb>H}TekeUvGac5fsa|9~Z z_gzIp4P=XpMem)`%z(Qqk&_`_3datlS@(PRIJktKRFoXC1X@}kSU<2ySAH=02aZzu z2|;4Zbq4J?+myVxpN-OXz5YjmwWfcm1I~jA$`v=Au|6DO=lfe<@5Y~x;%uG zbSPcYoq`}K-Q6A1-Q8W%0@B?rC7seK0#edl&+_^Ho~r~9gCD%8Y+n@IG(+1h_g;9?q^l(f5$u0#)i0=9jO|gvY(2Dm`iTX z-V@3w7()%_LPohdPfMu@n1U0`|yI9(a4yi z1Z_T3#KO^0VNI*2o4}K05uDd=BlAvhPla00hIZL<~ zHp{j6J%?L7E#2?me4;HhcQvHwz*hR8LL(Y~GELFWPL7j1iWWa5=ff1$Gcu~4!m&9r zTJO}?vb#BNzjAo#LEt2ucehr8gTX#PaRZws&_FQ4{x_!h_Izf<-BI?PX76E%igHrO zLVZJ>tKlaWIA3`N{O;EaAt#`7DB-1@_uaY+3l+#np)9em71jB!Udg0v-SRj56G*5W zNpV9D-KnwZRngM?U||vFXK_VuwdF>FKT*)geh>Z2^{Y-LJx=!uE{A%NGKYuJH zx;<4~T3a3!F%aobef5aZ=FS?Jv*l!WnVf!YBxR$~BX_r;@IOn5+I4@}Kor_|seuPg z$3)D?e zjTiN=zY&x|ec#q-c_kh@C^W2>u7<>>D8CZv{lO{2e^kvg^D}9~cqc`M(+ zP@kZdzS{Crl&2vHT2hj;`+NfhDKND0o(Hdf}qv*+w1&I>@<#ao&xvq zAR68ceis*3fx_rl1-AVjQejzlWG%ieCu&dyc7?#%;I!#3IA!nWZ8;o9%~o!{F#JQ~ zQl=E&i^_wXJ_NGyuX3jbXd~4E?Sc+p=y=y)*JA8~e7w#IfvYuw#)}Ji95k#N8={G= z>}@+-YnG+Fo2N!7hqw*&T$V-grg@Rc3nT)nMgvS{_SehQdxI9Pt>=^bn~NH+Epu}t z^wdE2;zXF>;v9-K-sl@k>;DZ4@1Ne?rBifJ(wt(99lLTgy5Qgo<Be4t#-X^7Lc{o4)GD$o@ z?3I1=`fBng3KStIqC!cs?d|Oj-%#*puN!o{9VnsCJ!$hRPEJ1G^xnOGy6>n&)3sia znG9)k3b-;(D7EK~Wn5vC6=12Wl*(7TFR;X|f~*Edz^mz#Kf-dcl-ULo5q z$mOVHkzeFr;nEqUop+E)e5gh0(2kS$PwtzWt0UoBqjjMw%HZ^Z^(z=?>yMHa_OojR zmAf;EQc4<_rXDu%P!m-vpQK#_a!tPcd9Yyksyi5QIR>n}eSWs*-t;p;a%0q+{z-U# zO(Q5!4lGyM_2lt%F7%uTx_upw^@RWEJaxgmVu3*`FWcRJ7rXQF6VV$(QK&QgPGm<# zEf43x`DoaS%OUogc52cVYCD>d?%1&fEyjK}J1dx=AJ}gFjJf}MS z^(uxmnX@^2qV`f2l={-!2rQm@S!1uMId8Ib8${i{*sjC4zR%EUfR0i(%q@oI8Ren z6kg?(GT8B!M!W`1p-8LmF2@#<8x~@OP)CJIy*nN48TZP(+tvXE~r4 zA(2{@$4XRG#QuSuDA&LW5o*SYe!eM%E6b}(G*?XB(bNz`ve`~DLia&Y4M zNVVm@(qsuQCDP6A;GGk+6?dlX(|8Eq={Y)dYx~=FO4BGEL8(l-2DNT|_z-`71Z_C> z-yQ6*JQ_-xqUqgktLSAV3aqgNi+@yyA3X`O_I49-_&|O#-7adh zHCkZSy1Hf3*9j1wBgKMSuT~&nX@#1YC94`QbT(+aq1D2HJg{8P$l-e7uelopheA zkJONP5neAB4N9Khgfx6*Gwbr=axC^x>sY`wH}*nz!P`&WOH5u$jsngS(_EW;TRt4lMqF}a#}G#T~}zt1P(?;pfD7g+eV%iaOG^*(+I_hNgP>TfX2pf?Ph6^yV0kNIVrABz)W6p-n zTY}Fb%;vv~1h5x9VzWzbr+G;6k~bz!`5T?C8Yz-NA z+2}C{N5W}^LYa4(#X@S)G9;{s0Lmy{gu99>D{b-Klby*(48c*||tVzKXg`40x z@mXje&7{bj9VAzQ63E93V!OGv^2|JD7!(}+Zd?>C?YbLoD%4ZxJvoNU-4T* zD9)%UelpCv?+dd$yk+$ZMMI#p0i%|YH<9+~__9eMB%cC%qdwR%3;~vVX;4052&WK> z=?2!rtO-{X>U&7d5gk<0sc*l*L!pS+8um=)9jVsbOp4UT4N0a)tGP7wiC1(L2Y4Sn zW^YKQ?7UaK_V{3N8e*nems9*$g;eo?itF#>A>XA0w9u99Sm6FWkkAw;qC z#0)w%{{txH)E3uXS-Vz=j#zJ!BOcxsR0=Dcwwjx_Hg2X8E}V+&0w>-^~vYTg&<= z$KAsT3DaYV%*0#@eU`fRfQ$qT$@R@Q^m;!|UEjwMpn4C+M$ge<*ZnvTPeT*^zq(CCK)|4L(vXyw_Il4^w#3E&#(EHK;$yG5~t<;+^Bl#2?>VWNR;SD20ps*>wV#Xdetc* z_EH@MZiw>Lde%=6PAv~i*DTkiG|0r3&KEEzVk(w0!cM!&bcdxg?@XUSr{_DyP3JSi zc{>*6gldh53K(_q5EMc*TA@4|-}$ZNxr|WNc4?HRrtAX6QZA%ff~1lKJ55GCp-Xet z^ijW>z~1TM=?qrWqb(ebY9H0ljOdE(cbNt8@#a}pJ9X*W(8S}$#6fXM~k0la#WpZE}PN?lFwb$Z*6?J z_N+SFb3Z1%bUC!H9USi?3QaL4a*cz;nc(%g(4f{8q+D2H6LUYQvadBR`I2K7)*H{E z9UJ!9&S;(&bFERuI|%z*;NJ`*QWL`EcUNQkSAUeG39~Wlih>j)IWIZ8JIvr}uMf^= zJX%?~IbCY5nXg1BMm$l>~# zCHQS46P_}i5n+%XcgoSFanHIdC*)VU@ptva)vR_fCf^_G(nLW<^fhy%naLOYZgc(l zZJvZRacqZf@z+r7g2<3*xTH~ABKp`WpV$}~l!{t>9H1&|yglhbPmQ~AK=S{_l0t_>`40zvIiTr*IvhVIzY#!ppg7CHrHnnsy4@ zX;seQ(LAm1G%euS5nowdWcBcTWqof`!mWiVQ2x)-=%c6Oabs&)G-Ei)lHzF=i4{#G z4`ne7!Qe4@Za~^_-!Ej#fO`cK>DPa!>O1Qa?f7r8$4oBNV~*f84f?Eaun12>xqdIM#J@F} zErU5?Z0{8-NvS|4Bf$`UY-aTo+>5IuuHS$vN*nW)H5uerJ}dDLVxJ!zAT8jE^ycytk+`I%8{)Fpx{Jf<#t0%1Cb(+S?t8BYwh5tsOvH`iG||!q%f0fEfl8hRtoLV{gaI? z9M7*D?Te?fld?fc((<(YN0^PK=bdMf%)chBkRmr16aI=U3QNU*5y<$>LM2hxqbtE6 zNn`CWl!tga^?1E}Fdt^Z8pQAlBx<9R>>n=h`RjLdeyF-oo%(PzwEGUv`uK*{l{S8U z+S4)xN*s2F8~SXK)Cq=rFexu`W?09oj4q{VgHSu)yi+suY44L|O{JOxaxL;tOVGqV z-COuH{Bz*%Fk3ZjL06nO7k8=+wIKkv#bGwrw&cq{*7_MXoNLw%|2nJF1?Tb0Y-2G0 zHg2CcmyAB-mSCC0I#}D1crOrZKp<1W(9a4{tduY+N0gDZu<(;c69|(==W(lgmIJwb z5{;&iLPA21!NAHr>)o+AQ_lX{ZQW_mcjEeloUN751a4x+85$X^xhhe z-Ch}+#0qSi=hPou9URKGF^x09%%+_aAqZmVdtTRsiX>c(3Z-KU8uTqu-+3CYdOn== z?>g}3?0;%@MG`iWK&PGWy;OQwzDFZcJwc3fBAU(~NlEJk6Wn2gmnRdKA z-`V9Z^6(wiiAo-E&fFeQh;}-e&z7s^S&OD>GggnOs!M)cnn5l~O*x*VBK-(MT0P{l z+oJww|M}MEKCI#MMN><0W=`kB+lnE$5XGc>5$hQW0u<(!+N6PTk|gM{y_@Q1xqa#U zul(MepRL_u$7)Da8yoT(Oci8_QyF-}u9RzmXO@23^G5(;8ln``iia%?Mf0Lq_hRZp zr#jiT8geVriJ$5x(YE_NG(+M`;1zlu{hM|ye-S*m99P5(6S0k{ETv)I6{o(N;Eko4j8bZDklUas9g}F=W?z@E+AP&LAs5tEI%Lnak7mAfb893fA zv+uP~tU!L{`6mh!`3k&F{xKuyb8KxOt6s=LlgH$A_r2aZ-_bjVZ_O=uP9_0Fv2FvP zs^2=onXTblQl(Qyn9X_eHEFNIjOpR&E#C|+iRLp?VBUt#NgM|0Oy$Cnm@FUZfCc%3X^v-cgt@XtQs@ z#eC0}be}Z#v0!8?z!{BFW#-b9b!?4J79hMMW}Bn;Tqu%U*cU>ku+~^8;`w{tHY43= z@~)RA%4Vd?q%LuLPB5rqv0!N}{O6^&t^4E0e_D^6EL^>R&J>4x$*~vlmk(O`YGeYd z+Mn@aayoZG3N_qD2?w%yBj~I3%!J;4jn?IR&FIDMNMFkBJx3%XT0@DuwHJLbaPev- zMtXL!v15dCA4R0b33^1){(m-NBVC%zN*|!OVTrEhsU9NP2boH;Jn&xoRDM`fti4Ws z)IuRp>?#c!W7Z4aCF0$9RU!N3@&x-S1ezon;Nz&5LZl^Yec2x-GONC9*{qgkw* z7?o!(nFpQOb7me|%39ELlP6wzHN-~c9O*WoCMp$Am{>~E(!S5sfTU3&s)7>FkD)A- zngc^KW=8l3Gg2dGoSEw{_qqEzW!zkK3(Qe#hx$~B`12pXH_Uy+g4!8?LNmkamlee6 z3v#*oW5nI?FHD5=Rl2`I&S4E3aXMTkrcMVNpCh=b*E}z1dL9}J1>5=LRgF|*3pz^h zoIZg`TE~PPrT*|u*Hm?WL2z2)P-Ta;=C^ih3~olvkP>W*^#TT>l`Z4p9XGC6YmZc= zqHzpYe@9 z6;(i|Ib`A?!_t=fPHY&)Q`OJYeNYu#rHB%$)sW8(>CJcb>djYCY`4^#C8&q=6cEl& zCWT#+u;~VL_55j#xHmj^uN%^UcvT0JR)q5-KIzDOfhI4X=n8dqyGprACLNF^M4QM*;f8M7PtUh+7 zkIN)P2&W|$kCShCf@s8kzO8Cza@8LprZjV$@*I7qCfaI&=y)^!vaH2^4f;_v690I> zF+Cleu?_HodN*0zI}|7rh5InT{Zb#=o!pf=iDzbrF|6;x85fC#kEsme!g5i#h2N6F z%<$l8he#~)G9nx8nMUc6_jl&$q!PHIF0_gnqq%?K8%Hswnu^H$6)dQM^_!UW5bX*F ze$o-;+0rYclWQ{ays>Q};ziX9-o#!@m^=jjW}{L9tZ*`wFC+_T?D zfSLW-u_)1GmEO&F`#udnotMhS!7DcU6zu}7phloV1RNm81>>Lkb7#_v_hNkUoe+|u zceyY>uSOhnNLvg=V?AB-VXC?`o@fr5erOa~8bWy2Xpk^KsjmStx);M=D)|gq*re%Q z?%Nql@HDMR?7ZGY?<`gvSV#^k2!CW2Gok4la8m=7HEgE;_;}x>N!&8AVMX>wP(dYq znEbEj0_F}?Zer|0XR0p1mP>rp`J`Snb7ovst^_AGB0ZtZz_Oy>7CAW;57ZOZI-xM*8I7c&b29zs&fxoo z8n1)u-v9| zE|+8ZDO9IFTeE4)9TjhwD!(;s`U7%BQe5Axh4%Fri-`zCsSg3ef7i~-c#Pj`mqHTtIAY8ufX;XS_ku!|xr5oWMa&{Tn7xPhVI zm(t)D`flgj-f03&VR#EZfyYw1ie`U(tzuHZj#$9{r!eUq(jTMd{kB<-x)zL8&er*q z#9RNIt4?&wv(hs2N70x$l!Gm3fO-G+8#f&_t=gia4qZC&A_K`3ZGNKl!D|l9dxqE!zI~1E40+^d!DHrt z#?ReXNDDqWW4(?1QN+JL7rqS|(^%a2OfT#|&WMay)a0VPtUYQdmyxOZ5a`<4RI8`L zQRw2AZ6l(;bd7{4LzS9*lKzOk+n)Z)^>#b=7$aa`ANh@p5S|u(`{LryJ2$@Eo7ff! z$Ni*RDE_q@CVJi`Lpsrbt}MKAqy~!&&#%Oz(a;PiFI^pEK<6k9WP)56{t>bK@XFmc zr$Q{kDjWO!lSVA%|gQ=rOnhS z1rJ@mx%Y#n3%XK^vGj~o8I^QAB8;sw9+85;Z~J4sx`XpCYt@wWAW1#+{2AAyH{g2v zIFRaA+*ae7Q=p?d(1gK`iXD+YC52YD6pGatjSSI(rV5VKw^?E_#6x(gYr@)9rC7|W z3GA;$3Bza;ZT0(K2&J|7ZeE3bVVESvfThZqsbt?aX(^&mW+|(}ZD~-M)1fL<&Wy4m z-FX&40whQM3=F?pr&f(zF;!$aL&MS-u|{Cf5ZXh-7JkUxqEX?qv|%(|cR!%bwy+i6 zn#qVUBBVU~wr*LkibzM1FY9>PJxKGfvT+h{moJZ=#>S||O`z{nk^zb!J3IUL2o&0` zo5JkNLKR={Mi0@QV>AxZz*x<3F{4NB3Jz0DXWO%ZWAbEdh@w~VR!US)MBj9t;~wewc;63D}k6hTQZ$^r|gwn z;R@RMAZpc2Ld>c1x#f6Tt$aY&>-sWQw$9^opewJBvsm% za~D}}R%4CIdDLYWM&)fYlGLH+C{?T5ccJ;S@RK>)I$xo{gorHEdNlbN z_$^)Cpg|l3O>Yize+WbVLc@!1ctP$D6)enZG*?$g)z;nOY~RcEaoT&6!sAIzr5&1z z<_F4A4@oMSp7%e-_AU)kz*caeb94;fiP2x%xMyQo?ikT6e6{RXV}i!&Y}{Zsm9PSr z1>+1Qzx!IIW{Pn57>o-zdP0YC0+(y zM&XcJw5sB|)vF@b+mk=NV3C+vNI9PsaG;;+w!{0;wM7l}?E6UHs~}2q`?fNya;bu9 zVy*Q%+!x7pUm&WN&h;6_j)esr+JG$A=z2CtBAckFd-%z8qqD-?2}NW zj4s)g>FV+H;U|(33Em=ljGuu4O;{Tr(^kr(3!>V&ZVwdb1Jv+Y)5${MH2n4y>4(vGptWbJMMbqEAS@940%%j+hmX(X>C2ydX}r zsZPGL^?R3q8BIqZ7IX&G?;{3DeZM|lEx%uvCD8cKsRC&yr!Q>Yj{W+fSP$E9D#$aU zp)hf$Th?b_q3S0QF9k<){^{v-%e=N2Yx2USEFJTl4;#is+DUR-J#68{&g^# zMW7cyyrgOk%}6;FJ6|?u2YW!!m!I5J4ikBO^>tW2GO7KEJ71nI0udf`Dq^Y%mMXsq zs`qOo#_UwdPMm$&6#dX=LpWx`L%V)0wHCimBpi#i26ESv^u6~idHDZ>GQ}0 z$=@EB<_4a^jw=l`7vCsY>Njm<*P=#-;r9(v?_7}{7h~7AApf)SN}3qI!P{<}4JWKj zUd?o`K6OtzV`dg9vsrBe9XW$z8=zCFJ)3?vi<#`k>XWsh#ep>37?jyBMpKXG7%9HK zVRb&_`9Vd_;yIF@P-;?j;@TfL$lG zZ|}{FicXrdFHHMJjxS`+3Eia z8i4|ppC(WI6NC*L?~8wA%gX8tF>JqXunbG13_^`AB-^Ammwdvh$jbX%uT_OZZk~}A z5^EzoDG!GlY-8LOoJ>(6blT)V*&mho) z|7)LuHPX&*ShcsAaGFTb&Sa>I$Wj&fZ9RHx&|;__!_iaxuz`wT{`$shrmpXiY50iD zVqKIrS?A793nT4>RqO~VC$lc?ain(Z!s`efsgrn}Twh&m3<-)rwwsy9^~`}to-?zK zia5Fi;sau>;1jt0tvnzx)I3 zALORXf2bmaM8M!`UQecYgqd%I8N|ZZqejMAz!Fr5T$W5!%{B)A^B#ivhYEHd8}FE1 z-V~!;WG#Gtn*@8}w<)rd)k8P5QWa!5>E1BKvL2*ys83H*PZ0d|&AxrCYDi>Is=nm8 zQJ+9)&@L7;!%osH^oi4bmbPn8mqOY_OQ&5N9{KE`B26N@(>eJWPTbRLsaJHJS_L_?$mpVBzxwc|Q|JcE8P z@%EwQAO^d&R5?_kOSLOK3nOGZ{=|?^$Unx2JJyIBJ@CDW_}$qQdm*3sh*J&64CE%L z0@4VLJWt-iSS^^FpQ5y17u?yFCC8h^#97E$N7)fft(96!*$^(FMPjC8`KoNoYI14! zPe#Gb_1-QUrLy#}^ttl#EAxtU?9zk{yC?y`FnLiv2P8M&4WK6Zg8RVZdZ~*XSC_e3 z=y18=UQ_*&<)LZTW{4`h_1ZC;n!rscZeX58rogsJw|JWAojwhWxS1U-9@NZe46Ikt zmfa{*x?G9M+A(+UEuqVq2S^(A5x~AbId1NmKVqH9dhU%J@NgzQl>&YSgiU3Wtuvw zFE``_HU6E?ztHCW7v=4H6skI18Gv{rf7%jy3VhKQ1Oy?r^)Y;{-lKtJi2f94BE(MG ztUP2%NMcRHzU4NVW6;;AlKF-=B~;+}Rbyd$FOUHh+df>%=0ZXQ?;zOUqCB6|_7KIi|a z{Cmwy&epx>Jh-!kN!<=q^z4TyQ|$qS@hm>p3&;lHJO2OfZ~b;tIv~w=it(_mE{kY?H^>DFKE4gS&7IQ%gLGr=a9}x!EYp zbB|B=Boj{j!lOrUxufMpedQa4<7yz1i+Jz=l&|K$`H`S)&!^{hHz((K=T9g&jA6X! zr(eJ57{g3#QQj;Q$J9aA*CXDBZ*@iPwzm)N%9~>gw(xoUd+wH?CFS_3=3JPG$4(y} zd&YI_Ko~J}*XQxH_Qah5=?rNmgw@{vK>oHUx`I6<{@Y69$)z|%W=-#wy!5Dq#^bNk z8s(15Z&(qI@R_Da;b9Au&;nC^@1<(CHLQ3l-+j_g(o&eex|If(*>-0tYw$Zu8 z2_60X6H(czI%|&CbQOD~-6csZT4I@+qA_R4_06X{CCWczs2f^*C=yZ0O`jV*D35sG z-lq&wJ?N)zGGL0_tmOiO#!tet0bez}3dPim2wuS?+sGH$4 z(@P*0aG??gJuo6$ckCuqI2*=zz^VFWQAl*Kn?_=fy>i`ew>=7H+X@;I)IkUnaS zNq?1^o_iy$`!+4zbE;#pccsp&eS@;MwemyqSpt9aB~v4X8Ws!rwpNe|!Oywp8CwWx zASQj^VA0p{q`01Cf!L`yw|I{5U2~a=-d{|eOeSBm;(afu0&+$qoht1kdPlzLdVz}^CWMXHtba&Pl;6Ax!Hb%vUVA>;CCfn$zbIwC_8xl>K_ef$p2B=E$`4> zz{H%sFWdm<>wX4khkzlvds1QN9SIb04Syi|emFz}8m|s!{9_lZTAf`uP6$ra8+p!r zP*U(G53)D^GSmFgZ9TtWeaDn9gdI=3SkhQ@6zUu|mU!HITJ4Y`gj0rEMuKpwY8}5* zI(w6}4rM$!pe>Je&bS_@E-v~TM_v$xu#n7B(@%OOFgG3nQiPDGynFy+ml?FOt?R^# zISl#B8@lcu9O^Ini}7yF?UScmrKdd2hZOm$8lCSQet?#ZrWakp+hMiNAHUS}bmFcf z706$5ic~wr{tJZP8bGoV)Vo_22A7D#cyrwG{Ri{A_0b zh=;9|%{`9CBR5Tj!J2MOj_G-8&bm9Y|%vFg$h zn$QZJw#iM*-3{2+W8mF#-rP4GE+@7f;Jtsg?k)a7hrc-RCVgiXMRASyZHV-=n*O&> z_vA!JMt>nFp;9BeannOY>Q@SLmtL7^bwuk;)ynYq$=(^kg`}q5XDVy|+@qq0Oh%DN z#TSyw^#m(`oaYN!&iDD+cNd@$82{b}W88w|%g;9LxzBd@0~+RAk5?ITxeu;W%K3kL zPUbesmV7jXb08ZAbfmKD=*`w>lN8}qZG)}b@C6(l3DH9mNQ7vWSMq5TUT>nXB`3%V zb$*D4wB&KVEzOGQ+nUWlHuUl1RFMd0Nqj6YgJyo`MS~Lf#VP?=4J z?6u2>`G6`8bRSq-T;hDNW9~vK#dZDT%Ic>U4UwOIs>Q@eaVI)3Qbx3?A2-Kj>e`r5L;cRC2v~z0 z(viQZuQ|StlzN6QgcOBddr`h&sRb=rgMd~z^g@mWU=((TPUe5)=`Rxc$3*d5SUfk& z*^l&yCVml6_z*i;7x!3VJUXpk6T(KhEGjN_NugT}n|H)UKFC`dO9va#l{r&K?Gp0I zhMAON6_}QXbRWb+syAP=bDY0Z*n!VP&fHtusOo)BZk_Yve)zoz^Z}Nyne4-+mtF!C zr_E9yV2$>Grk?Q^k}L945FUd2b(l!BDf67GNh~t{L;T271B~E}%_6Fm7H)OnwR=%5 z3-hF3XsiT6Dmw{&ch9^7uVPmL$EgxBOw*}zM#w+XEswHEUjMmG`m9VDw&#+O z@cE=OVvD=-jvxwIQ_!q^q1BoGWp%yqcVCVsEV3m{+}QtxoNXC+Y*;R}gJkv2e{H>c zIFT?}hwa-Bc{^!HnT6vYLnoP#fNIulMv*7Hrg{_eEEK{KKQ=}q5uOCo^Cd)c5+DWcr z>*&e(2sMibQB>AF{_ zrJ|(6>i~wfw7kp%d}2BfjRp8at`tgAjSo@S?p)yt^IYfDJ?q)8G8%$LHHo_5d#~k( zmpKX5^OBIc&8N**{g1E5VixcQ2}h60c+OvA@au3fbMz5?wX%eu=*Pc$cSl&aFS_FZchL*_VgI8FTWg`%9{(@y z2MQIXD{ewV4@Uya-;FgdzOeqo)l6(MDrJS>h#Ww&uQ@y0U~ z2O+`0KAk7gq>}j3$EZcv1XWi%s6>PfZxBF#C{22#dtxJeuwhfrrS8#ty1ckJ0dZi^ z{=Vf4oBTBd;e7Smj~6oBi(aM(P!UcXB6&B$N!K%~UqU8D3n`p(bK~`|;)5b>Z|Jak zv;@~7^}cA|p}zmF#KD8;ee7sLX04jp>WI98K(>Ed1s6l;lB5fDKBRY$t;^zw0v6{B zhZ9inp3YLc06i{5@DozaAk{w_FBZ$Ym}^+%C-DeHXqf^?g_4^Ob@sm>)TlRj!vo#D zWLc|DFxyQ=+@B4&rR$`}aL6iF^K#&5a44}Kp)XxH3hDf$L+X!z8)9-&E10Mm&4{TT zvU+tnIuSBnX}*o$cp-WM06Sna04=Dq)o%eQ2MBAKve-gk$A0|_ylDURpS4j^;8O3H z@jBlb5lra@N3JqlBsW3e(fWaQJxab4WT8XYlHXxvoBN*@94xgyL1L&^=M|R$Z6lHH(U^OndlQNDU=I?h~V@f{YY^;}LPe%DqTwszHpJ>KmOq$aL8k`!W! ziKo+xCwb*~A8C3vO4=YW-gcC_zdfK}=3t8OF6+%yc_`nsJZ1LXSJ6EI_-ES@hAz@YUbk zXfG27@psk9>42#aJ0c}fifmm@GjReD-wIan?Qv~17C~=8>)Kz7PR$S z#1QX0bG`Hh>B@sIdxfdBVCzd^(vbp`$<9n={+AC80BR~#N6odNic$exEXH`xT^;wa z$1Q8mLO#A7>&pZoiLuDe+H0Ad2GnJrVz@DbuLr84^l*Y_UQ0>Rb&XRVU;_T_yp zWBn6n5%jZ`)VWSfe29~$ef<-wNlMxo!U`3>2&?RI%T3}IFFs1V6#6H+1ksPIV>7H{ zFzzZJi`YJ}Ugrqx6v3?i>NBBbWc*ca&|B|(plY+ytN?Vj8$2kwj2_N*ezHEe-*CZT zq(+UA;6!eICoc(+k(k$Wzk1)D;RrPmDKzu_-4JbWn{*83BOlCma(m0&wWg8^PZl}K zEUAk2*FwefBqDK~AuX~he6pF9_n`%q0fD?VflUeL9vuWMhC#J;b;-%eg$Gq#Yy0sJ zSkm1m{*F)%lI^3|DoNgkS-6jft z%-t8bRmNh4kI>^Z_`Gy84q48m2quOZMF3J-?tZ>wu{974SV?PxUL5vlWM1Up(L)R{ ze4D+*b4UHmOT?aqV74EsZtf@ESwbYf6ZLdZ!VhtKDo1VFOTB%fqn}~t0KB~pfKk%y z^6ALW&o3uI?$)q|>4EpY+bhoCbFBoBZp0KX2lsLD+&%NR=|T^tU^?H4|7^Jt-8RiA z!opib+osX7OfFQL} zfq1kN;P29~uqcG05x+%naEKC4haTHTHUIdq=x33kS*`Cs5KCI6Ql?r0epF-OHeis1 zALx3S|9oj-|MKUP%~HMO{q^w+!n*3^BLT3g6+|`c03zWHC~bXYKA8qzmrtOdO)`qs zysvel@!71r!7pRJ49w?7IeMlzJ4cezF2GW@`>dnLv!l?d(d=*i} ziLo)5qt9TlUoYI`2J1ydfJF$^_1L4_{__VT&`;J2Fa!yBoxUd|CORK?Tx$xTL6IT} zCxwGC0KEG}a;UJlm|jRo4}hKOT@DxWKv?}3z>AJgce>Rqy+&U&G1t+ox{qJ`Zk3%T1q&`CMklxYo;G zV$dZJld>S)gaZuKcfBs3ms$V)`*%4BqFD_bL=Q<=n}x%2R{nUc13Un$J3onRI_t;V z+uMlna0S3_=kT~Nw^?hyUAYg2R%QZHPXN~}@a@Gs1O?aj9qa|GG_JrLL-naBRA_Sa z4-OuHx0O1$H3h&(D)REVK|yd}_)qV|qvTy&8i4_oli+lJ{5d`ju%#VIX=&frIy@cK z`Q?Xzcik^lF3GoDZME3^5j;FGF%S6FD%#rR;AO~jt+^G1pb?i|9?X|DH!tifKx=mc zu3&OpT<*NCr!1HO$|t(F+0a4gaAXWP*x28gbXxxc=8USXZ6!da2*AU`TXg$DqoAYD z=ZXZ@13DrOF0Q5k9h4*>!6M8<6Eb+6ORxT|mZ1`IemCe1@c@w!STM6?YVR2} zYQ?BOXV-fucfSWFNJB*>-W!TO4}5QUmz1B8{>CflFsc4c`|n`!dTb}FU97sUTb%EV z-GMkccV5p|hbdM1det=tmr?WE)?n;8E@m##;Yy1V3kwSYw|y>nEQtqbp8zPa!fktQ zY_M`bWTL?5aySY8&I3e5OMsX`dG$&OOu^1XPSp#`rjN`XAN*@US{mNv-i!nmwH(D0 z9hM7pFCq+vPY*MM!(v(tIK`P%J)aVwWg`x(cRm3Al@quu96oZPr=^tywg?7!eRaHE z^)kre;>%i%@9lYJqQmNCAMxtxI>kcmGQtaa`}XYxU}buIY4VH3e*H>hh|y?& z@mHaLV4zf^&NM^FSMXGAu?87V5W}w;uvJ0uQu_3GPcI;#116Gy=gYT7yYts=$no9GK|9J?(;vS!G~(IImyl^+&#;CnM|LkB8DWfFrXf>1B30YM3ttPomji znaocT`(ok(o>zat6^kl(Fm{(ib(_`JWV|3fZinANfSg`LXxl6@G&F<)fz*3YBq6TB z!2sX4!8Q9`LA1m3T2)Fa0u1xAFRrX;_xJUE2i~1Pz@r`jYOp~31=gnhkX$x@35kHm zcknhf8?Dm=0s_E>w;jag-wkG{`4c%gc`Hz(st2+`z;C_09k9$JVq@n)YUEq|I|({I zJ}rRIQjJBqaRRFgIHxi&ad&_Nk3&RcefbHh@+}w!qEf1?04x@lN&9{+$^qnn-OMsoW}z6q_!`CAXE5T*oaX)ln{4w{=++LmNeCrEE!Ba+w*K6kS#> ztxC~-3~h2PVbr2!$&Q@YN1fgCc$~-i4^F?%!)HF<&*$@gy`S&*clx{zbw^Z;z~@K@ zc>yykp)(kw>e$%BiHZF~bt=s-$^%^7-F-)jIJhh;&#N+;d=r<;m6eexjb5^imT>Lu ze2Fq`Mn+yujahi_6LUmFftP#ttrgcul6(Jp6W{u<3QJr2X`IO{Tw+^0xPy&h~I716&`+mT(Xxi ztz;pbV0e;^8L>LI*9;z}q*ya*^%~3koKsXxy#Now)I~;(CEJZ2X05fy=Cf^WkF`bj z>vs3_%sn)CzkdX;QO-SsNzjD{4{Sw(9qr;Df>#eQRS&3AakBUHm{uD$h;fsk2!dD^ z)?*jm_sj}$kaX-AyU;Fs_^xtLFZAcO>c>Q2agk$Q#t%RIgn?H3K0&wzm!D3*>2rOo zylWIEDt;Byva_=@`pxOR!#FIhK9cT(1v!%Nf6BZcwp4ZGd;N;=qVHh=de(`}pb3)<}4F8s^7rpks zTm7!BGE1>c4TKre_4QIsO(Hr93JQR-SSDXW4F zJao5~Ilu>K7`~tVEkV^wr33N*ec?0X@q5o>w+$N-h&?FZI{u{X$sr{(zA05=K%YIx zNZ6@mV8CWbGoQ47{9`Q^H{kXWg%?)?v$kIQ<|1-?|Ka&wPmq*!1(}h?AC&)Q-}-y^ z9N6p(kp``(V}@DHV<;{}sSj`eU6G54Hl;4&w^l^*ar}}}Q=4#$MIqeu&MFQw3jA20 zQFxDv5Q?H=k5znY3#%H?lo>^jDS?zn!*)I+JzYJmE^pI`Wqru_5NLQnJ^N9A_b{5B zyt1-Jv17Z6xOhPJ{>-`2s(GxPvY*jk$@P4TsO5yU84m4FY%{4>FUKoGi?> zO3uj4EChBHUAp9epD%ILH8co>kTUVkoja~+Nt38D@Q%9}RN;$L;4dslfi}hOoUez? zqiL;M=l;z#Pvm-7te-=EN4DX{OpJ2@5G1i6vI}!_)%n)tCPIVqZ~L557XHIYTU(oO z3x-keQZgfJG$tPx7bkvRng>75JUj&fWS;BAGYL#3^~{ER8bex3S9c4p0xXSYdu}b= zb(qB5>(~Dz?yRn9jPvR06iveWtV#&fXGm`l`vSvou`dz>lq4F+T)yj$>?c+7pXul4 z7mfQM%l^&2=2mwXO-<#xmR+GUdS8XGt1G}xYinp&0$=7XvpwF>{fr@Rd_I43rg+EX zc%~?CJpaGocitk}42U+B`2aZN({Sp{S3XG*zrEyc$y&q3&J#KURFN$N1Gs70o$b(F z((Mfaj)t)#vh{_ru_%6>xlfJ09+5|fWpas@*ItGuf2Mz56hQuMH*cWe0%#T>iS=ej zwmmJ~U3cVG!}09wix7Wp9L37$k@TI4h_V{_*&M1Kzda_04ftOc*s>uqLRj7;s(ralBvLa&6^MV)9)~1!5L#g@D+QvY_ujHR!jO z9s|{rRX8=G8)17No-^^ZPF<)by|UxVah#Y{wZzj-5u-yEuwtj3$)g_Zk@Y))v}H(hIw2}o9aaR0vJzJ0~#qPk>K z7OYOfFuP?G0YelES4a{X0_FPy;YoqvY;TgkvRdMVtl*mbrLF;_@7?oKOE1A>wqq%z z8(8AoQ4<%9AHf_krQ9-hzUa&mQLOU#czS+8f&Yj1-CHXnE{dB4g6r4Ufy}p{s19An z8m;yj1@b=1HsVJPw2`C4n&^oOCo}`gaGvE@@U(R$H{K37cj>aRr@O)aCBgfcSJdFlTeohZj(iFu;*6Xe0osc{?JI zR!0qbW8v+rtgKwqSZQ?{XgiQQzfdEl3;V!~W6DC$soGMiBN8ALj$VQkX6z?Bz|@_+ACPe8DJ%2u-P?%#wvk4oxk(13I5Q!rl;`{C z$`_Y?l9AyH=IOW+g-69rDnAI>1z>w-e|u!tZ5aaP5Lr?*IygAoh>0~b44OYUa;Uy) z=GeWgxAZ@DzobcdZ-#W}dR}dW(<#~l5!#^E6y>!|ukWAe4!!X$8Xf(Wqiw>Yb2aNF zU#8=r^RSDuvu24<7n?~s>pvgweUGuUq;~*Iu#6KY@-X5Q;-d~H3Utn7zdy=qLv$II z)^%{53%^iN5g@zx!LU(QPOi!#hfb|nOq)UVAl4im9na|PjISCO>kthJY`H4(464i@ z91JEW5v`t+b+3k$4i(zxxLR1mpS3F_s=#vLZOgZ(wXAZvYs7JlrMCG=sZtVntV7r} zOz=hsjR<}fZg*+jKTsKsb_iNcT*OEK%ES7H4<~XICEdoqAPFI~SWA}N&b6m`?)ntE z!3xY;2XbWm8?_F;zpqIEl1BxwP~DCfR?zCoA=@O9#0cl zO0tYTo1$W!n<}DS(uVhAF{qDs19-N3Z$!zM#adHzs$;LMP)9g?s_$`aH7YiWs@X%n zj8m=u@Urf!_ssL{8WfEhT)Rgphqi>8SiHH0qNp{si_9oWlG3cB Date: Tue, 6 Aug 2024 14:05:53 -0400 Subject: [PATCH 3/4] comparison table for sprott dynamics --- .../sprott_small_perturbation.py | 158 ++++++++++++++++++ 1 file changed, 158 insertions(+) create mode 100644 DYMAG_solver_version/sprott_small_perturbation.py diff --git a/DYMAG_solver_version/sprott_small_perturbation.py b/DYMAG_solver_version/sprott_small_perturbation.py new file mode 100644 index 0000000..83534f2 --- /dev/null +++ b/DYMAG_solver_version/sprott_small_perturbation.py @@ -0,0 +1,158 @@ +import nolds +import numpy as np +import matplotlib.pyplot as plt +import networkx as nx +import torch +from PDE_layer import PDE_layer +import phate +import time +import argparse +from data_module import RandomDataset +from torch_geometric.loader import DataLoader +from torch_geometric.data import Data, Dataset +import pandas as pd +import os +import multiprocessing as mp +import numpy as np + +def compute_lyapunov_for_node(traj): + # Calculate the maximum Lyapunov exponent for a given trajectory + return nolds.lyap_r(traj) + +def compute_lyapunov_for_feature(x, num_nodes, traj_ind): + lyap_max = 0 + # Calculate Lyapunov exponents for each node in the feature + for node_ind in range(num_nodes): + traj = x[:, node_ind, traj_ind] + l = compute_lyapunov_for_node(traj) + if l > lyap_max: + lyap_max = l + return lyap_max + +def measure_lyapunov(x): + # Assume that x is outputs + x = x.detach().cpu().numpy() + num_steps, n_nodes, num_features = x.shape + + # Prepare arguments for parallel processing + args = [(x, n_nodes, traj_ind) for traj_ind in range(num_features)] + + # Use multiprocessing to parallelize across features + with mp.Pool() as pool: + lyap_max_list = pool.starmap(compute_lyapunov_for_feature, args) + + return lyap_max_list + +def get_time_diffs(x1, x2): + # x1 and x2 have shape (num_steps, n_nodes, num_features) + diff = x1-x2 + # want the norm to be matrix norm, resulting in final shape of num_steps + norm = np.linalg.norm(diff, axis=(1,2)) + return norm + +if __name__ == '__main__': + argparser = argparse.ArgumentParser() + argparser.add_argument('--n_nodes', type=int, default=25, help='Number of nodes') + argparser.add_argument('--num_reps_per_n_node', type=int, default=5, help='Number of repetitions per number of nodes') + argparser.add_argument('--num_graphs', type=int, default=10, help='Number of graphs') + argparser.add_argument('--signal_type', type=str, default='dirac', help='Type of signal to use') + argparser.add_argument('--sampling_interval', type=float, default=0.2, help='Sampling interval') + argparser.add_argument('--final_t', type=float, default=20, help='Final time') + argparser.add_argument('--step_size', type=float, default=0.01, help='Step size') + argparser.add_argument('--edge_addition', type=int, default=2, help='Number of edges to add to the graph') + args = argparser.parse_args() + + # generate one graph + n_nodes = args.n_nodes + seed = 24 + torch.manual_seed(seed) + np.random.seed(seed) + signal_type = args.signal_type + + graph_params_bounds = {'n': (n_nodes, n_nodes), 'p': (.1, .5)} + # get dataset + data_unperturbed = RandomDataset(random_graph_model='er', + num_graphs=1, + graph_params_bounds=graph_params_bounds, + node_features=signal_type) + + # create a data object with the same graph but add a random edge + G = data_unperturbed[0] + data_unperturbed = data_unperturbed[0] + edge_index = G.edge_index + edge_addition = args.edge_addition + added_edges = 0 + while added_edges < edge_addition: + node0, node1 = torch.randint(0, n_nodes, (2,)) + edge_in_index = (node0 == edge_index[0]) & (node1 == edge_index[1]) + if edge_in_index.any(): + continue + edge_index = torch.cat([edge_index, torch.tensor([[node0, node1], [node1, node0]])], dim=1) + print(f'added edge between {node0} and {node1}') + added_edges += 1 + data_perturbed = Data(x=G.x, edge_index=edge_index) + + pde_layer_sprott = PDE_layer(dynamics='sprott', n_largest_graph = n_nodes,sampling_interval = args.sampling_interval, final_t = args.final_t, step_size = args.step_size) + pde_layer_heat = PDE_layer(dynamics='heat', n_largest_graph = n_nodes,sampling_interval = args.sampling_interval, final_t = args.final_t, step_size = args.step_size) + + batch = torch.tensor([0 for _ in range(n_nodes)], dtype=torch.long) + sprott_perturbed = pde_layer_sprott(data_perturbed.x, data_perturbed.edge_index, batch) + sprott_unperturbed = pde_layer_sprott(data_unperturbed.x, data_perturbed.edge_index, batch) + heat_perturbed = pde_layer_heat(data_perturbed.x, data_perturbed.edge_index, batch) + heat_unperturbed = pde_layer_heat(data_unperturbed.x, data_perturbed.edge_index, batch) + + # get lyapunov exponents + print('getting lyapunov for sprott perturbed') + lyap_sprott_perturbed = measure_lyapunov(sprott_perturbed) + print('getting lyapunov for sprott unperturbed') + lyap_sprott_unperturbed = measure_lyapunov(sprott_unperturbed) + print('getting lyapunov for heat perturbed') + lyap_heat_perturbed = measure_lyapunov(heat_perturbed) + print('getting lyapunov for heat unperturbed') + lyap_heat_unperturbed = measure_lyapunov(heat_unperturbed) + + # output has shape (num_steps, n_nodes, num_features) + deltas_sprott = get_time_diffs(sprott_perturbed, sprott_unperturbed) + deltas_heat = get_time_diffs(heat_perturbed, heat_unperturbed) + + # print out deltas sprott and heat every 5 steps + print('sprott deltas every second:') + deltas_sprott_sparse = deltas_sprott[::5] + print(deltas_sprott_sparse) + print('heat deltas every second:') + deltas_heat_sparse = deltas_heat[::5] + print(deltas_heat_sparse) + + + # plot the results + plt.figure() + plt.plot(deltas_sprott, label='Sprott') + plt.plot(deltas_heat, label='Heat') + plt.legend() + plt.xlabel('Time step') + plt.ylabel('Norm of difference') + plt.title('Difference between perturbed and unperturbed graphs') + plt.show() + + # make a 2x2 subplot and plot the solution over time at node 0 feature 0 + plt.figure() + plt.subplot(2,2,1) + plt.plot(sprott_perturbed[:, 0, 0]) + plt.title('Sprott perturbed') + plt.subplot(2,2,2) + plt.plot(sprott_unperturbed[:, 0, 0]) + plt.title('Sprott unperturbed') + plt.subplot(2,2,3) + plt.plot(heat_perturbed[:, 0, 0]) + plt.title('Heat perturbed') + plt.subplot(2,2,4) + plt.plot(heat_unperturbed[:,0,0]) + plt.title('Heat unperturbed') + plt.show() + + # print out lyapunov exponents + print('Lyapunov exponents:') + print('Sprott perturbed:', lyap_sprott_perturbed) + print('Sprott unperturbed:', lyap_sprott_unperturbed) + print('Heat perturbed:', lyap_heat_perturbed) + print('Heat unperturbed:', lyap_heat_unperturbed) From f1daae76d72eb467959c0487a026b98527f79f31 Mon Sep 17 00:00:00 2001 From: Charles Xu Date: Tue, 6 Aug 2024 14:06:29 -0400 Subject: [PATCH 4/4] fixed example for DYMAG --- DYMAG_solver_version/DYMAG.py | 7 +++++-- DYMAG_solver_version/visualize_trajectories.py | 18 +++++++++++++----- 2 files changed, 18 insertions(+), 7 deletions(-) diff --git a/DYMAG_solver_version/DYMAG.py b/DYMAG_solver_version/DYMAG.py index dd5493c..5146cc0 100644 --- a/DYMAG_solver_version/DYMAG.py +++ b/DYMAG_solver_version/DYMAG.py @@ -78,12 +78,15 @@ def reset_parameters(self): if __name__ == '__main__': # test the model num_nodes = 10 - num_features = 5 + num_features = 100 x = torch.randn(num_nodes, num_features) edge_index = torch.tensor([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9], [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]], dtype=torch.long) batch_index = torch.tensor([0, 0, 0, 0, 0, 1, 1, 1, 1, 1], dtype=torch.long) - model = DYMAG(num_features, 5, heat_derivative_func) + model = DYMAG(num_features, 1) print(model(x, edge_index, batch_index).shape) print(model) + # get the number of trainable parameters for the model + print(sum(p.numel() for p in model.parameters() if p.requires_grad)) + import pdb; pdb.set_trace() diff --git a/DYMAG_solver_version/visualize_trajectories.py b/DYMAG_solver_version/visualize_trajectories.py index a1fe076..e23f566 100644 --- a/DYMAG_solver_version/visualize_trajectories.py +++ b/DYMAG_solver_version/visualize_trajectories.py @@ -16,7 +16,8 @@ n_nodes = 25 p = 0.2 num_traj = 5 - dynamics = 'sprott' + dynamics = 'heat' + small_perturb = True if dynamics == 'sprott': phate_knn = 5 final_t = 100 @@ -38,7 +39,14 @@ G = nx.erdos_renyi_graph(n_nodes, p) edge_index = torch.tensor(list(G.edges)).t().contiguous() edge_index = torch.cat([edge_index, edge_index[[1, 0]]], dim=1) - x = torch.randn(n_nodes, num_traj) * 10 + if small_perturb: + x = torch.randn(n_nodes, 1) * 2 + # replicate this n_traj times + x = x.repeat(1, num_traj) + # add noise + x = x + torch.randn(n_nodes, num_traj) + else: + x = torch.randn(n_nodes, num_traj) # zero center the data x = x - x.mean(dim=0) batch = torch.tensor([0 for _ in range(n_nodes)], dtype=torch.long) @@ -164,11 +172,11 @@ plt.title(f'PHATE of {dynamics} dynamics') # Save the plot - plt.savefig(f'figs/phate_plot_{dynamics}_{final_t}_{sampling_interval}_{phate_knn}_{seed}.png') + plt.savefig(f'figs/phate_plot_{dynamics}_{final_t}_{sampling_interval}_{phate_knn}_{seed}_{small_perturb}.png') # Show the plot plt.show() plt.close() - print(f'saved to figs/phate_plot_{dynamics}_{final_t}_{sampling_interval}_{phate_knn}_{seed}.png') + print(f'saved to figs/phate_plot_{dynamics}_{final_t}_{sampling_interval}_{phate_knn}_{seed}_{small_perturb}.png') # Plot the results but color by the trajectory fig = plt.figure() @@ -196,6 +204,6 @@ plt.title(f'PHATE of {dynamics} dynamics') # Save the plot - plt.savefig(f'figs/phate_plot_traj_color_{dynamics}_{final_t}_{sampling_interval}_{phate_knn}_{seed}.png') + plt.savefig(f'figs/phate_plot_traj_color_{dynamics}_{final_t}_{sampling_interval}_{phate_knn}_{seed}_{small_perturb}.png') # Show the plot plt.show()