diff --git a/docs/image/lewis.png b/docs/image/lewis.png new file mode 100644 index 0000000..10ec91a Binary files /dev/null and b/docs/image/lewis.png differ diff --git a/docs/image/lewis.svg b/docs/image/lewis.svg new file mode 100644 index 0000000..267644a --- /dev/null +++ b/docs/image/lewis.svg @@ -0,0 +1,110 @@ + + + + + + + + + +x0 + +x0 + + + +x1 + +x1 + + + +x0->x1 + + +0.40 + + + +x4 + +x4 + + + +x1->x4 + + +0.20 + + + +y + +y + + + +x1->y + + +0.50 + + + +x2 + +x2 + + + +x2->x1 + + +-0.60 + + + +x3 + +x3 + + + +x3->y + + +-0.10 + + + +x4->y + + +-0.50 + + + +x5 + +x5 + + + +x5->x3 + + +0.50 + + + +x5->x4 + + +-0.70 + + + diff --git a/docs/reference/index.rst b/docs/reference/index.rst index 1f53a4c..3d8f32e 100644 --- a/docs/reference/index.rst +++ b/docs/reference/index.rst @@ -31,5 +31,6 @@ API Reference missingness_lingam abic_lingam causal_effect + lewis utils tools diff --git a/docs/reference/lewis.rst b/docs/reference/lewis.rst new file mode 100644 index 0000000..ea2cc1d --- /dev/null +++ b/docs/reference/lewis.rst @@ -0,0 +1,8 @@ +.. module:: lingam + +LEWIS +===== + +.. autoclass:: LEWIS + :members: + :inherited-members: \ No newline at end of file diff --git a/docs/tutorial/index.rst b/docs/tutorial/index.rst index b084255..5f513ec 100644 --- a/docs/tutorial/index.rst +++ b/docs/tutorial/index.rst @@ -46,3 +46,4 @@ Contents: high_dim_direct_lingam missingness_lingam abic_lingam + lewis diff --git a/docs/tutorial/lewis.rst b/docs/tutorial/lewis.rst new file mode 100644 index 0000000..071b54a --- /dev/null +++ b/docs/tutorial/lewis.rst @@ -0,0 +1,705 @@ +Causal Explanations with LEWIS +============================== + +This notebook implements **LEWIS**, a causal explanation framework for +black-box machine learning models. Using probabilistic contrastive +counterfactuals, we compute **Necessity (Nec)**, **Sufficiency (Suf)**, +and **Necessity-and-Sufficiency (NeSuf)** scores to quantify causal +feature importance without relying on model internals. This notebook +follows the methodology proposed by Galhotra et al. (SIGMOD 2021). + +**Reference:** + +* Sainyam Galhotra, Romila Pradhan, Babak Salimi (2021). Explaining Black-Box Algorithms Using Probabilistic Contrastive Counterfactuals. SIGMOD ’21: International Conference on Management of Data, Virtual Event, China, June 20-25, 2021. + +.. code-block:: python + + import numpy as np + import pandas as pd + import operator + import graphviz + import matplotlib.pyplot as plt + + from sklearn.preprocessing import KBinsDiscretizer + from sklearn.model_selection import train_test_split + from sklearn.ensemble import RandomForestClassifier + + import lingam + from lingam.utils import make_dot + + print([np.__version__, pd.__version__, graphviz.__version__, lingam.__version__]) + + np.set_printoptions(precision=3, suppress=True) + np.random.seed(0) + + +.. parsed-literal:: + + ['1.26.4', '2.3.0', '0.20.3', '1.12.1'] + + +Test data +--------- + +We create test data consisting of 7 variables. + +.. code-block:: python + + m = np.array([ + [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [ 0.4, 0.0,-0.6, 0.0, 0.0, 0.0, 0.0], + [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0], + [ 0.0, 0.2, 0.0, 0.0, 0.0,-0.7, 0.0], + [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [ 0.0, 0.5, 0.0,-0.1,-0.5, 0.0, 0.0], + ]) + + generate_error = lambda p: np.random.uniform(-p, p, size=1000) + + error_vars = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5] + params = [0.5 * np.sqrt(12 * v) for v in error_vars] + e = np.array([generate_error(p) for p in params]) + + cols = [f"x{i}" for i in range(len(m) - 1)] + ["y"] + X = np.linalg.pinv(np.eye(len(m)) - m) @ e + df = pd.DataFrame(X.T, columns=cols) + + display(make_dot(m, labels=cols)) + df.head() + + + +.. image:: ../image/lewis.svg + + + + +.. raw:: html + +
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
x0x1x2x3x4x5y
00.119568-0.1825000.763061-0.373908-0.315997-0.326320-0.202553
10.527104-0.954103-0.0585820.0696400.319117-0.495717-0.818257
20.2517180.0074410.0567200.154034-0.146964-1.0567110.510195
30.1099410.922016-0.6110970.6805210.069137-0.3612320.605658
4-0.187007-1.3462110.2573020.447058-0.904446-0.655983-1.317686
+
+
+ + +Discretization +-------------- + +LEWIS assumes that all features have finite, discrete domains. +Therefore, continuous-valued features must be discretized before +computing LEWIS scores. Each continuous feature is binned into a small +number of ordered categories (e.g., Low < Medium < High), enabling +contrastive value pairs :math:`(x, x′)` and well-defined computation of +Necessity, Sufficiency, and Necessity-and-Sufficiency scores. + +.. code-block:: python + + kbd = KBinsDiscretizer(n_bins=2, encode="ordinal", strategy="uniform", subsample=None) + df["y"] = kbd.fit_transform(df["y"].values.reshape(-1, 1)) + + feature_names = [f"x{i}" for i in range(len(m) - 1)] + for name in feature_names: + kbd = KBinsDiscretizer(n_bins=5, encode="ordinal", strategy="uniform", subsample=None) + df[name] = kbd.fit_transform(df[name].values.reshape(-1, 1)) + + df.head() + + + + +.. raw:: html + +
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
x0x1x2x3x4x5y
02.02.04.01.02.01.00.0
13.01.02.02.02.01.00.0
23.02.02.02.02.00.01.0
32.03.01.03.02.01.01.0
42.01.03.03.01.01.00.0
+
+
+ + +Model Training and Prediction with RandomForestClassifier +--------------------------------------------------------- + +We train a black-box prediction model using a RandomForestClassifier. +The model learns a non-linear decision function by aggregating +predictions from multiple decision trees trained on bootstrapped samples +of the data. After training, the model is used to generate predictions +(and class probabilities) for each instance. These predicted outcomes +serve as the decision variable used in subsequent LEWIS analysis, while +the internal structure of the model remains opaque. + +.. code-block:: python + + X = df[feature_names] + y = df["y"] + + X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, random_state=0) + + clf = RandomForestClassifier(max_depth=5, random_state=0) + clf.fit(X_train, y_train) + + y_pred = clf.predict(X_test) + X_test["y"] = y_pred + X_test.head() + + + + +.. raw:: html + +
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
x0x1x2x3x4x5y
9932.02.04.02.01.02.00.0
8591.00.04.03.02.01.00.0
2984.03.04.01.04.01.00.0
5532.01.04.02.00.04.00.0
6723.02.02.01.01.02.01.0
+
+
+ + +Computing LEWIS Scores +---------------------- + +We compute LEWIS explanation scores—**Necessity (Nec)**, **Sufficiency +(Suf)**, and **Necessity-and-Sufficiency (NeSuf)**—for each selected +feature using a causal, model-agnostic approach. Based on the given +causal graph, we first identify an adjustment set that satisfies the +**backdoor criterion** for the relationship between the target feature +and the outcome. These variables are used only for probability +estimation and are not part of the user-facing explanation. + +For an ordinal feature with ordered values +:math:`(x^{(1)} \prec \dots \prec x^{(L)})`, we compute Nec, Suf, and +NeSuf for **all ordered value pairs** :math:`(x, x')` such that +:math:`(x \succ x')`. After evaluating all pairs, we select the pair +that **maximizes NeSuf** as the most causally responsive contrast. The +corresponding Nec and Suf values for this same pair are then reported, +ensuring that all three scores are aligned to a single, interpretable +intervention. + +.. code-block:: python + + backdoor = { + 'x0': [], + 'x1': [], + 'x2': [], + 'x3': ['x5'], + 'x4': ['x1', 'x5'], + 'x5': [], + } + +.. code-block:: python + + def max_nesuf(score_dic, score, name): + pre_score = [0, 0, 0] + if name in score_dic.keys(): + pre_score = score_dic[name] + if score[2] > pre_score[2]: # nesuf + score_dic[name] = score + return score_dic + +.. code-block:: python + + score_dic = {} + for name in feature_names: + val = {} + for v in df[name].unique(): + df_temp = df[df[name] == v] + vc = df_temp["y"].value_counts().reindex([0, 1], fill_value=0) + val[v] = vc[1] * 1.0 / (vc[0] + vc[1]) + sorted_val = sorted(val.items(), key=operator.itemgetter(1)) + + for i, (before, _) in enumerate(sorted_val): + for after, _ in sorted_val[i + 1 :]: + lewis = lingam.LEWIS() + score = lewis.get_scores( + X_test, + x_names=[name], + x_values=[after], + x_prime_values=[before], + o_name="y", + c_names=backdoor[name]) + score_dic = max_nesuf(score_dic, score, name) + print(f"{name}: {before}->{after}: score={score}") + + + +.. parsed-literal:: + + x0: 0.0->1.0: score=(0.28876449588651437, 0.15671962820931234, 0.11307290998254543) + x0: 0.0->2.0: score=(0.37054073719681124, 0.22722774483181898, 0.16394438035923292) + x0: 0.0->3.0: score=(0.308951540247843, 0.1725738099335212, 0.12451167156864434) + x0: 0.0->4.0: score=(0.3041121895598342, 0.16868933442058193, 0.12170902996578553) + x0: 1.0->2.0: score=(0.11497772655799332, 0.08361171323457248, 0.05087147037668749) + x0: 1.0->3.0: score=(0.028383066149784867, 0.018800605652118815, 0.011438761586098911) + x0: 1.0->4.0: score=(0.02157891947822522, 0.014194218923715932, 0.008636119983240098) + x0: 2.0->3.0: score=(0.0, 0.0, 0.0) + x0: 2.0->4.0: score=(0.0, 0.0, 0.0) + x0: 3.0->4.0: score=(0.0, 0.0, 0.0) + x1: 0.0->1.0: score=(1.0, 0.058517831098751374, 0.058517831098751374) + x1: 0.0->2.0: score=(1.0, 0.3863071221311797, 0.3863071221311797) + x1: 0.0->3.0: score=(1.0, 0.6916588068976024, 0.6916588068976024) + x1: 0.0->4.0: score=(1.0, 0.610516293905489, 0.610516293905489) + x1: 1.0->2.0: score=(0.8485199269018905, 0.348163036815634, 0.3277892910324283) + x1: 1.0->3.0: score=(0.9153949454338188, 0.6724938577835782, 0.633140975798851) + x1: 1.0->4.0: score=(0.9041502549843325, 0.5863079313025591, 0.5519984628067376) + x1: 2.0->3.0: score=(0.44147733206212025, 0.4975643286375128, 0.3053516847664227) + x1: 2.0->4.0: score=(0.36724518905144576, 0.3653442623497988, 0.22420917177430932) + x1: 3.0->4.0: score=(0.0, 0.0, 0.0) + x2: 4.0->3.0: score=(0.37781207464998834, 0.12785476560870568, 0.10561677393784988) + x2: 4.0->2.0: score=(0.5960755908256189, 0.31071620571720965, 0.2566728201473194) + x2: 4.0->0.0: score=(0.6241556609044814, 0.3496612983303619, 0.2888441281382971) + x2: 4.0->1.0: score=(0.6962790808358484, 0.48269334060441804, 0.39873768641480434) + x2: 3.0->2.0: score=(0.350799987082434, 0.2096685653922427, 0.1510560462094695) + x2: 3.0->0.0: score=(0.3959311587667224, 0.2543229315200742, 0.18322735420044717) + x2: 3.0->1.0: score=(0.5118501873958845, 0.406857207955013, 0.2931209124769545) + x2: 2.0->0.0: score=(0.06951813121731897, 0.056500809878571416, 0.03217130799097767) + x2: 2.0->1.0: score=(0.2480748569145521, 0.2495012015568321, 0.142064866267485) + x2: 0.0->1.0: score=(0.19189704999930096, 0.2045580893963687, 0.10989355827650732) + x3: 0.0->3.0: score=(0.0, 0.0, 0.07576064369254501) + x3: 0.0->2.0: score=(0.24019885850208111, 0.0, 0.06756468929060586) + x3: 0.0->4.0: score=(0.4375471879801787, 0.33627411063087126, 0.2394474103512867) + x3: 0.0->1.0: score=(0.5215173421020096, 0.2077705378411908, 0.2560619882491983) + x3: 3.0->2.0: score=(5.8380840876851114e-05, 0.0, 0.0016793883085030473) + x3: 3.0->4.0: score=(0.06087756395956861, 0.18097315391122015, 0.1735621093691839) + x3: 3.0->1.0: score=(0.4050632226945426, 0.2552560187517776, 0.1901766872670955) + x3: 2.0->4.0: score=(0.08389050626406856, 0.26072855163571107, 0.17188272106068087) + x3: 2.0->1.0: score=(0.39325969742974465, 0.27799227073061517, 0.18849729895859244) + x3: 4.0->1.0: score=(0.009866759685432035, 0.3007171788408852, 0.016614577897911605) + x4: 4.0->3.0: score=(0.0, 0.0002437752006358201, 0.01417600649350649) + x4: 4.0->2.0: score=(0.0, 0.43864694304558277, 0.3479847269469221) + x4: 4.0->1.0: score=(0.36470213549225855, 0.6477582075236524, 0.5560876613241783) + x4: 4.0->0.0: score=(0.8540240160406432, 0.6193767844978033, 0.624207928623127) + x4: 3.0->2.0: score=(0.0, 0.3662085477202382, 0.33380872045341564) + x4: 3.0->1.0: score=(0.34372689656444605, 0.6086744374259367, 0.5419116548306718) + x4: 3.0->0.0: score=(0.8449920368334224, 0.6353532663565686, 0.6100319221296207) + x4: 2.0->1.0: score=(0.0, 0.408105331786863, 0.21286296955041636) + x4: 2.0->0.0: score=(0.24400727130132333, 0.5372555633357295, 0.27622320167620484) + x4: 1.0->0.0: score=(0.17773403145390043, 0.252226191821934, 0.07080119345612482) + x5: 1.0->0.0: score=(0.0, 0.0, 0.0) + x5: 1.0->2.0: score=(0.2233835942176238, 0.08602778414260957, 0.06622184994605693) + x5: 1.0->3.0: score=(0.5867549758794374, 0.42466177557054957, 0.3268930922716893) + x5: 1.0->4.0: score=(0.6595932839891901, 0.5795249365522417, 0.44610254432152086) + x5: 0.0->2.0: score=(0.4502876323702336, 0.1594757752886144, 0.13348733208373778) + x5: 0.0->3.0: score=(0.7074927867488406, 0.4708967005285375, 0.3941585744093701) + x5: 0.0->4.0: score=(0.7590499242327206, 0.613314857297541, 0.5133680264592017) + x5: 2.0->3.0: score=(0.46789042692929883, 0.37050797119720974, 0.2606712423256323) + x5: 2.0->4.0: score=(0.5616797256969114, 0.5399476525078897, 0.3798806943754639) + x5: 3.0->4.0: score=(0.17625937121629423, 0.2691689069247335, 0.11920945204983158) + + +Visualization of LEWIS Scores Sorted by NeSuf +--------------------------------------------- + +Finally, we visualize the LEWIS scores using a horizontal bar chart, +where features are **sorted in descending order of NeSuf**. NeSuf +represents the strongest indicator of causal importance, capturing how +often a feature change would flip the decision outcome in both +directions. By sorting features by NeSuf, the plot highlights which +features are most causally influential overall. For each feature, the +corresponding **Nec and Suf values for the same maximizing contrast** +are displayed alongside NeSuf, enabling a consistent and interpretable +comparison. + +.. code-block:: python + + def plot_scores_barh_featurecol( + df, + feature_col="feature", + value_cols=("Nec", "Suf", "NeSuf"), + title="LEWIS Scores", + reverse_for_image=True, + ): + needed = [feature_col, *value_cols] + missing = [c for c in needed if c not in df.columns] + if missing: + raise ValueError(f"df is missing columns {missing}. Required columns: {needed}") + + df_plot = ( + df.loc[:, needed] + .dropna(subset=[feature_col]) + .drop_duplicates(subset=[feature_col]) + .set_index(feature_col) + .loc[:, list(value_cols)] + ) + + if reverse_for_image: + df_plot = df_plot.iloc[::-1] + + y = np.arange(len(df_plot.index)) + h = 0.22 + + fig, ax = plt.subplots(figsize=(3.0, 3.6), dpi=120) + + bars1 = ax.barh( + y + h, + df_plot[value_cols[0]].values, + height=h, + color="#4C72B0", + edgecolor="black", + linewidth=0.6, + label=value_cols[0], + ) + bars2 = ax.barh( + y, + df_plot[value_cols[1]].values, + height=h, + color="#F0E442", + edgecolor="black", + linewidth=0.6, + label=value_cols[1], + ) + bars3 = ax.barh( + y - h, + df_plot[value_cols[2]].values, + height=h, + color="#2CA02C", + edgecolor="black", + linewidth=0.6, + label=value_cols[2], + ) + + ax.set_yticks(y) + ax.set_yticklabels(df_plot.index) + ax.set_xlim(0, 1.0) + ax.set_xticks([0.0, 0.5, 1.0]) + ax.xaxis.tick_top() + ax.xaxis.set_label_position("top") + ax.set_xlabel(title) + ax.legend(loc="lower right", fontsize=7, frameon=True) + + def add_labels(bars): + for b in bars: + w = b.get_width() + yy = b.get_y() + b.get_height() / 2 + ax.text(w + 0.01, yy, f"{w:.2f}", va="center", ha="left", fontsize=7) + + add_labels(bars1) + add_labels(bars2) + add_labels(bars3) + + ax.grid(False) + plt.tight_layout() + plt.show() + +.. code-block:: python + + nec = [] + suf = [] + nesuf = [] + for name in feature_names: + score = score_dic[name] + nec.append(score[0]) + suf.append(score[1]) + nesuf.append(score[2]) + score_df = pd.DataFrame() + score_df['feature'] = feature_names + score_df['Nec'] = nec + score_df['Suf'] = suf + score_df['NeSuf'] = nesuf + + plot_scores_barh_featurecol(score_df.sort_values(['NeSuf'], ascending=False)) + + + +.. image:: ../image/lewis.png + + diff --git a/examples/LEWIS.ipynb b/examples/LEWIS.ipynb new file mode 100644 index 0000000..8bdff33 --- /dev/null +++ b/examples/LEWIS.ipynb @@ -0,0 +1,865 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "2b36f2e5-b5cc-4584-ae89-15e751ba6d3a", + "metadata": {}, + "source": [ + "# Causal Explanations with LEWIS\n", + "This notebook implements **LEWIS**, a causal explanation framework for black-box machine learning models. Using probabilistic contrastive counterfactuals, we compute **Necessity (Nec)**, **Sufficiency (Suf)**, and **Necessity-and-Sufficiency (NeSuf)** scores to quantify causal feature importance without relying on model internals. This notebook follows the methodology proposed by Galhotra et al. (SIGMOD 2021).\n", + "\n", + "**Reference:**\n", + "* Sainyam Galhotra, Romila Pradhan, Babak Salimi (2021). Explaining Black-Box Algorithms Using Probabilistic Contrastive Counterfactuals. SIGMOD '21: International Conference on Management of Data, Virtual Event, China, June 20-25, 2021." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "129b1ab6-37cb-4201-be5a-e7713d7f0c76", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['1.26.4', '2.3.0', '0.20.3', '1.12.1']\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import operator\n", + "import graphviz\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from sklearn.preprocessing import KBinsDiscretizer\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn.ensemble import RandomForestClassifier\n", + "\n", + "import lingam\n", + "from lingam.utils import make_dot\n", + "\n", + "print([np.__version__, pd.__version__, graphviz.__version__, lingam.__version__])\n", + "\n", + "np.set_printoptions(precision=3, suppress=True)\n", + "np.random.seed(0)" + ] + }, + { + "cell_type": "markdown", + "id": "5fb857e3-ba10-46a7-887e-545222bd1648", + "metadata": {}, + "source": [ + "## Test data\n", + "We create test data consisting of 7 variables." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "cd0cd8f6-6295-4e4e-9088-bbad46b5a070", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "x0\n", + "\n", + "x0\n", + "\n", + "\n", + "\n", + "x1\n", + "\n", + "x1\n", + "\n", + "\n", + "\n", + "x0->x1\n", + "\n", + "\n", + "0.40\n", + "\n", + "\n", + "\n", + "x4\n", + "\n", + "x4\n", + "\n", + "\n", + "\n", + "x1->x4\n", + "\n", + "\n", + "0.20\n", + "\n", + "\n", + "\n", + "y\n", + "\n", + "y\n", + "\n", + "\n", + "\n", + "x1->y\n", + "\n", + "\n", + "0.50\n", + "\n", + "\n", + "\n", + "x2\n", + "\n", + "x2\n", + "\n", + "\n", + "\n", + "x2->x1\n", + "\n", + "\n", + "-0.60\n", + "\n", + "\n", + "\n", + "x3\n", + "\n", + "x3\n", + "\n", + "\n", + "\n", + "x3->y\n", + "\n", + "\n", + "-0.10\n", + "\n", + "\n", + "\n", + "x4->y\n", + "\n", + "\n", + "-0.50\n", + "\n", + "\n", + "\n", + "x5\n", + "\n", + "x5\n", + "\n", + "\n", + "\n", + "x5->x3\n", + "\n", + "\n", + "0.50\n", + "\n", + "\n", + "\n", + "x5->x4\n", + "\n", + "\n", + "-0.70\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "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", + " \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", + "
x0x1x2x3x4x5y
00.119568-0.1825000.763061-0.373908-0.315997-0.326320-0.202553
10.527104-0.954103-0.0585820.0696400.319117-0.495717-0.818257
20.2517180.0074410.0567200.154034-0.146964-1.0567110.510195
30.1099410.922016-0.6110970.6805210.069137-0.3612320.605658
4-0.187007-1.3462110.2573020.447058-0.904446-0.655983-1.317686
\n", + "
" + ], + "text/plain": [ + " x0 x1 x2 x3 x4 x5 y\n", + "0 0.119568 -0.182500 0.763061 -0.373908 -0.315997 -0.326320 -0.202553\n", + "1 0.527104 -0.954103 -0.058582 0.069640 0.319117 -0.495717 -0.818257\n", + "2 0.251718 0.007441 0.056720 0.154034 -0.146964 -1.056711 0.510195\n", + "3 0.109941 0.922016 -0.611097 0.680521 0.069137 -0.361232 0.605658\n", + "4 -0.187007 -1.346211 0.257302 0.447058 -0.904446 -0.655983 -1.317686" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m = np.array([\n", + " [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],\n", + " [ 0.4, 0.0,-0.6, 0.0, 0.0, 0.0, 0.0],\n", + " [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],\n", + " [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0],\n", + " [ 0.0, 0.2, 0.0, 0.0, 0.0,-0.7, 0.0],\n", + " [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],\n", + " [ 0.0, 0.5, 0.0,-0.1,-0.5, 0.0, 0.0],\n", + "])\n", + "\n", + "generate_error = lambda p: np.random.uniform(-p, p, size=1000)\n", + "\n", + "error_vars = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]\n", + "params = [0.5 * np.sqrt(12 * v) for v in error_vars]\n", + "e = np.array([generate_error(p) for p in params])\n", + "\n", + "cols = [f\"x{i}\" for i in range(len(m) - 1)] + [\"y\"]\n", + "X = np.linalg.pinv(np.eye(len(m)) - m) @ e\n", + "df = pd.DataFrame(X.T, columns=cols)\n", + "\n", + "display(make_dot(m, labels=cols))\n", + "df.head()" + ] + }, + { + "cell_type": "markdown", + "id": "fc408161-1188-4d6b-9d2c-a2b5b8d5a710", + "metadata": {}, + "source": [ + "## Discretization\n", + "LEWIS assumes that all features have finite, discrete domains. Therefore, continuous-valued features must be discretized before computing LEWIS scores. Each continuous feature is binned into a small number of ordered categories (e.g., Low < Medium < High), enabling contrastive value pairs $(x, x′)$ and well-defined computation of Necessity, Sufficiency, and Necessity-and-Sufficiency scores." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "2e7c3d6b-addb-4534-bdfc-7d18c92f1c91", + "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", + " \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", + "
x0x1x2x3x4x5y
02.02.04.01.02.01.00.0
13.01.02.02.02.01.00.0
23.02.02.02.02.00.01.0
32.03.01.03.02.01.01.0
42.01.03.03.01.01.00.0
\n", + "
" + ], + "text/plain": [ + " x0 x1 x2 x3 x4 x5 y\n", + "0 2.0 2.0 4.0 1.0 2.0 1.0 0.0\n", + "1 3.0 1.0 2.0 2.0 2.0 1.0 0.0\n", + "2 3.0 2.0 2.0 2.0 2.0 0.0 1.0\n", + "3 2.0 3.0 1.0 3.0 2.0 1.0 1.0\n", + "4 2.0 1.0 3.0 3.0 1.0 1.0 0.0" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "kbd = KBinsDiscretizer(n_bins=2, encode=\"ordinal\", strategy=\"uniform\", subsample=None)\n", + "df[\"y\"] = kbd.fit_transform(df[\"y\"].values.reshape(-1, 1))\n", + "\n", + "feature_names = [f\"x{i}\" for i in range(len(m) - 1)]\n", + "for name in feature_names:\n", + " kbd = KBinsDiscretizer(n_bins=5, encode=\"ordinal\", strategy=\"uniform\", subsample=None)\n", + " df[name] = kbd.fit_transform(df[name].values.reshape(-1, 1))\n", + "\n", + "df.head()" + ] + }, + { + "cell_type": "markdown", + "id": "50a9ae11-ce51-4c38-aea3-6bc596144d31", + "metadata": {}, + "source": [ + "## Model Training and Prediction with RandomForestClassifier\n", + "We train a black-box prediction model using a RandomForestClassifier. The model learns a non-linear decision function by aggregating predictions from multiple decision trees trained on bootstrapped samples of the data. After training, the model is used to generate predictions (and class probabilities) for each instance. These predicted outcomes serve as the decision variable used in subsequent LEWIS analysis, while the internal structure of the model remains opaque." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "28fd0259-813e-47de-b02a-5b4062fb66b6", + "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", + " \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", + "
x0x1x2x3x4x5y
9932.02.04.02.01.02.00.0
8591.00.04.03.02.01.00.0
2984.03.04.01.04.01.00.0
5532.01.04.02.00.04.00.0
6723.02.02.01.01.02.01.0
\n", + "
" + ], + "text/plain": [ + " x0 x1 x2 x3 x4 x5 y\n", + "993 2.0 2.0 4.0 2.0 1.0 2.0 0.0\n", + "859 1.0 0.0 4.0 3.0 2.0 1.0 0.0\n", + "298 4.0 3.0 4.0 1.0 4.0 1.0 0.0\n", + "553 2.0 1.0 4.0 2.0 0.0 4.0 0.0\n", + "672 3.0 2.0 2.0 1.0 1.0 2.0 1.0" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "X = df[feature_names]\n", + "y = df[\"y\"]\n", + "\n", + "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, random_state=0)\n", + "\n", + "clf = RandomForestClassifier(max_depth=5, random_state=0)\n", + "clf.fit(X_train, y_train)\n", + "\n", + "y_pred = clf.predict(X_test)\n", + "X_test[\"y\"] = y_pred\n", + "X_test.head()" + ] + }, + { + "cell_type": "markdown", + "id": "28cccbe5-5f3f-40b3-8269-fb4e06083af4", + "metadata": {}, + "source": [ + "## Computing LEWIS Scores\n", + "We compute LEWIS explanation scores—**Necessity (Nec)**, **Sufficiency (Suf)**, and **Necessity-and-Sufficiency (NeSuf)**—for each selected feature using a causal, model-agnostic approach. Based on the given causal graph, we first identify an adjustment set that satisfies the **backdoor criterion** for the relationship between the target feature and the outcome. These variables are used only for probability estimation and are not part of the user-facing explanation.\n", + "\n", + "For an ordinal feature with ordered values $(x^{(1)} \\prec \\dots \\prec x^{(L)})$, we compute Nec, Suf, and NeSuf for **all ordered value pairs** $(x, x')$ such that $(x \\succ x')$. After evaluating all pairs, we select the pair that **maximizes NeSuf** as the most causally responsive contrast. The corresponding Nec and Suf values for this same pair are then reported, ensuring that all three scores are aligned to a single, interpretable intervention." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "4f8d2d6c-d758-4eb7-bfee-a5026f2a8f93", + "metadata": {}, + "outputs": [], + "source": [ + "backdoor = {\n", + " 'x0': [],\n", + " 'x1': [],\n", + " 'x2': [],\n", + " 'x3': ['x5'],\n", + " 'x4': ['x1', 'x5'],\n", + " 'x5': [],\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "4b4368a6-f845-4672-a8a4-17a92cc727f6", + "metadata": {}, + "outputs": [], + "source": [ + "def max_nesuf(score_dic, score, name):\n", + " pre_score = [0, 0, 0]\n", + " if name in score_dic.keys():\n", + " pre_score = score_dic[name]\n", + " if score[2] > pre_score[2]: # nesuf\n", + " score_dic[name] = score\n", + " return score_dic" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "54085376-f2b5-4571-aebe-fc7f3cb2544a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x0: 0.0->1.0: score=(0.28876449588651437, 0.15671962820931234, 0.11307290998254543)\n", + "x0: 0.0->2.0: score=(0.37054073719681124, 0.22722774483181898, 0.16394438035923292)\n", + "x0: 0.0->3.0: score=(0.308951540247843, 0.1725738099335212, 0.12451167156864434)\n", + "x0: 0.0->4.0: score=(0.3041121895598342, 0.16868933442058193, 0.12170902996578553)\n", + "x0: 1.0->2.0: score=(0.11497772655799332, 0.08361171323457248, 0.05087147037668749)\n", + "x0: 1.0->3.0: score=(0.028383066149784867, 0.018800605652118815, 0.011438761586098911)\n", + "x0: 1.0->4.0: score=(0.02157891947822522, 0.014194218923715932, 0.008636119983240098)\n", + "x0: 2.0->3.0: score=(0.0, 0.0, 0.0)\n", + "x0: 2.0->4.0: score=(0.0, 0.0, 0.0)\n", + "x0: 3.0->4.0: score=(0.0, 0.0, 0.0)\n", + "x1: 0.0->1.0: score=(1.0, 0.058517831098751374, 0.058517831098751374)\n", + "x1: 0.0->2.0: score=(1.0, 0.3863071221311797, 0.3863071221311797)\n", + "x1: 0.0->3.0: score=(1.0, 0.6916588068976024, 0.6916588068976024)\n", + "x1: 0.0->4.0: score=(1.0, 0.610516293905489, 0.610516293905489)\n", + "x1: 1.0->2.0: score=(0.8485199269018905, 0.348163036815634, 0.3277892910324283)\n", + "x1: 1.0->3.0: score=(0.9153949454338188, 0.6724938577835782, 0.633140975798851)\n", + "x1: 1.0->4.0: score=(0.9041502549843325, 0.5863079313025591, 0.5519984628067376)\n", + "x1: 2.0->3.0: score=(0.44147733206212025, 0.4975643286375128, 0.3053516847664227)\n", + "x1: 2.0->4.0: score=(0.36724518905144576, 0.3653442623497988, 0.22420917177430932)\n", + "x1: 3.0->4.0: score=(0.0, 0.0, 0.0)\n", + "x2: 4.0->3.0: score=(0.37781207464998834, 0.12785476560870568, 0.10561677393784988)\n", + "x2: 4.0->2.0: score=(0.5960755908256189, 0.31071620571720965, 0.2566728201473194)\n", + "x2: 4.0->0.0: score=(0.6241556609044814, 0.3496612983303619, 0.2888441281382971)\n", + "x2: 4.0->1.0: score=(0.6962790808358484, 0.48269334060441804, 0.39873768641480434)\n", + "x2: 3.0->2.0: score=(0.350799987082434, 0.2096685653922427, 0.1510560462094695)\n", + "x2: 3.0->0.0: score=(0.3959311587667224, 0.2543229315200742, 0.18322735420044717)\n", + "x2: 3.0->1.0: score=(0.5118501873958845, 0.406857207955013, 0.2931209124769545)\n", + "x2: 2.0->0.0: score=(0.06951813121731897, 0.056500809878571416, 0.03217130799097767)\n", + "x2: 2.0->1.0: score=(0.2480748569145521, 0.2495012015568321, 0.142064866267485)\n", + "x2: 0.0->1.0: score=(0.19189704999930096, 0.2045580893963687, 0.10989355827650732)\n", + "x3: 0.0->3.0: score=(0.0, 0.0, 0.07576064369254501)\n", + "x3: 0.0->2.0: score=(0.24019885850208111, 0.0, 0.06756468929060586)\n", + "x3: 0.0->4.0: score=(0.4375471879801787, 0.33627411063087126, 0.2394474103512867)\n", + "x3: 0.0->1.0: score=(0.5215173421020096, 0.2077705378411908, 0.2560619882491983)\n", + "x3: 3.0->2.0: score=(5.8380840876851114e-05, 0.0, 0.0016793883085030473)\n", + "x3: 3.0->4.0: score=(0.06087756395956861, 0.18097315391122015, 0.1735621093691839)\n", + "x3: 3.0->1.0: score=(0.4050632226945426, 0.2552560187517776, 0.1901766872670955)\n", + "x3: 2.0->4.0: score=(0.08389050626406856, 0.26072855163571107, 0.17188272106068087)\n", + "x3: 2.0->1.0: score=(0.39325969742974465, 0.27799227073061517, 0.18849729895859244)\n", + "x3: 4.0->1.0: score=(0.009866759685432035, 0.3007171788408852, 0.016614577897911605)\n", + "x4: 4.0->3.0: score=(0.0, 0.0002437752006358201, 0.01417600649350649)\n", + "x4: 4.0->2.0: score=(0.0, 0.43864694304558277, 0.3479847269469221)\n", + "x4: 4.0->1.0: score=(0.36470213549225855, 0.6477582075236524, 0.5560876613241783)\n", + "x4: 4.0->0.0: score=(0.8540240160406432, 0.6193767844978033, 0.624207928623127)\n", + "x4: 3.0->2.0: score=(0.0, 0.3662085477202382, 0.33380872045341564)\n", + "x4: 3.0->1.0: score=(0.34372689656444605, 0.6086744374259367, 0.5419116548306718)\n", + "x4: 3.0->0.0: score=(0.8449920368334224, 0.6353532663565686, 0.6100319221296207)\n", + "x4: 2.0->1.0: score=(0.0, 0.408105331786863, 0.21286296955041636)\n", + "x4: 2.0->0.0: score=(0.24400727130132333, 0.5372555633357295, 0.27622320167620484)\n", + "x4: 1.0->0.0: score=(0.17773403145390043, 0.252226191821934, 0.07080119345612482)\n", + "x5: 1.0->0.0: score=(0.0, 0.0, 0.0)\n", + "x5: 1.0->2.0: score=(0.2233835942176238, 0.08602778414260957, 0.06622184994605693)\n", + "x5: 1.0->3.0: score=(0.5867549758794374, 0.42466177557054957, 0.3268930922716893)\n", + "x5: 1.0->4.0: score=(0.6595932839891901, 0.5795249365522417, 0.44610254432152086)\n", + "x5: 0.0->2.0: score=(0.4502876323702336, 0.1594757752886144, 0.13348733208373778)\n", + "x5: 0.0->3.0: score=(0.7074927867488406, 0.4708967005285375, 0.3941585744093701)\n", + "x5: 0.0->4.0: score=(0.7590499242327206, 0.613314857297541, 0.5133680264592017)\n", + "x5: 2.0->3.0: score=(0.46789042692929883, 0.37050797119720974, 0.2606712423256323)\n", + "x5: 2.0->4.0: score=(0.5616797256969114, 0.5399476525078897, 0.3798806943754639)\n", + "x5: 3.0->4.0: score=(0.17625937121629423, 0.2691689069247335, 0.11920945204983158)\n" + ] + } + ], + "source": [ + "score_dic = {}\n", + "for name in feature_names:\n", + " val = {}\n", + " for v in df[name].unique():\n", + " df_temp = df[df[name] == v]\n", + " vc = df_temp[\"y\"].value_counts().reindex([0, 1], fill_value=0)\n", + " val[v] = vc[1] * 1.0 / (vc[0] + vc[1])\n", + " sorted_val = sorted(val.items(), key=operator.itemgetter(1))\n", + "\n", + " for i, (before, _) in enumerate(sorted_val):\n", + " for after, _ in sorted_val[i + 1 :]:\n", + " lewis = lingam.LEWIS()\n", + " score = lewis.get_scores(\n", + " X_test,\n", + " x_names=[name],\n", + " x_values=[after],\n", + " x_prime_values=[before],\n", + " o_name=\"y\", \n", + " c_names=backdoor[name])\n", + " score_dic = max_nesuf(score_dic, score, name)\n", + " print(f\"{name}: {before}->{after}: score={score}\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "40967ab8-5750-4545-b5d8-eada914241e4", + "metadata": {}, + "source": [ + "## Visualization of LEWIS Scores Sorted by NeSuf\n", + "Finally, we visualize the LEWIS scores using a horizontal bar chart, where features are **sorted in descending order of NeSuf**. NeSuf represents the strongest indicator of causal importance, capturing how often a feature change would flip the decision outcome in both directions. By sorting features by NeSuf, the plot highlights which features are most causally influential overall. For each feature, the corresponding **Nec and Suf values for the same maximizing contrast** are displayed alongside NeSuf, enabling a consistent and interpretable comparison." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "8a6d2678-38f0-43eb-b278-6441d39cf6e6", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_scores_barh_featurecol(\n", + " df,\n", + " feature_col=\"feature\",\n", + " value_cols=(\"Nec\", \"Suf\", \"NeSuf\"),\n", + " title=\"LEWIS Scores\",\n", + " reverse_for_image=True,\n", + "):\n", + " needed = [feature_col, *value_cols]\n", + " missing = [c for c in needed if c not in df.columns]\n", + " if missing:\n", + " raise ValueError(f\"df is missing columns {missing}. Required columns: {needed}\")\n", + "\n", + " df_plot = (\n", + " df.loc[:, needed]\n", + " .dropna(subset=[feature_col])\n", + " .drop_duplicates(subset=[feature_col])\n", + " .set_index(feature_col)\n", + " .loc[:, list(value_cols)]\n", + " )\n", + "\n", + " if reverse_for_image:\n", + " df_plot = df_plot.iloc[::-1]\n", + "\n", + " y = np.arange(len(df_plot.index))\n", + " h = 0.22 \n", + "\n", + " fig, ax = plt.subplots(figsize=(3.0, 3.6), dpi=120)\n", + "\n", + " bars1 = ax.barh(\n", + " y + h,\n", + " df_plot[value_cols[0]].values,\n", + " height=h,\n", + " color=\"#4C72B0\",\n", + " edgecolor=\"black\",\n", + " linewidth=0.6,\n", + " label=value_cols[0],\n", + " )\n", + " bars2 = ax.barh(\n", + " y,\n", + " df_plot[value_cols[1]].values,\n", + " height=h,\n", + " color=\"#F0E442\",\n", + " edgecolor=\"black\",\n", + " linewidth=0.6,\n", + " label=value_cols[1],\n", + " )\n", + " bars3 = ax.barh(\n", + " y - h,\n", + " df_plot[value_cols[2]].values,\n", + " height=h,\n", + " color=\"#2CA02C\",\n", + " edgecolor=\"black\",\n", + " linewidth=0.6,\n", + " label=value_cols[2],\n", + " )\n", + "\n", + " ax.set_yticks(y)\n", + " ax.set_yticklabels(df_plot.index)\n", + " ax.set_xlim(0, 1.0)\n", + " ax.set_xticks([0.0, 0.5, 1.0])\n", + " ax.xaxis.tick_top()\n", + " ax.xaxis.set_label_position(\"top\")\n", + " ax.set_xlabel(title)\n", + " ax.legend(loc=\"lower right\", fontsize=7, frameon=True)\n", + "\n", + " def add_labels(bars):\n", + " for b in bars:\n", + " w = b.get_width()\n", + " yy = b.get_y() + b.get_height() / 2\n", + " ax.text(w + 0.01, yy, f\"{w:.2f}\", va=\"center\", ha=\"left\", fontsize=7)\n", + "\n", + " add_labels(bars1)\n", + " add_labels(bars2)\n", + " add_labels(bars3)\n", + "\n", + " ax.grid(False)\n", + " plt.tight_layout()\n", + " plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "d3bbef13-c1b6-49c3-ba3c-b42ad9c8a182", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVsAAAGjCAYAAACCDraZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAABJ0AAASdAHeZh94AABNJUlEQVR4nO3deVxU5ds/8M8Asg67IDsugChgaiKKJpnljyVZFEnFBcwsQS01Lb9qgKJPpuFWVFaI8jX3lMRd1NI0NVFxsBJEAlFhFFlGB2S5f3/4MI/jgMLAnFm43q/XvIL73OfMdWy8PHOf+74OjzHGQAghRKG0lB0AIYR0BJRsCSGEA5RsCSGEA5RsCSGEA5RsCSGEA5RsCSGEA5RsCSGEA5RsCSGEA5RsCSGEA5RsCSGEA5RsCSGEA5RsCSGEA5RsO4DU1FTweDz8+eefzfYpKCgAj8dr9vX5558DAAIDA2Fubo7n6xddvnwZPB4Pzs7OMsc+ceIEeDweNm7cCAA4deoUeDwedu/eLdXv2rVrCA8Ph7OzM/T19WFvb4+33noLGzZsaNF57t+/H35+frC2toahoSG6d++OiIgIHD58uEX7E6JIOsoOgKiW8ePHIzAwUKa9X79+AIChQ4fi0KFDEAgE8PLykmz//fffoaOjg8LCQty+fRsODg5S2xr3bc7Zs2cxfPhwODk54b333oONjQ2Kiorwxx9/YN26dZg1a9YL4169ejXmz58PPz8/LFy4EIaGhsjLy8Px48exfft2+Pv7t+rPgZD2RsmWSOnfvz8mTpzY7PbGhHnmzBmZZBsYGIgTJ07gzJkzGDdunGTbmTNnYGlpiV69ejV73OXLl8PU1BQXL16EmZmZ1LbS0tIXxlxXV4dly5bhrbfewtGjR2W2v2z/9tTQ0IAnT55AX1+fs/ck6oGGEUirDBw4ELq6upKr1Ua///47hg0bhoEDB0pta2howB9//AFfX1/weLxmj3vz5k14eHjIJFoAsLa2fmFM9+/fR2VlJYYMGdLk9uf3r66uRnx8PNzc3KCvrw9bW1uMHj0aN2/elPR59OgR5s2bB0dHR+jp6aFnz55YvXq1zPAJj8fDzJkzsXXrVnh4eEBPT08ybFFcXIypU6eiS5cu0NPTg4eHB1JSUmTi27BhAzw8PGBoaAhzc3MMGDAAP/300wvPmagfurIlUh4/foz79+/LtJuZmUFHRwf6+vp49dVXcebMGcm2oqIiFBUVwdfXF+Xl5Thw4IBk27Vr11BZWfnCIQQAcHZ2xrlz5yAQCODp6dmqmK2trWFgYID9+/dj1qxZsLCwaLZvfX093n77bWRmZmLcuHH48MMPUVVVhWPHjkEgEKBHjx5gjCE4OBgnT57Eu+++i759++LIkSOYP38+iouLsWbNGqljnjhxAjt37sTMmTPRuXNndO3aFSUlJRg0aJAkGVtZWeHQoUN49913UVlZiY8++ggA8P3332P27NkIDw/Hhx9+iOrqamRnZ+P8+fOYMGFCq/4ciIpjRONt2rSJAWAXL15sts+tW7cYgGZf586dk/SdP38+A8Bu377NGGNs27ZtTF9fn9XU1LCDBw8ybW1tVllZyRhj7KuvvmIA2O+//y7Z/+TJkwwA27Vrl6Tt6NGjTFtbm2lra7PBgwezBQsWsCNHjrAnT5606Bw/++wzBoAZGRmxgIAAtnz5cnbp0iWZfikpKQwAS0pKktnW0NDAGGNs3759DABLTEyU2h4eHs54PB7Ly8uTtAFgWlpaLCcnR6rvu+++y2xtbdn9+/el2seNG8dMTU3Z48ePGWOMhYSEMA8PjxadI1FvlGw7gNYk2+nTp7Njx47JvCoqKiR909PTGQC2bds2xhhjM2fOZEOGDGGMMfbw4UPG4/HY0aNHGWOMjR8/XpKIGzWVbBlj7MKFCywsLIwZGhpKkryVlRVLT09v0Xn+9NNPbOjQoUxLS0uyf79+/dj169clfYKCgljnzp1ZbW1ts8eZPn261D8Yjc6dO8cAsA0bNkjaALDhw4dL9WtoaGBmZmZs+vTpTCgUSr0a/1+cOXOGMcbYlClTmKmpKbtw4UKLzpGoL0q2HUBrku2qVateerz79+8zHo/HZs6cyRhjrF+/fmzBggWS7R4eHiwuLo4xxpiTkxN77bXXpPZvLtk2qqmpYRcuXGALFy5k+vr6rFOnTjJXji9SUVHBjh49yiZMmMAAsB49ejCxWMwYY8zd3V3yD0Nz/t//+3/M0dFRpr28vJwBYB9//LGkDQCbOnWqVL+SkpIXfksAwH7++WfGGGPXr19n9vb2DABzcXFhMTExkkRMNAvdICOtZmlpCXd3d5w5cwYikQjZ2dnw9fWVbPf19cWZM2dw+/ZtFBYWvnS89nm6urrw9vbGihUr8M0336C2tha7du1q8f4mJiZ46623sHXrVkyZMgU3b97E+fPnWxVDaxgYGEj93tDQAACYOHEijh071uSr8WZer1698M8//2D79u0YOnQo9uzZg6FDhyIuLk5h8RLloBtkRC5Dhw5FSkoKjh49ivr6eplku23bNpw6dUrSV14DBgwAANy9e1fu/Tdv3izZv0ePHjh//jxqa2vRqVOnJvdxdnbG8ePHUVVVBWNjY0n733//Ldn+IlZWVjA2NkZ9fT3efPPNl8ZoZGSEd955B++88w6ePHmC0aNHY/ny5Vi4cCFNIdMgdGVL5DJ06FDU19dj9erVcHV1hZWVlWSbr68vRCIRkpOToaWlJZWIm3Py5EmZaVUAcPDgQQBAz549m9338ePHOHfuXJPbDh06JLX/mDFjcP/+fXz11VcyfRvfPzAwEPX19TJ91qxZAx6Ph4CAgBeei7a2NsaMGYM9e/ZAIBDIbBcKhZKfHzx4ILVNV1cXvXv3BmMMtbW1L3wfol7oyrYDSUlJaXLp6ocffij5OSsrC//9739l+vTo0QODBw+W/N54tXru3DlERUVJ9XVzc0Pnzp1x7tw5eHl5NTl39nmzZs3C48ePERYWBnd3dzx58gRnz57Fjh070LVrV0RHRze77+PHj+Hr64tBgwbB398fjo6OKC8vx759+3D69GmEhoZKVsBNnjwZW7Zswdy5c3HhwgW89tprePToEY4fP46YmBiEhIRg1KhRGD58OBYtWoSCggK88sorOHr0KNLT0/HRRx+hR48eLz2fzz//HCdPnoSPjw/ee+899O7dG2VlZcjKysLx48dRVlYGABg5ciRsbGwwZMgQdOnSBX/99Re++uorBAUFSV1VEw2g5DFjwoHGG2TNvYqKil469WvKlCkyx7Wzs2MA2MaNG2W2BQcHMwBsxowZMtuaukF26NAhNnXqVObu7s74fD7T1dVlLi4ubNasWaykpOSF51dbW8u+//57FhoaypydnZmenh4zNDRk/fr1Y6tWrZKaCcEYY48fP2aLFi1i3bp1Y506dWI2NjYsPDyc3bx5U9KnqqqKzZkzh9nZ2bFOnToxV1dXtmrVKsn0sEYAWGxsbJNxlZSUsNjYWObo6Ch5nxEjRkj9eX333Xds2LBhzNLSkunp6bEePXqw+fPnS83+IJqBx1gT390IIYS0KxqzJYQQDlCyJYQQDlCyJYQQDlCyJYQQDlCyJYQQDlCyJYQQDlCyJYQQDqhEsq2pqcEnn3wCOzs7GBgYwMfHB8eOHWvRvsXFxYiIiICZmRlMTEwQEhKC/Px8BUdMFEnez0N8fHyTD6uk+gLqSyQSIS4uDv7+/rCwsACPx0NqamqL9y8vL8f06dNhZWUFIyMjDB8+HFlZWYoL+AVUYrluVFQUdu/ejY8++giurq5ITU1FYGAgTp48+cIiJiKRCMOHD0dFRQX+85//oFOnTlizZg38/Pxw5coVWFpacngWpL3I+3lo9M0334DP50t+19bWVmS4RIHu37+PpUuXwsnJCa+88oqkuFFLNDQ0ICgoCFevXsX8+fPRuXNnJCcn4/XXX8elS5fg6uqquMCbouwlbOfPn5epoyoWi1mPHj3Y4MGDX7jvypUrGQCpwst//fUX09bWZgsXLlRYzERx2vJ5iIuLYwCYUChUdJiEI9XV1ezu3buMMcYuXrzIALBNmza1aN8dO3bILAsvLS1lZmZmbPz48YoI94WUPoywe/duaGtrY/r06ZI2fX19vPvuuzh37hyKiopeuK+3tze8vb0lbe7u7hgxYgR27typ0LiJYrTl89CIMYbKysomq4gR9aKnpwcbGxu59t29eze6dOmC0aNHS9qsrKwQERGB9PR01NTUtFeYLaL0ZHv58mW4ubnBxMREqn3gwIEAgCtXrjS5X0NDA7KzsyX1Tp/f9+bNm6iqqmr3eIliyft5eFb37t1hamoKY2NjTJw4ESUlJYoIlai4y5cvo3///tDSkk5zAwcOxOPHj3Hjxg1O41F6sr179y5sbW1l2hvb7ty50+R+ZWVlqKmpkWtforrk/TwAgLm5OWbOnInvvvsOu3fvxrRp07Bjxw689tprqKysVFjMRDW15bOkCEq/QSYWi6GnpyfT3ngHWSwWN7sfALn2JapL3s8DIF2XF3haKHzgwIGIjIxEcnIyPv300/YNlqi0tnyWFEHpV7YGBgZNjp1UV1dLtje3HwC59iWqS97PQ3MmTJgAGxsbHD9+vF3iI+qjvT9LbaX0ZGtra9vk86Ua2+zs7Jrcz8LCAnp6enLtS1SXvJ+HF3F0dJQ8GYF0HIr4LLWF0pNt3759cePGDZkxtcanofbt27fJ/bS0tODl5YU///xTZtv58+fRvXt3eqyIGpL389AcxhgKCgqknpFGOoa+ffsiKytL8rTjRufPn4ehoSHc3Nw4jUfpyTY8PBz19fXYuHGjpK2mpgabNm2Cj48PHB0dAQCFhYWSp5s+u+/FixelEu4///yDEydOYOzYsdycAGlXbfk8PPsgxUbffPMNhEIh/P39FRs4Uaq7d+/i77//lnpIZnh4OEpKSvDzzz9L2u7fv49du3Zh1KhRTY7nKhTnM3ubMHbsWKajo8Pmz5/PvvvuO+br68t0dHTYr7/+Kunj5+fHng+3srKS9ejRg1lbW7MvvviCrVmzhjk6OjI7OztWWlrK9WmQdiLv58HAwIBFRUWxL7/8kn399dds/PjxjMfjsb59+7JHjx5xfRqknWzYsIEtW7aMzZgxgwFgo0ePZsuWLWPLli1j5eXljDHGpkyZwgCwW7duSfarq6tjgwYNYnw+nyUkJLCvv/6aeXh4MGNjY/b3339zfh4qkWzFYjH7+OOPmY2NDdPT02Pe3t7s8OHDUn2a+svFGGNFRUUsPDycmZiYMD6fz95++22Wm5vLVehEAeT9PEybNo317t2bGRsbs06dOjEXFxf2ySefsMrKSi7DJ+3M2dm52QeRNibXppItY4yVlZWxd999l1laWjJDQ0Pm5+fHLl68yP1JMHrgIyGEcELpY7aEENIRULIlhBAOULIlhBAOULIlhBAOULIlhBAOULIlhBAOULIlhBAOqFWyLS0tRXx8PEpLS5UdClEB9HkgjdThs6BWyVYoFCIhIaHJNfCk46HPA2mkDp8FtUq2hBCirijZEkIIByjZEkIIB5T+DLLWaCwoff36dSVHQlRBXl6e1H9Jx9WYE1T5wZ5qlWwbHz0cERGh5EiIKgkNDVV2CERF3LhxA4MHD1Z2GE1Sq2Tr5eUFANi3bx9cXFyUHA0hRFXk5eUhNDRUkiNUkVol28bHWLi4uMDDw0PJ0RBCVA3nj7ppBbpBRgghHKBkSwghHKBkSwghHKBkSwghHKBkSwghHKBkSwghHKBkSwghHKBkSwghHKBkSwghHKBkSwghHKBkSwghHKBkSwghHKBkSwghHFCrql+NPD09lR0CIYS0ilomW7/J62Hc2UnZYRBCVETV/UL8umW2ssN4IRpGIIRopG+++Qb9+/dHp06dEB8f32y/hoYGfPTRRzAzM0OXLl2wZs0aqe2HDh2Ci4sLjIyMEBISgocPH8oVDyVbQohGsrW1RXx8PMaMGfPCft9++y1OnTqFGzdu4MyZM1i9ejUyMzMBAKWlpRg/fjzWr18PoVAIMzMzzJ4t3xU0JVtCiEYKDQ1FcHAwzMzMXtgvLS0NH3/8MaytreHq6or33nsPW7ZsAQDs3bsXAwYMQGBgIAwNDREfH49du3ZBLBa3Oh5KtoSQDu369evo06eP5HcvLy/k5OQ0ua1bt27o1KkTbt682er3oWRLCOnQRCIRTExMJL+bmJhAJBI1ue357a1ByZYQ0qHx+XxUVlZKfq+srASfz29y2/PbW4OSLSGkQ+vduzeuXbsm+V0gEEie3v38toKCAtTW1qJHjx6tfh9KtoQQjVRXV4fq6mrU19dL/fy8iRMnYvXq1RAKhcjLy8P333+PyZMnAwDCwsJw8eJFHD58GI8fP0ZCQgLGjh0LAwODVsdDyZYQopESExNhYGCAH374AcuXL4eBgQHS0tJw+vRpqWGAGTNmwM/PD66urvD19cXcuXMxYsQIAIC1tTV++uknxMbGonPnznjw4AHWr18vVzw8xhhrlzPjQE5ODjw9PWkFGSFESuMKsmeHAFSNWi7XVfVleYS0lp4eDzU1anPdQ+SgsGR79+5drFu3DufPn8eff/4JkUiEkydP4vXXX2/zsc9lDkSvnq2/G0iIqjJzOAE1+pKpchq/9aoyhY3Z/vPPP1i5ciWKi4vh5eWlqLchRCMJhUIEBQXByMgIPXv2lCwfbUpqaipcXV3B5/PRq1cvyYR7xhji4+Ph6OgIMzMzvPfee3jy5AlXp0Ceo7Bk++qrr+LBgwe4ceMG5s6dq6i3IUQjxcbGwsbGBkKhEKtWrUJERATKyspk+h04cABr1qxBeno6qqqqsH//flhYWAB4moR37dqFP/74A0VFRbh37x6WLl3K9amQ/9XqZCsWi+Hu7g53d3ep9cFlZWWwtbWFr68v6uvrYWxsLPmfTghpOZFIhH379iEhIQGGhoYIDg6Gl5cX0tPTZfouXboUSUlJ6N27N3g8HlxcXGBubg7gaSJ+//33YW9vD2NjY3z66adITU3l+GxIo1YnWwMDA2zevBl5eXlYtGiRpD02NhYVFRVITU2FtrZ2mwMrLS1FTk6O1CsvL6/NxyVE1eXm5oLP58PBwUHS9ux6/Ub19fXIysqCQCCAo6MjunfvjsTERKmx3+d/Li4uRkVFheJPgsiQ6waZj48PFixYgJUrVyIsLAwlJSXYvn071q5dCzc3t3YJLDk5GQkJCe1yLELUSXPr8R88eCDVVlJSgrq6Ohw9ehTXrl1DeXk5Ro4cCWdnZ0yaNAn+/v748ssvERoaClNTU/zP//wPAODRo0cwNTXl7HzIU3LPRoiPj0dGRgamTJkCkUgEPz8/ues8NiUmJgZjx46VasvLy0NoaGi7vQchqqil6/EbVzEtWLAAZmZmMDMzw/vvv4+DBw9i0qRJmDp1KoqKiuDn54e6ujrMmzcPx44dQ5cuXTg7F/J/5L5Bpquri5SUFNy6dQtVVVXYtGkTeDxeuwVmbW0NDw8PqZeLi0u7HZ8QVeXq6gqRSITi4mJJW1OT9c3NzWFnZyf19+7Zn7W0tJCQkICCggLcvn0bHh4e6N+/f7sM85HWa9NshCNHjgAAqqurkZub2y4BEdLR8fl8hISEIC4uDmKxGBkZGcjOzkZISIhM36ioKHzxxReoqqrC7du3sXHjRgQFBQEA7t+/j/z8fDDGkJOTg3nz5iEuLo7r0yH/S+5km52djaVLlyI6Ohr9+vXDtGnTaOCdkHaSnJyMO3fuwNLSEnPnzsWOHTtgYWGBrVu3Sl3hxsXFwdbWFg4ODhg0aBAmTJiAiRMnAng6V3fkyJGSZ2fNmzcPAQEByjqlDk+u2gi1tbXw8fHBw4cPkZ2djVu3bsHb2xuRkZFISUmR6b97926MHTu2zSvIGleJ0AoyomloBVnbNOYGjauNkJiYiCtXriAzMxPGxsbo06cPPvvsMyxevBjh4eEIDAyU9AMgmbKSlpaGM2fOAAAWL14sd9CDR1yQe19CVJGOHo2jarpWX9lmZWXBx8cHM2bMkCo1Vl9fj8GDB6O4uBg5OTkwMzN74Q0zef4Vb/zXy2W5C/Tt9Vu9PyGqShAloCvbNlCHK9tWj9n2798ftbW1MjUdtbW1ceHCBRQXF0ueZskYa/ZFCGke1UbQPFQ8nBAVRLURNA8lW0JUDNVG0EyUbAlRMVQbQTNRsiVExTRXG0EkEkm1PV8b4cSJE9iyZQv++9//AgD8/f3x7bff4t9//0V5eblUbQTCPUq2hKgYeWsjdO3aVVIbAQCmTp2KiIgI+Pn5wdPTE2+++SY6depEtRGUhJItISqGaiNoJkq2hKgYqo2gmSjZEqKCqDaC5pGrNoKy0AoyoqloBVnbaOQKMkIIIa0n95MalClvET2LjGgWWwdbZYdAFEwtk60qf1UghJCm0DACIYRwgJItIYRwgJItIYRwgJItIYRwgJItIYRwgJItIYRwgJItIYRwgJItIYRwgJItIYRwgJItIYRwgJItIYRwgJItIYRwgJItIYRwgJItIYRwgJItIYRwgJItIYRwgJItIYRwQC2f1ODp6ansEAjpsOzsHVF8u1DZYagdtUy2fpPXw7izk7LDIKRDykgKlfpdKBQiKioKp06dgoODA5KTkzFixAiZ/QoKCvD+++/j/PnzMDIywowZM7B48WIAQGpqKqZNmwZ9/f97avb169fh5KQ5f89pGIEQ0iaxsbGwsbGBUCjEqlWrEBERgbKyMpl+s2bNgpOTE4RCIc6cOYPk5GQcOXJEsv3111+HSCSSvDQp0QKUbAkhbSASibBv3z4kJCTA0NAQwcHB8PLyQnp6ukzfgoICREREoFOnTujWrRuGDh2K69evKyFq5aBkSwiRW25uLvh8PhwcHCRtXl5eyMnJkekbGxuLHTt2oKamBrm5ufjjjz8wfPhwyfY//vgDlpaW6N27N7799ltO4ueSWo7ZEkJUg0gkgomJiVSbiYkJHjx4INP3tddew7fffgsjIyPU19dj2bJl6Nu3LwDAz88PAoEATk5OuHjxIsLCwmBlZYUxY8ZwcRqcoCtbQojc+Hw+KisrpdoqKyvB5/Ol2urr6+Hv74+oqChUV1cjPz8fW7duxS+//AIA6NatG7p27QotLS34+Phg9uzZ+Pnnnzk7Dy5QsiWEyM3V1RUikQjFxcWSNoFAAA8PD6l+ZWVluH37NmbMmAEdHR1069YNQUFByMzMbPK4WlpaYIwpNHauUbIlhMiNz+cjJCQEcXFxEIvFyMjIQHZ2NkJCQqT6WVlZwcnJCd9//z0aGhpQVFSEAwcOwMvLCwBw+PBhCIVCAEBWVhbWr1+P4OBgzs9HkSjZEkLaJDk5GXfu3IGlpSXmzp2LHTt2wMLCAlu3bpW6wt29ezd++uknmJubY+DAgQgMDMTUqVMBAMeOHYOHhweMjIwwbtw4fPLJJxg3bpyyTkkheEyNrtVzcnLg6elJixoIUaKMpFCV+4rfmBuaGsJQFXRlSwghHFDLqV+/bpmt7BAIgZ4eDzU1qnWFxwU7e0dlh6CWOE227733Hn744QcEBQUhIyND7uOcyxyIXj35L+9IiAKZOZyQ+jrd0hoBwNNaAMuXL8fdu3fh6OiIjIwM9OjRAwcOHMDy5cuRk5MjGb9cuXIlOnXqxNVpEQXhbBjhzz//RGpqqlShCUI0SUtrBBw4cABr1qxBeno6qqqqsH//flhYWAB4Okc1Pj4e9+7dw9WrV3Hx4kWsWrWK61MhCsDJlS1jDLNnz8bkyZObnVdHiDprrBGQn58vUyMgOjpaqu/SpUuRlJSE3r17AwBcXFwk28aPHy/52cDAAJMmTcL+/fu5OQmiUHJd2YrFYri7u8Pd3R1isVjSXlZWBltbW/j6+qK+vl7SnpaWBoFAgOXLl7c9YkJUUEtrBNTX1yMrKwsCgQCOjo7o3r07EhMTm727/9tvv6ns3XXSOnJd2RoYGGDz5s0YMmQIFi1ahKSkJABPv0ZVVFQgNTUV2traAICqqip88skn+M9//gMbG5sWv0dpaalkknOjvLw8ecIlROFaWiOgpKQEdXV1OHr0KK5du4by8nKMHDkSzs7OmDRpklTfPXv2IDMzE1evXlV4/ETx5B5G8PHxwYIFC7By5UqEhYWhpKQE27dvx9q1a+Hm5ibpt3TpUhgYGGDOnDmtOn5ycjISEhLkDY8QTrW0RoCBgQEAYMGCBTAzM4OZmRnef/99HDx4UCrZnjx5EjNmzMDBgwdhbW2t+BMgCtemMdv4+HhkZGRgypQpEIlE8PPzw+zZ/zct68aNG1i3bh22bdsGPT29Vh07JiYGY8eOlWrLy8tDaGhoW0ImRCGerRFgb28P4GmNgMmTJ0v1Mzc3h52dHXg8nqTt2Z8B4Pz584iIiMCuXbswYMAAxQdPONGmZKurq4uUlBR4e3tDX18fmzZtkvrgfPjhh/D19ZWrTJq1tTX9i07UxrM1AjZs2IDMzMwmawQAQFRUFL744gv069cPFRUV2Lhxo+TxMNeuXcOoUaPw448/4vXXX+f4LIgitXnqV+NjLaqrq5GbmytpP3HiBA4fPowPP/wQBQUFklddXR3EYjEKCgpkvnYRos5aWiMgLi4Otra2cHBwwKBBgzBhwgRMnDgRAJCUlIQHDx5gwoQJ4PP54PP5CAgIUNYpkXbUptoI2dnZ8Pb2RmRkJK5cuYL79+/j2rVrMDU1RWpqqsyUl+etWbMGH330UYvfr3H9My1qIKrg+UUNRHnUoTaC3MMItbW1iIqKgp2dHdatW4dbt27B29sbc+bMQUpKCt544w3s3btXZr/p06fD2dkZixYtkpRXI4QQjcfk9NlnnzEej8dOnDghaUtMTGQA2IEDB5rdz9nZmQUFBcn1ngKBgAGgF71U5mXrYCvXZ5m0r8bcIBAIlB1Ks+S6ss3KysKKFSswc+ZMqQe2ffrpp0hPT8d7772HnJwcmJmZyXP4l3JZ7gJ9e1r2S5RPECVQdghETch1g6x///6ora3F+vXrpdq1tbVx4cIFFBcXN5toCwoK2lSEhhBVJxQKERQUBCMjI/Ts2fOFS9RTU1Ph6uoKPp+PXr164ebNmwCe1k/w9fWFqakp7OzsMHfuXNTW1nJ1CkQBqJ4tIe2MCtKQpqhlPVtCVBUVpCHNoStbQtoRFaQhzaErW0LaERWkIc2hK1tC2pG8BWm6du0qKUjzrMaCNPv376fl62qOki0h7ejZgjSNmlrV1JqCNDt37qSCNBqAki0h7ejZgjRisRgZGRkvLUhTVVWF27dvY+PGjQgKCgJABWk0ESVbQtoZFaQhTWlTIRquNRaboBVkRFUIogRUjEYFaHQhGmXKW0SPxyGqwdbBVtkhEDWhlslWlf/1IoSQptCYLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcEAtn9Tg6emp7BAIUUl29o4ovl2o7DBIE9Qy2fpNXg/jzk7KDoMQlZORFKrsEEgzaBiBEA0mFAoRFBQEIyMj9OzZE5mZmU328/DwkDwync/nQ1tbG7NmzZJsv3fvHsLDw2FqagpLS0ssXLiQq1PQGGp5ZUsIaZnY2FjY2NhAKBTi+PHjiIiIQG5uLiwsLKT65eTkSH6uqamBjY0NxowZI2kLDg7GuHHjsGXLFvB4POTl0ROuW4uSLSEaSiQSYd++fcjPz4ehoSGCg4Ph5eWF9PR0REdHN7vf/v37YWJiAj8/PwDAoUOHoKenh7lz50r6eHl5KTx+TUPDCIRoqNzcXPD5fDg4OEjavLy8pK5im5KWloaJEyeCx+MBAC5cuICuXbsiICAAnTt3xogRI/DXX38pNHZNRMmWEA0lEolgYmIi1WZiYgKRSNTsPg8ePMChQ4cwadIkSVtxcTG2b9+O2bNn486dO/D390doaCjq6+sVFrsmomRLiIbi8/morKyUaqusrASfz292n+3bt6Nv375wd3eXtBkYGGDo0KEICAiArq4uPv74Y9y7d4/GbVuJki0hGsrV1RUikQjFxcWSNoFAAA8Pj2b3SUtLk7qqBZ7Oa28cUiDyo2RLiIbi8/kICQlBXFwcxGIxMjIykJ2djZCQkCb75+bmIisrC+PHj5dqHz16NK5evYrjx4+jvr4ea9euhY2NDVxcXLg4DY1ByZYQDZacnIw7d+7A0tISc+fOxY4dO2BhYYGtW7fKXOGmpaXB398fnTt3lmq3tLTE7t27MWvWLJiZmWHv3r3Yu3cvtLW1uTwVtcdjjDFlB9FSOTk58PT0pBVkhDQjIykUavRXut005oaXDZMok1rOs/11y2xlh0A6MD09HmpqVDOh2dk7KjsE0gyFJtvU1NRmJ0/fvXsXNjY2ch33XOZA9OrZ/B1VQhTJzOFEh7x6JG3DyZjt0qVLkZaWJvUyMzPj4q0J4URLaxAATy9CXF1dwefz0atXL9y8eRMAUFJSglGjRsHa2pru/msgToYRAgICMGDAAC7eihClaGkNggMHDmDNmjVIT0+XJNrGPlpaWggMDERsbCwCAgKUcRpEgeS6shWLxXB3d4e7uzvEYrGkvaysDLa2tvD19ZVZXVJVVUUrTohGaqxBkJCQIFOD4HlLly5FUlISevfuDR6PBxcXF5ibmwMArKysMGPGDPTt25fjMyBckCvZGhgYYPPmzcjLy8OiRYsk7bGxsaioqEBqaqrUtJDhw4fDxMRE8kHMzc196XuUlpYiJydH6kUrVogqamkNgvr6emRlZUEgEMDR0RHdu3dHYmIijf92EHIPI/j4+GDBggVYuXIlwsLCUFJSgu3bt2Pt2rVwc3MDABgaGiIqKkqSbC9duoSkpCT4+voiKysLjo7N3zlNTk5GQkKCvOERwpnmahA8ePBAqq2kpAR1dXU4evQorl27hvLycowcORLOzs4yq7aI5mnTPNsnT55gwIABEIlEEIlE6N27N06ePPnCwf0zZ85g2LBhmD59Or799ttm+5WWlkIoFEq15eXlITQ0lGYjEKV6fjbC5cuXMWLECJSVlUnaZs2aBT09PaxevVrS9vDhQ1hYWODUqVOS8oVffvkl/vzzT2zbtk3S7969e7C1taUr3lbQ+Hm2urq6SElJgbe3N/T19bFp06aX3kUdOnQofHx8cPz48Rf2s7a2hrW1dVvCI4QTz9YgsLe3B/C0BsHkyZOl+pmbm8POzk7q7wjNOug42jz168iRIwCA6urqFo3FAoCjo6PUVQAh6qw1NQiioqLwxRdfoKqqCrdv38bGjRsRFBQk2V5dXY2amhqZn4n6a1Oyzc7OxtKlSxEdHY1+/fph2rRpqKioeOl++fn5sLKyastbE6JSWlqDIC4uDra2tnBwcMCgQYMwYcIETJw4UbLdwMAAXbt2lfzcs2dPrk+FKIjcY7a1tbXw8fHBw4cPkZ2djVu3bsHb2xuRkZFISUkB8HSi9/NJ9eDBgwgKCsLs2bOxbt26Vr1n47gMjdkSZaIVZKpHo8dsExMTceXKFWRmZsLY2Bh9+vTBZ599hsWLFyM8PByBgYHw9fVFv379MGDAAJiamiIrKwspKSlwdHTEf/7zn/Y8D0IIUW1MDpcuXWI6Ojps1qxZUu11dXXM29ub2dnZsYcPH7JFixaxvn37MlNTU9apUyfm5OTEZsyYwe7duyfP2zKBQMAA0IteTb54OjxO3sfZyVauzy9RnMbcIBAIlB1Ks9SyxKLLchfo2+srOxyiYgRRApmv90KhEFFRUTh16hQcHByQnJyMESNGyOwbFRWFbdu2oVOnTgAAZ2dnyaKEkpISTJs2DefPn4dQKKQhBBWkDsMIVDycaLRnaxasWrUKERERzc6EWbJkiWTO+LOrvxprFmzZsoWrsIkGUst6toS0RGPNgvz8fJmaBc2V/mxKY82Ce/fuKTBaounoypZorJbWLGi0Zs0aWFpawtfXF7/++itXYZIOgpIt0VjN1SwQiUQyfT/88EPk5eXh7t27iI2NRXBwMP7991+uQiUdACVborH4fD4qKyul2iorK8Hny87R7tevH8zNzaGrq4vIyEgMHjwYR48e5SpU0gFQsiUa69maBY1aerdaS0uLZh2QdkXJlmis1tQs2LNnDx49eoS6ujrs2LEDp0+fxptvvinZTjULSFtRsiUaraU1C9asWQM7OztYWloiKSkJ+/btQ/fu3SXbqWYBaSua+kU0mpWVFQ4ePCjTHhkZicjISMnvZ86ceeFxaEiBtBVd2RJCCAfU8so2bxE9i4zIsnWwVXYIhDRLLZOtKq9/JoSQptAwAiGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcICSLSGEcEAtn9Tg6emp7BAIaXd29o4ovl2o7DCIgqhlsvWbvB7GnZ2UHQYh7SojKVSmTSgUIioqCqdOnYKDgwOSk5MxYsQImX4eHh74999/Jb+LxWLExMRgw4YNAIBDhw5h1qxZuHv3Lt58802kpqbC3NxcYedCZNEwAiEqLDY2FjY2NhAKhVi1ahUiIiJQVlYm0y8nJwcikQgikQgPHjyAiYkJxowZAwAoLS3F+PHjsX79egiFQpiZmWH27Nlcn0qHR8mWEBUlEomwb98+JCQkwNDQEMHBwfDy8kJ6evoL99u/fz9MTEzg5+cHANi7dy8GDBiAwMBAGBoaIj4+Hrt27YJYLObiNMj/omRLiIrKzc0Fn8+Hg4ODpM3Lyws5OTkv3C8tLQ0TJ04Ej8cDAFy/fh19+vSRbO/WrRs6deqEmzdvKiZw0iRKtoSoKJFIBBMTE6k2ExMTiESiZvd58OABDh06hEmTJrXpOKT9UbIlREXx+XxUVlZKtVVWVoLP5ze7z/bt29G3b1+4u7u36Tik/VGyJURFubq6QiQSobi4WNImEAjg4eHR7D5paWlSV7UA0Lt3b1y7dk3ye0FBAWpra9GjR4/2D5o0i5ItISqKz+cjJCQEcXFxEIvFyMjIQHZ2NkJCQprsn5ubi6ysLIwfP16qPSwsDBcvXsThw4fx+PFjJCQkYOzYsTAwMODiNMj/omRLiApLTk7GnTt3YGlpiblz52LHjh2wsLDA1q1bZa5w09LS4O/vj86dO0u1W1tb46effkJsbCw6d+6MBw8eYP369VyeBgHAY4wxZQfRUjk5OfD09KRFDUQjZSSFQo3+OqqUxtzwsmEWZaIrW0II4YBaLtf9dQutfumI9PR4qKnR3Cs/O3tHZYdAFEhhyTYzMxNbt27FmTNncPv2bdjY2OCNN97AsmXLYGtr26Zjn8sciF49adpKR2PmcELma3ZLawc0KigoQK9evRAZGYkffvgBAMAYw+LFi5GSkoLq6moMHToU3333Hezs7BR6PqRjUdgwwieffIJTp04hLCwM69evx7hx47Bz507069cP9+7dU9Tbkg6mpbUDGs2ZMwf9+/eXavv555+RlpaG8+fPo6SkBJaWlpg3b56iQycdjMKubJOSkjB06FBoaf1fPvf394efnx+++uorJCYmKuqtSQfRWDsgPz9fpnZAdHS0TP8jR46AMYa33noLt2/flrQXFBTgtddeg5PT05uuERERWLhwIWfnQTqGVl/ZisViuLu7w93dXaqQRVlZGWxtbeHr64v6+noMGzZMKtECwLBhw2BhYYG//vqr7ZGTDq81tQOePHmC+fPn48svv5TZFh4ejhs3buDWrVsQi8XYtm0bRo4cqdDYScfT6itbAwMDbN68GUOGDMGiRYuQlJQE4OnXuYqKCqSmpkJbW7vJfRtLwD0/D7AppaWlEAqFUm15eXmtDZdosObW/D948ECmb1JSEgIDA5tcNWVjY4OBAweie/fu0NbWRp8+fZCcnKywuEnHJNcwgo+PDxYsWICVK1ciLCwMJSUl2L59O9auXQs3N7dm91u7di2ePHmCd95556XvkZycjISEBHnCIx1ES9f8FxcXIyUlBVlZWU0eJyEhAdevX0dpaSn4fD4WLlyIKVOm4Oeff1ZY7KTjkXvMNj4+HhkZGZgyZQpEIhH8/PxeWJD4t99+Q0JCAiIiIvDGG2+89PgxMTEYO3asVFteXh5CQ0PlDZlomGdrB9jb2wN4Wjtg8uTJUv0uXryIoqIiuLi4AHh6RdzQ0ICCggIcP34cV69exbhx42BlZQUAmDZtGoYMGcLtyRCNJ3ey1dXVRUpKCry9vaGvr49NmzZJ6mc+7++//0ZYWBg8PT0l021extraGtbW1vKGRzqAZ2sHbNiwAZmZmU3WDggICMCtW7ckv69evRp3796VLFkdMGAAdu7cibFjx4LP5yMlJQVeXl6cngvRfG2a+nXkyBEAQHV1NXJzc5vsU1RUhJEjR8LU1BQHDx6EsbFxW96SECktqR2gp6cHGxsbyYvP58PAwACWlpYAnk5T7Nq1K3r16oUuXbogOzsbP/74ozJPi2gguWsjZGdnw9vbG5GRkbhy5Qru37+Pa9euwdTUVNLnwYMHGDp0KMrKynDmzBm4urq2KdjG9c+0qKFjampRAyGABtdGqK2tRVRUFOzs7LBu3TqkpqaipKQEc+bMkfR59OgRAgMDUVxcjIMHD7Y50RJCiDqTa8w2MTERV65cQWZmJoyNjdGnTx989tlnWLx4McLDwxEYGIjIyEhcuHABU6dOxV9//SU1t5bP57fpRtfgERfk3pcoHk+HB1bX/legzk5tW+ZNiDK1ehghKysLPj4+mDFjhlRNzPr6egwePBjFxcXIyclB3759pZ5j/yxnZ2cUFBS0OtjGrwouy12gb6/f6v0JNwRRAvq6TzilkcMI/fv3R21trUzxYW1tbVy4cAHFxcUwMzNDQUEBGGNNvuRJtER9CYVCBAUFwcjICD179kRmZuYL+xcUFMDAwADTpk2Tak9NTYWDgwNMTEwQHR2NJ0+eKDJsQtoV1bMlCtcexWKuXbuGOXPmYO/evSgqKkJRURGWLVum6NAJaTeUbIlCNRaLSUhIkCkW05Rni8U866effsKYMWPg7e0NU1NTLF68GFu2bOHiFAhpF5RsiUK1V7GY69evo0+fPlLHKCwshEgkUkzghLQzSrZEoZorFtNUknxRsZjnj9P4MyVboi7U8rE4RH20V7GY54/T+PPzxyFEVVGyJQrVXsVievfujWvXrkn6CwQCODk5UbIlaoOGEYhCPVssRiwWIyMj44XFYq5cuYIrV67ggw8+QFhYGHbs2AEAmDBhAvbs2YNLly6hoqICy5cvl0nYhKgySrZE4dqjWIyXlxeSkpIQHBwMBwcH2NnZYfHixco8LUJaRe5CNMpAK8jUA60gI1xThxVkajlmm7eIHo+jymwdqIYBIc9Ty2Sryv96EUJIU2jMlhBCOEDJlhBCOEDJlhBCOEDJlhBCOEDJlhBCOEDJlhBCOEDJlhBCOEDJlhBCOEDJlhBCOEDJlhBCOEDJlhBCOEDJlhBCOEDJlhBCOEDJlhBCOEDJlhBCOEDJlhBCOEDJlhBCOKCWT2rw9PRUdghEhdnZO6L4dqGywyBEilomW7/J62Hc2UnZYRAVlZEUquwQCJFBwwikQxAKhQgKCoKRkRF69uyJzMzMJvtFRUVBT08PfD4ffD5f6ll3Bw4cgK+vL0xNTWFnZ4e5c+eitraWq1Mgao6SLekQYmNjYWNjA6FQiFWrViEiIgJlZWVN9l2yZAlEIhFEIhFycnIk7ZWVlYiPj8e9e/dw9epVXLx4EatWreLqFIiaU8thBEJaQyQSYd++fcjPz4ehoSGCg4Ph5eWF9PR0REdHt/g448ePl/xsYGCASZMmYf/+/YoImWggurIlGi83Nxd8Ph8ODg6SNi8vL6mr1metWbMGlpaW8PX1xa+//trscX/77TepYQZCXoSSLdF4IpEIJiYmUm0mJiYQiUQyfT/88EPk5eXh7t27iI2NRXBwMP7991+Zfnv27EFmZibmzp2rsLiJZqFkSzQen89HZWWlVFtlZSX4fL5M3379+sHc3By6urqIjIzE4MGDcfToUak+J0+exIwZM7B//35YW1srNHaiOSjZEo3n6uoKkUiE4uJiSZtAIGjREICWlhYYY5Lfz58/j4iICOzcuRMDBgxQSLxEM1GyJRqPz+cjJCQEcXFxEIvFyMjIQHZ2NkJCQmT67tmzB48ePUJdXR127NiB06dP48033wQAXLt2DaNGjcKPP/6I119/neOzIOqOki3pEJKTk3Hnzh1YWlpi7ty52LFjBywsLLB161apK9w1a9bAzs4OlpaWSEpKwr59+9C9e3cAQFJSEh48eIAJEyZI5uEGBAQo65SImuGxZ78jqbicnBx4enrSCjLyQhlJoVCjjzVpB425oaXDQ8qglvNsf90yW9khKISeHg81NZQk2srO3lHZIRAiQ2HJ9rfffsPq1atx+fJlCIVCmJmZoW/fvliyZAmGDBnSpmOfyxyIXj1l7ySrOzOHE3RFRoiGUtiY7Y0bN6ClpYUPPvgAX3/9NT7++GPcu3cPw4YNw+HDhxX1thqlpev5582bhx49esDY2Bh9+vRBRkaGZFtJSQlGjRoFa2tr8Hg8rkInhDxHYVe206ZNw7Rp06TaYmJi0L17d6xduxb+/v6KemuN8ex6/uPHjyMiIgK5ubmwsLCQ6mdsbIxDhw7BxcUFv/76K8LCwnD58mV069YNWlpaCAwMRGxsLN3MIUSJWn1lKxaL4e7uDnd3d4jFYkl7WVkZbG1t4evri/r6+ib3NTQ0hJWVFcrLy+UOuKNoXM+fkJAgs57/efHx8XBzc4OWlhaGDx+O3r17IysrCwBgZWWFGTNmoG/fvhyfASHkWa2+sjUwMMDmzZsxZMgQLFq0CElJSQCeXoVVVFQgNTUV2trakv6VlZV48uQJ7t+/jy1btkAgEOA///nPS9+ntLQUQqFQqi0vL6+14aqt1q7nb/Tw4UMIBAL07t1b0SESQlpBrmEEHx8fLFiwACtXrkRYWBhKSkqwfft2rF27Fm5ublJ9IyIicOTIEQCArq4u3n//fSxZsuSl75GcnIyEhAR5wtMIza3nf/DgQbP7NDQ0IDo6GmPGjEGvXr0UHSIhpBXkHrONj49HRkYGpkyZApFIBD8/P8yeLTsl6/PPP8e8efNQVFSEzZs348mTJ6irq3vp8WNiYjB27Fiptry8PISGhsobslppzXr+RjExMaioqMCOHTsUHR4hpJXkTra6urpISUmBt7c39PX1sWnTpibvdj87Vjhx4kT0798fUVFR2L179wuPb21t3aGLfDy7nt/e3h7A0/X8kydPbrL/ggULcOnSJZw4cQJ6enpchkoIaYE2Tf1qHB6orq5Gbm7uS/vr6uoiODgYP//8s9TNNSKrNev5ExMTkZGRgcOHD8PY2Fhme3V1NWpqamR+JoRwR+5km52djaVLlyI6Ohr9+vXDtGnTUFFR8dL9xGIxGGOoqqqS9607jJau51+yZAlu3rwJZ2dnyZr9rVu3SrYbGBiga9eukp979uzJ9akQ0uHJVRuhtrYWPj4+ePjwIbKzs3Hr1i14e3sjMjISKSkpAJ7OJnh+GKC8vBx9+vQBABQWtv5R043rn2kFGSHkWRpbGyExMRFXrlxBZmamZNXSZ599hsWLFyM8PByBgYEICAiAg4MDfHx8YG1tjcLCQmzatAl37tyhGziEkI6HtdKlS5eYjo4OmzVrllR7XV0d8/b2ZnZ2duzhw4fsq6++YkOHDmWdO3dmOjo6zMrKio0aNYr99ttvrX1LCYFAwADQ65mXrYOt3H+ehGiKxtwgEAiUHUqz1LLEostyF+jb6ys7HJUgiBJIDT0IhUJERUXh1KlTcHBwQHJyMkaMGCGz37x587Bv3z6UlpaiW7duWLFiBd5++23J9hs3bmDmzJk4e/YsjIyMsGTJEsycOZOTcyKktdRhGIGKh2uYZ+sprFq1ChERESgrK5Pp11hPoaKiAuvWrcPEiRNx69YtAE9nLAQEBGDKlCkoKytDbm6u5GkFhBD5qGU9W9K0xnoK+fn5MvUUoqOjpfrGx8dLfn62nkK3bt2wadMm+Pr6IjIyEsDTKXvPr2YjhLQOXdlqkPaqp3DhwgVYWFhg8ODBsLa2RmhoqNTDEgkhrUfJVoM0V09BJBI1u09T9RSKi4uxefNmrF+/HoWFhXBycmp25RohpGVoGEGDtFc9BQMDA4SFhcHb2xsAEBcXBysrK4jFYhgYGCgmeEI0HF3ZapBn6yk0etHd2cZ6Cr/88otUPQVPT0+pOhc8Ho+e8kBIG1Gy1SDtVU9h4sSJ+OWXX3DlyhXU1tZi2bJlGD58OF3VEtIGlGw1THvUU+jVqxe+/vprhIaGwsrKCrm5udi8ebOyTokQjUCLGtTc84saCOmIaFEDIYQQAGo6GyFvUcd5FtnL2DrYKjsEQkgLqGWyVeWvCoQQ0hQaRiCEEA5QsiWEEA6o5TACIZqksrISDx8+RH19vbJDUWna2towNzdX26JIlGwJUaInT57gzp07YIxBV1eXVuo1gzGGx48fo7q6Gvr6+tDV1VV2SK1GyZYQJRIKhWCMwcnJCUZGRsoOR6U9evQIhYWFEAqFsLe3V3Y4rUZjtoQo0ZMnT6Cjo0OJtgWMjIygo6ODJ0+eKDsUuVCyJUSJGGPQ1tZWdhhqQ0tLS21XTFKyJYSoDXUe06ZkSwghHKBkSwghHKBkS4iKsXdwkhRsV8TL3sGpRXF07doVzs7OqK2tlbR98MEHUg8LJS1HU78IUTF3iovw9tx9Cjt+RlJoi/tWVVVh06ZNmD59usLi6SjoypYQ0qw5c+ZgxYoVUle3jXbv3g0PDw9YWFggODgYpaWlkm0nTpzAgAEDYGJiAldXV5w+fZrLsFUSJVtCSLOGDx8OJycnpKamSrVfuHABH330EbZv346SkhK4u7sjJiYGAJCfn4/Q0FDEx8fj4cOHyMzMhK0tlQKlYQRCyAvFxcVh2rRpiIqKkrSlpKQgJiYGXl5eAJ4+ZsnCwgJ1dXXYtm0bRo0ahbfffhsA4OTUsjFiTUdXtoSQFxoxYgTs7e2lnkNXWFiI5cuXw8zMDGZmZnB0dISOjg7u3buH27dvo1u3bkqMWDVRsiWEvFRcXJzU2K29vT2WLVuG8vJyyUssFsPBwQGOjo4oKChQbsAqSC2Traenp0KnxtCr/aYPEc3w1ltvwcbGBvv27QMAREdH46uvvsLVq1cBAGVlZUhPTwcAjB8/Hvv378fBgwfR0NCAoqIi3Lx5U1mhqwy1HLP1m7wexp3pL7syPT99SCgUIioqCqdOnYKDgwOSk5MxYsQImf3i4uKQkpKCiooKdOnSBQsXLsTUqVMBACtWrMCKFSskfevq6tCpUydUVVUp9FxUjZ29Y6umZ8lzfHnExcXB398fAODr64vVq1dj8uTJuHXrFiwsLBAREYGQkBB069YNe/bswfz58/HOO+/A1tYWKSkp6NGjR3uehtpRy0eZU7JVvoykUKmCIBERETA2NsaGDRtw/PhxREdHIzc3FxYWFlL75ebmws7ODkZGRrhx4wb8/Pxw9OhRyY2WZ82YMQNisVjmTrgmyc/PBwB0795dyZGoh+b+vOhR5qRDEIlE2LdvHxISEmBoaIjg4GB4eXlJvlY+y9XVVVJOsLGoyK1bt2T6PXnyBDt37sSkSZMUGzwhHKFkS9osNzcXfD4fDg4OkjYvLy/k5OQ02f/zzz+HkZER3NzcYG9vjzfffFOmz4EDB2BoaIjhw4crLG5CuETJlrSZSCSSeS6UiYkJRCJRk/0//fRTiEQi/PHHHxgzZkyTjzhJS0tDZGQktLToI0o0A32SSZvx+XxUVlZKtVVWVoLP5ze7D4/Hg4+PD+7cuYONGzdKbSsrK8OBAwcwefJkhcRLiDJQsiVt5urqCpFIhOLiYklbS29U1NXVIS8vT6pt586d8PT0RO/evds9VkKUhZItaTM+n4+QkBDExcVBLBYjIyMD2dnZCAkJken7/fffo7y8HA0NDTh58iS2bt2KN954Q6pPWloaXdUSjUPJlrSL5ORk3LlzB5aWlpg7dy527NgBCwsLbN26VeoKNyMjAz169ICpqSlmzpyJ1atXS9bQA0+n9ly8eBHjx49XxmkQojBquaiBqB4rKyscPHhQpj0yMhKRkZGS35uaDvas7t27q+3TUwl5EbqyJYS0u927d8Pe3h58Pl+qzm1HppZXtr9uma3sEDijp8dDTY3qLfKTd8knebmuznb4t/Cuwo7v7GSLgn/vtKjvb7/9hgULFuCvv/6Cjo4OXnnlFfz4448vreq1YMECbNq0CSNHjmyPkDWCwpNteXk5FixYgL179+Lx48cYOHAgvvzyS/Tv31/uY57LHIhePZufVqRJzBxOSC2LbWkNgnnz5mHfvn0oLS1Ft27dsGLFCsnY6B9//IHp06ejsLAQurq6CAgIwNdff/3CqVqEO/8W3kX57Tde3lFOZg4nWtSvoqICISEh+PHHHxEaGorHjx/j2LFj0NbWfum+hYWFNJvkOQodRmhoaEBQUBB++uknzJw5E1988QVKS0vx+uuvIzc3V5FvrbFiY2NhY2MDoVCIVatWISIiAmVlZTL9jI2NcejQIVRUVGDdunWYOHGiZFmsi4sLDh06hPLychQUFKChoQEJCQlcnwpRcTdu3ICenh5Gjx4NLS0t8Pl8hIWFwcnJCVFRUUhMTJT0TU1NlawE5PP5qK+vR8+ePTFw4EBlha9yFJpsd+/ejbNnzyI1NRVxcXGIjY3FqVOnoK2tjbi4OEW+tUZqTQ2C+Ph4uLm5QUtLC8OHD0fv3r2RlZUFAOjcuTPs7e0BAIwxaGlpUQk8IsPNzQ1PnjzBtGnTcOzYMZmFK81pXDn4zz//4MKFC4oMUa3IlWzFYjHc3d3h7u4OsVgsaS8rK4OtrS18fX1RX1+P3bt3o0uXLhg9erSkj5WVFSIiIpCeno6ampq2n0EH0toaBI0ePnwIgUAg9bWusLAQZmZm4PP52LNnD2bOnKmwuIl6MjU1xW+//YaamhpMmjQJVlZWmDhxYocredle5Eq2BgYG2Lx5M/Ly8rBo0SJJe2xsLCoqKpCamgptbW1cvnwZ/fv3l1nfPnDgQDx+/Bg3btxo9j1KS0uRk5Mj9Xp+pVFH09oaBMDToZzo6GiMGTMGvXr1krQ7OTmhvLwcJSUlWLhwIT0nijTJ09MTaWlpuHfvHs6ePYuzZ89i+fLlyg5LLcl9g8zHxwcLFizAypUrERYWhpKSEmzfvh1r166Fm5sbAODu3bsYNmyYzL6NT9q8c+dOk3VMgaeT5GkcUZo8NQhiYmJQUVGBHTt2NLnd2toa/v7+mDBhAn3lIy/06quvYvTo0RAIBHB2dpb6VltSUqLEyNRDm8Zs4+Pj4eHhgSlTpiAmJgZ+fn6YPfv/pmWJxWLo6enJ7Kevry/Z3pyYmBgIBAKpV+MjOTqq1tYgWLBgAS5duoRffvmlyf8PjZqqT0DI33//jTVr1uDOnafTxG7cuIH9+/dj4MCBeOWVV3DgwAFUVlYiPz8fP/74o5KjVX1tmvqlq6uLlJQUeHt7Q19fH5s2bZIUhAaeDjc0NS5bXV0t2d4ca2trWFtbtyU8jfNsDYINGzYgMzOz2RoEiYmJyMjIwOnTp2FsbCy1LSMjA66urnBzc8Pdu3exaNEimfoERHmcnWxbPD1L3uO3hLGxMc6ePYsvvvgClZWVsLS0RHh4OD799FPU19fjyJEjcHBwQK9evTB+/Hj8/vvvCotZE7R5nu2RI0cAPE2gubm5UpOdbW1tcfeu7OTsxjY7O7u2vn2Hk5ycjClTpsDS0hIODg5SNQhWrFghuVm2ZMkS6OrqwtnZWbLvd999h8jISNy7dw+zZ89GSUkJTE1NERAQgC+++EJZp0Se09IFB4pmb2+PXbt2Nbt9z549zW5To6dtcaZNyTY7OxtLly5FdHQ0rly5gmnTpuHatWswNTUFAPTt2xenT59GQ0OD1E2y8+fPw9DQUDK2S1qupTUIXvRhnzZtGqZNm6aQ+AghTZN7zLa2thZRUVGws7PDunXrkJqaipKSEsyZM0fSJzw8HCUlJfj5558lbffv38euXbswatSoF44jEkKIJpH7yjYxMRFXrlxBZmYmjI2N0adPH3z22WdYvHgxwsPDERgYiPDwcAwaNAjR0dG4fv06OnfujOTkZNTX17dppsHgEepx15ynwwOra9vXqZaOrxFCVJtcyTYrKwsrVqzAzJkzpR7I9+mnnyI9PR3vvfcecnJyYGZmhoMHD2L+/PlYv349xGIxvL29kZqaip49e8odtMtyF+jb68u9P1cEUQIauyKEAJBzGKF///6ora3F+vXrpdq1tbVx4cIFFBcXw8zMDABgbm6OH374Affv38ejR49w6tQpDBgwoM2BqyuhUIigoCAYGRmhZ8+eyMzMbLLf7t27MWjQIOjr6yMqKkpm+7179xAeHg5TU1NYWlpi4cKFCo6cENIWalliUZ09W0jm+PHjiIiIQG5uLiwsLKT6WVhY4OOPP8bZs2ebLDQTHByMcePGYcuWLeDxeDRPlhAVR8mWQ42FZPLz82UKyURHR0v1bZz3mpeXJ5NsDx06BD09PcydO1fS1txKPEKIaqAnNXBI3kIyz7tw4QK6du2KgIAAdO7cGSNGjMBff/3V3uESohSffPIJLCws8Oqrryo7lHZFyZZD8hSSaUpxcTG2b9+O2bNn486dO/D390doaCjq6+vbM1yiJHaOduDxeAp72Tm2bDFR165d4ezsjNraWknbBx98gPj4+Jfu+91338HNzQ18Ph92dnaIiIho0XsWFhYiOTkZeXl5uHTpUov2URc0jMAheQrJNMXAwABDhw5FQEAAAODjjz9GYmIi8vLy2jTLg6iGu7fvwjPVU2HHF0QJWty3qqoKmzZtwvTp01u8z6lTp7B06VIcOXIEnp6eKC0txS+//NKifQsLC9GlSxeZexiagK5sOdTaQjLN8fT0lKpBQYiizJkzBytWrJC6um20e/dueHh4wMLCAsHBwZIHO168eBGvvfYaPD2f/oNhbW0ttWKxa9euOHPmjOT3xqc+nD59Gm+99Rby8/PB5/Mxb948BZ8dtyjZcujZQjJisRgZGRnNFpKpr69HdXU16urqpH4GgNGjR+Pq1as4fvw46uvrsXbtWtjY2MDFxYXrUyIabvjw4XByckJqaqpU+4ULF/DRRx9h+/btKCkpgbu7O2JiYgA8Lb+6f/9+JCYm4vz5800m6qa89tprOHToELp37w6RSIQvv/yyvU9HqSjZciw5ORl37tyBpaUl5s6dK1VI5tkr3LS0NBgYGGDJkiX473//CwMDA8kznywtLbF7927MmjULZmZm2Lt3L/bu3duiB/ER0lpxcXEyV7cpKSmIiYmBl5cXOnXqhCVLliA9PR11dXUYNmwYtm/fjt9//x0jRoyAtbW11PPKOioas+VYSwvJREVFNbmYodHw4cNpBgLhxIgRI2Bvb4/NmzdL2goLC5GWliZVLU5HRwf37t2Dg4MDRo0ahVGjRqGurg779u3DhAkTMGDAAPj7+yvjFFSCWibbvEXqMYHf1oHqGhDNEBcXh/fff1+yPN/e3h7Lli2TmuvdFB0dHYSHh+Pzzz+HQCCAv78/jIyMZJ7y0BGGwNRyGEEgeFpzQNVfd4pUoy4pIW311ltvwcbGRvK0lOjoaHz11Ve4evUqgKcPe218ynN6ejp27dqFiooKMMZw9OhR5OTkSB5r/sorr2Dnzp2or6/H8ePHcerUKWWcEufUMtkSQrgXFxcnWc3o6+uL1atXY/LkyTAxMUH//v0lT2owNzfHN998g27dusHU1BRz5szB+vXrJc8jTEhIwOXLl2FmZoYff/yxyRvEmojH1KgsVU5ODjw9PeWaLkWIKsrPzwcAdO/eXdJm52iHu7dln3DSXmwdbNX2W1dTf16AeuQGtRyzJUSTqWsiJC9GwwiEEMIBSraEEMIBSraEKJka3TZROnX+s6JkS4gSaWlpoa6uTrIUmzSv8c/p2Sd1qxO6QUaIEpmamuLevXu4efMmdHR0qMBQMxhjqKurQ0NDA0xNTZUdjlwo2RKiRObm5gCAiooKNDQ0KDka1cXj8aCnpwdTU1PJn5m6oWRLiJKZm5urbQIhLaeegx+EEKJmKNkSQggHKNkSQggH1GrMtqamBsDTx3sTQkijxpzQmCNUkVol26KiIgBAaGiocgMhhKikoqIi9O/fX9lhNEmtqn6Vl5fj119/haOjI/T09JQdDiFERdTU1KCoqAh+fn4wMzNTdjhNUqtkSwgh6opukBFCCAco2RJCCAco2RJCCAco2RJCCAco2RJCCAco2RJCCAco2RJCCAco2RJCCAco2RJCCAf+P8buVWUJVvrwAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "nec = []\n", + "suf = []\n", + "nesuf = []\n", + "for name in feature_names:\n", + " score = score_dic[name]\n", + " nec.append(score[0])\n", + " suf.append(score[1])\n", + " nesuf.append(score[2])\n", + "score_df = pd.DataFrame()\n", + "score_df['feature'] = feature_names\n", + "score_df['Nec'] = nec\n", + "score_df['Suf'] = suf\n", + "score_df['NeSuf'] = nesuf\n", + "\n", + "plot_scores_barh_featurecol(score_df.sort_values(['NeSuf'], ascending=False))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5392beaf-5021-40af-b929-f6bc27f9bf82", + "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.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/lingam/__init__.py b/lingam/__init__.py index 3a9e108..190808a 100644 --- a/lingam/__init__.py +++ b/lingam/__init__.py @@ -27,6 +27,7 @@ from .missingness_lingam import mLiNGAM from .multi_group_resit import MultiGroupRESIT from .abic_lingam import ABICLiNGAM +from .lewis import LEWIS __all__ = [ "ICALiNGAM", @@ -56,6 +57,7 @@ "mLiNGAM", "MultiGroupRESIT", "ABICLiNGAM", + "LEWIS", ] __version__ = "1.12.1" diff --git a/lingam/lewis.py b/lingam/lewis.py new file mode 100644 index 0000000..4350348 --- /dev/null +++ b/lingam/lewis.py @@ -0,0 +1,252 @@ +""" +Python implementation of the LiNGAM algorithms. +The LiNGAM Project: https://sites.google.com/view/sshimizu06/lingam +""" + +import numpy as np +import pandas as pd +from sklearn.ensemble import RandomForestClassifier + + +class LEWIS(object): + """LEWIS explainer for computing necessity, sufficiency, and necessity-and-sufficiency scores. [1]_ [2]_ + + References + ---------- + .. [1] Sainyam Galhotra, Romila Pradhan, Babak Salimi (2021). + Explaining Black-Box Algorithms Using Probabilistic Contrastive Counterfactuals. + SIGMOD '21: International Conference on Management of Data, Virtual Event, China, June 20-25, 2021. + .. [2] https://sainyamgalhotra.github.io/lewis.zip + """ + + def __init__(self, epsilon=1e-10, random_state=0): + """Initialize LEWIS explainer. + + Parameters + ---------- + epsilon : float, optional (default=1e-10) + Small constant to avoid division by zero. + random_state : int, optional (default=0) + Random state for reproducibility. + """ + self._model_map = {} + self._epsilon = epsilon + self._random_state = random_state + + def _get_prob( + self, df, conditional_names, conditional_values, target_names, target_values + ): + """ + Estimate conditional probability. + + Parameters + ---------- + df : pandas.DataFrame + Input dataframe. + conditional_names : list of str + List of feature names for conditioning. + conditional_values : list + List of values corresponding to the conditional features. + target_names : list of str + List of target column names (typically one element). + target_values : list + List of values that the target should match. + + Returns + ------- + prob : float + Estimated conditional probability. + """ + # Check model cache + cache_key = ( + tuple(conditional_names), + tuple(target_names), + tuple(target_values), + ) + if cache_key in self._model_map: + model = self._model_map[cache_key] + prediction = model.predict_proba([conditional_values])[0][1] + return float(prediction) + + # Label target values + mask = np.ones(len(df), dtype=bool) + for t, v in zip(target_names, target_values): + mask &= df[t].astype(int) == int(v) + new_lst = mask.astype(int) + count = new_lst.sum() + + # If no conditions, simply return the occurrence rate of target values + if len(conditional_names) == 0: + return count * 1.0 / df.shape[0] + + # If all targets have the same value, probability is deterministic + if len(list(set(new_lst))) == 1: + if new_lst[0] == 1: + return 1 + else: + return 0 + + X = df[conditional_names].values + + model = RandomForestClassifier(random_state=self._random_state) + model.fit(X, new_lst) + self._model_map[cache_key] = model + prediction = model.predict_proba([conditional_values])[0][1] + return float(prediction) + + def get_scores( + self, + df, + x_names, + x_values, + x_prime_values, + o_name, + k_names=[], + k_values=[], + c_names=[], + ): + """Compute LEWIS explanation scores (Necessity, Sufficiency, and Necessity-and-Sufficiency). + + Parameters + ---------- + df : pandas.DataFrame + Input data frame. + x_names : list of str + Name of the attribute (or set of attributes) X under causal evaluation. + x_values : list + Target (intervened) value x of X, typically representing an improved or alternative value. + x_prime_values : list + Baseline or contrastive value x' of X against which x is compared. + o_name : str + Name of the outcome variable produced by the black-box model. + k_names : list of str, optional (default=[]) + Names of contextual variables defining the conditioning set K (used for global, local, or contextual explanations). + k_values : list, optional (default=[]) + Values corresponding to k_names, forming the concrete context k. + c_names : list of str, optional (default=[]) + Names of adjustment variables C satisfying the backdoor criterion. + + Returns + ------- + necessity : float + The necessity score. + sufficiency : float + The sufficiency score. + necessity_and_sufficiency : float + The necessity and sufficiency score. + """ + # Check parameters + if not isinstance(df, pd.DataFrame): + raise ValueError("df must be a pandas DataFrame") + if not isinstance(x_names, list) or not all(isinstance(v, str) for v in x_names): + raise ValueError("x_names must be a list of strings") + if not isinstance(x_values, list): + raise ValueError("x_values must be a list") + if not isinstance(x_prime_values, list): + raise ValueError("x_prime_values must be a list") + if not isinstance(o_name, str): + raise ValueError("o_name must be a string") + if not isinstance(k_names, list) or not all(isinstance(v, str) for v in k_names): + raise ValueError("k_names must be a list of strings") + if not isinstance(k_values, list): + raise ValueError("k_values must be a list") + if not isinstance(c_names, list) or not all(isinstance(v, str) for v in c_names): + raise ValueError("c_names must be a list of strings") + if len(x_names) != len(x_values): + raise ValueError( + f"x_names ({len(x_names)}) and x_values ({len(x_values)}) must have the same length" + ) + if len(x_names) != len(x_prime_values): + raise ValueError( + f"x_names ({len(x_names)}) and x_prime_values ({len(x_prime_values)}) must have the same length" + ) + if len(k_names) != len(k_values): + raise ValueError( + f"k_names ({len(k_names)}) and k_values ({len(k_values)}) must have the same length" + ) + + c_names_ = [v for v in c_names if v not in k_names] + + if len(c_names_) == 0: + c_values = [tuple()] + elif len(c_names_) == 1: + c_values = [(val,) for val in df[c_names_[0]].unique()] + else: + c_values = list(df.groupby(c_names_).groups.keys()) + + p_o_doxk = 0 + p_o_doxpk = 0 + nesuf = 0 + self._model_map = {} + for c_value in c_values: + cxk_names = k_names + c_names_ + x_names + cxk_values = k_values + list(c_value) + x_values + cxpk_values = k_values + list(c_value) + x_prime_values + + # P[o|cxk] + p_o_cxk = self._get_prob(df, cxk_names, cxk_values, [o_name], [1]) + + # P[o|cx'k] + if p_o_cxk > self._epsilon: + p_o_cxpk = self._get_prob(df, cxk_names, cxpk_values, [o_name], [1]) + else: + p_o_cxpk = 0 + continue + + if len(c_names_) > 0: + # P[c,k] + p_ck = self._get_prob(df, k_names, k_values, c_names_, c_value) + + # P[c|xk] + xk_names = k_names + x_names + xk_values = k_values + x_values + p_c_xk = self._get_prob(df, xk_names, xk_values, c_names_, c_value) + + # P[c|x'k] + xpk_values = k_values + x_prime_values + p_c_xpk = self._get_prob(df, xk_names, xpk_values, c_names_, c_value) + else: + p_ck = 1 + p_c_xpk = 1 + p_c_xk = 1 + + # P[o'|cx'k] + p_op_cxpk = 1 - p_o_cxpk + + p_o_doxk += p_o_cxk * p_c_xpk # P[o|do(x),k] + p_o_doxpk += p_op_cxpk * p_c_xk # P[o|do(x'),k] + nesuf += (p_o_cxk - p_o_cxpk) * p_ck # P[o|do(x),k]-P[o|do(x'),k] + + # P[o|k] + if len(k_names) > 0: + p_o_k = self._get_prob(df, k_names, k_values, [o_name], [1]) + else: + p_o_k = df[df[o_name] == 1].shape[0] * 1.0 / df.shape[0] + + xk_names = k_names + x_names + xk_values = k_values + x_values + xpk_values = k_values + x_prime_values + + # P[o|xk], P[o'|xk] + p_o_xk = self._get_prob(df, xk_names, xk_values, [o_name], [1]) + p_op_xk = 1 - p_o_xk + + # P[o|x'k], P[o'|x'k] + p_o_xpk = self._get_prob(df, xk_names, xpk_values, [o_name], [1]) + p_op_xpk = 1 - p_o_xpk + + if p_o_xk > self._epsilon: + nec = (p_o_doxpk - p_op_xk) * 1.0 / p_o_xk + else: + nec = 0.0 + + if p_op_xpk > self._epsilon: + suf = (p_o_doxk - p_o_xpk) * 1.0 / p_op_xpk + else: + suf = 0.0 + + return ( + float(np.clip(nec, 0.0, 1.0)), + float(np.clip(suf, 0.0, 1.0)), + float(np.clip(nesuf, 0.0, 1.0)), + ) diff --git a/tests/test_lewis.py b/tests/test_lewis.py new file mode 100644 index 0000000..38b1232 --- /dev/null +++ b/tests/test_lewis.py @@ -0,0 +1,217 @@ +import pandas as pd +from lingam.lewis import LEWIS + + +def test_get_scores_success(): + df = pd.DataFrame( + { + "x1": [0, 1, 0, 1, 0, 1, 0, 1, 0, 1], + "x2": [0, 0, 0, 0, 0, 1, 1, 1, 1, 1], + "x3": [0, 3, 0, 3, 3, 3, 0, 3, 2, 1], + "x4": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + "x5": [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + "y": [1, 0, 1, 0, 1, 0, 1, 0, 1, 0], + } + ) + explainer = LEWIS() + explainer.get_scores(df, ["x3"], [0], [2], "y") + explainer.get_scores(df, ["x3"], [3], [0], "y") + explainer.get_scores(df, ["x3"], [3], [0], "y", [], [], ["x1"]) + explainer.get_scores(df, ["x3"], [3], [0], "y", [], [], ["x1", "x2"]) + explainer.get_scores(df, ["x3"], [3], [0], "y", ["x2"], [1], ["x1"]) + explainer.get_scores(df, ["x1"], [0], [1], "x4") + explainer.get_scores(df, ["x1"], [0], [1], "x5") + + +def test_get_score_invalid_inputs(): + df = pd.DataFrame({"X": [0, 1, 0], "Y": [1, 0, 1]}) + explainer = LEWIS() + + # df must be a pandas DataFrame + try: + explainer.get_scores( + df.values, + x_names=["X"], + x_values=[1], + x_prime_values=[0], + o_name="Y", + ) + except ValueError: + pass + else: + raise AssertionError + + # x_names must be a list of strings + try: + explainer.get_scores( + df, + x_names=1, + x_values=[1], + x_prime_values=[0], + o_name="Y", + ) + except ValueError: + pass + else: + raise AssertionError + + # x_values must be a list + try: + explainer.get_scores( + df, + x_names=["X"], + x_values=1, + x_prime_values=[0], + o_name="Y", + ) + except ValueError: + pass + else: + raise AssertionError + + # x_prime_values must be a list + try: + explainer.get_scores( + df, + x_names=["X"], + x_values=[1], + x_prime_values=1, + o_name="Y", + ) + except ValueError: + pass + else: + raise AssertionError + + # o_name must be a string + try: + explainer.get_scores( + df, + x_names=["X"], + x_values=[1], + x_prime_values=[0], + o_name=1, + ) + except ValueError: + pass + else: + raise AssertionError + + # k_names must be a list of strings + try: + explainer.get_scores( + df, + x_names=["X"], + x_values=[1], + x_prime_values=[0], + o_name="Y", + k_names=1, + k_values=[], + c_names=[], + ) + except ValueError: + pass + else: + raise AssertionError + + # k_values must be a list + try: + explainer.get_scores( + df, + x_names=["X"], + x_values=[1], + x_prime_values=[0], + o_name="Y", + k_names=[], + k_values=1, + c_names=[], + ) + except ValueError: + pass + else: + raise AssertionError + + # c_names must be a list of strings + try: + explainer.get_scores( + df, + x_names=["X"], + x_values=[1], + x_prime_values=[0], + o_name="Y", + k_names=[], + k_values=[], + c_names=1, + ) + except ValueError: + pass + else: + raise AssertionError + + # x_names and x_values must have the same length + try: + explainer.get_scores( + df, + x_names=["X", "Z"], + x_values=[1], + x_prime_values=[0, 1], + o_name="Y", + k_names=[], + k_values=[], + c_names=[], + ) + except ValueError: + pass + else: + raise AssertionError + + # x_names and x_prime_values must have the same length + try: + explainer.get_scores( + df, + x_names=["X", "Z"], + x_values=[1, 2], + x_prime_values=[0], + o_name="Y", + k_names=[], + k_values=[], + c_names=[], + ) + except ValueError: + pass + else: + raise AssertionError + + # k_names and k_values must have the same length + try: + explainer.get_scores( + df, + x_names=["X", "Z"], + x_values=[1, 2], + x_prime_values=[0, 1], + o_name="Y", + k_names=["K1"], + k_values=[], + c_names=[], + ) + except ValueError: + pass + else: + raise AssertionError + + # c_names must be a list of strings + try: + explainer.get_scores( + df, + x_names=["X", "Z"], + x_values=[1, 2], + x_prime_values=[0, 1], + o_name="Y", + k_names=[], + k_values=[], + c_names=["C1", 1], + ) + except ValueError: + pass + else: + raise AssertionError