App
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 
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 */