/staggbin/getsplitresmask.py

#!/usr/bin/env python

import glob
import os
#from matplotlib import pyplot
import numpy
import frealignlib
import sys
import subprocess

def getnptcls(paramfilename):
	f=open(paramfilename,'r')
	lines=f.readlines()
	f.close()
	nptcls=int(lines[-1].split()[-1])
	return nptcls

def getrange(paramfilename):
	f=open(paramfilename,'r')
	f.readline()
	datarange=f.readline()
	f.close()
	
	words=datarange.split()
	
	if len(words) >= 3:
		start=words[1]
		end=words[2]
	elif len(words) == 2:
		start=words[1]
		end=None

	return start,end

def sortf(x,y):
	xa=int(x.split('_')[-1])
	ya=int(y.split('_')[-1])
	return xa-ya
		
paramfilename=sys.argv[1]
###change this in the future, make as an option
innerdiam=0
outerdiam=0
apix=float(sys.argv[2])
sets=glob.glob('set_*')

basedir=os.getcwd()
fscdict={}
nptcls=[]
respoint5=[]
respoint143=[]
rmeaspoint5=[]
rmeaspoint143=[]
resout=open('p_vs_res.m.txt','w')
resout.write('Start\tStop\tNPtcls\tFSC_0.5\tFSC_0.143\t1/FSC_0.5\t1/FSC_0.143\n')

sets.sort(cmp=sortf)
for set in sets:
	os.chdir(set)
	print set
	ns, res, fsc = frealignlib.parseCombine('frealign.combine.out')
	nptcls.append(getnptcls(paramfilename))
	start,end=getrange(paramfilename)
	fscdict[nptcls[-1]]=(ns,res,fsc)

	
	e='even'
	eroot=e.split('.')[0]
	eout=eroot+'.m.mrc'
	o='odd'
	oroot=o.split('.')[0]
	oout=oroot+'.m.mrc'
	
	ecommand='e2proc3d.py --process=mask.auto3d:radius=90:nshells=7:threshold=0.002:nshellsgauss=3 %s %s' % (e,eout)
	ocommand='e2proc3d.py --process=mask.auto3d:radius=90:nshells=7:threshold=0.002:nshellsgauss=3 %s %s' % (o,oout)
	
	print ecommand
	os.system(ecommand)
	print ocommand
	os.system(ocommand)
	
	fsc='%s_%s.eotest' % (eout, oout)
	fsccommand='proc3d %s %s fsc=%s' % (eout, oout, fsc)
	print fsccommand
	os.system(fsccommand)
	
	f=open(fsc)
	lines=f.readlines()
	f.close()
	
	freq=[]
	corr=[]
	for line in lines:
		words=line.split()
		freq.append(int(words[0]))
		corr.append(float(words[1]))
	
	for n in range(len(freq)):
		if corr[n] < 0.5:
			freq0_5=(freq[n]+freq[n-1])/2.0
			res0_5=len(freq)*2*apix/freq0_5
			respoint5.append(res0_5)
			break
			
	for n in range(len(freq)):
		if corr[n] < 0.143:
			freq0_143=(freq[n]+freq[n-1])/2.0
			res0_143=len(freq)*2*apix/freq0_143
			respoint143.append(res0_143)
			break
	

	resout.write('%s\t%s\t%d\t%f\t%f\t%f\t%f\n'% (start,end,nptcls[-1],respoint5[-1],respoint143[-1],1/respoint5[-1],1/respoint143[-1]))
	#os.system('rm working.mrc')
	#os.system('rm even')
	os.system('rm even.m.mrc')
	#os.system('rm odd')
	os.system('rm odd.m.mrc')
	os.system('rm weights')
	os.chdir(basedir)

resout.close()