← Back to branch summary

~i-martividal/uvmultifit/trunk-1

« back to all changes in this revision

Viewing changes to uvmultifit.py

  • Committer: I. Marti-Vidal
  • Date: 2020-10-02 17:18:30 UTC
  • Revision ID: i.marti-vidal@uv.es-20201002171830-d32ffzvcbp1y0sue
Updated the documentation with new Closure Test

Show diffs side-by-side

added added

removed removed

Lines of Context:
uvmultifit.py
411
411
- **EXAMPLE 8**: Simultaneous self-calibration and source fitting. We fit a source flux 
412
412
  density and perform Global Fringe Fitting (i.e., we fit for phases, delays, and delay
413
413
  rates) in one shot. Let's fit stokes RR and assume that there are 3 antennas (for 
414
 
  simplicity):
 
414
  simplicity), where "EB" and "FD" are going to be referred to, say, "LA":
415
415
 
416
416
  >>> stokes = 'RR'
417
417
  >>> model = ['delta']
418
 
  >>> phase_gains = {1:'p[0] + 6.2832*(p[1]*(nu-nu0) + p[2]*t)',
419
 
  >>>                2:'p[3] + 6.2832*(p[4]*(nu-nu0) + p[5]*t)'}
 
418
  >>> phase_gains = {'nuG': {'EB':'p[0] + 6.2832*p[1]*(nu-nu0)',
 
419
  >>>                        'FD':'p[3] + 6.2832*p[4]*(nu-nu0)'},
 
420
  >>>                 'tG': {'EB':'p[2]*t', 'FD':'p[5]*t'}}
 
421
  >>>
420
422
  >>> var = '0,0,p[6]'
421
423
 
422
 
  Good initial estimates for the phase gains can be obtained using the "QuinnFF" method.
 
424
  Notice that "nuG" contains the part of the gains that depends on frequency ("nu"; being "nu0"
 
425
  the lowest frequency in the dataset) and "tG" contains the part of the gains that depends on
 
426
  time. Good initial estimates for the phase gains can be obtained using the "QuinnFF" method.
 
427
 
 
428
|
 
429
|
 
430
 
 
431
- **EXAMPLE 9**: Fit of a dataset with corrupted phases (and amplitudes) to a simple model 
 
432
  source with non-zero closures. If the structure is a double source:
 
433
 
 
434
  >>> model = ['delta','delta']
 
435
  >>> var = ['0,0,1','p[1],p[2],p[3]']
 
436
  >>> phase_closure_Wgt = -1.0
 
437
  >>> amp_closure_Wgt = -1.0
 
438
 
 
439
  In this case, we **only** want to use the closure quantities in the fit (since the antenna
 
440
  gains have been corrupted), so we set "phase_closure_Wgt" and "amp_closure_Wgt" both to a 
 
441
  *negative number*. If we would like to use closures **and** visibilities simultaneously, we 
 
442
  should set these parameters to positive numbers (which would be the *relative weights between
 
443
  visibilities and closures*). Notice that (since the closures are insenstive to the absolute 
 
444
  flux scale and the absolute source position), one of the components has no position offset 
 
445
  and a normalized flux density. Hence, the parameters of the second source will be referred 
 
446
  to those of the first one.
 
447
 
 
448
 
423
449
 
424
450
 
425
451
=============
589
615
import gc
590
616
import re
591
617
import shutil
592
 
if False:
 
618
 
 
619
## Set to False when regenerating the documentation:
 
620
if True:
 
621
 
593
622
  from taskinit import gentools
594
623
  from clearcal_cli import clearcal_cli as clearcal
595
624
  ms = gentools(['ms'])[0]
871
900
            UVMultiFit will only use closure information in the fit (i.e., the 
872
901
            contribution from the visibilities will be turned off).
873
902
 
874
 
**trispec_refants** : `list'
 
903
**trispec_refants** : `list`
875
904
  List with the indiced of the antennas for which we have a reliable and robust amplitude
876
905
  calibration. If the phase closures are being used in the fit, UVMultiFit will add the 
877
906
  contribution of the amplitudes of the trispectra associated to these antennas. This 
4217
4246
 
4218
4247
    sys.stdout.write('\n Initial Chi Square: %.4e\n'%CurrChi) ; sys.stdout.flush()
4219
4248
 
4220
 
 #   HOLD = raw_input('HOLD')
 
4249
#    HOLD = raw_input('HOLD')
4221
4250
 
4222
4251
 
4223
4252
    controlIter = 0
4336
4365
    if self.doClos and not self.onlyClos:
4337
4366
      self.residuals(self.par2[2,:],-1,dof=False,doClos=False,onlyClos=False)
4338
4367
    elif self.onlyClos:
4339
 
      self.residuals(self.par2[2,:],-1,dof=False,doClos=True,onlyClos=True)
 
4368
      self.residuals(self.par2[2,:],-1,dof=False,doClos=True,onlyClos=True,doDump=True)
4340
4369
 
4341
4370
 
4342
4371
  #  self.getPar2(mode=1)
4435
4464
#  COMPUTE RESIDUALS FOR A MODEL REALIZATION
4436
4465
#
4437
4466
 
4438
 
  def residuals(self,p,mode=-1,dof=True,doClos=False,onlyClos=True):
 
4467
  def residuals(self,p,mode=-1,dof=True,doClos=False,onlyClos=True,doDump=False):
4439
4468
    """ Compute the residuals, fixed model, and/or covariance matrix and Chi square.
4440
4469
 
4441
4470
This method has a wide and flexible usage, depending on the value of 'mode'
4456
4485
  Whether closure quantities are taken into account to compute the Hessian and Gradient.
4457
4486
onlyClos : `bool`
4458
4487
  Whether ONLY closure quantities are used.
 
4488
doDump : `bool`
 
4489
  If True, writes the observed and (final) model closure quantities (if they are used 
 
4490
  in the fit) to a file.
4459
4491
 
4460
4492
The so-called 'output array' is the data that will be saved into the measurement set(s)
4461
4493
then the "writeModel" method of the parent UVMultiFit instance is called."""
4541
4573
           self.varbuffer[0][0,i,:nnu] = tempvar[i]
4542
4574
      #     print self.varbuffer[0][0,i,:nnu]
4543
4575
 
4544
 
         self.Ccompmodel(spw,nui,0,0)
 
4576
         self.Ccompmodel(spw,nui,0,0,0,0)
4545
4577
 
4546
4578
     self.imod[0] = modbackup
4547
4579
     return 0
4671
4703
           for i in range(len(tempvar)): 
4672
4704
             self.varbuffer[j+1][midx,i,0] = tempvar[i]
4673
4705
 
4674
 
     ChiSq += self.Ccompmodel(spw,nui,mode,doClos,onlyClos)
 
4706
     ChiSq += self.Ccompmodel(spw,nui,mode,doClos,onlyClos,doDump)
4675
4707
 
4676
4708
     if mode ==-2:
4677
4709
       if nui == -1: