From e0a78842702ab534b62fe165566b5d8c587295c6 Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@sheffield.ac.uk>
Date: Wed, 6 Aug 2014 22:15:28 +0100
Subject: [PATCH] hmc shortcut

---
 GPy/inference/optimization/__init__.py |   2 +-
 GPy/inference/optimization/hmc.py      | 106 +++++++++++++++++++++++--
 2 files changed, 99 insertions(+), 9 deletions(-)

diff --git a/GPy/inference/optimization/__init__.py b/GPy/inference/optimization/__init__.py
index 04e245e3..1590568f 100644
--- a/GPy/inference/optimization/__init__.py
+++ b/GPy/inference/optimization/__init__.py
@@ -1,3 +1,3 @@
 from scg import SCG
 from optimization import *
-from hmc import HMC,Gmodel
+from hmc import HMC,HMC_shortcut
diff --git a/GPy/inference/optimization/hmc.py b/GPy/inference/optimization/hmc.py
index a5b6fe19..93c0e5e3 100644
--- a/GPy/inference/optimization/hmc.py
+++ b/GPy/inference/optimization/hmc.py
@@ -48,13 +48,103 @@ class HMC:
     def _computeH(self,):
         return self.model.objective_function()+self.p.size*np.log(2*np.pi)/2.+np.log(np.linalg.det(self.M))/2.+np.dot(self.p, np.dot(self.Minv,self.p[:,None]))/2.
 
-class Gmodel:
-    def __init__(self,):
-        self.cov = np.array([[1., 0.99],[0.99, 1.]])
-        self.param_array = np.random.rand(2)
+class HMC_shortcut:
+    def __init__(self,model,M=None,stepsize_range=[1e-6, 1e-1],groupsize=5, Hstd_th=[1e-3, 20.]):
+        self.model = model
+        self.stepsize_range = np.log10(stepsize_range)
+        self.p = np.empty_like(model.optimizer_array.copy())
+        self.groupsize = groupsize
+        self.Hstd_th = Hstd_th
+        if M is None:
+            self.M = np.eye(self.p.size)
+        else:
+            self.M = M
+        self.Minv = np.linalg.inv(self.M)
 
-    def grad(self,):
-        return -np.dot(np.linalg.inv(self.cov),self.param_array[:,None])
+    def sample(self, m_iters=1000, hmc_iters=20):
+        thetas = np.empty((m_iters,self.p.size))
+        ps = np.empty((m_iters,self.p.size))
+        for i in xrange(m_iters):
+            # sample a stepsize from the uniform distribution
+            stepsize = np.exp10(np.random.rand()*(self.stepsize_range[1]-self.stepsize_range[0])+self.stepsize_range[0])
+            self.p[:] = np.random.multivariate_normal(np.zeros(self.p.size),self.M)
+            H_old = self._computeH()
+            p_old = self.p.copy()
+            theta_old = self.model.optimizer_array.copy()
+            #Matropolis
+            self._update(hmc_iters, stepsize)
+            H_new = self._computeH()
+
+            if H_old>H_new:
+                k = 1.
+            else:
+                k = np.exp(H_old-H_new)
+            if np.random.rand()<k:
+                thetas[i] = self.model.optimizer_array
+                ps[i] = self.p
+            else:
+                thetas[i] = theta_old
+                ps[i] = p_old
+                self.model.optimizer_array = theta_old
+        return thetas, ps
+
+    def _update(self, hmc_iters, stepsize):
+        theta_buf = np.empty((2*hmc_iters+1,self.model.optimizer_array.size))
+        p_buf = np.empty((2*hmc_iters+1,self.p.size))
+        H_buf = np.empty((2*hmc_iters+1,))
+        # Set initial position
+        theta_buf[hmc_iters] = self.model.optimizer_array
+        p_buf[hmc_iters] = self.p
+        H_buf[hmc_iters] = self._computeH()
+
+        reversal = []
+        pos = 1
+        for i in xrange(hmc_iters):
+            self.p[:] += -self.stepsize/2.*self.model._transform_gradients(self.model.objective_function_gradients())
+            self.model.optimizer_array = self.model.optimizer_array + self.stepsize*np.dot(self.Minv, self.p)
+            self.p[:] += -self.stepsize/2.*self.model._transform_gradients(self.model.objective_function_gradients())
+
+            theta_buf[hmc_iters+pos] = self.model.optimizer_array
+            p_buf[hmc_iters+pos] = self.p
+            H_buf[hmc_iters+pos] = self._computeH()
+
+            if i<self.groupsize:
+                pos += 1
+                continue
+            else:
+                if len(reversal)==0:
+                    Hlist = range(pos,pos-self.groupsize,-1)
+                    if self._testH(H_buf[Hlist]):
+                        pos += 1
+                    else:
+                        # Reverse the trajectory for the 1st time
+                        reversal.add(pos)
+                        pos = -1
+                        self.model.optimizer_array = theta_buf[hmc_iters]
+                        self.p[:] = -p_buf[hmc_iters]
+                else:
+                    Hlist = range(pos,pos+self.groupsize)
+                    if self._testH(H_buf[Hlist]):
+                        pos += -1
+                    else:
+                        # Reverse the trajectory for the 2nd time
+                        r = (hmc_iters - i)%((reversal[0]-pos)*2)
+                        if r>(reversal[0]-pos):
+                            pos_new = 2*reversal[0] - r - pos
+                        else:
+                            pos_new = 2*pos + r - reversal[0]
+                        self.model.optimizer_array = theta_buf[hmc_iters+pos_new]
+                        self.p[:] = p_buf[hmc_iters+pos_new] # the sign of momentum might be wrong!
+                        break
+
+
+    def _testH(self, Hlist):
+        Hstd = np.std(Hlist)
+        if Hstd<self.Hstd_th[0] or Hstd>self.Hstd_th[1]:
+            return False
+        else:
+            return True
+
+    def _computeH(self,):
+        return self.model.objective_function()+self.p.size*np.log(2*np.pi)/2.+np.log(np.linalg.det(self.M))/2.+np.dot(self.p, np.dot(self.Minv,self.p[:,None]))/2.
 
-    def objective_function(self,):
-        return np.dot(self.param_array, np.dot(np.linalg.inv(self.cov),self.param_array[:,None]))/2.