Calculate SOAP potential for ranking MHC2 peptides.¶
Directory: SOAP/examples/mhc2/
Preprocess PDB if necessary and generate sequence files, example file: SOAP/examples/mhc2/pdb_X_2.2A_0.25rfree.pir. Please see the module
sequenceson how to generate different subsets of PDBs.Prepare Decoys. The original decoy files generated using modeller, mhc2_original.tar.gz and the prepared decoys, mhc2.tar.gz can be found in the example directory.
import SOAP.decoys
import shutil
import os
import sys
import glob
import pickle
# define the paths for the decoys; assume they are found under the directory
# containing this script
workingPath = os.path.abspath(os.path.dirname(sys.argv[0]))
originalDecoyPath = os.path.join(workingPath,'mhc2_original')
preparedDecoyPath = os.path.join(workingPath,'mhc2')
decoySetName='mhc2'
#copy the files from originalDecoyPath to target path
"""The pre-prepared decoys directory should look like the following:
root:
decoySet1: could use the PDB code of the receptor here
Decoy files
[*base*]: if this file is present in this directory, it should be the receptor file shared by all decoys, and the final decoys are generated by combining this
rmsd.pickle
decoySet2
...
"""
fl=os.listdir(originalDecoyPath)
for f in fl:
if not os.path.isdir(os.path.join(originalDecoyPath,f)):
continue
shutil.rmtree(os.path.join(preparedDecoyPath,f), ignore_errors=True)
os.makedirs(os.path.join(preparedDecoyPath,f))
os.chdir(os.path.join(preparedDecoyPath,f))
for d in glob.glob(os.path.join(originalDecoyPath,
f, 'selected_peptides','*')):
shutil.copy(d, '.')
shutil.move('pMHC_'+f+'.pdb', 'base.pdb')
# Touch
open('needattachtobase', 'w')
rl=open(os.path.join(originalDecoyPath,f,'all_scores.txt')).read().split('\n')
rmsddict={}
scoredict={}
for item in rl:
if len(item)<10:
continue
il=item.split(',')
rmsddict[il[0].strip()[:-4]]=float(il[2])
scoredict[il[0].strip()[:-4]]=float(il[1])
pickle.dump(rmsddict,open('rmsd.pickle','wb'))
pickle.dump(scoredict,open('score.pickle','wb'))
# the decoys are now in the format that SOAP can process;
# use SOAP.decoys to build the DecoySet
do=SOAP.decoys.DecoySet(dsname=decoySetName, sourcedir=preparedDecoyPath)
do.build()
Run SOAP script to select models.
Select the best models with distance features alone, only considering interface atom pairs.
from __future__ import print_function
from SOAP import *
#define recovery fucntions
rfeature=[['d',0,20,0.05]]
par={'uniform':4,'featurebins':rfeature[0],'distribute':'lownumberpriority','firstpoint':2.25}
slo={'key':'sacnfflex',
'valueset':{
'sacnfflex':{'sftype':'sacnf','par':par,'parvalue':[1,1,1,1]},
'sacnf52':{'sftype':'sacnf','par':[2.75,3.75],'parvalue':[1,1,1,1]},
'sacnf53':{'sftype':'sacnf','par':[1.75,2.75,3.75],'parvalue':[1,1,1,1,1]},
'sacnf62':{'sftype':'sacnf','par':[2.25,3.75],'parvalue':[1,1,1,1]},
'sacnf63':{'sftype':'sacnf','par':[2.25,3.75,4.75],'parvalue':[1,1,1,1,1]},
'sacnf64':{'sftype':'sacnf','par':[1.75,2.75,3.75,4.75],'parvalue':[1,1,1,1,1,1]},
'sacnf73':{'sftype':'sacnf','par':[2.25,3.75,5.25],'parvalue':[1,1,1,1,1]},
'sacnf74':{'sftype':'sacnf','par':[2.75,3.75,4.75,5.75],'parvalue':[1,1,1,1,1,1]},
'sacnf75':{'sftype':'sacnf','par':[1.75,2.75,3.75,4.75,5.75],'parvalue':[1,1,1,1,1,1,1]},
'sacnf84':{'sftype':'sacnf','par':[1.75,3.25,4.75,6.25],'parvalue':[1,1,1,1,1,1]},
'sacnf104':{'sftype':'sacnf','par':[1.75,3.75,5.75,7.75],'parvalue':[1,1,1,1,1,1]},
'sacnf105':{'sftype':'sacnf','par':[1.75,3.75,5.25,6.75,8.25],'parvalue':[1,1,1,1,1,1,1]},
'sacnf124':{'sftype':'sacnf','par':[2.75,4.75,6.75,8.75],'parvalue':[1,1,1,1,1,1]},
'sacnf125':{'sftype':'sacnf','par':[1.75,3.75,5.25,7.75,9.75],'parvalue':[1,1,1,1,1,1,1]},
'sacnf154':{'sftype':'sacnf','par':[2.75,5.75,8.75,11.75],'parvalue':[1,1,1,1,1,1]},
'sacnf156':{'sftype':'sacnf','par':[2.75,4.75,6.75,8.75,10.75,12.75],'parvalue':[1,1,1,1,1,1,1,1]},
'sacnf158':{'sftype':'sacnf','par':[2.75,4.25,5.75,7.25,8.75,10.25,11.75,13.25],'parvalue':[1,1,1,1,1,1,1,1,1,1]}
}}
#use more recovery function models, not recommended unless you are sure what this means and think this will help
if 0:
for valuetype in ['aon','aoo','acn','aco','1oo','1cn','0oo']:
if valuetype.startswith('a'):
gl=['e','o','l','f','p','s']
elif valuetype.startswith('1'):
gl=['e','o','f','p']
elif valuetype.startswith('0'):
gl=['l','s']
for gentype in gl:
splinetype='s'+partype+valuetype+converttype+gentype
slo[splinetype+'2']={'sftype':splinetype,'par':[2.25,3.75]}
slo[splinetype+'3']={'sftype':splinetype,'par':[2.25,3.75,4.75]}
slo[splinetype+'4']={'sftype':splinetype,'par':[1.75,2.75,3.75,4.75]}
#
ni=40
initvalues=list(np.arange(0,ni)/10.0+1)
initvalues=np.array(initvalues)
#scoring table definitions using raw stats from PDB and recovery functions
sfeature=[rfeature[0],'a158','as158']
pm=['','npend']
ref1={'type':'sf','features':rfeature,'sftype':slo,'par':slo,'parvalue':slo,'ratio':[1.0]}
scaledsp1={'type':'scaledsp','pdbset':'X2_2.2A_0.25rfree','features':sfeature,'genmethod':'cs1','pm':pm,'refs':[ref1],'ratio':[1.0]}
#define the benchmark scheme
bmtype={'type':'dsscore','dslist':['mhc2'],'criteria':'dcg','combine':'scoresum','bm':'cs1'}
#define the continuous parameter (model parameter) to be optimzied
search1={'object':ref1,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'dfire','values':initvalues}}
search2={'object':ref1,'key':'par','pos':[0,1,2,3],'InitialGenerator':{'type':'dfire','values':initvalues}}
#define the models to be searched
dsearch2={'object':rfeature[0],'key':2,'valueset':[20,5.6,5.8,6,7]}#[6,5.4,5,4.8,4.6,4.4,4.2]#15, 7,6,5,4.5,4,3.5,3#,10,8,7,6,5#15,12,10, 8,7,15,12,10, 8,7,6,
dsearch4={'object':par,'key':'uniform','valueset':[7,6,5,4,3]}
dsearch5={'object':rfeature[0],'key':3,'valueset':[0.05]}#0.1,0.025,0.05 #,0.25,0.1,0.05 #0.05,0.1,0.25,
dsearch7={'object':pm,'key':0,'valueset':['','pp5','pp1','ks0.2','ks0.5']}#'gps0.5','gps0.2'
dsearch8={'object':pm,'key':1,'valueset':['npend','nbsum']}
dsearch9={'object':[scaledsp1,scaledsp1],'key':['genmethod','pdbset'],'valueset':[['cs1','X2_2.2A_0.25rfree']]}#,['cs1','X2_2.2A_0.25rfree_30'],['cs1','X2_2.2A_0.25rfree_60'],['cs1','X2_2.2A_0.25rfree_95'],['bs20dsp','X_2.2A_0.25rfree'],['bs20dsp','X_2.2A_0.25rfree_30'],['bs20dsp','X_2.2A_0.25rfree_60'],['bs15dsp','X_2.2A_0.25rfree'],['bs10dsp','X_2.2A_0.25rfree']]}#,,'cs1'
#optimization and sampling parameters
inner={'improved':2,'maxloop':100,'minloop':2}
outer={'improved':4,'maxloop':5023}
td=list(np.sqrt(np.linspace(1,float(10)**2,ni)))
td=np.array(td)
tune_interval=200
sampleschedule={'inner':inner,'outer':outer}
ssm={'sm':'mcp','reset_betweenruns':2,'blockupdate':True, 'exchange_method':1,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':200,'add_bias':False,'temperature_distribution':td}
ssm0={'sm':'mcs','reset_betweenruns':2,'blockupdate':False,'using_globalbest':True,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':201,'add_bias':False,'temperature_distribution':td}
ssm2={'sm':'mcp','reset_betweenruns':2,'blockupdate':False, 'exchange_method':1,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':200,'add_bias':False,'temperature_distribution':td}
ssm20={'sm':'mcs','reset_betweenruns':2,'blockupdate':True,'using_globalbest':True,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':201,'add_bias':False,'temperature_distribution':td}
ssm1={'sm':'mca','reset_betweenruns':2,'blockupdate':True, 'using_globalbest':True,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':201,'add_bias':False,'temperature_distribution':np.zeros(ni)+2}
ssm3={'sm':'powell','blockupdate':False}
sml=[ssm20,ssm2,ssm0,ssm,ssm0,ssm,ssm0,ssm, ssm1]
#define the final model
model1={'scorers':[scaledsp1,ref1],'bmtype':bmtype,'searches':[search1,search2], 'runsperscorer':ni,
'dsearches':[dsearch2,dsearch5,dsearch4,dsearch7,dsearch8,dsearch9],'sml':sml,'cvk':2,'repeat':1,'fold':3}#,dsearch2,dsearch5,dsearch6,dsearch7 #,'testperc':0.33
#test the model on local host
if 0:
so=scorer(model=convert2old(model1))
print(so.assess_ideal())
print(so.assess_worst())
print(so.assess_sascore(slevel='top1_rmsdbbif_mean_'))
pdb.set_trace()
#model selection
spl=spss(model=model1)
#setting up the parameters for model selection
spl.currentcutoffpercinitial=0.0
spl.currentcutoffpercratio=0.0
spl.maximumround=10
#search for the best model
spl.find_best_par()
spl.log()
Select the best models with distance and accessibility features, multiple recovery functions, only considering interface atom pairs.
from __future__ import print_function
from SOAP import *
import os
import profile
import sys
import pdb
import numpy as np
#recovery function
ml=[]
rfeature=[['d',0,6,0.2]]
par={'object':rfeature, 'index1':0,'index2':2,'typelist':[['fixbin',0,1],['fixbin',1,2],['fixbin',1,3],['numafter',[1,3],2]],'parindex':0}
sfeature=[rfeature[0],'a158','as158']
pm=['','npend']
par={'uniform':5,'featurebins':rfeature[0],'distribute':'lownumberpriority','firstpoint':2.25}
slo={'key':'sacnfflex',
'valueset':{
'sacnfflex':{'sftype':'sacnf','par':par,'parvalue':[1,1,1,1]},
'psn3':{'sftype':'psn','par':range(1,15,5),'parvalue':range(3)},
'psn4':{'sftype':'psn','par':range(1,15,4),'parvalue':range(4)},
'psn5':{'sftype':'psn','par':range(1,15,3),'parvalue':range(5)},
'psn8':{'sftype':'psn','par':range(0,15,2),'parvalue':range(8)},
}}
#different step method
ni=40
ssl={}
initvalues=list(np.arange(0,ni)/10.0+1)
initvalues=np.array(initvalues)
bmtype={'type':'dsscore','dslist':['mhc2'],'criteria':'dcg','combine':'scoresum','bm':'cs1'}
blocks=[]
kk=0
for apn in [36]:#57
print(apn)
noc=np.load(runenv.libdir+'ap'+str(apn)+'.npy').max()+1
print(noc)
scorers=[]
searches=[]
for i in range(noc):
print(i)
sbm=['cs1','ap'+str(apn)+'n',str(i)]
ref1={'type':'sf','features':sfeature,'sftype':slo,'par':slo,'parvalue':slo,'ratio':[1.0+0.00001*i],'bm':sbm}
scaledsp1={'type':'scaledsp','pdbset':'X2_2.2A_0.25rfree_60','features':sfeature,'genmethod':'cs1','pm':pm,'refs':[ref1],'ratio':ref1['ratio'],'bm':sbm}
search1={'object':ref1,'key':'parvalue','pos':[0,1,2,3,4,5,6],'InitialGenerator':{'type':'dfire','values':initvalues}}
#search2={'object':ref1,'key':'par','pos':[0,1,2,3,4,5,6],'InitialGenerator':{'type':'dfire','values':initvalues}}
#search3={'object':ref1,'key':'ratio','pos':[0],'InitialGenerator':[[1.0+0.00001*i] for j in range(40)]}
scorers.append(ref1)
scorers.append(scaledsp1)
searches.append(search1)
#searches.append(search2)
#searches.append(search3)
blocks.append([kk,kk+1,kk+2])
kk=kk+3
ssl['full']=[scorers,searches]
#define the final mode
#sampling and searching parameters
inner={'improved':2,'maxloop':10,'minloop':1}
outer={'improved':4,'maxloop':200}
#inner={'improved':1,'maxloop':2,'minloop':1}
#outer={'improved':2,'maxloop':2}
td=list(np.sqrt(np.linspace(1,float(10)**2,ni)))
td=np.array(td)
sampleschedule={'inner':inner,'outer':outer}
ssm={'sm':'mcp','reset_betweenruns':2,'blockupdate':True, 'exchange_method':1,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':200,'add_bias':False,'temperature_distribution':td}
ssm0={'sm':'mcs','reset_betweenruns':2,'blockupdate':False,'using_globalbest':True,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':201,'add_bias':False,'temperature_distribution':td}
ssm2={'sm':'mcp','reset_betweenruns':2,'blockupdate':False, 'exchange_method':1,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':200,'add_bias':False,'temperature_distribution':td}
ssm20={'sm':'mcs','reset_betweenruns':2,'blockupdate':True,'using_globalbest':True,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':201,'add_bias':False,'temperature_distribution':td}
ssm1={'sm':'mca','reset_betweenruns':2,'blockupdate':True, 'using_globalbest':True,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':201,'add_bias':False,'temperature_distribution':np.zeros(ni)+2}
sml=[ssm20,ssm2,ssm0,ssm,ssm0,ssm,ssm0,ssm, ssm1]
#sml=[ssm1]
#sml=[ssm]
dsearch6={'object':ssm,'key':'temperature_distribution','valueset':[np.sqrt(np.linspace(1,float(10)**2,ni)),
np.sqrt(np.linspace(1,float(20)**2,ni)),
np.sqrt(np.linspace(1,float(30)**2,ni)),
np.sqrt(np.linspace(2,float(20)**2,ni)),
np.sqrt(np.linspace(3,float(30)**2,ni))]}
dsearch5={'object':ssm,'key':'add_bias','valueset':[False]}
dsearch4={'object':ssm,'key':'reset_betweenruns','valueset':[2.0]}#,1.56,1.57,1.58,1.59
dsearch3={'object':ssm,'key':'stepmethod','valueset':['mxmp2.0']}
dsearch2={'object':ssm,'key':'sm','valueset':['mcp']}
ssmodel={'index':'full','valueset':ssl}
dsearch1={'object':ssmodel,'key':'index','valueset':['full']}
dsearch6={'object':slo,'key':'key','valueset':slo['valueset'].keys()}
dsearch8={'object':inner,'key':'improved','valueset':[1,2,3]}
dsearch9={'object':inner,'key':'maxloop','valueset':[10,30,100]}
dsearch10={'object':inner,'key':'minloop','valueset':[2,5,7,10]}
dsearch11={'object':outer,'key':'improved','valueset':[4,6,8]}
dsearch12={'object':ssm0,'key':'sm','valueset':['mcs','mca']}
dsearch13={'object':ssm0,'key':'using_globalbest','valueset':[True,False]}
dsearch14={'object':ssm0,'key':'temperature_distribution','valueset':[np.zeros(2*ni)+2,td]}
dsearch15={'object':ssm,'key':'exchange_method','valueset':[1]}
dsearch16={'object':[ssm2,ssm20],'key':['blockupdate','blockupdate'],
'valueset':[[True,True],[False,False]]}
dsearch17={'object':[ssm2,ssm20,ssm,ssm0,ssm1],'key':['blockupdate','blockupdate','blockupdate','blockupdate','blockupdate'],
'valueset':[[True,True,True,True,True],[False,False,False,False,False],[True,True,False,False,False]]}
dsearch21={'object':par,'key':'uniform','valueset':[8,5]}
dsearch32={'object':rfeature[0],'key':2,'valueset':[15, 12,10,8,6]}#,10,8,7,6,5#15,12,10, 8,7,15,12,10, 8,7,6,
dsearch34={'object':par,'key':'uniform','valueset':[7,5,4]} #7,6,5,4,3
dsearch35={'object':rfeature[0],'key':3,'valueset':[0.05,0.1,0.5]}#0.025,0.05 #,0.25,0.1,0.05 #0.05,0.1,0.25,
dsearch39={'object':[scaledsp1,scaledsp1],'key':['genmethod','pdbset'],'valueset':[['cs1','X2_2.2A_0.25rfree'],['bs20dsp','X_2.2A_0.25rfree'],['bs10dsp','X_2.2A_0.25rfree'],['bs2dsp','X_2.2A_0.25rfree'],['bs5dsp','X_2.2A_0.25rfree']]}#,,'cs1'['bs20dsp','X_2.2A_0.25rfree'],
dsearch7={'object':pm,'key':0,'valueset':['','pp5','pp1','ks0.2']}#'gps0.6','gps0.3',,'ks0.2','ks0.4'
dsearch8={'object':pm,'key':1,'valueset':['nbsum']}
scorerlist=ssl['full'][0]
#scorerlist=bm['scorers']
rfeature=[['b',0,1.0,0.1]]
sfeature=[rfeature[0],'a158']
pm=['pp5','','npend']
sratio=[3.8]
scaledsp2={'type':'scaledsp','pdbset':'X_2.2A_0.25rfree','features':sfeature,'genmethod':'','pm':'nbsum','refs':[],'ratio':sratio}
ref2={'type':'sf','features':rfeature,'sftype':'bins','par':range(10),'parvalue':[1 for i in range(10)],'bm':'','ratio':sratio}
search2={'object':ref2,'key':'parvalue','pos':range(0,10),'InitialGenerator':{'type':'dfire','values':list(np.zeros(40))}}
search3={'object':scaledsp2,'key':'ratio','pos':[0],'InitialGenerator':[[3.8] for i in range(40)]}
scorerlist=scorerlist+[scaledsp2,ref2]
searchlist=ssl['full'][1]+[search2,search3]
blocks.append([kk,kk+1])
model1={'scorers':scorerlist,'bmtype':bmtype,'searches':searchlist,#'thread':4,
'dsearches':[dsearch32,dsearch34,dsearch35,dsearch39],'sml':sml,'cvk':2,'fold':3,'repeat':1,'initialpars':'all','runsperscorer':40}
if 0:
so=scorer(model=convert2old(model1))
print(so.assess_model())
opt=optimizer(scorer=so)
opt.get_initial_value()
print(so.assess(opt.initialvalue))
opt.optimize()
pdb.set_trace()
spso=sps(modellist=[convert2old(model1)])
#pdb.set_trace()
spso.cv()#spl.eval_allpars()
logpath=spso.cvlist[0].logpath
elif 0:
so=scorer(model=convert2old(model1))
print(so.assess_model())
opt=optimizer(scorer=so)
opt.get_initial_value()
optres=opt.optimize()
print(so.assess(opt.initialvalue))
spso=sps(modellist=[convert2old(model1)])
pdb.set_trace()
spso.cv()#spl.eval_allpars()
logpath=spso.cvlist[0].logpath
elif 0:
spl=spss(model=model1)
spl.eval_allpars()
pdb.set_trace()
spl=spss(model=model1) # set of statistical potentials/discrete searches
spl.currentcutoffpercinitial=0.0# no randomness on discrete searches
spl.currentcutoffpercratio=0.0
spl.maximumround=10
#spl.eval_allpars()
spl.find_best_par()#find the best parameters on discrete and continuous variables.
spl.log()
pdb.set_trace()
Select the best models with orientation feature, considering both the interface and the introchain atom pairs.
from SOAP import *
import os
import profile
import sys
import pdb
import numpy as np
rfeature=[['dt',0,6,0.2]]
rrfeature=[['dt',0,6,6]]
runenv.otherrun=False
par={'uniform':4,'featurebins':rfeature[0],'distribute':'lownumberpriority','firstpoint':2.25}
slo={'key':'sacnfflex',
'valueset':{
'sacnfflex':{'sftype':'sacnf','par':par,'parvalue':[1,1,1,1]},
'psn3':{'sftype':'psn','par':range(1,15,5),'parvalue':range(3)},
'psn4':{'sftype':'psn','par':range(1,15,4),'parvalue':range(4)},
'psn5':{'sftype':'psn','par':range(1,15,3),'parvalue':range(5)},
'psn8':{'sftype':'psn','par':range(0,15,2),'parvalue':range(8)},
}}
sfeature=[rfeature[0],'g6#180','gs6#180','h12#180#-180','t306','ts306']
pm=['','nbsum']
ni=40
initvalues=list(np.arange(0,ni)/10.0+1)
initvalues=np.array(initvalues)
ratio1=[1.0]
#protein-peptide interactions
ref1={'type':'sf','features':rfeature,'sftype':slo,'par':slo,'parvalue':slo,'ratio':ratio1}
ref2={'type':'sf','features':[rrfeature[0],'g6#180'],'sftype':'spline','par':[0,60,120,180],'parvalue':[1,1,1,1],'ratio':ratio1}
ref3={'type':'sf','features':[rrfeature[0],'gs6#180'],'sftype':'spline','par':[0,60,120,180],'parvalue':[1,1,1,1],'ratio':ratio1}
ref4={'type':'sf','features':[rrfeature[0],'h12#180#-180'],'sftype':'spline','par':[-180,-60,60,180],'parvalue':[1,1,1,1],'ratio':ratio1}
scaledsp1={'type':'scaledsp','pdbset':'X_2.2A_0.25rfree','features':sfeature,'genmethod':'bs20dsp','pm':pm,'refs':[ref1,ref2,ref3,ref4],'ratio':ratio1}
ratio2=[1.04]
#peptide intra-atom interactions
ref21={'type':'sf','features':'dt30#6','sftype':slo,'par':slo,'parvalue':slo,'ratio':ratio2,'bm':'bs2dsp'}
ref22={'type':'sf','features':'dt1#6g6#180','sftype':'spline','par':[0,60,120,180],'parvalue':[1,1,1,1],'ratio':ratio2,'bm':'bs2dsp'}
ref23={'type':'sf','features':'dt1#6gs6#180','sftype':'spline','par':[0,60,120,180],'parvalue':[1,1,1,1],'ratio':ratio2,'bm':'bs2dsp'}
ref24={'type':'sf','features':'dt1#6h12#180#-180','sftype':'spline','par':[-180,-60,60,180],'parvalue':[1,1,1,1],'ratio':ratio2,'bm':'bs2dsp'}
scaledsp2={'type':'scaledsp','pdbset':'X_2.2A_0.25rfree','features':'dt30#6g6#180gs6#180h12#180#-180t306ts306','genmethod':'bs2dsp','bm':'bs2dsp','pm':pm,'refs':[ref21,ref22,ref23,ref24],'ratio':ratio2}
ref31={'type':'sf','features':'r20','sftype':'bins','par':[],'parvalue':[0]*21,'ratio':[1.0],'bm':''}
#define the benchmark scheme
bmtype={'type':'dsscore','dslist':['mhc2'],'criteria':'dcg','combine':'scoresum','bm':'cs1'}
#define the search positions
search1={'object':ref1,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'dfire','values':initvalues}}
search2={'object':ref1,'key':'par','pos':[0,1,2],'InitialGenerator':{'type':'dfire','values':initvalues}}
#define the search positions
search3={'object':ref2,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'sin','values':np.arange(0,ni)/10.0}}
search4={'object':ref3,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'sin','values':np.arange(0,ni)/10.0}}
search5={'object':ref4,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'sin','values':np.arange(0,ni)/10.0}}
search21={'object':ref21,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'dfire','values':initvalues}}
search22={'object':ref21,'key':'par','pos':[0,1,2,3],'InitialGenerator':{'type':'dfire','values':initvalues}}
#define the search positions
search23={'object':ref22,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'sin','values':np.arange(0,ni)/10.0}}
search24={'object':ref23,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'sin','values':np.arange(0,ni)/10.0}}
search25={'object':ref24,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'sin','values':np.arange(0,ni)/10.0}}
search9={'object':scaledsp2,'key':'ratio','pos':[0],'InitialGenerator':[[0.99] for i in range(40)]}
search31={'object':ref31,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'random'}}
search32={'object':ref31,'key':'ratio','pos':[0],'InitialGenerator':[[0.99] for i in range(ni)]}
#discrete search options
dsearch1={'object':[scaledsp1,bmtype],'key':['genmethod','bm'],'valueset':[['bs1','bs1'],['bs2','bs2'],['bs3','bs3']
,['bs4','bs4'],['bs5','bs5'],['bs6','bs6'],['bs7','bs7'],['bs8','bs8'],['bs9','bs9'],['bs10','bs10']
,['bs10','bs1'],['bs10','bs2'],['bs10','bs3'],['bs10','bs4']]}
dsearch2={'object':ref1,'key':'features','valueset':['dt800#20','dt200#20','dt40#20','dt600#15','dt300#15','dt150#15','dt30#15','dt100#10']}
dsearch3={'object':ref2,'key':'features','valueset':['g24#180','g12#180','g6#180']}
dsearch4={'object':ref3,'key':'features','valueset':['gs24#180','gs12#180','gs6#180']}
dsearch5={'object':ref4,'key':'features','valueset':['h48#360','h24#360','gs12#360']}
dsearch6={'object':ref1,'key':'par','valueset':[[1,2,8,15],[1,3,8,15],[2,3,8,15],range(1,16,2),range(1,16),range(0,16)]}
dsearch7={'object':ref2,'key':'par','valueset':[[0,60,120,180],range(0,181,30),range(0,181,20),range(0,181,10)]}
dsearch8={'object':ref3,'key':'par','valueset':[[0,60,120,180],range(0,181,30),range(0,181,20),range(0,181,10)]}
dsearch10={'object':scaledsp1,'key':'features','valueset':['dt30#6g6#180gs6#180h12#180#-180t306ts306','dt12#6g6#180gs6#180h12#180#-180t306ts306']}
dsearch11={'object':scaledsp1,'key':'pdbset','valueset':['X_2.2A','X_2.2A_nopsi']}
dsearch1={'object':[scaledsp1,bmtype],'key':['genmethod','bm'],'valueset':[['bs1dsp','bs1dsp'],['bs2dsp','bs2dsp'],['bs4dsp','bs4dsp'],['bs10dsp','bs10dsp']]}#'bs4dsp','bs5dsp','bs9dsp',,'bs2dsp'
dsearch2={'object':rfeature[0],'key':2,'valueset':[15, 12,10,8,6,5]}#,10,8,7,6,5#15,12,10, 8,7,15,12,10, 8,7,6,
dsearch6={'object':slo,'key':'key','valueset':slo['valueset'].keys()}#slo['valueset'].keys()'psn5','psn3','psn4','psn8'
dsearch7={'object':pm,'key':0,'valueset':['','pp5','pp1','ks0.2']}#'gps0.6','gps0.3',,'ks0.2','ks0.4'
dsearch8={'object':pm,'key':1,'valueset':['nbsum']}
dsearch9={'object':[scaledsp1,scaledsp1],'key':['genmethod','pdbset'],'valueset':[['cs1','X2_2.2A_0.25rfree'],['bs20dsp','X_2.2A_0.25rfree']]}#,,'cs1'['bs20dsp','X_2.2A_0.25rfree'],
dsearch9={'object':scaledsp1,'key':'genmethod','valueset':['bs20dsp','bs15dsp','bs10dsp','bs4dsp','bs2dsp']}#,,'cs1'['bs20dsp','X_2.2A_0.25rfree'],
dsearch2={'object':[rfeature[0],rrfeature[0],rrfeature[0]],'key':[2,2,3],'valueset':[[6,6,6],[8,8,8]]}#,,[10,10,10],[12,12,12][5.4,5.4,5.4],[4.4,4.4,4.4],10,8,7,6,5#15,12,10, 8,7,15,12,10, 8,7,6,
dsearch4={'object':par,'key':'uniform','valueset':[12,8,6,4]} #7,6,5,4,3,[12,10,8,7,6,4]
dsearch5={'object':rfeature[0],'key':3,'valueset':[0.2]}#0.025,0.05 #,0.25,0.1,0.05 #0.05,0.1,0.25,
#sampling and seaching parameters
inner={'improved':2,'maxloop':11,'minloop':1}
outer={'improved':4,'maxloop':2001}
td=list(np.sqrt(np.linspace(1,float(10)**2,ni)))
td=np.array(td)
tune_interval=203
sampleschedule={'inner':inner,'outer':outer}
ssm={'sm':'mcp','reset_betweenruns':2,'blockupdate':False, 'exchange_method':1,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':tune_interval,'add_bias':False,'temperature_distribution':td}
ssm0={'sm':'mcs','reset_betweenruns':1,'blockupdate':False,'using_globalbest':True,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp1.0','tune_interval':tune_interval,'add_bias':False,'temperature_distribution':td}
ssm2={'sm':'mcp','reset_betweenruns':2,'blockupdate':False, 'exchange_method':1,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':tune_interval,'add_bias':False,'temperature_distribution':td}
ssm20={'sm':'mcs','reset_betweenruns':2,'blockupdate':False,'using_globalbest':True,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':tune_interval,'add_bias':False,'temperature_distribution':td}
ssm1={'sm':'mca','reset_betweenruns':2,'blockupdate':False, 'using_globalbest':True,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':tune_interval,'add_bias':False,'temperature_distribution':np.zeros(ni)+2}
ssm3={'sm':'powell','blockupdate':False}
sml=[ssm20,ssm2,ssm0,ssm,ssm0,ssm,ssm0,ssm, ssm1]
#sml=[ssm3,ssm1]
#sml=[ssm0]
dsearch99={'object':[ssm,ssm2],'key':['blockupdate','blockupdate'],'valueset':[[False]*2]}#,
model1={'scorers':[scaledsp1,ref1,ref2,ref3,ref4,scaledsp2,ref21,ref22,ref23,ref24,ref31],'bmtype':bmtype,'searches':[search1,search2,search3,search4,search5,search21,search22,search23,search24,search25,search9,search31,search32], 'runsperscorer':ni,#number of replicas per scorers
#'initialmodelpath':inipath,
'dsearches':[dsearch4,dsearch2,dsearch9],#[dsearch4,dsearch2,dsearch9],,dsearch2,dsearch9
'sml':sml,#search schedule
'cvk':2,#number of repeation of cross validation
'repeat':1,#number of set of replicas
'fold':4}#fold of cross validation
#'testperc':0.00}#,dsearch2,dsearch5,dsearch6,dsearch7 #,'testperc':0.33
#
#test code in local host
if 0:
#sml=[ssm3,ssm1]
#dsearch6,dsearch1,dsearch2,dsearch7,dsearch9
#define the final model
#define the final model
so=scorer(model=convert2old(model1))
#opt=optimizer(scorer=so)
#opt.get_initial_value()
#optres=opt.optimize()
#profile.run('opt.optimize()')
so.assess_model()
pdb.set_trace()
spl=spss(model=model1)
#spl.eval_allpars()
pdb.set_trace()
spl.currentcutoffpercinitial=0.0
spl.currentcutoffpercratio=0.0
spl.maximumround=30
spl.find_best_par()
spl.log()
if 0:#load old result, can be used to detect sampling or other problems
logpath='/bell3/gqdong/statpot/results/mhc2/dcg/runs/57295-nan,nan_nan+nan_nan+nan_2312.403_0.000__nan+nan_nan+nan_0.000_0.000 /'
bestmodel=pickle.load(open(os.path.join(logpath,'bestmodel.pickle'), 'rb'))
for j in range(len(allsearches[i]['object'][allsearches[i]['key']])):
allsearches[i]['object'][allsearches[i]['key']][j]=bestmodel['searches'][k]['object'][allsearches[i]['key']][j]
spl=spss(model=model1) # set of statistical potentials/discrete searches
spl.currentcutoffpercinitial=0.0# no randomness on discrete searches
spl.currentcutoffpercratio=0.0
spl.maximumround=10
#spl.eval_allpars()
spl.find_best_par()#find the best parameters on discrete and continuous variables.
spl.log()
pdb.set_trace()
Run SOAP script to calculate the optimal statistical potential using the best model found in the last step.
from SOAP import *
import os
import profile
import sys
import pdb
import numpy as np
rfeature=[['dt',0,6,0.2]]
rrfeature=[['dt',0,6,6]]
runenv.otherrun=False
par={'uniform':10,'featurebins':rfeature[0],'distribute':'lownumberpriority','firstpoint':2.25}
slo={'key':'sacnfflex',
'valueset':{
'sacnfflex':{'sftype':'sacnf','par':par,'parvalue':[1,1,1,1]},
'psn3':{'sftype':'psn','par':range(1,15,5),'parvalue':range(3)},
'psn4':{'sftype':'psn','par':range(1,15,4),'parvalue':range(4)},
'psn5':{'sftype':'psn','par':range(1,15,3),'parvalue':range(5)},
'psn8':{'sftype':'psn','par':range(0,15,2),'parvalue':range(8)},
}}
sfeature=[rfeature[0],'g6#180','gs6#180','h12#180#-180','t306','ts306']
pm=['','nbsum']
ni=40
initvalues=list(np.arange(0,ni)/10.0+1)
initvalues=np.array(initvalues)
ratio1=[1.0]
#decoyset(dsname='ppd4s',sourcedir='/bell2/gqdong/rawdecoyfiles/ppd/ppd4s/').update_ppdt_pir()
#protein-peptide interactions
ref1={'type':'sf','features':rfeature,'sftype':slo,'par':slo,'parvalue':slo,'ratio':ratio1}
ref2={'type':'sf','features':[rrfeature[0],'g6#180'],'sftype':'spline','par':[0,60,120,180],'parvalue':[1,1,1,1],'ratio':ratio1}
ref3={'type':'sf','features':[rrfeature[0],'gs6#180'],'sftype':'spline','par':[0,60,120,180],'parvalue':[1,1,1,1],'ratio':ratio1}
ref4={'type':'sf','features':[rrfeature[0],'h12#180#-180'],'sftype':'spline','par':[-180,-60,60,180],'parvalue':[1,1,1,1],'ratio':ratio1}
scaledsp1={'type':'scaledsp','pdbset':'X2_2.2A_0.25rfree','features':sfeature,'genmethod':'cs1','pm':pm,'refs':[ref1,ref2,ref3,ref4],'ratio':ratio1}
ratio2=[1.04]
#peptide intra-atom interactions
ref21={'type':'sf','features':'dt30#6','sftype':slo,'par':slo,'parvalue':slo,'ratio':ratio2,'bm':'bs2dsp'}
ref22={'type':'sf','features':'dt1#6g6#180','sftype':'spline','par':[0,60,120,180],'parvalue':[1,1,1,1],'ratio':ratio2,'bm':'bs2dsp'}
ref23={'type':'sf','features':'dt1#6gs6#180','sftype':'spline','par':[0,60,120,180],'parvalue':[1,1,1,1],'ratio':ratio2,'bm':'bs2dsp'}
ref24={'type':'sf','features':'dt1#6h12#180#-180','sftype':'spline','par':[-180,-60,60,180],'parvalue':[1,1,1,1],'ratio':ratio2,'bm':'bs2dsp'}
scaledsp2={'type':'scaledsp','pdbset':'X_2.2A_0.25rfree','features':'dt30#6g6#180gs6#180h12#180#-180t306ts306','genmethod':'bs2dsp','bm':'bs2dsp','pm':pm,'refs':[ref21,ref22,ref23,ref24],'ratio':ratio2}
ref31={'type':'sf','features':'r20','sftype':'bins','par':[],'parvalue':[0]*21,'ratio':[1.0],'bm':''}
#define the benchmark scheme
bmtype={'type':'dsscore','dslist':['mhc2'],'criteria':'dcg','combine':'scoresum','bm':'cs1'}
#top10_rmsd10_irmsd4 10xtop10_rmsd10_irmsd4+10xtop10_rmsd5_irmsd2+ +100xtop10_rmsd5_irmsd2+top10_rmsd+top10_irmsd top10_rmsd10_irmsd4+top20_rmsd+
#top1000_rlrank__rmsdallif1FIRST+top1000_rlrank__rmsdallif2FIRST+top1000_rlrank__rmsdallif3FIRST
#3xtop1_rmsdallif_mean_+2xtop3_rmsdallif_mean_
#define the search positions
search1={'object':ref1,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'dfire','values':initvalues}}
search2={'object':ref1,'key':'par','pos':[0,1,2],'InitialGenerator':{'type':'dfire','values':initvalues}}
#define the search positions
search3={'object':ref2,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'sin','values':np.arange(0,ni)/10.0}}
search4={'object':ref3,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'sin','values':np.arange(0,ni)/10.0}}
search5={'object':ref4,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'sin','values':np.arange(0,ni)/10.0}}
search21={'object':ref21,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'dfire','values':initvalues}}
search22={'object':ref21,'key':'par','pos':[0,1,2,3],'InitialGenerator':{'type':'dfire','values':initvalues}}
#define the search positions
search23={'object':ref22,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'sin','values':np.arange(0,ni)/10.0}}
search24={'object':ref23,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'sin','values':np.arange(0,ni)/10.0}}
search25={'object':ref24,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'sin','values':np.arange(0,ni)/10.0}}
search9={'object':scaledsp2,'key':'ratio','pos':[0],'InitialGenerator':[[0.99] for i in range(40)]}
search31={'object':ref31,'key':'parvalue','pos':[0,1,2,3],'InitialGenerator':{'type':'random'}}
#search32={'object':ref31,'key':'ratio','pos':[0],'InitialGenerator':[[0.99] for i in range(ni)]}
#discrete search options
dsearch1={'object':[scaledsp1,bmtype],'key':['genmethod','bm'],'valueset':[['bs1','bs1'],['bs2','bs2'],['bs3','bs3']
,['bs4','bs4'],['bs5','bs5'],['bs6','bs6'],['bs7','bs7'],['bs8','bs8'],['bs9','bs9'],['bs10','bs10']
,['bs10','bs1'],['bs10','bs2'],['bs10','bs3'],['bs10','bs4']]}
dsearch2={'object':ref1,'key':'features','valueset':['dt800#20','dt200#20','dt40#20','dt600#15','dt300#15','dt150#15','dt30#15','dt100#10']}
dsearch3={'object':ref2,'key':'features','valueset':['g24#180','g12#180','g6#180']}
dsearch4={'object':ref3,'key':'features','valueset':['gs24#180','gs12#180','gs6#180']}
dsearch5={'object':ref4,'key':'features','valueset':['h48#360','h24#360','gs12#360']}
dsearch6={'object':ref1,'key':'par','valueset':[[1,2,8,15],[1,3,8,15],[2,3,8,15],range(1,16,2),range(1,16),range(0,16)]}
dsearch7={'object':ref2,'key':'par','valueset':[[0,60,120,180],range(0,181,30),range(0,181,20),range(0,181,10)]}
dsearch8={'object':ref3,'key':'par','valueset':[[0,60,120,180],range(0,181,30),range(0,181,20),range(0,181,10)]}
#dsearch9={'object':ref4,'key':'par','valueset':[[0,120,240,360],range(0,361,60),range(0,361,30),range(0,361,10)]}
dsearch10={'object':scaledsp1,'key':'features','valueset':['dt30#6g6#180gs6#180h12#180#-180t306ts306','dt12#6g6#180gs6#180h12#180#-180t306ts306']}
dsearch11={'object':scaledsp1,'key':'pdbset','valueset':['X_2.2A','X_2.2A_nopsi']}
dsearch1={'object':[scaledsp1,bmtype],'key':['genmethod','bm'],'valueset':[['bs1dsp','bs1dsp'],['bs2dsp','bs2dsp'],['bs4dsp','bs4dsp'],['bs10dsp','bs10dsp']]}#'bs4dsp','bs5dsp','bs9dsp',,'bs2dsp'
dsearch2={'object':rfeature[0],'key':2,'valueset':[15, 12,10,8,6,5]}#,10,8,7,6,5#15,12,10, 8,7,15,12,10, 8,7,6,
#dsearch3={'object':[rfeature[0],sfeature[0]],'key':[3,3],'valueset':[[0.025,0.025],[0.05,0.05],[0.1,0.1],[0.25,0.25],[0.5,0.5]]}
#dsearch4={'object':sfeature[0],'key':2,'valueset':[20,17,15,12,10,8,6]}
dsearch6={'object':slo,'key':'key','valueset':slo['valueset'].keys()}#slo['valueset'].keys()'psn5','psn3','psn4','psn8'
dsearch7={'object':pm,'key':0,'valueset':['','pp5','pp1','ks0.2']}#'gps0.6','gps0.3',,'ks0.2','ks0.4'
dsearch8={'object':pm,'key':1,'valueset':['nbsum']}
dsearch9={'object':[scaledsp1,scaledsp1],'key':['genmethod','pdbset'],'valueset':[['cs1','X2_2.2A_0.25rfree'],['bs20dsp','X_2.2A_0.25rfree']]}#,,'cs1'['bs20dsp','X_2.2A_0.25rfree'],
dsearch2={'object':[rfeature[0],rrfeature[0],rrfeature[0]],'key':[2,2,3],'valueset':[[6,6,6],[8,8,8]]}#,,[10,10,10],[12,12,12][5.4,5.4,5.4],[4.4,4.4,4.4],10,8,7,6,5#15,12,10, 8,7,15,12,10, 8,7,6,
#[6.4,6.4,6.4],[6.2,6.2,6.2],[5.6,5.6,5.6],,[5,5,5]
#dsearch3={'object':[rfeature[0],sfeature[0]],'key':[3,3],'valueset':[[0.025,0.025],[0.05,0.05],[0.1,0.1],[0.25,0.25],[0.5,0.5]]}
dsearch4={'object':par,'key':'uniform','valueset':[10]} #7,6,5,4,3,[12,10,8,7,6,4]
dsearch5={'object':rfeature[0],'key':3,'valueset':[0.2]}#0.025,0.05 #,0.25,0.1,0.05 #0.05,0.1,0.25,
#dsearch5={'object':rfeature[0],'key':3,'valueset':[0.05,0.1,0.5]}#0.025,0.05 #,0.25,0.1,0.05 #0.05,0.1,0.25,
#dsearch9={'object':bmtype,'key':'bm','valueset':['ss1dsp','bs1dsp','bs2dsp','bs10dsp']}#,
#sampling and seaching parameters
inner={'improved':2,'maxloop':11,'minloop':1}
outer={'improved':4,'maxloop':2001}
td=list(np.sqrt(np.linspace(1,float(10)**2,ni)))
td=np.array(td)
tune_interval=202
sampleschedule={'inner':inner,'outer':outer}
ssm={'sm':'mcp','reset_betweenruns':2,'blockupdate':False, 'exchange_method':1,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':tune_interval,'add_bias':False,'temperature_distribution':td}
ssm0={'sm':'mcs','reset_betweenruns':1,'blockupdate':False,'using_globalbest':True,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp1.0','tune_interval':tune_interval,'add_bias':False,'temperature_distribution':td}
ssm2={'sm':'mcp','reset_betweenruns':2,'blockupdate':False, 'exchange_method':1,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':tune_interval,'add_bias':False,'temperature_distribution':td}
ssm20={'sm':'mcs','reset_betweenruns':2,'blockupdate':False,'using_globalbest':True,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':tune_interval,'add_bias':False,'temperature_distribution':td}
ssm1={'sm':'mca','reset_betweenruns':2,'blockupdate':False, 'using_globalbest':True,
'sample_schedule':sampleschedule,
'stepmethod':'mxmp2.0','tune_interval':tune_interval,'add_bias':False,'temperature_distribution':np.zeros(ni)+2}
ssm3={'sm':'powell','blockupdate':False}
sml=[ssm20,ssm2,ssm0,ssm,ssm0,ssm,ssm0,ssm, ssm1]
#sml=[ssm3,ssm1]
#sml=[ssm0]
dsearch99={'object':[ssm,ssm2],'key':['blockupdate','blockupdate'],'valueset':[[False]*2]}#,
model1={'scorers':[scaledsp1,ref1,ref2,ref3,ref4,scaledsp2,ref21,ref22,ref23,ref24,ref31],'bmtype':bmtype,'searches':[search1,search2,search3,search4,search5,search21,search22,search23,search24,search25,search9,search31], 'runsperscorer':ni,#number of replicas per scorers
#'initialmodelpath':inipath,
'dsearches':[dsearch4],#[dsearch4,dsearch2,dsearch9],,dsearch2,dsearch9
'sml':sml,#search schedule
'cvk':0,#number of repeation of cross validation
'repeat':4,#number of set of replicas
'fold':0,#fold of cross validation
'testperc':0.00}#,dsearch2,dsearch5,dsearch6,dsearch7 #,'testperc':0.33
#
spl=spss(model=model1) # set of statistical potentials/dicrete searches
spl.currentcutoffpercinitial=0.0# no randomness on dicrete searches
spl.currentcutoffpercratio=0.0
spl.maximumround=10
spl.eval_allpars()
#spl.find_best_par()#find the best parameters on discrete and continuous variables.
spl.log()
pdb.set_trace()
Write out the optimal potential in HDF5 or lib format for use in Modeller, IMP or other packages.
"""
Write out the SOAP tables for use in Modeller and IMP.
Also write our raw scores for comparison with scores from IMP or Modeller.
"""
from __future__ import print_function
from SOAP.scorer import *
import pdb
#the path to the optimal potential
logpath=runenv.basedir+'results/mhc2/dcg/runs/57466-nan,nan_nan+nan_nan+nan_2297.581_0.000__nan+nan_nan+nan_0.000_0.000 /'
print(logpath)
bm=pickle.load(open(logpath+'cvmodel.pickle', 'rb'))
so=scorer(model=bm['model'])
#assess the model using the best parameter (necessary for writing the optimal potential)
print(so.assess(bm['allresult'][0]['bestpar'],slevel='dcg'))
print(so.assess_model())
#write out the potential in HDF5 format; the path to the HDF5 files are printed out in log
so.write_potential(filetype='hdf5')
#write out detailed score for comparison with other implementation of SOAP_mhc2
np.save('allscore.npy',so.score)
jdso=pickle.load(open(runenv.basedir+'Decoys/mhc2/mhc2.pickle', 'rb'))
rl=[', '.join([dso.dnlist[i],str(so.score[i]),str(so.scorearray[9,i]*so.ratioarray[9]),str(sum(so.scorearray[4:9,i]*so.ratioarray[4:9])),str(sum(so.scorearray[:4,i])+so.scorearray[10,i])]+map(str,list(so.scorerlist[-1].idist[i]))) for i in range(len(dso.dnlist))]
open('scores','w').write('\n'.join(rl))
print(os.getcwd())
#enter debug mode, so that you can examine the detail scores if necessary
pdb.set_trace()