Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/imagerhelpers/_gclean.py: 63%
323 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
1########################################################################3
2# _gclean.py
3#
4# Copyright (C) 2021,2022,2023
5# Associated Universities, Inc. Washington DC, USA.
6#
7# This script is free software; you can redistribute it and/or modify it
8# under the terms of the GNU Library General Public License as published by
9# the Free Software Foundation; either version 2 of the License, or (at your
10# option) any later version.
11#
12# This library is distributed in the hope that it will be useful, but WITHOUT
13# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15# License for more details.
16#
17# You should have received a copy of the GNU Library General Public License
18# along with this library; if not, write to the Free Software Foundation,
19# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20#
21# Correspondence concerning AIPS++ should be adressed as follows:
22# Internet email: casa-feedback@nrao.edu.
23# Postal address: AIPS++ Project Office
24# National Radio Astronomy Observatory
25# 520 Edgemont Road
26# Charlottesville, VA 22903-2475 USA
27#
28import os
29import asyncio
30from functools import reduce
31import copy
32import numpy as np
33import shutil
34import time
35import subprocess
37from casatasks.private.imagerhelpers.imager_return_dict import ImagingDict
38from casatasks import deconvolve, tclean, imstat
40###
41### import check versions
42###
43_GCV001 = True
44_GCV002 = True
45_GCV003 = True
46_GCV004 = True
49# from casatasks.private.imagerhelpers._gclean import gclean
50class gclean:
51 '''gclean(...) creates a stream of convergence records which indicate
52 the convergence quaility of the tclean process. The initial record
53 describes the initial dirty image.
54 It is designed for use with the interactive clean GUI, but it could
55 be used independently. It can be used as a regular generator:
56 for rec in gclean( vis='refim_point_withline.ms', imagename='test', imsize=512, cell='12.0arcsec',
57 specmode='cube', interpolation='nearest', nchan=5, start='1.0GHz', width='0.2GHz',
58 pblimit=-1e-05, deconvolver='hogbom', niter=500, cyclefactor=3, scales=[0, 3, 10] ):
59 # use rec to decide when to stop, for example to check stopcode or peak residual:
60 # if (rec[0] > 1) or (min(rec[1][0][0]['peakRes']) < 0.001):
61 # break
62 print(rec)
63 or as an async generator:
64 async for rec in gclean( vis='refim_point_withline.ms', imagename='test', imsize=512, cell='12.0arcsec',
65 specmode='cube', interpolation='nearest', nchan=5, start='1.0GHz', width='0.2GHz',
66 pblimit=-1e-05, deconvolver='hogbom', niter=500, cyclefactor=3, scales=[0, 3, 10] ):
67 # use rec to decide when to stop
68 print(rec)
71 See also: __next__(...) for a description of the returned rec
73 TODO: do we need to preserve any hidden state between tclean calls for the iterbotsink and/or synthesisimager tools?
74 '''
76 def _tclean( self, *args, **kwargs ):
77 """ Calls tclean records the arguments in the local history of tclean calls.
79 The full tclean history for this instance can be retrieved via the cmds() method."""
80 arg_s = ', '.join( map( lambda a: self._history_filter(len(self._exe_cmds), None, repr(a)), args ) )
81 kw_s = ', '.join( map( lambda kv: self._history_filter(len(self._exe_cmds), kv[0], "%s=%s" % (kv[0],repr(kv[1]))), kwargs.items()) )
82 if len(arg_s) > 0 and len(kw_s) > 0:
83 parameters = arg_s + ", " + kw_s
84 else:
85 parameters = arg_s + kw_s
86 self._exe_cmds.append( "tclean( %s )" % parameters )
87 self._exe_cmds_per_iter[-1] += 1
88 return tclean( *args, **kwargs )
90 def _deconvolve( self, *args, **kwargs ):
91 """ Calls deconvolve records the arguments in the local history of deconvolve calls.
93 The full deconvolve history for this instance can be retrieved via the cmds() method."""
94 arg_s = ', '.join( map( lambda a: self._history_filter(len(self._exe_cmds), None, repr(a)), args ) )
95 kw_s = ', '.join( map( lambda kv: self._history_filter(len(self._exe_cmds), kv[0], "%s=%s" % (kv[0],repr(kv[1]))), kwargs.items()) )
96 if len(arg_s) > 0 and len(kw_s) > 0:
97 parameters = arg_s + ", " + kw_s
98 else:
99 parameters = arg_s + kw_s
100 self._exe_cmds.append( "deconvolve( %s )" % parameters )
101 self._exe_cmds_per_iter[-1] += 1
102 return deconvolve( *args, **kwargs )
104 def _remove_tree( self, directory ):
105 if os.path.isdir(directory):
106 shutil.rmtree(directory)
107 self._exe_cmds.append( f'''shutil.rmtree( {repr(directory)} )''' )
108 self._exe_cmds_per_iter[-1] += 1
110 def cmds( self, history=False ):
111 """ Returns the history of all tclean calls for this instance. If ``history``
112 is set to True then the full history will be returned, otherwise the commands
113 executed for generating the latest result are returned.
114 """
116 if history:
117 return self._exe_cmds
118 else:
119 if self._exe_cmds_per_iter[-1] > 0:
120 # Return the last N commands
121 return self._exe_cmds[-self._exe_cmds_per_iter[-1]:]
122 else:
123 # If convergence is hit, no commands were run so return nothing
124 return []
127 def update( self, msg ):
128 """ Interactive clean parameters update.
130 Args:
131 msg: dict with possible keys 'niter', 'cycleniter', 'nmajor', 'threshold', 'cyclefactor'
133 Returns:
134 stopcode : Stop code in case of error (-1 on error, 0 if no error), int
135 stopdesc : Exception error message, str
136 """
137 if 'niter' in msg:
138 try:
139 self._niter = int(msg['niter'])
140 if self._niter < -1:
141 return -1, f"niter must be >= -1"
142 except ValueError as err:
143 return -1, "niter must be an integer"
145 if 'cycleniter' in msg:
146 try:
147 self._cycleniter = int(msg['cycleniter'])
148 if self._cycleniter < -1:
149 return -1, f"cycleniter must be >= -1"
150 except ValueError:
151 return -1, "cycleniter must be an integer"
153 if 'nmajor' in msg:
154 try:
155 self._nmajor = int(msg['nmajor'])
156 if self._nmajor < -1:
157 return -1, f"nmajor must be >= -1"
158 except ValueError as e:
159 return -1, "nmajor must be an integer"
161 if 'threshold' in msg:
162 try:
163 self._threshold = float(msg['threshold'])
164 if self._threshold < 0:
165 return -1, f"threshold must be >= 0"
166 except ValueError:
167 if isinstance(msg['threshold'], str) and "jy" in msg['threshold'].lower():
168 self._threshold_to_float(msg['threshold']) # Convert str to float
169 else:
170 return -1, f"threshold must be a number, or a number with units (Jy/mJy/uJy)"
172 if 'cyclefactor' in msg:
173 try:
174 self._cyclefactor = float(msg['cyclefactor'])
175 if self._cyclefactor <= 0:
176 return -1, f"cyclefactor must be > 0"
177 except ValueError:
178 return -1, "cyclefactor must be a number"
180 return 0, ""
184 def _threshold_to_float(self, msg=None):
185 # Convert threshold from string to float if necessary
186 if msg is not None:
187 if isinstance(msg, str):
188 if "mJy" in msg:
189 self._threshold = float(msg.replace("mJy", "")) / 1e3
190 elif "uJy" in msg:
191 self._threshold = float(msg.replace("uJy", "")) / 1e6
192 elif "Jy" in msg:
193 self._threshold = float(msg.replace("Jy", ""))
194 else:
195 if isinstance(self._threshold, str):
196 if "mJy" in self._threshold:
197 self._threshold = float(self._threshold.replace("mJy", "")) / 1e3
198 elif "uJy" in self._threshold:
199 self._threshold = float(self._threshold.replace("uJy", "")) / 1e6
200 elif "Jy" in self._threshold:
201 self._threshold = float(self._threshold.replace("Jy", ""))
204 def __init__( self, vis, imagename, field='', spw='', timerange='', uvrange='', antenna='', scan='', observation='', intent='', datacolumn='corrected', imsize=[100], cell=[ ],
205 phasecenter='', stokes='I', startmodel='', specmode='cube', reffreq='', nchan=-1, start='', width='', outframe='LSRK', veltype='radio', restfreq='', interpolation='linear',
206 perchanweightdensity=True, gridder='standard', wprojplanes=int(1), mosweight=True, psterm=False, wbawp=True, conjbeams=False, usepointing=False, pointingoffsetsigdev=[ ],
207 pblimit=0.2, deconvolver='hogbom', smallscalebias=0.0, niter=0, threshold='0.1Jy', nsigma=0.0, cycleniter=-1, nmajor=1, cyclefactor=1.0, minpsffraction=0.05,
208 maxpsffraction=0.8, scales=[], restoringbeam='', pbcor=False, nterms=int(2), weighting='natural', robust=float(0.5), npixels=0, gain=float(0.1), pbmask=0.2, sidelobethreshold=3.0,
209 noisethreshold=5.0, lownoisethreshold=1.5, negativethreshold=0.0, smoothfactor=float(1.0), minbeamfrac=0.3, cutthreshold=0.01, growiterations=75, dogrowprune=True,
210 minpercentchange=-1.0, verbose=False, fastnoise=True, savemodel='none', usemask='user', mask='', parallel=False, history_filter=lambda index, arg, history_value: history_value ):
212 self._vis = vis
213 self._imagename = imagename
214 self._imsize = imsize
215 self._cell = cell
216 self._phasecenter = phasecenter
217 self._stokes = stokes
218 self._startmodel = startmodel
219 self._specmode = specmode
220 self._reffreq = reffreq
221 self._nchan = nchan
222 self._start = start
223 self._width = width
224 self._outframe = outframe
225 self._veltype = veltype
226 self._restfreq = restfreq
227 self._interpolation = interpolation
228 self._perchanweightdensity = perchanweightdensity
229 self._gridder = gridder
230 self._wprojplanes = wprojplanes
231 self._mosweight = mosweight
232 self._psterm = psterm
233 self._wbawp = wbawp
234 self._conjbeams = conjbeams
235 self._usepointing = usepointing
236 self._pointingoffsetsigdev = pointingoffsetsigdev
237 self._pblimit = pblimit
238 self._deconvolver = deconvolver
239 self._smallscalebias = smallscalebias
240 self._niter = niter
241 self._threshold = threshold
242 self._cycleniter = cycleniter
243 self._minpsffraction = minpsffraction
244 self._maxpsffraction = maxpsffraction
245 self._nsigma = nsigma
246 self._nmajor = nmajor
247 self._cyclefactor = cyclefactor
248 self._scales = scales
249 self._restoringbeam = restoringbeam
250 self._pbcor = pbcor
251 self._nterms = nterms
252 self._exe_cmds = [ ]
253 self._exe_cmds_per_iter = [ ]
254 self._history_filter = history_filter
255 self._finalized = False
256 self._field = field
257 self._spw = spw
258 self._timerange = timerange
259 self._uvrange = uvrange
260 self._antenna = antenna
261 self._scan = scan
262 self._observation = observation
263 self._intent = intent
264 self._datacolumn = datacolumn
265 self._weighting = weighting
266 self._robust = robust
267 self._npixels = npixels
268 self._gain = gain
269 self._pbmask = pbmask
270 self._sidelobethreshold = sidelobethreshold
271 self._noisethreshold = noisethreshold
272 self._lownoisethreshold = lownoisethreshold
273 self._negativethreshold = negativethreshold
274 self._smoothfactor = smoothfactor
275 self._minbeamfrac = minbeamfrac
276 self._cutthreshold = cutthreshold
277 self._growiterations = growiterations
278 self._dogrowprune = dogrowprune
279 self._minpercentchange = minpercentchange
280 self._verbose = verbose
281 self._fastnoise = fastnoise
282 self._savemodel = savemodel
283 self._parallel = parallel
284 self._usemask = usemask
285 self._mask = mask
286 self.global_imdict = ImagingDict()
287 self.current_imdict = ImagingDict()
288 self._major_done = 0
289 self.hasit = False # Convergence flag
290 self.stopdescription = '' # Convergence flag
291 self._initial_mask_exists = False
292 self._convergence_result = (None,None,None,None,None,{ 'chan': None, 'major': None })
293 # ^^^^ ^^^^ ^^^^ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^----->>> convergence info
294 # | | | | +----->>> Number of global iterations remaining for current run (niterleft)
295 # | | | +---------->>> Number of major cycles remaining for current run (nmajorleft)
296 # | | +---------------->>> major cycles done for current run (nmajordone)
297 # | +------------------>>> tclean stopcode
298 # +----------------------->>> tclean stopdescription
300 # Convert threshold from string to float, interpreting units.
301 # XXX : We should ideally use quantities, but we are trying to
302 # stick to "public API" funtions inside _gclean
303 self._threshold_to_float()
306 def __add_per_major_items( self, tclean_ret, major_ret, chan_ret ):
307 '''Add meta-data about the whole major cycle, including 'cyclethreshold'
308 '''
310 if 'cyclethreshold' in tclean_ret:
312 rdict = dict( major=dict( cyclethreshold=[tclean_ret['cyclethreshold']] if major_ret is None else (major_ret['cyclethreshold'] + [tclean_ret['cyclethreshold']]) ),
313 chan=chan_ret )
314 else:
315 rdict = dict( major=dict( cyclethreshold=major_ret['cyclethreshold'].append(tclean_ret['cyclethreshold']) ),
316 chan=chan_ret )
318 return rdict
321 def _calc_deconv_controls(self, imdict, niterleft=0, threshold=0, cycleniter=-1):
322 """
323 Calculate cycleniter and cyclethreshold for deconvolution.
324 """
326 use_cycleniter = niterleft #niter - imdict.returndict['iterdone']
328 if cycleniter > -1 : # User is forcing this number
329 use_cycleniter = min(cycleniter, use_cycleniter)
331 psffrac = imdict.returndict['maxpsfsidelobe'] * self._cyclefactor
332 psffrac = max(psffrac, self._minpsffraction)
333 psffrac = min(psffrac, self._maxpsffraction)
335 # TODO : This assumes the default field (i.e., field=0);
336 # This won't work for multiple fields.
337 cyclethreshold = psffrac * imdict.get_peakres()
338 cyclethreshold = max(cyclethreshold, threshold)
340 return int(use_cycleniter), cyclethreshold
343 def __update_convergence(self):
344 """
345 Accumulates the per-channel/stokes summaryminor keys across all major cycle calls so far.
347 The "iterDone" key will be replaced with "iterations", and for the "iterations" key,
348 the value in the returned cummulative record will be a rolling sum of iterations done
349 for tclean calls so far, one value per minor cycle.
350 For example, if there have been two clean calls, and in the first call channel 0 had
351 [1] iteration in 1 minor cycle, and for the second call channel 0 had [6, 10, 9, 1]
352 iterations in 4 minor cycles), then the resultant "iterations" key for channel 0 would be:
353 [1, 7, 17, 26, 27]
354 """
356 keys = ['modelFlux', 'iterDone', 'peakRes', 'stopCode', 'cycleThresh']
358 # Grab tuples of keys of interest
359 outrec = {}
360 for nn in range(self.global_imdict.nchan):
361 outrec[nn] = {}
362 for ss in range(self.global_imdict.nstokes):
363 outrec[nn][ss] = {}
364 for key in keys:
365 # Replace iterDone with iterations
366 if key == 'iterDone':
367 # Maintain cumulative sum of iterations per entry
368 outrec[nn][ss]['iterations'] = np.cumsum(self.global_imdict.get_key(key, stokes=ss, chan=nn))
369 # Replace iterDone with iterations
370 #outrec[nn][ss]['iterations'] = self.global_imdict.get_key(key, stokes=ss, chan=nn)
371 else:
372 outrec[nn][ss][key] = self.global_imdict.get_key(key, stokes=ss, chan=nn)
374 return outrec
377 def _check_initial_mask(self):
378 """
379 Check if a mask from a previous run exists on disk or not.
380 """
382 if self._usemask == 'user' and self._mask == '':
383 maskname = self._imagename + '.mask'
385 if os.path.exists(maskname):
386 self._initial_mask_exists = True
387 else:
388 self._initial_mask_exists = False
390 def _fix_initial_mask(self):
391 """
392 If on start up, no user mask is provided, then flip the initial mask to
393 be all zeros for interactive use.
394 """
396 from casatools import image
397 ia = image()
399 if self._usemask == 'user' and self._mask == '':
400 maskname = self._imagename + '.mask'
402 # This means the mask was newly created by deconvolve, so flip it
403 if os.path.exists(maskname) and self._initial_mask_exists is False:
404 ia.open(maskname)
405 ia.set(0.0)
406 ia.close()
408 def _update_peakres(self):
409 if self._deconvolver == 'mtmfs':
410 residname = self._imagename + '.residual.tt0'
411 else:
412 residname = self._imagename + '.residual'
414 maskname = self._imagename + '.mask'
415 if not os.path.exists(maskname):
416 maskname = ''
418 peakres = imstat(imagename=residname, mask=maskname)['max']
419 if len(maskname) > 0:
420 masksum = imstat(imagename=maskname)['sum']
421 else:
422 masksum = []
424 if len(peakres) > 0:
425 peakres = peakres[0]
426 else:
427 peakres = None
429 if len(masksum) > 0:
430 masksum = masksum[0]
431 else:
432 masksum = None
434 return peakres, masksum
436 def __next__( self ):
437 """ Runs tclean and returns the (stopcode, convergence result) when executed with the python builtin next() function.
439 The returned convergence result is a nested dictionary:
440 {
441 channel id: {
442 stokes id: {
443 summary key: [values, one per minor cycle]
444 },
445 },
446 }
448 See also: gclean.__update_convergence(...)
449 """
451 # vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv------------>>> is done to produce an initial dirty image
452 if self._niter < 1 and self._convergence_result[2] is not None:
453 self._convergence_result = ( f'nothing to run, niter == {self._niter}',
454 self._convergence_result[1],
455 self._major_done,
456 self._nmajor,
457 self._niter,
458 self._convergence_result[5] )
459 return self._convergence_result
460 else:
461 ### CALL SEQUENCE:
462 ### tclean(niter=0),deconvolve(niter=0),tclean(niter=100),deconvolve(niter=0),tclean(niter=100),tclean(niter=0,restoration=True)
463 self._exe_cmds_per_iter.append(0)
464 try:
465 if self._convergence_result[1] is None:
467 # initial call to tclean(...) creates the initial dirty image with niter=0
468 tclean_ret = self._tclean( vis=self._vis, mask=self._mask, imagename=self._imagename, imsize=self._imsize, cell=self._cell,
469 phasecenter=self._phasecenter, stokes=self._stokes, startmodel=self._startmodel, specmode=self._specmode,
470 reffreq=self._reffreq, gridder=self._gridder, wprojplanes=self._wprojplanes, mosweight=self._mosweight,
471 psterm=self._psterm, wbawp=self._wbawp, conjbeams=self._conjbeams, usepointing=self._usepointing,
472 interpolation=self._interpolation, perchanweightdensity=self._perchanweightdensity,
473 nchan=self._nchan, start=self._start, width=self._width, veltype=self._veltype, restfreq=self._restfreq,
474 outframe=self._outframe, pointingoffsetsigdev=self._pointingoffsetsigdev, pblimit=self._pblimit,
475 deconvolver=self._deconvolver, smallscalebias=self._smallscalebias, cyclefactor=self._cyclefactor,
476 scales=self._scales, restoringbeam=self._restoringbeam, pbcor=self._pbcor, nterms=self._nterms,
477 field=self._field, spw=self._spw, timerange=self._timerange, uvrange=self._uvrange, antenna=self._antenna,
478 scan=self._scan, observation=self._observation, intent=self._intent, datacolumn=self._datacolumn,
479 weighting=self._weighting, robust=self._robust, npixels=self._npixels, interactive=False, niter=0,
480 gain=self._gain, calcres=True, calcpsf=True, restoration=False, parallel=self._parallel, fullsummary=True)
482 # Check if a mask from a previous run exists on disk
483 self._check_initial_mask()
485 deconv_ret = self._deconvolve(imagename=self._imagename, startmodel=self._startmodel,
486 deconvolver=self._deconvolver, restoration=False,
487 threshold=self._threshold, niter=0,
488 nsigma=self._nsigma, fullsummary=True, fastnoise=self._fastnoise, usemask=self._usemask,
489 mask=self._mask, pbmask=self._pbmask, sidelobethreshold=self._sidelobethreshold, noisethreshold=self._noisethreshold,
490 lownoisethreshold=self._lownoisethreshold, negativethreshold=self._negativethreshold, smoothfactor=self._smoothfactor,
491 minbeamfrac=self._minbeamfrac, cutthreshold=self._cutthreshold, growiterations=self._growiterations,
492 dogrowprune=self._dogrowprune, verbose=self._verbose)
494 # If no mask from a previous run exists, over-write the ones with zeros for the default mask
495 self._fix_initial_mask()
497 self.current_imdict.returndict = self.current_imdict.merge(tclean_ret, deconv_ret)
498 self.global_imdict.returndict = self.current_imdict.returndict
500 ## Initial call where niterleft and nmajorleft are same as original input values.
501 self.hasit, self.stopdescription = self.global_imdict.has_converged(self._niter, self._threshold, self._nmajor)
503 self.current_imdict.returndict['stopcode'] = self.hasit
504 self.current_imdict.returndict['stopDescription'] = self.stopdescription
505 self._major_done = 0
506 else:
507 # Reset convergence every time, since we return control to the GUI after a single major cycle
508 self.current_imdict.returndict['iterdone'] = 0.
510 # Mask can be updated here...
511 # Check for mask update - peakres + masksum
512 _peakres, _masksum = self._update_peakres()
514 self.hasit, self.stopdescription = self.global_imdict.has_converged(self._niter, self._threshold, self._nmajor, masksum=_masksum, peakres=_peakres)
516 #self.global_imdict.returndict['stopcode'] = self.hasit
517 #self.global_imdict.returndict['stopDescription'] = self.stopdescription
519 #self.current_imdict.returndict['stopcode'] = self.hasit
520 #self.current_imdict.returndict['stopDescription'] = self.stopdescription
522 # Has not, i.e., not converged
523 if self.hasit ==0 :
524 use_cycleniter, cyclethreshold = self._calc_deconv_controls(self.current_imdict, self._niter, self._threshold, self._cycleniter)
526 # Run the minor cycle
527 deconv_ret = self._deconvolve(imagename=self._imagename, startmodel=self._startmodel,
528 deconvolver=self._deconvolver, restoration=False,
529 threshold=cyclethreshold, niter=use_cycleniter, gain=self._gain, fullsummary=True)
531 # Run the major cycle
532 tclean_ret = self._tclean( vis=self._vis, imagename=self._imagename, imsize=self._imsize, cell=self._cell,
533 phasecenter=self._phasecenter, stokes=self._stokes, specmode=self._specmode, reffreq=self._reffreq,
534 gridder=self._gridder, wprojplanes=self._wprojplanes, mosweight=self._mosweight, psterm=self._psterm,
535 wbawp=self._wbawp, conjbeams=self._conjbeams, usepointing=self._usepointing, interpolation=self._interpolation,
536 perchanweightdensity=self._perchanweightdensity, nchan=self._nchan, start=self._start,
537 width=self._width, veltype=self._veltype, restfreq=self._restfreq, outframe=self._outframe,
538 pointingoffsetsigdev=self._pointingoffsetsigdev, pblimit=self._pblimit, deconvolver=self._deconvolver,
539 smallscalebias=self._smallscalebias, cyclefactor=self._cyclefactor, scales=self._scales,
540 restoringbeam=self._restoringbeam, pbcor=self._pbcor, nterms=self._nterms, field=self._field,
541 spw=self._spw, timerange=self._timerange, uvrange=self._uvrange, antenna=self._antenna,
542 scan=self._scan, observation=self._observation, intent=self._intent, datacolumn=self._datacolumn,
543 weighting=self._weighting, robust=self._robust, npixels=self._npixels, interactive=False,
544 niter=0, restart=True, calcpsf=False, calcres=True, restoration=False, threshold=self._threshold,
545 nsigma=self._nsigma, cycleniter=self._cycleniter, nmajor=1, gain=self._gain,
546 sidelobethreshold=self._sidelobethreshold, noisethreshold=self._noisethreshold,
547 lownoisethreshold=self._lownoisethreshold, negativethreshold=self._negativethreshold,
548 minbeamfrac=self._minbeamfrac, growiterations=self._growiterations, dogrowprune=self._dogrowprune,
549 minpercentchange=self._minpercentchange, fastnoise=self._fastnoise, savemodel=self._savemodel,
550 maxpsffraction=self._maxpsffraction,
551 minpsffraction=self._minpsffraction, parallel=self._parallel, fullsummary=True )
553 # Replace return dict with new return dict
554 # The order of the dicts into merge is important.
555 self.current_imdict.returndict = self.current_imdict.merge(tclean_ret, deconv_ret)
557 # Append new return dict to global return dict
558 self.global_imdict.returndict = self.global_imdict.concat(self.global_imdict.returndict, self.current_imdict.returndict)
559 self._major_done = self.current_imdict.returndict['nmajordone']
561 ## Decrement count for the major cycle just done...
562 self.__decrement_counts()
564 # Use global imdict for convergence check
565 if deconv_ret['stopcode'] == 7: ## Tell the convergence checker that the mask is zero and iterations were skipped
566 self.hasit, self.stopdescription = self.global_imdict.has_converged(self._niter, self._threshold, self._nmajor, masksum=0)
567 else:
568 self.hasit, self.stopdescription = self.global_imdict.has_converged(self._niter, self._threshold, self._nmajor)
571 self.global_imdict.returndict['stopcode'] = self.hasit
572 self.global_imdict.returndict['stopDescription'] = self.stopdescription
574 if not self.hasit:
575 # If we haven't converged, run deconvolve to update the mask
576 self._deconvolve(imagename=self._imagename, startmodel=self._startmodel, deconvolver=self._deconvolver, restoration=False, threshold=self._threshold, niter=0,
577 nsigma=self._nsigma, fullsummary=True, fastnoise=self._fastnoise, usemask=self._usemask, mask='', pbmask=self._pbmask,
578 sidelobethreshold=self._sidelobethreshold, noisethreshold=self._noisethreshold, lownoisethreshold=self._lownoisethreshold,
579 negativethreshold=self._negativethreshold, smoothfactor=self._smoothfactor, minbeamfrac=self._minbeamfrac, cutthreshold=self._cutthreshold,
580 growiterations=self._growiterations, dogrowprune=self._dogrowprune, verbose=self._verbose)
583 if len(self.global_imdict.returndict) > 0 and 'summaryminor' in self.global_imdict.returndict and sum(map(len,self.global_imdict.returndict['summaryminor'].values())) > 0:
584 self._convergence_result = ( self.global_imdict.returndict['stopDescription'] if 'stopDescription' in self.global_imdict.returndict else '',
585 self.global_imdict.returndict['stopcode'] if 'stopcode' in self.global_imdict.returndict else 0,
586 self._major_done,
587 self._nmajor,
588 self._niter,
589 self.__add_per_major_items( self.global_imdict.returndict,
590 self._convergence_result[5]['major'],
591 self.__update_convergence()))
592 else:
593 self._convergence_result = ( f'tclean returned an empty result',
594 self._convergence_result[1],
595 self._major_done,
596 self._nmajor,
597 self._niter,
598 self._convergence_result[5] )
599 except Exception as e:
600 self._convergence_result = ( str(e),
601 -1,
602 self._major_done,
603 self._nmajor,
604 self._niter,
605 self._convergence_result[5] )
606 return self._convergence_result
608 return self._convergence_result
610 def __decrement_counts( self ):
611 ## Update niterleft and nmajorleft now.
612 if self.hasit == 0: ##If not yet converged.
613 if self._nmajor != -1: ## If -1, don't touch it.
614 self._nmajor = self._nmajor - 1
615 if self._nmajor<0: ## Force a floor
616 self._nmajor=0
617 self._niter = self._niter - self.current_imdict.get_key('iterdone')
618 if self._niter<0: ## This can happen when we're counting niter across channels in a single minor cycle set, and it crosses the total.
619 self._niter=0 ## Force a floor
620 else:
621 return ##If convergence has been reached, don't try to decrement further.
624 def __reflect_stop( self ):
625 ## if python wasn't hacky, you would be able to try/except/raise in lambda
626 time.sleep(1)
627 try:
628 return self.__next__( )
629 except StopIteration:
630 raise StopAsyncIteration
632 async def __anext__( self ):
633 ### asyncio.run cannot be used here because this is called
634 ### within an asyncio loop...
635 loop = asyncio.get_event_loop( )
636 result = await loop.run_in_executor( None, self.__reflect_stop )
637 return result
639 def __iter__( self ):
640 return self
642 def __aiter__( self ):
643 return self
645 def __split_filename( self, path ):
646 return os.path.splitext(os.path.basename(path))
648 def __default_mask_name( self ):
649 imgparts = self.__split_filename( self._imagename )
650 return f'{imgparts[0]}.mask'
652 def mask(self):
653 #return self.__default_mask_name() if self._mask == '' else self._mask
654 return f'{self._imagename}.mask' if self._mask == '' else self._mask
656 def reset(self):
657 #if not self._finalized:
658 # raise RuntimeError('attempt to reset a gclean run that has not been finalized')
659 self._finalized = False
660 self._convergence_result = ( None,
661 self._convergence_result[1],
662 self._major_done,
663 self._nmajor,
664 self._niter,
665 self._convergence_result[5] )
667 def restore(self):
668 """ Restores the final image, and returns a path to the restored image. """
669 deconv_ret = self._deconvolve(imagename=self._imagename,
670 deconvolver=self._deconvolver,
671 restoration=True, niter=0,
672 fullsummary=True)
674 return { "image": f"{self._imagename}.image" }
676 def has_next(self):
677 return not self._finalized