diff --git a/docs/source/userguide/examples/GraphicalModels/sparse-gaussian-precision-matrix.ipynb b/docs/source/userguide/examples/GraphicalModels/sparse-gaussian-precision-matrix.ipynb index 9dd3294..b5f37ae 100644 --- a/docs/source/userguide/examples/GraphicalModels/sparse-gaussian-precision-matrix.ipynb +++ b/docs/source/userguide/examples/GraphicalModels/sparse-gaussian-precision-matrix.ipynb @@ -603,46 +603,6 @@ "plt.show()" ] }, - { - "cell_type": "code", - "execution_count": 16, - "id": "8d797c5e", - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAA64AAAEiCAYAAAD564/4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABcA0lEQVR4nO3deVyU5f4//tfMADPsi+ymgrvmQqEirpkomJV2KtMslVzatAxbtE/uFUfrlEczrU5up9M3T4vtomYuJ8NdK01NTdRUcEFAQAZm5vr90Y+pEa77ZtBJbu7X8zzmcfJ+33Pf94zy4r7u7W0QQggQERERERER1VHG670BREREREREREo4cCUiIiIiIqI6jQNXIiIiIiIiqtM4cCUiIiIiIqI6jQNXIiIiIiIiqtM4cCUiIiIiIqI6jQNXIiIiIiIiqtM4cCUiIiIiIqI6jQNXIiIiIiIiqtM4cKUaWbZsGQwGA3Jyctx6n8FgwIwZMzyyTUSkDfU9B2bMmAGDweDWe3JycmAwGLBs2TLPbBSRDt1yyy245ZZbrvdmXFNazoq4uDiMGjXKrffUJk9JP3Q3cDUYDDV6bdy48XpvKhGR008//YR77rkHTZo0gcViQcOGDdGvXz8sWLDgem8aEZGiyoPfstfWrVtrvKyff/4ZM2bMcPtAuqe9+eabmhxcEmmJQQghrvdG/JXee+89lz+vWLEC69atw7///W+X6f369UNUVNRfuWl1mt1uR0VFBcxms1tHwsrKyuDl5QUvLy8Pbh1R/fb999+jT58+aNy4MUaOHIno6GicPHkSW7duxdGjR3HkyJHrvYmKDAYDpk+fXm/PutpsNthsNlgslhq/RwgBq9UKb29vmEwmD24d0fW3bNkypKenY9asWYiPj69ST0tLQ3h4eI2W9dFHH+Hee+/Fhg0bqpxdLS8vBwD4+Phc9Ta7q127dggPD7/mJz5ycnIQHx+PpUuXun328nqzWq0wGo3w9vau8Xtqk6ekH7obTTzwwAMuf966dSvWrVtXZfqVSktL4efn58lNu2ZKSkrg7+9/TZdpMplqtXPF4CG6ei+99BKCg4OxY8cOhISEuNTOnj17fTZKozyRj7U5OGcwGJiPpDsDBgxAp06dPLb86zFgrQ+EECgrK4Ovr+81Xa7ZbHb7PTzZQUp0d6lwTdxyyy1o164ddu3ahV69esHPzw/PP/88APm9WtVdx19QUICJEyeiUaNGMJvNaN68OebMmQOHw6G6DXFxcbj99tuxdu1aJCQkwGKxoG3btvjkk09c5qu8/GbTpk147LHHEBkZiRtuuMFZX716NXr27Al/f38EBgZi4MCB2L9/f5X1HTx4EEOGDEFERAR8fX3RqlUr/N///V+V9fz50pydO3ciNTUV4eHh8PX1RXx8PB566CGX5Vb3fe3ZswcDBgxAUFAQAgIC0Ldv3yqXCVWub8uWLcjIyEBERAT8/f1x11134dy5c6rfH1F9cvToUdx4441VBq0AEBkZ6fLnpUuX4tZbb0VkZCTMZjPatm2LRYsWVXlfZcZs3LgRnTp1gq+vL9q3b+88W/DJJ5+gffv2sFgsSExMxJ49e1zeP2rUKAQEBODXX39Famoq/P39ERsbi1mzZqEmF/KcOnUKDz30EKKiomA2m3HjjTdiyZIlNfo+DAYDxo8fj//85z9o1aqVcxs3b97sMl/lvVI///wz7r//foSGhqJHjx7O+nvvvYfExET4+voiLCwMQ4cOxcmTJ6usb9u2bbjtttsQGhoKf39/dOjQAf/85z+rrOfP1q1bhx49eiAkJAQBAQFo1aqV8/cIIL9v7dtvv3VmdkhICAYNGoQDBw5U+7mOHDmCUaNGISQkBMHBwUhPT0dpaWmNvkOiuuqDDz5AYmIiAgMDERQUhPbt2zt/3pYtW4Z7770XANCnT58qt3ddeY/rxo0bYTAY8N///hczZ85Ew4YNERgYiHvuuQeFhYWwWq2YOHEiIiMjERAQgPT0dFitVpftqUmmxsXFYf/+/di0aZNzm/68HTXdHywoKMCoUaMQHByMkJAQjBw5EgUFBTX63ir3mzZv3oyHH34YDRo0QFBQEEaMGIGLFy9W2d7bb78da9asceb/W2+95da2OhwO/POf/3T+noiIiEBaWhp27tzpsp4/7xtXVFRg5syZaNGiBSwWCxo0aIAePXpg3bp1znmqy1ObzYbZs2ejWbNmMJvNiIuLw/PPP1/l76ryc3333Xfo0qULLBYLmjZtihUrVtToO6S6j4c0JC5cuIABAwZg6NCheOCBB9y+bLi0tBS9e/fGqVOn8PDDD6Nx48b4/vvvMWXKFJw5cwbz5s1TXcbhw4dx33334ZFHHsHIkSOxdOlS3HvvvcjKykK/fv1c5n3ssccQERGBadOmoaSkBADw73//GyNHjkRqairmzJmD0tJSLFq0CD169MCePXsQFxcHAPjxxx/Rs2dPeHt7Y9y4cYiLi8PRo0fxxRdf4KWXXqp2286ePYv+/fsjIiICkydPRkhICHJycqoMrK+0f/9+9OzZE0FBQXj22Wfh7e2Nt956C7fccgs2bdqEpKQkl/knTJiA0NBQTJ8+HTk5OZg3bx7Gjx+PlStXqn5/RPVFkyZNkJ2djX379qFdu3aK8y5atAg33ngj7rzzTnh5eeGLL77AY489BofDgccff9xl3iNHjuD+++/Hww8/jAceeACvvvoq7rjjDixevBjPP/88HnvsMQBAZmYmhgwZgkOHDsFo/ON4p91uR1paGrp27Yq5c+ciKysL06dPh81mw6xZs6TbmJeXh65duzoHoBEREVi9ejVGjx6NoqIiTJw4UfU72bRpE1auXIknnngCZrMZb775JtLS0rB9+/Yq39G9996LFi1a4OWXX3YOql966SVMnToVQ4YMwZgxY3Du3DksWLAAvXr1wp49e5wHCdatW4fbb78dMTExePLJJxEdHY0DBw7gyy+/xJNPPlnttu3fvx+33347OnTogFmzZsFsNuPIkSPYsmWL4mf65ptvMGDAADRt2hQzZszA5cuXsWDBAnTv3h27d+92ZnalIUOGID4+HpmZmdi9ezf+9a9/ITIyEnPmzFH9/oiuh8LCQpw/f95lmsFgQIMGDQD8/vM2bNgw9O3b1/nv+MCBA9iyZQuefPJJ9OrVC0888QTmz5+P559/Hm3atAEA5//LZGZmwtfXF5MnT8aRI0ewYMECeHt7w2g04uLFi5gxYwa2bt2KZcuWIT4+HtOmTXO+tyaZOm/ePEyYMAEBAQHOg/6V+4013R8UQmDQoEH47rvv8Mgjj6BNmzZYtWoVRo4c6dZ3PH78eISEhGDGjBk4dOgQFi1ahOPHjzsH8ZUOHTqEYcOG4eGHH8bYsWPRqlUrt/ZdR48ejWXLlmHAgAEYM2YMbDYb/ve//2Hr1q3Ss+ozZsxAZmYmxowZgy5duqCoqAg7d+7E7t27q+zX/tmYMWOwfPly3HPPPZg0aRK2bduGzMxMHDhwAKtWrXKZ98iRI7jnnnswevRojBw5EkuWLMGoUaOQmJiIG2+80a3vkuogoXOPP/64uPJr6N27twAgFi9eXGV+AGL69OlVpjdp0kSMHDnS+efZs2cLf39/8csvv7jMN3nyZGEymcSJEycUt6tJkyYCgPj444+d0woLC0VMTIy46aabnNOWLl0qAIgePXoIm83mnH7p0iUREhIixo4d67Lc3NxcERwc7DK9V69eIjAwUBw/ftxlXofDUWU9x44dE0IIsWrVKgFA7NixQ/FzXPl9DR48WPj4+IijR486p50+fVoEBgaKXr16VVlfSkqKy3Y89dRTwmQyiYKCAsX1EtUna9euFSaTSZhMJpGcnCyeffZZsWbNGlFeXl5l3tLS0irTUlNTRdOmTV2mVWbM999/75y2Zs0aAUD4+vq65MFbb70lAIgNGzY4p40cOVIAEBMmTHBOczgcYuDAgcLHx0ecO3fOOf3KHBg9erSIiYkR58+fd9mmoUOHiuDg4Go/w58BEADEzp07ndOOHz8uLBaLuOuuu5zTpk+fLgCIYcOGubw/JydHmEwm8dJLL7lM/+mnn4SXl5dzus1mE/Hx8aJJkybi4sWLLvP+OZcq11Pp9ddfFwBcvoMrHTt2TAAQS5cudU5LSEgQkZGR4sKFC85pP/zwgzAajWLEiBFV1vfQQw+5LPOuu+4SDRo0kK6T6Hqp/J1e3ctsNjvne/LJJ0VQUJDL/syVPvzwwyp5VKl3796id+/ezj9v2LBBABDt2rVzycthw4YJg8EgBgwY4PL+5ORk0aRJE5dpNc3UG2+80WXdlWq6P/jpp58KAGLu3LnOeWw2m+jZs2eVrKhO5XecmJjo8lnnzp0rAIjPPvvMOa0y/7Oysmq1rd9++60AIJ544okq2/HnbLxy37hjx45i4MCBip/jyjzdu3evACDGjBnjMt/TTz8tAIhvv/22yufavHmzc9rZs2eF2WwWkyZNUlwvaQMvFZYwm81IT0+v9fs//PBD9OzZE6GhoTh//rzzlZKSArvdXuWSturExsbirrvucv658pKPPXv2IDc312XesWPHutyDum7dOhQUFGDYsGEu6zeZTEhKSsKGDRsAAOfOncPmzZvx0EMPoXHjxi7LVHoIU+XZiC+//BIVFRWqnwX4/ezM2rVrMXjwYDRt2tQ5PSYmBvfffz++++47FBUVubxn3LhxLtvRs2dP2O12HD9+vEbrJKoP+vXrh+zsbNx555344YcfMHfuXKSmpqJhw4b4/PPPXeb98z1KlWc3evfujV9//RWFhYUu87Zt2xbJycnOP1de8XDrrbe65EHl9F9//bXKto0fP97535VnUMvLy/HNN99U+1mEEPj4449xxx13QAjhkk+pqakoLCzE7t27Vb+T5ORkJCYmOv/cuHFjDBo0CGvWrIHdbneZ95FHHnH58yeffAKHw4EhQ4a4rD86OhotWrRw5uOePXtw7NgxTJw4scpl2jXJx88++6xGt4YAwJkzZ7B3716MGjUKYWFhzukdOnRAv3798PXXX1d5z5Wfq2fPnrhw4UKVHCWqKxYuXIh169a5vFavXu2sh4SEoKSkxOXS0WthxIgRLg8ISkpKghCiyu1NSUlJOHnyJGw2m3OaO5lanZruD3799dfw8vLCo48+6nyvyWTChAkT3Pqs48aNc/msjz76KLy8vKpkSHx8PFJTU2u1rR9//LHzoXtXUsvG/fv34/DhwzX+PJXbnZGR4TJ90qRJAICvvvrKZXrbtm3Rs2dP558jIiLQqlWran9/kfbwUmGJhg0bXtVN/ocPH8aPP/6IiIiIaus1eaBK8+bNqwRAy5YtAfx+f1R0dLRz+pVP6asMhVtvvbXaZQcFBQH4Y0dU7fLDK/Xu3Rt33303Zs6ciddffx233HILBg8ejPvvv196M/65c+dQWlqKVq1aVam1adMGDocDJ0+edLmU48rBdGhoKABUuV+DqL7r3LkzPvnkE5SXl+OHH37AqlWr8Prrr+Oee+7B3r170bZtWwDAli1bMH36dGRnZ1e537GwsBDBwcHOP1/581VZa9SoUbXTr/y5MxqNLgehANeMqs65c+dQUFCAt99+G2+//Xa189QkH1u0aFFlWsuWLVFaWopz586p5qMQotplAHDu9B09ehSA+/l433334V//+hfGjBmDyZMno2/fvvjb3/6Ge+65x+VS6z+rPBgny8c1a9ZUebCUUj5WZjxRXdKlSxfFhzM99thj+O9//4sBAwagYcOG6N+/P4YMGYK0tLSrWq87WedwOFBYWOi8fNmdTK1OTfcHjx8/jpiYGAQEBLjUq8sEJVfmWkBAAGJiYqpkcnVPd67pth49ehSxsbEuB9lqYtasWRg0aBBatmyJdu3aIS0tDQ8++CA6dOggfc/x48dhNBrRvHlzl+nR0dEICQmpciLjyr9r4Pds5H5j/cCBq4S7T1a78gi/w+FAv3798Oyzz1Y7f+XO3bVy5fZWHuX/97//7bIDV+lqn9hmMBjw0UcfYevWrfjiiy+wZs0aPPTQQ/jHP/6BrVu3Vgne2pI9yVjoq4sTkZOPjw86d+6Mzp07o2XLlkhPT8eHH36I6dOn4+jRo+jbty9at26N1157DY0aNYKPjw++/vprvP7661XO/sl+vjz5c1e5DQ888ID03i2lnZjaqC4fDQYDVq9eXe1nvdr88vX1xebNm7FhwwZ89dVXyMrKwsqVK3Hrrbdi7dq116z9DfOR6pvIyEjs3bsXa9aswerVq7F69WosXboUI0aMwPLly2u93NpmnbuZWp2/en+wpqrbz/X0tvbq1QtHjx7FZ599hrVr1+Jf//oXXn/9dSxevBhjxoxRfG9NWzEyF+s3DlzdFBoaWuUJb+Xl5Thz5ozLtGbNmqG4uBgpKSm1XteRI0cghHD5Yf3ll18AoMpDOq7UrFkzAL//ElDahsqzJfv27avVNnbt2hVdu3bFSy+9hPfffx/Dhw/HBx98UG0ARUREwM/PD4cOHapSO3jwIIxGY5Wjn0QkV3nmojJ/vvjiC1itVnz++ecuR50rL3291hwOB3799VeXnRm1jIqIiEBgYCDsdvtV5WN1l5r98ssv8PPzk54tqNSsWTMIIRAfH6+4I1aZo/v27XN7W41GI/r27Yu+ffvitddew8svv4z/+7//w4YNG6pdVpMmTQBAmo/h4eHXvI0PUV3k4+ODO+64A3fccQccDgcee+wxvPXWW5g6dWq1V6J5kjuZKtuumu4PNmnSBOvXr0dxcbHLwbPqMkHJ4cOH0adPH+efi4uLcebMGdx2222q763ptjZr1gxr1qxBfn6+22ddw8LCkJ6ejvT0dBQXF6NXr16YMWOGdODapEkTOBwOHD582OUhXHl5eSgoKHBmJ+kD73F1U7Nmzarcn/r2229XOeM6ZMgQZGdnY82aNVWWUVBQ4HL/hMzp06ddnpZWVFSEFStWICEhodqzqH+WmpqKoKAgvPzyy9Xeg1rZUiYiIgK9evXCkiVLcOLECZd5lI5OXbx4sUo9ISEBAKo8nrySyWRC//798dlnn7lcspKXl4f3338fPXr04OVtRNXYsGFDtT+Plff+VF5KVnmk+c/zFhYWYunSpR7btjfeeMP530IIvPHGG/D29kbfvn2rnd9kMuHuu+/Gxx9/XO0Bs5q2u8rOzna5F/bkyZP47LPP0L9/f9Uzmn/7299gMpkwc+bMKt+rEAIXLlwAANx8882Ij4/HvHnzqhywVMrH/Pz8KtPU8jEmJgYJCQlYvny5y7r27duHtWvX1mink0jrKn/2KhmNRucVGJU/O5UHcGraJuZquJOp/v7+1W5TTfcHb7vtNthsNpdWO3a7HQsWLHBrm99++22X/b5FixbBZrNhwIABqu+t6bbefffdEEJg5syZVeZTysYr/34DAgLQvHlzaS4CcGbfld04XnvtNQDAwIEDpe+l+odnXN00ZswYPPLII7j77rvRr18//PDDD1izZg3Cw8Nd5nvmmWfw+eef4/bbb3c+hrukpAQ//fQTPvroI+Tk5FR5z5VatmyJ0aNHY8eOHYiKisKSJUuQl5dXo53QoKAgLFq0CA8++CBuvvlmDB06FBEREThx4gS++uordO/e3bnDOX/+fPTo0QM333wzxo0bh/j4eOTk5OCrr77C3r17q13+8uXL8eabb+Kuu+5Cs2bNcOnSJbzzzjsICgpS3MF68cUXnf0NH3vsMXh5eeGtt96C1WrF3LlzVT8XkR5NmDABpaWluOuuu9C6dWuUl5fj+++/x8qVKxEXF+d8kFz//v2dZysefvhhFBcX45133kFkZGSVq0KuBYvFgqysLIwcORJJSUlYvXo1vvrqKzz//POKZz3//ve/Y8OGDUhKSsLYsWPRtm1b5OfnY/fu3fjmm2+qHfhdqV27dkhNTXVphwOg2h2pKzVr1gwvvvgipkyZgpycHAwePBiBgYE4duwYVq1ahXHjxuHpp5+G0WjEokWLcMcddyAhIQHp6emIiYnBwYMHsX///mp37oDf7+PavHkzBg4ciCZNmuDs2bN48803ccMNN7j0kb3SK6+8ggEDBiA5ORmjR492tsMJDg6utn84kdasXr0aBw8erDK9W7duaNq0KcaMGYP8/HzceuutuOGGG3D8+HEsWLAACQkJzrNtCQkJMJlMmDNnDgoLC2E2m519Vq81dzI1MTERixYtwosvvojmzZsjMjISt956a433B++44w50794dkydPRk5ODtq2bYtPPvmkRg+A+rPy8nL07dvX2cLszTffRI8ePXDnnXeqvrem29qnTx88+OCDmD9/Pg4fPoy0tDQ4HA7873//Q58+fVwe2vdnbdu2xS233ILExESEhYVh586d+Oijj6TzA0DHjh0xcuRIvP322ygoKEDv3r2xfft2LF++HIMHD3Y5u0w68Fc+wrgukrXDufHGG6ud3263i+eee06Eh4cLPz8/kZqaKo4cOVLlkd9C/N6SZsqUKaJ58+bCx8dHhIeHi27duolXX3212jYWf9akSRMxcOBAsWbNGtGhQwdhNptF69atxYcffugyX+Xjz2VtaTZs2CBSU1NFcHCwsFgsolmzZmLUqFEubSSEEGLfvn3irrvuEiEhIcJisYhWrVqJqVOnVllPZTuc3bt3i2HDhonGjRsLs9ksIiMjxe23315luaimfdDu3btFamqqCAgIEH5+fqJPnz4uLTmUPlflo+2reww+UX21evVq8dBDD4nWrVuLgIAA4ePjI5o3by4mTJgg8vLyXOb9/PPPRYcOHYTFYhFxcXFizpw5YsmSJS4/v0L8kTFXAiAef/xxl2mVrVteeeUV57SRI0cKf39/cfToUdG/f3/h5+cnoqKixPTp04Xdbq+yzCtzIC8vTzz++OOiUaNGwtvbW0RHR4u+ffuKt99+W/X7qNzG9957T7Ro0UKYzWZx0003VcmFyrYKsrY0H3/8sejRo4fw9/cX/v7+onXr1uLxxx8Xhw4dcpnvu+++E/369ROBgYHC399fdOjQQSxYsKDKeiqtX79eDBo0SMTGxgofHx8RGxsrhg0b5tJiorp2OEII8c0334ju3bsLX19fERQUJO644w7x888/1+hzXZnTRHWFUjucP/8cfPTRR6J///4iMjJS+Pj4iMaNG4uHH35YnDlzxmV577zzjmjatKkwmUwu+wSydjg13Xeq7merppmam5srBg4cKAIDAwUAl+2o6f7ghQsXxIMPPiiCgoJEcHCwePDBB8WePXvcaoezadMmMW7cOBEaGioCAgLE8OHDXVpsCSHPf3e21WaziVdeeUW0bt1a+Pj4iIiICDFgwACxa9cul/X8ed/4xRdfFF26dBEhISHC19dXtG7dWrz00ksuy70yT4UQoqKiQsycOVPEx8cLb29v0ahRIzFlyhRRVlZWo8915b8L0i6DELxbuS6Ki4tDu3bt8OWXX17vTSEiqmLUqFH46KOPUFxc/Jev22Aw4PHHH3e5TJmISM+WLVuG9PR07NixQ/HJzURaxntciYiIiIiIqE7jwJWIiIiIiIjqNA5ciYiIiIiIqE7jwLWOysnJ4f2tdN1s3rwZd9xxB2JjY2EwGPDpp5+qvmfjxo24+eabYTab0bx5cyxbtqzKPAsXLkRcXBwsFguSkpKwffv2a7/x9JdYtmzZdbm/Ffij7Q7R1WLWUX0xatQoCCF4fytVq75kHQeuRFRFSUkJOnbsiIULF9Zo/mPHjmHgwIHo06cP9u7di4kTJ2LMmDEu7UJWrlyJjIwMTJ8+Hbt370bHjh2RmpqKs2fPeupjEBEpYtYRkR7Ul6zjU4WJSJHBYMCqVaswePBg6TzPPfccvvrqK+zbt885bejQoSgoKEBWVhYAICkpCZ07d3aeKXM4HGjUqBEmTJiAyZMne/QzEBGpYdYRkR5oOet4xpVIJ6xWK4qKilxeVqv1miw7OzsbKSkpLtNSU1ORnZ0N4PeG6Lt27XKZx2g0IiUlxTkPEdG1wKwjIj3QY9Z51XTG1mNja7WCvQuO1up9RHphsfjW6n3u/kwObTgOM2fOdJk2ffp0zJgxo1br/7Pc3FxERUW5TIuKikJRUREuX76Mixcvwm63VzvPwYMHr3r915LS92pQONT344uHFJdbHuiQ1i45TkproWgmrdl/OC+teQX6yLelmbe0BgCFjhxpLSZfvj22QPmvlIs+x6W1YGNTaU1A/r0BQLHjN2kt9FKMtGawyy82svua5Cs0KG4OrOYyac3rh8vSmrF9A2ntiHW9tNbIt7O05lvqL60BgMMs/5wXjUektWCj/N+Ad4FdWvOKDlTcHhlmnWcofa83NJXnx5dPHfDE5hDVG9yv81zW1XjgSkR1i9IgqjpTpkxBRkaGyzSz2XwNt4iI6Npj1hGRHjDr1HHgSqRRRqPKqZ8rmM1mjwVadHQ08vLyXKbl5eUhKCgIvr6+MJlMMJlM1c4THR3tkW0iovqBWUdEesCsU8d7XIk0ymBw7+VJycnJWL/e9XLGdevWITk5GQDg4+ODxMREl3kcDgfWr1/vnIeIqDrMOiLSA2adOp5xJdIoowcPOxUXF+PIkT/ucTt27Bj27t2LsLAwNG7cGFOmTMGpU6ewYsUKAMAjjzyCN954A88++yweeughfPvtt/jvf/+Lr776yrmMjIwMjBw5Ep06dUKXLl0wb948lJSUID093XMfhIg0j1lHRHrArFPHgSuRRhncvKTEHTt37kSfPn2cf668h2LkyJFYtmwZzpw5gxMnTjjr8fHx+Oqrr/DUU0/hn//8J2644Qb861//QmpqqnOe++67D+fOncO0adOQm5uLhIQEZGVlVbmxn4joz5h1RKQHzDp1Ne7jyqcKE3lGbZ8+d9PEG9yaf888+dNX6Q98qjCfKsynCtetpwoz6zyDTxUm8gzu13kOz7gSaZS7T58jItIiZh0R6QGzTp3HB64JE+RHZgGekSWqLaOn78zXKaVfHELh5F/751spLnfHkP9Ja5Et5UdZiyLzpTWvBPkZtc3n35XWvtnwubQGAK0a3iit3dnyOWnNLuRnGxv8Jr90yNZYvi0Gofyb3Ab5WUxbsPzMsvGy/MygOHFJWvMKVz6S7lcqX6fdzyatWY2l0lpzS1/5+1AsrRkqlC+oMp0vkdb8G8vPxpWIM9KaJUR+5jhAcWvkmHWeoXRW9bdfy6U17tcReQazTh3PuBJpFI/MEZEeMOuISA+Ydeo4cCXSKE8+fY6IqK5g1hGRHjDr1HHgSqRRPDJHRHrArCMiPWDWqePAlUijTCbeC0FE9R+zjoj0gFmnjgNXIo3iJSVEpAfMOiLSA2adOg5ciTTKk42qiYjqCmYdEekBs04dB65EGsUjc0SkB8w6ItIDZp26Gg9clfpyqfX0UqL0XvYCI5LjTfye8eOLh6Q1tV6tSjr/t6e09vPz+6W1g0XrpLXWQf2ktf9s/Le0VnpZ3qMRAGb0mSetmS7Le5EuOTZeWnugZIa0drlRgbQW6JD3EwWAX4u2SWtRfheltaYFHaW1wpYV8u05ofxr0xEp7/Na1kL+vkL7r9JaVEVbac3HEiit2UKU+7gazfIQKcdZaU0IeQ9cL4NFcZ21wazzjC+fOiCtcb+O6K/HrFPHM65EGsVG1USkB8w6ItIDZp06DlyJNIpH5ohID5h1RKQHzDp1HLgSaRTvhSAiPWDWEZEeMOvUceBKpFF8+hwR6QGzjoj0gFmnjgNXIo3ikTki0gNmHRHpAbNOHQeuRBrFe/iJSA+YdUSkB8w6dRy4EmmUkZeUEJEOMOuISA+YdequycCVPV6J/np8+pxnlAc6pLUdQ/4nrSn1aVXT9uUbpbWfp/4srTlOyvtp/v3uD6S1p1feq7g9uwpWSGvJZfL3pjd/U1oznSuT1rzLzNLa3vKPpDUAcCj0FI3YHC2t7ezypbTWAfLPeLrhXsXtiYK8P+zZcnnfzEbGLtLaAUeWtNbSkCatfXN6rrQGACGWKGmtq3hQWhNKO1cG+c8PatnilVn31+N+HdFfj1mnjmdciTSKR+aISA+YdUSkB8w6dRzbE2mUl8no1stdCxcuRFxcHCwWC5KSkrB9+3bpvLfccgsMBkOV18CBA53zjBo1qko9LU1+toiICGDWEZE+MOvU8YwrkUZ58sjcypUrkZGRgcWLFyMpKQnz5s1DamoqDh06hMjIyCrzf/LJJygvL3f++cKFC+jYsSPuvdf1ksu0tDQsXbrU+WezWX6JKBERwKwjIn1g1qnjGVcijTIajG693PHaa69h7NixSE9PR9u2bbF48WL4+flhyZIl1c4fFhaG6Oho52vdunXw8/OrEnBms9llvtDQ0Fp/fiLSB2YdEekBs04dB65EGmU0Gtx6Wa1WFBUVubysVmuV5ZaXl2PXrl1ISUn507qMSElJQXZ2do227d1338XQoUPh7+/vMn3jxo2IjIxEq1at8Oijj+LChQtX9yUQUb3HrCMiPWDWqePAlUij3A24zMxMBAcHu7wyMzOrLPf8+fOw2+2IinJ94mhUVBRyc3NVt2v79u3Yt28fxowZ4zI9LS0NK1aswPr16zFnzhxs2rQJAwYMgN0ufyosERGzjoj0gFmnjve4EmmU0c1O1VOmTEFGRobLNE/ci/Duu++iffv26NLFtb3H0KFDnf/dvn17dOjQAc2aNcPGjRvRt2/fa74dRFQ/MOuISA+Ydeo8PnBV68tV235g7AVGemc0unfBhNlsrlGghYeHw2QyIS8vz2V6Xl4eoqPlfTEBoKSkBB988AFmzZqlup6mTZsiPDwcR44cqVM7c5ccJ6W1yJY3SGs/P79fcblKvVoV3ze7rbS2/a6N0trlmIvSWnq/0YrrbBHcW1o7F3xGWgvcEyCtGZuGSGsl63+V1s53ypHWAKBrlLzfqMFULq3lFPworRVXnJPWTEYfxe0JCY2X1mK9EqU1Y7G8/2m4b5z8fRXybekWM1JeBHC06DtpzVQoX3BxtPx7NRnkzVqVvzk5Zl3dwv06Is9g1qnjpcJEGuXuJSU15ePjg8TERKxfv945zeFwYP369UhOTlZ874cffgir1YoHHnhAdT2//fYbLly4gJiYmBpvGxHpD7OOiPSAWaeOA1cijTIaDG693JGRkYF33nkHy5cvx4EDB/Doo4+ipKQE6enpAIARI0ZgypQpVd737rvvYvDgwWjQoIHL9OLiYjzzzDPYunUrcnJysH79egwaNAjNmzdHampq7b8EIqr3mHVEpAfMOnW8x5VIozzZ7+u+++7DuXPnMG3aNOTm5iIhIQFZWVnOG/tPnDhR5ZKWQ4cO4bvvvsPatWurLM9kMuHHH3/E8uXLUVBQgNjYWPTv3x+zZ89mf0MiUsSsIyI9YNap48CVSKPcvRfCXePHj8f48eOrrW3cuLHKtFatWkEIUe38vr6+WLNmzbXcPCLSCWYdEekBs04dB65EGmVw8zIRIiItYtYRkR4w69Rx4EqkUZ68pISIqK5g1hGRHjDr1HHgSqRRDDgi0gNmHRHpAbNO3XUfuCr15mIvMCI5o4EPBfeEUMjzoygyX1o7WLROcbk/T/1ZWlPq1aqky6pbpLXtfX6Q1mKDblJcboVCc1C7KJPWfMJN8mUeL5LWAnrHSWtdA5R7kfo7IqS1g0mrpbW/XZgsrf0WfFBaa1zQTnF7vjrzurQ2oKL6e4sAoPxsibTWoLP834ehwi6thVxqIK0BQIutN0tr27t/Jq3FCXk/2lDRXHGdtcGs0xbu1xHVDrNO3XUfuBJR7XiZGHBEVP8x64hID5h16jhwJdIoXlJCRHrArCMiPWDWqePAlUijPP3YdCKiuoBZR0R6wKxTx4ErkUYZ+dh0ItIBZh0R6QGzTh0HrkQaxUtKiEgPmHVEpAfMOnUcuBJpFAOOiPSAWUdEesCsU8eBK5FG8bHpRKQHzDoi0gNmnbo6PXBlLzAiOR6Z8wz7D+elNa8Ef2mtdVA/xeU6Tsr7bW6/a6O0ptSrVUmXJzpKa/unH1B87+ngfdJauFcbac0RY5EvtFDe//WCf460FmJoIV8mgMvGC9JaM/9bpTXDZYWfn7nyZRZ3P6u4PZ173SOtOX6V/xvwDveTL7TUJi2V+8mXabDJ++oCQPHPJ+W1zuekNW9DgHyd5UJxnbXBrKs/uF9HJMesU1enB65EJGc0KO+UEhHVB8w6ItIDZp06DlyJNIqPTSciPWDWEZEeMOvUceBKpFEmHpkjIh1g1hGRHjDr1HHgSqRRRiMDjojqP2YdEekBs04dB65EGsV7IYhID5h1RKQHzDp1HLgSaRTvhSAiPWDWEZEeMOvUceBKpFHeJp/rvQlERB7HrCMiPWDWqdPswPWv7gWmtk6iv5qnG1UvXLgQr7zyCnJzc9GxY0csWLAAXbp0qXbeZcuWIT093WWa2WxGWdkfvTuFEJg+fTreeecdFBQUoHv37li0aBFatFDu0/lX8wqU/+LYfP5dae0/G/+tuNy/3/2BtHY55qK0tr3PD9KaUq9WJTfOlPdiBYAfJ+yR1o413imtBflESWsNWsn/niPPyfuCOswV0hoAlAXlS2uhjjBpTZjky71hQh9prShQ/ncFAIXWY9JaRGgjae3z0rnS2m3mF6S1fJu8J2+ov/LvtIbDk6Q1//DW0poPgqQ1g92huM7aYNbpA/frSO+Ydep4TppIo4xGk1svd6xcuRIZGRmYPn06du/ejY4dOyI1NRVnz56VvicoKAhnzpxxvo4fP+5Snzt3LubPn4/Fixdj27Zt8Pf3R2pqqksIEhFdiVlHRHrArFPHgSuRRpkMJrde7njttdcwduxYpKeno23btli8eDH8/PywZMkS6XsMBgOio6Odr6ioP86+CSEwb948vPDCCxg0aBA6dOiAFStW4PTp0/j0009r+xUQkQ4w64hID5h16jhwJdIoo8Ho1stqtaKoqMjlZbVaqyy3vLwcu3btQkpKyh/rMhqRkpKC7Oxs6fYUFxejSZMmaNSoEQYNGoT9+/c7a8eOHUNubq7LMoODg5GUlKS4TCIiZh0R6QGzTh0HrkQa5e4lJZmZmQgODnZ5ZWZmVlnu+fPnYbfbXY6sAUBUVBRyc3Or3ZZWrVphyZIl+Oyzz/Dee+/B4XCgW7du+O233wDA+T53lklEBDDriEgfmHXqNPtwJiK9c7ff15QpU5CRkeEyzWw2X5NtSU5ORnJysvPP3bp1Q5s2bfDWW29h9uzZ12QdRKRPzDoi0gNmnToOXIk0yuTmjflms7lGgRYeHg6TyYS8vDyX6Xl5eYiOjq7Rury9vXHTTTfhyJEjAOB8X15eHmJiYlyWmZCQUMNPQER6xKwjIj1g1qnjpcJEGuXuvRA15ePjg8TERKxfv945zeFwYP369S5H35TY7Xb89NNPzjCLj49HdHS0yzKLioqwbdu2Gi+TiPSJWUdEesCsU1cvz7h6oheY2nvZC4z+au4+Ct0dGRkZGDlyJDp16oQuXbpg3rx5KCkpcfb0GjFiBBo2bOi8l2LWrFno2rUrmjdvjoKCArzyyis4fvw4xowZA+D3J9NNnDgRL774Ilq0aIH4+HhMnToVsbGxGDx4sMc+R22UN/OW1r7Z8Lm0Vnq5XHG5T6+8V1pL7zdaWosNukla2z9d3sNTrVerkg4L5Ovc+M810pqvIUJaM/xWLK3ZFbbFER6oUAW8DfIesAabkNaMZfK12kvlPV5L/PKkNUC5l22p2SattfSV/6L3KpD/24o8e4O09qW/vDcsANwRNU1ayyuT/9tqZuwtrdl+kffV9e4q/7tSwqwj7teRHjDr1NXLgSuRHrh7L4Q77rvvPpw7dw7Tpk1Dbm4uEhISkJWV5bwJ/8SJEzAa/zjad/HiRYwdOxa5ubkIDQ1FYmIivv/+e7Rt29Y5z7PPPouSkhKMGzcOBQUF6NGjB7KysmCxWDz2OYhI+5h1RKQHzDp1BiGE/JD0n5SVXfbYRvyVrubInBIemaPaslh8a/W+/x56xq35h7R6pVbr0Zviy/IzRi9uGCStHfvtnOJyfbzlv5CUzrj2CHpEWjMVyc8MXs0ZVyW1PePq92vt7kxxxCmfcS02yp9eGFghP/vpdaFqy4BKSmdcc5scl9YAwMvoI60FGORnR09clrcPaH25p7RmO1sqrX3p/09pDVA+43rUsUlaUzrjaj8g//nx7Sr//EqYdZ7B/Tpl3K+j2uJ+nefwjCuRRrnbfJqISIuYdUSkB8w6dRy4EmmUJy8pISKqK5h1RKQHzDp1HLgSaZSXSf4QISKi+oJZR0R6wKxTx4ErkUa52++LiEiLmHVEpAfMOnUcuBJpFC8pISI9YNYRkR4w69TpbuDKXmBUX7jTfJpqrtCRI621anijtDajzzzF5e4qWCGttQiWP6W1wih/wu3p4H3S2o8T9khrSn1a1dzyZKq0tveRndKaXcj7phpahUprJ23bFLengY88ex3e8p8RYZc/UP/Cpl+ktaD05orbY7EFSWtKfWVbF3ST1hxWeR9X65lL0lpoh0hpDQAMdoe0FmZpIq0V4KR8nc1iFddZG8w6UsL9OqovmHXqdDdwJaoveGSOiPSAWUdEesCsU8eBK5FGMeCISA+YdUSkB8w6dRy4EmmUgQFHRDrArCMiPWDWqePAlUijeGSOiPSAWUdEesCsU8eBK5FGGcGAI6L6j1lHRHrArFPHgSuRRvHIHBHpAbOOiPSAWaeOA1cijWLAEZEeMOuISA+Ydeo4cP0Ttb5cV9MPjOha4038nhGTL/85v7Plc9Ka6bJNcbnJZfdKa+eCz0hrdlEmrYV7tZHWjjWW91Td+M810hqg3KtVUdswaclYLu8ZWuZVIn9fhWf62llv8JHWwobJ+9yWolBxuQb5x1TmI/95FqFmac0RJ+/x2s3wqOIqHfLWuggrbSytXfI/L61d9D8trUWileL2yDDrqLa4X0dawqxTx4ErkUbxXggi0gNmHRHpAbNOHQeuRBplNHjmTBQRUV3CrCMiPWDWqePAlUijvIzySx2JiOoLZh0R6QGzTh0HrkQaxXshiEgPmHVEpAfMOnUcuBJpFO+FICI9YNYRkR4w69TxYmoijTIaTG693LVw4ULExcXBYrEgKSkJ27dvl877zjvvoGfPnggNDUVoaChSUlKqzD9q1CgYDAaXV1pamtvbRUT6wqwjIj1g1qnjwJVIozwZcCtXrkRGRgamT5+O3bt3o2PHjkhNTcXZs2ernX/jxo0YNmwYNmzYgOzsbDRq1Aj9+/fHqVOnXOZLS0vDmTNnnK//9//+X60/PxHpA7OOiPSAWafOIIQQNZmxrOyyRzekPlPqE6bWY4zqP4vFt1bvO1Sw1q35W4X0r/G8SUlJ6Ny5M9544w0AgMPhQKNGjTBhwgRMnjxZ9f12ux2hoaF44403MGLECAC/H5krKCjAp59+6tZ2/9XKC+Q9RQt8Tkhrqw7MUVxuevM3pTXbPnlfTJ9wP2nNESOvXfA6LK35GaKkNQAw/6zQjFShV2vCE02ltV0LDkhrBXZ5DkbYW8u3BYDNR76tVnFRWvOvaCCtlXvLe+caDd6K22O+IN8eW4j8wRsOk/xXcZEjR1oTkDdjDbXXvkelwSbfHq8Cee9YW6j8M/qE+NdqW5h1nsH9utrjfh0p4X6d5/CMK5FGGQwmt15WqxVFRUUuL6vVWmW55eXl2LVrF1JSUpzTjEYjUlJSkJ2dXaNtKy0tRUVFBcLCXAc5GzduRGRkJFq1aoVHH30UFy5cuLovgYjqPWYdEekBs04dB65EGmWEya1XZmYmgoODXV6ZmZlVlnv+/HnY7XZERbmelYuKikJubm6Ntu25555DbGysS0impaVhxYoVWL9+PebMmYNNmzZhwIABsNvlZ4yIiJh1RKQHzDp1fKowkUa526h6ypQpyMjIcJlmNpuv5SYBAP7+97/jgw8+wMaNG2GxWJzThw4d6vzv9u3bo0OHDmjWrBk2btyIvn37XvPtIKL6gVlHRHrArFPHgSuRRrl7Y77ZbK5RoIWHh8NkMiEvL89lel5eHqKjoxXf++qrr+Lvf/87vvnmG3To0EFx3qZNmyI8PBxHjhzhzhwRSTHriEgPmHXqeKkwkUZ56ulzPj4+SExMxPr1653THA4H1q9fj+TkZOn75s6di9mzZyMrKwudOnVSXc9vv/2GCxcuICYmpsbbRkT6w6wjIj1g1qnjwJVIowwwufVyR0ZGBt555x0sX74cBw4cwKOPPoqSkhKkp6cDAEaMGIEpU6Y4558zZw6mTp2KJUuWIC4uDrm5ucjNzUVxcTEAoLi4GM888wy2bt2KnJwcrF+/HoMGDULz5s2Rmpp67b4UIqp3mHVEpAfMOnW8VJhIo2rTfLqm7rvvPpw7dw7Tpk1Dbm4uEhISkJWV5byx/8SJEzAa/zjutWjRIpSXl+Oee+5xWc706dMxY8YMmEwm/Pjjj1i+fDkKCgoQGxuL/v37Y/bs2R65H4OI6g9mHRHpAbNOHfu4XmfsBUa17fd1pvQHt+aP8etYq/XoTV7pz9Jag9/k/U+tuZcUl2uJDZLWhL/8GGLF8SL5QhWumTG1kvdbNfxWLH8jAPvlCvkqW4RKaxVm+fsSJ7SR1rbO3ymt+Qr5+gDAq0C+zvwQ+dMSg60NpbVii/xx/n6GSMXtEZD3cfUqlteK/eU9Z02Q90b1Msh7+aodka8Q8p7FZhEg354im7RWHiJfn58lUHF7ZJh1nsH9Os/gfh1xv85zeMaVSKPcvUyEiEiLmHVEpAfMOnUcuBJplIG3qBORDjDriEgPmHXqOHAl0igDDNd7E4iIPI5ZR0R6wKxTx4ErkWbxyBwR6QGzjoj0gFmnhgNXIo3ikTki0gNmHRHpAbNOHQeuRBplMPDIHBHVf8w6ItIDZp06DlyJNItH5ohID5h1RKQHzDo17ONah7EXmD7Utt/XhcvH3Jq/gW98rdajN2VlZbV6X7E4pVgPKouW1krW/yqtBfSOk9Yu+OdIa5HnGklr9uJyaQ0AHI3lPTytXvLen5ccJ6U1f6P883d9opO0tmX+FmkNALwh39ZyyHvg+kDeV9e3RP4zWR5gV9wen2J5OwNhku+UGM/J/93ZGst7tV4W56U1syFEWgMAu5D/O/AyyL8Dk8I/n8ve8u882Ff+b0AJs84zuF/31+N+nT5wv85zeMaVSKN4LwQR6QGzjoj0gFmnjgNXIs3ivRBEpAfMOiLSA2adGg5ciTSKR+aISA+YdUSkB8w6dRy4EmmUgUfmiEgHmHVEpAfMOnUcuBJpFo/MEZEeMOuISA+YdWo4cCXSKB6ZIyI9YNYRkR4w69Rx4EqkUbwXgoj0gFlHRHrArFPHPq4apdQLDGA/MC2pbb+vS5cvuDV/oG+DWq1Hby6XlUprBiE/GmpwKEfpD8UfSWvnL+dIa12jRkprfoZIac2roEJaswcqH7M8adsmrRkN8u/gBsj7sQov+S/kYkOutNb9ie7SGgD8+I8j0prxsrzn6vnA49JamFXeG6/YovxzF3DYLK3ZWvpLa8Z9hfL3tZO/r1wUyJcJH2kNALwM8v6wPkes0lppM/n3KiCvhfrKewsrYdZ5Bvfr6hbu19Uf3K/zHJ5xJdIoHpkjIj1g1hGRHjDr1HHgSqRVCme+iIjqDWYdEekBs04VvyEijTK4+T93LVy4EHFxcbBYLEhKSsL27dsV5//www/RunVrWCwWtG/fHl9//bVLXQiBadOmISYmBr6+vkhJScHhw4fd3i4i0hdmHRHpAbNOHQeuRBplgNGtlztWrlyJjIwMTJ8+Hbt370bHjh2RmpqKs2fPVjv/999/j2HDhmH06NHYs2cPBg8ejMGDB2Pfvn3OeebOnYv58+dj8eLF2LZtG/z9/ZGamoqysrKr+h6IqH5j1hGRHjDr1PHhTBrFm/jrj9rexF9aVuzW/H6WgBrPm5SUhM6dO+ONN94AADgcDjRq1AgTJkzA5MmTq8x/3333oaSkBF9++aVzWteuXZGQkIDFixdDCIHY2FhMmjQJTz/9NACgsLAQUVFRWLZsGYYOHerWZ/EkPpyJD2fiw5nq1sOZmHWewf26uoX7dfUH9+s8l3U840qkUZ46MldeXo5du3YhJSXFOc1oNCIlJQXZ2dnVvic7O9tlfgBITU11zn/s2DHk5ua6zBMcHIykpCTpMomIAGYdEekDs04dH85EpFHu3t9gtVphtbqeRTGbzTCbXc8QnT9/Hna7HVFRUS7To6KicPDgwWqXnZubW+38ubm5znrlNNk8RETVYdYRkR4w69Rx4KpRapeMKF1ywstN6gv3Ai4zMxMzZ850mTZ9+nTMmDHjGm6T9hU7fpPWbJBfWvdrkfzyWgBwCPkllF2jHpTW/B0R0tplo/yy1bKgfGnN26B8eVEDH+VL1mRsBoe0Zs6X17xD5dujdCkwAHSY1Fxa2/fiIWktrFx+OTBM8p8tiyFMcXvKW9qkNS+YpDVji1BpzSEuSWuBJeHSmsGufPl6RbB8e87F/yqtBSNOWvN21O4SOWXMOqr/uF9HzDp1HLgSaVWN7k7/w5QpU5CRkeEy7cqjcgAQHh4Ok8mEvLw8l+l5eXmIjo6udtnR0dGK81f+f15eHmJiYlzmSUhIcO+DEJG+MOuISA+Ydap4jyuRRhmEcOtlNpsRFBTk8qou4Hx8fJCYmIj169c7pzkcDqxfvx7JycnVbktycrLL/ACwbt065/zx8fGIjo52maeoqAjbtm2TLpOICGDWEZE+MOvU8YwrkVa5eWTOHRkZGRg5ciQ6deqELl26YN68eSgpKUF6ejoAYMSIEWjYsCEyMzMBAE8++SR69+6Nf/zjHxg4cCA++OAD7Ny5E2+//TYAwGAwYOLEiXjxxRfRokULxMfHY+rUqYiNjcXgwYM990GISPuYdUSkB8w6VRy4EmmVBwPuvvvuw7lz5zBt2jTk5uYiISEBWVlZzpvwT5w4AaPxjws2unXrhvfffx8vvPACnn/+ebRo0QKffvop2rVr55zn2WefRUlJCcaNG4eCggL06NEDWVlZsFgsnvsgRKR9zDoi0gNmnSr2ca2neBO/dtS235f1krzfaHXMgfKejfSHc6W/SGuKD2cqrP3DmW4MHSCt+TuipLXLRvkDmMpE7R/OZKzlMU2zQf6AIaWHM5WEynuGWirkPUyB2j+cyWFWuFPGKH9Ahs1H/jkAwAGlhzPJe7warfJfxWU+8ocz+ZbIv5+reThTvr36J00CQLAxTlpTejiT2Z9ZV5dwv05buF+nHdyv8xyecSXSKIUHuBIR1RvMOiLSA2adOg5cibSqZhdLEBFpG7OOiPSAWaeKlwrrEC83qVtqe0lJ+cUSt+b3CVW+5JJ+Z8uTX5ZpC/aW1n6zblVcbsTm6h85DwAGk/yy1WNJ+6S1Zv63SmvedvllqWqXkDq8a/fA+cs4J63ZhPx3iMEgv2Q1pCiyVtsCAO1eaCWt/TBPnnX5Bvnl4g0qWiiu0+Ej/+6MVvnl4qXmAmnN1yHvHVtqPC+t+Zcp95xVahkoFHrZWr3k2eN/Uf7vzhQbpLw9Esw6z+B+Xf3B/bq6hft1nsN2OERERERERFSn8VJhIo0y8JISItIBZh0R6QGzTh0HrkRaxXwjIj1g1hGRHjDrVHHgSqRVDDgi0gNmHRHpAbNOFQeuRFrFS0qISA+YdUSkB8w6VRy4EmmUgflGRDrArCMiPWDWqePAlUirGHBEpAfMOiLSA2adKvZxJRfsBfbXq22/L1uuvN9odbyiA2u1Hr2xny6S13zl/Ua9CsoVl7vT70tpLafgR2ntb4bJ0prBT37sUakPp7FM3k8UAIRCn1frDT7Smle5Qoc1hWWWWUqltVKRJ18mgLDyeGlNKPSj7ThRnnU/vnJYWiszFytuj++v8u+9vJm8x6nPJfm22gPlf89lIl9aCyhU7puq1JfY+5T878QWI8+sMmOhtBbsK+9lrIRZ5xncr9MH7tf99bhf5zk840qkVbwXgoj0gFlHRHrArFPFgSuRRvFeCCLSA2YdEekBs06dwnVdRERERERERNcfz7gSaZWDh+aISAeYdUSkB8w6VRy4EmkULykhIj1g1hGRHjDr1HHgSqRVvImfiPSAWUdEesCsU8WBK5FWMd+ISA+YdUSkB8w6VezjSjXGXmCeUdt+X45jF92a3xgfWqv16E35xRJpTZyQ91grbVmhuFxfREhr3194U1prGpwkX+jcC9LSDRP6SGvinHKeX9j0i7QWNuwmaa3cYpXXIO9/GlgSJq0p9aMFACjUL3jJ+7E2KG8urXV4poW09v387xU3J6AiUlor8Zb/fVmFvP9pqJBnr5C3FoZ3vvK/ydNBh6S12KJW0lpFqLyXr/dFeT9jU6xyX1kZZp1ncL+OuF/nGdyv8xw+VZhIo4QQbr08JT8/H8OHD0dQUBBCQkIwevRoFBfLByn5+fmYMGECWrVqBV9fXzRu3BhPPPEECgtdd9wNBkOV1wcffOCxz0FEdROzjoj0gFmnjpcKE2mV43pvwO+GDx+OM2fOYN26daioqEB6ejrGjRuH999/v9r5T58+jdOnT+PVV19F27Ztcfz4cTzyyCM4ffo0PvroI5d5ly5dirS0NOefQ0JCPPlRiKguYtYRkR4w61Rx4EqkUaIOPDb9wIEDyMrKwo4dO9CpUycAwIIFC3Dbbbfh1VdfRWxsbJX3tGvXDh9//LHzz82aNcNLL72EBx54ADabDV5ef8RSSEgIoqOjPf9BiKjOYtYRkR4w69TxUmEirRLCvZcHZGdnIyQkxBluAJCSkgKj0Yht27bVeDmFhYUICgpyCTcAePzxxxEeHo4uXbpgyZIlHr00hojqKGYdEekBs04Vz7gSaZS7R+asViusVteH55jNZpjN5lpvQ25uLiIjXR9G4+XlhbCwMOTm5tZoGefPn8fs2bMxbtw4l+mzZs3CrbfeCj8/P6xduxaPPfYYiouL8cQTT9R6e4lIe5h1RKQHzDp1PONKpFUO4dYrMzMTwcHBLq/MzMxqFz158uRqb6L/8+vgwYNX/RGKioowcOBAtG3bFjNmzHCpTZ06Fd27d8dNN92E5557Ds8++yxeeeWVq14nEWkMs46I9IBZp4pnXIk0yt3LK6ZMmYKMjAyXabKjcpMmTcKoUaMUl9e0aVNER0fj7NmzLtNtNhvy8/NV72G4dOkS0tLSEBgYiFWrVsHb21tx/qSkJMyePRtWq/WqjiYSkbYw64hID5h16jhwpRpT6uml1AtM7b1US24+fc6dy0ciIiIQESHvO1opOTkZBQUF2LVrFxITEwEA3377LRwOB5KS5P1Hi4qKkJqaCrPZjM8//xwWi0V1XXv37kVoaKjnd+QU2oZ6hct7swWeUI7T0w33Smsmo7wvZuOCdtJacfez0lpRoLwfXIlfnrQGAEHp8h6npZD3G/UxBEprfpD3Ny0PkPd/VepvCgAWg7wHbAOrvB9rmVnek1epV2u3J7opbs+2+bulNZuQ980MMN4grRlL7NKaUp/b/BDly7oCDVUfsuFcbrG8B6whSL4zUhEm/7es0HJWGbOOyCO4X1fHMOtUceBKpFF14eEdbdq0QVpaGsaOHYvFixejoqIC48ePx9ChQ51Pnjt16hT69u2LFStWoEuXLigqKkL//v1RWlqK9957D0VFRSgqKgLwe7CaTCZ88cUXyMvLQ9euXWGxWLBu3Tq8/PLLePrpp6/nxyWi64BZR0R6wKxTx4ErkUYJW91o+PWf//wH48ePR9++fWE0GnH33Xdj/vz5znpFRQUOHTqE0tJSAMDu3budT6Zr3tz1zN6xY8cQFxcHb29vLFy4EE899RSEEGjevDlee+01jB079q/7YERUJzDriEgPmHXqDKKGw/uyMvklTkS8pKT2LBb55adKynaecm89nRrWaj16U15QIq2ZSmzSmtKllQBwuuERae1Y0U5praf9QWmteL/8UmFbb/mlNyV2lUuFTY2kNaFwLZPSpcJGyC8vtcMzlwp7WeXPHyzzkV8qbIP8993VXCpcLoqkNV+j/FJqc4n8GLPSpcKFPsoZ4W3wl9YCTsov87LH+sm3x0u+PTW5dKw6zDrP4H4dKeF+Xe1xv85zeMaVSKPqQqNqIiJPY9YRkR4w69Rx4EqkVXXjihIiIs9i1hGRHjDrVHHgSqRRdeEmfiIiT2PWEZEeMOvUceBKpFW8pISI9IBZR0R6wKxTxYcz0V9C6SZ/vd/gX9ub+Ev/d9yt+f16NqnVevTm0uUL0prfRYUHDPkqd6kU3vIHBVlN8of2bDqzSFrrHH2PtFZYflpaC/KJktYAIMwu739qULiUyaTwgKryBvLvzqtYvlBxWv6wLAAobyl/CJWXkD8MyPSr/DsXjQKktcve8vcBQNITN0trP849LF+uRf6wKCWWXxSKzYMV33vZmC+t7b6wUlrrbRgtX6hB/nAmU2yQ4vbIMOs8g/t1dDW4XyfH/TrP4RlXIo3iJSVEpAfMOiLSA2adOg5cibSKN/ETkR4w64hID5h1qjhwJdIoPjadiPSAWUdEesCsU8eBK5FW8ZISItIDZh0R6QGzThUHrkQaxSNzRKQHzDoi0gNmnToOXIm0igFHRHrArCMiPWDWqeLAlUij+PQ5ItIDZh0R6QGzTh37uNJ1p/deYLXt91X4xUG35g++o3Wt1qM3l7f+Jq2Z/OS9SMvkrU8BAGfLD0hrseZEac3nuFVac1TYpTVTqLyHaWmYTVoDAHOZ/L1KhI+8V63xsnxbhUne+9Pmq/wrygh5/1xjmfwRjRUW+XdQLuS9Wm1C+XdhsLWhtNbhWfk/kh3zf5LWfOx+0pqhQv4ZDSpH74v8zklrwZfCpTV7gPyYt1Do42rxrd2/K2adZ3C/jjyF+3Xcr/MUnnEl0ijeC0FEesCsIyI9YNap48CVSKOEnQ2/iKj+Y9YRkR4w69Rx4EqkVQ4GHBHpALOOiPSAWaeKA1cijRJ2XlJCRPUfs46I9IBZp44DVyKNEjwyR0Q6wKwjIj1g1qnjwJVIo3gvBBHpAbOOiPSAWadO3ruAiOo2h8O9l4fk5+dj+PDhCAoKQkhICEaPHo3i4mLF99xyyy0wGAwur0ceecRlnhMnTmDgwIHw8/NDZGQknnnmGdhsym1ciKgeYtYRkR4w61TxjCtdd0o9vfTeC0xJXXls+vDhw3HmzBmsW7cOFRUVSE9Px7hx4/D+++8rvm/s2LGYNWuW889+fn/0qrTb7Rg4cCCio6Px/fff48yZMxgxYgS8vb3x8ssve+yzAICxfQNpzWosldYK7b8qLreRsYt8ncXyX0DlZ0ukNe9weX/Pz0vnSmstfZOlNQBoXdBNXvSR900t95X3uS31vyitBZ6Ufw5jUZl8WwAYW4TK12kukNb8LwVJa5f8C6W1AOMNittz2XJJWlPq1dr5ifbS2g/z5FknvBWOP9uUd2x8DRHSmsNb3o/VbpTvaHjLPz5Qyz6uzDoibeF+Xe0w69Rx4EqkUXXhkpIDBw4gKysLO3bsQKdOnQAACxYswG233YZXX30VsbGx0vf6+fkhOjq62tratWvx888/45tvvkFUVBQSEhIwe/ZsPPfcc5gxYwZ8fHw88nmIqO5h1hGRHjDr1PFSYSKNEg6HWy9PyM7ORkhIiDPcACAlJQVGoxHbtm1TfO9//vMfhIeHo127dpgyZQpKS/84m5mdnY327dsjKirKOS01NRVFRUXYv3//tf8gRFRnMeuISA+Ydep4xpVIq9w8Mme1WmG1Wl2mmc1mmM3mWm9Cbm4uIiMjXaZ5eXkhLCwMubm50vfdf//9aNKkCWJjY/Hjjz/iueeew6FDh/DJJ584l/vncAPg/LPScomoHmLWEZEeMOtU8YwrkUYJh3DrlZmZieDgYJdXZmZmtcuePHlylZvsr3wdPHiw1ts+btw4pKamon379hg+fDhWrFiBVatW4ehRfd/fQkRVMeuISA+Ydep4xpVIo9y9F2LKlCnIyMhwmSY7Kjdp0iSMGjVKcXlNmzZFdHQ0zp496zLdZrMhPz9fep9DdZKSkgAAR44cQbNmzRAdHY3t27e7zJOXlwcAbi2XiLSPWUdEesCsU8eBK5FGuXt/gzuXj0RERCAiQv7E0UrJyckoKCjArl27kJiYCAD49ttv4XA4nKFVE3v37gUAxMTEOJf70ksv4ezZs85LVtatW4egoCC0bdu2xsslIu1j1hGRHjDr1PFSYSKNEhV2t16e0KZNG6SlpWHs2LHYvn07tmzZgvHjx2Po0KHOJ8+dOnUKrVu3dh5pO3r0KGbPno1du3YhJycHn3/+OUaMGIFevXqhQ4cOAID+/fujbdu2ePDBB/HDDz9gzZo1eOGFF/D4449f1b0bRKQ9zDoi0gNmnTqecaU6jb3A5OrCY9OB358iN378ePTt2xdGoxF333035s+f76xXVFTg0KFDzqfL+fj44JtvvsG8efNQUlKCRo0a4e6778YLL7zgfI/JZMKXX36JRx99FMnJyfD398fIkSNd+oN5yhHremmtuaWvtBZVoXzE8IAjS1oL942T1hp0Vlhuqbyf5m3mF6Q1r4Jy+TIBOKzyugiV/4IpcuRIa74GeX9cW2N5H1cb5H1jAcAh5I1DfR1h0po9UL7MULs8W4wlyjsLpfKPAh+7vKjUq7XjRPn2ZM+XP+UxqChYvjEAyi3yDCk0n5XWwvLll3XZ/a/9bgWzjqj+4H6dHLNOHQeuRBrlqUehuyssLEyxKXVcXByE+KOpdqNGjbBp0ybV5TZp0gRff/31NdlGItIuZh0R6QGzTh0HrkRaZRfq8xARaR2zjoj0gFmnigNXIo2qK0fmiIg8iVlHRHrArFPHgSuRRtWVeyGIiDyJWUdEesCsU8eBK5FG8cgcEekBs46I9IBZp44DVyKt4r0QRKQHzDoi0gNmnSoOXIk0ikfmiEgPmHVEpAfMOnUcuJJm6b0XmLAx4DyhkW9nac2KYmnNx6LQGBRAS0OatGaskL/PoNBkvNxPXsu3HZDWIs/eIF8hAOsZeW9UR5xCj1ch3x4vg7yH6WVxXr5MKPdNDSwJl9aK/eXLNQijtGYxyfu/CpNBcXssv8hrhsbyn1nhLd8epV6tyU8kSWs7F/ws3xgAfvJWrfAJbCit2QPk34HV+7J8mVBocquAWUekD9yvY9ap4cCVSKN4Ez8R6QGzjoj0gFmnjgNXIo3ikTki0gNmHRHpAbNOHQeuRBrFI3NEpAfMOiLSA2adOg5ciTSKR+aISA+YdUSkB8w6dRy4EmmUw6r80BoiovqAWUdEesCsU8eBK5FG8cgcEekBs46I9IBZp44DVyKN4r0QRKQHzDoi0gNmnToOXKleqm0vMLX31iU8MucZvqX+0pqhQkhrthB5DQC+OT1XWusWM1JaC7nUQL49NpO0Fuov/3f+pb98WwAgtEOktNbN8Ki05msLldYc8jalMBtCpLVyIe8pCwAGu/x79y+T92M1WuU/PwaHvLFufkiu4vYENW8kX26ZwmVgCj/PQUXB0ppSr9ZOE9rK1wdg94JD0pp3qfx7tXpbpTWTway4ztpg1hER9+sI4MCVSLMYcESkB8w6ItIDZp06DlyJNIqXlBCRHjDriEgPmHXqOHAl0igemSMiPWDWEZEeMOvUceBKpFEMOCLSA2YdEekBs06dwuMyiKguE3aHWy9Pyc/Px/DhwxEUFISQkBCMHj0axcXF0vlzcnJgMBiqfX344YfO+aqrf/DBBx77HERUNzHriEgPmHXqeMaVSKPqypG54cOH48yZM1i3bh0qKiqQnp6OcePG4f333692/kaNGuHMmTMu095++2288sorGDBggMv0pUuXIi0tzfnnkJCQa779RFS3MeuISA+Ydeo4cCXSqLpwE/+BAweQlZWFHTt2oFOnTgCABQsW4LbbbsOrr76K2NjYKu8xmUyIjo52mbZq1SoMGTIEAQEBLtNDQkKqzEtE+sKsIyI9YNapMwghlJsP/v/Kyi7XeiVEWqLUD8wTvcAsFt9avW9X0qtuzZ+47elarUfJkiVLMGnSJFy8eNE5zWazwWKx4MMPP8Rdd92luoxdu3ahU6dO2LJlC7p16+acbjAYEBsbC6vViqZNm+KRRx5Beno6DAbDNf8cf2Ytlmed6VSJtOaIVP573HZ5hbRmNgVIay223CytFf98UlprODxJWrOFKffaNCj88nSY5b1j4VDo/ekl/+4MCnetGAwK6wNggo+8VibfHqXPcc7+k7QWaKr6S7um7CiX1nwNEdKaAzZpze+sfH1lkcrf3c0TWklrP2X+Iq2ZSuXbYwv0ltZ8QuU9kpUw6zyD+3WkF9yvq7m6nnU840qkUe5eUmK1WmG1Wl2mmc1mmM3Kgxglubm5iIyMdJnm5eWFsLAw5Obm1mgZ7777Ltq0aeMSbgAwa9Ys3HrrrfDz88PatWvx2GOPobi4GE888UStt5eItIdZR0R6wKxTx4czEWmUsAu3XpmZmQgODnZ5ZWZmVrvsyZMnS2+0r3wdPHjwqj/D5cuX8f7772P06NFValOnTkX37t1x00034bnnnsOzzz6LV1555arXSUTawqwjIj1g1qnjGVcijXJY5ZfrVWfKlCnIyMhwmSY7Kjdp0iSMGjVKcXlNmzZFdHQ0zp51vVbRZrMhPz+/RvcwfPTRRygtLcWIESNU501KSsLs2bNhtVqv6mgiEWkLs46I9IBZp44DVyKNcveSEncuH4mIiEBEhPyeu0rJyckoKCjArl27kJiYCAD49ttv4XA4kJQkv8ey0rvvvos777yzRuvau3cvQkNDuSNHpDPMOiLSA2adOg5ciTSqLjw2vU2bNkhLS8PYsWOxePFiVFRUYPz48Rg6dKjzyXOnTp1C3759sWLFCnTp0sX53iNHjmDz5s34+uuvqyz3iy++QF5eHrp27QqLxYJ169bh5ZdfxtNPX/sHERBR3casIyI9YNap48CVSKPqwmPTAeA///kPxo8fj759+8JoNOLuu+/G/PnznfWKigocOnQIpaWlLu9bsmQJbrjhBvTv37/KMr29vbFw4UI89dRTEEKgefPmeO211zB27FiPfx4iqluYdUSkB8w6dWyHQ3QFrTw2fUvDqW7N3/3U7FqtR2/YDoftcNgOp261w2HWeQb360gvuF9Xf/CMK9EVlELsrw4/JQ5RN47M1TcXjUekNf/G8kFLORRGEAC6igelNVNhhbS2vftn0lpx53PSmn94a2ktr+yAtAYAYZYm8lppY2lNmOS92JQG5w6jfIBpOnJJWgOAc/G/SmvhJvl34H2qVFqL9ZMP6ESx/O8KADb7LZfWevnIjyw7vOXfXaFZ/m/LJ7ChtOZdqnxcWmlw2n5KS2ntx1cOS2uGimufS8w6Iroa3K+rPzhwJdIoR80uliAi0jRmHRHpAbNOHQeuRBpl55E5ItIBZh0R6QGzTh0HrkQaxUtKiEgPmHVEpAfMOnUcuBJpFC8pISI9YNYRkR4w69Rx4EqkUTwyR0R6wKwjIj1g1qnjwJVIoxhwRKQHzDoi0gNmnToOXIk0ipeUEJEeMOuISA+Ydeo4cCVyQ13qBcYjc54RbJT/PZaIM9KaEHbF5QqjvE9ncXS5tBYnEqU1b4O8N6oPgqS1Zsbe0hoAFOCktHbJ/7y0FpIbLK05fOW/bgx+JmmttJny9xqMOGnNaiyRvzFG/t0p/V0ZgrwVt6d34WhpzRYg/w7sRpu0FpYfLX9fgHxbrd5WaQ0A/M4ZpTWlXq0dnmkhrW2Zv0Vaa4BAxe2RYdYRkadwv05bOHAl0qgKh/IOPRFRfcCsIyI9YNap48CVSKN4ZI6I9IBZR0R6wKxTx4ErkUYx4IhID5h1RKQHzDp1HLgSaRRv4iciPWDWEZEeMOvUceBKpFE8MkdEesCsIyI9YNap48CVSKPsDDgi0gFmHRHpAbNOHQeuRBrFI3NEpAfMOiLSA2adOg5cia6R2vYCO/jO6Vqtj/dCeIZ3gfxx9JaQBtKal8GivGCD/BeSSeG9oaK5fJHl8n8DBrt8fbZf8qU1AAhtFiutXfSX/3u1hfrIaxb59lSIImlNQLk9gLfDV1rzKqyQ1opDC6W1gIt+0lpFmPwzAgAM8r6qQqHmfUm+SLu//Fe11fuytGYymOULBWALVOhXWyH/+1Lq1dr9ie7SGrOOiLSE+3V1DweuRBrFI3NEpAfMOiLSA2adOg5ciTSKAUdEesCsIyI9YNapM17vDSCi2nEI4dbLU1566SV069YNfn5+CAkJqdF7hBCYNm0aYmJi4Ovri5SUFBw+fNhlnvz8fAwfPhxBQUEICQnB6NGjUVxc7IFPQER1GbOOiPSAWaeOA1cijbILh1svTykvL8e9996LRx99tMbvmTt3LubPn4/Fixdj27Zt8Pf3R2pqKsrKypzzDB8+HPv378e6devw5ZdfYvPmzRg3bpwnPgIR1WHMOiLSA2adOl4qTKRRdeWSkpkzZwIAli1bVqP5hRCYN28eXnjhBQwaNAgAsGLFCkRFReHTTz/F0KFDceDAAWRlZWHHjh3o1KkTAGDBggW47bbb8OqrryI2Vv7wICKqX5h1RKQHzDp1PONKpFF15ZISdx07dgy5ublISUlxTgsODkZSUhKys7MBANnZ2QgJCXGGGwCkpKTAaDRi27Ztf/k2E9H1w6wjIj1g1qnjGVcijXL3yJzVaoXVanWZZjabYTYrt8y41nJzcwEAUVFRLtOjoqKctdzcXERGRrrUvby8EBYW5pyHiPSBWUdEesCsU1fjgavFIu+VR0TKatvTS8lTjo/dmn/GjBnOyz8qTZ8+HTNmzKgy7+TJkzFnzhzF5R04cACtW7d2axu0wCs6UFoLuJoFK7R5VekMes15d639J4lEK3lR3v5U5TPKv/Or4i//vRWMIPn7FH7dmdTWGSv/i1Z8r69KH2AJH6UvXU3tVokGCn9fzDrt4H4dUe0x664PnnEl0okpU6YgIyPDZZrsqNykSZMwatQoxeU1bdq0VtsRHR0NAMjLy0NMTIxzel5eHhISEpzznD171uV9NpsN+fn5zvcTEVWHWUdEeqDHrOPAlUgn3Ll8JCIiAhERER7Zjvj4eERHR2P9+vXOQCsqKsK2bducT7BLTk5GQUEBdu3ahcTERADAt99+C4fDgaSkJI9sFxHVD8w6ItIDPWYdH85ERFflxIkT2Lt3L06cOAG73Y69e/di7969Lr25WrdujVWrVgEADAYDJk6ciBdffBGff/45fvrpJ4wYMQKxsbEYPHgwAKBNmzZIS0vD2LFjsX37dmzZsgXjx4/H0KFD+ZRNIroumHVEpAd1OusEEdFVGDlypABQ5bVhwwbnPADE0qVLnX92OBxi6tSpIioqSpjNZtG3b19x6NAhl+VeuHBBDBs2TAQEBIigoCCRnp4uLl269Bd9KiIiV8w6ItKDupx1hv9/5URERERERER1Ei8VJiIiIiIiojqNA1ciIiIiIiKq0zhwJSIiIiIiojqNA1ciIiIiIiKq0zhwJSIiIiIiojqNA1ciIiIiIiKq0zhwJSIiIiIiojqNA1ciIiIiIiKq0zhwJSIiIiIiojqNA1ciIiIiIiKq0zhwJSIiIiIiojqNA1ciIiIiIiKq0/4/wD/Jt4HtGD8AAAAASUVORK5CYII=", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "fig, axes = plt.subplots(1, 3, figsize=(9.6, 3.0))\n", - "cmap = cm.PiYG\n", - "\n", - "sns.heatmap(prec, vmin=-1, vmax=1, cmap=cmap, ax=axes[0])\n", - "axes[0].set_xticks([])\n", - "axes[0].set_yticks([])\n", - "axes[0].set_title('True precision')\n", - "\n", - "sns.heatmap(prec_sample, vmin=-1, vmax=1, cmap=cmap, ax=axes[1])\n", - "axes[1].set_xticks([])\n", - "axes[1].set_yticks([])\n", - "axes[1].set_title('Sample precision')\n", - "\n", - "sns.heatmap(prec_scope, vmin=-1, vmax=1, cmap=cmap, ax=axes[2])\n", - "axes[2].set_xticks([])\n", - "axes[2].set_yticks([])\n", - "axes[2].set_title('Estimated precision matrix')\n", - "\n", - "plt.tight_layout()\n", - "plt.show()" - ] - }, { "cell_type": "markdown", "id": "0172cdbf", diff --git a/docs/source/userguide/quickstart.py b/docs/source/userguide/quickstart.py new file mode 100644 index 0000000..f66e164 --- /dev/null +++ b/docs/source/userguide/quickstart.py @@ -0,0 +1,22 @@ +from sklearn.datasets import make_regression +## generate data +n, p, s = 100, 10, 3 +X, y, true_coefs = make_regression(n_samples=n, n_features=p, n_informative=s, coef=True, random_state=0) +print("X shape:", X.shape) +print("Y shape:", y.shape) +import jax.numpy as jnp +def objective_function(coefs): + return jnp.linalg.norm(y - X @ coefs) +from scope import ScopeSolver +scope_solver = ScopeSolver( + dimensionality=p, ## there are p parameters + sparsity=s, ## we want to select s variables +) +scope_solver.solve(objective_function) +import numpy as np +est_support_set = scope_solver.get_support() +print("Estimated effective predictors:", est_support_set) +print("True effective predictors:", np.nonzero(true_coefs)[0]) +est_coefs = scope_solver.get_estimated_params() +print("Estimated coefficient:", np.around(est_coefs, 2)) +print("True coefficient:", np.around(true_coefs, 2)) \ No newline at end of file diff --git a/docs/source/userguide/quickstart.rst b/docs/source/userguide/quickstart.rst index c90c29d..f04ae84 100644 --- a/docs/source/userguide/quickstart.rst +++ b/docs/source/userguide/quickstart.rst @@ -1,86 +1,160 @@ :parenttoc: True -A Quick Example +Quick Example for Beginner ====================== Introduction --------------------------------------- +----------------- -Here, we will take linear regression as an example, introduce the basic usage of :ref:`scope `. to solve a variables selection problem. +.. Here, we will take linear regression as an start-up example to introduce the basic usage of :ref:`scope `. -Suppose we collect :math:`n` independent observations for a response variable and :math:`p` explanatory variables, say :math:`y \in R^n` and :math:`X \in R^{n\times p}`. Let :math:`\epsilon_1, \ldots, \epsilon_n` be i.i.d zero-mean random noises and :math:`\epsilon = (\epsilon_1, \ldots, \epsilon_n)`, the linear model has a form: - -.. math:: - y=X \beta^{*} +\epsilon. +Presume we collect a dataset with :math:`n=100` independent observations on a interesting response variable and :math:`p=10` predictive variables: .. code-block:: python from sklearn.datasets import make_regression ## generate data - n, p, k= 10, 5, 3 - x, y, true_params = make_regression(n_samples=n, n_features=p, n_informative=k, coef=True) + n, p, s = 100, 10, 3 + X, y, true_coefs = make_regression(n_samples=n, n_features=p, n_informative=s, coef=True, random_state=0) + print("X shape:", X.shape) + >>> X shape: (100, 10) + print("Y shape:", y.shape) + >>> Y shape: (100,) + +We are interested in explaining the variation of the response variable with :math:`p` predictive variable, but in practice, the underlying true parameters ``true_coefs`` is unknown. Our only knowledge is that only :math:`s=3` predictive variables would influence the response variable, and the relationship between the predictive variables and the response variable has a form of linear model: + +.. math:: + y=X \theta^{*} +\epsilon. + +where + +- :math:`y \in R^n` and :math:`X \in R^{n\times p}` records the observations on the response variable and the predictive variables, respectively; + +- :math:`\epsilon = (\epsilon_1, \ldots, \epsilon_n)` is a vector that contact :math:`n` i.i.d zero-mean random noises; + +- :math:`\theta^{*}` is an unknown coefficient vector to be estimated. + +With ``skscope``, users can well estimate the :math:`\theta^{*}` and find out the predictive variables that is influential of the response variable, so called variables selection in literature. +To see this, a present a step-by-step procedure below. +A Step-by-step procedure +------------------------------- -write the objective function with ``JAX`` package ------------------------------------------------------- +Sparsity-constrained optimization perspective +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Formalize the problem as a sparse convex optimization problem: +We first formulate the variables selection problem as an optimization problem: .. math:: - arg\min_{\beta \in R^p}L(\beta) := ||y-X \beta||^{2} s.t. ||\beta||_0 \leq s, + \arg\min_{\theta \in R^p} ||y-X \theta||^{2} s.t. ||\theta||_0 \leq s, -where objective function :math:`L(\beta)` is our optimization objective and the parameter, denoted as ``params``, is a vector which SCOPE needs to optimize and is imposed by sparse-constrain. +where -The length of ``params`` is actually the dimension of the optimization problem so denoted as ``dimensionality``. The number of nonzero parameters is the sparsity level and denote as ``sparsity``. +- :math:`\theta` is a coefficient vector that needs to be optimized -From the perspective of variable selection, each parameter corresponds to a variable, and the nonzero parameters correspond to the selected variables. +- :math:`||y-X \theta||^{2}` is the objective function that quantitatively measures the quality of the coefficient vector :math:`\theta`; + +- :math:`||\theta||_0` counts the number of non-zero element in :math:`\theta`. And thus, the sparsity constraint :math:`||\theta||_0 \leq s` reflects the prior knowledge that the number of effective predictive variable is less than :math:`s`. + +The intuitive explanation for the optimization problem is: finding a small subset of predictors, so that the resulting linear model is expected to have the most desirable prediction accuracy. Therefore, it fully leverages information the observed data and our prior knowledge. -In the above example, :math:`\beta` is the parameters, so ``dimensionality`` is :math:`p` and ``sparsity`` is :math:`s`. -``JAX`` is Autograd version of ``Numpy`` for high-performance machine learning research. It is a Python library that provides a NumPy-compatible multidimensional array API and automatic differentiation. As the usage of ``JAX`` is similar to ``Numpy``, we will not introduce it here. For more information, please refer to `JAX `__ . +Solve SCO via ``skscope`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -In the above example, we can define the objective function :math:`L(\beta)` as ``custom_objective(params)``: +The above optimization problem is a sparsity-constrained optimization (SCO) problem. Since ``skscope`` is designed for general sparsity-constrained optimization, it can definitively solve the SCO problem. More importantly, ``skscope`` is extremely easy to use --- it can solve the SCO problem so long as the objective function :math:`||y-X \theta||^{2}` is properly implemented. Next, for the above example, the objective function is programmed into a python function ``objective_function(coefs)``: .. code-block:: python import jax.numpy as jnp + def objective_function(coefs): + return jnp.linalg.norm(y - X @ coefs) - def custom_objective(params): - return jnp.sum( - jnp.square(y - x @ params) - ) - +Notice that, the first line imports the only required ``jax`` library [*]_, which is a Python library that provides a NumPy-compatible multidimensional array API and automatic differentiation. As the usage of ``jax`` is very similar to ``numpy``, we can simply understand ``jax.numpy`` is equivalent to ``numpy`` [*]_. -initialize ``ScopeSolver`` and solve the problem ------------------------------------------------------------- +Moreover, we can see that the programmed ``objective_function`` is exactly match to the mathematics function :math:`||y-X \theta||^{2}`. -Those concepts are introduced in the previous section. +Then, we will input ``objective_function`` into one solver in ``skscope`` to get the practice solution of the SCO problem. -- ``dimensionality`` is the number of parameters and must be offered. -- ``sparsity`` is the sparsity level, int or list of int. If it is an int, the solver will use the given sparsity level, otherwise, the solver will search the best sparsity level from the given list which will be introduced in the next section. +.. The length of ``coefs`` is actually the dimension of the optimization problem so denoted as ``dimensionality``. The number of nonzero parameters is the sparsity level and denote as ``sparsity``. -.. code-block:: python +.. From the perspective of variable selection, each parameter corresponds to a variable, and the nonzero parameters correspond to the selected variables. - solver = ScopeSolver( - dimensionality=p, ## there are p parameters - sparsity=k ## we want to select k variables +Here, we import the ``ScopeSolver`` [*]_ and properly configure it. + +.. code-block:: python + from skscope import ScopeSolver + scope_solver = ScopeSolver( + dimensionality=p, ## there are p parameters + sparsity=s, ## we want to select s variables ) -``solve`` is the main function of SCOPE, it takes the objective function as optimization objective and commands the algorithm to begin the optimization process. +In our configuration, we set: + +- ``dimensionality``, i.e., is the number of parameters and must be offered. + +and + +- ``sparsity``, i.e., the desired sparsity level. + +In the above example, :math:`\beta` is the parameters, so ``dimensionality`` is :math:`p` and ``sparsity`` is :math:`s`. + +Those concepts are introduced in the previous section. + +With ``scope_solver`` and ``objective_function``, we use one-line command to solve the SCO: .. code-block:: python - solver.solve(custom_objective) + scope_solver.solve(objective_function) + + +``solve`` is the main method of solver in ``skscope``, it takes the objective function as optimization objective and commands the algorithm to conduct the optimization process. + +Get solutions +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Since the solvers in ``skscope`` are all coupled with necessary functions to extract, we can also use one line to get the result. For example, we may interested in obtaining the effective variables, which can be extract by using +the ``get_support`` method. The code is below. + +.. code-block:: python + import numpy as np + est_support_set = scope_solver.get_support() + print("Estimated effective predictors:", est_support_set) + >>> Estimated effective predictors: [3 5 6] + print("True effective predictors:", np.nonzero(true_coefs)[0]) + >>> True effective predictors: [3 5 6] +We can see that the estimated effective predictive variables are exactly the true one, which reflects the accuracy of the solver in ``skscope``. -get results ------------------------------------------------------- +Else, we may interested in the regression coefficient: -- ``get_estimated_params`` returns the optimized parameters. -- ``get_support`` returns the index of selected variables (nonzero parameters). +- ``get_estimated_params`` returns the optimized coefficient. .. code-block:: python - beta = solver.get_estimated_params() - support_set = solver.get_support() + est_coefs = scope_solver.get_estimated_params() + print("Estimated coefficient:", np.around(est_coefs, 2)) + >>> Estimated coefficient: [ 0. 0. 0. 82.19 0. 88.31 70.05 0. 0. 0. ] + print("True coefficient:", np.around(true_coefs, 2)) + >>> True coefficient: [ 0. 0. 0. 82.19 0. 88.31 70.05 0. 0. 0. ] + +For the output, we see that the estimated coefficient approaches to the underlying true coefficient. + +Further reading +--------------------------- + +- `JAX library `__ + +- A bunch of `machine learning methods `__ implemented on the ``skscope`` + +- More `advanced features <../feature/index.html>`__ implemented in ``skscope`` + +Footnotes +--------------------------- + +.. [*] For simplicity, we just introduce the purpose of ``JAX`` library. For more information, please refer to `JAX `__. + +.. [*] If you know nothing about ``numpy``, we can turn to `this introduction `__. + +.. [*] We skip the algorithmic detail about ``scopeSolver``. Please turn the paper "sparsity-constrained optimization via splicing iteration" if your are interested in. \ No newline at end of file