From be57ad58461f1e1a9d1334c48d398e4a27e151f3 Mon Sep 17 00:00:00 2001 From: Joseph Date: Mon, 24 Oct 2022 12:13:47 +0100 Subject: [PATCH 01/27] Change default parameter handling --- markov_builder/MarkovChain.py | 120 ++++++++++++++++++++----------- markov_builder/example_models.py | 6 +- 2 files changed, 81 insertions(+), 45 deletions(-) diff --git a/markov_builder/MarkovChain.py b/markov_builder/MarkovChain.py index cacf68e..9da9c77 100644 --- a/markov_builder/MarkovChain.py +++ b/markov_builder/MarkovChain.py @@ -22,24 +22,24 @@ class MarkovChain(): """ def __init__(self, states: list = [], state_attributes_class: - MarkovStateAttributes = None, seed: int = None, - name: str = None): + MarkovStateAttributes = None, seed: int = None, name: str = + None, open_state: str = None, rates: list = None, + rate_dictionary: dict = None, auxiliary_expression: str = + None, auxiliary_symbol: str = None, shared_variables_dict: dict + = None, auxiliary_params_dict: dict = None): # Initialise the graph representing the states. Each directed edge has # a `rate` attribute which is a string representing the transition rate # between the two nodes self.graph = nx.DiGraph() - if states: - self.graph.add_nodes_from(states) self.rates = set() # Initialise a random number generator for simulation. Optionally, a # seed can be specified. self.rng = default_rng(seed) self.name = name - self.shared_rate_variables = set() - + self.shared_variables = {} self.rate_expressions = {} self.default_values = {} @@ -53,7 +53,24 @@ def __init__(self, states: list = [], state_attributes_class: self.auxiliary_variable = 'auxiliary_expression' if not issubclass(self.state_attributes_class, MarkovStateAttributes): - raise Exception("state_attirbutes_class must be a dataclass") + raise Exception("state_attributes_class must be a subclass of MarkovStateAttributes") + + if states: + for state in states: + self.add_state(state) + + if rates: + for r in rates: + self.add_both_transitions(*r) + if shared_variables_dict: + self.parameterise_rates(rate_dictionary, shared_variables_dict) + + if states and open_state and rates and rate_dictionary and auxiliary_expression and\ + auxiliary_symbol and shared_variables_dict and auxiliary_params_dict: + open_state = self.get_state_symbol(open_state) + self.define_auxiliary_expression(sp.sympify(auxiliary_expression.format(open_state)), + auxiliary_symbol, + auxiliary_params_dict) def mirror_model(self, prefix: str, new_rates: bool = False) -> None: """ Duplicate all states and rates in the model such that there are two identical components. @@ -404,14 +421,14 @@ def sample_trajectories(self, no_trajectories: int, time_range: list = [0, 1], # default missing values to those in self.default_values param_dict = {param: param_dict[param] if param in param_dict - else self.default_values[param] + else {**self.default_values, **self.shared_variables}[param] for param in param_list} else: param_dict = self.default_values - if starting_distribution is None: - starting_distribution = np.around(np.array([no_trajectories] * no_nodes) / no_nodes) - starting_distribution[0] += no_trajectories - starting_distribution.sum() + if starting_distribution is None: + starting_distribution = np.around(np.array([no_trajectories] * no_nodes) / no_nodes) + starting_distribution[0] += no_trajectories - starting_distribution.sum() distribution = starting_distribution @@ -554,7 +571,7 @@ def substitute_rates(self, rates_dict: dict): d['label'] = d['rate'] d['rate'] = str(sp.sympify(d['rate']).subs(rates_dict)) - def parameterise_rates(self, rate_dict: dict, shared_variables: list = []) -> None: + def parameterise_rates(self, rate_dict: dict, shared_variables: dict = {}) -> None: """Define a set of parameters for the transition rates. Parameters declared as @@ -569,14 +586,6 @@ def parameterise_rates(self, rate_dict: dict, shared_variables: list = []) -> No """ - # Check that shared_variables is a list (and not a string!) - if isinstance(shared_variables, str): - raise TypeError("shared_variables is a string but must be a list") - - for v in shared_variables: - if v in self.reserved_names: - raise Exception('name %s is reserved', v) - # Validate rate dictionary for r in rate_dict: if r not in self.rates: @@ -584,13 +593,21 @@ def parameterise_rates(self, rate_dict: dict, shared_variables: list = []) -> No rate_expressions = {} default_values_dict = {} + for r in self.rates: if r in rate_dict: - if len(rate_dict[r]) == 2: + default_values = [] + dummy_variables = [] + if len(rate_dict[r]) == 1: + expression = rate_dict[r][0] + elif len(rate_dict[r]) == 2: expression, dummy_variables = rate_dict[r] - default_values = [] - else: + elif len(rate_dict[r]) == 3: expression, dummy_variables, default_values = rate_dict[r] + else: + raise ValueError(f"Rate dictionary was malformed. \ + Entry with key {r} ({rate_dict[r]}) should be a tuple/list of length 1-3") + if len(dummy_variables) < len(default_values): raise ValueError("dummy variables list and default values list have mismatching lengths.\ Lengths {} and {}".format(len(dummy_variables), len(default_values))) @@ -598,10 +615,11 @@ def parameterise_rates(self, rate_dict: dict, shared_variables: list = []) -> No expression = sp.sympify(expression) for symbol in expression.free_symbols: - variables = list(dummy_variables) + list(shared_variables) - if str(symbol) not in variables: - raise Exception( - f"Symbol, {symbol} was not found in dummy variables or shared_variables, {variables}.") + symbol = str(symbol) + variables = list(dummy_variables) + list(shared_variables.keys()) + if symbol not in variables: + # Add symbol to shared variables dictionary + shared_variables[symbol] = None subs_dict = {u: f"{r}_{u}" for i, u in enumerate(dummy_variables)} rate_expressions[r] = sp.sympify(expression).subs(subs_dict) @@ -609,13 +627,14 @@ def parameterise_rates(self, rate_dict: dict, shared_variables: list = []) -> No for u, v in zip(dummy_variables, default_values): new_var_name = f"{r}_{u}" if new_var_name in default_values_dict: - raise Exception(f"A parameter with label {new_var_name} is already present in the model.") + raise ValueError(f"A parameter with label {new_var_name} is already present in the model.") default_values_dict[new_var_name] = v self.rate_expressions = {**self.rate_expressions, **rate_expressions} - self.default_values = {**self.default_values, **default_values_dict} + self.default_values = default_values_dict - self.shared_rate_variables = self.shared_rate_variables.union(shared_variables) + for key in shared_variables.keys(): + self.default_values[key] = shared_variables[key] def get_parameter_list(self) -> List[str]: """ @@ -630,7 +649,7 @@ def get_parameter_list(self) -> List[str]: for symbol in self.rate_expressions[r].free_symbols: rates.add(str(symbol)) - rates = rates.union([str(sym) for sym in self.shared_rate_variables]) + rates = rates.union([str(sym) for sym in self.shared_variables]) return sorted(rates) @@ -703,9 +722,18 @@ def generate_myokit_model(self, name: str = "", comp.add_variable(parameter) if parameter in self.default_values: comp[parameter].set_rhs(self.default_values[parameter]) + elif parameter in self.auxiliary_variables: + comp[parameter].set_rhs(self.auxiliary_variables[parameter]) + elif parameter in self.shared_variables: + comp[parameter].set_rhs(self.shared_variables[parameter]) for rate in self.rates: - comp.add_variable(rate) + try: + comp.add_variable(rate) + except myokit.DuplicateName: + # Variable has already been added + pass + if rate in self.rate_expressions: expr = self.rate_expressions[rate] comp[rate].set_rhs(str(expr)) @@ -724,7 +752,7 @@ def generate_myokit_model(self, name: str = "", var.set_state_value(0) # All but one states start with 0 occupancy. If they were all 0 nothing - # would happen! + # would happen. comp[states[-1]].set_state_value(1) # Write equation for eliminated state using the fact that the state @@ -739,12 +767,14 @@ def generate_myokit_model(self, name: str = "", comp[self.get_state_symbol(eliminate_state)].set_rhs(rhs_str) # Add auxiliary equation if required - if self.auxiliary_expression is not None: + if self.auxiliary_expression is not None and self.auxiliary_variable: comp.add_variable(self.auxiliary_variable) comp[self.auxiliary_variable].set_rhs(str(self.auxiliary_expression)) + return model - def define_auxiliary_expression(self, expression: sp.Expr, label: str = None, default_values: dict = {}) -> None: + def define_auxiliary_expression(self, expression: sp.Expr, label: str = + None, default_values: dict = {}) -> None: """Define an auxiliary output variable for the model. :param expression: A sympy expression defining the auxiliary variable @@ -763,12 +793,6 @@ def define_auxiliary_expression(self, expression: sp.Expr, label: str = None, de raise Exception() state_symbols = [self.get_state_symbol(state) for state in self.graph.nodes()] - for symbol in expression.free_symbols: - if str(symbol) not in state_symbols: - if self.shared_rate_variables is None: - self.shared_rate_variables = set(str(symbol)) - else: - self.shared_rate_variables.add(str(symbol)) for symbol in default_values: symbol = sp.sympify(symbol) @@ -777,9 +801,21 @@ def define_auxiliary_expression(self, expression: sp.Expr, label: str = None, de if symbol in self.default_values: raise Exception() - self.default_values = {**self.default_values, **default_values} + for symbol in expression.free_symbols: + if str(symbol) not in state_symbols: + symbol = str(symbol) + if symbol in default_values.keys(): + self.shared_variables[symbol] = default_values[symbol] + elif symbol in self.shared_variables: + continue + else: + self.shared_variables[symbol] = np.NaN + self.auxiliary_expression = expression + # TODO check these variables are not elsewhere in the model + self.auxiliary_variables = default_values + def as_latex(self, state_to_remove: str = None, include_auxiliary_expression: bool = False, column_vector=True, label_order: list = None) -> str: """Creates a LaTeX expression describing the Markov chain, its parameters and diff --git a/markov_builder/example_models.py b/markov_builder/example_models.py index db9c65b..9380db3 100644 --- a/markov_builder/example_models.py +++ b/markov_builder/example_models.py @@ -63,7 +63,7 @@ def construct_four_state_chain(): 'k_3': positive_rate_expr + ((0.0873, 8.91E-3),), 'k_4': negative_rate_expr + ((5.15E-3, 0.003158),)} - mc.parameterise_rates(rate_dictionary, shared_variables=('V',)) + mc.parameterise_rates(rate_dictionary) open_state = mc.get_state_symbol('O') auxiliary_expression = sp.sympify(f"g_Kr * {open_state} * (V - E_Kr)") @@ -122,7 +122,7 @@ def construct_wang_chain(): 'b_1': negative_rate_expr + ((0.006497, 0.03268),) } - mc.parameterise_rates(rate_dictionary, shared_variables=('V',)) + mc.parameterise_rates(rate_dictionary) open_state = mc.get_state_symbol('O') @@ -222,7 +222,7 @@ def construct_kemp_model(): 'b_h': positive_rate_expr + ((2.70e-01, 1.58e-02),), } - mc.parameterise_rates(rate_dictionary, shared_variables=('V',)) + mc.parameterise_rates(rate_dictionary) open_state = mc.get_state_symbol('O') From 9d952fc0a477a009501a00fea3b4e8463a4b65df Mon Sep 17 00:00:00 2001 From: Joseph Date: Mon, 24 Oct 2022 12:15:35 +0100 Subject: [PATCH 02/27] Remove unnecessary line from test --- tests/test_MarkovChain.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_MarkovChain.py b/tests/test_MarkovChain.py index 2ce6fa0..fa74892 100644 --- a/tests/test_MarkovChain.py +++ b/tests/test_MarkovChain.py @@ -80,7 +80,6 @@ def test_construct_examples(self): name), show_parameters=True) nx.drawing.nx_agraph.write_dot(mc.graph, "%s_dotfile.dot" % name) - nx.drawing.nx_agraph.write_dot(mc.graph, "%s_dotfile.dot" % name) def test_parameterise_rates(self): """ From 2297bf97459a07d763a8e2665dc3630c884f359c Mon Sep 17 00:00:00 2001 From: Joseph Date: Mon, 24 Oct 2022 12:20:24 +0100 Subject: [PATCH 03/27] Add first 4 models. Add tests to generate .mmt and graph files --- markov_builder/models/__init__.py | 0 .../models/thirty_models/__init__.py | 7 +++ markov_builder/models/thirty_models/model0.py | 37 +++++++++++ markov_builder/models/thirty_models/model1.py | 32 ++++++++++ markov_builder/models/thirty_models/model2.py | 35 +++++++++++ markov_builder/models/thirty_models/model3.py | 40 ++++++++++++ .../thirty_models/tests/test_30_models.py | 63 +++++++++++++++++++ 7 files changed, 214 insertions(+) create mode 100644 markov_builder/models/__init__.py create mode 100644 markov_builder/models/thirty_models/__init__.py create mode 100644 markov_builder/models/thirty_models/model0.py create mode 100644 markov_builder/models/thirty_models/model1.py create mode 100644 markov_builder/models/thirty_models/model2.py create mode 100644 markov_builder/models/thirty_models/model3.py create mode 100644 markov_builder/models/thirty_models/tests/test_30_models.py diff --git a/markov_builder/models/__init__.py b/markov_builder/models/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/markov_builder/models/thirty_models/__init__.py b/markov_builder/models/thirty_models/__init__.py new file mode 100644 index 0000000..f97907f --- /dev/null +++ b/markov_builder/models/thirty_models/__init__.py @@ -0,0 +1,7 @@ +from .model0 import model_00 +from .model1 import model_01 +from .model2 import model_02 +from .model3 import model_03 + + +__all__ = ['model_00', 'model_01', 'model_02', 'model_03'] diff --git a/markov_builder/models/thirty_models/model0.py b/markov_builder/models/thirty_models/model0.py new file mode 100644 index 0000000..fb40244 --- /dev/null +++ b/markov_builder/models/thirty_models/model0.py @@ -0,0 +1,37 @@ +import numpy as np + +from markov_builder import MarkovChain +from markov_builder.rate_expressions import negative_rate_expr, positive_rate_expr + + +class model_00(MarkovChain): + description = "" + states = ('O', 'C', 'I', 'IC') + rates = [('O', 'C', 'k_2', 'k_1'), + ('I', 'IC', 'k_2', 'k_1'), + ('O', 'I', 'k_3', 'k_4'), + ('C', 'IC', 'k_3', 'k_4')] + + open_state = 'O' + shared_variables_dict = {'V': np.NaN} + + rate_dictionary = {'k_1': positive_rate_expr + ((2.26E-4, 6.99E-2),), + 'k_2': negative_rate_expr + ((3.44E-5, 5.460E-2),), + 'k_3': positive_rate_expr + ((0.0873, 8.91E-3),), + 'k_4': negative_rate_expr + ((5.15E-3, 0.003158),)} + auxiliary_expression = "g_Kr * {} * (V - E_Kr)" + auxiliary_symbol = 'I_Kr' + + auxiliary_params_dict = {'g_Kr': 0.1524, + 'E_Kr': -88 + } + + def __init__(self): + super().__init__(states=self.states, + open_state=self.open_state, + rates=self.rates, + rate_dictionary=self.rate_dictionary, + auxiliary_expression=self.auxiliary_expression, + auxiliary_symbol=self.auxiliary_symbol, + shared_variables_dict=self.shared_variables_dict, + auxiliary_params_dict=self.auxiliary_params_dict) diff --git a/markov_builder/models/thirty_models/model1.py b/markov_builder/models/thirty_models/model1.py new file mode 100644 index 0000000..124942b --- /dev/null +++ b/markov_builder/models/thirty_models/model1.py @@ -0,0 +1,32 @@ +from numpy import NaN + +from markov_builder.MarkovChain import MarkovChain +from markov_builder.rate_expressions import negative_rate_expr, positive_rate_expr + + +class model_01(MarkovChain): + description = "" + states = ('O', 'C') + rates = [('O', 'C', 'k_12', 'k_21')] + + open_state = 'O' + shared_variables_dict = {'V': NaN} + rate_dictionary = {'k_12': positive_rate_expr + ((2.26E-4, 6.99E-2),), + 'k_21': negative_rate_expr + ((3.45e-5, 0.05462),)} + + auxiliary_expression = "g_Kr * {} * (V - E_Kr)" + auxiliary_symbol = 'I_Kr' + + auxiliary_params_dict = {'g_Kr': 0.1524, + 'E_Kr': -88 + } + + def __init__(self): + super().__init__(states=self.states, + open_state=self.open_state, + rates=self.rates, + rate_dictionary=self.rate_dictionary, + auxiliary_expression=self.auxiliary_expression, + auxiliary_symbol=self.auxiliary_symbol, + shared_variables_dict=self.shared_variables_dict, + auxiliary_params_dict=self.auxiliary_params_dict) diff --git a/markov_builder/models/thirty_models/model2.py b/markov_builder/models/thirty_models/model2.py new file mode 100644 index 0000000..544e2e9 --- /dev/null +++ b/markov_builder/models/thirty_models/model2.py @@ -0,0 +1,35 @@ +from numpy import NaN + +from markov_builder.MarkovChain import MarkovChain +from markov_builder.rate_expressions import negative_rate_expr, positive_rate_expr + + +class model_02(MarkovChain): + description = "" + states = ('O', 'C', 'I') + rates = [('C', 'O', 'a_m', 'b_m'), + ('O', 'I', 'b', 'a')] + + open_state = 'O' + shared_variables_dict = {'V': NaN} + rate_dictionary = {'a_m': positive_rate_expr + ((2.26E-4, 6.99E-2),), + 'b_m': negative_rate_expr + ((3.44E-5, 5.460E-2),), + 'b': positive_rate_expr + ((0.0873, 8.91E-3),), + 'a': negative_rate_expr + ((5.15E-3, 0.003158),)} + + auxiliary_expression = "g_Kr * {} * (V - E_Kr)" + auxiliary_symbol = 'I_Kr' + + auxiliary_params_dict = {'g_Kr': 0.1524, + 'E_Kr': -88 + } + + def __init__(self): + super().__init__(states=self.states, + open_state=self.open_state, + rates=self.rates, + rate_dictionary=self.rate_dictionary, + auxiliary_expression=self.auxiliary_expression, + auxiliary_symbol=self.auxiliary_symbol, + shared_variables_dict=self.shared_variables_dict, + auxiliary_params_dict=self.auxiliary_params_dict) diff --git a/markov_builder/models/thirty_models/model3.py b/markov_builder/models/thirty_models/model3.py new file mode 100644 index 0000000..94dfd61 --- /dev/null +++ b/markov_builder/models/thirty_models/model3.py @@ -0,0 +1,40 @@ +from numpy import NaN + +from markov_builder.MarkovChain import MarkovChain +from markov_builder.rate_expressions import negative_rate_expr, positive_rate_expr + + +class model_03(MarkovChain): + description = "" + states = ('O', 'C', 'I', 'IC') + rates = [('O', 'C', 'b1', 'a1'), + ('I', 'IC', 'a3', 'b3'), + ('O', 'I', 'a4', 'b4'), + ('C', 'IC', 'a4', 'b4')] + + open_state = 'O' + shared_variables_dict = {'V': NaN} + rate_dictionary = {'a1': positive_rate_expr + ((2.26E-4, 6.99E-2),), + 'b1': negative_rate_expr + ((3.45E-5, 5.462E-2),), + 'a3': negative_rate_expr + ((5.15e-3, 0.03158),), + 'b3': ('(a3 * a1)/b1', (), ()), + 'a4': positive_rate_expr + ((0.08730, 8.91e-3),), + 'b4': negative_rate_expr + ((5.15e-3, 0.03158),) + } + + auxiliary_expression = "g_Kr * {} * (V - E_Kr)" + auxiliary_symbol = 'I_Kr' + + auxiliary_params_dict = {'g_Kr': 0.1524, + 'E_Kr': -88 + } + + def __init__(self): + super().__init__(states=self.states, + open_state=self.open_state, + rates=self.rates, + rate_dictionary=self.rate_dictionary, + auxiliary_expression=self.auxiliary_expression, + auxiliary_symbol=self.auxiliary_symbol, + shared_variables_dict=self.shared_variables_dict, + auxiliary_params_dict=self.auxiliary_params_dict) diff --git a/markov_builder/models/thirty_models/tests/test_30_models.py b/markov_builder/models/thirty_models/tests/test_30_models.py new file mode 100644 index 0000000..bdfd107 --- /dev/null +++ b/markov_builder/models/thirty_models/tests/test_30_models.py @@ -0,0 +1,63 @@ +#!/usr/bin/env python3 + +import logging +import os +import unittest + +import myokit +import networkx as nx + +from markov_builder.models.thirty_models import ( + model_00, + model_01, + model_02, + model_03, +) + + +class TestThirtyModels(unittest.TestCase): + + def setUp(self): + """Run by unittest before the tests in this class are performed. + + Create an output directory (if it doesn't already exist). An + alternative output path can be used by setting the + MARKOVBUILDER_TEST_OUTPUT environment variable + + """ + test_output_dir = os.environ.get('MARKOVBUILDER_TEST_OUTPUT', os.path.join( + os.path.dirname(os.path.realpath(__file__)), self.__class__.__name__)) + if not os.path.exists(test_output_dir): + os.makedirs(test_output_dir) + self.output_dir = test_output_dir + logging.info("outputting to " + test_output_dir) + + self.models = [model_00, model_01, model_02, model_03] + + def test_generate_myokit(self): + for model in self.models: + name = model.__name__ + logging.debug(f"intiating {name}") + mc = model() + myokit_model = mc.generate_myokit_model() + myokit.save(os.path.join(self.output_dir, f"{model.__name__}.mmt"), myokit_model) + return + + def test_visualise_graphs(self): + for model in self.models: + name = model.__name__ + logging.debug(f"intiating {name}") + mc = model() + + mc.draw_graph(os.path.join(self.output_dir, f"{name}_graph.html"), + show_parameters=False) + mc.draw_graph(os.path.join(self.output_dir, f"{name}_graph_with_parameters.html"), + show_parameters=True) + + nx.drawing.nx_agraph.write_dot(mc.graph, os.path.join(self.output_dir, + "%s_dotfile.dot" % name)) + + +if __name__ == "__main__": + logging.getLogger().setLevel(logging.DEBUG) + unittest.main() From ae770abc3fe8d34913ced5a8ae704d78d3d0e858 Mon Sep 17 00:00:00 2001 From: Joseph Date: Tue, 25 Oct 2022 11:56:01 +0100 Subject: [PATCH 04/27] Add tests for connectedness as reversibility. Add models 4,5 --- markov_builder/MarkovChain.py | 26 +++++++++- .../models/thirty_models/__init__.py | 4 +- markov_builder/models/thirty_models/model4.py | 52 +++++++++++++++++++ markov_builder/models/thirty_models/model5.py | 50 ++++++++++++++++++ .../thirty_models/tests/test_30_models.py | 20 ++++++- 5 files changed, 148 insertions(+), 4 deletions(-) create mode 100644 markov_builder/models/thirty_models/model4.py create mode 100644 markov_builder/models/thirty_models/model5.py diff --git a/markov_builder/MarkovChain.py b/markov_builder/MarkovChain.py index 9da9c77..a151a5b 100644 --- a/markov_builder/MarkovChain.py +++ b/markov_builder/MarkovChain.py @@ -485,6 +485,18 @@ def get_equilibrium_distribution(self, param_dict: dict = None) -> Tuple[List[st ss = np.append(ss, 1 - ss.sum()) return labels, ss + def is_connected(self) -> bool: + """Checks if the graph is strongly connected that is, if each state can be + reached from every other state. This function returns true even if all + transition rates are 0. + + :return: A bool which is true if the graph is strongly connected and + false otherwise + + """ + + return nx.algorithms.components.is_strongly_connected(self.graph) + def is_reversible(self) -> bool: """Checks symbolically whether or not the Markov chain is reversible for any set of non-zero transition rate values. @@ -498,7 +510,7 @@ def is_reversible(self) -> bool: # Digraph must be strongly connected in order for the chain to be # reversible. In other words it must be possible to transition from any # state to any other state in some finite number of transitions - if not nx.algorithms.components.is_strongly_connected(self.graph): + if not self.is_connected(): return False undirected_graph = self.graph.to_undirected(reciprocal=False, as_view=True) @@ -512,6 +524,13 @@ def is_reversible(self) -> bool: forward_rate_list = [sp.sympify(self.graph.get_edge_data(frm, to)['rate']) for frm, to in iter] backward_rate_list = [sp.sympify(self.graph.get_edge_data(frm, to)['rate']) for to, frm in iter] + logging.debug(self.rate_expressions) + + # Substitute in expressions + forward_rate_list = [rate.subs(self.rate_expressions) for rate in forward_rate_list] + backward_rate_list = [rate.subs(self.rate_expressions) for rate in + backward_rate_list] + logging.debug("Rates moving forwards around the cycle are: %s", forward_rate_list) logging.debug("Rates moving backwards around the cycle are: %s", backward_rate_list) @@ -630,7 +649,10 @@ def parameterise_rates(self, rate_dict: dict, shared_variables: dict = {}) -> No raise ValueError(f"A parameter with label {new_var_name} is already present in the model.") default_values_dict[new_var_name] = v - self.rate_expressions = {**self.rate_expressions, **rate_expressions} + rate_expressions = {rate: expr.subs(rate_expressions) for rate, expr in + rate_expressions.items()} + + self.rate_expressions = rate_expressions self.default_values = default_values_dict for key in shared_variables.keys(): diff --git a/markov_builder/models/thirty_models/__init__.py b/markov_builder/models/thirty_models/__init__.py index f97907f..2ab1bf1 100644 --- a/markov_builder/models/thirty_models/__init__.py +++ b/markov_builder/models/thirty_models/__init__.py @@ -2,6 +2,8 @@ from .model1 import model_01 from .model2 import model_02 from .model3 import model_03 +from .model4 import model_04 +from .model5 import model_05 -__all__ = ['model_00', 'model_01', 'model_02', 'model_03'] +__all__ = ['model_00', 'model_01', 'model_02', 'model_03', 'model_04', 'model_05'] diff --git a/markov_builder/models/thirty_models/model4.py b/markov_builder/models/thirty_models/model4.py new file mode 100644 index 0000000..88dc1da --- /dev/null +++ b/markov_builder/models/thirty_models/model4.py @@ -0,0 +1,52 @@ +from numpy import NaN + +from markov_builder.MarkovChain import MarkovChain + + +class model_04(MarkovChain): + description = "" + states = ('O', 'C', 'I', 'IC') + rates = [('O', 'C', 'b1', 'a1'), + ('I', 'IC', 'a3', 'b3'), + ('O', 'I', 'a4', 'b4'), + ('C', 'IC', 'a4', 'b4')] + + open_state = 'O' + shared_variables_dict = {'V': NaN, + 'p1': 2.26e-4, + 'p2': 0.06990, + 'p3': 3.45e-5, + 'p4': 0.05462, + 'p5': 0.08730, + 'p6': 8.91e-3, + 'p7': 5.15e-3, + 'p8': 0.03158, + 'p9': 5.15e-3, + 'p10': 0.03158, + 'p11': 0.15240, + } + + rate_dictionary = {'a1': ('p1 * exp(p2*V)',), + 'b1': ('p3 * exp(-p4*V)',), + 'a4': ('p5 * exp(p6*V)',), + 'b4': ('p7 * exp(-p8*V)',), + 'a3': ('p9 * exp(-p10*V)',), + 'b3': ('(a3*a1)/b1',) + } + + auxiliary_expression = "g_Kr * {} * (V - E_Kr)" + auxiliary_symbol = 'I_Kr' + + auxiliary_params_dict = {'g_Kr': 0.1524, + 'E_Kr': -88 + } + + def __init__(self): + super().__init__(states=self.states, + open_state=self.open_state, + rates=self.rates, + rate_dictionary=self.rate_dictionary, + auxiliary_expression=self.auxiliary_expression, + auxiliary_symbol=self.auxiliary_symbol, + shared_variables_dict=self.shared_variables_dict, + auxiliary_params_dict=self.auxiliary_params_dict) diff --git a/markov_builder/models/thirty_models/model5.py b/markov_builder/models/thirty_models/model5.py new file mode 100644 index 0000000..3df3bf3 --- /dev/null +++ b/markov_builder/models/thirty_models/model5.py @@ -0,0 +1,50 @@ +from numpy import NaN + +from markov_builder.MarkovChain import MarkovChain + + +class model_05(MarkovChain): + description = "" + states = ('O', 'C', 'I', 'IC') + rates = [('O', 'C', 'bm', 'am'), + ('I', 'IC', 'bm', 'am'), + ('O', 'I', 'b1', 'a1'), + ('C', 'IC', 'b2', 'a2')] + + open_state = 'O' + shared_variables_dict = {'V': NaN, + 'p1': 2.26e-4, + 'p2': 0.06990, + 'p3': 3.45e-5, + 'p4': 0.05462, + 'p5': 0.08730, + 'p6': 8.91e-3, + 'p7': 5.15e-3, + 'p8': 0.03158, + 'p9': 5.15e-3, + 'p10': 0.03158, + 'p11': 0.15240, + } + + rate_dictionary = {'am': ('p1 * exp(p2*V)',), + 'bm': ('p3 * exp(-p4*V)',), + 'b1': ('p5 * exp(p6*V)',), + 'a1': ('p7 * exp(-p8*V)',), + 'b2': ('p9 * exp(p10*V)',), + 'a2': ('(b2*a1)/b1',) + } + + auxiliary_expression = "p11 * {} * (V - E_Kr)" + auxiliary_symbol = 'I_Kr' + + auxiliary_params_dict = {} + + def __init__(self): + super().__init__(states=self.states, + open_state=self.open_state, + rates=self.rates, + rate_dictionary=self.rate_dictionary, + auxiliary_expression=self.auxiliary_expression, + auxiliary_symbol=self.auxiliary_symbol, + shared_variables_dict=self.shared_variables_dict, + auxiliary_params_dict=self.auxiliary_params_dict) diff --git a/markov_builder/models/thirty_models/tests/test_30_models.py b/markov_builder/models/thirty_models/tests/test_30_models.py index bdfd107..84beb09 100644 --- a/markov_builder/models/thirty_models/tests/test_30_models.py +++ b/markov_builder/models/thirty_models/tests/test_30_models.py @@ -12,6 +12,8 @@ model_01, model_02, model_03, + model_04, + model_05, ) @@ -32,7 +34,7 @@ def setUp(self): self.output_dir = test_output_dir logging.info("outputting to " + test_output_dir) - self.models = [model_00, model_01, model_02, model_03] + self.models = [model_00, model_01, model_02, model_03, model_04, model_05] def test_generate_myokit(self): for model in self.models: @@ -57,6 +59,22 @@ def test_visualise_graphs(self): nx.drawing.nx_agraph.write_dot(mc.graph, os.path.join(self.output_dir, "%s_dotfile.dot" % name)) + def test_connected(self): + for model in self.models: + name = model.__name__ + logging.debug(f"intiating {name}") + + mc = model() + self.assertTrue(mc.is_connected()) + + def test_reversible(self): + for model in self.models: + name = model.__name__ + logging.debug(f"intiating {name}") + + mc = model() + self.assertTrue(mc.is_reversible()) + if __name__ == "__main__": logging.getLogger().setLevel(logging.DEBUG) From 0b6f9d884ac296994b7c34ae4df6f059b935a61d Mon Sep 17 00:00:00 2001 From: Joseph Date: Tue, 1 Nov 2022 13:33:22 +0000 Subject: [PATCH 05/27] Fix draw_graph handling of transition rates --- markov_builder/MarkovChain.py | 15 +++-- markov_builder/models/thirty_models/model0.py | 4 +- markov_builder/models/thirty_models/model1.py | 2 +- markov_builder/models/thirty_models/model2.py | 13 ++-- markov_builder/models/thirty_models/model3.py | 18 +++--- markov_builder/models/thirty_models/model5.py | 2 +- .../thirty_models/tests/test_30_models.py | 60 +++++++++++++++++++ 7 files changed, 88 insertions(+), 26 deletions(-) diff --git a/markov_builder/MarkovChain.py b/markov_builder/MarkovChain.py index a151a5b..171bc99 100644 --- a/markov_builder/MarkovChain.py +++ b/markov_builder/MarkovChain.py @@ -560,7 +560,7 @@ def draw_graph(self, filepath: str = None, show_options: bool = if len(self.rate_expressions) == 0: raise Exception() else: - data['label'] = str(self.rate_expressions[data['rate']]) + data['label'] = str(sp.sympify(data['rate']).subs(self.rate_expressions)) nt = pyvis.network.Network(directed=True) nt.from_nx(self.graph) @@ -750,11 +750,14 @@ def generate_myokit_model(self, name: str = "", comp[parameter].set_rhs(self.shared_variables[parameter]) for rate in self.rates: - try: - comp.add_variable(rate) - except myokit.DuplicateName: - # Variable has already been added - pass + free_symbols = sp.parse_expr(rate).free_symbols + for symb in free_symbols: + if symb not in comp.variables(): + try: + comp.add_variable(str(symb)) + except myokit.DuplicateName: + # Variable has already been added + pass if rate in self.rate_expressions: expr = self.rate_expressions[rate] diff --git a/markov_builder/models/thirty_models/model0.py b/markov_builder/models/thirty_models/model0.py index fb40244..4c76d29 100644 --- a/markov_builder/models/thirty_models/model0.py +++ b/markov_builder/models/thirty_models/model0.py @@ -16,9 +16,9 @@ class model_00(MarkovChain): shared_variables_dict = {'V': np.NaN} rate_dictionary = {'k_1': positive_rate_expr + ((2.26E-4, 6.99E-2),), - 'k_2': negative_rate_expr + ((3.44E-5, 5.460E-2),), + 'k_2': negative_rate_expr + ((3.45E-5, 5.460E-2),), 'k_3': positive_rate_expr + ((0.0873, 8.91E-3),), - 'k_4': negative_rate_expr + ((5.15E-3, 0.003158),)} + 'k_4': negative_rate_expr + ((5.15E-3, 0.03158),)} auxiliary_expression = "g_Kr * {} * (V - E_Kr)" auxiliary_symbol = 'I_Kr' diff --git a/markov_builder/models/thirty_models/model1.py b/markov_builder/models/thirty_models/model1.py index 124942b..76f4edf 100644 --- a/markov_builder/models/thirty_models/model1.py +++ b/markov_builder/models/thirty_models/model1.py @@ -7,7 +7,7 @@ class model_01(MarkovChain): description = "" states = ('O', 'C') - rates = [('O', 'C', 'k_12', 'k_21')] + rates = [('O', 'C', 'k_21', 'k_12')] open_state = 'O' shared_variables_dict = {'V': NaN} diff --git a/markov_builder/models/thirty_models/model2.py b/markov_builder/models/thirty_models/model2.py index 544e2e9..da4fcf1 100644 --- a/markov_builder/models/thirty_models/model2.py +++ b/markov_builder/models/thirty_models/model2.py @@ -7,15 +7,16 @@ class model_02(MarkovChain): description = "" states = ('O', 'C', 'I') - rates = [('C', 'O', 'a_m', 'b_m'), - ('O', 'I', 'b', 'a')] + rates = [('O', 'C', 'bm', 'am'), + ('O', 'I', 'ah', 'bh')] open_state = 'O' shared_variables_dict = {'V': NaN} - rate_dictionary = {'a_m': positive_rate_expr + ((2.26E-4, 6.99E-2),), - 'b_m': negative_rate_expr + ((3.44E-5, 5.460E-2),), - 'b': positive_rate_expr + ((0.0873, 8.91E-3),), - 'a': negative_rate_expr + ((5.15E-3, 0.003158),)} + rate_dictionary = {'am': positive_rate_expr + ((2.26E-4, 6.99E-2),), + 'bm': negative_rate_expr + ((3.45E-5, 5.462E-2),), + 'ah': positive_rate_expr + ((0.08730, 8.91e-3),), + 'bh': negative_rate_expr + ((5.15e-3, 0.03158),), + } auxiliary_expression = "g_Kr * {} * (V - E_Kr)" auxiliary_symbol = 'I_Kr' diff --git a/markov_builder/models/thirty_models/model3.py b/markov_builder/models/thirty_models/model3.py index 94dfd61..92fa2c2 100644 --- a/markov_builder/models/thirty_models/model3.py +++ b/markov_builder/models/thirty_models/model3.py @@ -7,19 +7,17 @@ class model_03(MarkovChain): description = "" states = ('O', 'C', 'I', 'IC') - rates = [('O', 'C', 'b1', 'a1'), - ('I', 'IC', 'a3', 'b3'), - ('O', 'I', 'a4', 'b4'), - ('C', 'IC', 'a4', 'b4')] + rates = [('O', 'C', 'bm', 'am'), + ('I', 'IC', 'bm', 'am'), + ('O', 'I', 'ah', 'bh'), + ('C', 'IC', 'ah', 'bh')] open_state = 'O' shared_variables_dict = {'V': NaN} - rate_dictionary = {'a1': positive_rate_expr + ((2.26E-4, 6.99E-2),), - 'b1': negative_rate_expr + ((3.45E-5, 5.462E-2),), - 'a3': negative_rate_expr + ((5.15e-3, 0.03158),), - 'b3': ('(a3 * a1)/b1', (), ()), - 'a4': positive_rate_expr + ((0.08730, 8.91e-3),), - 'b4': negative_rate_expr + ((5.15e-3, 0.03158),) + rate_dictionary = {'am': positive_rate_expr + ((2.26E-4, 6.99E-2),), + 'bm': negative_rate_expr + ((3.45E-5, 5.462E-2),), + 'ah': positive_rate_expr + ((0.08730, 8.91e-3),), + 'bh': negative_rate_expr + ((5.15e-3, 0.03158),), } auxiliary_expression = "g_Kr * {} * (V - E_Kr)" diff --git a/markov_builder/models/thirty_models/model5.py b/markov_builder/models/thirty_models/model5.py index 3df3bf3..106513c 100644 --- a/markov_builder/models/thirty_models/model5.py +++ b/markov_builder/models/thirty_models/model5.py @@ -37,7 +37,7 @@ class model_05(MarkovChain): auxiliary_expression = "p11 * {} * (V - E_Kr)" auxiliary_symbol = 'I_Kr' - auxiliary_params_dict = {} + auxiliary_params_dict = {'E_Kr': -88} def __init__(self): super().__init__(states=self.states, diff --git a/markov_builder/models/thirty_models/tests/test_30_models.py b/markov_builder/models/thirty_models/tests/test_30_models.py index 84beb09..cfd8d56 100644 --- a/markov_builder/models/thirty_models/tests/test_30_models.py +++ b/markov_builder/models/thirty_models/tests/test_30_models.py @@ -7,6 +7,11 @@ import myokit import networkx as nx +import myokit as mk +import numpy as np + +import matplotlib.pyplot as plt + from markov_builder.models.thirty_models import ( model_00, model_01, @@ -75,6 +80,61 @@ def test_reversible(self): mc = model() self.assertTrue(mc.is_reversible()) + def test_output(self): + mmt_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), + 'models_myokit') + Erev = -88 + + comparison_plot_dir = os.path.join(self.output_dir, 'comparison_plots') + + if not os.path.exists(comparison_plot_dir): + os.makedirs(comparison_plot_dir) + + for i, model in enumerate(self.models): + name = model.__name__ + logging.debug(f"intiating {name}") + + mc = model() + mk_protocol = mk.load_protocol(os.path.join(mmt_dir, + 'simplified-staircase.mmt')) + + mk_model = mk.load_model(os.path.join(mmt_dir, + f"model-{i}.mmt")) + + sim = mk.Simulation(mk_model, mk_protocol) + sim.set_tolerance(1e-9, 1e-9) + sim.set_constant('nernst.EK', Erev) + + tmax = mk_protocol.characteristic_time() + times = np.linspace(0, tmax, int(tmax) + 1) + sim.pre(1000) + + log = sim.run(tmax, log_times=times, log=['ikr.IKr']) + + mk_IKr = np.array(log['ikr.IKr']) + + generated_mk_model = mc.generate_myokit_model() + + sim = mk.Simulation(generated_mk_model, mk_protocol) + sim.set_tolerance(1e-9, 1e-9) + sim.set_constant('markov_chain.E_Kr', Erev) + sim.pre(1000) + + log = sim.run(tmax, log_times=times, log=['markov_chain.I_Kr']) + gen_mk_IKr = np.array(log['markov_chain.I_Kr']) + + fig = plt.figure(figsize=(12, 9)) + fig.gca().plot(times[:-1], mk_IKr, label='original myokit simulation') + fig.gca().plot(times[:-1], gen_mk_IKr, label='generated myokit simulation') + + fig.gca().legend() + + fig.savefig(os.path.join(comparison_plot_dir, + f"{name}_myokit_comparison")) + + error = np.sqrt(np.mean((gen_mk_IKr - mk_IKr)**2)) + self.assertLess(error, 1e-2) + if __name__ == "__main__": logging.getLogger().setLevel(logging.DEBUG) From d50282716d794537ccf828c42c1f5a6fdc258ebf Mon Sep 17 00:00:00 2001 From: Joseph Date: Tue, 1 Nov 2022 13:35:45 +0000 Subject: [PATCH 06/27] Add more models --- .../models/thirty_models/__init__.py | 5 +- markov_builder/models/thirty_models/model6.py | 57 +++++++++++++++++++ markov_builder/models/thirty_models/model7.py | 48 ++++++++++++++++ markov_builder/models/thirty_models/model8.py | 54 ++++++++++++++++++ .../thirty_models/tests/test_30_models.py | 19 ++++--- 5 files changed, 174 insertions(+), 9 deletions(-) create mode 100644 markov_builder/models/thirty_models/model6.py create mode 100644 markov_builder/models/thirty_models/model7.py create mode 100644 markov_builder/models/thirty_models/model8.py diff --git a/markov_builder/models/thirty_models/__init__.py b/markov_builder/models/thirty_models/__init__.py index 2ab1bf1..9b1b43a 100644 --- a/markov_builder/models/thirty_models/__init__.py +++ b/markov_builder/models/thirty_models/__init__.py @@ -4,6 +4,9 @@ from .model3 import model_03 from .model4 import model_04 from .model5 import model_05 +from .model6 import model_06 +from .model7 import model_07 +from .model8 import model_08 -__all__ = ['model_00', 'model_01', 'model_02', 'model_03', 'model_04', 'model_05'] +__all__ = ['model_00', 'model_01', 'model_02', 'model_03', 'model_04', 'model_05', 'model_06', 'model_07', 'model_08'] diff --git a/markov_builder/models/thirty_models/model6.py b/markov_builder/models/thirty_models/model6.py new file mode 100644 index 0000000..73fa696 --- /dev/null +++ b/markov_builder/models/thirty_models/model6.py @@ -0,0 +1,57 @@ +from numpy import NaN + +from markov_builder.MarkovChain import MarkovChain + + +class model_06(MarkovChain): + description = "" + states = ('O', 'C', 'I', 'IC') + rates = [('O', 'C', 'a1', 'b1'), + ('O', 'I', 'a2', 'b2'), + ('I', 'IC', 'a3', 'b3'), + ('C', 'IC', 'a4', 'b4')] + + open_state = 'O' + shared_variables_dict = {'V': NaN, + 'p1': 2.26e-4, + 'p2': 0.06990, + 'p3': 3.45e-5, + 'p4': 0.05462, + 'p5': 0.08730, + 'p6': 8.91e-3, + 'p7': 5.15e-3, + 'p8': 0.03158, + 'p9': 5.15e-3, + 'p10': 0.03158, + 'p11': 5.15e-3, + 'p12': 0.03158, + 'p13': 5.15e-3, + 'p14': 0.03158, + 'p15': 0.15240 + } + + rate_dictionary = { + 'b1': ('p1 * exp(p2*V)',), + 'a1': ('p3 * exp(-p4*V)',), + 'a2': ('p5 * exp(p6*V)',), + 'b2': ('p7 * exp(-p8*V)',), + 'b3': ('p9 * exp(p10*V)',), + 'a3': ('p11 * exp(-p12*V)',), + 'a4': ('p13 * exp(p14*V)',), + 'b4': ('(b3*b2*a1*a4)/(b1*a2*a3)',) + } + + auxiliary_expression = "p15 * {} * (V - E_Kr)" + auxiliary_symbol = 'I_Kr' + + auxiliary_params_dict = {'E_Kr': -88} + + def __init__(self): + super().__init__(states=self.states, + open_state=self.open_state, + rates=self.rates, + rate_dictionary=self.rate_dictionary, + auxiliary_expression=self.auxiliary_expression, + auxiliary_symbol=self.auxiliary_symbol, + shared_variables_dict=self.shared_variables_dict, + auxiliary_params_dict=self.auxiliary_params_dict) diff --git a/markov_builder/models/thirty_models/model7.py b/markov_builder/models/thirty_models/model7.py new file mode 100644 index 0000000..fdfba85 --- /dev/null +++ b/markov_builder/models/thirty_models/model7.py @@ -0,0 +1,48 @@ +from numpy import NaN + +from markov_builder.MarkovChain import MarkovChain + + +class model_07(MarkovChain): + description = "" + states = ('O', 'C1', 'C2', 'I') + rates = [ + ('O', 'C1', 'bm*2', 'am'), + ('C1', 'C2', 'bm', 'am*2'), + ('O', 'I', 'a1', 'b1') + ] + + open_state = 'O' + shared_variables_dict = {'V': NaN, + 'p1': 2.26e-4, + 'p2': 0.06990, + 'p3': 3.45e-5, + 'p4': 0.05462, + 'p5': 0.08730, + 'p6': 8.91e-3, + 'p7': 5.15e-3, + 'p8': 0.03158, + 'p9': 0.15240, + } + + rate_dictionary = { + 'am': ('p1 * exp(p2*V)',), + 'bm': ('p3 * exp(-p4*V)',), + 'a1': ('p5 * exp(p6*V)',), + 'b1': ('p7 * exp(-p8*V)',), + } + + auxiliary_expression = "p9 * {} * (V - E_Kr)" + auxiliary_symbol = 'I_Kr' + + auxiliary_params_dict = {'E_Kr': -88} + + def __init__(self): + super().__init__(states=self.states, + open_state=self.open_state, + rates=self.rates, + rate_dictionary=self.rate_dictionary, + auxiliary_expression=self.auxiliary_expression, + auxiliary_symbol=self.auxiliary_symbol, + shared_variables_dict=self.shared_variables_dict, + auxiliary_params_dict=self.auxiliary_params_dict) diff --git a/markov_builder/models/thirty_models/model8.py b/markov_builder/models/thirty_models/model8.py new file mode 100644 index 0000000..bb2c3fe --- /dev/null +++ b/markov_builder/models/thirty_models/model8.py @@ -0,0 +1,54 @@ +from numpy import NaN + +from markov_builder.MarkovChain import MarkovChain + + +class model_08(MarkovChain): + description = "" + states = ('O', 'Y1', 'Y2', 'Y4') + rates = [ + ('O', 'Y2', 'k43', 'k34'), + ('O', 'Y4', 'k56', 'k65'), + ('Y1', 'Y2', 'k12', 'k21') + ] + + open_state = 'O' + shared_variables_dict = {'V': NaN, + 'p1': 2.26e-4, + 'p2': 0.06990, + 'p3': 3.45e-5, + 'p4': 0.05462, + 'p5': 0.08730, + 'p6': 8.91e-3, + 'p7': 5.15e-3, + 'p8': 0.03158, + 'p9': 0.08730, + 'p10': 8.91e-3, + 'p11': 5.15e-3, + 'p12': 0.03158, + 'p13': 0.15240 + } + + rate_dictionary = { + 'k12': ('p1 * exp( p2 * V)',), + 'k21': ('p3 * exp(-p4 * V)',), + 'k34': ('p5 * exp( p6 * V)',), + 'k43': ('p7 * exp(-p8 * V)',), + 'k56': ('p9 * exp( p10 * V)',), + 'k65': ('p11 * exp(-p12 * V)',) + } + + auxiliary_expression = "p13 * {} * (V - E_Kr)" + auxiliary_symbol = 'I_Kr' + + auxiliary_params_dict = {'E_Kr': -88} + + def __init__(self): + super().__init__(states=self.states, + open_state=self.open_state, + rates=self.rates, + rate_dictionary=self.rate_dictionary, + auxiliary_expression=self.auxiliary_expression, + auxiliary_symbol=self.auxiliary_symbol, + shared_variables_dict=self.shared_variables_dict, + auxiliary_params_dict=self.auxiliary_params_dict) diff --git a/markov_builder/models/thirty_models/tests/test_30_models.py b/markov_builder/models/thirty_models/tests/test_30_models.py index cfd8d56..cf67a8a 100644 --- a/markov_builder/models/thirty_models/tests/test_30_models.py +++ b/markov_builder/models/thirty_models/tests/test_30_models.py @@ -4,14 +4,12 @@ import os import unittest +import matplotlib.pyplot as plt import myokit -import networkx as nx - import myokit as mk +import networkx as nx import numpy as np -import matplotlib.pyplot as plt - from markov_builder.models.thirty_models import ( model_00, model_01, @@ -19,13 +17,16 @@ model_03, model_04, model_05, + model_06, + model_07, + model_08, ) class TestThirtyModels(unittest.TestCase): def setUp(self): - """Run by unittest before the tests in this class are performed. + """This is run by unittest before the tests in this class are performed. Create an output directory (if it doesn't already exist). An alternative output path can be used by setting the @@ -39,7 +40,8 @@ def setUp(self): self.output_dir = test_output_dir logging.info("outputting to " + test_output_dir) - self.models = [model_00, model_01, model_02, model_03, model_04, model_05] + self.models = [model_00, model_01, model_02, model_03, model_04, + model_05, model_06, model_07, model_08] def test_generate_myokit(self): for model in self.models: @@ -95,8 +97,9 @@ def test_output(self): logging.debug(f"intiating {name}") mc = model() - mk_protocol = mk.load_protocol(os.path.join(mmt_dir, - 'simplified-staircase.mmt')) + mk_protocol_filename = os.path.join(mmt_dir, + 'simplified-staircase.mmt') + mk_protocol = mk.load_protocol(mk_protocol_filename) mk_model = mk.load_model(os.path.join(mmt_dir, f"model-{i}.mmt")) From beacf8fc03759bfeddfb4b202525a991f4c0be9c Mon Sep 17 00:00:00 2001 From: Joseph Date: Fri, 28 Apr 2023 13:52:43 +0100 Subject: [PATCH 07/27] Fix get_transition_matrix. Fix reversibility check. Add show_html switch to draw_graph --- markov_builder/MarkovChain.py | 31 +++++++++++++++++++------------ markov_builder/__init__.py | 2 +- tests/test_MarkovChain.py | 6 +++--- 3 files changed, 23 insertions(+), 16 deletions(-) diff --git a/markov_builder/MarkovChain.py b/markov_builder/MarkovChain.py index 171bc99..62347ca 100644 --- a/markov_builder/MarkovChain.py +++ b/markov_builder/MarkovChain.py @@ -267,9 +267,13 @@ def get_transition_matrix(self, use_parameters: bool = False, else: edge = self.graph.get_edge_data(current_state, incident_state) if edge is not None: - row.append(edge["rate"]) + rate = edge["rate"] + if isinstance(rate, str): + row.append(edge["rate"]) + else: + row.append(edge['rate'][0]) else: - row.append(0) + row.append(sp.sympify('0')) matrix.append(row) matrix = sp.Matrix(matrix) @@ -285,7 +289,7 @@ def get_transition_matrix(self, use_parameters: bool = False, if label_order is None: return list(self.graph.nodes), matrix else: - if len(label_order) != self.graph.nodes(): + if len(label_order) != len(self.graph.nodes()): raise Exception("Not all states accounted for in label order") permutation = [label_order.index(n) for n in self.graph.nodes] @@ -498,7 +502,8 @@ def is_connected(self) -> bool: return nx.algorithms.components.is_strongly_connected(self.graph) def is_reversible(self) -> bool: - """Checks symbolically whether or not the Markov chain is reversible for any set of non-zero transition rate values. + """Checks symbolically whether or not the Markov chain is reversible for any set of non-zero + transition rate values. We assume that all transition rates are always non-zero and follow Colquhoun et al. (2004) https://doi.org/10.1529/biophysj.103. @@ -520,9 +525,9 @@ def is_reversible(self) -> bool: cycle.append(cycle[0]) logging.debug("Checking cycle {}".format(cycle)) - iter = list(zip(cycle, itertools.islice(cycle, 1, None))) - forward_rate_list = [sp.sympify(self.graph.get_edge_data(frm, to)['rate']) for frm, to in iter] - backward_rate_list = [sp.sympify(self.graph.get_edge_data(frm, to)['rate']) for to, frm in iter] + iterator = list(zip(cycle, itertools.islice(cycle, 1, None))) + forward_rate_list = [sp.sympify(self.graph.get_edge_data(frm, to)['rate']) for frm, to in iterator] + backward_rate_list = [sp.sympify(self.graph.get_edge_data(frm, to)['rate']) for to, frm in iterator] logging.debug(self.rate_expressions) @@ -538,19 +543,21 @@ def is_reversible(self) -> bool: logging.debug("Not all rates were specified.") return False - forward_rate_product = sp.prod(forward_rate_list) - backward_rate_product = sp.prod(backward_rate_list) - if(forward_rate_product - backward_rate_product).evalf() != 0: + forward_rate_product = sp.prod(forward_rate_list).subs(self.rate_expressions) + backward_rate_product = sp.prod(backward_rate_list).subs(self.rate_expressions) + if (forward_rate_product - backward_rate_product).evalf() != 0: return False return True def draw_graph(self, filepath: str = None, show_options: bool = - False, show_rates: bool = False, show_parameters: bool = False): + False, show_rates: bool = False, show_parameters: bool = False, + show_html=False): """Visualise the graph as a webpage using pyvis. :param filepath: An optional filepath to save the file to. If this is None, will be opened as a webpage instead. :param show_options: Whether or not the options menu should be displayed on the webpage :param show_parameters: Whether or not we should display the transition rates instead of their labels + :param show_html: Whether or not to open the outputted html file in the browser """ for _, _, data in self.graph.edges(data=True): @@ -569,7 +576,7 @@ def draw_graph(self, filepath: str = None, show_options: bool = nt.show_buttons() if filepath is not None: nt.save_graph(filepath) - else: + if show_html: nt.show('Markov_builder_graph.html') def substitute_rates(self, rates_dict: dict): diff --git a/markov_builder/__init__.py b/markov_builder/__init__.py index 8a1c82c..bb8cd4d 100644 --- a/markov_builder/__init__.py +++ b/markov_builder/__init__.py @@ -16,7 +16,7 @@ finally: # Always manually delete frame # https://docs.python.org/2/library/inspect.html#the-interpreter-stack - del(frame) + del frame # # Version info diff --git a/tests/test_MarkovChain.py b/tests/test_MarkovChain.py index fa74892..aaef5db 100644 --- a/tests/test_MarkovChain.py +++ b/tests/test_MarkovChain.py @@ -23,7 +23,7 @@ def setUp(self): """ test_output_dir = os.environ.get('MARKOVBUILDER_TEST_OUTPUT', os.path.join( - os.path.dirname(os.path.realpath(__file__)), self.__class__.__name__)) + 'test_output', self.__class__.__name__)) if not os.path.exists(test_output_dir): os.makedirs(test_output_dir) self.output_dir = test_output_dir @@ -162,11 +162,11 @@ def test_assert_reversibility_using_cycles(self): mc = example_models.construct_non_reversible_chain() logging.debug("graph is %s", mc.graph) - assert(not mc.is_reversible()) + self.assertFalse(mc.is_reversible()) logging.info("Checking reversibility of non-reversible chain") mc.add_open_trapping() logging.debug("graph is %s", mc.graph) - assert(not mc.is_reversible()) + self.assertFalse(mc.is_reversible()) def test_equate_rates(self): """ From 051a36fbea687b642fefab4264a1fa12a3e22ad2 Mon Sep 17 00:00:00 2001 From: Joseph Date: Wed, 3 May 2023 11:27:34 +0100 Subject: [PATCH 08/27] Add models 11 and 30. Remove duplicated test file --- .flake8 | 12 +- .../models/thirty_models/__init__.py | 5 +- .../models/thirty_models/model11.py | 66 ++++++++ .../models/thirty_models/model30.py | 116 ++++++++++++++ .../thirty_models/tests/test_30_models.py | 144 ------------------ setup.py | 5 +- 6 files changed, 197 insertions(+), 151 deletions(-) create mode 100644 markov_builder/models/thirty_models/model11.py create mode 100644 markov_builder/models/thirty_models/model30.py delete mode 100644 markov_builder/models/thirty_models/tests/test_30_models.py diff --git a/.flake8 b/.flake8 index b6b5463..c65fe30 100644 --- a/.flake8 +++ b/.flake8 @@ -1,10 +1,14 @@ [flake8] max_line_length = 120 ignore = - W391 # allow empty line at end of file - W503 # break before binary operator - allow either style - W504 # break after binary operator - allow either style - E226 # missing whitespace around arithmetic operator +# allow empty line at end of file + W391 +# break before binary operator - allow either style + W503 +# break after binary operator - allow either style + W504 +# missing whitespace around arithmetic operator + E226 exclude= .git, venv, diff --git a/markov_builder/models/thirty_models/__init__.py b/markov_builder/models/thirty_models/__init__.py index 9b1b43a..633b8c5 100644 --- a/markov_builder/models/thirty_models/__init__.py +++ b/markov_builder/models/thirty_models/__init__.py @@ -7,6 +7,9 @@ from .model6 import model_06 from .model7 import model_07 from .model8 import model_08 +from .model11 import model_11 +from .model30 import model_30 -__all__ = ['model_00', 'model_01', 'model_02', 'model_03', 'model_04', 'model_05', 'model_06', 'model_07', 'model_08'] +__all__ = ['model_00', 'model_01', 'model_02', 'model_03', 'model_04', + 'model_05', 'model_06', 'model_07', 'model_08', 'model_11', 'model_30'] diff --git a/markov_builder/models/thirty_models/model11.py b/markov_builder/models/thirty_models/model11.py new file mode 100644 index 0000000..13e2e5f --- /dev/null +++ b/markov_builder/models/thirty_models/model11.py @@ -0,0 +1,66 @@ +from numpy import NaN + +from markov_builder.MarkovChain import MarkovChain + + +class model_11(MarkovChain): + description = "" + states = ('O', 'I', 'IC1', 'IC2', 'C2', 'C1') + rates = [ + ('O', 'I', 'ah', 'bh'), + ('C1', 'O', 'a1', 'b1'), + ('IC1', 'I', 'a3', 'b3'), + ('IC2', 'IC1', 'a4', 'b4'), + ('C2', 'C1', 'a2', 'b2'), + ('C2', 'IC2', 'ah', 'bh'), + ('C1', 'IC1', 'ah', 'bh') + ] + + open_state = 'O' + shared_variables_dict = {'V': NaN, + 'p1': 2.26e-4, + 'p2': 0.06990, + 'p3': 3.45e-5, + 'p4': 0.05462, + 'p5': 0.08730, + 'p6': 8.91e-3, + 'p7': 5.15e-3, + 'p8': 0.03158, + 'p9': 2.26e-4, + 'p10': 0.06990, + 'p11': 3.45e-5, + 'p12': 0.05462, + 'p13': 0.08730, + 'p14': 8.91e-3, + 'p15': 5.15e-3, + 'p16': 0.03158, + 'p17': 0.15240 + } + + rate_dictionary = { + 'a1': ('p1 * exp( p2 * V)',), + 'b1': ('p3 * exp(-p4 * V)',), + 'ah': ('p5 * exp( p6 * V)',), + 'bh': ('p7 * exp(-p8 * V)',), + 'a2': ('p9 * exp( p10 * V)',), + 'b2': ('p11 * exp(-p12 * V)',), + 'a3': ('p13 * exp( p14 * V)',), + 'b3': ('(a3*b1)/a1',), + 'a4': ('p15 * exp(p16* V)',), + 'b4': ('(a4*b2)/a2',) + } + + auxiliary_expression = "p17 * {} * (V - E_Kr)" + auxiliary_symbol = 'I_Kr' + + auxiliary_params_dict = {'E_Kr': -85} + + def __init__(self): + super().__init__(states=self.states, + open_state=self.open_state, + rates=self.rates, + rate_dictionary=self.rate_dictionary, + auxiliary_expression=self.auxiliary_expression, + auxiliary_symbol=self.auxiliary_symbol, + shared_variables_dict=self.shared_variables_dict, + auxiliary_params_dict=self.auxiliary_params_dict) diff --git a/markov_builder/models/thirty_models/model30.py b/markov_builder/models/thirty_models/model30.py new file mode 100644 index 0000000..10b094a --- /dev/null +++ b/markov_builder/models/thirty_models/model30.py @@ -0,0 +1,116 @@ +from numpy import NaN + +from markov_builder.MarkovChain import MarkovChain + + +class model_30(MarkovChain): + description = "" + states = ('O', 'I', 'IC1', 'IC2', 'C2', 'C1', 'C3', 'IC3', 'IC4', 'C4') + rates = [ + ('C1', 'O', 'a3', 'b3'), + ('C4', 'C3', 'a12', 'b12'), + ('C3', 'C2', 'a1', 'b1'), + ('IC1', 'I', 'a10', 'b10'), + ('IC2', 'IC1', 'a9', 'b9'), + ('IC3', 'IC2', 'a8', 'b8'), + ('IC4', 'IC3', 'a11', 'b11'), + ('C2', 'C1', 'a2', 'b2'), + ('O', 'I', 'a4', 'b4'), + ('C2', 'IC2', 'a6', 'b6'), + ('C1', 'IC1', 'a5', 'b5'), + ('C4', 'IC4', 'a13', 'b13'), + ('C3', 'IC3', 'a7', 'b7') + ] + + open_state = 'O' + shared_variables_dict = {'V': NaN, + 'p1': 2.26e-4, + 'p2': 0.06990, + 'p3': 3.45e-5, + 'p4': 0.05462, + 'p5': 0.08730, + 'p6': 8.91e-3, + 'p7': 5.15e-3, + 'p8': 0.03158, + 'p9': 2.26e-4, + 'p10': 0.06990, + 'p11': 3.45e-5, + 'p12': 0.05462, + 'p13': 0.08730, + 'p14': 8.91e-3, + 'p15': 5.15e-3, + 'p16': 0.03158, + 'p17': 0.08730, + 'p18': 8.91e-3, + 'p19': 5.15e-3, + 'p20': 0.03158, + 'p21': 0.06990, + 'p22': 3.45e-5, + 'p23': 0.05462, + 'p24': 0.08730, + 'p25': 8.91e-3, + 'p26': 5.15e-3, + 'p27': 0.03158, + 'p28': 0.08730, + 'p29': 8.91e-3, + 'p30': 5.15e-3, + 'p31': 0.03158, + 'p32': 5.15e-3, + 'p33': 0.03158, + 'p34': 0.03158, + 'p35': 8.91e-3, + 'p36': 5.15e-3, + 'p37': 0.03158, + 'p38': 0.08730, + 'p39': 8.91e-3, + 'p40': 5.15e-3, + 'p41': 0.03158, + 'p42': 5.15e-3, + 'p43': 0.03158, + 'p44': 0.03158, + 'p45': 0.15240 + } + + rate_dictionary = { + 'a1': ('p1 * exp( p2 * V)',), + 'b1': ('p3 * exp(-p4 * V)',), + 'a2': ('p5 * exp( p6 * V)',), + 'b2': ('p7 * exp(-p8 * V)',), + 'a3': ('p9 * exp( p10 * V)',), + 'b3': ('p11 * exp(-p12* V)',), + 'a4': ('p13 * exp( p14 * V)',), + 'b4': ('p15 * exp(-p16 * V)',), + 'a8': ('p17 * exp(p18* V)',), + 'b8': ('p19 * exp(-p20 * V)',), + 'a9': ('p21 * exp(p22* V)',), + 'b9': ('p23 * exp(-p24 * V)',), + 'a10': ('p25 * exp(p26* V)',), + 'b10': ('p27 * exp(-p28 * V)',), + 'a5': ('p29 * exp(p30 * V)',), + 'b5': ('a10*a5*b4*(b3)/(a3*a4*b10)',), + 'a6': ('p31 * exp(p32 * V)',), + 'b6': ('a9*a6*b5*(b2)/(a2*a5*b9)',), + 'a7': ('p33 * exp(p34 * V)',), + 'b7': ('a8*a7*b1*(b6)/(a1*a6*b8)',), + 'a11': ('p35 * exp(p36* V)',), + 'b11': ('p37 * exp(-p38 * V)',), + 'a12': ('p39 * exp(p40* V)',), + 'b12': ('p41 * exp(-p42 * V)',), + 'a13': ('p43 * exp(p44 * V)',), + 'b13': ('a13*a11*b7*b12/(a12*a7*b11)',), + } + + auxiliary_expression = "p45 * {} * (V - E_Kr)" + auxiliary_symbol = 'I_Kr' + + auxiliary_params_dict = {'E_Kr': -85} + + def __init__(self): + super().__init__(states=self.states, + open_state=self.open_state, + rates=self.rates, + rate_dictionary=self.rate_dictionary, + auxiliary_expression=self.auxiliary_expression, + auxiliary_symbol=self.auxiliary_symbol, + shared_variables_dict=self.shared_variables_dict, + auxiliary_params_dict=self.auxiliary_params_dict) diff --git a/markov_builder/models/thirty_models/tests/test_30_models.py b/markov_builder/models/thirty_models/tests/test_30_models.py deleted file mode 100644 index cf67a8a..0000000 --- a/markov_builder/models/thirty_models/tests/test_30_models.py +++ /dev/null @@ -1,144 +0,0 @@ -#!/usr/bin/env python3 - -import logging -import os -import unittest - -import matplotlib.pyplot as plt -import myokit -import myokit as mk -import networkx as nx -import numpy as np - -from markov_builder.models.thirty_models import ( - model_00, - model_01, - model_02, - model_03, - model_04, - model_05, - model_06, - model_07, - model_08, -) - - -class TestThirtyModels(unittest.TestCase): - - def setUp(self): - """This is run by unittest before the tests in this class are performed. - - Create an output directory (if it doesn't already exist). An - alternative output path can be used by setting the - MARKOVBUILDER_TEST_OUTPUT environment variable - - """ - test_output_dir = os.environ.get('MARKOVBUILDER_TEST_OUTPUT', os.path.join( - os.path.dirname(os.path.realpath(__file__)), self.__class__.__name__)) - if not os.path.exists(test_output_dir): - os.makedirs(test_output_dir) - self.output_dir = test_output_dir - logging.info("outputting to " + test_output_dir) - - self.models = [model_00, model_01, model_02, model_03, model_04, - model_05, model_06, model_07, model_08] - - def test_generate_myokit(self): - for model in self.models: - name = model.__name__ - logging.debug(f"intiating {name}") - mc = model() - myokit_model = mc.generate_myokit_model() - myokit.save(os.path.join(self.output_dir, f"{model.__name__}.mmt"), myokit_model) - return - - def test_visualise_graphs(self): - for model in self.models: - name = model.__name__ - logging.debug(f"intiating {name}") - mc = model() - - mc.draw_graph(os.path.join(self.output_dir, f"{name}_graph.html"), - show_parameters=False) - mc.draw_graph(os.path.join(self.output_dir, f"{name}_graph_with_parameters.html"), - show_parameters=True) - - nx.drawing.nx_agraph.write_dot(mc.graph, os.path.join(self.output_dir, - "%s_dotfile.dot" % name)) - - def test_connected(self): - for model in self.models: - name = model.__name__ - logging.debug(f"intiating {name}") - - mc = model() - self.assertTrue(mc.is_connected()) - - def test_reversible(self): - for model in self.models: - name = model.__name__ - logging.debug(f"intiating {name}") - - mc = model() - self.assertTrue(mc.is_reversible()) - - def test_output(self): - mmt_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), - 'models_myokit') - Erev = -88 - - comparison_plot_dir = os.path.join(self.output_dir, 'comparison_plots') - - if not os.path.exists(comparison_plot_dir): - os.makedirs(comparison_plot_dir) - - for i, model in enumerate(self.models): - name = model.__name__ - logging.debug(f"intiating {name}") - - mc = model() - mk_protocol_filename = os.path.join(mmt_dir, - 'simplified-staircase.mmt') - mk_protocol = mk.load_protocol(mk_protocol_filename) - - mk_model = mk.load_model(os.path.join(mmt_dir, - f"model-{i}.mmt")) - - sim = mk.Simulation(mk_model, mk_protocol) - sim.set_tolerance(1e-9, 1e-9) - sim.set_constant('nernst.EK', Erev) - - tmax = mk_protocol.characteristic_time() - times = np.linspace(0, tmax, int(tmax) + 1) - sim.pre(1000) - - log = sim.run(tmax, log_times=times, log=['ikr.IKr']) - - mk_IKr = np.array(log['ikr.IKr']) - - generated_mk_model = mc.generate_myokit_model() - - sim = mk.Simulation(generated_mk_model, mk_protocol) - sim.set_tolerance(1e-9, 1e-9) - sim.set_constant('markov_chain.E_Kr', Erev) - sim.pre(1000) - - log = sim.run(tmax, log_times=times, log=['markov_chain.I_Kr']) - gen_mk_IKr = np.array(log['markov_chain.I_Kr']) - - fig = plt.figure(figsize=(12, 9)) - fig.gca().plot(times[:-1], mk_IKr, label='original myokit simulation') - fig.gca().plot(times[:-1], gen_mk_IKr, label='generated myokit simulation') - - fig.gca().legend() - - fig.savefig(os.path.join(comparison_plot_dir, - f"{name}_myokit_comparison")) - - error = np.sqrt(np.mean((gen_mk_IKr - mk_IKr)**2)) - self.assertLess(error, 1e-2) - - -if __name__ == "__main__": - logging.getLogger().setLevel(logging.DEBUG) - unittest.main() diff --git a/setup.py b/setup.py index ce11e52..bbfbe9d 100644 --- a/setup.py +++ b/setup.py @@ -55,8 +55,9 @@ 'plotly>=5.3', 'symengine>=0.8', 'sympy>=1.8', - 'pyvis==0.1.9', - 'myokit==1.33.0' + 'pyvis>=0.1.9', + 'myokit>=1.33.9', + 'pygraphviz>=1.10' ], extras_require={ 'test': [ From 5310b55946095b3de6968101329f7d182f3b8677 Mon Sep 17 00:00:00 2001 From: Joseph Date: Wed, 3 May 2023 11:39:03 +0100 Subject: [PATCH 09/27] Lint and isort --- markov_builder/MarkovChain.py | 2 +- tests/test_MarkovChain.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/markov_builder/MarkovChain.py b/markov_builder/MarkovChain.py index fa185ed..f70a99b 100644 --- a/markov_builder/MarkovChain.py +++ b/markov_builder/MarkovChain.py @@ -513,7 +513,7 @@ def is_connected(self) -> bool: def is_reversible(self) -> bool: """Checks symbolically if the Markov chain is reversible for any set of non-zero transition rate values. - + We assume that all transition rates are always non-zero and follow Colquhoun et al. (2004) https://doi.org/10.1529/biophysj.103. diff --git a/tests/test_MarkovChain.py b/tests/test_MarkovChain.py index 72331b8..b1bb69b 100644 --- a/tests/test_MarkovChain.py +++ b/tests/test_MarkovChain.py @@ -24,7 +24,7 @@ def setUp(self): """ test_output_dir = os.environ.get('MARKOVBUILDER_TEST_OUTPUT', os.path.join('test_output', - os.path.dirname(os.path.abspath(__file__)), + os.path.dirname(os.path.abspath(__file__)), self.__class__.__name__)) if not os.path.exists(test_output_dir): From 57415dcbd75442590519741aba2df2a791fbcf90 Mon Sep 17 00:00:00 2001 From: Joseph Date: Wed, 3 May 2023 11:57:01 +0100 Subject: [PATCH 10/27] Update test_parameterise_rates --- markov_builder/MarkovChain.py | 2 +- tests/test_MarkovChain.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/markov_builder/MarkovChain.py b/markov_builder/MarkovChain.py index f70a99b..598e05b 100644 --- a/markov_builder/MarkovChain.py +++ b/markov_builder/MarkovChain.py @@ -618,7 +618,7 @@ def parameterise_rates(self, rate_dict: dict, shared_variables: dict = {}) -> No V is the membrane voltage (a variable shared between transition rates). :param rate_dict: A dictionary with a 2-tuple containing an expression and dummy variables for each rate. - :param shared_variables: A list of variables that may be shared between transition rates + :param shared_variables: A dictionary of variables that may be shared between transition rates """ diff --git a/tests/test_MarkovChain.py b/tests/test_MarkovChain.py index b1bb69b..4eef5ef 100644 --- a/tests/test_MarkovChain.py +++ b/tests/test_MarkovChain.py @@ -101,7 +101,7 @@ def test_parameterise_rates_no_default(self): 'k_4': negative_rate_expr, } - mc.parameterise_rates(rate_dictionary, shared_variables=('V',)) + mc.parameterise_rates(rate_dictionary, shared_variables={'V': 'V'}) def test_myokit_output(self): """ From f199b6f5d808e1dbb94285a57e03cd5cd8b51961 Mon Sep 17 00:00:00 2001 From: Joseph Date: Wed, 3 May 2023 12:04:57 +0100 Subject: [PATCH 11/27] Add test_30_models.py --- tests/test_30_models.py | 179 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 179 insertions(+) create mode 100644 tests/test_30_models.py diff --git a/tests/test_30_models.py b/tests/test_30_models.py new file mode 100644 index 0000000..addd323 --- /dev/null +++ b/tests/test_30_models.py @@ -0,0 +1,179 @@ +#!/usr/bin/env python3 + +import itertools +import logging +import os +import unittest + +import matplotlib.pyplot as plt +import myokit +import myokit as mk +import networkx as nx +import numpy as np +import sympy as sp + +from markov_builder.models.thirty_models import ( + model_00, + model_01, + model_02, + model_03, + model_04, + model_05, + model_06, + model_07, + model_08, + model_11, + model_30, +) + + +class TestThirtyModels(unittest.TestCase): + + def setUp(self): + """This is run by unittest before the tests in this class are performed. + + Create an output directory (if it doesn't already exist). An + alternative output path can be used by setting the + MARKOVBUILDER_TEST_OUTPUT environment variable + + """ + test_output_dir = os.environ.get('MARKOVBUILDER_TEST_OUTPUT', os.path.join( + 'test_output', self.__class__.__name__)) + if not os.path.exists(test_output_dir): + os.makedirs(test_output_dir) + self.output_dir = test_output_dir + logging.info("outputting to " + test_output_dir) + + self.models = [ + model_00, + model_01, model_02, model_03, model_04, + model_05, model_06, model_07, model_08, + model_11, model_30 + ] + + self.model_indices = [0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 30] + + def test_generate_myokit(self): + for model in self.models: + name = model.__name__ + logging.debug(f"initiating {name}") + mc = model() + myokit_model = mc.generate_myokit_model() + myokit.save(os.path.join(self.output_dir, f"{model.__name__}.mmt"), myokit_model) + return + + def test_visualise_graphs(self): + for model in self.models: + name = model.__name__ + logging.debug(f"initiating {name}") + mc = model() + + mc.draw_graph(os.path.join(self.output_dir, f"{name}_graph.html"), + show_parameters=False) + mc.draw_graph(os.path.join(self.output_dir, f"{name}_graph_with_parameters.html"), + show_parameters=True) + + nx.drawing.nx_agraph.write_dot(mc.graph, os.path.join(self.output_dir, + "%s_dotfile.dot" % name)) + + def test_connected(self): + for model in self.models: + name = model.__name__ + logging.debug(f"initiating {name}") + + mc = model() + self.assertTrue(mc.is_connected()) + + def test_reversible(self): + for model in self.models: + name = model.__name__ + logging.debug(f"initiating {name}") + + mc = model() + if not mc.is_reversible(): + undirected_graph = mc.graph.to_undirected(reciprocal=False, as_view=True) + cycle_basis = nx.cycle_basis(undirected_graph) + + for cycle in cycle_basis: + iterator = list(zip(cycle, itertools.islice(cycle, 1, None))) + forward_rate_list = [sp.sympify(mc.graph.get_edge_data(frm, to)['rate']) for frm, to in iterator] + backward_rate_list = [sp.sympify(mc.graph.get_edge_data(frm, to)['rate']) for to, frm in iterator] + + logging.debug(mc.rate_expressions) + + # Substitute in expressions + forward_rate_list = [rate.subs(mc.rate_expressions) for + rate in forward_rate_list] + backward_rate_list = [rate.subs(mc.rate_expressions) for rate in + backward_rate_list] + + forward_rate_product = sp.prod(forward_rate_list) + backward_rate_product = sp.prod(backward_rate_list) + if (forward_rate_product - backward_rate_product).evalf() != 0: + logging.error("Rates moving forwards around the cycle are: %s", forward_rate_list) + logging.error("Rates moving backwards around the cycle are: %s", backward_rate_list) + + self.assertTrue(mc.is_reversible()) + + def test_output(self): + mmt_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), + 'models_myokit') + Erev = -88 + + comparison_plot_dir = os.path.join(self.output_dir, 'comparison_plots') + + if not os.path.exists(comparison_plot_dir): + os.makedirs(comparison_plot_dir) + + for i, model in zip(self.model_indices, self.models): + if i is None or model is None: + continue + name = model.__name__ + logging.debug(f"initiating {name}") + + mc = model() + mk_protocol_filename = os.path.join(mmt_dir, + 'simplified-staircase.mmt') + mk_protocol = mk.load_protocol(mk_protocol_filename) + + mk_model = mk.load_model(os.path.join(mmt_dir, + f"model-{i}.mmt")) + + sim = mk.Simulation(mk_model, mk_protocol) + sim.set_tolerance(1e-9, 1e-9) + sim.set_constant('nernst.EK', Erev) + + tmax = mk_protocol.characteristic_time() + times = np.linspace(0, tmax, int(tmax) + 1) + sim.pre(1000) + + log = sim.run(tmax, log_times=times, log=['ikr.IKr']) + + mk_IKr = np.array(log['ikr.IKr']) + + generated_mk_model = mc.generate_myokit_model() + + sim = mk.Simulation(generated_mk_model, mk_protocol) + sim.set_tolerance(1e-9, 1e-9) + sim.set_constant('markov_chain.E_Kr', Erev) + sim.pre(1000) + + log = sim.run(tmax, log_times=times, log=['markov_chain.I_Kr']) + gen_mk_IKr = np.array(log['markov_chain.I_Kr']) + + fig = plt.figure(figsize=(12, 9)) + fig.gca().plot(times[:-1], mk_IKr, label='original myokit simulation') + fig.gca().plot(times[:-1], gen_mk_IKr, label='generated markov_builder simulation') + + fig.gca().legend() + + fig.savefig(os.path.join(comparison_plot_dir, + f"{name}_myokit_comparison")) + + error = np.sqrt(np.mean((gen_mk_IKr - mk_IKr)**2)) + self.assertLess(error, 1e-2) + + +if __name__ == "__main__": + logging.getLogger().setLevel(logging.DEBUG) + unittest.main() From 79f72958dbf35aeba8147059c8106bf01b8f7dac Mon Sep 17 00:00:00 2001 From: Joseph Date: Wed, 3 May 2023 13:53:42 +0100 Subject: [PATCH 12/27] Add myokit files for testing --- .../model_00_myokit_comparison.png | Bin 0 -> 192672 bytes .../model_01_myokit_comparison.png | Bin 0 -> 174954 bytes .../model_02_myokit_comparison.png | Bin 0 -> 179771 bytes tests/models_myokit/model-0.mmt | 54 ++++++++ tests/models_myokit/model-1.mmt | 44 ++++++ tests/models_myokit/model-10.mmt | 64 +++++++++ tests/models_myokit/model-11.mmt | 82 +++++++++++ tests/models_myokit/model-12.mmt | 77 +++++++++++ tests/models_myokit/model-13.mmt | 88 ++++++++++++ tests/models_myokit/model-14.mmt | 57 ++++++++ tests/models_myokit/model-15.mmt | 69 ++++++++++ tests/models_myokit/model-16.mmt | 66 +++++++++ tests/models_myokit/model-17.mmt | 60 +++++++++ tests/models_myokit/model-18.mmt | 72 ++++++++++ tests/models_myokit/model-19.mmt | 60 +++++++++ tests/models_myokit/model-2.mmt | 52 +++++++ tests/models_myokit/model-20.mmt | 74 ++++++++++ tests/models_myokit/model-21.mmt | 88 ++++++++++++ tests/models_myokit/model-22.mmt | 77 +++++++++++ tests/models_myokit/model-23.mmt | 107 +++++++++++++++ tests/models_myokit/model-24.mmt | 58 ++++++++ tests/models_myokit/model-25.mmt | 76 +++++++++++ tests/models_myokit/model-26.mmt | 58 ++++++++ tests/models_myokit/model-27.mmt | 81 +++++++++++ tests/models_myokit/model-28.mmt | 102 ++++++++++++++ tests/models_myokit/model-29.mmt | 85 ++++++++++++ tests/models_myokit/model-3.mmt | 57 ++++++++ tests/models_myokit/model-30.mmt | 127 ++++++++++++++++++ tests/models_myokit/model-4.mmt | 60 +++++++++ tests/models_myokit/model-5.mmt | 60 +++++++++ tests/models_myokit/model-6.mmt | 66 +++++++++ tests/models_myokit/model-7.mmt | 57 ++++++++ tests/models_myokit/model-8.mmt | 61 +++++++++ tests/models_myokit/model-9.mmt | 61 +++++++++ tests/models_myokit/simplified-staircase.mmt | 35 +++++ 35 files changed, 2235 insertions(+) create mode 100644 tests/models_myokit/comparison_plots/model_00_myokit_comparison.png create mode 100644 tests/models_myokit/comparison_plots/model_01_myokit_comparison.png create mode 100644 tests/models_myokit/comparison_plots/model_02_myokit_comparison.png create mode 100644 tests/models_myokit/model-0.mmt create mode 100644 tests/models_myokit/model-1.mmt create mode 100644 tests/models_myokit/model-10.mmt create mode 100644 tests/models_myokit/model-11.mmt create mode 100644 tests/models_myokit/model-12.mmt create mode 100644 tests/models_myokit/model-13.mmt create mode 100644 tests/models_myokit/model-14.mmt create mode 100644 tests/models_myokit/model-15.mmt create mode 100644 tests/models_myokit/model-16.mmt create mode 100644 tests/models_myokit/model-17.mmt create mode 100644 tests/models_myokit/model-18.mmt create mode 100644 tests/models_myokit/model-19.mmt create mode 100644 tests/models_myokit/model-2.mmt create mode 100644 tests/models_myokit/model-20.mmt create mode 100644 tests/models_myokit/model-21.mmt create mode 100644 tests/models_myokit/model-22.mmt create mode 100644 tests/models_myokit/model-23.mmt create mode 100644 tests/models_myokit/model-24.mmt create mode 100644 tests/models_myokit/model-25.mmt create mode 100644 tests/models_myokit/model-26.mmt create mode 100644 tests/models_myokit/model-27.mmt create mode 100644 tests/models_myokit/model-28.mmt create mode 100644 tests/models_myokit/model-29.mmt create mode 100644 tests/models_myokit/model-3.mmt create mode 100644 tests/models_myokit/model-30.mmt create mode 100644 tests/models_myokit/model-4.mmt create mode 100644 tests/models_myokit/model-5.mmt create mode 100644 tests/models_myokit/model-6.mmt create mode 100644 tests/models_myokit/model-7.mmt create mode 100644 tests/models_myokit/model-8.mmt create mode 100644 tests/models_myokit/model-9.mmt create mode 100644 tests/models_myokit/simplified-staircase.mmt diff --git a/tests/models_myokit/comparison_plots/model_00_myokit_comparison.png b/tests/models_myokit/comparison_plots/model_00_myokit_comparison.png new file mode 100644 index 0000000000000000000000000000000000000000..3fd3dbf7c486e21a082614fa9806a36cb5904c8e GIT binary patch literal 192672 zcmeFZXH=AD7X>=TL=%l7c95cIET|Ag1gRP+!GfZQA}tXSrAqIhpE05|1r?QML<9sG z5$VkWjEXQKO}aADn-Di?7@mqJTyYBsW*S%qVMRCgeKJR(X*=O&4o{RgnHRsF{ zp2c7==4kz_c7VZ{JCDJb<^Syr{N|jFks$tWo1^*>$Ah*f99@33H)HJm)$z2It)tbc zW6Pb*>>W+4k&@Z4{G_AfX$J*qY3sjVAZ2TBF8w-Ejg1fa?)1<44h)9WO8W0K z*BDi2#xw>)OKpd)Yj}5q+g06`-1bQm*)1ijLVSHgb<}*n-ZAZB+;aI<8VinX-}}>X zdEtO624W zEsZXnX#c%e>Yb;GW%6})Qm!R$*8QxqN|PvEE&Sk;zWBKq9MXf!`Qqoyr8UdH{^IY9 zGhff2_tlrb*#6ab-~P{gmY@Clw=ZA7xV-(VpZ@0;KQEuQ;(vaz>hc*Z_ZL6&f8PG> z7k|aNt@`&e{#lTJSL2_x_%}BG-6?-Z!@qgs@1?*W{#}iKSL3f0__xA-84CZd#=ooa z*9!cLMZXM%e^=w*)%a@#{_Wbn426GJMGYw zXKde^sGqXh*rO+PAj_@4^-5@HXyL_w0^C?zn9RZtL#yf8O@wX-D}!)d26wv6yN_Ubqya=7;59zrTDH*P3iomu4;EGr`rr zMMU>sKf%AYMU{WEJp2F7m;Zl<+2uda|1QC0&mVsHA=$bvZSA^sFNfAL7?)k+btUcX z?Io`5GV@r-c(msqetP@)x-2J4t9K7q%$hyh^!e{|lAXG$tb)LrUkDeLwJ7-Yg6QH4 zqRvMu67>)6|Mv8^U-0Mqe@(M#5RUSha$C4~@dLru&`w% zw`BclQ!ULsHZbn**G`-H-z#J+7rNx;Ax7|z492^F;Nba!E1s;l=l;bfzJH(YIr{nCBeUH_$-0t@QZ1_k zpU+j?U{QI05iZ70SMplScheZ}-{>4TFu$;{Fc}ebB;Cen;gTihN@D<*4E@#3v6gK_#~aY;!4?u{j=DI&Q0`i6_o zo;@2drBBADiPKLxv3SkbwN?S`mUHe_=^zBx&F$$GatWr zA`jWScI~Mxb5JWQ;aytFJgYECe6BTmnd$e0) z<4I(|`|5BX~ z!|%AG9>cvxAJ{&(CD*Qf;o5UvEPA10 z^p2XVSFaxT_$~2d<^9SPvimrs)4;-Rgh8P5$u%B>zUw{*ohE zj;9h0GdFiuJy^E))|QJ>Mp^JiMoWQZh9JgeA1a03f78Pr6ZM7C&6Xp5YTh%8m z@o1{#J?#i#+GEoB6PtzQY(yHeU4lv@6#E8iQ8Tto=IxF-P=;F8j^B%T4U}C|b^T_f zJv!)E>+74`!X+jJ*ur=0I7Rdt`N4{V(b)VoVtuy{8ybrA<@3BfKRumESNQJXuRA}s zv=ma(+h~yH@7z;cxQQilu*E_2RUeW|nze!F=a*|XnU{w17tJ7D>htH%cXXmP%CM8# z8*@GIZZnV9UI%n_7ux0y>TNjkaCsi?F;4Hn(JiOmt_oBhddGj+N8Rp7jj*lOm)0G; zS8zexdudlSGYD6yfhE6lm}#sNs~s^bcL_s8@{5YOU2o6L%|F;=b=us+ivdgq0?*~&h&AS~6XO}+WNsT*>(jZ*A* zaSHlGAnMbGO!KM6!+Sqhg8pUScVw)8XLIy$Ic5y+uFiY ziX9|7EqA3g^_vK7I;K|dSiVE9HE{DYeR1ytYv#*g8Wnye> z9PTqU$_;CpSgozZEl^hJD`mQAD&{d=(h`P8d-DXY7XIg(*>C&KqHq~SmS@;;MxF|I z&y#C(HQi!UpRtFXo!!So!3qeo$yz6dEke5^zRRK_PLMZb%Udkh2p6XS(|06&Bp&d+rpyrdPfx@W7H|GT_%mD1i^zxH>@rk;3k2`+o$d?t`6+IIJMHym$FtIrXMIFZvfNUaoiRP4niDAz0A9x|UFr z$Yk&KKQ6o)m6T}o8gB6a)LU_IUxTTVJo3RvZFSz_Y9p5h8_rYTo@{IFaQRKzian2) zDmb2u@yb(~`1E|FNG|URSGZJpxZZWJ^6>7$b2I%yO}w|&rrV19jFm?X2@7bbQ6er1 zShgg+-hOvXsnhbDHh~!r4AKL4oSRWN!sQH)wWR}JxHey0a|h3c_G#Xw)%7kst~gJy zr!GB6Q_SsnI`@}?q2Z6>yb%_mn&$Y}$lA8?k_bgPZTSzsWnkw!3khgg@v@)$Oji4- z2Z<@a<(V=Vx2?X+HP>o#+#Ec|>zV?SO ztnhb@OQemIF*Vjzu4qX>39x7hfQMgnkcxYdqQbRGif}xHwTrf4H=IQ7ckHU#So`=? zNM1Ew%bDnBH!;^f%0Xu0or5uueD{y*pUbyhRQBNUu2p~7eqeo8V* zv(h&_{q&pcfW9opj``Sm`;aiOf6zWhBHNW~E3nw1 z_4UPmwyJ7Xl2MqQlC-pRpki+Z6^(Gmw^~$Z4wwk6FF*S^Qgvzfgp>5OUFyhB(by(s z(V9W>QzM*clr5n)fQJHGpYgg3ROtl!3j53tzb|XPWtp^YaGMR%cKGRM-<9FX9`K3m zs(RozTRP`QG2~SHCMJ^m*icUxH}@VgS{Y)G7ECumZzE>{kyrT562pfq-YA+1t$(@* zhVIQ@pZ{LCo|zMDmuQmjy~$&!n=0h%8=F?DjI~D}!2QN^;59eR%tW41S_19kw-CC}VJ2ib?y3PK04#HAUkQ*&N zxfXdL_eh?nTh4fIHvbA!Ee=hLvV8)v{bQK=xE%Tr$?|@2=XbwuI{xA^ZL?z?Z)3Q0 zji?8v2&heXK6GoILpA1`>_5JswsRjOHD4;Gode{xm_0DzjGiHV6L_?Fk$ zdnaDs5G)Rr4q9U3x#{TRlUG!H(CutO=jBS_nVvsst)@oX!Wb9Ne@~>c;)&95q4kH8 z{Bgn3KpJnW5t4FLh8jvAaio4jo7Yl%{xrHaF5bT}mmci*OMm<&v~BU?#YJbPP2cRS zn_!gfOiO=#qrqZvp9wp9B-c0S&rdOD2Lc9}8S=0*3IgOBkF?edotPNmmZjS^jeIVa zE@%BNATDNF;Pi8;R{=&`ixNkFmI9vg@3P4oZM}hLk&5aM-Yun-U#;NO zd6oMvvN6=9!KvCiKN69%)h2ha%HLE}(b=D1k6x0Y5VrrN;U*tdz_CSGeZSU{uvj)V zz=n2(i@1o7Yh982;)edr?`Wz;A=YIn0Z*W)8!BGz^c?jPj|f(q_;_h_t0EK*=xkv; z^cbH?Kk;`7e7Daa&Bl*ieZ4W$|B{pLsx4SKKVi9Py06Xn(U}_1;?x_l9lV=%kGu zTRggJv+{U%3iq`P`)+k+X67dEiII_^^n9bvdh8yBDK>9Z)vhx>wg?+`-t5I2-Z4Ew zSreoEE3fS;JSU*CG}PGrR6MKy^XFqPF0bIjSEMzoy0ncFX)_Ji>m{yD9e)@2`}acg z^LdkfGR@D8wqYN0#=4k6e$sjJrScsAfi7`ncgKldmwZ2s0O20i`b1Z>U3T?$tirom zA-j4qLU+Y9cc-wmHC0Rse?fT-s!1_tA%WfT;Z1P?&#%|+kV1a*M(xvha5R-Jq=Q2N zmyncY5t?G?$t~pinM_FJni_d~d)M0+Z#1io05q}dVDB#&I^)QachN-b7r*OZB9ghC z2~3SpDI=^Z66;b`Z_n?=Rg`S9d(sBD$e=5way~sHbT_O}7g_Jp{-{F|~boIG^Xm>#I5s}+s5!nYyYD56|L&eW?p6Bcg(+;6FhOrz(My*nXQtT!i)a8HHct*W z1^N2kxpQafmXrH4oo7d?`gl{*WbFFke0$XcgLu|5ltlvCM1WuuN2*Mr$12(fknuT3 zjUj*hix&39Eq3+PdhKO{7HO=XbE1ctBD4 zP+9^r&pnj;K9G-#e*9&-g4H@HsX!Ea1cD|1^D|pdqO0*2Q&tdtj_SQ`YCNBpz02kJ z20vz_qT&TzMT_brqwj?_ugVRn#hr3c_b(!#il~KAp5c9N#+*x6R%#XjNG}u=r0l^~ zB5am?i8fWvXj>GPoqAeY8gwiTT~^sL$ln>0Lv`g{lr8W~7o*LXQc-O? zdkhU#GAihRJKFZ_fSE2aUgJa8HtDH}9O)aKZsfMIoawyS^Wzmw2?+vdmq1uuAN9~T z^T4l;jec(Fa6+;|{fS(=NT4?k(M3DGpgM0T_1)uB2i4SG15d2h4qu8--PoVY@)t^> zheU~~-iJ5I1+I?iO93M$|M|trVoz~~ysjC~XU#uu*j0^muP@g`?O`5_*R6M|HZsc^ zy5lExNle-OZZ#sv7pwUA^>AMM9ltxcTQl)A5vm&t)@bixF57Z)9(~o#q{$Ln?}ITP ze#ZW7{(YIv+VadLMvgi~V^O>b^YiR-)Dr)bCp#Oy-Te2@o132kYbhaTUiEeTQR3|%m4{zF z@;iQL+ndIIEt%8rpA4wl);#(FMfpnMsm&D7LS1NdPBF7zX{bF$2R${+n`#nbQRB0@ zDIFY98Pg;`H}WvrngdS<&>9=Gjx(P8@dcI2(9NO(b?MTjw_OMKDa!)c9&YWrFjw?e zawPioo6Je{SNmHXiY>7oLfMkusFFc=2-iDm@zT&!!x@+CvB!&0G>z{$ppJ$o8w1$< zu;;O|(&($RO9kgF{qgpmwjPl1|@5!OO5H17k6v#^7q5NyJ3SUKRt5aSUozqxZSyCci?UGRTC(zNoAQLZKL4hLS|eSZ`-&6c()jX7 z-_rfp)ZN%k;(#g}`!m5oy*v_Xz9thzcd4v}E)dOf1;_Sznj&^;;8)ZUD}$D-2rde` zeyfM>wGeHEcv_YdpraMg5oGMEPQuLuwe@Ap{($Uo43H;V*#$UKXp=Ee$?(#O6a5ZZ z;7{5s6AtGLyxrYlN*NqE{rdW?DoY<*)odoBW?-x*y(!SvyFYAZfmHa+W+e7>Z2ucu zP8~3WX2Q?MzzNZ4&2$j{2Tb}2sry?BTxEI?;IV4RlFjW!nyLVZeqAl@(VbkD>t5%j>S$WL!WNn zkf4OE7@XV-WS*HM-8VXJ+V<(#K7>_Ba^vGu<;G0?jMKjkbX7~^F)yHviau%qAB_QB z%Q#2XN61=MU3JU5Z(ks&d9@=4@z+=Gl%3^ln_oSB#`iD2eK>`bN4477trxMp2(;8L zCR*s`28|u43PRa#SUPFb3ojm_qg=4L%B(@ZmlLEKXzR&YCX}5V;Z&thk*AJUuUPi^ z(_e%kd7JspN&y9Il-he!>@uCKoCQ&3Bf3O0D zxHgyZ%=)skpPg)d@-;QP$J_EB<;zWDe8HyPo|R~vt6-WdtuD|w3?e|00 zCO&M*{$Wf(NaOWajeUUhxd6%CpX)D(xvw?#sxb@g8FF)(8|mEmL(&H4DZ{>r%%+I^ z$sSvs{rfK)`?M)GT||>DwNhC7miW}C3wFgD)72!ey^=THzki;9(x>gjHu(dCs$*H3 z*Qkza?YXt(RFtYgi{0|PHUV8=8b1IUqGyIiw{5j;G(5YYC;QTfSS~Y)x09PYopI>F z|41Ki2xW4|n58d}ym;0u%%_5SWvb~U$0}&YqBtapB6loyi+5lIbK5-UHFa*Rrx7dtbBSzH4-8pET?*09xhg~NOgUFH z#cE6whpyG8dQAxV^Vh7frX>Cu}x0vd#a^nK}N$Qstwbq9}4O^n8rP6E0$ zhdY#B1iQ69DYNy~hCRCrh8hkq^QXrAO+iowfJ*CeYAX#=8H?MjS!b>_DK>{;75e2Y z%$bjtAn%IYS??1QopbI1s?1aie}mq`e-PAuv6Gh@TXh~Xo@s|q5|k?d<=JMJX0DXT z0}cC1cs(~MqxsUG^^!OOeDA4pZ#IK=kl3;PH?4Zkww?Jk9sQTh`LYfl)Qyp#Nx9QK0CuCf?MrjPrqw z>(*U_z`PaQ?e)S_38jIG*RNLt7_UDa%rpP|@nte_&A=F1i$FqYvSR-7r}JA66%4m~ zO-@W_T-*00F_DGF_e7#OJJKkrtBoBce*XOV4kwV@W+s?%FioID7|btYa`xajtjj$`6NUe7|(w z9N{hRIQ=nTwy*)*bCIRj;-`9#qE8a7g%

0&>1nljG*EPnT|ME5OrzC964qi)Gb@ z8FNLJ?o;BbQF|<;xam4*8{o^3!jAj2hB8&?4kY_i(U8Tp+et|5DhxYp>}Gb886jh? ztsZo6#n7RZY^mbIKcOh3$=az>>Nj zN`iFoZJK7NnLS4vJdat~OG~r&5$*Us5*RNbRBEltkJzI(yq&##)J!pa8_>_|^x&SU1 z=4)NXhXxxHk*$aU6dD$SKys6*x<6lYI-~B3R_@HMYBZSbP5D0VMi9C4Ct9S79eWn+ zQxYU8%H4?6J3^ozWP1J3viZMRC(#_MI(cfM%Y;|UaUwMeqR5@inhEnkGn5<&bvpdm z2OZON33{lb9tV3Lw4u=1TqG+WF0<~Cnudd z@Ya8jxU73d-!6d!<}>nKd~@KnwywGVTL-pZ%V@o2DY2%FH`(3{;+5o=297LJfaD@s zL-&pt9W>gzZrYg`I+Vx?T|svx)12KjX#+-LNw$&?9*iAzwy$r>$T$ZnE=@gICJ_jp zbHKgGU&y0j(2{sb^cb(idM+oM`W^)tNDEosMUDz+BL`aBEQqQTRlwE&hNGfd$8pLZ zZ|G}A9p24!u(#)Pwf$Z$KY*`6$1i@^#LH>yh>niVKwpz;B$&xb_udFmz&lcJODnyRnfwb;K5{=x3U`6az}7YczM z3tVuy{t<%N;6ET$#XsIQ?ab!|Ux?+0NT7(hb7z@N#)`4wK5=o+PdkdZfQ>AGej#O# zA%C>I_kzS7UT^XUn454kWcA)%J$XRm(uH;`(Iq71anK-ID)A=zxi3}*6)ZYzD993Z zG28=bMMA>#Lh8B7(zCTBHKjLsdlcE65)_Q`8qz1$M_YAr;7G2!6X6zQltbr3R+odv zSb~gN*J8@g;^-G!-1F*evNv+rW`muTcCXgO7C^r_BLt{CyJ)aVSp^N__}LB;f?Q1PeKBm5K;lkquvY?`PLZBE-qM>Q*4q3AkkC}t` z+PH|T-feIkB>#m%G%aG_ft(fI)I#2z_#B$=-!#!GSB4Z6hgJ8)dL$oz8GNAHr$4Nz zzGcI>X+@k4IpqQ`AfE~o4;~G}%?tstoF;Y@EKjsg>mfjO@D=EsvCae}4AwAriinG` zE|Y>G(d9_B!T=4#;gNKacqI8wg%cpK9a=W#PbcJ;8amgQZE_!Izp>5Xgqjs^QXE1? zJp#B43LL>!yhCE~h9j5r?g5-^!agSF!+CZHz8|ZgX$U5;!<5*5l1$4U{`#CWH|j@; zRX6b-H5&_G8KktCWIIw!i#?f4=l(D>wrXlb0i(dV2z6k!NVCpiGuKNr1DghYcz(hB zw=u8`+WBcfRYlZiVP)}TWs*u>zPvczm#>N^`q)HSTi~Xo20@^v! znH;k_vwv3;LG>2PR{_*P z9Av8~n*eW^zPh@OV4gn*&MhvNCC;3p>w}GLu|B`g z6+LRfdh8_&_Rq-M=>9Its>F-9;vexzi-|0B5fU(_N9f%BZo{5XWhA{I*aYrZ1pWs_ z>|&-v+r9iE&lq&K+*j+F$GhhrDI*+7?Z({$@%|iw%V4PG;E1^B9gx-@ zKw#p@s^<6Kf4efRFL)5Xk)kJOr|&**WugUl8!_30yP!SBH-8uM03sq$CnrR*k?VIC z?rotq?}nm_4TJ;Hn4g}BlD-E8av3@tcO~fNFcN07hwH--en|iISByh}k&*JmhVkXz zsxH7>{FGL-%(j5-J(NOPDI5&)5X^?gGhLiVTS8e>!V7@>fGnzQ$GeGzmq^1lPBOk)gQ&;O|W?dSet0{2-(Vg zpPmZ@LFzl<5l;%Gi6>Xxsy<`WT-AwB@}v?0v?`x`UOgS2vGvD#>oYg`OihraXo@|ogU0zJ zx>3^m$c{j+ur!v^FW5>e$cjv@7M5`}n0`@GT7V{%E%W#B`c-)gh)E;}M!dQ_B7^)n zqgp5zYn>>4r7<ckgj zpAqv9lyCMWE|Dxjbw+0;b(3a_xh^;LLin1Wn$KXo-@fu6XDTC(0UFLsB|RfU9wgO5 z)P%j_VCDCJ{5Vpk{i=<&ZfumdLN0lv!5(a@S%yJ-`P(Dud;b2)GqX+uTNi;)gn{a0 zP+Liqs%oiH|7$XrjU9o1h-CC@$`aqDfUwK>{A#@^KC8o2l{dck3r`-z{| zj7tH1zqqoJQ`a=*oTiDl&i&E(f1ZNh5^26fZR|Bz$ybQO3J6$F7c&@<@5>%y`J^EO1JEC=6*#!F{OvX3rG{Zt%+kMQ(qEk{}O+-VE*@ZthhT zaFC@a;}nzAmy?_R*I@tu+R5O!FX+m+U9%T&ytMD0wr}vAfU?A{bWWE^*A2RH{M3R@u{$IFx^}Vx|F9D#(H(#*|6`=bht=*%|G#{do9?`hzAGjZp%Mj-V49&T_&d%dq+PeP%bc)J5|w^;Jut)#E3io zH;eNA^8UxCCT+y)az5rthW{y@-T(dX;FI6-QMYc=^oP3pHgo@ctHc*?y&RHqK|JKI zZ~S5a@BQV=(8i~B{KffyaTv&JiFWn?@AYON@ znjZu6$KxLHHo%A!&YggBDeg_1>g5)AoO6a?JFyCKeDjVw{J^EJu|_QRCo*AjoE8BW zd#qx}FaLX$8N1TVJpMcqC+_qYmmm6Q`U?IDXD9@Pbi#vJ2X0tI*a{5W~UzarD)7{{5Bit-g^^-?h?C?)N3&?K&>W`&g zy)FDoWc&6fJGbdw`etE}#MHtMXJ*Y3S+wYzuhl>Oyy)wJ{nLLFUjE&xv)h-ixVGhh zOUOoEMjFi~pN_`&eZpDy?H~ShHj5k`UR-9}Wo5bGEcr$v(8W4>dW(q< zAzv$bZuIi`iY)qnKf=c-op%xwHxMoa=d#aVR57iRDTg<>q9B-pfzUa6R2XLUK+utt zq6}IldL~_|91{(;!6<9XR;zdOND-&z56-%-h7|B;luutHM?nsU3Qj2rkZRLhrB6@D zF_%zXGoSotAm<9GcEj^Pm$r2C@t;MIj33zV13$bPO@Rm))JrG}5iiH8V3a82V&Ex-6g% z1Dq@Z5muXtcu0$ylcPl082q9iD5~{CkS^A)ZNHHl2!y+U3NVq; zWEMh1laf#blk$S>?nk>5O!$l($?%Zz?qp>t5Ul}eJb-QOOOVk$2Z?xU0HpV zX-&4d&xCb_Al2r!!KSrmu44DYk0AeTD*7{Hpy#Y6i!{>Ofws|^-!HvP3~RiX04{5H zN}wvPu_Vi8=T=qT9~Z=qD)NAd*PszOo!gWG>T(l=%ny2P#m6_lBr~IXeqT&X%*%Ap z&36IfpJq9OaHv&=?a0a4(g4n^D}#mtu!qFQ$DcaQhU+w#>DjVg3k_put^l6;6=FtS z!UeU=XtzUtC_U}Bi0{N)laXgj(#9Op{usV9W^uB+x^-CH(zWH<;b2@ta%1(=Rudq> zU3p`b%+;h5411D4w<^-WCNt(~S6i1dPa4wT$&kkEt<^`_>RupsZcGfss_vzwxcm6= z4a62A$X zHY^8kRfj}iQ6ejhYuZEWFBs(;W?BN6#)NZMRqSnE%=@zGTyC@ z$j=pL1LJV0%OR)7CUnUfJ=-8HHitw_94B&(lrhr*#R4U`7%4;qk~p^Fg2?=~a0l&k?8)mSAt7>=o9P~O4;`w^Jb|{u4qXC^ zj1+b8o0dpreheO|h*ep1@ZO&Ie%IhtcHq4Vr_B_ZN5pu~1n!BG!3p#OOpa~r^FTV8 z=r&jN3j$|5$D2t!di9hG(h_MWd%Z=7RDKlM1&N`5j&pw&+5uBYk(W#+I}g(kiF-8( zB2yp$i;}W?h4T?Jw+s)PLU!=7gSV#$!M$H`_?-B}a+>7G>C&w^qB{C=-!V{LK_o&@ zw859Sgm2|R?`D1Uuu6UZXfIk#KWc$}rrLapQ2(PJJ$m$Xkf;I*5@N(?$JWKaYDFyY zlm1-TTFAmG>}eZ?{@N1^mOm&185;UIVbYI2JrI$gNOK-9j(1ME6=42@58i3xQ3~m{ zoVH!# z;M^MAF$r!~k$g%tme6^W@-ysG(om9_)TdW*$udLbC-L^rrV(JP?8bZ2s{_Z8%Zc!T zF0V%IQq{fkjuB+3rQ}d4W*t0uZ~^%=uP?ekTFnuHAcEE_q14Oi8&N8ll#?mTktY&* z2v5(l7dJ}Y=&m^jAD?fU>~*0Yz5s(uDXnPrgs}o1BB!5BOD@+|RNt}VsouCYlIQO} zT7^6}EMPPF08#C%odTeet{g(tQnl=vfR0D7{8UAzK^@uaMlfGus0iiT#a~=ir3j)0 z`k&&X;K<*S^)MPWWAn{D{F#L}I|VAmnlNc}48cerYnt$+jy?~}Y+`yKynl+XUcItl zk^M5v(}ZR|Qh5<@RL(XT`6<^DMZSorJ5$s#{eFBQb5hKDbnHQ5qCeAnFoCwWTzpc| z$RN{4l*FkX%XqY3-cVc$q2pC`R$EBZ00&7h4V~lf&_usO7#(XaxJZjebal_TkH?Ep z50$1vhYnc}M(RRRBDjJ_v}UgS`!nyz7wBG310vXcTZ!r&+K^r*xqdz8hH2nM5bM?n zrpSIfyAc7&#ZuX7&=#wbDIDg)&r38(2y+yxP)H+V7nV%=Edw@bOh; zPEp?#!ko?sTkt2eLL|O9JC9LYLcL{Q-ndOq2H02f1Rr^PQipUQ@@N(kxk!8k)eN#? zQ#S?yK3}UkJ}%A|&9O>CbUoVaX0lzA!cQ9xp+(K&MvUNy@-_WuzkLS_7wo>MU!7aC zH(2+k)!F#*z$#z#_1)*$ynqWl=9OxnwC(74w7u&@`||WJ+LDa zR3hc68wmobmSM}ah z_HgEViK!K)r0(m>k8eOJBH;O`m;Q`~P96p5QIZR6^DM-yp7}Gu5ish*7vX>8(Q&Ok5el-ziW5g^2#P}Arx8{)N$qVdSKY-A8&MZv9_e<1 zKsE-ICSa)yA~O(dI?WUpH1Gk2W)VxoVNX>aRk(oOufAtE?PAitY@-cI=L}K^_-lXHe((nk=fbVX&9;=>Hhxg%ttCIqy$lk zud6?(t9vyS_!2oU`>7WZXH{dIKt(BnmghW3!VE_l(;rU)UU#o1#Ar~tNbn_0){`@s z`1SY>^Wd;6jy(}UWOKdn9~7`XRlwgp2k*y=tB_Bj{=ewDE;@yD&>FNjJyMCRd| zyTTB8BG4RXIQ#6_vBMrQW$Vxc5-WJx^cfYs>({N*0-JNRZ;ZB2LaE?NjaBG;?z>D9 zP)FBtr?hY9L2*;l3a1Fl|60hExgM~5sag%*AhlIZ%^vdzvaY$5&n^4l#oymoOL|j6 zmN2|erLSr|z`R#xN`<^!1@iN7RV}pjH7M#y!}U&HK_pS#g-X<@9FS>}aQ}WV?QBaU zGzTFVY4|lP$fdZt7k5+pevIPD(#yaD18UeChzVJDTOgUwaCcUCOk6jxd8sfA0psoa z(a)oi`+~1~*&h9DDB|A!#~l86Q;3TBP3_nan$pv+rmB%!;u)44i6e-cq6>Q7j}W<@ASV zG%#VoB1}L=s4UPL@AKZieLGSfHFdQzz8ga&Zh9B8thb@abHLQvo-UEf?w(HupGU32 zoxfZ!X34+{lL6qdELTOM`N?#F*PW6X8_`v{=&i|zk}4NpXXsFW`uQ9h8>|jw(oij= zoct0mf)m&Usm_LgyadZ4EUOdsLsGYb9E4)L2=*D%+&1Zr8v{r!0RVZ}GqY?4Of!(@ zIpWb785xbT<`GVM-6##I9HMoo5`(^Q((rwOw3EL`Y0!fpAzcO3nErSoxr z-x1Wl#gKXHT1}6noY+ea?Nko+qsZnaZMEd9G4`1pthIxT#X%=X)N-}Q6xJjIUttOX zE7gJI0J^%c)Jk&z1^$G#b$WxIq8!FJU9m%&SBj7RoWc&1rj(OtcTnLK8IVu`46^SY z9d%^4rQjQFD!o!xpvU4Sm`qBdW*t-{_35$tsS;-dw|elZR1QG^toNfC!`;Z5T$Cp% z%ieooHe_qxI+jIhKUO;Z>F_A}S(#JAJ$N+E-5hZBrOMzxE5vu<_b|X5Mo1 zs*QmoqZx*WlP9y^r~0a_v_JX`eG>=w&J-Q)ZBWeXkLcTC(MKkPkrrv6ns$M99)`Cr z-MS^PGlaYY(dfE#JP6!1Veqo~KeMHLeSK!A6Z@cX!Ihj#vK9?!w}h$Y^^73(K5(&? z5@(aJ6t8<<$I^g5$45QE*7^~U1ck9GgYjE4_0WURxIrq1z_|Jnql)&)sIiLJW6-0C z_fS{kk4_xMM3G(>-$je|!jeT8j4*PYk`M4K#{*??>p3m!4;WnygcxM)&_Ll%)J=~( z!b$*5hIMZperHPX5Y|ZzwE~m~I-z6m_g~q7KO?wZVL?Gq2zIc^+w^;JAYQTE8>yjm zfBtLOOa~bFFbrMr>BIB%QF-8uMO3=Us`(%SrULIJ{H8b4Wj7TE&(4EU_j98pB*=!# zpV}7Q68EyLSToE=RpQRmhsGYo=^ox5XyWI`hTHF=9v&cwcplsyftPJb-BCL@*Z{H4 z;XtWS`%YCfk`12k6r^Jv z*i0A3Nv|(*MkE-wH$&u;Ox(`XS`$&xF69KjcmSo{fV?VnBZlKm#D_%(D%|Iu*C+rG zijbuNUg1<*n2|1l7TRjjgAWi=R~H{Tdjt;YTuX;RSCUW2S{{;AX-H}X0=_CxPT*i_ z8nIx(KkT+(k<&C}`ra_rSTcOUv1OkiOT?#e6mUU+3NN2Go^ReD-m?NvjxH@O6_;sV z<%(&#WE2iLi@sK)Z)9%n*>Lo6Tw^LV&3YMK*c}!lOU@Nk{qcmd?qHG^tZOMnQRDS( z!b-!fB)TRpufRy@p=+mUMG?LJKS*!64~l?A)JpfGVyS?}d^1+<#-`oE_81b|>fQt7WHjD% zbIRG)+VBzwxov>Fa6h4r_8ggI>O+LN`c zAJV(Gr@2+f|4w*6AWs8me*le;v!RJgp`mKzb!YK#P>*D6u*2jv3n3svM+(5Pf|k+1 zYok%tV`meivpH3d%+5So*Ocya>b1iduzcl|i{v#zOOMOi`ly#X8Y3@*LT{*)g*aRFP#{sX|vGW#hC4r=^O`zAG^_@zB{5 zsjm<-LE_Yf%uRU3_Q+ca8qSyB9w2M1p{5ra{6);us8LJV*F}^?Rb!f`t?ea3d6xAR z4<{f1-ISd=EceWsW}#60<;E1kW_GRYbLX0+8$Q-$g&jB zwDIm_fDw^FiF3_pzxwSwO%AK*D(lf5JnIg4GpujWc{)x>FIzY@IS#r9b1|kt4N0n& z(ao*t@D`jf@i0jHe%z~7Ch8P-Te_z}7kC{#Z#ZrzbRW_Sh8sH~AOIg_Df6Ez;2xU35Xj{%6>} zlFQEKRb6Hdd*wP_Y%idP=P`uF1ilwrjL(^WV6>L2qkImF(!&1y5emd09UJDap2_=@Vsy$T1t)a%Na4+sF z0En;e{`8AlU0%AJx@FFkoKKK;Y z*q5@73jufe#xcIj0tMa3+@)ld7WK`W;7}xPMuZxf3Ha??x{)DzOn%faRpm_(NA{?d zh$9-eA?AiWWrSIzjIxe9y`zylzEjpCG<&c0;Fzt=Cox9tt$%p2bz26)#twe<+s7BL zDicg=H!dk`Pi3=epXDW649jp#W5dI`9^LjDtl!l9pK5!f!Xai`U}n-U34!Mfyj^_^ zuvM$$oD629l7E6`*KRpJ*x0hUqP$FujprA>$Y0_8nd!}d1UtJyK-AUO8Sfa|7*!sH zB1}6Ztv#II!MYazqJ@*5LfHDo#*KGwrLBPhz}?WXOg(j4F1#Uj35*YaZ_tYO+i-kH zq^wP8e{M@8Qq&fUra5x;^JuDr5vRfr)vz9@8b{hV#-*U zcm#5ffJjb1$|KT(*kT%i@C3jr*-@ybMvX8vL=}O4INbr$A%P%be8Td-e-EMPRM+QQ zKZG)U7%#s0yV~3?-z_^&>bJDuyR;7_mr@1I15Iz zfM(%|H4sH(#~tAWv7l4kGTf~x14q_UMQ1BIWa7#yE&j<$WO=1IC?c(m(qSBsN8`)E zqSZ(aq$*ozUC!FJqos(B&q?uK!3J%%m|DWjMXOrxtbj^ebcB;j7;mA5;Hh zxa6;d$GB{%_3ga3Q;lw2gE#ieY_v>Teod`U-{+N%I>Xu$-cWrGU|nU(BOCBEwiQ_z zFEIZ4>h~o5%j)1wkbgoE(!>FeuPxuV(yG-iV$tCCtHA+bK( z0FIir$hwih4rn$Oe!)i*W0*Y9@abtQEz}4W*FEr)Vr;qa_eMT5_qUHv#(5sr5&fKm z*mT%Xyh%%*A_NbBy=A#BEw8S*Ez5a%jU4reA4?Xd z1*xHJC253_*DBjuYMwT)twHXxmKCp#r@N^=lvygBcc>?b&uN{eHF1k z^wvG66<@Wt4>kGl5+9c0K#_T9p?5_Ok)?a@EDi7hQ7?o|(6I(Ih6O3S;?!*JPId$C z=^2NvRDRbIwR5ATp^i>jd+tn;)7$kJtu!-(BLn0*#sVl7%mT<&z16)BW0_*Wwhs(- zSm|y{p26n_^+vH72sS`KW#qQ=R0?`)y#;3OH-0~&eK z#AcFB87sirV%HxoMd2Z!)ku&`}U9gwRf%}AdEoXuF?1o-yX9*i)S z6kL9H@`rWe4bv1tL3@wXx2c56;^DUyf>{!-Pe(T-kNzmk@d1rgg2UEToI;5dqGPQ> zO0Vq00tzX)p03JB7JUv;v8ez>5HIFu!m#;|B>^V!2p&3zuhJs>myVH9!(HLD7E0cL zE45tgvWB?j-NfyJxVEZ1nzn>gTD_5e=;Zz=7n>|gLh}YJnffcDsohMe75p#VfJj~VK)g%%cx*K@pUi?5%3TI8JD)xK3ASYslk{Di3? z8gKoxRaK8iJ;tHcEL)V(2*Wj%r@eVqt;4NA+F?2Q$uNQZsaMEr1&>sv=T*k#V#m4= zAW!9>pI|{(ZtRr(YP;px*uR-E(t3b^D1%*=c2M)bbLq zWd~-l+ARlD&~_Lntqu%Xb&)nUQQPipWG^B9TVMn)b8=iw)5vgbBsfo?UKW9;)K`o| zqop@3@pn2vZxcqdYwa*pNr#`#X`Bh{#h+flET3f?Wq=XLmv!x3RWV@6h>$3! z4j3$FPa@1^{Jn7OL+UB>T6BSC0W|Wex%wba`UA--)QR(R(o_>B_uFp7rxr4E2}HX4}w*!a}FtR zaG^9iM2#ELe52D0eS!nnr5A~COom=nNM^hsrO!Xhh@I-{`Zg#Qrb@fFT5Bo;%R3xB zaXg4W35Nzww>(LKAwxU~Nl4VDIPix~2hcFFBoYXiygHDMha&jIClt{gu@ETN2UBA( zyFr}cwM7+fpuvw}+{>5bdq@TQ+ni`*iK3l=2Aw}dp5S9RBz3TAl+G_A^SmOR_OS<4 zD;;c-*tROLTQf$69;&#L1gp z_(U!t-`gZ;|3))w2vkK94**myXSM{Dl*$>sGRzc*9ykE6`wSh)0h;EeiP4j2^_Mrgal}E z{;_{#fh8jpjeXN-4}KECt$>Hg2}ml;EY}fZ(>ZO0;6N7ASzW%qD3ygIfG#mSy@SXc zvVsxSi$YY{yAB5>(MJ-=rm$1piFgIFUvObLxb57FS)SLO8tIu?kr&8$#y{9lWjK>v z`THooQKV_CD!tjDX-p5#elE?htp%-3!WOi0Q%r;ecuPr0kUJl1=pll>LkGD_kk_3^ zG8BJ^oH{>LG}pe|;bt5(6xZMaS(3mFpO$|+sBjLoN11F(gn?38!+(X866!%<(d>&3&( zf~#)EvoD!7U#^<^VtxmMFc9RjFW2Tb?oM-rK9Y8HNd~~9d2V$dZn7-jhes2lB*MYh zbqtTepeGySfyba0SzQxW_|Jb8$=WtD#Bnp(UQ#%Qg;@yrO6O<5oO~7L`pY}#l0+B1 zH|sU&#CAlG(7M-V;(SkeQZ$L;pk<{~uqbd))Q?v>f}z~H?wtoMA0`j@$9?(wOijK= z>wn*&hvrs@whbG;J;kPfmpSqJ#jIORSSdR+8~k}LoSQ_)8M`8%lf%H$mI^Z-zdwlO z)VLXplLFqezWXb&I=aWN3hA2P?tdWO@G#@D46lzjnV->e-GGjunv;ZsRJqrfoAMXQ-P}D3+rM+=5+cEb@2xba{8`grNV6Ugm50jF*^Gr8NNb=ohL9kXj`z-LR|C% z7vlD6ko-w}i^l{IO*6SyP>oNVpIQ9s75^yGeK>ImrJv^aXk+_zD#ICc4D&c-7O--% z!ntdtV!C$3gdQsm-{$l?i$=(BuXBt!bQmG+d$M%Nn&{{VLIw0kSvNaB6A17ptK79<;6X zMr|eup;Qv!3$VUsSr4ljO~UGJJL~|mlCQmxEZ>a0*CeU*4xUY%Ubd_m^m>WT`OHl) z$kOytG)a^Fbes;c$%v#|i^P<`AvS{sd(nBfrj|5d%6#d!Oi-ffG!iSfWq~50t zjP#@n4K*Er&5{!M%g*}R>$`BV2>bi}Y>?B}3(y0IfaR@`RHECag9GfRgugmdH2-hM z1Tgxy^MYspG;Mp^5FIMUMUpK<65@J(87?nN3y*4I`8lhX_w(hwRMrzu(9F z{=VPq`k(9Rdb)CaKJRhA@A0}{_x;Y@sYUL4kPLjd$TkSJj-%W?D#<`(-LGO`gRY{O zbrNQ6(b|8uP55IU>Jx^A3JBB)thFQw=oyGhgF?NRuoYj6m4PWqSWKhZ=YMOp%MK%N z9~HI2ydrI*fm|qXt09_@I>Q1TMO8 z_5ig4)whCUECY9%5?myNEdv93R9^thI%=%h2bDf%Eio!~pf)e{0`BUA_6Gris6!4= zHo(8?U|zfo=^0c-rWh?MM4%BG>_VI!^r$h1_tQ?6x?}EVC3n*c^ z_M#v{XgYUw&~35sv9K5;W{+#*v`?LGPwNt=II=qZ!9LToe!SA)#-9imwv zj~@*Z(RLyF)XkmY51H^cMr{`X`Ky3-43Ge)hvJA+&`wnyI;Et1{~1(J02DIf=ARNI0)C)9!vm34u|GX(4^2T*`&P#h7XQ@ZIJjL)w<>1#VjM6|qW zboaET8#FNwex^Zb*G@F_8a0SYp9D^ITlzFq$5x=f{a#sM?xAe}jw)Hq!WZMa`ZB?& zf=L|wy29@ZK*0-c)$~g$fd-=3X&Sb3W8!ZUVT>pR)dxt3dF%9)3b10m03hlw_94Z> z!0ykAL$Q0Fez`c>RVYp~&#WVI*XmF9Ssxmp?FV*Tyc_$9=+$L*$rsQiwHkE+>aa#U za4#qfqh&(0Sv4TYbb#31tIuyLmK*x_H0t^iFWU+d9rus<8Z}d>VrFOyD9CL z*M{E#Y8CbH%+m>9M6ftI)N&xI8MWUbqF9MC5f>e`6ygGZQcBc-f&R4tTBbz24XhLDm1ub~1e1N%bDmyNXMKcGrI#HQ2> zf!h_RA8ukJtk05vmT9RyUIt)cJap5Xkzb?VPHt z7N5&;^hZ^{h-hP<4iz@SNsDfqu0x{Md-i?X`52Z7c}LBQ;L212t>iF@&*Ed> z%hJBeMNsJvaMzSvIiXq$I^=BR5s6^uUwt5bJnQgfcRLNJ@$=Epp&Dsx!eK)v5p|eD z#u8e8A%^8YRD$PG9u0LhMjdO>Ee4}+9#BC7z!z92D91~U`px-J32`H$z2sB6&E7BI z=}75#&xP`m(JAiU_kfU{fWwP=Q{L4-5KtXd!MR98NV+W7m$4K5Hx8cwkE^ZRr?^gEUgRu_kL^cI>h3h|dZa$L&ib zwsbx|m)X&i8}Dv=a;Q9a5yt~%)>cXP^T8=5)$-^Y?qeZa+c}(gu~6V{>;s}3LFA@z z3kDI9uJX2u6nOWEAlT+&7Bp=47h8Nu;U|nh|BfYnoLByQE$byox$UQ=&xuYKU2eSW ze!k1^iDVt1{22qbO}mNQm(|;PzKt#}i4t|5*aprYp<@&n!3W%81^}y%^)xQm-fFINjN`sxl{eYO~t_BHI z{=|?_^567Tke9cm+dZ^hlbMJ6!*8SLX~7ku53M4T$nkNRb1+c*e&^P-PCmT15N^(QTFO+wY~A zsSwL!N+s|Z6w zPAjBOWiu5(hzc(2>75nb*SP;!e=!m_^l{fyVA=_Sq;ZLQ9}nM{&i=eiWVL`g2n<35 zb-_-e%NFRPln0A>K^fYj3rWDwLLN9ok`{kLzY}GTJ+59DSVP(h1#s>jo zlHeXfl9dM;d~&=qZT(Gdk}_JTG0i}NbK$~+_jE*?D@T+To5$(fFDRnuN|n{B%3c-%MRB(KqsCN z4hVWMQ78MBtbP0TAmHU~y0@w8Po>xi!C4j9>rp_A{Tn>+}LM=CSU7 z9aym)7pIgrhPm|{L;=P@|4VmhlZ69Jb|uJX{dCvu-@=)tR;atwNT?naG&kypGG5F; z!ohMuTz+U?lnBu+9X^kld#8 z)QUKW3Z7#IiwW?|b&>5ONogGO72R0%EgkKyoJALeR_pNpP|+7JkzatN_xmr+?yelS z8gTuBynBEs!`htkc5aFQ(Ptiv%(~!gab}1vF$N%q(^z<0CA++l{LHEIy zf;wZ)P8U|Q&%CLOqL>BJ>uCIv>|}JH{RtmWPsL(z-QGW4UmdN>TVc6LJ&g)gO9`z;Jaqu6mvk%)U9hEi7taJC=5O3;Ds#)W$?tA?y?XyjI&;6R5 z)d4!Ej?~5Yyq#^=DRI=yRCtszBq?`oJps0Dv=+C%P8UjtYG~qqH0m6;D*sbRrMY`n zA^96Dyg1{3?NrB?2kaB6Hm)rPd#ncKuXl2vdT?$^PA2vkn%V6wIqJ|aZz}fz04XV9 zI3o)TA^Y+OzqReil96L{s|BW>t}7L@tz{Q}a+8v}3v9-(Bmb>6$jHNOl`<9FbvF*| z_Hl6UyqZ^-uzKmt|N1=)I0thfrgRs;1fm%b`;kDTLlZ^Ps?*o2X$DDm-n6d;?;_3Z;Rcfj$)M!Q}1qZV4CHbxl8r zj<9<i=UljsOp6 zz*7PE93^xe@&-VT4wYY>6nDNi76ILn|9&wRKLW6rMQ|~y%*!Y~z1ms6auu~oLX90! z*(ITz0#!ky0u=<1VyEg+MYwr-~hK1ezZn^@Qb{fcP z#qb z9-#8}0p+tc0`R_3bU$8(ZX4>u4*@G)4t2&xO^yn_x}bU)Llg-|z}+jWxYT@t+9r}z z02@%lKmc#i``yRfiB3yK*WBbr-DPinT|q5nHP~AJEkli+#21_Vaei$OQ6fG;8jc_F z1QFCf39jlW|NKgtdhY;K=X7p8o(U-+*$Z!J3DyG{n$JAw5+f+njS> z20n4%Qt8Gsq6dJY<-_P+_g+>dTR*Kx{{U$R>6NIr?un z03LDRsl9)SYO63)KPAZEIz<7Y+l7!M16mZPK{s3eN*yrCs1VwR_}BFGs6`31m$^=b zdOxB34Agi7>j@R@_aWK^Qg=|3o;zm%3Hl7(BWugWTQz_`uAK~y<}d%9)LdH|Dhs*P z0KRbugCKrnj*5rODM@s0B9tiiLRJ`oc4ts^0BX$!fdX)*h!Ot%+M(wHG4)ZD(?RVy zGjvrdBrXL(tJB-f(a#|GwYKJKC;i;I3?Y{@Ko*#*At51r3sjITH~I521Cg(_ z%ceq5nI|DX3{>ar53|J7wwDph02MSBI-VmZN6Nkrl4fSXHfyT1@jjz4S&;08KYKX< zdN=%SgrvorD(JDJBra=WVnQGnol8yO03k()2@mj?LoGa@M{YDIioc-!s1bUI(BX`0 z>#b$-JQx|WrB~Mtr|u(=R+{f``mV#S1~)+3?EAk$S4|Oo+2gu|;CO9HWmM$Rhf3lQ z0}Po6#3)9%2&$-Vht#myJY0MdIoagz%DW0keSxXc+7dNVAeDU7V~K#?-Y!}b2&j^Y z=IegLa06@p-TF?QO+f`p#Q-6|4{ASi?Z*+b1LfU1Qm^(x+=qHD6992kS%q{c1ASal z5;J0=qv|GY3LeIg7-=v;J9al7@?d4U^qrIbb*^Jl;B)V z%mg^y2=OkPf7u`bgG)-zrLHF$Bj{^*mnSR%Zn)395!eI}^fh!Dgltba6njFSH~{7; z^XH)hcs&KQE1b;GYC}A?DzyC|_Q_34$cCYs?+8ksTZ(4S*PEbaFp%Tu%UcrCZ1f3A zvGD`MyuZQ$HT>(3|LdP#gi?r+`uor`??Hmqv>G)DMuZ!91F8#y#$jJ0Kxez+hE|Xd zprNIO^Fzf#XxH*L{G3Evx`a}T(Vho9eB&6jP_ZZwn`(4`>LZle9Q%u)eGui;?pf3b zicmdggR=5kET4)&6{crd*a?e*-}j~pab@}-JB5AZL{ATpcR7;tKt6{KlVoWVlVSe( zHSI20f?T@b`j%hCgnfv8ge>tTR4SpsQN(CAjw*jqi7;XgbreBP1u$F#EZ9HPWUXTX zYC$ScQPPN=3%a`vCAtvT2$dhAEPo}GHWuhg0N0IhVI9lMI~QA8oc|uW3e0|A)OaN0 z8ts~Z<{>C_N>O+$dQBdZSdkxFh`bqs3kAN5j>vo@+%8BcDhK;PaJtcLT)>3Dl>v3S z0yTOqb(*)Fjes=k5l|lnRM0UMRNoa5H>g>vLU*em5+3wOlL78H%4WjpDioFlUk$Q{ z)CiB#%34>ay7ZXG`JXv4z8f$HSigW|WMUu_;su5K=gh>Tl&ja>&Q3_l0ww1j2x}h@ zI!tqm35)~i5jOG<5r;vUDB{H+S99c_L!{7&oUkK673;aRuFE7{g5)T!1=QrJoNkmG zwzbPtAgqh9AkqO^(7OlPKg9Hb4EuqHQPOP(#5iwA4M?EeC8Ci>w%Gv@w;D2GG>FSj zU~cOnmJak%XEy|T>tu?GOX{XOKOi>@&CF4X__g zK;CW#;1NyD3;2A@F#?!^NlGosMo3Y;1tWwC3w_W{R~t8{qCPo6j7?Bkif#=;H9B0- zEasNMI!fE2K#yK;YzkaXG*7~TvZX^AH$Pb4)Fu`KVX#pu=!BaaNYk0^KVtf*o0of& z@+r@Tx3{jC;F)CPDe!`^jGwz-t&5k!HWck*5HA&uO)j-Fm^LP{HAP9Oql~aT20{NJS?lwMY3wHpGUhg!C-|w25Noghgos=< zAOR*ArQH0d$yQUOm0$#NwX%oS;V-S7@88>g&TI&yhz<28M23)mm;mGJ1 z@i&`W*HNOOe>WdUL?3s)Y0P&ysO9Z?YGI4~>*@IrAlv1jn;QY>`h($JLsUVEkDeX` zm~BcS2)B;=ZF`?5Ys-UxxfcLxsEn?IBincQ=m?18|JX!$o#^uo@MY_5bPyX!tV7$4 z7q_n#Bib14b%>ClOxG7e_$i~wpF1~@*l)b&^fqBZ!`npU{$tzr1I|J*?sY1BWbO>X zfo!+jE(xk^N5XBH+l<~n<=%i>xA2zb7~8d&bt4=)z6piXu>|=jTjgaj23cfcw6*gu zF*VBV%s_*c+>*3_BbYa1=TGDnn!>tbn4?AsC&!(1GZWqzi_Oc9bM=q21HC?I7c z&m-7;-ObHjn_{9_DyXU!wl5pV)#jm-v}D9&vb$#_q_M+6X=z@Uy)8{}m%e6SRe+mm zDyT1s8vchZPu9D+N!D%dRKO9C)sCw{-v@#e0SaA+ms|&oo@^eB+dScnt?zR@M651w z+)NedFW&uHb2E=87s%S$NSoHP7M-Rtp zscyjkTf39@YH`=&aG!eJUoi1Fa~Of-3fmI zq6sqwOkZ_SGaPNTpwwrI!oj25=x6SviW?ZFJHLjT`+0yW$IfCXtQ!*3I?BjyKOCWW zrG>s2yRdz3zzYZjA8I-HFkONeAAO5%z6n<3D8Z_2)040WdsnXrM}ZFwE&OT(_*&O3 z@AhY*U9nA~W`9t1=JP(YK!a3vw9SN9?%GVD!lDPJ;&qY;-@~k#pOlo$FB?n2yr$mZ z97meG-gSdH%CcG>AXy8C(z9}-l5(yonLolt2}Pd&G2vj|KMl>0SO{2F4}(D(G;bKK z%`68jZNlz-!#W3e0E!KbU|Zkl?Gp%>YjsnEtc6BEN{!uqTecd0?F6&6#d06mJ_C~F zT-%_UJS?18X*dnPLTy`<0F_4{+v}V9G%uVHSS25&k_DzN-$)V$@a=`I?|L@ddGlMp zq=FzPvs*YE60}80tgyPA7(1g4C*b_1lI>7Teehn3IT*6tvDRKky1GAe4paQsBtDY> zoO;Qq1`X%RhMl>%DHt#1;Q5ylpJg9kU{a4~XGSA#be4w9(Avkn9^;;EVUe2dC19-L zpFdGfxhFh=b-06P@1Y(h55JQ~o=8^~ad2AMJxl$Zz@vJ3$btJCY0%A11=ID!+}1MD zvFzU0S+$r?TVaNu{?u+!B!RTpZ4JTMWTULrPw|^*{?)Xq~u^l)#B7}LGKW1q~P9Dj|kAlJt zS^2tk43!x%0L!;$rW+Irr&@}fkzrPdcvmQEiuFm0pvf_Jtau|Kqv)#>G2^`_QKWjY z%t?AO+r_YvW2t&55F3~9HoMiVK_HS1CtE>U#IcKh13%x=DrldswvHAW91iC8)*W*e zNzke!)l8%Bc01j>;u%Jf!##ZJpnO%u;U1m~5LOYhtqq*5)izDeYBENZKG>j<4Q?8i&JJFXh| zyy|ycU{3L%W#B(FCrYb-(=s^pa27j^;tDGBUZwg{BAbEAeL*A1L6P?ER0vHLh2cv@LIULy1LeQDg_ z(28OEII+}(_)15D2TrRL11I*?dIsb0)+T!p-AIFE-eT$eys&?USGVf&Y-pTu`Gp*l zPa`^i>LYW^yT5o!t7G3s409DbT0T$E#V=icDPJA-W%*>TrE5ho-ZzYr++XUf;uu+7 z@%Z78m3hxMKbI;3x|pl}J^hJ3#Iz6h49|iYneZ$sim3JC&)mIg)H#1M^_IFw%l~x} zfjbdAdZ6GS#H?Oo+%_Yk{8rFg%AouLIk_sXROL7K>ic%1K(;# z^)9EqUX)Qh68PeCA-vy0Yc%_&$v^A{6soHl6;9-niWOc%)ZPjXek3kcGLwyyT-OIz zuz$mHR5hNL`17`a5uiB(IVyr4-nO~QfV~xX?crgrZglkX75>#j`Qs`s-|2TAzM6gx zuQT=|iEpAh=zBx8$9?tq+uwaozWwPt6mzB3>zFYmnVF$ve~agQvdlYLa!y%ux#9|^ z%6P0rRu}!=iDR()Cm8HI1&Zy)I&RsSbAXAupAdTs`#%)X<>vRNcYMTM=EFY`9t=%h ztbQNp#JAvVR7zYugkh0U6wa~FytU$`BRtrzvog*e_zv$lb7lV4KPlJ9q*Drp>A#ER zbIrXDKbB*vRgvmMe@j}U%iyHV!Kc-0QX&?&^ckASbj|tiSRpQ0oaG7W&u~?~o%vDm ze`rPUPV$ugrsodlF6Ijl_Dt}J>$2-=w_qjz>6{4WGA~!b)f=&~8#Z-jRn-=$GvwZs znbBF1`6tblr8cfOBr_s|F~l#5d)hGBZkisq!eD4AC3-xEb#G-5hls3nxmeuI>PfCU zp)1%5PcBtGjrQMWOq$4<^J-p`Rmtg1XEuNmwnm3JxFTR=o?gw4s?)lhztsMMN!!Wd z!4F)v4UfzQA$mFM;&AJ%+i2ozU0$(bm-8YXT4mH?@ngb^B)&#mwta6vsOgVOS9p`B z1x0g}-&0H+jti&rgtB{&4Xy1N6rv1jCi|HUFH(zi%igd?JtC{PbKUb^u)$0RgkKcouRQ|Mh15^E( zwD99OcXcajMd)wIm@CT1v9^|5v%KzTyp?}J{h0B(Nr_O|G%5VSUJ0v@c0N8=@Vl!G zD<@XH*S9-O-z-HXHI{BKX%8RzLtmdBWl(5pa@Mu4vi&^1l~q@qXVD_;$fa-I@pP^4 z#lG#-SJUKTd%seY`^+R0 z$GOhajx^2}D?YHa)N5*`drMoD9X9KB5i%S$YclJK_B=7QGru4%wYWt1BF9Fu z;q<1zVTsOG`^fyr%WGPBSReVsj-s#^>|-4~Rf>Y5h97_C+k{7q=uC~9k3qLQB^A>mt910Iwzn{J;O}^w0@#!?QJREG*g8XK{AEiefhr*&|)a zbu-ij`U6&s?U74PWBWmwJ$d9s@35(HcgM5mrv++JGYqb z#IU#RtXJIaj#3=cH~TJ7n{|R1J9n?xb>?oBKnqxhOOC)2OPxLv?9F z_1w1y*=OV(=llbg+_(+#aJZOSDY31~ZH5b-oU*OGZE>T-Ru!a_gEI51S@WV@gm~7d zzlwFV`P@1XLcJ9o76p>inLlI1#2GjmMOn~zI~}sXKb+D?eXB6jxxmn5HR$gQ3qyX4 zwt|vV^k|91bbiL%j2s(%aJ|%pSAJG`83}=-{2^BBb}6D#rFh+CUsKq^dHK+BIYy%^ zYN~jK7FXYDdp<=&s?b$cL#6SJV2B?{KNOp=3dd_OFgy+kK`cu&``ezIuC`F%VLJ1` zeD)J zVMfitoTuL(hL>OrXG)Fv-PO&^1@c_r0r&E@qAAMZjac;V7SU?gT8s^AUY9O&ATwj{NJ zz#czH#xHFYdcm$f)KagzK16bbJ91-1&usO%vQ-1!qNY)lRKGA1BtCXkjyTkomMX_| zbipzNJb8b`{zV{suuQGS-Jh)Ar7hdg%7b;8$l)6KJ5fP;Np`UKELx!m%%)7glw%s! z3vGNA0*wl6#Px9{`lCZb2b?BOPIlv_ADO+7Vvc(D@Z-Zd&I1JjsSp*pr0bS3rdqCQ zP(44pFj(Q$X4v#K7+YE7U+I3~)7`67=<#2BlYgN9!zV9X<5D&3WCE|gLoDkVbYdJy zsISgDr;2tP6vsxTJt#h=VFCMZd)P{)jKKrXF8bHagL$G_M5Q}0-U~hA_O|Go0otL% z&c5KtvJuzB9^}am zHP!|x+0~zTK6`smUbtY~Wp8D~>m{Mi?91SPUDi`~xPB77`1(}>kc$3FL=>fg8k|t; zhTg-eqsv$-BErq#nU zbqfAteVys~Fv5#9J;gK5R{D8@fDGmq_Bk_U{MqkE^@Ej~nE0h9T5zulnYhm5Xx$Wk z7dLzmx|i2kuC^XVg>Na>OS2AbY;w(GI@NPh^FI?bBIACpE3>Ryp6SoWB@ zUQhQM<`g4$qP?e4(yFE(><2vtq!)@?Rm&{x%+EHE;u8fj{a<&Gspg+PkR~9ds&xNo z+jCh=F}*y`u_h*Ua{VB$A&SVjtL8p5DO%1^txc0t1|AhtI!kR0QVYsfQiU|7x)Tu) z1%+;aPoi5MvXs1g_IEb^^ktt_XDXM5HP-lyRkXuI9*AdS9g$iR~S}w<7yg>VDfRRcdMS;Z#;dE3dej zDpqAHpP(;EtO~X9PHN<0P{@+R7XO{Ro@x(X_gj+YK^y&#_9@0{G|BD@fUriUIbyqZd`wMXm6xyWD2f=n1b!HIwH_tkd<>OziehG+iMpohQIG`w zq+%JP3!$&sejje#qaUC2+AP7lSSIQ1RVnK>Ir(Dg)V!-lpU-%Hd}FDqgJbY2vGwoe z5R0qZW#B~}bgamgmHOHW`Dljb`kDW_g?K0rNwz!B?%Bt#rXQCOW%fPVC;HCKcLg)m zTyu%gW6d+0zSdAZRYBvpt#9o|OI7t2zEtmGX?XEXuM%0cW(KuHup~H!1GSD)H$tze z{z#O4Zaa1?vCC1$r0`l(qSFNHcr%Z@mwX%}y<=f)^d2{M6WvAM<7s*<*#uyVH)>1q ztFlVAn};h=j=by%8)=Tae3HD+Pe`DfT9o+!&8n7UI*$}?t66l17$T;ciCmf>*JVf2Vmmy>KLjZ2(UgLqShQx+@hD7pT@1My~Q z8kb}>OI0x z>~X9qeynWISVsX#_ihOOm6_Cl&atB}6>WnBG}_s4@4OZ=L-nfetr61h9J!jB_hj3t#^c3h<1 zcf+YQ$Zcy8?p`GRIH1sagii6q0kVrvgKkok=2?C9;y+99eA=&Vd)Us$U$qU~ z&#UX;>O+Xe65ZYiv$2o5ju8<2mSmLNT>CT%?I+_mxKRA2egW}Za-_XaR6Fs+ct`k} z#pm_>IhkKAYk0ZmS`kOuCssKqJ*#ob>3gkjqfbEM8y{@5Q%_J`k*&^?VMk%I<8~&4 z@wd-A5?W3|JvAVu^)f8 z?#uIHwl5->mW2a6cEtOF=>Zz!oOmZ9&@4c6>gjX5)zRcP2Gd>Yk-b<-Jo!^mdp*(B z;9%^CS?*m8qv{;%$uAq4?fj-@1M3bd;Rv01S)8I#7x96*jzY5Al9frdtoU!IxaeQ- zzs)=!QqiR8U)kc8Z6U-ezNsnXKX_8o~@??Ofox7QmTj*>)1xbhMrA%rW^^#hZ6j;&^q5?OE|9jXg~IokfS+b%PkWoM*l5Xq-w0W#igBOASBEHuzmQq2fP3 zhj%siwCogxDhz*@TX@Ff_?tcVCR)cNBt!k;?iHH*3;U0~{Yx8olY2~$JW86mzwah` zQzs?xXx3gyXUxq9=|Om-(g-QpakNz9=XBl)R#kiy_Gc=g#s7@W7L=iy&R|^{n7N0o z4yKPvD7Zg!}^>2=3#nT)F?TE$p;(jz6?A4+U*5i%U87v&3UMBvH~;)ir231tPA!YveB?$@&DSXH(L= zQgr5k8YZgxfI>-u4PSQtOwcI^Z53y|X&%Y6l~}!3UaI8iab>$jvs|YWL6QA+bjEL% zYdmD4@?{$7vnMu5d?np_se_4Bn1r|n(s)Cywv|y)ap5auqpnyr`k1ke5YlKK{1edt|tD;h;n4t?#2kd|H&RI6wgFODGusLW|brM0K%JUv|f)IR+> z*YI0{O7(Mm5k2#}`K6(VRP#G^mLXZPt<51JmpfwIaGUI zq&*gQ1LZ#ynbr2}^gxS7?HhA|0N+6^CJ%w-J#;6Poqa@1OQSghqQK(W9>J{oh19YV z-bE@S{gUbrF?h_`Fnpzs?O>syT4gXvCMPDb?|o%qMQ@B4Rh{@^eApeSu2&^~w`R&{ zz1x`Gc&jqpl0MkMnb^ymW_Z0d)iYIy+>F5xGZ$f-(ECWETt}VzJX`HvS-W{^hE}uN z&u6G<^ojjEyG4B9=3Qdx3$NUDjB!xP7M1xd?jfYX&}rXRp#po2(uNnGs>4EeV=EJ9 zXPF<*sjBNA`^4?dFmjt$r}z_%l6cRbK#|I9{5#HQx8L3Eliw#Q|L8!6ZR%j7oW5@c zo+k_z)c+BhEk4QbtPz|QZYATgI&cqvg1GYEo6_F(>L)S($Z1u%*M^XX%A51EDspW=FAe>UaH&dXV(MKV&WF>8_q zwz!-{lQ1OMdw0KW{vZZv7EL8Pqta?GFNMw=IQMUEsXk8{89&?M@aylbC?4R@J)jji zgQLA?D{8lt$YotaJ5pLz?wqh(qvJ50aL;*08{3?MKc!zQY&n!78*@R`b-=hy>k=_m5XOO zxg1HUj;|0O@*meP$ulXQDMusqONR8Wi_sBsgPA zs~7#-^Q8yJn@VM=9P9JB9OlAvCbLddT+`&C(CgvRB$ur(5x%8wXOJ#7l5G8N!Y5(# zws#-m=Y3`kVh(3@+u!(IIx?XeVEVT-KVueS;VE7N(bFDsP4XNL`tF8n!5Fn1ZpW$- z=*9>lZ_Y+}a@g__0}TK9vfJFah_YYFJKDcIC2cKfObS{#jyXnG$4f6?++P}TSe}k- z-(#|i?vjV;$aXRUjpNmm#$B%>$%IdyxDyz8vyQDG{nnMz&~L3#J+!-tzvm{43|rdI zjjydWN8JP75;mq{MC^Jw@PxQX%jH{ zYCC^k$=_UTB{E!6F2mQfC%VX$PNpjCM*LF0*Z{{hijX}56=OY1QeUpFI&g`T&fQx+ zHT;Oe9{wsSEgXW;~%ysvPDTY18nl_JKlDn37cH; zXMb{KJLVg9`PXn$0><%uqH3z{=s4SeF(OL3N1MFsNn%>My39m=m@6M0^t?lf6D2DIqn!d3ub*V9hK==(CdXK--dedW>tOyP zpJ&eMDs4uIfb?MW*FdTthm~XtDnA#Dag~`FzYH2b^ZV$hYKjxYSiOGe?OKI5U7FIQ=h<7`qzbU!6U4zj#<^cw}gC z*rt!AP>g+`^sg0mZsyxwUW1H(%-Engo6<8u3wbHM7xxt%w6%YC=QEL-_ngcS5s<1b z%Lt_DINs0D8J}pO?K9GTtdHHv@Lpa&n}{-{_c47iH5a~dDd*BT}s@x}7S z_D58Vd8KBHR4-p0e*M6L&c~CLkH7MJTNx**r3+m0X<+5IKx{QO`PU=iu!@kg{_oW< zvy&3Jc)?89$*YDp6~4I-o{hMB>5}5hyN>9Ol%D?V+j`ut9EYU=!MA2MdE@sPi|6pg zP1!#``<=uUu*=wu`Avq1H64z+FD;SDQ{X00T(YWZV@InI=eXNw%vo)<;=>PCp zOfM-Bp()1uwOda##P($oN4wn}5Xf zPAA`JVNRLy*{3-R+W2FivfNY@(|?ao{(V;S;a-);p^Dl=Zs+Fo7CyY48tHk)INX0I z=X8Sg->@)ICyZirX~N)XKjqzQpj}Qs&9ld=_SnIcasn5&dY#2+1(5Wd|fV zm;E=-Zh*<4pGD59AFWLPD}xFjm%?hFc5Mtw=sWg{VNeiadAxFvkVEZT`}cCu7cCE*5C8pXyTt9syj1Arp8wJ z**dxEbDuLuvpf^&cn|fGy}93`5}If6PPEU)A9d0Ty%VFj>QI$E*Tq$Ars9+J+xI8N z-iFTwE`Ffg$KFLjNKcOhcCCA;MD+A=`LlBArG10e_gYo@p6eG@EsxuAoxu>#h^6s; zQq|8-s>`g{ebLLT<#7iunfT)Bx8v2XC00z#i~7x)rl#jBDrz%n{Aq)!o+VqXR^R)V zp~uUK!^ur=A;DNS=;gn!TcZ=cOJIrTSgcm$)dvyO1%E=Qi&f<6FMPLA9jBfPt%B$412rBy_bULhI|Ht6eC8lS@>FITws~6O@buuh&u7lI$S|GJ&eB}l$rgr^jTQr!K0kOw-20)^l`P84 zm5OJs`*5zI?3YEQoEVC|0tIJ19|d#zad~Gx3m?=wGX2!oi&)~0>z7f-Ki1V|bfP;n zDfX3og4G+A%d7Cb!C*$+r;3^Wtn=qhx&NomRXYBawwQ`{y9%l1`9cf3$Yjf$rI%;& zdi)O>OhupMNfpH3kCjMl`n3Dw`1pUBiGK`DBB~yxX;;jZtcu7;IA1zzS}stpBD^%V z9AzSe6X5i`lzfc}dqU#0vc%Fzn1FdppWmo*Dlc>-_J}9$8Vj$}XDTELNbYu^qLr(v z@Z6zzH%>SVp zFAee?mFwu+eUW5M1$`c0*Q4+0Z1VE?a6cAT1dB(Wa=A2Q8M>PLFO;8|yOdS8IQ6@o zZ@Bqn(aGmmKDDar#a>egV--N5)=oIM}?@XoxEqF=CbofDnXUR za%_4JCwCW2v>!jTGSb?wnO^6>7;gt20Jc`YEDCIbFa$sGtp9R z224tZ`rjV53(Tg5-I-Iz&Rx;2xI(hh>7$iVbdKvm_CzYLzLrJvV~SBzD*e%YF7_SG z$x_9?B6LSI4f45|4Yktzn`uTvE&@F@I^e~B^ zEv`aVe43Q|l<0ovnkt1dlc{Nh3@n;~d6uJyFl@5vA(zc^l>RZ_T38*V-y<>(Hyr>R z(c@^kJe4$9eTsACvJa&M*VBYXw`Q6M?QfRlxk#Ev0tq96zqvmYdh4bl`lrrHI(^~7 zSwY*YFGH!shHbQ(8PQ?em8?&le0HTP%yq!X+A?MF`Kwhemz}R}JuXYd*`_=Q_-LBa zu=iunNq$%PvkUQ@)wmD=*-kgS$(QkavDJN)!Q*e8=i9LJ_$BLQ)0&vd%7%q_*3jvt z`v-S1WrPtF;x;pL^nICFN$&^BRUQwCzg;?aY^7huEKevcwi554Zu@YOd{oBZaQ_~a zXF{Q(JOh3)%@Q|MBy}xo^*A*z zGVJc4tf3y&Z?OE|c}c7M8W)!VU@v z+vbUl7cCG6b}!8e4mSE)U%C6M?Md)$iXV-FoTOjwsl9akIo;Q?Fj)QT7o#%2WV6SO z@KF+GGVxU8QR(K=+h-Ip1<8(k?I#}hzM?;hRc$Wu@XdoZBg+Xce_yqpH~js6PfFm6 zrXz;KqY~7eP^m)dI+&3YJJwC(&&gz>=lF-%`j5Img>S9_?TNCXC@ymy1!a?0ZV6rb z`F=1$qWAk7H;vzWU6EqwX6Sm#!oh{Fz8nATE%*27(U1G*7VMJxbp#AeY#W?(zZ5S& z;TnA&zH-7omiMMWEw621KUIQYM@NTYmse!pzax}}6?2_fT+jgNmMtgkD>=DO~RO)EL<(CuXKySTrO%WcrVH{#}W)nh#V#7bcims|;nNbcbuX2x*=b>5XIa z9+PUT-J2PSXDPPA+n%*3gxVJd0cqx=i!bNZZpCs6*UR6u6?@ZE*pWtPj)=sX!b~t%;K?OVO z+vC|Juxj*D>$}#dT%!cre{CPjjWL(H-XXVgF^(a<)zK$MZXXso`U%^0M#Zb?#5>+w zF?bQK6pisVAL+}c6wyW9&+UJ`Ou*;o=T}Y{tkQK%tSFxqzGL$;^mm#P*^zx`Uf*Dl zX%oJvbhdcpNt4|qPj9BxivC&We}lQ`J4~tfE6e|vX%@Ut=sMVAI6PJE_aY&;A#2Rj zLq}k4>h?89i!Ws}LOFTX?L>@c{nXUwc7J^Fx?%t0D1H5$DIBfLlD>iOlO;S=WjW1aZJp+AI@Eh%lc4;-Pxm&3C==MJ!p!zQ}ylk(}UPv{-vOPTySJ#q z@2dJ-xnuN2c=SyepR?_k#SCmg_n+{-nA>E^nF*`;Nng2s7z2SSn%o#a*=9KHWgyjt zecVP_w0iTj(oBhbaz*r&2Crf5Qf}o)YkTEkue4W2 zI<2}+oq;FC4Q>5Q90wmbg;qvWO3*0S<5@Wr{i0in`1xMal&F+)%W&Q|{USb($w|m5 zicwpA5p&~E&o~28gYL=cn(F+5&%Lj~ddds`1jUfte3ecd*pndit%(2zX zG%VoN*|D70>)PSF7}{U>3_Uv5=63FdeMY9PLFYS0ip~sI)RcuF`w<&d2 zCS%S1 z@2uyl$(bAPpqY0YWU#JcFz~MnzU1QSVyx?@8PDO{{GHS9iMA+by=jWLiFGVfe^o@O zfW|k~%x=+tTEe~l%gGDn|1nn@7M2j6pTtt8K7IX=hi|&LZ`R7^4^+X+9S>Yeu&ssS zM-^yR8{X0jMSMTn5uESq9d9vt(6IZ+HF9HrA*1UT=Tr@Mp5D`-B#-ZvsUjrY7xO^4gM@v|Z zbfMWRT(i%tJoik$_($*a5e@DK0Xs)5!9qXsSYIoG5A>cc+oEn!2j8bs@j1$&IxuG% zTcy4#nDq6U#;6WR+4&(wF&XdT*-M9Uha$hGH z?1e)7>Eq&mCf~A>&JX9CMb>-ba<3gT79H-lwyyiLL;l`W5Xe1i^teZ_YXgx>u&o15 zn2K_}fKEFd1 zMmYlXLIvVNb5$7<`1!EKi1I8YN%+>u$28aS&(^8lz;U-9j)`N9 z?(cf23wx;gW0s}$zDJ1(ao7DkdE6{XT;i5kUqkP^l-#|HZIWj~D9xyM(6CrEOUhX! znq(Yl38z{8aCp0SDjD^C=Dk-t9_oSjj5*BS`|(G=ZhAkg0nw$3S4#6YGbj85x$bwo zeEF>(}7W1UO1l-$tX;J#Oqewd-^6nphX#L@nZCmGh{@4m-Flh*__`;H>)c>u=;tVxGuT z4fxr4TSlig<4h})h$C2qMD*^tDbFu^KMo20;2fOzomt^uKD{6~SazuDTTmBy&7^_N zd1{*Z250G)3(F1_Pq0<<3Vyn01i~ccCK*Y-g5x-I;y)%$t>ZGtY1n?Kn<`X`X1vwc z3f?J+Hxr+d8*m);92R0pF^K9X4o@j53Fx8L10C}BH>f!Eo(Lf+1%(fS;uSx&w!Z(aqEnIt*Yk3#k;Mp+&+ z2IZ(=n+RckXXk;b>`WII@m8_&MQg={VaoqU)q6lSnLY30D68x$iXEh?Ac%+%=~8zs zh*T-ki$tj+V(2Yc7Nsw}7o|u?dJiZi0--~Kv?xIWgd!!BK$8D`1Ac#>|DLnQJqsr9 zd+(i@XP%im520+tZydfE*WUJ0gf2K^PG0?6F*|>#&xrKrFWKyzAJB$H3oI79LFbT3 zD**b|fP>!};f_yfy9I=@%xa#^ylp7KHKK&ij+9DYpQN<;9t4)&e4^50r=3GDTiSnq z!K{3!eBhnl4VV&R};E8|y)GifSNc=BJa^>C5v+YM3lnyS!_=wLqe3_T}CSff&0$mFjk)m`O zfqMWrndUzC$N2F4%9ZGxp9xtgsqmg>7r2v`A|y5W<+P&5Cnx0LEy#vT#W%f4E_|r< z=Y&Fzl~)ny-9zzdcH^(+?l3?4Cdzb1yY>O|Mb>lK`umJu8RViF+gQDedpq@|lKy=B zmax5j+aswp|6578geB$uy+9*^E;!3aa!b8Di};zS0%|^!EzmAx zBru4-p}u}Bf9bAXNmv)?oTy=oY_c61hV8_rzbhBg>K1~!SAXRe->ScXV-w7TyY?8* ze|zVmmyg|8^cA=pYPMBFmeeG~<9CrR{hgbt0QXF17G(?flNYpGqsmL<<_gKV<9+(1 zKUpu%xd;Z*?%rDabU^FibhbY1fAf~y$f9Jv;MbO(YGwnCRcJzZj%KO%oMn4$uRr3p zm~$uC1}2K(&#y}ua?8)s?ZbJ9(^{#?z55(A~)k>UmQ2FMrt$)}%xaRX@r^w*kT-lkk<`P$ZAzTO)%K zU0YQz>RQ@}pJ4=1fm5?|=~YuucGJ#VYEmVyl%ymGS_)d3e}TXev(L6(X1Fi8`W=ch zdY22&!5^Ke$m>CDwM9qad=--3hGB%Et{olCTRSKWxyJ@j>6P`5UW6{0&}+y`ja+@AC7oIeL}wgVk3!Cr8~3UBbt0+b zIaElKb9Zn3k<3@^??6i08m^{YzFxCgl9~7Rt(e7XSYkWY$mjJuW}n9^z{|E8jFgaAA|Y9{N|w_Faq|U{IhwmVikIHB$+k zdW%LAsbJ-Dv8hsqXFB zCxLE|n<-aLTdtu$zg9Y`hJNeV!C1*4-C8pe$MevQ3RseQwh5Yvq9*s7uQ#*@R&n^A zt*z~~9b$aj6SO^V1EN!XYN6WVqZc$~^vG4sL5?)u?m{ULlgj_nBewEq=lOYekDer= zTdT*pGbPy`4xde85l-7^)z{$A(AcvHZAP>^tYTqdTbrPNdB>@;r$FlALP3;n7hKvkx#R~ zAb{*}iI{zt!~KHEYv1yf7>vcN?ZiZ^(MCy0XvlGT;XDk~&6n!8_@){h@u_gg_`ct+ z;2kxr#7h3d{3FQd7e+aVzkh%LYx3V3NN*dE@x2~ZCTxQqblOQjI4}9~#NpT%qXs0# zP@fbgdFW5-l_LGAstEeuX5Qa-3?k<4aang8(9hW_FG|S?LM~QosK14!0BdqLyj$$+ zg&Pi26(g?wI|bOy_F1Ey*4}n69PP>N7rVl?w_SyR3sFc7`s>Y*MN&>>xQ?MWg?MAn zk?p^Z5GYY}r+cfLRS}a7BJ9<<2=9gb1bckHs3SC>{+X1}s+M7)6UX`N5B>-;*{#5# znQ^{XtKtQf>#&^NF1S`()fO)2useoIC4qCnLW)}aPYfleWB%|ik=MBJ#(QwknQIUE zkUNfp_^6^>yS=s?B;w&)>d+|V6fUYw|L0oWti47E8=2$gjH}K>_m%HLW}o3ZrS%z@%Va) z0m;c8{?u#U7LJQIfUIdqi`P?OM$4P7?u1%GILVcwiPhbTdsgJYOH1I`ZCVx_lqtCM{? zg<)o5v6(CqxCn;}St@mqcQ||c&qq4Kj(Z1Po$;nxeOOU0>-B0&&djJxgNR!Cm!^=* z+LT6aJt*iG z&7!1{Hs|Tvf>|=eP=HWBrT@7+~0k^Ps2O(6bZELH#7ccr$8EZuOzQlz(9bdSE&L^{Q zvA4cbrkKQVa=zGmwXCcHQADFIR^^^cRux@113(4k%o zJ^1S%|8Th0$;*SJ3;Z}hsq8iQ<2kp-?z*E(Tez4(5#DP)gpgv3f}w9ie?N^gU*p1% z%Du=ewOZ`S?<<_ZWIyK=lsEp~U&h`jdOjYkOW5hEHYh4_y2&4hnz2M^fK_6eE_$!~ z^R5>-^7kBe(l1>FZc~U6Lhqikk0AfvSTw`#&rVEClx}u4@iCTx3cDeQTJ3u9%=I*M z?pbY)*{jZYe*iNecs|Os1NfJ>Fe=O{y>HqP4E%ld8MTweT8$W-Pj)h->uJb3xrYaV zzZ_si%(=?$P;j}S;E6A*v6^Vo2a7Pe-1ecs zG68-G{2=2Y=O9?bXhHBAs%_w_uK<64ICZIKp#uYc%ycH6@$G0kM|sYbRGcbeAki3^ z;qvI3->~vb3e@+V8{00q!})h@>1CooA2ZmCjI6)Gk&=$RzJI2|7-tGyLpZt|>J>{) zti)Ex-ECS}1Fuhj0u$@m?nu8(UkLUJ46EEbU=LYY{qJK0X91h}y}A%Z%BgWBvFBrP zDQm^D46fN0Z0K2ti$+NiFzeQCJ ze-UMaUFPw3x#qWF`a%n`_90!9yxZ#F#YAK(r&A7lTLGFzUm8#;_e8$oc01yDiRt@K zL)m*8r0CB44*c1ZTQ*L}(BQd}T@8WDbO~rN0;4uIG7TJ$bqk(|`j~}w00-mXeaZoORaw54U3eYW5?jTc@leu zaz@S2mVEe2PEH(mv`X*5ZRU$@@PijwmFzui|JQGU47khQg#5k-TVHi z7M4weqF64ecAX30Rf2)9Wc_DiKeTW`PPq9$^7mjpjz#eF z;`$TS3lyVW<=#9KXFG=;F;?!kf4w>Lp!R(!a5l7@|2-suSpak1q!FDSXboT}tLICC z^ANqAa!{PpIU`@K%q@YUukuFN1D}x=l|`&$dGr8ERTVoYWxQhc)^dMiqL#L3P3ZHf zo3>aWMcw*MY@3vF8uCha%4q7&L9mWU@G(!p=ez_Y$${9eV^Gdro>l*Up#^e1&Ho{` z`d%iuZQ>1zmTRUPYsWx+qHkgS=%B(*bL|l3)}mo5Bn~Z zargsM&f>SMG=duwD`Q28?Z(>PUOf!o@c|T(!hM&`(Tc?ql5D!-m&fQJ^?eWA@P&bU zjy#Wa)av9N{@0tb3fEUTg(IZ>E|y~x<3Lvf9d=5_6uemTf0YZi_)vuWlERXUIjHj_ z0F1pFCxBMW(%&dQb(m}48AQ}uLr3N3duQ5{(0ijd96hwX9%xv8k|f>+tE2M3)!n_l z)s&Xwffh`%z7hwOE~o6DGLAk7>0&9Uoj%13qy|SPIJOZoNr=xxJ5PX*J#e;IW?fq3=upCjN$_Be7AbJv2E1 z@!@6cZ4(A2)sZoY;12$_Nmnn4SU#c&n--PRGS(Am8yuOnz04}Te**dT_UaQT_@KRt znlOb<>Vje?u;Dt`3rn&0R$ao?VNnzr9_q;t;?I%-_VdVUJ{?6@V# zpw`D49Il=F#q&LsB%TZUn+vs3m<*3mkz8)763%soDuu*Dr~`!uh_)dcfjX`1_T<;bQklmp=-xb{{NIL z0vyvX!1-r~A*Wgd`Fgrq@TH=mn}4UO`4md$EXw@^zs{7=#IS4T45%uBoy~rP6|j1|M+h0q>voF7dc*R5{yVH`-@s8@{_KqLidw{?~9>vD(${ zLe~B(>bdT-;(pT(qQF62neZ9(u(sH&TNwQnVm(=FNaSL;-I$%GwSJC7?(RZG*=4B* z&v+*YBoUC&NJ7>Hg6kul{al8&_j=i@qICUA=_z4kq9qjbGoEiOOq{lL+IE1#RA_nn z8kUfjY07-@5b(Zt^wfL>JG5V<2p6rV$7M+j+|R}``ZH&l8OJ%J9JFPi6L~gmun$@5 zWoJ_S?(DgYOSSZR`Wbk1OJb#4x#dYw{&sG$9SXl?IyX)}ObvJLMoCt|PGNzf6F^jz7loCU`lq@+rPIB&v0R zRf^LEwr`yWx>Lx)inNQVi$%=o&;bhY--~$q?(|s@ttaAXYW{yi|NqO}{6gp;*8h4r z7O>0>wsfy$ZBWrX6dlJ2++g?v2pbiJXYK>DAjTo@8>yRnY2^LgHnH+;XE^~Qk{*w^ z1U@X|-qQ9rL%4UHXnED@5V1XXMUR(yQ?92Llq9CF$@P3J(;q0jQ|c;|Wb$5#2l>bE zK3PWjUEJ;|1S10-5uGcWle7@o zauN+Q>IW|gQX8CEUwrd3$N%0Tst{~3>{Z^DhVBV_@7*Bpbxp&!dkz=?Rq%H?&OG#| zpN8Euk#IC+ZxZZ0_Us{WbOG-R7juxsCdpgV%8NUHX(iYl{rJ`)1+N>8^++hsKP zRwf+{24JW&oLdS9s|_o6zi_uT>LdNp{t@%?q>joif7MK^l7lU(48b<(5mMfhlp zKH94>C^9Zg3{_6x=<36xxdn7tv;bHDm45v?5tq(+F&pgh)_1amxf11-Xvfg^R`-)k zgcI)PL0IR2?%dFx2h(QQl^f~%%f`~jy%7?QfLwd$MO3`F7lb3nea5um77)B$jImYi z`t?A__Ze{0oTvZnFy&O#ZiB*c$v!o1URrQSM-`fB4@u2E{Z@AQ=FOYd<25eOEa@7c z=$Va`#laq4hMIDe%UbOu!TYD|vm zut1WcTggq^$&di41c_fYJk+pp-e3iYH{)O~Q>`h3cb<-Xx(JB`v>P^2hDVx^$YH|W z2wonX8A7vR$*m;hqcU}F`_La&I`>6lWuGh}8SY$mO|9v-VOKd#a1(h^u(HFCCzN-; z$P(7ufTx%1c`Q0Od4&_iAV4v&TJL*_V3-cmA?)NNFVLRrA^`o|XTvj}-@OBh0q7W;Bsj|f22DfDbDUfDmz@)?*nFkpqoSz)MtOQ0 z(5k(Wu?_Z(@I=hV7)Z?HJxyj*bnpcf@Yyls>=EFe6nWm*e&6{riiw*$CQ|I8oDcI~$83kS^w8HJ%re9RcI8 z+aRseblP0{)PTAe&{jJv$;I>m9ngAT`ON13E1?qF=an8X zN;{!aH}6cZ-x7i`v>3?QDnL>RCr_GIJq<$cqLHDSjmv(_QOjMHv=_c}>XzQ?BuMRr z4Y}1KODZhC5g`_PAF?zfNI&ytrIaa;&7xdwumfursfeV5eJK!u zgr2Qi3rQsL*x$nV%yUtGDgp#(7o`jF9~^=wzI@N#TX2-;IG4%+ziEEHAh0}&6f&f~ zyzGgV&*w#gq1|nepe!51iMN6NGW+TnuF8dlo#yI&haYM;xs zOdL1CWdJQxdL?EwJmG|`1@+7CWC*bnCrkPYBKl6d*;%-t)?)mo8h9eB`w-4iik{A! z{@cgPVc%^TniUs@t`b*e&Zl1BEQ*b+q zcoyZY&m&*}ccvwM5l?f1^keXN7N&e-9~>E$xq8pWiGX3hUHe(D2)G!lD?ZVV0+c15nJJvuESlbHwyYJMuqN;PRR#h| zI#*(s?x$4!Fios1>X

{#9y^NH;|Ms>eF%N8_eI$cccH8WWS?yg|wzUlu&sUtH$ z#GMASRNn5;R$5yBj91#Io`8hdfp7Ki+DKjY71~&u*qOyabKYaR(zuVP_yg<}+UA3d zGd7jm{gU80TKGG0S}r5OA5xy%IFaUhL$ZdDxQf$ZurnG+f7E-+m~nv44G{aSacPjN z74M1{Q=TCKVUY@4N$O2EhZKLxT!jBbfB#8fR`VYpM>s!g08-+V#c`%Su1)%`2<)m{ zl5s%_IVxokyfo7t#QH|>lVZaGffa$7IYN$GZe6C+mLW;@DM+}x>G1(8yHS8MWi$P}KvWji4$Fyqw0#%?KS!Pn_cn2ZiezGC z=Ro}C&Q9!okPdxwNixZNyp#S4tljH}a%!24iaX1UxuC&-kyLrPFxeTHP*Uu@dY_l*fM) z^J3!gmfqb}{t)w6ne0AwwxmZ6x)P})DO%7JS&I_8i~dkxF}}QP2>Oy<34?Z7J5V4H zPvBAX{j^W9UamJ`I`m@e*AF;~7syHopnZ}S;GrC_4`@VME|~kzx8@cluo7ykzz;3s z$)+C4Ks-vaxM8(3LHyA_!19D3SHn=!R*j(HqgkLMfVTgbdJfV)+PQ4lT9Pc_u{b*X z5eHVb7(dAu@X~OuxE1rZq^r>BN;)qh&MeOP0<@1hRZibdgp9?hJ~q>}l%KPAanVyx z?@fH1x=Y4HO_d{NVT)&?b|6p)nR4P87XgrNpQiynl^o!*q6O}eWrjEgF7S|8&M1i zB}9Bw_$dqt8u)i?tN@+d$xYNIxvQ4xN=~H*?b^BQF$qU`qJpvwx%3AwcyJCg-a6V* z*kv^d$B$;T&I4Y(x*ZQi+CCVEd|3ewqZ0etXP*$axH^23iHQz+nimsFl9Q4szTq=7 zGjJ4baSkbzX~I8MVLvzpfgeH=3z7h72B*_^i7xYbAC>AL0E47&BJpCtopT-^Ux@{E z?h5X;^!22D^NGr_*O=Uzz>MY$fL6nq%gn2&%K$VcWPwxYn~hE!%FjD{xO$^yXP;2A z0Q9`hg@{zyk^cUEBBG+YdegJ#w85eRA9Yq04gjPg${k3cuXrFgJE%8n6|eB>M1~4bBKFiMw*vKR0N=~CiJsGus^Qyp zO`IekU$t9P3Dj0L=eU7(iR=ssMYCnLnzdJaA%2dpM!3An8xsqA+M4@2=uNHGoY~2fK(jq4CkLLM1L<=qr2>)o<+1AJEW`ovnh7L927Y@W&e&0b`0suSkwfk& zd0gy~bquzgWSBDj%_nW+;dPm`>FxJOh~^QvtSt0M$5Gjc1oCAjrge~C1ibQk+JetF zD4?uU0AAemM4dLc9?vJ^@~N@0QJV-VB2@<-1J@;#dNkm3QNwC-A_??lUaQ-c+YW>x z!;$RrOvRbS#*5|d^sgcC>5(923N1S@@PMN7dMqec<>n&k&YZ*2HH+YCC zi12L+02gy)$7-jsb3hEtcb(MU-fq<3>7sl&1aw^>BDTHJ;XHoc)Sgs%@ZUhQPAK)S zfX|MAes`IP1nz-N!`E6;TIx65H(c$n+}`#cWAYQjetc^uevcaI*ib$cZRhas=fYBn zon>+8tr?-=vjX9KFWAI%>Vj=#j(v&6PIUV2bV5l5p-Jf^9c@?$l!sJ~7~feIwB{k)%2hXXN|ys1{{cEo$A! z+0IU3ZZm)u14uEeQ(@RMx0l_=w!+^ z(*YA`^tBT6lW1^6?qe^g7=0v1|7u%-1vwv%g64#a3FnYorKg4a62fHK--3vkk5*p&%D*y-Lv18fRs(;A*0_O*66k=@ami|LL1K*b z!<_j73HYplP*jD+`R)?DVfGLKUjL z=}HlXjJ0DlQ^&J1Oqu@}aQnP$ulThAIJX|4w`Z~N7wJ3jFqk(B(44HJ0)&P_o3=ul zWNjABm%R-e0N@D)&_b#1VV?dJUiheVDcomwV-aX8s*ukmv-aAHWLRoC(d6C;pd|f zT1zxb{GW~QrzHc@U)7Qq)U|D}6b=4^ZdsrNd+>gzVC5GpY>X@VJrp_r-&5;A`Onn( zUe9)zV>(FCA-@Cz+JcZt!hsR4j1!e)q@5>sXJJ6|Is{rinES%Wfv$~sHFm_JnHtuSQ4a%w?j}e3oig=CqOC&1DC4=*M|RxI#>q zK}jKL{ms4mG+Rcw_(_5DwtuMl@h4c$Uhii;_8A<|k+=|T^5~a*jnD7&;HdjarX%^_ zh{^ywS+g%EO13^bB!D6$}R#61PTO4v_>1j09=8TG&8D~aGmHSqu1{!0QSuFJs9iI4D1)@SY`ze9l2J%k+?IAuwlHzWWlC7TpPkDNgS7tH!L z<2)D?JZu{_v`kg60w85;m#EAKVLP9vMjBB9*Xw3KYyqQVhsX;lwu8^jvWM~AW_TVr zoaO~+VFn5nm@rR;A@P7Hk;)VcX^%Te-LO#8ZM;|z<$(QQbQR!&4BguI5HW^zSJ}%n z6rWd~%Q&t;-|medD$}*Z(!0YR>7f|^OL1Ix_Keb1+VeKvC4+W1#nQY@E#F$Z=*Px+K~0_fGV zdM8`6Tc1cOvPWF2*{!Ypnft~@bQL#v_5!fk4NO%|AN_2`9>72`3|t5<2PN_| zck8b12sTq;SI?QPSNsNipfcl!I|h~)q)R{e?{2QiKXm>Idfi3Whk7C@8oJDAM;PdA zoQY!2O&xdyarmJ7rujNZT6F#d)%eoyeIW5?L0g`EB9jb$)cjx75>0-p!d2C#$j2YW z#yhX|v|RTZU;=J~f6v(S3C{P+h^khD$}{vyjaIgxd*I3&mZdGV&~Plqq|3Osdq4k> zkM|QnG|FpQ)6R;0l~VVx3x00|X!kQZc?+O$=bkC=zF&fAWDv(!-n482Cqcb8dTLF~ z-+MoD$vq+UO5kwPVnMLeTHR9WZPFr#2oh?b{~M~-vVIDhsOQ(z^n_gXm7DdTWTVmz zs%08uZg7yHe)a_IG)+F{Hj7^i@#p{B)QFltJppKJN$B0wK)$;@H@7Lx zZ2S0yQ^ z39SoMlC=-ScKSd3qVvg(+PQz>K^DXWCMKzB^2YLGMy74C(XhXN5-$kaH*hj>+7Ws- z*D-QoHlo2Ftyj+-8q$$|d1C%SRF2lI@971MAQ;y^vzmY zv~vR8U?t_@k#>EJ;M~HOCj#HZQAtpO6IciuRZVF>hh4s*NdAq_Ra?8Gkz&{BEy$qO zEN>pfbS{g2F=@!<^)nNz+m^kVfe7e-|r|ZB-RHO9Z?Yr6g zurq65f~R}~pJY5HP7A(wzw+Bsp+J?f@Hs#d=bBxInkXW$@isQa+G-wcX*_Jil*3H# z-pP_Tn!3%D>t;Yk0aj;rZx072v7ATQ#6aGf=xT6#s$&i;j(AK_{bpRk36>jIrHl4i`+MB0G1Psf!} zrbzPo)sx=#E=#mrno~s)lXxa|Rv9G|wyo30;6ZYUz z{FjVA#od^Zw3H}gN91fg^wzvJ)VBtXVfL!f?44W|-c9gXLzwrm%Aie@t&Ex*fPJn$!3IuzqJA+`6S##IIt&wPbm23`ui_{rH|gdq2sVE=)wv zY3~pcZu)5zWu8?b2X;B}v!nCg^3pfI#iV{imLFhZCeG&b85S*`_q&Vffq^62MufeA zY*S-$y8uIW*1zNY$^h>Ou$>J-TqsD0WL&_BKz{pPJMw;AFu*W!rNOZal^VostJ?8q6$A-;Ii~z>i6l5+cd4`Z8EKR?+dh z^Fipu2Rls^;X`zz{Fvz#N7G^vN$|vtyZxXNCplq%3KSh*#H^P^Qc#wD%8cG>*qpza zlhMG*_(Gb>_(q;`aE0#N^1Ihx#?{j;E2|dS(|f?hURI}MTv6!UTovd_X37W${2kbo z)|{~>!%Atw2iuth=!>s>SiPkT!f294_~1+$AdUau7e7;d9{`c})4g~D@_|E5E*%^6 zP)om!Uc`L^OF^yF^Pq(EwNP}_-<*4-S{@MP;*yFl&`5@_c56Rt3cggZlT-pOw;3g! zjQ-fY@2LeqUmR?+&k69J;O^PI7YS9^Tesps!7!{>;EtY){EHpME;Ez_*;55gt36;hoBJfK5gSy{pr8{bQU&&VsqUt{_0kXd8;qJiYrAA{i7XL0Pw(v&ft>Gxvk3{bb~@;7ZOu; zo~5ngVx^0Cpe^?LBsr{x^ftX06i!)TkQKek6oEi(z)eChaMb_Rvnal@R$2+gMYn0r<<;JBxlSv>XkicstaMG=4+Rnk8^GPO0qcn zc4L`0L!WHF^%Yj?Kv6Hi_SP+G6RMnPh%Eag_*KB9L)O(F1zHGzRqZY81&^I&>Z2J- ze}(Xj`LF5yqb=b+CI+I{ad~-@z=%${ZUS_M)+PTs?(kk7YA8b+RqsOg`)sjGNxv|; zUUV7?y_Z#|Y)Wlvt2?SiH9V}y(60gw zkK=8W5=5>9*DR)6o;1al%$0aDiBD@z`Y~Xy&Rerm5J+BB+7ht|xNR$&BBz<*{cPM8 z@Wa7SKdZdk+#o_rFU;v-Fz7dcYIS!8%Lv67GFjC2He3&3wl&Zc{4=vSpF%n0(J14t zROS*eFdPC$;!^y&j}kWKn(zdXq^RN=PrL_wmHm3y7=7i=%J)W2UZ(Gl7sl)4Yz(D` z8Qb#}I`<)6Zcmq$7zwANUFCsOF=!i@3Aj~I0(=R)(0x!)jqKsxN5;@LWP|A zK5XwEQv4QeFCA7Z9;NF)maN8hdO~-s@JuM^mD9Y;=skOIoRvHPUn;J zHmgu3Mfb-=~3n&7lFjY#{F4YI!G{ytB z1|4wWmETDA!h-3jEf*3IwZOb$)iRF3v1*UPmkll0?X!1UUmdqb^w!>(X&`OXkQO%4 z8`CisP0yZkhpSg3g$0otUut&vd*8)^kuq7a3xV299G%iipDqcaCmFtbe-CZC6qP|q z+>fS*(72(WRe2md8qR^DePKt^pPKD$)hn$iLm);;xwOOL&G9Oz!Of*^6Rv93Vo zXcd^7=#)o${%Vjtjmx9k2Qf@Jze&IK)VygvtNhr*FilaHNd*2>G`{i1Qtei{wG-ik zQlr276FswrRlol2{(r&FPI5#%{hD^^FQ!{+T{8=BuODDiaUQt?R_O&O-oQwGpi*vr zT5LaAlw_Mg+NPY^0rw`)DGB4n(eI5zMrR+5Cg%XC*{c@(*bE0e#KuWIAK!^PIF$TD zskBIAuCxmDOS-qYeqFuxP0I<{m+jwQ+=5^8n=gQWu6GV`>#e)Kbly$VHM^0%EiM?| znk`S0(4y|%aoJQ&9hGoJPI~ zfv#lm+z=k4UO_a9h_)53rAHx>2;vaQpO|HI6@+IMtsuWd(%kHzx378a`-^uQX7yYo zojMD}=4-^3wkY^iSJZp!rCMC;K(#r$kEM_5wdB<@*)uKP)X^7F?X^Ws!IzbhUQc(M z&%Kr+wn63~KllQpxZJVt&}_EiENzNQ z94$CR`O!06f|owL4SM$5T=Vl{Z}!PPqqsyQ$~mizmw(sJ$M%}6X|MdM+o8aQzsuND z#%q3&`62?Vi-Mkn&@V~H=aA>f$=kZY5VMrqj$SsPHP@pjsY~#Xx|C2mu!H@4RL!T^zBXzc5UWp!Rzk z@zkcTX39y@Z7GOi02@lkl6~EeMUS8XTv;Ts2_U$Tegl<;S!UllkX4itWyR>*wB{?lsRIXEIF!ucgdw*6I5s) zSJ(O!Wo|=P3dYQCib9o=ArZr;pg;IPHS%Zg4O~QyCd_tD>i45?|hesg6t#O1V8UCvwXpsR{G zfrIIA&u1x_$0g`rBT~K=s7D|^m@?wS0ht$LqKPC@VCs1=>ta8Ws`;q;U=-_afTW`H zvfRfxo^F$&+Y#yV-6qo|AJ1h<3KQd)FI|dCyUSuLn-IFR<`dKTGhu5$8yzWGzw-PM zE~$npvb?;|aIb0OQoO3wkF4(eYChT1Tf?s;MUYFrvhJZeogUd}aXvZ4xEPf$Hl?5=r> z|6IT0($D9ADQOZVxMZ|ySXwHXsa=wu_$){rM>Ci7efnU0L#!n9?ZfL_8Cgty$kFhs z_VJB{9nn!Ok(pcZz?}&~q)Bst339o_+H(M;q}5x0{oZSQJ_+1b8)#2(!=9D(#kh4jD!f3){b>pv^V65)>{yqiZS{n7sy6=}lHtZ`X zT-O8G4vqWS3A?Mvh?YpIN9U_07DAO+)8?CA?PY;%-D`F09)y%?`UQ77ZP8_C2p_@> zD+BfSOvDqj=jgGi?{8<&qUm0yPw{2`0Rsn{gXDC6`iHc@u4yR`LbpG>dwamsiD#Hp z*+4$8Ys&7a!+zmiz-f^ORM~EcO`*k~wkrcFGN>X9}J@ zmik4O4ReE~-fd1oWZw13>BOT`A8=GhN&0qD&067Rd*JFEs4=yU@Wnsd*3j79T$?J6 zP?A`x-}V9Vllk;HCW_V2#zK=2md*6k&%kyDnb6A!5Qm}pw-_aFXPd=edF(&;gBQzk zqtsoYqQpuUTLm!NB%)t$Eu@{}$ldaDcQX5vzHlk6WNEwfluw=!d260v?5I`uVSPg8 zuBqvmT9fwVfC-E2U!pClZQH#*jw!SyAM0+tCeF>!h8uN*liIhP?p<}6j7+xf=(Xfc zhEV3t{M#we-Uu0QzN51;Ct}R-2{v3ybZPJ8?RlB5dgd+K3*35+I!9e2l+b@0M`6b<6f}@*iS;)5=q*djDauWdlch-E<^f zBS}f-y7!`T+gA^z6CUaGNTb1;8toFk$cnZ$pW={jc+VGwlMem*uCv=z!9tr+3K;pHk@}qbEMcygDwdjYjVBwJJ1r4 zV(fUoE^`}6yAoc(yI8$RXqBvF4MC)@P>i;gOyMm@{m*!WB5k-elPRv@ZK)!JUHQp# zEXoxZo>4Jl)@Xhb3t0Wzp`-q$Gn88`akp>jA$rBMNXQSML`v;~u$h%g7L+NZ_?`Yv zQ)A!b$^yd`m3NG0DA^Id*NFm(W6$Yj9MRkzg<&S1Zm36;eEnj#gR+4+4;^l0N|7B$ za#yGxWdcESC0m=a9Q-%B)&`Hc?66M`BbU)?$!;^3T|V<63lr>>HyQ%^gaJI z(y)X|;DC_rjCrIp#4U6DGsD*a28nq;!BJ)QisGb{g4JV#w$J)987Lea?eR|!2o$n! zrcm!YQTC=12Wsb{<-4EN=J=B3jt(wGcWbk)sLfV;WbKV55&9(*#zvkcR zNT5B4x5aWs-4CU%Bu{;R|C&ler8H#ehf5TDYNf1g@U<@CCQnwUZ{+V;USh1MO9wMw zczc9s@DuuOndkFuF=+2^eV!bac~V@R9R?#UEHuf0oKI-a^gdKH#*%d~l z3V5Fhro9hJh7BIa#~aIgcGk*$?1+oWr_1N4$)N@MjsjLGDQl}Kgrn}N{izn!^&+`K z!zck8)T*A#g}XL!+S+_UKhu$>^d>Q?!|J*hh^V`WP+uam3vQb)Oe!n{pr;>dK4vw` zNb5*Ml2^oFb*`g%;wVWV?=dUa!ERXa|Sv=J4CP{-$}T*kZ;`C=#YN09=K?+p6zN-g$G*MJlVB7IC15&8#%@z4dPGU^`?1MP$unkY z462O7Wd8ljUy(k@2bI!}0P+dS{mlu_{WmS?571W65p5z|mf7*?fFk{sP~+7lo}H|G zvYxZ2>B3mOt%^iSN(J`75!;O@^{Fq)v=^>G>SG-BeR1XajY{C4cvPo8q*{h_HkotF zR~2$)#4xd*11Ud8w6(L1lc%DB>t9R(|CEIQvskYe-JCE2d}xZ$G=mQfO!Ee_7)D)z zCjhe;U@vw-2CJ{k)OeHf3{$V9PQddT~6|d#d z1lTK~uTXpkbLwTrlWI|NiJYgQ9w$!7ZlOhbK0y1KNRmgVMLUL{&UyJCRH(3mD8<4e zN#ABXkR8^uUOr7Pb)ZHuE58E;i*)9+*DoVuI$dVY_J(IW)Bj%ZU*vHa15=zQU{!`C zJjmzd=0N6*B{ye`|Lg$h-r`Vv1ZIGNue-LtTO*BK61eIJB@z4u@D$<#XX{@+$E7L4 zs(o(R6!}sB$d!YEW^+zGm~fJ$FyMUkhad26{q^SL$u%*g#kgN5A|GozK`9<7UyC#= z6MrtKozBkQv;l|Bi1EIw7o&VM`sS(?fsP+6_I-6gr&{Y-)tO) z%nWXLV`L?Wl0vzU%g)`=TUho9ya1Nx(^B9r=bxr`LxZp%>5lpX3+aQz%1Q2OLIzK&$a>?<6G?0lMyDR+M=y4PtG0(SjYD3kQB$&->ylTYVZ#YNpQ8yq zgBBK*D9>#c>#p`n=*n?$zMDCS)!^P*!c9^rgoVTt7-R{>U~OX|nMtPNS}>3VU?dGz zSFS!|egp=NVC+p}jUAI;`1jarnMIx1#3x+r1pSfyL-!-4!aAvF7g@A# z$}qa;#3j#q9sChE?9rAxZD-?z*%bl4i(66ZF=*8}=e?zIZnHk^eb>jR6moWC-IMj^ zpSd6Mo)`)A?zp;UeH)aEcP)N5hFX=8`#8_vJr960xWi!WzLSJ8*@-y{b;UJX!J*pb ztcWiU}R=5jnSWkvwH%NMl}cqTU%-1j+B zIHXV}xa>noWJd?rI)g0f*L!3bi{W@2xL}`gPkIum%DCtXXUFzSod#`T9gJ!S{%IpU zFuyX<)-slHTaV*Qu)@nv>t9n7NnFV@cLzTK9$oG)U{whfHTVrygUgylXAJXd-I0it z&8W@x(VR9gi!e*16Ypov?W7n5?tEwbxlu?pnDO=*PZBIua*mv{qI^Y&+f_kp*{|{ z37mR3G=Y(Dg`+$sh)-I^soC6v$1o6YykKd(I-^PcQVGKJsd()r#{P zQPT3s>K-)1ad0mZp|Q?EIR=VH`HHS{&7PY8em!5^b4zWJw^)GGR3OIWMiJ^6R5=Uq zemL~k8?NR0T8G67iR-;EHDz&7$(?g~EXW#K@wV3-w9rRQIoRbnF+aN9I$+x<{4liP zgM;5H8F^}Fa{&W8>80jscksU-V0tx7Vgi&j3(FZ{Ywy0Af9V{`@-F0SG?mu!_m&a7 z+T9N=r1&SrU>F$8o&t*NDFWRlXTtVsSckAi0f)cbi%KBp_qHx^Bpau3`@koNfCs+p zgMVUW|C~9X!D;@-dv(3unrqNaF4l((s6n0TyH)+9nD5<5nBUS~yL@w>fv_E7>{lw7 z3g~)7@&%oZkd@hAO51r zg>8Q`4QS8x@eq)y$lHT?RORRX)oTDAR2J|&pK>0x*+T>rX_(HPtUh?d5e$MJ*-0-6 zS}HD1(L2E8r6v-aVDa<4x~h;EqssU?mVtxS#_Hh+4J{8_tcgt;rb@NUINny zQ8k2uZ zlr<$}pAqqjY&7h%J#WxVV$LU{ z(4f2EIGc`LwmW+8;6KV=xgn>imEKtRV!Fb0GbOTcMt>;U*|JV}KfKO7DN)wu0a;5E z;S`UrwWz1&QM0yK4|LrNa95tZasNNYcz$CR3cAuiYr_JU^gi6QJiOnmI-$^_b05ysY7Rj6K=)b)Bm^W`f~NF;%K~uJ$Ad zb*A)xVjzV>3e>(IpR5bn&@(b1`eam-rEYLY_0i-g>4p125|MXcbwu6#@g`uCwS z+=UUH2$96cZG3d_kpzF1o6|k(VVBF68GiZItIl%%d6Tfoe zCmfA=Ct-RgSz@$TW0HA#)2YXO(RO_&*_`zWO%%yCJFhC)j~;q#5-MEB;A-3DD&(a5 zKbduBi&j<+09V`jE73dARoA=0S#;LASlZktVYlFqb1Jj3NP=&mi}|wwX?0)h~gY)qR!`Y*DUyFR5pW@w2N;c!s*?s%L9-KT@AmioN|<1 z6K1bWX}WWxsxA2wk!@0Dl^&f^;fZlyniJ#3N0`sAQEa7Q0XF1JkFaD4;4Xs?c+yft*{9eQ2R#Hhk%R+^hz zkHq{=jY(tv?NA4Ytq<%6Xs|8wZ_#bgdR-+zK+k*3qp8RCH`q#Py6zD={(Wj-eVsXxGyDL@n=+!l0)Uyd7#NMx0Bd zj1tv905|!b-B(`U(H8Ef7!G@-cd}7l8V+P2P^zE_xJfB49PcHflf*gA-<{VyVZIaM zG{s7ZVl-eLW&Ijq0%co_y#e`d;?FA{WHo^1kitYaYKg&DoO;L+EC&`Fia*o6jN