From b9b721fa2b6a35fe8c155e24b6f8a1318e44483c Mon Sep 17 00:00:00 2001
From: Alec Bills <abills@andrew.cmu.edu>
Date: Thu, 31 Oct 2024 07:50:25 -0700
Subject: [PATCH 01/13] add coupled variable to expression tree and
 discretisation

---
 src/pybamm/__init__.py                        |  1 +
 src/pybamm/discretisations/discretisation.py  |  5 ++
 .../expression_tree/coupled_variable.py       | 54 +++++++++++++++++++
 3 files changed, 60 insertions(+)
 create mode 100644 src/pybamm/expression_tree/coupled_variable.py

diff --git a/src/pybamm/__init__.py b/src/pybamm/__init__.py
index 68529156e3..3de52e5724 100644
--- a/src/pybamm/__init__.py
+++ b/src/pybamm/__init__.py
@@ -39,6 +39,7 @@
 from .expression_tree.parameter import Parameter, FunctionParameter
 from .expression_tree.scalar import Scalar
 from .expression_tree.variable import *
+from .expression_tree.coupled_variable import *
 from .expression_tree.independent_variable import *
 from .expression_tree.independent_variable import t
 from .expression_tree.vector import Vector
diff --git a/src/pybamm/discretisations/discretisation.py b/src/pybamm/discretisations/discretisation.py
index af4bd2edd6..7255e4923a 100644
--- a/src/pybamm/discretisations/discretisation.py
+++ b/src/pybamm/discretisations/discretisation.py
@@ -938,6 +938,11 @@ def _process_symbol(self, symbol):
             if symbol._expected_size is None:
                 symbol._expected_size = expected_size
             return symbol.create_copy()
+        
+        elif isinstance(symbol, pybamm.CoupledVariable):
+            new_symbol = self.process_symbol(symbol.children[0])
+            return new_symbol
+
         else:
             # Backup option: return the object
             return symbol
diff --git a/src/pybamm/expression_tree/coupled_variable.py b/src/pybamm/expression_tree/coupled_variable.py
new file mode 100644
index 0000000000..2b1cd61be5
--- /dev/null
+++ b/src/pybamm/expression_tree/coupled_variable.py
@@ -0,0 +1,54 @@
+import pybamm
+
+from pybamm.type_definitions import DomainType
+
+
+class CoupledVariable(pybamm.Symbol):
+    """
+    A node in the expression tree representing a variable whose equation is set by a different model or submodel.
+
+    
+    Parameters
+    ----------
+    name : str
+        The variable's name. If the 
+    """
+    def __init__(
+        self,
+        name: str,
+        domain: DomainType = None,
+    ) -> None:
+        super().__init__(name, domain=domain)
+
+    
+    def _evaluate_for_shape(self):
+        """
+        Returns the scalar 'NaN' to represent the shape of a parameter.
+        See :meth:`pybamm.Symbol.evaluate_for_shape()`
+        """
+        return pybamm.evaluate_for_shape_using_domain(self.domains)
+
+
+    def create_copy(self):
+        """See :meth:`pybamm.Symbol.new_copy()`."""
+        new_input_parameter = CoupledVariable(
+            self.name, self.domain, expected_size=self._expected_size
+        )
+        return new_input_parameter
+
+    @property
+    def children(self):
+        return self._children
+
+    @children.setter
+    def children(self, expr):
+        self._children = expr
+
+
+    def set_coupled_variable(self, symbol, expr):
+        if self == symbol:
+            symbol.children = [expr,]
+        else:
+            for child in symbol.children:
+                self.set_coupled_variable(child, expr)
+        symbol.set_id()
\ No newline at end of file

From f30749187ae2c777e8c21ff5b45031610e735c31 Mon Sep 17 00:00:00 2001
From: Alec Bills <abills@andrew.cmu.edu>
Date: Thu, 31 Oct 2024 11:13:46 -0700
Subject: [PATCH 02/13] add test; add coupledvariable dict to model

---
 src/pybamm/models/base_model.py               | 24 ++++++
 .../test_coupled_variable.py                  | 76 +++++++++++++++++++
 2 files changed, 100 insertions(+)
 create mode 100644 tests/unit/test_expression_tree/test_coupled_variable.py

diff --git a/src/pybamm/models/base_model.py b/src/pybamm/models/base_model.py
index f6f47acc55..dfe3e838a7 100644
--- a/src/pybamm/models/base_model.py
+++ b/src/pybamm/models/base_model.py
@@ -56,6 +56,7 @@ def __init__(self, name="Unnamed model"):
         self._boundary_conditions = {}
         self._variables_by_submodel = {}
         self._variables = pybamm.FuzzyDict({})
+        self._coupled_variables = {}
         self._summary_variables = []
         self._events = []
         self._concatenated_rhs = None
@@ -182,6 +183,29 @@ def boundary_conditions(self):
     def boundary_conditions(self, boundary_conditions):
         self._boundary_conditions = BoundaryConditionsDict(boundary_conditions)
 
+    @property
+    def coupled_variables(self):
+        """Returns a dictionary mapping strings to expressions representing variables needed by the model but whose equations were set by other models."""
+        return self._coupled_variables
+    
+    @coupled_variables.setter
+    def coupled_variables(self, coupled_variables):
+        for name, var in coupled_variables.items():
+            if (
+                isinstance(var, pybamm.CoupledVariable)
+                and var.name != name
+                # Exception if the variable is also there under its own name
+                and not (var.name in coupled_variables and coupled_variables[var.name] == var)
+            ):
+                raise ValueError(
+                    f"Coupled variable with name '{var.name}' is in coupled variables dictionary with "
+                    f"name '{name}'. Names must match."
+                )
+        self._coupled_variables = coupled_variables
+
+    def list_coupled_variables(self):
+        list(self._coupled_variables.keys())
+
     @property
     def variables(self):
         """Returns a dictionary mapping strings to expressions representing the model's useful variables."""
diff --git a/tests/unit/test_expression_tree/test_coupled_variable.py b/tests/unit/test_expression_tree/test_coupled_variable.py
new file mode 100644
index 0000000000..48bf51c37c
--- /dev/null
+++ b/tests/unit/test_expression_tree/test_coupled_variable.py
@@ -0,0 +1,76 @@
+#
+# Tests for the CoupledVariable class
+#
+
+import pytest
+
+import numpy as np
+
+import pybamm
+
+def combine_models(list_of_models):
+    model = pybamm.BaseModel()
+    
+    for submodel in list_of_models:
+        model.coupled_variables.update(submodel.coupled_variables)
+        model.variables.update(submodel.variables)
+        model.rhs.update(submodel.rhs)
+        model.algebraic.update(submodel.algebraic)
+        model.initial_conditions.update(submodel.initial_conditions)
+        model.boundary_conditions.update(submodel.boundary_conditions)
+    
+    for name, coupled_variable in model.coupled_variables.items():
+        if name in model.variables:
+            for sym in model.rhs.values():
+                coupled_variable.set_coupled_variable(sym, model.variables[name])
+            for sym in model.algebraic.values():
+                coupled_variable.set_coupled_variable(sym, model.variables[name])
+    return model
+
+
+class TestCoupledVariable:
+    def test_coupled_variable(self):
+        model_1 = pybamm.BaseModel()
+        model_1_var_1 = pybamm.CoupledVariable("a")
+        model_1_var_2 = pybamm.Variable("b")
+        model_1.rhs[model_1_var_2] = -0.2 * model_1_var_1
+        model_1.variables["b"] = model_1_var_2
+        model_1.coupled_variables["a"] = model_1_var_1
+        model_1.initial_conditions[model_1_var_2] = 1.0
+
+        model_2 = pybamm.BaseModel()
+        model_2_var_1 = pybamm.Variable("a")
+        model_2_var_2 = pybamm.CoupledVariable("b")
+        model_2.rhs[model_2_var_1] = - 0.2 * model_2_var_2
+        model_2.variables["a"] = model_2_var_1
+        model_2.coupled_variables["b"] = model_2_var_2
+        model_2.initial_conditions[model_2_var_1] = 1.0
+
+        model = combine_models([model_1, model_2])
+
+        params = pybamm.ParameterValues({})
+        geometry = {}
+
+        # Process parameters
+        params.process_model(model)
+        params.process_geometry(geometry)
+
+        # mesh and discretise
+        submesh_types = {}
+        var_pts = {}
+        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
+
+
+        spatial_methods = {}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+
+        # solve
+        solver = pybamm.CasadiSolver()
+        t = np.linspace(0, 10, 1000)
+        solution = solver.solve(model, t)
+
+        np.testing.assert_almost_equal(solution["a"].entries, solution["b"].entries, decimal=10)
+
+
+

From 65a45aff486dc5492fb9a684b90542ebe62ac640 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 31 Oct 2024 18:18:40 +0000
Subject: [PATCH 03/13] style: pre-commit fixes

---
 src/pybamm/discretisations/discretisation.py     |  2 +-
 src/pybamm/expression_tree/coupled_variable.py   | 14 +++++++-------
 src/pybamm/models/base_model.py                  |  6 ++++--
 .../test_coupled_variable.py                     | 16 +++++++---------
 4 files changed, 19 insertions(+), 19 deletions(-)

diff --git a/src/pybamm/discretisations/discretisation.py b/src/pybamm/discretisations/discretisation.py
index 7255e4923a..2ca8c87649 100644
--- a/src/pybamm/discretisations/discretisation.py
+++ b/src/pybamm/discretisations/discretisation.py
@@ -938,7 +938,7 @@ def _process_symbol(self, symbol):
             if symbol._expected_size is None:
                 symbol._expected_size = expected_size
             return symbol.create_copy()
-        
+
         elif isinstance(symbol, pybamm.CoupledVariable):
             new_symbol = self.process_symbol(symbol.children[0])
             return new_symbol
diff --git a/src/pybamm/expression_tree/coupled_variable.py b/src/pybamm/expression_tree/coupled_variable.py
index 2b1cd61be5..14fd8fbcdd 100644
--- a/src/pybamm/expression_tree/coupled_variable.py
+++ b/src/pybamm/expression_tree/coupled_variable.py
@@ -7,12 +7,13 @@ class CoupledVariable(pybamm.Symbol):
     """
     A node in the expression tree representing a variable whose equation is set by a different model or submodel.
 
-    
+
     Parameters
     ----------
     name : str
-        The variable's name. If the 
+        The variable's name. If the
     """
+
     def __init__(
         self,
         name: str,
@@ -20,7 +21,6 @@ def __init__(
     ) -> None:
         super().__init__(name, domain=domain)
 
-    
     def _evaluate_for_shape(self):
         """
         Returns the scalar 'NaN' to represent the shape of a parameter.
@@ -28,7 +28,6 @@ def _evaluate_for_shape(self):
         """
         return pybamm.evaluate_for_shape_using_domain(self.domains)
 
-
     def create_copy(self):
         """See :meth:`pybamm.Symbol.new_copy()`."""
         new_input_parameter = CoupledVariable(
@@ -44,11 +43,12 @@ def children(self):
     def children(self, expr):
         self._children = expr
 
-
     def set_coupled_variable(self, symbol, expr):
         if self == symbol:
-            symbol.children = [expr,]
+            symbol.children = [
+                expr,
+            ]
         else:
             for child in symbol.children:
                 self.set_coupled_variable(child, expr)
-        symbol.set_id()
\ No newline at end of file
+        symbol.set_id()
diff --git a/src/pybamm/models/base_model.py b/src/pybamm/models/base_model.py
index dfe3e838a7..85111af7f7 100644
--- a/src/pybamm/models/base_model.py
+++ b/src/pybamm/models/base_model.py
@@ -187,7 +187,7 @@ def boundary_conditions(self, boundary_conditions):
     def coupled_variables(self):
         """Returns a dictionary mapping strings to expressions representing variables needed by the model but whose equations were set by other models."""
         return self._coupled_variables
-    
+
     @coupled_variables.setter
     def coupled_variables(self, coupled_variables):
         for name, var in coupled_variables.items():
@@ -195,7 +195,9 @@ def coupled_variables(self, coupled_variables):
                 isinstance(var, pybamm.CoupledVariable)
                 and var.name != name
                 # Exception if the variable is also there under its own name
-                and not (var.name in coupled_variables and coupled_variables[var.name] == var)
+                and not (
+                    var.name in coupled_variables and coupled_variables[var.name] == var
+                )
             ):
                 raise ValueError(
                     f"Coupled variable with name '{var.name}' is in coupled variables dictionary with "
diff --git a/tests/unit/test_expression_tree/test_coupled_variable.py b/tests/unit/test_expression_tree/test_coupled_variable.py
index 48bf51c37c..53056e9b25 100644
--- a/tests/unit/test_expression_tree/test_coupled_variable.py
+++ b/tests/unit/test_expression_tree/test_coupled_variable.py
@@ -2,15 +2,15 @@
 # Tests for the CoupledVariable class
 #
 
-import pytest
 
 import numpy as np
 
 import pybamm
 
+
 def combine_models(list_of_models):
     model = pybamm.BaseModel()
-    
+
     for submodel in list_of_models:
         model.coupled_variables.update(submodel.coupled_variables)
         model.variables.update(submodel.variables)
@@ -18,7 +18,7 @@ def combine_models(list_of_models):
         model.algebraic.update(submodel.algebraic)
         model.initial_conditions.update(submodel.initial_conditions)
         model.boundary_conditions.update(submodel.boundary_conditions)
-    
+
     for name, coupled_variable in model.coupled_variables.items():
         if name in model.variables:
             for sym in model.rhs.values():
@@ -41,7 +41,7 @@ def test_coupled_variable(self):
         model_2 = pybamm.BaseModel()
         model_2_var_1 = pybamm.Variable("a")
         model_2_var_2 = pybamm.CoupledVariable("b")
-        model_2.rhs[model_2_var_1] = - 0.2 * model_2_var_2
+        model_2.rhs[model_2_var_1] = -0.2 * model_2_var_2
         model_2.variables["a"] = model_2_var_1
         model_2.coupled_variables["b"] = model_2_var_2
         model_2.initial_conditions[model_2_var_1] = 1.0
@@ -60,7 +60,6 @@ def test_coupled_variable(self):
         var_pts = {}
         mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
 
-
         spatial_methods = {}
         disc = pybamm.Discretisation(mesh, spatial_methods)
         disc.process_model(model)
@@ -70,7 +69,6 @@ def test_coupled_variable(self):
         t = np.linspace(0, 10, 1000)
         solution = solver.solve(model, t)
 
-        np.testing.assert_almost_equal(solution["a"].entries, solution["b"].entries, decimal=10)
-
-
-
+        np.testing.assert_almost_equal(
+            solution["a"].entries, solution["b"].entries, decimal=10
+        )

From 4d638c9ebb984065c71c653c2f10d41bddcf53ad Mon Sep 17 00:00:00 2001
From: Alec Bills <abills@andrew.cmu.edu>
Date: Thu, 31 Oct 2024 11:21:20 -0700
Subject: [PATCH 04/13] pre-commit

merge
---
 src/pybamm/expression_tree/coupled_variable.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/pybamm/expression_tree/coupled_variable.py b/src/pybamm/expression_tree/coupled_variable.py
index 14fd8fbcdd..b85589cbe0 100644
--- a/src/pybamm/expression_tree/coupled_variable.py
+++ b/src/pybamm/expression_tree/coupled_variable.py
@@ -11,7 +11,9 @@ class CoupledVariable(pybamm.Symbol):
     Parameters
     ----------
     name : str
-        The variable's name. If the
+        name of the node
+    domain : iterable of str
+        list of domains that this variable is valid over
     """
 
     def __init__(

From 0d4f12d4e45ff1feef68e413f59d413abf1c2e7d Mon Sep 17 00:00:00 2001
From: Alec Bills <abills@andrew.cmu.edu>
Date: Thu, 31 Oct 2024 11:40:22 -0700
Subject: [PATCH 05/13] Trigger CI


From c259bb0e6046dc5fd591e39d487ba6cf340b2915 Mon Sep 17 00:00:00 2001
From: Alec Bills <abills@andrew.cmu.edu>
Date: Thu, 31 Oct 2024 12:36:51 -0700
Subject: [PATCH 06/13] add tests for coverage; valentin comments

---
 .../expression_tree/coupled_variable.py       | 11 +++++-----
 src/pybamm/models/base_model.py               |  2 +-
 .../test_coupled_variable.py                  | 20 +++++++++++++++++++
 3 files changed, 26 insertions(+), 7 deletions(-)

diff --git a/src/pybamm/expression_tree/coupled_variable.py b/src/pybamm/expression_tree/coupled_variable.py
index b85589cbe0..04d03d2792 100644
--- a/src/pybamm/expression_tree/coupled_variable.py
+++ b/src/pybamm/expression_tree/coupled_variable.py
@@ -13,7 +13,7 @@ class CoupledVariable(pybamm.Symbol):
     name : str
         name of the node
     domain : iterable of str
-        list of domains that this variable is valid over
+        list of domains that this coupled variable is valid over
     """
 
     def __init__(
@@ -31,11 +31,9 @@ def _evaluate_for_shape(self):
         return pybamm.evaluate_for_shape_using_domain(self.domains)
 
     def create_copy(self):
-        """See :meth:`pybamm.Symbol.new_copy()`."""
-        new_input_parameter = CoupledVariable(
-            self.name, self.domain, expected_size=self._expected_size
-        )
-        return new_input_parameter
+        """Creates a new copy of the coupled variable."""
+        new_coupled_variable = CoupledVariable(self.name, self.domain)
+        return new_coupled_variable
 
     @property
     def children(self):
@@ -46,6 +44,7 @@ def children(self, expr):
         self._children = expr
 
     def set_coupled_variable(self, symbol, expr):
+        """Sets the children of the coupled variable to the expression passed in expr. If the symbol is not the coupled variable, then it searches the children of the symbol for the coupled variable. The coupled variable will be replaced by its first child (symbol.children[0], which should be expr) in the discretisation step."""
         if self == symbol:
             symbol.children = [
                 expr,
diff --git a/src/pybamm/models/base_model.py b/src/pybamm/models/base_model.py
index 85111af7f7..f7e8f70f32 100644
--- a/src/pybamm/models/base_model.py
+++ b/src/pybamm/models/base_model.py
@@ -206,7 +206,7 @@ def coupled_variables(self, coupled_variables):
         self._coupled_variables = coupled_variables
 
     def list_coupled_variables(self):
-        list(self._coupled_variables.keys())
+        return list(self._coupled_variables.keys())
 
     @property
     def variables(self):
diff --git a/tests/unit/test_expression_tree/test_coupled_variable.py b/tests/unit/test_expression_tree/test_coupled_variable.py
index 53056e9b25..3e60c412e5 100644
--- a/tests/unit/test_expression_tree/test_coupled_variable.py
+++ b/tests/unit/test_expression_tree/test_coupled_variable.py
@@ -7,6 +7,8 @@
 
 import pybamm
 
+import pytest
+
 
 def combine_models(list_of_models):
     model = pybamm.BaseModel()
@@ -72,3 +74,21 @@ def test_coupled_variable(self):
         np.testing.assert_almost_equal(
             solution["a"].entries, solution["b"].entries, decimal=10
         )
+
+        assert set(model.list_coupled_variables()) == set(["a", "b"])
+
+    def test_create_copy(self):
+        a = pybamm.CoupledVariable("a")
+        b = a.create_copy()
+        assert a == b
+
+    def test_setter(self):
+        model = pybamm.BaseModel()
+        a = pybamm.CoupledVariable("a")
+        coupled_variables = {"a": a}
+        model.coupled_variables = coupled_variables
+        assert model.coupled_variables == coupled_variables
+
+        with pytest.raises(ValueError, match="Coupled variable with name"):
+            coupled_variables = {"b": a}
+            model.coupled_variables = coupled_variables

From 969a8796c3adb2bfe31ff79a921933ae29fd1623 Mon Sep 17 00:00:00 2001
From: Pradyot Ranjan <99216956+prady0t@users.noreply.github.com>
Date: Tue, 5 Nov 2024 02:53:15 +0530
Subject: [PATCH 07/13] Using `tempfile` for test_symbol_visualise to remove
 flaky test. (#4544)

* Using tempfile for test_symbol_visualise to remove falky test

Signed-off-by: Pradyot Ranjan <99216956+pradyotRanjan@users.noreply.github.com>

* style: pre-commit fixes

* Using tmp_path fixture

Signed-off-by: Pradyot Ranjan <99216956+pradyotRanjan@users.noreply.github.com>

* style: pre-commit fixes

---------

Signed-off-by: Pradyot Ranjan <99216956+pradyotRanjan@users.noreply.github.com>
Co-authored-by: Pradyot Ranjan <99216956+pradyotRanjan@users.noreply.github.com>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Eric G. Kratz <kratman@users.noreply.github.com>
---
 .../unit/test_expression_tree/test_symbol.py  | 24 +++++++++----------
 1 file changed, 11 insertions(+), 13 deletions(-)

diff --git a/tests/unit/test_expression_tree/test_symbol.py b/tests/unit/test_expression_tree/test_symbol.py
index a61d86cbe0..d7374e8da7 100644
--- a/tests/unit/test_expression_tree/test_symbol.py
+++ b/tests/unit/test_expression_tree/test_symbol.py
@@ -3,8 +3,6 @@
 #
 
 import pytest
-import os
-from tempfile import TemporaryDirectory
 
 import numpy as np
 from scipy.sparse import csr_matrix, coo_matrix
@@ -391,17 +389,17 @@ def test_symbol_repr(self):
             pybamm.grad(c).__repr__(),
         )
 
-    def test_symbol_visualise(self):
-        with TemporaryDirectory() as dir_name:
-            test_stub = os.path.join(dir_name, "test_visualize")
-            test_name = f"{test_stub}.png"
-            c = pybamm.Variable("c", "negative electrode")
-            d = pybamm.Variable("d", "negative electrode")
-            sym = pybamm.div(c * pybamm.grad(c)) + (c / d + c - d) ** 5
-            sym.visualise(test_name)
-            assert os.path.exists(test_name)
-            with pytest.raises(ValueError):
-                sym.visualise(test_stub)
+    @pytest.fixture(scope="session")
+    def test_symbol_visualise(self, tmp_path):
+        temp_file = tmp_path / "test_visualize.png"
+        c = pybamm.Variable("c", "negative electrode")
+        d = pybamm.Variable("d", "negative electrode")
+        sym = pybamm.div(c * pybamm.grad(c)) + (c / d + c - d) ** 5
+        sym.visualise(str(temp_file))
+        assert temp_file.exists()
+
+        with pytest.raises(ValueError):
+            sym.visualise(str(temp_file.with_suffix("")))
 
     def test_has_spatial_derivatives(self):
         var = pybamm.Variable("var", domain="test")

From 106c249d7beb890c9c6e9a6633a43e1bf9e51886 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Mon, 4 Nov 2024 21:47:42 +0000
Subject: [PATCH 08/13] chore: update pre-commit hooks (#4564)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

updates:
- [github.com/astral-sh/ruff-pre-commit: v0.7.1 → v0.7.2](https://github.com/astral-sh/ruff-pre-commit/compare/v0.7.1...v0.7.2)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
---
 .pre-commit-config.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index dfc47aad43..60cd57e988 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -4,7 +4,7 @@ ci:
 
 repos:
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: "v0.7.1"
+    rev: "v0.7.2"
     hooks:
       - id: ruff
         args: [--fix, --show-fixes]

From f4e4955e41932faacde6394de98f92646e036cfc Mon Sep 17 00:00:00 2001
From: Alec Bills <48105066+aabills@users.noreply.github.com>
Date: Mon, 4 Nov 2024 14:40:32 -0800
Subject: [PATCH 09/13] add reaction heating (#4557)

* add reaction heating

* style: pre-commit fixes

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Eric G. Kratz <kratman@users.noreply.github.com>
---
 src/pybamm/models/submodels/thermal/base_thermal.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/pybamm/models/submodels/thermal/base_thermal.py b/src/pybamm/models/submodels/thermal/base_thermal.py
index 42d90f1bcf..384cb1723a 100644
--- a/src/pybamm/models/submodels/thermal/base_thermal.py
+++ b/src/pybamm/models/submodels/thermal/base_thermal.py
@@ -143,8 +143,12 @@ def _get_standard_coupled_variables(self, variables):
             phase_names = ["primary ", "secondary "]
 
         if self.options.electrode_types["negative"] == "planar":
-            Q_rxn_n = pybamm.FullBroadcast(
-                0, ["negative electrode"], "current collector"
+            i_n = variables["Lithium metal total interfacial current density [A.m-2]"]
+            eta_r_n = variables["Lithium metal interface reaction overpotential [V]"]
+            Q_rxn_n = pybamm.PrimaryBroadcast(
+                i_n * eta_r_n / self.param.n.L,
+                ["negative electrode"],
+                "current collector",
             )
             Q_rev_n = pybamm.FullBroadcast(
                 0, ["negative electrode"], "current collector"

From b0382e6f1e08ef28987be852fe7fad2f94cee750 Mon Sep 17 00:00:00 2001
From: "Eric G. Kratz" <kratman@users.noreply.github.com>
Date: Mon, 4 Nov 2024 20:28:32 -0500
Subject: [PATCH 10/13] Cleanup formatting in error messages (#4565)

* Remove some unneeded links

* Cleaning up some formatting

* Update src/pybamm/discretisations/discretisation.py

* Update src/pybamm/expression_tree/operations/convert_to_casadi.py
---
 CONTRIBUTING.md                                    |  6 +++---
 src/pybamm/discretisations/discretisation.py       | 12 +++++-------
 src/pybamm/expression_tree/averages.py             |  8 ++++----
 src/pybamm/expression_tree/binary_operators.py     | 12 ++++++------
 .../operations/convert_to_casadi.py                |  6 ++----
 src/pybamm/expression_tree/unary_operators.py      |  4 ++--
 src/pybamm/expression_tree/vector.py               |  4 +---
 src/pybamm/meshes/one_dimensional_submeshes.py     |  4 ++--
 src/pybamm/meshes/scikit_fem_submeshes.py          | 10 +++++-----
 src/pybamm/models/base_model.py                    | 14 +++++---------
 src/pybamm/simulation.py                           | 14 ++++++--------
 src/pybamm/solvers/casadi_algebraic_solver.py      |  8 +++-----
 12 files changed, 44 insertions(+), 58 deletions(-)

diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index fd3426bd3b..eb510f7054 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -44,7 +44,7 @@ You now have everything you need to start making changes!
 
 ### B. Writing your code
 
-6. PyBaMM is developed in [Python](https://www.python.org)), and makes heavy use of [NumPy](https://numpy.org/) (see also [NumPy for MatLab users](https://numpy.org/doc/stable/user/numpy-for-matlab-users.html) and [Python for R users](https://www.rebeccabarter.com/blog/2023-09-11-from_r_to_python)).
+6. PyBaMM is developed in [Python](https://www.python.org), and makes heavy use of [NumPy](https://numpy.org/).
 7. Make sure to follow our [coding style guidelines](#coding-style-guidelines).
 8. Commit your changes to your branch with [useful, descriptive commit messages](https://chris.beams.io/posts/git-commit/): Remember these are
    publicly visible and should still make sense a few months ahead in time.
@@ -116,8 +116,8 @@ PyBaMM provides a utility function `import_optional_dependency`, to check for th
 
 Optional dependencies should never be imported at the module level, but always inside methods. For example:
 
-```
-def use_pybtex(x,y,z):
+```python
+def use_pybtex(x, y, z):
     pybtex = import_optional_dependency("pybtex")
     ...
 ```
diff --git a/src/pybamm/discretisations/discretisation.py b/src/pybamm/discretisations/discretisation.py
index 2ca8c87649..3d9579ff9c 100644
--- a/src/pybamm/discretisations/discretisation.py
+++ b/src/pybamm/discretisations/discretisation.py
@@ -500,8 +500,8 @@ def check_tab_conditions(self, symbol, bcs):
 
         if domain != "current collector":
             raise pybamm.ModelError(
-                f"""Boundary conditions can only be applied on the tabs in the domain
-            'current collector', but {symbol} has domain {domain}"""
+                "Boundary conditions can only be applied on the tabs in the domain "
+                f"'current collector', but {symbol} has domain {domain}"
             )
         # Replace keys with "left" and "right" as appropriate for 1D meshes
         if isinstance(mesh, pybamm.SubMesh1D):
@@ -893,11 +893,9 @@ def _process_symbol(self, symbol):
                 y_slices = self.y_slices[symbol]
             except KeyError as error:
                 raise pybamm.ModelError(
-                    f"""
-                    No key set for variable '{symbol.name}'. Make sure it is included in either
-                    model.rhs or model.algebraic in an unmodified form
-                    (e.g. not Broadcasted)
-                    """
+                    f"No key set for variable '{symbol.name}'. Make sure it is included in either "
+                    "model.rhs or model.algebraic in an unmodified form "
+                    "(e.g. not Broadcasted)"
                 ) from error
             # Add symbol's reference and multiply by the symbol's scale
             # so that the state vector is of order 1
diff --git a/src/pybamm/expression_tree/averages.py b/src/pybamm/expression_tree/averages.py
index 5fa6c5f00f..11538ea153 100644
--- a/src/pybamm/expression_tree/averages.py
+++ b/src/pybamm/expression_tree/averages.py
@@ -251,8 +251,8 @@ def z_average(symbol: pybamm.Symbol) -> pybamm.Symbol:
     # Symbol must have domain [] or ["current collector"]
     if symbol.domain not in [[], ["current collector"]]:
         raise pybamm.DomainError(
-            f"""z-average only implemented in the 'current collector' domain,
-            but symbol has domains {symbol.domain}"""
+            "z-average only implemented in the 'current collector' domain, "
+            f"but symbol has domains {symbol.domain}"
         )
     # If symbol doesn't have a domain, its average value is itself
     if symbol.domain == []:
@@ -285,8 +285,8 @@ def yz_average(symbol: pybamm.Symbol) -> pybamm.Symbol:
     # Symbol must have domain [] or ["current collector"]
     if symbol.domain not in [[], ["current collector"]]:
         raise pybamm.DomainError(
-            f"""y-z-average only implemented in the 'current collector' domain,
-            but symbol has domains {symbol.domain}"""
+            "y-z-average only implemented in the 'current collector' domain, "
+            f"but symbol has domains {symbol.domain}"
         )
     # If symbol doesn't have a domain, its average value is itself
     if symbol.domain == []:
diff --git a/src/pybamm/expression_tree/binary_operators.py b/src/pybamm/expression_tree/binary_operators.py
index 1d630887b2..42f2f74e24 100644
--- a/src/pybamm/expression_tree/binary_operators.py
+++ b/src/pybamm/expression_tree/binary_operators.py
@@ -36,7 +36,7 @@ def _preprocess_binary(
     # Check both left and right are pybamm Symbols
     if not (isinstance(left, pybamm.Symbol) and isinstance(right, pybamm.Symbol)):
         raise NotImplementedError(
-            f"""BinaryOperator not implemented for symbols of type {type(left)} and {type(right)}"""
+            f"BinaryOperator not implemented for symbols of type {type(left)} and {type(right)}"
         )
 
     # Do some broadcasting in special cases, to avoid having to do this manually
@@ -389,9 +389,9 @@ def _binary_jac(self, left_jac, right_jac):
             return left @ right_jac
         else:
             raise NotImplementedError(
-                f"""jac of 'MatrixMultiplication' is only
-             implemented for left of type 'pybamm.Array',
-             not {left.__class__}"""
+                f"jac of 'MatrixMultiplication' is only "
+                "implemented for left of type 'pybamm.Array', "
+                f"not {left.__class__}"
             )
 
     def _binary_evaluate(self, left, right):
@@ -1541,8 +1541,8 @@ def source(
 
     if left.domain != ["current collector"] or right.domain != ["current collector"]:
         raise pybamm.DomainError(
-            f"""'source' only implemented in the 'current collector' domain,
-            but symbols have domains {left.domain} and {right.domain}"""
+            "'source' only implemented in the 'current collector' domain, "
+            f"but symbols have domains {left.domain} and {right.domain}"
         )
     if boundary:
         return pybamm.BoundaryMass(right) @ left
diff --git a/src/pybamm/expression_tree/operations/convert_to_casadi.py b/src/pybamm/expression_tree/operations/convert_to_casadi.py
index 274fd95154..af53b9acd8 100644
--- a/src/pybamm/expression_tree/operations/convert_to_casadi.py
+++ b/src/pybamm/expression_tree/operations/convert_to_casadi.py
@@ -231,8 +231,6 @@ def _convert(self, symbol, t, y, y_dot, inputs):
 
         else:
             raise TypeError(
-                f"""
-                Cannot convert symbol of type '{type(symbol)}' to CasADi. Symbols must all be
-                'linear algebra' at this stage.
-                """
+                f"Cannot convert symbol of type '{type(symbol)}' to CasADi. Symbols must all be "
+                "'linear algebra' at this stage."
             )
diff --git a/src/pybamm/expression_tree/unary_operators.py b/src/pybamm/expression_tree/unary_operators.py
index aa90fd6f4c..f41897e2de 100644
--- a/src/pybamm/expression_tree/unary_operators.py
+++ b/src/pybamm/expression_tree/unary_operators.py
@@ -992,8 +992,8 @@ def __init__(self, name, child, side):
         if side in ["negative tab", "positive tab"]:
             if child.domain[0] != "current collector":
                 raise pybamm.ModelError(
-                    f"""Can only take boundary value on the tabs in the domain
-                'current collector', but {child} has domain {child.domain[0]}"""
+                    "Can only take boundary value on the tabs in the domain "
+                    f"'current collector', but {child} has domain {child.domain[0]}"
                 )
         self.side = side
         # boundary value of a child takes the primary domain from secondary domain
diff --git a/src/pybamm/expression_tree/vector.py b/src/pybamm/expression_tree/vector.py
index 6dc358afb0..e9067a4ffd 100644
--- a/src/pybamm/expression_tree/vector.py
+++ b/src/pybamm/expression_tree/vector.py
@@ -29,9 +29,7 @@ def __init__(
             entries = entries[:, np.newaxis]
         if entries.shape[1] != 1:
             raise ValueError(
-                f"""
-                Entries must have 1 dimension or be column vector, not have shape {entries.shape}
-                """
+                f"Entries must have 1 dimension or be column vector, not have shape {entries.shape}"
             )
         if name is None:
             name = f"Column vector of length {entries.shape[0]!s}"
diff --git a/src/pybamm/meshes/one_dimensional_submeshes.py b/src/pybamm/meshes/one_dimensional_submeshes.py
index 8f27049411..027c6d0421 100644
--- a/src/pybamm/meshes/one_dimensional_submeshes.py
+++ b/src/pybamm/meshes/one_dimensional_submeshes.py
@@ -297,8 +297,8 @@ def __init__(self, lims, npts, edges=None):
 
         if (npts + 1) != len(edges):
             raise pybamm.GeometryError(
-                f"""User-suppled edges has should have length (npts + 1) but has length
-                {len(edges)}.Number of points (npts) for domain {spatial_var.domain} is {npts}.""".replace(
+                "User-suppled edges has should have length (npts + 1) but has length "
+                f"{len(edges)}.Number of points (npts) for domain {spatial_var.domain} is {npts}.".replace(
                     "\n                ", " "
                 )
             )
diff --git a/src/pybamm/meshes/scikit_fem_submeshes.py b/src/pybamm/meshes/scikit_fem_submeshes.py
index ba624c7f48..2d769f3ca2 100644
--- a/src/pybamm/meshes/scikit_fem_submeshes.py
+++ b/src/pybamm/meshes/scikit_fem_submeshes.py
@@ -92,8 +92,8 @@ def read_lims(self, lims):
         # check coordinate system agrees
         if spatial_vars[0].coord_sys != spatial_vars[1].coord_sys:
             raise pybamm.DomainError(
-                f"""spatial variables should have the same coordinate system,
-                but have coordinate systems {spatial_vars[0].coord_sys} and {spatial_vars[1].coord_sys}"""
+                "spatial variables should have the same coordinate system, "
+                f"but have coordinate systems {spatial_vars[0].coord_sys} and {spatial_vars[1].coord_sys}"
             )
         return spatial_vars, tabs
 
@@ -360,9 +360,9 @@ def __init__(self, lims, npts, y_edges=None, z_edges=None):
             # check that npts equals number of user-supplied edges
             if npts[var.name] != len(edges[var.name]):
                 raise pybamm.GeometryError(
-                    f"""User-suppled edges has should have length npts but has length {len(edges[var.name])}.
-                     Number of points (npts) for variable {var.name} in
-                     domain {var.domain} is {npts[var.name]}."""
+                    f"User-supplied edges has should have length npts but has length {len(edges[var.name])}. "
+                    f"Number of points (npts) for variable {var.name} in "
+                    f"domain {var.domain} is {npts[var.name]}."
                 )
 
             # check end points of edges agree with spatial_lims
diff --git a/src/pybamm/models/base_model.py b/src/pybamm/models/base_model.py
index f7e8f70f32..6cb3af2fd7 100644
--- a/src/pybamm/models/base_model.py
+++ b/src/pybamm/models/base_model.py
@@ -1110,7 +1110,7 @@ def check_ics_bcs(self):
         for var in self.rhs.keys():
             if var not in self.initial_conditions.keys():
                 raise pybamm.ModelError(
-                    f"""no initial condition given for variable '{var}'"""
+                    f"no initial condition given for variable '{var}'"
                 )
 
     def check_variables(self):
@@ -1132,11 +1132,9 @@ def check_variables(self):
         for var in all_vars:
             if var not in vars_in_keys:
                 raise pybamm.ModelError(
-                    f"""
-                    No key set for variable '{var}'. Make sure it is included in either
-                    model.rhs or model.algebraic, in an unmodified form
-                    (e.g. not Broadcasted)
-                    """
+                    f"No key set for variable '{var}'. Make sure it is included in either "
+                    "model.rhs or model.algebraic, in an unmodified form "
+                    "(e.g. not Broadcasted)"
                 )
 
     def check_no_repeated_keys(self):
@@ -1550,9 +1548,7 @@ def check_and_convert_bcs(self, boundary_conditions):
                 # Check types
                 if bc[1] not in ["Dirichlet", "Neumann"]:
                     raise pybamm.ModelError(
-                        f"""
-                        boundary condition types must be Dirichlet or Neumann, not '{bc[1]}'
-                        """
+                        f"boundary condition types must be Dirichlet or Neumann, not '{bc[1]}'"
                     )
 
         return boundary_conditions
diff --git a/src/pybamm/simulation.py b/src/pybamm/simulation.py
index da0ac08316..ab22972741 100644
--- a/src/pybamm/simulation.py
+++ b/src/pybamm/simulation.py
@@ -519,14 +519,12 @@ def solve(
                     dt_eval_max = np.max(np.diff(t_eval))
                     if dt_eval_max > np.nextafter(dt_data_min, np.inf):
                         warnings.warn(
-                            f"""
-                            The largest timestep in t_eval ({dt_eval_max}) is larger than
-                            the smallest timestep in the data ({dt_data_min}). The returned
-                            solution may not have the correct resolution to accurately
-                            capture the input. Try refining t_eval. Alternatively,
-                            passing t_eval = None automatically sets t_eval to be the
-                            points in the data.
-                            """,
+                            f"The largest timestep in t_eval ({dt_eval_max}) is larger than "
+                            f"the smallest timestep in the data ({dt_data_min}). The returned "
+                            "solution may not have the correct resolution to accurately "
+                            "capture the input. Try refining t_eval. Alternatively, "
+                            "passing t_eval = None automatically sets t_eval to be the "
+                            "points in the data.",
                             pybamm.SolverWarning,
                             stacklevel=2,
                         )
diff --git a/src/pybamm/solvers/casadi_algebraic_solver.py b/src/pybamm/solvers/casadi_algebraic_solver.py
index 2dd6f2d341..b139199f8c 100644
--- a/src/pybamm/solvers/casadi_algebraic_solver.py
+++ b/src/pybamm/solvers/casadi_algebraic_solver.py
@@ -146,11 +146,9 @@ def _integrate(self, model, t_eval, inputs_dict=None, t_interp=None):
                 )
             else:
                 raise pybamm.SolverError(
-                    f"""
-                    Could not find acceptable solution: solver terminated
-                    successfully, but maximum solution error ({casadi.mmax(casadi.fabs(fun))})
-                    above tolerance ({self.tol})
-                    """
+                    "Could not find acceptable solution: solver terminated "
+                    f"successfully, but maximum solution error ({casadi.mmax(casadi.fabs(fun))}) "
+                    f"above tolerance ({self.tol})"
                 )
 
         # Concatenate differential part

From 4a505629d9497723fc2642735f884e066f155876 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Tue, 5 Nov 2024 11:46:49 +0000
Subject: [PATCH 11/13] docs: add sensitivities notebook (#4559)

* docs: add sensitivities notebook

* docs: add sens notebook to toctree

* add install line

---------

Co-authored-by: Eric G. Kratz <kratman@users.noreply.github.com>
---
 docs/source/examples/index.rst                |   1 +
 .../sensitivities_and_data_fitting.ipynb      | 327 ++++++++++++++++++
 2 files changed, 328 insertions(+)
 create mode 100644 docs/source/examples/notebooks/parameterization/sensitivities_and_data_fitting.ipynb

diff --git a/docs/source/examples/index.rst b/docs/source/examples/index.rst
index 3f77578ef5..6ddaf5867e 100644
--- a/docs/source/examples/index.rst
+++ b/docs/source/examples/index.rst
@@ -87,6 +87,7 @@ The notebooks are organised into subfolders, and can be viewed in the galleries
     notebooks/parameterization/change-input-current.ipynb
     notebooks/parameterization/parameter-values.ipynb
     notebooks/parameterization/parameterization.ipynb
+    notebooks/parameterization/sensitivities_and_data_fitting.ipynb
 
 .. nbgallery::
     :caption: Simulations and Experiments
diff --git a/docs/source/examples/notebooks/parameterization/sensitivities_and_data_fitting.ipynb b/docs/source/examples/notebooks/parameterization/sensitivities_and_data_fitting.ipynb
new file mode 100644
index 0000000000..09b562c1ca
--- /dev/null
+++ b/docs/source/examples/notebooks/parameterization/sensitivities_and_data_fitting.ipynb
@@ -0,0 +1,327 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "bd9929be",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Note: you may need to restart the kernel to use updated packages.\n"
+     ]
+    }
+   ],
+   "source": [
+    "%pip install \"pybamm[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
+    "# import dependencies\n",
+    "import pybamm\n",
+    "import numpy as np\n",
+    "import matplotlib.pylab as plt\n",
+    "import scipy.optimize"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b1223d98",
+   "metadata": {},
+   "source": [
+    "# Sensitivities and data fitting using PyBaMM\n",
+    "\n",
+    "PyBaMM input parameters [`pybamm.InputParameter`](https://docs.pybamm.org/en/stable/source/api/expression_tree/input_parameter.html) can be used to run many simulations with varying parameters. Here we will demonstrate PyBaMM's ability to calculate the senstivities of model outputs with respect to input parameters. \n",
+    "\n",
+    "To be more specific, given a model output $f(a)$, where $a$ is an input parameter, we wish to calculate $\\frac{\\partial f}{\\partial a}(a)$.\n",
+    "\n",
+    "First we will demonstrate using a toy model, given by the equations\n",
+    "\n",
+    "$$\\frac{dy}{dt} = a y$$\n",
+    "\n",
+    "with a scalar state variable $y$ and a scalar parameter $a$, and initial conditions\n",
+    "\n",
+    "$$y(0) = 1$$\n",
+    "\n",
+    "We will also define a model output given by $f = y^2$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "25970cdf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# setup a simple test model\n",
+    "model = pybamm.BaseModel(\"name\")\n",
+    "y = pybamm.Variable(\"y\")\n",
+    "a = pybamm.InputParameter(\"a\")\n",
+    "model.rhs = {y: a * y}\n",
+    "model.initial_conditions = {y: 1}\n",
+    "model.variables = {\"y squared\": y**2}\n",
+    "\n",
+    "solver = pybamm.IDAKLUSolver(rtol=1e-10, atol=1e-10)\n",
+    "t_eval = np.linspace(0, 1, 80)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9f3d61bf",
+   "metadata": {},
+   "source": [
+    "Note that we have used the [`pybamm.IDAKLUSolver`](https://docs.pybamm.org/en/stable/source/api/solvers/idaklu_solver.html) solver for this example, this is currently the recommended solver for calculating sensitivities in PyBaMM.\n",
+    "\n",
+    "We can solve the model using a specific value of $a = 1$. However, now we will also calculate the forward sensitivities of the model by setting the argument `calculate_sensitivies=True`. Note that this argument will also accept a list of input parameters to calculate the sensitivities for, but setting it to `True` will calculate the sensitivities for **all** input parameters of the model "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "1b1781a9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "solution = solver.solve(model, [0, 1], inputs={\"a\": 1}, calculate_sensitivities=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4d6f176e",
+   "metadata": {},
+   "source": [
+    "We can now access the solution as normal, and the sensitivities using the syntax: `solution[output_name].sensitivities[input_parameter_name]`\n",
+    "\n",
+    "Note that a current restriction to the sensitivity calculation is that it will only return the sensitivities at the values of `t_eval` used to solve the model. Any interpolation between these values will have to be done manually"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "bf0a2d9c",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAHWCAYAAAARl3+JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABZNUlEQVR4nO3deViU5f4/8PcMMMM+CLIKCIIriuKaS6lpLqnZnmVmtlqe0+KpY55TlqdFW34dO+XR6phaarZ81WxRU3MtVxQFFwRBZBEQlBnWYZi5f38MTJGorPM888z7dV1zFcMz8n4cuflwryohhAAREREROTS11AGIiIiIqOVY1BEREREpAIs6IiIiIgVgUUdERESkACzqiIiIiBSARR0RERGRArCoIyIiIlIAFnVERERECuAqdYC2ZLFYkJeXBx8fH6hUKqnjEFEjCSFQWlqKsLAwqNX83bOp2PYROa6WtH+KLury8vIQEREhdQwiaqbs7GyEh4dLHcPhsO0jcnzNaf8UXdT5+PgAsP7F+Pr6SpyGiBrLYDAgIiLC9j1MTcO2j8hxtaT9U3RRVzfs4Ovry4aNyAFx6LB52PYROb7mtH+crEJERESkACzqiIiIiBSARR0RERGRArCoIyIiIlIAFnVERERECsCijoiIiEgBWNQRERERKQCLOiIiIiIFYFFHREREpAAs6oiIiIgUgEUdERERkQKwqCMiIiJSABZ1RERERArAoo6IiIhIAVjUEZFdfH0oG2/+eBIpuXqpoxAR2dXcdcexeEc69BWmNv06rm36pxMR1Vp3NAf7My4hqr0XenbQSR2HiMguisuM+PJgNgDgocEd2/RrsaeOiNpclcmMI1klAIAhMe2lDUNEZEfJtaMTnQK94OPu1qZfi0UdEbW5xKzLqDZbEKpzR1SAp9RxiIjsJjnHWtT1ssMIBYs6Impzv50tAgAMjgmASqWSOA0Rkf3U9dSxqCMiRfjtbDEADr0SkfOpK+riw/3a/GuxqCOiNlVaZcLx2uGHwTEBEqchIrKfi6VGXNBXQaUC4sJ82/zrsagjojZ16NwlmC0CUQGe6ODnIXUcIiK7qdvCKSbQG17att9whEUdEbWp39KtQ6/spSMiZ1M3ShFvp22cWNQRUZuqm083mPPpiMjJJOeWAIDd9uZkUUdEbeZyeTVOXjAAAAZ3Yk8dETmX3xdJsKgjIge3P8PaS9cl2BuBPlqJ0xAR2U+hoQoFBiPUKqCHHRZJACzqiKgN/Vq7Px23MiEiZ1PXSxcb5A1PjX1OZWVRR0Rt5tfaRRJDY1nUEZFzOW47ScLPbl+TRR0RtYnckkpkFpXDRa3CoE7+UschIrKr30+SsM/QK8CijojayK/p1qHX+HAdfNv4EGsiIjkRQvxe1NnhJIk6LOqIqE3UFXXDOPR6hd27d2PSpEkICwuDSqXChg0brnrtzJkzoVKpsGjRIrvlI6KWKTAYcbHUCBe1Cj1C2VNHRA5MCMH5dNdQXl6O3r17Y/Hixde8bv369di/fz/CwsLslIyIWsPxnBIAQOcgb3hoXOz2de2zHIOInMqZgjIUlRnh7qZGQqSf1HFkZ/z48Rg/fvw1r8nNzcVf//pXbNmyBRMmTLBTMiJqDSm2+XT22Z+uDnvqiKjV7a0deh0YHQCtq/1+S1UKi8WCadOm4cUXX0RcXJzUcYioiY7bedPhOrIt6qKioqBSqa54zJo1S+poRHQdv8+n4ykSzfH222/D1dUVzzzzTKOuNxqNMBgM9R5EJA0hhK2nzl7Hg9WR7fDroUOHYDabbR+npKTglltuwT333CNhKiK6HpPZggMZnE/XXImJifjggw9w5MgRqFSqRr1mwYIFmD9/fhsnI6LGuKCvQlFZNVzVKnS34yIJQMY9dYGBgQgJCbE9fvjhB8TExGD48OFSRyOia0jKLkF5tRntPN3QPcS+DZoS7NmzB4WFhYiMjISrqytcXV2RlZWFv/3tb4iKimrwNXPnzoVer7c9srOz7RuaiGzqNh3uEuwDdzf7Tj+RbU/dH1VXV2PVqlWYPXv2NX9zNRqNMBqNto85BEFkf3vSrEOvQ2PbQ61uXE8T/W7atGkYPXp0vefGjh2LadOmYcaMGQ2+RqvVQqvl2bpEcpCcWwLA/oskAAcp6jZs2ICSkhI8/PDD17yOQxBE0tubdhEAcGNnDr1eTVlZGdLT020fZ2ZmIikpCf7+/oiMjERAQP25iG5ubggJCUHXrl3tHZWImigpuwQA0DvCz+5fW7bDr3+0bNkyjB8//rp7NXEIgkha+kqTrUEb1jlQ2jAydvjwYSQkJCAhIQEAMHv2bCQkJGDevHkSJyOilrBYBI5nW4dfe0ewp+4KWVlZ2LZtG9atW3fdazkEQSStfWeLYRFAp0AvdPDzkDqObI0YMQJCiEZff+7cubYLQ0StJqOoDKXGGni4uaBrsI/dv77se+qWL1+OoKAgbr5J5AD21A693sReOiJyQknZv2867Opi/xJL1kWdxWLB8uXLMX36dLi6yr5Tkcjp7eV5r0TkxJKyLwOQZugVkHlRt23bNpw/fx6PPPKI1FGI6DrOF1cgq7gCrmoVbojhpsNE5HyO1fbU9YloJ8nXl3X315gxY5o074SIpLMn3Tr02jeyHby1sm5aiIhaXZXJjFMXrFupsaeOiBza3tr96YZxKxMickIn8vSosQi099ZKtlCMRR0RtViN2WI775X70xGRM0qyDb3qGn3EX2tjUUdELXYsRw9DVQ10Hm6ID/eTOg4Rkd3V7dHZR4JNh+uwqCOiFtt9xjqfblhse7jwaDAickLHJDxJog6LOiJqsV21Rd1NXTj0SkTOp7jMiPOXKgBA0tEKFnVE1CIlFdU4nlMCALipCzcdJiLnczzHOp8uJtALOg83yXKwqCOiFtmbXgSLALoEeyNUx6PBiMj5HJXB0CvAoo6IWqhuPh2PBiMiZ1W3SCKBRR0ROSohBHafsW5lwqFXInJGQghZLJIAWNQRUQucKShDvqEKWlc1Bkb7Sx2HiMjuzhVXQF9pgsZVjW4hvpJmYVFHRM1WN/R6Q6cAuLu5SJyGiMj+6nrp4sJ8oXGVtqxiUUdEzbY7rW4rEw69EpFzksOmw3VY1BFRs1RU1+BAxiUAwHAWdUTkpFjUEZHD23e2GNVmC8LbeSAm0EvqOEREdmesMeNkngEAizoicmA7U61DryO6Bkp2eDURkZRO5BlQbbagnacbIv09pY7Doo6Imk4IgZ1nCgEAI7oESZyGiEgaR7IuAwD6dWwni19uWdQRUZNlFJUj+1IlNC5qDIkNkDoOEZEkjpy3FnV9O7aTOIkVizoiarK6odeB0f7w1LhKnIaIyP6EEEis66mLZFFHRA5qZ2rt0GtXrnolIueUW1KJAoMRrmoV4sP9pI4DgEUdETVRZbUZBzKtW5mwqCMiZ3XkfAkAoEeYLzw08th8nUUdETXJvowiVNdY0MHPAzGB3lLHISKSRN0iib4yGXoFWNQRURPtOM2tTIiI6ubTyWWRBMCijoiaQAiBX05b59Pd3I1bmRCRc6qorsHJC9ZNh/uxqCMiR5RWWIbckkpoXdUYEtNe6jhERJI4nqOH2SIQ4uuOMJ271HFsWNQRUaPV9dINjgmQzcRgIiJ7S5TZpsN1WNQRUaNx6JWI6PdFEgmRftIG+RMWdUTUKPoKk+2305FdWdQRkXMSQthOkpDTfDqARR0RNdLutIswWwRig7wRIYODq4mIpJBZVI7LFSZoXNWIC9NJHaceFnVE1Cg7OPRKRGQbsYjvoIPGVV5llLzSEJEsmS0CO89Y96fj0CsROTO5Dr0CLOqIqBGO5ZTgUnk1fNxd0T9Kfg0ZEZG9HMkqASCvTYfrsKgjouv65ZR16PWmLoFwc2GzQUTOSV9pwpnCUgDyOh6sDltnIrqubacKAACju3PolYicV1J2CYQAIv09EeijlTrOFVjUEdE15VyuwOn8UqhVwIguLOqIyHn9cdNhOWJRR0TXVLfhcL+O7dDOSyNxGiIi6SRmXQIA9JXZpsN1WNQR0TVtq51PN6p7sMRJlGP37t2YNGkSwsLCoFKpsGHDBtvnTCYT5syZg169esHLywthYWF46KGHkJeXJ11gIoLJbLEtkhgQ7S9tmKtgUUdEV1VmrMH+s8UAOJ+uNZWXl6N3795YvHjxFZ+rqKjAkSNH8Morr+DIkSNYt24dUlNTcdttt0mQlIjqnMgzoNJkhs7DDV2CfKSO0yBXqQMQkXztTStCtdmCjgGeiAn0ljqOYowfPx7jx49v8HM6nQ5bt26t99xHH32EgQMH4vz584iMjLRHRCL6k0OZ1qHXAVHtoFarJE7TMPbUEdFVba9d9TqqWzBUKnk2Ys5Ar9dDpVLBz89P6ihETuuAraiT59ArwJ46IroKs0VgR6p1Ph2HXqVTVVWFOXPm4P7774evr2+D1xiNRhiNRtvHBoPBXvGInILFInC4dpGEXOfTAeypI6KrSMq+jKKyulMk5NuIKZnJZMK9994LIQSWLFly1esWLFgAnU5ne0RERNgxJZHypV8sQ0mFCe5uavQM00kd56pY1BFRg7aetPbSjewaJLtDq51BXUGXlZWFrVu3XrWXDgDmzp0LvV5ve2RnZ9sxKZHyHawdek2IaCfr9pDDr0TUoK0n8wEAt/TgVib2VlfQpaWlYceOHQgICLjm9VqtFlqt/Ha3J1KKQ+fkP/QKsKgjogZkXCzD2YvlcHNRYXjXQKnjKE5ZWRnS09NtH2dmZiIpKQn+/v4IDQ3F3XffjSNHjuCHH36A2WxGfr61wPb394dGww2gieytbuXrIBZ1RORotp60rnq9oVMAfN3dJE6jPIcPH8bIkSNtH8+ePRsAMH36dLz22mvYuHEjAKBPnz71Xrdjxw6MGDHCXjGJCNajEvP0VXBVq5Ag05Mk6rCoI6IrbKvdyoRDr21jxIgREEJc9fPX+hwR2Vfd0GtcBx08NfIum+Q72w9Abm4uHnzwQQQEBMDDwwO9evXC4cOHpY5FpGjFZUbbodWjeTQYETm5g5nW9nBgVDuJk1yfbEvOy5cvY+jQoRg5ciQ2bdqEwMBApKWloV07+f+lEjmy7acLYRFAzw6+CPPzkDoOEZGkDmZaj0qU86bDdWRb1L399tuIiIjA8uXLbc9FR0dLmIjIOdTNp2MvHRE5u+IyI85eLAfgGEWdbIdfN27ciP79++Oee+5BUFAQEhIS8Omnn0odi0jRKqvN2JN2EQAwpkeIxGmIiKR16Jx16LVzkDfaecl/5blsi7qMjAwsWbIEnTt3xpYtW/DUU0/hmWeewcqVK6/6GqPRCIPBUO9BRI2368xFVJksiPD3QPdQH6njEBFJylH2p6sj2+FXi8WC/v3746233gIAJCQkICUlBUuXLsX06dMbfM2CBQswf/58e8YkUpSfT1j3QxvTIwQqlUriNERE0qor6uS+P10d2fbUhYaGokePHvWe6969O86fP3/V1/CoHKLmM5kttq1MxsZx6JWInFu5sQYn8qwjfo4wnw6QcU/d0KFDkZqaWu+5M2fOoGPHjld9DY/KIWq+AxmXYKiqQYCXBv06cpU5ETm3xKzLMFsEOvh5OMxOALLtqXv++eexf/9+vPXWW0hPT8eaNWvwySefYNasWVJHI1KkLSd+P+vVRc2hVyJybvsyrFuZ3NDp2mcvy4lsi7oBAwZg/fr1+PLLL9GzZ0+8/vrrWLRoEaZOnSp1NCLFsViEbSuTMXHcyoSIaL+tqHOMoVdAxsOvADBx4kRMnDhR6hhEinc8V498QxW8NC4YEtNe6jhERJIqM9bgeI4eADA4hj11RORA6oZeR3QLgrubi8RpiIikdejcJZgtAhH+Hghv5yl1nEZjUUfk5IQQ2JxiLerGcdUrEdHvQ6/RjtNLB7CoI3J6qQWlyCwqh8ZVjZHdgqSOQ0Qkuf1nrUWdIw29AizqiJzepmRrL91NnQPhrZX1NFsiojZnqDIhOdc6n86RVr4CLOqInF7d0Ov4nhx6JSI6lHkJFgF0DPB0mP3p6rCoI3JiGRfLkFpQCle1CqO7cysTIqK6+XSDHayXDmBRR+TUNtX20g2OCYDO003iNERE0qvbdNjR5tMBLOqInNrvQ6+hEichIpKevsJkO+/V0ebTASzqiJxWzuUKJOfqoVbxFAkiIgA4eO4ShAA6tfdCsK+71HGajEUdkZOqW/U6IMof7b21EqchIpLevtqtTG5wwKFXgEUdkdP6MfkCAGBCPIdeiYiA3+fTOeLQK8Cijsgp5ZZUIim7BCoVMI5bmRARoaSiGqfz6+bT+UucpnlY1BE5oU21vXQDovwR5ON480aIiFrb/gzrfLrYIG+HbRdZ1BE5IdvQay8OvRIRAX8479VBe+kAFnVETievpBJHz1uHXnmKBBGRVd0iicGd2kucpPlY1BE5mZ/qhl47+iPIAZfsExG1tsLSKqQWlAJwzE2H67CoI3IydUXdrb3YS0dEBAC/phcBAOLCfOHvpZE4TfOxqCNyInkllThSN/TK+XRERACAvWnWoddhnR136BVgUUfkVH48/vuqV0fcLZ2IqLUJIbA3/SIAYFgsizoichA/HM8DAEzihsNERACAsxfLUGAwQuOqxoAox135CrCoI3Ia54srcCzHetbruJ4s6oiIAGBvmnU+3YCodnB3c5E4TcuwqCNyEj8kW3vpBscEINCHZ70SEQHA3tpFEsNiAyVO0nIs6oicxA/HrPPpJsaHSZyEiEgeTGYL9mdcAgDc6OCLJAAWdURO4ezFMpy8YICrWoVxcdzKhIgIAI5ll6DMWIN2nm7oEeordZwWY1FH5ATqeumGdW6Pdg68BxMRUWuqG3odEtsearVK4jQtx6KOSOGEEPi+dtUrh16JiH5Xt0jC0bcyqcOijkjhTl0oRXphGTSuaoyJC5Y6DhGRLJRWmXA0uwQAizoichAbj1l76W7uGgRfdzeJ0xARycOBjEswWwQ6Bngiwt9T6jitgkUdkYJZLALf1xZ1t/Xh0CsRUZ3ftzJRRi8dwKKOSNGOnL+M3JJKeGtdcXO3IKnjUK3du3dj0qRJCAsLg0qlwoYNG+p9XgiBefPmITQ0FB4eHhg9ejTS0tKkCUukUCzqiMihfJdk7aUbExfs8DulK0l5eTl69+6NxYsXN/j5d955B//5z3+wdOlSHDhwAF5eXhg7diyqqqrsnJRImfL1VUgvLINKBQyJUU5R5yp1ACJqGzVmC35Ktm5lcltvDr3Kyfjx4zF+/PgGPyeEwKJFi/Dyyy9j8uTJAIDPP/8cwcHB2LBhA6ZMmWLPqESKVNdLF99BB52ncuYas6eOSKF+PVuM4vJqBHhpMFRBwwtKl5mZifz8fIwePdr2nE6nw6BBg7Bv374GX2M0GmEwGOo9iOjqdp25CAC4sbPjHw32RyzqiBTqu6O5AIBbe4XCzYXf6o4iPz8fABAcXH/7meDgYNvn/mzBggXQ6XS2R0RERJvnJHJUZovAnjRrUTeiK4s6IpK5iuoabD5hLQBuT+DQq9LNnTsXer3e9sjOzpY6EpFsHcspQUmFCb7urugT4Sd1nFbFoo5IgbaeLEBFtRkR/h7oG9lO6jjUBCEh1rN5CwoK6j1fUFBg+9yfabVa+Pr61nsQUcN2pv4+9OqqsFEMZd0NEQEANtQOvd7RpwNUKsc/z9CZREdHIyQkBNu3b7c9ZzAYcODAAQwePFjCZETKUDefbrjChl4Brn4lUpyiMiN2155nODmhg8RpqCFlZWVIT0+3fZyZmYmkpCT4+/sjMjISzz33HN544w107twZ0dHReOWVVxAWFobbb79dutBEClBcZsTxnBIAwPAuLOqISOZ+OJYHs0UgPlyHmEBvqeNQAw4fPoyRI0faPp49ezYAYPr06VixYgX+/ve/o7y8HE888QRKSkowbNgwbN68Ge7u7lJFJlKEPWlFEALoHuqLYF/lfT+xqCNSmA21Gw7f3oe9dHI1YsQICCGu+nmVSoV//etf+Ne//mXHVETKZxt6VWAvHcA5dUSKkllUjqTsErioVZjEDYeJiGwsFoHdZ5S5lUkdFnVECrK+doHE0Nj2CPTRSpyGiEg+UvL0KC6vhrfWFf06KnNXABZ1RAphsQisO5IDALirL4deiYj+qG4rk6GxAYrdkF2Zd0XkhA5nXUbO5Up4a10xpkfD+5kRETmrXbah1yCJk7QdFnVEClHXS3drrxB4aFwkTkNEJB8lFdU4ev4yAOUukgBY1BEpQpXJjB+PXwAA3Nk3XOI0RETysietCBYBdAn2Rpifh9Rx2gyLOiIF2HqyAKXGGnTw88DAKH+p4xARyYrStzKpI9ui7rXXXoNKpar36Natm9SxiGSpbuj1zr4doFbzWDAiojoWi3CK+XSAzDcfjouLw7Zt22wfu7rKOi6RJApLq2zHgt3BY8GIiOpJydPjYqkRXhoX9I9S5lYmdWRdJbm6uiIkhKv4iK5l/ZFcmC0CfSP90InHghER1bPtVCEA4MbOgdC6KnsRmWyHXwEgLS0NYWFh6NSpE6ZOnYrz589LHYlIVoQQ+DbROvR6d78IidMQEcnP9lMFAIBR3ZU99ArIuKdu0KBBWLFiBbp27YoLFy5g/vz5uPHGG5GSkgIfH58GX2M0GmE0Gm0fGwwGe8UlksSxHD3SCsvg7qbGxN6hUschIpKVC/pKnMgzQKUCRnZjUSeZ8ePH2/4/Pj4egwYNQseOHfH111/j0UcfbfA1CxYswPz58+0VkUhy3yZmAwDGxYXA191N4jRERPLyy2nr0GtChB/aeyv/6ERZD7/+kZ+fH7p06YL09PSrXjN37lzo9XrbIzs7244JieyrymTGxqQ8ABx6JSJqyPba+XSjugdLnMQ+HKaoKysrw9mzZxEaevUhJq1WC19f33oPIqXaerIAhqoahOncMSQmQOo4RESyUlltxq/p1p0BRrOok9YLL7yAXbt24dy5c/jtt99wxx13wMXFBffff7/U0Yhk4ZvaBRJ39Qvn3nRERH+yN70IxhoLOvh5oEuwc+wMINs5dTk5Obj//vtRXFyMwMBADBs2DPv370dgoLJ3gyZqjNySSuxJs26meXc/HgtGRPRndateR3cPgkrlHL/4yraoW7t2rdQRiGTr28M5EAK4oZM/OgZ4SR2HiEhWLBZhWyThLPPpABkPvxJRwywWgW9qV73eN4ALJIiI/iwlT4/C2lMkBnVynvOwWdQROZjfzhYj53IlfNxdMb4n96YjIvqzulMkbuqi/FMk/ohFHZGD+eqwtZducp8wuLs5T2NFRNRYv58i4TxDrwCLOiKHUlJRjS0n8gEA9/WPlDgNEZH8/PEUiRFdnWtxJYs6Igey4Wguqmss6B7qi54duA8jEdGf1W047CynSPwRizoiByGEwNpDtQsk+oc7zRJ9IqKm2HrSOYdeARZ1RA4jKbsEp/NLoXVV444E7k1HRPRnhioTfjtrPUViXM8QidPYH4s6Igex9qC1l25Cr1DoPN0kTkNEJD87ThfCZBaIDfJGTKBznCLxRyzqiBxAaZUJG4/lAQCmDOQCCSKihtQtJBsX53y9dACLOiKHsPFYHipNZsQEemFAVDup4xARyU6VyYwdp63HJ45lUUdEcvXlwfMAgPsHRnKBBBFRA/akFaHSZEYHPw+n3R2ARR2RzKXk6pGSa4DGRY07+3KBBBFRQzanWIdex8QFO+0vvyzqiGRu9YEsAMDYniHw99JInIaISH5MZgu2n7ZuZeKsQ68AizoiWTNUmfBdknWBxIODuECCiKghBzMvoaTChAAvDQZE+UsdRzIs6ohk7LujuaioNiM2yBsDo523oSIiupa6Va+juwfDRe2cQ68Aizoi2RJCYPUB6wKJqYO4QIKIqCEWi/h9KxMn3HD4j1jUEcnUkfOXcTq/FO5uXCBBRHQ1x3JKUGAwwlvriiGxAVLHkRSLOiKZWr3f2kt3W+8w6Dx4ggQRUUM21/bSjegaCK2ri8RppMWijkiGLpVX44fkCwCAqYM6SpyGiEiehBDYksKh1zos6ohk6OvD2aiusaBnB1/Eh+ukjkNEJEupBaU4V1wBjYsaI7oGSR1HcizqiGTGbBFYtd+6N91DN0RxgQQR0VX8eNw6ojG8ayC8ta4Sp5Eeizoimdl1phA5lyuh83DDpN5hUschCZjNZrzyyiuIjo6Gh4cHYmJi8Prrr0MIIXU0ItkQQtiKuonxoRKnkQeWtUQy8/k+ay/dvf3D4aFx7km/zurtt9/GkiVLsHLlSsTFxeHw4cOYMWMGdDodnnnmGanjEcnCqQulyCgqh8ZVjVHdg6WOIwss6ohkJKu4HLvOXATABRLO7LfffsPkyZMxYcIEAEBUVBS+/PJLHDx4UOJkRPLxw3HraTsjOfRqw+FXIhlZfeA8hACGdwlEVHsvqeOQRIYMGYLt27fjzJkzAIBjx45h7969GD9+fIPXG41GGAyGeg8iJRNC4MfkuqFXTlOpw9KWSCYqq8346lA2AGDaDeylc2YvvfQSDAYDunXrBhcXF5jNZrz55puYOnVqg9cvWLAA8+fPt3NKIumcyDMgq7gC7m5q3NyNq17rsKeOSCY2JOVCX2lChL8HRrKRcmpff/01Vq9ejTVr1uDIkSNYuXIl3nvvPaxcubLB6+fOnQu9Xm97ZGdn2zkxkX19Xzv0enO3IHhx6NWGfxNEMiCEwIpfzwEApg+OcuoDqQl48cUX8dJLL2HKlCkAgF69eiErKwsLFizA9OnTr7heq9VCq9XaOyaRJOqveuXQ6x+xp45IBvZnXEJqQSk83FxwT/8IqeOQxCoqKqBW12+eXVxcYLFYJEpEJB/Hc/TIuVwJDzcXjOSGw/Wwp45IBlb8lgkAuLNvB57zSpg0aRLefPNNREZGIi4uDkePHsX777+PRx55ROpoRJKrWyAxqnsQt336ExZ1RBLLuVyBrScLAADTh0RJG4Zk4cMPP8Qrr7yCp59+GoWFhQgLC8OTTz6JefPmSR2NSFLccPjaWNQRSeyL/VmwCGBobAC6BPtIHYdkwMfHB4sWLcKiRYukjkIkK0ezS5BbUgkvjQvPem0A59QRSaiiugZfHjgPAHh4SLTEaaghLi4c3iGSi7peutE9guHuxu/NP2NRRySh/0vMgaGqBh0DPDGK25jIEs9bJZIHs0Xg+2PWrUwm9OLQa0OaXNRVVlYiNzf3iudPnDjRKoGInIXFIrC8dhuTGUOioOY2Jm3u5MmTeO+995CXZ/3BsH79+uu+RqXi+0IkB/vOFqOw1AidhxuHXq+iSUXdt99+i86dO2PChAmIj4/HgQMHbJ+bNm1aq4cjUrJdZy4io6gcPlpX3M1tTOxi/vz5uOWWWzB//nwcO3YMW7dulToSETXShiRrh9KE+FBoXDnQ2JAm/a288cYbSExMRFJSEpYvX45HH30Ua9asAcAhCqKm+uxX6zYmUwZG8DBqO/Hz80Pv3r2xdOlSLF++HMePH2/yn2EymbB161bs2bMHxcXFbZCSiP6symTG5pR8AMDtfTpInEa+mvSTxGQyITg4GADQr18/7N69G3fccQfS09M5REHUBGcKSrEnrQhqFfDQ4Cip4ziNW265BYB1SPXf//433nvvvSb/GXfeeSdCQ0Oxbt06tGvXDhUVFejVqxc2b97c2nGJqNa2UwUoM9agg58H+ndsJ3Uc2WpST11QUFC932z9/f2xdetWnDp1qlm/8RI5q//tyQAAjI0LQYS/p8RpnMfdd98NAJg5cyYuXryIF198scl/xvnz5/HJJ58gPDwcaWlp+Mc//oH4+PjWjkpEf7DhqHUe7OQ+YZx/fA1NKuq++OILBAXVn5yo0Wjw5ZdfYteuXa0ajEipCkurbA3UYzd2kjiNcxo/fjxuvfVWvPbaaygvL2/Sa93d3QFY277q6mrMmjULe/fubYuYRATgcnk1dqYWAgDuSODQ67U0qagLDw9HSEhIveeys7MBAEOHDm29VEQK9vlvWag2W9A30g/9OIwgicmTJ+PAgQMIDg7GkCFDsHTp0kafq/rMM8/g0qVLuOuuuzBz5kwsW7YMRUVFbZyYyHn9mHwBNRaBHqG+6MwN2q+pxctHunXrhnnz5qGioqI18hApWkV1DVYdyAIAPM5eOkm5uLhgwoQJeP755/Hyyy+jR48e+P777696/dmzZwEAU6dOhb+/P+bMmYObbroJp0+fxrfffmuv2EROZ8NR66pX9tJdX4uLuq1bt2LLli3o3LkzVqxY0QqRiJTr/xJzUFJhQqS/J8bEhVz/BdQmxo0bh8jISDzwwAM4fvw4PvzwQ6xevRobNmzAc8891+BrZs6ciejoaAwePBhPPvkkFi9ejJiYGLz88sucU0fURrIvVeBw1mWoVMCk3mFSx5G9Fu+jMGTIEBw4cACff/45/vnPf+LDDz/EokWLcOONN7ZGPiLFMFsElu21bmPy6LBouHCyr92dPXsWMTExWLhwIXr16nXFEWDLli1Dt27dGnxt3Z52b731Fg4dOoTc3Fxs3LgR27dvR1RUFNLT09s8P5Gz+a52b7ohMQEI0blLnEb+Wm1zrIceegh33303Fi5ciPHjx2PcuHF49913ER3N8yyJAODnE/k4V1wBnYcb7ukfLnUcpzRz5kykp6cjJCQE8fHx9R46nQ4A8NNPP13zz/j666+RlJRk+/jnn3/G6tWr2zI2kVMSQmB97dDrZO5N1yitviXzmDFj8Nhjj2H9+vXo0aMH/v73v6OsrKy1vwyRQxFCYOku65yshwZ3hKeGmw1LYevWrcjMzMSkSZNQWFiI3NxcvPHGG/D390dsbCwAoFOna891dHd3x8mTJ20fjxkzBikpKW2am8gZncgz4OzFcmhd1RjXk9NVGqPFP1mWLl2KQ4cO4dChQzh16hTUajV69uyJmTNnonfv3li7di169OiBdevWoX///q2RmcjhHMi8hGM5emhc1Zg+JErqOE6vJb1ty5Ytw3333YcRI0agT58+SE5O5ubrRG3g28QcAMDoHsHwdXeTOI1jaHFP3Ztvvgm9Xo+HHnoIO3bsQElJCRITE7F48WI88cQT+OWXXzBz5kw8/PDDLfo6CxcuhEqluuokZiI5+7i2l+6efuFo762VOA21pLctLi4OiYmJuPHGG3Hu3Dl07NgRmzZtaquoRE7JWGO2nfV6Tz9OV2msFvfU1e1Tdy2PPvooXnnllWZ/jUOHDuHjjz/mCjNySKfzDdiRehFqFbcxkYuW9rZpNBrce++9uPfee9swJZHz2n6qECUVJoT4uuPGzoFSx3EYrT6nriFBQUH45ZdfmvXasrIyTJ06FZ9++inateNGreR4PtltPRJsfM9QRLX3kjgNAU3rbRNC2DkdEX1z2NphdGffDtwpoAnsMltbpVJh+PDhzXrtrFmzMGHCBIwePRpvvPFGKycjals5lyuwMcl6JNgTN7GXTk4a29vW2JMmiKh1FBiqsOvMRQDA3Rx6bRJZL8Fbu3Ytjhw5gkOHDjXqeqPRCKPRaPvYYDC0VTSiRvl0dwZqLAJDYwPQO8JP6jhERLK37kguLALo37EdOgV6Sx3Hodhl+LU5srOz8eyzz2L16tW2A7SvZ8GCBdDpdLZHREREG6ckurqiMiPWHrIOITw9IlbiNERE8ieEsA29cj/PppNtUZeYmIjCwkL07dsXrq6ucHV1xa5du/Cf//wHrq6uMJvNV7xm7ty50Ov1tkdjFnEQtZXlv2bCWGNB73AdhsQESB2HiEj2jpy/jIyicni4uWBCPI8FayrZDr+OGjUKycnJ9Z6bMWMGunXrhjlz5lxxvA8AaLVaaLXcLoKkV1plwuf7sgAAT42I5T5mRESN8M1h695043uFwFsr2xJFtmT7N+bj44OePXvWe87LywsBAQFXPE8kN6v2n0dpVQ1ig7wxpkew1HGIiGSvoroGPxy/AAC4px+nTzWHbIdfiRxVZbUZ/9tj3cZk5vAYqLkcn4joujan5KPMWIMIfw8MivaXOo5Dkm1PXUN27twpdQSi61pz8DyKy6sR4e+ByX04J4SIqDHqhl7v7hvBX4abiT11RK2oymTGJ7utR4I9NTwWbi78FiMiup5zReXYl1EMlQq4q18HqeM4LP7EIWpF3ybmoMBgRKjOnQ0TEVEjfXnwPABgeJdAhLfzlDiN42JRR9RKTGYLluy09tI9eVMnaF2vXKFNRET1GWvM+CbROvT6wMBIidM4NhZ1RK1k/ZFc5JZUor23FlPYMBERNcqWEwW4VF6NEF933NwtSOo4Do1FHVErMJkt+HBHGgDgiZui4e7GXjoiosZYc8C6p+e9AyLgynnILcK/PaJWsP5oLrIvVSLAS4MHb+godRwiIodw9mIZ9mdcgloFTBnAvelaikUdUQvVmC1YvCMdAPDk8E7w1DjUTkFERJJZW7tAYmTXIIT5eUicxvGxqCNqofVHc5FVXMFeOiKiJqgymfFt3QKJQZyH3BpY1BG1wB976Z64ib10RESNteVEPi5XmBCmc8eIrlwg0RpY1BG1wLqjuThXXAF/9tIRETXJ6gPWodf7BkTChSdItAoWdUTNVF1jwX+2W1e8zhzeCV5a9tIRETVGemEpDmZegotahfu4QKLVsKgjaqZvE3OQc9m6L920G6KkjkNE5DBW7bf20t3cLQghOneJ0ygHizqiZjDWmPHRL9ZeuqdHxMBDw33piIgao8xYY1sg8dBgTltpTSzqiJph7cFs5OmrEOLrzlVb1CZyc3Px4IMPIiAgAB4eHujVqxcOHz4sdSyiFvu/xByUGWsQE+iFYbHtpY6jKJwERNREldVm24rXWTfH8vQIanWXL1/G0KFDMXLkSGzatAmBgYFIS0tDu3btpI5G1CIWi8DKfecAANOHREGl4gKJ1sSijqiJPt93DoWlRnTw88C9/cOljkMK9PbbbyMiIgLLly+3PRcdHS1hIqLWsTe9CBkXy+GtdcWdfdl+tjYOvxI1gaHKhCW7zgIAnr+lC7Su7KWj1rdx40b0798f99xzD4KCgpCQkIBPP/30qtcbjUYYDIZ6DyI5WvnbOQDA3f3C4c0dA1odizqiJvjfnkyUVJgQG+SNOxI6SB2HFCojIwNLlixB586dsWXLFjz11FN45plnsHLlygavX7BgAXQ6ne0REcEtIkh+sorL8UtqIQAukGgrLOqIGqm4zIhlezIAAH+7pQs3y6Q2Y7FY0LdvX7z11ltISEjAE088gccffxxLly5t8Pq5c+dCr9fbHtnZ2XZOTHR9n+/LghDA8C6B6BToLXUcRWJRR9RIi3ecRXm1Gb066DCuZ4jUcUjBQkND0aNHj3rPde/eHefPn2/weq1WC19f33oPIjkpN9bg68PWXzYeHhIlbRgFY1FH1Ag5lyuwan8WAODFsV25Yova1NChQ5GamlrvuTNnzqBjRw5ZkWNafzQXpVU1iArwxPAugVLHUSwWdUSN8P7WM6g2WzAkJgA3dua+StS2nn/+eezfvx9vvfUW0tPTsWbNGnzyySeYNWuW1NGImkwIYVsgMW1wFNScutJmWNQRXcepCwasP5oLAJgzrht76ajNDRgwAOvXr8eXX36Jnj174vXXX8eiRYswdepUqaMRNdnutCKkFZbBU+OCe7gNVJviemKi63h3SyqEACbEh6J3hJ/UcchJTJw4ERMnTpQ6BlGLfbrbusBsyoBI+Lq7SZxG2dhTR3QNBzKK8cvpQriqVXhhTFep4xAROZQTeXrsTS+Ci1qFGUOjpI6jeCzqiK7CYhF466dTAID7BkQgur2XxImIiBzL//ZkAgBu7RWKCH9PidMoH4s6oqv4IfkCjuXo4aVxwXOju0gdh4jIoeSVVOL7Y3kAgMdv5DF39sCijqgBxhoz3tl8GgDw5PAYBPpoJU5ERORYVvx2DjUWgRs6+SM+3E/qOE6BRR1RA77Yl4Wcy5UI8tHiMf6GSUTUJKVVJnx5wLpZ9hM3dZI4jfNgUUf0JyUV1fjwl3QAwAtjusJTw0XiRERN8dWhbJQaaxAb5I0RXYKkjuM0WNQR/cl/tqdDX2lC12Af3NWPeyoRETWFyWzBZ3utCyQevzGamw3bEYs6oj/IuFiGz/edAwC8PLE7XNgYERE1yU/JF5Cnr0J7by0m9+kgdRynwqKO6A8WbDqNGovAzd2CcGNnnk9IRNQUFovAkp1nAQDTB3eEu5uLxImcC4s6olq/nS3C1pMFcFGr8I9bu0kdh4jI4Ww/XYjT+aXw1rriocFRUsdxOizqiACYLQKv/2DdaPjBQZGIDfKROBERkWMRQuCjHdZFZg/e0BE6Tx4JZm8s6ohgXal16oIBvu6ueJYbDRMRNdmv6cU4ll0Crasajw7jVlBSYFFHTk9fYcJ7P6cCAJ6/pQv8vTQSJyIicjwf7UgDANw/MJIbtkuERR05vQ+2p+FSeTU6B3njwRs6Sh2HiMjhJGZdwv6MS3BzUXGzYQmxqCOnll5YatvCZN6kHnBz4bcEEVFTfVS7YfudCeEI8/OQOI3z4k8wclpCCMz//iRqLAKjuwdzCxMiomZIydVjR+pFqFXAUyNipI7j1FjUkdPacqIAe9KKoHFR4+UJ3aWOQ0TkkP6709pLNzE+DFHtvSRO49xY1JFTqqw24/UfTgIAnhzeiQ0REVEzpBeWYlNKPgBg1shYidMQizpySkt2piO3pBId/Dzw9Ag2REREzfHvbWkQAhjTIxhdQ7i/p9RY1JHTySoux9LdGQCAlyd0h4eGx9gQETXV6XwDfjx+AYB1OyiSHos6cipCCLy28QSqaywYFtse43qGSB2JiMgh/XvrGQDAhF6h6B7qK3EaAljUkZPZcqIAO1Ivws1FhfmT46BSqaSORETkcFJy9dhyogAqFfDc6M5Sx6FaLOrIaZQba/Cv708AAJ68KQYxgd4SJyIickx1vXSTe4ehczDn0smFbIu6JUuWID4+Hr6+vvD19cXgwYOxadMmqWORA/vPL2nI01chvJ0HV2kRETXT0fOXsf10IVzUKp6VLTOyLerCw8OxcOFCJCYm4vDhw7j55psxefJknDhxQupo5IBS80uxbE8mAGD+bXFcHEFE1Ezv1/bS3ZHQAdHcDkpWXKUOcDWTJk2q9/Gbb76JJUuWYP/+/YiLi5MoFTkii0XgH+uTUWMRuKVHMEZ1D5Y6EhGRQzp07hL2pBXBVa3Cs6M4l05uZFvU/ZHZbMY333yD8vJyDB48+KrXGY1GGI1G28cGg8Ee8Ujmvjx0HolZl+GlccH82/gLARFRcwgh8N6WVADAPf0jEOHvKXEi+jPZDr8CQHJyMry9vaHVajFz5kysX78ePXr0uOr1CxYsgE6nsz0iIiLsmJbkqLC0Cgs3nQYA/G1MVx40TUTUTDtTL+JA5iVoXNX4682clyxHsi7qunbtiqSkJBw4cABPPfUUpk+fjpMnT171+rlz50Kv19se2dnZdkxLcvT6D6dQWlWDXh10mD4kSuo4REQOyWwRtl+QZwyJ4i/IMiXr4VeNRoPYWOtvA/369cOhQ4fwwQcf4OOPP27weq1WC61Wa8+IJGO/nC7A98fyoFYBb93RCy5q7klHRNQc647kILWgFL7urnhqRIzUcegqZN1T92cWi6XenDmiqykz1uDl9SkAgEeHRaNXuE7iREREjqnKZLateJ01MhZ+nhqJE9HVyLanbu7cuRg/fjwiIyNRWlqKNWvWYOfOndiyZYvU0cgBvLv5NPL0VYj098TsW7pKHYeIyGGt/O0cLuirEKZz5zQWmZNtUVdYWIiHHnoIFy5cgE6nQ3x8PLZs2YJbbrlF6mgkc4lZl/H5/iwA1mFX7klHRNQ8JRXVWLwjHQAwe0xXuLuxPZUz2RZ1y5YtkzoCOaAqkxlz/u84hADu6ReOYZ3bSx2JiMhh/XfnWRiqatAtxAd3JHSQOg5dh0PNqSO6nv9sT0N6YRkCfbT454TuUschInJYuSWVWPHbOQDAnHHduNjMAbCoI8VIztHj490ZAIA3bu/JybxERC3w9qbTqK6x4IZO/hjRNVDqONQILOpIEaprLHjx22MwWwQmxodibFyI1JGIiBzW4XOXsPFYHlQq4OUJPaBSsZfOEbCoI0X4aEc6TueXwt9Lw6PAiIhawGIRmP+9daP/KQMi0LMDt4RyFCzqyOGl5Optq7P+NTkOAd7cgJqIqLm+PZKD5Fw9fLSu+NsYbgnlSFjUkUMz1pgx++skmC0CE+JDMTE+TOpIRK1u4cKFUKlUeO6556SOQgpXWmXCO5tTAQDPjOqM9vwl2aGwqCOH9sG2NJwpKEN7bw1en9xT6jhEre7QoUP4+OOPER8fL3UUcgKLd5xFUZkR0e29uNGwA2JRRw7ryPnLWLrrLADgjdt7wd+Lq11JWcrKyjB16lR8+umnaNeundRxSOGyisvx2d5MAMDLE7pD48oSwdHwHSOHVFFdg9lfJcEigNv7hGFcT652JeWZNWsWJkyYgNGjR1/zOqPRCIPBUO9B1FRv/ngK1WYLbuzcHjd3C5I6DjWDbE+UILqWN388hXPFFQjVuWM+h11JgdauXYsjR47g0KFD1712wYIFmD9/vh1SkVLtOF2In08WwEWtwryJ3MLEUbGnjhzOjtRCrD5wHgDw3j29ofNwkzgRUevKzs7Gs88+i9WrV8Pd3f2618+dOxd6vd72yM7OtkNKUorKajPmbUwBADw6LBqdg30kTkTNxZ46ciiXyqvx92+PAwBmDI3C0Fie7UrKk5iYiMLCQvTt29f2nNlsxu7du/HRRx/BaDTCxeX3g9W1Wi20Wq5SpOZZvCMd2ZcqEapzx7OjOksdh1qARR05DCEE5vzfcVwsNSI2yBtzxnWTOhJRmxg1ahSSk5PrPTdjxgx069YNc+bMqVfQEbXE2Ytl+Hi3dcHZq5N6wEvLssCR8d0jh/HlwWxsPVkANxcVPpjSB+5u/MFGyuTj44OePevPFfXy8kJAQMAVzxM1lxACr2xIgcksMLJrII9XVADOqSOHcPZiGV7/wXpszd/HdkNcGI+tISJqiY3H8vDb2WJoXdWYf1tPLo5QAPbUkewZa8x4du1RVJrMGBobgEeHRUsdicjudu7cKXUEUhB9pQmv/3AKAPDXm2MRGeApcSJqDeypI9l7Z3MqUnINaOfphv93Tx+o1fxtkoioJd7ZfBpFZUZ0CvTC4zd1kjoOtRIWdSRrO04XYlntDufv3t0bIbrrb+9ARERXt+9ssW1bqDdu7wmtK+cnKwWLOpKtQkMVXvjmGADg4SFRGN0jWOJERESOrbLajJfWWbeFun9gJIbEcFsoJWFRR7Jktgg8uzYJxeXV6B7qi5fGc/sSIqKW+ve2M8gqrkCIrzvm3sp2VWlY1JEsffhLGvZlFMNT44KPHkjg9iVERC10LLsE/9uTAQB4686e8HXnaTxKw6KOZOe39CJ8sD0NAPDmHT0RE+gtcSIiIsdWXWPB3789DosAJvcJw83dOJ1FiVjUkawUllbh2a+SIARwb/9w3JEQLnUkIiKHt3hHOlILShHgpcGrk+KkjkNthEUdyUaN2YJnvjyKi6VGdAn2xvzbuHM+EVFLncjT47870wEAr90WB38vjcSJqK2wqCPZ+Pe2M9ifcQmeGhf8d2o/eGg4j46IqCWqTGY8/1USTGaBsXHBmBgfKnUkakMs6kgWfjldgMU7rIdKv31XPGKDOI+OiKil3tuSijMFZWjvrcFbd/TiUWAKx6KOJHe+uALPf/X7fnSTeodJnIiIyPH9drYIy361bt7+9l3xCPDWSpyI2hqLOpJUZbUZT65KhL7ShD4RfvjHrd2ljkRE5PAMVSa88PUxCAHcPzACo7pztaszYFFHkhFC4J8bknHqggHtvTVY8mBfaFz5T5KIqKVe23gCefoqRPp74uUJPaSOQ3bCn6AkmVX7s7DuSC5c1Cp8eH9fhOo8pI5EROTwNiVfwLojuVCrgPfv7Q0vravUkchOWNSRJA5kFGP+9ycBAHPGdcXgmACJExEROb6cyxWY83/Ws11nDo9B/yh/iRORPbGoI7vLK6nE06uPoMYicFvvMDx+YyepIxEROTxT7V6fhqoa9I7ww3Oju0gdieyMRR3ZVZXJjCe+OIzi8mr0CPXF23fFc4k9EVEr+PfWMzhyvgQ+Wld8OCWBc5SdEN9xshshBOb833Gk5Brg76XBJw9xg2EiotawJ+0iluyy7vW58K54RAZ4SpyIpMCijuzmvzvP4rukPLiqVVj8QF+Et2OjQ0TUUoWlVXi+9szsBwZFYgJPjXBaLOrILraeLMB7P6cCsJ49yIURREQtZ7EIzP7qGIrKqtEtxAfzJnL7EmfGoo7a3KkLBjy39iiEAKbd0BEP3tBR6khERIrwwfY07E0vgoebCz56IAHubpzS4sxY1FGbKiytwqMrDqG82owhMQGYN4m/RRIRtYZfThfgg+1pAIA3bu+J2CAfiROR1FjUUZupMpnx+OeJyNNXoVOgF5ZM7Qc3F/6TIyJqqXNF5XhubRIA6wjIXf3CpQ1EssCfsNQmLBaBv319DMeyS+Dn6YbPpg+AztNN6lhERA6voroGM1clwlBVg76RfniF8+ioFos6ahNvbzmNH5MvwM1FhY8f7Ieo9l5SRyIicnhCCMxdl4zT+aVo763Bf6f24350ZMN/CdTqVu3Pwse7MgAA79wdj0GduNKViKg1rPjtHL5LyoOLWoWPHuiLEJ271JFIRljUUav65XQB5n2XAgCYfUsX3JHAeR5ERK1hT9pFvPHjKQDA3PHdcAN/YaY/YVFHreZYdglmrT4KiwDu7heOv94cK3UkIiJFSC8sw9Orj8BsEbgzoQMeHRYtdSSSIRZ11CrOFZXjkRWHUGky48bO7fHWHb14pisRUSu4XF6NR1ceQmlVDfp3bIcFd7F9pYaxqKMWKyozYvrygygur0ZcmC+WPMiJu0REraG6xoKnViciq7gCHfw8sHRaP2hducEwNYw/ealFSqtMeHj5QWQVVyC8nQeWzxgAb62r1LGIiByeEAKvbkzB/oxL8NK4YNnD/dHeWyt1LJIx2RZ1CxYswIABA+Dj44OgoCDcfvvtSE1NlToW/UGVyYwnv0hESq4B/l4afP7IQAT5cCUWEVFr+HRPBr48mA2VCvjP/QnoFuIrdSSSOdkWdbt27cKsWbOwf/9+bN26FSaTCWPGjEF5ebnU0QiA2SIw++sk/Ha2GF4aF6ycMRCdAr2ljkVEpAjfJeXirZ9OAwD+eWt3jOoeLHEicgSyHSfbvHlzvY9XrFiBoKAgJCYm4qabbpIoFQHWIYF/rEvGT8n50Lio8elD/dErXCd1LCIiRfg1vQgvfHMMAPDI0GiudKVGk21R92d6vR4A4O/vf9VrjEYjjEaj7WODwdDmuZyNEAJv/ngKXx3OhloFLJrSB0Ni20sdi4hIEU7k6fHkF4kwmQUmxIfi5QndudKVGk22w69/ZLFY8Nxzz2Ho0KHo2bPnVa9bsGABdDqd7REREWHHlM7ho1/S8b+9mQCAhXfG49ZeoRInIiJShuxLFXh4+SGUGWtwQyd/vH9vb6jVLOio8RyiqJs1axZSUlKwdu3aa143d+5c6PV62yM7O9tOCZ3D//Zk4P9tPQMAeGViD9w7gEUzEVFrqNsa6mKpEd1CfPDxtP7cuoSaTPbDr3/5y1/www8/YPfu3QgPv/aRU1qtFlotl3u3hVX7s2zH0zw3ujPneBARtRJ9hQnTlh1ExsVyhOncsXzGAOg83KSORQ5ItkWdEAJ//etfsX79euzcuRPR0SwipPJtYg5e3mA9z3Xm8Bg8O6qzxImIiJShzFiDh1ccxKkLBrT31mLVY4MQqvOQOhY5KNkOv86aNQurVq3CmjVr4OPjg/z8fOTn56OyslLqaE5lw9FcvPitdRXWw0OiMGdcV07aJWpj3KfTOVSZzHhs5SEcPV8CnYcbVj3GraGoZWRb1C1ZsgR6vR4jRoxAaGio7fHVV19JHc1pbDyWh9lfJ0EI4P6BkZg3sQcLOiI74D6dylddY8FTqxKxP+MSvLWu+PyRgdxcmFpM1sOvJJ0fjufh+a+SYBHAlAERePP2nlyFRWQn3KdT2aprLPjLmiPYkXoR7m5qLJveH70j/KSORQog26KOpLPxmLWgM1sE7ukXjrfu6MWCjkhC19unk3t0Og5jjRmzVh/BtlOF0Liq8fG0/hjUKUDqWKQQsh1+JWlsOJqL59YehdkicG//cLx9VzwLOiIJNWafTu7R6RiqTGbM/CIR204VQutqPY1neJdAqWORgrCoI5tvDmdj9tfWIdf7B0Zg4Z0s6Iik1ph9OrlHp/xVmcx48ovEPwy5DmBBR62Ow68EAPhifxZeqd22ZOqgSLw+mXPoiKTW2H06uUenvFVU1+DJLxKxJ60IHm4uWPZwfwyJ4fGK1PpY1BH+tyfDtrHwjKFRXOVKJDHu06kc+goTHll5CIlZl+GpccFnDw/ADZxDR22ERZ0TE0Jg0bY0fLA9DQDw9IgYvDiW+9ARSW3WrFlYs2YNvvvuO9s+nQCg0+ng4cGNaR1FoaEK05YdRGpBKXzdXbF8xgD069jwYhei1sCizkkJIfD6D6fw2a+ZAIC/3dIFf7k5lgUdkQwsWbIEADBixIh6zy9fvhwPP/yw/QNRk2UVl+PBZQeQfakSQT5afP4o96GjtseizgnVmC2Yuy4Z3yTmAABem9QDDw/l8A6RXHCfTsd2Ms+A6csP4mKpER0DPPHFI4MQGeApdSxyAizqnEyVyYy/rDmKbacKoFYB79zdG3f3u/oEbCIiarzdZy7i6dVHUGasQfdQX6x8ZACCfNyljkVOgkWdEzFUmfDYysM4mHkJGlc1Fj/QF7f0CJY6FhGRInx58Dxe3pACs0VgULQ/PnmoP3QeblLHIifCos5J5Our8PDygzidXwofrSv+N527mBMRtQaLReDdn1OxZOdZAMCdCR2w8K54aFy5FSzZF4s6J5BeWIqHlh1Enr4KgT5arJgxAHFhOqljERE5vCqTGS98cww/HL8AAHh2VGc8N7ozF52RJFjUKdzBzEt4/PPD0Fea0CnQCytnDESEPyfsEhG1VF5JJZ78IhHJuXq4uaiw4M54zlEmSbGoU7CNx/LwwtfHUG22ICHSD8umD4C/l0bqWEREDu/QuUt4alUiisqq0c7TDYun9uUpESQ5FnUKJITAkl1n8c7mVADA2LhgfDAlAe5uLhInIyJyfGsOnMerG1NgMgt0C/HBpw/15wgIyQKLOoWprrHg5Q3J+PqwdQ+6R4ZG458TusOF57gSEbVIlcmMf/1wEmsOnAcATOgVinfviYenhj9KSR74L1FBSiqqMXNVIvZnXIJaBcybyE2FiYhaQ1ZxOZ5efQQn8gxQqYAXxnTF0yNiuCCCZIVFnUKkF5bh8c8PI7OoHF4aF3z0QF+M7BYkdSwiIoe3OeUCXvzmOEqNNWjn6YZ/39cHI7qyfSX5YVGnALvOXMRf1hxBaVUNOvh54H/T+6N7KM8YJCJqieoaCxZsOoXlv54DAPTr2A4f3p+AMD8PaYMRXQWLOgcmhMCyvZl466dTsAigf8d2WDqtH9p7a6WORkTk0NILy/DcV0eRkmsAADxxUye8OLYr3Fy4oTDJF4s6B1VlMuMf65Kx7mguAODufuF4846e0LpyhSsRUXMJIbDqwHm8+eNJVJks8PN0wzt3xWNMXIjU0Yiui0WdA8orqcTMVYk4nqOHi1qFf97aHTOGRnHCLhFRCxSVGTHn2+PYfroQAHBj5/Z4757eCPZ1lzgZUeOwqHMwv6UX4S9fHsWl8mr4ebrhvw/0xZBYbnhJRNQSm5Iv4JXvUlBUVg2NixpzxnfDjCFRUHM7KHIgLOochBACn+zOwNubT8MigLgwXyx9sB83vCQiaoGiMiPmfZeCn5LzAQBdg32waEofLjYjh8SizgHoK0148Ztj+PlkAQDgrr7W+XM8IYKIqHmEENh4LA+vbTyByxUmuKhVeHpEDP5ycyznJpPDYlEncyfy9Hh69RFkFVdA46LGK5N64MFBkZw/R0TUTDmXK/Dqdydsc+e6h/ri3bvj0bODTuJkRC3Dok6mhBBYc/A85n9/EtU1FnTw88B/p/ZF7wg/qaMRETkkk9mCZXsz8cG2NFSazHBzUWHWyFg8PSIWGlduVUKOj0WdDJVWmfCP9Sn4/lgeAODmbkF4/97e8PPUSJyMiMgxHTp3Cf9cn4wzBWUAgIHR/njz9p7oHOwjcTKi1sOiTmaOZZfgmbVHkVVcAVe1Cn8f1xWPDevEFVhERM1wQV+JhZtO47sk6y/J/l4a/OPW7rirbwdOYyHFYVEnExaLwP/2ZuCdzamosQh08PPAf+5PQL+O7aSORkTkcCqrzfhkdwaW7EpHlckClQq4r38E5ozrhnZeHPUgZWJRJwMFhir87etj2JteBAC4tVcIFtwRD52nm8TJiIgci8Ui8P3xPLy96TTy9FUArEcovjopDr3CuRCClI1FncS2nMjHnP87jpIKE9zd1Jg3MQ73D4zgsAARURMIIbA7rQhvbzqNkxes57WG6dwx99bumBgfyjaVnAKLOomUVpnw+g8n8fXhHABAzw6+WHRfAmKDvCVORkTkWJKyS/D2ptPYl1EMAPDWuuLJmzrhsRs7wUPDPefIebCok8DBzEuY/XUSci5XQqUCnripE/52S1cuqSciaoLkHD0+2H4G205Z95vTuKgxbXBHzBoZC3/OmyMnxKLOjqpMZry7JRWf/ZoJIYDwdh54/94+GBjtL3U0IiKHcSy7BB9sT8MvtZsHq1XAHQnheP6Wzghvx6MTyXmxqLOTI+cv44VvjiHjYjkA4N7+4XhlYg/4uHMxBBFRYyRmXcZHv6RhR+pFANZi7vY+HTDr5ljEBHLqChGLujZWWW3G//s5Fctqe+eCfLRYeFcv3NwtWOpoRESyJ4TAztSLWLLrLA5mXgJQW8wldMBfRsaiE4s5IhsWdW1of0Yx5q5LRmaRtXfuzoQOmDepB0+GICK6jiqTGd8l5eJ/ezKRVmg9BcLNRYXb+3TA0yNjEd3eS+KERPLDoq4N6CtNWLjpNL48eB4AEOyrxYI72TtHRHQ9+foqrNqfhTUHz+NSeTUA62rW+wdG4JFh0QjVeUickEi+WNS1IiEENqXk47WNJ1BYagQA3D8wEi+N7wadB+fOERE1RAiBfRnFWLU/Cz+fKECNRQCw7jM3Y2g07hsYAV/OPya6LhZ1rSS3pBLzNqRge+1qrOj2XlhwZy/c0ClA4mRERPJUXGbEuiO5WHvoPM7WLiIDgIHR/pgxJAq39AiGqwu3eiJqLBZ1LWQyW7BsbyY+2JaGSpMZbi4qPDU8Bk+PjIW7Gze9JCL6oxqzBXvSivBtYg5+PpkPk9naK+elccHtCR3w4A0d0T3UV+KURI6JRV0L7M8oxisbUmyTeAdG+ePNO3qic7CPxMmIiOTlZJ4B64/m4LukPNv0FACID9dhyoBI3NYnDN5a/kgiagl+BzVDvr4Kb/10ChuP5QEAArw0mHtrd9zVtwPPFyQiqpVVXI4fjl/AxqQ8pBaU2p7399Jgcp8w3NMvAj3C2CtH1FpY1DWBscaMZXszsfiXdJRXm6FSAVMHReKFMV25TQkREYDzxRX4MfkCNqVcwPEcve15jYsao7oH4faEDhjZNYjHIhK1ARZ1jSCEwLZThXjjx5PIKq4AAPSN9MO/JvdEzw46idMREUlHCIFTF0rx88l8bE7Jx+n833vk1CpgSEx73NY7DGN7hnAXAKI2xqLuOk7mGfDmTyfxa3oxAOuJEC+N74bb+3SAWs2hViJyPlUmM/ZlFGPH6UJsO1mAPH2V7XNqFXBDpwDc2isUY+NCEOijlTApkXORdVG3e/duvPvuu0hMTMSFCxewfv163H777Xb52oWGKry/9Qy+OpwNIQCNqxqPDYvGrJGx8OJkXiKyg8WLF+Pdd99Ffn4+evfujQ8//BADBw60ew4hBDKKyrHnzEXsOnMRv50thrHGYvu8u5saw2LbY0xcCEZ3D4a/F6ejEElB1tVJeXk5evfujUceeQR33nmnfb6msQaf7M7AJ7szUGkyAwAmxodizrhuiPD3tEsGIqKvvvoKs2fPxtKlSzFo0CAsWrQIY8eORWpqKoKCgtr86xcYqrDvbDF+TS/Cb2eLkVtSWe/zoTp3jOgahNHdgzAkpj08NNzCiUhqsi7qxo8fj/Hjx9vla1XXWLD20Hn8Z3saisqsR9P0jfTDPyd0R7+O/nbJQERU5/3338fjjz+OGTNmAACWLl2KH3/8EZ999hleeumlVv96uSWVOJhZjIOZl3EgoxgZReX1Pq9xUaN/VDvc1CUQI7sGoUuwN1f7E8mMrIu6pjIajTAaf9//yGAwNOp1204W4PU/LIKICvDEnHHdMK5nCBstIrK76upqJCYmYu7cubbn1Go1Ro8ejX379l1xfXPbPgCYuy4Zu1IL682LAwCVCugZpsOQmAAMiW2PgVH+7I0jkjlFFXULFizA/Pnzm/y6fEMVsoor0N5bi2dHd8aUARFw49E0RCSRoqIimM1mBAcH13s+ODgYp0+fvuL65rZ9gLWHLk9fBRe1Cj3DfDEgyh83dArAgCh/6Dy5WpXIkSiqqJs7dy5mz55t+9hgMCAiIuK6r7tvQASMNRZMGRDBRRBE5HCa2/YBwF9GxmLmTZ3QJ9IPnhq2f0SOTFHfwVqtFlpt05fPu7mo8eiw6DZIRETUdO3bt4eLiwsKCgrqPV9QUICQkJArrm9u2wcAA6M5Z5hIKTjGSEQkMxqNBv369cP27dttz1ksFmzfvh2DBw+WMBkRyZmse+rKysqQnp5u+zgzMxNJSUnw9/dHZGSkhMmIiNrW7NmzMX36dPTv3x8DBw7EokWLUF5eblsNS0T0Z7Iu6g4fPoyRI0faPq6bMzJ9+nSsWLFColRERG3vvvvuw8WLFzFv3jzk5+ejT58+2Lx58xWLJ4iI6qiEEELqEG3FYDBAp9NBr9fD19dX6jhE1Ej83m0Z/v0ROa6WfP9yTh0RERGRArCoIyIiIlIAFnVERERECsCijoiIiEgBWNQRERERKQCLOiIiIiIFYFFHREREpAAs6oiIiIgUgEUdERERkQKwqCMiIiJSABZ1RERERArgKnWAtlR3rK3BYJA4CRE1Rd33rIKPpm5TbPuIHFdL2j9FF3WlpaUAgIiICImTEFFzlJaWQqfTSR3D4bDtI3J8zWn/VELBvwpbLBbk5eXBx8cHKpXqmtcaDAZEREQgOzsbvr6+dkrY9nhfjkep99aU+xJCoLS0FGFhYVCrOUukqdj2WSn13pR6X4By781e7Z+ie+rUajXCw8Ob9BpfX19F/UOqw/tyPEq9t8beF3vomo9tX31KvTel3heg3Htr6/aPvwITERERKQCLOiIiIiIFYFFXS6vV4tVXX4VWq5U6SqvifTkepd6bUu/L0Sn5fVHqvSn1vgDl3pu97kvRCyWIiIiInAV76oiIiIgUgEUdERERkQKwqCMiIiJSAMUWdYsXL0ZUVBTc3d0xaNAgHDx48JrXf/PNN+jWrRvc3d3Rq1cv/PTTT/U+L4TAvHnzEBoaCg8PD4wePRppaWlteQtX1ZR7+/TTT3HjjTeiXbt2aNeuHUaPHn3F9Q8//DBUKlW9x7hx49r6Nq7QlPtasWLFFZnd3d3rXSOX96wp9zVixIgr7kulUmHChAm2a+Tyfu3evRuTJk1CWFgYVCoVNmzYcN3X7Ny5E3379oVWq0VsbCxWrFhxxTVN/d6lKym1/VNq2wew/QMcp/2TddsnFGjt2rVCo9GIzz77TJw4cUI8/vjjws/PTxQUFDR4/a+//ipcXFzEO++8I06ePClefvll4ebmJpKTk23XLFy4UOh0OrFhwwZx7Ngxcdttt4no6GhRWVlpr9sSQjT93h544AGxePFicfToUXHq1Cnx8MMPC51OJ3JycmzXTJ8+XYwbN05cuHDB9rh06ZK9bkkI0fT7Wr58ufD19a2XOT8/v941cnjPmnpfxcXF9e4pJSVFuLi4iOXLl9uukcP7JYQQP/30k/jnP/8p1q1bJwCI9evXX/P6jIwM4enpKWbPni1OnjwpPvzwQ+Hi4iI2b95su6apf190JaW2f0pt+4Rg+1fHUdo/Obd9iizqBg4cKGbNmmX72Gw2i7CwMLFgwYIGr7/33nvFhAkT6j03aNAg8eSTTwohhLBYLCIkJES8++67ts+XlJQIrVYrvvzyyza4g6tr6r39WU1NjfDx8RErV660PTd9+nQxefLk1o7aJE29r+XLlwudTnfVP08u71lL369///vfwsfHR5SVldmek8P79WeNadj+/ve/i7i4uHrP3XfffWLs2LG2j1v690XKbf+U2vYJwfbvahyh/ZNb26e44dfq6mokJiZi9OjRtufUajVGjx6Nffv2Nfiaffv21bseAMaOHWu7PjMzE/n5+fWu0el0GDRo0FX/zLbQnHv7s4qKCphMJvj7+9d7fufOnQgKCkLXrl3x1FNPobi4uFWzX0tz76usrAwdO3ZEREQEJk+ejBMnTtg+J4f3rDXer2XLlmHKlCnw8vKq97yU71dzXe/7rDX+vpydUts/pbZ9ANu/a1FK+2fPtk9xRV1RURHMZjOCg4PrPR8cHIz8/PwGX5Ofn3/N6+v+25Q/sy00597+bM6cOQgLC6v3j2fcuHH4/PPPsX37drz99tvYtWsXxo8fD7PZ3Kr5r6Y599W1a1d89tln+O6777Bq1SpYLBYMGTIEOTk5AOTxnrX0/Tp48CBSUlLw2GOP1Xte6verua72fWYwGFBZWdkq/76dnVLbP6W2fQDbv6tRUvtnz7bPtcVpyWEsXLgQa9euxc6dO+tNqp0yZYrt/3v16oX4+HjExMRg586dGDVqlBRRr2vw4MEYPHiw7eMhQ4age/fu+Pjjj/H6669LmKz1LFu2DL169cLAgQPrPe+I7xeRlJTU9gFs/+o40ntmL4rrqWvfvj1cXFxQUFBQ7/mCggKEhIQ0+JqQkJBrXl/336b8mW2hOfdW57333sPChQvx888/Iz4+/prXdurUCe3bt0d6enqLMzdGS+6rjpubGxISEmyZ5fCeteS+ysvLsXbtWjz66KPX/Tr2fr+a62rfZ76+vvDw8GiVfwfOTqntn1LbPoDtX0OU1v7Zs+1TXFGn0WjQr18/bN++3facxWLB9u3b6/1m80eDBw+udz0AbN261XZ9dHQ0QkJC6l1jMBhw4MCBq/6ZbaE59wYA77zzDl5//XVs3rwZ/fv3v+7XycnJQXFxMUJDQ1sl9/U0977+yGw2Izk52ZZZDu9ZS+7rm2++gdFoxIMPPnjdr2Pv96u5rvd91hr/DpydUts/pbZ9ANu/hiit/bNr29ekZRUOYu3atUKr1YoVK1aIkydPiieeeEL4+fnZlnxPmzZNvPTSS7brf/31V+Hq6iree+89cerUKfHqq682uKTfz89PfPfdd+L48eNi8uTJkm1p0pR7W7hwodBoNOLbb7+ttwS8tLRUCCFEaWmpeOGFF8S+fftEZmam2LZtm+jbt6/o3LmzqKqqku19zZ8/X2zZskWcPXtWJCYmiilTpgh3d3dx4sSJevcu9XvW1PuqM2zYMHHfffdd8bxc3q+6LEePHhVHjx4VAMT7778vjh49KrKysoQQQrz00kti2rRptuvrlvW/+OKL4tSpU2Lx4sUNLuu/1t8XXZ9S2z+ltn3NuTe2f9K+Z3Ju+xRZ1AkhxIcffigiIyOFRqMRAwcOFPv377d9bvjw4WL69On1rv/6669Fly5dhEajEXFxceLHH3+s93mLxSJeeeUVERwcLLRarRg1apRITU21x61coSn31rFjRwHgiserr74qhBCioqJCjBkzRgQGBgo3NzfRsWNH8fjjj0vyQ7Qp9/Xcc8/Zrg0ODha33nqrOHLkSL0/Ty7vWVP/LZ4+fVoAED///PMVf5ac3q8dO3Y0+G+r7n6mT58uhg8ffsVr+vTpIzQajejUqVO9/afqXOvvixpHqe2fUts+Idj+1XGE9k/ObZ9KCCGa1rdHRERERHKjuDl1RERERM6IRR0RERGRArCoIyIiIlIAFnVERERECsCijoiIiEgBWNQRERERKQCLOiIiIiIFYFFHREREpAAs6oiIiIgUgEUdOaznn38ed955p9QxiIjsim0fXQ2LOnJYBw8eRP/+/aWOQURkV2z76Gp49is5nOrqanh5eaGmpsb23KBBg7B//34JUxERtS22fXQ9rlIHIGoqV1dX/Prrrxg0aBCSkpIQHBwMd3d3qWMREbUptn10PSzqyOGo1Wrk5eUhICAAvXv3ljoOEZFdsO2j6+GcOnJIR48eZaNGRE6HbR9dC4s6ckhJSUls2IjI6bDto2thUUcOKTk5GX369JE6BhGRXbHto2thUUcOyWKxIDU1FXl5edDr9VLHISKyC7Z9dC0s6sghvfHGG1ixYgU6dOiAN954Q+o4RER2wbaProX71BEREREpAHvqiIiIiBSARR0RERGRArCoIyIiIlIAFnVERERECsCijoiIiEgBWNQRERERKQCLOiIiIiIFYFFHREREpAAs6oiIiIgUgEUdERERkQKwqCMiIiJSABZ1RERERArw/wG4APAEUxDXOQAAAABJRU5ErkJggg==",
+      "text/plain": [
+       "<Figure size 640x480 with 2 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "fig, axs = plt.subplots(1, 2)\n",
+    "axs[0].plot(t_eval, solution[\"y squared\"](t_eval))\n",
+    "axs[0].set_ylabel(r\"$y^2$\")\n",
+    "axs[0].set_xlabel(r\"$t$\")\n",
+    "axs[1].plot(solution.t, solution[\"y squared\"].sensitivities[\"a\"])\n",
+    "axs[1].set_ylabel(r\"$\\frac{dy^2}{da}$\")\n",
+    "axs[1].set_xlabel(r\"$t$\")\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "495bcf05",
+   "metadata": {},
+   "source": [
+    "## Sensitivities for the DFN model\n",
+    "\n",
+    "We can do the same for the DFN model included in PyBaMM. We will setup the DFN model using \"Current function\" as an input parameter. This is the parameter we wish to calculate the sensitivities with respect to."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "e9119add",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# now lets do the same for the DFN model\n",
+    "\n",
+    "# load model\n",
+    "model = pybamm.lithium_ion.DFN()\n",
+    "\n",
+    "# load parameter values and process model and geometry\n",
+    "param = model.default_parameter_values\n",
+    "\n",
+    "# we want to calculate the sensitivities of the \"Current function\" parameter, so set\n",
+    "# this an an input parameter\n",
+    "param.update({\"Current function [A]\": \"[input]\"})\n",
+    "\n",
+    "solver = pybamm.IDAKLUSolver(rtol=1e-3, atol=1e-6)\n",
+    "\n",
+    "sim = pybamm.Simulation(model, parameter_values=param, solver=solver)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f198fbfe",
+   "metadata": {},
+   "source": [
+    "We can now evaluate the senstivities of, for example, the \"Terminal voltage\" output of the model with respect to the input parameter \"Current function\"."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "1d794537",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAGwCAYAAAC5ACFFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABsoElEQVR4nO3deViUVf8G8HsGGPZddgFBVERxT8Q1g1yw1PRnoeSShmaapWZpqaW+aWaZaZaZuaVmZVppqeG+IbiACiIKmMgyICIM+zbP7w90ihCcgRmGgftzXXO98mxzH4fX+Xae85wjEgRBABERERE9lljbAYiIiIgaMxZLRERERLVgsURERERUCxZLRERERLVgsURERERUCxZLRERERLVgsURERERUC31tB2gK5HI50tLSYG5uDpFIpO04REREpARBEJCXlwdnZ2eIxTX3H7FYUoO0tDS4urpqOwYRERHVwd27d9GyZcsa97NYUgNzc3MAlX/ZFhYWWk5DREREypDJZHB1dVV8j9eExZIaPLr1ZmFhwWKJiIhIxzxpCA0HeBMRERHVgsUSERERUS1YLBERERHVgsUSERERUS1YLBERERHVgsUSERERUS3qNHWAXC5HQkICMjMzIZfLq+zr37+/WoIRERERNQYqF0vnz5/HuHHjcOfOHQiCUGWfSCRCRUWF2sIRERERaZvKxdJrr72GHj164I8//oCTkxPXQiMiIqImTeVi6datW9izZw+8vLw0kYeIiIioUVF5gLefnx8SEhI0kYWIiIio0VG5Z+mNN97A3LlzIZVK4evrCwMDgyr7O3XqpLZwRERERNomEv47SvsJxOLqnVEikQiCIDTbAd4ymQyWlpbIzc1V60K6DwpKkVNUBidLIxgZ6KntukRERKT897fKPUu3b9+uVzBS3nPrziA1pwj7Xu+Nrm7W2o5DRETULKlcLLm7u2siBz2GpbEBUnOKkFtUpu0oREREzVadJqVMTEzEmjVrEBcXBwDw8fHBm2++idatW6s1XHNnaVw5HiynkMUSERGRtqj8NNzhw4fh4+ODyMhIdOrUCZ06dUJERAQ6dOiAsLAwTWRstlqYGwIAsvJLtJyEiIio+VK5Z2n+/PmYPXs2Pv7442rb3333XTz77LNqC9fc2T8sljJkxVpOQkRE1Hyp3LMUFxeHKVOmVNs+efJkXL9+XS2hqJKLlTEAIDWnSMtJiIiImi+ViyU7OztER0dX2x4dHQ17e3t1ZKKHXG1MAADJ2YVaTkJERNR8qXwbLjQ0FFOnTkVSUhJ69+4NADh79ixWrlyJOXPmqD1gc9bKtrJYupNVqJjHioiIiBqWyj1LixYtwuLFi7Fu3ToMGDAAAwYMwJdffokPP/wQCxcu1ERGAEB2djZCQkJgYWEBKysrTJkyBfn5+bWes3HjRjz99NOwsLCASCRCTk5OtWNatWoFkUhU5fXf8Vja4mZrArEIyCspR2YeB3kTERFpg8rFkkgkwuzZs5GSkoLc3Fzk5uYiJSUFb775pkZ7PkJCQhAbG4uwsDAcOHAAp06dwtSpU2s9p7CwEEOGDMF7771X63FLly5Fenq64vXGG2+oM3qdGerroVULUwBAvDRPy2mIiIiapzrNs/SIubm5unLUKi4uDocOHcKFCxfQo0cPAMC6desQFBSETz/9FM7Ozo8976233gIAnDhxotbrm5ubw9HRUek8JSUlKCn5p6dHJpMpfa6q2jtaIOleAW5IZejf1k5j70NERESPp1TPUrdu3fDgwQMAQNeuXdGtW7caX5oQHh4OKysrRaEEAIGBgRCLxYiIiKj39T/++GPY2tqia9euWLVqFcrLy2s9fsWKFbC0tFS8XF1d652hJt6OlQVpbJrmCjIiIiKqmVI9SyNGjIChoaHizw090FgqlVZ70k5fXx82NjaQSqX1uvasWbPQrVs32NjY4Ny5c1iwYAHS09OxevXqGs9ZsGBBlcHsMplMYwWTb0tLAMC1lFyNXJ+IiIhqp1Sx9MEHHyj+/OGHH6rtzefPn4+VK1fWesyjJVU05d9FT6dOnSCRSDBt2jSsWLFCUSD+l6GhYY371K1zSysAQFJWAbILSmFjKmmQ9yUiIqJKKo9Z8vT0xIULF2Bra1tle05ODrp164akpCSlrzV37lxMmjTpie/n6OiIzMzMKtvLy8uRnZ2t0lgjZfj5+aG8vBx///032rVrp9Zr14W1qQRe9mZIyMzHpTsP8KyPg7YjERERNSsqF0t///03Kioqqm0vKSlBSkqKSteys7ODnd2TBy37+/sjJycHly5dQvfu3QEAx44dg1wuh5+fn0rv+STR0dEQi8WNaoLNHu7WSMjMx8W/s1ksERERNTCli6Xff/9d8efDhw/D0tJS8XNFRQWOHj0KDw8P9aZ7qH379hgyZAhCQ0OxYcMGlJWVYebMmQgODlY8CZeamoqAgABs374dPXv2BFA51kkqlSIhIQEAcO3aNZibm8PNzQ02NjYIDw9HREQEBg4cCHNzc4SHh2P27Nl4+eWXYW1trZG21EUvT1vsvnAXZxKytB2FiIio2VG6WBo5ciSAynmWJk6cWGWfgYEBWrVqhc8++0yt4f5t586dmDlzJgICAiAWizF69GisXbtWsb+srAzx8fEoLPxnaZANGzZgyZIlip/79+8PANiyZQsmTZoEQ0ND7N69Gx9++CFKSkrg4eGB2bNnN7qZyPt4tQBQ+URcVn4JWpg1zHgpIiIiAkSCIAiqnODh4YELFy6gRYsWmsqkc2QyGSwtLZGbmwsLCwuNvEfQF6dxPV2Gz8Z0xujuLTXyHkRERM2Jst/fKs/gffv2bRZKWhD4cKzS4dj6TZVAREREqlG5WJo1a1aV21+PfPnll4oZs0n9BneoLJZO3bqHgpLaJ80kIiIi9VG5WPrll1/Qp0+fatt79+6NPXv2qCUUVefjZAGPFqYoLpPjr+vsXSIiImooKhdL9+/fr/Ik3CMWFhbIyuLTWpoiEokwvHPlk3+/XErVchoiIqLmQ+ViycvLC4cOHaq2/eDBg/D09FRLKHq8/+veEiIRcCYhC7ezCrQdh4iIqFlQeVLKOXPmYObMmbh37x6eeeYZAMDRo0fx2WefYc2aNerOR//iamOCge3scexGJnaev4OFz/loOxIREVGTp3KxNHnyZJSUlOCjjz7CsmXLAACtWrXC119/jQkTJqg9IFX1ci83HLuRiZ8vpWDOoLYwkaj8ERIREZEKVJ5n6d/u3bsHY2NjmJmZqTOTzmmIeZYeqZALGPjpCSRnF2LRcz6Y0lczs6YTERE1dRqbZ+nf7Ozsmn2h1ND0xCJMf7o1AGDjqUSUlFdfp4+IiIjUR+V7OB4eHhCJRDXuT0pKqlcgerJR3Vyw9ugtpOcW46cLdzHev5W2IxERETVZKhdL/514sqysDFFRUTh06BDmzZunrlxUC0N9PUx/ujUW/xaLL44mYFS3ljA15NglIiIiTVD5G/bNN9987Pb169fj4sWL9Q5Eygl+yg3fnbmNO/cLsen0bbwZ2EbbkYiIiJqkeo1Z+rehQ4fil19+Udfl6Akk+mK8PagdgMqxS1n5JVpORERE1DSprVjas2cPbGxs1HU5UsIwXyd0bmmJgtIKrDoUr+04RERETZLKt+G6du1aZYC3IAiQSqW4d+8evvrqK7WGo9qJxSIsfr4DRn99Dj9evIsXn2qJ7u4sWImIiNRJ5WJp5MiRVX4Wi8Wws7PD008/DW9vb3XlIiV1d7fGiz1a4qeLKVj4ayz2z+wDfT21dRgSERE1e0oVS3PmzMGyZctgamqKgQMHwt/fHwYGBprORkqaP7Q9DsdmIC5dho2nk/D6017ajkRERNRkKNUFsW7dOuTn5wMABg4ciAcPHmg0FKnGxlSCRQ/XiVtz5BYSMvO1nIiIiKjpUKpnqVWrVli7di0GDRoEQRAQHh4Oa2vrxx7bv39/tQYk5Yzu5oL9V9Jw8uY9zNtzBT9P8+ftOCIiIjVQam24X3/9Fa+99hoyMzMhEolQ0ykikQgVFc1v+Y2GXBuuNqk5RRjy+SnklZRj3uB2mDGQt+OIiIhqota14UaOHAmpVAqZTAZBEBAfH48HDx5Ue2VnZ6utAaQ6FytjfDC8AwDg87CbiEnN1XIiIiIi3afSfRozMzMcP34cHh4esLS0fOyLtGt0NxcM7uCAcrmAN3dHoai0+fX0ERERqZPKg1oGDBgAfX2uQ9ZYiUQirBjVCfbmhki8V4D//XFd25GIiIh0GkcAN0E2phJ8/lIXiETAzohkHIqRajsSERGRzmKx1ET18WqBqf09AQDv7LmCu9mFWk5ERESkm1gsNWFvD2qHrm5WkBWXY+YPUSgtl2s7EhERkc5RuViaPHky8vLyqm0vKCjA5MmT1RKK1MNAT4x1Y7vC0tgAV+7mYNXhG9qOREREpHNULpa2bduGoqKiatuLioqwfft2tYQi9WlpbYJV/9cJAPDt6ds4Gpeh5URERES6ReliSSaTITc3F4IgIC8vDzKZTPF68OAB/vzzT9jb22syK9XRoA6OeKVPKwDAnJ+uIOUBxy8REREpS+k5AKysrCASiSASidC2bdtq+0UiEZYsWaLWcKQ+84d643JyDq7czcGMXVH4eZo/JPocskZERPQkShdLx48fhyAIeOaZZ/DLL7/AxsZGsU8ikcDd3R3Ozs4aCUn1Z6ivhy/HdsVz687gyt0cLP8zDh8+nO2biIiIaqbU2nD/dufOHbi6ukIsZq/EI41lbThlHI3LwJRtFwEAX4V0Q5Cvk5YTERERaYey398qT8Xt7u6OnJwcREZGIjMzE3J51cfRJ0yYoHpaajAB7R3w2oDW2HAyEe/suQpvR3N42plpOxYREVGjpXLP0v79+xESEoL8/HxYWFhAJBL9czGRqFkupqtLPUsAUF4hx7hNEYi8nQ1vR3Pse70PjCV62o5FRETUoJT9/lb5XtrcuXMxefJk5OfnIycnBw8ePFC8NFkoZWdnIyQkBBYWFrCyssKUKVOQn59f6/FvvPEG2rVrB2NjY7i5uWHWrFnIzc2tclxycjKGDRsGExMT2NvbY968eSgvL9dYOxoDfT0xvhzbFS3MDHFDmof3912DijUzERFRs6FysZSamopZs2bBxMREE3lqFBISgtjYWISFheHAgQM4deoUpk6dWuPxaWlpSEtLw6effoqYmBhs3boVhw4dwpQpUxTHVFRUYNiwYSgtLcW5c+ewbds2bN26FYsXL26IJmmVvYURvhzXFXpiEfZGpWJXZLK2IxERETVKKt+GGzVqFIKDg/Hiiy9qKlM1cXFx8PHxwYULF9CjRw8AwKFDhxAUFISUlBSln8L7+eef8fLLL6OgoAD6+vo4ePAgnnvuOaSlpcHBwQEAsGHDBrz77ru4d+8eJBKJUtfVtdtw/7bhZCI+PngDEj0xfn7NH51drbQdiYiIqEFobID3sGHDMG/ePFy/fh2+vr4wMDCosn/48OGqp32C8PBwWFlZKQolAAgMDIRYLEZERAReeOEFpa7z6C9DX19fcV1fX19FoQQAgwcPxvTp0xEbG4uuXbs+9jolJSUoKSlR/CyTyerSrEZhWn9PXL7zAH9dz8DrOy9j/xt9YWOqXJFIRETUHKhcLIWGhgIAli5dWm2fSCRCRUVF/VP9h1QqrTY7uL6+PmxsbCCVSpW6RlZWFpYtW1bl1p1UKq1SKAFQ/FzbdVesWNFkJuAUiUT49MXOGL7uDP6+X4g3d0dh6ys9oScWPflkIiKiZkDlMUtyubzGl6qF0vz58xWzgtf0unGj/ou/ymQyDBs2DD4+Pvjwww/rfb0FCxYgNzdX8bp79269r6lNFkYG2DC+O4wMxDh9Kwufh93UdiQiIqJGQ+WepX8rLi6GkZFRnc+fO3cuJk2aVOsxnp6ecHR0RGZmZpXt5eXlyM7OhqOjY63n5+XlYciQITA3N8e+ffuq3DZ0dHREZGRkleMzMjIU+2piaGgIQ0PDWt9X13g7WmDl6E54c3c0vjyegM6uVnjWx+HJJxIRETVxKvcsVVRUYNmyZXBxcYGZmRmSkpIAAIsWLcJ3332n0rXs7Ozg7e1d60sikcDf3x85OTm4dOmS4txjx45BLpfDz8+vxuvLZDIMGjQIEokEv//+e7XCzt/fH9euXatSiIWFhcHCwgI+Pj4qtaUpGNHFBZN6twIAzPkxGrezCrQbiIiIqBFQuVj66KOPsHXrVnzyySdVnhbr2LEjNm3apNZwj7Rv3x5DhgxBaGgoIiMjcfbsWcycORPBwcGKJ+FSU1Ph7e2t6Cl6VCgVFBTgu+++g0wmg1QqhVQqVdwuHDRoEHx8fDB+/HhcuXIFhw8fxsKFCzFjxowm13OkrPeC2qOHuzXySsoxfcclFJY27TmniIiInkTlYmn79u3YuHEjQkJCoKf3z6zPnTt3Vsv4oprs3LkT3t7eCAgIQFBQEPr27YuNGzcq9peVlSE+Ph6FhYUAgMuXLyMiIgLXrl2Dl5cXnJycFK9HY4z09PRw4MAB6Onpwd/fHy+//DImTJjw2MHrzYVEX4z1Id0UE1Yu2MsJK4mIqHlTeZ4lY2Nj3LhxA+7u7jA3N8eVK1fg6emJ69evo2fPnrXOqt1U6fI8SzWJSLqPcZsiUCEXsGR4B0x8eHuOiIioqdDYcic+Pj44ffp0te179uypcV4i0j1+nrZYMNQbALDswHVc/Lv5rflHREQE1OFpuMWLF2PixIlITU2FXC7H3r17ER8fj+3bt+PAgQOayEhaMqWvB6Lu5uCPq+mYvvMy/nijL+wt6v70IxERkS5SuWdpxIgR2L9/P44cOQJTU1MsXrwYcXFx2L9/P5599llNZCQtEYlE+GR0J7R1MMO9vBJM33kZpeVybcciIiJqUCoVS+Xl5Vi6dCk8PDwQFhaGzMxMFBYW4syZMxg0aJCmMpIWmRrq45vxPWBupI9Ldx5g2YHr2o5ERETUoFQqlvT19fHJJ5+gvJyPkzcnHi1MsealLgCA78/fwc8XdXvGciIiIlWofBsuICAAJ0+e1EQWasQC2jvgrcA2AID3f43BtZRcLSciIiJqGCoP8B46dCjmz5+Pa9euoXv37jA1Na2yf/jw4WoLR43LrGfaICY1F0fiMvHajkv4fWYf2Jo1z8k7iYio+VB5niWxuObOKJFIpPJiuk1BU5xnqSay4jKM+PIsbmcVoHdrW2yf3BP6eip3UBIREWmdxuZZksvlNb6aY6HU3FgYGeCb8d1hItHDucT7+ORwvLYjERERaZRKxVJZWRn09fURExOjqTykA9o6mOPTMZ0BABtPJWH/lTQtJyIiItIclYolAwMDuLm5sQeJEOTrhNcGtAYAzNtzBdfTZFpOREREpBkq34Z7//338d577yE7m8tfNHfzBrdDvzYtUFwmx9TvL+JBQam2IxEREamdygO8u3btioSEBJSVlcHd3b3a03CXL19Wa0Bd0JwGeP9XTmEpRqw/izv3C9HHyxbbXuGAbyIi0g3Kfn+rPHXAyJEj65OLmhgrEwk2ju+BF746i7MJ9/HxwRtY+JyPtmMRERGpjco9S1Rdc+5ZeuRQTDpe21HZq/j5S53xQteWWk5ERERUO41NHUD0OEM6OuGNZ7wAAPN/ucYZvomIqMlQuVgSi8XQ09Or8UXN1+zAtgjwtkdJuRzTvr+IrPwSbUciIiKqN5XHLO3bt6/Kz2VlZYiKisK2bduwZMkStQUj3SMWi/B5cBeMXH8WSfcK8PrOy9j5qh8MOOCbiIh0mNrGLO3atQs//vgjfvvtN3VcTqdwzFJVCZn5GLn+LPJLyjHR3x1LRnTUdiQiIqJqGnzMUq9evXD06FF1XY50mJe9Gda81AUAsC38Dn66cFe7gYiIiOpBLcVSUVER1q5dCxcXF3VcjpqAQB8HzA5sCwBY+GsMopIfaDkRERFR3ag8Zsna2hoikUjxsyAIyMvLg4mJCXbs2KHWcKTb3njGC7FpufjregZe23EJ+2f2hb2FkbZjERERqUTlMUvbtm2r8rNYLIadnR38/PxgbW2t1nC6gmOWapZfUo4X1p/Frcx8dHOzwg9Te8FQn09NEhGR9in7/c1JKdWAxVLtbmcVYPiXZ5BXXI6xPd2wYpSvtiMRERGpf4D3rVu3MHbsWMhk1VeXz83Nxbhx45CUlFS3tNSkebQwxdqxXSESAT9EJmNnxB1tRyIiIlKa0sXSqlWr4Orq+tjKy9LSEq6urli1apVaw1HTMbCdPeYNbgcA+PD3WFz4O1vLiYiIiJSjdLF08uRJjBkzpsb9L774Io4dO6aWUNQ0TR/QGsN8nVBWIWD6jstIzy3SdiQiIqInUrpYSk5Ohr29fY37W7Rogbt3OZ8O1UwkEmHVmE7wdjRHVn4Jpm6/hKLSCm3HIiIiqpXSxZKlpSUSExNr3J+QkMDBzfREJhJ9fDuhB2xMJbiWmot3frkKPmNARESNmdLFUv/+/bFu3boa969duxb9+vVTSyhq2lxtTPBVSDfoi0XYfyUNX52ouQgnIiLSNqWLpQULFuDgwYP4v//7P0RGRiI3Nxe5ubmIiIjA6NGjcfjwYSxYsECTWakJ6eVpiyUjOgAAVh2Ox1+xUi0nIiIiejyV5lk6cOAAJk+ejPv371fZbmtri02bNmH48OFqD6gLOM9S3S3+LQbbw+/AVKKHX17vDW9H/v0REVHD0NiklEVFRTh06BASEhIgCALatm2LQYMGwcTEpN6hdRWLpborq5BjwneRCE+6D1cbY/w2oy9sTCXajkVERM0AZ/BuQCyW6udBQSlGrD+L5OxC9PK0wfdT/GCgp5Y1nomIiGqk9hm8iTTF2lSCTRN7wMxQH+eTsrFkf6y2IxERESnoTLGUnZ2NkJAQWFhYwMrKClOmTEF+fn6tx7/xxhto164djI2N4ebmhlmzZiE3N7fKcSKRqNpr9+7dmm4O/UdbB3OseakLRCJgx/lkfH+eS6IQEVHjoDPFUkhICGJjYxEWFoYDBw7g1KlTmDp1ao3Hp6WlIS0tDZ9++iliYmKwdetWHDp0CFOmTKl27JYtW5Cenq54jRw5UoMtoZoE+jgolkRZ8nsswhPvP+EMIiIizdOJMUtxcXHw8fHBhQsX0KNHDwDAoUOHEBQUhJSUFDg7Oyt1nZ9//hkvv/wyCgoKoK+vD6CyZ2nfvn0qFUglJSUoKSlR/CyTyeDq6soxS2ogCALe+jEav0WnwdrEAL/N6As32+b78AAREWmOxsYs6enpITMzs9r2+/fvQ09PT9XLKSU8PBxWVlaKQgkAAgMDIRaLERERofR1Hv1lPCqUHpkxYwZatGiBnj17YvPmzU+cUXrFihWwtLRUvFxdXVVrENVIJBJh5ehO6NzSEg8KyxC6/SLyS8q1HYuIiJoxlYulmgqJkpISSCSaeeRbKpVWW5dOX18fNjY2kEqVm8wwKysLy5Ytq3brbunSpfjpp58QFhaG0aNH4/XXX691pnKgcoLOR5Ny5ubmck08NTMy0MM343vA3twQ8Rl5eGt3NOTyRt8BSkRETZT+kw+ptHbtWgCV/+W/adMmmJmZKfZVVFTg1KlT8Pb2VunN58+fj5UrV9Z6TFxcnErXfByZTIZhw4bBx8cHH374YZV9ixYtUvy5a9euKCgowKpVqzBr1qwar2doaAhDQ8N656KaOVoaYeOEHnjxm3AcicvA6rCbePvheCYiIqKGpHSx9PnnnwOo7FnasGFDlVtuEokErVq1woYNG1R687lz52LSpEm1HuPp6QlHR8dqt/7Ky8uRnZ0NR0fHWs/Py8vDkCFDYG5ujn379sHAwKDW4/38/LBs2TKUlJSwINKyLq5WWDnaF7N/vIIvjyegraM5hndWbnwaERGRuihdLN2+fRsAMHDgQOzbtw9WVlb1fnM7OzvY2dk98Th/f3/k5OTg0qVL6N69OwDg2LFjkMvl8PPzq/E8mUyGwYMHw9DQEL///juMjIye+F7R0dGwtrZmodRIvNC1JW5I8/DNySTM+/kKPGxN4dvSUtuxiIioGVFpzFJZWRmSk5ORnp6uqTyP1b59ewwZMgShoaGIjIzE2bNnMXPmTAQHByuehEtNTYW3tzciIyMBVBZKgwYNQkFBAb777jvIZDJIpVJIpVJUVFQAAPbv349NmzYhJiYGCQkJ+Prrr7F8+XK88cYbDdo+qt07g73xjLc9SsrlCN1+ERmyYm1HIiKiZkSlYsnAwADFxdr5otq5cye8vb0REBCAoKAg9O3bFxs3blTsLysrQ3x8PAoLCwEAly9fRkREBK5duwYvLy84OTkpXo8GZBsYGGD9+vXw9/dHly5d8M0332D16tX44IMPtNJGejw9sQhfBHeBl70ZpLJihG6/iKLSCm3HIiKiZkLleZaWL1+OmzdvYtOmTdUewW+uuDZcw0i+X4gR68/gQWEZgnwd8eXYbhCLRdqORUREOkpjC+m+8MILOHr0KMzMzODr6wtTU9Mq+/fu3Vu3xDqMxVLDufB3NkK+jUBphRwzB3rxCTkiIqozZb+/Ve4asrKywujRo+sVjqiunmplgxWjfDH358on5DztTDGqW0ttxyIioiZM5WJpy5YtmshBpLTR3VsiKSsf648nYv4v1+BqY4KnWtloOxYRETVRKs/gvXnzZsU0AkTaMvfZdhja0RGlFXJM+/4Sku8XajsSERE1USoXSytWrICXlxfc3Nwwfvx4bNq0CQkJCZrIRlQjsViE1S92ga+LJbILSjF52wXIisu0HYuIiJoglYulW7duITk5GStWrICJiQk+/fRTtGvXDi1btsTLL7+siYxEj2Us0cOmiT3gaGGEhMx8zNh5GeUVcm3HIiKiJkblp+H+rbCwEKdPn8YPP/yAnTt3QhAElJc3vxXi+TScdsWk5mLMhnAUlVVggr87lo7oqO1IRESkA5T9/la5Z+mvv/7Ce++9h969e8PW1hYLFiyAtbU19uzZg3v37tUrNFFddHSxxJrgLhCJgO3hd7Dt3N/ajkRERE2Iyj1LYrEYdnZ2mDt3LqZOnaqWNeJ0HXuWGodvTiZixcEbEIuAzZOewtPt7LUdiYiIGjGN9SytXr0affr0wSeffIIOHTpg3Lhx2LhxI27evFmvwET1NbW/J17s0RJyAZi5Kwrx0jxtRyIioiagXmOWrl27hpMnT+LYsWM4cOAA7O3tkZKSos58OoE9S41Habkc47+LQMTtbLS0NsavM/qghZmhtmMREVEjpLGeJQAQBAGXL19GWFgYDh8+jOPHj0Mul8POzq7OgYnUQaIvxoaXu6OVrQlSHhRh2veXUFzGRXeJiKjuVC6Wnn/+edja2qJnz57YuXMn2rZti23btiErKwtRUVGayEikEmtTCb6b9BQsjPRx6c4DzP/lKurRgUpERM2cysudeHt7Y9q0aejXrx8sLS01kYmo3lrbmeHrl7tj4uZI/BqdBk87M8wKaKPtWEREpIPqNWaJKnHMUuP1Q2QyFuy9BgBY81IXjOzqouVERETUWGh0zBKRrhjb0w1T+3sCAObtuYLwxPtaTkRERLqGxRI1efOHeCPI1xFlFQKmfX8RCZmcUoCIiJTHYomavEeL7nZzs4KsuByTtlzAvbwSbcciIiIdwWKJmgUjAz1smviUYkqBKdsuoLC0+a1jSEREqqtTsZSYmIiFCxdi7NixyMzMBAAcPHgQsbGxag1HpE42phJsfaUnrE0McDUlF7N+iEKFnM83EBFR7VQulk6ePAlfX19ERERg7969yM/PBwBcuXIFH3zwgdoDEqlTqxam2DSxByT6YhyJy8TS/bGcg4mIiGqlcrE0f/58/O9//0NYWBgkEoli+zPPPIPz58+rNRyRJnR3t8Gal7oAALaF38F3Z25rNxARETVqKhdL165dwwsvvFBtu729PbKystQSikjTgnyd8F6QNwDgoz/jcPBaupYTERFRY6VysWRlZYX09OpfLFFRUXBx4YR/pDtC+3lifC93CALw1o/RuHTngbYjERFRI6RysRQcHIx3330XUqkUIpEIcrkcZ8+exdtvv40JEyZoIiORRohEInzwvA8CvO1RUi5H6PaL+DurQNuxiIiokVG5WFq+fDm8vb3h6uqK/Px8+Pj4oH///ujduzcWLlyoiYxEGqOvJ8a6cV3h62KJ7IJSvLL1Ah4UlGo7FhERNSJ1XhsuOTkZMTExyM/PR9euXdGmTfNdpJRrw+m+TFkxXvjqHFJzitDD3Ro7XvWDkYGetmMREZEGKfv9zYV01YDFUtNwMyMPo78+h7zicgzr5IR1wV0hFou0HYuIiDRE2e9vfVUvPGfOnMduF4lEMDIygpeXF0aMGAEbGxtVL02kVW0dzPHN+O6YuDkSf1xNR0trYywY2l7bsYiISMtU7lkaOHAgLl++jIqKCrRr1w4AcPPmTejp6cHb2xvx8fEQiUQ4c+YMfHx8NBK6sWHPUtOy93IK5vx0BQCwbGRHjO/lruVERESkCcp+f6s8wHvEiBEIDAxEWloaLl26hEuXLiElJQXPPvssxo4di9TUVPTv3x+zZ8+uVwOItGVUt5aY82xbAMDi32LwJ+dgIiJq1lTuWXJxcUFYWFi1XqPY2FgMGjQIqampuHz5MgYNGtRsJqlkz1LTIwgC3v81BrsikiHRE2PzpKfQt00LbcciIiI10ljPUm5urmLx3H+7d+8eZDIZgMqJK0tL+fg16S6RSIRlIzoiyNcRpRVyTP3+Iq7czdF2LCIi0oI63YabPHky9u3bh5SUFKSkpGDfvn2YMmUKRo4cCQCIjIxE27Zt1Z2VqEHpiUX4/KUu6OvVAoWlFZi0JRIJmfnajkVERA1M5WLpm2++QUBAAIKDg+Hu7g53d3cEBwcjICAAGzZsAAB4e3tj06ZNag2anZ2NkJAQWFhYwMrKClOmTEF+fu1fXNOmTUPr1q1hbGwMOzs7jBgxAjdu3KhyTHJyMoYNGwYTExPY29tj3rx5KC8vV2t20l2G+nrYML47Ore0xIPCMkz4LgJpOUXajkVERA2ozvMs5efnIykpCQDg6ekJMzMztQb7r6FDhyI9PR3ffPMNysrK8Morr+Cpp57Crl27ajxn48aN8Pb2hpubG7Kzs/Hhhx8iOjoat2/fhp6eHioqKtClSxc4Ojpi1apVSE9Px4QJExAaGorly5crnY1jlpq+7IJSjNlwDon3CtDazhQ/v9YbNqYSbcciIqJ6aFKTUsbFxcHHxwcXLlxAjx49AACHDh1CUFAQUlJS4OzsrNR1rl69is6dOyMhIQGtW7fGwYMH8dxzzyEtLQ0ODg4AgA0bNuDdd9/FvXv3IJEo92XIYql5SM0pwv99fQ7pucXo3NISO0N7wcxQ5anKiIiokdDYAG8AuHjxIt555x0EBwdj1KhRVV6aEB4eDisrK0WhBACBgYEQi8WIiIhQ6hoFBQXYsmULPDw84Orqqriur6+volACgMGDB0MmkyE2NrbGa5WUlEAmk1V5UdPnYmWM76f0hLWJAa6k5OK17y+hpLxC27GIiEjDVC6Wdu/ejd69eyMuLg779u1DWVkZYmNjcezYMVhaWmoiI6RSKezt7ats09fXh42NDaRSaa3nfvXVVzAzM4OZmRkOHjyIsLAwRY+RVCqtUigBUPxc23VXrFgBS0tLxetR8UVNn5e9Oba80hMmEj2cScjCnB+voELe6DtniYioHlQulpYvX47PP/8c+/fvh0QiwRdffIEbN27gxRdfhJubm0rXmj9/PkQiUa2v/w7IVlVISAiioqJw8uRJtG3bFi+++CKKi4vrdc0FCxYgNzdX8bp79269rke6pYurFTaO7wEDPRH+uJaORb/FQAfuZhMRUR2pPOAiMTERw4YNAwBIJBIUFBRAJBJh9uzZeOaZZ7BkyRKlrzV37lxMmjSp1mM8PT3h6OhYbW6n8vJyZGdnw9HRsdbzH/X+tGnTBr169YK1tTX27duHsWPHwtHREZGRkVWOz8jIAIBar2toaAhDQ8Na35eatr5tWmDNS10x84fL2BWRDFtTCeYOaqftWEREpAEqF0vW1tbIy8sDUDmbd0xMDHx9fZGTk4PCwkKVrmVnZwc7O7snHufv74+cnBxcunQJ3bt3BwAcO3YMcrkcfn5+Sr+fIAgQBAElJSWK63700UfIzMxU3OYLCwuDhYVFs1nXjupuWCcn5BR1xPv7YrDuWAJsTCV4pY+HtmMREZGaqXwbrn///ggLCwMAjBkzBm+++SZCQ0MxduxYBAQEqD0gALRv3x5DhgxBaGgoIiMjcfbsWcycORPBwcGKJ+FSU1Ph7e2t6ClKSkrCihUrcOnSJSQnJ+PcuXMYM2YMjI2NERQUBAAYNGgQfHx8MH78eFy5cgWHDx/GwoULMWPGDPYckVJC/Nwx9+E6ckv2X8evUalaTkREROqmcs/Sl19+qRjz8/7778PAwADnzp3D6NGjsXDhQrUHfGTnzp2YOXMmAgICIBaLMXr0aKxdu1axv6ysDPHx8YreLSMjI5w+fRpr1qzBgwcP4ODggP79++PcuXOKXiQ9PT0cOHAA06dPh7+/P0xNTTFx4kQsXbpUY+2gpmfmM17ILizFlrN/4+2fr8DS2AADve2ffCIREekEnZhnqbHjPEsklwuY+/MV7ItKhZGBGDum+KFHKxttxyIiolpobJ4lPT29xy6ke//+fejp6al6OaImQSwW4ZP/64SB7exQXCbH5K0XcEPK+beIiJoClYulmjqiSkpKlJ7xmqgpMtAT46uQ7ujhbg1ZcTkmfBeJu9mqPfRARESNj9Jjlh6NDxKJRNi0aVOVteAqKipw6tQpeHt7qz8hkQ4xlujhu4lP4aWN4bghzcPL30Vgz2u9YWfOBwaIiHSV0mOWPDwqH4m+c+cOWrZsWeWWm0QiQatWrbB06VKVHuVvKjhmif4rQ1aM/9twDnezi+DjZIHd03rBwshA27GIiOhfNLaQ7sCBA7F3715YW1vXO2RTwWKJHufvrAL834ZwZOWXoKeHDbZP7gkjA47rIyJqLDQ2wPv48eMslIiU0KqFKbZNfgrmhvqIvJ2NmbuiUF4h13YsIiJSkVI9S3PmzFH6gqtXr65XIF3EniWqTUTSfUzYHImScjn+r3tLrPq/ThCJRNqORUTU7Cn7/a3UAO+oqCil3pRfAETV+Xna4stx3fDajkvYcykFNqYSvBfUXtuxiIhISZyUUg3Ys0TK+PniXczbcxUAMH+oN14b0FrLiYiImjeNjVn6t5SUFKSkpNTnEkTNxpgernj/YY/Sxwdv4McLyVpOREREylC5WJLL5Vi6dCksLS3h7u4Od3d3WFlZYdmyZZDLOXiVqDah/T0VPUoL9l7DoRiplhMREdGTqFwsvf/++/jyyy/x8ccfIyoqClFRUVi+fDnWrVuHRYsWaSIjUZPy7pB2eKmHK+QCMGt3FM4lZmk7EhER1ULlMUvOzs7YsGEDhg8fXmX7b7/9htdffx2pqalqDagLOGaJVFVeIceMXZdxODYDZob62D21Fzq6WGo7FhFRs6KxMUvZ2dmPXdbE29sb2dnZql6OqFnS1xPji+Cu8Pe0RX5JOSZujkTSvXxtxyIiosdQuVjq3Lkzvvzyy2rbv/zyS3Tu3FktoYiaAyMDPWyc0B0dXSxwv6AU47+LRHpukbZjERHRf6h8G+7kyZMYNmwY3Nzc4O/vDwAIDw/H3bt38eeff6Jfv34aCdqY8TYc1UdWfgnGbAjH7awCuNua4IfQXnC2MtZ2LCKiJk9jt+EGDBiA+Ph4vPDCC8jJyUFOTg5GjRqF+Pj4ZlkoEdVXCzND7HjVD642xrhzvxDBG88jLYc9TEREjQUnpVQD9iyROqTmFGHsxvNIzi6Em40JfpjaCy7sYSIi0hiN9Sx5eXnhww8/xK1bt+oVkIiqcrEyxu6pveBua4Lk7EIEbwxHKnuYiIi0TuViacaMGfjjjz/Qrl07PPXUU/jiiy8glXJiPSJ1cP5XwXQ3uwjBG8OR8qBQ27GIiJo1lYul2bNn48KFC7hx4waCgoKwfv16uLq6YtCgQdi+fbsmMhI1K06WlQVTK0XBdB53s1kwERFpi1rGLJ0/fx7Tp0/H1atXUVFRoY5cOoVjlkgTpLnFCN4Yjr/vFypu0bnamGg7FhFRk9EgC+lGRkbirbfewgsvvICbN29izJgx9bkcEf2Lo6URdk/1h0cLU6TmsIeJiEhbVC6Wbt68iQ8++ABt27ZFnz59EBcXh5UrVyIjIwO7d+/WREaiZsvR0gg/hPZiwUREpEUq34YTi8V46qmnMG7cOAQHB8PBwUFT2XQGb8ORpmXIijF243kkZRXAxcoYP4T2gpstb8kREdWHst/fKhdLt27dQps2beodsClhsUQN4d8Fk/PDW3QsmIiI6k5jY5ZYKBFph4OFEXZP7QVPO1OkPRz8fed+gbZjERE1efUa4E1EDcvewgi7Q3uhtaJgOs+CiYhIw1gsEekYewsj/DC1smBKzy3GS9+cx99ZLJiIiDSFxRKRDrI3rxyz5GVvBqmssofpNgsmIiKNYLFEpKPszA3xQ2gvtFEUTOEsmIiINECpp+HmzJmj9AVXr15dr0C6iE/DkTbdyyvBuG/P41ZmPhwsKgsoTzszbcciImr0lP3+1lfmYlFRUUq9qUgkUi4dEamNnbkhfpjaC+O+PY+bGfkI3nj+4VNzLJiIiNRBLWvDNXfsWaLGICu/BCHfRiA+Iw/2Dwuo1iyYiIhq1CBrwxFR49HCzBC7Qv3QzsEcmXklGLvxPBLv5Ws7FhGRzqtTsXTx4kW88847CA4OxqhRo6q8NCU7OxshISGwsLCAlZUVpkyZgvz82r8Ipk2bhtatW8PY2Bh2dnYYMWIEbty4UeUYkUhU7cU17khX2T4smLwdKwum4I3nkZDJgomIqD5ULpZ2796N3r17Iy4uDvv27UNZWRliY2Nx7NgxWFpaaiIjACAkJASxsbEICwvDgQMHcOrUKUydOrXWc7p3744tW7YgLi4Ohw8fhiAIGDRoECoqKqoct2XLFqSnpyteI0eO1Fg7iDTN1swQO1+tLJjuKQqmPG3HIiLSWSqPWerUqROmTZuGGTNmwNzcHFeuXIGHhwemTZsGJycnLFmyRO0h4+Li4OPjgwsXLqBHjx4AgEOHDiEoKAgpKSlwdnZW6jpXr15F586dkZCQgNatWwOo7Fnat2+fSgVSSUkJSkpKFD/LZDK4urpyzBI1KtkFpRj37XnckOahhZkhfgj1QxsHc23HIiJqNDQ2ZikxMRHDhg0DAEgkEhQUFEAkEmH27NnYuHFj3RPXIjw8HFZWVopCCQACAwMhFosRERGh1DUKCgqwZcsWeHh4wNXVtcq+GTNmoEWLFujZsyc2b96MJ9WPK1asgKWlpeL13+sRNQY2phLsCu2F9k4WyMovwdhvz+NWBnuYiIhUpXKxZG1tjby8yn9wXVxcEBMTAwDIyclBYWGhetM9JJVKYW9vX2Wbvr4+bGxsIJVKaz33q6++gpmZGczMzHDw4EGEhYVBIpEo9i9duhQ//fQTwsLCMHr0aLz++utYt25drddcsGABcnNzFa+7d+/WvXFEGmRjKsGuV/3g42SBrPxSjP32PG6yYCIiUonKxVL//v0RFhYGABgzZgzefPNNhIaGYuzYsQgICFDpWvPnz3/sAOt/v/47IFtVISEhiIqKwsmTJ9G2bVu8+OKLKC4uVuxftGgR+vTpg65du+Ldd9/FO++8g1WrVtV6TUNDQ1hYWFR5ETVW1qYS7HzVDx2cHxZMG1kwERGpQuUxS9nZ2SguLoazszPkcjk++eQTnDt3Dm3atMHChQthbW2t9LXu3buH+/fv13qMp6cnduzYgblz5+LBgweK7eXl5TAyMsLPP/+MF154Qan3Ky0thbW1NTZt2oSxY8c+9pg//vgDzz33HIqLi2FoaKjUdTnPEumCnMJShGyKQGyaDLYPb9G1c+QYJiJqvtQ6g/e/2djYKP4sFosxf/78uiUEYGdnBzs7uyce5+/vj5ycHFy6dAndu3cHABw7dgxyuRx+fn5Kv58gCBAEocrg7P+Kjo6GtbW10oUSka6wMqnsYXr5uwjEpMow9tvz2DzpKXRxtdJ2NCKiRk3lYgkA5HI5EhISkJmZCblcXmVf//791RLs39q3b48hQ4YgNDQUGzZsQFlZGWbOnIng4GDFk3CpqakICAjA9u3b0bNnTyQlJeHHH3/EoEGDYGdnh5SUFHz88ccwNjZGUFAQAGD//v3IyMhAr169YGRkhLCwMCxfvhxvv/222ttA1BhYmUiwY4ofJmyOxNWUXIzdeB5fhXTDQG/7J59MRNRMqVwsnT9/HuPGjcOdO3eqPTUmEomqzWGkLjt37sTMmTMREBAAsViM0aNHY+3atYr9ZWVliI+PVwwyNzIywunTp7FmzRo8ePAADg4O6N+/P86dO6cYLG5gYID169dj9uzZEAQBXl5eWL16NUJDQzXSBqLGwMqk8hbc9B2XcPpWFl7dfhEfj/LFmB58qpOI6HFUHrPUpUsXtG3bFkuWLIGTk1O1xXM1OTFlY8UxS6SLSsvlmP/LVeyNSgUAvD2oLWYM9OKC2ETUbGhszNKtW7ewZ88eeHl51SsgEWmXRF+Mz17sDHsLI2w4mYhP/7qJDFkJPhzeAXpiFkxERI+oPHWAn58fEhISNJGFiBqYSCTC/KHe+OB5H4hEwPfn72DGzssoLtPM7XQiIl2kcs/SG2+8gblz50IqlcLX1xcGBgZV9nfq1Elt4YioYbzSxwN25oaY8+MVHIqVYsJ3kfh2Qg9Ymhg8+WQioiZO5TFLYnH1ziiRSARBEDQ6wLsx45glairCE+9j6vaLyCspR1sHM2x9pSecrYy1HYuISCOU/f5WuVi6c+dOrfvd3d1VuVyTwGKJmpK4dBkmbo5EZl4JnCyNsG1yT7TlArxE1ARprFii6lgsUVOT8qAQEzdHIvFeASyM9LFp4lPo6WHz5BOJiHSIWoul33//HUOHDoWBgQF+//33Wo8dPny46ml1HIslaooeFJRiyrYLuJycA4m+GGuDu2BIRydtxyIiUhu1FktisRhSqRT29vaPHbOkuBjHLLFYoialqLQCb/wQhSNxGRCJgKXDO2C8fyttxyIiUgtlv7+VmjpALpcrZr2Wy+U1vppjoUTUlBlL9LDh5W4Y29MNggAs+i0Wqw7fqDZ7PxFRU6byPEtE1Lzo64mx/IWOmB3YFgCw/ngi3tlzFWUV8iecSUTUNNRpId0LFy7g+PHjj11Id/Xq1WoJRkSNh0gkwpuBbWBvYYj3913Dz5dScC+/BF+FdIOJpE7/jBAR6QyV/5Vbvnw5Fi5ciHbt2sHBwaHKOlJcU4qoaRvb0w12ZoaY+cNlnIi/h7HfRmDzxB6wNTPUdjQiIo1ReeoABwcHrFy5EpMmTdJQJN3DAd7U3Fy68wBTtl1ATmEZPFqYYtsrPeFma6LtWEREKlHrAO8qJ4jF6NOnT73CEZFu6+5ujT2v9YaLlTFuZxVg1NfnEJOaq+1YREQaoXKxNHv2bKxfv14TWYhIh3jZm2Hv673h7WiOrPwSvPRNOE7fuqftWEREaqfybTi5XI5hw4bh5s2b8PHxqbaQ7t69e9UaUBfwNhw1Z7LiMkzbfgnhSfehLxbh0zGdMbKri7ZjERE9kcZuw82aNQvHjx9H27ZtYWtrC0tLyyovImpeLIwMsHXyUxjWyQnlcgFv/RiNb08laTsWEZHaqNyzZG5ujt27d2PYsGGayqRz2LNEBMjlAv73Rxw2n70NAJjS1wPvB7WHWMynZImocdJYz5KNjQ1at25dr3BE1PSIxSIseq493gvyBgB8d+Y23vwxGiXlnNmfiHSbysXShx9+iA8++ACFhYWayENEOkwkEmFq/9b4/KXO0BeLsP9KGl7ZcgF5xWXajkZEVGcq34br2rUrEhMTIQgCWrVqVW2A9+XLl9UaUBfwNhxRdadu3sP0HZdQUFqB9k4W2PrKU3CwMNJ2LCIiBWW/v1WewXvkyJH1yUVEzUT/tnb4cZo/Jm2JRFy6DMO/PIOvQrqju7u1tqMREalEpWKpvLwcIpEIkydPRsuWLTWViYiaiI4ulvhlem9M2XYRCZn5CN4YjiXDO2Kcn5u2oxERKU2lMUv6+vpYtWoVysvLNZWHiJoYd1tT/DqjD4Z0cERZhYD39l3D/F+ucuA3EekMlQd4P/PMMzh58qQmshBRE2VmqI+vX+6GeYPbQSQCdl+4i5e+OY/03CJtRyMieiKVxywNHToU8+fPx7Vr19C9e3eYmppW2T98+HC1hSOipkMkEmHGQC90cLbArB+iEH03B8+vO4P147rBz9NW2/GIiGqk8tNwYnHNnVEikQgVFc2va51PwxGp5s79Akz7/hJuSPOgLxZh4bD2mNi7FUQiTmBJRA1HY5NSyuXyGl/NsVAiItW525pi7+u98XxnZ5TLBXy4/zrm/nwFxWX8N4SIGh+Vi6V/Ky4uVlcOImpmTCT6WBvcpXJJFBGw93Iq/m/DOaQ84IS3RNS4qFwsVVRUYNmyZXBxcYGZmRmSkioXzFy0aBG+++47tQckoqZLJBIhtL8ndkzxg42pBDGpMjy/7gzOJWRpOxoRkYLKxdJHH32ErVu34pNPPoFEIlFs79ixIzZt2qTWcETUPPT2aoHfZ/ZBRxcLPCgsw8vfReDbU0lQcUglEZFGqFwsbd++HRs3bkRISAj09PQU2zt37owbN26oNRwRNR8trU2w57XeGNXNBXIB+OjPOMzaHY3CUs7rRkTapXKxlJqaCi8vr2rb5XI5ysq4WCYR1Z2RgR4+G9MZS4Z3UCzEO+qrc7hzv0Db0YioGVO5WPLx8cHp06erbd+zZw+6du2qllBE1HyJRCJM7N0Ku0J7oYWZBDekeXh+3RmciM/UdjQiaqaULpYmT56MvLw8LF68GDNnzsTKlSshl8uxd+9ehIaG4qOPPsLixYs1FjQ7OxshISGwsLCAlZUVpkyZgvz8fKXOFQQBQ4cOhUgkwq+//lplX3JyMoYNGwYTExPY29tj3rx5XM6FqBHo6WGDA2/0QxdXK8iKy/HK1gtYfzyB45iIqMEpXSxt27YNRUVFGDFiBPbv348jR47A1NQUixcvRlxcHPbv349nn31WY0FDQkIQGxuLsLAwHDhwAKdOncLUqVOVOnfNmjWPneyuoqICw4YNQ2lpKc6dO4dt27Zh69atGi36iEh5jpZG+HFaL4zt6QpBAFYdjsf0HZeRX8L/oCGihqP0DN5isRhSqRT29vaazlRNXFwcfHx8cOHCBfTo0QMAcOjQIQQFBSElJQXOzs41nhsdHY3nnnsOFy9ehJOTE/bt24eRI0cCAA4ePIjnnnsOaWlpcHBwAABs2LAB7777Lu7du1flab/acAZvIs3bFZGMD36PQVmFAC97M3wzvjta25lpOxYR6TCNzOCdl5cHmUxW60sTwsPDYWVlpSiUACAwMBBisRgRERE1nldYWIhx48Zh/fr1cHR0fOx1fX19FYUSAAwePBgymQyxsbE1XrekpKRB2k1E/xjn54bdU/3hYGGIhMx8jPzyLI5cz9B2LCJqBlQqltq2bQtra+vHvqysrGBtba2RkI/r0dLX14eNjQ2kUmmN582ePRu9e/fGiBEjarzuvwslAIqfa7vuihUrYGlpqXi5uroq2xQiqofu7tbY/0ZfPNXKGnkl5Xh1+0V8HnYTcjnHMRGR5uircvCePXtgY2OjtjefP38+Vq5cWesxcXFxdbr277//jmPHjiEqKqpO59dmwYIFmDNnjuJnmUzGgomogdibG2Hnq73wvz+uY3v4HXxx9BZiUnPxeXAXWBgZaDseETVBKhVLffr0UeuYpblz52LSpEm1HuPp6QlHR0dkZlZ9bLi8vBzZ2dmPvb0GAMeOHUNiYiKsrKyqbB89ejT69euHEydOwNHREZGRkVX2Z2RUduvXdF0AMDQ0hKGhYa25iUhzJPpiLB3REb4ulnj/1xgcvZGJEV+excbx3dHGwVzb8YioiVGpWFI3Ozs72NnZPfE4f39/5OTk4NKlS+jevTuAymJILpfDz8/vsefMnz8fr776apVtvr6++Pzzz/H8888rrvvRRx8hMzNTUQSGhYXBwsICPj4+9WkaETWAMT1c0c7RHK99fwm3swowcv1ZfDqmM4b6Omk7GhE1IUqPWXJ3d6+yvElDat++PYYMGYLQ0FBERkbi7NmzmDlzJoKDgxVPwqWmpsLb21vRU+To6IiOHTtWeQGAm5sbPDw8AACDBg2Cj48Pxo8fjytXruDw4cNYuHAhZsyYwZ4jIh3RqaUV9r/RF708bVBQWoHpOy9j5aEbKK+QazsaETURShdLt2/fhq2trSaz1Grnzp3w9vZGQEAAgoKC0LdvX2zcuFGxv6ysDPHx8SgsLFT6mnp6ejhw4AD09PTg7++Pl19+GRMmTMDSpUs10QQi0hBbM0PsmOKHV/tW/ofQ1ycSMWL9WcSk5mo5GRE1BUrPs0Q14zxLRI3H/itpWPhrDHKLyqAnFuHVvh54K7AtjCXa6RknosZLI/MsERE1ds93dsaROQPwXCcnVMgFfHMqCYPXnMLZhCxtRyMiHcViiYiaHDtzQ3w5rhs2TegBJ0sjJGcXImRTBOb9fAU5haXajkdEOkapYsnGxgZZWZX/VfZoQV0iosYu0McBf83ujwn+7hCJgJ8vpSBw9Sn8cTWdC/ISkdKUKpZKS0sVS3ps27YNxcXFGg1FRKQu5kYGWDqiI36e5g8vezNk5Zdgxq7LCN1+Cem5RdqOR0Q6QKkB3s8++ywyMjLQvXt3bNu2DS+99BKMjY0fe+zmzZvVHrKx4wBvIt1QUl6Br44n4qsTCSirEGBmqI93h3ojpKcbxGKRtuMRUQNT6wDvHTt2ICgoCPn5+RCJRMjNzcWDBw8e+yIiaqwM9fUw+9m2+GNWP3R1s0J+STkW/RqDF78JR0JmvrbjEVEjpfLUAR4eHrh48aJW51xqbNizRKR7KuQCdpy/g08O3UBBaQUkemLMfMYLrw1oDYk+n30hag6U/f7mPEtqwGKJSHel5hRh4b5rOB5/DwDQ1sEMH4/uhG5u1lpORkSaptF5lk6ePInnn38eXl5e8PLywvDhw3H69Ok6hyUi0hYXK2NsnvQUvgjuAhtTCW5m5GP01+ewZH8sCkrKtR2PiBoBlYulHTt2IDAwECYmJpg1axZmzZoFY2NjBAQEYNeuXZrISESkUSKRCCO6uODInAEY1c0FggBsOfs3Bn1+CsfjM7Udj4i0TOXbcO3bt8fUqVMxe/bsKttXr16Nb7/9FnFxcWoNqAt4G46oaTl18x7e23cNKQ8qpxYY2cUZi57zga0ZF9gmako0NmbJ0NAQsbGx8PLyqrI9ISEBHTt2bJZzMLFYImp6CkvLsfqvm9h89jbkAmBtYoDFz/tgZBcXiEScZoCoKdDYmCVXV1ccPXq02vYjR47A1dVV1csRETVKJhJ9LHzOB/te7wNvR3M8KCzD7B+vYOKWC7ibXajteETUgPRVPWHu3LmYNWsWoqOj0bt3bwDA2bNnsXXrVnzxxRdqD0hEpE2dXa2w/42+2HgqCV8cvYVTN+9h0Oen8PbgdpjUuxX0OJklUZNXp6kD9u3bh88++0wxPql9+/aYN28eRowYofaAuoC34Yiah8R7+Viw9xoib2cDqCykVo72hbcj/39PpIs4z1IDYrFE1HzI5QJ2X7iLFX/GIa+kHPpiEaY/3RozBnrByEBP2/GISAUanWeJiKi5EotFGOfnhiNzB2CQjwPK5QLWHUtA0NrTih4nImpaWCwREdWBg4URNk7ogQ0vd4OduSGS7hXgxW/C8f6+a5AVl2k7HhGpEYslIqJ6GNLRCUfmDMDYnpVPA++MSMag1afwV6xUy8mISF1YLBER1ZOlsQFWjOqEXaF+aGVrAqmsGFO/v4TXd15CZl7zm3uOqKmpd7FUUVGB6OhoPHjwQB15iIh0Vu/WLXDorf6Y/nRr6IlF+POaFIGfncRPF+6Cz9IQ6S6Vi6W33noL3333HYDKQmnAgAHo1q0bXF1dceLECXXnIyLSKUYGenh3iDd+n9kHvi6WkBWX451friJkUwT+zirQdjwiqgOVi6U9e/agc+fOAID9+/fj9u3buHHjBmbPno33339f7QGJiHRRB2dL7Hu9N94Pag8jAzHOJd7H4DWnsOFkIsor5NqOR0QqULlYysrKgqOjIwDgzz//xJgxY9C2bVtMnjwZ165dU3tAIiJdpa8nRmh/T/z11gD09WqBknI5Pj54AyPWn0VMaq624xGRklQulhwcHHD9+nVUVFTg0KFDePbZZwEAhYWF0NPjhGxERP/lZmuC76f0xKdjOsPS2ACxaTKMWH8WK/6MQ1FphbbjEdETqFwsvfLKK3jxxRfRsWNHiEQiBAYGAgAiIiLg7e2t9oBERE2BSCTC/3VviSNzBuD5zs6okAv45lQShnxxCucSsrQdj4hqUaflTvbs2YO7d+9izJgxaNmyJQBg27ZtsLKyapbrw3G5EyJS1dG4DCz8NQbpuZVTC7zYoyXeD/KBpYmBlpMRNR8NsjZccXExjIyM6np6k8FiiYjqIq+4DJ8ejsf283cgCEALM0MsGd4BQb6OEIlE2o5H1ORpbG24iooKLFu2DC4uLjAzM0NSUhIAYNGiRYopBYiI6MnMjQywZERH7HnNH172ZsjKL8GMXZcRuv0S0nOLtB2PiB5SuVj66KOPsHXrVnzyySeQSCSK7R07dsSmTZvUGo6IqDno7m6DP2b1xZsBbWCgJ8KRuAw8u/oUvj9/B3I5J7Mk0jaVi6Xt27dj48aNCAkJqfL0W+fOnXHjxg21hiMiai4M9fUw+9m2+GNWP3R1s0J+STkW/RqDF78JR0JmvrbjETVrKhdLqamp8PLyqrZdLpejrIwrbRMR1UdbB3Psea03lgzvAFOJHi7eeYCgL05j3dFbKC3nZJZE2qByseTj44PTp09X275nzx507dpVLaGIiJozPbEIE3u3wl9zBmBgOzuUVsjxWdhNPL/uDKKSuQ4nUUPTV/WExYsXY+LEiUhNTYVcLsfevXsRHx+P7du348CBA5rISETULLlYGWPzpKfw+5U0LN1/HfEZeRj19TkM7+yM0H6e6Ohiqe2IRM2Cyj1LI0aMwP79+3HkyBGYmppi8eLFiIuLw/79+xWzeWtCdnY2QkJCYGFhASsrK0yZMgX5+crdxxcEAUOHDoVIJMKvv/5aZZ9IJKr22r17twZaQESkOpFIhBFdXHBkzgCM7tYSggD8Fp2G59adQcim8zh58x7qMQMMESlB5Z4lAOjXrx/CwsLUnaVWISEhSE9PR1hYGMrKyvDKK69g6tSp2LVr1xPPXbNmTa1zlmzZsgVDhgxR/GxlZaWOyEREamNtKsFnL3bGK31a4dvTSThwNR1nE+7jbMJ9eDuaI7SfJ57v7AyJvsr/DUxET1CvSSkbSlxcHHx8fHDhwgX06NEDAHDo0CEEBQUhJSUFzs7ONZ4bHR2N5557DhcvXoSTkxP27duHkSNHKvaLRKJq21TFSSmJqKGl5hRh85nb2B2ZjIKH68s5WhjhlT6tMNbPDRZGnAmc6EnUOimljY0NsrIq1y6ytraGjY1NjS9NCA8Ph5WVlaJQAoDAwECIxWJERETUeF5hYSHGjRuH9evXw9HRscbjZsyYgRYtWqBnz57YvHnzE7u0S0pKIJPJqryIiBqSi5UxFj3ng3MLAvDuEG/YmxtCKivGioM30HvFMXz0x3Wk5XBiSyJ1UOo23Oeffw5zc3PFnxt6Gn6pVAp7e/sq2/T19WFjYwOpVFrjebNnz0bv3r1rXa9u6dKleOaZZ2BiYoK//voLr7/+OvLz8zFr1qwaz1mxYgWWLFmiekOIiNTM0tgA059ujcl9W+H36DR8ezoJNzPy8e3p29hy9m88/3AwuI8ze72J6kqrt+Hmz5+PlStX1npMXFwc9u7di23btiE+Pr7KPnt7eyxZsgTTp0+vdt7vv/+OuXPnIioqCmZmZgCUu+W2ePFibNmyBXfv3q3xmJKSEpSUlCh+lslkcHV15W04ItI6QRBw4uY9bDyZhPCk+4rt/dq0wNT+nujr1YLrzhE9pOxtOJUHeOvp6SE9Pb1aT8/9+/dhb2+PiooKpa81d+5cTJo0qdZjPD094ejoiMzMzCrby8vLkZ2dXePttWPHjiExMbHaYO3Ro0ejX79+OHHixGPP8/Pzw7Jly1BSUgJDQ8PHHmNoaFjjPiIibRKJRBjYzh4D29njWkouNp5Owp/X0nH6VhZO38pCeycLTO3vgec6OcNAj4PBiZShcs+SWCx+7G2xtLQ0tG7dGkVF6r9H/miA98WLF9G9e3cAwF9//YUhQ4bUOMBbKpUqxlk94uvriy+++ALPP/88PDw8HvteH330ET777DNkZ2crnY8DvImoMbubXYjNZ2/jxwt3UfhwMLiTpREm9/FAcE9XmHMwODVTau9ZWrt2LYDK/2rZtGmT4tYWAFRUVODUqVPw9vauR+SatW/fHkOGDEFoaCg2bNiAsrIyzJw5E8HBwYpCKTU1FQEBAdi+fTt69uwJR0fHx/Y6ubm5KQql/fv3IyMjA7169YKRkRHCwsKwfPlyvP322xppBxGRNrjamOCD5zvgzYA22BmRjK3n/kZ6bjE++jMOa4/ewjg/N0zq0wpOlsbajkrUKCldLH3++ecAKu+Hb9iwocoiuhKJBK1atcKGDRvUn/ChnTt3YubMmQgICIBYLMbo0aMVBRwAlJWVIT4+HoWFhUpf08DAAOvXr8fs2bMhCAK8vLywevVqhIaGaqIJRERaZWUiwYyBXni1nwd+i0rDxtNJSMjMxzenkvDdmdsY3qVyMHh7J/aQE/2byrfhBg4ciL1798La2lpTmXQOb8MRkS6SywWcuJmJb04mIeL2P0MP+re1w7T+nujd2paDwalJU/b7WycmpWzsWCwRka67cjcHG08n4eC1dMgffit0cLbA1P6eCPJ14mBwapLUWizNmTMHy5Ytg6mpKebMmVPrsatXr1Y9rY5jsURETcXd7EJ8d6ZyMHhRWeVgcBcrY7zSpxWCe7rBzLBOq2QRNUpqLZYGDhyIffv2wcrKCgMHDqz5YiIRjh07VrfEOozFEhE1NQ8KSrEz4g62nvsbWfmlAABzI32E+LnjlT6t4GBhpOWERPXH23ANiMUSETVVxWUV+DUqFRtPJyHpXgEAwEBPhBFdXDC1vyfaOphrOSFR3WmsWNqxYwdGjRoFExOTeodsKlgsEVFTJ5cLOHYjExtPJSHy738Ggz/dzg5T+3vC35ODwUn3aKxYsrOzQ1FREYYPH46XX34ZgwcPrjKNQHPEYomImpOo5Af49nQSDsVIFYPBO7pYYGr/1gjq6Ah9DgYnHaGxYqm8vByHDh3CDz/8gN9++w0mJiYYM2YMQkJC0Lt373oH10UsloioObpzvwDfnbmNny7eRXGZHEDlYPApfT3w0lOuMOVgcGrkGmTMUmFhIfbt24ddu3bhyJEjaNmyJRITE+t6OZ3FYomImrMHBaX4/vwdbDv3N+4XVA4GtzDSx8u93DGpdyvYczA4NVINNsA7KysLu3fvxoYNGxAXF6fSQrpNBYslIqLKweB7L6di0+kkJGVVDgaX6IkxsmvlzOBtOBicGhmNFkuPepR27tyJo0ePwtXVFWPHjkVISIjG1odrzFgsERH9Qy4XcCQuAxtPJeHinQeK7c9422Nqf0/4edhwMDg1ChorloKDg3HgwAGYmJjgxRdfREhICPz9/esdWJexWCIierxLdx7g21NJOHxdikffNp1bWiK0vyeGdOBgcNIuZb+/VR59p6enh59++olPwRER0RN1d7dG9/HdcTurAN+dScLPF1NwJSUXM3dFwdXGGFP6eODFp1xhIuFgcGq8OCmlGrBniYhIOffzS/D9+TvYHn4H2Q8Hg1saG2B8L3dM7N0KduaGWk5IzYlab8OtXbsWU6dOhZGREdauXVvrsbNmzVI9rY5jsUREpJqi0gr8cjkFm04n4e/7hQAAib4Ygzs44lkfBwxoawdLYwMtp6SmTq3FkoeHBy5evAhbW1t4eHjUfDGRCElJSXVLrMNYLBER1U2FXEDY9QxsPJWIy8k5iu36YhF6etggsL0DAts7wM2Wq0aQ+nFtuAbEYomIqP6i7+bgcKwUR65n4FZmfpV97RzMEdDeHoE+DujS0gpiMZ+mo/rTWLG0dOlSvP3229XWhisqKsKqVauwePHiuiXWYSyWiIjU6++sAhyJy8CRuAxc+PsBKuT/fFW1MDNEgLc9Atrbo2+bFhwcTnWmsWJJT08P6enpsLe3r7L9/v37sLe356SULJaIiNQqt7AMJ25mIux6Bk7G30NeSblin6G+GH29WiDQxwEB3vacLZxUorGpAwRBeOxkYleuXIGNjY2qlyMiIqqVpYkBRnRxwYguLigtl+PC39kIu17Z65TyoAhHb2Ti6I1MAJVzOAW2d0CgjwO8Hc05+SWphdI9S9bW1hCJRIrq69+/gBUVFcjPz8drr72G9evXayxsY8WeJSKihicIAuIz8nA0rrLXKfpuTpX9LlbGCHw4zsnPwxYSfU6ASVWp/Tbctm3bIAgCJk+ejDVr1sDS0lKxTyKRoFWrVs12Jm8WS0RE2peZV4xjcZk4EpeJMwn3UFwmV+wzM9THgLZ2CPSxx8B29rAykWgxKTUWGhuzdPLkSfTu3RsGBpz/4hEWS0REjUtRaQXOJmThSFwGjt7IxL28EsU+PbEIPdytFbfrPFqYajEpaZNaiyWZTKa4iEwmq/XY5lgssFgiImq85HIBV1NzceThOKcb0rwq+1vbmSLQp3I+p25u1tDjtATNhlqLpX8/AScWix87YO7RwG8+DcdiiYioMbubXYijcRk4EpeJ80n3Uf6vaQlsTCUY2M4ez/rYo18bO5gaclqCpkytxdLJkyfRp08f6Ovr4+TJk7UeO2DAANXT6jgWS0REuklWXIZTN+/hyPUMHLuRCVnxP9MSSPTE8G9t+7DXyR5OlsZaTEqawBm8GxCLJSIi3VdWIcfFvx/gaFwGwuIycOfhmnWPdHC2QGB7Bzzr44AOzhaclqAJ0FixdOjQIZiZmaFv374AgPXr1+Pbb7+Fj48P1q9fD2tr6/ol10EsloiImhZBEJB4Lx9h1zNxJC4Dl5Mf4N/flo4WRgj0sUdAewf4e9rCyEBPe2GpzjRWLPn6+mLlypUICgrCtWvX0KNHD8ydOxfHjx+Ht7c3tmzZUu/wuobFEhFR05aVX4LjNyoLp1M3s1BU9s/4XBOJHvq3sUOgjwMGtrODrZmhFpOSKjRWLJmZmSEmJgatWrXChx9+iJiYGOzZsweXL19GUFAQpFJpvcPrGhZLRETNR3FZBcKT7iuersuQ/TMtgUgEdHezVjxd19rOlLfrGjGNLXcikUhQWFh5H/fIkSOYMGECAMDGxuaJ0woQERHpOiMDPQxsVzm55f9GdkRMqkyx6G9smgwX7zzAxTsP8PHBG2hla6KYz6mHuzX09TiLuC5SuWdp+PDhKC0tRZ8+fbBs2TLcvn0bLi4u+OuvvzBz5kzcvHlTU1kbLfYsERERAKTlFCmmJQhPvI/Sin9mETeV6KGDiyU6uVjCt6UlOrW0gruNCcSc10lrNHYbLjk5Ga+//jru3r2LWbNmYcqUKQCA2bNno6KiAmvXrq1fch3EYomIiP4rv6Qcp2/eQ1hcBo7fyMSDwrJqx5gb6cPXxbLy1dISnVys4GpjzFt3DYRTBzQgFktERFSbCnnl03VXU3JxLSUHV1NzcT1NhpJyebVjLY0N0KllZQHVqaUlfFtawdnSiAWUBmi0WJLL5UhISEBmZibk8qofdP/+/VVPq+NYLBERkarKKuS4lZGPa6k5lUVUai5upOdVuXX3iK2pBB0fFU8ulbfwHCwMWUDVk8aKpfPnz2PcuHG4c+cO/nsqlzthsURERHVXWi7HzYy8h8VTZREVL82rsiTLI3bmhv8a/2SJji6WsDc30kJq3aWxYqlLly5o27YtlixZAicnp2pVraWlZd0SP0F2djbeeOMN7N+/H2KxGKNHj8YXX3wBMzOzGs95+umnqy3PMm3aNGzYsEHxc3JyMqZPn47jx4/DzMwMEydOxIoVK6Cvr/yDgiyWiIhIU4rLKnBDmld5++5hD9TNjDw8pn6Co4XRw7FPlUWUr4sl532qhcamDrh16xb27NkDLy+vegVUVUhICNLT0xEWFoaysjK88sormDp1Knbt2lXreaGhoVi6dKniZxMTE8WfKyoqMGzYMDg6OuLcuXNIT0/HhAkTYGBggOXLl2usLURERMoyMtBDF1crdHG1UmwrKq3A9fRcXEvJxdXUyv9NuJcPqawY0uvFCLueoTjWxcr44dgnS8VgcisTiRZaortU7ll65pln8M4772DIkCGaylRNXFwcfHx8cOHCBfTo0QNA5bIrQUFBSElJgbOz82PPe/rpp9GlSxesWbPmsfsPHjyI5557DmlpaXBwcAAAbNiwAe+++y7u3bsHiUS5Xyb2LBERkbYVlJQjNk2Gqyk5uPawgErKKnjssW42JlV6oDq6WMLCyKCBE2ufxnqW3njjDcydOxdSqRS+vr4wMKj6l9upUyfV0z5BeHg4rKysFIUSAAQGBkIsFiMiIgIvvPBCjefu3LkTO3bsgKOjI55//nksWrRI0bsUHh4OX19fRaEEAIMHD8b06dMRGxuLrl27PvaaJSUlKCn5Z8ZWTsZJRETaZmqoj54eNujpYaPYJisuQ2yqrMog8jv3C5GcXfn642q64ljPFqZVBpF3cLGEmaHKZUKTpPLfwujRowEAkydPVmwTiUQQBEFjA7ylUins7e2rbNPX14eNjU2ty6uMGzcO7u7ucHZ2xtWrV/Huu+8iPj4ee/fuVVz334USAMXPtV13xYoVWLJkSV2bQ0RE1CAsjAzg39oW/q1tFdtyCksRkyrD1dScytt4KblIzSlCUlYBkrIK8PuVNACVS7e0tjOrMojcx8kSxpLmt2iwysXS7du31fbm8+fPx8qVK2s9Ji4urs7Xnzp1quLPvr6+cHJyQkBAABITE9G6des6X3fBggWYM2eO4meZTAZXV9c6X4+IiKihWJlI0LdNC/Rt00Kx7X5+Ca6l5iImNVfRA5WeW4yEzHwkZOZjb1QqAEAsAtrYmyuKJ18XS7R3soCRQdMuoFQultzd3dX25nPnzsWkSZNqPcbT0xOOjo7IzMyssr28vBzZ2dlwdHRU+v38/PwAAAkJCWjdujUcHR0RGRlZ5ZiMjMpBcbVd19DQEIaGfLqAiIiaBlszQzzdzh5Pt/vnLk5mXvE/xVNKLq6k5CIrvwTxGXmIz8jDnkspAAB9sQhtHMzhZW+G1namaG1nhtZ2ZvC0M20yRVSdbkZ+//332LBhA27fvo3w8HC4u7tjzZo18PDwwIgRI5S+jp2dHezs7J54nL+/P3JycnDp0iV0794dAHDs2DHI5XJFAaSM6OhoAICTk5Piuh999BEyMzMVt/nCwsJgYWEBHx8fpa9LRETU1NibG+EZbyM84105PEUQBGTISnA1JaeyiHpYSGUXlCIuXYa49Krjd0WiyifxHhVOj4qo1vamsDPTrQk1VX4a7uuvv8bixYvx1ltv4aOPPkJMTAw8PT2xdetWbNu2DcePH9dI0KFDhyIjIwMbNmxQTB3Qo0cPxdQBqampCAgIwPbt29GzZ08kJiZi165dCAoKgq2tLa5evYrZs2ejZcuWirmXKioq0KVLFzg7O+OTTz6BVCrF+PHj8eqrr6o0dQCfhiMiouZIEASk5RYjNrXyybvEzHwk3stH4r0C5BZVXwvvEXMj/WpFlJe9KdxsTCHRFzdYfo1NSunj44Ply5dj5MiRMDc3x5UrV+Dp6YmYmBg8/fTTyMrKqnf4x8nOzsbMmTOrTEq5du1axaSUf//9Nzw8PHD8+HE8/fTTuHv3Ll5++WXExMSgoKAArq6ueOGFF7Bw4cIqfyF37tzB9OnTceLECZiammLixIn4+OOPOSklERFRHQmCgPsFpUi6V1BZPD0sopKyCnA3u/CxE2oCgJ5YBHcbk2o9Ua3tzDQyN5TGiiVjY2PcuHED7u7uVYqlW7duoVOnTigqKqp3eF3DYomIiEg5xWUVuHO/sFoRlZiZj4LSmp+o3/t6b3Rzs1ZrFo3Ns+Th4YHo6OhqA70PHTqE9u3bq56UiIiImg0jAz20czRHO0fzKtsfjYmqvI2XX6VXKi23GB62plpKXIdiac6cOZgxYwaKi4shCAIiIyPxww8/YMWKFdi0aZMmMhIREVETJxKJ4GhpBEdLI/TxalFlX2FpOUwk2psgU+V3fvXVV2FsbIyFCxeisLAQ48aNg7OzM7744gsEBwdrIiMRERE1Y9oslIA6jFn6t8LCQuTn51ebXbu54ZglIiIi3aPs97fKz+cVFRWhsLAQAGBiYoKioiKsWbMGf/31V93TEhERETVSKhdLI0aMwPbt2wEAOTk56NmzJz777DOMGDECX3/9tdoDEhEREWmTysXS5cuX0a9fPwDAnj174OjoiDt37mD79u1Yu3at2gMSERERaZPKxVJhYSHMzSsf9/vrr78watQoiMVi9OrVC3fu3FF7QCIiIiJtUrlY8vLywq+//oq7d+/i8OHDGDRoEAAgMzOTg5uJiIioyVG5WFq8eDHefvtttGrVCn5+fvD39wdQ2cvUtWtXtQckIiIi0qY6TR0glUqRnp6Ozp07QyyurLciIyNhYWEBb29vtYds7Dh1ABERke7R2HInAODo6AhHR8cq23r27FmXSxERERE1airfhiMiIiJqTlgsEREREdWCxRIRERFRLVgsEREREdVCu8v4NhGPHiiUyWRaTkJERETKevS9/aSJAVgsqUFeXh4AwNXVVctJiIiISFV5eXmwtLSscX+d5lmiquRyOdLS0mBubg6RSKS268pkMri6uuLu3bvNbv6m5tr25tpugG1vjm1vru0G2PbG0nZBEJCXlwdnZ2fFvJGPw54lNRCLxWjZsqXGrm9hYaH1Xyhtaa5tb67tBtj25tj25tpugG1vDG2vrUfpEQ7wJiIiIqoFiyUiIiKiWrBYasQMDQ3xwQcfwNDQUNtRGlxzbXtzbTfAtjfHtjfXdgNsu661nQO8iYiIiGrBniUiIiKiWrBYIiIiIqoFiyUiIiKiWrBYIiIiIqoFi6VGbP369WjVqhWMjIzg5+eHyMhIbUeqlw8//BAikajKy9vbW7G/uLgYM2bMgK2tLczMzDB69GhkZGRUuUZycjKGDRsGExMT2NvbY968eSgvL2/optTq1KlTeP755+Hs7AyRSIRff/21yn5BELB48WI4OTnB2NgYgYGBuHXrVpVjsrOzERISAgsLC1hZWWHKlCnIz8+vcszVq1fRr18/GBkZwdXVFZ988ommm/ZET2r7pEmTqv0ODBkypMoxutj2FStW4KmnnoK5uTns7e0xcuRIxMfHVzlGXb/fJ06cQLdu3WBoaAgvLy9s3bpV082rlTJtf/rpp6t97q+99lqVY3Sx7V9//TU6deqkmFzR398fBw8eVOxvqp/5k9rdJD9vgRql3bt3CxKJRNi8ebMQGxsrhIaGClZWVkJGRoa2o9XZBx98IHTo0EFIT09XvO7du6fY/9prrwmurq7C0aNHhYsXLwq9evUSevfurdhfXl4udOzYUQgMDBSioqKEP//8U2jRooWwYMECbTSnRn/++afw/vvvC3v37hUACPv27auy/+OPPxYsLS2FX3/9Vbhy5YowfPhwwcPDQygqKlIcM2TIEKFz587C+fPnhdOnTwteXl7C2LFjFftzc3MFBwcHISQkRIiJiRF++OEHwdjYWPjmm28aqpmP9aS2T5w4URgyZEiV34Hs7Owqx+hi2wcPHixs2bJFiImJEaKjo4WgoCDBzc1NyM/PVxyjjt/vpKQkwcTERJgzZ45w/fp1Yd26dYKenp5w6NChBm3vvynT9gEDBgihoaFVPvfc3FzFfl1t+++//y788ccfws2bN4X4+HjhvffeEwwMDISYmBhBEJruZ/6kdjfFz5vFUiPVs2dPYcaMGYqfKyoqBGdnZ2HFihVaTFU/H3zwgdC5c+fH7svJyREMDAyEn3/+WbEtLi5OACCEh4cLglD5RSwWiwWpVKo45uuvvxYsLCyEkpISjWavq/8WDHK5XHB0dBRWrVql2JaTkyMYGhoKP/zwgyAIgnD9+nUBgHDhwgXFMQcPHhREIpGQmpoqCIIgfPXVV4K1tXWVdr/77rtCu3btNNwi5dVULI0YMaLGc5pK2zMzMwUAwsmTJwVBUN/v9zvvvCN06NChynu99NJLwuDBgzXdJKX9t+2CUPnl+eabb9Z4TlNpuyAIgrW1tbBp06Zm9ZkLwj/tFoSm+XnzNlwjVFpaikuXLiEwMFCxTSwWIzAwEOHh4VpMVn+3bt2Cs7MzPD09ERISguTkZADApUuXUFZWVqXN3t7ecHNzU7Q5PDwcvr6+cHBwUBwzePBgyGQyxMbGNmxD6uj27duQSqVV2mlpaQk/P78q7bSyskKPHj0UxwQGBkIsFiMiIkJxTP/+/SGRSBTHDB48GPHx8Xjw4EEDtaZuTpw4AXt7e7Rr1w7Tp0/H/fv3FfuaSttzc3MBADY2NgDU9/sdHh5e5RqPjmlM/y78t+2P7Ny5Ey1atEDHjh2xYMECFBYWKvY1hbZXVFRg9+7dKCgogL+/f7P5zP/b7kea2ufNhXQboaysLFRUVFT5RQIABwcH3LhxQ0up6s/Pzw9bt25Fu3btkJ6ejiVLlqBfv36IiYmBVCqFRCKBlZVVlXMcHBwglUoBAFKp9LF/J4/26YJHOR/Xjn+3097evsp+fX192NjYVDnGw8Oj2jUe7bO2ttZI/voaMmQIRo0aBQ8PDyQmJuK9997D0KFDER4eDj09vSbRdrlcjrfeegt9+vRBx44dFbnU8ftd0zEymQxFRUUwNjbWRJOU9ri2A8C4cePg7u4OZ2dnXL16Fe+++y7i4+Oxd+9eALrd9mvXrsHf3x/FxcUwMzPDvn374OPjg+jo6Cb9mdfUbqBpft4slqjBDB06VPHnTp06wc/PD+7u7vjpp5+0/o88NYzg4GDFn319fdGpUye0bt0aJ06cQEBAgBaTqc+MGTMQExODM2fOaDtKg6up7VOnTlX82dfXF05OTggICEBiYiJat27d0DHVql27doiOjkZubi727NmDiRMn4uTJk9qOpXE1tdvHx6dJft68DdcItWjRAnp6etWemsjIyICjo6OWUqmflZUV2rZti4SEBDg6OqK0tBQ5OTlVjvl3mx0dHR/7d/Jony54lLO2z9bR0RGZmZlV9peXlyM7O7tJ/V0AgKenJ1q0aIGEhAQAut/2mTNn4sCBAzh+/Dhatmyp2K6u3++ajrGwsND6f3DU1PbH8fPzA4Aqn7uutl0ikcDLywvdu3fHihUr0LlzZ3zxxRdN/jOvqd2P0xQ+bxZLjZBEIkH37t1x9OhRxTa5XI6jR49WuSes6/Lz85GYmAgnJyd0794dBgYGVdocHx+P5ORkRZv9/f1x7dq1Kl+mYWFhsLCwUHT/NnYeHh5wdHSs0k6ZTIaIiIgq7czJycGlS5cUxxw7dgxyuVzxj46/vz9OnTqFsrIyxTFhYWFo166d1m9DqSIlJQX379+Hk5MTAN1tuyAImDlzJvbt24djx45Vu02ort9vf3//Ktd4dIw2/114UtsfJzo6GgCqfO662PbHkcvlKCkpadKf+eM8avfjNInPWyvDyumJdu/eLRgaGgpbt24Vrl+/LkydOlWwsrKq8vSArpk7d65w4sQJ4fbt28LZs2eFwMBAoUWLFkJmZqYgCJWP2bq5uQnHjh0TLl68KPj7+wv+/v6K8x89bjpo0CAhOjpaOHTokGBnZ9fopg7Iy8sToqKihKioKAGAsHr1aiEqKkq4c+eOIAiVUwdYWVkJv/32m3D16lVhxIgRj506oGvXrkJERIRw5swZoU2bNlUen8/JyREcHByE8ePHCzExMcLu3bsFExMTrU8dUFvb8/LyhLffflsIDw8Xbt++LRw5ckTo1q2b0KZNG6G4uFhxDV1s+/Tp0wVLS0vhxIkTVR6XLiwsVByjjt/vR49Tz5s3T4iLixPWr1+v9cfIn9T2hIQEYenSpcLFixeF27dvC7/99pvg6ekp9O/fX3ENXW37/PnzhZMnTwq3b98Wrl69KsyfP18QiUTCX3/9JQhC0/3Ma2t3U/28WSw1YuvWrRPc3NwEiUQi9OzZUzh//ry2I9XLSy+9JDg5OQkSiURwcXERXnrpJSEhIUGxv6ioSHj99dcFa2trwcTERHjhhReE9PT0Ktf4+++/haFDhwrGxsZCixYthLlz5wplZWUN3ZRaHT9+XABQ7TVx4kRBECqnD1i0aJHg4OAgGBoaCgEBAUJ8fHyVa9y/f18YO3asYGZmJlhYWAivvPKKkJeXV+WYK1euCH379hUMDQ0FFxcX4eOPP26oJtaotrYXFhYKgwYNEuzs7AQDAwPB3d1dCA0NrfYfALrY9se1GYCwZcsWxTHq+v0+fvy40KVLF0EikQienp5V3kMbntT25ORkoX///oKNjY1gaGgoeHl5CfPmzasy744g6GbbJ0+eLLi7uwsSiUSws7MTAgICFIWSIDTdz7y2djfVz1skCILQcP1YRERERLqFY5aIiIiIasFiiYiIiKgWLJaIiIiIasFiiYiIiKgWLJaIiIiIasFiiYiIiKgWLJaIiIiIasFiiYiIiKgWLJaIiIiIasFiiYioFrNnz8aoUaO0HYOItIjFEhFRLSIjI9GjRw9txyAiLeLacEREj1FaWgpTU1OUl5crtvn5+eH8+fNaTEVE2qCv7QBERI2Rvr4+zp49Cz8/P0RHR8PBwQFGRkbajkVEWsBiiYjoMcRiMdLS0mBra4vOnTtrOw4RaRHHLBER1SAqKoqFEhGxWCIiqkl0dDSLJSJisUREVJNr166hS5cu2o5BRFrGYomIqAZyuRzx8fFIS0tDbm6utuMQkZawWCIiqsH//vc/bN26FS4uLvjf//6n7ThEpCWcZ4mIiIioFuxZIiIiIqoFiyUiIiKiWrBYIiIiIqoFiyUiIiKiWrBYIiIiIqoFiyUiIiKiWrBYIiIiIqoFiyUiIiKiWrBYIiIiIqoFiyUiIiKiWrBYIiIiIqrF/wNHTA6OTX7COQAAAABJRU5ErkJggg==",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "solution = sim.solve(\n",
+    "    [0, 3600], inputs={\"Current function [A]\": 0.15652}, calculate_sensitivities=True\n",
+    ")\n",
+    "plt.plot(\n",
+    "    solution.t, solution[\"Terminal voltage [V]\"].sensitivities[\"Current function [A]\"]\n",
+    ")\n",
+    "\n",
+    "plt.xlabel(r\"$t$\")\n",
+    "plt.ylabel(\"sensitivities of Terminal voltage wrt Current fuction\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d84c83fb",
+   "metadata": {},
+   "source": [
+    "## Sensitivities and data fitting\n",
+    "\n",
+    "Sensitivities are often used to aid data fitting by providing a means to calculate the gradient of the function to be minimised. Take for example the data fitting exercise we introduced in the previous notebook. Once again we will generate some fake data for fitting, like so:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "33d95c39",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "t_eval = np.linspace(0, 3600, 100)\n",
+    "data = sim.solve([0, 3600], inputs={\"Current function [A]\": 0.2222})[\n",
+    "    \"Terminal voltage [V]\"\n",
+    "](t_eval)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d5045b57",
+   "metadata": {},
+   "source": [
+    "Now we will contruct a function to minimise, but here we will return both the value of the function, and its gradient with respect to the input parameter \"Current function\". Note that our objective function is the sum of squared different between the vector $\\mathbf{f}$, the simulated \"Terminal voltage\", and $\\mathbf{d}$, the vector of fake data, given by\n",
+    "\n",
+    "$$\\mathcal{O}(a) = \\sum_{i=0}^{i=N} (f_i(a) - d_i)^2$$ \n",
+    "\n",
+    "where $a$ is the parameter to be optimised (in this case \"Current function\"), $f_i$ is each element of the vector $\\mathbf{f}$, and $d_i$ is each element of $\\mathbf{d}$. We wish to also find the gradient of this function wrt the parameter $a$, which is:\n",
+    "\n",
+    "$$\\frac{\\partial \\mathcal{O}}{\\partial a}(a) = 2 \\sum_{i=0}^{i=N} (f_i(a) - d_i) \\frac{\\partial f_i}{\\partial a} $$ \n",
+    "\n",
+    "Using these equations, we will define a function that takes in as an argument $a$, and returns $(\\mathcal{O}(a), \\frac{\\partial \\mathcal{O}}{\\partial a}(a))$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "ffad2bc0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def sum_of_squares_jac(parameters):\n",
+    "    sol = sim.solve(\n",
+    "        [0, 3600],\n",
+    "        t_interp=t_eval,\n",
+    "        inputs={\"Current function [A]\": parameters[0]},\n",
+    "        calculate_sensitivities=True,\n",
+    "    )\n",
+    "    term_voltage = sol[\"Terminal voltage [V]\"].data\n",
+    "    term_voltage_sens = sol[\"Terminal voltage [V]\"].sensitivities[\n",
+    "        \"Current function [A]\"\n",
+    "    ]\n",
+    "\n",
+    "    f = np.sum((term_voltage - data) ** 2)\n",
+    "    g = 2 * np.sum((term_voltage - data) * term_voltage_sens)\n",
+    "    print(\n",
+    "        f\"evaluating function and jacobian for p = {parameters[0]}, \\tloss = {f}, grad = {g}\"\n",
+    "    )\n",
+    "    return f, g"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fdcae8ac",
+   "metadata": {},
+   "source": [
+    "We can then use this function along with an optimisation algorithm to recover the value of the Current function that was used to generate the data. In this case we will use the `scipy.optimize` module again. This module allows the use of a function in the form given above to perform the minimisation, using both the value of the objective function and its gradient to find the minimum value of $a$ in the least number of steps.\n",
+    "\n",
+    "Once again, we will place bounds on \"Current function [A]\" between $(0.01, 0.6)$, and use a random starting value $x_0$ between these bounds."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "44f52a7e",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "starting parameter is 0.4035076613514513\n",
+      "evaluating function and jacobian for p = 0.4035076613514513, \tloss = 0.358632833514908, grad = 3.7278265436340283\n",
+      "evaluating function and jacobian for p = 0.01, \tloss = 0.8765036816401419, grad = -9.691403058800336\n",
+      "evaluating function and jacobian for p = 0.23970725060294026, \tloss = 0.0036706826105448887, grad = 0.41303112361463107\n",
+      "evaluating function and jacobian for p = 0.21929734316448973, \tloss = 0.00010628748342624681, grad = -0.07293742235816741\n",
+      "evaluating function and jacobian for p = 0.22236059911277223, \tloss = 2.5627985467376454e-07, grad = 0.0035066000930420284\n",
+      "evaluating function and jacobian for p = 0.22222008304298907, \tloss = 7.758859239380127e-09, grad = 2.935984691530423e-05\n",
+      "evaluating function and jacobian for p = 0.22221889660489277, \tloss = 7.74141508243804e-09, grad = -1.183255640750795e-08\n",
+      "recovered parameter is 0.22221889660489277\n"
+     ]
+    }
+   ],
+   "source": [
+    "bounds = (0.01, 0.6)\n",
+    "x0 = np.random.uniform(low=bounds[0], high=bounds[1])\n",
+    "\n",
+    "print(\"starting parameter is\", x0)\n",
+    "res = scipy.optimize.minimize(sum_of_squares_jac, [x0], bounds=[bounds], jac=True)\n",
+    "print(\"recovered parameter is\", res.x[0])"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "env",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

From 9a479cf1da0fbd8cf026c6ca577b2369dc59efe2 Mon Sep 17 00:00:00 2001
From: "Eric G. Kratz" <kratman@users.noreply.github.com>
Date: Wed, 6 Nov 2024 12:19:25 -0500
Subject: [PATCH 12/13] Make ParameterValues.pop return a value (#4571)

* Fix pop

* Remove some text
---
 src/pybamm/parameters/parameter_values.py | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/src/pybamm/parameters/parameter_values.py b/src/pybamm/parameters/parameter_values.py
index 0a0e49cd8f..003f8beca2 100644
--- a/src/pybamm/parameters/parameter_values.py
+++ b/src/pybamm/parameters/parameter_values.py
@@ -1,6 +1,3 @@
-#
-# Parameter values for a simulation
-#
 import numpy as np
 import pybamm
 import numbers
@@ -184,7 +181,7 @@ def items(self):
         return self._dict_items.items()
 
     def pop(self, *args, **kwargs):
-        self._dict_items.pop(*args, **kwargs)
+        return self._dict_items.pop(*args, **kwargs)
 
     def copy(self):
         """Returns a copy of the parameter values. Makes sure to copy the internal

From f05cae2ecda531cde049ccfc69d470cdb845bf98 Mon Sep 17 00:00:00 2001
From: Martin Robinson <martinjrobins@gmail.com>
Date: Thu, 7 Nov 2024 16:28:03 +0000
Subject: [PATCH 13/13] bug: use direct casadi bspline function for 1D & 2D
 cubic interp (#4572)

* bug: use direct casadi bspline function for 1D cubic interp

* bug: use direct bspline func for 2d cubic interp
---
 .../operations/convert_to_casadi.py           | 27 +++++++++++++++++++
 1 file changed, 27 insertions(+)

diff --git a/src/pybamm/expression_tree/operations/convert_to_casadi.py b/src/pybamm/expression_tree/operations/convert_to_casadi.py
index af53b9acd8..6b61b35263 100644
--- a/src/pybamm/expression_tree/operations/convert_to_casadi.py
+++ b/src/pybamm/expression_tree/operations/convert_to_casadi.py
@@ -7,6 +7,7 @@
 import casadi
 import numpy as np
 from scipy import special
+from scipy import interpolate
 
 
 class CasadiConverter:
@@ -165,6 +166,18 @@ def _convert(self, symbol, t, y, y_dot, inputs):
                             # for some reason, pybamm.Interpolant always returns a column vector, so match that
                             test = test.T
                         return test
+                    elif solver == "bspline":
+                        bspline = interpolate.make_interp_spline(
+                            symbol.x[0], symbol.y, k=3
+                        )
+                        knots = [bspline.t]
+                        coeffs = bspline.c.flatten()
+                        degree = [bspline.k]
+                        m = len(coeffs) // len(symbol.x[0])
+                        f = casadi.Function.bspline(
+                            symbol.name, knots, coeffs, degree, m
+                        )
+                        return f(converted_children[0])
                     else:
                         return casadi.interpolant(
                             "LUT", solver, symbol.x, symbol.y.flatten()
@@ -176,6 +189,20 @@ def _convert(self, symbol, t, y, y_dot, inputs):
                             symbol.y.ravel(order="F"),
                             converted_children,
                         )
+                    elif solver == "bspline" and len(converted_children) == 2:
+                        bspline = interpolate.RectBivariateSpline(
+                            symbol.x[0], symbol.x[1], symbol.y
+                        )
+                        [tx, ty, c] = bspline.tck
+                        [kx, ky] = bspline.degrees
+                        knots = [tx, ty]
+                        coeffs = c
+                        degree = [kx, ky]
+                        m = 1
+                        f = casadi.Function.bspline(
+                            symbol.name, knots, coeffs, degree, m
+                        )
+                        return f(casadi.hcat(converted_children).T).T
                     else:
                         LUT = casadi.interpolant(
                             "LUT", solver, symbol.x, symbol.y.ravel(order="F")