diff --git a/examples/lanczos/Lanczos.ipynb b/examples/lanczos/Lanczos.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..1c3504f44a58a8b0a92d93b3e947134b1a0f46eb --- /dev/null +++ b/examples/lanczos/Lanczos.ipynb @@ -0,0 +1,133 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 34, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def step(v1,v0,A,beta):\n", + " w=np.dot(A,v1)\n", + " alpha=np.real(np.dot(np.conj(w),v1))\n", + " w-=alpha*v1+beta*v0\n", + " beta=np.linalg.norm(w)\n", + " w/=beta\n", + " return w,alpha,beta" + ] + }, + { + "cell_type": "code", + "execution_count": 59, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "N=800\n", + "M=200\n", + "\n", + "#mat=(np.random.rand(N,N)-0.5)+1j*(np.random.rand(N,N)-0.5)\n", + "#mat+=np.conj(mat.T)\n", + "\n", + "#evs=np.linalg.eigvalsh(mat)" + ] + }, + { + "cell_type": "code", + "execution_count": 60, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "vs=np.zeros((2,N),dtype=complex)\n", + "vs[1]=np.random.rand(N)\n", + "vs[1]/=np.linalg.norm(vs[1])\n", + "\n", + "alphas=np.zeros(M+1)\n", + "betas=np.zeros(M+1)\n", + "\n", + "for i in xrange(M):\n", + " vs[i%2],alphas[i],betas[i+1]=step(vs[(i+1)%2],vs[i%2],mat,betas[i])\n", + "alphas[-1]=np.real(np.dot(np.conj(vs[i%2]),np.dot(mat,vs[i%2])))" + ] + }, + { + "cell_type": "code", + "execution_count": 61, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3UAAAJKCAYAAACLVWzsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3X+MHeme1/dPzT133b5wzXp7Hcp3tCHk3CRIYyb0XWZp\nK0TnINApox5BUAISiiUnivgrcq9CSNhFGTz3L7IzioR7kMhm4Y+bNRKgQLLotOIewtAOG9nywBqa\nHgxSerhRxri3oW/vmmjH1hqe/GGfnjrVT51Tp6pO1fM89X5JpTvn3OPqp36dU9/6Ps/3iYwxAgAA\nAAD46Y22GwAAAAAAKI+gDgAAAAA8RlAHAAAAAB4jqAMAAAAAjxHUAQAAAIDHCOoAAAAAwGO9thsQ\nRRFzKgAAAADoNGNMVPbfOpGpM8awsLAssAwGA+u1NBgMGm3HeDxWv9+fakO/39d4PG5t37jYJp+W\n0WhkPbeSJGn8OCZJosFgoCRJNB6PNR6PNRqNNBgMNBqNaj2mt27dqq3d2fOv15t+ftr2+WhrYxzH\niuPYqXYWOS9u3brl/PXele8k2zXUxPWQ9521tra2tO8LFpZlLFW1nqkDsLhz585Z319ZWWm0HRsb\nG5Kkjz76SM+fP9fKyopu3rx5+n4bXGyTT168eGF9//nz5422Y2NjY+qYbW9v6yd/8id1cHBw+t7k\nv106ttnzb39/X8fHx1OfOTg40EcffdRau23XyNHRkR49ejT1ubbbaZM9L5IkmTonJPfa3eXvpCau\nh83NTR0cHEydB3Ec6+nTp1PntIvfF0CdCOqAEra3t7W1taUXL17o3Llz2tzcbPSHwvYj1u/3dfPm\nzcbaMJG9yXKBi23yhSsPDLK2tracv3mfSJ9/w+FQ9+7dO/OZpoPkrOw1MhwOrZ9ru53zuPIQYh7b\nd1LbvyNNKXI9fPHFF0qSpNS+WOQhxY0bN3TlypWg9ze6i6AOWJALGYP0j9jh4aHiOO7Mk18sl0sP\nDNKWffOeF9RU5WqQnOVLO7Py2r2/v6/hcOjszbsLvyN1K3IN5R2vzz//XJ999tnp60X3RdGHFMfH\nx6dBpe/7GzjDgf6jBvDJaDQyks4sSZK03TR4Zjwem9FoZAaDgRmNRmY8HrfdJGPMq3YlSWIGg4FJ\nksSJdvl63Y3HY9Pv96fa3O/3ndinabZ2xnFs1tbWnDs/02zt7vV6zu9vX8/nqmzH6/z587Xvi7z9\n27X9Db+8jolKx1Rk6oAF+dLdB25z+Um9bTxb2a5RdXE1gziPL+Opsu189uyZF2OSfBjDaNPV3xHb\n9fDkyRPt7++f+WyVfWH7vrB5+PCh0xldYBGRqaHaSqUGRJFpuw3AIpIk0ccff2x9/+7duy20CD7y\n5TyyBZ/9fl+3b99u/AZoe3vb+eAoFL6cn1l5Y7YGg4F2d3ebb1COvP27urrauTFfyzrX0t8XtmA/\nK45jXb58WRcuXOjU/oc7oiiSqTClAZk6YEG+ZgzgFl+e1LtUoMTFDGIZPhTI8OX8zPJlbKDtd6TX\n63VyzFde9cqjo6NKWbT094Xt4VTW4eGhDg8PT193Zf8jHAR1wIJ86U4Ft/ly8+nqzb3L3Vdn8aXd\nvpyfWcsKEOrma7fRZWii62/2b+zt7enk5GTmv3G1WqYPD4XQkioD8upYRKEUAB3kSwENVws6uNqu\neXxpty/np0260M/a2pqJ49j57RgMBtbz4uLFi04XqlmGJq6RooVU0outcFC22NWtW7eWWvzK5+sS\n84lCKQDgH18yvq52N3Y1gziPL+325fy0SXe7S5LEi0nV8zKjJycnneuO2cQ1UrSQSlq2e+be3t7p\n+xOffPKJXr58OfWZ9Di9q1ev6v79+6dZtuzrzc1NScrNxOV1h09nFMv8jXn/psg6Qz8vfUBQB6AT\nXOyy4sMk6a7e3PvaPdCndocwhtGXILpIkOFqd8C6NXGN5HX5TAdo89g+mw7oJp+ZFfTZgsDsutOB\n4T/4B//A2pb0WMwyf2Pev5n3+uDgQJ9++mmlwLCJ4LOpAHfRddZ2LVdJ89WxiO6XAJaMLivh8fWY\n0u5m5XWzW11dda5bY7rb6MWLF+d2B/Rh/5fR1rmW3v+rq6sLd8/s+pKdbzA7X+Sir+M4PtN1uuo6\nm/gbZdY5Ob+lat0vmdIA3nMxAwO3+FqeHbP5OsWBj+329RqyFabp9XpTWYa2puiYJW9/Z4U6BULb\n10iRaplAnZIk0c7OTqUpDZwI6l4FqgAAAADQRdXmqXujzqaUZQyL78t4vK1+/9uSotOl3/+2xuPt\npf7d0SiZ+puTJUmutb5PurqMx9sajRINBkONRsnSzwHOE44xy/KXUK6hwWBo3Y6LF3/EufN5PN5W\nklzTYDDU6uqPWtudXc6f/8bU6zi+rLW17zi3bWX2RdvfOenjkSTXNB5vn3nv1q33T1+vrX1HcXx5\n6nj0el+f+TqOL5/5N7OuO9u9V5m/Me/fzHudPe9YFl+S5JqqolAKatHWBMW+DIJvSttdUV2dg8vV\nCo4+cvUYS+2f/yEL5RryqcrkopNnS9KXX3459dpWsTFdjdGHa8SV75y8wlaz2pDtRrq+vq4HDx7k\nvp5cT7MKuKSvO1shq0X/RpF/U+T1nTt3ZnZzXvR1HMeSpgu6VF1nE3+jzDonx3RnZ0dVONH9su02\nhKjpm5vhcHj6g5g2GAy0u7u7tL/r6ziPZbD98DU9VsTl49H2GI1QuHqMXTj/Q2e7QfWtrHnR4Kjt\n89kmvf9tk5WXEcex80Geq985TfHlt2vRALat4HPZf6PMOifHNIqqdb9UlSordSyvmoA6tVE5qq0J\ndX2tyLYMLkxqnDeB7mAwaKwNWC5Xj7EL539Z2QmMffj+8vm7t0iVybfeesvpY2Lb/9nqg2UW2wTb\nbXP1Oweo2+uYqHRMRffLALXRFbKtrjmuzKHlQrcvF7qi+jQHF8px9Ri7cP6X4UrXskW11eW+DtnJ\nyW1ZoM8//1yfffbZ6WvXjklet7tsF7hFudhd09XvHMA1BHUBauPmps3gqu0JnF25KXPhhy+UcTfI\n5+oxduH8L8PX4MjXIDrLdj6fP3/+zNg0F4+J7bfvnXfeqTShdlY2yLNNMr3sQM/V7xzANQR1AWrr\n5qbt4KotrtyUufDD50rm1GcuZH1ncfUYu3D+l+FrcORrEJ1lO5+fPHmi/f39M599+PChhsOhk9fl\nRPZ3OD3OqY4g7+DgQB988MFU0Lvsh5iufucAzqnSd7OORYypq53PYx185FJ///RYkSRJOOae4dqt\nxsfz39exgCGfq3nHJIRtTV8ja2trJo7jyuPwJJnV1VWnxuEBPlLFMXVUvwyUL9WSQtD1ylyoD+dS\n9/hctTOEapg2PlfIXFTdmTypmfPX9R4NQBlVq1/S/TJQXe0K2QZfu33BPb52xUN5PnctmzePmmvF\nRYrKHpO9vT2dnJyc+VwI1+Wi3TVt4w2zDg4OdOPGDV25cmUpAVdI5xpQJzJ1QA3IjKIOZOrgq5DP\n3bxtW11dXVrg4gpbNnbRCpt1Z+5CPtfQbWTqAAeQGUUdyPrWiy5azQk5y2y7Lnu9no6Pj3Xv3j1J\n4WaK5lXYLDIJet2Zu7xzzYdCNsAyEdQB6CQXb/h97ornGrpoNSuUapg22evSFsi4OOXBsszrdmtT\nZwCcd66dnJwEH2QDs9D9EkDn+FycAsXQRatZtmsqjuPWJ65ehuFweBo8pA0GA+3u7jbfoJalu2gW\nydxJ1bquFg0k19bWdOnSJace3AGz0P0SABbkytyCWB6fuwO6mEWeJ5vNmhTZePTo0elnQsme5GWK\n9vf3O9n9r+nMXdFCNo8fPw7y/APyENQB6Byfb/hRjK/dAX3uNpq+uU+SZOqGWgrnwUmXx9jNU6Sr\nalaZ8yJ7rtmy8tnv82VX5QTa9kbbDQCApvl6w++y7e1tJUmi4XCoJEm0vb3dans2NzfV7/en3vOh\n6MysLLJPQn5wsrGxodu3bytJEg0GA62ururly5dTn/HxmNVlY2NDd+/e1e7urr73ve+duQ5tvvji\ni9LfH7ZrPe+7fBJ4f/zxx/rJn/zJ1r+ngDqRqQPQOVSZrJeL2SVfi86EEgyF/uAknSnKG2Pn2zFb\nhqKZu88//1yfffbZ6etFvj9s1/rR0dGZTHEWmTuEhqAOQOf4esPvKlfHKPo41UgowVDeg5P19XUl\nSeLVeMF5QjlmyzJvzJ1tQvNFvz9sk6g3XZUTaBtBHYLjQpEBF9qA2Xy84XdVKNklF4SSRbY9OLFN\nXB3CjXQox6wJtvPiyZMn2t/fP/PZKt8fZcf2kbmDzwjqEBQXuoG50AagSb5kKnx42BJSFjn74CRJ\nEiczulXlBbBbW1v68MMPnT3X2mI7L2xBXdXvj7bn01sWH77H0BJjTKvLqyYA9RiNRkbSmSVJkk61\nYWI8HpvRaGQGg4EZjUZmPB433gaEbzwem36/P3W+9/t9p843H9oYusFgYP1uHAwGbTetVpxri7Ht\nrziOzdraWq2/XePx2CRJYgaDgVldXbWei9lldXXVqd9Pzq2wvY6JSsdUZOoQFBe6gbnQBomMIZrj\nQ3bJ1XF/XeJLRrcqzrXFNDXHYdXM3d7eni5fvqwLFy6cZsgkNZo149zCLAR1CIoLNw0utEFy58uf\nriLV+LL/XB+j6MrDli6zjT2L41hHR0dBTdrNuba4puc4LDPm7vDwUIeHh6ev9/b2Tt9Pv5cO/K5e\nvar79++ffn9nXy96vuedWw8fPjy9hqr+DfiLoK4BvtyUhcCFAesutEFy48aCbGE17L/6uPKwpYxQ\nfkOaysi0zedzzQVN/XaVydylpYO59Hvp9z/55JOpOQyzrxfN/uWdWycnJ6cZxTJ/Y17wWSQ49WGd\nrra7tu/0Kn0361gU+Jg6+j83L91vPkmSVva1C21wYWyfC23wGfuvPr5+F/va7iJCPb9DPmZNaOu8\nKDPmru4ljmMTx3Hu+MK1tbUz/38df6PX61V67cs6XW335PtBqjamjqBuyUL90YL7XLix6EphhGVh\n/9XLhYctiwr5NyTk89vHc80VLvx22drgyhLHsfnOd75jBoOBuXjxYuvtYalnSZLESNWCusi8Cqxa\nE0WRebU9AAAAANBFkYwxUdl//UadTSnrVcYwzGU0SiRFZ5YkudZ625a5jMfb6ve/PbXN/f63NR5v\nt942Fs4DX9rpartYmltC/g2xnd9xfFlra9/RYDDUaJQEc66Px9sajZLgtqtLy3i8rSS5psFgqCS5\npvF4e+q9tbXvKI4vT53Pvd7XZ74uuwwGw9M2Za+huv4GS7NLklxTVRRKWTJXimY0zZXKi2iXD6Xu\nJXfPV1/2n89cL0IS8m9IVwqnUPAoDHkVftPvbW9vn5mE/sGDB9bXk/M9XVgljmNJ9iIsaZOiO3kT\n3y/6N3q93lRxlUVf+7JOV9s9+U7f2dlRFU50v2y7DcuWvci7cFM2HA5PKzGlDQYD7e7uNt8gYAbO\n126y3Wz3+33dvn3bqe/orvyGJEmijz/+2Pr+3bt3W2hRPULdLlRnu7alsw860gHCot9R8/7GvOCz\nyGtf1ulquyff6VFUrfslQR2Wgh8x+ITztZs47m4J9eFKqNuFZnTloQ5UOaij+yWWIuQuQwgP52s3\nuTCXY1mudxstI9T53ULdLjQjr9snkEVQh6VgLBB8wvlaPx+CDl9vtkMdo2V7uBLHsY6OjjQcDp09\nj+bhoRGAJtD9EmiADze4QF18GqvmQzuzQu42mu5qVsd4IlfQhQ7APIypQyEEFe3x9cYRKMunoMPH\nm+2ujNHy6TwCgKoYU4e5Qu2q4wtXy+UDy+LTWDUfx6v42m10UT6dR2XwsBVAnQjqOoCgol0u3Zhw\nE4Em+Bx0+HCNdGWMls/n0Tw8bAVQN4K6DnApqOgiV25MuIkIA0HH8vhyjXSlsI+v51ERPGwFUDeC\nug5wJahoS9s3wa7cmHATsZi2z5u8NhF0LI9P10i22+j29raSJHHqfK3Kdh6tr69ra2tLH374odfb\nycNWAHUjqOsAV4KKNrhwE+zKDa4rNxEuBktZLpw3Nj4HHT5w5RpZlKvnax3S51FI29n1h60A6kdQ\n1wGuBBVtcOUm2IUbXBduIny5KXPlvMnyNejwhQvXSBmunq91C2k7u/ywFcByENR1hAtBRRu4Cf6K\nCzcRvtyUuXre+Bp0+MLXya9dPV/rFtJ2hty1dBl86OEBtI2gDkHjJvgrLmRsfbkpc/W8cSEwD1n2\nGplMfv3o0aPTz7iYWc47X/f3950ORheVt53Pnj3zcjxhqF1L68a+AQoyxrS6vGoCsBzj8dj0+30j\n6XTp9/tmPB633bROGo1GU8disiRJ0nbTprh83ozHY5MkiRkMBiZJEifaFCqfz9der+fk+VuFbTvj\nODZxHHu/rb6ca21g36ArXsdEpWMqMnUImgvZKXzFl0yTy+eNr12pfew+5UtmOXu+7u/v6/j4eOoz\nLnZzXpTtujw6OprKpEp+bmveufbw4cOgsq1l+HIdAm0jqEPwfL0JDpHLwVJWF0rGN8XX7lOudsO1\nSZ+vw+FQ9+7dO/OZEAKE7HU5HA6tn/NtW/POtZOTk9Nj6cM1sww+XYdAmwjqgJb4mLmog49Btq9B\niSt8KZCT5WvhlC4FCKFsq+1cyzo4ONCNGzd05coVZ8+9ZfClhwfQNoI6NKqrgUwWQYJfXA5KfLim\nfO0+5WvhlKIBggvnb1WhbGv2XNvb29PJycmZzx0fH58Gq3t7e7p8+bIuXLjg7LVfB596eABSi7/L\nVQbk1bGIQimd4XLxiaa5MvB7PB6b0WhkBoOBGY1GrR0LV9qRZzAYWI/XYDBotV2+XFOunO9V+bQd\n6YI6Fy9etLb74sWLzl5ziwhxW/POtVlLv983t27dcvq7FAhdld9lVSyUQlCHxvh0Q7RsLgQJrgQE\nrrRjFlfPXVfbleXDMS4i77p1PWAoEiD4eDxsQtlW2zVTZDl//vzU6ziOzdra2un5SdAHLFeV3+Wq\nQd3M7pdRFP2YpP9Z0r/xulH/kzFmK4qif1/S/yjpN0j6vqT/1BjzLy3//vuSnkn6V5J+3RjzE7P+\nHsLmaxesZXBh4LcrXQpdaccseWM61tfXWy2e4ss1FUr3KV/HbxXtovjee+8535V3nlDGphWpaGrz\n5ZdfTr0+PDzU4eHh6etPPvlEL1++PH2d7cJ59epV3b9///QcyL7e3NyUJO/PE5/50OW+y9r8XZ43\npu7XJf1Xxpi/H0XRb5T096Io+puS/oKkP26M+TtRFP3nkv4bSX/a8u+NpKEx5ge1thpeciGQcYUL\nA79dCQhcaccstqBkfX1dd+7caXVcpE/XVAjVRH0dv1V0vNbjx4+dHy84T0hj0+ZNTl5GOqCT5gd9\ntiBw8u/S7y0SGNoCxUX3dzawKRJ81vF328Z4fPe1+bscvcr2FfxwFP1vkv6cpP/FGPPDr9/7MUl3\njTFvWT7/TyX9TmNM7uOlKIrMIm2Av2xfRv1+X7dv3+7kl9H29narmYskSfTxxx9b3797927n2rEo\nF9ptu6biOHbyBjXN5++C9HWbFzC89dZbevPNN529ecw7d21WV1edzmjNs8i2TvT7fV2/ft25ACB9\n7k2K9qSDq/Pnz5/J1LWh1+tNBYLzXtv296wAzPb9kV1nHMeSpoNP22eqBqPS7MCx7nXmZWzT16mL\n7XZlnU39jZ/92Z+dOvd++Id/WL/1t/7W3HNtcn5HUSRjTFTgMrMr2k9T0r8l6f+R9E1J/5ekP/j6\n/T8u6VnOv/lc0iNJf1fSH8v5zML9VeGv9ID2JEnoz98iV8Y52dqRHQfi4nniwrhIY6avqbW1NRPH\ncevHdB5fxgLOk7cd2XFNrh0D2zW3srIyd7yWD9dl1rLGprmw7dnf01u3bpXaVheW7P7u9Xq511CZ\nIjJFluzfnPc6juMz37eLrqPMOqtuR1vtdmGdTbU7W7Dpa1/72sx/Mzm/pQYKpUj6jXoVmP1Hr1//\ne5J2Xr/3pyX9i5x/d/n1/16S9Pcl/YeWzyzxKw/ALK4E2SEFJqurq63d+PkSLPlacCTLFjBkb05d\nOC/y2p6+9tfW1ha6cXT1urRJb+vq6urC22lbXAzyjJn/XbpogODSMvkey/v+YGHxeUmSxEjVgrq5\n3S+jKPq6pLGk/90Y82ct//+/K+nnjTG/a856bkn6/4wx/0PmfSPdSr0zfL0AAAAAQIh2Xy8T363U\n/XJe9ctI0l+U9I/SAV0URZeMMf88iqI3JP13kv685d9+Q9LXjDH/Moqi3yBpJOm7tr9jzPtl2w8g\nMMPh8LSAQdpgMNDu7m7zDZohPcYlb6xDk2PsiowfcmG8UNGiDz6O5yo6hsvFbStyPs/iw3hOqZmx\naT7ui0kBqAcPHuS+nhTzmrX/Fh1TV2Z/z1tnkTF1gEuSJNHOTrV1zMzURVH0uyX9n5L29Co9KEl/\nStK/I+m/fP36rxlj/tTrz39L0s8ZYzaiKPq3Jf3115/pSfpLxpg/Y/kbZl62EKiTq+WAXWiXC23I\nuyl28SY4LS8YvXjxot5+++1G2l00WMreRLVRpKRIwZE0nwqpLFql0MVtq6Paos+BTbaybVW+7Isy\nFg0Mba+z+7tIADb5TSgSfGY/U0cwWiRwXMY653G13S6s09V2T34D3n333UqZuoWqXy4DQR2a5GrV\nPRfa5UIb8tqR9wXo0k1R0SzZstu9aLA00WbQXDS7tba2pkuXLjn3QCarTMbLxYcW8zJai3IhS1xU\n3due5eJ3WJtmBYZ532NVe29UDUbnBY7LWOe86peuttuVdbra7kn188aqXy5redUENG08HpvRaOTc\nIO9lc7WQhAvtcqENE0UKG7R9zLKKVtlz4ZjOW5ouAlF032UrNLparCKtTPVFVwuQ1FFwxIeqkjbz\nCpCUWdbW1jr5O7wol36b2uZK1WoshyoWSiGo66Aufym4UobexXa50IZF2uVipcT0jV+2pHEb7V6k\nOmPbQV4dAYOrAUKZbXOtWmZW2akCfD2GWXUEedmHFF35HV5Ul+9ZbFypWo36EdRhYV1+6uXqtrvQ\nLhfasEi7XP+Bd6Xdy5jLatntLjuPmg/nRaiZu7qyV10O8iT3g/m2EMigCwjqsDBXMzJNcPWJnwvt\ncqENRdvlYvCZ5XK768iKLfvms4551Fy9SQ4xc2fM/MCmTJbYlwA3a96+KPKQwpdtBVCPqkEdhVI6\nKK8oQZOl19uUHRw9GaDaNhfa5UIb5rUrb9B8k1Umi/Kh3XVUOWyqAEyI1RhDqZZps6yqki4WlZkn\nuy+Ojo706NGjuf+uK7/LAFS5UApBXQe5UuUQ7nNhioMsV6pMLsrldtdR6a+Jm89lVGN04TwpUy3T\nl0qgWaEew0UVDebbfugDoDkEdSjF1YwM3OFq8F/0Zsi1J9w+tbvMjfdbb72lN998s9Ego44AwbXg\nqOh5srKyoufPn5++duHaLCPEY1jUosG8i5lmAPUhqFuAi1kHwFUud9Mt0q2x6hxGyxBCu/NuPl2b\n0LxoxsvF4KjMdkh+dkvMKhPkuXgMF1WmGy5BHhAW5qkryNUiEICrfCmok1dl0vXCEq5WG51nkWkS\n2tyWKiX3XTp3ylYCDeX3rWxRH5eOYVFFpkSZd8xv3brF3HeAp0ShlGJczjp0GdnTfG3vG1+uGdsT\n7l6vp5cvX56+dvHJfV731uvXr+v+/ftOXxPZ7ttPnjzR/v7+mc+1nXWcl/XJZnhsXDh3yhbZCCFz\nl2a7Znw5hosqMgbXJpsxz2bzrl696vz3C9BVZOoK8iXr0CVkT/O5sG9caENRRZ7mu5gBKzKHnKv7\nPM2XrGPZaRJc3I5Q5rpbVNlj6FvmblmTu/d6vanX2XkAyfQB7RGZumJ8yTp0Cccknyv7xseCOsPh\nUPfu3TvzfttZoyJcOe6LsmVQfBjv43MFwjLj7kLL3EndmRKijkqhRWR7OdiuY0lTvUiy2b8msoFt\n92QBlqFqpq5XZ2Nctrm5qYODgzNdnW7evNliq7rtxYsX1vfndaXpAlf2zcbGxpkfStd/TM+dO2d9\nf2VlpeGWLM6V476oyfHP3oCmuwhOvntdOley7c4Ljk5OTk4fFLiyHelrs2hgc3x87Nx2VFX0GKYd\nHBzoxo0bzge42e/feUFetutlGemATpIODw+n/sbe3t7p+xOffPLJ1L/Lvt7b25vbBVSaHSimj5Ht\nfM/+jTLB56Kvl/E3aHc390Vt30FV0nx1LGqo+6UxZ7tt0K2gXb502WqDq/vGhy6ZeW30oVuRq8d9\nUb5uR9Euby5uR5mCIi5uR1Vd6ppapPu2i4utC2gcx4W7iRY5v4uss+rrZfwN2t3NfTH5DpKqdb/s\nVFAHt/gQILTF1X3jy826r2PVXD3ui/J5DHORCoQXL150+uFA0cDG9e0oq0yA69uYuzzpbV9bW5t7\nc8nCwuLGkiSJkQIYU/dqewAAAACgi6qNqXujzqaU9SpjyMLS/DIeb6vf/7ak6HTp97+t8Xi79ba5\nuIxGydS+mixJcq31ts1aBoOhtd2DwbD1ts1afD0/fW13ke3w8fwvuh1ra9/RaJRoMBhqNEq8O16z\ntj9JrmkwGGp19Ufn7gcp0vnz3/D+/F103yTJNd269f7p67W17yiOL0/thzi+fOa9Xu/rM1+zsLDM\nX5LkmqpyIlPXdhvQXb5WG2xL3txqrleS8/U4+9puyc/KqTbp7djb29PJycmZz/hQWbXIdmTnfPPh\n2l5UmWqZEyFWD53Hdh1LmnpvfX1dDx48sL62FXSxVdiUNPMzNpPjYfsbRdZZ9fUy/gbt7ua+mHzX\nvvvuu6qSqSOoQ6f5XP6+LT7erPsajHJ+uiUvyPbtZn+Ria19eICwqDLTQWT5MGWHK7K/GdkgcF6g\naAvast/fVYPPMq+X8Tdodzf3xeQ+quqUBgR16DSfMyGucH2KgwnbjcWy51KqivPTLbaHA3lPXF07\nl9Js25HN0k24OEdfnWz7oszUAD4cd5/5+DARWBTz1AEV+DR/oYvBk+2GyNX5r+bN6+Viu306P+dx\n8fxdVJESz/0jAAAgAElEQVQ50Q4ODvTRRx85vW3Z7VhZWdHR0dHUvIITLs7RVyfbvlhfX9edO3cW\n6qLpy9x3vrLNmQogo0rpzDqWV02AC8bjsfPzeC2DD/MXulrq3pcpDrJ8arcP5+c8rp6/Vfk8dUNW\n0WkQQin9P0+ZaRHSS3petdD3FYB6vI6JSsdUZOpa4OITa18yF8vgwxPAra2tM0+NXcgIvHjxwvq+\nrRuXS3xqtw/n5zyunr9VnTt3zvr+yspKwy2pLpuxyiumcnx8HHTmbmJeZn+ew8PDqTFgIe8rAG4g\nqGuYq8FTqDddoXA1CPH1ptbXdvvK1fO3qrzusevr60qSxKkHd0WkA5kixVS68huRDXhthTvm6cq+\nAtAegrqGuRo8hXrTFQpXgxBfx3z52m5fuXr+VlV0PJYLD+4WZbtGbB4+fKjhcOhV8FpGNmNepoJm\nV/YVgHYQ1DXM1eAp1JuuULgahNhuan2oSuZruyU3u2/P4+r5W4fszX6SJE4+uFtUkaIwUviFVPKU\n6Z7Z1X0FoBkEdQ1zNXgK+aYrBC4HIb6O+fKx3a52357H5fO3bq4+uCtj0cClqxUgy3TP9DHQB+A2\n5qlrmMuTIDMPDLrMhwwY89a5L5QJym3SvxF5hVTSXPlta0ORffXWW2/pzTffdPo7B0BzmHzcQwRP\ngFtcftiSNhwOT7tvpQ0GA+3u7jbfIJwRygTl8xQppCKFEcxWlbevspOch3BeACiPycc95GO3r6b4\nkC1BeFwtYJTlavdtfCWUCcrnKVpIpStTIMxi21fZgE4K47wA0B6COjjDlfFCBJbV+Lj/fBkHxdhX\nP6Qf3OVlV107txZVtJBK2sHBgd577z3vvh+qso0pffLkifb39898lgqZAMoiqIMzXMiWuBJY+srX\n/edLBqxLBUdC4cu5VUaZCpCPHz/Wo0ePTl/78P1QB1uVVFtQR4VMAGW90XYDgAkXsiWzAkvM5+v+\n29zcVL/fn3rP1QzYxsaG7t69q93dXd29e9fbG77t7W0lSaLhcKgkSbS9vd12k5bCp3Orio2NDd2+\nfVtJkmgwGGh1ddX6uez3+aRiZujnQZbtvMjy4bsTgDvI1MEZLjzRdiGwLMrFbo4+7b80MmDN8jWj\nW0aXzq15mbuVlRXrd0EXx91lz4u8Cpmuf3cCcAdBHZzhwnghFwLLIly9KfZl/9lQwKg5LnS1blL2\n3JpkKV16IFM3WzB7dHQ01fXSJuTzICt9XuRVyNzf32eMHYBCCOrgDBeeaLsQWBbh6k2xL/sP7fI1\no1sHVx/ILIMtmC0y7q6LxUJs3529Xq+TWUwA5RDUwSltZ0tcCCyLcPWm2Jf9h3b5nNGtytUHMk0o\nWjGzi8VCik6FcePGjc7P+wfAjsnHAQ/lddVJkkR3795toUXhcXHMYih8mex9GZhA/itFM3ddnMA8\n7zxJ68o1A3QFk48DHUQ3x+XqUhe5NnQ5o9vlLGVW0WIhXeyCmHeepJG5A5BGpg7w1Pb2didvipvg\nUyaUjKJfupylnCfvurN9zrXrsG5Fs5hp/X5f169f1/379/k+ADxEpg7oqLbHH4bM1TGLWWQU/dPl\nLOU8th4INl0opFJ0/GHawcGBPvjgA3355Zen7+3t7eny5cu6cOFC0PsLAJk6ADjDl0ydL+3EbGRb\nv5LugVAkkOlKlrNM5s4mjmOCPMBRZOoAoGa+jFn0JaOIfGRbp82bwDyLyqGLOTw81OHh4enrbCbv\n6tWrdN8EPEWmDgAsfBizSKbOfxzD2dLXYV4hlbfeektvvvlmpwIRW8B7/vz5qa6XZfR6Pb18+fL0\ntS2zJ2kqs5wNBJsIDMluI0RVM3UEdcjFlybgNopu+I8pDorLC4CzwUxXroHsg6f19XXduXOnchfN\nWeI4lqSpbF82EJwXGNqCPql4oPjs2TM9ffp0qg11BJ+Lvl7G36Dd3dwXk/trgjosBTeL0whwy2Pf\nLZcPGUXkI1NX3CLZqa7uv/T3gS34cYEt6JMWCxTnKbPORV8v42/Q7m7ui8n99bvvvlspqJMxptXl\nVRPgmtFoZCSdWZIkabtpjRuPx6bf70/th36/b8bjcdtNO2M8HpvRaGQGg4EZjUatt9GnfYd2uXbu\nNoVrZDHj8dgkSWIGg4FJksRcuXLF+lt18eLFzp1LNun9tba2ZuI4tu4vFhaWdpckSYwkYyrEVE5k\n6l5tDwAAAAB0UbXul2/U2ZSyjGFxbRmNEknRmSVJrrXetqaXwWBo3ReDwbD1trl+zHzZd0WW8Xhb\no1GiwWCo0SjReLzdeptCWVw8d1n8WMbjbfX737aeP5xLxfZfklzTYDBUklzTrVvvn75eW/uO4vjy\n1H6M48tn3uv1vj7zNQsLy/wlSa6pKqY0gJUvJd2bcO7cOev7KysrDbdkNhfL2/uy7+ah7PxyuXju\ntolxqMVlS/3nVcjswoTlZaSnkLCxjdmVdKZAy4MHD6yvbeP66hiDFMexvvWtb+mb3/xmbuGURdfp\n63gsF/6GL+t0td2T++udnR1VQVAHq+wPZZcLMPgS4LoYQPmy7+bZ2to6U1GuK/NjNcHFc7ctPEBY\nXDowySs8c3JyclpllP1ZXF7Qt8i+s1XpTAeBiwaKtvuRqsFnmdfL+Bu0u5v7oq77ayfG1LXdBrjL\nlSfWPlQYdLViqQ/7bh6fys67cs0swtVztw1Uw6ymyITlEvsTgHuqTmlApg7OcumJ9bwuKi5wNbvq\nw76bx5dMkkvXzCJcPXfbQFfUauiOCaCryNTBWTyxhit8ySRxzfiPY1ivvP2Z5uK1DKB7yNQhWDyx\nhit8ySRxzfgvlHGorrDtz6yDgwPduHFDV65cIXMHwFsEdXCWL13e0A0+dCMN6ZrxcWxgHXx5gOCL\not0xj4+PKaQCwGt0v4SzfOnyBrgilGsmlO2Ae4p0x5SktbU1Xbp0qXMPFQC0h+6XCBZPrIHFhHLN\nMIXEtK5mLZehSHdMSXr8+LEePXp0+prsHQDXkakDgJK42V4On6aQWDaylvVLT7Oyv7+v4+PjQv9u\ndXWVcXcAloZMHQC0wJfpA3wMPEMaG1gVWcv6pcfH2q7jlZUVa3Ehxt0BcBlBHVCCjzfKqJcPN9u+\nBJ5ZVID8ChVNl8vWZfno6Giq66UNFTMBuIagDliQrzfKqJcPN9s+BJ42eWMDpVeFLrr0MIWs5fJl\nK9vavuNt0pm7vb09Xb58WRcuXOjMuQnALQR1wIJ8vVFGvXy42fYh8MxT5Ea7Cw9TyFo2L/tQoci4\nu8PDQx0eHp6+7sK5CcAtBHXAgny+UUZ9fLjZ9iHwLKqrD1PIWrZj3ri7eeieCaBpBHXAgny6UXZx\n7J+LbSrDh+kD8gLP9fV17wKCLj9MIWvZrjKZO4nCKgCaRVAHLMiHDI3k5o2fi22qInuz7Rpb4Lm+\nvq47d+54dwzyHqbs7+9rOBx6E5zWoatZyzbVkbnj+ABYJuapA0pIz3PkYoZGetU16+OPP7a+f/fu\n3RZa5Gab6uRDFtLXY2C7ke71enr58uXp667M38Y8fu1L/wY8e/ZMT58+nRpTZ3Px4kW9/fbbzn43\nAGgX89QBLXA9QyO52V3NxTbVxZcspK/HoEgXuK5kQ3zqAh4qW5fYed0zT05O6I4JYGneaLsBAJbD\nxRs/F9tUl1ld4lwyrxtjkiTa3t5uuFXFbGxs6O7du9rd3dWVK1esn3n48KHz21HV5uam+v3+1Htx\nHOvo6Cj4bXdV+tz83ve+d+b4ZB0cHOi9995TkiQcMwC1IFMHBMrFsX8utqkuvmTAbMeg1+t5V9Qh\nLzjtQjYkm7WcdP9LT5gd6rb7IHt89vb2dHJycuZzjx8/5pgBqA1j6oCAuTj2z8U21cGnsWpFuoqt\nrq46XY69aLEK17ejDnnnXhe23Qd5x8eGYwZ0V9UxdQR1DvChuAKA2WxBhg+FO/KKbqS5uh3p4DQv\nG5Lm6nZU5fMx7ALbd8PKysrcLD7HDOgWgjrP+XojCOAsH7OQRbMIrmcQQtmOMrq87b7IfjccHR1N\ndb3MwzEDuoOgznO+ddkio2jHvoGND+dFmTm3XHzwVHY7rl+/rvv37zt9jOYps+1xHOvy5cu6cOGC\nt9vtM44ZgCymNPCcL8UVXCnX7uJNsiv7Bm7x5bwoMlVAlotTB5Tdjg8++EBffvnl6Xt7e3ve3TiX\n2fbDw8OpedUODg706aefeh/g+qKOY+bjuQpgeWZm6qIoWpF0T9I5ST8k6ReMMT8dRdGPSPorkn6L\npO9L+iPGmF+x/Ptrkv6spK9J+gvGmJ+xfIZMnQeZOhfa6WpXVRf2Ddzj63lRNIPg+kTKZTIhNj5m\nR8pu+/nz56cC3Oy2X716laBvSeo4X22Z5+wxW/S1q8fYxQe8QFVLzdQZY55HUfR7jDG/FkVRT9Iv\nRlH0uyX9AUl/0xjzQRRFf1LST71e0g37mqQ/J+n3SXoi6dMoiv6GMeZx2caGyJcS7y5kFGfNA9bm\nl7kL+wbu8fW8KJpBcH3qgDKZEBsfsyNltz0d0Elnt/2TTz7Ry5cvT19n90WRAEESN+MWdZyvtsxz\n9pgt+tp2vkvTx7COwHGRdU6m8Jh3Xdbdzqrtbmudvra7S/uiru/BwmPqoij6hl5l7f4zSX9N0sAY\n88tRFMWSdo0xvy3z+auSbhljrr1+/VOSZIz57zOf63SmTvKjuIILWYe8Cm+DwUC7u7uNtMHGhX0D\n94RyXhTNIKytrenSpUvO3qzbtiObmSrDlsmT3Apc6spaztPr9aYCguzrOI4laepmfF42sKuBYlPH\nbFG2YzjvuJc5L+b9mzbaWUe721inr+3u0r6Y9Dp79913l1soJYqiNyT9kqS+pD9vjPlvoyg6McZc\nfP3/R5J+MHmd+nf/iaTEGPPHXr++Lul3GWNuZj7X+aDOBy50fXT1JtmFfeO7ELvShHReFJk6IFui\n3cVtzT5AW19f1507d2q9cS4SuLRxfqe33ZbpqCPArUMdN1ChBIrzjhmAcCRJop2dnWaqX0ZR9Jsk\n7Uj6aUl/PR3ERVH0A2PMj2Q+/x9LulYkqJNupd4Zvl4AAAAAIES7r5eJ7zZT/dIY86tRFG1L+nFJ\nvxxFUWyMOYyi6LKkI8s/eSLpx1Kvf0zSF/Z1v1+8xeg0H7qqYjGuZmDrFkrmbpGJlF0vpmLTRnbE\nhUyeNH/bF+0CF7IyGcW6u+kWyTxX7UoGoBmvMnXV1jGv+uWPSnppjPmVKIrO61Wm7ruSEknHxpif\neT1W7oeNMdlCKT1J/0TS75X0zyQ9lPRHs4VS6H4JNMfFbo6ujpWsW0jBa5mJlH0MYKV2gjxX9pUt\naHjw4EHhoK9IV8muWlY33VnHbNHXtmPsyhikOI71rW99S9/85jcba6er47Fc+Bu+rNPVdjcypi6K\not8u6XuS3ni9/Lwx5sPXUxr8VUn/plJTGkRR9C1JP2eM2Xj973+/vprS4C8aY/6M5W8Q1AENcDVT\nFFKwM0vIwWvRog4hHNN5QV5dgcvq6qquXLnizMMXm0UDiElV50WygV0OFF3I4Np6x0iqLXAsu85s\nL50m2llHu9tYp6/t7tK+mJzPVac0KDymblkI6oBmuBo8uRps1s3V/V+XIsVUfOyOOc+8m8k6snsh\nXg8TVTNLXQoUXQjyACwPQR2AQlzOFHVhrGRXglcpP4BNC3XbberowhlK8N8EFwPFZSDIA8JCUAc4\nyMWxa6FninzQheBV6lZ3zDLKBHkhZjl9UiVQXGY33bR+v6/r168vZVJjAMtHUAc4xtWMjKvtsnEx\nKF6GkLezSHdMF7LELkjvq/39fR0fH8/8vKvXLfI10U1XOjvfYKjZvJC/O9FdBHWAY1zOiPmQKfIp\n+KyiK9spuX1NuIYsZ3c1UWm1jonZ2/5+6tJ3J7qFoA5wjMtj13zQlQCgK9sp2W/CQs0g1IEsJ6Rm\ngrxFxwbaunguGhjaAkVpes6+Wf8mL5udrRi7yDqbaHdb6/S13V3aF5PfP4I6wDFdullfhq4ExV3Z\nzol5N6g8abfL+z5ZW1vTpUuXnMqgYLnmXUPZrpfLkv07bczrNU+X5jhz4W/4sk5X213XPHW9sv8Q\ngN3m5qYODg7OdA2ZjKHAbOfOnbO+v7Ky0nBLlqsr2zmxsbFxGnQkSXJmsvKDgwN99NFHBCYZtu+T\nOI719OnTqX04+f/Zf+FKX0OSvXjLnTt35nbdrSobOGaDrUVf27KP8/7NPGXW2US721inr+3u0r6Y\n/P5V5USmTiJTBwAAAKCrqnW/dCJTR+9LAGm2J9CuDdavQ5En7SF2S6SLcjVd67qL8uZ12Vy0a1lT\nXTyLmIyhK1IxFnBdkiTa2am2DieCOgBIS3czshXZCKWrWbY7VZIkZ7pOhdgtkS7K1XSt6y7KK9Jl\nc5GJ2W0Pntocg7SxsZFbiGnRdfo6HsuFv+HLOl1t9+T3b6diVEdQB8BpW1tbnQh0JOnFixfW958/\nf95wS5Zrctxcn17DVXnj7I6OjjQcDoPKZqNe2SCvjHfeeadyYJh+nZ2zr+i/mWxH3vfJoutsot1t\nrNPXdndpX9T1++fEmLq22wC/MQlpeT7suy51Netyt0QfzkWXUE0UAMJSdUoDMnXwWshd85bNl32X\n19Vsf38/uKxEVzMwvpyLLqGaKAAgjaAOXutS17y6+bLvbIFOr9fT8fHxaQYvlAAg241okoEJvXS9\nL+eiq7rSbRcAkO+NthsAVOHLzcz29raSJNFwOFSSJNre3m67Sd7su42NDd2+fVtJkmgwGGh1dXVp\nc7y4YGNjQ3fv3tXu7q4uXbp0Zs6bkLZ1wpdz0VUUTkHXuPibCrSNTB285sPNjKtdy3zYdxPprmZ5\nY+xCDAC6Euz4dC66iGqi6BJXf1OBtpGpg9c2NzfV7/en3nPtZmZW17I2+bDvbOaNsQvpqW1Xgh1f\nz0VXZLPZSZLo+vXr2traCu6aAFz9TQXaRqYOXvOhNLqr2RYf9p1Nl8bYdSUD4+u56JKuzO0IuPqb\nCrSNKQ2AJetymfplSZdz39/f1/Hx8ZnPhLJ/bZMF379/P/jS/0xxUB7fOQgZ5zdCxZQGgOO6km1p\nUpExdg8fPgxiGoAuZmC6sp3LQiYDIeM3FbDrfFDH02AsG13Llitv3NnJyUlw3TG7Uvq/K9u5LF0Z\ni4lu4jcVsOt0UMfTYDQlnW1xmY8POWxPbbNCCQi6koHpynYuC5kMhM6X31SgSZ0O6ngaXI2PAQDy\n+fqQI/vUdm9vTycnJ2c+98UXXyhJEq/P165kYLqynctiy2Ssr69ra2tLH374obfnPwAgX6eDOp4G\nl+drAIB8Pj/kSD+1zRtE//nnn+uzzz47fe3j+dqVDExXtnOZujgWEwC6rNNBHU+Dy/M5AIBdKA85\nbAHB+fPn9eWXX059zsfztSsZmLwxM5K8z7a2ge9rAAhfp4M6ngaXF0oA0CbXuq+G8pDDFhA8efJE\n+/v7Zz7rY4XMrmRgsmNmQt7WZeP7GgDC1+mgjgpK5YUSALTFxRvUkB5yZAOCJEmsQZ3vFTK7lIHp\n0rbWLe/7en9/37uHGgAAu04HdRIVlMoKKQBog4s3qCE/5ChaIfPGjRu6cuWKNze5XcrAdGlb62Y7\n/3u9no6Pj71+qAEA+ErngzqUE3IA0ARXb1BtXd5CGMNUtEKmbze5XcrA0DugvOz5v7+/r+Pj46nP\ntP1QCQBQTWSMabcBUWTabgPQtLwKjUmS6O7duy206CxbF9F+v6/bt297f+OXt/9tn3PleNjYjlGv\n19PLly9PX4dyzGzbGsexLl++rAsXLgQVwC7bcDg8fXiR9tZbb+nNN9/0/iEOAPgoiiIZY6Ky/55M\nnYNcK6CB+vnQfdXFLqJ1KdIdU3K/kEqXMjDZbX327JmePn2qR48enX7Gh+yqC/KyniFM+wEAXUVQ\n5xgXC2igfj50X3W1i2gdigRDkh+FVNJdZvMyMCEcM+nsfITpgE4KJ4BdtpCn/QCArqL7pWN86Jbn\nOjKd9ejSuWh7mGLj+rbnHbPV1VWvCsAUkRfAXrx4UW+//XZQ27oM29vbhab9YH8CQDPofhmYkLMj\nTSDTWR8fuojWpWghFde7Y3apymFeF0IfsqsuKDPtx97eHmMYAcBRZOoc06XsyDKw/+qVfZo/CehC\nz4QWKaTiagGS9DHL61a6tramS5cueX0Mi2ZXQ8xSLkPR/ZlGoRoAqA+ZusB0KTuyDD5nOl3sNmqb\n4qALmdCi89q5ON6oyBi7x48fe19gpMw0FWSa8hXdn2mHh4c6PDw8fZ3dv1evXtX9+/envtOk6YdC\nts9wTABgcWTqHGTLjvAjV4yvmTpfpg/wdf+Wkb4O825wXS8BX3TqBsn/jNYi2zrR7/d1/fp1ggqL\nMvszKzu9RhzHkjQVCNo+MyswLBMoLvp6Geuk3f6vk3b7v05X2z353amaqZMxptXlVROAeozHY9Pv\n942k06Xf75vxeNx202YajUZTbZ4sSZK03bQpg8HA2s7BYNB205Yq7/icP3/e6XPNdj2srKxYtyW9\nxHFs1tbWzGAwMKPRyKltymPb1iJL9hhmt/3WrVtmNBrlvvZh35RRdn/WvfR6vZmv4zg2cRwv9G/a\nWCft9n+dtNv/dbra7sm9gyRjKsRUZOoQHB8znXnd5AaDgXZ3d5tvUI4uZerSbJlUWwl4yb19kb0e\njo6OzkwFMI8vY6eKjCdcVDaTNC+zVMeTX1f2b3p/TuYFTGfZAAD1SJJEOzs7lTJ1TgR1rwJVAAAA\nAOiiAAqlkKhD1/kypk6a//Te1XbXLfQKmYtwdWzavHM1L9vqmiJj0+aNRVvG8Zi3f8uMqQOALnqV\nqau2DicydW23AXCBj91Gu9odUwpnwvIypextsgGSi102s9fY+vq67ty5U3nbXVSkm+gyg7zJ/n3w\n4MGZKVHyPlMkMCwTKC76ehnrpN3+r5N2+79OV9s9eQD87rvv+t/9su02ACjHl7GAy1KkQubFixf1\n9ttvOxPY2DQxdsrFIE9aPNsUSmbJh+NhCwwXDRTLvF7GOmm3/+uk3f6v09V2Tx7iV61+SVAHoLQu\nZ+qyfO6OmdVEkOfqvlgkqLDtmzqe/LbB1S60ANAVBHVAoFycjDzLp7GAyxZKd0ybZY1N831uPMne\nbVoq/+S37Ni0OmSPY1evZQCoU9H7OYI6IEA+BUu2zEZXn/YX6Y4ZQtfUZYxNc/X8bsOiY9OKBIJl\nhRB4A0BbFrmfI6gDAuRrt0afgtFlyzuGod4k19FlM9R904RZgWBdXWhdHYcHAK5a5H6ualDnxJQG\nAKa9ePHC+v7z588bbslitra2zmRrDg4O9NFHH3Xu5m9zc1MHBwdT+6PX6+n4+Pi0uMzk/wth32xs\nbExtR5kgL9R904Ts/s+qowvt4eHh1L/hGAHAbE3ezxHUAQ46d+6c9f2VlZWGW7IYX4PRZZjc6M6a\nDy7kgHdWkFdkbryQ900bZh2Psl1oDw4OdOPGDbKrAJCjyfs5gjrAQbYsT7/fPx1L4ypfg9FlSd9I\n503/0JWAN70vihaVefjwoYbDIQHDEtgye++8887Ck9Kns6t7e3t0zwSAlCbv5wjqAAdlszy+TEbu\nazDahLyAd39/v3OBS5EspiSdnJzQHbNBZQLvtGz3TFuQJ8n5qr4AUJcm7+colAKgVrYS79y02W+S\nsxUKu1pUJuTpIHxW93yFtqkYssVXrl692tnquQC6jeqXAOCJIuPKuhq4dGU6CJ8tOi6yDNucfPOC\nPonsHwD/EdQBHeHDZOQ2vrZ72fLG2BG4+DulR5eU6Z5ZhyITsc8LBOe9tgWKi66jjb9Bu9kXIbW7\nS/ticl9UNaiTMabV5VUTgOUZj8dmNBqZwWBgRqORGY/HbTdpYePx2PT7fSPpdOn3+85vi6/tbsJo\nNJraL5MlSZK2m9a6vPPm1q1b3l/LIRmPxyZJEjMYDMza2pqJ49h6Tre99Hq9hV7HcXxmWxZdRxt/\ng3azL0Jqd5f2xeS+SJIxFWIqMnUIWiiTYfuaufC13U3IOzevX7/OmCIVK7nv47Ucsnlj8GxZNQDA\nq/uinZ0d/7tfvgpUAQAAAKCLqnW/fKPOppRlDEsIy3i8rdEo0WAw1GiUaDzebr1Ng8FQUnRmGQyG\nrbdtkWU0SqzbkSTXWm9biO1mf7m3hHIts0wv4/G2kuSaBoOhkuSabt16//T12tp3FMeXp453r/f1\nqddxfPnMZ1hYWFh8W5LkmqpinjrUwtaVzIV5pUKZDNvX+d98bXdbXrx4YX2/KxOUzxLKtYxptknQ\n02zdcB88eDA1ZYqkmd0+s8VV5r22dRNddB1t/A3azb4Iqd1d2heT+6KdnR1VQVCHWmxtbZ2pgnZw\ncKCPPvqo1aAulKDC18nIfW13Wwhc8oVyLWMx84K+9Ocm5gWC815nA8Uy62jjb9Bu9kVI7e7Svqjr\nvsiJMXVttwHVuVyencmw4Qtbxjtbnr2rhVMk+806RWUAACGoOqUBmTrUwuUMQ9EnvUDbspnNSVey\nR48enX7GhW7NbUlfy652+QYAoA1OFEqB/zY3N9Xv96feo2sUbLa3t5UkiYbDoZIk0fb2dttNcsrG\nxobu3r2r3d1dXbp06Uz590m35q6b1eUbALqC31RMkKlDLRg71bzt7W1tbW151fWM7MpiKJySj30D\noOv4TUUaQR1qQzfH5vj6Re5qQR1XudytuW3sGwBdx28q0uh+CXjI165nZFcWQ7fmfOwbAF3HbyrS\nyNR5wMdudlguX7/Iya4shm7N+dg3ALqO31SkEdQ5ztdudlguX7/ImWtscXRrzpfdN5OCATwAA9AF\n/KYijaDOcfSXho2vX+RkV6ojc2/HAzAAXcNvKtKYfNxxLk/qjXYxqXr32AKXfr+v27dvd/7YJ0mi\njz/+2Pr+3bt3W2gRAADFMfl44HztZoflo1te95C5z+frOFMAAOpA9UvHUeENwASBSz4egAEAuoxM\nnRTxlx0AACAASURBVOPoL40uYJxYMQQu+XwdZwoAQB0YUwcExMfgiHFixbGvZsuOM11fX9f9+/e9\nuh4AAN1UdUzdzKAuiqIVSfcknZP0Q5J+wRjz01EU/YikvyLpt0j6vqQ/Yoz5Fcu//76kZ5L+laRf\nN8b8hOUzBHVADXy94afAxWIokFOMr9cDAKCblhrUvf4D3zDG/FoURT1JvyjpT0j6A5L+hTHmgyiK\n/qSki8aYn7L8238q6ceNMT+YsX6COjTKx2xWEb4GR1R4xTL4ej0AALpp6dUvjTG/9vo/f0jS1ySd\n6FVQN3j9/vck7Uo6E9RN2li2cUDdQp7LytciGowTwzL4ej0AAFDG3OqXURS9EUXR35f0y5L+tjHm\nM0m/2Rjzy68/8suSfnPOPzeS/o8oiv5uFEV/rJYWAxXMKgnvO1+DIyq8Yhl8vR4AACijSKbuX0v6\nHVEU/SZJO1EU/Z7M/2+iKMrrP/kfGGOeRlF0SdLfjKLoHxtj/k72Q++///7pfw+HQw2HwwU2ASgu\n5Kf3vlb/o8JrNaF2J67Kdj3EcayjoyMNh0P2FQCgVbu7u7UOMyk8pYEx5lejKNqW9OOSfjmKotgY\ncxhF0WVJRzn/5unr//3nURT9r5J+QtLMoA5YppCf3vscHDGRejkhdyeuKns9PHv2TE+fPtWjR49O\nP8O+AgC0JZvI+u53v1tpffOqX/6opJfGmF+Joui8pB1J35WUSDo2xvxMFEU/JemHs4VSoij6hqSv\nGWP+ZRRFv0HSx5K+a4z5OPM5CqWgMVTE8wPZp2IoBlIc+woA4LJlF0q5LOl7URS9oVfj737eGPO3\noih6JOmvRlH0X+j1lAavG/MtST9njNmQFEv661EUTf7OX8oGdEDTfM5mdQXZp+JC7k5ct7x99cUX\nXyhJEh4gAAC8NjOoM8b8Q0nfsbz/A0m/z/L+P5O08fq/P5f0O+ppJlAfuvq5bVYxG47btJC7E9ct\nb199/vnn+uyzz05f8wABAOCjudUvAaBJZJ+Ko3JocbZ9df78eX355ZdT74VSDRcA0C2FC6UAQBPI\nPhVHd+LibPvqyZMn2t/fP/PZhw8fUiETAOCVmYVSGmkAhVIApFDMBk3JK56SxrkHAGhC1UIpBHUA\nnLO9vT2VUVlfX9f9+/cpZoFa2R4g2KyururKlSucewCApSGog7MoS486kLnDMqUfIOzt7enk5GTm\n5+M41uXLl3XhwgW+1wAAtSGog5O4EXdDCIE184uhKUW6Y2YR5AEA6rDseeo6J4SbYBdQlr59ocz3\nRjXMxfAdVt7m5qYODg7mdsdMOzw81OHh4enrg4MDffrpp3QXBgA0iqAuJZSbYBdwI96+UAJrqmEW\nx3dYNdkKmfv7+zo+Pl5oHQcHB/rggw+mpkrY29ubyuZdvXqVoA8AUCuCupRQboJdwI14+0IJrG3Z\nkziOdXR0RNn5DL7DqtvY2DjdV0ULqWRl577LZvM++eQTvXz58vT1vKCPIBAAMA9BXUooN8EusN2I\nMylys0IJrLPZk2fPnunp06d69OjR6WfIRr3Cd1i98s69dIBWRjqgk+YHfdnXti6e8wJBAkMACBuF\nUlIoyFCvbFl6JkVuVqjFarhO87Fvli/9vWYL8s6fP38mU7cM2b/T6/WmAr95r20FXiRNjcdcNFAs\nEkjW/Td8WSft9n+dtNv/dbra7slDtqqFUmSMaXV51QQ3jMdj0+/3jaTTpd/vm/F43HbTsETj8diM\nRiMzGAzMaDQK6niPx2OTJIkZDAYmSZIgtm0wGExdo5NlMBi03bTW8R3WvOw1duvWrTPHwIcljmMT\nx/HUe71er9bXy/gbvqyTdvu/Ttrt/zpdbffkd1qSMRViKjJ1GWSXuiXUbFbIyEbNxndY++Zl87JZ\nMwBAtyVJop2dHf/nqXsVqAIAAABAFwUwT51DiTpvMBdVPYbDoe7du3fm/cFgoN3d3eYbhELmZULI\ntn6F7wr3ZbOr6+vrevDgwczXd+7cmephsOiYOgCAO15l6qqtw4mgDothLqr6hFIhsmvSZeeTJJmq\nhClRxn+C7wo/pM/not55552FA8HJa9uDkDiOJWlmN9Gqr5fxN3xZJ+32f5202/91utruSXX4nYpR\nnRPdL9tug28YU1QfxtT5Ly/bevHiRb399tudzk7xXYE8trGXkkoHikVeL+Nv+LJO2u3/Omm3/+t0\ntd2Tse9Vq18S1HmILoP1orCE3/ICl7SuBup8VwAA4IeqQR3dLz1El8F6len6BHfYJrrP6mp3zLzv\nimfPnilJEsbZAQAQCII6D9luYif9cYGumQQjk2zr3t6eTk5Oznzuiy++6FwgY/uuiONYT58+nRqH\nyDg7AAD8RvdLT9FlELDL6455/vx5ffnll6evu9IlM/tdcXR0dKawjMQ4OwAA2sSYOgCFdaG8va34\nTTagm1hdXdWVK1eC3Rc2FJYBAMA9jKkDUEhXyttnu2OurKzoyZMn2t/fP/PZ4+Pj0wAnxH1hkzfO\n7uTk5HRf7O3t6fLly7pw4QJBHgAAHiBTB3REl8vbF6mQOflc6PvCFtzPE8cxQR4AAEtEpg5AIS9e\nvLC+//z584Zb0rwiFTIl6eHDhxoOh0EHLkULy6QdHh5OTaRKJg8AALcQ1AEd0eWpMLKBzP7+vo6P\nj898ritdENPTeBTNYqZlg7yDgwN9+umnun///ul4zatXr069Dmn/AQDgGrpfAhmhFhOxdbvrSgXI\nrDJdEEPdV2X2hU22GE2v19PLly9PX2e7cGaDPoJAAECXUf0SqFHogQ9TYXwlvS+KdEGUwq2Wmd4X\nz54909OnT6cyccuQDfrmBYGbm5uSFOQDFwAACOqAGnW5mEiXlemC2O/3df369SCzS20EefPEcSxJ\nU+0g+wcACAVBHVCjvDm8BoOBdnd3m28QGlG2C2JXJjSfF+TlzQPYtDqyfwSCAIA2ENTBK66PVyNT\n1111ZafSXTRDDRCy3XjX19d1586dqaA4G1C5yJb9q2MsoDQ7UKz6mq6oABAegjp4w4fxaj60Ec1I\nBy551TLnKZIpCuW8sgV6Dx48yA2S52XVXFUkGyjNDhSrvq6jK6orwWhX1km7/V8n7fZ/na62e3Iv\nQFAHb/iSBaOYCLJswX4dXQ5DDvKyZgV9RYJAWyCDfHUEik0Eo11ZJ+32f5202/91utruSfLg3Xff\nJaiDHxivBp8V6XJYVcjFVxZle7giaWYXWV+zfwCAbkuSRDs7O/4HdRJBHQAAAICuqtb9sldnU8oi\nUdcNjFdDyGzn9zIyRV3qslm3edm/OsYCttVlCADgr1eZumrrcCJT13Yb0BzGqyFkixYLqQMPRpZr\nkbGARQLFOl7X0RXVlWC0K+uk3f6vk3b7v05X282YOgCVuT7FRGiWNd9behoFjmE3LRp8uhCMdmmd\ntNv/ddJu/9fparsnCQ6qXwIohe6w7VtG8RWOIQAA/iGogyQyLssW4v71ZYqJrqljEnSOIQAAfqka\n1DlRKAXV2DIuk//2PfBwQaj798WLF9b3nz9/3nBLkLaxsTF1XpUJ8h4+fKjhcBjMAwgAADAbmboA\nkHFZrlD3b6jbFbp0kLe/v6/j4+OZn6c7JgAA7quaqXujzsagHWRclivU/bu5ual+vz/1Xr/fPx3k\nCzdtbGzo7t272t3d1fe+970zxzDr4OBAH330UUOtAwAAbaD7ZQDOnTtnfX9lZaXhloQp1P07ydww\nxYS/ssdwb29PJycnZz7n+wMIAAAwG90vA0AVw+Vi/8IXeV1qmfIAAAC3USgFZFyWjP0LX2xuburg\n4GDqAUSv19Px8bHu3bsnKYwiPwAAYBqZOgAISJFCKhTDAQDALcxTBwCwGg6Hpxm6tMFgoN3d3eYb\nBAAArOh+CQCwyivy8+zZMyVJohcvXjDODgCAABDUASVsb29ra2sruJviULerq2xj7OI41tOnT/Xo\n0aPT9xhnBwCA3wjq0CofgwhbNcwQbopD3a4usxX5OTo6mgropK/msuM4AwDgJ8bUoTW+ThWQVzbe\n9+IToW4XpjHODgAA91QdU/dGnY0BFrG1tTUV0ElfZQxc9uLFC+v7vk/wHOp2YVreOLv9/X0Nh0Ml\nSaLt7e2GWwUAAKqg++UcPnYP9IWvQUTeTfHKykrDLalXqNuFacxlBwBAeAjqZmCM0XL5GkTYbor7\n/b5u3rzZYquqC3W7MC07zs42l93BwYFu3LihK1eu8DALAAAPMKZuBsYYLZevY+qk6QmeV1ZWdPPm\nTefbXESo24V8eWPs0ny5LgEA8BWTjy8RBQWWjyACaFfew6us1dVVMncAACwJk48vka/dA32ysbER\nxM1hqGMvQ90ufMXW7daGMXcAALiLoG4Gn8cYcTPenFDHXoa6XZhWZIxdFmPuALiEex6A7pdz+dg9\n0Oexaj4KdexlqNuF2WzfH/Pw/QKgLdzzwDd5DyHofrlkPnYPnDX/m2/b4gNfp2aYJ9TtwmxlM3fv\nvfceT8oBNI57HvhkVi+oqgjqAsTNeLNCHXuZt13Pnj1TkiTcvAcs/TCraObu8ePHevTo0enrvb09\nXb58WRcuXOA8AbA03PPAJ7MeQlRFUBegUIMMV/k89nIW23bFcaynT59O3bwzzi5sRTN32Ruow8ND\nHR4enr7OBnlXr17V/fv3eTgAoBLueeCTZT6EYExdgOhf3rzs2Mv19fUgbliz23V0dDQV0E0wzq47\nbN8vKysrC/8g9Xo9vXz58vR1HMdk9gAsjHse+GRWvYKdnR3mqcNZPhZ4kcKoYBXyDwxzN0IqHuxX\nYQvyJHn//QCgfr7e86B7Zt0jvvvuuwR1CEMowVDIVSPzto2JqbutTMXMRcVxLElTXTqzgR9dOgEA\nrst7CFG1+iVBHZwRSjAUcjbLdvOe7UbnYyCO6tI/Us+ePdPTp0+nArAmzOvSaQv6JLJ/AID2EdQh\nGKEEQ6Fns9I373lFM3wLxFG/eUFeNgBrgi3okxbL/s17bQsUF11HG3+DdrMvQmo3+8L/dndpX9Q1\nTx1BHZwRSqauS9msvED84sWLevvtt70PYFEfWzGhBw8etJrZKyJ77c57bQsUF11HG3+DdrMvQmo3\n+8L/dndpXzCmDsEJZUyd1J1sVl4gnubrMUSz5mX2bD+UAACEIJjqlxJBHQAAAICuqtb90onJx0nU\nIWR52ay1tTVdunTJ+wIN6QzL3t6eTk5OznyG7phYhkW7dBbpJgMAQNNeZeqqrcOJTF3bbeiCEOZ/\n85WtW6ntZjKEbop0x4RLZgV9kzLSkhYq6MI4ELf+hi/rpN3+r5N2+79OV9vNmDoUFtJYNV8VnazZ\n93F2RecrC6USKMIzLxCc9zobKJZZRxt/g3azL0JqN/vC/3Z3aV80Mk9dFEUrku5JOifphyT9gjHm\np6Mo+sOS3pf02yS9Y4z5pZx/f03Sn5X0NUl/wRjzM5bPENQtWShVJUMSctXIIt0x03jAAAAAuq5q\nUDdzTJ0x5nkURb/HGPNrURT1JP1iFEW/W9I/lPSHJP3sjIZ9TdKfk/T7JD2R9GkURX/DGPO4bGNR\nzosXL6zvP3/+vOGWLC7UbqPnzp2zvn9ycnIa7E2yXb5t78bGxmmbi3THPDg40HvvvRfkcQYAAGjC\n3EIpxphfe/2fP6RXGbcfGGP+sfQqopzhJyT938aY77/+7F+W9AclEdQ1LC+AWFlZabgli7F15fM1\n0Mna3NzUwcHBzG6KBwcHunHjhtfdFItspyQ9fvx4qjvq3t7e1ATQPm47AABAU+YGdVEUvSHplyT1\nJf15Y8w/KrjuNyX9v6nXX0j6XQu3EJXZbqz7/f5pv19XbW1tnQkGDg4O9NFHH3l/gz9p/7xuisfH\nx15n7rLbmTdnXzZrfHh4ODXImCAPAAAgX5FM3b+W9DuiKPpNknaiKBoaY3YLrJuBco7I3linB2W6\nzOduo0WU6aboY+YuvZ227OvKysrcYzovyLt69aru379P900AANBJheepM8b8ahRF25J+p6TdAv/k\niaQfS73+Mb3K1p3x/vvvn/73cDjUcDgs2iwUlL6x9oWv3UbLKNpNMbTM3axKoLNkg7xPPvlkqjyw\nLbMnaWrcHoEgAABoy+7urnZ3d2tb37zqlz8q6aUx5leiKDovaUfSd40xf+v1//+3Jf0JY8zfs/zb\nnqR/Iun3Svpnkh5K+qPZQilUv0Serk3FkK4amddNMSuEqQGKToNQRZF5Y+I4npn9swWB0uxAcdHX\ndazT1/MAAIAuW/aUBr9d0vckvfF6+XljzIdRFP0hSVuSflTSr0p6ZIz5/VEUfUvSzxljNl7/+9+v\nr6Y0+IvGmD9j+RsEdS3wpapkds4oH7qN1qFMoONzwJs+zrYJoNvgy6SliwanrgajvgbR7Ivw2s2+\n8L/d7Av/292lfTG5B68a1MkY0+ryqglo0ng8Nv1+3+jVuEcjyfT7fTMej9tuGlLG47FJksQMBgOz\nuro6dbzylrW1NTMajcxgMDCj0cjbY5re9rW1NRPHcaHtZzm79Hq9qddxHJ/Zn9nPLPp6Gev0td3s\nC//bzb7wv93sC//b3aV9MbkHl2RMhZhqZqauCWTqmufzZOS+ZBjrVjRzly064nP2Lm1eJi+brQIA\nAPBFkiTa2dlZXvfLJkRRZF4FqgAAAADQRdW6XxaufrlMJOqa5Wumztd2L0PZqQFCydzNYxuLKU1X\n3VxfX9eDBw8KZ/98GVMHAAD88ipTV20dTgR1aJavk5G/eBH2vHWLKDs1QCiTt8+TN4XHrO3OBoLp\noM/2el6gWOZ11XUWCU5dDUZ9DaLZF+G1m33hf7vZF/63u0v7YnIPvlMxqnOi+2XbbegiH6tKkqmb\nrei4u4sXL+rtt9/u1JjELpkXnLoYjPoaRLMv2v8bvqyTdvu/Ttrt/zpdbffkHnypUxo0gaDOHa4X\nIenavHVlLDrXHfsPAACgfQR1qIUvAVPeWCmXg9G2FM3ckekEAABoV9WgjjF1kPQqKMre/Ls4/io7\nVsoWuEz+26V2tyE77m5vb08nJydnPvf8effGJAIAAITkjbYbADf4WoRkVjCKV4Hd3bt3tbu7q3fe\necf6mZWVlYZbBQAAgDoR1C1oe3tbSZJoOBwqSRJtb2+33aRanDt3zvq+6zf8vgajbdjc3FS/3596\nL45jHR0dBXc+AwAAdAndLxcQclc/X6c5yAtG9/f3NRwOGWOXku2OOSl/n54GIZTzGQAAoEsolLKA\n0Evq+zjNgS3Qts3/4VrBFxeEfj4DAAD4gkIpDQq9q5+tCEmSJE5Xlcxmn2xl/F0s+OKC0M9nAACA\nriCoW4Cv487K8KmraToYHQ6Hunfv3pnPEKic1aXzGQAAIGQUSlmArdCED+POyvC1qiSBSnEUTgEA\nAAgDmboFZLv6+TLurAxfu+bZCr6kAxVXu5G2gcIpAAAAYSCoW1B23FmofM14EagsJn0+J0kytZ8k\nxiMCAAD4gO6XsPK5q2l6wu1Lly7p8PBw6v/3oRtpG3zNzgIAAHQdmTpY2bqarq+va2trSx9++KE3\n3RgJVIrzNTsLAADQdQR1yJXumudTNcw0ApXiGI8IAADgJ4I6FDKrGqbLN/kEKsUxHhEAAMBPBHUo\nxNdujAQqi6FwCgAAgH8olIJCfO7GSOGUcnwN5AEAALqGoA6F+FwNM41ApTifA3kAAIAuofslCgml\nGiaBSnG28Yj9fl/r6+tKkkQvXrzw5rgDANAF29vb2tra4je6gwjqUFgI1TApnFJcXiB/584d7447\nAACh8/XeDPWIjDHtNiCKTNttwOKSJNHHH39sff/u3bsttKi47e3tM4VT0uPs+v2+bt++zReghc/H\nHQCAkPEb7bcoimSMicr+e8bUoRSfx6ZROKU8n487AAAh4ze62wjqUEooY9P4AlxMKMcdAIDQ8Bvd\nbQR1KCWUaph8AS7GdtzTYxKTJNH29nZLrQMAoLtCuTdDORRKQSm2Iho3b970bhwahVMWw2TuAAC4\nKZR7M5RDoRR0HoVTymNQNgAAQHUUSoEztre3lSSJd93wKJxSHmMSAQAA2kf3S9QilLlRCFIWkzcm\ncX9/n+6rAAAADSGoQy22tramAjrpqwyXTzf0eUHKs2fPlCSJXrx4QaCSYhuT2Ov1dHx8rHv37kny\nM7gHAADwCUEdahFKhiuvcArFQOyyg7L39/d1fHw89Rkfg3sAAACfENShFqFMDWCrHHV0dDQV0EkE\nKmkbGxun+2E4HJ5m6NJ8C+4BAAB8QqGUinwtDlK3kOZGSRdOuXv3ri5cuGD9HIHKWfPG2HX5GgEA\nAFgWMnUVhFIcpA4hz40SShayCYyxAwAAaB7z1FXAHF3dYAve4zjW5cuXdeHCBQqnZKTn/bONsZO4\nRgAAANKqzlNHpq6CUIqDYLZsFnIyQTmFU+wYYwcAANAsxtRVQLe82UIab8gE5eUwxg4AAGD5yNRV\nYBs/5GtxkLqFPN6QDG1xjLEDAABYPsbUVZQePxRScZCqQh5vmLdtq6urunLlCmPsMhhjBwAAMFvV\nMXUEdViKvLFUg8FAu7u7zTeoRrYsZK/X08uXL09f9/t93b59m8AuI++8eOutt/Tmm2/qxYsXBMUA\nAKBzKJQCJ4U83jBbOMWWfWJycru88+Lzzz/XZ599dvp6b2+P6qIAAAAFEdRhKUIfb0iFx3Js58X5\n8+f15ZdfTn3u8PBwqhgN4+4AAADy0f0SS9OV8YaMsVtM9rx48uSJ9vf35/679P68evWq7t+/T3dN\nAAAQBMbUAS1jjF01eUHxLNn9y2TwAADAZwR1gAOo8FieLSiuyhbkSdLW1tZpdi+b7Vv0dR3rJPgE\nAAASQR3gnLwxdhcvXtTbb7/NzbxFOih+9uyZnj59emaC9yriOJakqXVms32Lvq5jnU0En8sIRptY\np6/tZl/43272hf/tZl/43+4u7YvJPWHVoE7GmFaXV00AwjEajYykmUu/3zfj8bjtpjprPB6bJEnM\nYDAwq6urc/dnKEscxyaO46n3er1era+X8TdoN/sipHazL/xvN/vC/3Z3aV9M7gklGVMhpiJTh0Zt\nb29PPb0IMWNVtDsh3TGLKTJmEQAAwFdJkmhnZ8f/7pevAlUAAAAA6KJq3S/fqLMpZRnD0oVlNEok\nRWeWJLnWetva2O7z578x9brf/7bG4+3W2+vjMh5vK0muaTAYam3tO4rjy1P7No4vn3mv1/t6pdd1\nrJOFhYWFhYWFJUmuqSomH0djXrx4YX0/9Em6i064fXBwoBs3bjC3XQnpyeAl+xyJkqbeW19f14MH\nD0q/rrpOW0GYZRR0WUaBlybW6Wu72Rf+t5t94X+72Rf+t7tL+6Lf7+vmzZva2dlRFU50v2y7DWhG\n3nxkXRhbVmbCbeZeC18TwWfdwWhT6/S13ewL/9vNvvC/3ewL/9vdpX1x8+bNWqpfEtShMbaCF12d\nlLvMhNtd3VcAAAChI6iDV2yZiS4GKWUn3F5dXaV7JgAAQGAI6gBPpQPc/f19HR8fL/Tv6Z4JAAAQ\nBoI6x3RhHjbUr2zmLq3f7+v69eu6f/8+5x8AAIBHCOocwpgxVJHO3NkqIxaRrapJNg8AAMB9BHUO\n6XJ1R9SvavdMG4I8AAAA91QN6pinrkZdnYcNy5Gee62O7pnSq3lT0tm/yfoI7AAAAPxFUFejc+fO\nWd9fWVlpuCUIzSTomtU90zah+TxMeA4AAOA/ul/WiDF1i6OwTHnZ6SHW19d1586dStk8umcCAAA0\njzF1jmEetuIIgutXR7GVNII8AACA5SOog7coLLN8dQd5BN0AAAD1I6iDt4bDoe7du3fm/cFgoN3d\n3eYb1AF1VNRcW1vTpUuX6DILAABQE6pfwlsUlmleHRU1Hz9+rEePHp2+3tvbo4smAMApjNlH1xDU\noTWbm5s6ODg4M6bu5s2bLbaqO4pU1LTJTtHBNAkAAJfYHlry24TQ0f0SraKwjFvmjcFbWVkpNO/i\n6uoq0yQAAFrBmH34iDF1AJYmG3QfHR1Ndb0sguIqAIAmMWYfPmJMHYClSY/Bk8qNw2OCcwBAkxiz\njy6aGdRFUbQi6Z6kc5J+SNIvGGN+OoqiPyzpfUm/TdI7xphfyvn335f0TNK/kvTrxpifqK/pAJpW\ndhze8fHx6VNTxjUAAJaJMfvoorndL6Mo+oYx5teiKOpJ+kVJf0LSv5D0ryX9rKT/ekZQ908l/bgx\n5gcz1k/3S8BjZaZJYMwdAGCZGLMP3zQ2pi6Kom/oVdbuhjHmH71+729rflD3O40xuXd5BHVAOMp0\nz+z3+7p+/bru379P6WkAANBJSx9TF0XRG5J+6f9v745j4zzv+4B/H+saUurmzlGMnNIFaHfogNVu\nEAFDJqEYyKCDz4GCtcGGdgUMeMMQDOhEGgOKLIGRWf6rWbIONWUsGOJsCOIuQ4EtXafDZKZBZWSA\ntSSLZleOu7XEArSdHHXKHHeYxIbJsz9IySR1pI48knfv8fMBCPN9dXzuhyd5+fJ7z/M+T5JOkk/f\nDnQDqkl+u5Ty/ST/stb6md2VCTTB5umZg4zcLS0t5ZOf/GRu3ry54dz69gAA2NpORup+JMkLST5a\na720du5eI3Unaq3XSikPJvlSkrla61c2vWbiR+psgMlhtdsNzhNTNAGAw+PAVr+stX63lNJL8leT\nXBrwZ66t/fdPSilfTPK+JF/Z/Lpz587d+X52djazs7ODljX2bIC5MwLwZNnNyN1t6xdXeeWVV3Li\nxIncf//9mZqayunTp03XBAAa69KlS3u6xca2I3WllHckWam1vlFKOZrVkbqna61fXvv330nyy7XW\n/9rnZ48lOVJr/dNSyg8nWVz72cVNr5vokTobYA6uXwC2x9lk6fe/8dGjRzdMvRxEq9XKysrKneN2\nu33P0JdkwwcGgiEAMC72daGUUspPJflckvvWvj5fa/1UKeVDSRaSvCPJd5NcqbV+oJTyriSfnD1r\niAAAE8pJREFUqbWeKaX8pST/fq2pVpJfr7X+Sp/3mOhQZwPMwQnAh8PmFclOnTqV559/fldTNLfS\nL/Ql2bD1wk6D4W6C4k6P96NNdTe/TXU3v011N79NdTe/zXGt+/aHysOGutRaR/q1WsLkeuSRR2pW\nF4zZ8NXtdkdd2tiZmZnp21czMzOjLo19duHChdrtduvMzEw9fvx43/8fHPRXq9Xa9rjdbtd2u72j\nnxlFm+pufpvqbn6b6m5+m+pufpvjWnen06kXLlyoSWodIlMNvFDKfpn0kTpTCgdnpI5kuMVVAACa\nptvt5oUXXjiYfer2SymlrgZVAACAw+iAVr/cTxM8UMcObX7eam5uzogmfZ/Du3z5cm7dupU333wz\n165du+fzcsn2z9QBAIzC6kjdcG2MxUjdqGsAmm270Hf7w4EkQwXD3QTFnR7vR5vqbn6b6m5+m+pu\nfpvqbn6b41r37ceyPvjBDzZ/+uWoawC4VzDcaVDczfF+tKnu5rep7ua3qe7mt6nu5rc5rnXfnpW2\nr1saHAShDgAAOMyGDXX37WUxAAAAHCyhDgAAoMGEOgAAgAYT6gAAABpMqBuBXq+Xbreb2dnZdLvd\n9Hq9UZc0tvQVAABsbyw2Hz9Mer1ennjiiSwtLd05d/t7m2xvpK8AAHav1+tlYWEhy8vLmZqayvz8\nvL+hJpQtDQ5Yt9vN4uJi3/MXL14cQUXjS18BAOxOvw/Hb290LdiNH1saNMzy8nLf87du3TrgSsaf\nvgIA2J2FhYUNgS5ZnfF0/vz5EVXEfhLqDtjU1FTf89PT0wdcyfjTVwAAu+PD8cNFqDtg8/Pz6XQ6\nG851Op3Mzc2NqKLxpa8AAHbHh+OHi4VSDtjtOcznz5/PrVu3Mj09nbm5OXOb+9BXAAC7Mz8/n6Wl\npbueqfPh+GSyUAoAAEygXq/nw/GGGHahFKEOAABghKx+CQAAcIgJdQAAAA0m1I2BXq+Xbreb2dnZ\ndLvd9Hq9UZcEAAA0hNUvR6zX6+WJJ57YsDLR7e89yAoAANyLkboRW1hY2BDoktVQd/78+RFVNN6M\nagIAwEZG6kZseXm57/lbt24dcCXjz6gmAADczUjdiE1NTfU9Pz09fcCVjD+jmgAAcDehbsTm5+fT\n6XQ2nOt0OpmbmxtRRePLqCYAANzN9MsRuz1t8Pz587l161amp6czNzdnOmEfRjUBAOBupdY62gJK\nqaOugWbo90xdp9PJM888IwQDANBYpZTUWstuf95IHY1hVBMAAO5mpA4AAGCEhh2ps1AKAABAgwl1\nY8gG2wAAwKA8UzdmbLC9c71eLwsLC1leXs7U1FTm5+f1FQAAh4ZQN2a222BbULmbEAwAwGFn+uWY\nscH2zmwXggEA4DAQ6saMDbZ3RggGAOCwE+rGzPz8fDqdzoZznU4nc3NzI6povAnBAAAcdp6pGzM2\n2N6Z+fn5LC0tbZiCKQQDAHCY2Hycxuv1ekIwAACNNezm40IdAADACA0b6jxT1wA2I98Z/QUAwGHi\nmboxZx+2ndFfAMBmvV4vCwsLWV5eztTUVObn5/1dwEQx/XLMdbvdLC4u9j1/8eLFEVQ03vQXALBe\nvw98O51OnnnmmUMZ7ATc8TTs9EsjdWPOPmw7o78AgPUWFhY2BLpkdRbP+fPnD12YMaNpcnmmbszZ\nh21n9BcAsJ4PfN+yXcCl2YS6MWcz8p3RXwDAej7wfYuAO7lMvxxzNiPfGf0FAKw3Pz+fpaWlu56p\nO4wf+Aq4k8tCKQ3lIdfB6SsAONx6vZ4PfGPRmHFmoZRDyEOug9NXAMCZM2fc92NG0yQzUtdAlu0f\nnL4CAGDcDTtSZ6GUBvKQ6+D0FQAAk06oayAPuQ5OXwEAMOmEugaybP/g+vVVu93O9evXMzs7m263\nm16vN6LqAABgeBZKaSAPuQ5uc1+9+eabuXbtWq5cuXLnNRZOAQCgySyUMiEs2z8YC6cAADBubGmA\nZft3wMIpAABvMTAwGYS6CbCwsLAh0CWroe78+fMuyk22Wjjl6tWrmZ2d9csMADg0DAxMDqFuAhh9\nGtz8/HyWlpY2/PJqtVq5ceNGXnzxxSR+mQEAh4OBgckh1E0Ay/YPbvPCKVevXs2NGzc2vGZpaSmP\nP/54Hn74YSN3AMDEMjAwOYS6CdBv9Gn9sv2CyUZnzpy50xezs7N3RujWM3IHAEw6AwOTQ6ibAJbt\n372tfpmtt7S0lI9//OMeIgYAJkq/gQF7HzeTLQ0mkGX7B9fvAeF+pqenN0xFaLfbOXHiRO6//34h\nDwBorF6vt2Hv41OnTuWll17yQfYBG3ZLA6FuAm01pfCBBx7Ie97zHhfoJut/mfV7xm4QQh4A0HT9\nPuzudDp55pln/F2zz4Q67rLVSN16LtD++v0y2zxKN4jNIe/06dN3feqVxJROAGBsbPU35PHjxy0g\nt8+EOu4y6JRCF2h/m6chXL9+fcPzibvRarWysrJy57jdbidJXn/99Q3ntguC9zruFxR32sYo3kPd\n+mKS6tYXza9bXzS/bn2x+/d4+eWX88Ybb2z7N82xY8dy5MiR3HfffWm1Wjl79myS5Nlnn83Kykpa\nrVYeeuihvPrqq7s+3o82D+I9dtPm2bNnc+7cuaFDXWqtI/1aLYG9duHChdrtduvMzEx94IEHapJt\nvzqdTr1w4cKoyx5LFy5cqJ1O5559uNdfrVZrR8ftdru22+2h2hjFe6hbX0xS3fqi+XXri+bXrS+G\ne4/dfJVShm7jINocx7pbrVZ96qmnapJah8hURuoOgUGmYybJyZMn8+CDD274NMfo3ar1o3e3Vxdd\nP8oGAAC7cfz48dy4caP50y9XgyoAAMBhNNz0y7HYp85A3f6zwuP+utdI3iDP1AEAcPisjtQN18ZY\njNSNuobD5iBXeDysoa/fni+XL1++c3x7U8+dBMF7HfcLijttYxTvoW59MUl164vm160vml+3vtjb\n9zhy5Ei+//3vZztri3xs+5qd2o82D+I9dtpmq9XKk08+maeffrr50y9HXcNhdFArPBrZG9y9guC9\njjcHxd20MYr3ULe+mKS69UXz69YXza97kDbf/va3Z3Fx8c4KhI888ki+853vjH3d4/Aem/tuXFaR\nHIf3GOXql0IdSQbfBmEYnU4njz322D2X2hX8AID9YoNtxpFQx545iBUejx49mps3b945NroHAByk\nrVYF73a7uXjx4ggqguFD3VgslMJ4OHPmzIYAtR8hb32gS7Ih0CWrc7zXv8crr7wi5AEAe2Z5ebnv\n+Vu3dra2AIyTbUNdKWU6yYtJppK8Lcl/qLV+rJTyqSQfTPJnSZaS/L1a63f7/PyjSX4tyZEkz9Va\n/+ke188+2mnI2zzqtheEPABgL01NTfU9Pz09fcCVwN655/TLUsqxWuv/K6W0kvznJL+c5GiSL9da\nf1BK+USS1Fo/uunnjiT570n+RpI/TvK1JL9Ya31t0+tMv2yo7Rb26Bf6Nk+93Av9ntMT9ACArXim\njnF0YM/UlVKOZXXU7vFa6zfXnf9Qkr9Va31s0+tPJ3mq1vro2vFHk6TW+olNrxPqJlS/0Pf8889v\n+CW6F6N7m8OiX8wAwHY2/40yNzfn7wZGat9DXSnlviTfSNJJ8ula60c2/ft/TPKFWuu/2XT+byfp\n1lo/vHb8WJK/Vmud2/Q6oe4Q2eno3m4dP348Dz/8sJE7AADG3r4vlFJr/UGS95ZSfiTJC6WU2Vrr\npbU3fzLJn20OdLd/dNAizp07d+f72dnZzM7ODvqjNMzm5/Q226vFWW7cuJEXX3wxiefwAAAYL5cu\nXcqlS5f2rL0dbWlQSvl4kpu11n9WSvm7ST6c5GdqrXctF1RKOZXk3Lrplx9L8oPNi6UYqWM79wp5\nu3lOz/RMAADGyb5OvyylvCPJSq31jVLK0SQvJHk6yQ8l+dUkM7XW/73Fz7ayulDKzyT5X0m+Ggul\nMKRBntMbhOmZAACMi/0OdT+V5HNJ7lv7+nyt9VOllN/P6hYH31l76Uu11l8qpbwryWdqrWfWfv4D\neWtLg8/WWn+lz3sIdQxlfdC7evVqbty4saOfN3IHAMAoHdjql/tFqGMv9VumeBBG7gAAGJV9XygF\nmuR2GNvpYivrF1a5HQgFOwAAmsBIHRNvN9MzjdwBwOHR6/WysLCQ5eVl935GwvRL2IHdTM/0zB0A\nTK5+fxu493PQhDrYod2M3J08eTIPPvigT/AAYMJ0u90sLi72PX/x4sURVMRh5Jk62KH1G6APOnL3\n2muv5cqVK3eOPXcHAJNheXm57/lbt+7ahhnGlpE6Dr3dbonguTsAaL6tRurc5zlIpl/CHuo3cjc9\nPX3PT+va7XZOnDiR+++/3y9/AGiQfvf+VquVlZWVO8eesWO/CXWwx9aP3E1PT+f69esbpl4OYnPI\nO336dF566SXP5AHAGBpk1o6RO/aTUAf7bLcbmq/X7xO/xx57bEPQ2xz8dno8Pz+fJBuWZB62zYN4\nD3Xri0mqW180v2590fy6h23z5ZdfzhtvvLHtfb3dbuf9739/FhcXs7KyklarlYceeiivvvrqro/P\nnj2bJHn22Wf3rM2DeI+mtDmudZ89ezbnzp0bOtSl1jrSr9USYLxduHChdrvdOjMzU48fP16TDP11\n9OjRDcetVmuo43a7Xdvt9p62eRDvoW59MUl164vm160vml/3XrQ5qq9SSiPfoyltjmPdrVarPvXU\nUzVJrUNkKiN1sEN7MXIHAADJ6tTeGzduNH/65WpQBQAAOIwmYJ86A3U03foHrN98881cu3Ytr7/+\n+p1/b216pu7o0aO5efPmKEoFAGCMrI7UDdfGfXtTChxuZ86cycWLF3Pp0qV84xvfyHPPPZdut5uZ\nmZl0u908+eSTG44/8pGPpNPpbGij1WoNddxut9Nut/e0zYN4D3Xvb5vqbn6b6m5+m+pufptHjhzZ\ncHzs2LG7zu2FUna/TsYo36MpbR7Ee+y0zda6BVaGet9xmH456hpgFDZvnXDq1Klcvnx518dzc3NJ\nsqdtHsR7qFtfTFLd+qL5deuL5td9UH3xta99beQrJ47DezSlzXGte69WvxTqAAAARmjYUGf6JQAA\nQIMJdQAAAA0m1AEAADSYUAcAANBgQh0AAECDCXUAAAANJtQBAAA0mFAHAADQYEIdAABAgwl1AAAA\nDSbUAQAANJhQBwAA0GBCHQAAQIMJdQAAAA0m1AEAADSYUAcAANBgQh0AAECDCXUAAAANJtQBAAA0\nmFAHAADQYEIdAABAgwl1AAAADSbUAQAANJhQBwAA0GBCHQAAQIMJdQAAAA0m1AEAADSYUAcAANBg\nQh0AAECDCXUAAAANJtQBAAA0mFAHAADQYEIdAABAgwl1AAAADSbUAQAANJhQBwAA0GBCHQAAQIMJ\ndQAAAA0m1AEAADSYUAcAANBgQh0AAECDCXUAAAANJtQBAAA0mFAHAADQYEIdAABAgwl10HCXLl0a\ndQnQeK4jGI5rCEZLqIOGcyOF4bmOYDiuIRgtoQ4AAKDBhDoAAIAGK7XW0RZQymgLAAAAGLFaa9nt\nz4481AEAALB7pl8CAAA0mFAHAADQYCMNdaWUR0spv1dK+f1Syj8eZS3QFKWUb5VSXimlXCmlfHXt\n3NtLKV8qpfyPUspiKeUvjLpOGBellH9VSvl2KeV3153b8poppXxs7b70e6WUR0ZTNYyPLa6hc6WU\nP1q7F10ppXxg3b+5hmCdUsq7Sym/U0p5tZRytZQyv3Z+z+5FIwt1pZQjSZ5N8miSn0zyi6WUvzKq\neqBBapLZWuvJWuv71s59NMmXaq1/OcmX146BVf86q/ea9fpeM6WUn0zyC1m9Lz2a5F+UUsxq4bDr\ndw3VJP987V50stb6nxLXEGzhe0n+Ua31oSSnkvzDtdyzZ/eiUV5k70vyB7XWb9Vav5fk3yb52RHW\nA02yeXWkv5nkc2vffy7Jzx1sOTC+aq1fSfJ/Np3e6pr52SRfqLV+r9b6rSR/kNX7FRxaW1xDyd33\nosQ1BHeptb5ea/1va9//3ySvJfnR7OG9aJSh7keT/OG64z9aOwdsryb57VLK10spH147985a67fX\nvv92kneOpjRojK2umXdl9X50m3sTbG2ulPJyKeWz66aNuYZgG6WUH0tyMsl/yR7ei0YZ6uylALvz\n07XWk0k+kNXh+7++/h/r6j4lri8Y0ADXjOsJ7vbpJD+e5L1JriX51W1e6xqCJKWUP5fk3yV5otb6\np+v/bdh70ShD3R8nefe643dnYyIF+qi1Xlv7758k+WJWh+O/XUppJ0kp5USS66OrEBphq2tm873p\nL66dA9aptV6va5I8l7emhrmGoI9Syg9lNdB9vtb6m2un9+xeNMpQ9/UkP1FK+bFSytuy+jDgb42w\nHhh7pZRjpZQ/v/b9Dyd5JMnvZvXaeXztZY8n+c3+LQBrtrpmfivJ3ymlvK2U8uNJfiLJV0dQH4y1\ntT9Ab/tQVu9FiWsI7lJKKUk+m+SbtdZfW/dPe3Yvau1tyYOrta6UUs4meSHJkSSfrbW+Nqp6oCHe\nmeSLq78b0kry67XWxVLK15P8Rinl7yf5VpKfH12JMF5KKV9IMpPkHaWUP0zyT5J8In2umVrrN0sp\nv5Hkm0lWkvzS2kgEHFp9rqGnksyWUt6b1Slh/zPJP0hcQ7CFn07yWJJXSilX1s59LHt4LyquMwAA\ngOaybwgAAECDCXUAAAANJtQBAAA0mFAHAADQYEIdAABAgwl1AAAADSbUAQAANJhQBwAA0GD/HzdM\nIvHqcSRGAAAAAElFTkSuQmCC\n", + "text/plain": [ + "<matplotlib.figure.Figure at 0x11097ee90>" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(figsize=(15,10))\n", + "\n", + "t=np.diag(alphas[:],0)+np.diag(betas[1:],1)+np.diag(betas[1:],-1)\n", + "\n", + "for i in xrange(2,M):\n", + " plt.plot(i*np.ones(i),np.linalg.eigvalsh(t[:i,:i]),'ok')\n", + "\n", + "\n", + "for i in xrange(20):\n", + " plt.axhline(evs[-(1+i)])\n", + "plt.ylim(1.01*evs[-1],1.01*evs[-20]);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "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.8" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/exercises/ex03_skeleton/trising_skeleton.cpp b/exercises/ex03_skeleton/trising_skeleton.cpp new file mode 100644 index 0000000000000000000000000000000000000000..700f6cae1f477568a78306804910343d3e06da08 --- /dev/null +++ b/exercises/ex03_skeleton/trising_skeleton.cpp @@ -0,0 +1,56 @@ +// ---- CQP FS16, Ex 3.1 Skeleton code ---- +// ---------------------------------------- +#include <valarray> +#include <vector> +#include <iostream> +#include <cmath> +#include <complex> + +typedef std::valarray< std::complex<double> > state_vector_t; +typedef std::vector<double> param_vector_t; +typedef unsigned state_t; + +// Implement the time evolution of state x by n*dt here. Use the split +// operator scheme we discussed in the exercise class. +void evolve(unsigned l, const param_vector_t &j, const param_vector_t &hx, double dt, unsigned n, state_vector_t& x) +{ +// ........ +} + +void printMagnetization(unsigned l, const state_vector_t& v, std::ostream& outputStream) +{ + std::vector<double> res(l,0.0); + for(int r=0; r<l; r++){ + double m(0); + // Measure magnetization at site r here. To be completed. + outputStream << m << "\t\t"; + } + outputStream << std::endl; +} + +int main(){ + // ======== input + unsigned l = 15; // length of chain + param_vector_t j(l-1,0.0); // Ising coupling + param_vector_t hx(l,0.4); // transverse magnetic field + double tmax = 20.; // maximum time reached in time evolutin + double tmeas = 0.5; // time after which a measurment is performed + double dt = 0.1; // time step used for the evolution + unsigned nmeas = tmax/tmeas; // number of measurements + unsigned nstep = tmeas/dt; // number of steps between measurements + // ======== input + + // Starting configuration + state_t initialInd = 1<<7; + state_vector_t x(0., 1 << l); + x[initialInd] = 1; + printMagnetization(l,x,std::cout); + + // Do time evolution + for(int n=1; n<=nmeas; ++n){ + evolve(l, j, hx, dt, nstep, x); + printMagnetization(l,x,std::cout); + } + + return 0; +} diff --git a/exercises/ex03_skeleton/trising_skeleton.py b/exercises/ex03_skeleton/trising_skeleton.py new file mode 100644 index 0000000000000000000000000000000000000000..0abce03467185ef06ae19b651e477df280650d4b --- /dev/null +++ b/exercises/ex03_skeleton/trising_skeleton.py @@ -0,0 +1,76 @@ +# ---- CQP FS16, Ex3.1 skeleton ---- +# ---------------------------------- +import numpy as np +import matplotlib.pyplot as plt + +def evolve(l, j, hx, v, dt, n): + """ Evolve the state v in time by n * dt + + Parameters: + l length of chain + j Ising coupling for each neighbors, 1-d numpy array with length at least l-1 + hx transverse field for each site, 1-d numpy array with length at least l + v the state at time t, complex 1-d numpy array with length at least 1 << l + dt timestep + n number of consecutive timesteps + + Return: + v the state at time t+n*dt + """ + # Implement the time evolution of the state v here. + # Use the split-operator scheme we discussed in the + # exercise class. + return v + + +def magnetization(v, l): + """ Measure the magnetization per site for the state v + + Parameters: + v the state to be measured, complex 1-d numpy array with length at least 1 << l + l length of chain + + Return: + m the magnetization at each site, 1-d numpy array with length l + """ + # Implement the measurement of the magnetization here. + return m + + + +l = 9 # length of chain +dim = 1 << l +j = 1.00*np.ones((l)) # Ising coupling +hx = 0.40*np.ones((l-1)) # transverse field strength Gamma + +# Starting configuration +ind = int('000010000',base=2) # specify spin configuration: 0 = down, 1 = up +v = np.zeros((dim),dtype=complex) +v[ind] = 1. # start time evolution with a single basis state + +# Parameters for time evolution +tend = 10. # duration of the time evolution +tmeas = 0.5 # time after which a measurment is performed +dt = 0.1 # time step used for the evolution +nmeas = int((tend)/tmeas) # number of measurements +nstep = int((tmeas)/dt) # number of steps between measurements + +# measure magnetization +m=np.zeros((nmeas,l)); +m[0,:] = magnetization(v, l) + +# Do time evolution +for t in range(1,nmeas): + print(t*tmeas) + v=evolve(l,j,hx,v,dt,nstep) + m[t,:] = magnetization(v, l) + +# plot magnetization +plt.pcolor(np.array(range(l+1)),np.linspace(0.,tend,nmeas),m, + vmin=-1.0, vmax=1.0) +plt.xlabel('site') +plt.ylabel('time') +plt.axis([0,l,0,tend]) +plt.title("Magnetisation, J = "+str(j[0])+" , hx = "+str(hx[0])) +plt.colorbar() +plt.show()