Task #1892: optional Normalize Variance action added that can fix higher rates...

Task #1892: optional Normalize Variance action added that can fix higher rates of false positives at times the sky temperature is high

.gitattributes

@@ -398,6 +398,7 @@ CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/frequencyselectionaction.h
 CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/fringestopaction.h -text
 CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/imageraction.h -text
 CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/iterationaction.h -text
+CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/normalizevarianceaction.h -text
 CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/plotaction.h -text
 CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/quickcalibrateaction.h -text
 CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/rawappenderaction.h -text
@@ -642,6 +643,7 @@ CEP/DP3/AOFlagger/src/strategy/actions/foreachmsaction.cpp -text
 CEP/DP3/AOFlagger/src/strategy/actions/frequencyselectionaction.cpp -text
 CEP/DP3/AOFlagger/src/strategy/actions/fringestopaction.cpp -text
 CEP/DP3/AOFlagger/src/strategy/actions/imageraction.cpp -text
+CEP/DP3/AOFlagger/src/strategy/actions/normalizevarianceaction.cpp -text
 CEP/DP3/AOFlagger/src/strategy/actions/plotaction.cpp -text
 CEP/DP3/AOFlagger/src/strategy/actions/slidingwindowfitaction.cpp -text
 CEP/DP3/AOFlagger/src/strategy/actions/spatialcompositionaction.cpp -text

CEP/DP3/AOFlagger/include/AOFlagger/quality/defaultstatistics.h

@@ -101,6 +101,25 @@ class DefaultStatistics : public Serializable
 			return *this;
 		}
 		
+		DefaultStatistics ToSinglePolarization() const
+		{
+			if(_polarizationCount == 1)
+				return *this;
+			
+			DefaultStatistics singlePol(1);
+			for(unsigned p=0;p<_polarizationCount;++p)
+			{
+				singlePol.rfiCount[0] += rfiCount[p];
+				singlePol.count[0] += count[p];
+				singlePol.sum[0] += sum[p];
+				singlePol.sumP2[0] += sumP2[p];
+				singlePol.dCount[0] += dCount[p];
+				singlePol.dSum[0] += dSum[p];
+				singlePol.dSumP2[0] += dSumP2[p];
+			}
+			return singlePol;
+		}
+		
 		virtual void Serialize(std::ostream &stream) const
 		{
 			SerializeToUInt32(stream, _polarizationCount);

CEP/DP3/AOFlagger/include/AOFlagger/quality/statisticscollection.h

@@ -127,6 +127,11 @@ class StatisticsCollection : public Serializable
 			loadBaseline<false>(qualityData);
 		}
 		
+		void LoadTimeStatisticsOnly(QualityTablesFormatter &qualityData)
+		{
+			loadTime<false>(qualityData);
+		}
+		
 		void Add(QualityTablesFormatter &qualityData)
 		{
 			loadTime<true>(qualityData);

CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/action.h

@@ -52,6 +52,7 @@ namespace rfiStrategy {
 		FringeStopActionType,
 		ImagerActionType,
 		IterationBlockType,
+		NormalizeVarianceActionType,
 		PlotActionType,
 		QuickCalibrateActionType,
 		RawAppenderActionType,

CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/normalizevarianceaction.h

+/***************************************************************************
+ *   Copyright (C) 2008 by A.R. Offringa   *
+ *   offringa@astro.rug.nl   *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ *   This program is distributed in the hope that it will be useful,       *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU General Public License for more details.                          *
+ *                                                                         *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+#ifndef RFISTRATEGYNORMALIZEVARIANCEACTION_H
+#define RFISTRATEGYNORMALIZEVARIANCEACTION_H
+
+#include <map>
+
+#include <AOFlagger/strategy/actions/action.h>
+
+#include <AOFlagger/strategy/control/artifactset.h>
+
+namespace rfiStrategy {
+
+	/**
+		@author A.R. Offringa <offringa@astro.rug.nl>
+	*/
+	class NormalizeVarianceAction : public Action {
+		public:
+			NormalizeVarianceAction() :
+				_medianFilterSizeInS(60.0*15.0)
+			{
+			}
+			virtual ~NormalizeVarianceAction()
+			{
+			}
+			virtual std::string Description()
+			{
+				return "Normalize variance over time";
+			}
+			virtual void Perform(ArtifactSet &artifacts, ProgressListener &progress);
+
+			virtual ActionType Type() const { return NormalizeVarianceActionType; }
+
+		private:
+			double _medianFilterSizeInS;
+			
+			void correctData(std::vector<Image2DPtr> &data, size_t timeStep, double stddev);
+			void correctDataUpTo(std::vector<Image2DPtr> &data, size_t &dataTimeIndex, double time, const std::vector<double> &observationTimes, double stddev);
+	};
+}
+
+#endif

CEP/DP3/AOFlagger/src/CMakeLists.txt

@@ -116,6 +116,7 @@ set(STRATEGY_ACTION_FILES
   strategy/actions/frequencyselectionaction.cpp
   strategy/actions/fringestopaction.cpp
   strategy/actions/imageraction.cpp
+  strategy/actions/normalizevarianceaction.cpp
   strategy/actions/plotaction.cpp
   strategy/actions/slidingwindowfitaction.cpp
   strategy/actions/spatialcompositionaction.cpp

CEP/DP3/AOFlagger/src/strategy/actions/normalizevarianceaction.cpp

+/***************************************************************************
+ *   Copyright (C) 2008 by A.R. Offringa   *
+ *   offringa@astro.rug.nl   *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ *   This program is distributed in the hope that it will be useful,       *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU General Public License for more details.                          *
+ *                                                                         *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+#include <AOFlagger/strategy/actions/normalizevarianceaction.h>
+
+#include <map>
+
+#include <AOFlagger/util/progresslistener.h>
+
+#include <AOFlagger/quality/defaultstatistics.h>
+#include <AOFlagger/quality/qualitytablesformatter.h>
+#include <AOFlagger/quality/statisticalvalue.h>
+#include <AOFlagger/quality/statisticscollection.h>
+
+#include <AOFlagger/strategy/algorithms/medianwindow.h>
+#include <AOFlagger/strategy/imagesets/msimageset.h>
+#include <AOFlagger/quality/statisticsderivator.h>
+
+namespace rfiStrategy {
+	
+void NormalizeVarianceAction::Perform(ArtifactSet &artifacts, ProgressListener &progress)
+{
+	if(!artifacts.HasImageSet())
+		throw std::runtime_error("Normalize variance called without image set");
+	ImageSet *imageSet = artifacts.ImageSet();
+	MSImageSet *msImageSet = dynamic_cast<MSImageSet*>(imageSet);
+	if(msImageSet == 0)
+		throw std::runtime_error("Normalize variance actions needs measurement set");
+	std::string filename = msImageSet->Reader()->Set().Location();
+	QualityTablesFormatter qtables(filename);
+	StatisticsCollection statCollection(msImageSet->Reader()->Set().GetPolarizationCount());
+	statCollection.LoadTimeStatisticsOnly(qtables);
+	statCollection.IntegrateTimeToOneChannel();
+	
+	TimeFrequencyData &original = artifacts.OriginalData();
+	const std::vector<double> &observationTimes = artifacts.MetaData()->ObservationTimes();
+	size_t
+		width = original.ImageWidth();
+	
+	std::vector<Image2DPtr> data;
+	for(unsigned img=0;img<original.ImageCount();++img)
+		data.push_back(Image2D::CreateCopy(original.GetImage(img)));
+		
+	// Calculate all stddevs
+	const std::map<double, DefaultStatistics> &statMap = statCollection.TimeStatistics();
+	std::map<double, double> stddevs;
+	for(std::map<double, DefaultStatistics>::const_iterator i = statMap.begin();
+			i != statMap.end(); ++i)
+	{
+		double stddev = StatisticsDerivator::GetStatisticAmplitude(
+			QualityTablesFormatter::DStandardDeviationStatistic,
+			i->second.ToSinglePolarization(), 0);
+		stddevs.insert(std::pair<double, double>(i->first, stddev));
+	}
+	
+	// Add the first half of the window
+	const double halfWindowTime = _medianFilterSizeInS * 0.5;
+	MedianWindow<double> window;
+	std::map<double, double>::const_iterator windowRightSideIterator = stddevs.begin();
+	const double startTime = windowRightSideIterator->first;
+	do {
+		window.Add(windowRightSideIterator->second);
+		++windowRightSideIterator;
+	} while(windowRightSideIterator->first - startTime < halfWindowTime);
+	
+	// Add the second half, and start correcting the data
+	size_t dataTimeIndex = 0;
+	while(windowRightSideIterator != stddevs.end() &&
+		windowRightSideIterator->first - startTime < _medianFilterSizeInS)
+	{
+		correctDataUpTo(data, dataTimeIndex, windowRightSideIterator->first, observationTimes, window.Median());
+		window.Add(windowRightSideIterator->second);
+		++windowRightSideIterator;
+	}
+	
+	// Slide window until right side hits end
+	std::map<double, double>::const_iterator windowLeftSideIterator = stddevs.begin();
+	const double endTime = stddevs.rbegin()->first;
+	while(windowRightSideIterator != stddevs.end() && windowRightSideIterator->first < endTime)
+	{
+		correctDataUpTo(data, dataTimeIndex, windowRightSideIterator->first, observationTimes, window.Median());
+		
+		window.Add(windowRightSideIterator->second);
+		window.Remove(windowLeftSideIterator->second);
+		
+		++windowRightSideIterator;
+		++windowLeftSideIterator;
+	}
+	
+	// Slide until window center hits end
+	while(windowRightSideIterator->first + halfWindowTime < endTime)
+	{
+		correctDataUpTo(data, dataTimeIndex, windowLeftSideIterator->first + _medianFilterSizeInS, observationTimes, window.Median());
+		window.Remove(windowLeftSideIterator->second);
+		++windowLeftSideIterator;
+	}
+	
+	while(dataTimeIndex < width)
+	{
+		correctData(data, dataTimeIndex, window.Median());
+		++dataTimeIndex;
+	}
+	
+	// Replace images
+	for(unsigned img=0;img<original.ImageCount();++img)
+		original.SetImage(img, data[img]);
+}
+
+void NormalizeVarianceAction::correctDataUpTo(std::vector<Image2DPtr> &data, size_t &dataTimeIndex, double rightSideTime, const std::vector<double> &observationTimes, double stddev)
+{
+	size_t width = (*data.begin())->Width();
+	double halfWindowWidth = _medianFilterSizeInS*0.5;
+	while(dataTimeIndex < width &&
+		observationTimes[dataTimeIndex] + halfWindowWidth < rightSideTime)
+	{
+		correctData(data, dataTimeIndex, stddev);
+		++dataTimeIndex;
+	}
+}
+
+void NormalizeVarianceAction::correctData(std::vector<Image2DPtr> &data, size_t timeStep, double stddev)
+{
+	AOLogger::Debug << "Scaling timestep " << timeStep << " with " << stddev << '\n';
+	num_t oneOverStddev = 1.0 / stddev;
+	
+	for(std::vector<Image2DPtr>::iterator i=data.begin();i!=data.end();++i)
+	{
+		Image2DPtr image = *i;
+		
+		for(unsigned y=0;y<image->Height();++y)
+			image->SetValue(timeStep, y, image->Value(timeStep, y) * oneOverStddev);
+	}
+}
+
+} // end of namespace

CEP/DP3/AOFlagger/src/strategy/control/actionfactory.cpp

@@ -42,6 +42,7 @@
 #include <AOFlagger/strategy/actions/fringestopaction.h>
 #include <AOFlagger/strategy/actions/imageraction.h>
 #include <AOFlagger/strategy/actions/iterationaction.h>
+#include <AOFlagger/strategy/actions/normalizevarianceaction.h>
 #include <AOFlagger/strategy/actions/plotaction.h>
 #include <AOFlagger/strategy/actions/quickcalibrateaction.h>
 #include <AOFlagger/strategy/actions/rawappenderaction.h>
@@ -85,6 +86,7 @@ const std::vector<std::string> ActionFactory::GetActionList()
 	list.push_back("Fringe stopping recovery");
 	list.push_back("Image");
 	list.push_back("Iteration");
+	list.push_back("Normalize variance");
 	list.push_back("Phase adapter");
 	list.push_back("Plot");
 	list.push_back("Quickly calibrate");
@@ -149,6 +151,8 @@ Action *ActionFactory::CreateAction(const std::string &action)
 		return new ImagerAction();
 	else if(action == "Iteration")
 		return new IterationBlock();
+	else if(action == "Normalize variance")
+		return new NormalizeVarianceAction();
 	else if(action == "Phase adapter")
 		return new Adapter();
 	else if(action == "Plot")