Tutorials: Chromosome Chromosome Optimization

Import the necessary packages

import numpy as np
import os
import subprocess
import sys
import copy

Important parameters

We defined the important parameters in the next blocks.

run_number = 1 #The index of the current iteration
N_replicas = 1 #The number of the total replicas

nAtom  = 60642 #The number of the atoms in our system
nIdeal = 2490-1 #The number of the ideal parameters
ncompt = 3 #The number of the compartments
nChrom = 22 #The number of the chromosomes including inter-chromosomal interactions, excluding the sex chromosome
nInter = int(nChrom*(nChrom-1)/2)
ncv    = int(nIdeal + ncompt*(ncompt+1)/2 + nInter)

chop_ideals             = 1000

## Adam optimizer parameters
m_dw                    = np.loadtxt('adam_chr_chr_param/%02d/mdw.txt'%(run_number-1))
v_dw                    = np.loadtxt('adam_chr_chr_param/%02d/vdw.txt'%(run_number-1))
m_db                    = np.loadtxt('adam_chr_chr_param/%02d/mdb.txt'%(run_number-1))
v_db                    = np.loadtxt('adam_chr_chr_param/%02d/vdb.txt'%(run_number-1))
beta1                   = 0.9
beta2                   = 0.999
epsilon                 = 1e-8
t                       = int(np.loadtxt('adam_chr_chr_param/%02d/t.txt'%(run_number-1)))

start_cv                = 1 #The starting index of the parameters that we hope to optimize
end_cv                  = 2726 #The ending index of the parameters that we hope to optimize

old_iter                = run_number-1

cvInd   = np.zeros((ncv, ), dtype=float)
irun    = 0

We normalize the contact probabilities firstly, and people can have their own path for the contact probs.

for replica in range(1,N_replicas+1,1): #We show the situation when the N_replicas = 1
    #If simulation didn't complete
    if os.path.exists("contact_prob/contact_prob.txt"):
        cvInd   += np.loadtxt("contact_prob/contact_prob.txt")
        irun    += np.loadtxt("contact_prob/nframes.txt")
cvInd /= irun

We load the experimental results and scale them to the same magnitude of the simulated contact probabilities.

expt        = np.loadtxt("expt_constraints_HFF_100KB.txt")
expt        = cvInd[0]/expt[0]*expt #We scale the experimental contacts to the same magnitude of simulated results

grad        = -cvInd + expt

Adam optimization

## START TO DO THE ADAM OPTIMIZATION
## momentum beta 1
# *** weights *** #
m_dw        = beta1*m_dw + (1-beta1)*grad
# *** biases *** #
m_db        = beta1*m_db + (1-beta1)*grad
## rms beta 2
# *** weights *** #
v_dw        = beta2*v_dw + (1-beta2)*(grad**2)
# *** biases *** #
v_db        = beta2*v_db + (1-beta2)*grad

We save the updated parameters in Adam for the next iteration

subprocess.call(["mkdir -p adam_chr_chr_param/%02d"%(run_number)],shell=True,stdout=subprocess.PIPE)
np.savetxt('adam_chr_chr_param/%02d/mdw.txt'%(run_number), m_dw.reshape((-1,1)), fmt='%15.12e')
np.savetxt('adam_chr_chr_param/%02d/vdw.txt'%(run_number), v_dw.reshape((-1,1)), fmt='%15.12e')
np.savetxt('adam_chr_chr_param/%02d/mdb.txt'%(run_number), m_db.reshape((-1,1)), fmt='%15.12e')
np.savetxt('adam_chr_chr_param/%02d/vdb.txt'%(run_number), v_db.reshape((-1,1)), fmt='%15.12e')
np.savetxt('adam_chr_chr_param/%02d/t.txt'%(run_number), np.array([t+1]).reshape((-1,1)), fmt='%d')

## bias correction
m_dw_corr   = m_dw/(1-beta1**t)
m_db_corr   = m_db/(1-beta1**t)
v_dw_corr   = v_dw/(1-beta2**t)
v_db_corr   = v_db/(1-beta2**t)

dalpha     = m_dw_corr/(np.sqrt(v_dw_corr)+epsilon)

Start to update the parameters

#We set all the default learning rates as 0.001
eta1           = 0.001
eta2           = 0.001
eta3           = 0.001

#Load the old potential
ideal_old         = np.loadtxt('potential/%02d/ideal_param_file.txt'%(run_number-1), usecols=[1])
compt_old         = np.loadtxt('potential/%02d/compt_param_file.txt'%(run_number-1), usecols=[2])
inter_old         = np.loadtxt('potential/%02d/interchr_param_file.txt'%(run_number-1), usecols=[2])

ideal_new = copy.deepcopy(ideal_old)
compt_new = copy.deepcopy(compt_old)
inter_new = copy.deepcopy(inter_old)

ideal_new[2:chop_ideals+1]          -= dalpha[1:chop_ideals]*eta1
compt_new[np.array([0,1,2,4,5,7])]  -= dalpha[nIdeal:nIdeal+ncompt*(ncompt+1)//2]*eta2 #The 6 positions represent the 6 type-type interactions
inter_new                           -= dalpha[nIdeal+ncompt*(ncompt+1)//2:]*eta3

Save all the updated parameters as the next iteration.

#Save the new potential
subprocess.call(["mkdir -p potential/%02d"%(run_number)],shell=True,stdout=subprocess.PIPE)

ideal_idx = np.loadtxt('potential/%02d/ideal_param_file.txt'%(run_number-1), usecols=[0], dtype='int')
compt_idx = np.loadtxt('potential/%02d/compt_param_file.txt'%(run_number-1), usecols=[0,1], dtype='int')
inter_idx = np.loadtxt('potential/%02d/interchr_param_file.txt'%(run_number-1), usecols=[0,1], dtype='int')

np.savetxt('potential/%02d/ideal_param_file.txt'%run_number,
           np.hstack((ideal_idx.reshape((-1,1)), ideal_new.reshape((-1,1)))), fmt='%d %.6f')
np.savetxt('potential/%02d/compt_param_file.txt'%run_number,
           np.hstack((compt_idx, compt_new.reshape((-1,1)))), fmt='%d %d %.6f')
np.savetxt('potential/%02d/interchr_param_file.txt'%run_number,
           np.hstack((inter_idx, inter_new.reshape((-1,1)))), fmt='%d %d %.6f')