mirror of
https://github.com/wassname/simpeg.git
synced 2026-08-02 13:00:17 +08:00
218 lines
20 KiB
Plaintext
218 lines
20 KiB
Plaintext
{
|
|
"metadata": {
|
|
"name": ""
|
|
},
|
|
"nbformat": 3,
|
|
"nbformat_minor": 0,
|
|
"worksheets": [
|
|
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"collapsed": false,
|
|
"input": [
|
|
"import SimPEG\n",
|
|
"from SimPEG.mesh import TensorMesh\n",
|
|
"from SimPEG.regularization import Regularization\n",
|
|
"import SimPEG.inverse as inverse\n",
|
|
"from SimPEG.inverse import Minimize\n",
|
|
"import numpy as np"
|
|
],
|
|
"language": "python",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"prompt_number": 63
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"collapsed": false,
|
|
"input": [
|
|
"from SimPEG.forward.LinearProblem import example as LinExample\n",
|
|
"\n",
|
|
"prob, m_true = LinExample(1000)\n",
|
|
"M = prob.mesh\n",
|
|
"\n",
|
|
"reg = Regularization(M)\n",
|
|
"opt = inverse.InexactGaussNewton(maxIter=20)\n",
|
|
"inv = inverse.Inversion(prob,reg,opt,beta0=1e-4)\n",
|
|
"m0 = np.zeros_like(m_true)\n",
|
|
"\n",
|
|
"mrec = inv.run(m0)\n",
|
|
"\n",
|
|
"# plt.figure(1)\n",
|
|
"# for i in range(prob.G.shape[0]):\n",
|
|
"# plt.plot(prob.G[i,:])\n",
|
|
"\n",
|
|
"plt.figure(2)\n",
|
|
"\n",
|
|
"plt.plot(M.vectorCCx, m_true, 'b-')\n",
|
|
"plt.plot(M.vectorCCx, mrec, 'r-')\n"
|
|
],
|
|
"language": "python",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"output_type": "stream",
|
|
"stream": "stdout",
|
|
"text": [
|
|
"====================== SimPEG Inversion ======================\n",
|
|
" # beta phi_d phi_m f norm(dJ) #LS\n",
|
|
"--------------------------------------------------------------\n",
|
|
" 0 1.00e-04 9.53e+02 0.00e+00 9.53e+02 8.95e+03 0\n",
|
|
" 1 1.00e-04 6.17e-01 2.65e+04 3.27e+00 1.58e+01 0\n",
|
|
" 2 1.00e-04 9.22e-02 2.30e+04 2.39e+00 1.33e+02 0\n",
|
|
" 3 1.00e-04 5.11e-02 2.29e+04 2.34e+00 1.02e+00 0"
|
|
]
|
|
},
|
|
{
|
|
"output_type": "stream",
|
|
"stream": "stdout",
|
|
"text": [
|
|
"\n",
|
|
" 4 1.00e-04 4.08e-02 2.28e+04 2.32e+00 1.38e+01 0\n",
|
|
"------------------------- STOP! -------------------------\n",
|
|
"1 : |fc-fOld| = 1.4678e-02 <= tolF*(1+|fStop|) = 9.5423e+01\n",
|
|
"1 : |xc-xOld| = 5.5359e-02 <= tolX*(1+|x0|) = 1.0000e-01\n",
|
|
"1 : |g| = 1.3766e+01 <= tolG*(1+|fStop|) = 9.5423e+01\n",
|
|
"0 : |g| = 1.3766e+01 <= 1e3*eps = 1.0000e-02\n",
|
|
"0 : iter = 4\t <= maxIter\t = 20\n",
|
|
"========================= DONE! =========================\n",
|
|
"\n"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"output_type": "pyout",
|
|
"prompt_number": 16,
|
|
"text": [
|
|
"[<matplotlib.lines.Line2D at 0x110646dd0>]"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"output_type": "display_data",
|
|
"png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVHX+P/DXUQyvIYjX9RZoCgJCImBcvV9Ia3/6LbG0\n1drFS2qabfnNXbXatL6rlWRk33L9tl7yWqG2KdYCpol4x9JUAm29ggRqigLz+f3xydFxZpj79bye\nj8c8ambOnPPmeHjNh8/5nM9RhBACRESkGvVcXQARETkXg5+ISGUY/EREKsPgJyJSGQY/EZHKMPiJ\niFTG5uCfMGECWrdujfDwcIPv5+TkwM/PD1FRUYiKisLrr79u6yaJiMgGPrauYPz48Zg6dSrGjRtn\ndJnk5GRkZWXZuikiIrIDm1v8iYmJ8Pf3r3MZXiNGROQ+HN7HrygKdu/ejcjISMycORNFRUWO3iQR\nEdXB4cH/0EMP4eeff0ZBQQFCQ0Mxffp0R2+SiIjqIuyguLhYhIWFmVxOo9GIVq1aiaqqKr33goOD\nBQA++OCDDz4seAQHB1uc2Q5v8V+8eFHbx79582ZERETA19dXb7mioiIIIfgQAnPnznV5De7y4L7g\nvuC+qPthTfe5zaN60tLSkJubi7KyMnTo0AHz589HdXU1ACA9PR0bNmxAZmYmfHx8EBERgUWLFtm6\nSSIisoHNwb9mzZo6358yZQqmTJli62aIiMhOeOWuG0pJSXF1CW6D++IO7os7uC9sowghhKuLAOSw\nTzcphYjIY1iTnWzxExGpDIOfiEhlGPxERCrD4CciUhkGPxGRyjD4iYhUhsFPRKQyDH4iIpVh8BMR\nqQyDn4hIZRj8REQqw+AnIlIZBj8Rkcow+ImIVIbBT0SkMgx+IiKVYfATEakMg5+ISGUY/EREKsPg\nJyJSGQY/EZHK+Li6AFIZIYADB4DiYiA4GIiMBBTF1VURqQqDn5zn4EFgwgTg11+BHj2AwkIgMBBY\nvhwIDXV1dUSqwa4eco7sbGDQIGDmTODHH4HPPgNOnJBfBCkpwJ49rq6QSDUUIYRwdREAoCgK3KQU\nsrfCQqBfPxn2CQn672/dCjz7LPDdd0Dnzk4vj8iTWZOdDH5yrBs3ZD/+nDnA2LHGl3vrLeDLL4Fv\nvgHq8Q9RInNZk538DSPHeu01IDy87tAHgBdeAKqrgY8/dk5dRCrGFj85zvHjQGKi7Opp08b08vv2\nAY8+Ks8BNG3q+PqIvABb/ORe5swBZs0yL/QBIDoaSE4GFi92bF1EKscWPznG/v3AiBHAyZNA48bm\nf+7ECSA+HigpAZo0cVh5RN6CLX5yHwsWAC+9ZFnoA8CDD8qRPytWOKQsImKLnxyhuBjo3Vu22q3p\nq9+9W54MPnmSI3yITGCLn9zDkiXywixrT9A+/DDg5wd8/bV96yIiAGzxk71duQI88ABw6BDQoYP1\n61m6FNi5E/j0U/vVRuSFXNLinzBhAlq3bo3w8HCjy8yePRtBQUHo1asXjh8/busmyZ2tXStH5tgS\n+gAwZgzw1VfA5cv2qYuItGwO/vHjx+Orr74y+v7evXuxc+dO7Nu3D7NmzcKsWbNs3SS5s3/8Q3bz\n2MrfHxg6FFi3zvZ1EZEOm4M/MTER/v7+Rt/Pz8/HqFGjEBAQgLS0NBw7dszWTZK7OnZMntgdMsQ+\n63viCWD9evusi4i0HD4t8969ezH2rsv1W7ZsiaKiIgQHBzt602QljQb49ls5g4Ilunz4D4jkcSjK\ntc9hVc93CBIKxmPP+gu4FWDmRWD3CA7mvG9E93J48Ash9E48KEZuvDFv3jzt/6ekpCAlJcWBlZEx\nx44Bw4YBsbHmf0YRGny6axX+HJWN02/Yq5KGmN0kFd/P2YSs9pMt/nRZGdCpE5CVZa96iFwvJycH\nOTk5Nq3DLqN6SkpKMHz4cBQWFuq9l5GRgZqaGsyYMQMAEBwcjKKiIv1COKrHbRw+DIwbJ/9rtp07\ngSlTgCNH7FtMVpacwsGKA/3LL4H33pP/JfJWbjmOPzY2Fhs3bsTly5exevVqhISEOHqTZCMhrLgb\n4rp1sk/e3gYNknfuKi21/7qJVMrmrp60tDTk5uairKwMHTp0wPz581H9W+dweno6YmJikJCQgOjo\naAQEBGDlypU2F02OZXHw19YCGzYAeXn2L6ZhQ6B/fzm009TUzvdQFPmzEJEum4N/zZo1JpdZuHAh\nFi5caOumyEksDv6dO4G2bYGuXR1TUGoqsGULg5/ITjhlA+mxOPjXrgUef9xh9WDYMGD7dsuHGRGR\nQQx+0mNR8Gs08l66o0Y5rqC2beW4zN27LfoYW/xEhjH4SY9FYVlQALRoAXTp4rB6ANzp7rEQg59I\nH4Of9FjU4t+8Wd5wxdGGDJHdPRaweGQSkUow+EmPRcGflQUMH+7QegDcmd//0iWzP8KuHiLDGPyk\nx+zgLykBLlyw7BJfa/n4AElJwL//7fhtEXk5Bj/pMTv4N28GHnkEqF/f4TUBAAYMsOjmLGzxExnG\n4Cc9Zge/s7p5buvfn8FPZAcMftJjVvBfuQLk5wMDBzqlJgBAjx7Ar7/KqZ+JyGoMftJjVvB/8w0Q\nF2f9fXWtoSgWtfrZ4icyjMFPeswK/u3b5QRqzta/v/zSMQODn8gwBj/pMSv4s7Od281zW1KSnBuI\niU5kNQY/6TEZ/MXFwNWrQHi402rSCg4GamqA06dNLsoWP5FhDH7SYzL4s7Pl0Mp6Ljh8FAVISJD3\nhjQDg59IH4Of9JgV/K7o5rktMVF295jAKRuIDGPwk546g7+2Vp5cdWXwJySYHfxs8RPpY/CTnjqD\nf/9+OU1yu3ZOrUlHRARw9qy8m7oJDH4ifQx+0lNn8Lu6mweQ8/bExQG7dtW5GLt6iAxj8JMetw9+\nQPbzmzjBy64eIsMY/GSQweC/dk129SQnO70ePWb28xORPgY/6THa4s/NBaKjgSZNnF6TnthYoLAQ\nuH7d6CJs8RMZxuAnPUaD3126eQCgUSMgNBQ4cMDoIgx+IsMY/KSnzuB3xfw8xsTGyhlCicgiDH7S\nYzD4//Mf4OJFICrKJTUZFBcH7Nlj9G22+IkMY/CTHoPBv2MH0K+f8+62ZQ4TLX4GP5FhDH7SYzD4\n3al//7YuXeSNWc6dc3UlRB6FwU969FrJGo17Br+iyO4eI61+tviJDGPwkx69Fv+RI0Dz5kDnzq4q\nybjY2Dr7+Rn8RPoY/KRHL/jdsbV/m4kWPxHpY/CTHoPB707DOO8WEyOvJq6p0XuLXT1EhjH4SY9O\n8N+4AXz3HZCS4sqSjGveHGjfHjh61ODbDH4ifQx+0qMT/N9+K6dB9vNzaU11MjKsk109RIYx+EmP\nTvC7czfPbUYu5GJXD5FhDH7SoxP827e774nd22JigH37XF0Fkcdg8JMebfBfvAiUlMhgdWdhYUBR\nkbyY6y5s8RMZZnPw5+XlISQkBF27dkVGRobe+zk5OfDz80NUVBSioqLw+uuv27pJcjBt8H/9tTyp\n6+Pj6pLqdt99MvwPHdJ5mcFPZJjNv9HTp0/HsmXL0KlTJwwePBhpaWkIDAzUWSY5ORlZWVm2boqc\nRBv827e7f//+bdHRQEEBEB/v6kqI3J5NLf7KykoAQFJSEjp16oRBgwYh38DoCsFml0cRAlAg3PvC\nrXtFR+v187PFT2SYTcFfUFCA7t27a5+HhoZizz2jKxRFwe7duxEZGYmZM2eiqKjIlk2SEwgBtL96\nDGjQQE6E5gkMBD/A4CcyxOEndx966CH8/PPPKCgoQGhoKKZPn+7oTZKNhAB6XtwODB7sOYPhQ0Pl\nPQOuXNG+5CmlEzmbTX38vXv3xosvvqh9/v3332PIkCE6yzRr1kz7/8888wxeeeUV3Lx5E76+vnrr\nmzdvnvb/U1JSkOKuV4t6OSGAiIvbgUHPuLoU8/n4AD17ylsx/nbcsKuHvFFOTg5ycnJsWodNwe/3\n29WceXl56NixI7KzszF37lydZS5evIhWrVpBURRs3rwZERERBkMf0A1+ch3l1k10K/0W6LfK1aVY\n5nZ3z10NBgY/eZt7G8Xz58+3eB02j+p55513kJ6ejurqakybNg2BgYFYtmwZACA9PR0bNmxAZmYm\nfHx8EBERgUWLFtm6SXKwwB934axfD3T193d1KZaJjga2btU+ZVcPkWGKcJMhN4qicPSPmzj6yMs4\nccYX/++I5S0Jlzp2DBg+HDh1CoCct230aKPztxF5BWuyk1fukp52hdtwtK2HjN+/24MPApcuAb/8\non2JbQkifQx+0nXxIpqWFqO4pZtP02BI/fpAVJScnx/s6iEyhsFPunbswPmQvtDUb+DqSqxz13h+\njuohMozBT7q2b8fZHoM8t7V8z4VcDH4ifQx+ukMIGfyh3hH8HvszEDkYg5/uOHoUaNIEV1oGu7oS\n63XpAlRWAqWlANjiJzKEwU93bNsGDByof7N1T6IoQK9ewP79nvszEDkYg5/u+PJLYOhQAB4c/IC2\nu4cnd4kMY/CTVFkp+8b79/fsFj9wZ25+MPiJDGHwk7R9O5CQADRp4h3Bz64eIqMY/CRt3QqkpgKA\n5wd/p05AVRV8Ss+zxU9kAIOfAI1G9u97S/ArChAdjYbf72fwExnA4CfZH96qFdC5MwAvCH5ABv9R\n/TtyERGDnwBgyxbgkUe0T70i+Hv1gu/RfWzxExnA4Ced/n3AS4I/Ohq+hfsgNEx+onsx+NXu7Fng\n9GmgTx/tS14R/O3bA0KgTe1ZV1dC5HYY/Gr3+eeyte9z52ZsXhH8ioJb4dEIq2I/P9G9GPxqt2ED\nMHKkzkteEfwAboZHI/wmg5/oXgx+NSstBQ4eBAbp3m3LW4L/VkQ0Im4x+InuxeBXs88/BwYPBho1\n0nnZa4I/vBfCbu3nvA1E92Dwq9nGjXrdPID3BL+mTTtUK/cBZ864uhQit8LgV6tffgG++w4YNkzv\nLW8JfkUBjjTQvSMXETH41SsrC+jXD2jaVO8tbwl+ADhyH4Of6F4MfrX69FPgv/7L4FveEvyKAhxm\ni59ID4NfjS5cAPbsAR57zODb3hT8Rxr0ksHPE7xEWgx+NVqzBnj0UaBxY4Nve0vwA0BZvVZAs2bA\nTz+5uhQit8HgV6N//hMYO9bo294S/NqfIZrdPUR3Y/CrzdGjwKVLQEqK0UW8JfiB33p4GPxEOhj8\navPJJ8CTTwL16xtdxFuCX3uzdQY/kQ4f04uQ17h5E1ixAvj22zoX85bzoNovr169gAMH5J3G6rGt\nQ8TfAjXZtAmIiAAefLDOxbylxQ/89iXWooV8nDzp6nKI3AKDX00yM4FJk0wu5i3Br+3qAdjdQ3QX\nBr9aHD0KnDoFjBhhclFvCn6t6Ghg/36X1ULkThj8avH3vwNTpgANGphc1FuCH2CLn8gQntxVg59/\nlnPzFBWZtbi3BL/Oz/DQQ/LeA7W1dY5oIlIDBr8tLl4Edu8GTpyQU/9WVQG+vkCbNkCPHkBcHPC7\n37m6SuDtt4Hx4wF/f7MW96bg17b4mzeX/y4//giEhrq0LiJXszn48/LykJ6ejpqaGkybNg1Tp07V\nW2b27NlYu3Yt/P39sWrVKnTv3t3wysaPl0PvunUDBgxwz/QpKgJWrgTWrwfOnQMefhgICZGPhg1l\n+J8/D/zf/wHp6UBQEDB6NDBhggwfZ7twQdZy+LDZH/GW4AfuGZp6u7tHzcF/4wbw1VfAtm3AkSNA\ncbF8rX59oHVroHNnuZ9iYuRFfgZmb/VoN24AZ8/K39Fff5XPr1+XP3/DhvLRrJlsJLRta3RaE09n\nc/BPnz4dy5YtQ6dOnTB48GCkpaUhMDBQ+/7evXuxc+dO7Nu3D9u2bcOsWbOwZcsWwyuLiZEHY2am\nDNXHHgPGjZMHoCuTSAggLw944w3ZXfDEE8BHHwG9e9fdbVBTA+TkyOANCpLhP3u2HFroLK++Cvzh\nD0D79mZ/RAjvGO6ud8jcDv5x41xSj0uVlwNvvQV8/DHQsycwdCgwZow8Lps0kcfqxYtyAMC+ffKv\nxDFjgPh44Pe/Bx5/3DUNF2vduAHk58sGT2GhHNxw8qQM+3btZKg3bSqDvVEjeY1HVZV8VFTIfXH+\nPHDffXIfdesmH927yyHRISGe3WUobFBRUSEiIyO1z6dOnSq2bNmis8ySJUvE22+/rX0eFBRkcF06\npWg0Qpw+LcTbbwsRGipEjx5CLF8uxK1btpRrOY1GiK1bhYiPF6JLFyE++kiImzetW9eZM0JMmiRE\ny5ZCvPuuENXV9q3VkBMnhAgMFKKszKKPzZkjxLx5DqrJiUpLhQgIuOuFnBwh+vRxWT0uodEI8fHH\nQrRqJUR6uhCnTpn/2cpKIdatE2LUKCH8/IR44gkh/vUvIWpqHFevtWprhThwQIg33xRiwAAhmjYV\nIi5O/s69/74QO3cKcemS3B/m0miEKC8XYt8+IVauFOIvf5H7oksXuf6EBCGef16IVauE+Okny9Zt\nR9bEuE3Bn52dLUaPHq19npmZKebMmaOzzFNPPSW2bdumfR4bGytOGTj4jBav0QiRnS1E//5CtG8v\n/2F/+cWWsk2rqZEHfGSkEOHhQqxZY7+gLiwUol8/IXr1EuLoUfus0xCNRojhw4V44w2LPzpnjhDz\n5zugJicrKxPC3/+uFyorhWjc2Dlfuu7gyhUhRo8WIiJCiEOHbFvX5ctCLF0qRO/eQrRrJ8RLLwlx\n7Jh96rTWzz/LBmFammxQPfigEFOmCPHFF/Lf2pHKy4XYsUOIhQuFGDlSiDZthPjd72Qt778vf7dr\nax1bw2+sCX6Hn9wV8stF5zXFSLeN4ZGGCoABAAYgShzA8y8vxpCXgrBSGYeMetNRojxgt1obiFtI\nE6vxomYhrsAPC+rNx1blEYix9QDjk1laKAwQO/Cs+F+8FpaMRfVexNvKC6hV7PtP8XvNJszTnELv\nL9fj1l8t+2xtLfDuu3YtxyX0DrP77wc6dgR++EH+ue5B0tKADRvMX76luIQttUNxWInEtHp7UBXd\nyMYKAgBMBjAZIeIHjHtrBZ58sy9OozM+qfcHrFOeQKXi2K6gJuIakkUOBohs9BfZaIVL+Ebpj6+V\ngdihLMCZXzoByyAfDucPoP9vDwBCIBhFSFi7E4mf5iFeLII/fsEuJQE7lSR8qyThIKLs/ntu4iJ8\noxRxbypboLKyEikpKTh48CAAYOrUqRgyZAhSU1O1y2RkZKCmpgYzZswAAAQHB6PIwLBCRVEwZ85c\n7fPk5BQkJ6cY3vB//oN67y1BvRUfQ/TtD82MFyBiYq39MYCKCtRb/hHqZbwDERIKzayXIPr2c/x5\nhZIS1P/TBKCqCrXLPwG6dLHPei9ehE9MFGpXr4OIT7BqFT4+nn+Ct7wcCA6WtxfWGjtWnjN65hlX\nlWWVvn2B//7vOidVveP8efj0T4bmiTRo/jrPcf+QNTVQtm9DvU9WQPk6G2LgYGgeGQExeIh9zmNV\nVUHJ3wMl599Qcv8N5eABiN4xEP0HQgwYCBEZ5d797OfOQfl2J5Sduai3Mw/4z88QfR6GSEiCSEqG\n6BUtzyFYKDc3B7m5Odrnr78+X69xbZKtf2ZERkaK3NxcUVxcLLp16yZKS0t13s/Pzxfx8fGirKxM\nrFq1SqSmphpcj1WlXLkizwN06iT74desEaKiwrzP3rgh+yvHjJH9l2PGyD5CZ6utFeKdd4Ro0UKI\nzEzb+wlramQf5yuv2Kc+D1ZeLv9pdWRkCPHssy6pxxZJSfIUhUnl5bJ78rXXHF6TjrIyIT78UIhH\nHxWiWTPZvz5zpuwyPX5c/r7VpbJSdkd98okQM2YIkZws+9FjYoT485/l7+q1a075URymtFSITZuE\nmD5diKgo+fP16yf7VXNyTO8jI6zJTpta/PLbJxcTJ05EdXU1pk2bhmnTpmHZMvm3Vnp6OgDg5Zdf\nxtq1axEQEICVK1ciJCREbz2Kolj+rXVbTY2cgGz5cmDXLiAyUo7gCA0FWraUZ++vXwcuXwaOHwcO\nHQK++w4IC5NDLceMAe4aieQSx47J1mirVnLkRdu2lq9DCGDGDDmSITtbNttVrKJCjk6sqLjrxf37\n5SinwkIXVWWdxEQ5qCwxsY6FqquBQYPkqJ2333bdn2xVVfL6lvx8eYvPH36QFxG2aCFHBjVpIlu6\n16/LUTaXLsnaO3WSv5NRUfIRFwf4+bnmZ3CGigo5U25eHpCbC3z/vbzQMDlZPvr0kfvKBGuy0+bg\ntxebgv9u16/LnXn4sAzTy5eBq1flDmzeXA7JiogAEhKAgADbt2dP1dXA3/4mh7MuWSKH0Jn7yysE\n8Je/yCt08/I8a+idg1RUyCyprLzrxepqeSHb2bMeFSrx8XI0Znx8HQvNnCmP+S1b3K8LpLZWDtGu\nrJRhf+uWHErZpIlsdLVo4fl9i7a6elU2SHNz5ePQISA8XA4b79lTNma7dNH7MmDwe4uCAnkxW6tW\nwKJFsvVTlytX5Dw8x4/LX/rWrZ1Tp5urrAQ6dJC7R0dSkvySHDjQJXVZo08fYPFi+V+D1q6V14js\n2+d+DRqyzo0b8i+mgwfl48ABee/oVq3kl2ZQENCsGZS1ay3OTnX3Bbir3r3lt/1HHwHDhslv+z/9\nSQZVs2Z3lrtyRV5FvGABMHiwbCV46ZWG1jDagOzTR7asPCj467yHzJkzwHPPyatxGfreo1EjeVa/\nb987r9XWAiUl8q+m4mLZslm71uJVM/jdlY8PMHGi7I9euxZ4/33g6aflP3RAgOzHKC6Wgf/ppyb6\nANTLYEOoTx/gww+dXostjAa/RiNHKM2YIfuHybvVry+HqgE2DUlm8Lu7hg1l4D/9tDx/UVQkxyfe\nf7+8bNzX19UVui2dSdruFhcnp8/woFsxGi31gw9k3/Cf/+z0mshzMfg9SePG8mQP2aZNG3li98QJ\nOfeKBzAY/OfPA3/9qxzMoPIRXGQZz2juEFnBaIsfkN09e/Y4tR5baDQGzlm8+CLwxz96zJcXuQ8G\nP3mtOkcH3j7B6yH0ZkzNzZXDdufMcVlN5LkY/OTVjLb44+I8Kvh1unpqauQonsWLzbrAh+heDH7y\nWnV29fTsKcdEX73q1JqspRP8K1bIC55GjnRlSeTBGPzkters6rnvPnlh3N69TqvHFtrgv34dmDcP\nePNNXulKVmPwk1er84LGhx+WI2I8gDb4MzKA2Fj5ILISx4CR16qzqweQUzcsXuy0emyh0QA+V38B\n/v53YOdOV5dDHo4tfvJaJntC4uNlV8+tW06pxxYaDRDwv2/K+1Bz+CbZiC1+8mp1tvibNwe6dpUT\nmz38sNNqssb91ZfR7NMPgcLDri6FvABb/OS1THb1ALK7Jy/PKfXY4ukrGbgx+PdyriYiGzH4Sd2S\nk90/+K9dw7hrS3FtMufjIftg8JPXMqvFn5Ag79pWU+OUmqzy4YfYfV9f1Hbp5upKyEsw+MlrmRX8\nLVsC7dvLO7a5o5s3gcWLsaTJbE+ZSJQ8AA8louRkOfeNO/rnP4GwMBypH8XgJ7vhoURey6wWPyBP\n8Lpj8NfWyit0Z8/2pFsHkAfgoURey+wZDfr2lcHvbv38GzbIrqikJAY/2RUPJfJqZrX4W7cGOnd2\nr3l7hJD3Up49G1AUBj/ZFQ8l8loWzWE2cCCQne2wWiz21Vfyct3UVAAedZdI8gA8lMhreXTwL1gA\nvPyyNu0Z/GRPPJTI65nV3ZOYKId0Xrni8HpM2rULOHsWePxx7UsMfrInHkrk9cwK/kaN5FTHOTmO\nLse0BQvk/XTvuoE6g5/siYcS0W2DBgHbt7u2hsOHgQMHgD/8QedlgzdbJ7ISg5+8mtlj+QFg8GDg\nX/+y4AMOsHAh8PzzQMOGOi/r3WydyAY8lMirWRT8ERFyLP+xYw6tyahTp+QJ5okT9d5iVw/ZEw8l\notsUBRgxAvjiC9ds/3/+B5g0Cbj/fr23GPxkTzyUyKtZ1OIHZPBnZTmsHqPOnQPWrwemTdN7Swj5\nYB8/2QuDn7yaxcGfnCy7ei5ccFhNBi1eDIwdK6douMft0Gfwk70w+Inudt998iTvli3O22ZZGbB8\nuRzCaQC7ecjeeDiRV7O4xQ84v5//3XeBUaPkfQEMYPCTvSlCuHLs2h2KosBNSiEv4usrL8b19bXg\nQ5WVQMeOQEkJ4O/vqNKkigqgSxc5QVxQkMFFqqrkfeGrqhxbCnkma7KT7Qjyeha3J/z8gAEDgE2b\nHFKPjqVLgWHDjIY+wBY/2Z/Vh9PVq1fx6KOPomPHjnjsscdw7do1g8t17twZERERiIqKQkxMjNWF\nElnDqq4eABgzBli92u716Pj1V2DJEjn1ch141S7Zm9XBn5mZiY4dO+LkyZNo3749PvjgA4PLKYqC\nnJwcHDx4EHvdab5zoroMGyanTjh/3nHbyMyUd/8KCalzMV61S/Zm9eG0d+9ePPPMM/D19cWECROQ\nn59vdFn23ZOrWN3ib9QIeOwxx7X6KyuBt94C5s0zuSi7esjefEwvYlhBQQG6d+8OAOjevbvR1ryi\nKOjXrx8eeOABTJgwASNGjLB2k0QWUxR5B0OLTu7+JvCBZ9F70QT863cz7d7XErZuERqHDMXeoz2A\no3Uv++uvDH6yrzqDf+DAgbhg4EKWv/3tb2a34nft2oW2bdvi2LFjGD58OGJiYtCmTRuDy867q/WT\nkpKClJQUs7ZBZMyf/gRs3Wrlh8XD6HLdFz9+8G9837qf3Wq6v+oShn21FLMH7kepmeePn33Wbpsn\nD5eTk4McG6cPt3o458iRIzFnzhxERUVh//79WLBgATZs2FDnZ2bOnImQkBD88Y9/1C+EwznJHb3/\nvpyjf906+63z+edl/82SJfZbJ6mWU4dzxsbGYvny5bhx4waWL1+OuLg4vWWuX7+Oq1evAgBKS0ux\nbds2DBkyxNpNEjnfU08BO3bIO2LZw/HjwKpVwCuv2Gd9RFawOvgnTZqEM2fOoFu3bjh79iwm/jaV\n7Llz55D62w2iL1y4gMTERERGRmL06NF44YUX0KFDB/tUTuQM998vb4qyaJHt6xICmDoVmDMHaN3a\n9vURWYntpsLvAAAHgklEQVRX7hKZcu4cEBYGnDgBBAZav55Nm4C5c4GDB3Vuq0hkC2uyk8FPZI5J\nk4BmzeQQTGtUVMgbvXzyCcBBC2RHDH4iRzl/HggPB/LzgeBgyz8/dqycCuK99+xfG6ka5+ohcpS2\nbYEXXgBmzLD8irC1a+UXxptvOqY2Igsx+InMNXOmnLFzxQrzP1NYCDz3HLBmDdCkiaMqI7IIzzAR\nmcvXV07h0LcvEB0tu37qcvo0MHy4nG+/Vy/n1EhkBrb4iSwRFgZkZABDhwJFRcaX+/FHoF8/+VfC\nmDHOq4/IDAx+IkuNHi2HZcbHA59/rtvnX1srb6OYkCCnWzZw83QiV+OoHiJr5eYCkyfLCdzi44Fb\nt4BvvpG3UMzIAB56yNUVkgpwOCeRs2k0wJ49wOHDQIMGQGys6b5/Ijti8BMRqQzH8RMRkUkMfiIi\nlWHwExGpDIOfiEhlGPxERCrD4CciUhkGPxGRyjD4iYhUhsFPRKQyDH4iIpVh8BMRqQyDn4hIZRj8\nREQqw+AnIlIZBj8Rkcow+ImIVIbBT0SkMgx+IiKVYfATEakMg5+ISGUY/EREKsPgJyJSGQY/EZHK\nMPiJiFSGwU9EpDJWB//69evRo0cP1K9fHwcOHDC6XF5eHkJCQtC1a1dkZGRYuzkiIrITq4M/PDwc\nn332GZKSkupcbvr06Vi2bBl27NiBpUuXoqyszNpNqkZOTo6rS3Ab3Bd3cF/cwX1hG6uDv3v37njw\nwQfrXKayshIAkJSUhE6dOmHQoEHIz8+3dpOqwYP6Du6LO7gv7uC+sI1D+/gLCgrQvXt37fPQ0FDs\n2bPHkZskIiITfOp6c+DAgbhw4YLe62+88QaGDx/usKKIiMiBhI1SUlLE/v37Db5XUVEhIiMjtc+f\ne+45sWXLFoPLBgcHCwB88MEHH3xY8AgODrY4t+ts8ZtLCGHwdT8/PwByZE/Hjh2RnZ2NuXPnGlz2\n1KlT9iiFiIhMsLqP/7PPPkOHDh2wZ88epKamYujQoQCAc+fOITU1VbvcO++8g/T0dAwYMACTJ09G\nYGCg7VUTEZHVFGGsuU5ERF7JqVfumnMx1+zZsxEUFIRevXrh+PHjzizPqUzti1WrVqFnz57o2bMn\nxowZgxMnTrigSucw9yK/goIC+Pj4YNOmTU6szrnM2RcFBQXo3bs3QkJCkJKS4twCncjUvrhx4wae\nfvppREVFITk5GV988YULqnS8CRMmoHXr1ggPDze6jMW5afFZARtERkaK3NxcUVJSIrp16yZKS0t1\n3s/Pzxfx8fHi8uXLYvXq1SI1NdWZ5TmVqX2xe/duUVFRIYQQYsWKFeKpp55yRZlOYWpfCCFETU2N\n6Nu3r0hNTRUbNmxwQZXOYWpfaDQaERYWJrKzs4UQwuC+8ham9kVmZqaYNGmSEEKIkpISERQUJDQa\njStKdai8vDxx4MABERYWZvB9a3LTaS1+cy7mys/Px6hRoxAQEIC0tDQcO3bMWeU5lTn7ok+fPtqT\n46mpqcjNzXV6nc5g7kV+GRkZGDVqFFq2bOnsEp3GnH2xb98+REREYMCAAQDgtefMzNkXfn5+uHr1\nKqqrq1FeXo7GjRtDURRXlOtQiYmJ8Pf3N/q+NbnptOA352KuvXv3IjQ0VPu8ZcuWKCoqclaJTmPp\nhW0ffvih1143Yc6+OHv2LL744gtMmjQJALzylxswb19s27YNiqIgMTERw4cPx7Zt25xdplOYsy/S\n0tJQW1uLwMBAJCQkYNWqVc4u0y1Yk5t2Gc5pL0IIvaGh3vpLbq4dO3Zg5cqV2L17t6tLcZnnn38e\nCxcuhKIoBo8RNamqqsKhQ4ewY8cOXL9+HQMHDsTRo0fRqFEjV5fmdO+99x58fHxw/vx5FBYWIjU1\nFadPn0a9euqadNia3HTaHurdu7fOSYfvv/8ecXFxOsvExsbihx9+0D4vLS1FUFCQs0p0GnP2BQAc\nOXIEEydORFZWFpo3b+7MEp3GnH2xf/9+jB49Gg888AA2btyIyZMnIysry9mlOpw5+6JPnz4YOnQo\n2rRpg6CgIERHRyMvL8/ZpTqcOfsiLy8PTz75JBo3bozY2Fi0a9fOqwdBGGNNbjot+O++mKukpATZ\n2dmIjY3VWSY2NhYbN27E5cuXsXr1aoSEhDirPKcyZ1+cOXMGI0eOxKpVq9ClSxdXlOkU5uyLn376\nCcXFxSguLsaoUaOQmZmJESNGuKJchzJnX8TFxSE3NxfXr19HeXk5Dh48iPj4eFeU61Dm7Iv+/ftj\n8+bN0Gg0+Omnn1BeXq7TPaQW1uSmU7t6bl/MVV1djWnTpiEwMBDLli0DAKSnpyMmJgYJCQmIjo5G\nQEAAVq5c6czynMrUvnj11VdRXl6OiRMnAgAaNGiAvXv3urJkhzG1L9TE1L5o0aIFxo8fj+joaLRs\n2RKvvvoqmjZt6uKqHcPUvhg9ejR++OEH7b549913XVyxY6SlpSE3NxdlZWXo0KED5s+fj+rqagDW\n5yYv4CIiUhl1nQUhIiIGPxGR2jD4iYhUhsFPRKQyDH4iIpVh8BMRqQyDn4hIZRj8REQq8/8Beea0\nop6JcSoAAAAASUVORK5CYII=\n",
|
|
"text": [
|
|
"<matplotlib.figure.Figure at 0x10f6877d0>"
|
|
]
|
|
}
|
|
],
|
|
"prompt_number": 16
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"collapsed": false,
|
|
"input": [
|
|
"x = np.array([1,2,3])\n",
|
|
"l = 1*np.ones(x.shape)\n",
|
|
"u = 2*np.ones(x.shape)\n",
|
|
"\n",
|
|
"activeSet = np.logical_or(x ==l, x == u)\n",
|
|
"\n",
|
|
"import scipy.sparse as sp\n",
|
|
"shape = (x.size, np.sum(activeSet))\n",
|
|
"print shape\n",
|
|
"print (np.arange(shape[1]), np.where(activeSet)[0])\n",
|
|
"print (np.ones(shape[1]), (np.arange(shape[1]), np.where(activeSet)[0]))\n",
|
|
"rI = sp.csr_matrix((np.ones(shape[1]), (np.arange(shape[1]), np.where(activeSet)[0])), shape=shape)\n",
|
|
"print rI.todense()\n",
|
|
"\n"
|
|
],
|
|
"language": "python",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"output_type": "stream",
|
|
"stream": "stdout",
|
|
"text": [
|
|
"(3, 2)\n",
|
|
"(array([0, 1]), array([0, 1]))\n",
|
|
"(array([ 1., 1.]), (array([0, 1]), array([0, 1])))\n",
|
|
"[[ 1. 0.]\n",
|
|
" [ 0. 1.]\n",
|
|
" [ 0. 0.]]\n"
|
|
]
|
|
}
|
|
],
|
|
"prompt_number": 61
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"collapsed": false,
|
|
"input": [
|
|
"\n",
|
|
"\n",
|
|
"class ProjectedGradient(Minimize):\n",
|
|
" name = 'InexactGaussNewton'\n",
|
|
"\n",
|
|
" maxIterCG = 10\n",
|
|
" tolCG = 1e-5\n",
|
|
" \n",
|
|
" lower = -0.5\n",
|
|
" upper = 1. # ensure these are vectors that are the same size.\n",
|
|
" \n",
|
|
" explore = False\n",
|
|
"\n",
|
|
" def projection(self, p):\n",
|
|
" \"\"\"Make sure we are feasible.\"\"\"\n",
|
|
" return np.median(np.c_[l,x,u],axis=1)\n",
|
|
" \n",
|
|
" def activeSet(self, x):\n",
|
|
" \"\"\"If we are on a bound\"\"\"\n",
|
|
" return np.logical_or(x == self.lower, x == self.upper)\n",
|
|
" \n",
|
|
" def bindingSet(self, x):\n",
|
|
" \"\"\"\n",
|
|
" If we are on a bound and the negative gradient points away from the feasible set. \n",
|
|
" \n",
|
|
" Optimality condition. (Satisfies Kuhn-Tucker) MoreToraldo91\n",
|
|
" \n",
|
|
" \"\"\"\n",
|
|
" bind_up = np.logical_and(x == self.lower, self.g >= 0)\n",
|
|
" bind_low = np.logical_and(x == self.upper, self.g <= 0)\n",
|
|
" return np.logical_or(bind_up, bind_low)\n",
|
|
" \n",
|
|
" def findSearchDirection(self):\n",
|
|
" self.aSet_prev = self.activeSet(self.xc)\n",
|
|
" if self.explorePG or not self.exploreCG:\n",
|
|
" p = -self.g\n",
|
|
" else:\n",
|
|
" Z = I_something\n",
|
|
" def reduceHess(v):\n",
|
|
" # Z is tall and skinny\n",
|
|
" Z.T*self.H(Z*v)\n",
|
|
" \n",
|
|
" p, info = sp.linalg.cg(reduceHess, -Z.T*self.g, tol=self.tolCG, maxiter=self.maxIterCG)\n",
|
|
" p = Z*p # bring up to full size\n",
|
|
" aSet_after = self.activeSet(self.xc+p)\n",
|
|
" return p\n",
|
|
" \n",
|
|
" def doEndIteration(self, xt):\n",
|
|
" aSet = self.activeSet(self.xc)\n",
|
|
" bSet = self.bindingSet(self.xc)\n",
|
|
" # implement 3.8\n",
|
|
" self.f_last - self.f <= self.eta_2 * max_decrease # where max decrease\n",
|
|
" # if true go to CG\n",
|
|
" # don't do too many steps of PG in a row.\n",
|
|
" \n",
|
|
" self.explorePG = not np.all(aSet == self.aSet_prev) # explore proximal gradient\n",
|
|
" self.exploreCG = np.all(aSet == bSet) # explore conjugate gradient"
|
|
],
|
|
"language": "python",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"prompt_number": 64
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"collapsed": false,
|
|
"input": [],
|
|
"language": "python",
|
|
"metadata": {},
|
|
"outputs": []
|
|
}
|
|
],
|
|
"metadata": {}
|
|
}
|
|
]
|
|
} |