diff --git a/pymbar/__init__.py b/pymbar/__init__.py index 56de3338..94287079 100644 --- a/pymbar/__init__.py +++ b/pymbar/__init__.py @@ -28,6 +28,10 @@ __maintainer__ = "Levi N. Naden, Michael R. Shirts and John D. Chodera" __email__ = "levi.naden@choderalab.org,michael.shirts@colorado.edu,john.chodera@choderalab.org" +from pymbar.logging_utils import setup_logging + +setup_logging() + from pymbar import timeseries, testsystems, confidenceintervals from pymbar.mbar import MBAR from pymbar.other_estimators import bar, bar_zero, exp, exp_gauss diff --git a/pymbar/logging_utils.py b/pymbar/logging_utils.py new file mode 100644 index 00000000..04fef9b0 --- /dev/null +++ b/pymbar/logging_utils.py @@ -0,0 +1,146 @@ +""" +Utilities for our custom logging features +""" +from contextlib import contextmanager +import logging + + +class ForceEmitLogger(logging.Logger): + def log(self, level, msg, *args, **kwargs): + """ + Do what ``logging.Logger`` does, but offer to override all + level thresholds if ``force_emit=True`` is passed to the + ``.log()`` call. + """ + if not isinstance(level, int): + if logging.raiseExceptions: + raise TypeError("level must be an integer") + else: + return + if self.isEnabledFor(level): + self._log(level, msg, args, **kwargs) + elif kwargs.pop("force_emit", None): + with self.temporary_level(level): + self._log(level, msg, args, **kwargs) + + def debug(self, msg, *args, **kwargs): + """ + Do what ``logging.Logger`` does, but offer to override all + level thresholds if ``force_emit=True`` is passed to the + ``.log()`` call. + """ + if self.isEnabledFor(logging.DEBUG): + self._log(logging.DEBUG, msg, args, **kwargs) + elif kwargs.pop("force_emit", None): + with self.temporary_level(logging.DEBUG): + self._log(logging.DEBUG, msg, args, **kwargs) + + def info(self, msg, *args, **kwargs): + """ + Do what ``logging.Logger`` does, but offer to override all + level thresholds if ``force_emit=True`` is passed to the + ``.log()`` call. + """ + if self.isEnabledFor(logging.INFO): + self._log(logging.INFO, msg, args, **kwargs) + elif kwargs.pop("force_emit", None): + with self.temporary_level(logging.INFO): + self._log(logging.INFO, msg, args, **kwargs) + + def warning(self, msg, *args, **kwargs): + """ + Do what ``logging.Logger`` does, but offer to override all + level thresholds if ``force_emit=True`` is passed to the + ``.log()`` call. + """ + if self.isEnabledFor(logging.WARNING): + self._log(logging.WARNING, msg, args, **kwargs) + elif kwargs.pop("force_emit", None): + with self.temporary_level(logging.INFO): + self._log(logging.INFO, msg, args, **kwargs) + + def error(self, msg, *args, **kwargs): + """ + Do what ``logging.Logger`` does, but offer to override all + level thresholds if ``force_emit=True`` is passed to the + ``.log()`` call. + """ + if self.isEnabledFor(logging.ERROR): + self._log(logging.ERROR, msg, args, **kwargs) + elif kwargs.pop("force_emit", None): + with self.temporary_level(logging.ERROR): + self._log(logging.ERROR, msg, args, **kwargs) + + def critical(self, msg, *args, **kwargs): + """ + Do what ``logging.Logger`` does, but offer to override all + level thresholds if ``force_emit=True`` is passed to the + ``.log()`` call. + """ + if self.isEnabledFor(logging.CRITICAL): + self._log(logging.CRITICAL, msg, args, **kwargs) + elif kwargs.pop("force_emit", None): + with self.temporary_level(logging.CRITICAL): + self._log(logging.CRITICAL, msg, args, **kwargs) + + @contextmanager + def temporary_level(self, level=logging.INFO): + """ + Context manager to change the logging level for a block. + """ + logger_level = self.getEffectiveLevel() + handler_levels = [h.level for h in self.handlers] + self.setLevel(level) + try: + yield + finally: + self.setLevel(logger_level) + for handler, handler_level in zip(self.handlers, handler_levels): + handler.level = handler_level + + +class PerLevelFormatter(logging.Formatter): + """ + Adapted from https://stackoverflow.com/a/14859558 + """ + + FORMATS = { + logging.ERROR: "ERROR! %(message)s", + logging.WARNING: "WARNING: %(message)s", + logging.INFO: "%(message)s", + logging.DEBUG: "Debug: %(message)s", + } + + def __init__(self, fmt="%(levelName)d: %(message)s", datefmt=None, style="%", **kwargs): + super().__init__(fmt=fmt, datefmt=datefmt, style=style, **kwargs) + + def format(self, record): + + # Save the original format configured by the user + # when the logger formatter was instantiated + format_orig = self._style._fmt + self._style._fmt = self.FORMATS.get(record.levelno, self._style._fmt) + # Call the original formatter class to do the grunt work + result = super().format(record) + # Restore the original format configured by the user + self._style._fmt = format_orig + + return result + + +def setup_logging(level=logging.WARNING): + """ + Basic configuration for the project, with a custom formatter + """ + logging.setLoggerClass(ForceEmitLogger) + logger = logging.getLogger(__name__.split(".", 1)[0]) + handler = logging.StreamHandler() + formatter = PerLevelFormatter() + handler.setFormatter(formatter) + logger.addHandler(handler) + configure_logging_level(level) + + +def configure_logging_level(level): + logger = logging.getLogger(__name__.split(".", 1)[0]) + logger.setLevel(level) diff --git a/pymbar/mbar.py b/pymbar/mbar.py index 5048d6ad..1ca0e863 100644 --- a/pymbar/mbar.py +++ b/pymbar/mbar.py @@ -216,8 +216,9 @@ def __init__( K, N = np.shape(u_kn) - if verbose: - logger.info("K (total states) = {:d}, total samples = {:d}".format(K, N)) + logger.info( + "K (total states) = {:d}, total samples = {:d}".format(K, N), force_emit=verbose + ) if np.sum(self.N_k) != N: raise ParameterError( @@ -285,9 +286,8 @@ def __init__( logger.warning(dedent(msg[1:])) # Print number of samples from each state. - if self.verbose: - logger.info("N_k = ") - logger.info(self.N_k) + logger.info("N_k = ", force_emit=verbose) + logger.info(self.N_k, force_emit=verbose) # Determine list of k indices for which N_k != 0 self.states_with_samples = np.where(self.N_k != 0)[0] @@ -295,8 +295,9 @@ def __init__( # Number of states with samples. self.K_nonzero = self.states_with_samples.size - if verbose: - logger.info("There are {:d} states with samples.".format(self.K_nonzero)) + logger.info( + "There are {:d} states with samples.".format(self.K_nonzero), force_emit=verbose + ) # Initialize estimate of relative dimensionless free energy of each state to zero. # Note that f_k[0] will be constrained to be zero throughout. @@ -306,8 +307,7 @@ def __init__( # If an initial guess of the relative dimensionless free energies is # specified, start with that. if initial_f_k is not None: - if self.verbose: - logger.info("Initializing f_k with provided initial guess.") + logger.info("Initializing f_k with provided initial guess.", force_emit=verbose) # Cast to np array. initial_f_k = np.array(initial_f_k, dtype=np.float64) # Check shape @@ -317,20 +317,19 @@ def __init__( ) # Initialize f_k with provided guess. self.f_k = initial_f_k - if self.verbose: - logger.info(self.f_k) + logger.info(self.f_k, force_emit=verbose) # Shift all free energies such that f_0 = 0. self.f_k[:] = self.f_k[:] - self.f_k[0] else: # Initialize estimate of relative dimensionless free energies. self._initializeFreeEnergies(verbose, method=initialize) - if self.verbose: - logger.info( - "Initial dimensionless free energies with method {:s}".format(initialize) - ) - logger.info("f_k = ") - logger.info(self.f_k) + logger.info( + "Initial dimensionless free energies with method {:s}".format(initialize), + force_emit=verbose, + ) + logger.info("f_k = ", force_emit=verbose) + logger.info(self.f_k, force_emit=verbose) if solver_protocol is None: solver_protocol = ({"method": None},) @@ -348,13 +347,10 @@ def __init__( self.Log_W_nk = mbar_solvers.mbar_log_W_nk(self.u_kn, self.N_k, self.f_k) # Print final dimensionless free energies. - if self.verbose: - logger.info("Final dimensionless free energies") - logger.info("f_k = ") - logger.info(self.f_k) - - if self.verbose: - logger.info("MBAR initialization complete.") + logger.info("Final dimensionless free energies", force_emit=verbose) + logger.info("f_k = ", force_emit=verbose) + logger.info(self.f_k, force_emit=verbose) + logger.info("MBAR initialization complete.", force_emit=verbose) @property def W_nk(self): @@ -440,15 +436,16 @@ def compute_effective_sample_number(self, verbose=False): for k in range(self.K): w = np.exp(self.Log_W_nk[:, k]) N_eff[k] = 1 / np.sum(w ** 2) - if verbose: - logger.info( - "Effective number of sample in state {:d} is {:10.3f}".format(k, N_eff[k]) - ) - logger.info( - "Efficiency for state {:d} is {:6f}/{:d} = {:10.4f}".format( - k, N_eff[k], len(w), N_eff[k] / len(w) - ) - ) + logger.info( + "Effective number of sample in state {:d} is {:10.3f}".format(k, N_eff[k]), + force_emit=verbose, + ) + logger.info( + "Efficiency for state {:d} is {:6f}/{:d} = {:10.4f}".format( + k, N_eff[k], len(w), N_eff[k] / len(w) + ), + force_emit=verbose, + ) return N_eff @@ -1359,8 +1356,7 @@ def compute_entropy_and_enthalpy( >>> results = mbar.compute_entropy_and_enthalpy() """ - if verbose: - logger.info("Computing average energy and entropy by MBAR.") + logger.info("Computing average energy and entropy by MBAR.", force_emit=verbose) dims = len(np.shape(u_kn)) if dims == 3: @@ -1640,19 +1636,17 @@ def _initializeFreeEnergies(self, verbose=False, method="zeros"): * mean-reduced-potential: the mean reduced potential is used """ - if method == "zeros": # Use zeros for initial free energies. - if verbose: - logger.info("Initializing free energies to zero.") + logger.info("Initializing free energies to zero.", force_emit=verbose) self.f_k[:] = 0.0 elif method == "mean-reduced-potential": # Compute initial guess at free energies from the mean reduced # potential from each state - if verbose: - logger.info( - "Initializing free energies with mean reduced potential for each state." - ) + logger.info( + "Initializing free energies with mean reduced potential for each state.", + force_emit=verbose, + ) means = np.zeros([self.K], float) for k in self.states_with_samples: means[k] = self.u_kn[k, 0 : self.N_k[k]].mean() diff --git a/pymbar/other_estimators.py b/pymbar/other_estimators.py index ac878356..d912b773 100644 --- a/pymbar/other_estimators.py +++ b/pymbar/other_estimators.py @@ -262,7 +262,9 @@ def bar( # this data set is returning NAN -- will likely not work. Return 0, print a warning: # consider returning more information about failure logger.warning( - "BAR is likely to be inaccurate because of poor overlap. Improve the sampling, or decrease the spacing betweeen states. For now, guessing that the free energy difference is 0 with no uncertainty." + "BAR is likely to be inaccurate because of poor overlap. " + "Improve the sampling, or decrease the spacing betweeen states. " + "For now, guessing that the free energy difference is 0 with no uncertainty." ) if compute_uncertainty: result_vals["Delta_f"] = 0.0 @@ -276,8 +278,9 @@ def bar( while FUpperB * FLowerB > 0: # if they have the same sign, they do not bracket. Widen the bracket until they have opposite signs. # There may be a better way to do this, and the above bracket should rarely fail. - if verbose: - logger.info("Initial brackets did not actually bracket, widening them") + logger.info( + "Initial brackets did not actually bracket, widening them", force_emit=verbose + ) FAve = (UpperB + LowerB) / 2 UpperB = UpperB - max(abs(UpperB - FAve), 0.1) LowerB = LowerB + max(abs(LowerB - FAve), 0.1) @@ -303,8 +306,7 @@ def bar( if FNew == 0: # Convergence is achieved. - if verbose: - logger.info("Convergence achieved.") + logger.info("Convergence achieved.", force_emit=verbose) relative_change = 10 ** (-15) break @@ -321,19 +323,16 @@ def bar( # Check for convergence. if DeltaF == 0.0: # The free energy difference appears to be zero -- return. - if verbose: - logger.info("The free energy difference appears to be zero.") + logger.info("The free energy difference appears to be zero.", force_emit=verbose) break if iterated_solution: relative_change = abs((DeltaF - DeltaF_old) / DeltaF) - if verbose: - logger.info("relative_change = {:12.3f}".format(relative_change)) + logger.info("relative_change = {:12.3f}".format(relative_change), force_emit=verbose) if (iteration > 0) and (relative_change < relative_tolerance): # Convergence is achieved. - if verbose: - logger.info("Convergence achieved.") + logger.info("Convergence achieved.", force_emit=verbose) break if method == "false-position" or method == "bisection": @@ -349,18 +348,19 @@ def bar( message = "WARNING: Cannot determine bound on free energy" raise BoundsError(message) - if verbose: - logger.info("iteration {:5d}: DeltaF = {:16.3f}".format(iteration, DeltaF)) + logger.info( + "iteration {:5d}: DeltaF = {:16.3f}".format(iteration, DeltaF), force_emit=verbose + ) # Report convergence, or warn user if not achieved. if iterated_solution: if iteration < maximum_iterations: - if verbose: - logger.info( - "Converged to tolerance of {:e} in {:d} iterations ({:d} function evaluations)".format( - relative_change, iteration, nfunc - ) - ) + logger.info( + "Converged to tolerance of {:e} in {:d} iterations ({:d} function evaluations)".format( + relative_change, iteration, nfunc + ), + force_emit=verbose, + ) else: message = "WARNING: Did not converge to within specified tolerance. max_delta = {:f}, TOLERANCE = {:f}, MAX_ITS = {:d}".format( relative_change, relative_tolerance, maximum_iterations @@ -518,15 +518,13 @@ def bar( ) raise ParameterError(message) - if verbose: - logger.info("DeltaF = {:8.3f} +- {:8.3f}".format(DeltaF, dDeltaF)) + logger.info("DeltaF = {:8.3f} +- {:8.3f}".format(DeltaF, dDeltaF), force_emit=verbose) result_vals["Delta_f"] = DeltaF result_vals["dDelta_f"] = dDeltaF return result_vals else: - if verbose: - logger.info("DeltaF = {:8.3f}".format(DeltaF)) + logger.info("DeltaF = {:8.3f}".format(DeltaF), force_emit=verbose) result_vals["Delta_f"] = DeltaF return result_vals diff --git a/pymbar/pmf.py b/pymbar/pmf.py index 05214a09..b74cb6f1 100644 --- a/pymbar/pmf.py +++ b/pymbar/pmf.py @@ -228,8 +228,7 @@ def __init__(self, u_kn, N_k, verbose=False, mbar_options=None, timings=True, ** # self._random = np.random # self._seed = None - if self.verbose: - logger.info("PMF initialized") + logger.info("PMF initialized", force_emit=verbose) # TODO: see above about not storing np.random # @property @@ -1829,11 +1828,12 @@ def prob(x): if n % sample_every == 0: csamples[:, n // sample_every] = results["c"] logposteriors[n // sample_every] = results["logposterior"] - if n % print_every == 0 and verbose: + if n % print_every == 0: logger.info( "MC Step {:d} of {:d}".format(n, niterations), str(results["logposterior"]), str(bspline.c), + force_emit=verbose, ) # We now have a distribution of samples of parameters sampled according @@ -1843,8 +1843,7 @@ def prob(x): t_mc = 0 g_mc = None - if verbose: - logger.info("Done MC sampling") + logger.info("Done MC sampling", force_emit=verbose) if decorrelate: t_mc, g_mc, Neff = timeseries.detect_equilibration(logposteriors) @@ -1853,28 +1852,34 @@ def prob(x): ) equil_logp = logposteriors[t_mc:] g_mc = timeseries.statistical_inefficiency(equil_logp) - if verbose: - logger.info("Statistical inefficiency of log posterior is {:.3g}".format(g_mc)) + logger.info( + "Statistical inefficiency of log posterior is {:.3g}".format(g_mc), + force_emit=verbose, + ) g_c = np.zeros(len(c)) for nc in range(len(c)): g_c[nc] = timeseries.statistical_inefficiency(csamples[nc, t_mc:]) - if verbose: - logger.info("Time series for spline parameters are: {:s}".format(str(g_c))) + logger.info( + "Time series for spline parameters are: {:s}".format(str(g_c)), force_emit=verbose + ) maxgc = np.max(g_c) meangc = np.mean(g_c) guse = g_mc # doesn't affect the distribution that much indices = timeseries.subsample_correlated_data(equil_logp, g=guse) logposteriors = equil_logp[indices] csamples = (csamples[:, t_mc:])[:, indices] - if verbose: - logger.info("samples after decorrelation: {:d}".format(np.shape(csamples)[1])) + logger.info( + "samples after decorrelation: {:d}".format(np.shape(csamples)[1]), + force_emit=verbose, + ) self.mc_data["samples"] = csamples self.mc_data["logposteriors"] = logposteriors self.mc_data["mc_parameters"] = mc_parameters self.mc_data["acceptance_ratio"] = self.mc_data["naccept"] / niterations - if verbose: - logger.info("Acceptance rate: {:5.3f}".format(self.mc_data["acceptance_ratio"])) + logger.info( + "Acceptance rate: {:5.3f}".format(self.mc_data["acceptance_ratio"]), force_emit=verbose + ) self.mc_data["nequil"] = t_mc # the start of the "equilibrated" data set self.mc_data["g_logposterior"] = g_mc # statistical efficiency of the log posterior self.mc_data["g_parameters"] = g_c # statistical efficiency of the parametere diff --git a/pymbar/timeseries.py b/pymbar/timeseries.py index fd336634..0d77f13b 100644 --- a/pymbar/timeseries.py +++ b/pymbar/timeseries.py @@ -733,17 +733,16 @@ def subsample_correlated_data(A_t, g=None, fast=False, conservative=False, verbo # Compute the statistical inefficiency for the timeseries. if not g: - if verbose: - logger.info("Computing statistical inefficiency...") + logger.info("Computing statistical inefficiency...", force_emit=verbose) g = statistical_inefficiency(A_t, A_t, fast=fast) - if verbose: - logger.info("g = {:f}".format(g)) + logger.info("g = {:f}".format(g), force_emit=verbose) if conservative: # Round g up to determine the stride we can use to pick out regularly-spaced uncorrelated samples. stride = int(math.ceil(g)) - if verbose: - logger.info("conservative subsampling: using stride of {:d}".format(stride)) + logger.info( + "conservative subsampling: using stride of {:d}".format(stride), force_emit=verbose + ) # Assemble list of indices of uncorrelated snapshots. indices = range(0, T, stride) @@ -757,18 +756,19 @@ def subsample_correlated_data(A_t, g=None, fast=False, conservative=False, verbo if (n == 0) or (t != indices[n - 1]): indices.append(t) n += 1 - if verbose: - logger.info("standard subsampling: using average stride of {:f}".format(g)) + logger.info( + "standard subsampling: using average stride of {:f}".format(g), force_emit=verbose + ) # Number of samples in subsampled timeseries. N = len(indices) - if verbose: - logger.info( - "The resulting subsampled set has {:d} samples (original timeseries had {:d}).".format( - N, T - ) - ) + logger.info( + "The resulting subsampled set has {:d} samples (original timeseries had {:d}).".format( + N, T + ), + force_emit=verbose, + ) # Return the list of indices of uncorrelated snapshots. return indices diff --git a/pymbar/utils.py b/pymbar/utils.py index 26222d8f..1f67d646 100644 --- a/pymbar/utils.py +++ b/pymbar/utils.py @@ -25,6 +25,7 @@ from itertools import zip_longest import warnings + import numpy as np import numexpr @@ -390,7 +391,7 @@ def check_w_normalized(W, N_k, tolerance=1.0e-4): # ============================================================================================ # Exception classes -# ============================================================================================= +# ============================================================================================ class ParameterError(Exception):