remove $Id: tags after discussion on the mailign list: http://markmail.org/message...
[blender.git] / intern / itasc / WDLSSolver.cpp
1 /** \file itasc/WDLSSolver.cpp
2  *  \ingroup itasc
3  */
4 /*
5  * WDLSSolver.hpp.cpp
6  *
7  *  Created on: Jan 8, 2009
8  *      Author: rubensmits
9  */
10
11 #include "WDLSSolver.hpp"
12 #include "kdl/utilities/svd_eigen_HH.hpp"
13
14 namespace iTaSC {
15
16 WDLSSolver::WDLSSolver() : m_lambda(0.5), m_epsilon(0.1) 
17 {
18         // maximum joint velocity
19         m_qmax = 50.0;
20 }
21
22 WDLSSolver::~WDLSSolver() {
23 }
24
25 bool WDLSSolver::init(unsigned int nq, unsigned int nc, const std::vector<bool>& gc)
26 {
27         m_ns = std::min(nc,nq);
28     m_AWq = e_zero_matrix(nc,nq);
29     m_WyAWq = e_zero_matrix(nc,nq);
30     m_WyAWqt = e_zero_matrix(nq,nc);
31         m_S = e_zero_vector(std::max(nc,nq));
32         m_Wy_ydot = e_zero_vector(nc);
33         if (nq > nc) {
34                 m_transpose = true;
35             m_temp = e_zero_vector(nc);
36             m_U = e_zero_matrix(nc,nc);
37                 m_V = e_zero_matrix(nq,nc);
38             m_WqV = e_zero_matrix(nq,nc);
39         } else {
40                 m_transpose = false;
41             m_temp = e_zero_vector(nq);
42             m_U = e_zero_matrix(nc,nq);
43                 m_V = e_zero_matrix(nq,nq);
44             m_WqV = e_zero_matrix(nq,nq);
45         }
46     return true;
47 }
48
49 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)
50 {
51         double alpha, vmax, norm;
52         // Create the Weighted jacobian
53     m_AWq = A*Wq;
54         for (int i=0; i<Wy.size(); i++)
55                 m_WyAWq.row(i) = Wy(i)*m_AWq.row(i);
56
57     // Compute the SVD of the weighted jacobian
58         int ret;
59         if (m_transpose) {
60                 m_WyAWqt = m_WyAWq.transpose();
61                 ret = KDL::svd_eigen_HH(m_WyAWqt,m_V,m_S,m_U,m_temp);
62         } else {
63                 ret = KDL::svd_eigen_HH(m_WyAWq,m_U,m_S,m_V,m_temp);
64         }
65     if(ret<0)
66         return false;
67
68     m_WqV.noalias() = Wq*m_V;
69
70     //Wy*ydot
71         m_Wy_ydot = Wy.array() * ydot.array();
72     //S^-1*U'*Wy*ydot
73         e_scalar maxDeltaS = e_scalar(0.0);
74         e_scalar prevS = e_scalar(0.0);
75         e_scalar maxS = e_scalar(1.0);
76         e_scalar S, lambda;
77         qdot.setZero();
78         for(int i=0;i<m_ns;++i) {
79                 S = m_S(i);
80                 if (S <= KDL::epsilon)
81                         break;
82                 if (i > 0 && (prevS-S) > maxDeltaS) {
83                         maxDeltaS = (prevS-S);
84                         maxS = prevS;
85                 }
86                 lambda = (S < m_epsilon) ? (e_scalar(1.0)-KDL::sqr(S/m_epsilon))*m_lambda*m_lambda : e_scalar(0.0);
87                 alpha = m_U.col(i).dot(m_Wy_ydot)*S/(S*S+lambda);
88                 vmax = m_WqV.col(i).array().abs().maxCoeff();
89                 norm = fabs(alpha*vmax);
90                 if (norm > m_qmax) {
91                         qdot += m_WqV.col(i)*(alpha*m_qmax/norm);
92                 } else {
93                         qdot += m_WqV.col(i)*alpha;
94                 }
95                 prevS = S;
96         }
97         if (maxDeltaS == e_scalar(0.0))
98                 nlcoef = e_scalar(KDL::epsilon);
99         else
100                 nlcoef = (maxS-maxDeltaS)/maxS;
101     return true;
102 }
103
104 }