Coverage for /home/casatest/venv/lib/python3.12/site-packages/casatasks/private/task

1################################################

2# single dish + interfermeter join image reconstruction task

5################################################

7import platform

8import os

9import shutil

10import numpy

11import copy

12import time

14from casatasks import casalog

16from casatasks.private.imagerhelpers.imager_base import PySynthesisImager

17from casatasks.private.imagerhelpers.imager_parallel_continuum import PyParallelContSynthesisImager

18from casatasks.private.imagerhelpers.imager_parallel_cube import PyParallelCubeSynthesisImager

19from casatasks.private.imagerhelpers.input_parameters import ImagerParameters

20from .cleanhelper import write_tclean_history, get_func_params

21from .sdint_helper import *

22from casatools import table

23from casatools import synthesisimager,synthesisutils

25try:

26 from casampi.MPIEnvironment import MPIEnvironment

27 from casampi import MPIInterface

28 mpi_available = True

29except ImportError:

30 mpi_available = False

33# setup functions

34def setup_imagerObj(paramList=None):

35 """

36 setup imaging parameters

37 """

38 defaultconstructor = False

39 if paramList!=None:

40 if not isinstance(paramList, ImagerParameters):

41 raise RuntimeError("Internal Error: invalid paramList")

42 else:

43 defaultconstructor = True

45 if defaultconstructor:

46 return PySynthesisImager

47 else:

48 return PySynthesisImager(params=paramList)

51def setup_imager(imagename, specmode,calcres,calpsf,inparams):

52 """

53 Setup cube imaging for major cycles.

54 - Do initialization

55 - and run a major cycle

56 """

57 # create a local copy of input params dict so that it can be modified

58 locparams = copy.deepcopy(inparams)

60 # cube imaging setup

61 locparams['imagename']=imagename

62 locparams['specmode']='cube'

63 locparams['niter']=0

64 locparams['deconvolver']='hogbom'

66 #casalog.post("local inparams(msname) in setup_imager==",locparams['msname'])

67 params = ImagerParameters(**locparams)

69 ## Major cycle is either PySynthesisImager or PyParallelCubeSynthesisImager

70 imagertool = setup_imagerObj(params)

72 #self.imagertool = PySynthesisImager(params=params)

73 imagertool.initializeImagers()

74 imagertool.initializeNormalizers()

75 imagertool.setWeighting()

76 if 'psfphasecenter' in locparams:

77 psfphasecenter = locparams['psfphasecenter']

78 else:

79 psfphasecenter = ''

81 ## Extra one for psfphasecenter...

82 imagerInst=None

83 if((psfphasecenter != '') and (gridder=='mosaic')):

84 imagerInst = setup_imagerObj()

87 gridder = locparams['gridder']

89 if calpsf == True:

90 imagertool.makePSF()

91 imagertool.makePB()

92 if((psfphasecenter != '') and (gridder=='mosaic')):

93 casalog.post("doing with different phasecenter psf", "INFO")

94 imagertool.unlockimages(0)

95 psfParameters=paramList.getAllPars()

96 psfParameters['phasecenter']=psfphasecenter

97 psfParamList=ImagerParameters(**psfParameters)

98 psfimager=imagerInst(params=psfParamList)

99 psfimager.initializeImagers()

100 psfimager.setWeighting()

101 psfimager.makeImage('psf', psfParameters['imagename']+'.psf')

102

103 # can take out this since niter is fixed to 0

104 if locparams['niter'] >=0 :

105 ## Make dirty image

106 if calcres == True:

107 t0=time.time();

108 imagertool.runMajorCycle(isCleanCycle=False)

109 t1=time.time();

110 casalog.post("***Time for major cycle (calcres=T): "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_tclean");

111

112 ## In case of no deconvolution iterations....

113 #if locparams['niter']==0 and calcres==False:

114 # if savemodel != "none":

115 # imagertool.predictModel()

116

117 return imagertool

118

119def setup_deconvolver(imagename,specmode,inparams):

120 """

121 Cube or MFS minor cycles.

122 """

123 inparams['imagename']=imagename

124 params = ImagerParameters(**inparams)

125 deconvolvertool = setup_imagerObj(params)

126

127 ## Why are we initializing these ?

128 deconvolvertool.initializeImagers()

129 deconvolvertool.initializeNormalizers()

130 deconvolvertool.setWeighting()

131

132

133 ### These three should be unncessary. Need a 'makeimage' method for csys generation.

134 deconvolvertool.makePSF() ## Make this to get a coordinate system

135 #deconvolvertool.makeImage('psf', imagename+'.psf')

136 deconvolvertool.makePB() ## Make this to turn .weight into .pb maps

137

138 ## Initialize deconvolvers. ( Order is important. This cleans up a leftover tablecache image.... FIX!)

139 deconvolvertool.initializeDeconvolvers()

140 deconvolvertool.initializeIterationControl() # This needs to be run before runMajorCycle

141 deconvolvertool.runMajorCycle(isCleanCycle=False) ## Make this to make template residual images.

142

143 return deconvolvertool

144

145def setup_sdimaging(template='',output='', inparms=None, sdparms=None):

146 """

147 Make the SD cube Image and PSF

148

149 Option 1 : Use/Regrid cubes for the observed image and PSF

150 Option 2 : Make the SD image and PSF cubes using 'tsdimager's usage of the SD gridder option.

151

152 Currently, only Option 1 is supported.

153

154 """

155 sdintlib = SDINT_helper()

156 if 'sdpsf' in sdparms:

157 sdpsf = sdparms['sdpsf']

158 else:

159 raise RuntimeError("Internal Error: missing sdpsf parameter")

160

161 if 'sdimage' in sdparms:

162 sdimage = sdparms['sdimage']

163 else:

164 raise RuntimeError("Internal Error: missing sdimage parameter")

165 if 'pblimit' in inparms:

166 pblimit = inparms['pblimit']

167

168 if sdpsf !="":

169 ## check the coordinates of psf with int psf

170 sdintlib.checkpsf(sdpsf, template+'.psf')

171

172 ## Regrid the input SD image and PSF cubes to the target coordinate system.

173 sdintlib.regridimage(imagename=sdimage, template=template+'.residual', outfile=output+'.residual')

174 sdintlib.regridimage(imagename=sdimage, template=template+'.residual', outfile=output+'.image')

175

176 if sdpsf !="":

177 sdintlib.regridimage(imagename=sdpsf, template=template+'.psf', outfile=output+'.psf')

178 else:

179 ## Make an internal sdpsf image if the user has not supplied one.

180 casalog.post("Constructing a SD PSF cube by evaluating Gaussians based on the restoring beam information in the regridded SD Image Cube")

181 sdintlib.create_sd_psf(sdimage=output+'.residual', sdpsfname=output+'.psf')

182

183 ## Apply the pbmask from the INT image cube, to the SD cubes.

184 #TTB: Create *.mask cube

185

186 sdintlib.addmask(inpimage=output+'.residual', pbimage=template+'.pb', pblimit=pblimit)

187 sdintlib.addmask(inpimage=output+'.image', pbimage=template+'.pb', pblimit=pblimit)

188

189 sdintlib.deleteTmpFiles()

190

191

192

193def sdintimaging(

194 usedata,

195 ####### Single dish input data

196 sdimage,

197 sdpsf,

198 sdgain,

199 dishdia,

200 ####### Interferometer Data Selection

201 vis,#='',

202 selectdata,

203 field,#='',

204 spw,#='',

205 timerange,#='',

206 uvrange,#='',

207 antenna,#='',

208 scan,#='',

209 observation,#='',

210 intent,#='',

211 datacolumn,#='corrected',

212

213

214 ####### Image definition

215 imagename,#='',

216 imsize,#=[100,100],

217 cell,#=['1.0arcsec','1.0arcsec'],

218 phasecenter,#='J2000 19:59:28.500 +40.44.01.50',

219 stokes,#='I',

220 projection,#='SIN',

221 startmodel,#='',

222

223 ## Spectral parameters

224 specmode,#='mfs',

225 reffreq,#='',

226 nchan,#=1,

227 start,#='',

228 width,#='',

229 outframe,#='LSRK',

230 veltype,#='',

231 restfreq,#=[''],

232# sysvel,#='',

233# sysvelframe,#='',

234 interpolation,#='',

235# chanchunks,#=1,

236 perchanweightdensity, #=''

237 ##

238 ####### Gridding parameters

239 gridder,#='ft',

240 facets,#=1,

241 psfphasecenter,#='',

242

243 wprojplanes,#=1,

244

245 ### PB

246 vptable,

247 mosweight, #=True

248 aterm,#=True,

249 psterm,#=True,

250 wbawp ,#= True,

251# conjbeams ,#= True,

252 cfcache ,#= "",

253 usepointing, #=false

254 computepastep ,#=360.0,

255 rotatepastep ,#=360.0,

256 pointingoffsetsigdev ,#=0.0,

257

258 pblimit,#=0.01,

259# normtype,#='flatnoise',

260

261 ####### Deconvolution parameters

262 deconvolver,#='hogbom',

263 scales,#=[],

264 nterms,#=1,

265 smallscalebias,#=0.0

266

267 ### restoration options

268 restoration,

269 restoringbeam,#=[],

270 pbcor,

271

272 ##### Outliers

273# outlierfile,#='', ### RESTRICTION : No support for outlier fields for joint SD-INT imaging.

274

275 ##### Weighting

276 weighting,#='natural',

277 robust,#=0.5,

278 noise,#0.0Jy

279 npixels,#=0,

280# uvtaper,#=False,

281 uvtaper,#=[],

282

283

284 ##### Iteration control

285 niter,#=0,

286 gain,#=0.1,

287 threshold,#=0.0,

288 nsigma,#=0.0

289 cycleniter,#=0,

290 cyclefactor,#=1.0,

291 minpsffraction,#=0.1,

292 maxpsffraction,#=0.8,

293 interactive,#=False,

294 fullsummary,#=False,

295 nmajor,#=-1,

296

297 ##### (new) Mask parameters

298 usemask,#='user',

299 mask,#='',

300 pbmask,#='',

301 # maskthreshold,#='',

302 # maskresolution,#='',

303 # nmask,#=0,

304

305 ##### automask by multithresh

306 sidelobethreshold,#=5.0,

307 noisethreshold,#=3.0,

308 lownoisethreshold,#=3.0,

309 negativethreshold,#=0.0,

310 smoothfactor,#=1.0,

311 minbeamfrac,#=0.3,

312 cutthreshold,#=0.01,

313 growiterations,#=100

314 dogrowprune,#=True

315 minpercentchange,#=0.0

316 verbose, #=False

317 fastnoise, #=False

318

319 ## Misc

320

321 restart,#=True,

322

323 #savemodel,#="none",

324

325# makeimages,#="auto"

326 calcres,#=True,

327 calcpsf):#=True,

328

329 ####### State parameters

330 #parallel):#=False):

331

332

333 ##################################################

334 # copied from SDINT.do_reconstruct

335 #################################################

336 int_cube = imagename+'.int.cube'

337 sd_cube = imagename+'.sd.cube'

338 joint_cube = imagename+'.joint.cube'

339 joint_multiterm = imagename+'.joint.multiterm'

340

341 if specmode=='mfs':

342 decname = joint_multiterm

343 else:

344 decname = joint_cube

345

346 #####################################################

347 #### Sanity checks and controls

348 #####################################################

349

350 if interactive:

351 # Check for casaviewer, if it does not exist flag it up front for macOS

352 # since casaviewer is no longer provided by default with macOS.

353 try:

354 import casaviewer as __test_casaviewer

355 except:

356 if platform.system( ) == "Darwin":

357 casalog.post(

358 "casaviewer is no longer available for macOS, for more information see: http://go.nrao.edu/casa-viewer-eol Please restart by setting interactive=F",

359 "WARN",

360 "task_sdintimaging",

361 )

362 raise RuntimeError( "casaviewer is no longer available for macOS, for more information see: http://go.nrao.edu/casa-viewer-eol" )

363

364 ### Move these checks elsewhere ?

365 inpparams=locals().copy()

366 ###now deal with parameters which are not the same name

367 #casalog.post("current inpparams=",inpparams)

368 #casalog.post("inpparams.keys()=",inpparams.keys())

369 locvis=inpparams.pop('vis')

370 #casalog.post("LOCVIS====",locvis)

371 if type(locvis)==list:

372 llocvis = [v.lstrip() for v in locvis]

373 else:

374 llocvis = locvis.lstrip()

375 inpparams['msname']=llocvis

376 inpparams['timestr']= inpparams.pop('timerange')

377 inpparams['uvdist']= inpparams.pop('uvrange')

378 inpparams['obs']= inpparams.pop('observation')

379 inpparams['state']= inpparams.pop('intent')

380 inpparams['loopgain']=inpparams.pop('gain')

381 inpparams['scalebias']=inpparams.pop('smallscalebias')

382

383 sdparms={}

384 sdparms['sdimage']=inpparams['sdimage']

385 sdparms['sdpsf']=inpparams['sdpsf']

386 sdparms['sdgain']=inpparams['sdgain']

387

388 if usedata!='int': # check sd parameters

389

390 _myia = image()

391

392 if not os.path.exists(sdparms['sdimage']):

393 casalog.post( "Input image sdimage = '"+str(sdparms['sdimage'])+"' does not exist.", "WARN", "task_sdintimaging" )

394 return

395 else:

396 try:

397 _myia.open(sdparms['sdimage'])

398 except Exception as instance:

399 casalog.post( "Input image sdimage = '"+str(sdparms['sdimage'])+"' cannot be opened.", "WARN", "task_sdintimaging" )

400 casalog.post( str(instance), "WARN", "task_sdintimaging" )

401 return

402

403 mysummary = _myia.summary(list=False)

404 _myia.close()

405

406 try:

407 freqaxis_index = list(mysummary['axisnames']).index('Frequency')

408 except(ValueError):

409 casalog.post('The image '+sdparms['sdimage']+' has no frequency axis. Try adding one with ia.adddegaxis() .',

410 'WARN', 'task_sdintimaging')

411 return

412

413 if freqaxis_index != 3:

414 casalog.post('The image '+sdparms['sdimage']+' has its frequency axis on position '+str(freqaxis_index)+

415 ' whereas it should be in position 3 (counting from 0). Use task imtrans() with order=["r", "d", "s", "f"] to fix this.',

416 'WARN', 'task_sdintimaging')

417 return

418

419 if stokes!='I' and usedata=='sdint':

420 casalog.post('You have specified parameter stokes=\"'+str(stokes)+'\" but presently only stokes=\"I\" is supported when usedata=\"sdint\".',

421 'WARN', 'task_sdintimaging')

422 return

423

424

425 if sdparms['sdpsf']!='':

426 if not os.path.exists(sdparms['sdpsf']):

427 casalog.post( "Input image sdpsf = '"+str(sdparms['sdpsf'])+"' does not exist.", "WARN", "task_sdintimaging" )

428 return

429 else:

430 try:

431 _myia.open(sdparms['sdpsf'])

432 _myia.close()

433 except Exception as instance:

434 casalog.post( "Input image sdpsf = '"+str(sdparms['sdpsf'])+"' cannot be opened.", "WARN", "task_sdintimaging" )

435 casalog.post( str(instance), "WARN", "task_sdintimaging" )

436 return

437

438 if (sdparms['sdgain']*0!=0 or sdparms['sdgain']<=0):

439 casalog.post('Invalid sdgain: '+str(sdparms['sdgain']), 'WARN')

440 casalog.post("The sdgain parameter needs to be chosen as a number > 0 which represents the weight of the SD contribution relative to the INT contribution to the joint image.", "WARN", "task_sdintimaging")

441 return

442

443 if (dishdia*0!=0 or dishdia<=0):

444 casalog.post('Invalid dishdia: '+str(dishdia), 'WARN')

445 casalog.post("The dishdia parameter needs to provide the diameter (meters) of the SD telescope which produced the SD image.", "WARN", "task_sdintimaging")

446 return

447

448

449 if specmode=='cont':

450 specmode='mfs'

451 inpparams['specmode']='mfs'

452

453 # from sdint

454 # automatically decide if pb need to be applied

455 if gridder=='mosaic' or gridder=='awproject':

456 applypb = True

457 else:

458 applypb = False

459

460 if (deconvolver=="mtmfs") and (specmode!='mfs') and (specmode!='cube' or nterms!=1) and (specmode!='cubedata' or nterms!=1):

461 casalog.post( "The MSMFS algorithm (deconvolver='mtmfs') applies only to specmode='mfs' or specmode='cube' with nterms=1 or specmode='cubedata' with nterms=1.", "WARN", "task_sdintimaging" )

462 return

463

464 if(deconvolver=="mtmfs" and (specmode=='cube' or specmode=='cubedata') and nterms==1 ):

465 casalog.post( "The MSMFS algorithm (deconvolver='mtmfs') with specmode='cube', nterms=1 is currently not supported. Please use deconvolver='multiscale' instead for cubes.", "WARN", "task_sdintimaging" )

466 return

467

468 if(specmode=='mfs' and deconvolver!='mtmfs'):

469 casalog.post("Currently, only the multi-term MFS algorithm is supported for specmode=mfs. To make a single plane MFS image (while retaining the frequency dependence for the cube major cycle stage), please pick nterms=1 along with deconvolver=mtmfs. The scales parameter is still usable for multi-scale multi-term deconvolution","WARN","task_sdintimaging")

470 return;

471

472 if(usedata=='sd'):

473 casalog.post("The Single-Dish-Only mode of sdintimaging is better supported via the deconvolve task which supports spectral cube, mfs and multi-term mfs deconvolution in the image domain alone. The deconvolve task is the more appropriate version to use for stand-alone image-domain deconvolution, and will not have the bookkeeping overheads currently present in the sdintimaging task's sd-only mode. Please note that the 'sd' option of the sdintimaging task will be removed in a subsequent release. Please refer to the task deconvolve documentation for instructions on how to prepare image and psf cubes for the deconvolve task for all these modes.","WARN","task_sdintimaging");

474

475 if (nmajor < -1):

476 casalog.post("Negative values less than -1 for nmajor are reserved for possible future implementation", "WARN", "task_sdintimaging")

477 return

478

479# if parallel==True:

480# casalog.post("Cube parallelization (all major cycles) is currently not supported via task_sdintimaging. This will be enabled after a cube parallelization rework.")

481# return;

482

483 #####################################################

484 #### Construct ImagerParameters object

485 #####################################################

486

487 imager = None

488 paramList = None

489 deconvolvertool = None

490

491 # Put all parameters into dictionaries and check them.

492 ##make a dictionary of parameters that ImagerParameters take

493 defparm=dict(list(zip(ImagerParameters.__init__.__code__.co_varnames[1:], ImagerParameters.__init__.__defaults__)))

494

495

496 ###assign values to the ones passed to tclean and if not defined yet in tclean...

497 ###assign them the default value of the constructor

498 bparm={k: inpparams[k] if k in inpparams else defparm[k] for k in defparm.keys()}

499

500 ###default mosweight=True is tripping other gridders as they are not

501 ###expecting it to be true

502 if(bparm['mosweight']==True and bparm['gridder'].find("mosaic") == -1):

503 bparm['mosweight']=False

504

505 ## Two options have been removed from the interface. Hard-code them here.

506 bparm['normtype'] = 'flatnoise' ## Hard-code this since the pbcor steps assume it.

507 bparm['conjbeams']=False

508

509 #paramList=ImagerParameters(**bparm)

510

511 #paramList.printParameters()

512

513 if len(pointingoffsetsigdev)>0 and pointingoffsetsigdev[0]!=0.0 and usepointing==True and gridder.count('awproj')>1:

514 casalog.post("pointingoffsetsigdev will be used for pointing corrections with AWProjection", "WARN")

515

516 #=========================================================

517 ####set the children to load c++ libraries and applicator

518 ### make workers ready for c++ based mpicommands

519 cppparallel=False

520 if mpi_available and MPIEnvironment.is_mpi_enabled:

521 mint=MPIInterface.MPIInterface()

522 cl=mint.getCluster()

523 cl._cluster.pgc("from casatools import synthesisimager", False)

524 cl._cluster.pgc("si=synthesisimager()", False)

525

526 cl._cluster.pgc("si.initmpi()", False)

527 cppparallel=True

528 ###ignore chanchunk

529 bparm['chanchunks']=1

530

531 #################################################

532 #### start of more computing-intensive work #####

533 #################################################

534

535 synu = synthesisutils()

536

537 retrec={}

538

539 try:

540 mysdintlib = SDINT_helper()

541 ## Init major cycle elements

542 casalog.post("INT cube setup ....")

543 t0=time.time();

544 imager=setup_imager(int_cube, specmode, calcres, calcpsf, bparm)

545 mysdintlib.copy_restoringbeam(fromthis=int_cube+'.psf', tothis=int_cube+'.residual')

546

547 t1=time.time();

548 casalog.post("***Time for initializing imager (INT cube) : "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_sdintimaging");

549

550 ## Init minor cycle elements

551 if niter>0 or restoration==True:

552 casalog.post("Combined image setup ....")

553 t0=time.time();

554 deconvolvertool=setup_deconvolver(decname, specmode, bparm )

555

556 t1=time.time();

557 casalog.post("***Time for seting up deconvolver(s): "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_sdintimaging");

558

559 if usedata!='int':

560 ### debug (remove it later)

561 casalog.post("SD cube setup ....")

562 setup_sdimaging(template=int_cube, output=sd_cube, inparms=bparm, sdparms=sdparms )

563

564

565 ####now is the time to check estimated memory

566 # need to move to somewhere below???

567 imager.estimatememory()

568

569 ## Do sanity checks on INT and SD cubes

570 ### Frequency range of cube, data selection range. mtmfs reffreq.

571 ### nchan too small or too large

572 ### sumwt : flagged channels in int cubes

573 ### sd cube empty channels ? Weight image ?

574 validity, inpparams = mysdintlib.check_coords(intres=int_cube+'.residual', intpsf=int_cube+'.psf',

575 intwt = int_cube+'.sumwt',

576 sdres=sd_cube+'.residual', sdpsf=sd_cube+'.psf',

577 sdwt = '',

578 pars=inpparams)

579

580 if validity==False:

581 casalog.post('Exiting from the sdintimaging task due to inconsistencies found between the interferometer-only and singledish-only image and psf cubes. Please modify inputs as needed','WARN')

582 if imager != None:

583 imager.deleteTools()

584 if deconvolvertool != None:

585 deconvolvertool.deleteTools()

586 mysdintlib.deleteTmpFiles()

587 return

588

589 #### SDINT specific feathering....

590 ## Feather INT and SD residual images (feather in flat-sky. output has common PB)

591 casalog.post("Feathering INT and SD residual images...")

592 mysdintlib.feather_residual(int_cube, sd_cube, joint_cube, applypb, inpparams)

593 mysdintlib.feather_int_sd(sdcube=sd_cube+'.psf',

594 intcube=int_cube+'.psf',

595 jointcube=joint_cube+'.psf',

596 sdgain=sdgain,

597 dishdia=dishdia,

598 usedata=usedata,

599 chanwt = inpparams['chanwt'])

600

601 #print("Fitting for cube")

602 synu.fitPsfBeam(joint_cube)

603

604 ###############

605 ##### Placeholder code for PSF renormalization if needed

606 ##### Note : If this is enabled, we'll need to restrict the use of 'faceting' as .sumwt shape changes.

607 #mysdintlib.calc_renorm(intname=int_cube, jointname=joint_cube)

608 #mysdintlib.apply_renorm(imname=joint_cube+'.psf', sumwtname=joint_cube+'.sumwt')

609 #mysdintlib.apply_renorm(imname=joint_cube+'.residual', sumwtname=joint_cube+'.sumwt')

610 ###############

611

612 #casalog.post("feather_int_sd DONE")

613

614 if specmode=='mfs':

615 ## Calculate Spectral PSFs and Taylor Residuals

616 casalog.post("Calculate spectral PSFs and Taylor Residuals...")

617 mysdintlib.cube_to_taylor_sum(cubename=joint_cube+'.psf',

618 cubewt=int_cube+'.sumwt',

619 chanwt=inpparams['chanwt'],

620 mtname=joint_multiterm+'.psf',

621 nterms=nterms, reffreq=inpparams['reffreq'], dopsf=True)

622 mysdintlib.cube_to_taylor_sum(cubename=joint_cube+'.residual',

623 cubewt=int_cube+'.sumwt',

624 chanwt=inpparams['chanwt'],

625 mtname=joint_multiterm+'.residual',

626 nterms=nterms, reffreq=inpparams['reffreq'], dopsf=False)

627

628 #print("Fit for multiterm")

629 if(deconvolver=='mtmfs' and nterms==1): # work around file naming issue

630 os.system('rm -rf '+joint_multiterm+'tmp.psf')

631 os.system('ln -sf '+joint_multiterm+'.psf.tt0 '+joint_multiterm+'tmp.psf')

632 synu.fitPsfBeam(joint_multiterm+'tmp',nterms=nterms)

633 os.system('rm -rf '+joint_multiterm+'tmp.psf')

634 else:

635 synu.fitPsfBeam(joint_multiterm,nterms=nterms)

636

637 if niter>0 :

638 isit = deconvolvertool.hasConverged()

639 deconvolvertool.updateMask()

640

641 while ( not deconvolvertool.hasConverged() ):

642

643 t0=time.time();

644

645 #### Print out the peak of the residual image here to check !!!

646# if specmode=='mfs':

647# print('Max of joint residual before initminorcycle' + str(imstat(joint_multiterm+'.residual.tt0',verbose=False)['max'][0]))

648# else:

649# print('Max of joint residual before initminorcycle' + str(imstat(joint_cube+'.residual',verbose=False)['max'][0]))

650

651

652

653 deconvolvertool.runMinorCycle()

654

655# if specmode=='mfs':

656# print('Max of joint residual after minorcycle' + str(imstat(joint_multiterm+'.residual.tt0',verbose=False)['max'][0]))

657# else:

658# print('Max of joint residual after minorcycle' + str(imstat(joint_cube+'.residual',verbose=False)['max'][0]))

659

660

661 t1=time.time();

662 casalog.post("***Time for minor cycle: "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_sdintimaging");

663

664 ### sdint specific feathering steps HERE

665 ## Prepare the joint model cube for INT and SD major cycles

666 if specmode=='mfs':

667 ## Convert Taylor model coefficients into a model cube : int_cube.model

668 mysdintlib.taylor_model_to_cube(cubename=int_cube, ## output

669 mtname=joint_multiterm, ## input

670 nterms=nterms, reffreq=inpparams['reffreq'])

671 else:

672 ## Copy the joint_model cube to the int_cube.model

673 shutil.rmtree(int_cube+'.model',ignore_errors=True)

674 shutil.copytree(joint_cube+'.model', int_cube+'.model')

675 hasfile=os.path.exists(joint_cube+'.model')

676 #casalog.post("DEBUG: has joint cube .image===",hasfile)

677

678 if applypb==True:

679 ## Take the int_cube.model to flat sky.

680 if specmode=='cube':

681 ## Divide the model by the frequency-dependent PB to get to flat-sky

682 fdep_pb = True

683 else:

684 ## Divide the model by the common PB to get to get to flat-sky

685 fdep_pb = False

686 mysdintlib.modify_with_pb(inpcube=int_cube+'.model',

687 pbcube=int_cube+'.pb',

688 cubewt=int_cube+'.sumwt',

689 chanwt=inpparams['chanwt'],

690 action='div', pblimit=pblimit,freqdep=fdep_pb)

691

692 if usedata!="int":

693 ## copy the int_cube.model to the sd_cube.model

694 shutil.rmtree(sd_cube+'.model',ignore_errors=True)

695 shutil.copytree(int_cube+'.model', sd_cube+'.model')

696

697 if applypb==True:

698 ## Multiply flat-sky model with freq-dep PB

699 mysdintlib.modify_with_pb(inpcube=int_cube+'.model',

700 pbcube=int_cube+'.pb',

701 cubewt=int_cube+'.sumwt',

702 chanwt=inpparams['chanwt'],

703 action='mult', pblimit=pblimit, freqdep=True)

704

705 ## Major cycle for interferometer data

706 t0=time.time();

707 # print('Max of int residual before major cycle' + str(imstat(int_cube+'.residual',verbose=False)['max'][0]))

708 # print('Max of int model before major cycle' + str(imstat(int_cube+'.model',verbose=False)['max'][0]))

709

710 if usedata != "sd":

711 imager.runMajorCycle()

712 # track nmajor for the deconvolvertool.hasConverged() method

713 deconvolvertool.majorCnt = imager.majorCnt

714

715 # print('Max of int residual after major cycle' + str(imstat(int_cube+'.residual',verbose=False)['max'][0]))

716 t1=time.time();

717 casalog.post("***Time for major cycle: "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_tclean");

718

719 if usedata!="int":

720 ## Major cycle for Single Dish data (uses the flat sky cube model in sd_cube.model )

721 mysdintlib.calc_sd_residual(origcube=sd_cube+'.image',

722 modelcube=sd_cube+'.model',

723 residualcube=sd_cube+'.residual', ## output

724 psfcube=sd_cube+'.psf')

725

726 ## Feather the residuals

727 mysdintlib.feather_residual(int_cube, sd_cube, joint_cube, applypb, inpparams)

728 ###############

729 ##### Placeholder code for PSF renormalization if needed

730 #mysdintlib.apply_renorm(imname=joint_cube+'.residual', sumwtname=joint_cube+'.sumwt')

731 ###############

732

733 if specmode=='mfs':

734 ## Calculate Spectral Taylor Residuals

735 mysdintlib.cube_to_taylor_sum(cubename=joint_cube+'.residual',

736 cubewt=int_cube+'.sumwt',

737 chanwt=inpparams['chanwt'],

738 mtname=joint_multiterm+'.residual',

739 nterms=nterms, reffreq=inpparams['reffreq'], dopsf=False)

740

741# if specmode=='mfs':

742# print('Max of residual after feather step ' + str(imstat(joint_multiterm+'.residual.tt0',verbose=False)['max'][0]))

743# else:

744# print('Max of residual after feather step ' + str(imstat(joint_cube+'.residual',verbose=False)['max'][0]))

745

746

747 deconvolvertool.updateMask()

748

749 ## Get summary from iterbot

750 #if type(interactive) != bool:

751 #retrec=imager.getSummary();

752 retrec=deconvolvertool.getSummary(fullsummary);

753 retrec['nmajordone'] = imager.majorCnt

754 if calcres==True:

755 retrec['nmajordone'] = retrec['nmajordone'] + 1 ## To be consistent with tclean. Remove, when we can change the meaning of nmajordone to exclude the initial major cycles.

756

757 ## Restore images.

758 if restoration==True:

759 t0=time.time();

760 deconvolvertool.restoreImages()

761 t1=time.time();

762 casalog.post("***Time for restoring images: "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_tclean");

763 if pbcor==True:

764 #if applypb==True:

765 t0=time.time();

766 if specmode=='mfs':

767 mysdintlib.pbcor(imagename=decname+'.image.tt0' , pbimage=decname+'.pb.tt0' , cutoff=pblimit,outfile=decname+'.image.tt0.pbcor')

768 else:

769 mysdintlib.pbcor(imagename=joint_cube+'.image' , pbimage=int_cube+'.pb' , cutoff=pblimit,outfile=joint_cube+'.image.pbcor')

770

771 #imager.pbcorImages()

772 t1=time.time();

773 casalog.post("***Time for pb-correcting images: "+"%.2f"%(t1-t0)+" sec", "INFO3", "task_tclean");

774

775 ##close tools

776 # needs to deletools before concat or lock waits for ever

777 imager.deleteTools()

778 deconvolvertool.deleteTools()

779

780

781 finally:

782 if imager != None:

783 imager.deleteTools()

784 if(cppparallel):

785 ###release workers back to python mpi control

786 si=synthesisimager()

787 si.releasempi()

788

789 #clean up tmp files

790 mysdintlib.deleteTmpFiles()

791

792 # Write history at the end, when hopefully all temp files are gone from disk,

793 # so they won't be picked up. They need time to disappear on NFS or slow hw.

794 # Copied from tclean.

795 try:

796 params = get_func_params(sdintimaging, locals())

797 write_tclean_history(imagename, 'sdintimaging', params, casalog)

798 except Exception as exc:

799 casalog.post("Error updating history (logtable): {} ".format(exc),'WARN')

800

801 return retrec

802

803##################################################

Coverage for /home/casatest/venv/lib/python3.12/site-packages/casatasks/private/task_sdintimaging.py: 7%

324 statements