Line data Source code
1 : /*
2 : * PointingDirectionCache.cc
3 : *
4 : * Created on: Dec 1, 2016
5 : * Author: jjacobs
6 : */
7 :
8 : #include <msvis/MSVis/PointingDirectionCache.h>
9 :
10 : #include <casacore/casa/Exceptions/Error.h>
11 : #include <casacore/casa/BasicSL/String.h>
12 : #include <casacore/measures/Measures/MDirection.h>
13 : #include <casacore/ms/MeasurementSets/MSPointingColumns.h>
14 :
15 : #include <memory>
16 :
17 : using namespace casacore;
18 :
19 : namespace casa {
20 : namespace vi {
21 :
22 0 : PointingColumns::PointingColumns (const MSPointingColumns & pointingColumns)
23 0 : : pointingColumns_p (pointingColumns)
24 0 : {}
25 :
26 : Pointing
27 0 : PointingColumns::getPointingRow (int row, double /*targetTime*/, bool asMeasure) const
28 : {
29 0 : Pointing pointing;
30 :
31 : // See if the row's time is no earlier than 5 minutes before target
32 : // An optimization to avoid unnecessary direction measure conversions
33 : // for rows unlikely to be used (the conversions are fairly expensive)
34 :
35 0 : pointing.time = pointingColumns_p.time () (row);
36 0 : pointing.interval = pointingColumns_p.interval () (row);
37 0 : pointing.antennaId = pointingColumns_p.antennaId() (row);
38 0 : pointing.row = row;
39 0 : pointing.source = this;
40 :
41 0 : if (asMeasure){
42 :
43 : // Pointing is likely to be useful so get all the data.
44 :
45 0 : pointing.direction.reset (new MDirection (pointingColumns_p.directionMeas (row)));
46 : }
47 :
48 0 : return pointing;
49 0 : }
50 :
51 : int
52 0 : PointingColumns::nRows () const
53 : {
54 0 : return pointingColumns_p.nrow();
55 : }
56 :
57 :
58 0 : PointingDirectionCache::PointingDirectionCache (int nAntennas,
59 0 : const PointingSource & pointingSource)
60 0 : : antennaLevelCache_p (MinTimeEntries, MaxTimeEntries, nAntennas, this),
61 0 : antennaEarliestTime_p (nAntennas, -1.0),
62 0 : lastRowRead_p (-1),
63 0 : nFallbacks_p (0),
64 0 : nHits_p (0),
65 0 : pointingEofReached_p (false),
66 0 : pointingSource_p (pointingSource)
67 0 : {}
68 :
69 0 : PointingDirectionCache::~PointingDirectionCache ()
70 : {
71 : // printf ("--> nHits=%d, nFallbacks_p=%d, %%fallback=%f, unused=0%3.2f\n", nHits_p, nFallbacks_p,
72 : // 100.0 * nFallbacks_p / (nHits_p + nFallbacks_p), (nAdded_p - nHits_p) * 1.0f / nAdded_p);
73 0 : }
74 :
75 : void
76 0 : PointingDirectionCache::fillCache (int antenna, double time, bool flushAndRewind) const
77 : {
78 0 : if (flushAndRewind){
79 :
80 : // A time earlier than what was cached is needed: reload the cache from the
81 : // beginning. If done with any frequency on large files it could produce a
82 : // performance problem.
83 :
84 0 : antennaLevelCache_p.flushTimes ();
85 0 : lastRowRead_p = -1;
86 : }
87 :
88 0 : int nRows = pointingSource_p.nRows();
89 :
90 : while (true) {
91 :
92 : // Begin reading rows from the pointing table. Each row will be an entry for
93 : // an antenna at a particular time. The implementation will fail if the
94 : // timestamps for each antenna are not monotonically increasing.
95 :
96 0 : lastRowRead_p ++;
97 :
98 0 : if (lastRowRead_p >= nRows){
99 0 : pointingEofReached_p = true;
100 0 : return; // can't fill anymore
101 : }
102 :
103 : // Pull out the information from the current row of the pointing table and
104 : // add it to the cache.
105 :
106 0 : Pointing pointing = pointingSource_p.getPointingRow (lastRowRead_p, time, false);
107 0 : antennaLevelCache_p.addEntry (pointing);
108 :
109 0 : if (antennaEarliestTime_p [pointing.antennaId] < 0) {
110 0 : antennaEarliestTime_p [pointing.antennaId] = pointing.time;
111 : }
112 :
113 : // See if the current row satisfies the underlying cache read request.
114 : // If so, return.
115 :
116 0 : if (pointing.antennaId == antenna){
117 :
118 0 : if (pd_cache::timeMatch (time, pointing.time, pointing.interval)){
119 0 : return;
120 : }
121 :
122 0 : if (time < pointing.time){
123 : // Target time before first pointing time, so give up.
124 0 : return;
125 : }
126 : }
127 0 : }
128 : }
129 :
130 : std::pair <bool, casacore::MDirection>
131 0 : PointingDirectionCache::getPointingDirection (int antenna, double time, const MDirection & phaseCenter) const
132 : {
133 : using namespace pd_cache;
134 :
135 : // Try to fill the direction request from the cache.
136 :
137 0 : if (noDataForAntenna (antenna, time)){
138 0 : nFallbacks_p ++;
139 0 : return std::make_pair (true, phaseCenter); // antenna not in pointing subtable --> return default.
140 : }
141 :
142 : CacheAccessStatus status;
143 0 : const MDirection * direction = nullptr;
144 0 : MDirection result;
145 :
146 0 : std::tie (status, direction) = antennaLevelCache_p.getPointingDirection (antenna, time);
147 :
148 0 : if (status == CacheAccessStatus::Hit){
149 :
150 0 : result = * direction;
151 0 : nHits_p ++;
152 :
153 0 : } else if (status == CacheAccessStatus::MissInternal){
154 :
155 : // The cache contained values surrounding the requested value but not that
156 : // matched it. Use the phaseCenter as a reasonable fallback.
157 :
158 0 : result = phaseCenter;
159 0 : nFallbacks_p ++;
160 :
161 : } else {
162 :
163 : // The request failed so add more data to the cache.
164 :
165 0 : bool flushAndRewind = status == CacheAccessStatus::MissPrior;
166 :
167 0 : fillCache (antenna, time, flushAndRewind);
168 :
169 : // Try to fill the direction request again.
170 :
171 0 : std::tie (status, direction) = antennaLevelCache_p.getPointingDirection (antenna, time);
172 :
173 0 : if (status != CacheAccessStatus::Hit){
174 :
175 : // Miss after cache fill so use the phaseCenter as a fallback.
176 :
177 0 : result = phaseCenter;
178 0 : nFallbacks_p ++;
179 : } else {
180 0 : result = * direction;
181 0 : nHits_p ++;
182 : }
183 : }
184 :
185 0 : bool northPole = result.toString() == "00:00:00.000000 90.00.00.00000 J2000";
186 0 : if (northPole){
187 0 : printf ("Returning north pole: ant=%d, t=%f\n", antenna, time);
188 : }
189 :
190 0 : return std::make_pair (true, result);
191 0 : }
192 :
193 : const PointingSource *
194 0 : PointingDirectionCache::getPointingSource () const
195 : {
196 0 : return & pointingSource_p;
197 : }
198 :
199 : bool
200 0 : PointingDirectionCache::noDataForAntenna (int antenna, double time) const
201 : {
202 : // If the entire pointing subtable has been scanned once and the antenna
203 : // was not seen then return true. This is intended to prevent excessive
204 : // reading of the pointing subtable when an antenna is missing.
205 :
206 0 : bool noData = false;
207 :
208 0 : if (antennaEarliestTime_p [antenna] >= 0){
209 :
210 : // If there is an earliest time for the antenna, then see if the requested
211 : // time is before it.
212 :
213 0 : noData = time < antennaEarliestTime_p [antenna];
214 :
215 : } else {
216 :
217 : // No data has been seen (yet?) for the current antenna.
218 : // If the entire subtable has been scanned once already, then
219 : // there is no data for this antenna.
220 :
221 0 : noData = pointingEofReached_p;
222 : }
223 :
224 0 : return noData;
225 : }
226 :
227 : namespace pd_cache {
228 :
229 : bool
230 0 : timeMatch (double time, double rowsTime, double rowsInterval)
231 : {
232 : // The time matches if it falls within an interval around the pointing
233 : // table time.
234 :
235 0 : bool match = time >= rowsTime - rowsInterval &&
236 0 : time <= rowsTime + rowsInterval;
237 :
238 0 : return match;
239 : }
240 :
241 0 : AntennaLevelCache::AntennaLevelCache (int minTimes, int maxTimes, int nAntennas,
242 0 : const PointingDirectionCache * pdCache)
243 0 : : pointingDirectionCache_p (pdCache),
244 0 : timeLevelCache_p (nAntennas, TimeLevelCache (minTimes, maxTimes, this))
245 0 : {}
246 :
247 :
248 : void
249 0 : AntennaLevelCache::addEntry (Pointing & pointing)
250 : {
251 0 : timeLevelCache_p [pointing.antennaId].addEntry (pointing);
252 0 : }
253 :
254 : void
255 0 : AntennaLevelCache::flushTimes ()
256 : {
257 0 : for (auto & timeCache : timeLevelCache_p){
258 :
259 : // Flush the time cache for each antenna
260 :
261 0 : timeCache.flush ();
262 : }
263 0 : }
264 :
265 : std::pair<CacheAccessStatus,const casacore::MDirection *>
266 0 : AntennaLevelCache::getPointingDirection (int antenna, double time)
267 : {
268 0 : return timeLevelCache_p [antenna].getPointingDirection(time);
269 : }
270 :
271 : const PointingSource *
272 0 : AntennaLevelCache::getPointingSource () const
273 : {
274 0 : return pointingDirectionCache_p->getPointingSource ();
275 : }
276 :
277 0 : TimeLevelEntry::TimeLevelEntry (const Pointing & pointing,
278 0 : const TimeLevelCache * tlCache)
279 0 : : direction_p (std::move (pointing.direction)),
280 0 : row_p (pointing.row),
281 0 : timeCenter_p (pointing.time),
282 0 : timeLevelCache_p (tlCache),
283 0 : interval_p (pointing.interval)
284 : {
285 0 : bool northPole = direction_p && direction_p->toString() == "00:00:00.000000 90.00.00.00000 J2000";
286 0 : if (northPole){
287 0 : printf ("Returning north pole: ant=%d, t=%f\n", pointing.antennaId, pointing.time);
288 : }
289 0 : }
290 :
291 0 : TimeLevelEntry::TimeLevelEntry (const TimeLevelEntry & other)
292 0 : : direction_p (std::move (other.direction_p)),
293 0 : row_p (other.row_p),
294 0 : timeCenter_p (other.timeCenter_p),
295 0 : timeLevelCache_p (other.timeLevelCache_p),
296 0 : interval_p (other.interval_p)
297 0 : {}
298 :
299 0 : TimeLevelEntry::~TimeLevelEntry ()
300 : {
301 0 : }
302 :
303 : TimeLevelEntry &
304 0 : TimeLevelEntry::operator= (const TimeLevelEntry & other)
305 : {
306 0 : if (& other != this){
307 :
308 0 : direction_p = std::move (other.direction_p);
309 0 : other.direction_p = nullptr;
310 :
311 0 : row_p = other.row_p;
312 0 : timeCenter_p = other.timeCenter_p;
313 0 : interval_p = other.interval_p;
314 : }
315 :
316 0 : return * this;
317 : }
318 :
319 : const casacore::MDirection *
320 0 : TimeLevelEntry::getDirection () const{
321 :
322 : // If the direction measure is actually cached, then simply return it.
323 :
324 0 : if (! direction_p){
325 :
326 : // Refetch the row and actually get the direction this time. Put it in
327 : // the cache and return a pointer to it.
328 :
329 0 : Pointing p = timeLevelCache_p->getPointingSource()->getPointingRow (row_p, timeCenter_p, true);
330 0 : direction_p.reset (new MDirection (* p.direction.get()));
331 :
332 0 : }
333 :
334 0 : return direction_p.get();
335 : }
336 :
337 : double
338 0 : TimeLevelEntry::getInterval () const{
339 0 : return interval_p;
340 : }
341 :
342 : double
343 0 : TimeLevelEntry::getTime () const
344 : {
345 0 : return timeCenter_p;
346 : }
347 :
348 :
349 : bool
350 0 : operator== (double t, const TimeLevelEntry & tle)
351 : {
352 0 : return timeMatch (t, tle.getTime(), tle.getInterval());
353 : }
354 :
355 : bool
356 0 : operator== (const TimeLevelEntry & tle, double t)
357 : {
358 0 : return t == tle;
359 : }
360 :
361 : bool
362 0 : operator< (double t, const TimeLevelEntry & tle)
363 : {
364 0 : return t < tle.getTime() - tle.getInterval();
365 : }
366 :
367 : bool
368 0 : operator< (const TimeLevelEntry & tle, double t)
369 : {
370 0 : return t > tle.getTime() + tle.getInterval();
371 : }
372 :
373 :
374 0 : TimeLevelCache::TimeLevelCache (int minTimes, int maxTimes, const AntennaLevelCache * alCache)
375 0 : : antennaLevelCache_p (alCache),
376 0 : cache_p (),
377 0 : maxTimes_p (maxTimes),
378 0 : minTimes_p (minTimes)
379 0 : {}
380 :
381 :
382 : void
383 0 : TimeLevelCache::addEntry (Pointing & pointing)
384 : {
385 0 : if (cache_p.size() + 1 > (uint) maxTimes_p){
386 :
387 : // Compact the cache so that it contains minTimes_p elements.
388 :
389 0 : int nToErase = cache_p.size() - minTimes_p;
390 0 : cache_p.erase (cache_p.begin(), cache_p.begin() + nToErase);
391 : }
392 :
393 : // Add in the new entry on the end.
394 :
395 0 : if (! cache_p.empty() && pointing.time < cache_p.back().getTime()){
396 :
397 : // The pointing table is not increasing monotonically for this antenna.
398 : // As a kluge, flush the cache and continue. This could behave badly
399 : // in some situations.
400 :
401 0 : cache_p.clear();
402 : }
403 :
404 0 : cache_p.push_back (TimeLevelEntry (pointing, this));
405 0 : }
406 :
407 : void
408 0 : TimeLevelCache::flush ()
409 : {
410 : // Get rid of everything.
411 :
412 0 : cache_p.clear();
413 0 : }
414 :
415 : std::pair<CacheAccessStatus, const casacore::MDirection *>
416 0 : TimeLevelCache::getPointingDirection (double time)
417 : {
418 0 : if (cache_p.empty ()){
419 0 : return std::make_pair (CacheAccessStatus::MissPost, nullptr);
420 : }
421 :
422 : // Find the first element that is >= the requested time.
423 :
424 0 : Cache::iterator lowerBound = std::lower_bound (cache_p.begin(), cache_p.end(), time);
425 :
426 0 : if (lowerBound == cache_p.end()){
427 :
428 : // Handle the case where there is no lower bound
429 :
430 0 : TimeLevelEntry & tle = cache_p.back();
431 0 : double upperTime = tle.getTime() + tle.getInterval();
432 :
433 0 : if (time <= upperTime){
434 :
435 : // Last element in cache's interval contains the target
436 : // time: we found it
437 :
438 0 : return std::make_pair (CacheAccessStatus::Hit, tle.getDirection());
439 : }
440 :
441 : // Off the top end of existing data in the cache
442 :
443 0 : return std::make_pair (CacheAccessStatus::MissPost, nullptr);
444 : }
445 :
446 0 : if (time == * lowerBound){
447 :
448 : // We've found it
449 :
450 0 : return std::make_pair (CacheAccessStatus::Hit, lowerBound->getDirection());
451 :
452 : }
453 :
454 0 : if (lowerBound == cache_p.begin()){
455 :
456 : // Nothing below it so it's a miss
457 :
458 0 : return std::make_pair (CacheAccessStatus::MissPrior, nullptr);
459 : }
460 :
461 0 : if (time == * (lowerBound - 1)){
462 :
463 0 : return std::make_pair (CacheAccessStatus::Hit,
464 0 : (lowerBound - 1)->getDirection());
465 : }
466 :
467 : // Cache miss: must be between two cache entries (this could be handled
468 : // by interpolating between the two entries, but might indicate a flaw in
469 : // the algorithm or the data; for now just let the caller sort it out).
470 :
471 0 : return std::make_pair (CacheAccessStatus::MissInternal, nullptr);
472 : }
473 :
474 : const PointingSource *
475 0 : TimeLevelCache::getPointingSource () const
476 : {
477 0 : return antennaLevelCache_p->getPointingSource ();
478 : }
479 :
480 :
481 : } // end namespace pd_cache
482 :
483 : } // end namespace vi
484 : } // end namepsace casa
485 :
|
