• Operation343TreeDefinitionLanguage

Operation 343

Operation 343 is an attempt to perform a self calibration of a real WSRT data set using MeqTimba. We have had some success doing that using the Glish interface. Now that there is a Python interface, the "matrix343.g" script is rewritten and modified. It is now called "matrix343.py". This page serves as a kind of log of the issues that I have and their solutions (if available).

Importing MS data into TreeDefinitionLanguage

• During the tree definition (_define_forest) create any parameters that you want to get from a measurement set as Meq.Constant.I used Meq.Parm first, but that did not work.

• Hijack (copy) read_msvis_header.py and modify it to suit your needs

Setting solvable parameters in _tdl_job_something(mqs, parent)

   1 from Timba.TDL import *
   2 
   3 def set_MAB_node_state (mqs, node, fields_record):
   4     """helper function to set the state of a node specified by name or
   5     nodeindex""";
   6     rec = record(state=fields_record);
   7     if isinstance(node,str):
   8         rec.name = node;
   9     elif isinstance(node,int):
  10         rec.nodeindex = node;
  11     else:
  12         raise TypeError,'illegal node argumnent';
  13 # pass command to kernel
  14     mqs.meq('Node.Set.State',rec);
  15     pass
  16 
  17 
  18 
  19 def create_solver_defaults(num_iter=6, solvable=[]):
  20     solver_defaults=record()
  21     solver_defaults.num_iter     = num_iter
  22     solver_defaults.save_funklets= True
  23     solver_defaults.last_update  = True
  24 #See example in TDL/MeqClasses.py
  25     solver_defaults.solvable     = record(command_by_list=(record(name=solvable,
  26                                          state=record(solvable=True)),
  27                                          record(state=record(solvable=False))))
  28     return solver_defaults
  29 
  30 
  31 
  32 def _tdl_job_source_flux_fit_no_calibration(mqs, parent):
  33     msname          = '3C343.MS'
  34     inputrec        = create_inputrec(msname, tile_size=1500)
  35     outputrec       = create_outputrec()
  36 
  37     source_list = create_initial_source_model()
  38 
  39     solvables = []
  40     for source in source_list:
  41         solvables.append('stokes:I:'+source.name)
  42         solvables.append('stokes:Q:'+source.name)
  43         pass
  44 
  45     solver_defaults = create_solver_defaults(solvable=solvables)
  46 
  47     set_MAB_node_state(mqs, 'solver', solver_defaults)
  48 
  49     mqs.init(record(mandate_regular_grid=False),\
                    inputinit=inputrec, outputinit=outputrec)
  50 
  51 
  52     pass

MeqSum / MeqAdd

• If you want to add the results of one or more children, keeping the Cells layout, use MeqAdd.

• If you want to sum/reduce over all points in a VellSet, use MeqSum

Memory consumption and runtime

When solving for I and Q for both sources (12 coefficients total), 1 channel, 1438 timeslots per tile, the runtime per solver iteration is approximately 4 seconds. Total memory consumption is approximately 600 MB

6 channels: 39 seconds per iteration: memory consumption 2600 MB. This implies that every channel costs of the order of 400 MB with the current tiling and cache setup (100). 6.5s /channel/iteration

16 channels, 1438 timeslots, cachepolicy: 0, memory: 2100 MB, 1m50s per iteration. 120MB RAM per channel. 6.9s/channel/iteration

Phase solutions

solving for antenna based global phases (xx and yy): 2.57 seconds/timeslot with 16 channels, two sources, cachepolicy: 0. Memory consumption at start: 290 MB. Memory consumption near end: 1370 MB. Is this normal, or is something leaking 1MB of memory per solution interval?

Interesting effect: when solving for spatially dependent phases, the solutions fail horribly. I have also seen this in the Glish version of the reduction. We need to find out whether this is fundamental or not.

Tiled solutions: tile_size: 150 timeslots. Tiling in time: 1 timeslot, frequency: 1 timeslot. Tiling is specified in thee initrec of a MeqParm, for example: tiling=record(freq=51,time=30). If one does not want any tiling in frequency, omit the freq= parameter.

The total number of solvables per timeslot is 13*2 = 26. This setup results in the following timings:

When the solve matrix contains approximately 260 solvables, the matrix inversion starts to dominate the runtime for this setup.

Initrecords and things ending on _index

!!!! initrec fields ending on _index are 0-based in python and 1-based in Glish !!!!