<commit_before>import numpy as np
import helper
import glob
import sys
import argparse
from scipy.ndimage import filters
from pyemma.coordinates.clustering import AssignCenters
import itertools

"""
    Script that computes the absolute binding free energy of an MSM.
    It needs:
        1) Trajectory wildcard in order to compute cluster volume
        2) A discretized/clusterCenters.dat file with the cluster center
        3) A MSM_object.pkl obtained with pyemma in order to obtain the stationary distribution
        4) For the moment, it needs of a reweightingT, in order to do a histogram reweighting, but does not seem to work that well
"""

def assignNewTrajecories(trajs, clusterCenters):
    assign = AssignCenters(clusterCenters)
    dTrajs = assign.assign(trajs)
    return dTrajs

def expandTrajs(trajList):
    d = 0.5
    combinations = np.array([[0,1,0], [1,0,0], [0,0,1], [0,-1,0], [-1,0,0], [0,0,-1]])
    return list(trajList+combinations*d)

def parseArgs():
    parser = argparse.ArgumentParser(description="Script that computes delta G")
    parser.add_argument('trajectory', type=str)
    parser.add_argument('-r', '--reweight', type=float, default=1000)
    args = parser.parse_args()
    return args.trajectory, args.reweight

def writePDB(pmf_xyzg, title="clusters.pdb"):
    templateLine = "HETATM%s  H1  CLT L 502    %s%s%s  0.75%s           H\n"

    content = ""
    for i, line in enumerate(pmf_xyzg):
        number = str(i).rjust(5)
        x = ("%.3f"%line[0]).rjust(8)
        y = ("%.3f"%line[1]).rjust(8)
        z = ("%.3f"%line[2]).rjust(8)
        g = ("%.3f"%line[3]).rjust(8)

        content += templateLine%(number, x, y, z, g)
    f = open(title, 'w')
    f.write(content)
    f.close()

def calcQ(T, gpmf):
    Q = np.sum(np.exp(-gpmf/(kb*T)))
    return Q

def reweightProbabilities(T, Torig, origProb):
    """
        Torig: temperature at which simulations were run
    """
    kb = 0.0019872041
    gpmf = -kb*Torig*np.log(origProb)
    gpmf[gpmf == -np.inf] = np.inf
    gpmf -= gpmf.min() #It does not make any difference here

    print gpmf[gpmf == np.inf]
    print gpmf[gpmf == -np.inf]
    a = (T-Torig) / (kb * T * Torig)
    correction1 = np.exp(a * gpmf)
    beta = 1/kb/T
    betaOrig = 1/kb/Torig
    Qorig = np.sum(np.exp(-betaOrig*gpmf))
    Qnew = np.sum(np.exp(-beta*gpmf))
    correction2 = Qorig / Qnew

    return correction1 * correction2 * origProb


def main(trajWildcard, reweightingT=1000):
    #clusteringObject = helper.loadMSM('clustering_object.pkl')
    #allClusters = clusteringObject.clustercenters
    allClusters = np.loadtxt("discretized/clusterCenters.dat")

    MSMObject = helper.loadMSM('MSM_object.pkl')
    pi = MSMObject.stationary_distribution
    np.savetxt("stationaryDist_small.dat", pi)


    r = allClusters[MSMObject.connected_sets[0]]


    #filename = "output.txt"
    #data = np.loadtxt(filename)

    #r = data[:,0:3]
    #pi = data[:,3]

    d = 0.75

    originalFilenames = glob.glob(trajWildcard)

    originalCoordinates = []
    for i, originalFilename in enumerate(originalFilenames):
        trajOriginalCoordinates = list(np.loadtxt(originalFilename, usecols=(1,2,3)))<commit_msg>Add pseudocounts to count matrix to avoid large  errors due to disconnected clusters<commit_after>import numpy as np
import helper
import glob
import sys
import argparse
from scipy.ndimage import filters
from pyemma.coordinates.clustering import AssignCenters
from AdaptivePELE.pyemma_scripts import runMarkovChainModel as run
import itertools

"""
    Script that computes the absolute binding free energy of an MSM.
    It needs:
        1) Trajectory wildcard in order to compute cluster volume
        2) A discretized/clusterCenters.dat file with the cluster center
        3) A MSM_object.pkl obtained with pyemma in order to obtain the stationary distribution
        4) For the moment, it needs of a reweightingT, in order to do a histogram reweighting, but does not seem to work that well
"""

def assignNewTrajecories(trajs, clusterCenters):
    assign = AssignCenters(clusterCenters)
    dTrajs = assign.assign(trajs)
    return dTrajs

def expandTrajs(trajList):
    d = 0.5
    combinations = np.array([[0,1,0], [1,0,0], [0,0,1], [0,-1,0], [-1,0,0], [0,0,-1]])
    return list(trajList+combinations*d)

def parseArgs():
    parser = argparse.ArgumentParser(description="Script that computes delta G")
    parser.add_argument('trajectory', type=str)
    parser.add_argument('-r', '--reweight', type=float, default=1000)
    args = parser.parse_args()
    return args.trajectory, args.reweight

def writePDB(pmf_xyzg, title="clusters.pdb"):
    templateLine = "HETATM%s  H1  CLT L 502    %s%s%s  0.75%s           H\n"

    content = ""
    for i, line in enumerate(pmf_xyzg):
        number = str(i).rjust(5)
        x = ("%.3f"%line[0]).rjust(8)
        y = ("%.3f"%line[1]).rjust(8)
        z = ("%.3f"%line[2]).rjust(8)
        g = ("%.3f"%line[3]).rjust(8)

        content += templateLine%(number, x, y, z, g)
    f = open(title, 'w')
    f.write(content)
    f.close()

def calcQ(T, gpmf):
    Q = np.sum(np.exp(-gpmf/(kb*T)))
    return Q

def reweightProbabilities(T, Torig, origProb):
    """
        Torig: temperature at which simulations were run
    """
    kb = 0.0019872041
    gpmf = -kb*Torig*np.log(origProb)
    gpmf[gpmf == -np.inf] = np.inf
    gpmf -= gpmf.min() #It does not make any difference here

    print gpmf[gpmf == np.inf]
    print gpmf[gpmf == -np.inf]
    a = (T-Torig) / (kb * T * Torig)
    correction1 = np.exp(a * gpmf)
    beta = 1/kb/T
    betaOrig = 1/kb/Torig
    Qorig = np.sum(np.exp(-betaOrig*gpmf))
    Qnew = np.sum(np.exp(-beta*gpmf))
    correction2 = Qorig / Qnew

    return correction1 * correction2 * origProb


def main(trajWildcard, reweightingT=1000):
    #clusteringObject = helper.loadMSM('clustering_object.pkl')
    #allClusters = clusteringObject.clustercenters
    allClusters = np.loadtxt("discretized/clusterCenters.dat")

    MSMObject = helper.loadMSM('MSM_object.pkl')
    if len(allClusters) == MSMObject.stationary_distribution.size:
        pi = MSMObject.stationary_distribution
        r = allClusters[MSMObject.connected_sets[0]]
    else:
        ######################
        # Reconstruct stationary distribution with pseudocounts to ensure
        # connectivity
        counts = MSMObject.count_matrix_full
        counts += 1/float(counts.shape[0])
        trans = run.buildRevTransitionMatrix(counts)
        eiv, eic = run.getSortedEigen(trans)
        pi = run.getStationaryDistr(eic[:, 0])
        r = allClusters


    #filename = "output.txt"
    #data = np.loadtxt(filename)

    #r = data[:,0:3]
    #pi = data[:,3]

    d = 0.75

    originalFilenames = glob.glob(trajWildcard)

    originalCoordinates = []
    for i, originalFilename in enumerate(originalFilenames):
        trajOriginalCoordinates = list(np.loadtxt(originalFilename, usecols=(1,2,3)))
{'fraction_non_alphanumeric': 0.1294777729823047, 'fraction_numerical': 0.023162134944612285, 'mean_word_length': 3.231284236153378, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '664262', 'n_tokens_mistral': 2539, 'n_tokens_neox': 2438, 'n_words': 581}