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