Commits

Ville Suoranta authored 496bd6166f9 Merge
Merge pull request #453 in CASA/casa from feature/CAS-9004 to master

* commit 'd3c678955ca35fe50461fe831333c26ac61451a3': added some extra notes on phasecenter parameter Added TRACKFIELD option for mosaic Put some logging info about tracking added mosaicftnew to check for mosaic prevent using TRACKFIELD with mosaic Forgot to shift the second antenna beam by the same amount as the first. Fixed that now. Fixed an error for large angle shift for moving sources added the special case of phasecenter='TRACKFIELD' to track on the moving phasecenter of the FIELD subtables of the MSs. seriously let's get openmp on the mac added some description to phasecenter added support for tracking to mosaic and cleaned up and sped up some code made makepb shift beam is tracking a moving source Fixed an issue with the fact that ephemeris_id is an optional column Brought up interface and modified some the code to track ephemerides tables object (CAS-9004)
No tags

code/msvis/MSVis/VisBufferUtil.cc

Modified
32 32 #include <casa/Exceptions/Error.h>
33 33 #include <casa/Arrays/ArrayMath.h>
34 34 #include <casa/Arrays/MatrixMath.h>
35 35 #include <casa/Arrays/Array.h>
36 36 #include <casa/Arrays/Vector.h>
37 37 #include <casa/Arrays/Matrix.h>
38 38 #include <casa/Arrays/Cube.h>
39 39 #include <casa/OS/Timer.h>
40 40 #include <measures/Measures/UVWMachine.h>
41 41 #include <measures/Measures/MeasTable.h>
42 +#include <ms/MSSel/MSSelectionTools.h>
42 43 #include <ms/MeasurementSets/MSPointingColumns.h>
43 44 #include <msvis/MSVis/VisBufferUtil.h>
44 45 #include <msvis/MSVis/StokesVector.h>
45 46 #include <msvis/MSVis/VisibilityIterator.h>
46 47 #include <msvis/MSVis/VisibilityIterator2.h>
47 48 #include <msvis/MSVis/VisBuffer.h>
48 49 #include <ms/MeasurementSets/MSColumns.h>
49 50 #include <casa/iostream.h>
50 51 #include <iomanip>
52 +#ifdef _OPENMP
53 +#include <omp.h>
54 +#endif
51 55 using namespace std;
52 56
53 57 using namespace casacore;
54 58 namespace casa { //# NAMESPACE CASA - BEGIN
55 59
56 60 // <summary>
57 61 // </summary>
58 62
59 63 // <reviewed reviewer="" date="" tests="tMEGI" demos="">
60 64
632 636 const ROMSPointingColumns& mspc=vb.msColumns().pointing();
633 637 Int guessIndex=0;
634 638 for (uInt k=0; k <nTimes; ++k){
635 639 for (uInt a=0; a < nAnt; ++a){
636 640 std::ostringstream oss;
637 641 oss.precision(13);
638 642 oss << t[uniqIndx[k]] << "_" << a;
639 643 String key=oss.str();
640 644 //String key=String::toString(t[uniqIndx[k]])+String("_")+String::toString(a);
641 645 Int row=mspc.pointingIndex(a, t[uniqIndx[k]], guessIndex);
642 - cerr << "String "<< key << "pointing row "<< row << endl;
646 + //cerr << "String "<< key << "pointing row "<< row << endl;
643 647 timeAntIndex_p[oldMSId_p][key]=row > -1 ? cachedPointingDir_p[oldMSId_p].shape()[0] : -1;
644 648 guessIndex=row;
645 649 if(row >-1){
646 650 cachedPointingDir_p[oldMSId_p].resize(cachedPointingDir_p[oldMSId_p].nelements()+1, true);
647 651 cachedPointingDir_p[oldMSId_p][cachedPointingDir_p[oldMSId_p].nelements()-1]=mspc.directionMeas(row);
648 652 }
649 653
650 654 }
651 655 }
652 656
653 657 }
654 658 tim.show("After caching all ant pointings");
655 659 }
656 660
657 661 /////
658 662 // String index=String::toString(vb.time()(vbrow))+String("_")+String::toString(antid);
659 663 std::ostringstream oss;
660 664 oss.precision(13);
661 665 oss << vb.time()(vbrow) << "_" << antid ;
662 666 String index=oss.str();
663 667 Int rowincache=timeAntIndex_p[oldMSId_p][index];
664 - cerr << "key "<< index << " index " << rowincache << endl;
668 + //cerr << "key "<< index << " index " << rowincache << endl;
665 669 tim.show("retrieved cache");
666 670 if(rowincache <0)
667 671 return vb.phaseCenter();
668 672 return cachedPointingDir_p[oldMSId_p][rowincache];
669 673
670 674
671 675
672 676 }
673 677 MDirection VisBufferUtil::getPointingDir(const vi::VisBuffer2& vb, const Int antid, const Int vbrow){
674 678 Timer tim;
675 679 tim.mark();
676 680 ROMSColumns msc(vb.ms());
677 - //MDirection outdir;
678 681 if(oldMSId_p != vb.msId()){
679 - tim.mark();
680 - //cerr << "MSID: "<< oldMSId_p << " " << vb.msId() << endl;
681 682 oldMSId_p=vb.msId();
682 683 if(timeAntIndex_p.shape()(0) < (oldMSId_p+1)){
683 684 timeAntIndex_p.resize(oldMSId_p+1, true);
684 685 cachedPointingDir_p.resize(oldMSId_p+1, true);
685 686 }
686 687 if( timeAntIndex_p[oldMSId_p].empty()){
687 688 Vector<Double> tOrig;
688 689 msc.time().getColumn(tOrig);
689 690 Vector<Double> t;
690 691 rejectConsecutive(tOrig, t);
691 692 Vector<uInt> uniqIndx;
692 693 uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, t, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
693 694 uInt nAnt=msc.antenna().nrow();
694 695 const ROMSPointingColumns& mspc=msc.pointing();
695 - Int guessIndex=0;
696 + Vector<Double> tUniq(nTimes);
696 697 for (uInt k=0; k <nTimes; ++k){
697 - for (uInt a=0; a < nAnt; ++a){
698 + tUniq[k]= t[uniqIndx[k]];
699 + }
700 + Bool tstor, timcolstor, intcolstor, antcolstor;
701 + Double * tuniqptr=tUniq.getStorage(tstor);
702 + Int cshape=cachedPointingDir_p[oldMSId_p].shape()[0];
703 + cachedPointingDir_p[oldMSId_p].resize(cshape+nTimes*nAnt, True);
704 + Vector<Double> timecol;
705 + Vector<Double> intervalcol;
706 + Vector<Int> antcol;
707 + mspc.time().getColumn(timecol, True);
708 + mspc.interval().getColumn(intervalcol, True);
709 + mspc.antennaId().getColumn(antcol, True);
710 + Double *tcolptr=timecol.getStorage(timcolstor);
711 + Double *intcolptr=intervalcol.getStorage(intcolstor);
712 + Int * antcolptr=antcol.getStorage(antcolstor);
713 + Int npointrow=vb.ms().pointing().nrow();
714 +#pragma omp parallel for firstprivate(nTimes, tuniqptr, tcolptr, antcolptr, intcolptr, npointrow), shared(mspc)
715 + for (uInt a=0; a < nAnt; ++a){
716 +
717 + //Double wtime1=omp_get_wtime();
718 + Vector<Int> indices;
719 + Vector<MDirection> theDirs(nTimes);
720 + pointingIndex(tcolptr, antcolptr, intcolptr, npointrow, a, nTimes, tuniqptr, indices);
721 +
722 +#pragma omp critical
723 + {
724 + for (uInt k=0; k <nTimes; ++k){
725 +
726 +
698 727 std::ostringstream oss;
699 728 oss.precision(13);
700 - oss << t[uniqIndx[k]] << "_" << a;
729 + oss << tuniqptr[k] << "_" << a;
701 730 String key=oss.str();
702 - //String key=String::toString(t[uniqIndx[k]])+String("_")+String::toString(a);
703 - Int row=mspc.pointingIndex(a, t[uniqIndx[k]], guessIndex);
704 - //cerr << "String "<< key << " pointing row "<< row << endl;
705 - timeAntIndex_p[oldMSId_p][key]=row > -1 ? cachedPointingDir_p[oldMSId_p].shape()[0] : -1;
706 - guessIndex=row;
707 - if(row >-1){
708 - cachedPointingDir_p[oldMSId_p].resize(cachedPointingDir_p[oldMSId_p].nelements()+1, true);
709 - cachedPointingDir_p[oldMSId_p][cachedPointingDir_p[oldMSId_p].nelements()-1]=mspc.directionMeas(row);
731 +
732 + timeAntIndex_p[oldMSId_p][key]=indices[k] > -1 ? cshape : -1;
733 +
734 + if(indices[k] >-1){
735 +
736 + cachedPointingDir_p[oldMSId_p][cshape]=mspc.directionMeas(indices[k]);
737 + cshape+=1;
710 738 }
739 +
740 +
741 + }
742 + }//end critical
711 743
712 - }
744 +
713 745 }
714 -
746 +
747 + cachedPointingDir_p[oldMSId_p].resize(cshape, True);
715 748 }
716 - tim.show("After caching all ant pointings");
749 +
717 750 }
718 751
719 752 /////
720 753 // String index=String::toString(vb.time()(vbrow))+String("_")+String::toString(antid);
721 754 std::ostringstream oss;
722 755 oss.precision(13);
723 756 oss << vb.time()(vbrow) << "_" << antid ;
724 757 String index=oss.str();
725 758 Int rowincache=timeAntIndex_p[oldMSId_p][index];
726 - //cerr << "key "<< index << " index " << rowincache << endl;
759 + ///////TESTOO
760 + /* if(rowincache>=0){
761 + cerr << "msid " << oldMSId_p << " key "<< index << " index " << rowincache<< " " << cachedPointingDir_p[oldMSId_p][rowincache] << endl;
762 + }*/
763 + /////////////
727 764 //tim.show("retrieved cache");
728 765 if(rowincache <0)
729 766 return vb.phaseCenter();
730 767 return cachedPointingDir_p[oldMSId_p][rowincache];
731 768
732 769
733 770
734 771 }
735 772
773 + void VisBufferUtil::pointingIndex(Double*& timecol, Int*& antcol, Double*& intervalcol, const Int nrow, const Int ant, const Int ntimes, Double*& ptime, Vector<Int>& indices){
774 +
775 + indices.resize(ntimes);
776 +
777 + indices.set(-1);
778 + Int guessRow=0;
779 +
780 + for(Int pt=0; pt < ntimes; ++pt){
781 + //cerr << " " << guessRow ;
782 + for (Int k=0; k< 2; ++k){
783 + Int start=guessRow;
784 + Int end=nrow;
785 + if(k==1){
786 + start=0;
787 + end=guessRow;
788 + }
789 + for (Int i=start; i<end; i++) {
790 + if(ant == antcol[i]){
791 + Double halfInt=0.0;
792 + if(intervalcol[i]==0.0){
793 + Int counter=0;
794 + Int adder=1;
795 +
796 + while(!( (timecol[i+counter]!=timecol[i]))){
797 + counter=counter+adder;
798 + if(nrow <= i+counter){
799 + adder=-1;
800 + counter=0;
801 + }
802 + ////Could not find another point (interval is infinite) hence only 1 valid point
803 + if( (i+counter) < 0){
804 + indices[pt]=i;
805 + break;
806 + }
807 + }
808 + halfInt = abs(timecol[i+counter]-timecol[i])/2.0;
809 + }
810 + else{
811 + halfInt = intervalcol[i]/2.0;
812 + }
813 + if (halfInt>0.0) {
814 + if (timecol[i] >= ptime[pt] - halfInt && timecol[i] <= ptime[pt] + halfInt) {
815 + indices[pt]=i;
816 + guessRow=i;
817 + break;
818 + }
819 + } else {
820 + // valid for all times (we should also handle interval<0 -> timestamps)
821 + indices[pt]=i;
822 + guessRow=i;
823 + break;
824 + }
825 +
826 + }//if ant
827 + }//start end
828 +
829 + }//k
830 +
831 + }//pt
832 + //cerr << "ant " << ant << " indices " << indices << endl;
833 + }
834 +
736 835 MDirection VisBufferUtil::getPhaseCenter(const vi::VisBuffer2& vb, const Double timeo){
737 836 //Timer tim;
738 837
739 838 Double timeph = timeo > 0.0 ? timeo : vb.time()(0);
740 839 //MDirection outdir;
741 840 if(oldPCMSId_p != vb.msId()){
742 841 ROMSColumns msc(vb.ms());
743 842 //tim.mark();
744 843 //cerr << "MSID: "<< oldPCMSId_p << " " << vb.msId() << endl;
745 844 oldPCMSId_p=vb.msId();
791 890 if(fabs(timeph - (low->first)) < fabs(timeph - (upp->first))){
792 891 retval=low->second;
793 892 }
794 893 else{
795 894 retval=upp->second;
796 895 }
797 896
798 897 }
799 898 }
800 899 //tim.show("retrieved cache");
900 + //cerr << std::setprecision(12) <<"msid " << oldPCMSId_p << " time "<< timeph << " val " << retval.toString() << endl;
801 901
802 -
803 902 return retval;
804 903
805 904
806 905
807 906 }
808 907
908 + MDirection VisBufferUtil::getEphemDir(const vi::VisBuffer2& vb,
909 + const Double timeo){
910 +
911 + Double timeEphem = timeo > 0.0 ? timeo : vb.time()(0);
912 + ROMSColumns msc(vb.ms());
913 + const ROMSFieldColumns& msfc=msc.field();
914 + Int fieldId=vb.fieldId()(0);
915 + return msfc.ephemerisDirMeas(fieldId, timeEphem);
916 +
917 +
918 +
919 + }
809 920 //utility to reject consecutive similar value for sorting
810 921 void VisBufferUtil::rejectConsecutive(const Vector<Double>& t, Vector<Double>& retval){
811 922 uInt n=t.nelements();
812 923 if(n >0){
813 924 retval.resize(n);
814 925 retval[0]=t[0];
815 926 }
816 927 else
817 928 return;
818 929 Int prev=0;

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

Add shortcut