diff --git a/.gitignore b/.gitignore index bcfb3d73..793b6202 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,4 @@ *.pyc -SimPEG.sublime-project -SimPEG.sublime-workspace -docs/_build/ \ No newline at end of file +*.sublime-project +*.sublime-workspace +docs/_build/ diff --git a/SimPEG/forward/LinearProblem.py b/SimPEG/forward/LinearProblem.py index d30a5b4d..ef92957b 100644 --- a/SimPEG/forward/LinearProblem.py +++ b/SimPEG/forward/LinearProblem.py @@ -13,13 +13,13 @@ class LinearProblem(Problem): return self.G.dot(m) def J(self, m, v, u=None): - return G.dot(v) + return self.G.dot(v) def Jt(self, m, v, u=None): - return G.T.dot(v) + return self.G.T.dot(v) -if __name__ == '__main__': - N = 100 + +def example(N): h = np.ones(N)/N M = TensorMesh([h]) @@ -28,8 +28,6 @@ if __name__ == '__main__': p = -0.25 q = 0.25 - - g = lambda k: np.exp(p*jk[k]*M.vectorCCx)*np.cos(2*np.pi*q*jk[k]*M.vectorCCx) G = np.empty((nk, M.nC)) @@ -38,12 +36,6 @@ if __name__ == '__main__': G[i,:] = g(i) - - plt.figure(1) - for i in range(nk): - plt.plot(G[i,:]) - - m_true = np.zeros(M.nC) m_true[M.vectorCCx > 0.3] = 1. m_true[M.vectorCCx > 0.45] = -0.5 @@ -55,29 +47,29 @@ if __name__ == '__main__': d_obs = d_true + noise - # plt.figure(3) - # plt.plot(d_true,'-o') - # plt.plot(d_obs,'r-o') - - - - - prob = LinearProblem(M) prob.G = G prob.dobs = d_obs prob.std = np.ones_like(d_obs)*0.1 + return prob, m_true + + +if __name__ == '__main__': + + prob, m_true = example(100) + M = prob.mesh + reg = Regularization(M) - opt = InexactGaussNewton(maxIter=20) - inv = Inversion(prob,reg,opt,beta0=1e-4) - m0 = np.zeros_like(m_true) mrec = inv.run(m0) + plt.figure(1) + for i in range(prob.G.shape[0]): + plt.plot(prob.G[i,:]) plt.figure(2) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 6221246f..daf4aebf 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -154,6 +154,7 @@ class Minimize(object): self._iterLS = 0 self._STOP = np.zeros((5,1),dtype=bool) + x0 = self.projection(x0) # ensure that we start of feasible. self.x0 = x0 self.xc = x0 self.xOld = x0 @@ -296,17 +297,18 @@ class Minimize(object): :rtype: numpy.ndarray,bool :return: (xt, passLS) """ - # Armijo linesearch - descent = np.inner(self.g, p) + # Projected Armijo linesearch t = 1 iterLS = 0 while iterLS < self.maxIterLS: xt = self.projection(self.xc + t*p) ft = self.evalFunction(xt, return_g=False, return_H=False) - if ft < self.f + t*self.LSreduction*descent: + descent = np.inner(self.g, xt - self.xc) # this takes into account multiplying by t, but is important for projection. + if ft < self.f + self.LSreduction*descent: break iterLS += 1 t = self.LSshorten*t + # TODO: Check if t is tooo small. self._iterLS = iterLS return xt, iterLS < self.maxIterLS diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index 140b56fa..05b31685 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -8,8 +8,13 @@ import numpy as np import unittest import inspect -happiness = ['The test be workin!', 'You get a gold star!', 'Yay passed!', 'Happy little convergence test!', 'That was easy!', 'Testing is important.', 'You are awesome.', 'Go Test Go!', 'Once upon a time, a happy little test passed.', 'And then everyone was happy.'] -sadness = ['No gold star for you.','Try again soon.','Thankfully, persistence is a great substitute for talent.','It might be easier to call this a feature...','Coffee break?', 'Boooooooo :(', 'Testing is important. Do it again.'] +try: + import getpass + name = getpass.getuser()[0].upper() + getpass.getuser()[1:] +except Exception, e: + name = 'You' +happiness = ['The test be workin!', 'You get a gold star!', 'Yay passed!', 'Happy little convergence test!', 'That was easy!', 'Testing is important.', 'You are awesome.', 'Go Test Go!', 'Once upon a time, a happy little test passed.', 'And then everyone was happy.','Not just a pretty face '+name,'You deserve a pat on the back!','Well done '+name+'!', 'Awesome, '+name+', just awesome.'] +sadness = ['No gold star for you.','Try again soon.','Thankfully, persistence is a great substitute for talent.','It might be easier to call this a feature...','Coffee break?', 'Boooooooo :(', 'Testing is important. Do it again.',"Did you put your clever trousers on today?",'Just think about a dancing dinosaur and life will get better!','You had so much promise '+name+', oh well...', name.upper()+' ERROR!','Get on it '+name+'!', 'You break it, you fix it.'] class OrderTest(unittest.TestCase): """ diff --git a/notebooks/Bound Constraint Problem.ipynb b/notebooks/Bound Constraint Problem.ipynb new file mode 100644 index 00000000..f641fbc9 --- /dev/null +++ b/notebooks/Bound Constraint Problem.ipynb @@ -0,0 +1,218 @@ +{ + "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": [ + "[]" + ] + }, + { + "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": [ + "" + ] + } + ], + "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": {} + } + ] +} \ No newline at end of file