Python MBAR.computeExpectations方法代码示例

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computeExpectations [as 别名]
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)

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computeExpectations [as 别名]
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)

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computeExpectations [as 别名]
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)

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computeExpectations [as 别名]
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)

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computeExpectations [as 别名]
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)

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computeExpectations [as 别名]
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_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)
    
    x_kn = convert_xn_to_x_kn(x_n) ** 2.

    x_kn = x_kn.values  # Convert to numpy for MBAR
    x_kn[np.isnan(x_kn)] = 0.0  # Convert nans to 0.0

    mu, sigma = mbar.computeExpectations(x_kn)
    mu0 = test.analytical_x_squared()
    
    z = (mu0 - mu) / sigma
    eq(z / z_scale_factor, np.zeros(len(z)), decimal=0)

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computeExpectations [as 别名]
def test_exponential_mbar_xkn():
    """Exponential Distribution Test: can MBAR calculate E(x_kn)?"""
    test = exponential_distributions.ExponentialTestCase(rates)
    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)
    
    x_kn = convert_xn_to_x_kn(x_n)

    x_kn = x_kn.values  # Convert to numpy for MBAR
    x_kn[np.isnan(x_kn)] = 0.0  # Convert nans to 0.0

    mu, sigma = mbar.computeExpectations(x_kn)
    mu0 = test.analytical_means()
    
    z = (mu0 - mu) / sigma
    eq(z / z_scale_factor, np.zeros(len(z)), decimal=0)

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computeExpectations [as 别名]
  # thermodynamic state
  elif observe == 'potential energy':
    A_kn = U_kln

  # observable for estimation is the position
  elif observe == 'position': 
    A_kn = numpy.zeros([K,N_max], dtype = numpy.float64)
    for k in range(0,K):
      A_kn[k,0:N_k[k]] = x_kn[k,0:N_k[k]]  

  elif observe == 'position^2': 
    A_kn = numpy.zeros([K,N_max], dtype = numpy.float64)
    for k in range(0,K):
      A_kn[k,0:N_k[k]] = x_kn[k,0:N_k[k]]**2   

  (A_k_estimated, dA_k_estimated) = mbar.computeExpectations(A_kn)

  As_k_estimated = numpy.zeros([K],numpy.float64)
  dAs_k_estimated = numpy.zeros([K],numpy.float64)

  # 'standard' expectation averages

  ifzero = numpy.array(N_k != 0)

  for k in range(K):
    if (ifzero[k]):
      if (observe == 'position') or (observe == 'position^2'):
        As_k_estimated[k] = numpy.average(A_kn[k,0:N_k[k]])
        dAs_k_estimated[k]  = numpy.sqrt(numpy.var(A_kn[k,0:N_k[k]])/(N_k[k]-1))
      elif (observe == 'RMS displacement' ) or (observe == 'potential energy'):
        As_k_estimated[k] = numpy.average(A_kn[k,k,0:N_k[k]])

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computeExpectations [as 别名]
    def MBAR_analysis(self, debug = False):
	"""MBAR analysis for populations and BICePs score"""
	# load necessary data first
	self.load_data()

	# Suppose the energies sampled from each simulation are u_kln, where u_kln[k,l,n] is the reduced potential energy
	#   of snapshot n \in 1,...,N_k of simulation k \in 1,...,K evaluated at reduced potential for state l.
	self.K = self.nlambda   # number of thermodynamic ensembles
	# N_k[k] will denote the number of correlated snapshots from state k
	N_k = np.array( [len(self.traj[i]['trajectory']) for i in range(self.nlambda)] )
	nsnaps = N_k.max()
	u_kln = np.zeros( (self.K, self.K, nsnaps) )
	nstates = int(self.states)
	print 'nstates', nstates
	states_kn = np.zeros( (self.K, nsnaps) )

	# Get snapshot energies rescored in the different ensembles
	"""['step', 'E', 'accept', 'state', 'sigma_noe', 'sigma_J', 'sigma_cs', 'sigma_pf''gamma']
	[int(step), float(self.E), int(accept), int(self.state), int(self.sigma_noe_index), int(self.sigma_J_index), int(self.sigma_cs_H_index), int(self.sigma_cs_Ha_index), int(self.sigma_cs_N_index), int(self.sigma_cs_Ca_index), int(self.sigma_pf_index), int(self.gamma_index)]               	"""

	for n in range(nsnaps):

  		for k in range(self.K):
    			for l in range(self.K):
				if debug:
      					print 'step', self.traj[k]['trajectory'][n][0],
      				if k==l:
          				print 'E%d evaluated in model %d'%(k,k), self.traj[k]['trajectory'][n][1],
          				u_kln[k,k,n] = self.traj[k]['trajectory'][n][1]
          			state, sigma_noe_index, sigma_J_index, sigma_cs_H_index, sigma_cs_Ha_index, sigma_cs_N_index, sigma_cs_Ca_index, sigma_pf_index, gamma_index = self.traj[k]['trajectory'][n][3:] 	# IMPORTANT: make sure the order of these parameters is the same as the way they are saved in PosteriorSampler
          			print 'state, sigma_noe_index, sigma_J_index, sigma_cs_H_index, sigma_cs_Ha_index, sigma_cs_N_index, sigma_cs_Ca_index, sigma_pf_index, gamma_index', state, sigma_noe_index, sigma_J_index, sigma_cs_H_index, sigma_cs_Ha_index, sigma_cs_N_index, sigma_cs_Ca_index, sigma_pf_index, gamma_index 
          			states_kn[k,n] = state
          			sigma_noe = self.traj[k]['allowed_sigma_noe'][sigma_noe_index]
          			sigma_J = self.traj[k]['allowed_sigma_J'][sigma_J_index]
          			sigma_cs_H = self.traj[k]['allowed_sigma_cs_H'][sigma_cs_H_index]
          			sigma_cs_Ha = self.traj[k]['allowed_sigma_cs_Ha'][sigma_cs_Ha_index]
          			sigma_cs_N = self.traj[k]['allowed_sigma_cs_N'][sigma_cs_N_index]
          			sigma_cs_Ca = self.traj[k]['allowed_sigma_cs_Ca'][sigma_cs_Ca_index]
          			sigma_pf = self.traj[k]['allowed_sigma_pf'][sigma_pf_index]
          			u_kln[k,l,n] = self.sampler[l].neglogP(0, state, sigma_noe, sigma_J, sigma_cs_H, sigma_cs_Ha, sigma_cs_N, sigma_cs_Ca, sigma_pf, gamma_index)
				if debug:
					print 'E_%d evaluated in model_%d'%(k,l), u_kln[k,l,n]


	# Initialize MBAR with reduced energies u_kln and number of uncorrelated configurations from each state N_k.
 	# u_kln[k,l,n] is the reduced potential energy beta*U_l(x_kn), where U_l(x) is the potential energy function for state l,
	# beta is the inverse temperature, and and x_kn denotes uncorrelated configuration n from state k.
	# N_k[k] is the number of configurations from state k stored in u_knm
	# Note that this step may take some time, as the relative dimensionless free energies f_k are determined at this point.
	mbar = MBAR(u_kln, N_k)

	# Extract dimensionless free energy differences and their statistical uncertainties.
#	(Deltaf_ij, dDeltaf_ij) = mbar.getFreeEnergyDifferences()
	#(Deltaf_ij, dDeltaf_ij, Theta_ij) = mbar.getFreeEnergyDifferences(uncertainty_method='svd-ew')
	(Deltaf_ij, dDeltaf_ij, Theta_ij) = mbar.getFreeEnergyDifferences(uncertainty_method='approximate')
	#print 'Deltaf_ij', Deltaf_ij
	#print 'dDeltaf_ij', dDeltaf_ij
	beta = 1.0 # keep in units kT
	#print 'Unit-bearing (units kT) free energy difference f_1K = f_K - f_1: %f +- %f' % ( (1./beta) * Deltaf_ij[0,K-1], (1./beta) * dDeltaf_ij[0,K-1])
	self.f_df = np.zeros( (self.nlambda, 2) )  # first column is Deltaf_ij[0,:], second column is dDeltaf_ij[0,:]
	self.f_df[:,0] = Deltaf_ij[0,:]
	self.f_df[:,1] = dDeltaf_ij[0,:]

	# Compute the expectation of some observable A(x) at each state i, and associated uncertainty matrix.
        # Here, A_kn[k,n] = A(x_{kn})
        #(A_k, dA_k) = mbar.computeExpectations(A_kn)
        self.P_dP = np.zeros( (nstates, 2*self.K) )  # left columns are P, right columns are dP
	if debug:
        	print 'state\tP\tdP'
    	for i in range(nstates):
        	A_kn = np.where(states_kn==i,1,0)
        	(p_i, dp_i) = mbar.computeExpectations(A_kn, uncertainty_method='approximate')
        	self.P_dP[i,0:self.K] = p_i
        	self.P_dP[i,self.K:2*self.K] = dp_i
		print i
        	for p in p_i: print p,
        	for dp in dp_i: print dp,
		print 
	pops, dpops = self.P_dP[:,0:self.K], self.P_dP[:,self.K:2*self.K]

	# save results
	self.save_MBAR()

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computeExpectations [as 别名]
f_df[:,0] = Deltaf_ij[0,:]
f_df[:,1] = dDeltaf_ij[0,:]
print 'Writing %s...'%args.bayesfactorfile
savetxt(args.bayesfactorfile, f_df)
print '...Done.'


# Compute the expectation of some observable A(x) at each state i, and associated uncertainty matrix.
# Here, A_kn[k,n] = A(x_{kn})
#(A_k, dA_k) = mbar.computeExpectations(A_kn)
P_dP = np.zeros( (nstates, 2*K) )  # left columns are P, right columns are dP
if (1):
    print 'state\tP\tdP'
    for i in range(nstates):
        A_kn = np.where(states_kn==i,1,0)
        (p_i, dp_i) = mbar.computeExpectations(A_kn, uncertainty_method='approximate')
        P_dP[i,0:K] = p_i
        P_dP[i,K:2*K] = dp_i
        print i,
        for p in p_i: print p,
        for dp in dp_i: print dp,
        print

pops, dpops = P_dP[:,0:K], P_dP[:,K:2*K]

print 'Writing %s...'%args.popsfile
savetxt(args.popsfile, P_dP)
print '...Done.'

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computeExpectations [as 别名]
  
  # observable for estimation is the position
  elif observe == 'position':
    state_dependent = False
    A_kn = numpy.zeros([K,N_max], dtype = numpy.float64)
    for k in range(0,K):
      A_kn[k,0:N_k[k]] = x_kn[k,0:N_k[k]]

  # observable for estimation is the position^2
  elif observe == 'position^2':
    state_dependent = False
    A_kn = numpy.zeros([K,N_max], dtype = numpy.float64)
    for k in range(0,K):
      A_kn[k,0:N_k[k]] = x_kn[k,0:N_k[k]]**2

  (A_k_estimated, dA_k_estimated) = mbar.computeExpectations(A_kn, state_dependent = state_dependent)

  # need to additionally transform to get the square root
  if observe == 'RMS displacement':
    A_k_estimated = numpy.sqrt(A_k_estimated)
    # Compute error from analytical observable estimate.
    dA_k_estimated = dA_k_estimated/(2*A_k_estimated)

  As_k_estimated = numpy.zeros([K],numpy.float64)
  dAs_k_estimated = numpy.zeros([K],numpy.float64)

  # 'standard' expectation averages - not defined if no samples
  nonzeros = numpy.arange(K)[Nk_ne_zero]

  totaln = 0
  for k in nonzeros:

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computeExpectations [as 别名]

#.........这里部分代码省略.........
    #pertepsi *= 2
    #pertsig = numpy.array([sig_samp_space[1]])
    #pertsig *= 1.5
    #pertNk = numpy.array([energies.N_k[1]])
    #pertrdfs = rdfs[:,1,:]
 

    #Begin computing expectations
    Erdfs = numpy.zeros([Npert, basebins])
    dErdfs = numpy.zeros([Npert, basebins])
    eff_u_kln = numpy.zeros([energies.nstates, energies.nstates + Npert, energies.itermax])
    eff_u_kln[:energies.nstates, :energies.nstates, :] = energies.u_kln
    
    for l in xrange(Npert):
        #if not savedata or not os.path.isfile('lj%s/rdf%sfromn%s.npz'%(nstates-Npert+l, nstates-Npert+l, nstates)):
        if crunchset:
            filecheck = filestring % (nstates, pertname[l])
        else:
            filecheck = filestring % (fileobjs[0]+l, fileobjs[1]+l, fileobjs[2])
        if not savedata or not os.path.isfile(filecheck) or effnum:
            print "Working on state %d/%d" %(l+1, Npert)
            q = pertq[l]
            epsi = pertepsi[l]
            sig = pertsig[l]
            #u_kn_P = numpy.zeros([nstates,energies.itermax])
            u_kn_P = flamC12sqrt(epsi,sig)*energies.const_R_matrix + flamC6sqrt(epsi,sig)*energies.const_A_matrix + flamC1(q)*energies.const_q_matrix + flamC1(q)**2*energies.const_q2_matrix + energies.const_unaffected_matrix
            if effnum:
                eff_u_kln[:,energies.nstates+l,:] = u_kn_P
            else:
                for bin in xrange(basebins):
                    stdout.flush()
                    stdout.write('\rWorking on bin %d/%d'%(bin,basebins-1))
                    #Erdfs[l, bin], dErdfs[l, bin] = mbar.computePerturbedExpectation(u_kn_P, rdfs[bin,:,:])
                    Erdfs[l, bin], dErdfs[l, bin] = mbar.computeExpectations(rdfs[bin,:,:], u_kn=u_kn_P, compute_uncertainty=False)
                stdout.write('\n')
            if savedata:
                #savez('lj%s/rdf%sfromn%s.npz'%(nstates-Npert+l, nstates-Npert+l, nstates), Erdfs=Erdfs[l,:], dErdfs=dErdfs[l,:])
                if crunchset:
                    savez(filestring % (nstates, pertname[l]), Erdfs=Erdfs[l,:], dErdfs=dErdfs[l,:])
                else:
                    savez(filestring % (fileobjs[0]+l, fileobjs[1]+l, fileobjs[2]), Erdfs=Erdfs[l,:], dErdfs=dErdfs[l,:])
        else:
            #rdfdata = numpy.load('lj%s/rdf%sfromn%s.npz'%(nstates-Npert+l, nstates-Npert+l, nstates))
            if crunchset:
                rdfdata = numpy.load(filestring % (nstates, pertname[l]))
            else:
                rdfdata = numpy.load(filestring % (fileobjs[0]+l, fileobjs[1]+l, fileobjs[2]))
            Erdfs[l,:] = rdfdata['Erdfs']
            dErdfs[l,:] = rdfdata['dErdfs']

    #Effective number of samples
    if effnum:
        N_k = numpy.zeros(energies.nstates+Npert, dtype=numpy.int32)
        f_k = numpy.zeros(energies.nstates+Npert)
        N_k[:energies.nstates] = energies.N_k
        f_k[:mbar.K] = mbar.f_k
        eff_mbar = MBAR(eff_u_kln, N_k, verbose = True, initial_f_k=f_k, subsampling_protocol=[{'method':'L-BFGS-B','options':{'disp':True}}], subsampling=1)
        for l in xrange(Npert):
            w = numpy.exp(eff_mbar.Log_W_nk[:,mbar.K+l])
            neff = 1.0/numpy.sum(w**2)
            print "Effective number of samples for {0:s}: {1:f}".format(pertname[l], neff)
        pdb.set_trace()


    if bootstrap: #I spun this off into its own section, although it probably could have been folded into the previous
        from numpy.random import random_integers

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computeExpectations [as 别名]
  
  # observable for estimation is the position
  elif observe == 'position':
    state_dependent = False
    A_kn = numpy.zeros([K,N_max], dtype = numpy.float64)
    for k in range(0,K):
      A_kn[k,0:N_k[k]] = x_kn[k,0:N_k[k]]

  # observable for estimation is the position^2
  elif observe == 'position^2':
    state_dependent = False
    A_kn = numpy.zeros([K,N_max], dtype = numpy.float64)
    for k in range(0,K):
      A_kn[k,0:N_k[k]] = x_kn[k,0:N_k[k]]**2

  results = mbar.computeExpectations(A_kn, state_dependent = state_dependent)
  A_k_estimated = results['mu']
  dA_k_estimated = results['sigma']

  # need to additionally transform to get the square root
  if observe == 'RMS displacement':
    A_k_estimated = numpy.sqrt(A_k_estimated)
    # Compute error from analytical observable estimate.
    dA_k_estimated = dA_k_estimated/(2*A_k_estimated)

  As_k_estimated = numpy.zeros([K],numpy.float64)
  dAs_k_estimated = numpy.zeros([K],numpy.float64)

  # 'standard' expectation averages - not defined if no samples
  nonzeros = numpy.arange(K)[Nk_ne_zero]