wiki

#!/usr/bin/env python
import sys,os,time,re,datetime,smtplib

#------PPI_FACES.py---------------------------
#
# This script works with it's partner python, PPI_FACES.sh, and a MATLAB script 
# PPI_FACES.m (in Scripts folder on Munin) to run PPI analysis.  It
# uses the PPI "button" in SPM to first extract a timeseries from
# based on a user specified roi.  It takes in an order number for a group of subjects
# to set the order specific matrices that are used in the PPI analysis.  This is
# the "11 regressor model" that should only be used for Faces contrasts other
# than Faces > Shapes.  If you want Faces > Shapes please use PPI_ALLFACES.
# If you are interested in CARDS please use PPI_CARDS. 
# To minimize user error in entering contrast and order number specific matrices,
# Vanessa decided to hard code the matrices into the .m script, and they are
# set based on the order number variable specified by the user. The script then 
# moves into a single subject analysis, and uses the  extracted values 
# (under PPI.P, PPI.ppi, and PPI.Y) as regressors.  Additionally, 
# it uses the ART motion outliers from the subject's original BOLD 
# preprocessing for the task of interest as an additional regressor.  
# Output goes into the user specified folder under Analysis/SPM/Analyzed/(Subject)
# 

#########user section#########################
#user specific constants
username = "abc1"                               #your cluster login name (use what shows up in qstatall)
useremail = "email@name.com"                 #email to send job notices to
template_f = file("PPI_FACES.sh")            #job template location (on head node)
experiment = "NAME.01"                           #experiment name for qsub
nodes = 400                                     #number of nodes on cluster
maintain_n_jobs = 100                           #leave one in q to keep them moving through
min_jobs = 5                                    #minimum number of jobs to keep in q even when crowded
n_fake_jobs = 50                                #during business hours, pretend there are extra jobs to try and leave a few spots open
sleep_time = 20                                 #pause time (sec) between job count checks
max_run_time = 720                              #maximum time any job is allowed to run in minutes
max_run_hours=24	                        #maximum number of hours submission script can run
warning_time=18                                 #send out a warning after this many hours informing you that the deamon is still running
                                                #make job files  these are the lists to be traversed
                                                #all iterated items must be in "[ ]" separated by commas.
#experiment variables 
subnums = ["22222222_22222"]

runs = [1]                     # [ run01 ] range cuts the last number off any single runs should still be in [ ] or can be runs=range(1,2)

# ppi general variables
task = "faces"                 # MUST be either "faces" or "cards" or "rest" SCRIPT ONLY WORKING WITH FACES CURRENTLY!
outputfolder="mask_design_200"    # This is the name that you would like for the output folder under "Analyzed/Subject/PPI"
                               # You should append the number of subjects used to create the group mask "rAMY_cluster_141"
                               # Since this folder will be in a higher directory called "PPI" - you don't need PPI in the name!
ordernumber ="1"               # This is the order number that will specify the matrix weights for each subject
			       # [ condition#   sess   weight; condition#  sess  weight ] with as many conditions needed
			       # SPM.mat Thresh Variables
connum = "1"                   # Number of the contrast from single subject SPM.mat that we want to extract values from
conname = "One_gr_Two"    # Name of contrast of interest (Faces_gr_Shapes)
threshdesc = "none"            # This is either "FWE" or "none"
thresh = "1"                   # Threshold we want to use when creating VOI (1)

voxextent = "0"                # Voxel extent threshold for when creating VOI
fcontrast = "37"               # Index of F-contrast number used to adjust data GUI default is 0
numsess = "1"                  # Choice for the number of the session (should it be 1?)

# Mask with another contrast?
maskwithother = "no"           # Do you want to mask with another contrast? (yes/no)
otherconnumber = "1"           # If yes to mask_with_other, the contrast #
otherconthresh = "1"           # If yes to mask_with_other, the contrast threshold
otherinclusive = "0"           # If yes to mask_with_other, the type of mask to do (0-inclusive, 1-exclusive)

# VOI Variables
voiname = "mask_voi_200"       # Name that we want our resulting VOI to have (VOI_ is auto. prefixed)
                               # You should append the number of subjects used to create the group mask "rAMY_cluster_141"
voitype = "mask"               # The type of voi to make (box, sphere, or mask) Mask should
                               # be used if the user wants to do a cluster

# Sphere Variables             # ONLY USED IF voitype = "sphere"
spherecenterx ="22"            # Coordinate of center voxel of sphere
spherecentery ="-4"
spherecenterz ="20"
sphereradius = "5"             # Radius of sphere

# Cube Variables               # ONLY USED IF voitype = "box"
boxcenterx ="10"               # Coordinate of center voxel of box
boxcentery ="10"
boxcenterz ="10"
boxdimx = "2"                  # The dimensions of the box
boxdimy = "2"
boxdimz = "2"

# Mask Variables               # ONLY USED IF voitype = "mask"
maskthresh = ".5"              # Threshold for the mask to define the VOI
subincluded ="200s"            # The name of the subject-number-folder under ROI/PPI/Task/ with masks in it ("141s")
maskname ="mask_cluster.img"   # Tbe full name (including extension) of a mask to be used to define VOI
                               # Must be saved to ROI/PPI/Task_Name/#Subjects/
####!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!##########################################################################
 
def daemonize(stdin='/dev/null',stdout='/dev/null',stderr='/dev/null'):
	try:
		#try first fork
		pid=os.fork()
		if pid>0:
			sys.exit(0)
	except OSError, e:
		sys.stderr.write("for #1 failed: (%d) %s\n" % (e.errno,e.strerror))
		sys.exit(1)
	os.chdir("/")
	os.umask(0)
	os.setsid()
	try:
		#try second fork
		pid=os.fork()
		if pid>0:
			sys.exit(0)
	except OSError, e:
			sys.stderr.write("for #2 failed: (%d) %s\n" % (e.errno, e.strerror))
			sys.exit(1)
	for f in sys.stdout, sys.stderr: f.flush()
	si=file(stdin,'r')
	so=file(stdout,'a+')
	se=file(stderr,'a+',0)
	os.dup2(si.fileno(),sys.stdin.fileno())
	os.dup2(so.fileno(),sys.stdout.fileno())
	os.dup2(se.fileno(),sys.stderr.fileno())
	
 
 
start_dir = os.getcwd()
 
daemonize('/dev/null',os.path.join(start_dir,'daemon.log'),os.path.join(start_dir,'daemon.log'))
sys.stdout.close()
os.chdir(start_dir)
temp=time.localtime()
hour,minute,second=temp[3],temp[4],temp[5]
prev_hr=temp[3]
t0=str(hour)+':'+str(minute)+':'+str(second)
log_name=os.path.join(start_dir,'daemon.log')
log=file(log_name,'w')
log.write('Daemon started at %s with pid %d\n' %(t0,os.getpid()))
log.write('To kill this process type "kill %s" at the head node command line\n' % os.getpid())
log.close()
t0=time.time()
master_clock=0
 
#build allowed timedelta
kill_time_limit = datetime.timedelta(minutes=max_run_time)
 
 
def _check_jobs(username, kill_time_limit, n_fake_jobs):
#careful, looks like all vars are global
#see how many jobs we have  in
 
	#set number of jobs to maintain based on time of day.
	cur_time = datetime.datetime.now() #get current time  #time.localtime()  #get current time
	if (cur_time.weekday > 4) | (cur_time.hour < 8) | (cur_time.hour > 17):
		n_other_jobs = 0
	else: #its a weekday, fake an extra 6 jobs to leave 5 nodes open
		n_other_jobs = n_fake_jobs
 
	n_jobs = 0
	
	status = os.popen("qstat  -u '*'")
	status_list = status.readlines()
 
	for line in status_list:
		#are these active or q'd jobs?
		if (line.find(" r ") > -1):
			running = 1
		elif (line.find(" qw ") > -1):   #all following jobs are in queue not running
			running = 0
 
		#if job is mine
		if (line.find(username) > 0) & (line.find("interact.q") < 0):   #name is in the line, not including first spot
			n_jobs = n_jobs + 1
			if running == 1:   #if active job, check how long its been running and delete it if too long
				job_info = line.split()  #get job information
				start_date = job_info[5].split("/")  #split job start date
				start_time = job_info[6].split(":")  #split time from hours:minutes:seconds format
				started = datetime.datetime(int(start_date[2]), int(start_date[0]), int(start_date[1]),
							int(start_time[0]), int(start_time[1]), int(start_time[2]))
				if ((cur_time - started) > kill_time_limit) & (line.find("stalled") == -1):   #if the active job is over max run time, delete it
					os.system("qdel %s" % (job_info[0]))   #delete the run away job
					print("Job %s was deleted because it ran for more than the maximum time." % (job_info[0]))
 
		# if line starts " ###" and isnt an interactive job
		elif bool(re.match( "^\d+", line )) & (line.find("interact") < 0) & (line.find("(Error)") < 0):
			n_other_jobs = n_other_jobs + 1
	return n_jobs, n_other_jobs
		
#make a directory to write job files to and store the start directory
tmp_dir = str(os.getpid())
os.mkdir(tmp_dir)
 
#read in template
template = template_f.read()
template_f.close()
os.chdir(tmp_dir)
 
#for each subject
for subnum in subnums:
	#for each run
	for run in runs:
		
		#Check for changes in user settings
		user_settings=("/home/%s/user_settings.txt") % (username)
		if os.path.isfile(user_settings):
			f=file(user_settings)
			settings=f.readlines()
			f.close()
			for line in settings:
				exec(line)
 
		#define substitutions, make them in template 
		runstr = "%05d" %(run)  
		tmp_job_file = template.replace( "SUB_USEREMAIL_SUB", useremail )
		tmp_job_file = tmp_job_file.replace( "SUB_SUBNUM_SUB", str(subnum) )
		tmp_job_file = tmp_job_file.replace( "SUB_RUNNUM_SUB", str(run) )
		tmp_job_file = tmp_job_file.replace( "SUB_ATASK_SUB", str(task) )
		tmp_job_file = tmp_job_file.replace( "SUB_ORD_SUB", str(ordernumber) )
		tmp_job_file = tmp_job_file.replace( "SUB_NAMEOFVOI_SUB", str(voiname) )
		tmp_job_file = tmp_job_file.replace( "SUB_OUTPUTFOL_SUB", str(outputfolder) )
		tmp_job_file = tmp_job_file.replace( "SUB_CONNUM_SUB", str(connum) )
		tmp_job_file = tmp_job_file.replace( "SUB_CONNAME_SUB", str(conname) )
		tmp_job_file = tmp_job_file.replace( "SUB_THRESH_SUB", str(thresh) )
		tmp_job_file = tmp_job_file.replace( "SUB_DESCOFTHRESH_SUB", str(threshdesc) )
		tmp_job_file = tmp_job_file.replace( "SUB_VOXEXT_SUB", str(voxextent) )
		tmp_job_file = tmp_job_file.replace( "SUB_DATAADJUST_SUB", str(fcontrast) )
		tmp_job_file = tmp_job_file.replace( "SUB_NUMSESS_SUB", str(numsess) )
		tmp_job_file = tmp_job_file.replace( "SUB_OTHERMASK_SUB", str(maskwithother) )
		tmp_job_file = tmp_job_file.replace( "SUB_OTHERMASKTHRESH_SUB", str(otherconthresh) )
		tmp_job_file = tmp_job_file.replace( "SUB_OTHERMASKCON_SUB", str(otherconnumber) )
		tmp_job_file = tmp_job_file.replace( "SUB_OTHERMASKINCLUSIVE_SUB", str(otherinclusive) )
		tmp_job_file = tmp_job_file.replace( "SUB_VMASKTYPE_SUB", str(voitype) )
		tmp_job_file = tmp_job_file.replace( "SUB_CENTERSPHEREX_SUB", str(spherecenterx) )
		tmp_job_file = tmp_job_file.replace( "SUB_CENTERSPHEREY_SUB", str(spherecentery) )
		tmp_job_file = tmp_job_file.replace( "SUB_CENTERSPHEREZ_SUB", str(spherecenterz) )
		tmp_job_file = tmp_job_file.replace( "SUB_RADIUSSPHERE_SUB", str(sphereradius) )
		tmp_job_file = tmp_job_file.replace( "SUB_CENTERBOXX_SUB", str(boxcenterx) )
                tmp_job_file = tmp_job_file.replace( "SUB_CENTERBOXY_SUB", str(boxcentery) )
		tmp_job_file = tmp_job_file.replace( "SUB_CENTERBOXZ_SUB", str(boxcenterz) )
		tmp_job_file = tmp_job_file.replace( "SUB_DIMBOXX_SUB", str(boxdimx) )
		tmp_job_file = tmp_job_file.replace( "SUB_DIMBOXY_SUB", str(boxdimy) )
		tmp_job_file = tmp_job_file.replace( "SUB_DIMBOXZ_SUB", str(boxdimz) )
		tmp_job_file = tmp_job_file.replace( "SUB_SUBINCLUDED_SUB", str(subincluded) )
		tmp_job_file = tmp_job_file.replace( "SUB_THEMASKVOI_SUB", str(maskname) )
		tmp_job_file = tmp_job_file.replace( "SUB_MASKVOITHRESH_SUB", str(maskthresh) )
		
		#make fname and write job file to cwd
		tmp_job_fname = "_".join( ["PPI_FACES", subnum, runstr ] )
		tmp_job_fname = ".".join( [ tmp_job_fname, "job" ] )
		tmp_job_f = file( tmp_job_fname, "w" )
		tmp_job_f.write(tmp_job_file)
		tmp_job_f.close()
 
 
		#wait to submit the job until we have fewer than maintain in q
		n_jobs = maintain_n_jobs
		while n_jobs >= maintain_n_jobs: 
 
			#count jobs
			n_jobs, n_other_jobs = _check_jobs(username, kill_time_limit, n_fake_jobs)   #count jobs, delete jobs that are too old
 
			#adjust job submission by how may jobs are submitted
			#set to minimum number if all nodes are occupied
			#should still try to leave # open on weekdays
			if ((n_other_jobs+ n_jobs) > (nodes+1)): 
				n_jobs = maintain_n_jobs  - (min_jobs - n_jobs)
 
			if n_jobs >= maintain_n_jobs: 
				time.sleep(sleep_time)
			elif n_jobs < maintain_n_jobs:
				cmd = "qsub -v EXPERIMENT=%s %s"  % ( experiment, tmp_job_fname )
				dummy, f = os.popen2(cmd)
				time.sleep(sleep_time)
 
	#Check what how long daemon has been running
	t1=time.time()
	hour=(t1-t0)/3600
	log=file(log_name,'a+')
	log.write('Daemon has been running for %s hours\n' % hour)
	log.close()
	now_hr=time.localtime()[3]
	if now_hr>prev_hr:
		master_clock=master_clock+1
	prev_hr=now_hr
	serverURL="email.biac.duke.edu"
	if master_clock==warning_time:
		headers="From: %s\r\nTo: %s\r\nSubject: Daemon job still running!\r\n\r\n" % (useremail,useremail)
		text="""Your daemon job has been running for %d hours.  It will be killed after %d.
		To kill it now, log onto the head node and type kill %d""" % (warning_time,max_run_hours,os.getpid())
		message=headers+text
		mailServer=smtplib.SMTP(serverURL)
		mailServer.sendmail(useremail,useremail,message)
		mailServer.quit()
	elif master_clock==max_run_hours:
		headers="From: %s\r\nTo: %s\r\nSubject: Daemon job killed!\r\n\r\n" % (useremail,useremail)
		text="Your daemon job has been killed.  It has run for the maximum time alotted"
		message=headers+text
		mailServer=smtplib.SMTP(serverURL)
		mailServer.sendmail(useremail,useremail,message)
		mailServer.quit()
		ID=os.getpid()
		os.system('kill '+str(ID))
 
 
 
#wait for jobs to complete
#delete them if they run too long
n_jobs = 1
while n_jobs > 0:
	n_jobs, n_other_jobs = _check_jobs(username, kill_time_limit, n_fake_jobs)
	time.sleep(sleep_time)
 
 
#remove tmp job files move to start dir and delete tmpdir
#terminated jobs will prevent this from executing
#you will then have to clean up a "#####" directory with
# ".job" files written in it.
cmd = "rm *.job"
os.system(cmd)
os.chdir(start_dir)
os.rmdir(tmp_dir)