diff --git a/econml/__init__.py b/econml/__init__.py index 13b63dd15..99e8c4060 100644 --- a/econml/__init__.py +++ b/econml/__init__.py @@ -21,6 +21,7 @@ 'dowhy', 'utilities', 'federated_learning', + 'validate', '__version__'] from ._version import __version__ diff --git a/econml/iv/dr/_dr.py b/econml/iv/dr/_dr.py index c06df6278..5c1df1c94 100644 --- a/econml/iv/dr/_dr.py +++ b/econml/iv/dr/_dr.py @@ -2186,10 +2186,10 @@ def true_heterogeneity_function(X): est.fit(Y=y, T=T, Z=Z, X=X) >>> est.effect(X[:3]) - array([-1.74672..., 1.57225..., -1.58916...]) + array([-1.74672..., 1.57..., -1.58916...]) >>> est.effect_interval(X[:3]) - (array([-7.05230..., -6.78656..., -5.11344...]), - array([3.55885..., 9.93108..., 1.93512...])) + (array([-7.05230..., -6..., -5.11344...]), + array([3.55885..., 9.9..., 1.93512...])) """ def __init__(self, *, diff --git a/econml/tests/test_drtester.py b/econml/tests/test_drtester.py new file mode 100644 index 000000000..9e890c190 --- /dev/null +++ b/econml/tests/test_drtester.py @@ -0,0 +1,263 @@ +import unittest + +import numpy as np +import pandas as pd +import scipy.stats as st +from sklearn.ensemble import RandomForestClassifier, GradientBoostingRegressor + +from econml.validate.drtester import DRtester +from econml.dml import DML + + +class TestDRTester(unittest.TestCase): + + @staticmethod + def _get_data(num_treatments=1): + np.random.seed(576) + + N = 20000 # number of units + K = 5 # number of covariates + + # Generate random Xs + X_mu = np.zeros(5) # Means of Xs + # Random covariance matrix of Xs + X_sig = np.diag(np.random.rand(5)) + X = st.multivariate_normal(X_mu, X_sig).rvs(N) + + # Effect of Xs on outcome + X_beta = np.random.uniform(0, 5, K) + # Effect of treatment on outcomes + D_beta = np.arange(num_treatments + 1) + # Effect of treatment on outcome conditional on X1 + DX1_beta = np.array([0] * num_treatments + [3]) + + # Generate treatments based on X and random noise + beta_treat = np.random.uniform(-1, 1, (num_treatments + 1, K)) + D1 = np.zeros((N, num_treatments + 1)) + for k in range(num_treatments + 1): + D1[:, k] = X @ beta_treat[k, :] + np.random.gumbel(0, 1, N) + D = np.array([np.where(D1[i, :] == np.max(D1[i, :]))[0][0] for i in range(N)]) + D_dum = pd.get_dummies(D) + + # Generate Y (based on X, D, and random noise) + Y_sig = 1 # Variance of random outcome noise + Y = X @ X_beta + (D_dum @ D_beta) + X[:, 1] * (D_dum @ DX1_beta) + np.random.normal(0, Y_sig, N) + Y = Y.to_numpy() + + train_prop = .5 + train_N = np.ceil(train_prop * N) + ind = np.array(range(N)) + train_ind = np.random.choice(N, int(train_N), replace=False) + val_ind = ind[~np.isin(ind, train_ind)] + + Xtrain, Dtrain, Ytrain = X[train_ind], D[train_ind], Y[train_ind] + Xval, Dval, Yval = X[val_ind], D[val_ind], Y[val_ind] + + return Xtrain, Dtrain, Ytrain, Xval, Dval, Yval + + def test_multi(self): + Xtrain, Dtrain, Ytrain, Xval, Dval, Yval = self._get_data(num_treatments=2) + + # Simple classifier and regressor for propensity, outcome, and cate + reg_t = RandomForestClassifier(random_state=0) + reg_y = GradientBoostingRegressor(random_state=0) + + cate = DML( + model_y=reg_y, + model_t=reg_t, + model_final=reg_y + ).fit(Y=Ytrain, T=Dtrain, X=Xtrain) + + # test the DR outcome difference + my_dr_tester = DRtester( + model_regression=reg_y, + model_propensity=reg_t, + cate=cate + ).fit_nuisance( + Xval, Dval, Yval, Xtrain, Dtrain, Ytrain + ) + dr_outcomes = my_dr_tester.dr_val_ + + ates = dr_outcomes.mean(axis=0) + for k in range(dr_outcomes.shape[1]): + ate_errs = np.sqrt(((dr_outcomes[:, k] - ates[k]) ** 2).sum() / + (dr_outcomes.shape[0] * (dr_outcomes.shape[0] - 1))) + + self.assertLess(abs(ates[k] - (k + 1)), 2 * ate_errs) + + res = my_dr_tester.evaluate_all(Xval, Xtrain) + res_df = res.summary() + + for k in range(3): + if k == 0: + with self.assertRaises(Exception) as exc: + res.plot_cal(k) + self.assertTrue(str(exc.exception) == 'Plotting only supported for treated units (not controls)') + else: + self.assertTrue(res.plot_cal(k) is not None) + + self.assertGreater(res_df.blp_pval.values[0], 0.1) # no heterogeneity + self.assertLess(res_df.blp_pval.values[1], 0.05) # heterogeneity + + self.assertLess(res_df.cal_r_squared.values[0], 0) # poor R2 + self.assertGreater(res_df.cal_r_squared.values[1], 0) # good R2 + + self.assertLess(res_df.qini_pval.values[1], res_df.qini_pval.values[0]) + + def test_binary(self): + Xtrain, Dtrain, Ytrain, Xval, Dval, Yval = self._get_data(num_treatments=1) + + # Simple classifier and regressor for propensity, outcome, and cate + reg_t = RandomForestClassifier(random_state=0) + reg_y = GradientBoostingRegressor(random_state=0) + + cate = DML( + model_y=reg_y, + model_t=reg_t, + model_final=reg_y + ).fit(Y=Ytrain, T=Dtrain, X=Xtrain) + + # test the DR outcome difference + my_dr_tester = DRtester( + model_regression=reg_y, + model_propensity=reg_t, + cate=cate + ).fit_nuisance( + Xval, Dval, Yval, Xtrain, Dtrain, Ytrain + ) + dr_outcomes = my_dr_tester.dr_val_ + + ate = dr_outcomes.mean(axis=0) + ate_err = np.sqrt(((dr_outcomes - ate) ** 2).sum() / + (dr_outcomes.shape[0] * (dr_outcomes.shape[0] - 1))) + truth = 1 + self.assertLess(abs(ate - truth), 2 * ate_err) + + res = my_dr_tester.evaluate_all(Xval, Xtrain) + res_df = res.summary() + + for k in range(2): + if k == 0: + with self.assertRaises(Exception) as exc: + res.plot_cal(k) + self.assertTrue(str(exc.exception) == 'Plotting only supported for treated units (not controls)') + else: + self.assertTrue(res.plot_cal(k) is not None) + + self.assertLess(res_df.blp_pval.values[0], 0.05) # heterogeneity + self.assertGreater(res_df.cal_r_squared.values[0], 0) # good R2 + self.assertLess(res_df.qini_pval.values[0], 0.05) # heterogeneity + + def test_nuisance_val_fit(self): + Xtrain, Dtrain, Ytrain, Xval, Dval, Yval = self._get_data(num_treatments=1) + + # Simple classifier and regressor for propensity, outcome, and cate + reg_t = RandomForestClassifier(random_state=0) + reg_y = GradientBoostingRegressor(random_state=0) + + cate = DML( + model_y=reg_y, + model_t=reg_t, + model_final=reg_y + ).fit(Y=Ytrain, T=Dtrain, X=Xtrain) + + # test the DR outcome difference + my_dr_tester = DRtester( + model_regression=reg_y, + model_propensity=reg_t, + cate=cate + ).fit_nuisance(Xval, Dval, Yval) + + dr_outcomes = my_dr_tester.dr_val_ + + ate = dr_outcomes.mean(axis=0) + ate_err = np.sqrt(((dr_outcomes - ate) ** 2).sum() / + (dr_outcomes.shape[0] * (dr_outcomes.shape[0] - 1))) + truth = 1 + self.assertLess(abs(ate - truth), 2 * ate_err) + + # use evaluate_blp to fit on validation only + blp_res = my_dr_tester.evaluate_blp(Xval) + + self.assertLess(blp_res.pvals[0], 0.05) # heterogeneity + + for kwargs in [{}, {'Xval': Xval}]: + with self.assertRaises(Exception) as exc: + my_dr_tester.evaluate_cal(kwargs) + self.assertTrue( + str(exc.exception) == "Must fit nuisance models on training sample data to use calibration test" + ) + + def test_exceptions(self): + Xtrain, Dtrain, Ytrain, Xval, Dval, Yval = self._get_data(num_treatments=1) + + # Simple classifier and regressor for propensity, outcome, and cate + reg_t = RandomForestClassifier(random_state=0) + reg_y = GradientBoostingRegressor(random_state=0) + + cate = DML( + model_y=reg_y, + model_t=reg_t, + model_final=reg_y + ).fit(Y=Ytrain, T=Dtrain, X=Xtrain) + + # test the DR outcome difference + my_dr_tester = DRtester( + model_regression=reg_y, + model_propensity=reg_t, + cate=cate + ) + + # fit nothing + for func in [my_dr_tester.evaluate_blp, my_dr_tester.evaluate_cal, my_dr_tester.evaluate_qini]: + with self.assertRaises(Exception) as exc: + func() + if func.__name__ == 'evaluate_cal': + self.assertTrue( + str(exc.exception) == "Must fit nuisance models on training sample data to use calibration test" + ) + else: + self.assertTrue(str(exc.exception) == "Must fit nuisances before evaluating") + + my_dr_tester = my_dr_tester.fit_nuisance( + Xval, Dval, Yval, Xtrain, Dtrain, Ytrain + ) + + for func in [ + my_dr_tester.evaluate_blp, + my_dr_tester.evaluate_cal, + my_dr_tester.evaluate_qini, + my_dr_tester.evaluate_all + ]: + with self.assertRaises(Exception) as exc: + func() + if func.__name__ == 'evaluate_blp': + self.assertTrue( + str(exc.exception) == "CATE predictions not yet calculated - must provide Xval" + ) + else: + self.assertTrue(str(exc.exception) == + "CATE predictions not yet calculated - must provide both Xval, Xtrain") + + for func in [ + my_dr_tester.evaluate_cal, + my_dr_tester.evaluate_qini, + my_dr_tester.evaluate_all + ]: + with self.assertRaises(Exception) as exc: + func(Xval=Xval) + self.assertTrue( + str(exc.exception) == "CATE predictions not yet calculated - must provide both Xval, Xtrain") + + cal_res = my_dr_tester.evaluate_cal(Xval, Xtrain) + self.assertGreater(cal_res.cal_r_squared[0], 0) # good R2 + + my_dr_tester = DRtester( + model_regression=reg_y, + model_propensity=reg_t, + cate=cate + ).fit_nuisance( + Xval, Dval, Yval, Xtrain, Dtrain, Ytrain + ) + qini_res = my_dr_tester.evaluate_qini(Xval, Xtrain) + self.assertLess(qini_res.pvals[0], 0.05) diff --git a/econml/validate/__init__.py b/econml/validate/__init__.py new file mode 100644 index 000000000..a1c2adad7 --- /dev/null +++ b/econml/validate/__init__.py @@ -0,0 +1,11 @@ +# Copyright (c) PyWhy contributors. All rights reserved. +# Licensed under the MIT License. + +""" +A suite of validation methods for CATE models. +""" + +from .drtester import DRtester + + +__all__ = ['DRtester'] diff --git a/econml/validate/drtester.py b/econml/validate/drtester.py new file mode 100644 index 000000000..1948bed80 --- /dev/null +++ b/econml/validate/drtester.py @@ -0,0 +1,594 @@ +from typing import Tuple, Union, List + +import numpy as np +import pandas as pd +import scipy.stats as st +from sklearn.model_selection import check_cv +from sklearn.model_selection import cross_val_predict, StratifiedKFold, KFold +from statsmodels.api import OLS +from statsmodels.tools import add_constant + +from .results import CalibrationEvaluationResults, BLPEvaluationResults, QiniEvaluationResults, EvaluationResults +from .utils import calculate_dr_outcomes, calc_qini_coeff + + +class DRtester: + + """ + Validation tests for CATE models. Includes the best linear predictor (BLP) test as in Chernozhukov et al. (2022), + the calibration test in Dwivedi et al. (2020), and the QINI coefficient as in Radcliffe (2007). + + **Best Linear Predictor (BLP)** + + Runs ordinary least squares (OLS) of doubly robust (DR) outcomes on DR outcome predictions from the CATE model + (and a constant). If the CATE model captures true heterogeneity, then the OLS estimate on the CATE predictions + should be positive and significantly different than 0. + + **Calibration** + + First, units are binned based on out-of-sample defined quantiles of the CATE predictions (s(Z)). Within each bin + (k), the absolute difference between the mean CATE prediction and DR outcome is calculated, along with the + absolute difference in the mean CATE prediction and the overall ATE. These measures are then summed across bins, + weighted by a probability a unit is in each bin. + + .. math:: + + \\mathrm{Cal}_G = \\sum_k \\pi(k) |{\\E[s(Z) | k] - \\E[Y^{DR} | k]}| + + \\mathrm{Cal}_O = \\sum_k \\pi(k) |{\\E[s(Z) | k] - \\E[Y^{DR}]}| + + The calibration r-squared is then defined as + + .. math:: + + \\mathcal{R^2}_C = 1 - \\frac{\\mathrm{Cal}_G}{\\mathrm{Cal}_O} + + The calibration r-squared metric is similar to the standard R-square score in that it can take any value + less than or equal to 1, with scores closer to 1 indicating a better calibrated CATE model. + + **QINI** + + Units are ordered by predicted CATE values and a running measure of the average treatment effect in each cohort is + kept as we progress through ranks. The QINI coefficient is then the area under the resulting curve, with a value + of 0 interpreted as corresponding to a model with randomly assigned CATE coefficients. All calculations are + performed on validation dataset results, using the training set as input. + + More formally, the QINI curve is given by the following function: + + .. math:: + + \\tau_{QINI}(q) = \\mathrm{Cov}(Y^{DR}(g,p), \\mathbb{1}\\{\\hat{\\tau}(Z) \\geq \\hat{\\mu}(q)\\}) + + Where :math:`q` is the desired quantile, :math:`\\hat{\\mu}` is the quantile function, and :math:`\\hat{\\tau}` is + the predicted CATE function. + :math:`Y^{DR}(g,p)` refers to the doubly robust outcome difference (relative to control) for the given observation. + + The QINI coefficient is then given by: + + .. math:: + + QINI = \\int_0^1 \\tau_{QINI}(q) dq + + Parameters + ---------- + model_regression: estimator + Nuisance model estimator used to fit the outcome to features. Must be able to implement `fit' and `predict' + methods + + model_propensity: estimator + Nuisance model estimator used to fit the treatment assignment to features. Must be able to implement `fit' + method and either `predict' (in the case of binary treatment) or `predict_proba' methods (in the case of + multiple categorical treatments). + + n_splits: integer, default 5 + Number of splits used to generate cross-validated predictions + + References + ---------- + + + [Chernozhukov2022] V. Chernozhukov et al. + Generic Machine Learning Inference on Heterogeneous Treatment Effects in Randomized Experiments + arXiv preprint arXiv:1712.04802, 2022. + ``_ + + [Dwivedi2020] R. Dwivedi et al. + Stable Discovery of Interpretable Subgroups via Calibration in Causal Studies + arXiv preprint arXiv:2008.10109, 2020. + ``_ + + + [Radcliffe2007] N. Radcliffe + Using control groups to target on predicted lift: Building and assessing uplift model. + Direct Marketing Analytics Journal (2007), pages 14–21. + """ + + def __init__( + self, + *, + model_regression, + model_propensity, + cate, + cv: Union[int, List] = 5 + ): + self.model_regression = model_regression + self.model_propensity = model_propensity + self.cate = cate + self.cv = cv + + def get_cv_splitter(self, random_state: int = 123): + """ + Generates splitter object for cross-validation. + Checks if the cv object passed at initialization is a splitting mechanism or a number of folds and returns + appropriately modified object for use in downstream cross-validation. + + Parameters + ---------- + random_state: integer + seed for splitter, default is 123 + + Returns + ------ + None, but adds attribute 'cv_splitter' containing the constructed splitter object. + """ + splitter = check_cv(self.cv, [0], classifier=True) + if splitter != self.cv and isinstance(splitter, (KFold, StratifiedKFold)): + splitter.shuffle = True + splitter.random_state = random_state + self.cv_splitter = splitter + + def get_cv_splits(self, vars: np.array, T: np.array): + """ + Generates splits for cross-validation, given a set of variables and treatment vector. + + Parameters + ---------- + vars: (n x k) matrix or vector of length n + Features used in nuisance models + T: vector of length n_val + Treatment assignment vector. Control status must be minimum value. It is recommended to have + the control status be equal to 0, and all other treatments integers starting at 1. + + Returns + ------ + list of folds of the data, on which to run cross-validation + """ + if not hasattr(self, 'cv_splitter'): + self.get_cv_splitter() + + all_vars = [var if np.ndim(var) == 2 else var.reshape(-1, 1) for var in vars if var is not None] + to_split = np.hstack(all_vars) + folds = self.cv_splitter.split(to_split, T) + + return list(folds) + + def fit_nuisance( + self, + Xval: np.array, + Dval: np.array, + yval: np.array, + Xtrain: np.array = None, + Dtrain: np.array = None, + ytrain: np.array = None, + ): + """ + Generates nuisance predictions and calculates doubly robust (DR) outcomes either by (1) cross-fitting in the + validation sample, or (2) fitting in the training sample and applying to the validation sample. If Xtrain, + Dtrain, and ytrain are all not None, then option (2) will be implemented, otherwise, option (1) will be + implemented. In order to use the `evaluate_cal' method then Xtrain, Dtrain, and ytrain must all be specified. + + Parameters + ---------- + Xval: (n_val x k) matrix or vector of length n + Features used in nuisance models for validation sample + Dval: vector of length n_val + Treatment assignment of validation sample. Control status must be minimum value. It is recommended to have + the control status be equal to 0, and all other treatments integers starting at 1. + yval: vector of length n_val + Outcomes for the validation sample + Xtrain: (n_train x k) matrix or vector of length n, default ``None`` + Features used in nuisance models for training sample + Dtrain: vector of length n_train, default ``None'' + Treatment assignment of training sample. Control status must be minimum value. It is recommended to have + the control status be equal to 0, and all other treatments integers starting at 1. + ytrain: vector of length n_train, defaul ``None`` + Outcomes for the training sample + + Returns + ------ + self, with added attributes for the validation treatments (Dval), treatment values (tmts), + number of treatments excluding control (n_treat), boolean flag for whether training data is provided + (fit_on_train), doubly robust outcome values for the validation set (dr_val), and the DR ATE value (ate_val). + If training data is provided, also adds attributes for the doubly robust outcomes for the training + set (dr_train) and the training treatments (Dtrain) + """ + + self.Dval = Dval + + # Unique treatments (ordered, includes control) + self.treatments = np.sort(np.unique(Dval)) + + # Number of treatments (excluding control) + self.n_treat = len(self.treatments) - 1 + + # Indicator for whether + self.fit_on_train = (Xtrain is not None) and (Dtrain is not None) and (ytrain is not None) + + if self.fit_on_train: + # Get DR outcomes in training sample + reg_preds_train, prop_preds_train = self.fit_nuisance_cv(Xtrain, Dtrain, ytrain) + self.dr_train_ = calculate_dr_outcomes(Dtrain, ytrain, reg_preds_train, prop_preds_train) + + # Get DR outcomes in validation sample + reg_preds_val, prop_preds_val = self.fit_nuisance_train(Xtrain, Dtrain, ytrain, Xval) + self.dr_val_ = calculate_dr_outcomes(Dval, yval, reg_preds_val, prop_preds_val) + else: + # Get DR outcomes in validation sample + reg_preds_val, prop_preds_val = self.fit_nuisance_cv(Xval, Dval, yval) + self.dr_val_ = calculate_dr_outcomes(Dval, yval, reg_preds_val, prop_preds_val) + + # Calculate ATE in the validation sample + self.ate_val = self.dr_val_.mean(axis=0) + + return self + + def fit_nuisance_train( + self, + Xtrain: np.array, + Dtrain: np.array, + ytrain: np.array, + Xval: np.array + ) -> Tuple[np.array, np.array]: + """ + Fits nuisance models in training sample and applies to generate predictions in validation sample. + + Parameters + ---------- + Xtrain: (n_train x k) matrix + Training sample features used to predict both treatment status and outcomes + Dtrain: array of length n_train + Training sample treatment assignments. Should have integer values with the lowest-value corresponding to + the control treatment. It is recommended to have the control take value 0 and all other treatments be + integers starting at 1 + ytrain: array of length n_train + Outcomes for training sample + Xval: (n_train x k) matrix + Validation sample features used to predict both treatment status and outcomes + + Returns + ------- + 2 (n_val x n_treatment + 1) arrays corresponding to the predicted outcomes under treatment status and predicted + treatment probabilities, respectively. Both evaluated on validation set. + """ + + # Fit propensity in treatment + model_propensity_fitted = self.model_propensity.fit(Xtrain, Dtrain) + # Predict propensity scores + prop_preds = model_propensity_fitted.predict_proba(Xval) + + # Possible treatments (need to allow more than 2) + n = Xval.shape[0] + reg_preds = np.zeros((n, self.n_treat + 1)) + for i in range(self.n_treat + 1): + model_regression_fitted = self.model_regression.fit( + Xtrain[Dtrain == self.treatments[i]], ytrain[Dtrain == self.treatments[i]]) + reg_preds[:, i] = model_regression_fitted.predict(Xval) + + return reg_preds, prop_preds + + def fit_nuisance_cv( + self, + X: np.array, + D: np.array, + y: np.array + ) -> Tuple[np.array, np.array]: + """ + Generates nuisance function predictions using k-folds cross validation. + + Parameters + ---------- + X: (n x k) matrix + Features used to predict treatment/outcome + D: array of length n + Treatment assignments. Should have integer values with the lowest-value corresponding to the + control treatment. It is recommended to have the control take value 0 and all other treatments be integers + starting at 1 + y: array of length n + Outcomes + + Returns + ------- + 2 (n x n_treatment + 1) arrays corresponding to the predicted outcomes under treatment status and predicted + treatment probabilities, respectively. + """ + + splits = self.get_cv_splits([X], D) + prop_preds = cross_val_predict(self.model_propensity, X, D, cv=splits, method='predict_proba') + + # Predict outcomes + # T-learner logic + N = X.shape[0] + reg_preds = np.zeros((N, self.n_treat + 1)) + for k in range(self.n_treat + 1): + for train, test in splits: + model_regression_fitted = self.model_regression.fit(X[train][D[train] == self.treatments[k]], + y[train][D[train] == self.treatments[k]]) + reg_preds[test, k] = model_regression_fitted.predict(X[test]) + + return reg_preds, prop_preds + + def get_cate_preds( + self, + Xval: np.array, + Xtrain: np.array = None + ): + """ + Generates predictions from fitted CATE model. If Xtrain is None, then the predictions are generated using + k-folds cross-validation on the validation set. If Xtrain is specified, then the CATE is assumed to have + been fitted on the training sample (where the DR outcomes were generated using k-folds CV), and then applied + to the validation sample. + + Parameters + ---------- + Xval: (n_val x n_treatment) matrix + Validation set features to be used to predict (and potentially fit) DR outcomes in CATE model + Xtrain (n_train x n_treatment) matrix, defaul ``None`` + Training set features used to fit CATE model + + Returns + ------- + None, but adds attribute cate_preds_val_ for predicted CATE values on the validation set and, if training + data is provided, attribute cate_preds_train_ for predicted CATE values on the training set + """ + base = self.treatments[0] + vals = [self.cate.effect(X=Xval, T0=base, T1=t) for t in self.treatments[1:]] + self.cate_preds_val_ = np.stack(vals).T + + if Xtrain is not None: + trains = [self.cate.effect(X=Xtrain, T0=base, T1=t) for t in self.treatments[1:]] + self.cate_preds_train_ = np.stack(trains).T + + def evaluate_cal( + self, + Xval: np.array = None, + Xtrain: np.array = None, + n_groups: int = 4 + ) -> CalibrationEvaluationResults: + """ + Implements calibration test as in [Dwivedi2020] + + Parameters + ---------- + Xval: (n_val x n_treatment) matrix, default ``None`` + Validation sample features for CATE model. If not specified, then `fit_cate' method must already have been + implemented + Xtrain: (n_train x n_treatment) matrix, default ``None`` + Training sample features for CATE model. If not specified, then `fit cate' method must already have been + implemented (with Xtrain specified) + n_groups: integer, default 4 + Number of quantile-based groups used to calculate calibration score. + + Returns + ------- + CalibrationEvaluationResults object showing the results of the calibration test + """ + if not hasattr(self, 'dr_train_'): + raise Exception("Must fit nuisance models on training sample data to use calibration test") + + # if CATE is given explicitly or has not been fitted at all previously, fit it now + if (not hasattr(self, 'cate_preds_train_')) or (not hasattr(self, 'cate_preds_val_')): + if (Xval is None) or (Xtrain is None): + raise Exception('CATE predictions not yet calculated - must provide both Xval, Xtrain') + self.get_cate_preds(Xval, Xtrain) + + cal_r_squared = np.zeros(self.n_treat) + df_plot = pd.DataFrame() + for k in range(self.n_treat): + cuts = np.quantile(self.cate_preds_train_[:, k], np.linspace(0, 1, n_groups + 1)) + probs = np.zeros(n_groups) + g_cate = np.zeros(n_groups) + se_g_cate = np.zeros(n_groups) + gate = np.zeros(n_groups) + se_gate = np.zeros(n_groups) + for i in range(n_groups): + # Assign units in validation set to groups + ind = (self.cate_preds_val_[:, k] >= cuts[i]) & (self.cate_preds_val_[:, k] <= cuts[i + 1]) + # Proportion of validations set in group + probs[i] = np.mean(ind) + # Group average treatment effect (GATE) -- average of DR outcomes in group + gate[i] = np.mean(self.dr_val_[ind, k]) + se_gate[i] = np.std(self.dr_val_[ind, k]) / np.sqrt(np.sum(ind)) + # Average of CATE predictions in group + g_cate[i] = np.mean(self.cate_preds_val_[ind, k]) + se_g_cate[i] = np.std(self.cate_preds_val_[ind, k]) / np.sqrt(np.sum(ind)) + + # Calculate group calibration score + cal_score_g = np.sum(abs(gate - g_cate) * probs) + # Calculate overall calibration score + cal_score_o = np.sum(abs(gate - self.ate_val[k]) * probs) + # Calculate R-square calibration score + cal_r_squared[k] = 1 - (cal_score_g / cal_score_o) + + df_plot1 = pd.DataFrame({'ind': np.array(range(n_groups)), + 'gate': gate, 'se_gate': se_gate, + 'g_cate': g_cate, 'se_g_cate': se_g_cate}) + df_plot1['tmt'] = self.treatments[k + 1] + df_plot = pd.concat((df_plot, df_plot1)) + + self.cal_res = CalibrationEvaluationResults( + cal_r_squared=cal_r_squared, + df_plot=df_plot, + treatments=self.treatments + ) + + return self.cal_res + + def evaluate_blp( + self, + Xval: np.array = None, + Xtrain: np.array = None + ) -> BLPEvaluationResults: + """ + Implements the best linear predictor (BLP) test as in [Chernozhukov2022]. `fit_nusiance' method must already + be implemented. + + Parameters + ---------- + Xval: (n_val x k) matrix, default ``None'' + Validation sample features for CATE model. If not specified, then `fit_cate' method must already have been + implemented + Xtrain: (n_train x k) matrix, default ``None'' + Training sample features for CATE model. If specified, then CATE is fitted on training sample and applied + to Xval. If specified, then Xtrain, Dtrain, Ytrain must have been specified in `fit_nuisance' method (and + vice-versa) + + Returns + ------- + BLPEvaluationResults object showing the results of the BLP test + """ + + if not hasattr(self, 'dr_val_'): + raise Exception("Must fit nuisances before evaluating") + + # if CATE is given explicitly or has not been fitted at all previously, fit it now + if not hasattr(self, 'cate_preds_val_'): + if Xval is None: # need at least Xval + raise Exception('CATE predictions not yet calculated - must provide Xval') + self.get_cate_preds(Xval, Xtrain) + + if self.n_treat == 1: # binary treatment + reg = OLS(self.dr_val_, add_constant(self.cate_preds_val_)).fit() + params = [reg.params[1]] + errs = [reg.bse[1]] + pvals = [reg.pvalues[1]] + else: # categorical treatment + params = [] + errs = [] + pvals = [] + for k in range(self.n_treat): # run a separate regression for each + reg = OLS(self.dr_val_[:, k], add_constant(self.cate_preds_val_[:, k])).fit(cov_type='HC1') + params.append(reg.params[1]) + errs.append(reg.bse[1]) + pvals.append(reg.pvalues[1]) + + self.blp_res = BLPEvaluationResults( + params=params, + errs=errs, + pvals=pvals, + treatments=self.treatments + ) + + return self.blp_res + + def evaluate_qini( + self, + Xval: np.array = None, + Xtrain: np.array = None, + percentiles: np.array = np.linspace(5, 95, 50) + ) -> QiniEvaluationResults: + """ + Calculates QINI coefficient for the given model as in Radcliffe (2007), where units are ordered by predicted + CATE values and a running measure of the average treatment effect in each cohort is kept as we progress + through ranks. The QINI coefficient is then the area under the resulting curve, with a value of 0 interpreted + as corresponding to a model with randomly assigned CATE coefficients. All calculations are performed on + validation dataset results, using the training set as input. + + Parameters + ---------- + Xval: (n_val x k) matrix, default ``None'' + Validation sample features for CATE model. If not specified, then `fit_cate' method must already have been + implemented + Xtrain: (n_train x k) matrix, default ``None'' + Training sample features for CATE model. If specified, then CATE is fitted on training sample and applied + to Xval. If specified, then Xtrain, Dtrain, Ytrain must have been specified in `fit_nuisance' method (and + vice-versa) + percentiles: one-dimensional array, default ``np.linspace(5, 95, 50)'' + Array of percentiles over which the QINI curve should be constructed. Defaults to 5%-95% in intervals of + 5%. + + Returns + ------- + QiniEvaluationResults object showing the results of the QINI fit + """ + if not hasattr(self, 'dr_val_'): + raise Exception("Must fit nuisances before evaluating") + + if (not hasattr(self, 'cate_preds_train_')) or (not hasattr(self, 'cate_preds_val_')): + if (Xval is None) or (Xtrain is None): + raise Exception('CATE predictions not yet calculated - must provide both Xval, Xtrain') + self.get_cate_preds(Xval, Xtrain) + + if self.n_treat == 1: + qini, qini_err = calc_qini_coeff( + self.cate_preds_train_, + self.cate_preds_val_, + self.dr_val_, + percentiles + ) + qinis = [qini] + errs = [qini_err] + else: + qinis = [] + errs = [] + for k in range(self.n_treat): + qini, qini_err = calc_qini_coeff( + self.cate_preds_train_[:, k], + self.cate_preds_val_[:, k], + self.dr_val_[:, k], + percentiles + ) + + qinis.append(qini) + errs.append(qini_err) + + pvals = [st.norm.sf(abs(q / e)) for q, e in zip(qinis, errs)] + + self.qini_res = QiniEvaluationResults( + params=qinis, + errs=errs, + pvals=pvals, + treatments=self.treatments + ) + + return self.qini_res + + def evaluate_all( + self, + Xval: np.array = None, + Xtrain: np.array = None, + n_groups: int = 4 + ) -> EvaluationResults: + """ + Implements the best linear prediction (`evaluate_blp'), calibration (`evaluate_cal') and QINI coefficient + (`evaluate_qini') methods. + + Parameters + ---------- + Xval: (n_cal x k) matrix, default ``None'' + Validation sample features for CATE model. If not specified, then `fit_cate' method must already have been + implemented + Xtrain: (n_train x k) matrix, default ``None'' + Training sample features for CATE model. If not specified, then `fit_cate' method must already have been + implemented + + Returns + ------- + EvaluationResults object summarizing the results of all tests + """ + + if (not hasattr(self, 'cate_preds_train_')) or (not hasattr(self, 'cate_preds_val_')): + if (Xval is None) or (Xtrain is None): + raise Exception('CATE predictions not yet calculated - must provide both Xval, Xtrain') + self.get_cate_preds(Xval, Xtrain) + + blp_res = self.evaluate_blp() + cal_res = self.evaluate_cal(n_groups=n_groups) + qini_res = self.evaluate_qini() + + self.res = EvaluationResults( + blp_res=blp_res, + cal_res=cal_res, + qini_res=qini_res + ) + + return self.res diff --git a/econml/validate/results.py b/econml/validate/results.py new file mode 100644 index 000000000..0b3edaeaa --- /dev/null +++ b/econml/validate/results.py @@ -0,0 +1,245 @@ +import numpy as np +import pandas as pd + +from typing import List + + +class CalibrationEvaluationResults: + """ + Results class for calibration test. + + Parameters + ---------- + cal_r_squared: list or numpy array of floats + Sequence of calibration R^2 values + + df_plot: pandas dataframe + Dataframe containing necessary data for plotting calibration test GATE results + + treatments: list or numpy array of floats + Sequence of treatment labels + """ + def __init__( + self, + cal_r_squared: np.array, + df_plot: pd.DataFrame, + treatments: np.array + ): + self.cal_r_squared = cal_r_squared + self.df_plot = df_plot + self.treatments = treatments + + def summary(self) -> pd.DataFrame: + """ + Constructs dataframe summarizing the results of the calibration test. + + Parameters + ---------- + None + + Returns + ------- + pandas dataframe containing summary of calibration test results + """ + + res = pd.DataFrame({ + 'treatment': self.treatments[1:], + 'cal_r_squared': self.cal_r_squared, + }).round(3) + return res + + def plot_cal(self, tmt: int): + """ + Plots group average treatment effects (GATEs) and predicted GATEs by quantile-based group in validation sample. + + Parameters + ---------- + tmt: integer + Treatment level to plot + + Returns + ------- + matplotlib plot with predicted GATE on x-axis and GATE (and 95% CI) on y-axis + """ + if tmt == 0: + raise Exception('Plotting only supported for treated units (not controls)') + + df = self.df_plot + df = df[df.tmt == tmt].copy() + rsq = round(self.cal_r_squared[np.where(self.treatments == tmt)[0][0] - 1], 3) + df['95_err'] = 1.96 * df['se_gate'] + fig = df.plot( + kind='scatter', + x='g_cate', + y='gate', + yerr='95_err', + xlabel='Group Mean CATE', + ylabel='GATE', + title=f"Treatment = {tmt}, Calibration R^2 = {rsq}" + ) + + return fig + + +class BLPEvaluationResults: + """ + Results class for BLP test. + + Parameters + ---------- + params: list or numpy array of floats + Sequence of estimated coefficient values + + errs: list or numpy array of floats + Sequence of estimated coefficient standard errors + + pvals: list or numpy array of floats + Sequence of estimated coefficient p-values + + treatments: list or numpy array of floats + Sequence of treatment labels + """ + def __init__( + self, + params: List[float], + errs: List[float], + pvals: List[float], + treatments: np.array + ): + self.params = params + self.errs = errs + self.pvals = pvals + self.treatments = treatments + + def summary(self): + """ + Constructs dataframe summarizing the results of the BLP test. + + Parameters + ---------- + None + + Returns + ------- + pandas dataframe containing summary of BLP test results + """ + res = pd.DataFrame({ + 'treatment': self.treatments[1:], + 'blp_est': self.params, + 'blp_se': self.errs, + 'blp_pval': self.pvals + }).round(3) + return res + + +class QiniEvaluationResults: + """ + Results class for QINI test. + + Parameters + ---------- + params: list or numpy array of floats + Sequence of estimated QINI coefficient values + + errs: list or numpy array of floats + Sequence of estimated QINI coefficient standard errors + + pvals: list or numpy array of floats + Sequence of estimated QINI coefficient p-values + + treatments: list or numpy array of floats + Sequence of treatment labels + """ + def __init__( + self, + params: List[float], + errs: List[float], + pvals: List[float], + treatments: np.array + ): + self.params = params + self.errs = errs + self.pvals = pvals + self.treatments = treatments + + def summary(self): + """ + Constructs dataframe summarizing the results of the QINI test. + + Parameters + ---------- + None + + Returns + ------- + pandas dataframe containing summary of QINI test results + """ + res = pd.DataFrame({ + 'treatment': self.treatments[1:], + 'qini_est': self.params, + 'qini_se': self.errs, + 'qini_pval': self.pvals + }).round(3) + return res + + +class EvaluationResults: + """ + Results class for combination of all tests. + + Parameters + ---------- + cal_res: CalibrationEvaluationResults object + Results object for calibration test + + blp_res: BLPEvaluationResults object + Results object for BLP test + + qini_res: QiniEvaluationResults object + Results object for QINI test + """ + def __init__( + self, + cal_res: CalibrationEvaluationResults, + blp_res: BLPEvaluationResults, + qini_res: QiniEvaluationResults + ): + self.cal = cal_res + self.blp = blp_res + self.qini = qini_res + + def summary(self): + """ + Constructs dataframe summarizing the results of all 3 tests. + + Parameters + ---------- + None + + Returns + ------- + pandas dataframe containing summary of all test results + """ + res = self.blp.summary().merge( + self.qini.summary(), + on='treatment' + ).merge( + self.cal.summary(), + on='treatment' + ) + return res + + def plot_cal(self, tmt: int): + """ + Plots group average treatment effects (GATEs) and predicted GATEs by quantile-based group in validation sample. + + Parameters + ---------- + tmt: integer + Treatment level to plot + + Returns + ------- + matplotlib plot with predicted GATE on x-axis and GATE (and 95% CI) on y-axis + """ + return self.cal.plot_cal(tmt) diff --git a/econml/validate/utils.py b/econml/validate/utils.py new file mode 100644 index 000000000..1ed45225f --- /dev/null +++ b/econml/validate/utils.py @@ -0,0 +1,94 @@ +from typing import Tuple + +import numpy as np + + +def calculate_dr_outcomes( + D: np.array, + y: np.array, + reg_preds: np.array, + prop_preds: np.array +) -> np.array: + """ + Calculates doubly-robust (DR) outcomes using predictions from nuisance models + + Parameters + ---------- + D: vector of length n + Treatment assignments. Should have integer values with the lowest-value corresponding to the + control treatment. It is recommended to have the control take value 0 and all other treatments be integers + starting at 1 + y: vector of length n + Outcomes + reg_preds: (n x n_treat) matrix + Outcome predictions for each potential treatment + prop_preds: (n x n_treat) matrix + Propensity score predictions for each treatment + + Returns + ------- + Doubly robust outcome values + """ + + # treat each treatment as a separate regression + # here, prop_preds should be a matrix + # with rows corresponding to units and columns corresponding to treatment statuses + dr_vec = [] + d0_mask = np.where(D == 0, 1, 0) + y_dr_0 = reg_preds[:, 0] + (d0_mask / np.clip(prop_preds[:, 0], .01, np.inf)) * (y - reg_preds[:, 0]) + for k in np.sort(np.unique(D)): # pick a treatment status + if k > 0: # make sure it is not control + dk_mask = np.where(D == k, 1, 0) + y_dr_k = reg_preds[:, k] + (dk_mask / np.clip(prop_preds[:, k], .01, np.inf)) * (y - reg_preds[:, k]) + dr_k = y_dr_k - y_dr_0 # this is an n x 1 vector + dr_vec.append(dr_k) + dr = np.column_stack(dr_vec) # this is an n x n_treatment matrix + + return dr + + +def calc_qini_coeff( + cate_preds_train: np.array, + cate_preds_val: np.array, + dr_val: np.array, + percentiles: np.array +) -> Tuple[float, float]: + """ + Helper function for QINI coefficient calculation. See documentation for "evaluate_qini" method + for more details. + + Parameters + ---------- + cate_preds_train: (n_train x n_treatment) matrix + Predicted CATE values for the training sample. + cate_preds_val: (n_val x n_treatment) matrix + Predicted CATE values for the validation sample. + dr_val: (n_val x n_treatment) matrix + Doubly robust outcome values for each treatment status in validation sample. Each value is relative to + control, e.g. for treatment k the value is Y(k) - Y(0), where 0 signifies no treatment. + percentiles: one-dimensional array + Array of percentiles over which the QINI curve should be constructed. Defaults to 5%-95% in intervals of 5%. + + Returns + ------- + QINI coefficient and associated standard error. + """ + qs = np.percentile(cate_preds_train, percentiles) + toc, toc_std, group_prob = np.zeros(len(qs)), np.zeros(len(qs)), np.zeros(len(qs)) + toc_psi = np.zeros((len(qs), dr_val.shape[0])) + n = len(dr_val) + ate = np.mean(dr_val) + for it in range(len(qs)): + inds = (qs[it] <= cate_preds_val) # group with larger CATE prediction than the q-th quantile + group_prob = np.sum(inds) / n # fraction of population in this group + toc[it] = group_prob * ( + np.mean(dr_val[inds]) - ate) # tau(q) = q * E[Y(1) - Y(0) | tau(X) >= q[it]] - E[Y(1) - Y(0)] + toc_psi[it, :] = np.squeeze( + (dr_val - ate) * (inds - group_prob) - toc[it]) # influence function for the tau(q) + toc_std[it] = np.sqrt(np.mean(toc_psi[it] ** 2) / n) # standard error of tau(q) + + qini_psi = np.sum(toc_psi[:-1] * np.diff(percentiles).reshape(-1, 1) / 100, 0) + qini = np.sum(toc[:-1] * np.diff(percentiles) / 100) + qini_stderr = np.sqrt(np.mean(qini_psi ** 2) / n) + + return qini, qini_stderr diff --git a/notebooks/CATE validation.ipynb b/notebooks/CATE validation.ipynb new file mode 100644 index 000000000..57b20444b --- /dev/null +++ b/notebooks/CATE validation.ipynb @@ -0,0 +1,543 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + }, + "pycharm": { + "name": "#%%\n" + }, + "scrolled": true + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/shap/utils/_clustering.py:35: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + " def _pt_shuffle_rec(i, indexes, index_mask, partition_tree, M, pos):\n", + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/shap/utils/_clustering.py:54: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + " def delta_minimization_order(all_masks, max_swap_size=100, num_passes=2):\n", + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/shap/utils/_clustering.py:63: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + " def _reverse_window(order, start, length):\n", + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/shap/utils/_clustering.py:69: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + " def _reverse_window_score_gain(masks, order, start, length):\n", + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/shap/utils/_clustering.py:77: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + " def _mask_delta_score(m1, m2):\n", + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/shap/links.py:5: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + " def identity(x):\n", + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/shap/links.py:10: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + " def _identity_inverse(x):\n", + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/shap/links.py:15: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + " def logit(x):\n", + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/shap/links.py:20: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + " def _logit_inverse(x):\n", + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/shap/utils/_masked_model.py:362: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + " def _build_fixed_single_output(averaged_outs, last_outs, outputs, batch_positions, varying_rows, num_varying_rows, link, linearizing_weights):\n", + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/shap/utils/_masked_model.py:384: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + " def _build_fixed_multi_output(averaged_outs, last_outs, outputs, batch_positions, varying_rows, num_varying_rows, link, linearizing_weights):\n", + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/shap/maskers/_tabular.py:185: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + " def _single_delta_mask(dind, masked_inputs, last_mask, data, x, noop_code):\n", + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/shap/maskers/_tabular.py:196: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + " def _delta_masking(masks, x, curr_delta_inds, varying_rows_out,\n", + "/opt/anaconda3/envs/cate_test/lib/python3.9/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n", + " from .autonotebook import tqdm as notebook_tqdm\n", + "The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n", + "The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import scipy.stats as st\n", + "from sklearn.ensemble import RandomForestClassifier, GradientBoostingRegressor\n", + "\n", + "from econml.metalearners import TLearner\n", + "from econml.dml import DML\n", + "\n", + "from econml.validate.drtester import DRtester" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + }, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "## Generate Data\n", + "\n", + "Setting with 5 continuous covariates, 2 categorical treatments (and a control), and a continous outcome. Treatment assignment is based on the covariates and random noise. The outcome is based on covariates, treatment, and random noise. For the second treatment, the effect is conditional on the second covariate (X[:, 1])." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "np.random.seed(123)\n", + "\n", + "N = 20000 # number of units\n", + "K = 5 # number of covariates\n", + "num_treatments = 2 # number of treatments (excluding control)\n", + "\n", + "# Generate random Xs\n", + "X_mu = np.zeros(5) # Means of Xs\n", + "# Random covariance matrix of Xs\n", + "X_sig = np.diag(np.random.rand(5))\n", + "X = st.multivariate_normal(X_mu, X_sig).rvs(N)\n", + "\n", + "# Effect of Xs on outcome\n", + "X_beta = np.random.uniform(0, 5, K)\n", + "# Effect of treatment on outcomes\n", + "D_beta = np.array([0, 1, 2])\n", + "# Effect of treatment on outcome conditional on X1\n", + "DX1_beta = np.array([0, 0, 3])\n", + "\n", + "# Generate treatments based on X and random noise\n", + "beta_treat = np.random.uniform(-1, 1, (num_treatments + 1, K))\n", + "D1 = np.zeros((N, num_treatments + 1))\n", + "for k in range(num_treatments + 1):\n", + " D1[:, k] = X @ beta_treat[k, :] + np.random.gumbel(0, 1, N)\n", + "D = np.array([np.where(D1[i, :] == np.max(D1[i, :]))[0][0] for i in range(N)])\n", + "D_dum = pd.get_dummies(D)\n", + "\n", + "# Generate Y (based on X, D, and random noise)\n", + "Y_sig = 1 # Variance of random outcome noise\n", + "Y = X @ X_beta + (D_dum @ D_beta) + X[:, 1] * (D_dum @ DX1_beta) + np.random.normal(0, Y_sig, N)\n", + "Y = Y.to_numpy()\n", + "\n", + "# Split into training/validation samples\n", + "train_prop = .5\n", + "train_N = np.ceil(train_prop * N)\n", + "ind = np.array(range(N))\n", + "train_ind = np.random.choice(N, int(train_N), replace=False)\n", + "val_ind = ind[~np.isin(ind, train_ind)]\n", + "\n", + "Xtrain, Dtrain, Ytrain = X[train_ind], D[train_ind], Y[train_ind]\n", + "Xval, Dval, Yval = X[val_ind], D[val_ind], Y[val_ind]" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + }, + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "## Validation\n", + "We use a simple (un-tuned) random forest classifier to predict treatment and an un-tuned gradient boosted regressor to predict outcomes. In the true data, there is heterogeneity by the second covariate (X[:, 1]) for the second treatment (D = 2) but not for the first (D = 1). Therefore, we should see a more significant blp_est, more significant qini coefficient, and higher calibration r-squared for the second treatment. We can also plot the calibration for each treatment. There should be a strong positive relationship for treatment 2 but a flat relationship for treatment 1.\n", + "\n", + "We test two CATE models: a double machine learning (DML) estimator and a simple T-learner." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "model_regression = GradientBoostingRegressor(random_state=0)\n", + "model_propensity = RandomForestClassifier(random_state=0)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "`sparse` was renamed to `sparse_output` in version 1.2 and will be removed in 1.4. `sparse_output` is ignored unless you leave `sparse` to its default value.\n", + "The final model has a nonzero intercept for at least one outcome; it will be subtracted, but consider fitting a model without an intercept if possible.\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "est_t = TLearner(models=model_regression)\n", + "est_dm = DML(model_y=model_regression, model_t=model_propensity, model_final=model_regression)\n", + "\n", + "est_t.fit(Ytrain, Dtrain, X=Xtrain)\n", + "est_dm.fit(Ytrain, Dtrain, X=Xtrain)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### DML" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
treatmentblp_estblp_seblp_pvalqini_estqini_seqini_pvalcal_r_squared
01-0.1370.1420.335-0.0150.0210.242-5.506
121.2090.0950.0000.3730.0240.0000.090
\n", + "
" + ], + "text/plain": [ + " treatment blp_est blp_se blp_pval qini_est qini_se qini_pval \\\n", + "0 1 -0.137 0.142 0.335 -0.015 0.021 0.242 \n", + "1 2 1.209 0.095 0.000 0.373 0.024 0.000 \n", + "\n", + " cal_r_squared \n", + "0 -5.506 \n", + "1 0.090 " + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Initialize DRtester and fit/predict nuisance models\n", + "dml_tester = DRtester(\n", + " model_regression=model_regression, \n", + " model_propensity=model_propensity, \n", + " cate=est_dm\n", + ").fit_nuisance(Xval, Dval, Yval, Xtrain, Dtrain, Ytrain)\n", + "\n", + "res_dml = dml_tester.evaluate_all(Xval, Xtrain)\n", + "res_dml.summary()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHHCAYAAABDUnkqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8pElEQVR4nO3dd3xUVf7/8fcQkkkoAQIkBAyEotKLQdgQaRo2AhtWV4RFpamAGpQvrFJXAijFAqJIEaTpIs1FVxcWFlkigrAoJixKE2kqhE4qEJI5vz/8ZXRIIROSTLi8no/HPPSeOefOZ06GyTv3njtjM8YYAQAAWEQZTxcAAABQlAg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AEqVCRMmyGazubSFhoZqwIABzu0lS5bIZrPp66+/LuHqXMXFxclmsykuLs6jdQBwRbiBW2w2W4FuJfVmP2fOHC1ZsqREHutG7d27VxMmTNDRo0dL/LFTU1MVGxur+++/XwEBAbLZbEU2b8nJyZo4caJatGihChUqyM/PT02bNtWoUaN04sSJInkMTyuNr7NOnTq5/Jvz8/NT8+bNNXPmTDkcjnzHHj9+XDVr1pTNZtPrr7+eZ7/9+/dr5MiRatmypSpWrKjg4GB1797d46Hyt7IDZm63HTt2XHd8dpi+9ubr65tr/4ULF6pRo0by9fXV7bffrlmzZuW575UrVyo8PFzly5dX5cqV1a5dO/3nP/8p9HNFwZX1dAG4ubz//vsu2++99542btyYo71Ro0YlUs+cOXNUrVo1l7/qS6u9e/dq4sSJ6tSpk0JDQ0v0sc+ePatJkyapdu3aatGiRZGFz8OHDysyMlLHjx/Xww8/rMGDB8vHx0f/+9//tHDhQn300Uc6ePDgDT/OgQMHVKaM5/4Wy+t11qFDB126dEk+Pj4eqeu2227T1KlTJf3yM/7ggw80fPhwnTlzRpMnT851zIULF9S1a1elpaWpQ4cOGjlypEJCQtS7d+8cfd99910tXLhQDz30kJ555hklJSXpnXfe0e9+9zutX79ekZGRxfr83PHcc8/p7rvvdmlr0KBBgcfPnTtXFSpUcG57eXnl6PPOO+/oqaee0kMPPaQRI0boiy++0HPPPaf09HSNGjXKpe+ECRM0adIk9ezZUwMGDNDVq1f17bff6ueff3bzmaFQDHADYmJiTEFeRmlpacXy+E2aNDEdO3Ysln0XtdWrVxtJZvPmzSX+2JcvXzYnT540xhjz1VdfGUlm8eLFN7TPq1evmhYtWphy5cqZL774Isf9SUlJZuzYsW7vNzY29rqvqcWLFxtJ5quvvnJ7/w6Hw6Snp7s1pjS+zjp27GiaNGni0nbp0iVTp04dU7FiRZOZmZljzOXLl02HDh2Mv7+/2b59u7l8+bLp3r27sdvt5vPPP8/R/+uvvzYpKSkubWfPnjXVq1c3ERERRfuECmnz5s1Gklm9enWhxme/3s6cOZNvv/T0dFO1alXTvXt3l/ZHH33UlC9f3pw/f97Ztn37dmOz2cyMGTMKVRNuHKelUOQ6deqkpk2bateuXerQoYPKlSunsWPHSpKuXLmi2NhYNWjQQHa7XSEhIRo5cqSuXLniso/Fixfr3nvvVWBgoOx2uxo3bqy5c+e69AkNDdV3332nzz//3HkouVOnTpJ+XZOxdetWPffcc6pevboqV66sIUOGKCMjQxcvXlS/fv1UpUoVValSRSNHjpQxxmX/DodDM2fOVJMmTeTr66ugoCANGTJEFy5cyFHHH/7wB23dulVt2rSRr6+v6tWrp/fee8/ZZ8mSJXr44YclSZ07dy7x03d2u101atQo0n3+/e9/1+7duzVu3Djdc889Oe739/d3OXrwxRdf6OGHH1bt2rWdP/vhw4fr0qVL132sa9fcZEtPT9eQIUNUtWpV+fv7q1+/fnn+fDZs2KDWrVvLz89P77zzjqQbf53lteZm9erVCgsLk5+fn6pVq6bHHnssx1/sAwYMUIUKFfTzzz/rgQceUIUKFVS9enU9//zzysrKuu6c5MbX11d33323UlJSdPr0aZf7jDHq37+/du/erY0bN+p3v/ud7Ha71qxZoy5duuiBBx7Q3r17XcaEhYW5HM2QpKpVq6p9+/bat29foWosTikpKcrMzCzUWGOMkpOTc7wPZNu8ebPOnTunZ555xqU9JiZGaWlpWrt2rbNt5syZqlGjhoYNGyZjjFJTUwtVEwqP01IoFufOnVPXrl315z//WY899piCgoLkcDjUo0cPbd26VYMHD1ajRo20Z88evfHGGzp48KA+/vhj5/i5c+eqSZMm6tGjh8qWLatPP/1UzzzzjBwOh2JiYiT98gby7LPPqkKFCho3bpwkKSgoyKWOZ599VjVq1NDEiRO1Y8cOzZ8/X5UrV9aXX36p2rVra8qUKVq3bp1ee+01NW3aVP369XOOHTJkiJYsWaKBAwfqueee05EjR/T2228rPj5e27Ztk7e3t7PvoUOH1LNnTz3xxBPq37+/Fi1apAEDBigsLExNmjRRhw4d9Nxzz+mtt97S2LFjnaft8jt9d+XKFaWkpBRovqtVq1agfkXpk08+kST17du3QP1Xr16t9PR0Pf3006patap27typWbNm6aefftLq1asLVcPQoUNVuXJlTZgwQQcOHNDcuXN17NgxZ+jIduDAAfXp00dDhgzRoEGDdOedd0oqutfZb2W/Zu6++25NnTpVp06d0ptvvqlt27YpPj5elStXdvbNyspSVFSU2rZtq9dff12fffaZpk+frvr16+vpp58u1JwcPXpUNpvN5XEkaeTIkdqwYYM2btzocvrGx8dHf//739WzZ0917dpVO3bsUHBwcL6PkZiYWKDX3NWrV5WUlFSgugMCAm7o1OPAgQOVmpoqLy8vtW/fXq+99ppat25d4PH16tVTamqqypcvrwceeEDTp093+TnHx8dLUo59hoWFqUyZMoqPj9djjz0mSdq0aZPatWunt956Sy+//LLOnTunGjVqaNy4cRo6dGihnyPc4NHjRrjp5XZaqmPHjkaSmTdvnkv7+++/b8qUKZPjFMa8efOMJLNt2zZnW26nDaKioky9evVc2vI6XZB92iIqKso4HA5ne3h4uLHZbOapp55ytmVmZprbbrvNZT9ffPGFkWSWLVvmst/169fnaK9Tp46RZLZs2eJsO336tLHb7eYvf/mLs83d01LZz6EgN3cU1WmpVq1amUqVKhW4f24/06lTpxqbzWaOHTvmbMvttFSdOnVM//79ndvZcxMWFmYyMjKc7a+++qqRZP7xj3+4jJVk1q9fX6Ca3HmdZZ8Syf6ZZmRkmMDAQNO0aVNz6dIlZ79//vOfRpIZP368s61///5Gkpk0aZLLPlu1amXCwsJyPNa1OnbsaBo2bGjOnDljzpw5Y/bv329eeOEFIynHqZOitGXLFmOz2cyLL7543b7Z81OQ25EjRwpVz7Zt28xDDz1kFi5caP7xj3+YqVOnmqpVqxpfX1/zzTffXHf8zJkzzdChQ82yZcvMhx9+aIYNG2bKli1rbr/9dpOUlOTsFxMTY7y8vHLdR/Xq1c2f//xnY4wx58+fN5JM1apVTYUKFcxrr71mVq5cae6///5c3xdRPDhyg2Jht9s1cOBAl7bVq1erUaNGatiwoc6ePetsv/feeyX9cti3Xbt2kiQ/Pz/n/UlJSbp69ao6duyoDRs2KCkpSZUqVSpQHU888YTLX/Bt27bV9u3b9cQTTzjbvLy81Lp1a+3atcul1kqVKqlLly4utWYfpt+8ebMeeeQRZ3vjxo3Vvn1753b16tV155136vDhwwWqMzdRUVHauHFjoccXt+TkZFWsWLHA/X/7M01LS9OlS5fUrl07GWMUHx+v2rVru13D4MGDXY6gPf300xo7dqzWrVunHj16ONvr1q2rqKiofGu6kddZtq+//lqnT5/WhAkTXK626d69uxo2bKi1a9dq4sSJLmOeeuopl+327dvnWKCfl/3796t69eoubT169NDChQvdqrugTp8+rUceeUR169bVyJEjr9u/RYsWBX4NF/a0abt27ZzvG9Ivz79nz55q3ry5xowZo/Xr1+c7ftiwYS7bDz30kNq0aaNHH31Uc+bM0ejRoyUp34Xjvr6+ztOr2aegzp07pxUrVjgXavfs2VPNmjXTyy+/rCFDhhTquaLgCDcoFrVq1crxRvD9999r3759Od6Ms/12jcC2bdsUGxur7du3Kz093aWfO790rv2FmT0uJCQkR/tv12p8//33SkpKUmBg4HVrze1xJKlKlSo51n+4Izg4+LqnBzzJ39/frfB2/PhxjR8/Xp988kmOeSnoqYtr3X777S7bFSpUUHBwcI7L7evWrZvr+KJ6nWU7duyYJDlPe/1Ww4YNtXXrVpc2X1/fHP8e3HndhIaGasGCBXI4HPrhhx80efJknTlzJs/LmG9EWlqa/vCHPyglJUVbt27NsRYnN1WqVCmSK6qysrJ05swZl7aAgIA8w0aDBg30xz/+UWvWrFFWVlauVz7l55FHHtFf/vIXffbZZ85w4+fnp4yMjFz7X7582RmUs//r7e2tnj17OvuUKVNGvXv3VmxsrI4fP16oMI+CI9ygWPz2L+JsDodDzZo104wZM3Idkx04fvjhB913331q2LChZsyYoZCQEPn4+GjdunV64403rvsZHr+V15tabu3mNwsJHQ6HAgMDtWzZslzHX/sLKa/HMXksTiyIS5cuFfiXflEvFi6Ihg0bKj4+Xj/++GOOsHitrKwsdenSRefPn9eoUaPUsGFDlS9fXj///LMGDBjg1s+0MHJ7PRbl66yw3P2le63y5cu7hIeIiAjdddddGjt2rN56660bLc8pIyNDf/rTn/S///1PGzZsUNOmTQs87vz58wXqW7169Tzn48cff8wRUDdv3uxc2J2bkJAQZWRkKC0tTf7+/gWq4drxv609ODhYWVlZOn36tMsfPRkZGTp37pxq1qwp6ZfQ5evrq8qVK+d4PtnjLly4QLgpZoQblJj69etr9+7duu+++3J8Au1vffrpp7py5Yo++eQTlzeAzZs35+ib335utNbPPvtMERERuf5iLAx3a125cmWOU3t5uZEQVVjR0dFavny5/va3v2nMmDH59t2zZ48OHjyopUuXuizavtHTbt9//706d+7s3E5NTdXJkyfVrVu3644tjtdZnTp1JP2ygDn7dGu2AwcOOO8vLs2bN9djjz2md955R88//3yR/AJ1OBzq16+fNm3apFWrVqljx44FHvvll1+6/Hzyc+TIkTw//6lGjRo5XistWrTId3+HDx+Wr69vgY4wXcsYo6NHj6pVq1bOtpYtW0r65dTjb19fX3/9tRwOh/P+MmXKqGXLlvrqq6+UkZHhcnQp+0Mt8zp6jaJDuEGJ6dWrl9atW6cFCxZo8ODBLvddunRJDodD5cuXd/6189tf2ElJSVq8eHGOfZYvX14XL14sllrnzJmjl156SVOmTHG5LzMzU6mpqTmuRrme8uXLS1KB6y3ta2569uypqVOnavLkyerUqZPCw8Nd7k9JSdG0adM0efLkXH+mxhi9+eabN1TD/PnzNXDgQOe6m7lz5yozM1Ndu3a97tjieJ21bt1agYGBmjdvnh5//HHZ7XZJ0r/+9S/t27dP48ePL8jTuiEjR47Ue++9pxkzZmjmzJk3vL9nn31WK1eu1DvvvKM//elPbo0tqjU3vr6+eZ7eOnPmTI6wsHv3bn3yySfq2rWryxVYx48fV3p6uho2bJjv+Llz5+rMmTO6//77nW333nuvAgICNHfuXJdwM3fuXJUrV07du3d3tvXu3Vs7duzQ0qVLNWjQIEm/nLpatmyZGjdu7DzKg+JDuEGJ6du3r1atWqWnnnpKmzdvVkREhLKysrR//36tWrXK+Tkkv//97+Xj46Po6GgNGTJEqampWrBggQIDA3Xy5EmXfYaFhWnu3Ll6+eWX1aBBAwUGBub4i7kwOnbsqCFDhmjq1KlKSEjQ73//e3l7e+v777/X6tWr9eabb7qcTy+Ili1bysvLS6+88oqSkpJkt9udn7GSm6Jec/P222/r4sWLzr8eP/30U/3000+SfvkFlr2+JPtS5sWLF+f7yc/e3t5as2aNIiMj1aFDB/Xq1UsRERHy9vbWd999pw8++EBVqlTR5MmT1bBhQ9WvX1/PP/+8fv75Z/n7++vvf//7Da1Jkn45JXDfffepV69eOnDggObMmaN77rnHZTFxXorjdebt7a1XXnlFAwcOVMeOHdWnTx/npeChoaEaPnz4DT3fgmjcuLG6deumd999Vy+++KKqVq1a6H3NnDlTc+bMUXh4uMqVK6e//e1vLvc/+OCDztCem6Jac5Of3r17y8/PT+3atVNgYKD27t2r+fPnq1y5cpo2bZpL3379+unzzz93CbR16tRR79691axZM/n6+mrr1q1asWKFWrZs6bLw18/PTy+99JJiYmL08MMPKyoqSl988YX+9re/afLkyQoICHD2HTJkiN59913FxMTo4MGDql27tt5//30dO3ZMn376abHOB/4/T12mBWvI61Lwaz85NVtGRoZ55ZVXTJMmTYzdbjdVqlQxYWFhZuLEiS6XXX7yySemefPmxtfX14SGhppXXnnFLFq0KMclo4mJiaZ79+6mYsWKRpLzct28PsE2r08j7d+/vylfvnyOeufPn2/CwsKMn5+fqVixomnWrJkZOXKkOXHihLNPnTp1cr30tmPHjjkuH16wYIGpV6+e8fLyKvFPK86+JDq322/ndNasWXleOp2bCxcumPHjx5tmzZqZcuXKGV9fX9O0aVMzZswY56ciG2PM3r17TWRkpKlQoYKpVq2aGTRokNm9e3eOy9LduRT8888/N4MHDzZVqlQxFSpUMI8++qg5d+5cjrF5XRp9o6+zay8Fz7Zy5UrTqlUrY7fbTUBAgHn00UfNTz/95NInr9dcQT6h2Zj8/53FxcUZSSY2Nva6+8lP9uXqBXndeMqbb75p2rRpYwICAkzZsmVNcHCweeyxx8z333+fo2/2x1T81pNPPmkaN25sKlasaLy9vU2DBg3MqFGjTHJycq6PN3/+fHPnnXcaHx8fU79+ffPGG2+4fNxEtlOnTpn+/fubgIAAY7fbTdu2bQv8bwo3zmaMB07WAyi1evXqpaNHj2rnzp2eLgUACoXTUgCcjDGKi4vLcfoBAG4mHLkBAACWwhdnAgAASyHcAAAASyHcAAAASyHcAAAAS7nlrpZyOBw6ceKEKlasWGwf3Q8AAIqWMUYpKSmqWbOmyydP5+aWCzcnTpy47pf8AQCA0unHH3/Ubbfdlm+fWy7cVKxYUdIvk1OYb4oFAAAlLzk5WSEhIc7f4/m55cJN9qkof39/wg0AADeZgiwpYUExAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMJNEUnPyFTo6LUKHb1W6RmZni4HAIASV1p+FxJuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBugBKWnpGp0NFrFTp6rdIzMj1dDgBYDuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYikfDzZYtWxQdHa2aNWvKZrPp448/zrf/mjVr1KVLF1WvXl3+/v4KDw/Xhg0bSqZYAABwU/BouElLS1OLFi00e/bsAvXfsmWLunTponXr1mnXrl3q3LmzoqOjFR8fX8yVAgCAm0VZTz54165d1bVr1wL3nzlzpsv2lClT9I9//EOffvqpWrVqVcTVAQCAm5FHw82NcjgcSklJUUBAQJ59rly5oitXrji3k5OTS6I0AADgITf1guLXX39dqamp6tWrV559pk6dqkqVKjlvISEhJVghAAAoaTdtuPnggw80ceJErVq1SoGBgXn2GzNmjJKSkpy3H3/8sQSrBAAAJe2mPC21YsUKPfnkk1q9erUiIyPz7Wu322W320uoMgAA4Gk33ZGb5cuXa+DAgVq+fLm6d+/u6XIAAEAp49EjN6mpqTp06JBz+8iRI0pISFBAQIBq166tMWPG6Oeff9Z7770n6ZdTUf3799ebb76ptm3bKjExUZLk5+enSpUqeeQ5AACA0sWjR26+/vprtWrVynkZ94gRI9SqVSuNHz9eknTy5EkdP37c2X/+/PnKzMxUTEyMgoODnbdhw4Z5pH4AAFD6ePTITadOnWSMyfP+JUuWuGzHxcUVb0EAAOCmd9OtuQEAAMgP4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAbkLpGZkKHb1WoaPXKj0j09PlAKUK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4aYYHD2b7ukSAAC4ZRFuisDF9AwNWrrLud3trS/Ub+FOJaVf9WBVAADcmgg3ReC55Qna/sNZl7Zth87q2eXxHqoIAIBbF+HmBh0+k6ot35+R45r2LGO05fszOnI2zSN1AQBwqyLc3KBj5/NfX3P0HOEGAICSRLi5QXUCyuV7f2jV8iVUCQAAkAg3N6xe9QrqcHv1HBPpZbOpw+3VVbca4QYAgJJEuCkCs/q0Unj9ai5tEQ2qaVafVh6qCACAWxfhpghUKuetBf3DnNvrnmuv955oo0rlvD1YFQAAtybCTTEIrZb/OhwAAFB8CDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDeABx09m+7pEgDAcgg3QAm6mJ6hQUt3Obe7vfWF+i3cqaT0qx6sCgCshXADlKDnlido+w9nXdq2HTqrZ5fHe6giALAej4abLVu2KDo6WjVr1pTNZtPHH3983TFxcXG66667ZLfb1aBBAy1ZsqTY6wSKwuEzqdry/Rk5rmnPMkZbvj+jI2fTPFIXAFiNR8NNWlqaWrRoodmzZxeo/5EjR9S9e3d17txZCQkJ+r//+z89+eST2rBhQzFXCty4Y+fzX19z9BzhBgCKQllPPnjXrl3VtWvXAvefN2+e6tatq+nTp0uSGjVqpK1bt+qNN95QVFRUcZUJFIk6AeXyvT+0avkSqgQAit/Rs+lqXNPfI499U6252b59uyIjI13aoqKitH37dg9VBBRcveoV1OH26jn+0XnZbOpwe3XVrUa4AXDzKk0XTNxU4SYxMVFBQUEubUFBQUpOTtalS5dyHXPlyhUlJye73ABPmdWnlcLrV3Npi2hQTbP6tPJQRQBQNErTBRM3VbgpjKlTp6pSpUrOW0hIiKdLwi2sUjlvLegf5txe91x7vfdEG1Uq5+3BqgDgxpS2CyZuqnBTo0YNnTp1yqXt1KlT8vf3l5+fX65jxowZo6SkJOftxx9/LIlSgQIJrZb/OhwAuBmUtgsmPLqg2F3h4eFat26dS9vGjRsVHh6e5xi73S673V7cpQEAcMsqbRdMePTITWpqqhISEpSQkCDpl0u9ExISdPz4cUm/HHXp16+fs/9TTz2lw4cPa+TIkdq/f7/mzJmjVatWafjw4Z4oHwAAqPRdMOHRcPP111+rVatWatXql8WUI0aMUKtWrTR+/HhJ0smTJ51BR5Lq1q2rtWvXauPGjWrRooWmT5+ud999l8vAAQDwsNJ0wYRHT0t16tRJxpg878/t04c7deqk+Hg+qh4AgNIk+4KJxuN/+WDddc+153NuAACAdXjyggnCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBS3ws2qVauUkZHh3P7pp5/kcDic2+np6Xr11VeLrjoAAAA3uRVu+vTpo4sXLzq3GzdurKNHjzq3U1JSNGbMmKKqDQAAwG1uhRtjTL7bAAAAnsaaGwAAYCmEGwAAYCll3R2wYcMGVapUSZLkcDi0adMmffvtt5Lksh4HAADAE9wON/3793fZHjJkiMu2zWa7sYoAAABugFvh5reXfQMAAJRGbq25efzxx5WSklJctQAAANwwt8LN0qVLdenSpeKqBQAA4Ibd0OfcAAAAlDZuLyhOSUmRr69vvn38/f0LXRAAAMCNcDvc3HHHHXneZ4yRzWZTVlbWDRUFAABQWG6Hmw8//FABAQHFUQsAAMANczvcREREKDAwsDhqAQAAuGFF/vULnJICAACe5Fa4qVOnjry8vHK97+DBgxo5cqRuu+22IikMAACgMNwKN0eOHFHVqlWd2+np6Vq8eLHat2+vxo0ba8uWLRoxYkSRFwkAAFBQbq+5kaQdO3bo3Xff1erVq1W7dm3t27dPmzdvVvv27Yu6PgAAALe4deRm+vTpatKkiXr27KkqVapoy5Yt2rNnj2w2m8sRHQAAAE9x68jNqFGjNGrUKE2aNCnPtTcAAACe5NaRm5deekmrV69W3bp1NWrUKH377bfFVRcAAEChuBVuxowZo4MHD+r9999XYmKi2rZtqxYtWsgYowsXLhRXjQAAAAVWqM+56dixo5YuXaqTJ0/qmWee0V133aUOHTqoXbt2mjFjRlHXCAAAUGA39CF+/v7+GjJkiHbu3Kndu3erbdu2mjZtWlHVBgAA4Da3FhRfunRJmzZt0h/+8AdJv5ymunLlyq87K1tWP/zwQ9FWCAAA4Aa3ws3SpUu1du1aZ7h5++231aRJE/n5+UmSDhw4oJo1a2r48OFFXykAAEABuHVaatmyZRo8eLBL2wcffKDNmzdr8+bNevXVV7Vq1aoiLRAAAMAdboWbQ4cOqVmzZs5tX19flSnz6y7atGmjvXv3Fl11AAAAbnLrtNTFixdd1ticOXPG5X6Hw+FyPwAAQElz68jNbbfdlu8H9/3vf//jW8EBAIBHuRVuunXrpvHjx+vy5cs57rt06ZImTpyo7t27F1lxAAAA7nLrtNTYsWO1atUq3XnnnRo6dKjuuOMOSb9cJfX2228rMzNTY8eOLZZCAQAACsKtcBMUFKQvv/xSTz/9tEaPHi1jjCTJZrOpS5cumjNnjoKCgoqlUAAAgIJwK9xIUt26dbV+/XqdP39ehw4dkiQ1aNBAAQEBRV4cAACAu9wON9kCAgLUpk2boqwFAADght3Qd0sBAACUNoQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKWU9XYBVlPMpq6PTunu6DAAAbnkcuQEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAGAm9zRs+meLgEoVQg3AHCTuZieoUFLdzm3u731hfot3Kmk9KserAooPQg3AHCTeW55grb/cNalbduhs3p2ebyHKgJKF4+Hm9mzZys0NFS+vr5q27atdu7cmW//mTNn6s4775Sfn59CQkI0fPhwXb58uYSqBQDPOnwmVVu+PyPHNe1ZxmjL92d05GyaR+oCShOPhpuVK1dqxIgRio2N1TfffKMWLVooKipKp0+fzrX/Bx98oNGjRys2Nlb79u3TwoULtXLlSo0dO7aEKwcAzzh2Pv/1NUfPEW4Aj4abGTNmaNCgQRo4cKAaN26sefPmqVy5clq0aFGu/b/88ktFRETokUceUWhoqH7/+9+rT58+1z3aAwBWUSegXL73h1YtX0KVAKWXx8JNRkaGdu3apcjIyF+LKVNGkZGR2r59e65j2rVrp127djnDzOHDh7Vu3Tp169atRGoGAE+rV72COtxePcebt5fNpg63V1fdaoQbwGPfLXX27FllZWUpKCjIpT0oKEj79+/Pdcwjjzyis2fP6p577pExRpmZmXrqqafyPS115coVXblyxbmdnJxcNE8AADxkVp9WembZN9r2m0XFEQ2qaVafVh6sCig9PL6g2B1xcXGaMmWK5syZo2+++UZr1qzR2rVr9dJLL+U5ZurUqapUqZLzFhISUoIVA0DRq1TOWwv6hzm31z3XXu890UaVynl7sCqg9PDYkZtq1arJy8tLp06dcmk/deqUatSokeuYF198UX379tWTTz4pSWrWrJnS0tI0ePBgjRs3TmXK5MxqY8aM0YgRI5zbycnJBBwAlhJaLf91OMCtxmNHbnx8fBQWFqZNmzY52xwOhzZt2qTw8PBcx6Snp+cIMF5eXpIkY0yuY+x2u/z9/V1uAADAujx25EaSRowYof79+6t169Zq06aNZs6cqbS0NA0cOFCS1K9fP9WqVUtTp06VJEVHR2vGjBlq1aqV2rZtq0OHDunFF19UdHS0M+QAAIBbm0fDTe/evXXmzBmNHz9eiYmJatmypdavX+9cZHz8+HGXIzV//etfZbPZ9Ne//lU///yzqlevrujoaE2ePNlTTwEAAJQyHg03kjR06FANHTo01/vi4uJctsuWLavY2FjFxsaWQGUAAOBmdFNdLQUAAHA9hBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGAphBsAAGApHv9uKeBWU86nrI5O6+7pMgDAsjhyAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALKWspwsAAADWUM6nrI5O6+7pMjhyAwAArIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALMXj4Wb27NkKDQ2Vr6+v2rZtq507d+bb/+LFi4qJiVFwcLDsdrvuuOMOrVu3roSqBQAApV1ZTz74ypUrNWLECM2bN09t27bVzJkzFRUVpQMHDigwMDBH/4yMDHXp0kWBgYH68MMPVatWLR07dkyVK1cu+eIBAECp5NFwM2PGDA0aNEgDBw6UJM2bN09r167VokWLNHr06Bz9Fy1apPPnz+vLL7+Ut7e3JCk0NLQkSwYAAKWcx05LZWRkaNeuXYqMjPy1mDJlFBkZqe3bt+c65pNPPlF4eLhiYmIUFBSkpk2basqUKcrKysrzca5cuaLk5GSXGwAAsC6PhZuzZ88qKytLQUFBLu1BQUFKTEzMdczhw4f14YcfKisrS+vWrdOLL76o6dOn6+WXX87zcaZOnapKlSo5byEhIUX6PAAAQOni8QXF7nA4HAoMDNT8+fMVFham3r17a9y4cZo3b16eY8aMGaOkpCTn7ccffyzBigEAQEnz2JqbatWqycvLS6dOnXJpP3XqlGrUqJHrmODgYHl7e8vLy8vZ1qhRIyUmJiojI0M+Pj45xtjtdtnt9qItHgAAlFoeO3Lj4+OjsLAwbdq0ydnmcDi0adMmhYeH5zomIiJChw4dksPhcLYdPHhQwcHBuQYbAABw6/HoaakRI0ZowYIFWrp0qfbt26enn35aaWlpzqun+vXrpzFjxjj7P/300zp//ryGDRumgwcPau3atZoyZYpiYmI89RQAAEAp49FLwXv37q0zZ85o/PjxSkxMVMuWLbV+/XrnIuPjx4+rTJlf81dISIg2bNig4cOHq3nz5qpVq5aGDRumUaNGeeopAACAUsZmjDGeLqIkJScnq1KlSkpKSpK/v7+nywGAQknPyFTj8RskSXsnRamcj0f/VgWKnTu/v2+qq6UAAACuh3ADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAshXADAAAspaynCwAAuK+cT1kdndbd02UApRJHbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKWU9XQBJc0YI0lKTk72cCUAAKCgsn9vZ/8ez88tF25SUlIkSSEhIR6uBAAAuCslJUWVKlXKt4/NFCQCWYjD4dCJEydUsWJF2Ww2T5fjMcnJyQoJCdGPP/4of39/T5fjUcyFK+bDFfPhivn4FXPhqrjnwxijlJQU1axZU2XK5L+q5pY7clOmTBnddtttni6j1PD39+cf5f/HXLhiPlwxH66Yj18xF66Kcz6ud8QmGwuKAQCApRBuAACApRBublF2u12xsbGy2+2eLsXjmAtXzIcr5sMV8/Er5sJVaZqPW25BMQAAsDaO3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3FjY7NmzFRoaKl9fX7Vt21Y7d+7Mt//q1avVsGFD+fr6qlmzZlq3bl0JVVr83JmLBQsWqH379qpSpYqqVKmiyMjI687dzcbd10a2FStWyGaz6YEHHijeAkuYu/Nx8eJFxcTEKDg4WHa7XXfcccct++9FkmbOnKk777xTfn5+CgkJ0fDhw3X58uUSqrb4bNmyRdHR0apZs6ZsNps+/vjj646Ji4vTXXfdJbvdrgYNGmjJkiXFXmdJcXc+1qxZoy5duqh69ery9/dXeHi4NmzYUDLFGljSihUrjI+Pj1m0aJH57rvvzKBBg0zlypXNqVOncu2/bds24+XlZV599VWzd+9e89e//tV4e3ubPXv2lHDlRc/duXjkkUfM7NmzTXx8vNm3b58ZMGCAqVSpkvnpp59KuPLi4e58ZDty5IipVauWad++vfnjH/9YMsWWAHfn48qVK6Z169amW7duZuvWrebIkSMmLi7OJCQklHDlxcPd+Vi2bJmx2+1m2bJl5siRI2bDhg0mODjYDB8+vIQrL3rr1q0z48aNM2vWrDGSzEcffZRv/8OHD5ty5cqZESNGmL1795pZs2YZLy8vs379+pIpuJi5Ox/Dhg0zr7zyitm5c6c5ePCgGTNmjPH29jbffPNNsddKuLGoNm3amJiYGOd2VlaWqVmzppk6dWqu/Xv16mW6d+/u0ta2bVszZMiQYq2zJLg7F9fKzMw0FStWNEuXLi2uEktUYeYjMzPTtGvXzrz77rumf//+lgo37s7H3LlzTb169UxGRkZJlVii3J2PmJgYc++997q0jRgxwkRERBRrnSWtIL/MR44caZo0aeLS1rt3bxMVFVWMlXlGQeYjN40bNzYTJ04s+oKuwWkpC8rIyNCuXbsUGRnpbCtTpowiIyO1ffv2XMds377dpb8kRUVF5dn/ZlGYubhWenq6rl69qoCAgOIqs8QUdj4mTZqkwMBAPfHEEyVRZokpzHx88sknCg8PV0xMjIKCgtS0aVNNmTJFWVlZJVV2sSnMfLRr1067du1ynro6fPiw1q1bp27dupVIzaWJVd9Hi4rD4VBKSkqJvJfecl+ceSs4e/assrKyFBQU5NIeFBSk/fv35zomMTEx1/6JiYnFVmdJKMxcXGvUqFGqWbNmjjetm1Fh5mPr1q1auHChEhISSqDCklWY+Th8+LD+85//6NFHH9W6det06NAhPfPMM7p69apiY2NLouxiU5j5eOSRR3T27Fndc889MsYoMzNTTz31lMaOHVsSJZcqeb2PJicn69KlS/Lz8/NQZaXD66+/rtTUVPXq1avYH4sjN0A+pk2bphUrVuijjz6Sr6+vp8spcSkpKerbt68WLFigatWqebqcUsHhcCgwMFDz589XWFiYevfurXHjxmnevHmeLs0j4uLiNGXKFM2ZM0fffPON1qxZo7Vr1+qll17ydGkoRT744ANNnDhRq1atUmBgYLE/HkduLKhatWry8vLSqVOnXNpPnTqlGjVq5DqmRo0abvW/WRRmLrK9/vrrmjZtmj777DM1b968OMssMe7Oxw8//KCjR48qOjra2eZwOCRJZcuW1YEDB1S/fv3iLboYFeb1ERwcLG9vb3l5eTnbGjVqpMTERGVkZMjHx6dYay5OhZmPF198UX379tWTTz4pSWrWrJnS0tI0ePBgjRs3TmXK3Dp/Q+f1Purv739LH7VZsWKFnnzySa1evbrEjoDfOq+6W4iPj4/CwsK0adMmZ5vD4dCmTZsUHh6e65jw8HCX/pK0cePGPPvfLAozF5L06quv6qWXXtL69evVunXrkii1RLg7Hw0bNtSePXuUkJDgvPXo0UOdO3dWQkKCQkJCSrL8IleY10dERIQOHTrkDHmSdPDgQQUHB9/UwUYq3Hykp6fnCDDZwc/cYl9daNX30RuxfPlyDRw4UMuXL1f37t1L7oGLfckyPGLFihXGbrebJUuWmL1795rBgwebypUrm8TERGOMMX379jWjR4929t+2bZspW7asef31182+fftMbGyspS4Fd2cupk2bZnx8fMyHH35oTp486bylpKR46ikUKXfn41pWu1rK3fk4fvy4qVixohk6dKg5cOCA+ec//2kCAwPNyy+/7KmnUKTcnY/Y2FhTsWJFs3z5cnP48GHz73//29SvX9/06tXLU0+hyKSkpJj4+HgTHx9vJJkZM2aY+Ph4c+zYMWOMMaNHjzZ9+/Z19s++FPyFF14w+/btM7Nnz7bUpeDuzseyZctM2bJlzezZs13eSy9evFjstRJuLGzWrFmmdu3axsfHx7Rp08bs2LHDeV/Hjh1N//79XfqvWrXK3HHHHcbHx8c0adLErF27toQrLj7uzEWdOnWMpBy32NjYki+8mLj72vgtq4UbY9yfjy+//NK0bdvW2O12U69ePTN58mSTmZlZwlUXH3fm4+rVq2bChAmmfv36xtfX14SEhJhnnnnGXLhwoeQLL2KbN2/O9b0g+/n379/fdOzYMceYli1bGh8fH1OvXj2zePHiEq+7uLg7Hx07dsy3f3GyGXOLHTcEAACWxpobAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbwOISExM1bNgwNWjQQL6+vgoKClJERITmzp2r9PR0T5eXp9DQUNlsNq1YsSLHfU2aNJHNZtOSJUtKvrBcJCYm6tlnn1W9evVkt9sVEhKi6OjoHN8zJElTp06Vl5eXXnvtNWdb9nPN6zZgwABJyvP+3OYIuJXxreCAhR0+fFgRERGqXLmypkyZombNmslut2vPnj2aP3++atWqpR49euQ69urVq/L29i7hil2FhIRo8eLF+vOf/+xs27FjhxITE1W+fHkPVvaro0ePOuf4tddeU7NmzXT16lVt2LBBMTEx2r9/v0v/RYsWaeTIkVq0aJFeeOEFSdJXX32lrKwsSdKXX36phx56SAcOHJC/v78kuXyj9OLFi3X//fe77LNy5crF+AyBm1Cxf8EDAI+Jiooyt912m0lNTc31fofD4fx/SWbOnDkmOjralCtXzvldWnPmzDH16tUz3t7e5o477jDvvfeec8yRI0eMJBMfH+9su3DhgpFkNm/ebIz59fto/vnPf5pmzZoZu91u2rZte90vZa1Tp44ZPXq0sdvt5vjx4872QYMGmWeffdZUqlTJ5Xt7Lly4YJ544glTrVo1U7FiRdO5c2eTkJDgvP/QoUOmR48eJjAw0JQvX960bt3abNy4McdjTp482QwcONBUqFDBhISEmHfeeSffOrt27Wpq1aqV6xxf+/1KcXFxplatWiYjI8PUrFnTbNu2LceY7PnK7buZJJmPPvoo33oAGMNpKcCizp07p3//+9+KiYnJ8yiHzWZz2Z4wYYIefPBB7dmzR48//rg++ugjDRs2TH/5y1/07bffasiQIRo4cKA2b97sdj0vvPCCpk+frq+++krVq1dXdHS0rl69mu+YoKAgRUVFaenSpZKk9PR0rVy5Uo8//niOvg8//LBOnz6tf/3rX9q1a5fuuusu3XfffTp//rwkKTU1Vd26ddOmTZsUHx+v+++/X9HR0Tp+/LjLfqZPn67WrVsrPj5ezzzzjJ5++mkdOHAg1/rOnz+v9evX5znH1x5RWbhwofr06SNvb2/16dNHCxcuzPf5AygkT6crAMVjx44dRpJZs2aNS3vVqlVN+fLlTfny5c3IkSOd7ZLM//3f/7n0bdeunRk0aJBL28MPP2y6detmjHHvyM2KFSucfc6dO2f8/PzMypUr86y/Tp065o033jAff/yxqV+/vnE4HGbp0qWmVatWxhjjcuTmiy++MP7+/uby5csu+6hfv36+R16aNGliZs2a5fKYjz32mHPb4XCYwMBAM3fu3FzH//e//811jnOTlJRk/Pz8nEeT4uPjTYUKFUxKSopLv+sdufH19XX+/LJvx44du+7jA7cSjtwAt5idO3cqISFBTZo00ZUrV1zua926tcv2vn37FBER4dIWERGhffv2uf244eHhzv8PCAjQnXfeWaD9dO/eXampqdqyZYsWLVqU61Gb3bt3KzU1VVWrVlWFChWctyNHjuiHH36Q9MuRm+eff16NGjVS5cqVVaFCBe3bty/HkZvmzZs7/99ms6lGjRo6ffp0rrUZYwr03CVp+fLlql+/vlq0aCFJatmyperUqaOVK1cWeB+S9MYbbyghIcHlVrNmTbf2AVgdC4oBi2rQoIFsNluOUyr16tWT5LpINZu7i3TLlPnl76Pf/pK/3qkmd5UtW1Z9+/ZVbGys/vvf/+qjjz7K0Sc1NVXBwcGKi4vLcV/2qaHnn39eGzdu1Ouvv64GDRrIz89PPXv2VEZGhkv/axdR22w2ORyOXGu7/fbbZbPZciwazs3ChQv13XffqWzZX992HQ6HFi1apCeeeOK647PVqFFDDRo0KHB/4FbEkRvAoqpWraouXbro7bffVlpaWqH20ahRI23bts2lbdu2bWrcuLEkqXr16pKkkydPOu9PSEjIdV87duxw/v+FCxd08OBBNWrUqEB1PP744/r888/1xz/+UVWqVMlx/1133aXExESVLVtWDRo0cLlVq1bNWfeAAQP04IMPqlmzZqpRo4aOHj1aoMfPS0BAgKKiojR79uxc5/jixYuSpD179ujrr79WXFycyxGXuLg4bd++vUDhCEDBceQGsLA5c+YoIiJCrVu31oQJE9S8eXOVKVNGX331lfbv36+wsLB8x7/wwgvq1auXWrVqpcjISH366adas2aNPvvsM0m/HP353e9+p2nTpqlu3bo6ffq0/vrXv+a6r0mTJqlq1aoKCgrSuHHjVK1aNT3wwAMFeh6NGjXS2bNnVa5cuVzvj4yMVHh4uB544AG9+uqruuOOO3TixAmtXbtWDz74oFq3bq3bb79da9asUXR0tGw2m1588cU8j8i4Y/bs2YqIiFCbNm00adIkNW/eXJmZmdq4caPmzp2rffv2aeHChWrTpo06dOiQY/zdd9+thQsXunzuTX4uXryoxMREl7aKFSuWmkvjgVLB04t+ABSvEydOmKFDh5q6desab29vU6FCBdOmTRvz2muvmbS0NGc/5XGZcX6XghtjzN69e014eLjx8/MzLVu2NP/+979zXVD86aefmiZNmhgfHx/Tpk0bs3v37nzrzl5QnJdrLwVPTk42zz77rKlZs6bx9vY2ISEh5tFHH3VeRn7kyBHTuXNn4+fnZ0JCQszbb79tOnbsaIYNG5bvY7Zo0cJ5WXxeTpw4YWJiYkydOnWMj4+PqVWrlunRo4fZvHmzuXLliqlatap59dVXcx37yiuvmMDAQJORkWGMuf6C4txuU6dOzbc+4FZjM8aNFXEA4Ka4uDh17txZFy5c4MPmAJQI1twAAABLIdwAAABL4bQUAACwFI7cAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAAS/l/9UL239OF3uQAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "res_dml.cal.plot_cal(1)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi0AAAHHCAYAAABz3mgLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3NUlEQVR4nO3deXwV1f3/8fcFkhuyQghhkbAjS9jXhlQiFaUQQW0RRUVEvuASUFxAoL+K2K+NuGFlU5StrkAVq2KhSBNABAUMFmRHQCyEgEASEkgg9/z+4JtbLllvCLk58no+HvN4MOeemfncyXDvOzNnJg5jjBEAAEAlV8XXBQAAAJQGoQUAAFiB0AIAAKxAaAEAAFYgtAAAACsQWgAAgBUILQAAwAqEFgAAYAVCCwAAsAKhBYDPPPPMM3I4HB5tjRs31n333eeeX7BggRwOhzZt2lTB1XlKTk6Ww+FQcnKyT+sArmaEFhTJ4XCUaqqoD/FZs2ZpwYIFFbKty7V9+3Y988wzOnDgQIVve+PGjRo9erSio6MVFBSkhg0bavDgwdq9e/dlrzsjI0NTpkxRhw4dFBwcrOrVq6tt27Z66qmndPjw4XKo3vcq43F2/fXXe/yfq169utq3b69XX31VLper2GV//PFH1a9fXw6HQy+99FKR/Xbu3Knx48erY8eOCgkJUb169RQfH++zsOhyufTCCy+oSZMmCggIUPv27fX++++XevlTp05p1KhRql27toKCgtS7d299++23BfqdPn1aY8eOVYMGDeR0OtW6dWvNnj27PN8KylE1XxeAyuvtt9/2mP/rX/+qlStXFmhv3bp1hdQza9YsRUREePwWXllt375dU6ZM0fXXX6/GjRtX6LanTp2qdevW6fbbb1f79u2VmpqqGTNmqHPnztqwYYPatm1bpvX+8MMP6tOnj3788UfdfvvtGjVqlPz9/fXvf/9bc+fO1dKlS8slGO3atUtVqvju96mijrNevXrpzJkz8vf390ldDRo0UGJioiTp+PHjeu+99/TYY4/p2LFjeu655wpd5uTJk+rXr5+ysrLUq1cvjR8/XlFRUbrjjjsK9H3rrbc0d+5c/f73v9fDDz+s9PR0vfHGG/rVr36l5cuXq0+fPlf0/V3qD3/4g55//nmNHDlS3bp109///nfdddddcjgcuvPOO4td1uVyKT4+Xt99953GjRuniIgIzZo1S9dff702b96sFi1aSJLy8vLUt29fbdq0SQkJCWrRooVWrFihhx9+WCdPntSkSZMq4q3CGwYopYSEBFOaQyYrK+uKbD86OtrExcVdkXWXtyVLlhhJJikpqcK3vW7dOpOTk+PRtnv3buN0Os3dd99dpnWeO3fOdOjQwQQGBpq1a9cWeD09Pd1MmjTJ6/VOnjy5xGNq/vz5RpLZuHGj1+t3uVwmOzvbq2Uq43EWFxdnoqOjPdrOnDljGjVqZEJCQsz58+cLLHP27FnTq1cvExoaatavX2/Onj1r4uPjjdPpNKtXry7Qf9OmTSYzM9Oj7fjx46Z27domNja2fN9QCX766Sfj5+dnEhIS3G0ul8tcd911pkGDBoW+34stWrTISDJLlixxt6WlpZkaNWqYIUOGuNsWL15sJJm5c+d6LP/73//eBAQEmKNHj5bTO0J5IbSg1AoLLfkfpps2bTLXXXedqV69unn00UeNMRc+NJ9++mnTrFkz4+/vbxo0aGDGjRtnzp4967GOefPmmd69e5vatWsbf39/07p1azNr1iyPPo0aNTKSPKb8L5b8L7W1a9eaMWPGmIiICBMWFmZGjRplcnJyzMmTJ83QoUNNjRo1TI0aNcy4ceOMy+XyWH9eXp6ZNm2aadOmjXE6nSYyMtKMGjXKnDhxokAd8fHxZu3ataZbt27G6XSaJk2amIULF7r75Ndz6eSLAHOxzp07m86dO5dp2Q8++MBIMs8991yp+q9Zs8YMGjTIREVFuX/2Y8eOLRAgCgstjRo1MsOGDXPP5+/P1atXm1GjRpnw8HATEhJihg4dWuTPZ/ny5aZLly7G6XSaadOmGWMu/zhLSkoq9Oe4ePFi07lzZxMQEGBq1apl7r77bvPTTz959Bk2bJgJCgoyP/30k7nllltMUFCQiYiIME888USJX8DGFB5ajDFm0KBBRpI5fPiwR7vL5TJ33HGHCQsLM19//bW7PScnx9x8882mZs2a5vvvvy9xu8YY87vf/c6Eh4eXqm95mTlzppFUoMb33nvP/X+9OLfffrupU6eOycvL82gfNWqUCQwMdH8GjRkzxkgq8ItW/i8dc+bMKYd3g/LE5SFctp9//ln9+vXTnXfeqXvuuUd16tSRy+XSwIED9eWXX2rUqFFq3bq1tm7dqmnTpmn37t36+OOP3cvPnj1b0dHRGjhwoKpVq6ZPP/1UDz/8sFwulxISEiRJr776qsaMGaPg4GD94Q9/kCTVqVPHo44xY8aobt26mjJlijZs2KA5c+aoRo0a+uqrr9SwYUP9+c9/1ueff64XX3xRbdu21b333ute9oEHHtCCBQs0fPhwPfLII9q/f79mzJihlJQUrVu3Tn5+fu6+e/fu1aBBgzRixAgNGzZM8+bN03333acuXbooOjpavXr10iOPPKLXXntNkyZNcl8+K+4yWk5OjjIzM0u1vyMiIkrV72LGGB09elTR0dFeLytJn3zyiSRp6NChpeq/ZMkSZWdn66GHHlKtWrX0zTffaPr06frpp5+0ZMmSMtUwevRo1ahRQ88884x27dql2bNn6+DBg+4Bsvl27dqlIUOG6IEHHtDIkSPVsmVLSeV3nF0s/5jp1q2bEhMTdfToUf3lL3/RunXrlJKSoho1arj75l+K6NGjh1566SV98cUXevnll9WsWTM99NBDZdonBw4ckMPh8NiOJI0fP14rVqzQypUr1a1bN3e7v7+/PvzwQw0aNEj9+vXThg0bVK9evWK3kZqaWqpj7ty5c0pPTy9V3eHh4cVeAkxJSVFQUFCB/zPdu3d3v/7rX/+62OU7d+5cYBvdu3fXnDlztHv3brVr1045OTmqWrVqgUt+gYGBkqTNmzdr5MiRpXpPqCC+Tk2wR1FnWiSZ119/3aP97bffNlWqVCnwG9Hrr79uJJl169a52wo7fd+3b1/TtGlTj7aiTtvn/ybet29fjzMoMTExxuFwmAcffNDddv78edOgQQOP9axdu9ZIMu+++67HepcvX16gPf838TVr1rjb0tLSjNPpNE888YS7zdvLQ0WdnSlsKou333670NPgpdWpUycTFhZW6v6F/UwTExONw+EwBw8edLd5c6alS5cuJjc3193+wgsvGEnm73//u8eykszy5ctLVZM3x9mlZ1pyc3NNZGSkadu2rTlz5oy732effWYkmaefftrdNmzYMCPJPPvssx7r7NSpk+nSpUuBbV0qLi7OtGrVyhw7dswcO3bM7Ny504wbN85IMvHx8SUuX1Zr1qwxDofD/PGPfyyxb/7+Kc20f//+YtcVHx9f4OdizIVLz5LMhAkTil0+KCjI3H///QXaly1b5nF8vPzyy4WeuZkwYYKRZG6++eYS3jUqGmdacNmcTqeGDx/u0bZkyRK1bt1arVq10vHjx93tv/nNbyRJSUlJ6tmzpySpevXq7tfT09N17tw5xcXFacWKFUpPT1dYWFip6hgxYoTHb9w9evTQ+vXrNWLECHdb1apV1bVrV23evNmj1rCwMN14440etXbp0kXBwcFKSkrSXXfd5W5v06aNrrvuOvd87dq11bJlS/3www+lqrMwffv21cqVK8u8fHF27typhIQExcTEaNiwYWVaR0ZGhkJCQkrd/+KfaVZWls6cOaOePXvKGKOUlBQ1bNjQ6xpGjRrlccbroYce0qRJk/T5559r4MCB7vYmTZqob9++xdZ0OcdZvk2bNiktLU3PPPOMAgIC3O3x8fFq1aqVli1bpilTpngs8+CDD3rMX3fddQUGthdl586dql27tkfbwIEDNXfuXK/qLq20tDTdddddatKkicaPH19i/w4dOpT6GK5bt26xr585c0ZOp7NAe/5+PnPmTLksf9ddd+nZZ5/V/fffr5kzZ6pFixb65z//qVmzZpVqO6h4hBZctmuuuabA6dU9e/Zox44dBT5k86Wlpbn/vW7dOk2ePFnr169Xdna2Rz9vvkwu/SLMXy4qKqpA+8mTJz1qTU9PV2RkZIm1FrYdSapZs6bHOr1Vr169Ek/Tl0Vqaqri4+MVFhamv/3tb6patWqZ1hMaGupVKPvxxx/19NNP65NPPimwX0p7CeFS+Xd85AsODla9evUK3FbepEmTQpcvr+Ms38GDByXJffnpYq1atdKXX37p0RYQEFDg/4M3x03jxo315ptvyuVyad++fXruued07Ngxj8BUXrKysnTzzTcrMzNTX375pYKDg0tcpmbNmuV2h1H16tWVk5NToP3s2bPu18tj+bp16+qTTz7R0KFDddNNN0m6cKxPnz5dw4YNK9X7RsUitOCyFfYB4nK51K5dO73yyiuFLpMfJPbt26cbbrhBrVq10iuvvKKoqCj5+/vr888/17Rp00p8BsXFivpCLqzdGONRa2RkpN59991Cl7/0i6ao7Vy8Tm+dOXOm1F/mJf2Wmi89PV39+vXTqVOntHbtWtWvX7/M9bVq1UopKSk6dOhQgRB4qby8PN144406ceKEnnrqKbVq1UpBQUH6z3/+o/vuu8+rn2lZFHY8ludxVlZlDYz5goKCPEJBbGysOnfurEmTJum111673PLccnNz9bvf/U7//ve/tWLFilLfIp+bm6sTJ06Uqm/t2rWL3R/16tVTUlKSjDEeZ0+PHDkiSSUey/Xq1XP3vVhhy/fq1Us//PCDtm7dqqysLHXo0MH9zKFrr722VO8HFYfQgiuiWbNm+u6773TDDTcUeOLpxT799FPl5OTok08+8TiDkZSUVKBvceu53Fq/+OILxcbGlvgbXGl5W+uiRYsKXGIrSmnC0dmzZzVgwADt3r1bX3zxhdq0aeNVPZcaMGCA3n//fb3zzjuaOHFisX23bt2q3bt3a+HChR6DnS/38teePXvUu3dv9/zp06d15MgR9e/fv8Rlr8Rx1qhRI0kXBv7mX/bMt2vXLvfrV0r79u11zz336I033tCTTz5Zpktul3K5XLr33nu1atUqLV68WHFxcaVe9quvvvL4+RRn//79xT6/qGPHjnrrrbe0Y8cOj2P366+/dr9enI4dO2rt2rVyuVweg3G//vprBQYGFggjVatW9VjnF198IUkV/mwalIwn4uKKGDx4sP7zn//ozTffLPDamTNnlJWVJem/v31e/EWcnp6u+fPnF1guKChIp06duiK15uXl6U9/+lOB186fP1+mbQYFBUlSqZfNH9NSmqkkeXl5uuOOO7R+/XotWbJEMTExXtd/qUGDBqldu3Z67rnntH79+gKvZ2Zmuu+2KexnaozRX/7yl8uqYc6cOTp37px7fvbs2Tp//rz69etX4rJX4jjr2rWrIiMj9frrr3tcivjHP/6hHTt2KD4+vsR1XK7x48fr3LlzRZ7R9NaYMWO0aNEizZo1S7/73e+8WjZ/TEtpppLOFt5yyy3y8/Nzjy2RLvzsXn/9dV1zzTXu8XDShbMnO3fu9Dg2Bg0apKNHj+qjjz5ytx0/flxLlizRgAEDCh3vku/YsWOaOnWq2rdvT2iphDjTgiti6NChWrx4sR588EElJSUpNjZWeXl52rlzpxYvXqwVK1aoa9euuummm+Tv768BAwbogQce0OnTp/Xmm28qMjKywOndLl26aPbs2frf//1fNW/eXJGRkQV+wy2LuLg4PfDAA0pMTNSWLVt00003yc/PT3v27NGSJUv0l7/8RYMGDfJqnR07dlTVqlU1depUpaeny+l06je/+U2R42bKc0zLE088oU8++UQDBgzQiRMn9M4773i8fs8997j/nX/L7vz584t90rCfn58++ugj9enTR7169dLgwYMVGxsrPz8/ff/993rvvfdUs2ZNPffcc2rVqpWaNWumJ598Uv/5z38UGhqqDz/88LLG/EgXLj/ccMMNGjx4sHbt2qVZs2bp17/+tccg3KJciePMz89PU6dO1fDhwxUXF6chQ4a4b3lu3LixHnvssct6v6XRpk0b9e/fX2+99Zb++Mc/qlatWmVe16uvvqpZs2YpJiZGgYGBBY6b2267zR3GC1OeY1oaNGigsWPH6sUXX9S5c+fUrVs3ffzxx1q7dq3effddj0tLEydO1MKFCz3O3gwaNEi/+tWvNHz4cG3fvt39RNy8vLwCg6Pj4uIUExOj5s2bKzU1VXPmzNHp06f12Wef+fTJzCiCr25bgn2Ke7hcYXJzc83UqVNNdHS0cTqdpmbNmqZLly5mypQpJj093d3vk08+Me3btzcBAQGmcePGZurUqWbevHkFbo1MTU018fHxJiQkpNCHy136xNT822mPHTvm0Z7/oK9LzZkzx3Tp0sVUr17dhISEmHbt2pnx48d7PLgr/+Fll4qLiytwm+ybb75pmjZtaqpWrVqhD5fLvw29qOli06dPL/IW4cKcPHnSPP3006Zdu3YmMDDQBAQEmLZt25qJEyeaI0eOuPtt377d9OnTxwQHB5uIiAgzcuRI89133xlJZv78+e5+ZXm4XM2aNU1wcLC5++67zc8//1xg2aJuAb7c46yoh8stWrTIdOrUyTidThMeHl7sw+UuVZonAhtT/P+z5ORkI8lMnjy5xPUUJ/+27KKmkm5TLm95eXnmz3/+s2nUqJHx9/c30dHR5p133imy7kvrO3HihBkxYoSpVauWCQwMNHFxcYU+Vfmxxx4zTZs2NU6n09SuXdvcddddZt++fVfqbeEyOYy5jNGDAKw2ePBgHThwQN98842vSwGAEnF5CLhKGWOUnJxc4DIAAFRWnGkBAABWYJQRAACwAqEFAABYgdACAACsQGgBAABWsPruIZfLpcOHDyskJOSKPeIdAACUL2OMMjMzVb9+fa8e4md1aDl8+HCJf7wNAABUTocOHVKDBg1K3d/q0BISEiLpwpsODQ31cTUAAKA0MjIyFBUV5f4eLy2rQ0v+JaHQ0FBCCwAAlvF2aAcDcQEAgBUILQAAwAqEFgAAYAVCCwAAsAKhBQAAWIHQAgAArEBoAQAAViC0AAAAKxBaAACAFQgtAADACoQWAABgBUILAACwAqEFAABYodKElueff14Oh0Njx471dSkAAKASqhShZePGjXrjjTfUvn17X5cCAAAkZeeeV+MJy9R4wjJl5573dTmSKkFoOX36tO6++269+eabqlmzpq/LAQAAlZTPQ0tCQoLi4+PVp0+fEvvm5OQoIyPDYwIAAFeHar7c+AcffKBvv/1WGzduLFX/xMRETZky5QpXBQAAKiOfnWk5dOiQHn30Ub377rsKCAgo1TITJ05Uenq6ezp06NAVrhIAAFQWPjvTsnnzZqWlpalz587utry8PK1Zs0YzZsxQTk6Oqlat6rGM0+mU0+ms6FIBAEAl4LPQcsMNN2jr1q0ebcOHD1erVq301FNPFQgsAADg6uaz0BISEqK2bdt6tAUFBalWrVoF2gEAAHx+9xAAAEBp+PTuoUslJyf7ugQAAFBJcaYFAABYgdACAACsQGgBAABWILQAAAArEFoAAIAVCC0AAKBYB45n+7oESYQWAABwiVPZuRq5cLN7vv9ra3Xv3G+Unn3Oh1URWgAAwCUeeX+L1u877tG2bu9xjXk/xUcVXUBoAQAAbj8cO601e47JdUl7njFas+eY9h/P8kldEqEFAABc5OCJ4sevHPiZ0AIAACqBRuGBxb7euFZQBVVSEKEFAAC4Na0drF4tahcICFUdDvVqUVtNIggtAACgkpg+pJNimkV4tMU2j9D0IZ18VNEFhBYAAOAhLNBPbw7r4p7//JHr9NcR3RUW6OfDqggtAACgBI0jih/nUlEILQAAwAqEFgAAYAVCCwAAsAKhBQAAWIHQAgAArEBoAQAAViC0AAAAKxBaAACAFQgtAADACoQWAABgBUILAACwAqEFAABYgdACAACsUM3XBQAAgMon0L+aDjwf7+syPHCmBQAAWIHQAgAArEBoAQAAViC0AAAAKxBaAACAFQgtAADACoQWAABgBUILAACwAqEFAABYgdACAACsQGgBAABWILQAAAArEFoAAIAVCC0AAMAKhBYAAGAFQgsAALACoQUAAFiB0AIAAKxAaAEAAFYgtAAAACsQWgAAgBUILQAAwAqEFgAAYAVCCwAAsAKhBQAAWIHQAgAArEBoAQAAViC0AAAAKxBaAACAFQgtAADACoQWAABgBUILAACwAqEFAABYgdACAACsQGgBAABWILQAAAArEFoAAIAVCC0AAMAKhBYAAGAFQgsAALACoQUAAFiB0AIAAKzg09Aye/ZstW/fXqGhoQoNDVVMTIz+8Y9/+LIkAABQSfk0tDRo0EDPP/+8Nm/erE2bNuk3v/mNbrnlFn3//fe+LAsAAFRCDmOM8XURFwsPD9eLL76oESNGlNg3IyNDYWFhSk9PV2hoaAVUBwAALldZv7+rXcGavJKXl6clS5YoKytLMTExhfbJyclRTk6Oez4jI6OiygMAAD7m84G4W7duVXBwsJxOpx588EEtXbpUbdq0KbRvYmKiwsLC3FNUVFQFVwsAAHzF55eHcnNz9eOPPyo9PV1/+9vf9NZbb2n16tWFBpfCzrRERUVxeQgAAIuU9fKQz0PLpfr06aNmzZrpjTfeKLEvY1oAALBPWb+/fX556FIul8vjbAoAAIDk44G4EydOVL9+/dSwYUNlZmbqvffeU3JyslasWOHLsgAAQCXk09CSlpame++9V0eOHFFYWJjat2+vFStW6MYbb/RlWQAAoBLyaWiZO3euLzcPAAAsUunGtAAAABSG0AIAAKxAaAEAAFYgtAAAACsQWgAAgBUILQAAwAqEFgAAYAVCCwAAsAKhBQAAWIHQAgAArEBoAQAAViC0AAAAKxBaAACAFQgtAADACoQWAABgBUILAACwAqEFAABYgdACAACsQGgBAABWILQAAAArEFoAAIAVCC0AAMAKhBYAAGAFQgsAALACoQUAAFiB0AIAAKxAaAEAAFYgtAAAACsQWgAAgBUILQAAwAqEFgAAYAVCCwAAsAKhBQB8LDv3vBpPWKbGE5YpO/e8r8sBKi1CCwAAsAKhBQAAWIHQAgAArEBoAQAAViC0AAAAKxBaAACAFQgtAADACoQWAABgBUILAACwAqEFACqRA8ezfV0CUGkRWgDAh05l52rkws3u+f6vrdW9c79RevY5H1YFVE6EFgDwoUfe36L1+457tK3be1xj3k/xUUVA5UVoAQAf+eHYaa3Zc0yuS9rzjNGaPce0/3iWT+oCKitCCwD4yMETxY9fOfAzoQW4GKEFAHykUXhgsa83rhVUQZUAdiC0AICPNK0drF4tahf4IK7qcKhXi9pqEkFoAS5GaAEAH5o+pJNimkV4tMU2j9D0IZ18VBFQeRFaAMCHwgL99OawLu75zx+5Tn8d0V1hgX4+rAqonAgtAFCJNI4ofpwLcDUjtAAAACsQWgAAgBUILQAAwAqEFgAAYAVCCwAAsAKhBQAAWIHQAgAArEBoAQAAVvAqtCxevFi5ubnu+Z9++kku13//qHp2drZeeOGF8qsOAK4Cgf7VdOD5eB14Pl6B/tV8XQ5QaXkVWoYMGaJTp06559u0aaMDBw645zMzMzVx4sTyqg0AAMDNq9BijCl2HgAA4EphTAsAALACoQUAAFjB6xFfK1asUFhYmCTJ5XJp1apV2rZtmyR5jHcBAAAoTw7jxcCUKlVKPjHjcDiUl5d3WUWVVkZGhsLCwpSenq7Q0NAK2SYAALg8Zf3+9upMy8W3NwMAAFQkr8a03H///crMzLxStQAAABTJq9CycOFCnTlz5krVAgAAUKTLek4LAABARfH6lufMzExlZGQUO5VWYmKiunXrppCQEEVGRurWW2/Vrl27vC0JAABcBby+5fnaa68t8jVjjFd3D61evVoJCQnq1q2bzp8/r0mTJummm27S9u3bFRQU5G1pAADgF8zr0PK3v/1N4eHh5bLx5cuXe8wvWLBAkZGR2rx5s3r16lUu2wAAAL8MXoeW2NhYRUZGXolalJ6eLklFhqKcnBzl5OS45725FAUAAOxW7o/xL+uD5Vwul8aOHavY2Fi1bdu20D6JiYkKCwtzT1FRUZdTKgAAsIhXoaVRo0aqWrVqoa/t3r1b48ePV4MGDcpUSEJCgrZt26YPPvigyD4TJ05Uenq6ezp06FCZtgUAAOzj1eWh/fv3e8xnZ2dr0aJFmjdvntavX6+uXbvq8ccf97qI0aNH67PPPtOaNWuKDT1Op1NOp9Pr9QMAAPt5PaZFkjZs2KC33npLS5YsUcOGDbVjxw4lJSXpuuuu82o9xhiNGTNGS5cuVXJyspo0aVKWcgAAwFXAq8tDL7/8sqKjozVo0CDVrFlTa9as0datW+VwOFSrVi2vN56QkKB33nlH7733nkJCQpSamqrU1FSeugsAAArw6q88V6tWTU899ZSeffZZj7Etfn5++u6779SmTRvvNu5wFNo+f/583XfffSUuz195BgDAPmX9/vbqTMuf/vQnLVmyRE2aNNFTTz2lbdu2eV3oxYwxhU6lCSwAAODq4lVomThxonbv3q23335bqamp6tGjhzp06CBjjE6ePHmlagQAACjbc1ri4uK0cOFCHTlyRA8//LA6d+6sXr16qWfPnnrllVfKu0YAAADvxrQUZ9u2bZo7d67effddpaWllccqS8SYFgAA7FPW72+vbnk+c+aMVq1apZtvvlnShctFFz9Wv1q1atq3b583qwQAACgVr0LLwoULtWzZMndomTFjhqKjo1W9enVJ0q5du1S/fn099thj5V8pUMlk555Xm6dXSJK2P9tXgf5leuwRAKCUvBrT8u6772rUqFEebe+9956SkpKUlJSkF154QYsXLy7XAgEAACQvQ8vevXvVrl0793xAQICqVPnvKrp3767t27eXX3UAAAD/x6vz2adOnfIYw3Ls2DGP110ul8frAAAA5cWrMy0NGjQo9oFy//73v8v8V54BAACK41Vo6d+/v55++mmdPXu2wGtnzpzRlClTFB8fX27FAQAA5PPq8tCkSZO0ePFitWzZUqNHj9a1114r6cJdQzNmzND58+c1adKkK1IoAAC4unkVWurUqaOvvvpKDz30kCZMmKD859I5HA7deOONmjVrlurUqXNFCgUAAFc3rx8s0aRJEy1fvlwnTpzQ3r17JUnNmzdXeHh4uRcHAACQr8xPwwoPD1f37t3LsxYAAIAilekPJgIAAFQ0QgsAALACoQUAAFiB0AIAAKxAaAEAAFYgtAAAACsQWgAAgBUILQAAwAqEFgAAYAVCCwAAsAKhBQAAWIHQApSDA8ezfV0CAPziEVqAMjiVnauRCze75/u/tlb3zv1G6dnnfFgVAPyyEVqAMnjk/S1av++4R9u6vcc15v0UH1UEAL98hBbASz8cO601e47JdUl7njFas+eY9h/P8kldAPBLR2gBvHTwRPHjVw78TGgBgCuB0AJ4qVF4YLGvN64VVEGVAMDVhdACeKlp7WD1alG7wH+eqg6HerWorSYRhBYAuBIILUAZTB/SSTHNIjzaYptHaPqQTj6qCAB++QgtQBmEBfrpzWFd3POfP3Kd/jqiu8IC/XxYFQD8shFagHLQOKL4cS4AgMtHaAEAAFYgtAAAACsQWgAAgBUILQAAwAqEFgAAYAVCCwAAsAKhBQAAWIHQAgAArEBoAQAAViC0AAAAKxBaAACAFQgtAADACoQWAABgBUILAACwAqEFAABYgdACAACsQGgBAABWILQAAAArEFoAAIAVCC0AAMAK1XxdAGCrQP9qOvB8vK/LAICrBmdaAACAFQgtAADACoQWAABgBUILAACwAqEFAABYgdACAACsQGgBAABWILQAAAArEFoAAIAVCC0AAMAKhBYAAGAFQgsAALACoQUAAFiB0AIAAKxAaAEAAFbwaWhZs2aNBgwYoPr168vhcOjjjz/2ZTkAAKAS82loycrKUocOHTRz5kxflgEAACxQzZcb79evn/r16+fLEgAAgCV8Glq8lZOTo5ycHPd8RkaGD6sBAAAVyaqBuImJiQoLC3NPUVFRvi4JAABUEKtCy8SJE5Wenu6eDh065OuSAABABbHq8pDT6ZTT6fR1GQAAwAesOtMCAACuXj4903L69Gnt3bvXPb9//35t2bJF4eHhatiwoQ8rAwAAlY1PQ8umTZvUu3dv9/zjjz8uSRo2bJgWLFjgo6oAAEBl5NPQcv3118sY48sSAACAJRjTAgAArEBoAQAAViC0AAAAKxBaAACAFQgtAADACoQWAABgBUILAACwAqEFAABYgdACAACsQGgBAABWILQAAAArEFoAAIAVCC0AAMAKhBYAAGAFQgsAALACoQUAAFiB0AIAAKxAaAEAAFYgtBQiO/e8Gk9YpsYTlik797yvywEAACK0AAAASxBaAACAFQgtAADACoQWAABgBUILAACwAqEFAABYgdACAACsQGgBAABWILQAAAArEFoAAIAVCC0AAMAKhBYAAGAFQgsAALACoQUAAFiB0AIAAKxAaAEAAFYgtAAAACsQWgAAgBUILQAAwAqEFgAAYAVCCwAAsAKhBQAAWIHQAgAArEBoAQAAViC0AAAAKxBaAACAFQgtJThwPNvXJQAAABFaCjiVnauRCze75/u/tlb3zv1G6dnnfFgVAAAgtFzikfe3aP2+4x5t6/Ye15j3U3xUEQAAkAgtHn44dlpr9hyT65L2PGO0Zs8x7T+e5ZO6AAAAocXDwRPFj1858DOhBQAAXyG0XKRReGCxrzeuFVRBlQAAgEsRWi7StHawerWoXWCnVHU41KtFbTWJILQAAOArhJZLTB/SSTHNIjzaYptHaPqQTj6qCAAASISWAsIC/fTmsC7u+c8fuU5/HdFdYYF+PqwKAAAQWkrQOKL4cS4AAKBiEFoAAIAVCC0AAMAKhBYAAGAFQgsAALACoQUAAFiB0AIAAKxAaAEAAFYgtAAAACsQWgAAgBUILQAAwAqEFgAAYAVCCwAAsAKhBQAAWIHQAgAArEBoAQAAViC0AAAAK1SK0DJz5kw1btxYAQEB6tGjh7755htflwQAACoZn4eWRYsW6fHHH9fkyZP17bffqkOHDurbt6/S0tJ8XRoAAKhEfB5aXnnlFY0cOVLDhw9XmzZt9PrrryswMFDz5s3zdWkAAKAS8Wloyc3N1ebNm9WnTx93W5UqVdSnTx+tX7++QP+cnBxlZGR4TAAA4Org09By/Phx5eXlqU6dOh7tderUUWpqaoH+iYmJCgsLc09RUVEVVSoAAPAxn18e8sbEiROVnp7ung4dOuTrkgAAQAWp5suNR0REqGrVqjp69KhH+9GjR1W3bt0C/Z1Op5xOZ0WVBwAAKhGfnmnx9/dXly5dtGrVKneby+XSqlWrFBMT48PKAABAZePTMy2S9Pjjj2vYsGHq2rWrunfvrldffVVZWVkaPny4r0sDAACViM9Dyx133KFjx47p6aefVmpqqjp27Kjly5cXGJxbkQL9q+nA8/E+2z4AACjIYYwxvi6irDIyMhQWFqb09HSFhob6uhwAAFAKZf3+turuIQAAcPUitAAAACsQWgAAgBUILQAAwAqEFgAAYAVCCwAAsAKhBQAAWIHQAgAArEBoAQAAViC0AAAAKxBaAACAFQgtAADACoQWAABgBUILAACwQjVfF3A5jDGSLvyJawAAYIf87+387/HSsjq0ZGZmSpKioqJ8XAkAAPBWZmamwsLCSt3fYbyNOZWIy+XS4cOHFRISIofD4etyfCIjI0NRUVE6dOiQQkNDfV3OLwL7tPyxT8sf+7R8sT/LX3H71BijzMxM1a9fX1WqlH6kitVnWqpUqaIGDRr4uoxKITQ0lP9o5Yx9Wv7Yp+WPfVq+2J/lr6h96s0ZlnwMxAUAAFYgtAAAACsQWizndDo1efJkOZ1OX5fyi8E+LX/s0/LHPi1f7M/ydyX2qdUDcQEAwNWDMy0AAMAKhBYAAGAFQgsAALACoQUAAFiB0GKBmTNnqnHjxgoICFCPHj30zTffFNt/yZIlatWqlQICAtSuXTt9/vnnFVSpPbzZpwsWLJDD4fCYAgICKrDaym3NmjUaMGCA6tevL4fDoY8//rjEZZKTk9W5c2c5nU41b95cCxYsuOJ12sTbfZqcnFzgGHU4HEpNTa2Ygiu5xMREdevWTSEhIYqMjNStt96qXbt2lbgcn6VFK8s+LY/PUkJLJbdo0SI9/vjjmjx5sr799lt16NBBffv2VVpaWqH9v/rqKw0ZMkQjRoxQSkqKbr31Vt16663atm1bBVdeeXm7T6ULT3Q8cuSIezp48GAFVly5ZWVlqUOHDpo5c2ap+u/fv1/x8fHq3bu3tmzZorFjx+p//ud/tGLFiitcqT283af5du3a5XGcRkZGXqEK7bJ69WolJCRow4YNWrlypc6dO6ebbrpJWVlZRS7DZ2nxyrJPpXL4LDWo1Lp3724SEhLc83l5eaZ+/fomMTGx0P6DBw828fHxHm09evQwDzzwwBWt0ybe7tP58+ebsLCwCqrObpLM0qVLi+0zfvx4Ex0d7dF2xx13mL59+17ByuxVmn2alJRkJJmTJ09WSE22S0tLM5LM6tWri+zDZ6l3SrNPy+OzlDMtlVhubq42b96sPn36uNuqVKmiPn36aP369YUus379eo/+ktS3b98i+19tyrJPJen06dNq1KiRoqKidMstt+j777+viHJ/kThGr5yOHTuqXr16uvHGG7Vu3Tpfl1NppaenS5LCw8OL7MNx6p3S7FPp8j9LCS2V2PHjx5WXl6c6dep4tNepU6fIa9Wpqale9b/alGWftmzZUvPmzdPf//53vfPOO3K5XOrZs6d++umniij5F6eoYzQjI0NnzpzxUVV2q1evnl5//XV9+OGH+vDDDxUVFaXrr79e3377ra9Lq3RcLpfGjh2r2NhYtW3btsh+fJaWXmn3aXl8llr9V56BihATE6OYmBj3fM+ePdW6dWu98cYb+tOf/uTDyoALWrZsqZYtW7rne/bsqX379mnatGl6++23fVhZ5ZOQkKBt27bpyy+/9HUpvxil3afl8VnKmZZKLCIiQlWrVtXRo0c92o8ePaq6desWukzdunW96n+1Kcs+vZSfn586deqkvXv3XokSf/GKOkZDQ0NVvXp1H1X1y9O9e3eO0UuMHj1an332mZKSktSgQYNi+/JZWjre7NNLleWzlNBSifn7+6tLly5atWqVu83lcmnVqlUeafViMTExHv0laeXKlUX2v9qUZZ9eKi8vT1u3blW9evWuVJm/aByjFWPLli0co//HGKPRo0dr6dKl+te//qUmTZqUuAzHafHKsk8vVabP0ssaxosr7oMPPjBOp9MsWLDAbN++3YwaNcrUqFHDpKamGmOMGTp0qJkwYYK7/7p160y1atXMSy+9ZHbs2GEmT55s/Pz8zNatW331Fiodb/fplClTzIoVK8y+ffvM5s2bzZ133mkCAgLM999/76u3UKlkZmaalJQUk5KSYiSZV155xaSkpJiDBw8aY4yZMGGCGTp0qLv/Dz/8YAIDA824cePMjh07zMyZM03VqlXN8uXLffUWKh1v9+m0adPMxx9/bPbs2WO2bt1qHn30UVOlShXzxRdf+OotVCoPPfSQCQsLM8nJyebIkSPuKTs7292Hz1LvlGWflsdnKaHFAtOnTzcNGzY0/v7+pnv37mbDhg3u1+Li4sywYcM8+i9evNhce+21xt/f30RHR5tly5ZVcMWVnzf7dOzYse6+derUMf379zfffvutD6qunPJvt710yt+Hw4YNM3FxcQWW6dixo/H39zdNmzY18+fPr/C6KzNv9+nUqVNNs2bNTEBAgAkPDzfXX3+9+de//uWb4iuhwvalJI/jjs9S75Rln5bHZ6nj/zYOAABQqTGmBQAAWIHQAgAArEBoAQAAViC0AAAAKxBaAACAFQgtAADACoQWAABgBUILAACwAqEFsEhqaqoeffRRNW/eXAEBAapTp45iY2M1e/ZsZWdn+7q8IjVu3FgOh0MffPBBgdeio6PlcDi0YMGCii+sEKmpqRozZoyaNm0qp9OpqKgoDRgwoMDfoZGkxMREVa1aVS+++KK7Lf+9FjXdd999klTk64XtIwAXVPN1AQBK54cfflBsbKxq1KihP//5z2rXrp2cTqe2bt2qOXPm6JprrtHAgQMLXfbcuXPy8/Or4Io9RUVFaf78+brzzjvdbRs2bFBqaqqCgoJ8WNl/HThwwL2PX3zxRbVr107nzp3TihUrlJCQoJ07d3r0nzdvnsaPH6958+Zp3LhxkqSNGzcqLy9PkvTVV1/p97//vXbt2qXQ0FBJ8vhL1vPnz9dvf/tbj3XWqFHjCr5DwHLl8kcIAFxxffv2NQ0aNDCnT58u9HWXy+X+tyQza9YsM2DAABMYGGgmT55sjDFm1qxZpmnTpsbPz89ce+215q9//at7mf379xtJJiUlxd128uRJI8kkJSUZY/77N3E+++wz065dO+N0Ok2PHj1K/CNyjRo1MhMmTDBOp9P8+OOP7vaRI0eaMWPGmLCwMI+/WXLy5EkzYsQIExERYUJCQkzv3r3Nli1b3K/v3bvXDBw40ERGRpqgoCDTtWtXs3LlygLbfO6558zw4cNNcHCwiYqKMm+88Uaxdfbr189cc801he7jkydPeswnJyeba665xuTm5pr69eubdevWFVgmf39duqwxF35GS5cuLbYeAJ64PARY4Oeff9Y///lPJSQkFHlWwuFweMw/88wzuu2227R161bdf//9Wrp0qR599FE98cQT2rZtmx544AENHz5cSUlJXtczbtw4vfzyy9q4caNq166tAQMG6Ny5c8UuU6dOHfXt21cLFy6UJGVnZ2vRokW6//77C/S9/fbblZaWpn/84x/avHmzOnfurBtuuEEnTpyQJJ0+fVr9+/fXqlWrlJKSot/+9rcaMGCAfvzxR4/1vPzyy+ratatSUlL08MMP66GHHtKuXbsKre/EiRNavnx5kfv40jMgc+fO1ZAhQ+Tn56chQ4Zo7ty5xb5/AOXA16kJQMk2bNhgJJmPPvrIo71WrVomKCjIBAUFmfHjx7vbJZmxY8d69O3Zs6cZOXKkR9vtt99u+vfvb4zx7kzLBx984O7z888/m+rVq5tFixYVWX+jRo3MtGnTzMcff2yaNWtmXC6XWbhwoenUqZMxxnicaVm7dq0JDQ01Z8+e9VhHs2bNij1TEh0dbaZPn+6xzXvuucc973K5TGRkpJk9e3ahy3/99deF7uPCpKenm+rVq7vP/qSkpJjg4GCTmZnp0a+kMy0BAQHun1/+dPDgwRK3D1ytONMCWOybb77Rli1bFB0drZycHI/Xunbt6jG/Y8cOxcbGerTFxsZqx44dXm83JibG/e/w8HC1bNmyVOuJj4/X6dOntWbNGs2bN6/QsyzfffedTp8+rVq1aik4ONg97d+/X/v27ZN04UzLk08+qdatW6tGjRoKDg7Wjh07Cpxpad++vfvfDodDdevWVVpaWqG1GS/+4P3777+vZs2aqUOHDpKkjh07qlGjRlq0aFGp1yFJ06ZN05YtWzym+vXre7UO4GrCQFzAAs2bN5fD4ShwaaNp06aSPAd35vN2cGuVKhd+h7n4y7ukSz7eqlatmoYOHarJkyfr66+/1tKlSwv0OX36tOrVq6fk5OQCr+VfonnyySe1cuVKvfTSS2revLmqV6+uQYMGKTc316P/pYOPHQ6HXC5XobW1aNFCDoejwGDbwsydO1fff/+9qlX770eoy+XSvHnzNGLEiBKXz1e3bl01b9681P2Bqx1nWgAL1KpVSzfeeKNmzJihrKysMq2jdevWWrdunUfbunXr1KZNG0lS7dq1JUlHjhxxv75ly5ZC17Vhwwb3v0+ePKndu3erdevWparj/vvv1+rVq3XLLbeoZs2aBV7v3LmzUlNTVa1aNTVv3txjioiIcNd933336bbbblO7du1Ut25dHThwoFTbL0p4eLj69u2rmTNnFrqPT506JUnaunWrNm3apOTkZI8zJMnJyVq/fn2pQg+AsuFMC2CJWbNmKTY2Vl27dtUzzzyj9u3bq0qVKtq4caN27typLl26FLv8uHHjNHjwYHXq1El9+vTRp59+qo8++khffPGFpAtna371q1/p+eefV5MmTZSWlqb/9//+X6HrevbZZ1WrVi3VqVNHf/jDHxQREaFbb721VO+jdevWOn78uAIDAwt9vU+fPoqJidGtt96qF154Qddee60OHz6sZcuW6bbbblPXrl3VokULffTRRxowYIAcDof++Mc/FnkGxRszZ85UbGysunfvrmeffVbt27fX+fPntXLlSs2ePVs7duzQ3Llz1b17d/Xq1avA8t26ddPcuXM9nttSnFOnTik1NdWjLSQkpNLcAg5UOr4eVAOg9A4fPmxGjx5tmjRpYvz8/ExwcLDp3r27efHFF01WVpa7n4q4nba4W56NMWb79u0mJibGVK9e3XTs2NH885//LHQg7qeffmqio6ONv7+/6d69u/nuu++KrTt/IG5RLr3lOSMjw4wZM8bUr1/f+Pn5maioKHP33Xe7b5fev3+/6d27t6levbqJiooyM2bMMHFxcebRRx8tdpsdOnRw3/5dlMOHD5uEhATTqFEj4+/vb6655hozcOBAk5SUZHJyckytWrXMCy+8UOiyU6dONZGRkSY3N9cYU/JA3MKmxMTEYusDrmYOY7wYfQbgqpacnKzevXvr5MmTPAQNQIVjTAsAALACoQUAAFiBy0MAAMAKnGkBAABWILQAAAArEFoAAIAVCC0AAMAKhBYAAGAFQgsAALACoQUAAFiB0AIAAKxAaAEAAFb4/2z23HHxLl6gAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "res_dml.cal.plot_cal(2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### T-Learner" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + }, + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
treatmentblp_estblp_seblp_pvalqini_estqini_seqini_pvalcal_r_squared
01-0.1850.1110.096-0.0440.0220.023-2.747
120.7160.0600.0000.3710.0250.0000.626
\n", + "
" + ], + "text/plain": [ + " treatment blp_est blp_se blp_pval qini_est qini_se qini_pval \\\n", + "0 1 -0.185 0.111 0.096 -0.044 0.022 0.023 \n", + "1 2 0.716 0.060 0.000 0.371 0.025 0.000 \n", + "\n", + " cal_r_squared \n", + "0 -2.747 \n", + "1 0.626 " + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Initialize DRtester and fit/predict nuisance models\n", + "t_tester = DRtester(\n", + " model_regression=model_regression, \n", + " model_propensity=model_propensity, \n", + " cate=est_t\n", + ").fit_nuisance(Xval, Dval, Yval, Xtrain, Dtrain, Ytrain)\n", + "\n", + "res_t = t_tester.evaluate_all(Xval, Xtrain)\n", + "res_t.summary()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHHCAYAAABDUnkqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABAo0lEQVR4nO3deVxUZf//8feoMOACigpuCC65kSJqGnq73WGk3ljdd2m2aGZmZWVZuVWa9Su1xWxxyXKtNLXM6qu3ZiaZppkGZbkrLqW4y6qgzPX7owdzO7LIKDBwfD0fj3nkuc51znzOxWnmzdmwGWOMAAAALKKMpwsAAAAoTIQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbAABgKYQbACXKiy++KJvN5tIWGhqq+++/3zk9Z84c2Ww2bd68uZircxUbGyubzabY2FiP1gHAFeEGbrHZbAV6FdeH/dSpUzVnzpxiea+rtW3bNr344ovav39/sb93amqqxo4dq1tuuUUBAQGy2WyFNm7JyckaN26cwsPDVbFiRfn6+ur666/XiBEjdPjw4UJ5D08riftZly5dXP6f8/X1VYsWLTR58mQ5HI58lz148KBq1aolm82mN954I89+O3bs0PDhw9WyZUtVqlRJNWvWVM+ePT0eKi/2888/67HHHlNYWJgqVKigunXrqnfv3tq1a1eBlr90HC9+eXl55bnc3r175ePjk2vIvtJ1ovCU83QBKF0++ugjl+l58+Zp1apVOdqbNm1aLPVMnTpV1apVc/mtvqTatm2bxo0bpy5duig0NLRY3/vEiRN66aWXVLduXYWHhxda+Ny3b5+ioqJ08OBB3XnnnXrooYfk7e2t3377TTNnztQXX3xR4C+Z/OzcuVNlynjud7G89rNOnTrp7Nmz8vb29khdderU0fjx4yX9/TOeP3++nnrqKR0/flyvvPJKrsucPn1a3bt3V1pamjp16qThw4crODhYffr0ydH3ww8/1MyZM/Wf//xHjz76qJKSkvT+++/rxhtv1IoVKxQVFVWk21cQEydO1Pr163XnnXeqRYsWSkxM1HvvvadWrVpp48aNuv766/Nd/rnnntODDz7o0paWlqaHH35YN998c57LPfXUUypXrpwyMjIKbZ0oRAa4CkOGDDEF2Y3S0tKK5P3DwsJM586di2TdhW3x4sVGklmzZk2xv/e5c+fMkSNHjDHG/Pzzz0aSmT179lWt8/z58yY8PNyUL1/e/PDDDznmJyUlmdGjR7u93rFjx152n5o9e7aRZH7++We31+9wOEx6erpby5TE/axz584mLCzMpe3s2bMmJCTEVKpUyVy4cCHHMufOnTOdOnUyfn5+ZsOGDebcuXOmZ8+exm63m++//z5H/82bN5uUlBSXthMnTpjq1aubDh06FO4GXaH169ebjIwMl7Zdu3YZu91u7rnnnita50cffWQkmU8++STX+StWrDDe3t7m+eefL/B+eLl1onBxWgqFrkuXLrr++uu1ZcsWderUSeXLl9fo0aMlSRkZGRo7dqwaNmwou92u4OBgDR8+PMdvP7Nnz9Y///lPBQYGym63q1mzZpo2bZpLn9DQUP3xxx/6/vvvnYd8u3TpIul/12SsW7dOTzzxhKpXr67KlStr8ODByszM1JkzZ9SvXz9VqVJFVapU0fDhw2WMcVm/w+HQ5MmTFRYWJh8fHwUFBWnw4ME6ffp0jjr+9a9/ad26dWrbtq18fHxUv359zZs3z9lnzpw5uvPOOyVJXbt2LfbTd3a7XTVq1CjUdX7++ef69ddf9dxzz+kf//hHjvl+fn4uRw9++OEH3Xnnnapbt67zZ//UU0/p7Nmzl32vS6+5yZaenq7BgweratWq8vPzU79+/fL8+axcuVJt2rSRr6+v3n//fUlXv5/ldc3N4sWL1bp1a/n6+qpatWq699579ddff7n0uf/++1WxYkX99ddfuu2221SxYkVVr15dzzzzjLKysi47Jrnx8fHRDTfcoJSUFB07dsxlnjFG/fv316+//qpVq1bpxhtvlN1u15IlS9StWzfddttt2rZtm8syrVu3VsWKFV3aqlatqo4dO2r79u1XVGNha9++fY4jZ9ddd53CwsKuuMb58+erQoUKuvXWW3PMO3/+vIYOHaqhQ4eqQYMGhbJOFD5OS6FInDx5Ut27d9ddd92le++9V0FBQXI4HOrVq5fWrVunhx56SE2bNtXWrVv11ltvadeuXVq6dKlz+WnTpiksLEy9evVSuXLl9PXXX+vRRx+Vw+HQkCFDJEmTJ0/W448/rooVK+q5556TJAUFBbnU8fjjj6tGjRoaN26cNm7cqBkzZqhy5cr68ccfVbduXb366qtavny5Xn/9dV1//fXq16+fc9nBgwdrzpw5GjBggJ544gklJCTovffeU1xcnNavX+9y7nzPnj264447NHDgQPXv31+zZs3S/fffr9atWyssLEydOnXSE088oXfeeUejR492nrbL7/RdRkaGUlJSCjTe1apVK1C/wvTVV19Jku67774C9V+8eLHS09P1yCOPqGrVqtq0aZPeffdd/fnnn1q8ePEV1fDYY4+pcuXKevHFF7Vz505NmzZNBw4ccIaObDt37lTfvn01ePBgDRo0SI0bN5ZUePvZxbL3mRtuuEHjx4/X0aNH9fbbb2v9+vWKi4tT5cqVnX2zsrIUHR2tdu3a6Y033tC3336rN998Uw0aNNAjjzxyRWOyf/9+2Ww2l/eRpOHDh2vlypVatWqVbrjhBme7t7e3Pv/8c91xxx3q3r27Nm7cqJo1a+b7HomJiQXa586fP6+kpKQC1R0QEFBopx6NMTp69KjCwsLcXvb48eNatWqV+vTpowoVKuSYP3nyZJ0+fVrPP/+8lixZUijrRBHw8JEjlHK5nZbq3LmzkWSmT5/u0v7RRx+ZMmXK5DiFMX36dCPJrF+/3tmW22mD6OhoU79+fZe2vE4XZJ+2iI6ONg6Hw9keGRlpbDabefjhh51tFy5cMHXq1HFZzw8//JDrIeQVK1bkaA8JCTGSzNq1a51tx44dM3a73Tz99NPONndPS2VvQ0Fe7iis01IRERHG39+/wP1z+5mOHz/e2Gw2c+DAAWdbbqelQkJCTP/+/Z3T2WPTunVrk5mZ6Wx/7bXXjCTz5ZdfuiwryaxYsaJANbmzn61Zs8blZ5qZmWkCAwPN9ddfb86ePevs93//939GkhkzZoyzrX///kaSeemll1zWGRERYVq3bp3jvS7VuXNn06RJE3P8+HFz/Phxs2PHDvPss88aSaZnz56XXf5KrV271thsNvPCCy9ctm/2+BTklZCQUGg1Zp8CmjlzptvLvvvuu0aSWb58eY55R44cMZUqVTLvv/++Mabgp0fzWyeKBkduUCTsdrsGDBjg0rZ48WI1bdpUTZo00YkTJ5zt//znPyVJa9asUfv27SVJvr6+zvlJSUk6f/68OnfurJUrVyopKUn+/v4FqmPgwIEuv8G3a9dOGzZs0MCBA51tZcuWVZs2bbRlyxaXWv39/dWtWzeXWrMP069Zs0Z33323s71Zs2bq2LGjc7p69epq3Lix9u3bV6A6cxMdHa1Vq1Zd8fJFLTk5WZUqVSpw/4t/pmlpaTp79qzat28vY4zi4uJUt25dt2t46KGHXI6gPfLIIxo9erSWL1+uXr16Odvr1aun6OjofGu6mv0s2+bNm3Xs2DG9+OKL8vHxcbb37NlTTZo00bJlyzRu3DiXZR5++GGX6Y4dO+a4QD8vO3bsUPXq1V3aevXqpZkzZ7pVd0EdO3ZMd999t+rVq6fhw4dftn94eHiB9+HCOm26Y8cODRkyRJGRkerfv7/by8+fP1/Vq1dXt27dcswbMWKE6tevn+Ni4atZJ4oG4QZFonbt2jnOg+/evVvbt2/P8WGc7eJrBNavX6+xY8dqw4YNSk9Pd+nnzpfOpV+Y2csFBwfnaL/4Wo3du3crKSlJgYGBl601t/eRpCpVquS4/sMdNWvWvOzpAU/y8/NzK7wdPHhQY8aM0VdffZVjXAp66uJS1113nct0xYoVVbNmzRy329erVy/X5QtrP8t24MABSXKe9rpYkyZNtG7dOpc2Hx+fHP8/uLPfhIaG6oMPPpDD4dDevXv1yiuv6Pjx4y7BqrCkpaXpX//6l1JSUrRu3boc1+LkpkqVKoVyR1VWVpaOHz/u0hYQEJDjMyYxMVE9e/aUv7+/PvvsM5UtW9at99m3b582bNigxx57TOXKuX49bty4UR999JFWr17t1umz/NaJosNIo0hc/BtxNofDoebNm2vSpEm5LpMdOPbu3aubbrpJTZo00aRJkxQcHCxvb28tX75cb7311mWf4XGxvD7ccms3F11Q7HA4FBgYqE8++STX5S/9QsrrfcwlFym74+zZswX+0i/si4ULokmTJoqLi9OhQ4dyhMVLZWVlqVu3bjp16pRGjBihJk2aqEKFCvrrr790//33u/UzvRK57Y+FuZ9dKXe/fC9VoUIFl/DQoUMHtWrVSqNHj9Y777xzteU5ZWZm6t///rd+++03rVy58rK3V1+83KlTpwrUt3r16nmOx6FDh3IE1DVr1jgv7Jb+DqPdu3fXmTNn9MMPP6hWrVoFet+LzZ8/X5J0zz335Jg3fPhwdezYUfXq1XOG5+yjukeOHNHBgwdz/SUnv3Wi6BBuUGwaNGigX3/9VTfddFOOJ9Be7Ouvv1ZGRoa++uorlw+LNWvW5Oib33quttZvv/1WHTp0yPWL8Uq4W+vChQtznNrLy9WEqCsVExOjBQsW6OOPP9aoUaPy7bt161bt2rVLc+fOdblo+2pPu+3evVtdu3Z1TqempurIkSPq0aPHZZctiv0sJCRE0t8XMGefbs22c+dO5/yi0qJFC9177716//339cwzz1zRqb5LORwO9evXT6tXr9aiRYvUuXPnAi/7448/uvx88pOQkJDn859q1KiRY18JDw93/vvcuXOKiYnRrl279O2336pZs2YFrvFi8+fPV4MGDXTjjTfmmHfw4EEdOHAg16OAvXr1kr+/v86cOePWOlF0CDcoNr1799by5cv1wQcf6KGHHnKZd/bsWTkcDlWoUMH529vFX9hJSUmaPXt2jnVWqFAh1w+Uwqh16tSpevnll/Xqq6+6zLtw4YJSU1Nz3I1yOdl3SRS03pJ+zc0dd9yh8ePH65VXXlGXLl0UGRnpMj8lJUUTJkzQK6+8kuvP1Bijt99++6pqmDFjhgYMGOC87mbatGm6cOGCunfvftlli2I/a9OmjQIDAzV9+nQ98MADstvtkqT//ve/2r59u8aMGVOQzboqw4cP17x58zRp0iRNnjz5qtf3+OOPa+HChXr//ff173//261lC+uaGx8fnzxPb2VlZalPnz7asGGDvvzyyxz74cWOHDmipKQkNWjQIMeTguPi4rR9+3a98MILuS47Y8aMHKcuv/vuO7377rt644031KRJkxzLXG6dKDqEGxSb++67T4sWLdLDDz+sNWvWqEOHDsrKytKOHTu0aNEi53NIbr75Znl7eysmJkaDBw9WamqqPvjgAwUGBurIkSMu62zdurWmTZum//f//p8aNmyowMDAHL8xX4nOnTtr8ODBGj9+vOLj43XzzTfLy8tLu3fv1uLFi/X222/rjjvucGudLVu2VNmyZTVx4kQlJSXJbrc7n7GSm8K+5ua9997TmTNnnH8S4euvv9aff/4p6e8vsOzrS7JvZZ49e3a+T3728vLSkiVLFBUVpU6dOql3797q0KGDvLy89Mcff2j+/PmqUqWKXnnlFTVp0kQNGjTQM888o7/++kt+fn76/PPPr+qaJOnv0x433XSTevfurZ07d2rq1Kn6xz/+4XIxcV6KYj/z8vLSxIkTNWDAAHXu3Fl9+/Z13goeGhqqp5566qq2tyCaNWumHj166MMPP9QLL7ygqlWrXvG6Jk+erKlTpyoyMlLly5fXxx9/7DL/9ttvz/fW5sK65iY/Tz/9tL766ivFxMTo1KlTOWq89957nf8eNWqU5s6dm+tRouxT0HmdPsrtycLZgbdz585q06ZNjvmXWyeKkOdu1IIV5HUr+KVPTs2WmZlpJk6caMLCwozdbjdVqlQxrVu3NuPGjTNJSUnOfl999ZVp0aKF8fHxMaGhoWbixIlm1qxZOW4ZTUxMND179jSVKlUykpy36+Z1i2b2bcbHjx93ae/fv7+pUKFCjnpnzJhhWrdubXx9fU2lSpVM8+bNzfDhw83hw4edfUJCQnK99bZz5845bh/+4IMPTP369U3ZsmWL/WnF2bdE5/a6eEyzb1vN7dbp3Jw+fdqMGTPGNG/e3JQvX974+PiY66+/3owaNcr5VGRjjNm2bZuJiooyFStWNNWqVTODBg0yv/76a47b0t25Ffz77783Dz30kKlSpYqpWLGiueeee8zJkydzLJvXrdFXu59deit4toULF5qIiAhjt9tNQECAueeee8yff/7p0ievfa4gT2g2Jv//z2JjY40kM3bs2MuuJz/Zt6sXZL/xlOxHT+T1ulj29lxad1ZWlqldu7Zp1aqVW++d363gV7pOFA6bMR44WQ+gxOrdu7f279+vTZs2eboUALginJYC4GSMUWxsbI5D+wBQmnDkBgAAWAp/OBMAAFgK4QYAAFgK4QYAAFgK4QYAAFjKNXe3lMPh0OHDh1WpUqUie3Q/AAAoXMYYpaSkqFatWpf946XXXLg5fPjwZf/IHwAAKJkOHTqkOnXq5Nvnmgs3lSpVkvT34Pj5+Xm4GgAAUBDJyckKDg52fo/n55oLN9mnovz8/Ag3AACUMgW5pIQLigEAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgpJeuYFhY5cptCRy5SeecHT5QAAcM0i3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEvxaLhZu3atYmJiVKtWLdlsNi1duvSyy2RkZOi5555TSEiI7Ha7QkNDNWvWrKIvFgAAlArlPPnmaWlpCg8P1wMPPKB///vfBVqmd+/eOnr0qGbOnKmGDRvqyJEjcjgcRVwpAAAoLTwabrp3767u3bsXuP+KFSv0/fffa9++fQoICJAkhYaGFlF1AACgNCpV19x89dVXatOmjV577TXVrl1bjRo10jPPPKOzZ896ujQAAFBCePTIjbv27dundevWycfHR1988YVOnDihRx99VCdPntTs2bNzXSYjI0MZGRnO6eTk5OIqFwAAeECpOnLjcDhks9n0ySefqG3bturRo4cmTZqkuXPn5nn0Zvz48fL393e+goODi7lqAABQnEpVuKlZs6Zq164tf39/Z1vTpk1ljNGff/6Z6zKjRo1SUlKS83Xo0KHiKhcAAHhAqQo3HTp00OHDh5Wamups27Vrl8qUKaM6derkuozdbpefn5/LCwAAWJdHw01qaqri4+MVHx8vSUpISFB8fLwOHjwo6e+jLv369XP2v/vuu1W1alUNGDBA27Zt09q1a/Xss8/qgQcekK+vryc2AQAAlDAeDTebN29WRESEIiIiJEnDhg1TRESExowZI0k6cuSIM+hIUsWKFbVq1SqdOXNGbdq00T333KOYmBi98847HqkfAACUPB69W6pLly4yxuQ5f86cOTnamjRpolWrVhVhVQAAoDQrVdfcAAAAXA7hBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWIpHw83atWsVExOjWrVqyWazaenSpQVedv369SpXrpxatmxZZPUBAIDSx6PhJi0tTeHh4ZoyZYpby505c0b9+vXTTTfdVESVAQCA0qqcJ9+8e/fu6t69u9vLPfzww7r77rtVtmxZt472AAAA6yt119zMnj1b+/bt09ixYwvUPyMjQ8nJyS4vAABgXaUq3OzevVsjR47Uxx9/rHLlCnbQafz48fL393e+goODi7hKAADgSaUm3GRlZenuu+/WuHHj1KhRowIvN2rUKCUlJTlfhw4dKsIqAQCAp3n0mht3pKSkaPPmzYqLi9Njjz0mSXI4HDLGqFy5cvrmm2/0z3/+M8dydrtddru9uMsFAAAeUmrCjZ+fn7Zu3erSNnXqVH333Xf67LPPVK9ePQ9VBgAAShKPhpvU1FTt2bPHOZ2QkKD4+HgFBASobt26GjVqlP766y/NmzdPZcqU0fXXX++yfGBgoHx8fHK0AwCAa5dHw83mzZvVtWtX5/SwYcMkSf3799ecOXN05MgRHTx40FPlAQCAUshmjDGeLqI4JScny9/fX0lJSfLz8yu09aZnXlCzMSslSdteilZ571Jzxg8AgBLPne/vUnO3FAAAQEEQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQbgAAgKUQborA/hPpni4BAIBrFuGmEJxJz9SguVuc0z3e+UH9Zm5SUvp5D1YFAMC1iXBTCJ5YEK8Ne0+4tK3fc0KPL4jzUEUAAFy7CDdXad/xVK3dfVyOS9qzjNHa3ceVcCLNI3UBAHCtItxcpQOn8r++Zv9Jwg0AAMWJcHOVQgLK5zs/tGqFYqoEAABIhJurVr96RXW6rnqOgSxrs6nTddVVrxrhBgCA4kS4KQTv9o1QZINqLm0dGlbTu30jPFQRAADXLsJNIfAv76UP+rd2Ti9/oqPmDWwr//JeHqwKAIBrE+GmCIRWy/86HAAAUHQINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFI8Gm7Wrl2rmJgY1apVSzabTUuXLs23/5IlS9StWzdVr15dfn5+ioyM1MqVK4unWAAAUCp4NNykpaUpPDxcU6ZMKVD/tWvXqlu3blq+fLm2bNmirl27KiYmRnFxcUVcKQAAKC3KefLNu3fvru7duxe4/+TJk12mX331VX355Zf6+uuvFRERUcjVAQCA0sij4eZqORwOpaSkKCAgIM8+GRkZysjIcE4nJycXR2kAAMBDSvUFxW+88YZSU1PVu3fvPPuMHz9e/v7+zldwcHAxVggAAIpbqQ038+fP17hx47Ro0SIFBgbm2W/UqFFKSkpyvg4dOlSMVQIAgOJWKk9Lffrpp3rwwQe1ePFiRUVF5dvXbrfLbrcXU2UAAMDTSt2RmwULFmjAgAFasGCBevbs6elyAABACePRIzepqanas2ePczohIUHx8fEKCAhQ3bp1NWrUKP3111+aN2+epL9PRfXv319vv/222rVrp8TEREmSr6+v/P39PbINAACgZPHokZvNmzcrIiLCeRv3sGHDFBERoTFjxkiSjhw5ooMHDzr7z5gxQxcuXNCQIUNUs2ZN52vo0KEeqR8AAJQ8Hj1y06VLFxlj8pw/Z84cl+nY2NiiLQgAAJR6pe6aGwAAgPwQbgAAgKUQbgAAgKW4FW4WLVqkzMxM5/Sff/4ph8PhnE5PT9drr71WeNUBAAC4ya1w07dvX505c8Y53axZM+3fv985nZKSolGjRhVWbQBKgPTMCwoduUyhI5cpPfOCp8sBgMtyK9xcemdTfnc6AQAAeALX3AAAAEsh3AAAAEtx+yF+K1eudP6pA4fDodWrV+v333+XJJfrcQAAADzB7XDTv39/l+nBgwe7TNtstqurCAAA4Cq4FW4uvu0bAACgJHLrmpsHHnhAKSkpRVULAADAVXMr3MydO1dnz54tqloAAACu2lU95wYAACBbSXnop9sXFKekpMjHxyffPn5+fldcEAAAwNVwO9w0atQoz3nGGNlsNmVlZV1VUQAAAFfK7XDz2WefKSAgoChqAQAAuGpuh5sOHTooMDCwKGoBAAC4aoX+5xc4JQUAADzJrXATEhKismXL5jpv165dGj58uOrUqVMohQEAAFwJt8JNQkKCqlat6pxOT0/X7Nmz1bFjRzVr1kxr167VsGHDCr1IAACAgnL7mhtJ2rhxoz788EMtXrxYdevW1fbt27VmzRp17NixsOsDAABwi1tHbt58802FhYXpjjvuUJUqVbR27Vpt3bpVNpvN5YgOAACAp7h15GbEiBEaMWKEXnrppTyvvQEAAPAkt47cvPzyy1q8eLHq1aunESNG6Pfffy+qugAAAK6IW+Fm1KhR2rVrlz766CMlJiaqXbt2Cg8PlzFGp0+fLqoaAQAACuyKnnPTuXNnzZ07V0eOHNGjjz6qVq1aqVOnTmrfvr0mTZpU2DUCAAAU2FU9xM/Pz0+DBw/Wpk2b9Ouvv6pdu3aaMGFCYdUGAADgNrcuKD579qxWr16tf/3rX5L+Pk2VkZHxv5WVK6e9e/cWboUAAABucCvczJ07V8uWLXOGm/fee09hYWHy9fWVJO3cuVO1atXSU089VfiVAgAAFIBbp6U++eQTPfTQQy5t8+fP15o1a7RmzRq99tprWrRoUaEWCAAA4A63ws2ePXvUvHlz57SPj4/KlPnfKtq2batt27YVXnUAAABucuu01JkzZ1yusTl+/LjLfIfD4TIfAACguLl15KZOnTr5Prjvt99+46+CAwAAj3Ir3PTo0UNjxozRuXPncsw7e/asxo0bp549exZacQAAAO5y67TU6NGjtWjRIjVu3FiPPfaYGjVqJOnvu6Tee+89XbhwQaNHjy6SQgEAAArCrXATFBSkH3/8UY888ohGjhwpY4wkyWazqVu3bpo6daqCgoKKpFAAAICCcCvcSFK9evW0YsUKnTp1Snv27JEkNWzYUAEBAYVeHAAAgLvcDjfZAgIC1LZt28KsBQAA4Kpd1d+WAgAAKGkINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwBgcemZFxQ6cplCRy5TeuYFT5cDFDnCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBSPhpu1a9cqJiZGtWrVks1m09KlSy+7TGxsrFq1aiW73a6GDRtqzpw5RV4nAAAoPTwabtLS0hQeHq4pU6YUqH9CQoJ69uyprl27Kj4+Xk8++aQefPBBrVy5sogrBQAApcUV/+HMwtC9e3d17969wP2nT5+uevXq6c0335QkNW3aVOvWrdNbb72l6OjooioTAACUIqXqmpsNGzYoKirKpS06OlobNmzIc5mMjAwlJye7vAAAgHWVqnCTmJiooKAgl7agoCAlJyfr7NmzuS4zfvx4+fv7O1/BwcHFUSoAAPCQUhVursSoUaOUlJTkfB06dMjTJQEAgCLk0Wtu3FWjRg0dPXrUpe3o0aPy8/OTr69vrsvY7XbZ7fbiKA8AAJQAperITWRkpFavXu3StmrVKkVGRnqoIgAAUNJ4NNykpqYqPj5e8fHxkv6+1Ts+Pl4HDx6U9PcppX79+jn7P/zww9q3b5+GDx+uHTt2aOrUqVq0aJGeeuopT5QPAABKII+Gm82bNysiIkIRERGSpGHDhikiIkJjxoyRJB05csQZdCSpXr16WrZsmVatWqXw8HC9+eab+vDDD7kNHAAAOHn0mpsuXbrIGJPn/NyePtylSxfFxcUVYVUAAKA0K1XX3AAAAFwO4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFhKqfrDmSVZee9y2j+hp6fLAADgmseRGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAUOj2n0j32HsTbgAAwFU7k56pQXO3OKd7vPOD+s3cpKT088VeC+EGAABctScWxGvD3hMubev3nNDjC+KKvRbCDQAAuCr7jqdq7e7jclzSnmWM1u4+roQTacVaD+EGAABclQOn8r++Zv9Jwg0AAChFQgLK5zs/tGqFYqrkb4QbAABwVepXr6hO11XPESrK2mzqdF111atGuAEAAKXMu30jFNmgmktbh4bV9G7fiGKvhXADoMA8+dwKACWbf3kvfdC/tXN6+RMdNW9gW/mX9yr2Wgg3APJUkp5bAaB0Ca2W/3U4RYlwAyBPJem5FQBQUIQbALkqac+tAICCItwAyFVJe24FABQU4QZArkracysAoKAINwByVdKeWwEABUW4AZCnkvTcCgAoKMINgDyVpOdWAEBBEW4AFJgnn1sBAAVFuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZCuAEAAJZSIsLNlClTFBoaKh8fH7Vr106bNm3Kt//kyZPVuHFj+fr6Kjg4WE899ZTOnTtXTNUCAICSzOPhZuHChRo2bJjGjh2rX375ReHh4YqOjtaxY8dy7T9//nyNHDlSY8eO1fbt2zVz5kwtXLhQo0ePLubKAQBASeTxcDNp0iQNGjRIAwYMULNmzTR9+nSVL19es2bNyrX/jz/+qA4dOujuu+9WaGiobr75ZvXt2/eyR3sAAMC1waPhJjMzU1u2bFFUVJSzrUyZMoqKitKGDRtyXaZ9+/basmWLM8zs27dPy5cvV48ePXLtn5GRoeTkZJcXAFyr9p9I93QJQJHzaLg5ceKEsrKyFBQU5NIeFBSkxMTEXJe5++679dJLL+kf//iHvLy81KBBA3Xp0iXP01Ljx4+Xv7+/8xUcHFzo2wEAJdWZ9EwNmrvFOd3jnR/Ub+YmJaWf92BVQNHy+Gkpd8XGxurVV1/V1KlT9csvv2jJkiVatmyZXn755Vz7jxo1SklJSc7XoUOHirliAPCcJxbEa8PeEy5t6/ec0OML4jxUEVD0ynnyzatVq6ayZcvq6NGjLu1Hjx5VjRo1cl3mhRde0H333acHH3xQktS8eXOlpaXpoYce0nPPPacyZVzzmt1ul91uL5oNAIASbN/xVK3dfTxHe5YxWrv7uBJOpKletQoeqAwoWh49cuPt7a3WrVtr9erVzjaHw6HVq1crMjIy12XS09NzBJiyZctKkowxRVcsAJQyB07lf33N/pNpxVQJULw8euRGkoYNG6b+/furTZs2atu2rSZPnqy0tDQNGDBAktSvXz/Vrl1b48ePlyTFxMRo0qRJioiIULt27bRnzx698MILiomJcYYcAIAUElA+3/mhVTlqA2vyeLjp06ePjh8/rjFjxigxMVEtW7bUihUrnBcZHzx40OVIzfPPPy+bzabnn39ef/31l6pXr66YmBi98sorntoEACiR6levqE7XVde63cfluKi9rM2mDg2rcUoKlmUz19i5nOTkZPn7+yspKUl+fn6eLgco8dIzL6jZmJWSpG0vRau8t8d/J4IbktLP69FPftH6iy4q7nRddb3bN0L+5b08WBmsqCg/L9z5/i51d0sBAArOv7yXPujf2jm9/ImOmjewLcEGlka4AYBrSGi1/K/DAayAcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACylnKcLAAAA1lDeu5z2T+jp6TI4cgMAAKyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACylRISbKVOmKDQ0VD4+PmrXrp02bdqUb/8zZ85oyJAhqlmzpux2uxo1aqTly5cXU7UAAKAkK+fpAhYuXKhhw4Zp+vTpateunSZPnqzo6Gjt3LlTgYGBOfpnZmaqW7duCgwM1GeffabatWvrwIEDqly5cvEXDwAAShyPh5tJkyZp0KBBGjBggCRp+vTpWrZsmWbNmqWRI0fm6D9r1iydOnVKP/74o7y8vCRJoaGhxVkyAAAowTx6WiozM1NbtmxRVFSUs61MmTKKiorShg0bcl3mq6++UmRkpIYMGaKgoCBdf/31evXVV5WVlVVcZQMAgBLMo0duTpw4oaysLAUFBbm0BwUFaceOHbkus2/fPn333Xe65557tHz5cu3Zs0ePPvqozp8/r7Fjx+bon5GRoYyMDOd0cnJy4W4EAAAoUUrEBcXucDgcCgwM1IwZM9S6dWv16dNHzz33nKZPn55r//Hjx8vf39/5Cg4OLuaKAQBAcfJouKlWrZrKli2ro0ePurQfPXpUNWrUyHWZmjVrqlGjRipbtqyzrWnTpkpMTFRmZmaO/qNGjVJSUpLzdejQocLdCAAAUKJ4NNx4e3urdevWWr16tbPN4XBo9erVioyMzHWZDh06aM+ePXI4HM62Xbt2qWbNmvL29s7R3263y8/Pz+UFAACsy+OnpYYNG6YPPvhAc+fO1fbt2/XII48oLS3NefdUv379NGrUKGf/Rx55RKdOndLQoUO1a9cuLVu2TK+++qqGDBniqU0ALK28dzntn9BT+yf0VHlvj99gCQCX5fFPqj59+uj48eMaM2aMEhMT1bJlS61YscJ5kfHBgwdVpsz/MlhwcLBWrlypp556Si1atFDt2rU1dOhQjRgxwlObAAAAShCbMcZ4uojilJycLH9/fyUlJXGKCsA1IT3zgpqNWSlJ2vZSNEfgUCq58/3t8dNSAAAAhYlwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIUnOQGAxWX/CQ3gWsGRGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCnlPF1AcTPGSJKSk5M9XAkAACio7O/t7O/x/Fxz4SYlJUWSFBwc7OFKAACAu1JSUuTv759vH5spSASyEIfDocOHD6tSpUqy2WxKTk5WcHCwDh06JD8/P0+X53GMhyvGwxXj4YrxcMV4uGI8XF3teBhjlJKSolq1aqlMmfyvqrnmjtyUKVNGderUydHu5+fHzncRxsMV4+GK8XDFeLhiPFwxHq6uZjwud8QmGxcUAwAASyHcAAAAS7nmw43dbtfYsWNlt9s9XUqJwHi4YjxcMR6uGA9XjIcrxsNVcY7HNXdBMQAAsLZr/sgNAACwFsINAACwFMINAACwFMINAACwFEuGmylTpig0NFQ+Pj5q166dNm3alG//xYsXq0mTJvLx8VHz5s21fPlyl/nGGI0ZM0Y1a9aUr6+voqKitHv37qLchELlznh88MEH6tixo6pUqaIqVaooKioqR//7779fNpvN5XXLLbcU9WYUGnfGY86cOTm21cfHx6XPtbR/dOnSJcd42Gw29ezZ09mntO4fa9euVUxMjGrVqiWbzaalS5dedpnY2Fi1atVKdrtdDRs21Jw5c3L0cffzqKRwdzyWLFmibt26qXr16vLz81NkZKRWrlzp0ufFF1/MsW80adKkCLei8Lg7HrGxsbn+v5KYmOjS71rZP3L7XLDZbAoLC3P2Kcz9w3LhZuHChRo2bJjGjh2rX375ReHh4YqOjtaxY8dy7f/jjz+qb9++GjhwoOLi4nTbbbfptttu0++//+7s89prr+mdd97R9OnT9dNPP6lChQqKjo7WuXPnimuzrpi74xEbG6u+fftqzZo12rBhg4KDg3XzzTfrr7/+cul3yy236MiRI87XggULimNzrpq74yH9/TTNi7f1wIEDLvOvpf1jyZIlLmPx+++/q2zZsrrzzjtd+pXG/SMtLU3h4eGaMmVKgfonJCSoZ8+e6tq1q+Lj4/Xkk0/qwQcfdPlCv5L9raRwdzzWrl2rbt26afny5dqyZYu6du2qmJgYxcXFufQLCwtz2TfWrVtXFOUXOnfHI9vOnTtdtjcwMNA571raP95++22XcTh06JACAgJyfHYU2v5hLKZt27ZmyJAhzumsrCxTq1YtM378+Fz79+7d2/Ts2dOlrV27dmbw4MHGGGMcDoepUaOGef31153zz5w5Y+x2u1mwYEERbEHhcnc8LnXhwgVTqVIlM3fuXGdb//79za233lrYpRYLd8dj9uzZxt/fP8/1Xev7x1tvvWUqVapkUlNTnW2lef/IJsl88cUX+fYZPny4CQsLc2nr06ePiY6Odk5f7fiWFAUZj9w0a9bMjBs3zjk9duxYEx4eXniFeUhBxmPNmjVGkjl9+nSefa7l/eOLL74wNpvN7N+/39lWmPuHpY7cZGZmasuWLYqKinK2lSlTRlFRUdqwYUOuy2zYsMGlvyRFR0c7+yckJCgxMdGlj7+/v9q1a5fnOkuKKxmPS6Wnp+v8+fMKCAhwaY+NjVVgYKAaN26sRx55RCdPnizU2ovClY5HamqqQkJCFBwcrFtvvVV//PGHc961vn/MnDlTd911lypUqODSXhr3D3dd7rOjMMa3NHM4HEpJScnx2bF7927VqlVL9evX1z333KODBw96qMLi0bJlS9WsWVPdunXT+vXrne3X+v4xc+ZMRUVFKSQkxKW9sPYPS4WbEydOKCsrS0FBQS7tQUFBOc5zZktMTMy3f/Z/3VlnSXEl43GpESNGqFatWi7/A95yyy2aN2+eVq9erYkTJ+r7779X9+7dlZWVVaj1F7YrGY/GjRtr1qxZ+vLLL/Xxxx/L4XCoffv2+vPPPyVd2/vHpk2b9Pvvv+vBBx90aS+t+4e78vrsSE5O1tmzZwvl/7/S7I033lBqaqp69+7tbGvXrp3mzJmjFStWaNq0aUpISFDHjh2VkpLiwUqLRs2aNTV9+nR9/vnn+vzzzxUcHKwuXbrol19+kVQ4n8+l1eHDh/Xf//43x2dHYe4f19xfBUfBTZgwQZ9++qliY2NdLqK96667nP9u3ry5WrRooQYNGig2NlY33XSTJ0otMpGRkYqMjHROt2/fXk2bNtX777+vl19+2YOVed7MmTPVvHlztW3b1qX9Wto/kLv58+dr3Lhx+vLLL12uMenevbvz3y1atFC7du0UEhKiRYsWaeDAgZ4otcg0btxYjRs3dk63b99ee/fu1VtvvaWPPvrIg5V53ty5c1W5cmXddtttLu2FuX9Y6shNtWrVVLZsWR09etSl/ejRo6pRo0auy9SoUSPf/tn/dWedJcWVjEe2N954QxMmTNA333yjFi1a5Nu3fv36qlatmvbs2XPVNRelqxmPbF5eXoqIiHBu67W6f6SlpenTTz8t0AdOadk/3JXXZ4efn598fX0LZX8rjT799FM9+OCDWrRoUY7TdpeqXLmyGjVqZLl9Iy9t27Z1buu1un8YYzRr1izdd9998vb2zrfv1ewflgo33t7eat26tVavXu1sczgcWr16tctv3xeLjIx06S9Jq1atcvavV6+eatSo4dInOTlZP/30U57rLCmuZDykv+/+efnll7VixQq1adPmsu/z559/6uTJk6pZs2ah1F1UrnQ8LpaVlaWtW7c6t/Va3D+kvx+fkJGRoXvvvfey71Na9g93Xe6zozD2t9JmwYIFGjBggBYsWODyeIC8pKamau/evZbbN/ISHx/v3NZrcf+QpO+//1579uwp0C9GV7V/FMplySXIp59+aux2u5kzZ47Ztm2beeihh0zlypVNYmKiMcaY++67z4wcOdLZf/369aZcuXLmjTfeMNu3bzdjx441Xl5eZuvWrc4+EyZMMJUrVzZffvml+e2338ytt95q6tWrZ86ePVvs2+cud8djwoQJxtvb23z22WfmyJEjzldKSooxxpiUlBTzzDPPmA0bNpiEhATz7bffmlatWpnrrrvOnDt3ziPb6A53x2PcuHFm5cqVZu/evWbLli3mrrvuMj4+PuaPP/5w9rmW9o9s//jHP0yfPn1ytJfm/SMlJcXExcWZuLg4I8lMmjTJxMXFmQMHDhhjjBk5cqS57777nP337dtnypcvb5599lmzfft2M2XKFFO2bFmzYsUKZ5/LjW9J5u54fPLJJ6ZcuXJmypQpLp8dZ86ccfZ5+umnTWxsrElISDDr1683UVFRplq1aubYsWPFvn3ucnc83nrrLbN06VKze/dus3XrVjN06FBTpkwZ8+233zr7XEv7R7Z7773XtGvXLtd1Fub+YblwY4wx7777rqlbt67x9vY2bdu2NRs3bnTO69y5s+nfv79L/0WLFplGjRoZb29vExYWZpYtW+Yy3+FwmBdeeMEEBQUZu91ubrrpJrNz587i2JRC4c54hISEGEk5XmPHjjXGGJOenm5uvvlmU716dePl5WVCQkLMoEGDSsX/jNncGY8nn3zS2TcoKMj06NHD/PLLLy7ru5b2D2OM2bFjh5FkvvnmmxzrKs37R/atu5e+sre/f//+pnPnzjmWadmypfH29jb169c3s2fPzrHe/Ma3JHN3PDp37pxvf2P+vlW+Zs2axtvb29SuXdv06dPH7Nmzp3g37Aq5Ox4TJ040DRo0MD4+PiYgIMB06dLFfPfddznWe63sH8b8/ZgMX19fM2PGjFzXWZj7h80YY9w/3gMAAFAyWeqaGwAAAMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINYHGJiYkaOnSoGjZsKB8fHwUFBalDhw6aNm2a0tPTPV1enkJDQ2Wz2fTpp5/mmBcWFiabzaY5c+YUf2G5SExM1OOPP6769evLbrcrODhYMTExOf72lCSNHz9eZcuW1euvv+5sy97WvF7333+/JOU5P7cxAq5l5TxdAICis2/fPnXo0EGVK1fWq6++qubNm8tut2vr1q2aMWOGateurV69euW67Pnz5+Xl5VXMFbsKDg7W7NmzdddddznbNm7cqMTERFWoUMGDlf3P/v37nWP8+uuvq3nz5jp//rxWrlypIUOGaMeOHS79Z82apeHDh2vWrFl69tlnJUk///yzsrKyJEk//vij/vOf/2jnzp3y8/OTJPn6+jqXnz17tm655RaXdVauXLkItxAoha7ojzYAKBWio6NNnTp1TGpqaq7zHQ6H89+SzNSpU01MTIwpX7688++JTZ061dSvX994eXmZRo0amXnz5jmXSUhIMJJMXFycs+306dNGklmzZo0x5n9/g+b//u//TPPmzY3dbjft2rVz+eO0uQkJCTEjR440drvdHDx40Nk+aNAg8/jjjxt/f3+Xv+V0+vRpM3DgQFOtWjVTqVIl07VrVxMfH++cv2fPHtOrVy8TGBhoKlSoYNq0aWNWrVqV4z1feeUVM2DAAFOxYkUTHBxs3n///Xzr7N69u6ldu3auY3z69GmX6djYWFO7dm2TmZlpatWqZdavX59jmezxunRZY/7+GX3xxRf51gPAGE5LARZ18uRJffPNNxoyZEieRzlsNpvL9Isvvqjbb79dW7du1QMPPKAvvvhCQ4cO1dNPP63ff/9dgwcP1oABA7RmzRq363n22Wf15ptv6ueff1b16tUVExOj8+fP57tMUFCQoqOjNXfuXElSenq6Fi5cqAceeCBH3zvvvFPHjh3Tf//7X23ZskWtWrXSTTfdpFOnTkmSUlNT1aNHD61evVpxcXG65ZZbFBMTo4MHD7qs580331SbNm0UFxenRx99VI888oh27tyZa32nTp3SihUr8hzjS4+ozJw5U3379pWXl5f69u2rmTNn5rv9AK6Qp9MVgKKxceNGI8ksWbLEpb1q1aqmQoUKpkKFCmb48OHOdknmySefdOnbvn17M2jQIJe2O++80/To0cMY496Rm08//dTZ5+TJk8bX19csXLgwz/pDQkLMW2+9ZZYuXWoaNGhgHA6HmTt3romIiDDGGJcjNz/88IPx8/Mz586dc1lHgwYN8j3yEhYWZt59912X97z33nud0w6HwwQGBppp06bluvxPP/2U6xjnJikpyfj6+jqPJsXFxZmKFSualJQUl36XO3Lj4+Pj/Pllvw4cOHDZ9weuJRy5Aa4xmzZtUnx8vMLCwpSRkeEyr02bNi7T27dvV4cOHVzaOnTooO3bt7v9vpGRkc5/BwQEqHHjxgVaT8+ePZWamqq1a9dq1qxZuR61+fXXX5WamqqqVauqYsWKzldCQoL27t0r6e8jN88884yaNm2qypUrq2LFitq+fXuOIzctWrRw/ttms6lGjRo6duxYrrUZYwq07ZK0YMECNWjQQOHh4ZKkli1bKiQkRAsXLizwOiTprbfeUnx8vMurVq1abq0DsDouKAYsqmHDhrLZbDlOqdSvX1+S60Wq2dy9SLdMmb9/P7r4S/5yp5rcVa5cOd13330aO3asfvrpJ33xxRc5+qSmpqpmzZqKjY3NMS/71NAzzzyjVatW6Y033lDDhg3l6+urO+64Q5mZmS79L72I2mazyeFw5FrbddddJ5vNluOi4dzMnDlTf/zxh8qV+9/HrsPh0KxZszRw4MDLLp+tRo0aatiwYYH7A9cijtwAFlW1alV169ZN7733ntLS0q5oHU2bNtX69etd2tavX69mzZpJkqpXry5JOnLkiHN+fHx8ruvauHGj89+nT5/Wrl271LRp0wLV8cADD+j777/XrbfeqipVquSY36pVKyUmJqpcuXJq2LChy6tatWrOuu+//37dfvvtat68uWrUqKH9+/cX6P3zEhAQoOjoaE2ZMiXXMT5z5owkaevWrdq8ebNiY2NdjrjExsZqw4YNBQpHAAqOIzeAhU2dOlUdOnRQmzZt9OKLL6pFixYqU6aMfv75Z+3YsUOtW7fOd/lnn31WvXv3VkREhKKiovT1119ryZIl+vbbbyX9ffTnxhtv1IQJE1SvXj0dO3ZMzz//fK7reumll1S1alUFBQXpueeeU7Vq1XTbbbcVaDuaNm2qEydOqHz58rnOj4qKUmRkpG677Ta99tpratSokQ4fPqxly5bp9ttvV5s2bXTddddpyZIliomJkc1m0wsvvJDnERl3TJkyRR06dFDbtm310ksvqUWLFrpw4YJWrVqladOmafv27Zo5c6batm2rTp065Vj+hhtu0MyZM12ee5OfM2fOKDEx0aWtUqVKJebWeKBE8PRFPwCK1uHDh81jjz1m6tWrZ7y8vEzFihVN27Ztzeuvv27S0tKc/ZTHbcb53QpujDHbtm0zkZGRxtfX17Rs2dJ88803uV5Q/PXXX5uwsDDj7e1t2rZta3799dd8686+oDgvl94KnpycbB5//HFTq1Yt4+XlZYKDg80999zjvI08ISHBdO3a1fj6+prg4GDz3nvvmc6dO5uhQ4fm+57h4eHO2+LzcvjwYTNkyBATEhJivL29Te3atU2vXr3MmjVrTEZGhqlatap57bXXcl124sSJJjAw0GRmZhpjLn9BcW6v8ePH51sfcK2xGePGFXEA4KbY2Fh17dpVp0+f5mFzAIoF19wAAABLIdwAAABL4bQUAACwFI7cAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAAS/n/L+zTgW40TfcAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "res_t.cal.plot_cal(1)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAHHCAYAAACRAnNyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA2j0lEQVR4nO3de5yN9f7//+fCnE+MMSTjXGEYh3HY02wjIZvZdNhSUkludBhiq0h9cmjvGmWXyikV0kFJRbUrNr4zyCaHKOVclDDO5sgMs96/P/rN2pYZY4bhemse99tt3W5zva/Ta611zaznel/v6xqXMcYIAADAQhWcLgAAAOBcCCoAAMBaBBUAAGAtggoAALAWQQUAAFiLoAIAAKxFUAEAANYiqAAAAGsRVAAAgLUIKgAuubFjx8rlcnm11a1bV/fdd59n+q233pLL5dK6desuc3XeUlNT5XK5lJqa6mgdAH5HUIFcLleJHpfrD/fUqVP11ltvXZZ9XazNmzdr7Nix2r1792Xf99q1azV48GBFR0crKChItWvXVu/evbV9+/aL3nZGRobGjRun5s2bKzg4WAEBAWratKlGjhypffv2lUH1zrPxOLvhhhu8fucCAgIUExOjl19+WW63u9h1f/31V9WsWVMul0v/+te/zrnc1q1bNWLECLVo0UIhISG66qqrlJiY6FhAdLvdeuGFF1SvXj35+/srJiZG77//fqm2sWTJEt14440KCwtTSEiIYmNjNXfuXM/8I0eOaMKECUpISFC1atVUuXJl/elPf/Ja5mzffvutevbsqfDwcAUGBqpp06Z69dVXL/h54sJVcroAOO+dd97xmn777be1ePHiQu2NGze+LPVMnTpVERERXt+2bbV582aNGzdON9xwg+rWrXtZ9/38889r5cqVuv322xUTE6O0tDRNnjxZrVq10urVq9W0adML2u7PP/+szp0769dff9Xtt9+uQYMGydfXV99//71mzJih+fPnl0kY2rZtmypUcO670rmOs4SEBJ04cUK+vr6O1FWrVi0lJydLkg4fPqw5c+bo73//uw4dOqRnn322yHWOHTumbt26KTs7WwkJCRoxYoSioqJ0xx13FFr2zTff1IwZM/S3v/1NDz/8sNLT0zV9+nT96U9/0sKFC9W5c+dL+vzO9tRTT2n8+PEaOHCg2rRpo08//VR33XWXXC6X7rzzzvOuP2vWLA0YMEBdunTRc889p4oVK2rbtm3as2ePZ5lVq1bpqaeeUvfu3fV///d/qlSpkj7++GPdeeednt/hM/3nP/9Rjx491LJlSz399NMKDg7WTz/9pN9++63Mnz9KwABnSUpKMiU5NLKzsy/J/qOjo02HDh0uybbL2rx584wkk5KSctn3vXLlSpObm+vVtn37duPn52f69u17Qds8deqUad68uQkMDDQrVqwoND89Pd08+eSTpd7umDFjzntMzZo1y0gya9euLfX23W63ycnJKdU6Nh5nHTp0MNHR0V5tJ06cMHXq1DEhISHm9OnThdY5efKkSUhIMKGhoWbVqlXm5MmTJjEx0fj5+Zlly5YVWn7dunUmMzPTq+3w4cOmWrVqJj4+vmyf0Hn89ttvxsfHxyQlJXna3G63ad++valVq1aRz/dMu3btMgEBAeaRRx4pdrmff/7Z7N6926vN7XabG2+80fj5+ZmsrCxPe3p6uqlevbq59dZbTX5+/gU8K5Q1ggoKKSqoFPwBXbdunWnfvr0JCAgwQ4cONcb8/ody9OjRpkGDBsbX19fUqlXLPP744+bkyZNe25g5c6bp2LGjqVatmvH19TWNGzc2U6dO9VqmTp06RpLXo+DDpOCDbMWKFWbIkCEmIiLChIWFmUGDBpnc3Fxz7Ngxc88995jKlSubypUrm8cff9y43W6v7efn55uJEyeaJk2aGD8/PxMZGWkGDRpkjh49WqiOxMREs2LFCtOmTRvj5+dn6tWrZ2bPnu1ZpqCesx9OhJYztWrVyrRq1eqC1v3ggw+MJPPss8+WaPnly5ebXr16maioKM97P2zYsEKhoaigUqdOHdOvXz/PdMHruWzZMjNo0CATHh5uQkJCzD333HPO92fhwoUmNjbW+Pn5mYkTJxpjLv44S0lJKfJ9/PDDD02rVq2Mv7+/qVq1qunbt6/57bffvJbp16+fCQoKMr/99pu5+eabTVBQkImIiDCPPvroeT90jSk6qBhjTK9evYwks2/fPq92t9tt7rjjDhMWFma++eYbT3tubq7561//aqpUqWJ+/PHH8+7XGGNuu+02Ex4eXqJly8qUKVOMpEI1zpkzx/O7XpyRI0caX19fc/z4cWOMMZmZmYV+54vz6quvGknm+++/97RNmzbNSDKbN282xhiTlZVFYHEYp35QYkeOHFG3bt1055136u6771b16tXldrvVs2dPff311xo0aJAaN26sTZs2aeLEidq+fbsWLFjgWX/atGmKjo5Wz549ValSJX3++ed6+OGH5Xa7lZSUJEl6+eWXNWTIEAUHB+upp56SJFWvXt2rjiFDhqhGjRoaN26cVq9erddff12VK1fWf//7X9WuXVvPPfecvvzyS02YMEFNmzbVvffe61n3gQce0FtvvaX+/fvrkUce0a5duzR58mRt2LBBK1eulI+Pj2fZnTt3qlevXhowYID69eunmTNn6r777lNsbKyio6OVkJCgRx55RK+++qqefPJJz6mx4k6R5ebmKjMzs0Svd0RERImWO5MxRgcOHFB0dHSp15Wkzz77TJJ0zz33lGj5efPmKScnRw899JCqVq2qNWvWaNKkSfrtt980b968C6ph8ODBqly5ssaOHatt27Zp2rRp+uWXXzyDXAts27ZNffr00QMPPKCBAwfquuuuk1R2x9mZCo6ZNm3aKDk5WQcOHNArr7yilStXasOGDapcubJn2fz8fHXt2lXt2rXTv/71Ly1ZskQvvviiGjRooIceeuiCXpPdu3fL5XJ57UeSRowYoUWLFmnx4sVq06aNp93X11cff/yxevXqpW7dumn16tW66qqrit1HWlpaiY65U6dOKT09vUR1h4eHF3t6b8OGDQoKCir0O9O2bVvP/D//+c/nXH/JkiVq1KiRvvzySz3++OPau3evqlSpoqSkJI0bN+68pxbT0tIkef+uLVmyRKGhodq7d69uueUWbd++XUFBQbrnnns0ceJE+fv7n/d5o4w5nZRgn3P1qEgyr732mlf7O++8YypUqFDom89rr71mJJmVK1d62orqmu/ataupX7++V9u5uuQLvnF37drV61tTXFyccblc5sEHH/S0nT592tSqVctrOytWrDCSzHvvvee13YULFxZqL/jGvXz5ck/bwYMHjZ+fn3n00Uc9baU99XOuXpiiHhfinXfeMZLMjBkzLmj9li1bmrCwsBIvX9R7mpycbFwul/nll188baXpUYmNjTV5eXme9hdeeMFIMp9++qnXupLMwoULS1RTaY6zs3tU8vLyTGRkpGnatKk5ceKEZ7l///vfRpIZPXq0p61fv35GknnmmWe8ttmyZUsTGxtbaF9n69Chg2nUqJE5dOiQOXTokNm6dat5/PHHjSSTmJh43vUv1PLly43L5TJPP/30eZcteH1K8ti1a1ex20pMTCz0vhjz+2llSeaJJ54odv3Q0FBTpUoV4+fnZ55++mnz0UcfmbvuuqtE6x45csRERkaa9u3be7XHxMSYwMBAExgYaIYMGWI+/vhjM2TIECPJ3HnnncVuE5cGPSooMT8/P/Xv39+rbd68eWrcuLEaNWqkw4cPe9pvvPFGSVJKSoquv/56SVJAQIBnfnp6uk6dOqUOHTpo0aJFSk9PV1hYWInqGDBggNc363bt2mnVqlUaMGCAp61ixYpq3bq11q9f71VrWFiYunTp4lVrbGysgoODlZKSorvuusvT3qRJE7Vv394zXa1aNV133XX6+eefS1RnUbp27arFixdf8PrF2bp1q5KSkhQXF6d+/fpd0DYyMjIUEhJS4uXPfE+zs7N14sQJXX/99TLGaMOGDapdu3apaxg0aJBXz9ZDDz2kJ598Ul9++aV69uzpaa9Xr566du1abE0Xc5wVWLdunQ4ePKixY8d6fZtOTExUo0aN9MUXXxQajPnggw96Tbdv377Q4PRz2bp1q6pVq+bV1rNnT82YMaNUdZfUwYMHddddd6levXoaMWLEeZdv3rx5iY/hGjVqFDv/xIkT8vPzK9Re8DqfOHGi2PWzsrLkdrs1fvx4jRw5UpL0t7/9TUePHtUrr7yiJ598ssjj2e12q2/fvjp+/LgmTZpUaJs5OTl68MEHPVf53HbbbcrLy9P06dP1zDPP6Jprrim2LpQtggpK7Oqrry50JcSOHTu0ZcuWQn9YCxw8eNDz88qVKzVmzBitWrVKOTk5XsuV5gPk7A+/gvWioqIKtR87dsyr1vT0dEVGRp631qL2I0lVqlTx2mZpXXXVVeftgr8QaWlpSkxMVFhYmD766CNVrFjxgrYTGhpaqiD266+/avTo0frss88KvS4lPT1wtrM/BIKDg3XVVVcVugS8Xr16Ra5fVsdZgV9++UWSPKeWztSoUSN9/fXXXm3+/v6Ffh9Kc9zUrVtXb7zxhtxut3766Sc9++yzOnTo0CU55ZCdna2//vWvyszM1Ndff63g4ODzrlOlSpUyuzIoICBAubm5hdpPnjzpmX++9bOzs9WnTx+v9j59+mjhwoXasGGDEhISCq03ZMgQLVy4UG+//baaN29eaJsF2zjTXXfdpenTp2vVqlUElcuMoIISK+qPhtvtVrNmzfTSSy8VuU5BePjpp5/UqVMnNWrUSC+99JKioqLk6+urL7/8UhMnTjzvPSLOdK4P4aLajTFetUZGRuq9994rcv2zP1zOtZ8zt1laJ06cKPEH+Pm+jRZIT09Xt27ddPz4ca1YsUI1a9a84PoaNWqkDRs2aM+ePYWC39ny8/PVpUsXHT16VCNHjlSjRo0UFBSkvXv36r777ivVe3ohijoey/I4u1AXGhILBAUFeQWB+Ph4tWrVSk8++WSZ3scjLy9Pt912m77//nstWrSoxJez5+Xl6ejRoyVatlq1asW+HldddZVSUlJkjPHqJd2/f78knfdYrlmzpnbs2FFofFHBl5GiwuG4ceM0depUjR8/vsixWDVr1tSPP/5Yqm3i0iKo4KI0aNBA3333nTp16lTozqNn+vzzz5Wbm6vPPvvMq6ciJSWl0LLFbedia12yZIni4+PP+02tpEpb69y5cwudPjuXkgSikydPqkePHtq+fbuWLFmiJk2alKqes/Xo0UPvv/++3n33XY0aNarYZTdt2qTt27dr9uzZXgOWL/bU1o4dO9SxY0fPdFZWlvbv36/u3bufd91LcZzVqVNH0u+DdwtOaRbYtm2bZ/6lEhMTo7vvvlvTp0/XY489dkGn087mdrt17733aunSpfrwww/VoUOHEq/73//+1+v9Kc6uXbuKvb9QixYt9Oabb2rLli1ex+4333zjmV+c2NhY7dixQ3v37lX9+vU97QU3JTz7y8eUKVM0duxYDRs2zHOqqKhtLl68WHv37vXqRTvXNnHpcWdaXJTevXtr7969euONNwrNO3HihLKzsyX971vmmR++6enpmjVrVqH1goKCdPz48UtSa35+vv7xj38Umnf69OkL2mdQUJAklXjdgjEqJXmcT35+vu644w6tWrVK8+bNU1xcXKnrP1uvXr3UrFkzPfvss1q1alWh+ZmZmZ6rZIp6T40xeuWVVy6qhtdff12nTp3yTE+bNk2nT59Wt27dzrvupTjOWrdurcjISL322mtepym++uorbdmyRYmJiefdxsUaMWKETp06dc6ey9IaMmSI5s6dq6lTp+q2224r1boFY1RK8jhfr+DNN98sHx8fTZ061dNmjNFrr72mq6++2jO+Tfq9l2Xr1q1ex0bBDe3OHL/jdrs1a9YshYeHKzY21tM+d+5cPfLII+rbt2+xr2Pv3r0LbVP6/UZ5lSpV0g033FDsc0LZo0cFF+Wee+7Rhx9+qAcffFApKSmKj49Xfn6+tm7dqg8//FCLFi1S69atddNNN8nX11c9evTQAw88oKysLL3xxhuKjIz0dPMWiI2N1bRp0/TPf/5TDRs2VGRkZKFvsheiQ4cOeuCBB5ScnKyNGzfqpptuko+Pj3bs2KF58+bplVdeUa9evUq1zRYtWqhixYp6/vnnlZ6eLj8/P914443nHAdTlmNUHn30UX322Wfq0aOHjh49qnfffddr/t133+35ueDy2lmzZhV7x18fHx998skn6ty5sxISEtS7d2/Fx8fLx8dHP/74o+bMmaMqVaro2WefVaNGjdSgQQM99thj2rt3r0JDQ/Xxxx9fdNd4Xl6eOnXqpN69e2vbtm2aOnWq/vznP3sNpD2XS3Gc+fj46Pnnn1f//v3VoUMH9enTx3N5ct26dfX3v//9op5vSTRp0kTdu3fXm2++qaefflpVq1a94G29/PLLmjp1quLi4hQYGFjouLn11ls9AbwoZTlGpVatWho2bJgmTJigU6dOqU2bNlqwYIFWrFih9957z+u00ahRozR79myvXpqbb75ZnTp1UnJysg4fPqzmzZtrwYIF+vrrrzV9+nTPQN01a9bo3nvvVdWqVdWpU6dCp3+vv/56T49My5Ytdf/992vmzJk6ffq0OnTooNTUVM2bN0+jRo26qFOruEBOXW4EexV3w7ei5OXlmeeff95ER0cbPz8/U6VKFRMbG2vGjRtn0tPTPct99tlnJiYmxvj7+5u6deua559/3sycObPQZYxpaWkmMTHRhISEFHnDt7PvXFpw6euhQ4e82gtuvnW2119/3cTGxpqAgAATEhJimjVrZkaMGOF1M62CG4qdrUOHDoUuaX3jjTdM/fr1TcWKFS/rDd8KLhk/1+NMkyZNOuflvEU5duyYGT16tGnWrJkJDAw0/v7+pmnTpmbUqFFm//79nuU2b95sOnfubIKDg01ERIQZOHCg+e6774wkM2vWLM9yF3LDtypVqpjg4GDTt29fc+TIkULrnuty3Ys9zs51w7e5c+eali1bGj8/PxMeHl7sDd/OVpI78xpT/O9ZamqqkWTGjBlz3u0Up+AS6nM9zndJcVnLz883zz33nKlTp47x9fU10dHR5t133z1n3WfXl5mZaYYOHWpq1KhhfH19TbNmzQqtf77bApx5rBrz+9+0sWPHmjp16hgfHx/TsGFDzw0Fcfm5jLmIkYEArgi9e/fW7t27tWbNGqdLAYBS4dQP8AdnjFFqamqhLn4AuBLQowIAAKzFVT8AAMBaBBUAAGAtggoAALAWQQUAAFjrir7qx+12a9++fQoJCblkt10HAABlyxijzMxM1axZUxUqFN9nckUHlX379p33H6cBAAA77dmzR7Vq1Sp2mSs6qISEhEj6/YmGhoY6XA0AACiJjIwMRUVFeT7Hi3NFB5WC0z2hoaEEFQAArjAlGbbBYFoAAGAtggoAALAWQQUAAFjLmqAyfvx4uVwuDRs2zOlSAACAJawIKmvXrtX06dMVExPjdCkAAMAijgeVrKws9e3bV2+88YaqVKnidDkAAMAijgeVpKQkJSYmqnPnzuddNjc3VxkZGV4PAADwx+XofVQ++OADffvtt1q7dm2Jlk9OTta4ceMucVUAAMAWjvWo7NmzR0OHDtV7770nf3//Eq0zatQopaenex579uy5xFUCAAAnuYwxxokdL1iwQLfeeqsqVqzoacvPz5fL5VKFChWUm5vrNa8oGRkZCgsLU3p6OnemBQDgClGaz2/HTv106tRJmzZt8mrr37+/GjVqpJEjR543pAAAgD8+x4JKSEiImjZt6tUWFBSkqlWrFmoHAADlk+NX/QAAAJyLVf89OTU11ekSAACARehRAQAAheTknVbdJ75Q3Se+UE7eacfqIKgAAABrEVQAAIC1CCoAAMBaBBUAAGAtggoAALAWQQUAAFiLoAIAAKxFUAEAANYiqAAAAGsRVAAAgLUIKgAAwFoEFQAAYC2CCgAAsBZBBQAAFGv34RzH9k1QAQAAXo7n5Gng7PWe6e6vrtC9M9YoPefUZa+FoAIAALw88v5GrfrpsFfbyp2HNeT9DZe9FoIKAADw+PlQlpbvOCT3We35xmj5jkPadTj7stZDUAEAAB6/HC1+PMruIwQVAADgkDrhgcXOr1s16DJV8juCCgAA8KhfLVgJ11QrFBAqulxKuKaa6kUQVAAAgIMm9WmpuAYRXm3xDSM0qU/Ly14LQQUAAHgJC/TRG/1iPdNfPtJebw9oq7BAn8teC0EFAAAUq25E8eNWLiWCCgAAsBZBBQAAWIugAgAArEVQAQAA1iKoAAAAaxFUAACAtQgqAADAWgQVAABgLYIKAACwFkEFAABYi6ACAACsRVABAADWquR0AQAAwD6BvpW0e3yi02XQowIAAOxFUAEAANYiqAAAAGsRVAAAgLUIKgAAwFoEFQAAYC2CCgAAsBZBBQAAWIugAgAArEVQAQAA1iKoAAAAaxFUAACAtQgqAADAWgQVAABgLYIKAACwFkEFAABYi6ACAACsRVABAADWIqgAAABrEVQAAIC1CCoAAMBaBBUAAGAtggoAALAWQQUAAFiLoAIAAKxFUAEAANYiqAAAAGsRVAAAgLUIKgAAwFoEFQAAYC2CCgAAsBZBBQAAWIugAgAArOVoUJk2bZpiYmIUGhqq0NBQxcXF6auvvnKyJAAAYBFHg0qtWrU0fvx4rV+/XuvWrdONN96om2++WT/++KOTZQEAAEu4jDHG6SLOFB4ergkTJmjAgAHnXTYjI0NhYWFKT09XaGjoZagOAABcrNJ8fle6TDWdV35+vubNm6fs7GzFxcU5XQ4AALCA40Fl06ZNiouL08mTJxUcHKz58+erSZMmRS6bm5ur3Nxcz3RGRsblKhMAADjA8at+rrvuOm3cuFHffPONHnroIfXr10+bN28uctnk5GSFhYV5HlFRUZe5WgAAcDlZN0alc+fOatCggaZPn15oXlE9KlFRUYxRAQDgCnJFjlEp4Ha7vcLImfz8/OTn53eZKwIAAE5xNKiMGjVK3bp1U+3atZWZmak5c+YoNTVVixYtcrIsAABgCUeDysGDB3Xvvfdq//79CgsLU0xMjBYtWqQuXbo4WRYAALCEo0FlxowZTu4eAABYzvGrfgAAAM6FoAIAAKxFUAEAANYiqAAAAGsRVAAAgLUIKgAAwFoEFQAAYC2CCgAAsBZBBQAAWIugAgAArEVQAQAA1iKoAAAAaxFUAACAtQgqAADAWgQVAABgLYIKAACwFkEFAABYi6ACAACsRVABAADWIqgAAABrEVQAAIC1CCoAYJmcvNOq+8QXqvvEF8rJO+10OYCjCCoAAMBaBBUAAGAtggoAALAWQQUAAFiLoAIAAKxFUAEAANYiqAAAAGsRVAAAgLUIKgAAwFoEFQAAYC2CCgBYbPfhHKdLABxFUAEAixzPydPA2es9091fXaF7Z6xRes4pB6sCnENQAQCLPPL+Rq366bBX28qdhzXk/Q0OVQQ4i6ACAJb4+VCWlu84JPdZ7fnGaPmOQ9p1ONuRugAnEVQAwBK/HC1+PMruIwQVlD8EFQCwRJ3wwGLn160adJkqAexBUAEAS9SvFqyEa6oV+sNc0eVSwjXVVC+CoILyh6ACABaZ1Kel4hpEeLXFN4zQpD4tHaoIcBZBBQAsEhboozf6xXqmv3ykvd4e0FZhgT4OVgU4h6ACABarG1H8uBXgj46gAgAArEVQAQAA1iKoAAAAaxFUAACAtQgqAADAWgQVAABgLYIKAACwFkEFAABYi6ACAACsVcnpAgAA3gJ9K2n3+ESnywCsQI8KAACwFkEFAABYi6ACAACsRVABAADWIqgAAABrEVQAAIC1CCoAAMBaBBUAAGAtggoAALAWQQUAAFiLoAIAAKxFUAEAANYiqAAAAGsRVAAAgLUIKgAAwFoEFQAAYC2CCgAAsBZBBQAAWIugAgAArFWqoPLhhx8qLy/PM/3bb7/J7XZ7pnNycvTCCy+UeHvJyclq06aNQkJCFBkZqVtuuUXbtm0rTUkAAOAPrFRBpU+fPjp+/LhnukmTJtq9e7dnOjMzU6NGjSrx9pYtW6akpCStXr1aixcv1qlTp3TTTTcpOzu7NGUBAIA/qEqlWdgYU+x0aS1cuNBr+q233lJkZKTWr1+vhISEi9o2AAC48pUqqFxq6enpkqTw8PAi5+fm5io3N9cznZGRcVnqAgAAzrBmMK3b7dawYcMUHx+vpk2bFrlMcnKywsLCPI+oqKjLXCUAALicSt2jsmjRIoWFhUn6PVwsXbpUP/zwgyR5jV8praSkJP3www/6+uuvz7nMqFGjNHz4cM90RkYGYQUAgD8wlynFQJMKFc7fAeNyuZSfn1+qIgYPHqxPP/1Uy5cvV7169Uq8XkZGhsLCwpSenq7Q0NBS7RMAADijNJ/fpepROfNS5LJgjNGQIUM0f/58paamliqkAACAP75SjVG5//77lZmZWWY7T0pK0rvvvqs5c+YoJCREaWlpSktL04kTJ8psHwAA4MpVqlM/FStW1P79+xUZGVk2O3e5imyfNWuW7rvvvvOuz6kfAACuPJfs1M/F3jflUm8PAAD8sZT6qp/MzEz5+/sXuwy9GwAAoCyUOqhce+2155xnjLmgq34AlJ2cvNNqMnqRJGnzM10V6GvVfR0BoFRK/Rfso48+OuedYwEAAMpSqYNKfHx8mQ2mBQAAKE6Z30Kf0z4AAKCslCqo1KlTRxUrVixy3vbt2zVixAjVqlWrTAoDAAAoVVDZtWuXqlat6pnOycnRrFmz1L59ezVp0kTLly/3+l88AAAAF+OCLgdYvXq13nzzTc2bN0+1a9fWli1blJKSovbt25d1fQAAoBwrVY/Kiy++qOjoaPXq1UtVqlTR8uXLtWnTJrlcLq+eFgAAgLJQqh6VkSNHauTIkXrmmWfOOVYFAACgrJSqR+Uf//iH5s2bp3r16mnkyJH64YcfLlVdAAAApQsqo0aN0vbt2/XOO+8oLS1N7dq1U/PmzWWM0bFjxy5VjQAAoJy6oPuodOjQQbNnz9b+/fv18MMPq1WrVkpISND111+vl156qaxrBAAA5dRF3fAtNDRUDzzwgNasWaPvvvtO7dq10/jx48uqNgAAUM6VajDtiRMntHTpUv31r3+V9PupoNzc3P9trFIl/fTTT2VbIQAAKLdKFVRmz56tL774whNUJk+erOjoaAUEBEiStm3bppo1a+rvf/972VcKAADKnVKd+nnvvfc0aNAgr7Y5c+YoJSVFKSkpeuGFF/Thhx+WaYEAAKD8KlVQ2blzp5o1a+aZ9vf3V4UK/9tE27ZttXnz5rKrDgAAlGulOvVz/PhxrzEphw4d8prvdru95gMAAFyMUvWo1KpVq9ibvH3//ff892QAAFBmShVUunfvrtGjR+vkyZOF5p04cULjxo1TYmJimRUH4OLsPpzjdAkAcFFcxhhT0oUPHDigFi1ayNfXV4MHD9a1114r6ferfSZPnqzTp09rw4YNql69+iUr+EwZGRkKCwtTenq6QkNDL8s+AZsdz8lT0nsbtPKnw562hGuqaVKflgoL9HGwMgD4n9J8fpcqqEjSrl279NBDD2nx4sUqWNXlcqlLly6aOnWq6tevf+GVlxJBBfB274w1+nrHIbnPaKvocim+YYTeHtDWsboA4Eyl+fwu1WBaSapXr54WLlyoo0ePaufOnZKkhg0bKjw8/MKqBVAmfj6UpeU7DhVqzzdGy3cc0q7D2aoXEeRAZQBw4UodVAqEh4erbVu+oQG2+OVo8eNRdh8hqAC48lzU//oBYI864YHFzq9blZAC4MpDUAH+IOpXC1bCNdUK/VJXdLmUcE01elMAXJEIKsAfyKQ+LRXXIMKrLb5hhCb1aelQRQBwcQgqwB9IWKCP3ugX65n+8pH2entAWy5NBnDFIqgAf2B1I4oftwIAtiOoAAAAaxFUAACAtQgqAADAWgQVAABgLYIKAACwFkEFAABYi6ACAACsRVABAADWIqgAAABrEVQAAIC1CCoAAMBaBBUAAGAtggoAALAWQQUAAFiLoAIAAKxFUAEAANaq5HQBAMpWoG8l7R6f6HQZAFAm6FEBAADWIqgAAABrEVQAAIC1CCoAAMBaBBUAAGAtggoAALAWQQUAAFiLoAIAAKxFUAEAANYiqAAAAGsRVAAAgLUIKgAAwFoEFQAAYC2CCgAAsBZBBQAAWIugAgAArEVQAQAA1iKoAAAAaxFUAACAtQgqAADAWgQVAABgLYIKAACwFkEFAABYy9Ggsnz5cvXo0UM1a9aUy+XSggULnCwHAABYxtGgkp2drebNm2vKlClOlgEAACxVycmdd+vWTd26dXOyBAAAYDHGqAAAAGs52qNSWrm5ucrNzfVMZ2RkOFgNAAC41K6oHpXk5GSFhYV5HlFRUU6XBAAALqErKqiMGjVK6enpnseePXucLgkAAFxCV9SpHz8/P/n5+TldBgAAuEwcDSpZWVnauXOnZ3rXrl3auHGjwsPDVbt2bQcrAwAANnA0qKxbt04dO3b0TA8fPlyS1K9fP7311lsOVQUAAGzhaFC54YYbZIxxsgQAAGCxK2owLQAAKF8IKgAAwFoEFQAAYC2CCgAAsBZBBQAAWIugAgAArEVQAQAA1iKoAAAAaxFUAACAtQgqAADAWgQVAABgLYIKAACwFkEFAABYi6ACAACsRVABAADWIqgAAABrEVQAAIC1CCoAAMBaBBUAAGAtggoAALAWQQUAAFiLoAIAAKxFUAEAANYiqAAAAGsRVAAAgLUIKgAAwFoEFQAAYC2CShFy8k6r7hNfqO4TXygn77TT5QAAUG4RVAAAgLUIKgAAwFoEFQAAYC2CCgAAsBZBBQAAWIugAgAArEVQAQAA1iKoAAAAaxFUAACAtQgqAADAWgQVAABgLYIKAACwFkEFAABYi6ACAACsRVABAADWIqgAAABrEVQAAIC1CCoAAMBaBBUAAGAtggoAALAWQQUAAFiLoAIAAKxFUAEAANYiqAAAAGsRVAAAgLUIKgAAwFoEFQAAYC2CynnsPpzjdAkAAJRbBJWzHM/J08DZ6z3T3V9doXtnrFF6zikHqwIAoHwiqJzlkfc3atVPh73aVu48rCHvb3CoIgAAyi+Cyhl+PpSl5TsOyX1We74xWr7jkHYdznakLgAAyiuCyhl+OVr8eJTdRwgqAABcTgSVM9QJDyx2ft2qQZepEgAAIBFUvNSvFqyEa6oVelEqulxKuKaa6kUQVAAAuJwIKmeZ1Kel4hpEeLXFN4zQpD4tHaoIAIDyi6BylrBAH73RL9Yz/eUj7fX2gLYKC/RxsCoAAMongsp51I0oftwKAAC4dAgqAADAWgQVAABgLYIKAACwFkEFAABYi6ACAACsRVABAADWsiKoTJkyRXXr1pW/v7/atWunNWvWOF0SAACwgONBZe7cuRo+fLjGjBmjb7/9Vs2bN1fXrl118OBBp0sDAAAOczyovPTSSxo4cKD69++vJk2a6LXXXlNgYKBmzpzpdGkAAMBhjgaVvLw8rV+/Xp07d/a0VahQQZ07d9aqVasKLZ+bm6uMjAyvBwAA+ONyNKgcPnxY+fn5ql69uld79erVlZaWVmj55ORkhYWFeR5RUVGXq1QAAOAAx0/9lMaoUaOUnp7ueezZs8fpkgAAwCVUycmdR0REqGLFijpw4IBX+4EDB1SjRo1Cy/v5+cnPz+9ylQcAABzmaI+Kr6+vYmNjtXTpUk+b2+3W0qVLFRcX52BlAADABo72qEjS8OHD1a9fP7Vu3Vpt27bVyy+/rOzsbPXv39/p0gAAgMMcDyp33HGHDh06pNGjRystLU0tWrTQwoULCw2wBQAA5Y/jQUWSBg8erMGDBztdBgAAsMwVddUPAAAoXwgqAADAWgQVAABgLYIKAACwFkEFAABYi6ACAACsRVABAADWIqgAAABrEVQAAIC1CCoAAMBaVtxC3zaBvpW0e3yi02UAAFDu0aMCAACsRVABAADWIqgAAABrEVQAAIC1CCoAAMBaBBUAAGAtggoAALAWQQUAAFiLoAIAAKxFUAEAANYiqAAAAGsRVAAAgLUIKgAAwFoEFQAAYC2CCgAAsFYlpwu4GMYYSVJGRobDlQAAgJIq+Nwu+BwvzhUdVDIzMyVJUVFRDlcCAABKKzMzU2FhYcUu4zIliTOWcrvd2rdvn0JCQuRyuZwu54qTkZGhqKgo7dmzR6GhoU6XUy7xHtiB98F5vAfOu5zvgTFGmZmZqlmzpipUKH4UyhXdo1KhQgXVqlXL6TKueKGhofxhcBjvgR14H5zHe+C8y/UenK8npQCDaQEAgLUIKgAAwFoElXLMz89PY8aMkZ+fn9OllFu8B3bgfXAe74HzbH0PrujBtAAA4I+NHhUAAGAtggoAALAWQQUAAFiLoAIAAKxFUCnHpkyZorp168rf31/t2rXTmjVrnC6pXFm+fLl69OihmjVryuVyacGCBU6XVK4kJyerTZs2CgkJUWRkpG655RZt27bN6bLKnWnTpikmJsZzk7G4uDh99dVXTpdVro0fP14ul0vDhg1zuhRJBJVya+7cuRo+fLjGjBmjb7/9Vs2bN1fXrl118OBBp0srN7Kzs9W8eXNNmTLF6VLKpWXLlikpKUmrV6/W4sWLderUKd10003Kzs52urRypVatWho/frzWr1+vdevW6cYbb9TNN9+sH3/80enSyqW1a9dq+vTpiomJcboUDy5PLqfatWunNm3aaPLkyZJ+/79JUVFRGjJkiJ544gmHqyt/XC6X5s+fr1tuucXpUsqtQ4cOKTIyUsuWLVNCQoLT5ZRr4eHhmjBhggYMGOB0KeVKVlaWWrVqpalTp+qf//ynWrRooZdfftnpsuhRKY/y8vK0fv16de7c2dNWoUIFde7cWatWrXKwMsA56enpkn7/kIQz8vPz9cEHHyg7O1txcXFOl1PuJCUlKTEx0euzwQZX9D8lxIU5fPiw8vPzVb16da/26tWra+vWrQ5VBTjH7XZr2LBhio+PV9OmTZ0up9zZtGmT4uLidPLkSQUHB2v+/Plq0qSJ02WVKx988IG+/fZbrV271ulSCiGoACj3kpKS9MMPP+jrr792upRy6brrrtPGjRuVnp6ujz76SP369dOyZcsIK5fJnj17NHToUC1evFj+/v5Ol1MIQaUcioiIUMWKFXXgwAGv9gMHDqhGjRoOVQU4Y/Dgwfr3v/+t5cuXq1atWk6XUy75+vqqYcOGkqTY2FitXbtWr7zyiqZPn+5wZeXD+vXrdfDgQbVq1crTlp+fr+XLl2vy5MnKzc1VxYoVHauPMSrlkK+vr2JjY7V06VJPm9vt1tKlSzkvjHLDGKPBgwdr/vz5+n//7/+pXr16TpeE/5/b7VZubq7TZZQbnTp10qZNm7Rx40bPo3Xr1urbt682btzoaEiR6FEpt4YPH65+/fqpdevWatu2rV5++WVlZ2erf//+TpdWbmRlZWnnzp2e6V27dmnjxo0KDw9X7dq1HaysfEhKStKcOXP06aefKiQkRGlpaZKksLAwBQQEOFxd+TFq1Ch169ZNtWvXVmZmpubMmaPU1FQtWrTI6dLKjZCQkEJjs4KCglS1alUrxmwRVMqpO+64Q4cOHdLo0aOVlpamFi1aaOHChYUG2OLSWbdunTp27OiZHj58uCSpX79+euuttxyqqvyYNm2aJOmGG27wap81a5buu+++y19QOXXw4EHde++92r9/v8LCwhQTE6NFixapS5cuTpcGS3AfFQAAYC3GqAAAAGsRVAAAgLUIKgAAwFoEFQAAYC2CCgAAsBZBBQAAWIugAgAArEVQAQAA1iKoAFeQtLQ0DR06VA0bNpS/v7+qV6+u+Ph4TZs2TTk5OU6Xd05169aVy+XSBx98UGhedHS0XC6XNXfjTUtL05AhQ1S/fn35+fkpKipKPXr08PrfWAWSk5NVsWJFTZgwwdNW8FzP9Si46+255hf1GgHlGbfQB64QP//8s+Lj41W5cmU999xzatasmfz8/LRp0ya9/vrruvrqq9WzZ88i1z116pR8fHwuc8XeoqKiNGvWLN15552ettWrVystLU1BQUEOVvY/u3fv9rzGEyZMULNmzXTq1CktWrRISUlJ2rp1q9fyM2fO1IgRIzRz5kw9/vjjkqS1a9cqPz9fkvTf//5Xf/vb37Rt2zaFhoZKktf/EZo1a5b+8pe/eG2zcuXKl/AZAlcgA+CK0LVrV1OrVi2TlZVV5Hy32+35WZKZOnWq6dGjhwkMDDRjxowxxhgzdepUU79+fePj42OuvfZa8/bbb3vW2bVrl5FkNmzY4Gk7duyYkWRSUlKMMcakpKQYSebf//63adasmfHz8zPt2rUzmzZtKrb2OnXqmCeeeML4+fmZX3/91dM+cOBAM2TIEBMWFmZmzZrltd8BAwaYiIgIExISYjp27Gg2btzomb9z507Ts2dPExkZaYKCgkzr1q3N4sWLC+3z2WefNf379zfBwcEmKirKTJ8+vdg6u3XrZq6++uoiX+Njx455Taemppqrr77a5OXlmZo1a5qVK1cWWqfg9Tp7XWN+f4/mz59fbD0AjOHUD3AFOHLkiP7zn/8oKSnpnL0PLpfLa3rs2LG69dZbtWnTJt1///2aP3++hg4dqkcffVQ//PCDHnjgAfXv318pKSmlrufxxx/Xiy++qLVr16patWrq0aOHTp06Vew61atXV9euXTV79mxJUk5OjubOnav777+/0LK33367Dh48qK+++krr169Xq1at1KlTJx09elTS7/95unv37lq6dKk2bNigv/zlL+rRo4d+/fVXr+28+OKLat26tTZs2KCHH35YDz30kLZt21ZkfUePHtXChQvP+Rqf3dMxY8YM9enTRz4+PurTp49mzJhR7PMHcIGcTkoAzm/16tVGkvnkk0+82qtWrWqCgoJMUFCQGTFihKddkhk2bJjXstdff70ZOHCgV9vtt99uunfvbowpXY/KBx984FnmyJEjJiAgwMydO/ec9depU8dMnDjRLFiwwDRo0MC43W4ze/Zs07JlS2OM8epRWbFihQkNDTUnT5702kaDBg2K7RGJjo42kyZN8trn3Xff7Zl2u90mMjLSTJs2rcj1v/nmmyJf46Kkp6ebgIAATy/Phg0bTHBwsMnMzPRa7nw9Kv7+/p73r+Dxyy+/nHf/QHlCjwpwBVuzZo02btyo6Oho5ebmes1r3bq11/SWLVsUHx/v1RYfH68tW7aUer9xcXGen8PDw3XdddeVaDuJiYnKysrS8uXLNXPmzCJ7U7777jtlZWWpatWqCg4O9jx27dqln376SdLvPSqPPfaYGjdurMqVKys4OFhbtmwp1KMSExPj+dnlcqlGjRo6ePBgkbWZUvwj+ffff18NGjRQ8+bNJUktWrRQnTp1NHfu3BJvQ5ImTpyojRs3ej1q1qxZqm0Af3QMpgWuAA0bNpTL5Sp02qJ+/fqSvAdoFijtANUKFX7/3nLmB/b5TueUVqVKlXTPPfdozJgx+uabbzR//vxCy2RlZemqq65SampqoXkFp18ee+wxLV68WP/617/UsGFDBQQEqFevXsrLy/Na/uwBxC6XS263u8jarrnmGrlcrkIDZosyY8YM/fjjj6pU6X9/Qt1ut2bOnKkBAwacd/0CNWrUUMOGDUu8PFAe0aMCXAGqVq2qLl26aPLkycrOzr6gbTRu3FgrV670alu5cqWaNGkiSapWrZokaf/+/Z75GzduLHJbq1ev9vx87Ngxbd++XY0bNy5RHffff7+WLVumm2++WVWqVCk0v1WrVkpLS1OlSpXUsGFDr0dERISn7vvuu0+33nqrmjVrpho1amj37t0l2v+5hIeHq2vXrpoyZUqRr/Hx48clSZs2bdK6deuUmprq1ROSmpqqVatWlSjoACg5elSAK8TUqVMVHx+v1q1ba+zYsYqJiVGFChW0du1abd26VbGxscWu//jjj6t3795q2bKlOnfurM8//1yffPKJlixZIun3Xpk//elPGj9+vOrVq6eDBw/q//7v/4rc1jPPPKOqVauqevXqeuqppxQREaFbbrmlRM+jcePGOnz4sAIDA4uc37lzZ8XFxemWW27RCy+8oGuvvVb79u3TF198oVtvvVWtW7fWNddco08++UQ9evSQy+XS008/fc6ektKYMmWK4uPj1bZtWz3zzDOKiYnR6dOntXjxYk2bNk1btmzRjBkz1LZtWyUkJBRav02bNpoxY4bXfVWKc/z4caWlpXm1hYSEWHO5NmAFpwfJACi5ffv2mcGDB5t69eoZHx8fExwcbNq2bWsmTJhgsrOzPcvpHJe+Fnd5sjHGbN682cTFxZmAgADTokUL85///KfIwbSff/65iY6ONr6+vqZt27bmu+++K7bugsG053L25ckZGRlmyJAhpmbNmsbHx8dERUWZvn37ei5t3rVrl+nYsaMJCAgwUVFRZvLkyaZDhw5m6NChxe6zefPmnku1z2Xfvn0mKSnJ1KlTx/j6+pqrr77a9OzZ06SkpJjc3FxTtWpV88ILLxS57vPPP28iIyNNXl6eMeb8g2mLeiQnJxdbH1DeuIwpxQgyAOVaamqqOnbsqGPHjnFjMgCXBWNUAACAtQgqAADAWpz6AQAA1qJHBQAAWIugAgAArEVQAQAA1iKoAAAAaxFUAACAtQgqAADAWgQVAABgLYIKAACwFkEFAABY6/8DQu9UKqAXqh0AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "res_t.cal.plot_cal(2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.9.17" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file