From e17e539bceee9f39caff5dd2260883facd7e2e1e Mon Sep 17 00:00:00 2001
From: Niu <mniu@wpga-pmurmac.lan>
Date: Sat, 2 Aug 2014 14:06:24 +0100
Subject: [PATCH 1/7] initial implementation of hmc

---
 GPy/inference/optimization/__init__.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/GPy/inference/optimization/__init__.py b/GPy/inference/optimization/__init__.py
index 1a8f043b..04e245e3 100644
--- a/GPy/inference/optimization/__init__.py
+++ b/GPy/inference/optimization/__init__.py
@@ -1,2 +1,3 @@
 from scg import SCG
 from optimization import *
+from hmc import HMC,Gmodel

From 18a1513edb6b3ad9fe23d7e039eb755d19a764cb Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@sheffield.ac.uk>
Date: Tue, 5 Aug 2014 16:48:35 +0100
Subject: [PATCH 2/7] add hmc.py

---
 GPy/inference/optimization/hmc.py | 63 +++++++++++++++++++++++++++++++
 1 file changed, 63 insertions(+)
 create mode 100644 GPy/inference/optimization/hmc.py

diff --git a/GPy/inference/optimization/hmc.py b/GPy/inference/optimization/hmc.py
new file mode 100644
index 00000000..a8bb0405
--- /dev/null
+++ b/GPy/inference/optimization/hmc.py
@@ -0,0 +1,63 @@
+"""HMC implementation"""
+
+import numpy as np
+
+
+class HMC:
+    def __init__(self,model,M=None,stepsize=1e-1):
+        self.model = model
+        self.stepsize = stepsize
+        self.p = np.empty_like(model.optimizer_array.copy())
+        if M is None:
+            self.M = np.eye(self.p.size)
+        else:
+            self.M = M
+        self.Minv = np.linalg.inv(self.M)
+
+    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):
+            #Gibbs
+            self.p[:] = np.random.multivariate_normal(np.ones(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)
+            H_new = self._computeH()
+
+            k = np.exp(H_old-H_new)
+            print k
+            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):
+        for i in xrange(hmc_iters):
+            g = self.p.copy()
+            g[:] = 1e-2
+#            self.p[:] += self.stepsize/2.*self.model.grad()[:,0]#*-self.model._transform_gradients(self.model.objective_function_gradients())
+            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[:,None])[:,0]
+            self.p[:] += self.stepsize/2.*-self.model._transform_gradients(self.model.objective_function_gradients())
+            #self.model.optimizer_array = self.model.optimizer_array - self.stepsize*self.model._transform_gradients(self.model.objective_function_gradients())
+
+    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.optimizer_array = np.random.rand(2)
+
+    def grad(self,):
+        return -np.dot(np.linalg.inv(self.cov),self.optimizer_array[:,None])
+
+    def objective_function(self,):
+        return np.dot(self.optimizer_array, np.dot(np.linalg.inv(self.cov),self.optimizer_array[:,None]))/2.

From 0a1e7aface8f50f4655bb7f47bca056510e22cea Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@sheffield.ac.uk>
Date: Wed, 6 Aug 2014 11:01:08 +0100
Subject: [PATCH 3/7] debug hmc code

---
 GPy/inference/optimization/hmc.py | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/GPy/inference/optimization/hmc.py b/GPy/inference/optimization/hmc.py
index a8bb0405..783aadc7 100644
--- a/GPy/inference/optimization/hmc.py
+++ b/GPy/inference/optimization/hmc.py
@@ -28,11 +28,14 @@ class HMC:
             H_new = self._computeH()
 
             k = np.exp(H_old-H_new)
+            print H_old,H_new
             print k
-            if np.random.rand()<k:
+            a = np.random.rand()
+            if a<k:
                 thetas[i] = self.model.optimizer_array
                 ps[i] = self.p
             else:
+                print '###'
                 thetas[i] = theta_old
                 ps[i] = p_old
                 self.model.optimizer_array = theta_old
@@ -40,13 +43,10 @@ class HMC:
 
     def _update(self, hmc_iters):
         for i in xrange(hmc_iters):
-            g = self.p.copy()
-            g[:] = 1e-2
-#            self.p[:] += self.stepsize/2.*self.model.grad()[:,0]#*-self.model._transform_gradients(self.model.objective_function_gradients())
-            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[:,None])[:,0]
-            self.p[:] += self.stepsize/2.*-self.model._transform_gradients(self.model.objective_function_gradients())
-            #self.model.optimizer_array = self.model.optimizer_array - self.stepsize*self.model._transform_gradients(self.model.objective_function_gradients())
+            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)
+            print self.model.objective_function()
+            self.p[:] += -self.stepsize/2.*self.model._transform_gradients(self.model.objective_function_gradients())
 
     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.
@@ -54,10 +54,10 @@ class HMC:
 class Gmodel:
     def __init__(self,):
         self.cov = np.array([[1., 0.99],[0.99, 1.]])
-        self.optimizer_array = np.random.rand(2)
+        self.param_array = np.random.rand(2)
 
     def grad(self,):
-        return -np.dot(np.linalg.inv(self.cov),self.optimizer_array[:,None])
+        return -np.dot(np.linalg.inv(self.cov),self.param_array[:,None])
 
     def objective_function(self,):
-        return np.dot(self.optimizer_array, np.dot(np.linalg.inv(self.cov),self.optimizer_array[:,None]))/2.
+        return np.dot(self.param_array, np.dot(np.linalg.inv(self.cov),self.param_array[:,None]))/2.

From 23120c53e2670de8b9c17fa6a75e87cb60be24ef Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@sheffield.ac.uk>
Date: Wed, 6 Aug 2014 11:10:34 +0100
Subject: [PATCH 4/7] remove prints

---
 GPy/inference/optimization/hmc.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/GPy/inference/optimization/hmc.py b/GPy/inference/optimization/hmc.py
index 783aadc7..859877f2 100644
--- a/GPy/inference/optimization/hmc.py
+++ b/GPy/inference/optimization/hmc.py
@@ -28,14 +28,11 @@ class HMC:
             H_new = self._computeH()
 
             k = np.exp(H_old-H_new)
-            print H_old,H_new
-            print k
             a = np.random.rand()
             if a<k:
                 thetas[i] = self.model.optimizer_array
                 ps[i] = self.p
             else:
-                print '###'
                 thetas[i] = theta_old
                 ps[i] = p_old
                 self.model.optimizer_array = theta_old
@@ -45,7 +42,6 @@ class HMC:
         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)
-            print self.model.objective_function()
             self.p[:] += -self.stepsize/2.*self.model._transform_gradients(self.model.objective_function_gradients())
 
     def _computeH(self,):

From d11a2c482c3547096a7ded079357ec3b20a173f0 Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@sheffield.ac.uk>
Date: Wed, 6 Aug 2014 16:33:35 +0100
Subject: [PATCH 5/7] correct the initial distribution of p

---
 GPy/inference/optimization/hmc.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/GPy/inference/optimization/hmc.py b/GPy/inference/optimization/hmc.py
index 859877f2..a5b6fe19 100644
--- a/GPy/inference/optimization/hmc.py
+++ b/GPy/inference/optimization/hmc.py
@@ -18,8 +18,7 @@ class HMC:
         thetas = np.empty((m_iters,self.p.size))
         ps = np.empty((m_iters,self.p.size))
         for i in xrange(m_iters):
-            #Gibbs
-            self.p[:] = np.random.multivariate_normal(np.ones(self.p.size),self.M)
+            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()
@@ -27,9 +26,11 @@ class HMC:
             self._update(hmc_iters)
             H_new = self._computeH()
 
-            k = np.exp(H_old-H_new)
-            a = np.random.rand()
-            if a<k:
+            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:

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 6/7] 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.

From c3bb7b28a137ebb8c006c1b0f605fa663df5f890 Mon Sep 17 00:00:00 2001
From: Zhenwen Dai <z.dai@sheffield.ac.uk>
Date: Thu, 7 Aug 2014 14:28:46 +0100
Subject: [PATCH 7/7] debug HMC shortcut

---
 GPy/core/parameterization/parameter_core.py | 19 ++++++
 GPy/inference/optimization/hmc.py           | 69 ++++++++++++---------
 2 files changed, 57 insertions(+), 31 deletions(-)

diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py
index e359409e..e8ef186f 100644
--- a/GPy/core/parameterization/parameter_core.py
+++ b/GPy/core/parameterization/parameter_core.py
@@ -888,6 +888,25 @@ class Parameterizable(OptimizationHandlable):
             self._param_array_ = np.empty(self.size, dtype=np.float64)
         return self._param_array_
 
+    @property
+    def unfixed_param_array(self):
+        """
+        Array representing the parameters of this class.
+        There is only one copy of all parameters in memory, two during optimization.
+
+        !WARNING!: setting the parameter array MUST always be done in memory:
+        m.param_array[:] = m_copy.param_array
+        """
+        if self.__dict__.get('_param_array_', None) is None:
+            self._param_array_ = np.empty(self.size, dtype=np.float64)
+                    
+        if self.constraints[__fixed__].size !=0:
+            fixes = np.ones(self.size).astype(bool)
+            fixes[self.constraints[__fixed__]] = FIXED
+            return self._param_array_[fixes]
+        else:
+            return self._param_array_
+
     @param_array.setter
     def param_array(self, arr):
         self._param_array_ = arr
diff --git a/GPy/inference/optimization/hmc.py b/GPy/inference/optimization/hmc.py
index 93c0e5e3..8c65cdf0 100644
--- a/GPy/inference/optimization/hmc.py
+++ b/GPy/inference/optimization/hmc.py
@@ -15,13 +15,12 @@ class HMC:
         self.Minv = np.linalg.inv(self.M)
 
     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))
+        params = np.empty((m_iters,self.p.size))
         for i in xrange(m_iters):
             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()
+            params[i] = self.model.unfixed_param_array
             #Matropolis
             self._update(hmc_iters)
             H_new = self._computeH()
@@ -31,13 +30,10 @@ class HMC:
             else:
                 k = np.exp(H_old-H_new)
             if np.random.rand()<k:
-                thetas[i] = self.model.optimizer_array
-                ps[i] = self.p
+                params[i] = self.model.unfixed_param_array
             else:
-                thetas[i] = theta_old
-                ps[i] = p_old
                 self.model.optimizer_array = theta_old
-        return thetas, ps
+        return params
 
     def _update(self, hmc_iters):
         for i in xrange(hmc_iters):
@@ -49,9 +45,9 @@ class HMC:
         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 HMC_shortcut:
-    def __init__(self,model,M=None,stepsize_range=[1e-6, 1e-1],groupsize=5, Hstd_th=[1e-3, 20.]):
+    def __init__(self,model,M=None,stepsize_range=[1e-6, 1e-1],groupsize=5, Hstd_th=[1e-5, 3.]):
         self.model = model
-        self.stepsize_range = np.log10(stepsize_range)
+        self.stepsize_range = np.log(stepsize_range)
         self.p = np.empty_like(model.optimizer_array.copy())
         self.groupsize = groupsize
         self.Hstd_th = Hstd_th
@@ -62,14 +58,13 @@ class HMC_shortcut:
         self.Minv = np.linalg.inv(self.M)
 
     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))
+        params = 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])
+            stepsize = np.exp(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()
+            params[i] = self.model.unfixed_param_array
             theta_old = self.model.optimizer_array.copy()
             #Matropolis
             self._update(hmc_iters, stepsize)
@@ -80,13 +75,10 @@ class HMC_shortcut:
             else:
                 k = np.exp(H_old-H_new)
             if np.random.rand()<k:
-                thetas[i] = self.model.optimizer_array
-                ps[i] = self.p
+                params[i] = self.model.unfixed_param_array
             else:
-                thetas[i] = theta_old
-                ps[i] = p_old
                 self.model.optimizer_array = theta_old
-        return thetas, ps
+        return params
 
     def _update(self, hmc_iters, stepsize):
         theta_buf = np.empty((2*hmc_iters+1,self.model.optimizer_array.size))
@@ -99,31 +91,43 @@ class HMC_shortcut:
 
         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())
+        i=0
+        while i<hmc_iters:
+            self.p[:] += -stepsize/2.*self.model._transform_gradients(self.model.objective_function_gradients())
+            self.model.optimizer_array = self.model.optimizer_array + stepsize*np.dot(self.Minv, self.p)
+            self.p[:] += -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()
+            i+=1
 
             if i<self.groupsize:
                 pos += 1
                 continue
             else:
                 if len(reversal)==0:
-                    Hlist = range(pos,pos-self.groupsize,-1)
+                    Hlist = range(hmc_iters+pos,hmc_iters+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]
+                        reversal.append(pos)
+                        if hmc_iters-i>pos:
+                            pos = -1
+                            i += pos
+                            self.model.optimizer_array = theta_buf[hmc_iters]
+                            self.p[:] = -p_buf[hmc_iters]
+                        else:
+                            pos_new = pos-hmc_iters+i
+                            self.model.optimizer_array = theta_buf[hmc_iters+pos_new]
+                            self.p[:] = -p_buf[hmc_iters+pos_new]
+                            break
                 else:
-                    Hlist = range(pos,pos+self.groupsize)
+                    Hlist = range(hmc_iters+pos,hmc_iters+pos+self.groupsize)
+#                    print Hlist
+#                    print self._testH(H_buf[Hlist])
+
                     if self._testH(H_buf[Hlist]):
                         pos += -1
                     else:
@@ -132,14 +136,17 @@ class HMC_shortcut:
                         if r>(reversal[0]-pos):
                             pos_new = 2*reversal[0] - r - pos
                         else:
-                            pos_new = 2*pos + r - reversal[0]
+                            pos_new = pos + r
                         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!
+#                        print reversal[0],pos,pos_new
+#                        print H_buf
                         break
 
-
     def _testH(self, Hlist):
         Hstd = np.std(Hlist)
+#        print Hlist
+#        print Hstd
         if Hstd<self.Hstd_th[0] or Hstd>self.Hstd_th[1]:
             return False
         else: