ampl_iv.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gprACN0E2fRM"
      },
      "source": [
        "# Simpel AMPL implementation of IV Extrapolation approach"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 219,
      "metadata": {
        "id": "YZotGiKU3RZY"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Please provide a valid license UUID. You can use a free https://ampl.com/ce license.\n"
          ]
        },
        {
          "data": {
            "application/vnd.jupyter.widget-view+json": {
              "model_id": "cdb5d32e335e478e9f9d51e96bf5d57f",
              "version_major": 2,
              "version_minor": 0
            },
            "text/plain": [
              "VBox(children=(Output(), HBox(children=(Text(value='', description='License UUID:', style=DescriptionStyle(des…"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Google Colab & Kaggle integration\n",
        "from amplpy import AMPL, tools\n",
        "ampl = tools.ampl_notebook(\n",
        "    modules=[\"highs\"], # modules to install\n",
        "    license_uuid=\"default\") # license to use"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BZ5JZD0FzXAg"
      },
      "source": [
        "# Diet Model\n",
        "\n",
        "Problem: minimize linear mapping over parameters $(\\theta_0, \\theta_1)$ subject to the constraint that parameters also give identified parameters.\n",
        "\n",
        "Sets:\n",
        "- `THETA`: Parameter space.\n",
        "- `THETA0`: parameter space of the model for $m_0$\n",
        "- `THETA1`: parameter space of the model for $m_1$\n",
        "- `IDENTIF`: set of identified IV-like estimands.\n",
        "\n",
        "Parameters:\n",
        "- `iv {IDENTIF}`: value of identified IV parameter(s).\n",
        "- `gamma {THETA}`: $\\Gamma^*_{dk}(b_{dk})$ for $d\\in\\{0,1\\}$ and $k=1,\\cdots,K_d$ - contribution to objective by each parameter.\n",
        "\n",
        "Variables:\n",
        "- `t0 {THETA0}`: value of parameters for $m_0$\n",
        "- `t1 {THETA1}`: value of parameters for $m_1$\n",
        "\n",
        "Objective:\n",
        "- `beta_star`: linear map giving $\\overline{\\beta}_{fd}^*$: sum over terms in `Gamma`.\n",
        "\n",
        "Constraints:\n",
        "- `iv {IV_identified}`: ensure that parameters result in identified IV-like estimands.\n",
        "\n",
        "The problem is then modeled as follows:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 220,
      "metadata": {
        "id": "4oBfIvMzbnvk"
      },
      "outputs": [],
      "source": [
        "%%ampl_eval\n",
        "reset;\n",
        "set NUTR;\n",
        "set FOOD;\n",
        "\n",
        "param cost {FOOD} > 0;\n",
        "param f_min {FOOD} >= 0;\n",
        "param f_max {j in FOOD} >= f_min[j];\n",
        "\n",
        "param n_min {NUTR} >= 0;\n",
        "param n_max {i in NUTR} >= n_min[i];\n",
        "\n",
        "param amt {NUTR,FOOD} >= 0;\n",
        "\n",
        "var Buy {j in FOOD} >= f_min[j], <= f_max[j];\n",
        "\n",
        "minimize Total_Cost:  sum {j in FOOD} cost[j] * Buy[j];\n",
        "\n",
        "subject to Diet {i in NUTR}:\n",
        "n_min[i] <= sum {j in FOOD} amt[i,j] * Buy[j] <= n_max[i];"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DVPF1DQixgSl"
      },
      "source": [
        "## Prepare `pandas.DataFrame` objects with all the data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 221,
      "metadata": {
        "id": "-L5H4Wzoxd36"
      },
      "outputs": [],
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "food_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"BEEF\", 3.59, 2, 10),\n",
        "        (\"CHK\", 2.59, 2, 10),\n",
        "        (\"FISH\", 2.29, 2, 10),\n",
        "        (\"HAM\", 2.89, 2, 10),\n",
        "        (\"MCH\", 1.89, 2, 10),\n",
        "        (\"MTL\", 1.99, 2, 10),\n",
        "        (\"SPG\", 1.99, 2, 10),\n",
        "        (\"TUR\", 2.49, 2, 10),\n",
        "    ],\n",
        "    columns=[\"FOOD\", \"cost\", \"f_min\", \"f_max\"],\n",
        ").set_index(\"FOOD\")\n",
        "\n",
        "# Create a pandas.DataFrame with data for n_min, n_max\n",
        "nutr_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"A\", 700, 20000),\n",
        "        (\"C\", 700, 20000),\n",
        "        (\"B1\", 700, 20000),\n",
        "        (\"B2\", 700, 20000),\n",
        "        (\"NA\", 0, 50000),\n",
        "        (\"CAL\", 16000, 24000),\n",
        "    ],\n",
        "    columns=[\"NUTR\", \"n_min\", \"n_max\"],\n",
        ").set_index(\"NUTR\")\n",
        "\n",
        "amt_df = pd.DataFrame(\n",
        "    np.matrix(\n",
        "        [\n",
        "            [60, 8, 8, 40, 15, 70, 25, 60],\n",
        "            [20, 0, 10, 40, 35, 30, 50, 20],\n",
        "            [10, 20, 15, 35, 15, 15, 25, 15],\n",
        "            [15, 20, 10, 10, 15, 15, 15, 10],\n",
        "            [928, 2180, 945, 278, 1182, 896, 1329, 1397],\n",
        "            [295, 770, 440, 430, 315, 400, 379, 450],\n",
        "        ]\n",
        "    ),\n",
        "    columns=food_df.index.tolist(),\n",
        "    index=nutr_df.index.tolist(),\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RtOMAi0KzobG"
      },
      "source": [
        "## Pandas dataframes\n",
        "\n",
        "Cost, minimum, and maximum amount of each food:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 222,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 332
        },
        "id": "dlA9U_2Vzg28",
        "outputId": "64ebdbf5-9980-4b2c-a03c-87624cbd4d63"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>cost</th>\n",
              "      <th>f_min</th>\n",
              "      <th>f_max</th>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>FOOD</th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>BEEF</th>\n",
              "      <td>3.59</td>\n",
              "      <td>2</td>\n",
              "      <td>10</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>CHK</th>\n",
              "      <td>2.59</td>\n",
              "      <td>2</td>\n",
              "      <td>10</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>FISH</th>\n",
              "      <td>2.29</td>\n",
              "      <td>2</td>\n",
              "      <td>10</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>HAM</th>\n",
              "      <td>2.89</td>\n",
              "      <td>2</td>\n",
              "      <td>10</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>MCH</th>\n",
              "      <td>1.89</td>\n",
              "      <td>2</td>\n",
              "      <td>10</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>MTL</th>\n",
              "      <td>1.99</td>\n",
              "      <td>2</td>\n",
              "      <td>10</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>SPG</th>\n",
              "      <td>1.99</td>\n",
              "      <td>2</td>\n",
              "      <td>10</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>TUR</th>\n",
              "      <td>2.49</td>\n",
              "      <td>2</td>\n",
              "      <td>10</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "      cost  f_min  f_max\n",
              "FOOD                    \n",
              "BEEF  3.59      2     10\n",
              "CHK   2.59      2     10\n",
              "FISH  2.29      2     10\n",
              "HAM   2.89      2     10\n",
              "MCH   1.89      2     10\n",
              "MTL   1.99      2     10\n",
              "SPG   1.99      2     10\n",
              "TUR   2.49      2     10"
            ]
          },
          "execution_count": 222,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "food_df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6QHGKGRl1Sgw"
      },
      "source": [
        "Minimum and maximum amounts of each nutrient:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 223,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 269
        },
        "id": "-nSsA-qBzjz5",
        "outputId": "5a1c70b5-7e69-4dba-c595-8be75a514c61"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>n_min</th>\n",
              "      <th>n_max</th>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>NUTR</th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>A</th>\n",
              "      <td>700</td>\n",
              "      <td>20000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>C</th>\n",
              "      <td>700</td>\n",
              "      <td>20000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>B1</th>\n",
              "      <td>700</td>\n",
              "      <td>20000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>B2</th>\n",
              "      <td>700</td>\n",
              "      <td>20000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>NA</th>\n",
              "      <td>0</td>\n",
              "      <td>50000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>CAL</th>\n",
              "      <td>16000</td>\n",
              "      <td>24000</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "      n_min  n_max\n",
              "NUTR              \n",
              "A       700  20000\n",
              "C       700  20000\n",
              "B1      700  20000\n",
              "B2      700  20000\n",
              "NA        0  50000\n",
              "CAL   16000  24000"
            ]
          },
          "execution_count": 223,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "nutr_df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mJ_3H1QA1hvY"
      },
      "source": [
        "Amount of each nutrient in each food:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 224,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 300
        },
        "id": "GQD03caEzlX4",
        "outputId": "254e49e0-f961-4f5f-8767-c3e855d15ae5"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>BEEF</th>\n",
              "      <th>CHK</th>\n",
              "      <th>FISH</th>\n",
              "      <th>HAM</th>\n",
              "      <th>MCH</th>\n",
              "      <th>MTL</th>\n",
              "      <th>SPG</th>\n",
              "      <th>TUR</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>A</th>\n",
              "      <td>60</td>\n",
              "      <td>8</td>\n",
              "      <td>8</td>\n",
              "      <td>40</td>\n",
              "      <td>15</td>\n",
              "      <td>70</td>\n",
              "      <td>25</td>\n",
              "      <td>60</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>C</th>\n",
              "      <td>20</td>\n",
              "      <td>0</td>\n",
              "      <td>10</td>\n",
              "      <td>40</td>\n",
              "      <td>35</td>\n",
              "      <td>30</td>\n",
              "      <td>50</td>\n",
              "      <td>20</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>B1</th>\n",
              "      <td>10</td>\n",
              "      <td>20</td>\n",
              "      <td>15</td>\n",
              "      <td>35</td>\n",
              "      <td>15</td>\n",
              "      <td>15</td>\n",
              "      <td>25</td>\n",
              "      <td>15</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>B2</th>\n",
              "      <td>15</td>\n",
              "      <td>20</td>\n",
              "      <td>10</td>\n",
              "      <td>10</td>\n",
              "      <td>15</td>\n",
              "      <td>15</td>\n",
              "      <td>15</td>\n",
              "      <td>10</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>NA</th>\n",
              "      <td>928</td>\n",
              "      <td>2180</td>\n",
              "      <td>945</td>\n",
              "      <td>278</td>\n",
              "      <td>1182</td>\n",
              "      <td>896</td>\n",
              "      <td>1329</td>\n",
              "      <td>1397</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>CAL</th>\n",
              "      <td>295</td>\n",
              "      <td>770</td>\n",
              "      <td>440</td>\n",
              "      <td>430</td>\n",
              "      <td>315</td>\n",
              "      <td>400</td>\n",
              "      <td>379</td>\n",
              "      <td>450</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "     BEEF   CHK  FISH  HAM   MCH  MTL   SPG   TUR\n",
              "A      60     8     8   40    15   70    25    60\n",
              "C      20     0    10   40    35   30    50    20\n",
              "B1     10    20    15   35    15   15    25    15\n",
              "B2     15    20    10   10    15   15    15    10\n",
              "NA    928  2180   945  278  1182  896  1329  1397\n",
              "CAL   295   770   440  430   315  400   379   450"
            ]
          },
          "execution_count": 224,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "amt_df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tYA-jeGUxujz"
      },
      "source": [
        "## Load the data from the `pandas.DataFrame` objects"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 225,
      "metadata": {
        "id": "6CPg9YzCxuTT"
      },
      "outputs": [],
      "source": [
        "# 1. Send the data from \"amt_df\" to AMPL and initialize the indexing set \"FOOD\"\n",
        "ampl.set_data(food_df, \"FOOD\")\n",
        "# 2. Send the data from \"nutr_df\" to AMPL and initialize the indexing set \"NUTR\"\n",
        "ampl.set_data(nutr_df, \"NUTR\")\n",
        "# 3. Set the values for the parameter \"amt\" using \"amt_df\"\n",
        "ampl.get_parameter(\"amt\").set_values(amt_df)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6NnyUc2qxzIM"
      },
      "source": [
        "# Solve with HiGHS"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 226,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "11ojqR5JyVqB",
        "outputId": "7cd7ae32-e9b8-45f0-d232-82f23832f1c8"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "HiGHS 1.5.1: HiGHS 1.5.1: optimal solution; objective 119.9897589\n",
            "4 simplex iterations\n",
            "0 barrier iterations\n",
            "Objective is: 119.98975893599335\n"
          ]
        }
      ],
      "source": [
        "# Specify the solver to use (e.g., HiGHS)\n",
        "ampl.option[\"solver\"] = \"highs\"\n",
        "# Solve\n",
        "ampl.solve()\n",
        "# Stop if the model was not solved\n",
        "assert ampl.get_value(\"solve_result\") == \"solved\"\n",
        "# Get objective entity by AMPL name\n",
        "totalcost = ampl.get_objective('Total_Cost')\n",
        "# Print it\n",
        "print(\"Objective is:\", totalcost.value())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rrwelmdwygur"
      },
      "source": [
        "## Reassign data for specific instances and resolve"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 227,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "plPCDhdhjivI",
        "outputId": "9f77f712-afcb-4bea-eaee-b1f0471bdaff"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Increased costs of beef and ham.\n",
            "HiGHS 1.5.1: HiGHS 1.5.1: optimal solution; objective 144.0120033\n",
            "0 simplex iterations\n",
            "0 barrier iterations\n",
            "New objective value: 144.01200332502077\n"
          ]
        }
      ],
      "source": [
        "# Reassign data - specific instances\n",
        "cost = ampl.get_parameter('cost')\n",
        "cost.set_values({'BEEF': 5.01, 'HAM': 4.55})\n",
        "print(\"Increased costs of beef and ham.\")\n",
        "\n",
        "# Resolve and display objective\n",
        "ampl.solve()\n",
        "# Stop if the model was not solved\n",
        "assert ampl.get_value(\"solve_result\") == \"solved\"\n",
        "print(\"New objective value:\", totalcost.value())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-Mj4AP9Qy60F"
      },
      "source": [
        "### Reassign data for all instances and resolve"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 228,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "j3j84nkVy6SD",
        "outputId": "ab75fc49-701a-4543-b7b6-51bdd606b775"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Updated all costs.\n",
            "HiGHS 1.5.1: HiGHS 1.5.1: optimal solution; objective 164.10625\n",
            "4 simplex iterations\n",
            "0 barrier iterations\n",
            "New objective value: 164.10625\n"
          ]
        }
      ],
      "source": [
        "# Reassign data - all instances\n",
        "cost.set_values(\n",
        "    {\n",
        "        \"BEEF\": 3,\n",
        "        \"CHK\": 5,\n",
        "        \"FISH\": 5,\n",
        "        \"HAM\": 6,\n",
        "        \"MCH\": 1,\n",
        "        \"MTL\": 2,\n",
        "        \"SPG\": 5.01,\n",
        "        \"TUR\": 4.55,\n",
        "    }\n",
        ")\n",
        "print(\"Updated all costs.\")\n",
        "\n",
        "# Resolve and display objective\n",
        "ampl.solve()\n",
        "# Stop if the model was not solved\n",
        "assert ampl.get_value(\"solve_result\") == \"solved\"\n",
        "print(\"New objective value:\", totalcost.value())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zmmo0fWhzFPi"
      },
      "source": [
        "## Get the values of the variable Buy in a pandas.DataFrame object"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 229,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 300
        },
        "id": "5BFTCHMyzEOg",
        "outputId": "3cb92934-9213-41ed-baa5-d05b84a1c54d"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>Buy.val</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>BEEF</th>\n",
              "      <td>10.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>CHK</th>\n",
              "      <td>2.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>FISH</th>\n",
              "      <td>2.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>HAM</th>\n",
              "      <td>6.851974</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>MCH</th>\n",
              "      <td>10.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>MTL</th>\n",
              "      <td>10.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>SPG</th>\n",
              "      <td>6.765351</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>TUR</th>\n",
              "      <td>2.000000</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "        Buy.val\n",
              "BEEF  10.000000\n",
              "CHK    2.000000\n",
              "FISH   2.000000\n",
              "HAM    6.851974\n",
              "MCH   10.000000\n",
              "MTL   10.000000\n",
              "SPG    6.765351\n",
              "TUR    2.000000"
            ]
          },
          "execution_count": 229,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "df = ampl.get_variable('Buy').get_values().to_pandas()\n",
        "df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E9p--jVyzH7D"
      },
      "source": [
        "## Get the values of an expression into a pandas.DataFrame object"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 230,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 300
        },
        "id": "X-Zj8cWsyzOD",
        "outputId": "6f9381fe-c329-4aab-91ed-69041e1d3170"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>100*Buy[j]/(Buy[j].ub)</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>BEEF</th>\n",
              "      <td>100.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>CHK</th>\n",
              "      <td>20.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>FISH</th>\n",
              "      <td>20.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>HAM</th>\n",
              "      <td>68.519737</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>MCH</th>\n",
              "      <td>100.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>MTL</th>\n",
              "      <td>100.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>SPG</th>\n",
              "      <td>67.653509</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>TUR</th>\n",
              "      <td>20.000000</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "      100*Buy[j]/(Buy[j].ub)\n",
              "BEEF              100.000000\n",
              "CHK                20.000000\n",
              "FISH               20.000000\n",
              "HAM                68.519737\n",
              "MCH               100.000000\n",
              "MTL               100.000000\n",
              "SPG                67.653509\n",
              "TUR                20.000000"
            ]
          },
          "execution_count": 230,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "df = ampl.get_data('{j in FOOD} 100*Buy[j]/Buy[j].ub').to_pandas()\n",
        "df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The n+1 Bernstein basis polynomials of degree $n$ are defined as \n",
        "$$ b_{v,n}(x) = {n\\choose v} x^v(1-x)^{n-v} $$\n",
        "for $v=0, \\cdots, n$.\n",
        "\n",
        "A linear combinations of Bernstein basis polynomoials\n",
        "$$ B_n(x) = \\sum_{v=0}^n\\beta_v b_{v,n}(x)$$\n",
        "is called a Bernstein polynomial."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 231,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Define some of the functions that we need\n",
        "from math import comb\n",
        "\n",
        "def bernstein_basis(n, v, x):\n",
        "    \"\"\"Bernstein polynomial basis of degree n and index v at point x\"\"\"\n",
        "    return comb(n, v) * x ** v * (1 - x) ** (n - v)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 232,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Weights for target\n",
        "def w_target_late(d, u, u_lo, u_hi):\n",
        "    \"\"\"Target weight w(d,u): LATE(u_lo, u_hi)\n",
        "    Args:\n",
        "        d (np.int): value of treatment\n",
        "        u (np.int): latent probability to get treatment\n",
        "        u_lo (np.float): lower bound of the LATE\n",
        "        u_hi (np.float): upper bound of the LATE\n",
        "    \"\"\"\n",
        "    # Return 1 divided by u_hi - u_lo if u_lo < u < u_hi, 0 otherwise\n",
        "    if u_lo < u < u_hi: w = 1 / (u_hi - u_lo)\n",
        "    else: w = 0\n",
        "\n",
        "    if d == 1: return w\n",
        "    else: return -w"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 233,
      "metadata": {},
      "outputs": [],
      "source": [
        "import scipy.integrate as integrate\n",
        "\n",
        "def gamma_star0_late(m0, weight, u_lo, u_hi):\n",
        "    \"\"\"Gamma*0: the first term of the gamma* function for binary IV, no covariates\n",
        "    Args:\n",
        "        m0 (function): MTE function for d = 0\n",
        "        weight (function): weights corresponding to target d = 0\n",
        "    \"\"\"\n",
        "    def func(u): return m0(u) * weight(0, u, u_lo, u_hi)\n",
        "    return integrate.quad(func, 0, 1)[0]\n",
        "\n",
        "def gamma_star1_late(m1, weight, u_lo, u_hi):\n",
        "    \"\"\"Gamma*0: the first term of the gamma* function for binary IV, no covariates\n",
        "    Args:\n",
        "        m1 (function): MTE function for d = 1\n",
        "        weight (function): weights corresponding to target d = 1\n",
        "    \"\"\"\n",
        "    def func(u): return m1(u) * weight(1, u, u_lo, u_hi)\n",
        "    return integrate.quad(func, 0, 1)[0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 234,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Compute the contributions of parameters for gamma_star(basis_function)\n",
        "def bern_2_0(u): return bernstein_basis(2, 0, u)\n",
        "def bern_2_1(u): return bernstein_basis(2, 1, u)\n",
        "def bern_2_2(u): return bernstein_basis(2, 2, u)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 235,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "-0.27999999999985575"
            ]
          },
          "execution_count": 235,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "gamma_star0_late(bern_2_0, w_target_late, u_lo=0.2, u_hi=0.8)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 236,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "0.27999999999985575"
            ]
          },
          "execution_count": 236,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "gamma_star1_late(bern_2_0, w_target_late, u_lo=0.2, u_hi=0.8)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 237,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Compute the LATEs from the example\n",
        "def m0_dgp(u): \n",
        "    return 0.6 * bern_2_0(u) + 0.4 * bern_2_1(u) + 0.3 * bern_2_2(u)\n",
        "\n",
        "def m1_dgp(u):\n",
        "    return 0.75 * bern_2_0(u) + 0.5 * bern_2_1(u) + 0.25 * bern_2_2(u)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 238,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGiCAYAAADA0E3hAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABcFUlEQVR4nO3dd3yV5cH/8c/Jyd5kJ8wkhB0CBoGALGUIPg5aC446sT+pVlSqLbgQ5/PYaq1WtNaBbQWp4qyIRBmyBBlh7ySEkZBFNmSe3x9XiMSgkpDknJN836/X9YJzc98nV3IR8uWaFpvNZkNERETECbjYuwIiIiIi50vBRURERJyGgouIiIg4DQUXERERcRoKLiIiIuI0FFxERETEaSi4iIiIiNNQcBERERGnoeAiIiIiTkPBRURERJxGo4PLN998w5VXXklUVBQWi4WPP/74Z59ZtWoViYmJeHp6EhMTw2uvvdaUuoqIiEg71+jgUlpaSkJCAn/729/O6/60tDQmTZrEiBEj2Lp1Kw899BAzZsxg8eLFja6siIiItG+WCzlk0WKx8NFHH3HNNdf86D1//OMf+fTTT9mzZ0/dtenTp7Nt2zbWr1/f1A8tIiIi7ZBrS3+A9evXM378+HrXJkyYwJtvvkllZSVubm4NnikvL6e8vLzudU1NDfn5+QQHB2OxWFq6yiIiItIMbDYbxcXFREVF4eLSPNNqWzy4ZGVlER4eXu9aeHg4VVVV5ObmEhkZ2eCZZ599lrlz57Z01URERKQVHDlyhE6dOjXLe7V4cAEa9JKcGZ36sd6T2bNnM3PmzLrXhYWFdOnShf379xMUFNRyFZWfVVlZyYoVKxgzZsw5e8uk9agtHIfawrGoPRxHfn4+PXr0wM/Pr9nes8WDS0REBFlZWfWuZWdn4+rqSnBw8Dmf8fDwwMPDo8H1oKCgH31GWkdlZSXe3t4EBwfrHwQ7U1s4DrWFY1F7OJ7mnObR4vu4JCUlkZycXO/asmXLGDRokP5CiYiISKM0OriUlJSQkpJCSkoKYJY7p6SkkJGRAZhhnptvvrnu/unTp3P48GFmzpzJnj17eOutt3jzzTd54IEHmuczEBERkXaj0UNFmzZtYsyYMXWvz8xFueWWW5g/fz6ZmZl1IQYgOjqaJUuWcP/99/PKK68QFRXFSy+9xC9/+ctmqL6IiIi0J40OLqNHj+antn6ZP39+g2ujRo1iy5Ytjf1QIiIi2Gw2qqqqqK6uPq/7KysrcXV15fTp0+f9jDSN1WrF1dW1VbcqaZVVRSIiIk1RUVFBZmYmZWVl5/2MzWYjIiKCI0eOaO+vVuDt7U1kZCTu7u6t8vEUXERExCHV1NSQlpaG1WolKioKd3f38woiNTU1lJSU4Ovr22ybnklDNpuNiooKcnJySEtLIy4urlW+3gouIiLikCoqKqipqaFz5854e3uf93M1NTVUVFTg6emp4NLCvLy8cHNz4/Dhw3Vf85amFhUREYem8OHYWrt99LdBREREnIaCi4iIiDgNBRcREZE2aP78+QQGBtq7Gs1OwUVERMQBZGRkcOWVV+Lj40NISAgzZsygoqLC3tVyOFpVJCIiYmfV1dVcccUVhIaGsmbNGvLy8rjllluw2Wy8/PLL9q6eQ1GPi4iIOAWbzUZZRdV5lVMV1ed97/mUn9ox/lxGjx7NPffcw3333UeHDh0IDw/n9ddfp7S0lNtuuw0/Pz9iY2P54osvAHP48O7du/n3v//NwIEDGTt2LM8//zz/+Mc/KCoqOq+POX/+fLp06YK3tzeTJ08mLy+vwT1PPfUUYWFh+Pn5cccddzBr1iwGDBhQ9+e33nor11xzDXPnziUsLAx/f3/uvPNOh+r5UY+LiIg4hVOV1fR57Eu7fOzdT0zA271xPzLfeecd/vCHP7Bx40YWLVrEb3/7Wz7++GMmT57MQw89xF/+8hduuukmMjIyWL9+Pf369SMqKqru+QkTJlBeXs7mzZvrnRF4Lhs2bOD222/nmWee4Re/+AVLly5lzpw59e559913efrpp5k3bx7Dhw/nvffe4/nnnyc6OrrefV9//TWenp6sWLGC9PR0brvtNkJCQnj66acb9fm3FPW4iIiItICEhAQeeeQR4uLimD17Nl5eXoSEhPCb3/yGuLg4HnvsMfLy8ti+fTtZWVmEh4fXe75Dhw64u7uTlZX1sx/rr3/9KxMmTGDWrFn06NGDGTNmMGHChHr3vPzyy0ybNo3bbruNHj168NhjjxEfH9/gvdzd3Xnrrbfo27cvV1xxBU888QQvvfQSNTU1F/YFaSbqcREREafg5WZl9xMTfva+mpoaiouK8fP3a7bN0bzcrI1+pn///nW/t1qtBAcH1wsKZ4JKdnY2wDmPM7DZbOd1zMGePXuYPHlyvWtJSUksXbq07vW+ffu466676t0zePBgli9fXu9aQkJCvZ2Kk5KSKCkp4ciRI3Tt2vVn69LSFFxERMQpWCyW8xquqampocrdire7q1133XVzc6v32mKx1Lt2JpDU1NQQERHBhg0b6t1/8uRJKisrG/TEnMv5zsH5YQhqzNwdRzmwUkNFIiIidpaUlMTOnTvJzMysu7Zs2TI8PDxITEz82ef79OnDt99+W+/aD1/37NmTjRs31ru2adOmBu+1bds2Tp06Ve99fH196dSp03l9Li1NwUVERMTOxo8fT58+fbjpppvYunUrX3/9NQ888AC/+c1v8Pf3/9nnZ8yYwdKlS3nuuefYv38/f/vb3+oNEwHcc889vPnmm7zzzjscOHCAp556iu3btzfoSamoqGDatGns3r2bL774gjlz5vC73/3OYc6McoxaiIiItGNWq5XPP/8cT09Phg8fzpQpU7jmmmv485//fF7PDx06lDfeeIOXX36ZAQMGsGzZMh555JF699x4443Mnj2bBx54gIsuuoi0tDRuvfXWBic6X3bZZcTFxTFy5EimTJnClVdeyeOPP95cn+oF0xwXERGRZrZy5coG19LT0xtcO3uOSZcuXfjvf//b5I95++23c/vtt9e79vvf/77e60cffZRHH3207vW4cePo3r17g/eaO3cuc+fObXJdWpKCi4iISDtQVlbGa6+9xoQJE7BarSxcuJCvvvqK5ORke1etUTRUJCIi4uAmTpyIr6/vOcszzzxzXu9hsVhYsmQJI0aMIDExkc8++4zFixczduzYFq5981KPi4iIiIN744036q30OVtQUNB5vYeXlxdfffXVT94zf/78xlat1Sm4iIiIOLiOHTvauwoOQ0NFIiIi4jQUXERERMRpKLiIiIiI01BwEREREaeh4CIiIiJOQ8FFRESkDZo/fz6BgYH2rkazU3ARERFxAPfeey+JiYl4eHgwYMAAe1fHYSm4iIiIOACbzcbtt9/O1KlT7V0Vh6bgIiIizsFmg4rS8yuVZed/7/mUsw5DPB+jR4/mnnvu4b777qNDhw6Eh4fz+uuvU1paym233Yafnx+xsbF88cUXdc+89NJL3H333cTExDTpyzN//ny6dOmCt7c3kydPJi8vr8E9Tz31FGFhYfj5+XHHHXcwa9aser07t956K9dccw1z584lLCwMf39/7rzzTioqKppUp5agnXNFRMQ5VJbBM1E/e5sLENjcH/uh4+Du06hH3nnnHf7whz+wceNGFi1axG9/+1s+/vhjJk+ezEMPPcRf/vIXbrrpJjIyMvD29r6g6m3YsIHbb7+dZ555hl/84hcsXbqUOXPm1Lvn3Xff5emnn2bevHkMHz6c9957j+eff57o6Oh693399dd4enqyYsUK0tPTue222wgJCeHpp5++oDo2F/W4iIiItICEhAQeeeQR4uLimD17Nl5eXoSEhPCb3/yGuLg4HnvsMfLy8ti+ffsFf6y//vWvTJgwgVmzZtGjRw9mzJjBhAkT6t3z8ssvM23aNG677TZ69OjBY489Rnx8fIP3cnd356233qJv375cccUVPPHEE7z00kvU1NRccD2bg3pcRETEObh5m56Pn1FTU0NRcTH+fn64uDTT/8/dGt8j0r9//7rfW61WgoOD6wWF8PBwALKzsy+4env27GHy5Mn1riUlJbF06dK61/v27eOuu+6qd8/gwYNZvnx5vWsJCQn1eoCSkpIoKSnhyJEjdO3a9YLreqEUXERExDlYLOc3XFNTA27V5t7mCi5N4ObmVu+1xWKpd81isQA0S0+G7Tzn4Jz5mI197lzP2ouGikRERJxcnz59+Pbbb+td++Hrnj17snHjxnrXNm3a1OC9tm3bxqlTp+q9j6+vL506dWrGGjedelxEREQcwMGDBykpKSErK4tTp06RkpICmFDi7u7+k8/OmDGDYcOG8dxzz3HNNdewbNmyesNEAPfccw+/+c1vGDRoEMOGDWPRokVs3769wSqmiooKpk2bxiOPPMLhw4eZM2cOv/vd75pv2O0CKbiIiIg4gDvuuINVq1bVvR44cCAAaWlpdOvW7SefHTp0KG+88QZz5szh8ccfZ+zYsTzyyCM8+eSTdffceOONpKam8sADD3D69GmmTJnCrbfe2qAX5rLLLiMuLo6RI0dSXl7Oddddx+OPP95sn+eFUnARERFpZitXrmxwLT09vcG1s+eYnOuZxrj99tu5/fbb6137/e9/X+/1o48+yqOPPlr3ety4cXTv3r3Be82dO5e5c+deUH1aioKLiIhIO1BWVsZrr73GhAkTsFqtLFy4kK+++ork5GR7V61RHGPASkRERH7UxIkT8fX1PWd55plnzus9LBYLS5YsYcSIESQmJvLZZ5+xePFixo4d28K1b17qcREREXFwb7zxRr2VPmcLCgo6r/fw8vLiq6+++sl75s+f39iqtToFFxEREQfXsWNHe1fBYWioSEREHFpjNkmT1tfa7aPgIiIiDunMLrNlZWV2ron8lDPt88OdgluKhopERMQhWa1WAgMD687y8fb2Pq9t52tqaqioqOD06dMOs2laW2Sz2SgrKyM7O5vAwECsVmurfFwFFxERcVgRERFA4w4itNlsnDp1Ci8vL4c5X6ctCwwMrGun1qDgIiIiDstisRAZGUlYWBiVlZXn9UxlZSXffPMNI0eObLXhi/bKzc2t1XpazlBwERERh2e1Ws/7B6TVaqWqqgpPT08FlzZIg38iIiLiNBRcRERExGkouIiIiIjTUHARERERp6HgIiIiIk5DwUVERESchoKLiIiIOA0FFxEREXEaCi4iIiLiNBRcRERExGkouIiIiIjTUHARERERp6HgIiIiIk5DwUVEREScRpOCy7x584iOjsbT05PExERWr179k/e/++67JCQk4O3tTWRkJLfddht5eXlNqrCIiIi0X40OLosWLeK+++7j4YcfZuvWrYwYMYKJEyeSkZFxzvvXrFnDzTffzLRp09i1axfvv/8+3333HXfccccFV15ERETal0YHlxdeeIFp06Zxxx130Lt3b1588UU6d+7Mq6++es77v/32W7p168aMGTOIjo7mkksu4c4772TTpk0XXHkRERFpX1wbc3NFRQWbN29m1qxZ9a6PHz+edevWnfOZYcOG8fDDD7NkyRImTpxIdnY2H3zwAVdcccWPfpzy8nLKy8vrXhcVFQFQWVlJZWVlY6oszezM11/tYH9qC8ehtnAsag/H0RJt0KjgkpubS3V1NeHh4fWuh4eHk5WVdc5nhg0bxrvvvsvUqVM5ffo0VVVVXHXVVbz88ss/+nGeffZZ5s6d2+D6ihUr8Pb2bkyVpYUkJyfbuwpSS23hONQWjkXtYX9lZWXN/p6NCi5nWCyWeq9tNluDa2fs3r2bGTNm8NhjjzFhwgQyMzN58MEHmT59Om+++eY5n5k9ezYzZ86se11UVETnzp0ZM2YMwcHBTamyNJPKykqSk5MZN24cbm5u9q5Ou6a2cBxqC8ei9nAcLbEQp1HBJSQkBKvV2qB3JTs7u0EvzBnPPvssw4cP58EHHwSgf//++Pj4MGLECJ566ikiIyMbPOPh4YGHh0eD625ubvpL6CDUFo5DbeE41BaORe1hfy3x9W/U5Fx3d3cSExMbdL8lJyczbNiwcz5TVlaGi0v9D2O1WgHTUyMiIiJyvhq9qmjmzJm88cYbvPXWW+zZs4f777+fjIwMpk+fDphhnptvvrnu/iuvvJIPP/yQV199ldTUVNauXcuMGTMYPHgwUVFRzfeZiIiISJvX6DkuU6dOJS8vjyeeeILMzEz69evHkiVL6Nq1KwCZmZn19nS59dZbKS4u5m9/+xu///3vCQwM5NJLL+X//u//mu+zEBERkXahSZNz77rrLu66665z/tn8+fMbXLvnnnu45557mvKhREREROrorCIRERFxGgouIiIi4jQUXERERMRpKLiIiIiI01BwEREREaeh4CIiIiJOQ8FFREREnIaCi4iIiDgNBRcRERFxGgouIiIi4jScKri4bHkHyvLtXQ0RERGxE6cKLtblj8Ofe8CiX8PeJVBdae8qiYiISCtq0iGL9mIL7QNFe2DPZ6Z4h0D/KTDgBoiIt3f1REREpIU5VXCpuuVzqMyEbQth+yIozYFv55kSHm8CTPyvwDfU3lUVERGRFuBUQ0UARPSDCU/DzD1w/SLofRVY3eHEDvhyNrzQCxZeD7s/haoKe9dWREREmpFT9bjUY3WDnpebUpYPOxdDygI4vgX2LTHFqwP0u9b0xEQNBIvF3rUWERGRC+C8weVs3kEw+Dem5OwzAWb7IijOhO/+YUpoL0i4HvpPBf9Ie9dYREREmsD5hop+TmhPGDcX7t8Fv15selxcPSFnL3w1B/7SB/79S9jxAVSesndtRUREpBHaRo/LubhYoftYU04Xwq6PTU/MkW/h4FemeARAv8kw4EbodLGGkkRERBxc2w0uZ/MMgMRbTMk7ZFYlbXsPCo/A5vmmBHc3Q0kJ10FAJ3vXWERERM6h7Q0V/ZzgWLj0Ebh3O9z8qQkrbt6QdxCWPwl/6Qf/vBq2LYKKMnvXVkRERM7SPnpczsXFBWJGmTLpT2b5dMoCOLwGUlea8rkv9LnGrErqkmSeEREREbtpv8HlbB5+MPBGU06mm96WbQvM71P+bUpg1++HkoKi7V1jERGRdkldCD/UoRuM/iPMSIHbvoCBvwZ3Pyg4DKv+F14aAG9Pgq3/hvJiO1dWRESkfVFw+TEWC3QdBle/Ag/sh1/8A2JGAxY4vBY+udsc+Pjh/zPDSjU1dq6wiIhI26ehovPh7m0Oc+w/BQqPmhVJ2xaaCb3bF5ni38kMIw24wUwAFhERkWanHpfGCugEIx+A322CaV9B4m1mP5iio7D6z/DyRfDGONj0NpwqsHdtRURE2hT1uDSVxQKdLzbl8v81ZyNtW2g2tju60ZQv/gi9/wcSboDYMWZTPBEREWkyBZfm4OYJ/X5hSnEWbP+PWVqds8cc/rhzMfhGQMJUE2LCetm7xiIiIk5JQ0XNzS8Chs+Au9bD/1sJg+8EryAoyYK1f4V5Q+D1MbDxH+ZUaxERETlv6nFpKRYLRA00ZfxTcOBLSFlofj2+xZQvH4Iel5sJvd3HgtXN3rUWERFxaAourcHVHXpfaUppLux43wwlZW2HPZ+a4hMK8VNgwPUQEW/vGouIiDgkBZfW5hMCQ39rStZOM6F3+yIozYFvXzElPN4EmPgp4Btq7xqLiIg4DM1xsaeIfjDhaZi5B65fBH2uBqs7nNhhhpFe6AULrjPnKFWV27u2IiIidqceF0dgdYOel5tSlm9WIW1bCMc2w/4vTPHqAP2uNfNhogaaOTQiIiLtjIKLo/EOgsG/MSVnn5kLs30RFGfCd/8wJbSXCTD9p5pVTCIiIu2EhoocWWhPGDcX7t8Fv15selxcPSFnLyQ/Bi/0hn9fa3poKk/bu7YiIiItTj0uzsDFapZLdx8Lpwth10dmafWRb+FgsikeAdBvMgy4ETpdrKEkERFpkxRcnI1nACTeakreITMXZtt7UHgENs83Jbi7OfCx/3UQ2Nm+9RUREWlGGipyZsGxcOkjcO92uPlTSLge3LzNqdXLn4IX4+Gdq2DbIqgotXdtRURELph6XNoCFxeIGWXKpD+Z5dPbFkL6akhbZcrnvtDnGjOpt0uSeUZERMTJKLi0NR5+MPBGU06mm96WbQvM71P+bUpgV9M7M+B66NDNzhUWERE5f/pvd1vWoRuM/iPMSIHbvoCBN4G7HxQchlX/C39NgLevgK3/hvJie9dWRETkZym4tAcWC3QdBlf/DR7YD7/4B8SMASxweA18cjf8uQd8eCekroSaGnvXWERE5Jw0VNTeuHtD/ymmFB41m9ulLDATere/Z0pAZ7O53YAbzARgERERB6Eel/YsoBOM+D38bhNM+woG3W72gyk8Aqv/DC9fBG+Mg01vw6kCe9dWREREPS6CGUrqfLEpE56FfUtML8yhr+HoRlO++CP0ugJL/FSwaShJRETsQ8FF6nPzhH6/MKU4C7b/x4SYnD2w60Ncd33IBNdAXLw2w8BfQ1gve9dYRETaEQ0VyY/zi4DhM+Cu9fD/VsLgO7F5BeFZVYB1/cswbwi8PgY2/sOcai0iItLCFFzk51ksEDUQJj1H1b072Rg9g5oeE8HFFY5vgSUPmFVJi26CfV9AdaW9aywiIm2UhoqkcazuZAYOonrSY7iUF8DODyDlXcjaAXs+NcUnFOKnmA3uIuLtXWMREWlDFFyk6XxDYehvTcnaaY4Z2L4ISnPg21dMCY83y6rjf2XuFxERuQAaKpLmEdEPJjwNM/fA9e9B76vAxQ1O7IAvZ8MLvWDh9bDnM6iqsHdtRUTESanHRZqX1Q16TjSlLB92Ljarko5vMcus9y0BryCIv9b0xEQOMHNoREREzoOCi7Qc7yAY/BtTsveawx63LYKSLNj4uimhvc1cmP5TzSomERGRn6ChImkdYb1g3BNw/y64cTH0+yVYPcz+MMmPwQu94d/Xmh6aytP2rq2IiDgo9bhI67K6QtxYU04VwK6PzFDS0Y1wMNkUjwCzAd6AG6DTxRpKEhGROgouYj9egTDoNlNyD5pVSdveg6KjsPltU4K7Q8L1kHCdOVtJRETaNQ0ViWMI6Q6XPQr37YCbP4H+14Gbtzm1evmT8Jd+8M+rzRyZijJ711ZEROxEPS7iWFxcIGa0KeV/ht2fQMpCOLwGUlea8rkf9L0aEm6ArsM0lCQi0o4ouIjj8vAzBzkO/DWcTDfDSNsWmt9v/bcpgV3NXJiE66BDNztXWEREWpqGisQ5dOgGo2fBjBS47QsTZtz9oOAwrHwW/poAb19hwkx5sb1rKyIiLUTBRZyLxWKGh65+BR7YD5NfN8NKWMxw0id3mwMfP7wTUldBTY29aywiIs2oScFl3rx5REdH4+npSWJiIqtXr/7J+8vLy3n44Yfp2rUrHh4exMbG8tZbbzWpwiJ13L0hYaqZzHv/Trj0UbMKqbIMtr8H/7wK/tofvn4S8g7Zu7YiItIMGj3HZdGiRdx3333MmzeP4cOH8/e//52JEyeye/duunTpcs5npkyZwokTJ3jzzTfp3r072dnZVFVVXXDlReoEdIKRD8CI38PR78zeMDs/hMIjsPrPpnQeYpZW951slmKLiIjTaXRweeGFF5g2bRp33HEHAC+++CJffvklr776Ks8++2yD+5cuXcqqVatITU0lKCgIgG7dujWpsvkl5QQHN+lRaS8sFug82JTL/xf2fW5WJR36Go5sMGXpLOh1hVmVFDsGXKz2rrWIiJynRgWXiooKNm/ezKxZs+pdHz9+POvWrTvnM59++imDBg3iueee41//+hc+Pj5cddVVPPnkk3h5eZ3zmfLycsrLy+teFxUVATD2xbUMiI1kTM9QxvQMoXeEHxYthW1VlZWV9X51bFboeZUpxZm47PwAl+3vYcndZ44W2LkYm28ENfG/oib+Ogjtae8KN4pztUXbprZwLGoPx9ESbdCo4JKbm0t1dTXh4eH1roeHh5OVlXXOZ1JTU1mzZg2enp589NFH5Obmctddd5Gfn/+j81yeffZZ5s6d2+C6Ddh2tJBtRwt58euDBLjb6NvBlB7+Ntz1H+dWk5ycbO8qNEEsdHqIwOA0OuetodPJ9biXZGFd/zLW9S9z0juGjKBLONZhKJWuvvau7HlzzrZom9QWjkXtYX9lZc2/YajFZrPZzvfm48eP07FjR9atW0dSUlLd9aeffpp//etf7N27t8Ez48ePZ/Xq1WRlZREQEADAhx9+yLXXXktpaek5e13O1ePSuXNnth9IZ3tODcv35bD2UB6nK79fMeLh6sLQmCDG9AhhdM9QOgaeuzdHLkxlZSXJycmMGzcONzc3e1fnwlSVYzmYbHphDiZjsVUDYLO6Y4ubQE3/67DFXApWx/w821RbODm1hWNReziOvLw8IiMjKSwsxN/fv1nes1E9LiEhIVit1ga9K9nZ2Q16Yc6IjIykY8eOdaEFoHfv3thsNo4ePUpcXFyDZzw8PPDw8GhwPaqDL/Hdg7kxKZrTldWsT81jxd5svt6TzbGCU6zan8uq/bnw3730DPdjTK8wLusdxsDOgbhatfK7Obm5uTn/PwhubhA/2ZSSHNjxPqQswHJiB5a9n+Gy9zPwCYX4KTDgeoiIt3eNz6lNtEUbobZwLGoP+2uJr3+jgou7uzuJiYkkJyczefLkuuvJyclcffXV53xm+PDhvP/++5SUlODra7rf9+/fj4uLC506Nf3QPE83K2N6hjGmZxhzr7Kx/0QJX+89wYq92Ww+fJJ9J4rZd6KY11YdIsDLjVE9Qrm0VxijeoTSwce9yR9X2ijfUEi6y5SsHWZC7/ZFUJoD375iSkS8mdAb/ytzv4iItLpGryqaOXMmN910E4MGDSIpKYnXX3+djIwMpk+fDsDs2bM5duwY//znPwG44YYbePLJJ7ntttuYO3cuubm5PPjgg9x+++0/Ojm3sSwWCz0j/OgZ4cddo7tTUFbBqv05LN+bzcp9ORSequTTbcf5dNtxXCwwoHMgl/YKY0yvMPpE+muCr9QXEQ+Xx8O4uXDwK0h5F/YtNYEmazYkPwpx483S6h6Xg6uCsIhIa2l0cJk6dSp5eXk88cQTZGZm0q9fP5YsWULXrl0ByMzMJCMjo+5+X19fkpOTueeeexg0aBDBwcFMmTKFp556qvk+ix8I9Hbn6gEduXpAR6qqa9h6pIDle7NZsTebvVnFbMkoYEtGAX9etp9wfw/Tc9MrjEu6h+DjoeObpJbVDXpONKUs36xESnkXjm+FfUtM8QqC+GvNeUmRA3Tgo4hIC2vU5Fx7KSoqIiAggNzcXIIvcCOX4wWnWLHPhJi1B/M4VVld92fuVhcGRwcxumcoY3qFERPio96YH6isrGTJkiVMmjSp/Y4dZ+81AWb7f6DkrPleob3NXJj+U8EvosWrobZwHGoLx6L2cBx5eXmEhITYb3JuWxAV6MWNQ7py45CunK6sZkNaPiv2ZrNiXzaH88pYczCXNQdzeerzPXQJ8mZMz1BG9wojKSYYTzettxYgrBeMfxIumwOpK2HbAtjzX8jZA8mPwVePQ+xlJsT0vALcPO1dYxGRNqPdBZezebpZGdUjlFE9Qplj60Nabikr9uWwYm82G9LyyMgv4531h3ln/WE8XF1Iig2umxDcJdjb3tUXe7O6QtxYU04VwK6PYNtCszvvwWRTPAOg7y9gwI3QaZCGkkRELlC7Di5ns1gsxIT6EhPqy7RLoiktr2LNwVxW7sth5b5sMgtP1/4+hznsIibEh9E9wxjdM5TB0UHqjWnvvAJh0G2m5B40AWbbe1B0FDa/bUpwdzOhN+E6c7aSiIg0moLLj/DxcGVC3wgm9I3AZrOx70QxK/eZlUqbD58kNbeU1Nw03lqbhpeblWGxwYzuGcronmF0DlJvTLsW0h0uexTGPAzp35il1Xs+hbyDsPxJWP4UxIwyvTC9/secci0iIudFweU8WCwWekX40yvCn+mjYik6XcnaA7W9MfuzOVFUztd7s/l6bzawi5hQH0b3UG9Mu+fiAjGjTSn/M+z+xISYw2vM3JjUleDuB32vNiGmS5KGkkREfoaCSxP4e7oxMT6SifGR2Gw29mQWs3K/2TNm8+GTpOaUkppjemM83VxIiglmdE+z+V23EB97V1/swcMPBv7alPw0s7ldygIoOAxb/21Kh25mg7uE66BDV3vXWETEISm4XCCLxUKfKH/6RPlz1+ju5+yNWbEvhxX7cgDoFuxtJgT3DCUpJgQvnQzZ/gRFw+hZMPIPkLHeBJjdH8PJdFj5jCndRpj5MH2uBg/nOfBRRKSlKbg0sx/2xuzNMnNjVu3PZlP6SdLzykivXank7urCkOigupVN3cN8tW9Me+LiAt2GmzLpObOkOuVdSPsG0lebsuQBE14SrjdhxkVnbolI+6bg0oIsFgu9I/3pHenPb0fHUny6knWH8li1P4dV+3I4VnCK1QdyWX3A7BsTFeDJqJ4mxAzrHoK/pzZOajfcfSBhqikFR2D7e2Y+TP6h2hVKCyGgsxlGSrgegmPtXWMREbtQcGlFfp5u9VYqHcopqe2NyWFDWj7HC0+zcOMRFm48gtXFwkVdAhkZZ4aV+kUF4OKi3ph2IbAzjHwQRjwAR78zvTA7P4LCI/DNn0zpPBRL/BRcq7UiSUTaFwUXO7FYLHQP86N7mB93jIjhVEU136bl8U1tb0xqbinfpZ/ku/STPJ+8nyAfdy7pHsLIHqGMjAshzF+7sbZ5Fgt0HmzK5f9rzkZKWQCHlsORb3E98i2XW9ywVC+DgTdCzBhw0ZwpEWnbFFwchJe7tW5XXq6EI/llfHPAhJh1h/LIL62oO+EaoHekPyN7hDAqLpTEbh3wcNUPrDbNzQv6/dKUokzY8R9sW9/FmrsPdn1oil8U9J9iDnwM7WnvGouItAgFFwfVOci77kylyuoathw+yTcHcvhmfy47jhWyJ7OIPZlF/H1VKl5uVobGBJnemB6hOhyyrfOPhOH3UnXxb1m3+BUu8TuKdddiKD4Oa180pWOimQvT75fgHWTvGouINBsFFyfgZnVhSEwwQ2KCeXAC5JWUs+ZgLqv257D6QC45xfWXXHcM9GJkjxBGxoUyLDaEAG9N8m2TLBYKvGOomfA7rJc/A/u/NJN4DyyDY5tN+fIh6DnR7A/T/TKw6u+CiDg3BRcnFOzrwdUDOnL1gI51G+B9cyCH1Qdy+C7tJMcKTtVN8nWxQELnQEbEhTKqRwgJnQJxtWpJbZvj6gF9rjKlJAd2vG9Orc7aYXbs3f0J+IRCfO1QUkQ/e9dYRKRJFFyc3Nkb4E0fFUtZRRUb0vL5prY35mB2CVszCtiaUcBLXx/Az8OVpNhgRtRO8u0arJ182xzfUEi6y5SsHWZZ9fZFUJoD375iSkS8OWYg/lfgE2LvGouInDcFlzbG2931+0m+wPGCU6w+YELMmoO5FJRVsmz3CZbtPgFAlyBvLokLYWRcCEmxIQR4aSihTYmIh8vjYdxcOPiVWZW07wsTaJbOgmWPQNx40wsTNwFc3e1dYxGRn6Tg0sZFBXox9eIuTL24C9U1NnYdL2T1ATM/Zsvhk2Tkl7FgQwYLNmR8P6zUPYRL4kIZ2CUQNw0rtQ1WNzPXpedEKMuHnYvN/jDHt5pl1vuWgFeQ6YEZcD1EDtCBjyLikBRc2hGri4X+nQLp3ymQu8d0p6S8ig2pebW79+ZwKKf0+2Gl5QfxcbcyNCaYS+JCGBEXQmyozsxpE7yDYPBvTMneU7sz7yIoyYKNfzclrI9ZldR/CvhF2LvGIiJ1FFzaMV8PVy7rHc5lvcMBM6y05qA5gmDtwVzySyv4em82X+/NBiAywJOkmCB8ii0MLiknsoOGlZxeWG8Y9wRc+hikrjQTevf8F7J3Q/Kj8NUc6D7WhJiek8BNGx+KiH0puEidqEAvpgzqzJRBnampsbE7s4g1B3NZcyCXjen5ZBae5sOtxwEr//q/VfSK8OOS7iFcEhfC4OggvN3118lpWV0hbqwppwpg10dmPszRjWZ59YFl4Blg9oVJuAE6DdJQkojYhX7SyDm5uFjo1zGAfh0DmD4qltOV1XyXns+qfdl8sSWNY2UW9mYVszermDfWpOFmtXBRlw5c0j2E4XEh9O8YoGXXzsorEAbdZkrugdqhpPeg6BhsesuU4DgzF6b/dRDQ0d41FpF2RMFFzounm5URcaEM7RZIfPVBhoway3eHC1lTu1rpWMEpNqTlsyEtn+eT9+Pn6UpSTDDDu4cwvHsIsaHazdcphcTBZY/BmIchfbXphdn9KeQdgK+fgK+fhJjRZlVSr/8Bdx36KCItS8FFmiTYx50rE6K4MiEKm81Gel4Zaw7msvZALusO5VJ0uqresusIf8/aEGPCTLgOiXQuLlYTUGJGw6Q/mw3tti2Ew2shdYUp7n7Q9xqzP0yXoRpKEpEWoeAiF8xisRAd4kN0iA83De1KdY2NnccKTZA5mMumwyfJKjrN4i1HWbzlKADdw3wZHmtCzJCYYO0f40w8/eGim0zJTzPDSNsWQsFh2PovUzpEmwm9CddBh672rrGItCEKLtLsrC4WEjoHktDZLLs+XVnNpvSTrD1kgsyOY4UczC7hYHYJ76w/jIsF4jsFMiw2mOGxIQzq1gFPN5127RSComHMbBj1R8hYVzuU9AmcTIOVz5jSbYQZSup9FXhoSb2IXBgFF2lxnm5WLokzq48ACsoq+DY1j7UH81h7KJfUnFK2HSlg25ECXl15CHdXFxK7dGBYbDDDuofQv1OANsJzdC4u0O0SUyb9CfZ8ZkJM2jdmbkz6avj8AXOWUsL1Jsy4qE1FpPEUXKTVBXq7c3m/SC7vFwlAZuEp1tWGmHUH88gqOs361DzWp+bxfPJ+fD1cGRwdZIJMbAi9IvxwcdH8CYfl7mOGiBKug4Ij5pykbQsh72DtCqWFENC59p7rITjW3jUWESei4CJ2FxngxS8TO/HLxE7YbDYO5ZSyrjbErE/No/BUJcv3ZrO8diO8Dt5uJMUGkxQbwrDYYGJCtGLJYQV2hpEPwIjfw9HvzDEDOz+CwiPwzZ9M6TzULK3uO9nsFSMi8hMUXMShWCwWuof50j3Ml5uTulFdY2NPZhHrDuWy9mAe36Xnc7KskiU7sliyIwuAcH8PhsWGkBQbzLDYYDp10JJch2OxQOfBplz+v7D3c9Pzcmg5HPnWlC/+aJZUD7geYsaYlUwiIj+g4CIOzXrWRnj/b2QsFVU1bD9awLpDeaw7lMuWwwWcKCrno63H+GjrMcCceD0sNtj0ysQEE6al147FzQvirzWlKPP7oaScvbDzA1P8IqH/VDOpN7SnvWssIg5EwUWcirurC4O6BTGoWxAzLovjdGU1mw+fZH1tkNl2tJCM/DIy8st477sjAMSE+pggExPC0Jgggn097PxZSB3/SLjkPhh+LxzfAikLTXApzoS1L5rSMdHMhen3S3NApIi0awou4tQ83ax1u/NCT0rKq/guLd9M7j2Ux87jhaTmlJKaU8q/v80AoGe4H0mxwQyNCWJIdDAdfNzt+0mIGUrqmGjKhKdh/1ITYg4sg2ObTfnyIeg50WxwF3uZOV9JRNodfedLm+Lr4cqYXmGM6RUGQGFZJRvS8lh3KI9vU/PYm1XMvhOmzF+XjsUCvSL8GRoTRFJMMEOigwnw1mZ4duXqAX2uNqUkG3a8b0LMiR1mj5jdn4BPGPSfYoaSwvvau8Yi0ooUXKRNC/B2Y3zfCMb3jQAgr6ScDWn5rD9kViwdzC5hT2YRezKLeHutCTJ9Iv0ZGhPM0JhgBkcHaVdfe/INg6S7TcncbubCbP8PlGbD+r+ZEtHf9MLEXws+IfausYi0MAUXaVeCfT2YFB/JpHizh0xOcTnfppremPWpeaTmlLLreBG7jhfx5po0LBboG+XP0GgTZC5WkLGfyP6mjHsCDiTDtgWwbylkbYel22HZwxA3wfTCxI0HVw0BirRFCi7SroX6edQdFglwouh0bZDJZ0NqHqm5pew8VsTOY0W8URtk6vXIdAvS0FJrs7pBr0mmlObBzsUmxBzfCvs+N8U7GPpda0JMZIIOfBRpQxRcRM4S7u/J1QM6cvWAjsC5g8wPe2R6R/gzpHai75DoIE32bU0+wTDk/5mSvcccM7B9EZScgI1/NyWsjwkw8VPAL9zeNRaRC6TgIvITfhhksotO823tHJkNaWZoaXdmEbtr58gA9IrwY0h0EENq58iEaPl16wjrDeOfhMvmQOoKs0vv3iWQvRuWPQLJc6D7WLPBXY+J4Kb9fUSckYKLSCOE+XtyVUIUV9UOLWUXn2ZDaj4b0vLYkJrPgewS9mYVszermHfWHwYgNtSHITGmN2ZoTDDh2hCvZVldIW6cKadOwq6PzKqkoxvhwJemeAZ8P5TUMVFDSSJORMFF5AKE+XnWmyOTW1LOd2n5bEjLr1t+fSinlEM5pSzYYPaR6RrszZDoIAbXDi116uCls5ZailcHGHS7KbkHag95XARFR2HTm6aE9DAb3CVcB/5R9q6xiPwMBReRZhTi68HE+Egm1q5aOllawXfpJshsTMtn1/FCDueVcTivjP9sOgpAVIAng2uHli7uFkRsqA6NbBEhcXDZYzDmEUj/xsyH2f0p5O6Hr+fC109A7BhIuAF6XQHuOvNKxBEpuIi0oA4+7vX2kSk6Xcnm9JO1QSaP7UcLOV54mo9TjvNxynEAgn3cGRwdxMXdghgcHUTvSH+sLgoyzcbFBWJGm3LF82ZDu5QFcHitOfTx0HJw94O+15j9YboM1VCSiANRcBFpRf6ebvV29i2rqGJrRgEbUvPYmJ7P1owC8kor+GJnFl/sNKdf+3m4ktitAxd3C2JIdBDxnQLwcNXJyc3Cww8G/tqU/DTY9p4ZTio4DFv/ZUqH6O+Hkjp0tXeNRdo9BRcRO/J2dz3rrCUor6pmx9HCuqGlzYdPUlxexcp9OazclwOAh6sLCZ0DSewSgK3AwsjyKjq4aS+ZCxYUDWNmw6g/QsY6M6F398dwMg1WPmNKtxEmxPS5Gjx87V1jkXZJwUXEgXi4WutOv757DFTX2NiTWcTG2iDzXXo+eaUVda/ByutPL6dPlL8ZWqp9NtRPS7CbzMUFul1iyqTnYM9nZigp7RtIX23Kkgehz1VmVVLXS+xdY5F2RcFFxIFZXSz06xhAv44B3H5JNDabjUM5pWbC76Fcvtl7nPxyS93uvmf2kokJ8WFQ7fDSxd2C6BrsrQm/TeHuY4aIEq6DgiOw/T3TE5N/qHaF0kII6IJLv1/hUx5h79qKtAsKLiJOxGKx0D3Ml+5hvlw7MJIlS44wcPilbD1axKb0k3yXns++E8Wk5paSmltat3Ip1M+Di7t1YFBXE2R6R/rhanWx82fjZAI7w8gHYcQDcGSjOWZg54dQmIF17fOMBWqKPoCBN5qJvZ4B9q6xSJuk4CLi5CIDPOkS4le3u29hWSWbM/LZmHaSTen5bD9aSE5xOUt2ZLFkh5nw6+1u5aIuHUjsanplBnYJxMdD/xycF4sFugwx5fL/hb2fU7N1AZbU5bgc3QBHN8AXf4Be/2OGkmJGg4smU4s0F/1LJdLGBHi7cWmvcC7tZc7lOV1ZzfajhXyXbubIbD58kuLTVaw5mMuag7mAGZLqE+lfF2QGdeugHX7Ph5sXxF9Lda+rWf7Ju4wNy8O6fRHk7oOdH5jiFwUJU83+MKE97F1jEaen4CLSxnm6WRkcbfaEAaipsbE/u5jv0k2PzHdp+RwvPM2OY4XsOFbI/HXpAHQO8mJQ1yASu3ZgULcO9Ajzw0X7yfyo024dqEm6EeuImXB8i5kLs+N9KD4Oa/5iSsdE0wvT75dmV18RaTQFF5F2xsXFQq8If3pF+HPTULMvyfGCU2w6fJLN6fl8l36SvVlFHMk/xZH8Y3y09RgAfp6uXNSlA4O6miGmAV0C8XbXPyENWCwmoHRMhAlPw/6lZlXSgWQ4ttmUpbOh5yQTYmIvM+crich50XeLiBAV6MVVgV51h0cWn65ka0aBCTOHzcZ4xaerWLU/h1X7zX4yZw8vnSlRgV72/DQcj6uH2fOlz9VQkm16YFIWwImdZo+Y3R+DbzjE/8qEmPC+9q6xiMNTcBGRBvw83RjZI5SRPUIBqKquYU9mMZsP57Pp8Em2HD55zuGlqABPLjoryPSO9MdNq5cM3zBIutuUzO0mwOz4D5ScgPV/MyUywcyFif8V+ATbu8YiDknBRUR+lqvVhfhOAcR3CuDW4dHA98NLWw6fZPPhk+zOLOJ44WmOb8/kv9szAfB0cyGhU2BdkBnYpQNBPu72/FQcQ2R/U8Y/aYaQUt6F/V9C5jZTlj0MPS43u/TGjQdXfc1EzlBwEZEm+eHwUml5FduOFrA5/SRbMk6yJaOAwlOVbEgzp2OfERPiw0VdO3BRlw5c1DWQuDC/9nuIpNUNek0ypTQPdi42ISYzBfb+1xTvYOh3rRlKikzQgY/S7im4iEiz8PFwZVhsCMNizblLNTU2UnNL2FzbI7Mlo4CD2SV1m+N9sNlsjufn4UpC58DaMBPIwM4dCPBuh2cv+QTDkP9nSvYeM5S0/T9QkgUb/25KWB8TYOKngF+4vWssYhcKLiLSIlxcLHQP86N7mB9TL+4CQEFZBVszCtiSYcJMypECisvr7ykD0D3Ml4F1YaYDcWG+7WspdlhvM4x02RxIXWF6YfYugezdsOwRSJ4D3S8zIabHRHDTnjvSfii4iEirCfR2Z0yvMMb0CgPMpN99J4rZklHA1sNmiCk9r4yD2SUczC7h/R/2ynQJZGCXDgzsEkigdzuY92F1hbhxppw6Cbs+MvvDHN0IB5aZ4hlg9oUZcKNZgq2hJGnjFFxExG5crS70jQqgb1RA3Z4yeSXldb0yWzJOsu1I4Tl7ZWJCfBhwJsh0DqRXRBs/f8mrAwy63ZTcA7VDSYug6BhsesuUkB5mQm//qRDQ0d41FmkRCi4i4lCCfT0Y2yecsX3MHI4zvTJnwkxKRkHdPJnU3FI+3GI2yPNysxLfKYCBtfNkBnYJbLvHFoTEwdg5cOkjkPaNCTF7PoPc/fD1XPj6CYgdY5ZW97oC3L3tXWORZqPgIiIO7exemV/X9sqcLK0g5agZXtp6pICUDDNXZmNaPhvPWsEUFeBpemU6m51++0UF4OXehg48dLGagBI7Bk4XmQ3tUhZCxjo4tNwUD39zWnXCDdBlqIaSxOkpuIiI0+ng486YnmGM6WnmytTU2DiUU8LWIwVszShga8ZJ9p8oNvvKnHUqttXFQq8IPwZ0DmRA50AGdgkkJqSNTPz19IeLbjYlPxW2vQfbFkJBBmz5pylBMWYoKeE6COxi7xqLNEmTgsu8efP405/+RGZmJn379uXFF19kxIgRP/vc2rVrGTVqFP369SMlJaUpH1pEpAEXFwtx4X7EhfsxZVBnwOwrs/1oISlHTJBJOVJAdnE5u44Xset4Ee9uyADMGUwJnUyQSagNNKF+Hvb8dC5cUAyMeQhGzYLDa02A2fWxCTQrnjal2wizKqn3VeDha+8ai5y3RgeXRYsWcd999zFv3jyGDx/O3//+dyZOnMju3bvp0uXHE3xhYSE333wzl112GSdOnLigSouI/BwfD1eSYoNJijVb59tsNjILT5NypMCUjAK2HzNnMP1w4m/HQK/aIBNAQqdA+nUMwMfDCTuoXVwgeoQpE58z82C2LYC01ZBeWz5/wJylNOB66HqJeUbEgTX6O/GFF15g2rRp3HHHHQC8+OKLfPnll7z66qs8++yzP/rcnXfeyQ033IDVauXjjz/+yY9RXl5OeXl53euioiIAKisrqaysbGyVpRmd+fqrHexPbdF4oT6ujOsVwrheZpO8quoa9meXsO1oIduPFrHtaAEHc0o5VnCKYwWn+HyHObrAxQJxYb707xRA/44B9O/kT1yYb905TE7RFi4e0PdaUwqP4LLjP7hsfw/LyTQTZrYtwBbQmZr4KdT0vw46RNu7xk3mFO3RTrREG1hsNpvtfG+uqKjA29ub999/n8mTJ9ddv/fee0lJSWHVqlXnfO7tt99m3rx5rF+/nqeeeoqPP/74J4eKHn/8cebOndvg+oIFC/D21ux4EWk5p6vgSKmF9BLIKLFwuMRCYUXDOTBuLjY6+UAXXxtdfGx09bUR4ulkc19tNoJKD9A5fw0dT27AreZU3R/l+fQgI+gSjncYQpVVp35L05SVlXHDDTdQWFiIv79/s7xno3pccnNzqa6uJjy8/lbT4eHhZGVlnfOZAwcOMGvWLFavXo2r6/l9uNmzZzNz5sy610VFRXTu3JkxY8YQHKwTU+2psrKS5ORkxo0bh5tbO9yW3YGoLVpPVtFpdhwtYvuxQrYfK2THsSKKT1eRVgxpxd8nFX9PV/p3CiC+oz/9O5pfnWNJ9n1QeYqq/Utw2b4IS9pKgkv3E1y6nwGZC7H1nERN/+uxdRthVjI5OH1vOI68vLxmf88mDdpafvBfCpvN1uAaQHV1NTfccANz586lR48e5/3+Hh4eeHg0nBzn5uamv4QOQm3hONQWLa9zsBudg/2YlGA2daupsZGWV8q2IwVsq50zs+tYAUWnq1hzMI81B7//xzrc34P+nQJJ6BRA/06B9O8U4Ji7/rq5wYDrTCnKNJvbpSzAkrsPy67FuOxaDH5RkDDVLK0OPf9/0+1F3xv21xJf/0YFl5CQEKxWa4Pelezs7Aa9MADFxcVs2rSJrVu38rvf/Q6AmpoabDYbrq6uLFu2jEsvvfQCqi8i0vpcXCzEhvoSG+rLLy7qRGVlJZ/+dwkxAy9hV1YJ244UsONYIftPFHOiqJzk3SdI3v39ooQuQd5mvkxtmOnXMQBfR5r86x8Jl9wHw++F41vM3jA73ofi47DmL6Z0HGQm9Pb7pdnVV6SVNOo7xd3dncTERJKTk+vNcUlOTubqq69ucL+/vz87duyod23evHksX76cDz74gOho5538JSJyNlcX6NfRn4Hdgus2yiurqGLX8SK2HSlg+9FCth8tID2vjIx8U/673Uz+tVjMEQZnemTiOwbQJ8ofb3c7hxmLxZx/1DERJjwN+5eaXXoPJMOxTaYsnQ09J5qzkmIvM+cribSgRv8NmzlzJjfddBODBg0iKSmJ119/nYyMDKZPnw6Y+SnHjh3jn//8Jy4uLvTr16/e82FhYXh6eja4LiLS1ni7u3JxtyAu7hZUd62wrJIdxwrZdrSA7UcL2HmsiGMFpziUU8qhnFI+2mqOMDArmfyIr+2Z6dcxgD6R/ni62WmOiauHWTbd52ooyTY9MCkL4MRO2P2JKT5h0H+K2R8mvK996iltXqODy9SpU8nLy+OJJ54gMzOTfv36sWTJErp2Nf/DyMzMJCMjo9krKiLSFgR4u3FJXAiXxIXUXcstKWfH0cK6XpntxwrJKS5n34li9p0o5oPaU7KtLhbiwnyJ7xhAfG3PTG97hBnfMEi625TM7SbA7PgPlGbD+r+ZEtHf9MLE/wp8tKhCmk+jlkPbS1FREQEBAeTm5mpVkZ1VVlayZMkSJk2apElvdqa2cBwt0RYnik6bMHOskB1HzZyZ3JKKBvf9MMzYrWemutIMIaW8C/u/hJra/TtcXKHH5eaogbjx4NryE5P1veE48vLyCAkJsd9yaBERaR3h/p6E9/GsOyXbZrPVLssuZMcxU3bWhpm9WcXszSrm/bN6ZrqH+tKvdkl2v9aYM2N1g16TTCnNg50fmJ6YzBTY+19TvINND0zC9RCZ4GSb3oijUHAREXECFouFyAAvIgO8GN83Avg+zGw/akKMCTNF5JZ8P8y0eMuZ5yE21Jd+USbInAkz/p4t0CPhEwxD7jTlxG6zM+/2/0DJCdjwmilhfc2qpPgp4NdwVarIj1FwERFxUmeHmQlnhZkTReVnBZlCdh4v5ERROQezSziYXcLHKcfr3qNbsDd9owLo29GfflEB9I3yJ9i3GQ+ZDO8D45+Cyx6HQ8tNiNm7BLJ3wbJHIHkOdB9rQkyPieDmDBv2iT0puIiItCEWi4WIAE8iAr4fZgLILj7NruNF7Kwdatp13KxmSs8rIz2vrO5cJoDIAE/6RvmbQBPlT9+OAUQFeJ5zo9HzZnWFHuNNOXUSdn5oTq0++h0c+NIUz0CzL8yAG8wSbA0lyTkouIiItANhfp6E9fRkTM+wumsnSyvYdbyIXcdNmNl9vIjU3FIyC0+TWXiar/Zk193bwduNPmeHmSh/okN8sbo0IVx4dYCLp5mSe8DMhdm+CIqOwaY3TQnpYebCJFwH/lHN8SWQNkLBRUSknerg495gaXbx6Ur2ZBaz67iZL7PreCEHs0s4WVbJ2oN5rD3rOANPNxd6RfjX9c70ifKnZ7gfXu6NWNEUEgdj58Clj0DaKrNL757PIHc/fD0Xlj8JMaPN0upeV4CbDnxs7xRcRESkjp+nG4Ojgxgc/f2meacrqzlwosSEmeOmZ2ZPZjGnKqtJqT2r6QyX2knAfaL86RPpX/frz86bcbFC7KWmnC6C3R+bnpiM9WZuzKHl4OEPfSeboaTOQzSU1E4puIiIyE/ydLOaDe86BdRdq66xkZ5XWjfUtPt4EbuPF5FXWsGB7BIOZJfwyVmTgCP8PekT5U/vSD/6RAbQO9KPbsE+uJxrqMnTHy662ZT8VNj2numJKcyALe+YEhRjDntMmAqBXVrjyyAOQsFFREQazXrWQZNXJZg5KDabjZzi8rowsyezmN2ZRaTllpJVdJqsotMs3/v9vBlvdys9I/zoHelP70jTM9Mrwg+fsw+cDIqBMQ/BqFlweK2Z0LvrYxNoVjxlSrcRphem91Xg4dvKXwlpbQouIiLSLCwWC2H+noT5ezKm1/eTgEvKq9iXVcSu40XsySxid2YxezOLKKuoZmtGAVszCs56D+ga5F0XZkzxo2OgF5boERA9AiY+Z+bBbFsAad9A+mpTPn8A+lyNpd+vwFZjh6+AtAYFFxERaVG+Hq4kdg0isev382aqa2yk5ZayO7M2zNSGmuzi8rol2l/szKq738/Tld4RJsT0ivSnd+REel4/Ba+yY7BtkQkx+amwbQGu2xYwzj0EF98dMPBGCI61x6ctLUTBRUREWp3VxUL3MF+6h30/1ASQV1LOnsxi9pwJNJlFHMwuofh0FRvT89mYnl93r8UC0cE+9Iq8lF59rmaI20H6Zv8Xn4Of4l2eC2ueN6VLklla3XeymT8jTk3BRUREHEawrweXxHnUW6JdUVXDwewS9mYVsTfrTKgpJreknNTcUlJzS1my40zvzBUEuV/ONa4bmOq9kbiSTbhkrDerk774I/T+HzMfJnqUWckkTkfBRUREHJq7q4tZVh1Vv7ckp7icfVnF7M0yQWZvVhEHTpSQXwFvVQzjrbJhhJPPZOsafmldTVzVMdjxPux4n1Oe4ZT2uhb/oTfjHtHLTp+ZNIWCi4iIOKVQPw9C/er3zlRW13Agq5BFS1fjFdmdA9lhfJoZxWuFV9Lfksq11m+4yrqOwNMn8Ep5BVJeYY+1B9tCrqAk9iq6depIzwgzGficS7XF7hRcRESkzXCzuhAX5stFITYmjYvDzc2cfl14qpL9J4rZm3U1fzmei9/hrxhUsJRLSKF39X56n9hPedbfSK65iEerR7LFdSCxEYH0DPejZ4QfPcP96BHhR0hzHkApTaLgIiIibV6AlxsXdwvi4m5BQFcgEZvtD2Qdz6D0u4V0OPA+waUH+R/rBv7HuoFsWyAfHR/O4iMjec/Wue59gn3c6VEbZnqE+9Ej3Je4cD8CvNzs9rm1NwouIiLSLlksFiI7doWOs8D2R8jaDikLse34D2Fledzp+jl3un7OYfc4PmYU7xQPIq/Un/WpeaxPzav3XpEBnsSF+9EjzJcetaEmLsy3/mZ60iz0FRUREbFYIDIBIhOwjHsCDiabs5L2L6VrxQHu5QAzvN6hqPNl7AydxGouYm/OKfZnFXO89jTtzMLTfLM/p97bdgz0oke4rwkytT003cN88XbXj9+m0ldORETkbK7u5iTqXldAaR7s/ABS3sWSuY2Aw0sZfngpw72DIX4KXH49RYG9OZBdwv4TJew/UVxbSsgpLudYwSmOFZxixb76gaZTB6+6Xpm42l+7q4fmvOgrJCIi8mN8gmHInaac2G126N22CEqzYcOrsOFV/MP6kjjgBhL7T4HB3x/4mF9awYETxezPLjG/nijmwIkS8korOHryFEdPnqp3dhN830MTF+5H9zDfukDj56k5NGcouIiIiJyP8D4w/im47HE4tNyEmL1LIHsXLHsYkh+DuHFml96eEwny8WBITDBDYoLrvU1eSTkHs0vYn13CwdremQPZJeSW/HgPTYS/J3Hh5lDLuHBfuoeacBPk496KXwDHoOAiIiLSGFZX6DHelFMnYeeH5tTqo9/B/qWmeAZC/LWQcAN0vMjMoakV7OtBsK9Hg0BzsrSCA9klHMg2PTMHa39/oqi87nTt1Qdy6z0T5ONO91BfYmt7Zs6UqABPLJa2uQ+NgouIiEhTeXWAi6eZknvATOjdvgiKjsF3b5gS0hMGXA/9p4J/1I++VQcfdwZHBzE4Oqje9cJTlRzMLuFgdnFtmDGh5ujJU+SXVrCxtP4ZTgDe7lZiQ32JDfWpF2i6BPng7urSIl+K1qLgIiIi0hxC4mDsHLj0EUhbBSkLYc9nkLsPvnocvn4CYkabXpheV4C793m9bYCXG4ldO5DYtUO962UVVaTmlNaGmtqSU0J6billFdXsOFbIjmOF9Z6xuljoGuRNbJhvXbCJDfMlNsSXAG/nmEej4CIiItKcXKwQe6kpp4tg98cmxGSsM3NjDi0HD3/oew0MuBE6D6k3lHS+vN1d6dcxgH4dA+pdr6yu4XBeGQezSziUU8Kh2kBzKLuE0orquoMpkzlR77kQX3dizoSZUF9ian/t1MEbqwMdf6DgIiIi0lI8/eGim03JT4Vt75n5MAUZsOWfpgTFmAm9CddBYJeff8+f4WZ1qRsaOpvNZuNEUbkJMzmmhyY1p5RDOSVkFp4mt6SC3JJ8NqbVH3Zyt7rQLcSbmBATZmLOhBo79dIouIiIiLSGoBgY8xCMmgWH15oAs+tjE2hWPG1KtxGmF6b3leDh+7Nv2RgWi4WIAE8iAjwZ3j2k3p+VlFeRVhtiUnNKOFT7+7TcUsqramr3qClp8J7BPu7EhPoQHVIbaEJ8iAn1adG5NAouIiIircnFBaJHmDLxOTMPZtsCSPsG0leb8vnvoc/VMOAG6DrcPNOCfD1cie8UQHyn+sNONTU2jhWcqg00paTm1v6aU0pW0WnySivIK63gu/ST9T9FC3QO8ibKy9bsdVVwERERsRcPX7PiaMD1Zvho2yJIeRdOptVudrfADB+dGUoKimnV6rm4WOgc5E3nIG9G96z/Z6XlVaTVzpdJrQ02abW/L62o5nBeGWnlZc1eJwUXERERRxDYBUY9CCMfgCMbzNLqXR+ZQLPq/0zpMsyEnD7XmPkzduTjce7JwTabjZziclJzS9l+6Bh3vti8H1fBRURExJFYLNBlqCkT/w/2fm56YQ6tMCuTMtbBkj+YeTADrofoUWYlk4OwWCyE+XsS5u9JXADc2czvr+AiIiLiqNy8zA688ddC0XGzuV3KAsjdDzv+Y4p/R7O53YAbzF4ybZxzb58nIiLSXvhHwSX3w90b4Y7lMGgaeAaYXXrXvAB/GwRvjIXv3jRHEbRR6nERERFxJhYLdEo0ZcIzsP8Ls8Hdwa/MeUlHv4Ols6HXJLO0OmaMOV+pjWg7n4mIiEh74+YJfSebUnzCDB2lLIDs3WZi766PwDcc+k8xISast71rfME0VCQiItIW+IXDsHvgt+vgzm9gyHTwDoaSE7DuZZg3FP4+Cjb8HUrz7F3bJlOPi4iISFtisUBkginjnoSDyaYXZv9SyEwx5cuHoccEM6E3bjxYneOARVBwERERabtc3c1J1L2uML0sOz8wS6szt8He/5riHQzxU8zS6oj+TTrwsTUpuIiIiLQHPsEw5E5TTuw2AWb7f6A0Gza8akpYX9ML038K+IbZu8bnpDkuIiIi7U14H5jwNMzcAze8b3bitbpD9i5Y9jA83wsWTDWHQFaV27u29ajHRUREpL2yukKP8aacOgk7F5ul1cc2mTkx+5eCZ6DZAC/hBuh4kd2HkhRcREREBLw6wMV3mJKzv/aQx0VQfBy+e8OUkJ5mLkz/qWZDPDvQUJGIiIjUF9oDxj4O9++EX38I8b8CV0/I3QdfPQ5/6Qv/+gXs+AAqT7Vq1dTjIiIiIufmYoXul5lyutDMedm2EDLWw6GvTfHwNxvgDbgROg9u8aEkBRcRERH5eZ4BkHiLKXmHYNt7phRmwJZ3TAmKrR1Kug4CO7dINRRcREREpHGCY+HSh2H0bDi8xkzo3f0J5B+C5U/B8qchegSW6Kua/UMruIiIiEjTuLhA9EhTJv0J9nxqdulNXw1p3+C6d1Wzf0gFFxEREblwHr5m87oBN8DJw7B9EbZv/wXsbNYPo1VFIiIi0rw6dIVRf6Bq2opmf2sFFxEREWkZLbDCSMFFREREnIaCi4iIiDgNBRcRERFxGgouIiIi4jQUXERERMRpKLiIiIiI01BwEREREaeh4CIiIiJOQ8FFREREnIaCi4iIiDgNBRcRERFxGgouIiIi4jSaFFzmzZtHdHQ0np6eJCYmsnr16h+998MPP2TcuHGEhobi7+9PUlISX375ZZMrLCIiIu1Xo4PLokWLuO+++3j44YfZunUrI0aMYOLEiWRkZJzz/m+++YZx48axZMkSNm/ezJgxY7jyyivZunXrBVdeRERE2hfXxj7wwgsvMG3aNO644w4AXnzxRb788kteffVVnn322Qb3v/jii/VeP/PMM3zyySd89tlnDBw48Jwfo7y8nPLy8rrXRUVFAFRWVlJZWdnYKkszOvP1VzvYn9rCcagtHIvaw3G0RBs0KrhUVFSwefNmZs2aVe/6+PHjWbdu3Xm9R01NDcXFxQQFBf3oPc8++yxz585tcH3FihV4e3s3psrSQpKTk+1dBamltnAcagvHovawv7KysmZ/z0YFl9zcXKqrqwkPD693PTw8nKysrPN6j+eff57S0lKmTJnyo/fMnj2bmTNn1r0uKiqic+fOjBkzhuDg4MZUWZpZZWUlycnJjBs3Djc3N3tXp11TWzgOtYVjUXs4jry8vGZ/z0YPFQFYLJZ6r202W4Nr57Jw4UIef/xxPvnkE8LCwn70Pg8PDzw8PBpcd3Nz019CB6G2cBxqC8ehtnAsag/7a4mvf6OCS0hICFartUHvSnZ2doNemB9atGgR06ZN4/3332fs2LGNr6mIiIi0e41aVeTu7k5iYmKDccPk5GSGDRv2o88tXLiQW2+9lQULFnDFFVc0raYiIiLS7jV6qGjmzJncdNNNDBo0iKSkJF5//XUyMjKYPn06YOanHDt2jH/+85+ACS0333wzf/3rXxk6dGhdb42XlxcBAQHN+KmIiIhIW9fo4DJ16lTy8vJ44oknyMzMpF+/fixZsoSuXbsCkJmZWW9Pl7///e9UVVVx9913c/fdd9ddv+WWW5g/f/6FfwYiIiLSbjRpcu5dd93FXXfddc4/+2EYWblyZVM+hIiIiEgDOqtIREREnIaCi4iIiDgNBRcRERFxGgouIiIi4jQUXERERMRpKLiIiIiI01BwEREREaeh4CIiIiJOQ8FFREREnIaCi4iIiDgNBRcRERFxGgouIiIi4jQUXERERMRpKLiIiIiI01BwEREREaeh4CIiIiJOQ8FFREREnIaCi4iIiDgNBRcRERFxGgouIiIi4jQUXERERMRpKLiIiIiI01BwEREREaeh4CIiIiJOQ8FFREREnIaCi4iIiDgNBRcRERFxGgouIiIi4jQUXERERMRpKLiIiIiI01BwEREREaeh4CIiIiJOQ8FFREREnIaCi4iIiDgNBRcRERFxGgouIiIi4jQUXERERMRpKLiIiIiI01BwEREREaeh4CIiIiJOQ8FFREREnIaCi4iIiDgNBRcRERFxGgouIiIi4jQUXERERMRpKLiIiIiI01BwEREREaeh4CIiIiJOQ8FFREREnIaCi4iIiDgNBRcRERFxGgouIiIi4jQUXERERMRpKLiIiIiI01BwEREREaeh4CIiIiJOQ8FFREREnIaCi4iIiDgNBRcRERFxGgouIiIi4jQUXERERMRpKLiIiIiI01BwEREREafRpOAyb948oqOj8fT0JDExkdWrV//k/atWrSIxMRFPT09iYmJ47bXXmlRZERERad8aHVwWLVrEfffdx8MPP8zWrVsZMWIEEydOJCMj45z3p6WlMWnSJEaMGMHWrVt56KGHmDFjBosXL77gyouIiEj74trYB1544QWmTZvGHXfcAcCLL77Il19+yauvvsqzzz7b4P7XXnuNLl268OKLLwLQu3dvNm3axJ///Gd++ctfnvNjlJeXU15eXve6sLAQgPz8/MZWV5pZZWUlZWVl5OXl4ebmZu/qtGtqC8ehtnAsag/Hcebnts1ma743tTVCeXm5zWq12j788MN612fMmGEbOXLkOZ8ZMWKEbcaMGfWuffjhhzZXV1dbRUXFOZ+ZM2eODVBRUVFRUVFpA+XQoUONiRs/qVE9Lrm5uVRXVxMeHl7venh4OFlZWed8Jisr65z3V1VVkZubS2RkZINnZs+ezcyZM+teFxQU0LVrVzIyMggICGhMlaWZFRUV0blzZ44cOYK/v7+9q9OuqS0ch9rCsag9HEdhYSFdunQhKCio2d6z0UNFABaLpd5rm83W4NrP3X+u62d4eHjg4eHR4HpAQID+EjoIf39/tYWDUFs4DrWFY1F7OA4Xl+ZbxNyodwoJCcFqtTboXcnOzm7Qq3JGRETEOe93dXUlODi4kdUVERGR9qxRwcXd3Z3ExESSk5PrXU9OTmbYsGHnfCYpKanB/cuWLWPQoEGaNCUiIiKN0ui+m5kzZ/LGG2/w1ltvsWfPHu6//34yMjKYPn06YOan3HzzzXX3T58+ncOHDzNz5kz27NnDW2+9xZtvvskDDzxw3h/Tw8ODOXPmnHP4SFqX2sJxqC0ch9rCsag9HEdLtIXFZmv8GqV58+bx3HPPkZmZSb9+/fjLX/7CyJEjAbj11ltJT09n5cqVdfevWrWK+++/n127dhEVFcUf//jHuqAjIiIicr6aFFxERERE7EFnFYmIiIjTUHARERERp6HgIiIiIk5DwUVERESchsMEl3nz5hEdHY2npyeJiYmsXr36J+9ftWoViYmJeHp6EhMTw2uvvdZKNW37GtMWH374IePGjSM0NBR/f3+SkpL48ssvW7G2bVtjvy/OWLt2La6urgwYMKBlK9iONLYtysvLefjhh+natSseHh7Exsby1ltvtVJt27bGtsW7775LQkIC3t7eREZGctttt5GXl9dKtW27vvnmG6688kqioqKwWCx8/PHHP/tMs/zsbrZTjy7Ae++9Z3Nzc7P94x//sO3evdt277332nx8fGyHDx8+5/2pqak2b29v27333mvbvXu37R//+IfNzc3N9sEHH7RyzduexrbFvffea/u///s/28aNG2379++3zZ492+bm5mbbsmVLK9e87WlsW5xRUFBgi4mJsY0fP96WkJDQOpVt45rSFldddZVtyJAhtuTkZFtaWpptw4YNtrVr17ZirdumxrbF6tWrbS4uLra//vWvttTUVNvq1attffv2tV1zzTWtXPO2Z8mSJbaHH37YtnjxYhtg++ijj37y/ub62e0QwWXw4MG26dOn17vWq1cv26xZs855/x/+8Adbr1696l278847bUOHDm2xOrYXjW2Lc+nTp49t7ty5zV21dqepbTF16lTbI488YpszZ46CSzNpbFt88cUXtoCAAFteXl5rVK9daWxb/OlPf7LFxMTUu/bSSy/ZOnXq1GJ1bI/OJ7g0189uuw8VVVRUsHnzZsaPH1/v+vjx41m3bt05n1m/fn2D+ydMmMCmTZuorKxssbq2dU1pix+qqamhuLi4WU8CbY+a2hZvv/02hw4dYs6cOS1dxXajKW3x6aefMmjQIJ577jk6duxIjx49eOCBBzh16lRrVLnNakpbDBs2jKNHj7JkyRJsNhsnTpzggw8+4IorrmiNKstZmutnd5NOh25Oubm5VFdXNzikMTw8vMHhjGdkZWWd8/6qqipyc3OJjIxssfq2ZU1pix96/vnnKS0tZcqUKS1RxXajKW1x4MABZs2axerVq3F1tfu3dpvRlLZITU1lzZo1eHp68tFHH5Gbm8tdd91Ffn6+5rlcgKa0xbBhw3j33XeZOnUqp0+fpqqqiquuuoqXX365NaosZ2mun91273E5w2Kx1Htts9kaXPu5+891XRqvsW1xxsKFC3n88cdZtGgRYWFhLVW9duV826K6upobbriBuXPn0qNHj9aqXrvSmO+LmpoaLBYL7777LoMHD2bSpEm88MILzJ8/X70uzaAxbbF7925mzJjBY489xubNm1m6dClpaWk6dsZOmuNnt93/WxYSEoLVam2QlrOzsxskszMiIiLOeb+rqyvBwcEtVte2riltccaiRYuYNm0a77//PmPHjm3JarYLjW2L4uJiNm3axNatW/nd734HmB+eNpsNV1dXli1bxqWXXtoqdW9rmvJ9ERkZSceOHQkICKi71rt3b2w2G0ePHiUuLq5F69xWNaUtnn32WYYPH86DDz4IQP/+/fHx8WHEiBE89dRT6qFvRc31s9vuPS7u7u4kJiaSnJxc73pycjLDhg075zNJSUkN7l+2bBmDBg3Czc2txera1jWlLcD0tNx6660sWLBA48bNpLFt4e/vz44dO0hJSakr06dPp2fPnqSkpDBkyJDWqnqb05Tvi+HDh3P8+HFKSkrqru3fvx8XFxc6derUovVty5rSFmVlZbi41P9RZ7Vage//ty+to9l+djdqKm8LObO87c0337Tt3r3bdt9999l8fHxs6enpNpvNZps1a5btpptuqrv/zJKq+++/37Z7927bm2++qeXQzaSxbbFgwQKbq6ur7ZVXXrFlZmbWlYKCAnt9Cm1GY9vih7SqqPk0ti2Ki4ttnTp1sl177bW2Xbt22VatWmWLi4uz3XHHHfb6FNqMxrbF22+/bXN1dbXNmzfPdujQIduaNWtsgwYNsg0ePNhen0KbUVxcbNu6datt69atNsD2wgsv2LZu3Vq3NL2lfnY7RHCx2Wy2V155xda1a1ebu7u77aKLLrKtWrWq7s9uueUW26hRo+rdv3LlStvAgQNt7u7utm7dutleffXVVq5x29WYthg1apQNaFBuueWW1q94G9TY74uzKbg0r8a2xZ49e2xjx461eXl52Tp16mSbOXOmraysrJVr3TY1ti1eeuklW58+fWxeXl62yMhI24033mg7evRoK9e67VmxYsVP/vvfUj+7LTab+spERETEOdh9jouIiIjI+VJwEREREaeh4CIiIiJOQ8FFREREnIaCi4iIiDgNBRcRERFxGgouIiIi4jQUXERERMRpKLiIiIiI01BwEREREaeh4CIiIiJO4/8DXIZ5Z+12A7kAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Plot m0_dgp and m1_dgp over the unit interval\n",
        "import matplotlib.pyplot as plt\n",
        "u = np.linspace(0, 1, 100)\n",
        "plt.plot(u, m0_dgp(u), label='m0_dgp')\n",
        "plt.plot(u, m1_dgp(u), label='m1_dgp')\n",
        "plt.grid()\n",
        "plt.legend()\n",
        "plt.xlim(0, 1)\n",
        "plt.ylim(0, 1)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 239,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "0.04591666666666566"
            ]
          },
          "execution_count": 239,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Compute implied LATE\n",
        "a = integrate.quad(lambda u: m1_dgp(u) * w_target_late(1, u, 0.35, 0.9), 0, 1)[0]\n",
        "b = integrate.quad(lambda u: m0_dgp(u) * w_target_late(0, u, 0.35, 0.9), 0, 1)[0]\n",
        "a + b"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Diet Model\n",
        "\n",
        "Problem: minimize linear mapping over parameters $(\\theta_0, \\theta_1)$ subject to the constraint that parameters also give identified parameters.\n",
        "\n",
        "Sets:\n",
        "- `THETA`: Parameter space.\n",
        "- `THETA0`: parameter space of the model for $m_0$\n",
        "- `THETA1`: parameter space of the model for $m_1$\n",
        "- `IDENTIF`: set of identified IV-like estimands.\n",
        "\n",
        "Parameters:\n",
        "- `iv {IDENTIF}`: value of identified IV parameter(s).\n",
        "- `gamma {THETA}`: $\\Gamma^*_{dk}(b_{dk})$ for $d\\in\\{0,1\\}$ and $k=1,\\cdots,K_d$ - contribution to objective by each parameter.\n",
        "\n",
        "Variables:\n",
        "- `t0 {THETA0}`: value of parameters for $m_0$\n",
        "- `t1 {THETA1}`: value of parameters for $m_1$\n",
        "\n",
        "Objective:\n",
        "- `beta_star`: linear map giving $\\overline{\\beta}_{fd}^*$: sum over terms in `Gamma`.\n",
        "\n",
        "Constraints:\n",
        "- `iv {IV_identified}`: ensure that parameters result in identified IV-like estimands.\n",
        "\n",
        "The problem is then modeled as follows:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 240,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Compute the contribution for K_0 = 3 and K_1 = 3\n",
        "gamma_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"theta_01\", gamma_star0_late(bern_2_0, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_02\", gamma_star0_late(bern_2_1, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_03\", gamma_star0_late(bern_2_2, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_11\", gamma_star1_late(bern_2_0, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_12\", gamma_star1_late(bern_2_1, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_13\", gamma_star1_late(bern_2_2, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "    ],\n",
        "    columns=[\"THETA\", \"gamma\"],\n",
        ").set_index(\"THETA\")\n",
        "\n",
        "# Compute the contribution for K_0 = 3 and K_1 = 3\n",
        "gamma_ident_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"theta_01\", gamma_star0_late(bern_2_0, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "        (\"theta_02\", gamma_star0_late(bern_2_1, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "        (\"theta_03\", gamma_star0_late(bern_2_2, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "        (\"theta_11\", gamma_star1_late(bern_2_0, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "        (\"theta_12\", gamma_star1_late(bern_2_1, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "        (\"theta_13\", gamma_star1_late(bern_2_2, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "    ],\n",
        "    columns=[\"THETA\", \"gamma_ident\"],\n",
        ").set_index(\"THETA\")\n",
        "\n",
        "# Join gamma_df and gamma_ident_df\n",
        "gamma_df = gamma_df.join(gamma_ident_df)\n",
        "\n",
        "a = integrate.quad(lambda u: m1_dgp(u) \n",
        "    * w_target_late(1, u, 0.35, 0.6), 0, 1)[0]\n",
        "\n",
        "b = integrate.quad(lambda u: m0_dgp(u) \n",
        "    * w_target_late(0, u, 0.35, 0.6), 0, 1)[0]\n",
        "\n",
        "iv_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"iv_1\", a+b)\n",
        "    ],\n",
        "    columns = [\"IDENTIF\", \"iv\"]\n",
        ").set_index(\"IDENTIF\")\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 241,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>iv</th>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>IDENTIF</th>\n",
              "      <th></th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>iv_1</th>\n",
              "      <td>0.079417</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "               iv\n",
              "IDENTIF          \n",
              "iv_1     0.079417"
            ]
          },
          "execution_count": 241,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "iv_df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 242,
      "metadata": {},
      "outputs": [],
      "source": [
        "%%ampl_eval\n",
        "reset;\n",
        "set THETA;\n",
        "set IDENTIF;\n",
        "\n",
        "param gamma {THETA};\n",
        "param gamma_ident {THETA};\n",
        "param iv {IDENTIF};\n",
        "\n",
        "var Theta_val {j in THETA} >= 0, <= 1;\n",
        "\n",
        "maximize beta: sum {j in THETA} gamma[j] * Theta_val[j];\n",
        "\n",
        "subject to Identified {i in IDENTIF}:\n",
        "sum {j in THETA} gamma_ident[j] * Theta_val[j] == iv[i];\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 243,
      "metadata": {},
      "outputs": [],
      "source": [
        "ampl.set_data(gamma_df, \"THETA\")\n",
        "ampl.set_data(iv_df, \"IDENTIF\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 244,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "HiGHS 1.5.1: HiGHS 1.5.1: optimal solution; objective 0.3608654721\n",
            "0 simplex iterations\n",
            "0 barrier iterations\n",
            "Objective is: 0.3608654721269454\n",
            "argmax:           Theta_val.val\n",
            "theta_01       1.000000\n",
            "theta_02       0.000000\n",
            "theta_03       0.000000\n",
            "theta_11       0.000000\n",
            "theta_12       0.265017\n",
            "theta_13       1.000000\n"
          ]
        }
      ],
      "source": [
        "# Specify the solver to use (e.g., HiGHS)\n",
        "ampl.option[\"solver\"] = \"highs\"\n",
        "# Solve\n",
        "ampl.solve()\n",
        "\n",
        "# Stop if the model was not solved\n",
        "assert ampl.get_value(\"solve_result\") == \"solved\"\n",
        "# Get objective entity by AMPL name\n",
        "beta_hi = ampl.get_objective('beta')\n",
        "# Print it\n",
        "print(\"Objective is:\", beta_hi.value())\n",
        "argmax = ampl.get_variable('Theta_val').get_values().to_pandas()\n",
        "print(\"argmax:\", argmax)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 245,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Code up constant splines\n",
        "def cs_basis(u_lo, u_hi, u):\n",
        "    \"\"\"Constant spline basis function (no covariates)\"\"\"\n",
        "    if u_lo <= u < u_hi: return 1\n",
        "    else: return 0\n",
        "\n",
        "cs_0_35 = lambda u: cs_basis(0, 0.35, u)\n",
        "cs_35_60 = lambda u: cs_basis(0.35, 0.60, u)\n",
        "cs_60_90 = lambda u: cs_basis(0.60, 0.90, u)\n",
        "cs_90_1 = lambda u: cs_basis(0.90, 1, u)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 246,
      "metadata": {},
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "C:\\Users\\budde\\AppData\\Local\\Temp\\ipykernel_8704\\3842910436.py:10: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
            "  If increasing the limit yields no improvement it is advised to analyze \n",
            "  the integrand in order to determine the difficulties.  If the position of a \n",
            "  local difficulty can be determined (singularity, discontinuity) one will \n",
            "  probably gain from splitting up the interval and calling the integrator \n",
            "  on the subranges.  Perhaps a special-purpose integrator should be used.\n",
            "  return integrate.quad(func, 0, 1)[0]\n",
            "C:\\Users\\budde\\AppData\\Local\\Temp\\ipykernel_8704\\3842910436.py:19: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
            "  If increasing the limit yields no improvement it is advised to analyze \n",
            "  the integrand in order to determine the difficulties.  If the position of a \n",
            "  local difficulty can be determined (singularity, discontinuity) one will \n",
            "  probably gain from splitting up the interval and calling the integrator \n",
            "  on the subranges.  Perhaps a special-purpose integrator should be used.\n",
            "  return integrate.quad(func, 0, 1)[0]\n",
            "C:\\Users\\budde\\AppData\\Local\\Temp\\ipykernel_8704\\3842910436.py:10: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
            "  If increasing the limit yields no improvement it is advised to analyze \n",
            "  the integrand in order to determine the difficulties.  If the position of a \n",
            "  local difficulty can be determined (singularity, discontinuity) one will \n",
            "  probably gain from splitting up the interval and calling the integrator \n",
            "  on the subranges.  Perhaps a special-purpose integrator should be used.\n",
            "  return integrate.quad(func, 0, 1)[0]\n",
            "C:\\Users\\budde\\AppData\\Local\\Temp\\ipykernel_8704\\3842910436.py:19: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
            "  If increasing the limit yields no improvement it is advised to analyze \n",
            "  the integrand in order to determine the difficulties.  If the position of a \n",
            "  local difficulty can be determined (singularity, discontinuity) one will \n",
            "  probably gain from splitting up the interval and calling the integrator \n",
            "  on the subranges.  Perhaps a special-purpose integrator should be used.\n",
            "  return integrate.quad(func, 0, 1)[0]\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "(              gamma  gamma_ident\n",
              " THETA                           \n",
              " theta_0_1  0.000000          0.0\n",
              " theta_0_2 -0.454545         -1.0\n",
              " theta_0_3 -0.545455          0.0\n",
              " theta_0_4  0.000000          0.0\n",
              " theta_1_1  0.000000          0.0\n",
              " theta_1_2  0.454545          1.0\n",
              " theta_1_3  0.545455          0.0\n",
              " theta_1_4  0.000000          0.0,\n",
              "                iv\n",
              " IDENTIF          \n",
              " iv_1     0.079417)"
            ]
          },
          "execution_count": 246,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Repeat for constant spline basis\n",
        "gamma_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"theta_0_1\", gamma_star0_late(cs_0_35, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_0_2\", gamma_star0_late(cs_35_60, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_0_3\", gamma_star0_late(cs_60_90, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_0_4\", gamma_star0_late(cs_90_1, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_1_1\", gamma_star1_late(cs_0_35, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_1_2\", gamma_star1_late(cs_35_60, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_1_3\", gamma_star1_late(cs_60_90, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_1_4\", gamma_star1_late(cs_90_1, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "    ],\n",
        "    columns=[\"THETA\", \"gamma\"],\n",
        ").set_index(\"THETA\")\n",
        "\n",
        "# Compute the contribution for K_0 = 3 and K_1 = 3\n",
        "gamma_ident_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"theta_0_1\", gamma_star0_late(cs_0_35, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "        (\"theta_0_2\", gamma_star0_late(cs_35_60, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "        (\"theta_0_3\", gamma_star0_late(cs_60_90, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "        (\"theta_0_4\", gamma_star0_late(cs_90_1, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "\n",
        "        (\"theta_1_1\", gamma_star1_late(cs_0_35, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "        (\"theta_1_2\", gamma_star1_late(cs_35_60, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "        (\"theta_1_3\", gamma_star1_late(cs_60_90, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "        (\"theta_1_4\", gamma_star1_late(cs_90_1, w_target_late, u_lo=0.35, u_hi=0.6)),\n",
        "    ],\n",
        "    columns=[\"THETA\", \"gamma_ident\"],\n",
        ").set_index(\"THETA\")\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "# Join gamma_df and gamma_ident_df\n",
        "gamma_df = gamma_df.join(gamma_ident_df)\n",
        "\n",
        "a = integrate.quad(lambda u: m1_dgp(u) \n",
        "    * w_target_late(1, u, 0.35, 0.6), 0, 1)[0]\n",
        "\n",
        "b = integrate.quad(lambda u: m0_dgp(u) \n",
        "    * w_target_late(0, u, 0.35, 0.6), 0, 1)[0]\n",
        "\n",
        "iv_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"iv_1\", a+b)\n",
        "    ],\n",
        "    columns = [\"IDENTIF\", \"iv\"]\n",
        ").set_index(\"IDENTIF\")\n",
        "\n",
        "gamma_df, iv_df\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 247,
      "metadata": {},
      "outputs": [],
      "source": [
        "\n",
        "%%ampl_eval\n",
        "reset;\n",
        "set THETA;\n",
        "set IDENTIF;\n",
        "\n",
        "param gamma {THETA};\n",
        "param gamma_ident {THETA};\n",
        "param iv {IDENTIF};\n",
        "\n",
        "var Theta_val {j in THETA} >= 0, <= 1;\n",
        "\n",
        "maximize beta: sum {j in THETA} gamma[j] * Theta_val[j];\n",
        "\n",
        "subject to Identified {i in IDENTIF}:\n",
        "sum {j in THETA} gamma_ident[j] * Theta_val[j] == iv[i];\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 248,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "HiGHS 1.5.1: HiGHS 1.5.1: optimal solution; objective 0.5815530303\n",
            "0 simplex iterations\n",
            "0 barrier iterations\n",
            "Objective is: 0.5815530303030305\n",
            "argmax:            Theta_val.val\n",
            "theta_0_1       0.000000\n",
            "theta_0_2       0.000000\n",
            "theta_0_3       0.000000\n",
            "theta_0_4       0.000000\n",
            "theta_1_1       0.000000\n",
            "theta_1_2       0.079417\n",
            "theta_1_3       1.000000\n",
            "theta_1_4       0.000000\n"
          ]
        }
      ],
      "source": [
        "ampl.set_data(gamma_df, \"THETA\")\n",
        "ampl.set_data(iv_df, \"IDENTIF\")\n",
        "\n",
        "# Specify the solver to use (e.g., HiGHS)\n",
        "ampl.option[\"solver\"] = \"highs\"\n",
        "# Solve\n",
        "ampl.solve()\n",
        "\n",
        "# Stop if the model was not solved\n",
        "assert ampl.get_value(\"solve_result\") == \"solved\"\n",
        "# Get objective entity by AMPL name\n",
        "beta_hi = ampl.get_objective('beta')\n",
        "# Print it\n",
        "print(\"Objective is:\", beta_hi.value())\n",
        "argmax = ampl.get_variable('Theta_val').get_values().to_pandas()\n",
        "print(\"argmax:\", argmax)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 249,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Weights for IV-like estimands\n",
        "def s_iv_slope(z, exp_z, cov_d_z):\n",
        "    \"\"\"IV-like specification s(d,z): IV slope\n",
        "    Args:\n",
        "        z (np.int): value of the instrument\n",
        "        exp_z (np.float): expected value of the instrument\n",
        "        cov_d_z (np.float): covariance between treatment and instrument\n",
        "    \"\"\"\n",
        "    return (z - exp_z) / cov_d_z"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 250,
      "metadata": {},
      "outputs": [],
      "source": [
        "def compute_iv_slope(m1, m0):\n",
        "    \"\"\"IV-like specification gamma_0(d,z): IV slope\n",
        "    Args:\n",
        "        d (np.int): value of treatment\n",
        "        z (np.int): value of the instrument\n",
        "        exp_z (np.float): expected value of the instrument\n",
        "        cov_d_z (np.float): covariance between treatment and instrument\n",
        "\n",
        "        z takes 3 values: 0, 1, 2\n",
        "        P(z=k): 0.5, 0.4, 0.1\n",
        "        P(d=1|z=k): 0.35, 0.6, 0.7\n",
        "    \"\"\"\n",
        "\n",
        "    exp_z = 0 * 0.5 + 1 * 0.4 + 2 * 0.1\n",
        "    ed = 0.35 * 0.5 + 0.6 * 0.4 + 0.7 * 0.1\n",
        "    ez = 0.5 * 0 + 0.4 * 1 + 0.1 * 2\n",
        "    edz = 0.6 * 0.4 + 2 * 0.7 * 0.1    \n",
        "    \n",
        "    cov_d_z = edz - ed * ez\n",
        "\n",
        "    def check(u): return 1 if u >= 0.35 else 0\n",
        "    def weight(u): return s_iv_slope(1, 0, exp_z, cov_d_z) * check(u)\n",
        "    a0 = integrate.quad(lambda u: m1(u) * weight(u), 0, 1)[0]\n",
        "\n",
        "    def check(u): return 1 if u >= 0.6 else 0\n",
        "    def weight(u): return s_iv_slope(1, 1, exp_z, cov_d_z) * check(u)\n",
        "    a1 = integrate.quad(lambda u: m1(u) * weight(u), 0, 1)[0]\n",
        "\n",
        "    def check(u): return 1 if u >= 0.7 else 0\n",
        "    def weight(u): return s_iv_slope(1, 2, exp_z, cov_d_z) * check(u)\n",
        "    a2 = integrate.quad(lambda u: m1(u) * weight(u), 0, 1)[0]\n",
        "\n",
        "    # Repeat for d = 0\n",
        "    def check(u): return 1 if u < 0.35 else 0\n",
        "    def weight(u): return s_iv_slope(0, 0, exp_z, cov_d_z) * check(u)\n",
        "    b0 = integrate.quad(lambda u: m0(u) * weight(u), 0, 1)[0]\n",
        "\n",
        "    def check(u): return 1 if u < 0.6 else 0\n",
        "    def weight(u): return s_iv_slope(0, 1, exp_z, cov_d_z) * check(u)\n",
        "    b1 = integrate.quad(lambda u: m0(u) * weight(u), 0, 1)[0]\n",
        "\n",
        "    def check(u): return 1 if u < 0.7 else 0\n",
        "    def weight(u): return s_iv_slope(0, 2, exp_z, cov_d_z) * check(u)\n",
        "    b2 = integrate.quad(lambda u: m0(u) * weight(u), 0, 1)[0]\n",
        "\n",
        "    return a0, a1, a2, b0, b1, b2\n",
        "    # return 0.5 * (a0 + b0) + 0.4 * (a1 + b1) + 0.1 * (a2 + b2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 251,
      "metadata": {},
      "outputs": [],
      "source": [
        "def gamma_0_iv_slope(m0, weight, supp_z, f_z, prop_z):\n",
        "\n",
        "    \"\"\"\n",
        "    Args:\n",
        "        m0 (function): MTE function for d = 0\n",
        "        weight (function): weights corresponding to IV slope parameter\n",
        "        supp_z (np.array): support of z\n",
        "        f_z (np.array): density of z\n",
        "        prop_z (np.array): propensity score given z\n",
        "    \n",
        "    Returns:\n",
        "        val (float): value of gamma_0 at m0\n",
        "    \"\"\"\n",
        "    # Compute expected value of z\n",
        "    ez = np.sum(supp_z * f_z)\n",
        "    ed = np.sum(prop_z * f_z)\n",
        "    edz = np.sum(supp_z * prop_z * f_z)\n",
        "\n",
        "    # Compute covariance between d and z\n",
        "    cov_dz = edz - ed * ez\n",
        "\n",
        "    # Define function to be integrated over\n",
        "    def func(u, z): \n",
        "        if prop_z[np.where(supp_z == z)[0][0]] < u: return m0(u) * weight(z, ez, cov_dz)\n",
        "        else : return 0\n",
        "\n",
        "    # Integrate func over u in [0,1] for every z in supp_z\n",
        "    return np.sum([integrate.quad(func, 0, 1, args=(z,))[0] * f_z[i] \n",
        "         for i, z in enumerate(supp_z)])\n",
        "\n",
        "def gamma_1_iv_slope(m1, weight, supp_z, f_z, prop_z):\n",
        "\n",
        "    \"\"\"\n",
        "    Args:\n",
        "        m1 (function): MTE function for d = 1\n",
        "        weight (function): weights corresponding to IV slope parameter\n",
        "        supp_z (np.array): support of z\n",
        "        f_z (np.array): density of z\n",
        "        prop_z (np.array): propensity score given z\n",
        "    \n",
        "    Returns:\n",
        "        val (float): value of gamma_0 at m1\n",
        "    \"\"\"\n",
        "    # Compute expected value of z\n",
        "    ez = np.sum(supp_z * f_z)\n",
        "    ed = np.sum(prop_z * f_z)\n",
        "    edz = np.sum(supp_z * prop_z * f_z)\n",
        "\n",
        "    # Compute covariance between d and z\n",
        "    cov_dz = edz - ed * ez\n",
        "\n",
        "    # Define function to be integrated over\n",
        "    def func(u, z): \n",
        "        if prop_z[np.where(supp_z == z)[0][0]] > u: return m1(u) * weight(z, ez, cov_dz)\n",
        "        else : return 0\n",
        "\n",
        "    # Integrate func over u in [0,1] for every z in supp_z\n",
        "    return np.sum([integrate.quad(func, 0, 1, args=(z,))[0] * f_z[i] \n",
        "         for i, z in enumerate(supp_z)])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 252,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "(-0.4251001872659388, 0.49873595505617896, 0.07363576779024017)"
            ]
          },
          "execution_count": 252,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# get index of elements in supp_z which are equal to 0\n",
        "supp_z = np.array([0,1,2])\n",
        "f_z = np.array([0.5, 0.4, 0.1])\n",
        "prop_z = np.array([0.35, 0.6, 0.7])\n",
        "\n",
        "a, b = gamma_0_iv_slope(m0_dgp, s_iv_slope, supp_z, f_z, prop_z), gamma_1_iv_slope(m1_dgp, s_iv_slope, supp_z, f_z, prop_z)\n",
        "a, b, a+b"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 273,
      "metadata": {},
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "C:\\Users\\budde\\AppData\\Local\\Temp\\ipykernel_8704\\3842910436.py:10: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
            "  If increasing the limit yields no improvement it is advised to analyze \n",
            "  the integrand in order to determine the difficulties.  If the position of a \n",
            "  local difficulty can be determined (singularity, discontinuity) one will \n",
            "  probably gain from splitting up the interval and calling the integrator \n",
            "  on the subranges.  Perhaps a special-purpose integrator should be used.\n",
            "  return integrate.quad(func, 0, 1)[0]\n",
            "C:\\Users\\budde\\AppData\\Local\\Temp\\ipykernel_8704\\3842910436.py:19: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
            "  If increasing the limit yields no improvement it is advised to analyze \n",
            "  the integrand in order to determine the difficulties.  If the position of a \n",
            "  local difficulty can be determined (singularity, discontinuity) one will \n",
            "  probably gain from splitting up the interval and calling the integrator \n",
            "  on the subranges.  Perhaps a special-purpose integrator should be used.\n",
            "  return integrate.quad(func, 0, 1)[0]\n",
            "C:\\Users\\budde\\AppData\\Local\\Temp\\ipykernel_8704\\1355938745.py:28: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
            "  If increasing the limit yields no improvement it is advised to analyze \n",
            "  the integrand in order to determine the difficulties.  If the position of a \n",
            "  local difficulty can be determined (singularity, discontinuity) one will \n",
            "  probably gain from splitting up the interval and calling the integrator \n",
            "  on the subranges.  Perhaps a special-purpose integrator should be used.\n",
            "  return np.sum([integrate.quad(func, 0, 1, args=(z,))[0] * f_z[i]\n",
            "C:\\Users\\budde\\AppData\\Local\\Temp\\ipykernel_8704\\1355938745.py:58: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
            "  If increasing the limit yields no improvement it is advised to analyze \n",
            "  the integrand in order to determine the difficulties.  If the position of a \n",
            "  local difficulty can be determined (singularity, discontinuity) one will \n",
            "  probably gain from splitting up the interval and calling the integrator \n",
            "  on the subranges.  Perhaps a special-purpose integrator should be used.\n",
            "  return np.sum([integrate.quad(func, 0, 1, args=(z,))[0] * f_z[i]\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "(              gamma   gamma_ident\n",
              " THETA                            \n",
              " theta_0_1  0.000000  0.000000e+00\n",
              " theta_0_2 -0.454545 -8.426966e-01\n",
              " theta_0_3 -0.545455 -1.573034e-01\n",
              " theta_0_4  0.000000  2.775558e-17\n",
              " theta_1_1  0.000000 -3.330669e-16\n",
              " theta_1_2  0.454545  8.426966e-01\n",
              " theta_1_3  0.545455  1.573034e-01\n",
              " theta_1_4  0.000000  0.000000e+00,\n",
              "                iv\n",
              " IDENTIF          \n",
              " iv_1     0.073636)"
            ]
          },
          "execution_count": 273,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Repeat for IV slope coefficient (Figure 2 in the paper)\n",
        "supp_z = np.array([0,1,2])\n",
        "f_z = np.array([0.5, 0.4, 0.1])\n",
        "prop_z = np.array([0.35, 0.6, 0.7])\n",
        "\n",
        "# Cotribution to target parameter\n",
        "gamma_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"theta_0_1\", gamma_star0_late(cs_0_35, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_0_2\", gamma_star0_late(cs_35_60, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_0_3\", gamma_star0_late(cs_60_90, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_0_4\", gamma_star0_late(cs_90_1, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_1_1\", gamma_star1_late(cs_0_35, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_1_2\", gamma_star1_late(cs_35_60, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_1_3\", gamma_star1_late(cs_60_90, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "        (\"theta_1_4\", gamma_star1_late(cs_90_1, w_target_late, u_lo=0.35, u_hi=0.9)),\n",
        "    ],\n",
        "    columns=[\"THETA\", \"gamma\"],\n",
        ").set_index(\"THETA\")\n",
        "\n",
        "# Contribution to identified parameter\n",
        "gamma_ident_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"theta_0_1\", gamma_0_iv_slope(cs_0_35, s_iv_slope, supp_z, f_z, prop_z)),\n",
        "        (\"theta_0_2\", gamma_0_iv_slope(cs_35_60, s_iv_slope, supp_z, f_z, prop_z)),\n",
        "        (\"theta_0_3\", gamma_0_iv_slope(cs_60_90, s_iv_slope, supp_z, f_z, prop_z)),\n",
        "        (\"theta_0_4\", gamma_0_iv_slope(cs_90_1, s_iv_slope, supp_z, f_z, prop_z)),\n",
        "\n",
        "        (\"theta_1_1\", gamma_1_iv_slope(cs_0_35, s_iv_slope, supp_z, f_z, prop_z)),\n",
        "        (\"theta_1_2\", gamma_1_iv_slope(cs_35_60, s_iv_slope, supp_z, f_z, prop_z)),\n",
        "        (\"theta_1_3\", gamma_1_iv_slope(cs_60_90, s_iv_slope, supp_z, f_z, prop_z)),\n",
        "        (\"theta_1_4\", gamma_1_iv_slope(cs_90_1, s_iv_slope, supp_z, f_z, prop_z)),\n",
        "    ],\n",
        "    columns=[\"THETA\", \"gamma_ident\"],\n",
        ").set_index(\"THETA\")\n",
        "\n",
        "# Join gamma_df and gamma_ident_df\n",
        "gamma_df = gamma_df.join(gamma_ident_df)\n",
        "\n",
        "a, b = gamma_0_iv_slope(m0_dgp, s_iv_slope, supp_z, f_z, prop_z), gamma_1_iv_slope(m1_dgp, s_iv_slope, supp_z, f_z, prop_z)\n",
        "\n",
        "iv_df = pd.DataFrame(\n",
        "    [\n",
        "        (\"iv_1\", a+b)\n",
        "    ],\n",
        "    columns = [\"IDENTIF\", \"iv\"]\n",
        ").set_index(\"IDENTIF\")\n",
        "\n",
        "gamma_df, iv_df\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 254,
      "metadata": {},
      "outputs": [],
      "source": [
        "%%ampl_eval\n",
        "reset;\n",
        "set THETA;\n",
        "set IDENTIF;\n",
        "\n",
        "param gamma {THETA};\n",
        "param gamma_ident {THETA};\n",
        "param iv {IDENTIF};\n",
        "\n",
        "var Theta_val {j in THETA} >= 0, <= 1;\n",
        "\n",
        "maximize beta: sum {j in THETA} gamma[j] * Theta_val[j];\n",
        "subject to Identified {i in IDENTIF}:\n",
        "sum {j in THETA} gamma_ident[j] * Theta_val[j] == iv[i];"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 255,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "HiGHS 1.5.1:   Warning code 1 for call Highs_addRows(lp(), AccConstraints.lb.size(), AccConstraints.lb.data(), AccConstraints.ub.data(), AccConstraints.coeffs.size(), AccConstraints.starts.data(), AccConstraints.indices.data(), AccConstraints.coeffs.data())\n",
            "HiGHS 1.5.1: optimal solution; objective 0.5003247475\n",
            "0 simplex iterations\n",
            "0 barrier iterations\n",
            "Objective is: 0.5003247474747358\n",
            "argmax:            Theta_val.val\n",
            "theta_0_1       0.000000\n",
            "theta_0_2       0.099286\n",
            "theta_0_3       0.000000\n",
            "theta_0_4       0.000000\n",
            "theta_1_1       0.000000\n",
            "theta_1_2       0.000000\n",
            "theta_1_3       1.000000\n",
            "theta_1_4       0.000000\n"
          ]
        }
      ],
      "source": [
        "\n",
        "ampl.set_data(gamma_df, \"THETA\")\n",
        "ampl.set_data(iv_df, \"IDENTIF\")\n",
        "\n",
        "# Specify the solver to use (e.g., HiGHS)\n",
        "ampl.option[\"solver\"] = \"highs\"\n",
        "# Solve\n",
        "ampl.solve()\n",
        "\n",
        "# Stop if the model was not solved\n",
        "assert ampl.get_value(\"solve_result\") == \"solved\"\n",
        "# Get objective entity by AMPL name\n",
        "beta_hi = ampl.get_objective('beta')\n",
        "# Print it\n",
        "print(\"Objective is:\", beta_hi.value())\n",
        "argmax = ampl.get_variable('Theta_val').get_values().to_pandas()\n",
        "print(\"argmax:\", argmax)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 256,
      "metadata": {},
      "outputs": [],
      "source": [
        "%%ampl_eval\n",
        "reset;\n",
        "set THETA;\n",
        "set IDENTIF;\n",
        "\n",
        "param gamma {THETA};\n",
        "param gamma_ident {THETA};\n",
        "param iv {IDENTIF};\n",
        "\n",
        "var Theta_val {j in THETA} >= 0, <= 1;\n",
        "\n",
        "minimize beta: sum {j in THETA} gamma[j] * Theta_val[j];\n",
        "subject to Identified {i in IDENTIF}:\n",
        "sum {j in THETA} gamma_ident[j] * Theta_val[j] == iv[i];"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 257,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "HiGHS 1.5.1:   Warning code 1 for call Highs_addRows(lp(), AccConstraints.lb.size(), AccConstraints.lb.data(), AccConstraints.ub.data(), AccConstraints.coeffs.size(), AccConstraints.starts.data(), AccConstraints.indices.data(), AccConstraints.coeffs.data())\n",
            "HiGHS 1.5.1: optimal solution; objective -0.4208873737\n",
            "0 simplex iterations\n",
            "0 barrier iterations\n",
            "Objective is: -0.4208873737373858\n",
            "argmax:            Theta_val.val\n",
            "theta_0_1       0.000000\n",
            "theta_0_2       0.000000\n",
            "theta_0_3       1.000000\n",
            "theta_0_4       0.000000\n",
            "theta_1_1       0.000000\n",
            "theta_1_2       0.274048\n",
            "theta_1_3       0.000000\n",
            "theta_1_4       0.000000\n"
          ]
        }
      ],
      "source": [
        "\n",
        "ampl.set_data(gamma_df, \"THETA\")\n",
        "ampl.set_data(iv_df, \"IDENTIF\")\n",
        "\n",
        "# Specify the solver to use (e.g., HiGHS)\n",
        "ampl.option[\"solver\"] = \"highs\"\n",
        "# Solve\n",
        "ampl.solve()\n",
        "\n",
        "# Stop if the model was not solved\n",
        "assert ampl.get_value(\"solve_result\") == \"solved\"\n",
        "# Get objective entity by AMPL name\n",
        "beta_hi = ampl.get_objective('beta')\n",
        "# Print it\n",
        "print(\"Objective is:\", beta_hi.value())\n",
        "argmax = ampl.get_variable('Theta_val').get_values().to_pandas()\n",
        "print(\"argmax:\", argmax)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 296,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "(array([0., 0., 0., 0., 1., 1., 1., 1.]),\n",
              " ['theta0_0',\n",
              "  'theta0_1',\n",
              "  'theta0_2',\n",
              "  'theta0_3',\n",
              "  'theta1_0',\n",
              "  'theta1_1',\n",
              "  'theta1_2',\n",
              "  'theta1_3'])"
            ]
          },
          "execution_count": 296,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "d0 = [\"theta0_\" + str(i) for i in range(4)]\n",
        "d1 = [\"theta1_\" + str(i) for i in range(4)]\n",
        "d = d0 + d1\n",
        "\n",
        "a = np.zeros(4)\n",
        "b = np.ones(4)\n",
        "\n",
        "d, a, b\n",
        "\n",
        "c = np.concatenate((a, b))\n",
        "c, d\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 304,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "'THETA'"
            ]
          },
          "execution_count": 304,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "df.index.rename(\"THETA\", inplace=True)\n",
        "df.index.name"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "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.10.9"
    },
    "vscode": {
      "interpreter": {
        "hash": "24209f11fee719204716071d008026089e3b3353af3d540e4dfc4f270a7a7039"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}