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)