vgvassilev / clad

clad -- automatic differentiation for C/C++
GNU Lesser General Public License v3.0
277 stars 118 forks source link

incorrect values for jacobian (#2) #473

Open FROL256 opened 2 years ago

FROL256 commented 2 years ago

Input:

struct CamInfo
{
  float projM[16];
  float width;
  float height;
};

void VertexShader(const CamInfo& u, float vx, float vy, float vz, 
                               float output[2])
{
  const float W    =   vx * u.projM[3] + vy * u.projM[7] + vz * u.projM[11] + u.projM[15]; 
  const float xNDC =  (vx * u.projM[0] + vy * u.projM[4] + vz * u.projM[ 8] + u.projM[12])/W;
  const float yNDC = -(vx * u.projM[1] + vy * u.projM[5] + vz * u.projM[ 9] + u.projM[13])/W;
  output[0] = (xNDC*0.5f + 0.5f)*u.width;
  output[1] = (yNDC*0.5f + 0.5f)*u.height; 
}

...
auto h_jac = clad::jacobian(VertexShader, "vx,vy,vz");

The output jacobian is always zero (well, the generated code is obviously wrong, some-thing is missed). If i perform some manual inlining:

static inline void VertexShader(const CamInfo& u, float vx, float vy, float vz, 
                                float output[2])
{
  const float W    =   vx * u.projM[3] + vy * u.projM[7] + vz * u.projM[11] + u.projM[15]; 
  output[0] = (((vx * u.projM[0] + vy * u.projM[4] + vz * u.projM[ 8] + u.projM[12])/W)*0.5f + 0.5f)*u.width;
  output[1] = ((-(vx * u.projM[1] + vy * u.projM[5] + vz * u.projM[ 9] + u.projM[13])/W)*0.5f + 0.5f)*u.height; 
}
...
auto h_jac = clad::jacobian(VertexShader, "vx,vy,vz");

It works better, but unfortunately jacobian is still numerically wrong (

If i split VertexShader in two functions VS_X and VS_Y, clad produces correct gradients for them. I have made a branch for this bug here: https://github.com/FROL256/clad_autodiff_examples/tree/jacobian_test1

The first col is printed from gradients, the second one is from jacobian:

345.373 393.846
345.373 0

69.6805 118.154
69.6805 0

30.2959 78.7692
30.2959 0
grimmmyshini commented 2 years ago

Hi @FROL256, thanks for reporting this issue! I can also reproduce this issue, the outputs are 0 because the output is never assigned to the jacobian. Could you also please send the arguments you executed the function with?

grimmmyshini commented 2 years ago

Some diagnostic info and my thoughts on this:

The minimal example to reproduce is:

void minimal(double x, double y, double out[2]) {
  double temp1 = x + y;
  double temp2 = (x * y);

  out[0] = temp1;
  out[1] = temp2;
}

this generates code where the output jacobian is not even assigned anything:

void minimal_jac(double x, double y, double out[2], double *jacobianMatrix) {
    double _d_temp1 = 0;
    double _t0;
    double _t1;
    double _d_temp2 = 0;
    double temp1 = x + y;
    _t1 = x;
    _t0 = y;
    double temp2 = (_t1 * _t0);
    out[0] = temp1;
    out[1] = temp2;
    {
        double _r0 = _d_temp2 * _t0;
        double _r1 = _t1 * _d_temp2;
    }
}

However this works perfectly well:

void minimal(double x, double y, double out[2]) {
  out[0] = x + y;
  out[1] = x * y;
}

results in

void minimal_jac(double x, double y, double out[2], double *jacobianMatrix) {
    double _t0;
    double _t1;
    out[0] = x + y;
    _t1 = x;
    _t0 = y;
    out[1] = x * y;
    {
        double _r0 = 1 * _t0;
        jacobianMatrix[2UL] += _r0;
        double _r1 = _t1 * 1;
        jacobianMatrix[3UL] += _r1;
    }
    {
        jacobianMatrix[0UL] += 1;
        jacobianMatrix[1UL] += 1;
    }
}

Even though they are virtually the same function. Also, no jacobian test encapsulates the first case (i.e. all tests are super trivial where expressions are directly being assigned to the output like in the second example here.). This makes me think if supporting that is even possible with the current setup....

FROL256 commented 2 years ago

Could you also please send the arguments you executed the function with?

Sure, everything is here.

https://github.com/FROL256/clad_autodiff_examples/blob/jacobian_test1/SourceFile.cpp

vgvassilev commented 1 year ago

Is that still an issue with the new version of clad?

FROL256 commented 1 year ago

didn't work with it for a while, srry (