10 #define VP_DEBUG_MODE 45
13 #include <sot/core/solver-hierarchical-inequalities.hh>
34 SolverHierarchicalInequalities::ConstraintList constraintH;
35 std::vector<bubMatrix> Jes;
36 std::vector<bubVector> ees;
37 std::vector<bubMatrix> Jis;
38 std::vector<bubVector> eiInfs;
39 std::vector<bubVector> eiSups;
40 std::vector<ConstraintMem::BoundSideVector>
bounds;
45 bubVector ee, eiInf, eiSup;
46 ConstraintMem::BoundSideVector eiBoundSide;
49 if (bs !=
"variable") {
50 cerr <<
"!! '" << bs <<
"'" << endl;
55 cerr <<
"!! '" << bs <<
"'" << endl;
60 sotRotationComposedInExtenso Qh(nJ);
61 std::deque<SolverHierarchicalInequalities *> solvers;
63 for (
unsigned int level = 0;; ++level) {
68 }
else if (bs !=
"level") {
69 cerr <<
"!! '" << bs <<
"'" << endl;
73 if (bs !=
"equalities") {
74 cerr <<
"!! '" << bs <<
"'" << endl;
81 for (
int i = 0;
i < me; ++
i) {
82 for (
int j = 0; j < nJ; ++j) off >> Je(
i, j);
88 if (bs !=
"inequalities") {
89 cerr <<
"!! '" << bs <<
"'" << endl;
96 eiBoundSide.resize(mi);
98 for (
int i = 0;
i < mi; ++
i) {
99 for (
int j = 0; j < nJ; ++j) off >> Ji(
i, j);
101 eiBoundSide[
i] = ConstraintMem::BOUND_VOID;
104 eiBoundSide[
i] = (ConstraintMem::BoundSideType)(
105 eiBoundSide[
i] | ConstraintMem::BOUND_INF);
106 eiInf(
i) = atof(number.c_str());
112 eiBoundSide[
i] = (ConstraintMem::BoundSideType)(
113 eiBoundSide[
i] | ConstraintMem::BOUND_SUP);
114 eiSup(
i) = atof(number.c_str());
120 struct timeval t0, t1;
123 <<
"/* ----------------------------------------------------- */"
126 <<
"/* ----------------------------------------------------- */"
129 <<
"/* ----------------------------------------------------- */"
131 sotDEBUG(5) <<
"ee" << level <<
" = " << (MATLAB)ee << endl;
132 sotDEBUG(5) <<
"eiInf" << level <<
" = " << (MATLAB)eiInf << endl;
133 sotDEBUG(5) <<
"eiSup" << level <<
" = " << (MATLAB)eiSup << endl;
134 sotDEBUG(5) <<
"Je" << level <<
" = " << (MATLAB)Je << endl;
135 sotDEBUG(5) <<
"Ji" << level <<
" = " << (MATLAB)Ji << endl;
136 gettimeofday(&t0, NULL);
141 eiInfs.push_back(eiInf);
142 eiSups.push_back(eiSup);
143 bounds.push_back(eiBoundSide);
145 sotDEBUG(1) <<
"--- Level " << level << std::endl;
146 SolverHierarchicalInequalities *solver =
147 new SolverHierarchicalInequalities(nJ, Qh, Rh, constraintH);
148 solver->initConstraintSize(Je.size1() + Ji.size1());
149 if (solvers.size() == 0)
150 solver->setInitialConditionVoid();
152 SolverHierarchicalInequalities *solverPrec = solvers.back();
153 solver->setInitialCondition(solverPrec->u0, solverPrec->rankh);
155 solver->recordInitialConditions();
156 solvers.push_back(solver);
158 solver->solve(Je, ee, Ji, eiInf, eiSup, eiBoundSide);
160 gettimeofday(&t1, NULL);
161 dtsolver = (t1.tv_sec - t0.tv_sec) * 1000. +
162 (t1.tv_usec - t0.tv_usec + 0.) / 1000.;
163 sotDEBUG(1) <<
"u" << level <<
" = " << (MATLAB)solver->u0 << endl;
164 cout <<
"dtSOLV_" << level <<
" = " << dtsolver <<
"%ms" << endl;
166 solver->computeDifferentialCondition();
172 for (
unsigned int repet = 0; repet < 1; ++repet) {
173 cout <<
"Repet = " << repet << std::endl;
174 struct timeval t0, t1;
176 constraintH.resize(0);
179 for (
unsigned int level = 0; level < solvers.size(); ++level) {
180 gettimeofday(&t0, NULL);
181 sotDEBUG(1) <<
"--- Level " << level << std::endl;
182 SolverHierarchicalInequalities *solver = solvers[level];
184 solver->setInitialConditionVoid();
186 SolverHierarchicalInequalities *solverPrec = solvers[level - 1];
187 solver->setInitialCondition(solverPrec->u0, solverPrec->rankh);
190 solver->recordInitialConditions();
191 #ifdef WITH_WARM_START
194 gettimeofday(&t1, NULL);
195 dtsolver = (t1.tv_sec - t0.tv_sec) * 1000. +
196 (t1.tv_usec - t0.tv_usec + 0.) / 1000.;
197 sotDEBUG(1) <<
"dtWS_" << level <<
" = " << dtsolver <<
"%ms"
199 cout <<
" " << dtsolver;
200 gettimeofday(&t0, NULL);
202 solver->solve(Jes[level], ees[level], Jis[level], eiInfs[level],
203 eiSups[level],
bounds[level], solver->getSlackActiveSet());
204 sotDEBUG(1) <<
"u" << level <<
" = " << (MATLAB)solver->u0 << endl;
205 gettimeofday(&t1, NULL);
206 dtsolver = (t1.tv_sec - t0.tv_sec) * 1000. +
207 (t1.tv_usec - t0.tv_usec + 0.) / 1000.;
208 cout <<
" " << dtsolver;
209 sotDEBUG(1) <<
"dtSOLV_" << level <<
" = " << dtsolver <<
"%ms"
211 gettimeofday(&t0, NULL);
213 solver->computeDifferentialCondition();
218 gettimeofday(&t1, NULL);
219 dtsolver = (t1.tv_sec - t0.tv_sec) * 1000. +
220 (t1.tv_usec - t0.tv_usec + 0.) / 1000.;
221 cout <<
" " << dtsolver;
222 sotDEBUG(1) <<
"dtREC_" << level <<
" = " << dtsolver <<
"%ms"
225 std::cout << std::endl;
228 for (
unsigned int level = 0; level < solvers.size(); ++level) {
229 delete solvers[level];
230 solvers[level] = NULL;
236 void deparse(std::vector<bubMatrix> Jes, std::vector<bubVector> ees,
237 std::vector<bubMatrix> Jis, std::vector<bubVector> eiInfs,
238 std::vector<bubVector> eiSups,
239 std::vector<ConstraintMem::BoundSideVector> bounds) {
241 cout <<
"variable size " << Jes[0].size2() << endl;
243 for (
unsigned int i = 0;
i < Jes.size(); ++
i) {
244 bubMatrix &Je = Jes[
i];
245 bubMatrix &Ji = Jis[
i];
246 bubVector &ee = ees[
i];
247 bubVector &eiInf = eiInfs[
i];
248 bubVector &eiSup = eiSups[
i];
249 ConstraintMem::BoundSideVector &boundSide =
bounds[
i];
255 <<
"equalities " << ee.size() << endl;
257 for (
unsigned int i = 0;
i < ee.size(); ++
i) {
258 for (
unsigned int j = 0; j < Je.size2(); ++j) cout << Je(
i, j) <<
" ";
259 cout <<
"\t" << ee(
i) << endl;
262 unsigned int nbIneq = 0;
263 for (
unsigned int i = 0;
i < boundSide.size(); ++
i) {
264 if (boundSide[
i] & ConstraintMem::BOUND_INF) nbIneq++;
265 if (boundSide[
i] & ConstraintMem::BOUND_SUP) nbIneq++;
268 cout << endl <<
"inequalities " << nbIneq << endl;
269 if (eiInf.size() > 0)
270 for (
unsigned int i = 0;
i < eiInf.size(); ++
i) {
271 if (boundSide[
i] & ConstraintMem::BOUND_INF) {
272 for (
unsigned int j = 0; j < Ji.size2(); ++j)
273 cout << -Ji(
i, j) <<
" ";
275 <<
" X " << -eiInf(
i) << endl;
277 if (boundSide[
i] & ConstraintMem::BOUND_SUP) {
278 for (
unsigned int j = 0; j < Ji.size2(); ++j) cout << Ji(
i, j) <<
" ";
280 <<
" X " << eiSup(
i) << endl;
284 sotDEBUG(15) << endl << endl <<
"end" << endl;
289 std::ifstream off(
filename.c_str());
295 bubVector ee, eiInf, eiSup;
296 ConstraintMem::BoundSideVector eiBoundSide;
299 if (bs !=
"variable") {
300 cerr <<
"!! '" << bs <<
"'" << endl;
305 cerr <<
"!! '" << bs <<
"'" << endl;
311 std::vector<bubMatrix> Jes;
312 std::vector<bubVector> ees;
313 std::vector<bubMatrix> Jis;
314 std::vector<bubVector> eiInfs;
315 std::vector<bubVector> eiSups;
316 std::vector<ConstraintMem::BoundSideVector>
bounds;
318 for (
unsigned int level = 0;; ++level) {
322 }
else if (bs !=
"level") {
323 cerr <<
"!! '" << bs <<
"'" << endl;
327 if (bs !=
"equalities") {
328 cerr <<
"!! '" << bs <<
"'" << endl;
335 for (
int i = 0;
i < me; ++
i) {
336 for (
int j = 0; j < nJ; ++j) off >> Je(
i, j);
341 if (bs !=
"inequalities") {
342 cerr <<
"!! '" << bs <<
"'" << endl;
349 eiBoundSide.resize(mi);
351 for (
int i = 0;
i < mi; ++
i) {
352 for (
int j = 0; j < nJ; ++j) off >> Ji(
i, j);
354 eiBoundSide[
i] = ConstraintMem::BOUND_VOID;
357 eiBoundSide[
i] = (ConstraintMem::BoundSideType)(
358 eiBoundSide[
i] | ConstraintMem::BOUND_INF);
359 eiInf(
i) = atof(number.c_str());
365 eiBoundSide[
i] = (ConstraintMem::BoundSideType)(
366 eiBoundSide[
i] | ConstraintMem::BOUND_SUP);
367 eiSup(
i) = atof(number.c_str());
377 eiInfs.push_back(eiInf);
378 eiSups.push_back(eiSup);
379 bounds.push_back(eiBoundSide);
388 void randBound(ConstraintMem::BoundSideVector &M,
const unsigned int row) {
390 for (
unsigned int i = 0;
i < row; ++
i) {
391 double c = ((
rand() + 0.0) / RAND_MAX * 2) - 1.;
393 M[
i] = ConstraintMem::BOUND_INF;
395 M[
i] = ConstraintMem::BOUND_SUP;
399 void randTest(
const unsigned int nJ,
const bool enableSolve[]) {
400 bubVector eiInf0(nJ), ee0(1), eiSup0(nJ);
401 bubMatrix Ji0(nJ, nJ), Je0(1, nJ);
402 ConstraintMem::BoundSideVector bound0(nJ, ConstraintMem::BOUND_INF);
403 Ji0.assign(bub::identity_matrix<double>(nJ));
404 eiInf0.assign(bub::zero_vector<double>(nJ));
405 Je0.assign(bub::zero_matrix<double>(1, nJ));
410 bubVector ee1, eiInf1, eiSup1;
412 ConstraintMem::BoundSideVector bound1;
414 sotDEBUG(15) <<
"ee1 = " << (MATLAB)ee1 << std::endl;
415 randVector(eiInf1, 3);
416 sotDEBUG(15) <<
"eiInf1 = " << (MATLAB)eiInf1 << std::endl;
417 randVector(eiSup1, 3);
418 sotDEBUG(15) <<
"eiSup1 = " << (MATLAB)eiSup1 << std::endl;
422 randMatrix(Ji1, 3, nJ);
423 sotDEBUG(15) <<
"Ji1 = " << (MATLAB)Ji1 << std::endl;
425 randMatrix(xhi, 1, 3);
427 bub::project(Je1, bub::range(2, 3), bub::range(0, nJ)) = bub::prod(xhi, Je1);
428 sotDEBUG(15) <<
"Je1 = " << (MATLAB)Je1 << std::endl;
433 bubVector ee2, eiInf2, eiSup2;
435 ConstraintMem::BoundSideVector bound2;
437 sotDEBUG(15) <<
"ee2 = " << (MATLAB)ee2 << std::endl;
438 randVector(eiInf2, 3);
439 sotDEBUG(15) <<
"eiInf2 = " << (MATLAB)eiInf2 << std::endl;
440 randVector(eiSup2, 3);
441 sotDEBUG(15) <<
"eiSup2 = " << (MATLAB)eiSup2 << std::endl;
443 randMatrix(Je2, 3, nJ);
444 sotDEBUG(15) <<
"Je2 = " << (MATLAB)Je2 << std::endl;
445 randMatrix(Ji2, 3, nJ);
446 sotDEBUG(15) <<
"Ji2 = " << (MATLAB)Ji2 << std::endl;
448 bubVector ee3, eiInf3, eiSup3;
450 ConstraintMem::BoundSideVector bound3;
452 sotDEBUG(15) <<
"ee3 = " << (MATLAB)ee3 << std::endl;
453 randVector(eiInf3, 3);
454 sotDEBUG(15) <<
"eiInf3 = " << (MATLAB)eiInf3 << std::endl;
455 randVector(eiSup3, 3);
456 sotDEBUG(15) <<
"eiSup3 = " << (MATLAB)eiSup3 << std::endl;
458 randMatrix(Je3, 3, nJ);
459 sotDEBUG(15) <<
"Je3 = " << (MATLAB)Je3 << std::endl;
460 randMatrix(Ji3, 3, nJ);
461 sotDEBUG(15) <<
"Ji3 = " << (MATLAB)Ji3 << std::endl;
463 bubVector ee4(nJ), eiInf4(1), eiSup4(1);
464 bubMatrix Je4(nJ, nJ), Ji4(1, nJ);
465 ConstraintMem::BoundSideVector bound4(1, ConstraintMem::BOUND_INF);
466 Je4.assign(bub::identity_matrix<double>(nJ));
467 ee4.assign(bub::zero_vector<double>(nJ));
468 Ji4.assign(bub::zero_matrix<double>(1, nJ));
469 eiInf4.assign(bub::zero_vector<double>(1));
470 eiSup4.assign(bub::zero_vector<double>(1));
473 std::vector<bubMatrix> Jes;
474 std::vector<bubVector> ees;
475 std::vector<bubMatrix> Jis;
476 std::vector<bubVector> eiInfs;
477 std::vector<bubVector> eiSups;
478 std::vector<ConstraintMem::BoundSideVector>
bounds;
480 if (enableSolve[0]) {
484 eiInfs.push_back(eiInf0);
485 eiSups.push_back(eiSup0);
489 if (enableSolve[1]) {
493 eiInfs.push_back(eiInf1);
494 eiSups.push_back(eiSup1);
497 if (enableSolve[2]) {
501 eiInfs.push_back(eiInf2);
502 eiSups.push_back(eiSup2);
505 if (enableSolve[3]) {
509 eiInfs.push_back(eiInf3);
510 eiSups.push_back(eiSup3);
513 if (enableSolve[4]) {
517 eiInfs.push_back(eiInf4);
518 eiSups.push_back(eiSup4);
523 sotRotationComposedInExtenso Qh(nJ);
525 SolverHierarchicalInequalities::ConstraintList constraintH;
526 SolverHierarchicalInequalities solver(nJ, Qh, Rh, constraintH);
527 solver.initConstraintSize((40 + 6 + 6 + 6 + 40) + 2);
531 if (enableSolve[0]) solver.solve(Je0, ee0, Ji0, eiInf0, eiSup0, bound0);
532 sotDEBUG(15) <<
"/* ----------------------------------------------------- */"
534 sotDEBUG(15) <<
"/* ----------------------------------------------------- */"
536 sotDEBUG(15) <<
"/* ----------------------------------------------------- */"
538 if (enableSolve[1]) solver.solve(Je1, ee1, Ji1, eiInf1, eiSup1, bound1);
539 sotDEBUG(15) <<
"/* ----------------------------------------------------- */"
541 sotDEBUG(15) <<
"/* ----------------------------------------------------- */"
543 sotDEBUG(15) <<
"/* ----------------------------------------------------- */"
545 if (enableSolve[2]) solver.solve(Je2, ee2, Ji2, eiInf2, eiSup2, bound2);
546 sotDEBUG(15) <<
"/* ----------------------------------------------------- */"
548 sotDEBUG(15) <<
"/* ----------------------------------------------------- */"
550 sotDEBUG(15) <<
"/* ----------------------------------------------------- */"
552 if (enableSolve[3]) solver.solve(Je3, ee3, Ji3, eiInf3, eiSup3, bound3);
553 sotDEBUG(15) <<
"/* ----------------------------------------------------- */"
555 sotDEBUG(15) <<
"/* ----------------------------------------------------- */"
557 sotDEBUG(15) <<
"/* ----------------------------------------------------- */"
559 if (enableSolve[4]) solver.solve(Je4, ee4, Ji4, eiInf4, eiSup4, bound4);
575 for (
int i = 0;
i < 10; ++
i)
578 parseTest(
"/home/nmansard/src/StackOfTasks/tests//t.txt");