Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / intern / itasc / WSDLSSolver.cpp
1 /* $Id: WSDLSSolver.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 "WSDLSSolver.hpp"
9 #include "kdl/utilities/svd_eigen_HH.hpp"
10 #include <cstdio>
11
12 namespace iTaSC {
13
14 WSDLSSolver::WSDLSSolver() :
15         m_ns(0), m_nc(0), m_nq(0)
16
17 {
18         // default maximum speed: 50 rad/s
19         m_qmax = 50.0;
20 }
21
22 WSDLSSolver::~WSDLSSolver() {
23 }
24
25 bool WSDLSSolver::init(unsigned int _nq, unsigned int _nc, const std::vector<bool>& gc)
26 {
27         if (_nc == 0 || _nq == 0 || gc.size() != _nc)
28                 return false;
29         m_nc = _nc;
30         m_nq = _nq;
31         m_ns = std::min(m_nc,m_nq);
32     m_AWq = e_zero_matrix(m_nc,m_nq);
33     m_WyAWq = e_zero_matrix(m_nc,m_nq);
34     m_U = e_zero_matrix(m_nc,m_nq);
35         m_S = e_zero_vector(std::max(m_nc,m_nq));
36     m_temp = e_zero_vector(m_nq);
37     m_V = e_zero_matrix(m_nq,m_nq);
38     m_WqV = e_zero_matrix(m_nq,m_nq);
39         m_Wy_ydot = e_zero_vector(m_nc);
40         m_ytask = gc;
41     return true;
42 }
43
44 bool WSDLSSolver::solve(const e_matrix& A, const e_vector& Wy, const e_vector& ydot, const e_matrix& Wq, e_vector& qdot, e_scalar& nlcoef)
45 {
46         unsigned int i, j, l;
47         e_scalar N, M;
48
49         // Create the Weighted jacobian
50     m_AWq = (A*Wq).lazy();
51         for (i=0; i<m_nc; i++)
52                 m_WyAWq.row(i) = Wy(i)*m_AWq.row(i);
53
54     // Compute the SVD of the weighted jacobian
55         int ret = KDL::svd_eigen_HH(m_WyAWq,m_U,m_S,m_V,m_temp);
56     if(ret<0)
57         return false;
58
59         m_Wy_ydot = Wy.cwise() * ydot;
60     m_WqV = (Wq*m_V).lazy();
61         qdot.setZero();
62         e_scalar maxDeltaS = e_scalar(0.0);
63         e_scalar prevS = e_scalar(0.0);
64         e_scalar maxS = e_scalar(1.0);
65         for(i=0;i<m_ns;++i) {
66                 e_scalar norm, mag, alpha, _qmax, Sinv, vmax, damp;
67                 e_scalar S = m_S(i);
68                 bool prev;
69                 if (S < KDL::epsilon)
70                         break;
71                 Sinv = e_scalar(1.)/S;
72                 if (i > 0) {
73                         if ((prevS-S) > maxDeltaS) {
74                                 maxDeltaS = (prevS-S);
75                                 maxS = prevS;
76                         }
77                 }
78                 N = M = e_scalar(0.);
79                 for (l=0, prev=m_ytask[0], norm=e_scalar(0.); l<m_nc; l++) {
80                         if (prev == m_ytask[l]) {
81                                 norm += m_U(l,i)*m_U(l,i);
82                         } else {
83                                 N += std::sqrt(norm);
84                                 norm = m_U(l,i)*m_U(l,i);
85                         }
86                         prev = m_ytask[l];
87                 }
88                 N += std::sqrt(norm);
89                 for (j=0; j<m_nq; j++) {
90                         for (l=0, prev=m_ytask[0], norm=e_scalar(0.), mag=e_scalar(0.); l<m_nc; l++) {
91                                 if (prev == m_ytask[l]) {
92                                         norm += m_WyAWq(l,j)*m_WyAWq(l,j);
93                                 } else {
94                                         mag += std::sqrt(norm);
95                                         norm = m_WyAWq(l,j)*m_WyAWq(l,j);
96                                 }
97                                 prev = m_ytask[l];
98                         }
99                         mag += std::sqrt(norm);
100                         M += fabs(m_V(j,i))*mag;
101                 }
102                 M *= Sinv;
103                 alpha = m_U.col(i).dot(m_Wy_ydot);
104                 _qmax = (N < M) ? m_qmax*N/M : m_qmax;
105                 vmax = m_WqV.col(i).cwise().abs().maxCoeff();
106                 norm = fabs(Sinv*alpha*vmax);
107                 if (norm > _qmax) {
108                         damp = Sinv*alpha*_qmax/norm;
109                 } else {
110                         damp = Sinv*alpha;
111                 }
112                 qdot += m_WqV.col(i)*damp;
113                 prevS = S;
114         }
115         if (maxDeltaS == e_scalar(0.0))
116                 nlcoef = e_scalar(KDL::epsilon);
117         else
118                 nlcoef = (maxS-maxDeltaS)/maxS;
119     return true;
120 }
121
122 }