bullet3/test/InverseDynamics/test_invdyn_jacobian.cpp

// Test of kinematic consistency: check if finite differences of velocities, accelerations
// match positions

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <iostream>

#include <gtest/gtest.h>

#include "Bullet3Common/b3Random.h"

#include "CloneTreeCreator.hpp"
#include "CoilCreator.hpp"
#include "DillCreator.hpp"
#include "RandomTreeCreator.hpp"
#include "BulletInverseDynamics/MultiBodyTree.hpp"
#include "MultiBodyTreeDebugGraph.hpp"

using namespace btInverseDynamics;

#if (defined BT_ID_HAVE_MAT3X) && (defined BT_ID_WITH_JACOBIANS)
// minimal smart pointer to make this work for c++2003
template <typename T>
class ptr {
    ptr();
    ptr(const ptr&);
public:
    ptr(T* p) : m_p(p) {};
    ~ptr() { delete m_p; }
    T& operator*() { return *m_p; }
    T* operator->() { return m_p; }
    T*get() {return m_p;}
    const T*get() const {return m_p;}
    friend bool operator==(const ptr<T>& lhs, const ptr<T>& rhs) { return rhs.m_p == lhs.m_p; }
    friend bool operator!=(const ptr<T>& lhs, const ptr<T>& rhs) { return !(rhs.m_p == lhs.m_p);
}

private:
    T* m_p;
};

void calculateDotJacUError(const MultiBodyTreeCreator& creator, const int nloops,
                           double* max_error) {
    // tree1 is used as reference to compute dot(Jacobian)*u from acceleration(dot(u)=0)
    ptr<MultiBodyTree> tree1(CreateMultiBodyTree(creator));
    ASSERT_TRUE(0x0 != tree1);
    CloneTreeCreator clone(tree1.get());
    // tree2 is used to compute dot(Jacobian)*u using the calculateJacobian function
    ptr<MultiBodyTree> tree2(CreateMultiBodyTree(clone));
    ASSERT_TRUE(0x0 != tree2);

    const int ndofs = tree1->numDoFs();
    const int nbodies = tree1->numBodies();
    if (ndofs <= 0) {
        *max_error = 0;
        return;
    }

    vecx q(ndofs);
    vecx u(ndofs);
    vecx dot_u(ndofs);
    vecx zero(ndofs);
    setZero(zero);

    double max_lin_error = 0;
    double max_ang_error = 0;

    for (int loop = 0; loop < nloops; loop++) {
        for (int i = 0; i < q.size(); i++) {
            q(i) = b3RandRange(-B3_PI, B3_PI);
            u(i) = b3RandRange(-B3_PI, B3_PI);
        }

        EXPECT_EQ(0, tree1->calculateKinematics(q, u, zero));
        EXPECT_EQ(0, tree2->calculatePositionAndVelocityKinematics(q, u));
        EXPECT_EQ(0, tree2->calculateJacobians(q, u));

        for (int idx = 0; idx < nbodies; idx++) {
            vec3 tmp1, tmp2;
            vec3 diff;
            EXPECT_EQ(0, tree1->getBodyLinearAcceleration(idx, &tmp1));
            EXPECT_EQ(0, tree2->getBodyDotJacobianTransU(idx, &tmp2));
            diff = tmp1 - tmp2;
            double lin_error = maxAbs(diff);

            if (lin_error > max_lin_error) {
                max_lin_error = lin_error;
            }

            EXPECT_EQ(0, tree1->getBodyAngularAcceleration(idx, &tmp1));
            EXPECT_EQ(0, tree2->getBodyDotJacobianRotU(idx, &tmp2));
            diff = tmp1 - tmp2;
            double ang_error = maxAbs(diff);
            if (ang_error > max_ang_error) {
                max_ang_error = ang_error;
            }
        }
    }
    *max_error = max_ang_error > max_lin_error ? max_ang_error : max_lin_error;
}

void calculateJacobianError(const MultiBodyTreeCreator& creator, const int nloops,
                            double* max_error) {
    // tree1 is used as reference to compute the Jacobian from velocities with unit u vectors.
    ptr<MultiBodyTree> tree1(CreateMultiBodyTree(creator));
    ASSERT_TRUE(0x0 != tree1);
    // tree2 is used to compute the Jacobians using the calculateJacobian function
    CloneTreeCreator clone(tree1.get());
    ptr<MultiBodyTree> tree2(CreateMultiBodyTree(clone));
    ASSERT_TRUE(0x0 != tree2);

    const int ndofs = tree1->numDoFs();
    const int nbodies = tree1->numBodies();

    if (ndofs <= 0) {
        *max_error = 0;
        return;
    }

    vecx q(ndofs);
    vecx zero(ndofs);
    setZero(zero);
    vecx one(ndofs);

    double max_lin_error = 0;
    double max_ang_error = 0;

    for (int loop = 0; loop < nloops; loop++) {
        for (int i = 0; i < q.size(); i++) {
            q(i) = b3RandRange(-B3_PI, B3_PI);
        }

        EXPECT_EQ(0, tree2->calculatePositionKinematics(q));
        EXPECT_EQ(0, tree2->calculateJacobians(q));

        for (int idx = 0; idx < nbodies; idx++) {
            mat3x ref_jac_r(3, ndofs);
            mat3x ref_jac_t(3, ndofs);
            ref_jac_r.setZero();
            ref_jac_t.setZero();
            // this re-computes all jacobians for every body ...
            // but avoids std::vector<eigen matrix> issues
            for (int col = 0; col < ndofs; col++) {
                setZero(one);
                one(col) = 1.0;
                EXPECT_EQ(0, tree1->calculatePositionAndVelocityKinematics(q, one));
                vec3 vel, omg;
                EXPECT_EQ(0, tree1->getBodyLinearVelocity(idx, &vel));
                EXPECT_EQ(0, tree1->getBodyAngularVelocity(idx, &omg));
                setMat3xElem(0, col, omg(0), &ref_jac_r);
                setMat3xElem(1, col, omg(1), &ref_jac_r);
                setMat3xElem(2, col, omg(2), &ref_jac_r);
                setMat3xElem(0, col, vel(0), &ref_jac_t);
                setMat3xElem(1, col, vel(1), &ref_jac_t);
                setMat3xElem(2, col, vel(2), &ref_jac_t);
            }

            mat3x jac_r(3, ndofs);
            mat3x jac_t(3, ndofs);
            mat3x diff(3, ndofs);

            EXPECT_EQ(0, tree2->getBodyJacobianTrans(idx, &jac_t));
            EXPECT_EQ(0, tree2->getBodyJacobianRot(idx, &jac_r));
            sub(ref_jac_t,jac_t,&diff);
            double lin_error = maxAbsMat3x(diff);
            if (lin_error > max_lin_error) {
                max_lin_error = lin_error;
            }
            sub(ref_jac_r, jac_r,&diff);
            double ang_error = maxAbsMat3x(diff);
            if (ang_error > max_ang_error) {
                max_ang_error = ang_error;
            }
        }
    }
    *max_error = max_ang_error > max_lin_error ? max_ang_error : max_lin_error;
}

void calculateVelocityJacobianError(const MultiBodyTreeCreator& creator, const int nloops,
                                    double* max_error) {
    // tree1 is used as reference to compute the velocities directly
    ptr<MultiBodyTree> tree1(CreateMultiBodyTree(creator));
    ASSERT_TRUE(0x0 != tree1);
    // tree2 is used to compute the velocities via jacobians
    CloneTreeCreator clone(tree1.get());
    ptr<MultiBodyTree> tree2(CreateMultiBodyTree(clone));
    ASSERT_TRUE(0x0 != tree2);

    const int ndofs = tree1->numDoFs();
    const int nbodies = tree1->numBodies();

    if (ndofs <= 0) {
        *max_error = 0;
        return;
    }

    vecx q(ndofs);
    vecx u(ndofs);

    double max_lin_error = 0;
    double max_ang_error = 0;

    for (int loop = 0; loop < nloops; loop++) {
        for (int i = 0; i < q.size(); i++) {
            q(i) = b3RandRange(-B3_PI, B3_PI);
            u(i) = b3RandRange(-B3_PI, B3_PI);
        }

        EXPECT_EQ(0, tree1->calculatePositionAndVelocityKinematics(q, u));
        EXPECT_EQ(0, tree2->calculatePositionKinematics(q));
        EXPECT_EQ(0, tree2->calculateJacobians(q));

        for (int idx = 0; idx < nbodies; idx++) {
            vec3 vel1;
            vec3 omg1;
            vec3 vel2;
            vec3 omg2;
            mat3x jac_r2(3, ndofs);
            mat3x jac_t2(3, ndofs);

            EXPECT_EQ(0, tree1->getBodyLinearVelocity(idx, &vel1));
            EXPECT_EQ(0, tree1->getBodyAngularVelocity(idx, &omg1));
            EXPECT_EQ(0, tree2->getBodyJacobianTrans(idx, &jac_t2));
            EXPECT_EQ(0, tree2->getBodyJacobianRot(idx, &jac_r2));
            omg2 = jac_r2 * u;
            vel2 = jac_t2 * u;

            double lin_error = maxAbs(vel1 - vel2);
            if (lin_error > max_lin_error) {
                max_lin_error = lin_error;
            }
            double ang_error = maxAbs(omg1 - omg2);
            if (ang_error > max_ang_error) {
                max_ang_error = ang_error;
            }
        }
    }
    *max_error = max_ang_error > max_lin_error ? max_ang_error : max_lin_error;
}

// test nonlinear terms: dot(Jacobian)*u (linear and angular acceleration for dot_u ==0)
// from Jacobian calculation method and pseudo-numerically using via the kinematics method.
TEST(InvDynJacobians, JacDotJacU) {
    const int kNumLevels = 5;
#ifdef B3_USE_DOUBLE_PRECISION
	const double kMaxError = 1e-12;
#else
    const double kMaxError = 5e-5;
#endif
    const int kNumLoops = 20;
    for (int level = 0; level < kNumLevels; level++) {
        const int nbodies = BT_ID_POW(2, level);
        CoilCreator coil(nbodies);
        double error;
        calculateDotJacUError(coil, kNumLoops, &error);
        EXPECT_GT(kMaxError, error);
        DillCreator dill(level);
        calculateDotJacUError(dill, kNumLoops, &error);
        EXPECT_GT(kMaxError, error);
    }

    const int kRandomLoops = 100;
    const int kMaxRandomBodies = 128;
    for (int loop = 0; loop < kRandomLoops; loop++) {
        RandomTreeCreator random(kMaxRandomBodies);
        double error;
        calculateDotJacUError(random, kNumLoops, &error);
        EXPECT_GT(kMaxError, error);
    }
}

// Jacobians: linear and angular acceleration for dot_u ==0
// from Jacobian calculation method and pseudo-numerically using via the kinematics method.
TEST(InvDynJacobians, Jacobians) {
    const int kNumLevels = 5;
#ifdef B3_USE_DOUBLE_PRECISION
	const double kMaxError = 1e-12;
#else
	const double kMaxError = 5e-5;
#endif
	const int kNumLoops = 20;
    for (int level = 0; level < kNumLevels; level++) {
        const int nbodies = BT_ID_POW(2, level);
        CoilCreator coil(nbodies);
        double error;
        calculateJacobianError(coil, kNumLoops, &error);
        EXPECT_GT(kMaxError, error);
        DillCreator dill(level);
        calculateDotJacUError(dill, kNumLoops, &error);
        EXPECT_GT(kMaxError, error);
    }
    const int kRandomLoops = 20;
    const int kMaxRandomBodies = 16;
    for (int loop = 0; loop < kRandomLoops; loop++) {
        RandomTreeCreator random(kMaxRandomBodies);
        double error;
        calculateJacobianError(random, kNumLoops, &error);
        EXPECT_GT(kMaxError, error);
    }
}

// test for jacobian*u == velocity
TEST(InvDynJacobians, VelocitiesFromJacobians) {
    const int kRandomLoops = 20;
    const int kMaxRandomBodies = 16;
    const int kNumLoops = 20;
#ifdef B3_USE_DOUBLE_PRECISION
	const double kMaxError = 1e-12;
#else
	const double kMaxError = 5e-5;
#endif
	for (int loop = 0; loop < kRandomLoops; loop++) {
        RandomTreeCreator random(kMaxRandomBodies);
        double error;
        calculateVelocityJacobianError(random, kNumLoops, &error);
        EXPECT_GT(kMaxError, error);
    }
}
#endif

int main(int argc, char** argv) {
	b3Srand(1234);
    ::testing::InitGoogleTest(&argc, argv);
    return RUN_ALL_TESTS();
}