{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["# Introduction to a numpy API for ONNX: CustomClassifier\n", "\n", "This notebook shows how to write python classifier using similar functions as numpy offers and get a class which can be inserted into a pipeline and still be converted into ONNX."]}, {"cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [{"data": {"text/html": ["
run previous cell, wait for 2 seconds
\n", ""], "text/plain": [""]}, "execution_count": 2, "metadata": {}, "output_type": "execute_result"}], "source": ["from jyquickhelper import add_notebook_menu\n", "add_notebook_menu()"]}, {"cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": ["%load_ext mlprodict"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## A custom binary classifier\n", "\n", "Let's imagine a classifier not that simple about simple but not that complex about predictions. It does the following:\n", "* compute the barycenters of both classes,\n", "* determine an hyperplan containing the two barycenters of the clusters,\n", "* train a logistic regression on both sides.\n", "\n", "Some data first..."]}, {"cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [{"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnQ0lEQVR4nO3df3Qc5bkf8O9raaVdW1pD8SY3tbHWxGmKw6EIy4SU5OYQG0LSJIAPAfRHTogXcE6vayBpLwQX/jE9TQ65QLjcRuZW1Nw2EpSQQnvbi8E4PSm94coydluwQ0pAAlGSXRssbFmyV+jpH6OVZndndmdnZ+ad2fl+ztkje7U/3v2h95n3fZ/3GSUiICKi+FmiuwFERKQHAwARUUwxABARxRQDABFRTDEAEBHFVLvuBjRixYoVks1mdTeDiChSDhw4cFREMpXXRyoAZLNZjI6O6m4GEVGkKKXGra7nFBARUUwxABARxRQDABFRTEVqDcBKsVjExMQEZmZmdDfFF8lkEqtWrUIikdDdFCJqMZEPABMTE+ju7kY2m4VSSndzPCUiOHbsGCYmJrBmzRrdzSGiFhP5KaCZmRmcc845Ldf5A4BSCuecc07Ljm6Ioq5QAPbvN35GUeQDAICW7PxLWvm1EUXZ8DDQ0wNccYXxc3hYd4sa1xIBgIgoSIUCkMsB09PA5KTxM5eL3kiAAUCTt956C5/97Gexdu1a3HDDDThz5ozuJhGRQ2NjQEdH+XWJhHF9lDAAaHLnnXfijjvuwBtvvIGzzz4bg4ODuptERA5ls0DlMVuxaFwfJbEMAF4u3Nx777146KGHFv6/Y8cO/OQnP6l5HxHBvn37cN111wEAvv3tb+OZZ55pvjFEFIhMBhgcBFIpIJ02fg4OGtdHSeTTQBs1PGzM1XV0GBF8cBDo73f/eFu2bMHmzZtx++23Y25uDk888QT27duHiy66yPL2Q0ND+NjHPoazzjoL7e3G279q1Sq8++677htBRIHr7wc2bTKmfbLZ6HX+QMwCgHnhZnrauC6XMz5Etx9eNpvFOeecg4MHD+IPf/gDent70dPTg0OHDtne5+jRo+6ejIhCJZOJZsdfEqsAUFq4KXX+wOLCTTMf4s0334zdu3fj97//PbZs2YITJ07gC1/4guVth4aGcP755+P48eOYnZ1Fe3s7JiYmsHLlSvcNICJyIVYBwK+Fm2uvvRb33nsvisUihoaG0NbWVnMEAACXX345fv7zn+PGG2/E448/jquvvrq5RhARNShWi8B+Ldx0dHTg8ssvx/XXX4+2tjZH9/nRj36EBx54AGvXrsWxY8eQy+WaawQRUYNiNQIA/Fm4mZubw8svv4ynnnrK8X3OO+88jIyMNP/kREQuxWoEUJLJABs2eNP5Hz58GGvXrsXGjRvxqU99qvkHJCIKSOxGAF5bt24d3nzzTd3NICJqWCxHAERExABARBRbDABERDHFAEBEFFMMAJo88sgjWLt2LZRSLA1BRFowAGhy2WWXYe/evejp6dHdFCKKqXgGAA/rQbspBw0Avb29yEateDgRtZT47QPwuB60m3LQ69atc/18REReiVcA8KEetJty0EREYRCvAOBTPehGy0FzBEBEYRCvAOBTPWg35aCJiHSL1yKwT/Wg3ZSDfvjhh7Fq1SpMTEzgwgsvxM0339xUG4iIGhWvEQDgSz1oN+Wgt2/fju3btzf93EREbsVrBFDiYT1oloOmuPIwm5o00TYCUEqdC+CvAHwcgAB4VETqJ9CHDMtBUxx5nE1NmugcAcwC+L6IrANwKYA/UUq5So8REU8bFiat/NoomszZ1JOTxs9cjiOBKNIWAETkPRF5Zf7fJwAcAbCy0cdJJpM4duxYS3aUIoJjx44hmUzqbgrRglI2tVkpm5qiJRSLwEqpLIBeAH9n8btbAdwKAKtXr666bymTptCihx/JZBKrVq3S3QyiBT5lU5MG2gOAUqoLwNMAbheRDyt/LyKPAngUAPr6+qoO8xOJBNasWeN7O4nIUMqmzuWMI/9i0ZNsatJAawBQSiVgdP4/E5Ff6GwLETnnQzY1aaAzC0gBGARwREQe0NUOInInk2HHH3U6s4AuA/AtAF9SSh2av3xVY3uIiGJF2whARF4CoHQ9PxFR3MVzJ7BXuBWSiCKMAcCt4WGgpwe44grj5/Cw7hYRETWEAcANboUkohbAAOAGt0ISUQtgAHCDWyGJqnBJLHoYANzw6cQyRFHldkmMQUMvFaUian19fTI6Oqq7GYsKBW6FpNgrFIxO33yq7VQKGB+v/WfBktLBUUodEJG+yus5AmiGhyeWIYoqN0tizKMIBwYAImqKmyUx5lGEAwMAETXFzZIY8yjCgQGgFXFlLfKi9hH29xtz/nv3Gj/rzeUzjyIcGAB08POvmzuUIy+qH2GjS2KNBg3yHrOAguZn6oPbdAwKDX6E5AdmAYWB36kPXFmLvKA/wqhNNZG3GACC5PdfN1fWIi/Ij9CLqSYGkGhjAAiS33/dXFmLvKA+Qi8Go1Fdq6BFXAMIWmkNwHw2ba9Xv7hDOfL8/gj37zc67snJxevSaWNBdsMGZ+3jWkV02K0BaD0pfCwFcTZtnqw18vz+CJsdjJZmM80BoDSbya9edHAKSIeolJDgBG/LanaqictNrYEBgKxxgrflNZOHz+Wm1sA1AKrGCV5yiMtN0cA1AHLOLi2VE7xUgctN0cYpIKrW1VV+9A8Y/+/q8u85ud4QCfyYWgsDAFU7edKY8jFLJo3r/cD1hkjgx9R6uAZA1YJcA+B6QyTwY4o21gIi54JM8WD9okhMq/Bjak0MAGQtqFq9MU8oj8q0SlQ+pigE0zBhACB7QWxYi3FCeZTOixuFjykqwTRMuAYQBXFIto7Da6zQbD0eHcL6MXGNojauAURVXA5rolIew0NRmVYxC+vHxDUKdxgAwixKcwTUsChMq0RFFINpGDAAhBkPa1qejvPituJCKYOpOywFEWY8rNHKyXy3mznxyvsEWU7Bz1NS6xZEpfVWo3UEoJR6TCmVV0q9qrMdocXDGm2cLL24WZ7RuaQThxnFsK5RhJXWLCCl1B8DOAngr0Tkgnq3ZxZQNp7f7IBfv5OMEjdZJ7ozVaKYdUTeCGUWkIj8CsD7OtsQCXE+rNFwyOxk6cXN8ozuJR3OKFKl0C8CK6VuVUqNKqVGC600VqX6NM1ZOOko7W7T1WW/wKq7A/ZyRrEVF5LjKPQBQEQeFZE+EenLxPEIOM40HTI76SitbpPLAevX2w9WwrCk40XWUVy2psSB9p3ASqksgL/mGkAE+T03r3nSvJEsoK4uo/N30tQoL+noXscgd0K5BkAhVm+MH8RhoOZDZrulF/NbU7rNyZPOBytRXtLRvY5B3tKdBjoM4NcAPq2UmlBK5XS2h+bV69yDnJvXsVPKRqEA3Hef9Vuje34/KHF5nXGhOwuoX0Q+ISIJEVklIoM62xN5XqzMOencgz4MDMEh8/AwsHo1cM891m9NGOb3gxCX1xkX3AncKrza4lnq3M2TvKXOvfRXHrPDwFJMnJmp/p35rYnLTtS4vM444BpAK/BySsZJ5x6zw0CrAU+J1VsThvl9v9M0w/I6qTkMAK3AyykZp517iObmnarbKdrcwComAuGNe0zTJMdEJDKX9evXC1nI50VSKRFg8ZJKGdc385gjI809RogMDRlvyfLlxs+hocZuUPp1Oi2STIrs3BnOt8aPr4InjWqh71IUARgViz5Ve6feyIUBoAZzD2XZw8VX3U7RYa8ZhX5sZMSIYeaXkk4b19fi22urG3kpCHYBgFNArSKCUzJBqTtD5nAKLQrz3m7W532bMopD+dGIYwAIs0ZX8qLQQ2lQt1NsoaymRtfnfe2juWss9BgAwooreZ6p2yl6kNUUpuJojQwGfe2jWyiwtirttYAaEZtaQCy44ou6NXhcFunx/SxbPhYP8v2rVnpzEgmj82+lU5D5xYfPm7WAooRDZ1/UnSFzMYXm+zS3aSQoPT343X3Dno4yfN/SwbWpxliN/H0cXnIEEEYcAQSnyaMtX8+yZfE9OIUU/mFyHD96LONpXxrlCqUtw+rvPpEA2tubHl5yBBAlMdtpq40H6yy+TnNbjASLSODjM2OeJ9M0OvgJ05pHy7Aa+ReLvmZRMQCEFYfO/vJo7iaTMe5mlst5FKstoksCRYwhq3VGkPkJPrHbcm7m8QfPABBmTOv0j0frLIWCMTgzGxz06CBtfiQ4m0hhEmmcQgpbMIijyDgeZXh9pM7Ufh9VjvyTSesRgYdZVAwAFE8ezd34vV5f2NSP89rGsQl70YNxPAljJPjgg/WPC/w4Umd+gs/MI/+33wZ27/Z1KpjloClWFhc7M8gMDlanKDb4x+V3qvvYGPBhZwbvzCy2q7sbuPji2vczH6mX1hRzOaOMczP9B1P7A5DJLH5IPtfe5giAYqPqiBjNr7P4vV5v1eHOztbvcP06Umd+ggY+TgUzDZRiwe/MWj/TKN3spYry6yXv2aWBcgqIWopdx+TkRGfNMI/aveZmFqB0pN7kDFfNx2fHH30MANQyzCUZ0qcL+MsdY/jy1iyQydjOXXd1GVkyYT+SddPh8tSNVE/NNQClVFop9UmL6y/0r0lEzpXSHI8cWVz0/MrkMH4z04NL7zHKJ2B42HLuOpcD1q9vPksmzJuimElMtdiuASilrgfwEIA8gASAm0Rk//zvXhGROnkI3uMaAJmZj/hnZoAlS4Bl0wWMowdLYT35XZoi6uoyOv9m58h9LwRH5AE3pSDuBrBeRC4C8B0A/14pdW3p8bxvIsWCR4fLhQJw55YCPjO9H4nJAk6fNjrzLMZQrJzZNKW/lI6IT55sPkuGm6Io6moFgDYReQ8ARGQEwOUA/qVSajuA6KQOUXh4uDPpw13GNM8LuALj6MENGEYqBVzS9grSOFF+Y4tEdS/y2cOyKSrMU1C2ItnoFmR1nsj5aaG/BfDJiuu6AbwI4LTd/fy88JzAEebl2crzeZmreKwppOQfdR6Wj5IVzwGI3H+/5Qlvmz2NchhOwB7YKXe9PGkwzxMcODR6UngAFwJYa3F9AsC9dvfz88IAEGFuz1bu8LGOIy1/u3V39XN0dhoXm86m2X6t2SDilFU7AwtAXnbYYYiaMWQXAGpNAT0DYLNSqq10hVLq4wD+HYBveD0SoRbnZQ0Bi8dKp4r43G2XVD/H6dPGxWaSvtksmSCKttrNnB08aCx8m3k+BeX1QofbeTNOGfmiVgBYD+A8AIeUUl9SSt0GYATArwFcEkTjqIU0UUOg6m/f4rHU4CBw/vnl13d0AJ2d5Q/mwyS9n6mWdv3vrl3ANdcAU1Plt/e8SqjXCx1uDgRYf9o/VsMC8wXAbQDmAEwAWFXv9n5eOAXkIy/neD18npqzD3aPlc+L7NwpkkxWrwckkyJ79kRmysFq5qy725jVsnppTmZnGprR8WPKppF5s0aeP6jvcATBxRrAWQB2ATgE4EoYewL+D4Av2d3H7wsDgE9Cuijnuu+xuiMg0t4u0tERutdZi9VL6ew0goD5umXLjLjm5vHqvqd+LHQ47aydrh2F9DscFm4CwJsA/jmAdtN1F8HIDhq2u5+fFwYAH4R4Uc71uvHISPVr6uioPmwOyesssesTK/vfgYH6H5ndY7l+T/0+uq41mnPyYkP6HQ4LuwBQaw3gj0XkxyIya5ouOiQi/xjAPm8moEi7sCSzW3C9btzVVb7FFzAeqN1+g5hutaa5Kxeat26tvZxS67Fcv6duFzoqFxusFh9qNdjJ2lGIv8OhZxUVwnrhCMAHmo6enB5Qupp9sBoBlNJBw3KUaHoD8nmRc5N56cOIrEDecdPcpoYGlbpaNS2zbVv1NI3T71+tLwxHAHWh0SmgIC4ArgLwOoA3ANxV7/YMAD4JrEcofzqn07UNzz7YdQiluZOAXqetijfgvW9ukymk5AMslymk5AYMeblFwvKxyt5TP6Z37NZhKj+TPXu82R8S8Hc4akIXAAC0AfgdjFTTDgD/C8C6WvdhAPBRQBkUgW9equwQNGeKFA7nZbaz/A2Yq+gYp5CSc5N5t5ukG3t//Vo8tYpElZd02ggAHu4QZxaQNbsAoPOUkJcAeENE3hSRMwCeAHC1xvbEW0B1gwObrrXboeXh6zxyBHj8ceOnE8PDwLW9Yzh5uvwNqKysWEQCf7ljzFUTG9pu4Wc1O6vFhkrFItDb6905Jln7unFWUSGIC4DrAPxb0/+/BeARi9vdCmAUwOjq1at9io8UlDBN1zZzwLhtW/lr2Lat/nOlUiIrkJcp1J4amfPgDXH02rwsz2GlchRWWgOwmqbh0buvEMIpIEcBwHzhFFBrCMN0beXMx8CA8/7n8GHrvvull+zvY+5rb8CQTCElx5E2poNqdYx+CiIaV3bs7Oi1CGMA+ByAPab//wDAD2rdhwGgdejsB+zWJ7u7a/S/8w0uHM7LD79fnrFj3mpg13dXPucK5OXzncbjlW7w/p4ReWWPu7l/1xoZykStMw97+wIUxgDQDmOz2RosLgJ/ptZ9GADIC/XWJ6sOgueHC6dTy2UGHTKDRFnGjtMD6Fojn9Lvli0LcBDQyAjASUpnmHBncJnQBQCjTfgqgN/CyAbaUe/2DADkhXoZimXT4HVuPIVU2Uig3hS6Xe7+30+Ujyo6Onw+cM3nRR5+WGTp0hov3nRbJymdYTnSDtNCU0jYBQCdWUAQkf8mIv9ARD4pIv9KZ1uoNVltPDVnynR1Vd+nbGesVdqS+bZIIIsx6/tatMEqUeX//dkw/m+x/OxmZ84Y5Z59MTwMrFwJbN8OnDpV8YIsXkCd9wBAuHbecmewY1oDAFFT6tQ0dlJeYd8+YGCgRhZinXTGZYki/pDM2mYw1q1kXCjggodyWIppnIVJLMU0HkMOK1Aou41ntfCPHAFuusno6Cslk9YpmE5TOt2c28EPXp57otVZDQvCeuEUEC2oM8fb6CxAzfXChTWA9MIawHGkpdhhPG8zdcxkZEQ+6l5edqPjSMulbSPG7ZpJV7J6HVZ1pM27cmvd12lKZxiEIdUsRGAzBaSM30VDX1+fjI6O6m4G1VEoGKPtbNanPTmFgnE4bS74lkoZh/TzT7h/v3HUPTm5eJN02tgXtmGDizbO3+FoVxbvvANkMYaze7PG72weqF4b7F7LKaTw3MA4Nm9G9esEgO5uYHYWuPtuozJcjTaUtd/qsUoq3r9a78HC8/j+QTcp7O0LkFLqgIj0Vf3CKiqE9cIRQPgFknzhYANTraNvz9ro1Shk/nE+6jb2BUwOzBdJ2727uvB/5cXpOQ5qpT7Vyl+lloAwZgE1emEACLfAki8cPpHVLIBnbWyiDbaPV5raKd2pXuffSDaOVXs7OkSefNJZ5c0w5tSHsU0hZRcAuAhMngks+cJhwRurckCetdHhAzk5aXyhYGT8fPABgKNHF+vznDixeKNly5y1y+7FWL1nu3cD119ff9U6jOfkDWObIohrAOQZB1Pz3j9hg3O8R48UcG3vGH5zOoujyCy08cAB4OTJBtcEPHixw8PAc98exk+LORTRgWXtM2hfIuVZLMmkkao0MwPccYf9PL6TNtR6z6xeUzIJKBXgh+pA4F+06LNbA+AIgDx1991Gn9FsYUdHGq3+ODyMFet78OISI9/+puQwUinjgHv9+toHk1WZmA2V3bRWKAB3bingp0UjDXQ5JtE+expSmcI4MwNccomx4Ds+DuzcufgmJxLGSMRpG2q9Z1ajmrY2YElFN6E7p555/t6xmhcK64VrAOFVuR66c2fIpmYt5sBnO1Py+kt5x2fQslpnLRzOy2u7TTV9GvDKnrxsTe6W4yif6688P4Akk7XP6OLVXLjVOkEyGb5dtdzp2zBwEZj8Eta/x7J+0SZz6LXdI3Jed3kZBnNCUSPZRHVT9CsWeudSKZlEd1WHXxUAgnwzrVatw5hTH8Y2hRgDAPlmZKQ6YcXLsvJuVHbOTw9Y9+Qf3j9QdTpGc39rl3FqdSIrwHgfzk3m5bmdI9bDiO5uI/smkajq9CfRbWwu073JillALYcBgHwzMFDdEeocAdgdtU8OVBw1ls4RbLrhFFJGsKjzWFansgVEbsGATKPTOLKvlXtacSku7ZLJh3dHp9QyRYpdAOAiMDWlUDASUyo9+KC+hAy7NcLXL67Iybz44qobJrsT2Lzm4MKKr9Va7394sIAN2I/06fLaPLdgF3bhu0jiNNI4AVU6xeLBg0B7e802t8tHSN/41cU3jac3pADU/lYS1VHqbM0ZeV1dRt9ai5+79GvWAstkyp+w4oZqZhpz37gG6OjAktkzwOAg+jdtwlXPjGEMWXzyrb1I35EDlizBmx/NYUtiEM8k+9F5ooBHsK3q/L4LHb9dMbVly4C5OZ/TpYhsWA0LwnrhFFD4uFkADqJchOM1QtMNi4mkzKCj/MUkEouNTSaN0gvmuftEh7yyJy//Zfue6sVbwCi+ls9bz5Mlk8ZcEqd5yGfgFBD5odF0+EJhcaPr5KTxM5fzptKxWd0duKXE/k2bgPFxfPDUXlyjnsU0UmU3k2JxsbEzM0YRNhNVPINeHMTX/olNQ+64w3gztm41NnN1dhrF3FIp4LHHgCuvDO7I38uy0tQSGACoaU7KHZQEuYfHdhrdVEZAenrwu117sR8bcKSzFx2oU/feTm8vVCJRfl17O/C97y3+f+tW4J13gBdfNLYer10bXGfM0glkxWpYENYLp4Cir96Uke/JLxYNmEJKVnbkpa1N5CFsk7n5tMw5QOaWtNXM3pFEYrGxQ0PGtM6yZcZPu3mnoM9XG9aNGhQYcAqIwqDWlFG9g1TzDIbr2QyLIUgRCXzizBjO/qiAWzAIBSxe5j4ybpRKGZdt24wyDMuWGT8ff3xxiNHfD7z9NvDLXxo/7Sq/BTEHVuc1s3QCAcwCIg36+42p98pzi5T6xVJGUS5n3K4UHHI5ox87dcqoT5ZKGck1g4MVfW2tFCOLFKEEihhDFlmM4Qw6sBQWxdbm5ox0zvPPB+691/7xK7OMKlmlTZU6Y7/WAniKRLLBEQBpUTk/X+sgtfKguVg0+jPzAfTRI/NDgl27ID09mL3cmN//cNdw9Ujh7rshySQmkcYppLAFgziKDMaQtV8D6Ow0yoVaNb4ROjpjDwrXEVpzEd1qXiisF64BtK5a09TmcgwrUF63BxD5TnJIZjuNMguVqZhTSMl53UbBt5e2lc+9H/rmTlnZkS97zhswJFNI+VuPR1cdG+4udi/odRuPgaUgKOzs+sVScFgss7BMptEpN2NAViAvU7Avs3Ac6YWAUXW7VEr2PZmvOk/6CuRl38adIslk+WkavcTOOFxqfR5hWERv8vvCAEDhVPHFtvuev7xlwLJq5n24S45juW0AmEJqYdRQdbt0Wt7fM2L7t/30QF4+3zmyMIKI2EEfOVXv6N7BOai1ts8BBgAKH9MXey6Vkjd2Di12/KVIcPiwsVu28jB9/jKNTplCsiowTKJrobpn6ajeagQg+by/5w6mcHPyQev8Mnj03HYBgIvAFLhCATj4fAFiWtlV09P4xD05rF9dwP/8Z/P5oF/8IrBuHfD1rwOnT1s+1hl04D7swCmkMIk0TrelcCsGsBH70INxPLWkH+k0MJXK4OA264VQX88dTOHm5IPWuYju8xeR5wSmQJQyM195xaiO8Nkl+/HM1BVYjsmF20wijevwFJ7FNdapmBZOIYUejAMAPt0xhreXZPHOzOIfZioFPPMM0Ns7//fqsAodTzsbE4180H5WMPSifTXwnMCkTWmD18aNwHe/a3yXX53KIoHqfHwAKKLD6mHmbzRfbiGZxGxHCt9NDOJMOoOpVAZX3bMBH3Zmqm5+9tmmvxWHKZzMnIyJRj5oHSW6ff4icgTQQnQcoNRrwJEjxtG31QzODRjGY8ihiAQSKGILBvEiNmEcPdYjgGQSePZZ4NxzjZz8bBYFZBaeEvD+qF37e0rBCPsH3WT77EYA2hd2G7lwEXhR5TnBd+7UnKZsUQdnaMh27Xbhcm4yL0/fNSLnJvMLi7AL+fqlxa/SickdvCieKlYTprWGGmwWgTkCiCBzWYTpaaNKQUWVYu/nq2sdgRQKkJUroYrFhaskkUBP27t4ZyaDFSggC+OEKkdh3Le722hzqYxD1cOXrujqWjjad/piwn4w13LMX0jL2hykm90IgLWAIsaqZo6VZsvLlHWie2v/ge974CAuN3X+AIBiERcnDuIyHMMgcjiDDnTgDLa2DeLzf9GPNWuMm/X2Gj+rSujUq6lTQxN3pUbVK+JEocZF4IixygqzYldexrKcScWV5qqc61cXMHuTffXKQgF44AHrNiwtHscgcliKaZyFSSzFNB5P5JBBAddcA1x/PUvTRx7zZSNNSwBQSn1TKfWaUmpOKVW9MEG2rGqJVUqljJOylwqplZg79tWrgfvuAz7cVV6D+cNdw7hzSwGfmd6PxGQBH58Zw9QZ+z/wsTHgSGcvzlRk7kiiA7d8/6zqjJ5EAg/eNhZoNWTyESuNRpvVwoDfFwDnA/g0gP8OoM/p/bgIbDAvdHZ0GOckSaeNtdKdO43Tz5oXhJ8eyMv7e4yF1sqaN5W7Yz9qS8gppOQDLJcppORmDFjvoD18WGRkRAqHjTIJpSJqJ7BMTmG+dk4+L3MVuxhnO43ibI3uqucaY4hx5T30EMZSEAwA7lVmAZn/be5zb5zvmM8sXV5WGgEQ6cOIfFBRH8eqmuY/bR8wOvLSH/i2bWUR5qVtQ5JKiazpMmrnPD1g6qUrOofJgaFQnkSemsQIHWqRDQAAbgUwCmB09erVPr09jfHru+7F41aWTq48ei8VR7Otj2NRTfO5nSPltXksevDC4bx92yteWCMHjKzJQ9Q8uwDg2xqAUmqvUupVi8vVjTyOiDwqIn0i0pcJQVaBX+fW9upxzVOypTNcmRWRQBZjAICjyGALBnEKKcx1p42NVhULeulUEV/eml3cBXnypOWi34qTY/abJCt2UFrV3rE71wbXGIl8ZBUVgrogYlNAfh2Nev24pSPsNV3VR/hzqZT82V35siPwpwdMR+j1Ds99eBNqTfFwBEDUPLAaaPP8Ohr1+nFLR9hP7stgdqC8jogaHMT3/nUG4+PAU08ZhdK+sNl0hG51eG7mcW2SeudIZ00eIv9o2QmslLoWwJ8DyAA4DuCQiHy53v107wT2q0Kk75UnLbbGNr1506Pttvv3G9Nek4tFQZFOG/FnwwbPn44olux2ArMURINKHWciYaQ7e7Xr3fZxfej5wlTqOExtIWpVLAftkXozJJ4+rk8rzs1MOdkt1rrFKR4ifTgCCCsfD43dPrSfNb84xUPkH44AosbH/Ec3R931Fmu9aFPQ59ogijtWAw0rn2us9PcbBRudHnWX4pF51NBsxVEi0osjgLAKYHK8kaNuq3h0+rRRrp+IookBwG/NrJr6teLsomnmeJRKGdctWQKsX89yzkRRxQDgJy+yeHyaHHfTtP5+4MAB4wxkwOI5QFjOmSiaGAD8UCgAzz/v76ppk81z27STJ42SQWaszUMUTQwAXisdWm/eXH3OxpD0lM0kGPH8H0StgwHAS+ZD66mp6t+HpKdsphPnxi2i1sEA4CW7E/YuWxaqnrLZTtzntWkiCgj3AXjJ6tA6mQR+8QugtzcUnX9Jo/sAKmUyoXo5ROQCRwBesjq0fuwx4MorQ9lbOk0w8rr+DxGFAwOA11psfsSvM6ARkX4sBke2WKqZqDWwGBw1jOfjJWptDABkizn/RK2NAYBsMeefqLUxDZRqajZdlIjCiwGA6mLOP1Fr4hQQEVFMMQAQEcUUAwARUUwxABARxRQDABFRTDEAEBHFFAMAEUUPS9R6ggGAiKKFJWo9wwBARNFhPu3q5KTxM5fjSMAlBgAiig6WqPUUAwARRQdL1HqKAYCIooMlaj3FYnBEFC0sUesZLQFAKXU/gK8DOAPgdwC+IyLHdbSFiCKIJWo9oWsK6AUAF4jIhQB+C+AHmtpBRBRbWgKAiDwvIrPz/30ZwCo/n8/PPSPcj0JEURWGReAtAP7G7pdKqVuVUqNKqdGCi17Wzz0j3I9CRFGmRMSfB1ZqL4A/svjVDhF5dv42OwD0AdgsDhrS19cno6OjjttQKBgd8/T04nWpFDA+3vz0oZ+PTUTkJaXUARHpq7zet0VgEdlUp0E3AfgagI1OOn83SntGzJ10ac9Is520n49NRBQEXVlAVwH4UwBfFJFTfj2Pn3tGuB+FiKJO1xrAIwC6AbyglDqklBrw40n83DPC/ShEFHW+rQH4odE1gJJCwb89I34+NhGRFwJfAwgTP/eMcD8KEUVVGNJAiYhIAwYAIqKYYgAgIoopBgAiophiACAiiqlIpYEqpQoAxnW3o4YVAI7qboQmfO3xxNceDT0iUpWvGKkAEHZKqVGrXNs44Gvna4+bVnjtnAIiIoopBgAiophiAPDWo7oboBFfezzxtUcY1wCIiGKKIwAiophiACAiiikGAI8ppe5XSv1GKfW/lVL/SSl1lu42BUUp9U2l1GtKqTmlVKTT45xQSl2llHpdKfWGUuou3e0JklLqMaVUXin1qu62BEkpda5S6pdKqcPz3/XbdLepGQwA3nsBwAUiciGA3wL4geb2BOlVAJsB/Ep3Q/ymlGoD8BcAvgJgHYB+pdQ6va0K1G4AV+luhAazAL4vIusAXArgT6L8uTMAeExEnheR2fn/vgxglc72BElEjojI67rbEZBLALwhIm+KyBkATwC4WnObAiMivwLwvu52BE1E3hORV+b/fQLAEQAr9bbKPQYAf20B8De6G0G+WAngHdP/JxDhjoAap5TKAugF8Heam+JaLM4I5jWl1F4Af2Txqx0i8uz8bXbAGC7+LMi2+c3JaydqdUqpLgBPA7hdRD7U3R63GABcEJFNtX6vlLoJwNcAbJQW22hR77XHyLsAzjX9f9X8ddTilFIJGJ3/z0TkF7rb0wxOAXlMKXUVgD8F8A0ROaW7PeSb/QA+pZRao5TqAHAjgP+suU3kM6WUAjAI4IiIPKC7Pc1iAPDeIwC6AbyglDqklBrQ3aCgKKWuVUpNAPgcgP+qlNqju01+mV/o3wZgD4yFwP8oIq/pbVVwlFLDAH4N4NNKqQmlVE53mwJyGYBvAfjS/N/3IaXUV3U3yi2WgiAiiimOAIiIYooBgIgophgAiIhiigGAiCimGACIiGKKAYCoAfPVIN9SSv29+f+fPf//rFLqOaXUcaXUX+tuJ5ETDABEDRCRdwD8FMAP56/6IYBHRWQMwP0wcsSJIoEBgKhxDwK4VCl1O4DPA/gxAIjIiwBOaGwXUUNYC4ioQSJSVEr9CwDPAbhSRIq620TkBkcARO58BcB7AC7Q3RAitxgAiBqklLoIwBUwzgh1h1LqE3pbROQOAwBRA+arQf4URh34t2Es/P5Yb6uI3GEAIGrMLQDeFpEX5v//bwCcr5T6olLqfwB4CsDG+QqZX9bWSiIHWA2UiCimOAIgIoopBgAiophiACAiiikGACKimGIAICKKKQYAIqKYYgAgIoqp/w8IIAhOsw7PfQAAAABJRU5ErkJggg==\n", "text/plain": ["
"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["from sklearn.datasets import make_classification\n", "from pandas import DataFrame\n", "\n", "X, y = make_classification(200, n_classes=2, n_features=2, n_informative=2,\n", " n_redundant=0, n_clusters_per_class=2, hypercube=False)\n", "\n", "df = DataFrame(X)\n", "df.columns = ['X1', 'X2']\n", "df['y'] = y\n", "ax = df[df.y == 0].plot.scatter(x=\"X1\", y=\"X2\", color=\"blue\", label=\"y=0\")\n", "df[df.y == 1].plot.scatter(x=\"X1\", y=\"X2\", color=\"red\", label=\"y=1\", ax=ax);"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Split into train and test as usual."]}, {"cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": ["from sklearn.model_selection import train_test_split\n", "X_train, X_test, y_train, y_test = train_test_split(X, y)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The model..."]}, {"cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [{"data": {"text/plain": ["array([1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0,\n", " 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0,\n", " 1, 1, 1, 0, 0, 0], dtype=int64)"]}, "execution_count": 6, "metadata": {}, "output_type": "execute_result"}], "source": ["import numpy\n", "from sklearn.base import ClassifierMixin, BaseEstimator\n", "from sklearn.linear_model import LogisticRegression\n", "\n", "class TwoLogisticRegression(ClassifierMixin, BaseEstimator):\n", " \n", " def __init__(self):\n", " ClassifierMixin.__init__(self)\n", " BaseEstimator.__init__(self)\n", " \n", " def fit(self, X, y, sample_weights=None):\n", " if sample_weights is not None:\n", " raise NotImplementedError(\"weighted sample not implemented in this example.\")\n", " \n", " # Barycenters\n", " self.weights_ = numpy.array([(y==0).sum(), (y==1).sum()])\n", " p1 = X[y==0].sum(axis=0) / self.weights_[0]\n", " p2 = X[y==1].sum(axis=0) / self.weights_[1]\n", " self.centers_ = numpy.vstack([p1, p2])\n", " \n", " # A vector orthogonal\n", " v = p2 - p1\n", " v /= numpy.linalg.norm(v)\n", " x = numpy.random.randn(X.shape[1])\n", " x -= x.dot(v) * v\n", " x /= numpy.linalg.norm(x)\n", " self.hyperplan_ = x.reshape((-1, 1))\n", " \n", " # sign\n", " sign = ((X - p1) @ self.hyperplan_ >= 0).astype(numpy.int64).ravel()\n", " \n", " # Trains models\n", " self.lr0_ = LogisticRegression().fit(X[sign == 0], y[sign == 0])\n", " self.lr1_ = LogisticRegression().fit(X[sign == 1], y[sign == 1])\n", "\n", " return self\n", " \n", " def predict_proba(self, X):\n", " sign = self.predict_side(X).reshape((-1, 1))\n", " prob0 = self.lr0_.predict_proba(X)\n", " prob1 = self.lr1_.predict_proba(X)\n", " prob = prob1 * sign - prob0 * (sign - 1)\n", " return prob\n", " \n", " def predict(self, X):\n", " prob = self.predict_proba(X)\n", " return prob.argmax(axis=1)\n", "\n", " def predict_side(self, X):\n", " return ((X - self.centers_[0]) @ self.hyperplan_ >= 0).astype(numpy.int64).ravel()\n", " \n", " \n", "model = TwoLogisticRegression()\n", "model.fit(X_train, y_train)\n", "model.predict(X_test)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Let's compare the model a single logistic regression. It shouuld be better. The same logistic regression applied on both sides is equivalent a single logistic regression and both half logistic regression is better on its side."]}, {"cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [{"data": {"text/plain": ["(0.5, 0.64)"]}, "execution_count": 7, "metadata": {}, "output_type": "execute_result"}], "source": ["from sklearn.metrics import accuracy_score\n", "lr = LogisticRegression().fit(X_train, y_train)\n", "accuracy_score(y_test, lr.predict(X_test)), accuracy_score(y_test, model.predict(X_test))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["However, this is true in average but not necessarily true for one particular datasets. But that's not the point of this notebook."]}, {"cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [{"data": {"text/plain": ["array([[-0.27222367, -0.16954845],\n", " [ 0.06570281, -0.17501428]])"]}, "execution_count": 8, "metadata": {}, "output_type": "execute_result"}], "source": ["model.centers_"]}, {"cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [{"data": {"text/plain": ["array([[0.01617249],\n", " [0.99986922]])"]}, "execution_count": 9, "metadata": {}, "output_type": "execute_result"}], "source": ["model.hyperplan_"]}, {"cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [{"data": {"text/plain": ["(array([[ 1.27524081, -0.16215767]]), array([[-0.25198847, -0.58704473]]))"]}, "execution_count": 10, "metadata": {}, "output_type": "execute_result"}], "source": ["model.lr0_.coef_, model.lr1_.coef_"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Let's draw the model predictions. Colored zones indicate the predicted class, green line indicates the hyperplan splitting the features into two. A different logistic regression is applied on each side."]}, {"cell_type": "code", "execution_count": 10, "metadata": {"scrolled": false}, "outputs": [{"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+v0lEQVR4nO3deXxU1fn48c+5s2SSQMIkBDBsEcSCIqJYgxQXRCvma8Gli9raBVtaf63aWm1dailarbYurfqtikvtV63aiiKtQF3qUosEUalCtSIaQCIQQiAQksxyz++PyUxmJjOTyWRm7kzmeb9efdXMnXvnZIDznPU5SmuNEEKIwmNYXQAhhBDWkAAghBAFSgKAEEIUKAkAQghRoCQACCFEgbJbXYC+GDykQldVj7K6GEII0Sd+v2LrRgfhiy6VgtETvNhsmV+J+fF77+7SWldFv55XAaCqehQ3PLrc6mIIIUSfbNrg4MaLqmjf3z3o4iox+c61TYw/3Jvxzz//6NGbY70uQ0BCCJFhVdV+/FH1vN8XeN1KEgCEECLDytwmCxa24CwyKS41cRYFfi5zm5aWK6+GgIQQIl/NmNPO5NpOmhptVFX7La/8QQKAEEJkTZnbzImKP0iGgIQQokBJABBCiAIlAUAIIQqUBAAhhChQEgCEEKJASQAQQogCJQFACCEKlAQAIYQoUBIAhBCiQEkAEEKIAiUBQAghCpQEACGEKFASAIQQokBJABBCiAIlAUAIIQqUBAAhhChQEgCEEKJASQAQQogCJQFACCEKlAQAIYQoUBIAhBBp1dpisGmDg9YWqV5ynd3qAgghBo5VK4pZfJ0bmwP8XliwsIUZc9qtLpaIw7IQrZQarZR6SSn1H6XUBqXUpVaVRQjRf60tBouvc+PpNGjfb+DpNFi8yC09gRxm5Z+MD/ix1vowYDrwfaXUYRaWRwjRD02NNmyOyNds9sDrIjdZFgC01p9qrd/q+u99wHvASKvKI4Ton6pqP35v5Gt+X+B1kZtyom+mlKoBjgLqY1xboJRaq5Rau69ld9bLJoRITpnbZMHCFpxFJsWlJs6iwM9lbtPqook4LJ8EVkoNApYAP9Rat0Zf11ovBhYDjDtsis5y8YQQfTBjTjuTaztparRRVe2Xyj/HWRoAlFIOApX/o1rrp6wsixAiPcrcplT8ecLKVUAKeAB4T2t9m1XlEEKIQmXlHMDngAuAk5VS67r+V2dheYQQoqBYNgSktX4NUFZ9vhBCFLqcWAUkhBAi+yQAdClqaaZiwzqKWpqtLooQQmSF5ctAc8GYFUupve4KTIcDw+ulfuEtbJkzz+piCSFERhV8D6CopZna667A3tmBc/8+7J0d1C66XHoCQogBr+ADQGnjVkxHZAIT0+6gtHGrRSUSQojsKPgA0FY9GsMbmcDE8Hlpqx5tUYmEECI7Cj4AdLorqV94C74iF57SwfiKXNQvvIVOd6XVRRNiwEj1kBg5XCazZBIY2DJnHjtqZ1LauJW26tFS+QuRRqkeEiOHy2SehNUune5Kdh8+VSp/IdIo1UNi5HCZ7JBvUwiRMakeEiOHy2SHBAAhRMakekiMHC6THRIAhBAZk+ohMXK4THbIJHAeKmpplgnrAtLaYuT1ASupHhIjh8tkngSANMlWpSxpKwrLQFkJk+ohMXK4TGbJEFAajFmxlLl105l10fnMrZvOmJXPZORzJG1FYZGVMCLT8upvkrFjI6W/OZXS35xqdVFCslkpS9qKwmLFShjZeFVY8nYIKBgE2q543tpyBCvlzo7Qa8FKOd1DQZK2orBkeyXMQBluEsnL2wAQFKs3kM2gkM1KOZi2onbR5Zh2B4bPK2krBrDgSpjFi9zY7IHKP1MrYcKHm+gMvLZ4kZvJtZ0pf16+T14XgrwPALFks3eQ7UpZ0lYUlmythAkNN3V2vxYcbkrlM6U3kR8GZAAICu8dZDIYZLtS7nRXSsVfQLKxEiadw02Z6E2IzBjQASBcoonjdAQHqZRFPkvncFO6exMicwomACSSrZ5CLpJNZSIoXcNNksYhf0gAiBLdUxjIAUE2lYlo6RhuyubktegfpbW2ugxJO6SyWN92Wo1lnz+QgkFRSzNz66ZjD1u+6itysWz5aukJiLSQVUC54/yjR7+ptT4m+nXpAfRBruw9SId4m8cysX9BFCZJ45D7JACkwOq9B+ngLSnFFtb6B7B1duAtKc3K58vcw8Airf38JAEgTfKtd+A40Ia/qAh7Z/dSDb+zCMeBtox/tsw9DCyy5j9/SQBIs3xZURTYqawiX1Qq42klwnMnBdNn1C66nB21M6UnkIdkzX9+kwCQQbm8osiqtBLZzJ2Ub/JxGEXW/Oc3CQBZlGu9AyvSSkhCu9jydRglX9f852OwzQTJ+WqRYFprq1Nbd7or2X341Ky1voM9D1+RC0/pYHxFroJPaJfPef/z8ejGVSuKuaRuBDdeVMUldSNYtbLY6iJZRnoAFnJ2+Chp8+K8fhYelz0negXZIAntIuX7MEo+Hd0ocxaRJABYpLphL1PXbEcbCmVq1tWOoDHD+YpyieRO6pavwyjh8mXNf74H23TL/T7mAOTs8DF1zXbsfo3Da2L3a6bWb8fZ4Yt7T64MGYn0y8dhlHw1EIJtOkkPwAIlbV60ocDfnYZDG4qSNi8eV+9/JLm8ukikJp+GUfKZ5CmKJAHAAgdKHSgzMgeTMjUHSh1x7kgs11YXidRYNYxSaCtiJNh2kwBgAY/LzrraEUytj5wDSKb13xsJBumRbKXYn8ozFyrefF1+2l/5MmeRaZYGAKXUg8AZwE6t9WQry5JtjWPL2TW8lJI2LwdKHWmp/KPJUFFqkq0U+1N55kLFKytihNWTwA8Bcywug2U8Ljt7KoszUvnHki8TyUUtzVRsWEdRS3PWPzvZNfn9WbufK+v+QytiwgRXxIjCYGkA0Fq/Cuy2sgyFKlcDwZgVS5lbN51ZF53P3LrpjFn5TFY/P9lKsT+VZ65UvLIiRuT8HIBSagGwAKCqJOeLm3dyKbV1LiSKS7ZSTPS+3sb2c6XizcSKmFyY1xDJy/kaVWu9GFgMgRPBLC5OQbAqtXUuJIpLtlKM9771q4t6HdvPpaWI6VwRkwvzGqJvLD8SUilVA/wtmUlgq4+ELGTZCAa5dExlKquAAC6pGxGYVO3iLDK5Y/n2mM8YSK3l1hajT7+7yC45ElL0EMxF1NsqJGeHjyFXzox4XyYCglUpqmOJt0wwutIOf9+mDY4+pRkYSEsRJcVCfrJ6GehjwEnAUKXUJ8BCrfUDVpapUMTMRTS2POn3xZtA7m9gyOVEcb0NceTK2L4VCvl3z2eWBgCt9XlWfn6hCs9FFExHMbV+O7uGl0b0BJJ9X7h07D3IxURxrS0G917nxptgzXwuje1nWyH/7vlMhoDyRLLDNclINhdRf3MWwcDZmfzikyV4OyOP0Iw1xFHIaQYK+XfPVxIA8kCywzXJSjYXUaZyFuVbIGhtMVj6QBnRZyj74gxxDKSx/b4q5N89H0kAyHGpDMP0JtlcRJnKWZRvKSqaGm3YneD1hL+qOfPC1ryp7AbSiiORPhIAclw6hmFiSTYXUbZzFmUiGPS38os1welwamafcyBNJcwsWZ8v4pEAkOPSPQwTzuOyJ1WhJ/u+dEh3MEi28itqaY678iifJzgl4ZtIRAJAjstk6uhc19+homQrvzErllJ73RWYDgeGN7D3YMuceRHPytcJzv6sz8/GsFGiwCsyb+DXIgNANoZh8kFfewfJVH59yT+UjxOcqa7Pz8awUTKBV2SW1emgC46zw8eQ5vaE5//Gku3U0bkumdTWyVR+ofxDYYL5hwaCVM4bzka66vDA69y/D3tnB7WLLrckBXghk9oki9K9nFMExBsqSmbsvq16NIY3MkoYPi9t1aMzX/As6evwVTbSOuRC4j8hASBrMrGcU8QWHhBOBY6tc7P57KdjVn7pzD+Uy2PmfRm+ykZah0IIvPlAap4sydRyTtG7clcLU5afFPFa+BxCOvIPDaQx82ysesqlxH8DRSqNA8vTQfdFPqeDdnb4OGXZpkAPoIvPpnhh7ngJABbr73LTbKRCjpUq2+t08afH1uA82J2Wz4iWyz0aESlW4yC8UXPO7KmSDtpKhbycM1cF8ys5r5/VrzTXVo2ZH/A4uffcA0xbVJ2RjV3ZWPWUi4n/8k2slWzTr/0h2m7DdDgxvF6GQsxWgtQ+WSTLOXNHsmmukwkIVdV+fJ7I17IxZu7Ay0bvON7MoY1dknIi+2I1Dgy/D+X3QWegVTIGamLdKzVQlmVzV62IrS8T8snsPVi/ugitFRB4ls2mMzZmfuwvLueAx4kDL/N5gF1UUWw3c+LgFUk5YY1YjYNoOviXM4rsAxAFJzQhHyY4IZ9IrL0HwTXzPq8ikC1UYdhhcm1n3OekasucefzpsTXMsT/HWDbzBIHjNPrS22htMdi0wZHWNf3B52Z674CILdg48BW58JQOxucs6rG3RUWnsu0iTVFRcNKRXykYBBqbJ2M3H8TD4NA1ewaPQlz7fjVrVQWBbYS6T72NTLbQ5UhIa0WvZBte/1rECqstnR0Nse6TACAKTjon5IeVbsNnRt6XqaMQI3sbAYZdJ9XbyHRSODkS0nrhE+rRAWHX7Kktse6R/pkoGHs73GxsnszeDjeNY8t5Ye54Xp81mhfmjk95R3a5q4WLa6/BaWunxLEPp62di6ddzkH3z05z6cNa2WGCvY1U7rUleW8yUkk5ITKr013J7sOnJlxlJT0AURBeaajjrjW/xG748Jl2Lq69hhPGrkjLhPwJY1dw5PDV7GwbybDSbZS7Ao2tdB98059WdjZa6PmaMbWQSQ9A5L3wln2863et+SUefzEHvIPx+Iu5s/6GuO9PRbmrhQmV60OVfyzJJLBLpD+t7Gy10MvcJuMP90rlnyekByDyWryWfbidbSOxGz48YY1du+FjZ9vIhBV2JiUKAol6Cv1pZUsLXUSTACDyVnjL3uOHoTSxavVZHDvkn7jK94feF2ui1mfacdnb2Ng8OWLYJhf0tvegPzt08/FMA5E5CYeAlFJlSqnxMV6fkrkiCZFYcMjno5ZJ2I3AgshzeYzNjGWFeQZnrXyb6s17Q++PNVF7yrgnuezvS/j5Sw/y7WUv8urm09NSpnQOKwH9GjISojdxk8Eppb4M/BbYCTiAb2qt3+i69pbW+uhsFTIon5PBifSIHPJx4DcN3HovmxlLCd1r2mMl2tvb4WZn20hc9jYu+/sSPP7i0DWnrZ37585OqSeQzDBUJqTjzGRRGM4/enTMZHCJegBXA9O01lOBbwEPK6XO6roWc1eZEOFSPf0sno69g3i9/kzK/Pu7JnNdKAWHqPfxErnGMdbO3uBEbYevNNRzCArOCfRVNiaY45HegeivRHMANq31pwBa6zVKqVnA35RSo4mTV0KIoHSfflbdsJcj6zdymnlGKA/OE5yH09bJ2Z+9m9JV+yLeb/jMuDt7480JDCvd1udy5cIEc/ScwUBIyCZporMjUQ9gX/j4f1cwOAmYBxye4XKJPBaebM3hNbH7NVPrt6fcEwg+z2H6KaeVEtp5kAsZShM+087YIRtRRlSbRIGj0x+zBxJz81btNSlV2OkMJumw9vu3c8nn3fxqfgmXfN7NqpXFvd+UoqKWZio2rEv7Ob5jVixlbt10Zl10PnPrpjNm5TNpfb7olqgH8D2ihnq01vuUUnOAqzJaKpHX0n36WazneXEwwXifGbVPM8LXhGkzsJndrV0NnLSyAdMWuwcSb/NWXwWDyZ31N0TMAWSy9R+cy4gud/SqKIDF18Kxb3yBcldLWucMMnU6Wazc9rWLLmdH7UzpCWRAon+NS4F7lFK3aq39AEqp4cCtwETguswXT+SjdCRb6+15JUYbP5xzGa7y/Rzo6HndZoJCYzPjp3sud7WkpaJOVzBJRrwJ570dbtY2noBNRe7sDR+OSia1dTIyWUnLYfHZlSgATAN+BaxTSl0KHAFcBvwa+HoWyibyVH+TrUW3cGM979+1w0Jr/SOuK1B+DQrsYcPfmT5/OV3BJJFYLfw762+gzTOIB9++Cpvy0+4rjbgn3nBUf9JUZLKS7s9h8TJv0Hdx/zVorVuA73VV/i8AjcB0rfUn2SqcyF+pnn4Wr4Xb2/Max5Zj95gc8dZOTJvC7uvZA/HaDYY0t+ftaWyxJpxtys/9b12Nz3SFvVNTbG/Dr21JD0eV/ubUUOAd9NPfJZw87k8l3ZtUD4tPdkhKgkSkRPsAhgA3A7XAT4A6YDZwqdb6H9kqYDjZB9A/wTNwc7UC3Nvh5tvLXkxpfb6zw8cpyzYFTvnqogGfXaE0bB5XxtiPWtO2KskKsb4fh9GJ3fDS7hsUeq3Yvp/vTPslx1S/mnSvJFbgnXbXZXHfP2blMz0q6XTMAQT1paIuamlmbt30wJBUF1+Ri2XLV0fcm6l5i3wQbx9AolrgLeD3wPe11j7gOaXUVOD3SqnNWuvzMlNUkQnpXpaZCf1ZUhlrothnV7w7bTgtlcWc+PeGpI6AzAXxJnljTThfeNSveODtyDUZfm3rU+Ufb2jp/uvjB94tVzwfkW8+3a3pvhwWn8yQlEwux5bob/8J0cM9Wut1wAyl1HcyWiqRVn05A9dK/VlSGWui2O7TGKaJw2emdVVSJvW2qzjWhHOJc3+vq5DiBRVILfAG5xA6gc4s7EhO1CNIZkhKJpdjSzQHEHesX2t9X2aKIzIh3csyM6U/Syo9LjvrjxrGkWt3hNYuK2DyW028clpNWlclZUq8lviRw1f36AmE/9zbKqTegkp/9zKkOqEcXanHq+R7G7pJZt4gk/MW+Sx3/vWLjEn3ssy+StT6jNafJZWtFS58doXDFxnoHD4zbUdAplP0nMzOtpGMUNsZxi4aqGEXVUkPgcVbhZRMUEn3XoZklptGV+qbzjyX8Usf71HJJzt0E30EYnSrPtXJ5YHO0n8BXZvKfgfYgPu11jdZWZ6BKp1n4PZVKonSUl1SeaDUgYreENwV6PZUFqe0KilTYs3JuD3reM93OB6cOPEwnwd42jyzX7uKkx3eiQ68VTRR0tz/7ypWMIhVqR/6xEOBnltUJd+XoZve5g16CxKFyLJ/BUopG/C/wKnAJ8AbSqllWuv/WFWmgSzVZZn9keyQRrr0Fug8LrvlFT8EktodWb8Ru9k9J3NU/Xa01tghlNX0QS7kM0e/0q/vqi/DO8HAm6kFA8FgMKS5HW16Er43WMmne+imL5PLhcDKIyGPBT7UWn+ktfYAjxPIMyQyxOOys6eyOGuVYLD1GS7VrJvJStdh78nYuvdgXvxoHlv3Hpz0Pa801PHblbdxwIzcsKUBVGSSXbvdw2z3yn6Vsa95j9KdxymWWEOS0YKVfHDoxlfkwlM6GF+RS4Zu0sjK5tBIYGvYz58Q2HMQQSm1AFgAUFVifetNJM+qRGm9tfT7MicRzz1vXM2KD78W+vl/JjzCgmNuTHhPsEdUZu7HQWSrVgHRe3IMHT+jaV/0ZV4lGwsGYvXUgvs0TNegHuPzMnSTOTlfo2qtFwOLIbARzOLiiD6wIlFab9JxeMvWvQd3Vf7dLfZnN36NmWOe5bBh/457X7BHtMtfxXwe4EEuxIuDEqONf9cOA8jYPE2y8yrZWjAQa0hy4+Sq7p/fvYvSd+8KzRvI0E1mWBkAtgHhA3mjul4TA0g2E6X1Jt6cxMFD3qfDVxq3fNGrdT5ojn0i6rUv/R+XTr8ybkAJ7xE9wXm8yClMMN4PJbUDLJ+o9rjsbD64jHEfdh+puXlcWcKyRH8/ye44j+6pxeq5pSuBnYjNygDwBjBBKXUwgYr/XOB8C8sjMiQbidKSEWtFjNbwo5VP47B5YvYIoidEVx01ngrVzDG8QQMHs4uqrncqfKYz4SR3dI+o1RzEjNqnIw6wb6KKnYxkGNsoJ/vfmbPDx9iPWyPywI/9qJWNk6tiVubR38/mg8sY+3FmUm5kIhgUem4gywKA1tqnlPoB8HcCy0Af1FpvSHRPhzJ5sqKZyQdKmNiRuYMuxMAUa07Ca7oAhdcsAiJXKcXaQV27dgvH8AiX8mecYSeTBfW2bj9Rjyg4PGVTfnymg28ffQNzJvwlzd9CYn2ZA4j1/Yz7cG8geGR4x3l/spkGFXJuoCBL5wC01suB5cm+v9Hp4dGhTdi14vqtY0JB4H1XO+tLDsQMDKlck+cl97x8E90C9/gdGEpHJFcLr8BjVYZOPCighMC69Ae5kBc5JdQTSGaSO1aPKDRB7N9PDQ00UMPdaxeB0sw55Mk0fQPxOTt8lLd0YPOYKH9kJtB4cwCxvp9o2dpxHu9s5HiBQXIDBeT8JHA4DWgFPmXw1pxzGT3hB3yw+01uXHUuPtOL3djH1TMe59CKaQApXevLPdd/PDwiCF07egs+pSMCVLzXE92TD88LXs/VALW3w80a5WJH1SaO8ajQtRPGrqDk4OdZ64RJbYO4a0lkYlvPQW9Sf/Aa/B7FlKgJ0ddHwUs1MKsBjgsmSjE0I0Y9yd6Dd6A/PomLRy2LqNyD5avZM5zypokRrf7wstuaJ3Ou+Th3831eG2Xwek0n/2z4Gfe9eQ3HjXoxo0No1Q17OWr1pwRP1dQK/ApMu5FwMjqZ5ZxWp9yI11OQ3EABeRUAlFIoDOyGg0mVxwHwXvPr+EwvJn58ZuDnYIWdyrW+3BMMQgBvbbwL33u3BK6FBah4rye6Jx+el44A2uPaq2clDFAOrbht/XDG6UF4XPa4AeqVhjru2PYFfBfUgc3D01pxwycjQ8+7ueueF7XinM9/iyXP/SHQIzjoTfjGqSwxfDyjFdczhqFdyxVXjVLM+ZqJxwZOP7z4x0AQeHNUJx984/v4DI3jxOsZtnU0RAVQL6CH7KfoX7ejt36Oi2uvYdhnnooo+7Vtxdytn+Xfo9o58xt0fc4vqHh4csSQ0iZjH+8525jkKWW8OTji30cqAXSzfy/v2j/FO7I7qCkNfkPxh9Mref0gP5M8TiZ2xP6coVHLOZ8+toQPHG2csEVx3FYigkcu9Ha33jmT9SUHOGqPk9lRG8zqh3fysP85Juz2h/5eDnR5FQBGDBrHCZO+yKTK40J/QJMqj8NuOPCZRASGVK/J85J7XjoCaI9rhpE4QJma3Xt2s+Bfe6hfeAtvjd/aI0CVD72Eu+pG4PvszWDzgOHH57ex1qmZ2AHrSw7gUxpTgQ+Nfdw/uH/ubHa2jaT+4DUsMXyha+tLDjBxbCW7hpfy1JBddNj2og2NR8PKmiKmNfp5/PhifMa+QM8UM3BPV4UT/CytAO2hc9Qb0HAKd9bfwLwjlkSUY0vpdpRh8nJNoPL3G+DR4Br7XGhIaVdLE9dMa+4KDnu4/7UyKiqGJQyGiXpxu1qa+Om0ZjwT4YYTuoMawKpRmkuO3IXX0Nh1S9znwRh2DR9PSZuXN6t8XDm+EZ8Chwm/+qg6FKRyrbf7RKWi/OQKvv5cJ6ZrEPXDOzn16xpvw++xb7kvonECgQbKe82vR9Q9iV5P9VqmnoeTyJ2HXfIqALhsJczraoEGHVoxjatnPB7zl0/lmjwvuedNqjwOm+FAm2DLeICy4/f6cZow+wMP9s7AeO2rj9/L01H3NDXasDmAhpPA7wTtAdPJ8KbxUPQpNXuGY6vYDYYfu1ZMPlASGpP3exTPaBUIDF3XILA88RDKcbAXnwYDg7KKkfxjnsEhtOMw93dVlN33AIGhHdOGCWA6A2UiMM8wvGk89uFbQ581yVOKQ+3lpIZAD8OjCfy+rh2hCelPDzRHBIcGfytfWbaPd48axt+O6IwIKB+ylykdDtZXRAa8YICK9byXa7oDwKs14DV63hcdQAOvV+Jx2Xl7SHPomtfQvD3Ew/jdgefFvq844bVU7unL814eZ1A9NxC8/nDwPrx6d+Ca38+mp7/DUbsDQ0Hvu9q5MSqgjL74tagerSNBbze5a5l8HpUcSgx5FQDiObRiWtwuWyrX5Hm9X9tVPxP90AvYal5BN5zILvtEDp3TnpEAddNBv6TtsWuY/YEnVEGZdgdH76/ocU+r8uP3Ap8cB398EWpexr7lOI6ddnVolY0x+l+YY1ZxTul/mDike+PWxI5irt86JubQQfS1CaVOhjbs5eI12zlmNLw6Bg4qrWDokMh7rvxoEjdu/RL+j04NlInARPGxuoPRYc8bbxazrnYEn63fzspHFK+O0VQUVTNu8AZKmr04PH5O2KJw+nUoOMxqAJtfc+TaHVzQBI98jdAQ1QUv7eW4rXtp/HwFT1T2DGolbd4ezzuxITDPZhpwUGkldr27x32TD5RgjxEkM3EtG88L7j2Y5AG7bon5vJgB5TensqmiGd/QjphBI5VrmXweRKzsDYl7JGQuGnfYFH3Do0kvGhIZ0tpicEndCDyd3amknEUmdyzfnvAs2VQle+Rf0KqVxSxe5MZmB78PLp52OUcOX53ycZOxxDqC0mdTvDB3fI8J01c3n95jN3S8zWLhm6iGbm8LrbE3TA2mZs3IQEv9pIawiegusSapfTbFnedW8/YQT0RQC5b/jYN06Hmf3QZvHVfN7uEloWGlYICasscRKtc7Q7xpHZdPdC0XnpfrCzKSeZ65GFM3ahtRJACIPtu0wcGNF1XRvr87ABSXmlx9TxPjD/cmuDN1fT2DtrXFoKnRRlW1n4Pun83G5sn8/KUHOeDtnjgtcezjulnzmVC5vs/lGdLcznEvbcXh7Q54XofB67NGs6ey53LZ3vIP7e1ws69lKDU0YLgD2UCjA4zfCGxcs+k4zTkCLfjwa4nKVL15b4/UE7E2bSXKDprr50ynq3y5HKCSed7D9zW9rzv1pOjfK/f+xETOq6ruGmYJ4/cFXs+UviYEK3ObEb2RYaXbGOJv4TDeDx220p/EdH3NmZNoN/QrDXXsXD2Fe/VFeHFQbBzgw8PLe6yxN20Gb3yuGndzO4f+JzBebffpuMGgtzIlkyI80XGi4T2UXDxnOp1prSd2FMfcBxPv9VSvZep5D3ua2mLdJwFA9FmZ22TBwpaIYZYFC1syMvwTrq8JwcLXgE/a3sDHjKedEhx4WaDuZkTtvyl3taSUHTRdh+zs7XDzeP2lbNQTu84BaAcTDt2wGyM6O6jfpNXtYtdBg9h8iJuSNi9lLe1MfqupR1bNZMvUW+bUeDuDy1s6cvqc6Xw5B9tq8k2IlEyu7eSy25oBTc1EX8Yr//4IVQamxkkrAA8Z83lh+Dieb5jba3bQeAEiHYfs7GwbyTj1ER6coYNgAPzKhjKjcvCHxYNgxb2nspjto8riZ9XsZ2UXr6ejIafPmc6Xc7CtJt+E6LNVK4pZfJ0bmwP83kDrf8ac9t5vtEjsygDMlmIer7+UKeZ6GvyBYaHoZG69pY/ub/K2YaXb2GOW46Ij4nW79mPaFbaw841NuxGzAksmq2aq4vV0Wt0uS8+Z7o3V52DnCwkAok9aWwwWX+cOrADqDLy2eJGbybWdOdMLCE4Aj+1wU+5qiVsZ2HfY2WhOjDiD96/GGaGdt70daRkdHC486leMr3gv4VBS9KTkpO0NvKmOxaftaKCdYuyGl/8cXcnkt3f2KLMVFVi8no5V50wnw8pzsPOJfBuiTxred6CiDhK12aGp0ZYTASCid9L+YqjFHl0ZrD9qGCe//W/smBFn8E7wvx+aGE50oDrQIzjcvXYRxfY2/NrGlUf/iNnulREVZsSkpN9k02fcHPLBHmymxklguMdpdPDKnBrayovwOY2cqcBi9SqsOGe6L3K9fLlAvhGRtFUrirl3kRuvJ3LdSaZXACWrZ++kONRiZ2zkYSuBYSEgrNheHPzk8CtDrfdER1rGCg6gaPcN4lwe46dv/B92uwdDm6yrHcGu4aU9JiUPfa9nL0HbFA5fIJDmQwWWzuGmTMj18lnNykPhRR4JVq5ej0H3SnONo8jMygqgZITSQIQJb7EHJ009LnvMYaFS234mj17NkOb2QHrkOAeqV9HEJM9/GeLvWYFP5188xDcpoR2nzx86VL28pSMwDxFG0XM9f/QwT3iZhUg3+VslkhKqXDu7XytyaX506y6mHOexrFzhYu1PiLfWP9YY8cYxVZywcgvaINRyP2X4Xzn5+BdooIbB7l0csfUjjnhmJ6YBHzOe+ep+nrKdQ7uvlN9xCRdzV4/P0oYKbNCKkzpZA35bIBTIOLXIJvmbJpISq3LVGmom+mLfECZ8V24mewrh+xPsZluvh9CHD7Gs2nESl/z7zzjwQ1cRj3r9UzAUpqFQ5nq2jhlEzcf7UIDNBPDxkG0+M2YuoWlXNRevvzfmpixlalrdrkDAWf0pNjOy5e83YM3xI2l1u6TyF1klf9tEUlLd/JXtJaMz5rQzubaT/TdfGnM1TvQqHI/LThNVLHv3Qr7H0oi1+IYG5dfYusbsg5V/OK00k53vMKj0jR5lCayV727VBwPO2A9bAjt5wyZ3dx00KN1fRa9yPY2DyDz5UxdJC1auybbmrVoyWuY2OShGfp94qQF2to1kmzEKp9n3oSzDDKw599pjT6e9NmsMe4Z1Z5YMbtQK7uS1qvJNZ5oEkb9kElj0SZnbZPzh3qQq8FiTssElo9kWnhrA4TWx+zVHrt5Jx95BuOxtNJoHcR8XBlrtdB8/mogG3j06cCBLW3kRH00oj7j/ownlEZV/uOBEdEmbF2dH78No6RTru5havz3r5RDWkx6AyJjeksZla24AYu8GPmCWctvy3/GmcQzD9A6+wwMRQzyawNGIphE4JjGUZ0cFxvXXTxvGlkPcofdvOOYgGiZU4G5up6WymLbyorjlsbIFLmkSRJD8aYuMSTRvkMzcQHiAAPoVLGIt+3TgZRMT8JlOqtneIx+P0mAqjTLpquwres2z01ZelLDiB+sTlUmaBBEkAUBkVKx5g2TmBsIDhKcdlKFwFOm4waKopTmUKtp+/7k9yhFc9nnk6p0cMEtx4GU+D7CLKgAaqMFJ5ByAAuxdsWbyW01sH1WWlo1FVrfAJU2CCJI/cZFx0bn5Y+0pCE8nEStA4AefNzBAEwwWVTRR2riVIe+9y7Rbr8Nvd2B4vbw+ZRSbR1f1yOC5a3gpL08/lNtfv5k39PRQ5Q+wiyrm8wAPciEmBqW0RQ4HpbGCzoUWeD7sMs5n+bLCKndLJgas3uYGwgPEUJqooSF0iAsEgkXVkmeY+8BlaLsde9v+QGu968jIY9Zu48trX2afoyy0F+Bc/XhozP1EdTbz6d7AFVyV/wTn8SKncAyr+asxD3tYJZ3OCjpXWuCSJiEz8mmFlawCElkXnBtwFpkUl5o4o9JJBAPEd7iXrYzmBWazmTFcxS8ZShNDvE2ccd9l2Ds7cHRV/uG8OBhNIwe8g/H4i3ls9Q+ZumZHaNWL0/TxkDGfH065mmJ75EFJuxhK5REfsm76CHw2g06bHZ/NSHsF3Ti2nBfmjuf1WaN5Ye74nK0gRCRnhy+UKiTe9VxaYRUsryNOY1/Cv8ia8HH6GXMq4+4pKHObPDTnds595nIU4OoaB7qBa7mWG7jBczXtOHFG5dAPcuClgZrQz+PUR/gxsNP9GdqAmWUvc2vUOdlOWyenHfIXHtt+HvP0pYxTH/GRHsd5/JYTiH2Qe6qkBZ5fkmnZWz2/E6+8R8ARsd4jPQCRFWNWLGVu3XRmXXQ+c+umU7VkWUTlX9TSTMWGdRS1NFPU0syXV1zdM1EaUEwHV/MrHFETthpoZRAHKI6Y3AX4SI/DRuTKIWVqDHd7jGRvVwOBVM+N5hhe859EozmGO+tvYG+HG1GYkm3Z58L8DvQsr4pT10vzQ2RUa4tB2/stfOm6KwJj9F3j9DNuuJz5peewwzeCe868ny8tvRTT4cDwevjojC9i2h3YPLF35npxcDNX8DNuxIsDl93LJb7beIujaaCG3aoCGx6K7J34TDvn1f6WdQyPOeZ+wtgVHDl8dcSE8cbmyXHPAUj2zGAxsCTbss+V+Z1Y5Y1FAoDImOBSzmONrVzY6aQ8bMjGi4PhbVtox+ArT/wQO93BYcKSRxM+14GX+/gu9/FdDrF/xGajhk8Z3n1dtXP7nHPo8JWGKvVG4q96KXe1RFTsic4BEIWpLy37XFhhFau8scgQkEi71haDd153cm/XUs7328fhIHLZT3CcvoYGvET+IwrmydeAt7gE02bDb3fgKR2Mx1HMAtv9tJVW0lpUyegFn2GPsyry2TYfHb5SJlSuj6jYk82tH+8cAGn9F65gy95nU3gdBj6bStiyt/och+jyaoi5e1J6ACJtWlsMXnyyhKUPlGHYwNsZGMUPX2NvuOyYHb6Icfro4BDkLSnlzZ9ex6czZwOEJpCPp4qJjU2hZaPPPFAWcV86WuuxhoZEYcuFln1fhJf33ec2vxvrPbn9G4iQbObN6U34ap5OdyXQdVzkde6uSr9nFrUnOI9/Ok/mzlvf5I1t43n61gkU201afZU8cebvuGDpxdg6OyPuNPx+Pp05O/QZwf8vI3JjWV/OAOiL6KEhIfJt5VawvF6IuQ41f36TAhEr/83H7zl45NYhWcupn8iYFUupXXQ52mZD+f3U/+JW1teeFTgusjPWiKKmqFijTTh7oZ3O46YwBbhj9vbQ7+lyn86yBfWMX/Ioh99/B6bDieHzUr/wllCln0hvZwAI6+TLjthCJX8iOSQ8/423U6FNjcMFHW1dreoM59SP1bKPvj795z/C5u9uTEy/9ofU33ESNscIhnb23LXrcGq+u3A3JYM1NRO7h3qi00N0uiv5z7cvYdM5X01YhnjinQEgrJNPO2ILlQSAHBEz/w0Kf1vP94bnzenvZwZb4ZNXP0XtdVd0LcUMtL63zJkX8f7Ghz/A8Ef2JA2/j0P3vcM57eu4l+/gwYkTD/N5gKecX+Gks9q4Z2FF0r2XTndlnyp+kZuszngqkiN/EjkiVoK0eMLz5kSLNVcQb8w+2Ntwe5r4WP8Eu7d7KWbtosvZUTsz9P7WFoPljwzm+zE+cwh7eMD4MU5/eyid8h/tF3L8PUdw2UWHpXwimE2bVBv7canYv2s0dVHPA9lFFG3iadpK8/J7MdtbM/YxubQjVsRnyZ+EUupLwC+AScCxWuu1VpQjl8RKkBZJ4yrRmP5AKxpg0wZHREUfXqn7PHDmt1tZMORRTrr18oiWfXDMPlgxH84W2nHgDMuFb9odlDZuDQWApkYb64uOwuNzUhS2C9dvd+AdXAZFdsIX86giO66tn2BzHBY362dvqo39jKyqYPAQN0olPp7Ltv0DGDw84XsEaK1prawEvkvTkt9k7HNyZUesSMyqULweOBu416LPzznRh6d4PYE5AKcr0OK/4Md7qJnkparaz/rVRVxSNyI0rHLZjzfy2ZGbeGrRNDye7tb2S7/38BBXRGyyql10OfW3ncRwu43hnVtooIYGanosxTR8XtqqR4d6D56SsezwHc43eCiUMtmGn3/99FZaJ07G8Pa83z55ZMKsn71xKX9Slb9InlKKMpeDXVWjM/o5ubIjViRmyZ+G1vo9QP5hR4k+PAV6noIVPVdwLo9x7Y0Xolx2NnoC6+uf4DwAatjc45Qr0+7gc2/8ke+13RUxXr/Adj8P2S/EtDtCK3CGr/5nxLyAeeYdfG/phayyncwobwMnXjGEKeeUAlC/8BZqF10ecb/zYHfcE8F609pi4Bus0H6FkjojrZRSoDK/BzTf1s0XIqV179uFM/bhSr0MXJ5oCEgptQBYADB0xMhpdyxfndYyZHp9fbqfv2mDgxsvqqJ9v8FQmtjM2IgK/gDFjGUzu6iKed3nLAIF9s7OiHuuv/ptPju7IzRXADC3bnooxz6Ar8jFo39aw9YDw2L+PvFWEfX1OwgOZS175n2GDp1EVbWf0vLE99m2f9Drc0W3Tdt2sO3uH1hdDJEl8x57/02t9THRr2esGaCUekEptT7G/+b1fnc3rfVirfUxWutjBrsr0lrGVSuKuaRuBDdeVMUldSNYtbI4558fPldQQwMenBHXvTio4WOgewfuAYrpLBmMr8jFhm9fgumIvMdWYuekSRvpdFey+/CpdLorKW3ciumIHK817Q5GHNjM+MO9MSvy8PvDlbnNHve0thhs2uCgtSXyr2B4D0drhdaKpkYbpjXp1FOyeesnnHL2+Rw9q475F1+Ox5NwckcIy2QsAGitT9FaT47xv2cy9Zl9EV7RtO838HQaLF7k7lEh5drzww9T2VEypsc5tqVFHqZ8qwJH12ErTxd9heuvfpuX7/0Ty5avZtM5X+0xXm/ze0Ot/qC26tExx/Wj35eKRIExtBoqjFLg9SY3XLir2cZb7xSzq9nW+5sz5Be/vp2LvnUBb720nPLyMh7+y1OWlUWIRAo2GVysiia4QiXXnz9jTjt3LN/Ognth1TW34Cty4SkNtPDrF97CKRe7uHP5di79dTOX3baLmtmDQy3zTncl9Qt73hPdak/2fX3VW2CMtRpKa3A4eh+qfHJZGVNO+Axnfb2GKSd8hiV/Lev1nkRuvP0u7v7Dw6Gfr7/lDu75wyMJ79Fa8+rra5h3+qkAnHf2XJY//49+lUOITLFqGehZwJ1AFfCsUmqd1vq0bJaht3Npc/35wZ20TYfPZdnJn+sx9r5+dVFoSWj0Bqwtc+axo3Zmrztuk31fX/R2IHz4aiilNEppqqr9GL38Td3VbOOSq0bR3mHQ3jVtcfGVozhxxn8ZWpnad/61L53FBf/vR1z0rQswTZOnnl3BM488wPFnfDHm+++7/WaqKisoHzwYuz1Q4OoRI2jcvjOlzxci06xaBfQ08LQVnx0UveyyLytUcuH54aJ3z8baVRy9ASvZHbfp3pmbTGAMroYqH2wyZoK318rftv0DtmwrxuHQocofwG7XbNnmZGhlanmTxowaScWQct7Z8B47dzUz5bBJjBlZzT//9mTce5p3Sx4ikT8Kel1W9LLLdFfOiZ7fW96d/uitlW2lZANjmdvEbuheK/+gMSM9PeYJfD7FmJGxTxVL1gVfPoc/LXmGnU27+OoXz2Tf/jbqzv1GzPfed/vNfOaQcezdtw+fz4fdbqdx+3aqRwzrVxmEyJSCDgDQMylZNp4/ZsXSXvPu9Ec6hp8yuTw2E4F3aKWfO2/6hIuvHIXdrvH5FHfe9EnKwz9BZ3x+Nr/67f/i9Xm577c3Y7PZEvYAAI6f/lmeWfE853zhdB57ahmnnzKrX2UQIlMKPgBkW1FLM7VR5+NG593pr/4OP4WnlMhU+ulMBN5zvtDKiTP+y5ZtTsaM9PS78gdwOh3MnP5ZyssGY7MlN4H/i5/8iAsv/Qk33HYnUw6fyAVfOrvf5RAiEyQAZFlofX3YBqvovDvpkGorO5n5g1w2tNKf8ph/LKZpsnbdOzx0561J31MzZjQvPv1Y2sogRKYU7DJQq2RyfX20WBuwepPp5bH55P2Nmzj65DpOnFHL+IPHWl0cIdJOegBZFlxfH503J1dy4MeaP/ClcflqPpk4YTzrXl5pdTGEyBgJACno7wqeTKyvT5fg/ME9CyvwdQUC0wfr64ssO4YyHsn/I0T/SADoo3St4MnWyVeprOaZXNuJUprgCKHfr/JqHkAIkRwJAH2QjRU86ZTqap6mRht2J3jDltDnyj4CIUT6yCRwEopamqnYsA73++tjZsgsbdxqUcni608yukynsRBC5AYJAL0Ys2Ipc+umM+ui8zn+RxdiC1u+CZlbwdNf/VnNE55xtLjUxFlkZiyNxUC0+P/+xNGz6nCPP0JSQ4icJkNACcQa8vHb7PiKijDtzpxbwROuv634TKfJyBTVvBtjWyPmyGp0ZXrPj0jW9GlHMefkEznj/PmWfL4QyZIAkECsTVt+VzGv/fpuvGXlObeCJ1w6ktFlOk1Gf8RaAeRYtpySqxaCww5eHwduWoT3C3Upf8aNt9+Fe0g5F33rAiCQDrqqsoLvfetrCe+bcviklD9TiGySAJBAvE1beyZOztmKP1y+tuJToZp3U3LVQlRHB3TF65IrF9I6Y3rKPYFU0kFPnDA+1V9BiKyTAJBArm/aSkYut+LTydjWGGj5h0/R2O0Y2xrxpxgAUkkHLUQ+kQDQi1zetJVOmcz+mQ3myGrwRh0c7PMFXu+HvqaDlh6AyCcSAJKQrU1bVslG9s9M05UVHLhpESVXLgS7HXyBOYD+TgSnkg5aiHwhAaDA5Xv2z3DeL9TROmN6WlcBpZIO+t6HHuWO+x5kR1MzM//nHE496Xju+NWifpdFiHSTAFDgcvn0sFToyoqUx/xjSSUd9He/+VW++82vpq0MQmSKbAQrcLLrNz5JBy0GOukBFLhsHl6fDtnMACrpoMVAJwFAFNR+ASFENwkAAiic/QJCiG4yByCEEAVKAoAQQhQoCQBCJLDlk20cN+csq4sBwKgjjrW6CGKAkQAgBpQ1Teu47d37WdO0zuqi9JnP5+v9TUKkkUwCi7zR2xLQNU3rOPP57+Dxe3DanCw99T6OrZra7881TT+XXvUL1ry9joOGD+OGa37C9y6/mleW/RmATR9vZv6lV/DKsj8z5YTTOLPuNF545Z8Uu1zcd/vNjKsZw67m3Vx27fV80vgpADf+7KdMP+Yobvrd7/l4y1YatnzCqOqDmH38DP72/Iu07tvPp9t38uUzz+Cnl1wUUZ79bQf46ncvYU9rK16vl59ddjF1p57Mlk+28aX5FzF92tGhsj567x0Uu1z9/g7EwCQ9ADFgvLZ9LR6/Bz8mHtPLa9vXpuW5mxq28O0LzuX1lUspLxvMO/95j7JBg3j3P+8D8OiSpXz1nHmh95cNHsSqFU/znQvO46pf3gzAldffxEXzL+AfSx/nj7+/nUuvXhh6/383bmLpw/fxwO9+DcBb/17P//3v7by2fAlLlz/H2+9siCiPq8jJw3f/lleW/Zm/PvogP7vxFrTWMcu6bOXzafkOxMAkPQAxYMwccQxOmxOP6cVpOJg54pi0PHfsqJEccdhEAI6cfBhbP2nkgq+czaNPLuWGa67g6WdX8uJTj4Xe/8UvnA7AOV84natvCFTqr/yrnv9++FHoPfv2t7G/7QAAp58yK6KVftLnjqPCPQSAL5w2m9VvvsVRUw4PXddac/2tv2PVmjcxDINPd+xk567muGUtVM4OHyVtXg6UOvC4pKqLRb4VMWAcWzWVpafex2vb1zJzxDFpGf4BcDqdof+2GTY6/J3MnXMqv77jHk447limTj4sVGEDKKV6/LepTZ5f8iiuoqIezy8pLo74Ofx+AEXkz3955lmam1t4+ZkncDgcTDnhNDo7O+OWtRBVN+xl6prtaEOhTM262hE0ji23ulg5R4aAxIBybNVULjvi22mr/ONxFRVx8vEz+PHPf8n555wZce2pZ1eG/v+zRx0JwKyZx7H4j38KvSc4fBTLy/96nZY9e2nv6ODZ5/9B7bSjIq637tvP0MoKHA4H/3x9DVu3FW4rPxZnh4+pa7Zj92scXhO7XzO1fjvODplkjyYBQIgUfWne/2AYipOPnxHx+p69rXyu7mzufehRbrzmJwDc/POrWPfuBj5XdzbTT5vHg3/6c9znHj1lMl//fz9iZt05zJ1zasTwT/Bz163fwIzTz+Lxp5dx6PiD0//L5bGSNi/aiOw1aUNR0uaNc0fhUsHJo3ww7rAp+oZHl1tdDJElE4wWDp7wmdDP2UwEl4w773uI1n37uOayi0OvTTnhNF5a+jiVFe6UnvmnJ5fy9voN/OYX16SrmDFt2raDbXf/IKOfYRVnh49Tlm3C7u+u23w2xQtzxxfsXMC8x95/U2vdY1KsML8NIfrpa9+7lI+3bGXZIw9YXRQRxeOys652BFPrI+cACrXyT0S+EZEXcq31/8g9v4v5+juv/r1fzz3/i2dy/hfP7NczBDSOLWfX8FJZBdQL+VZETtNa91gVI/pHaw164Gd+9bjsUvH3wpJJYKXUb5RS7yul3lFKPa2UGmJFOURu69A29u1pIZ/mqXKd1prWDi+epq1WF0XkAKvC4/PAVVprn1LqZuAq4KcWlUXkqEZzEDTtxrWrCdW60+riDAzaxNO0lebl91pdEpEDLAkAWuvnwn5cDXwxmfv8fsWmDY6MnVrV2mLIqVg5xK8Mtuoy0FB695etLo4QA04uDJDNB56Id1EptQBYEPhpDDdeVIXfGzi3dsac9rQVYtWKYhZf58bmICPPF0KIXJOxOQCl1AtKqfUx/jcv7D3XAD7g0XjP0Vov1lofE1jDWkX7fgNPp8HiRW5aW9JT/NYWg8XXufF0Ghl5vhBC5KKM9QC01qckuq6U+iZwBjBbpzDLZ7NDU6MtLUM1TY02bA4gLG1KOp8vUlf6m1OtLoIQA5YlO4GVUnOA24ATtdZNyd83VENN10/ahHffBW8aEnw47HDEEaDCmvzahHWNYO7o//MHlKHALqsLkWPkO+lJvpOerPxOxmqtq6JftCoAfAgUAc1dL63WWn8v6wXphVJqbazt04VMvpOe5DvpSb6TnnLxO7FqFdAhVnyuEEKIbjLLKYQQBUoCQGKLrS5ADpLvpCf5TnqS76SnnPtO8iodtBBCiPSRHoAQQhQoCQBCCFGgJAD0QjKX9qSU+pJSaoNSylRK5dSytmxSSs1RSv1XKfWhUupKq8uTC5RSDyqldiql1ltdllyglBqtlHpJKfWfrn8zl1pdpnASAHr3PDBZaz0F+IBA5tJCtx44G3jV6oJYRSllA/4XOB04DDhPKXWYtaXKCQ8Bc6wuRA7xAT/WWh8GTAe+n0t/TyQA9EJr/ZzWOrjbeDUwysry5AKt9Xta6/9aXQ6LHQt8qLX+SGvtAR4H5vVyz4CntX4V2G11OXKF1vpTrfVbXf+9D3gPGGltqbpJAOib+cAKqwshcsJIIPxUlU/IoX/YIvcopWqAo4B6i4sSkgvpoC2nlHoBGBHj0jVa62e63tNr5tKBJJnvRAiRHKXUIGAJ8EOtdavV5QmSAEDmM5fmo96+E8E2YHTYz6O6XhMiglLKQaDyf1Rr/ZTV5QknQ0C96Mpc+hNgrtb6gNXlETnjDWCCUupgpZQTOBdYZnGZRI5RSingAeA9rfVtVpcnmgSA3t0FDAaeV0qtU0rdY3WBrKaUOksp9QlwHPCsUurvVpcp27oWBvwA+DuBib0/a603WFsq6ymlHgNeBz6jlPpEKXWh1WWy2OeAC4CTu+qPdUqpOqsLFSSpIIQQokBJD0AIIQqUBAAhhChQEgCEEKJASQAQQogCJQFACCEKlAQAIfqgK7vjx0qpiq6f3V0/1yilViql9iil/mZ1OYVIhgQAIfpAa70VuBu4qeulm4DFWusG4DcE1nwLkRckAAjRd7cD05VSPwRmArcAaK1fBPZZWC4h+kRyAQnRR1prr1LqCmAl8HmttdfqMgmRCukBCJGa04FPgclWF0SIVEkAEKKPlFJTgVMJnPD0I6XUQdaWSIjUSAAQog+6sjveTSCv+xYCE7+3WFsqIVIjAUCIvvkOsEVr/XzXz78HJimlTlRK/RP4CzC7KxPmaZaVUogkSDZQIYQoUNIDEEKIAiUBQAghCpQEACGEKFASAIQQokBJABBCiAIlAUAIIQqUBAAhhChQ/x9SuD4g259VzgAAAABJRU5ErkJggg==\n", "text/plain": ["
"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["import matplotlib.pyplot as plt\n", "\n", "def draw_line(ax, v, p0, rect, N=50, label=None, color=\"black\"):\n", " x1, x2, y1, y2 = rect\n", " v = v / numpy.linalg.norm(v) * (x2 - x1)\n", " points = [p0 + v * ((i * 2. / N - 2) + (x1 - p0[0]) / v[0]) for i in range(0, N * 4 + 1)]\n", " arr = numpy.vstack(points)\n", " arr = arr[arr[:, 0] >= x1]\n", " arr = arr[arr[:, 0] <= x2]\n", " arr = arr[arr[:, 1] >= y1]\n", " arr = arr[arr[:, 1] <= y2]\n", " ax.plot(arr[:, 0], arr[:, 1], '.', label=label, color=color)\n", "\n", "def zones(ax, model, X):\n", " r = (X[:, 0].min(), X[:, 0].max(), X[:, 1].min(), X[:, 1].max())\n", " h = .02 # step size in the mesh\n", " xx, yy = numpy.meshgrid(numpy.arange(r[0], r[1], h), numpy.arange(r[2], r[3], h))\n", " Z = model.predict(numpy.c_[xx.ravel(), yy.ravel()])\n", " Z = Z.reshape(xx.shape)\n", " return ax.pcolormesh(xx, yy, Z, cmap=plt.cm.Paired)\n", "\n", "fig, ax = plt.subplots(1, 1)\n", "zones(ax, model, X)\n", "df[df.y == 0].plot.scatter(x=\"X1\", y=\"X2\", color=\"blue\", label=\"y=0\", ax=ax)\n", "df[df.y == 1].plot.scatter(x=\"X1\", y=\"X2\", color=\"red\", label=\"y=1\", ax=ax);\n", "rect = (df.X1.min(), df.X1.max(), df.X2.min(), df.X2.max())\n", "draw_line(ax, model.centers_[1] - model.centers_[0], model.centers_[0],\n", " rect, N=100, label=\"hyperplan\", color=\"green\")\n", "ax.legend();"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Conversion to ONNX = second implementation\n", "\n", "The conversion fails as expected because there is no registered converter for this new model."]}, {"cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["MissingShapeCalculator\n", "---\n", "Unable to find a shape calculator for type ''.\n", "It usually means the pipeline being converted contains a\n", "transformer or a predictor with no corresponding converter\n", "implemented in sklearn-onnx. If the converted is implemented\n", "in another library, you need to register\n", "the converted so that it can be used by sklearn-onnx (function\n", "update_registered_converter). If the model is not yet covered\n", "by sklearn-onnx, you may raise an issue to\n", "https://github.com/onnx/sklearn-onnx/issues\n", "to get the converter implemented or even contribute to the\n", "project. If the model is a custom model, a new converter must\n", "be implemented. Examples can be found in the gallery.\n", "\n"]}], "source": ["from skl2onnx import to_onnx\n", "one_row = X_train[:1].astype(numpy.float32)\n", "try:\n", " to_onnx(model, one_row)\n", "except Exception as e:\n", " print(e.__class__.__name__)\n", " print(\"---\")\n", " print(e)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Writing a converter means implementing the prediction methods with ONNX operators. That's very similar to learning a new mathematical language even if this language is very close to *numpy*. Instead of having a second implementation of the predictions, why not having a single one based on ONNX? That way the conversion to ONNX would be obvious. Well do you know ONNX operators? Not really... Why not using then numpy functions implemented with ONNX operators? Ok! But how?"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## A single implementation with ONNX operators\n", "\n", "A classifier needs two pethods, `predict` and `predict_proba` and one graph is going to produce both of them. The user need to implement the function producing this graph, a decorator adds the two methods based on this graph."]}, {"cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": ["from mlprodict.npy import onnxsklearn_class\n", "from mlprodict.npy.onnx_variable import MultiOnnxVar\n", "import mlprodict.npy.numpy_onnx_impl as nxnp\n", "import mlprodict.npy.numpy_onnx_impl_skl as nxnpskl\n", "\n", "\n", "@onnxsklearn_class('onnx_graph')\n", "class TwoLogisticRegressionOnnx(ClassifierMixin, BaseEstimator):\n", " \n", " def __init__(self):\n", " ClassifierMixin.__init__(self)\n", " BaseEstimator.__init__(self)\n", " \n", " def fit(self, X, y, sample_weights=None):\n", " if sample_weights is not None:\n", " raise NotImplementedError(\"weighted sample not implemented in this example.\")\n", " \n", " # Barycenters\n", " self.weights_ = numpy.array([(y==0).sum(), (y==1).sum()])\n", " p1 = X[y==0].sum(axis=0) / self.weights_[0]\n", " p2 = X[y==1].sum(axis=0) / self.weights_[1]\n", " self.centers_ = numpy.vstack([p1, p2])\n", " \n", " # A vector orthogonal\n", " v = p2 - p1\n", " v /= numpy.linalg.norm(v)\n", " x = numpy.random.randn(X.shape[1])\n", " x -= x.dot(v) * v\n", " x /= numpy.linalg.norm(x)\n", " self.hyperplan_ = x.reshape((-1, 1))\n", " \n", " # sign\n", " sign = ((X - p1) @ self.hyperplan_ >= 0).astype(numpy.int64).ravel()\n", " \n", " # Trains models\n", " self.lr0_ = LogisticRegression().fit(X[sign == 0], y[sign == 0])\n", " self.lr1_ = LogisticRegression().fit(X[sign == 1], y[sign == 1])\n", "\n", " return self\n", " \n", " def onnx_graph(self, X):\n", " h = self.hyperplan_.astype(X.dtype)\n", " c = self.centers_.astype(X.dtype)\n", "\n", " sign = ((X - c[0]) @ h) >= numpy.array([0], dtype=X.dtype)\n", " cast = sign.astype(X.dtype).reshape((-1, 1))\n", "\n", " prob0 = nxnpskl.logistic_regression( # pylint: disable=E1136\n", " X, model=self.lr0_)[1]\n", " prob1 = nxnpskl.logistic_regression( # pylint: disable=E1136\n", " X, model=self.lr1_)[1]\n", " prob = prob1 * cast - prob0 * (cast - numpy.array([1], dtype=X.dtype))\n", " label = nxnp.argmax(prob, axis=1)\n", " return MultiOnnxVar(label, prob)"]}, {"cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [{"data": {"text/plain": ["TwoLogisticRegressionOnnx()"]}, "execution_count": 14, "metadata": {}, "output_type": "execute_result"}], "source": ["model = TwoLogisticRegressionOnnx()\n", "model.fit(X_train, y_train)"]}, {"cell_type": "code", "execution_count": 14, "metadata": {"scrolled": false}, "outputs": [{"name": "stderr", "output_type": "stream", "text": ["C:\\Python395_x64\\lib\\site-packages\\xgboost\\compat.py:36: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.\n", " from pandas import MultiIndex, Int64Index\n"]}, {"data": {"text/plain": ["array([1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0,\n", " 1, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0,\n", " 1, 1, 1, 0, 0, 0], dtype=int64)"]}, "execution_count": 15, "metadata": {}, "output_type": "execute_result"}], "source": ["model.predict(X_test.astype(numpy.float32))"]}, {"cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [{"data": {"text/plain": ["array([[0.44604164, 0.55395836],\n", " [0.5958315 , 0.40416852],\n", " [0.41722754, 0.5827725 ],\n", " [0.5319096 , 0.46809047],\n", " [0.47805768, 0.5219424 ]], dtype=float32)"]}, "execution_count": 16, "metadata": {}, "output_type": "execute_result"}], "source": ["model.predict_proba(X_test.astype(numpy.float32))[:5]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["It works with double too."]}, {"cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [{"data": {"text/plain": ["array([[0.44604165, 0.55395835],\n", " [0.5958315 , 0.4041685 ],\n", " [0.41722751, 0.58277249],\n", " [0.53190958, 0.46809042],\n", " [0.47805765, 0.52194235]])"]}, "execution_count": 17, "metadata": {}, "output_type": "execute_result"}], "source": ["model.predict_proba(X_test.astype(numpy.float64))[:5]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["And now the conversion to ONNX."]}, {"cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": ["onx = to_onnx(model, X_test[:1].astype(numpy.float32),\n", " options={id(model): {'zipmap': False}})"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Let's check the output."]}, {"cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [{"data": {"text/plain": ["{'label': array([1, 0, 1, 0, 1], dtype=int64),\n", " 'probabilities': array([[0.44604164, 0.55395836],\n", " [0.5958315 , 0.40416852],\n", " [0.41722754, 0.5827725 ],\n", " [0.5319096 , 0.46809047],\n", " [0.47805768, 0.5219424 ]], dtype=float32)}"]}, "execution_count": 19, "metadata": {}, "output_type": "execute_result"}], "source": ["from mlprodict.onnxrt import OnnxInference\n", "\n", "oinf = OnnxInference(onx)\n", "oinf.run({'X': X_test[:5].astype(numpy.float32)})"]}, {"cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": []}], "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}, "language_info": {"codemirror_mode": {"name": "ipython", "version": 3}, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.5"}}, "nbformat": 4, "nbformat_minor": 4}