Fix bugs after rebase

pymc3/step_methods/hmc/base_hmc.py

@@ -2,33 +2,32 @@
 from pymc3.step_methods import arraystep
 from .quadpotential import quad_potential, QuadPotentialDiagAdapt
 from pymc3.step_methods.hmc import integration
-from pymc3.theanof import inputvars
+from pymc3.theanof import inputvars, floatX
+from pymc3.tuning import guess_scaling
 import numpy as np
-import theano
 
 
 class BaseHMC(arraystep.GradientSharedStep):
     default_blocked = True
 
     def __init__(self, vars=None, scaling=None, step_scale=0.25, is_cov=False,
-                 model=None, blocked=True, use_single_leapfrog=False,
-                 potential=None, integrator="leapfrog", **theano_kwargs):
+                 model=None, blocked=True, potential=None,
+                 integrator="leapfrog", **theano_kwargs):
         """Superclass to implement Hamiltonian/hybrid monte carlo
 
         Parameters
         ----------
         vars : list of theano variables
         scaling : array_like, ndim = {1,2}
-            Scaling for momentum distribution. 1d arrays interpreted matrix diagonal.
+            Scaling for momentum distribution. 1d arrays interpreted matrix
+            diagonal.
         step_scale : float, default=0.25
             Size of steps to take, automatically scaled down by 1/n**(1/4)
         is_cov : bool, default=False
-            Treat scaling as a covariance matrix/vector if True, else treat it as a
-            precision matrix/vector
-        model : pymc3 Model instance.  default=Context model
-        blocked: Boolean, default True
-        use_single_leapfrog: Boolean, will leapfrog steps take a single step at a time.
-            default False.
+            Treat scaling as a covariance matrix/vector if True, else treat
+            it as a precision matrix/vector
+        model : pymc3 Model instance
+        blocked: bool, default=True
         potential : Potential, optional
             An object that represents the Hamiltonian with methods `velocity`,
             `energy`, and `random` methods.

pymc3/step_methods/hmc/nuts.py

@@ -90,7 +90,7 @@ class NUTS(BaseHMC):
     def __init__(self, vars=None, Emax=1000, target_accept=0.8,
                  gamma=0.05, k=0.75, t0=10, adapt_step_size=True,
                  max_treedepth=10, on_error='summary',
-                 adapt_mass_matrix=True, early_max_treedepth=8,
+                 early_max_treedepth=8,
                  **kwargs):
         R"""
         Parameters
@@ -153,7 +153,7 @@ def __init__(self, vars=None, Emax=1000, target_accept=0.8,
         This is usually achieved by setting the `tune` parameter if
         `pm.sample` to the desired number of tuning steps.
         """
-        super(NUTS, self).__init__(vars, use_single_leapfrog=True, **kwargs)
+        super(NUTS, self).__init__(vars, **kwargs)
 
         self.Emax = Emax
 
@@ -168,7 +168,6 @@ def __init__(self, vars=None, Emax=1000, target_accept=0.8,
         self.log_step_size_bar = 0
         self.m = 1
         self.adapt_step_size = adapt_step_size
-        self.adapt_mass_matrix = adapt_mass_matrix
         self.max_treedepth = max_treedepth
         self.early_max_treedepth = early_max_treedepth
 
@@ -218,7 +217,7 @@ def astep(self, q0):
 
         self.m += 1
 
-        if self.tune and self.adapt_mass_matrix:
+        if self.tune:
             self.potential.adapt(q, q_grad)
 
         stats = {

pymc3/step_methods/hmc/quadpotential.py

@@ -128,11 +128,17 @@ def __init__(self, n, initial_mean, initial_diag=None, initial_weight=0,
         self._n_samples = 0
         self.adaptation_window = adaptation_window
 
-    def velocity(self, x):
-        return self._var_theano * x
+    def velocity(self, x, out=None):
+        return np.multiply(self._var, x, out=out)
 
-    def energy(self, x):
-        return 0.5 * x.dot(self._var_theano * x)
+    def energy(self, x, velocity=None):
+        if velocity is not None:
+            return 0.5 * x.dot(velocity)
+        return 0.5 * x.dot(self._var * x)
+
+    def velocity_energy(self, x, v_out):
+        self.velocity(x, out=v_out)
+        return 0.5 * scipy.linalg.blas.ddot(x, v_out)
 
     def random(self):
         vals = floatX(normal(size=self._n))