Coverage for /home/casatest/venv/lib/python3.12/site-packages/casatasks/private/imagerhelpers/imager_base.py: 43%

1import os 

2import math 

3import shutil 

4import string 

5import time 

6import re 

7import copy 

8from typing import TYPE_CHECKING 

9 

10 

11from casatools import ( 

12 synthesisimager, 

13 synthesisdeconvolver, 

14 synthesisnormalizer, 

15 iterbotsink, 

16 ctsys, 

17 table, 

18 image, 

19) 

20from casatasks import casalog 

21from casatasks.private.imagerhelpers.summary_minor import SummaryMinor 

22if TYPE_CHECKING: 

23 from casatasks.private.imagerhelpers.input_parameters import ImagerParameters 

24 

25ctsys_hostinfo = ctsys.hostinfo 

26_tb = table() 

27_ia = image() 

28 

29""" 

30A set of helper functions for tclean. 

31 

32Summary... 

33  

34""" 

35 

36############################################# 

37class PySynthesisImager: 

38 

39 def __init__(self,params: 'ImagerParameters'): 

40 ################ Tools 

41 self.initDefaults() 

42 

43 # Check all input parameters, after partitioning setup. 

44 

45 # Selection Parameters. Dictionary of dictionaries, indexed by 'ms0','ms1',... 

46 self.allselpars = params.getSelPars() 

47 # Imaging/Deconvolution parameters. Same for serial and parallel runs 

48 self.alldecpars = params.getDecPars() 

49 self.allimpars = params.getImagePars() 

50 self.allgridpars = params.getGridPars() 

51 self.allnormpars = params.getNormPars() 

52 self.weightpars = params.getWeightPars() 

53 # Iteration parameters 

54 self.iterpars = params.getIterPars() ## Or just params.iterpars 

55 

56 # CFCache params 

57 self.cfcachepars = params.getCFCachePars() 

58 ## Number of fields ( main + outliers ) 

59 self.NF = len(self.allimpars.keys()) 

60 self.stopMinor = {} ##[0]*self.NF 

61 for immod in range(0, self.NF): 

62 self.stopMinor[str(immod)] = 1.0 

63 ## Number of nodes. This gets set for parallel runs 

64 ## It can also be used serially to process the major cycle in pieces. 

65 self.NN = 1 

66 ## for debug mode automask incrementation only 

67 self.ncycle = 0 

68 # For the nmajor parameter 

69 self.majorCnt = 0 

70 

71 # isvalid = self.checkParameters() 

72 # if isvalid==False: 

73 # casalog.post('Invalid parameters') 

74 

75 ############################################# 

76 # def checkParameters(self): 

77 # # Copy the imagename from impars to decpars, for each field. 

78 # for immod in range(0,self.NF): 

79 # self.alldecpars[str(immod)]['imagename'] = self.allimpars[str(immod)]['imagename'] 

80 # return True 

81 

82 ############################################# 

83 def makeCFCache(self, exists): 

84 # Make the CFCache and re-load it. The following calls become 

85 # NoOps (in SynthesisImager.cc) if the gridder is not one 

86 # which uses CFCache. 

87 if exists: 

88 casalog.post("CFCache already exists") 

89 else: 

90 self.dryGridding(); 

91 self.fillCFCache(); 

92 self.reloadCFCache(); 

93 

94 

95 def initializeImagers(self): 

96 

97 ## Initialize the tool for the current node 

98 self.SItool = synthesisimager() 

99 

100 ## casalog.post('impars ', self.allimpars['0']['specmode'], 'frame', self.allimpars['0']['outframe']) 

101 ## Send in selection parameters for all MSs in the list. 

102 for mss in sorted((self.allselpars).keys()): 

103 # if(self.allimpars['0']['specmode']=='cubedata'): 

104 # self.allselpars[mss]['outframe']='Undefined' 

105 self.SItool.selectdata(self.allselpars[mss]) 

106 # self.SItool.selectdata( **(self.allselpars[mss]) ) 

107 

108 ## For each image-field, define imaging parameters 

109 # nimpars = copy.deepcopy(self.allimpars) 

110 # for fld in range(0,self.NF): 

111 # self.SItool.defineimage( **( nimpars[str(fld)] ) ) 

112 

113 # If cfcache directory already exists, assume that it is 

114 # usable and is correct. makeCFCache call then becomes a 

115 # NoOp. 

116 

117 cfCacheName='' 

118 exists=False 

119 if(self.allgridpars['0']['gridder'].startswith('awpr') ): 

120 cfCacheName=self.allgridpars['0']['cfcache']; 

121 if (cfCacheName == ''): 

122 cfCacheName = self.allimpars['0']['imagename'] + '.cf' 

123 self.allgridpars['0']['cfcache']= cfCacheName 

124 exists = (os.path.exists(cfCacheName) and os.path.isdir(cfCacheName)); 

125 else: 

126 cfCacheName='' 

127 exists=True 

128 

129 for fld in range(0,self.NF): 

130 # casalog.post("self.allimpars=",self.allimpars,"\n") 

131 

132 # print(f'####allimpars={self.allimpars[str(fld)]} \n allgridpars={self.allgridpars[str(fld)]}') 

133 self.SItool.defineimage( 

134 self.allimpars[str(fld)], self.allgridpars[str(fld)] 

135 ) 

136 

137 ###for cases when synthesisnormalizer is setup in c++ send the normalizer info 

138 ###all images have the same normtype etc..so first one is good enough 

139 self.SItool.normalizerinfo(self.allnormpars["0"]) 

140 ###commenting this out so that tuneSelect is done after weighting 

141 ###CAS-11687 

142 # For cube imaging: align the data selections and image setup 

143 

144 #if self.allimpars['0']['specmode'] != 'mfs' and self.allimpars['0']['specmode'] != 'cubedata': 

145 # self.SItool.tuneselectdata() 

146 ###For cubes create cfcache ahead of each partition trying 

147 ### to create it as it is not multiprocess safe 

148 if(("cube" in self.allimpars['0']['specmode']) or ("awphpg" in self.allgridpars['0']['gridder'])): 

149 self.makeCFCache(exists); 

150 ### Warning about awp2/mosaic not having conjbeam thus will not be correct on first major cycles 

151 ## CAS-14146 : Krishna : Moved this warning to task_tclean.py along with the other warnings there. 

152 #if( ("mfs" in self.allimpars['0']['specmode']) and ("mtmfs" in self.allimpars['0']['deconvolver']) and (self.allgridpars['0']['gridder'] in ['awp2', 'mosaic']) ): 

153 # casalog.post( 

154 # "You may consider using specmode=mvc with "+self.allgridpars['0']['gridder'] 

155 # +" as this gridder does not use conjbeams \n thus need a couple of major cycle to converge to the correct answer", 

156 # "WARN" 

157 # ) 

158 

159 

160 ############################################# 

161 

162 def initializeDeconvolvers(self): 

163 for immod in range(0, self.NF): 

164 self.SDtools.append(synthesisdeconvolver()) 

165 self.SDtools[immod].setupdeconvolution(decpars=self.alldecpars[str(immod)]) 

166 

167 ############################################# 

168 ## Overloaded by ParallelCont 

169 def initializeNormalizers(self): 

170 for immod in range(0, self.NF): 

171 self.PStools.append(synthesisnormalizer()) 

172 normpars = self.allnormpars[str(immod)] 

173 self.PStools[immod].setupnormalizer(normpars=normpars) 

174 

175 ############################################# 

176 

177 def initializeIterationControl(self): 

178 # note that in CASA5 this is casac.synthesisiterbot 

179 self.IBtool = iterbotsink() 

180 itbot = self.IBtool.setupiteration(iterpars=self.iterpars) 

181 

182 ############################################# 

183 def estimatememory(self): 

184 # casalog.post("MEMORY usage ", self.SItool.estimatememory(), type(self.SItool.estimatememory())) 

185 # griddermem=0 

186 if self.SItool != None: 

187 griddermem = self.SItool.estimatememory() 

188 deconmem = 0 

189 for immod in range(0, self.NF): 

190 ims = self.allimpars[str(immod)]["imsize"] 

191 if type(ims) == int: 

192 ims = [ims, ims] 

193 if len(ims) == 1: 

194 ims.append(ims[0]) 

195 # casalog.post('shape', self.allimpars[str(immod)]['imsize'], len(ims) ) 

196 # casalog.post("DECON mem usage ", self.SDtools[immod].estimatememory(ims)) 

197 if len(self.SDtools) > immod: 

198 if self.SDtools != None: 

199 deconmem += self.SDtools[immod].estimatememory(ims) 

200 availmem = ctsys_hostinfo()["memory"]["available"] 

201 if (deconmem + griddermem) > 0.8 * availmem: 

202 casalog.post( 

203 "Memory available " 

204 + str(availmem) 

205 + " kB is very close to amount of required memory " 

206 + str(deconmem + griddermem) 

207 + " kB", 

208 "WARN", 

209 ) 

210 else: 

211 casalog.post( 

212 "Memory available " 

213 + str(availmem) 

214 + " kB and required memory " 

215 + str(deconmem + griddermem) 

216 + " kB", 

217 "INFO2", 

218 ) 

219 

220 ############################################ 

221 def restoreImages(self): 

222 for immod in range(0, self.NF): 

223 self.SDtools[immod].restore() 

224 

225 ############################################# 

226 def pbcorImages(self): 

227 for immod in range(0, self.NF): 

228 self.SDtools[immod].pbcor() 

229 

230 def getSummary(self,fullsummary,fignum=1): 

231 summ = self.IBtool.getiterationsummary() 

232 casalog.post('getSummary call: fullsummary='+str(fullsummary)) 

233 if ('stopcode' in summ): 

234 summ['stopDescription'] = self.getStopDescription(summ['stopcode']) 

235 if ('summaryminor' in summ): 

236 summ['summaryminor'] = SummaryMinor.convertMatrix(summ['summaryminor'],fullsummary) 

237 #self.plotReport( summ, fignum ) 

238 

239 return summ 

240 

241 ############################################# 

242 def deleteImagers(self): 

243 if self.SItool != None: 

244 self.SItool.done() 

245 self.SItool = None 

246 

247 def deleteDeconvolvers(self): 

248 for immod in range(0, len(self.SDtools)): 

249 self.SDtools[immod].done() 

250 

251 def deleteNormalizers(self): 

252 for immod in range(0, len(self.PStools)): 

253 self.PStools[immod].done() 

254 

255 def deleteIterBot(self): 

256 if self.IBtool != None: 

257 self.IBtool.done() 

258 

259 def deleteCluster(self): 

260 # casalog.post('no cluster to delete') 

261 return 

262 

263 def deleteWorkDir(self): 

264 # No .workdirectory to delete 

265 return 

266 

267 def initDefaults(self): 

268 # Reset globals/members 

269 self.NF = 1 

270 self.stopMinor = {"0": 1.0} # Flag to call minor cycle for this field or not. 

271 self.NN = 1 

272 self.SItool = None 

273 self.SDtools = [] 

274 self.PStools = [] 

275 self.IBtool = None 

276 

277 ############################################# 

278 

279 def deleteTools(self): 

280 self.deleteImagers() 

281 self.deleteDeconvolvers() 

282 self.deleteNormalizers() 

283 self.deleteIterBot() 

284 self.deleteWorkDir() 

285 self.initDefaults() 

286 self.deleteCluster() 

287 

288 ############################################# 

289 

290 def getStopDescription(self, stopflag): 

291 stopreasons = [ 

292 "iteration limit", # 1 

293 "threshold", # 2 

294 "force stop", # 3 

295 "no change in peak residual across two major cycles", # 4 

296 "peak residual increased by more than 3 times from the previous major cycle", # 5 

297 "peak residual increased by more than 3 times from the minimum reached", # 6 

298 "zero mask", # 7 

299 "any combination of n-sigma and other valid exit criterion", # 8 

300 "reached nmajor", # 9 

301 ] 

302 if stopflag > 0: 

303 return stopreasons[stopflag - 1] 

304 return None 

305 

306 def hasConverged(self): 

307 # Merge peak-res info from all fields to decide iteration parameters 

308 time0 = time.time() 

309 self.IBtool.resetminorcycleinfo() 

310 for immod in range(0, self.NF): 

311 initrec = self.SDtools[immod].initminorcycle() 

312 # print('INIT Minor cycle dict {}'.format(initrec)) 

313 self.IBtool.mergeinitrecord(initrec, immod) 

314 

315 # # Run interactive masking (and threshold/niter editors) 

316 # self.runInteractiveGUI2() 

317 

318 # Check with the iteration controller about convergence. 

319 reachedNmajor = ( 

320 self.iterpars["nmajor"] >= 0 and self.majorCnt >= self.iterpars["nmajor"] 

321 ) 

322 stopflag = self.IBtool.cleanComplete(reachedMajorLimit=reachedNmajor) 

323 if stopflag > 0: 

324 casalog.post( 

325 "Reached global stopping criterion : " 

326 + self.getStopDescription(stopflag), 

327 "INFO", 

328 ) 

329 

330 # revert the current automask to the previous one 

331 # if self.iterpars['interactive']: 

332 for immod in range(0, self.NF): 

333 if self.alldecpars[str(immod)]["usemask"].count("auto") > 0: 

334 prevmask = self.allimpars[str(immod)]["imagename"] + ".prev.mask" 

335 if os.path.isdir(prevmask): 

336 # Try to force rmtree even with an error as an nfs mounted disk gives an error 

337 # shutil.rmtree(self.allimpars[str(immod)]['imagename']+'.mask') 

338 shutil.rmtree( 

339 self.allimpars[str(immod)]["imagename"] + ".mask", 

340 ignore_errors=True, 

341 ) 

342 # For NFS mounted disk it still leave .nfs* file(s) 

343 if os.path.isdir( 

344 self.allimpars[str(immod)]["imagename"] + ".mask" 

345 ): 

346 import glob 

347 

348 if glob.glob( 

349 self.allimpars[str(immod)]["imagename"] + ".mask/.nfs*" 

350 ): 

351 for item in os.listdir(prevmask): 

352 src = os.path.join(prevmask, item) 

353 dst = os.path.join( 

354 self.allimpars[str(immod)]["imagename"] 

355 + ".mask", 

356 item, 

357 ) 

358 if os.path.isdir(src): 

359 shutil.move(src, dst) 

360 else: 

361 shutil.copy2(src, dst) 

362 shutil.rmtree(prevmask) 

363 else: 

364 shutil.move( 

365 prevmask, 

366 self.allimpars[str(immod)]["imagename"] + ".mask", 

367 ) 

368 casalog.post( 

369 "[" 

370 + str(self.allimpars[str(immod)]["imagename"]) 

371 + "] : Reverting output mask to one that was last used ", 

372 "INFO", 

373 ) 

374 casalog.post( 

375 "***Time taken in checking hasConverged " + str(time.time() - time0), 

376 "INFO3", 

377 ) 

378 return stopflag > 0 

379 

380 ############################################# 

381 def updateMask(self): 

382 # Setup mask for each field ( input mask, and automask ) 

383 maskchanged = False 

384 time0 = time.time() 

385 for immod in range(0, self.NF): 

386 maskchanged = maskchanged | self.SDtools[immod].setupmask() 

387 

388 # Run interactive masking (and threshold/niter editors), if interactive=True 

389 maskchanged = maskchanged | self.runInteractiveGUI2() 

390 

391 time1 = time.time() 

392 casalog.post("Time to update mask " + str(time1 - time0) + "s", "INFO3") 

393 ## Return a flag to say that the mask has changed or not. 

394 return maskchanged 

395 

396 ############################################# 

397 def runInteractiveGUI2(self): 

398 maskchanged = False 

399 forcestop = True 

400 if self.iterpars["interactive"] == True: 

401 self.stopMinor = self.IBtool.pauseforinteraction() 

402 # casalog.post("Actioncodes in python : " , self.stopMinor) 

403 

404 for akey in self.stopMinor: 

405 if self.stopMinor[akey] < 0: 

406 maskchanged = True 

407 self.stopMinor[akey] = abs(self.stopMinor[akey]) 

408 

409 # Check if force-stop has happened while savemodel != "none". 

410 # If so, warn the user that unless the Last major cycle has happened, 

411 # the model won't have been written into the MS, and to do a 'predict' run. 

412 forcestop = True 

413 for akey in self.stopMinor: 

414 forcestop = forcestop and self.stopMinor[akey] == 3 

415 

416 if self.iterpars["savemodel"] != "none": 

417 if forcestop == True: 

418 self.predictModel() 

419 # if self.iterpars['savemodel'] == "modelcolumn": 

420 # wstr = "Saving model column" 

421 # else: 

422 # wstr = "Saving virtual model" 

423 # casalog.post("Model visibilities may not have been saved in the MS even though you have asked for it. Please check the logger for the phrases 'Run (Last) Major Cycle' and '" + wstr +"'. If these do not appear, then please save the model via a separate tclean run with niter=0,calcres=F,calcpsf=F. It will pick up the existing model from disk and save/predict it. Reason for this : For performance reasons model visibilities are saved only in the last major cycle. If the X button on the interactive GUI is used to terminate a run before this automatically detected 'last' major cycle, the model isn't written. However, a subsequent tclean run as described above will predict and save the model. ","WARN") 

424 

425 # casalog.post('Mask changed during interaction : ', maskchanged) 

426 return maskchanged or forcestop 

427 

428 ############################################# 

429 def makePSF(self): 

430 self.makePSFCore() 

431 divideInPython = ( 

432 self.allimpars["0"]["specmode"] == "mfs" 

433 or self.allimpars["0"]["deconvolver"] == "mtmfs" 

434 ) 

435 

436 ### Gather PSFs (if needed) and normalize by weight 

437 for immod in range(0, self.NF): 

438 # for cube normalization is done in C++ 

439 if divideInPython: 

440 self.PStools[immod].gatherpsfweight() 

441 self.PStools[immod].dividepsfbyweight() 

442 # continuum A style gridders need their .weight divided by sumwt except for awphpg  

443 if(("awphpg" not in self.allgridpars['0']['gridder'])): 

444 self.PStools[immod].divideweightbysumwt() 

445 self.check_psf(immod) 

446 

447 def check_psf(self, immod): 

448 if self.SDtools != []: 

449 if immod <= len(self.SDtools) - 1: 

450 self.SDtools[immod].checkrestoringbeam() 

451 

452 ############################################# 

453 def calcVisAppSens(self): 

454 

455 return self.SItool.apparentsens() 

456 

457 ############################################# 

458 

459 def runMajorCycle(self, isCleanCycle=True): 

460 # @param isCleanCycle: true if this is part of the major/minor cleaning loop, false if this is being used for some other purpose (such as generating the residual image) 

461 if self.IBtool != None: 

462 lastcycle = self.IBtool.cleanComplete(lastcyclecheck=True) > 0 

463 else: 

464 lastcycle = True 

465 

466 divideInPython = ( 

467 self.allimpars["0"]["specmode"] == "mfs" 

468 or self.allimpars["0"]["deconvolver"] == "mtmfs" 

469 ) 

470 ##norm is done in C++ for cubes 

471 if not divideInPython: 

472 self.runMajorCycleCore(lastcycle) 

473 if isCleanCycle: 

474 self.majorCnt += 1 

475 if self.IBtool != None: 

476 self.IBtool.endmajorcycle() 

477 return 

478 

479 for immod in range(0, self.NF): 

480 self.PStools[immod].dividemodelbyweight() 

481 self.PStools[immod].scattermodel() 

482 self.runMajorCycleCore(lastcycle) 

483 if isCleanCycle: 

484 self.majorCnt += 1 

485 

486 if self.IBtool != None: 

487 self.IBtool.endmajorcycle() 

488 ### Gather residuals (if needed) and normalize by weight 

489 for immod in range(0, self.NF): 

490 self.PStools[immod].gatherresidual() 

491 self.PStools[immod].divideresidualbyweight() 

492 self.PStools[immod].multiplymodelbyweight() 

493 

494 ############################################# 

495 def predictModel(self): 

496 for immod in range(0, self.NF): 

497 self.PStools[immod].dividemodelbyweight() 

498 self.PStools[immod].scattermodel() 

499 

500 self.predictModelCore() 

501 ###return the model images back to whatever state they were 

502 for immod in range(0, self.NF): 

503 self.PStools[immod].multiplymodelbyweight() 

504 

505 ## self.SItool.predictmodel(); 

506 

507 ############################################# 

508 def dryGridding(self): 

509 self.SItool.drygridding(**(self.cfcachepars)) 

510 

511 ############################################# 

512 ## Overloaded for parallel runs 

513 def fillCFCache(self): 

514 cfcName = self.allgridpars["0"]["cfcache"] 

515 cflist = [] 

516 if not (cfcName == ""): 

517 cflist = [f for f in os.listdir(cfcName) if re.match(r"CFS*", f)] 

518 # cflist = ["CFS_0_0_CF_0_0_0.im"]; 

519 self.cfcachepars["cflist"] = cflist 

520 

521 # self.SItool.fillcfcache(**(self.cfcachepars), self.allgridpars['0']['gridder'],cfcName); 

522 

523 self.SItool.fillcfcache( 

524 cflist, 

525 self.allgridpars["0"]["gridder"], 

526 cfcName, 

527 self.allgridpars["0"]["psterm"], 

528 self.allgridpars["0"]["aterm"], 

529 self.allgridpars["0"]["conjbeams"], 

530 ) 

531 

532 ############################################# 

533 def reloadCFCache(self): 

534 self.SItool.reloadcfcache() 

535 

536 ############################################# 

537 def makePB(self): 

538 ###for cube standard gridder pb is made in c++ with psf 

539 if not ( 

540 "stand" in self.allgridpars["0"]["gridder"] 

541 and "cube" in self.allimpars["0"]["specmode"] 

542 ): 

543 self.makePBCore() 

544 for immod in range(0, self.NF): 

545 self.PStools[immod].normalizeprimarybeam() 

546 

547 ############################################# 

548 def makePBCore(self): 

549 self.SItool.makepb() 

550 

551 ############################################# 

552 def checkPB(self): 

553 """Checks for common problem cases in the .pb image""" 

554 if self.SItool is None: 

555 # Seems to be None for specmode='mfs', parallel=True 

556 return 

557 

558 import numpy as np 

559 

560 facetIdx = 0 # TODO iterate over facets 

561 imagename = self.SItool.getImageName(facetIdx, "PB") 

562 _ia.open(imagename) 

563 # Case 1: non-zeroes on edge of .pb 

564 pixelVals = _ia.getregion().copy() 

565 pixelVals[1:-2][ 

566 1:-2 

567 ] = 0 # zero out everything that isn't at the edge of 'right ascension' and 'declination' indexes 

568 if pixelVals.max() > 0: 

569 idx = np.unravel_index([pixelVals.argmax()], pixelVals.shape) 

570 idx = [ 

571 x[0] for x in idx 

572 ] # (array([296]), array([147]), array([0]), array([0])) --> [296, 147, 0, 0] 

573 casalog.post( 

574 f"Warning! Non-zero values at the edge of the .pb image can cause unexpected aliasing effects! (found value {pixelVals.max()} at index {idx})", 

575 "WARN", 

576 ) 

577 # release the image 

578 _ia.close() 

579 _ia.done() 

580 

581 ############################################# 

582 def makeSdImage(self): 

583 self.makeSdImageCore() 

584 if TYPE_CHECKING: 

585 from casatools.synthesisnormalizer import synthesisnormalizer 

586 for immod in range(0, self.NF): 

587 normalizer: synthesisnormalizer = self.PStools[immod] 

588 normalizer.gatherresidual() 

589 normalizer.divideresidualbyweight(singledish=True) 

590 

591 ############################################# 

592 def makeSdPSF(self): 

593 self.makeSdPSFCore() 

594 for immod in range(0, self.NF): 

595 self.PStools[immod].gatherresidual() 

596 self.PStools[immod].dividepsfbyweight() 

597 

598 ############################################# 

599 def makeSdImageCore(self): 

600 self.SItool.makesdimage() 

601 

602 ############################################# 

603 def makeSdPSFCore(self): 

604 self.SItool.makesdpsf() 

605 

606 ############################################# 

607 def makeImage( 

608 self, imagetype="observed", image="", compleximage="", imagefieldid=0 

609 ): 

610 """ 

611 This should replace makeSDImage, makeSDPSF and makePSF 

612 etc in the long run 

613 But for now you can do the following images i.e string recognized by type 

614 "observed", "model", "corrected", "psf", "residual", "singledish-observed", 

615 "singledish", "coverage", "holography", "holography-observed" 

616 For holography the FTmachine should be SDGrid and the baselines 

617 selected should be those that are pointed up with the antenna which is rastering. 

618 """ 

619 self.SItool.makeimage(imagetype, image, compleximage, imagefieldid) 

620 

621 ############################################# 

622 

623 ## Overloaded for parallel runs 

624 def setWeighting(self): 

625 ## Set weighting parameters, and all pars common to all fields. 

626 self.SItool.setweighting(**(self.weightpars)) 

627 ##moved the tuneselect after weighting so that 

628 ##the weight densities use all the data selected CAS-11687 

629 ###For cube imaging: align the data selections and image setup 

630 ### the tuneSelect is now done in C++ CubeMajorCycleAlgorith.cc 

631 # if self.allimpars['0']['specmode'] != 'mfs' and self.allimpars['0']['specmode'] != 'cubedata': 

632 # self.SItool.tuneselectdata() 

633 

634 # casalog.post("get set density from python") 

635 # self.SItool.getweightdensity() 

636 # self.SItool.setweightdensity() 

637 

638 ############################################# 

639 ## Overloaded for parallel runs 

640 def makePSFCore(self): 

641 self.SItool.makepsf() 

642 

643 ############################################# 

644 ## Overloaded for parallel runs 

645 def runMajorCycleCore(self, lastcycle): 

646 controldict = {"lastcycle": lastcycle} 

647 if ("0" in self.alldecpars) and ("usemask" in self.alldecpars["0"]): 

648 controldict["usemask"] = self.alldecpars["0"]["usemask"] 

649 self.SItool.executemajorcycle(controls=controldict) 

650 

651 ############################################# 

652 ## Overloaded for parallel runs 

653 def predictModelCore(self): 

654 self.SItool.predictmodel() 

655 

656 ############################################# 

657 

658 def runMinorCycle(self): 

659 return self.runMinorCycleCore() 

660 

661 ############################################# 

662 

663 def runMinorCycleCore(self): 

664 

665 # Set False for release packages. 

666 # Only set this to True for testing and debugging automask in parallel mode 

667 # since in parallel mode, runtime setting of the enviroment variable 

668 # currently does not work. 

669 # False = disable always save intermediate images mode 

670 alwaysSaveIntermediateImages = False 

671 

672 # Get iteration control parameters 

673 iterbotrec = self.IBtool.getminorcyclecontrols() 

674 

675 # TT debug - comment out after debugging.... 

676 casalog.post("Minor Cycle controls : " + str(iterbotrec)) 

677 

678 self.IBtool.resetminorcycleinfo() 

679 

680 # 

681 # Run minor cycle 

682 self.ncycle += 1 

683 retval = False 

684 for immod in range(0, self.NF): 

685 if self.stopMinor[str(immod)] < 3: 

686 

687 # temporarily disable the check (=> always save the intermediate images 

688 if alwaysSaveIntermediateImages or ( 

689 "SAVE_ALL_RESIMS" in os.environ 

690 and os.environ["SAVE_ALL_RESIMS"] == "true" 

691 ): 

692 resname = self.allimpars[str(immod)]["imagename"] + ".residual" 

693 tempresname = ( 

694 self.allimpars[str(immod)]["imagename"] 

695 + ".inputres" 

696 + str(self.ncycle) 

697 ) 

698 if os.path.isdir(resname): 

699 shutil.copytree(resname, tempresname) 

700 modname = self.allimpars[str(immod)]["imagename"] + ".model" 

701 tempmodname = ( 

702 self.allimpars[str(immod)]["imagename"] 

703 + ".inputmod" 

704 + str(self.ncycle) 

705 ) 

706 if os.path.isdir(modname): 

707 shutil.copytree(modname, tempmodname) 

708 

709 exrec = self.SDtools[immod].executeminorcycle(iterbotrecord=iterbotrec) 

710 

711 # casalog.post('.... iterdone for ', immod, ' : ' , exrec['iterdone']) 

712 retval = retval or exrec["iterdone"] > 0 

713 self.IBtool.mergeexecrecord(exrec, immod) 

714 if alwaysSaveIntermediateImages or ( 

715 "SAVE_ALL_AUTOMASKS" in os.environ 

716 and os.environ["SAVE_ALL_AUTOMASKS"] == "true" 

717 ): 

718 maskname = self.allimpars[str(immod)]["imagename"] + ".mask" 

719 tempmaskname = ( 

720 self.allimpars[str(immod)]["imagename"] 

721 + ".autothresh" 

722 + str(self.ncycle) 

723 ) 

724 if os.path.isdir(maskname): 

725 shutil.copytree(maskname, tempmaskname) 

726 

727 # Some what duplicated as above but keep a copy of the previous mask 

728 # for interactive automask to revert to it if the current mask 

729 # is not used (i.e. reached deconvolution stopping condition). 

730 ## no longer needed as of CAS-9386 for cubes. 

731 # if self.iterpars['interactive'] and self.alldecpars[str(immod)]['usemask']=='auto-thresh': 

732 if self.alldecpars[str(immod)]["usemask"].count("auto") > 0: 

733 maskname = self.allimpars[str(immod)]["imagename"] + ".mask" 

734 prevmaskname = ( 

735 self.allimpars[str(immod)]["imagename"] + ".prev.mask" 

736 ) 

737 if os.path.isdir(maskname): 

738 if os.path.isdir(prevmaskname): 

739 shutil.rmtree(prevmaskname) 

740 shutil.copytree(maskname, prevmaskname) 

741 return retval 

742 

743 ############################################# 

744 def runMajorMinorLoops(self): 

745 self.runMajorCycle(isCleanCycle=False) 

746 while not self.hasConverged(): 

747 self.runMinorCycle() 

748 self.runMajorCycle() 

749 

750 ############################################# 

751 

752 def plotReport(self, summ={}, fignum=1): 

753 

754 if not ( 

755 "summaryminor" in summ 

756 and "summarymajor" in summ 

757 and "threshold" in summ 

758 and summ["summaryminor"].shape[0] == 6 

759 ): 

760 casalog.post( 

761 "Cannot make summary plot. Please check contents of the output dictionary from tclean." 

762 ) 

763 return summ 

764 

765 import matplotlib.pyplot as pl 

766 from numpy import max as amax 

767 

768 # 0 : iteration number (within deconvolver, per cycle) 

769 # 1 : peak residual 

770 # 2 : model flux 

771 # 3 : cyclethreshold 

772 # 4 : deconvolver id 

773 # 5 : subimage id (channel, stokes..) 

774 

775 pl.ioff() 

776 

777 fig, ax = pl.subplots(nrows=1, ncols=1, num=fignum) 

778 pl.clf() 

779 minarr = summ["summaryminor"] 

780 if minarr.size == 0: 

781 casalog.post("Zero iteration: no summary plot is generated.", "WARN") 

782 else: 

783 iterlist = minarr[0, :] 

784 eps = 0.0 

785 peakresstart = [] 

786 peakresend = [] 

787 modfluxstart = [] 

788 modfluxend = [] 

789 itercountstart = [] 

790 itercountend = [] 

791 peakresstart.append(minarr[1, :][0]) 

792 modfluxstart.append(minarr[2, :][0]) 

793 itercountstart.append(minarr[0, :][0] + eps) 

794 peakresend.append(minarr[1, :][0]) 

795 modfluxend.append(minarr[2, :][0]) 

796 itercountend.append(minarr[0, :][0] + eps) 

797 for ii in range(0, len(iterlist) - 1): 

798 if iterlist[ii] == iterlist[ii + 1]: 

799 peakresend.append(minarr[1, :][ii]) 

800 peakresstart.append(minarr[1, :][ii + 1]) 

801 modfluxend.append(minarr[2, :][ii]) 

802 modfluxstart.append(minarr[2, :][ii + 1]) 

803 itercountend.append(iterlist[ii] - eps) 

804 itercountstart.append(iterlist[ii] + eps) 

805 

806 peakresend.append(minarr[1, :][len(iterlist) - 1]) 

807 modfluxend.append(minarr[2, :][len(iterlist) - 1]) 

808 itercountend.append(minarr[0, :][len(iterlist) - 1] + eps) 

809 

810 # pl.plot( iterlist , minarr[1,:] , 'r.-' , label='peak residual' , linewidth=1.5, markersize=8.0) 

811 # pl.plot( iterlist , minarr[2,:] , 'b.-' , label='model flux' ) 

812 # pl.plot( iterlist , minarr[3,:] , 'g--' , label='cycle threshold' ) 

813 

814 pl.plot(itercountstart, peakresstart, "r.--", label="peak residual (start)") 

815 pl.plot( 

816 itercountend, 

817 peakresend,