CDPP / AMDA_Kernel

Commit e336270ad14534c6d0e3e782b8ca3fe7c84beb27

Authored by Benjamin Renard
1 parent c2168efe
Exists in master and in 41 other branches AllLegends, CDF_ISTP, CDF_ISTP_HOPEFULLY_ALL_WORKING, FER_10400, FER_10412, FER_10765, FER_11000, FER_11100, FER_11242, FER_11277, FER_11500, FER_11537, FER_11588, FER_11651, FER_19836, FER_5344, FER_6031, FER_6966, FER_7318, FER_7618, FER_7857, FER_8455, FER_9629, FER_CR41, FER_CR42, FER_CR43, FER_TFCat, Histo1D2D, JUNO_JEDI, MAZ_10936, MAZ_6018, MAZ_9016_New, MAZ_FIX_9016, SpeasyGet, amdadev, cu2022, das2, density_maps, develop, evol_9042, ubuntu

Fix correlation & covariance

Showing 6 changed files with 306 additions and 422 deletions   Show diff stats
Inline Side-by-side
src/ExternLib/StatisticFunctions/AbstractFunc.hh
  Diff comments   View file @e336270
... ... @@ -204,7 +204,86 @@ private:
204 204  
205 205 protected:
206 206 OutputElemType _nanVal;
207   -};
  207 +};
  208 +
  209 +template <typename InputElemType, typename OutputElemType>
  210 +class Abstract2ParamsFunc : public AbstractFuncBase {
  211 +public:
  212 + /**
  213 + * @brief Constructor.
  214 + * @details Create the ParamData type of the input ParamData.
  215 + */
  216 + Abstract2ParamsFunc(Process& pProcess, TimeIntervalListSPtr pTimeIntervalList, ParamDataSpec<InputElemType>& param1Input, ParamDataSpec<InputElemType>& param2Input, double windowtime)
  217 + : AbstractFuncBase(pProcess, pTimeIntervalList, windowtime),
  218 + _param1Input(param1Input),
  219 + _param2Input(param2Input),
  220 + _paramOutput(new ParamDataSpec<OutputElemType>) {
  221 + _paramDataOutput=_paramOutput;
  222 + }
  223 +
  224 + virtual ~Abstract2ParamsFunc() {
  225 + }
  226 +
  227 + virtual void pushData(double time, InputElemType& elem1, InputElemType& elem2) = 0;
  228 +
  229 + virtual OutputElemType compute() = 0;
  230 +
  231 + /**
  232 + * @overload Operation::write(ParamDataIndexInfo &pParamDataIndexInfo)
  233 + */
  234 +
  235 + void write(ParamDataIndexInfo &pParamDataIndexInfo) {
  236 + if ((pParamDataIndexInfo._nbDataToProcess > 0)) {
  237 + for (unsigned int _index = pParamDataIndexInfo._startIndex ;
  238 + _index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess;
  239 + ++_index)
  240 + {
  241 + double crtTime = _param1Input.getTime(_index);
  242 + InputElemType crtVal1 = _param1Input.get(_index);
  243 + InputElemType crtVal2 = _param2Input.get(_index);
  244 + if (needToChangeTarget(crtTime)) {
  245 + _paramOutput->pushTime(getTarget());
  246 + _paramOutput->push(compute());
  247 + pushData(crtTime, crtVal1, crtVal2);
  248 + nextTarget();
  249 + bool skip = false;
  250 + while (!skip && needToChangeTarget(crtTime)) {
  251 + _paramOutput->pushTime(getTarget());
  252 + _paramOutput->push(compute());
  253 + skip = nextTarget();
  254 + }
  255 + }
  256 + else {
  257 + pushData(crtTime, crtVal1, crtVal2);
  258 + if (needInit()) {
  259 + init();
  260 + }
  261 + }
  262 + }
  263 + }
  264 + if (pParamDataIndexInfo._timeIntToProcessChanged || pParamDataIndexInfo._noMoreTimeInt) {
  265 + if (!needInit()) {
  266 + do {
  267 + if (inInt(getTarget())) {
  268 + _paramOutput->pushTime(getTarget());
  269 + _paramOutput->push(compute());
  270 + }
  271 + } while (nextTarget());
  272 + }
  273 + }
  274 + }
  275 +
  276 + double getInputParamSampling() {
  277 + return _param1Input.getMinSampling();
  278 + }
  279 +
  280 +public:
  281 + ParamDataSpec<InputElemType>& _param1Input;
  282 +
  283 + ParamDataSpec<InputElemType>& _param2Input;
  284 +
  285 + ParamDataSpec<OutputElemType>* _paramOutput;
  286 +};
208 287  
209 288 template <typename InputElemType, typename OutputElemType>
210 289 class ClassicAbstractFunc : public AbstractFunc<InputElemType,OutputElemType> {
... ... @@ -293,15 +372,70 @@ public:
293 372 _targets.push_back(time);
294 373 }
295 374  
296   - virtual void resetFunc() {
297   - _mem.clear();
298   - _targets.clear();
299   - }
  375 + virtual void resetFunc() {
  376 + _mem.clear();
  377 + _targets.clear();
  378 + }
  379 +
  380 +protected:
  381 + std::list<double> _targets;
  382 +
  383 + std::list<std::pair<double,InputElemType> > _mem;
  384 +};
  385 +
  386 +template <typename InputElemType, typename OutputElemType>
  387 +class Sm2ParamsAbstractFunc : public Abstract2ParamsFunc<InputElemType, OutputElemType> {
  388 +public:
  389 + Sm2ParamsAbstractFunc(Process& pProcess, TimeIntervalListSPtr pTimeIntervalList, ParamDataSpec<InputElemType>& param1Input, ParamDataSpec<InputElemType>& param2Input, double windowtime)
  390 + : Abstract2ParamsFunc<InputElemType,OutputElemType>(pProcess, pTimeIntervalList, param1Input, param2Input, windowtime) {
  391 + }
  392 +
  393 + virtual ~Sm2ParamsAbstractFunc() {
  394 + }
  395 +
  396 + virtual void init() {
  397 + if (!_targets.empty()) {
  398 + Abstract2ParamsFunc<InputElemType,OutputElemType>::setTarget(_targets.front());
  399 + Abstract2ParamsFunc<InputElemType,OutputElemType>::setNeedInit(false);
  400 + _targets.pop_front();
  401 + }
  402 + }
  403 +
  404 + virtual bool nextTarget() {
  405 + if (!_targets.empty()) {
  406 + bool res = Abstract2ParamsFunc<InputElemType,OutputElemType>::setTarget(_targets.front());
  407 + _targets.pop_front();
  408 + while (!_mem.empty() && !Abstract2ParamsFunc<InputElemType,OutputElemType>::inWindow(_mem.front().first)) {
  409 + _mem.pop_front();
  410 + }
  411 + return res;
  412 + }
  413 + return false;
  414 + }
  415 +
  416 + virtual bool needToChangeTarget(double crtTime) {
  417 + return !Abstract2ParamsFunc<InputElemType,OutputElemType>::needInit() &&
  418 + !Abstract2ParamsFunc<InputElemType,OutputElemType>::inWindow(crtTime) && !_targets.empty() ;
  419 + }
  420 +
  421 + virtual double getSampling() {
  422 + return Abstract2ParamsFunc<InputElemType,OutputElemType>::getInputParamSampling();
  423 + }
  424 +
  425 + virtual void pushData(double time, InputElemType& elem1, InputElemType& elem2) {
  426 + _mem.push_back(std::make_pair(time, std::make_pair(elem1,elem2)));
  427 + _targets.push_back(time);
  428 + }
  429 +
  430 + virtual void resetFunc() {
  431 + _mem.clear();
  432 + _targets.clear();
  433 + }
300 434  
301 435 protected:
302 436 std::list<double> _targets;
303 437  
304   - std::list<std::pair<double,OutputElemType> > _mem;
  438 + std::list<std::pair<double,std::pair<InputElemType,InputElemType > > > _mem;
305 439 };
306 440  
307 441 } /* namespace StatisticFunctions */
... ...
src/ExternLib/StatisticFunctions/CorrelationFunctions.hh
  Diff comments   View file @e336270
... ... @@ -33,9 +33,6 @@ namespace AMDA {
33 33 namespace Parameters {
34 34 namespace StatisticFunctions {
35 35  
36   -#define AVERAGE_TIME 1200 // (seconds)
37   -#define MAX_GAP_SIZE 3600 // (seconds)
38   -
39 36 enum COEFS {
40 37 COVARIANCE = 1,
41 38 PAERSON = 2,
... ... @@ -53,150 +50,6 @@ namespace AMDA {
53 50 {"4", COEFS::KENDALL},
54 51 };
55 52  
56   - template <typename InputElemType, typename OutputElemType>
57   - class AbstractCorrelationFunc : public AbstractFuncBase {
58   - public:
59   -
60   - /**
61   - * @brief Constructor.
62   - * @details Create the ParamData type of the input ParamData.
63   - */
64   - AbstractCorrelationFunc(Process& pProcess, TimeIntervalListSPtr pTimeIntervalList, ParamDataSpec<InputElemType>& firstParamInput, ParamDataSpec<InputElemType>& secondParamInput, double windowtime, std::string correlationType)
65   - : AbstractFuncBase(pProcess, pTimeIntervalList, windowtime),
66   - _correlationType(correlationType),
67   - _firstParamInput(firstParamInput),
68   - _secondParamInput(secondParamInput),
69   - _paramOutput(new ParamDataSpec<OutputElemType>) {
70   - _paramDataOutput = _paramOutput;
71   - }
72   -
73   - virtual ~AbstractCorrelationFunc() {
74   - }
75   -
76   - virtual OutputElemType computeCorrelation(std::list<std::pair<double, std::pair<InputElemType, InputElemType>>>&mem, OutputElemType& nanVal, std::string type) = 0;
77   -
78   - virtual void init() {
79   - AbstractCorrelationFunc<InputElemType, OutputElemType>::setTarget(AbstractCorrelationFunc<InputElemType, OutputElemType>::getIntStartTime());
80   - AbstractCorrelationFunc<InputElemType, OutputElemType>::setNeedInit(false);
81   - }
82   -
83   - virtual bool nextTarget() {
84   - double target = AbstractCorrelationFunc<InputElemType, OutputElemType>::getTarget() + AbstractCorrelationFunc<InputElemType, OutputElemType>::getWindowTime();
85   - bool res = AbstractCorrelationFunc<InputElemType, OutputElemType>::setTarget(target);
86   - while (!_mem.empty() && !AbstractCorrelationFunc<InputElemType, OutputElemType>::inWindow(_mem.front().first)) {
87   - _mem.pop_front();
88   - }
89   - return res;
90   - }
91   -
92   - virtual bool needToChangeTarget(double crtTime) {
93   - return !AbstractCorrelationFunc<InputElemType, OutputElemType>::needInit() && !AbstractCorrelationFunc<InputElemType, OutputElemType>::inWindow(crtTime);
94   - }
95   -
96   - virtual double getSampling() {
97   - return AbstractCorrelationFunc<InputElemType, OutputElemType>::getWindowTime();
98   - }
99   -
100   - virtual void pushData(double time, std::pair<InputElemType, InputElemType>& elem) {
101   - _mem.push_back(std::make_pair(time, elem));
102   - }
103   -
104   - virtual void resetFunc() {
105   - _mem.clear();
106   - }
107   -
108   - void pushSecondParamData(ParamDataIndexInfo &pParamDataIndexInfo) {
109   - for (unsigned int _index = pParamDataIndexInfo._startIndex;
110   - _index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess;
111   - ++_index) {
112   - double time = _secondParamInput.getTime(_index);
113   - InputElemType val_ = _secondParamInput.get(_index);
114   - _secondParamInputData.push_back(std::pair<double, InputElemType> (time, val_));
115   - }
116   - }
117   -
118   - virtual InputElemType getValue(std::vector<std::pair<double, InputElemType> >& input, double time) = 0;
119   -
120   - /**
121   - * @overload Operation::write(ParamDataIndexInfo &pParamDataIndexInfo)
122   - */
123   -
124   - void write(ParamDataIndexInfo &pParamDataIndexInfo) {
125   -
126   - if ((pParamDataIndexInfo._nbDataToProcess > 0)) {
127   - if (pParamDataIndexInfo._startIndex == 0) {
128   - _nanVal = _firstParamInput.get(0);
129   - _nanVal << NotANumber();
130   - }
131   - for (unsigned int _index = pParamDataIndexInfo._startIndex;
132   - _index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess;
133   - ++_index) {
134   - double crtTime = _firstParamInput.getTime(_index);
135   - InputElemType firstVal = _firstParamInput.get(_index);
136   - // get the second element
137   - InputElemType secondVal = getValue(_secondParamInputData, crtTime);
138   - std::pair<InputElemType, InputElemType> crtVal(firstVal, secondVal);
139   -
140   - if (needToChangeTarget(crtTime)) {
141   - _paramOutput->pushTime(getTarget());
142   - _paramOutput->push(compute());
143   - pushData(crtTime, crtVal);
144   - nextTarget();
145   - bool skip = false;
146   - while (!skip && needToChangeTarget(crtTime)) {
147   - _paramOutput->pushTime(getTarget());
148   - _paramOutput->push(compute());
149   - skip = nextTarget();
150   - }
151   - } else {
152   - pushData(crtTime, crtVal);
153   - if (needInit()) {
154   - init();
155   - }
156   - }
157   - }
158   - }
159   - if (pParamDataIndexInfo._timeIntToProcessChanged || pParamDataIndexInfo._noMoreTimeInt) {
160   - if (!needInit()) {
161   - do {
162   - if (inInt(getTarget())) {
163   - _paramOutput->pushTime(getTarget());
164   - _paramOutput->push(compute());
165   - }
166   - } while (nextTarget());
167   - }
168   - }
169   -
170   - }
171   -
172   - double getInputParamSampling() {
173   - return _firstParamInput.getMinSampling();
174   - }
175   -
176   - OutputElemType compute() {
177   - return computeCorrelation(_mem, AbstractCorrelationFunc<InputElemType, OutputElemType>::_nanVal, _correlationType);
178   - }
179   -
180   -
181   - protected:
182   - OutputElemType _nanVal;
183   -
184   - std::string _correlationType;
185   -
186   - std::list<std::pair<double, std::pair<InputElemType, InputElemType>> > _mem;
187   -
188   - private:
189   - ParamDataSpec<InputElemType>& _firstParamInput;
190   -
191   - ParamDataSpec<InputElemType>& _secondParamInput;
192   -
193   - ParamDataSpec<OutputElemType>* _paramOutput;
194   -
195   - std::vector<std::pair<double, InputElemType> > _secondParamInputData;
196   -
197   -
198   - };
199   -
200 53 /**
201 54 *
202 55 * @param pProcess
... ... @@ -207,74 +60,40 @@ namespace AMDA {
207 60 * @param type
208 61 */
209 62 template <typename InputElemType, typename OutputElemType>
210   - class Correlation : public AbstractCorrelationFunc<InputElemType, OutputElemType> {
  63 + class Correlation : public Sm2ParamsAbstractFunc<InputElemType, OutputElemType> {
211 64 public:
212 65  
213 66 Correlation(Process & pProcess, TimeIntervalListSPtr pTimeIntervalList, ParamDataSpec<InputElemType>& firstParamInput,
214 67 ParamDataSpec<InputElemType>& secondParamInput, double windowtime, std::string correlationType) :
215   - AbstractCorrelationFunc<InputElemType, OutputElemType> (pProcess, pTimeIntervalList, firstParamInput, secondParamInput, windowtime, correlationType) {
  68 + Sm2ParamsAbstractFunc<InputElemType, OutputElemType> (pProcess, pTimeIntervalList, firstParamInput,
  69 + secondParamInput, windowtime), _correlationType(correlationType) {
216 70  
217 71 }
218 72  
219 73 virtual ~Correlation() {
220 74 }
221   -
222   - InputElemType getValue(std::vector<std::pair<double, InputElemType> >& input, double time) {
223   - double min_t = time - AVERAGE_TIME / 2.;
224   - double max_t = time + AVERAGE_TIME / 2.;
225   - std::vector<std::pair<double, InputElemType> > values_for_mean;
226   - InputElemType nanVal;
227   - nanVal << NotANumber();
228   - std::pair<double, InputElemType> prev_value(NAN, nanVal);
229   - std::pair<double, InputElemType> next_value(NAN, nanVal);
230   - InputElemType value = nanVal;
231   - for (auto it = input.begin(); it != input.end(); ++it) {
232   - if (it->first == time) {
233   - value = it->second;
234   - return value;
235   - break;
236   - } else if (isNAN(it->second))
237   - continue;
238   - else if (it->first > max_t) {
239   - next_value = *it;
240   - break;
241   - } else if (it->first < min_t) {
242   - prev_value = *it;
243   - } else {
244   - values_for_mean.push_back(*it);
245   - }
246   - }
247   - if (!values_for_mean.empty()) {
248   - //Compute mean
249   - InputElemType sum = 0;
250   - for (auto it = values_for_mean.begin(); it != values_for_mean.end(); ++it) {
251   - sum += it->second;
252   - }
253   - value = sum / (InputElemType) values_for_mean.size();
254   - } else {
255   - if (!isNAN(prev_value.first) && !isNAN(next_value.first) && (next_value.first - prev_value.first <= MAX_GAP_SIZE)) {
256   - //Compute interpolated value
257   - value = prev_value.second + (time - prev_value.first) / (next_value.first - prev_value.first) * (next_value.second - prev_value.second);
258   - }
259   - }
260   -
261   - return value;
262   - }
263 75  
264   - OutputElemType computeCorrelation(std::list<std::pair<double, std::pair<InputElemType, InputElemType>>>&mem, OutputElemType& nanVal, std::string type) {
265   - OutputElemType result = nanVal;
266   - if (mem.empty()) {
  76 + virtual OutputElemType compute() {
  77 + return computeCorrelation(Abstract2ParamsFunc<InputElemType, OutputElemType>::_paramOutput);
  78 + }
  79 +
  80 + double computeCorrelation(ParamDataSpec<double>* /*out*/) {
  81 + double result = 0.;
  82 + result << NotANumber();
  83 + if (Sm2ParamsAbstractFunc<InputElemType, OutputElemType>::_mem.empty()) {
267 84 return result;
268 85 }
269 86 std::list<std::pair<InputElemType, InputElemType>> list;
270   - for (typename std::list<std::pair<double, std::pair < InputElemType, InputElemType>>>::iterator it = mem.begin(); it != mem.end(); ++it) {
  87 + for (typename std::list<std::pair<double, std::pair < InputElemType, InputElemType>>>::iterator it = Sm2ParamsAbstractFunc<InputElemType, OutputElemType>::_mem.begin(); it != Sm2ParamsAbstractFunc<InputElemType, OutputElemType>::_mem.end(); ++it) {
  88 + if (isNAN(it->second.first) || isNAN(it->second.second))
  89 + continue;
271 90 list.push_back(it->second);
272 91 }
273   - if (coefsToStr.find(type) == coefsToStr.end()) {
274   - BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("StatisticFunctions::CorrelationFunction unknown correlation type " + type));
  92 + if (coefsToStr.find(_correlationType) == coefsToStr.end()) {
  93 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("StatisticFunctions::CorrelationFunction unknown correlation type " + _correlationType));
275 94 }
276 95  
277   - switch (coefsToStr[type]) {
  96 + switch (coefsToStr[_correlationType]) {
278 97 case COEFS::COVARIANCE:
279 98 getCovariance(list, result);
280 99 break;
... ... @@ -288,90 +107,30 @@ namespace AMDA {
288 107 getKendall(list, result);
289 108 break;
290 109 default:
291   - BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("StatisticFunctions::CorrelationFunction unknown correlation type :" + type));
  110 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("StatisticFunctions::CorrelationFunction unknown correlation type :" + _correlationType));
292 111 }
293 112  
294 113 return result;
295 114 }
296 115  
297   - };
298   -
299   - template <typename InputElemType, typename OutputElemType, typename dataType>
300   - class Correlation1D : public AbstractCorrelationFunc<InputElemType, OutputElemType> {
301   - public:
302   -
303   - Correlation1D(Process & pProcess, TimeIntervalListSPtr pTimeIntervalList, ParamDataSpec<InputElemType>& firstParamInput,
304   - ParamDataSpec<InputElemType>& secondParamInput, double windowtime, std::string correlationType) :
305   - AbstractCorrelationFunc<InputElemType, OutputElemType> (pProcess, pTimeIntervalList, firstParamInput, secondParamInput, windowtime, correlationType) {
306   -
307   - }
308   -
309   - virtual ~Correlation1D() {
310   - }
311   -
312   - InputElemType getValue(std::vector<std::pair<double, InputElemType> >& input, double time) {
313   - double min_t = time - AVERAGE_TIME / 2.;
314   - double max_t = time + AVERAGE_TIME / 2.;
315   - std::vector<std::pair<double, InputElemType> > values_for_mean;
316   - InputElemType nanVal;
317   - nanVal << NotANumber();
318   - std::pair<double, InputElemType> prev_value(NAN, nanVal);
319   - std::pair<double, InputElemType> next_value(NAN, nanVal);
320   - InputElemType value = nanVal;
321   - for (auto it = input.begin(); it != input.end(); ++it) {
322   - if (it->first == time) {
323   - value = it->second;
324   - return value;
325   - break;
326   - } else if (isNAN(it->second))
327   - continue;
328   - else if (it->first > max_t) {
329   - next_value = *it;
330   - break;
331   - } else if (it->first < min_t) {
332   - prev_value = *it;
333   - } else {
334   - values_for_mean.push_back(*it);
335   - }
336   - }
337   - if (!values_for_mean.empty()) {
338   - //Compute mean
339   - InputElemType sum = values_for_mean[0].second;
340   - for(int i= 0; i < sum.size();++i){
341   - for (int j=1 ; j <values_for_mean.size(); ++j ){
342   -
343   - sum[i] = sum[i] + values_for_mean[j].second[i];
344   - }
345   - value[i] = sum[i] / (float) values_for_mean.size();
346   - }
347   -
348   - } else {
349   - if (!isNAN(prev_value.first) && !isNAN(next_value.first) && (next_value.first - prev_value.first <= MAX_GAP_SIZE)) {
350   - //Compute interpolated value
351   - for(int i= 0; i < prev_value.second.size();++i)
352   - value[i] = prev_value.second[i] + (time - prev_value.first) / (next_value.first - prev_value.first) * (next_value.second[i] - prev_value.second[i]);
353   - }
354   - }
355   -
356   - return value;
357   - }
358   -
359   - OutputElemType computeCorrelation(std::list<std::pair<double, std::pair<InputElemType, InputElemType>>>&mem, OutputElemType& nanVal, std::string type) {
360   - OutputElemType result = nanVal;
361   - if (mem.empty()) {
  116 + std::vector<double> computeCorrelation(ParamDataSpec<std::vector<double> >* /*out*/) {
  117 + std::vector<double> result;
  118 + result.resize(result.size());
  119 + result << NotANumber();
  120 + if (Sm2ParamsAbstractFunc<InputElemType, OutputElemType>::_mem.empty()) {
362 121 return result;
363 122 }
364   - int n_ = mem.begin()->second.first.size();
  123 + int n_ = Sm2ParamsAbstractFunc<InputElemType, OutputElemType>::_mem.begin()->second.first.size();
365 124 for(int i = 0; i < n_; ++i){
366   - std::list<std::pair<dataType, dataType>> list;
367   - for (typename std::list<std::pair<double, std::pair < InputElemType, InputElemType>>>::iterator it = mem.begin(); it != mem.end(); ++it) {
368   - list.push_back(std::make_pair(it->second.first[i], it->second.second[i]));
  125 + std::list<std::pair<double, double>> list;
  126 + for (typename std::list<std::pair<double, std::pair < InputElemType, InputElemType>>>::iterator it = Sm2ParamsAbstractFunc<InputElemType, OutputElemType>::_mem.begin(); it != Sm2ParamsAbstractFunc<InputElemType, OutputElemType>::_mem.end(); ++it) {
  127 + list.push_back(std::make_pair((double)it->second.first[i], (double)it->second.second[i]));
369 128 }
370   - if (coefsToStr.find(type) == coefsToStr.end()) {
371   - BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("StatisticFunctions::CorrelationFunction unknown correlation type " + type));
  129 + if (coefsToStr.find(_correlationType) == coefsToStr.end()) {
  130 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("StatisticFunctions::CorrelationFunction unknown correlation type " + _correlationType));
372 131 }
373 132  
374   - switch (coefsToStr[type]) {
  133 + switch (coefsToStr[_correlationType]) {
375 134 case COEFS::COVARIANCE:
376 135 getCovariance(list, result[i]);
377 136 break;
... ... @@ -385,7 +144,7 @@ namespace AMDA {
385 144 getKendall(list, result[i]);
386 145 break;
387 146 default:
388   - BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("StatisticFunctions::CorrelationFunction unknown correlation type :" + type));
  147 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("StatisticFunctions::CorrelationFunction unknown correlation type :" + _correlationType));
389 148 }
390 149 list.clear();
391 150 }
... ... @@ -393,6 +152,9 @@ namespace AMDA {
393 152 return result;
394 153 }
395 154  
  155 + private:
  156 + std::string _correlationType;
  157 +
396 158 };
397 159  
398 160 }
... ...
src/ExternLib/StatisticFunctions/CorrelationProcess.cc
  Diff comments   View file @e336270
... ... @@ -34,41 +34,31 @@ using namespace log4cxx;
34 34 namespace AMDA {
35 35 namespace Parameters {
36 36  
37   - CorrelationProcess::CorrelationProcess(Parameter &parameter) : SingleParamProcess_CRTP(parameter) {
  37 + CorrelationProcess::CorrelationProcess(Parameter &parameter) : MultiParamProcess_CRTP(parameter) {
38 38 _type = "";
39 39 }
40 40  
41   - CorrelationProcess::CorrelationProcess(const CorrelationProcess& pProcess, Parameter &parameter) : SingleParamProcess_CRTP(pProcess, parameter) {
  41 + CorrelationProcess::CorrelationProcess(const CorrelationProcess& pProcess, Parameter &parameter) : MultiParamProcess_CRTP(pProcess, parameter),
  42 + _firstParamName(pProcess._firstParamName), _secondParamName(pProcess._secondParamName) {
42 43 _type = "";
43 44 }
44 45  
45 46 CorrelationProcess::~CorrelationProcess() {
46   - if (_secondInputParam != nullptr)
47   - _secondInputParam->closeConnection(this);
48 47 }
49 48  
50   - void CorrelationProcess::establishConnection() {
51   -
52   -
53   - if (_attributList.size() < 2) {
54   - BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_PROCESS_ERR) << AMDA::ex_msg(std::string("CorrelationProcess::parse require 3 attributes'")));
55   - }
56   -
57   - //Imf parameter
58   - ParameterCreatorFromExpression creator(_parameter.getParameterManager());
59   - std::string secondParamExpression;
60   -
61   - if(_attributList[0][0] == '$' || _attributList[0][0] == '#')
62   - secondParamExpression=_attributList[0];
63   - else
64   - secondParamExpression = "$" + _attributList[0];
65   - _secondInputParam = creator.getOneParameterFromExpression(_parameter,secondParamExpression ,true);
66   - if (_secondInputParam == nullptr) {
67   - BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_PROCESS_ERR) << AMDA::ex_msg(std::string("CorrelationProcess::parse cannot retrieve imf param")));
68   - }
69   - _secondInputParam->openConnection(this);
  49 + void CorrelationProcess::parse() {
  50 + ParameterCreatorFromExpression creator(_parameter.getParameterManager());
  51 + ParameterSPtr lParameter = creator.getOneParameterFromExpression(_parameter,_expression, isUserProcess());
  52 + _paramNameList[lParameter->getId()].first = lParameter;
  53 + _firstParamName = lParameter->getId();
  54 + lParameter = creator.getOneParameterFromExpression(_parameter,_attributList[0], isUserProcess());
  55 + _paramNameList[lParameter->getId()].first = lParameter;
  56 + _secondParamName = lParameter->getId();
  57 + }
70 58  
71   - SingleParamProcess::establishConnection();
  59 + void CorrelationProcess::establishConnection() {
  60 + parse();
  61 + MultiParamProcess::establishConnection();
72 62 }
73 63  
74 64 TimeStamp CorrelationProcess::init() {
... ... @@ -77,17 +67,19 @@ namespace AMDA {
77 67 if(_type.empty() && _attributList.size() >=3)
78 68 _type = _attributList[2];
79 69  
80   - TimeStamp time = _parameterInput->init(this, _timeIntervalList);
81   - _paramInput = _parameterInput->getParamData(this).get();
  70 + TimeStamp timeStamp = MultiParamProcess::init();
  71 +
  72 + ParamData* lfirstParamInput = _paramNameList[_firstParamName].first->getParamData(this).get();
  73 +
  74 + ParamData* lSecondParamInput = _paramNameList[_secondParamName].first->getParamData(this).get();
  75 +
  76 + StatisticCorrelationCreator lCreator(*this, _timeIntervalList, *lfirstParamInput, *lSecondParamInput, _windowtime, _type);
82 77  
83   - _secondInputParam->init(this, _timeIntervalList);
84   - ParamData* lSecondParamInput = _secondInputParam->getParamData(this).get();
  78 + _operation = lCreator.getOperation();
  79 + _paramData = ParamDataSPtr(_operation->getParamOutput());
  80 + _paramData->setMinSampling(lfirstParamInput->getMinSampling());
85 81  
86   - StatisticCorrelationCreator lCreator(*this, _timeIntervalList, *_paramInput, *lSecondParamInput, _secondInputParam, _windowtime, _type);
87   - _operation = lCreator.getOperation();
88   - _paramData = ParamDataSPtr(_operation->getParamOutput());
89   - _paramData->setMinSampling(_paramInput->getMinSampling());
90   - return time;
  82 + return timeStamp;
91 83 }
92 84  
93 85 // Covariance
... ...
src/ExternLib/StatisticFunctions/CorrelationProcess.hh
  Diff comments   View file @e336270
... ... @@ -14,13 +14,13 @@
14 14 #ifndef CORRELATIONPROCESS_HH
15 15 #define CORRELATIONPROCESS_HH
16 16  
17   -#include "SingleParamProcess.hh"
  17 +#include "MultiParamProcess.hh"
18 18 #include "ParamInfo.hh"
19 19  
20 20 namespace AMDA {
21 21 namespace Parameters {
22 22  
23   - class CorrelationProcess : public SingleParamProcess_CRTP<CorrelationProcess> {
  23 + class CorrelationProcess : public MultiParamProcess_CRTP<CorrelationProcess> {
24 24 public:
25 25 CorrelationProcess(Parameter &parameter);
26 26 CorrelationProcess(const CorrelationProcess& pProcess, Parameter &parameter);
... ... @@ -38,19 +38,21 @@ namespace AMDA {
38 38  
39 39  
40 40 protected:
  41 + /**
  42 + * @brief If the expression is not a Single parameter,
  43 + * it must ask the creation of a parameter responsible of the formula calculation.
  44 + */
  45 + void parse();
41 46  
42   - // void getSecondParamData();
43   -
  47 +
44 48 std::string _type;
45 49  
  50 + private:
46 51 double _windowtime;
47 52  
48   - private:
49   - /**
50   - * @brief list of param name intput
51   - * @detail this list must be ordered
52   - */
53   - ParameterSPtr _secondInputParam;
  53 + std::string _firstParamName;
  54 +
  55 + std::string _secondParamName;
54 56  
55 57 };
56 58  
... ...
src/ExternLib/StatisticFunctions/StatisticCorrelationCreator.hh
  Diff comments   View file @e336270
... ... @@ -30,8 +30,8 @@ namespace AMDA {
30 30 /**
31 31 * @brief Constructor.
32 32 */
33   - StatisticCorrelationCreator(Process& pProcess, TimeIntervalListSPtr pTimeIntervalList, ParamData& paramInput, ParamData& secondParamData, ParameterSPtr secondInputParam, double windowtime, std::string type_)
34   - : _process(pProcess), _timeIntervalList(pTimeIntervalList), _paramData(paramInput), _secondParamData(secondParamData), _secondInputParam(secondInputParam),_windowtime(windowtime), _type(type_) {
  33 + StatisticCorrelationCreator(Process& pProcess, TimeIntervalListSPtr pTimeIntervalList, ParamData& paramInput, ParamData& secondParamData, double windowtime, std::string type_)
  34 + : _process(pProcess), _timeIntervalList(pTimeIntervalList), _paramData(paramInput), _secondParamData(secondParamData), _windowtime(windowtime), _type(type_) {
35 35 _paramData.accept(*this);
36 36 }
37 37  
... ... @@ -172,41 +172,20 @@ namespace AMDA {
172 172  
173 173 template <typename Type>
174 174 void createOperation() {
175   - _operation = new StatisticFunctions::Correlation<Type, Type>(_process, _timeIntervalList, dynamic_cast<ParamDataSpec<Type>&> (_paramData),
  175 + _operation = new StatisticFunctions::Correlation<Type, double>(_process, _timeIntervalList, dynamic_cast<ParamDataSpec<Type>&> (_paramData),
176 176 dynamic_cast<ParamDataSpec<Type>&> (_secondParamData), _windowtime, _type);
177   - pushSecondParamData<Type>();
178   -
179   - }
180   -
181   - template <typename Type>
182   - void pushSecondParamData(){
183   - try {
184   - ParamDataIndexInfo lParamDataIndexInfo;
185   - lParamDataIndexInfo = _secondInputParam->getAsync(&_process).get();
186   - while ((!lParamDataIndexInfo._noMoreTimeInt && !lParamDataIndexInfo._timeIntToProcessChanged) || (lParamDataIndexInfo._nbDataToProcess > 0)) {
187   - reinterpret_cast<StatisticFunctions::AbstractCorrelationFunc<Type, Type>*> (_operation)->pushSecondParamData(lParamDataIndexInfo);
188   - if (lParamDataIndexInfo._timeIntToProcessChanged || lParamDataIndexInfo._noMoreTimeInt)
189   - break;
190   - lParamDataIndexInfo = _secondInputParam->getAsync(&_process).get();
191   - }
192   - } catch (...) {
193   - throw;
194   - }
195   - }
  177 + }
196 178  
197 179 template <typename Type>
198 180 void createOperation1D() {
199   - _operation = new StatisticFunctions::Correlation1D<std::vector<Type>, std::vector<Type>, Type>(_process, _timeIntervalList, dynamic_cast<ParamDataSpec<std::vector<Type>>&> (_paramData),
  181 + _operation = new StatisticFunctions::Correlation<std::vector<Type>, std::vector<double>>(_process, _timeIntervalList, dynamic_cast<ParamDataSpec<std::vector<Type>>&> (_paramData),
200 182 dynamic_cast<ParamDataSpec<std::vector<Type>>&> (_secondParamData), _windowtime, _type);
201   - pushSecondParamData<std::vector<Type>>();
202   -
203 183 }
204 184  
205 185 Process &_process;
206 186 TimeIntervalListSPtr _timeIntervalList;
207 187 ParamData &_paramData;
208 188 ParamData &_secondParamData;
209   - ParameterSPtr _secondInputParam;
210 189 double _windowtime;
211 190 std::string _type;
212 191 Operation *_operation;
... ...
src/ExternLib/StatisticFunctions/Toolbox.hh
  Diff comments   View file @e336270
... ... @@ -191,25 +191,29 @@ namespace AMDA {
191 191 }
192 192  
193 193 template <typename Type>
194   - std::pair<Type, Type> getMean(std::list<std::pair<Type, Type>> &list) {
195   - std::pair<Type, Type> result(0, 0);
  194 + std::pair<double, double> getMean(std::list<std::pair<Type, Type>> &list) {
  195 + std::pair<double, double> result(0, 0);
196 196 std::pair<int, int> counter(0, 0);
197 197  
198 198 for (auto elt : list) {
199 199 if (!isNAN(elt.first)) {
200   - result.first += elt.first;
  200 + result.first += (double)elt.first;
201 201 counter.first += 1;
202 202 }
203 203 if (!isNAN(elt.second)) {
204   - result.second += elt.second;
  204 + result.second += (double)elt.second;
205 205 counter.second += 1;
206 206 }
207 207 }
208 208 if (counter.first != 0)
209 209 result.first /= counter.first;
  210 + else
  211 + result.first = NAN;
210 212  
211 213 if (counter.second != 0)
212 214 result.second /= counter.second;
  215 + else
  216 + result.second = NAN;
213 217  
214 218 return result;
215 219 }
... ... @@ -244,34 +248,28 @@ namespace AMDA {
244 248 }
245 249  
246 250 template <typename Type>
247   - std::vector<Type>rankify(std::vector<Type>& X) {
248   -
249   - int N = X.size();
250   -
251   - // Rank Vector
252   - std::vector<Type> Rank_X(N);
253   -
254   - for (int i = 0; i < N; i++) {
  251 + std::vector<double>rankify(std::vector<Type>& X) {
  252 + std::vector<double> Rank_X;
  253 + if (X.empty()) {
  254 + return Rank_X;
  255 + }
  256 + Rank_X.resize(X.size());
  257 + int i = 0;
  258 + for (typename std::vector<Type>::iterator it1 = X.begin(); it1 != X.end(); ++it1) {
255 259 int r = 1, s = 1;
256 260  
257   - // Count no of smaller elements
258   - // in 0 to i-1
259   - for (int j = 0; j < N; j++) {
260   - if (X[j] < X[i]) r++;
261   - if (X[j] == X[i] && i != j) s++;
  261 + for (typename std::vector<Type>::iterator it2 = X.begin(); it2 != X.end(); ++it2) {
  262 + if (it1 == it2) {
  263 + continue;
  264 + }
  265 + if (*it2 < *it1) r++;
  266 + if (*it1 == *it2) s++;
262 267 }
263 268  
264   - // Count no of smaller elements
265   - // in i+1 to N-1
266   - /**
267   - for (int j = i + 1; j < N; j++) {
268   - if (X[j] < X[i]) r++;
269   - if (X[j] == X[i]) s++;
270   - }*/
271   -
272 269 // Use Fractional Rank formula
273 270 // fractional_rank = r + (n-1)/2
274   - Rank_X[i] = (Type) r + (s - 1) * 0.5;
  271 + Rank_X[i] = (double) r + (s - 1) * 0.5;
  272 + ++i;
275 273 }
276 274  
277 275 // Return Rank Vector
... ... @@ -279,107 +277,124 @@ namespace AMDA {
279 277 }
280 278  
281 279 template <typename Type>
282   - bool getCovariance(std::list<std::pair<Type, Type>> &list, Type &result) {
  280 + bool getCovariance(std::list<std::pair<Type, Type>> &list, double &result) {
  281 + result = NAN;
283 282 if (list.empty()) {
284 283 return false;
285 284 }
286   - std::pair<Type, Type> mean = getMean(list);
287   - result = 0;
  285 + std::pair<double, double> mean = getMean(list);
  286 + if (isNAN(mean.first) || isNAN(mean.second)) {
  287 + return false;
  288 + }
288 289 int counter = 0;
  290 + double sum = 0.;
289 291 for (auto elt : list) {
290 292 if (!isNAN(elt.first) && !isNAN(elt.second)) {
291   - result += (elt.first - mean.first)*(elt.second - mean.second);
  293 + sum += ((double)elt.first - mean.first)*((double)elt.second - mean.second);
292 294 counter += 1;
293 295 }
294 296 }
295 297 if (counter != 0)
296   - result /= counter;
  298 + result = sum / counter;
297 299 return true;
298 300 }
299 301  
300 302 template <typename Type>
301   - bool getPearson(std::list<std::pair<Type, Type>> &list, Type &result) {
  303 + bool getPearson(std::list<std::pair<Type, Type>> &list, double &result) {
  304 + result = NAN;
302 305 if (list.empty()) {
303 306 return false;
304 307 }
305   - result = 0;
306 308 int counter = 0;
307   - Type sum1 = 0, sum2 = 0;
308   - Type sum1Sq = 0, sum2Sq = 0;
309   - Type sum12 = 0;
  309 + double sum1 = 0, sum2 = 0;
  310 + double sum1Sq = 0, sum2Sq = 0;
  311 + double sum12 = 0;
  312 +
310 313 for (auto elt : list) {
311 314 if (!isNAN(elt.first) && !isNAN(elt.second)) {
312   - sum1 += elt.first;
313   - sum2 += elt.second;
314   - sum1Sq += elt.first * elt.first;
315   - sum2Sq += elt.second * elt.second;
316   - sum12 += elt.first * elt.second;
  315 + sum1 += (double)elt.first;
  316 + sum2 += (double)elt.second;
  317 + sum1Sq += (double)elt.first * (double)elt.first;
  318 + sum2Sq += (double)elt.second * (double)elt.second;
  319 + sum12 += (double)elt.first * (double)elt.second;
317 320 counter += 1;
318 321 }
319 322 }
320   - if (counter > 0)
  323 + if (counter > 1) {
321 324 result = (counter * sum12 - sum1 * sum2) /
322 325 (std::sqrt((counter * sum1Sq - sum1 * sum1)*(counter * sum2Sq - sum2 * sum2)));
  326 + }
323 327 return true;
324 328 }
325 329  
326 330 template <typename Type>
327   - bool getSpearman(std::list<std::pair<Type, Type>> &list, Type &result) {
  331 + bool getSpearman(std::list<std::pair<Type, Type>> &list, double &result) {
  332 + result = NAN;
  333 +
328 334 if (list.empty()) {
329 335 return false;
330 336 }
331 337  
332   - std::vector<Type> X, Y;
  338 + std::vector<double> X, Y;
333 339 for (auto elt : list) {
334 340 if (!isNAN(elt.first) && !isNAN(elt.second)) {
335   - X.push_back(elt.first);
336   - Y.push_back(elt.second);
  341 + X.push_back((double)elt.first);
  342 + Y.push_back((double)elt.second);
337 343 }
338 344 }
339   - int n = X.size();
340   - if (n == 0)
341   - return true;
342   - std::vector<Type> rank_X = rankify(X);
343   - std::vector<Type> rank_Y = rankify(Y);
344   - std::list<std::pair<Type, Type>> rankList;
  345 +
  346 + if (X.empty()) {
  347 + return false;
  348 + }
  349 +
  350 + std::vector<double> rank_X = rankify(X);
  351 + std::vector<double> rank_Y = rankify(Y);
  352 + std::list<std::pair<double, double>> rankList;
345 353  
346   - for (int i = 0; i < n; i++)
  354 + for (int i = 0; i < X.size(); i++)
347 355 rankList.push_back(std::make_pair(rank_X[i], rank_Y[i]));
348 356  
349   - result = 0;
350   - getPearson(rankList, result);
351   - return true;
  357 + return getPearson(rankList, result);
352 358 }
353 359  
354 360 template <typename Type>
355   - bool getKendall(std::list<std::pair<Type, Type>> &list, Type &result) {
  361 + bool getKendall(std::list<std::pair<Type, Type>> &list, double &result) {
  362 + result = NAN;
  363 +
356 364 if (list.empty()) {
357 365 return false;
358 366 }
359   - result = 0;
360   - std::vector<Type> X, Y;
  367 +
  368 + std::vector<double> X, Y;
361 369 for (auto elt : list) {
362 370 if (!isNAN(elt.first) && !isNAN(elt.second)) {
363   - X.push_back(elt.first);
364   - Y.push_back(elt.second);
  371 + X.push_back((double)elt.first);
  372 + Y.push_back((double)elt.second);
365 373 }
366 374 }
367   - int n = X.size();
368   - if (n == 0)
369   - return true;
370   - std::vector<Type> rank_x = rankify(X);
371   - std::vector<Type> rank_y = rankify(Y);
372   - Type nc = 0;
373   - for (int i = 0; i < n; i++){
374   - for (int j=i;j<n;j++){
375   - if((i != j )&& ((rank_x[i] > rank_x[j] && rank_y[i] > rank_y[j]) ||
376   - (rank_x[i] < rank_x[j] && rank_y[i] < rank_y[j])))
377   - nc += 1;
  375 +
  376 + if (X.empty()) {
  377 + return false;
  378 + }
  379 +
  380 + std::vector<double> rank_x = rankify(X);
  381 + std::vector<double> rank_y = rankify(Y);
  382 +
  383 + long nc = 0;
  384 + long nd = 0;
  385 + for (int i = 0; i < rank_x.size()-1; i++){
  386 + for (int j = i+1; j<rank_x.size(); j++){
  387 + if (((rank_x[i] > rank_x[j]) && (rank_y[i] > rank_y[j])) ||
  388 + ((rank_x[i] < rank_x[j]) && (rank_y[i] < rank_y[j])))
  389 + ++nc;
  390 + if (((rank_x[i] > rank_x[j]) && (rank_y[i] < rank_y[j])) ||
  391 + ((rank_x[i] < rank_x[j]) && (rank_y[i] > rank_y[j])))
  392 + ++nd;
  393 +
378 394 }
379 395  
380 396 }
381   - result = (Type) (4*nc - n*(n-1))/(Type) (n * (n - 1));
382   -
  397 + result = ((double)(nc-nd)/(0.5*rank_x.size()*(rank_x.size()-1)));
383 398 return true;
384 399 }
385 400  
... ...