Use logging instead of print (#1416)

pymc3/__init__.py

@@ -25,3 +25,10 @@
 
 from . import glm
 from .data import *
+
+import logging
+_log = logging.getLogger('pymc3')
+if not logging.root.handlers:
+    _log.setLevel(logging.INFO)
+    handler = logging.StreamHandler()
+    _log.addHandler(handler)

pymc3/distributions/multivariate.py

@@ -8,8 +8,9 @@
 import theano.tensor as tt
 
 from scipy import stats
-from theano.tensor.nlinalg import det, matrix_inverse, trace, eigh
+from theano.tensor.nlinalg import det, matrix_inverse, trace
 
+import pymc3 as pm
 from . import transforms
 from .distribution import Continuous, Discrete, draw_values, generate_samples
 from ..model import Deterministic
@@ -20,6 +21,7 @@
 __all__ = ['MvNormal', 'MvStudentT', 'Dirichlet',
            'Multinomial', 'Wishart', 'WishartBartlett', 'LKJCorr']
 
+
 def get_tau_cov(mu, tau=None, cov=None):
     """
     Find precision and standard deviation
@@ -165,7 +167,6 @@ def logp(self, value):
         mu = self.mu
 
         d = S.shape[0]
-        n = value.shape[0]
 
         X = value - mu
 
@@ -297,7 +298,7 @@ def random(self, point=None, size=None):
     def logp(self, x):
         n = self.n
         p = self.p
-        
+
         if x.ndim==2:
             x_sum = x.sum(axis=0)
             n_sum = n * x.shape[0]
@@ -307,16 +308,16 @@ def logp(self, x):
 
         return bound(
             factln(n_sum) + tt.sum(x_sum * tt.log(p) - factln(x_sum)),
-            tt.all(x >= 0), 
-            tt.all(x <= n), 
+            tt.all(x >= 0),
+            tt.all(x <= n),
             tt.eq(tt.sum(x_sum), n_sum),
             tt.isclose(p.sum(), 1),
             n >= 0)
 
 
 def posdef(AA):
     try:
-        fct = np.linalg.cholesky(AA)
+        np.linalg.cholesky(AA)
         return 1
     except np.linalg.LinAlgError:
         return 0
@@ -347,7 +348,7 @@ def perform(self, node, inputs, outputs):
         try:
             z[0] = np.array(posdef(x), dtype='int8')
         except Exception:
-            print('Failed to check if positive definite', x)
+            pm._log.exception('Failed to check if positive definite', x)
             raise
 
     def infer_shape(self, node, shapes):
@@ -488,9 +489,9 @@ def WishartBartlett(name, S, nu, is_cholesky=False, return_cholesky=False, testv
 
     c = tt.sqrt(ChiSquared('c', nu - np.arange(2, 2 + n_diag), shape=n_diag,
                            testval=diag_testval))
-    print('Added new variable c to model diagonal of Wishart.')
+    pm._log.info('Added new variable c to model diagonal of Wishart.')
     z = Normal('z', 0, 1, shape=n_tril, testval=tril_testval)
-    print('Added new variable z to model off-diagonals of Wishart.')
+    pm._log.info('Added new variable z to model off-diagonals of Wishart.')
     # Construct A matrix
     A = tt.zeros(S.shape, dtype=np.float32)
     A = tt.set_subtensor(A[diag_idx], c)

pymc3/model.py

@@ -3,6 +3,7 @@
 import theano.tensor as tt
 from theano.tensor.var import TensorVariable
 
+import pymc3 as pm
 from .memoize import memoize
 from .theanof import gradient, hessian, inputvars
 from .vartypes import typefilter, discrete_types, continuous_types
@@ -81,7 +82,7 @@ def get_named_nodes(graph):
 
 
 def _get_named_nodes(graph, nodes):
-    if graph.owner == None:
+    if graph.owner is None:
         if graph.name is not None:
             nodes.update({graph.name: graph})
     else:
@@ -287,11 +288,11 @@ def Var(self, name, dist, data=None):
                 var = TransformedRV(name=name, distribution=dist, model=self,
                                     transform=dist.transform)
                 if self.verbose:
-                    print('Applied {transform}-transform to {name}'
-                          ' and added transformed {orig_name} to model.'.format(
-                              transform=dist.transform.name,
-                              name=name,
-                              orig_name='{}_{}_'.format(name, dist.transform.name)))
+                    pm._log.info('Applied {transform}-transform to {name}'
+                                 ' and added transformed {orig_name} to model.'.format(
+                                    transform=dist.transform.name,
+                                    name=name,
+                                    orig_name='{}_{}_'.format(name, dist.transform.name)))
                 self.deterministics.append(var)
                 return var
         elif isinstance(data, dict):

pymc3/plots.py

@@ -1,8 +1,9 @@
 import numpy as np
 from scipy.stats import kde, mode
-from .stats import autocorr, quantiles, hpd
 from numpy.linalg import LinAlgError
 import matplotlib.pyplot as plt
+import pymc3 as pm
+from .stats import quantiles, hpd
 
 __all__ = ['traceplot', 'kdeplot', 'kde2plot',
            'forestplot', 'autocorrplot', 'plot_posterior']
@@ -69,7 +70,7 @@ def traceplot(trace, varnames=None, transform=lambda x: x, figsize=None,
     if ax is None:
         fig, ax = plt.subplots(n, 2, squeeze=False, figsize=figsize)
     elif ax.shape != (n, 2):
-        print('traceplot requires n*2 subplots')
+        pm._log.warning('traceplot requires n*2 subplots')
         return None
 
     for i, v in enumerate(varnames):
@@ -168,14 +169,14 @@ def kde2plot_op(ax, x, y, grid=200):
 
 def kdeplot(data, ax=None):
     if ax is None:
-        f, ax = subplots(1, 1, squeeze=True)
+        f, ax = plt.subplots(1, 1, squeeze=True)
     kdeplot_op(ax, data)
     return ax
 
 
 def kde2plot(x, y, grid=200, ax=None):
     if ax is None:
-        f, ax = subplots(1, 1, squeeze=True)
+        f, ax = plt.subplots(1, 1, squeeze=True)
     kde2plot_op(ax, x, y, grid)
     return ax
 

pymc3/sampling.py

@@ -4,7 +4,7 @@
 from numpy.random import randint, seed
 from numpy import shape, asarray
 
-from . import backends
+import pymc3 as pm
 from .backends.base import merge_traces, BaseTrace, MultiTrace
 from .backends.ndarray import NDArray
 from .model import modelcontext, Point
@@ -69,7 +69,7 @@ def assign_step_methods(model, step=None, methods=(NUTS, HamiltonianMC, Metropol
         if var not in assigned_vars:
             selected = max(methods, key=lambda method: method._competence(var))
             if model.verbose:
-                print('Assigned {0} to {1}'.format(selected.__name__, var))
+                pm._log.info('Assigned {0} to {1}'.format(selected.__name__, var))
             selected_steps[selected].append(var)
 
     # Instantiate all selected step methods
@@ -262,7 +262,7 @@ def _choose_backend(trace, chain, shortcuts=None, **kwds):
         return NDArray(**kwds)
 
     if shortcuts is None:
-        shortcuts = backends._shortcuts
+        shortcuts = pm.backends._shortcuts
 
     try:
         backend = shortcuts[trace]['backend']

pymc3/stats.py

@@ -9,6 +9,7 @@
 
 from scipy.stats.distributions import pareto
 
+import pymc3 as pm
 from .backends import tracetab as ttab
 
 __all__ = ['autocorr', 'autocov', 'dic', 'bpic', 'waic', 'loo', 'hpd', 'quantiles',
@@ -114,7 +115,7 @@ def waic(trace, model=None, n_eff=False):
     model : PyMC Model
         Optional model. Default None, taken from context.
     n_eff: bool
-        if True the effective number parameters will be returned. 
+        if True the effective number parameters will be returned.
         Default False
 
     Returns
@@ -131,8 +132,8 @@ def waic(trace, model=None, n_eff=False):
 
     vars_lpd = np.var(log_py, axis=0)
     if np.any(vars_lpd > 0.4):
-        warnings.warn("""For one or more samples the posterior variance of the 
-        log predictive densities exceeds 0.4. This could be indication of 
+        warnings.warn("""For one or more samples the posterior variance of the
+        log predictive densities exceeds 0.4. This could be indication of
         WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details
         """)
     p_waic = np.sum(vars_lpd)
@@ -155,7 +156,7 @@ def loo(trace, model=None, n_eff=False):
     model : PyMC Model
         Optional model. Default None, taken from context.
     n_eff: bool
-        if True the effective number parameters will be computed and returned. 
+        if True the effective number parameters will be computed and returned.
         Default False
 
     Returns
@@ -189,7 +190,7 @@ def loo(trace, model=None, n_eff=False):
     if np.any(pareto_fit[0] > 0.5):
         warnings.warn("""Estimated shape parameter of Pareto distribution
         is for one or more samples is greater than 0.5. This may indicate
-        that the variance of the Pareto smoothed importance sampling estimate 
+        that the variance of the Pareto smoothed importance sampling estimate
         is very large.""")
 
     # Calculate expected values of the order statistics of the fitted Pareto
@@ -404,7 +405,7 @@ def quantiles(x, qlist=(2.5, 25, 50, 75, 97.5), transform=lambda x: x):
         return dict(zip(qlist, quants))
 
     except IndexError:
-        print("Too few elements for quantile calculation")
+        _log.warning("Too few elements for quantile calculation")
 
 
 def df_summary(trace, varnames=None, stat_funcs=None, extend=False, include_transformed=False,
@@ -438,7 +439,7 @@ def df_summary(trace, varnames=None, stat_funcs=None, extend=False, include_tran
         addition to, rather than in place of, the default statistics.
         This is only meaningful when `stat_funcs` is not None.
     include_transformed : bool
-        Flag for reporting automatically transformed variables in addition to 
+        Flag for reporting automatically transformed variables in addition to
         original variables (defaults to False).
     alpha : float
         The alpha level for generating posterior intervals. Defaults
@@ -545,7 +546,7 @@ def summary(trace, varnames=None, alpha=0.05, start=0, batches=100, roundto=3,
       The number of digits to round posterior statistics.
 
     include_transformed : bool
-      Flag for summarizing automatically transformed variables in addition to 
+      Flag for summarizing automatically transformed variables in addition to
       original variables (defaults to False).
 
     tofile : None or string

pymc3/tuning/starting.py

@@ -6,6 +6,7 @@
 from scipy import optimize
 import numpy as np
 from numpy import isfinite, nan_to_num, logical_not
+import pymc3 as pm
 from ..vartypes import discrete_types, typefilter
 from ..model import modelcontext, Point
 from ..theanof import inputvars
@@ -51,10 +52,10 @@ def find_MAP(start=None, vars=None, fmin=None, return_raw=False,
     kwargs["disp"] = model.verbose > 1
 
     if disc_vars and kwargs["disp"]:
-        print("Warning: vars contains discrete variables. MAP " +
-              "estimates may not be accurate for the default " +
-              "parameters. Defaulting to non-gradient minimization " +
-              "fmin_powell.")
+        pm._log.warning("Warning: vars contains discrete variables. MAP " +
+                        "estimates may not be accurate for the default " +
+                        "parameters. Defaulting to non-gradient minimization " +
+                        "fmin_powell.")
         fmin = optimize.fmin_powell
 
     if fmin is None: