Python MBAR.computePerturbedExpectation方法代码示例

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computePerturbedExpectation [as 别名]
      A_kn[k,0:N_k[k]] = (x_kn[k,0:N_k[k]] - O_extra[nth])**2 # observable is the squared displacement

  # observable is the potential energy, a 3D array since the potential energy is a function of 
  # thermodynamic state
  elif observe == 'potential energy':
    A_kn = unew_kln[:,nth,:]/beta

  # position and position^2 can use the same observables  
  # observable for estimation is the position
  elif observe == 'position': 
    A_kn = A_kn_all['position']

  elif observe == 'position^2': 
    A_kn = A_kn_all['position^2']

  (A_k_estimated, dA_k_estimated) = mbar.computePerturbedExpectation(unew_kln[:,nth,:],A_kn)

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

  A_k_error = A_k_estimated - A_k_analytical[observe][nth]

  print "Analytical estimator of %s is" % (observe)
  print A_k_analytical[observe][nth]

  print "MBAR estimator of the %s is" % (observe)
  print A_k_estimated

  print "MBAR estimators differ by X standard deviations"

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

#.........这里部分代码省略.........
        except:
            raise
        nkbin, nkiter = rdfk.shape
        rdfs[:,k,:nkiter] = rdfk
  
    #Do halfs
    if sliceset is 1:
        darange = xrange(0,25)
        nrange = 25
    elif sliceset is 2:
        darange = xrange(25,Nparm)
        nrange = 25
    #Do fifths
    elif sliceset is 3:
        darange = xrange(0,10)
        nrange = 10
    elif sliceset is 4:
        darange = xrange(10,20)
        nrange = 10
    elif sliceset is 5:
        darange = xrange(20,30)
        nrange = 10
    elif sliceset is 6:
        darange = xrange(30,40)
        nrange = 10
    elif sliceset is 7:
        darange = xrange(40,Nparm)
        nrange = Nparm-40
    #Do Everything
    else:
        darange = xrange(Nparm)
        nrange = Nparm
    #Generate the rmins for each lj combo
    #                    epsi   sig
    rmins = numpy.zeros([Nparm, Nparm])
    q = 0
    run_start_time = time.time()
    number_of_iterations = Nparm*nrange
    iteration = 0
    gr = (numpy.sqrt(5)-1)/2.0
    #Set up units for closest point
    Uplotx = units.Quantity(plotx, units.nanometer)
    zeroR = numpy.zeros(plotx.shape)
    OWsig = 3.15061 * units.angstrom
    OWepsi = 0.6364 * units.kilojoules_per_mole
    for iepsi in darange:
        epsi = epsi_range[iepsi]
        for isig in xrange(Nparm):
            if not os.path.isfile('esq_%s/RDFns%iE%iS%i.npz' %(spacename, nstates, iepsi, isig)):
                initial_time = time.time()
                iteration += 1
                print "Working on e %i and s %i" % (iepsi, isig)
                sig = sig_range[isig]
                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
                #Estimate zero-th order RDF
                guessepsi = geo(epsi*units.kilojoules_per_mole,OWepsi)
                guesssig = geo(sig*units.nanometer, OWsig)
                rdfguess = exp0(Uplotx, guessepsi, guesssig)
                depart0 = plotx[numpy.isclose(rdfguess,zeroR)].max()
                #Find nearest value in array
                #rcurrent = numpy.argmin(numpy.abs(sig*0.75-plotx))
                rcurrent = numpy.argmin(numpy.abs(depart0-plotx))
                foundmin = False
                #Estimate current min
                Ecurrent = mbar.computePerturbedExpectation(u_kn_P, rdfs[rcurrent,:,:], compute_uncertainty=False)
                if numpy.isclose(Ecurrent, 0):
                    #Case where sigma already is g=0
                    while foundmin is False:
                        rnew = rcurrent + 1
                        Enew = mbar.computePerturbedExpectation(u_kn_P, rdfs[rnew,:,:], compute_uncertainty=False)
                        if numpy.isclose(Enew, 0):
                            #Have not found where g(r)=0 and g(r+dr)>0
                            rcurrent = rnew
                        else:
                            foundmin = True
                else:
                    #Case where sigma is not g=0
                    while foundmin is False:
                        rnew = rcurrent - 1
                        Enew = mbar.computePerturbedExpectation(u_kn_P, rdfs[rnew,:,:], compute_uncertainty=False)
                        #Reversed logic of previous since we are searching backwards
                        if not numpy.isclose(Enew, 0):
                            #Have not found where g(r)=0 and g(r+dr)>0
                            rcurrent = rnew
                        else:
                            foundmin = True
                #rmins[iepsi,isig] = rcurrent
                laptime = time.clock()
                # Show timing statistics. copied from Repex.py, copywrite John Chodera
                final_time = time.time()
                elapsed_time = final_time - initial_time
                estimated_time_remaining = (final_time - run_start_time) / (iteration) * (number_of_iterations - iteration)
                estimated_total_time = (final_time - run_start_time) / (iteration) * (number_of_iterations)
                estimated_finish_time = final_time + estimated_time_remaining
                print "Iteration took %.3f s." % elapsed_time
                print "Estimated completion in %s, at %s (consuming total wall clock time %s)." % (str(datetime.timedelta(seconds=estimated_time_remaining)), time.ctime(estimated_finish_time), str(datetime.timedelta(seconds=estimated_total_time)))
                savez('esq_%s/RDFns%iE%iS%i.npz' %(spacename, nstates, iepsi, isig), rmins=numpy.array([plotx[rnew]]))
            else:
                rmins[iepsi,isig] = numpy.load('esq_%s/RDFns%iE%iS%i.npz' %(spacename, nstates, iepsi, isig))['rmins']
    savez('LJEffHS.npz', rmins=rmins)

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

#.........这里部分代码省略.........
            #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
        Nboot = bootstrapcount
        ErdfsB = numpy.zeros([Nboot, Npert, basebins])
        run_start_time = time.time()
        number_of_iterations = Npert*Nboot
        iteration = 0
        for l in xrange(Npert): #Go through each bootstrap
            print "Bootstraping 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
            for Nb in xrange(Nboot): #Bootstrap
                initial_time=time.time()
                iteration += 1
                u_kn_B = numpy.zeros(u_kn_P.shape)
                rdfsB = numpy.zeros(rdfs.shape)
                for k in xrange(nstates): #Roll random numbers
                    N_k = energies.N_k[k]
                    samplepool = random_integers(0,N_k-1,N_k)
                    for i in xrange(len(samplepool)):
                        u_kn_B[k,i] = u_kn_P[k,samplepool[i]] #Have to do this since slicing the first index alone reverses the return order for the last 2 indicies
                        rdfsB[:,k,i] = rdfs[:,k,samplepool[i]]
                #for bin in xrange(basebins):
                for bin in xrange(60): #Cut down for time
                    stdout.flush()
                    #stdout.write('\rWorking on bootstrap %d/%d and bin %d/%d'%(Nb,Nboot-1,bin,basebins-1))
                    stdout.write('\rWorking on bootstrap %d/%d and bin %d/%d'%(Nb,Nboot-1,bin,60-1))
                    ErdfsB[Nb, l, bin] = mbar.computePerturbedExpectation(u_kn_B, rdfsB[bin,:,:], compute_uncertainty=False)[0] #Dont compute error for speed
                laptime = time.clock()
                # Show timing statistics. copied from Repex.py, copywrite John Chodera
                final_time = time.time()
                elapsed_time = final_time - initial_time
                estimated_time_remaining = (final_time - run_start_time) / (iteration) * (number_of_iterations - iteration)
                estimated_total_time = (final_time - run_start_time) / (iteration) * (number_of_iterations)
                estimated_finish_time = final_time + estimated_time_remaining
                print "Iteration took %.3f s." % elapsed_time
                print "Estimated completion in %s, at %s (consuming total wall clock time %s)." % (str(datetime.timedelta(seconds=estimated_time_remaining)), time.ctime(estimated_finish_time), str(datetime.timedelta(seconds=estimated_total_time)))
        savez('BootstrapedIonRDF.npz', Erdfs = ErdfsB)
                    

    if singleparms is None and not crunchset:
        f,a = plt.subplots(Npert, 2)
        for l in xrange(Npert):
            Nl = pertNk[l]
            rdftraj = pertrdfs[:,l,:Nl]
            if Npert == 1:
                a[0].plot(plotx, Erdfs[l,:], '-k')
                #a[0].plot(plotx, Erdfs[l,:] + dErdfs[l,:], '--k')
                #a[0].plot(plotx, Erdfs[l,:] - dErdfs[l,:], '--k')
                a[1].plot(plotx, rdftraj.sum(axis=1)/float(Nl), '-k')
                alims0 = a[0].get_ylim()
                alims1 = a[1].get_ylim()
                a[0].set_ylim(0, numpy.amax(numpy.array([alims0[1], alims1[1]])))
                a[1].set_ylim(0, numpy.amax(numpy.array([alims0[1], alims1[1]])))
            else:
                a[l,0].plot(plotx, Erdfs[l,:], '-k')
                #a[l,0].plot(plotx, Erdfs[l,:] + dErdfs[l,:], '--k')
                #a[l,0].plot(plotx, Erdfs[l,:] - dErdfs[l,:], '--k')
                a[l,1].plot(plotx, rdftraj.sum(axis=1)/float(Nl), '-k')
                alims0 = a[l,0].get_ylim()
                alims1 = a[l,1].get_ylim()
                a[l,0].set_ylim(0, numpy.amax(numpy.array([alims0[1], alims1[1]])))
                a[l,1].set_ylim(0, numpy.amax(numpy.array([alims0[1], alims1[1]])))
    elif crunchset:
        f,a = plt.subplots(Npert)
        for l in xrange(Npert):
            a[l].plot(plotx, Erdfs[l,:], '-k')
            a[l].set_ylim([0,a[l].get_ylim()[1]])
            print "Peak for %s at %f" % (pertname[l], plotx[numpy.argmax(Erdfs[l,:])])
    else:
        f,a=plt.subplots(1,1)
        a.plot(plotx, Erdfs[0,:], '-k')
        
    plt.show()
    pdb.set_trace()