I've seen a few instances of people suffering in silence (or worse, proclaiming that MeqTrees are too slow/too memory-hungry for their needs) when their trees are suboptimal by orders of magnitude... It's a (mostly social) problem that worries me and I have no idea how to solve it yet.
What this means is - don't suffer in silence! Ask a question! A MeqTrees guru may be able to help you.
For example, Oleg pointed out that the following script could be sped up (tests showed by a factor of about 20) by changing the order of operations (which allows us to eliminate the time-expensive Compounder node). Trace out the sequence of operations when the `pre_rotate' flag is set to True (the script will run very quickly) vs what is done when `pre-rotate' is set to False.
1 # a script to generate a sequence of CLAR beam images as a
2 # function of time.
3 # The original script is available as Timba/WH/contrib/AGW/MG_AGW_clar_beam.py
4 #
5
6 from Timba.TDL import *
7 from Timba.Meq import meq
8 from Timba.Meq import meqds
9
10 from make_multi_dim_request import *
11
12 import Meow.Bookmarks
13
14 # setup a few bookmarks
15 Settings.forest_state = record(bookmarks=[
16 Meow.Bookmarks.PlotPage("CLAR Beam",["beam_rot"],["AzEl"])
17 ]);
18
19 # Timba.TDL.Settings.forest_state is a standard TDL name.
20 # This is a record passed to Set.Forest.State.
21 Settings.forest_state.cache_policy = 100;
22
23 # get position of field centre RA DEC
24 TDLCompileMenu('Field Centre RA and DEC',
25 TDLOption('fc_ra','RA of field centre (radians)',[0,1,2,3],more=float),
26 TDLOption('fc_dec','DEC of field centre (radians)',[0,1,2,3],more=float),
27 );
28
29 # get telescope location for use in simulation
30 TDLCompileOption('telescope','Telescope Site',['DRAO','VLA'])
31
32 # do fast or slow calculation
33 TDLCompileOption('pre_rotate','Pre-Rotate L and M frame (fast method)',[True, False])
34
35
36
37 def _define_forest (ns):
38 """define_forest() is a standard TDL name. When a forest script is
39 loaded by, e.g., the browser, this method is automatically called to
40 define the forest. The 'ns' argument is a NodeScope object in which
41 the forest is to be defined, usually this is simply the global scope.
42 """;
43
44 # we build up the mathematical expression of a CLAR power
45 # pattern using the formulae
46 # log16 = (-1.0) * log(16.0)
47 # L,M give direction cosines wrt field centre in
48 # AzEl coordinate frame
49 # L_gain = (L * L) / (width * width)
50 # sin_factor = sqr(sin(field centre elevation))
51 # M_gain = (sin_factor * M * M ) / (width * width)
52 # power pattern = exp(log16 * (L_gain + M_gain))
53
54 # constant used in beam calculation
55 ns.ln_16 << Meq.Constant(-2.7725887)
56
57 # CLAR HPBW of 3 arcmin at zenith = 0.00087266 radians
58 ns.HPBW << Meq.Constant(0.00087266)
59
60 # beam is centred in L,M grid
61 ns.centre << Meq.Constant(0.0)
62
63 # create a MeqComposer containing ra dec position
64 ns.RADec <<Meq.Composer(fc_ra, fc_dec)
65
66 # we create an AzEl node with an Observatory name
67 ns.AzEl << Meq.AzEl(radec=ns.RADec, observatory=telescope)
68
69 # create a Parallactic angle node
70 pa = ns.ParAngle << Meq.ParAngle(radec=ns.RADec, observatory = telescope)
71 # rotation matrix to go from AzEl to Ra,Dec L,M
72 ns.P << Meq.Matrix22(Meq.Cos(pa),-Meq.Sin(pa),Meq.Sin(pa),Meq.Cos(pa))
73
74 # get the elevation
75 ns.El << Meq.Selector(ns.AzEl, index=1)
76
77 # get sine of elevation - used to get CLAR beam broadening
78 ns.sine_el << Meq.Sin(ns.El)
79
80 # square this sine value
81 ns.sine_el_sq << Meq.Sqr(ns.sine_el)
82
83 # create L and M axis nodes
84 laxis = ns.laxis << Meq.Grid(axis=2);
85 maxis = ns.maxis << Meq.Grid(axis=3);
86 # attempt to de-rotate the beam in AzEl coordinates to sky coordinates
87 if pre_rotate:
88 ns.lm_pre_rot << Meq.Composer(laxis,maxis) # returns an lm 2-vector
89
90 ns.rot_lm << Meq.MatrixMultiply(ns.P,ns.lm_pre_rot); # rotated lm
91 ns.l_rot << Meq.Selector(ns.rot_lm,index=0)
92 ns.m_rot << Meq.Selector(ns.rot_lm,index=1)
93
94 ns.l_sq << Meq.Sqr(ns.l_rot - ns.centre)
95 ns.m_sq << Meq.Sqr(ns.m_rot - ns.centre)
96 ns.m_sq_sin << Meq.Multiply(ns.m_sq, ns.sine_el_sq)
97
98 # Add l and m gains
99 ns.l_and_m_sq << Meq.Add(ns.l_sq, ns.m_sq_sin)
100 ns.beam_rot << Meq.Exp((ns.l_and_m_sq * ns.ln_16) / Meq.Sqr(ns.HPBW))
101
102 else:
103 ns.l_sq << Meq.Sqr(laxis - ns.centre)
104 ns.m_sq << Meq.Sqr(maxis - ns.centre)
105 ns.m_sq_sin << Meq.Multiply(ns.m_sq, ns.sine_el_sq)
106
107 # Add l and m gains
108 ns.l_and_m_sq << Meq.Add(ns.l_sq, ns.m_sq_sin)
109 ns.exp_gain << Meq.Exp((ns.l_and_m_sq * ns.ln_16) / Meq.Sqr(ns.HPBW))
110
111 # attempt to de-rotate the beam in AzEl coordinates to sky coordinates
112 ns.lm_pre_rot << Meq.Composer(laxis,maxis) # returns an lm 2-vector
113
114 ns.rot_lm << Meq.MatrixMultiply(ns.P,ns.lm_pre_rot); # rotated lm
115 ns.l_rot << Meq.Selector(ns.rot_lm,index=0)
116 ns.m_rot << Meq.Selector(ns.rot_lm,index=1)
117 ns.lm_rot << Meq.Composer(Meq.Grid(axis=0),Meq.Grid(axis=1),ns.l_rot,ns.m_rot)
118
119 ns.resampler << Meq.Resampler(ns.exp_gain,dep_mask = 0xff)
120 ns.beam_rot << Meq.Compounder(children=[ns.lm_rot,ns.resampler],common_axes=[hiid('l'),hiid('m')])
121
122 ##################################################
123 def _test_forest (mqs,parent,wait=False):
124 """test_forest() is a standard TDL name. When a forest script is
125 loaded by, e.g., the browser, and the "test" option is set to true,
126 this method is automatically called after define_forest() to run a
127 test on the forest. The 'mqs' argument is a meqserver proxy object.
128 """;
129
130 # any large time range will do: we observe the changes in the beam
131 # pattern in timesteps of 2400s, or 40 min
132 delta_t = 2400.0
133 # Approx start of 2001
134 t0 = 4485011731.14 - delta_t
135 t1 = t0 + delta_t
136
137 # here, any frequency range will do
138 f0 = 800.0
139 f1 = 1300.0
140
141 m_range = [-0.0025,0.0025];
142 l_range = [-0.0025,0.0025];
143 lm_num = 31;
144 counter = 0
145 for i in range(12):
146 t0 = t0 + delta_t
147 t1 = t0 + delta_t
148 request = make_multi_dim_request(counter=counter, dom_range = [[f0,f1],[t0,t1],l_range,m_range], nr_cells = [1,1,lm_num,lm_num])
149 counter = counter + 1
150 # execute request
151 mqs.meq('Node.Execute',record(name='beam_rot',request=request),wait=False);
152 ##################################################
153
154 if __name__=='__main__':
155 ns=NodeScope()
156 _define_forest(ns)
157 ns.Resolve()
158 print "Added %d nodes" % len(ns.AllNodes())