search

rjhogan / Adept-2

Combined array and automatic differentiation library in C++
http://www.met.reading.ac.uk/clouds/adept/
Apache License 2.0
163 stars 29 forks source link

Computing multiple Jacobians from the same recording #13

Closed VictorCeballos closed 5 years ago

VictorCeballos commented 5 years ago

Hello, I am having problems when computing different Jacobians (using different sets of independent variables) from the same recording.

Consider the following code to be differentiated:

    adept::Stack stack;
    adept::aVector3 u = {1.0, 2.0, 3.0};
    adept::aReal h = 5;
    stack.new_recording();
    adept::aVector3 v = adept::Vector3{2,3,4} * u + 9 * h;

then, if we compute the derivative of v wrt to h, the following code is working just fine:

    int n = 1; // independent
    int m = 3; // dependent

    adept::Real jac[3] = {0.0};
    //stack.independent(u);
    stack.independent(h);
    stack.dependent(v);
    stack.jacobian(jac);
    print_jacobian(m, n, jac); // jac = {9.0, 9.0, 9.0}

However, if we first compute the full Jacobian, then clear the independent variables and add only h, the result is not right:

    int n = 3+1; // independent
    int m = 3; // dependent

    adept::Real jac[12] = {0.0};
    stack.independent(u);
    stack.independent(h);
    stack.dependent(v);
    stack.jacobian(jac);
    print_jacobian(m, n, jac); // jac = { 2.0, 0.0, 0.0, 9.0, 
                               //         0.0, 3.0, 0.0, 9.0, 
                               //         0.0, 0.0, 4.0, 9.0 } RIGHT

    stack.clear_independents();
    stack.clear_dependents();

    n = 1;
    m = 3;

    adept::Real jac1[3] = {0.0};
    //stack.independent(u);
    stack.independent(h);
    stack.dependent(v);
    stack.jacobian(jac1);
    print_jacobian(m, n, jac1); // jac1 = {9.0, 9.0, 25.0} WRONG!

Am I doing something wrong? Is this a problem with the clear_independents/clear_dependents functionality? This happens to me on Ubuntu 18.04.2 with GCC 7.4.0

Thanks!

-- For reference

void print_jacobian(int m, int n, adept::Real* jac)
{
    for (int i = 0; i < m; ++i)
    {
        for (int j = 0; j < n; ++j)
        {
            printf("\t%.3f", jac[m*j+i]);
        }
        printf("\n");
    }
}
rjhogan commented 5 years ago

Hi, Strange - it works for me on OpenSUSE with GCC 6.3.0, 7.3.0 and 8.3.0. The full code I ran is:

include "adept_arrays.h"

include "iostream"

void print_jacobian(int m, int n, adept::Real jac) { for (int i = 0; i < m; ++i) { for (int j = 0; j < n; ++j) { printf("\t%.3f", jac[mj+i]); } printf("\n"); } }

int main() { adept::Stack stack; adept::aVector3 u = {1.0, 2.0, 3.0}; adept::aReal h = 5; stack.new_recording(); adept::aVector3 v = adept::Vector3{2,3,4} u + 9 h; int n = 3+1; // independent int m = 3; // dependent

adept::Real jac[12] = {0.0}; stack.independent(u); stack.independent(h); stack.dependent(v); stack.jacobian(jac); print_jacobian(m, n, jac); // jac = { 2.0, 0.0, 0.0, 9.0, // 0.0, 3.0, 0.0, 9.0, // 0.0, 0.0, 4.0, 9.0 } RIGHT stack.clear_independents(); stack.clear_dependents();

n = 1; m = 3;

adept::Real jac1[3] = {0.0}; //stack.independent(u); stack.independent(h); stack.dependent(v);

std::cout << stack; stack.print_statements();

stack.jacobian(jac1); print_jacobian(m, n, jac1); }

This produces the following output:

    2.000   0.000   0.000   9.000
    0.000   3.000   0.000   9.000
    0.000   0.000   4.000   9.000

Automatic Differentiation Stack (address 0x7ffc1101f690): Currently attached - thread safe Recording status: Recording is ON 3 statements (1048576 allocated) and 6 operations (1048576 allocated) 7 gradients currently registered and a total of 7 needed (current index 7) Gradient list has no gaps Computation status: 0 gradients assigned (0 allocated) Jacobian size: 3x1 Independent indices: 3 Dependent indices: 4 5 6 Parallel Jacobian calculation not available 1: d[4] = + 2d[0] + 9d[3] 2: d[5] = + 3d[1] + 9d[3] 3: d[6] = + 4d[2] + 9d[3] 9.000 9.000 9.000

You can see that the Jacobian is correct, and the list of differential statements correctly matches your operations. What output do you get when these stack print statements are included?

Robin.

VictorCeballos commented 5 years ago

Hi, I'm sorry I didn't do further testings before I opened the issue. It seems that when I compile this simple code on its own I do obtain the expected results:

    2.000   0.000   0.000   9.000
    0.000   3.000   0.000   9.000
    0.000   0.000   4.000   9.000
Automatic Differentiation Stack (address 0x7ffcf0ceb5e0):
   Currently attached - thread safe
   Recording status:
      Recording is ON
      3 statements (1048576 allocated) and 6 operations (1048576 allocated)
      4 gradients currently registered and a total of 4 needed (current index 4)
      Gradient list has no gaps
   Computation status:
      0 gradients assigned (0 allocated)
      Jacobian size: 3x1
      Independent indices: 0
      Dependent indices:   1 2 3
      Parallel Jacobian calculation not available
1: d[1] =  + 2*d[3] + 9*d[0]
2: d[2] =  + 3*d[4] + 9*d[0]
3: d[3] =  + 4*d[5] + 9*d[0]
    9.000
    9.000
    9.000

However, when I compile this code as part of a larger project, the jacobians start getting wrong values:

    2.000   0.000   0.000   9.000
    0.000   3.000   0.000   9.000
    0.000   0.000   4.000   9.000
Automatic Differentiation Stack (address 0x7ffd4d887c90):
   Currently attached - thread safe
   Recording status:
      Recording is ON
      3 statements (1048576 allocated) and 6 operations (1048576 allocated)
      4 gradients currently registered and a total of 4 needed (current index 4)
      Gradient list has no gaps
   Computation status:
      0 gradients assigned (0 allocated)
      Jacobian size: 3x1
      Independent indices: 0
      Dependent indices:   1 2 3
      Parallel Jacobian calculation not available
1: d[1] =  + 2*d[3] + 9*d[0]
2: d[2] =  + 3*d[4] + 9*d[0]
3: d[3] =  + 4*d[5] + 9*d[0]
    9.000
    9.000
    25.000

I do not understand what I am doing wrong, as even if I'm compiling alongside other files in the project, this is the only code that gets executed (in a main.cpp).

I have now tried to run this code from within the application, which is a Qt application, and it crashes.

However, I have managed to create a minimal Qt C++ application that displays the same problem (the wrong Jacobians):

main.cpp


#include <QApplication>
#include "mainwindow.h"

int main(int argc, char argv[]) { QApplication a(argc, argv); MainWindow w = new MainWindow(); w->show(); return a.exec(); }


> mainwindow.cpp

ifndef MAINWINDOW_H

define MAINWINDOW_H

include

class MainWindow : public QMainWindow { Q_OBJECT

public:

MainWindow();

void test();

};

endif // MAINWINDOW_H


> mainwindow.cpp

include "mainwindow.h"

include

include

include "adept_arrays.h"

include "iostream"

MainWindow::MainWindow() { QPushButton* button = new QPushButton;

QHBoxLayout *layout = new QHBoxLayout;
layout->addWidget(button);

QWidget *window = new QWidget();
window->setLayout(layout);

setCentralWidget(window);

QObject::connect(button, &QPushButton::clicked, this, &MainWindow::test);

}

void print_jacobian(int m, int n, adept::Real jac) { for (int i = 0; i < m; ++i) { for (int j = 0; j < n; ++j) { printf("\t%.3f", jac[mj+i]); } printf("\n"); } }

void MainWindow::test() { adept::Stack stack; adept::aVector3 u = {1.0, 2.0, 3.0}; adept::aReal h = 5; stack.new_recording(); adept::aVector3 v = adept::Vector3{2,3,4} u + 9 h;

int n = 3+1; // independent
int m = 3;   // dependent

adept::Real jac[12] = {0.0};
stack.independent(u);
stack.independent(h);
stack.dependent(v);
stack.jacobian(jac);
print_jacobian(m, n, jac);

stack.clear_independents();
stack.clear_dependents();

n = 1;
m = 3;

adept::Real jac1[3] = {0.0};
//stack.independent(u);
stack.independent(h);
stack.dependent(v);

std::cout << stack;
stack.print_statements();

stack.jacobian(jac1);
print_jacobian(m, n, jac1);

}

Sometimes it gets even worse:
2.000   221.000 0.000   9.000
0.000   3.000   0.000   9.000
0.000   221.000 4.000   9.000

Automatic Differentiation Stack (address 0x7ffdfa775c30): Currently attached - thread safe Recording status: Recording is ON 3 statements (1048576 allocated) and 6 operations (1048576 allocated) 4 gradients currently registered and a total of 4 needed (current index 4) Gradient list has no gaps Computation status: 0 gradients assigned (0 allocated) Jacobian size: 3x1 Independent indices: 0 Dependent indices: 1 2 3 Parallel Jacobian calculation not available 1: d[1] = + 2d[3] + 9d[0] 2: d[2] = + 3d[4] + 9d[0] 3: d[3] = + 4d[5] + 9d[0] 9.000 672.000 25.000 2.000 221.000 4.000 9.000 0.000 6.000 0.000 9.000 0.000 221.000 8.000 9.000 Automatic Differentiation Stack (address 0x7ffdfa775c30): Currently attached - thread safe Recording status: Recording is ON 3 statements (1048576 allocated) and 6 operations (1048576 allocated) 4 gradients currently registered and a total of 4 needed (current index 4) Gradient list has no gaps Computation status: 0 gradients assigned (0 allocated) Jacobian size: 3x1 Independent indices: 0 Dependent indices: 1 2 3 Parallel Jacobian calculation not available 1: d[1] = + 2d[3] + 9d[0] 2: d[2] = + 3d[4] + 9d[0] 3: d[3] = + 4d[5] + 9d[0] 9.000 672.000 41.000 free(): invalid next size (fast)


I'm sorry to raise this issue, as it seems a bit tangential to Adept...

I've done some further testing and it seems that Forward-mode differentiation is behaving ok.
I might be able to add autodiff to my application in this way. However, I cannot seem to be able to differentiate wrt a vector. Is this something that Adept can do?

Specifically, this fails to compile:
`adept::set_gradients(u, 3, u_ad);`

Here's the full code:

include

include "adept_arrays.h"

int main() { adept::Stack stack; adept::aVector3 u = {1.0, 2.0, 3.0}; adept::aReal h = 5; stack.new_recording(); adept::aVector3 v = adept::Vector3{2,3,4} u + 9 h;

// FORWARD
{
    std::cout << "FORWARD" << std::endl;

    h.set_gradient(1.0);
    stack.forward();
    adept::Vector dv_dh = v.get_gradient();
    std::cout << dv_dh << std::endl;
    stack.clear_gradients();

    //adept::Vector u_ad = {1.0, 1.0, 1.0};
    //adept::set_gradients(u, 3, u_ad);

    // compile error:

    // error: no matching function for call to ‘set_gradients(adept::aVector3&, int, adept::Vector&)’
    //    adept::set_gradients(u, 3, u_ad);
    // candidate: template<class Type> void adept::set_gradients(adept::Active<Type>*, adept::Index, const Type*)
    //    void set_gradients(Active<Type>* a, Index n, const Type* data)
    // mismatched types ‘adept::Active<Type>*’ and ‘adept::FixedArray<double, true, 3>’

    //stack.forward();
    //adept::Vector dv_du = v.get_gradient();
    //std::cout << dv_du << std::endl;
    //stack.clear_gradients();

    // alternative: this works...

    u[0].set_gradient(1.0);
    stack.forward();
    adept::Vector dv_du0 = v.get_gradient();
    std::cout << dv_du0 << std::endl;
    stack.clear_gradients();

    // but...

    //u[1].set_gradient(1.0);

    // run-time error:

    // terminate called after throwing an instance of 'adept::gradient_out_of_range'
    //   what():  Gradient index out of range: probably aReal objects have been created after a set_gradient(s) call

    //stack.forward();
    //adept::Vector dv_du1 = v.get_gradient();
    //std::cout << dv_du1 << std::endl;
    //stack.clear_gradients();

    //u[2].set_gradient(1.0);
    //stack.forward();
    //adept::Vector dv_du2 = v.get_gradient();
    //std::cout << dv_du2 << std::endl;
    //stack.clear_gradients();
}

// REVERSE
{
    std::cout << "REVERSE" << std::endl;

    // does not compile

    //adept::Vector3 v_ad = {1.0, 1.0, 1.0};
    //adept::set_gradients(v, 3, v_ad);
    //stack.reverse();
    //adept::Vector3 dv_du;
    //adept::get_gradients(u, dv_du);
    //adept::Real dv_dh = h.get_gradient();

    //std::cout << dv_du << std::endl;
    //std::cout << dv_dh << std::endl;
}

}



Thanks!
rjhogan commented 5 years ago

I've had a go at compiling a Qt application but it looks like a bit of a hassle and I don't have time to figure it out. Can you make a failing application that doesn't depend on Qt?

VictorCeballos commented 5 years ago

Hi, While trying to come up with a simpler application I've realised that simply adding a print statement before the Adept code does the trick:

main.cpp


#include <iostream>
#include <adept_arrays.h>

void print_jacobian(int m, int n, adept::Real jac) { for (int i = 0; i < m; ++i) { for (int j = 0; j < n; ++j) { printf("\t%.3f", jac[mj+i]); } printf("\n"); } }

int main(int argc, char **argv) { std::cout << "Hello World!" << std::endl;

adept::Stack stack;
adept::aVector3 u = {1.0, 2.0, 3.0};
adept::aReal h = 5;
stack.new_recording();
adept::aVector3 v = adept::Vector3{2,3,4} * u + 9 * h;

int n = 3+1; // independent
int m = 3;   // dependent

adept::Real jac[12] = {0.0};
stack.independent(u);
stack.independent(h);
stack.dependent(v);
stack.jacobian(jac);
print_jacobian(m, n, jac);

stack.clear_independents();
stack.clear_dependents();

n = 1;
m = 3;

adept::Real jac1[3] = {0.0};
//stack.independent(u);
stack.independent(h);
stack.dependent(v);

std::cout << stack;
stack.print_statements();

stack.jacobian(jac1);
print_jacobian(m, n, jac1);

return 0;

}


Can you confirm that this is failing for you as well?

Also, I've been working on this further and I've seen that if I use native arrays or std::vectors of adoubles (instead of FixedArrays) everything works as expected. So it seems this is a problem with Adept arrays and not the computation of the Jacobians.
rjhogan commented 5 years ago

Hi,

At work we have the following versions of g++:

gnu/4.8.5 gnu/5.3.0(old) gnu/6.3.0(default) gnu/8.3.0 gnu/4.9.1 gnu/6.1.0 gnu/7.3.0(jun19:new)

I can get it to fail with 4.8.5 only - the others are fine. But this is enough to convince me that there is probably a bug in Adept so I will investigate. The failing run outputs this line:

  Independent indices: 0

which implies that the independent variable has been indexed wrongly somehow.

I'll let you know if I figure it out - will probably take over a week to get on to it properly.

Robin

On 15/05/19 12:59, Victor Ceballos wrote:

Hi, While trying to come up with a simpler application I've realised that simply adding a print statement before the Adept code does the trick:

main.cpp

include

include

void print_jacobian(int m, int n, adept::Real jac) { for (int i = 0; i < m; ++i) { for (int j = 0; j < n; ++j) { printf("\t%.3f", jac[mj+i]); } printf("\n"); } }

int main(int argc, char **argv) { std::cout << "Hello World!" << std::endl;

adept::Stack stack;
adept::aVector3 u = {1.0, 2.0, 3.0};
adept::aReal h = 5;
stack.new_recording();
adept::aVector3 v = adept::Vector3{2,3,4} * u + 9 * h;

int n = 3+1; // independent
int m = 3;   // dependent

adept::Real jac[12] = {0.0};
stack.independent(u);
stack.independent(h);
stack.dependent(v);
stack.jacobian(jac);
print_jacobian(m, n, jac);

stack.clear_independents();
stack.clear_dependents();

n = 1;
m = 3;

adept::Real jac1[3] = {0.0};
//stack.independent(u);
stack.independent(h);
stack.dependent(v);

std::cout << stack;
stack.print_statements();

stack.jacobian(jac1);
print_jacobian(m, n, jac1);

return 0;

}

Can you confirm that this is failing for you as well?

Also, I've been working on this further and I've seen that if I use native arrays or std::vectors of adoubles (instead of FixedArrays) everything works as expected. So it seems this is a problem with Adept arrays and not the computation of the Jacobians.

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/rjhogan/Adept-2/issues/13?email_source=notifications&email_token=ADC2EWJVS3C3UU2FWURQ2ULPVP3JXA5CNFSM4HMQBFG2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODVONVPY#issuecomment-492624575, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ADC2EWMFSXNC7H6LQRS2PBDPVP3JXANCNFSM4HMQBFGQ.

VictorCeballos commented 5 years ago

Hi, No problem - I'm glad we've managed to reproduce the problem.

I've noticed that the indices are not consistent across different applications. For the failing application, I get:

      Independent indices: 0
      Dependent indices:   1 2 3
      Parallel Jacobian calculation not available
1: d[1] =  + 2*d[3] + 9*d[0]
2: d[2] =  + 3*d[4] + 9*d[0]
3: d[3] =  + 4*d[5] + 9*d[0]

Shouldn't this be instead?:

1: d[1] =  + 2*d[4] + 9*d[0]
2: d[2] =  + 3*d[5] + 9*d[0]
3: d[3] =  + 4*d[6] + 9*d[0]

Or something along these lines? It seems some indices get used twice.

Strangely, these expressions are the same even when I comment out the printing line (hence obtaining correct Jacobians). They are the expressions I've been getting all along.

rjhogan commented 5 years ago

Hi, Are you using the 2.0.5 tarball, or the latest git snapshot? The tarball is missing this fix: https://github.com/rjhogan/Adept-2/commit/be51a8e28c847872d58b666c84ab7810803bbc99 that solves a problem with initializing FixedArrays from initializer_lists, which could have led to the indices not being registered properly. I also found that problems occur when I compiled the adept library with one version of the compiler and then tried to compile/link my application using another version. Less likely as a problem, but did you define ADEPT_SUPPORT_HUGE_ARRAYS? (safer not to) Do any of these help? Robin.

VictorCeballos commented 5 years ago

Hi, I was using the 2.0.5 tarball from the website. Using the latest version from git solves the problem. Thanks for the help!