Python pymbar.MBAR类代码示例

def test_mbar_computeEntropyAndEnthalpy():

    """Can MBAR calculate f_k, <u_k> and s_k ??"""

    for system_generator in system_generators:
        name, test = system_generator()
        x_n, u_kn, N_k_output, s_n = test.sample(N_k, mode='u_kn')
        eq(N_k, N_k_output)
        mbar = MBAR(u_kn, N_k)
        f_ij, df_ij, u_ij, du_ij, s_ij, ds_ij = mbar.computeEntropyAndEnthalpy(u_kn)

        fa = test.analytical_free_energies()
        ua = test.analytical_observable('potential energy')
        sa = test.analytical_entropies()

        fa_ij = np.array(np.matrix(fa) - np.matrix(fa).transpose())
        ua_ij = np.array(np.matrix(ua) - np.matrix(ua).transpose())
        sa_ij = np.array(np.matrix(sa) - np.matrix(sa).transpose())

        z = convert_to_differences(f_ij,df_ij,fa)
        eq(z / z_scale_factor, np.zeros(np.shape(z)), decimal=0)
        z = convert_to_differences(u_ij,du_ij,ua)
        eq(z / z_scale_factor, np.zeros(np.shape(z)), decimal=0)
        z = convert_to_differences(s_ij,ds_ij,sa)
        eq(z / z_scale_factor, np.zeros(np.shape(z)), decimal=0)

def test_mbar_computePMF():

    """ testing computePMF """

    name, test = generate_ho()
    x_n, u_kn, N_k_output, s_n = test.sample(N_k, mode='u_kn')
    mbar = MBAR(u_kn,N_k)
    #do a 1d PMF of the potential in the 3rd state:
    refstate = 2
    dx = 0.25
    xmin = test.O_k[refstate] - 1
    xmax = test.O_k[refstate] + 1
    within_bounds = (x_n >= xmin) & (x_n < xmax)
    bin_centers = dx*np.arange(np.int(xmin/dx),np.int(xmax/dx)) + dx/2
    bin_n = np.zeros(len(x_n),int)
    bin_n[within_bounds] = 1 + np.floor((x_n[within_bounds]-xmin)/dx)
    # 0 is reserved for samples outside the domain.  We will ignore this state
    range = np.max(bin_n)+1
    results = mbar.computePMF(u_kn[refstate,:], bin_n, range, uncertainties = 'from-specified', pmf_reference = 1)
    f_i = results['f_i']
    df_i = results['df_i']
    f0_i = 0.5*test.K_k[refstate]*(bin_centers-test.O_k[refstate])**2
    f_i, df_i = f_i[2:], df_i[2:] # first state is ignored, second is zero, with zero uncertainty
    normf0_i = f0_i[1:] - f0_i[0] # normalize to first state
    z = (f_i - normf0_i) / df_i
    eq(z / z_scale_factor, np.zeros(len(z)), decimal=0)

    def update_logZ_with_mbar(self):
        """
        Use MBAR to update logZ estimates.
        """
        if not self.ncfile:
            raise Exception("Cannot update logZ using MBAR since no NetCDF file is storing history.")

        if not self.sampler.update_scheme == 'global-jump':
            raise Exception("Only global jump is implemented right now.")

        if not self.ncfile:
            raise Exception("Must have a storage file attached to use MBAR updates")

        # Extract relative energies.
        if self.verbose:
            print('Updating logZ estimate with MBAR...')
        initial_time = time.time()
        from pymbar import MBAR
        #first = int(self.iteration / 2)
        first = 0
        u_kn = np.array(self.ncfile.variables['u_k'][first:,:]).T
        [N_k, bins] = np.histogram(self.ncfile.variables['state_index'][first:], bins=(np.arange(self.sampler.nstates+1) - 0.5))
        mbar = MBAR(u_kn, N_k)
        Deltaf_ij, dDeltaf_ij, Theta_ij = mbar.getFreeEnergyDifferences(compute_uncertainty=True, uncertainty_method='approximate')
        self.logZ[:] = -mbar.f_k[:]
        self.logZ -= self.logZ[0]
        final_time = time.time()
        elapsed_time = final_time - initial_time
        self._timing['MBAR time'] = elapsed_time
        if self.verbose:
            print('MBAR time    %8.3f s' % elapsed_time)

def test_mbar_computeOverlap():

    # tests with identical states, which gives analytical results.

    d = len(N_k)
    even_O_k = 2.0*np.ones(d)
    even_K_k = 0.5*np.ones(d)
    even_N_k = 100*np.ones(d)
    name, test = generate_ho(O_k = even_O_k, K_k = even_K_k)
    x_n, u_kn, N_k_output, s_n = test.sample(even_N_k, mode='u_kn')
    mbar = MBAR(u_kn, even_N_k)

    overlap_scalar, eigenval, O = mbar.computeOverlap()

    reference_matrix = np.matrix((1.0/d)*np.ones([d,d]))
    reference_eigenvalues = np.zeros(d)
    reference_eigenvalues[0] = 1.0
    reference_scalar = np.float64(1.0)

    eq(O, reference_matrix, decimal=precision)
    eq(eigenval, reference_eigenvalues, decimal=precision)
    eq(overlap_scalar, reference_scalar, decimal=precision)

    # test of more straightforward examples
    for system_generator in system_generators:
        name, test = system_generator()
        x_n, u_kn, N_k_output, s_n = test.sample(N_k, mode='u_kn')
        mbar = MBAR(u_kn, N_k)
        overlap_scalar, eigenval, O = mbar.computeOverlap()

        # rows of matrix should sum to one
        sumrows = np.array(np.sum(O,axis=1))
        eq(sumrows, np.ones(np.shape(sumrows)), decimal=precision)
        eq(eigenval[0], np.float64(1.0), decimal=precision)

def test_mbar_getWeights():

    """ testing getWeights """

    for system_generator in system_generators:
        name, test = system_generator()
        x_n, u_kn, N_k_output, s_n = test.sample(N_k, mode='u_kn')
        mbar = MBAR(u_kn, N_k)
        # rows should be equal to zero
        W = mbar.getWeights()
        sumrows = np.sum(W,axis=0)
        eq(sumrows, np.ones(len(sumrows)), decimal=precision)

def test_mbar_computeExpectationsInner():

    """Can MBAR calculate general expectations inner code (note: this just tests completion)"""

    for system_generator in system_generators:
        name, test = system_generator()
        x_n, u_kn, N_k_output, s_n = test.sample(N_k, mode='u_kn')
        eq(N_k, N_k_output)
        mbar = MBAR(u_kn, N_k)
        A_in = np.array([x_n, x_n ** 2, x_n ** 3])
        u_n = u_kn[:2,:]
        state_map = np.array([[0,0],[1,0],[2,0],[2,1]],int)
        [A_i, d2A_ij] = mbar.computeExpectationsInner(A_in, u_n, state_map)

def test_mbar_computeExpectations_position_differences():

    """Can MBAR calculate E(x_n)??"""

    for system_generator in system_generators:
        name, test = system_generator()
        x_n, u_kn, N_k_output, s_n = test.sample(N_k, mode='u_kn')
        eq(N_k, N_k_output)
        mbar = MBAR(u_kn, N_k)
        mu_ij, sigma_ij = mbar.computeExpectations(x_n, output = 'differences')
        mu0 = test.analytical_observable(observable = 'position')
        z = convert_to_differences(mu_ij, sigma_ij, mu0)
        eq(z / z_scale_factor, np.zeros(np.shape(z)), decimal=0)

def test_mbar_computeExpectations_position_averages():

    """Can MBAR calculate E(x_n)??"""

    for system_generator in system_generators:
        name, test = system_generator()
        x_n, u_kn, N_k_output, s_n = test.sample(N_k, mode='u_kn')
        eq(N_k, N_k_output)
        mbar = MBAR(u_kn, N_k)
        mu, sigma = mbar.computeExpectations(x_n)
        mu0 = test.analytical_observable(observable = 'position')

        z = (mu0 - mu) / sigma
        eq(z / z_scale_factor, np.zeros(len(z)), decimal=0)

def test_exponential_mbar_xkn_squared():
    """Exponential Distribution Test: can MBAR calculate E(x_kn^2)"""
    test = exponential_distributions.ExponentialTestCase(rates)
    x_kn, u_kln, N_k_output = test.sample(N_k)

    eq(N_k, N_k_output)

    mbar = MBAR(u_kln, N_k)

    mu, sigma = mbar.computeExpectations(x_kn ** 2)
    mu0 = test.analytical_x_squared()

    z = (mu0 - mu) / sigma
    eq(z / z_scale_factor, np.zeros(len(z)), decimal=0)

def test_exponential_mbar_xkn_squared():
    """Harmonic Oscillators Test: can MBAR calculate E(x_kn^2)??"""
    test = harmonic_oscillators.HarmonicOscillatorsTestCase(O_k, k_k)
    x_kn, u_kln, N_k_output = test.sample(N_k)

    eq(N_k, N_k_output)

    mbar = MBAR(u_kln, N_k)

    mu, sigma = mbar.computeExpectations(x_kn ** 2)
    mu0 = test.analytical_x_squared()

    z = (mu0 - mu) / sigma
    eq(z / z_scale_factor, np.zeros(len(z)), decimal=0)

def test_mbar_computeEffectiveSampleNumber():
    """ testing computeEffectiveSampleNumber """

    for system_generator in system_generators:
        name, test = system_generator()
        x_n, u_kn, N_k_output, s_n = test.sample(N_k, mode='u_kn')
        eq(N_k, N_k_output)
        mbar = MBAR(u_kn, N_k)

        # one mathematical effective sample numbers should be between N_k and sum_k N_k
        N_eff = mbar.computeEffectiveSampleNumber()
        sumN = np.sum(N_k)
        assert all(N_eff > N_k)
        assert all(N_eff < sumN)

def test_exponential_mbar_free_energies():
    """Exponential Distribution Test: can MBAR calculate correct free energy differences?"""
    test = exponential_distributions.ExponentialTestCase(rates)
    x_kn, u_kln, N_k_output = test.sample(N_k, mode='u_kln')
    eq(N_k, N_k_output)

    mbar = MBAR(u_kln, N_k)
    fe, fe_sigma = mbar.getFreeEnergyDifferences()
    fe, fe_sigma = fe[0,1:], fe_sigma[0,1:]

    fe0 = test.analytical_free_energies()
    fe0 = fe0[1:] - fe0[0]

    z = (fe - fe0) / fe_sigma
    eq(z / z_scale_factor, np.zeros(len(z)), decimal=0)

def test_mbar_computeExpectations_potential():

    """Can MBAR calculate E(u_kn)??"""

    for system_generator in system_generators:
        name, test = system_generator()
        x_n, u_kn, N_k_output, s_n = test.sample(N_k, mode='u_kn')
        eq(N_k, N_k_output)
        mbar = MBAR(u_kn, N_k)
        mu, sigma = mbar.computeExpectations(u_kn, state_dependent = True)
        mu0 = test.analytical_observable(observable = 'potential energy')
        print(mu)
        print(mu0)
        z = (mu0 - mu) / sigma
        eq(z / z_scale_factor, np.zeros(len(z)), decimal=0)

def test_harmonic_oscillators_mbar_free_energies():
    """Harmonic Oscillators Test: can MBAR calculate correct free energy differences?"""
    test = harmonic_oscillators.HarmonicOscillatorsTestCase(O_k, k_k)
    x_kn, u_kln, N_k_output = test.sample(N_k)

    eq(N_k, N_k_output)

    mbar = MBAR(u_kln, N_k)
    fe, fe_sigma = mbar.getFreeEnergyDifferences()
    fe, fe_sigma = fe[0,1:], fe_sigma[0,1:]

    fe0 = test.analytical_free_energies()
    fe0 = fe0[1:] - fe0[0]

    z = (fe - fe0) / fe_sigma
    eq(z / z_scale_factor, np.zeros(len(z)), decimal=0)

def test_harmonic_oscillators_mbar_free_energies():
    """Harmonic Oscillators Test: can MBAR calculate correct free energy differences?"""
    test = harmonic_oscillators.HarmonicOscillatorsTestCase(O_k, k_k)
    x_n, u_kn, origin = test.sample(N_k)
    u_ijn, N_k_output = convert_ukn_to_uijn(u_kn)
    
    eq(N_k, N_k_output.values)

    mbar = MBAR(u_ijn.values, N_k)
    fe, fe_sigma = mbar.getFreeEnergyDifferences()
    fe, fe_sigma = fe[0], fe_sigma[0]

    fe0 = test.analytical_free_energies()

    z = (fe - fe0) / fe_sigma
    z = z[1:]  # First component is undetermined.
    eq(z / z_scale_factor, np.zeros(len(z)), decimal=0)