From 206a9fc110835c1ce11af29d164926df86522188 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9=20Offringa?= <offringa@astron.nl>
Date: Sat, 17 Mar 2012 22:40:20 +0000
Subject: [PATCH] Task #1892: added variance normalization to strategy files,
made the action thread safe
---
.../actions/normalizevarianceaction.h | 23 ++++-
.../strategy/control/strategyreader.h | 1 +
.../strategy/control/strategywriter.h | 1 +
.../AOFlagger/strategy/control/types.h | 3 +-
.../actions/normalizevarianceaction.cpp | 83 ++++++++++++-------
.../src/strategy/control/strategyreader.cpp | 20 +++--
.../src/strategy/control/strategywriter.cpp | 12 ++-
7 files changed, 103 insertions(+), 40 deletions(-)
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/normalizevarianceaction.h b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/normalizevarianceaction.h
index 25c8f1bed93..63976ae26a6 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/normalizevarianceaction.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/normalizevarianceaction.h
@@ -22,6 +22,8 @@
#include <map>
+#include <boost/thread/mutex.hpp>
+
#include <AOFlagger/strategy/actions/action.h>
#include <AOFlagger/strategy/control/artifactset.h>
@@ -34,23 +36,36 @@ namespace rfiStrategy {
class NormalizeVarianceAction : public Action {
public:
NormalizeVarianceAction() :
- _medianFilterSizeInS(60.0*15.0)
- {
- }
- virtual ~NormalizeVarianceAction()
+ _medianFilterSizeInS(60.0*15.0), _isInitialized(false)
{
}
+ virtual ~NormalizeVarianceAction();
virtual std::string Description()
{
return "Normalize variance over time";
}
virtual void Perform(ArtifactSet &artifacts, ProgressListener &progress);
+ virtual void Initialize()
+ {
+ clean();
+ }
+ virtual void Finish()
+ {
+ clean();
+ }
virtual ActionType Type() const { return NormalizeVarianceActionType; }
+ double MedianFilterSizeInS() const { return _medianFilterSizeInS; }
+ void SetMedianFilterSizeInS(double filterSize) { _medianFilterSizeInS = filterSize; }
private:
double _medianFilterSizeInS;
+ boost::mutex _mutex;
+ bool _isInitialized;
+ std::map<double, double> _stddevs;
+ void initializeStdDevs(ArtifactSet &artifacts);
+ void clean();
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);
};
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/control/strategyreader.h b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/control/strategyreader.h
index 57165808e52..f40fcdb72ed 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/control/strategyreader.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/control/strategyreader.h
@@ -77,6 +77,7 @@ class StrategyReader {
class Action *parseFringeStopAction(xmlNode *node);
class Action *parseImagerAction(xmlNode *node);
class Action *parseIterationBlock(xmlNode *node);
+ class Action *parseNormalizeVarianceAction(xmlNode *node);
class Action *parsePlotAction(xmlNode *node);
class Action *parseQuickCalibrateAction(xmlNode *node);
class Action *parseRawAppenderAction(xmlNode *node);
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/control/strategywriter.h b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/control/strategywriter.h
index 2c09c84dbf8..7586f279643 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/control/strategywriter.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/control/strategywriter.h
@@ -92,6 +92,7 @@ namespace rfiStrategy {
void writeFringeStopAction(const class FringeStopAction &action);
void writeImagerAction(const class ImagerAction &action);
void writeIterationBlock(const class IterationBlock &action);
+ void writeNormalizeVarianceAction(const class NormalizeVarianceAction &action);
void writePlotAction(const class PlotAction &action);
void writeQuickCalibrateAction(const class QuickCalibrateAction &action);
void writeRawAppenderAction(const class RawAppenderAction &action);
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/control/types.h b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/control/types.h
index 7e3747c54ab..3ed997f4807 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/control/types.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/control/types.h
@@ -83,7 +83,8 @@ class Observatorium;
// 3.4 : Renamed ThresholdAction to SumThresholdAction
// 3.5 : Added the AbsThresholdAction
// 3.6 : Added the DirectionProfileAction and the EigenValueVerticalAction.
-#define STRATEGY_FILE_FORMAT_VERSION 3.6
+// 3.7 : Added the NormalizeVarianceAction
+#define STRATEGY_FILE_FORMAT_VERSION 3.7
// The earliest format version which can be read by this version of the software
#define STRATEGY_FILE_FORMAT_VERSION_REQUIRED 3.4
diff --git a/CEP/DP3/AOFlagger/src/strategy/actions/normalizevarianceaction.cpp b/CEP/DP3/AOFlagger/src/strategy/actions/normalizevarianceaction.cpp
index 06d4b38c8b3..022346b4e6e 100644
--- a/CEP/DP3/AOFlagger/src/strategy/actions/normalizevarianceaction.cpp
+++ b/CEP/DP3/AOFlagger/src/strategy/actions/normalizevarianceaction.cpp
@@ -33,20 +33,58 @@
#include <AOFlagger/quality/statisticsderivator.h>
namespace rfiStrategy {
+
+NormalizeVarianceAction::~NormalizeVarianceAction()
+{
+}
+
+void NormalizeVarianceAction::initializeStdDevs(ArtifactSet &artifacts)
+{
+ // The thread that calls this function first will initialize the
+ // std dev. When a new measurement set is read, Initialize or Finalize
+ // will be called, causing a clean().
+ boost::mutex::scoped_lock lock(_mutex);
+ if(!_isInitialized)
+ {
+ 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();
+ _isInitialized = true;
+
+ // Calculate all stddevs
+ const std::map<double, DefaultStatistics> &statMap = statCollection.TimeStatistics();
+ _stddevs.clear();
+ std::map<double, double>::iterator pos = _stddevs.begin();
+ for(std::map<double, DefaultStatistics>::const_iterator i = statMap.begin();
+ i != statMap.end(); ++i)
+ {
+ double stddev = StatisticsDerivator::GetStatisticAmplitude(
+ QualityTablesFormatter::DStandardDeviationStatistic,
+ i->second.ToSinglePolarization(), 0);
+ pos = _stddevs.insert(pos, std::pair<double, double>(i->first, stddev));
+ }
+ }
+}
+
+void NormalizeVarianceAction::clean()
+{
+ boost::mutex::scoped_lock lock(_mutex);
+ _isInitialized = false;
+ _stddevs.clear(); // frees a bit of memory.
+}
+
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();
+ initializeStdDevs(artifacts);
TimeFrequencyData &original = artifacts.OriginalData();
const std::vector<double> &observationTimes = artifacts.MetaData()->ObservationTimes();
@@ -57,22 +95,10 @@ void NormalizeVarianceAction::Perform(ArtifactSet &artifacts, ProgressListener &
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();
+ std::map<double, double>::const_iterator windowRightSideIterator = _stddevs.begin();
const double startTime = windowRightSideIterator->first;
do {
window.Add(windowRightSideIterator->second);
@@ -81,7 +107,7 @@ void NormalizeVarianceAction::Perform(ArtifactSet &artifacts, ProgressListener &
// Add the second half, and start correcting the data
size_t dataTimeIndex = 0;
- while(windowRightSideIterator != stddevs.end() &&
+ while(windowRightSideIterator != _stddevs.end() &&
windowRightSideIterator->first - startTime < _medianFilterSizeInS)
{
correctDataUpTo(data, dataTimeIndex, windowRightSideIterator->first, observationTimes, window.Median());
@@ -90,9 +116,9 @@ void NormalizeVarianceAction::Perform(ArtifactSet &artifacts, ProgressListener &
}
// 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)
+ 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());
@@ -136,7 +162,6 @@ void NormalizeVarianceAction::correctDataUpTo(std::vector<Image2DPtr> &data, siz
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)
diff --git a/CEP/DP3/AOFlagger/src/strategy/control/strategyreader.cpp b/CEP/DP3/AOFlagger/src/strategy/control/strategyreader.cpp
index ca82c6f80f4..cf4441cf48f 100644
--- a/CEP/DP3/AOFlagger/src/strategy/control/strategyreader.cpp
+++ b/CEP/DP3/AOFlagger/src/strategy/control/strategyreader.cpp
@@ -43,6 +43,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>
@@ -280,6 +281,8 @@ Action *StrategyReader::parseAction(xmlNode *node)
newAction = parseImagerAction(node);
else if(typeStr == "IterationBlock")
newAction = parseIterationBlock(node);
+ else if(typeStr == "NormalizeVarianceAction")
+ newAction = parseNormalizeVarianceAction(node);
else if(typeStr == "PlotAction")
newAction = parsePlotAction(node);
else if(typeStr == "QuickCalibrateAction")
@@ -533,7 +536,7 @@ Action *StrategyReader::parseFourierTransformAction(xmlNode *)
return newAction;
}
-class Action *StrategyReader::parseFrequencyConvolutionAction(xmlNode *node)
+Action *StrategyReader::parseFrequencyConvolutionAction(xmlNode *node)
{
FrequencyConvolutionAction *newAction = new FrequencyConvolutionAction();
newAction->SetConvolutionSize(getInt(node, "convolution-size"));
@@ -541,14 +544,14 @@ class Action *StrategyReader::parseFrequencyConvolutionAction(xmlNode *node)
return newAction;
}
-class Action *StrategyReader::parseFrequencySelectionAction(xmlNode *node)
+Action *StrategyReader::parseFrequencySelectionAction(xmlNode *node)
{
FrequencySelectionAction *newAction = new FrequencySelectionAction();
newAction->SetThreshold(getDouble(node, "threshold"));
return newAction;
}
-class Action *StrategyReader::parseFringeStopAction(xmlNode *node)
+Action *StrategyReader::parseFringeStopAction(xmlNode *node)
{
FringeStopAction *newAction = new FringeStopAction();
newAction->SetFitChannelsIndividually(getBool(node, "fit-channels-individually"));
@@ -559,13 +562,13 @@ class Action *StrategyReader::parseFringeStopAction(xmlNode *node)
return newAction;
}
-class Action *StrategyReader::parseImagerAction(xmlNode *)
+Action *StrategyReader::parseImagerAction(xmlNode *)
{
ImagerAction *newAction = new ImagerAction();
return newAction;
}
-class Action *StrategyReader::parseIterationBlock(xmlNode *node)
+Action *StrategyReader::parseIterationBlock(xmlNode *node)
{
IterationBlock *newAction = new IterationBlock();
newAction->SetIterationCount(getInt(node, "iteration-count"));
@@ -574,6 +577,13 @@ class Action *StrategyReader::parseIterationBlock(xmlNode *node)
return newAction;
}
+Action *StrategyReader::parseNormalizeVarianceAction(xmlNode *node)
+{
+ NormalizeVarianceAction *newAction = new NormalizeVarianceAction();
+ newAction->SetMedianFilterSizeInS(getDouble(node, "median-filter-size-in-s"));
+ return newAction;
+}
+
Action *StrategyReader::parseSetFlaggingAction(xmlNode *node)
{
SetFlaggingAction *newAction = new SetFlaggingAction();
diff --git a/CEP/DP3/AOFlagger/src/strategy/control/strategywriter.cpp b/CEP/DP3/AOFlagger/src/strategy/control/strategywriter.cpp
index 91db62daf33..3405d3cc65d 100644
--- a/CEP/DP3/AOFlagger/src/strategy/control/strategywriter.cpp
+++ b/CEP/DP3/AOFlagger/src/strategy/control/strategywriter.cpp
@@ -40,6 +40,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>
@@ -161,6 +162,9 @@ namespace rfiStrategy {
case IterationBlockType:
writeIterationBlock(static_cast<const IterationBlock&>(action));
break;
+ case NormalizeVarianceActionType:
+ writeNormalizeVarianceAction(static_cast<const NormalizeVarianceAction&>(action));
+ break;
case PlotActionType:
writePlotAction(static_cast<const PlotAction&>(action));
break;
@@ -399,7 +403,13 @@ namespace rfiStrategy {
writeContainerItems(action);
}
- void StrategyWriter::writePlotAction(const class PlotAction &action)
+ void StrategyWriter::writeNormalizeVarianceAction(const NormalizeVarianceAction &action)
+ {
+ Attribute("type", "NormalizeVarianceAction");
+ Write<double>("median-filter-size-in-s", action.MedianFilterSizeInS());
+ }
+
+ void StrategyWriter::writePlotAction(const PlotAction &action)
{
Attribute("type", "PlotAction");
Write<int>("plot-kind", action.PlotKind());
--
GitLab