Including Clad in cmpad

[bradbell]

I am trying to install clad on fedora38.

The current version of clang++ on fedora38 is 16.0.6, but the install instructions
https://github.com/vgvassilev/clad#how-to-install
say that version 16.0x is not supported. Does this mean I cannot use clad ?

[vgvassilev]

Hi @bradbell,

For me the readme currently says: At the moment, LLVM/Clang 5.0.x - 16.0.x are supported.. You should be able to use clad v1.2 with clang 16.0.6. If that does not work, you can install earlier version of clang and install earlier version of Clad.

Best, Vassil

[bradbell]

My mistake. I have gotten some simple clad examples to run on fedora38. My next question is how to store a clad gradient object in a class member variable ?

The gradient object is create (in the examples) using statements of the form:

   auto f_dx    = clad::gradient(f);

The problem is, if I want f_dx to be a member variable (so it does not need to be created each time), I need to declare the type of f_dx before the assignment (and I do not know it). Using typeid(f_dx).name() I get something like this:

typeid( f_dx ).name() = N4clad12CladFunctionIPFvPdNS_9array_refIvEEENS_8NoObjectELb1EEE

which is not helpful.

[vgvassilev]

Can you show the concrete example? We create the code of the function only once and it is quite cheap to copy the CladFunction around. If that is not sufficient, you can forward declare the gradient function and use it as a standalone function.

[bradbell]

I am trying to include clad in the following project:
https://cmpad.readthedocs.io/en/latest/

The only real work to do this is the create a class, derived from cmpad::gradient, that computes the gradient using clad; see
https://cmpad.readthedocs.io/en/latest/gradient.html

Here is an example of what that looks like for adolc:
https://cmpad.readthedocs.io/en/latest/adolc_gradient.hpp.html

During the setup function, I want the clad version of the gradient to store an object in the derived class member variable. That way I can separate the setup cost from the evaluation cost. The problem is I do not know to declare the variable that will hold the the clad gradient ?

[vgvassilev]

@vaithak, could you help @bradbell out with this?

[vaithak]

Sure, will look into this.

[bradbell]

@vaithak The project cmpad is under construction and hence is not yet up to my standards for documentation.

I think it would be easier to just answer my question above. In case that you would like to try to implement
include/cmpad/clad/gradient.hpp
see
https://github.com/bradbell/cmpad/tree/main/include/cmpad
https://cmpad.readthedocs.io/en/latest/get_package.sh.html

[vaithak]

Hi @bradbell, there are two options. The first option, as @vgvassilev suggested, is to use auto f_dx = and f_dx.execute(...) at the same place, the reason being that clad works at compile time to compute the function for gradient or derivative. So, the first statement won't incur any extra cost other than the CladFunction object creation, which is not expensive.

Suppose you have to store the generated function object in a member variable. In that case, the types are different for forward mode and reverse mode, and both depend on the type of the function to be differentiated.
For ex, if the function to be differentiated is

double f (double x, double y) {
  return x*y;
}

The corresponding type for the forward mode member variable (f_dx) will be:
clad::CladFunction <double (*) (double, double), clad::NoObect, false>
i.e., the first template param is the type of the function to be differentiated.

Similarly, for reverse mode, the type will be:
clad::CladFunction <void (*) (double, double, clad::array_ref<void>, clad::array_ref<void>), clad::NoObject, true>
i.e., the first template param is similar to the function to be differentiated but return type is void, along with added parameters types of clad::array_ref<void> for each input param type.

[bradbell]

Thanks for the quick response. I suspect this means that clads rates (speeds) will not be much different if one includes the setup time or not. For cppad_jit and cppadcg, the source code for the derivatives is computed at run time and then compiled and linked in,
so their seutp times are very large.

This brings up an issue. Does clad include all of the if-then branching in the function, so one does not have to re-compile for different argument values ?

[vaithak]

Yup, the generated function will be same irrespective of the value of the argument.

[bradbell]

My next issue is using a function object inside of the function that is being differentiated; Here is my example that fails / runs when I use (do not use) a function object:

# include "clad/Differentiator/Differentiator.h"
# include <iostream>

// USE_FUN_OBJ
// The progrogam crashes (runs) if USE_FUN_OBJ is 1 (0) 
# define USE_FUN_OBJ 1

//
// f_obj
struct fun_obj {
   double operator()(double *x, size_t n)
   {  double ret = 0;
      for (size_t i = 0; i < n; i++) {
         ret += (i+1) * x[i];
      }
      return ret;
   }
} f_obj;

//
// f
# if USE_FUN_OBJ
double f(double *x, size_t n)
{  return f_obj(x, n);
}
# else
double f(double *x, size_t n)
{  double ret = 0;
   for (size_t i = 0; i < n; i++) {
      ret += (i+1) * x[i];
   }
   return ret;
}
# endif

// main
int main() {
   //
   // n
   size_t n = 3;
   //
   // x
   double x[] = {1.0, 2.0, 3.0};
   //
   // fx
   double fx = f(x, n);
   std::cout << "f(x, n) = " << fx << "\n";
   //
   // fx
   fx = f_obj(x, n);
   std::cout << "f_obj(x, n) = " << fx << "\n";
   //
   // f_dx
   auto f_dx    = clad::gradient(f, "x");
   //
   // dx, dx_ref
   double dx[] = {0.0, 0.0, 0.0};
   clad::array_ref<double> dx_ref(dx, n);
   //
   // dx, dx_ref
   f_dx.execute(x, n, dx_ref);
   for(size_t i = 0; i < n; ++i)
      std::cout << "df/dx[" << i << "] = " << dx[i] << "\n";
   //
   return 0;
}

[vaithak]

Hi @bradbell, I tried running this locally to figure out how to fix this. It turns out, differentiation for the call expression to a function object inside a function is not supported in clad (cc: @vgvassilev for confirmation).

One way to fix this is by directly passing the function object in clad::gradient expression (clad::gradient(f_obj, "x")), along with this, another change is required to make it work, pass an empty object of the struct while calling execute method. So, the overall change in main will look like:

auto f_dx = clad::gradient(f_obj, "x");
...
...
fun_obj temp_obj;
f_dx.execute(x, n, &temp_obj, dx_ref);
...

[vgvassilev]

Hi @bradbell, I tried running this locally to figure out how to fix this. It turns out, differentiation for the call expression to a function object inside a function is not supported in clad (cc: @vgvassilev for confirmation).

One way to fix this is by directly passing the function object in clad::gradient expression (clad::gradient(f_obj, "x")), along with this, another change is required to make it work, pass an empty object of the struct while calling execute method. So, the overall change in main will look like:

auto f_dx = clad::gradient(f_obj, "x");
...
...
fun_obj temp_obj;
f_dx.execute(x, n, &temp_obj, dx_ref);
...

Hi @bradbell, I tried running this locally to figure out how to fix this. It turns out, differentiation for the call expression to a function object inside a function is not supported in clad (cc: @vgvassilev for confirmation).

That should work, can you open an issue with a minimal reproducing example?

[bradbell]

I have verified that your suggested changes work.
Thanks.

[bradbell]

I am still having trouble using function objects. Here is my next example that demonstrates the problem:

# include "clad/Differentiator/Differentiator.h"
# include <iostream>
# include <vector>

// CASE
// Case 1 and 2 run while case 3 crashes,
// The error message for case 3 makes it seem like case 2 should fail. 
// To be specific, it starts with:
//    no matching constructor for initialization of 'fun_1_class'
# define CASE 3

// FUN_CLASS, USE_CLAD
# if CASE == 1
# define FUN_CLASS fun_1_class
# define USE_CLAD  1
# endif
# if CASE == 2
# define FUN_CLASS fun_2_class
# define USE_CLAD 0
# endif
# if CASE == 3
# define FUN_CLASS fun_2_class
# define USE_CLAD 1
# endif


//
// fun_1_class
class fun_1_class {
private:
   // n_
   const size_t n_;
   //
public:
   // ctor
   fun_1_class(size_t n)
   : n_(n)
   { }
   //
   // operator
   double operator()(const double* x)
   {  double ret = 0;
      for (size_t i = 0; i < n_; i++) {
         ret += (i+1) * x[i];
      }
      return ret;
   }
};

//
// fun_2_class
class fun_2_class {
private:
   //
   // fun_1_
   fun_1_class fun_1_;
   //
public:
   // ctor
   fun_2_class(size_t n)
   : fun_1_(n)
   { }
   //
   // operator
   double operator()(const double* x)
   {  return fun_1_(x);
   }
};

//
// main
int main() {
   //
   // n
   size_t n = 4;
   //
   // x
   std::vector<double> x(n);
   for(size_t i = 0; i < n; ++i)
      x[i] = double(i + 1);
   //
   // fun
   FUN_CLASS fun(n);
   //
   // fx
   double fx = fun(x.data());
   std::cout << "fun(x.data()) = " << fx << "\n";
# if USE_CLAD
   //
   // f_dx
   auto f_dx    = clad::gradient(fun);
   //
   // dx, dx_ref
   std::vector<double> dx(n);
   for(size_t j = 0; j < n; ++j)
      dx[j] = 0.0;
   clad::array_ref<double> dx_ref(dx.data(), n);
   //
   // dx, dx_ref
   FUN_CLASS not_used(n);
   f_dx.execute(x.data(), &not_used, dx_ref);
   for(size_t i = 0; i < n; ++i)
      std::cout << "df/dx[" << i << "] = " << dx[i] << "\n";
# endif
   //
   return 0;
}

[vaithak]

This is a similar issue as earlier; the assertion error is different, which, as you noted, is related to the improper instantiation of a class object in the generated code.
The context is also slightly different; in the earlier case, a function was calling into a function object, but in this case, one function object calls the other. I also tried replicating your scenario (of cmpad), i.e., calling the algo function object from a struct object of clad, which is ultimately passed for differentiation but got the same error.

But the crux of the issue is the same, i.e., there is an issue in differentiating call expressions to a function object.

I created an issue for this: #625.

[bradbell]

I have subscribed to issue #625 and will hold off on adding clad to cmpad for now.

Thanks.

[bradbell]

So that you understand the context of my difficulty, I have pushed a clad branch with my current attempt at creating a gradient for cmpad. Here is the source code:
https://github.com/bradbell/cmpad/blob/clad/include/cmpad/clad/gradient.hpp

[bradbell]

I can confirm that the example in issue #625 now runs on my Fedora 38 system using the following commands:

clad>clad.sh issue_625.cpp
clang++ -Xclang -add-plugin -Xclang clad -Xclang -load -Xclang /home/bradbell/repo/cmpad.git/cpp/build/prefix/lib64/clad.so -I /home/bradbell/repo/cmpad.git/cpp/build/prefix/include -std=c++11 issue_625.cpp -o issue_625
clad.sh.sh: OK
clad>./issue_625
f(x, n) = 14
df/dx[0] = 1
df/dx[1] = 2
df/dx[2] = 3

Unfortunately, when I try to compile the code (CASE 3) above I get the following:

clad>mv discussion_618 discussion_618.cpp
clad>clad.sh discussion_618.cpp
clang++ -Xclang -add-plugin -Xclang clad -Xclang -load -Xclang /home/bradbell/repo/cmpad.git/cpp/build/prefix/lib64/clad.so -I /home/bradbell/repo/cmpad.git/cpp/build/prefix/include -std=c++11 discussion_618.cpp -o discussion_618
discussion_618.cpp:65:11: error: no matching constructor for initialization of 'fun_1_class'
   double operator()(const double* x)
          ^
...

[vaithak]

Hi @bradbell, this is a current assumption in clad. While differentiating a function containing a class object, it assumes that a default constructor for this class exists.
The reason for this is that at some places at the very start of the generated derivative function, it may declare placeholder objects for the user-defined class (placeholder because they are initialized using copy initialization later in the generated code).

For example, CASE 3 will have a statement like fun_2_class _t2;. This will only work if fun_2_class contains a default constructor. So, just adding a default constructor for fun_2_class will fix this.

For your use case, would it be possible to define a default constructor for your custom-defined class or give default values to the parameters of the already-defined constructor?

[bradbell]

I created the example below (which uses the default constructor for the function objects) and seems to have the same problem.

# include "clad/Differentiator/Differentiator.h"
# include <iostream>
# include <vector>

// CASE
// Case 1 and 2 run while case 3 crashes,
// The error message for case 3 makes is not much help. To be specific:
// ... Differentiator/Array.h:177:16: error: read-only variable is not assignable
//   m_arr[i] = arr[i];
//
# define CASE 3

// FUN_CLASS, USE_CLAD
# if CASE == 1
# define FUN_CLASS fun_1_class
# define USE_CLAD  1
# endif
# if CASE == 2
# define FUN_CLASS fun_2_class
# define USE_CLAD 0
# endif
# if CASE == 3
# define FUN_CLASS fun_2_class
# define USE_CLAD 1
# endif


//
// fun_1_class
class fun_1_class {
private:
   // n_
   size_t n_;
   //
public:
   // setup
   void setup(size_t n)
   {  n_ = n; }
   //
   // operator
   double operator()(const double* x)
   {  double ret = 0;
      for (size_t i = 0; i < n_; i++) {
         ret += (i+1) * x[i];
      }
      return ret;
   }
};

//
// fun_2_class
class fun_2_class {
private:
   //
   // fun_1_
   fun_1_class fun_1_;
   //
public:
   // setup
   void setup(size_t n)
   {  fun_1_.setup(n); }
   //
   // operator
   double operator()(const double* x)
   {  return fun_1_(x);
   }
};

//
// main
int main() {
   //
   // n
   size_t n = 4;
   //
   // x
   std::vector<double> x(n);
   for(size_t i = 0; i < n; ++i)
      x[i] = double(i + 1);
   //
   // fun
   FUN_CLASS fun;
   fun.setup(n);
   //
   // fx
   double fx = fun(x.data());
   std::cout << "fun(x.data()) = " << fx << "\n";
# if USE_CLAD
   //
   // f_dx
   auto f_dx    = clad::gradient(fun);
   //
   // dx, dx_ref
   std::vector<double> dx(n);
   for(size_t j = 0; j < n; ++j)
      dx[j] = 0.0;
   clad::array_ref<double> dx_ref(dx.data(), n);
   //
   // dx, dx_ref
   FUN_CLASS not_used;
   f_dx.execute(x.data(), &not_used, dx_ref);
   for(size_t i = 0; i < n; ++i)
      std::cout << "df/dx[" << i << "] = " << dx[i] << "\n";
# endif
   //
   return 0;
}

[vaithak]

This is a separate issue related to improper handling of const parameters. I have created a fix for this: #632, and verified that the above program works after this fix.
Also, thanks for taking the time and effort to report these issues. Really appreciated 👍🏼

[bradbell]

I am having trouble figuring out how to use cmake to build a simple project that uses clad. The test_clad.sh test script below tires to do this. When I run it, I get the message:

[ 50%] Building CXX object CMakeFiles/xam_1.dir/xam_1.cpp.o
/usr/bin/clang++  -I/home/bradbell/trash/work/clad.git/include -Xclang -add-plugin clad -load clad.git/build/lib64/clad.so -MD -MT CMakeFiles/xam_1.dir/xam_1.cpp.o -MF CMakeFiles/xam_1.dir/xam_1.cpp.o.d -o CMakeFiles/xam_1.dir/xam_1.cpp.o -c /home/bradbell/trash/work/xam_1.cpp
clang-16: error: no such file or directory: 'clad'
clang-16: error: no such file or directory: 'clad.git/build/lib64/clad.so'
clang-16: warning: -load: 'linker' input unused [-Wunused-command-line-argument]

test_clad.sh

#! /usr/bin/env bash
set -e -u
# -----------------------------------------------------------------------------
# bash function that echos and executes a command
echo_eval() {
	echo $*
	eval $*
}
# -----------------------------------------------------------------------------
#
# work
if [ ! -e work ]
then
   mkdir work
else
   read -p 'The work sub-directory already exists,  use it [y/n] ?' use
   if [ "$use"  != 'y' ]
   then
      echo 'aborting'
      exit 1
   fi
fi 
cd work
#
# clad.git
if [ ! -d clad.git ]
then
   git clone https://github.com/vgvassilev/clad.git clad.git
fi
cd clad.git
#
# work/clad.git/build
if [ ! -d build ]
then
   mkdir build
fi
cd build
cmake -S .. -B . -D LLVM_EXTERNAL_LIT="$(which lit)"
make -j $(nproc) 
#
# work
cd ../..
#
# clad_so
clad_so=$(find clad.git/build -name 'clad.so')
clad_so="$(pwd)/$clad_so"
#
# xam_1.cpp
cat << EOF > xam_1.cpp
#include "clad/Differentiator/Differentiator.h"
#include <iostream>

double f(double x, double y) { return x * y; }

int main() {
   // f_dx
   auto f_dx = clad::differentiate(f, "x");
   //
   // x, y, f_dx_value
   double x           = 3.0;
   double y           = 4.0;
   double f_dx_value = f_dx.execute(x, y);
   //
   // ok
   bool ok = f_dx_value == y;
   //
   if( ok )
      return 0;
   return 1;
}
EOF
#
# xam_1
echo_eval clang++ \
   -Xclang -add-plugin -Xclang clad \
   -Xclang -load -Xclang $clad_so \
   -I clad.git/include \
   -std=c++11 \
   xam_1.cpp -o xam_1
#
# xam_1
if ! ./xam_1
then
   echo 'test_clad.sh: Error'
   exit 1
fi
#
# CMakeLists.txt
cat << EOF > CMakeLists.txt
CMAKE_MINIMUM_REQUIRED(VERSION 3.5)
PROJECT(example_1)
INCLUDE_DIRECTORIES( \${PROJECT_SOURCE_DIR}/clad.git/include )
ADD_COMPILE_OPTIONS(
   -Xclang -add-plugin -Xclang clad
   -Xclang -load -Xclang $clad_so
)
ADD_EXECUTABLE(xam_1 xam_1.cpp)
EOF
#
# work/build
if [ ! -d build ]
then
   mkdir build
fi
cd build
cmake -S .. -B . \
   -D CMAKE_CXX_COMPILER=clang++ \
   -D CMAKE_VERBOSE_MAKEFILE=YES
make
#
# xam_1
if ! ./xam_1
then
   echo 'test_clad.sh: Error'
   exit 1
fi
#
echo 'test_clad.sh: OK'
exit 0

[vgvassilev]

Would that help you

clad/cmake/modules/AddCladBenchmark.cmake

Line 20 in bba8a47

function(CB_ADD_GBENCHMARK benchmark)

Here we attach the plugin to the invocation.

[bradbell]

So how would I modify my CMakeLists.txt file in the example above; i.e.,

CMAKE_MINIMUM_REQUIRED(VERSION 3.5)
PROJECT(example_1)
INCLUDE_DIRECTORIES( \${PROJECT_SOURCE_DIR}/clad.git/include )
ADD_COMPILE_OPTIONS(
   -Xclang -add-plugin -Xclang clad
   -Xclang -load -Xclang $clad_so
)
ADD_EXECUTABLE(xam_1 xam_1.cpp)

where $clad_so gets replaced by the location of the clad.so file ?

@vgvassilev
The following line was missing and has been added to test_clad.sh above:

clad_so="$(pwd)/$clad_so"

[vgvassilev]

Does it work in the end?

[bradbell]

Here is a modification to CMakeLists.txt that works for me:

CMAKE_MINIMUM_REQUIRED(VERSION 3.5)
PROJECT(example_1)
INCLUDE_DIRECTORIES( \${PROJECT_SOURCE_DIR}/clad.git/include )
ADD_EXECUTABLE(xam_1 xam_1.cpp)
TARGET_COMPILE_OPTIONS(xam_1 PUBLIC -fplugin=$clad_so )

[bradbell]

It it does not work. When I execute the script, the build with

echo_eval clang++ \
   -Xclang -add-plugin -Xclang clad \
   -Xclang -load -Xclang $clad_so \
   -I clad.git/include \
   -std=c++11 \
   xam_1.cpp -o xam_1

works but using cmake and make does not. The final output is:

[ 50%] Building CXX object CMakeFiles/xam_1.dir/xam_1.cpp.o
/usr/bin/clang++  -I/home/bradbell/trash/work/clad.git/include -Xclang -add-plugin clad -load /home/bradbell/trash/work/clad.git/build/lib64/clad.so -MD -MT CMakeFiles/xam_1.dir/xam_1.cpp.o -MF CMakeFiles/xam_1.dir/xam_1.cpp.o.d -o CMakeFiles/xam_1.dir/xam_1.cpp.o -c /home/bradbell/trash/work/xam_1.cpp
clang-16: error: no such file or directory: 'clad'
clang-16: warning: -load: 'linker' input unused [-Wunused-command-line-argument]
clang-16: warning: /home/bradbell/trash/work/clad.git/build/lib64/clad.so: 'linker' input unused [-Wunused-command-line-argument]

[vgvassilev]

Can you use:

echo_eval clang++ -fplugin=$clad_so
-I clad.git/include
-std=c++11
xam_1.cpp -o xam_1

[bradbell]

The solution to this problem is above; see
#618 (reply in thread)

[bradbell]

I have another case that is giving an unexpected compile error. To be specific:

 50%] Building CXX object CMakeFiles/example.dir/example.cpp.o
/usr/bin/clang++  -I/home/bradbell/trash/work/clad.git/include -g -fplugin=/home/bradbell/trash/work/clad.git/build/lib64/clad.so -MD -MT CMakeFiles/example.dir/example.cpp.o -MF CMakeFiles/example.dir/example.cpp.o.d -o CMakeFiles/example.dir/example.cpp.o -c /home/bradbell/trash/work/example.cpp
error: invalid operands to binary expression ('size_t' (aka 'unsigned long') and 'std::vector<double>')
In file included from /home/bradbell/trash/work/example.cpp:1:
...

Note that the first line in example.cpp is:

# include <vector>

If you change use_clad to have value '0' , instead of '1' (in the script below) the program runs fine. Below is the script that demonstrates this problem. I called this script clad_two.sh (because I have a clad_one.sh that works. )

#! /usr/bin/env bash
set -e -u
# -----------------------------------------------------------------------------
# use_clad
# This flag can have the value '0' or '1'
# If it is '0', the example.cpp program below works fine.
# If it is '1', the example.cpp program below crashes.
use_clad='1'
# -----------------------------------------------------------------------------
# bash function that echos and executes a command
echo_eval() {
	echo $*
	eval $*
}
# -----------------------------------------------------------------------------
#
# program
program='clad_two.sh'
#
# work
if [ ! -e work ]
then
   mkdir work
else
   r=''
   while [ "$r" != 'u' ] && [ "$r" != 'r' ] && [ "$r" != 'a' ] 
   do
      read -p 'The work directory exists:  use, remove, abort [u/r/a] ?' r
   done
   if [ "$r"  == 'a' ]
   then
      echo "$program: aborting"
      exit 1
   elif [ "$r" == 'r' ]
   then
      echo_eval rm -rf work
      mkdir work
   fi
fi 
cd work
#
# work/clad.git
if [ ! -d clad.git ]
then
   git clone https://github.com/vgvassilev/clad.git clad.git
fi
cd clad.git
#
# work/clad.git/build
if [ ! -d build ]
then
   mkdir build
fi
cd build
cmake -S .. -B . -D LLVM_EXTERNAL_LIT="$(which lit)"
make -j $(nproc) 
#
# work
cd ../..
#
# clad_so
clad_so=$(find clad.git/build -name 'clad.so' )
clad_so="$(pwd)/$clad_so"
#
# work/example.cpp
cat << EOF > example.cpp
# include <vector>
# include <clad/Differentiator/Differentiator.h>
//
# define USE_CLAD $use_clad
//
// det_of_minor
template <class Scalar>
Scalar det_of_minor(
   const std::vector<Scalar>&    a  ,
   size_t                        n  ,
   size_t                        m  ,
   std::vector<size_t>&          r  ,
   std::vector<size_t>&          c  )
{  assert( a.size() == n * n );
   assert( r.size() == n + 1 );
   assert( c.size() == n + 1 );
   // END PROTOTYPE
   //
   // R0 = R(0)
   size_t R0 = r[n];
   assert( R0 < n );
   //
   // Cj = C(0)
   size_t Cj = c[n];
   assert( Cj < n );
   //
   //
   // check if this is a 1 by 1 minor
   if( m == 1 ) return a[ R0 * n + Cj ];
   //
   // detM
   // initialize determinant of the minor M
   Scalar detM(0.0);
   //
   // sign
   // initialize sign of factor for next sub-minor
   int sign = 1;
   //
   // r
   // remove row with index 0 in M from all the sub-minors of M
   r[n] = r[R0];
   //
   // C(j-1)
   // initial index in c for previous column of the minor M
   size_t Cj1 = n;
   //
   // for each column of M
   for(size_t j = 0; j < m; j++)
   {
      // M[0,j] = A[ R0, Cj ]
      // element with index (0, j) in the minor M
      assert( Cj < n );
      Scalar M0j = a[ R0 * n + Cj ];
      //
      // remove column with index j in M to form next sub-minor S of M
      c[Cj1] = c[Cj];
      //
      // detS
      // compute determinant of S, the sub-minor of M with
      // row R(0) and column C(j) removed.
      Scalar detS = det_of_minor(a, n, m - 1, r, c);
      //
      // restore column with index j in represenation of M as a minor of A
      c[Cj1] = Cj;
      //
      // detM
      // include this sub-minor term in the summation
      if( sign > 0 )
         detM = detM + M0j * detS;
      else
         detM = detM - M0j * detS;
      //
      // advance to next column of M
      Cj1  = Cj;
      Cj   = c[Cj];
      sign = - sign;
   }
   //
   // r
   // restore row zero to the minor representation for M
   r[n] = R0;
   //
   // return the determinant of the minor M
   return detM;
}
//
// wrapper
struct wrapper {
   // n_
   size_t n_;
   //
   // ell_
   size_t ell_;
   //
   // r_, c_
   std::vector<size_t> r_, c_;
   //
   // setup
   void setup(size_t n)
   {  n_ = n; 
      ell_ = size_t( std::sqrt( double(n) ) );
      if( ell_ * ell_ != n )
         ++ell_;
      r_.resize(ell_ + 1);
      c_.resize(ell_ + 1);
      for(size_t i = 0; i < ell_; ++i)
      {  r_[i] = i + 1;
         c_[i] = i + 1;
      }
      r_[ell_] = 0;
      c_[ell_] = 0;
   }
   //
   // operator()
   double operator()(const double* x)
   {  std::vector<double> vx(n_);
      for(size_t j = 0; j < n_; ++j)
         vx[j] = x[j];
      return det_of_minor(vx, ell_, ell_, r_, c_);
   }
};
//
// main
int main() {
   //
   // ok
   bool ok = true;
   //
   // n_arg
   size_t n_arg = 4;
   //
   // x
   std::vector<double> x(n_arg);
   for(size_t i = 0; i < n_arg; ++i)
      x[i] = double(i + 1);
   //
   // fun
   wrapper fun;
   fun.setup(n_arg);
   //
   // fx
   double fx    = fun(x.data());
   double check = x[0] * x[3] - x[1] * x[2];
   ok &= fx == check;
   //
# if USE_CLAD   
   {  //
      // f_dx
      auto f_dx    = clad::gradient(fun);
      //
      // dx, dx_ref
      std::vector<double> dx(n_arg);
      for(size_t j = 0; j < n_arg; ++j)
         dx[j] = 0.0;
      clad::array_ref<double> dx_ref(dx.data(), n_arg);
      //
      // dx, dx_ref
      wrapper not_used;
      f_dx.execute(x.data(), &not_used, dx_ref);
      for(size_t i = 0; i < n_arg; ++i)
         ok &= dx[i] == double(i+ 1);
   }
# endif
   if( ok )
      return 0;
   return 1;
}
EOF
#
# work/CMakeLists.txt
cat << EOF > CMakeLists.txt
CMAKE_MINIMUM_REQUIRED(VERSION 3.5)
SET(CMAKE_VERBOSE_MAKEIFLE TRUE)
PROJECT(example_1)
SET(CMAKE_BUILD_TYPE DEBUG)
INCLUDE_DIRECTORIES( 
   \${PROJECT_SOURCE_DIR}/clad.git/include 
)
ADD_EXECUTABLE(example example.cpp)
TARGET_COMPILE_OPTIONS(example PUBLIC -fplugin=$clad_so )
EOF
#
# work/build
if [ ! -d build ]
then
   mkdir build
fi
cd build
if [ -e CMakeCache.txt ]
then
   rm CMakeCache.txt
fi
cmake -S .. -B . \
   -D CMAKE_CXX_COMPILER=clang++ \
   -D CMAKE_VERBOSE_MAKEFILE=YES
#
# work/build/example
make clean
make example
if ! ./example
then
   echo "$program: work/build/example: Error"
   exit 1
else
   echo ' -----------------------------------'
   echo "$program: work/build/example: OK"
   echo ' -----------------------------------'
fi
#
echo "$program: OK"
exit 0

[vaithak]

Currently, clad has very limited support for handling stl vectors, especially when the param to be differentiated is itself a vector or an index of it. I tried some basic examples with the vector class, and it seems to be failing; Cc: @vgvassilev for confirmation.

I also saw your library's structure and if I am not wrong, this function implementation (det_of_minor) is a generic function for all implementations, right? If so, I guess it won't be possible to modify the function's interface to take a C-style array instead of the vector class to make this work with clad.

[vgvassilev]

cc: @parth-07 on the level of std::vector support.

[bradbell]

My goal in cmpad is to have a single set of source code for algorithms that can be differentiated by many packages and to compare the speed of those algorithms. I also want cmpad to provide an example implementation of derivative calculations for many packages. I want to make it as easy as possible to add a package and to make it possible to test any subset of packages. For example, if a clad developer wanted to improve the speed of the clad implementation, she or he would only need to attach clad to cmpad to work on that (not worry about the other packages).

I could have made the interface use raw C vectors together with a size argument, but I did not do this because indexing off the end of the vector is a very common problem. CppAD::vector acts very much like std::vector and checks for indexing off the end of the vector. So I decided to use cmpad::vector as either CppAD::vector or std::vector, when one chooses CppAD::vector (and NDEBUG is defined) it checks for indexing errors.

I think what is needed on the clad side is better documentation of what operations can be differentiated. This would expand with time as more and more operations are handled. I suspect that this is determined by what LLVM operations are currently handled by clad ? If so, I would embed that documentation in the source code that handles the operations .

[bradbell]

It has been a while so I thought I would check and see if clad really supports function objects, demos/Functor seems to indicate that it does.

/*
Expect that this  should work for all of 0 or 1 choices of USE_FUN_OBJ and FUN_OBJ_N_CONST.

1. If USE_FUN_OBJ is 0, this program executes without error.
2. If USE_FUN_OBJ is 1 and FUN_OBJ_N_CONST is 1,
   program exits with messge 'cannot assign to non-static data member 'n_' ...
3. If USE_FUN_OBJ is 1 and FUN_OBJ_N_CONST is 0,
   program exits with a segementation fault.
-------------------------------------------------------------------------------
*/
# include "clad/Differentiator/Differentiator.h"
# include <iostream>
//
// USE_FUN_OBJ
# define USE_FUN_OBJ 1
//
// FUN_OBJ_N_CONST
# define FUN_OBJ_N_CONST 0
//
// f(x, n)
double f(double *x, size_t n)
{  double ret = 0;
   for (size_t i = 0; i < n; i++) {
      ret += (i+1) * x[i];
   }
   return ret;
}

// f_obj
struct f_obj {
   // n_
# if FUN_OBJ_N_CONST
   const size_t n_;
# else
   size_t n_;
# endif
   //
   // ctor
   f_obj(size_t n) : n_(n)
   { }
   //
   // f_obj(x)
   double operator()(double* x)
   {  return f(x, n_); }
};

int main() {
   // x, n, F
   double x[] = {1.0, 2.0, 3.0};
   size_t n   = sizeof(x) / sizeof(x[0]);
   f_obj F(n);
   //
   // f(x, n)
   std::cout << "f(x, n) = " << F(x) << "\n";
   //
   // dx
# if USE_FUN_OBJ
   auto f_dx = clad::gradient(F, "x");
# else
   auto f_dx = clad::gradient(f, "x");
# endif
   double dx[n];
   clad::array_ref<double> dx_ref(dx, n);
   f_dx.execute(x, n, dx_ref);
   //
   // ok
   bool ok    = true;
   for(size_t i = 0; i < n; ++i)
   {  std::cout << "dx[" << i << "] = " << dx[i] << "\n";
      ok &= dx[i] == i+1;
   }
   if( ok )
   {  std::cout << "temp: OK\n";
      return 0;
   }
   std::cout << "temp: Error\n";
   return 1;
}