compile ok mais segmentations

Hacene SI HADJ MOHAND
1 parent 5adfc136
Showing 1 changed file with 185 additions and 52 deletions Show diff stats
src/ExternLib/DataFiltering/DataFiltering.hh
@@ -16,6 +16,7 @@
 #include <list>
 #include <c++/4.8.2/bits/stl_vector.h>
+#include <algorithm>
 #include "Parameter.hh"
 #include "ParamData.hh"
@@ -26,10 +27,15 @@
 namespace AMDA {
     namespace Parameters {
+        /**
+         * Spike filtering using 
+         * An algorithm for finding spurious points in turbulent signals, Computers in Physics, V. 7 No. 5, pp. 599 - 607
+         */
         namespace DataFiltering {
             template <typename ElemType>
             struct FilteringContainer {
+                ~FilteringContainer(){}
                 // size of real points(not nan and not spike)  
                 int _size = 0;
                 int _nFilteredPoints = 0;
@@ -43,6 +49,10 @@ namespace AMDA {
             template <typename ElemType>
             struct FilteringContainer1D {
+                FilteringContainer1D(){}
+                ~FilteringContainer1D(){
+                    delete &_size;
+                }
                 // size of real points(not nan and not spike)  
                 std::vector<int> _size;
                 std::vector<int> _nFilteredPoints;
@@ -50,7 +60,7 @@ namespace AMDA {
                 std::vector<ElemType> _sumsq;
                 std::vector<ElemType> _ave;
                 std::vector<ElemType> _sig;
-                std::vector<std::list<ElemType> > _values;
+                std::list<std::vector<ElemType> > _values;
                 std::list<double> _times;
             };
@@ -80,7 +90,7 @@ namespace AMDA {
                         double crtTime = _paramInput.getTime(_index);
                         ElemType inputElt = _paramInput.get(_index);
-                        // initialiser et remplir container
+                        // initialiser et remplir le container container
                         if (_container._size < _nPoints) {
                             _container._values.push_back(inputElt);
@@ -92,72 +102,90 @@ namespace AMDA {
                                 _container._size += 1;
                             }
                         } else {
+                            // traitement une fois le container rempli 
                             _container._ave = (ElemType) _container._sum / _container._size;
                             _container._sig = std::sqrt((_container._sumsq - _container._sum * _container._sum / _container._size) / _container._size);
                             // filter Elements
                             for (auto it = _container._values.begin(); it != _container._values.end(); ++it) {
                                 ElemType crt_val = *it;
-                                if(isNAN(crt_val))
+                                if (isNAN(crt_val))
                                     continue;
                                 if (std::abs(crt_val - _container._ave) > _factor * _container._sig) {
                                     _container._sum -= crt_val;
                                     _container._sumsq -= crt_val*crt_val;
                                     _container._size -= 1;
                                     _container._nFilteredPoints += 1;
-                                     *it << NotANumber();
+                                    *it << NotANumber();
                                 }
                             }
                             _paramOutput->pushTime(_container._times.front());
                             _paramOutput->getDataList().push_back(_container._values.front());
                             _container._values.pop_front();
                             _container._times.pop_front();
-                            if(!isNAN(_container._values.front())){
-                            _container._size -= 1;
-                            _container._sum -= _container._values.front();
-                            _container._sumsq -= _container._values.front() * _container._values.front();
+                            if (!isNAN(_container._values.front())) {
+                                _container._size -= 1;
+                                _container._sum -= _container._values.front();
+                                _container._sumsq -= _container._values.front() * _container._values.front();
                             }
-                            
+
                             _container._values.push_back(inputElt);
                             _container._times.push_back(crtTime);
-                             if (!isNAN(inputElt)) {
+                            if (!isNAN(inputElt)) {
                                 _container._sum += inputElt;
                                 _container._sumsq += inputElt*inputElt;
                                 _container._size += 1;
                             }
-                        }   
-                            // last value filter and write
-                            if (_index == pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess - 1) {
+                        }
+                        // last value filter and write
+                        if (_index == pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess - 1) {
                             _container._ave = (ElemType) _container._sum / _container._size;
                             _container._sig = std::sqrt((_container._sumsq - _container._sum * _container._sum / _container._size) / _container._size);
-                              for (auto it = _container._values.begin(); it != _container._values.end(); ++it) {
+                            for (auto it = _container._values.begin(); it != _container._values.end(); ++it) {
                                 ElemType crt_val = *it;
-                                if (std::abs(crt_val - _container._ave) > _factor * _container._sig){ 
+                                if (std::abs(crt_val - _container._ave) > _factor * _container._sig) {
                                     *it << NotANumber();
                                     _container._nFilteredPoints += 1;
                                 }
-                              }
-                                
-                                std::list<double>::iterator itTimes = _container._times.begin();
-                                for (ElemType elt : _container._values) {
-                                    _paramOutput->pushTime(*itTimes);
-                                    _paramOutput->getDataList().push_back(elt);
-                                    ++itTimes;
-                                }
                             }
-                        
+
+                            std::list<double>::iterator itTimes = _container._times.begin();
+                            for (ElemType elt : _container._values) {
+                                _paramOutput->pushTime(*itTimes);
+                                _paramOutput->getDataList().push_back(elt);
+                                ++itTimes;
+                            }
+                        }
+
                     }
                 }
             private:
-                ParamDataSpec<ElemType>& _paramInput;
-
-                ParamDataSpec<ElemType>* _paramOutput;
-
+                /**
+                 * input data
+                 */
+                ParamDataSpec<ElemType >& _paramInput;
+
+                /**
+                 * oupout filtered data 
+                 */
+                ParamDataSpec<ElemType >* _paramOutput;
+
+                /**
+                 * factor in term of sigma used in the filtering   
+                 */
                 double _factor;
+
+                /**
+                 Number of point used in the fitering algorithm
+                 */
                 int _nPoints;
+
+                /**
+                 * Data container used during filtering 
+                 */
                 FilteringContainer<ElemType> _container;
             };
@@ -172,50 +200,155 @@ namespace AMDA {
                 _paramInput(paramInput),
                 _paramOutput(new ParamDataSpec<std::vector<ElemType> >()),
                 _factor(factor),
-                _nPoints(nPoints) {
+                _nPoints(nPoints),
+                _container(new FilteringContainer1D<ElemType> ()){
                     _paramDataOutput = _paramOutput;
                 }
                 virtual ~DataFiltering1D() {
-
                 }
                 void write(ParamDataIndexInfo &pParamDataIndexInfo) {
+                    
                     // init container 
-                    std::vector<std::list<ElemType> > outputList;
+                   
                     int size = _paramInput.get(pParamDataIndexInfo._startIndex).size();
-                    std::vector<int> vecInt(size, 0);
-                    std::vector<ElemType> vecElement(size, 0);
-                    _container._size = vecInt;
-                    _container._nFilteredPoints = vecInt;
-                    _container._ave = vecElement;
-                    _container._sum = vecElement;
-                    _container._sumsq = vecElement;
-                    _container._sig = vecElement;
-                    for (int j = 0; j < size; ++j) {
-                        std::list<ElemType> list;
-                        _container._values.push_back(list);
-                        outputList.push_back(list);
-                    }
-                    for (int i = 0; i < size; ++i) {
+                    std::vector<int> vectInt(size,0);
+                    _container->_size = vectInt;
+                    /**
+                    _container->_nFilteredPoints.resize(size,0);
+                    _container->_ave.resize(size,0);
+                    _container->_sum.resize(size,0);
+                    _container->_sumsq.resize(size,0);
+                    _container->_sig.resize(size,0);
+                 */
+                    for (unsigned int _index = pParamDataIndexInfo._startIndex;
+                            _index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess;
+                            ++_index) {
+                        double crtTime = _paramInput.getTime(_index);
+                        std::vector<ElemType> inputElt = _paramInput.get(_index);
-                    }
-                    
-                    
-                    
-                
-                }
+                         _paramOutput->pushTime(crtTime);
+                         _paramOutput->getDataList().push_back(inputElt);
+                                 /**
+                        // initialiser et remplir container
+                        if (*std::min_element(std::begin(_container._size), std::end(_container._size)) < _nPoints) {
+
+                            _container._values.push_back(inputElt);
+                            _container._times.push_back(crtTime);
+                            for (int i = 0; i < size; i++) {
+                                if (!isNAN(inputElt[i])) {
+                                    _container._sum[i] += inputElt[i];
+                                    _container._sumsq[i] += inputElt[i]*inputElt[i];
+                                    _container._size[i] += 1;
+                                }
+                            }
+                        } else {
+                            for (int i = 0; i < size; i++) {
+                                _container._ave[i] = (ElemType) _container._sum[i] / _container._size[i];
+                                _container._sig[i] = std::sqrt((_container._sumsq[i] - _container._sum[i] * _container._sum[i] / _container._size[i]) / _container._size[i]);
+                            }
+
+                            // filter Elements
+                            for (auto elt = _container._values.begin(); elt != _container._values.end(); ++elt) {
+                                int i;
+                                for (auto it = elt->begin(); it != elt->end(); ++it) {
+                                    i = it - elt->begin();
+                                    ElemType crt_val = *it;
+                                    if (isNAN(crt_val)) {
+                                        continue;
+                                    }
+
+                                    if (std::abs(crt_val - _container._ave[i]) > _factor * _container._sig[i]) {
+                                        _container._sum[i] -= crt_val;
+                                        _container._sumsq[i] -= crt_val*crt_val;
+                                        _container._size[i] -= 1;
+                                        _container._nFilteredPoints[i] += 1;
+                                        *it << NotANumber();
+                                    }
+                                }
+                            }
+                            _paramOutput->pushTime(_container._times.front());
+                            _paramOutput->getDataList().push_back(_container._values.front());
+                            _container._values.pop_front();
+                            _container._times.pop_front();
+                            for (int i = 0; i < size; i++) {
+                                if (!isNAN(_container._values.front()[i])) {
+                                    _container._size[i] -= 1;
+                                    _container._sum[i] -= _container._values.front()[i];
+                                    _container._sumsq[i] -= _container._values.front()[i] * _container._values.front()[i];
+                                }
+                                if (!isNAN(inputElt[i])) {
+                                    _container._sum[i] += inputElt[i];
+                                    _container._sumsq[i] += inputElt[i] * inputElt[i];
+                                    _container._size[i] += 1;
+                                }
+                            }
+                            _container._values.push_back(inputElt);
+                            _container._times.push_back(crtTime);
+                        }
+                        // last value filter and write
+                        if (_index == pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess - 1) {
+                               for (int i = 0; i < size; i++) {
+                            _container._ave[i] = (ElemType) _container._sum[i] / _container._size[i];
+                            _container._sig[i] = std::sqrt((_container._sumsq[i] - _container._sum[i] * _container._sum[i] / _container._size[i]) / _container._size[i]);
+                               }
+                            for (auto elt = _container._values.begin(); elt != _container._values.end(); ++elt) {
+                                int i;
+                                for (auto it = elt->begin(); it != elt->end(); ++it) {
+                                    i = it - elt->begin();
+                                    ElemType crt_val = *it;
+                                    if (isNAN(crt_val)) {
+                                        continue;
+                                    }
+
+                                    if (std::abs(crt_val - _container._ave[i]) > _factor * _container._sig[i]) {
+                                        _container._sum[i] -= crt_val;
+                                        _container._sumsq[i] -= crt_val*crt_val;
+                                        _container._size[i] -= 1;
+                                        _container._nFilteredPoints[i] += 1;
+                                        *it << NotANumber();
+                                    }
+                                }
+                            }
+                            std::list<double>::iterator itTimes = _container._times.begin();
+                            for (auto elt : _container._values) {
+                                _paramOutput->pushTime(*itTimes);
+                                _paramOutput->getDataList().push_back(elt);
+                                ++itTimes;
+                            }
+                        }
+                                  * */
+                    }
+                }
             private:
+                /**
+                 * input data
+                 */
                 ParamDataSpec<std::vector<ElemType> >& _paramInput;
+                /**
+                 * oupout filtered data 
+                 */
                 ParamDataSpec<std::vector<ElemType> >* _paramOutput;
+                /**
+                 * factor in term of sigma used in the filtering   
+                 */
                 double _factor;
+
+                /**
+                 Number of point used in the fitering algorithm
+                 */
                 int _nPoints;
-                FilteringContainer1D<ElemType> _container;
+
+                /**
+                 * Data container used during filtering 
+                 */
+                FilteringContainer1D<ElemType>* _container;
             };
         }