Source code for SimUtils.CC3DPipeline

import os,gzip,string,time,shutil
import numpy as np
import Image
#~ from CC3DTools.CC3DReaders import *
from .AnalysisUtils import *
from .ImageUtils import *
from .Readers import *
from .Cell import *
from multiprocessing import Pool
#from Compactness import getCompactness
#from DrawCells import drawCells

## @package SimUtils.CC3DPipeline 
# Utility functions for a simulation pipeline (as used on the cluster).

[docs]def makeImages(id,trange,inpath,outpath,cm='default.ctb',gzipped=False,timestamp=False,label=False,scale=1,bc=None,fontsize=6,fieldname=None,border=True): """ Make images for a single simulation simulation :param id: simulation identifier :type id: str :param trange: list of mcs to draw :param inpath: path containing data files :type inpath: str :param outpath: path to save images to :type outpath: str :param cm: file containing the colormap :type cm: str :param gzipped: data is gzipped :type gzipped: bool :param timestamp: add time stamp to the image :type timestamp: bool :param label: add id as label to the image :type label: bool :param scale: scaling of the image :type scale: number :param bc: color of cell boundaries (r,g,b) :param fontsize: size of the fonts used for label and time stamp; font size will be multiplied by scale. :type fontsize: int :param fieldname: name of chemical field :type fieldname: str :param border: cut of border pixels :type border: bool """ t0 = time.time() tlen = len(str(trange[-1])) colormap = readColorMap(cm) for t in trange: outfile = id+'_'+string.zfill(str(t),tlen)+'.png' im = makeImage(id,inpath,t,colormap,timestamp,label,scale,bc,fontsize,border,gzipped,fieldname) im.save(outpath+outfile) print str(len(trange)) + ' images drawn in '+str(time.time()-t0)+' seconds'
def makeImagesMC(id,trange,inpath,outpath,cm='default.ctb',gzipped=False,timestamp=False,label=False,scale=1,bc=None,fontsize=6,fieldname=None,border=True,ncores=8): """ Make images for a single simulation simulation, using multiple processors :param id: simulation identifier :type id: str :param trange: list of mcs to draw :param inpath: path containing data files :type inpath: str :param outpath: path to save images to :type outpath: str :param cm: file containing the colormap :type cm: str :param gzipped: data is gzipped :type gzipped: bool :param timestamp: add time stamp to the image :type timestamp: bool :param label: add id as label to the image :type label: bool :param scale: scaling of the image :type scale: number :param bc: color of cell boundaries (r,g,b) :param fontsize: size of the fonts used for label and time stamp; font size will be multiplied by scale. :type fontsize: int :param fieldname: name of chemical field :type fieldname: str :param border: cut of border pixels :type border: bool :param ncores: number of cores that may be used :type ncores: int """ #~ print os.getcwd() t0 = time.time() tlen = len(str(trange[-1])) colormap = readColorMap(cm) jobs = [[id,inpath,t,postfix,colormap,outpath+'/'+id+'_'+string.zfill(str(t),tlen)+'.png', \ timestamp,label,scale,bc,fontsize,border,gzipped,fieldname] for t in trange] pool = Pool(processes=ncores) pool.map(makeImageMC,jobs) print str(len(trange)) + ' images drawn in '+str(time.time()-t0)+' seconds' def makeImageMC(args): """ Wrapper function for :func:`~SimUtils.CC3DPipeline.makeImage` to work with multiprocessing :param args: list of all arguments of :func:`~SimUtils.CC3DPipeline.makeImage` """ [id,inpath,t,colormap,outfile,timestamp,label,scale,bc,fontsize,border,gzipped,fieldname] = args im = makeImage(id,inpath,t,colormap,timestamp,label,scale,bc,fontsize,border,gzipped,fieldname) im.save(outfile) def makeImage(id,inpath,t,colormap,timestamp=False,label=False,scale=1,bc=None,fontsize=6,border=True,gzipped=True,fieldname=None): """ Draw morphology for one timestep :param id: simulation identifier :param inpath: path containing data files :param t: time step :param outpath: path to save images to :param colormap: dictionary with cell types as keys and colors (r,g,b) as values :param timestamp: add time stamp to the image :type timestamp: bool :param label: add id as label to the image :type label: bool :param scale: scaling of the image :type scale: number :param bc: color of cell boundaries (r,g,b) :param fontsize: size of the fonts used for label and time stamp; font size will be multiplied by scale. :type fontsize: int :param border: cut of border pixels :type border: bool :param gzipped: data is gzipped :type gzipped: bool :param fieldname: name of chemical field :type fieldname: str :return: image object .. seealso:: :func:`~SimUtils.ImageUtils.drawCells`, :func:`~SimUtils.ImageUtils.addTimeStamp`, :func:`~SimUtils.ImageUtils.addLabel` """ sigma = readSigma(id,t,inpath,gzipped,border) types = readTau(id,t,inpath,gzipped,border) if fieldname is not None: field = readChemField(id,t,inpath,fieldname,gzipped,border) else: field = np.zeros_like(sigma) (nx,ny) = sigma.shape im = drawCells(sigma,types,field,colormap,scale=2,nx=nx,ny=ny,bc=bc) if timestamp: fs = int(fontsize*(nx/200.0)*scale) addTimeStamp(im,str(t),fs,fc=bc) if label: fs = int(1.25*fontsize*(nx/200.0)*scale) addLabel(im,str(id),fs,fc=bc) return im def getNodesAndDrawOnMorph(id,trange,inpath,outpath,closedisk=5,thinloop=10,d=6,cm='default.ctb',gzipped=False): """ Calculate nodes for a simulations and draw nodes on the morphology :param id: simulation identifier :param trange: list of time steps :param inpath: path containing data files :param outpath: path to save images to :param closedisk: size of the disk used for closing. :param thinloop: number of thinning iterations. :param d: minimal distance between two branchpoints. :param cm file: containing the colormap :param gzipped: if True, data is expected to be gzipped, and stored in inpath/id/ .. seealso:: :func:`~SimUtils.AnalysisUtils.getNodes`, :func:`~SimUtils.AnalysisUtils.drawNodesColored` """ t0 = time.time() tlen = len(str(trange[-1])) colormap = readColorMap(cm) postfix = '.data' if gzipped: postfix += '.gz' inpath += '/'+id+'/' for t in trange: im = np.loadtxt(inpath+'/'+id+'_TF_'+str(t)+postfix) im = im>0 im = im[1:im.shape[0]-1,1:im.shape[1]-1] nodes,skel = getNodes(im,closedisk=closedisk,thinloop=thinloop) endpoints = getEndPoints(skel,nodes) for ep in endpoints: nodes.append(ep) im = makeImage(id,inpath,t,postfix,colormap,scale=2) im = im.rotate(-90) colormapNodes = dict((i,(0,0,0)) for i in range(len(nodes))) nodedict = dict((i,node) for i,node in enumerate(nodes)) #~ colormapNodes = {i:(0,0,0) for i in range(len(nodes))} #~ nodedict = {i:node for i,node in enumerate(nodes)} outfile = id+'_'+string.zfill(str(t),tlen) drawNodesColored(nodedict,im,outpath+'/'+outfile,colormapNodes,scale=2,r=2) #~ im.save(outpath+outfile) print str(len(trange)) + ' images drawn in '+str(time.time()-t0)+' seconds' def getNodesAndSkelOnMorph(id,trange,inpath,outpath,closedisk=5,thinloop=10,d=6,cm='default.ctb',gzipped=False): """ Calculate nodes for a simulations and draw nodes on the morphology together with the skeleton. :param id: simulation identifier :param trange: list of time steps :param inpath: path containing data files :param outpath: path to save images to :param closedisk: size of the disk used for closing. :param thinloop: number of thinning iterations. :param d: minimal distance between two branchpoints. :param cm file: containing the colormap :param gzipped: if True, data is expected to be gzipped, and stored in inpath/id/ .. seealso:: :func:`~SimUtils.AnalysisUtils.getNodes`, :func:`~SimUtils.AnalysisUtils.drawNodesColored` """ t0 = time.time() tlen = len(str(trange[-1])) colormap = readColorMap(cm) postfix = '.data' if gzipped: postfix += '.gz' inpath += '/'+id+'/' for t in trange: im = np.loadtxt(inpath+'/'+id+'_TF_'+str(t)+postfix) im = im>0 im = im[1:im.shape[0]-1,1:im.shape[1]-1] nodes,skel = getNodes(im,closedisk=closedisk,thinloop=thinloop,d=d) endpoints = getEndPoints(skel,nodes) for ep in endpoints: nodes.append(ep) edges = getEdges(nodes,skel) nnodes = copy.deepcopy(nodes) con = dict((i,0) for i,n in enumerate(nodes)) for e in edges: con[e[0]] += 1 colormapNodes = dict((i,(0,0,0)) for i in range(len(nodes))) for n in con: if con[n] < 2: colormapNodes[n] = (150,150,150) if con[n] == 2: colormapNodes[n] = (0,255,0) nodes = nnodes im = makeImage(id,inpath,t,postfix,colormap,scale=2,bc=(180,180,180)) im = im.rotate(-90) nodedict = dict((i,node) for i,node in enumerate(nodes)) #~ colormapNodes = {i:(0,0,0) for i in range(len(nodes))} #~ nodedict = {i:node for i,node in enumerate(nodes)} outfile = id+'_'+string.zfill(str(t),tlen) im = drawSkelOnMorph(skel,im,None,scale=2,save=False) drawNodesColored(nodedict,im,outpath+'/'+outfile,colormapNodes,scale=2,r=2) #~ im.save(outpath+outfile) print str(len(trange)) + ' images drawn in '+str(time.time()-t0)+' seconds' def getConnectedEdgesOnMorph(id,trange,inpath,outpath,closedisk=5,thinloop=10,d=6,cm='default.ctb',gzipped=False): """ Calculate nodes for a simulations and draw nodes on the morphology :param id: simulation identifier :param trange: list of time steps :param inpath: path containing data files :param outpath: path to save images to :param closedisk: size of the disk used for closing. :param thinloop: number of thinning iterations. :param d: minimal distance between two branchpoints. :param cm file: containing the colormap :param gzipped: if True, data is expected to be gzipped, and stored in inpath/id/ .. seealso:: :func:`~SimUtils.AnalysisUtils.getNodes`, :func:`~SimUtils.AnalysisUtils.drawNodesColored` """ t0 = time.time() tlen = len(str(trange[-1])) colormap = readColorMap(cm) postfix = '.data' if gzipped: postfix += '.gz' inpath += '/'+id+'/' for t in trange: im = np.loadtxt(inpath+'/'+id+'_TF_'+str(t)+postfix) im = im>0 im = im[1:im.shape[0]-1,1:im.shape[1]-1] nodes,skel = getNodes(im,closedisk=closedisk,thinloop=thinloop,d=d) endpoints = getEndPoints(skel,nodes) for ep in endpoints: nodes.append(ep) nodes = mergeNodes(nodes,d) con = getNofEdges(nodelist,skel,w) for n,node in enumerate(copy.deepcopy(nodes)): if con[n] == 2: nodes.remove(node) colormapNodes = dict((i,(0,0,0)) for i in range(len(nodes))) im = makeImage(id,inpath,t,postfix,colormap,scale=2,bc=(180,180,180)) im = im.rotate(-90) im = drawSkelOnMorph(skel,im,None,scale=2,save=False) nodedict = dict((i,node) for i,node in enumerate(nodes)) outfile = id+'_'+string.zfill(str(t),tlen) drawNodesColored(nodedict,im,outpath+'/'+outfile,colormapNodes,scale=2,r=2) print str(len(trange)) + ' images drawn in '+str(time.time()-t0)+' seconds'
[docs]def getCompactnessForSim(id,trange,inpath,outpath=None,gzipped=False,border=True): """ Calculate compactness for one simulation :param id: simulation identifier :param trange: range of mcs to calculate :param inpath: path containing data files :param gzipped: if True, data is expected to be gzipped, and stored in inpath/id/ .. seealso:: :func:`~SimUtils.AnalysisUtils.getCompactness` """ t0 = time.time() if outpath is None: outpath = inpath if gzipped: f = gzip.open(outpath+'/'+id+'/'+id+'_compactness.data.gz','w') else: f = open(outpath+'/'+id+'_compactness.data','w') data = np.zeros((len(trange),2)) for i,t in enumerate(trange): sigma = readSigma(id,t,inpath,gzipped,border) data[i,:] = [t,getCompactness(getLCC(sigma))] f.write('#MCS\tcompactness\n') np.savetxt(f,data) f.close() print 'compactness calculated for '+str(len(trange))+' morphologies in '+str(time.time()-t0)+' seconds'
def getNodesAndEdges(sigma,merge=True,purge=False,closedisk=3,thinloop=8,d=10,w=2,nodedist=10): """ Calculate nodes and edges for single morphology :param sigma: CPM grid :param merge: merge nodes that are too close :type merge: bool :param purge: purge nodes that connect two nodes :type purge: bool :param closedisk: size of the disk used for closing. :param thinloop: number of thinning iterations. :param d: minimal distance between two branchpoints (used in node detection) :param w: node width that will be cut out of the first image :param nodedist: minimal distance betweeen to nodes (used in merging of nodes) :return: dictionary with nodes {nodeid : (x,y)}, list of edges [(source,target)], skeleton image .. seealso:: :func:`~SimUtils.AnalysisUtils.getNodes`, :func:`~SimUtils.AnalysisUtils.getEndPoints`, :func:`~SimUtils.AnalysisUtils.getEdges`, :func:`~SimUtils.AnalysisUtils.mergeNodesAndEdges`, :func:`~SimUtils.AnalysisUtils.purgeNodes`, """ nodes,skel = getNodes(sigma>0,closedisk=closedisk,thinloop=thinloop,d=d) endpoints = getEndPoints(skel,nodes,d=d) for ep in endpoints: nodes[max(nodes.keys())+1] = ep #~ nodes.append(ep) edges = getEdges(nodes,skel,w=w) if merge: nodes,edges = mergeNodesAndEdges(nodes,edges,d=nodedist) if purge: nodes,edges = purgeNodes(nodes,edges) return nodes,edges,skel def getTypesOnTipsForSim(id,trange,inpath,typenames,gzipped=False,border=True): """ Count occurance of cell typs on sprout end points for a simulation :param id: simulation id :param trange: list of MCS :param inpath: path to simulation data :param typename: dict with {typeID : typename} :param gzipped: data is gzipped :type gzipped: bool :param border: CPM field is surrounded by a 1 pixel thick border cell :type border: bool .. seealso:: :func:`~SimUtils.CC3DPipeline.getNodesAndEdges`, :class:`~SimUtils.AnalysisUtils.TipFinder` """ if gzipped: f = gzip.open(inpath+'/'+id+'/'+id+'_TIPS.data.gz','w') else: f = open(inpath+'/'+id+'_TIPS.data','w') f.write('#MCS\t'+'\t'.join([str(typenames[i]) for i in typenames])) for t in trange: #~ print t sigma = readSigma(id,t,inpath,gzipped=gzipped,border=border) tau = readTau(id,t,inpath,gzipped=gzipped,border=border) nodes,edges,skel = getNodesAndEdges(sigma,closedisk=3,thinloop=8,d=5) tf = TipFinder(sigma,tau,nodes,edges) tips = tf.findAllCellsOnTips() tiptypes = tips.values() f.write('\n'+str(t)+'\t'+'\t'.join([str(tiptypes.count(type)) for type in typenames])) f.close() def doMorphometry(id,trange,inpath,gzipped=False,save=True,closedisk=3,thinloop=8,d=10): """ Calculate morphometrics with default settings for one simulation. All morphometrics are collected and save in a file inpath/id_morph.data.gz :param id: simulation identifier :param trange: range of mcs to calculate :param inpath: path containing data files :param gzipped: if True, data is expected to be gzipped, and stored in inpath/id/ .. seealso:: :func:`~SimUtils.CC3DPipeline.getNodesAndEdges` """ t0 = time.time() data = np.zeros((len(trange),9)) for i,t in enumerate(trange): if gzipped: im = np.loadtxt(inpath+'/'+id+'/'+id+'_TF_'+str(t)+'.data.gz') else: im = np.loadtxt(inpath+'/'+id+'_TF_'+str(t)+'.data') im = im>0 im = im[1:im.shape[0]-1,1:im.shape[1]-1] c = getCompactness(getLCC(im)) meshes = getMeshes(im) nodes,edges,skel = getNodesAndEdges(im) degree = np.bincount((np.array(edges)).flatten()) nEP = np.sum(degree==1) nBP = np.sum(degree>2) nEdges = len(edges) if len(meshes) > 0: data[i,:] = [t,c,len(nodes),nEdges,len(meshes),np.mean(meshes),np.std(meshes),nEP,nBP] else: data[i,:] = [t,c,len(nodes),nEdges,0,0,0,nEP,nBP] if save: if gzipped: f = gzip.open(inpath+'/'+id+'/'+id+'_morph.data.gz','w') else: f = gzip.open(inpath+'/'+id+'_morph.data.gz','w') f.write('#MCS\tcompactness\tnnodes\tnnedges\tnmeshes\tmeshAavg\tmeshAstd\tnEP\tnBP\n') np.savetxt(f,data) f.close() else: print '#MCS\tcompactness\tnnodes\tnnedges\tnmeshes\tmeshAavg\tmeshAstd\tnEP\tnBP\n' return data print 'mophometrics calculated for '+str(len(trange))+' morphologies in '+str(time.time()-t0)+' seconds' def calcOrderParameter(id,trange,inpath,radii,border=True,gzipped=False,weighted=False,outpath=None): """ Function to call :func:`~SimUtils.CC3DPipeline.getOrderParameterForSim`, needed for compatibility """ getOrderParameterForSim(id,trange,inpath,radii,border,gzipped,weighted,outpath)
[docs]def getOrderParameterForSim(id,trange,inpath,radii,border=True,gzipped=False,weighted=False,outpath=None): """ Calculate orderparameters for one simulation. All order parameters are collected and save in a file inpath/id_orderparameter.data :param id: simulation identifier :param trange: range of mcs to calculate :param inpath: path containing data files :param radii: list of radii for wich the order parameter is calculates :param border: remove border pixels from data :param gzipped: if True, data is expected to be gzipped, and stored in inpath/id/, and the output file will be gzipped. .. seealso:: :func:`~SimUtils.AnalysisUtils.getOrderParameter` """ if outpath is None: outpath = inpath t0 = time.time() if gzipped: f = gzip.open(outpath+'/'+id+'/'+id+'_orderparameter.data.gz','w') else: f = open(outpath+'/'+id+'_orderparameter.data','w') data = np.zeros((len(trange),len(radii)+1)) f.write('#MCS') for r in radii: f.write('\t'+str(r)) f.write('\n') for i,t in enumerate(trange): sigma = readSigma(id,t,inpath,gzipped,border) (nx,ny) = sigma.shape data[i,0] = t angles = getAngleField(sigma) for j,r in enumerate(radii): #~ data[i,j+1] = getOrderParameter(cells,sigma,r,weighted) data[i,j+1] = getOrderParameterFromGrid(sigma,angles,r,weighted) np.savetxt(f,data) f.close() print 'Order parameter calculated for '+str(len(trange))+' simulations and '+str(len(radii))+' radii in '+str(time.time()-t0)+' seconds'
def calcOrderParameterExtra(id,inpath,radii,border=True): """ Calculate extra orderparameters for one simulation. All order parameters are collected and added to a file inpath/id_orderparameter.data.gz. This function expects all data to be gzipped. :param id: simulation identifier :param trange: range of mcs to calculate :param inpath: path containing data files :param radii: list of radii for wich the order parameter is calculates :param border: remove border pixels from data .. seealso:: :func:`~SimUtils.AnalysisUtils.getOrderParameter` """ t0 = time.time() # read existing data dataE = np.loadtxt(inpath+'/'+id+'_orderparameter.data.gz',skiprows=1) f = gzip.open(inpath+'/'+id+'_orderparameter.data.gz','r') head = f.readline() f.close() trange = dataE[:,0] data = np.zeros((len(trange),len(radii)+dataE.shape[1])) data[:,0:dataE.shape[1]] = dataE co = dataE.shape[1] f = gzip.open(inpath+'/'+id+'_orderparameter.data.gz','w') f.write(head.rstrip('\n')) for r in radii: f.write('\t'+str(r)) f.write('\n') for i,t in enumerate(trange): #~ print t cells = readCellDictFromSigma(inpath+'/'+id+'_CF_'+str(int(t))+'.data.gz',inpath+'/'+id+'_TF_'+str(int(t))+'.data.gz',ignore=[1]) sigma = np.loadtxt(inpath+'/'+id+'_CF_'+str(int(t))+'.data.gz') (nx,ny) = sigma.shape if border: sigma[0:nx,0] = 0 sigma[0:nx,ny-1] = 0 sigma[0,0:ny] = 0 sigma[nx-1,0:ny] = 0 #~ data[i,0] = t for j,r in enumerate(radii): data[i,j+co] = getOrderParameter(cells,sigma,r) np.savetxt(f,data) f.close() print 'Order parameter calculated for '+str(len(trange))+' and '+str(len(radii))+' radii in '+str(time.time()-t0)+' seconds' ## Copy simulation data to scratch # Copy and unpack data in scratch to resume post-processing on cluster. # @param id simulation identifier # @param source source folder def copyDataBack(id,source,gunzip=False): """ Copy simulation data back to the scratch drive :param id: simulation id :type id: str :param source: directory with data :type source: str :param gunzip: extract gzipped data :type gunzip: bool """ wd = os.getcwd() # copy data from home to node os.chdir('/scratch/') os.system('cp '+source+'/data_'+id+'.tar ./') os.system('tar -xf data_'+id+'.tar') # extract gzipped data and move to scratch if gunzip: os.chdir('/scratch/'+id) os.system('gunzip *.data.gz') os.system('mv *.data /scratch/') os.chdir('/scratch/') os.system('rm -r '+id) # change working directory to original working directory os.chdir(wd) ## Copy simulation data to home # Copy gzipped simulation data to home before analysis etc. # @param id simulation identifier # @param target target folder def copyData(id,target): """ Gzip data and copy back to home to secure the data before post-processing. :param id: simulation id :param target: directory to copy data to """ wd = os.getcwd() os.chdir('/scratch/') os.mkdir(id) os.system('cp '+id+'_*.data '+id) os.chdir('/scratch/'+id) print os.getcwd() print 'gzip -q '+id+'*.data' os.system('gzip -q '+id+'*.data') os.chdir('/scratch/') print os.getcwd() print 'tar -cf '+target+'/data_'+id+'.tar '+id os.system('tar -cf '+target+'/data_'+id+'.tar '+id) os.system('rm -r '+id) os.chdir(wd) print os.getcwd() ## Move simulation data to home # Text data (txt, data) is gzipped and tarred, images are stored in a gzipped tarball. # @param id simulation identifier # @param target target folder # @param data copy .data files to target folder # @param images copy .png files to target folder # @param txt copy .txt files to target folder def moveData(id,target,data=False,images=False,txt=False): """ Gzip simulation results and copy to home. Different types of data are packed in different ways: * data: all .data files are gzipped and moved to a folder `id`, the folder is tarred on moved to `target`. * images: all .png files are moved to a folder `id`, the folder is gzipped and moved to `target`. * txt: all .txt files are gzipped and moved to a folder `id`, the folder is tarred on moved to `target`. :param id: simulation id :param target: target directory :param data: copy data to target :param images: copy images to target :param txt: copy txt files to target """ wd = os.getcwd() os.chdir('/scratch/') if data: if not os.path.isdir(id): os.makedirs(id) print 'gzip -q '+id+'_*.data' os.system('gzip -q '+id+'_*.data') print 'mv '+id+'_*.data.gz '+id os.system('mv '+id+'_*.data.gz '+id) print 'tar -cf '+target+'/data_'+id+'.tar '+id os.system('tar -cf '+target+'/data_'+id+'.tar '+id) os.system('rm -r '+id) if images: os.mkdir(id) print 'mv '+id+'_*.png '+id os.system('mv '+id+'_*.png '+id) print 'tar -zcf '+target+'/images_'+id+'.tar.gz '+id os.system('tar -zcf '+target+'/images_'+id+'.tar.gz '+id) os.system('rm -r '+id) if txt: os.mkdir(id) print 'gzip -q '+id+'_*.txt' os.system('gzip -q '+id+'_*.txt') print 'mv '+id+'_*.txt.gz '+id os.system('mv '+id+'_*.txt.gz '+id) print 'tar -cf '+target+'/txt_'+id+'.tar '+id os.system('tar -cf '+target+'/txt_'+id+'.tar '+id) os.system('rm -r '+id) os.chdir(wd) def getPIFfromData(id,t,inpath,piffile,typemap,gzipped=False): """ Generate PIF from data file :param id: simulation id :param t: simulation time step :param inpath: path to data files :param piffile: file to save pif to :param typemap: map of cell type names and numbers {typeID : typeName} :param gzipped: data is gzipped """ # load cellfield (sigma) postfix = '.data' if gzipped: postfix += '.gz' sigma = np.loadtxt(inpath+'/'+id+'_CF_'+str(t)+postfix,dtype=np.int) types = np.loadtxt(inpath+'/'+id+'_TF_'+str(t)+postfix,dtype=np.int) # open piffile f = open(piffile,'w') # iterate over grid and write pif for i in range(0,np.shape(sigma)[1]): for j in range(0,np.shape(sigma)[0]): if sigma[j][i] > 0: id = str(sigma[j][i]) type = typemap[types[j][i]] #~ x = str(sigma.shape[0]-i-1) x = str(i) y = str(sigma.shape[1]-j-1) #~ y = str(sigma.shape[1]-j-1) f.write(' '.join([id,type,x,x,y,y])+' 0 0\n') #~ f.write(str(sigma[j][i])+' '+str(typemap[types[j][i]])+' '+str(i)+' '+str(i)+' '+str(j)+' '+str(j)+ ' 0 0\n') f.close() def calculateClusters(simid,trange,indir,outdir,r,th=10,minlabsize=50,opendisk=1,mincellsize=.25,gzipped=False,border=False): """ Calculate clusters for a simulation. :param simid: simulation identifier :param trange: list of timesteps :param indir: data directory :param outdir: data to save output to :param r: radius used for the relative director (:func:`~SimUtils.AnalysisUtils.calculateClustersForSim`) :param th: threshold used for the cluster detection algorithm (func:`~SimUtils:AnalysisUtils:getCellClusters`) :param minlabsize: minimum size of an area to be used for cluster detection (func:`~SimUtils:AnalysisUtils:getCellClusters`) :param opendisk: disk size used for opening during cluster detection algorithm (func:`~SimUtils:AnalysisUtils:getCellClusters`) :param mincellsize: minumum fraction of a cell that must bo in a cluster to be added to that cluster func:`~SimUtils:AnalysisUtils:getCellClusters`) .. seealso:: :func:`~SimUtils.CC3DTools.calculateClustersForSim`, :func:`~SimUtils.AnalysisUtils.getRelativeDirField`, :func:`~SimUtils.AnalysisUtils.getCellClusters`, :func:`~SimUtils.AnalysisUtils.calcClusterInterfaces` """ calculateClustersForSim([simid,trange,indir,outdir,r,th,minlabsize,opendisk,mincellsize,gzipped]) def saveTrajectories(simid,trange,indir,outdir,gzipped=False,border=False,tdict=None): if tdict is None: tdict = getTrajectories(simid,trange,indir,gzipped=gzipped,border=border) data = np.zeros((len(trange),2*len(tdict)+1)) header = '#MCS' data[:,0] = trange for i,id in enumerate(tdict.keys()): data[:,2*i+1:2*i+3] = tdict[id] header += '\t'+str(id)+' X\t'+str(id)+' Y' f = open(outdir+simid+'_traj.data','w') f.write(header+'\n') np.savetxt(f,data) f.close() def drawTrajectoriesForSim(simid,trange,indir,outdir,cmtypes,gzipped=False,border=False,tdict=None,bg=(255,255,255),scale=5,fontpath='/usr/share/fonts/liberation/LiberationSans-Regular.ttf'): sigma = readSigma(simid,trange[0],indir,gzipped,border) tau = readTau(simid,trange[0],indir,gzipped,border) cm = dict((id,cmtypes[np.unique(tau[np.where(sigma==id)])[0]]) for id in np.unique(sigma)) if tdict is None: tdict = getTrajectories(simid,trange,indir,gzipped=gzipped,border=border) drawTrajectories(tdict,outdir+simid+'_traj.png',cm,scale=scale,crop=True,offset=10,bg=bg) def getTrajectories(simid,trange,indir,gzipped=False,border=False): """ Calculate trajactories for all cells in a simulation :param simid: simulation identifier :param trange: list of timesteps :param indir: data directory :param gzipped: set true for gzipped data :param border: set true if data contains a border cell .. seealso:: :func:`~SimUtils.AnalysisUtils.getCoM` """ tdict = {} pool = Pool(processes=7) traj = pool.map(_getCOMforT, [[simid,t,indir,gzipped,border] for t in trange]) for i,t in enumerate(trange): if traj[i][0] != t: print 'this should not happen' exit() for id,com in traj[i][1].iteritems(): tdict.setdefault(id,np.zeros((len(trange),2))) tdict[id][i] = com return tdict def _getCOMforT(args): [simid,t,indir,gzipped,border] = args com = {} sigma = readSigma(simid,t,indir,gzipped=gzipped,border=border) for i in np.unique(sigma): if i == 0: continue com[i] = getCoM(np.where(sigma==i)) return (t,com) #~ def calculateClustersForSim(args): #~ """ Function that actually calculates and save data. This function is designed so it works well with the :mod:`multiprocess` module; instead of multiple arguments all arguments are presented in one list. For user convinience, :func:`~SimUtils.CC3DTools.calculateClusters` is created that uses seperate arguments. #~ """ #~ [simid,trange,indir,outdir,r,th,minlabsize,opendisk,mincellsize,gzipped] = args #~ cellTCdict = {} #~ for t in trange: #~ print simid,t #~ sigma = readSigma(simid,t,indir,gzipped=gzipped) #~ cells = dict((cell,ClusterCell(cell,np.where(sigma==cell))) for cell in np.unique(sigma[sigma>0])) #~ cells = {cell : ClusterCell(cell,np.where(sigma==cell)) for cell in np.unique(sigma[sigma>0])} #~ diffield = getRelativeDirField(sigma,cells,r) #~ clusters = getCellClusters(diffield,sigma,th,minlabsize,opendisk,mincellsize) #~ for id,c in clusters.iteritems(): #~ for cid in c.cells: #~ cells[cid].cluster = id #~ calcClusterInterfaces(sigma,clusters,cells) #~ for cid,cell in cells.iteritems(): #~ if cell.cluster in clusters: #~ cluster = clusters[cell.cluster] #~ else: #~ cluster = None #~ cellTCdict.setdefault(cid,ClusterCellTC(cid)).addTimeStep(t,cell.com,cell.l,cluster) #~ sdata = None #~ for cid,cellTC in cellTCdict.iteritems(): #~ savedata = cellTC.getData() #~ if sdata is None: #~ sdata = savedata #~ else: #~ sdata = np.row_stack((sdata,savedata)) #~ header = cellTC.header #~ fmt = '%d\t%d\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%d\t%d\t%d\t%.1f\t%.1f' #~ f = open(outdir+simid+'_clusters.data','w') #~ f.write(header) #~ np.savetxt(f,sdata,fmt=fmt) #~ f.close() #~ print 'saved data for '+simid def calcClusters(id,trange,indir,outdir,r,th,minlabsize,opendisk,mincellsize,gzipped=False,border=False): """ Function to call :func:`~SimUtils.CC3DPipeline.getClustersForSim`, needed for compatibility """ getClustersForSim(id,trange,indir,outdir,r,th,minlabsize,opendisk,mincellsize,gzipped,border)
[docs]def getClustersForSim(id,trange,indir,outdir,r,th,minlabsize,opendisk,mincellsize,gzipped=False,border=False): cellTCdict = {} t0 = time.time() for t in trange: #~ print 'find clusters for mcs =',t sigma = readSigma(id,t,indir,gzipped=gzipped,border=border) #~ sigma = readGrid(indir+'/'+id+'/'+id+'_CF_'+str(t)+'.data.gz') diffield = getRelativeDirField(sigma,r) #~ diffield = np.loadtxt(indir+'/'+id+'/'+id+'_reldir_r='+str(r)+'_'+str(t)+'.data.gz') clusters = getCellClusters(diffield,sigma,th,minlabsize,opendisk,mincellsize) cellids = np.unique(sigma) # remove ECM (sigma=0) from cellids for cid,cluster in clusters.iteritems(): csz = cluster.getClusterSize() # save data per cells for cellid in cluster.cells: p = np.where(sigma==cellid) cellTCdict.setdefault(cellid,ClusterCellTC(cellid)).addTimeStep(t,p,cid,csz) cellids = np.delete(cellids,np.where(cellids==cellid)) # save data for cells not in a cluster for cellid in cellids: if cellid == 0: continue p = np.where(sigma==cellid) cellTCdict.setdefault(cellid,ClusterCellTC(cellid)).addTimeStep(t,p,np.nan,np.nan) # write data to file per cell for cellid,celltc in cellTCdict.iteritems(): data = np.column_stack((celltc.time,calcMSDTransForCellTC(celltc.com), calcMSDRotForCellTC(celltc.laxis),celltc.clusterId,celltc.clusterSize)) f = open(outdir+'/'+id+'_cell_'+str(int(cellid))+'_cluster+MSD.data','w') f.write('#MCS\tMSDTrans\tMSDRot\tclusterId\tclusterSize\n') np.savetxt(f,data) f.close() print 'clusters calculated for '+str(len(trange))+' morphologies in '+str(time.time()-t0)+' seconds'
[docs]def createPBSScripts(runid,joblist,command,time,ncores=8,ppn=8,path='batchScripts/',xserver=False,python27=False): """ Create a set of PBS scripts to run a simulation on lisa. For each job in joblist a single line command is added to the script:: python command jobid > log/jobid.out 2> log/jobid.err & When the number of commands is equal to ppn, a new PBS script is created. :param runid: identifier for the scripts :param joblist: list of job identifiers :param command: command that runs the simulation :param time: requested walltime on the cluster :param ncores: types of nodes requested (either, 8, 12 or 16) :param ppn: number of processers per node that will be used :param path: location where pbs scripts are saved :param xserver: add lines to pbs script to mimic xserver :param python27: add lines to pbs script to use python 2.7 instead of the default python 2.6 .. seealso:: :func:`~SimUtils.CC3DPipeline.createPBS`, :func:`~SimUtils.CC3DPipeline.addCommandToPBS`, :func:`~SimUtils.CC3DPipeline.finishPBS` """ nscripts = 0 if ppn > ncores: ppn = ncores for i,jobid in enumerate(joblist): if (i%ppn == 0): if nscripts > 0: finishPBS(fn) fn = path+'/run_'+runid+'_part_'+str(nscripts) createPBS(fn,time,ncores,ppn,xserver,python27) nscripts += 1 addCommandToPBS(fn,command+' '+jobid,'log/'+jobid) finishPBS(fn)
def createPBS(filename,time,ncores=None,ppn=None,xserver=False,python27=False): """ Create a set of PBS script, set job settings and add generic commands:: #PBS -S /bin/bash #PBS -lnodes=1:ncores:ppn=ppn #PBS -lwalltime=time :param runid: identifier for the scripts :param joblist: list of job identifiers :param command: command that runs the simulation :param time: requested walltime on the cluster :param ncores: types of nodes requested (either, 8, 12 or 16) :param ppn: number of processers per node that will be used :param path: location where pbs scripts are saved :param xserver: add lines to pbs script to mimic xserver :param python27: add lines to pbs script to use python 2.7 instead of the default python 2.6 """ f = open(filename,'w') f.write('#PBS -S /bin/bash\n') f.write('#PBS -lnodes=1') if (ncores == 8) or (ncores == 12) or (ncores == 16): f.write(':cores'+str(ncores)) if ppn: f.write(':ppn='+str(ppn)) f.write('\n') f.write('#PBS -lwalltime='+str(time)+'\n') if xserver: f.write('export DISPLAY=:1\n') f.write('( Xvfb $DISPLAY -auth /dev/null & )\n') if python27: f.write('module load python/2.7.2 \n') f.write('cd $HOME\n') f.close() def addCommandToPBS(filename,command,log): """ Add single line command to existing PBS script:: python command > log.out 2> log/j.err &. :param filename: path to existing pbs script :param command: complete command (including arguments) :param log: path to and name of log files """ f = open(filename,'a') f.write(command+' > '+log+'.out 2> '+log+'.err &\n') f.close() def finishPBS(filename): """ Finish pbs file """ f = open(filename,'a') f.write('wait\n') f.close() def packScripts(runid,simlist,scriptdir,target='batchScripts'): cwd = os.getcwd() os.chdir(scriptdir) command = 'tar -cf scripts_'+runid+'.tar' for sim in simlist: command += ' '+str(sim)+'.xml' #~ command += ' '+str(sim)+'-*.xml' os.system(command) os.chdir(cwd) if os.path.isfile(target+'/scripts_'+runid+'.tar'): os.remove(target+'/scripts_'+runid+'.tar') shutil.move(scriptdir+'/scripts_'+runid+'.tar',target+'/') def packPBS(runid,target='batchScripts'): cwd = os.getcwd() os.chdir(target) command = 'tar -cf jobs_'+runid+'.tar run_'+runid+'*' os.system(command) os.chdir(cwd)