Implemented functionality of ChiSquaredFrequency class

Core/Algorithms/inc/ChiSquaredFrequency.h

@@ -24,8 +24,18 @@ public:
 
     virtual double calculateChiSquared(const OutputData<double> *p_simulation_data=0);
 
+    void setCutoff(double cutoff) {
+        if (cutoff>=0.0 && cutoff<=1.0) m_cutoff = cutoff;
+    }
+    double getCutoff() { return m_cutoff; }
+
     //! return output data which contains chi^2 values
     virtual OutputData<double > *createChi2DifferenceMap() const;
+protected:
+    virtual void initWeights();
+    OutputData<complex_t> *mp_real_ft;
+    OutputData<complex_t> *mp_simulation_ft;
+    double m_cutoff;
 };
 
 

Core/Algorithms/inc/IChiSquaredModule.h

@@ -52,7 +52,7 @@ public:
     virtual OutputData<double > *createChi2DifferenceMap() const=0;
 
 protected:
-    void initWeights();
+    virtual void initWeights();
     OutputData<double> *mp_real_data;
     OutputData<double> *mp_simulation_data;
     OutputData<double> *mp_weights;

Core/Algorithms/src/ChiSquaredFrequency.cpp

@@ -2,11 +2,16 @@
 
 ChiSquaredFrequency::ChiSquaredFrequency(const OutputData<double>& real_data)
 : IChiSquaredModule(real_data)
+, mp_real_ft(0)
+, mp_simulation_ft(0)
+, m_cutoff(1.0)
 {
 }
 
 ChiSquaredFrequency::~ChiSquaredFrequency()
 {
+    delete mp_real_ft;
+    delete mp_simulation_ft;
 }
 
 double ChiSquaredFrequency::calculateChiSquared(
@@ -15,12 +20,10 @@ double ChiSquaredFrequency::calculateChiSquared(
     if (p_simulation_data!=0) {
         setSimulationData(*p_simulation_data);
     }
-    if (mp_simulation_data==0) {
-        throw LogicErrorException("No simulation data present for calculating chi squared.");
-    }
+
     double result = 0.0;
-    size_t data_size = mp_real_data->getAllocatedSize();
     initWeights();
+    size_t data_size = mp_weights->getAllocatedSize();
     OutputData<double> *p_difference = createChi2DifferenceMap();
     mp_weights->resetIndex();
     p_difference->resetIndex();
@@ -33,17 +36,56 @@ double ChiSquaredFrequency::calculateChiSquared(
 
 OutputData<double>* ChiSquaredFrequency::createChi2DifferenceMap() const
 {
-    OutputData<double > *p_difference = mp_simulation_data->clone();
-    p_difference->setAllTo(0.0);
+    OutputData<double> *p_difference = mp_weights->clone();
 
-    mp_simulation_data->resetIndex();
-    mp_real_data->resetIndex();
     p_difference->resetIndex();
-    while (mp_real_data->hasNext()) {
-        double value_simu = mp_simulation_data->next();
-        double value_real = mp_real_data->next();
-        double squared_difference = mp_squared_function->calculateSquaredDifference(value_real, value_simu);
+    mp_simulation_ft->resetIndex();
+    mp_real_ft->resetIndex();
+
+    while (p_difference->hasNext()) {
+        complex_t diff = mp_simulation_ft->next() - mp_real_ft->next();
+        double squared_difference = std::norm(diff);
         p_difference->next() = squared_difference;
     }
 
-    return p_difference;}
+    return p_difference;
+}
+
+void ChiSquaredFrequency::initWeights()
+{
+    if (mp_simulation_data==0) {
+        throw LogicErrorException("No simulation data present for calculating chi squared.");
+    }
+    delete mp_real_ft;
+    delete mp_simulation_ft;
+    mp_real_ft = new OutputData<complex_t>();
+    mp_simulation_ft = new OutputData<complex_t>();
+    fourierTransform(*mp_real_data, mp_real_ft);
+    fourierTransform(*mp_simulation_data, mp_simulation_ft);
+    delete mp_weights;
+    mp_weights = new OutputData<double>();
+    size_t rank = mp_simulation_ft->getDimension();
+    int *n_dims = new int[rank];
+    for (size_t i=0; i<rank; ++i) {
+        n_dims[i] = mp_simulation_ft->getAxis(i)->getSize();
+    }
+    mp_weights->setAxisSizes(rank, n_dims);
+    delete n_dims;
+    mp_weights->resetIndex();
+    size_t nbr_rows = mp_weights->getAllSizes()[0];
+    size_t row_length = mp_weights->getAllSizes()[1];
+    size_t row_number = 0;
+    size_t counter = 0;
+    while (mp_weights->hasNext()) {
+        double weight = 0.0;
+        size_t column_index = counter%row_length;
+        int shift = (column_index>=row_length/2 ? -(int)row_length : 0);
+        int centered_column_index = (int)column_index - shift;
+        if (row_number<m_cutoff*nbr_rows && centered_column_index < m_cutoff*row_length/2.0) {
+            weight = 1.0;
+        }
+        mp_weights->next() = weight;
+        ++counter;
+        if (counter%row_length==0) ++row_number;
+    }
+}