refactoring kwargs, avoid global q, use class attribute instead

classical/liouville.py

@@ -104,7 +104,7 @@ class LiouvillePowersNR(LiouvillePowers):
         H(q, p) = p^2/2 + V(q) and non-recursive Chebyshev polynomial
         expansion or Newtonian polynomial interpolation """
 
-    def __init__(self, nterms, system, recursion=None, sign=None, **kwargs):
+    def __init__(self, nterms, system, recursion=None, **kwargs):
         """ the kinetic energy must be explicit and H(q, p) = p^2/2 + V(q) """
         from combinatorics import cnk_chebyshev_T2, cnk_newton_1
 
@@ -115,6 +115,7 @@ class LiouvillePowersNR(LiouvillePowers):
         assert np.isclose(self.shift, 0.0)
 
         if recursion == 'chebyshev':
+            sign = kwargs.get('sign')
             assert sign is not None
             signed = True if sign == -1 else False
             vfunc = np.vectorize(cnk_chebyshev_T2)
@@ -128,7 +129,7 @@ class LiouvillePowersNR(LiouvillePowers):
                          momentum=system.p, function=system.p,
                          recursion='taylor', **kwargs)
         self.ilnp = self.Lf
-        self.ilnq = [q] + [e*self.scale for e in self.ilnp[:-1]]
+        self.ilnq = [self.q] + [e*self.scale for e in self.ilnp[:-1]]
         self.ilnqexpr = np.inner(self.ilnq, self.cnk)
         self.ilnpexpr = np.inner(self.ilnp, self.cnk)
         self.get_val_float = self.get_val_float_recurs
@@ -146,8 +147,8 @@ class LiouvillePowersNR(LiouvillePowers):
             and p, evaluates (iL)^n q using the identity (iL)^n q == (iL)^{n-1} p
             and finally computes the values of the Chebyshev/Newton polynomials """
         repl = {self.q: qval, self.p: pval}
-        ilnprec = [float(expr.xreplace(repl)) for expr in self.ilnp]
-        ilnqrec = [float(qval)] + [e*self.scale for e in ilnprec[:-1]]
-        qvec = np.inner(ilnqrec, self.cnk)
-        pvec = np.inner(ilnprec, self.cnk)
+        ilnprec = [pval] + [expr.xreplace(repl) for expr in self.ilnp[1:]]
+        ilnqrec = [qval] + [e*self.scale for e in ilnprec[:-1]]
+        qvec = np.inner(np.array(ilnqrec, dtype=float), self.cnk)
+        pvec = np.inner(np.array(ilnprec, dtype=float), self.cnk)
         return (qvec, pvec)