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");