Skip to content

Commit 6e0399d

Browse files
pasabanovpasabanov
committed
Unstable "Not-a-knot" cubic spline
1 parent d5c73a9 commit 6e0399d

File tree

3 files changed

+82
-149
lines changed

3 files changed

+82
-149
lines changed

‎ALFI/ALFI/spline/cubic.h‎

Lines changed: 75 additions & 112 deletions
Original file line numberDiff line numberDiff line change
@@ -115,66 +115,37 @@ namespace alfi::spline {
115115
coeffs = util::spline::simple_spline<Number,Container>(X, Y, 3);
116116
return;
117117
}
118-
119-
const Container<Number> dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
120-
121-
Container<Number> B(n);
122-
118+
const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
123119
Container<Number> diag(n - 2), right(n - 2);
124-
diag[0] = ((dX[0] + dX[1]) * (dX[0] + 2 * dX[1])) / dX[1];
125-
for (SizeT i = 1; i < n - 3; ++i) {
126-
diag[i] = 2 * (dX[i] + dX[i+1]);
127-
}
128-
diag[n-3] = ((dX[n-2] + dX[n-3]) * (dX[n-2] + 2 * dX[n-3])) / dX[n-3];
129120
for (SizeT i = 0; i < n - 2; ++i) {
121+
diag[i] = 2 * (dX[i] + dX[i+1]);
130122
right[i] = 3 * (dY[i+1]/dX[i+1] - dY[i]/dX[i]);
131123
}
132-
133-
const auto m_1 = dX[1] / diag[0];
134-
diag[1] -= m_1 * (dX[1] - dX[0]*dX[0]/dX[1]);
135-
right[1] -= m_1 * right[0];
136-
for (SizeT i = 2; i < n - 3; ++i) {
137-
const auto m = dX[i] / diag[i-1];
138-
diag[i] -= m * dX[i];
139-
right[i] -= m * right[i-1];
124+
const auto C = util::linalg::tridiag_solve(
125+
dX[0] - dX[1], 2*dX[0] + dX[1], dX[0] / (dX[0]+dX[1]) * right[0],
126+
2*dX[n-2] + dX[n-3], dX[n-2] - dX[n-3], dX[n-2] / (dX[n-2]+dX[n-3]) * right[n-3],
127+
dX.begin(), diag.begin(), dX.begin() + 1, right.begin(), n
128+
);
129+
for (SizeT i = 0, index = 0; i < n - 1; ++i) {
130+
coeffs[index++] = (C[i+1] - C[i]) / (3*dX[i]);
131+
coeffs[index++] = C[i];
132+
coeffs[index++] = dY[i]/dX[i] - dX[i] * ((2*C[i] + C[i+1]) / 3);
133+
coeffs[index++] = Y[i];
140134
}
141-
const auto m_n_3 = (dX[n-3] - dX[n-2]*dX[n-2]/dX[n-3]) / diag[n-4];
142-
diag[n-3] -= m_n_3 * dX[n-3];
143-
right[n-3] -= m_n_3 * right[n-4];
144-
145-
B[n-2] = right[n-3] / diag[n-3];
146-
for (SizeT i = n - 3; i >= 2; --i) {
147-
B[i] = (right[i-1] - dX[i] * B[i+1]) / diag[i-1];
148-
}
149-
B[1] = (right[0] - (dX[1] - dX[0]*dX[0]/dX[1]) * B[2]) / diag[0];
150-
B[0] = B[1] + dX[0]/dX[1] * (B[1] - B[2]);
151-
B[n-1] = B[n-2] + dX[n-2]/dX[n-3] * (B[n-2] - B[n-3]);
152-
153-
Container<Number> A(n - 1);
154-
for (SizeT i = 1; i < A.size() - 1; ++i) {
155-
A[i] = (B[i+1] - B[i]) / (3*dX[i]);
156-
}
157-
A[0] = A[1];
158-
A[n-2] = A[n-3];
159-
Container<Number> C(n - 1);
160-
for (SizeT i = 0; i < C.size(); ++i) {
161-
C[i] = dY[i]/dX[i] - dX[i] * ((2*B[i] + B[i+1]) / 3);
162-
}
163-
util::spline::merge_coeffs(coeffs, {A, B, C, Y});
164135
},
165136
[&](const typename Types::Periodic&) {
166137
if (n <= 2) {
167138
coeffs = util::spline::simple_spline<Number,Container>(X, Y, 3);
168139
return;
169140
}
170141

171-
const Container<Number> dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
142+
const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
172143

173-
Container<Number> B(n);
144+
Container<Number> C(n);
174145
if (n == 3) {
175-
B[0] = 3 * (dY[0]/dX[0] - dY[1]/dX[1]) / (dX[0] + dX[1]);
176-
B[1] = -B[0];
177-
B[2] = B[0];
146+
C[0] = 3 * (dY[0]/dX[0] - dY[1]/dX[1]) / (dX[0] + dX[1]);
147+
C[1] = -C[0];
148+
C[2] = C[0];
178149
} else {
179150
Container<Number> diag(n - 1), right(n - 1);
180151
for (SizeT i = 0; i < n - 2; ++i) {
@@ -204,38 +175,35 @@ namespace alfi::spline {
204175
diag[n-2] -= last_row[n-3] / diag[n-3] * last_col[n-3];
205176
right[n-2] -= last_row[n-3] / diag[n-3] * right[n-3];
206177

207-
B[n-1] = right[n-2] / diag[n-2];
208-
B[n-2] = (right[n-3] - last_col[n-3] * B[n-1]) / diag[n-3];
178+
C[n-1] = right[n-2] / diag[n-2];
179+
C[n-2] = (right[n-3] - last_col[n-3] * C[n-1]) / diag[n-3];
209180
for (SizeT i = n - 3; i >= 1; --i) {
210-
B[i] = (right[i-1] - dX[i]*B[i+1] - last_col[i-1]*B[n-1]) / diag[i-1];
181+
C[i] = (right[i-1] - dX[i]*C[i+1] - last_col[i-1]*C[n-1]) / diag[i-1];
211182
}
212-
B[0] = B[n-1];
183+
C[0] = C[n-1];
213184
}
214185

215-
Container<Number> A(n - 1);
216-
for (SizeT i = 0; i < A.size(); ++i) {
217-
A[i] = (B[i+1] - B[i]) / (3*dX[i]);
186+
for (SizeT i = 0, index = 0; i < n - 1; ++i) {
187+
coeffs[index++] = (C[i+1] - C[i]) / (3*dX[i]);
188+
coeffs[index++] = C[i];
189+
coeffs[index++] = dY[i]/dX[i] - dX[i] * ((2*C[i] + C[i+1]) / 3);
190+
coeffs[index++] = Y[i];
218191
}
219-
Container<Number> C(n - 1);
220-
for (SizeT i = 0; i < C.size(); ++i) {
221-
C[i] = dY[i]/dX[i] - dX[i] * ((2*B[i] + B[i+1]) / 3);
222-
}
223-
util::spline::merge_coeffs(coeffs, {A, B, C, Y});
224192
},
225193
[&](const typename Types::ParabolicEnds&) {
226194
coeffs = compute_coeffs(X, Y, typename Types::FixedThird{0, 0});
227195
},
228196
[&](const typename Types::Clamped& c) {
229-
const Container<Number> dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
197+
const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
230198

231-
Container<Number> B(n);
199+
Container<Number> C(n);
232200
if (n == 2) {
233-
B[0] = (3*dY[0]/dX[0] - c.df_n - 2*c.df_1) / dX[0];
234-
B[1] = (c.df_1 + 2*c.df_n - 3*dY[0]/dX[0]) / dX[0];
201+
C[0] = (3*dY[0]/dX[0] - c.df_n - 2*c.df_1) / dX[0];
202+
C[1] = (c.df_1 + 2*c.df_n - 3*dY[0]/dX[0]) / dX[0];
235203
} else if (n == 3) {
236-
B[1] = (3*(dY[1]/dX[1] - dY[0]/dX[0]) + (c.df_1 - c.df_n)) / (dX[0] + dX[1]);
237-
B[0] = (3*(dY[0]/dX[0] - c.df_1)/dX[0] - B[1]) / 2;
238-
B[2] = (3*(c.df_n - dY[1]/dX[1])/dX[1] - B[1]) / 2;
204+
C[1] = (3*(dY[1]/dX[1] - dY[0]/dX[0]) + (c.df_1 - c.df_n)) / (dX[0] + dX[1]);
205+
C[0] = (3*(dY[0]/dX[0] - c.df_1)/dX[0] - C[1]) / 2;
206+
C[2] = (3*(c.df_n - dY[1]/dX[1])/dX[1] - C[1]) / 2;
239207
} else {
240208
Container<Number> diag(n - 2), right(n - 2);
241209
diag[0] = 1.5*dX[0] + 2*dX[1];
@@ -253,66 +221,60 @@ namespace alfi::spline {
253221
right[i] -= m * right[i-1];
254222
}
255223

256-
B[n-2] = right[n-3] / diag[n-3];
224+
C[n-2] = right[n-3] / diag[n-3];
257225
for (SizeT i = n - 3; i >= 1; --i) {
258-
B[i] = (right[i-1] - dX[i] * B[i+1]) / diag[i-1];
226+
C[i] = (right[i-1] - dX[i] * C[i+1]) / diag[i-1];
259227
}
260-
B[n-1] = 1.5*(c.df_n/dX[n-2] - dY[n-2]/dX[n-2]/dX[n-2]) - 0.5*B[n-2];
261-
B[0] = 1.5*(dY[0]/dX[0]/dX[0] - c.df_1/dX[0]) - 0.5*B[1];
228+
C[n-1] = 1.5*(c.df_n/dX[n-2] - dY[n-2]/dX[n-2]/dX[n-2]) - 0.5*C[n-2];
229+
C[0] = 1.5*(dY[0]/dX[0]/dX[0] - c.df_1/dX[0]) - 0.5*C[1];
262230
}
263231

264-
Container<Number> A(n - 1);
265-
for (SizeT i = 0; i < A.size(); ++i) {
266-
A[i] = (B[i+1] - B[i]) / (3*dX[i]);
267-
}
268-
Container<Number> C(n - 1);
269-
for (SizeT i = 0; i < C.size(); ++i) {
270-
C[i] = dY[i]/dX[i] - dX[i] * ((2*B[i] + B[i+1]) / 3);
232+
for (SizeT i = 0, index = 0; i < n - 1; ++i) {
233+
coeffs[index++] = (C[i+1] - C[i]) / (3*dX[i]);
234+
coeffs[index++] = C[i];
235+
coeffs[index++] = dY[i]/dX[i] - dX[i] * ((2*C[i] + C[i+1]) / 3);
236+
coeffs[index++] = Y[i];
271237
}
272-
util::spline::merge_coeffs(coeffs, {A, B, C, Y});
273238
},
274239
[&](const typename Types::FixedSecond& s) {
275-
const Container<Number> dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
240+
const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
276241

277-
Container<Number> B(n);
278-
B[0] = s.d2f_1 / 2;
279-
B[n-1] = s.d2f_n / 2;
242+
Container<Number> C(n);
243+
C[0] = s.d2f_1 / 2;
244+
C[n-1] = s.d2f_n / 2;
280245
if (n > 2) {
281246
Container<Number> diag(n - 2), right(n - 2);
282247
for (SizeT i = 0; i < n - 2; ++i) {
283248
diag[i] = 2 * (dX[i] + dX[i+1]);
284249
right[i] = 3 * (dY[i+1]/dX[i+1] - dY[i]/dX[i]);
285250
}
286-
right[0] -= B[0] * dX[0];
287-
right[n-3] -= B[n-1] * dX[n-2];
251+
right[0] -= C[0] * dX[0];
252+
right[n-3] -= C[n-1] * dX[n-2];
288253
for (SizeT i = 1; i < n - 2; ++i) {
289254
const auto m = dX[i] / diag[i-1];
290255
diag[i] -= m * dX[i];
291256
right[i] -= m * right[i-1];
292257
}
293-
B[n-2] = right[n-3] / diag[n-3];
258+
C[n-2] = right[n-3] / diag[n-3];
294259
for (SizeT i = n - 3; i >= 1; --i) {
295-
B[i] = (right[i-1] - dX[i] * B[i+1]) / diag[i-1];
260+
C[i] = (right[i-1] - dX[i] * C[i+1]) / diag[i-1];
296261
}
297262
}
298263

299-
Container<Number> A(n - 1);
300-
for (SizeT i = 0; i < A.size(); ++i) {
301-
A[i] = (B[i+1] - B[i]) / (3*dX[i]);
302-
}
303-
Container<Number> C(n - 1);
304-
for (SizeT i = 0; i < C.size(); ++i) {
305-
C[i] = dY[i]/dX[i] - dX[i] * ((2*B[i] + B[i+1]) / 3);
264+
for (SizeT i = 0, index = 0; i < n - 1; ++i) {
265+
coeffs[index++] = (C[i+1] - C[i]) / (3*dX[i]);
266+
coeffs[index++] = C[i];
267+
coeffs[index++] = dY[i]/dX[i] - dX[i] * ((2*C[i] + C[i+1]) / 3);
268+
coeffs[index++] = Y[i];
306269
}
307-
util::spline::merge_coeffs(coeffs, {A, B, C, Y});
308270
},
309271
[&](const typename Types::FixedThird& t) {
310-
const Container<Number> dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
272+
const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
311273

312-
Container<Number> B(n);
274+
Container<Number> C(n);
313275
if (n == 2) {
314-
B[0] = -dX[0] * (t.d3f_1 + t.d3f_n) / 8;
315-
B[1] = -B[0];
276+
C[0] = -dX[0] * (t.d3f_1 + t.d3f_n) / 8;
277+
C[1] = -C[0];
316278
} else {
317279
Container<Number> diag(n), right(n);
318280
diag[0] = 1;
@@ -338,36 +300,37 @@ namespace alfi::spline {
338300
diag[n-1] -= m_n_1 * dX[n-2];
339301
right[n-1] -= m_n_1 * right[n-2];
340302

341-
B[n-1] = right[n-1] / diag[n-1];
303+
C[n-1] = right[n-1] / diag[n-1];
342304
for (SizeT i = n - 2; i >= 1; --i) {
343-
B[i] = (right[i] - dX[i] * B[i+1]) / diag[i];
305+
C[i] = (right[i] - dX[i] * C[i+1]) / diag[i];
344306
}
345-
B[0] = right[0] + B[1];
307+
C[0] = right[0] + C[1];
346308
}
347309

348-
Container<Number> A(n - 1);
349-
for (SizeT i = 1; i < A.size() - 1; ++i) {
350-
A[i] = (B[i+1] - B[i]) / (3*dX[i]);
310+
Container<Number> D(n - 1);
311+
for (SizeT i = 1; i < D.size() - 1; ++i) {
312+
D[i] = (C[i+1] - C[i]) / (3*dX[i]);
351313
}
352314
if (n == 2) {
353-
A[0] = (t.d3f_1 + t.d3f_n) / 12;
315+
D[0] = (t.d3f_1 + t.d3f_n) / 12;
354316
} else {
355-
A[0] = t.d3f_1 / 6;
356-
A[n-2] = t.d3f_n / 6;
317+
D[0] = t.d3f_1 / 6;
318+
D[n-2] = t.d3f_n / 6;
357319
}
358-
Container<Number> C(n - 1);
359-
for (SizeT i = 0; i < C.size(); ++i) {
360-
C[i] = dY[i]/dX[i] - dX[i] * ((2*B[i] + B[i+1]) / 3);
320+
for (SizeT i = 0, index = 0; i < n - 1; ++i) {
321+
coeffs[index++] = D[i];
322+
coeffs[index++] = C[i];
323+
coeffs[index++] = dY[i]/dX[i] - dX[i] * ((2*C[i] + C[i+1]) / 3);
324+
coeffs[index++] = Y[i];
361325
}
362-
util::spline::merge_coeffs(coeffs, {A, B, C, Y});
363326
},
364327
[&](const typename Types::NotAKnotStart&) {
365328
if (n <= 4) {
366329
coeffs = util::spline::simple_spline<Number,Container>(X, Y, 3);
367330
return;
368331
}
369332

370-
const Container<Number> dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
333+
const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
371334

372335
{ // first four points
373336
const auto h1 = dX[0], h2 = dX[1];
@@ -414,7 +377,7 @@ namespace alfi::spline {
414377
return;
415378
}
416379

417-
const Container<Number> dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
380+
const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
418381

419382
{ // last four points
420383
const auto h1 = dX[n-4], h2 = dX[n-3];

‎ALFI/ALFI/util/linalg.h‎

Lines changed: 7 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -64,43 +64,22 @@ namespace alfi::util::linalg {
6464
return X;
6565
}
6666

67-
// WRONG CODE
68-
// Tried to implement this for "Not-a-knot" spline, because it had the third coefficient in the first equation,
69-
// but decided that it would be better to derive and equation with only two coefficients.
70-
// See https://github.com/ALFI-lib/alfi-lib.github.io/pull/8
71-
template <typename Number = DefaultNumber, template <typename> class Container = DefaultContainer, typename I1, typename I2, typename I3, typename I4>
67+
template <typename Number = DefaultNumber, template <typename> class Container = DefaultContainer>
7268
Container<Number> tridiag_solve(
73-
Number a11, Number a12, Number a13, Number r1,
74-
Number ann2, Number ann1, Number ann, Number rn,
75-
I1& lower_begin, const I1& lower_end,
76-
I2& diag_begin, const I2& diag_end,
77-
I3& upper_begin, const I3& upper_end,
78-
I4& right_begin, const I4& right_end
69+
Number a11, Number a12, Number r1,
70+
Number ann1, Number ann, Number rn,
71+
auto lower, auto diag, auto upper, auto right,
72+
SizeT n
7973
) {
80-
const auto n = right_end - right_begin + 2;
81-
assert(n == (lower_end - lower_begin + 2));
82-
assert(n == (diag_end - diag_begin + 2));
83-
assert(n == (upper_end - upper_begin + 2));
84-
auto lower = lower_begin;
85-
auto diag = diag_begin;
86-
auto upper = upper_begin;
87-
auto right = right_begin;
8874
const auto m_1 = *lower / a11;
8975
*diag -= m_1 * a12;
90-
*upper -= m_1 * a13;
9176
*right -= m_1 * r1;
9277
++lower, ++diag, ++upper, ++right;
93-
for (; lower < lower_end; ++lower, ++diag, ++upper, ++right) {
78+
for (SizeT i = 3; i < n; ++i, ++lower, ++diag, ++upper, ++right) {
9479
const auto m = *lower / *(diag - 1);
9580
*diag -= m * *(upper - 1);
9681
*right -= m * *(right - 1);
9782
}
98-
if (ann2 != 0) {
99-
const auto m_n = ann2 / *(lower - 1);
100-
ann1 -= m_n * *(diag - 1);
101-
ann -= m_n * *(upper - 1);
102-
rn -= m_n * *(right - 1);
103-
}
10483
const auto m_n = ann1 / *(diag - 1);
10584
ann -= m_n * *(upper - 1);
10685
rn -= m_n * *(right - 1);
@@ -109,7 +88,7 @@ namespace alfi::util::linalg {
10988
for (auto i = n - 2; i > 0; --i) {
11089
X[i] = (*--right - *--upper * X[i+1]) / *--diag;
11190
}
112-
X[0] = (r1 - a12 * X[1] - a13 * X[2]) / a11;
91+
X[0] = (r1 - a12 * X[1]) / a11;
11392
return X;
11493
}
11594
}

‎ALFI/ALFI/util/spline.h‎

Lines changed: 0 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -100,13 +100,4 @@ namespace alfi::util::spline {
100100
<< ") is too big (bigger than 4). Returning an empty array..." << std::endl;
101101
return {};
102102
}
103-
104-
template <typename Number = DefaultNumber, template <typename> class C1 = DefaultContainer, template <typename> class C2 = DefaultContainer>
105-
void merge_coeffs(auto& coeffs, const C1<C2<Number>>& separate) {
106-
for (SizeT index = 0, i = 0; index < coeffs.size(); ++i) {
107-
for (SizeT k = 0; k < separate.size(); ++k) {
108-
coeffs[index++] = separate[k][i];
109-
}
110-
}
111-
}
112103
}

0 commit comments

Comments
 (0)