clad icon indicating copy to clipboard operation
clad copied to clipboard

Clarifications about peruse of local variables.

Open nestor013 opened this issue 6 months ago • 2 comments

Hi, First thanks for this great tool !

Now when trying jacobian I found that I cannot use intermediate variables inside the "function". See the example below, it is a bit contrived here, but in more useful cases I could need some intermediate variables. As seen in the dump of h_1, _d_a3 is set to 0. Note that, it seems to work when using differentiate, and gradient.

Is this a limitation, an issue, or am I missing something ?

Here is the code:

#include "clad/Differentiator/Differentiator.h"

void h(double a, double b, double output[]) {
  output[0] = a * a * a;
  output[1] = a * a * a + b * b * b;
  output[2] = 2 * (a + b);
}

void h_1(double a, double b, double output[]) {
  double a3 = a * a * a;

  output[0] = a3;
  output[1] = a3 + b * b * b;
  output[2] = 2 * (a + b);
}

int main() {
  auto h_jac = clad::jacobian(h);
  auto h1_jac = clad::jacobian(h_1);
    
  h_jac.dump();
  h1_jac.dump();
}

And its output:

The code is: 
void h_jac(double a, double b, double output[], double *jacobianMatrix) {
    output[0] = a * a * a;
    output[1] = a * a * a + b * b * b;
    output[2] = 2 * (a + b);
    {
        jacobianMatrix[4UL] += 2 * 1;
        jacobianMatrix[5UL] += 2 * 1;
    }
    {
        jacobianMatrix[2UL] += 1 * a * a;
        jacobianMatrix[2UL] += a * 1 * a;
        jacobianMatrix[2UL] += a * a * 1;
        jacobianMatrix[3UL] += 1 * b * b;
        jacobianMatrix[3UL] += b * 1 * b;
        jacobianMatrix[3UL] += b * b * 1;
    }
    {
        jacobianMatrix[0UL] += 1 * a * a;
        jacobianMatrix[0UL] += a * 1 * a;
        jacobianMatrix[0UL] += a * a * 1;
    }
}

The code is: 
void h_1_jac(double a, double b, double output[], double *jacobianMatrix) {
    double _d_a3 = 0;
    double a3 = a * a * a;
    output[0] = a3;
    output[1] = a3 + b * b * b;
    output[2] = 2 * (a + b);
    {
        jacobianMatrix[4UL] += 2 * 1;
        jacobianMatrix[5UL] += 2 * 1;
    }
    {
        jacobianMatrix[3UL] += 1 * b * b;
        jacobianMatrix[3UL] += b * 1 * b;
        jacobianMatrix[3UL] += b * b * 1;
    }
}

nestor013 avatar Aug 24 '24 11:08 nestor013