Lagrange-Hardcoded-Binary.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.00000000e+00  1.01325000e+05  2.50000000e+03  3.00000000e-01\n",
      "  7.00000000e-01  1.00000000e+00 -2.20583626e+05 -1.94101859e+05\n",
      " -1.11878832e+02  2.14055264e+05]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAD8CAYAAAC/1zkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8VfX9x/HXJzuMECBhhRGQMMKGCLgVEBAHLhQnQ6VatbZqW0db7M8ua6vW1tbSypQhOAqKi+FkhxFIwgrICAQSCAkjIet+f3+cE71iLgnk5p57cz/Px+M+knzvGZ8jyDvfMz5XjDEopZRS3hDidAFKKaXqDw0VpZRSXqOhopRSyms0VJRSSnmNhopSSimv0VBRSinlNRoqSimlvEZDRSmllNdoqCillPKaMKcL8LW4uDiTmJjodBlKKRVQ1q9ff8QYE1/dckEXKomJiaSmpjpdhlJKBRQR2VuT5fT0l1JKKa/RUFFKKeU1GipKKaW8RkNFKaWU12ioKKWU8hoNFaWUUl6joaKUUsprNFSUUioI/OfL3Xy2PbfO96OhopRS9VzeiRJe/GQ7n23TUFFKKVVLc9bso7TCxfiLE+t8XxoqSilVj5WUVzBr9V6GdGtBp/hGdb4/DRWllKrHPkjL4cjJEiZe0tEn+9NQUUqpesoYw9QV39ClZSMu6dzcJ/usVaiIyBgRyRARl4ikuI2Hi8gMEdkiIltF5Gm390aKyHYRyRKRp9zGO4rIGhHZKSJviUiEPR5p/5xlv5/ots7T9vh2ERlRm2NRSqn6Zu03+WQcPM7ESzoiIj7ZZ21nKunAzcCXZ4yPASKNMb2AAcCPRCRRREKB14BrgGTgDhFJttd5AXjZGJMEHAPus8fvA44ZYzoDL9vLYa83FugBjAT+aW9fKaUUMHXFNzRtEM6N/RJ8ts9ahYoxZqsxZntVbwENRSQMiAZKgePAQCDLGLPbGFMKzANGixWhQ4C37fVnADfa34+2f8Z+f6i9/GhgnjGmxBjzDZBlb18ppYJexsFCPs08zN2DOxAV7rvft+vqmsrbwCkgB9gH/MUYkw8kAPvdlsu2x5oDBcaY8jPGcV/Hfr/QXt7Ttn5ARCaJSKqIpObl5dX+6JRSys/9+ePtxESFc/9lnXy632o/+VFElgKtqnjrWWPMQg+rDQQqgDZAU+AreztVndQzZxnnPNf5/qAxU4ApACkpKVUuo5RS9cXKrCN8sSOPZ0d1p0l0uE/3XW2oGGOGncd27wQ+NsaUAbkisgJIwZpZtHNbri1wEDgCxIpImD0bqRwHawbSDsi2T6c1AfLdxs/cllJKBS2Xy/Cnj7eREBvNPRd18Pn+6+r01z5giFgaAoOBbcA6IMm+0ysC60L7ImOMAT4DbrXXHwdUzoIW2T9jv7/cXn4RMNa+O6wjkASsraPjUUqpgPBheg6bswt5/OouPr2WUqm2txTfJCLZwEXAYhH5xH7rNaAR1t1h64BpxpjN9izkEeATYCsw3xiTYa/zS+BxEcnCumbyhj3+BtDcHn8ceArAXm8+kAl8DDxsjKmozfEopVQgKymv4MVPttOtVWOf3vHlTqxf+oNHSkqKSU1NdboMpZTyun9+nsWfP97OjIkDuaJLvFe3LSLrjTEp1S2nT9QrpVQ9cLCgmL8vy2J4ckuvB8q50FBRSql64PeLt+Iyhl9fl1z9wnVIQ0UppQLc1zuPsHhLDg9f1Zl2zRo4WouGilJKBbDScheTF6XTvlkDJl3u2wcdq6KhUkPGGFbuOsLpMr3BTCnlP17/Yhe78k7x3A3JjtxCfCYNlRpav/cYd/5nDe+n6fOVSin/kJV7gn8sz+L6Pm0Y0q2l0+UAGio1NqBDU7q0bMSMVXsIttuwlVL+x+UyPPXOFhpEhjL5emcvzrvTUKkhEeHeixJJP3CcDfsKnC5HKRXkZq/dR+reY/zq2mTiGkU6Xc63NFTOwU39EmgcFcaMlXucLkUpFcQOFhTzwkfbuCwpjlv6O/PkvCcaKuegYWQYYwa048MtOeQeP+10OUqpIORyGX7x9mZcxvCHm3r57BMda0pD5Rzde1EHKoxh9pp9TpeilApCs9fs5eusI/zq2mTHn0mpiobKOUqMa8iVXeKZs3YfpeUup8tRSgWRPUdO8YcPt3F5l3juGNiu+hUcoKFyHsZdnEjeiRK9vVgp5TMVLsOTC9IIDxX+fEtvvzvtVUlD5TxcnhRPcusY/rZsp85WlFI+8foXu0jde4zfju5BqyZRTpfjkYbKeQgJEX4+oiv78ouYn7rf6XKUUvXcpv0FvLxkB9f1bs2Nff3rbq8zaaicpyu7xpPSoSmvLtuprVuUUnXmZEk5j83bSMuYKH7vh3d7nUlD5TyJCL8Y2Y3cEyX63IpSqs48tyiD/flFvHx7X5pEhztdTrU0VGphYMdmXNElnn99sYvjp8ucLkcpVc8s3HSAt9dn8/BVnRnYsZnT5dSIhkot/XxEVwqKynj9811Ol6KUqkd2553kmXe3MKBDU34yNMnpcmpMQ6WWeiY04aZ+Cfz3q2/Yd7TI6XKUUvXA6bIKHp6zkfCwEP5+Rz/CQwPnn+rAqdSPPXVNN8JChd8tznS6FKVUPfD8B5lszTnOS7f1oU1stNPlnBMNFS9oGRPFw1d15tPMw3y984jT5SilAtjCTQeYvWYfP7qik998Rsq50FDxkvsu7Uj7Zg347fsZlFXoA5FKqXO3/dAJnnpnCxcmNuXJ4V2dLue8aKh4SVR4KM9e252duSeZtWqv0+UopQLM8dNlPPjmehpFhfHanf0D6jqKu8Cs2k8NT27J5V3i+eun28kpLHa6HKVUgHC5DE/MT2N/fhH/vKs/LWL8tw1LdTRUvEhE+N3onlQYw+SFGU6Xo5QKEK99lsWSzMM8M6o7FyYGxvMonmioeFn75g14bGgXPs08zCcZh5wuRynl55ZkHuavS3ZwY982TLgk0elyak1DpQ7cf1lHurVqzOSFGZwsKXe6HKWUn9p5+AQ/e2sTvRKa8Cc/bmd/LjRU6kB4aAh/vLkXh0+c5s8fb3O6HKWUHyosKuOBmalEhYcy5d4BRIWHOl2SV2io1JF+7Zsy7qJEZq7ay6pdR50uRynlR8oqXDw8ZwMHCor59z39ad0ksB5wPBsNlTr0i5FdSWzegF+8k8YpPQ2mlAKMMUxelMHXWUf4w029GNAhsC/Mn0lDpQ41iAjjxTF9yD5WzJ8+0tNgSimYtmIPc9bs46ErL2BMin9+znxtaKjUsQsTmzHxko7MWr2XFVnawkWpYLZ822F+tziTkT1a8fMAfWK+OhoqPvDzEV3pFNeQny9Io7BIP3dFqWC0JbuQR+ZspEebJrx0ex9CQgL/Tq+qaKj4QFR4KC/f3pfcEyU8878tGGOcLkkp5UP784uYMH0dTRtE8Mb4FBpEhDldUp3RUPGRPu1i+dnVXVi8OYd3NxxwuhyllI8UFpUxftpaSssrmD7hQlo0DtwWLDWhoeJDD15xAQM7NuM3C9PZe/SU0+UoperY6bIKHpiZyv78Yqbcm0JSy8ZOl1TnahUqIjJGRDJExCUiKW7jESIyTUS2iEiaiFzp9t4AezxLRF4V+xFSEWkmIktEZKf9tak9LvZyWSKyWUT6u21rnL38ThEZV5tj8YXQEOHl2/sSEiI8Nm8TpeXaIl+p+qq8wsVP5m5k3d58/nJbHwZ3au50ST5R25lKOnAz8OUZ4w8AGGN6AVcDfxWRyn39C5gEJNmvkfb4U8AyY0wSsMz+GeAat2Un2esjIs2AycAgYCAwuTKI/FlCbDQv3NKbTfsLePETvc1YqfrIGMOvF2bwaeZhJl+XzA192jhdks/UKlSMMVuNMdureCsZKxgwxuQCBUCKiLQGYowxq4x1tXomcKO9zmhghv39jDPGZxrLaiDW3s4IYIkxJt8YcwxYwncB5ddG9WrNuIs68J+vvmFJ5mGny1FKednLS3Ywd+0+Hr7qAsZf0tHpcnyqrq6ppAGjRSRMRDoCA4B2QAKQ7bZctj0G0NIYkwNgf21hjycA+6tYx9N4QHjm2u70TIjhyQVpZB8rcrocpZSX/Per3by6PIvbU9oF7Kc31ka1oSIiS0UkvYrX6LOsNhXrH/lU4BVgJVAOVHVjdnX313pap8bbEpFJIpIqIql5eXnV7M43IsNCee3O/rhchofnbKSkvMLpkpRStfTWun38bvFWru3Vmj/c3KtedB0+V9WGijFmmDGmZxWvhWdZp9wY8zNjTF9jzGggFtiJFTRt3RZtCxy0vz9sn9bC/pprj2djzXLOXMfTeFX1TDHGpBhjUuLj46s7ZJ/p0LwhL47pTdr+An77fqbT5SilamHx5hyefncLl3eJ5+Xb+xJaTx9urE6dnP4SkQYi0tD+/mqg3BiTaZ/WOiEig+27vu4FKsNpEVB5B9e4M8bvte8CGwwU2tv5BBguIk3tC/TD7bGAMrJnax668gLmrNnH/HX7q19BKeV3lmYe5rF5G+nfvimv392fiLDgfVqjVo91ishNwN+BeGCxiGwyxozAuh7yiYi4gAPAPW6rPQRMB6KBj+wXwJ+A+SJyH7APGGOPfwiMArKAImACgDEmX0SeB9bZy/2fMSa/NsfjlCeHd2VLdiG/WphOt9aN6d021umSlFI19MWOPH48ewM92sQwdcKF9fpp+ZqQYGsZkpKSYlJTU50u4wfyT5Vy/d+/xhjDwkcuJb5xpNMlKaWqsXLXESZMW8cF8Y2Y+8BgmjQId7qkOiMi640xKdUtF7xzND/TrGEE/75nAMeKyvjRrFS9cK+Un1u16yj3TU+lQ/MGvHn/oHodKOdCQ8WP9Exowl9v68OGfQU88266Np5Uyk+t2nWUidPX0bZpNLPvH0yzhhFOl+Q3NFT8zKherfnpsCTe2ZDNf77a7XQ5SqkzuAfKnAcG66nqMwT3FSU/9ZMhSew8fJI/frSNDs0bMqJHK6dLUkoBX+88wv0z19GuaQMNFA90puKHQkKEv4zpQ5+2sTw2byOb9hc4XZJSQW/5tsNMnLGOxOYNmTtJA8UTDRU/FR0Ryn/HpRDfOJL7Z6xjf762clHKKR+n5/CjWevp2rIx8yYNJq6RBoonGip+LK5RJNPGD6SswjB+2lqOnSp1uiSlgs4767N5eM5GereNZfYDg4htoBflz0ZDxc91btGIKfcMYP+xYibOWEdRabnTJSkVNKat+IYnFqQxuFMzZk4cSEyU3jZcHQ2VADCoU3NeHduXtP0FPDJnI2UV+uFeStUlYwyvLN3Bb9/PZESPlkwdfyENI/W+pprQUAkQI3u25nc39mL5tlyeemcLLpc+w6JUXahwGSYvyuCVpTu5pX9bXruzP5FhoU6XFTA0egPInYPac+RkCS8t2UHjqDAmX58clK21laorJeUVPP5WGou35DDp8k48NbIbIUHabfh8aagEmEeHdKawuIw3vv6GRpFhPDki+D4ESKm6cPx0GQ/OWs/KXUd5dlR3Hri8k9MlBSQNlQAjIvzq2u4UlZbzj8+yaBgZxkNXXuB0WUoFtJzCYiZMW0dW7kleuq0PN/dvW/1KqkoaKjVVdhq+eAH63Q3Nnf1HXET43Y29KCqt4IWPtxERFsJ9lwbX52Ar5S2ZB48zYfpaTpVUMH3CQC5NinO6pICmoVJTB9bDylfh65eg4xWQMhG6XQuhztxiGGo/dV9a7uL5DzIJEZhwiQaLUufiix15PDx7A40iw1jw4EV0bx3jdEkBT+/+qqnES+BnGTDkV5C/GxaMg5d7wLLnoWCfIyWFh4bw6h39GNGjJb99P5MZK/c4UodSgWjmqj1MmLaWds0a8N7DF2ugeIl+SNf5cFVA1lJInQo77E8wThpuzV6SroYQ395+WFru4pE5G/g08zC/uS6ZiXoqTCmPyitc/G7xVqav3MPQbi149Y5++gxKDdT0Q7o0VGqrYB9smGm9Th6GJu2g/zjofw809l134dJyF4/N28hH6Yf4xciu/PjKzj7bt1KBorCojEfmbuCrnUe479KOPDOqO6F6y3CNaKh4UGcfJ1xRBts/tGYvuz+HkDDoOsqavXS8AkLq/kxjeYWLJxaksXDTQX4ypDM/u7qLPseilC0r9yQPzEwl+1gRz4/uydiB7Z0uKaDUNFR0zuctoeGQPNp6Hd0F66fBxtmwdRE06wQDJkDfu6Bh8zorISw0hJdu60tUWCivLs/iREk5v742WR/eUkFv+bbDPDZ3ExFhIcx5YDAXJjZzuqR6S2cqdanstBUqqVNh3yoIjYDkG63ZS/vBUEezCJfL8PziTKat2MNN/RL48629CQ/VezJU8HG5DH9fnsUry3aQ3DqGKfemkBAb7XRZAUlnKv4gPAp632a9Dmdas5e0ebBlPsR3t8Klz+0Q1cSruw0JEX5zXTLNG0bwl093UFhcxmt39ic6QvsXqeBx/HQZj7+VxtKth7m5XwJ/uLkXUeH6/0Bd05mKr5WegvR3rNnLwY0Q3gB63mIFTEJ/r+/uzdV7+fXCdPq2i+WNcRfSrKF+FoSq/zIOFvLj2Rs4cKyYX13bnXEXJ+r1xVrSC/UeOB4q7g5ssGYvW96GsiJo3dcKl163QkRDr+3m4/QcHpu3iTax0UyfcCEdmntv20r5m/nr9vPrhenENgjntTv7k6LXT7xCQ8UDvwqVSqcLYfN8a/aSmwmRMdD7dkiZAC17eGUX6/fmc/+MVEJE+O+4FPq1b+qV7SrlL06VlPObhRm8syGbSzo3529j++nH/nqRhooHfhkqlYyB/WutcMl4DypKoN1ga/aSPNq6RlMLu/NOMn7aOg4dP81fxvThhj5tvFS4Us7amnOcR+ZsYPeRUzw6JInHhibp8ydepqHigV+HiruifNg0G1KnQf4uiG5q3ZI8YALEnf+DjfmnSnlw1nrW7snnsaFJ/HRYkp5rVgHLGMOba/bx/AeZNIkO52+39+XiztoQsi5oqHgQMKFSyeWCPV9as5dti8FVXuuGliXlFTzzbjrvbMjm2l6teXFMbxpE6I2AKrAcPVnCL9/ZwtKth7m8Szwv3dZHT3fVIQ0VDwIuVNydOAwbZ8H6GVC4Dxq1hH73wIBxEHtuTwcbY5jy5W5e+HgbXVo2Zso9KbRv3qCOClfKu77ckccTC9IoLCrjqWu6Mf7iRH3It45pqHgQ0KFSyVUBWcus2cvOT6xrMefZ0PKLHXk8OmcDIsLf7+jH5V3i67BwpWqnqLScP364jVmr99KlZSP+Nrafdhf2EQ0VD+pFqLgr2O/W0PIQxLSFAePPqaHlniOn+NGs9ezIPcFPh3bh0SGd9bc+5Xc27DvGE/PT2HP0FPdd0pEnR3TVhxl9SEPFg3oXKpUqymD7R3ZDy8/OuaFlUWk5z76XznsbD3B5l3heub2vPiip/MLpsgpeWrKD/361m9ZNovnLmD5cdEHd9dBTVdNQ8aDehoq7o7tg/XTY+CYU59e4oaUxhrlr9/PcogyaNYzglbF9GdxJ/+dVzlm/N5+fL9jM7iOnuGNge54Z1Y3GUc582mqw01DxIChCpVLZadj6vt3QcmWNG1qmHyjk0bkb2XvUuuf/0SGdCdOGlMqHTpaU8+LH25i5ei9tmkTzwi299bPjHaah4kFQhYq73K3WMy9pc6HkeLUNLU+WlPObhem8u+EAKR2a8vLtfWnXTO8OU3Vv2dbD/Op/6Rw6fppxFyXy5IiuNNJPZnSchooHQRsqlUpPQfq7dkPLDdU2tHxvYza/+V8GLmOYfEMPxgxoqw9LqjpxoKCY3y7K4NPMw3Rt2Zg/3tKL/tpOyG9oqHgQ9KHi7uBGa/ayZcH3G1r2vAUiG327WPaxIp6Yn8aab/K5Orklv7+pJy0a165ljFKVSstdTF/5Da8s3YnLGH4yNIn7L+1ERJiecvUnGioeaKhUoQYNLV0uwxtff8OLn24nOjyU525I5sa+CTprUbXy5Y48nns/g915pxjWvQWTr++hp1n9VE1DpVa/CojIiyKyTUQ2i8h7IhLr9t7TIpIlIttFZITb+Eh7LEtEnnIb7ygia0Rkp4i8JSIR9nik/XOW/X5idftQ5yiqCQx8AB5aCRM/tW5F3jAT/nUxvDEC0t4ipKKEBy7vxIc/uYwL4hvys7fSuH9GKgcKip2uXgWgPUdOMWlmKvdOXYvLZZg6PoX/jrtQA6UeqNVMRUSGA8uNMeUi8gKAMeaXIpIMzAUGAm2ApUAXe7UdwNVANrAOuMMYkyki84F3jTHzROR1IM0Y8y8R+THQ2xjzoIiMBW4yxtzuaR/GmIqz1awzlRoqyodNc6zZyxkNLSuaXcC0Fd/w1093IAJPDO/K+IsTtSusqlZhcRn/WL6T6Sv3EB4awsNXdeb+yzoSGaYPMfo7n5/+EpGbgFuNMXeJyNMAxpg/2u99AjxnL/qcMWaEPf60PfYnIA9oZQfURZXLVa5rjFklImHAISAeeKqqfRhjVp2tTg2Vc2QMfFPZ0PIDt4aWE9jf4ip+/cEOPt+eR8+EGJ4f3VM/p0VVqaS8gtmr9/H35TspKC5jzIC2PDm8Ky1i9NpcoHDiM+onAm/Z3ycAq93ey7bHAPafMT4IaA4UGGPKq1g+oXIdO3AK7eXPtg/lLSLQ6Qrr5d7QcsF42jVswbR+97C82zU881khN/1zJbeltOWXI7vRXLvFKqxrce9vPshfPt3O/vxiLuncnKev6U7PhB/exq7qh2pDRUSWAlU1kXrWGLPQXuZZoByYXblaFcsbqr6GY86y/Nm2dbZ1vkdEJgGTANq3P7duvspN45Zw+ZNw6c++bWgpK15mqHmJFRcM472QETy7wfBx+iF+MjSJey9K1Dt4gpQxhiWZh3lpyQ62HTpBcusYZk7sxWVJcXpzRz1XbagYY4ad7X0RGQdcBww1351LywbauS3WFjhof1/V+BEgVkTC7NmK+/KV28q2T381AfKr2ceZxzAFmALW6a+zHY+qgZBQ6DLcetkNLcM2zGTMySXc2KwNi0KH86fFFzJr9V6evqYbI3q00n9IgoQxhs935PHK0p2k7S+gY1xD/ja2L9f3bqNNSoNEbS/UjwReAq4wxuS5jfcA5vDdRfRlQBLW7GIHMBQ4gHWh/k5jTIaILADecbtQv9kY808ReRjo5Xah/mZjzG2e9qEX6h1yRkNLI6F8HTqQfxddSXHCJfzimmQGaR+xessYw/Jtuby6bCdp2YUkxEbz2NAkbu6foC1+6gmfXKgXkSwgEjhqD602xjxov/cs1nWWcuCnxpiP7PFRwCtAKDDVGPN7e7wTMA9oBmwE7jbGlIhIFDAL6Ic1QxlrjNl9tn2cjYaKD9gNLc3GN5HifPbTilllV5HT8WYeGDmQ3m1jq9+GCgjlFS4Wb8nhX5/vYtuhE7RrFs0jV3Xmpn5t9dRnPaMPP3qgoeJDdkNL17o3CNm/ilLC+LBiINsSxjDq2pvo3U7vFAtUp0rKeXt9Nv/9ejf784vp3KIRD15xAaP7tiFcZyb1koaKBxoqDsndSumaNzCb5hJZcZIdrgTWNL+RbiMeIKVrol5zCRA5hcXMWLmXOWv2cvx0Of3ax/LQFRcwrHtLvWZSz2moeKCh4rDSUxRvXMDxr6fQ8kQGxSaCFdFX0uDiBxh06dX6AKUfMsawenc+M1ft4dPMwxhjGNmzFfdd2okBHXS2GSw0VDzQUPEfJfvWs/fT12iX/QHRlLBdOpHX7S76XHMfjWP0HyunFRSV8u6GA8xdu4+duSdpEh3O2AvbcffgDtpOJQhpqHigoeJ/yk8dY/uSqTRKn0mH8j2cNNGkx40k/qqHuKDnIKfLCyoVLsPKXUd4e302H6UforTcRZ92sdw1sD3X92lDdIS2UwlWGioeaKj4MWPI2rCMwi+n0LNgOZFSxtaw7hxLvpuew8cR06ix0xXWW9sPnWDhpgO8t/EAOYWniYkKY3TfBMYObEePNvr0u9JQ8UhDJTAUHj3E1o//Tdtd82jrOsgx04iNzUYRfdF9pAwYqHcYecHuvJN8lH6I99MOsu3QCUJDhMuS4rh1QFuGdW9JVLjOStR3NFQ80FAJLMblYve6jyha+R+6FX5JOBWspSdZHW6jw0VjGJTUSh+uqyFjDJk5x1mamctH6TlsO3QCgP7tY7mxXwKjerUmTnu2KQ80VDzQUAlcpQU57Fv6b5pum0Pz8sPkmSYslCEcThpL/959uKxLvH6W+RlOlZSzatdRPt+Ry/KtuRwsPI0IDGjflFG9WjOyZyvaxEY7XaYKABoqHmio1AOuCkq3L6Hgy9eJy/kCMHxe0Ye3zDBOdRjCJUmtuCwpjuTWMUH37ERZhYvN2QWs2nWUFVlHSd2bT1mFITo8lEuT4rg6uSVDurXQGYk6ZxoqHmio1DMF+3Gtn0F56gwiinPJlThmlV7JWxVXUdagBRcmNmNQp+ZcmNiU7q1j6t21mOOny0jbX0DqnmOk7s1n474Cikqt9nfdW8dweVIcV3SJZ0BiU/0gLFUrGioeaKjUU2c0tHRJKOmNLmF66RDeK+yMIYTIsBB6JTShb7tYeiTEkNy6CRfENwyYazL5p0rZmnOcrTnHyTh4nLTsAnbnnQIgRKBbqxhSEptyUafmDOrUnGYNIxyuWNUnGioeaKgEAbuhJRvfhOJ8ypsksqPtrXwUNoQVOZBx8Dgl5S4AIsJC6BTXkKSWjekc34jEuAa0b9aADs0b0rRBuM/bx5wqKedAQTF7jxax9+gp9hw9RVbuSbJyT3LkZOm3y7VoHEnvtrH0bdfE+to+lpiocJ/WqoKLhooHGipBpLwEMhdZs5d9KyE0ApJHU95/PLuje5ORc5ytOSfIyj3JztwT7M8v/t7qUeEhtG4STesmUcQ3jqR5w0iaN4ogJjqcmKgwYqLCiQoPJSo8hKjwUMJDhRCxXgaocLkodxlKy10Ul1ZQXFbByZJyCovLKCwuo6CojCMnSzhyspTc46fJKTxNYXHZ92poHBVGUotGdG7RiKQWjeneOoburRvrJ2sqn9NQ8UBDJUjlboXUaZA2D0oKIb4bpEyE3rdDtNWK/3RZBfvzi6xZQn4ROQXF5Bw/TU5BMUdPlXL0ZCknS8qr2VHNRYSFEN8DoG2CAAAQOElEQVQokrhGEcQ3jrQCLDaKhNhoOjRvSIdmDYh1YLakVFU0VDzQUAlypUWQ8a41ezmwHsKiodctVsC06Q/V/AN+uqyC46fLOF5czonTZRSXVVBS5uJ0WQVlLoMxhgqXIUSE0BDrFRkWQnR4KFERoTSKDKNJdLg9ywnRwFABQ0PFAw0V9a2Dm6xw2fI2lJ2C1n2scOl5K0Q2cro6pfyKhooHGirqB04Xwub51umx3AyIaAx9bocBE6BVT6erU8ovaKh4oKGiPDIG9q+F9dMg/V2oKIF2g6zZS/JoCNcnz1Xw0lDxQENF1UhRPqTNtU6PHc2C6KbQ9y5r9hLX2enqlPI5DRUPNFTUOTEG9nxlhcvW98FVDh0vt2YvXa+FMH3AUAWHmoaKdt9T6mxErBDpeDmcOAyb3oTU6bBgPDRsAf3vgf7joGkHpytVyi/oTEWpc+WqgF3LrdnLjo+t2UzS1dbsJWk4hGiPLVX/6ExFqboSEmqFSNLVUJgNG2bC+hkwdyzEtIUB46DfPRDT2ulKlfI5nako5Q0VZdasJXWqNYuRUOg2ypq9dLwSQgKjaaVSnuhMRSlfCg2H7tdbr6O7YMMMq6Hl1vehaUdImWDdPdYwzulKlapTOlNRqq6Ul1ihkjoV9q74tqElKROh/UXVtoRRyp/oTEUpp4VFQq9brVfuNuuhyk1zYcuCKhtaKlUf6ExFKV+qZUNLpZyiMxWl/FFEA+h3t/U6uMmavWxeYF1/0YaWqh7QmYpSTjt9HLbMh3VTtaGl8lvapsUDDRXlt4yB7HXWqbHKhpZtB1qzlx43akNL5SgNFQ80VFRAOLOhZVSsdUtyygSIS3K6OhWENFQ80FBRAUUbWio/oRfqlaoPztrQMt5qBzNgHDRNdLhQpSw6U1Eq0FTV0LLzsO8aWobq74rK+3SmolR95amh5bw7ICbBasXf/x6IaeN0pSoI6UxFqfqgqoaWXa+xZi+drtKGlqrWdKaiVDBxb2iZvxvWT7ceqNz2gXW9ZcAE64FLbWip6litfn0RkRdFZJuIbBaR90Qk1h5vLiKfichJEfnHGesMEJEtIpIlIq+KWH0pRKSZiCwRkZ3216b2uNjLZdn76e+2rXH28jtFZFxtjkWpeqNZJ7j6/+DxrXDLG9YpsaWT4aXu8PZ9sGeFdR1GqTpQ2znxEqCnMaY3sAN42h4/DfwaeLKKdf4FTAKS7NdIe/wpYJkxJglYZv8McI3bspPs9RGRZsBkYBAwEJhcGURKKb5raDnhQ/jxGutU2M4lMH0UvDYIVr8OxQVOV6nqmVqFijHmU2NMuf3jaqCtPX7KGPM1Vrh8S0RaAzHGmFXGupgzE7jRfns0MMP+fsYZ4zONZTUQa29nBLDEGJNvjDmGFXCVAaWUcteiG1zzAjyxDUa/ZvUW+/iX8Ndu8L+HIXu9zl6UV3jzmspE4K1qlkkAst1+zrbHAFoaY3IAjDE5ItLCbZ39VazjaVwp5Ymnhpab3oRWva3ZTK8x2tBSnbdqZyoislRE0qt4jXZb5lmgHJhd3eaqGKvu1yNP69R4WyIySURSRSQ1Ly+vmt0pFSTa9IXr/2bNXq79q/X8ywc/tWYvHzwOh9KdrlAFoGpnKsaYYWd7375Afh0w1FR/f3I29ikyW1vgoP39YRFpbc9SWgO5buu0q2KdbODKM8Y/93AMU4ApYN1SXE2NSgWXqBi48H5Iue+7hpYb34TUN7ShpTpntb37ayTwS+AGY0xRdcvbp7dOiMhg+66ve4GF9tuLgMo7uMadMX6vfRfYYKDQ3s4nwHARaWpfoB9ujymlzocItBsIN71uzV5G/AGK8+F/D1qzl4+fgSM7na5S+blaPfwoIllAJHDUHlptjHnQfm8PEANEAAXAcGNMpoikANOBaOAj4FFjjBGR5sB8oD2wDxhjjMm3w+cfWBfhi4AJxphUex8TgWfsff/eGDOtupr14UelzkFVDS0TL7NmL92u04aWQUS7FHugoaLUeTqZCxtnWQ9WFuzThpZBRkPFAw0VpWpJG1oGJW3TopSqG9rQUp2FzlSUUrWnDS3rPZ2pKKV85wcNLWdY11+0oWXQ0ZmKUqpulJdYd4ylToW9KyA0ArrfYM1eOlxs3cKsAobOVJRSzqpsaNnrVsjdZrWE2TQX0t+GuK5WuPS5HaK1D2x9ojMVpZTvlBZBxruQOg0OpEJYNPS8xQqYhP46e/FjOlNRSvkf94aWOWlWuGyerw0t6xGdqSilnHX6OGxZYF17OZwOEY2h921WwLTq6XR1yqYPP3qgoaKUnzIGslOtcMl4F8pPa0NLP6Kh4oGGilIBoCgf0uZZAXN0J0TFQt+7IGUCxCU5XV1Q0lDxQENFqQBiDOz52q2hZZk2tHSIXqhXSgU+Eeh4mfU6mWt9zsv6afD2BG1o6ad0pqKUCiwul1tDy4+0oaWP6ExFKVU/hYRA0jDrVZgNG2bBBm1o6S90pqKUCnwV5W4NLZdpQ8s6oDMVpVTwCA2D7tdZr28bWr6pDS0doDMVpVT99G1Dy2mw92sICYfk0drQ8jzpTEUpFdzcG1rmbbfCJW2OW0PLCdBnrDa09DKdqSilgkdpEWS8Z117OZAKYVFuDS0H6OzlLHSmopRSZ4poAP3usl7fa2g5G1r1cmto2djpSgOWzlSUUsHNY0PLCVbQKEDbtHikoaKUqlKVDS0vtBta3hT0DS01VDzQUFFKVavKhpZ3Wrcmx3dxujpHaKh4oKGilKoxjw0tJ0C364OqoaVeqFdKqdry2NByot3Q8m4YMF4bWrrRmYpSSp2LKhtaDrUbWo6otw0tdaailFJ14XsNLQ/Ahpl2Q8s7oXEbqxV//3uDtqGlzlSUUqq2PDa0nACdhtSLhpY6U1FKKV+ptqHleOh7NzSKd7rSOqczFaWUqgtVNrS8wW5oeUnAtYTRmYpSSjnJY0PLdyCuixUu9bChpc5UlFLKVwK4oaXOVJRSyt9U1dByy4J61dBSZypKKeWkkhNWp+TUaXB4C0Q0shtaTvSrhpbapsUDDRWllF/y84aWGioeaKgopfyeHza01FDxQENFKRUw/KihZU1DpVaPeYrIiyKyTUQ2i8h7IhJrj18tIutFZIv9dYjbOgPs8SwReVXEut1BRJqJyBIR2Wl/bWqPi71clr2f/m7bGmcvv1NExtXmWJRSyu9UNrQcMw0ez4Shk6Fgr9XQ8uVkWPocHNvjdJXfU9veAUuAnsaY3sAO4Gl7/AhwvTGmFzAOmOW2zr+ASUCS/Rppjz8FLDPGJAHL7J8BrnFbdpK9PiLSDJgMDAIGApMrg0gppeqdRi3gssfhJ2lw9zvQbhCs+Bv8rS+8eQtsW2y1i3FYrULFGPOpMabyKFYDbe3xjcaYg/Z4BhAlIpEi0hqIMcasMtZ5t5nAjfZyo4EZ9vczzhifaSyrgVh7OyOAJcaYfGPMMayAqwwopZSqn0JCoPMwGDsbfpoOV/wSDmdYDS1f6QWf/wmOH6x+O3VVnhe3NRH4qIrxW4CNxpgSIAHIdnsv2x4DaGmMyQGwv7awxxOA/VWs42n8B0RkkoikikhqXl7eOR2UUkr5rSYJcNXTVrjcPhtadIfP/wgv94R5d0HWUqtVvw9V+/CjiCwFWlXx1rPGmIX2Ms8C5cDsM9btAbwADK8cqmI71d0p4GmdGm/LGDMFmALWhfpq9qeUUoHlew0tv7Fa8W+YZTW0jO1gXdj3UUPLamcqxphhxpieVbwqA2UccB1wl3G7lUxE2gLvAfcaY3bZw9nYp8hsbYHKedph+7QW9tdct3XaVbGOp3GllApezTrCsOesC/u3ToXY9tYF/Ze6wyfP1vnua3v310jgl8ANxpgit/FYYDHwtDFmReW4fVrrhIgMtu/6uhdYaL+9COuiPvZX9/F77bvABgOF9nY+AYaLSFP7Av1we0wppVRYpNVXbPwH8PBaGPgANGlX/Xq1VKvnVEQkC4gEjtpDq40xD4rIr7DuBNvptvhwY0yuiKQA04ForGswjxpjjIg0B+YD7YF9wBhjTL4dPv/AughfBEwwxqTa+58IPGNv//fGmGnV1azPqSil1LnThx890FBRSqlz55OHH5VSSil3GipKKaW8RkNFKaWU12ioKKWU8hoNFaWUUl6joaKUUsprNFSUUkp5TdA9pyIiecDeWmwiDqu1f7AItuMFPeZgocd8bjoYY6ptHhZ0oVJbIpJakweA6otgO17QYw4Wesx1Q09/KaWU8hoNFaWUUl6joXLupjhdgI8F2/GCHnOw0GOuA3pNRSmllNfoTEUppZTXaKhUQURGish2EckSkaeqeD9SRN6y318jIom+r9K7anDMj4tIpohsFpFlItLBiTq9qbpjdlvuVhEx9mcBBbSaHLOI3Gb/WWeIyBxf1+htNfi73V5EPhORjfbf71FO1OktIjJVRHJFJN3D+yIir9r/PTaLSH+vFmCM0ZfbCwgFdgGdgAggDUg+Y5kfA6/b348F3nK6bh8c81VAA/v7h4LhmO3lGgNfAquBFKfr9sGfcxKwEWhq/9zC6bp9cMxTgIfs75OBPU7XXctjvhzoD6R7eH8U1gckCjAYWOPN/etM5YcGAlnGmN3GmFJgHjD6jGVGAzPs798GhtqfUBmoqj1mY8xn5ruPjF4NtPVxjd5Wkz9ngOeBPwOnfVlcHanJMT8AvGaMOQZgjMn1cY3eVpNjNkCM/X0T4KAP6/M6Y8yXQP5ZFhkNzDSW1UCsiLT21v41VH4oAdjv9nO2PVblMsaYcqAQaO6T6upGTY7Z3X1Yv+kEsmqPWUT6Ae2MMR/4srA6VJM/5y5AFxFZISKrRWSkz6qrGzU55ueAu0UkG/gQeNQ3pTnmXP9/Pydh3tpQPVLVjOPMW+RqskwgqfHxiMjdQApwRZ1WVPfOeswiEgK8DIz3VUE+UJM/5zCsU2BXYs1GvxKRnsaYgjqura7U5JjvAKYbY/4qIhcBs+xjdtV9eY6o03+/dKbyQ9lAO7ef2/LD6fC3y4hIGNaU+WzTTX9Xk2NGRIYBzwI3GGNKfFRbXanumBsDPYHPRWQP1rnnRQF+sb6mf7cXGmPKjDHfANuxQiZQ1eSY7wPmAxhjVgFRWD2y6qsa/f9+vjRUfmgdkCQiHUUkAutC/KIzllkEjLO/vxVYbuwrYAGq2mO2TwX9GytQAv08O1RzzMaYQmNMnDEm0RiTiHUd6QZjTKoz5XpFTf5u/w/rpgxEJA7rdNhun1bpXTU55n3AUAAR6Y4VKnk+rdK3FgH32neBDQYKjTE53tq4nv46gzGmXEQeAT7BunNkqjEmQ0T+D0g1xiwC3sCaImdhzVDGOldx7dXwmF8EGgEL7HsS9hljbnCs6Fqq4THXKzU85k+A4SKSCVQAPzfGHHWu6tqp4TE/AfxHRH6GdRpofCD/kigic7FOX8bZ14kmA+EAxpjXsa4bjQKygCJgglf3H8D/7ZRSSvkZPf2llFLKazRUlFJKeY2GilJKKa/RUFFKKeU1GipKKaW8RkNFKaWU12ioKKWU8hoNFaWUUl7z/5xsxVZgYGmRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from autograd import grad\n",
    "from pycalphad import Database, Model, variables as v\n",
    "from pycalphad.tests.datasets import ALFE_TDB\n",
    "from pycalphad.codegen.sympydiff_utils import build_functions\n",
    "\n",
    "dbf = Database(ALFE_TDB)\n",
    "mod = Model(dbf, ['AL', 'FE'], 'LIQUID')\n",
    "state_variables = [v.N, v.P, v.T]\n",
    "site_fracs = mod.site_fractions\n",
    "build_output = build_functions(mod.GM, tuple(state_variables + site_fracs),\n",
    "                               include_grad=True, include_hess=False)\n",
    "cf, gf, hf = build_output.func, build_output.grad, build_output.hess\n",
    "\n",
    "def toy_objective(x):\n",
    "    # N,P,T,y1,y2,np,mu_a,mu_b,tau_1,tau_2\n",
    "    inp = np.ascontiguousarray(x[:, :5])\n",
    "    obj = cf(*x[:, :5])\n",
    "    lagrange = x[:, 6]*(0.3 - x[:, 5]*x[:, 3]) + x[:, 7]*(0.7 - x[:, 5]*x[:, 4]) + \\\n",
    "               x[:, 8]*(2500 - x[:, 2]) + x[:, 9]*(1 - (x[:, 3]+x[:, 4]))\n",
    "    return x[:,5] * obj + lagrange\n",
    "\n",
    "def toy_gradient(x):\n",
    "    grad = gf(*x[:5])\n",
    "    obj = cf(*x[:5])\n",
    "    out = np.zeros(x.shape[0])\n",
    "    out[:5] = [x[5]*g for g in grad]\n",
    "    out[5] = obj\n",
    "    out[6] = 0.3 - x[5]*x[3]\n",
    "    out[5] += -x[6]*x[3]\n",
    "    out[3] += -x[6]*x[5]\n",
    "    \n",
    "    out[7] = 0.7 - x[5]*x[4]\n",
    "    out[5] += -x[7]*x[4]\n",
    "    out[4] += -x[7]*x[5]\n",
    "    \n",
    "    out[8] = 2500 - x[2]\n",
    "    out[2] += -x[8]\n",
    "\n",
    "    out[9] = 1 - (x[3]+x[4])\n",
    "    out[3] += -x[9]\n",
    "    out[4] += -x[9]\n",
    "    return out\n",
    "\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "delta = 1e-3\n",
    "result = fsolve(toy_gradient, [1.00000000e+00, 1.01325000e+05, 2.50000000e+03, 7.01402806e-01,\n",
    "       2.98597194e-01, 3.00000000e-01, 0.00000000e+00, 0.00000000e+00, 0, 0])\n",
    "print(result)\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = np.linspace([1,1e5, 2500, 1e-4, 1-1e-4, 1, 0, 0, 0, 0], [1, 1e5, 2500, 1-1e-4, 1e-4, 1, 0, 0, 0,0], 100)\n",
    "y = toy_objective(x)\n",
    "plt.plot(x[:, 3], y)\n",
    "plt.plot([0, 1], [result[7], result[6]])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('fixed_total_moles', array([0.3, 0.7]))\n"
     ]
    }
   ],
   "source": [
    "from pycalphad.core.problem import Problem\n",
    "from pycalphad.core.composition_set import CompositionSet\n",
    "from pycalphad.codegen.callables import build_phase_records\n",
    "\n",
    "conds = {v.N: 1, v.P: 1e5, v.T: 2500, v.X('FE'): 0.7}\n",
    "comps = [v.Species('AL'), v.Species('FE')]\n",
    "\n",
    "prx = build_phase_records(dbf, comps, ['LIQUID'], conds,\n",
    "                         models={'LIQUID': mod}, build_gradients=True)\n",
    "prx = prx['LIQUID']\n",
    "compset = CompositionSet(prx)\n",
    "prob = Problem([compset], comps, {str(key): value for key, value in conds.items()}, [0, 0])\n",
    "\n",
    "def pycalphad_obj(x):\n",
    "    obj = prob.objective(x[:8])\n",
    "    c = prob.cl[2:] - prob.constraints(x[:8])[2:]\n",
    "    lagrange = np.sum(-c * x[8:])\n",
    "    return obj+lagrange\n",
    "\n",
    "def pycalphad_grad(x):\n",
    "    out = np.zeros(x.shape[0])\n",
    "    res = prob.gradient(x[:8])\n",
    "    c = prob.cl[2:] - prob.constraints(x[:8])[2:]\n",
    "    out[:8] = res\n",
    "    out[:8] += -np.dot(prob.jacobian(x[:8])[2:].T, x[8:])\n",
    "    out[8:] = -c\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00000000e+00,  0.00000000e+00, -3.18868759e+01,  7.70156224e+03,\n",
       "       -3.80257214e+03, -1.99247053e+05,  8.95791582e-02,  6.10420842e-01,\n",
       "       -0.00000000e+00,  5.55111512e-17])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = [1.00000000e+00, 1.01325000e+05, 2.50000000e+03, 7.01402806e-01,\n",
    "       2.98597194e-01, 3.00000000e-01, 0.00000000e+00, 0.00000000e+00, 0, 0]\n",
    "pycalphad_grad(np.array(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00000000e+00,  0.00000000e+00, -3.18868759e+01,  7.70156224e+03,\n",
       "       -3.80257214e+03, -1.99247053e+05,  8.95791582e-02,  6.10420842e-01,\n",
       "        0.00000000e+00,  0.00000000e+00])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "toy_gradient(np.array(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.00000000e+00  1.01325000e+05  2.50000000e+03  3.00000000e-01\n",
      "  7.00000000e-01  1.00000000e+00 -2.20583626e+05 -1.94101859e+05\n",
      " -1.11878832e+02  2.14055264e+05]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAD8CAYAAAC/1zkdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8VfX9x/HXJzuMECBhhRGQMMKGCLgVEBAHLhQnQ6VatbZqW0db7M8ua6vW1tbSypQhOAqKi+FkhxFIwgrICAQSCAkjIet+f3+cE71iLgnk5p57cz/Px+M+knzvGZ8jyDvfMz5XjDEopZRS3hDidAFKKaXqDw0VpZRSXqOhopRSyms0VJRSSnmNhopSSimv0VBRSinlNRoqSimlvEZDRSmllNdoqCillPKaMKcL8LW4uDiTmJjodBlKKRVQ1q9ff8QYE1/dckEXKomJiaSmpjpdhlJKBRQR2VuT5fT0l1JKKa/RUFFKKeU1GipKKaW8RkNFKaWU12ioKKWU8hoNFaWUUl6joaKUUsprNFSUUioI/OfL3Xy2PbfO96OhopRS9VzeiRJe/GQ7n23TUFFKKVVLc9bso7TCxfiLE+t8XxoqSilVj5WUVzBr9V6GdGtBp/hGdb4/DRWllKrHPkjL4cjJEiZe0tEn+9NQUUqpesoYw9QV39ClZSMu6dzcJ/usVaiIyBgRyRARl4ikuI2Hi8gMEdkiIltF5Gm390aKyHYRyRKRp9zGO4rIGhHZKSJviUiEPR5p/5xlv5/ots7T9vh2ERlRm2NRSqn6Zu03+WQcPM7ESzoiIj7ZZ21nKunAzcCXZ4yPASKNMb2AAcCPRCRRREKB14BrgGTgDhFJttd5AXjZGJMEHAPus8fvA44ZYzoDL9vLYa83FugBjAT+aW9fKaUUMHXFNzRtEM6N/RJ8ts9ahYoxZqsxZntVbwENRSQMiAZKgePAQCDLGLPbGFMKzANGixWhQ4C37fVnADfa34+2f8Z+f6i9/GhgnjGmxBjzDZBlb18ppYJexsFCPs08zN2DOxAV7rvft+vqmsrbwCkgB9gH/MUYkw8kAPvdlsu2x5oDBcaY8jPGcV/Hfr/QXt7Ttn5ARCaJSKqIpObl5dX+6JRSys/9+ePtxESFc/9lnXy632o/+VFElgKtqnjrWWPMQg+rDQQqgDZAU+AreztVndQzZxnnPNf5/qAxU4ApACkpKVUuo5RS9cXKrCN8sSOPZ0d1p0l0uE/3XW2oGGOGncd27wQ+NsaUAbkisgJIwZpZtHNbri1wEDgCxIpImD0bqRwHawbSDsi2T6c1AfLdxs/cllJKBS2Xy/Cnj7eREBvNPRd18Pn+6+r01z5giFgaAoOBbcA6IMm+0ysC60L7ImOMAT4DbrXXHwdUzoIW2T9jv7/cXn4RMNa+O6wjkASsraPjUUqpgPBheg6bswt5/OouPr2WUqm2txTfJCLZwEXAYhH5xH7rNaAR1t1h64BpxpjN9izkEeATYCsw3xiTYa/zS+BxEcnCumbyhj3+BtDcHn8ceArAXm8+kAl8DDxsjKmozfEopVQgKymv4MVPttOtVWOf3vHlTqxf+oNHSkqKSU1NdboMpZTyun9+nsWfP97OjIkDuaJLvFe3LSLrjTEp1S2nT9QrpVQ9cLCgmL8vy2J4ckuvB8q50FBRSql64PeLt+Iyhl9fl1z9wnVIQ0UppQLc1zuPsHhLDg9f1Zl2zRo4WouGilJKBbDScheTF6XTvlkDJl3u2wcdq6KhUkPGGFbuOsLpMr3BTCnlP17/Yhe78k7x3A3JjtxCfCYNlRpav/cYd/5nDe+n6fOVSin/kJV7gn8sz+L6Pm0Y0q2l0+UAGio1NqBDU7q0bMSMVXsIttuwlVL+x+UyPPXOFhpEhjL5emcvzrvTUKkhEeHeixJJP3CcDfsKnC5HKRXkZq/dR+reY/zq2mTiGkU6Xc63NFTOwU39EmgcFcaMlXucLkUpFcQOFhTzwkfbuCwpjlv6O/PkvCcaKuegYWQYYwa048MtOeQeP+10OUqpIORyGX7x9mZcxvCHm3r57BMda0pD5Rzde1EHKoxh9pp9TpeilApCs9fs5eusI/zq2mTHn0mpiobKOUqMa8iVXeKZs3YfpeUup8tRSgWRPUdO8YcPt3F5l3juGNiu+hUcoKFyHsZdnEjeiRK9vVgp5TMVLsOTC9IIDxX+fEtvvzvtVUlD5TxcnhRPcusY/rZsp85WlFI+8foXu0jde4zfju5BqyZRTpfjkYbKeQgJEX4+oiv78ouYn7rf6XKUUvXcpv0FvLxkB9f1bs2Nff3rbq8zaaicpyu7xpPSoSmvLtuprVuUUnXmZEk5j83bSMuYKH7vh3d7nUlD5TyJCL8Y2Y3cEyX63IpSqs48tyiD/flFvHx7X5pEhztdTrU0VGphYMdmXNElnn99sYvjp8ucLkcpVc8s3HSAt9dn8/BVnRnYsZnT5dSIhkot/XxEVwqKynj9811Ol6KUqkd2553kmXe3MKBDU34yNMnpcmpMQ6WWeiY04aZ+Cfz3q2/Yd7TI6XKUUvXA6bIKHp6zkfCwEP5+Rz/CQwPnn+rAqdSPPXVNN8JChd8tznS6FKVUPfD8B5lszTnOS7f1oU1stNPlnBMNFS9oGRPFw1d15tPMw3y984jT5SilAtjCTQeYvWYfP7qik998Rsq50FDxkvsu7Uj7Zg347fsZlFXoA5FKqXO3/dAJnnpnCxcmNuXJ4V2dLue8aKh4SVR4KM9e252duSeZtWqv0+UopQLM8dNlPPjmehpFhfHanf0D6jqKu8Cs2k8NT27J5V3i+eun28kpLHa6HKVUgHC5DE/MT2N/fhH/vKs/LWL8tw1LdTRUvEhE+N3onlQYw+SFGU6Xo5QKEK99lsWSzMM8M6o7FyYGxvMonmioeFn75g14bGgXPs08zCcZh5wuRynl55ZkHuavS3ZwY982TLgk0elyak1DpQ7cf1lHurVqzOSFGZwsKXe6HKWUn9p5+AQ/e2sTvRKa8Cc/bmd/LjRU6kB4aAh/vLkXh0+c5s8fb3O6HKWUHyosKuOBmalEhYcy5d4BRIWHOl2SV2io1JF+7Zsy7qJEZq7ay6pdR50uRynlR8oqXDw8ZwMHCor59z39ad0ksB5wPBsNlTr0i5FdSWzegF+8k8YpPQ2mlAKMMUxelMHXWUf4w029GNAhsC/Mn0lDpQ41iAjjxTF9yD5WzJ8+0tNgSimYtmIPc9bs46ErL2BMin9+znxtaKjUsQsTmzHxko7MWr2XFVnawkWpYLZ822F+tziTkT1a8fMAfWK+OhoqPvDzEV3pFNeQny9Io7BIP3dFqWC0JbuQR+ZspEebJrx0ex9CQgL/Tq+qaKj4QFR4KC/f3pfcEyU8878tGGOcLkkp5UP784uYMH0dTRtE8Mb4FBpEhDldUp3RUPGRPu1i+dnVXVi8OYd3NxxwuhyllI8UFpUxftpaSssrmD7hQlo0DtwWLDWhoeJDD15xAQM7NuM3C9PZe/SU0+UoperY6bIKHpiZyv78Yqbcm0JSy8ZOl1TnahUqIjJGRDJExCUiKW7jESIyTUS2iEiaiFzp9t4AezxLRF4V+xFSEWkmIktEZKf9tak9LvZyWSKyWUT6u21rnL38ThEZV5tj8YXQEOHl2/sSEiI8Nm8TpeXaIl+p+qq8wsVP5m5k3d58/nJbHwZ3au50ST5R25lKOnAz8OUZ4w8AGGN6AVcDfxWRyn39C5gEJNmvkfb4U8AyY0wSsMz+GeAat2Un2esjIs2AycAgYCAwuTKI/FlCbDQv3NKbTfsLePETvc1YqfrIGMOvF2bwaeZhJl+XzA192jhdks/UKlSMMVuNMdureCsZKxgwxuQCBUCKiLQGYowxq4x1tXomcKO9zmhghv39jDPGZxrLaiDW3s4IYIkxJt8YcwxYwncB5ddG9WrNuIs68J+vvmFJ5mGny1FKednLS3Ywd+0+Hr7qAsZf0tHpcnyqrq6ppAGjRSRMRDoCA4B2QAKQ7bZctj0G0NIYkwNgf21hjycA+6tYx9N4QHjm2u70TIjhyQVpZB8rcrocpZSX/Per3by6PIvbU9oF7Kc31ka1oSIiS0UkvYrX6LOsNhXrH/lU4BVgJVAOVHVjdnX313pap8bbEpFJIpIqIql5eXnV7M43IsNCee3O/rhchofnbKSkvMLpkpRStfTWun38bvFWru3Vmj/c3KtedB0+V9WGijFmmDGmZxWvhWdZp9wY8zNjTF9jzGggFtiJFTRt3RZtCxy0vz9sn9bC/pprj2djzXLOXMfTeFX1TDHGpBhjUuLj46s7ZJ/p0LwhL47pTdr+An77fqbT5SilamHx5hyefncLl3eJ5+Xb+xJaTx9urE6dnP4SkQYi0tD+/mqg3BiTaZ/WOiEig+27vu4FKsNpEVB5B9e4M8bvte8CGwwU2tv5BBguIk3tC/TD7bGAMrJnax668gLmrNnH/HX7q19BKeV3lmYe5rF5G+nfvimv392fiLDgfVqjVo91ishNwN+BeGCxiGwyxozAuh7yiYi4gAPAPW6rPQRMB6KBj+wXwJ+A+SJyH7APGGOPfwiMArKAImACgDEmX0SeB9bZy/2fMSa/NsfjlCeHd2VLdiG/WphOt9aN6d021umSlFI19MWOPH48ewM92sQwdcKF9fpp+ZqQYGsZkpKSYlJTU50u4wfyT5Vy/d+/xhjDwkcuJb5xpNMlKaWqsXLXESZMW8cF8Y2Y+8BgmjQId7qkOiMi640xKdUtF7xzND/TrGEE/75nAMeKyvjRrFS9cK+Un1u16yj3TU+lQ/MGvHn/oHodKOdCQ8WP9Exowl9v68OGfQU88266Np5Uyk+t2nWUidPX0bZpNLPvH0yzhhFOl+Q3NFT8zKherfnpsCTe2ZDNf77a7XQ5SqkzuAfKnAcG66nqMwT3FSU/9ZMhSew8fJI/frSNDs0bMqJHK6dLUkoBX+88wv0z19GuaQMNFA90puKHQkKEv4zpQ5+2sTw2byOb9hc4XZJSQW/5tsNMnLGOxOYNmTtJA8UTDRU/FR0Ryn/HpRDfOJL7Z6xjf762clHKKR+n5/CjWevp2rIx8yYNJq6RBoonGip+LK5RJNPGD6SswjB+2lqOnSp1uiSlgs4767N5eM5GereNZfYDg4htoBflz0ZDxc91btGIKfcMYP+xYibOWEdRabnTJSkVNKat+IYnFqQxuFMzZk4cSEyU3jZcHQ2VADCoU3NeHduXtP0FPDJnI2UV+uFeStUlYwyvLN3Bb9/PZESPlkwdfyENI/W+pprQUAkQI3u25nc39mL5tlyeemcLLpc+w6JUXahwGSYvyuCVpTu5pX9bXruzP5FhoU6XFTA0egPInYPac+RkCS8t2UHjqDAmX58clK21laorJeUVPP5WGou35DDp8k48NbIbIUHabfh8aagEmEeHdKawuIw3vv6GRpFhPDki+D4ESKm6cPx0GQ/OWs/KXUd5dlR3Hri8k9MlBSQNlQAjIvzq2u4UlZbzj8+yaBgZxkNXXuB0WUoFtJzCYiZMW0dW7kleuq0PN/dvW/1KqkoaKjVVdhq+eAH63Q3Nnf1HXET43Y29KCqt4IWPtxERFsJ9lwbX52Ar5S2ZB48zYfpaTpVUMH3CQC5NinO6pICmoVJTB9bDylfh65eg4xWQMhG6XQuhztxiGGo/dV9a7uL5DzIJEZhwiQaLUufiix15PDx7A40iw1jw4EV0bx3jdEkBT+/+qqnES+BnGTDkV5C/GxaMg5d7wLLnoWCfIyWFh4bw6h39GNGjJb99P5MZK/c4UodSgWjmqj1MmLaWds0a8N7DF2ugeIl+SNf5cFVA1lJInQo77E8wThpuzV6SroYQ395+WFru4pE5G/g08zC/uS6ZiXoqTCmPyitc/G7xVqav3MPQbi149Y5++gxKDdT0Q7o0VGqrYB9smGm9Th6GJu2g/zjofw809l134dJyF4/N28hH6Yf4xciu/PjKzj7bt1KBorCojEfmbuCrnUe479KOPDOqO6F6y3CNaKh4UGcfJ1xRBts/tGYvuz+HkDDoOsqavXS8AkLq/kxjeYWLJxaksXDTQX4ypDM/u7qLPseilC0r9yQPzEwl+1gRz4/uydiB7Z0uKaDUNFR0zuctoeGQPNp6Hd0F66fBxtmwdRE06wQDJkDfu6Bh8zorISw0hJdu60tUWCivLs/iREk5v742WR/eUkFv+bbDPDZ3ExFhIcx5YDAXJjZzuqR6S2cqdanstBUqqVNh3yoIjYDkG63ZS/vBUEezCJfL8PziTKat2MNN/RL48629CQ/VezJU8HG5DH9fnsUry3aQ3DqGKfemkBAb7XRZAUlnKv4gPAp632a9Dmdas5e0ebBlPsR3t8Klz+0Q1cSruw0JEX5zXTLNG0bwl093UFhcxmt39ic6QvsXqeBx/HQZj7+VxtKth7m5XwJ/uLkXUeH6/0Bd05mKr5WegvR3rNnLwY0Q3gB63mIFTEJ/r+/uzdV7+fXCdPq2i+WNcRfSrKF+FoSq/zIOFvLj2Rs4cKyYX13bnXEXJ+r1xVrSC/UeOB4q7g5ssGYvW96GsiJo3dcKl163QkRDr+3m4/QcHpu3iTax0UyfcCEdmntv20r5m/nr9vPrhenENgjntTv7k6LXT7xCQ8UDvwqVSqcLYfN8a/aSmwmRMdD7dkiZAC17eGUX6/fmc/+MVEJE+O+4FPq1b+qV7SrlL06VlPObhRm8syGbSzo3529j++nH/nqRhooHfhkqlYyB/WutcMl4DypKoN1ga/aSPNq6RlMLu/NOMn7aOg4dP81fxvThhj5tvFS4Us7amnOcR+ZsYPeRUzw6JInHhibp8ydepqHigV+HiruifNg0G1KnQf4uiG5q3ZI8YALEnf+DjfmnSnlw1nrW7snnsaFJ/HRYkp5rVgHLGMOba/bx/AeZNIkO52+39+XiztoQsi5oqHgQMKFSyeWCPV9as5dti8FVXuuGliXlFTzzbjrvbMjm2l6teXFMbxpE6I2AKrAcPVnCL9/ZwtKth7m8Szwv3dZHT3fVIQ0VDwIuVNydOAwbZ8H6GVC4Dxq1hH73wIBxEHtuTwcbY5jy5W5e+HgbXVo2Zso9KbRv3qCOClfKu77ckccTC9IoLCrjqWu6Mf7iRH3It45pqHgQ0KFSyVUBWcus2cvOT6xrMefZ0PKLHXk8OmcDIsLf7+jH5V3i67BwpWqnqLScP364jVmr99KlZSP+Nrafdhf2EQ0VD+pFqLgr2O/W0PIQxLSFAePPqaHlniOn+NGs9ezIPcFPh3bh0SGd9bc+5Xc27DvGE/PT2HP0FPdd0pEnR3TVhxl9SEPFg3oXKpUqymD7R3ZDy8/OuaFlUWk5z76XznsbD3B5l3heub2vPiip/MLpsgpeWrKD/361m9ZNovnLmD5cdEHd9dBTVdNQ8aDehoq7o7tg/XTY+CYU59e4oaUxhrlr9/PcogyaNYzglbF9GdxJ/+dVzlm/N5+fL9jM7iOnuGNge54Z1Y3GUc582mqw01DxIChCpVLZadj6vt3QcmWNG1qmHyjk0bkb2XvUuuf/0SGdCdOGlMqHTpaU8+LH25i5ei9tmkTzwi299bPjHaah4kFQhYq73K3WMy9pc6HkeLUNLU+WlPObhem8u+EAKR2a8vLtfWnXTO8OU3Vv2dbD/Op/6Rw6fppxFyXy5IiuNNJPZnSchooHQRsqlUpPQfq7dkPLDdU2tHxvYza/+V8GLmOYfEMPxgxoqw9LqjpxoKCY3y7K4NPMw3Rt2Zg/3tKL/tpOyG9oqHgQ9KHi7uBGa/ayZcH3G1r2vAUiG327WPaxIp6Yn8aab/K5Orklv7+pJy0a165ljFKVSstdTF/5Da8s3YnLGH4yNIn7L+1ERJiecvUnGioeaKhUoQYNLV0uwxtff8OLn24nOjyU525I5sa+CTprUbXy5Y48nns/g915pxjWvQWTr++hp1n9VE1DpVa/CojIiyKyTUQ2i8h7IhLr9t7TIpIlIttFZITb+Eh7LEtEnnIb7ygia0Rkp4i8JSIR9nik/XOW/X5idftQ5yiqCQx8AB5aCRM/tW5F3jAT/nUxvDEC0t4ipKKEBy7vxIc/uYwL4hvys7fSuH9GKgcKip2uXgWgPUdOMWlmKvdOXYvLZZg6PoX/jrtQA6UeqNVMRUSGA8uNMeUi8gKAMeaXIpIMzAUGAm2ApUAXe7UdwNVANrAOuMMYkyki84F3jTHzROR1IM0Y8y8R+THQ2xjzoIiMBW4yxtzuaR/GmIqz1awzlRoqyodNc6zZyxkNLSuaXcC0Fd/w1093IAJPDO/K+IsTtSusqlZhcRn/WL6T6Sv3EB4awsNXdeb+yzoSGaYPMfo7n5/+EpGbgFuNMXeJyNMAxpg/2u99AjxnL/qcMWaEPf60PfYnIA9oZQfURZXLVa5rjFklImHAISAeeKqqfRhjVp2tTg2Vc2QMfFPZ0PIDt4aWE9jf4ip+/cEOPt+eR8+EGJ4f3VM/p0VVqaS8gtmr9/H35TspKC5jzIC2PDm8Ky1i9NpcoHDiM+onAm/Z3ycAq93ey7bHAPafMT4IaA4UGGPKq1g+oXIdO3AK7eXPtg/lLSLQ6Qrr5d7QcsF42jVswbR+97C82zU881khN/1zJbeltOWXI7vRXLvFKqxrce9vPshfPt3O/vxiLuncnKev6U7PhB/exq7qh2pDRUSWAlU1kXrWGLPQXuZZoByYXblaFcsbqr6GY86y/Nm2dbZ1vkdEJgGTANq3P7duvspN45Zw+ZNw6c++bWgpK15mqHmJFRcM472QETy7wfBx+iF+MjSJey9K1Dt4gpQxhiWZh3lpyQ62HTpBcusYZk7sxWVJcXpzRz1XbagYY4ad7X0RGQdcBww1351LywbauS3WFjhof1/V+BEgVkTC7NmK+/KV28q2T381AfKr2ceZxzAFmALW6a+zHY+qgZBQ6DLcetkNLcM2zGTMySXc2KwNi0KH86fFFzJr9V6evqYbI3q00n9IgoQxhs935PHK0p2k7S+gY1xD/ja2L9f3bqNNSoNEbS/UjwReAq4wxuS5jfcA5vDdRfRlQBLW7GIHMBQ4gHWh/k5jTIaILADecbtQv9kY808ReRjo5Xah/mZjzG2e9qEX6h1yRkNLI6F8HTqQfxddSXHCJfzimmQGaR+xessYw/Jtuby6bCdp2YUkxEbz2NAkbu6foC1+6gmfXKgXkSwgEjhqD602xjxov/cs1nWWcuCnxpiP7PFRwCtAKDDVGPN7e7wTMA9oBmwE7jbGlIhIFDAL6Ic1QxlrjNl9tn2cjYaKD9gNLc3GN5HifPbTilllV5HT8WYeGDmQ3m1jq9+GCgjlFS4Wb8nhX5/vYtuhE7RrFs0jV3Xmpn5t9dRnPaMPP3qgoeJDdkNL17o3CNm/ilLC+LBiINsSxjDq2pvo3U7vFAtUp0rKeXt9Nv/9ejf784vp3KIRD15xAaP7tiFcZyb1koaKBxoqDsndSumaNzCb5hJZcZIdrgTWNL+RbiMeIKVrol5zCRA5hcXMWLmXOWv2cvx0Of3ax/LQFRcwrHtLvWZSz2moeKCh4rDSUxRvXMDxr6fQ8kQGxSaCFdFX0uDiBxh06dX6AKUfMsawenc+M1ft4dPMwxhjGNmzFfdd2okBHXS2GSw0VDzQUPEfJfvWs/fT12iX/QHRlLBdOpHX7S76XHMfjWP0HyunFRSV8u6GA8xdu4+duSdpEh3O2AvbcffgDtpOJQhpqHigoeJ/yk8dY/uSqTRKn0mH8j2cNNGkx40k/qqHuKDnIKfLCyoVLsPKXUd4e302H6UforTcRZ92sdw1sD3X92lDdIS2UwlWGioeaKj4MWPI2rCMwi+n0LNgOZFSxtaw7hxLvpuew8cR06ix0xXWW9sPnWDhpgO8t/EAOYWniYkKY3TfBMYObEePNvr0u9JQ8UhDJTAUHj3E1o//Tdtd82jrOsgx04iNzUYRfdF9pAwYqHcYecHuvJN8lH6I99MOsu3QCUJDhMuS4rh1QFuGdW9JVLjOStR3NFQ80FAJLMblYve6jyha+R+6FX5JOBWspSdZHW6jw0VjGJTUSh+uqyFjDJk5x1mamctH6TlsO3QCgP7tY7mxXwKjerUmTnu2KQ80VDzQUAlcpQU57Fv6b5pum0Pz8sPkmSYslCEcThpL/959uKxLvH6W+RlOlZSzatdRPt+Ry/KtuRwsPI0IDGjflFG9WjOyZyvaxEY7XaYKABoqHmio1AOuCkq3L6Hgy9eJy/kCMHxe0Ye3zDBOdRjCJUmtuCwpjuTWMUH37ERZhYvN2QWs2nWUFVlHSd2bT1mFITo8lEuT4rg6uSVDurXQGYk6ZxoqHmio1DMF+3Gtn0F56gwiinPJlThmlV7JWxVXUdagBRcmNmNQp+ZcmNiU7q1j6t21mOOny0jbX0DqnmOk7s1n474Cikqt9nfdW8dweVIcV3SJZ0BiU/0gLFUrGioeaKjUU2c0tHRJKOmNLmF66RDeK+yMIYTIsBB6JTShb7tYeiTEkNy6CRfENwyYazL5p0rZmnOcrTnHyTh4nLTsAnbnnQIgRKBbqxhSEptyUafmDOrUnGYNIxyuWNUnGioeaKgEAbuhJRvfhOJ8ypsksqPtrXwUNoQVOZBx8Dgl5S4AIsJC6BTXkKSWjekc34jEuAa0b9aADs0b0rRBuM/bx5wqKedAQTF7jxax9+gp9hw9RVbuSbJyT3LkZOm3y7VoHEnvtrH0bdfE+to+lpiocJ/WqoKLhooHGipBpLwEMhdZs5d9KyE0ApJHU95/PLuje5ORc5ytOSfIyj3JztwT7M8v/t7qUeEhtG4STesmUcQ3jqR5w0iaN4ogJjqcmKgwYqLCiQoPJSo8hKjwUMJDhRCxXgaocLkodxlKy10Ul1ZQXFbByZJyCovLKCwuo6CojCMnSzhyspTc46fJKTxNYXHZ92poHBVGUotGdG7RiKQWjeneOoburRvrJ2sqn9NQ8UBDJUjlboXUaZA2D0oKIb4bpEyE3rdDtNWK/3RZBfvzi6xZQn4ROQXF5Bw/TU5BMUdPlXL0ZCknS8qr2VHNRYSFEN8DoG2CAAAQOElEQVQokrhGEcQ3jrQCLDaKhNhoOjRvSIdmDYh1YLakVFU0VDzQUAlypUWQ8a41ezmwHsKiodctVsC06Q/V/AN+uqyC46fLOF5czonTZRSXVVBS5uJ0WQVlLoMxhgqXIUSE0BDrFRkWQnR4KFERoTSKDKNJdLg9ywnRwFABQ0PFAw0V9a2Dm6xw2fI2lJ2C1n2scOl5K0Q2cro6pfyKhooHGirqB04Xwub51umx3AyIaAx9bocBE6BVT6erU8ovaKh4oKGiPDIG9q+F9dMg/V2oKIF2g6zZS/JoCNcnz1Xw0lDxQENF1UhRPqTNtU6PHc2C6KbQ9y5r9hLX2enqlPI5DRUPNFTUOTEG9nxlhcvW98FVDh0vt2YvXa+FMH3AUAWHmoaKdt9T6mxErBDpeDmcOAyb3oTU6bBgPDRsAf3vgf7joGkHpytVyi/oTEWpc+WqgF3LrdnLjo+t2UzS1dbsJWk4hGiPLVX/6ExFqboSEmqFSNLVUJgNG2bC+hkwdyzEtIUB46DfPRDT2ulKlfI5nako5Q0VZdasJXWqNYuRUOg2ypq9dLwSQgKjaaVSnuhMRSlfCg2H7tdbr6O7YMMMq6Hl1vehaUdImWDdPdYwzulKlapTOlNRqq6Ul1ihkjoV9q74tqElKROh/UXVtoRRyp/oTEUpp4VFQq9brVfuNuuhyk1zYcuCKhtaKlUf6ExFKV+qZUNLpZyiMxWl/FFEA+h3t/U6uMmavWxeYF1/0YaWqh7QmYpSTjt9HLbMh3VTtaGl8lvapsUDDRXlt4yB7HXWqbHKhpZtB1qzlx43akNL5SgNFQ80VFRAOLOhZVSsdUtyygSIS3K6OhWENFQ80FBRAUUbWio/oRfqlaoPztrQMt5qBzNgHDRNdLhQpSw6U1Eq0FTV0LLzsO8aWobq74rK+3SmolR95amh5bw7ICbBasXf/x6IaeN0pSoI6UxFqfqgqoaWXa+xZi+drtKGlqrWdKaiVDBxb2iZvxvWT7ceqNz2gXW9ZcAE64FLbWip6litfn0RkRdFZJuIbBaR90Qk1h5vLiKfichJEfnHGesMEJEtIpIlIq+KWH0pRKSZiCwRkZ3216b2uNjLZdn76e+2rXH28jtFZFxtjkWpeqNZJ7j6/+DxrXDLG9YpsaWT4aXu8PZ9sGeFdR1GqTpQ2znxEqCnMaY3sAN42h4/DfwaeLKKdf4FTAKS7NdIe/wpYJkxJglYZv8McI3bspPs9RGRZsBkYBAwEJhcGURKKb5raDnhQ/jxGutU2M4lMH0UvDYIVr8OxQVOV6nqmVqFijHmU2NMuf3jaqCtPX7KGPM1Vrh8S0RaAzHGmFXGupgzE7jRfns0MMP+fsYZ4zONZTUQa29nBLDEGJNvjDmGFXCVAaWUcteiG1zzAjyxDUa/ZvUW+/iX8Ndu8L+HIXu9zl6UV3jzmspE4K1qlkkAst1+zrbHAFoaY3IAjDE5ItLCbZ39VazjaVwp5Ymnhpab3oRWva3ZTK8x2tBSnbdqZyoislRE0qt4jXZb5lmgHJhd3eaqGKvu1yNP69R4WyIySURSRSQ1Ly+vmt0pFSTa9IXr/2bNXq79q/X8ywc/tWYvHzwOh9KdrlAFoGpnKsaYYWd7375Afh0w1FR/f3I29ikyW1vgoP39YRFpbc9SWgO5buu0q2KdbODKM8Y/93AMU4ApYN1SXE2NSgWXqBi48H5Iue+7hpYb34TUN7ShpTpntb37ayTwS+AGY0xRdcvbp7dOiMhg+66ve4GF9tuLgMo7uMadMX6vfRfYYKDQ3s4nwHARaWpfoB9ujymlzocItBsIN71uzV5G/AGK8+F/D1qzl4+fgSM7na5S+blaPfwoIllAJHDUHlptjHnQfm8PEANEAAXAcGNMpoikANOBaOAj4FFjjBGR5sB8oD2wDxhjjMm3w+cfWBfhi4AJxphUex8TgWfsff/eGDOtupr14UelzkFVDS0TL7NmL92u04aWQUS7FHugoaLUeTqZCxtnWQ9WFuzThpZBRkPFAw0VpWpJG1oGJW3TopSqG9rQUp2FzlSUUrWnDS3rPZ2pKKV85wcNLWdY11+0oWXQ0ZmKUqpulJdYd4ylToW9KyA0ArrfYM1eOlxs3cKsAobOVJRSzqpsaNnrVsjdZrWE2TQX0t+GuK5WuPS5HaK1D2x9ojMVpZTvlBZBxruQOg0OpEJYNPS8xQqYhP46e/FjOlNRSvkf94aWOWlWuGyerw0t6xGdqSilnHX6OGxZYF17OZwOEY2h921WwLTq6XR1yqYPP3qgoaKUnzIGslOtcMl4F8pPa0NLP6Kh4oGGilIBoCgf0uZZAXN0J0TFQt+7IGUCxCU5XV1Q0lDxQENFqQBiDOz52q2hZZk2tHSIXqhXSgU+Eeh4mfU6mWt9zsv6afD2BG1o6ad0pqKUCiwul1tDy4+0oaWP6ExFKVU/hYRA0jDrVZgNG2bBBm1o6S90pqKUCnwV5W4NLZdpQ8s6oDMVpVTwCA2D7tdZr28bWr6pDS0doDMVpVT99G1Dy2mw92sICYfk0drQ8jzpTEUpFdzcG1rmbbfCJW2OW0PLCdBnrDa09DKdqSilgkdpEWS8Z117OZAKYVFuDS0H6OzlLHSmopRSZ4poAP3usl7fa2g5G1r1cmto2djpSgOWzlSUUsHNY0PLCVbQKEDbtHikoaKUqlKVDS0vtBta3hT0DS01VDzQUFFKVavKhpZ3Wrcmx3dxujpHaKh4oKGilKoxjw0tJ0C364OqoaVeqFdKqdry2NByot3Q8m4YMF4bWrrRmYpSSp2LKhtaDrUbWo6otw0tdaailFJ14XsNLQ/Ahpl2Q8s7oXEbqxV//3uDtqGlzlSUUqq2PDa0nACdhtSLhpY6U1FKKV+ptqHleOh7NzSKd7rSOqczFaWUqgtVNrS8wW5oeUnAtYTRmYpSSjnJY0PLdyCuixUu9bChpc5UlFLKVwK4oaXOVJRSyt9U1dByy4J61dBSZypKKeWkkhNWp+TUaXB4C0Q0shtaTvSrhpbapsUDDRWllF/y84aWGioeaKgopfyeHza01FDxQENFKRUw/KihZU1DpVaPeYrIiyKyTUQ2i8h7IhJrj18tIutFZIv9dYjbOgPs8SwReVXEut1BRJqJyBIR2Wl/bWqPi71clr2f/m7bGmcvv1NExtXmWJRSyu9UNrQcMw0ez4Shk6Fgr9XQ8uVkWPocHNvjdJXfU9veAUuAnsaY3sAO4Gl7/AhwvTGmFzAOmOW2zr+ASUCS/Rppjz8FLDPGJAHL7J8BrnFbdpK9PiLSDJgMDAIGApMrg0gppeqdRi3gssfhJ2lw9zvQbhCs+Bv8rS+8eQtsW2y1i3FYrULFGPOpMabyKFYDbe3xjcaYg/Z4BhAlIpEi0hqIMcasMtZ5t5nAjfZyo4EZ9vczzhifaSyrgVh7OyOAJcaYfGPMMayAqwwopZSqn0JCoPMwGDsbfpoOV/wSDmdYDS1f6QWf/wmOH6x+O3VVnhe3NRH4qIrxW4CNxpgSIAHIdnsv2x4DaGmMyQGwv7awxxOA/VWs42n8B0RkkoikikhqXl7eOR2UUkr5rSYJcNXTVrjcPhtadIfP/wgv94R5d0HWUqtVvw9V+/CjiCwFWlXx1rPGmIX2Ms8C5cDsM9btAbwADK8cqmI71d0p4GmdGm/LGDMFmALWhfpq9qeUUoHlew0tv7Fa8W+YZTW0jO1gXdj3UUPLamcqxphhxpieVbwqA2UccB1wl3G7lUxE2gLvAfcaY3bZw9nYp8hsbYHKedph+7QW9tdct3XaVbGOp3GllApezTrCsOesC/u3ToXY9tYF/Ze6wyfP1vnua3v310jgl8ANxpgit/FYYDHwtDFmReW4fVrrhIgMtu/6uhdYaL+9COuiPvZX9/F77bvABgOF9nY+AYaLSFP7Av1we0wppVRYpNVXbPwH8PBaGPgANGlX/Xq1VKvnVEQkC4gEjtpDq40xD4rIr7DuBNvptvhwY0yuiKQA04ForGswjxpjjIg0B+YD7YF9wBhjTL4dPv/AughfBEwwxqTa+58IPGNv//fGmGnV1azPqSil1LnThx890FBRSqlz55OHH5VSSil3GipKKaW8RkNFKaWU12ioKKWU8hoNFaWUUl6joaKUUsprNFSUUkp5TdA9pyIiecDeWmwiDqu1f7AItuMFPeZgocd8bjoYY6ptHhZ0oVJbIpJakweA6otgO17QYw4Wesx1Q09/KaWU8hoNFaWUUl6joXLupjhdgI8F2/GCHnOw0GOuA3pNRSmllNfoTEUppZTXaKhUQURGish2EckSkaeqeD9SRN6y318jIom+r9K7anDMj4tIpohsFpFlItLBiTq9qbpjdlvuVhEx9mcBBbSaHLOI3Gb/WWeIyBxf1+htNfi73V5EPhORjfbf71FO1OktIjJVRHJFJN3D+yIir9r/PTaLSH+vFmCM0ZfbCwgFdgGdgAggDUg+Y5kfA6/b348F3nK6bh8c81VAA/v7h4LhmO3lGgNfAquBFKfr9sGfcxKwEWhq/9zC6bp9cMxTgIfs75OBPU7XXctjvhzoD6R7eH8U1gckCjAYWOPN/etM5YcGAlnGmN3GmFJgHjD6jGVGAzPs798GhtqfUBmoqj1mY8xn5ruPjF4NtPVxjd5Wkz9ngOeBPwOnfVlcHanJMT8AvGaMOQZgjMn1cY3eVpNjNkCM/X0T4KAP6/M6Y8yXQP5ZFhkNzDSW1UCsiLT21v41VH4oAdjv9nO2PVblMsaYcqAQaO6T6upGTY7Z3X1Yv+kEsmqPWUT6Ae2MMR/4srA6VJM/5y5AFxFZISKrRWSkz6qrGzU55ueAu0UkG/gQeNQ3pTnmXP9/Pydh3tpQPVLVjOPMW+RqskwgqfHxiMjdQApwRZ1WVPfOeswiEgK8DIz3VUE+UJM/5zCsU2BXYs1GvxKRnsaYgjqura7U5JjvAKYbY/4qIhcBs+xjdtV9eY6o03+/dKbyQ9lAO7ef2/LD6fC3y4hIGNaU+WzTTX9Xk2NGRIYBzwI3GGNKfFRbXanumBsDPYHPRWQP1rnnRQF+sb6mf7cXGmPKjDHfANuxQiZQ1eSY7wPmAxhjVgFRWD2y6qsa/f9+vjRUfmgdkCQiHUUkAutC/KIzllkEjLO/vxVYbuwrYAGq2mO2TwX9GytQAv08O1RzzMaYQmNMnDEm0RiTiHUd6QZjTKoz5XpFTf5u/w/rpgxEJA7rdNhun1bpXTU55n3AUAAR6Y4VKnk+rdK3FgH32neBDQYKjTE53tq4nv46gzGmXEQeAT7BunNkqjEmQ0T+D0g1xiwC3sCaImdhzVDGOldx7dXwmF8EGgEL7HsS9hljbnCs6Fqq4THXKzU85k+A4SKSCVQAPzfGHHWu6tqp4TE/AfxHRH6GdRpofCD/kigic7FOX8bZ14kmA+EAxpjXsa4bjQKygCJgglf3H8D/7ZRSSvkZPf2llFLKazRUlFJKeY2GilJKKa/RUFFKKeU1GipKKaW8RkNFKaWU12ioKKWU8hoNFaWUUl7z/5xsxVZgYGmRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.optimize import fsolve\n",
    "\n",
    "delta = 1e-3\n",
    "result = fsolve(pycalphad_grad, [1.00000000e+00, 1.01325000e+05, 2.50000000e+03, 7.01402806e-01,\n",
    "       2.98597194e-01, 3.00000000e-01, 0.00000000e+00, 0.00000000e+00, 0, 0])\n",
    "print(result)\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = np.linspace([1,1e5, 2500, 1e-4, 1-1e-4, 1, 0, 0, 0, 0], [1, 1e5, 2500, 1-1e-4, 1e-4, 1, 0, 0, 0,0], 100)\n",
    "y = [pycalphad_obj(z) for z in x]\n",
    "plt.plot(x[:, 3], y)\n",
    "plt.plot([0, 1], [result[7], result[6]])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00000000e+00,  0.00000000e+00,  8.52651283e-14,  2.38651410e-09,\n",
       "        2.85217538e-09, -5.50062396e-09,  1.16573418e-15, -1.11022302e-15,\n",
       "       -0.00000000e+00, -0.00000000e+00])"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pycalphad_grad(np.array(result))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00000000e+00,  0.00000000e+00,  8.52651283e-14,  2.38651410e-09,\n",
       "        2.85217538e-09, -5.50062396e-09,  1.11022302e-15, -1.11022302e-15,\n",
       "        0.00000000e+00,  0.00000000e+00])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "toy_gradient(np.array(result))"
   ]
  }
 ],
 "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}