Deriv.hh
4.21 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
/*
* Deriv.hh
*
* Created on: Dec 3, 2012
* Author: f.casimir
*/
#ifndef Deriv_HH_
#define Deriv_HH_
#include "Parameter.hh"
#include "ParamData.hh"
#include "DataTypeMath.hh"
#include "Operation.hh"
namespace AMDA {
namespace Parameters {
/**
* @class Deriv
* @brief It is responsible to compute the deriv data.
* @details This class implement the interface Operation.
* @code
* 3 point Lagrangian interpolation (points {x0,y0; x1,y1; x2,y2}):
* df/dx = y0*(2x-x1-x2)/(x01*x02) + y1*(2x-x0-x2)/(x10*x12)+y2*(2x-x0-x1)/(x20*x21),
* where
* x01= x0-x1; x02=x0-x2, x12=x1-x2...
* @endcode
*/
template<class TParamData>
class Deriv : public Operation {
public:
/**
* @brief Element type of paramData
*/
typedef typename TParamData::ElementType ElementType;
/**
* @brief Constructor.
* @details Create the ParamData type of the input ParamData.
*/
Deriv(Process& pProcess, TParamData& paramInput)
: Operation(pProcess),
_paramInput(paramInput),
_paramOutput(new TParamData()),
_firstPointComputed(false) {
_paramDataOutput=_paramOutput;
}
/**
* @overload Operation::write(ParamDataIndexInfo &pParamDataIndexInfo)
*/
void write(ParamDataIndexInfo &pParamDataIndexInfo) {
//Store data in a working list
for (unsigned int index = pParamDataIndexInfo._startIndex;
index< pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess; index++) {
_data.push_back(std::make_pair(_paramInput.getTime(index), _paramInput.getDataList()[index]));
}
ElementType res;
typename std::list<std::pair<double, ElementType> >::iterator itData0, itData1, itData2;
if (!_firstPointComputed && (_data.size() >= 3)) {
//First point calculation
itData0 = _data.begin();
itData1 = std::next(itData0, 1);
itData2 = std::next(itData0, 2);
res = computeDeriv(*itData0, *itData1, *itData2);
_paramOutput->getDataList().push_back(res);
_paramOutput->pushTime((*itData0).first);
_firstPointComputed = true;
}
while (_data.size() >= 3) {
//Calculation for other points than the first one or the last one
itData0 = _data.begin();
itData1 = std::next(itData0, 1);
itData2 = std::next(itData0, 2);
res = computeDeriv(*itData0, *itData1, *itData2);
_paramOutput->getDataList().push_back(res);
_paramOutput->pushTime((*itData1).first);
_lastPopData = *itData0;
_data.pop_front();
}
if (pParamDataIndexInfo._timeIntToProcessChanged || pParamDataIndexInfo._noMoreTimeInt) {
//Finish process for this interval
if (!_firstPointComputed) {
//Only one or two points => cannot compute deriv => set NaN values
while (!_data.empty()) {
res << NotANumber();
_paramOutput->getDataList().push_back(res);
_paramOutput->pushTime((*_data.begin()).first);
_data.pop_front();
}
}
else {
//Last point calculation
itData1 = _data.begin();
itData2 = std::next(itData1, 1);
res = computeDeriv(_lastPopData, *itData1, *itData2);
_paramOutput->getDataList().push_back(res);
_paramOutput->pushTime((*itData2).first);
_data.clear();
}
}
}
ElementType computeDeriv(std::pair<double, ElementType> data0, std::pair<double, ElementType> data1, std::pair<double, ElementType> data2) {
ElementType y0 = data0.second;
ElementType y1 = data1.second;
ElementType y2 = data2.second;
double x0 = data0.first;
double x1 = data1.first;
double x2 = data2.first;
double x = x1;
double x02 = x0 -x2;
double x01 = x0 -x1;
double x10 = x1 -x0;
double x12 = x1 -x2;
double x20 =x2 - x0 ;
double x21 = x2 -x1;
double _2x = 2*x;
return y0*(_2x-x1-x2)/(x01*x02) + y1*(_2x-x0-x2)/(x10*x12)+y2*(_2x-x0-x1)/(x20*x21);
}
/**
* @overload Operation::reset(double pStartTime, double pTimeInt)
* @brief reset static data to process another TimeInterval
*/
virtual void reset() {
Operation::reset();
_data.clear();
_firstPointComputed = false;
}
private:
/**<
* @brief It is the channel of data derived
*/
TParamData &_paramInput;
/**<
* @brief It is the channel of the data derived
*/
TParamData *_paramOutput;
bool _firstPointComputed;
std::pair<double, ElementType> _lastPopData;
std::list<std::pair<double, ElementType> > _data;
};
} /* namespace Parameters */
} /* namespace AMDA */
#endif /* Deriv_HH_ */