wy36101299/optimal design3.ipynb

## optimal design3.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "import sympy as sp\n",
    "import math\n",
    "import numpy as np\n",
    "import datetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhoAAACHCAYAAABQ3vCSAAAYKGlDQ1BJQ0MgUHJvZmlsZQAAWIWV\nWQdUVEuT7nsnzzCEIeecc06Sc5ScBGHIQxKGJCBRRBFFBUWCoCIiIIiJKKIgKioigiJmFBFR9CEK\niAjsJeh7+/7ds2f7nL7zTXVV3a+7q0PNAMDLRY2OjoCZAYiMiqM7WRgLeXh6CeFGAQTwgAwUgCo1\nIDbayMHBFvyvZW4Y0UbKkPyqr/9d738sLIFBsQEAQA4I9g+MDYhE8GUA0DwB0fQ4ADADiFw0MS56\nFX9HMBsdIQgAFr+KQ9Yx3yr2X8dKazouTiYINgUAz0Cl0kMAYFz1L5QQEIL4YYxG2ihRgbQoRDUb\nwfoBodRAAHh6EB25yMhtq3gawVL+//AT8t98+v/xSaWG/MHrfVkreFNabHQENen/ORz/d4mMiP/9\nDhGkMoTSLZ1W+4yMW034NptVzIDgjih/+80IpiC4lxa4pr+Kn4XGW7pu6E8FxJogYwY4AIBBINXU\nBsHIWMIc8eGuRhtYhUpfs0X0YXtanJXLBvanb3Pa8A8nBMWaOf/GoUFWths+90RF2P/GFcE0cysE\nI5EGX04OdXFf5wn3JNDc7BHMiOCB2HBnmw39V8mhJva/dejxTqucxRD8PZhu7rSug+KKjP3dL5RC\nAHWNAxeCDeNCXSzXbVEeQbEetr+5BQaZmq1zQAUGRblucEYh0WXstGGbEx3hsKGPqgiKsHBaH2fU\n+dgE59+2g3FIgK2PA+ptGNXaYZ0/ai46zsFlnRsaDWyBCTAFQiAeqf5gGwgDtP6plink23qLOaAC\nOggBQUB+Q/Lbwn2tJQp5OoNk8BlBQSD2j53xWmsQSEDkS3+k6095ELzWmrBmEQ7eIzgSzYPWR+ui\nbZGnIVJV0Fpo7d92Qky/34o1w5piLbHmWOk/PAIQ1hFIpQPaf8r+tsS8xzzEvMU8xoxingIbpDUI\n6fMqw6g/PXMD79a8bHz3pWXR/8VcCNiBUcTOfKN3/oj15G8dtATCWh1tjNZD+CPc0RxoHiCPVkN6\nYoQ2QPqmjkj/yTD+D4u/x/Lf71vl988+bsgZZRjVN1j4/+Fv8kfr315M/jFGgcinzb81UXtQl1C3\nUV2oO6gOVAsQQl1DtaL6UFdX8Z9IeLcWCb/f5rTGLRzxQ/uto3RWaVLp13+8nbrBgL423yAuaHvc\n6oIw2RadRKeFhMYJGSE7cpCQVVSAgpyQipKyBgCr+/v69vHNaW3fhjge/C0LSwFAUxAR3vhbFjQM\nQPtLZEsj/i2T2IWEPBqAO34B8fSEdRl69YEBRMCErAxuIABEgRTSJxWgAXSBITAD1mAzcAGeYCsy\n6qEgEmGdCHaATJAD8sBBcASUguPgFKgBDeAiaAEdoAvcAvfAAHgMniOxMQ4+gWkwBxYhCMJBZIgV\n4oYEIXFIFlKBtCB9yAyyhZwgT8gPCoGioHhoB7QTyoMKoFLoJFQLXYDaoC7oDvQQegq9gSahGegn\njIIZYDaYH5aAFWEt2Ai2gV1gHzgEjoGT4Ww4Hy6GK+F6uBnugu/Bj+FR+BM8iwIoEooDJYySR2mh\nTFCbUV6oYBQdlYbaiypCVaLOodqRuR5CjaKmUAtoLJoVLYSWR+LTEu2KDkDHoNPQ+9Cl6Bp0M7oH\nPYR+g55GL2PIGD6MLEYHY4XxwIRgEjE5mCJMNaYJcxNZUeOYOSwWy4GVxGoia9MTG4ZNwe7DlmMb\nsdexD7Fj2FkcDseNk8Xp4TbjqLg4XA6uBFePu4YbxI3jfuBJeEG8Ct4c74WPwmfhi/B1+E78IH4C\nv0hgJogTdAibCYGEJMIBQhWhnfCAME5YJLIQJYl6RBdiGDGTWEw8R7xJfEH8RiKRREjaJEcSjZRB\nKiadJ/WS3pAWGCgMMgwmDN4M8Qz5DGcYrjM8ZfhGJpMlyIZkL3IcOZ9cS75BfkX+wcjKqMBoxRjI\nmM5YxtjMOMj4hYnAJM5kxLSVKZmpiOkS0wOmKWYCswSzCTOVOY25jLmN+QnzLAsrizLLZpZIln0s\ndSx3WD5QcBQJihklkJJNOUW5QRljRbGKspqwBrDuZK1ivck6zoZlk2SzYgtjy2NrYOtnm2ansKux\nu7FvZy9jv8o+yoHikOCw4ojgOMBxkWOY4ycnP6cRZxBnLuc5zkHOeS5eLkOuIK69XI1cj7l+cgtx\nm3GHcx/ibuF+yYPmkeFx5EnkqeC5yTPFy8aryxvAu5f3Iu8zPphPhs+JL4XvFF8f3yy/AL8FfzR/\nCf8N/ikBDgFDgTCBwwKdApOCrIL6gjTBw4LXBD8KsQsZCUUIFQv1CE0L8wlbCscLnxTuF14UkRRx\nFckSaRR5KUoU1RINFj0s2i06LSYoZie2Q+ys2DNxgriWeKj4UfHb4vMSkhLuErslWiQ+SHJJWkkm\nS56VfCFFljKQipGqlHokjZXWkg6XLpcekIFl1GVCZcpkHsjCshqyNNly2YdyGDltuSi5Srkn8gzy\nRvIJ8mfl3yhwKNgqZCm0KHxRFFP0UjykeFtxWUldKUKpSum5MkXZWjlLuV15RkVGJUClTOWRKlnV\nXDVdtVX1q5qsWpBahdqIOqu6nfpu9W71JQ1NDbrGOY1JTTFNP81jmk+02LQctPZp9WpjtI2107U7\ntBd0NHTidC7q/KUrrxuuW6f7YZPkpqBNVZvG9ET0qHon9Ub1hfT99E/ojxoIG1ANKg3eGooaBhpW\nG04YSRuFGdUbfTFWMqYbNxnPm+iYpJpcN0WZWpjuNe03o5i5mpWavTIXMQ8xP2s+baFukWJx3RJj\naWN5yPKJFb9VgFWt1bS1pnWqdY8Ng42zTanNW1sZW7ptux1sZ21XaPfCXtw+yr5lM9hstblw80sH\nSYcYhyuOWEcHxzLH907KTjucbjuzOvs61znPuRi7HHB57irlGu/a7cbk5u1W6zbvbupe4D7qoeiR\n6nHPk8eT5tnqhfNy86r2mt1ituXIlnFvde8c72EfSZ/tPne28myN2HrVl8mX6nvJD+Pn7lfn94u6\nmVpJnfW38j/mPx1gEnA04FOgYeDhwMkgvaCCoIlgveCC4A8heiGFIZOhBqFFoVM0E1op7WuYZdjx\nsPnwzeFnwlci3CMaI/GRfpFtUZSo8KiebQLbtm97GC0bnRM9GqMTcyRmmm5Dr46FYn1iW+PYkKtO\nX7xU/K74Nwn6CWUJPxLdEi9tZ9ketb0vSSYpN2ki2Tz5dAo6JSCle4fwjswdb1KNUk+mQWn+ad3p\nounZ6eMZFhk1mcTM8Mz7WUpZBVnfd7rvbM/mz87IHttlsetsDmMOPefJbt3dx/eg99D29Oeq5pbk\nLu8N3Hs3TymvKO/XvoB9d/cr7y/ev5IfnN9/QONAxUHswaiDw4cMDtUUsBQkF4wV2hU2HxY6vPfw\n9yO+R+4UqRUdP0o8Gn90tNi2uLVErORgya/S0NLHZcZljcf4juUemy8PLB+sMKw4d5z/eN7xnydo\nJ0ZOWpxsrpSoLDqFPZVw6n2VW9Xt01qna6t5qvOql85EnRmtcarpqdWsra3jqztwFj4bf3ay3rt+\noMG0ofWc/LmTjRyNeefB+fjzHy/4XRi+aHOx+5LWpXOXxS8fa2Jt2tsMNSc1T7eEtoy2erY+bLNu\n627XbW+6onDlTIdwR9lV9qsHOomd2Z0r15KvzV6Pvj7VFdI11u3b/fyGx41HPY49/TdtbvbeMr91\n47bR7Wu9er0dd3TutN3VuttyT+Nec596X9N99ftN/Rr9zQ80H7QOaA+0P9z0sHPQYLBryHTo1iOr\nR/ce2z9+OOw6PPLE+8noSODIh6cRT78+S3i2+DzjBebF3pfML4te8b2qfC39unFUY/TqG9M3fW+d\n3z4fCxj79C723a/x7Pfk90UTghO1H1Q+dEyaTw583PJx/FP0p8WpnM8sn499kfpy+S/Dv/qmPabH\nv9K/rszs+8b97cx3te/dsw6zr+Yi5xbn9/7g/lGzoLVw+6f7z4nFxF+4X8VL0kvtyzbLL1YiV1ai\nqXTq2lUAhVQ4OBiAmTMAkD0BYEXyOCLjev61UVDQatqxqktGbjGbkNtWIeiHKJAHVAPDcCQ8hgpC\nzaDzMEqYUWw5LgxvSpAgMpJgBhSZhVGWyYqZznKS8pJNgN2f4yIXmtuP5zqfIH+uwFchH+F7ojpi\npyXYJDOkJmTsZRvlGRUCFC8pLaroqsaqHVfv0XijuaDNoMOjK7NJS89U397AyzDUKME4x6TItMas\n3fyuxTPLD1bzNmhbZjs+e8nNyg46jsZOVs72Lk6urm7u7h4enp5eXl5bvLy9fLy2evi6+TlR7fzN\nA/QD1YNkggVDWENxoYu0L2Fvwh9F3EZW5dlt5dH7Y5Lo1FijOO64L/FdCUcTt223ThJNWkp+ktK4\nY0+qX5pmOiOytq5kFmSF7tTLZs3+sKszp3B36J5NuRy5S3noffr7Gw5oHbx4aKlQ8LDsEYUipaPK\nxaolaqXqZerHNMp1KsyPB50oPjlyir3K6LRPddSZ5Jqc2kN1ZWdP1zc2tJ270Th4/vNF4UvRlwea\npVsiWovbmtsfXJnoWO7kuKZ83a2roPtDj+XNslv3b7/pnb6LvSfeZ3E/sD/2QcSA60PNQYEh4tDC\no7HH94evPWkf6Xh67VnX884XjS8PvYp4bTzKPTrzZuBt21jNu7Lxg+93TSR9iJz0+2j3SXWKMvXp\n860vVX/lTId9tZ9R+ybyXXrWZ67zh9LC4Z+vf3EveSxXraysxgkgAV7kluiE5Dr14D0kCW2DrsO8\ncBY8g4pG/UDvwQhjbmLjcAq4b/huQjkxlRTI4EF2ZvRg8meOZ8mj1LAOsP3gkOT04SrkfsBL5rPl\n3yfQL0QWdhQ5JDogTpIwk0yQqpZ+KPNdjlleSkFNUVtJW1lVRVpVQI1ZHVL/rjGOnFa92m06tbql\nm/L0UvTDDLYY2hsZG2uaKJiKmfGYM1tgLRYtp63GrUds+mw77c7bV24udMh2jHWiOju46LvKunG5\nY9y/erzw7PW6tOW4d65P7FYfXzM/OSor9Yf/y4CuwKqgPcHhIXahijQW2rewx+HNEcWRSVEe2zSi\nKdGTMdfohbGBcerxmPjhhNOJ8dvNktiSxpIvpmTssE/lS/2Y1p6+PyMs0ynLFIkMnV0aOUq7ZfeI\n5wru5c6j7CPtR+9fyp878PXgzKGFQtxhriNSRZpHTYsdSraUhpTRj6WW764oOH7sxJmTrZWDpxZO\nS1d7n8mraap9VrdcL9xgdi60cf/5lgtfLqlf3tX0sIXcqtdGay+5cq9jpVP9WsT1qq4XN1h6DG/S\nbuXdruvtvTN5j9ynet+rP+tB/cCTQeyQ2iPfx9nDVU96Rt4/Iz5XfOHyMulVxevbo/Nvlcfo7y6N\nz0zIfQiZPPXx9RTvZ48vx/6a/prwTX6WMk9cgH9++nVlmbYx/0TACeSAJZLxHAV3ISxkAR2CxmA9\n+CSKjNqFxqELMBKY69hAHAV3B7+HYE8UJC6QHjG0kk8zljAVMB9gKaCUsp5ma2bv5XjFucBN4ZHn\nNeej8u8QOCp4Tqhb+JHIuOhnsRnxaeTWNCLVLX1aZqesl5yiPCQ/qFClmKhkrSykPK/Sr1qtlqbu\nqiGnCWuOaDVoZ+m46croLm0a0Duln2hgbShoOGvUZ3zaJM3UzUzBHGP+wuKy5V4rP2sNG5LNqG2T\nXa69D7JTYByeOtY7ZTg7u4i7zLn2upW6h3voepI8n3ud3ZLsbenD4fNu6wXfdD8bKid1zP9cQHKg\neRBL0PPg6pCYUF0amtYfdjTcL0I64mtka1TGNvNoQnRfzD66dSw+9mZcVrx+/GJCS2LsdoXtk0lV\nyb4pPCmPduSnWqbBaZ3paRmWmfyZi1mjO3uzL+wqy8neHbnHLVd/r0QeOW9234v9N/LrDhw+mHko\nsYBeGH0YuRYUxRyNKY4uiSqllfkdcy63rrA97nMi6WR55c1TX06zV2uesa1xqnWs23I2pf5yw2Kj\nxfnCC68vyV5OaOpqIbU6txW3P+8QvhrRefU6a1dY940e3ptxt/p7Je6k3n3UJ3M/q39swO3h8FDA\no9nhPSM8TxueG74YfpUxavfW+d3B9/OTh6dufnWZf7o6/+u/w60WLJKdnjYHwO0wAM7aCM4HQLwG\nOT82AeBABsBFG8DcJQC6Gg0gb6k/54cAMELOjp2gCtxEdg8ssn9YQeHQfqgRyfW+w5ywLuwL74Rr\n4H74G4oHZYQKRR1EMvC3aBJaA01F70e3oScw7BgzTDySdY1gGbBG2ETsOewHnAjOF1eBe4UXwYfi\nz+OXCLaEE4Q5ogOxgUQmRZEGGbQYTpJJ5ATyGKMjYxeTClMVMzfzQRYCyy4KTMlixbDmsjGzlbCL\ns1/iMOUY4dzGheeq4jbhfsezm1ee9zFfOr8c/wuBfEFTwSWhduFkEX1RjOgDsWPi4RJ6khTJj1I9\n0pUyWbJBcrby2goKiopK+squKhGqO5Etv0ljSHNOm1/HQjdhU63eawNuQ3ejEuPXplJm8ea3LHms\nQqyP2By1TbAztFux79q8zyHMkeaU7Xze5Z0bj7uzR75n3xayt6NP0dYRPyaqqr9FgGtgYFB68NmQ\nDzTlsMzwoUgpJPKexWjSi2J/xLsn1Cd+TuJMVkox3uGZmp7WlkHIDM26n62xq3I305603Ik8o33Z\n+5vyRw8yHrIvOH9Y7cjNo/bF90sty26VO1b8ONFb2Vl1sfpoTXIdrX7LOaPz7BfeXGpoSm/Z2uZ1\nZcfVlmsL3do9kbf29pbcrepr7O8ceDg48Rj/RO/p/uffX3mNNo2RxqkT7R/xU5JfwF/lXwVmir/z\nzTbPRy6o//z1q3nZd23/EAM2IAYUgQ7wFsJDCpALlAxVIpn+V5gbNobD4SPwdfgTkrObIKdJOaoP\ntYiWRXuj89Fd6FmMDIaKKcY8xJKwFthd2B4cFmeNO4AbwYvhY/E3CXyERMIwUZt4gkQkJZImGDwY\n7pNNyB2MWozNTBpMbcwGzLeQHPUpJYgyw5rFxsZWw27A/pQjnpONs5nLixvmrufx5CXwdvDFInP9\nQeCMIE1IQWhOuEvkgKivmKo4UfydRLdklVSudKyMn6yjnJn8JgVNRXUlDWUdFWNVO7Ut6lEaOZrV\nWg+0l3XVNm3TO6c/Y6htlG08ZCpplm7+3FLXqtx6ydbBrtD+7uZfjgpOgc4VLs+QOfb2OOn5cYu6\n906fIV9xv1hqh/9yoF5QanBXKIHmFnY6fD7SNurUtl8xnvTWOO74HQnPtislpSRfTfmZqpOWmd6f\nKZKVtHNol1JO/u4vufZ76/IW9xvm7zjQdHC2wLSw8gihiH50pES/9NQxfPm2iuETeierT7FV5VZj\nz+TXCtRdrrdrGGvcfoF08fhltaa7LX6ts+17OviuNl1z74K7m3pot/hu99/JuKfW97H/1MCWQZah\na48DnoCRsmfaz1++3P1aZfT1233vdMenJiom7T/OTu35vPCX1fSurxdm+r99+L4yxzWv+sNlYcfP\nusWPS1rLR9bmXxq4gHRQC4bAMiSNzH4GVA+NwFhYDfaHD8FdyC1CFOWGykVdRX1DS6P90CXoIQwT\nxhaTh7mHJWOdsCXYtzh5XCruAV4Cn4l/Q7AgXCSKEctIHKQjDBwMpWRBcjWjEmMHky3Ta+S+wcRS\nT7GjfGUtYTNlm2E/xeHGSeLs4krm1uCe42njTeez4udE5vqq4BEhOnIDURPlEUMjZ8+YxFPJQakH\nSGb+WPaV3Cf5X4oUJTlla2RFF6p1qn/WFNJy1y7QGdzEruejX2ewaORgXGdKMIs0f2JpbXXLxtZ2\nxJ7mABzLnDe5vHXL9zDwnN1ywYfuq+E3418WKBvUECITWh0mEV4bqRjVFm0RMxIbGY9NqNxulPQ6\nZXsqNi0/gy2zZKdIdkOO7u77uQF50L6z+d4HsYfKCgUOHynCHU0snij1Khsq96j4fqK2MqgKd3pv\n9VyNR23TWbb6uIbhRu3zFRcxl6IuP222aGlrU2qv6xC7WnaN8Xpq18cb7j09t1Run7xDuZt9b/5+\nRP+7AZ+HT4fcHz0Zdnly56nKs4Lnn17qv8p//fKN/Nu0sYFx0ffbJ+5Pin5M+HR9avmL8l/W055f\nPWfsv236LjqLm3071z6f8UPvx/RC5k/KzxOLhMWYxae/jH+V/PqwpLm0a+nRsugybblheXpFbWX7\nypXV+Y8NVlVZOz4gBmMAMK9WVr5JAIArAGDp0MrKYuXKytIpJMl4AcD1iPX/dtbOGmYAjt1YRbeS\nxzL+/R/LfwGxasSrzZ1JGAAAAZ1pVFh0WE1MOmNvbS5hZG9iZS54bXAAAAAAADx4OnhtcG1ldGEg\neG1sbnM6eD0iYWRvYmU6bnM6bWV0YS8iIHg6eG1wdGs9IlhNUCBDb3JlIDUuNC4wIj4KICAgPHJk\nZjpSREYgeG1sbnM6cmRmPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5LzAyLzIyLXJkZi1zeW50YXgt\nbnMjIj4KICAgICAgPHJkZjpEZXNjcmlwdGlvbiByZGY6YWJvdXQ9IiIKICAgICAgICAgICAgeG1s\nbnM6ZXhpZj0iaHR0cDovL25zLmFkb2JlLmNvbS9leGlmLzEuMC8iPgogICAgICAgICA8ZXhpZjpQ\naXhlbFhEaW1lbnNpb24+NTM4PC9leGlmOlBpeGVsWERpbWVuc2lvbj4KICAgICAgICAgPGV4aWY6\nUGl4ZWxZRGltZW5zaW9uPjEzNTwvZXhpZjpQaXhlbFlEaW1lbnNpb24+CiAgICAgIDwvcmRmOkRl\nc2NyaXB0aW9uPgogICA8L3JkZjpSREY+CjwveDp4bXBtZXRhPgrD7+IRAAA9DUlEQVR4Ae2dCdxV\n0/f/tynzmMxDMlSoRJSvKKJUivRFSQkVSjJlppApEZlCyJjMSYaIkCEhmWUWQmaZh/Pb7/X97/Pf\n9zznDs9z732ec+9d6/W6955hn332/uxzz/7stdZee4nAilFRBBQBRUARUAQUAUWgCAgsWYQ8NUtF\nQBFQBBQBRUARUAQEASUa+iAoAoqAIqAIKAKKQNEQUKJRNGg1Y0VAEVAEFAFFQBFQoqHPgCKgCCgC\nioAioAgUDQElGkWDVjNWBBQBRUARUAQUASUa+gwoAoqAIqAIKAKKQNEQUKJRNGg1Y0VAEVAEFAFF\nQBFQoqHPgCKgCCgCioAioAgUDQElGkWDVjNWBBQBRUARUAQUASUa+gwoAoqAIqAIKAKKQNEQUKJR\nBGj/+OMP8+eff+ac808//WT++uuvnNOnS8h9yYuPiiKgCCgCioAikAQElk5CIaJlePvtt83vv/9u\nfvnlF7PVVluZ1VdfPZok5/2nn37afPLJJ7Hp11xzTbPllluaDTfc0Cy5ZP6ci/Jyv2effda0b9/e\n7L777rH39Q9++OGH5sQTTzQdO3Y0gwYN8k9Va/u7774zjz32mPnggw/M3LlzzbXXXpsXbtW6eUxi\n2u+zzz4zG2ywgVluueViUughRUARUAQUgUpAIP/etYAo/fjjj2bs2LFm1KhR5r333jOvvPKKGT16\ntPn5559rfJd//vlHtAsjRowwN954o2yjbSDPqVOnmmHDhpl58+bVOP/ohX///beU+eOPP46eit2H\nBN1///3mlltuiT2fy0HWxZswYYJgVr9+fcnv008/zeXSgqR59913zeeff56S18svv2wOOeQQc+ut\nt6Yc1x1FQBFQBBSBCkOA1VuTIjNmzAh23XXXwHb8gR0RB1ZDEPTu3Tt44IEH8i7isssuGwwZMiTM\n599//w2+//77oEuXLsE666wTHi/ExvLLLx9cd911OWVlzR2B1WpIPRcsWJDTNdFEltwEVjMTPPro\no4E1m0h+5FtbcsEFFwT33ntvyu0s0QnOO++84J133kk5rjuKgCKgCCgClYVAYjQaP/zwgxk3bpxZ\nddVVzcYbb2wsMTArrLCCWXrppY0lGnnTP/wXfFliiSXMaqutJmr9L7/80rz22mv+aWMfAzHdFMJ3\nIiXjyE69evXMGmusYTB9YGZwgmYEU4wlRO5Qxt/ffvtN6rLyyiubTTbZxJBvbQjaoYceekjw8u+H\nOeqUU04xjRs39g/LNmaVxYsXG7RNKoqAIqAIKALljUAifTSikGMWsBqC6OG89zEvfPXVV2a33XYz\nzZs3D/P7+uuvjdUOSGfIwW7duoUkwGpBpGOlA+X6zTff3Gy99dYG4pJJyG/FFVc0EAhIT6dOncLk\np512muGDQDjmzJljlllmGYOvyiqrrCJk4+CDDw7T+xuknz17tpSVe+AI2rVrV/Pggw9KXlZjYxo0\naGC++OIL8eHAtLLXXntJFpg33njjDdOyZUvxUWEf2XPPPc3aa68t23xBCPA7+eijj4Q44FdCvhC0\ne+65x7z44ovGai7MTTfdZLp3724gO6THfNSkSRPTunVryYuy4T/C59dffzWbbrqpadOmjRBLyMdz\nzz1nrFbHbLvttoKTI3+uDq5A1HO99dYzzZo1c4f0VxFQBBQBRSChCCRGo0HHuu666xYVJjrKKVOm\nyOfOO+80J598sjiCXn/99eF98TU44ogjxOdgxx13lA6PfSeXXXaZueuuu2QEjx8JeUS1JS6t+734\n4ovNzTffLNegoRgwYIB55pln3GkhMjvvvLPs09nis4GWY6eddjKLFi0y/fv3D9PGbTitB79u+/XX\nX5frKCOChoZy33DDDWEWHDvrrLOMNX0Y8ODamTNnmiOPPDJMwwZ+MhAKznMN5YfcsM2xpZZaStrO\n3ZtrIDZnn3221Jt9fGIuv/xyqTcOuJA7iMjQoUM5LcL14Mkxrmf/mmuuMccdd5xLIr/gUQzimXIT\n3VEEFAFFQBEoDAK2s0iE2E4lsJ1KYEepgdUoSJnsaDto1apVYGuadxnJw3aQ4peBb8bdd98d7LHH\nHsEOO+wQ2E4wzN86MAaDBw8WHxEO4utgtQrB7bffLmmsdiOwJgHZxreC621HLPvuy/fRuO+++4KV\nVlop+Oabb9zp4Pjjjw+sRiPc9zesViCwGoXwEPe32pJwP24DHw1rLgmeeuqp8DS+Edtvv31gNQvh\nMcoNBr5QfjurJ7CaETmMjwf+LJaQyT71pv6Uywnn77jjDtnFl4bro2JNKsHAgQMFS85Z7UxgiVtg\nHXzDpFarEViNUGBnyITHrKYosBqZgOuRMWPGBGuttVZgCVOYxjqfBs8//3y4rxuKgCKgCCgCyUUg\nMRoNTA89evQQM8T06dMNMxmYsYBKvlCC3wd+GXx69uwpo23MCoyuGXEzA4T74t+AzwN+I6j4Se+0\nFrbzNG3btpW0lkDIOdKlE9v5Sx6M+knHh/zSCXlyDf4axxxzjGgOKEN1BQ1R3LRSTCVR2WWXXcKp\nsJg97OMq/iGkw0yB9OvXT375wsxxwAEHhPtxG9TTkqvwFBqKpk2bhiYoTlhCJtoQNDi+0B6UH6Et\nuB8aGidbbLGFmFzcvv4qAoqAIqAIJBeBRPlo0OljZsCGf8kll5h99tlHTBhMcy2G2Nkm0oGiisck\nYjUc4teAnwCdvROcVFu0aCGdL/E2cH7EBwLfBnwkMglECROAnx++Be3atYu9DF+RkSNHCsGaP3++\n+EJYDYu58MILY9NX9yD+JFGxWo3ooXA/jpiEJ+2GTwDwwYAUgWtUvv32W/HpiB5nv1jtG3cvPaYI\nKAKKgCJQuwgkimjgXLlw4ULp8HEARMPAMTraYgsaFZwTcTLEwXPvvfdOuSX+HThp4j/gSBAJ8PlA\ncBJl1E6AMV+YfUHe0fzSBRHDH2L8+PHinImj6qRJk8yxxx4rfhxoUmpbGjZsKPE5/PtCnKgrWhdw\ncQLp4OP7tLhz/BKvBM0RhNIX8M4kzDwqtv9OpvvrOUVAEVAEFIGaI5AY0wlVwEGSAFqo3VGXT5w4\n0TADhFkfTnAozKa2d2mz/TJrgpkS3ItInkQKtX4NElXTv/aqq64SjQSzICAIHTp0kNN0thAMhJkX\nkKSoELSKjnLy5Mkpp4466qiUfbfDjA60KwgzPyBZEBW0J5kEbUI0aFam9P45yu7E11BwDGKFGemi\niy5ySWSGy+mnny7766+/vjjMYnKCRECOEGa3+DNXmF0CCSNiqRNICnkzy8UXPw3aJma94BTrBBNb\n37595dlwx/RXEVAEFAFFIJkILGXV9COTUjQIBqp3OiRiZ7zwwgsSJdSfxkgETc6ddNJJORUbjQNp\nISyYPSAKmDH4MAPj1VdflRkVzpRB2HB8KCAXCFoM0jKDYqONNjKYM5jJwSicmRp0lHSYfNBeQI4g\nC5gDmJrKzBG0GdYBVEwHkJpTTz3VbLPNNkJuopXgWmbBQGrwU+D++DL06tUrmlT233rrLbkHfhVP\nPvmkefjhh6UTJmw702mpByYhtCLWsdMwhRWNBCSAuoMPdSItfiiQPe7JNFy0O2hR+Fx55ZWiwXjk\nkUckzzPPPFPIBHWGnEByaDuuY7YOfheUBdzBCpMNM2sgipAz6slsF66/4oorDKYV62Aqs14oB7FU\naK/hw4cLgUMbstlmmxk0LNaRV/LAtySqHYkFSQ8qAoqAIqAI1BkCMp2hzu4ec2P8M2baKZask0En\nE7X3ozUgDaaVXAStA+nTBYdCgxGnliccuvOr6Ny5c+icSIdKfnTWmFogRdOmTZOOlM6RchNoixE+\nHSwOkAidq5vSijMj8SXiBP8QOl86XmJRsJ1Nm0E+3A+NDPeGEDmZNWuW+JHYmRsSfp1pqUwvJdYG\nJOX999+XpBCPRo0aSd0gQwh1o47UlTrzIaAaJMk3dxDUDFKBtLeaIUgLacEe4X6QBCdvvvmmrMmC\nPwrEwQlkwpmU0JTwgQxxL8rhBDKDhiOu3Vwa/VUEFAFFQBFIBgKJIxrJgEVLoQgoAoqAIqAIKAKF\nQCBRPhqFqJDmoQgoAoqAIqAIKALJQUCJRnLaQkuiCCgCioAioAiUHQJKNMquSbVCioAioAgoAopA\nchBQopGcttCSKAKKgCKgCCgCZYeAEo2ya1KtkCKgCCgCioAikBwElGgkpy20JIqAIqAIKAKKQNkh\noESj7JpUK6QIKAKKgCKgCCQHASUayWkLLYkioAgoAoqAIlB2CCjRKLsm1QopAoqAIqAIKALJQUCJ\nRnLaQkuiCCgCioAioAiUHQJKNMquSbVCioAioAgoAopAchBQopGcttCSKAKKgCKgCCgCZYeAEo2y\na1KtkCKgCCgCioAikBwElGgkpy20JIqAIqAIKAKKQNkhoESj7JpUK6QIKAKKgCKgCCQHASUayWkL\nLYkioAgoAoqAIlB2CCjRKLsm1QopAoqAIqAIKALJQUCJRnLaQkuiCCgCioAioAiUHQJKNMquSbVC\nioAioAgoAopAchBQopGcttCSKAKKgCKgCCgCZYeAEo2ya1KtkCKgCDgEnnjiCXPwwQebHXbYwTz4\n4IPusP4qAjkj8P7778sztM4665iTTz7ZfPXVVzlfmy7h1VdfbTp27Gh23XVX8+abb6ZLVjbHlWiU\nTVNqRRQBRcBHYOTIkeamm26Sz1FHHWX69OljFi1a5CfRbUUgIwJffPGFOeGEE8w+++xjPv30U7Pm\nmmsanquffvop43WZTo4YMcKsuuqqZurUqaZz587mwAMPNG+99VamS0r+nBKNkm/C9BX48ccfzZ9/\n/pk+QYWc+fvvv833339fIbXVajoEJk2aZFq0aCG7/fr1M2uttZa59dZb3emi/y5cuLDo9yjlGwRB\nYL755ptEV+G7776Td+i6665r6tWrZzp16mRefPFF8/XXX9e43D/88IN59dVXzbLLLmuOOOIIs8km\nm5jrr7++xvmVwoVKNEqhldKUkYcdVjx58uQqKRYvXmwuv/xyM3/+/CrnMh2YPn26OfvsszMlKblz\nzzzzjLnmmmtKrtyFKvAvv/xixo0bZ954441CZVnr+WD2OP3006vc9/fff5e25eUflauuusrsvffe\ncvjzzz+vVdJ95513mtdffz1aJN33EHjttdcMZDDJstFGG5mTTjrJrL/++lJMSMIKK6wgpIMDEIYH\nHngg9tnq37+/+eyzz6pUb8iQIeaQQw6R4xCtb7/91jRu3LhKunI6oESjRFuTB/6MM84wH330kTy0\nd9xxR1gTRgo33HCD+fjjj8M/SHgyy0b79u3N008/be6+++4sKUvj9Jw5c8yAAQNMly5dqhT43Xff\nNbNmzapyPO7AK6+8Yt55553wFPtcXwry4Ycfmrlz58oLMlpe6l8KHSL2bEaTmD98YVTIiPDMM8+s\n8lLv0KGD2XTTTSU5z3PDhg3N4MGD/cuLsk3nc//995uWLVum5A/xv+CCC+TjP0spiSI7Lj2/f/zx\nR+Rs6u7zzz8vZqLUo8nd413zn//8p0oBeW/l+r+scnGBD6yyyiqmXbt2ZsMNNxQyMWPGDLP99tuL\nCYVbbb311uL7M2XKlCp3pn5Dhw6tov3YYostTNOmTSX9U089JUSje/fuVa4vpwNKNEq0NXFQQmux\n5JJLGlj3yiuvHNbkn3/+MWgmBg4caFZfffXweC4bqAdh6KjyGHGUutx3332GP3Hz5s3Dqjz88MPy\nAqAjokPIJtj16ch4OSyxxBJmmWWWMfvvv3+Vji1bPnV1Hm1G165dTaNGjaQImNPQEBx22GFml112\nERJSV2XL9b48lzvvvLMkv/baa8PLaA9ICCaSW265JTzuNiDd1JVO+K677jKYE3MVNEF0etURbPeQ\n/FatWpnVVltNLuW/CtHFhHf44YeLrb9Xr16SLl3emF0Y5W6++eYykPj333+lrR555JHwEsgH9edD\nOkbQaG6SLhAmHHTHjBljtttuu7C4Z511ljnxxBMFq5tvvjk8noQNnoWJEyeaTz75xBx33HFmueWW\nk2ItvfTS5tRTTxVtW9Scgu8FzwEkl3dyVGbPnm1mzpxpbrzxxpT3dzRdOewr0SjRVuSFxsuITnTa\ntGnSkbiqvPDCC9IZwpxrIvyJGDlikuEPVqrCC56X/jHHHJNSBUbA3bp1q5a6Eo9zXo68ZE455RQx\nWUFUki48C19++aXZbbfdwqJCUOm4jz322PBYqWwMGzZMnsuo5z9tgyc/znu+vPfee6J5gijSWdN+\nuQqaHqfizvUazFP8N+lkIKQIZh2OkRfE//zzzzfrrbeeGTVqVNpsObfxxhubnj17mrXXXls6M9qN\nOvqy++67S514DzCqptNLukD6MBlAdH3ZcccdzR577BH61fjn6nob356lllpKtFE8Fz5hbWg1ZXvt\ntZe5+OKLU4rJfwxnT7QzmEd84Z19++23m4suusjUr19f3if++XLbXrrcKlQp9cFk0qxZM9OgQQOz\n/PLLp1QbTQZe9tXVZrhM0JKg9n388cdlKpcbCbvzpfJ73XXXySiYF7YvTZo0MXyqM2pCdfrf//7X\nbLXVVn5Wid+mczvyyCPNGmusEZaVbbQApSi0G8885glMKU623HJLeWYZ8R966KFyGP+kvn37SieP\nihqH4G222cZdktMvI9jqCNNp+b9ATJ3QKT377LOhHZ8ZByuttJKYPV2a6C/k8LHHHks5jFYKrQxm\nu2233VbO7bfffmbQoEEp6ZK+Q8eLWSEqPJPgDQnxNbTRdLW9j9Zz9OjR8t+HzH3wwQeiJfPfr7xv\nIRuYSjbYYIOwiLwvILg8r+4/RxsOHz5c6sh7nAFjbZj0wkLVwYZqNOoA9HxuiTmDEdH48ePFCY7R\nzoUXXpiSJSN53/b5119/yT6jHq6FRSNPPvmkdLjYISEWqPKcMBojH3xBSlV4We+7774FKT4vA16A\nl112mZkwYUJB8syWCQ6FkAJiQKDCZ1SOOv7444/Pdml4ns4WopEEwU8E9TijVlT8vJSp3znnnJNz\n8VBVI7ygo4KGifydYE/nWI8ePcSWjunlgAMOcKdz+sUeXx3BATVKZiAIjGCdKSWX/BhBP/rooylJ\n8QkiD98MSAKwZHZEbQhTPKlP27ZtDU7WV155peH9wX4uQrvRLv77yb8OQu9mCvnHi7XNfxpfN4gh\nJsa3335bNC3813gfYjZFe4apC38MPkx35b3rC1pS2gES6As+RGiIMWU7YSCHxhitKvnxf2D6bDmL\najRKrHV5yaDJwO6MzZA/gC/Y+7DZ+o5okAZGeqhjmVGC3bp3797ykMPK6bh4WaACdILZpU2bNgUx\nnXA//oCUORehPIz48pFff/3V8Nlpp53yyUaupXNDbQrJoGy8nMDd99vI+yZeBthz6WggF9h92eZl\nRPtxf3DEydV1gqSBTDrPeJcVGhuwdyp8d7y2f5kZgl/FvHnz5JcXOiM9RraMXFGlUwemn+IkiYmA\nlzPmBV8TQ7nx+MecgDc/I0Ve2k6YguhP586XYK244orS7i7/XH5Rl0eJBj5UfJygXSFuwiWXXOIO\nVflFS+lGwJzEQRsNI8cc2eL4pZdeKm2P1oSZVajv6cDo4AotP//8s4zEeYcwqse/h3cQ7xw665df\nfll8LvDBYIBCG/k+GJSHNsOfBHNBnKB1osOuDY0GzyW+aJifILtoWfhvQSTQUPB/at26tfjVxJU1\negyzHu+cqEAifKLB88unkkSJRom2NiNsXqy5CkSCqHZ0PLyMiCuAihd1H6aWdIKzknPCS5cm23Fe\nLnR26UYx0euxZ+ZLNBjJY58vhDi7Op0a5aIzY3R+7rnninYjarqi00cblC12B8F/eOlCDH3BL8aR\nQToV8iMdpIkXNCYDRzJQwzJCxJxQHc2Af79s23QacQGFsFVDen2Ja2NIBLihfkaoz5577inPL50V\nH0gGRAS/IM4zcqajx8cmTjj/22+/SRp3HjMCUwd5ftJ1ZC5t9BeNAG2Gw6YT/mOMwJ977jl3SH4Z\nzbrZLP4J/EMgSD6p8M+zja8GDspoWoihkIvwLOM4iUPp2LFjw0vQdNE5MiBA+OW/zP/amVbCxHaD\njh/cMgkEZbPNNpOAUtF0EAl8TyAPToMCYeR+kBxHKiCLkFwIJW3pi9OQ5vv/Js9ou/j38bfRNsS9\nK/kfU19mlSC0HwMwBmm8H33TnJ9fum2eRwZU0cEf6fHHQNtVqaJEowRbvqbTERkdEosA7QZmE9eB\nZoIgbh54pvRx59AGRNW9cemSeozRIk6V7oXOqBWSwIiUl0uUaDCi44VDx0UHjf8AL186wAULFoQe\n6IzaGNVHr2fUzqgUcojgTIiGiQ7Ot/9yjk6NfCEaxRLs5rfddpthRMtLm1EfjnuMCDFP+RJHNFBL\nE88C3wQ6cuzVEBT2IRxOGBXTgdOZfmy1OYw077nnntjRH50COENGfEEbhKq7ukQDMsR/gjo5oc1p\nt6ivBB1RHNGgzNkELRjtzag52u5x10J2zjvvPOn8sev7Habv4Mu17ENucUiNIxrE3KGeCHWlg/U1\nQhyH2EKC3LPOMSe0mcPVOd0yqwLxByM8o+meR55VplunE7RZufpooPmII8DRvHHI9HFz5yE7TH1H\nIHOYSTALUXY/Zgv34VmAJDHAwIk6rv15th0Bc/fQ3/8hoESjBJ8EHJJwcMzVLupXkZe7U0czjZUZ\nFOleClznm2D8fKq7zQv73nvvzemypHnOjxw5Urz96Wx9lScdL5qHqND54QvAOTo+iBa/7KPd8YUR\ncFR42bsXOufAjc4mzgSCatcRkmg+hdqng4FUUX46JgIYUadoJ5XufhAzVx9UyxANVNO+CYBrIW9o\nNhDMeBAuVPRxwrPvO+O5NNwLh9HqCkQuOiOFETmdLipxX9LVO65z9q/D9AYRwEeKtvTJq5/O34aI\n8h9EU+PjhSYLjQiaDn9mDOTFd0T183KdKseYRkobRrVpnEtXP9eGpMGcl+5/SrqoDwPXILQZfhjp\nBC0D5D0Xwf8qF98G6hkn1NNpVnDiRbvhiJOfHr82/mf8pxk8cF+myGPa84Xr0830Q/tUyfL/DZyV\njEKJ1Z1RMS/kuD8QL+t0wgiQOd2M0Bn18Kfi5RgXeIqRHaMLp6JPl2cux2H6jAooc/TDKCJ6rKYa\nG78svMwYWVGHTBJ3npELOLnRCepoHGndCBKyAJYQtriXMi9vOgU6E4iE+6UTgIT4n7gXPaYTtBi0\nEyNHOlzU2Qht5Qdny1Q3riuEUEdeotSFjpx82ade0U/c/cATJ2aeKUbo5OHIER0WhBfBTALxRZil\nAcZR9TUdESr8uE4Sx0Qkrk3kRIYvronWhXtEj7Ef979z9+VZd+YBdzvKjNp8pjVDQhjRKvB/YNsJ\n+/jdOKGzBQ+0R8x2QktDGjp3tIwQVwYNPtl66KGHpGxx5JV8/brwLNIO/jG3na5+PJNoG3gG0WSh\n3UJo23QmLkngffGf4QMm+Uq69nH1cL/pngccPyEOCP4vEDTS8pxCLhDek/wHMQtBStDA0b5xvhgc\nR3sZFchlNhIavabc9lWjUYItyrRNFouKE1TujDxZEdBNxeSlhF8GHedE67zFbAbO4Y2PKhC1NiMm\nHA7daBBbLoSGl1G+QqfLyLA2hVEVH5ywon9+1KSYNj62qm5e2jjngRmjDjpRXqLY++nQGcXwQucl\nRAfBi50ROS8mVNluRFTIutF2jJqImcBsBe7JS4xy4oeRru2jZcCMAZGh3JA5J3Tg+K84HxZ8POjY\nUBtHTTPumnx+mZ3AbBCIwGmnnSZY0gkwq4IQ1NF1HmgT0h500EFVZiDQkePnEKcKp4yFIMZ+XSE1\n1RGeH7QQ+D4h+ChAANAIga2vffDbBNKBnwN151nDxwFNCv4P/jV09KTlf4pjqOvAIOw8LwSSc//7\n6pQ7W1qeIe7FPagjJhDuT0eMhpX75iKQILSbECfq5gvPCcQKQsV/F7IJKXFRNP20hdhmAIHDJ/9j\nzIMMJCAQkF9nkmPg1t76wTiTmiN4vnbJlYX3Kj4eUeF/W4w2id4nyftKNJLcOjFlo4Oj88jUceOM\nREfqP9yo+XmJYfd2gn2f0RICs3dOUezzZ2fEUJ0peVyXJKGjxb4eDdjFqIWZDkRebGhV9HS0aAyc\nfwHaEDoHN3uAFz8vUl7wvNBRmbK8M6poOsxCC/dnXj1EoH///qJNYdSIhoOZJ9HZJenuTzkxe0BM\neJk6gVS89NJLMqMIIgPxYDRdrOBsdBSQWjpP7OXOzwMCxQwG8HXCi55ZDJAMXvDOx8Wd54VPJxU3\nBRIfBMhyISWThjDuPpSbaZGOaPCs8ZyBP+1HOzrx24RyYw51fhuQXpxGfcFcgjMmwrNJR4mDNySA\n/zfYYpJy2iL/2ny3IQhEIoXo0J4jR44UgkH98NfJ5ADr3xtyhXmBAQDX+QLB572D0zPknncdTtDF\nIho4c0NsIEsMxDBRYtZi8ODIHc/m0UcfLcWE4DJQIzic/8xyErIFAXMOsa5eaDP4+DOI3LmK+rUj\nOZUSQMCOvgM7xS2wL+HAEo2MJbYjxcCq+FLS2JFgYEdX8rGsXM7x647xa00Ccpy09qUS2FkC4bGU\nzEpkh3rYl3dgVb4pJaaefr3ZJm0m4RrbyQW2Iwxsx5w1faa8cjnn2ou0rrzpymgDVAXWeS02W0u0\nAtthp5wjP/Jy+bl9foshLn9w45ljPx3mVqsSWA99SWOnxAb2xZ9SJGs6CCw5rvJcWmIYWF8jyTfl\ngjx2rFNqYKePVisHS4wCa2YLbOefch11jv7f/ASc45NJwCyaxj2TPJeuPTPlkc85V35XVvcfintu\nrCklsKP+2NtZkhFYLWqVcy5/l2+6Z6TKhXkc8O9hzTkZ/9vWMTmwA7XADtiq3NH638hzGT1htZ+B\n1YoGlkRFT1XUvvpolAitxOmJCHWs54DpJJOg7WDkwIjI/nklKdoJ1H18GIUg/Lpj/Dp/AdSi+Elg\nj3TH5IIS+6LOjMLc7BBXfOrk15tt0mYSrsEvwdm2s6XPlFcu51x7kdaVN+6ejMgwAWHOwTxhX/op\n2aPBYrTJaNh2RHKO/MjL5ef2+S2GuPwxw/HMsR+HOaYCprjiqEiZo9oONBksTEZ8DL+saGiYPkgQ\nO/ItlKAJizqCZssbLRIjcJwFHd5cQ3mj/zc/L87xySTULZrGPZM8l649M+WRzzlXfldWysPHbwvy\nf9LOaEEbgXkLDQ8mEV/Q3uDgSjuiwXLi8nf58lvsOvn3QDvJMxp3T/yI0ILy7sWkjAbJCU69OPqO\nGDHCHZJfS2LMTOuXgy+c70uTkqhCdpYaaaVC6lrS1URliekDVThBmzIJf3zUfzg0Mq0UO2eugvqS\nKXiQlUKroXMtQyHTEXkPdTXOZ+k8wgt5v9rOi2eCqYXYz+mcsXtHX/yooXkWmF2Eej6pwnNKrAhI\ntf9x5aXjx2RARFEnEGk6AUyFhOKujUBP7t5xv5gYMG1h58cnKm6aadx15XTMDtXFPEYbQigaWvNk\n1ATLfxGzHe8b52eSVAxwgsVERVtiZmGwx/uROmG2gnxA6PF18wXSjH8GhCpafz9dJWyLDr4SKlrq\ndWTUhvczL1LfiSxdvWDTjCzw08jVpk9eOL85L+liODqmK28xjzuPfXwDKnFkwciaFx6jR5yFS1Gw\n3UMocLbzX9oQSLQ62MyZghg3Gq2L+jItF6dQ7Pkq8QgweGLE7/zE4lPV7VGcsNGWQTacdhinXqZ4\n8xwyiMG/DUdqtEu+EA2XwRrvnUoXJRpl/AQwsoiObrNVtybXZMszCeepF1JdPJJQ9kKUgfrziare\nC5F3beSRqfyYi2jXpLUtZU5amWqjrapzj6Rj5J47fp3Qpu5/5I7HtTPEJCnE15W9rn6VaNQV8npf\nRUARUAQUAUWgAhDI7H1UAQBoFRUBRUARUAQUAUWgeAgo0SgetpqzIqAIKAKKgCJQ8Qgo0aj4R0AB\nUAQUAUVAEVAEiodA4SadF6+MmrMioAgoAkVDgIioTBMmfgxTMaOLZVX3xsyEIT+mbjJllynFLuJn\ndfPS9IpAOSCgRKMcWlHroAgoAjVCgPDrrFEBKWBKOL8EmXLh6KubKXkwlZiFxyAuBFNj/RWmubr1\nM6qbp6ZXBEodATWdlHoLavkVAUWgxgiwdgqLu6F1IFAd2ojbbrutxvkR64Z8WMiPSJF9+vSROAzE\nAVFRBCoVASUaldryWu+MCNBhEEyIjkel5giAIdEUkyoEWTr88MNlkTzKyMq5flhsyo92IiosMc65\nqBBfwV+mneBOaE3QdKgoApWKgBKNSm15rXdaBFhlk9Uc6WBcYB4IB+SD5aT5EF46V3HX8Btdi4R7\nufNESnQBgMgb34FPP/3ULF68WG5FuGPSlorccccdEkGRVS19odMl4qKLtOifq+1tSAEr8WIysQu2\niamDqI9OiKxLhEdCnPttR51YW4YVaX1h5V0i8rLkOEJUSbQb1YnO6+en24pAOSCgPhrl0Ipah5wQ\noCMnBHm2aH0sCMVS29jqWXSJES4L2rFsu1t6ns6fDinb2iGTJk0yzz33nHRmLHyGs6FdbVXK+8QT\nT5hnn31W7oH9HmLDMtPuPIs4XXjhhWbTTTeVBaoId8zaLXRwSRdC2bO+Cv4PDiOI2jvvvCNLhBNK\nnzUgqrMOTzHrjLMmhAB8MaWgiSAEP88L66uwngVt4ZZgp06QDNbzYCnz6JLnlJW6s7jauHHjUsKm\nF7MemrcikEQEVKORxFbRMhUFAVZYhGxkE1Z87dixY+i8h30dotG5c2cZqbI2Ayp28kMtnk7sUtnm\n+uuvlw6Va1n/wi2It2jRInPvvffKfSA0+++/v+Q1fPhwM3/+fMmSNVpQz0NI7BLVpkWLFiWxdsZP\nP/1kJkyYIITJkQwqBNGAdEEyWNEzCWYptEasocL6QbQDRAMtBJokBKLZo0cP0SrNnj1bjrkvtBe9\nevUyaG6iwsJbEEXICGQKp1AVRaBiEbCqWhVFoCIQ6Nu3b3DcccdlrKvVMARNmjRJSWNJBQsdBJMn\nTw6P2+XK5ZglLuExf8MuHS35TJ06NTz8+uuvB1atLvu2ww3sCpaBVdeH560TYmCXqg7s6FmO2RVA\nA7u4cni+VDYsUQqs5ieYM2dObJGnTJkSbLbZZgHp6lqsKSqwBC6YO3euFOWjjz4K7AJtgdVEpBTt\nvvvuCzp16pRyjB07sySwPh6B1YKF5ywxDazpJFi4cKEcAwfrpxKe1w1FoNIQUNNJxVLMyqv4mmuu\nKSPQTDVnJH7UUUelJMGsgR0/l1Vz3YXOAbBt27bmpZdektE7ZhNs/ggqeFZ1XLBggbvENGvWTMw6\nmB0GDBggx/FlQAvj/DTQEMQt4BRmkucGK/eiwcG8RB24H34hrB6MaYEl57MJ/g7ULVM8Cla8dHXK\nll8+5/EDAX/uVb9+fTGJEOMCcxgmKbQrlkCISYe6Y+5AY7XDDjuk3JYlwjF9cS1mFSc8GzwXlqCI\nDwf5oun68MMPQ3NK8+bNzUUXXSSXUB7Sgq0ltFnNeO4++qsIlDICSjRKufW07NVCYOeddzYPPfSQ\nxDVIdyHOmVGhs8V84QRnRpxB6WCcs6g7534xk9gRrSwVDnFAMIUMGzbMNG7cWAjHmDFjQvMM5/HB\noCPCTwNZe+21ZQl0TCYsSU2H2ahRo9CHQxIV+OvYY481dKrUGYfYDh06CLm48cYbZdpmz549s94R\n0wlL00eXzc56YRESTJw40UAgGjRoYE488UQxg9SrV09MO1ZbZA488EBjtVzGalnER6Z9+/ZS56jv\nyEYbbWS4DnPX0UcfHZaUOmJCoW3Ii+ehS5cu8nGJIDPOGfTnn382Z599tmFaLWWLOsq6a/RXESgn\nBJRolFNral3yRgC7upsxkC4zZn9AWPCrSNdR0NkyuuWXjoxO6sorrzS9e/eWGQzkjT+AE/LEro/W\nALs/goPhMcccY/bYYw8ZjUM4WrduLQSka9eu7tKC/eK4iE/CYYcdJgSIUTgdI06bdIxoABCO4WOB\nD8ree+9trEkqLAMzMxjNo7mBHKUTOt5iR8uEFOAbAblDI3PGGWcIkYMQ0DZgiUDoBg0alK6oOR1H\n6wRJhFS0a9cu7TVoQMCM5yYJRCxtQfWEIlBABJRoFBBMzSo5CKDmPuSQQ8y8efNSCmVto2b06NEy\nCnVOl9afooqqPOUib4frcQ5kFEtgpmydBaNbOh8EU8Jbb70l0yWjHRsaEzQd1qfDWD8NSU9nNHjw\n4NBUgjkC4gLJiSMaEBU0D5QRc8cuu+wiJgzJzPvChEOkyqjgsMoHoTw4RM6aNUv2Tz75ZPnlC1JE\n1EtiSbRs2TI8zgadLXgyUyPT7B7MM3T+mcT6NkhdXZqrrrrK4ESLdoLZHuCJc+nAgQPF3OHSuV9M\nHPvuu6/sgiuOnwTRQtBmOAEvPghmqWKapmjbXLRCrmz6qwiUAwJKNMqhFbUOVRBgaiIzNayTn3Qw\nJMCEgU8EnZQvUXu8f87fpqM677zzTLdu3cR8ARGgg8q3Y5o2bZqUE9U+HRGzM+hEL7nkEsn/hBNO\n8IshaVMO/L8dFzbbdZpxaThmnR3TnQqPQyL8wFPhCbuBZoMPHX9NBW1GJiJCvs6Pxd3D+Tm4fX6Z\nWkqI7zjx2xWNQ1y9MfFYR0/RzmASgyChQcJvQ0URUAQKg4ASjcLgqLkkDAE7q0FKhE+DEzQL+EH0\n69fPHcr5l/gPTIOkY+vfv79cR2AmRs0ciwoaD5xPswlaCMqEecLZ8bHhjx8/3syYMUNG79E81l13\n3egh2XcEIPZkDgfxKcFJEVMC9cWnBUHzAPEg/2zOoJAHnB/xiyCvdBofHF8zmVa4LxoePvkKjqDT\np083u+++u2TFPvXBdIWGCEJ6/vnnC8nAzDV27FiZShxHTPItC9ejzYFQVse5uBD31TwUgbpCQONo\n1BXyet9EIuD7TbgC4p9ArASCb9HpYp7ggwkG3wuEhbT22msv6bjYZ2RMx+yvm4HvAvl3796dJGJC\nwfkSx0OuJ09mt/izGriHEzpLTBOZfABc2pr8nnLKKRIMjHJSHic4thLEjHtnE5whqSNxRpgxk04I\nfoavRzEFXww7JdnYacWizXL3gtgxe4TOHu2Pb8LB3IbTLjNn4oQZKXFCm6Tz1/HTg8m5554rzw7P\nlYoiUAkIqEajElpZ6ygIZDMpkGjPPfc0BGYaMmRIiBrTEekcuH7o0KHyy0lMJphoEDoznCmffvpp\nIRxoUhgdE+UTFT5OiUSIhEwwUqajI5gTI//TTz9d8nAzWAgchRCJknyYLsnInnujZfBnPUjCAn0R\nBZNOGD+GVq1aSZnxZcCUQxmjMzG4rdN6+EVgGinanCgxYWYGDprOBAIZwzmS6Khu2q+fT77bmM2Y\nHYRJCcdP1lwhzDjEgqiuOL7SFsw4cUIaNEZRbQbmFYhRdMouvioPPvigOPBmMwVxD8wzaDQwHfkh\nzd399VcRKEcElGiUY6tqnWIRwByAWSCTQDRsoCyJXeFMH02bNpUOAhJBx4hDJ8IUUKf+Z6YIUyCx\n7zthVgqzSOjwEJwbMakgaC3OOeccIRzs4+BIh+hCXHMM8wKOlTi2QlSIMrrTTjtxqiiCyQDnWUxM\nEBxiXUC6fMdJd2M6a/wb4rQW1BGNAHj7HTN1w6xCZ+uEjr0YJIP8IQdgjhMtjruUB01LOmdMyOJE\nO+UUcugTDUxAzEbChyYqkDDq1aZNm+ip2H3amVgdhGdXUQQqBYEl7Cjtf+7WlVJjrWfFIsAMBUbR\n2VTczLDAEZJw4BALlVQE0BKgmYCAoGkZNWpUlWBeLERGx45DaykIJJIZNhBD58vhtFUsqIYWBqIR\nnWWD6QvtBE7CucqLL74ocTcgNOkcbnPNS9MpAqWAgPpolEIraRkLggAag2wkgxsRR4LpoYy2Vaoi\ngLkHjQbmFLQDmEiiHSbaHLQd0dVNq+ZW90fo+CFOffr0EZKBRsutdYK5BO0DsS/QxvjC84EZJBpJ\n1k8T3cYvAx8YCEs6R9noNbqvCJQ6AqrRKPUW1PIXBQFGscSrYMSqWo2aQ0y8D3xN4pxsa55r4a4k\nJggB2tC+OGFqK8HIMLMx0wjHVbcYnkvDrBWIFjNWVBQBRSAzAko0MuOjZxUBRaCMESA2Cr4cvjBb\nKG7Zdz8NGg+cRpm5oqIIKAKZEVCikRkfPasIKAKKgCKgCCgCeSCgPhp5gKeXKgKKgCKgCCgCikBm\nBJRoZMZHzyoCioAioAgoAopAHggo0cgDPL1UEVAEFAFFQBFQBDIjoAG7MuOjZxUBRUARqDYChDEn\nsBtB3JgiS9RXApMxlTabo2m1b6YXKAIJR0CdQRPeQHVRvKuvvlo88VkrgngIucSeqItyJv2ejz32\nmGGBNCJtEtGT9UN+/PHHsNgdO3aUiJU6cyGEpOQ3IBjE1fjzzz9lYba77rrLsMAfodoJ8c6Kt1dc\ncUVRI7yWPIhagbJDQE0nZdekNa8QQWIHDRpk6CBZRIzfm266qeYZVviVRJsk9DXxODp06CBBoM46\n6ywJ633ppZfKWhsamLe8HhLavG3btrJ+CiSTabBEGiWq6MCBA02TJk2yhsEvL0S0NoqAMUo09CkI\nESDSIYGqiJLIAloHH3ywEI4wQYltsN4GwZZYJ6QuhDgLbm0SommiPnf7RNIEYxZmq0uBBCVZIGLF\nXuU1n/rzfEEmWA8FYdl5t8IrEUDRZrVo0SK8BSv8ouFQUQQqCQH10aik1s5SVzpmVhr1hRdpqcp1\n111nfv31VwmsVKwVTx02LMhFlEgCQLnVSTE7udDcRJ5cf/31RYvBNXQ2LNrGCqK1LYQPZzEwFjzr\n1q2bWbRokYTFJnqnK29tlynufixkhumhXr16spquW5AuLm1dHCNM/YQJEwS/0047zRDoq0uXLvK8\nUR4CgbFeik800GypKAKVhkDtv+UqDeESq2+21U1LqTqsPoq9nFVWfWH0ifaGVT3jlj7308Zt01F/\n88034UqspGE11n322ce0a9cuJBpu9VeWEmcVVnxeXEfOSJdPXQhlP+OMM4TsoDEgvPauu+5aZcGw\nuiibf09WrWUdlb59+4btxFoht9xyS5gMc0R1BCJ4++23yyU4Z7L+TSaB7CxcuFDatX79+oaF+Vi1\ndpdddjHsUzZW1YXQImixnLz22muyfgq+GU78VWHdMf1VBModATWdlHsLV3j9IE5rrbWWoDBlyhSD\njwRLueOLwoJY1RG0FqzSCZkYPXp0yqV0JqwMG6ehgNjgBArxqGtTCYWGFI0YMULW8YBosMw52hdM\nO0kTtEAQAtabgayB4b///iuhv2lLFjuj489FnnnmGZkFQmffvn17IQFPPvlkxksnT54sC8c57RMa\noLlz58o1aCvQCq222mpV8qCM06dPTyEeVRLpAUWgQhBI3pulQoBPajVRB/sS9xL1z5fS9nbbbSed\nFsvA10Swv6MFqS4mzESAaGyxxRY1uW3Br6EsrFCKvwjE56+//jKXXXZZOCov+A0LlCH+Q8zQwdmy\nf//+5tRTT5XOnFlS2WTx4sVCEtGAQBY233xzM27cOCEbmAzTCffbfvvt5RoW2EOjMWzYsHTJw+Mv\nvfSSQfvCgm0qikClI6Cmk0p/Arz6owpmCiZObI0aNTLvvvuu6d69u5eitDcZyeOoh0mlJkLHzGfW\nrFnVuhw/A8hJUlaBXX755UWjA8FgZtHYsWPFmTEp5UsHLtgzg8M3P6A5gERkE8wwEBVWZXUC4ejd\nu7fMEDn88MPd4ZRftBas0lqd/wGEdObMmYIpWjO0XQ0aNEjJV3cUgUpCQDUaFdDajOIYzR955JHm\n2muvlRovWLDAHHbYYXK8devW5uGHH5bRPiYBXqo4t0E0Bg8eXAEIFaeKH3zwgdjvGdniEMiS48RX\niJN58+aZJ554QswCcedzOQZhaNWqlbQp7T1nzhy5bOTIkeEx2n/ZZZcV/xL8CaZNmyb7mE+SLmBI\nDAoIMfLmm2+aVVddVZ7rbGWHaOCf4ftQcA2+FRCQTIIjL7jxwbcnk/zwww/iQ8JMGWZtodlgxpGK\nIlDJCKhGowJaH4c1Xs7Mihg/frw59NBDxSmRlyyjWEZ2RDCkA+rZs6d8Cg0L/g1Mm81FGEEutdRS\nuSRNdBpmIfDJJjgbQuzwJyGwU019JVDtM2to+PDhBidGZpMgBFxjRD106FDpmF150GxgQqipTJw4\n0UyaNCmrQyX5E5ysEBExcd7lObr//vtl2vLdd9+dkykLguAcNv36gjVEMJ3g+IvWh/8Ns3KYVYLp\n5qSTTop9RtFcDRkyJF12elwRqEgElGhUQLPjIc8H/wI6IbQWOAFCMpiC50aIxYICZz3UyY0bN87p\nFkypJbBRpQid0wknnCA+HPlGCcW8wJRLyCXEAw0KbcxUXzrKQgqkyM2iyZZvTWb3xOXpfCaoF7NH\nCCgHQcv2bOGoy/P+3XffhTNYyL9Tp05xtwmPEXQLwoZTKnL55ZcLaeL/BOFQUQQUgewIKNHIjlHZ\npOjVq5dEqWT0zIuZEWmxSQbgMRODDrC2hdkmcbNAClEOVOT4BmDDz1fQLDALolCy5ZZbSofYr18/\n06NHDyGWhSYZlDVXjU2h6uXyIZ4GzzJTciEA+JpAADCjpBPSEY/DBdZy6fgfMP04naBt8qdHox3i\n87ENAqdEIx1qelwRSEVAiUYqHmW9x4sYZzp8BfDRICpks2bNasVJEQdTVNC5SDrHvFyu9dMwuwJn\nwWIIo+NcR/PFuH+2POlYnd8F0yz79OkjUzGzXVed82+//bYh5HYugtmkefPmuSRNm+bOO+80mDLw\nfUDQZBBnA7McHb8fGCtdJjNmzDAHHXRQutNVjjMVmvgn+LkgzNIpB7OeVEa/FIFaQkCJRi0BXde3\nwVSCLRp1M7EHiAtw7rnnikYjl+l6+ZafUSdaAF/oAHlx44zXtGlTcUbFObFly5ZVIpT61+W6zei1\nukQDcoIWBCKRSehsCtXh0DZ0oM6clem+uZz7+uuvJV4Izr44QDINlLU2mM5ZKBMG5cB8EW3TaPnw\nFYFkRDUJ0XTZ9jGn4RcBsXBEw8XPoI3BDnHhymk/39cFPyQ0XGgonHA9s1Buvvlmd0iuxzzjSORn\nn32Woinh+cDfoxBtj0aM+hC4zS9rWBjdUATKBAElGmXSkJmqwXLVeL9PnTpV1N3EUGC5ajQMdD44\nguIUWAz1uisXKmo6Cl+i+/65YmzTCdFJECuEDpIROb4jmBqYgsh5yolDJVqfqDBVkfTMIiCENwGg\ncLKEJOUjlAWyR7kgXzjl1kRYc+Pll1+WmSSUc7/99pPOFYdN2hpTWcOGDc2+++5bk+yrXIMGIRct\nQpULa3AAExXtAoGgHdDOEYANUoB5CB8MhJgXkGjIA1O0fSEoGe3KrCoCbeFQypTnbbfdNkwGKSLC\n66hRo+QYfijkB8HAVwNfl65du+bkABtmmmYDPxrai3pArlUUgXJFQIlGubasVy9IBquxIryUEezM\nbdq0kc4V50Gc4opJNOSmdfyFVoVZCpgU6JCY8siUUiJMIpxnRg7mAKZOEpHSF0bnHIcIMDrG9IR/\nBVE28xF8B+hoHnnkEYlyWVOiwawK4jdAgnBWdKPyzp07S10ff/xxKWahiEY+da7utUxLhRQTmZWl\n1tGQoBFguqtbxIw8WUvEkcHoPZjeDUFh5Vw0PZCUMWPGCNlwaQmw5RNHfD8gJxAMtE5oMmhv32/D\nXVvd3wEDBojp0g9RXt08NL0iUAoILGFVjv/TOZZCabWMNUKAF7Kb2kfnQ0cbVXvT2VU34mWNClOL\nF6FiP+ecc0Rj427LOh905JAKBPW5C1HOPjjRoTFFMboWCaPaOBOAfz3OtZmcC7lHnODDwvojkD7I\nS02FDpQPI39n/qHj9adw+uWt6X1q4zoifqJlwuzjhGeZmSPUh+fYD95FGsxGhA3HPIimLiqYWTCH\n0JZoKqJxNcAO05lP9iAYPDcQE+4JSXdCWSAjkDf8naor5J1Pe1f3fppeEagLBFSjUReo1/I9UTtH\nZ0dALEqlwykkXHQumcSZQ6Ikg2syzWrgPM6KjsCwn6vA9VnHg9E47ZKP+ATD5UPnyKfUxPlg+OWO\ne5b98xBFzBHEDIkTfCHiCIhL68iZ2+cXIlCslWOVZPhI63a5IrBkuVZM66UIgMDs2bNlBJwrGtdc\nc41ETM01vUvHyJaRctQvwJ3P9OucApkuWQgnw0z3KqVz+KygoQDXXAXz0xFHHJFr8rzSobnC/8l3\nMM0rQ71YEShTBJRolGnDarX+F40SvwRmLOQq2P9rItwDH46axFbAnIWDYlIWXatJ/YtxDXEyIGHg\nmqvgyInvUW3I/PnzxfEWh2CcS1UUAUUgHgH10YjHRY+WAQL4U2DjRx1ebBW1s+Pj2Bc1U5UBlHVS\nBTCFaOAvkUTTD2YafDoginGmtjoBTW+qCCQQASUaCWwULZIioAgoAoqAIlAuCKjppFxaUuuhCCgC\nioAioAgkEAElGglsFC2SIqAIKAKKgCJQLggo0SiXltR6KAKKgCKgCCgCCURAiUYCG0WLpAgoAoqA\nIqAIlAsCSjTKpSW1HoqAIqAIKAKKQAIRUKKRwEbRIikCioAioAgoAuWCgBKNcmlJrYcioAgoAoqA\nIpBABJRoJLBRtEiKgCKgCCgCikC5IKBEo1xaUuuhCCgCioAioAgkEAElGglsFC2SIqAIKAKKgCJQ\nLggo0SiXltR6KAKKgCKgCCgCCURAiUYCG0WLpAgoAoqAIqAIlAsCSjTKpSW1HoqAIqAIKAKKQAIR\nUKKRwEbRIikCioAioAgoAuWCgBKNcmlJrYcioAgoAoqAIpBABJRoJLBRtEiKgCKgCCgCikC5IPB/\n4GANvb4Qd+oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Image(filename='/Users/wy/Desktop/beales_function.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class GoldSearch(object):\n",
    "\n",
    "    def  __init__(self):\n",
    "        self.l = 10**-5\n",
    "        self.alpha = (math.sqrt(5)-1)/2.\n",
    "\n",
    "    def g_lambda(self, a, b):\n",
    "        return a+(1-self.alpha)*(b-a)\n",
    "\n",
    "    def g_mu(self, a, b):\n",
    "        return a+self.alpha*(b-a)\n",
    "\n",
    "    def goldSearch(self, a, b,lambda_k,mu_k,function,k = 1):\n",
    "        # step1\n",
    "        if (b - a) < self.l:\n",
    "            return (a+b)/2.\n",
    "\n",
    "        if function(lambda_k) > function(mu_k):\n",
    "            # step2\n",
    "            a = lambda_k\n",
    "            b = b\n",
    "            lambda_k = mu_k\n",
    "            mu_k = self.g_mu(a,b)\n",
    "            k = k+1\n",
    "            return self.goldSearch(a,b,lambda_k,mu_k,function,k)\n",
    "\n",
    "        elif function(lambda_k) <= function(mu_k):\n",
    "            # step3\n",
    "            a = a\n",
    "            b = mu_k\n",
    "            mu_k = lambda_k\n",
    "            lambda_k = self.g_lambda(a,b)\n",
    "            k = k+1\n",
    "            return self.goldSearch(a,b,lambda_k,mu_k,function,k)\n",
    "GoldSearch = GoldSearch()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gradient(f):\n",
    "    return [sp.lambdify((x1,x2), f.diff(x, 1), 'numpy') for x in [x1,x2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Fletcher_Reeves(f,xj):\n",
    "    lambda_j = sp.symbols('lambda_j')\n",
    "    e = 10**-5\n",
    "    sj = np.array(map(lambda fun : fun( xj[0],xj[1] ),gradient(f)))*(-1)\n",
    "    i = 1\n",
    "    while np.linalg.norm(sj) > e:\n",
    "        i = i+1\n",
    "        tmp = xj+lambda_j*sj\n",
    "        new_f = f.subs([(x1,tmp[0]),(x2,tmp[1])])\n",
    "        lambdaJ = GoldSearch.goldSearch(a,b,GoldSearch.g_lambda(a,b),GoldSearch.g_mu(a,b),sp.lambdify(lambda_j , new_f))\n",
    "        xj_1 = xj+lambdaJ*sj\n",
    "        sj_1 = np.array(map(lambda fun : fun( xj_1[0],xj_1[1] ),gradient(f)))*(-1)\n",
    "        beta_j = np.dot(sj_1.T,sj_1)/np.dot(sj.T,sj)\n",
    "        sj_1 = sj_1+beta_j*sj\n",
    "        sj = sj_1\n",
    "        xj = xj_1\n",
    "    return xj_1,i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 3.00000314  0.50000078]\n",
      "24\n",
      "0:00:00.916073\n"
     ]
    }
   ],
   "source": [
    "a = -5\n",
    "b = 5\n",
    "x1,x2 = sp.symbols('x1,x2')\n",
    "f = (1.5-x1*(1-x2))**2 + (2.25-x1*(1-x2**2))**2 + (2.625-x1*(1-x2**3))**2\n",
    "# 初始點\n",
    "xj = np.array([1,1])\n",
    "start = datetime.datetime.now()\n",
    "xj_1,i = Fletcher_Reeves(f,xj)\n",
    "end = datetime.datetime.now()\n",
    "print xj_1\n",
    "print i\n",
    "print end - start"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Fletcher_Reeves\n",
    "初始點 (1,1)  \n",
    "GoldSearch interval -5 ~ 5  \n",
    "e = 10**-5  \n",
    "number of iterations : 24  \n",
    "run time : 0.91s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def DFP(f,xi):\n",
    "    lambda_i = sp.symbols('lambda_i')\n",
    "    e = 10**-3\n",
    "    gradient_f = (np.array(map(lambda fun : fun( xi[0],xi[1] ),gradient(f)))).reshape(2,1)\n",
    "    Bi = np.identity(2)\n",
    "    i = 0\n",
    "    while abs(np.linalg.norm(gradient_f)) > e:\n",
    "        i = i+1\n",
    "        si = (np.dot(Bi,gradient_f)*(-1)).reshape(1,2)[0]\n",
    "        tmp = xi+lambda_i*si\n",
    "        new_f = f.subs([(x1,tmp[0]),(x2,tmp[1])])\n",
    "        lambdaI = GoldSearch.goldSearch(a,b,GoldSearch.g_lambda(a,b),GoldSearch.g_mu(a,b),sp.lambdify(lambda_i , new_f))\n",
    "        xi_1 = xi+lambdaI*si\n",
    "        gradient_f_1 = (np.array(map(lambda fun : fun( xi_1[0],xi_1[1] ),gradient(f)))).reshape(2,1)\n",
    "        if abs(np.linalg.norm(gradient_f_1)) > e:\n",
    "            gi = (gradient_f_1 - gradient_f).reshape(1,2)[0]\n",
    "            Mi = (np.dot(si.reshape(2,1),si.reshape(2,1).T))*lambdaI/np.dot(si.T,gi)\n",
    "            Ni = np.dot(np.dot(Bi,gi).reshape(2,1),np.dot(Bi,gi).T.reshape(1,2))*(-1)/np.dot(np.dot(gi.T,Bi),gi)\n",
    "            Bi = Bi+Mi+Ni\n",
    "            xi = xi_1\n",
    "            gradient_f = (np.array(map(lambda fun : fun( xi[0],xi[1] ),gradient(f)))).reshape(2,1)\n",
    "        else:\n",
    "            return xi_1,i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 3.00002223  0.49998904]\n",
      "8\n",
      "0:00:00.343704\n"
     ]
    }
   ],
   "source": [
    "a = -5\n",
    "b = 5\n",
    "x1,x2 = sp.symbols('x1,x2')\n",
    "f = (1.5-x1*(1-x2))**2 + (2.25-x1*(1-x2**2))**2 + (2.625-x1*(1-x2**3))**2\n",
    "xi = np.array([1,1])\n",
    "\n",
    "start = datetime.datetime.now()\n",
    "xi_1,i = DFP(f,xi)\n",
    "end = datetime.datetime.now()\n",
    "print xi_1\n",
    "print i\n",
    "print end - start"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# DFP\n",
    "初始點 (1,1)  \n",
    "GoldSearch interval -5 ~ 5  \n",
    "e = 10**-5  \n",
    "number of iterations : 8  \n",
    "run time : 0.34s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def BFGS(f,xi):\n",
    "    lambda_i = sp.symbols('lambda_i')\n",
    "    e = 10**-3\n",
    "    gradient_f = (np.array(map(lambda fun : fun( xi[0],xi[1] ),gradient(f)))).reshape(2,1)\n",
    "    Bi = np.identity(2)\n",
    "    i = 0\n",
    "    while abs(np.linalg.norm(gradient_f)) > e:\n",
    "        i = i+1\n",
    "        si = (np.dot(Bi,gradient_f)*(-1)).reshape(1,2)[0]\n",
    "        tmp = xi+lambda_i*si\n",
    "        new_f = f.subs([(x1,tmp[0]),(x2,tmp[1])])\n",
    "        lambdaI = GoldSearch.goldSearch(a,b,GoldSearch.g_lambda(a,b),GoldSearch.g_mu(a,b),sp.lambdify(lambda_i , new_f))\n",
    "        xi_1 = xi+lambdaI*si\n",
    "        gradient_f_1 = (np.array(map(lambda fun : fun( xi_1[0],xi_1[1] ),gradient(f)))).reshape(2,1)\n",
    "        if abs(np.linalg.norm(gradient_f_1)) > e:\n",
    "            gi = (gradient_f_1 - gradient_f).reshape(1,2)[0]\n",
    "            di = xi_1-xi\n",
    "            Mi = ((1 + np.dot(np.dot(gi.T,Bi),gi)/np.dot(di.T,gi))*np.dot(di.reshape(2,1),di.reshape(1,2)))/np.dot(di.T,gi)\n",
    "            Ni = np.dot(np.dot(di.reshape(2,1),gi.reshape(1,2)),Bi)*(-1)/np.dot(di.T,gi)\n",
    "            Qi = np.dot(np.dot(Bi,gi).reshape(2,1),di.reshape(1,2))*(-1)/np.dot(di.T,gi)\n",
    "            Bi = Bi+Mi+Ni+Qi\n",
    "            xi = xi_1\n",
    "            gradient_f = (np.array(map(lambda fun : fun( xi[0],xi[1] ),gradient(f)))).reshape(2,1)\n",
    "        else:\n",
    "            return xi_1,i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 3.00002222  0.49998909]\n",
      "8\n",
      "0:00:00.389444\n"
     ]
    }
   ],
   "source": [
    "a = -5\n",
    "b = 5\n",
    "x1,x2 = sp.symbols('x1,x2')\n",
    "f = (1.5-x1*(1-x2))**2 + (2.25-x1*(1-x2**2))**2 + (2.625-x1*(1-x2**3))**2\n",
    "xi = np.array([1,1])\n",
    "\n",
    "start = datetime.datetime.now()\n",
    "xi_1,i = BFGS(f,xi)\n",
    "end = datetime.datetime.now()\n",
    "print xi_1\n",
    "print i\n",
    "print end - start"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BFGS\n",
    "初始點 (1,1)  \n",
    "GoldSearch interval -5 ~ 5  \n",
    "e = 10**-5  \n",
    "number of iterations : 8  \n",
    "run time : 0.38s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.000000\n",
      "         Iterations: 56\n",
      "         Function evaluations: 107\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 3.00002489,  0.50000749])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import fmin\n",
    "\n",
    "def fun(X):\n",
    "    return (1.5-X[0]*(1-X[1]))**2 + (2.25-X[0]*(1-X[1]**2))**2 + (2.625-X[0]*(1-X[1]**3))**2\n",
    "\n",
    "fmin(fun,np.array([1,1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# scipy python做科學計算的lib\n",
    "出處 ： http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin.html#scipy.optimize.fmin  \n",
    "Minimize a function using the downhill simplex algorithm.  \n",
    "This algorithm only uses function values, not derivatives or second derivatives.  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}