Skip to content
Merged
Show file tree
Hide file tree
Changes from all 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
71 changes: 32 additions & 39 deletions openpathsampling/numerics/wham.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,10 +7,11 @@
# from several files.
import pandas as pd
import numpy as np

import sys
import logging
logger = logging.getLogger(__name__)


class WHAM(object):
"""
Weighted Histogram Analysis Method
Expand Down Expand Up @@ -59,13 +60,13 @@ def __init__(self, tol=1e-10, max_iter=1000000, cutoff=0.05,
@property
def float_format(self):
"""Float output format. Example: 10.8 (default)"""
return lambda x : "{:" + self._float_format + "f}".format(x)
return lambda x: "{:" + self._float_format + "f}".format(x)

@float_format.setter
def float_format(self, value):
self._float_format = value

def load_files(self,fnames): # pragma: no cover
def load_files(self, fnames): # pragma: no cover
"""Load a file or files into the internal structures.

Requires either pandas or something else with pandas-like read_table
Expand All @@ -87,23 +88,22 @@ def load_files(self,fnames): # pragma: no cover
try:
for fname in fnames:
frames.append(pd.read_table(fname, index_col=0, sep=" ",
usecols=[0,1], header=None))
usecols=[0, 1], header=None))
except TypeError:
frames.append(pd.read_table(fnames, index_col=0, sep=" ",
usecols=[0,1], header=None))
usecols=[0, 1], header=None))
fnames = [fnames]
df = pd.concat(frames, axis=1)
df.columns=fnames
df.columns = fnames
return df


def prep_reverse_cumulative(self, df, cutoff=None, tol=None):
"""
Created cleaned dataframe for further analysis.

This version assumes that the input is the result of a reversed
cumulative histogram. That means that it cleans leading input where
the initial
the initial

Parameters
----------
Expand Down Expand Up @@ -131,7 +131,7 @@ def prep_reverse_cumulative(self, df, cutoff=None, tol=None):
hist_max = df.max(axis=0)
raw_cutoff = cutoff*hist_max
cleaned_df = df.apply(
lambda s : [x if x > raw_cutoff[s.name] else 0.0 for x in s]
lambda s: [x if x > raw_cutoff[s.name] else 0.0 for x in s]
)

if self.interfaces is not None:
Expand All @@ -140,10 +140,10 @@ def prep_reverse_cumulative(self, df, cutoff=None, tol=None):
if type(self.interfaces) is not pd.Series:
self.interfaces = pd.Series(data=self.interfaces,
index=df.columns)
greater_almost_equal = lambda a, b : (a >= b or
abs(a - b) < 10e-10)
greater_almost_equal = lambda a, b: (a >= b or
abs(a - b) < 10e-10)
cleaned_df = cleaned_df.apply(
lambda s : [
lambda s: [
(
s.iloc[i]
if greater_almost_equal(s.index[i],
Expand All @@ -155,18 +155,18 @@ def prep_reverse_cumulative(self, df, cutoff=None, tol=None):
)
else:
# clear duplicates of leading values
test_f = lambda val1, val2, val_max : (
test_f = lambda val1, val2, val_max: (
abs(val1 - val2) > tol or abs(val1 - val_max) > tol
)
cleaned_df = cleaned_df.apply(
lambda s : [
s.iloc[i] if test_f(s.iloc[i], s.iloc[i+1], s.max()) else 0.0
lambda s: [
s.iloc[i]
if test_f(s.iloc[i], s.iloc[i+1], s.max()) else 0.0
for i in range(len(s)-1)
] + [s.iloc[-1]]
)
return cleaned_df


def unweighting_tis(self, cleaned_df):
"""
Calculates the "unweighting" values for each histogram.
Expand All @@ -189,13 +189,12 @@ def unweighting_tis(self, cleaned_df):
unweighting values for the input dataframe
"""
unweighting = cleaned_df.copy().applymap(
lambda x : 1.0 if x > 0.0 else 0.0
lambda x: 1.0 if x > 0.0 else 0.0
)
return unweighting


def sum_k_Hk_Q(self, cleaned_df):
"""Sum over histograms for each histogram bin. Length is n_bins
r"""Sum over histograms for each histogram bin. Length is n_bins

Called ``sum_hist`` in other codes, or :math:`\sum_k H_k(Q)` in F&S.
This is the sum over histograms of values for a given histogram bin.
Expand All @@ -212,7 +211,6 @@ def sum_k_Hk_Q(self, cleaned_df):
"""
return cleaned_df.sum(axis=1)


def n_entries(self, cleaned_df):
"""List of counts of entries per histogram. Length is n_hists

Expand All @@ -231,7 +229,6 @@ def n_entries(self, cleaned_df):
"""
return cleaned_df.sum(axis=0)


def weighted_counts_tis(self, unweighting, n_entries):
"""
Product of unweighting and n_entries.
Expand All @@ -249,14 +246,13 @@ def weighted_counts_tis(self, unweighting, n_entries):
pandas.DataFrame
weighted counts matrix, size n_hists by n_dims
"""
weighted_counts = unweighting.apply(lambda s : [x * n_entries[s.name]
for x in s])
weighted_counts = unweighting.apply(lambda s: [x * n_entries[s.name]
for x in s])
return weighted_counts


def generate_lnZ(self, lnZ, unweighting, weighted_counts,
sum_k_Hk_Q, tol=None):
"""
def generate_lnZ(self, lnZ, unweighting, weighted_counts, sum_k_Hk_Q,
tol=None):
r"""
Perform the WHAM iteration to estimate ln(Z_i) for each histogram.

Parameters
Expand Down Expand Up @@ -286,12 +282,11 @@ def generate_lnZ(self, lnZ, unweighting, weighted_counts,
diff = self.tol + 1 # always start above the tolerance
iteration = 0
hists = weighted_counts.columns
bins = weighted_counts.index
# TODO: probably faster if we make wc this a sparse matrix
wc = weighted_counts.values
unw = unweighting.values
lnZ_old = pd.Series(data=lnZ, index=hists)
Z_new = pd.Series(index=hists)
Z_new = pd.Series(index=hists, dtype='float64')
sum_k_Hk_byQ = sum_k_Hk_Q.values
while diff > tol and iteration < self.max_iter:
Z_old = np.exp(lnZ_old)
Expand Down Expand Up @@ -325,7 +320,7 @@ def generate_lnZ(self, lnZ, unweighting, weighted_counts,
# multiplication, which numpy can do very quickly.
#############################################################

w_i = unw[:,i]
w_i = unw[:, i]

# numerator: w_{i,Q} * sum_k_Hk_byQ
numerator_byQ = np.multiply(w_i, sum_k_Hk_byQ)
Expand All @@ -351,7 +346,6 @@ def generate_lnZ(self, lnZ, unweighting, weighted_counts,
self.convergence = (iteration, diff)
return lnZ_old


def get_diff(self, lnZ_old, lnZ_new, iteration):
"""Calculate the difference for this iteration.

Expand All @@ -372,7 +366,7 @@ def get_diff(self, lnZ_old, lnZ_new, iteration):
difference between old and new to use for convergence testing
"""
# get error
diff=0
diff = 0
diff = sum(abs(lnZ_old - lnZ_new))
# check status (mainly for debugging)
if (iteration % self.sample_every == 0): # pragma: no cover
Expand All @@ -382,7 +376,6 @@ def get_diff(self, lnZ_old, lnZ_new, iteration):
logger.debug("lnZnew = " + str(lnZ_new))
return diff


def output_histogram(self, lnZ, sum_k_Hk_Q, weighted_counts):
"""Recombine the data into a joined histogram

Expand All @@ -404,7 +397,8 @@ def output_histogram(self, lnZ, sum_k_Hk_Q, weighted_counts):
"""
Z = np.exp(lnZ)
Z0_over_Zi = Z.iloc[0] / Z
output = pd.Series(index=sum_k_Hk_Q.index, name="WHAM")
output = pd.Series(index=sum_k_Hk_Q.index, name="WHAM",
dtype='float64')
for val in sum_k_Hk_Q.index:
sum_w_over_Z = sum([
weighted_counts.loc[val, hist_i] * Z0_over_Zi[hist_i]
Expand Down Expand Up @@ -438,7 +432,7 @@ def guess_lnZ_crossing_probability(self, cleaned_df):
pandas.Series
initial guess for values of ln(Z_i) for each histogram
"""
df = cleaned_df.apply(lambda s : s/s.max())
df = cleaned_df.apply(lambda s: s/s.max())
# pandas magic, see http://stackoverflow.com/questions/18327624
first_nonzero = df.apply(lambda s: s[s == 1.0].index[0])
df = df.loc[first_nonzero]
Expand All @@ -464,7 +458,7 @@ def check_cleaned_overlaps(self, cleaned_df):
RuntimeError
if the input doesn't have enough of an overlap
"""
for col_idx in range(1,len(cleaned_df.columns)):
for col_idx in range(1, len(cleaned_df.columns)):
col = cleaned_df.columns[col_idx]
prev_col = cleaned_df.columns[col_idx - 1]

Expand Down Expand Up @@ -503,7 +497,7 @@ def wham_bam_histogram(self, input_df):
n_entries = self.n_entries(cleaned)
unweighting = self.unweighting_tis(cleaned)
weighted_counts = self.weighted_counts_tis(unweighting,
n_entries)
n_entries)
try:
lnZ = self.generate_lnZ(guess, unweighting, weighted_counts,
sum_k_Hk_Q)
Expand Down Expand Up @@ -540,7 +534,6 @@ def parsing(parseargs): # pragma: no cover
return opts, args


import sys, os
if __name__ == "__main__": # pragma: no cover
opts, args = parsing(sys.argv[1:])
wham = WHAM(tol=opts.tol, max_iter=opts.max_iter, cutoff=opts.cutoff)
Expand All @@ -558,5 +551,5 @@ def parsing(parseargs): # pragma: no cover
wham_hist = wham.wham_bam_histogram(df)
print(wham_hist.to_string(
header=False,
float_format=lambda x : "{:10.8f}".format(x)
float_format=lambda x: "{:10.8f}".format(x)
))
18 changes: 7 additions & 11 deletions openpathsampling/tests/test_resampling_statistics.py
Original file line number Diff line number Diff line change
@@ -1,22 +1,16 @@
import pandas as pd
from pandas.util.testing import assert_frame_equal
from nose.tools import (assert_equal, assert_not_equal,
assert_almost_equal, raises)
from nose.plugins.skip import Skip, SkipTest
from .test_helpers import (
true_func, assert_equal_array_array, make_1d_traj, data_filename,
assert_items_equal
)

from pandas.testing import assert_frame_equal
from nose.tools import assert_equal
from .test_helpers import assert_items_equal
import openpathsampling as paths
from openpathsampling.high_level.interface_set import GenericVolumeInterfaceSet

import logging
logging.getLogger('openpathsampling.initialization').setLevel(logging.CRITICAL)
logging.getLogger('openpathsampling.ensemble').setLevel(logging.CRITICAL)
logging.getLogger('openpathsampling.storage').setLevel(logging.CRITICAL)
logging.getLogger('openpathsampling.netcdfplus').setLevel(logging.CRITICAL)


class TestResamplingStatistics(object):
# NOTE: we test the mean_df and std_df functions within this
def setup(self):
Expand Down Expand Up @@ -90,6 +84,7 @@ def test_percentile(self):
check_dtype=False
)


class TestBlockResampling(object):
def setup(self):
self.samples = list(range(100))
Expand Down Expand Up @@ -129,7 +124,8 @@ def test_n_blocks_initialization(self):
assert_items_equal(resampler_2.unassigned, list(range(96, 100)))

def test_n_per_block_initialization(self):
resampler = paths.numerics.BlockResampling(self.samples, n_per_block=10)
resampler = paths.numerics.BlockResampling(self.samples,
n_per_block=10)
assert_equal(resampler.n_total_samples, 100)
assert_equal(resampler.n_blocks, 10)
assert_equal(resampler.n_per_block, 10)
Expand Down
Loading