Simbody 3.7
Loading...
Searching...
No Matches
GCVSPLUtil.h
Go to the documentation of this file.
1#ifndef SimTK_SIMMATH_GCVSPL_UTIL_H_
2#define SimTK_SIMMATH_GCVSPL_UTIL_H_
3
4/* -------------------------------------------------------------------------- *
5 * Simbody(tm): SimTKmath *
6 * -------------------------------------------------------------------------- *
7 * This is part of the SimTK biosimulation toolkit originating from *
8 * Simbios, the NIH National Center for Physics-Based Simulation of *
9 * Biological Structures at Stanford, funded under the NIH Roadmap for *
10 * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. *
11 * *
12 * Portions copyright (c) 2008-12 Stanford University and the Authors. *
13 * Authors: Peter Eastman *
14 * Contributors: *
15 * *
16 * Licensed under the Apache License, Version 2.0 (the "License"); you may *
17 * not use this file except in compliance with the License. You may obtain a *
18 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
19 * *
20 * Unless required by applicable law or agreed to in writing, software *
21 * distributed under the License is distributed on an "AS IS" BASIS, *
22 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
23 * See the License for the specific language governing permissions and *
24 * limitations under the License. *
25 * -------------------------------------------------------------------------- */
26
27#include "SimTKcommon.h"
29
30// These are global functions.
31int SimTK_gcvspl_(const SimTK::Real *, const SimTK::Real *, int *, const SimTK::Real *, const SimTK::Real *, int *, int *,
32 int *, int *, SimTK::Real *, SimTK::Real *, int *, SimTK::Real *, int *);
33SimTK::Real SimTK_splder_(int *, int *, int *, SimTK::Real *, const SimTK::Real *, const SimTK::Real *, int *, SimTK::Real *, int);
34
35namespace SimTK {
36
45public:
46 static void gcvspl(const Vector& x, const Vector& y, const Vector& wx, Real wy, int m, int md, Real val, Vector& c, Vector& wk, int& ier);
47 template <int K>
48 static void gcvspl(const Vector& x, const Vector_<Vec<K> >& y, const Vector& wx, Vec<K> wy, int m, int md, Real val, Vector_<Vec<K> >& c, Vector& wk, int& ier);
49 static Real splder(int derivOrder, int degree, Real t, const Vector& x, const Vector& coeff);
50 template <int K>
51 static Vec<K> splder(int derivOrder, int degree, Real t, const Vector& x, const Vector_<Vec<K> >& coeff);
52};
53
54template <int K>
55void GCVSPLUtil::gcvspl(const Vector& x, const Vector_<Vec<K> >& y, const Vector& wx, Vec<K> wy, int degree, int md, Real val, Vector_<Vec<K> >& c, Vector& wk, int& ier) {
56 SimTK_APIARGCHECK_ALWAYS(degree > 0 && degree%2==1, "GCVSPLUtil", "gcvspl", "degree must be positive and odd");
57 SimTK_APIARGCHECK_ALWAYS(y.size() >= x.size(), "GCVSPLUtil", "gcvspl", "y is shorter than x");
58 SimTK_APIARGCHECK_ALWAYS(wx.size() >= x.size(), "GCVSPLUtil", "gcvspl", "wx and x must be the same size");
59 SimTK_APIARGCHECK_ALWAYS(x.hasContiguousData(), "GCVSPLUtil", "gcvspl", "x must have contiguous storage (i.e. not be a view)");
60 SimTK_APIARGCHECK_ALWAYS(wk.hasContiguousData(), "GCVSPLUtil", "gcvspl", "wk must have contiguous storage (i.e. not be a view)");
61
62 // Create various temporary variables.
63
64 int m = (degree+1)/2;
65 int n = x.size();
66 int ny = y.size();
67 Vector yvec(ny*K);
68 int index = 0;
69 for (int j = 0; j < K; ++j)
70 for (int i = 0; i < ny; ++i)
71 yvec[index++] = y[i][j];
72 int nc = n*K;
73 Vector cvec(nc);
74 wk.resize(6*(m*n+1)+n);
75 int k = K;
76
77 // Invoke GCV.
78
79 SimTK_gcvspl_(&x[0], &yvec[0], &ny, &wx[0], &wy[0], &m, &n, &k, &md, &val, &cvec[0], &n, &wk[0], &ier);
80 if (ier != 0) {
81 SimTK_APIARGCHECK_ALWAYS(n >= 2*m, "GCVSPLUtil", "gcvspl", "Too few data points");
82 SimTK_APIARGCHECK_ALWAYS(ier != 2, "GCVSPLUtil", "gcvspl", "The values in x must be strictly increasing");
83 SimTK_APIARGCHECK_ALWAYS(ier == 0, "GCVSPLUtil", "gcvspl", "GCVSPL returned an error code");
84 }
85 c.resize(n);
86 index = 0;
87 for (int j = 0; j < K; ++j)
88 for (int i = 0; i < n; ++i)
89 c[i][j] = cvec[index++];
90}
91
92template <int K>
93Vec<K> GCVSPLUtil::splder(int derivOrder, int degree, Real t, const Vector& x, const Vector_<Vec<K> >& coeff) {
94 assert(derivOrder >= 0);
95 assert(t >= x[0] && t <= x[x.size()-1]);
96 assert(x.size() == coeff.size());
97 assert(degree > 0 && degree%2==1);
98 assert(x.hasContiguousData());
99
100 // Create various temporary variables.
101
102
103 Vec<K> result;
104 int m = (degree+1)/2;
105 int n = x.size();
106 int interval = (int) ceil(n*(t-x[0])/(x[n-1]-x[0]));
107
108 const int MaxCheapM = 32;
109 Real qbuf[2*MaxCheapM];
110
111 Real *q = qbuf; // tentatively
112 if (m > MaxCheapM)
113 q = new Real[2*m]; // use heap instead; don't forget to delete
114
115 int offset = (int) (&coeff[1][0]-&coeff[0][0]);
116
117 // Evaluate the spline one component at a time.
118
119 for (int i = 0; i < K; ++i)
120 result[i] = SimTK_splder_(&derivOrder, &m, &n, &t, &x[0], &coeff[0][i], &interval, q, offset);
121
122 if (m > MaxCheapM)
123 delete[] q;
124
125 return result;
126}
127
128} // namespace SimTK
129
130#endif // SimTK_SIMMATH_GCVSPL_UTIL_H_
131
132
#define SimTK_APIARGCHECK_ALWAYS(cond, className, methodName, msg)
Definition ExceptionMacros.h:224
SimTK::Real SimTK_splder_(int *, int *, int *, SimTK::Real *, const SimTK::Real *, const SimTK::Real *, int *, SimTK::Real *, int)
int SimTK_gcvspl_(const SimTK::Real *, const SimTK::Real *, int *, const SimTK::Real *, const SimTK::Real *, int *, int *, int *, int *, SimTK::Real *, SimTK::Real *, int *, SimTK::Real *, int *)
Includes internal headers providing declarations for the basic SimTK Core classes,...
This is the header file that every Simmath compilation unit should include first.
#define SimTK_SIMMATH_EXPORT
Definition SimTKmath/include/simmath/internal/common.h:64
This class provides entry points for using the GCVSPL algorithm in terms of SimTK data types.
Definition GCVSPLUtil.h:44
static void gcvspl(const Vector &x, const Vector &y, const Vector &wx, Real wy, int m, int md, Real val, Vector &c, Vector &wk, int &ier)
static Real splder(int derivOrder, int degree, Real t, const Vector &x, const Vector &coeff)
bool hasContiguousData() const
Definition MatrixBase.h:837
This is a fixed-length column vector designed for no-overhead inline computation.
Definition Vec.h:184
int size() const
Definition VectorBase.h:396
VectorBase & resize(int m)
Definition VectorBase.h:451
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition Assembler.h:37
SimTK_Real Real
This is the default compiled-in floating point type for SimTK, either float or double.
Definition SimTKcommon/include/SimTKcommon/internal/common.h:606