Task #9787: DPPP gaincal: consider all flags in solint

CEP/DP3/DPPP/include/DPPP/GainCal.h

@@ -95,14 +95,6 @@ namespace LOFAR {
       // Invert solution (for applying it)
       casa::Cube<casa::DComplex> invertSol(const casa::Cube<casa::DComplex>& sol);
 
-      // Counts the number of antennas with non-flagged data,
-      // Set a map for the used antennas in iS, returns the number of antennas
-      void setAntennaMaps (const casa::Bool* flag, uint freqCell);
-
-      // Remove rows and colums corresponding to antennas with too much
-      // flagged data from vis and mvis
-      void removeDeadAntennas ();
-
       // Fills the matrices itsVis and itsMVis
       void fillMatrices (casa::Complex* model, casa::Complex* data,
                          float* weight, const casa::Bool* flag);
@@ -153,11 +145,10 @@ namespace LOFAR {
 
       uint             itsMaxIter;
       double           itsTolerance;
-      bool             itsPropagateSolutions; // Not used currently, TODO: use this
+      bool             itsPropagateSolutions;
       uint             itsSolInt;  // Time cell size
       uint             itsNChan;   // Frequency cell size
       uint             itsNFreqCells;
-      uint             itsMinBLperAnt;
 
       uint             itsTimeSlotsPerParmUpdate;
       uint             itsConverged;

CEP/DP3/DPPP/src/GainCal.cc

@@ -83,7 +83,6 @@ namespace LOFAR {
         itsSolInt        (parset.getInt(prefix + "solint", 1)),
         itsNChan         (parset.getInt(prefix + "nchan", 0)),
         itsNFreqCells    (0),
-        itsMinBLperAnt   (parset.getInt(prefix + "minblperant", 4)),
         itsTimeSlotsPerParmUpdate
                          (parset.getInt(prefix + "timeslotsperparmupdate", 500)),
         itsConverged     (0),
@@ -329,7 +328,7 @@ namespace LOFAR {
         // Start new solution interval
 
         for (uint freqCell=0; freqCell<itsNFreqCells; freqCell++) {
-          setAntennaMaps(flag, freqCell);
+          iS[freqCell].clearStationFlagged();
           iS[freqCell].resetVis();
         }
       }
@@ -475,58 +474,6 @@ namespace LOFAR {
       }
     }
 
-    void GainCal::setAntennaMaps (const Bool* flag, uint freqCell) {
-      uint nCr=info().ncorr();
-      uint nCh=info().nchan();
-      uint nBl=info().nbaselines();
-
-      iS[freqCell].clearStationFlagged();
-
-      casa::Vector<casa::uInt> dataPerAntenna; // nAnt
-      dataPerAntenna.resize(info().antennaNames().size());
-
-      // TODO: implement smarter graph algorithm here that requires
-      // less than O(n_ant^3) worst case operations.
-      bool flaggedAntsAdded=true;
-      uint loopcount=0;
-      while (flaggedAntsAdded) {
-        dataPerAntenna=0;
-        for (uint bl=0;bl<nBl;++bl) {
-          uint ant1=info().getAnt1()[bl];
-          uint ant2=info().getAnt2()[bl];
-          if (ant1==ant2) {
-            continue;
-          }
-          uint chmax=min((freqCell+1)*itsNChan, nCh);
-          for (uint ch=freqCell*itsNChan;ch<chmax;++ch) {
-            for (uint cr=0;cr<nCr;++cr) {
-              if (!(flag[bl*nCr*nCh + ch*nCr + cr]
-                    || iS[freqCell].getStationFlagged()[ant1]
-                    || iS[freqCell].getStationFlagged()[ant2]
-                   )) {
-                dataPerAntenna(ant1)++;
-                dataPerAntenna(ant2)++;
-              }
-            }
-          }
-        }
-
-        flaggedAntsAdded=false;
-        for (uint ant=0; ant<info().antennaNames().size(); ++ant) {
-          if (dataPerAntenna(ant)>nCr*itsMinBLperAnt) {
-            iS[freqCell].getStationFlagged()[ant]=false; // Index in stefcal numbering
-          } else { // Not enough data
-            //cout<<"flagging station "<<ant<<", "<<dataPerAntenna(ant)<<endl;
-            if (!iS[freqCell].getStationFlagged()[ant]) {
-              iS[freqCell].getStationFlagged()[ant]=true;
-              flaggedAntsAdded=true;
-            }
-          }
-        }
-        loopcount++;
-      }
-    }
-
     void GainCal::stefcal () {
       itsTimerSolve.start();
 

CEP/DP3/DPPP/src/StefCal.cc

@@ -285,7 +285,13 @@ namespace LOFAR {
             vis_p=&_vis(IPosition(6,0,1,time,ch,1,st1)); for (uint st2=0;st2<_nSt;++st2) { t(3) += conj(_z(st2+_nSt*ch+_nSt*_nChan*time,3)) * vis_p[st2]; }// itsVis(IPosition(6,st2,1,time,ch,1,st1))
           }
         }
-        DComplex invdet= 1./(w(0) * w (3) - w(1)*w(2));
+        DComplex invdet= w(0) * w (3) - w(1)*w(2);
+        if (abs(invdet)==0) {
+          _stationFlagged[st1] = true;
+          _g(st1,0) = 0;
+          continue;
+        }
+        invdet= 1./invdet;
         _g(st1,0) = invdet * ( w(3) * t(0) - w(1) * t(2) );
         _g(st1,1) = invdet * ( w(3) * t(1) - w(1) * t(3) );
         _g(st1,2) = invdet * ( w(0) * t(2) - w(2) * t(0) );
@@ -333,6 +339,7 @@ namespace LOFAR {
         //cout<<", w="<<ww<<"       ";
         if (ww==0) {
           _stationFlagged[st1%_nSt]=true;
+          _g(st1,0)=0;
           continue;
         }
         _g(st1,0)=tt/ww;