Commits

Takeshi Nakazato authored 0fa342f675a 
CAS-13842 use MSColumn object which is held by VI/VB2
cffd9dc77f7
master

casatools/src/code/msvis/MSVis/VisBufferUtil.cc

Modified
38 38 #include <casa/Arrays/Cube.h>
39 39 #include <casacore/casa/Utilities/Sort.h>
40 40 #include <casa/OS/Timer.h>
41 41 #include <casa/OS/Path.h>
42 42 #include <measures/Measures/UVWMachine.h>
43 43 #include <measures/Measures/MeasTable.h>
44 44 #include <ms/MSSel/MSSelectionTools.h>
45 45 #include <ms/MeasurementSets/MSPointingColumns.h>
46 46 #include <msvis/MSVis/VisBufferUtil.h>
47 47 #include <msvis/MSVis/StokesVector.h>
48 +#include <msvis/MSVis/ViImplementation2.h>
48 49 #include <msvis/MSVis/VisibilityIterator.h>
49 50 #include <msvis/MSVis/VisibilityIterator2.h>
50 51 #include <msvis/MSVis/VisBuffer.h>
51 52 #include <ms/MeasurementSets/MSColumns.h>
52 53 #include <casa/iostream.h>
53 54 #include <fstream>
54 55 #include <iomanip>
55 56 #ifdef _OPENMP
56 57 #include <omp.h>
57 58 #endif
58 59 using namespace std;
59 60 ​
60 61 using namespace casacore;
61 62 namespace casa { //# NAMESPACE CASA - BEGIN
62 63 ​
63 -// <summary>
64 +// <summary>
64 65 // </summary>
65 66 ​
66 67 // <reviewed reviewer="" date="" tests="tMEGI" demos="">
67 68 ​
68 69 // <prerequisite>
69 70 // </prerequisite>
70 71 //
71 72 // <etymology>
72 73 // </etymology>
73 74 //
74 -// <synopsis>
75 -// </synopsis>
75 +// <synopsis>
76 +// </synopsis>
76 77 //
77 78 // <example>
78 79 // <srcblock>
79 80 // </srcblock>
80 81 // </example>
81 82 //
82 83 // <motivation>
83 84 // </motivation>
84 85 //
85 86 // <todo asof="">
86 87 // </todo>
87 88 ​
88 89 ​
89 90 VisBufferUtil::VisBufferUtil(): oldMSId_p(-1), oldPCMSId_p(-1), timeAntIndex_p(0), cachedPointingDir_p(0), cachedPhaseCenter_p(0){};
90 91 ​
91 92 ​
92 93 // Construct from a VisBuffer (sets frame info)
93 94 VisBufferUtil::VisBufferUtil(const VisBuffer& vb): oldMSId_p(-1), oldPCMSId_p(-1), timeAntIndex_p(0), cachedPointingDir_p(0), cachedPhaseCenter_p(0) {
94 95 ​
95 - // The nominal epoch
96 + // The nominal epoch
96 97 MEpoch ep=vb.msColumns().timeMeas()(0);
97 98 ​
98 99 // The nominal position
99 100 String observatory;
100 101 MPosition pos;
101 102 if (vb.msColumns().observation().nrow() > 0) {
102 103 observatory = vb.msColumns().observation().telescopeName()
103 104 (vb.msColumns().observationId()(0));
104 105 }
105 - if (observatory.length() == 0 ||
106 + if (observatory.length() == 0 ||
106 107 !MeasTable::Observatory(pos,observatory)) {
107 108 // unknown observatory, use first antenna
108 109 pos=vb.msColumns().antenna().positionMeas()(0);
109 110 }
110 -
111 +​
111 112 // The nominal direction
112 113 MDirection dir=vb.phaseCenter();
113 114 ​
114 115 // The nominal MeasFrame
115 116 mframe_=MeasFrame(ep, pos, dir);
116 117 ​
117 118 }
118 119 ​
119 120 // Construct from a VisBuffer (sets frame info)
120 121 VisBufferUtil::VisBufferUtil(const vi::VisBuffer2& vb): oldMSId_p(-1), oldPCMSId_p(-1), timeAntIndex_p(0), cachedPointingDir_p(0), cachedPhaseCenter_p(0) {
183 184 false, false);
184 185 retval = !uvwMachine.isNOP();
185 186 dphase.resize(vb.nRows());
186 187 dphase.set(0.0);
187 188 if(uvw.nelements() ==0)
188 189 uvw=vb.uvw();
189 190 for (rownr_t row=0; row< vb.nRows(); ++row){
190 191 Vector<Double> eluvw(uvw.column(row));
191 192 uvwMachine.convertUVW(dphase(row), eluvw);
192 193 }
193 -
194 +​
194 195 return retval;
195 196 }
196 197 ​
197 198 // Set the visibility buffer for a PSF
198 199 void VisBufferUtil::makePSFVisBuffer(VisBuffer& vb) {
199 200 CStokesVector coh(Complex(1.0), Complex(0.0), Complex(0.0), Complex(1.0));
200 201 vb.correctedVisibility()=coh;
201 202 }
202 203 ​
203 204 ​
204 -Bool VisBufferUtil::interpolateFrequency(Cube<Complex>& data,
205 - Cube<Bool>& flag,
205 +Bool VisBufferUtil::interpolateFrequency(Cube<Complex>& data,
206 + Cube<Bool>& flag,
206 207 const VisBuffer& vb,
207 - const Vector<Float>& outFreqGrid,
208 - const MS::PredefinedColumns whichCol,
208 + const Vector<Float>& outFreqGrid,
209 + const MS::PredefinedColumns whichCol,
209 210 const MFrequency::Types freqFrame,
210 211 const InterpolateArray1D<Float,Complex>::InterpolationMethod interpMethod){
211 212 ​
212 213 Cube<Complex> origdata;
213 214 // Convert the visibility frequency to the frame requested
214 215 Vector<Double> visFreqD;
215 216 convertFrequency(visFreqD, vb, freqFrame);
216 217 //convert it to Float
217 218 Vector<Float> visFreq(visFreqD.nelements());
218 219 convertArray(visFreq, visFreqD);
219 -
220 +​
220 221 //Assign which column is to be regridded to origdata
221 222 if(whichCol==MS::MODEL_DATA){
222 223 origdata.reference(vb.modelVisCube());
223 224 }
224 225 else if(whichCol==MS::CORRECTED_DATA){
225 226 origdata.reference(vb.correctedVisCube());
226 227 }
227 228 else if(whichCol==MS::DATA){
228 229 origdata.reference(vb.visCube());
229 230 }
239 240 //interpolate the data and the flag to the output frequency grid
240 241 InterpolateArray1D<Float,Complex>::
241 242 interpolate(data,flag, outFreqGrid,visFreq,flipdata,flipflag,interpMethod);
242 243 flipdata.resize();
243 244 //reflip the data and flag to be in the same order as in Visbuffer output
244 245 swapyz(flipdata,data);
245 246 data.resize();
246 247 data.reference(flipdata);
247 248 flipflag.resize();
248 249 swapyz(flipflag,flag);
249 - flag.resize();
250 + flag.resize();
250 251 flag.reference(flipflag);
251 252 ​
252 253 return true;
253 254 ​
254 255 }
255 256 void VisBufferUtil::getFreqRange(Double& freqMin, Double& freqMax, vi::VisibilityIterator2& vi, MFrequency::Types freqFrame){
256 257 vi.originChunks();
257 258 vi.origin();
258 259 ​
259 260 Double freqEnd=0.0;
260 261 Double freqStart=C::dbl_max;
261 262 vi::VisBuffer2* vb=vi.getVisBuffer();
262 263 for (vi.originChunks(); vi.moreChunks();vi.nextChunk())
263 264 {
264 265 for (vi.origin(); vi.more();vi.next())
265 266 {
266 267 Double localmax, localmin;
267 - IPosition localmaxpos(1,0);
268 + IPosition localmaxpos(1,0);
268 269 IPosition localminpos(1,0);
269 270 Vector<Double> freqs=vb->getFrequencies(0, freqFrame);
270 271 if(freqs.nelements() ==0){
271 272 throw(AipsError("Frequency selection error" ));
272 273 }
273 274 minMax(localmin,localmax,localminpos, localmaxpos, freqs);
274 275 //localmax=max(freqs);
275 276 //localmin=min(freqs);
276 277 //freqEnd=max(freqEnd, localmax);
277 278 //freqStart=min(freqStart, localmin);
278 279 ​
279 - Int nfreq = freqs.nelements();
280 + Int nfreq = freqs.nelements();
280 281 Vector<Int> curspws = vb->spectralWindows();
281 282 // as the vb row 0 is used for getFrequencies, the same row 0 is used here
282 - Vector<Double> chanWidths = vi.subtableColumns().spectralWindow().chanWidth()(curspws[0]);
283 - // freqs are in channel center freq so add the half the width to the values to return the edge frequencies
283 + Vector<Double> chanWidths = vi.subtableColumns().spectralWindow().chanWidth()(curspws[0]);
284 + // freqs are in channel center freq so add the half the width to the values to return the edge frequencies
284 285 if (nfreq==1) {
285 286 freqEnd=max(freqEnd, localmax+fabs(chanWidths[0]/2.0));
286 287 freqStart=min(freqStart, localmin-fabs(chanWidths[0]/2.0));
287 288 }
288 289 else {
289 290 freqEnd=max(freqEnd, localmax+fabs(chanWidths[localmaxpos[0]]/2.0));
290 291 freqStart=min(freqStart, localmin-fabs(chanWidths[localminpos[0]]/2.0));
291 292 }
292 -
293 +​
293 294 }
294 295 }
295 296 freqMin=freqStart;
296 297 freqMax=freqEnd;
297 298 }
298 299 ​
299 300 Bool VisBufferUtil::getFreqRangeFromRange(casacore::Double& outfreqMin, casacore::Double& outfreqMax, const casacore::MFrequency::Types inFreqFrame, const casacore::Double infreqMin, const casacore::Double infreqMax, vi::VisibilityIterator2& vi, casacore::MFrequency::Types outFreqFrame){
300 301 ​
301 302 ​
302 303 if(inFreqFrame==outFreqFrame){
303 304 outfreqMin=infreqMin;
304 305 outfreqMax=infreqMax;
305 306 return True;
306 307 }
307 308 ​
308 -
309 +​
309 310 vi.originChunks();
310 311 vi.origin();
311 312 try{
312 313 outfreqMin=C::dbl_max;
313 314 outfreqMax=0;
314 315 vi::VisBuffer2* vb=vi.getVisBuffer();
315 316 MSColumns msc(vi.ms());
316 317 // The nominal epoch
317 318 MEpoch ep=msc.timeMeas()(0);
318 -
319 +​
319 320 // The nominal position
320 321 String observatory;
321 322 MPosition pos;
322 323 if (msc.observation().nrow() > 0) {
323 324 observatory = msc.observation().telescopeName()
324 325 (msc.observationId()(0));
325 326 }
326 327 if (observatory.length() == 0 ||
327 328 !MeasTable::Observatory(pos,observatory)) {
328 329 // unknown observatory, use first antenna
329 330 pos=msc.antenna().positionMeas()(0);
330 331 }
331 -
332 +​
332 333 // The nominal direction
333 334 MDirection dir=vb->phaseCenter();
334 335 MeasFrame mFrame(ep, pos, dir);
335 336 // The conversion engine:
336 - MFrequency::Convert toNewFrame(inFreqFrame,
337 + MFrequency::Convert toNewFrame(inFreqFrame,
337 338 MFrequency::Ref(outFreqFrame, mFrame));
338 -
339 +​
339 340 for (vi.originChunks(); vi.moreChunks();vi.nextChunk())
340 341 {
341 342 for (vi.origin(); vi.more();vi.next()){
342 343 //assuming time is fixed in visbuffer
343 344 mFrame.resetEpoch(vb->time()(0)/86400.0);
344 -
345 +​
345 346 // Reset the direction (ASSUMES phaseCenter is constant in the VisBuffer)
346 347 mFrame.resetDirection(vb->phaseCenter());
347 348 Double temp=toNewFrame(infreqMin).getValue().getValue();
348 349 if(temp < outfreqMin)
349 350 outfreqMin = temp;
350 -
351 +​
351 352 temp=toNewFrame(infreqMax).getValue().getValue();
352 353 if(temp > outfreqMax)
353 - outfreqMax = temp;
354 + outfreqMax = temp;
354 355 }
355 356 }
356 357 }
357 358 catch(...){
358 359 //Could not do a conversion
359 360 return False;
360 -
361 +​
361 362 }
362 363 //cerr << "min " << outfreqMin << " max " << outfreqMax << endl;
363 364 return True;
364 -
365 +​
365 366 }
366 367 ​
367 -void VisBufferUtil::convertFrequency(Vector<Double>& outFreq,
368 - const VisBuffer& vb,
368 +void VisBufferUtil::convertFrequency(Vector<Double>& outFreq,
369 + const VisBuffer& vb,
369 370 const MFrequency::Types freqFrame){
370 371 Int spw=vb.spectralWindow();
371 372 MFrequency::Types obsMFreqType=(MFrequency::Types)(vb.msColumns().spectralWindow().measFreqRef()(spw));
372 373 ​
373 374 // The input frequencies (a reference)
374 375 Vector<Double> inFreq(vb.frequency());
375 376 ​
376 377 // The output frequencies
377 378 outFreq.resize(inFreq.nelements());
378 379 ​
392 393 ​
393 394 // Reset the timestamp (ASSUMES TIME is constant in the VisBuffer)
394 395 mframe_.resetEpoch(vb.time()(0)/86400.0);
395 396 ​
396 397 // Reset the direction (ASSUMES phaseCenter is constant in the VisBuffer)
397 398 mframe_.resetDirection(vb.msColumns().field().phaseDirMeasCol()(vb.fieldId())(IPosition(1,0)));
398 399 ​
399 400 // cout << "Frame = " << mframe_ << endl;
400 401 ​
401 402 // The conversion engine:
402 - MFrequency::Convert toNewFrame(obsMFreqType,
403 + MFrequency::Convert toNewFrame(obsMFreqType,
403 404 MFrequency::Ref(newMFreqType, mframe_));
404 405 ​
405 406 // Do the conversion
406 407 for (uInt k=0; k< inFreq.nelements(); ++k)
407 408 outFreq(k)=toNewFrame(inFreq(k)).getValue().getValue();
408 -
409 +​
409 410 }
410 411 else{
411 412 // The requested frame is the same as the observed frame
412 413 outFreq=inFreq;
413 414 }
414 415 ​
415 416 }
416 417 ​
417 -void VisBufferUtil::convertFrequency(Vector<Double>& outFreq,
418 - const vi::VisBuffer2& vb,
418 +void VisBufferUtil::convertFrequency(Vector<Double>& outFreq,
419 + const vi::VisBuffer2& vb,
419 420 const MFrequency::Types freqFrame){
420 421 Int spw=vb.spectralWindows()(0);
421 422 MFrequency::Types obsMFreqType=(MFrequency::Types)(MSColumns(vb.ms()).spectralWindow().measFreqRef()(spw));
422 423 ​
423 -
424 424 ​
425 - // The input frequencies
425 +​
426 + // The input frequencies
426 427 Vector<Int> chanNums=vb.getChannelNumbers(0);
427 -
428 +​
428 429 Vector<Double> inFreq(chanNums.nelements());
429 430 Vector<Double> spwfreqs=MSColumns(vb.ms()).spectralWindow().chanFreq().get(spw);
430 431 for (uInt k=0; k < chanNums.nelements(); ++k){
431 432 ​
432 433 inFreq[k]=spwfreqs[chanNums[k]];
433 434 }
434 435 ​
435 436 // The output frequencies
436 437 outFreq.resize(inFreq.nelements());
437 438 ​
451 452 ​
452 453 // Reset the timestamp (ASSUMES TIME is constant in the VisBuffer)
453 454 mframe_.resetEpoch(vb.time()(0)/86400.0);
454 455 ​
455 456 // Reset the direction (ASSUMES phaseCenter is constant in the VisBuffer)
456 457 mframe_.resetDirection(vb.phaseCenter());
457 458 ​
458 459 // cout << "Frame = " << mframe_ << endl;
459 460 ​
460 461 // The conversion engine:
461 - MFrequency::Convert toNewFrame(obsMFreqType,
462 + MFrequency::Convert toNewFrame(obsMFreqType,
462 463 MFrequency::Ref(newMFreqType, mframe_));
463 464 ​
464 465 // Do the conversion
465 466 for (uInt k=0; k< inFreq.nelements(); ++k)
466 467 outFreq(k)=toNewFrame(inFreq(k)).getValue().getValue();
467 -
468 +​
468 469 }
469 470 else{
470 471 // The requested frame is the same as the observed frame
471 472 outFreq=inFreq;
472 473 }
473 474 ​
474 475 //cerr << std::setprecision(9) << " infreq " << inFreq[152] << " " << outFreq[152] << " vb freq " << vb.getFrequencies(0, freqFrame)[152] << endl;
475 476 ​
476 477 }
477 478 ​
478 - void VisBufferUtil::toVelocity(Vector<Double>& outVel,
479 - const VisBuffer& vb,
479 + void VisBufferUtil::toVelocity(Vector<Double>& outVel,
480 + const VisBuffer& vb,
480 481 const MFrequency::Types freqFrame,
481 482 const MVFrequency restFreq,
482 483 const MDoppler::Types veldef){
483 484 ​
484 485 // The input frequencies (a reference)
485 486 Vector<Double> inFreq(vb.frequency());
486 487 ​
487 488 // The output velocities
488 489 outVel.resize(inFreq.nelements());
489 490 ​
490 491 // Reset the timestamp (ASSUMES TIME is constant in the VisBuffer)
491 492 mframe_.resetEpoch(vb.time()(0)/86400.0);
492 -
493 +​
493 494 // Reset the direction (ASSUMES phaseCenter is constant in the VisBuffer)
494 495 //mframe_.resetDirection(vb.phaseCenter());
495 496 mframe_.resetDirection(vb.msColumns().field().phaseDirMeasCol()(vb.fieldId())(IPosition(1,0)));
496 -
497 +​
497 498 // The frequency conversion engine:
498 499 Int spw=vb.spectralWindow();
499 500 MFrequency::Types obsMFreqType=(MFrequency::Types)(vb.msColumns().spectralWindow().measFreqRef()(spw));
500 501 ​
501 502 MFrequency::Types newMFreqType=freqFrame;
502 503 if (freqFrame==MFrequency::N_Types)
503 504 // Don't convert frame
504 505 newMFreqType=obsMFreqType;
505 506 ​
506 - MFrequency::Convert toNewFrame(obsMFreqType,
507 + MFrequency::Convert toNewFrame(obsMFreqType,
507 508 MFrequency::Ref(newMFreqType, mframe_));
508 509 ​
509 510 // The velocity conversion engine:
510 511 MDoppler::Ref dum1(MDoppler::RELATIVISTIC);
511 512 MDoppler::Ref dum2(veldef);
512 513 MDoppler::Convert dopConv(dum1, dum2);
513 514 ​
514 515 // Cope with unspecified rest freq
515 516 MVFrequency rf=restFreq;
516 517 if (restFreq.getValue()<=0.0)
619 620 MDoppler eh2 = dopConv(eh);
620 621 outVel(k)=eh2.getValue().get().getValue();
621 622 }
622 623 ​
623 624 ​
624 625 ​
625 626 ​
626 627 }
627 628 ​
628 629 MDirection VisBufferUtil::getPointingDir(const VisBuffer& vb, const Int antid, const Int vbrow){
629 -
630 +​
630 631 Timer tim;
631 632 tim.mark();
632 633 const MSColumns& msc=vb.msColumns();
633 634 //cerr << "oldMSId_p " << oldMSId_p << " vb " << vb.msId() << endl;
634 635 if(vb.msId() < 0)
635 636 throw(AipsError("VisBuffer is not attached to an ms so cannot get pointing "));
636 637 if(oldMSId_p != vb.msId()){
637 638 oldMSId_p=vb.msId();
638 639 if(timeAntIndex_p.shape()(0) < (oldMSId_p+1)){
639 640 timeAntIndex_p.resize(oldMSId_p+1, true);
661 662 Vector<Int> antcol;
662 663 mspc.time().getColumn(timecol, True);
663 664 mspc.interval().getColumn(intervalcol, True);
664 665 mspc.antennaId().getColumn(antcol, True);
665 666 Double *tcolptr=timecol.getStorage(timcolstor);
666 667 Double *intcolptr=intervalcol.getStorage(intcolstor);
667 668 Int * antcolptr=antcol.getStorage(antcolstor);
668 669 Int npointrow=msc.pointing().nrow();
669 670 #pragma omp parallel for firstprivate(nTimes, tuniqptr, tcolptr, antcolptr, intcolptr, npointrow), shared(mspc)
670 671 for (uInt a=0; a < nAnt; ++a){
671 -
672 +​
672 673 //Double wtime1=omp_get_wtime();
673 674 Vector<ssize_t> indices;
674 675 Vector<MDirection> theDirs(nTimes);
675 676 pointingIndex(tcolptr, antcolptr, intcolptr, npointrow, a, nTimes, tuniqptr, indices);
676 -
677 +​
677 678 #pragma omp critical
678 679 {
679 680 for (uInt k=0; k <nTimes; ++k){
680 -
681 -
681 +​
682 +​
682 683 // std::ostringstream oss;
683 684 // oss.precision(13);
684 685 // oss << tuniqptr[k] << "_" << a;
685 686 // String key=oss.str();
686 687 std::pair<double, int> key=make_pair(t[uniqIndx[k]],a);
687 688 timeAntIndex_p[oldMSId_p][key]=indices[k] > -1 ? cshape : -1;
688 -
689 +​
689 690 if(indices[k] >-1){
690 -
691 +​
691 692 cachedPointingDir_p[oldMSId_p][cshape]=mspc.directionMeas(indices[k]);
692 693 cshape+=1;
693 694 }
694 -
695 -
695 +​
696 +​
696 697 }
697 698 }//end critical
698 -
699 -
699 +​
700 +​
700 701 }
701 -
702 +​
702 703 cachedPointingDir_p[oldMSId_p].resize(cshape, True);
703 704 }
704 -
705 +​
705 706 }
706 707 ​
707 708 /////
708 709 // String index=String::toString(vb.time()(vbrow))+String("_")+String::toString(antid);
709 710 // std::ostringstream oss;
710 711 // oss.precision(13);
711 712 // oss << vb.time()(vbrow) << "_" << antid ;
712 713 // String index=oss.str();
713 714 //cerr << "index "<< index << endl;
714 715 ////////////
781 782 mspc.interval().getColumn(intervalcol, True);
782 783 mspc.antennaId().getColumn(antcol, True);
783 784 Double *tcolptr=timecol.getStorage(timcolstor);
784 785 Double *intcolptr=intervalcol.getStorage(intcolstor);
785 786 Int * antcolptr=antcol.getStorage(antcolstor);
786 787 Int npointrow=vb.ms().pointing().nrow();
787 788 //ofstream myfile;
788 789 //myfile.open ("POINTING.txt", ios::trunc);
789 790 #pragma omp parallel for firstprivate(nTimes, tuniqptr, tcolptr, antcolptr, intcolptr, npointrow, timeType, timeUnit, pos, dirframe), shared(mspc)
790 791 for (uInt a=0; a < nAnt; ++a){
791 -
792 +​
792 793 //Double wtime1=omp_get_wtime();
793 794 Vector<ssize_t> indices;
794 795 //Vector<MDirection> theDirs(nTimes);
795 796 pointingIndex(tcolptr, antcolptr, intcolptr, npointrow, a, nTimes, tuniqptr, indices);
796 -
797 +​
797 798 #pragma omp critical
798 799 {
799 800 MEpoch timenow(Quantity(tuniqptr[0], timeUnit),timeType);
800 801 MeasFrame mframe(timenow, pos);
801 802 MDirection::Convert cvt(MDirection(), MDirection::Ref(dirframe, mframe));
802 803 for (uInt k=0; k <nTimes; ++k){
803 -
804 -
804 +​
805 +​
805 806 /*std::ostringstream oss;
806 807 oss.precision(13);
807 808 oss << tuniqptr[k] << "_" << a;
808 809 String key=oss.str();
809 810 */
810 811 // myfile <<std::setprecision(13) << tuniqptr[k] << "_" << a << " index "<< indices[k] << "\n";
811 812 pair<double, int> key=make_pair(tuniqptr[k],a);
812 813 timeAntIndex_p[oldMSId_p][key]=indices[k] > -1 ? cshape : -1;
813 -
814 +​
814 815 if(indices[k] >-1){
815 816 timenow=MEpoch(Quantity(tuniqptr[k], timeUnit),timeType);
816 817 mframe.resetEpoch(timenow);
817 818 cachedPointingDir_p[oldMSId_p][cshape]=cvt(mspc.directionMeas(indices[k]));
818 -
819 +​
819 820 cshape+=1;
820 821 }
821 -
822 -
822 +​
823 +​
823 824 }
824 825 }//end critical
825 -
826 -
826 +​
827 +​
827 828 }
828 829 cachedPointingDir_p[oldMSId_p].resize(cshape, True);
829 830 }
830 -
831 +​
831 832 }
832 833 ​
833 834 /////
834 835 // String index=String::toString(vb.time()(vbrow))+String("_")+String::toString(antid);
835 836 /* std::ostringstream oss;
836 837 oss.precision(13);
837 838 oss << vb.time()(vbrow) << "_" << antid ;
838 839 String index=oss.str();
839 840 */
840 841 pair<double, int> index=make_pair(vb.time()(vbrow),antid);
841 842 rowincache=timeAntIndex_p[oldMSId_p].at(index);
842 843 ​
843 844 //tim.show("retrieved cache");
844 845 }///if usepointing
845 846 if(rowincache <0)
846 847 return getPhaseCenter(vb);
847 848 return cachedPointingDir_p[oldMSId_p][rowincache];
848 849 ​
849 850 ​
850 851 ​
851 852 }
852 -
853 +​
853 854 void VisBufferUtil::pointingIndex(Double*& timecol, Int*& antcol, Double*& intervalcol, const rownr_t nrow, const Int ant, const Int ntimes, Double*& ptime, Vector<ssize_t>& indices){
854 -
855 +​
855 856 indices.resize(ntimes);
856 -
857 +​
857 858 indices.set(-1);
858 -
859 +​
859 860 for(Int pt=0; pt < ntimes; ++pt){
860 861 //cerr << " " << guessRow ;
861 -
862 +​
862 863 /*for (Int k=0; k< 2; ++k){
863 864 Int start=guessRow;
864 865 Int end=nrow;
865 866 if(k==1){
866 867 start=0;
867 868 end=guessRow;
868 869 }
869 870 */
870 871 Double nearval=1e99;
871 872 ssize_t nearestIndx=-1;
874 875 for (ssize_t i=start; i<end; i++) {
875 876 if(intervalcol[i]<=0.0 && ant==antcol[i]){
876 877 if(abs(timecol[i]-ptime[pt]) < nearval){
877 878 nearestIndx=i;
878 879 nearval=abs(timecol[i]-ptime[pt]);
879 880 }
880 881 }
881 882 }
882 883 indices[pt]=nearestIndx;
883 884 ​
884 -
885 +​
885 886 for (ssize_t i=start; i<end; i++) {
886 887 if(ant == antcol[i]){
887 888 Double halfInt=0.0;
888 889 Bool done=False;
889 890 if(intervalcol[i]<=0.0){
890 891 // if(abs(timecol[inx[i]]-ptime[pt]) < nearval){
891 892 // nearestIndx=inx[i];
892 893 // }
893 894 ssize_t counter=0;
894 895 ssize_t adder=1;
895 896 done=False;
896 897 // while(!( (timecol[i+counter]!=timecol[i]))){
897 898 while(!done){
898 899 counter=counter+adder;
899 900 if((ssize_t)nrow <= i+counter){
900 - adder=-1;
901 + adder=-1;
901 902 counter=0;
902 903 }
903 904 ////Could not find another point (interval is infinite) hence only 1 valid point
904 905 if( (i+counter) < 0){
905 906 cerr << "HIT BREAK " << endl;
906 907 indices[pt]=i;
907 908 break;
908 909 }
909 910 if( (antcol[i+counter]==ant && timecol[i+counter] != timecol[i]) ){
910 911 done=True;
911 912 }
912 913 }
913 -
914 +​
914 915 //if(ant==12 && abs(timecol[i+counter]-timecol[i]) > 10.0){
915 916 //cerr << "i " << i << " counter " << counter << " done " << done << " adder " << adder << "ant count "<< antcol[i+counter] << " diff " << abs(timecol[i+counter]-timecol[i]) << endl;
916 917 // }
917 918 halfInt = abs(timecol[i+counter]-timecol[i])/2.0;
918 919 }
919 920 else{
920 921 halfInt = intervalcol[i]/2.0;
921 922 }
922 923 if (halfInt>0.0) {
923 -
924 +​
924 925 if ((timecol[i] >= (ptime[pt] - halfInt)) && (timecol[i] <= (ptime[pt] + halfInt))) {
925 926 ////TESTOO
926 927 //if(ant==12){
927 928 // cerr << "timecol " << timecol[i] << " halfInt " << halfInt << " TEST " << timecol[nearestIndx] << " inx " << i << " " << nearestIndx << endl;
928 929 //}
929 930 indices[pt]=abs(timecol[i]-ptime[pt]) < nearval ? i : nearestIndx;
930 931 ////////TESTOO
931 932 if(indices[pt] > 4688000){
932 933 cerr << indices[pt] << " timecol " << timecol[i] << " halfInt " << halfInt << " TEST " << timecol[nearestIndx] << " nearval " << nearval << " inx " << i << " " << nearestIndx << endl;
933 -
934 +​
934 935 }
935 936 ///////////////////
936 937 break;
937 938 }
938 939 } else {
939 940 // valid for all times (we should also handle interval<0 -> timestamps)
940 941 cerr << "JUMPY " << i << " ant " << ant << " halfint " << halfInt << " done "<< done << endl;
941 942 indices[pt]=i;
942 943 break;
943 944 }
944 945 ​
945 946 }//if ant
946 947 }//start end
947 -
948 +​
948 949 //}//k
949 -
950 +​
950 951 }//pt
951 -
952 +​
952 953 //cerr << "ant " << ant << " indices " << indices << endl;
953 954 }
954 955 ​
955 956 MDirection VisBufferUtil::getPhaseCenter(const vi::VisBuffer2& vb, const Double timeo){
956 957 //Timer tim;
957 -
958 - Double timeph = timeo > 0.0 ? timeo : vb.time()(0);
958 +​
959 + Double timeph = timeo > 0.0 ? timeo : vb.time()(0);
959 960 //MDirection outdir;
960 961 if(oldPCMSId_p != vb.msId()){
961 962 MSColumns msc(vb.ms());
962 963 //tim.mark();
963 964 //cerr << "MSID: "<< oldPCMSId_p << " " << vb.msId() << endl;
964 965 oldPCMSId_p=vb.msId();
965 966 if(cachedPhaseCenter_p.shape()(0) < (oldPCMSId_p+1)){
966 967 cachedPhaseCenter_p.resize(oldPCMSId_p+1, true);
967 968 cachedPhaseCenter_p[oldPCMSId_p]=map<Double, MDirection>();
968 969 }
969 970 if( cachedPhaseCenter_p[oldPCMSId_p].empty()){
970 971 Vector<Double> tOrig;
971 972 msc.time().getColumn(tOrig);
972 973 Vector<Int> fieldId;
973 974 msc.fieldId().getColumn(fieldId);
974 975 Vector<Double> t;
975 976 Vector<Int> origindx;
976 977 rejectConsecutive(tOrig, t, origindx);
977 978 Vector<uInt> uniqIndx;
978 -
979 +​
979 980 uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, t, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
980 981 //cerr << "ntimes " << nTimes << " " << uniqIndx << "\n origInx " << origindx << endl;
981 982 const MSFieldColumns& msfc=msc.field();
982 983 for (uInt k=0; k <nTimes; ++k){
983 984 //cerr << t[uniqIndx[k]] << " " << fieldId[origindx[uniqIndx[k]]] << endl;
984 985 //cerr << msfc.phaseDirMeas(fieldId[origindx[uniqIndx[k]]], t[uniqIndx[k]]) << endl;
985 986 //cerr << "size " << cachedPhaseCenter_p[oldPCMSId_p].size() << endl;
986 987 String ephemIfAny=msfc.ephemPath(fieldId[origindx[uniqIndx[k]]]);
987 988 if(ephemIfAny=="" || !Table::isReadable(ephemIfAny, False)){
988 989 (cachedPhaseCenter_p[oldPCMSId_p])[t[uniqIndx[k]]]=msfc.phaseDirMeas(fieldId[origindx[uniqIndx[k]]], t[uniqIndx[k]]);
989 990 }
990 991 else{
991 992 Vector<MDirection> refDir(msfc.referenceDirMeasCol()(fieldId[origindx[uniqIndx[k]]]));
992 993 (cachedPhaseCenter_p[oldPCMSId_p])[t[uniqIndx[k]]]=getEphemBasedPhaseDir(vb, ephemIfAny, refDir(0), t[uniqIndx[k]]);
993 994 }
994 -
995 +​
995 996 }
996 -
997 +​
997 998 ​
998 999 }
999 1000 //tim.show("After caching all phasecenters");
1000 1001 }
1001 1002 //tim.mark();
1002 1003 MDirection retval;
1003 1004 auto it=cachedPhaseCenter_p[oldPCMSId_p].find(timeph);
1004 1005 if(it != cachedPhaseCenter_p[oldPCMSId_p].end()){
1005 1006 retval=it->second;
1006 1007 }
1018 1019 retval=low->second;
1019 1020 }
1020 1021 else{
1021 1022 retval=upp->second;
1022 1023 }
1023 1024 ​
1024 1025 }
1025 1026 }
1026 1027 //tim.show("retrieved cache");
1027 1028 //cerr << std::setprecision(12) <<"msid " << oldPCMSId_p << " time "<< timeph << " val " << retval.toString() << endl;
1028 -
1029 +​
1029 1030 return retval;
1030 1031 ​
1031 1032 ​
1032 1033 ​
1033 1034 }
1034 1035 ​
1035 1036 ​
1036 1037 MDirection VisBufferUtil::getEphemBasedPhaseDir(const vi::VisBuffer2& vb, const String& ephemPath, const MDirection&refDir, const Double t){
1037 - MEpoch ep(Quantity(t, "s"), ROMSColumns(vb.ms()).timeMeas()(0).getRef());
1038 + MEpoch ep(Quantity(t, "s"), vb.getVi()->getImpl()->getEpoch().getRef());
1038 1039 mframe_.resetEpoch(ep);
1039 1040 if(!Table::isReadable(ephemPath, False))
1040 1041 return refDir;
1041 1042 MeasComet mcomet(Path(ephemPath).absoluteName());
1042 1043 mframe_.set(mcomet);
1043 1044 MDirection::Ref outref1(MDirection::AZEL, mframe_);
1044 1045 MDirection tmpazel=MDirection::Convert(MDirection(MDirection::COMET), outref1)();
1045 1046 MDirection::Types outtype=(MDirection::Types) refDir.getRef().getType();
1046 1047 MDirection::Ref outref(outtype, mframe_);
1047 1048 MDirection outdir=MDirection::Convert(tmpazel, outref)();
1048 1049 MVDirection mvoutdir(outdir.getAngle());
1049 1050 MVDirection mvrefdir(refDir.getAngle());
1050 1051 //copying what ROMSFieldColumns::extractDirMeas does
1051 1052 mvoutdir.shift(mvrefdir.getAngle(), True);
1052 1053 return MDirection(mvoutdir, outtype);
1053 1054 }
1054 -
1055 - MDirection VisBufferUtil::getEphemDir(const vi::VisBuffer2& vb,
1055 +​
1056 + MDirection VisBufferUtil::getEphemDir(const vi::VisBuffer2& vb,
1056 1057 const Double timeo){
1057 1058 ​
1058 - Double timeEphem = timeo > 0.0 ? timeo : vb.time()(0);
1059 - MSColumns msc(vb.ms());
1060 - const MSFieldColumns& msfc=msc.field();
1059 + Double timeEphem = timeo > 0.0 ? timeo : vb.time()(0);
1060 + const MSFieldColumns &msfc = vb.subtableColumns().field();
1061 1061 Int fieldId=vb.fieldId()(0);
1062 1062 MDirection refDir(Quantity(0, "deg"), Quantity(0, "deg"),(MDirection::Types)msfc.ephemerisDirMeas(fieldId, timeEphem).getRef().getType());
1063 1063 //Now do the parallax correction
1064 1064 return getEphemBasedPhaseDir(vb, msfc.ephemPath(fieldId), refDir, timeEphem);
1065 1065 ​
1066 1066 ​
1067 1067 }
1068 -
1068 +​
1069 1069 //utility to reject consecutive similar value for sorting
1070 1070 void VisBufferUtil::rejectConsecutive(const Vector<Double>& t, Vector<Double>& retval){
1071 1071 uInt n=t.nelements();
1072 1072 if(n >0){
1073 1073 retval.resize(n);
1074 1074 retval[0]=t[0];
1075 1075 }
1076 1076 else
1077 1077 return;
1078 1078 Int prev=0;
1090 1090 uInt n=t.nelements();
1091 1091 if(n >0){
1092 1092 retval.resize(n);
1093 1093 indx.resize(n);
1094 1094 retval[0]=t[0];
1095 1095 indx[0]=0;
1096 1096 }
1097 1097 else
1098 1098 return;
1099 1099 Int prev=0;
1100 - for (uInt k=1; k < n; ++k){
1100 + for (uInt k=1; k < n; ++k){
1101 1101 if(t[k] != retval(prev)){
1102 1102 ++prev;
1103 1103 //retval.resize(prev+1, true);
1104 1104 retval[prev]=t[k];
1105 1105 //indx.resize(prev+1, true);
1106 1106 indx[prev]=k;
1107 1107 }
1108 1108 }
1109 1109 retval.resize(prev+1, true);
1110 1110 indx.resize(prev+1, true);
1111 -
1111 +​
1112 1112 }
1113 1113 ​
1114 1114 // helper function to swap the y and z axes of a Cube
1115 1115 void VisBufferUtil::swapyz(Cube<Complex>& out, const Cube<Complex>& in)
1116 1116 {
1117 1117 IPosition inShape=in.shape();
1118 1118 uInt nxx=inShape(0),nyy=inShape(2),nzz=inShape(1);
1119 - //resize breaks references...so out better have the right shape
1119 + //resize breaks references...so out better have the right shape
1120 1120 //if references is not to be broken
1121 1121 if(out.nelements()==0)
1122 1122 out.resize(nxx,nyy,nzz);
1123 1123 Bool deleteIn,deleteOut;
1124 1124 const Complex* pin = in.getStorage(deleteIn);
1125 1125 Complex* pout = out.getStorage(deleteOut);
1126 1126 uInt i=0, zOffset=0;
1127 1127 for (uInt iz=0; iz<nzz; ++iz, zOffset+=nxx) {
1128 1128 Int yOffset=zOffset;
1129 1129 for (uInt iy=0; iy<nyy; ++iy, yOffset+=nxx*nzz) {
1130 - for (uInt ix=0; ix<nxx; ++ix){
1130 + for (uInt ix=0; ix<nxx; ++ix){
1131 1131 pout[i++] = pin[ix+yOffset];
1132 1132 }
1133 1133 }
1134 1134 }
1135 1135 out.putStorage(pout,deleteOut);
1136 1136 in.freeStorage(pin,deleteIn);
1137 1137 }
1138 1138 ​
1139 1139 // helper function to swap the y and z axes of a Cube
1140 1140 void VisBufferUtil::swapyz(Cube<Bool>& out, const Cube<Bool>& in)
1150 1150 for (uInt iz=0; iz<nzz; iz++, zOffset+=nxx) {
1151 1151 Int yOffset=zOffset;
1152 1152 for (uInt iy=0; iy<nyy; iy++, yOffset+=nxx*nzz) {
1153 1153 for (uInt ix=0; ix<nxx; ix++) pout[i++] = pin[ix+yOffset];
1154 1154 }
1155 1155 }
1156 1156 out.putStorage(pout,deleteOut);
1157 1157 in.freeStorage(pin,deleteIn);
1158 1158 }
1159 1159 ​
1160 -
1160 +​
1161 1161 ​
1162 1162 ​
1163 1163 } //# NAMESPACE CASA - END
1164 1164 ​

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut