Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / intern / itasc / WDLSSolver.cpp
1 /* $Id: WDLSSolver.cpp 20749 2009-06-09 11:27:30Z ben2610 $
2  * WDLSSolver.hpp.cpp
3  *
4  *  Created on: Jan 8, 2009
5  *      Author: rubensmits
6  */
7
8 #include "WDLSSolver.hpp"
9 #include "kdl/utilities/svd_eigen_HH.hpp"
10
11 namespace iTaSC {
12
13 WDLSSolver::WDLSSolver() : m_lambda(0.5), m_epsilon(0.1) 
14 {
15         // maximum joint velocity
16         m_qmax = 50.0;
17 }
18
19 WDLSSolver::~WDLSSolver() {
20 }
21
22 bool WDLSSolver::init(unsigned int nq, unsigned int nc, const std::vector<bool>& gc)
23 {
24         m_ns = std::min(nc,nq);
25     m_AWq = e_zero_matrix(nc,nq);
26     m_WyAWq = e_zero_matrix(nc,nq);
27     m_U = e_zero_matrix(nc,nq);
28         m_S = e_zero_vector(std::max(nc,nq));
29     m_temp = e_zero_vector(nq);
30     m_V = e_zero_matrix(nq,nq);
31     m_WqV = e_zero_matrix(nq,nq);
32         m_Wy_ydot = e_zero_vector(nc);
33     return true;
34 }
35
36 bool WDLSSolver::solve(const e_matrix& A, const e_vector& Wy, const e_vector& ydot, const e_matrix& Wq, e_vector& qdot, e_scalar& nlcoef)
37 {
38         double alpha, vmax, norm;
39         // Create the Weighted jacobian
40     m_AWq = A*Wq;
41         for (int i=0; i<Wy.size(); i++)
42                 m_WyAWq.row(i) = Wy(i)*m_AWq.row(i);
43
44     // Compute the SVD of the weighted jacobian
45         int ret = KDL::svd_eigen_HH(m_WyAWq,m_U,m_S,m_V,m_temp);
46     if(ret<0)
47         return false;
48
49     m_WqV = (Wq*m_V).lazy();
50
51     //Wy*ydot
52         m_Wy_ydot = Wy.cwise() * ydot;
53     //S^-1*U'*Wy*ydot
54         e_scalar maxDeltaS = e_scalar(0.0);
55         e_scalar prevS = e_scalar(0.0);
56         e_scalar maxS = e_scalar(1.0);
57         e_scalar S, lambda;
58         qdot.setZero();
59         for(int i=0;i<m_ns;++i) {
60                 S = m_S(i);
61                 if (S <= KDL::epsilon)
62                         break;
63                 if (i > 0 && (prevS-S) > maxDeltaS) {
64                         maxDeltaS = (prevS-S);
65                         maxS = prevS;
66                 }
67                 lambda = (S < m_epsilon) ? (e_scalar(1.0)-KDL::sqr(S/m_epsilon))*m_lambda*m_lambda : e_scalar(0.0);
68                 alpha = m_U.col(i).dot(m_Wy_ydot)*S/(S*S+lambda);
69                 vmax = m_WqV.col(i).cwise().abs().maxCoeff();
70                 norm = fabs(alpha*vmax);
71                 if (norm > m_qmax) {
72                         qdot += m_WqV.col(i)*(alpha*m_qmax/norm);
73                 } else {
74                         qdot += m_WqV.col(i)*alpha;
75                 }
76                 prevS = S;
77         }
78         if (maxDeltaS == e_scalar(0.0))
79                 nlcoef = e_scalar(KDL::epsilon);
80         else
81                 nlcoef = (maxS-maxDeltaS)/maxS;
82     return true;
83 }
84
85 }