From 991a3e5fbc947fd06e7a36c6dad847d274a1854f Mon Sep 17 00:00:00 2001
From: seogi <sgkang09@gmail.com>
Date: Tue, 3 Jun 2014 16:20:36 -0700
Subject: [PATCH] Working projected gauss newton CG

---
 SimPEG/Optimization.py | 114 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 114 insertions(+)

diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py
index 662e30da..c9747ed8 100644
--- a/SimPEG/Optimization.py
+++ b/SimPEG/Optimization.py
@@ -872,3 +872,117 @@ class NewtonRoot(object):
                 break
 
         return x
+
+class ProjectedGNCG(BFGS, Minimize, Remember):
+
+    def __init__(self, **kwargs):
+        Minimize.__init__(self, **kwargs)
+
+    name = 'Projected GNCG'
+
+    maxIterCG = 5
+    tolCG = 1e-1
+
+    lower = -np.inf
+    upper = np.inf
+
+    def _startup(self, x0):
+        # ensure bound vectors are the same size as the model
+        if type(self.lower) is not np.ndarray:
+            self.lower = np.ones_like(x0)*self.lower
+        if type(self.upper) is not np.ndarray:
+            self.upper = np.ones_like(x0)*self.upper
+
+
+    @Utils.count
+    def projection(self, x):
+        """projection(x)
+
+            Make sure we are feasible.
+
+        """
+        return np.median(np.c_[self.lower,x,self.upper],axis=1)
+
+    @Utils.count
+    def activeSet(self, x):
+        """activeSet(x)
+
+            If we are on a bound
+
+        """
+        return np.logical_or(x <= self.lower, x >= self.upper)
+
+    @property
+    def approxHinv(self):
+        """
+            The approximate Hessian inverse is used to precondition CG.
+
+            Default uses BFGS, with an initial H0 of *bfgsH0*.
+
+            Must be a scipy.sparse.linalg.LinearOperator
+        """
+        _approxHinv = getattr(self,'_approxHinv',None)
+        if _approxHinv is None:
+            M = sp.linalg.LinearOperator( (self.xc.size, self.xc.size), self.bfgs, dtype=self.xc.dtype )
+            return M
+        return _approxHinv
+
+    @approxHinv.setter
+    def approxHinv(self, value):
+        self._approxHinv = value        
+
+    @Utils.timeIt
+    def findSearchDirection(self):
+
+        """
+            findSearchDirection()
+            Finds the search direction based on either CG or steepest descent.
+        """
+        Active = self.activeSet(self.xc)
+        temp = sum((np.ones_like(self.xc.size)-Active))
+        allBoundsAreActive =  temp == self.xc.size
+
+
+        if allBoundsAreActive:
+            Hinv = SolverICG(self.H, M=self.approxHinv, tol=self.tolCG, maxiter=self.maxIterCG)
+            p = Hinv * (-self.g)
+            return p
+        else:
+
+            
+            delx = np.zeros(self.g.size)
+            resid = -(1-Active) * self.g
+
+            # Begin CG iterations.
+            cgiter = 0
+            cgFlag = 0
+            normResid0 = norm(resid)
+
+            while cgFlag == 0:
+
+                cgiter = cgiter + 1
+                dc = (1-Active)*(self.approxHinv*resid)
+                rd = np.dot(resid, dc)
+
+                #  Compute conjugate direction pc.
+                if cgiter == 1:
+                    pc = dc
+                else:
+                    betak = rd / rdlast
+                    pc = dc + betak * pc
+
+                #  Form product Hessian*pc.
+                Hp = self.H*pc
+                Hp = (1-Active)*Hp
+
+                #  Update delx and residual.
+                alphak = rd / np.dot(pc, Hp)
+                delx = delx + alphak*pc
+                resid = resid - alphak*Hp
+                rdlast = rd
+
+                if np.logical_or(norm(resid)/normResid0 <= self.tolCG, cgiter == self.maxIterCG):
+                    cgFlag = 1
+                # End CG Iterations
+
+            return delx