App

- Main Page
- Related Pages
- Namespaces
- Classes
- Files

- File List
- File Members

CubicSpline.cpp
Go to the documentation of this file.
1
2 #include "CubicSpline.hpp"
3
4 namespace app {
5
6 CubicSpline::CubicSpline(double const * x, double const * y, unsigned count, double y1_0, double y1_n) {
7 if (count < 2) { initOK = 0; return; }
8 valsCount = count;
9 lastVal = count - 1;
10 jsav = 0; cor = 0;
11
12 if (x == NULL) { initOK = 0; return; }
13 if (y == NULL) { initOK = 0; return; }
14
15 xVals = (double *) malloc((unsigned long int)(valsCount * sizeof(double))); // ANDROID_SPECIFIC přetypování v malloc
16 if (xVals == NULL) { initOK = 0; return; }
17 memcpy(xVals, x, (unsigned long int)(count * sizeof(double))); // xVals = x;
18
19 yVals = (double *) malloc((unsigned long int)(valsCount * sizeof(double))); // ANDROID_SPECIFIC přetypování v malloc
20 if (yVals == NULL) { initOK = 0; free(xVals); return; }
21 memcpy(yVals, y, (unsigned long int)(count * sizeof(double))); // yVals = y; // ANDROID_SPECIFIC přetypování v malloc
22
23 y2Vals = (double *) malloc((unsigned long int)(valsCount * sizeof(double))); // ANDROID_SPECIFIC přetypování v malloc
24 if (y2Vals == NULL) { initOK = 0; free(xVals); free(yVals); return; }
25
26 ascnd = (xVals[lastVal] >= xVals[0]);
27 dj = pow(count, 0.25);
28 dj = (dj > 1) ? 1 : dj; //MIN(1,(int)pow((Doub)n,0.25)); // WTF???
29 if ( y2ValsInit(xVals, yVals, y1_0, y1_n) ) { initOK = 0; free(xVals); free(yVals); free(y2Vals); return; }
30 initOK = 1;
31 }
32
33 CubicSpline::~CubicSpline() {
34 if (xVals != NULL) { free(xVals); }
35 if (yVals != NULL) { free(yVals); }
36 if (y2Vals != NULL) { free(y2Vals); }
37 }
38
39 double CubicSpline::getValueAt(double x) {
40 if (initOK) { // Při inicializaci nenastal žádný problém. Díky tomuto zmizí testování v ostatních metodách
41 if (x <= xVals[0]) {
42 return yVals[0];
43 } else if (x >= xVals[lastVal]) {
44 return yVals[lastVal];
45 }
46 int jlo = cor ? hunt(x) : locate(x); // Prom je zbytečná, lze rovnou volat s tímto funkci
47 /*
48 if (jlo < 0) {
49 return 0; // Proč? Co to znamená? NESMYSL - nemůže vrátit menší, zajistěno při návratu
50 }
51 */
52 return interpolate(jlo, x);
53 } else {
54 return 0;
55 }
56 }
57
58 int CubicSpline::locate(const double x) {
59 /* NELZE - pojistka v initOK
60 if (valsCount < 2) {
61 return -1; // throw("locate size error");
62 }
63 */
64 int jl = 0;
65 int ju = lastVal; //
66 int jm = 0;
67 while (ju - jl > 1) {
68 jm = (ju + jl) >> 1;
69 if (x >= xVals[jm] == ascnd) {
70 jl = jm;
71 } else {
72 ju = jm;
73 }
74 }
75 cor = abs(jl - jsav) > dj ? 0 : 1;
76 jsav = jl;
77 int retVal = valsCount - 2;
78 retVal = (retVal < jl) ? retVal : jl;
79 return (retVal > 0) ? retVal : 0; // je kontrola na zápornost nutná? Asi ne.. Ale okrajová optimalizace
80 }
81
82 int CubicSpline::hunt(const double x) {
83 /* NELZE
84 if (valsCount < 2) {
85 return -1; // throw("hunt size error");
86 }
87 */
88 int jl = jsav;
89 int jm, ju;
90 int inc = 1;
91 if (jl < 0 || jl > lastVal) {
92 jl = 0;
93 ju = lastVal;
94 } else {
95 if (x >= xVals[jl] == ascnd) {
96 while (1) {
97 ju = jl + inc;
98 if (ju >= lastVal) {
99 ju = lastVal;
100 break;
101 } else if (x < xVals[ju] == ascnd)
102 break;
103 else {
104 jl = ju;
105 inc += inc;
106 }
107 }
108 } else {
109 ju = jl;
110 while (1) {
111 jl = jl - inc;
112 if (jl <= 0) {
113 jl = 0;
114 break;
115 } else if (x >= xVals[jl] == ascnd)
116 break;
117 else {
118 ju = jl;
119 inc += inc;
120 }
121 }
122 }
123 }
124 while (ju - jl > 1) {
125 jm = (ju + jl) >> 1;
126 if (x >= xVals[jm] == ascnd)
127 jl = jm;
128 else
129 ju = jm;
130 }
131 cor = abs(jl - jsav) > dj ? 0 : 1;
132 jsav = jl;
133 int retVal = valsCount - 2;
134 retVal = (retVal < jl) ? retVal : jl;
135 return (retVal > 0) ? retVal : 0; // MAX(0,MIN(valsCount-2,jl)); // MAX(0,MIN(n-mm,jl-((mm-2)>>1))); // n = count, mm = kernelWidth
136 }
137
138 int CubicSpline::y2ValsInit(const double *xv, const double *yv, double y1_0, double y1_n) {
139 // Jako optimalizace doplnena lastVal, pro pochopeni viz kniha.
140 int i, k;
141 double p, qn, sig, un;
142 double * u = (double *) malloc ((unsigned long int)((lastVal) * sizeof(double))); // ANDROID_SPECIFIC přetypování v malloc
143 if (u == NULL) { return 1; }
144 if (y1_0 > 0.99e99) { // Natural spline
145 y2Vals[0] = u[0] = 0.0;
146 } else {
147 y2Vals[0] = -0.5;
148 u[0] = (3.0 / (xv[1] - xv[0]))
149 * ((yv[1] - yv[0]) / (xv[1] - xv[0]) - y1_0);
150 }
151 for (i = 1; i < lastVal; i++) {
152 sig = (xv[i] - xv[i - 1]) / (xv[i + 1] - xv[i - 1]);
153 p = sig * y2Vals[i - 1] + 2.0;
154 y2Vals[i] = (sig - 1.0) / p;
155 u[i] = (yv[i + 1] - yv[i]) / (xv[i + 1] - xv[i])
156 - (yv[i] - yv[i - 1]) / (xv[i] - xv[i - 1]);
157 u[i] = (6.0 * u[i] / (xv[i + 1] - xv[i - 1]) - sig * u[i - 1]) / p;
158 }
159 if (y1_n > 0.99e99) { // Natural spline
160 qn = un = 0.0;
161 } else {
162 qn = 0.5;
163 un = (3.0 / (xv[lastVal] - xv[valsCount - 2]))
164 * (y1_n - (yv[lastVal] - yv[valsCount - 2]) / (xv[lastVal] - xv[valsCount - 2]));
165 }
166 y2Vals[lastVal] = (un - qn * u[valsCount - 2]) / (qn * y2Vals[valsCount - 2] + 1.0);
167 for (k = valsCount - 2; k >= 0; k--) {
168 y2Vals[k] = y2Vals[k] * y2Vals[k + 1] + u[k];
169 }
170 free(u);
171 return 0;
172 }
173
174 double CubicSpline::interpolate(int jl, double x) const{
175 int klo = jl, khi = jl + 1;
176 double y, h, b, a;
177 h = xVals[khi] - xVals[klo];
178 // if (h == 0.0) { return 0; } //throw("Bad input to routine splint"); // Toto by nikdy nemelo nastat, protoze mam vzdy min. 2 body a vzdy vracim jl o 2 mensi, nez je max...
179 a = (xVals[khi] - x) / h;
180 b = (x - xVals[klo]) / h;
181 y = a * yVals[klo] + b * yVals[khi]
182 + ((a * a * a - a) * y2Vals[klo] + (b * b * b - b) * y2Vals[khi]) * (h * h)
183 / 6.0;
184 return y;
185 }
186
187 } /* namespace app */