Use SAMS t0 in analysis if available
@@ -1402,9 +1402,8 @@ class MultiStateSamplerAnalyzer(PhaseAnalyzer):
# Correct for potentially-changing log weights
if has_log_weights:
- u_n[iteration] += - np.sum(log_weights[states_slice, iteration])
- # JDC: Don't include changing normalization of expanded ensemble so as not to pull in unconverged log weight region
- #+ (n_replicas * logsumexp(-f_l[:] + log_weights[:, iteration]))
+ u_n[iteration] += - np.sum(log_weights[states_slice, iteration]) \
+ + (n_replicas * logsumexp(-f_l[:] + log_weights[:, iteration]))
logger.debug("Done.")
return u_n
@@ -1963,17 +1962,29 @@ class MultiStateSamplerAnalyzer(PhaseAnalyzer):
n_uncorrelated_iterations : int
"""
u_n = self.get_effective_energy_timeseries(energies, replica_state_indices)
+
+ # For SAMS, if there is a second-stage start time, use only the asymptotically optimal data
+ t0 = 1 # discard minimization frame
+ try:
+ iteration = len(u_n)
+ data = self._reporter.read_online_analysis_data(iteration, 't0')
+ t0 = max(t0, int(data['t0'][0]))
+ logger.debug('t0 found; using initial t0 = {} instead of 1'.format(t0))
+ except Exception as e:
+ # No t0 found
+ pass
+
# Discard equilibration samples.
# TODO: if we include u_n[0] (the energy right after minimization) in the equilibration detection,
# TODO: then number_equilibrated is 0. Find a better way than just discarding first frame.
- i_t, g_i, n_effective_i = utils.get_equilibration_data_per_sample(u_n[1:])
+ i_t, g_i, n_effective_i = utils.get_equilibration_data_per_sample(u_n[t0:])
n_effective_max = n_effective_i.max()
i_max = n_effective_i.argmax()
n_equilibration = i_t[i_max]
g_t = g_i[i_max]
equilibration_data = [n_equilibration, g_t, n_effective_max]
- # Discard also minimization frame.
- equilibration_data[0] += 1
+ # Account for initially discarded frames
+ equilibration_data[0] += t0
self._equilibration_data = tuple(equilibration_data)
logger.debug('Equilibration data: {}'.format(equilibration_data))
return self._equilibration_data
