"
+ ],
+ "text/plain": [
+ "cue A B C D X Y\n",
+ "mean rating 9.84 1.47 0.35 0.81 4.97 8.42"
+ ]
+ },
+ "execution_count": 3,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "cue_names = [\"A\", \"B\", \"C\", \"D\", \"X\", \"Y\"]\n",
+ "\n",
+ "participant_means = ratings.groupby([\"ppt\", \"cue\"]).rating.mean().unstack()[cue_names]\n",
+ "observed = participant_means.mean()\n",
+ "observed.round(2).to_frame(\"mean rating\").T"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b68011b3",
+ "metadata": {},
+ "source": [
+ "Participants rated the blocked cue X lower than its control Y, so learning about X was blocked.\n",
+ "But X was not rated low: its mean rating is close to 5, the middle of the scale.\n",
+ "\n",
+ "## The models\n",
+ "\n",
+ "We use the two learning rules of the [blocking example](blocking.ipynb).\n",
+ "$V_i$ is the associative strength of cue $i$, $R$ is the outcome, and $\\alpha$ is the learning rate, which is the same for all cues.\n",
+ "\n",
+ "**Bush and Mosteller (1951).** Each cue that is present learns from its own error, independently of the other cues:\n",
+ "\n",
+ "$$\n",
+ "\\Delta V_i = \\alpha \\, (R - V_i)\n",
+ "$$\n",
+ "\n",
+ "**Rescorla and Wagner (1972).** All cues that are present share one error, the difference between the outcome and their summed prediction:\n",
+ "\n",
+ "$$\n",
+ "\\Delta V_i = \\alpha \\left( R - \\sum_{j \\in \\text{present}} V_j \\right)\n",
+ "$$\n",
+ "\n",
+ "Spicer et al. write the learning rate as the product of a cue salience and an outcome learning rate, and collapse the two into one parameter, $G$, which is our $\\alpha$.\n",
+ "\n",
+ "Participants gave ratings, not associative strengths.\n",
+ "Following Spicer et al., who in turn follow Gluck and Bower (1988), the rating $P_i$ of cue $i$ is a logistic function of its associative strength:\n",
+ "\n",
+ "$$\n",
+ "P_i = \\frac{10}{1 + e^{-\\theta \\, (V_i - \\beta)}}\n",
+ "$$\n",
+ "\n",
+ "$\\beta$ is the associative strength that is rated 5, the middle of the scale, and $\\theta$ sets how steeply the rating rises around it.\n",
+ "\n",
+ "Before training, every cue has the same starting strength $V_0$:\n",
+ "\n",
+ "$$\n",
+ "V_i = V_0 \\quad \\text{for all cues } i\n",
+ "$$\n",
+ "\n",
+ "As usually applied, both models start at $V_0 = 0$.\n",
+ "The model has four parameters: $\\alpha$, $\\beta$, $\\theta$, and, when it is free, $V_0$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d34c13cd",
+ "metadata": {},
+ "source": [
+ "We extend the model function of the [blocking example](blocking.ipynb) in three ways:\n",
+ "\n",
+ "- on the first trial, all six associative strengths are set to the starting strength,\n",
+ "- the model predicts a rating from the summed strength of the cues present,\n",
+ "- on test trials there is no feedback, so the strengths are not updated.\n",
+ "\n",
+ "Like Spicer et al., we fit ratings divided by 10, so that they are on the same scale as associative strength."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "24e40a59",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-10-01T08:52:15.474443Z",
+ "iopub.status.busy": "2026-10-01T08:52:15.474443Z",
+ "iopub.status.idle": "2026-10-01T08:52:15.478487Z",
+ "shell.execute_reply": "2026-10-01T08:52:15.478487Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "from functools import partial\n",
+ "\n",
+ "from cpm.generators import Parameters, Value, Wrapper\n",
+ "from cpm.models.learning import DeltaRule, SeparableRule\n",
+ "from cpm.models.utils import Nominal\n",
+ "\n",
+ "\n",
+ "def model(parameters, trial, rule):\n",
+ " values = np.asarray(parameters.values).copy()\n",
+ " if trial.trial == 1:\n",
+ " values[:] = parameters.initial_strength\n",
+ " cues = np.array([trial.cue_1, trial.cue_2], dtype=int)\n",
+ " present = Nominal(target=cues[cues > 0], bits=6)\n",
+ "\n",
+ " prediction = np.sum(values * present)\n",
+ " rating = 1 / (1 + np.exp(-parameters.theta * (prediction - parameters.beta)))\n",
+ " if trial.phase < 3: # no feedback at test\n",
+ " update = rule(weights=values, feedback=[trial.outcome], input=present, alpha=parameters.alpha)\n",
+ " update.compute()\n",
+ " values += update.weights.flatten()\n",
+ "\n",
+ " return {\"values\": values, \"prediction\": prediction, \"rating\": rating, \"dependent\": np.array([rating])}\n",
+ "\n",
+ "\n",
+ "rules = {\"Bush and Mosteller\": SeparableRule, \"Rescorla and Wagner\": DeltaRule}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "898423ec",
+ "metadata": {},
+ "source": [
+ "Next, we write the experiment as a table, with one row per trial, as in the [blocking example](blocking.ipynb).\n",
+ "The cues are numbered A = 1, B = 2, C = 3, D = 4, X = 5 and Y = 6.\n",
+ "Each cue is tested once: learning stops at test, so a second test block would give the same predictions.\n",
+ "\n",
+ "The dataset does not include the order in which each participant saw the training trials.\n",
+ "We use one fixed order instead; with learning rates as small as the ones we find below, the order does not change the fitted ratings."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "8547045f",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-10-01T08:52:15.481093Z",
+ "iopub.status.busy": "2026-10-01T08:52:15.480088Z",
+ "iopub.status.idle": "2026-10-01T08:52:15.487963Z",
+ "shell.execute_reply": "2026-10-01T08:52:15.487963Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "
phase
\n",
+ "
cue_1
\n",
+ "
cue_2
\n",
+ "
outcome
\n",
+ "
trial
\n",
+ "
\n",
+ " \n",
+ " \n",
+ "
\n",
+ "
52
\n",
+ "
2
\n",
+ "
2
\n",
+ "
6
\n",
+ "
1
\n",
+ "
53
\n",
+ "
\n",
+ "
\n",
+ "
53
\n",
+ "
2
\n",
+ "
3
\n",
+ "
4
\n",
+ "
0
\n",
+ "
54
\n",
+ "
\n",
+ "
\n",
+ "
54
\n",
+ "
3
\n",
+ "
1
\n",
+ "
0
\n",
+ "
0
\n",
+ "
55
\n",
+ "
\n",
+ "
\n",
+ "
55
\n",
+ "
3
\n",
+ "
2
\n",
+ "
0
\n",
+ "
0
\n",
+ "
56
\n",
+ "
\n",
+ "
\n",
+ "
56
\n",
+ "
3
\n",
+ "
3
\n",
+ "
0
\n",
+ "
0
\n",
+ "
57
\n",
+ "
\n",
+ "
\n",
+ "
57
\n",
+ "
3
\n",
+ "
4
\n",
+ "
0
\n",
+ "
0
\n",
+ "
58
\n",
+ "
\n",
+ "
\n",
+ "
58
\n",
+ "
3
\n",
+ "
5
\n",
+ "
0
\n",
+ "
0
\n",
+ "
59
\n",
+ "
\n",
+ "
\n",
+ "
59
\n",
+ "
3
\n",
+ "
6
\n",
+ "
0
\n",
+ "
0
\n",
+ "
60
\n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ],
+ "text/plain": [
+ " phase cue_1 cue_2 outcome trial\n",
+ "52 2 2 6 1 53\n",
+ "53 2 3 4 0 54\n",
+ "54 3 1 0 0 55\n",
+ "55 3 2 0 0 56\n",
+ "56 3 3 0 0 57\n",
+ "57 3 4 0 0 58\n",
+ "58 3 5 0 0 59\n",
+ "59 3 6 0 0 60"
+ ]
+ },
+ "execution_count": 5,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "number = {cue: i + 1 for i, cue in enumerate(cue_names)}\n",
+ "\n",
+ "\n",
+ "def trials(phase, stimuli, outcomes):\n",
+ " return [\n",
+ " {\"phase\": phase, \"cue_1\": number[s[0]], \"cue_2\": number[s[1]] if len(s) > 1 else 0, \"outcome\": outcome}\n",
+ " for s, outcome in zip(stimuli, outcomes)\n",
+ " ]\n",
+ "\n",
+ "\n",
+ "stage_1 = trials(1, [\"A\", \"B\", \"C\"], [1, 0, 0]) * 12\n",
+ "stage_2 = trials(2, [\"AX\", \"BY\", \"CD\"], [1, 1, 0]) * 6\n",
+ "test = trials(3, cue_names, [0] * 6)\n",
+ "design = pd.DataFrame(stage_1 + stage_2 + test)\n",
+ "design[\"trial\"] = np.arange(1, len(design) + 1)\n",
+ "design.tail(8)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "389f76a0",
+ "metadata": {},
+ "source": [
+ "## Fitting the mean ratings\n",
+ "\n",
+ "Spicer et al. fitted one set of parameters to the six mean ratings of the group, by minimising the sum of squared errors (SSE) between the predicted and observed ratings.\n",
+ "To do this with [`FminBound`](#cpm.optimisation.FminBound), we treat the whole group as a single participant: the `observed` column holds the mean rating of each cue on its test trial, and nothing on the training trials, where no ratings were made.\n",
+ "\n",
+ "The squared-error loss of [`Distance.SSE`](#cpm.optimisation.minimise.Distance.SSE) does not accept missing values, so we write a loss function that passes it the test trials only."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "88023cb0",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-10-01T08:52:15.489970Z",
+ "iopub.status.busy": "2026-10-01T08:52:15.489970Z",
+ "iopub.status.idle": "2026-10-01T08:52:15.494486Z",
+ "shell.execute_reply": "2026-10-01T08:52:15.494486Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "from cpm.optimisation import FminBound, minimise\n",
+ "\n",
+ "design[\"ppt\"] = \"group\"\n",
+ "design[\"observed\"] = np.nan\n",
+ "test_trials = design.phase.eq(3).to_numpy()\n",
+ "design.loc[test_trials, \"observed\"] = observed[cue_names].to_numpy() / 10\n",
+ "\n",
+ "\n",
+ "def test_sse(predicted, observed, **kwargs):\n",
+ " rated = ~np.isnan(np.ravel(observed))\n",
+ " return minimise.Distance.SSE(np.ravel(predicted)[rated], np.ravel(observed)[rated])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "46200890",
+ "metadata": {},
+ "source": [
+ "The parameters get bounds but no informative priors, since we minimise the SSE alone.\n",
+ "The starting strength is a fixed value, or a free parameter between −1 and 1, as in Spicer et al."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "2ea03fe6",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-10-01T08:52:15.496384Z",
+ "iopub.status.busy": "2026-10-01T08:52:15.496384Z",
+ "iopub.status.idle": "2026-10-01T08:52:15.499285Z",
+ "shell.execute_reply": "2026-10-01T08:52:15.499285Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "def make_parameters(initial_strength):\n",
+ " return Parameters(\n",
+ " alpha=Value(value=0.1, lower=0, upper=1, prior=\"uniform\"),\n",
+ " beta=Value(value=0.5, lower=0, upper=1, prior=\"uniform\"),\n",
+ " theta=Value(value=5, lower=0, upper=100, prior=\"uniform\"),\n",
+ " initial_strength=initial_strength,\n",
+ " values=np.zeros(6),\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7efb49d4",
+ "metadata": {},
+ "source": [
+ "The fit can end in a local minimum, so we start the optimiser from a grid of starting values and keep the best fit, as Spicer et al. did.\n",
+ "The error surface of the Bush–Mosteller model is flat along some directions, so we also tighten the stopping rule of the optimiser.\n",
+ "The function below fits both models, and returns their best-fitting parameters, their SSE, and the $R^2$ of the predicted against the observed mean ratings, which Spicer et al. used to judge how well a model fits."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "b509acac",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-10-01T08:52:15.501292Z",
+ "iopub.status.busy": "2026-10-01T08:52:15.501292Z",
+ "iopub.status.idle": "2026-10-01T08:52:15.505161Z",
+ "shell.execute_reply": "2026-10-01T08:52:15.505161Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "from itertools import product\n",
+ "\n",
+ "\n",
+ "def fit_models(initial_strength, guesses):\n",
+ " summary, predictions = {}, {}\n",
+ " for name, rule in rules.items():\n",
+ " wrapper = Wrapper(model=partial(model, rule=rule), parameters=make_parameters(initial_strength), data=design)\n",
+ " fit = FminBound(\n",
+ " model=wrapper,\n",
+ " data=design,\n",
+ " ppt_identifier=\"ppt\",\n",
+ " minimisation=test_sse,\n",
+ " initial_guess=guesses,\n",
+ " number_of_starts=len(guesses),\n",
+ " # passed on to scipy.optimize.fmin_l_bfgs_b\n",
+ " approx_grad=True,\n",
+ " factr=10,\n",
+ " pgtol=1e-12,\n",
+ " )\n",
+ " fit.optimise(display=False)\n",
+ "\n",
+ " best = {key: float(value) for key, value in fit.parameters[0].items()}\n",
+ " wrapper.reset(parameters=best)\n",
+ " wrapper.run()\n",
+ " predictions[name] = 10 * wrapper.dependent[test_trials].ravel()\n",
+ " r_squared = np.corrcoef(predictions[name], observed[cue_names])[0, 1] ** 2\n",
+ " summary[name] = {\"initial_strength\": float(np.asarray(initial_strength)), **best, \"SSE\": fit.fit[0][\"fun\"], \"R²\": r_squared}\n",
+ " return pd.DataFrame(summary).T, pd.DataFrame(predictions, index=cue_names)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f528749f",
+ "metadata": {},
+ "source": [
+ "Spicer et al. show the fits as the predicted mean rating of each cue, against the distribution of the participants' ratings.\n",
+ "We draw the same figure: the grey dots are the mean ratings of individual participants, the black line is the mean of the group, and the two models are drawn on either side of it, so that the distance to the line is the error of the model."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "708ab22d",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-10-01T08:52:15.507169Z",
+ "iopub.status.busy": "2026-10-01T08:52:15.507169Z",
+ "iopub.status.idle": "2026-10-01T08:52:15.513761Z",
+ "shell.execute_reply": "2026-10-01T08:52:15.513761Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "SURFACE, INK, MUTED, GRID = \"#fcfcfb\", \"#0b0b0b\", \"#52514e\", \"#e6e5e1\"\n",
+ "styles = { # colour and shape, so that the models are not told apart by colour alone\n",
+ " \"Bush and Mosteller\": {\"color\": \"#2a78d6\", \"marker\": \"o\", \"offset\": -0.15},\n",
+ " \"Rescorla and Wagner\": {\"color\": \"#eb6834\", \"marker\": \"D\", \"offset\": 0.15},\n",
+ "}\n",
+ "\n",
+ "\n",
+ "def plot_fit(predictions, title):\n",
+ " positions = np.arange(len(cue_names))\n",
+ " jitter = np.random.default_rng(0).uniform(-0.07, 0.07, participant_means.shape)\n",
+ " fig, ax = plt.subplots(figsize=(7.5, 4.4), facecolor=SURFACE)\n",
+ " ax.set_facecolor(SURFACE)\n",
+ "\n",
+ " violins = ax.violinplot([participant_means[cue] for cue in cue_names], positions=positions, widths=0.8, showextrema=False)\n",
+ " for body in violins[\"bodies\"]:\n",
+ " body.set(facecolor=\"#eeede9\", edgecolor=\"none\", alpha=1)\n",
+ " ax.scatter((positions + jitter).ravel(), participant_means.to_numpy().ravel(), s=10, color=\"#a9a8a2\", linewidths=0, zorder=2, label=\"participants\")\n",
+ " ax.hlines(observed[cue_names], positions - 0.32, positions + 0.32, color=INK, linewidth=2.5, zorder=3, label=\"observed mean\")\n",
+ " for name, style in styles.items():\n",
+ " ax.scatter(\n",
+ " positions + style[\"offset\"], predictions[name], s=80, color=style[\"color\"], marker=style[\"marker\"],\n",
+ " edgecolors=SURFACE, linewidths=2, zorder=4, label=name,\n",
+ " )\n",
+ "\n",
+ " ax.set_xticks(positions, cue_names, fontsize=11, color=INK)\n",
+ " ax.set_yticks(range(0, 11, 2))\n",
+ " ax.set_ylim(-0.5, 10.5)\n",
+ " ax.set_xlim(-0.6, len(cue_names) - 0.4)\n",
+ " ax.set_xlabel(\"cue\", color=MUTED)\n",
+ " ax.set_ylabel(\"mean causal rating\", color=MUTED)\n",
+ " ax.tick_params(length=0, colors=MUTED)\n",
+ " ax.yaxis.grid(True, color=GRID, linewidth=0.8)\n",
+ " ax.set_axisbelow(True)\n",
+ " ax.spines[[\"top\", \"right\", \"left\"]].set_visible(False)\n",
+ " ax.spines[\"bottom\"].set_color(GRID)\n",
+ " ax.set_title(title, loc=\"left\", color=INK, fontsize=12, pad=30)\n",
+ " legend = ax.legend(\n",
+ " ncol=4, loc=\"lower left\", bbox_to_anchor=(0, 1.0), frameon=False, fontsize=9, labelcolor=MUTED,\n",
+ " handletextpad=0.3, columnspacing=1.2, borderaxespad=0.3,\n",
+ " )\n",
+ " legend.legend_handles[0].set_sizes([30]) # the participant dots are too small to read in the legend\n",
+ " plt.tight_layout()\n",
+ " plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5dbb52da",
+ "metadata": {},
+ "source": [
+ "The values that Spicer et al. report for their fits (Tables S2 and S3 of their supplementary materials) are below, for comparison.\n",
+ "They are rounded to two decimal places."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "id": "769939ba",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-10-01T08:52:15.515768Z",
+ "iopub.status.busy": "2026-10-01T08:52:15.515768Z",
+ "iopub.status.idle": "2026-10-01T08:52:15.520153Z",
+ "shell.execute_reply": "2026-10-01T08:52:15.520153Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "reported = pd.DataFrame(\n",
+ " [\n",
+ " [\"zero\", \"Bush and Mosteller\", 0.00, 0.04, 0.21, 13.41, 0.24, 0.71],\n",
+ " [\"zero\", \"Rescorla and Wagner\", 0.00, 0.01, 0.04, 75.81, 0.24, 0.71],\n",
+ " [\"free\", \"Bush and Mosteller\", 0.45, 0.02, 0.49, 30.52, 0.06, 0.93],\n",
+ " [\"free\", \"Rescorla and Wagner\", 0.43, 0.04, 0.41, 19.31, 0.00, 1.00],\n",
+ " ],\n",
+ " columns=[\"starting strength\", \"model\", \"initial_strength\", \"alpha\", \"beta\", \"theta\", \"SSE\", \"R²\"],\n",
+ ").set_index([\"starting strength\", \"model\"])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a98eaed4",
+ "metadata": {},
+ "source": [
+ "## All cues start at zero\n",
+ "\n",
+ "We first fit the models as they are usually applied, with every cue starting at zero.\n",
+ "There are three free parameters, so the grid has $2^3 = 8$ starting points."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "id": "b1571d35",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-10-01T08:52:15.522161Z",
+ "iopub.status.busy": "2026-10-01T08:52:15.522161Z",
+ "iopub.status.idle": "2026-10-01T08:52:39.246352Z",
+ "shell.execute_reply": "2026-10-01T08:52:39.246352Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plot_fit(zero_predictions, \"All cues start at zero\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "19dc7479",
+ "metadata": {},
+ "source": [
+ "This is Figure 1 of Spicer et al.: both models fit equally poorly, with an SSE of 0.24 and an $R^2$ of 0.71.\n",
+ "The Rescorla–Wagner model predicts blocking, since it rates X lower than Y, but the difference is small, and neither model captures the low rating of B or the high rating of Y.\n",
+ "\n",
+ "The problem is the starting strength.\n",
+ "After the B− trials of Stage 1, B still has an associative strength of zero, just like the new cue Y.\n",
+ "In Stage 2, B and Y are then trained together on BY+ trials, so both models predict the same rating for them, although participants rated B 1.5 and Y 8.4.\n",
+ "The best the models can do is to learn very little, and use a steep response function to place B, X and Y near the middle of the scale.\n",
+ "The fitted Rescorla–Wagner model takes this to the extreme: its learning rate is close to zero, and $\\theta$ is at its upper bound of 100, so that the response function is almost a step.\n",
+ "Its parameters differ from the reported ones for this reason; the fit is equally good.\n",
+ "\n",
+ "## Starting from uncertainty\n",
+ "\n",
+ "A rat that has never encountered a cue has not learned anything about it, so its associative strength should start at zero.\n",
+ "A participant who sees a new food, however, does not know whether it causes stomach ache or not.\n",
+ "Spicer et al. proposed that the associative strength of a novel cue should start at an intermediate value, which reflects this uncertainty, and they made the starting strength a free parameter.\n",
+ "If it represents uncertainty, the best-fitting value should be near the middle of the range of associative strengths.\n",
+ "\n",
+ "We make the starting strength a free parameter, with a uniform prior between −1 and 1, and add it to the grid of starting points."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "id": "eac403fd",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-10-01T08:52:39.481325Z",
+ "iopub.status.busy": "2026-10-01T08:52:39.481325Z",
+ "iopub.status.idle": "2026-10-01T08:54:02.425098Z",
+ "shell.execute_reply": "2026-10-01T08:54:02.425098Z"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
\n",
+ "
\n",
+ "
initial_strength
\n",
+ "
alpha
\n",
+ "
beta
\n",
+ "
theta
\n",
+ "
SSE
\n",
+ "
R²
\n",
+ "
\n",
+ " \n",
+ " \n",
+ "
\n",
+ "
cpm
\n",
+ "
Bush and Mosteller
\n",
+ "
0.417
\n",
+ "
0.006
\n",
+ "
0.429
\n",
+ "
87.679
\n",
+ "
0.060
\n",
+ "
0.928
\n",
+ "
\n",
+ "
\n",
+ "
Rescorla and Wagner
\n",
+ "
0.428
\n",
+ "
0.043
\n",
+ "
0.410
\n",
+ "
19.317
\n",
+ "
0.001
\n",
+ "
0.999
\n",
+ "
\n",
+ "
\n",
+ "
Spicer et al. (2021)
\n",
+ "
Bush and Mosteller
\n",
+ "
0.450
\n",
+ "
0.020
\n",
+ "
0.490
\n",
+ "
30.520
\n",
+ "
0.060
\n",
+ "
0.930
\n",
+ "
\n",
+ "
\n",
+ "
Rescorla and Wagner
\n",
+ "
0.430
\n",
+ "
0.040
\n",
+ "
0.410
\n",
+ "
19.310
\n",
+ "
0.000
\n",
+ "
1.000
\n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ],
+ "text/plain": [
+ " initial_strength alpha beta \\\n",
+ "cpm Bush and Mosteller 0.417 0.006 0.429 \n",
+ " Rescorla and Wagner 0.428 0.043 0.410 \n",
+ "Spicer et al. (2021) Bush and Mosteller 0.450 0.020 0.490 \n",
+ " Rescorla and Wagner 0.430 0.040 0.410 \n",
+ "\n",
+ " theta SSE R² \n",
+ "cpm Bush and Mosteller 87.679 0.060 0.928 \n",
+ " Rescorla and Wagner 19.317 0.001 0.999 \n",
+ "Spicer et al. (2021) Bush and Mosteller 30.520 0.060 0.930 \n",
+ " Rescorla and Wagner 19.310 0.000 1.000 "
+ ]
+ },
+ "execution_count": 13,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "initial_strength = Value(value=0.5, lower=-1, upper=1, prior=\"uniform\")\n",
+ "grid = [(0.05, 0.25), (0.25, 0.75), (1, 3), (0.25, 0.75)] # alpha, beta, theta, initial_strength\n",
+ "free_summary, free_predictions = fit_models(initial_strength, np.array(list(product(*grid))))\n",
+ "\n",
+ "pd.concat({\"cpm\": free_summary, \"Spicer et al. (2021)\": reported.loc[\"free\"]}).round(3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "id": "e51295f5",
+ "metadata": {
+ "execution": {
+ "iopub.execute_input": "2026-10-01T08:54:02.428172Z",
+ "iopub.status.busy": "2026-10-01T08:54:02.428172Z",
+ "iopub.status.idle": "2026-10-01T08:54:02.521726Z",
+ "shell.execute_reply": "2026-10-01T08:54:02.521726Z"
+ },
+ "tags": [
+ "thumbnail"
+ ]
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtIAAAGuCAYAAAC9YuGNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAACfjklEQVR4nOzdd5hTVfoH8O+9N20mUzK9wMzQhyKIqICKAqICoquigIIIWFYsu+5vbVhQEV27rq5iV5QiYMG2YsMGFlBZQFCGJjMDTJ/JMJOe3PP7IyRMMum56e/neXg0mZvk5La899z3vIezWi0MhBBCCCGEkKDwsW4AIYQQQgghiYgCaUIIIYQQQkJAgTQhhBBCCCEhoECaEEIIIYSQEFAgTQghhBBCSAgokCaEEEIIISQEFEgTQgghhBASAgqkCSGEEEIICQEF0oQQQgghhISAAukYGTp0KK677rpYNyPpXXfddSgt7RHy66dMmYIpU6ZI2KLQPf30Mxg27Hjk5ORizJgxsW4OIYQQkvKSNpDeuXMnZs++AscddxwKC4swcOAgXHDBhXjxxRddlnv88Sfw8ccfR6QNmzZtwkMPPQStVhuR94+UXbt24aGHHkJ1dXWsmxIQvV6Phx56CBs2bIh1UyJm/fqvcM8992D06FFYsuQ53HPPPbFuEgHw+eef46GHHop1MwghhMRIUgbSmzZtwrhx47Fjxw7MmTMHjz32GK64YjZ4nsfzz7/gsuyTTz6J//73vxFqx2Y8/PAjaG9v7/a3X375Bc8880xEPjdcu3ZV4eGHH0FNTU2smxIQg8GAhx9+BBs3bpT8vdeuXYu1a9dK/r7B+u6778DzPJ599llcdtllOOecc2LdJALg88+/wMMPPxLrZhBCCIkRWawbEAmPP/4EsrKy8PXXX0Gj0bj8rampKeKfr9PpoFarfS6jVCoj3o5oYIzBaDQiLS0t1k2JCIVCEesmAACam5uQlpbmtz2iKMJsNkOlUkWpZZFntVohimLcbItIS/ZjihBCkklS9kj/+eefGDRoYLcgGgAKCgqc/5+drYFOp8PKlW8hO1uD7GyNM2+5pqYG//znzTjxxJNQVFSMXr1644or5nRLd1ixYgWyszXYuHEj/vnPm9G3bz8MGjQYDz30EBYuXAgAGDbseOf7O17vniPteJ+ffvoJd955J/r06YuSklLMmjULzc3NLp8piiIeeughVFYORHFxCc477zzs2rUr4Lzrd955F2ecMRY9evREz55lOOWUU/H888872zFnzhwAwHnnne9styNtYujQoZg+fQa+/HI9xo4dh6KiYrz++usAAK1WiwULFmDw4CEoKCjE8OEn4Kmn/g1RFJ2fXV1djexsDZ555j94/fWlOP744SgoKMS4cePx669burV17dr3MXLkKBQWFmH06FPw0Ucf4brrrsPQoUOd79enT18AwMMPP+Jsr/vt9sOHD2PmzJkoLe2BPn364q677obNZvO7rjzlSL/44osYNWo0iotLUF5egbFjx+Htt9/2+T5msxkPPvggzjhjLMrKylFSUopJkybju+++89uG7GwNli9fAZ1O5/x+K1ascP7tlltuxZo1azBq1GgUFBTiyy+/dH7nG264Af369UdBQSFGjRqNZcuWdXt/k8mEf/3rXxg+/AQUFBRi8OAhWLjwHphMpoDWz+jRp+B//9uKs88+B0VFxRg6dBheffW1kL5/1/1jyZIlzv1j165dIb3Hyy+/jGHDjkdxcQkuvPAiHDx4EIwxPProoxg0aDCKiopx2WWXobW1rdt3++KLLzBp0mSUlJSiR4+emDZtOv744w/n36+77jq8/PLLzu3g+OcgiiKWLFmCUaNGo7CwCP369cdNN/0DbW1al8/xdUwRQgiJb0nZI11WVoaff/4Zv//+OwYPHux1uZdeehF/+9vfceKJIzB37lwAQO/evQEAW7b8D5s3b8LFF09FaWkP1NTU4NVXX8V5552HTZs2IT093eW9br75FuTn5+H222+DTqfH2Wefhb179+Gdd97BQw/9C3l5eQCA/Px8n22/9dbboNFosGDB7aiursHzzz8PufxWLF167If1vvsW4emnn8bkyZMwYcIE/PbbDkydejGMRqPfdfPVV1/jqquuwtixY3HfffcBAHbvrsJPP23Cddddh1NPPQ3z51+LF154ETfffDMqKwcAACorK53vsWfPHlx11VWYN28e5syZg/79+0Gv12PKlCk4fLgO8+bNQ8+ePbF58yYsWrQIDQ31ePjhh13a8fbbb6OzsxPz5s0Dx3F4+umnMXv2bGzbthVyuRwA8Nlnn2HevHkYMmQw7r33Hmi1Wtx4499QWlrqfJ/8/Hw8+eST+Oc//4nzzjsPf/nL+QCAIUOGOJex2WyYOvVinHjiiVi8eDG++eYbPPvss+jduzeuvvoqv+usq6VL38Btt92OCy64APPnz4fJZMSOHTvxyy+/YNq0aV5f19HRgTffXIZLLrkYc+bMQWdnJ5YtW4apUy/GV1+tx7Bhw7y+9qWXXsTSpUvx669b8J//2NOBRo4c5fz7d999h7Vr1+Kvf70Gubl5KC8vR2NjI84662xwHIe//vUa5OXl48svv8CNN/4NHR0duP766wHYg71LL70MP/30E+bOnYMBAyrx++87sWTJEuzbtxcrV670u060Wi2mTZuGiy66EJdccjHWrn0f//znP6FQyDF79uyQvv+KFStgNBoxd+5cKJUK5OTkBP0ea9asgcViwbXX/hVtbW14+ulnMHfuPJxxxhnYuHEj/vGPm7B//368+OJLWLjwbjz33HPO165atQrz51+HCRMmYNGi+6DXG/Dqq69i4sRJ2LDhO1RUVGDevHmoq6vH119/jZdech17AQA33fQPrFy5ErNmzcK1116L6upqvPzyy9i+fTs+//wz534OeD6mCCGEJACr1cKS7d/nn3/OsrKyWFZWFhs/fjy788472WeffcYMBn23ZQsLC9k111zd7fmOjiPdnvvxxx+YWq1my5cvcz73xhtLmVqtZmeddRYzmYwuyz/55JNMrVazffv2dnuvQYMGunyu432mTJnCLBaz8/nbbruVZWVlsZaWZma1Wtjhw4dYdnY2mzFjusv7PfDAYqZWqz1+l67/brnlZlZSUtKtrV3/vfPOO0ytVrOvv/7KY7vVajX79NNPXZ5/6KF/scLCQlZVtcvl+bvvvptlZWWxAwf+ZFarhe3bt5ep1WpWVlbGmpqanMt9+OEHTK1Ws48//sj53MiRJ7MBAwYwrbbN+dw333zN1Go1GzRooPO5hoZ6plar2eLF93dr7zXXXM3UajX7178edHn+lFNOYaeddprffWnixHPYxInnOB9Pnz6NnXTSSUHvkyaTken1OpfnmpubWO/evdm1117r9/XXXHM1Kyws7Pa8Wq1mmZmZbMeO31yenz9/Puvbty9rbGxwef6KK65gpaUlzv17xYrlLDMzk23Y8J3Lci+99BJTq9Xs++83+l0/arWa/fvfTzmf0+t1bPTo0axXr17OYy7Q7+/YP0pKSlh9fV1I69DxHhUVFc7jxmq1sIULFzK1Ws1Gjx7FjEaD8/k5c+awnJwcptN1MqvVwrTaNlZaWsquv/56l886fPgQKy0tcXn+H//4B1Or1d3Wy4YN3zG1Ws3eemuly/Offvppt+e9HVP0j/7RP/pH/+L/X1Kmdpx55nh88cUXmDx5Mnbs2IGnn34aU6dOxcCBg/DJJ58E9B5d8xMtFgtaW1vRp08fZGdnY9u2bd2WnzPnCgiCEHbb586dC47jnI9POeUU2Gw21NbWAgC+/fZbWK1WXHXV1S6v++tfrw3o/bOzs6HT6fD111+H3MaKigqcddYEl+fef/99nHLKKdBoNGhpaXH+GzduHGw2G3744QeX5adOnYqcHI3z8SmnnAoAOHDgAACgrq4OO3f+jksvvRQZGRnO5caMGYMhQ7zfZfDmyiuvdHl86qmnOD8rGNnZ2Th8+JDHNBRfBEFw5viKoojW1jbYbDaccMIJHvenYJx22mkYOHCg8zFjDB9++CEmTZoExpjL9pgw4Uy0tx9xfub777+PyspKDBgwwGW5sWPPAAB8953/SigymQzz5s1zPlYoFJg3bx6ampqwdevWkL7/X/5yfre7N8G+x4UXXojs7Gzn45NOOhEAMH36DMhkMpfnzWYzDh+uAwB8/fXXaG9vxyWXXOyyTgRBwIknnhRQdZj3338f2dlZGD9+vMt7nHDCcGRkZHRbr56OKUIIIfEvKVM7AODEE0dgxYrlMJvN+O23Hfj444+xZMkSXHHFHGzcuMEl8PDEYDDgySefxIoVK3H48GEwxpx/O3LkSLflKyoqJGl3WVlPl8eOPG9HCb2aGntA3adPb5flcnNzPOaEu7v66quxdu37uPjiS1BaWoozzxyPiy66CGeddVbAbfT0Xfft248dO3Y685XduQ/y7NnT9Xs6gmrH93RcOLh/TwDo3bsPtm8PPPhUqVTdgjKNRhNSWcJ//OMf+Oabb3HmmWeiT58+OPPMMzFt2iUYPXq039euXLkSzz77LHbv3gOLxeJ8Ptx9x/31zc3NaG9vx9KlS7F06VKPr3Fsj3379qOqqsrrdmtu9j84t6SkuNvg2n797O9XXV2Dk08+GUBw39/bOgnmPdz3saysLABAjx49PD7v2B/27dsPADj//L94bINjeV/27duP9vYj6NvXc4qG+3qV6vxBCCEkupI2kHZQKBQ48cQROPHEEejXry+uv/4GvP/++1iwYIHP1916621YsWIFrr/+Opx88khkZWWB4zhceeWVLoPnHFQqaUbY87znXu2ugXw4CgoKsHHjBqxfvx5ffPElvvjiCyxfvgKXXXYpXnjhBf9vAHisJiCKIsaPH4+bbrrJ42scgZWDt957qb5nIJ8VisrKSvzyy8/49NPPsH79l/jwww/xyiuv4Pbbb8Odd97p9XWrV6/Gddddj/POm4K///3vyM8vgCAIePLJJ3HgwJ9htcl9ezj2zxkzpuOyy2Z6fM1xxw1xLjtkyGA8+OC/PC7Xs2fok9l0Fez393Q8Bfse3ra7v33Psf5eeulFFBYWdVtOJvO/P4miiIKCAudgRHf5+Xkuj6lCByGEJKakD6S7OuGEEwAA9fUNzue6plF09eGHH2DmzMvw4IMPOp8zGo0ea0J74+Wtw1JeXgYA2L//T/Tq1cv5fGtra8A9rAqFApMnT8bkyZMhiiL++c+b8frrr+PWW29D3759Qmp37969odPpMH78uOBf7EFZ2bHv6e7PP/e7PPa2DSNFrVbj4oun4uKLp8JsNuPyy2fj8cefwD//+U+vZec++OAD9OrVC8uXL3dpbyQm88jPz0dmZiZsNtHv9ujduzd27NiBcePGhrwe6+rqu5V83Lt3HwCgoqIcgDTfP1rr0DHgOD+/wO/687bKevfujW+++QajR4+iIJkQQpJYUuZIf/fddx57Nj///AsAcBkRn56e7jE45nmh23u8+OJLAZVMO/be9sAimODbn7Fjx0Imk+HVV191ef6ll14K6PWtra0uj3med1a4MJvt5c4cAVEw7b7ooguxefNmfPnl+m5/02q1sFqtAb8XAJSUlGDw4MFYtWoVOjs7nc9v3LgRO3f+7rKsI1CRcj17477+FAoFKisrwRhzSTVw57jT0HWf+uWXX7B582bJ2ygIAv7yl/Px4Ycf4vfff+/2967lFC+66EIcPnwYS5e+0W05g8EAnU7n9/OsVqtLuTaz2YzXX38d+fn5GD58OABpvn+01uGECWciKysLTz75hMdt2nX9OY5x94vYiy66EDabDY8++li311ut1oSb7ZQQQohnSdkjfdttt0Ov1+P8889D//4DYLGYsWnTZrz33nsoLy/HrFmznMsOHz4c33zzLZ599lkUF5egV68KnHTSSZg0aSJWrVqNrKwsVFYOxM8/b8Y333yL3NzcgNvhCCIWL34AF188FTKZHJMnT/I7WYsvhYWFmD9/Pp599llceumlmDDhLOzYsQNffvkl8vLy/PYq/u1vf0NbWxvOOOMMlJaWora2Fi+++BKGDh3qLHE3dOhQCIKAf//7aRw5cgQKhRJjx57hUoPb3d///nd88sk6zJgxAzNnzsTw4cOh1+vw+++/44MPPsRvv213lgAM1D33LMRll83ExIkTMWvWLGi1Wrz00ssYPHiwS3CdlpaGgQMH4r331qJv337IycnB4MGDfJY+DNWFF16EoqIijBo1CoWFBaiq2o2XX34ZEyeeg8zMTK+vmzRpIj766CPMmjUL55wzEdXV1XjttdcwcOBA6HSdXl8Xqvvuuw8bNmzAhAlnYc6cK1BZORBtbW3Ytm0bvvnmG1RXHwAAXHrppVi79n383//9HzZs2IDRo0fBZrNh9+49WLt2Ld577z2MGHGCz88qKSnBv//9NGpqatCvXz+8995a/Pbbb3j66aedJd6k+P7RWoeOIPqvf70WZ5wxFlOnTkV+fj4OHjyIzz77DKNHj8bjj9sDZMcxfvvtt2PChAngeQGXXHIxxowZg3nz5uHJJ5/Eb7/9hjPPHA+5XI59+/bh/fc/wMMPP4wLL7xAsjYTQgiJjaQMpB94YDHWrn0fn3/+BZYufQNmsxk9e/bE1VdfhVtvvdVlUN6//vUgbrrpJjzwwIMwGAyYOfMynHTSSXj44YchCALWrHkbJpMJo0aNwgcfvI+pU6cG3I4TTxyBu+++C6+99jq+/PJLiKKI7du3hRVIA8D99y9Cenoa3njjTXzzzbc4+eSTsXbte5g4cRKUSt8z2k2fPh1Ll76BV155Fe3t7SgqKsTUqVNxxx0LwPP2GxRFRUV46qmn8OSTT+LGG/8Gm82Gjz/+yGcgnZ6ejk8++S+eeOJJvP/++1i1ahUyMzPRr18/3HHHgoAGaLmbPHkyXn31VTz88MO4775F6Nu3L55/fglWrnwLu3btcln2P/95Brfeas9TNpvNWLDg9ogE0vPmzcPbb6/Bc889B51Oh9LSUlx77bW49dZbfL5u1qxZaGhoxNKlr2P9+q9QWVmJl19+Ce+//35EpjYvLCzEV199hUceeRQfffQxXnnlVeTm5mLgwIFYtGiRczme57Fy5Qo899wSrFq1Ch9//DHS0tLQq1cvXHfd/G657Z5oNBo8//zzuO222/DGG2+isLAAjz/+GObOnSPp94/mOpw2bRqKi4vx1FP/xn/+8wxMJjNKSkpw6qmnuFyI/+Uv5+Paa/+Kd999D6tXrwFjDJdccjEA4N//fgrDhw/H66+/jvvvXwyZTIby8jJMnz4do0eP8vbRhBBCEghntVqkH91Fok6r1aKiohfuvvtuv0FdohszZgzy8vLxwQfvx7opKW/KlCloaWnFTz/9GOumEEIIIVGXlDnSyc5gMHR7zjHF9+mnj4l2cyLGYrF0y63esGEDfvttB8aMSZ7vSQghhJDElJSpHcnuvffew8qVK3H22edArVbjp59+wjvvvIMzzzwzoHrGieLw4cO44IILMWPGdBQXl2DPnt147bXXUVRUhKuuutL/GxBCCCGERBAF0gloyJAhEAQZnn76aXR0dKCwsBDXXTcfd999d6ybJimNRoPhw4fjzTeXobm5Genp6TjnnHNw3333BTXokxBCCCEkEihHmhBCCCGEkBBQjjQhhBBCCCEhoECaEEIIIYSQEFAgTQghhBBCSAgokCaEEEIIISQEFEgTQgghhBASAgqkCSGEEEIICQEF0oQQQgghhISAAmlCCCGEEEJCQIE0IYQQQgghIaBAmhBCCCGEkBBQIE0IIYQQQkgIKJAmhBBCCCEkBBRIE0IIIYQQEgIKpAkhhBBCCAkBBdKEEEIIIYSEgALpJPLHH79j/rXXBLTsvx58AJ999mmEW0R27tyBuXNmx7oZIdm8eRNuuH5+rJtBgnTD9fOxefOmmH3+rbfcjG++/ipmnx+Kruvsm6+/wq233BzjFgVGbG+KdROi6rln/4Olr78W62YAAA78+SemT7s41s0gcYAC6QTV2NiI6dMuhk6ncz43aNBgvPDiywG9/s677sbEiZMi1TwX9917D/7734+j8lmEJIL77r0HMy+bgdmXz8IVs2fhn/93E3788YdYNyvidu7cgenTLsbCu+9yed5iseDKeXO6ndNCsWbNajz66MNhvUciMG9Zj/b7psG8ZX3Y79V1f5w3dw7uvWch9u3bK0ErU9NPP/6Iv15zlctzq95aienTLkZjY6PzuV9//QVz58yGaLNFu4lEQhRIJyCr1RrrJpA4Q/tEdFTVmfDNH52oqjOF/V6zLp+NZctX4I03l+Pyy6/AM08/jaamRv8vTHBpaWloampEXd1h53M//7wZ2dnZMWyVNKJ1HJq3rIdu6b2AxQTd0nslCaYd++PLL7+C/v374/HHHpOgpdKxJVCwOXjIELS3t+PQoYPO53bu3IEePXri9507jj23YwcGDx4CXhBi0UyvEmldxwNZrBsQ7zo7O9DR0YHMzExkZGSG9V43XD8fEyachU2bfkJ9fT0GDBiA666/Ebm5uVi+7E388MP36OzsRF5ePqbPmIFTTjkVgP0AfOzRRzBz5iysXbsWGk02mprst/TmX/tXAMBfr70Wubm5eOzRR7D0jWUAAKvFgnfffQcbN25Ae3s7CgoKccONf0OfPn1w37334OSRIzFlynnO959x6WVY+967YAw4++yzMW36DHAch+amJjz//BIcOHAAomjDgAGVuOrqa1BYWAjAfrtNkMlgNBiwZcuvyMnJxV+vvRZDhhyHN99Yij/++AO7d1dh1VtvYdCgQbjzrrvx8Ucf4r///S90uk5kZmZi6sWXYMKEs8Jav7Gi1Wrx2muvYOeOnVAoFDjjjDMwfcalELqcHNet+8Tjum1saMCLLz6PvXv3ged59OzZA3cvvBdKpRJGgwErVizHL7/8AovFjOHDT8CVV16FdLUajY2NuPGG63Dd9TfgvXffhdFowKmnngaD0Yjrr7/B+bnvv78WO3f8hrvuvgeMMaxb9wk+/+xTaLVa9OrVG1df81f07NkTANDS0oLnlzyHPXt2o7i4BKNGj/b5vW+4fj7OOvscbN70Ew4ePIhBgwbj73+/CatWvYWNGzcgKysLN9z4N1RWDgRgDzLefedtbNy4ATqdDpWVA3HNX+37LYCAjoHZV8zBO2+vgclkwplnTsDls6+QdFt6s3m/Hk+sa8YfXQLoQSVK3Dw5HyP7pIf13hzHYcSJJ0KtTsfhw4dRUFCIb77+Cv/973/x2ONPOJe79ZabMWXKFIwbf6bP/QYA6urqcNedC1BbW4vevfvgb3+/Cfn5+R4/P9z1/um6T/DBB+/DZDLh7LPPCej7nnHGWHz99deYOXMWAOCbr7/GuPFnYsXyZc7lrFYr1qxehQ0bNsBsNuO4447DVVddjazsbDDGsGLFcnz7zTcwm03QaDS4Ys5c2Gw2rH3vPTAmYvbl9vdetnyF333fl2CPw5dfiWzqgTOIFo8GO6LN/hiAYsSEsN9fJpdj7Lhx+PDDD3Ckvd25vn2tP1/n8+3btmHVqrdw+PAhKBQKTD53Ci66aCoA4LvvvsXa995FW1sbysrKMe/Kq9CnTx8A9l7yfv364cCBA6iq2oWb/vF/3dr6zDNPY+eO32A0GlFcXILZV1yB444b6vF7hfM7BgA6nQ4vvvg8tm/bBo1Gg3N83NHNyspCWVkZdu6wB88mkwm1tbW44oo52LFzJ8aNPxMAsHPnTpwxdmxA32Xduk/w4dHj7JxzJuLXX391ng8c54tRo0fj03WfgOM4XHjRVEyZcp7z9d9/vxFr33sXzc3NKCkpwdx5VzrPzZ7W9UknnexvVyFHUSDthdlsxs+bf0R9/bFek+LiUpw88hQoFIqQ33f9+vW48667kJ9fgJdffgn/eeZp3HvfIlT06oXz/3IBMjMy8ONPP+LZ/zyDvn36orCoCABgMBhxoLoa/376GQBAe3s7brzhOrzw4ktQq9UA7D96Xa1YsRx//PEH7rzrbhQXl6Du8GHIFXKP7TIYjPhz/37859klaG5uwgOL70dhURHGjRsPkTGcd/75GDLkOFitVrzw/BK8+MLzWHjPvc7X//jD97jt9gX4+99vwtr312LJc8/iuSUv4Io5c7F//35n0A4Ahw8fxqpVb+GRRx9Djx49odVq0d6uDXmdOtTW1uLgwYP+FwxDz549UVZW5vLcM0//GxqNBs89twQdnR146F8PQqlSYepUe/6cr3X71lsrUVRcgjvuvBsAsG/fXgi8/UbRkuefg8ALePyJJyEIAl54YQleffUV/O3vNzk/+5eff8bDjzwKmUyGgwdrcf+iRbj6qquhOBpQbfjuW1x49Efr888/w9dfrcftC+5AYWERPvvsUzzy8EN46ql/QyaX45mnn0JBYRFeevlVNDc34V8PPuh3ffzw/fe4fcEdSE9Lw8KFd+Guu+7AzFmX48orr8I777yNl196EY8/8RQA+63N/fv34/7FDyIzIwMr31qJp//9JBbd/wAABHQMHKytxdPPPIvGxkbcseA2nDBihPOHLlI279dj/tJDsIquz/9RZ8L8pYfwwtweYQXToiji119/gdlsRq9evQN6ja/9BrBv99tuWwBNTg6eePxRrF71Fm648W8e3yuc9b7jt9/w1ltv4a677kafPn3w9ttrUFtb47f948aNx/3334dLZ1wKrVaLffv2Ye7ceS6B9Ptr38Ovv/6KxYsfQEZGBl544Xk888zTuHvhPdi+fRu+37gBjzz6GHJzc9Hc1ASzxYLS0lJcNHUqDhz4E7fdtsD5Xv72fV+CPQ4jqVsQ7SBhMG02mfDV+vXIzMyCOiMDgO/119jU5PV8/uef+/Hoow/jxr/9HSeddDLMJhMOHjoEAPj995145eWXsOCOOzFgQCU+++xT/OvBxXjmmWeRfvQ37ZtvvsaCBXeib79+sJjN2PTTTy5tHTp0qPN898l/P8aTTzyO55a8gLS0tG7fK5zfMQB4/bVXodfp8dySF2AymfDoI77Th4YMOQ47d+7EORMnoWrXLvTt2xfDhh2Pd955GwCg1+lw4MABXHe048PXd/ntt+1Ys3oV7rprIXr16oV3330HBw/WunzewYO1OEM5Fi+8+DKqqnbhgcX348QTT0JxcTG2bPkVy958A7fdfgd69eqFn3/ejEcefghPP/MsMjMzPa5rEjhK7fDCPYgGgPr6w/j55x/Det9zJk5Ejx49oVQqcfnls7Fz5w60tLTg9NPPQHZ2NnhBwGmnjUFpaQ9U7a5yvo4xEbNmXQ6lUunsdfKFMYYvv/wCV8yZg5KSUnAch9IePVBQUOhleRGzLp8NpVKJHj16YuKkydjw3bcAgMLCQpxwwggoFAqkp6dj6tSLsWvXHxDFY5HFCSfYf1h5QcD4cePR1NSEjo4Oj5/F8zwYswe+ZpO9N6miolcQa9Gz5cuXY9KkyRH9t3z5cpfPbG1pwY4dv+GKOXOhSktDQUEhpk69BN9+83VA61aQyaBta0NTUyNkMhkqKwdCJpfjSHs7Nv20CVddfQ3UajVUKhVmzLgUP/zwg0s+3bRp06FWq6FUKtG3bz/k5eXi519+BmD/EWtubsbIkaMAAJ99ug7TZ1yKkpJSCIKAc8+dArPZjD1796C5uRl//PEHZs++wtnOs8/x37t4zsSJyM/PR7pajRNOGIGMjEyMGjUavCDg1FNPQ21tLawWCxhj+Oyzz3DFnLnIycmBTC7HpZdehl27qtDc3AwAfo8BgOHSy2ZCoVCgZ8+eGDCgEvv37w9pXwnGE+uauwXRDlbR/vdQrFyxAnPnzMbsy2fh8ccew9SLLwk4vcHbfuNwzsRJKCwqgkKhwJjTz/C5nsJZ7xs2fIfTTz8dAyorIZPLMW36DCiVKr/tt5+LCrBt+zZ88+03OPXUU7sFtN999y0uvvgS5BcUQJWWhivmzMX27dvQ2toKQRBgNltwsLYWVqsV+QUFKC0t9fp5vvZ9X0I5DiPFaxDtcDSYDjXNw7k/zp6FjRs34pZbb3XeVfO1/nydz7/84gucetoYjB59CmQyGdLVagwYMAAA8N133+H008/A4MFDIJPJMGXKeVCrM7BlyxZnm8aMOR39+vcHx3HOzoGuxo8/E+lqNWQyGf5ywYVgjKG6utrj9wvnd0y02fDDD9/j0ksvg1qtRm5uLv7ylwt8rs8hxx2H33/fCcDeyTVo8BDkFxSA43g0NNTj9z9+R0aGGhUVFX6/y8YNG5zrQiaX4+JLpnXb1zIzM3H++X+BTCbDkCHHoaCgAAcO/Hl0+32Kv/zlAvTp0wc8z2PUqNHo0aMH/rfl14DXNfGOeqQ96Ozs6BZEO9TXHUZnZ0fIaR4F+QXO/9doNJDL5WhtbcGPP/6Ar9Z/iZaWVnAcYDQa0XHkiHPZtLQ0Z89zII4cOQKTyYSSEu8/Ll3J5QqXH/GCggK0trba36u9Ha+//hp27foDer0egH1wkNFgcPYcaDQa52uVKvsPqcFgcF7tdlVcXIwbbrwRn326Ds8veQ79+w/A5ZfPRq/egfXGxZOW1hbI5QqX719YVISWlhbnY1/rdvbsK/D2mtVYfP8icByHsePG45JLpqGxqQmMibjxhutcPo/nOWi1Wudj99v1Z5wxFt99+w1OO20Mvv32W4waNdp5wm1qasJ/nnkafJeeS6vVipaWFshkMo/t9Kfr8gqlEtka18eMMZjMZtj0ephMRtx7z0Jw3LHXy2QytLQ0Iz8/Hx9//JHfY6Drj4dSpYTBYPDbxnBU1Zlc0jk8+aPOhN31JgwoDu7HZ+asWc67NPV1dXjkkYehVqsDSo/wtt84tq3L8ahUwmj0vp7CWe9tbW0YPGSI828ymQw5Occ+25dx48/E1199herqA7jppu637VtaWlFQeGwfzM3NtZ8vW1pw3HFDMX3GDKxe/RYOHjyEocOG4orZc5y96O587fu+hHocSs1vEO0QRs+0Y39sbWnBI488jOrqagwaNBiA7/U3aNBgr+fz5uYmDDz6Hu5aW1pc9h3AHuy2tB7bJr7WqyiKWL3qLfz44w9ob28Hx3EwGAzo6DjicflwfsdsNpvzgs3ZNj/nx8GDh+DIkQ7U1tbg9993Yuasy48+Pxg7d+zAwYMHMXjwEHAc5/e7eD7Oclw+Lztb4/JYqVTBePQ4bWpqxFtvrcSaNaudf7fZbM7fISDy+3Ayo0DaA289qQ6dHaEH0k3Nx8oVtbe3w2KxwGa14e01a3DvvfehV+/e4Hket95yM1iX13Fdow8AvNtjd1lZWVAqlaivr+t2wHlisZjR3t7uDIyam5uduasrV66AyWzCI488hqzsbBz480/cdtstLu3zxb3tAHDqqafh1FNPg9lkwurVq/Cf/zyDJ558KsB3jB95uXmwWMzQarXOk3BTUyPy8vKcy/hat9nZ2bj6Gnuee011NRYvvh/l5eUYWDkQHMfjxZde8djL5Rj5zfGuN5XGnH4G1qxZjdaWFny/cSNuuukfx9qal4e5c6/E8BNO6PZ+zc3NHtsplYzMTCiVSvzroYfQo0f3vNRdf/zh9xiIhTqtJeDlgg2kuyouKcEJI0bg119/wdlnnwOVKg1ms2sA3zVw87bfjB59SlCfG+56z8nJQXPTsXOa1WpFW5vW+wu6OPXU0/DG0tdRVFSEPn37ulQzAIC8vFw0NTahf397D6a2rQ0WiwW5R4+tiRMnYeLESdDrdHj55Zfw2uuvYsGCOz2eG33t+77k5+WFdBxKSWxvgm7ZYv9BtPMFNuiWLYas7zDw2f4vht3l5uXh2vnX4d57FmLkyFHIzc31u/68nc/z8wtQX1fn9XOa3LZ5U1Mj8nKPnTt9rdeNGzdg48aNuOvuhSgpKQHHcZg39wow5nnvDed3LCszE4IgQ3NTk/M839zsu/RgRkYGevWqwJYtW1BTU4t+/foDsA9E3PHbbzh48CDGn3lmQN8lJycHLV3OxzabDW1tbQG03C4vLx+TJp+Lc86Z6HWZSO7DyY7WnAeeelG7yvDzd1++/OJzHD50CGaTCSuWL8OgQYOhN+jB8zyysrLAGMNXX633m2eYlZVlv0VUX+/x7xzHYcKEs/DmG2+gvq4OjDEcPnTIa1UAjuOxcsVymE0mHD50CJ99ug5jTj8DAKA36KFUKJGuVqOjowNvv70mqO+crcl2aefhQ4ewfds2mE0myGQyqNLSXAbmheryyy/Hp5+ui+i/yy+/3OUzc/PyMGTIcVj25hswGo1obmrCe+++i7FjxzmX8bVuf/jhezQ3NYExhnS1GjzPQxAEaHJycPLIk/Hqq6/gyNHeQW1bGzZv8l0fOD8/HwMHDsLzzy+x3+I77lj+8MRJk7F69SocPpqjqNfr8fPPm2EwGJCfn4/KyoFY0aWdX37xedjbxIHneZx99jl48403nAF6R0cHfvj+e3tbQjgGoqFE4zt/NtjlvGlsbMT/tmxBebn9Nm+vXr3Q0NCAP/74HTabDR988D46O49d4Hvbb4IV7no/bcwYbNiwAXv27IbVYsE777wNk8kY0GvT0tJw732L8H//57lm8+lnjMXatfbBUUaDAW+8sRRDhw5Dbm4u9u7di6qqXbBaLFAoFFCqVBB4+/fP1mjQ3NTsUnnA177vS6jHoZT47AKoZy8E+AC3Ly9APXthSEG0Q58+fTBkyBCsfe9dAL7Xn6/z+YSzzsb332/E5k2bYLPZoNfpsHv3bgDAGaefgY0bN2DXrl2w2WxYt+4TdHR04oQRIwJqo8FggEwmQ1ZmJqxWK955e43P7RnO7xgvCDjl1FOxevUq6HQ6tLa24qMPP/D7uiFDjsN/P/4Iffr0gfxo6tLgwYOxbdtWHDhwwDmY0N93OW3MGGzcuBH79u21D9p+9x2YTIFXDpo4aRI++vAD7N+3z36X0GTC9u3b/N6RIYGhHmkPMjIyUVxc6jG9o7ikNKzqHePHn4mnn34K9fX16N9/AP5+0z+Qk5OD0aNH4+ab/wm5XIYzzhjrHE3rjUKpxLRp0/Cvfz0Aq9WKq6/+K3JyXXueZ10+2377d/EidHR0oLDQXrXDU550WpoKvXr3xo03Xg/GGCacdbYzGJw+/VI89+x/MG/uHOTl5eK8887Hzz9vDvg7T5lyHpY89yzmzpmNgQMHYebMWUdvyR4Ex3GoqOiF62+4MeD386asrKzbQMBouOmmf+DVV1/BDdfPt+ekjjkdf7ngQufffa3b/fv34803lkKn00GtVuPMM890jpa+4Ya/Yc3qVbhjwe3o7OxAdnY2Tj31NIwcNcpne84YOxZLnnsWF1401eVuwKRJk8HzPB5//FG0tLRApUrDwIEDnSfzm276B55/fgmuvvpKlJSUYvz4M7F+/ZeSraeZM2fhgw8/wP2L7oVWq0VmZiaOO24oTj3tNAwffkLQx0A0VJYoMahE6TO9Y1CpMqTe6BXLl2HVW28BANTqdIwcNRqXXDINgL2H+vLLZ+PJJx6HKDJMPvdc9Ox5bN/2td8EI9z1PmzY8Zhx6aV44vHHYDabcfbZ56CsrDzg1/ft28/r3y668CKYjEbcfdcdsFgsGDLkOOcAP4NBjzffeAMNDfUQBBkGDBiAa4720J8y+hRs3LABV181D4wxLH1jmd9935dQj0MpOdI0/KZ38ALUcxdJUr1j6tSLsWjRvbjgwot8rj+r1er1fN6nTx/cfMutWL1qFZ577j9QqVSYfO4UDBgwAIOHDMG8K6/CC88/h7Y2LcrLy3DnXXcFnMI4duw4/LZ9O66/fj7S0tIxZcoUlzuB7sL9Hbvyyqvw4gvP4/rr5iMnx161Y9++fT5fM+S44/Dxxx/h7C49wUVFxZDL5cjOznJWPfH3XYYNOx7Tpk3DY48+ArPZjHPOmYiSkhK/A2UdTjrpZFgsFrz44vNoaGiAXC5H3379cPVVgU3gRnzjrFZLrO+exiWz2Yyff/4R9V1qnRaXlOLkk0Ov2nHD9fMxZ+485+CveOEoceUom0cIOcZb1Q4AkPEIu2oHIYHymSstYRBN4pvVYsGVV87FnXctxMCBse9wSHXUI+2FQqHAaaeNRWdnhz0nWoI60oSQxDOyTzpemNsjYnWkCQmU155pCqKT3qZNP+GE4SdAZAyrVr2FzMxM9OvbN9bNIqBA2q+MDAqgCUl1I/ukY/UN5dhdb0Kd1oISjTyswYWEhKpbME1BdEr47ttv8fyS58CYfQzFbbcvCDi1g0QWpXYQQgghCca8ZT10yxZDPXshBdGExBAF0oQQQkgCEtubwqrOQQgJHwXShBBCCCGEhIDqSBNCCCGEEBICCqQJIYQQQggJAQXShBBCCCGEhIACaUIIIYQQQkJAgTQhhBBCCCEhoECaEEIIIYSQEFAgTQghhBBCSAgokCaEEEIIISQEFEgTQgghhBASAgqkCSGEEEIICQEF0oQQQgghhISAAmlCCCGEEEJCQIE0IYQQQgghIaBAmhBCCCGEkBBQIE0IIYQQQkgIZLH88N9/34kPP/wAf+7fj7a2Ntxy620YOXKU8++MMaxZvQrr138JnU6PgQMrcfU1f0VJSWkMW00IIYQQQkiMe6RNJhN6VfTCVVdd4/HvH3zwPtat+wTX/PVa/Ouhh6BUqvDgA4thNpuj3NLUIbY3xboJhBBCCCEJIaaB9AknjMCll83EyFGjuv2NMYZP/vsxpl58CU4+eSQqKnrhxhv/hra2Nvz88+YYtDb5mbesR/t902Desj7WTSGEEEIIiXsxTe3wpbGxAVqtFsOGDnM+l65Wo1+//thdVYXTThsTw9bFhiiK+PGHDWhoqANjDDzPQ63OAGMMOl0nGGNur+DA8xxEUXQ+I5crcO6UCyCTuW5685b10C29FxBt9v8CUIyYEOmvlFIYY9i7dzdaWpqQl1eAfv0GAEC35ziOi3FLCSEkdJ7OdXRei09dt1Vubj6aGhvQ0tIEQZChf/9KiIyhatfvsFotLttQoVBApUpDRUVv9OtfmdLbN24Daa1WCwDI1mhcns/WZDv/5knXoFEq3QNU6fA873MHZIw5v9Omn75Hff1h599EUURHxxEf784giq5tt1jM+PWXTRg56lTn53YNou1v3D2Y7toOqaTagbd3TxV++20rAODQwVqwo+vT5TnGnAE2IYQkiq6/k/v27u52rusb4nktmX8nIhlbBMp9WzlYLBbn8w5d22symWAymbB9+//AGAt5+4YrkvsHzweWtBG3gXSompvqJX9PmUzWrQc3XIIgQ7o602UnqKozoU5rQYlGjsoSJQD7TsLzPPS6DrS2Nof9uUVFJTjp5NHeg2gHt2Da2Q59J2xWS9jtsFgssNls/hdMIl0vggCgsbEegOuJtO7wQWiys6LYKkIICZ9SqXT+rjQ21rn8rbGxHsXFRUG/pyiKST0mqus6ixX3bRXae4S2fcMV6f2jsCiwwhZxG0hrjvZEt2u1yMnJcT7frm1Hr169vL4uv6BY8rYwxmAxm2A2m+Ae+IQqLT3DeQBt3q/HE+ua8Uedyfn3QSVK3Dw5HyP7pIPjOKSlZyAnJxf19aHv9EVFJTjl1NMhCAIAH0G0g4dgOj09I6xgmuM4yBUqqDOyY34Cibb8lhY0NTU6Hzv28YaGBudzJaU9I7IPE0JIpIiiCL3u2N1RjSbH5bymcbuz3BVjDLW1NdBq26DR5KCsrNz528DzPPLyi5Lyt8JiMcNk1Me6Gd22VWjvoZGmMUHieR6anDzIZPKYfL5D3AbShYVF0Gg0+G3Hb+jVuzcAQK/XY+/ePThn4kSvrwu0Kz5YQlo6lKo0WK0WWCxm2KyWkG/LyBVKZzs379dj/tJDsLplTfxRZ8L8pYfwwtweGNknHTzP44QTTsTPP/+E5mZ7zzTP81CpVDCbzbBarT4/MycnN7gg2sFLMK3THYEYYI8yx3GQyeSQyRWQyeRJeVIMRO/efSDarNBqtdBoNCgrKwcAyGRytLe3Izcvn3IJCSEJx2Zz/f1xnNvcz3We1NbWoKpqF4BjnQrl5RXOvzPGnL9bycRqiY+e9q7bKisrG21tLWhr00ImE1BeXgFRZDhwYD9sNpvzt4kxBoVCAYVCidLSUp/bN9KsFjMUCmXMPh+IcSBtNBhQX38sFaOxsREH/vwTGRkZyC8owLlTzsN7776DkuISFBYWYtXqt5CTk4OTTx4Zk/ZyHAe5XAG5XHE0Z9gGq9UCm9WKmppqHDx4KKD3GTX6FOf/P7GuuVsQ7WAV7X9ffYN9J7XZGCwWEdnZuS7L1dQcRGtbq8/PPPvsic6TkdjeBN2yxf6DaAfRBt2yxZD1HQY+uwAcx6G1VYvdR09+7jiOQ1l5BXpVVECQycDzAgWHsK+X8vIKlx8JAOjTpx8USlWMWkUIIeER3X5LvJ3rPNFq29wea11eJ4pWxHGfX0hEUex28REr7tvK0x3/Pn36RLlVgbMe7dSMZYwR071z3/59WHTfvc7Hb76xFAAwduw43HDj33DBBRfCZDTixRdfgF6vw8CBA3HnXQuhUChi1OJjOI6DIMggCDJACby39iM88sgjfl+XlpbmTM+oqjO5pHN48kedCbvrTRhQrERmVhYunz0XBoMh6PZ+/vnXWLfuE6hUKvDZBVDPXhhYjzQA8ALUsxeCzy4AABiNRsyZMw9btmzx+pIFC27HHXfcEXQ7CSGEJJZwxrv4SwMJ9M5nIrFKMM6IHGO1WiCXxy4ujGkgPWTIcVjz9rte/85xHGZcehlmXHpZFFsVmkAvhjIyMpz/X6cN7GCq01owoFjpfH0ogfSWLVswc+YsrFy5AiqVylmNw28wzQtQz13kXN5oNGLmzFk+g2hCCCGpI5xg118aSDIOSo+X3uhkYbNaYxpIx3RCllTU2dnp/P8STWAJ8l2X6/r6YK1fvx4zZ86C0WgEYM95Vs9dBPBe8s+8BNHr19OELYQQQuBMcwyVI7Vg2LDjUV5e0e0WvU20xUWZOClRIC2tWK9Pzmq1JNceGiO1tbU4ePBgQMsOGzYMarUaADDjuRqf6R2DSpVYfb39Cr3jyBFs2vRj2G3Ny8vHsOOH+x546BZEi6KIqqoqnzW8u+rZsyfKysrCbmsyMRr1MJuM3Z5XqdIpR5oQkpBsNit0nb7mMwhfRqYmYoUEoo0xho4jWkhVASzSfFVViRccxyMzSxO7z6dAOra8Ve0AABkPZ9UOADDodbBYfOdUB0qQyZHepQSfSzDtFkQTaVAgTQhJNhazCQaDLqKfkZaeEdNb91KyB9Jt/heMEzU11c6qKgBQWTkwoEGk0ZaZlROzAD85LvES2Mg+6Xhhbg8MKnEt3zKoROkSRDPGYJGgXA7HcVAoVFAp01yed6Z5yJUUREdbfF3cE0JIwGxuaR2MMdTUVGP79q2oqamWJC0jnNSReJNoaSruVVUOHDgg2XZNFtQjHUd21x+b2dAxuLArRzBtNhmCnq6b5wUolSrI5AqfV21ie5OzOgeRlslogMnUfaCoKk0d8zqYhBASCr2uw6UKRSR6MOVyBdLSM/wvmABE0YbOjvZYNyNg7tvTId56pjOzNOC42PQNJ1dxxgQ3oFjpMYB2sPcmKyGXK2Axm2AMYFYkjuOgUqX7DaAdKIgmhBASKPeqGv7qQkvxGYktsW5BOqqoHDhwAKYuqYlSbFdpxW69UmpHAuI4DgqlfZptgPN6G43nBagzsiFXKONucEBK8rIJaMsQQhIRYwyMud4d1Why3B5rwv4cMYkqdyTab7Gjqor7RC2xmhY8HlGPdAITBAH1DfUep1flOB7p6sykGelMCCEkvngqOxbM9ODBEEUxKaYKtwfSHBKlaodDpLarFDiOT92ZDUn4WltbXB47breo0tIpiI473rqkE6uHghBCAM+DAIOZHjzYz0qGQBoAeJ5PuAGUkdquUoh1rEORVoLLy3PNadZoNOA4DjJZYJO9kOihcJkQkkyiOX13Mk0VHs8XBJGouhJpsV6f1COd4Pr1GwAAaGyoQ3Z2NsrKyiHI5AmXh5USaJMQQpJINAcBupfZS2S8IAMkKGcbCbW1NR7TReMZL8Q2lKVAOsFxHIf+/SvRs0cPZ53pWN/mIN54jqQ5irAJIQkomukJYhJNqy0TZJBmajXpRaLqSqQJMQ6kKeJKEl17oKk3Oj553Sq0uQghCUYUxaje9o/250USH8epHZGouhJJHMfFvPOQeqSTRNcTTLKcbAghhMSnWAyWsw84TPywxTGOqetENvEinqtzeCKLg1TWxN8jCQC4zHTIgpz1kESJl4OdUjsIIYkmFoP/RFtyBNIAIJPHZyAdz9U5PImHwgqU2pEEGGMu9Tw91fYkcYziaEJIgonF4L9kGnAokyli3YSkIJPHfj1SIJ0E3ANnURSTbErV5OD99hNF0oSQxBKrHulkwfN80vSux0o8pHUAFEgnBYu5+/hfiyVexwQTd7E/DRBCSOAYYzHLkU4m8jjoTU1k8bL+KJBOcIyJzrJ3XVnMJhp0GHe8zWwY3VYQQkg4GGMx+X1JpsodQHykJSSyeFl/FEgnOLPZc1F3xlhcDmRIZd7vQFEkTQhJHLHsGU6mXmme5+NisFwikssVcZHWAVAgnfA8pXUE8jcSP+LkXEAIIQGJ5RicZBv/I1coY92EhBRP640y3RMQYwx79+5GS3Mj1Go1AKC9XQuNJgdlZeXgOA6MMezfvxednTrk5RegX78BcXP1lqq8l7mj7UIISRxSDvpjjKG2tgZabZvLb5jXz06iHmng2IC5WKWsMMZQU1ONurrDAICSklKUl1d43QbBbq9I4Lj4GqgZPy0hAdu7dze2b9vS7fmGhgYAQHl5BWpra1BVtQsAcOhQLQCgf//K6DWSdEfxMiEkCUgZzHb9rer6G+b1s5OsR5rjOMjlCphjdAe5trYGu3dXOR93dFQ5a0l7Wz6Y7RUJckX8pHUAlNqRkFpamrz+TavVHv1vm8vzrS3NkWwSCYiXCVni6IRACCG+uM9bEC733yrHb5g3yZbaAcQ2TcF9/duf0wI41lu9fftW1NRUgzEW9PaKBIU8ftI6AAqkE1JeXoHXv2k0mqP/zXF5PjcvP5JNIgGgeJkQkugYk3bmXPffKsdvmK/PT6bKHQDA8wJ4XojJZ7uvf/tzGgDHep8bGhpQVbULtbU1QW8vqQmCDLwQm3XlDaV2JKB+/QYAAJoa65GZmQkAaG9vh0ajQVlZOQA4/3ukowMFBUXO15BY8hRJU3RNCEkcUvcIO36rtFqty2+Y7zZYk6raBcdxkCsUMBkNUf/ssrJyMMZccqSPbZPuvc9Dhw5z/n+g20tK8TTI0IGzWi3JdWmXQkwmg98DLy1NHZc7XipijKHjiOuJieM4ZGZ17xEghJB4FMjvTqSp0tKhUKhi2gapiaKIzg5tVD/T38DBmppqZz40AFRWDox6PrS7zKycuEuHpB7pBBbIraB4uwVCCCEkccXDYL94aIPUeJ6HIJPDFsX5H/wNHAzlbkEkyeKodnRXFEgnMCGAIDlWeVekO88ngPg7KRBCiDfxMNgvHtoQCXK5IqqBtKfUja6BtKN6R6x7oR3iZUpwdxRIJxjGGPbsqUJN9Z8AgOLiYq91HAVBiMurN3IMbR5CSKJgjMVFHWdRtIExlnS/b3K5AkaDLmqfp9HkOHui7Y81Ufvs4HFxmxdPgXSC2bt3N37b/j/n4/Z2LQDPdRyFON3pUln3wvvJ9UNACEle8RBEA/aAPhkDaY6zB4vWKPVKRzJ1gxcEyOVKWCwmSVJx5HJ53G5vCqQTjKca0p2dnVBnZHV7nueouiEhhBBpxFNuss1mBc/H563+cMjliqgF0pFK3RBkcqSnZ4DjOCgUSuj1nWGnrMjiNK0DoDrSCcdTDem8/EIIgqzbP46nzRt/XK+o4/QCmxBCugk3N9nTBB+hiqegXkrRDBil3B4OXYNowB6sp6dnhH2HPF7TOgDqkU44/foNsO/8R3OkKyp6U43oRMIBYO5PEEJI/LOJ4c1oKOX00rY4STORGsdxEGQy2KzSzR7pjdTTfbsH0WJ7E/jsAmcwHWrPtEwWv2kdAPVIJxyO4zBgwECcdfZknHX2ZPQfMDCudzDiirYUISRRhdsLLOX00qKE05THG5ksOr3SUm4P9yDavGU92u+bBvOW9QDC65mO57QOgAJpQqLMLZSmiyBCSAIQxfCn5pZyemkp2hOvopXGINX28BRE65beC1hM0C29N+xgOp7TOgBK7SAkpiiMJoQkApsEPcBSV4kQbTYIsuQLY3ie91DhSXpSbA+vQbQj9Ua02R8DUIyYEHSaB8/z4ON8vBdNEU5IFHV2trvcHpXJ5EhXZ8awRYQQ4p/JaIDJFNupwd2pVOlQKJNrqnAAqK2txd69e6I6OUso8vLyMez44c7J4boF0V3xAtRzF0ExYgIA+8DV7du2oqWl2ednyOQKKJUq9OzZE2VlZZJ/Bykk36UcIXGMo9QOQkgCkqJHWmrJOuBw+fLlePjhR2LdDJ8mTJiAlStXBBZEA916pgVBwKDBQzBz5iysX7/e7+ctWHA77rjjDsnaL6X47i8nhBBCSMzF47Tc8Rjcp4IRI0Zg5coVUKnsdwP8BtEOR4NpR860SqXCypUrMGLEiEg3OaIokCYkhqg/mhAS7xgTwZgY62Z0I9psSTvgMJ7NmDHDGUSL7U3QLVvsP4h2EG3QLVsMsd0+uZxKpcL06dMj1dSooBxpQqJIpzviUh9ULlcgLT0jhi0ihBDfrFYL9LqOWDfDI3VGFgQhubJUa2trcfDgQRgMuriceCYrKxsnnnRy4GkdXXnIlf71l59x5Eh7t0UFQQZVWjoAxHWONAXShESRTtfhMoCEAmlCSLyLx4GGDqq0dCgUyTfgEAAMBh0sZlOsm+GR32odnrgF0Ywxn9U7FEoVVKp0ydsuNUrtICSKuqdyUHIHISS+xXMuss0afz22UhF4IdZN8MpmtUCv73Sm1ihGTIB67iLAW5uDDKKB+P7+XSXX/RBCEg3F0YSQOBePAw0d4jnIDxcvxEcgyRhDbW0NtNo2aDQ5KCsrB8dxzmDa0TPtCJK79UyHEETbXxYf398fCqQJIYQQ4pF9BsH4G2joIIr2AYdcEpYS5eOkR7a2tgZVVbsAAA0NDQCA8vIKAPAfTIcYRAPx8/39odQOQgghhHiUCD2+idDGUNgvDiJ/gcAYQ01NNbZv34qamupulVC02ja3x1qXx17TPOTKkINojuMT5uIornukRZsNa95egw3ffQetVovc3ByMHTceF198ScKsYEJ8ob2YEBLPEiFItdmskMnksW6G5DiOA8/zECM88YyvHmcA0GhynM/bH2u6vYennmlZ32HgswsABBdEA4j7acG7iutA+v0P3scXn3+GG274G3qWlWH/vn1YsuRZpKen49xzp8S6eYQQQkhS61quM17ZrFZAGetWREY0AmlPPc5dA+mysnLn8xqNxvnYnXswHWoQDVAgLZndVVU46aSTMeLEEwEAhYWF2Pj9BuzduzfGLSMkRHQnhRCSIBhjCdMjnbx50pEPKLOzNS49ztnZ2R4HGHYNrr1xD6ZDCaIBgKNAWhoDKiux/ssvcPjwYZSWluLAgQOo2rULV8yZ6/U1ohi/gyIIgVvuGWOM9llCSFxKhCAaOBbwJ8rgtKDE6OLAX7qHLzarBTrdEcjlSlgsppAmleHAxfy3MdCLmLgOpC+88CIY9Hr83z/+fvT2hohLL5uJ008/w+trmpvqo9hCQoIjl8uds0EBgMGgR0fHkRi2iBBCPBMEAXJ5YuQea9taYh54RQLP81AoFBH9jPZ2rdvjdgDuAw61AQfSgH2Mm8mmD7lNHR3tMd+ehUWlAS0X14H0jz/+gI0bN+DvN/0DZT3LcODAn1i69HXk5ORg3LjxHl+TX1Ac5VYSEjiDQedyiystLR3KBJi5iRCSeowGHaxB3pIPlrcaxcFKV2ckxCx4wbJZrTAYOiPy3o51r9e7BryOwYSeBhhKtb38ydbkuXQ6xbO4DqSXL3sTF1x4EU47bQwAoLyiAk3NzXh/7XteA+lESlAnqcf9dOMYlU0IIfEkWvnR4aQQdCXabMl5Lo1gMNl13QNAZmYmSkt7uAwmdB9gKNX28kcQhITZnnHdSpPJBN7tSofn+W41DglJVLQnE0LiEWNiVH5r/dUoDpQo2sCSMLUjkoPu3Nd9eroa5eUV4DgOHMehvLwCw4Yd73zO02tC3V7+JNLA0bgOpE888SS899672PLrr2hsbMTmTZvw8Ucf4eSRo2LdNEIIISRpWaNU9k6jyXF7rAn5vawJMjgyGJEMKENZ91JuL28SaTIWAOCsVkvcdooZDAasXvUWNm/ehPb2I8jNzcFpp43BJZdMgyxBBkAQ0pVe1+GScyhXKJGWpo5hiwghpDuDvhMWiznin3Ms5/ZYCkGoQZRCoYIqLfnypDs7tBEZeBfKupdye3nDCwIyMrIlfc9IiutAmpBkQ4E0ISQRdHRoEy5VItECsEDpOo8kTClCKchkcqSrM2PdjIDFdWoHIUmPLmMJIXEmUfONRZsNjCVeu/1JpMlJpJBo3zexWktIguseN1MkTQiJL9HKj46ERG67NzyXWqFaon3fxGotIQmPAmdCSHwLdjrneJLIbfcm0Xpow5Vo3zeu60gTknQojiaExDHGmGS9utGavKOrpOyRTrDAMlyJ9n0pkCYkhiiuJoTEE1EUJcszjtbkHV2Jog2iKCZcMOZLMn2XQCTa902s1hKS8NxCZ5pciBASR6RMjYjW5B3uki29I9ECy3BxlCNNCPGGwmZCSDyzShiERmPyDk+k/A7xINEmKAkHzyfed6XUDkKiyb1DOjatIISQbqTMjwaAsrJyAHCZvCMarFYrGGMJF5D5wvNCStSS5nkh1k0IGgXShEQVpXYQQuKTPVCT7pzEcRzKyysinhftjjERomiDICRPiMMLFEjHK0rtICSKKGwmhMSrZEqJSKbvAgBCAgaYoRCExPueFEgTEk3dImkKrQkh8SGZgs9k+i6AvUc6FSTi96RAmpCocg2cKbODEBIPRFGEaLPFuhmSsR3Nk04WyZSm4guldhBCfOp+Xk+eEz0hJHElWw8ukFzfieO4hOytDYYgyBJygCgF0oREFfPxiBBCYiOZgk6HZPtOsiTvlZbJ5LFuQkiC3iq33XoLPF8wcFAo5CgqLsG4ceNw3HFDw28dIUnE423GJLr1SAhJTIwxWC3JFXQCgNViAVMlTxk8QSYHzKaofBZjDDU11airOwwAKCkpRXl5RUTXpSBLzAuFoHukhw8fjoaGBiiVKgwZchyGDDkOKpUKDQ316Nu3H7RtbVh8//34+efNkWgvIYQQQiQkddm7eOEog5csotljW1tbg927q9DR0YGOjg7s3l2F2tqaiH5mouaBB93qIx0dOO/8v+CSS6a5PP/uu++gqakJdy+8B2tWr8K777yDk08eKVlDCUl83X+okmkwDCEkMVkt5lg3IWKsFkvCBmjuOI6DIMiiUk/afXp3+3PaiNUEl8nkCXvnIOge6R9/+AFjThvT7fnTTj0NP/7wg/3/x4zB4cOHwm8dIUmAMYY9e6qw6acfUFNT3S14pmCaEBIrjDFYYpBL7Egd2L59q8fzolQs1uS6SJDLFVH5HPfp3e3PaSL2ebIofa9ICPoyTaGQo2p3FYpLSlyer9pdBYXCftuBiSxqG5uQeLd3725s37bF5bloz/RFCCGeiKINTBSj/rm1tTWoqtoFAGhoaAAQmfOiaLNBFEXwfHLUVpDJFYBRH/HPKSsrB2PMJUc6klO8y+WJOdAQCCGQnjT5XLz80ovYv28f+vbrBwDYt3cvvvpqPS66aCoAYOu2rejVq5ekDSUkUbW0NLk87n57jAFIzFtahJDEFqtBhu6pA5FMG7BazFAoVRF572jjeT4q6R0cx6GiohcqKnpF9HMAR1pH4l7ocFarJej7KRs2fIdP163D4cP2K5XS0lJMnnwuxpx+OgDAbDIBHAeFgnqlCdmzp8qlR7qycqDLD0ZGpiZpeksIIYmls7M96hOxMMawbdtWNDU1Op/LyMhEjx49UFZWLnmurCCTQ63OlPQ9Y8lsNsFo0IX1Howx1NbWoK2t9egzHHJycpzr3/F3rbYNWVkaaLWtOHLkCLKzszF06PGS/malpWckdBZDSIE0ISRwjDHs3bsbzc2NyMzI6PZDoc7IhpDkhfYJIfFHFEV0dmij/rk1NdXOtA537h0NUsnM0iR0r2dXjDF0HOk+GDAY3raBY/372kaFhYU4/vgTwvp8B47jkJGpSdiBhkAIqR0OVosF7UeOdMutyi8oCLtRhCQTjuPQv38levXqDYO+08MSdC1LCIm+WFXr8FQR4tjfIpPiYbVYIFcoJX/fWOA4DnK5EhZL6DWlvW0Dx/r3tY3a24+E/Lnu5HJlQgfRQAiBdF3dYTy/ZAmqqqrc/mLP81y95m1pWkZIsvEyKp2qdhBCYsESYiDtuPUfKo0mxznAsPvfNCG/ry8WizlpAmkAUCj9B9K+tpO3bcCYCMaYz22UnZ0VfIO9UCTBNgk6kF7y3LPgeQEL7rgDOZoceJnmkBDixuvPDsXRhJAoE0UxqAFrcoUSCoXSpSazzWaF2WyCJcjZ9rpWhGCMIS0tHTzPQaPJiVhlCKvVAsaSaJZDQeZx0GGg28mxntva2qDX69DZab9b2tjYiNraGufftVotsrKyj+ZIdyA7OwtDhx4vyXeQyeTgkyCtMegc6dmXz8TDjzyKHj16RqpNhCQls8kIo4eyRYk+0IIQkngCHbAmCDKkpWf4HFwmiiIM+s6oTBQSjmQ711osZme6YDjbafv2rS69z0VFxRg2TJpg2Zd0dWZUZ2uMlKAz73v27ImOIx2RaAshSc3bLTZK7SCERFsg+dGCIEO6OtMlOKuqM+GbPzpRVXesd5PneaSrM+N+BsFkm8FRJpOD54Wwt5P75CuRnHjFQRBkSRFEAyH0SO/47TesWvUWLps5E+XlFd2qDaSnp0vaQEKShdGoh9lk7Pa8SpWeNDVOCSHxjzERHUe0fpfrWppz8349nljXjD+6BGaDSpS4eXI+Rvax/+7HqgpI4Lij1TuSI70DACxmEwSZPKztdKzUnRYajSYiJQjdJdPdgaAD6RnTL3G81O0vNNiQEF8MBp3HXEKlMg1KVVoMWkQISUWBpHXIFUqkpakB2IOz+UsPwephAkQZD7wwt4czSDPodWFVk4i0ZArgALjkfUdrO3WtMe3Iaw8m8OYFAWp1VtJc0AR9H+beexdFoh2EJD9vqR002pAQEkWBpDh0rabwxLpmj8EZAFhF+99X32AfnBZINYlYslrMSRVIdw1Go7Wdwp3eXaVMS5ogGgghkB48ZEgk2kFI0qMcaUJIrDHGYLX6nxbckUdbVWdySRPw5I86E3bXmzCgWBn3edIWqwWqJKre4RDN7RTO9O6CIIOQJLnRDgGtyerqAygrKwfP86iuPuBz2WjMy05IIvIaMFMgTQiJkkB6o7sGmXVa/0G3Y7kBxUrn693Pd+GmA0iGMdisFsiSqFcakG47BcK9xnQwgxOVquTqjQYCDKRvu/UWvPTyq8jOzsZtt94Ce360p5VPOdKEeOMthYN6pAkh0RLIJCxdz0klmsB6D7su5+mcFm46gJQsFnPSBdJSbadAdK0x7RicGAiZTJ40lTq6CiiQfva5JcjKynL+PyEkBJTaQQiJoUDTOgD7JB6CIENliRKDSpQ+0wYGlSqdvZzeakmHkw4gtWSbnAWAZNspEBzHoby8Iujtp1QlZ1W3gOpIFxQUOne4pqYm5ObmoaCg0OVfbm4empqaItpYQhIZ5UgTQmIpmDrK5i4Vhm6enA+Zl2hBxgM3T8o/9jqT50AuFrWKvWFH0zuSjRTbKVLkCmW3csnJIugJWRbdd59zKsmu9Ho9Ft13nxRtIiTpMMYokCaExFQgaR3OZc0miKK9BMTIPul4YW4PDCpRuiwzqETpUlJNFEWvlSDKyspRWTkQRUXFqKwcGLGpwANlsSRfIC3FdooMDkpl8pZ4DWHYJutWQRoAOjo6oFIpPfyFEOILBdKEkEgLJq3DwaDvRLo6ExzHYWSfdKy+oRy7602o01pQopE70wQc7++YrtqTUNMBIsVqNYOx9KRK7wAQ9naKBKVK5XPq8kQXcCD9+GOPHv0/Ds899yzk8mMJ46Ioorq6GgMGVErdPkKSgu9gmSVdvh4hJL6EMj22zWaFXteBtPQMZyA0oFjpEpgB9hjAoO/0mHcbN9U6PLTLZrVCJk++wW9A8NspUjieh0KR3DP3BhxIH5v6myEtLQ0KxbERrzKZDP0HDMCECWdJ3T5CkgJjXqrkH1sC3WcLJYQQaQST1tGVzWZFZ4cWcrkSCqVr/WGbzQqzyeQzTSCeqnW4s1fvSM5AuqtAtlOkqFTJ1+vvLuBA+vobbgQAFBQW4vzz/wKVKrmvMAiRkr/0DXuPdJQaQwhJKaGkdbizWI4FYsHUH46nah3ukjW9A7Bvc72uI6q9z+4EQZaU5e7cBZ20Mm3adAqiCQlSIIE0IYREQrhBtLtgzlfxVK3DHWMspoFmJHEcB6UqtgP8VGnJeZHiLqQ5In/68Uf88OP3aGluhtXquhM+8ujjkjSMkGRCgTQhJFZCTeuQQqiTd0SLxWJO2l5TxwQoUl9IBUIuV8T9dPFSCbpH+pNP/oslS56FJluDP//8E3379UdGRiYaGhowfPgJkWgjIQnPbyAt+suhJoSQ4DHGQhpoKBVHtY5hw45HeXlF3PVQWi3mpO7IiNUkKMk6+YonQV8ufP7Zp/jrtfMxZszp+Oabr3HBBRegqKgYq1e95bG+NCHE/2DDZD6RE0JiJxa9kYnEkd6RrL3SgiBArlDCYo7eQEOFMrnL3bkL+ps2NzejstJe5k6hUMBgMAIAzhg7Ft9/v1Ha1hGSJCi1gxASC7HsjU4U1iScnKWrqE6GwnFQKlNrHF3QgbRGo3H2POfnF2DPnt0AgMbGRlAsQIhn/gNpSu0ghEhLimodqcBiTe70Dp7noYhScKtUqsBxqdMbDYSQ2nHccUPxyy+/oHfvPhg3fjzeWPo6fvrxR+zfvw8jR46SvIGtLS1YvmI5tv5vC0wmM4qLi3H9DTegb99+kn8WIZHiLwc6mU/ihJDYsNmsdG4JABNFiKItqQfHKZUqmE0m2OcsiAyO45J+8hVPOKvVEtRaFUURjDEIggAA+P77jaiqqkJJSQnOPutsSYubd3Z24vbbbsGQIcfhnHMmIisrC3X1dSgqKkZxcbFkn0NIpHV2tEMUbV7/LpPJka7OjGKLCCHJzmjQw2w2xroZCUGpTIt5ubhIMxkNMJkMEXt/pSo95dI6gCADaZvNhrXvvYvxZ05AXl5eJNsFAFixfBmqqqpw/+IHIv5ZhERSx5E2nz1DvCAgIyM7ii0ihCS7zg4txDitCBRvU4enwjmYMREdR9oRiV5pjuOQkamJu6os0RDUfQxBEPDBBx/gjLHjItQcV7/88guOHz4cTz7xOH7/fSdyc/NwzsSJOOuss72+Jl5PGiR1McYCKH/HaN8lhEhGFG1xfU6Jt6nDRZsNVqs16atNyBWKiFTwkMuVAf3WJZJA94WgE4KGDh2K33/ficLCwqAbFazGxgZ88flnmHLe+bho6lTs27sXr7/2GmQyGcaNG+/xNc1N9RFvFyHB8jcbqCjaaN8lhEhGEATIJUy1lFo8Th2ubWuO64sPqSiVSkl7jhljOHJEK9n7xYvCotKAlgs6kB5+wglYuWI5ampq0KdPH6jc8mFOOvnkYN/SK1Fk6Nu3L2bOnAUA6N27D2pqa/HF5597DaTzCyh3msQXUbRBr+vwuQzHccjLL0rJ22KEEOkZ9J1xPf21RpPj7Im2P9bErjFHpaeroUpTx7oZEWc06GG1SlcWUaFQITMrx/+CSSroQPrVV14GAPz34488/JXD6jVvh9smp5wcDXr27OnyXM8ePbDpp5+8vibZb8uQxBPorIUcx9H+SwgJm2OSkXgWj1OH22xWcByX9B0aSpUK1k7pAmllik3A4i7oQHr1mnci0Q6PKisH4vDhwy7PHa6rQ0FBQdTaQEi4xABrRNtrSafuyYgQIg2bNb6DaODY1OGxTufoijEG0WaDIEveMngAIAgyCIJMkostQSYHf7SKW6qK61/tKeedjz17duO9995FfV0dNm7YgPVffoGJkybFummEBCzQwRdMTJ5BGoSQ2KFJWEKXKutOrlBK8j4Kid4nkQVdRzrafv31F6xcsQL19XUoLCzElPPO91m1g5B4YzIZYDL6r92ZlqaW7ORGCEldnZ3tEG3e69YT7wSZDGp1VqybEXGMMXQcafO/oC8ch8wULXnXVdwH0oQkukAnRUjVYvaEEOnYawVrY92MhJaZlZMSwaFe3wmrJfRcablcibT05B+c6U9cp3YQkgxYUDnShBASOmsC5EfHO1uK9ObL5YrwXq8I7/XJggJpQiKMcqQJIdES79U6EkGqrEOZLJw64xwEIbkHZQYqoLWg1+sDfsP09PSQG0NIMgquagchhIQuESp2xLtUCaQ5joNMJg9pgKVMLk+J9JdABBRIz5t7BQB/K4xB6jrShCSDQHuak2lqVUJI9CVC/ehEkEoXI0KogXRYvdnJJaBA+t57F0W6HYQkLcqRJoREQypMbx0NjIkQRTElJhmRyWQwhfg6YhfQmhg8ZEik20FIUgqml1mkHGlCSBhE6o2WjCjaUiKQ5nkB9oyDwH9/OI4/+joChDCzoYPJZEJzc1O3EcIVFb3CbRMhSSO4XmYGxhjlnRFCQmITU6PaRDTYbNaUSF/gOA6CTAZbEOkdyT7zY7CCXhtH2tuxZMlz+N///ufx75QjTcgxweY9UyBNCAkVTcIiHdGWOmkygiAEF0in+JTg7oK+b7F06evQ6XT410MPQaFQ4K677sYNN96IkpJi3H77gki0kZCExYLMWaQ8aUJIqChHWjpiCvXuC0GmaQS7fLILukd6x47fcNttC9C3bz/wPIf8ggIMO/54pKelY+3a9zDixBMj0U5CEpIYQo80IYQEizGWUsFfpKXSuuSD7GEOdvlkF3SPtMlkQlZ2NgBArVbjyJEjAIDy8nL8+ed+aVtHSIILtoc52B5sQggB6CJcaoyxlFmnwQ4c5LjkH4QZjKDXRmlpKQ4fPgTAPrDwyy8+R2tLCz7/4nPk5ORI3kBCElkoOdKEEBKsVOpBjZZUSZXhOC7gsTk8L9A4HjdBp3ZMPncKtG1tAIBp06bjwQcfwIYNGyCTyXDDDTdK3kBCElmw034HOgsiIYR0lcgX4Ywx1NbWQKttg0aTg7Ky8rgI1ux3FFMjjYHnhYAm80mFkoDB4qxWS1hHn8lkwqFDh5Cfn4+srCyp2kVIUtDrOoKaNUquUCItTR3BFhFCkpHZZITRqI91M0JSU1ONqqpdzseVlQNRXl7hcVleECCXK2GxmCJepSQtTQ25QhnRz4gXgf5WyeVKpKXTb1RXYV9ayGUy8BxHVymEeECpHYSQaEjkc4dW2+b2WOtxOUEmh1qdBaVSBbU6C0KE6zwn8joNFhdgDMfxsb9TEG+CL3/3+mv4av2XAOw1K++9dyFuv/1WXDf/WuzcuUPyBhKSyGiwISEkGhI56NNoctwea7otI8jkSE/PcKZ8cByH9PSMiAbTibxOg8UHOICQBhp2F/Qa+emnH1HRqxcA4Jdff0FjYxOe+vczmHLeeVj11kqp20dIQqMeaUJINLAgpniON2Vl5aisHIiiomJUVg5EWVm5y9/dg2ixvQlA5IPpRF6nQQswJz0ectfjTdCBdEdHh/Pq8X9btuCUU05BaWkpzhx/JmpqaiRvICGJKpTySRRIE0JCEsfnDsYYamqqsX37VtTUVHc7z3Ech/LyCgwbdjzKyytcgjX3INq8ZT3a75sG85b1ztdGLJiO43UqtUADZAqkuws6kM7OzsbBg7UQbTZs3boVw4YdDwAwmU2UJ02Ii+BPwjSzISEk2dTW1qCqahcaGhpQVbULtbWBdbp5CqJ1S+8FLCbolt4bnWA6RQQcSIMCaXdBl78bN/5MPPXkE8jJyQHHAUOHDQMA7NmzB6WlPSRvICGJKtTeZcYYXfUTQpKGp8GE3qpyOHgNoh31skWb/TEAxYgJzmBar++ELYhKSSRI9NPUTdCB9PTpM1BeVo7mlmaccsqpkMvtV4A8z+PCiy6SvIGEJCpvNaT91UylQJoQksi6nuOysrLR2alz+XtDQz2+/LIBACAIAnJycjFs2PHOu9p+g2gHCqYlQ6kdoQu7jjQhxDOrxQK9vqPb8/5qpqozsiAIQV/jEkJSmMGgg8Vs8rnMoUOHcOjQ4Yi3paGhATW11UG9JkeTg379+iMvLx/Djh8OQbBPhOI1iO6KF6CeuwiKERMAADabDdu3bUVLS7PPz+zRoxQ9eni/k65QqqBSpQf1PRKV1WqBXtf998od/T51F/TaeOftNT7/fsm06SE3hpBk4i3f2d9tThpwSAgJViC5q6tXv42n/v1MFFoTmhEjRmDduk+CC6KBbj3TgiBg0OAhmDz5XGzZssXry/7vH3/HP//5D69/T6184EC/ayqtk8AEHUhv3rzZ5bHNZkVjYyN4XkBxcREF0oQc5S0g1mhy0NDQ0OWxJqDXEUKIV0kQ38yYMQMqlQqAvcSdbtli/0G0g2iDbtliyPoOA59dAJVKhenTp/sMpP1KoTSG1Pmm0gs6kH70sce7PafX67HkuWcxcuRISRpFSDLwFhA7aqRqtVpoNJpuNVMpkCaEBCvQCTXi2erVqzF37hyoVCrw2QVQz14YWI80YE/vmL0QfHYBAMBoNGLNGt930P2+ZQoF0oFG0im0RgImWY50TXU1HnnkITy35AUp3o6QhGc06GE2G4N+nVKVDqVSFYEWEUKSlcVihkHf6XOZaOVIhyOecqTT1ZmQpUhJPZvNCl3nEb/LZWRmg+eFKLQocUiWMa7X66HX66V6O0ISXujl76iWNCEkOIFUU+jRo4fPwDFeGI16Z9UOR3DsNZh2C6IZYzAa9ejbtzf69u0dVjtoOmwSiKAD6U8++a/rE4yhra0N3333HYYPP0GqdhGS8EKeXpZSOwghQUqmXkKb1QK9vtN/MO0hiJay9B1NMucB/Tx1E3Qg/d+PP3Z5zPMcsrKyMHbcOFx00VTJGkZIogu1Z5lypAkhwUq2+r5+g+kIB9EcxyXdOiWRQXWkCYmQzs52iLZjvSf+JmJxkMnkSFdnRrOphJAk0NnRDjHQKhcJwuPkLMsWQz17YcSCaAAQBBnUGVmSvV+8s1os2LVrh9/fJ3VGtjN/ndhRVW1CIsWtZ7m2tsY5EYuj/J2naXKpR5oQEgpeEJIukPbUM+0ocQdEJogG7Osylezfvzeg3yfSXUiB9L59e/HjDz+gubkZVqvV5W+33HqbJA0jJNG5B8T+JmJxvo6S0AghIRB4AVb/iyUc92A60kE0YF+XqcS9wom33yfSXdCZ9N9/vxF333UXDh46iM2bN8Fqs6K2thY7dvyG9PTUmEqTkEC4B9IaTY7LY71eh5qa6u490NQjTQgJQTz3ojLGUFNTje3bt3o+7/nhCKYdr4tkEA0AfIpNg52Xl+/y2H2iMOJd0HvK2vfexZy5czFp0mRcMXsW5s27EoWFRXjpxReQk5Pj/w0ISQGefiQcE68cPnwIHR0d6OjocN5KoynCCSHhEuI4+As0tc0Xm9UCne4I5HIlLBaTyxgUqaVaHnCfPv1gNhu9ThRGvAu6R7qhoQEjRpwIAJDJZDAZTeA4DlPOOx9ffvmF5A0kJDF1D4Y5jkN5eUW3Ozdardb+iqM9Ntu2/Q979lRRQE0ICQrP8+DipGSbew90W1ury98d571giTYbTEZ9RINoXhBSrmKH4/dp2LDjUV5ekXLfPxxBX76q1WoYDQYAQG5uLmpqa1BeUQG9TgeTySx5AwlJRL6CYI0mx9kjY3+sAeC5x6Z//8rINZIQknRkggwWMfa/xe7ns8LCQpe/x3PqgCyOe/ZJ/Al6bxk0aDC2b9+G8ooKjD7lVCx9/TXs2LEDv23fhqFDh0aijYQkHF+BtOOWmfstNPfBiK0tzQAF0oSQIAgyOSyW2AfS7uczjuNQWTkwIVIHhBSZFpxII+hA+qqrrobZYk/unzr1YsgEAVVVVRg1ajSmXnyJ5A0kJBH5CqQdt9Dc8wPde6pz3QZ/EEKIPzJZfPSmdr/zluPxvOcQaJ39aIiXdRhJjDHs3bsbLS1NyMsrQO/efQJ6HSV8dBf03pKReWyiCJ7ncSHNZkhIdyGkN3ftqS4sLEa/fgMkbhQhJNnxvACe5yGKoc2sKhVvd968kWIwohTs+dHxkWceSXv37sb2bVsAAIcO1oKJIoqLi/y/kCLpboIOpLds+RU8z2P48BNcnt+2bStEUcQJJ4yQrHGEJKpQakF37alOV2fRYA9CSEhkMjnMZlNM2+Dtzps3gdbZjzRZiqR1tLQ0uTxubW0OLJAm3QR92bVyxXKPV7pMZFi5YrkkjSIk0YVdcYMqdhBCQpSIOb7udfZjNRhRJlPE5HOjLS+vwOVxbm6gqYTUweMu6B7purp69OxZ1u350h49UF9fL0mjCEl4YQbCNLshISRUidirGmwqSCRwHJcy9aMdqYOtLc3IzctHn959odMdiXGrElPQgXR6ejoaPZSyqa+vg1KplKxhhCSycANhqiFNCAkVx3GQyRWwxkH1jkAFmwoSCTK5ImVS6jiOs5dXPVoZShQDq8udIqsnKEGndpx88slYuvQ1l97n+ro6LHvzDZx00smSNo6QREWpHYSQWJInYK90rKX2Ogs0QqZI2h1ntVqC+sXW63R48MEHsH//PuTm5gEAWltbMHDgINxy621Qq9URaSghicRo0IU12EepSoNSmSZhiwghqYQxho4jbf4XJADsPbQZmZqU6ZF2x5iIjiNav8tlZuWk7DryJuhAGrAfoNu3b0P1gQNQKBQor6jA4MFDItE+QhKSQd8Z1qQICqUKKlW6/wUJIcQLva4DVqsl1s1ICHKFEmlpqdsRGOiFV1Z2bhRak1hCCqRj5f2172HlyhU499wpmDvvylg3hxCv9PrOsPITFQolVCl8UieEhM9iMcOg74x1MxJCujozIQdpSiWwQJpDVnaOn2VST8JUHd+7dy+++OILVFTEbiACIYEKN0eaUqQJIeGyB4Z0G94fe7WO5J/N0JdA0jUoo8OzhAikjQYD/vPMv3Ht/PlQqzNi3RxC/KPyd4SQGOM4DnJ5atRFDodcrqS8XyCASJnWkScJEUi/8uorOGHEiRg27PhYN4WQgITfI02BNCEkfHIFBdL+0Dqy8xtG08WGR3F/L+P77zfiz/378dDDjwS0vKdZFwmJtvADaZH2ZUJI2DiOB8fxYIzOJ57wvACO4+l8C8AeSvv+7Uql9cTzgfU1hxRIi6KI+vp6HGlvh+h2cEpZvaO5uRlLX38Ndy+8B4oArxibm2h2RRJ7SmV4twqtFgvty4QQSQiCALk8dQfS+WIyGaHX62LdjLigUCh8Bo8WizmlfpcKi0oDWi7oqh27d+/GM08/haamZnS/cuGwes3bwbydT5s3b8Ljjz3qsmFFUQTHceA4DitXrgLvNp1nKl0tkfjEGIOusz2s9+A4DuqMbIlaREjqYoxh3749zqmQ+/btn3K3qEVRhJ6mf/ZInZEFjkuILNeI0+s7INq8z3AoCDKkpafOOLVAe6SDDqRvveVmlJSWYPr0S5GTk9MtpyZdwglZDAYDmpqaXJ57fsmzKC3tgQsuvAjl5eWSfRYhUpFmIgQqM0RIIKxWi0sHCs8LkMmO3Wzds6cK27dtcT4edvwI+9TIKUan64CNakq7kMkVSE+hwNAff3XHaX15FnRqR319HW6++RYUl5REoj0u0tLSugXLSqUKmZmZFESTuCXNQEEGxljK9ZwREiyzyejy4y9XKF0C6ZYW186Y1pZmIAUDaYVCCQMF0i4UCmWsmxBf/PzecFS1w6Og72f069cf9fWpkyNDSLCo4gYh0eNeKtL9+MvLK3B5nJuXH/E2xSOZTE4X5l1wPJ/ytaPd+a/aEZVmJJyg96LJk8/Fm28uhVbbhvLyCggy1xzliopeUrXNo/sW3R/R9yckXFKNjmdMBMcJ/hckJJW5X7e6BdL9+g0AAGeOtONxqrHXlFbCbDbGuilxQaGg2tHu/K4PWl8eBZ0jPWP6JZ7eBvazmbSDDQlJRFJNy6vOyKIeE0L86OxohygeGyAlCDKoM7Ji2KL4ZbPZwh4IHQzGGGpra6DVtkGjyUFZWXm3YC2QZSIhI1MT8GCyVGE06mE2eb/QUqrSoFSmRbFFiSHoX+lnn1sSiXYQkjSkSu2gFBFC/HM/TmhWUO8EQYAgyGCzWaPyebW1Naiq2gUAaGhoAACUl1cEvYzUZDI5BdEe+LuAoR58z4IOpAsKCiPRDkKSBgXShEST23FCx41PCoUSBkN0Ammtts3tsbZbkBzIMlJTKFQRff9E5W8wIQ029Czk+8YHa2vR3NwMq9X1gDzp5JPDbhQhiUyyHGmqiU6IT4yx7j3SFEj7JJMrwBn1UVlPGk2Os5fZ/lgT0jJS4ngegoxS5jyhHOnQBL03NTTU4/HHHkVNTQ1cp5O0r2DKkSapjnqkCYkdOm58i+agw7Iye5larVYLjUbjfBzsMlJSyGmQoVeU2hGSoAPp119/DQWFhVh4z3248Ybr8K+HHkFnRwfefPMNzL7iiki0kZCEIlVPslQ924QkK2/HCNVg902uiE4gzXEcyssrfKZqBLKMlORUO9orypEOTdDZ9nt278aMGZciK8s+rSbPcRg4aBBmzpyF1197LRJtJCShSNUjJlLPGiE+eTvW6CLUN8egw1RDgwx98xtIU460R0HvUaIoIk1lL3+SlZWJ1rZWAEB+QQEOHz4kbesISUCUI01IdHg7RphIF6H+pOKsfqn4nYPhN1CmHmmPgr4kLSsrx4HqAygsKkK/fv3x4QcfQCaT48svv0BRUVEk2khIQhEptYOQqPB210ZkImgqI99kcgVg1KdMlROO4yHI5LFuRnyj1I6QBN0jPfXii52302bMuBSNjQ2495678b8tWzBv3lWSN5CQRCLlQCepAnJCklXXiVi6ors5/nEcB4VcEetmRA3NZOgf5UiHJuiZDT3p7OiAOiODVjJJeVLPHJaZlUPHFSFeGPSdsFjM3Z5XKFRQpaXHoEWJxWazQtd5JNbNiIqMzGzwPN2n8IUxho4jbR7/xnEcMrNyotyixBBy1n19XR22bv0fzCYTMjIzpWwTIQlL6nQM6pUmxDtvx4e3nmriShBk4IXkDy4FmYyC6ABwHOej44Y6dLwJOke6o6MDTz35BHbu3AGAwzP/+Q+Kiorx/PNLkKFW44o5c6VvJSEJQurAl4kikAI/dIQEizEGm5eAmQLpwCnkShht+lg3I6IUchpkGLiu84N0eZbujHoVdI/0G0tfhyAIWPL8i1Aqj+VXnXrqqdi69X+SNo6QRCP1D7jIKCAgxBPGRK8D5URRpIlZAiRP+jxpzj6wkgTEW8BMgbR3QQfS27Ztw6zLZyMvL8/l+ZKSEjQ1NUvWMEISkdSDnCi1gxDPbDbfF5k2mzVKLUlsHM9DlsTVLOQKBQWBQaBAOnhBB9Imk9GlJ9qhs7MTcnnqFXgnpCupA18KpAnxzF+g7C/QJsck82x/yd/jLjFvATMF0l4FHUgPGjQI3377rfMxBw6iKOKDD97HkCHHSdo4QhKN5KkdFAwQ4pHN6ieQtlqi1JLEJ5PJk7LHkeP4lJzBMRxee6RpsKFXQe9hsy6/Aovvvw/79+2D1WrF8uXLUFtbi87OTix+4MEINJGQxMCY9HmZomgDYywpf+QICRVjLIAeaSsdOwHiOHsescVsinVTJEVpHcHzFjDTevQupDrSep0On366DgeqD8BkNKJ37z6YOGkycnKoxiBJXTarFTqd9DVZMzI14PmQK1USknQsFjMM+k6/y6nVWRBk1CMZCKvVAr2uI9bNkJQ6IxsCVT0KitGgh9ls7Pa8UpkGpSotBi2KfyGdYdLVaky9+BKp20JIQvNWiitcos1GgTQhXVgDTNuwWi0USAdIEGTgOF7yWvixwvMCBdEh8NrzTD3SXoV0hjGbzaiprkb7kXYw0bVD+6STT5akYYQkmkjlM9tEG2RI3lH1hASDMQarJbBA2mIxUy9agDiOg1yu8NgbmYhokGGI/FTtYIxh797daGlpQl5eAfr1G5DyaR9BB9Jb//c/PPvsM+jo8HQLiMPqNW9L0CxCEk/keqSpjBchDqLNFnCvqSjaIIo2mtUuQLIkCqSpdnRo/JW/27t3N7Zv2wIAOHSwFgDQv39ldBoXp4IOpF977RWMPuVUXHLJNGg0mgg0iZDEFKmAl8p4EXKMxWIObnkz9UoHShCEpEjvoLSO0MlkMqjS1N2ed6zPlpYml+dbW5oBCqSD097ejvPOO5+CaEK6kGImNcYYamtroNW2QaPJQVlZOTiOo8odhBzFGIPF4r2yhKdjiNI7Amev3iGXtHqHt/NaJMnklAoXKp4XoFB4vwjJyytw9kQDQG5efjSaFdeCDqRHjT4Fv+/cgeLi4ki0h5CEJMUsarW1Naiq2gUAaGhoAACUl1cAsN+ipnqoJNXZrFafF6zejiGbzUrHT4DkEpfB83VeixTKj46cfv0GALD3ROfm5Tsfp7KgzyxXXXU1nnzicfzxxx8or6jodvvk3HOnSNY4QhKFFIG0Vtvm9ljr/MGxWSkQIMTsozca8H4MWcwmCGl0/ATCfp7hAEhTE9/XeS0SOI6nnPgI4jjOnhOd4ukcXQV9Zvl+4wZs374Ncrkcv/++E+hSvJvjKJAmqUmKQFqjyXH22NgfayR9f0ISmSiKsPrJj/Z2DJktZihV6ZQeFQBHeoe/dR0oX+e1SJDJk3OWRhK/gg6k33rrLUybPgMXXngR1bYlBEdnWfMzXXEgysrKAdh7bDQajfMxQIE0IYEMMvR6DDEGi8UMhUIZySYmDZlMukDa13ktEmQyyo8m0RV0IG21WnHqqadREE3IUVLVj+Y4DuXlFR5ve4qiCFEU6bgjKYkxBksAZdl8HUMWs4kC6QBJGYz62iaRQIE0ibagf5XHjhuHH374PhJtISQhWW3eJ4dgjKGmphrbt29FTU11WJU9bAHO5kZIsrHZrBDF8Eqy2WxWurMTIJ6PTJ6xlOdDTwSZjNI6SNQF3SMtiiI+/OB9bNu6FRUeBhvOmTtPssYRkgh8pXVIOWLdarNCDupRI6nHLFEVCbPZhDQadBgQmVwOs0naGvaRruBBvdEkFoI+o9TWVKN37972/6+tcfsrXQmS1MIYg9VHT7GUI9atFguYiupJk9QSyCDDQFnMJqho0GFAZIIcZkg7y2GkK3hQIE1iIehA+t777o9EOwhJSP5uFUs5Yp0xe540zdhFUomnmsbhTPJhsZigUKikbmbSEWTS99xHtoIHR2XvSEzQPS5CwuCrNxqQfsS61WqhQJqkDMaYx9rR4aQImM0USAeC4zgIgkzSvPJIVvCQUX40iREKpAkJg79bzlKPWLdazFAqKQggqcFms4J5GGQYToqAaLPRTIcBEmTSBtKRrOARiR50QgJBtbQICZFos4VdSSBYNpsVjEX3MwmJFW9TVWs0OW6PNZK8L3EViYuNSFXuEATKjyaxQZdwhITIYpVmAFSwrBYL5FQPlyQ5dnQSFU/CTRGw0EyHAZFFoJc3UpU7KOWNxAoF0oSEKJCZ1iL1uRRIk2Tn6/gKN0XAUW1HLleE2ryUwHE8eJ6X9M5bJCp3CALlR5PYodQOQkIgijbJZjQMltVqofQOkvSCKXkXSrqAVCX1kp3U6R3hpuV4Qr3RJJaoR5qQEFjMsf0RtlgsNN0xSVr+6rO7CyVdwH5BSnXZ/REEmaR33yJRuYMGjpJYor2PkCDZczelH6wUTG1ci9lEgXSEMMawd+9utLQ0IS+vAP36DaBgK8qCCaIBoK2t1e1xm99AmjEGm81Kk3j4wUscpEaicofUbSQkGLT3ERKkSFXrCKZXzWazQhRtNAFBBOzduxvbt20BABw6WAsA6N+/MpZNSjlWS3CBdHeBVYKwWi0USPsR/2kTHHieslRJ7NDeR0iQPE0QIQVPg3B8iXV6SbJqaWlyedzUWA+DQQczlUyLGqst2EDa9Y5BoHcQbEH2fKcijovvGQMFQaA7RiSmKJAmJAiMsYjVoA12EI7ZbJKsBis5Ji+vwOVxVlYWLGYT1R6OElEUPU7C4oljkKHBoHd53v1Y8sZms9ExFAA+jnul47ltJDVQagchQYhEbrRDsINwGBOphFcE9O3bHyajvtt2EMXYVGlJNcHMpNc1HQoAMjMzUVraI6gBbJQn7Z/AC5BufkNpCXHcW05SAwXShASIMQazKXKBdCiDcCxmEwXSEmOMedwOjDEwJoLj6EZeJAUTSLunQ6Wnq4MexGaz2SiQ9iOee33juW0kNdAvAiEBcgzwiydWqyVm9ayTla8a3dGeEj4VBbM/S1GTmI4f/+K51zf+B0OSZBfXPdJr176HzZt+wqFDh6BQKDCgshKXz5qN0h49Yt00koLMJmOsm+CR2WyEKk0d62YkDV/BMhNFgH63IyqYi1UpahLH28VxPOLitCoGx3F0h4jEXFwH0r/v3ImJEyehb79+sNlEvLVyBR544H48+dTTUKlUsW4eSSGiaAu6tm20mM0mKFVp9IMiEV890jQwLbIYY0H1+ktRk5gCaf8clTvibV3FczURkjo4q9WSML8MR9rbcfXVV+K+Rfdj8OAhsW4OSSEGgy7oqg3BTLASLqUyDUpVWkTeO9UYjXqYTUaP2y8tTQ2Fki7iI0UURXR2aL3+3dcxFc7xlpmVQyXU/NDrOuKuM0EuVyItne7GkdiK6x5pd3q9vcRRRkam12Uoh5FIjYliSKXPQpm2OFRmsxEyuYKCAQk4ep09bb9+/SvpHBNB/gYa+jqmwjnebDYr9W76EY/pHRzP0fFIIibQiX4SJpAWRRFLl76OysqBKC/3ngfX3FQfxVaRVCCTySCTBX+oeJpgJVKBNGMM2rZm2GjgVNgc29vT9tN1duCITRubhqUAnuehUHivQuPrmArneGtrbaa0HT8EQYBcHl/VTTo7jkAUtbFuBklShUWlAS2XMIH0q6+8jNraGty/+EGfy+UXFEepRSQVMCZC13kkpNdqNDnOnjH7Y41ErfJMoVAgXZ1FvdJhMhkNsFhMHrefOiMTCgWldkSKxWyCyWTw+ndfx1Q4x1t2dg5kVEbSJ6vVAqNBF9RrIp3elq3JhSAkTBhDklRC7IGvvvIytmz5FYsWLUZeXp7PZQPtiickECZj6HWjpagoEAzGGGw2KxQKZUQ/J9lxvP2H3tP24zmezjEx5OuYCut44+i3w59Q7spFOr1NJpPRIGsSc3E92JAxhtdefQWbN2/GfYsWoaQksG52QqTAmIiOI+0A4vYQ6YbjeGRkZlOvdBhMRoPXXlFVmpouVCLIaNDDbI5+mUkarOsfYwwdR9r8L9jF9u1bXe4SFBUVY9iw4yVqEYes7MCmgickkuK6R/rVV17Gxo0bcNttC5CmSoO2zX4Qp6enQ6GkHzMSWSaTEYkURAP24N9iNlFliSAwxrB37260tDQhL6/AZ08mXZ5Elq/Sg8eWCTxdINBlKT/aP0fN5kC2kUMk09voDgKJF3EdSH/++WcAgPvuu8fl+euvvwHjxp8ZiyaRFCGKYtxOwMILAuRyJSwWk8dZ2UwmI+QKJfVKB2jv3t3Yvm0LAODQwVrYbFaUlpR4XpjWaUQFEtAGky4Q6LLBBIepjOd52GyBr6tIprdRIE3iRVwH0mvefjfWTSApKl6DaEEmR3p6BjiOg0KhhF7fCZtbbVfGRJjNRiiVdKs6EC0tTS6P21pbvAbSdHESWWIAAW0w1TkCXZZ6pAPD8TwQRGEgKSbM8dkWQuIA7YmEuBFF0W+eZiwCqq5BtKMN6ekZEGTdS1I5JhQh/uXlFbg8zs31NaCZAulIYqL/fVajyXF7rAl72UACeALwcTSwL57aQlJbXPdIExILJqPngWZyhRIKhdKl3JLNZoXZbAppwpZguAfRYnsT+OwCZzDt3jPNGIPJZIBKlR7RdiWDfv0GAABaW5qRm5ePiopeXst8URgdOYyxgFIsgkkXCHRZRpN6BCSeeoHjqS0ktcV11Q5Cos1ms3arGy0IMqSlZ/jMyRNFEQZ9p9+Z2ULhHkSbt6yHbtliqGcvhGLEBAD2IMRTmkdGpoZyCYNksZhh0Hd6/Fu6OhMyD3cASPj8TQ8eaTRNuH++jo1oo2ORxAv6hSWkC/feaEGQIV2d6RKMVtWZ8M0fnaiqO9YLzfM80tWZkk8O4DGIXnovYDFBt/RemLesB+A9zcNb7zrxzldKDKXLRE6se4VpwKF/8XShQfWjSbyg1A5CjrJaLbC69eimdQliN+/X44l1zfijSwA9qESJmyfnY2SfdHAch7T0DMl61bwG0eLR0T6izf4YgGLEBI9pHhaLCQqbkmb/CoKvgC7WwV4yE8XYTm8viiJ4XohpG+JdPOUl83z8BPUktcXPUUFIDDHGYDToXZ6TK5TOnujN+/WYv/SQSxANAH/UmTB/6SFs3m9/Lc/zkMvDr3HuN4h2OBpM++qZpl7p4Nh8BHSxDvaSma/1Hg2eSkkSV5HKS2aMoaamGtu3b0VNTXWAd34okCbxgbqpSEqpra3FwYMHuz1vsVhgdpvNbtToU5z//8S6Zli9dEZaRfvfV99gH8hktdqwefPPIbcxLy8fw44f7j+IdvDQM61SpWP7tq1oaWkGAKhU6RC8TPHbs2dPlJWVhdzeZMIYg83qPc/dGoEceGIX60DWRoG0X5FK7Qh2KnGO4+MqzYSkNgqkSUpZvnw5Hn74Eb/LpaWlob6+DoA9J9q9J9rdH3Um7K43YUCxEplZWbh89lwYDMH3BI8YMQLr1n0CQbDfYvYbRDu4BdOCIGDQ4CGYPPlcbNmyxedLFyy4HXfccUfQbU1GomjzmSsr2mxHUwDoZp7UYh3IRmKgcDLieF7yFKdgaoPb20BBNIkf9GtAiAcZGRnO/6/TWnwseUzX5bq+PhgzZsyASmWf3ltsb4Ju2WL/QbSDaINu2WKI7fYJRlQqFaZPnx5SO1KVxWKWZBkSHFEUYz7Yz34RRYNJ/YlEnnQwtcEj1QZCQkV7IyEedHYeK/FUogmsxFLX5bq+PhirV6+G0WifDIbPLoB69kIg0AFQvAD17IXgs+0TjBiNRqxZsyakdqQixlhA9cAtZprsRmrx0hscL+2IZ5FIqSgrK0dl5UAUFRWjsnKg36nEKa2DxBOqI01SinuOtMlo6Fapw2HU6FOQmZkFAJjxXI3P9I5BpUqsvt5+8u84cgSbNv0YchsdOdJBpXfwAtRzFznrSttsNpccaQe5QgGFQuXyHOVI25nNxm4DTr1JS8+AXK6IcItSh8Ggi/ikRoFQKFU0iZEfBr0OFktst5VCoYIqjbYTiQ8USJOU5Wnyla7kCiXS0tQAjlXt8DTgUMYDL8ztgZF97Cd2KX5oAq7aAXQLor1NzuKQkZlNZb7cMMbQ2dEecHoBzwtQZ2RRz5gEGGPo7GyXPO+WMYba2hpotW3QaHJQVlbud3vxgoCMjGxJ25FsjEY9zCZjTNugVKVBqUyLaRsIcaDUDpKSPJW7c2cxmyAe/XEf2ScdL8ztgUElrqXtBpUoXYJoURQl6a2xWS3Q6zudKQSKEROgnruoe5pHkEE0ABipHF43ZrMpqBxdUbTBSrnSkhBFW0TqczsqQTQ0NKCqahdqa2v8t+XoYFLiXTxMhBIPbSDEgap2kJRktVoCyoc06DuRrs4Ex3EY2Scdq28ox+56E+q0FpRo5BhQfCywZoxJOn2uI5h29Ew7gmVnz3QIQTQAWC1m2KxWr+XwUg0TRZhMwV9cGI0GyOQK6pUOk9US2GDeYAVbCeJYe8xQKFV+l0tV8bC/x0MbCHGgX1KSchhjMBkDy4W12azQ6zqQlp7hLHk2oFjpEkAD9p5og75T8sFKXoPpZYuhnr0w6CDawWjUOy8QUp3JZADcBg8GkhbAmAiz2Ui3mMPAGAu4CkqwqRoaTY6zJrH9sSagz7FQIO0TH8Y5I5R0G89toB5pEj8okCYpp2vKRiBsNis6O7SQy5VQKF2n27bZrDCbTBEdfOMpmJb1HeaszhFsEA3Y2221WlJ+wJwo2mD2MMgt0AkiTEYjFHJlxGZ8S3aiaAt4tshgJ+1wVH7QarXQaDR+K0E42GxWiKKNxhF4Ec6+Huw29N4G6gAg8YPO/iSlMMZCuo0PABaLCbrOIzjS3oqOI2040t4KXeeRqIxgd8+ZDieIdjAZ9Slfxs1k9DxoylNagGcMJnNsB14lMk8XMd4Evk3sOI5DeXkFhg07HuXlFT57PnlBgFKVDt5RKScOKojEq3DuYgW7Db23gUIXEj9obyQpxWySpgZwoO/BGENNTTW2b9+KmprqsD7bPZgOJ4gGHAMjU3fAnK+Boe4TRGRnZ3vdjmaTKeUvSELBmBhUybtgJ+0IlCCTQ63OglKpglqdBUEmh8VM29SbcILYSG1DQmKJUjtIyhBDHFQWDqluZTrYrBbodEcglythsZgghjmtssmohzxFB8z5CuLc0wIYYz62I4PFYupWn5v4Zja5rn+O43wGr6GmavjiXmaS4zikp2dAr++050orlH7egQRDim3IcVxKnq9I/KI60iRlGA16mKN8G3779q0uA56KiooxbNjxUW2DP0pVOpQpOLiqs6M94Pxcf9tRkMmgVmdJ3sZkZa/brYVMroBC4WHcgdkU8Qla3INosb3JJWXKaNBBlaamoM2DjiPamE3pzvMCMjKp1jeJH9QjTVKCKIpRC6IdI9Pb2lqh17tWB/F1K1OqEe3BMpsMUCiUKRUwiKLoEkQ7UnDq6g4DAEpKSl3yat0rQIiiDT/99IPLsoyxlFqH4bBYzFBnZDsr4XQlCDKkpcmgVKZ5rITT9fiy45CTk+Ps3fS1HZ2f4WnCoy6VcDiOgypNTQNyvbDfPYjdZxMSTyiQJinBHMWUjq7pHF0VFhb6vJUpdRpIoBhjKVfGzb0nura2Brt3Vzkfd3RUOQerAa63pBkT0djY2G3ZQZka57TuxDtRFLulE1XVHavNXnl00iOe55GuzoRe1+ESTHs6vhobj13k+NqOgO9ZQ3VL7wUAZzAtk8npAsmDWK4PqpBD4g0F0iTp2XujozcK331kugPH8T5/gEKdQEIKZpMRCoUqZQIG9/KHnrZZ1/XvCMbKyyuwfftWj8syUQQokPaL444FYpv36/HEumb8UXfs+BxUosTNk/Mxsk86OI5DWnoGOju0zr97O77sFSC6d5N23Y6+gmgAHoNpGnTYHcfzQHjDM0L/7BQ5R5HEQZd2JOlFOy/afWT6sec1Qb0umiPa7b3SqVTyyzU48rTNvK1/b8syD0EccWWz2ZxVHzbv12P+0kMuQTQA/FFnwvylh7B5vz0tiud5yOVKZ/qNe7qUg0aj8bhtGBPBGPMfRDscDabNW9YD8D8IMhWFMylLuCiQJvGGeqRJUmNMhNkU3UDakQbQ1tYGwH5b2JHzHMjrpKxKEAx7r3Sq5Eq7fseysnIwxlxya72t/2CWJce4D057Yl0zrF7Gq1lF+99X32BfrwqlEvv27XZJ6cjIyEB6enq344sxhurqAzAdrQrS2NiItrY2VPTq6z+IdvDSM50ax4Z/sazjTLMaknhDVTtIUjOZDDAZo1vyLpGlpakhT4GSXxazCQaDTtL3TFdnQiaTS/qeiaC2thYHDx70uYxjKvCxY8cBsOdET3uuxu97v3NjOQYU2/fHV195EY1Nx3LTc3Ny0bdvP4+v27dvL1qPDkYsKSnF5ZfPgUxm7zfyG0R3xQtQz10ExYgJAOwpQTt37kRnZ6fPl/Xs2RNlZWX+3z9BmU1GGI2e7wxEWqqco0jioB5pkrQYY91q1RLfzGZjavxIRaBnMVV7K5cvX46HH37E73IFBQXYu3cPAKBOG9gkQnVaizOQfujhx9HU1BR0+x555BFnEC22N0G3bHFgQTRg75lethiyvsPAZxeA53ls3Pg9FixY4PNlCxbcjjvuuCPotiaKWA74o8GGJN7QHkmSltVqiVmt00Rls9lgs1r9L5jgIhP0pmYgHaiuvbglmsB67rsu568X2JvVq1fDeHQqeD67AOrZCwE+wEGhvAD17IXO+tJGoxFr1qwJqR3JJKZVO1L0gpXELwqkSdKK9IQOycrsZdrspBJEQlvXad6rqw+guvqAlynfKUvOF4PBgO3bfwMAVJYoMajE952PQaVKZ2/09u3bYTCElqK1ZcsWzJw5yxlMK0ZMgHruIv/BtFtah9FoxMyZs7Bly5aQ2pFMYpkjHcvPJsQTypEmSUkURZeSWSQ4mVk5Sd3zYzabYAwwR7qmptpjXXAAqKwc6CytlpaekZKTd3jLkRZFEUaDzuVio0ePnhg0eAiAY1U7PA04lPHAC3N7YGSfdADAH7/vxKFDvvOw/cnLy8ew44c7a337zJV2C6JtNhu2b9uKlpZmAPZgTpWW7nFCGSD5c6QZE9FxRBuTz072cxNJPBRIk6REgwzDk+wDegz6Tlgs5oCWdZ8evKuuU4UrlCqoVOmStTGR2Ww26HVHPJaNy8jUOANQf3WkAWkvigMqgecWRDPGoNd3wmZ1zeu2TxiT5TWYTmaMMXQc8VzPO5I4jkNmlufyooTECg02JEkp0CCJeGaxmJM2kGZMDGr/cJ8e3PVvGuf/W8wmKJVpKd9b5iuIBuwXMenqTHAch5F90rH6hnLsrj82s6EjnQOwB2wGfWi50R7bZrVAr+90BtOOYNkZTAcYRAP2AF+nOwK1OhN8oDnXSYLjOHAcH/UxKJTWQeIR9UiTpCOKNnR2tMe6GQkvWW+hGvQ6WILIA2eMoba2BlqtFtnZ2QCA9vZ2Z63vruso1XulbTYr9LoOvxOYCIIMaekZPntzRVGEQd/pMj24VDz2TC9bDPXshQEF0V1xnH0q81SbHr6zsx2iLbrTGwoyOdTqzKh+JiH+UCBNkk4sa5wmk2TM+TWZjDBFeN9I9rQYb6xWC/S6TgQz6FIuV0KhVEIQjt0ctdmsMJtMQV3shMI9mBbbm5zVOQINoh04jjsaTKfOTV69rgPWANePVORyBdLSM6L6mYT4Q4E0STqxOMEnI7lCibQ0daybIQnGGEwmQ9RmuVSq0lNolkh7KlC4KRixmIrbPZgGgg+iu0qlSXkMBl3UKyOl+h0fEp8o4YgkFcYYBdESsVqSYz2Kog16XUdUp4o3GfUw6Dshislfx9xsNkqSxxztIBo4ljPt+OxwgmjAfhGfKmU3YzFVN00PTuJR6tyHIilBDHTGMuIXYyJEUUzYqgT2mS2NMJliU73FarWgs6MdSpUKCoUq6XqnY71+pWKzWqDTHYFcroTFYgo779dwtOSfQqmSqIXxKRYzDNKshiQeUSBNkkokBialMpvNCp5PrDxpxhjMZiPMJmNMejndWgOT0QCzyQSlUgV5kqR7MMZgMuphTpLeV9Fmg8kmXe680ah3BtPJsL094WPwvahHmsQjCqRJUrFFeRR5shNtNiBGKZ/eJvrwRhQZrFYTLBYLEPMA2jOO4yCTKyCTKcDzwQUi8TLJB2MMRoOOSkz6YTIZwBiDUpWcJRGpR5oQOwqkSVKJdjmmZBfLVJnly5fj4YcfidnnS2HEiBGYMWMGVq9eHfbU0gsW3I477rhDopaFhoLo4JjN9rz8ZAymY5HylWzrkCQHurwjSSXRBncxxlBTU43t27eipqY6DlIRXCXa+nSXlpaGgoICpKWlRf2zJ0yYgHXrPsH8+ddi3bpPMGHChKi3QUqMMRiNegqig+RIM0o29slRohfYchxPgTSJSxRIk6TBGIv6TFvhqq2tQVXVLjQ0NKCqahdqa2ti3SQXiTh4U61WY968udiwYQPq6+uwd+8e1NfXYcOGDZg3by7U6siX9JswYQJWrlwBlco+4EylUmHlyhUJHUybzaaUqUghNZPJkDT55F1Fs1c6UQc9k+RHdaRJ0mCMoeNIW6ybEZTt27e6TD9dVFSMYcOOj2GLusvKzo3J59bW1qK2thaiaIPVYgmorKFGo8Gw44dD4WNCFLPZhO3btkKr1QbUjn379qK1rdX5ODcnF3379vO6fF5ePoYdP9w5013XiT5sNhu2b9uKlpZm3x/KcZDJ5JDL5OB4e09cLHOk7ZOtdMTks5OJOiMrqSZtiWbNfpqMhcSr5DmiScqLt7SIQGg0OS6BtEajiV1jvGCMRf2WqiiKKCzMR44mM+D0EkGQIV2d6dLWqjoT6rQWlGjkqCyxB9cKhRInnjQSel1HQFVeiosLUVW1y/m4snIgyssrPLfBz9TTgiBg+AkjgqpVzPMC5AolFIrYVE9x5EXH0rFp2tug0eSgZ88yHDxY63zsPlW7r9f6WjbSDAYd1OqspElRiObgPxpoSOIVBdIkaSRiIF1WVg4A0Gq10Gg0zsepymq1wmw2whpCHm5alwB28349nljXjD/qjt1OH1SixM2T8zGyTzo4jkNaegY6O7R+3zfQbeQxiF56LyDa7P8FoBgxwT6ddHpGwMG0KNpgMuphMuohlyugUKqi2qtpMZtinivvSIECgIaGBrS1taKxsdH5GIDXixv31/paNtJEmw1WqwVyeWKVlPQmuqkdQtQ+i5BgUGoHSRqiaENnR3usmxGQeOol8yczKyeibWOMwWa1wmQyhFwHvOt05pv36zF/6SFYPcR+Mh54YW4PjOxjn2bYoNfBYgkud5UxhurqatTUHIDVakVubi5GjDgZ6owsj0G0Ey9APXcRFCMmON8n1Fn0BJkcSqUqKtNRd3a0xzxX3j0FSqlUwmQ6tt3cU6Icx1draytaW1tcymLGOn1KkMmgVmfF7POlJMXU8IFKpenXSWJJiB7pTz9dh48+/ABarRYVFb1w5ZVXoV///rFuFok78RmIehJPvWSxZLPZYDTqYLOGN5FO15zoJ9Y1ewyiAcAq2v+++gZ7r7JCqQw6kK6trcGePVXOxzwvQ7rac0+0izB7pruyWS3QWy2QyeRQpakj1jPIGIt5EA10T4HKzs529kjb/65xWb7r8dX9vTQen48Wm9Uak3SpSKDBhoQkQNWOH77/Hm++sRSXTJuORx55DBUVFXjwwcVob0+MnkcSPYn0w6TVtrk91samIQGI1Hq1WMzQdbaHHUQDcKY6VNWZXNI5PPmjzoTd9SaX1wWj67YrKirBKaee7rzt7DWIdjgaTJu3rAcAZzAthNjTZrVa0NnZHrEBX/FSBaesrByVlQNRVFSMysqBGDr0eJfH7uk27seXg0wmS/n0KSlFM92Co1kNSZyK+x7pjz/+CBMmnIXx488EAFzz12uxZcsWfP3Velx40dQYt47EE3vAxwGIfbaSKIr47bdtaG9vR1ZWNjIzM3HgwJ9ec00bGxvxxRefOR9zHI++ffuiV6/eMb1AiFQvkNVqkeyWcNf1U6cNLKCs01owoNjei71x43cwGAzO95LJZLAeDe41mhyYzSbo9fbpo7vm4efk5OKUU093VufwG0Q7eOmZ1umOhDahEGPQ6zqgzsh2tkUq8XRxyhiDTtcJrbYNhw4dQklJCbKyslBbW4Pdu6ucy/iSm5sbV98p0XEcB47jAh6f4imljTGG7du3orW1FTKZDOXlvVBRUeGynXieakiT+BXXgbTVYsH+/ftw4UUXOZ/jeR5Dhw3D7t27Pb4m1oNiSGzxPAdRjH0g/dtv25y3npuaGtHU1OhzefeeP8ZE7N27B4IgxDTlg+P4iBxTUtbU7fojXqIJrGe363Jm87GBjYwx+xTjR7V1KXvnrry8l0uJO92yxf6DaAfRBt2yxZD1HQY+uwAcx0EuV8Jk0wf2eg/MZiOUSuknnuH5yOwDwegaLAOAyWTCnj3BlePLyMjA0KGxLy0pCLKjNe9jf56SAsfzYAFeAHpKaWtra0VTUxMAe6rXnj1V4HnO5bzHxcE+SFJPoB1JcR1IH+nogCiK0GRrXJ7XZGfj8KFDHl/T3FQfhZaReOXoIYk1qVKPWltbUVxcIsl7hYIxE3Q66QcTSb2drEfzhStLlBhUovSZ3jGoVOnsjda2tbkMRAtGTc0B9O7TD4IggM8ugHr2wsB6pAH7wMPZC531pRljMOh1YaVomEwmdByRPuUtHo6p1taWsN8jLS3deachliJ1TMVKMPuH+3ZsbW3FEQ/7rPt5j5lM0HUmzzojiaGwqDSg5eI6kA5FfkFxrJtACHLz8lF32PPFXjCKS3ogN69Qghaljpsn5/us2nHzpHzn4z8P7Av5c9raWvHjDxuc6R2Oahx+g2kP1TsYY8jKzgm5LcmuuKSns9cy9PegYynW3LdjcUkPyOTybudK2lYkkcR1IJ2VmQme56Ft17o8r21v9zrymkb2kngwevQYbPrpe7S1tSInJwfZmlzsrvrdpfeT43gIgoC0tDSYTKajKQb22708z2Pw4KHo378y5r2BiWZkn3S8MLeHzzrSDscddzyqD+x3bheO46BQKGA2W8BxQH5+AQxGI3Sd9jQC19vLHLTaNtTVHUaPHj3tr/UXTLsF0Y7PpG3sW//+lQCAmuo/YTKZoFQqUVbeCxyAAwf2Q6frPHpBAnQ9hhhjUCiU6D9gIB1LccCxDVpbmpGbl49+/QagX78B+OnHjWhuboQgyGhbkYQT93Wk77xjAfr164crr7oagP2H7PrrrsWkSZNpsCEhxKfd9cdmNnSkc0RDIHWkCSGEJL647pEGgPPOOx/PPfcf9OnbF/369ccn//0YJpMJ445W8SCEEG8GFCujGkA7dOuZpiCaEEKSUtz3SAPAp+s+wYdHJ2Tp1as35l15Jfr3HxDrZhFCiE/mLeuhW7YY6tkLKYgmhJAklBCBNCGEJCqxvclZnYMQQkhyoUCaEEIIIYSQEFCJC0IIIYQQQkJAgTQhhBBCCCEhoECaEEIIIYSQEFAgTQghhBBCSAgokCaEEEIIISQEFEgTQgghhBASgrif2TAYjDEwRtX8CCGEEEJIeDiOA8dxPpdJukC6uak+1s0ghBBCCCEJLr+g2G8gnVQTslCPNCGEEEIIkULK9UgH8oUJIYQQQgiRAg02JIQQQgghJARJ1SOd6m695Z+orq7GovsXY9CgwbFuDnGzZs1qvPP2GudjuVyOwsJCjBt/Js4//y/gebqujTe//PwzPv10Hfbt2wej0Yjc3Fwcf/zxOO/8v6C0tDTWzUtpXY8njuOgUqUhPz8fgwcPxsRJk9GzZ88Yt5A4PPjAYjQ01OOJJ/8NuVzufH7/vn24884FmDt3HiZNPjeGLSQOi+67B1qtFo899gRkXbYVADz5xOPYs2c3nnrqaajS0mLUwvhDgXSSqK2tQXV1NQBg44YNFEjHKYVCgXvuXQQAMJtN2LljB1auWA4mirjwoqkxbh3pasXyZfjgg/cxevQpuHb+fGRlZaGhoQFff/UV/v3Uk3j0scdj3cSU1/V4MhoMqKmpxpdffoH167/E/OuuxxlnjI1xCwkAXHX1Nbj5n/+Hte+9i+kzLgUAiDYbXnrpRfTu3QfnTJwU4xYSh2uuuRa33nozPvjwA1x88SXO57f+73/46acfccutt1EQ7YYC6SSxYcMGcByPwYMH46effsS8K6+CTEabN95wHI8BAwY4Hx933FDU1NRg06ZNFEjHkS1bfsUHH7yPiy++BDMuvcz5/ODBQzB+/Jn49ddfYtg64uB+PA07/nicM3ESHn7oX3jh+SWorKxEUVFxDFtIAKC4uBgXTZ2K9959B2PGnI7SHj2w7tN1OHDgTzz08CN0Ny6OlPbogQsvmor33n0XY8aMQVFRMcxmM1599RWcdPLJGDlyVKybGHdo700CjDF8v3EjjjvuOJx3/vno6OjA1q3/i3WzSIDS0tJgs9li3QzSxccffYTsbA0uvmSax7+feOJJUW4RCZRCocCVV14Fq9WK9evXx7o55KgLL7gQhYWFePnll9Dc3IzVq97CpMnnonfvPrFuGnHz/+3da1BU5wHG8WcDu0OKCworLpVYjc6wKkyMxrtEjJIGjAomRk211XhJazJp0qbOdHpR2zrRmU4+VU0CqNXEqClqtK2ZEVCEmijRmIrlUhW5KEhEEQdZLsv2g+iIuegcTc7u9v+b2S/nPWfmOezs7DPvvpw3LTVNPXs6lJmRIUnauSNLV640aMGCRSYn800U6QBQWlqqL76o07iEBD3yyBDZ7XYVFOSbHQtfw+PxyOPxqLm5WZ8WFurw4U80atQos2Ohk8fjUWlpieLj4/lVx0/FPPSQIiIi9N+yUrOjoFOw1aqFixbr5MkiLfv97xQaGqqZncs84FuCrVYtXvxTHT/+mbKy/qbduz/UzFmzFRkZaXY0n8S3RAAoKMiX1WrTyBEjFRwcrJGjRiv/YJ7czc2sZfIxLS1uzZ71XJdjY8aMVWpqmkmJcLurV6+qra1NDofD7Ci4B5GRDjU0NJgdA7eIi4tXXFy8iopO6JVXXtWDfD/5rEGDBytxwhPatvV99ev3sJKfSjY7ks9iRtrPeTweffLxIT069FF9LzRUkjRu3Di1tLToyJHDJqfD7Ww2m95YtVpvrFqtP/xxpebNf0HHj3+mt95+y+xouB3PpPdzXkm8h76kuqpKxcXFslgsOnmyyOw4uIMbEzxTpkzVA0FBJqfxXRRpP/f558fV2Niox4Y9pqamJjU1NalPnx+oR48eKigoMDsebmOxPKD+/Qeof/8BcrlcSkmZrGefnaED+3NVWVlpdjxIstvtslptunjxotlRcA/q6+vVvXt3s2Ogk9frVXr6O4qOduqFBQuVm5ujsrIys2PhG9xY2sYSt2/GX8fP3SjLa9eukdau6TLW2NioK1euKDw83IxouEu9O593W11VpT59+picBkFBQXK5YlV04t/yeDwKYibG71RVVerSpUtKTJxgdhR0OnBgv0pKirVs+QoNHDhI+QcPKiP9Ha1atZrZTvg1ZqT9WEtLiz4tPKLhw0do2fIVXV4/f/U1eTweHTr0L7Nj4g6qOmei7WF2k5PghslPT1FDQ4N27Mj6yvFjx45+x4lwt1pbW7U+M1NWq1VPTJxkdhzo+v8dvLt5k8aPT9SgQYNlsVi0cNFiVVZWaO9He82OB9wTZqT9WGHhEbndbiWnpGjw4Lgvje/+cJcK8vOVzI5RPsPr7bj5c2Z7e5vOnDmjrKwsxcTEsImODxk6dJimTkvVB9u3qbqqSmPHjpM9zK66ujrtz83VtWvXNHToMLNj/t+79fPkdjersrJSOdn7dOHCBS156WVFRUWZnBCStHnzJknSnLk/vnmsb9++eio5Rdu3bdXo0WMUERFhVjzgnlCk/VhBQb4cDsdXlmhJGj8+URs3blBtba2cTjYl8AWtra367W9+Len6EoLIyEglPP64Zsx4jnVoPmbOnLmKjY3VR3v3at26NXK7W65vET5kiKZOmWp2PKjr5ykkJEQ9e0YpLi5er/9qqXr3ZotwX1Bc/B/lHdivF1+8vjvorWbOnKWPDx3SXzdu0Gu/+KVJCYF7Y2lvb/OaHQIAAADwN6yRBgAAAAygSAMAAAAGUKQBAAAAAyjSAAAAgAEUaQAAAMAAijQAAABgAEUaAAAAMIAiDQAAABhAkQYAAAAMoEgDAAAABlCkAQAAAAOCzQ4AALg/Ojo6tGfPbmXv26f6+osKD++upKQkxbpcWrF8mTZs3KTQ0FBJ0tnyci1d+rr+smadoqKiJEklxcXasuU9nT59WmFhdg0fMVLPP/8jhYSEmHlbAOCzKNIAECC2bHlPOdnZ+sm8eXK5Bqrh8mWdO3/urq6tra3VypV/0qzZs/WzJUvU2Nio9ZkZWp+ZoSUvvfwtJwcA/8TSDgAIAM3Nzdr7z39ozty5SkycIKfTKdfAgZo4cdJdXb9r5w4lJCRo8uSnFR39fcXGujR//gLl5eWptbX1W04PAP6JGWkACADnqqvV1tam+Lh4Q9dXVJxVRUWF8vPzbznqldfbobq6OsXExNyfoAAQQCjSABAAbDbb145ZLBZJktfrvXms3ePpco7b7dakpCeVkpzypesdDsd9SgkAgYUiDQABwBkdLZvNphNFJzSxV68uY2Fh4ZKkhsuX1a1bN0nS2bPlXc7p1+9hnauukjM6+rsJDAABgDXSABAAbDabpqWm6d3Nm5WXd0C1tbUqKytTbk62nE6nIiMd2v7BNtXUnNexo0f19z27u1w/LTVVpaWlysxI19nyctXUnFdh4RFlZqSbdEcA4Pss7e1t3jufBgDwdR0dHdq1c4dycrJ16dJl9ejRXUlP/lBpadNVUlKijPS3VVNTqwED+is5ebLefPPPXR5/d+rUKW19f4vKykrl9UpOZy+NHjNW06c/Y/KdAYBvokgDAAAABrC0AwAAADCAIg0AAAAYQJEGAAAADKBIAwAAAAZQpAEAAAADKNIAAACAARRpAAAAwACKNAAAAGAARRoAAAAwgCINAAAAGECRBgAAAAygSAMAAAAG/A9FgVTOz+X5kgAAAABJRU5ErkJggg==",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plot_fit(free_predictions, \"Starting strength is a free parameter\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "12a108d0",
+ "metadata": {},
+ "source": [
+ "This is Figure 2 of Spicer et al.\n",
+ "With a free starting strength, the Rescorla–Wagner model fits the mean ratings almost perfectly ($R^2 = 1.00$), and its best-fitting parameters are the reported ones.\n",
+ "The best-fitting starting strength is 0.43, an intermediate value, as the uncertainty account predicts.\n",
+ "\n",
+ "Both models now separate B from Y: B loses associative strength on the B− trials of Stage 1, while Y enters Stage 2 at the starting strength.\n",
+ "Only the Rescorla–Wagner model also separates X from Y.\n",
+ "At the end of Stage 1, A predicts stomach ache, so the summed prediction of A and X on AX+ trials is close to the outcome, and X barely learns.\n",
+ "X ends the experiment close to where it started, near $\\beta$, so it is rated close to 5: the model is as uncertain about X as it was about any new food.\n",
+ "The Bush–Mosteller model has no summed error, so X and Y learn in the same way, and it predicts the same rating for both.\n",
+ "\n",
+ "The Bush–Mosteller parameters differ from the reported ones, because the model predicts almost the same ratings for a range of parameter values.\n",
+ "Our fit has an SSE of 0.060, and the reported parameters give 0.061.\n",
+ "\n",
+ "## Summary and next steps\n",
+ "\n",
+ "You fitted two learning models to the mean ratings of a blocking experiment, with a logistic response function and a starting strength that can be a free parameter, and recreated Figures 1 and 2 of Spicer et al. (2021).\n",
+ "When all the test cues are fitted, the Rescorla–Wagner model only explains the experiment better than the Bush–Mosteller model if novel cues start at an intermediate associative strength.\n",
+ "\n",
+ "To take this further:\n",
+ "\n",
+ "- Spicer et al. go on to fit the redundancy effect, and a manipulation of the outcome base rate (Jones et al., 2019), with the same models; their data and code are linked from the [supplementary materials](https://osf.io/7u6re/).\n",
+ "- To fit a model to each participant, rather than to the mean of a group, see [Tutorial 1: Build and fit your first model](../tutorials/first-model.ipynb).\n",
+ "\n",
+ "This example is based on a first version by Zhaotong Jia.\n",
+ "\n",
+ "## References\n",
+ "\n",
+ "Bush, R. R., & Mosteller, F. (1951). A mathematical model for simple learning. *Psychological Review*, 58(5), 313–323. https://doi.org/10.1037/h0054388\n",
+ "\n",
+ "Gluck, M. A., & Bower, G. H. (1988). From conditioning to category learning: An adaptive network model. *Journal of Experimental Psychology: General*, 117(3), 227–247. https://doi.org/10.1037/0096-3445.117.3.227\n",
+ "\n",
+ "Jones, P. M., Zaksaite, T., & Mitchell, C. J. (2019). Uncertainty and blocking in human causal learning. *Journal of Experimental Psychology: Animal Learning and Cognition*, 45(1), 111–124. https://doi.org/10.1037/xan0000185\n",
+ "\n",
+ "Rescorla, R. A., & Wagner, A. R. (1972). A theory of Pavlovian conditioning: Variations in the effectiveness of reinforcement and nonreinforcement. In A. H. Black & W. F. Prokasy (Eds.), *Classical conditioning II: Current research and theory* (pp. 64–99). Appleton-Century-Crofts.\n",
+ "\n",
+ "Spicer, S. G., Wills, A. J., Jones, P. M., Mitchell, C. J., & Dome, L. (2021). Representing uncertainty in the Rescorla-Wagner model: Blocking, the redundancy effect, and outcome base rate. *Open Journal of Experimental Psychology and Neuroscience*, 1, 14–21. https://doi.org/10.46221/ojepn.2021.6623"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "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.12.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/docs/examples/index.md b/docs/examples/index.md
index 9eca123..505df62 100644
--- a/docs/examples/index.md
+++ b/docs/examples/index.md
@@ -19,6 +19,18 @@ Why does a cue learned in compound with a pretrained cue gain so little strength
{bdg-success}`Beginner` {bdg-light}`15 min`
:::
+:::{grid-item-card} Fitting causal ratings in a blocking experiment
+:link: blocking-ratings
+:link-type: doc
+:img-top: ../_static/thumbnails/blocking-ratings.png
+
+{bdg-primary-line}`Associative learning`
+^^^
+Fit both learning rules to the ratings of a human blocking experiment, and recreate the finding of Spicer et al. (2021) that novel cues need an intermediate starting strength.
++++
+{bdg-warning}`Intermediate` {bdg-light}`25 min`
+:::
+
:::{grid-item-card} Model-based and model-free control in the two-step task
:link: two-step-task
:link-type: doc
@@ -49,6 +61,7 @@ Estimate how well confidence ratings track the accuracy of perceptual decisions.
:caption: Associative learning
blocking
+blocking-ratings
```
```{toctree}
diff --git a/test/datasets/test_datasets.py b/test/datasets/test_datasets.py
index 13a943d..0dadce6 100644
--- a/test/datasets/test_datasets.py
+++ b/test/datasets/test_datasets.py
@@ -2,7 +2,7 @@
import pytest
import pandas as pd
from unittest.mock import patch, mock_open
-from cpm.datasets import load_csv, load_bandit_data, load_risky_choices, load_two_step_data
+from cpm.datasets import load_csv, load_bandit_data, load_risky_choices, load_two_step_data, load_blocking_data
@pytest.fixture
@@ -94,3 +94,27 @@ def test_load_two_step_data_columns():
expected = {"ppt", "trial", "s1", "s2", "choice", "stimuli_left", "points", "timeout_1", "timeout_2"}
assert expected.issubset(data.columns)
assert data.ppt.nunique() > 1
+
+
+@patch("cpm.datasets.base.load_csv")
+def test_load_blocking_data(mock_load_csv):
+ # Mock load_csv to return a DataFrame
+ mock_load_csv.return_value = pd.DataFrame({"col1": [1, 3], "col2": [2, 4]})
+
+ # Call the function
+ result = load_blocking_data()
+
+ # Assertions
+ mock_load_csv.assert_called_once_with("blocking_ratings.csv")
+ assert isinstance(result, pd.DataFrame)
+ assert result.equals(pd.DataFrame({"col1": [1, 3], "col2": [2, 4]}))
+
+
+def test_load_blocking_data_means():
+ # The packaged file must reproduce the mean ratings reported by Spicer et al. (2021)
+ data = load_blocking_data()
+ assert {"ppt", "block", "trial", "cue", "rating", "rt"}.issubset(data.columns)
+ assert data.ppt.nunique() == 41
+ means = data.groupby(["ppt", "cue"]).rating.mean().unstack().mean()
+ reported = {"A": 9.84, "B": 1.47, "C": 0.35, "D": 0.81, "X": 4.97, "Y": 8.42}
+ assert means.round(2).to_dict() == reported