Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Configure logging with support for verbose flags #391

Open
wants to merge 5 commits into
base: pymbar4
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions pymbar/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,10 @@
__maintainer__ = "Levi N. Naden, Michael R. Shirts and John D. Chodera"
__email__ = "[email protected],[email protected],[email protected]"

from pymbar.utils import configure_logging

configure_logging()

from pymbar import timeseries, testsystems, confidenceintervals
from pymbar.mbar import MBAR
from pymbar.other_estimators import bar, bar_zero, exp, exp_gauss
Expand Down
68 changes: 40 additions & 28 deletions pymbar/mbar.py
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@
DataError,
logsumexp,
check_w_normalized,
log_level,
)

logger = logging.getLogger(__name__)
Expand Down Expand Up @@ -217,7 +218,8 @@ def __init__(
K, N = np.shape(u_kn)

if verbose:
logger.info("K (total states) = {:d}, total samples = {:d}".format(K, N))
with log_level(__name__, logging.INFO):
logger.info("K (total states) = {:d}, total samples = {:d}".format(K, N))

if np.sum(self.N_k) != N:
raise ParameterError(
Expand Down Expand Up @@ -286,8 +288,9 @@ def __init__(

# Print number of samples from each state.
if self.verbose:
logger.info("N_k = ")
logger.info(self.N_k)
with log_level(__name__, logging.INFO):
logger.info("N_k = ")
logger.info(self.N_k)

# Determine list of k indices for which N_k != 0
self.states_with_samples = np.where(self.N_k != 0)[0]
Expand All @@ -296,7 +299,8 @@ 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))
with log_level(__name__, logging.INFO):
logger.info("There are {:d} states with samples.".format(self.K_nonzero))

# Initialize estimate of relative dimensionless free energy of each state to zero.
# Note that f_k[0] will be constrained to be zero throughout.
Expand All @@ -307,7 +311,8 @@ def __init__(
# specified, start with that.
if initial_f_k is not None:
if self.verbose:
logger.info("Initializing f_k with provided initial guess.")
with log_level(__name__, logging.INFO):
logger.info("Initializing f_k with provided initial guess.")
# Cast to np array.
initial_f_k = np.array(initial_f_k, dtype=np.float64)
# Check shape
Expand All @@ -318,19 +323,21 @@ def __init__(
# Initialize f_k with provided guess.
self.f_k = initial_f_k
if self.verbose:
logger.info(self.f_k)
with log_level(__name__, logging.INFO):
logger.info(self.f_k)
# 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)
with log_level(__name__, logging.INFO):
logger.info(
"Initial dimensionless free energies with method {:s}".format(initialize)
)
logger.info("f_k = ")
logger.info(self.f_k)

if solver_protocol is None:
solver_protocol = ({"method": None},)
Expand All @@ -349,12 +356,14 @@ def __init__(

# Print final dimensionless free energies.
if self.verbose:
logger.info("Final dimensionless free energies")
logger.info("f_k = ")
logger.info(self.f_k)
with log_level(__name__, logging.INFO):
logger.info("Final dimensionless free energies")
logger.info("f_k = ")
logger.info(self.f_k)

if self.verbose:
logger.info("MBAR initialization complete.")
with log_level(__name__, logging.INFO):
logger.info("MBAR initialization complete.")

@property
def W_nk(self):
Expand Down Expand Up @@ -441,14 +450,15 @@ def compute_effective_sample_number(self, verbose=False):
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)
with log_level(__name__, logging.INFO):
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)
)
)
)

return N_eff

Expand Down Expand Up @@ -1360,7 +1370,8 @@ def compute_entropy_and_enthalpy(

"""
if verbose:
logger.info("Computing average energy and entropy by MBAR.")
with log_level(__name__, logging.INFO):
logger.info("Computing average energy and entropy by MBAR.")

dims = len(np.shape(u_kn))
if dims == 3:
Expand Down Expand Up @@ -1640,19 +1651,20 @@ 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.")
with log_level(__name__, logging.INFO):
logger.info("Initializing free energies to zero.")
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."
)
with log_level(__name__, logging.INFO):
logger.info(
"Initializing free energies with mean reduced potential for each state."
)
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()
Expand Down
35 changes: 22 additions & 13 deletions pymbar/other_estimators.py
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@
# =============================================================================================
import logging
import numpy as np
from pymbar.utils import ParameterError, ConvergenceError, BoundsError, logsumexp
from pymbar.utils import ParameterError, ConvergenceError, BoundsError, logsumexp, log_level

logger = logging.getLogger(__name__)

Expand Down Expand Up @@ -277,7 +277,8 @@ def bar(
# 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")
with log_level(__name__, logging.INFO):
logger.info("Initial brackets did not actually bracket, widening them")
FAve = (UpperB + LowerB) / 2
UpperB = UpperB - max(abs(UpperB - FAve), 0.1)
LowerB = LowerB + max(abs(LowerB - FAve), 0.1)
Expand All @@ -304,7 +305,8 @@ def bar(
if FNew == 0:
# Convergence is achieved.
if verbose:
logger.info("Convergence achieved.")
with log_level(__name__, logging.INFO):
logger.info("Convergence achieved.")
relative_change = 10 ** (-15)
break

Expand All @@ -322,18 +324,21 @@ def bar(
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.")
with log_level(__name__, logging.INFO):
logger.info("The free energy difference appears to be zero.")
break

if iterated_solution:
relative_change = abs((DeltaF - DeltaF_old) / DeltaF)
if verbose:
logger.info("relative_change = {:12.3f}".format(relative_change))
with log_level(__name__, logging.INFO):
logger.info("relative_change = {:12.3f}".format(relative_change))

if (iteration > 0) and (relative_change < relative_tolerance):
# Convergence is achieved.
if verbose:
logger.info("Convergence achieved.")
with log_level(__name__, logging.INFO):
logger.info("Convergence achieved.")
break

if method == "false-position" or method == "bisection":
Expand All @@ -350,17 +355,19 @@ def bar(
raise BoundsError(message)

if verbose:
logger.info("iteration {:5d}: DeltaF = {:16.3f}".format(iteration, DeltaF))
with log_level(__name__, logging.INFO):
logger.info("iteration {:5d}: DeltaF = {:16.3f}".format(iteration, DeltaF))

# 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
with log_level(__name__, logging.INFO):
logger.info(
"Converged to tolerance of {:e} in {:d} iterations ({:d} function evaluations)".format(
relative_change, iteration, nfunc
)
)
)
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
Expand Down Expand Up @@ -519,14 +526,16 @@ def bar(
raise ParameterError(message)

if verbose:
logger.info("DeltaF = {:8.3f} +- {:8.3f}".format(DeltaF, dDeltaF))
with log_level(__name__, logging.INFO):
logger.info("DeltaF = {:8.3f} +- {:8.3f}".format(DeltaF, dDeltaF))
result_vals["Delta_f"] = DeltaF
result_vals["dDelta_f"] = dDeltaF
return result_vals

else:
if verbose:
logger.info("DeltaF = {:8.3f}".format(DeltaF))
with log_level(__name__, logging.INFO):
logger.info("DeltaF = {:8.3f}".format(DeltaF))
result_vals["Delta_f"] = DeltaF
return result_vals

Expand Down
31 changes: 19 additions & 12 deletions pymbar/pmf.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
import numpy as np
import pymbar

from pymbar.utils import kln_to_kn, ParameterError, DataError, logsumexp
from pymbar.utils import kln_to_kn, ParameterError, DataError, logsumexp, log_level
from pymbar import timeseries

# bunch of imports needed for doing newton optimization of B-splines
Expand Down Expand Up @@ -229,7 +229,8 @@ def __init__(self, u_kn, N_k, verbose=False, mbar_options=None, timings=True, **
# self._seed = None

if self.verbose:
logger.info("PMF initialized")
with log_level(__name__, logging.INFO):
logger.info("PMF initialized")

# TODO: see above about not storing np.random
# @property
Expand Down Expand Up @@ -1830,11 +1831,12 @@ def prob(x):
csamples[:, n // sample_every] = results["c"]
logposteriors[n // sample_every] = results["logposterior"]
if n % print_every == 0 and verbose:
logger.info(
"MC Step {:d} of {:d}".format(n, niterations),
str(results["logposterior"]),
str(bspline.c),
)
with log_level(__name__, logging.INFO):
logger.info(
"MC Step {:d} of {:d}".format(n, niterations),
str(results["logposterior"]),
str(bspline.c),
)

# We now have a distribution of samples of parameters sampled according
# to the posterior.
Expand All @@ -1844,7 +1846,8 @@ def prob(x):
g_mc = None

if verbose:
logger.info("Done MC sampling")
with log_level(__name__, logging.INFO):
logger.info("Done MC sampling")

if decorrelate:
t_mc, g_mc, Neff = timeseries.detect_equilibration(logposteriors)
Expand All @@ -1854,27 +1857,31 @@ 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))
with log_level(__name__, logging.INFO):
logger.info("Statistical inefficiency of log posterior is {:.3g}".format(g_mc))
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)))
with log_level(__name__, logging.INFO):
logger.info("Time series for spline parameters are: {:s}".format(str(g_c)))
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]))
with log_level(__name__, logging.INFO):
logger.info("samples after decorrelation: {:d}".format(np.shape(csamples)[1]))

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"]))
with log_level(__name__, logging.INFO):
logger.info("Acceptance rate: {:5.3f}".format(self.mc_data["acceptance_ratio"]))
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
Expand Down
Loading