Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_tclean.py: 77%
258 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
1################################################
2# Refactored Clean task
3#
4# v1.0: 2012.10.05, U.R.V.
5#
6################################################
8import platform
9import os
10import shutil
11import numpy
12import copy
13import filecmp
14import time
15import pdb
17from casatasks import casalog
19from casatasks.private.imagerhelpers.imager_base import PySynthesisImager
20from casatasks.private.imagerhelpers.input_parameters import saveparams2last
21from casatasks.private.imagerhelpers.imager_parallel_continuum import PyParallelContSynthesisImager
22from casatasks.private.imagerhelpers.imager_parallel_cube import PyParallelCubeSynthesisImager
23from casatasks.private.imagerhelpers.imager_mtmfs_via_cube import PyMtmfsViaCubeSynthesisImager
24from casatasks.private.imagerhelpers.input_parameters import ImagerParameters
25from casatasks.private.imagerhelpers.imager_return_dict import ImagingDict
26from .cleanhelper import write_tclean_history, get_func_params
27from casatools import table
28from casatools import image
29from casatools import synthesisutils
30from casatools import synthesisimager
33try:
34 from casampi.MPIEnvironment import MPIEnvironment
35 from casampi import MPIInterface
37 mpi_available = True
38except ImportError:
39 mpi_available = False
41# if you want to save tclean.last.* from python call of tclean uncomment the decorator
42#@saveparams2last(multibackup=True)
44def tclean(
45 ####### Data Selection
46 vis, # ='',
47 selectdata,
48 field, # ='',
49 spw, # ='',
50 timerange, # ='',
51 uvrange, # ='',
52 antenna, # ='',
53 scan, # ='',
54 observation, # ='',
55 intent, # ='',
56 datacolumn, # ='corrected',
57 ####### Image definition
58 imagename, # ='',
59 imsize, # =[100,100],
60 cell, # =['1.0arcsec','1.0arcsec'],
61 phasecenter, # ='J2000 19:59:28.500 +40.44.01.50',
62 stokes, # ='I',
63 projection, # ='SIN',
64 startmodel, # ='',
65 ## Spectral parameters
66 specmode, # ='mfs',
67 reffreq, # ='',
68 nchan, # =1,
69 start, # ='',
70 width, # ='',
71 outframe, # ='LSRK',
72 veltype, # ='',
73 restfreq, # =[''],
74 # sysvel,#='',
75 # sysvelframe,#='',
76 interpolation, # ='',
77 perchanweightdensity, # =''
78 ##
79 ####### Gridding parameters
80 gridder, # ='ft',
81 facets, # =1,
82 psfphasecenter, # ='',
83 wprojplanes, # =1,
84 ### PB
85 vptable,
86 mosweight, # =True
87 aterm, # =True,
88 psterm, # =True,
89 wbawp, # = True,
90 conjbeams, # = True,
91 cfcache, # = "",
92 usepointing, # =false
93 computepastep, # =360.0,
94 rotatepastep, # =360.0,
95 pointingoffsetsigdev, # =[10.0],
96 pblimit, # =0.01,
97 normtype, # ='flatnoise',
98 ####### Deconvolution parameters
99 deconvolver, # ='hogbom',
100 scales, # =[],
101 nterms, # =1,
102 smallscalebias, # =0.0
103 fusedthreshold, # =0.0
104 largestscale, # =-1
105 ### restoration options
106 restoration,
107 restoringbeam, # =[],
108 pbcor,
109 ##### Outliers
110 outlierfile, # ='',
111 ##### Weighting
112 weighting, # ='natural',
113 robust, # =0.5,
114 noise, # 0.0Jy
115 npixels, # =0,
116 # uvtaper,#=False,
117 uvtaper, # =[],
118 ##### Iteration control
119 niter,#=0,
120 gain,#=0.1,
121 threshold,#=0.0,
122 nsigma,#=0.0
123 cycleniter,#=0,
124 cyclefactor,#=1.0,
125 minpsffraction,#=0.1,
126 maxpsffraction,#=0.8,
127 interactive,#=False,
128 nmajor,#=-1,
129 fullsummary,#=False,
131 ##### (new) Mask parameters
132 usemask, # ='user',
133 mask, # ='',
134 pbmask, # ='',
135 # maskthreshold,#='',
136 # maskresolution,#='',
137 # nmask,#=0,
138 ##### automask by multithresh
139 sidelobethreshold, # =5.0,
140 noisethreshold, # =3.0,
141 lownoisethreshold, # =3.0,
142 negativethreshold, # =0.0,
143 smoothfactor, # =1.0,
144 minbeamfrac, # =0.3,
145 cutthreshold, # =0.01,
146 growiterations, # =100
147 dogrowprune, # =True
148 minpercentchange, # =0.0
149 verbose, # =False
150 fastnoise, # =False
151 ## Misc
152 restart, # =True,
153 savemodel, # ="none",
154 # makeimages,#="auto"
155 calcres, # =True,
156 calcpsf, # =True,
157 psfcutoff, # =0.35
158 ####### State parameters
159 parallel,
160): # =False):
162 #####################################################
163 #### Sanity checks and controls
164 #####################################################
165 ### Move these checks elsewhere ?
166 inpparams=locals().copy()
167# saveinputs(inpparams)
168 ###now deal with parameters which are not the same name
169 inpparams['msname']= inpparams.pop('vis')
170 inpparams['timestr']= inpparams.pop('timerange')
171 inpparams['uvdist']= inpparams.pop('uvrange')
172 inpparams['obs']= inpparams.pop('observation')
173 inpparams['state']= inpparams.pop('intent')
174 inpparams['loopgain']=inpparams.pop('gain')
175 inpparams['scalebias']=inpparams.pop('smallscalebias')
176 #
177 # Force chanchunks=1 always now (CAS-13400)
178 inpparams['chanchunks']=1
180 if specmode=='cont':
181 specmode='mfs'
182 inpparams['specmode']='mfs'
183# if specmode=='mfs' and nterms==1 and deconvolver == "mtmfs":
184# casalog.post( "The MTMFS deconvolution algorithm (deconvolver='mtmfs') needs nterms>1.Please set nterms=2 (or more). ", "WARN", "task_tclean" )
185# return
187 # Force chanchunks=1 always now (CAS-13400)
188 inpparams["chanchunks"] = 1
190 if specmode == "cont":
191 specmode = "mfs"
192 inpparams["specmode"] = "mfs"
193 # if specmode=='mfs' and nterms==1 and deconvolver == "mtmfs":
194 # casalog.post( "The MTMFS deconvolution algorithm (deconvolver='mtmfs') needs nterms>1.Please set nterms=2 (or more). ", "WARN", "task_tclean" )
195 # return
197 if deconvolver == "mtmfs" and (specmode == "cube" or specmode == "cubedata"):
198 casalog.post(
199 "The MSMFS algorithm (deconvolver='mtmfs') with specmode='cube' is not supported",
200 "WARN",
201 "task_tclean",
202 )
203 return
205 if (
206 (specmode == "cube" or specmode == "cubedata" or specmode == "cubesource")
207 ) and (parallel == False and mpi_available and MPIEnvironment.is_mpi_enabled):
208 casalog.post(
209 "When CASA is launched with mpi, the parallel=False option has no effect for 'cube' imaging for gridder='mosaic','wproject','standard' and major cycles are always executed in parallel.\n",
210 "WARN",
211 "task_tclean",
212 )
213 # casalog.post( "Setting parameter parallel=False with specmode='cube' when launching CASA with mpi has no effect except for awproject.", "WARN", "task_tclean" )
216 ## Part of CAS-13814, checking for the only options compatible with mtmfs_via_cube.
217 if specmode=="mvc":
218 if deconvolver != 'mtmfs' or nterms<=1 :
219 casalog.post("The specmode='mvc' option requires the deconvolver to be 'mtmfs' and 'nterms>1.",
220 "WARN",
221 "task_tclean")
222 return
224 ## Part of CAS-13814, moving warnings about pbcor and widebandpbcor from the C++ code to here, for better access to user-settings.
225 if specmode=='mfs' and deconvolver=='mtmfs' and gridder in ['standard','mosaic'] and pbcor==True:
226 casalog.post("For specmode='mfs' and deconvolver='mtmfs', the option of pbcor=True divides each restored Taylor coefficient image by the pb.tt0 image. This correction ignores the frequency-dependence of the primary beam and does not correct for PB spectral index. It is scientifically valid only for small fractional bandwidths. For more accurate wideband primary beam correction (if needed), please use one of the following options : (1) specmode='mvc' with gridder='standard' or 'mosaic' with pbcor=True, (2) conjbeams=True and wbawp=True with gridder='awproject' and pbcor=True.",
227 "WARN",
228 "task_tclean")
229 if perchanweightdensity == False and weighting == "briggsbwtaper":
230 casalog.post(
231 "The briggsbwtaper weighting scheme is not compatable with perchanweightdensity=False.",
232 "WARN",
233 "task_tclean",
234 )
236 if (specmode == "mfs" or specmode == "cont") and weighting == "briggsbwtaper":
237 casalog.post(
238 "The briggsbwtaper weighting scheme is not compatable with specmode='mfs' or 'cont'.",
239 "WARN",
240 "task_tclean",
241 )
242 return
244 if npixels != 0 and weighting == "briggsbwtaper":
245 casalog.post(
246 "The briggsbwtaper weighting scheme is not compatable with npixels != 0.",
247 "WARN",
248 "task_tclean",
249 )
250 return
252 if facets > 1 and parallel == True:
253 casalog.post(
254 "Facetted imaging currently works only in serial. Please choose pure W-projection instead.",
255 "WARN",
256 "task_tclean",
257 )
259 if nmajor < -1:
260 casalog.post(
261 "Negative values less than -1 for nmajor are reserved for possible future implementation",
262 "WARN",
263 "task_tclean",
264 )
265 return
267 ## CAS-13814
268 if (specmode == "mvc" and gridder == 'awproject' and (conjbeams==True or wbawp==False) ):
269 casalog.post(
270 "specmode='mvc' requires frequency-dependent primary beams to be used during cube gridding. Please set conjbeams=False and wbawp=True for the awproject gridder.",
271 "WARN",
272 "task_tclean",
273 )
274 return
276 #CAS-13814
277 if (specmode == "mvc" and gridder == 'mosaic' and conjbeams==True):
278 casalog.post(
279 "specmode='mvc' requires frequency-dependent primary beams to be used during cube gridding. Please set conjbeams=False with the mosaic gridder.",
280 "WARN",
281 "task_tclean",
282 )
283 return
286 #####################################################
287 #### Construct ImagerParameters object
288 #####################################################
290 imager = None
291 paramList = None
293 # Put all parameters into dictionaries and check them.
294 ##make a dictionary of parameters that ImagerParameters take
296 defparm = dict(
297 list(
298 zip(
299 ImagerParameters.__init__.__code__.co_varnames[1:],
300 ImagerParameters.__init__.__defaults__,
301 )
302 )
303 )
304 ###assign values to the ones passed to tclean and if not defined yet in tclean...
305 ###assign them the default value of the constructor
306 bparm = {k: inpparams[k] if k in inpparams else defparm[k] for k in defparm.keys()}
308 ###default mosweight=True is tripping other gridders as they are not
309 ###expecting it to be true
310 if bparm["mosweight"] == True and bparm["gridder"].find("mosaic") == -1:
311 bparm["mosweight"] = False
313 if specmode == "mfs":
314 bparm["perchanweightdensity"] = False
316 # deprecation message
317 if usemask == "auto-thresh" or usemask == "auto-thresh2":
318 casalog.post(
319 usemask
320 + " is deprecated, will be removed in CASA 5.4. It is recommended to use auto-multithresh instead",
321 "WARN",
322 )
324 # paramList.printParameters()
326 if (
327 len(pointingoffsetsigdev) > 0
328 and pointingoffsetsigdev[0] != 0.0
329 and usepointing == True
330 and gridder.count("awproj") > 1
331 ):
332 casalog.post(
333 "pointingoffsetsigdev will be used for pointing corrections with AWProjection",
334 "WARN",
335 )
336 # elif usepointing==True and pointingoffsetsigdev[0] == 0:
337 # casalog.post("pointingoffsetsigdev is set to zero which is an unphysical value, will proceed with the native sky pixel resolution instead". "WARN")
339 ##pcube may still need to be set to True for some combination of ftmachine etc...
340 # =========================================================
341 concattype = ""
342 pcube = False
343 if parallel == True and specmode != "mfs":
344 if specmode != "mfs":
345 pcube = False
346 parallel = False
347 else:
348 pcube = True
349 # =========================================================
350 ####set the children to load c++ libraries and applicator
351 ### make workers ready for c++ based mpicommands
352 cppparallel = False
353 if (
354 mpi_available
355 and MPIEnvironment.is_mpi_enabled
356 and specmode != "mfs"
357 and not pcube
358 ):
359 mint = MPIInterface.MPIInterface()
360 cl = mint.getCluster()
361 cl._cluster.pgc("from casatools import synthesisimager", False)
362 cl._cluster.pgc("si=synthesisimager()", False)
363 cl._cluster.pgc("si.initmpi()", False)
364 cppparallel = True
365 ###ignore chanchunk
366 bparm["chanchunks"] = 1
368 if interactive:
369 # catch non operational case (parallel cube tclean with interative=T)
370 if pcube:
371 casalog.post(
372 "Interactive mode is not currently supported with parallel apwproject cube CLEANing, please restart by setting interactive=F",
373 "WARN",
374 "task_tclean",
375 )
376 return False
378 # Check for casaviewer, if it does not exist flag it up front for macOS
379 # since casaviewer is no longer provided by default with macOS. Returning
380 # False instead of throwing an exception results in:
381 #
382 # RuntimeError: No active exception to reraise
383 #
384 # from tclean run from casashell.
385 try:
386 import casaviewer as __test_casaviewer
387 except:
388 if platform.system( ) == "Darwin":
389 casalog.post(
390 "casaviewer is no longer available for macOS, for more information see: http://go.nrao.edu/casa-viewer-eol Please restart by setting interactive=F",
391 "WARN",
392 "task_tclean",
393 )
394 raise RuntimeError( "casaviewer is no longer available for macOS, for more information see: http://go.nrao.edu/casa-viewer-eol" )
396 #casalog.post('parameters {}'.format(bparm))
397 paramList=ImagerParameters(**bparm)
399 ## Setup Imager objects, for different parallelization schemes.
400 imagerInst = PySynthesisImager
401 if specmode == "mvc":
402 imager = PyMtmfsViaCubeSynthesisImager(params=paramList)
403 imagerInst = PyMtmfsViaCubeSynthesisImager
404 elif parallel == False and pcube == False:
405 imager = PySynthesisImager(params=paramList)
406 imagerInst = PySynthesisImager
407 elif parallel == True:
408 imager = PyParallelContSynthesisImager(params=paramList)
409 imagerInst = PySynthesisImager
410 elif pcube == True:
411 imager = PyParallelCubeSynthesisImager(params=paramList)
412 imagerInst = PyParallelCubeSynthesisImager
413 # virtualconcat type - changed from virtualmove to virtualcopy 2016-07-20
414 # using ia.imageconcat now the name changed to copyvirtual 2019-08-12
415 concattype = "copyvirtual"
416 else:
417 casalog.post("Invalid parallel combination in doClean.", "ERROR")
418 return
420 retrec = {}
422 try:
423 # if (1):
424 # pdb.set_trace()
425 ## Init major cycle elements
426 t0 = time.time()
427 imager.initializeImagers()
429 # Construct the CFCache for AWProject-class of FTMs. For
430 # other choices the following three calls become NoOps.
431 # imager.dryGridding();
432 # imager.fillCFCache();
433 # imager.reloadCFCache();
435 imager.initializeNormalizers()
436 imager.setWeighting()
437 t1 = time.time()
438 casalog.post(
439 "***Time for initializing imager and normalizers: "
440 + "%.2f" % (t1 - t0)
441 + " sec",
442 "INFO3",
443 "task_tclean",
444 )
446 ## Init minor cycle elements
447 if niter > 0 or restoration == True:
448 t0 = time.time()
449 imager.initializeDeconvolvers()
450 t1 = time.time()
451 casalog.post(
452 "***Time for initializing deconvolver(s): "
453 + "%.2f" % (t1 - t0)
454 + " sec",
455 "INFO3",
456 "task_tclean",
457 )
459 ####now is the time to check estimated memory
460 imager.estimatememory()
462 if niter > 0:
463 t0 = time.time()
464 imager.initializeIterationControl()
465 t1 = time.time()
466 casalog.post(
467 "***Time for initializing iteration controller: "
468 + "%.2f" % (t1 - t0)
469 + " sec",
470 "INFO3",
471 "task_tclean",
472 )
474 ## Make PSF
475 if calcpsf == True:
476 t0 = time.time()
478 imager.makePSF()
479 if (psfphasecenter != "") and ("mosaic" in gridder):
480 ###for some reason imager keeps the psf open delete it and recreate it afterwards
481 imager.deleteTools()
482 mytb = table()
483 psfname = (
484 bparm["imagename"] + ".psf.tt0"
485 if (os.path.exists(bparm["imagename"] + ".psf.tt0"))
486 else bparm["imagename"] + ".psf"
487 )
488 mytb.open(psfname)
489 miscinf = mytb.getkeyword("miscinfo")
490 iminf = mytb.getkeyword("imageinfo")
491 # casalog.post ('miscinfo {} {}'.format(miscinf, iminf))
492 mytb.done()
493 casalog.post("doing with different phasecenter psf")
494 imager.unlockimages(0)
495 psfParameters = paramList.getAllPars()
496 psfParameters["phasecenter"] = psfphasecenter
497 psfParamList = ImagerParameters(**psfParameters)
498 psfimager = imagerInst(params=psfParamList)
499 psfimager.initializeImagers()
500 psfimager.setWeighting()
501 psfimager.makeImage("psf", psfParameters["imagename"] + ".psf")
502 psfimager.deleteTools()
503 mytb.open(psfname, nomodify=False)
504 mytb.putkeyword("imageinfo", iminf)
505 mytb.putkeyword("miscinfo", miscinf)
506 mytb.done()
507 mysu=synthesisutils()
508 mysu.fitPsfBeam(imagename=bparm['imagename'],
509 nterms=(bparm['nterms'] if deconvolver=="mtmfs" else 1),
510 psfcutoff=bparm['psfcutoff'])
511 imager = PySynthesisImager(params=paramList)
512 imager.initializeImagers()
513 imager.initializeNormalizers()
514 imager.setWeighting()
515 ###redo these as we destroyed things for lock issues
516 ## Init minor cycle elements
517 if niter > 0 or restoration == True:
518 imager.initializeDeconvolvers()
519 if niter > 0:
520 imager.initializeIterationControl()
522 t1 = time.time()
523 if specmode != "mfs" and ("stand" in gridder):
524 casalog.post(
525 "***Time for making PSF and PB: " + "%.2f" % (t1 - t0) + " sec",
526 "INFO3",
527 "task_tclean",
528 )
529 else:
530 casalog.post(
531 "***Time for making PSF: " + "%.2f" % (t1 - t0) + " sec",
532 "INFO3",
533 "task_tclean",
534 )
536 imager.makePB()
538 t2 = time.time()
539 if specmode == "mfs" and ("stand" in gridder):
540 casalog.post(
541 "***Time for making PB: " + "%.2f" % (t2 - t1) + " sec",
542 "INFO3",
543 "task_tclean",
544 )
546 if gridder in ["mosaic", "awproject"]:
547 imager.checkPB()
549 if niter >= 0:
551 ## Make dirty image
552 if calcres == True:
553 t0 = time.time()
554 imager.runMajorCycle(isCleanCycle=False)
555 t1 = time.time()
556 casalog.post(
557 "***Time for major cycle (calcres=T): "
558 + "%.2f" % (t1 - t0)
559 + " sec",
560 "INFO3",
561 "task_tclean",
562 )
564 ## In case of no deconvolution iterations....
565 if niter == 0 and calcres == False:
566 if savemodel != "none":
567 imager.predictModel()
569 # CAS-13960 : Construct return dict for niter=0 case
570 # If residual image does not exist, summaryminor will not be
571 # populated.
572 if niter==0:
573 id = ImagingDict()
574 retrec = id.construct_residual_dict(paramList)
576 ## Do deconvolution and iterations
577 if niter>0 :
578 t0=time.time();
579 isit = imager.hasConverged()
580 imager.updateMask()
581 # if((type(usemask)==str) and ('auto' in usemask)):
582 # isit = imager.hasConverged()
583 isit = imager.hasConverged()
585 t1 = time.time()
586 casalog.post(
587 "***Time to update mask: " + "%.2f" % (t1 - t0) + " sec",
588 "INFO3",
589 "task_tclean",
590 )
591 while not isit:
592 t0 = time.time()
593 ### sometimes after automasking it does not do anything
594 doneMinor = imager.runMinorCycle()
595 t1 = time.time()
596 casalog.post(
597 "***Time for minor cycle: " + "%.2f" % (t1 - t0) + " sec",
598 "INFO3",
599 "task_tclean",
600 )
602 t0 = time.time()
603 if doneMinor:
604 imager.runMajorCycle()
605 t1 = time.time()
606 casalog.post(
607 "***Time for major cycle: " + "%.2f" % (t1 - t0) + " sec",
608 "INFO3",
609 "task_tclean",
610 )
612 imager.updateMask()
613 t2 = time.time()
614 casalog.post(
615 "***Time to update mask: " + "%.2f" % (t2 - t1) + " sec",
616 "INFO3",
617 "task_tclean",
618 )
619 isit = imager.hasConverged() or (not doneMinor)
621 ## Get summary from iterbot
622 #if type(interactive) != bool:
623 retrec=imager.getSummary(fullsummary);
625 if savemodel!='none' and (interactive==True or usemask=='auto-multithresh' or nsigma>0.0):
626 paramList.resetParameters()
627 if parallel and specmode == "mfs":
628 # For parallel mfs, also needs to reset the parameters for each node
629 imager.resetSaveModelParams(paramList)
630 imager.initializeImagers()
631 imager.predictModel()
633 ## Restore images.
634 if restoration == True:
635 t0 = time.time()
636 imager.restoreImages()
637 t1 = time.time()
638 casalog.post(
639 "***Time for restoring images: " + "%.2f" % (t1 - t0) + " sec",
640 "INFO3",
641 "task_tclean",
642 )
643 if pbcor == True:
644 t0 = time.time()
645 imager.pbcorImages()
646 t1 = time.time()
647 casalog.post(
648 "***Time for pb-correcting images: "
649 + "%.2f" % (t1 - t0)
650 + " sec",
651 "INFO3",
652 "task_tclean",
653 )
654 ######### niter >=0 end if
656 finally:
657 ##close tools
658 # needs to deletools before concat or lock waits for ever
659 if imager != None:
660 imager.deleteTools()
661 if cppparallel:
662 ###release workers back to python mpi control
663 si = synthesisimager()
664 si.releasempi()
665 if pcube:
666 casalog.post("running concatImages ...")
667 casalog.post(
668 "Running virtualconcat (type=%s) of sub-cubes" % concattype,
669 "INFO2",
670 "task_tclean",
671 )
672 imager.concatImages(type=concattype)
673 # CAS-10721
674 # if niter>0 and savemodel != "none":
675 # casalog.post("Please check the casa log file for a message confirming that the model was saved after the last major cycle. If it doesn't exist, please re-run tclean with niter=0,calcres=False,calcpsf=False in order to trigger a 'predict model' step that obeys the savemodel parameter.","WARN","task_tclean")
677 # Write history at the end, when hopefully all .workdirectory, .work.temp, etc. are gone
678 # from disk, so they won't be picked up. They need time to disappear on NFS or slow hw.
679 try:
680 params = get_func_params(tclean, locals())
681 write_tclean_history(imagename, "tclean", params, casalog)
682 except Exception as exc:
683 casalog.post("Error updating history (logtable): {} ".format(exc), "WARN")
685 return retrec
688##################################################