casa: parallel/pimager.py Source File

Go to the documentation of this file.
00001 from taskinit import *
00002 from parallel_go import *
00003 from cleanhelper import *
00004 from parallel.parallel_cont import imagecont
00005 from simple_cluster import simple_cluster
00006 from odict import *
00007 from casac import casac
00008 import numpy as np
00009 import random
00010 import string
00011 import time
00012 import os
00013 import shutil
00014 import pdb
00015 import copy
00016 import glob
00017 class pimager():
00018     def __init__(self, cluster=''):
00019         self.msinfo=odict()
00020         self.engineinfo=odict()
00021         self.spw=''
00022         self.field=''
00023         self.phasecenter=''
00024         self.ftmachine='ft'
00025         self.wprojplanes=32
00026         self.facets=1
00027         self.imsize=[512, 512]
00028         self.cell=['1arcsec', '1arcsec']
00029         self.weight='natural'
00030         self.robust=0.0
00031         self.stokes='I'
00032         self.visinmem=False
00033         self.numthreads=-1
00034         self.gain=0.1
00035         self.npixels=0
00036         self.uvtaper=False
00037         self.outertaper=[]
00038         self.timerange=''
00039         self.uvrange=''
00040         self.baselines=''
00041         self.scan=''
00042         self.observation=''
00043         self.pbcorr=False
00044         self.minpb=0.2
00045         self.cyclefactor=1.5
00046         self.c=cluster
00047         if self.c == '' :
00048             # Until we move to the simple cluster
00049             self.c = simple_cluster.getCluster()._cluster
00050         os.environ['IPYTHONDIR']='./i_serpiante'  
00051         shutil.rmtree(os.environ['IPYTHONDIR'], True)
00052     def __del__(self):
00053         #print 'ipythondir', os.environ['IPYTHONDIR']
00054         shutil.rmtree(os.environ['IPYTHONDIR'], True)
00055     @staticmethod
00056     def maxouterpsf(psfim='', image=''):
00057         ia.open(psfim)
00058         stat=ia.statistics()
00059         csys=ia.coordsys()
00060         ia.open(image)
00061         beampix=np.fabs(ia.restoringbeam()['major']['value']/qa.convert(csys.increment(type='direction', format='q')['quantity']['*1'], ia.restoringbeam()['major']['unit'])['value'])
00062         beampix=int(np.ceil(beampix))
00063         blc=copy.deepcopy(stat['maxpos'])
00064         trc=copy.deepcopy(stat['maxpos'])
00065         blc[0:2]=blc[0:2]-2*beampix
00066         trc[0:2]=trc[0:2]+2*beampix
00067         blcpeak=copy.deepcopy(blc)
00068         trcpeak=copy.deepcopy(trc)
00069         blcpeak[2:4]=0
00070         trcpeak[2:4]=0
00071         #print 'blc trc peak', blcpeak, trcpeak
00072         blc[0:2]=0
00073         trc[0:2]=ia.shape()[0]-1;
00074         maxplane=rg.box(blc=blc.tolist(),trc=trc.tolist())
00075         ia.open(psfim)
00076         ib=ia.subimage(region=maxplane)
00077         csys2=ib.coordsys()
00078         boxpeak=rg.wbox(blc=["%dpix"%x for x in blcpeak], trc=["%dpix"%x for x in trcpeak], csys=csys2.torecord())
00079         outerreg=rg.complement(boxpeak)
00080         statout=ib.statistics(region=outerreg)
00081         ia.done()
00082         ib.done()
00083         return np.max(statout['max'], np.fabs(statout['min']))
00084     @staticmethod
00085     def averimages(outimage='outimage', inimages=[]):
00086         if((type(inimages)==list) and (len (inimages)==0)):
00087             return False
00088         if(os.path.exists(outimage)):
00089             shutil.rmtree(outimage)
00090         if(type(inimages)==list):
00091             shutil.copytree(inimages[0], outimage)
00092             ia.open(outimage)        
00093             for k in range(1, len(inimages)) :
00094                 ia.calc('"'+outimage+'" +  "'+inimages[k]+'"')
00095             ia.calc('"'+outimage+'"'+'/'+str(len(inimages)))
00096             ia.done()
00097         elif(type(inimages)==str):
00098             shutil.copytree(inimages, outimage)
00099         else:
00100             return False
00101         return True  
00102     @staticmethod
00103     def weightedaverimages(outimage='outimage', inimages=[], wgtimages=[]):
00104         if((type(inimages)==list) and (len (inimages)==0)):
00105             return False
00106         if(len(wgtimages) != len (inimages)):
00107             raise Exception, 'number of images and weight images are different'
00108         if(os.path.exists(outimage)):
00109             shutil.rmtree(outimage)
00110         shutil.rmtree( '__sumweight_image', ignore_errors=True)
00111         if(type(inimages)==list):
00112             shutil.copytree(inimages[0], outimage)
00113             if(len(inimages)==1):
00114                 return True
00115             shutil.copytree(wgtimages[0], '__sumweight_image')
00116             ia.open(outimage)  
00117             ia.calc('"'+inimages[0]+'"*"'+wgtimages[0]+'"')
00118             ib,=gentools(['ia'])
00119             ib.open( '__sumweight_image')
00120             for k in range(1, len(inimages)) :
00121                 ia.calc('"'+outimage+'" +  "'+inimages[k]+'"*"'+wgtimages[k]+'"')
00122                 ib.calc('"'+'__sumweight_image' + '" +  "'+wgtimages[k]+'"')
00123             ib.done()
00124             ia.calc('"'+outimage+'"'+'/"__sumweight_image"')
00125             ia.done()
00126             shutil.rmtree( '__sumweight_image')
00127         elif(type(inimages)==str):
00128             shutil.copytree(inimages, outimage)
00129         else:
00130             return False
00131         return True
00132   
00133     @staticmethod
00134     def findchansel(msname='', spw='*', field='*', numpartition=1, freqrange=[0, 1e12]):
00135         im,ms=gentools(['im', 'ms'])
00136         if(field==''):
00137             field='*'
00138         fieldid=0
00139         fieldids=[]
00140         if(type(field)==str):
00141             fieldids=ms.msseltoindex(vis=msname, field=field)['field']
00142         elif(type(field)==int):
00143             fieldids=[field]
00144         elif(type(field)==list):
00145             fieldids=field
00146         a={}
00147         a['freqstart']=1.7977e308
00148         fid=0
00149         while((a['freqstart']> 1.0e100) and (fid < len(fieldids))):
00150             a=im.advisechansel(msname=msname, getfreqrange=True, fieldid=fieldids[fid], spwselection=spw)
00151             fid+=1
00152                          
00153         freqrange[0]=a['freqstart']
00154         freqrange[1]=a['freqend']
00155         thesel=ms.msseltoindex(msname, spw=spw)
00156         spwids=thesel['spw']
00157         nchans=np.sum(thesel['channel'][:,2]-thesel['channel'][:,1]+1)
00158         
00159         if(nchans < 1):
00160             return []
00161         tb.open(msname)
00162         spectable=string.split(tb.getkeyword('SPECTRAL_WINDOW'))
00163         if(len(spectable) ==2):
00164             spectable=spectable[1]
00165         else:
00166             spectable=msname+"/SPECTRAL_WINDOW"
00167         tb.open(spectable)
00168         elfreq=[]
00169         nspw=tb.nrows()
00170         for k in range(nspw):
00171             elfreq.append(tb.getcol('CHAN_FREQ', k, 1).flatten())
00172         tb.done()
00173         beginfreq=1e12
00174         endfreq=0.0
00175         ####number of channel for each partition
00176         nch=np.zeros(numpartition)
00177         nch[:]=int(nchans/numpartition)
00178         nch[0:(nchans%numpartition)] +=1
00179         sel=['-1']*numpartition
00180         row=0
00181         nextstart=thesel['channel'][row,1]-1
00182         nextspw=thesel['channel'][row,0]
00183         nextend=thesel['channel'][row,2]
00184         for k in range(numpartition):
00185             startch=nextstart+1
00186             spwsel=nextspw
00187             endch=nextstart
00188             for jj in range(int(nch[k])):
00189                 nextstart=endch+1
00190                 chanval=endch if (endch > -1) else 0
00191                 beginfreq=beginfreq if(beginfreq < elfreq[spwsel][chanval]) else elfreq[spwsel][chanval]
00192                 endfreq=endfreq if(endfreq > elfreq[spwsel][chanval]) else elfreq[spwsel][chanval]
00193                 endch +=1
00194                 if(endch>nextend):
00195                     row+=1
00196                     nextend=thesel['channel'][row,2]
00197                     if(startch <= (endch-1)):
00198                         sel[k]=(sel[k]+','+str(spwsel)+':'+str(startch)+'~'+str(endch-1)) if(sel[k] != '-1') else (str(spwsel)+':'+str(startch)+'~'+str(endch-1))
00199                     startch=thesel['channel'][row,1] 
00200                     endch=startch
00201                     nextstart=endch+1
00202                     spwsel=thesel['channel'][row,0]
00203                     nextspw=spwsel
00204             if(startch <= (endch)):
00205                 sel[k]=(sel[k]+','+str(spwsel)+':'+str(startch)+'~'+str(endch)) if(sel[k] != '-1') else (str(spwsel)+':'+str(startch)+'~'+str(endch))
00206             if(startch==endch):
00207                 nextstart=nextstart-1
00208         #freqrange[0]=beginfreq
00209         #freqrange[1]=endfreq
00210         #print  'FREQRANGE', freqrange
00211 
00212         return sel
00213                 
00214 
00215     @staticmethod
00216     def findchanselcont(msname='', spwids=[], numpartition=1, beginfreq=0.0, endfreq=1e12, freqrange=[0, 1e12]):
00217         numproc=numpartition
00218         spwsel=[]
00219         startsel=[]
00220         nchansel=[]
00221         for k in range(numproc):
00222             spwsel.append([])
00223             startsel.append([])
00224             nchansel.append({})
00225 
00226         tb.open(msname)
00227         spectable=string.split(tb.getkeyword('SPECTRAL_WINDOW'))
00228         if(len(spectable) ==2):
00229             spectable=spectable[1]
00230         else:
00231             spectable=msname+"/SPECTRAL_WINDOW"
00232         tb.open(spectable)
00233         channum=tb.getcol('NUM_CHAN')
00234         nspw=tb.nrows()
00235         if(len(spwids)==0):
00236             spwids=range(nspw)
00237         freqs=[]
00238         for k in range(nspw):
00239             freqs.append(tb.getcol('CHAN_FREQ', k, 1).flatten())
00240 
00241         allfreq=freqs[spwids[0]]
00242         flatspw=np.ndarray((channum[spwids[0]]), int)
00243         flatspw[:]=spwids[0]
00244         flatchan=np.array(range(channum[spwids[0]]))
00245         for k in range(1, len(spwids)) :
00246             allfreq=np.append(allfreq, freqs[spwids[k]])
00247             tpflatspw=np.ndarray((channum[spwids[k]]), int)
00248             tpflatspw[:]=spwids[k]
00249             flatspw=np.append(flatspw, tpflatspw)
00250             flatchan=np.append(flatchan, np.array(range(channum[spwids[k]])))
00251          ##number of channels in the ms
00252         #pdb.set_trace()
00253         sortind=np.argsort(allfreq)
00254         allfreq=np.sort(allfreq)
00255         numchanperms=len(allfreq)
00256         #print 'number of channels', numchanperms
00257         minfreq=np.min(allfreq)
00258         if(minfreq > endfreq):
00259             return -1, -1, -1
00260         minfreq=minfreq if minfreq > beginfreq else beginfreq
00261         maxfreq=np.max(allfreq)
00262         if(maxfreq < beginfreq):
00263             return -1, -1, -1
00264         maxfreq=maxfreq if maxfreq < endfreq else endfreq
00265         ###modify the beginfreq and endfreq
00266         freqrange[0]=minfreq
00267         freqrange[1]=maxfreq
00268         startallchan=0
00269         while (minfreq > allfreq[startallchan]):
00270             startallchan=startallchan+1
00271         if(startallchan > (allfreq.size -1) ):
00272             #return -1
00273             return -1, -1, -1
00274         endallchan=allfreq.size -1
00275         while((maxfreq < allfreq[endallchan]) and (endallchan > 0)):
00276             endallchan=endallchan-1
00277         if(endallchan < (startallchan+1)):
00278             ##fail
00279             return -1, -1, -1
00280         totchan=endallchan-startallchan+1
00281         bandperproc=(allfreq[endallchan]-allfreq[startallchan])/numproc
00282         chanperproc=totchan/numproc
00283         extrachan=0
00284         if((totchan%numproc) > 0):
00285             chanperproc=chanperproc+1
00286             extrachan=1
00287         actualchan=0
00288 
00289         spwcounter=0
00290         for k in range(numproc):
00291             spwsel[k].append(flatspw[sortind[actualchan]])
00292             startsel[k].append(flatchan[sortind[actualchan]])
00293             for j in range(chanperproc):
00294                 if(((actualchan%chanperproc) > 0) and (flatspw[sortind[actualchan]] !=flatspw[sortind[actualchan-1]])):
00295                     nchansel[k].update({flatspw[sortind[actualchan-1]]:(flatchan[sortind[actualchan-1]]-startsel[k][spwcounter]+1)})
00296                     if(not np.any(np.array(spwsel[k])==flatspw[sortind[actualchan]])):
00297                         spwcounter=spwcounter+1
00298                         spwsel[k].append(flatspw[sortind[actualchan]])
00299                         startsel[k].append(flatchan[sortind[actualchan]])
00300                     else:
00301                         ###we are back on a spw that is already selected
00302                         nchansel[k].pop(flatspw[sortind[actualchan]])
00303                         spwcounter=0
00304                         while (spwsel[k][spwcounter] != flatspw[sortind[actualchan]]):
00305                             spwcounter+=1
00306                 actualchan=actualchan+1
00307                 if(actualchan==totchan):
00308                     break
00309             #pdb.set_trace()
00310             if(len(nchansel[k]) != len(spwsel[k])):
00311                 nchansel[k].update({spwsel[k][spwcounter]:(flatchan[sortind[actualchan-1]]-startsel[k][spwcounter]+1)})
00312             if(actualchan==totchan):
00313                 break
00314             spwcounter=0
00315         retchan=[]
00316         for k in range(numproc):
00317             retchan.append([])
00318             for j in range(len(nchansel[k])):
00319                 retchan[k].append(nchansel[k][spwsel[k][j]])
00320 
00321         return spwsel, startsel, retchan
00322 
00323     @staticmethod
00324     def findchanselold(msname='', spwids=[], numpartition=1, beginfreq=0.0, endfreq=1e12, continuum=True):
00325         numproc=numpartition
00326         spwsel=[]
00327         startsel=[]
00328         nchansel=[]
00329         for k in range(numproc):
00330             spwsel.append([])
00331             startsel.append([])
00332             nchansel.append([])
00333         tb.open(msname)
00334         spectable=string.split(tb.getkeyword('SPECTRAL_WINDOW'))
00335         if(len(spectable) ==2):
00336             spectable=spectable[1]
00337         else:
00338             spectable=msname+"/SPECTRAL_WINDOW"
00339         tb.open(spectable)
00340         channum=tb.getcol('NUM_CHAN')
00341         nspw=tb.nrows()
00342         if(len(spwids)==0):
00343             spwids=range(nspw)
00344         freqs=[]
00345         for k in range(nspw):
00346             freqs.append(tb.getcol('CHAN_FREQ', k, 1).flatten())
00347 
00348         allfreq=freqs[spwids[0]]
00349         for k in range(1, len(spwids)) :
00350             allfreq=np.append(allfreq, freqs[spwids[k]])
00351          ##number of channels in the ms
00352         allfreq=np.sort(allfreq)
00353         numchanperms=len(allfreq)
00354         #print 'number of channels', numchanperms
00355         minfreq=np.min(allfreq)
00356         if(minfreq > endfreq):
00357             return -1, -1, -1
00358         minfreq=minfreq if minfreq > beginfreq else beginfreq
00359         maxfreq=np.max(allfreq)
00360         if(maxfreq < beginfreq):
00361             return -1, -1, -1
00362         maxfreq=maxfreq if maxfreq < endfreq else endfreq
00363         startallchan=0
00364         while (minfreq > allfreq[startallchan]):
00365             startallchan=startallchan+1
00366         if(startallchan > (allfreq.size -1) ):
00367             #return -1
00368             return -1, -1, -1
00369         endallchan=allfreq.size -1
00370         while((maxfreq < allfreq[endallchan]) and (endallchan > 0)):
00371             endallchan=endallchan-1
00372         if(endallchan < (startallchan+1)):
00373             ##fail
00374             return -1, -1, -1
00375         totchan=endallchan-startallchan+1
00376         bandperproc=(allfreq[endallchan]-allfreq[startallchan])/numproc
00377         chanperproc=totchan/numproc
00378         extrachan=0
00379         if((totchan%numproc) > 0):
00380             chanperproc=chanperproc+1
00381             extrachan=1
00382         startthissel=startallchan
00383         #endthissel=chanperproc+startthissel-1
00384         #if(endthissel > (len(allfreq))-1):
00385         #    enthissel=(len(allfreq))-1
00386         lowfreq=allfreq[startthissel]
00387         upfreq=lowfreq+bandperproc
00388         #pdb.set_trace()
00389         for k in range(numproc):
00390             donelow=False;
00391             doneup=False
00392             #while ((not donelow) or (not doneup)):
00393             for uu in range(1):
00394                 #####################
00395                 for j in range(nspw):
00396                     #print lowfreq, upfreq, np.min(freqs[j]), np.max(freqs[j])
00397                     ind=np.argsort(freqs[j])
00398                     if((lowfreq > freqs[j][ind[channum[j]-1]]) or 
00399                             (upfreq < freqs[j][ind[0]])):
00400                         #this spectral window is not in
00401                         continue
00402                     elif((lowfreq >= freqs[j][ind[0]]) and 
00403                        (upfreq <= freqs[j][ind[channum[j]-1]])):
00404                         ###whole selected band is in this spw
00405                         donelow=True
00406                         doneup=True
00407                         spwsel[k].append(j)
00408                         cc=0
00409                         if(freqs[j][channum[j]-1] > freqs[j][0]):
00410                             while(freqs[j][ind[cc]] < lowfreq):
00411                                 cc=cc+1
00412                         else:
00413                             while(freqs[j][ind[cc]] > upfreq):
00414                                 cc=cc+1
00415                         cc= (cc-extrachan) if (cc-extrachan) > 0 else 0 
00416                         startsel[k].append(ind[cc])
00417                         nchansel[k].append(chanperproc)
00418                     elif((lowfreq <= freqs[j][ind[0]]) and 
00419                          (upfreq >= freqs[j][ind[channum[j]-1]])):
00420                         ###this spw is within the whole range
00421                         spwsel[k].append(j)
00422                         startsel[k].append(0)
00423                         nchansel[k].append(channum[j])
00424                     elif((lowfreq >= freqs[j][ind[0]]) and 
00425                          (lowfreq < freqs[j][ind[channum[j]-1]]) and
00426                          (upfreq > freqs[j][ind[channum[j]-1]])):
00427                         ###lowbound in this spw but not upbound
00428                         donelow=True
00429                         spwsel[k].append(j)
00430                         matchlow=0
00431                         while(lowfreq > freqs[j][ind[matchlow]]):
00432                             matchlow=matchlow+1
00433                         matchlow=matchlow-1 if matchlow > 0 else 0
00434                         if(freqs[j][channum[j]-1] > freqs[j][0]):
00435                             #positive inc
00436                             startsel[k].append(ind[matchlow])
00437                             nchansel[k].append(channum[j]-ind[matchlow])
00438                         else:
00439                             #neg inc
00440                             startsel[k].append(0)
00441                             nchansel[k].append(ind[matchlow])
00442                     elif((upfreq > freqs[j][ind[0]]) and 
00443                          (upfreq <= freqs[j][ind[channum[j]-1]]) and
00444                          (lowfreq < freqs[j][ind[0]])):
00445                         ###upbound in this spw but not lowbound
00446                         doneup=True
00447                         spwsel[k].append(j)    
00448                         matchup=0
00449                         while(upfreq > freqs[j][ind[matchup]]):
00450                             matchup=matchup+1
00451                         matchup=matchup-1 if matchup > 0 else 0
00452                         if(freqs[j][channum[j]-1] > freqs[j][0]):
00453                             #positive inc
00454                             startsel[k].append(0)
00455                             nchansel[k].append(ind[matchup]+1)
00456                         else:
00457                             #neg inc
00458                             startsel[k].append(ind[matchup])
00459                             nchansel[k].append(channum[j]-ind[matchup]+1)
00460                     else:
00461                         print 'Huh ?'
00462             #startthissel=endthissel
00463             #endthissel=chanperproc+startthissel-1
00464             #if(endthissel > (len(allfreq))-1):
00465             #    enthissel=(len(allfreq))-1
00466             startthissel=startthissel+chanperproc-extrachan
00467             lowfreq=upfreq
00468             #lowfreq=upfreq if (upfreq < allfreq[startthissel]) else allfreq[startthissel]
00469             upfreq=lowfreq+bandperproc
00470         if(continuum):
00471             #no need to send same data to different processors
00472             for k in range (numproc-1):
00473                 for j in range(k+1, numproc):
00474                     for kk in range(len(spwsel[k])):
00475                         spwtopop=[]
00476                         for jj in range(len(spwsel[j])):
00477                             if(spwsel[k][kk] == spwsel[j][jj]):
00478                                 while(startsel[j][jj] < (startsel[k][kk]+nchansel[k][kk])):
00479                                     startsel[j][jj]=startsel[j][jj]+1
00480                                     nchansel[j][jj]=nchansel[j][jj]-1
00481                                 if(nchansel[j][jj] <=0):
00482                                     spwtopop.append(jj)
00483                         for jj in range(len(spwtopop)):
00484                             startsel[j].pop(spwtopop[jj])
00485                             spwsel[j].pop(spwtopop[jj])
00486                             nchansel[j].pop(spwtopop[jj])
00487 
00488         return spwsel, startsel, nchansel
00489 
00490     @staticmethod
00491     def copyimage(inimage='', outimage='', init=False, initval=0.0):
00492         ia.fromimage(outfile=outimage, infile=inimage, overwrite=True)
00493         ia.open(outimage)
00494         if(init):
00495             ia.set(initval)
00496         else:
00497             ####ib=iatool.create()
00498             ####ib.open(inimage)
00499             ####arr=ib.getchunk()
00500             ####ib.done()
00501             ia.insert(inimage, locate=[0,0,0,0])
00502         ia.done()
00503     @staticmethod
00504     def regridimage(outimage='', inimage='', templateimage=''):
00505         ia.open(templateimage)
00506         csys=ia.coordsys()
00507         shp=ia.shape()
00508         ia.open(inimage)
00509         ia.regrid(outfile=outimage, shape=shp, csys=csys.torecord(), axes=[0,1], overwrite=True)
00510         ia.done()
00511     def setupcommonparams(self, spw='*', field='*', phasecenter='', 
00512                           stokes='I', ftmachine='ft', wprojplanes=64, facets=1, 
00513                           imsize=[512, 512], pixsize=['1arcsec', '1arcsec'], weight='natural',
00514                           robust=0.0, npixels=0, gain=0.1,  uvtaper=False, outertaper=[], 
00515                           timerange='', uvrange='', baselines='', scan='', observation='', 
00516                           visinmem=False, pbcorr=False, minpb=0.2, numthreads=1, cyclefactor=1.5):
00517         self.spw=spw
00518         self.field=field
00519         self.phasecenter=phasecenter
00520         self.stokes=stokes
00521         self.ftmachine=ftmachine
00522         self.wprojplanes=wprojplanes
00523         self.facets=facets
00524         self.imsize=imsize
00525         self.cell=pixsize
00526         self.weight=weight
00527         self.robust=robust
00528         self.npixels=npixels
00529         self.gain=gain
00530         self.uvtaper=uvtaper
00531         self.outertaper=outertaper
00532         self.timrange=timerange
00533         self.uvrange=uvrange
00534         self.baselines=baselines
00535         self.scan=scan
00536         self.observation=observation
00537         self.visinmem=visinmem
00538         self.pbcorr=pbcorr
00539         self.minpb=minpb
00540         self.numthreads=numthreads
00541         self.cyclefactor=cyclefactor
00542     
00543     def pcontmt(self, msname=None, imagename=None, imsize=[1000, 1000], 
00544               pixsize=['1arcsec', '1arcsec'], phasecenter='', 
00545               field='', spw='*', stokes='I', ftmachine='ft', wprojplanes=128, facets=1, 
00546               majorcycles=-1, cyclefactor=1.5, niter=1000, npercycle=100, gain=0.1, threshold='0.0mJy', alg='clark', scales=[0], weight='natural', robust=0.0, npixels=0,  uvtaper=False, outertaper=[], timerange='', uvrange='', baselines='', scan='', observation='', pbcorr=False, minpb=0.2, 
00547               contclean=False, visinmem=False, interactive=False, maskimage='lala.mask',
00548               numthreads=1, savemodel=False, nterms=2):
00549         if(nterms==1):
00550             self.pcont(msname=msname, imagename=imagename, imsize=imsize, 
00551               pixsize=pixsize, phasecenter=phasecenter, field=field, spw=spw, stokes=stokes, ftmachine=ftmachine, wprojplanes=wprojplanes, facets=facets,  
00552               majorcycles=majorcycles, cyclefactor=cyclefactor, niter=niter, npercycle=npercycle, gain=gain, 
00553               threshold=threshold, alg=alg, scales=scales, weight=weight, robust=robust, 
00554               npixels=npixels,  uvtaper=uvtaper, outertaper=outertaper, 
00555               timerange=timerange, uvrange=uvrange, baselines=baselines, scan=scan, 
00556               observation=observation, pbcorr=pbcorr, minpb=minpb, contclean=contclean, 
00557               visinmem=visinmem, interactive=interactive, maskimage=maskimage,
00558               numthreads=numthreads, savemodel=savemodel)
00559         else:
00560             self.setupcommonparams(spw=spw, field=field, phasecenter=phasecenter, 
00561                                    stokes=stokes, ftmachine=ftmachine, wprojplanes=wprojplanes, 
00562                                    facets=facets, imsize=imsize, pixsize=pixsize, weight=weight, 
00563                                    robust=robust, npixels=npixels, gain=gain, uvtaper=uvtaper,
00564                                    outertaper=outertaper, timerange=timerange, uvrange=uvrange, 
00565                                    baselines=baselines, scan=scan, observation=observation, 
00566                                    visinmem=visinmem, pbcorr=pbcorr, minpb=minpb, numthreads=numthreads, 
00567                                    cyclefactor=cyclefactor)
00568             dc=casac.deconvolver()
00569             ia=casac.image()
00570             niterpercycle=niter/majorcycles if(majorcycles >0) else niter
00571             if(niterpercycle == 0):
00572                 niterpercycle=niter
00573                 majorcycles=1
00574             self.setupcluster()
00575             spwids=ms.msseltoindex(vis=msname, spw=self.spw)['spw']
00576             timesplit=0
00577             timeimage=0
00578             elimageroot=imagename
00579             elmask=maskimage
00580             owd=os.getcwd()
00581             fullpath=lambda a: owd+'/'+a if ((len(a) !=0) and a[0] != '/') else a
00582             imagename=fullpath(elimageroot)
00583             maskimage=fullpath(elmask) 
00584             models=[]
00585             psfs=[]
00586             residuals=[]
00587             restoreds=[]
00588             sumwts=[]
00589             ntaylor=nterms
00590             npsftaylor = 2 * nterms - 1
00591             for tt in range(0, nterms):
00592                 models.append(imagename+'.model.tt'+str(tt))
00593                 residuals.append(imagename+'.residual.tt'+str(tt));
00594                 restoreds.append(imagename+'.image.tt'+str(tt));
00595                 sumwts.append(imagename+'.sumwt.tt'+str(tt));
00596             #tempmodel=owd+'/tempmodel'
00597             #shutil.rmtree(tempmodel, True)
00598             for tt in range(0,npsftaylor):
00599               psfs.append(imagename+'.psf.tt'+str(tt));
00600             if(not contclean):
00601                 toberemoved=glob.glob(imagename+'*')
00602                 ##keep masks
00603                 tbr=copy.deepcopy(toberemoved)
00604                 for k in range(len(toberemoved)):
00605                     if(string.find(toberemoved[k], '.mask') > 0):
00606                            tbr.remove(toberemoved[k])
00607                 toberemoved=tbr
00608                 for rem in toberemoved:
00609                     print "Removing ", rem
00610                     shutil.rmtree(rem, True) 
00611             ###num of cpu 
00612             numcpu=self.numcpu
00613             time1=time.time()
00614             ###spw and channel selection
00615             #spwsel,startsel,nchansel=self.findchanselcont(msname, spwids, numcpu)
00616 
00617             #print 'SPWSEL ', spwsel, startsel, nchansel
00618             freqrange=[0.0, 0.0]
00619             spwsel=self.findchansel(msname, spw, field, numcpu, freqrange=freqrange)
00620             minfreq=freqrange[0]
00621             maxfreq=freqrange[1]
00622 
00623             out=range(numcpu) 
00624         
00625             c=self.c
00626             ####
00627             #the default working directory is somewhere 
00628             owd=os.getcwd()
00629             self.c.pgc('import os')
00630             self.c.pgc('os.chdir("'+owd+'")')
00631             ##################### 
00632             #c.pgc('casalog.filter()')
00633             c.pgc('from  parallel.parallel_cont import *')
00634 
00635             band=maxfreq-minfreq
00636             ###need to get the data selection for each process here
00637             ## this algorithm is a poor first try
00638             if(minfreq <0 ):
00639                 minfreq=0.0
00640             freq='"%s"'%(str((minfreq+band/2.0))+'Hz')
00641             band='"%s"'%(str((band*1.1))+'Hz')
00642             ###define image names
00643             imlist=[]
00644             char_set = string.ascii_letters
00645             cfcachelist=[]
00646             for k in range(numcpu):
00647                 substr=self.workingdirs[k]+'/Temp_'+str(k)+'_'+string.join(random.sample(char_set,8), sep='')
00648                 while (os.path.exists(substr+'.model.tt0')):
00649                     substr=self.workingdirs[k]+'/Temp_'+str(k)+'_'+string.join(random.sample(char_set,8), sep='')
00650                    ###############
00651                 imlist.append(substr)
00652         ##continue clean or not
00653             if(contclean and os.path.exists(models[0])):
00654                 for k in range(numcpu):
00655                     for tt in range(0, nterms):
00656                         self.copyimage(inimage=models[tt], outimage=imlist[k]+'.model.tt'+str(tt), init=False)
00657             else:
00658                 for k in range(len(imlist)):
00659                     for tt in range(nterms):
00660                         shutil.rmtree(imlist[k]+'.model.tt'+str(tt), True)
00661             intmask=0
00662             donewgt=[False]*numcpu
00663          ### do one major cycle at the end
00664             donemaj=False
00665             maj=0;
00666             maxresid=-1.0
00667             while(not donemaj):
00668             
00669                 casalog.post('Starting Gridding for major cycle '+str(maj)) 
00670                 for k in range(numcpu):
00671                     imnam='"%s"'%(imlist[k])
00672                     runcomm='a.imagecont(msname='+'"'+msname+'", field="'+str(field)+'", spw="'+str(spwsel[k])+'", freq='+freq+', band='+band+', imname='+imnam+', nterms='+str(nterms)+', scales='+str(scales)+')'
00673                     print 'command is ', runcomm
00674                     out[k]=c.odo(runcomm,k)
00675                 over=False
00676                 printkounter=0
00677                 while(not over):
00678                     time.sleep(5)
00679                     overone=True
00680                     for k in range(numcpu):
00681                         cj=c.check_job(out[k],False)
00682                         overone=(overone and cj)
00683                     #print 'maj, cj, dwght', maj, cj, donewgt[k]                    
00684                     over=overone
00685                 self.combineimages(rootnames=imlist, nterms=nterms, outputrootname=imagename)
00686                 if((maskimage == '') or (maskimage==[])):
00687                     maskimage=imagename+'.mask'
00688                     ia.removefile(maskimage)
00689                 if (interactive and (intmask==0)):
00690                     if(maj==0):
00691                         ia.removefile(maskimage)
00692                     retdraw=im.drawmask(residuals[0],maskimage, niter=niterpercycle, npercycle=npercycle, threshold=threshold);
00693                     intmask=retdraw['stat']
00694                     ####The user may want to change these in the gui
00695                     threshold=retdraw['threshold']
00696                     niterpercycle=retdraw['niter']
00697                     npercycle=retdraw['npercycle']
00698                     print 'intmask', intmask
00699                     if(intmask==1):
00700                         interactive=False
00701                         im.done()
00702                     if(intmask==2):
00703                         interactive=False
00704                         break;
00705                 else:
00706                     if (maj==0):
00707                         if(os.path.exists(maskimage)):
00708                             self.regridimage(outimage='__lala.mask', inimage=maskimage, templateimage=residuals[0]);
00709                             shutil.rmtree(maskimage, True)
00710                             shutil.move('__lala.mask', maskimage)
00711                         else:
00712                             print 'DOING a full image mask'
00713                             self.copyimage(inimage=residuals[0], outimage=maskimage, init=True, initval=1.0)
00714             #########Things to be done after first major cycle only
00715            #########
00716                 if(maj==0):
00717                     dc.mtopen(ntaylor=nterms, scalevector=scales, psfs=psfs)
00718                     if(not contclean or (not os.path.exists(model))):
00719                         for tt in range(nterms):
00720                             self.copyimage(inimage=residuals[tt], outimage=models[tt], 
00721                                            init=True, initval=0.0) 
00722                 ##for mosaic stuff needs to be done here
00723                 
00724                     if((self.weight=='uniform') or (self.weight=='briggs')):
00725                         c.pgc('wtgrid=a.getweightgrid(msname="'+msname+'")')
00726                         sumweight=c.pull('wtgrid', 0)[0]
00727                         for jj in range(1,numcpu):
00728                             sumweight += c.pull('wtgrid', jj)[jj]
00729                             print 'Max min var of weight dens', np.max(sumweight), np.min(sumweight), np.var(sumweight)
00730                         c.push(wtgrid=sumweight)
00731                         c.pgc('a.setweightgrid(msname="'+msname+'", weight=wtgrid)')
00732                 newthresh=threshold
00733                 if(majorcycles <= 1):
00734                     ia.open(residuals[0])
00735                     residstat=ia.statistics()
00736                     ia.done()
00737                     maxresid=np.max([residstat['max'], np.fabs(residstat['min'])])
00738                     psfoutermax=self.maxouterpsf(psfim=psfs[0], image=restoreds[0])
00739                     newthresh=psfoutermax*cyclefactor*maxresid
00740                     oldthresh=qa.convert(qa.quantity(threshold, "Jy"), "Jy")['value']
00741                     if(newthresh < oldthresh):
00742                         newthresh=oldthresh
00743                     newthresh=qa.tos(qa.quantity(newthresh, "Jy"))
00744             ####no need to do this in last major cycle
00745                 needclean=((maj < majorcycles if (majorcycles >1) else True)  and 
00746                            ((maxresid > qa.convert(qa.quantity(threshold),'Jy')['value']) if(majorcycles <2) else True) 
00747                            and (niterpercycle > 1))
00748                 donemaj = not needclean
00749 
00750                 if(needclean):
00751                     elniter=npercycle if(interactive) else niterpercycle
00752                     retval = dc.mtclean(residuals=residuals, models=models, 
00753                                         niter=elniter, gain=gain, threshold=newthresh, 
00754                                         displayprogress=True, mask=maskimage)
00755                     print 'mtclean return', retval
00756                     if(majorcycles <=1):
00757                         niterpercycle=niterpercycle-retval['iterations']
00758                     maxresid=retval['maxresidual']
00759                     ia.open(models[0])
00760                     print 'min max of model', ia.statistics()['min'], ia.statistics()['max']
00761                     ia.done()
00762                     for k in range(len(imlist)):
00763                         for tt in range(nterms):
00764                             ia.open(imlist[k]+'.model.tt'+str(tt))
00765                             ia.insert(infile=models[tt], locate=[0,0,0,0])
00766                             ia.done()
00767                     maj +=1
00768             ia.open(imlist[numcpu-1]+'.image.tt0')
00769             beam=ia.restoringbeam()
00770             ia.done()
00771             dc.mtrestore(models=models, residuals=residuals, images=restoreds,
00772                          bmaj=beam['major'], bmin=beam['minor'], bpa=beam['positionangle'])
00773             alphaname = imagename+'.alpha'
00774             betaname = imagename+'.beta'
00775             dc.mtcalcpowerlaw(images=restoreds, residuals=residuals,
00776                               alphaname=alphaname, 
00777                               betaname=betaname,
00778                               threshold='0.001Jy', calcerror=True)
00779             dc.done
00780             c.pgc('del a')
00781 
00782 
00783     def combineimages(self, rootnames=[], nterms=2, outputrootname=''):
00784         combmodels=[]
00785         combpsfs=[]
00786         combresiduals=[]
00787         combwts=[]
00788         combimages=[]
00789         ncpu = len(rootnames)
00790         for tt in range(0,nterms):  
00791             combmodels.append(outputrootname+'.model.tt'+str(tt))
00792             combresiduals.append(outputrootname+'.residual.tt'+str(tt))
00793             combwts.append(outputrootname+'.sumwt.tt'+str(tt))
00794             combimages.append(outputrootname+'.image.tt'+str(tt))
00795             if not os.path.exists( combmodels[tt] ):
00796                 self.copyimage(inimage=rootnames[0]+'.model.tt'+ str(tt), outimage=combmodels[tt], init=True)
00797             if not os.path.exists( combresiduals[tt] ):
00798                 self.copyimage(inimage=rootnames[0]+'.residual.tt'+ str(tt), outimage=combresiduals[tt], init=True)
00799             if not os.path.exists( combwts[tt] ):
00800                 self.copyimage(inimage=rootnames[0]+'.sumwt.tt'+ str(tt), outimage=combwts[tt], init=True)
00801             if not os.path.exists( combimages[tt] ):
00802                 self.copyimage(inimage=rootnames[0]+'.image.tt'+ str(tt), outimage=combimages[tt], init=True)
00803             for chunk in range(0,ncpu):
00804                 chunkresidual=rootnames[chunk]+'.residual.tt'+str(tt)
00805                 ia.open(combresiduals[tt])
00806                 ia.calc( '"'+combresiduals[tt] + '"+"' + chunkresidual + '"*"' + rootnames[chunk]+'.sumwt.tt0"' )
00807                 ###note the above is multiplying by weight of tt0 
00808                 ia.close()
00809                 ia.open(combwts[tt])
00810                 ia.calc( '"'+combwts[tt] + '"+"' + rootnames[chunk]+'.sumwt.tt'+str(tt)+'"')
00811                 ia.close()
00812             ## Normalize residuals by number of chunks.
00813             ia.open( combresiduals[tt] )
00814             ia.calc( '"'+combresiduals[tt] + '"/"' + combwts[0]+'"' )
00815             ia.done()
00816         for tt in range(0,2*nterms-1):  
00817             combpsfs.append(outputrootname+'.psf.tt'+str(tt));
00818             if not os.path.exists( combpsfs[tt] ):
00819                 self.copyimage(inimage=rootnames[0]+'.psf.tt'+ str(tt), outimage=combpsfs[tt], init=True)
00820             ## Add together all chan-chunk images.
00821             for chunk in range(0,ncpu):
00822                 chunkpsf=rootnames[chunk]+'.psf.tt'+str(tt)
00823                 ia.open(combpsfs[tt])
00824                 ia.calc( '"'+combpsfs[tt] + '"+"' + chunkpsf + '"*"' + rootnames[chunk]+'.sumwt.tt0"' )
00825                 ia.close()
00826             ## Normalize PSFs 
00827             ia.open( combpsfs[tt] )
00828             ia.calc('"'+ combpsfs[tt] + '"/"' + combwts[0] +'"')
00829             ia.done()
00830         #done
00831     
00832     def pcont(self, msname=None, imagename=None, imsize=[1000, 1000], 
00833               pixsize=['1arcsec', '1arcsec'], phasecenter='', 
00834               field='', spw='*', stokes='I', ftmachine='ft', wprojplanes=128, facets=1, 
00835               hostnames='',  
00836               numcpuperhost=1, majorcycles=-1, cyclefactor=1.5, niter=1000, npercycle=100, gain=0.1, threshold='0.0mJy', alg='clark', scales=[0], weight='natural', robust=0.0, npixels=0,  uvtaper=False, outertaper=[], timerange='', uvrange='', baselines='', scan='', observation='', pbcorr=False, minpb=0.2, 
00837               contclean=False, visinmem=False, interactive=False, maskimage='lala.mask',
00838               numthreads=1, savemodel=False,
00839               painc=360., pblimit=0.1, dopbcorr=True, applyoffsets=False, cfcache='cfcache.dir',
00840               epjtablename=''):
00841 
00842         """
00843         msname= measurementset
00844         imagename = image
00845         imsize = list of 2 numbers  [nx,ny] defining image size in x and y
00846         pixsize = list of 2 quantities   ['sizex', 'sizey'] defining the pixel size e.g  ['1arcsec', '1arcsec']
00847         phasecenter = an integer or a direction string   integer is fieldindex or direction e.g  'J2000 19h30m00 -30d00m00'
00848         field = field selection string ...msselection style
00849         spw = spw selection string ...msselection style
00850         stokes= string e.g 'I', 'IV'
00851         ftmachine= the ftmachine to use ...'ft', 'wproject' etc
00852         wprojplanes is an interger that is valid only of ftmachine is 'wproject', 
00853         facets= integer do split image facet, 
00854         hostnames= list of strings ..empty string mean localhost
00855         numcpuperhos = integer ...number of processes to launch on each host
00856         majorcycles= integer number of CS major cycles to do, 
00857         cyclefactor= if majorcycles=-1. This determines number of majorcycles
00858         niter= integer ...total number of clean iteration 
00859         threshold=quantity ...residual peak at which to stop deconvolving
00860         alg= string  possibilities are 'clark', 'hogbom', 'msclean'
00861         weight= string e.g 'natural', 'briggs' or 'radial'
00862         robust= float  valid for 'briggs' only  e.g 0,0
00863          
00864         scales = scales to use when using alg='msclean'
00865         contclean = boolean ...if False the imagename.model is deleted if its on 
00866         disk otherwise clean will continue from previous run
00867         visinmem = load visibility in memory for major cycles...make sure totalmemory  available to all processes is more than the MS size
00868         interactive boolean ...get a viewer to draw mask on
00869         maskimage  a prior mask image to limit clean search
00870         numthreads  number of threads to use in each engine
00871         savemodel   when True the model is saved in the MS header for selfcal
00872         painc = Parallactic angle increment in degrees after which a new convolution function is computed (default=360.0deg)
00873         cfcache = The disk cache directory for convolution functions
00874         pblimit = The fraction of the peak of the PB to which the PB corrections are applied (default=0.1)
00875         dopbcorr = If true, correct for PB in the major cycles as well
00876         applyoffsets = If true, apply antenna pointing offsets from the pointing table given by epjtablename 
00877         epjtablename = Table containing antenna pointing offsets
00878         """
00879 
00880 
00881         niterpercycle=niter/majorcycles if(majorcycles >0) else niter
00882         if(niterpercycle == 0):
00883             niterpercycle=niter
00884             majorcycles=1
00885         num_ext_procs=0
00886         self.setupcommonparams(spw=spw, field=field, phasecenter=phasecenter, 
00887                                stokes=stokes, ftmachine=ftmachine, wprojplanes=wprojplanes, 
00888                                facets=facets, imsize=imsize, pixsize=pixsize, weight=weight, 
00889                                robust=robust, npixels=npixels, gain=gain, uvtaper=uvtaper,
00890                                outertaper=outertaper, timerange=timerange, uvrange=uvrange, 
00891                                baselines=baselines, scan=scan, observation=observation, 
00892                                visinmem=visinmem, pbcorr=pbcorr, minpb=minpb, numthreads=numthreads, 
00893                                cyclefactor=cyclefactor)
00894         
00895         self.setupcluster(hostnames,numcpuperhost, num_ext_procs)
00896       
00897         if(spw==''):
00898             spw='*'
00899         if(field==''):
00900             field='*'
00901         spwids=ms.msseltoindex(vis=msname, spw=spw)['spw']
00902         timesplit=0
00903         timeimage=0
00904         elimageroot=imagename
00905         elmask=maskimage
00906         owd=os.getcwd()
00907         fullpath=lambda a: owd+'/'+a if ((len(a) !=0) and a[0] != '/') else a
00908         imagename=fullpath(elimageroot)
00909         maskimage=fullpath(elmask)
00910             
00911         model=imagename+'.model' if (len(elimageroot) != 0) else (owd+'/elmodel')
00912         tempmodel=owd+'/tempmodel'
00913         shutil.rmtree(tempmodel, True)
00914         if(not contclean):
00915             print "Removing ", model, 'and', imagename+'.image'
00916             shutil.rmtree(model, True)
00917             shutil.rmtree(imagename+'.image', True)
00918             shutil.rmtree(imagename+'.residual', True)
00919 
00920         ###num of cpu 
00921         numcpu=self.numcpu
00922         time1=time.time()
00923         ###spw and channel selection
00924         #spwsel,startsel,nchansel=self.findchanselcont(msname, spwids, numcpu)
00925 
00926         #print 'SPWSEL ', spwsel, startsel, nchansel
00927         freqrange=[0.0, 0.0]
00928         spwsel=self.findchansel(msname, spw, field, numcpu, freqrange=freqrange)
00929         minfreq=freqrange[0]
00930         maxfreq=freqrange[1]
00931 
00932         out=range(numcpu) 
00933         
00934         c=self.c
00935        ####
00936         #the default working directory is somewhere 
00937         owd=os.getcwd()
00938         self.c.pgc('import os')
00939         self.c.pgc('os.chdir("'+owd+'")')
00940         ##################### 
00941         c.pgc('casalog.filter()')
00942         c.pgc('from  parallel.parallel_cont import *')
00943 
00944         band=maxfreq-minfreq
00945         ###need to get the data selection for each process here
00946         ## this algorithm is a poor first try
00947         if(minfreq <0 ):
00948             minfreq=0.0
00949         freq='"%s"'%(str((minfreq+band/2.0))+'Hz')
00950         band='"%s"'%(str((band*1.1))+'Hz')
00951         ###define image names
00952         imlist=[]
00953         char_set = string.ascii_letters
00954         cfcachelist=[]
00955         for k in range(numcpu):
00956             substr=self.workingdirs[k]+'/Temp_'+str(k)+'_'+string.join(random.sample(char_set,8), sep='')
00957             while (os.path.exists(substr+'.model')):
00958                 substr=self.workingdirs[k]+'/Temp_'+str(k)+'_'+string.join(random.sample(char_set,8), sep='')
00959             ###############
00960             imlist.append(substr)
00961             cfcachelist.append(cfcache+'_'+str(k));
00962         #####
00963         ##continue clean or not
00964         if(contclean and os.path.exists(model)):
00965             for k in range(numcpu):
00966                 self.copyimage(inimage=model, outimage=imlist[k]+'.model', init=False)
00967         else:
00968             for k in range(len(imlist)):
00969                 shutil.rmtree(imlist[k]+'.model', True)
00970         intmask=0
00971         donewgt=[False]*numcpu
00972         ### do one major cycle at the end
00973         donemaj=False
00974         maj=0;
00975         maxresid=-1.0
00976         while(not donemaj):
00977             
00978             casalog.post('Starting Gridding for major cycle '+str(maj)) 
00979             for k in range(numcpu):
00980                 imnam='"%s"'%(imlist[k])
00981                 c.odo('a.cfcache='+'"'+str(cfcachelist[k])+'"',k);
00982                 c.odo('a.painc='+str(painc),k);
00983                 c.odo('a.pblimit='+str(pblimit),k);
00984                 c.odo('a.dopbcorr='+str(dopbcorr),k);
00985                 c.odo('a.applyoffsets='+str(applyoffsets),k);
00986                 c.odo('a.epjtablename='+'"'+str(epjtablename)+'"',k);
00987                 runcomm='a.imagecont(msname='+'"'+msname+'", field="'+str(field)+'", spw="'+str(spwsel[k])+'", freq='+freq+', band='+band+', imname='+imnam+')'
00988                 print 'command is ', runcomm
00989                 out[k]=c.odo(runcomm,k)
00990             over=False
00991             printkounter=0
00992             while(not over):
00993                 time.sleep(5)
00994                 #printkounter +=1
00995                 overone=True
00996                 for k in range(numcpu):
00997                     #print 'k', k, out[k]
00998                     cj=c.check_job(out[k],False)
00999                     overone=(overone and cj)
01000                     #print 'maj, cj, dwght', maj, cj, donewgt[k]
01001                     if((maj==0) and cj and (not donewgt[k])):
01002                         print 'doing weight', 'a.getftweight(msname="'+msname+'", wgtimage="'+imlist[k]+'.wgt''")'
01003                         out[k]=c.odo('a.getftweight(msname="'+msname+'", wgtimage="'+imlist[k]+'.wgt''")', k)
01004                         donewgt[k]=True
01005                         while(not c.check_job(out[k],False)):
01006                             time.sleep(1)
01007                     
01008                     #if((printkounter==10) and not(c.check_job(out[k],False))):
01009                     #    print 'job ', k, 'is waiting'
01010                     #    printkounter=0
01011                 over=overone
01012             residual=imagename+'.residual'
01013             psf=imagename+'.psf'
01014             fluxim=imagename+'.flux'
01015             coverim=imagename+'.flux.pbcoverage'
01016             psfs=range(len(imlist))
01017             residuals=range(len(imlist))
01018             restoreds=range(len(imlist))
01019             weightims=range(len(imlist))
01020             fluxims=range(len(imlist))
01021             coverims=range(len(imlist))
01022             for k in range (len(imlist)):
01023                 psfs[k]=imlist[k]+'.psf'
01024                 residuals[k]=imlist[k]+'.residual'
01025                 restoreds[k]=imlist[k]+'.image'
01026                 fluxims[k]=imlist[k]+'.flux'
01027                 weightims[k]=imlist[k]+'.wgt'
01028                 coverims[k]=imlist[k]+'.flux.pbcoverage'
01029             self.weightedaverimages(residual, residuals, weightims)
01030             if((maskimage == '') or (maskimage==[])):
01031                 maskimage=imagename+'.mask'
01032                 ia.removefile(maskimage)
01033             if (interactive and (intmask==0)):
01034                 if(maj==0):
01035                     ia.removefile(maskimage)
01036                 retdraw=im.drawmask(imagename+'.residual',maskimage, niter=niterpercycle, npercycle=npercycle, threshold=threshold);
01037                 intmask=retdraw['stat']
01038                     ####The user may want to change these in the gui
01039                 threshold=retdraw['threshold']
01040                 niterpercycle=retdraw['niter']
01041                 npercycle=retdraw['npercycle']
01042                 print 'intmask', intmask
01043                 if(intmask==1):
01044                     interactive=False
01045                     im.done()
01046                 if(intmask==2):
01047                     interactive=False
01048                     break;
01049             else:
01050                 if (maj==0):
01051                     if(os.path.exists(maskimage)):
01052                         self.regridimage(outimage='__lala.mask', inimage=maskimage, templateimage=residual);
01053                         shutil.rmtree(maskimage, True)
01054                         shutil.move('__lala.mask', maskimage)
01055                     else:
01056                         print 'DOING a full image mask'
01057                         self.copyimage(inimage=residual, outimage=maskimage, init=True, initval=1.0);
01058            #########Things to be done after first major cycle only
01059            #########
01060             if(maj==0):
01061     #            copyimage(inimage=residual, outimage='lala.mask', init=True, initval=1.0)
01062                 if(not contclean or (not os.path.exists(model))):
01063                     self.copyimage(inimage=residual, outimage=model, 
01064                               init=True, initval=0.0)     
01065                 self.weightedaverimages(psf, psfs, weightims)
01066                 if(self.ftmachine=='mosaic'):
01067                     self.averimages(fluxim, fluxims)
01068                     self.averimages(coverim, coverims)
01069                     ##OR the mask and the fluximage into the mask
01070                     elim=casac.image()
01071                     elim.open(maskimage)
01072                     elim.calc('iif(mask("'+fluxim+'"), "'+maskimage+'", 0)')
01073                     elim.done()
01074                 if((self.weight=='uniform') or (self.weight=='briggs')):
01075                     c.pgc('wtgrid=a.getweightgrid(msname="'+msname+'")')
01076                     sumweight=c.pull('wtgrid', 0)[0]
01077                     for jj in range(1,numcpu):
01078                         sumweight += c.pull('wtgrid', jj)[jj]
01079                         print 'Max min var of weight dens', np.max(sumweight), np.min(sumweight), np.var(sumweight)
01080                     c.push(wtgrid=sumweight)
01081                     c.pgc('a.setweightgrid(msname="'+msname+'", weight=wtgrid)')
01082             newthresh=threshold
01083             if(majorcycles <= 1):
01084                 ia.open(residual)
01085                 residstat=ia.statistics()
01086                 maxresid=np.max(residstat['max'], np.fabs(residstat['min']))
01087                 psfoutermax=self.maxouterpsf(psfim=psf, image=restoreds[0])
01088                 newthresh=psfoutermax*cyclefactor*maxresid
01089                 oldthresh=qa.convert(qa.quantity(threshold, "Jy"), "Jy")['value']
01090                 if(newthresh < oldthresh):
01091                     newthresh=oldthresh
01092                 newthresh=qa.tos(qa.quantity(newthresh, "Jy"))
01093             ####no need to do this in last major cycle
01094             needclean=((maj < majorcycles if (majorcycles >1) else True)  and 
01095                        ((maxresid > qa.convert(qa.quantity(threshold),'Jy')['value']) if(majorcycles <2) else True) 
01096                        and (niterpercycle > 1))
01097             donemaj = not needclean
01098             if(needclean):
01099                 casalog.post('Deconvolving for major cycle '+str(maj))  
01100                 #incremental clean...get rid of tempmodel
01101                 shutil.rmtree(tempmodel, True)
01102                 elniter=npercycle if(interactive) else niterpercycle
01103                 rundecon='retval=a.cleancont(alg="'+str(alg)+'", thr="'+str(newthresh)+'", scales='+ str(scales)+', niter='+str(elniter)+',psf="'+psf+'", dirty="'+residual+'", model="'+tempmodel+'", mask="'+str(maskimage)+'")'
01104                 print 'Deconvolution command', rundecon
01105                 out[0]=c.odo(rundecon,0)
01106                 over=False
01107                 while(not over):
01108                     time.sleep(5)
01109                     over=c.check_job(out[0],False)
01110                 retval=c.pull('retval', 0)[0]
01111                 if(majorcycles <=1):
01112                     niterpercycle=niterpercycle-retval['iterations']
01113                 maxresid=retval['maxresidual']
01114                 ###incremental added to total 
01115                 ia.open(model)
01116                 ia.calc('"'+model+'" +  "'+tempmodel+'"')
01117                 ia.done()
01118                 ia.open(tempmodel)
01119             #arr=ia.getchunk()
01120                 imminmax= abs(ia.statistics()['min'])
01121                 imminmax=max(imminmax, ia.statistics()['max'])
01122                 print 'min max of incrmodel', ia.statistics()['min'], ia.statistics()['max']
01123                 ia.done()
01124                 if(imminmax == 0.0):
01125                     print 'Threshold reached'
01126                     break
01127                 ia.open(model)
01128             #arr=ia.getchunk()
01129                 print 'min max of model', ia.statistics()['min'], ia.statistics()['max']
01130                 ia.done()
01131                 for k in range(len(imlist)):
01132                     ia.open(imlist[k]+'.model')
01133                 #ia.putchunk(arr)
01134                     ia.insert(infile=model, locate=[0,0,0,0])
01135                     ia.done()
01136                 maj +=1
01137                 #if(majorcycles > 1):
01138                 #   donemaj=(maj >= majorcycles)
01139             
01140         
01141         restored=imagename+'.image'
01142         #ia.open(restored)
01143         #for k in range(1, len(imlist)) :
01144         #    ia.calc('"'+restored+'" +  "'+imlist[k]+'.image"')
01145         #ia.calc('"'+restored+'"'+'/'+str(len(imlist)))
01146         #ia.done()
01147         shutil.rmtree(imagename+'.image', True)
01148         self.weightedaverimages(restored, restoreds, weightims)
01149         #shutil.move(restored,  imagename+'.image')
01150         time2=time.time()
01151         ###Clean up
01152         for k in range(len(imlist)):
01153             shutil.rmtree(imlist[k]+'.model', True)
01154             shutil.rmtree(imlist[k]+'.psf', True)
01155             shutil.rmtree(imlist[k]+'.residual', True)
01156             shutil.rmtree(imlist[k]+'.image', True)
01157             shutil.rmtree(imlist[k]+'.wgt', True)
01158         if(savemodel):
01159             myim=casac.imager()
01160             myim.selectvis(vis=msname, spw=spw, field=field)
01161             myim.defineimage()
01162             myim.setoptions(ftmachine=ftmachine,wprojplanes=wprojplanes)
01163             myim.ft(model)
01164             myim.done()
01165         print 'Time to image is ', (time2-time1)/60.0, 'mins'
01166         casalog.post('Time to image is '+str((time2-time1)/60.0)+ ' mins')
01167         c.pgc('del a')
01168         #c.stop_cluster()
01169 
01170     def pcube_try(self, msname=None, imagename='elimage', imsize=[1000, 1000], 
01171               pixsize=['1arcsec', '1arcsec'], phasecenter='', 
01172               field='', spw='*', ftmachine='ft', wprojplanes=128, facets=1, 
01173               hostnames='', 
01174               numcpuperhost=1, majorcycles=1, niter=1000, gain=0.1, threshold='0.0mJy', alg='clark', scales=[0],
01175               mode='channel', start=0, nchan=1, step=1, restfreq='', weight='natural', 
01176               robust=0.0, npixels=0, 
01177               imagetilevol=100000,
01178               contclean=False, pbcorr=False, chanchunk=1, visinmem=False, maskimage='' , numthreads=-1,
01179               painc=360., pblimit=0.1, dopbcorr=True, applyoffsets=False, cfcache='cfcache.dir',
01180               epjtablename=''): 
01181 
01182         """
01183         msname= measurementset
01184         imagename = image
01185         imsize = list of 2 numbers  [nx,ny] defining image size in x and y
01186         pixsize = list of 2 quantities   ['sizex', 'sizey'] defining the pixel size e.g  ['1arcsec', '1arcsec']
01187         phasecenter = an integer or a direction string   integer is fieldindex or direction e.g  'J2000 19h30m00 -30d00m00'
01188         field = field selection string ...msselection style
01189         spw = spw selection string ...msselection style
01190         ftmachine= the ftmachine to use ...'ft', 'wproject' etc
01191         wprojplanes is an interger that is valid only of ftmachine is 'wproject', 
01192         facets= integer do split image facet, 
01193         hostnames= list of strings ..empty string mean localhost
01194         numcpuperhos = integer ...number of processes to launch on each host
01195         majorcycles= integer number of CS major cycles to do, 
01196         niter= integer ...total number of clean iteration 
01197         threshold=quantity string ...residual peak at which to stop deconvolving
01198         alg= string  possibilities are 'clark', 'hogbom', 'multiscale' and their 'mf'
01199         scales= list of scales in pixel for multiscale clean e.g [0, 3, 10]
01200         mode= channel definition, can be 'channel', 'frequency', 'velocity'
01201         start = first channel in the definition spec of mode, can be int, freq or vel quantity
01202         step = channel width specified in the definition of mode
01203         restfreq= what 'rest frequency' to use to calculate velocity from frequency
01204                                 empty string '' implies use the first restfreq in SOURCE of ms 
01205         weight= type of weight to apply
01206         contclean = boolean ...if False the imagename.model is deleted if its on 
01207         disk otherwise clean will continue from previous run
01208         chanchunk = number of channel to process at a go per process...careful not to 
01209        go above total memory available
01210        visinmem = load visibility in memory for major cycles...make sure totalmemory  available to all processes is more than the MS size
01211         painc = Parallactic angle increment in degrees after which a new convolution function is computed (default=360.0deg)
01212         cfcache = The disk cache directory for convolution functions
01213         pblimit = The fraction of the peak of the PB to which the PB corrections are applied (default=0.1)
01214         dopbcorr = If true, correct for PB in the major cycles as well
01215         applyoffsets = If true, apply antenna pointing offsets from the pointing table given by epjtablename 
01216         epjtablename = Table containing antenna pointing offsets
01217         """
01218 
01219         if(spw==''):
01220             spw='*'
01221         if(field==''):
01222             field='*'
01223         spwids=ms.msseltoindex(vis=msname, spw=spw)['spw']
01224         ###num of cpu per node
01225         numcpu=numcpuperhost
01226         time1=time.time()
01227         self.setupcommonparams(spw=spw, field=field, phasecenter=phasecenter, 
01228                                stokes=stokes, ftmachine=ftmachine, wprojplanes=wprojplanes, 
01229                                facets=facets, imsize=imsize, pixsize=pixsize, weight=weight, 
01230                                robust=robust, npixels=npixels, gain=gain,   
01231                                visinmem=visinmem, pbcorr=pbcorr, numthreads=numthreads)
01232         self.setupcluster(hostnames,numcpuperhost, 3)
01233         numcpu=self.numcpu
01234         ##Start an slave for my async use for cleaning up etc here
01235         buddy_id=[numcpu, numcpu+1, numcpu+2]
01236         self.c.push(numcpu=numcpu, targets=buddy_id) 
01237         #####################
01238         ###set the working directory here...
01239         owd=os.getcwd()
01240         self.c.pgc('import os')
01241         self.c.pgc('os.chdir("'+owd+'")')
01242         elimageroot=imagename
01243         #elmask=maskimage
01244         #fullpath=lambda a: owd+'/'+a if ((len(a) !=0) and a[0] != '/') else a
01245         #imagename=fullpath(elimageroot)
01246         #maskimage=fullpath(elmask)           
01247         model=imagename+'.model' if (len(elimageroot) != 0) else (owd+'/elmodel')
01248         if(not contclean or (not os.path.exists(model))):
01249             shutil.rmtree(model, True)
01250             shutil.rmtree(imagename+'.image', True)
01251             ##create the cube
01252             im.selectvis(vis=msname, spw=spw, field=field, writeaccess=False)
01253             im.defineimage(nx=imsize[0], ny=imsize[1], cellx=pixsize[0], celly=pixsize[1], 
01254                            phasecenter=phasecenter, mode=mode, spw=spwids.tolist(), nchan=nchan, step=step, start=start, restfreq=restfreq)
01255             im.setoptions(imagetilevol=imagetilevol) 
01256             print 'making model image (', model, ') ...'
01257             im.make(model)
01258             print 'model image (', model, ') made'
01259             im.done()
01260         #print 'LOCKS ', tb.listlocks()
01261         ia.open(model)
01262         csys=ia.coordsys()
01263         elshape=ia.shape()
01264         ia.done()
01265         if((maskimage != '') and (os.path.exists(maskimage))):
01266             ia.open(maskimage)
01267             maskshape=ia.shape()
01268             ia.done()
01269             if(np.any(maskshape != elshape)):
01270                 newmask=maskimage+'_regrid'
01271                 self.regridimage(outimage=newmask, inimage=maskimage, templateimage=model);
01272                 maskimage=newmask
01273         ###as image will have conversion to LSRK...need to get original stuff
01274         #originsptype=csys.getconversiontype('spectral', showconversion=False)
01275         #csys.setconversiontype(spectral=originsptype)
01276         fstart=csys.toworld([0,0,0,0],'n')['numeric'][3]
01277         fstep=csys.toworld([0,0,0,1],'n')['numeric'][3]-fstart
01278         fend=fstep*(nchan-1)+fstart
01279         
01280         #print 'LOCKS2 ', tb.listlocks()
01281         imepoch=csys.epoch()
01282         imobservatory=csys.telescope()
01283         shutil.rmtree(imagename+'.image', True)
01284         shutil.rmtree(imagename+'.residual', True)
01285         shutil.copytree(model, imagename+'.image')
01286         shutil.copytree(model, imagename+'.residual')
01287 
01288         out=range(numcpu)  
01289         self.c.pgc('from  parallel.parallel_cont import *')
01290         #spwlaunch='"'+spw+'"' if (type(spw)==str) else str(spw)
01291         #fieldlaunch='"'+field+'"' if (type(field) == str) else str(field)
01292         #pslaunch='"'+phasecenter+'"' if (type(phasecenter) == str) else str(phasecenter)
01293         #launchcomm='a=imagecont(ftmachine='+'"'+ftmachine+'",'+'wprojplanes='+str(wprojplanes)+',facets='+str(facets)+',pixels='+str(imsize)+',cell='+str(pixsize)+', spw='+spwlaunch +',field='+fieldlaunch+',phasecenter='+pslaunch+',weight="'+weight+'")'
01294         #print 'launch command', launchcomm
01295         #self.c.pgc(launchcomm)
01296         ###set some common parameters
01297         self.c.pgc('a.imagetilevol='+str(imagetilevol))
01298         self.c.pgc('a.visInMem='+str(visinmem))
01299         self.c.pgc('a.painc='+str(painc))
01300         self.c.pgc('a.cfcache='+'"'+str(cfcache)+'"')
01301         self.c.pgc('a.pblimit='+str(pblimit));
01302         self.c.pgc('a.dopbcorr='+str(dopbcorr));
01303         self.c.pgc('a.applyoffsets='+str(applyoffsets));
01304         self.c.pgc('a.epjtablename='+'"'+str(epjtablename)+'"');
01305 
01306         tb.clearlocks()
01307         #print 'LOCKS3', tb.listlocks()
01308         chancounter=0
01309         nchanchunk=nchan/chanchunk if (nchan%chanchunk) ==0 else nchan/chanchunk+1
01310         spwsel,startsel, nchansel=imagecont.findchanselLSRK(msname=msname, spw=spwids, 
01311                                                       field=field, 
01312                                                       numpartition=nchanchunk, 
01313                                                       beginfreq=fstart, endfreq=fend, chanwidth=fstep)
01314         print 'spwsel', spwsel, 'startsel', startsel,'nchansel', nchansel
01315         #print 'startsel', startsel
01316         #print  'nchansel', nchansel
01317         imnam='"%s"'%(imagename)
01318         donegetchan=np.array(range(nchanchunk),dtype=bool)
01319         doneputchan=np.array(range(nchanchunk),dtype=bool)
01320         readyputchan=np.array(range(nchanchunk), dtype=bool)
01321         cpudoing=np.array(range(nchanchunk), dtype=int)
01322         donegetchan.setfield(False,bool)
01323         doneputchan.setfield(False,bool)
01324         readyputchan.setfield(False, bool)
01325         chanind=np.array(range(numcpu), dtype=int)
01326         self.c.push(readyputchan=readyputchan, targets=buddy_id)
01327         #c.push(doneputchan=doneputchan, targets=buddy_id)
01328         buddy_is_ready=[True, True, True]
01329         buddy_ref=[False, False, False]
01330         cleanupcomm=['', '', '']
01331         cleanupcomm[2]='a.cleanupmodelimages(readyputchan=readyputchan,  imagename='+imnam+', nchanchunk='+str(nchanchunk)+', chanchunk='+str(chanchunk)+')'
01332         cleanupcomm[1]='a.cleanupresidualimages(readyputchan=readyputchan,  imagename='+imnam+', nchanchunk='+str(nchanchunk)+', chanchunk='+str(chanchunk)+')'
01333         cleanupcomm[0]='a.cleanuprestoredimages(readyputchan=readyputchan,  imagename='+imnam+', nchanchunk='+str(nchanchunk)+', chanchunk='+str(chanchunk)+')'
01334         def gen_command(ccounter):
01335             return 'a.imagechan_new(msname='+'"'+msname+'", start='+str(startsel[ccounter])+', numchan='+str(nchansel[ccounter])+', field="'+str(field)+'", spw='+str(spwsel[ccounter])+', cubeim='+imnam+', imroot='+imnam+',imchan='+str(ccounter)+',chanchunk='+str(chanchunk)+',niter='+str(niter)+',alg="'+alg+'", scales='+str(scales)+', majcycle='+str(majorcycles)+', thr="'+str(threshold)+'", mask="'+maskimage+'")'
01336 
01337         #while(chancounter < nchanchunk):
01338         chanind.setfield(-1, int)
01339         for k in range(numcpu):
01340             if(chancounter < nchanchunk):
01341                 chanind[k]=chancounter
01342                 runcomm=gen_command(chancounter)
01343                 print 'command is ', runcomm
01344                 out[k]=self.c.odo(runcomm,k)
01345                 chancounter=chancounter+1
01346         while(chancounter < nchanchunk):
01347                 over=False
01348                 while(not over):
01349                 #############loop waiting for a chunk of work
01350                     time.sleep(1)
01351                     for bud in range(3):
01352                         if(buddy_is_ready[bud]):
01353                             print 'SENDING ', cleanupcomm[bud]
01354                             self.c.push(readyputchan=readyputchan, targets=buddy_id[bud])
01355                         #c.push(doneputchan=doneputchan, targets=buddy_id)
01356                             buddy_ref[bud]=self.c.odo(cleanupcomm[bud], buddy_id[bud])
01357                         buddy_is_ready[bud]=self.c.check_job(buddy_ref[bud], False)
01358                         #print 'buddy_ready', bud, buddy_is_ready[bud]
01359                 #if(buddy_is_ready):
01360                 #    doneputchan=c.pull('doneputchan', buddy_id)[buddy_id]
01361                     overone=True
01362                     for k in range(numcpu):
01363                         overone=(overone and self.c.check_job(out[k],False))
01364                         if((chanind[k] > -1) and self.c.check_job(out[k],False) and 
01365                            (not readyputchan[chanind[k]])):
01366                             readyputchan[chanind[k]]=True      
01367                             if(chancounter < nchanchunk):
01368                                 chanind[k]=chancounter
01369                                 runcomm=gen_command(chancounter)
01370                                 print 'command is ', runcomm
01371                                 print 'processor ', k
01372                                 out[k]=self.c.odo(runcomm,k)
01373                                 chancounter+=1
01374                             overone=(overone and self.c.check_job(out[k],False))
01375                     over=overone
01376                ############
01377         time2=time.time()
01378         print 'Time to image is ', (time2-time1)/60.0, 'mins'
01379         ##sweep the remainder channels in case they are missed
01380         for bud in range(3):
01381             while(not buddy_is_ready[bud]):
01382                 buddy_is_ready[bud]=self.c.check_job(buddy_ref[bud], False)
01383             #doneputchan=c.pull('doneputchan', buddy_id)[buddy_id] 
01384             self.c.push(readyputchan=readyputchan, targets=buddy_id[bud])
01385             buddy_ref[bud]=self.c.odo(cleanupcomm[bud], buddy_id[bud])
01386         for bud in range(3):
01387             while(not buddy_is_ready[bud]):
01388                 buddy_is_ready[bud]=self.c.check_job(buddy_ref[bud], False)
01389         time2=time.time()
01390         print 'Time to image after cleaning is ', (time2-time1)/60.0, 'mins'
01391         #self.c.stop_cluster()
01392 
01393 #################
01394     def pcube_driver(self, msname=None, imagename='elimage', imsize=[1000, 1000], 
01395               pixsize=['1arcsec', '1arcsec'], phasecenter='', 
01396               field='', spw='*', ftmachine='ft', wprojplanes=128, facets=1, 
01397               hostnames='', 
01398               numcpuperhost=1, majorcycles=-1, cyclefactor=1.5, niter=1000, npercycle=100, gain=0.1, threshold='0.0mJy', alg='clark', scales=[0],
01399               mode='channel', start=0, nchan=1, step=1, restfreq='', stokes='I', weight='natural', 
01400               robust=0.0, npixels=0,uvtaper=False, outertaper=[], timerange='', uvrange='',baselines='', scan='', observation='',  pbcorr=False,  minpb=0.2,
01401               imagetilevol=100000,
01402               contclean=False, chanchunk=1, visinmem=False, maskimage='' , numthreads=1,  savemodel=False,
01403               painc=360., pblimit=0.1, dopbcorr=True, applyoffsets=False, cfcache='cfcache.dir',
01404               epjtablename='', interactive=False):
01405         """
01406         Drive pcube when interactive is on and depending on cyclefactor or not
01407         """
01408         if(interactive==False):
01409             self.pcube(msname=msname, imagename=imagename, imsize=imsize, 
01410               pixsize=pixsize, phasecenter=phasecenter, 
01411               field=field, spw=spw, ftmachine=ftmachine, wprojplanes=wprojplanes, facets=facets, 
01412               hostnames=hostnames, 
01413               numcpuperhost=numcpuperhost, majorcycles=majorcycles, cyclefactor=cyclefactor, niter=niter, gain=gain, threshold=threshold, alg=alg, scales=scales,
01414               mode=mode, start=start, nchan=nchan, step=step, restfreq=restfreq, stokes=stokes, weight=weight, 
01415               robust=robust, npixels=npixels,uvtaper=uvtaper, outertaper=outertaper, timerange=timerange, uvrange=uvrange, baselines=baselines, scan=scan, 
01416                   observation=observation,  pbcorr=pbcorr, minpb=minpb, imagetilevol=imagetilevol,
01417               contclean=contclean, chanchunk=chanchunk, visinmem=visinmem, maskimage=maskimage , numthreads=numthreads,  savemodel=savemodel,
01418               painc=painc, pblimit=pblimit, dopbcorr=dopbcorr, applyoffsets=applyoffsets, cfcache=cfcache, epjtablename=epjtablename)
01419             return
01420         ###interactive is true
01421         ###lets get the 0th iteration niter=0, majorcycles=1, maskimage='', savemodel=False
01422         self.pcube(msname=msname, imagename=imagename, imsize=imsize, pixsize=pixsize, phasecenter=phasecenter, 
01423               field=field, spw=spw, ftmachine=ftmachine, wprojplanes=wprojplanes, facets=facets, hostnames=hostnames, 
01424               numcpuperhost=numcpuperhost, majorcycles=1, cyclefactor=cyclefactor, niter=0, gain=gain, threshold=threshold, alg=alg, scales=scales,
01425               mode=mode, start=start, nchan=nchan, step=step, restfreq=restfreq, stokes=stokes, weight=weight, 
01426               robust=robust, npixels=npixels,uvtaper=uvtaper, outertaper=outertaper, timerange=timerange, uvrange=uvrange, baselines=baselines, scan=scan, 
01427                   observation=observation,  pbcorr=pbcorr, minpb=minpb, imagetilevol=imagetilevol,
01428               contclean=contclean, chanchunk=chanchunk, visinmem=visinmem, maskimage='', numthreads=numthreads,  savemodel=False,
01429               painc=painc, pblimit=pblimit, dopbcorr=dopbcorr, applyoffsets=applyoffsets, cfcache=cfcache, epjtablename=epjtablename)
01430         if(maskimage==''):
01431             maskimage=imagename+'.mask'
01432         residual=imagename+'.residual'
01433         ###now deal with interactive and    
01434         if(cyclefactor > 0.00000001):
01435             ####after this conclean has to be true
01436             psf=imagename+'.psf'
01437             newthresh=threshold
01438             psfoutermax=self.maxouterpsf(psfim=psf, image=imagename+'.image')
01439             oldthresh=qa.convert(qa.quantity(threshold, "Jy"), "Jy")['value']
01440             ###have to get the niter back from pcube
01441             #myim,=gentools(['im'])
01442             #retval=myim.drawmask(imagename+'.residual', maskimage)
01443             #myim.done()
01444             ###for now if oldthresh is 0 that is niter is determines end we will do one interactive loop
01445             notdone=True
01446             while(notdone):
01447                 
01448      
01449                 ia.open(residual)
01450                 residstat=ia.statistics()
01451                 ia.done()
01452                 maxresid=np.max(residstat['max'], np.fabs(residstat['min']))
01453                 newthresh=psfoutermax*cyclefactor*maxresid
01454                 if(newthresh < oldthresh):
01455                     newthresh=oldthresh
01456                 myim=casac.imager()
01457                 retdraw=myim.drawmask(imagename+'.residual',maskimage, niter=niter, npercycle=npercycle, threshold=qa.tos(qa.quantity(newthresh, "Jy")));
01458                 myim.done()
01459                 retintval=retdraw['stat']
01460                     ####The user may want to change these in the gui
01461                 somethresh=qa.convert(qa.quantity(retdraw['threshold'], "Jy"), "Jy")['value']
01462                 if(somethresh < newthresh):
01463                     newthresh=somethresh
01464                 niter=retdraw['niter']
01465                 npercycle=retdraw['npercycle']
01466                 if(retintval==2):
01467                     notdone=False
01468                     break
01469                 if(retintval==1):
01470                     ###The user has decided to be done with interactive
01471                     notdone=False
01472                     newthresh=oldthresh
01473                 newthresh=qa.tos(qa.quantity(newthresh, "Jy"))
01474                 threshold=newthresh                
01475                 retval=self.pcube(msname=msname, imagename=imagename, imsize=imsize, pixsize=pixsize, phasecenter=phasecenter, 
01476                            field=field, spw=spw, ftmachine=ftmachine, wprojplanes=wprojplanes, facets=facets, hostnames=hostnames, 
01477                            numcpuperhost=numcpuperhost, majorcycles=majorcycles, cyclefactor=cyclefactor, niter=npercycle, gain=gain, threshold=threshold, alg=alg, scales=scales,
01478                            mode=mode, start=start, nchan=nchan, step=step, restfreq=restfreq, stokes=stokes, weight=weight, 
01479                            robust=robust, npixels=npixels,uvtaper=uvtaper, outertaper=outertaper, timerange=timerange, uvrange=uvrange, baselines=baselines, scan=scan, 
01480                            observation=observation,  pbcorr=pbcorr, minpb=minpb, imagetilevol=imagetilevol,
01481                            contclean=True, chanchunk=chanchunk, visinmem=visinmem, maskimage=maskimage , numthreads=numthreads,  savemodel=False,
01482                            painc=painc, pblimit=pblimit, dopbcorr=dopbcorr, applyoffsets=applyoffsets, cfcache=cfcache, epjtablename=epjtablename)
01483                 ####Here we should extract the remainder iterations
01484                 niter=niter-retval['iterations']
01485                 print 'retval', retval, retval['maxresidual'], oldthresh
01486                 if ((niter < 1) or (retval['maxresidual'] < oldthresh)) :
01487                     notdone=False
01488         else: 
01489             #using majorcycles
01490             npercycle=niter/majorcycles
01491             for k in range(majorcycles):
01492                 myim=casac.imager()
01493                 retdraw=myim.drawmask(imagename+'.residual',maskimage, niter=niter, npercycle=npercycle, threshold=threshold);
01494                 myim.done()
01495                 niter=retdraw['niter']
01496                 npercycle=retdraw['npercycle']
01497                 threshold=retdraw['threshold']
01498                 retval=retdraw['stat']
01499                 if(retval==2):
01500                     break
01501                 self.pcube(msname=msname, imagename=imagename, imsize=imsize, pixsize=pixsize, phasecenter=phasecenter, 
01502                            field=field, spw=spw, ftmachine=ftmachine, wprojplanes=wprojplanes, facets=facets, hostnames=hostnames, 
01503                            numcpuperhost=numcpuperhost, majorcycles=1, cyclefactor=0, niter=npercycle, gain=gain, threshold=threshold, alg=alg, scales=scales,
01504                            mode=mode, start=start, nchan=nchan, step=step, restfreq=restfreq, stokes=stokes, weight=weight, 
01505                            robust=robust, npixels=npixels,uvtaper=uvtaper, outertaper=outertaper, timerange=timerange, uvrange=uvrange, baselines=baselines, scan=scan, 
01506                            observation=observation,  pbcorr=pbcorr, minpb=minpb, imagetilevol=imagetilevol,
01507                            contclean=True, chanchunk=chanchunk, visinmem=visinmem, maskimage=maskimage , numthreads=numthreads,  savemodel=False,
01508                            painc=painc, pblimit=pblimit, dopbcorr=dopbcorr, applyoffsets=applyoffsets, cfcache=cfcache, epjtablename=epjtablename)
01509         if(savemodel):
01510             myim=casac.imager()
01511             myim.selectvis(vis=msname, spw=spw, field=field)
01512             myim.defineimage()
01513             myim.setoptions(ftmachine=ftmachine,wprojplanes=wprojplanes)
01514             myim.ft(imagename+'.model')
01515             myim.done()
01516             
01517                 
01518 ############################################
01519     def pcube(self, msname=None, imagename='elimage', imsize=[1000, 1000], 
01520               pixsize=['1arcsec', '1arcsec'], phasecenter='', 
01521               field='', spw='*', ftmachine='ft', wprojplanes=128, facets=1, 
01522               hostnames='', 
01523               numcpuperhost=1, majorcycles=-1, cyclefactor=1.5, niter=1000, gain=0.1, threshold='0.0mJy', alg='clark', scales=[0],
01524               mode='channel', start=0, nchan=1, step=1, restfreq='', stokes='I', weight='natural', 
01525               robust=0.0, npixels=0,uvtaper=False, outertaper=[], timerange='', uvrange='',baselines='', scan='', observation='',  pbcorr=False,  minpb=0.2, 
01526               imagetilevol=100000,
01527               contclean=False, chanchunk=1, visinmem=False, maskimage='' , numthreads=1,  savemodel=False,
01528               painc=360., pblimit=0.1, dopbcorr=True, applyoffsets=False, cfcache='cfcache.dir',
01529               epjtablename=''): 
01530 
01531         """
01532         msname= measurementset
01533         imagename = image
01534         imsize = list of 2 numbers  [nx,ny] defining image size in x and y
01535         pixsize = list of 2 quantities   ['sizex', 'sizey'] defining the pixel size e.g  ['1arcsec', '1arcsec']
01536         phasecenter = an integer or a direction string   integer is fieldindex or direction e.g  'J2000 19h30m00 -30d00m00'
01537         field = field selection string ...msselection style
01538         spw = spw selection string ...msselection style
01539         ftmachine= the ftmachine to use ...'ft', 'wproject' etc
01540         wprojplanes is an interger that is valid only of ftmachine is 'wproject', 
01541         facets= integer do split image facet, 
01542         hostnames= list of strings ..empty string mean localhost
01543         numcpuperhos = integer ...number of processes to launch on each host
01544         majorcycles= integer number of CS major cycles to do, 
01545         niter= integer ...total number of clean iteration 
01546         threshold=quantity string ...residual peak at which to stop deconvolving
01547         alg= string  possibilities are 'clark', 'hogbom', 'multiscale' and their 'mf'
01548         scales= list of scales in pixel for multiscale clean e.g [0, 3, 10]
01549         mode= channel definition, can be 'channel', 'frequency', 'velocity'
01550         start = first channel in the definition spec of mode, can be int, freq or vel quantity
01551         step = channel width specified in the definition of mode
01552         restfreq= what 'rest frequency' to use to calculate velocity from frequency
01553                                 empty string '' implies use the first restfreq in SOURCE of ms 
01554         weight= type of weight to apply
01555         contclean = boolean ...if False the imagename.model is deleted if its on 
01556         disk otherwise clean will continue from previous run
01557         chanchunk = number of channel to process at a go per process...careful not to 
01558        go above total memory available
01559        visinmem = load visibility in memory for major cycles...make sure totalmemory  available to all processes is more than the MS size
01560        numthreads number of threads to use per engine while gridding
01561        savemodel if True save the model in the ms header for self calibration
01562         painc = Parallactic angle increment in degrees after which a new convolution function is computed (default=360.0deg)
01563         cfcache = The disk cache directory for convolution functions
01564         pblimit = The fraction of the peak of the PB to which the PB corrections are applied (default=0.1)
01565         dopbcorr = If true, correct for PB in the major cycles as well
01566         applyoffsets = If true, apply antenna pointing offsets from the pointing table given by epjtablename 
01567         epjtablename = Table containing antenna pointing offsets
01568         """
01569 
01570         if(spw==''):
01571             spw='*'
01572         if(field==''):
01573             field='*'
01574         spwids=ms.msseltoindex(vis=msname, spw=spw)['spw']
01575         ###num of cpu per node
01576         numcpu=numcpuperhost
01577         time1=time.time()
01578         self.setupcommonparams(spw=spw, field=field, phasecenter=phasecenter, 
01579                                stokes=stokes, ftmachine=ftmachine, wprojplanes=wprojplanes, 
01580                                facets=facets, imsize=imsize, pixsize=pixsize, weight=weight, 
01581                                robust=robust, npixels=npixels, gain=gain, uvtaper=uvtaper,
01582                                outertaper=outertaper, timerange=timerange, uvrange=uvrange, 
01583                                baselines=baselines, scan=scan, observation=observation, 
01584                                visinmem=visinmem, pbcorr=pbcorr, minpb=minpb, 
01585                                numthreads=numthreads, 
01586                                cyclefactor=cyclefactor)
01587              
01588         self.setupcluster(hostnames,numcpuperhost, 0)
01589         numcpu=self.numcpu
01590         ####
01591         #the default working directory is somewhere 
01592         owd=os.getcwd()
01593         self.c.pgc('import os')
01594         self.c.pgc('os.chdir("'+owd+'")')
01595         #####################
01596         model=imagename+'.model' 
01597         if(not contclean or (not os.path.exists(model))):
01598             shutil.rmtree(model, True)
01599             failedmods=glob.glob(imagename+'*.model')
01600             for someim in failedmods:
01601                 shutil.rmtree(someim, True)
01602             shutil.rmtree(imagename+'.image', True)
01603             ##create the cube
01604             im.selectvis(vis=msname, spw=spw, field=field, writeaccess=False)
01605             im.defineimage(nx=imsize[0], ny=imsize[1], cellx=pixsize[0], celly=pixsize[1], 
01606                            phasecenter=phasecenter, mode=mode, spw=spwids.tolist(), nchan=nchan, step=step, start=start, restfreq=restfreq)
01607             im.setoptions(imagetilevol=imagetilevol) 
01608             #print 'making model image (', model, ') ...'
01609             im.make(model)
01610             print 'model image (', model, ') made'
01611             im.done()
01612         #print 'LOCKS ', tb.listlocks()
01613         
01614         ia.open(model)
01615         elshape=ia.shape()
01616         csys=ia.coordsys()
01617         fstart=csys.toworld([0,0,0,0],'n')['numeric'][3]
01618         fstep=csys.toworld([0,0,0,1],'n')['numeric'][3]-fstart
01619         fend=fstep*(nchan-1)+fstart
01620         ia.done()
01621         ###handle mask
01622         if((maskimage != '')):
01623             if(os.path.exists(maskimage)):
01624                 ia.open(maskimage)
01625                 maskshape=ia.shape()
01626                 ia.done()
01627                 if(np.any(maskshape != elshape)):
01628                     newmask=maskimage+'_regrid'
01629                     self.regridimage(outimage=newmask, inimage=maskimage, templateimage=model);
01630                     maskimage=newmask
01631             else:
01632                 shutil.copytree(model, maskimage)
01633                 ia.open(maskimage)
01634                 ia.set(0.0)
01635                 ia.done()
01636                 
01637         #print 'LOCKS2 ', tb.listlocks()
01638         imepoch=csys.epoch()
01639         imobservatory=csys.telescope()
01640         shutil.rmtree(imagename+'.image', True)
01641         shutil.rmtree(imagename+'.residual', True)
01642         ## don't need to copy anymore with concat
01643         #shutil.copytree(model, imagename+'.image')
01644         #shutil.copytree(model, imagename+'.residual')
01645 
01646         out=range(numcpu)  
01647         self.c.pgc('from  parallel.parallel_cont import *')
01648         ###set some common parameters
01649         self.c.pgc('a.imagetilevol='+str(imagetilevol))
01650         self.c.pgc('a.visInMem='+str(visinmem))
01651         self.c.pgc('a.painc='+str(painc))
01652         self.c.pgc('a.cfcache='+'"'+str(cfcache)+'"')
01653         self.c.pgc('a.pblimit='+str(pblimit));
01654         self.c.pgc('a.dopbcorr='+str(dopbcorr));
01655         self.c.pgc('a.applyoffsets='+str(applyoffsets));
01656         self.c.pgc('a.epjtablename='+'"'+str(epjtablename)+'"');
01657 
01658         tb.clearlocks()
01659         #print 'LOCKS3', tb.listlocks()
01660         chancounter=0
01661         #####
01662         if(contclean):
01663             imagecont.getallchanmodel(imagename , chanchunk)
01664         #####
01665         timemake=time.time()
01666         print 'time to get make cubes', timemake - time1 
01667         nchanchunk=nchan/chanchunk if (nchan%chanchunk) ==0 else nchan/chanchunk+1
01668         ###spw and channel selection
01669         #spwsel,startsel, nchansel=imagecont.findchanselLSRK(msname=msname, spw=spwids, 
01670         #                                             field=field, 
01671         #                                             numpartition=nchanchunk, 
01672         #                                             beginfreq=fstart, endfreq=fend, chanwidth=fstep)
01673         #print 'time to calc selection ', time.time()-timemake
01674         #print 'spwsel', spwsel, 'startsel', startsel,'nchansel', nchansel
01675         ##print 'startsel', startsel
01676         ##print  'nchansel', nchansel
01677         imnam='"%s"'%(imagename)
01678         donegetchan=np.array(range(nchanchunk),dtype=bool)
01679         doneputchan=np.array(range(nchanchunk),dtype=bool)
01680         readyputchan=np.array(range(nchanchunk), dtype=bool)
01681         cpudoing=np.array(range(nchanchunk), dtype=int)
01682         donegetchan.setfield(False,bool)
01683         doneputchan.setfield(False,bool)
01684         readyputchan.setfield(False, bool)
01685         chanind=np.array(range(numcpu), dtype=int)
01686         def gen_command(ccounter):
01687             startfreq=str(fstart+ccounter*chanchunk*fstep)+'Hz'
01688             widthfreq=str(fstep)+'Hz'
01689             imnchan=chanchunk
01690             if((ccounter == (nchanchunk-1)) and ((nchan%chanchunk) != 0)):
01691                 imnchan=nchan%chanchunk
01692                 
01693             return 'a.imagechan(msname='+'"'+msname+'", start='+str(startsel[ccounter])+', numchan='+str(nchansel[ccounter])+', field="'+str(field)+'", spw='+str(spwsel[ccounter])+', imroot='+imnam+',imchan='+str(ccounter)+',niter='+str(niter)+',alg="'+alg+'", scales='+str(scales)+', majcycle='+str(majorcycles)+', thr="'+str(threshold)+'", fstart="'+startfreq+'", width="'+widthfreq+'", chanchunk='+str(imnchan)+', mask="'+maskimage+'")'
01694         def gen_command2(ccounter):
01695             imnchan=chanchunk
01696             startchan=ccounter*chanchunk
01697             if((ccounter == (nchanchunk-1)) and ((nchan%chanchunk) != 0)):
01698                 imnchan=nchan%chanchunk
01699             return 'retval=a.imagechan_selfselect(msname='+'"'+msname+'", field="'+str(field)+'", spwids='+str(spwids.tolist())+', imroot='+imnam+',imchan='+str(ccounter)+',niter='+str(niter)+',alg="'+alg+'", scales='+str(scales)+', majcycle='+str(majorcycles)+', thr="'+str(threshold)+'", chanchunk='+str(imnchan)+', mask="'+maskimage+'", startchan='+str(startchan)+')'
01700         #while(chancounter < nchanchunk):
01701         chanind.setfield(-1, int)
01702         for k in range(numcpu):
01703             if(chancounter < nchanchunk):
01704                 #if((len(nchansel[chancounter])==0) or (len(spwsel[chancounter])==0) or  (len(startsel[chancounter])==0)):
01705                     ###no need to process this channel
01706                 #    doneputchan[chancounter]=True
01707                 #else:
01708 ##need to retab this 
01709                 runcomm=gen_command2(chancounter)
01710                 print 'command is ', runcomm
01711                 out[k]=self.c.odo(runcomm,k)
01712 ##############
01713                 chanind[k]=chancounter
01714                 chancounter=chancounter+1
01715                 time.sleep(1) ###just to avoid an i/o bottle neck
01716 ###############
01717         retval={}
01718         retval['converged']=True
01719         retval['iterations']=0
01720         retval['maxresidual']=0.0
01721         #print 'numcpuused', chancounter, nchanchunk
01722         ##reset numcpu in case less than available is used
01723         numcpu=copy.deepcopy(chancounter)
01724         over=False
01725         #while((chancounter < nchanchunk) and (not over)):
01726         while(not over):
01727             #over=False
01728             #while(not over):
01729             #############loop waiting for a chunk of work
01730                 time.sleep(1)
01731                 overone=True
01732                 for k in range(numcpu):
01733                     overone=(overone and ((type(out[k])==int) or (self.c.check_job(out[k],False))))
01734                     #print k,  'checjob' , ((type(out[k])==int) or (self.c.check_job(out[k],False))), 'chanind', chanind[k], 'chancounter', chancounter, 'readypu', readyputchan[chanind[k]]
01735                     if((chanind[k] > -1) and (not readyputchan[chanind[k]]) and ((type(out[k])==int) or self.c.check_job(out[k],False)) ):
01736                         readyputchan[chanind[k]]=True 
01737                         if(type(out[k]) !=int):
01738                             ##done for this channel
01739                             retvals=self.c.pull('retval', k)
01740                             #print "RETVALS", retvals, 'retval', retval
01741                             retval['iterations']=max(retval['iterations'], retvals[k]['iterations'])
01742                             retval['maxresidual']=max(retval['maxresidual'], retvals[k]['maxresidual'])
01743                             retval['converged']= retval['converged'] and retvals[k]['converged']
01744                         if(chancounter < nchanchunk):
01745                             #if((len(nchansel[chancounter])==0) or (len(spwsel[chancounter])==0) or  (len(startsel[chancounter])==0)):
01746                                 ###no need to process this channel
01747                              #   doneputchan[chancounter]=True
01748                             #else:
01749 ##########to be retabbed
01750                             runcomm=gen_command2(chancounter)
01751                             print 'command is ', runcomm
01752                             print 'processor ', k
01753                             out[k]=self.c.odo(runcomm,k)
01754                             time.sleep(1)
01755 ###################
01756                             chanind[k]=chancounter
01757                             chancounter+=1
01758                         overone=(overone and ((type(out[k])==int) or self.c.check_job(out[k],False)))
01759                         #print 'over', over, 'chancounter', chancounter, 'chanind', chanind
01760                 over=(overone) and (chancounter >= nchanchunk)               
01761                ############
01762         timebegrem=time.time()
01763         print 'Time to image is ', (timebegrem-time1)/60.0, 'mins'
01764         chans=(np.array(range(nchanchunk))*chanchunk).tolist()
01765         imnams=[imagename]*nchanchunk
01766 #        ia.open(imnams[0]+'0.image')
01767 #        rb=ia.restoringbeam()
01768 #        ia.done()
01769 #        ia.open(imagename+'.image')
01770 #        ia.setrestoringbeam(beam=rb)
01771 #        ia.done()
01772         #imagecont.putchanimage2(model , [imnams[k]+str(k)+'.model' for k in range(nchanchunk)], chans, doneputchan.tolist(), True)
01773         #imagecont.putchanimage2(imagename+'.residual' ,[imnams[k]+str(k)+'.residual' for k in range(nchanchunk)] , chans, doneputchan.tolist(), True)
01774         
01775         #imagecont.putchanimage2(imagename+'.image' , [imnams[k]+str(k)+'.image' for k in range(nchanchunk)], chans, doneputchan.tolist(), True)
01776         #self.concatimages(model,  [imnams[k]+str(k)+'.model' for k in range(nchanchunk)], csys)
01777         imagecont.concatimages(model,  [imnams[k]+str(k)+'.model' for k in range(nchanchunk)], csys)
01778         #self.concatimages(imagename+'.residual' ,[imnams[k]+str(k)+'.residual' for k in range(nchanchunk)], csys)
01779         imagecont.concatimages(imagename+'.residual' ,[imnams[k]+str(k)+'.residual' for k in range(nchanchunk)], csys)
01780         #self.concatimages(imagename+'.image' , [imnams[k]+str(k)+'.image' for k in range(nchanchunk)], csys)
01781         imagecont.concatimages(imagename+'.image' , [imnams[k]+str(k)+'.image' for k in range(nchanchunk)], csys)
01782         #self.concatimages(imagename+'.psf' , [imnams[k]+str(k)+'.psf' for k in range(nchanchunk)], csys)
01783         imagecont.concatimages(imagename+'.psf' , [imnams[k]+str(k)+'.psf' for k in range(nchanchunk)], csys)
01784         if(self.ftmachine=='mosaic'):
01785             #self.concatimages(imagename+'.flux' , [imnams[k]+str(k)+'.flux' for k in range(nchanchunk)], csys)
01786             imagecont.concatimages(imagename+'.flux' , [imnams[k]+str(k)+'.flux' for k in range(nchanchunk)], csys)
01787             #self.concatimages(imagename+'.flux.pbcoverage' , [imnams[k]+str(k)+'.flux.pbcoverage' for k in range(nchanchunk)], csys)
01788             imagecont.concatimages(imagename+'.flux.pbcoverage' , [imnams[k]+str(k)+'.flux.pbcoverage' for k in range(nchanchunk)], csys)
01789         time2=time.time()
01790         print 'Time to concat/cleanup', (time2- timebegrem)/60.0, 'mins'
01791         if(savemodel):
01792             myim=casac.imager()
01793             myim.selectvis(vis=msname, spw=spw, field=field)
01794             myim.defineimage()
01795             myim.setoptions(ftmachine=ftmachine,wprojplanes=wprojplanes)
01796             myim.ft(model)
01797             myim.done()
01798 
01799         time2=time.time()
01800         print 'Time to image after cleaning is ', (time2-time1)/60.0, 'mins'
01801         return retval
01802         #self.c.stop_cluster()
01803 
01804 ##############################
01805     def concatimages(self, outputimage='', imagelist='', csys=None, removeinfile=True):
01806         ncpu=self.numcpu
01807         if(type(imagelist) != list):
01808             imagelist=[imagelist]
01809         imagelist=np.array(imagelist)
01810         if(len(imagelist)==1):
01811             ib=ia.fromimage(ourfile=outputimage, infile=imagelist[0], overwrite=True)
01812         else:
01813             nim=len(imagelist)
01814             self.c.pgc('from  parallel.parallel_cont import *')
01815             if((nim/ncpu) > 1):
01816                 intermimages=['intermimages_'+str(k) for k in range(ncpu)]
01817                 subpart=np.array([nim/ncpu for k in range(ncpu)])
01818                 subpart[range(nim%ncpu)] +=1
01819                 infiles=[[]]*ncpu
01820                 out=range(ncpu)
01821                 #infiles[0]=imagelist(range(subpart[0]))
01822                 for k in  range(ncpu):
01823                     infiles[k]=imagelist[range(np.sum(subpart[0:k]), np.sum(subpart[0:(k+1)]))].tolist() 
01824                     commcon='imagecont.concatimages(cubeimage="'+intermimages[k]+'", inim='+str(infiles[k])+', removeinfile='+str(removeinfile)+')'
01825                     print 'concat command', commcon
01826                     out[k]=self.c.odo(commcon,k)
01827                 over=False
01828                 while(not over):
01829                     time.sleep(1)
01830                     overone=True
01831                     for k in range(ncpu):
01832                         overone=(overone and ((type(out[k])==int) or (self.c.check_job(out[k],False))))
01833                     over=overone
01834             else:
01835                 intermimages=[imagelist[k] for k in range(len(imagelist))]
01836             imagecont.concatimages(outputimage , intermimages, csys, removeinfile)
01837 #################
01838     def pcube_expt(self, msname=None, imagename='elimage', imsize=[1000, 1000], 
01839               pixsize=['1arcsec', '1arcsec'], phasecenter='', 
01840               field='', spw='*', ftmachine='ft', wprojplanes=128, facets=1, 
01841               hostnames='', 
01842               numcpuperhost=1, majorcycles=1, niter=1000, gain=0.1, threshold='0.0mJy', alg='clark', scales=[0],
01843               mode='channel', start=0, nchan=1, step=1, stokes='I', restfreq='', weight='natural', 
01844               robust=0.0, npixels=0, 
01845               imagetilevol=100000,
01846               contclean=False, chanchunk=1, visinmem=False, maskimage='' , pbcorr=False, numthreads=-1,
01847               painc=360., pblimit=0.1, dopbcorr=True, applyoffsets=False, cfcache='cfcache.dir',
01848               epjtablename=''): 
01849 
01850         """
01851         msname= measurementset
01852         imagename = image
01853         imsize = list of 2 numbers  [nx,ny] defining image size in x and y
01854         pixsize = list of 2 quantities   ['sizex', 'sizey'] defining the pixel size e.g  ['1arcsec', '1arcsec']
01855         phasecenter = an integer or a direction string   integer is fieldindex or direction e.g  'J2000 19h30m00 -30d00m00'
01856         field = field selection string ...msselection style
01857         spw = spw selection string ...msselection style
01858         ftmachine= the ftmachine to use ...'ft', 'wproject' etc
01859         wprojplanes is an interger that is valid only of ftmachine is 'wproject', 
01860         facets= integer do split image facet, 
01861         hostnames= list of strings ..empty string mean localhost
01862         numcpuperhos = integer ...number of processes to launch on each host
01863         majorcycles= integer number of CS major cycles to do, 
01864         niter= integer ...total number of clean iteration 
01865         threshold=quantity string ...residual peak at which to stop deconvolving
01866         alg= string  possibilities are 'clark', 'hogbom', 'multiscale' and their 'mf'
01867         scales= list of scales in pixel for multiscale clean e.g [0, 3, 10]
01868         mode= channel definition, can be 'channel', 'frequency', 'velocity'
01869         start = first channel in the definition spec of mode, can be int, freq or vel quantity
01870         step = channel width specified in the definition of mode
01871         restfreq= what 'rest frequency' to use to calculate velocity from frequency
01872                                 empty string '' implies use the first restfreq in SOURCE of ms 
01873         weight= type of weight to apply
01874         contclean = boolean ...if False the imagename.model is deleted if its on 
01875         disk otherwise clean will continue from previous run
01876         chanchunk = number of channel to process at a go per process...careful not to 
01877        go above total memory available
01878        visinmem = load visibility in memory for major cycles...make sure totalmemory  available to all processes is more than the MS size
01879         painc = Parallactic angle increment in degrees after which a new convolution function is computed (default=360.0deg)
01880         cfcache = The disk cache directory for convolution functions
01881         pblimit = The fraction of the peak of the PB to which the PB corrections are applied (default=0.1)
01882         dopbcorr = If true, correct for PB in the major cycles as well
01883         applyoffsets = If true, apply antenna pointing offsets from the pointing table given by epjtablename 
01884         epjtablename = Table containing antenna pointing offsets
01885         """
01886 
01887         if(spw==''):
01888             spw='*'
01889         if(field==''):
01890             field='*'
01891         spwids=ms.msseltoindex(vis=msname, spw=spw)['spw']
01892         ###num of cpu per node
01893         numcpu=numcpuperhost
01894         time1=time.time()
01895         self.spw=spw
01896         self.field=field
01897         self.phasecenter=phasecenter
01898         self.ftmachine=ftmachine
01899         self.wprojplanes=wprojplanes
01900         self.facets=facets
01901         self.imsize=imsize
01902         self.cell=pixsize
01903         self.weight=weight
01904         self.robust=robust
01905         self.npixels=npixels
01906         self.visinmem=visinmem
01907         self.numthreads=numthreads
01908         self.stokes=stokes
01909         self.gain=gain
01910         self.pbcorr=pbcorr
01911         self.setupcluster(hostnames,numcpuperhost, 0)
01912         numcpu=self.numcpu
01913         ####
01914         #the default working directory is somewhere 
01915         owd=os.getcwd()
01916         hostname=os.getenv('HOSTNAME')
01917         ###start one engine locally to cleanup when possible
01918         self.c.start_engine(hostname,1,owd)
01919         self.c.pgc('import os')
01920         self.c.pgc('os.chdir("'+owd+'")')
01921         #####################
01922         model=imagename+'.model' 
01923         if(not contclean or (not os.path.exists(model))):
01924             shutil.rmtree(model, True)
01925             shutil.rmtree(imagename+'.image', True)
01926             ##create the cube
01927             im.selectvis(vis=msname, spw=spw, field=field, writeaccess=False)
01928             im.defineimage(nx=imsize[0], ny=imsize[1], cellx=pixsize[0], celly=pixsize[1], 
01929                            phasecenter=phasecenter, mode=mode, spw=spwids.tolist(), nchan=nchan, step=step, start=start, restfreq=restfreq)
01930             im.setoptions(imagetilevol=imagetilevol) 
01931             #print 'making model image (', model, ') ...'
01932             im.make(model)
01933             print 'model image (', model, ') made'
01934             im.done()
01935         #print 'LOCKS ', tb.listlocks()
01936         
01937         ia.open(model)
01938         elshape=ia.shape()
01939         csys=ia.coordsys()
01940         fstart=csys.toworld([0,0,0,0],'n')['numeric'][3]
01941         fstep=csys.toworld([0,0,0,1],'n')['numeric'][3]-fstart
01942         fend=fstep*(nchan-1)+fstart
01943         ia.done()
01944         ###handle mask
01945         if((maskimage != '') and (os.path.exists(maskimage))):
01946             ia.open(maskimage)
01947             maskshape=ia.shape()
01948             ia.done()
01949             if(np.any(maskshape != elshape)):
01950                 newmask=maskimage+'_regrid'
01951                 self.regridimage(outimage=newmask, inimage=maskimage, templateimage=model);
01952                 maskimage=newmask
01953         #print 'LOCKS2 ', tb.listlocks()
01954         imepoch=csys.epoch()
01955         imobservatory=csys.telescope()
01956         shutil.rmtree(imagename+'.image', True)
01957         shutil.rmtree(imagename+'.residual', True)
01958         shutil.copytree(model, imagename+'.image')
01959         shutil.copytree(model, imagename+'.residual')
01960 
01961         out=range(numcpu)
01962         outhelper=0
01963         self.c.pgc('from  parallel.parallel_cont import *')
01964         self.c.odo('a=imagecont()', numcpu)
01965         ###set some common parameters
01966         self.c.pgc('a.imagetilevol='+str(imagetilevol))
01967         self.c.pgc('a.visInMem='+str(visinmem))
01968         self.c.pgc('a.painc='+str(painc))
01969         self.c.pgc('a.cfcache='+'"'+str(cfcache)+'"')
01970         self.c.pgc('a.pblimit='+str(pblimit));
01971         self.c.pgc('a.dopbcorr='+str(dopbcorr));
01972         self.c.pgc('a.applyoffsets='+str(applyoffsets));
01973         self.c.pgc('a.epjtablename='+'"'+str(epjtablename)+'"');
01974 
01975         tb.clearlocks()
01976         #print 'LOCKS3', tb.listlocks()
01977         chancounter=0
01978         #####
01979         if(contclean):
01980             imagecont.getallchanmodel(imagename , chanchunk)
01981         #####
01982         timemake=time.time()
01983         print 'time to get make cubes', timemake - time1 
01984         nchanchunk=nchan/chanchunk if (nchan%chanchunk) ==0 else nchan/chanchunk+1
01985         ###spw and channel selection
01986         #spwsel,startsel, nchansel=imagecont.findchanselLSRK(msname=msname, spw=spwids, 
01987         #                                             field=field, 
01988         #                                             numpartition=nchanchunk, 
01989         #                                             beginfreq=fstart, endfreq=fend, chanwidth=fstep)
01990         #print 'time to calc selection ', time.time()-timemake
01991         #print 'spwsel', spwsel, 'startsel', startsel,'nchansel', nchansel
01992         ##print 'startsel', startsel
01993         ##print  'nchansel', nchansel
01994         imnam='"%s"'%(imagename)
01995         donegetchan=np.array(range(nchanchunk),dtype=bool)
01996         doneputchan=np.array(range(nchanchunk),dtype=bool)
01997         readyputchan=np.array(range(nchanchunk), dtype=bool)
01998         cpudoing=np.array(range(nchanchunk), dtype=int)
01999         donegetchan.setfield(False,bool)
02000         doneputchan.setfield(False,bool)
02001         readyputchan.setfield(False, bool)
02002         chanind=np.array(range(numcpu), dtype=int)
02003         cleanup='a.cleanupimages(readyputchan=readyputchan, imagename='+imnam+',nchanchunk='+str(nchanchunk)+',chanchunk='+str(chanchunk)+')'
02004         def gen_command(ccounter):
02005             startfreq=str(fstart+ccounter*chanchunk*fstep)+'Hz'
02006             widthfreq=str(fstep)+'Hz'
02007             imnchan=chanchunk
02008             if((ccounter == (nchanchunk-1)) and ((nchan%chanchunk) != 0)):
02009                 imnchan=nchan%chanchunk
02010                 
02011             return 'a.imagechan(msname='+'"'+msname+'", start='+str(startsel[ccounter])+', numchan='+str(nchansel[ccounter])+', field="'+str(field)+'", spw='+str(spwsel[ccounter])+', imroot='+imnam+',imchan='+str(ccounter)+',niter='+str(niter)+',alg="'+alg+'", scales='+str(scales)+', majcycle='+str(majorcycles)+', thr="'+str(threshold)+'", fstart="'+startfreq+'", width="'+widthfreq+'", chanchunk='+str(imnchan)+', mask="'+maskimage+'")'
02012         def gen_command2(ccounter):
02013             imnchan=chanchunk
02014             if((ccounter == (nchanchunk-1)) and ((nchan%chanchunk) != 0)):
02015                 imnchan=nchan%chanchunk
02016             return 'a.imagechan_selfselect(msname='+'"'+msname+'", field="'+str(field)+'", spwids='+str(spwids.tolist())+', imroot='+imnam+',imchan='+str(ccounter)+',niter='+str(niter)+',alg="'+alg+'", scales='+str(scales)+', majcycle='+str(majorcycles)+', thr="'+str(threshold)+'", chanchunk='+str(imnchan)+', mask="'+maskimage+'")'
02017         #while(chancounter < nchanchunk):
02018         chanind.setfield(-1, int)
02019         for k in range(numcpu):
02020             if(chancounter < nchanchunk):
02021                 #if((len(nchansel[chancounter])==0) or (len(spwsel[chancounter])==0) or  (len(startsel[chancounter])==0)):
02022                     ###no need to process this channel
02023                 #    doneputchan[chancounter]=True
02024                 #else:
02025 ##need to retab this 
02026                 runcomm=gen_command2(chancounter)
02027                 print 'command is ', runcomm
02028                 out[k]=self.c.odo(runcomm,k)
02029 ##############
02030                 chanind[k]=chancounter
02031                 chancounter=chancounter+1
02032 ###############
02033         #print 'numcpuused', chancounter, nchanchunk
02034         ##reset numcpu in case less than available is used
02035         numcpu=copy.deepcopy(chancounter)
02036         over=False
02037         helpover=True
02038         #while((chancounter < nchanchunk) and (not over)):
02039         while(not over):
02040             #over=False
02041             #while(not over):
02042             #############loop waiting for a chunk of work
02043                 time.sleep(1)
02044                 overone=True
02045                 for k in range(numcpu):
02046                     overone=(overone and ((type(out[k])==int) or (self.c.check_job(out[k],False))))
02047                     #print k,  'checjob' , ((type(out[k])==int) or (self.c.check_job(out[k],False))), 'chanind', chanind[k], 'chancounter', chancounter, 'readypu', readyputchan[chanind[k]]
02048                     if((chanind[k] > -1) and (not readyputchan[chanind[k]]) and ((type(out[k])==int) or self.c.check_job(out[k],False)) ):
02049                         readyputchan[chanind[k]]=True      
02050                         if(chancounter < nchanchunk):
02051                             #if((len(nchansel[chancounter])==0) or (len(spwsel[chancounter])==0) or  (len(startsel[chancounter])==0)):
02052                                 ###no need to process this channel
02053                              #   doneputchan[chancounter]=True
02054                             #else:
02055 ##########to be retabbed
02056                             runcomm=gen_command2(chancounter)
02057                             print 'command is ', runcomm
02058                             print 'processor ', k
02059                             out[k]=self.c.odo(runcomm,k)
02060 ###################
02061                             chanind[k]=chancounter
02062                             chancounter+=1
02063                         overone=(overone and ((type(out[k])==int) or self.c.check_job(out[k],False)))
02064                         if(helpover and ((type(outhelper)==int) or self.c.check_job(outhelper, False))):
02065                             self.c.push(readyputchan=readyputchan, targets=numcpu)
02066                             outhelper=self.c.odo(cleanup, numcpu)
02067                         #print 'over', over, 'chancounter', chancounter, 'chanind', chanind
02068                         
02069                 over=(overone) and (chancounter >= nchanchunk)               
02070                ############
02071         timebegrem=time.time()
02072         print 'Time to image is ', (timebegrem-time1)/60.0, 'mins'
02073         while(not self.c.check_job(outhelper, False)):
02074             time.sleep(1)
02075         chans=(np.array(range(nchanchunk))*chanchunk).tolist()
02076         imnams=[imagename]*nchanchunk
02077         imagecont.putchanimage2(model , [imnams[k]+str(k)+'.model' for k in range(nchanchunk)], chans, doneputchan.tolist(), True)
02078         imagecont.putchanimage2(imagename+'.residual' ,[imnams[k]+str(k)+'.residual' for k in range(nchanchunk)] , chans, doneputchan.tolist(), True)
02079         imagecont.putchanimage2(imagename+'.image' , [imnams[k]+str(k)+'.image' for k in range(nchanchunk)], chans, doneputchan.tolist(), True)
02080         time2=time.time()
02081         print 'Time to concat/cleanup', (time2- timebegrem)/60.0, 'mins'
02082         
02083 
02084         time2=time.time()
02085         print 'Time to image after cleaning is ', (time2-time1)/60.0, 'mins'
02086         #self.c.stop_engine(numcpu)
02087 
02088 ##############################
02089 
02090     def pcontmultims(self, msnames=[], workdirs=[], imagename=None, imsize=[1000, 1000], 
02091                      pixsize=['1arcsec', '1arcsec'], phasecenter='', 
02092                      field='', spw='*', freqrange=['', ''],  stokes='I', ftmachine='ft', wprojplanes=128, facets=1, 
02093                      hostnames='', 
02094                      numcpuperhost=1, majorcycles=1, niter=1000, gain=0.1, threshold='0.0mJy', alg='clark', scales=[], weight='natural', robust=0.0, npixels=0, pbcorr=False, 
02095                      contclean=False, visinmem=False, maskimage='lala.mask',  numthreads=1,
02096                      painc=360., pblimit=0.1, dopbcorr=True, applyoffsets=False, cfcache='cfcache.dir',
02097                      epjtablename=''):
02098         """
02099         msnames=  list containing measurementset names
02100         imagename = image
02101         imsize = list of 2 numbers  [nx,ny] defining image size in x and y
02102         pixsize = list of 2 quantities   ['sizex', 'sizey'] defining the pixel size e.g  ['1arcsec', '1arcsec']
02103         phasecenter = an integer or a direction string   integer is fieldindex or direction e.g  'J2000 19h30m00 -30d00m00'
02104         field = field selection string ...msselection style
02105         spw = spw selection string ...msselection style
02106         freqrange= continuum image frequency bound e.g ['1GHz', '1.23GHz'] .
02107         stokes= string e.g 'I' , 'IV'
02108         ftmachine= the ftmachine to use ...'ft', 'wproject' etc
02109         wprojplanes is an interger that is valid only of ftmachine is 'wproject', 
02110         facets= integer do split image facet, 
02111         hostnames= list of strings ..empty string mean localhost
02112         numcpuperhos = integer ...number of processes to launch on each host
02113         majorcycles= integer number of CS major cycles to do, 
02114         niter= integer ...total number of clean iteration 
02115         threshold=quantity ...residual peak at which to stop deconvolving
02116         alg= string  possibilities are 'clark', 'hogbom', 'msclean'
02117         weight= string  possibilities 'natural', 'briggs', 'radial'
02118         robust= float  robust factor for briggs
02119         scales = scales to use when using alg='msclean'
02120         contclean = boolean ...if False the imagename.model is deleted if its on 
02121         disk otherwise clean will continue from previous run
02122         visinmem = load visibility in memory for major cycles...make sure totalmemory  available to all processes is more than the MS size
02123         maskimage an image on disk to limit clean search
02124         painc = Parallactic angle increment in degrees after which a new convolution function is computed (default=360.0deg)
02125         cfcache = The disk cache directory for convolution functions
02126         pblimit = The fraction of the peak of the PB to which the PB corrections are applied (default=0.1)
02127         dopbcorr = If true, correct for PB in the major cycles as well
02128         applyoffsets = If true, apply antenna pointing offsets from the pointing table given by epjtablename 
02129         epjtablename = Table containing antenna pointing offsets
02130         """
02131         time1=time.time()
02132         if len(msnames)==0:
02133             return    
02134         niterpercycle=niter/majorcycles
02135         if(niterpercycle == 0):
02136             niterpercycle=niter
02137             majorcycles=1
02138         num_ext_procs=0
02139         self.spw=spw
02140         self.field=field
02141         self.phasecenter=phasecenter
02142         self.stokes=stokes
02143         self.ftmachine=ftmachine
02144         self.wprojplanes=wprojplanes
02145         self.facets=facets
02146         self.imsize=imsize
02147         self.cell=pixsize
02148         self.weight=weight
02149         self.robust=robust
02150         self.npixels=npixels
02151         self.visinmem=visinmem
02152         self.gain=gain
02153         self.numthreads=numthreads
02154         self.pbcorr=pbcorr
02155         self.setupcluster(hostnames,numcpuperhost, num_ext_procs)
02156         model=imagename+'.model' if (len(imagename) != 0) else 'elmodel'
02157         if(not contclean):
02158             print "Removing ", model, 'and', imagename+'.image'
02159             shutil.rmtree(model, True)
02160             shutil.rmtree(imagename+'.image', True)
02161             shutil.rmtree(imagename+'.residual', True)
02162         shutil.rmtree('tempmodel', True)
02163         if(len(freqrange) < 2):
02164             freqrange=['','']
02165         if(freqrange[0]==''):
02166             freqrange[0]=0.0
02167         if(freqrange[1]==''):
02168             freqrange[1]=1e24
02169         if(type(freqrange[0])==str):
02170             freqrange[0]=qa.convert(qa.quantity(freqrange[0], 'Hz'), 'Hz')['value']
02171         if(type(freqrange[1])==str):
02172             freqrange[1]=qa.convert(qa.quantity(freqrange[1], 'Hz'), 'Hz')['value']
02173         numms=len(msnames) if (type(msnames)==list) else 1
02174         if(type(msnames) == str):
02175             msnames=[msnames]
02176         freqrange=self.findfreqranges(msnames=msnames, spw=spw, freqmin=freqrange[0], freqmax=freqrange[1])
02177         ## the above will return the freqmin, freqmax in the data if freqrange is default
02178         freq=(freqrange[0]+freqrange[1])/2.0
02179         band=abs(freqrange[1]-freqrange[0])
02180         self.makeconttempimages(imagename, self.numcpu, contclean)
02181         def gen_comm(msnames, field, freq, band, imname):
02182             spwsel=[]
02183             startsel=[]
02184             nchansel=[]
02185             for k in range(len(msnames)): 
02186                 msname=msnames[k]
02187                 spwsel.append(self.msinfo[msname]['spwids'].tolist())
02188                 startsel.append(self.msinfo[msname]['startsel'])
02189                 nchansel.append(self.msinfo[msname]['nchansel'])
02190             freqsel='"%fHz"'%freq
02191             bandsel='"%fHz"'%band
02192             return 'a.imagecontmultims(msnames='+str(msnames)+', start='+str(startsel)+', numchan='+str(nchansel)+', field="'+str(field)+'", spw='+str(spwsel)+', freq='+freqsel+', band='+ bandsel+', imname="'+imname+'")'
02193         ###major cycle
02194         if(numms < self.numcpu):
02195             self.numcpu=numms
02196         out=range(self.numcpu)
02197         
02198         for k in range(self.numcpu):
02199             if(not self.engineinfo.has_key(k)):
02200                 self.engineinfo[k]={}
02201             self.engineinfo[k]['msnames']=[]
02202         t_msnames=copy.deepcopy(msnames)
02203         t_msnames.reverse()
02204         counter=0
02205         for k in range(len(msnames)):
02206             self.engineinfo[counter]['msnames'].append(t_msnames.pop())
02207             counter=counter+1
02208             if(counter==self.numcpu):
02209                 counter=0
02210         for maj in range(majorcycles):
02211             myrec=self.msinfo
02212             ##initial bunch of launches
02213             for k in range(self.numcpu):
02214                 runcomm=gen_comm(msnames=self.engineinfo[k]['msnames'], 
02215                                  field=self.field, freq=freq, band=band, 
02216                                  imname=self.engineinfo[k]['imname'])
02217                 #print 'cpu', k,  'command is ', runcomm
02218                 out[k]=self.c.odo(runcomm,k)
02219             over=False
02220             while (not over):
02221                 over=True
02222                 time.sleep(1)
02223                 for k in range(self.numcpu):
02224                     over=(over and self.c.check_job(out[k], False))
02225             residual=imagename+'.residual'
02226             psf=imagename+'.psf'
02227             psfs=range(self.numcpu)
02228             residuals=range(self.numcpu)
02229             restoreds=range(self.numcpu)
02230             for k in range (self.numcpu):
02231                 psfs[k]=self.engineinfo[k]['imname']+'.psf'
02232                 residuals[k]=self.engineinfo[k]['imname']+'.residual'
02233                 restoreds[k]=self.engineinfo[k]['imname']+'.image'
02234             self.averimages(residual, residuals)
02235             if(maj==0):
02236                 self.averimages(psf, psfs)
02237                 if(os.path.exists(maskimage)):
02238                     self.regridimage(outimage='__lala.mask', inimage=maskimage, templateimage=residual);
02239                     shutil.rmtree(maskimage, True)
02240                     shutil.move('__lala.mask', maskimage)
02241                 else:
02242                     self.copyimage(inimage=residual, outimage=maskimage, init=True, initval=1.0);
02243                 if(not contclean or (not os.path.exists(model))):
02244                     self.copyimage(inimage=residual, outimage=model, 
02245                                    init=True, initval=0.0)
02246             if scales==[]:
02247                 scales=[0]
02248             self.incrementaldecon(alg=alg, scales=scales, residual=residual, model=model, niter=niterpercycle, psf=psf, mask=maskimage, thr=threshold, cpuid=0, imholder=self.engineinfo)
02249         #######
02250         restored=imagename+'.image'
02251         self.averimages(restored, restoreds)
02252         ###close major cycle loop
02253         #shutil.rmtree(imagename+'.image', True)
02254         #shutil.move(restored,  imagename+'.image')
02255         time2=time.time()
02256         ###Clean up
02257         for k in range(self.numcpu):
02258             imlist=self.engineinfo[k]['imname']
02259             shutil.rmtree(imlist+'.model', True)
02260             shutil.rmtree(imlist+'.residual', True)
02261             shutil.rmtree(imlist+'.image', True)
02262             shutil.rmtree(imlist+'.psf', True)
02263         print 'Time to image is ', (time2-time1)/60.0, 'mins'
02264         #self.c.stop_cluster()
02265    
02266 
02267     def pcontmultims2(self, msnames=[], workdirs=[], imagename=None, imsize=[1000, 1000], 
02268                      pixsize=['1arcsec', '1arcsec'], phasecenter='', 
02269                      field='', spw='*', freqrange=['', ''],  stokes='I', ftmachine='ft', wprojplanes=128, facets=1, 
02270                      hostnames='', 
02271                      numcpuperhost=1, majorcycles=1, niter=1000, gain=0.1, threshold='0.0mJy', alg='clark', scales=[0], weight='natural', robust=0.0, npixels=0, 
02272                      contclean=False, visinmem=False, maskimage='lala.mask',
02273                      painc=360., pblimit=0.1, dopbcorr=True, applyoffsets=False, cfcache='cfcache.dir',
02274                      epjtablename=''):
02275         """
02276         msnames=  list containing measurementset names
02277         imagename = image
02278         imsize = list of 2 numbers  [nx,ny] defining image size in x and y
02279         pixsize = list of 2 quantities   ['sizex', 'sizey'] defining the pixel size e.g  ['1arcsec', '1arcsec']
02280         phasecenter = an integer or a direction string   integer is fieldindex or direction e.g  'J2000 19h30m00 -30d00m00'
02281         field = field selection string ...msselection style
02282         spw = spw selection string ...msselection style
02283         freqrange= continuum image frequency bound e.g ['1GHz', '1.23GHz'] .
02284         stokes= string e.g 'I' , 'IV'
02285         ftmachine= the ftmachine to use ...'ft', 'wproject' etc
02286         wprojplanes is an interger that is valid only of ftmachine is 'wproject', 
02287         facets= integer do split image facet, 
02288         hostnames= list of strings ..empty string mean localhost
02289         numcpuperhos = integer ...number of processes to launch on each host
02290         majorcycles= integer number of CS major cycles to do, 
02291         niter= integer ...total number of clean iteration 
02292         threshold=quantity ...residual peak at which to stop deconvolving
02293         alg= string  possibilities are 'clark', 'hogbom', 'msclean'
02294         weight= string  possibilities 'natural', 'briggs', 'radial'
02295         robust= float  robust factor for briggs
02296         scales = scales to use when using alg='msclean'
02297         contclean = boolean ...if False the imagename.model is deleted if its on 
02298         disk otherwise clean will continue from previous run
02299         visinmem = load visibility in memory for major cycles...make sure totalmemory  available to all processes is more than the MS size
02300         maskimage an image on disk to limit clean search
02301         painc = Parallactic angle increment in degrees after which a new convolution function is computed (default=360.0deg)
02302         cfcache = The disk cache directory for convolution functions
02303         pblimit = The fraction of the peak of the PB to which the PB corrections are applied (default=0.1)
02304         dopbcorr = If true, correct for PB in the major cycles as well
02305         applyoffsets = If true, apply antenna pointing offsets from the pointing table given by epjtablename 
02306         epjtablename = Table containing antenna pointing offsets
02307         """
02308         time1=time.time()
02309         if len(msnames)==0:
02310             return    
02311         niterpercycle=niter/majorcycles
02312         if(niterpercycle == 0):
02313             niterpercycle=niter
02314             majorcycles=1
02315         num_ext_procs=0
02316         self.spw=spw
02317         self.field=field
02318         self.phasecenter=phasecenter
02319         self.stokes=stokes
02320         self.ftmachine=ftmachine
02321         self.wprojplanes=wprojplanes
02322         self.facets=facets
02323         self.imsize=imsize
02324         self.cell=pixsize
02325         self.weight=weight
02326         self.robust=robust
02327         self.npixels=npixels
02328         self.gain=gain
02329         self.visinmem=visinmem
02330         
02331         self.setupcluster(hostnames,numcpuperhost, num_ext_procs)
02332         model=imagename+'.model' if (len(imagename) != 0) else 'elmodel'
02333         if(not contclean):
02334             print "Removing ", model, 'and', imagename+'.image'
02335             shutil.rmtree(model, True)
02336             shutil.rmtree(imagename+'.image', True)
02337             shutil.rmtree(imagename+'.residual', True)
02338         shutil.rmtree('tempmodel', True)
02339         if(len(freqrange) < 2):
02340             freqrange=['','']
02341         if(freqrange[0]==''):
02342             freqrange[0]=0.0
02343         if(freqrange[1]==''):
02344             freqrange[1]=1e24
02345         if(type(freqrange[0])==str):
02346             freqrange[0]=qa.convert(qa.quantity(freqrange[0], 'Hz'), 'Hz')['value']
02347         if(type(freqrange[1])==str):
02348             freqrange[1]=qa.convert(qa.quantity(freqrange[1], 'Hz'), 'Hz')['value']
02349         numms=len(msnames) if (type(msnames)==list) else 1
02350         if(type(msnames) == str):
02351             msnames=[msnames]
02352         freqrange=self.findfreqranges(msnames=msnames, spw=spw, freqmin=freqrange[0], freqmax=freqrange[1])
02353         ## the above will return the freqmin, freqmax in the data if freqrange is default
02354         freq=(freqrange[0]+freqrange[1])/2.0
02355         band=abs(freqrange[1]-freqrange[0])
02356         self.makeconttempimages(imagename, 0,contclean)
02357         def gen_comm(msname, startsel, nchansel, field, spw, freq, band, imname):
02358             spwsel=str(self.msinfo[msname]['spwids'].tolist())
02359             freqsel='"%fHz"'%freq
02360             bandsel='"%fHz"'%band
02361             return 'a.imagecont(msname='+'"'+msname+'", start='+str(startsel)+', numchan='+str(nchansel)+', field="'+str(field)+'", spw='+spwsel+', freq='+freqsel+', band='+ bandsel+', imname="'+imname+'")'
02362         ###major cycle
02363         out=range(self.numcpu)  
02364         mscpuindex=range(self.numcpu) 
02365         msassigned=odict()
02366         for k in range(len(msnames)):
02367             msassigned[msnames[k]]=-1
02368         for maj in range(majorcycles):
02369             ###set the processed flags off
02370             processing=odict()
02371             processed=odict()
02372             mscounter=0
02373             for k in range(len(msnames)):
02374                 processing[msnames[k]]=False
02375                 processed[msnames[k]]=False
02376             myrec=self.msinfo
02377             ##initial bunch of launches
02378             for k in range(self.numcpu):
02379                 msname=myrec.keys()[k]
02380                 if((msassigned[msname]==k) or (msassigned[msname] <0)):
02381                     mscpuindex[k]=msname
02382                     if(mscounter < len(msnames)):
02383                         runcomm=gen_comm(msname=msname, startsel=myrec[msname]['startsel'], nchansel=myrec[msname]['nchansel'],
02384                                     field=self.field, spw=myrec[msname]['spwsel'], freq=freq, band=band, imname=myrec[msname]['imname'])
02385                         print 'cpu', k,  'command is ', runcomm
02386                         out[k]=self.c.odo(runcomm,k)
02387                         mscounter +=1
02388                         processing[msname]=True
02389                         msassigned[msname]=k
02390             while (not np.alltrue(processed.values())):
02391                 overone=True
02392                 time.sleep(1)
02393                 for k in range(self.numcpu):
02394                     if(processing[mscpuindex[k]] and self.c.check_job(out[k], False)):
02395                         processed[mscpuindex[k]]=True
02396                         processing[mscpuindex[k]]=False
02397                         found=False
02398                         mscounter=0
02399                         while(not found):
02400                             msname=myrec.keys()[mscounter]
02401                             if(not processed[msname] and ((msassigned[msname]==k) or (msassigned[msname] <0))):
02402                                 processing[msname]=True
02403                                 runcomm=gen_comm(msname=msname, startsel=myrec[msname]['startsel'], 
02404                                                  nchansel=myrec[msname]['nchansel'],
02405                                                  field=self.field, spw=myrec[msname]['spwsel'], freq=freq, band=band, 
02406                                                  imname=myrec[msname]['imname'])
02407                                 print 'cpu', k, 'command is ', runcomm
02408                                 out[k]=self.c.odo(runcomm,k)
02409                                 mscpuindex[k]=msname
02410                                 msassigned[msname]=k
02411                                 found=True
02412                             mscounter+=1
02413                             if(mscounter == len(msnames)):
02414                                 found=True
02415 #                    if(processing[mscpuindex[k]] and self.c.check_job(out[k],False) and (not processed[mscpuindex[k]])):
02416 #                        processed[mscpuindex[k]]=True
02417 #                        if(mscounter < len(msnames)):
02418 #                            msname=myrec.keys()[mscounter]
02419 #                            runcomm=gen_comm(msname=msname, startsel=myrec[msname]['startsel'], nchansel=myrec[msname]['nchansel'],
02420 #                                    field=self.field, spw=myrec[msname]['spwsel'], freq=freq, band=band, imname=myrec[msname]['imname'])
02421 #                            print 'cpu', k, 'command is ', runcomm
02422 #                            out[k]=self.c.odo(runcomm,k)
02423 #                            mscpuindex[k]=msname
02424 #                            mscounter +=1
02425 #                            processing[msname]=True
02426             residual=imagename+'.residual'
02427             psf=imagename+'.psf'
02428             psfs=range(len(msnames))
02429             residuals=range(len(msnames))
02430             restoreds=range(len(msnames))
02431             for k in range (len(msnames)):
02432                 psfs[k]=self.msinfo[msnames[k]]['imname']+'.psf'
02433                 residuals[k]=self.msinfo[msnames[k]]['imname']+'.residual'
02434                 restoreds[k]=self.msinfo[msnames[k]]['imname']+'.image'
02435             self.averimages(residual, residuals)
02436             if(maj==0):
02437                 self.averimages(psf, psfs)
02438                 if(os.path.exists(maskimage)):
02439                     self.regridimage(outimage='__lala.mask', inimage=maskimage, templateimage=residual);
02440                     shutil.rmtree(maskimage, True)
02441                     shutil.move('__lala.mask', maskimage)
02442                 else:
02443                     self.copyimage(inimage=residual, outimage=maskimage, init=True, initval=1.0);
02444                 if(not contclean or (not os.path.exists(model))):
02445                     self.copyimage(inimage=residual, outimage=model, 
02446                                    init=True, initval=0.0)
02447             self.incrementaldecon(alg=alg, scales=scales, residual=residual, model=model, niter=niterpercycle, psf=psf, mask=maskimage, thr=threshold, cpuid=0, imholder=self.msinfo)
02448         #######
02449         restored=imagename+'.image'
02450         self.averimages(restored, restoreds)
02451         ###close major cycle loop
02452         #shutil.rmtree(imagename+'.image', True)
02453         #shutil.move(restored,  imagename+'.image')
02454         time2=time.time()
02455         ###Clean up
02456         for k in range(len(msnames)):
02457             imlist=self.msinfo[msnames[k]]['imname']
02458             shutil.rmtree(imlist+'.model', True)
02459             shutil.rmtree(imlist+'.residual', True)
02460             shutil.rmtree(imlist+'.image', True)
02461             shutil.rmtree(imlist+'.psf', True)
02462         print 'Time to image is ', (time2-time1)/60.0, 'mins'
02463         #self.c.stop_cluster()
02464 
02465     def incrementaldecon(self, alg, scales,residual, model, niter, psf,  mask, thr, cpuid, imholder):
02466         ##############
02467             ia.open(residual)
02468             print 'Residual', ia.statistics() 
02469             ia.done()
02470             ########
02471             #incremental clean...get rid of tempmodel
02472             shutil.rmtree('tempmodel', True)
02473             rundecon='a.cleancont(alg="'+str(alg)+'", scales='+str(scales)+',niter='+str(niter)+',psf="'+psf+'", dirty="'+residual+'", model="'+'tempmodel'+'", mask='+'"'+mask+'", thr="'+str(thr)+'")'
02474             print 'Deconvolution command', rundecon
02475             out=self.c.odo(rundecon,cpuid)
02476             over=False
02477             while(not over):
02478                 time.sleep(1)
02479                 over=self.c.check_job(out,False)
02480             ###incremental added to total 
02481             ia.open(model)
02482             ia.calc('"'+model+'" +  "tempmodel"')
02483             ia.done()
02484             ia.open('tempmodel')
02485             #arr=ia.getchunk()
02486             print 'min max of incrmodel', ia.statistics()['min'], ia.statistics()['max']
02487             ia.done()
02488             ia.open(model)
02489             #arr=ia.getchunk()
02490             print 'min max of model', ia.statistics()['min'], ia.statistics()['max']
02491             ia.done()
02492             imkeys=imholder.keys()
02493             for k in range(len(imkeys)):
02494                 imlist=imholder[imkeys[k]]['imname']
02495                 ia.open(imlist+'.model')
02496                 #ia.putchunk(arr)
02497                 ia.insert(infile=model, locate=[0,0,0,0])
02498                 ia.done()
02499         #######
02500     def makeconttempimages(self, imname, numim=0, contclean=False): 
02501         myrec=self.msinfo
02502         msnames=myrec.keys()
02503         if(numim==0):
02504             numim=len(msnames)
02505         char_set=string.ascii_letters
02506         for k in range(numim):
02507             substr='Temp_'+str(k)+'_'+string.join(random.sample(char_set,8), sep='')
02508             substr=string.replace(substr, '/','_')
02509             while (os.path.exists(substr+'.model')):
02510                 substr='Temp_'+str(k)+'_'+string.join(random.sample(char_set,8), sep='')
02511                 substr=string.replace(substr, '/','_')
02512             if(numim==len(msnames)):
02513                 myrec[msnames[k]]['imname']=substr
02514             self.engineinfo[k]={'imname':substr}
02515             model=imname+'.model'
02516             if(contclean and os.path.exists(model)):
02517                 self.copyimage(inimage=model, outimage=substr+'.model', init=False)
02518             else:
02519                 shutil.rmtree(substr+'.model', True)
02520     def findfreqranges(self, msnames=[], spw=['*'], freqmin=0.0,  freqmax=1.e24):
02521         if(type(spw)==str):
02522             spw=[spw]
02523         if (len(spw)==1):
02524             for k in range(1,len(msnames)) :
02525                 spw.append(spw[0])
02526         retfreqmin=1.0e24
02527         retfreqmax=-1.0e24
02528         for k in range(len(msnames)):
02529             self.msinfo[msnames[k]]=odict()
02530             myrec=self.msinfo[msnames[k]]
02531             freqrange=[0,0]
02532             elselec=ms.msseltoindex(vis=msnames[k], spw=spw[k])
02533             channel=elselec['channel']
02534             myrec['spwids']=elselec['spw']
02535             spws, starts, nchans=self.findchanselcont(
02536                 msname=msnames[k], spwids=myrec['spwids'], 
02537                 numpartition=1, beginfreq=freqmin, 
02538                 endfreq=freqmax, freqrange=freqrange)
02539             myrec['spwsel']=spws[0]
02540             myrec['startsel']=starts[0]
02541             myrec['nchansel']=nchans[0]
02542             ##try to respect channel selection
02543             if(len(channel > 0)):
02544                 if(channel.shape[0]==len(spws[0])):
02545                     for elspw in range(len(spws[0])):
02546                         if(myrec['startsel'][elspw] < channel[elspw][1]):
02547                             myrec['startsel'][elspw]=channel[elspw][1]
02548                         if((myrec['nchansel'][elspw] > (channel[elspw][2]- myrec['startsel'][elspw]))):
02549                             myrec['nchansel'][elspw]=(channel[elspw][2]- myrec['startsel'][elspw])+1  
02550             retfreqmax=freqrange[1] if (freqrange[1] > retfreqmax) else retfreqmax
02551             retfreqmin=freqrange[0] if (freqrange[0] <  retfreqmin) else retfreqmin
02552         freqmin=retfreqmin
02553         freqmax=retfreqmax
02554         return [freqmin, freqmax]
02555 
02556     
02557 
02558 
02559 
02560     def setupcluster(self, hostnames='', numcpuperhost=1, num_ext_procs=0, workingdir=''):
02561         """
02562         This sets up the cluster with numcpuperhost on each of the hostnames
02563         working dir by default is  pwd or it can be a vector of 1 or nhostnames
02564         if  1 
02565         num_ext_procs get to be set on master host (i.e this host) and working dir for them is pwd
02566         num_ext_procs is basically to help the master process do stuff asynchronously
02567         """
02568         if(type(self.c) == str):
02569             self.c=cluster()
02570             hostname=os.getenv('HOSTNAME')
02571             wd=os.getcwd()
02572         #pdb.set_trace()
02573             if((workingdir=='') or (workingdir==['']) or (workingdir==[])):
02574                 workingdir=[wd]
02575             if(type(workingdir)==str):
02576                 workingdir=[workingdir]
02577             self.workingdirs=workingdir
02578             if((hostnames==[]) or (hostnames=='')): 
02579                 self.hostnames=[hostname]
02580             else:
02581                 self.hostnames=hostnames
02582             if(len(self.workingdirs)==1):
02583                 for k in range(1, len(hostnames)):
02584                     self.workingdirs.append(self.workingdirs[0])
02585             if(len(hostnames) != len(self.workingdirs)):
02586                 raise ValueError, "Length of hostnames and workingdirs do not match"
02587 
02588             for k in range(len(hostnames)):
02589                 hostname=hostnames[k]
02590                 print 'starting ', numcpuperhost, 'on', hostname, 'with wd', self.workingdirs[k]
02591                 self.c.start_engine(hostname,numcpuperhost,self.workingdirs[k])
02592             self.numcpu=len(hostnames)*numcpuperhost
02593             self.numextraprocs=num_ext_procs
02594             if(num_ext_procs >0):
02595                 hostname=os.getenv('HOSTNAME')
02596                 self.c.start_engine(hostname,num_ext_procs,wd)
02597         else:
02598             ###grab some useful info from the preset cluster
02599             self.numcpu=len(self.c.get_ids())-num_ext_procs
02600             self.hostnames=self.c.get_nodes()
02601             self.c.pgc('import os')
02602             self.c.pgc('wd=os.getcwd()')
02603             wdrec=self.c.pull('wd')
02604             self.workingdirs=['']*len(wdrec)
02605             for k in range(len(wdrec)):
02606                 #print 'WORKING DIR for proc ', k , ' is ' , wdrec[k]
02607                 self.workingdirs[k]=wdrec[k]            
02608         ###do the common stuff to all child
02609         #self.c.pgc('casalog.filter()')
02610         self.c.pgc('from  parallel.parallel_cont import *')
02611         spwlaunch='"'+self.spw+'"' if (type(self.spw)==str) else str(self.spw)
02612         fieldlaunch='"'+self.field+'"' if (type(self.field) == str) else str(self.field)
02613         pslaunch='"'+self.phasecenter+'"' if (type(self.phasecenter) == str) else str(self.phasecenter)
02614         launchcomm='a=imagecont(ftmachine='+'"'+self.ftmachine+'",'+'wprojplanes='+str(self.wprojplanes)+',facets='+str(self.facets)+',pixels='+str(self.imsize)+',cell='+str(self.cell)+', spw='+spwlaunch +',field='+fieldlaunch+',phasecenter='+pslaunch+',weight="'+self.weight+'", robust='+ str(self.robust)+', npixels='+str(self.npixels)+', stokes="'+self.stokes+'", numthreads='+str(self.numthreads)+', gain='+str(self.gain)+', uvtaper='+str(self.uvtaper)+', outertaper='+str(self.outertaper)+', timerange="'+str(self.timerange)+'"'+', uvrange="'+str(self.uvrange)+'"'+', baselines="'+str(self.baselines)+'"'+', scan="'+str(self.scan)+'"'+', observation="'+str(self.observation)+'"'+', pbcorr='+str(self.pbcorr)+', minpb='+str(self.minpb)+', cyclefactor='+str(self.cyclefactor)+')'
02615         print 'launch command', launchcomm
02616         self.c.pgc(launchcomm);
02617         self.c.pgc('a.visInMem='+str(self.visinmem));
02618         #c.pgc('a.painc='+str(self.painc));
02619         #c.pgc('a.cfcache='+'"'+str(self.cfcache)+'"');
02620         #c.pgc('a.pblimit='+str(self.pblimit));
02621         #c.pgc('a.dopbcorr='+str(self.dopbcorr));
02622         #c.pgc('a.applyoffsets='+str(self.applyoffsets));
02623         #c.pgc('a.epjtablename='+'"'+str(self.epjtablename)+'"');
02624 
02625     def pcontmultidiskms(self, msnames=[], workdirs=[], imagename=None, imsize=[1000, 1000], 
02626                      pixsize=['1arcsec', '1arcsec'], phasecenter='', 
02627                      field='', spw='*', ftmachine='ft', wprojplanes=128, facets=1, 
02628                      hostnames='', 
02629                      numcpuperhost=1, majorcycles=1, niter=1000, threshold='0.0mJy', alg='clark', weight='natural',
02630                      contclean=False, visinmem=False, maskimage='lala.mask',
02631                      painc=360., pblimit=0.1, dopbcorr=True, applyoffsets=False, cfcache='cfcache.dir',
02632                      epjtablename=''):
02633         if len(msnames)==0:
02634             return                    
02635         niterpercycle=niter/majorcycles
02636         if(niterpercycle == 0):
02637             niterpercycle=niter
02638             majorcycles=1
02639         c=cluster()
02640         hostname=os.getenv('HOSTNAME')
02641         wd=os.getcwd()
02642         owd=wd
02643         timesplit=0
02644         timeimage=0
02645         model=imagename+'.model' if (len(imagename) != 0) else 'elmodel'
02646         shutil.rmtree('tempmodel', True)
02647         if(not contclean):
02648             print "Removing ", model, 'and', imagename+'.image'
02649             shutil.rmtree(model, True)
02650             shutil.rmtree(imagename+'.image', True)
02651         ###num of cpu per node
02652         numcpu=numcpuperhost
02653         if((hostnames==[]) or (hostnames=='')): 
02654             hostnames=[hostname]
02655         if (len(msnames) != len(hostnames)):
02656             raise 'Number of MSs and hosts has to match for now' 
02657             return
02658         if (len(msnames) != len(workdirs)):
02659             raise 'Number of Working Directiories and hosts has to match for now' 
02660             return
02661         print 'Hosts ', hostnames
02662         time1=time.time()
02663         print 'output will be in directory', owd
02664         hostdata=odict()
02665         for k in range(len(hostnames)):
02666             hostname=hostnames[k]
02667             myrec=odict()
02668             myrec['msname']=workdirs[k]+'/'+msnames[k]
02669             myrec['spwids']=ms.msseltoindex(vis=myrec['msname'], spw=spw)['spw']
02670             myrec['spwsel'],myrec['startsel'],myrec['nchansel']=self.findchanselcont(msnames[k], myrec['spwids'], numcpuperhost)
02671             c.start_engine(hostname,numcpuperhost,workdirs[k])
02672             hostdata[hostname]=myrec
02673         ##start an engine locally for deconvolution
02674         c.start_engine(os.getenv('HOSTNAME'),1,owd)
02675         numcpu=numcpu*len(hostnames)
02676 
02677         print 'SPWSEL ', hostdata 
02678 
02679         out=range(numcpu)  
02680         c.pgc('casalog.filter()')
02681         c.pgc('from  parallel.parallel_cont import *')
02682         spwlaunch='"'+spw+'"' if (type(spw)==str) else str(spw)
02683         fieldlaunch='"'+field+'"' if (type(field) == str) else str(field)
02684         pslaunch='"'+phasecenter+'"' if (type(phasecenter) == str) else str(phasecenter)
02685         launchcomm='a=imagecont(ftmachine='+'"'+ftmachine+'",'+'wprojplanes='+str(wprojplanes)+',facets='+str(facets)+',pixels='+str(imsize)+',cell='+str(pixsize)+', spw='+spwlaunch +',field='+fieldlaunch+',phasecenter='+pslaunch+',weight="'+weight+'")'
02686         print 'launch command', launchcomm
02687         c.pgc(launchcomm);
02688         #c.pgc('a.visInMem='+str(visinmem));
02689         #c.pgc('a.painc='+str(painc));
02690         #c.pgc('a.cfcache='+'"'+str(cfcache)+'"');
02691         #c.pgc('a.pblimit='+str(pblimit));
02692         #c.pgc('a.dopbcorr='+str(dopbcorr));
02693         #c.pgc('a.applyoffsets='+str(applyoffsets));
02694         #c.pgc('a.epjtablename='+'"'+str(epjtablename)+'"');
02695 
02696         allfreq=np.array([])
02697         for msid in range(len(hostnames)): 
02698             msname=msnames[msid]
02699             myrec=hostdata.values()[msid]
02700             tb.open(msname)
02701             spectable=string.split(tb.getkeyword('SPECTRAL_WINDOW'))
02702             if(len(spectable) ==2):
02703                 spectable=spectable[1]
02704             else:
02705                 spectable=msname+"/SPECTRAL_WINDOW"
02706             tb.open(spectable)
02707             #freqs=tb.getcol('CHAN_FREQ')
02708             for k in range(len(myrec['spwids'])) :
02709                 allfreq=np.append(allfreq, tb.getcol('CHAN_FREQ',myrec['spwids'][k],1))
02710             tb.done()
02711         ##number of channels in the ms
02712         numchanperms=len(allfreq)
02713         #print 'number of channels', numchanperms
02714         minfreq=np.min(allfreq)
02715         band=np.max(allfreq)-minfreq
02716         ##minfreq=minfreq-(band/2.0);
02717         ###need to get the data selection for each process here
02718         ## this algorithm is a poor first try
02719         if(minfreq <0 ):
02720             minfreq=0.0
02721         freq='"%s"'%(str((minfreq+band/2.0))+'Hz')
02722         band='"%s"'%(str((band*1.1))+'Hz')
02723         ###define image names
02724         imlist=[]
02725         cfcachelist=[]
02726         for kk in range(len(hostnames)):
02727             myrec=hostdata.values()[kk]
02728             for jj in range(numcpuperhost) :
02729                 k=kk*numcpuperhost+jj
02730                 substr=msnames[kk]+'_spw'
02731                 for u in range(len(myrec['spwsel'][jj])):
02732                     if(u==0):
02733                         substr=substr+'_'+str(myrec['spwsel'][jj][u])+'_chan_'+str(myrec['startsel'][jj][u])
02734                     else:
02735                         substr=substr+'_spw_'+str(myrec['spwsel'][jj][u])+'_chan_'+str(myrec['startsel'][jj][u])
02736                 imlist.append(workdirs[kk]+'/'+substr)
02737                 cfcachelist.append(cfcache+'_'+str(k));
02738         #####
02739         ##continue clean or not
02740         if(contclean and os.path.exists(model)):
02741             for k in range(numcpu):
02742                 self.copyimage(inimage=model, outimage=imlist[k]+'.model', init=False)
02743         else:
02744             for k in range(len(imlist)):
02745                 shutil.rmtree(imlist[k]+'.model', True)
02746         for maj in range(majorcycles):
02747             for jj in range(len(hostnames)):
02748                 msname=msnames[jj]
02749                 for kk in range(numcpuperhost):
02750                     spwsel=hostdata[hostnames[jj]]['spwsel'][kk]
02751                     startsel=hostdata[hostnames[jj]]['startsel'][kk]
02752                     nchansel=hostdata[hostnames[jj]]['nchansel'][kk]
02753                     k=jj*numcpuperhost+kk
02754                     imnam='"%s"'%(imlist[k])
02755                     #c.odo('a.cfcache='+'"'+str(cfcachelist[k])+'"',k)
02756                     runcomm='a.imagecont(msname='+'"'+msname+'", start='+str(startsel)+', numchan='+str(nchansel)+', field="'+str(field)+'", spw='+str(spwsel)+', freq='+freq+', band='+band+', imname='+imnam+')'
02757                     print 'command is ', runcomm
02758                     out[k]=c.odo(runcomm,k)
02759             over=False
02760             while(not over):
02761                 time.sleep(5)
02762                 overone=True
02763                 for k in range(numcpu):
02764                     #print 'k', k, out[k]
02765                     overone=(overone and c.check_job(out[k],False)) 
02766                     #print 'overone', overone 
02767                 over=overone
02768             residual=imagename+'.residual'
02769             psf=imagename+'.psf'
02770             psfs=range(len(imlist))
02771             residuals=range(len(imlist))
02772             restoreds=range(len(imlist))
02773             for k in range (len(imlist)):
02774                 psfs[k]=imlist[k]+'.psf'
02775                 residuals[k]=imlist[k]+'.residual'
02776                 restoreds[k]=imlist[k]+'.image'
02777             self.averimages(residual, residuals)
02778             if(maj==0):
02779                     self.averimages(psf, psfs)
02780                     if(os.path.exists(maskimage)):
02781                         self.regridimage(outimage='__lala.mask', inimage=maskimage, templateimage=residual);
02782                         shutil.rmtree(maskimage, True)
02783                         shutil.move('__lala.mask', maskimage)
02784                     else:
02785                         self.copyimage(inimage=residual, outimage=maskimage, init=True, initval=1.0);
02786                     if(not contclean or (not os.path.exists(model))):
02787                             self.copyimage(inimage=residual, outimage=model, 
02788                                     init=True, initval=0.0)    
02789             #ia.open(residual)
02790             #for k in range(1, len(imlist)) :
02791             #    ia.open(residual)
02792             #    ia.calc('"'+residual+'" +  "'+imlist[k]+'.residual"')
02793             #ia.done()
02794             #############
02795             #   ia.open(imlist[k]+'.residual')
02796             #   print 'residual of ', imlist[k], ia.statistics()['min'], ia.statistics()['max']
02797             #   ia.done()
02798             ##############
02799             ia.open(residual)
02800             print 'Residual', ia.statistics() 
02801             ia.done()
02802             ########
02803             #incremental clean...get rid of tempmodel
02804             shutil.rmtree('tempmodel', True)
02805             rundecon='a.cleancont(alg="'+str(alg)+'", niter='+str(niterpercycle)+',psf="'+psf+'", dirty="'+residual+'", model="'+'tempmodel'+'", mask="'+str(maskimage)+'", thr="'+str(threshold)+'")'
02806             print 'Deconvolution command', rundecon
02807             out[0]=c.odo(rundecon,numcpu)
02808             over=False
02809             while(not over):
02810                 time.sleep(5)
02811                 over=c.check_job(out[0],False)
02812             ###incremental added to total 
02813             ia.open(model)
02814             ia.calc('"'+model+'" +  "tempmodel"')
02815             ia.done()
02816             ia.open('tempmodel')
02817             #arr=ia.getchunk()
02818             print 'min max of incrmodel', ia.statistics()['min'], ia.statistics()['max']
02819             ia.done()
02820             ia.open(model)
02821             #arr=ia.getchunk()
02822             print 'min max of model', ia.statistics()['min'], ia.statistics()['max']
02823             ia.done()
02824             for k in range(len(imlist)):
02825                 ia.open(imlist[k]+'.model')
02826                 #ia.putchunk(arr)
02827                 ia.insert(infile=model, locate=[0,0,0,0])
02828                 ia.done()
02829         #######
02830         restored=imagename+'.image'
02831         self.averimages(restored, restoreds)
02832         ###close major cycle loop
02833         #shutil.rmtree(imagename+'.image', True)
02834         #shutil.move(restored,  imagename+'.image')
02835         time2=time.time()
02836         ###Clean up
02837         #for k in range(len(imlist)):
02838         #   shutil.rmtree(imlist[k]+'.model', True)
02839         #  shutil.rmtree(imlist[k]+'.residual', True)
02840         # shutil.rmtree(imlist[k]+'.image', True)
02841         print 'Time to image is ', (time2-time1)/60.0, 'mins'
02842         #c.stop_cluster()
02843 
02844 
02845         return