Commits
David Mehringer authored e71238a4ce8 Merge
1 1 | ########################################################################3 |
2 2 | # _gclean.py |
3 3 | # |
4 - | # Copyright (C) 2021,2022 |
4 + | # Copyright (C) 2021,2022,2023 |
5 5 | # Associated Universities, Inc. Washington DC, USA. |
6 6 | # |
7 7 | # This script is free software; you can redistribute it and/or modify it |
8 8 | # under the terms of the GNU Library General Public License as published by |
9 9 | # the Free Software Foundation; either version 2 of the License, or (at your |
10 10 | # option) any later version. |
11 11 | # |
12 12 | # This library is distributed in the hope that it will be useful, but WITHOUT |
13 13 | # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
14 14 | # FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public |
22 22 | # Internet email: aips2-request@nrao.edu. |
23 23 | # Postal address: AIPS++ Project Office |
24 24 | # National Radio Astronomy Observatory |
25 25 | # 520 Edgemont Road |
26 26 | # Charlottesville, VA 22903-2475 USA |
27 27 | # |
28 28 | import os |
29 29 | import asyncio |
30 30 | from functools import reduce |
31 31 | import copy |
32 + | import numpy as np |
33 + | import shutil |
34 + | import time |
35 + | import subprocess |
32 36 | |
37 + | from casatasks.private.imagerhelpers.imager_return_dict import ImagingDict |
38 + | from casatasks import deconvolve, tclean, imstat |
39 + | |
40 + | ### |
41 + | ### import check versions |
42 + | ### |
33 43 | _GCV001 = True |
34 44 | _GCV002 = True |
35 45 | _GCV003 = True |
46 + | _GCV004 = True |
47 + | |
36 48 | |
37 49 | # from casatasks.private.imagerhelpers._gclean import gclean |
38 50 | class gclean: |
39 51 | '''gclean(...) creates a stream of convergence records which indicate |
40 52 | the convergence quaility of the tclean process. The initial record |
41 53 | describes the initial dirty image. |
42 54 | It is designed for use with the interactive clean GUI, but it could |
43 55 | be used independently. It can be used as a regular generator: |
44 56 | for rec in gclean( vis='refim_point_withline.ms', imagename='test', imsize=512, cell='12.0arcsec', |
45 57 | specmode='cube', interpolation='nearest', nchan=5, start='1.0GHz', width='0.2GHz', |
58 70 | |
59 71 | See also: __next__(...) for a description of the returned rec |
60 72 | |
61 73 | TODO: do we need to preserve any hidden state between tclean calls for the iterbotsink and/or synthesisimager tools? |
62 74 | ''' |
63 75 | |
64 76 | def _tclean( self, *args, **kwargs ): |
65 77 | """ Calls tclean records the arguments in the local history of tclean calls. |
66 78 | |
67 79 | The full tclean history for this instance can be retrieved via the cmds() method.""" |
68 - | from casatasks import tclean |
69 80 | arg_s = ', '.join( map( lambda a: self._history_filter(len(self._exe_cmds), None, repr(a)), args ) ) |
70 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()) ) |
71 82 | if len(arg_s) > 0 and len(kw_s) > 0: |
72 83 | parameters = arg_s + ", " + kw_s |
73 84 | else: |
74 85 | parameters = arg_s + kw_s |
75 86 | self._exe_cmds.append( "tclean( %s )" % parameters ) |
87 + | self._exe_cmds_per_iter[-1] += 1 |
76 88 | return tclean( *args, **kwargs ) |
77 89 | |
78 90 | def _deconvolve( self, *args, **kwargs ): |
79 91 | """ Calls deconvolve records the arguments in the local history of deconvolve calls. |
80 92 | |
81 93 | The full deconvolve history for this instance can be retrieved via the cmds() method.""" |
82 - | from casatasks import deconvolve |
83 94 | arg_s = ', '.join( map( lambda a: self._history_filter(len(self._exe_cmds), None, repr(a)), args ) ) |
84 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()) ) |
85 96 | if len(arg_s) > 0 and len(kw_s) > 0: |
86 97 | parameters = arg_s + ", " + kw_s |
87 98 | else: |
88 99 | parameters = arg_s + kw_s |
89 100 | self._exe_cmds.append( "deconvolve( %s )" % parameters ) |
101 + | self._exe_cmds_per_iter[-1] += 1 |
90 102 | return deconvolve( *args, **kwargs ) |
91 103 | |
92 - | def cmds( self ): |
93 - | """ Returns the history of all tclean calls for this instance. """ |
94 - | return self._exe_cmds |
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 |
109 + | |
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 + | """ |
115 + | |
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 [] |
125 + | |
95 126 | |
96 127 | def update( self, msg ): |
97 128 | """ Interactive clean parameters update. |
98 129 | |
99 130 | Args: |
100 - | msg: dict with possible keys 'niter', 'cycleniter', 'nmajor', 'threshold', 'cyclefactor' and 'mask' |
131 + | msg: dict with possible keys 'niter', 'cycleniter', 'nmajor', 'threshold', 'cyclefactor' |
132 + | |
133 + | Returns: |
134 + | stopcode : Stop code in case of error (-1 on error, 0 if no error), int |
135 + | stopdesc : Exception error message, str |
101 136 | """ |
102 137 | if 'niter' in msg: |
103 138 | try: |
104 139 | self._niter = int(msg['niter']) |
105 - | except ValueError: |
106 - | pass |
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" |
144 + | |
107 145 | if 'cycleniter' in msg: |
108 146 | try: |
109 147 | self._cycleniter = int(msg['cycleniter']) |
148 + | if self._cycleniter < -1: |
149 + | return -1, f"cycleniter must be >= -1" |
110 150 | except ValueError: |
111 - | pass |
151 + | return -1, "cycleniter must be an integer" |
152 + | |
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" |
160 + | |
112 161 | if 'threshold' in msg: |
113 - | self._threshold = msg['threshold'] |
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)" |
171 + | |
114 172 | if 'cyclefactor' in msg: |
115 173 | try: |
116 - | self._cyclefactor = int(msg['cyclefactor']) |
174 + | self._cyclefactor = float(msg['cyclefactor']) |
175 + | if self._cyclefactor <= 0: |
176 + | return -1, f"cyclefactor must be > 0" |
117 177 | except ValueError: |
118 - | pass |
119 - | if 'mask' in msg: |
120 - | self._mask = msg['mask'] |
121 - | |
122 - | def __init__( self, vis, imagename, field='', spw='', timerange='', uvrange='', antenna='', scan='', observation='', intent='', datacolumn='corrected', |
123 - | imsize=[100], cell=[ ], phasecenter='', stokes='I', startmodel='', specmode='cube', reffreq='', nchan=-1, start='', width='', outframe='LSRK', |
124 - | restfreq='', interpolation='linear', perchanweightdensity=True, gridder='standard', wprojplanes=int(1), mosweight=True, psterm=False, wbawp=True, |
125 - | conjbeams=False, usepointing=False, pointingoffsetsigdev=[ ], pblimit=0.2, deconvolver='hogbom', niter=0, threshold='0.1Jy', nsigma=0.0, |
126 - | cycleniter=-1, nmajor=1, cyclefactor=1.0, scales=[], restoringbeam='', pbcor=False, nterms=int(2), weighting='natural', robust=float(0.5), |
127 - | npixels=0, gain=float(0.1), sidelobethreshold=3.0, noisethreshold=5.0, lownoisethreshold=1.5, negativethreshold=0.0, minbeamfrac=0.3, |
128 - | growiterations=75, dogrowprune=True, minpercentchange=-1.0, fastnoise=True, savemodel='none', usemask='user', mask='', parallel=False, |
129 - | history_filter=lambda index, arg, history_value: history_value ): |
178 + | return -1, "cyclefactor must be a number" |
179 + | |
180 + | return 0, "" |
181 + | |
182 + | |
183 + | |
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", "")) |
202 + | |
203 + | |
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 ): |
211 + | |
130 212 | self._vis = vis |
131 213 | self._imagename = imagename |
132 214 | self._imsize = imsize |
133 215 | self._cell = cell |
134 216 | self._phasecenter = phasecenter |
135 217 | self._stokes = stokes |
136 218 | self._startmodel = startmodel |
137 219 | self._specmode = specmode |
138 220 | self._reffreq = reffreq |
139 221 | self._nchan = nchan |
140 222 | self._start = start |
141 223 | self._width = width |
142 224 | self._outframe = outframe |
143 - | self._interpolation = interpolation |
225 + | self._veltype = veltype |
144 226 | self._restfreq = restfreq |
227 + | self._interpolation = interpolation |
145 228 | self._perchanweightdensity = perchanweightdensity |
146 229 | self._gridder = gridder |
147 230 | self._wprojplanes = wprojplanes |
148 231 | self._mosweight = mosweight |
149 232 | self._psterm = psterm |
150 233 | self._wbawp = wbawp |
151 234 | self._conjbeams = conjbeams |
152 235 | self._usepointing = usepointing |
153 236 | self._pointingoffsetsigdev = pointingoffsetsigdev |
154 237 | self._pblimit = pblimit |
155 238 | self._deconvolver = deconvolver |
239 + | self._smallscalebias = smallscalebias |
156 240 | self._niter = niter |
157 241 | self._threshold = threshold |
158 242 | self._cycleniter = cycleniter |
243 + | self._minpsffraction = minpsffraction |
244 + | self._maxpsffraction = maxpsffraction |
159 245 | self._nsigma = nsigma |
160 246 | self._nmajor = nmajor |
161 247 | self._cyclefactor = cyclefactor |
162 248 | self._scales = scales |
163 249 | self._restoringbeam = restoringbeam |
164 250 | self._pbcor = pbcor |
165 251 | self._nterms = nterms |
166 252 | self._exe_cmds = [ ] |
253 + | self._exe_cmds_per_iter = [ ] |
167 254 | self._history_filter = history_filter |
168 255 | self._finalized = False |
169 256 | self._field = field |
170 257 | self._spw = spw |
171 258 | self._timerange = timerange |
172 259 | self._uvrange = uvrange |
173 260 | self._antenna = antenna |
174 261 | self._scan = scan |
175 262 | self._observation = observation |
176 263 | self._intent = intent |
177 264 | self._datacolumn = datacolumn |
178 265 | self._weighting = weighting |
179 266 | self._robust = robust |
180 267 | self._npixels = npixels |
181 268 | self._gain = gain |
269 + | self._pbmask = pbmask |
182 270 | self._sidelobethreshold = sidelobethreshold |
183 271 | self._noisethreshold = noisethreshold |
184 272 | self._lownoisethreshold = lownoisethreshold |
185 273 | self._negativethreshold = negativethreshold |
274 + | self._smoothfactor = smoothfactor |
186 275 | self._minbeamfrac = minbeamfrac |
276 + | self._cutthreshold = cutthreshold |
187 277 | self._growiterations = growiterations |
188 278 | self._dogrowprune = dogrowprune |
189 279 | self._minpercentchange = minpercentchange |
280 + | self._verbose = verbose |
190 281 | self._fastnoise = fastnoise |
191 282 | self._savemodel = savemodel |
192 283 | self._parallel = parallel |
193 284 | self._usemask = usemask |
194 285 | self._mask = mask |
195 - | |
286 + | self.global_imdict = ImagingDict() |
287 + | self.current_imdict = ImagingDict() |
196 288 | self._major_done = 0 |
197 - | self._convergence_result = (None,None,None,{ 'chan': None, 'major': None }) |
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 }) |
198 293 | # ^^^^ ^^^^ ^^^^ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^----->>> convergence info |
199 - | # | | | |
200 - | # | | +---------------->>> major cycles done for current run |
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) |
201 297 | # | +------------------>>> tclean stopcode |
202 - | # +----------------------->>> error message |
203 - | |
204 - | def __filter_convergence( raw ): |
205 - | ### |
206 - | ### this function filters out the pieces of the `raw` tclean 'summaryminor' |
207 - | ### return dictionary that we care about |
208 - | ### |
209 - | ### the first index in the `raw` dictionary is the channel axis |
210 - | ### each channel may have a number of polarity dictionaries |
211 - | ### |
212 - | keep_keys = [ 'modelFlux', 'iterDone', 'peakRes', 'stopCode', 'cycleThresh' ] |
213 - | ret = {} |
214 - | for channel_k,channel_v in raw[0].items( ): # 0: main field in multifield imaging TODO worry about other fields |
215 - | ret[channel_k] = {} |
216 - | for stokes_k,stokes_v in channel_v.items( ): |
217 - | ret[channel_k][stokes_k] = {} |
218 - | for summary_k in keep_keys: |
219 - | ret[channel_k][stokes_k][summary_k] = copy.deepcopy(stokes_v[summary_k]) |
220 - | return ret |
298 + | # +----------------------->>> tclean stopdescription |
299 + | |
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() |
304 + | |
221 305 | |
222 306 | def __add_per_major_items( self, tclean_ret, major_ret, chan_ret ): |
223 307 | '''Add meta-data about the whole major cycle, including 'cyclethreshold' |
224 308 | ''' |
309 + | |
225 310 | if 'cyclethreshold' in tclean_ret: |
226 - | return dict( major=dict( cyclethreshold=[tclean_ret['cyclethreshold']] if major_ret is None else (major_ret['cyclethreshold'] + [tclean_ret['cyclethreshold']]) ), |
311 + | |
312 + | rdict = dict( major=dict( cyclethreshold=[tclean_ret['cyclethreshold']] if major_ret is None else (major_ret['cyclethreshold'] + [tclean_ret['cyclethreshold']]) ), |
227 313 | chan=chan_ret ) |
228 314 | else: |
229 - | return dict( major=dict( cyclethreshold=major_ret['cyclethreshold'].append(tclean_ret['cyclethreshold']) ), |
315 + | rdict = dict( major=dict( cyclethreshold=major_ret['cyclethreshold'].append(tclean_ret['cyclethreshold']) ), |
230 316 | chan=chan_ret ) |
231 317 | |
232 - | |
233 - | def __update_convergence( cumm_sm, new_sm ): |
234 - | """Accumulates the per-channel/stokes subimage 'summaryminor' records from new_sm to cumm_sm. |
235 - | param cumm_sm: cummulative summary minor records : { chan: { stoke: { key: [values] } } } |
236 - | param new_sm: new summary minor records : { chan: { stoke: { key: [values] } } } |
318 + | return rdict |
319 + | |
320 + | |
321 + | def _calc_deconv_controls(self, imdict, niterleft=0, threshold=0, cycleniter=-1): |
322 + | """ |
323 + | Calculate cycleniter and cyclethreshold for deconvolution. |
324 + | """ |
325 + | |
326 + | use_cycleniter = niterleft #niter - imdict.returndict['iterdone'] |
327 + | |
328 + | if cycleniter > -1 : # User is forcing this number |
329 + | use_cycleniter = min(cycleniter, use_cycleniter) |
330 + | |
331 + | psffrac = imdict.returndict['maxpsfsidelobe'] * self._cyclefactor |
332 + | psffrac = max(psffrac, self._minpsffraction) |
333 + | psffrac = min(psffrac, self._maxpsffraction) |
334 + | |
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) |
339 + | |
340 + | return int(use_cycleniter), cyclethreshold |
237 341 | |
238 - | For most "keys", the resultant "values" will be a list, one value per minor cycle. |
342 + | |
343 + | def __update_convergence(self): |
344 + | """ |
345 + | Accumulates the per-channel/stokes summaryminor keys across all major cycle calls so far. |
239 346 | |
240 347 | The "iterDone" key will be replaced with "iterations", and for the "iterations" key, |
241 348 | the value in the returned cummulative record will be a rolling sum of iterations done |
242 349 | for tclean calls so far, one value per minor cycle. |
243 350 | For example, if there have been two clean calls, and in the first call channel 0 had |
244 351 | [1] iteration in 1 minor cycle, and for the second call channel 0 had [6, 10, 9, 1] |
245 352 | iterations in 4 minor cycles), then the resultant "iterations" key for channel 0 would be: |
246 353 | [1, 7, 17, 26, 27] |
247 354 | """ |
248 355 | |
249 - | ### substitute 'iterations' for 'iterDone' |
250 - | replace_iters_key = lambda x: ('iterations' if x[0] == 'iterDone' else x[0], x[1]) |
251 - | new_sm = { |
252 - | chan_k: { |
253 - | stokes_k: { |
254 - | k: v for k,v in map( replace_iters_key, stokes_v.items( ) ) |
255 - | } for stokes_k,stokes_v in chan_v.items() |
256 - | } for chan_k,chan_v in new_sm.items() |
257 - | } |
356 + | keys = ['modelFlux', 'iterDone', 'peakRes', 'stopCode', 'cycleThresh'] |
357 + | |
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) |
258 373 | |
259 - | if cumm_sm is None: |
260 - | return new_sm |
374 + | return outrec |
375 + | |
376 + | |
377 + | def _check_initial_mask(self): |
378 + | """ |
379 + | Check if a mask from a previous run exists on disk or not. |
380 + | """ |
381 + | |
382 + | if self._usemask == 'user' and self._mask == '': |
383 + | maskname = self._imagename + '.mask' |
384 + | |
385 + | if os.path.exists(maskname): |
386 + | self._initial_mask_exists = True |
387 + | else: |
388 + | self._initial_mask_exists = False |
389 + | |
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 + | """ |
395 + | |
396 + | from casatools import image |
397 + | ia = image() |
398 + | |
399 + | if self._usemask == 'user' and self._mask == '': |
400 + | maskname = self._imagename + '.mask' |
401 + | |
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() |
407 + | |
408 + | def _update_peakres(self): |
409 + | if self._deconvolver == 'mtmfs': |
410 + | residname = self._imagename + '.residual.tt0' |
411 + | else: |
412 + | residname = self._imagename + '.residual' |
413 + | |
414 + | maskname = self._imagename + '.mask' |
415 + | if not os.path.exists(maskname): |
416 + | maskname = '' |
417 + | |
418 + | peakres = imstat(imagename=residname, mask=maskname)['max'] |
419 + | if len(maskname) > 0: |
420 + | masksum = imstat(imagename=maskname)['sum'] |
421 + | else: |
422 + | masksum = [] |
423 + | |
424 + | if len(peakres) > 0: |
425 + | peakres = peakres[0] |
261 426 | else: |
262 - | def accumulate_tclean_records( cumm_subsm_rec, new_subsm_rec ): |
263 - | """ |
264 - | param cumm_subsm_rec: cummulative subimage 'summaryminor' record : { "key": [minor cycle values,...] } |
265 - | param new_subsm_rec: new subimage 'summaryminor' record : { "key": [minor cycle values,...] } |
266 - | """ |
267 - | curr_iters_sum = max(cumm_subsm_rec['iterations']) if 'iterations' in cumm_subsm_rec else 0 |
268 - | if 'iterations' in new_subsm_rec: |
269 - | iterations_tuple = reduce( lambda acc, v: (acc[0]+v, acc[1] + [acc[0]+v+curr_iters_sum]), new_subsm_rec['iterations'], (0,[]) ) |
270 - | new_subsm_rec['iterations'] = iterations_tuple[1] # just want the sum of iterations list |
271 - | return { key: cumm_subsm_rec[key] + new_subsm_rec[key] for key in new_subsm_rec.keys( ) } |
272 - | return { channel_k: { |
273 - | stokes_k: accumulate_tclean_records( cumm_sm[channel_k][stokes_k], stokes_v ) |
274 - | for stokes_k,stokes_v in channel_v.items( ) } for channel_k,channel_v in new_sm.items( ) |
275 - | } |
427 + | peakres = None |
428 + | |
429 + | if len(masksum) > 0: |
430 + | masksum = masksum[0] |
431 + | else: |
432 + | masksum = None |
433 + | |
434 + | return peakres, masksum |
276 435 | |
277 436 | def __next__( self ): |
278 437 | """ Runs tclean and returns the (stopcode, convergence result) when executed with the python builtin next() function. |
279 438 | |
280 439 | The returned convergence result is a nested dictionary: |
281 440 | { |
282 441 | channel id: { |
283 442 | stokes id: { |
284 443 | summary key: [values, one per minor cycle] |
285 444 | }, |
286 445 | }, |
287 446 | } |
288 447 | |
289 448 | See also: gclean.__update_convergence(...) |
290 449 | """ |
291 - | if self._finalized: |
292 - | self._convergence_result = ( f'iteration terminated', |
293 - | self._convergence_result[1], |
294 - | self._major_done, |
295 - | self._convergence_result[3] ) |
296 - | raise StopIteration |
297 - | # ensure that at least the initial tclean run |
450 + | |
298 451 | # vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv------------>>> is done to produce an initial dirty image |
299 452 | if self._niter < 1 and self._convergence_result[2] is not None: |
300 453 | self._convergence_result = ( f'nothing to run, niter == {self._niter}', |
301 454 | self._convergence_result[1], |
302 455 | self._major_done, |
303 - | self._convergence_result[3] ) |
456 + | self._nmajor, |
457 + | self._niter, |
458 + | self._convergence_result[5] ) |
304 459 | return self._convergence_result |
305 460 | else: |
306 - | ### |
307 461 | ### CALL SEQUENCE: |
308 462 | ### tclean(niter=0),deconvolve(niter=0),tclean(niter=100),deconvolve(niter=0),tclean(niter=100),tclean(niter=0,restoration=True) |
309 - | ### |
310 - | if self._convergence_result[1] is None: |
311 - | # initial call to tclean(...) creates the initial dirty image with niter=0 |
312 - | tclean_ret = self._tclean( vis=self._vis, mask=self._mask, imagename=self._imagename, imsize=self._imsize, cell=self._cell, |
313 - | phasecenter=self._phasecenter, stokes=self._stokes, startmodel=self._startmodel, specmode=self._specmode, |
314 - | reffreq=self._reffreq, gridder=self._gridder, wprojplanes=self._wprojplanes, mosweight=self._mosweight, |
315 - | psterm=self._psterm, wbawp=self._wbawp, conjbeams=self._conjbeams, usepointing=self._usepointing, |
316 - | interpolation=self._interpolation, restfreq=self._restfreq, perchanweightdensity=self._perchanweightdensity, |
317 - | nchan=self._nchan, start=self._start, width=self._width, outframe=self._outframe, |
318 - | pointingoffsetsigdev=self._pointingoffsetsigdev, pblimit=self._pblimit, deconvolver=self._deconvolver, |
319 - | cyclefactor=self._cyclefactor, scales=self._scales, restoringbeam=self._restoringbeam, pbcor=self._pbcor, |
320 - | nterms=self._nterms, field=self._field, spw=self._spw, timerange=self._timerange, uvrange=self._uvrange, |
321 - | antenna=self._antenna, scan=self._scan, observation=self._observation, intent=self._intent, |
322 - | datacolumn=self._datacolumn, weighting=self._weighting, robust=self._robust, npixels=self._npixels, |
323 - | interactive=False, niter=0, gain=0.000001, calcres=True, restoration=False, parallel=self._parallel, fullsummary=True ) |
324 - | self._deconvolve( imagename=self._imagename, niter=0, usemask=self._usemask, restoration=False, deconvolver=self._deconvolver ) |
325 - | self._major_done = 0 |
326 - | else: |
327 - | tclean_ret = self._tclean( vis=self._vis, imagename=self._imagename, imsize=self._imsize, cell=self._cell, phasecenter=self._phasecenter, |
328 - | stokes=self._stokes, specmode=self._specmode, reffreq=self._reffreq, |
463 + | self._exe_cmds_per_iter.append(0) |
464 + | try: |
465 + | if self._convergence_result[1] is None: |
466 + | |
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) |
481 + | |
482 + | # Check if a mask from a previous run exists on disk |
483 + | self._check_initial_mask() |
484 + | |
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) |
493 + | |
494 + | # If no mask from a previous run exists, over-write the ones with zeros for the default mask |
495 + | self._fix_initial_mask() |
496 + | |
497 + | self.current_imdict.returndict = self.current_imdict.merge(tclean_ret, deconv_ret) |
498 + | self.global_imdict.returndict = self.current_imdict.returndict |
499 + | |
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) |
502 + | |
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. |
509 + | |
510 + | # Mask can be updated here... |
511 + | # Check for mask update - peakres + masksum |
512 + | _peakres, _masksum = self._update_peakres() |
513 + | |
514 + | self.hasit, self.stopdescription = self.global_imdict.has_converged(self._niter, self._threshold, self._nmajor, masksum=_masksum, peakres=_peakres) |
515 + | |
516 + | #self.global_imdict.returndict['stopcode'] = self.hasit |
517 + | #self.global_imdict.returndict['stopDescription'] = self.stopdescription |
518 + | |
519 + | #self.current_imdict.returndict['stopcode'] = self.hasit |
520 + | #self.current_imdict.returndict['stopDescription'] = self.stopdescription |
521 + | |
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) |
525 + | |
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) |
530 + | |
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, |
329 534 | gridder=self._gridder, wprojplanes=self._wprojplanes, mosweight=self._mosweight, psterm=self._psterm, |
330 535 | wbawp=self._wbawp, conjbeams=self._conjbeams, usepointing=self._usepointing, interpolation=self._interpolation, |
331 - | restfreq=self._restfreq, perchanweightdensity=self._perchanweightdensity, nchan=self._nchan, start=self._start, |
332 - | width=self._width, outframe=self._outframe, pointingoffsetsigdev=self._pointingoffsetsigdev, pblimit=self._pblimit, |
333 - | deconvolver=self._deconvolver, cyclefactor=self._cyclefactor, scales=self._scales, restoringbeam=self._restoringbeam, |
334 - | pbcor=self._pbcor, nterms=self._nterms, field=self._field, spw=self._spw, timerange=self._timerange, |
335 - | uvrange=self._uvrange, antenna=self._antenna, scan=self._scan, observation=self._observation, intent=self._intent, |
336 - | datacolumn=self._datacolumn, weighting=self._weighting, robust=self._robust, npixels=self._npixels, interactive=False, |
337 - | niter=self._niter, restart=True, calcpsf=False, calcres=False, restoration=False, threshold=self._threshold, |
338 - | nsigma=self._nsigma, cycleniter=self._cycleniter, nmajor=self._nmajor, gain=self._gain, |
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, |
339 546 | sidelobethreshold=self._sidelobethreshold, noisethreshold=self._noisethreshold, |
340 547 | lownoisethreshold=self._lownoisethreshold, negativethreshold=self._negativethreshold, |
341 548 | minbeamfrac=self._minbeamfrac, growiterations=self._growiterations, dogrowprune=self._dogrowprune, |
342 - | minpercentchange=self._minpercentchange, fastnoise=self._fastnoise, savemodel=self._savemodel, maxpsffraction=1, |
343 - | minpsffraction=0, parallel=self._parallel, fullsummary=True ) |
344 - | self._deconvolve( imagename=self._imagename, niter=0, usemask=self._usemask, restoration=False, deconvolver=self._deconvolver ) |
345 - | self._major_done = tclean_ret['nmajordone'] if 'nmajordone' in tclean_ret else 0 |
346 - | |
347 - | if len(tclean_ret) > 0 and 'summaryminor' in tclean_ret and sum(map(len,tclean_ret['summaryminor'].values())) > 0: |
348 - | new_summaryminor_rec = gclean.__filter_convergence(tclean_ret['summaryminor']) |
349 - | self._convergence_result = ( None, |
350 - | tclean_ret['stopcode'] if 'stopcode' in tclean_ret else 0, |
351 - | self._major_done, |
352 - | self.__add_per_major_items( tclean_ret, |
353 - | self._convergence_result[3]['major'], |
354 - | gclean.__update_convergence( self._convergence_result[3]['chan'], |
355 - | new_summaryminor_rec ) ) ) |
356 - | else: |
357 - | self._convergence_result = ( f'tclean returned an empty result', |
358 - | self._convergence_result[1], |
549 + | minpercentchange=self._minpercentchange, fastnoise=self._fastnoise, savemodel=self._savemodel, |
550 + | maxpsffraction=self._maxpsffraction, |
551 + | minpsffraction=self._minpsffraction, parallel=self._parallel, fullsummary=True ) |
552 + | |
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) |
556 + | |
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'] |
560 + | |
561 + | ## Decrement count for the major cycle just done... |
562 + | self.__decrement_counts() |
563 + | |
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) |
569 + | |
570 + | |
571 + | self.global_imdict.returndict['stopcode'] = self.hasit |
572 + | self.global_imdict.returndict['stopDescription'] = self.stopdescription |
573 + | |
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) |
581 + | |
582 + | |
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, |
359 602 | self._major_done, |
360 - | self._convergence_result[3] ) |
603 + | self._nmajor, |
604 + | self._niter, |
605 + | self._convergence_result[5] ) |
606 + | return self._convergence_result |
607 + | |
361 608 | return self._convergence_result |
362 609 | |
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. |
622 + | |
623 + | |
363 624 | def __reflect_stop( self ): |
364 625 | ## if python wasn't hacky, you would be able to try/except/raise in lambda |
626 + | time.sleep(1) |
365 627 | try: |
366 628 | return self.__next__( ) |
367 629 | except StopIteration: |
368 630 | raise StopAsyncIteration |
369 631 | |
370 632 | async def __anext__( self ): |
371 633 | ### asyncio.run cannot be used here because this is called |
372 634 | ### within an asyncio loop... |
373 635 | loop = asyncio.get_event_loop( ) |
374 636 | result = await loop.run_in_executor( None, self.__reflect_stop ) |
375 637 | return result |
376 638 | |
377 639 | def __iter__( self ): |
378 640 | return self |
379 641 | |
380 642 | def __aiter__( self ): |
381 643 | return self |
382 644 | |
645 + | def __split_filename( self, path ): |
646 + | return os.path.splitext(os.path.basename(path)) |
647 + | |
648 + | def __default_mask_name( self ): |
649 + | imgparts = self.__split_filename( self._imagename ) |
650 + | return f'{imgparts[0]}.mask' |
651 + | |
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 |
655 + | |
383 656 | def reset(self): |
384 657 | #if not self._finalized: |
385 658 | # raise RuntimeError('attempt to reset a gclean run that has not been finalized') |
386 659 | self._finalized = False |
387 660 | self._convergence_result = ( None, |
388 661 | self._convergence_result[1], |
389 662 | self._major_done, |
390 - | self._convergence_result[3] ) |
663 + | self._nmajor, |
664 + | self._niter, |
665 + | self._convergence_result[5] ) |
391 666 | |
392 667 | def restore(self): |
393 668 | """ Restores the final image, and returns a path to the restored image. """ |
394 - | tclean_ret = self._tclean( vis=self._vis, imagename=self._imagename, imsize=self._imsize, cell=self._cell, |
395 - | phasecenter=self._phasecenter, stokes=self._stokes, specmode=self._specmode, |
396 - | reffreq=self._reffreq, gridder=self._gridder, wprojplanes=self._wprojplanes, mosweight=self._mosweight, |
397 - | psterm=self._psterm, wbawp=self._wbawp, conjbeams=self._conjbeams, usepointing=self._usepointing, |
398 - | interpolation=self._interpolation, restfreq=self._restfreq, perchanweightdensity=self._perchanweightdensity, nchan=self._nchan, |
399 - | start=self._start, width=self._width, outframe=self._outframe, pointingoffsetsigdev=self._pointingoffsetsigdev, |
400 - | pblimit=self._pblimit, deconvolver=self._deconvolver, cyclefactor=self._cyclefactor, scales=self._scales, |
401 - | restoringbeam=self._restoringbeam, pbcor=self._pbcor, nterms=self._nterms, field=self._field, spw=self._spw, |
402 - | timerange=self._timerange, uvrange=self._uvrange, antenna=self._antenna, scan=self._scan, |
403 - | observation=self._observation, intent=self._intent, datacolumn=self._datacolumn, weighting=self._weighting, |
404 - | robust=self._robust, npixels=self._npixels, gain=self._gain, sidelobethreshold=self._sidelobethreshold, |
405 - | noisethreshold=self._noisethreshold, lownoisethreshold=self._lownoisethreshold, |
406 - | negativethreshold=self._negativethreshold, minbeamfrac=self._minbeamfrac, growiterations=self._growiterations, |
407 - | dogrowprune=self._dogrowprune, minpercentchange=self._minpercentchange, fastnoise=self._fastnoise, |
408 - | savemodel=self._savemodel, nsigma=self._nsigma, interactive=False, |
409 - | niter=0, restart=True, calcpsf=False, calcres=False, restoration=True, |
410 - | parallel=self._parallel ) |
669 + | deconv_ret = self._deconvolve(imagename=self._imagename, |
670 + | deconvolver=self._deconvolver, |
671 + | restoration=True, niter=0, |
672 + | fullsummary=True) |
673 + | |
411 674 | return { "image": f"{self._imagename}.image" } |
412 675 | |
413 676 | def has_next(self): |
414 677 | return not self._finalized |