diff --git a/release_on_pip.sh b/release_on_pip.sh index bbeba024..e6858124 100644 --- a/release_on_pip.sh +++ b/release_on_pip.sh @@ -33,16 +33,21 @@ pip install --use-feature=2020-resolver --upgrade setuptools wheel twine auditwh pip install -e .['dev'] # Rebuild the .c files from the .pyc files -python setup.py develop +# python setup.py develop # Build the distribution python setup.py sdist bdist_wheel # Linux: -auditwheel repair --plat manylinux2014_x86_64 ./dist/*linux_x86_64.whl -w ./dist/ +# auditwheel repair --plat manylinux2014_x86_64 ./dist/*linux_x86_64.whl -w ./dist/ +# auditwheel repair --plat manylinux2014_x86_64 ./dist/*.whl -w ./dist/ + + # Upload the distro -twine upload dist/*manylinux*.whl dist/*tar.gz +# twine upload dist/*manylinux*.whl dist/*tar.gz +twine upload dist/*.whl dist/*tar.gz + # Windows (EDIT names to upload all things .whl and tat.gz in dist/ !!!): # twine upload dist/*manylinux*.whl dist/*tar.gz diff --git a/setup.py b/setup.py index 02cc133d..c18cba0b 100644 --- a/setup.py +++ b/setup.py @@ -59,7 +59,7 @@ def run(self): # Run the setup setup( name="tigramite", - version="5.0.0.0", + version="5.0.0.1", packages=["tigramite", "tigramite.independence_tests", "tigramite.toymodels"], license="GNU General Public License v3.0", description="Tigramite causal discovery for time series", diff --git a/tigramite/causal_effects.py b/tigramite/causal_effects.py index 69e92823..f92ea6cd 100644 --- a/tigramite/causal_effects.py +++ b/tigramite/causal_effects.py @@ -122,10 +122,10 @@ def __init__(self, # Construct internal graph from input graph depending on graph type # and hidden variables - (self.graph, self.graph_type, - self.tau_max, self.hidden_variables) = self._construct_graph( - graph=graph, graph_type=graph_type, - hidden_variables=hidden_variables) + # (self.graph, self.graph_type, + # self.tau_max, self.hidden_variables) = + self._construct_graph(graph=graph, graph_type=graph_type, + hidden_variables=hidden_variables) # print(self.graph.shape) self._check_graph(self.graph) @@ -212,20 +212,21 @@ def _construct_graph(self, graph, graph_type, hidden_variables): if graph_type in ['dag', 'admg']: - tau_max = 0 if graph.ndim != 2: raise ValueError("graph_type in ['dag', 'admg'] assumes graph.shape=(N, N).") # Convert to shape [N, N, 1, 1] with dummy dimension # to process as tsg_dag or tsg_admg with potential hidden variables self.graph = np.expand_dims(graph, axis=(2, 3)) + + # tau_max needed in _get_latent_projection_graph self.tau_max = 0 if len(hidden_variables) > 0: - graph = self._get_latent_projection_graph() # stationary=False) - graph_type = "tsg_admg" + self.graph = self._get_latent_projection_graph() # stationary=False) + self.graph_type = "tsg_admg" else: - graph = self.graph - graph_type = 'tsg_' + graph_type + # graph = self.graph + self.graph_type = 'tsg_' + graph_type elif graph_type in ['tsg_dag', 'tsg_admg']: if graph.ndim != 4: @@ -237,10 +238,10 @@ def _construct_graph(self, graph, graph_type, hidden_variables): self.tau_max = graph.shape[2] - 1 if len(hidden_variables) > 0: - graph = self._get_latent_projection_graph() #, stationary=False) - graph_type = "tsg_admg" + self.graph = self._get_latent_projection_graph() #, stationary=False) + self.graph_type = "tsg_admg" else: - graph_type = graph_type + self.graph_type = graph_type elif graph_type in ['stationary_dag']: # Currently only stationary_dag without hidden variables is supported @@ -260,16 +261,16 @@ def _construct_graph(self, graph, graph_type, hidden_variables): for varlag in self.X.union(self.Y).union(self.S): maxlag_XYS = max(maxlag_XYS, abs(varlag[1])) - tau_max = maxlag_XYS + statgraph_tau_max + self.tau_max = maxlag_XYS + statgraph_tau_max stat_graph = deepcopy(graph) # Construct tsg_graph - graph = np.zeros((self.N, self.N, tau_max + 1, tau_max + 1), dtype=' 1: + print("\n## Predicting target %s" % str(y)) + if pred_params is not None: + for key in list(pred_params): + print("%s = %s" % (key, pred_params[key])) + # Default value for pred_params + if pred_params is None: + pred_params = {} + # Check this is a valid target + if y not in self.model.fit_results: + raise ValueError("y = %s not yet fitted" % str(y)) + + # Transform the data if needed + a_transform = self.model.fit_results[y]['data_transform'] + if a_transform is not None: + intervention_data = a_transform.transform(X=intervention_data) + + + # Now iterate through interventions (and potentially S) + for index, dox_vals in enumerate(intervention_data): + # Construct XZS-array + intervention_array = dox_vals.reshape(1, lenX) + predictor_array = intervention_array + + predicted_vals = self.model.fit_results[y]['model'].predict( + X=predictor_array, **pred_params) + predicted_array[index, iy] = predicted_vals.mean() - return effect \ No newline at end of file + return predicted_array \ No newline at end of file diff --git a/tigramite/models.py b/tigramite/models.py index 6991641e..96fdd6e3 100644 --- a/tigramite/models.py +++ b/tigramite/models.py @@ -150,6 +150,7 @@ def get_general_fitted_model(self, # Initialize the fit results fit_results = {} for y in self.Y: + # Construct array of shape (var, time) with first entry being # a dummy, second is y followed by joint X and Z (ignore the notation in construct_array) array, xyz = \ @@ -164,9 +165,6 @@ def get_general_fitted_model(self, if self.data_transform is not None: array = self.data_transform.fit_transform(X=array.T).T - # Cache array for use in prediction - self.observation_array = array - self.xyz = xyz # Fit the model # Copy and fit the model a_model = deepcopy(self.model) @@ -179,8 +177,11 @@ def get_general_fitted_model(self, target_array = array[np.where(xyz==1)[0][0], :] a_model.fit(X=predictor_array, y=target_array) + # Cache the results fit_results[y] = {} + fit_results[y]['observation_array'] = array + fit_results[y]['xyz'] = xyz fit_results[y]['model'] = a_model # Cache the data transform fit_results[y]['data_transform'] = deepcopy(self.data_transform) @@ -254,13 +255,13 @@ def get_general_prediction(self, conditions_data = a_transform.transform(X=conditions_data) # Extract observational Z from stored array - z_indices = list(np.where(self.xyz==1)[0][1:]) - z_array = self.observation_array[z_indices, :].T + z_indices = list(np.where(self.fit_results[y]['xyz']==1)[0][1:]) + z_array = self.fit_results[y]['observation_array'][z_indices, :].T Tobs = len(z_array) if self.conditions is not None and conditions_data is not None: - s_indices = list(np.where(self.xyz==2)[0]) - s_array = self.observation_array[s_indices, :].T + s_indices = list(np.where(self.fit_results[y]['xyz']==2)[0]) + s_array = self.fit_results[y]['observation_array'][s_indices, :].T # Now iterate through interventions (and potentially S) for index, dox_vals in enumerate(intervention_data): @@ -293,133 +294,6 @@ def get_general_prediction(self, return predicted_array - # def get_general_prediction(self, - # intervention_data=None, - # conditions_data=None, - # pred_params=None, - # ): - # r"""Predict effect of intervention with fitted model. - - # Uses the model.predict() function of the sklearn model. - - # Parameters - # ---------- - # intervention_data : data object, optional - # Tigramite dataframe object with optional new mask. Only the - # values for X will be extracted. - # conditions_data : data object, optional - # Tigramite dataframe object with optional new mask. Only the - # values for conditions will be extracted. - # pred_params : dict, optional - # Optional parameters passed on to sklearn prediction function. - - # Returns - # ------- - # Results from prediction. - # """ - - # pred_dict = {} - # for y in self.Y: - # # Print message - # if self.verbosity > 1: - # print("\n## Predicting target %s" % str(y)) - # if pred_params is not None: - # for key in list(pred_params): - # print("%s = %s" % (key, pred_params[key])) - # # Default value for pred_params - # if pred_params is None: - # pred_params = {} - # # Check this is a valid target - # if y not in self.fit_results: - # raise ValueError("y = %s not yet fitted" % str(y)) - # # Construct the array form of the data - # # Check if we've passed a new dataframe object - # observation_array, xyz = \ - # self.dataframe.construct_array(X=self.X, Y=[y] + self.Z, Z=self.conditions, - # tau_max=self.tau_max, - # # mask=self.test_mask, - # mask_type=self.mask_type, - # cut_off=self.cut_off, - # verbosity=self.verbosity) - - # intervention_array = np.copy(observation_array) - # if intervention_data is not None: - # tmp_array, _ = intervention_data.construct_array(X=self.X, Y=[y] + self.Z, - # Z=self.conditions, - # tau_max=self.tau_max, - # mask_type=self.mask_type, - # cut_off=self.cut_off, - # verbosity=self.verbosity) - - # # Only replace X-variables in intervention_array (necessary if lags of - # # X are in Z...) - # for index in np.where(xyz==0)[0]: - # intervention_array[index] = tmp_array[index] - - # if self.conditions is not None and conditions_data is not None: - # tmp_array, _ = conditions_data.construct_array(X=self.X, Y=[y] + self.Z, - # Z=self.conditions, - # tau_max=self.tau_max, - # mask_type=self.mask_type, - # cut_off=self.cut_off, - # verbosity=self.verbosity) - - # # Only replace condition-variables in intervention_array - # # (necessary if lags of X are in Z...) - # for index in np.where(xyz==2)[0]: - # intervention_array[index] = tmp_array[index] - - # # Transform the data if needed - # a_transform = self.fit_results[y]['data_transform'] - # if a_transform is not None: - # intervention_array = a_transform.transform(X=intervention_array.T).T - # # Cache the test array - # self.intervention_array = intervention_array - # # Run the predictor, for Y only the Z-part is used, the first index is y - # predictor_indices = list(np.where(xyz==0)[0]) \ - # + list(np.where(xyz==1)[0][1:]) \ - # + list(np.where(xyz==2)[0]) - # predictor_array = intervention_array[predictor_indices, :].T - - # pred_dict[y] = self.fit_results[y]['model'].predict( - # X=predictor_array, **pred_params) - - # # print(pred_dict[y]) - # if self.conditions is not None and conditions_data is not None: - - # a_conditional_model = deepcopy(self.conditional_model) - - # # Fit Y|do(X) on S - # conditions_array = observation_array[list(np.where(xyz==2)[0])] - # target_array = pred_dict[y] # array[np.where(xyz==1)[0][0], :] - - # if a_transform is not None: - # conditions_array = a_transform.transform(X=conditions_array.T).T - # target_array = a_transform.transform(X=target_array.T).T - - # a_conditional_model.fit(X=conditions_array.T, y=target_array) - # self.fit_results[y]['conditional_model'] = a_conditional_model - - # # Now predict conditional causal effect for new conditions - # tmp_array, _ = conditions_data.construct_array(X=self.X, Y=[y] + self.Z, - # Z=self.conditions, - # tau_max=self.tau_max, - # mask_type=self.mask_type, - # cut_off=self.cut_off, - # verbosity=self.verbosity) - - # # Construct conditions array - # new_conditions_array = tmp_array[list(np.where(xyz==2)[0])] - - # if a_transform is not None: - # new_conditions_array = a_transform.transform(X=new_conditions_array.T).T - - # pred_dict[y] = a_conditional_model.predict( - # X=new_conditions_array.T, **pred_params) - - # return pred_dict - - def get_fit(self, all_parents, selected_variables=None, tau_max=None, diff --git a/tutorials/tigramite_tutorial_general_causal_effect_analysis.ipynb b/tutorials/tigramite_tutorial_general_causal_effect_analysis.ipynb index 53901f1d..3934e624 100644 --- a/tutorials/tigramite_tutorial_general_causal_effect_analysis.ipynb +++ b/tutorials/tigramite_tutorial_general_causal_effect_analysis.ipynb @@ -389,7 +389,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 2, "metadata": {}, "outputs": [ { @@ -409,7 +409,7 @@ "S = []\n", "M = []\n", "\n", - "hidden_variables = [(0, 0)]\n", + "hidden_variables = {(0, 0)}\n", "\n", "\n", "\n" @@ -438,7 +438,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 3, "metadata": {}, "outputs": [ { @@ -471,7 +471,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 4, "metadata": {}, "outputs": [ { @@ -546,7 +546,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 5, "metadata": {}, "outputs": [ { @@ -590,7 +590,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 6, "metadata": {}, "outputs": [ { @@ -625,7 +625,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 7, "metadata": {}, "outputs": [ { @@ -660,7 +660,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 8, "metadata": {}, "outputs": [ { @@ -678,7 +678,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 9, "metadata": {}, "outputs": [ { @@ -722,7 +722,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 10, "metadata": {}, "outputs": [ { @@ -740,7 +740,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 11, "metadata": {}, "outputs": [ { @@ -784,7 +784,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 12, "metadata": {}, "outputs": [ { @@ -837,7 +837,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 13, "metadata": {}, "outputs": [ { @@ -874,7 +874,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ @@ -904,7 +904,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 15, "metadata": { "scrolled": true }, @@ -964,7 +964,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 17, "metadata": {}, "outputs": [ { @@ -997,18 +997,14 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 18, "metadata": {}, "outputs": [ { - "ename": "NameError", - "evalue": "name 'intervention_data1' is not defined", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m/tmp/ipykernel_4532/2201758403.py\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mtrue_effect\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mintervention_data1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mY\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mintervention_data2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mY\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"True effect = %.2f\"\u001b[0m \u001b[0;34m%\u001b[0m\u001b[0mtrue_effect\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mNameError\u001b[0m: name 'intervention_data1' is not defined" + "name": "stdout", + "output_type": "stream", + "text": [ + "True effect = 0.75\n" ] } ], @@ -1037,16 +1033,16 @@ }, { "cell_type": "code", - "execution_count": 45, + "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 45, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" } @@ -1073,15 +1069,15 @@ }, { "cell_type": "code", - "execution_count": 46, + "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "y1 = [[5.00149575]]\n", - "y2 = [[4.36091033]]\n" + "y1 = [[0.73288076]]\n", + "y2 = [[-0.00799763]]\n" ] } ], @@ -1110,14 +1106,14 @@ }, { "cell_type": "code", - "execution_count": 47, + "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Causal effect = 0.64\n" + "Causal effect = 0.74\n" ] } ], @@ -1137,12 +1133,12 @@ }, { "cell_type": "code", - "execution_count": 48, + "execution_count": 22, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAe7ElEQVR4nO3dfZQcdZ3v8fcnDzohJOYKA0KCCbJsJBCchBGJWTkgaEiI4UF2YY+uRpFsRK/oHqIEXVRcr3jRKzeyGnEvPqy4PmcE5CEiwciTYZJMAAlIQFiTcMmQu4kZTDCB7/2japKm09PdM9PV3dP9eZ3TZ6qraqq+XUn6m6pffb+liMDMzJrXsFoHYGZmteVEYGbW5JwIzMyanBOBmVmTcyIwM2tyI2odQH8dfPDBMWnSpFqHYWY2pKxevfq5iGgttGzIJYJJkybR2dlZ6zDMzIYUSU/3tcyXhszMmpwTgZlZk3MiMDNrckNujKCQ3bt3s3HjRnbt2lXrUBpKS0sLEyZMYOTIkbUOxcwy1BCJYOPGjYwZM4ZJkyYhqdbhNISIYOvWrWzcuJEjjzyy1uGYWYYaIhHs2rXLSaDCJHHQQQfR3d1d61DMml7H2k1cfftjbN62k8PHjWLRrMmcPW18xbbfEIkAcBLIgI+pWe11rN3E4p89xM7dLwKwadtOFv/sIYCKJQMPFpuZ1bGrb39sbxLotXP3i1x9+2MV24cTQYUMHz6ctra2va+rrrqqz3U7Ojp45JFH9r6/4ooruOOOOwYdw7Zt2/ja177W79/7zGc+w5e+9KVB79/MKm/ztp39mj8QDXNpqNZGjRpFV1dXWet2dHQwd+5cpkyZAsCVV15ZkRh6E8HFF19cke2ZWe0dPm4Umwp86R8+blTF9tGUZwQdazcx86o7OfKyXzDzqjvpWLsps31ddtllTJkyheOPP55LL72Ue++9lxtvvJFFixbR1tbGE088wfz58/nJT34CJC00Lr/8cmbMmEF7eztr1qxh1qxZHHXUUSxduhSAnp4eTjvtNKZPn87UqVP5+c9/vndfTzzxBG1tbSxatAiAq6++mje+8Y0cf/zxfPrTn94b1+c//3kmT57M6aefzmOPVe4U08wqa9GsyYwaOfxl80aNHM6iWZMrto+mOyPIauBl586dtLW17X2/ePFi3va2t7Fs2TIeffRRJLFt2zbGjRvHvHnzmDt3Luedd17BbR1xxBHcd999fOxjH2P+/Pncc8897Nq1i2OPPZaFCxfS0tLCsmXLGDt2LM899xwnnXQS8+bN46qrruLhhx/ee2ayfPlyHn/8cVatWkVEMG/ePFauXMno0aP5wQ9+wNq1a9mzZw/Tp0/nhBNOGPBnN7Ps9H4v+a6hCio28DKYA1vo0tCePXtoaWnhAx/4AGeeeSZz584ta1vz5s0DYOrUqfT09DBmzBjGjBlDS0sL27ZtY/To0Vx++eWsXLmSYcOGsWnTJp599tn9trN8+XKWL1/OtGnTgORM4vHHH2fHjh2cc845HHDAAS/bn5nVp7Onja/oF3++prs0VI2Bl14jRoxg1apVvPOd76Sjo4MzzjijrN975StfCcCwYcP2Tve+37NnDzfccAPd3d2sXr2arq4uDj300IJV1RHB4sWL6erqoquriw0bNnDhhRcCvjXUzPZpukTQ1wBLJQdeevX09LB9+3bmzJnDNddcs/eMYcyYMezYsWPA292+fTuHHHIII0eOZMWKFTz99NMFtztr1iyuv/56enp6ANi0aRNbtmzh5JNPZtmyZezcuZMdO3Zw0003DfxDmtmQ13SXhhbNmvyyMQKozMBL/hjBGWecwSWXXMJZZ53Frl27iAi+8pWvAHDBBRdw0UUXsWTJkr2DxP3xrne9i3e84x20t7fT1tbG61//egAOOuggZs6cyXHHHcfs2bO5+uqrWb9+PTNmzADgwAMP5Hvf+x7Tp0/n/PPPp62tjYkTJ/KWt7xlUJ/dzIY2RUR2G5eeAnYALwJ7IqI9b/kpwM+BP6SzfhYRRe+lbG9vj/wH06xfv55jjjmm7LiyLtduJP09tma2v3r4zpG0Ov87uFc1zghOjYjniiz/TUSUN4paIVkPvJiZ9apGi4jBaroxAjOzaqpGi4jByjoRBLBc0mpJC/pYZ4akdZJulXRsoRUkLZDUKamzr26YWV7ialY+pmaDV807FQcq60QwMyKmA7OBD0k6OW/5GmBiRLwB+CrQUWgjEXFdRLRHRHtra+t+y1taWti6dau/uCqo93kELS0ttQ7FbEir5p2KA5XpGEFEbE5/bpG0DDgRWJmz/E8507dI+pqkg0uMKexnwoQJbNy40b3zK6z3CWVmNnBZ3alYSZklAkmjgWERsSOdfjtwZd46rwGejYiQdCLJGcrW/u5r5MiRfoqWmdWlarSIGKwszwgOBZalFawjgO9HxG2SFgJExFLgPOCDkvYAO4ELwtd3zKzB1Pudipklgoh4EnhDgflLc6avBa7NKgYzMyvNt4+amTW5pmsxYWbWX/VQGZwlJwIzsyKGQmXwYPnSkJlZEUOhMniwnAjMzIoYCpXBg+VEYGZWxFCoDB4sJwIzsyKq8fD4WvNgsZlZEUOhMniwnAjMzEqo98rgwfKlITOzJudEYGbW5HxpyMwaXqNXBg+WE4GZNbRmqAweLF8aMrOG1gyVwYPlRGBmDa0ZKoMHy4nAzBpaM1QGD5YTgZk1tGaoDB4sDxabWUNrhsrgwXIiMLOG1+iVwYPlS0NmZk3OicDMrMk5EZiZNTmPEZhZ3XOLiGw5EZhZXXOLiOz50pCZ1TW3iMieE4GZ1TW3iMieE4GZ1TW3iMieE4GZ1TW3iMieB4vNrK65RUT2nAjMrO65RUS2Mr00JOkpSQ9J6pLUWWC5JC2RtEHSg5KmZxmPmZntrxpnBKdGxHN9LJsNHJ2+3gR8Pf1pZmZVUutLQ2cB342IAO6XNE7SYRHxTI3jMrMKcmVwfcv6rqEAlktaLWlBgeXjgT/mvN+YznsZSQskdUrq7O7uzihUM8tCb2Xwpm07CfZVBnes3VTr0CyVdSKYGRHTSS4BfUjSyXnLVeB3Yr8ZEddFRHtEtLe2tmYRp5llxJXB9S/TRBARm9OfW4BlwIl5q2wEjsh5PwHYnGVMZlZdrgyuf5klAkmjJY3pnQbeDjyct9qNwHvSu4dOArZ7fMCssbgyuP5leUZwKHC3pHXAKuAXEXGbpIWSFqbr3AI8CWwAvglcnGE8ZlYDrgyuf5ndNRQRTwJvKDB/ac50AB/KKgYzqz1XBte/Wt8+amZNwJXB9c1N58zMmpwTgZlZk3MiMDNrch4jMLOS3CKisTkRmFlRfnh84/OlITMryi0iGp8TgZkV5RYRjc+JwMyKcouIxudEYGZFuUVE4/NgsZkV5RYRjc+JwMxKcouIxlby0pCkvy1nnpmZDU3ljBEsLnOemZkNQX1eGpI0G5gDjJe0JGfRWGBP1oGZWeW4MtiKKTZGsBnoBOYBq3Pm7wA+lmVQZlY5rgy2UvpMBBGxDlgnaRnwfES8CCBpOPDKKsVnZoNUrDLYicCgvDGC5UBu5cgo4I5swjGzSnNlsJVSTiJoiYie3jfp9AHZhWRmleTKYCulnETwvKTpvW8knQD4vxJmQ4Qrg62UcgrKPgr8WNLm9P1hwPmZRWRmFeXKYCulZCKIiAckvR6YDAh4NCJ2Zx6ZmVWMK4OtmHIqiw8APgFcEhEPAZMkzc08MjMzq4pyxgi+BfwFmJG+3wj8S2YRmZlZVZUzRnBURJwv6e8BImKnJGUcl5nlcGWwZamcRPAXSaOAAJB0FPBCplGZ2V6uDLaslXNp6NPAbcARkm4AfgV8PNOozGwvPzPYslas6dzMiLgHWAmcC5xEctfQJRHxXJXiM2t6rgy2rBU7I+jtOHpfRGyNiF9ExM1OAmbV5cpgy1qxRLBb0reACZKW5L/K3YGk4ZLWSrq5wLJTJG2X1JW+rhjIhzBrZK4MtqwVGyyeC5wOvJWXt6Hur0uA9STPMSjkNxHhugSzPrgy2LJWLBEsiohPSHptRHxnIBuXNAE4E/g88E8D2YaZuTLYslXs0tAcSSOBCwax/WtI7jB6qcg6MyStk3SrpGMLrSBpgaROSZ3d3d2DCMfMzPIVSwS3Ac8Bx0v6k6QduT9LbThtQ7ElIopdVloDTIyINwBfBToKrRQR10VEe0S0t7a2ltq1mZn1Q5+JICIWRcSrgF9ExNiIGJP7s4xtzwTmSXoK+AHwVknfy9vHn3qfdRARtwAjJR084E9jZmb9Vk730bMkTQSOjog70irjERGxo8TvLQYWQ3J3EHBpRLw7dx1JrwGejYiQdCJJYto6oE9iVsfcIsLqWclEIOkiYAHwauAoYAKwFDhtIDuUtBAgIpYC5wEflLSH5GE3F0REDGS7ZvXKLSKs3qnU966kLuBE4LcRMS2d91BETM0+vP21t7dHZ2dnLXZtNiAzr7qTTQWqgMePG8U9l721BhFZM5K0OiLaCy0rp9fQCxHxl5yNjSBtQGdmpblFhNW7chLBryVdDoyS9Dbgx8BN2YZl1jjcIsLqXTmJ4DKgG3gI+EfgFuBTWQZl1kjcIsLqXTl3Db0EfDN9mVk/uUWE1btyHkxjZoPkFhFWz8q5NGRmZg3MicDMrMkVe0LZTRS5TTQi5mUSkVkdcmWwNbJiYwRfqloUZnXMlcHW6PpMBBHx62oGYlavij083onAGkE5vYaOBr4ATAFaeudHxOsyjMusbrgy2BpdOYPF3wK+DuwBTgW+C/x7lkGZ1RNXBlujKycRjIqIX5E0qHs6Ij5D8hxjs6bgymBrdOUUlO2SNAx4XNKHgU3AIdmGZVY/XBlsja6cRPBR4ADgI8DnSM4G3pthTGZ1x5XB1sjK6TX0QDrZI+lC4MCIKPnMYjMzGxpKjhFI+r6ksZJGA48Aj0lalH1oZmZWDeUMFk9JzwDOJmlB/VrgH7IMyszMqqecMYKRkkaSJIJrI2K3JD+hzIYUt4gw61s5ZwTfAJ4CRgMrJU0EPEZgQ0Zvi4hN23YS7GsR0bF2U61DM6sLJRNBRCyJiPERMScST5MUlpkNCcVaRJhZeS0mruhj0ZUVjsUsE24RYVZcOZeGns95vQjMBiZlGJNZRblFhFlx5dQRfDn3vaQvATdmFpFZhS2aNfllbaTBLSLMcg3kmcUHAO48akOGW0SYFVfOGMFD7HtS2XCgFY8P2BDjFhFmfSvnjGBuzvQe4NmI2JNRPGZmVmXljBE8DSDpEJIH0xwuiYj4z6yDMzOz7JVzaWge8GXgcGALMBFYDxybbWhm+7gy2Cw75dw++jngJOD3EXEkcBpwT6ZRmeVwZbBZtspJBLsjYiswTNKwiFgBtJW7A0nDJa2VdHOBZZK0RNIGSQ9Kml5+6NYsXBlslq1yBou3SToQWAncIGkLyaBxuS4huZQ0tsCy2cDR6etNJM9GflM/tm1NwJXBZtkq54zgLODPwMeA24AngHeUs3FJE4AzgX8rsu3vpj2M7gfGSTqsnG1b83BlsFm2+kwEkv5K0syIeD4iXoqIPRHxHaALGFfm9q8BPg681Mfy8cAfc95vTOflx7JAUqekzu7u7jJ3bY3CD483y1axM4JrgB0F5v85XVaUpLnAlohYXWy1AvP2e9ZBRFwXEe0R0d7a2lpq19Zgzp42ni+cO5Xx40YhYPy4UXzh3Km+a8isQoqNEUyKiAfzZ0ZEp6RJZWx7JjBP0hyS+oOxkr4XEe/OWWcjcETO+wnA5jK2bU3GlcFm2Sl2RtBSZFnJi7MRsTgiJkTEJOAC4M68JABJ87r3pHcPnQRsj4hnSm3bzMwqp1gieEDSRfkzJV0IFLvcU5SkhZIWpm9vAZ4ENgDfBC4e6HbNzGxgFFH48cOSDgWWAX9h3xd/O/AK4JyI+L9ViTBPe3t7dHZ21mLXNgiuDDarLUmrI6K90LI+xwgi4lngzZJOBY5LZ/8iIu7MIEZrYL2Vwb1FYb2VwYCTgVkdKKfp3ApgRRVisQZVrDLYicCs9sopKDMbFFcGm9U3JwLLnCuDzeqbE4FlzpXBZvVtIM8sNusXPzPYrL45EVhVuDLYrH750pCZWZNzIjAza3JOBGZmTc5jBFYWt4gwa1xOBFaSW0SYNTZfGrKS/PB4s8bmRGAluUWEWWNzIrCS3CLCrLE5EVhJbhFh1tg8WGwluUWEWWNzIrCyuEWEWePypSEzsybnRGBm1uR8aagJuCrYzIpxImhwrgo2s1J8aajBuSrYzEpxImhwrgo2s1KcCBqcq4LNrBQnggbnqmAzK8WDxQ3OVcFmVooTQRNwVbCZFeNLQ2ZmTS6zRCCpRdIqSesk/U7SZwusc4qk7ZK60tcVWcVjZmaFZXlp6AXgrRHRI2kkcLekWyPi/rz1fhMRczOMw8zMisgsEUREAD3p25HpK7LaXyNziwgzy1KmYwSShkvqArYAv4yI3xZYbUZ6+ehWScf2sZ0FkjoldXZ3d2cZct3pbRGxadtOgn0tIjrWbqp1aGbWIDJNBBHxYkS0AROAEyUdl7fKGmBiRLwB+CrQ0cd2rouI9ohob21tzTLkuuMWEWaWtarcNRQR24C7gDPy5v8pInrS6VuAkZIOrkZMQ4VbRJhZ1rK8a6hV0rh0ehRwOvBo3jqvkaR0+sQ0nq1ZxTQUuUWEmWUtyzOCw4AVkh4EHiAZI7hZ0kJJC9N1zgMelrQOWAJckA4yW8otIswsa1neNfQgMK3A/KU509cC12YVQyNwiwgzy5pbTAwBbhFhZllyiwkzsybnRGBm1uR8aagKXBlsZvXMiSBjfni8mdU7XxrKmCuDzazeORFkzJXBZlbvnAgy5spgM6t3TgQZc2WwmdU7DxZnzJXBZlbvnAiqwJXBZlbPfGnIzKzJORGYmTU5JwIzsybnMYIyuEWEmTUyJ4IS3CLCzBqdLw2V4BYRZtbonAhKcIsIM2t0TgQluEWEmTU6J4IS3CLCzBqdB4tLcIsIM2t0TgRlcIsIM2tkvjRkZtbknAjMzJpcU1wacmWwmVnfGj4RuDLYzKy4hr805MpgM7PiGj4RuDLYzKy4hk8Ergw2Myuu4ROBK4PNzIrLLBFIapG0StI6Sb+T9NkC60jSEkkbJD0oaXql4zh72ni+cO5Uxo8bhYDx40bxhXOneqDYzCyV5V1DLwBvjYgeSSOBuyXdGhH356wzGzg6fb0J+Hr6s6JcGWxm1rfMzggi0ZO+HZm+Im+1s4DvpuveD4yTdFhWMZmZ2f4yHSOQNFxSF7AF+GVE/DZvlfHAH3Peb0zn5W9ngaROSZ3d3d2ZxWtm1owyTQQR8WJEtAETgBMlHZe3igr9WoHtXBcR7RHR3tramkGkZmbNqyp3DUXENuAu4Iy8RRuBI3LeTwA2VyMmMzNLZHnXUKukcen0KOB04NG81W4E3pPePXQSsD0inskqJjMz21+Wdw0dBnxH0nCShPOjiLhZ0kKAiFgK3ALMATYAfwbeV2qjq1evfk7S0wOM6WDguQH+bpbqNS6o39gcV/84rv5pxLgm9rVAEftdkm9Ykjojor3WceSr17igfmNzXP3juPqn2eJq+MpiMzMrzonAzKzJNVsiuK7WAfShXuOC+o3NcfWP4+qfpoqrqcYIzMxsf812RmBmZnmcCMzMmlxDJwJJV0t6NG1xvay3wK3AemdIeixth31ZFeL627Q190uS+rwVTNJTkh6S1CWps47iqvbxerWkX0p6PP353/pYryrHq9Tnr0Z79QHGdYqk7enx6ZJ0RZXiul7SFkkP97G8VserVFy1Ol5HSFohaX367/GSAutU9phFRMO+gLcDI9LpLwJfLLDOcOAJ4HXAK4B1wJSM4zoGmEzSdqO9yHpPAQdX8XiVjKtGx+t/Apel05cV+nOs1vEq5/OTFEneStJL6yTgt1X4sysnrlOAm6v19ylnvycD04GH+1he9eNVZly1Ol6HAdPT6THA77P+O9bQZwQRsTwi9qRv7yfpZZTvRGBDRDwZEX8BfkDSHjvLuNZHxGNZ7mMgyoyr6scr3f530unvAGdnvL9iyvn8tWivXos/l7JExErg/xVZpSbt6MuIqyYi4pmIWJNO7wDWs39X5ooes4ZOBHneT5JB85XVCrtGAlguabWkBbUOJlWL43VopD2o0p+H9LFeNY5XOZ+/Fseo3H3OUPLUwFslHZtxTOWq53+DNT1ekiYB04ABtfAvV5a9hqpC0h3Aawos+mRE/Dxd55PAHuCGQpsoMG/Q99SWE1cZZkbEZkmHAL+U9Gj6v5haxlX149WPzVT8eBVQzufP5BiVUM4+1wATI3lq4Bygg+TpgLVWi+NVjpoeL0kHAj8FPhoRf8pfXOBXBnzMhnwiiIjTiy2X9F5gLnBapBfX8mTSCrtUXGVuY3P6c4ukZSSn/4P6YqtAXFU/XpKelXRYRDyTnv5u6WMbFT9eBZTz+WvRXr3kPnO/TCLiFklfk3RwRNS6uVpdtqOv5fFS8njfnwI3RMTPCqxS0WPW0JeGJJ0BfAKYFxF/7mO1B4CjJR0p6RXABSTtsWtK0mhJY3qnSQa+C97dUGW1OF43Au9Np98L7HfmUsXjVc7nr0V79ZJxSXqNJKXTJ5L8+9+acVzlqMt29LU6Xuk+/w+wPiL+Vx+rVfaYVXtEvJovkvbWfwS60tfSdP7hwC05680hGZl/guQSSdZxnUOS0V8AngVuz4+L5O6Pdenrd/USV42O10HAr4DH05+vruXxKvT5gYXAwnRawL+myx+iyJ1hVY7rw+mxWUdy88SbqxTXfwDPALvTv18X1snxKhVXrY7X35Bc5nkw57trTpbHzC0mzMyaXENfGjIzs9KcCMzMmpwTgZlZk3MiMDNrck4EZmZNzonABkxSTxnrfFTSAdWIp4/9t6VVob3v56mCHVMlfUPSzLx5SyT9c877T0r61wrs6y4V6QprNlC+fdQGTFJPRBxYYp2nSO5xLrsaU9LwiHhxsPGl25qf7v/Dldhege13ASfkxitpLMm936eT3A9+JzAtIrYNcl93AZdGRMVabEsaEfsaM1qT8hmBDVrat/0uST9R8vyHG9KKx4+QFH2tkLQiXfftku6TtEbSj9N+Kr3PErhC0t3AxyWtytn+JEkPptMnSPp12lju9rTlRO//lr8oaZWk30t6S1pheyVwvpJ+8udLmi/p2vR3Jkr6lZJ+7r+S9Np0/rfT/9XfK+lJSef18bmPAX6fn7QiaU3wSeBakqKfK/KTgKTZkn6UdwxvSqe/LqlTSS/6z/ax756c6fMkfTudbpX0U0kPpK+ZBX53fnrsbyJp0ndg+vnXKHmew1k5x329pG+msSyXNCpd9sb0uN2n5LkfD6fzh6fvH0iX/2Oh+K3OVKNSzq/GfAE96c9TgO0k/U6GAfcBf5Mue4r0GQHAwSS9f0an7z9B8iXZu97Hc7bdBbwuZ71PASOBe4HWdP75wPXp9F3Al9PpOcAd6fR84Nqc7e59D9wEvDedfj/QkU5/G/hx+lmmkLR3LvT5/wl4f5Hjcx9wdx/LRgD/mXMsvg68O53urZwenn6u43M+Y3vusU+nzwO+nU5/P+fYv5akTUH+vueTVNK+OieWsTl/RhtIKlcnkTRrbEuX/SgnxodJK22Bq0h7+gMLgE+l068EOoEja/131a/iryHfdM7qxqqI2Ah7L5dMAu7OW+ckki/We9IWLq8g+bLs9cOc6R8Bf0fyJXN++poMHEfSWRSSL8rc/iq9zblWp/svZQZwbjr97yQPwOnVEREvAY9IOrSP358FvK/QAkkTSLqphqQDI+Jl4ykRsUfSbcA7JP0EOBP4eLr475S00R5B8pCSKSTtBspxOjAlPT4AYyWNiaSvfa5fRkRvL34B/0PSycBLJO2Mez/zHyKiK51eDUxS8qS/MRFxbzr/+ySNHSHp8XR8zlnUq0g6dv6hzPitBpwIrFJeyJl+kcJ/t0TyBfT3fWzj+ZzpHwI/lvQzICLicUlTgd9FxIwSMfS1/1JyB8xyP89+LX/TAfBxkXY8LeB/A58heerbp4FFkj5P8oVPRLSRfMYPkTwc5YGI2CHpSOBS4I0R8V/pJZ+WErHmLh8GzIiInX3E1Sv3WL8LaCUZ69idjuv0bjP/z3UUhVsg9xLw3yPi9hL7tzriMQLL2g6Sx+1B0rhrpqS/guTLVNJfF/qliHiC5Ivnn9l3pvAY0CppRvr7I1X6YSG5+893L0mXTki+DPPPYIo5FVhRaIGk2SQPz/ku8DngHElTIuKTEdGWJgFILvVMBy5i32ccS/IlvT09E5ndx/6flXSMpGEkzQJ7LSdpltYbS1v+LxbwKmBLmgROBSYWWzki/gvYoaTrJew7hgC3Ax9U0kYZSX+tpBus1TEnAsvadcCtklZERDfJ9en/SAd/7wdeX+R3fwi8m+QyEZE8gvE84IuS1pGMI7y5xP5XkFwq6ZJ0ft6yjwDvS2P5B2C/h4QXMRu4LX+mpBbgGuDiSDxPcsnn2vx1Ixlkvjnd1s3pvHXAWpKul9cD9/Sx/8vS37mTl18e+wjQng7UPkLSsbKUG9Lf6SRJiI+W8TsXAtdJuo/kLGB7Ov/fgEeANekA8jfwlYe659tHzQZA0hrgTRGxu9ax1ELuuIeSuozDIqI/idTqiBOBmfVbena1mOR/+08D89MzPhuCnAjMzJqcxwjMzJqcE4GZWZNzIjAza3JOBGZmTc6JwMysyf1/VF9Uc2mDfQgAAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgCklEQVR4nO3df5wcdZ3n8dc7IZghBHNCQEiAoMtGAsEhjkjM6kNXMBBjAoiC5+6BollOOX/cw2gAFz09Dzz09BAV2T0EVxR/rImgSCKImxVwIYEEAiESENaZ8CCBNSHBRBP43B9VHTpDd03PdFdXTc/7+Xj0Y6qraqo+XUn6k6pvfT6liMDMzKyeUUUHYGZm5eZEYWZmmZwozMwskxOFmZllcqIwM7NMexUdQB4OOOCAmDJlStFhmJkNGytXrnwqIibWWtaRiWLKlCmsWLGi6DDMzIYNSY/XW+ZLT2ZmlsmJwszMMjlRmJlZpkLHKCRdDcwFNkbEMTWWvwn4CfC7dNaPI+KzQ9nXzp076e3tZceOHUOM1vobO3YskydPZsyYMUWHYmY5Know+xrgCuDbGev8a0TMbXZHvb29jB8/nilTpiCp2c2NeBHB008/TW9vL0cccUTR4ZhZjgpNFBGxXNKUduxrx44dThItJIn999+fTZs2FR2K2Yi35N4+Llu6jg2bt3PIhC4Wzp7KqcdNatn2h8MYxUxJqyX9XNLR9VaStEDSCkkr6n15OUm0lo+nWfGW3NvHBT++n77N2wmgb/N2Lvjx/Sy5t69l+yh7orgHODwiXg18FVhSb8WIuCoieiKiZ+LEmjUjZmYd57Kl69i+87k95m3f+RyXLV3Xsn2UOlFExDMRsS2dvgkYI+mAgsMastGjR9Pd3b37demll9Zdd8mSJTz44IO731988cXccsstTcewefNmvv71rw/69z7zmc/wxS9+sen9m1lrbdi8fVDzh6LowexMkl4OPBkRIel4ksT2dMFhDVlXVxerVq1qaN0lS5Ywd+5cpk2bBsBnPzukm71epJIoPvjBD7Zke2ZWrEMmdNFXIykcMqGrZfso9IxC0veAO4GpknolnSvpPEnnpaucAayRtBq4HDgr2vRIviX39jHr0l9yxKKfMevSX7b0el9/ixYtYtq0aRx77LF8/OMf54477uCGG25g4cKFdHd388gjj3DOOefwox/9CEhalFx44YXMnDmTnp4e7rnnHmbPns0rX/lKrrzySgC2bdvGW97yFmbMmMH06dP5yU9+sntfjzzyCN3d3SxcuBCAyy67jNe+9rUce+yxfPrTn94d1+c//3mmTp3KiSeeyLp1rTuNNbPWWTh7Kl1jRu8xr2vMaBbOntqyfRR919O7B1h+Bcnts21VGRyqXPerDA4BTd1JsH37drq7u3e/v+CCCzjppJNYvHgxDz30EJLYvHkzEyZMYN68ecydO5czzjij5rYOPfRQ7rzzTj72sY9xzjnncPvtt7Njxw6OPvpozjvvPMaOHcvixYvZb7/9eOqppzjhhBOYN28el156KWvWrNl9ZrNs2TIefvhh7rrrLiKCefPmsXz5csaNG8f111/Pvffey65du5gxYwavec1rhvzZzSwfle+kPO96KvWlp6JkDQ41c/BrXXratWsXY8eO5f3vfz9ve9vbmDu3sZKRefPmATB9+nS2bdvG+PHjGT9+PGPHjmXz5s2MGzeOCy+8kOXLlzNq1Cj6+vp48sknX7SdZcuWsWzZMo477jggORN5+OGH2bp1K6eddhr77LPPHvszs/I59bhJLU0M/ZV6MLso7Rgcqthrr7246667eMc73sGSJUs4+eSTG/q9l7zkJQCMGjVq93Tl/a5du7juuuvYtGkTK1euZNWqVRx00EE1q9IjggsuuIBVq1axatUq1q9fz7nnngv49lczSzhR1FBvEKiVg0MV27ZtY8uWLcyZM4evfOUru884xo8fz9atW4e83S1btnDggQcyZswYbrvtNh5//PGa2509ezZXX30127ZtA6Cvr4+NGzfyxje+kcWLF7N9+3a2bt3KjTfeOPQPaWbDmi891bBw9tQ9xiigNYND/ccoTj75ZD7ykY8wf/58duzYQUTw5S9/GYCzzjqLD3zgA1x++eW7B7EH4z3veQ9vf/vb6enpobu7m1e96lUA7L///syaNYtjjjmGU045hcsuu4y1a9cyc+ZMAPbdd1++853vMGPGDM4880y6u7s5/PDDecMb3tDUZzez4UttuomorXp6eqL/g4vWrl3LUUcd1fA28i6J7xSDPa5m9mJl+L6RtDIiemot8xlFHXkPDpmZQX53WbaSxyjMzArUjhYczRpRiaITL7MVycfTrHntvMtyqEZMohg7dixPP/20v9xapPI8irFjxxYditmw1s67LIdqxIxRTJ48md7eXj8/oYUqT7gzs6HL6y7LVhoxiWLMmDF+EpuZlU47WnA0a8QkCjOzsir7XZYjZozCzMyGxonCzMwy+dKTmVmTylBZnScnCjOzJgyHyupm+dKTmVkThkNldbOcKMzMmjAcKqub5URhZtaE4VBZ3SwnCjOzJiycPZWuMaP3mFe2yupmeTDbzKwJw6GyullOFGZmTSp7ZXWzfOnJzMwyOVGYmVkmX3oysxGv0yurm1XoGYWkqyVtlLSmznJJulzSekn3SZrR7hjNrLNVKqv7Nm8neKGyesm9fUWHVhpFX3q6Bjg5Y/kpwJHpawHwjTbEZGYjyEiorG5WoYkiIpYD/5Gxynzg25H4DTBB0sHtic7MRoKRUFndrKLPKAYyCfh91fvedN6LSFogaYWkFX7cqZk1aiRUVjer7IlCNeZFrRUj4qqI6ImInokTJ+Yclpl1ipFQWd2sst/11AscWvV+MrChoFjMrAONhMrqZpU9UdwAnC/peuB1wJaIeKLgmMysw3R6ZXWzCk0Ukr4HvAk4QFIv8GlgDEBEXAncBMwB1gN/BN5bTKRmZiNXoYkiIt49wPIAPtSmcMzMrIayD2abmVnByj5GYWY2ILfgyJcThZkNa5UWHJXq6koLDsDJokV86cnMhjW34MifE4WZDWtuwZE/JwozG9bcgiN/ThRmNqy5BUf+PJhtZsOaW3Dkz4nCzIY9t+DIly89mZlZJicKMzPL5EtPZlY4V1aXmxOFmRXKldXl50tPZlYoV1aXnxOFmRXKldXl50RhZoVyZXX5OVGYWaFcWV1+Hsw2s0K5srr8nCjMrHCurC43X3oyM7NMThRmZpbJicLMzDJ5jMLMmuYWHJ3NicLMmuIWHJ2v0EtPkk6WtE7SekmLaix/k6Qtklalr4uLiNPM6nMLjs5X2BmFpNHA14CTgF7gbkk3RMSD/Vb914iY2/YAzawhbsHR+Yo8ozgeWB8Rj0bEn4HrgfkFxmNmQ+AWHJ2vyEQxCfh91fvedF5/MyWtlvRzSUfX25ikBZJWSFqxadOmVsdqZnW4BUfnKzJRqMa86Pf+HuDwiHg18FVgSb2NRcRVEdETET0TJ05sXZRmlunU4yZxyenTmTShCwGTJnRxyenTPZDdQYq866kXOLTq/WRgQ/UKEfFM1fRNkr4u6YCIeKpNMZpZA9yCo7MNeEYh6Z2NzBuCu4EjJR0haW/gLOCGfvt5uSSl08en8T7dgn2bmVmDGrn0dEGD8wYlInYB5wNLgbXADyLiAUnnSTovXe0MYI2k1cDlwFkR0f/ylJmZ5ajupSdJpwBzgEmSLq9atB+wqxU7j4ibgJv6zbuyavoK4IpW7MvM6nNltWXJGqPYAKwA5gErq+ZvBT6WZ1Bm1j6urLaB1E0UEbEaWC1pMfBsRDwHuwvlXtKm+MwsZ1mV1U4UBo2NUSwDqitnuoBb8gnHzNrNldU2kEYSxdiI2FZ5k07vk19IZtZOrqy2gTSSKJ6VNKPyRtJrAP9Xw6xDuLLaBtJIwd1HgR9KqhTDHQycmVtEZtZWlXEI3/Vk9QyYKCLibkmvAqaStN14KCJ25h6ZmbWNK6stSyOV2fsAnwQ+EhH3A1Mkue23mdkI0cgYxbeAPwMz0/e9wP/MLSIzMyuVRsYoXhkRZ0p6N0BEbK/0XzKzcnBlteWpkUTxZ0ldpC3AJb0S+FOuUZlZw1xZbXlr5NLTp4GbgUMlXQfcCnwi16jMrGF+ZrXlLasp4KyIuB1YDpwOnEBy19NH/DwIs/JwZbXlLeuMotIx9s6IeDoifhYRP3WSMCsXV1Zb3rLGKHZK+hYwuV+bcQAi4sP5hWVmjVo4e+oeYxTgymprraxEMRc4Efhr9mwzbmYl4spqy1tWolgYEZ+UdFhEXNu2iMxs0FxZbXnKGqOYI2kMybOszcxshMo6o7gZeAoYJ+kZkjueovIzIvZrQ3xmZlawumcUEbEwIl4K/Cwi9ouI8dU/2xijmZkVqJHusfMlHQ4cGRG3pFXae0XE1vzDMxsZ3ILDyqyR7rEfAH4EfDOdNRlYkmNMZiNKpQVH3+btBC+04Fhyb1/RoZkBjbXw+BAwC3gGICIeBg7MMyizkcQtOKzsGkkUf4qIP1feSNqLtEGgmTXPLTis7BpJFP8i6UKgS9JJwA+BG1uxc0knS1onab2kRTWWS9Ll6fL7qp/dbdYp3ILDyq6RRLEI2ATcD/wdcBPwqWZ3LGk08DXgFGAa8G5J0/qtdgpwZPpaAHyj2f2alc3C2VPpGjN6j3luwWFl0shdT88D/5C+Wul4YH1EPAog6XpgPvBg1TrzgW9HRAC/kTRB0sER8USLYzErjFtwWNk18uCivEwCfl/1vhd4XQPrTAKcKKyjuAWHlVkjl57yUutxqv0HyRtZJ1lRWiBphaQVmzZtajo4MzNLFJkoeoFDq95PBjYMYR0AIuKqiOiJiJ6JEye2NFAzs5Es6wl3N5JxG2xEzGty33cDR0o6AugjaT74n/utcwNwfjp+8Tpgi8cnrIxcWW2dLGuM4ot57jgidkk6H1gKjAaujogHJJ2XLr+S5A6rOcB64I/Ae/OMyWwoKpXVlaK5SmU14GRhHUHJDUWdpaenJ1asWFF0GDZCzLr0l/TVKI6bNKGL2xf9dQERmQ2epJUR0VNr2YB3PUk6EriEpNZhbGV+RLyiZRGaDWOurLZO18hg9rdICt12AW8Gvg38U55BmQ0nrqy2TtdIouiKiFtJLlM9HhGfIXmOtpnhymrrfI0U3O2QNAp4OB187sPdY812c2W1dbpGEsVHgX2ADwOfIzmbODvHmMyGHVdWWydrpNfT3enkNknnAvtGxDP5hmVmZmXRyBPuvitpP0njSBr2rZO0MP/QzMysDBoZzJ6WnkGcSlIAdxjwt3kGZWZm5dHIGMUYSWNIEsUVEbFTUudV6dmI5hYcZvU1ckbxTeAxYBywXNLhpM/PNusElRYcfZu3E7zQgmPJvX1Fh2ZWCgMmioi4PCImRcScSDxOUnhn1hEuW7pud5+miu07n+OypesKisisXBpp4XFxnUWfbXEsZoVwCw6zbI1cenq26vUcyXOsp+QYk1lbuQWHWbZG6ii+VP1e0hdJnhNh1hEWzp66R5twcAsOs2pDeWb2PoA7x1rHcAsOs2yNjFHczwtPuhsNTMTjE9Zh3ILDrL5GzijmVk3vAp6MiF05xWNmZiXTyBjF4wCSDiR5cNEhkoiIf887ODMzK14jl57mAV8CDgE2AocDa4Gj8w3NrHGurDbLTyO3x34OOAH4bUQcAbwFuD3XqMwGwZXVZvlqJFHsjIingVGSRkXEbUB3vmGZNc6V1Wb5amQwe7OkfYHlwHWSNpIMapuVgiurzfLVyBnFfOCPwMeAm4FHgLfnGZTZYLiy2ixfdROFpL+QNCsino2I5yNiV0RcC6wCJrQrQLOBLJw9la4xo/eY58pqs9bJOqP4CrC1xvw/psvMSuHU4yZxyenTmTShCwGTJnRxyenTfdeTWYtkjVFMiYj7+s+MiBWSpjSzU0kvA75P0lzwMeBdEfGHGus9RpKsngN2RURPM/u1zuXKarP8ZJ1RjM1Y1uzF30XArRFxJHBr+r6eN0dEt5OEmVkxshLF3ZI+0H+mpHOBlU3udz5wbTp9LcljVs3MrISyLj19FFgs6T28kBh6gL2B05rc70ER8QRARDyRtgepJYBl6TO6vxkRV9XboKQFwAKAww47rMnwrN1cWW1WXnUTRUQ8Cbxe0puBY9LZP4uIXzayYUm3AC+vseiiQcQ3KyI2pInkF5IeiojldeK9CrgKoKenJ2qtY+VUqayuFM1VKqsBJwuzEmikKeBtwG2D3XBEnFhvmaQnJR2cnk0cTNJDqtY2NqQ/N0paDBxPUvhnHSSrstqJwqx4jRTc5eEG4Ox0+mzgJ/1XkDRO0vjKNPBWYE3bIrS2cWW1WbkVlSguBU6S9DBwUvoeSYdIuild5yDg15JWA3eRXPa6uZBoLVeurDYrt6E8CrVpaZPBt9SYvwGYk04/Cry6zaFZAfzMarNyKyRRmFXzM6vNys2JwkrBldVm5VXUGIWZmQ0TThRmZpbJicLMzDJ5jMJawi04zDqXE4U1zS04zDqbLz1Z07JacJjZ8OdEYU1zCw6zzuZEYU1zCw6zzuZEYU1bOHsqXWNG7zHPLTjMOocHs61pbsFh1tmcKKwl3ILDrHP50pOZmWVyojAzs0y+9GSAK6vNrD4nCnNltZll8qUnc2W1mWVyojBXVptZJicKc2W1mWVyojBXVptZJg9mmyurzSyTE4UBrqw2s/p86cnMzDIVkigkvVPSA5Kel9STsd7JktZJWi9pUTtjNDOzRFFnFGuA04Hl9VaQNBr4GnAKMA14t6Rp7QnPzMwqChmjiIi1AJKyVjseWB8Rj6brXg/MBx7MPcBhyC04zCwvZR6jmAT8vup9bzqvJkkLJK2QtGLTpk25B1cmlRYcfZu3E7zQgmPJvX1Fh2ZmHSC3RCHpFklrarzmN7qJGvOi3soRcVVE9EREz8SJE4cW9DDlFhxmlqfcLj1FxIlNbqIXOLTq/WRgQ5Pb7EhuwWFmeSrzpae7gSMlHSFpb+As4IaCYyolt+AwszwVdXvsaZJ6gZnAzyQtTecfIukmgIjYBZwPLAXWAj+IiAeKiLfs3ILDzPJU1F1Pi4HFNeZvAOZUvb8JuKmNoQ1LbsFhZnlyC48O4RYcZpaXMo9RmJlZCThRmJlZJl96KglXVptZWTlRlEClsrpSNFeprAacLMyscL70VAKurDazMnOiKAFXVptZmTlRlIArq82szJwoSsCV1WZWZh7MLgFXVptZmTlRlIQrq82srHzpyczMMjlRmJlZJicKMzPL5DGKFnELDjPrVE4ULeAWHGbWyXzpqQXcgsPMOpkTRQu4BYeZdTInihZwCw4z62ROFC3gFhxm1sk8mN0CbsFhZp3MiaJF3ILDzDqVLz2ZmVkmJwozM8tUyKUnSe8EPgMcBRwfESvqrPcYsBV4DtgVET15xeTKajOz2ooao1gDnA58s4F13xwRT+UZjCurzczqK+TSU0SsjYjSlC27strMrL6yj1EEsEzSSkkL8tqJK6vNzOrL7dKTpFuAl9dYdFFE/KTBzcyKiA2SDgR+IemhiFheZ38LgAUAhx122KBiPWRCF301koIrq83McjyjiIgTI+KYGq9GkwQRsSH9uRFYDByfse5VEdETET0TJ04cVKyurDYzq6+0l54kjZM0vjINvJVkELzlTj1uEpecPp1JE7oQMGlCF5ecPt0D2WZmFHd77GnAV4GJwM8krYqI2ZIOAf4xIuYABwGLJVXi/G5E3JxXTK6sNjOrrZBEERGLSS4l9Z+/AZiTTj8KvLrNoZmZWT+lvfRkZmbl4ERhZmaZnCjMzCyTE4WZmWVSRBQdQ8tJ2gQ8PsRfPwDItbfUEDmuwXFcg+O4BqcT4zo8ImoWoXVkomiGpBV5dqkdKsc1OI5rcBzX4Iy0uHzpyczMMjlRmJlZJieKF7uq6ADqcFyD47gGx3ENzoiKy2MUZmaWyWcUZmaWyYnCzMwyjfhEIekySQ9Juk/SYkkT6qx3sqR1ktZLWtSGuN4p6QFJz0uqe7ubpMck3S9plaQVJYqr3cfrZZJ+Ienh9Od/qrNeW47XQJ9ficvT5fdJmpFXLIOM602StqTHZ5Wki9sQ09WSNkqq+RiBAo/VQHG1/Vil+z1U0m2S1qb/Fj9SY53WHrOIGNEvkudc7JVOfwH4Qo11RgOPAK8A9gZWA9NyjusoYCrwK6AnY73HgAPaeLwGjKug4/W/gUXp9KJaf47tOl6NfH6SLsk/BwScAPxbG/7sGonrTcBP2/X3Kd3nG4EZwJo6y9t+rBqMq+3HKt3vwcCMdHo88Nu8/36N+DOKiFgWEbvSt78BJtdY7XhgfUQ8GhF/Bq4H5ucc19qIWJfnPoaiwbjafrzS7V+bTl8LnJrz/rI08vnnA9+OxG+ACZIOLkFcbRfJ443/I2OVIo5VI3EVIiKeiIh70umtwFqg/8N0WnrMRnyi6Od9JFm4v0nA76ve9/LiP5iiBLBM0sr0ueFlUMTxOiginoDkHxJwYJ312nG8Gvn8RRyjRvc5U9JqST+XdHTOMTWizP/+Cj1WkqYAxwH/1m9RS49ZIQ8uajdJtwAvr7Hookif4S3pImAXcF2tTdSY1/R9xY3E1YBZEbFB0oHALyQ9lP5PqMi42n68BrGZlh+vGhr5/LkcowE0ss97SHr+bJM0B1gCHJlzXAMp4lg1otBjJWlf4J+Bj0bEM/0X1/iVIR+zEZEoIuLErOWSzgbmAm+J9AJfP73AoVXvJwMb8o6rwW1sSH9ulLSY5PJCU198LYir7cdL0pOSDo6IJ9JT7I11ttHy41VDI58/l2PUbFzVXzgRcZOkr0s6ICKKbIBXxLEaUJHHStIYkiRxXUT8uMYqLT1mI/7Sk6STgU8C8yLij3VWuxs4UtIRkvYGzgJuaFeM9UgaJ2l8ZZpkYL7mHRptVsTxugE4O50+G3jRmU8bj1cjn/8G4L+kd6ecAGypXDrL0YBxSXq5lDyoXtLxJN8RT+cc10CKOFYDKupYpfv8f8DaiPg/dVZr7TFr94h92V7AepJreavS15Xp/EOAm6rWm0Nyd8EjJJdg8o7rNJL/FfwJeBJY2j8ukrtXVqevB8oSV0HHa3/gVuDh9OfLijxetT4/cB5wXjot4Gvp8vvJuLOtzXGdnx6b1SQ3d7y+DTF9D3gC2Jn+3Tq3JMdqoLjafqzS/f4VyWWk+6q+t+bkeczcwsPMzDKN+EtPZmaWzYnCzMwyOVGYmVkmJwozM8vkRGFmZpmcKCw3krY1sM5HJe3Tjnjq7L87raqtvJ+nFna7lfRNSbP6zbtc0t9Xvb9I0tdasK9fKaOjr9lQ+fZYy42kbRGx7wDrPEZyj3fD1aySRkfEc83Gl27rnHT/57diezW2vwp4TXW8kvYjuff9RJL74X8JHBcRm5vc16+Aj0dEy9qnS9orXmiaaSOUzygsd2nf/l9J+pGSZ39cl1aMfpikIO42Sbel675V0p2S7pH0w7SfTeU5EhdL+jXwCUl3VW1/iqT70unXSPqXtOnf0rSdR+V/21+QdJek30p6Q1qd/FngTCXPEzhT0jmSrkh/53BJtyrp53+rpMPS+dekZwV3SHpU0hl1PvdRwG/7J7VIWj9cBFxBUhR1cf8kIekUST/odwxvTKe/IWmFkmcR/I86+95WNX2GpGvS6YmS/lnS3elrVo3fPSc99jeSNFDcN/389yh5lsf8quO+VtI/pLEsk9SVLnttetzuVPLMlzXp/NHp+7vT5X9XK34rmXZUEvo1Ml/AtvTnm4AtJP1mRgF3An+VLnuM9PkQwAEkfZfGpe8/SfIlWlnvE1XbXgW8omq9TwFjgDuAien8M4Gr0+lfAV9Kp+cAt6TT5wBXVG1393vgRuDsdPp9wJJ0+hrgh+lnmUbSurvW5//vwPsyjs+dwK/rLNsL+PeqY/EN4G/S6UrV+ej0cx1b9Rl7qo99On0GcE06/d2qY38YSRuI/vs+h6QS+WVVsexX9We0nqTydwpJI83udNkPqmJcQ1qpDFxK+kwHYAHwqXT6JcAK4Iii/676lf0aEU0BrRTuiohe2H05Zgrw637rnEDyxXt72kJnb5Iv04rvV03/AHgXyZfQmelrKnAMSVdYSL5Iq/vbVJqnrUz3P5CZwOnp9D+RPBypYklEPA88KOmgOr8/G3hvrQWSJpN0wg1J+0bEHuM5EbFL0s3A2yX9CHgb8Il08buUtEjfi+QhNtNI2jk04kRgWnp8APaTND6S5xpU+0VEVJ7FIOB/SXoj8DxJu+rKZ/5dRKxKp1cCU5Q8JXJ8RNyRzv8uSdNNSPprHVt1FvZSko6rv2swfiuAE4W1y5+qpp+j9t89kXxBvbvONp6tmv4+8ENJPwYiIh6WNB14ICJmDhBDvf0PpHpAr/rzvKilczpAPyHSbrU1/F/gMyRPDPw0sFDS50kSAhHRTfIZP0Ty8Jy7I2KrpCOAjwOvjYg/pJeUxg4Qa/XyUcDMiNheJ66K6mP9HmAiyVjLznRcqbLN/n+uXdRucV0h4L9FxNIB9m8l4jEKK9pWksc5QtJYbZakv4Dky1bSX9b6pYh4hOSL6e954UxjHTBR0sz098do4IfJVO+/vztIOqxC8mXZ/wwoy5uB22otkHQKyYOVvg18DjhN0rSIuCgiutMkAcmlpBnAB3jhM+5H8iW+JT2TOaXO/p+UdJSkUSSNHCuWkTSzq8TS3f8Xa3gpsDFNEm8GDs9aOSL+AGxV0rUUXjiGAEuB/6qkTTaS/lJJJ18rMScKK9pVwM8l3RYRm0iuj38vHZz+DfCqjN/9PvA3JJehiOTxnmcAX5C0mmQc4/UD7P82kksxqySd2W/Zh4H3prH8LfCih9hnOAW4uf9MSWOBrwAfjMSzJJeUrui/biSD4D9Nt/XTdN5q4F6SrqVXA7fX2f+i9Hd+yZ6X3z4M9KQDyQ+SdBwdyHXp76wgSZgPNfA75wJXSbqT5CxiSzr/H4EHgXvSAe5v4isbpefbY81yIOke4HURsbPoWIpQPe6ipC7l4IgYTKK1EnGiMLOWS8/OLiA5W3gcOCc9Y7RhyInCzMwyeYzCzMwyOVGYmVkmJwozM8vkRGFmZpmcKMzMLNP/B/r31c5xS9pvAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] @@ -1186,7 +1182,7 @@ }, { "cell_type": "code", - "execution_count": 49, + "execution_count": 23, "metadata": {}, "outputs": [], "source": [ @@ -1223,7 +1219,7 @@ }, { "cell_type": "code", - "execution_count": 50, + "execution_count": 24, "metadata": {}, "outputs": [ { @@ -1252,7 +1248,7 @@ }, { "cell_type": "code", - "execution_count": 51, + "execution_count": 25, "metadata": {}, "outputs": [], "source": [ @@ -1286,7 +1282,7 @@ }, { "cell_type": "code", - "execution_count": 52, + "execution_count": 26, "metadata": {}, "outputs": [ { @@ -1346,7 +1342,7 @@ }, { "cell_type": "code", - "execution_count": 55, + "execution_count": 27, "metadata": {}, "outputs": [], "source": [ @@ -1362,7 +1358,7 @@ }, { "cell_type": "code", - "execution_count": 56, + "execution_count": 28, "metadata": {}, "outputs": [ { @@ -1426,16 +1422,16 @@ }, { "cell_type": "code", - "execution_count": 57, + "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 57, + "execution_count": 29, "metadata": {}, "output_type": "execute_result" } @@ -1454,15 +1450,15 @@ }, { "cell_type": "code", - "execution_count": 58, + "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Causal effect for S = -1.00 is -0.71\n", - "Causal effect for S = 1.00 is 0.72\n" + "Causal effect for S = -1.00 is -0.65\n", + "Causal effect for S = 1.00 is 0.69\n" ] } ], @@ -1506,7 +1502,7 @@ }, { "cell_type": "code", - "execution_count": 59, + "execution_count": 31, "metadata": {}, "outputs": [ { @@ -1592,7 +1588,7 @@ }, { "cell_type": "code", - "execution_count": 60, + "execution_count": 32, "metadata": {}, "outputs": [], "source": [ @@ -1617,16 +1613,16 @@ }, { "cell_type": "code", - "execution_count": 61, + "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 61, + "execution_count": 33, "metadata": {}, "output_type": "execute_result" } @@ -1642,7 +1638,7 @@ }, { "cell_type": "code", - "execution_count": 62, + "execution_count": 34, "metadata": {}, "outputs": [ { @@ -1682,7 +1678,7 @@ }, { "cell_type": "code", - "execution_count": 65, + "execution_count": 35, "metadata": {}, "outputs": [ { @@ -1764,7 +1760,7 @@ }, { "cell_type": "code", - "execution_count": 66, + "execution_count": 36, "metadata": {}, "outputs": [], "source": [ @@ -1799,7 +1795,7 @@ }, { "cell_type": "code", - "execution_count": 67, + "execution_count": 37, "metadata": {}, "outputs": [ { @@ -1838,7 +1834,7 @@ }, { "cell_type": "code", - "execution_count": 68, + "execution_count": 38, "metadata": {}, "outputs": [ { @@ -1876,7 +1872,7 @@ }, { "cell_type": "code", - "execution_count": 70, + "execution_count": 39, "metadata": {}, "outputs": [ {