patr-schm / TinyAD

Automatic Differentiation in Geometry Processing Made Simple
MIT License
371 stars 17 forks source link

How to use k at run time in dense mode? #6

Closed pressureless closed 1 year ago

pressureless commented 1 year ago

The paper said "The restriction of having to choose k at compile time can be lifted by using TinyAD::Scalar<Eigen::Dynamic, double> at a run time cost."

I tried to modify angle.cc as follows but got an assertion error.

How can I use k at run time? Thank you.

void test(Eigen::VectorXd variable, Eigen::VectorXd y){
    // Choose autodiff scalar type for 3 variables
    using ADouble = TinyAD::Scalar<Eigen::Dynamic, double>;

    // Init a vector of active variables
    auto x = ADouble::make_active(variable);

    // Compute angle between the two vectors
    ADouble angle = acos(x.dot(y) / (x.norm() * y.norm()));

    // Retreive gradient and Hessian w.r.t. x
    Eigen::Vector3d g = angle.grad;
    Eigen::Matrix3d H = angle.Hess;

    TINYAD_INFO("angle: " << std::endl << angle.val);
    TINYAD_INFO("g: " << std::endl << g);
    TINYAD_INFO("H: " << std::endl << H);
}

int main()
{
    Eigen::VectorXd variable(3);
    variable << 0.0, -1.0, 1.0;
    Eigen::VectorXd constant(3);
    constant << 2.0, 3.0, 5.0;
    test(variable, constant);

    return 0;
}
image image
patr-schm commented 1 year ago

Hi, thanks for pointing this out. Dynamic mode was broken and I just pushed a fix.

Note, however, that dynamic mode (in its current implementation) is pretty experimental and way more restricted than static mode. Especially:

Here's a version of your code snippet that works now:

#include <TinyAD/Scalar.hh>

void test(Eigen::VectorXd& variables, Eigen::VectorXd& y)
{
    // Choose number of variables k at run time
    TINYAD_ASSERT_EQ(variables.size(), y.size());
    const int k_dynamic = variables.size();

    // Autodiff scalar type
    using ADouble = TinyAD::Scalar<Eigen::Dynamic, double>;

    // Init a vector of active variables. Derivatives will have dimension k_dynamic.
    Eigen::VectorX<ADouble> x = ADouble::make_active(variables);

    // Compute angle between the two vectors
    // Sadly, Eigen's dot(..) is not available in dynamic mode
    ADouble dot_xy = ADouble::make_passive(0.0, k_dynamic);
    ADouble norm_x = ADouble::make_passive(0.0, k_dynamic);
    for (int i = 0; i < k_dynamic; ++i)
    {
        dot_xy += x[i] * y[i];
        norm_x += x[i] * x[i];
    }
    norm_x = sqrt(norm_x);

    ADouble angle = acos(dot_xy / (norm_x * y.norm()));

    // Retreive gradient and Hessian w.r.t. x
    Eigen::Vector3d g = angle.grad;
    Eigen::Matrix3d H = angle.Hess;

    TINYAD_INFO("angle: " << std::endl << angle.val);
    TINYAD_INFO("g: " << std::endl << g);
    TINYAD_INFO("H: " << std::endl << H);
}

int main()
{
    Eigen::VectorXd variable(3);
    variable << 0.0, -1.0, 1.0;
    Eigen::VectorXd constant(3);
    constant << 2.0, 3.0, 5.0;
    test(variable, constant);

    return 0;
}
pressureless commented 1 year ago

Thank you!