Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/imagerhelpers/imager_base.py: 16%

1import os 

2import math 

3import shutil 

4import string 

5import time 

6import re 

7import copy 

8from typing import TYPE_CHECKING 

9 

10from casatools import ( 

11 synthesisimager, 

12 synthesisdeconvolver, 

13 synthesisnormalizer, 

14 iterbotsink, 

15 ctsys, 

16 table, 

17 image, 

18) 

19from casatasks import casalog 

20from casatasks.private.imagerhelpers.summary_minor import SummaryMinor 

21if TYPE_CHECKING: 

22 from casatasks.private.imagerhelpers.input_parameters import ImagerParameters 

23 

24ctsys_hostinfo = ctsys.hostinfo 

25_tb = table() 

26_ia = image() 

27 

28""" 

29A set of helper functions for tclean. 

30 

31Summary... 

32  

33""" 

34 

35############################################# 

36class PySynthesisImager: 

37 

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

39 ################ Tools 

40 self.initDefaults() 

41 

42 # Check all input parameters, after partitioning setup. 

43 

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

45 self.allselpars = params.getSelPars() 

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

47 self.alldecpars = params.getDecPars() 

48 self.allimpars = params.getImagePars() 

49 self.allgridpars = params.getGridPars() 

50 self.allnormpars = params.getNormPars() 

51 self.weightpars = params.getWeightPars() 

52 # Iteration parameters 

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

54 

55 # CFCache params 

56 self.cfcachepars = params.getCFCachePars() 

57 ## Number of fields ( main + outliers ) 

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

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

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

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

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

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

64 self.NN = 1 

65 ## for debug mode automask incrementation only 

66 self.ncycle = 0 

67 # For the nmajor parameter 

68 self.majorCnt = 0 

69 

70 # isvalid = self.checkParameters() 

71 # if isvalid==False: 

72 # casalog.post('Invalid parameters') 

73 

74 ############################################# 

75 # def checkParameters(self): 

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

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

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

79 # return True 

80 

81 ############################################# 

82 def makeCFCache(self, exists): 

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

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

85 # which uses CFCache. 

86 if exists: 

87 casalog.post("CFCache already exists") 

88 else: 

89 self.dryGridding(); 

90 self.fillCFCache(); 

91 self.reloadCFCache(); 

92 

93 

94 def initializeImagers(self): 

95 

96 ## Initialize the tool for the current node 

97 self.SItool = synthesisimager() 

98 

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

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

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

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

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

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

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

106 

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

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

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

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

111 

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

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

114 # NoOp. 

115 

116 cfCacheName='' 

117 exists=False 

118 if(self.allgridpars['0']['gridder'].startswith('awpr') or self.allgridpars['0']['gridder'].startswith('awph') ): 

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

120 if (cfCacheName == ''): 

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

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

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

124 else: 

125 cfCacheName='' 

126 exists=True 

127 

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

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

130 

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

132 self.SItool.defineimage( 

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

134 ) 

135 

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

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

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

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

140 ###CAS-11687 

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

142 

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

144 # self.SItool.tuneselectdata() 

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

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

147 if("cube" in self.allimpars['0']['specmode']): 

148 self.makeCFCache(exists); 

149 

150 ############################################# 

151 

152 def initializeDeconvolvers(self): 

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

154 self.SDtools.append(synthesisdeconvolver()) 

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

156 

157 ############################################# 

158 ## Overloaded by ParallelCont 

159 def initializeNormalizers(self): 

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

161 self.PStools.append(synthesisnormalizer()) 

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

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

164 

165 ############################################# 

166 

167 def initializeIterationControl(self): 

168 # note that in CASA5 this is casac.synthesisiterbot 

169 self.IBtool = iterbotsink() 

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

171 

172 ############################################# 

173 def estimatememory(self): 

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

175 # griddermem=0 

176 if self.SItool != None: 

177 griddermem = self.SItool.estimatememory() 

178 deconmem = 0 

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

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

181 if type(ims) == int: 

182 ims = [ims, ims] 

183 if len(ims) == 1: 

184 ims.append(ims[0]) 

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

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

187 if len(self.SDtools) > immod: 

188 if self.SDtools != None: 

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

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

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

192 casalog.post( 

193 "Memory available " 

194 + str(availmem) 

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

196 + str(deconmem + griddermem) 

197 + " kB", 

198 "WARN", 

199 ) 

200 else: 

201 casalog.post( 

202 "Memory available " 

203 + str(availmem) 

204 + " kB and required memory " 

205 + str(deconmem + griddermem) 

206 + " kB", 

207 "INFO2", 

208 ) 

209 

210 ############################################ 

211 def restoreImages(self): 

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

213 self.SDtools[immod].restore() 

214 

215 ############################################# 

216 def pbcorImages(self): 

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

218 self.SDtools[immod].pbcor() 

219 

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

221 summ = self.IBtool.getiterationsummary() 

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

223 if ('stopcode' in summ): 

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

225 if ('summaryminor' in summ): 

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

227 #self.plotReport( summ, fignum ) 

228 

229 return summ 

230 

231 ############################################# 

232 def deleteImagers(self): 

233 if self.SItool != None: 

234 self.SItool.done() 

235 self.SItool = None 

236 

237 def deleteDeconvolvers(self): 

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

239 self.SDtools[immod].done() 

240 

241 def deleteNormalizers(self): 

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

243 self.PStools[immod].done() 

244 

245 def deleteIterBot(self): 

246 if self.IBtool != None: 

247 self.IBtool.done() 

248 

249 def deleteCluster(self): 

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

251 return 

252 

253 def deleteWorkDir(self): 

254 # No .workdirectory to delete 

255 return 

256 

257 def initDefaults(self): 

258 # Reset globals/members 

259 self.NF = 1 

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

261 self.NN = 1 

262 self.SItool = None 

263 self.SDtools = [] 

264 self.PStools = [] 

265 self.IBtool = None 

266 

267 ############################################# 

268 

269 def deleteTools(self): 

270 self.deleteImagers() 

271 self.deleteDeconvolvers() 

272 self.deleteNormalizers() 

273 self.deleteIterBot() 

274 self.deleteWorkDir() 

275 self.initDefaults() 

276 self.deleteCluster() 

277 

278 ############################################# 

279 

280 def getStopDescription(self, stopflag): 

281 stopreasons = [ 

282 "iteration limit", # 1 

283 "threshold", # 2 

284 "force stop", # 3 

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

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

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

288 "zero mask", # 7 

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

290 "reached nmajor", # 9 

291 ] 

292 if stopflag > 0: 

293 return stopreasons[stopflag - 1] 

294 return None 

295 

296 def hasConverged(self): 

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

298 time0 = time.time() 

299 self.IBtool.resetminorcycleinfo() 

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

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

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

303 self.IBtool.mergeinitrecord(initrec) 

304 

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

306 # self.runInteractiveGUI2() 

307 

308 # Check with the iteration controller about convergence. 

309 reachedNmajor = ( 

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

311 ) 

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

313 if stopflag > 0: 

314 casalog.post( 

315 "Reached global stopping criterion : " 

316 + self.getStopDescription(stopflag), 

317 "INFO", 

318 ) 

319 

320 # revert the current automask to the previous one 

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

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

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

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

325 if os.path.isdir(prevmask): 

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

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

328 shutil.rmtree( 

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

330 ignore_errors=True, 

331 ) 

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

333 if os.path.isdir( 

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

335 ): 

336 import glob 

337 

338 if glob.glob( 

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

340 ): 

341 for item in os.listdir(prevmask): 

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

343 dst = os.path.join( 

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

345 + ".mask", 

346 item, 

347 ) 

348 if os.path.isdir(src): 

349 shutil.move(src, dst) 

350 else: 

351 shutil.copy2(src, dst) 

352 shutil.rmtree(prevmask) 

353 else: 

354 shutil.move( 

355 prevmask, 

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

357 ) 

358 casalog.post( 

359 "[" 

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

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

362 "INFO", 

363 ) 

364 casalog.post( 

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

366 "INFO3", 

367 ) 

368 return stopflag > 0 

369 

370 ############################################# 

371 def updateMask(self): 

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

373 maskchanged = False 

374 time0 = time.time() 

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

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

377 

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

379 maskchanged = maskchanged | self.runInteractiveGUI2() 

380 

381 time1 = time.time() 

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

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

384 return maskchanged 

385 

386 ############################################# 

387 def runInteractiveGUI2(self): 

388 maskchanged = False 

389 forcestop = True 

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

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

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

393 

394 for akey in self.stopMinor: 

395 if self.stopMinor[akey] < 0: 

396 maskchanged = True 

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

398 

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

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

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

402 forcestop = True 

403 for akey in self.stopMinor: 

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

405 

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

407 if forcestop == True: 

408 self.predictModel() 

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

410 # wstr = "Saving model column" 

411 # else: 

412 # wstr = "Saving virtual model" 

413 # 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") 

414 

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

416 return maskchanged or forcestop 

417 

418 ############################################# 

419 def makePSF(self): 

420 self.makePSFCore() 

421 divideInPython = ( 

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

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

424 ) 

425 

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

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

428 # for cube normalization is done in C++ 

429 if divideInPython: 

430 self.PStools[immod].gatherpsfweight() 

431 self.PStools[immod].dividepsfbyweight() 

432 self.check_psf(immod) 

433 

434 def check_psf(self, immod): 

435 if self.SDtools != []: 

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

437 self.SDtools[immod].checkrestoringbeam() 

438 

439 ############################################# 

440 def calcVisAppSens(self): 

441 

442 return self.SItool.apparentsens() 

443 

444 ############################################# 

445 

446 def runMajorCycle(self, isCleanCycle=True): 

447 # @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) 

448 if self.IBtool != None: 

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

450 else: 

451 lastcycle = True 

452 

453 divideInPython = ( 

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

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

456 ) 

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

458 if not divideInPython: 

459 self.runMajorCycleCore(lastcycle) 

460 if isCleanCycle: 

461 self.majorCnt += 1 

462 if self.IBtool != None: 

463 self.IBtool.endmajorcycle() 

464 return 

465 

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

467 self.PStools[immod].dividemodelbyweight() 

468 self.PStools[immod].scattermodel() 

469 self.runMajorCycleCore(lastcycle) 

470 if isCleanCycle: 

471 self.majorCnt += 1 

472 

473 if self.IBtool != None: 

474 self.IBtool.endmajorcycle() 

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

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

477 self.PStools[immod].gatherresidual() 

478 self.PStools[immod].divideresidualbyweight() 

479 self.PStools[immod].multiplymodelbyweight() 

480 

481 ############################################# 

482 def predictModel(self): 

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

484 self.PStools[immod].dividemodelbyweight() 

485 self.PStools[immod].scattermodel() 

486 

487 self.predictModelCore() 

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

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

490 self.PStools[immod].multiplymodelbyweight() 

491 

492 ## self.SItool.predictmodel(); 

493 

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

495 def dryGridding(self): 

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

497 

498 ############################################# 

499 ## Overloaded for parallel runs 

500 def fillCFCache(self): 

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

502 cflist = [] 

503 if not (cfcName == ""): 

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

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

506 self.cfcachepars["cflist"] = cflist 

507 

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

509 

510 self.SItool.fillcfcache( 

511 cflist, 

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

513 cfcName, 

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

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

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

517 ) 

518 

519 ############################################# 

520 def reloadCFCache(self): 

521 self.SItool.reloadcfcache() 

522 

523 ############################################# 

524 def makePB(self): 

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

526 if not ( 

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

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

529 ): 

530 self.makePBCore() 

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

532 self.PStools[immod].normalizeprimarybeam() 

533 

534 ############################################# 

535 def makePBCore(self): 

536 self.SItool.makepb() 

537 

538 ############################################# 

539 def checkPB(self): 

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

541 if self.SItool is None: 

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

543 return 

544 

545 import numpy as np 

546 

547 facetIdx = 0 # TODO iterate over facets 

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

549 _ia.open(imagename) 

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

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

552 pixelVals[1:-2][ 

553 1:-2 

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

555 if pixelVals.max() > 0: 

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

557 idx = [ 

558 x[0] for x in idx 

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

560 casalog.post( 

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

562 "WARN", 

563 ) 

564 # release the image 

565 _ia.close() 

566 _ia.done() 

567 

568 ############################################# 

569 def makeSdImage(self): 

570 self.makeSdImageCore() 

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

572 self.PStools[immod].gatherresidual() 

573 self.PStools[immod].divideresidualbyweight(singledish=True) 

574 

575 ############################################# 

576 def makeSdPSF(self): 

577 self.makeSdPSFCore() 

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

579 self.PStools[immod].gatherresidual() 

580 self.PStools[immod].dividepsfbyweight() 

581 

582 ############################################# 

583 def makeSdImageCore(self): 

584 self.SItool.makesdimage() 

585 

586 ############################################# 

587 def makeSdPSFCore(self): 

588 self.SItool.makesdpsf() 

589 

590 ############################################# 

591 def makeImage( 

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

593 ): 

594 """ 

595 This should replace makeSDImage, makeSDPSF and makePSF 

596 etc in the long run 

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

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

599 "singledish", "coverage", "holography", "holography-observed" 

600 For holography the FTmachine should be SDGrid and the baselines 

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

602 """ 

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

604 

605 ############################################# 

606 

607 ## Overloaded for parallel runs 

608 def setWeighting(self): 

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

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

611 ##moved the tuneselect after weighting so that 

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

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

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

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

616 # self.SItool.tuneselectdata() 

617 

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

619 # self.SItool.getweightdensity() 

620 # self.SItool.setweightdensity() 

621 

622 ############################################# 

623 ## Overloaded for parallel runs 

624 def makePSFCore(self): 

625 self.SItool.makepsf() 

626 

627 ############################################# 

628 ## Overloaded for parallel runs 

629 def runMajorCycleCore(self, lastcycle): 

630 controldict = {"lastcycle": lastcycle} 

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

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

633 self.SItool.executemajorcycle(controls=controldict) 

634 

635 ############################################# 

636 ## Overloaded for parallel runs 

637 def predictModelCore(self): 

638 self.SItool.predictmodel() 

639 

640 ############################################# 

641 

642 def runMinorCycle(self): 

643 return self.runMinorCycleCore() 

644 

645 ############################################# 

646 

647 def runMinorCycleCore(self): 

648 

649 # Set False for release packages. 

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

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

652 # currently does not work. 

653 # False = disable always save intermediate images mode 

654 alwaysSaveIntermediateImages = False 

655 

656 # Get iteration control parameters 

657 iterbotrec = self.IBtool.getminorcyclecontrols() 

658 

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

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

661 

662 self.IBtool.resetminorcycleinfo() 

663 

664 # 

665 # Run minor cycle 

666 self.ncycle += 1 

667 retval = False 

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

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

670 

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

672 if alwaysSaveIntermediateImages or ( 

673 "SAVE_ALL_RESIMS" in os.environ 

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

675 ): 

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

677 tempresname = ( 

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

679 + ".inputres" 

680 + str(self.ncycle) 

681 ) 

682 if os.path.isdir(resname): 

683 shutil.copytree(resname, tempresname) 

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

685 tempmodname = ( 

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

687 + ".inputmod" 

688 + str(self.ncycle) 

689 ) 

690 if os.path.isdir(modname): 

691 shutil.copytree(modname, tempmodname) 

692 

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

694 

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

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

697 self.IBtool.mergeexecrecord(exrec, immod) 

698 if alwaysSaveIntermediateImages or ( 

699 "SAVE_ALL_AUTOMASKS" in os.environ 

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

701 ): 

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

703 tempmaskname = ( 

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

705 + ".autothresh" 

706 + str(self.ncycle) 

707 ) 

708 if os.path.isdir(maskname): 

709 shutil.copytree(maskname, tempmaskname) 

710 

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

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

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

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

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

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

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

718 prevmaskname = ( 

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

720 ) 

721 if os.path.isdir(maskname): 

722 if os.path.isdir(prevmaskname): 

723 shutil.rmtree(prevmaskname) 

724 shutil.copytree(maskname, prevmaskname) 

725 return retval 

726 

727 ############################################# 

728 def runMajorMinorLoops(self): 

729 self.runMajorCycle(isCleanCycle=False) 

730 while not self.hasConverged(): 

731 self.runMinorCycle() 

732 self.runMajorCycle() 

733 

734 ############################################# 

735 

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

737 

738 if not ( 

739 "summaryminor" in summ 

740 and "summarymajor" in summ 

741 and "threshold" in summ 

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

743 ): 

744 casalog.post( 

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

746 ) 

747 return summ 

748 

749 import matplotlib.pyplot as pl 

750 from numpy import max as amax 

751 

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

753 # 1 : peak residual 

754 # 2 : model flux 

755 # 3 : cyclethreshold 

756 # 4 : deconvolver id 

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

758 

759 pl.ioff() 

760 

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

762 pl.clf() 

763 minarr = summ["summaryminor"] 

764 if minarr.size == 0: 

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

766 else: 

767 iterlist = minarr[0, :] 

768 eps = 0.0 

769 peakresstart = [] 

770 peakresend = [] 

771 modfluxstart = [] 

772 modfluxend = [] 

773 itercountstart = [] 

774 itercountend = [] 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

789 

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

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

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

793 

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

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

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

797 

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

799 pl.plot( 

800 itercountend, 

801 peakresend, 

802 "r.-", 

803 label="peak residual (end)", 

804 linewidth=2.5, 

805 ) 

806 pl.plot(itercountstart, modfluxstart, "b.--", label="model flux (start)") 

807 pl.plot( 

808 itercountend, modfluxend, "b.-", label="model flux (end)", linewidth=2.5 

809 ) 

810 pl.plot( 

811 iterlist, minarr[3, :], "g--", label="cycle threshold", linewidth=2.5 

812 ) 

813 maxval = amax(minarr[1, :]) 

814 maxval = max(maxval, amax(minarr[2, :]))