InitialWorkingExample

  # stations list  STATIONS = range(1,15);  # 3 sources  SOURCES = ('a','b','c');    #--- these are generated automatically from the station list  # list of ifrs as station pairs: each entry is two indices, (s1,s2)  IFRS   = [ (s1,s2) for s1 in STATIONS for s2 in STATIONS if s1<s2 ];  # list of ifr strings of the form "s1-s2"  QIFRS  = [ '-'.join(map(str,ifr)) for ifr in IFRS ];  # combined list: each entry is ("s1-s2",s1,s2)  QQIFRS = [ ('-'.join(map(str,ifr)),) + ifr for ifr in IFRS ];    # namespace of all node classes  Meq = Namespace('Meq');     # global node scope & repository  nsglob = NodeScope();       # init some node groups  ROOT = NodeGroup('root');  SOLVERS = NodeGroup('solvers');    #------- create nodes for instrumental model  # some handy aliases  ZERO = nsglob.zero() << Meq.Constant(value=0);  UNITY = nsglob.unity() << Meq.Constant(value=1);    PHASE_CENTER_RA  = nsglob.ra0() << Meq.Parm();  PHASE_CENTER_DEC = nsglob.dec0() << Meq.Parm();    STATION_UWV = {};  STATION_POS = {};    ARRAY_POS = nsglob.xyz0() << Meq.Composer(    nsglob.x0() << Meq.Parm(),    nsglob.y0() << Meq.Parm(),    nsglob.z0() << Meq.Parm() );      # create station-related nodes and branches  for s in STATIONS:    STATION_POS[s] = nsglob.xyz(s) << Meq.Composer(      nsglob.x(s) << Meq.Parm(),      nsglob.y(s) << Meq.Parm(),      nsglob.z(s) << Meq.Parm() );    STATION_UWV[s] = nsglob.stuwv(s) << Meq.UWV(children={      'xyz': STATION_POS[s],      'xyz0': ARRAY_POS,      'ra': PHASE_CENTER_RA,      'dec': PHASE_CENTER_DEC    });    # create per-source station gains    for (q,src) in enumerate(SOURCES):      nsglob.G(s,q=q) << Meq.Composer(        nsglob.Gxx(s,q=q) << Meq.Polar(Meq.Parm(),Meq.Parm()), ZERO,        ZERO, nsglob.Gyy(s,q=q) << Meq.Polar(Meq.Parm(),Meq.Parm()),      dims=[2,2]);    # alternative: single gain with no direction dependence    # nsglob.G(s) << Meq.Composer(    #    nsglob.Gxx(s) << Meq.Polar(Meq.Parm(),Meq.Parm()), ZERO,    #    ZERO, nsglob.Gyy(s) << Meq.Polar(Meq.Parm(),Meq.Parm()),    #  dims=[2,2]);    # this function returns a per-station gain node, given a set of qualifiers  def STATION_GAIN (s=s,q=q,**qual):      # **qual swallows any remaining qualifiers    return nsglob.G(s,q=q);    # note alternative for no direction dependence:    # def STATION_GAIN (s=s,**qual):      #   return nsglob.G(s);      #------- end of instrumental model        #------- create model for unpolarized point source  # References instrumental model: STATION_GAIN(s,**qual), STATION_UWV[s].  # Returns unqualified predict node, should be qualified with ifr string.  def makeUnpolarizedPointSource (ns,**qual):    ns.lmn() << Meq.LMN(children={        'ra':   ns.ra() << Meq.Parm(),        'dec':  ns.dec() << Meq.Parm(),        'ra0':  PHASE_CENTER_RA,        'dec0': PHASE_CENTER_DEC    });    ns.stokes_i() << Meq.Parm();    # create per-station term subtrees    for s in STATIONS:      ns.sdft(s) << Meq.MatrixMultiply(        STATION_GAIN(s,**qual),        Meq.StatPSDFT(ns.lmn(),STATION_UWV[s])      );    # create per-baseline predicters    for (q,s1,s2) in QQIFRS:      ns.predict(q) << Meq.Multiply(          ns.stokes_i(),          ns.dft(q) << Meq.DFT(ns.sdft(s1),ns.sdft(s2)),      );    return ns.predict;      #------- create peeling unit  # inputs: an unqualified input node, will be qualified with ifr string.  # predicters: list of unqualified predict nodes, will be qualified with ifr string.  # Returns unqualified output node, should be qualified with ifr string.  def peelUnit (inputs,predicters,ns):    for q in QIFRS:      # create condeq branch      ns.condeq(q) << Meq.Condeq(        ns.measured(q) << Meq.PhaseShift(children=inputs(q)),        ns.predicted(q) << Meq.Add(*[prd(q) for prd in predicters])      );      # create subtract branch      ns.subtract(q) << Meq.Subtract(ns.measured(q),ns.predicted(q));    # creates solver and sequencers    ns.solver() << Meq.Solver(*[ns.condeq(q) for q in QIFRS]);    for q in QIFRS:      ns.reqseq(q) << Meq.ReqSeq(ns.solver(),ns.subtract(q));    # returns root nodes of unit    return ns.reqseq;    # create source predictors, each in its own subscope  predicter = {};  for (q,src) in enumerate(SOURCES):    predicter[q] = makeUnpolarizedPointSource(nsglob.subscope('predict',src),q=q);    # create spigots  for q in QIFRS:    nsglob.spigot(q) << Meq.Spigot();      # chain peel units, by connecting outputs to inputs. First input  # is spigot.  inputs = nsglob.spigot;  for (q,src) in enumerate(SOURCES):    ns_pu = nsglob.subscope('peelunit',q);    inputs = peelUnit(inputs,predicter.values(),ns=ns_pu);    SOLVERS << ns_pu.solver();      # create sinks, connect them to output of last peel unit  for q in QIFRS:    ROOT << nsglob.sink(q) << Meq.Sink(inputs(q));      # create data collectors (this simply shows off the use of arbitrary node  # groupings)  ROOT << nsglob.solver_collect() << Meq.DataCollect(*SOLVERS.values());      # deliberately create an orphan branch. This checkes orphan collection.  # this whole branch should go away, and so should the UNITY node, which  # is not used anywhere  nsglob.orphan() << Meq.Add(Meq.Const(value=0),UNITY,Meq.Add(UNITY,ZERO));  # resolve all nodes  nsglob.resolve(ROOT);```