|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
#include "toonz/ikjacobian.h"
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
#include <stdlib.h></stdlib.h>
|
|
Toshihiro Shimizu |
890ddd |
#include <math.h></math.h>
|
|
Toshihiro Shimizu |
890ddd |
#include <assert.h></assert.h>
|
|
Toshihiro Shimizu |
890ddd |
#include <iostream></iostream>
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
#include "tstopwatch.h"
|
|
Toshihiro Shimizu |
890ddd |
using namespace std;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
inline bool NearZero(double x, double tolerance)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
return (fabs(x) <= tolerance);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
/*
|
|
Toshihiro Shimizu |
890ddd |
#ifdef _DYNAMIC
|
|
Toshihiro Shimizu |
890ddd |
const double BASEMAXDIST = 0.02;
|
|
Toshihiro Shimizu |
890ddd |
#else
|
|
Toshihiro Shimizu |
890ddd |
const double MAXDIST = 0.08;
|
|
Toshihiro Shimizu |
890ddd |
#endif
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
const double DELTA = 0.4;
|
|
Toshihiro Shimizu |
890ddd |
const long double LAMBDA = 2.0; // solo per DLS. ottimale : 0.24
|
|
Toshihiro Shimizu |
890ddd |
const double NEARZERO = 0.0000000001;
|
|
Toshihiro Shimizu |
890ddd |
*/
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
//*******************************************************
|
|
Toshihiro Shimizu |
890ddd |
// Class VectorRn
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
VectorRn VectorRn::WorkVector;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
double VectorRn::MaxAbs() const
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
double result = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
double *t = x;
|
|
Toshihiro Shimizu |
890ddd |
for (long i = length; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
if ((*t) > result) {
|
|
Toshihiro Shimizu |
890ddd |
result = *t;
|
|
Toshihiro Shimizu |
890ddd |
} else if (-(*t) > result) {
|
|
Toshihiro Shimizu |
890ddd |
result = -(*t);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
t++;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
return result;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
//*************************************************************************
|
|
Toshihiro Shimizu |
890ddd |
// MatrixRmn
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn MatrixRmn::WorkMatrix; // Temporary work matrix
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Fill the diagonal entries with the value d. The rest of the matrix is unchanged.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::SetDiagonalEntries(double d)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Shinya Kitaoka |
12c444 |
long diagLen = std::min(NumRows, NumCols);
|
|
Toshihiro Shimizu |
890ddd |
double *dPtr = x;
|
|
Toshihiro Shimizu |
890ddd |
for (; diagLen > 0; diagLen--) {
|
|
Toshihiro Shimizu |
890ddd |
*dPtr = d;
|
|
Toshihiro Shimizu |
890ddd |
dPtr += NumRows + 1;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Fill the diagonal entries with values in vector d. The rest of the matrix is unchanged.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::SetDiagonalEntries(const VectorRn &d)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Shinya Kitaoka |
12c444 |
long diagLen = std::min(NumRows, NumCols);
|
|
Toshihiro Shimizu |
890ddd |
assert(d.length == diagLen);
|
|
Toshihiro Shimizu |
890ddd |
double *dPtr = x;
|
|
Toshihiro Shimizu |
890ddd |
double *from = d.x;
|
|
Toshihiro Shimizu |
890ddd |
for (; diagLen > 0; diagLen--) {
|
|
Toshihiro Shimizu |
890ddd |
*dPtr = *(from++);
|
|
Toshihiro Shimizu |
890ddd |
dPtr += NumRows + 1;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Fill the superdiagonal entries with the value d. The rest of the matrix is unchanged.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::SetSuperDiagonalEntries(double d)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Shinya Kitaoka |
12c444 |
long sDiagLen = std::min(NumRows, (long)(NumCols - 1));
|
|
Toshihiro Shimizu |
890ddd |
double *to = x + NumRows;
|
|
Toshihiro Shimizu |
890ddd |
for (; sDiagLen > 0; sDiagLen--) {
|
|
Toshihiro Shimizu |
890ddd |
*to = d;
|
|
Toshihiro Shimizu |
890ddd |
to += NumRows + 1;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Fill the superdiagonal entries with values in vector d. The rest of the matrix is unchanged.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::SetSuperDiagonalEntries(const VectorRn &d)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Shinya Kitaoka |
12c444 |
long sDiagLen = std::min((long)(NumRows - 1), NumCols);
|
|
Toshihiro Shimizu |
890ddd |
assert(sDiagLen == d.length);
|
|
Toshihiro Shimizu |
890ddd |
double *to = x + NumRows;
|
|
Toshihiro Shimizu |
890ddd |
double *from = d.x;
|
|
Toshihiro Shimizu |
890ddd |
for (; sDiagLen > 0; sDiagLen--) {
|
|
Toshihiro Shimizu |
890ddd |
*to = *(from++);
|
|
Toshihiro Shimizu |
890ddd |
to += NumRows + 1;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Fill the subdiagonal entries with the value d. The rest of the matrix is unchanged.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::SetSubDiagonalEntries(double d)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Shinya Kitaoka |
12c444 |
long sDiagLen = std::min(NumRows, NumCols) - 1;
|
|
Toshihiro Shimizu |
890ddd |
double *to = x + 1;
|
|
Toshihiro Shimizu |
890ddd |
for (; sDiagLen > 0; sDiagLen--) {
|
|
Toshihiro Shimizu |
890ddd |
*to = d;
|
|
Toshihiro Shimizu |
890ddd |
to += NumRows + 1;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Fill the subdiagonal entries with values in vector d. The rest of the matrix is unchanged.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::SetSubDiagonalEntries(const VectorRn &d)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Shinya Kitaoka |
12c444 |
long sDiagLen = std::min(NumRows, NumCols) - 1;
|
|
Toshihiro Shimizu |
890ddd |
assert(sDiagLen == d.length);
|
|
Toshihiro Shimizu |
890ddd |
double *to = x + 1;
|
|
Toshihiro Shimizu |
890ddd |
double *from = d.x;
|
|
Toshihiro Shimizu |
890ddd |
for (; sDiagLen > 0; sDiagLen--) {
|
|
Toshihiro Shimizu |
890ddd |
*to = *(from++);
|
|
Toshihiro Shimizu |
890ddd |
to += NumRows + 1;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Set the i-th column equal to d.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::SetColumn(long i, const VectorRn &d)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(NumRows == d.GetLength());
|
|
Toshihiro Shimizu |
890ddd |
double *to = x + i * NumRows;
|
|
Toshihiro Shimizu |
890ddd |
const double *from = d.x;
|
|
Toshihiro Shimizu |
890ddd |
for (i = NumRows; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
*(to++) = *(from++);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Set the i-th column equal to d.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::SetRow(long i, const VectorRn &d)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(NumCols == d.GetLength());
|
|
Toshihiro Shimizu |
890ddd |
double *to = x + i;
|
|
Toshihiro Shimizu |
890ddd |
const double *from = d.x;
|
|
Toshihiro Shimizu |
890ddd |
for (i = NumRows; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
*to = *(from++);
|
|
Toshihiro Shimizu |
890ddd |
to += NumRows;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Sets a "linear" portion of the array with the values from a vector d
|
|
Toshihiro Shimizu |
890ddd |
// The first row and column position are given by startRow, startCol.
|
|
Toshihiro Shimizu |
890ddd |
// Successive positions are found by using the deltaRow, deltaCol values
|
|
Toshihiro Shimizu |
890ddd |
// to increment the row and column indices. There is no wrapping around.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::SetSequence(const VectorRn &d, long startRow, long startCol, long deltaRow, long deltaCol)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
long length = d.length;
|
|
Toshihiro Shimizu |
890ddd |
assert(startRow >= 0 && startRow < NumRows && startCol >= 0 && startCol < NumCols);
|
|
Toshihiro Shimizu |
890ddd |
assert(startRow + (length - 1) * deltaRow >= 0 && startRow + (length - 1) * deltaRow < NumRows);
|
|
Toshihiro Shimizu |
890ddd |
assert(startCol + (length - 1) * deltaCol >= 0 && startCol + (length - 1) * deltaCol < NumCols);
|
|
Toshihiro Shimizu |
890ddd |
double *to = x + startRow + NumRows * startCol;
|
|
Toshihiro Shimizu |
890ddd |
double *from = d.x;
|
|
Toshihiro Shimizu |
890ddd |
long stride = deltaRow + NumRows * deltaCol;
|
|
Toshihiro Shimizu |
890ddd |
for (; length > 0; length--) {
|
|
Toshihiro Shimizu |
890ddd |
*to = *(from++);
|
|
Toshihiro Shimizu |
890ddd |
to += stride;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// The matrix A is loaded, in into "this" matrix, based at (0,0).
|
|
Toshihiro Shimizu |
890ddd |
// The size of "this" matrix must be large enough to accomodate A.
|
|
Toshihiro Shimizu |
890ddd |
// The rest of "this" matrix is left unchanged. It is not filled with zeroes!
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::LoadAsSubmatrix(const MatrixRmn &A)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(A.NumRows <= NumRows && A.NumCols <= NumCols);
|
|
Toshihiro Shimizu |
890ddd |
int extraColStep = NumRows - A.NumRows;
|
|
Toshihiro Shimizu |
890ddd |
double *to = x;
|
|
Toshihiro Shimizu |
890ddd |
double *from = A.x;
|
|
Toshihiro Shimizu |
890ddd |
for (long i = A.NumCols; i > 0; i--) { // Copy columns of A, one per time thru loop
|
|
Toshihiro Shimizu |
890ddd |
for (long j = A.NumRows; j > 0; j--) { // Copy all elements of this column of A
|
|
Toshihiro Shimizu |
890ddd |
*(to++) = *(from++);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
to += extraColStep;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// The matrix A is loaded, in transposed order into "this" matrix, based at (0,0).
|
|
Toshihiro Shimizu |
890ddd |
// The size of "this" matrix must be large enough to accomodate A.
|
|
Toshihiro Shimizu |
890ddd |
// The rest of "this" matrix is left unchanged. It is not filled with zeroes!
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::LoadAsSubmatrixTranspose(const MatrixRmn &A)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(A.NumRows <= NumCols && A.NumCols <= NumRows);
|
|
Toshihiro Shimizu |
890ddd |
double *rowPtr = x;
|
|
Toshihiro Shimizu |
890ddd |
double *from = A.x;
|
|
Toshihiro Shimizu |
890ddd |
for (long i = A.NumCols; i > 0; i--) { // Copy columns of A, once per loop
|
|
Toshihiro Shimizu |
890ddd |
double *to = rowPtr;
|
|
Toshihiro Shimizu |
890ddd |
for (long j = A.NumRows; j > 0; j--) { // Loop copying values from the column of A
|
|
Toshihiro Shimizu |
890ddd |
*to = *(from++);
|
|
Toshihiro Shimizu |
890ddd |
to += NumRows;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
rowPtr++;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Calculate the Frobenius Norm (square root of sum of squares of entries of the matrix)
|
|
Toshihiro Shimizu |
890ddd |
double MatrixRmn::FrobeniusNorm() const
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
return sqrt(FrobeniusNormSq());
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Multiply this matrix by column vector v.
|
|
Toshihiro Shimizu |
890ddd |
// Result is column vector "result"
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::Multiply(const VectorRn &v, VectorRn &result) const
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(v.GetLength() == NumCols && result.GetLength() == NumRows);
|
|
Toshihiro Shimizu |
890ddd |
double *out = result.GetPtr(); // Points to entry in result vector
|
|
Toshihiro Shimizu |
890ddd |
const double *rowPtr = x; // Points to beginning of next row in matrix
|
|
Toshihiro Shimizu |
890ddd |
for (long j = NumRows; j > 0; j--) {
|
|
Toshihiro Shimizu |
890ddd |
const double *in = v.GetPtr();
|
|
Toshihiro Shimizu |
890ddd |
const double *m = rowPtr++;
|
|
Toshihiro Shimizu |
890ddd |
*out = 0.0f;
|
|
Toshihiro Shimizu |
890ddd |
for (long i = NumCols; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
*out += (*(in++)) * (*m);
|
|
Toshihiro Shimizu |
890ddd |
m += NumRows;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
out++;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Multiply transpose of this matrix by column vector v.
|
|
Toshihiro Shimizu |
890ddd |
// Result is column vector "result"
|
|
Toshihiro Shimizu |
890ddd |
// Equivalent to mult by row vector on left
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::MultiplyTranspose(const VectorRn &v, VectorRn &result) const
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(v.GetLength() == NumRows && result.GetLength() == NumCols);
|
|
Toshihiro Shimizu |
890ddd |
double *out = result.GetPtr(); // Points to entry in result vector
|
|
Toshihiro Shimizu |
890ddd |
const double *colPtr = x; // Points to beginning of next column in matrix
|
|
Toshihiro Shimizu |
890ddd |
for (long i = NumCols; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
const double *in = v.GetPtr();
|
|
Toshihiro Shimizu |
890ddd |
*out = 0.0f;
|
|
Toshihiro Shimizu |
890ddd |
for (long j = NumRows; j > 0; j--) {
|
|
Toshihiro Shimizu |
890ddd |
*out += (*(in++)) * (*(colPtr++));
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
out++;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Form the dot product of a vector v with the i-th column of the array
|
|
Toshihiro Shimizu |
890ddd |
double MatrixRmn::DotProductColumn(const VectorRn &v, long colNum) const
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(v.GetLength() == NumRows);
|
|
Toshihiro Shimizu |
890ddd |
double *ptrC = x + colNum * NumRows;
|
|
Toshihiro Shimizu |
890ddd |
double *ptrV = v.x;
|
|
Toshihiro Shimizu |
890ddd |
double ret = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
for (long i = NumRows; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
ret += (*(ptrC++)) * (*(ptrV++));
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
return ret;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Add a constant to each entry on the diagonal
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn &MatrixRmn::AddToDiagonal(double d) // Adds d to each diagonal entry
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Shinya Kitaoka |
12c444 |
long diagLen = std::min(NumRows, NumCols);
|
|
Toshihiro Shimizu |
890ddd |
double *dPtr = x;
|
|
Toshihiro Shimizu |
890ddd |
for (; diagLen > 0; diagLen--) {
|
|
Toshihiro Shimizu |
890ddd |
*dPtr += d;
|
|
Toshihiro Shimizu |
890ddd |
dPtr += NumRows + 1;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
return *this;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Aggiunge i temini del vettore alla diagonale
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn &MatrixRmn::AddToDiagonal(const VectorRn &v) // Adds d to each diagonal entry
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Shinya Kitaoka |
12c444 |
long diagLen = std::min(NumRows, NumCols);
|
|
Toshihiro Shimizu |
890ddd |
double *dPtr = x;
|
|
Toshihiro Shimizu |
890ddd |
const double *dv = v.x;
|
|
Toshihiro Shimizu |
890ddd |
for (; diagLen > 0; diagLen--) {
|
|
Toshihiro Shimizu |
890ddd |
*dPtr += *(dv++);
|
|
Toshihiro Shimizu |
890ddd |
dPtr += NumRows + 1;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
return *this;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn &MatrixRmn::MultiplyScalar(const MatrixRmn &A, double k, MatrixRmn &dst)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
long length = A.NumCols;
|
|
Toshihiro Shimizu |
890ddd |
double *dPtr = dst.x;
|
|
Toshihiro Shimizu |
890ddd |
for (long i = dst.NumCols; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
double *aPtr = A.x; // Points to beginning of row in A
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
for (long j = dst.NumRows; j > 0; j--) {
|
|
Toshihiro Shimizu |
890ddd |
*dPtr = *aPtr * k;
|
|
Toshihiro Shimizu |
890ddd |
dPtr++;
|
|
Toshihiro Shimizu |
890ddd |
aPtr++;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
aPtr += A.NumRows;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
return dst;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
// Multiply two MatrixRmn's
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn &MatrixRmn::Multiply(const MatrixRmn &A, const MatrixRmn &B, MatrixRmn &dst)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(A.NumCols == B.NumRows && A.NumRows == dst.NumRows && B.NumCols == dst.NumCols);
|
|
Toshihiro Shimizu |
890ddd |
long length = A.NumCols;
|
|
Toshihiro Shimizu |
890ddd |
double *bPtr = B.x; // Points to beginning of column in B
|
|
Toshihiro Shimizu |
890ddd |
double *dPtr = dst.x;
|
|
Toshihiro Shimizu |
890ddd |
for (long i = dst.NumCols; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
double *aPtr = A.x; // Points to beginning of row in A
|
|
Toshihiro Shimizu |
890ddd |
for (long j = dst.NumRows; j > 0; j--) {
|
|
Toshihiro Shimizu |
890ddd |
*dPtr = DotArray(length, aPtr, A.NumRows, bPtr, 1);
|
|
Toshihiro Shimizu |
890ddd |
dPtr++;
|
|
Toshihiro Shimizu |
890ddd |
aPtr++;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
bPtr += B.NumRows;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
return dst;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Multiply two MatrixRmn's, Transpose the first matrix before multiplying
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn &MatrixRmn::TransposeMultiply(const MatrixRmn &A, const MatrixRmn &B, MatrixRmn &dst)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(A.NumRows == B.NumRows && A.NumCols == dst.NumRows && B.NumCols == dst.NumCols);
|
|
Toshihiro Shimizu |
890ddd |
long length = A.NumRows;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
double *bPtr = B.x; // bPtr Points to beginning of column in B
|
|
Toshihiro Shimizu |
890ddd |
double *dPtr = dst.x;
|
|
Toshihiro Shimizu |
890ddd |
for (long i = dst.NumCols; i > 0; i--) { // Loop over all columns of dst
|
|
Toshihiro Shimizu |
890ddd |
double *aPtr = A.x; // aPtr Points to beginning of column in A
|
|
Toshihiro Shimizu |
890ddd |
for (long j = dst.NumRows; j > 0; j--) { // Loop over all rows of dst
|
|
Toshihiro Shimizu |
890ddd |
*dPtr = DotArray(length, aPtr, 1, bPtr, 1);
|
|
Toshihiro Shimizu |
890ddd |
dPtr++;
|
|
Toshihiro Shimizu |
890ddd |
aPtr += A.NumRows;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
bPtr += B.NumRows;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
return dst;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Multiply two MatrixRmn's. Transpose the second matrix before multiplying
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn &MatrixRmn::MultiplyTranspose(const MatrixRmn &A, const MatrixRmn &B, MatrixRmn &dst)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(A.NumCols == B.NumCols && A.NumRows == dst.NumRows && B.NumRows == dst.NumCols);
|
|
Toshihiro Shimizu |
890ddd |
long length = A.NumCols;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
double *bPtr = B.x; // Points to beginning of row in B
|
|
Toshihiro Shimizu |
890ddd |
double *dPtr = dst.x;
|
|
Toshihiro Shimizu |
890ddd |
for (long i = dst.NumCols; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
double *aPtr = A.x; // Points to beginning of row in A
|
|
Toshihiro Shimizu |
890ddd |
for (long j = dst.NumRows; j > 0; j--) {
|
|
Toshihiro Shimizu |
890ddd |
*dPtr = DotArray(length, aPtr, A.NumRows, bPtr, B.NumRows);
|
|
Toshihiro Shimizu |
890ddd |
dPtr++;
|
|
Toshihiro Shimizu |
890ddd |
aPtr++;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
bPtr++;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
return dst;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Solves the equation (*this)*xVec = b;
|
|
Toshihiro Shimizu |
890ddd |
// Uses row operations. Assumes *this is square and invertible.
|
|
Toshihiro Shimizu |
890ddd |
// No error checking for divide by zero or instability (except with asserts)
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::Solve(const VectorRn &b, VectorRn *xVec) const
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(NumRows == NumCols && NumCols == xVec->GetLength() && NumRows == b.GetLength());
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Copy this matrix and b into an Augmented Matrix
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn &AugMat = GetWorkMatrix(NumRows, NumCols + 1);
|
|
Toshihiro Shimizu |
890ddd |
AugMat.LoadAsSubmatrix(*this);
|
|
Toshihiro Shimizu |
890ddd |
AugMat.SetColumn(NumRows, b);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Put into row echelon form with row operations
|
|
Toshihiro Shimizu |
890ddd |
AugMat.ConvertToRefNoFree();
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Solve for x vector values using back substitution
|
|
Toshihiro Shimizu |
890ddd |
double *xLast = xVec->x + NumRows - 1; // Last entry in xVec
|
|
Toshihiro Shimizu |
890ddd |
double *endRow = AugMat.x + NumRows * NumCols - 1; // Last entry in the current row of the coefficient part of Augmented Matrix
|
|
Toshihiro Shimizu |
890ddd |
double *bPtr = endRow + NumRows; // Last entry in augmented matrix (end of last column, in augmented part)
|
|
Toshihiro Shimizu |
890ddd |
for (long i = NumRows; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
double accum = *(bPtr--);
|
|
Toshihiro Shimizu |
890ddd |
// Next loop computes back substitution terms
|
|
Toshihiro Shimizu |
890ddd |
double *rowPtr = endRow; // Points to entries of the current row for back substitution.
|
|
Toshihiro Shimizu |
890ddd |
double *xPtr = xLast; // Points to entries in the x vector (also for back substitution)
|
|
Toshihiro Shimizu |
890ddd |
for (long j = NumRows - i; j > 0; j--) {
|
|
Toshihiro Shimizu |
890ddd |
accum -= (*rowPtr) * (*(xPtr--));
|
|
Toshihiro Shimizu |
890ddd |
rowPtr -= NumCols; // Previous entry in the row
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
assert(*rowPtr != 0.0); // Are not supposed to be any free variables in this matrix
|
|
Toshihiro Shimizu |
890ddd |
*xPtr = accum / (*rowPtr);
|
|
Toshihiro Shimizu |
890ddd |
endRow--;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// ConvertToRefNoFree
|
|
Toshihiro Shimizu |
890ddd |
// Converts the matrix (in place) to row echelon form
|
|
Toshihiro Shimizu |
890ddd |
// For us, row echelon form allows any non-zero values, not just 1's, in the
|
|
Toshihiro Shimizu |
890ddd |
// position for a lead variable.
|
|
Toshihiro Shimizu |
890ddd |
// The "NoFree" version operates on the assumption that no free variable will be found.
|
|
Toshihiro Shimizu |
890ddd |
// Algorithm uses row operations and row pivoting (only).
|
|
Toshihiro Shimizu |
890ddd |
// Augmented matrix is correctly accomodated. Only the first square part participates
|
|
Toshihiro Shimizu |
890ddd |
// in the main work of row operations.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::ConvertToRefNoFree()
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
// Loop over all columns (variables)
|
|
Toshihiro Shimizu |
890ddd |
// Find row with most non-zero entry.
|
|
Toshihiro Shimizu |
890ddd |
// Swap to the highest active row
|
|
Toshihiro Shimizu |
890ddd |
// Subtract appropriately from all the lower rows (row op of type 3)
|
|
Shinya Kitaoka |
12c444 |
long numIters = std::min(NumRows, NumCols);
|
|
Toshihiro Shimizu |
890ddd |
double *rowPtr1 = x;
|
|
Toshihiro Shimizu |
890ddd |
const long diagStep = NumRows + 1;
|
|
Toshihiro Shimizu |
890ddd |
long lenRowLeft = NumCols;
|
|
Toshihiro Shimizu |
890ddd |
for (; numIters > 1; numIters--) {
|
|
Toshihiro Shimizu |
890ddd |
// Find row with most non-zero entry.
|
|
Toshihiro Shimizu |
890ddd |
double *rowPtr2 = rowPtr1;
|
|
Toshihiro Shimizu |
890ddd |
double maxAbs = fabs(*rowPtr1);
|
|
Toshihiro Shimizu |
890ddd |
double *rowPivot = rowPtr1;
|
|
Toshihiro Shimizu |
890ddd |
long i;
|
|
Toshihiro Shimizu |
890ddd |
for (i = numIters - 1; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
const double &newMax = *(++rowPivot);
|
|
Toshihiro Shimizu |
890ddd |
if (newMax > maxAbs) {
|
|
Toshihiro Shimizu |
890ddd |
maxAbs = *rowPivot;
|
|
Toshihiro Shimizu |
890ddd |
rowPtr2 = rowPivot;
|
|
Toshihiro Shimizu |
890ddd |
} else if (-newMax > maxAbs) {
|
|
Toshihiro Shimizu |
890ddd |
maxAbs = -newMax;
|
|
Toshihiro Shimizu |
890ddd |
rowPtr2 = rowPivot;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
// Pivot step: Swap the row with highest entry to the current row
|
|
Toshihiro Shimizu |
890ddd |
if (rowPtr1 != rowPtr2) {
|
|
Toshihiro Shimizu |
890ddd |
double *to = rowPtr1;
|
|
Toshihiro Shimizu |
890ddd |
for (long i = lenRowLeft; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
double temp = *to;
|
|
Toshihiro Shimizu |
890ddd |
*to = *rowPtr2;
|
|
Toshihiro Shimizu |
890ddd |
*rowPtr2 = temp;
|
|
Toshihiro Shimizu |
890ddd |
to += NumRows;
|
|
Toshihiro Shimizu |
890ddd |
rowPtr2 += NumRows;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
// Subtract this row appropriately from all the lower rows (row operation of type 3)
|
|
Toshihiro Shimizu |
890ddd |
rowPtr2 = rowPtr1;
|
|
Toshihiro Shimizu |
890ddd |
for (i = numIters - 1; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
rowPtr2++;
|
|
Toshihiro Shimizu |
890ddd |
double *to = rowPtr2;
|
|
Toshihiro Shimizu |
890ddd |
double *from = rowPtr1;
|
|
Toshihiro Shimizu |
890ddd |
assert(*from != 0.0);
|
|
Toshihiro Shimizu |
890ddd |
double alpha = (*to) / (*from);
|
|
Toshihiro Shimizu |
890ddd |
*to = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
for (long j = lenRowLeft - 1; j > 0; j--) {
|
|
Toshihiro Shimizu |
890ddd |
to += NumRows;
|
|
Toshihiro Shimizu |
890ddd |
from += NumRows;
|
|
Toshihiro Shimizu |
890ddd |
*to -= (*from) * alpha;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
// Update for next iteration of loop
|
|
Toshihiro Shimizu |
890ddd |
rowPtr1 += diagStep;
|
|
Toshihiro Shimizu |
890ddd |
lenRowLeft--;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Calculate the c=cosine and s=sine values for a Givens transformation.
|
|
Toshihiro Shimizu |
890ddd |
// The matrix M = ( (c, -s), (s, c) ) in row order transforms the
|
|
Toshihiro Shimizu |
890ddd |
// column vector (a, b)^T to have y-coordinate zero.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::CalcGivensValues(double a, double b, double *c, double *s)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
double denomInv = sqrt(a * a + b * b);
|
|
Toshihiro Shimizu |
890ddd |
if (denomInv == 0.0) {
|
|
Toshihiro Shimizu |
890ddd |
*c = 1.0;
|
|
Toshihiro Shimizu |
890ddd |
*s = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
} else {
|
|
Toshihiro Shimizu |
890ddd |
denomInv = 1.0 / denomInv;
|
|
Toshihiro Shimizu |
890ddd |
*c = a * denomInv;
|
|
Toshihiro Shimizu |
890ddd |
*s = -b * denomInv;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Applies Givens transform to columns i and i+1.
|
|
Toshihiro Shimizu |
890ddd |
// Equivalent to postmultiplying by the matrix
|
|
Toshihiro Shimizu |
890ddd |
// ( c -s )
|
|
Toshihiro Shimizu |
890ddd |
// ( s c )
|
|
Toshihiro Shimizu |
890ddd |
// with non-zero entries in rows i and i+1 and columns i and i+1
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::PostApplyGivens(double c, double s, long idx)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(0 <= idx && idx < NumCols);
|
|
Toshihiro Shimizu |
890ddd |
double *colA = x + idx * NumRows;
|
|
Toshihiro Shimizu |
890ddd |
double *colB = colA + NumRows;
|
|
Toshihiro Shimizu |
890ddd |
for (long i = NumRows; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
double temp = *colA;
|
|
Toshihiro Shimizu |
890ddd |
*colA = (*colA) * c + (*colB) * s;
|
|
Toshihiro Shimizu |
890ddd |
*colB = (*colB) * c - temp * s;
|
|
Toshihiro Shimizu |
890ddd |
colA++;
|
|
Toshihiro Shimizu |
890ddd |
colB++;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Applies Givens transform to columns idx1 and idx2.
|
|
Toshihiro Shimizu |
890ddd |
// Equivalent to postmultiplying by the matrix
|
|
Toshihiro Shimizu |
890ddd |
// ( c -s )
|
|
Toshihiro Shimizu |
890ddd |
// ( s c )
|
|
Toshihiro Shimizu |
890ddd |
// with non-zero entries in rows idx1 and idx2 and columns idx1 and idx2
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::PostApplyGivens(double c, double s, long idx1, long idx2)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(idx1 != idx2 && 0 <= idx1 && idx1 < NumCols && 0 <= idx2 && idx2 < NumCols);
|
|
Toshihiro Shimizu |
890ddd |
double *colA = x + idx1 * NumRows;
|
|
Toshihiro Shimizu |
890ddd |
double *colB = x + idx2 * NumRows;
|
|
Toshihiro Shimizu |
890ddd |
for (long i = NumRows; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
double temp = *colA;
|
|
Toshihiro Shimizu |
890ddd |
*colA = (*colA) * c + (*colB) * s;
|
|
Toshihiro Shimizu |
890ddd |
*colB = (*colB) * c - temp * s;
|
|
Toshihiro Shimizu |
890ddd |
colA++;
|
|
Toshihiro Shimizu |
890ddd |
colB++;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// ********************************************************************************************
|
|
Toshihiro Shimizu |
890ddd |
// Singular value decomposition.
|
|
Toshihiro Shimizu |
890ddd |
// Return othogonal matrices U and V and diagonal matrix with diagonal w such that
|
|
Toshihiro Shimizu |
890ddd |
// (this) = U * Diag(w) * V^T (V^T is V-transpose.)
|
|
Toshihiro Shimizu |
890ddd |
// Diagonal entries have all non-zero entries before all zero entries, but are not
|
|
Toshihiro Shimizu |
890ddd |
// necessarily sorted. (Someday, I will write ComputedSortedSVD that handles
|
|
Toshihiro Shimizu |
890ddd |
// sorting the eigenvalues by magnitude.)
|
|
Toshihiro Shimizu |
890ddd |
// ********************************************************************************************
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::ComputeSVD(MatrixRmn &U, VectorRn &w, MatrixRmn &V) const
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Shinya Kitaoka |
12c444 |
assert(U.NumRows == NumRows && V.NumCols == NumCols && U.NumRows == U.NumCols && V.NumRows == V.NumCols && w.GetLength() == std::min(NumRows, NumCols));
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
double temp = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
VectorRn &superDiag = VectorRn::GetWorkVector(w.GetLength() - 1); // Some extra work space. Will get passed around.
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Choose larger of U, V to hold intermediate results
|
|
Toshihiro Shimizu |
890ddd |
// If U is larger than V, use U to store intermediate results
|
|
Toshihiro Shimizu |
890ddd |
// Otherwise use V. In the latter case, we form the SVD of A transpose,
|
|
Toshihiro Shimizu |
890ddd |
// (which is essentially identical to the SVD of A).
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn *leftMatrix;
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn *rightMatrix;
|
|
Toshihiro Shimizu |
890ddd |
if (NumRows >= NumCols) {
|
|
Toshihiro Shimizu |
890ddd |
U.LoadAsSubmatrix(*this); // Copy A into U
|
|
Toshihiro Shimizu |
890ddd |
leftMatrix = &U;
|
|
Toshihiro Shimizu |
890ddd |
rightMatrix = &V;
|
|
Toshihiro Shimizu |
890ddd |
} else {
|
|
Toshihiro Shimizu |
890ddd |
V.LoadAsSubmatrixTranspose(*this); // Copy A-transpose into V
|
|
Toshihiro Shimizu |
890ddd |
leftMatrix = &V;
|
|
Toshihiro Shimizu |
890ddd |
rightMatrix = &U;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Do the actual work to calculate the SVD
|
|
Toshihiro Shimizu |
890ddd |
// Now matrix has at least as many rows as columns
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
CalcBidiagonal(*leftMatrix, *rightMatrix, w, superDiag);
|
|
Toshihiro Shimizu |
890ddd |
ConvertBidiagToDiagonal(*leftMatrix, *rightMatrix, w, superDiag);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// ************************************************ CalcBidiagonal **************************
|
|
Toshihiro Shimizu |
890ddd |
// Helper routine for SVD computation
|
|
Toshihiro Shimizu |
890ddd |
// U is a matrix to be bidiagonalized.
|
|
Toshihiro Shimizu |
890ddd |
// On return, U and V are orthonormal and w holds the new diagonal
|
|
Toshihiro Shimizu |
890ddd |
// elements and superDiag holds the super diagonal elements.
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::CalcBidiagonal(MatrixRmn &U, MatrixRmn &V, VectorRn &w, VectorRn &superDiag)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
assert(U.NumRows >= V.NumRows);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// The diagonal and superdiagonal entries of the bidiagonalized
|
|
Toshihiro Shimizu |
890ddd |
// version of the U matrix
|
|
Toshihiro Shimizu |
890ddd |
// are stored in the vectors w and superDiag (temporarily).
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Apply Householder transformations to U.
|
|
Toshihiro Shimizu |
890ddd |
// Householder transformations come in pairs.
|
|
Toshihiro Shimizu |
890ddd |
// First, on the left, we map a portion of a column to zeros
|
|
Toshihiro Shimizu |
890ddd |
// Second, on the right, we map a portion of a row to zeros
|
|
Toshihiro Shimizu |
890ddd |
const long rowStep = U.NumCols;
|
|
Toshihiro Shimizu |
890ddd |
const long diagStep = U.NumCols + 1;
|
|
Toshihiro Shimizu |
890ddd |
double *diagPtr = U.x;
|
|
Toshihiro Shimizu |
890ddd |
double *wPtr = w.x;
|
|
Toshihiro Shimizu |
890ddd |
double *superDiagPtr = superDiag.x;
|
|
Toshihiro Shimizu |
890ddd |
long colLengthLeft = U.NumRows;
|
|
Toshihiro Shimizu |
890ddd |
long rowLengthLeft = V.NumCols;
|
|
Toshihiro Shimizu |
890ddd |
while (true) {
|
|
Toshihiro Shimizu |
890ddd |
// Apply a Householder xform on left to zero part of a column
|
|
Toshihiro Shimizu |
890ddd |
SvdHouseholder(diagPtr, colLengthLeft, rowLengthLeft, 1, rowStep, wPtr);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
if (rowLengthLeft == 2) {
|
|
Toshihiro Shimizu |
890ddd |
*superDiagPtr = *(diagPtr + rowStep);
|
|
Toshihiro Shimizu |
890ddd |
break;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Apply a Householder xform on the right to zero part of a row
|
|
Toshihiro Shimizu |
890ddd |
SvdHouseholder(diagPtr + rowStep, rowLengthLeft - 1, colLengthLeft, rowStep, 1, superDiagPtr);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
rowLengthLeft--;
|
|
Toshihiro Shimizu |
890ddd |
colLengthLeft--;
|
|
Toshihiro Shimizu |
890ddd |
diagPtr += diagStep;
|
|
Toshihiro Shimizu |
890ddd |
wPtr++;
|
|
Toshihiro Shimizu |
890ddd |
superDiagPtr++;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
int extra = 0;
|
|
Toshihiro Shimizu |
890ddd |
diagPtr += diagStep;
|
|
Toshihiro Shimizu |
890ddd |
wPtr++;
|
|
Toshihiro Shimizu |
890ddd |
if (colLengthLeft > 2) {
|
|
Toshihiro Shimizu |
890ddd |
extra = 1;
|
|
Toshihiro Shimizu |
890ddd |
// Do one last Householder transformation when the matrix is not square
|
|
Toshihiro Shimizu |
890ddd |
colLengthLeft--;
|
|
Toshihiro Shimizu |
890ddd |
SvdHouseholder(diagPtr, colLengthLeft, 1, 1, 0, wPtr);
|
|
Toshihiro Shimizu |
890ddd |
} else {
|
|
Toshihiro Shimizu |
890ddd |
*wPtr = *diagPtr;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Form U and V from the Householder transformations
|
|
Toshihiro Shimizu |
890ddd |
V.ExpandHouseholders(V.NumCols - 2, 1, U.x + U.NumRows, U.NumRows, 1);
|
|
Toshihiro Shimizu |
890ddd |
U.ExpandHouseholders(V.NumCols - 1 + extra, 0, U.x, 1, U.NumRows);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Done with bidiagonalization
|
|
Toshihiro Shimizu |
890ddd |
return;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Helper routine for CalcBidiagonal
|
|
Toshihiro Shimizu |
890ddd |
// Performs a series of Householder transformations on a matrix
|
|
Toshihiro Shimizu |
890ddd |
// Stores results compactly into the matrix: The Householder vector u (normalized)
|
|
Toshihiro Shimizu |
890ddd |
// is stored into the first row/column being transformed.
|
|
Toshihiro Shimizu |
890ddd |
// The leading term of that row (= plus/minus its magnitude is returned
|
|
Toshihiro Shimizu |
890ddd |
// separately into "retFirstEntry"
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::SvdHouseholder(double *basePt,
|
|
Toshihiro Shimizu |
890ddd |
long colLength, long numCols, long colStride, long rowStride,
|
|
Toshihiro Shimizu |
890ddd |
double *retFirstEntry)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Calc norm of vector u
|
|
Toshihiro Shimizu |
890ddd |
double *cPtr = basePt;
|
|
Toshihiro Shimizu |
890ddd |
double norm = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
long i;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
double aa0 = *cPtr;
|
|
Toshihiro Shimizu |
890ddd |
double aa1 = *basePt;
|
|
Toshihiro Shimizu |
890ddd |
double aa2 = *retFirstEntry;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
for (i = colLength; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
norm += Square(*cPtr);
|
|
Toshihiro Shimizu |
890ddd |
cPtr += colStride;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
norm = sqrt(norm); // Norm of vector to reflect to axis e_1
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Handle sign issues
|
|
Toshihiro Shimizu |
890ddd |
double imageVal; // Choose sign to maximize distance
|
|
Toshihiro Shimizu |
890ddd |
if ((*basePt) < 0.0) {
|
|
Toshihiro Shimizu |
890ddd |
imageVal = norm;
|
|
Toshihiro Shimizu |
890ddd |
norm = 2.0 * norm * (norm - (*basePt));
|
|
Toshihiro Shimizu |
890ddd |
} else {
|
|
Toshihiro Shimizu |
890ddd |
imageVal = -norm;
|
|
Toshihiro Shimizu |
890ddd |
norm = 2.0 * norm * (norm + (*basePt));
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
norm = sqrt(norm); // Norm is norm of reflection vector
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
if (norm == 0.0) { // If the vector being transformed is equal to zero
|
|
Toshihiro Shimizu |
890ddd |
// Force to zero in case of roundoff errors
|
|
Toshihiro Shimizu |
890ddd |
cPtr = basePt;
|
|
Toshihiro Shimizu |
890ddd |
for (i = colLength; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
*cPtr = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
cPtr += colStride;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
*retFirstEntry = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
return;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
*retFirstEntry = imageVal;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Set up the normalized Householder vector
|
|
Toshihiro Shimizu |
890ddd |
*basePt -= imageVal; // First component changes. Rest stay the same.
|
|
Toshihiro Shimizu |
890ddd |
// Normalize the vector
|
|
Toshihiro Shimizu |
890ddd |
norm = 1.0 / norm; // Now it is the inverse norm
|
|
Toshihiro Shimizu |
890ddd |
cPtr = basePt;
|
|
Toshihiro Shimizu |
890ddd |
for (i = colLength; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
*cPtr *= norm;
|
|
Toshihiro Shimizu |
890ddd |
cPtr += colStride;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Transform the rest of the U matrix with the Householder transformation
|
|
Toshihiro Shimizu |
890ddd |
double *rPtr = basePt;
|
|
Toshihiro Shimizu |
890ddd |
for (long j = numCols - 1; j > 0; j--) {
|
|
Toshihiro Shimizu |
890ddd |
rPtr += rowStride;
|
|
Toshihiro Shimizu |
890ddd |
// Calc dot product with Householder transformation vector
|
|
Toshihiro Shimizu |
890ddd |
double dotP = DotArray(colLength, basePt, colStride, rPtr, colStride);
|
|
Toshihiro Shimizu |
890ddd |
// Transform with I - 2*dotP*(Householder vector)
|
|
Toshihiro Shimizu |
890ddd |
AddArrayScale(colLength, basePt, colStride, rPtr, colStride, -2.0 * dotP);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// ********************************* ExpandHouseholders ********************************************
|
|
Toshihiro Shimizu |
890ddd |
// The matrix will be square.
|
|
Toshihiro Shimizu |
890ddd |
// numXforms = number of Householder transformations to concatenate
|
|
Toshihiro Shimizu |
890ddd |
// Each Householder transformation is represented by a unit vector
|
|
Toshihiro Shimizu |
890ddd |
// Each successive Householder transformation starts one position later
|
|
Toshihiro Shimizu |
890ddd |
// and has one more implied leading zero
|
|
Toshihiro Shimizu |
890ddd |
// basePt = beginning of the first Householder transform
|
|
Toshihiro Shimizu |
890ddd |
// colStride, rowStride: Householder xforms are stored in "columns"
|
|
Toshihiro Shimizu |
890ddd |
// numZerosSkipped is the number of implicit zeros on the front each
|
|
Toshihiro Shimizu |
890ddd |
// Householder transformation vector (only values supported are 0 and 1).
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::ExpandHouseholders(long numXforms, int numZerosSkipped, const double *basePt, long colStride, long rowStride)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
// Number of applications of the last Householder transform
|
|
Toshihiro Shimizu |
890ddd |
// (That are not trivial!)
|
|
Toshihiro Shimizu |
890ddd |
long numToTransform = NumCols - numXforms + 1 - numZerosSkipped;
|
|
Toshihiro Shimizu |
890ddd |
assert(numToTransform > 0);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
if (numXforms == 0) {
|
|
Toshihiro Shimizu |
890ddd |
SetIdentity();
|
|
Toshihiro Shimizu |
890ddd |
return;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Handle the first one separately as a special case,
|
|
Toshihiro Shimizu |
890ddd |
// "this" matrix will be treated to simulate being preloaded with the identity
|
|
Toshihiro Shimizu |
890ddd |
long hDiagStride = rowStride + colStride;
|
|
Toshihiro Shimizu |
890ddd |
const double *hBase = basePt + hDiagStride * (numXforms - 1); // Pointer to the last Householder vector
|
|
Toshihiro Shimizu |
890ddd |
const double *hDiagPtr = hBase + colStride * (numToTransform - 1); // Pointer to last entry in that vector
|
|
Toshihiro Shimizu |
890ddd |
long i;
|
|
Toshihiro Shimizu |
890ddd |
double *diagPtr = x + NumCols * NumRows - 1; // Last entry in matrix (points to diagonal entry)
|
|
Toshihiro Shimizu |
890ddd |
double *colPtr = diagPtr - (numToTransform - 1); // Pointer to column in matrix
|
|
Toshihiro Shimizu |
890ddd |
for (i = numToTransform; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
CopyArrayScale(numToTransform, hBase, colStride, colPtr, 1, -2.0 * (*hDiagPtr));
|
|
Toshihiro Shimizu |
890ddd |
*diagPtr += 1.0; // Add back in 1 to the diagonal entry (since xforming the identity)
|
|
Toshihiro Shimizu |
890ddd |
diagPtr -= (NumRows + 1); // Next diagonal entry in this matrix
|
|
Toshihiro Shimizu |
890ddd |
colPtr -= NumRows; // Next column in this matrix
|
|
Toshihiro Shimizu |
890ddd |
hDiagPtr -= colStride;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Now handle the general case
|
|
Toshihiro Shimizu |
890ddd |
// A row of zeros must be in effect added to the top of each old column (in each loop)
|
|
Toshihiro Shimizu |
890ddd |
double *colLastPtr = x + NumRows * NumCols - numToTransform - 1;
|
|
Toshihiro Shimizu |
890ddd |
for (i = numXforms - 1; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
numToTransform++; // Number of non-trivial applications of this Householder transformation
|
|
Toshihiro Shimizu |
890ddd |
hBase -= hDiagStride; // Pointer to the beginning of the Householder transformation
|
|
Toshihiro Shimizu |
890ddd |
colPtr = colLastPtr;
|
|
Toshihiro Shimizu |
890ddd |
for (long j = numToTransform - 1; j > 0; j--) {
|
|
Toshihiro Shimizu |
890ddd |
// Get dot product
|
|
Toshihiro Shimizu |
890ddd |
double dotProd2N = -2.0 * DotArray(numToTransform - 1, hBase + colStride, colStride, colPtr + 1, 1);
|
|
Toshihiro Shimizu |
890ddd |
*colPtr = dotProd2N * (*hBase); // Adding onto zero at initial point
|
|
Toshihiro Shimizu |
890ddd |
AddArrayScale(numToTransform - 1, hBase + colStride, colStride, colPtr + 1, 1, dotProd2N);
|
|
Toshihiro Shimizu |
890ddd |
colPtr -= NumRows;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
// Do last one as a special case (may overwrite the Householder vector)
|
|
Toshihiro Shimizu |
890ddd |
CopyArrayScale(numToTransform, hBase, colStride, colPtr, 1, -2.0 * (*hBase));
|
|
Toshihiro Shimizu |
890ddd |
*colPtr += 1.0; // Add back one one as identity
|
|
Toshihiro Shimizu |
890ddd |
// Done with this Householder transformation
|
|
Toshihiro Shimizu |
890ddd |
colLastPtr--;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
if (numZerosSkipped != 0) {
|
|
Toshihiro Shimizu |
890ddd |
assert(numZerosSkipped == 1);
|
|
Toshihiro Shimizu |
890ddd |
// Fill first row and column with identity (More generally: first numZerosSkipped many rows and columns)
|
|
Toshihiro Shimizu |
890ddd |
double *d = x;
|
|
Toshihiro Shimizu |
890ddd |
*d = 1;
|
|
Toshihiro Shimizu |
890ddd |
double *d2 = d;
|
|
Toshihiro Shimizu |
890ddd |
for (i = NumRows - 1; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
*(++d) = 0;
|
|
Toshihiro Shimizu |
890ddd |
*(d2 += NumRows) = 0;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// **************** ConvertBidiagToDiagonal ***********************************************
|
|
Toshihiro Shimizu |
890ddd |
// Do the iterative transformation from bidiagonal form to diagonal form using
|
|
Toshihiro Shimizu |
890ddd |
// Givens transformation. (Golub-Reinsch)
|
|
Toshihiro Shimizu |
890ddd |
// U and V are square. Size of U less than or equal to that of U.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::ConvertBidiagToDiagonal(MatrixRmn &U, MatrixRmn &V, VectorRn &w, VectorRn &superDiag) const
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
// These two index into the last bidiagonal block (last in the matrix, it will be
|
|
Toshihiro Shimizu |
890ddd |
// first one handled.
|
|
Toshihiro Shimizu |
890ddd |
long lastBidiagIdx = V.NumRows - 1;
|
|
Toshihiro Shimizu |
890ddd |
long firstBidiagIdx = 0;
|
|
Toshihiro Shimizu |
890ddd |
//togliere
|
|
Toshihiro Shimizu |
890ddd |
double aa = w.MaxAbs();
|
|
Toshihiro Shimizu |
890ddd |
double bb = superDiag.MaxAbs();
|
|
Toshihiro Shimizu |
890ddd |
|
|
Shinya Kitaoka |
12c444 |
double eps = 1.0e-15 * std::max(w.MaxAbs(), superDiag.MaxAbs());
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
while (true) {
|
|
Toshihiro Shimizu |
890ddd |
bool workLeft = UpdateBidiagIndices(&firstBidiagIdx, &lastBidiagIdx, w, superDiag, eps);
|
|
Toshihiro Shimizu |
890ddd |
if (!workLeft) {
|
|
Toshihiro Shimizu |
890ddd |
break;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Get ready for first Givens rotation
|
|
Toshihiro Shimizu |
890ddd |
// Push non-zero to M[2,1] with Givens transformation
|
|
Toshihiro Shimizu |
890ddd |
double *wPtr = w.x + firstBidiagIdx;
|
|
Toshihiro Shimizu |
890ddd |
double *sdPtr = superDiag.x + firstBidiagIdx;
|
|
Toshihiro Shimizu |
890ddd |
double extraOffDiag = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
if ((*wPtr) == 0.0) {
|
|
Toshihiro Shimizu |
890ddd |
ClearRowWithDiagonalZero(firstBidiagIdx, lastBidiagIdx, U, wPtr, sdPtr, eps);
|
|
Toshihiro Shimizu |
890ddd |
if (firstBidiagIdx > 0) {
|
|
Toshihiro Shimizu |
890ddd |
if (NearZero(*(--sdPtr), eps)) {
|
|
Toshihiro Shimizu |
890ddd |
*sdPtr = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
} else {
|
|
Toshihiro Shimizu |
890ddd |
ClearColumnWithDiagonalZero(firstBidiagIdx, V, wPtr, sdPtr, eps);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
continue;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Estimate an eigenvalue from bottom four entries of M
|
|
Toshihiro Shimizu |
890ddd |
// This gives a lambda value which will shift the Givens rotations
|
|
Toshihiro Shimizu |
890ddd |
// Last four entries of M^T * M are ( ( A, B ), ( B, C ) ).
|
|
Toshihiro Shimizu |
890ddd |
double A;
|
|
Toshihiro Shimizu |
890ddd |
A = (firstBidiagIdx < lastBidiagIdx - 1) ? Square(superDiag[lastBidiagIdx - 2]) : 0.0;
|
|
Toshihiro Shimizu |
890ddd |
double BSq = Square(w[lastBidiagIdx - 1]);
|
|
Toshihiro Shimizu |
890ddd |
A += BSq; // The "A" entry of M^T * M
|
|
Toshihiro Shimizu |
890ddd |
double C = Square(superDiag[lastBidiagIdx - 1]);
|
|
Toshihiro Shimizu |
890ddd |
BSq *= C; // The squared "B" entry
|
|
Toshihiro Shimizu |
890ddd |
C += Square(w[lastBidiagIdx]); // The "C" entry
|
|
Toshihiro Shimizu |
890ddd |
double lambda; // lambda will hold the estimated eigenvalue
|
|
Toshihiro Shimizu |
890ddd |
lambda = sqrt(Square((A - C) * 0.5) + BSq); // Use the lambda value that is closest to C.
|
|
Toshihiro Shimizu |
890ddd |
if (A > C) {
|
|
Toshihiro Shimizu |
890ddd |
lambda = -lambda;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
lambda += (A + C) * 0.5; // Now lambda equals the estimate for the last eigenvalue
|
|
Toshihiro Shimizu |
890ddd |
double t11 = Square(w[firstBidiagIdx]);
|
|
Toshihiro Shimizu |
890ddd |
double t12 = w[firstBidiagIdx] * superDiag[firstBidiagIdx];
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
double c, s;
|
|
Toshihiro Shimizu |
890ddd |
CalcGivensValues(t11 - lambda, t12, &c, &s);
|
|
Toshihiro Shimizu |
890ddd |
ApplyGivensCBTD(c, s, wPtr, sdPtr, &extraOffDiag, wPtr + 1);
|
|
Toshihiro Shimizu |
890ddd |
V.PostApplyGivens(c, -s, firstBidiagIdx);
|
|
Toshihiro Shimizu |
890ddd |
long i;
|
|
Toshihiro Shimizu |
890ddd |
for (i = firstBidiagIdx; i < lastBidiagIdx - 1; i++) {
|
|
Toshihiro Shimizu |
890ddd |
// Push non-zero from M[i+1,i] to M[i,i+2]
|
|
Toshihiro Shimizu |
890ddd |
CalcGivensValues(*wPtr, extraOffDiag, &c, &s);
|
|
Toshihiro Shimizu |
890ddd |
ApplyGivensCBTD(c, s, wPtr, sdPtr, &extraOffDiag, extraOffDiag, wPtr + 1, sdPtr + 1);
|
|
Toshihiro Shimizu |
890ddd |
U.PostApplyGivens(c, -s, i);
|
|
Toshihiro Shimizu |
890ddd |
// Push non-zero from M[i,i+2] to M[1+2,i+1]
|
|
Toshihiro Shimizu |
890ddd |
CalcGivensValues(*sdPtr, extraOffDiag, &c, &s);
|
|
Toshihiro Shimizu |
890ddd |
ApplyGivensCBTD(c, s, sdPtr, wPtr + 1, &extraOffDiag, extraOffDiag, sdPtr + 1, wPtr + 2);
|
|
Toshihiro Shimizu |
890ddd |
V.PostApplyGivens(c, -s, i + 1);
|
|
Toshihiro Shimizu |
890ddd |
wPtr++;
|
|
Toshihiro Shimizu |
890ddd |
sdPtr++;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
// Push non-zero value from M[i+1,i] to M[i,i+1] for i==lastBidiagIdx-1
|
|
Toshihiro Shimizu |
890ddd |
CalcGivensValues(*wPtr, extraOffDiag, &c, &s);
|
|
Toshihiro Shimizu |
890ddd |
ApplyGivensCBTD(c, s, wPtr, &extraOffDiag, sdPtr, wPtr + 1);
|
|
Toshihiro Shimizu |
890ddd |
U.PostApplyGivens(c, -s, i);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// DEBUG
|
|
Toshihiro Shimizu |
890ddd |
// DebugCalcBidiagCheck( V, w, superDiag, U );
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// This is called when there is a zero diagonal entry, with a non-zero superdiagonal entry on the same row.
|
|
Toshihiro Shimizu |
890ddd |
// We use Givens rotations to "chase" the non-zero entry across the row; when it reaches the last
|
|
Toshihiro Shimizu |
890ddd |
// column, it is finally zeroed away.
|
|
Toshihiro Shimizu |
890ddd |
// wPtr points to the zero entry on the diagonal. sdPtr points to the non-zero superdiagonal entry on the same row.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::ClearRowWithDiagonalZero(long firstBidiagIdx, long lastBidiagIdx, MatrixRmn &U, double *wPtr, double *sdPtr, double eps)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
double curSd = *sdPtr; // Value being chased across the row
|
|
Toshihiro Shimizu |
890ddd |
*sdPtr = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
long i = firstBidiagIdx + 1;
|
|
Toshihiro Shimizu |
890ddd |
while (true) {
|
|
Toshihiro Shimizu |
890ddd |
// Rotate row i and row firstBidiagIdx (Givens rotation)
|
|
Toshihiro Shimizu |
890ddd |
double c, s;
|
|
Toshihiro Shimizu |
890ddd |
CalcGivensValues(*(++wPtr), curSd, &c, &s);
|
|
Toshihiro Shimizu |
890ddd |
U.PostApplyGivens(c, -s, i, firstBidiagIdx);
|
|
Toshihiro Shimizu |
890ddd |
*wPtr = c * (*wPtr) - s * curSd;
|
|
Toshihiro Shimizu |
890ddd |
if (i == lastBidiagIdx) {
|
|
Toshihiro Shimizu |
890ddd |
break;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
curSd = s * (*(++sdPtr)); // New value pops up one column over to the right
|
|
Toshihiro Shimizu |
890ddd |
*sdPtr = c * (*sdPtr);
|
|
Toshihiro Shimizu |
890ddd |
i++;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// This is called when there is a zero diagonal entry, with a non-zero superdiagonal entry in the same column.
|
|
Toshihiro Shimizu |
890ddd |
// We use Givens rotations to "chase" the non-zero entry up the column; when it reaches the last
|
|
Toshihiro Shimizu |
890ddd |
// column, it is finally zeroed away.
|
|
Toshihiro Shimizu |
890ddd |
// wPtr points to the zero entry on the diagonal. sdPtr points to the non-zero superdiagonal entry in the same column.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::ClearColumnWithDiagonalZero(long endIdx, MatrixRmn &V, double *wPtr, double *sdPtr, double eps)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
double curSd = *sdPtr; // Value being chased up the column
|
|
Toshihiro Shimizu |
890ddd |
*sdPtr = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
long i = endIdx - 1;
|
|
Toshihiro Shimizu |
890ddd |
while (true) {
|
|
Toshihiro Shimizu |
890ddd |
double c, s;
|
|
Toshihiro Shimizu |
890ddd |
CalcGivensValues(*(--wPtr), curSd, &c, &s);
|
|
Toshihiro Shimizu |
890ddd |
V.PostApplyGivens(c, -s, i, endIdx);
|
|
Toshihiro Shimizu |
890ddd |
*wPtr = c * (*wPtr) - s * curSd;
|
|
Toshihiro Shimizu |
890ddd |
if (i == 0) {
|
|
Toshihiro Shimizu |
890ddd |
break;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
curSd = s * (*(--sdPtr)); // New value pops up one row above
|
|
Toshihiro Shimizu |
890ddd |
if (NearZero(curSd, eps)) {
|
|
Toshihiro Shimizu |
890ddd |
break;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
*sdPtr = c * (*sdPtr);
|
|
Toshihiro Shimizu |
890ddd |
i--;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Matrix A is ( ( a c ) ( b d ) ), i.e., given in column order.
|
|
Toshihiro Shimizu |
890ddd |
// Mult's G[c,s] times A, replaces A.
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::ApplyGivensCBTD(double cosine, double sine, double *a, double *b, double *c, double *d)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
double temp = *a;
|
|
Toshihiro Shimizu |
890ddd |
*a = cosine * (*a) - sine * (*b);
|
|
Toshihiro Shimizu |
890ddd |
*b = sine * temp + cosine * (*b);
|
|
Toshihiro Shimizu |
890ddd |
temp = *c;
|
|
Toshihiro Shimizu |
890ddd |
*c = cosine * (*c) - sine * (*d);
|
|
Toshihiro Shimizu |
890ddd |
*d = sine * temp + cosine * (*d);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Now matrix A given in row order, A = ( ( a b c ) ( d e f ) ).
|
|
Toshihiro Shimizu |
890ddd |
// Return G[c,s] * A, replace A. d becomes zero, no need to return.
|
|
Toshihiro Shimizu |
890ddd |
// Also, it is certain the old *c value is taken to be zero!
|
|
Toshihiro Shimizu |
890ddd |
void MatrixRmn::ApplyGivensCBTD(double cosine, double sine, double *a, double *b, double *c,
|
|
Toshihiro Shimizu |
890ddd |
double d, double *e, double *f)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
*a = cosine * (*a) - sine * d;
|
|
Toshihiro Shimizu |
890ddd |
double temp = *b;
|
|
Toshihiro Shimizu |
890ddd |
*b = cosine * (*b) - sine * (*e);
|
|
Toshihiro Shimizu |
890ddd |
*e = sine * temp + cosine * (*e);
|
|
Toshihiro Shimizu |
890ddd |
*c = -sine * (*f);
|
|
Toshihiro Shimizu |
890ddd |
*f = cosine * (*f);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Helper routine for SVD conversion from bidiagonal to diagonal
|
|
Toshihiro Shimizu |
890ddd |
bool MatrixRmn::UpdateBidiagIndices(long *firstBidiagIdx, long *lastBidiagIdx, VectorRn &w, VectorRn &superDiag, double eps)
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
long lastIdx = *lastBidiagIdx;
|
|
Toshihiro Shimizu |
890ddd |
double *sdPtr = superDiag.GetPtr(lastIdx - 1); // Entry above the last diagonal entry
|
|
Toshihiro Shimizu |
890ddd |
while (NearZero(*sdPtr, eps)) {
|
|
Toshihiro Shimizu |
890ddd |
*(sdPtr--) = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
lastIdx--;
|
|
Toshihiro Shimizu |
890ddd |
if (lastIdx == 0) {
|
|
Toshihiro Shimizu |
890ddd |
return false;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
*lastBidiagIdx = lastIdx;
|
|
Toshihiro Shimizu |
890ddd |
long firstIdx = lastIdx - 1;
|
|
Toshihiro Shimizu |
890ddd |
double *wPtr = w.GetPtr(firstIdx);
|
|
Toshihiro Shimizu |
890ddd |
while (firstIdx > 0) {
|
|
Toshihiro Shimizu |
890ddd |
if (NearZero(*wPtr, eps)) { // If this diagonal entry (near) zero
|
|
Toshihiro Shimizu |
890ddd |
*wPtr = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
break;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
if (NearZero(*(--sdPtr), eps)) { // If the entry above the diagonal entry is (near) zero
|
|
Toshihiro Shimizu |
890ddd |
*sdPtr = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
break;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
wPtr--;
|
|
Toshihiro Shimizu |
890ddd |
firstIdx--;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
*firstBidiagIdx = firstIdx;
|
|
Toshihiro Shimizu |
890ddd |
return true;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// ******************************************DEBUG STUFFF
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
bool MatrixRmn::DebugCheckSVD(const MatrixRmn &U, const VectorRn &w, const MatrixRmn &V) const
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
// Special SVD test code
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn IV(V.getNumRows(), V.getNumColumns());
|
|
Toshihiro Shimizu |
890ddd |
IV.SetIdentity();
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn VTV(V.getNumRows(), V.getNumColumns());
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn::TransposeMultiply(V, V, VTV);
|
|
Toshihiro Shimizu |
890ddd |
IV -= VTV;
|
|
Toshihiro Shimizu |
890ddd |
double error = IV.FrobeniusNorm();
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn IU(U.getNumRows(), U.getNumColumns());
|
|
Toshihiro Shimizu |
890ddd |
IU.SetIdentity();
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn UTU(U.getNumRows(), U.getNumColumns());
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn::TransposeMultiply(U, U, UTU);
|
|
Toshihiro Shimizu |
890ddd |
IU -= UTU;
|
|
Toshihiro Shimizu |
890ddd |
error += IU.FrobeniusNorm();
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn Diag(U.getNumRows(), V.getNumRows());
|
|
Toshihiro Shimizu |
890ddd |
Diag.SetZero();
|
|
Toshihiro Shimizu |
890ddd |
Diag.SetDiagonalEntries(w);
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn B(U.getNumRows(), V.getNumRows());
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn C(U.getNumRows(), V.getNumRows());
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn::Multiply(U, Diag, B);
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn::MultiplyTranspose(B, V, C);
|
|
Toshihiro Shimizu |
890ddd |
C -= *this;
|
|
Toshihiro Shimizu |
890ddd |
error += C.FrobeniusNorm();
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
bool ret = (fabs(error) <= 1.0e-13 * w.MaxAbs());
|
|
Toshihiro Shimizu |
890ddd |
assert(ret);
|
|
Toshihiro Shimizu |
890ddd |
return ret;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
//=============================================================================
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
const double PI = 3.1415926535897932384626433832795028841972;
|
|
Toshihiro Shimizu |
890ddd |
const double RadiansToDegrees = 180.0 / PI;
|
|
Toshihiro Shimizu |
890ddd |
const double DegreesToRadians = PI / 180;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
const double Jacobian::DefaultDampingLambda = 1.1;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
const double Jacobian::PseudoInverseThresholdFactor = 0.001;
|
|
Toshihiro Shimizu |
890ddd |
const double Jacobian::MaxAngleJtranspose = 30.0 * DegreesToRadians;
|
|
Toshihiro Shimizu |
890ddd |
const double Jacobian::MaxAnglePseudoinverse = 5.0 * DegreesToRadians;
|
|
Toshihiro Shimizu |
890ddd |
const double Jacobian::MaxAngleDLS = 5.0 * DegreesToRadians;
|
|
Toshihiro Shimizu |
890ddd |
const double Jacobian::MaxAngleSDLS = 45.0 * DegreesToRadians;
|
|
Toshihiro Shimizu |
890ddd |
const double Jacobian::BaseMaxTargetDist = 3.4;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
Jacobian::Jacobian(IKSkeleton *skeleton, std::vector<tpointd> &targetPos)</tpointd>
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
Jacobian::skeleton = skeleton;
|
|
Toshihiro Shimizu |
890ddd |
nEffector = skeleton->getNumEffector();
|
|
Toshihiro Shimizu |
890ddd |
nJoint = skeleton->getNodeCount() - nEffector; //numero dei giunti meno gli effettori
|
|
Toshihiro Shimizu |
890ddd |
nRow = 2 * nEffector;
|
|
Toshihiro Shimizu |
890ddd |
nCol = nJoint;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
target = (targetPos);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
Jend.SetSize(nRow, nCol); // Matrice jacobiana
|
|
Toshihiro Shimizu |
890ddd |
Jend.SetZero();
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
Jtarget.SetSize(nRow, nCol); // Matrice jacobiana basta sulle posizioni dei targets (non usata)
|
|
Toshihiro Shimizu |
890ddd |
Jtarget.SetZero();
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
U.SetSize(nRow, nRow); // matrice U per il calcolo SVD
|
|
Toshihiro Shimizu |
890ddd |
w.SetLength(min(nRow, nCol));
|
|
Toshihiro Shimizu |
890ddd |
V.SetSize(nCol, nCol); // matrice V per il calcolo SVD
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
dS.SetLength(nRow); // (Posizione Target ) - (posizione End effector)
|
|
Toshihiro Shimizu |
890ddd |
dTheta.SetLength(nCol); // Cambiamenti degli angoli dei Joints
|
|
Toshihiro Shimizu |
890ddd |
dPreTheta.SetLength(nCol);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Usato nel: metodo del trasposto dello Jacobiano & DLS & SDLS
|
|
Toshihiro Shimizu |
890ddd |
dT.SetLength(nRow);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Usato nel Selectively Damped Least Squares Method
|
|
Toshihiro Shimizu |
890ddd |
dSclamp.SetLength(nEffector);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
Jnorms.SetSize(nEffector, nCol); // Memorizza le norme della matrice attiva J
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
DampingLambdaSqV.SetLength(nRow);
|
|
Toshihiro Shimizu |
890ddd |
diagMatIdentity.SetLength(nCol);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
Reset();
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
void Jacobian::Reset()
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
// Usato nel Damped Least Squares Method
|
|
Toshihiro Shimizu |
890ddd |
DampingLambda = DefaultDampingLambda;
|
|
Toshihiro Shimizu |
890ddd |
DampingLambdaSq = Square(DampingLambda);
|
|
Toshihiro Shimizu |
890ddd |
for (int i = 0; i < DampingLambdaSqV.GetLength(); i++)
|
|
Toshihiro Shimizu |
890ddd |
DampingLambdaSqV[i] = DampingLambdaSq;
|
|
Toshihiro Shimizu |
890ddd |
for (int i = 0; i < diagMatIdentity.GetLength(); i++)
|
|
Toshihiro Shimizu |
890ddd |
diagMatIdentity[i] = 1.0;
|
|
Toshihiro Shimizu |
890ddd |
//DampingLambdaSDLS = 1.5*DefaultDampingLambda;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
dSclamp.Fill(HUGE_VAL);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Calcola il vettore deltaS vector, dS, (l' errore tra end effector e target
|
|
Toshihiro Shimizu |
890ddd |
// Calcola le matrce jacobiana J
|
|
Toshihiro Shimizu |
890ddd |
void Jacobian::computeJacobian()
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
// Scorro tutto lo skeleton per trovare tutti gli end effectors
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
int numNode = skeleton->getNodeCount();
|
|
Toshihiro Shimizu |
890ddd |
for (int index = 0; index < numNode; index++) {
|
|
Toshihiro Shimizu |
890ddd |
IKNode *n = skeleton->getNode(index);
|
|
Toshihiro Shimizu |
890ddd |
int effectorCount = skeleton->getNumEffector();
|
|
Toshihiro Shimizu |
890ddd |
if (n->IsEffector()) {
|
|
Toshihiro Shimizu |
890ddd |
int i = n->getEffectorNum();
|
|
Toshihiro Shimizu |
890ddd |
const TPointD &targetPos = target[i];
|
|
Toshihiro Shimizu |
890ddd |
TPointD temp;
|
|
Toshihiro Shimizu |
890ddd |
// Calcolo i valori di deltaS (differenza tra end effectors e target positions.)
|
|
Toshihiro Shimizu |
890ddd |
temp = targetPos;
|
|
Toshihiro Shimizu |
890ddd |
TPointD a = n->GetS();
|
|
Toshihiro Shimizu |
890ddd |
temp -= n->GetS();
|
|
Toshihiro Shimizu |
890ddd |
if (i < effectorCount - 1) {
|
|
Toshihiro Shimizu |
890ddd |
temp.x = 100 * temp.x;
|
|
Toshihiro Shimizu |
890ddd |
temp.y = 100 * temp.y;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
dS.SetCouple(i, temp);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Find all ancestors (they will usually all be joints)
|
|
Toshihiro Shimizu |
890ddd |
// Set the corresponding entries in the Jacobians J, K.
|
|
Toshihiro Shimizu |
890ddd |
IKNode *m = skeleton->getParent(n);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
while (m) {
|
|
Toshihiro Shimizu |
890ddd |
int j = m->getJointNum();
|
|
Toshihiro Shimizu |
890ddd |
//assert(j>=0 && j<skeleton->GetNumJoint());</skeleton->
|
|
Toshihiro Shimizu |
890ddd |
int numnode = skeleton->getNodeCount();
|
|
Toshihiro Shimizu |
890ddd |
assert(0 <= i && i < nEffector && 0 <= j && j < (skeleton->getNodeCount() - skeleton->getNumEffector()));
|
|
Toshihiro Shimizu |
890ddd |
if (m->isFrozen()) {
|
|
Toshihiro Shimizu |
890ddd |
Jend.SetCouple(i, j, TPointD(0.0, 0.0));
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
} else {
|
|
Toshihiro Shimizu |
890ddd |
temp = m->GetS(); // joint pos.
|
|
Toshihiro Shimizu |
890ddd |
temp -= n->GetS(); // -(end effector pos. - joint pos.)
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
double tx = temp.x;
|
|
Toshihiro Shimizu |
890ddd |
temp.x = temp.y;
|
|
Toshihiro Shimizu |
890ddd |
temp.y = -tx;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
if (i < effectorCount - 1) {
|
|
Toshihiro Shimizu |
890ddd |
temp.x = 100 * temp.x;
|
|
Toshihiro Shimizu |
890ddd |
temp.y = 100 * temp.y;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
Jend.SetCouple(i, j, temp);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
m = skeleton->getParent(m);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// The delta theta values have been computed in dTheta array
|
|
Toshihiro Shimizu |
890ddd |
// Apply the delta theta values to the joints
|
|
Toshihiro Shimizu |
890ddd |
// Nothing is done about joint limits for now.
|
|
Toshihiro Shimizu |
890ddd |
void Jacobian::UpdateThetas()
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
// Update the joint angles
|
|
Toshihiro Shimizu |
890ddd |
for (int index = 0; index < skeleton->getNodeCount(); index++) {
|
|
Toshihiro Shimizu |
890ddd |
IKNode *n = skeleton->getNode(index);
|
|
Toshihiro Shimizu |
890ddd |
if (n->IsJoint()) {
|
|
Toshihiro Shimizu |
890ddd |
int i = n->getJointNum();
|
|
Toshihiro Shimizu |
890ddd |
n->AddToTheta(dTheta[i]);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
// Aggiorno le posizioni dei joint
|
|
Toshihiro Shimizu |
890ddd |
skeleton->compute();
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
bool Jacobian::checkJointsLimit()
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
bool clampingDetected = false;
|
|
Toshihiro Shimizu |
890ddd |
/*
|
|
Toshihiro Shimizu |
890ddd |
Node* n = skeleton->getNode(3);
|
|
Toshihiro Shimizu |
890ddd |
int indexJoint = n->getJointNum();
|
|
Toshihiro Shimizu |
890ddd |
double theta = n->getTheta();
|
|
Toshihiro Shimizu |
890ddd |
double upperLimit = PI;
|
|
Toshihiro Shimizu |
890ddd |
double lowerLimit = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
if(theta >upperLimit || theta
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
if(theta
|
|
Toshihiro Shimizu |
890ddd |
clampingDetected = true;
|
|
Toshihiro Shimizu |
890ddd |
double clampingVar = theta - upperLimit;
|
|
Toshihiro Shimizu |
890ddd |
for(int i=0;i<jactive->getNumRows();i++)</jactive->
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
double tmp = clampingVar*Jactive->Get(i,indexJoint);
|
|
Toshihiro Shimizu |
890ddd |
dS[i] -= tmp;
|
|
Toshihiro Shimizu |
890ddd |
Jactive->Set(i,indexJoint,0.0);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
n->setTheta(upperLimit);
|
|
Toshihiro Shimizu |
890ddd |
diagMatIdentity.Set(indexJoint, 0.0);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
}*/
|
|
Toshihiro Shimizu |
890ddd |
return clampingDetected;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
void Jacobian::ZeroDeltaThetas()
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
dTheta.SetZero();
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Find the delta theta values using inverse Jacobian method
|
|
Toshihiro Shimizu |
890ddd |
// Uses a greedy method to decide scaling factor
|
|
Toshihiro Shimizu |
890ddd |
void Jacobian::CalcDeltaThetasTranspose()
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
const MatrixRmn &J = Jend;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
J.MultiplyTranspose(dS, dTheta);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Scale back the dTheta values greedily
|
|
Toshihiro Shimizu |
890ddd |
J.Multiply(dTheta, dT); // dT = J * dTheta
|
|
Toshihiro Shimizu |
890ddd |
double alpha = Dot(dS, dT) / dT.NormSq();
|
|
Toshihiro Shimizu |
890ddd |
assert(alpha > 0.0);
|
|
Toshihiro Shimizu |
890ddd |
// Also scale back to be have max angle change less than MaxAngleJtranspose
|
|
Toshihiro Shimizu |
890ddd |
double maxChange = dTheta.MaxAbs();
|
|
Toshihiro Shimizu |
890ddd |
double beta = MaxAngleJtranspose / maxChange;
|
|
Toshihiro Shimizu |
890ddd |
dTheta *= min(alpha, beta);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
void Jacobian::CalcDeltaThetasPseudoinverse()
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn &J = const_cast<matrixrmn &="">(Jend);</matrixrmn>
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// costruisco matrice J1
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn J1;
|
|
Toshihiro Shimizu |
890ddd |
J1.SetSize(2, J.getNumColumns());
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
for (int i = 0; i < 2; i++)
|
|
Toshihiro Shimizu |
890ddd |
for (int j = 0; j < J.getNumColumns(); j++)
|
|
Toshihiro Shimizu |
890ddd |
J1.Set(i, j, J.Get(i, j));
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// COSTRUISCO VETTORI ds1 e ds2
|
|
Toshihiro Shimizu |
890ddd |
VectorRn dS1(2);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
for (int i = 0; i < 2; i++)
|
|
Toshihiro Shimizu |
890ddd |
dS1.Set(i, dS.Get(i));
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// calcolo dtheta1
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn U, V;
|
|
Toshihiro Shimizu |
890ddd |
VectorRn w;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
U.SetSize(J1.getNumRows(), J1.getNumRows());
|
|
Toshihiro Shimizu |
890ddd |
w.SetLength(min(J1.getNumRows(), J1.getNumColumns()));
|
|
Toshihiro Shimizu |
890ddd |
V.SetSize(J1.getNumColumns(), J1.getNumColumns());
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
J1.ComputeSVD(U, w, V);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Next line for debugging only
|
|
Toshihiro Shimizu |
890ddd |
assert(J1.DebugCheckSVD(U, w, V));
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Calculate response vector dTheta that is the DLS solution.
|
|
Toshihiro Shimizu |
890ddd |
// Delta target values are the dS values
|
|
Toshihiro Shimizu |
890ddd |
// We multiply by Moore-Penrose pseudo-inverse of the J matrix
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
double pseudoInverseThreshold = PseudoInverseThresholdFactor * w.MaxAbs();
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
long diagLength = w.GetLength();
|
|
Toshihiro Shimizu |
890ddd |
double *wPtr = w.GetPtr();
|
|
Toshihiro Shimizu |
890ddd |
dTheta.SetZero();
|
|
Toshihiro Shimizu |
890ddd |
for (long i = 0; i < diagLength; i++) {
|
|
Toshihiro Shimizu |
890ddd |
double dotProdCol = U.DotProductColumn(dS1, i); // Dot product with i-th column of U
|
|
Toshihiro Shimizu |
890ddd |
double alpha = *(wPtr++);
|
|
Toshihiro Shimizu |
890ddd |
if (fabs(alpha) > pseudoInverseThreshold) {
|
|
Toshihiro Shimizu |
890ddd |
alpha = 1.0 / alpha;
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn::AddArrayScale(V.getNumRows(), V.GetColumnPtr(i), 1, dTheta.GetPtr(), 1, dotProdCol * alpha);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn JcurrentPinv(V.getNumRows(), J1.getNumRows()); // pseudoinversa di J1
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn JProjPre(JcurrentPinv.getNumRows(), J1.getNumColumns()); // Proiezione di J1
|
|
Toshihiro Shimizu |
890ddd |
if (skeleton->getNumEffector() > 1) {
|
|
Toshihiro Shimizu |
890ddd |
// calcolo la pseudoinversa di J1
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn VD(V.getNumRows(), J1.getNumRows()); // matrice del prodotto V*w
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
double *wPtr = w.GetPtr();
|
|
Toshihiro Shimizu |
890ddd |
pseudoInverseThreshold = PseudoInverseThresholdFactor * w.MaxAbs();
|
|
Toshihiro Shimizu |
890ddd |
for (int j = 0; j < VD.getNumColumns(); j++) {
|
|
Toshihiro Shimizu |
890ddd |
double *VPtr = V.GetColumnPtr(j);
|
|
Toshihiro Shimizu |
890ddd |
double alpha = *(wPtr++); // elemento matrice diagonale
|
|
Toshihiro Shimizu |
890ddd |
for (int i = 0; i < V.getNumRows(); i++) {
|
|
Toshihiro Shimizu |
890ddd |
if (fabs(alpha) > pseudoInverseThreshold) {
|
|
Toshihiro Shimizu |
890ddd |
double entry = *(VPtr++);
|
|
Toshihiro Shimizu |
890ddd |
VD.Set(i, j, entry * (1.0 / alpha));
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn::MultiplyTranspose(VD, U, JcurrentPinv);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// calcolo la proiezione J1
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn::Multiply(JcurrentPinv, J1, JProjPre);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
for (int j = 0; j < JProjPre.getNumColumns(); j++)
|
|
Toshihiro Shimizu |
890ddd |
for (int i = 0; i < JProjPre.getNumRows(); i++) {
|
|
Toshihiro Shimizu |
890ddd |
double temp = JProjPre.Get(i, j);
|
|
Toshihiro Shimizu |
890ddd |
JProjPre.Set(i, j, -1.0 * temp);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
JProjPre.AddToDiagonal(diagMatIdentity);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
//task priority strategy
|
|
Toshihiro Shimizu |
890ddd |
for (int i = 1; i < skeleton->getNumEffector(); i++) {
|
|
Toshihiro Shimizu |
890ddd |
// costruisco matrice Jcurrent (Ji)
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn Jcurrent(2, J.getNumColumns());
|
|
Toshihiro Shimizu |
890ddd |
for (int j = 0; j < J.getNumColumns(); j++)
|
|
Toshihiro Shimizu |
890ddd |
for (int k = 0; k < 2; k++)
|
|
Toshihiro Shimizu |
890ddd |
Jcurrent.Set(k, j, J.Get(k + 2 * i, j));
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// costruisco il vettore dScurrent
|
|
Toshihiro Shimizu |
890ddd |
VectorRn dScurrent(2);
|
|
Toshihiro Shimizu |
890ddd |
for (int k = 0; k < 2; k++)
|
|
Toshihiro Shimizu |
890ddd |
dScurrent.Set(k, dS.Get(k + 2 * i));
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Moltiplico Jcurrent per la proiezione di J(i-1)
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn Jdst(Jcurrent.getNumRows(), JProjPre.getNumColumns());
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn::Multiply(Jcurrent, JProjPre, Jdst);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Calcolo la pseudoinversa di Jdst
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn UU(Jdst.getNumRows(), Jdst.getNumRows()), VV(Jdst.getNumColumns(), Jdst.getNumColumns());
|
|
Toshihiro Shimizu |
890ddd |
VectorRn ww(min(Jdst.getNumRows(), Jdst.getNumColumns()));
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
Jdst.ComputeSVD(UU, ww, VV);
|
|
Toshihiro Shimizu |
890ddd |
assert(Jdst.DebugCheckSVD(UU, ww, VV));
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn VVD(VV.getNumRows(), J1.getNumRows()); // matrice V*w
|
|
Toshihiro Shimizu |
890ddd |
VVD.SetZero();
|
|
Toshihiro Shimizu |
890ddd |
pseudoInverseThreshold = PseudoInverseThresholdFactor * ww.MaxAbs();
|
|
Toshihiro Shimizu |
890ddd |
double *wwPtr = ww.GetPtr();
|
|
Toshihiro Shimizu |
890ddd |
for (int j = 0; j < VVD.getNumColumns(); j++) {
|
|
Toshihiro Shimizu |
890ddd |
double *VVPtr = VV.GetColumnPtr(j);
|
|
Toshihiro Shimizu |
890ddd |
double alpha = 50 * (*(wwPtr++)); // elemento matrice diagonale
|
|
Toshihiro Shimizu |
890ddd |
for (int i = 0; i < VV.getNumRows(); i++) {
|
|
Toshihiro Shimizu |
890ddd |
if (fabs(alpha) > 100 * pseudoInverseThreshold) {
|
|
Toshihiro Shimizu |
890ddd |
double entry = *(VVPtr++);
|
|
Toshihiro Shimizu |
890ddd |
VVD.Set(i, j, entry * (1.0 / alpha));
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn JdstPinv(VV.getNumRows(), J1.getNumRows());
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn::MultiplyTranspose(VVD, UU, JdstPinv);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
VectorRn dTemp(J1.getNumRows());
|
|
Toshihiro Shimizu |
890ddd |
Jcurrent.Multiply(dTheta, dTemp);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
VectorRn dTemp2(dScurrent.GetLength());
|
|
Toshihiro Shimizu |
890ddd |
for (int k = 0; k < dScurrent.GetLength(); k++)
|
|
Toshihiro Shimizu |
890ddd |
dTemp2[k] = dScurrent[k] - dTemp[k];
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Moltiplico JdstPinv per dTemp2
|
|
Toshihiro Shimizu |
890ddd |
VectorRn dThetaCurrent(JdstPinv.getNumRows());
|
|
Toshihiro Shimizu |
890ddd |
JdstPinv.Multiply(dTemp2, dThetaCurrent);
|
|
Toshihiro Shimizu |
890ddd |
for (int k = 0; k < dTheta.GetLength(); k++)
|
|
Toshihiro Shimizu |
890ddd |
dTheta[k] += dThetaCurrent[k];
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Infine mi calcolo la pseudoinversa di Jcurrent e quindi la sua proiezione che servirà al passo successivo
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// calcolo la pseudoinversa di Jcurrent
|
|
Toshihiro Shimizu |
890ddd |
Jcurrent.ComputeSVD(U, w, V);
|
|
Toshihiro Shimizu |
890ddd |
assert(Jcurrent.DebugCheckSVD(U, w, V));
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn VD(V.getNumRows(), Jcurrent.getNumRows()); // matrice del prodotto V*w
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
double *wPtr = w.GetPtr();
|
|
Toshihiro Shimizu |
890ddd |
pseudoInverseThreshold = PseudoInverseThresholdFactor * w.MaxAbs();
|
|
Toshihiro Shimizu |
890ddd |
for (int j = 0; j < VVD.getNumColumns(); j++) {
|
|
Toshihiro Shimizu |
890ddd |
double *VPtr = V.GetColumnPtr(j);
|
|
Toshihiro Shimizu |
890ddd |
double alpha = *(wPtr++); // elemento matrice diagonale
|
|
Toshihiro Shimizu |
890ddd |
for (int i = 0; i < V.getNumRows(); i++) {
|
|
Toshihiro Shimizu |
890ddd |
if (fabs(alpha) > pseudoInverseThreshold) {
|
|
Toshihiro Shimizu |
890ddd |
double entry = *(VPtr++);
|
|
Toshihiro Shimizu |
890ddd |
VD.Set(i, j, entry * (1.0 / alpha));
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn::MultiplyTranspose(VD, U, JcurrentPinv);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// calcolo la proiezione Jcurrent
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn::Multiply(JcurrentPinv, Jcurrent, JProjPre);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
for (int j = 0; j < JProjPre.getNumColumns(); j++)
|
|
Toshihiro Shimizu |
890ddd |
for (int k = 0; k < JProjPre.getNumRows(); k++) {
|
|
Toshihiro Shimizu |
890ddd |
double temp = JProjPre.Get(k, j);
|
|
Toshihiro Shimizu |
890ddd |
JProjPre.Set(k, j, -1.0 * temp);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
JProjPre.AddToDiagonal(diagMatIdentity);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
//sw.stop();
|
|
Toshihiro Shimizu |
890ddd |
//std::ofstream os("C:\\buttami.txt", std::ios::app);
|
|
Toshihiro Shimizu |
890ddd |
//sw.print(os);
|
|
Toshihiro Shimizu |
890ddd |
//os.close();
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Scale back to not exceed maximum angle changes
|
|
Toshihiro Shimizu |
890ddd |
double maxChange = 10 * dTheta.MaxAbs();
|
|
Toshihiro Shimizu |
890ddd |
if (maxChange > MaxAnglePseudoinverse) {
|
|
Toshihiro Shimizu |
890ddd |
dTheta *= MaxAnglePseudoinverse / maxChange;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
void Jacobian::CalcDeltaThetasDLS()
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
const MatrixRmn &J = Jend;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn::MultiplyTranspose(J, J, U); // U = J * (J^T)
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
U.AddToDiagonal(DampingLambdaSqV);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Use the next four lines instead of the succeeding two lines for the DLS method with clamped error vector e.
|
|
Toshihiro Shimizu |
890ddd |
// CalcdTClampedFromdS();
|
|
Toshihiro Shimizu |
890ddd |
// VectorRn dTextra(2*nEffector);
|
|
Toshihiro Shimizu |
890ddd |
// U.Solve( dT, &dTextra );
|
|
Toshihiro Shimizu |
890ddd |
// J.MultiplyTranspose( dTextra, dTheta );
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Use these two lines for the traditional DLS method
|
|
Toshihiro Shimizu |
890ddd |
// gennaro
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
U.Solve(dS, &dT);
|
|
Toshihiro Shimizu |
890ddd |
J.MultiplyTranspose(dT, dTheta);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Scalo per non superare l'nagolo massimod i cambiamento
|
|
Toshihiro Shimizu |
890ddd |
double maxChange = 100 * dTheta.MaxAbs();
|
|
Toshihiro Shimizu |
890ddd |
if (maxChange > MaxAngleDLS) {
|
|
Toshihiro Shimizu |
890ddd |
dTheta *= MaxAngleDLS / maxChange;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
void Jacobian::CalcDeltaThetasDLSwithSVD()
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
const MatrixRmn &J = Jend;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
J.ComputeSVD(U, w, V);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// For Debug
|
|
Toshihiro Shimizu |
890ddd |
assert(J.DebugCheckSVD(U, w, V));
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Calculate response vector dTheta that is the DLS solution.
|
|
Toshihiro Shimizu |
890ddd |
// Delta target values are the dS values
|
|
Toshihiro Shimizu |
890ddd |
// We multiply by DLS inverse of the J matrix
|
|
Toshihiro Shimizu |
890ddd |
long diagLength = w.GetLength();
|
|
Toshihiro Shimizu |
890ddd |
double *wPtr = w.GetPtr();
|
|
Toshihiro Shimizu |
890ddd |
dTheta.SetZero();
|
|
Toshihiro Shimizu |
890ddd |
for (long i = 0; i < diagLength; i++) {
|
|
Toshihiro Shimizu |
890ddd |
double dotProdCol = U.DotProductColumn(dS, i); // Dot product with i-th column of U
|
|
Toshihiro Shimizu |
890ddd |
double alpha = *(wPtr++);
|
|
Toshihiro Shimizu |
890ddd |
alpha = alpha / (Square(alpha) + DampingLambdaSq);
|
|
Toshihiro Shimizu |
890ddd |
MatrixRmn::AddArrayScale(V.getNumRows(), V.GetColumnPtr(i), 1, dTheta.GetPtr(), 1, dotProdCol * alpha);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Scale back to not exceed maximum angle changes
|
|
Toshihiro Shimizu |
890ddd |
double maxChange = dTheta.MaxAbs();
|
|
Toshihiro Shimizu |
890ddd |
if (maxChange > MaxAngleDLS) {
|
|
Toshihiro Shimizu |
890ddd |
dTheta *= MaxAngleDLS / maxChange;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
void Jacobian::CalcDeltaThetasSDLS()
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
const MatrixRmn &J = Jend;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Compute Singular Value Decomposition
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
J.ComputeSVD(U, w, V);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Next line for debugging only
|
|
Toshihiro Shimizu |
890ddd |
assert(J.DebugCheckSVD(U, w, V));
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Calculate response vector dTheta that is the SDLS solution.
|
|
Toshihiro Shimizu |
890ddd |
// Delta target values are the dS values
|
|
Toshihiro Shimizu |
890ddd |
int nRows = J.getNumRows();
|
|
Toshihiro Shimizu |
890ddd |
int numEndEffectors = skeleton->getNumEffector(); // Equals the number of rows of J divided by three
|
|
Toshihiro Shimizu |
890ddd |
int nCols = J.getNumColumns();
|
|
Toshihiro Shimizu |
890ddd |
dTheta.SetZero();
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Calculate the norms of the 3-vectors in the Jacobian
|
|
Toshihiro Shimizu |
890ddd |
long i;
|
|
Toshihiro Shimizu |
890ddd |
const double *jx = J.GetPtr();
|
|
Toshihiro Shimizu |
890ddd |
double *jnx = Jnorms.GetPtr();
|
|
Toshihiro Shimizu |
890ddd |
for (i = nCols * numEndEffectors; i > 0; i--) {
|
|
Toshihiro Shimizu |
890ddd |
double accumSq = Square(*(jx++));
|
|
Toshihiro Shimizu |
890ddd |
accumSq += Square(*(jx++));
|
|
Toshihiro Shimizu |
890ddd |
accumSq += Square(*(jx++));
|
|
Toshihiro Shimizu |
890ddd |
*(jnx++) = sqrt(accumSq);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Clamp the dS values
|
|
Toshihiro Shimizu |
890ddd |
CalcdTClampedFromdS();
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Loop over each singular vector
|
|
Toshihiro Shimizu |
890ddd |
for (i = 0; i < nRows; i++) {
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
double wiInv = w[i];
|
|
Toshihiro Shimizu |
890ddd |
if (NearZero(wiInv, 1.0e-10)) {
|
|
Toshihiro Shimizu |
890ddd |
continue;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
wiInv = 1.0 / wiInv;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
double N = 0.0; // N is the quasi-1-norm of the i-th column of U
|
|
Toshihiro Shimizu |
890ddd |
double alpha = 0.0; // alpha is the dot product of dT and the i-th column of U
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
const double *dTx = dT.GetPtr();
|
|
Toshihiro Shimizu |
890ddd |
const double *ux = U.GetColumnPtr(i);
|
|
Toshihiro Shimizu |
890ddd |
long j;
|
|
Toshihiro Shimizu |
890ddd |
for (j = numEndEffectors; j > 0; j--) {
|
|
Toshihiro Shimizu |
890ddd |
double tmp;
|
|
Toshihiro Shimizu |
890ddd |
alpha += (*ux) * (*(dTx++));
|
|
Toshihiro Shimizu |
890ddd |
tmp = Square(*(ux++));
|
|
Toshihiro Shimizu |
890ddd |
alpha += (*ux) * (*(dTx++));
|
|
Toshihiro Shimizu |
890ddd |
tmp += Square(*(ux++));
|
|
Toshihiro Shimizu |
890ddd |
alpha += (*ux) * (*(dTx++));
|
|
Toshihiro Shimizu |
890ddd |
tmp += Square(*(ux++));
|
|
Toshihiro Shimizu |
890ddd |
N += sqrt(tmp);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// M is the quasi-1-norm of the response to angles changing according to the i-th column of V
|
|
Toshihiro Shimizu |
890ddd |
// Then is multiplied by the wiInv value.
|
|
Toshihiro Shimizu |
890ddd |
double M = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
double *vx = V.GetColumnPtr(i);
|
|
Toshihiro Shimizu |
890ddd |
jnx = Jnorms.GetPtr();
|
|
Toshihiro Shimizu |
890ddd |
for (j = nCols; j > 0; j--) {
|
|
Toshihiro Shimizu |
890ddd |
double accum = 0.0;
|
|
Toshihiro Shimizu |
890ddd |
for (long k = numEndEffectors; k > 0; k--) {
|
|
Toshihiro Shimizu |
890ddd |
accum += *(jnx++);
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
M += fabs((*(vx++))) * accum;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
M *= fabs(wiInv);
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
double gamma = MaxAngleSDLS;
|
|
Toshihiro Shimizu |
890ddd |
if (N < M) {
|
|
Toshihiro Shimizu |
890ddd |
gamma *= N / M; // Scale back maximum permissable joint angle
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Calculate the dTheta from pure pseudoinverse considerations
|
|
Toshihiro Shimizu |
890ddd |
double scale = alpha * wiInv; // This times i-th column of V is the psuedoinverse response
|
|
Toshihiro Shimizu |
890ddd |
dPreTheta.LoadScaled(V.GetColumnPtr(i), scale);
|
|
Toshihiro Shimizu |
890ddd |
// Now rescale the dTheta values.
|
|
Toshihiro Shimizu |
890ddd |
double max = dPreTheta.MaxAbs();
|
|
Toshihiro Shimizu |
890ddd |
double rescale = (gamma) / (gamma + max);
|
|
Toshihiro Shimizu |
890ddd |
dTheta.AddScaled(dPreTheta, rescale);
|
|
Toshihiro Shimizu |
890ddd |
/*if ( gamma
|
|
Toshihiro Shimizu |
890ddd |
dTheta.AddScaled( dPreTheta, gamma/max );
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
else {
|
|
Toshihiro Shimizu |
890ddd |
dTheta += dPreTheta;
|
|
Toshihiro Shimizu |
890ddd |
}*/
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Scale back to not exceed maximum angle changes
|
|
Toshihiro Shimizu |
890ddd |
double maxChange = dTheta.MaxAbs();
|
|
Toshihiro Shimizu |
890ddd |
if (maxChange > 100 * MaxAngleSDLS) {
|
|
Toshihiro Shimizu |
890ddd |
dTheta *= MaxAngleSDLS / (MaxAngleSDLS + maxChange);
|
|
Toshihiro Shimizu |
890ddd |
//dTheta *= MaxAngleSDLS/maxChange;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
void Jacobian::CalcdTClampedFromdS()
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
long len = dS.GetLength();
|
|
Toshihiro Shimizu |
890ddd |
long j = 0;
|
|
Toshihiro Shimizu |
890ddd |
for (long i = 0; i < len; i += 2, j++) {
|
|
Toshihiro Shimizu |
890ddd |
double normSq = Square(dS[i]) + Square(dS[i + 1]); //+Square(dS[i+2]);
|
|
Toshihiro Shimizu |
890ddd |
if (normSq > Square(dSclamp[j])) {
|
|
Toshihiro Shimizu |
890ddd |
double factor = dSclamp[j] / sqrt(normSq);
|
|
Toshihiro Shimizu |
890ddd |
dT[i] = dS[i] * factor;
|
|
Toshihiro Shimizu |
890ddd |
dT[i + 1] = dS[i + 1] * factor;
|
|
Toshihiro Shimizu |
890ddd |
//dT[i+2] = dS[i+2]*factor;
|
|
Toshihiro Shimizu |
890ddd |
} else {
|
|
Toshihiro Shimizu |
890ddd |
dT[i] = dS[i];
|
|
Toshihiro Shimizu |
890ddd |
dT[i + 1] = dS[i + 1];
|
|
Toshihiro Shimizu |
890ddd |
//dT[i+2] = dS[i+2];
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
void Jacobian::UpdatedSClampValue()
|
|
Toshihiro Shimizu |
890ddd |
{
|
|
Toshihiro Shimizu |
890ddd |
// Traverse skeleton to find all end effectors
|
|
Toshihiro Shimizu |
890ddd |
TPointD temp;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
int numNode = skeleton->getNodeCount();
|
|
Toshihiro Shimizu |
890ddd |
for (int i = 0; i < numNode; i++) {
|
|
Toshihiro Shimizu |
890ddd |
IKNode *n = skeleton->getNode(i);
|
|
Toshihiro Shimizu |
890ddd |
if (n->IsEffector()) {
|
|
Toshihiro Shimizu |
890ddd |
int i = n->getEffectorNum();
|
|
Toshihiro Shimizu |
890ddd |
const TPointD &targetPos = target[i];
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
// Compute the delta S value (differences from end effectors to target positions.
|
|
Toshihiro Shimizu |
890ddd |
// While we are at it, also update the clamping values in dSclamp;
|
|
Toshihiro Shimizu |
890ddd |
temp = targetPos;
|
|
Toshihiro Shimizu |
890ddd |
temp -= n->GetS();
|
|
Toshihiro Shimizu |
890ddd |
double normSi = sqrt(Square(dS[i]) + Square(dS[i + 1]));
|
|
Toshihiro Shimizu |
890ddd |
double norma = sqrt(temp.x * temp.x + temp.y * temp.y);
|
|
Toshihiro Shimizu |
890ddd |
double changedDist = norma - normSi;
|
|
Toshihiro Shimizu |
890ddd |
|
|
Toshihiro Shimizu |
890ddd |
if (changedDist > 0.0) {
|
|
Toshihiro Shimizu |
890ddd |
dSclamp[i] = BaseMaxTargetDist + changedDist;
|
|
Toshihiro Shimizu |
890ddd |
} else {
|
|
Toshihiro Shimizu |
890ddd |
dSclamp[i] = BaseMaxTargetDist;
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|
|
Toshihiro Shimizu |
890ddd |
}
|