Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_sdintimaging.py: 75%
337 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
1################################################
2# single dish + interfermeter join image reconstruction task
3#
4#
5################################################
7import platform
8import os
9import shutil
10import numpy
11import copy
12import time
14from casatasks import casalog
16from casatasks.private.imagerhelpers.imager_base import PySynthesisImager
17from casatasks.private.imagerhelpers.imager_parallel_continuum import PyParallelContSynthesisImager
18from casatasks.private.imagerhelpers.imager_parallel_cube import PyParallelCubeSynthesisImager
19from casatasks.private.imagerhelpers.input_parameters import ImagerParameters
20#from casatasks import imregrid
21from .cleanhelper import write_tclean_history, get_func_params
22from .sdint_helper import *
23from casatools import table
24from casatools import synthesisimager,synthesisutils
26try:
27 from casampi.MPIEnvironment import MPIEnvironment
28 from casampi import MPIInterface
29 mpi_available = True
30except ImportError:
31 mpi_available = False
33sdintlib = SDINT_helper()
34synu = synthesisutils()
36# setup functions
37def setup_imagerObj(paramList=None):
38 """
39 setup imaging parameters
40 """
41 defaultconstructor = False
42 if paramList!=None:
43 if not isinstance(paramList, ImagerParameters):
44 raise RuntimeError("Internal Error: invalid paramList")
45 else:
46 defaultconstructor = True
48 if defaultconstructor:
49 return PySynthesisImager
50 else:
51 return PySynthesisImager(params=paramList)
54def setup_imager(imagename, specmode,calcres,calpsf,inparams):
55 """
56 Setup cube imaging for major cycles.
57 - Do initialization
58 - and run a major cycle
59 """
60 # create a local copy of input params dict so that it can be
61 # modified
62 locparams = copy.deepcopy(inparams)
64 # cube imaging setup
65 locparams['imagename']=imagename
66 locparams['specmode']='cube'
67 locparams['niter']=0
68 locparams['deconvolver']='hogbom'
70 #casalog.post("local inparams(msname) in setup_imager==",locparams['msname'])
71 params = ImagerParameters(**locparams)
72 #params = ImagerParameters(msname=self.vis, field=self.field,spw=self.spw,
73 # imagename=imagename,
74 # imsize=self.imsize, cell=self.cell, phasecenter=self.phasecenter,
75 # weighting=self.weighting,
76 # gridder=self.gridder, pblimit=self.pblimit,wprojplanes=self.wprojplanes,
77 # specmode='cube',nchan=self.nchan,
78 # reffreq=self.reffreq, width=self.width, start=self.start, interpolation=self.interpolation,
79 # interactive=self.interactive,
80 # deconvolver='hogbom', niter=0,
81 # wbawp=True)
83 ## Major cycle is either PySynthesisImager or PyParallelCubeSynthesisImager
84 imagertool = setup_imagerObj(params)
86 #self.imagertool = PySynthesisImager(params=params)
87 imagertool.initializeImagers()
88 imagertool.initializeNormalizers()
89 imagertool.setWeighting()
90 if 'psfphasecenter' in locparams:
91 psfphasecenter = locparams['psfphasecenter']
92 else:
93 psfphasecenter = ''
95 ## Extra one for psfphasecenter...
96 imagerInst=None
97 if((psfphasecenter != '') and (gridder=='mosaic')):
98 imagerInst = setup_imagerObj()
101 gridder = locparams['gridder']
103 if calpsf == True:
104 imagertool.makePSF()
105 imagertool.makePB()
106 if((psfphasecenter != '') and (gridder=='mosaic')):
107 casalog.post("doing with different phasecenter psf", "INFO")
108 imagertool.unlockimages(0)
109 psfParameters=paramList.getAllPars()
110 psfParameters['phasecenter']=psfphasecenter
111 psfParamList=ImagerParameters(**psfParameters)
112 psfimager=imagerInst(params=psfParamList)
113 psfimager.initializeImagers()
114 psfimager.setWeighting()
115 psfimager.makeImage('psf', psfParameters['imagename']+'.psf')
117 # can take out this since niter is fixed to 0
118 if locparams['niter'] >=0 :
119 ## Make dirty image
120 if calcres == True:
121 t0=time.time();
122 imagertool.runMajorCycle(isCleanCycle=False)
123 t1=time.time();
124 casalog.post("***Time for major cycle (calcres=T): "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_tclean");
126 ## In case of no deconvolution iterations....
127 #if locparams['niter']==0 and calcres==False:
128 # if savemodel != "none":
129 # imagertool.predictModel()
131 sdintlib.copy_restoringbeam(fromthis=imagename+'.psf', tothis=imagename+'.residual')
132 return imagertool
134def setup_deconvolver(imagename,specmode,inparams):
135 """
136 Cube or MFS minor cycles.
137 """
138 #params = ImagerParameters(msname=self.vis, field=self.field,spw=self.spw,
139 # imagename=imagename,
140 # imsize=self.imsize, cell=self.cell, phasecenter=self.phasecenter,
141 # weighting=self.weighting,
142 # gridder=self.gridder, pblimit=self.pblimit,wprojplanes=self.wprojplanes,
143 # specmode=self.specmode,nchan=self.nchan,
144 # reffreq=self.reffreq, width=self.width, start=self.start, interpolation=self.interpolation,
145 # interactive=self.interactive,
146 # deconvolver=self.deconvolver, scales=self.scales,nterms=self.nterms,
147 # niter=self.niter,cycleniter=self.cycleniter, threshold=self.threshold,
148 # mask=self.mask)
149 inparams['imagename']=imagename
150 params = ImagerParameters(**inparams)
151 deconvolvertool = setup_imagerObj(params)
153 ## Why are we initializing these ?
154 deconvolvertool.initializeImagers()
155 deconvolvertool.initializeNormalizers()
156 deconvolvertool.setWeighting()
159 ### These three should be unncessary. Need a 'makeimage' method for csys generation.
160 deconvolvertool.makePSF() ## Make this to get a coordinate system
161 #deconvolvertool.makeImage('psf', imagename+'.psf')
162 deconvolvertool.makePB() ## Make this to turn .weight into .pb maps
164 ## Initialize deconvolvers. ( Order is important. This cleans up a leftover tablecache image.... FIX!)
165 deconvolvertool.initializeDeconvolvers()
166 deconvolvertool.initializeIterationControl() # This needs to be run before runMajorCycle
167 deconvolvertool.runMajorCycle(isCleanCycle=False) ## Make this to make template residual images.
169 return deconvolvertool
171def setup_sdimaging(template='',output='', inparms=None, sdparms=None):
172 """
173 Make the SD cube Image and PSF
175 Option 1 : Use/Regrid cubes for the observed image and PSF
176 Option 2 : Make the SD image and PSF cubes using 'tsdimager's usage of the SD gridder option.
178 Currently, only Option 1 is supported.
180 """
181 sdintlib = SDINT_helper()
182 if 'sdpsf' in sdparms:
183 sdpsf = sdparms['sdpsf']
184 else:
185 raise RuntimeError("Internal Error: missing sdpsf parameter")
187 if 'sdimage' in sdparms:
188 sdimage = sdparms['sdimage']
189 else:
190 raise RuntimeError("Internal Error: missing sdimage parameter")
191 if 'pblimit' in inparms:
192 pblimit = inparms['pblimit']
194 if sdpsf !="":
195 ## check the coordinates of psf with int psf
196 sdintlib.checkpsf(sdpsf, template+'.psf')
198 ## Regrid the input SD image and PSF cubes to the target coordinate system.
199 #imregrid(imagename=sdimage, template=template+'.residual',
200 # output=output+'.residual',overwrite=True,axes=[0,1])
201 sdintlib.regridimage(imagename=sdimage, template=template+'.residual', outfile=output+'.residual')
202 #imregrid(imagename=sdimage, template=template+'.residual',
203 # output=output+'.image',overwrite=True,axes=[0,1])
204 sdintlib.regridimage(imagename=sdimage, template=template+'.residual', outfile=output+'.image')
206 if sdpsf !="":
207 #imregrid(imagename=sdpsf, template=template+'.psf',
208 # output=output+'.psf',overwrite=True,axes=[0,1])
209 sdintlib.regridimage(imagename=sdpsf, template=template+'.psf', outfile=output+'.psf')
210 else:
211 ## Make an internal sdpsf image if the user has not supplied one.
212 casalog.post("Constructing a SD PSF cube by evaluating Gaussians based on the restoring beam information in the regridded SD Image Cube")
213 sdintlib.create_sd_psf(sdimage=output+'.residual', sdpsfname=output+'.psf')
215 ## Apply the pbmask from the INT image cube, to the SD cubes.
216 #TTB: Create *.mask cube
218 sdintlib.addmask(inpimage=output+'.residual', pbimage=template+'.pb', pblimit=pblimit)
219 sdintlib.addmask(inpimage=output+'.image', pbimage=template+'.pb', pblimit=pblimit)
223def feather_residual(int_cube, sd_cube, joint_cube, applypb, inparm):
225 if applypb==True:
226 ## Take initial INT_dirty image to flat-sky.
227 sdintlib.modify_with_pb(inpcube=int_cube+'.residual',
228 pbcube=int_cube+'.pb',
229 cubewt=int_cube+'.sumwt',
230 chanwt=inparm['chanwt'],
231 action='div',
232 pblimit=inparm['pblimit'],
233 freqdep=True)
235 ## Feather flat-sky INT dirty image with SD image
236 sdintlib.feather_int_sd(sdcube=sd_cube+'.residual',
237 intcube=int_cube+'.residual',
238 jointcube=joint_cube+'.residual', ## output
239 sdgain=inparm['sdgain'],
240 dishdia=inparm['dishdia'],
241 usedata=inparm['usedata'],
242 chanwt = inparm['chanwt'])
244 if applypb==True:
245 if inparm['specmode'].count('cube')>0:
246 ## Multiply the new JOINT dirty image by the frequency-dependent PB.
247 fdep_pb = True
248 else:
249 ## Multiply new JOINT dirty image by a common PB to get the effect of conjbeams.
250 fdep_pb = False
251 sdintlib.modify_with_pb(inpcube=joint_cube+'.residual',
252 pbcube=int_cube+'.pb',
253 cubewt=int_cube+'.sumwt',
254 chanwt=inparm['chanwt'],
255 action='mult',
256 pblimit=inparm['pblimit'],
257 freqdep=fdep_pb)
259def deleteTmpFiles():
260 if os.path.exists('tmp_intplane'):
261 os.system('rm -rf tmp_intplane')
262 if os.path.exists('tmp_sdplane'):
263 os.system('rm -rf tmp_sdplane')
264 if os.path.exists('tmp_jointplane'):
265 os.system('rm -rf tmp_jointplane')
268def sdintimaging(
269 usedata,
270 ####### Single dish input data
271 sdimage,
272 sdpsf,
273 sdgain,
274 dishdia,
275 ####### Interferometer Data Selection
276 vis,#='',
277 selectdata,
278 field,#='',
279 spw,#='',
280 timerange,#='',
281 uvrange,#='',
282 antenna,#='',
283 scan,#='',
284 observation,#='',
285 intent,#='',
286 datacolumn,#='corrected',
289 ####### Image definition
290 imagename,#='',
291 imsize,#=[100,100],
292 cell,#=['1.0arcsec','1.0arcsec'],
293 phasecenter,#='J2000 19:59:28.500 +40.44.01.50',
294 stokes,#='I',
295 projection,#='SIN',
296 startmodel,#='',
298 ## Spectral parameters
299 specmode,#='mfs',
300 reffreq,#='',
301 nchan,#=1,
302 start,#='',
303 width,#='',
304 outframe,#='LSRK',
305 veltype,#='',
306 restfreq,#=[''],
307# sysvel,#='',
308# sysvelframe,#='',
309 interpolation,#='',
310# chanchunks,#=1,
311 perchanweightdensity, #=''
312 ##
313 ####### Gridding parameters
314 gridder,#='ft',
315 facets,#=1,
316 psfphasecenter,#='',
318 wprojplanes,#=1,
320 ### PB
321 vptable,
322 mosweight, #=True
323 aterm,#=True,
324 psterm,#=True,
325 wbawp ,#= True,
326# conjbeams ,#= True,
327 cfcache ,#= "",
328 usepointing, #=false
329 computepastep ,#=360.0,
330 rotatepastep ,#=360.0,
331 pointingoffsetsigdev ,#=0.0,
333 pblimit,#=0.01,
334# normtype,#='flatnoise',
336 ####### Deconvolution parameters
337 deconvolver,#='hogbom',
338 scales,#=[],
339 nterms,#=1,
340 smallscalebias,#=0.0
342 ### restoration options
343 restoration,
344 restoringbeam,#=[],
345 pbcor,
347 ##### Outliers
348# outlierfile,#='', ### RESTRICTION : No support for outlier fields for joint SD-INT imaging.
350 ##### Weighting
351 weighting,#='natural',
352 robust,#=0.5,
353 noise,#0.0Jy
354 npixels,#=0,
355# uvtaper,#=False,
356 uvtaper,#=[],
359 ##### Iteration control
360 niter,#=0,
361 gain,#=0.1,
362 threshold,#=0.0,
363 nsigma,#=0.0
364 cycleniter,#=0,
365 cyclefactor,#=1.0,
366 minpsffraction,#=0.1,
367 maxpsffraction,#=0.8,
368 interactive,#=False,
369 fullsummary,#=False,
370 nmajor,#=-1,
372 ##### (new) Mask parameters
373 usemask,#='user',
374 mask,#='',
375 pbmask,#='',
376 # maskthreshold,#='',
377 # maskresolution,#='',
378 # nmask,#=0,
380 ##### automask by multithresh
381 sidelobethreshold,#=5.0,
382 noisethreshold,#=3.0,
383 lownoisethreshold,#=3.0,
384 negativethreshold,#=0.0,
385 smoothfactor,#=1.0,
386 minbeamfrac,#=0.3,
387 cutthreshold,#=0.01,
388 growiterations,#=100
389 dogrowprune,#=True
390 minpercentchange,#=0.0
391 verbose, #=False
392 fastnoise, #=False
394 ## Misc
396 restart,#=True,
398 #savemodel,#="none",
400# makeimages,#="auto"
401 calcres,#=True,
402 calcpsf):#=True,
404 ####### State parameters
405 #parallel):#=False):
408 ##################################################
409 # copied from SDINT.do_reconstruct
410 #################################################
411 int_cube = imagename+'.int.cube'
412 sd_cube = imagename+'.sd.cube'
413 joint_cube = imagename+'.joint.cube'
414 joint_multiterm = imagename+'.joint.multiterm'
416 if specmode=='mfs':
417 decname = joint_multiterm
418 else:
419 decname = joint_cube
421 #####################################################
422 #### Sanity checks and controls
423 #####################################################
425 if interactive:
426 # Check for casaviewer, if it does not exist flag it up front for macOS
427 # since casaviewer is no longer provided by default with macOS.
428 try:
429 import casaviewer as __test_casaviewer
430 except:
431 if platform.system( ) == "Darwin":
432 casalog.post(
433 "casaviewer is no longer available for macOS, for more information see: http://go.nrao.edu/casa-viewer-eol Please restart by setting interactive=F",
434 "WARN",
435 "task_sdintimaging",
436 )
437 raise RuntimeError( "casaviewer is no longer available for macOS, for more information see: http://go.nrao.edu/casa-viewer-eol" )
439 ### Move these checks elsewhere ?
440 inpparams=locals().copy()
441 ###now deal with parameters which are not the same name
442 #casalog.post("current inpparams=",inpparams)
443 #casalog.post("inpparams.keys()=",inpparams.keys())
444 locvis=inpparams.pop('vis')
445 #casalog.post("LOCVIS====",locvis)
446 if type(locvis)==list:
447 llocvis = [v.lstrip() for v in locvis]
448 else:
449 llocvis = locvis.lstrip()
450 inpparams['msname']=llocvis
451 inpparams['timestr']= inpparams.pop('timerange')
452 inpparams['uvdist']= inpparams.pop('uvrange')
453 inpparams['obs']= inpparams.pop('observation')
454 inpparams['state']= inpparams.pop('intent')
455 inpparams['loopgain']=inpparams.pop('gain')
456 inpparams['scalebias']=inpparams.pop('smallscalebias')
458 sdparms={}
459 sdparms['sdimage']=inpparams['sdimage']
460 sdparms['sdpsf']=inpparams['sdpsf']
461 sdparms['sdgain']=inpparams['sdgain']
463 if usedata!='int': # check sd parameters
465 _myia = image()
467 if not os.path.exists(sdparms['sdimage']):
468 casalog.post( "Input image sdimage = '"+str(sdparms['sdimage'])+"' does not exist.", "WARN", "task_sdintimaging" )
469 return
470 else:
471 try:
472 _myia.open(sdparms['sdimage'])
473 except Exception as instance:
474 casalog.post( "Input image sdimage = '"+str(sdparms['sdimage'])+"' cannot be opened.", "WARN", "task_sdintimaging" )
475 casalog.post( str(instance), "WARN", "task_sdintimaging" )
476 return
478 mysummary = _myia.summary(list=False)
479 _myia.close()
481 try:
482 freqaxis_index = list(mysummary['axisnames']).index('Frequency')
483 except(ValueError):
484 casalog.post('The image '+sdparms['sdimage']+' has no frequency axis. Try adding one with ia.adddegaxis() .',
485 'WARN', 'task_sdintimaging')
486 return
488 if freqaxis_index != 3:
489 casalog.post('The image '+sdparms['sdimage']+' has its frequency axis on position '+str(freqaxis_index)+
490 ' whereas it should be in position 3 (counting from 0). Use task imtrans() with order=["r", "d", "s", "f"] to fix this.',
491 'WARN', 'task_sdintimaging')
492 return
495 if sdparms['sdpsf']!='':
496 if not os.path.exists(sdparms['sdpsf']):
497 casalog.post( "Input image sdpsf = '"+str(sdparms['sdpsf'])+"' does not exist.", "WARN", "task_sdintimaging" )
498 return
499 else:
500 try:
501 _myia.open(sdparms['sdpsf'])
502 _myia.close()
503 except Exception as instance:
504 casalog.post( "Input image sdpsf = '"+str(sdparms['sdpsf'])+"' cannot be opened.", "WARN", "task_sdintimaging" )
505 casalog.post( str(instance), "WARN", "task_sdintimaging" )
506 return
508 if (sdparms['sdgain']*0!=0 or sdparms['sdgain']<=0):
509 casalog.post('Invalid sdgain: '+str(sdparms['sdgain']), 'WARN')
510 casalog.post("The sdgain parameter needs to be chosen as a number > 0 which represents the weight of the SD contribution relative to the INT contribution to the joint image.", "WARN", "task_sdintimaging")
511 return
513 if (dishdia*0!=0 or dishdia<=0):
514 casalog.post('Invalid dishdia: '+str(dishdia), 'WARN')
515 casalog.post("The dishdia parameter needs to provide the diameter (meters) of the SD telescope which produced the SD image.", "WARN", "task_sdintimaging")
516 return
519 if specmode=='cont':
520 specmode='mfs'
521 inpparams['specmode']='mfs'
523 # from sdint
524 # automatically decide if pb need to be applied
525 if gridder=='mosaic' or gridder=='awproject':
526 applypb = True
527 else:
528 applypb = False
530 if (deconvolver=="mtmfs") and (specmode!='mfs') and (specmode!='cube' or nterms!=1) and (specmode!='cubedata' or nterms!=1):
531 casalog.post( "The MSMFS algorithm (deconvolver='mtmfs') applies only to specmode='mfs' or specmode='cube' with nterms=1 or specmode='cubedata' with nterms=1.", "WARN", "task_sdintimaging" )
532 return
534 if(deconvolver=="mtmfs" and (specmode=='cube' or specmode=='cubedata') and nterms==1 ):
535 casalog.post( "The MSMFS algorithm (deconvolver='mtmfs') with specmode='cube', nterms=1 is currently not supported. Please use deconvolver='multiscale' instead for cubes.", "WARN", "task_sdintimaging" )
536 return
538 if(specmode=='mfs' and deconvolver!='mtmfs'):
539 casalog.post("Currently, only the multi-term MFS algorithm is supported for specmode=mfs. To make a single plane MFS image (while retaining the frequency dependence for the cube major cycle stage), please pick nterms=1 along with deconvolver=mtmfs. The scales parameter is still usable for multi-scale multi-term deconvolution","WARN","task_sdintimaging")
540 return;
542 if(usedata=='sd'):
543 casalog.post("The Single-Dish-Only mode of sdintimaging is better supported via the deconvolve task which supports spectral cube, mfs and multi-term mfs deconvolution in the image domain alone. The deconvolve task is the more appropriate version to use for stand-alone image-domain deconvolution, and will not have the bookkeeping overheads currently present in the sdintimaging task's sd-only mode. Please note that the 'sd' option of the sdintimaging task will be removed in a subsequent release. Please refer to the task deconvolve documentation for instructions on how to prepare image and psf cubes for the deconvolve task for all these modes.","WARN","task_sdintimaging");
545 if (nmajor < -1):
546 casalog.post("Negative values less than -1 for nmajor are reserved for possible future implementation", "WARN", "task_sdintimaging")
547 return
549# if parallel==True:
550# casalog.post("Cube parallelization (all major cycles) is currently not supported via task_sdintimaging. This will be enabled after a cube parallelization rework.")
551# return;
553 #####################################################
554 #### Construct ImagerParameters object
555 #####################################################
557 imager = None
558 paramList = None
559 deconvolvertool = None
561 # Put all parameters into dictionaries and check them.
562 ##make a dictionary of parameters that ImagerParameters take
563 defparm=dict(list(zip(ImagerParameters.__init__.__code__.co_varnames[1:], ImagerParameters.__init__.__defaults__)))
566 ###assign values to the ones passed to tclean and if not defined yet in tclean...
567 ###assign them the default value of the constructor
568 bparm={k: inpparams[k] if k in inpparams else defparm[k] for k in defparm.keys()}
570 ###default mosweight=True is tripping other gridders as they are not
571 ###expecting it to be true
572 if(bparm['mosweight']==True and bparm['gridder'].find("mosaic") == -1):
573 bparm['mosweight']=False
575 ## Two options have been removed from the interface. Hard-code them here.
576 bparm['normtype'] = 'flatnoise' ## Hard-code this since the pbcor steps assume it.
577 bparm['conjbeams']=False
579 #paramList=ImagerParameters(**bparm)
581 #paramList.printParameters()
583 if len(pointingoffsetsigdev)>0 and pointingoffsetsigdev[0]!=0.0 and usepointing==True and gridder.count('awproj')>1:
584 casalog.post("pointingoffsetsigdev will be used for pointing corrections with AWProjection", "WARN")
586 #=========================================================
587 ####set the children to load c++ libraries and applicator
588 ### make workers ready for c++ based mpicommands
589 cppparallel=False
590 if mpi_available and MPIEnvironment.is_mpi_enabled:
591 mint=MPIInterface.MPIInterface()
592 cl=mint.getCluster()
593 cl._cluster.pgc("from casatools import synthesisimager", False)
594 cl._cluster.pgc("si=synthesisimager()", False)
596 cl._cluster.pgc("si.initmpi()", False)
597 cppparallel=True
598 ###ignore chanchunk
599 bparm['chanchunks']=1
601 #################################################
602 #### start of more computing-intensive work #####
603 #################################################
605 retrec={}
607 try:
608 sdintlib = SDINT_helper()
609 ## Init major cycle elements
610 ### debug (remove it later)
611 casalog.post("INT cube setup ....")
612 t0=time.time();
613 imager=setup_imager(int_cube, specmode, calcres, calcpsf, bparm)
616 ##imager.initializeImagers()
618 ##imager.initializeNormalizers()
619 ##imager.setWeighting()
620 t1=time.time();
621 casalog.post("***Time for initializing imager (INT cube) : "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_sdintimaging");
623 ## Init minor cycle elements
624 if niter>0 or restoration==True:
625 ### debug (remove it later)
626 casalog.post("Combined image setup ....")
627 t0=time.time();
628 deconvolvertool=setup_deconvolver(decname, specmode, bparm )
629 #imager.initializeDeconvolvers()
630 t1=time.time();
631 casalog.post("***Time for seting up deconvolver(s): "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_sdintimaging");
633 if usedata!='int':
634 ### debug (remove it later)
635 casalog.post("SD cube setup ....")
636 setup_sdimaging(template=int_cube, output=sd_cube, inparms=bparm, sdparms=sdparms )
639 ####now is the time to check estimated memory
640 # need to move to somewhere below???
641 imager.estimatememory()
643 ## Do sanity checks on INT and SD cubes
644 ### Frequency range of cube, data selection range. mtmfs reffreq.
645 ### nchan too small or too large
646 ### sumwt : flagged channels in int cubes
647 ### sd cube empty channels ? Weight image ?
648 validity, inpparams = sdintlib.check_coords(intres=int_cube+'.residual', intpsf=int_cube+'.psf',
649 intwt = int_cube+'.sumwt',
650 sdres=sd_cube+'.residual', sdpsf=sd_cube+'.psf',
651 sdwt = '',
652 pars=inpparams)
654 if validity==False:
655 casalog.post('Exiting from the sdintimaging task due to inconsistencies found between the interferometer-only and singledish-only image and psf cubes. Please modify inputs as needed','WARN')
656 if imager != None:
657 imager.deleteTools()
658 if deconvolvertool != None:
659 deconvolvertool.deleteTools()
660 deleteTmpFiles()
661 return
663 #### inpparams now has a new parameter "chanwt" with ones and zeros to indicate chans
664 #### that have data from both INT and SD cubes (they are the 'ones'). This is to be used in
665 #### feathering and in the cube-to-taylor sum and modify_with_pb.
667 #### END MOVE THIS SECTION to setup_imager ....for sdint
669 #### SDINT specific feathering....
670 ## Feather INT and SD residual images (feather in flat-sky. output has common PB)
671 casalog.post("Feathering INT and SD residual images...")
672 feather_residual(int_cube, sd_cube, joint_cube, applypb, inpparams)
673 sdintlib.feather_int_sd(sdcube=sd_cube+'.psf',
674 intcube=int_cube+'.psf',
675 jointcube=joint_cube+'.psf',
676 sdgain=sdgain,
677 dishdia=dishdia,
678 usedata=usedata,
679 chanwt = inpparams['chanwt'])
681 #print("Fitting for cube")
682 synu.fitPsfBeam(joint_cube)
684 ###############
685 ##### Placeholder code for PSF renormalization if needed
686 ##### Note : If this is enabled, we'll need to restrict the use of 'faceting' as .sumwt shape changes.
687 #sdintlib.calc_renorm(intname=int_cube, jointname=joint_cube)
688 #sdintlib.apply_renorm(imname=joint_cube+'.psf', sumwtname=joint_cube+'.sumwt')
689 #sdintlib.apply_renorm(imname=joint_cube+'.residual', sumwtname=joint_cube+'.sumwt')
690 ###############
692 #casalog.post("feather_int_sd DONE")
694 if specmode=='mfs':
695 ## Calculate Spectral PSFs and Taylor Residuals
696 casalog.post("Calculate spectral PSFs and Taylor Residuals...")
697 sdintlib.cube_to_taylor_sum(cubename=joint_cube+'.psf',
698 cubewt=int_cube+'.sumwt',
699 chanwt=inpparams['chanwt'],
700 mtname=joint_multiterm+'.psf',
701 nterms=nterms, reffreq=inpparams['reffreq'], dopsf=True)
702 sdintlib.cube_to_taylor_sum(cubename=joint_cube+'.residual',
703 cubewt=int_cube+'.sumwt',
704 chanwt=inpparams['chanwt'],
705 mtname=joint_multiterm+'.residual',
706 nterms=nterms, reffreq=inpparams['reffreq'], dopsf=False)
708 #print("Fit for multiterm")
709 if(deconvolver=='mtmfs' and nterms==1): # work around file naming issue
710 os.system('rm -rf '+joint_multiterm+'tmp.psf')
711 os.system('ln -sf '+joint_multiterm+'.psf.tt0 '+joint_multiterm+'tmp.psf')
712 synu.fitPsfBeam(joint_multiterm+'tmp',nterms=nterms)
713 os.system('rm -rf '+joint_multiterm+'tmp.psf')
714 else:
715 synu.fitPsfBeam(joint_multiterm,nterms=nterms)
717 if niter>0 :
718 isit = deconvolvertool.hasConverged()
719 deconvolvertool.updateMask()
721 while ( not deconvolvertool.hasConverged() ):
723 t0=time.time();
725 #### Print out the peak of the residual image here to check !!!
726# if specmode=='mfs':
727# print('Max of joint residual before initminorcycle' + str(imstat(joint_multiterm+'.residual.tt0',verbose=False)['max'][0]))
728# else:
729# print('Max of joint residual before initminorcycle' + str(imstat(joint_cube+'.residual',verbose=False)['max'][0]))
733 deconvolvertool.runMinorCycle()
735# if specmode=='mfs':
736# print('Max of joint residual after minorcycle' + str(imstat(joint_multiterm+'.residual.tt0',verbose=False)['max'][0]))
737# else:
738# print('Max of joint residual after minorcycle' + str(imstat(joint_cube+'.residual',verbose=False)['max'][0]))
741 t1=time.time();
742 casalog.post("***Time for minor cycle: "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_sdintimaging");
744 ### sdint specific feathering steps HERE
745 ## Prepare the joint model cube for INT and SD major cycles
746 if specmode=='mfs':
747 ## Convert Taylor model coefficients into a model cube : int_cube.model
748 sdintlib.taylor_model_to_cube(cubename=int_cube, ## output
749 mtname=joint_multiterm, ## input
750 nterms=nterms, reffreq=inpparams['reffreq'])
751 else:
752 ## Copy the joint_model cube to the int_cube.model
753 shutil.rmtree(int_cube+'.model',ignore_errors=True)
754 shutil.copytree(joint_cube+'.model', int_cube+'.model')
755 hasfile=os.path.exists(joint_cube+'.model')
756 #casalog.post("DEBUG: has joint cube .image===",hasfile)
758 if applypb==True:
759 ## Take the int_cube.model to flat sky.
760 if specmode=='cube':
761 ## Divide the model by the frequency-dependent PB to get to flat-sky
762 fdep_pb = True
763 else:
764 ## Divide the model by the common PB to get to get to flat-sky
765 fdep_pb = False
766 sdintlib.modify_with_pb(inpcube=int_cube+'.model',
767 pbcube=int_cube+'.pb',
768 cubewt=int_cube+'.sumwt',
769 chanwt=inpparams['chanwt'],
770 action='div', pblimit=pblimit,freqdep=fdep_pb)
772 if usedata!="int":
773 ## copy the int_cube.model to the sd_cube.model
774 shutil.rmtree(sd_cube+'.model',ignore_errors=True)
775 shutil.copytree(int_cube+'.model', sd_cube+'.model')
777 if applypb==True:
778 ## Multiply flat-sky model with freq-dep PB
779 sdintlib.modify_with_pb(inpcube=int_cube+'.model',
780 pbcube=int_cube+'.pb',
781 cubewt=int_cube+'.sumwt',
782 chanwt=inpparams['chanwt'],
783 action='mult', pblimit=pblimit, freqdep=True)
785 ## Major cycle for interferometer data
786 t0=time.time();
787 # print('Max of int residual before major cycle' + str(imstat(int_cube+'.residual',verbose=False)['max'][0]))
788 # print('Max of int model before major cycle' + str(imstat(int_cube+'.model',verbose=False)['max'][0]))
790 if usedata != "sd":
791 imager.runMajorCycle()
792 # track nmajor for the deconvolvertool.hasConverged() method
793 deconvolvertool.majorCnt = imager.majorCnt
795 # print('Max of int residual after major cycle' + str(imstat(int_cube+'.residual',verbose=False)['max'][0]))
796 t1=time.time();
797 casalog.post("***Time for major cycle: "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_tclean");
799 if usedata!="int":
800 ## Major cycle for Single Dish data (uses the flat sky cube model in sd_cube.model )
801 sdintlib.calc_sd_residual(origcube=sd_cube+'.image',
802 modelcube=sd_cube+'.model',
803 residualcube=sd_cube+'.residual', ## output
804 psfcube=sd_cube+'.psf')
806 ## Feather the residuals
807 feather_residual(int_cube, sd_cube, joint_cube, applypb, inpparams )
808 ###############
809 ##### Placeholder code for PSF renormalization if needed
810 #sdintlib.apply_renorm(imname=joint_cube+'.residual', sumwtname=joint_cube+'.sumwt')
811 ###############
813 if specmode=='mfs':
814 ## Calculate Spectral Taylor Residuals
815 sdintlib.cube_to_taylor_sum(cubename=joint_cube+'.residual',
816 cubewt=int_cube+'.sumwt',
817 chanwt=inpparams['chanwt'],
818 mtname=joint_multiterm+'.residual',
819 nterms=nterms, reffreq=inpparams['reffreq'], dopsf=False)
821# if specmode=='mfs':
822# print('Max of residual after feather step ' + str(imstat(joint_multiterm+'.residual.tt0',verbose=False)['max'][0]))
823# else:
824# print('Max of residual after feather step ' + str(imstat(joint_cube+'.residual',verbose=False)['max'][0]))
827 deconvolvertool.updateMask()
829 ## Get summary from iterbot
830 #if type(interactive) != bool:
831 #retrec=imager.getSummary();
832 retrec=deconvolvertool.getSummary(fullsummary);
833 retrec['nmajordone'] = imager.majorCnt
834 if calcres==True:
835 retrec['nmajordone'] = retrec['nmajordone'] + 1 ## To be consistent with tclean. Remove, when we can change the meaning of nmajordone to exclude the initial major cycles.
837 ## Restore images.
838 if restoration==True:
839 t0=time.time();
840 deconvolvertool.restoreImages()
841 t1=time.time();
842 casalog.post("***Time for restoring images: "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_tclean");
843 if pbcor==True:
844 #if applypb==True:
845 t0=time.time();
846 if specmode=='mfs':
847 sdintlib.pbcor(imagename=decname+'.image.tt0' , pbimage=decname+'.pb.tt0' , cutoff=pblimit,outfile=decname+'.image.tt0.pbcor')
848 else:
849 sdintlib.pbcor(imagename=joint_cube+'.image' , pbimage=int_cube+'.pb' , cutoff=pblimit,outfile=joint_cube+'.image.pbcor')
851 #imager.pbcorImages()
852 t1=time.time();
853 casalog.post("***Time for pb-correcting images: "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_tclean");
855 ##close tools
856 # needs to deletools before concat or lock waits for ever
857 imager.deleteTools()
858 deconvolvertool.deleteTools()
861 finally:
862 if imager != None:
863 imager.deleteTools()
864 if(cppparallel):
865 ###release workers back to python mpi control
866 si=synthesisimager()
867 si.releasempi()
869 #clean up tmp files
870 deleteTmpFiles()
872 # Write history at the end, when hopefully all temp files are gone from disk,
873 # so they won't be picked up. They need time to disappear on NFS or slow hw.
874 # Copied from tclean.
875 try:
876 params = get_func_params(sdintimaging, locals())
877 write_tclean_history(imagename, 'sdintimaging', params, casalog)
878 except Exception as exc:
879 casalog.post("Error updating history (logtable): {} ".format(exc),'WARN')
881 return retrec
883##################################################