Add a backup central difference numerical differentiation method.

Clad can now build calls to a "forward" (single-arg) and "reverse" (multi-arg) numerical diff function. Numerical diff can be used in the following contexts:
- With clad forward and reverse mode for all in-built scalar and non scalar types.
- With clad and standalone for "forward" and "reverse" derivatives.
- Standalone for in-built scalar and non-scalar types.
- Standalone with support for user defined types as input (both value and pointer forms).
- Standalone to differentiate overloaded operators.
- Standalone to differentiate functors.

Numerical diff error estimates can be printed with the help of the '-fprint-num-diff-errors' flag, and numerical diff may be disabled using the '-DCLAD_NO_NUM_DIFF' flag during compilation.

demos/CustomTypeNumDiff.cpp

@@ -0,0 +1,149 @@
+//--------------------------------------------------------------------*- C++ -*-
+// clad - The C++ Clang-based Automatic Differentiator
+//
+// A demo, describing how users can numerically differentiate functions that
+// take user-defined types as parameters using clad.
+//
+// author: Garima Singh
+//----------------------------------------------------------------------------//
+
+// To run the demo please type:
+// path/to/clang  -Xclang -add-plugin -Xclang clad -Xclang -load -Xclang \
+// path/to/libclad.so  -I../include/ -x c++ -std=c++11 CustomTypeNumDiff.cpp
+//
+// A typical invocation would be:
+// ../../../../obj/Debug+Asserts/bin/clang  -Xclang -add-plugin -Xclang clad \
+// -Xclang -load -Xclang ../../../../obj/Debug+Asserts/lib/libclad.dylib     \
+// -I../include/ -x c++ -std=c++11 CustomTypeNumDiff.cpp
+
+// Necessary for clad to work include
+#include "clad/Differentiator/Differentiator.h"
+
+#include <iostream> // for std::* functions.
+
+// A user defined class to emulate fixed point decimal numbers.
+// It stores the significand and inverted exponent as integers.
+// for a number of the form ax10e-b
+// 'number' refers to 'a' and 'scale' refers to 'b'.
+class FixedPointFloat {
+  long int number;
+  unsigned scale;
+
+public:
+  FixedPointFloat(unsigned num, unsigned n) {
+    // If there are any trailing zeros to num and a non-zero scale value
+    // normalize the number so that we do not do unnecessary multiplications.
+    while (!num % 10 && n) {
+      num = std::floor(num / 10);
+      n--;
+    }
+    number = num;
+    scale = n;
+  }
+  // A function to get the double equivalent of the decimal number stored.
+  double getDouble(unsigned num, unsigned scale) {
+    return (double)num / std::pow(10, scale);
+  }
+  // A function to scale the value of a given number, we use this to
+  // facilitate evaluation of different operations on this class.
+  unsigned getScaled(unsigned num, unsigned scale) {
+    for (unsigned i = 0; i < scale; i++) {
+      num *= 10;
+    }
+    return num;
+  }
+  // Gets the double representation of the number represented by the
+  // calling object.
+  double getNum() { return getDouble(number, scale); }
+  // Some operator overloads for this class...
+  FixedPointFloat operator+(FixedPointFloat a) {
+    unsigned xa, yb;
+    // Here, let us understand this with an example.
+    // consider two numbers -> 0.005 and 0.00002 ~ 5x10e-3 and 2x10e-5
+    // converting them to FixedPointFloat will give us a pair like so:
+    // (5, 3) and (2, 5)
+    // To add the above numbers, we follow the following...
+    // Get the max scale value (i.e the minimum exponent value)
+    // in our case, it is 5.
+    unsigned maxScale = std::max(a.scale, scale);
+    // Scale both numbers with the remaining scale value.
+    // essentially, do the following:
+    // (5x10e-3 + 2x10e-5) * (10e5/10e5)
+    // (5x10e(-3 + 5) + 2x10e(-5+5)) * (1/10e5)
+    // => (500 + 2) * (1/10e5)
+    xa = getScaled(a.number, maxScale - a.scale); // = 500
+    yb = getScaled(number, maxScale - scale);     // = 2
+    // => (500 + 2) * (1/10e5) = 502x10e-5
+    return FixedPointFloat(xa + yb, maxScale); // = (502, 5)
+  }
+  // Similar to above, but here, we just subtract the number instead of adding
+  // it.
+  FixedPointFloat operator-(FixedPointFloat a) {
+    unsigned xa, yb;
+    unsigned maxScale = std::max(a.scale, scale);
+    xa = getScaled(a.number, maxScale - a.scale);
+    yb = getScaled(number, maxScale - scale);
+    return FixedPointFloat(xa - yb, maxScale);
+  }
+
+  FixedPointFloat operator*(FixedPointFloat a) {
+    unsigned xa, yb;
+    // This is fairly straight forward. Let us take the same example as before
+    // 0.005 and 0.00002 ~ 5x10e-3 and 2x10e-5
+    // This operation is equivalent to
+    // = 5x10e-3 x 2x10e-5
+    // = (5x2) x 10e(-5-3)
+    // = 10x10e-8 (which is eventually reformed to 1x10e-7)
+    return FixedPointFloat(a.number * number, a.scale + scale); // = (1, 7)
+  }
+};
+
+// This function is essential if we want to differentiate a function with
+// user-defined data types as arguments. This function describes the "rules"
+// of how to differentiate user-defined data types by specifying how to
+// introduce a small change (h) in the object of the user-defined data type.
+// Details on how to overload this function are provided in the docs.
+FixedPointFloat
+updateIndexParamValue(FixedPointFloat arg, std::size_t idx, std::size_t currIdx,
+                      int multiplier, numerical_diff::precision& h_val,
+                      std::size_t n = 0, std::size_t i = 0) {
+  if (idx == currIdx) {
+    // Here we just introduce an h of 0.000001 to all our FixedPointFloat
+    // numbers.
+    FixedPointFloat h(1, 5);
+    h_val = (h_val == 0) ? h.getNum() : h_val;
+    FixedPointFloat Fmultiplier(multiplier, 0);
+    return arg + Fmultiplier * h;
+  }
+  return arg;
+}
+
+// A simple multiplication function.
+// currently we can only numerically differentiate the following types:
+// - all scalar numeric types
+// - all types with basic arithmetic operators overloaded.
+// - all types that are implicitly convertible to some numeric type.
+double func(FixedPointFloat x, FixedPointFloat y) { return (x * y).getNum(); }
+
+int main() {
+  // Define some values which will reflect the derivative later.
+  double dx = 0, dy = 0;
+  // Define our inputs: 3x10e-3 and 7x10e-2 or 0.003 and 0.07.
+  FixedPointFloat x(3, 3), y(7, 2);
+  // Push the dx and dy values to a tape created by us.
+  // This is how we return the derivative with respect to all arguments.
+  // The order of being placed in this tape should be the same as the order of
+  // the arguments being passed to the function.
+  clad::tape<clad::array_ref<
+      double /*This should be the return value of the function you want to differentiate.*/>>
+      grad = {};
+  // Place the l-value reference of the variables in the tape.
+  grad.emplace_back(&dx);
+  grad.emplace_back(&dy);
+  // Finally, call the numerical diff method, keep the order of the arguments to
+  // the function in mind!
+  numerical_diff::central_difference(func, grad, /*printErrors=*/0, x, y);
+  // Finally print the results!
+  std::cout << "Result of df/dx is = " << dx << "\nResult of df/dx is = " << dy
+            << std::endl;
+}

include/clad/Differentiator/DerivativeBuilder.h

@@ -97,17 +97,36 @@ namespace clad {
     clang::NamespaceDecl* m_BuiltinDerivativesNSD;
     /// A reference to the model to use for error estimation (if any).
     std::unique_ptr<FPErrorEstimationModel> m_EstModel = nullptr;
+    clang::NamespaceDecl* m_NumericalDiffNSD;
+    /// A flag to keep track of whether error diagnostics are requested by user
+    /// for numerical differentiation.
+    bool m_PrintNumericalDiffErrorDiag = false;
     DeclWithContext cloneFunction(const clang::FunctionDecl* FD,
-                                  clad::VisitorBase VB,
-                                  clang::DeclContext* DC,
+                                  clad::VisitorBase VB, clang::DeclContext* DC,
                                   clang::Sema& m_Sema,
                                   clang::ASTContext& m_Context,
                                   clang::SourceLocation& noLoc,
                                   clang::DeclarationNameInfo name,
                                   clang::QualType functionType);
+    /// Looks for a suitable overload for a given function.
+    ///
+    /// \param[in] DNI The identification information of the function overload
+    /// to be found.
+    /// \param[in] CallArgs The call args to be used to resolve to the
+    /// correct overload.
+    /// \param[in] forCustomDerv A flag to keep track of which
+    /// namespace we should look in for the overloads.
+    /// \param[in] namespaceShouldExist A flag to enforce assertion failure
+    /// if the overload function namespace was not found. If false and
+    /// the function containing namespace was not found, nullptr is returned.
+    ///
+    /// \returns The call expression if a suitable function overload was found,
+    /// null otherwise.
     clang::Expr*
     findOverloadedDefinition(clang::DeclarationNameInfo DNI,
-                             llvm::SmallVectorImpl<clang::Expr*>& CallArgs);
+                             llvm::SmallVectorImpl<clang::Expr*>& CallArgs,
+                             bool forCustomDerv = true,
+                             bool namespaceShouldExist = true);
     bool noOverloadExists(clang::Expr* UnresolvedLookup,
                           llvm::MutableArrayRef<clang::Expr*> ARargs);
     /// Shorthand to issues a warning or error.
@@ -130,6 +149,19 @@ namespace clad {
     /// an in-built one (TaylorApprox) or one provided by the user.
     void
     SetErrorEstimationModel(std::unique_ptr<FPErrorEstimationModel> estModel);
+    /// Fuction to set the error diagnostic printing value for numerical
+    /// differentiation.
+    ///
+    /// \param[in] \c value The new value to be set.
+    void setNumDiffErrDiag(bool value) {
+      m_PrintNumericalDiffErrorDiag = value;
+    }
+    /// Function to return if clad should emit error information for numerical
+    /// differentiation.
+    ///
+    /// \returns The flag  that controls printing of error information for
+    /// numerical differentiation.
+    bool shouldPrintNumDiffErrs() { return m_PrintNumericalDiffErrorDiag; }
     ///\brief Produces the derivative of a given function
     /// according to a given plan.
     ///

include/clad/Differentiator/Differentiator.h

@@ -12,6 +12,7 @@
 #include "BuiltinDerivatives.h"
 #include "CladConfig.h"
 #include "FunctionTraits.h"
+#include "NumericalDiff.h"
 #include "Tape.h"
 
 #include <assert.h>
@@ -58,6 +59,15 @@ namespace clad {
     return val;
   }
 
+  /// Add value to the end of the tape, return the same value.
+  /// A specialization for clad::array_ref types to use in reverse mode.
+  template <typename T, typename U>
+  CUDA_HOST_DEVICE clad::array_ref<T> push(tape<clad::array_ref<T>>& to,
+                                           U val) {
+    to.emplace_back(val);
+    return val;
+  }
+
   /// Remove the last value from the tape, return it.
   template <typename T>
   CUDA_HOST_DEVICE T pop(tape<T>& to) {