clad icon indicating copy to clipboard operation
clad copied to clipboard

error: no member named 'estimate_error' in namespace 'clad'

Open frhack opened this issue 1 year ago • 19 comments

Hi I'm successfully using clad, but cannot use estimate_error:

Can someone help me?

C++17clad C++14clad C++11clad standard clad tutorial

#include "clad/Differentiator/Differentiator.h"

float function_fn1(float x, float y) {
  float z;
  z = x + y;
  return z;
}
#include "clad/Differentiator/Differentiator.h"
auto df = clad::estimate_error(function_fn1);
//auto df = clad::estimate_error(function_fn1);
// Print the generated code to standard output.
df.dump();
// Declare the necessary variables.
double x, y, d_x, d_y, final_error = 0;
// Finally call execute on the generated code.
df.execute(x, y, &d_x, &d_y, final_error);

final_error; input_line_20:2:18: error: no member named 'estimate_error' in namespace 'clad'

frhack avatar Jul 15 '22 16:07 frhack

Hi! Are you perhaps using clad in the cling environment? If so, the clad version built by default with cling is the latest release which does not have the error estimation abilities. If you would like, you can freely use it outside of cling or build the latest master of clad with cling. Hope this helps!

grimmmyshini avatar Jul 15 '22 18:07 grimmmyshini

Hi! Are you perhaps using clad in the cling environment? If so, the clad version built by default with cling is the latest release which does not have the error estimation abilities. If you would like, you can freely use it outside of cling or build the latest master of clad with cling. Hope this helps!

Hi thanks, yes I'm using in cling. Is there a way to use in cling ? I'm in hurry for a presentation and all my environment is in cling

More precisely I had setup a working mybinder environment with:

  • numpy
  • clad
  • xeus-cling
  • matplotlib
  • cmake
  • python=3.7

Everything is working except to estimate_error

frhack avatar Jul 15 '22 18:07 frhack

Hi @frhack, when is your presentation deadline. I am available to meet tomorrow afternoon Geneva time. Can you send me an email and we can schedule a zoom call to unblock you.

vgvassilev avatar Jul 15 '22 19:07 vgvassilev

Hi @frhack, when is your presentation deadline. I am available to meet tomorrow afternoon Geneva time. Can you send me an email and we can schedule a zoom call to unblock you.

Hi I have and exam monday and I prepared a notebook, it's almost finisced, I'd like to add a a floating point error estimation

Thank you very much I'll send you an email with the link to the notebook

frhack avatar Jul 15 '22 19:07 frhack

hi @vgvassilev I sent you an email to the cern address

frhack avatar Jul 16 '22 14:07 frhack

@frhack, I have not received it. Can you try joining [snip] I will be there...

vgvassilev avatar Jul 16 '22 15:07 vgvassilev

This is the code I would try:



double derivative_fn(double x) {
  return (exp(x)*(cos(3*x) + sin(3*x)/2 + (3*sin(x))/2)) / 
            pow(pow(cos(x),3) + pow(sin(x),3),2);
}

auto df = clad::estimate_error(derivative_fn);

// Print the generated code to standard output.
df.dump();

double pi =  2 * asin(1);
double d_x, final_error = 0;
double point_x =  pi/4;


// compute the error in point_x
df.execute(point_x, &d_x, final_error);

....finally output the error

frhack avatar Jul 16 '22 16:07 frhack

hi @vgvassilev is everything clear ?

If nested calls are a problem, we can use temporary variables to compute intermediary values

frhack avatar Jul 17 '22 07:07 frhack

Hi @frhack,

here is the output I get: The file with slight modifications (eg std:: prefix)

cat TT.cpp 
#include "clad/Differentiator/Differentiator.h"
double derivative_fn(double x) {
  return (std::exp(x)*(std::cos(3*x) + std::sin(3*x)/2 + (3*std::sin(x))/2)) /
    std::pow(std::pow(std::cos(x),3) + std::pow(std::sin(x),3),2);
}

extern "C" int printf(const char*,...);
int main () {
  auto df = clad::estimate_error(derivative_fn);

  // Print the generated code to standard output.
  df.dump();

  double pi =  2 * asin(1);
  double d_x, final_error = 0;
  double point_x =  pi/4;


  // compute the error in point_x
  df.execute(point_x, &d_x, final_error);
  printf("final_error='%f'\n", final_error);
}

The compiler invocation with clad master:

vvassilev@vv-nuc ~/workspace/builds/clad $ clang-12 -std=c++11 -Xclang -add-plugin -Xclang clad -Xclang -plugin-arg-clad -Xclang -fdump-derived-fn -Xclang -load -Xclang lib/clad.so -I ../../sources/clad/include TT.cpp -lm -lstdc++
void derivative_fn_grad(double x, clad::array_ref<double> _d_x, double &_final_error) {
    double _t0;
    double _t1;
    double _t2;
    double _ret_value0 = 0;
    _t1 = 3 * x;
    _t2 = 3 * x;
    _t0 = (std::cos(_t1) + std::sin(_t2) / 2 + (3 * std::sin(x)) / 2);
    double derivative_fn_return = (std::exp(x) * _t0) / std::pow(std::pow(std::cos(x), 3) + std::pow(std::sin(x), 3), 2);
    _ret_value0 = derivative_fn_return;
    goto _label0;
  _label0:
    {
        typename __gnu_cxx::__promote_2<double, int>::__type _r0 = std::pow(std::pow(std::cos(x), 3) + std::pow(std::sin(x), 3), 2);
        double _r1 = 1 / _r0;
        double _r2 = _r1 * _t0;
        double _r3 = _r2 * clad::custom_derivatives::exp_pushforward(x, 1.).pushforward;
        * _d_x += _r3;
        double _r4 = std::exp(x) * _r1;
        double _r5 = _r4 * clad::custom_derivatives::cos_pushforward(_t1, 1.).pushforward;
        double _r6 = _r5 * x;
        double _r7 = _r4 / 2;
        double _r8 = _r7 * clad::custom_derivatives::sin_pushforward(_t2, 1.).pushforward;
        double _r9 = _r8 * x;
        double _r10 = _r4 / 2;
        double _r11 = _r10 * std::sin(x);
    }
    double _delta_x = 0;
    _delta_x += std::abs(* _d_x * x * 1.1920928955078125E-7);
    _final_error += _delta_x + std::abs(1. * _ret_value0 * 1.1920928955078125E-7);
}

The output:

vvassilev@vv-nuc ~/workspace/builds/clad $ ./a.out 
The code is: 
void derivative_fn_grad(double x, clad::array_ref<double> _d_x, double &_final_error) {
    double _t0;
    double _t1;
    double _t2;
    double _ret_value0 = 0;
    _t1 = 3 * x;
    _t2 = 3 * x;
    _t0 = (std::cos(_t1) + std::sin(_t2) / 2 + (3 * std::sin(x)) / 2);
    double derivative_fn_return = (std::exp(x) * _t0) / std::pow(std::pow(std::cos(x), 3) + std::pow(std::sin(x), 3), 2);
    _ret_value0 = derivative_fn_return;
    goto _label0;
  _label0:
    {
        typename __gnu_cxx::__promote_2<double, int>::__type _r0 = std::pow(std::pow(std::cos(x), 3) + std::pow(std::sin(x), 3), 2);
        double _r1 = 1 / _r0;
        double _r2 = _r1 * _t0;
        double _r3 = _r2 * clad::custom_derivatives::exp_pushforward(x, 1.).pushforward;
        * _d_x += _r3;
        double _r4 = std::exp(x) * _r1;
        double _r5 = _r4 * clad::custom_derivatives::cos_pushforward(_t1, 1.).pushforward;
        double _r6 = _r5 * x;
        double _r7 = _r4 / 2;
        double _r8 = _r7 * clad::custom_derivatives::sin_pushforward(_t2, 1.).pushforward;
        double _r9 = _r8 * x;
        double _r10 = _r4 / 2;
        double _r11 = _r10 * std::sin(x);
    }
    double _delta_x = 0;
    _delta_x += std::abs(* _d_x * x * 1.1920928955078125E-7);
    _final_error += _delta_x + std::abs(1. * _ret_value0 * 1.1920928955078125E-7);
}

final_error='0.000001'

Let me know if I should run something else.

vgvassilev avatar Jul 17 '22 08:07 vgvassilev

@vgvassilev wonderful thanks, just one more thing... can you please repeat the calculations replacing every occurrence of "double" with "long double" in order to see how the error change?

frhack avatar Jul 17 '22 09:07 frhack

@frhack if you want to have a more detailed overview (this will increase the error from the default model, then you can split your function like so (and template it so that comparison becomes easier.)

template <typename T = double>
T derivative_fn(T x) {
  T c1 = std::exp(x);
  T m1 = std::cos(3 * x) + std::sin(3 * x) / 2 + (3 * std::sin(x)) / 2;
  T t = std::pow(std::pow(std::cos(x), 3) + std::pow(std::sin(x), 3), 2);
  T result = (c1 * m1) / t;
  return result;
}

and the main function will look like this then:

extern "C" int printf(const char*,...);

int main () {
  auto df = clad::estimate_error(derivative_fn<>);

  // Print the generated code to standard output.
  df.dump();

  double pi =  2 * asin(1);
  double d_x, final_error = 0;
  long double point_x =  pi/4;

  // compute the error in point_x
  df.execute(point_x, &d_x, final_error);
  printf("final_error='%.10f'\n", final_error);
  printf("Actual Error = %.10Lf\n",
         std::abs(derivative_fn<long double>(point_x) - derivative_fn<float>(point_x)));
}

Here you can compare your actual function output to the estimated error. If you have some ground truth value for the function, you can compare that too. Note: the estimated error for the default model will be higher as it provides an upper bound. Finally, executing the above will give you this:

The code is: 
void derivative_fn_grad(double x, clad::array_ref<double> _d_x, double &_final_error) {
    double _t0;
    double _d_c1 = 0;
    double _delta_c1 = 0;
    double _EERepl_c10;
    double _t1;
    double _t2;
    double _t3;
    double _t4;
    double _t5;
    double _t6;
    double _d_m1 = 0;
    double _delta_m1 = 0;
    double _EERepl_m10;
    double _t7;
    double _t8;
    double _t9;
    double _t10;
    double _t11;
    double _d_t = 0;
    double _delta_t = 0;
    double _EERepl_t0;
    double _t12;
    double _t13;
    double _t14;
    double _t15;
    double _d_result = 0;
    double _delta_result = 0;
    double _EERepl_result0;
    _t0 = x;
    double c1 = std::exp(_t0);
    _EERepl_c10 = c1;
    _t1 = x;
    _t2 = 3 * _t1;
    _t3 = x;
    _t4 = 3 * _t3;
    _t6 = x;
    _t5 = std::sin(_t6);
    double m1 = std::cos(_t2) + std::sin(_t4) / 2 + (3 * _t5) / 2;
    _EERepl_m10 = m1;
    _t7 = x;
    _t8 = std::cos(_t7);
    _t9 = x;
    _t10 = std::sin(_t9);
    _t11 = std::pow(_t8, 3) + std::pow(_t10, 3);
    double t = std::pow(_t11, 2);
    _EERepl_t0 = t;
    _t14 = c1;
    _t13 = m1;
    _t15 = (_t14 * _t13);
    _t12 = t;
    double result = _t15 / _t12;
    _EERepl_result0 = result;
    double derivative_fn_return = result;
    goto _label0;
  _label0:
    _d_result += 1;
    {
        double _r20 = _d_result / _t12;
        double _r21 = _r20 * _t13;
        _d_c1 += _r21;
        double _r22 = _t14 * _r20;
        _d_m1 += _r22;
        double _r23 = _d_result * -_t15 / (_t12 * _t12);
        _d_t += _r23;
        _delta_result += std::abs(_d_result * _EERepl_result0 * 1.1920928955078125E-7);
    }
    {
        double _grad4 = 0.;
        int _grad5 = 0;
        clad::custom_derivatives::std::pow_pullback(_t11, 2, _d_t, &_grad4, &_grad5);
        double _r12 = _grad4;
        double _grad0 = 0.;
        int _grad1 = 0;
        clad::custom_derivatives::std::pow_pullback(_t8, 3, _r12, &_grad0, &_grad1);
        double _r13 = _grad0;
        double _r14 = _r13 * clad::custom_derivatives::cos_pushforward(_t7, 1.).pushforward;
        * _d_x += _r14;
        int _r15 = _grad1;
        double _grad2 = 0.;
        int _grad3 = 0;
        clad::custom_derivatives::std::pow_pullback(_t10, 3, _r12, &_grad2, &_grad3);
        double _r16 = _grad2;
        double _r17 = _r16 * clad::custom_derivatives::sin_pushforward(_t9, 1.).pushforward;
        * _d_x += _r17;
        int _r18 = _grad3;
        int _r19 = _grad5;
        _delta_t += std::abs(_d_t * _EERepl_t0 * 1.1920928955078125E-7);
    }
    {
        double _r1 = _d_m1 * clad::custom_derivatives::cos_pushforward(_t2, 1.).pushforward;
        double _r2 = _r1 * _t1;
        double _r3 = 3 * _r1;
        * _d_x += _r3;
        double _r4 = _d_m1 / 2;
        double _r5 = _r4 * clad::custom_derivatives::sin_pushforward(_t4, 1.).pushforward;
        double _r6 = _r5 * _t3;
        double _r7 = 3 * _r5;
        * _d_x += _r7;
        double _r8 = _d_m1 / 2;
        double _r9 = _r8 * _t5;
        double _r10 = 3 * _r8;
        double _r11 = _r10 * clad::custom_derivatives::sin_pushforward(_t6, 1.).pushforward;
        * _d_x += _r11;
        _delta_m1 += std::abs(_d_m1 * _EERepl_m10 * 1.1920928955078125E-7);
    }
    {
        double _r0 = _d_c1 * clad::custom_derivatives::exp_pushforward(_t0, 1.).pushforward;
        * _d_x += _r0;
        _delta_c1 += std::abs(_d_c1 * _EERepl_c10 * 1.1920928955078125E-7);
    }
    double _delta_x = 0;
    _delta_x += std::abs(* _d_x * x * 1.1920928955078125E-7);
    _final_error += _delta_x + _delta_result + _delta_t + _delta_c1 + _delta_m1;
}

final_error='0.0000020599'
Actual Error = 0.0000000460

grimmmyshini avatar Jul 17 '22 10:07 grimmmyshini

very interesting @grimmmyshini thanks! I understand that this is an upper bound estimation

but why you use float to estimate the actual error ? Shouldn't we have to use double ?

frhack avatar Jul 17 '22 10:07 frhack

That's because by default we use the float machine Epsilon. This can obviously be changed if you decide to implement your own custom models with something else.

grimmmyshini avatar Jul 17 '22 10:07 grimmmyshini

@vgvassilev wonderful thanks, just one more thing... can you please repeat the calculations replacing every occurrence of "double" with "long double" in order to see how the error change?

cat TT.cpp 
#include "clad/Differentiator/Differentiator.h"

long double derivative_fn(long double x) {
  return (std::exp(x)*(std::cos(3*x) + std::sin(3*x)/2 + (3*std::sin(x))/2)) /
    std::pow(std::pow(std::cos(x),3) + std::pow(std::sin(x),3),2);
}

extern "C" int printf(const char*,...);
int main () {
  auto df = clad::estimate_error(derivative_fn);

  // Print the generated code to standard output.
  df.dump();

  long double pi =  2 * asin(1);
  long double d_x;
  double final_error = 0;
  long double point_x =  pi/4;


  // compute the error in point_x
  //derivative_fn_grad(point_x, &d_x, final_error);
  df.execute(point_x, &d_x, final_error);
  printf("final_error='%f'\n", final_error);
}
vvassilev@vv-nuc ~/workspace/builds/clad $ ./a.out 
The code is: 
void derivative_fn_grad(long double x, clad::array_ref<long double> _d_x, double &_final_error) {
    long double _t0;
    long double _t1;
    long double _t2;
    double _ret_value0 = 0;
    _t1 = 3 * x;
    _t2 = 3 * x;
    _t0 = (std::cos(_t1) + std::sin(_t2) / 2 + (3 * std::sin(x)) / 2);
    long double derivative_fn_return = (std::exp(x) * _t0) / std::pow(std::pow(std::cos(x), 3) + std::pow(std::sin(x), 3), 2);
    _ret_value0 = derivative_fn_return;
    goto _label0;
  _label0:
    {
        typename __gnu_cxx::__promote_2<long double, int>::__type _r0 = std::pow(std::pow(std::cos(x), 3) + std::pow(std::sin(x), 3), 2);
        long double _r1 = 1 / _r0;
        long double _r2 = _r1 * _t0;
        long double _r3 = _r2 * clad::custom_derivatives::std::exp_pushforward(x, 1.L).pushforward;
        * _d_x += _r3;
        long double _r4 = std::exp(x) * _r1;
        long double _r5 = _r4 * clad::custom_derivatives::std::cos_pushforward(_t1, 1.L).pushforward;
        long double _r6 = _r5 * x;
        long double _r7 = _r4 / 2;
        long double _r8 = _r7 * clad::custom_derivatives::std::sin_pushforward(_t2, 1.L).pushforward;
        long double _r9 = _r8 * x;
        long double _r10 = _r4 / 2;
        long double _r11 = _r10 * std::sin(x);
    }
    double _delta_x = 0;
    _delta_x += std::abs(* _d_x * x * 1.1920928955078125E-7);
    _final_error += _delta_x + std::abs(1. * _ret_value0 * 1.1920928955078125E-7);
}

final_error='nan'

It seems we get a nan here because we statically synthesize double &_final_error as part of the error estimate?

@grimmmyshini, does it make sense?

vgvassilev avatar Jul 17 '22 10:07 vgvassilev

Yeah, it doesn't make a lot of sense to run the error estimation in long double mode since the Epsilon does not change.

grimmmyshini avatar Jul 17 '22 13:07 grimmmyshini

@vgvassilev with reference to the presentation on youtube ...

EL is not the error due to Taylor truncation ? You say "linearization", is it something different ?

frhack avatar Jul 17 '22 16:07 frhack

Hey! The machine Epsilon (E subscript M) is not the linearization error (E subscript L). You can read more about the machine epsilon here. And the linearization error is due to the fact the Taylor series is approximated to the first order.

grimmmyshini avatar Jul 17 '22 17:07 grimmmyshini

Hey! The machine Epsilon (E subscript M) is not the linearization error (E subscript L). You can read more about the machine epsilon here. And the linearization error is due to the fact the Taylor series is approximated to the first order.

thanks, yes machine epsilon is clear, I was asking about E subscript L If it's the error due to Taylor serie truncation

frhack avatar Jul 17 '22 18:07 frhack

Oh yes, that's the linearization error!

grimmmyshini avatar Jul 17 '22 18:07 grimmmyshini

We have released clad v1.0 and the binder should be up to date. I will mark this issue as resolved. Please re-open if necessary.

vgvassilev avatar Nov 18 '22 07:11 vgvassilev