GCC Code Coverage Report


Directory: ./
File: unittest/finite-differences.cpp
Date: 2024-08-27 18:20:05
Exec Total Coverage
Lines: 0 96 0.0%
Branches: 0 386 0.0%

Line Branch Exec Source
1 //
2 // Copyright (c) 2016-2019 CNRS INRIA
3 //
4
5 #include "pinocchio/multibody/model.hpp"
6 #include "pinocchio/multibody/data.hpp"
7 #include "pinocchio/multibody/sample-models.hpp"
8 #include "pinocchio/algorithm/joint-configuration.hpp"
9 #include "pinocchio/algorithm/kinematics.hpp"
10 #include "pinocchio/algorithm/jacobian.hpp"
11
12 #include <iostream>
13 #include <boost/test/unit_test.hpp>
14 #include <boost/utility/binary.hpp>
15
16 using namespace pinocchio;
17 using namespace Eigen;
18
19 template<bool local>
20 Data::Matrix6x finiteDiffJacobian(
21 const Model & model, Data & data, const Eigen::VectorXd & q, const Model::JointIndex joint_id)
22 {
23 Data::Matrix6x res(6, model.nv);
24 res.setZero();
25 VectorXd q_integrate(model.nq);
26 VectorXd v_integrate(model.nv);
27 v_integrate.setZero();
28
29 forwardKinematics(model, data, q);
30 const SE3 oMi_ref = data.oMi[joint_id];
31
32 double eps = 1e-8;
33 for (int k = 0; k < model.nv; ++k)
34 {
35 // Integrate along kth direction
36 v_integrate[k] = eps;
37 q_integrate = integrate(model, q, v_integrate);
38
39 forwardKinematics(model, data, q_integrate);
40 const SE3 & oMi = data.oMi[joint_id];
41
42 if (local)
43 res.col(k) = log6(oMi_ref.inverse() * oMi).toVector();
44 else
45 res.col(k) = oMi_ref.act(log6(oMi_ref.inverse() * oMi)).toVector();
46
47 res.col(k) /= eps;
48
49 v_integrate[k] = 0.;
50 }
51
52 return res;
53 }
54
55 template<typename Matrix>
56 void filterValue(MatrixBase<Matrix> & mat, typename Matrix::Scalar value)
57 {
58 for (int k = 0; k < mat.size(); ++k)
59 mat.derived().data()[k] =
60 math::fabs(mat.derived().data()[k]) <= value ? 0 : mat.derived().data()[k];
61 }
62
63 template<typename JointModel_>
64 struct init;
65
66 template<typename JointModel_>
67 struct init
68 {
69 static JointModel_ run()
70 {
71 JointModel_ jmodel;
72 jmodel.setIndexes(0, 0, 0);
73 return jmodel;
74 }
75 };
76
77 template<typename Scalar, int Options>
78 struct init<pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options>>
79 {
80 typedef pinocchio::JointModelRevoluteUnalignedTpl<Scalar, Options> JointModel;
81
82 static JointModel run()
83 {
84 typedef typename JointModel::Vector3 Vector3;
85 JointModel jmodel(Vector3::Random().normalized());
86
87 jmodel.setIndexes(0, 0, 0);
88 return jmodel;
89 }
90 };
91
92 template<typename Scalar, int Options>
93 struct init<pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options>>
94 {
95 typedef pinocchio::JointModelRevoluteUnboundedUnalignedTpl<Scalar, Options> JointModel;
96
97 static JointModel run()
98 {
99 typedef typename JointModel::Vector3 Vector3;
100 JointModel jmodel(Vector3::Random().normalized());
101
102 jmodel.setIndexes(0, 0, 0);
103 return jmodel;
104 }
105 };
106
107 template<typename Scalar, int Options>
108 struct init<pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options>>
109 {
110 typedef pinocchio::JointModelPrismaticUnalignedTpl<Scalar, Options> JointModel;
111
112 static JointModel run()
113 {
114 typedef typename JointModel::Vector3 Vector3;
115 JointModel jmodel(Vector3::Random().normalized());
116
117 jmodel.setIndexes(0, 0, 0);
118 return jmodel;
119 }
120 };
121
122 template<typename Scalar, int Options>
123 struct init<pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options>>
124 {
125 typedef pinocchio::JointModelHelicalUnalignedTpl<Scalar, Options> JointModel;
126
127 static JointModel run()
128 {
129 typedef typename JointModel::Vector3 Vector3;
130 JointModel jmodel(Vector3::Random().normalized());
131
132 jmodel.setIndexes(0, 0, 0);
133 return jmodel;
134 }
135 };
136
137 template<typename Scalar, int Options>
138 struct init<pinocchio::JointModelUniversalTpl<Scalar, Options>>
139 {
140 typedef pinocchio::JointModelUniversalTpl<Scalar, Options> JointModel;
141
142 static JointModel run()
143 {
144 typedef typename JointModel::Vector3 Vector3;
145 JointModel jmodel(XAxis::vector(), YAxis::vector());
146
147 jmodel.setIndexes(0, 0, 0);
148 return jmodel;
149 }
150 };
151
152 template<typename Scalar, int Options, int axis>
153 struct init<pinocchio::JointModelHelicalTpl<Scalar, Options, axis>>
154 {
155 typedef pinocchio::JointModelHelicalTpl<Scalar, Options, axis> JointModel;
156
157 static JointModel run()
158 {
159 JointModel jmodel(static_cast<Scalar>(0.5));
160
161 jmodel.setIndexes(0, 0, 0);
162 return jmodel;
163 }
164 };
165
166 template<typename Scalar, int Options, template<typename, int> class JointCollection>
167 struct init<pinocchio::JointModelTpl<Scalar, Options, JointCollection>>
168 {
169 typedef pinocchio::JointModelTpl<Scalar, Options, JointCollection> JointModel;
170
171 static JointModel run()
172 {
173 typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
174 JointModel jmodel((JointModelRX()));
175
176 jmodel.setIndexes(0, 0, 0);
177 return jmodel;
178 }
179 };
180
181 template<typename Scalar, int Options, template<typename, int> class JointCollection>
182 struct init<pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection>>
183 {
184 typedef pinocchio::JointModelCompositeTpl<Scalar, Options, JointCollection> JointModel;
185
186 static JointModel run()
187 {
188 typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 0> JointModelRX;
189 typedef pinocchio::JointModelRevoluteTpl<Scalar, Options, 1> JointModelRY;
190 JointModel jmodel((JointModelRX()));
191 jmodel.addJoint(JointModelRY());
192
193 jmodel.setIndexes(0, 0, 0);
194 return jmodel;
195 }
196 };
197
198 template<typename JointModel_>
199 struct init<pinocchio::JointModelMimic<JointModel_>>
200 {
201 typedef pinocchio::JointModelMimic<JointModel_> JointModel;
202
203 static JointModel run()
204 {
205 JointModel_ jmodel_ref = init<JointModel_>::run();
206
207 JointModel jmodel(jmodel_ref, 1., 0.);
208 jmodel.setIndexes(0, 0, 0);
209
210 return jmodel;
211 }
212 };
213
214 struct FiniteDiffJoint
215 {
216 void operator()(JointModelComposite & /*jmodel*/) const
217 {
218 }
219
220 template<typename JointModel>
221 void operator()(JointModelBase<JointModel> & /*jmodel*/) const
222 {
223 typedef typename JointModel::ConfigVector_t CV;
224 typedef typename JointModel::TangentVector_t TV;
225 typedef typename LieGroup<JointModel>::type LieGroupType;
226
227 JointModel jmodel = init<JointModel>::run();
228 std::cout << "name: " << jmodel.classname() << std::endl;
229
230 typename JointModel::JointDataDerived jdata_ = jmodel.createData();
231 typedef JointDataBase<typename JointModel::JointDataDerived> DataBaseType;
232 DataBaseType & jdata = static_cast<DataBaseType &>(jdata_);
233
234 CV q = LieGroupType().random();
235 jmodel.calc(jdata.derived(), q);
236 SE3 M_ref(jdata.M());
237
238 CV q_int(q);
239 const Eigen::DenseIndex nv = jdata.S().nv();
240 TV v(nv);
241 v.setZero();
242 double eps = 1e-8;
243
244 Eigen::Matrix<double, 6, JointModel::NV> S(6, nv), S_ref(jdata.S().matrix());
245
246 for (int k = 0; k < nv; ++k)
247 {
248 v[k] = eps;
249 q_int = LieGroupType().integrate(q, v);
250 jmodel.calc(jdata.derived(), q_int);
251 SE3 M_int = jdata.M();
252
253 S.col(k) = log6(M_ref.inverse() * M_int).toVector();
254 S.col(k) /= eps;
255
256 v[k] = 0.;
257 }
258
259 BOOST_CHECK(S.isApprox(S_ref, eps * 1e1));
260 std::cout << "S_ref:\n" << S_ref << std::endl;
261 std::cout << "S:\n" << S << std::endl;
262 }
263 };
264
265 BOOST_AUTO_TEST_SUITE(BOOST_TEST_MODULE)
266
267 BOOST_AUTO_TEST_CASE(test_S_finit_diff)
268 {
269 boost::mpl::for_each<JointModelVariant::types>(FiniteDiffJoint());
270 }
271
272 BOOST_AUTO_TEST_CASE(test_jacobian_vs_finit_diff)
273 {
274 pinocchio::Model model;
275 pinocchio::buildModels::humanoidRandom(model);
276 pinocchio::Data data(model);
277
278 VectorXd q = VectorXd::Ones(model.nq);
279 q.segment<4>(3).normalize();
280 computeJointJacobians(model, data, q);
281
282 Model::Index idx =
283 model.existJointName("rarm2") ? model.getJointId("rarm2") : (Model::Index)(model.njoints - 1);
284 Data::Matrix6x Jrh(6, model.nv);
285 Jrh.fill(0);
286
287 getJointJacobian(model, data, idx, WORLD, Jrh);
288 Data::Matrix6x Jrh_finite_diff = finiteDiffJacobian<false>(model, data, q, idx);
289 BOOST_CHECK(Jrh_finite_diff.isApprox(Jrh, sqrt(1e-8)));
290
291 getJointJacobian(model, data, idx, LOCAL, Jrh);
292 Jrh_finite_diff = finiteDiffJacobian<true>(model, data, q, idx);
293 BOOST_CHECK(Jrh_finite_diff.isApprox(Jrh, sqrt(1e-8)));
294 }
295
296 BOOST_AUTO_TEST_SUITE_END()
297