quadpotential

examples/bfgstest.py

@@ -4,16 +4,18 @@
 random.seed(1)
 seterr(invalid = 'raise')
 
-n = 20
+n = 6
 m =2
-c0 = array([[4., .2],
-               [.2, 1.]])[:m, :m]
+c0 = array([[100**2., .3],
+               [.3, 1.5]])[:m, :m]
 p0 = linalg.inv(c0)
 C =linalg.cholesky(c0)
 
 data = dot(C, random.normal(size = (m, n))).T
+#data = array([[-1.0, 0],[0,0], [0., 1.]])
 
-v0 = 1e-4
+
+v0 = 1e-8
 hessgen = HessApproxGen( n ,v0)
 
 for x in data:
@@ -25,6 +27,6 @@
 a = array([1.,0])[:m]
 #print h.S.dot(a)
 #print h.C.dot(a)
-print h.S.dot(array([1,0.]))
-print h.S.dot(array([0,1.]))
-print C.dot(a)
+
+print h.Hdot(array([1,0.]))
+print h.Hdot(array([0,1.]))

mcex/step_methods/__init__.py

@@ -2,7 +2,6 @@
 from hmc import *
 from fast_hmc import *
 from hmc_lbfgs import *
-from hmc_vlbfgs import *
 from velocity_hmc import *
 from split_hmc import *
 from metropolis import *    

mcex/step_methods/hmc.py

@@ -4,10 +4,7 @@
 @author: johnsalvatier
 ''' 
 from numpy import floor
-from numpy.linalg import solve
-from scipy.linalg import cholesky, cho_solve
-
-
+from quadpotential import *
 from utils import *
 from ..core import * 
 
@@ -59,43 +56,4 @@ def step(state, q0, logp, dlogp):
         return state, metrop_select(mr, q, q0)
 
     return array_step(step, vars, [logp_dict, dlogp_dict])
-
-
-
-def quad_potential(C, is_cov):
-    if is_cov:
-        return QuadPotential(C)
-    else :
-        return QuadPotential_Inv(C) 
-
-class QuadPotential_Inv(object):
-    def __init__(self, A):
-        self.L = cholesky(A, lower = True)
-
-    def velocity(self, x ):
-        return cho_solve((self.L, True), x)
-
-    def random(self):
-        n = normal(size = self.L.shape[0])
-        return dot(self.L, n)
-
-    def energy(self, x):
-        L1x = solve(self.L, x)
-        return .5 * dot(L1x.T, L1x)
-
-
-class QuadPotential(object):
-    def __init__(self, A):
-        self.A = A
-        self.L = cholesky(A, lower = True)
-
-    def velocity(self, x):
-        return dot(self.A, x)
-
-    def random(self):
-        n = normal(size = self.L.shape[0])
-        return solve(self.L.T, n)
-
-    def energy(self, x):
-        return .5 * dot(x, dot(self.A, x))
 

mcex/step_methods/hmc_lbfgs.py

@@ -4,10 +4,6 @@
 @author: johnsalvatier
 ''' 
 from numpy import floor
-from numpy.linalg import solve
-from scipy.linalg import cholesky, cho_solve
-
-
 from utils import *
 from ..core import * 
 from lbfgs import *
@@ -41,21 +37,21 @@ def step(state, q0, logp, dlogp):
         z = np.random.normal(size = n)
         p = p0 = state.hess.C.dot(z)
         #use the leapfrog method
-        p = p - (e/2) * -dlogp(q) # half momentum update
+        p = p - (e/2) * dlogp(q) # half momentum update
 
         for i in range(nstep): 
             #alternate full variable and momentum updates
-            q = q + e * state.hess.Bdot(p)
+            q = q + -e * state.hess.Hdot(p)
             if i != nstep - 1:
-                p = p - e * -dlogp(q)
+                p = p - e *dlogp(q)
 
-        p = p - (e/2) * -dlogp(q)  # do a half step momentum update to finish off
+        p = p - (e/2) * dlogp(q)  # do a half step momentum update to finish off
 
         p = -p 
 
         def energy(d):
-            Cd = state.hess.C.dot(d)
-            return .5 * dot(Cd, Cd)
+            Sd = state.hess.S.Tdot(d)
+            return .5 * dot(Sd, Sd)
 
 
         mr = (-logp(q0)) + energy(p0) - ((-logp(q))  + energy(p))

mcex/step_methods/lbfgs.py

@@ -20,7 +20,7 @@ def update(self, x, lp, dlp):
         self.grads.put(dlp)
 
         d = zip(self.lps, self.xs, self.grads)
-        d.sort(key = lambda a: a[0])
+        #d.sort(key = lambda a: a[0])
 
         lbfgs = LBFGS(self.S0)
         if len(d) > 0:
@@ -62,10 +62,10 @@ def update(self, a, b):
 
 class LBFGS(object):
     def __init__(self, Var0):
-        C0 = sqrt(Var0)
+        S0 = sqrt(Var0)
 
-        self.C = I_abT_Product(C0)
-        self.S = I_abT_Product(1./C0)
+        self.C = I_abT_Product(1./S0)
+        self.S = I_abT_Product(S0)
 
 
     def update(self, s, y, sy):

mcex/step_methods/velocity_hmc.py

@@ -4,11 +4,11 @@
 @author: johnsalvatier
 ''' 
 from numpy import floor
-from numpy.linalg import solve
-from scipy.linalg import cholesky, cho_solve
 import random
 from __builtin__ import sum as sumb
 
+from quadpotential import *
+
 from utils import *
 from ..core import * 
 
@@ -88,42 +88,3 @@ def Cdot(g):
 
     return array_step(step, logp_d_dict, vars)
 
-
-
-def quad_potential(C, is_cov):
-    if is_cov:
-        return QuadPotential(C)
-    else :
-        return QuadPotential_Inv(C) 
-
-class QuadPotential_Inv(object):
-    def __init__(self, A):
-        self.L = cholesky(A, lower = True)
-
-    def velocity(self, x ):
-        return cho_solve((self.L, True), x)
-
-    def random(self):
-        n = normal(size = self.L.shape[0])
-        return dot(self.L, n)
-
-    def energy(self, x):
-        L1x = solve(self.L, x)
-        return .5 * dot(L1x.T, L1x)
-
-
-class QuadPotential(object):
-    def __init__(self, A):
-        self.A = A
-        self.L = cholesky(A, lower = True)
-
-    def velocity(self, x):
-        return dot(self.A, x)
-
-    def random(self):
-        n = normal(size = self.L.shape[0])
-        return solve(self.L.T, n)
-
-    def energy(self, x):
-        return .5 * dot(x, dot(self.A, x))
-