Project
Loading...
Searching...
No Matches
Spline1DHelperOld.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
16
17#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
18
19#include "Spline1DHelperOld.h"
20#include "GPUCommonLogger.h"
21#include "TMath.h"
22#include "TMatrixD.h"
23#include "TVectorD.h"
24#include "TDecompBK.h"
25#include <vector>
26
27#include "TRandom.h"
28#include "TMath.h"
29#include "TCanvas.h"
30#include "TNtuple.h"
31#include "TFile.h"
32#include "GPUCommonMath.h"
33#include <iostream>
34
35templateClassImp(o2::gpu::Spline1DHelperOld);
36
37using namespace o2::gpu;
38
39template <typename DataT>
40Spline1DHelperOld<DataT>::Spline1DHelperOld() : mError(), mSpline(), mFdimensions(0)
41{
42}
43
44template <typename DataT>
45int32_t Spline1DHelperOld<DataT>::storeError(int32_t code, const char* msg)
46{
47 mError = msg;
48 return code;
49}
50
51template <typename DataT>
52void Spline1DHelperOld<DataT>::getScoefficients(const typename Spline1D<double>::Knot& knotL, double u,
53 double& cSl, double& cDl, double& cSr, double& cDr)
54{
57
58 // F(u) = cSl * Sl + cSr * Sr + cDl * Dl + cDr * Dr;
59
60 u = u - knotL.u;
61 double v = u * double(knotL.Li); // scaled u
62 double vm1 = v - 1.;
63 double a = u * vm1;
64 double v2 = v * v;
65 cSl = v2 * (2 * v - 3.) + 1; // == 2*v3 - 3*v2 + 1
66 cDl = vm1 * a; // == (v2 - 2v + 1)*u
67 cSr = 1. - cSl;
68 cDr = v * a; // == (v2 - v)*u
69}
70
71template <typename DataT>
72void Spline1DHelperOld<DataT>::getDScoefficients(const typename Spline1D<double>::Knot& knotL, double u,
73 double& cSl, double& cDl, double& cSr, double& cDr)
74{
75 u = u - knotL.u;
76 double dv = double(knotL.Li);
77 double v = u * dv; // scaled u
78 double v2 = v * v;
79 cSl = 6 * (v2 - v) * dv;
80 cDl = 3 * v2 - 4 * v + 1.;
81 cSr = -cSl;
82 cDr = 3 * v2 - 2 * v;
83 // at v==0 : 0, 1, 0, 0
84 // at v==1 : 0, 0, 0, 1
85}
86
87template <typename DataT>
88void Spline1DHelperOld<DataT>::getDDScoefficients(const typename Spline1D<double>::Knot& knotL, double u,
89 double& cSl, double& cDl, double& cSr, double& cDr)
90{
91 u = u - knotL.u;
92 double dv = double(knotL.Li);
93 double v = u * dv; // scaled u
94 cSl = (12 * v - 6) * dv * dv;
95 cDl = (6 * v - 4) * dv;
96 cSr = -cSl;
97 cDr = (6 * v - 2) * dv;
98}
99
100template <typename DataT>
101void Spline1DHelperOld<DataT>::getDDScoefficientsLeft(const typename Spline1D<double>::Knot& knotL,
102 double& cSl, double& cDl, double& cSr, double& cDr)
103{
104 double dv = double(knotL.Li);
105 cSl = -6 * dv * dv;
106 cDl = -4 * dv;
107 cSr = -cSl;
108 cDr = -2 * dv;
109}
110
111template <typename DataT>
112void Spline1DHelperOld<DataT>::getDDScoefficientsRight(const typename Spline1D<double>::Knot& knotL,
113 double& cSl, double& cDl, double& cSr, double& cDr)
114{
115 double dv = double(knotL.Li);
116 cSl = 6 * dv * dv;
117 cDl = 2 * dv;
118 cSr = -cSl;
119 cDr = 4 * dv;
120}
121
122template <typename DataT>
123void Spline1DHelperOld<DataT>::getDDDScoefficients(const typename Spline1D<double>::Knot& knotL,
124 double& cSl, double& cDl, double& cSr, double& cDr)
125{
126 double dv = double(knotL.Li);
127 cSl = 12 * dv * dv * dv;
128 cDl = 6 * dv * dv;
129 cSr = -cSl;
130 cDr = cDl;
131}
132
133template <typename DataT>
134void Spline1DHelperOld<DataT>::approximateFunctionClassic(Spline1DContainer<DataT>& spline,
135 double xMin, double xMax, std::function<void(double x, double f[/*spline.getFdimensions()*/])> F)
136{
139
140 int32_t Ndim = spline.getYdimensions();
141 const int32_t nKnots = spline.getNumberOfKnots();
142
143 TMatrixD A(nKnots, nKnots);
144
145 A.Zero();
146
147 /*
148 const Spline1D::Knot& knot0 = spline.getKnot(i);
149 double x = (u - knot0.u) * knot0.Li; // scaled u
150 double cS1 = (6 - 12*x)*knot0.Li*knot0.Li;
151 double cZ0 = (6*x-4)*knot0.Li;
152 double cZ1 = (6*x-2)*knot0.Li;
153 // f''(u) = cS1*(f1-f0) + cZ0*z0 + cZ1*z1;
154 */
155
156 // second derivative at knot0 is 0
157 {
158 const typename Spline1D<DataT>::Knot& knot0 = spline.getKnot(0);
159 double cZ0 = (-4) * knot0.Li;
160 double cZ1 = (-2) * knot0.Li;
161 // f''(u) = cS1*(f1-f0) + cZ0*z0 + cZ1*z1;
162 A(0, 0) = cZ0;
163 A(0, 1) = cZ1;
164 }
165
166 // second derivative at knot nKnots-1 is 0
167 {
168 const typename Spline1D<DataT>::Knot& knot0 = spline.getKnot(nKnots - 2);
169 double cZ0 = (6 - 4) * knot0.Li;
170 double cZ1 = (6 - 2) * knot0.Li;
171 // f''(u) = cS1*(f1-f0) + cZ0*z0 + cZ1*z1;
172 A(nKnots - 1, nKnots - 2) = cZ0;
173 A(nKnots - 1, nKnots - 1) = cZ1;
174 }
175
176 // second derivative at other knots is same from the left and from the right
177 for (int32_t i = 1; i < nKnots - 1; i++) {
178 const typename Spline1D<DataT>::Knot& knot0 = spline.getKnot(i - 1);
179 double cZ0 = (6 - 4) * knot0.Li;
180 double cZ1_0 = (6 - 2) * knot0.Li;
181 // f''(u) = cS1*(f1-f0) + cZ0*z0 + cZ1*z1;
182
183 const typename Spline1D<DataT>::Knot& knot1 = spline.getKnot(i);
184 double cZ1_1 = (-4) * knot1.Li;
185 double cZ2 = (-2) * knot1.Li;
186 // f''(u) = cS2*(f2-f1) + cZ1_1*z1 + cZ2*z2;
187 A(i, i - 1) = cZ0;
188 A(i, i) = cZ1_0 - cZ1_1;
189 A(i, i + 1) = -cZ2;
190 }
191
192 A.Invert();
193
194 spline.setXrange(xMin, xMax);
195 DataT* parameters = spline.getParameters();
196
197 TVectorD b(nKnots);
198 b.Zero();
199
200 double uToXscale = (((double)xMax) - xMin) / spline.getUmax();
201
202 for (int32_t i = 0; i < nKnots; ++i) {
203 const typename Spline1D<DataT>::Knot& knot = spline.getKnot(i);
204 double u = knot.u;
205 double f[Ndim];
206 F(xMin + u * uToXscale, f);
207 for (int32_t dim = 0; dim < Ndim; dim++) {
208 parameters[(2 * i) * Ndim + dim] = f[dim];
209 }
210 }
211
212 for (int32_t dim = 0; dim < Ndim; dim++) {
213
214 // second derivative at knot0 is 0
215 {
216 double f0 = parameters[(2 * 0) * Ndim + dim];
217 double f1 = parameters[(2 * 1) * Ndim + dim];
218 const typename Spline1D<DataT>::Knot& knot0 = spline.getKnot(0);
219 double cS1 = (6) * knot0.Li * knot0.Li;
220 // f''(u) = cS1*(f1-f0) + cZ0*z0 + cZ1*z1;
221 b(0) = -cS1 * (f1 - f0);
222 }
223
224 // second derivative at knot nKnots-1 is 0
225 {
226 double f0 = parameters[2 * (nKnots - 2) * Ndim + dim];
227 double f1 = parameters[2 * (nKnots - 1) * Ndim + dim];
228 const typename Spline1D<DataT>::Knot& knot0 = spline.getKnot(nKnots - 2);
229 double cS1 = (6 - 12) * knot0.Li * knot0.Li;
230 // f''(u) = cS1*(f1-f0) + cZ0*z0 + cZ1*z1;
231 b(nKnots - 1) = -cS1 * (f1 - f0);
232 }
233
234 // second derivative at other knots is same from the left and from the right
235 for (int32_t i = 1; i < nKnots - 1; i++) {
236 double f0 = parameters[2 * (i - 1) * Ndim + dim];
237 double f1 = parameters[2 * (i)*Ndim + dim];
238 double f2 = parameters[2 * (i + 1) * Ndim + dim];
239 const typename Spline1D<DataT>::Knot& knot0 = spline.getKnot(i - 1);
240 double cS1 = (6 - 12) * knot0.Li * knot0.Li;
241 // f''(u) = cS1*(f1-f0) + cZ0*z0 + cZ1*z1;
242
243 const typename Spline1D<DataT>::Knot& knot1 = spline.getKnot(i);
244 double cS2 = (6) * knot1.Li * knot1.Li;
245 // f''(u) = cS2*(f2-f1) + cZ1_1*z1 + cZ2*z2;
246 b(i) = -cS1 * (f1 - f0) + cS2 * (f2 - f1);
247 }
248
249 TVectorD c = A * b;
250 for (int32_t i = 0; i < nKnots; i++) {
251 parameters[(2 * i + 1) * Ndim + dim] = c[i];
252 }
253 }
254}
255
256template <typename DataT>
257void Spline1DHelperOld<DataT>::approximateDataPoints(
258 Spline1DContainer<DataT>& spline,
259 double xMin, double xMax,
260 double vx[], double vf[], int32_t nDataPoints)
261{
263
264 spline.setXrange(xMin, xMax);
265 setSpline(spline, spline.getYdimensions(), xMin, xMax, vx, nDataPoints);
266 approximateFunction(spline.getParameters(), vf);
267}
268
269template <typename DataT>
270void Spline1DHelperOld<DataT>::approximateFunction(
271 Spline1DContainer<DataT>& spline, double xMin, double xMax, std::function<void(double x, double f[/*spline.getFdimensions()*/])> F,
272 int32_t nAuxiliaryDataPoints)
273{
275 setSpline(spline, spline.getYdimensions(), nAuxiliaryDataPoints);
276 approximateFunction(spline.getParameters(), xMin, xMax, F);
277 spline.setXrange(xMin, xMax);
278}
279
280template <typename DataT>
281void Spline1DHelperOld<DataT>::approximateFunctionGradually(
282 Spline1DContainer<DataT>& spline, double xMin, double xMax, std::function<void(double x, double f[/*spline.getFdimensions()*/])> F,
283 int32_t nAuxiliaryDataPoints)
284{
286 setSpline(spline, spline.getYdimensions(), nAuxiliaryDataPoints);
287 approximateFunctionGradually(spline.getParameters(), xMin, xMax, F);
288 spline.setXrange(xMin, xMax);
289}
290
291template <typename DataT>
292void Spline1DHelperOld<DataT>::approximateFunction(
293 DataT* Sparameters, double xMin, double xMax, std::function<void(double x, double f[])> F) const
294{
297 std::vector<double> vF(getNumberOfDataPoints() * mFdimensions);
298 double mUtoXscale = (((double)xMax) - xMin) / mSpline.getUmax();
299 for (int32_t i = 0; i < getNumberOfDataPoints(); i++) {
300 F(xMin + mUtoXscale * mDataPoints[i].u, &vF[i * mFdimensions]);
301 }
302 approximateFunction(Sparameters, vF.data());
303}
304
305template <typename DataT>
306void Spline1DHelperOld<DataT>::approximateFunctionGradually(
307 DataT* Sparameters, double xMin, double xMax, std::function<void(double x, double f[])> F) const
308{
311 std::vector<double> vF(getNumberOfDataPoints() * mFdimensions);
312 double mUtoXscale = (((double)xMax) - xMin) / mSpline.getUmax();
313 for (int32_t i = 0; i < getNumberOfDataPoints(); i++) {
314 F(xMin + mUtoXscale * mDataPoints[i].u, &vF[i * mFdimensions]);
315 }
316 approximateFunctionGradually(Sparameters, vF.data());
317}
318
319template <typename DataT>
320int32_t Spline1DHelperOld<DataT>::setSpline(
321 const Spline1DContainer<DataT>& spline, int32_t nFdimensions, int32_t nAuxiliaryDataPoints)
322{
323 // Prepare creation of a best-fit spline
324 //
325 // Data points will be set at all integer U (that includes all knots),
326 // plus at nAuxiliaryDataPoints points between the integers.
327 //
328 // nAuxiliaryDataPoints must be >= 2
329 //
330 // nAuxiliaryDataPoints==1 is also possible, but there must be at least
331 // one integer U without a knot, in order to get 2*nKnots data points in total.
332 //
333 // The return value is an error index, 0 means no error
334
335 int32_t ret = 0;
336
337 mFdimensions = nFdimensions;
338 int32_t nPoints = 0;
339 if (!spline.isConstructed()) {
340 ret = storeError(-1, "Spline1DHelperOld<DataT>::setSpline: input spline is not constructed");
341 mSpline.recreate(0, 2);
342 nAuxiliaryDataPoints = 2;
343 nPoints = 4;
344 } else {
345 std::vector<int32_t> knots;
346 for (int32_t i = 0; i < spline.getNumberOfKnots(); i++) {
347 knots.push_back(spline.getKnot(i).getU());
348 }
349 mSpline.recreate(0, spline.getNumberOfKnots(), knots.data());
350
351 nPoints = 1 + mSpline.getUmax() + mSpline.getUmax() * nAuxiliaryDataPoints;
352 if (nPoints < 2 * mSpline.getNumberOfKnots()) {
353 nAuxiliaryDataPoints = 2;
354 nPoints = 1 + mSpline.getUmax() + mSpline.getUmax() * nAuxiliaryDataPoints;
355 ret = storeError(-3, "Spline1DHelperOld::setSpline: too few nAuxiliaryDataPoints, increase to 2");
356 }
357 }
358
359 const int32_t nPar = 2 * mSpline.getNumberOfKnots(); // n parameters for 1D
360
361 mDataPoints.resize(nPoints);
362 double scalePoints2Knots = ((double)mSpline.getUmax()) / (nPoints - 1.);
363
364 for (int32_t i = 0; i < nPoints; ++i) {
365 DataPoint& p = mDataPoints[i];
366 double u = i * scalePoints2Knots;
367 int32_t iKnot = mSpline.getLeftKnotIndexForU(u);
368 const typename Spline1D<double>::Knot& knot0 = mSpline.getKnot(iKnot);
369 const typename Spline1D<double>::Knot& knot1 = mSpline.getKnot(iKnot + 1);
370 double l = knot1.u - knot0.u;
371 double s = (u - knot0.u) * knot0.Li; // scaled u
372 double s2 = s * s;
373 double sm1 = s - 1.;
374
375 p.iKnot = iKnot;
376 p.isKnot = 0;
377 p.u = u;
378 p.cS1 = s2 * (3. - 2. * s);
379 p.cS0 = 1. - p.cS1;
380 p.cZ0 = s * sm1 * sm1 * l;
381 p.cZ1 = s2 * sm1 * l;
382 }
383
384 const int32_t nKnots = mSpline.getNumberOfKnots();
385
386 mKnotDataPoints.resize(nKnots);
387
388 for (int32_t i = 0; i < nKnots; ++i) {
389 const typename Spline1D<double>::Knot& knot = mSpline.getKnot(i);
390 int32_t iu = (int32_t)(knot.u + 0.1f);
391 mKnotDataPoints[i] = iu * (1 + nAuxiliaryDataPoints);
392 mDataPoints[mKnotDataPoints[i]].isKnot = 1;
393 }
394
395 TMatrixDSym A(nPar);
396 A.Zero();
397
398 for (int32_t i = 0; i < nPoints; ++i) {
399 DataPoint& p = mDataPoints[i];
400 int32_t j = p.iKnot * 2;
401 A(j + 0, j + 0) += p.cS0 * p.cS0;
402 A(j + 1, j + 0) += p.cS0 * p.cZ0;
403 A(j + 2, j + 0) += p.cS0 * p.cS1;
404 A(j + 3, j + 0) += p.cS0 * p.cZ1;
405
406 A(j + 1, j + 1) += p.cZ0 * p.cZ0;
407 A(j + 2, j + 1) += p.cZ0 * p.cS1;
408 A(j + 3, j + 1) += p.cZ0 * p.cZ1;
409
410 A(j + 2, j + 2) += p.cS1 * p.cS1;
411 A(j + 3, j + 2) += p.cS1 * p.cZ1;
412
413 A(j + 3, j + 3) += p.cZ1 * p.cZ1;
414 }
415
416 // copy symmetric matrix elements
417
418 for (int32_t i = 0; i < nPar; i++) {
419 for (int32_t j = i + 1; j < nPar; j++) {
420 A(i, j) = A(j, i);
421 }
422 }
423
424 TMatrixDSym Z(nKnots);
425 mLSMmatrixSvalues.resize(nKnots * nKnots);
426 for (int32_t i = 0, k = 0; i < nKnots; i++) {
427 for (int32_t j = 0; j < nKnots; j++, k++) {
428 mLSMmatrixSvalues[k] = A(i * 2 + 1, j * 2);
429 Z(i, j) = A(i * 2 + 1, j * 2 + 1);
430 }
431 }
432
433 {
434 TDecompBK bk(A, 0);
435 bool ok = bk.Invert(A);
436
437 if (!ok) {
438 ret = storeError(-4, "Spline1DHelperOld::setSpline: internal error - can not invert the matrix");
439 A.Zero();
440 }
441 mLSMmatrixFull.resize(nPar * nPar);
442 for (int32_t i = 0, k = 0; i < nPar; i++) {
443 for (int32_t j = 0; j < nPar; j++, k++) {
444 mLSMmatrixFull[k] = A(i, j);
445 }
446 }
447 }
448
449 {
450 TDecompBK bk(Z, 0);
451 if (!bk.Invert(Z)) {
452 ret = storeError(-5, "Spline1DHelperOld::setSpline: internal error - can not invert the matrix");
453 Z.Zero();
454 }
455 mLSMmatrixSderivatives.resize(nKnots * nKnots);
456 for (int32_t i = 0, k = 0; i < nKnots; i++) {
457 for (int32_t j = 0; j < nKnots; j++, k++) {
458 mLSMmatrixSderivatives[k] = Z(i, j);
459 }
460 }
461 }
462
463 return ret;
464}
465
466template <typename DataT>
467int32_t Spline1DHelperOld<DataT>::setSpline(
468 const Spline1DContainer<DataT>& spline, int32_t nFdimensions, double xMin, double xMax, double vx[], int32_t nDataPoints)
469{
470 // Prepare creation of a best-fit spline
471 //
472 // Data points will be set at all integer U (that includes all knots),
473 // plus at nAuxiliaryDataPoints points between the integers.
474 //
475 // nAuxiliaryDataPoints must be >= 2
476 //
477 // nAuxiliaryDataPoints==1 is also possible, but there must be at least
478 // one integer U without a knot, in order to get 2*nKnots data points in total.
479 //
480 // The return value is an error index, 0 means no error
481
482 int32_t ret = 0;
483
484 mFdimensions = nFdimensions;
485 int32_t nPoints = nDataPoints;
486 if (!spline.isConstructed()) {
487 ret = storeError(-1, "Spline1DHelperOld<DataT>::setSpline: input spline is not constructed");
488 mSpline.recreate(0, 2);
489 } else {
490 std::vector<int32_t> knots;
491 for (int32_t i = 0; i < spline.getNumberOfKnots(); i++) {
492 knots.push_back(spline.getKnot(i).getU());
493 }
494 mSpline.recreate(0, spline.getNumberOfKnots(), knots.data());
495 }
496
497 mSpline.setXrange(xMin, xMax);
498
499 const int32_t nPar = 2 * mSpline.getNumberOfKnots(); // n parameters for 1D
500
501 mDataPoints.resize(nPoints);
502
503 for (int32_t i = 0; i < nPoints; ++i) {
504 DataPoint& p = mDataPoints[i];
505 double u = mSpline.convXtoU(vx[i]);
506 int32_t iKnot = mSpline.getLeftKnotIndexForU(u);
507 p.iKnot = iKnot;
508 p.isKnot = 0;
509 p.u = u;
510 const typename Spline1D<double>::Knot& knot0 = mSpline.getKnot(iKnot);
511 getScoefficients(knot0, u, p.cS0, p.cZ0, p.cS1, p.cZ1);
512 }
513
514 const int32_t nKnots = mSpline.getNumberOfKnots();
515
516 TMatrixDSym A(nPar);
517 A.Zero();
518
519 for (int32_t i = 0; i < nPoints; ++i) {
520 DataPoint& p = mDataPoints[i];
521 int32_t j = p.iKnot * 2;
522 A(j + 0, j + 0) += p.cS0 * p.cS0;
523 A(j + 1, j + 0) += p.cS0 * p.cZ0;
524 A(j + 2, j + 0) += p.cS0 * p.cS1;
525 A(j + 3, j + 0) += p.cS0 * p.cZ1;
526
527 A(j + 1, j + 1) += p.cZ0 * p.cZ0;
528 A(j + 2, j + 1) += p.cZ0 * p.cS1;
529 A(j + 3, j + 1) += p.cZ0 * p.cZ1;
530
531 A(j + 2, j + 2) += p.cS1 * p.cS1;
532 A(j + 3, j + 2) += p.cS1 * p.cZ1;
533
534 A(j + 3, j + 3) += p.cZ1 * p.cZ1;
535 }
536
537 for (int32_t iKnot = 0; iKnot < nKnots - 2; ++iKnot) {
538 const typename Spline1D<double>::Knot& knot0 = mSpline.getKnot(iKnot);
539 const typename Spline1D<double>::Knot& knot1 = mSpline.getKnot(iKnot + 1);
540 // const typename Spline1D<double>::Knot& knot2 = mSpline.getKnot(iKnot + 2);
541 /*
542 another way to calculate f(u):
543 T uu = T(u - knotL.u);
544 T v = uu * T(knotL.Li); // scaled u
545 T vm1 = v-1;
546 T v2 = v * v;
547 float cSr = 3*v2 - 2*v3;
548 float cSl = 1-cSr;
549 float cDl = (v3-2*v2+v)*knotL.L;
550 float cDr = (v3-v2)*knotL.L;
551 return cSl*Sl + cSr*Sr + cDl*Dl + cDr*Dr;
552 ()'v:
553 aSr = 6*v - 6*v2
554 aSl = -aSr
555 aDl = (3*v2-4*v+1)*knotL.L;
556 aDr = (3*v2-2*v)*knotL.L;
557 ()''v
558 bSr = 6 -12*v
559 bSl = -bSr
560 bDl = (6*v-4)*knotL.L;
561 bDr = (6*v-2)*knotL.L;
562 ()''u
563 dSr = (6 - 12*v)*knotL.Li*knotL.Li;
564 dSl = -dSr
565 dDl = (6*v-4)*knotL.Li;
566 dDr = (6*v-2)*knotL.Li;
567 */
568
569 double l0 = knot0.Li;
570 double l1 = knot1.Li;
571 /*
572 v1 = +6*l0*l0*s0 + 2*l0*z0 - 6*l0*l0*s1 + 4*l0*z1 ;
573 v2 = -6*l1*l1*s1 - 4*l1*z1 + 6*l1*l1*s2 - 2*l1*z2 ;
574 v2-v1 = -6*l1*l1*s1 - 4*l1*z1 + 6*l1*l1*s2 - 2*l1*z2 -6*l0*l0*s0 - 2*l0*z0 + 6*l0*l0*s1 - 4*l0*z1
575 = -6*l0*l0*s0 - 2*l0*z0 -6*(l1*l1-l0*l0)*s1 - 4*(l0+l1)*z1 + 6*l1*l1*s2 - 2*l1*z2
576 */
577 double c = 0.01;
578 double cS0 = c * -3 * l0 * l0;
579 double cZ0 = c * -l0;
580 double cS1 = c * -3 * (l1 * l1 - l0 * l0);
581 double cZ1 = c * -2 * (l0 + l1);
582 double cS2 = c * 3 * l1 * l1;
583 double cZ2 = c * -l1;
584
585 int32_t j = iKnot * 2;
586
587 A(j + 0, j + 0) += cS0 * cS0;
588 A(j + 1, j + 0) += cS0 * cZ0;
589 A(j + 2, j + 0) += cS0 * cS1;
590 A(j + 3, j + 0) += cS0 * cZ1;
591 A(j + 4, j + 0) += cS0 * cS2;
592 A(j + 5, j + 0) += cS0 * cZ2;
593
594 A(j + 1, j + 1) += cZ0 * cZ0;
595 A(j + 2, j + 1) += cZ0 * cS1;
596 A(j + 3, j + 1) += cZ0 * cZ1;
597 A(j + 4, j + 1) += cZ0 * cS2;
598 A(j + 5, j + 1) += cZ0 * cZ2;
599
600 A(j + 2, j + 2) += cS1 * cS1;
601 A(j + 3, j + 2) += cS1 * cZ1;
602 A(j + 4, j + 2) += cS1 * cS2;
603 A(j + 5, j + 2) += cS1 * cZ2;
604
605 A(j + 3, j + 3) += cZ1 * cZ1;
606 A(j + 4, j + 3) += cZ1 * cS2;
607 A(j + 5, j + 3) += cZ1 * cZ2;
608
609 A(j + 4, j + 4) += cS2 * cS2;
610 A(j + 5, j + 4) += cS2 * cZ2;
611
612 A(j + 5, j + 5) += cZ2 * cZ2;
613 }
614
615 for (int32_t iKnot = -1; iKnot < nKnots - 2; ++iKnot) {
616
617 const typename Spline1D<double>::Knot& knot1 = mSpline.getKnot(iKnot + 1);
618 /*
619 ()''u
620 dSr = (3 - 6*v)*knotL.Li*knotL.Li;
621 dSl = -dSr
622 dDl = (3*v-2)*knotL.Li;
623 dDr = (3*v-1)*knotL.Li;
624 */
625
626 double l1 = knot1.Li;
627 /*
628 v2 = -3*l1*l1*s1 - 2*l1*z1 + 3*l1*l1*s2 - l1*z2 ;
629 */
630 double c = 0.01;
631 double cS1 = c * -3 * (l1 * l1);
632 double cZ1 = c * -2 * (l1);
633 double cS2 = c * 3 * l1 * l1;
634 double cZ2 = c * -l1;
635
636 int32_t j = iKnot * 2;
637
638 A(j + 2, j + 2) += cS1 * cS1;
639 A(j + 3, j + 2) += cS1 * cZ1;
640 A(j + 4, j + 2) += cS1 * cS2;
641 A(j + 5, j + 2) += cS1 * cZ2;
642
643 A(j + 3, j + 3) += cZ1 * cZ1;
644 A(j + 4, j + 3) += cZ1 * cS2;
645 A(j + 5, j + 3) += cZ1 * cZ2;
646
647 A(j + 4, j + 4) += cS2 * cS2;
648 A(j + 5, j + 4) += cS2 * cZ2;
649
650 A(j + 5, j + 5) += cZ2 * cZ2;
651 }
652
653 {
654 int32_t iKnot = nKnots - 2;
655 const typename Spline1D<double>::Knot& knot0 = mSpline.getKnot(iKnot);
656 /*
657 ()''u
658 dSr = (3 - 6*v)*knotL.Li*knotL.Li;
659 dSl = -dSr
660 dDl = (3*v-2)*knotL.Li;
661 dDr = (3*v-1)*knotL.Li;
662 */
663
664 double l0 = knot0.Li;
665 /*
666 v1 = +3*l0*l0*s0 + l0*z0 - 3*l0*l0*s1 + 2*l0*z1 ;
667 */
668 double c = 0.01;
669 double cS0 = c * 3 * l0 * l0;
670 double cZ0 = c * l0;
671 double cS1 = c * -3 * l0 * l0;
672 double cZ1 = c * 2 * l0;
673
674 int32_t j = iKnot * 2;
675
676 A(j + 0, j + 0) += cS0 * cS0;
677 A(j + 1, j + 0) += cS0 * cZ0;
678 A(j + 2, j + 0) += cS0 * cS1;
679 A(j + 3, j + 0) += cS0 * cZ1;
680
681 A(j + 1, j + 1) += cZ0 * cZ0;
682 A(j + 2, j + 1) += cZ0 * cS1;
683 A(j + 3, j + 1) += cZ0 * cZ1;
684
685 A(j + 2, j + 2) += cS1 * cS1;
686 A(j + 3, j + 2) += cS1 * cZ1;
687
688 A(j + 3, j + 3) += cZ1 * cZ1;
689 }
690
691 // copy symmetric matrix elements
692
693 for (int32_t i = 0; i < nPar; i++) {
694 for (int32_t j = i + 1; j < nPar; j++) {
695 A(i, j) = A(j, i);
696 }
697 }
698
699 TMatrixDSym Z(nKnots);
700 mLSMmatrixSvalues.resize(nKnots * nKnots);
701 for (int32_t i = 0, k = 0; i < nKnots; i++) {
702 for (int32_t j = 0; j < nKnots; j++, k++) {
703 mLSMmatrixSvalues[k] = A(i * 2 + 1, j * 2);
704 Z(i, j) = A(i * 2 + 1, j * 2 + 1);
705 }
706 }
707
708 {
709 TDecompBK bk(A, 0);
710 bool ok = bk.Invert(A);
711
712 if (!ok) {
713 ret = storeError(-4, "Spline1DHelperOld::setSpline: internal error - can not invert the matrix");
714 A.Zero();
715 }
716 mLSMmatrixFull.resize(nPar * nPar);
717 for (int32_t i = 0, k = 0; i < nPar; i++) {
718 for (int32_t j = 0; j < nPar; j++, k++) {
719 mLSMmatrixFull[k] = A(i, j);
720 }
721 }
722 }
723
724 {
725 TDecompBK bk(Z, 0);
726 if (!bk.Invert(Z)) {
727 ret = storeError(-5, "Spline1DHelperOld::setSpline: internal error - can not invert the matrix");
728 Z.Zero();
729 }
730 mLSMmatrixSderivatives.resize(nKnots * nKnots);
731 for (int32_t i = 0, k = 0; i < nKnots; i++) {
732 for (int32_t j = 0; j < nKnots; j++, k++) {
733 mLSMmatrixSderivatives[k] = Z(i, j);
734 }
735 }
736 }
737
738 return ret;
739}
740
741template <typename DataT>
742void Spline1DHelperOld<DataT>::approximateFunction(
743 DataT* Sparameters,
744 const double DataPointF[/*getNumberOfDataPoints() x nFdim*/]) const
745{
747
748 const int32_t nPar = 2 * mSpline.getNumberOfKnots();
749 double b[nPar];
750 for (int32_t idim = 0; idim < mFdimensions; idim++) {
751 for (int32_t i = 0; i < nPar; i++) {
752 b[i] = 0.;
753 }
754 for (int32_t i = 0; i < getNumberOfDataPoints(); ++i) {
755 const DataPoint& p = mDataPoints[i];
756 double* bb = &(b[p.iKnot * 2]);
757 double f = (double)DataPointF[i * mFdimensions + idim];
758 bb[0] += f * p.cS0;
759 bb[1] += f * p.cZ0;
760 bb[2] += f * p.cS1;
761 bb[3] += f * p.cZ1;
762 }
763
764 const double* row = mLSMmatrixFull.data();
765
766 for (int32_t i = 0; i < nPar; i++, row += nPar) {
767 double s = 0.;
768 for (int32_t j = 0; j < nPar; j++) {
769 s += row[j] * b[j];
770 }
771 Sparameters[i * mFdimensions + idim] = (float)s;
772 }
773 }
774}
775
776template <typename DataT>
777void Spline1DHelperOld<DataT>::approximateFunctionGradually(
778 DataT* Sparameters, const double DataPointF[/*getNumberOfDataPoints() x nFdim */]) const
779{
782 copySfromDataPoints(Sparameters, DataPointF);
783 approximateDerivatives(Sparameters, DataPointF);
784}
785
786template <typename DataT>
787void Spline1DHelperOld<DataT>::copySfromDataPoints(
788 DataT* Sparameters, const double DataPointF[/*getNumberOfDataPoints() x nFdim*/]) const
789{
792 for (int32_t i = 0; i < mSpline.getNumberOfKnots(); ++i) { // set F values at knots
793 int32_t ip = mKnotDataPoints[i];
794 for (int32_t d = 0; d < mFdimensions; d++) {
795 Sparameters[2 * i * mFdimensions + d] = DataPointF[ip * mFdimensions + d];
796 }
797 }
798}
799
800template <typename DataT>
801void Spline1DHelperOld<DataT>::approximateDerivatives(
802 DataT* Sparameters, const double DataPointF[/*getNumberOfDataPoints() x nFdim*/]) const
803{
807
808 const int32_t nKnots = mSpline.getNumberOfKnots();
809 const int32_t Ndim = mFdimensions;
810 double b[nKnots * mFdimensions];
811 for (int32_t i = 0; i < nKnots * Ndim; i++) {
812 b[i] = 0.;
813 }
814
815 for (int32_t i = 0; i < getNumberOfDataPoints(); ++i) {
816 const DataPoint& p = mDataPoints[i];
817 if (p.isKnot) {
818 continue;
819 }
820 for (int32_t d = 0; d < Ndim; d++) {
821 double f = (double)DataPointF[i * Ndim + d];
822 b[(p.iKnot + 0) * Ndim + d] += f * p.cZ0;
823 b[(p.iKnot + 1) * Ndim + d] += f * p.cZ1;
824 }
825 }
826
827 const double* row = mLSMmatrixSvalues.data();
828 for (int32_t i = 0; i < nKnots; ++i, row += nKnots) {
829 double s[Ndim];
830 for (int32_t d = 0; d < Ndim; d++) {
831 s[d] = 0.;
832 }
833 for (int32_t j = 0; j < nKnots; ++j) {
834 for (int32_t d = 0; d < Ndim; d++) {
835 s[d] += row[j] * Sparameters[2 * j * Ndim + d];
836 }
837 }
838 for (int32_t d = 0; d < Ndim; d++) {
839 b[i * Ndim + d] -= s[d];
840 }
841 }
842
843 row = mLSMmatrixSderivatives.data();
844 for (int32_t i = 0; i < nKnots; ++i, row += nKnots) {
845 double s[Ndim];
846 for (int32_t d = 0; d < Ndim; d++) {
847 s[d] = 0.;
848 }
849 for (int32_t j = 0; j < nKnots; ++j) {
850 for (int32_t d = 0; d < Ndim; d++) {
851 s[d] += row[j] * b[j * Ndim + d];
852 }
853 }
854 for (int32_t d = 0; d < Ndim; d++) {
855 Sparameters[(2 * i + 1) * Ndim + d] = (float)s[d];
856 }
857 }
858}
859
860template <typename DataT>
861int32_t Spline1DHelperOld<DataT>::test(const bool draw, const bool drawDataPoints)
862{
863 using namespace std;
864
865 // input function F
866
867 const int32_t Ndim = 5;
868 const int32_t Fdegree = 4;
869 double Fcoeff[Ndim][2 * (Fdegree + 1)];
870
871 auto F = [&](double x, double f[]) -> void {
872 double cosx[Fdegree + 1], sinx[Fdegree + 1];
873 double xi = 0;
874 for (int32_t i = 0; i <= Fdegree; i++, xi += x) {
875 GPUCommonMath::SinCosd(xi, sinx[i], cosx[i]);
876 }
877 for (int32_t dim = 0; dim < Ndim; dim++) {
878 f[dim] = 0; // Fcoeff[0]/2;
879 for (int32_t i = 1; i <= Fdegree; i++) {
880 f[dim] += Fcoeff[dim][2 * i] * cosx[i] + Fcoeff[dim][2 * i + 1] * sinx[i];
881 }
882 }
883 };
884
885 TCanvas* canv = nullptr;
886 TNtuple* nt = nullptr;
887 TNtuple* knots = nullptr;
888
889 auto ask = [&]() -> bool {
890 if (!canv) {
891 return 0;
892 }
893 canv->Update();
894 LOG(info) << "type 'q ' to exit";
895 std::string str;
896 std::getline(std::cin, str);
897 return (str != "q" && str != ".q");
898 };
899
900 LOG(info) << "Test 1D interpolation with the compact spline";
901
902 int32_t nTries = 100;
903
904 if (draw) {
905 canv = new TCanvas("cQA", "Spline1D QA", 1000, 600);
906 nTries = 10000;
907 }
908
909 double statDf1 = 0;
910 double statDf2 = 0;
911 double statDf1D = 0;
912 double statN = 0;
913
914 int32_t seed = 1;
915
916 for (int32_t itry = 0; itry < nTries; itry++) {
917
918 // init random F
919 for (int32_t dim = 0; dim < Ndim; dim++) {
920 gRandom->SetSeed(seed++);
921 for (int32_t i = 0; i < 2 * (Fdegree + 1); i++) {
922 Fcoeff[dim][i] = gRandom->Uniform(-1, 1);
923 }
924 }
925
926 // spline
927
928 int32_t nKnots = 4;
929 const int32_t uMax = nKnots * 3;
930
931 Spline1D<DataT, Ndim> spline1;
932 int32_t knotsU[nKnots];
933
934 do { // set knots randomly
935 knotsU[0] = 0;
936 double du = 1. * uMax / (nKnots - 1);
937 for (int32_t i = 1; i < nKnots; i++) {
938 knotsU[i] = (int32_t)(i * du); // + gRandom->Uniform(-du / 3, du / 3);
939 }
940 knotsU[nKnots - 1] = uMax;
941 spline1.recreate(nKnots, knotsU);
942 if (nKnots != spline1.getNumberOfKnots()) {
943 LOG(info) << "warning: n knots changed during the initialisation " << nKnots
944 << " -> " << spline1.getNumberOfKnots();
945 continue;
946 }
947 } while (0);
948
949 std::string err = FlatObject::stressTest(spline1);
950 if (!err.empty()) {
951 LOG(info) << "error at FlatObject functionality: " << err;
952 return -1;
953 } else {
954 // LOG(info) << "flat object functionality is ok" ;
955 }
956
957 nKnots = spline1.getNumberOfKnots();
958 int32_t nAuxiliaryPoints = 1;
959 Spline1D<DataT, Ndim> spline2(spline1);
960 spline1.approximateFunction(0., TMath::Pi(), F, nAuxiliaryPoints);
961
962 // if (itry == 0)
963 {
964 TFile outf("testSpline1D.root", "recreate");
965 if (outf.IsZombie()) {
966 LOG(info) << "Failed to open output file testSpline1D.root ";
967 } else {
968 const char* name = "Spline1Dtest";
969 spline1.writeToFile(outf, name);
970 Spline1D<DataT, Ndim>* p = spline1.readFromFile(outf, name);
971 if (p == nullptr) {
972 LOG(info) << "Failed to read Spline1D from file testSpline1D.root ";
973 } else {
974 spline1 = *p;
975 }
976 outf.Close();
977 }
978 }
979
980 Spline1DHelperOld<DataT> helper;
981 helper.setSpline(spline2, Ndim, nAuxiliaryPoints);
982 helper.approximateFunctionGradually(spline2, 0., TMath::Pi(), F, nAuxiliaryPoints);
983
984 // 1-D splines for each dimension
985 Spline1D<DataT, 1> splines3[Ndim];
986 {
987 for (int32_t dim = 0; dim < Ndim; dim++) {
988 auto F3 = [&](double u, double f[]) -> void {
989 double ff[Ndim];
990 F(u, ff);
991 f[0] = ff[dim];
992 };
993 splines3[dim].recreate(nKnots, knotsU);
994 splines3[dim].approximateFunction(0., TMath::Pi(), F3, nAuxiliaryPoints);
995 }
996 }
997
998 double stepX = 1.e-2;
999 for (double x = 0; x < TMath::Pi(); x += stepX) {
1000 double f[Ndim];
1001 DataT s1[Ndim], s2[Ndim];
1002 F(x, f);
1003 spline1.interpolate(x, s1);
1004 spline2.interpolate(x, s2);
1005 for (int32_t dim = 0; dim < Ndim; dim++) {
1006 statDf1 += (s1[dim] - f[dim]) * (s1[dim] - f[dim]);
1007 statDf2 += (s2[dim] - f[dim]) * (s2[dim] - f[dim]);
1008 DataT s1D = splines3[dim].interpolate(x);
1009 statDf1D += (s1D - s1[dim]) * (s1D - s1[dim]);
1010 }
1011 statN += Ndim;
1012 }
1013 // LOG(info) << "std dev : " << sqrt(statDf1 / statN) ;
1014
1015 if (draw) {
1016 delete nt;
1017 delete knots;
1018 nt = new TNtuple("nt", "nt", "u:f:s");
1019 double drawMax = -1.e20;
1020 double drawMin = 1.e20;
1021 double stepX = 1.e-4;
1022 for (double x = 0; x < TMath::Pi(); x += stepX) {
1023 double f[Ndim];
1024 DataT s[Ndim];
1025 F(x, f);
1026 spline1.interpolate(x, s);
1027 nt->Fill(spline1.convXtoU(x), f[0], s[0]);
1028 drawMax = std::max(drawMax, std::max(f[0], (double)s[0]));
1029 drawMin = std::min(drawMin, std::min(f[0], (double)s[0]));
1030 }
1031
1032 nt->SetMarkerStyle(8);
1033
1034 {
1035 TNtuple* ntRange = new TNtuple("ntRange", "nt", "u:f");
1036 drawMin -= 0.1 * (drawMax - drawMin);
1037
1038 ntRange->Fill(0, drawMin);
1039 ntRange->Fill(0, drawMax);
1040 ntRange->Fill(uMax, drawMin);
1041 ntRange->Fill(uMax, drawMax);
1042 ntRange->SetMarkerColor(kWhite);
1043 ntRange->SetMarkerSize(0.1);
1044 ntRange->Draw("f:u", "", "");
1045 delete ntRange;
1046 }
1047
1048 nt->SetMarkerColor(kGray);
1049 nt->SetMarkerSize(2.);
1050 nt->Draw("f:u", "", "P,same");
1051
1052 nt->SetMarkerSize(.5);
1053 nt->SetMarkerColor(kBlue);
1054 nt->Draw("s:u", "", "P,same");
1055
1056 knots = new TNtuple("knots", "knots", "type:u:s");
1057 for (int32_t i = 0; i < nKnots; i++) {
1058 double u = spline1.getKnot(i).u;
1059 DataT s[Ndim];
1060 spline1.interpolate(spline1.convUtoX(u), s);
1061 knots->Fill(1, u, s[0]);
1062 }
1063
1064 knots->SetMarkerStyle(8);
1065 knots->SetMarkerSize(1.5);
1066 knots->SetMarkerColor(kRed);
1067 knots->SetMarkerSize(1.5);
1068 knots->Draw("s:u", "type==1", "same"); // knots
1069
1070 if (drawDataPoints) {
1071 for (int32_t j = 0; j < helper.getNumberOfDataPoints(); j++) {
1072 const typename Spline1DHelperOld<DataT>::DataPoint& p = helper.getDataPoint(j);
1073 if (p.isKnot) {
1074 continue;
1075 }
1076 DataT s[Ndim];
1077 spline1.interpolate(spline1.convUtoX(p.u), s);
1078 knots->Fill(2, p.u, s[0]);
1079 }
1080 knots->SetMarkerColor(kBlack);
1081 knots->SetMarkerSize(1.);
1082 knots->Draw("s:u", "type==2", "same"); // data points
1083 }
1084
1085 if (!ask()) {
1086 break;
1087 }
1088 } // draw
1089 }
1090 // delete canv;
1091 // delete nt;
1092 // delete knots;
1093
1094 statDf1 = sqrt(statDf1 / statN);
1095 statDf2 = sqrt(statDf2 / statN);
1096 statDf1D = sqrt(statDf1D / statN);
1097
1098 LOG(info) << "\n std dev for Compact Spline : " << statDf1 << " / " << statDf2;
1099 LOG(info) << " mean difference between 1-D and " << Ndim
1100 << "-D splines : " << statDf1D;
1101
1102 if (statDf1 < 0.05 && statDf2 < 0.06 && statDf1D < 1.e-20) {
1103 LOG(info) << "Everything is fine";
1104 } else {
1105 LOG(info) << "Something is wrong!!";
1106 return -2;
1107 }
1108 return 0;
1109}
1110
1111template class o2::gpu::Spline1DHelperOld<float>;
1112template class o2::gpu::Spline1DHelperOld<double>;
1113
1114#endif
i
int32_t i
Definition GPUCommonAlgorithm.h:443
bb
const int16_t bb
Definition GPUTPCGMMerger.cxx:1461
j
uint32_t j
Definition RawData.h:0
c
uint32_t c
Definition RawData.h:2
templateClassImp
templateClassImp(o2::gpu::Spline1DHelperOld)
Spline1DHelperOld.h
Definition of Spline1DHelperOld class.
A
Definition A.h:16
o2::gpu::FlatObject::isConstructed
bool isConstructed() const
Tells if the object is constructed.
Definition FlatObject.h:262
o2::gpu::FlatObject::stressTest
static std::string stressTest(T &obj)
Test the flat object functionality for a child class T.
Definition FlatObject.h:411
o2::gpu::Spline1DHelperOld::getDDScoefficientsLeft
static void getDDScoefficientsLeft(const typename Spline1D< double >::Knot &knotL, double &cSl, double &cDl, double &cSr, double &cDr)
Definition Spline1DHelperOld.cxx:101
o2::gpu::Spline1DHelperOld::Spline1DHelperOld
Spline1DHelperOld()
_____________ Constructors / destructors __________________________
Definition Spline1DHelperOld.cxx:40
o2::gpu::Spline1DHelperOld::approximateFunction
void approximateFunction(Spline1DContainer< DataT > &spline, double xMin, double xMax, std::function< void(double x, double f[])> F, int32_t nAuxiliaryDataPoints=4)
Create best-fit spline parameters for a given input function F.
Definition Spline1DHelperOld.cxx:270
o2::gpu::Spline1DHelperOld::getScoefficients
static void getScoefficients(const typename Spline1D< double >::Knot &knotL, double u, double &cSl, double &cDl, double &cSr, double &cDr)
Definition Spline1DHelperOld.cxx:52
o2::gpu::Spline1DHelperOld::getDScoefficients
static void getDScoefficients(const typename Spline1D< double >::Knot &knotL, double u, double &cSl, double &cDl, double &cSr, double &cDr)
Definition Spline1DHelperOld.cxx:72
o2::gpu::Spline1DHelperOld::test
static int32_t test(const bool draw=0, const bool drawDataPoints=1)
Test the Spline1D class functionality.
Definition Spline1DHelperOld.cxx:861
o2::gpu::Spline1DHelperOld::getDDScoefficientsRight
static void getDDScoefficientsRight(const typename Spline1D< double >::Knot &knotL, double &cSl, double &cDl, double &cSr, double &cDr)
Definition Spline1DHelperOld.cxx:112
o2::gpu::Spline1DHelperOld::approximateDataPoints
void approximateDataPoints(Spline1DContainer< DataT > &spline, double xMin, double xMax, double x[], double f[], int32_t nDataPoints)
Create best-fit spline parameters for a given input function F.
Definition Spline1DHelperOld.cxx:257
o2::gpu::Spline1DHelperOld::approximateFunctionGradually
void approximateFunctionGradually(Spline1DContainer< DataT > &spline, double xMin, double xMax, std::function< void(double x, double f[])> F, int32_t nAuxiliaryDataPoints=4)
Create best-fit spline parameters gradually for a given input function F.
Definition Spline1DHelperOld.cxx:281
o2::gpu::Spline1DHelperOld::copySfromDataPoints
void copySfromDataPoints(DataT *Fparameters, const double DataPointF[]) const
a tool for the gradual approximation: set spline values S_i at knots == function values
Definition Spline1DHelperOld.cxx:787
o2::gpu::Spline1DHelperOld::approximateFunctionClassic
void approximateFunctionClassic(Spline1DContainer< DataT > &spline, double xMin, double xMax, std::function< void(double x, double f[])> F)
Create classic spline parameters for a given input function F.
Definition Spline1DHelperOld.cxx:134
o2::gpu::Spline1DHelperOld::getDDScoefficients
static void getDDScoefficients(const typename Spline1D< double >::Knot &knotL, double u, double &cSl, double &cDl, double &cSr, double &cDr)
Definition Spline1DHelperOld.cxx:88
o2::gpu::Spline1DHelperOld::getDDDScoefficients
static void getDDDScoefficients(const typename Spline1D< double >::Knot &knotL, double &cSl, double &cDl, double &cSr, double &cDr)
Definition Spline1DHelperOld.cxx:123
o2::gpu::Spline1DHelperOld::setSpline
int32_t setSpline(const Spline1DContainer< DataT > &spline, int32_t nFdimensions, int32_t nAuxiliaryDataPoints)
_______________ Interface for a step-wise construction of the best-fit spline _______________________...
Definition Spline1DHelperOld.cxx:320
o2::gpu::Spline1DHelperOld::approximateDerivatives
void approximateDerivatives(DataT *Fparameters, const double DataPointF[]) const
Definition Spline1DHelperOld.cxx:801
o2::gpu::Spline1DHelperOld::getNumberOfDataPoints
int32_t getNumberOfDataPoints() const
number of data points
Definition Spline1DHelperOld.h:105
o2::gpu::Spline1DHelperOld::getDataPoint
const DataPoint & getDataPoint(int32_t ip) const
Definition Spline1DHelperOld.h:126
o2::gpu::Spline1D::readFromFile
static Spline1D * readFromFile(TFile &inpf, const char *name)
read a class object from the file
Definition Spline1D.h:170
x
GLint GLenum GLint x
Definition glcorearb.h:403
v
const GLdouble * v
Definition glcorearb.h:832
s1
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s1
Definition glcorearb.h:5034
name
GLuint const GLchar * name
Definition glcorearb.h:781
f
GLdouble f
Definition glcorearb.h:310
b
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
a
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
v2
GLfloat GLfloat GLfloat v2
Definition glcorearb.h:813
std
Defining DataPointCompositeObject explicitly as copiable.
Definition CcdbObjectInfo.h:121
o2::gpu::Spline1DContainer::Knot::u
DataT u
u coordinate of the knot i (an integer number in float format)
Definition Spline1DSpec.h:55
o2::gpu::Spline1DHelperOld::DataPoint
Helper structure for 1D spline construction.
Definition Spline1DHelperOld.h:43
o2::gpu::Spline1DHelperOld::DataPoint::isKnot
bool isKnot
is the point placed at a knot
Definition Spline1DHelperOld.h:50
o2::gpu::Spline1DHelperOld::DataPoint::iKnot
int32_t iKnot
index of the left knot of the segment
Definition Spline1DHelperOld.h:49
LOG
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
row
std::vector< int > row
Definition test_ctf_io_itsmft.cxx:48
str
const std::string str
Definition test_ransEncodeDecode.cxx:37
msg
uint64_t const void const *restrict const msg
Definition x9.h:153