{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Lecture 4: Hard and soft constraints\n", "\n", "Goal: to illustrate examples of constraints in optimization problems " ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "scrolled": true }, "outputs": [], "source": [ "## MP 574 Lecture 4: Optimization Constraints\n", "##\n", "%matplotlib inline\n", "import numpy as np\n", "import numpy.linalg as la\n", "import matplotlib.pyplot as plt\n", "from IPython.display import display, Image\n", "import matplotlib.image as mpimg\n", "from os.path import dirname, join as pjoin\n", "from scipy import signal\n", "import scipy.io as sio\n", "import scipy.optimize as opt\n", "import numpy.random as rnd\n", "\n", "font = {'weight' : 'normal',\n", " 'size' : 20}" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhcd33v8fd3tFq2tVmyvEiWZFuJt2y27KyEkITg0CYpe1IoSymhD4Sl5dLSp/dJW3rv7UJvubQNgZQGSEqBACUN1CFAFrJ5k+zEidfIsiTLm3bZ1i7N7/4xM46iSNZInpkzc+bzeh491swczfmOz8xHP33POb9jzjlERCT1BbwuQEREYkOBLiLiEwp0ERGfUKCLiPiEAl1ExCcyvVpxSUmJq6qq8mr1IiIpqb6+vsM5VzrZY54FelVVFXV1dV6tXkQkJZlZ81SPqeUiIuITCnQREZ9QoIuI+IQCXUTEJxToIiI+MW2gm9mDZtZmZq9O8biZ2T+ZWYOZ7TGz9bEvU0REphPNCP07wObzPH4rUBP+uhu4/8LLEhGRmZo20J1zzwJd51nkDuAhF7INKDSzxbEqUCSRDp48w4PPH6Gx/azXpYjMWCxOLFoKHB13uzV834mJC5rZ3YRG8SxbtiwGqxa5cMGg4/FXT/LdrU3sOBIau3z553DdyhJ+7+pKbllThpl5W6RIFGKxU3Syd/qkV81wzj3gnKt1ztWWlk565qpIwv3dEwf49H/s4kTvAH926yp+/cfX88V3XExj+1k++XA9//JUg9clikQlFiP0VqBi3O1y4HgMnlck7h7dfYxv/qaR371yGf/rjnUEAqHxycqF8/nk9cv5Hz96mf/7q0NcvGg+t6xd5HG1IucXixH6Y8CHw0e7XAX0Oufe1G4RSTZ7Wnv405/sYVN1MX9529pzYR6RmRHgb99zKZeWF/BHP3yJQ6fOeFSpSHSiOWzx+8BW4GIzazWzj5vZH5rZH4YX2QI0Ag3AvwKfilu1IjHSfmaIux+qp2ReDvd/cD3ZmZN/FHKzMnjg92rJy8nkEw/V0ds/kuBKRaI3bcvFOXfXNI874NMxq0gkAb725CE6+4Z49NPXsmBeznmXXVSQyzc+tJ73fmMr9//mMF+6dVWCqhSZGZ0pKmnneM8AP9x5lPfXVrB2SUFUP7OhspjbL1vCQ1ub6Dw7FN8CRWZJgS5p576nQ0etfOptK2f0c5+5sYaBkTEeeK4xHmWJXDAFuqSVYz0DPFIXGp0vLZwzo59duXBeaJT+YrNG6ZKUFOiSViKj80/PcHQe8ZkbaxgaHeOBZzVKl+SjQJe00drdz4/qjvKBjRUsmeHoPOLcKH1rMx0apUuSUaBL2vj3bS0EHXzqhtmNziPuCffSf7jz6PQLiySQAl3SwljQ8dPdrdxwUemsR+cRKxfOY1N1MT+pbyV01K5IclCgS1p4vqGDU6eHeM+G8pg833vXl9PY0ceulp6YPJ9ILCjQJS38pL6VgjlZ3LR6YUye752XLmZOVgY/2dUak+cTiQUFuvje6cERnth7ktsvW0JOZkZMnnNeTiab1y3iZy8fZ3BkLCbPKXKhFOjie/+95wRDo8GYtVsi3rO+nDODo/xq36mYPq/IbCnQxfd+Ut/KitK5XFYe3Wn+0bp6xQIWF+Sq7SJJQ4EuvtbU0Uddczfv3VAR86sOZQSMd69fyrOH2mk7PRjT5xaZDQW6+NpjLx/HDN51xdK4PP+715cTdPCzPboEgHhPgS6+9qt9p7iiopBFBblxef4VpfO4qGwev1YfXZKAAl1860TvAK8c6+XmNWVxXc/Nq8vY0dSli1+I5xTo4ltP7m8D4O2r4xzoa8oYCzqeOdQW1/WITEeBLr716/2nqFyQx8qF8+K6nsvLCymZl63DF8VzCnTxpb6hUV5s6OTm1WUxP7plokDAuGlVGb852M7waDCu6xI5HwW6+NJzr7UzPBbk5ji3WyJuXlPGmaFRdhzpSsj6RCajQBdf+tW+NgrmZLGxqigh67tuZQm5WQF+vV9tF/GOAl18ZyzoeOrAKW5ctZDMjMS8xedkZ3DdylJ+te+UptQVzyjQxXd2tXTT3T+SsHZLxNvXLORYzwAHTp5J6HpFIhTo4jtPHWgjM2Bcf1FJQtd746qyc+sX8YICXXzn+dc6WF9ZxPzcrISut3R+DmsW5/P8ax0JXa9IhAJdfKW7b5hXj/dy3crEjs4jrqspob65m4FhzZEuiadAF1/Z2tiJc3DtygWerP/alSUMjwXZ2aTDFyXxFOjiK883dDAvJ5NLyws9Wf/GqiKyMwK80KC2iySeAl185YWGDq5aXkxWgg5XnCgvO5P1lYU8r0AXDyjQxTeOdvXT3NnPtR71zyOuW1nC3uOn6eob9rQOST8KdPGNSJvDqx2iEZFfKC8e1ihdEiuqQDezzWZ20MwazOxLkzy+zMyeNrPdZrbHzN4Z+1JFzu/5hg4Wzs+J++yK07lkaQHzczPVR5eEmzbQzSwDuA+4FVgD3GVmayYs9j+BR5xzVwB3Al+PdaEi5xMMOl483Ml1K0viPrvidDIzAly9fIH66JJw0YzQNwENzrlG59ww8APgjgnLOCA//H0BcDx2JYpMb//JUM/a6/55xHU1JRztGqCls9/rUiSNRBPoS4Gj4263hu8b7y+BD5lZK7AF+MxkT2Rmd5tZnZnVtbe3z6JckclF2hvJEuiROp5r0PtcEieaQJ/s79eJ08ndBXzHOVcOvBN42Mze9NzOuQecc7XOudrS0tKZVysyhW2NXSwvmRu3i0HP1PKSuZTl57C9UScYSeJEE+itQMW42+W8uaXyceARAOfcViAXSI6hkvjeWNCxs6mLK5cXe13KOWbGpuoF7DjSpel0JWGiCfSdQI2ZVZtZNqGdno9NWKYFuAnAzFYTCnT9rSkJceDkac4MjrKpOnkCHWBTdTEnTw9ytGvA61IkTUwb6M65UeAe4AlgP6GjWfaa2ZfN7PbwYl8APmFmLwPfBz7qNCyRBIlc9u3Kam/mb5nKleFfMNuOdHpciaSLzGgWcs5tIbSzc/x99477fh9wbWxLE4nOjiNdlBfNYUnhHK9LeYOVpfMoystix5Eu3l9bMf0PiFwgnSkqKc05x44jXUnXbgEIBIxN1cW6cLQkjAJdUtrh9j46+4bPtTeSzabqBbR09XOiV310iT8FuqS07eH+9KYk659HRH7RaJQuiaBAl5S240gXC+fnULUgz+tSJrV6cT7zcjIV6JIQCnRJWc45tjeG+udez98ylYyAUVtVxHYFuiSAAl1SVmv3ACdPDyZt/zxiU3UxDW1n6Tg75HUp4nMKdElZkVFvsvbPIyK/cOp0nVGJMwW6pKwdRzopzMuixuP5z6dzydJCcrMCartI3CnQJWXVNXezYVkRgUBy9s8jsjMDXFpeSH1zt9eliM8p0CUldfUN09jex4aqIq9LiUptZRF7j5+mf3jU61LExxTokpIio93ayuTeIRpRW1XEWNDx8tFer0sRH1OgS0qqa+4iK8O4tLzA61Kisn5Z6C+J+mb10SV+FOiSkuqbulm3tIDcrAyvS4lKYV42NQvnUac+usSRAl1SztDoGHuO9VJbmRr984jaqiJ2NXcTDGpmaYkPBbqknFeP9TI8GmRDivTPIzZUFnN6cJTX2s56XYr4lAJdUk5dU6htsSHVRujheuvUR5c4UaBLyqlv7qZqQR6l83O8LmVGKhfkUTIvm/om9dElPhToklKcc9Q3d6dcuwVCF47eUFmkHaMSNwp0SSlNnf109g1TmyInFE1UW1lMS1c/bWcGvS5FfEiBLiklMsFVqh3hEhE5s3WXRukSBwp0SSn1zd0UzMliRWlyT8g1lXVLCsjJDJzbsSsSSwp0SSn1zd2sX1aY9BNyTSU0UVcB9S0KdIk9BbqkjN6BEV5rO5tyhytOtL6yiL3HTjM4MuZ1KeIzCnRJGbvDo9rIvCipav2yIobHguw9rom6JLYU6JIydrX0EDC4rKLQ61IuSOQX0q7mHo8rEb9RoEvK2NXczapF+czNyfS6lAtSOj+HZcV5uuCFxJwCXVLCWNDx0tGelO+fR2yoLKK+pRvnNFGXxI4CXVLCoVNnODs06ptAX19ZRPuZIVq7B7wuRXxEgS4pIdKeSPUdohHrl4X2A+zS4YsSQwp0SQm7WropmZdDRfEcr0uJiYvL5jM3O0NnjEpMRRXoZrbZzA6aWYOZfWmKZd5vZvvMbK+Z/Udsy5R0tyt8QpFZap5QNFFmRoDLKgp1gpHE1LSBbmYZwH3ArcAa4C4zWzNhmRrgz4BrnXNrgc/HoVZJU51nh2jq7PdN/zxiQ2UR+0+coX941OtSxCeiGaFvAhqcc43OuWHgB8AdE5b5BHCfc64bwDnXFtsyJZ3tagkdr+23QF9fWcRY0PHyUZ1gJLERTaAvBY6Ou90avm+8i4CLzOwFM9tmZpsneyIzu9vM6sysrr29fXYVS9rZ1dJNVoaxbmmB16XE1PqK8AlGartIjEQT6JM1LScePJsJ1AA3AHcB3zKzN53O55x7wDlX65yrLS0tnWmtkqbqm7tZs6SA3KwMr0uJqYK8LFYunKcTjCRmogn0VqBi3O1y4Pgky/yXc27EOXcEOEgo4EUuyMhYkJeP9qTs/OfT2bCsiF06wUhiJJpA3wnUmFm1mWUDdwKPTVjmUeBtAGZWQqgF0xjLQiU97Tt+mqHRoG+OP59ofWUhPf0jNHb0eV2K+MC0ge6cGwXuAZ4A9gOPOOf2mtmXzez28GJPAJ1mtg94Gviic64zXkVL+oj0l9dXpvaEXFOJ7OjV8egSC1HNcuSc2wJsmXDfveO+d8Afh79EYqa+uZulhXNYXOCPE4omWl4yj4I5Wexq6eZ9tRXT/4DIeehMUUlqu5q7We/T/jlAIGCsX1aoHaMSEwp0SVrHewY43jvIhmX+bLdEbKgs4tCps/QOjHhdiqQ4Bbokrdf75/4docPrE47t1vHocoEU6JK0djX3kJsVYPXifK9LiavLKgoJ2OtnxIrMlgJdklZ9SzeXlReSleHvt+ncnExWL87XkS5ywfz9SZGUNTgyxt5jvb6bv2UqGyqL2N3SzVhQJxjJ7CnQJSntae1lNOjSKtD7hsc4ePKM16VIClOgS1KKHMZ3hU/PEJ0osmNU86PLhVCgS1Kqb+6mumQuxXOzvS4lIcqL5lA6P0d9dLkgCnRJOs456pu7fDsh12TMjNrKIuqau7wuRVKYAl2SzuH2Prr7R6itSp9Ah1Af/WjXAKdOD3pdiqQoBboknfrwKHVDZbHHlSTWxqrQ661rUttFZkeBLkmnrqmborwsVpTO9bqUhFqzJJ85WRlqu8isKdAl6dQ1d7OhshizyS6W5V9ZGQEuryjUCF1mTYEuSaXj7BBHOvrYmGb984jaqiL2nThN39Co16VIClKgS1KJjE7TbYdoRG1VMWNBx0tHNa+LzJwCXZJKfXMX2ZkB1i0t8LoUT1yxrBAz7RiV2VGgS1LZ2dTNZeUF5GRmeF2KJ/Jzs1i1KF87RmVWFOiSNAaGx9h7vDftDlecqLayiF3N3YyOBb0uRVKMAl2SxsutPYyMubTdIRpRWxWaqOuAJuqSGVKgS9KITMiVLjMsTqU2fIKRrjMqM6VAl6Sxs6mLlQvnUZiXHhNyTWVp4RwWF+Sys0l9dJkZBbokhbGgo76p+9zp7+luY1UxO4504ZwueCHRU6BLUth/4jRnhka5arkCHeDK5cW0nRmiubPf61IkhSjQJSlsa+wEYFO1Ah3gyvD/w/YjnR5XIqlEgS5JYceRLpYV57G4YI7XpSSFFaXzWDA3m+1H1EeX6CnQxXPBoGNHU5dG5+OYGZuqi9neqECX6CnQxXOvtZ2lp3/kXJtBQq6sLuZYzwCt3eqjS3QU6OK5SJ/4yuoFHleSXDaF/z92qO0iUVKgi+e2H+licUEuFcXqn4+3atF88nMzFegSNQW6eMo5x/bGUP883S5oMZ1AINxHV6BLlKIKdDPbbGYHzazBzL50nuXea2bOzGpjV6L42ZGOPjrODqndMoUrqxdwpKOPNl04WqIwbaCbWQZwH3ArsAa4y8zWTLLcfOCzwPZYFyn+FRl96giXyW06dzy6RukyvWhG6JuABudco3NuGPgBcMcky/018PeAhhIStR1HuiiZl512F4SO1tol+czNzlAfXaISTaAvBY6Ou90avu8cM7sCqHDO/fx8T2Rmd5tZnZnVtbe3z7hY8RfnHFsPd6p/fh6ZGQFqq4rZ2qgzRmV60QT6ZJ+0czMGmVkA+CrwhemeyDn3gHOu1jlXW1paGn2V4kuNHX2cPD3INStKvC4lqV2zYgENbWc5pT66TCOaQG8FKsbdLgeOj7s9H1gHPGNmTcBVwGPaMSrTebGhA4BrVyrQzyfy//Pi4Q6PK5FkF02g7wRqzKzazLKBO4HHIg8653qdcyXOuSrnXBWwDbjdOVcXl4rFN15o6GRJQS5VC/K8LiWprVmcT2FeFi80qO0i5zdtoDvnRoF7gCeA/cAjzrm9ZvZlM7s93gWKP40FHVsbO7lmZYn659MIBIyrly/gxYYOzY8u55UZzULOuS3Algn33TvFsjdceFnid/uOn6Z3YIRrV+r482hcs7KEx189SVNnP9UlOiJIJqczRcUTL4T7wdohGp1rV4R+8b3QoD66TE2BLp54oaGDlQvnUZaf63UpKaG6ZC6LC3K1Y1TOS4EuCTc0OsbOpq5zo06ZnplxzYoSth7uJBhUH10mp0CXhNvd0sPgSJBrdLjijFy7cgHd/SPsO3Ha61IkSSnQJeFebOggYHDVco3QZ0LHo8t0FOiScC8c7uSSpQUUzMnyupSUUpafy4rSuTyv49FlCgp0Saje/hF2t3TzlhpN/TAbb6kpZXtjJ4MjY16XIklIgS4J9VxDO0EHN1ysQJ+NGy4uZWg0yDZN1iWTUKBLQj1zsJ2COVlcXlHodSkp6arlC8jJDPDMQc1WKm+mQJeECQYdzxxs5y01JWRm6K03G7lZGVyzYgHPHGzzuhRJQvpUScLsO3GajrND3HDxQq9LSWk3XLyQps5+mjr6vC5FkowCXRImMqp860Xqn1+IyP4HjdJlIgW6JMwzB9u5ZGkBpfNzvC4lpVUumEt1yVyeOaQ+uryRAl0Soqd/mF0t3Tq6JUbeelEpWw/r8EV5IwW6JMRzr3XocMUYetuqhQyNBnWtUXkDBbokxDMH2ynMy+LyiiKvS/GFK6uLyc0K8BsdvijjKNAl7saCjt8cauMtNaVkBHR1olgIHb5YwpMHTukqRnKOAl3irr65m46zw9yypszrUnzl7WvKONo1wP4TZ7wuRZKEAl3i7om9J8nOCKh/HmM3ry7DLPT/KwIKdIkz5xy/ePUk19WUMD9XsyvGUun8HDZWFivQ5RwFusTV3uOnOdYzwDvWqt0SD7esLePAyTM6a1QABbrE2RN7TxKwUHtAYu8daxcBartIiAJd4uqJvSfZWFXMgnk6OzQeKorzWLskX4EugAJd4qix/SyHTp1l87pFXpfia5vXLmJXSw+nTg96XYp4TIEucfPE3lPA620BiY/IL8xf7jvlcSXiNQW6xM0v9p7k0vIClhTO8boUX1u5cB7LS+byi1dPeF2KeEyBLnHR3NnHy0d7uHXdYq9L8T0z49ZLFrH1cCdtZ9R2SWcKdImLR3cfB+COy5d4XEl6+J3LlxJ08LOXNUpPZwp0iTnnHI++dIyrlher3ZIgNWXzWbc0n0d3H/O6FPGQAl1i7qWjPRzp6OPdV5R7XUpaedcV5bxyrJeGNs3tkq6iCnQz22xmB82swcy+NMnjf2xm+8xsj5k9aWaVsS9VUsWju4+Rkxlg8yU6uiWRbrtsMQGDn2qUnramDXQzywDuA24F1gB3mdmaCYvtBmqdc5cCPwb+PtaFSmoYGQvysz0nuHlNGfmauyWhFs7P5S01pTy6+zjBoKbUTUfRjNA3AQ3OuUbn3DDwA+CO8Qs45552zvWHb24D9Ld2mnr2UDtdfcO86/KlXpeSlt51xVKO9Qyws6nL61LEA9EE+lLg6LjbreH7pvJx4PHJHjCzu82szszq2tt1pRU/+unuYxTlZfFWTZXriVvWlpGXncGjL6ntko6iCfTJLjEz6d9zZvYhoBb4ymSPO+cecM7VOudqS0v1gfebnv5hfrXvFLddtoSsDO1v90Jediab1y3i5y+foH941OtyJMGi+dS1AhXjbpcDxycuZGY3A38O3O6cG4pNeZJKflTXytBokLs2LfO6lLT2u5uWcWZolP966U0fU/G5aAJ9J1BjZtVmlg3cCTw2fgEzuwL4JqEwb4t9mZLsgkHHv29vZmNVEasX53tdTlrbUFnEqkXzeWhrs643mmamDXTn3ChwD/AEsB94xDm318y+bGa3hxf7CjAP+JGZvWRmj03xdOJTz77WTnNnP793dZXXpaQ9M+PDV1ex/8Rp6pu7vS5HEigzmoWcc1uALRPuu3fc9zfHuC5JMQ9vbaZkXg6bNbNiUrjj8iX8zZb9PLytmdqqYq/LkQTRniu5YEe7+nnqYBt3baogO1NvqWQwNyeT92woZ8srJ2g/o11a6UKfPrlg39vegoF2hiaZD11VyciY44c7W7wuRRJEgS4XZGB4jEfqjvL2NWWaiCvJrFw4j2tXLuB721sYGQt6XY4kgAJdLsj3d7TQ1TfMx69b7nUpMomPX1fNid5BfrpLJxqlAwW6zNrgyBjffPYwV1YXs6laO96S0dsuXsi6pfl8/ZkGRjVK9z0Fuszaj+pbOXV6iM/eVON1KTIFM+Oet9XQ1NnPz/fo4hd+p0CXWRkeDfKNZw6zflkh16xY4HU5ch63rCnj4rL5/MvTDZqF0ecU6DIrP93dyrGeAT5zUw1mk033I8kiEDDuuXElDW1nefzVk16XI3GkQJcZGxkL8vVnDnNpeQE3XKRJ1lLBOy9ZzPLSufzzU69plO5jCnSZsYe3NtPc2c/nb9boPFVkBIzP3VTDgZNn+PGuVq/LkThRoMuMdJ4d4qu/PsT1F5XytosXel2OzMDtly1hQ2URf/+Lg5wZHPG6HIkDBbrMyD/88hADw2Pc+9urNTpPMWbGX9y2hs6+If75qQavy5E4UKBL1F491ssPdrbw4aurWLlwvtflyCxcWl7I+zaU8+0XjtDYftbrciTGFOgSFeccX/7ZPorysvnczTruPJV98R2ryMnM4K9/vk/zpfuMAl2i8u/bmtnR1MUX33ExBXOyvC5HLkDp/Bw+f3MNTx9s17VHfUaBLtNqaDvL/96yn7deVMqdGyum/wFJeh+7tprayiLufXQvR7v6vS5HYkSBLuc1PBrk8z/czZysDL7y3ku1I9QnMgLGVz9wOQ74wiMvM6Zj031BgS7n9bUnD/HqsdP8zbsvZWF+rtflSAxVFOfxV7evZUdTF9989rDX5UgMKNBlSk8faOP+Zw7zgdoKNq/TpeX86N3rl/JblyzmH395iK2HO70uRy6QAl0mtfd4L/f8xy5WL87n3tvWeF2OxImZ8X/efQlVJXP55MN1NLTpUMZUpkCXNznRO8Dvf2cnBXOyePCjG5mbE9W1xCVFFczJ4tsf3Uh2ZoCPfWcHHWd1DdJUpUCXN+jtH+H3v1NH39AYD35sI2Xqm6eFiuI8vvWRjbSdHuIPvlvH2aFRr0uSWVCgyzltZwb5wANbOdx2lq9/cD2rFuV7XZIk0OUVhXztzit45VgvH/zXbXT3DXtdksyQAl0AONrVz/u+sZWWrn4e/OhGrte0uGlp87pFfONDG9h/8gzv/+ZWTvYOel2SzIACXdjV0s177n+Rnv4RvvcHV3JdTYnXJYmH3r6mjO9+bBMnegd5z/0v8uqxXq9Lkigp0NNYMOi4/5nDvO8bW8nODPDIJ6/mimVFXpclSeDqFQv4/ieuYjQY5N1ff5HvvHBE876kAAV6mmrt7ucj397B3/3iAJvXLuK/P/sWLl6kGRTldZeUF/D4567nLTUl/OXP9vGJh+rVgkly5tVv3draWldXV+fJutNZ39Ao9z9zmAeeayRg8Be3reXOjRU6pV+m5Jzj2y808bePHyAjYHzqhhV84vrl5GZleF1aWjKzeudc7aSPKdDTQ+/ACD/c2cK3njtC25khfufyJfzpratYXDDH69IkRbR09vM3j+/n8VdPsqQgl09cv5z31VYwT+cpJJQCPU0559jT2stPdrXy4/pW+ofHuGp5MX+yeRXr1SuXWdrW2MlXnjhIfXM383MyeV9tBe/ZsJQ1i/P1l14CKNDTSP/wKLuae3j2tXa2vHKC1u4BsjMC3HbZEj52bRXrlhZ4XaL4xEtHe3jw+SNseeUEo0FH1YI8br1kMW+pKWH9siK1ZOLkggPdzDYDXwMygG855/52wuM5wEPABqAT+IBzrul8z6lAv3BnBkc40tHHgRNn2HfiNHtae9jT2sto0JGVYVy3soRbL1nMLWvKKMzL9rpc8anOs0P8ct8ptrxyghcPdzIWdGRnBLi8opB1SwtYvXg+qxfnU10yV9NIxMAFBbqZZQCHgLcDrcBO4C7n3L5xy3wKuNQ594dmdifwLufcB873vAr0kLGgY2QsyNBIkKHRMQZGxugbGqN/eJQzQ6OcHhihd2CE7r4R2s8O0n5miJO9g7R09dPd//qV2/OyM1i9OJ+NVcVctbyY2qpi9TYl4U4PjlDX1MW2xi52HOniwMnTDI4Ezz2+YG425cV5LM7PpXR+DqXzcyjKyyJ/ThYFc7KYn5tJXnYmc7Mzyc0KkJOZQU5WgKyMABkBtXPg/IEezSd+E9DgnGsMP9kPgDuAfeOWuQP4y/D3Pwb+xczMxaGf88jOo/zrc40z+pkLKWL8S3jD87jX/3HOhf8Fhwv960L3Bx0EnSPoHGPB179Ggo7RsSAzua5AUV4WpfNzKMvP5dZLFrOsOI+qBXmsWpTPsuI8AnrDi8fyc7O4cVUZN64qA0IDlqbO0F+RzV19HO0a4GhXP4fbz7LtSCc94wYl0wkYZGUEyAwYGeO+Ahb5Cs0eaUboi/D3hO+PPNG4j8n4T8yF9P9n+pOfvamG2y5bMuv1TSWaQOJnRToAAAVjSURBVF8KHB13uxW4cqplnHOjZtYLLAA6xi9kZncDdwMsW7ZsVgUX5mVRUzZvxj9nM/4vf8MPT/btuTdA6A3zxjdO5I2VYUYgELqdYa+/CTMDRnZmgMxAgKxMIzc8EsnNzGBuTiZzczLIy86kIDxyKZiTRXamThuQ1JIRMFaUzmNF6eSf2aHRMXoHRjg9MErvwAh9Q6P0D4/SNzTG4OgYQyNBBkfHGBkN/SU7MhZkNDwoGg2GBkTOOYJBGHPu3EDqzQOtkOkGaLPhZvHD8boubzSBPlkSTnwF0SyDc+4B4AEItVyiWPeb3LJ2Ebes1cUWRPwgJzODhfMzWKhz2mIimiFfKzD+ysDlwPGpljGzTKAA6IpFgSIiEp1oAn0nUGNm1WaWDdwJPDZhmceAj4S/fy/wVDz65yIiMrVpWy7hnvg9wBOEDlt80Dm318y+DNQ55x4D/g142MwaCI3M74xn0SIi8mZRHdfmnNsCbJlw373jvh8E3hfb0kREZCZ02ISIiE8o0EVEfEKBLiLiEwp0ERGf8Gy2RTNrB5pn+eMlTDgLNQ3oNacHveb0cCGvudI5N+lV3D0L9AthZnVTTU7jV3rN6UGvOT3E6zWr5SIi4hMKdBERn0jVQH/A6wI8oNecHvSa00NcXnNK9tBFROTNUnWELiIiEyjQRUR8IuUC3cw2m9lBM2swsy95XU88mFmFmT1tZvvNbK+ZfS58f7GZ/crMXgv/W+R1rbFkZhlmttvMfh6+XW1m28Ov94fh6Zt9w8wKzezHZnYgvK2vToNt/Efh9/SrZvZ9M8v123Y2swfNrM3MXh1336Tb1UL+KZxne8xs/YWsO6UCPXzB6vuAW4E1wF1mtsbbquJiFPiCc241cBXw6fDr/BLwpHOuBngyfNtPPgfsH3f774Cvhl9vN/BxT6qKn68Bv3DOrQIuI/TafbuNzWwp8Fmg1jm3jtB03Hfiv+38HWDzhPum2q63AjXhr7uB+y9kxSkV6Iy7YLVzbhiIXLDaV5xzJ5xzu8LfnyH0QV9K6LV+N7zYd4Hf8abC2DOzcuC3gG+FbxtwI6GLjoP/Xm8+cD2hawngnBt2zvXg420clgnMCV/ZLA84gc+2s3PuWd58xbaptusdwEMuZBtQaGaLZ7vuVAv0yS5YvdSjWhLCzKqAK4DtQJlz7gSEQh9Y6F1lMff/gD8BguHbC4Ae59xo+LbftvVyoB34drjN9C0zm4uPt7Fz7hjwD0ALoSDvBerx93aOmGq7xjTTUi3Qo7oYtV+Y2TzgJ8DnnXOnva4nXszst4E251z9+LsnWdRP2zoTWA/c75y7AujDR+2VyYT7xncA1cASYC6hlsNEftrO04np+zzVAj2aC1b7gpllEQrz7znn/jN896nIn2Phf9u8qi/GrgVuN7MmQm20GwmN2AvDf5qD/7Z1K9DqnNsevv1jQgHv120McDNwxDnX7pwbAf4TuAZ/b+eIqbZrTDMt1QI9mgtWp7xw//jfgP3OuX8c99D4i3F/BPivRNcWD865P3POlTvnqght06eccx8EniZ00XHw0esFcM6dBI6a2cXhu24C9uHTbRzWAlxlZnnh93jkNft2O48z1XZ9DPhw+GiXq4DeSGtmVpxzKfUFvBM4BBwG/tzreuL0Gq8j9GfXHuCl8Nc7CfWVnwReC/9b7HWtcXjtNwA/D3+/HNgBNAA/AnK8ri/Gr/VyoC68nR8Fivy+jYG/Ag4ArwIPAzl+287A9wntIxghNAL/+FTblVDL5b5wnr1C6AigWa9bp/6LiPhEqrVcRERkCgp0ERGfUKCLiPiEAl1ExCcU6CIiPqFAFxHxCQW6iIhP/H8U8Ro3JwE1iwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Let's create a signal and its FT (ie: data, with noise)\n", "N = 101\n", "rmax = 40\n", "r = np.linspace(-rmax,rmax,N)\n", "x = 1.0*(np.exp(-(r*r/100.0)))\n", "\n", "# Calculate a (discrete) Fourier transform, with some noise\n", "xn = x + 0.04*rnd.normal(loc=0.0, scale=1.0, size=N)\n", "d = np.fft.fftshift(np.fft.fft(np.fft.ifftshift(xn)))/np.sqrt(N)\n", "\n", "plt.plot(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Unconstrained problem\n", "\n", "Suppose we are reconstructing an image $\\mathbf{x}$ from some data $\\mathbf{d} = \\mathbf{Ax}$. In this case, $\\mathbf{A}$ calculates a (discrete) Fourier transform. Our unconstrained optimization problem is:\n", "\n", "$$\\hat{\\mathbf{x}} = \\arg \\min_{\\mathbf{x}} \\| \\mathbf{Ax} - \\mathbf{d} \\|_2^2 $$\n", "\n", "Importantly, even though our formulation is based on a matrix-vector multiplication where $\\mathbf{A}$ is the DFT matrix, in practice we calculate $\\mathbf{A x}$ by taking an FFT of $\\mathbf{x}$. " ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [], "source": [ "def func(x1):\n", " # Note that the matrix A is implemented efficiently using FFT\n", " d1 = np.fft.fftshift(np.fft.fft(np.fft.ifftshift(x1)))/np.sqrt(N) \n", " value = la.norm(d1-d)**2\n", " return value" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "`gtol` termination condition is satisfied.\n", "Number of iterations: 3, function evaluations: 306, CG iterations: 2, optimality: 6.89e-07, constraint violation: 0.00e+00, execution time: 0.15 s.\n", "\n", "Final value of the cost function:\n", "2.277815493144445e-12\n", "\n", "Error with respect to original noiseless image:\n", "0.10361507228999815\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3ic1ZX48e8d9d67rGJJLnK3hTG2AQMBbFOc0A1sIJAQdiHZbEISEpYUkl0ISQhhw4bwA0JbejXVNGPAXe6WZNkqVu+9t7m/P2ZGqIwsWZ7RFJ3P8/hB8847ozOMdHTn3HPvq7TWCCGEcH0GRwcghBDCNiShCyGEm5CELoQQbkISuhBCuAlJ6EII4SY8HfWNIyMjdUpKiqO+vRBCuKS9e/fWa62jrN3nsISekpJCdna2o769EEK4JKVUyVj3SclFCCHchCR0IYRwE5LQhRDCTUhCF0IINyEJXQgh3IQkdCGEcBOS0IUQwk1IQhdiDJ29/Ty34wTdfQOODkWICXHYwiIhnFnfgJF/fX4fW4/VEeznxYbFCY4OSYhxyQhdiBGMRs1PXz3I1mN1eBoUe040OjokISZERuhCDKG15r/ez+OtA5X89OLZ7CxqIPtEk6PDEmJCxh2hK6WeUkrVKqWOjHG/Uko9opQqUEodUkottX2YQtjWlqO1rPnjFo5Wtw47/vKeMp78qpibV6bwb2vSOCMlnPyaNlo6+xwUqRATN5GSy9PA2pPcvw7IMP+7Dfj76YclhP309A/w6005nGjo5I7/20dHTz8ABbVt/OadHFanR/KrSzNRSpGVEobWsLdUyi7C+Y2b0LXWXwAn+2neADyrTXYCoUqpOFsFKIStPbejhNLGTu48L53i+g7uefMw3X0D/ODFA/h7e/LQNYswGBQAS2aEmevoUnYRzs8WNfQEoGzI7XLzsaqRJyqlbsM0iicpKckG31qIU9PU0csjnx7n3FlR3HXxbLw9DTz08TGK6jvIq2rlyZuyiA72HTzfz9uD+QkhZMvEqHABtuhyUVaOaWsnaq0f11pnaa2zoqKs7s8uhF399dPjtPf0c88lcwG447x0zs6I5FB5CzevTOGCuTGjHrM8NZyDZS3Sjy6cni0SejkwY8jtRKDSBs8rhE0V13fw/M4SrluexKyYIAA8DIpHrlvCfRvmcfe6OVYfl5UcRu+AkcMVLVMZrhCnzBYJfRPwbXO3ywqgRWs9qtwihKN9nFtNv1Hz7xdkDDseFuDNt89KwdfLw+rjslLCAaQfXTi9cWvoSqkXgTVApFKqHPg14AWgtX4MeB9YDxQAncB37BWsEKejpKGTMH8vYobUyCciPMCb9OhA9hQ3mn4ThHBS4yZ0rfXGce7XwB02i0gIOylt7CQp3H9Sjz0jJYx3D1VhNOrBDhghnI0s/RfTRmljJ0kRAZN67Bkp4bR195Nf02bjqISwHUnoYlroHzBS0dRF8iRH6KvSI/HyUDy/c8wLrgvhcJLQxbRQ2dxNv1FPuuQSE+zLxuVJvLynjNKGThtHJ4RtSEIX00JpoykJJ0VMLqED3HleOp4eioc/OWarsISwKUnoYlooaewAmPQIHSA62JebVqbw5oEKjkktXTghSehiWiht6MTb00DsKbYsjnT7OWkEenvy54/ybRSZELYjCV1MC6WNncwI8zvtlsOwAG++e/ZMNufU8OiWApo6em0UoRCnTxK6mBZKGibfgz7SrWencnZGJH/cnM+K+z/l568doqG9xybPLcTpkCsWCbentaa0sZPlqeE2eb5AH0+eu/VMjla38sz2El7fW05FcxfP3rJcFh0Jh5IRunB7TZ19tPf0M8NGI3SLObHB3H/FAu7bMI+vCup54qsimz6/EKdKErpweyUNpg6XyS4qGs+1Z8xg7bxY/rg5n8PlsiOjcBxJ6MLtWXrQk0+jB/1klFI8cOUCIgJ8+OFL+wcvaSfEVJOELtyeZWWnrUsuQ4X6e/OXaxdTXN/By3vKxn+AEHYgCV24vZLGTmKCfcbc79xWzkqLICLAWxYdCYeRhC7c3ulsm3uq0qICKaxrn5LvJcRIktCF2ytt6CQpfHLb5p6qtOgAiuo6puR7CTGSJHTh1rr7Bqhu7bbbhOhIaVGBNHT0ygpS4RCS0IVbK28y77I4RSWXmVGmTwJF9VJ2EVNPErpwayUNp79t7qlIiwoEoLBWyi5i6klCF24hr6qVJ74swnSJ268dLGsG7LeoaKTEMH+8PQwyMSocQvZyEW7hpd2lPLOjhCBfT649IwmAE/Ud/OOLIi7MjCEi0GdK4vAwKFIjAyShC4eQEbpwC3Xm3Q5/+04uJQ0daK355ZuH8fYw8LsN86c0Ful0EY4iCV24hfq2XtKiAvA0KH708gFe2lPG9sIGfr5uDrEhp3dRi1M1MzKQksZOevuNU/p9hZCELtxCfXsPc2KD+e8rFrC/tJlfvnmY5SnhXL88acpjSYsOYMCoKW2UUbqYWpLQhVuoa+8hMtCbSxfGc8XSBHw8Dfz3FQscsj+5pdOlQDpdxBSTSVHh8rr7Bmjr7icqyDTx+aerFnHP+rlTNhE60kxL66JMjIopJiN04fLqzROikeYEbjAohyVzMF3RKDbYVxK6mHITSuhKqbVKqXylVIFS6m4r9ycppbYopfYrpQ4ppdbbPlQhrKtvNy2zj3RgEh9pZpR0uoipN25CV0p5AI8C64BMYKNSKnPEaf8JvKK1XgJcB/yvrQMVYiz1beYRepDzJHTLrosjFzoJYU8TGaEvBwq01kVa617gJWDDiHM0EGz+OgSotF2IQpzc1yUXbwdH8rW0qADauvsH++OFmAoTSegJwNBLsJSbjw31G+BGpVQ58D7wA2tPpJS6TSmVrZTKrqurm0S4Qow2sobuDNKiZU8XMfUmktCt9X2N/By5EXhaa50IrAeeU0qNem6t9eNa6yytdVZUVNSpRyuEFfXtvQT5eNr9ikSnYrB1USZGxRSaSEIvB2YMuZ3I6JLKrcArAFrrHYAvEGmLAIUYT117z2DLorOIDfYlLsSXD49UOToUMY1MJKHvATKUUqlKKW9Mk56bRpxTClwAoJSaiymhS01FTIm6th6nKreAqXXyxhXJbCtoIL9arjEqpsa4CV1r3Q/cCWwG8jB1s+Qope5TSl1uPu0nwPeUUgeBF4GbtUzviylS395DZJDzTIhaXL88CR9PA09vL3Z0KGKamNBKUa31+5gmO4ce+9WQr3OBVbYNTYiJqW/rITLd+Sp8YQHefGtJAm/sq+BnF88hLMD5/ugI9yIrRYVL6+kfoLW73+lKLhbfWZVKT7+RF3aXOjoUMQ1IQhcurcEJV4kONTs2iFXpETy3o4S+AdlOV9iXJHTh0pxxUdFI31mZSnVrN5tzqh0dinBzktCFS7MkdGdrWxzq/DnRhPl7sa2g3tGhCDcnCV24tPo25y65gKmFMSM6iOM1sshI2JckdOHS6lxghA6mrQAKZLMuYWeS0IVLq2vrIdDJlv1bkxEdSHNnHw0dvY4ORbgxSejCpdWbLz3n7NLNm3VJ2UXYkyR04dJMCd25yy3wdUKXzbqEPUlCFy6tvr3XJRJ6XIgvAd4eFNZKQhf2IwlduDRn3cdlJKUU6dGBFEhCF3YkCV24rN5+I82dfUQF+jo6lAlJiw7keK3svCjsRxK6cFkNHZZriTr/CB1MdfSa1h5au/scHYpwU5LQhctyhUVFQ6VHWS5LJ2UXYR+S0IVLGbowxxmvJXoyGTFBAByXhC7sRBK6cCkPfXyMVQ98xtZjdV+vEnWRhD4jzA9vD4OM0IXdTOgCF0I4i63H6qho7uKmp3aTGhkAuE4N3dPDQGpkgHS6CLuRhC5cRv+AkfzqNv5lRTLengae/KqYAG8P/L1d58c4PTqQI5Utjg5DuCnX+U0Q015xfQc9/UYWzwjlymWJXJQZQ0uXa3WMpEcH8v6RKrr7Bpx+/xnheiShC5eRW9UKQGZ8MABnzoxwZDiTkh4diNZQVNcx+DqEsBWZFBUuI7eqFW8PA2nm9j9XJHu6CHuShC5cRm5lK+nRgXh7uu6PbWpkAAYF+dWtjg5FuCHX/c0Q005eVavLlyl8vTxYnhrOa3vL6e4bcHQ4ws1IQhcuobatm/r2XubGuXZCB/jB+RnUtPbw6t5yR4ci3IwkdOEScivNE6JukNBXpkWwLDmMv28poLff6OhwhBuRhC5cQl6VaZdCd0joSil+eEEGlS3dvL5PRunCdiShC5eQW9VKQqgfIf5ejg7FJs7JiGTRjFAe3VJA34CM0oVtTCihK6XWKqXylVIFSqm7xzjnGqVUrlIqRyn1gm3DFNNdbmWLW9TPLZRS/OiCDMqbunhzf4WjwxFuYtyErpTyAB4F1gGZwEalVOaIczKAXwCrtNbzgB/ZIVYxTXX1DlBc734LcdbMjiIlwp+PcqodHYpwExMZoS8HCrTWRVrrXuAlYMOIc74HPKq1bgLQWtfaNkwxneXXtGHUkBkX5OhQbEopxYLE0MH5ASFO10QSegJQNuR2ufnYULOAWUqpbUqpnUqptdaeSCl1m1IqWymVXVdXN7mIxbSTZ1nyHxfi4Ehsb05sEBXNXXIVI2ETE0noysoxPeK2J5ABrAE2Ak8opUJHPUjrx7XWWVrrrKioqFONVUxTuZWtBPp4khjm5+hQbM7StXNURunCBiaS0MuBGUNuJwKVVs55W2vdp7UuBvIxJXghTovWmq8K6lmSFIrBYG1s4drmmMtIR2UrAGEDE0noe4AMpVSqUsobuA7YNOKct4DzAJRSkZhKMEW2DFRMT8dr2ymu7+CiebGODsUuYoN9CfX3GiwrCXE6xk3oWut+4E5gM5AHvKK1zlFK3aeUutx82magQSmVC2wBfqq1brBX0GL62HzE1AFyUWaMgyOxD6UUc2KDZGJU2MSE9kPXWr8PvD/i2K+GfK2BH5v/CWEzH+XWsCQplJhgX0eHYjdz44J5aXcZA0aNhxuWlcTUkZWiwmlVNHdxuKKFi9203GIxNzaYrr4BShs7HR2KcHGS0IXTsiy4cfuEbu50kTq6OF2S0IXT2pxTzayYQFIjAxwdil1lxARiUHBUEro4TZLQhVNq7Ohld3Gj24/OwXTRi5lRgeTKxKg4TZLQhVP6JK8Go3b/covFnNgg6UUXp00SunBKn+XVkhDqxzw325BrLHPjgilvki0AxOmRhC6c0uGKFrJSwlBqerTxzTWvGM2vlrKLmDxJ6MLpNHf2UtHc5RZXJ5oo6XQRtiAJXTidXHNSc6cLWownNtiXED/ZAkCcHknowulYLgg9nRK6UorZMUEcr2l3dCjChUlCF04nr6qN6CAfooJ8HB3KlEqPCeR4bTumnTSEOHWS0IXTya1qdbvLzU1EelQgLV191Lf3OjoU4aIkoQun0ttvpKC2bVqVWywyYgIBOF4rnS5iciShC6dyvLaNvgE9rTpcLNKjTQm9sFbq6GJyJKELp2KZEJ2OJZfYYF8CfTw5LgldTJIkdOFU8qra8PPyICXCvTfkskYpRVp0IAWS0MUkSUIXTiW3qoXZsUHT9kIPGdGBMkIXkyYJXTgNrTW5ldOzw8UiPTqQurYeWjplTxdx6iShC6dR0dxFa3f/tJwQtcgwT4wW1Emnizh1ktCF07BcKHm6j9ABqaOLSZGELpxGbmUrSpn2Bp+uEsP88fE0yBYAYlIkoQunkVvVQmpEAP7eno4OxWE8DIqZUYEU1ElCF6dOErpwGjmVrcydxuUWi4zoQBmhi0mRhC6cQktnH+VNXcyPD3F0KA6XHh1IRXMXnb39jg5FuBhJ6MIp5FS2AEybS86dTMbgFgAdDo5EuBpJ6MIp5JiX/EtCH9LpUtfG7uJGLvjz5/zklYMOjkq4guk7+yScypHKFuJCfIkInF57oFuTHBGAp0Hx8CfHKW3sRGvoHTA6OizhAiY0QldKrVVK5SulCpRSd5/kvKuUUloplWW7EMV0kFPZyjypnwPg7WlgZlQApY2d3HRWCreuTqWyuZs+SepiHOOO0JVSHsCjwIVAObBHKbVJa5074rwg4IfALnsEKtxXZ28/hXXtXLIgztGhOI2Hr13CgFGzIDGEV7PLGDBqKpq6SImcfpuWiYmbyAh9OVCgtS7SWvcCLwEbrJz3O+BBoNuG8YlpIK+qDa1hfoKM0C0y44NZkGj6/5Fs3nmypLHTkSEJFzCRhJ4AlA25XW4+NkgptQSYobV+92RPpJS6TSmVrZTKrqurO+VghXuSDpeTSwr3B6BUEroYx0QSurV9TAevYquUMgB/AX4y3hNprR/XWmdprbOioqImHqVwazkVrYT5exEX4uvoUJxSdJAPPp4GShukjVGc3EQSejkwY8jtRKByyO0gYD7wuVLqBLAC2CQTo2Isz+44wfuHqwZvH6lsYX5CCEpNzz3Qx2MwKJLC/SlpkBG6OLmJtC3uATKUUqlABXAdcL3lTq11CxBpua2U+hy4S2udbdtQhbt45NPjdPQMMD8+hNgQX47VtHHL6lRHh+XUkiP8peQixjXuCF1r3Q/cCWwG8oBXtNY5Sqn7lFKX2ztA4V46evqpb++lq2+AX7x5iGM1potCy5L/k0sKDzD3pOvxTxbT1oQWFmmt3wfeH3HsV2Ocu+b0wxLuqqzJNMpcnR7JVwX1/P49U/erTIieXHKEP529A9S19xAdJHMNwjpZKSqmVKm5DnzXxbPpGzCys6iRAO/peVHoU5EUYe50aeiUhC7GJHu5iCllqQOnRPjzwJUL8fE0kBkfjGGaXhR6opLNrYsyMSpORkboYkqVNXYS5OtJiJ8Xof7e/PPmMwj283J0WE4vMcwfpaQXXZycJHQxpUobO0kK9x9sUVyZHjnOIwSY9neJD/GThC5OSkouYkpZEro4daZedFlcJMYmCV1MGaNRU9bUJQl9kqQXXYxHErqYMjVt3fT2G5khCX1SkiL8qW/vpb1HLk0nrJOELqaMpWVRRuiTkxxuau0slU4XMQZJ6GLKWMoFktAnJ9nSi94odXRhnSR0MWXKGjsxKIgP9XN0KC7JsrhIetHFWCShiylT2thJXIgf3p7yYzcZwb5ehPl7ycSoGJP0oYspIy2Lpy8pIoAPjlRTWNeOUcPZ6ZH84IIMR4clnIQMlcSUKW2UlsXTtfGMGaRFBWDUUNncxaOfF9DbLxePFiYyQhdTorO3n/r2nsE6sJic65Yncd3yJAA+PFLF7c/v43BFC8uSwxwcmXAGMkIXU6KssQtAetBtaHlqBAC7ihscHIlwFpLQxZSQlkXbCw/wZlZMILuKGh0dinASktDFlJCEbh/LU8PJPtFI/4DU0YUkdDFFyho7CfTxJMxftsq1pTNTI+joHSC3qtXRoQgnIAldTInSxk5mDNk2V9jGmanhAFJ2EYAkdDFFSho6SAqXFaK2Fh3sS2pkgEyMCkASupgCLZ19FNV3MDdOLgRtD2emhrO7uJEBo3Z0KMLBJKELu9t9ohGt4ayZEY4OxS2dOTOc1u5+8qvbHB2KcDBJ6MLudhY14ONpYNGMUEeH4pakH11YSEIXdrezqIGlSWH4enk4OhS3lBDqR2KYn0yMCknowr5aOvvIrWplhZRb7GrFzAh2FDVIP/o0Jwld2JWlfr5iZrijQ3Fr35gbTUtXH3tONDk6FOFAktCFXUn9fGqcMysKH08DH+VWOzoU4UATSuhKqbVKqXylVIFS6m4r9/9YKZWrlDqklPpUKZVs+1CFK5L6+dTw9/bk7IxIPsqpQWtpX5yuxk3oSikP4FFgHZAJbFRKZY44bT+QpbVeCLwGPGjrQIXrkfr51LooM5aK5i7ZBmAam8gIfTlQoLUu0lr3Ai8BG4aeoLXeorW2XBdrJ5Bo2zCFK5L6+dS6YG40BgUf5dQ4OhThIBNJ6AlA2ZDb5eZjY7kV+MDaHUqp25RS2Uqp7Lq6uolHKVyS1M+nVkSgD1nJ4XyUKwl9uppIQre2m5LVIp1S6kYgC/ijtfu11o9rrbO01llRUVETj1K4FK01h8tb+CSvRurnU+yieTHkVbVSJheSnpYmktDLgRlDbicClSNPUkp9A7gHuFxr3WOb8IQr6R8w8vzOEtb99Usu+9tXVLd0c+MKmR+fShdmxgDIKH2amsg1RfcAGUqpVKACuA64fugJSqklwD+AtVrrWptHKZze4fIWfvHmIY5UtLIgIYTffXM+ly+KJ8RP9j+fSskRAcyJDWJzTjW3rk51dDhiio2b0LXW/UqpO4HNgAfwlNY6Ryl1H5Cttd6EqcQSCLxq3u+6VGt9uR3jFk7koY/y+duWAiICffjfG5aybn6s7HvuQJctiuePm/M5UNbMYpm/mFaUo3pWs7KydHZ2tkO+t7CdqpYuzrr/M9YviOX+KxbKiNwJtPf0c+6DW8iICeTF762QP65uRim1V2udZe0+WSkqTsveEtNS89vPTZNk7iQCfTz5wfnp7CxqZOsx6SabTiShi9Oyt6QJXy+DXLzCyVx/ZjIzwv34w4f5GOXCF9OGJHRxWvaVNLEoMRQvD/lRcibengbuumg2eVWtbDo4qilNuCn5LRST1t03QE5lK8uSwxwdirDisoXxzIsP5g8fHqWyucvR4YgpIAldTNqh8hb6jVoSupMyGBT3X7GA9p5+rn5sByUNHY4OSdiZJHQxaZYJ0SVJktCd1cLEUF783go6e01J/XiNXHfUnUlCF5O2t6SJmZEBhAd4OzoUcRLzE0J4+ftnoYHrHt/JiXoZqbsrSehiUrTW7CttYqmUW1zCrJggXr5tBUatufmfu2lol9053JEkdEFZYycbH995ShNnJxo6aezolfq5C5kZFcgTN51BVUs33302m67eAUeHJGxMErrg4U+Os6OogW0F9RN+jKV+LgndtSxLDuOv1y3hQFkzP3p5v1zdyM1IQp/mSho6eOtABQDHTmHCbG9JE0G+nqRHBdorNGEna+fHctdFs9mcU8PB8hZHhyNsaCK7LQo39rfPCvA0KOJCfMmvaR/zvI6efu58YR+h/t4sTQ4bvFaowSD7hLiifzkrmUc+Pc6b+8plAy83IiP0aay0oZM39lewcXkSy1PCOVY99gj9lewytuTX8cWxOu596wjF9R2ckSLlFlcV7OvFNzJjeOdQFX0DxsHjh8qb+fZTu+ns7XdgdGKyJKFPY49uKcDDoPjXNWnMig2iurWblq6+UecNGDVPbStmWXIY2f/5Db746Xk8duMybl4l+227siuWJNDY0cvWfNMGXlprfr0phy+O1ZF/kj/uwnlJQp+mKpq7eH1fOdcvTyIm2JfZMUEAVheebM6ppqyxi++dnYpSiqQIf9bOjyXQRyp2ruycWVFEBHjz5n7THMrmnBr2lzYDUNXS7cjQxCRJQp+E/gEj33hoKy/sKnV0KJO2raCefqPmxhVJAMyKNSX0fCsJ/Ykvi0gK9+fCzNgpjVHYl5eHgcsWxfNxXg1NHb08uPkoCaF+ALL3i4uShD4JRypbKaht56sC191rOq+qFT8vD1IjTV0q8SG+BPp4jqqj7y1pYl9pM7esSsFDJkDdzreWJNDbb+T7z++lqK6Dey/NxM/LQ0boLkoS+iRY+rWPunCdMa+qldmxQYNJWilFRkzgqBH6k18VEezrydVZM6w9jXBxCxNDmBkVwO7iRpYmhXLxvBjiQn2papERuiuShD4JloR+or6D7j7XW22nteZodRtz44KGHZ8dE0R+ddvgYpPypk4+PFLNDSuSCZB6uVtSSnHl0kQAfr52Dkop4kP8qGiWEborkoR+irr7BsguaSIh1A+jhsK6sXu3ncG+0iau+ccO2rq/7l6pbu2mubNv1FWGZsUE0dTZR317LwCv7S1HAzecmTSVIYsp9t2zU3n7jlWcOTMCgPhQX6qkhu6S3CahHyhr5t63jtDc2WvX77O3pInefiM3rUwGcGh7V2FdO89sP3HScx744Ci7ixvZVtAweCyvqhVgVEKfHft1p4vRqHk1u5zV6ZEkhvnbNnDhVHw8PVg0ZHFRXIgfde099PYbT/IoMZ7atm4Gpvjyf26R0D84XMW1/9jBcztL+PZTu2ntHt1LbSvbCurxNCiuzUrC28NgtStkMlq7+7jvnVzufevIhH8IHvroGL/elENZY6fV+3cVNbC7uBGAHYVf79OSV2WK2ZLALWbFfN3psqOogYrmLqmdT0Pxob5oDTWtUnaZrE0HK1n+X5+y7Pcfc8cL+3g1u2xK/kC6dELXWvOPrYX82wv7mBcfzJ+uXkRuZSvf+eceOnrss9JtW2EDi2aEEuLvRVp04GmP0LXWfHikigsf2so/txfz3M4Sfv9e7riPa+vu45O8GgC+PG59U62/bSkgMtCb5anhbC8cPkJPDPMj2Ndr2PmRgd6EB3hzrKaNV7LLCPb15KLMmNN4dcIVxYXYtnWxu2+Abz66jc+O1oy6r7O3f8IXsdZac8vTe/jV20dsEpe97C1p4q5XD7J4RigXzIlhT3EjP33t0JRc29WlE/rr+yq4/4OjrJ8fxwvfW8FVyxL5n42mneRufWYPPf2Tm7Dcc6LR6mNbuvo4XN7MqjRTrXF2TOBJl8tPxIOb87n9+X1EBPjw9h2ruHV1Kv/cdoInvyo+6eM259TQ02/Ex9PAF8dGt08eKGvmy+P1fO/smZw/J5rjte3UtplGXHlVraPKLWDudIkOJPtEEx8eqeabSxLw9fI4rdcnXE+8uRf9ZK2L/QMTH23uK2niQFkzz2wvGXa8u2+Acx78nPveHX8AA/DBkWo+O1rLa3vLp6QZ4Wh1K7/ZlHNKr7WssZPvP5dNXIgvT918Bn++ZhHb7z4fbw/DlFwtyqUT+gu7SsiIDuR/Ni4ZTDzrFsTxx6sWsrOokce3Fp3yc35xrI6rH9vB3z8vHHXfrqIGjBpWpkcCMDs2mMqW7kmXeLr7BnhuRwnr5sey6c5VLEwM5Z71c1k7L5bfv5fLh0eqxnzs2wcqmBHux4bF8WwrrB/1Q/e3z44T6u/FDSuSWWn+A7SjsIHuvgGK6zusJnTTawrieG07Pf1GrpFyy7QUH+oLQOUYrYu7ixuZ/5vNPPDB0Qklux1Fpk+H2wrqh81xfZxbQ317D8/sOIbuqc8AABQVSURBVMGRipPv+tjdN8D9H+QR7OtJZ+8AW60MYmyppauP257dy9PbT3B4nNgs2nv6+e4z2fT2G3nypjMGr+Tl6WEgNTJgShooXDahF9a1s6+0mauWJY7a8e+KpYlcsiCO/9lScEqX2+ofMA6WO/5vV+momtf2wgb8vDxYkmSaQJoda1qUM9Yo/UhFy0n/Km89Vkd7Tz8blyfh6WF6KwwGxcPXLWZhYii/fPPIsI2TLOraethWUM+GRQmcOyuatu7+Ydug5lS28EleLbesSiXQx5N58SEE+Xqyo7CBYzVtGDXMHVE/t7DU0efGBTMv3nrSF+7N39uTED8vqqy0Lg4YNb99JweF4rGthVz/xC5qx6m17yhsIMzfi36j5qPcr8sub+6vICbYh4gAb+59+8hJSy/PbD9BWWMXf924hBA/Lz48Uj35FzgOrTU/e+0gFeaS04Gy5gk95uevH+J4bRv/e8My0qOHbyudFh1AYZ39L/3nsgn9jX3lGJRppZs1916aibeHgXvfPjLhTfxfyS7nWE07G5fPoK6thw9GjJC3FdRzRmo4Pp6mTwOzY00Jz9rEaE5lC1f+fTsXP/wFP3vtINVWPr6+d6iKMH+vwRG0ha+XB/96bhqNHb3sLGoY9bh3D1Vi1PDNJfGsSo/AoBhWdnnk0+ME+Xpy08oUADwMihUzI9hR1DBmh4uFpTf96mWJKCUrQ6eruBBfqzX01/eVk1PZyh+uWshD1yziUHkz6x/5is/za60+T0dPPwfKmrn2jCQSw/x4/7Dpd6q+vYetx+r41pJE7l43l/2lzby2t9zqc9S39/C3zwo4f040582O5sLMGD7Jq7HbJONT206wOaeGX6ybQ2yw7+D+NhZGo2Z7Qf2wwdazO0p471AVd108m9UZkaOeMy0qkNLGzkmXgSdqQgldKbVWKZWvlCpQSt1t5X4fpdTL5vt3KaVSbB3oUANGzRv7KjhnVhTRwb5Wz4kN8eUnF83iy+P1vHto7NKFRVt3Hw99nM8ZKWH81zcXkBoZwNNDWgLfP1zF8dp2zpsdNXjMslx+5MRoS1cf//r8PsL8vblpZQpv7a9kzZ+28PKer/d+6e4b4JO8GtbOjxscnQ+1ZnYU/t4evH949Ejk7QOVZMYFkx4dRKi/NwsTQ/nyuCmh51S2sDmnhltWpRLi9/Wk58q0CEoaOvkkr5YAbw+Swq23Ii5NCuOxG5dy44rkcf+fCfeVEOpH5YhBSEdPP3/cnM+SpFAuWxjHFUsTefuO1YT5e3HzP/fwizcO0z6iGSG7pIl+o2ZVegTrF8SxraCels4+3jlYyYBRc8XSBK5YkkBWchgPfHiUls7R5cu/fnKcrr4Bfrl+LgDr5sfS1t3PtsKJX2Frog6Xt3D/+3lcmBnDratTWTwjdNQI/a0DFVz/xC7WPvwFXxyrY39pE79/L5cL5kRz+zlpVp83LSqQAaOmtMF6R5qtjJvQlVIewKPAOiAT2KiUyhxx2q1Ak9Y6HfgL8AdbBzrUjsIGqlq6B1e4jeXbZ6WwICGE+97NpaD25PWrv39eSH17L/95SSYGg+LbZyWzv7SZQ+XNVDZ3cffrh1iUGDIs0SmlmBUzvNPFaNT85JWDVDZ38egNS/j1ZfP45MfnsnhGKPe+nUORuY625Wgtnb0DXLowzmo8vl4enD8nms051cPqlCUNHRwoa2bD4vjBY+dkRHKgrJmWzr7B0fktq4dvbbsyzTRq+CSvhtmxQWNemEIpxdr5cXh7uuyHN2ED1pb/P7a1kLq2Hu69NHPw09vs2CDe+cFqvn/OTF7aU8rah78YVubcXliPl4ciKzmc9Qvi6BvQfJRbzRv7KpgXH8ysGNPP4n0b5tPc2cs/vhg+dzVg1Ly1v4INixMGyxirMyIJ9PHkQyuDnZG01lbbgPeWNPLoloJRx5/beQI/Lw/+dNUilFIsSQqltLFz2EW1Pz1aS5i/FwNGzbef2s31/28XMcG+/PmaRWP+XqWZr+xl7zr6RH5rlwMFWusirXUv8BKwYcQ5G4BnzF+/Blyg7Ph5/fV95QT5enLhOC11HgbF/VcsoLtvgLUPf8Hv380dNYFZ3dLNgx8e5Ykvi/nWkoTBBRZXLkvE39uDp74q5kcvH2DAqPnrdUvwGjGanh0bzLGar5fLP/ZFIZ/k1fDL9XNZlhwOQFKEP49sXIKPp4F73jSVgN49XEVEgDdnpoaPGf8lC+Jo7Ogd7CUHeGlPGUrBZYuGJPRZURg1PPFVEZtzarh19fDROcCsmEAiArzRGuaMUW4RwiIuxI/mzr7BC11Ut3Tz+BdFXLYonqVJwy9s4uvlwS/Wz+XV759FS2cfv38vb/C+nYUNLJ4Rip+3B4sSQ0gI9eMfXxRxuKJlWLk0Mz6YM1Mj2JI/fLIzt7KVtp5+zpn1dRnDx9M02Pkot3rcSdnfvpPL+X/+fFhPfXVLN997di9/3Jw/LMEOGDWf5tWyZk40If6m3x/L1Zwso/T+ASNfHa/nG3Nj2Pwf5/CztbNJDPPjf29YSqi/95hxzIwKALB7HX0iCT0BKBtyu9x8zOo5Wut+oAWIGHEOSqnblFLZSqnsurrJzVK3dffxwZEqLlsUP6GWuvkJIWy5aw1XLUvkyW3FnPPgFq5+bDu3PZvNd5/JZvUfPuOxrYWcPyeaey6ZO/i4YF8vrlyayFsHKtld3Mh9G+aTEhkw6vlnxwTS1NlHXVsPz2w/wYMf5nPpwji+sypl2HnRQb7cvW4OO4oaeH5nCZ/l1bJuQazVcovFmtnR+Hl58J657lhQ286TXxZz+aL4wdYygEUzQgny8eR/PisgyNeT71i58IRSihXmWv1Y9XMhLAY7XcwTo5sOVtDTb+QnF84a8zFZKeHcviaNT/Jq2HOikdbuPg5XtHCW+dOhUor1C2IpqG3Hw6C4fMinTDCVBfOqWmns+LoTZlexaQ5pxczh6WTd/FiaOvuGDXZGqmrp4v92lVDS0MktT5vWpvQPGPnhS/vp6h1AKXhnSG/4gbImGjp6hw0UFySG4GFQgwn9YHkzLV19nDs7Ch9PD/5tTTof//hcFiae/DJ+AT6exIX4UjhOpeB0TSShWxtpj/wMM5Fz0Fo/rrXO0lpnRUVFWXnI+D44XE13n3HccstQkYE+PHDlQjbdsZrzZkfjaTBQ2thJfk0rN61MYetPz+Oxf1lGZKDPsMfdtDIZg4LLF8VzxVLrk6+WidF73z7CrzflcFFmDH++ZpHVCcWNZySRlRzGrzbl0NU3wCUL4kedM5Sf9/Cyyz1vHsbP24N7Lx1e8fLyMLAy3fQDb210brHK/IuVKQldjCM+xNKLbiq7fHikmvkJwVYHNUPdsiqV6CAfHvjgKLuKGjFqOGtIMr5koeln/uyMSKKDhs9/WdqBdwxZBLezqIHUyABiRsyVnTs7Cl8vA+8eHnt+7MkvizFq+P0353O0uo07X9jHnz8+xu7iRv7rW/NZnhLOOwcrBz9df5Rbg5eHYs2QeTJ/b09mxwQNTox+nl+HQcHZ6aeev9KiAik8ha67yZjIFnrlwNCG5ERg5JInyznlSilPIAQY+0/naUgI8+OarESWJp36hW0XJIbwl2sXT/j89OggPv3JGhLD/Mbs+JgVY6qNbc6p4ZKFcTx87eJRZRkLg0Hx31cs4JJHviTEz7SCczzrF8Tx3uEqfv76YXYVN/LAFQtG/eEB+NaSRIrrO6yOzi2uXJZAoK/npP7fiellcHFRczc1rd3sK23mrovGHp1b+Hl78B8XzuIXbxzmT5vz8fE0DLb5AixKDOHW1alcYmXuaGFiCAHeHmwvrOeShXEMGDW7ixtZv2D0uf7enly2MJ6X95Rx2cJ4zhrRKdbS2ceLu0u5bGEcN65IRim4580jbMmv45qsRK5Ymkhn7wD/+dYR8qrayIwP5uPcGlbMjBi1gnpxUijvHKjEaNRsPVbH0qSwwZLMqUiLCuCNfRVore3WQTaREfoeIEMplaqU8gauAzaNOGcTcJP566uAz/REewVP0ar0SB68yvoI2B5SIwPGTNAAEYE+LEwM4dqsGfz1JMncYlZMEH+6ehG/2zBvQheMOG+OaSTy+r5ylqeEj7nYZ+38WD76j3PHHJ2DqfZ4+aJ4aUcU44oJ9kUp0+IiS+/42vkTu2LV1csSmRkVQH5NG8uSw4aVRpVS3Htp5qg6PJg+aZ45M2Jwm4q8qlZau/tHlVssfnVZJikR/tz5wr5RE7jP7TxBR+8At5m7Tm44M5mfXDiLszMi+e3l8wHTYMnDoHjnUCWFde0U1XVYnZdbMiOUtp5+dp9o5FB5C+fOmlx1IS06kLaefuraesY/eZLGTejmmvidwGYgD3hFa52jlLpPKXW5+bQngQilVAHwY2BUa6M7e/uOVfzhqoUnrYcPtWFxAuusjDqs8ff25Pw50Xh5KP77ivljzqILYUvengYiA32oau5m85FqZkYFkB5tfTHaSJ4eBn528WyAUWssxrMyLYLi+g4qm7vYZa6PnznT+ifZIF8v/vEvWXT3DXD78/sGe7y7+wZ4evsJzp0VReaQxXE/uCCD5249Ez9v0x+Y8ABvVqdH8s7BSj7KMf3R+sZcKwnd/AnjkU+PA6a5rcmwdLoU2LHTZUJXLdBavw+8P+LYr4Z83Q1cbdvQXIe9R7y/uWwet52TNuFfKCFsIT7El9yqVnKrWvn+OTNP6bEXz4vlkY1LTnk0a2mv3V7YwM6iBpIj/Ac3C7MmPTqQP1+ziNuf38dVf9/BzKgA2rv7qW/v5fZzrfeED3XZonjuevUgT35VxPyE4GHNBhYzIwMJ8vVke2EDkYHek15B/XXrYsfg67Q1aTZ2AdHBvoPtU0JMlfhQPw5XtDBg1BMut1gopbh8UfxJS4DWzIkNIjzAm6+O17G7uJEVqeOP8NfOj+PXl2WilKm9MLukiTWzo1gxxsh+qIvmxeDtYaC+vZcL51p/jQaDGvz9OycjatKfkmOCfQjw9rBrp4tcV0wIYZVlZBwf4suChJAp+Z4Gg+KsmRG8f6Sa3n7jmOWWkb6zKvWkDQFjCfb1Ys3sKD7KrTnpupbFM0L58ng9586eXP0cTH/k0qID7bq4SEboQgirLL3oF8+PndKJ9JXpEYP7tJw5xoSoLf3g/Az+bU3aqGvsDnXJwjjOmhkx6fq5RVpUIEV2XFwkI3QhhFUpEaaec2ttg/ZkqS/PCPcjwUpN29YWJIawIPHkn0DmxAbz4m0rTvt7pUUF8Ob+Cjp7+/H3tn36lRG6EMKq8+dE8/YdqzgjZWJlD1tJifAnPTqQ809zNOyMLBOj9hqlywhdCGGVwaCGXTx6qiil2HTnqnHXdLiitOivN+mab4d5CUnoQginY49yhDNIjvDngjnRJ93I63S45/81IYRwQj6eHjx58xl2e373+0wjhBDTlCR0IYRwE5LQhRDCTUhCF0IINyEJXQgh3IQkdCGEcBOS0IUQwk1IQhdCCDeh7HSluPG/sVJ1QMkkHx4J1NswHFcgr3l6kNc8PZzOa07WWlvdx9dhCf10KKWytdZZjo5jKslrnh7kNU8P9nrNUnIRQgg3IQldCCHchKsm9McdHYADyGueHuQ1Tw92ec0uWUMXQggxmquO0IUQQowgCV0IIdyEyyV0pdRapVS+UqpAKXW3o+OxB6XUDKXUFqVUnlIqRyn17+bj4Uqpj5VSx83/DXN0rLaklPJQSu1XSr1rvp2qlNplfr0vK6Xsc5kXB1FKhSqlXlNKHTW/12dNg/f4P8w/00eUUi8qpXzd7X1WSj2llKpVSh0Zcszq+6pMHjHns0NKqaWn871dKqErpTyAR4F1QCawUSmV6dio7KIf+InWei6wArjD/DrvBj7VWmcAn5pvu5N/B/KG3P4D8Bfz620CbnVIVPbzV+BDrfUcYBGm1+6277FSKgH4IZCltZ4PeADX4X7v89PA2hHHxnpf1wEZ5n+3AX8/nW/sUgkdWA4UaK2LtNa9wEvABgfHZHNa6yqt9T7z122YftETML3WZ8ynPQN80zER2p5SKhG4BHjCfFsB5wOvmU9xt9cbDJwDPAmgte7VWjfjxu+xmSfgp5TyBPyBKtzsfdZafwE0jjg81vu6AXhWm+wEQpVScZP93q6W0BOAsiG3y83H3JZSKgVYAuwCYrTWVWBK+kC04yKzuYeBnwFG8+0IoFlr3W++7W7v9UygDvinucz0hFIqADd+j7XWFcCfgFJMibwF2It7v88WY72vNs1prpbQlZVjbtt3qZQKBF4HfqS1bnV0PPailLoUqNVa7x162Mqp7vReewJLgb9rrZcAHbhRecUac914A5AKxAMBmEoOI7nT+zwem/6cu1pCLwdmDLmdCFQ6KBa7Ukp5YUrm/6e1fsN8uMbyccz831pHxWdjq4DLlVInMJXRzsc0Yg81fzQH93uvy4FyrfUu8+3XMCV4d32PAb4BFGut67TWfcAbwErc+322GOt9tWlOc7WEvgfIMM+Ke2OaUNnk4Jhszlw/fhLI01o/NOSuTcBN5q9vAt6e6tjsQWv9C611otY6BdN7+pnW+gZgC3CV+TS3eb0AWutqoEwpNdt86AIgFzd9j81KgRVKKX/zz7jlNbvt+zzEWO/rJuDb5m6XFUCLpTQzKVprl/oHrAeOAYXAPY6Ox06vcTWmj12HgAPmf+sx1ZU/BY6b/xvu6Fjt8NrXAO+av54J7AYKgFcBH0fHZ+PXuhjINr/PbwFh7v4eA78FjgJHgOcAH3d7n4EXMc0R9GEagd861vuKqeTyqDmfHcbUATTp7y1L/4UQwk24WslFCCHEGCShCyGEm5CELoQQbkISuhBCuAlJ6EII4SYkoQshhJuQhC6EEG7i/wOBwcKSnVdXSQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x0 = 0*xn\n", "res = opt.minimize(func, x0, method='trust-constr', tol=1e-12, options={'gtol': 1e-6, 'disp': True})\n", "plt.plot(res.x)\n", "\n", "print('\\nFinal value of the cost function:')\n", "print(res.fun)\n", "\n", "print('\\nError with respect to original noiseless image:')\n", "print(la.norm(res.x-x)**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Hard constraints\n", "\n", "Further, suppose we know a priori which pixels within the image correspond to tissues (a total of $N_T$ pixels where the image intensity is different from zero), and which pixels correspond to air (a total of $N_A$ pixels where the image intensity is zero). The total number of pixels in the image is $N = N_T + N_A$. This constraint can be expressed as $\\mathbf{Cx} = \\mathbf{0}$, where $\\mathbf{C}$ is a subsampling matrix that selects all the pixels in the air regions (and ignores the rest), and $\\mathbf{0}$ is a zero-vector of the appropriate length $N_A$. \n", "\n", "\n", "Putting all these elements together, we have the following constrained formulation:\n", "$$\\hat{\\mathbf{x}} = \\arg \\min_{\\mathbf{x}} \\| \\mathbf{Ax} - \\mathbf{d} \\|_2^2 $$\n", "$$\\hbox{such that } \\mathbf{Cx} = \\mathbf{0}$$" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import LinearConstraint\n", "linear_constraint = LinearConstraint([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])\n", "C = np.identity(N)\n", "C1 = C[0:20,:]\n", "C2 = C[N-20:N,:]\n", "C = np.concatenate((C1,C2),axis=0)\n" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "linear_constraint = LinearConstraint(C, np.zeros(40), np.zeros(40))" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "`gtol` termination condition is satisfied.\n", "Number of iterations: 4, function evaluations: 408, CG iterations: 3, optimality: 1.30e-08, constraint violation: 0.00e+00, execution time: 0.2 s.\n", "\n", "Final value of the cost function:\n", "0.04002062856627106\n", "\n", "Error with respect to original noiseless image:\n", "0.06380854701416211\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXRbd5n/8fcjyfu+x7vj2FmbtknctOlOaUtSIClbacrSQofCsBQow9CBGX5QYObADMMyFEqhUOBASxcohQa6N11oQpykSRo7i+MktuPdjvfd+v7+kOw4jh07tqQrXT2vc3JiSVfSo1z7k6+f+73fK8YYlFJKhT6H1QUopZTyDQ10pZSyCQ10pZSyCQ10pZSyCQ10pZSyCZdVb5yenm6KioqsenullApJO3fubDXGZEz1mGWBXlRURHl5uVVvr5RSIUlEjk/3mLZclFLKJjTQlVLKJjTQlVLKJjTQlVLKJjTQlVLKJjTQlVLKJjTQlVLKJjTQlZpG39AIv3n9GAPDo1aXotSsWHZikVLBbHjUzad+u4sXD7aQGBPBpgtzrS5JqRnpCF2pSYwx3P34Pl482ILTIZQfO2l1SUrNio7QlZrkO08f5PFddXz+2sXsONZO+XENdBUaZhyhi8gvRKRZRN6c5nERkR+KSJWI7BWR1b4vUynfevVwKxt+8ArVLT2n3f/H3XX85KUj3HJxAXe+tYSyohQONHbRNTBsUaVKzd5sWi4PAuvP8vgGoNT75w7gJ/MvSyn/GR5189Un36SyoYtP/W73+EHPmrY+/uOJ/VxUlMI9G1cgIlxUlIoxsEtH6SoEzBjoxpiXgfazbLIJ+LXx2AYki0i2rwpUytce3lFLdUsvt64rpLKhi2/8pYKRUTef/f1uROB7778Ql9Pzo3FhfjJOh7BTA12FAF/00HOB2gm367z3NUzeUETuwDOKp6CgwAdvrdS56R4Y5vvPHmLtwlS+tnEFURFO7n+5muNtfeyu6eCHm1eRlxI7vn1clIvl2YnsOHa2MY1SwcEXs1xkivvMVBsaY+43xpQZY8oyMqZcn10pv/rp1mraeof4yg3LEBG++LYlrCpI5tWqVt69OpeNF+Sc8ZyyohTeqO1gaMRtQcVKzZ4vAr0OyJ9wOw+o98HrKuVTDZ39/PzVajZekMMF+ckARDgd/PgDq7nzraXcs+m8KZ9XVpjKwLCb/fWdgSxXqXPmi0B/Eviwd7bLJUCnMeaMdotSVvvrvkYGht3cdd3i0+7PTorhrusWEx81dQeyrCgFQPvoKujNZtriQ8DrwBIRqROR20XkEyLyCe8mW4BqoAr4GfBJv1Wr1Dwcb+slIcpFYVrszBtPkJUYTUFqrPbRVdCb8aCoMWbzDI8b4FM+q0gpP6lp7yM/NRaRqQ77nF1ZYQpbD7VgjJnT85UKBD31X4WNmvY+ClLPbXQ+pqwolbbeIY619fm4KqV8RwNdhQW321B7sv+c2y1j1i1KA+APu+p8WZZSPqWBrsJCU/cAQyNu8uc4Ql+YHscNKxfwi1eP0t475OPqlPINDXQVFmq8rZK5tlwAPn/tYvqGR/np1iO+Kkspn9JAV2Ghpn3+gV6alcCNF+byq9eP0dw14KPKlPIdDXQVFmrb+3AI5CTHzOt1PndtKcOjhh+/pKN0FXw00FVYqGnvIyc5hkjX/L7lC9PiuKksj99tr+H1I214Zu0qFRw00FVYOD6PKYuTfeaaUhJjXGz+2Tbe+aNXeWxnHaNuDXZlPQ10FRZqfRjoOckxvPyvb+Fb7zqPwWE3//LoHr71VKVPXlup+dBAV7bXOzhCa8/QnKcsTiU20sUHLi7kmc9fyYfXFfKL146y9VCLz15fqbnQQFe2V3vSM8NlricVnY2I8OUblrE4K55/eXQPbT2DPn8PpWZLA13Z3nEfzEE/m+gIJz+4eRWdfcN86fF9eqBUWUYDXdlerQ/moM9kWXYiX9qwlOcqm/jrm41+ex+lzkYDXdleTXsfCdEukmIi/Po+H7m0iJgIJ+XHdN10ZQ0NdGV7Ne19FKbNbdncc+FwCMUZcVS39vj1fZSajga6sr35LJt7rooz4jnSooGurKGBrmxt1G2oa+/36ZTFsylOj6PuZD8Dw6MBeT+lJtJAV7bW1DXA0Kg7YCP0RZnxGHNqZo1SgaSBrmzNF6ssnovi9DgAbbsoS2igK1sYGB6laYolbY+39QJQmBoXkDqKMzzvU62Briygga5s4b6tR7j6v186bWQ8NOLmgVePkpscQ05ydEDqiI10kZMUTXVLb0DeT6mJNNCVLdS09dE/PMpdv3+D4VE3APe/fIRDTT1848YVuJyB+1bXmS7KKhroyhZaegaJjXSyp66TH71QRXVLDz98oYq3n5/NNUuzAlpLcUYc1S29ugSACjiX1QUo5QutPUOsK04jKTaCH71Yxd/ebCTK5eD/vXN5wGspTo+je3CElp5BMhMC0+pRCnSErmyirWeQ9PgovrZxBQsSoznY1M2Xb1hmSaAuyowH4Eiz9tFVYGmgq5DndhvaeodIi48kMTqCn35oDV982xLeX5ZvST3FGZ5A1yUAVKBpy0WFvI7+YUbdhvT4KADOy03ivNwky+rJTowmOsKhI3QVcLMaoYvIehE5KCJVInL3FI8XiMiLIrJbRPaKyA2+L1WpqbV6LyqRnhBlcSUeDodQnB6vI3QVcDMGuog4gXuBDcByYLOITD7S9O/AI8aYVcDNwI99XahS02nt9gZ6fKTFlZwyNtNFqUCazQh9LVBljKk2xgwBDwObJm1jgETv10lAve9KVOrsWrwj9Iz44Bihg6ePXnuyTxfpUgE1m0DPBWon3K7z3jfR14APikgdsAX4jE+qU2oWWnuGAMZ76MFgUUacLtKlAm42gT7VVQEmnzGxGXjQGJMH3AD8RkTOeG0RuUNEykWkvKVFr5CufKOtZxCXQ/x+RaJzsWhspoueMaoCaDaBXgdMnP+Vx5ktlduBRwCMMa8D0UD65BcyxtxvjCkzxpRlZGTMrWKlJmntGSQ1LhKHw79XJDoXC9PjEIE9dZ1Wl6LCyGwCfQdQKiILRSQSz0HPJydtUwO8FUBEluEJdB2Cq4Bo7RkKqnYLQFyUi2uWZPLYzloGR7SPrgJjxkA3xowAnwaeBirxzGbZLyL3iMhG72ZfAD4mInuAh4DbjC5koQKktWcwaKYsTnTrpUW09gzx1N4Gq0tRYWJWJxYZY7bgOdg58b6vTvi6ArjMt6UpNTut3YOUeE+3DyZXlKazKCOOB/9+jHetyvX7RaqV0lP/VUgzxtDaMxRUUxbHiAi3XVrE3rpOdtd2WF2OCgMa6CqkdQ2MMDTqDroe+ph3r84jIcrFg68ds7oUFQY00FVIaxs/7T94zhKdKC7KxU0X5bNlX8OUl8hTypc00FVIC8aTiib78LpCRo3hsZ11VpeibE4DXYW0sYW50uKCN9AL0+IoTI2lor7L6lKUzWmgq5DWGuQtlzElmfFUNetZo8q/NNBVSGvtHkQEUmODO9AXZcRztLWXEe8FrJXyBw10FdJaeoZIjY3E5Qzub+VFmfEMjbqpO9lvdSnKxoL7p0CpGYxdSzTYjZ34pG0X5U8a6CqkeU77D+52C5xafbFKV19UfqSBrkJaMC7MNZWkmAgyEqJ0hK78SgNdhbTWnsGgnrI4UUmGznRR/qWBrkJW39AIfUOjIdFyAU8f/UhLD7oQqfIXDXQVslq7g/8s0YlKMuPpHhihxXtRa6V8TQNdhaxgvDj02YwfGNW2i/ITDXQVUjr7hqnxXnh5fGGuEAn08amLOtNF+cmsLnChVLD4xlMV/GFXHZ++ppSMeE/vPFR66FmJUcRHuXSErvxGA12FlDdqO4iNdPHD5w8THeH5BTNUZrmICIu8B0aV8gdtuaiQMTA8SnVLDx+5rIgf3bKKKJeTzIQoIl2h822sUxeVP+kIXYWMQ03duA0sy07khpXZXFKcRvfAiNVlnZNFmXE8vquOroFhEqMjrC5H2UzoDG1U2Kts8Kwnviw7EfAcDF2YHmdlSeesxDvT5YiO0pUfaKCrkFHZ0E1spJPC1FirS5kzXaRL+ZMGugoZFQ1dLFmQgMMhVpcyZwWpsUQ6HVQ06NWLlO9poKuQYIyhsqFrvN0SqlxOB9cszeSxnZ4+ulK+pIGuQsKJjn66B0ZCPtABPn1NCd0DI/z678esLkXZjAa6CgmVDd0ALM9OsLiS+TsvN4lrlmbywKtH6R0MrVk6KrhpoKuQMDbDZcmC0B+hg2eUfrJvmN9uP251KcpGNNBVSKhs6KIwLZb4KHucOrG6IIXLS9K5/+WjDAyPWl2OsolZBbqIrBeRgyJSJSJ3T7PNTSJSISL7ReR3vi1ThbvKhi6W26B/PtFnrimhtWeQh/9RY3UpyiZmDHQRcQL3AhuA5cBmEVk+aZtS4N+Ay4wxK4DP+aFWFaZ6B0c43t5niwOiE11cnEZpZjxbD7VYXYqyidmM0NcCVcaYamPMEPAwsGnSNh8D7jXGnAQwxjT7tkwVzg40dmMMtgt08BwgPdDYbXUZyiZmE+i5QO2E23Xe+yZaDCwWkddEZJuIrJ/qhUTkDhEpF5HylhYdlajZOXXKf+jPcJls6YIEGjoH6OgbsroUZQOzCfSpTsubfFFEF1AKXA1sBn4uIslnPMmY+40xZcaYsoyMjHOtVYWpyoYuEqNd5CbHWF2Kzy31/taho3TlC7MJ9Dogf8LtPKB+im3+ZIwZNsYcBQ7iCXil5sUYw6tVrawuTEEkdE/5n86yBZ7fOg7oUgDKB2YT6DuAUhFZKCKRwM3Ak5O2eQJ4C4CIpONpwVT7slAVng419XC8rY/rly+wuhS/yEiIIjUuUkfoyidmDHRjzAjwaeBpoBJ4xBizX0TuEZGN3s2eBtpEpAJ4EfiiMabNX0Wr8PHM/kYArl2WaXEl/iEiLMlKoFIDXfnArM7SMMZsAbZMuu+rE742wF3eP0r5zDMVTawqSCYzMdrqUvxmaXYCD/+jllG3wRnCK0kq6+mZoipo1Xf0s+9Ep23bLWOWLUikf3iUmvY+q0tRIU4DXQWtZyuaALh+RZbFlfjX0mw9MKp8QwNdBa1nKhpZlBHHIu9l2+yqNDMBh+jURTV/GugqKHX2DbO9up3rV9i73QIQE+mkKD2OA406Qlfzo4GugtKLB5sZcRuuX27vdsuYZQsSdYSu5k0DXQWlZyubyEyI4oK8M044tqWlCxI43tanF7xQ86KBroLS3roOLipKDekLQp+LsSUADjbpKF3NnQa6CjrdA8PUtvfbcjGu6SwdXwJAA13NnQa6CjpjvWQ7Lpc7nbyUGOKjXOMrSyo1FxroKuiMhdrynPAJdBGhNCueQ9pyUfOgga6CTmVDF8mxESyw8en+UynJiOdIS4/VZagQpoGugk5FfRfLFiTacrncsynJjKe1Z0gvdqHmTANdBZVRt+FgU3dY9c/HlGR6zoitatZRupobDXQVVI629jIw7A6r/vkYDXQ1XxroKqjY+fqhM8lLiSXS5dBAV3Omga6CSmVDFy6HjI9Ww4nTIRSnx1GlB0bVHGmgq6BS0dBFSWY8US6n1aVYoiQzXkfoas400FVQqWzoYnkYHhAdU5IZz4mOfvqHRq0uRYUgDXQVNNp7h2jqGgzLGS5jSjLjMQadj67mRANdBY1TB0TDN9BLMz0HgzXQ1VxooKugUVEfvjNcxhSlx+IQnbqo5kYDXQWNyoYushKjSIuPsroUy0S5nBSmxWmgqznRQFdBo6KhK6zbLWMWZehMFzU3GugqKAwMj3K4uYcVYXiG6GQlmfEca+tlZNRtdSkqxGigq6BwuKmHUbdhRU6S1aVYriQznuFRw/H2PqtLUSFGA10Fhf31nQA6Quf0NV16B0f49t8O8MiOWourUqHAZXUBSgHsr+8iIcpFfkqs1aVYblFGHAB/euME33yqgtr2fpZkJXDTRfkWV6aC3axG6CKyXkQOikiViNx9lu3eKyJGRMp8V6IKB/vrO1mWnRg2F4U+m4Roz8U9tuxrxCnC5SXp1LT3YYyxujQV5GYMdBFxAvcCG4DlwGYRWT7FdgnAncB2Xxep7G3Ubahs6A7LJXOn8/Grivnk1Yv462ev5NplmfQPj9LSM2h1WSrIzablshaoMsZUA4jIw8AmoGLSdt8AvgP8i08rVLZ3tLWX/uFR7Z9P8JHLFo5/XZjmacHUtveRmRBel+VT52Y2LZdcYOIRmTrvfeNEZBWQb4z5y9leSETuEJFyESlvaWk552KVPZ06IKozXKaSn+o5rnC8TWe9qLObTaBP1dQcb+aJiAP4HvCFmV7IGHO/MabMGFOWkZEx+yqVrVXUdxHpdFCaFX5roM9GXkoMIlCj0xjVDGYT6HXAxMPreUD9hNsJwHnASyJyDLgEeFIPjKrp1J3s42TvqQshVzR0sXhBPBFOnUU7legIJwsSo6nREbqawWx+gnYApSKyUEQigZuBJ8ceNMZ0GmPSjTFFxpgiYBuw0RhT7peKVcj76IM7uOXn2xkedWOMYX99eK+BPhsFqbE6QlczmjHQjTEjwKeBp4FK4BFjzH4RuUdENvq7QGUvo25DdUsvlQ1d/PyVozR2DdDeO6T98xkUpMbqmaNqRrM6scgYswXYMum+r06z7dXzL0vZVUNnPyNuQ0K0i+8/d4gIp+cQjc5wObvCtFhaugfpHxolJjI8L8+nZqZNSxVQY22Dr29cQaTLwX/99QAi4X1Ri9kYm+mibRd1NhroKqDq2vsBuKgolS/fsIxRt2FhWhxxUboKxdmMzUXXQFdnoz9FKqBq2vtwOoTspGhuviiflw+1sChDpyvOpGB8LnqvxZWoYKaBrgKq9mQfOcnRuLxTFH/ywTUWVxQaUmIjSIhyUasjdHUW2nJRAVXT3qcrKs6BiJCvM13UDDTQVUDVtveNtw/UuSlM07no6uw00FXA9A2N0NozND5jQ52bgtRY6tr7GXXrMrpqahroKmBqvTNcNNDnpiAtlqFRN01dA1aXooKUBroKmLF2gbZc5qYw1TN1UVddVNPRQFcBU6uBPi9j/24600VNRwNdBUxNex9xkU5SYiOsLiUk5SRH43QIx9t1Lrqams5DVwFT295HfmosInrd0LlwOR3kJsdQ2dDNoaZu3MaQnRRDUoz+B6k8NNBVwNSe7Bs/hV3NzaKMOF440MwLB5oBWJ6dyJbPXmFxVSpYaKCrgDDGUNPexxWleqWq+fjmu1ay6/hJHCJsPdTMI+V1tPYMkh4fZXVpKghooKuAaOkZZGDYrQdE5yk3OYbc5BgAFiRF80h5HTuOtrNhZbbFlalgoAdFVUCcmoMeY3El9rEyN4noCAfbj7ZbXYoKEhroKiB0yqLvRbocrC5IYccxDXTloYGuAmLspKI8XZjLp9YuTKWioYuugWGrS1FBQANdBURtex+ZCVFER+jl03xp7cJUjIGdx05aXYoKAhroKiBqdJVFv1iVn0KEU7SPrgANdBUAxhiOtvZqoPtBTKST8/OS+cfRNqtLUUFAA135XW17P83dg1xYkGx1Kba0dmEqe+s66R8atboUZTENdOV327yjx4sXpllciT2tXZjKiNuwu0b76OFOA1353fbqdlLjIinN1ItB+8OawhQcgvbRlQa68r9t1W2sLUrF4dBFufwhMTqC5TmJbNc+etjTQFd+VXeyjxMd/VxSnGp1Kba2rjiNXcc76BkcsboUZSENdOVX26s9bYCLi7V/7k/XLstiaNTN1oMtVpeiLKSBrvxq+9E2kmMjWJKVYHUptramMIWU2AierWi0uhRloVkFuoisF5GDIlIlIndP8fhdIlIhIntF5HkRKfR9qSoUbatu1/55ALicDq5ZmsULB5oZHnVbXY6yyIyBLiJO4F5gA7Ac2CwiyydtthsoM8acDzwGfMfXharQ09DZT017n7ZbAuS65Vl0DYywQ2e7hK3ZjNDXAlXGmGpjzBDwMLBp4gbGmBeNMWNXrt0G5Pm2TBWKxvrnekA0MK5cnE6Uy8EzFU1Wl6IsMptAzwVqJ9yu8943nduBv071gIjcISLlIlLe0qIHb+xuW3UbidEuli5ItLqUsBAb6eLyknSerWjCGGN1OcoCswn0qZqfU363iMgHgTLgv6d63BhzvzGmzBhTlpGhlyKzq76hER4tr+XZiibWLkzFqf3zgLlueRYnOvqpbOi2uhRlgdlcgq4OyJ9wOw+on7yRiFwLfAW4yhgz6JvyVCjp6Bvi+88d5rGddfQMjlCcHsc/X11idVlh5a3LshDZx7MVTSzP0d+Mws1sAn0HUCoiC4ETwM3ALRM3EJFVwE+B9caYZp9XqYKaMYan9jXwtSf3c7JvmE0X5rB5bQFlhSmI6Og8kDISoliVn8zT+xu5860l+u8fZmYMdGPMiIh8GngacAK/MMbsF5F7gHJjzJN4WizxwKPeb6AaY8xGP9atgoQxhjsffoM/76lnZW4Sv/roWlbkJFldVlh79+o8/v2JN9l6qIWrl2RaXY4KILHq4ElZWZkpLy+35L2V7xxr7eXq/3mJ2y4t4t/fvgyXU89Vs9rQiJtr/3cr8VEu/vKZy/UcAJsRkZ3GmLKpHtOfPjUvu7xLtt68Nl/DPEhEuhx84frFVDR08Zd9DVaXowJIfwLVvOw8fpL4KBelmXpqfzB55/k5LF2QwHefOcjQiJ45Gi400NW87Krp4ML8ZJ2aGGQcDuFL65dyvK2P35fXzvwEZQsa6GrOegZHONjYxWq9tFxQunpJBmuLUvnes4eoatZ56eFAA13N2d7aDtwGVhemWF2KmoKI8J/vXolDhPf/dBuVDV1Wl6T8TANdzdnO454DoqvyNdCDVUlmPI98/BIinA42/2wb++o6rS5J+ZEGupqzXTUnKcmMJyk2wupS1FkUZ8TzyMfXERfp4gM/38bhJm2/2JUGupoTt9uwu7aDNQU6Og8FBWmxPHzHJURFOLntlzto7h6wuiTlBxroiuNtvbznJ3/nREf/rJ9T3dpLR98wqwv1gGioyE+N5Re3XkR77xC3P1hO35Bef9RuNNAVP3juMDuPn+S1qtZZP2fshKLVOkIPKSvzkvi/zavYX9/JnQ+9ocvs2owGepg71trLE2+cAOBQ4+x7q7trTpIY7WJRRry/SlN+cu3yLL74tqU8V9nEG7UdVpejfGg2qy0qG/vxS1W4nA4yE6I41Nwz7Xa9gyN89uHdJMZEUFaYyrbqdlYVpOg6ISHqA5cU8L3nDvHE7hOs0t+ybENH6GGstr2PP+w6wS1rC7ioKPWssx8e21nHc5XNvHCgmS//cR9HW3sp0/nnISsxOoLrlmXx570Np11U+s0Tndz+4A76h0YtrE7NlQZ6GPvxS0dwiPDxq4pZnJVAQ+cAnf3DZ2w36jb88rWjXJifzO7/uI4XvnAV996ymlsvKwp80cpnblyVS3vvEK8c9lwO0hjD1/+8n+cPNHNQpzaGJA30MNXYOcBjO2u56aI8spNiWJzl6YVPdYr485VNHGvr45+uWIiIUJwRz9vPzyYxWuefh7KrFmeQHBvBE7s9FyB78WAzO455DnY3nMOMJxU8NNDD1GtVrQyPGj50SREAi7M8qyUeajqzj/7Aq0fJTY5h/YoFgSxR+Vmky8E7zs/mmYpGugaG+c7fDpKVGAVAfafOUw9FGuhh6kBjF1EuB4sy4gDITY4hNtLJoUm/ar95opPtR9u59dJCXe/chm68MJeBYTeffWg3Bxq7+fINy4hyOWjs1BF6KNKf0DBV2dDNkgUJ4yHtcAilmfFnBPoDrx4lLtLJ+y8qsKJM5WdrClPIS4nhxYMtLM9O5J3n55CTHKMj9BClgR6mDjR2sXTB6RelKM1KOK3l0tw1wJ/31PO+snySYrRfbkciwrtW5QLwxfVLcDiE7KRo7aGHKA10mzvU1M1nHtp92jS0lu5BWnuGWLog8bRtl2Ql0NI9yMneIQD+sPsEI27Dh9YVBrRmFVifuGoRv7itjKsXZwCQnRRDg47QQ5IGeghr7hrg2Yqms27z7b8e4M976tlW3TZ+34FGz7rYS7Mnj9A9M10ONXVjjOHR8lrWFKbo2aA2Fxfl4pqlWYh4ThLLSY6mqWuAkVG9dF2o0UAPYd95+iAf+3U5zV1Tj6Yq6rt4/kAzAK9PDPQGT5982aQR+vhMl+Ye3qjt4EhLL+9bk+eP0lUQy06KwW2guXvQ6lLUOdJAD1EDw6M8/WYjAK9Os6jWvS9VER/lYkVOIq8fORXolY1dLEiMJiUu8rTts5OiSYhycbipm0d31hEd4eDt52f770OooJSdHA1Ag49muoyMuvncw7vHF3RT/qOBHqJeOthM9+AIIvDq4TMDvbqlhy37GvjQukKuW57F/vrO8bNADzR0n9FuAc8BstKsePbUdfLnPfVsOC+bBD15KOxkJ3kCvb7DN330vSc6eeKNeh545ehp94+Mutn4o1f51d+P+eR9lAZ6yHpyTz3p8ZHccF42Lx9uPWMZ1J+8dIRIp4OPXraQS4rTcBv4x9F2hkfdVDX3nHFAdMzirAT21HbQPTCi7ZYwlZ0UA0w/Qj/c1M2GH7zC4zvrZvV6Y8dvXjjQfNrB+deOtLG3rpMfPn+YgWFdO8YXNNBDUPfAMM9XNvP2ldlctSSD1p5BDkxY+vZERz9/3H2CzWsLyEiIYlVBMlEuB38/0kp1Sy9Do26WTTFCh1N99LyUGC4pTgvI51HBJTHaRVykc9qZLt98qpLKhi6+8OgevvLHfQyOnD2MXz/SRpTLQf/wKFsPNY/f/8TuE0Q6HbT1DvFoea1PP0O40kAPQc/sb2JwxM3GC3O4ojQdOL3tct9LRxCBO64sBiDK5WRNYQqvH2k7NcPlLCN0gPesztOlccOUiJCdHEPDFC2Xlw+1sPVQC3dvWMonrlrEb7fXcNN9r1PT1jflaw2NuCk/dpL3rskjJTaCLfs8x316B0f425uNvGdNHqsKkrn/lWqdVeMDswp0EVkvIgdFpEpE7p7i8SgR+b338e0iUuTrQtUpT+6pJzc5htUFKWQnxVCSGc/L3hXzGjr7+f2OWt67Jp+c5Jjx56wrTuNAYzd/r2oj0umg2EysytgAAAnWSURBVHvK/2RrF6Zy5zUl3HZpUSA+igpS2UnRZ7RcRt2G/9xSSUFqLB+5rIi7Nyzlpx9aQ3VrL+t/8DK/215zRutv34kO+odHuaI0nbetWMDzlU0MDI/yTEUj/cOjvHt1Lp+4ahG17f1s8R7kV3M3Y6CLiBO4F9gALAc2i8jySZvdDpw0xpQA3wO+7etClUdbzyCvVrWy8cKc8XnDV5Sm84+j7QwMj3LfS0dwG8Mnr1502vPWLfK0T5544wQlmfFETLMuS6TLwV3XLzljBowKLzlJZ57+/2h5LQcau/nS+qVEuZwAvG3FAp7+3JWsKkjmy3/cx0ce3EFn36klmMdmV128MI0NK7PpHRrllcOt/HF3PXkpMawpSOG6ZVksyojjvpeO6CXx5mk2VyxaC1QZY6oBRORhYBNQMWGbTcDXvF8/BvxIRMT4Ye/sq+uk/Hi7r182KF2Qn3zGNTv/vKeeUbdh4wU54/ddWZrBL187xlN7G3hoRy3vXZNHfmrsac87Py+ZmAgn/cOjU85wUWqi7ORoWnsGGRpxE+ly0Ds4wnefPcSawhRuWHn6qps5yTH85qMX8+vXj/GNpyr57rMHuWfTeYDn/IelCxJIiYvk0kVpJMVE8OvXj/FaVSufvLpkvK338SsX8a+P7+XbE1Z8tLN1i9KmbXvOx2wCPReYeMSiDrh4um2MMSMi0gmkAafNpxORO4A7AAoK5rbY09+PtPJffz0wp+eGmking7/cefl4X7ujb4j/e6GK1QXJp63DcnFxKhFO4at/ehO32/Cpt5Sc+VouB2VFKbxyuPWME4qUmiwnKQZjoKlrgPzUWJ7a20BL9yD33rJ6/DfDiRwO4bbLFnK4uYffba/hny4vJispip3HT7J5rednPcLp4PrlWTzqnR1z46pTg5JNq3L40YtV3Lf1SGA+oMW+eeN5lgX6VEfGJo+8Z7MNxpj7gfsBysrK5jR6v/XSIt5/Uf5cnhpSuvpHuPHHr/HFR/fw+D9fisvp4Nt/O0BH/zC/uXHlaT9UsZEu1hSmsK26nZvKzhydj1m3KI1XDrfqCF3NaOzkovqOfvJTY3l6fyO5yTFcVHT2yw7e+dZSHt9Vx/8+e5BbLi5kYNjNugmzpW5Ymc2jO+tYmZtESeap78Mol5Pn7rqKvqER/3ygIBMd4fTL684m0OuAiQmaB9RPs02diLiAJMAvfZHoCKff/jGCSXJsJF/fuILPPLSbn71ylIuKUnjoH7V87IqFLM8583/2a5dlset4x5Sj8zHvWZ3HiZP9XFSU6s/SlQ2MnVzU0DlA7+AIr1S18oGLC6YcnU+UlRjNRy5byH1bjzA06kbE0z8fc1lJOityEvno5UVnPDfS5SDSpcdu5mM2gb4DKBWRhcAJ4GbglknbPAncCrwOvBd4wR/983DzjvOzeWpvA9977hA5SdHkJEXzuWsXT7ntbZcW8Y7zc1jg/UGcSlZiNN9610p/lats5NTJRQNsPdTC0Iibt83yilWfuHIRv912nC37GlmRk0hS7KmzjSNdDp668wq/1KxmMcvFGDMCfBp4GqgEHjHG7BeRe0Rko3ezB4A0EakC7gLOmNqozp2I8I0bzyM20smxtj6+tnEFcVFT/x/scjrOGuZKnYu4KBeJ0S4aOvt5en8jKbERlBWevd0yJik2gn++2vOb4jo9OS2gZjNCxxizBdgy6b6vTvh6AHifb0tTABkJUdz3wTXsq+vker2mpwqgnOQYjrf1savmJOtXLDinSxDedmkRVc09vK/M/se7gsmsAl1Z65LiND0NXwVcdlI0rxxuZcRtZt1uGRMT6eS7N13gp8rUdPTUf6XUlLKTYxhxG2IjnVzuXWJCBTcNdKXUlHK8x2SuXpIRFjPL7EADXSk1pbGZLufablHW0UBXSk3pLUsz+dgVC7l+uQZ6qNCDokqpKaXGRfKVt09eh08FMx2hK6WUTWigK6WUTWigK6WUTWigK6WUTWigK6WUTWigK6WUTWigK6WUTWigK6WUTYhV16EQkRbg+Byfns6k65WGAf3M4UE/c3iYz2cuNMZkTPWAZYE+HyJSbowps7qOQNLPHB70M4cHf31mbbkopZRNaKArpZRNhGqg3291ARbQzxwe9DOHB7985pDsoSullDpTqI7QlVJKTaKBrpRSNhFygS4i60XkoIhUicjdVtfjDyKSLyIvikiliOwXkc96708VkWdF5LD37xSra/UlEXGKyG4R+Yv39kIR2e79vL8XkUira/QlEUkWkcdE5IB3X68Lg338ee/39Jsi8pCIRNttP4vIL0SkWUTenHDflPtVPH7ozbO9IrJ6Pu8dUoEuIk7gXmADsBzYLCJ2vKTKCPAFY8wy4BLgU97PeTfwvDGmFHjee9tOPgtUTrj9beB73s97Erjdkqr85wfA34wxS4EL8Hx22+5jEckF7gTKjDHnAU7gZuy3nx8E1k+6b7r9ugEo9f65A/jJfN44pAIdWAtUGWOqjTFDwMPAJotr8jljTIMxZpf36248P+i5eD7rr7yb/Qq40ZoKfU9E8oC3Az/33hbgGuAx7yZ2+7yJwJXAAwDGmCFjTAc23sdeLiBGRFxALNCAzfazMeZloH3S3dPt103Ar43HNiBZRLLn+t6hFui5QO2E23Xe+2xLRIqAVcB2IMsY0wCe0AcyravM574P/Cvg9t5OAzqMMSPe23bb18VAC/BLb5vp5yISh433sTHmBPA/QA2eIO8EdmLv/Txmuv3q00wLtUCXKe6z7bxLEYkHHgc+Z4zpsroefxGRdwDNxpidE++eYlM77WsXsBr4iTFmFdCLjdorU/H2jTcBC4EcIA5Py2EyO+3nmfj0+zzUAr0OyJ9wOw+ot6gWvxKRCDxh/ltjzB+8dzeN/Trm/bvZqvp87DJgo4gcw9NGuwbPiD3Z+6s52G9f1wF1xpjt3tuP4Ql4u+5jgGuBo8aYFmPMMPAH4FLsvZ/HTLdffZppoRboO4BS71HxSDwHVJ60uCaf8/aPHwAqjTH/O+GhJ4FbvV/fCvwp0LX5gzHm34wxecaYIjz79AVjzAeAF4H3ejezzecFMMY0ArUissR711uBCmy6j71qgEtEJNb7PT72mW27nyeYbr8+CXzYO9vlEqBzrDUzJ8aYkPoD3AAcAo4AX7G6Hj99xsvx/Nq1F3jD++cGPH3l54HD3r9Tra7VD5/9auAv3q+LgX8AVcCjQJTV9fn4s14IlHv38xNAit33MfB14ADwJvAbIMpu+xl4CM8xgmE8I/Dbp9uveFou93rzbB+eGUBzfm899V8ppWwi1FouSimlpqGBrpRSNqGBrpRSNqGBrpRSNqGBrpRSNqGBrpRSNqGBrpRSNvH/AU57ZHOay3TRAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "res = opt.minimize(func, x0, method='trust-constr', tol=1e-6, constraints=linear_constraint, options={'gtol': 1e-6, 'disp': True})\n", "plt.plot(res.x)\n", "\n", "print('\\nFinal value of the cost function:')\n", "print(res.fun)\n", "\n", "print('\\nError with respect to original noiseless image:')\n", "print(la.norm(res.x-x)**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Soft constraints\n", "\n", "The example above applied hard constraints to an image reconstruction problem, and forced the pixel values to be zero in the image regions known to be air. A related approach to pose the reconstruction problem is to apply a soft constraint. In this case, instead of forcing the pixel values to be zero in the air regions, we can encourage small pixel values in this regions. For example, one might pose the reconstruction as follows: \n", "\n", "$$\\hat{\\mathbf{x}} = \\arg \\min_{\\mathbf{x}} \\| \\mathbf{Ax} - \\mathbf{d} \\|_2^2 + \\lambda \\| \\mathbf{Cx}\\|_2^2 $$\n", "\n", "for some value of $\\lambda>0$. Note that the last term in our formulation will result in a penalty for non-zero values of the image in the air region, but it will not penalize non-zero values in the tissue region. Further, this soft formulation will lead to higher penalty for higher values in the air region, whereas the hard constraint formulation above directly forbids any non-zero values, regardless of the specifics. " ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [], "source": [ "def func2(x):\n", " lam = 5.0\n", " value = func(x) + lam*la.norm(C.dot(x))**2 \n", " return value" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "`gtol` termination condition is satisfied.\n", "Number of iterations: 9, function evaluations: 918, CG iterations: 14, optimality: 9.05e-07, constraint violation: 0.00e+00, execution time: 0.34 s.\n", "\n", "Final value of the cost function:\n", "0.03335052380816429\n", "\n", "Error with respect to original noiseless image:\n", "0.06488464617064679\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxddZ3/8dfn3pt93/elWbpSoG1aWhapbBZ0iisCDoI44oY4g+PI6Izj6Mz8Rh0XGHFBQNRRKqAiSGVfW2lp2tIt6ZKmbZImzdrs6733+/vj3oQ0TZo0vTfnLp/n49FHc+89N/dze5J3v/dzvud7xBiDUkqp4GezugCllFK+oYGulFIhQgNdKaVChAa6UkqFCA10pZQKEQ6rXjg9Pd0UFxdb9fJKKRWUtm/f3maMyZjsMcsCvbi4mMrKSqteXimlgpKIHJvqMW25KKVUiNBAV0qpEKGBrpRSIUIDXSmlQoQGulJKhQgNdKWUChEa6EopFSI00JWaQv+wk1+/eZTBEZfVpSg1I5adWKRUIBtxufn8b3bwyoFWEmMiuP7CPKtLUmpaOkJXagJjDPf8fg+vHGjFbhMqj560uiSlZkRH6EpN8J3nDvD7HQ38w1Xz2Xa0g8pjGugqOEw7QheRh0WkRUT2TvG4iMh9IlIjIrtFZLnvy1TKtzYdauPae9+gtrX3lPv/uLOBn7x6mJsvKuSuK8uoKE5h/4luugdHLKpUqZmbScvlEWDdGR6/Fij3/rkD+Mm5l6WU/4y43Hz9qb1UN3Xz+d/uHDvoWdfez78+uY+VxSl8c/0SRISVxakYAzt0lK6CwLSBbox5Heg4wybXA78yHluAZBHJ8VWBSvnahm311Lb2ceuaIqqbuvnWn6twutx88Xc7EYEffPRCHHbPr8aFBcnYbcJ2DXQVBHzRQ88D6sfdbvDe1zRxQxG5A88onsLCQh+8tFJnp2dwhB++cJBV81L5xvolREXYeeD1Wo6197OzrpP7blpGfkrs2PZxUQ4W5ySy7eiZxjRKBQZfzHKRSe4zk21ojHnAGFNhjKnIyJh0fXal/Opnr9XS3jfM165bhIjw5fcsYFlhMptq2vjg8jzWX5B72nMqilN4u76TYafbgoqVmjlfBHoDUDDudj7Q6IPvq5RPNXUN8OCmWtZfkMsFBckARNht/Phjy7nrynK+ef15kz6voiiVwRE3+xq75rJcpc6aLwL9KeDj3tkuq4EuY8xp7RalrPaXPScYHHFz99XzT7k/JymGu6+eT3zU5B3IiuIUAO2jq4A3k2mLjwJvAgtEpEFEPikinxGRz3g32QjUAjXAz4HP+a1apc7BsfY+EqIcFKXFTr/xOFmJ0RSmxmofXQW8aQ+KGmNumuZxA3zeZxUp5Sd1Hf0UpMYiMtlhnzOrKErhtYOtGGNm9Xyl5oKe+q/CRl1HP4WpZzc6H1VRnEp73zBH2/t9XJVSvqOBrsKC222oPzlw1u2WUWtK0wD4w44GX5allE9poKuw0NwzyLDTTcEsR+jz0uO4bmk2D286QkffsI+rU8o3NNBVWKjztkpm23IB+Ier5tM/4uJnrx32VVlK+ZQGugoLdR3nHujlWQm8/8I8fvnmUVq6B31UmVK+o4GuwkJ9Rz82gdzkmHP6Pn9/VTkjLsOPX9VRugo8GugqLNR19JObHEOk49x+5IvS4rihIp/fbq3jzcPteGbtKhUYNNBVWDh2DlMWJ/rCFeUkxji46edb+JsfbeKJ7Q243Brsynoa6Cos1Psw0HOTY3j9n97Nf37gPIZG3Pzj47v4z2eqffK9lToXGugq5PUNOWnrHZ71lMXJxEY6+NhFRTz/D+/i42uKeHjzEV472Oqz76/UbGigq5BXf9Izw2W2JxWdiYjw1esWMT8rnn98fBftvUM+fw2lZkoDXYW8Yz6Yg34m0RF27r1xGV39I3zl93v0QKmyjAa6Cnn1PpiDPp1FOYl85dqFvFjdzF/2nvDb6yh1JhroKuTVdfSTEO0gKSbCr6/ziYuLiYmwU3lU101X1tBAVyGvrqOforTZLZt7Nmw2oSQjjtq2Xr++jlJT0UBXIe9cls09WyUZ8Rxu1UBX1tBAVyHN5TY0dAz4dMrimZSkx9FwcoDBEdecvJ5S42mgq5DW3D3IsMs9ZyP00sx4jHlnZo1Sc0kDXYU0X6yyeDZK0uMAtO2iLKGBrkLC4IiL5kmWtD3W3gdAUWrcnNRRkuF5nVoNdGUBDXQVEn762mHWfvfVU0bGw043D206Ql5yDLnJ0XNSR2ykg9ykaGpb++bk9ZQaTwNdhYS69n4GRlzc/bu3GXG5AXjg9cMcbO7lW+9fgsM+dz/qOtNFWUUDXYWE1t4hYiPt7Gro4kcv11Db2st9L9fw3vNzuGJh1pzWUpIRR21rny4BoOacw+oClPKFtt5h1pSkkRQbwY9eqeHZvSeIctj4t79ZPOe1lKTH0TPkpLV3iMyEuWn1KAU6Qlchor13iPT4KL6xfgnZidEcaO7hq9ctsiRQSzPjATjcon10Nbc00FXQc7sN7X3DpMVHkhgdwc9uWcGX37OAj1YUWFJPSYYn0HUJADXXtOWigl7nwAgutyE9PgqA8/KSOC8vybJ6chKjiY6w6QhdzbkZjdBFZJ2IHBCRGhG5Z5LHC0XkFRHZKSK7ReQ635eq1OTavBeVSE+IsrgSD5tNKEmP1xG6mnPTBrqI2IH7gWuBxcBNIjLxSNO/AI8ZY5YBNwI/9nWhSk2lrccb6PGRFlfyjtGZLkrNpZmM0FcBNcaYWmPMMLABuH7CNgZI9H6dBDT6rkSlzqzVO0LPiA+METp4+uj1J/t1kS41p2YS6HlA/bjbDd77xvsG8Lci0gBsBL7gk+qUmoG23mGAsR56ICjNiNNFutScm0mgT3ZVgIlnTNwEPGKMyQeuA34tIqd9bxG5Q0QqRaSytVWvkK58o713CIdN/H5ForNROjrTRc8YVXNoJoHeAIyf/5XP6S2VTwKPARhj3gSigfSJ38gY84AxpsIYU5GRkTG7ipWaoK13iNS4SGw2/16R6GzMS49DBHY1dFldigojMwn0bUC5iMwTkUg8Bz2fmrBNHXAlgIgswhPoOgRXc6Ktdzig2i0AcVEOrliQyRPb6xlyah9dzY1pA90Y4wTuBJ4DqvHMZtknIt8UkfXezb4EfEpEdgGPArcZXchCzZG23qGAmbI43q0XF9PWO8wzu5usLkWFiRmdWGSM2YjnYOf4+74+7usq4BLflqbUzLT1DFHmPd0+kFxWnk5pRhyP/PUoH1iW5/eLVCulp/6roGaMoa13OKCmLI4SEW67uJjdDV3srO+0uhwVBjTQVVDrHnQy7HIHXA991AeX55MQ5eCRzUetLkWFAQ10FdTax077D5yzRMeLi3Jww8oCNu5pmvQSeUr5kga6CmqBeFLRRB9fU4TLGJ7Y3mB1KSrEaaCroDa6MFdaXOAGelFaHEWpsVQ1dltdigpxGugqqLUFeMtlVFlmPDUtetao8i8NdBXU2nqGEIHU2MAO9NKMeI609eH0XsBaKX/QQFdBrbV3mNTYSBz2wP5RLs2MZ9jlpuHkgNWlqBAW2L8FSk1j9FqigW70xCdtuyh/0kBXQc1z2n9gt1vgndUXa3T1ReVHGugqqAXiwlyTSYqJICMhSkfoyq800FVQa+sdCugpi+OVZehMF+VfGugqaPUPO+kfdgVFywU8ffTDrb3oQqTKXzTQVdBq6wn8s0THK8uMp2fQSav3otZK+ZoGugpagXhx6DMZOzCqbRflJxroKqh09Y9Q573w8tjCXEES6GNTF3Wmi/KTGV3gQqlA8a1nqvjDjgbuvKKcjHhP7zxYeuhZiVHERzl0hK78RgNdBZW36zuJjXRw30uHiI7wfMAMllkuIkKp98CoUv6gLRcVNAZHXNS29vKJS4r50c3LiHLYyUyIItIRPD/GOnVR+ZOO0FXQONjcg9vAopxErluaw+qSNHoGnVaXdVZKM+P4/Y4GugdHSIyOsLocFWKCZ2ijwl51k2c98UU5iYDnYOi89DgrSzprZd6ZLod1lK78QANdBY3qph5iI+0UpcZaXcqs6SJdyp800FXQqGrqZkF2AjabWF3KrBWmxhJpt1HVpFcvUr6nga6CgjGG6qbusXZLsHLYbVyxMJMntnv66Er5kga6CgrHOwfoGXQGfaAD3HlFGT2DTn7116NWl6JCjAa6CgrVTT0ALM5JsLiSc3deXhJXLMzkoU1H6BsKrlk6KrBpoKugMDrDZUF28I/QwTNKP9k/wm+2HrO6FBVCNNBVUKhu6qYoLZb4qNA4dWJ5YQqXlqXzwOtHGBxxWV2OChEzCnQRWSciB0SkRkTumWKbG0SkSkT2ichvfVumCnfVTd0sDoH++XhfuKKMtt4hNrxVZ3UpKkRMG+giYgfuB64FFgM3icjiCduUA/8MXGKMWQL8vR9qVWGqb8jJsY7+kDggOt5FJWmUZ8bz2sFWq0tRIWImI/RVQI0xptYYMwxsAK6fsM2ngPuNMScBjDEtvi1ThbP9J3owhpALdPAcIN1/osfqMlSImEmg5wH14243eO8bbz4wX0Q2i8gWEVk32TcSkTtEpFJEKltbdVSiZuadU/6Df4bLRAuzE2jqGqSzf9jqUlQImEmgT3Za3sSLIjqAcmAtcBPwoIgkn/YkYx4wxlQYYyoyMjLOtlYVpqqbukmMdpCXHGN1KT630PupQ0fpyhdmEugNQMG42/lA4yTb/MkYM2KMOQIcwBPwSp0TYwybatpYXpSCSPCe8j+VRdmeTx37dSkA5QMzCfRtQLmIzBORSOBG4KkJ2zwJvBtARNLxtGBqfVmoCk8Hm3s51t7PNYuzrS7FLzISokiNi9QRuvKJaQPdGOME7gSeA6qBx4wx+0TkmyKy3rvZc0C7iFQBrwBfNsa0+6toFT6e33cCgKsWZVpciX+ICAuyEqjWQFc+MKOzNIwxG4GNE+77+rivDXC3949SPvN8VTPLCpPJTIy2uhS/WZiTwIa36nG5DfYgXklSWU/PFFUBq7FzgD3Hu0K23TJqUXYiAyMu6jr6rS5FBTkNdBWwXqhqBuCaJVkWV+JfC3P0wKjyDQ10FbCerzpBaUYcpd7LtoWq8swEbKJTF9W500BXAamrf4SttR1csyS02y0AMZF2itPj2H9CR+jq3Gigq4D0yoEWnG7DNYtDu90yalF2oo7Q1TnTQFcB6YXqZjITorgg/7QTjkPSwuwEjrX36wUv1DnRQFcBaXdDJyuLU4P6gtBnY3QJgAPNOkpXs6eBrgJOz+AI9R0DIbkY11QWji0BoIGuZk8DXQWc0V5yKC6XO5X8lBjioxxjK0sqNRsa6CrgjIba4tzwCXQRoTwrnoPaclHnQANdBZzqpm6SYyPIDuHT/SdTlhHP4dZeq8tQQUwDXQWcqsZuFmUnhuRyuWdSlhlPW++wXuxCzZoGugooLrfhQHNPWPXPR5Vles6IrWnRUbqaHQ10FVCOtPUxOOIOq/75KA10da400FVACeXrh04nPyWWSIdNA13Nmga6CijVTd04bDI2Wg0ndptQkh5HjR4YVbOkga4CSlVTN2WZ8UQ57FaXYomyzHgdoatZ00BXAaW6qZvFYXhAdFRZZjzHOwcYGHZZXYoKQhroKmB09A3T3D0UljNcRpVlxmMMOh9dzYoGugoY7xwQDd9AL8/0HAzWQFezoYGuAkZVY/jOcBlVnB6LTXTqopodDXQVMKqbuslKjCItPsrqUiwT5bBTlBanga5mRQNdBYyqpu6wbreMKs3QmS5qdjTQVUAYHHFxqKWXJWF4huhEZZnxHG3vw+lyW12KCjIa6CogHGruxeU2LMlNsroUy5VlxjPiMhzr6Le6FBVkNNBVQNjX2AWgI3ROXdOlb8jJt5/dz2Pb6i2uSgUDh9UFKAWwr7GbhCgHBSmxVpdiudKMOAD+9PZx/uOZKuo7BliQlcANKwssrkwFuhmN0EVknYgcEJEaEbnnDNt9WESMiFT4rkQVDvY1drEoJzFsLgp9JgnRnot7bNxzArsIl5alU9fRjzHG6tJUgJs20EXEDtwPXAssBm4SkcWTbJcA3AVs9XWRKrS53Ibqpp6wXDJ3Kp++vITPrS3lL198F1ctymRgxEVr75DVZakAN5OWyyqgxhhTCyAiG4DrgaoJ230L+A7wjz6tUIW8I219DIy4tH8+zicumTf2dVGapwVT39FPZkJ4XZZPnZ2ZtFzygPFHZBq8940RkWVAgTHmz2f6RiJyh4hUikhla2vrWRerQtM7B0R1hstkClI9xxWOteusF3VmMwn0yZqaY808EbEBPwC+NN03MsY8YIypMMZUZGRkzLxKFdKqGruJtNsozwq/NdBnIj8lBhGo02mMahozCfQGYPzh9XygcdztBOA84FUROQqsBp7SA6NqKg0n+znZ986FkKuaupmfHU+EXWfRTiY6wk52YjR1OkJX05jJb9A2oFxE5olIJHAj8NTog8aYLmNMujGm2BhTDGwB1htjKv1SsQp6tz+yjZsf3MqIy40xhn2N4b0G+kwUpsbqCF1Na9pAN8Y4gTuB54Bq4DFjzD4R+aaIrPd3gSq0uNyG2tY+qpu6efCNI5zoHqSjb1j759MoTI3VM0fVtGZ0YpExZiOwccJ9X59i27XnXpYKVU1dAzjdhoRoBz988SARds8hGp3hcmZFabG09gwxMOwiJjI8L8+npqdNSzWnRtsG/75+CZEOG//vL/sRCe+LWszE6EwXbbuoM9FAV3OqoWMAgJXFqXz1ukW43IZ5aXHERekqFGcyOhddA12dif4WqTlV19GP3SbkJEVz48oCXj/YSmmGTlecTuHYXPQ+iytRgUwDXc2p+pP95CZH4/BOUfzJ366wuKLgkBIbQUKUg3odoasz0JaLmlN1Hf26ouIsiAgFOtNFTUMDXc2p+o7+sfaBOjtFaToXXZ2ZBrqaM/3DTtp6h8dmbKizU5gaS0PHAC63LqOrJqeBruZMvXeGiwb67BSmxTLsctPcPWh1KSpAaaCrOTPaLtCWy+wUpXqmLuqqi2oqGuhqztRroJ+T0X83nemipqKBruZMXUc/cZF2UmIjrC4lKOUmR2O3Ccc6dC66mpzOQ1dzpr6jn4LUWET0uqGz4bDbyEuOobqph4PNPbiNIScphqQY/Q9SeWigqzlTf7J/7BR2NTulGXG8vL+Fl/e3ALA4J5GNX7zM4qpUoNBAV3PCGENdRz+XleuVqs7Ff3xgKTuOncQmwmsHW3issoG23iHS46OsLk0FAA10NSdae4cYHHHrAdFzlJccQ15yDADZSdE8VtnAtiMdXLs0x+LKVCDQg6JqTrwzBz3G4kpCx9K8JKIjbGw90mF1KSpAaKCrOaFTFn0v0mFjeWEK245qoCsPDXQ1J0ZPKsrXhbl8atW8VKqauukeHLG6FBUANNDVnKjv6CczIYroCL18mi+tmpeKMbD96EmrS1EBQANdzYk6XWXRL5YVpBBhF+2jK0ADXc0BYwxH2vo00P0gJtLO+fnJvHWk3epSVADQQFd+V98xQEvPEBcWJltdSkhaNS+V3Q1dDAy7rC5FWUwDXfndFu/o8aJ5aRZXEppWzUvF6TbsrNM+erjTQFd+t7W2g9S4SMoz9WLQ/rCiKAWboH10pYGu/G9LbTurilOx2XRRLn9IjI5gcW4iW7WPHvY00JVfNZzs53jnAKtLUq0uJaStKUljx7FOeoecVpeiLKSBrvxqa62nDXBRifbP/emqRVkMu9y8dqDV6lKUhTTQlV9tPdJOcmwEC7ISrC4lpK0oSiElNoIXqk5YXYqy0IwCXUTWicgBEakRkXsmefxuEakSkd0i8pKIFPm+VBWMttR2aP98DjjsNq5YmMXL+1sYcbmtLkdZZNpAFxE7cD9wLbAYuElEFk/YbCdQYYw5H3gC+I6vC1XBp6lrgLqOfm23zJGrF2fRPehkm852CVszGaGvAmqMMbXGmGFgA3D9+A2MMa8YY0avXLsFyPdtmSoYjfbP9YDo3HjX/HSiHDaer2q2uhRlkZkEeh5QP+52g/e+qXwS+MtkD4jIHSJSKSKVra168CbUbaltJzHawcLsRKtLCQuxkQ4uLUvnhapmjDFWl6MsMJNAn6z5OelPi4j8LVABfHeyx40xDxhjKowxFRkZeimyUNU/7OTxynpeqGpm1bxU7No/nzNXL87ieOcA1U09VpeiLDCTS9A1AAXjbucDjRM3EpGrgK8BlxtjhnxTngomnf3D/PDFQzyxvYHeIScl6XF8dm2Z1WWFlSsXZSGyhxeqmlmcq5+Mws1MAn0bUC4i84DjwI3AzeM3EJFlwM+AdcaYFp9XqQKaMYZn9jTxjaf2cbJ/hOsvzOWmVYVUFKUgoqPzuZSREMWygmSe23eCu64s03//MDNtoBtjnCJyJ/AcYAceNsbsE5FvApXGmKfwtFjigce9P0B1xpj1fqxbBQhjDHdteJundzWyNC+JX96+iiW5SVaXFdY+uDyff3lyL68dbGXtgkyry1FzSKw6eFJRUWEqKysteW3lO0fb+lj7P69y28XF/Mt7F+Gw67lqVht2urnq+68RH+Xgz1+4VM8BCDEist0YUzHZY/rbp87JDu+SrTeuKtAwDxCRDhtfumY+VU3d/HlPk9XlqDmkv4HqnGw/dpL4KAflmXpqfyD5m/NzWZidwPeeP8CwU88cDRca6Oqc7Kjr5MKCZJ2aGGBsNuEr6xZyrL2f31XWT/8EFRI00NWs9Q45OXCim+V6abmAtHZBBquKU/nBCwepadF56eFAA13N2u76TtwGlhelWF2KmoSI8F8fXIpNhI/+bAvVTd1Wl6T8TANdzdr2Y54DossKNNADVVlmPI99ejURdhs3/XwLexq6rC5J+ZEGupq1HXUnKcuMJyk2wupS1BmUZMTz2KfXEBfp4GMPbuFQs7ZfQpUGupoVt9uws76TFYU6Og8GhWmxbLhjNVERdm77xTZaegatLkn5gQa64lh7Hx/6yV853jkw4+fUtvXR2T/C8iI9IBosClJjefjWlXT0DfPJRyrpH9brj4YaDXTFvS8eYvuxk2yuaZvxc0ZPKFquI/SgsjQ/if+9aRn7Gru469G3dZndEKOBHuaOtvXx5NvHATh4Yua91Z11J0mMdlCaEe+v0pSfXLU4iy+/ZyEvVjfzdn2n1eUoH5rJaosqhP341RocdhuZCVEcbOmdcru+ISdf3LCTxJgIKopS2VLbwbLCFF0nJEh9bHUhP3jxIE/uPM4y/ZQVMnSEHsbqO/r5w47j3LyqkJXFqWec/fDE9gZerG7h5f0tfPWPezjS1keFzj8PWonREVy9KIundzedclHpvce7+OQj2xgYdllYnZotDfQw9uNXD2MT4dOXlzA/K4GmrkG6BkZO287lNvxi8xEuLEhm579ezctfupz7b17OrZcUz33RymfevyyPjr5h3jjkuRykMYZ/f3ofL+1v4YBObQxKGuhh6kTXIE9sr+eGlfnkJMUwP8vTC5/sFPGXqps52t7P3102DxGhJCOe956fQ2K0zj8PZpfPzyA5NoInd3ouQPbKgRa2HfUc7G46ixlPKnBooIepzTVtjLgMt6wuBmB+lme1xIPNp/fRH9p0hLzkGNYtyZ7LEpWfRTpsvO/8HJ6vOkH34AjfefYAWYlRADR26Tz1YKSBHqb2n+gmymGjNCMOgLzkGGIj7Ryc8FF77/Euth7p4NaLi3S98xD0/gvzGBxx88VHd7L/RA9fvW4RUQ4bJ7p0hB6M9Dc0TFU39bAgO2EspG02oTwz/rRAf2jTEeIi7Xx0ZaEVZSo/W1GUQn5KDK8caGVxTiJ/c34uuckxOkIPUhroYWr/iW4WZp96UYryrIRTWi4t3YM8vauRj1QUkBSj/fJQJCJ8YFkeAF9etwCbTchJitYeepDSQA9xB5t7+MKjO0+ZhtbaM0Rb7zALsxNP2XZBVgKtPUOc7BsG4A87j+N0G25ZUzSnNau59ZnLS3n4tgrWzs8AICcphiYdoQclDfQg1tI9yAtVzWfc5tt/2c/TuxrZUts+dt/+E551sRfmTByhe2a6HGzuwRjD45X1rChK0bNBQ1xclIMrFmYh4jlJLDc5mubuQZwuvXRdsNFAD2Lfee4An/pVJS3dk4+mqhq7eWl/CwBvjg/0Jk+ffNGEEfrYTJeWXt6u7+Rwax8fWZHvj9JVAMtJisFtoKVnyOpS1FnSQA9SgyMuntt7AoBNUyyqdf+rNcRHOViSm8ibh98J9OoT3WQnRpMSF3nK9jlJ0SREOTjU3MPj2xuIjrDx3vNz/PcmVEDKSY4GoMlHM12cLjd/v2Hn2IJuyn800IPUqwda6BlyIgKbDp0e6LWtvWzc08Qta4q4enEW+xq7xs4C3d/Uc1q7BTwHyMqz4tnV0MXTuxq59rwcEvTkobCTk+QJ9MZO3/TRdx/v4sm3G3nojSOn3O90uVn/o008svnIFM9UZ0sDPUg9tauR9PhIrjsvh9cPtZ22DOpPXj1MpN3G7ZfMY3VJGm4Dbx3pYMTlpqal97QDoqPmZyWwq76TnkGntlvCVE5SDDD1CP1Qcw/X3vsGv9/eMKPvN3r85uX9LaccnN98uJ3dDV388KVD9A3p2uy+oIEehHoGR3ipuoX3Ls3h8gUZtPUOsX/c0rfHOwf4487j3LSqkIyEKJYVJhPlsPHXw23UtvYx7HKzaJIROrzTR89PiWF1SdqcvB8VWBKjHcRF2qec6fIfz1RT3dTNlx7fxdf+uIch55kX8nrzcDtRDhsDIy5eO9gydv+TO48T5bDR2T/Cb7fW+fQ9hCsN9CD0/L5mhpxu1l+Yy2Xl6cCpbZefvnoYEbjjXSUARDnsrChK4c3D7e/McDnDCB3gQ8vzdWncMCUi5CTH0DRJy+X1g628drCVe65dyGcuL+U3W+u44advUtfeP+n3Gna6qTx6kg+vyCclNoKNezzHffqGnDy79wQfWpHPxaVpPPBGLYMjusLjuZpRoIvIOhE5ICI1InLPJI9HicjvvI9vFZFiXxcaSPYe7+JjD27hl389asnrP7WrkbzkGJYXppCTFENZZjyve1fMa+oa4Hfb6vnwit0PmGIAAA7oSURBVAJyk2PGnrOmJI39J3r4a007kXYbJd5T/idaNS+Vu64o47aLi+firagAlZMUfVrLxeU2/NfGagpTY/nEJcXcc+1CfnbLCmrb+lh37+v8dmvdaa2/Pcc7GRhxcVl5Ou9Zks1L1c0Mjrh4vuoEAyMuPrAsjzuvKKO1Z4jHK+vn8i2GpGkvcCEiduB+4GqgAdgmIk8ZY6rGbfZJ4KQxpkxEbgS+DXzUHwWP93Z9Jxv3NFHb2sfR9j46+oYpTotlQXYii3MSWFOaTmlG3Nj82plo7Rmi8mgH5xckkzcuEAEGhl388MWDPLjpCHYRNte00zvk5PPvLpv0e4243AhMuQaKMYaX97fw1pEOTnQPcqJrEIddWJyTyJLcJC4pSycjIeqU57T3DrGppo073lUy9r4uK0/nt1vrGBxx8dNXD+M2hs+tLT3leWtK0+AFePLt45RlxhMxRU2RDht3X7NgJv9UKoTlJsWc0sYDeLyynv0nerj/5uVEOewAvGdJNkvzkvjyE7v46h/38HzVCe796DKSYj0H00dnV100L42YSAcbttXzxqE2/rizkfyUGFYUpiACywuT+elrtdy4qnDSn81hp5utR9p59UArGQlR3LK6iLgo/12fp2tghLfrO+kaGOG687JntI5R9+AIu+u72H28k73HuzjY3ItdhNgoOymxkdx99XzOy0vyW80wsysWrQJqjDG1ACKyAbgeGB/o1wPf8H79BPAjERHjhwsWut2Gl/a38PPXa3nraAeRdhtFabGUpMdRUZRCbVsfG/c08ehbnp5cTlI0a0rTKEyNJTMhmoyEKCLsgsNmQwS6B0bo6B+muWuQ1w61sct7SS4Rz/KiH16RT8+gk21HO9hc00Zz9xA3rizgn9Yt5N+f3sd3nzvAsNPNZ9eWUtvax6GWHvY0dLGzvpM9x7sYcblJjY0kIyGK8qwELi1L45KydKqbevjhiwfZ19hNpN1GVlIU2YnR9Ay6+OWbxxh2uslIiOLlL11+ykyTp3c14nIb1l+QO3bfu8oz+MXmozyzu4lHt9Xz4RX5FKTGnvLvdn5+MjERdgZGXJPOcFFqvJzkaNp6hxh2uol02OgbcvK9Fw6yoiiF65aeuupmbnIMv779In715lG+9Uw133vhAN+8/jzAc/7DwuwEUuIiubg0jaSYCH715lE217TxubVlY229O68o4/ZHKnng9Vo+t7Z0bLDS1T/Cfz9bzdO7mugdchJptzHscvPz12v57NpS/nZ1EdER9knfw5DTxSv7WzjeOciinASW5CYR5bCx93gX24+dZMjp5sZVBWQmeGb1DDvd/N+WY2zYVsehll5G0+vRkjT+9+ZlpMdHMeJy88jmo2zYVkd8dASZCVHERdqpauo+5TmFqbEszE5ABPqHXexu6OTjD7/FY59eQ1mm/07Um0mg5wHjPws1ABdNtY0xxikiXUAacMp8OhG5A7gDoLBwdos9/fDFg9z3cg15yTH8y3sX8dGVBadNrTPGUN8xwKaaNt441MrrB9to653+JIkLCpK5++r5XDQvlc01bfyusp47f7sTgPT4SFYUpXDrxcVcXOrpW3//hguJsNu496VD3PfyobGdGeWwsTQviY97RxGtvUO0dA+xpbadp3c1jr1eYWos3/3w+XxgWd4pI4ARl5vNNW3c9ott/PjVw3xl3UIAOvuH+d+Xa1hemHzKOiwXlaQSYRe+/qe9uN1m0k8MkQ4bFcUpvHGo7bQTipSaKDcpBmOguXuQgtRYntndRGvPEPffvHzST7w2m3DbJfM41NLLb7fW8XeXlpCVFMX2Yye5aZXndz3CbuOaxVk87p0d8/5l7wxK3r0gk7ULMvjucwfYUtvOf31gKYdbe/nK73fT1jvMB5flcc2SbC4p87QOv/f8Af7jmWq+8+wBFucmcmFBMkVpsdhtgt0m7Gvs5s+7GukePHX2jMMmON2eX1QRuP+VGm5aVch5eUnc99Ih6jr6WVmcwt1XzWd5UQoNJ/v5+p/28b77NvGFK8v4xeaj1LT0smpeKlEOG3Xt/XQPjrAwO4H3nZ/LssJkzs9LHvuEMupIWx8f+elf+fhDW3n8sxef9unfV2YS6JP1KyaOvGeyDcaYB4AHACoqKmY1ev/QinzKshLO+DFIRChMi+XmtEJuvsjzwzTsdNPWO0Rb7xAjLoPbGNxuQ0J0BClxEaTERp7yP/1FJWncdWU524+dJDMxmuK02NN+kO024TsfOp/FOYl0DoxQnhlPeVY8JenxRDpOr80Yw8HmXjbXtJESF8H7zs+d9ONlhN3G2gWZfHB5Hg9tOsLNqwopSI3l28/up3NghF+/f+kptcRGOlhRlMKW2g5uqDh9dD5qTWkabxxq0xG6mtboyUWNnQMUpMby3L4T5CXHsLL4zJcdvOvKcn6/o4Hvv3CAmy8qYnDEzZpxs6WuW5rD49sbWJqXRFnmOz+HIsLDt67kN1uP8d9/2c9V33+NIaeb8sx4fv7xCs7PTx7bdnlhCr/5u9VsqW3n5f0tvF3Xye+21TMw7qBqTISddedl8/5leSzKSaC6qYe9x7voHXJyYUEyywtT6B1y8uNXavj1lmO43IaF2Qn88vZVvKs8/ZTfr/PykvjM/23na3/cS1FaLA/dWsGVi7LO6t9zXnocv7x9FTf+bAu3PLSVxz+9hrT4qOmfeJZmEugNQMG42/lA4xTbNIiIA0gCOnxS4QRFaXEUpU1+QO9MIh02cpNjTjlQOB2H3cZF00zds9mE2y+dN6PvJyIsyE5gQfbMAvXL71nAxj1N/Pez+/nExcU8+lY9n7psHotzTx9hX7Uoix3HOqfs54Nn5srxkwOsLE6d0eur8DV6clFT1yB9Q07eqGnjYxcVTns8Kisxmk9cMo+fvnaYYZcbEU//fNQlZeksyU3k9kuLT3uuzSbcsqaYKxZl8V8bqylKjeWuK8unbKmsLkkbm1rrdLnpGXTi8g7UEmMiTnleZkI0l3sXHxuVkRDFdz9yAXddWc7R9j4uLk3HPsnMriW5Sfz5zst4/VArVy/OmrKe6SzJTeLhT6zkloe28qe3G2ecG2dDpmtzewP6IHAlcBzYBtxsjNk3bpvPA0uNMZ/xHhT9oDHmhjN934qKClNZWXmu9Ye8H7xwkHtfOkROUjQCvHD35ZMeDHK63LT1DpPt/UVU6lz0DTlZ8m/P8ZV1CylKi+Vzv9nBhjtWz+jchK7+ES77zst0DzpZkpvIM3ddNgcVB4+jbX0UTfKJf6ZEZLsxpmKyx6Y9dGuMcQJ3As8B1cBjxph9IvJNEVnv3ewhIE1EaoC7gdOmNqrZ+fTlJWQlRtHUNcg31i+Z8si+w27TMFc+ExflIDHaQVPXAM/tO0FKbAQVRWdut4xKio3gs2s9nxTX6MlppylOP7uZd2djRvN+jDEbgY0T7vv6uK8HgY/4tjQFnv74j25ezq76Tq7Ra3qqOZSbHMOx9n521J1k3ZKZTd0bddvFxdS09PKRioLpN1Y+47+JnMpnVhanat9bzbmcpGjeONSG0214z1kOJmIi7Xzvhgv8VJmaip76r5SaVE5yDE63ITbSzqXeJSZUYNNAV0pNKtd7TGbtgoxZz+xQc0sDXSk1qdFldM+23aKso4GulJrUuxdm8qnL5nHNYg30YKEHRZVSk0qNi+Rr711sdRnqLOgIXSmlQoQGulJKhQgNdKWUChEa6EopFSI00JVSKkRooCulVIjQQFdKqRChga6UUiFi2gtc+O2FRVqBY7N8ejoTrlcaBvQ9hwd9z+HhXN5zkTEmY7IHLAv0cyEilVNdsSNU6XsOD/qew4O/3rO2XJRSKkRooCulVIgI1kB/wOoCLKDvOTzoew4PfnnPQdlDV0opdbpgHaErpZSaQANdKaVCRNAFuoisE5EDIlIjIvdYXY8/iEiBiLwiItUisk9Evui9P1VEXhCRQ96/U6yu1ZdExC4iO0Xkz97b80Rkq/f9/k5EIq2u0ZdEJFlEnhCR/d59vSYM9vE/eH+m94rIoyISHWr7WUQeFpEWEdk77r5J96t43OfNs90isvxcXjuoAl1E7MD9wLXAYuAmEQnFS6o4gS8ZYxYBq4HPe9/nPcBLxphy4CXv7VDyRaB63O1vAz/wvt+TwCctqcp/7gWeNcYsBC7A895Ddh+LSB5wF1BhjDkPsAM3Enr7+RFg3YT7ptqv1wLl3j93AD85lxcOqkAHVgE1xphaY8wwsAG43uKafM4Y02SM2eH9ugfPL3oenvf6S+9mvwTeb02Fvici+cB7gQe9twW4AnjCu0movd9E4F3AQwDGmGFjTCchvI+9HECMiDiAWKCJENvPxpjXgY4Jd0+1X68HfmU8tgDJIpIz29cOtkDPA+rH3W7w3heyRKQYWAZsBbKMMU3gCX0g07rKfO6HwD8Bbu/tNKDTGOP03g61fV0CtAK/8LaZHhSROEJ4HxtjjgP/A9ThCfIuYDuhvZ9HTbVffZppwRboMsl9ITvvUkTigd8Df2+M6ba6Hn8RkfcBLcaY7ePvnmTTUNrXDmA58BNjzDKgjxBqr0zG2ze+HpgH5AJxeFoOE4XSfp6OT3/Ogy3QG4CCcbfzgUaLavErEYnAE+a/Mcb8wXt38+jHMe/fLVbV52OXAOtF5CieNtoVeEbsyd6P5hB6+7oBaDDGbPXefgJPwIfqPga4CjhijGk1xowAfwAuJrT386ip9qtPMy3YAn0bUO49Kh6J54DKUxbX5HPe/vFDQLUx5vvjHnoKuNX79a3An+a6Nn8wxvyzMSbfGFOMZ5++bIz5GPAK8GHvZiHzfgGMMSeAehFZ4L3rSqCKEN3HXnXAahGJ9f6Mj77nkN3P40y1X58CPu6d7bIa6BptzcyKMSao/gDXAQeBw8DXrK7HT+/xUjwfu3YDb3v/XIenr/wScMj7d6rVtfrhva8F/uz9ugR4C6gBHgeirK7Px+/1QqDSu5+fBFJCfR8D/w7sB/YCvwaiQm0/A4/iOUYwgmcE/smp9iuelsv93jzbg2cG0KxfW0/9V0qpEBFsLRellFJT0EBXSqkQoYGulFIhQgNdKaVChAa6UkqFCA10pZQKERroSikVIv4/kBV5NOL5x+MAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x0 = 0*xn\n", "res = opt.minimize(func2, x0, method='trust-constr', tol=1e-12, options={'gtol': 1e-6, 'disp': True})\n", "plt.plot(res.x)\n", "\n", "print('\\nFinal value of the cost function:')\n", "print(res.fun)\n", "\n", "print('\\nError with respect to original noiseless image:')\n", "print(la.norm(res.x-x)**2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }