Enable Clad in ADBench

ADBench/run-all.ps1

@@ -29,52 +29,52 @@ https://github.com/awf/ADBench
 
 #>
 
-param(# Which build to test.  
+param(# Which build to test.
       # Builds should leave a script file 'cmake-vars-$buildtype.ps1' in the ADBench directory.
       # which sets $bindir to the build directory,
       # And if only some D,K are valid for GMM, sets $gmm_d_vals, and $gmm_k_vals
       [string]$buildtype="Release",
-      
-      # Estimated time of accurate result achievement. 
-      # A runner cyclically reruns measured function until total time becomes more than that value. 
+
+      # Estimated time of accurate result achievement.
+      # A runner cyclically reruns measured function until total time becomes more than that value.
       # Supported only by the benchmark runner-based tools
       # (those with ToolType cpp, dotnet, julia, or python).
       [double]$minimum_measurable_time = 0.5,
 
-      # Maximum number of times to run the function for timing                
-      [int]$nruns_f=10,        
-      
+      # Maximum number of times to run the function for timing
+      [int]$nruns_f=10,
+
       # Maximum number of times to run the jacobian for timing
-      [int]$nruns_J=10,        
-      
+      [int]$nruns_J=10,
+
       # How many seconds to wait before we believe we have accurate timings
-      [double]$time_limit=10,  
-      
-      # Kill the test after this many seconds 
+      [double]$time_limit=10,
+
+      # Kill the test after this many seconds
       [double]$timeout=300,
 
       # Kill the test if it consumes more than this many gigabytes of RAM.
       [double]$max_memory_amount_in_gb=[double]::PositiveInfinity,
-      
-      # Where to store the ouput, defaults to tmp/ in the project root 
+
+      # Where to store the ouput, defaults to tmp/ in the project root
       [string]$tmpdir="",
-      
+
       # Repeat tests, even if output file exists
       [switch]$repeat,
-      
+
       # Repeat only failed tests
       [switch]$repeat_failures,
-      
+
       # List of tools to run
       [string[]]$tools=@(),
-      
+
       # Don't delete produced jacobians even if they're accurate
-      [switch]$keep_correct_jacobians, 
-      
+      [switch]$keep_correct_jacobians,
+
       # GMM D values to try.  Must be a subset of the list of
-      # compiled values in ADBench/cmake-vars-$buildtype.ps1                         
+      # compiled values in ADBench/cmake-vars-$buildtype.ps1
       [int[]]$gmm_d_vals_param,
-      
+
       # GMM K values to run.  As above.
       [int[]]$gmm_k_vals_param,
 
@@ -209,7 +209,7 @@ function run_command ($indent, $outfile, $timeout, $cmd) {
                 $mem = [math]::round($mem / (1024 * 1024 * 1024), 2)
                 Write-Host "${indent}Killed due to consuming $mem GB of operating memory"
                 Store-NonFatalError "Process killed due to consuming $mem GB of operating memory`n[$cmd $args]"
-                
+
                 $status = [RunCommandStatus]::OutOfMemory
                 break
             }
@@ -327,16 +327,16 @@ Class Tool {
         This will create a Tool:
         - called "Finite"
         - run from binary executables
-        - runs all four tests 
+        - runs all four tests
         - does not do GMM in separate FULL and SPLIT modes
         - doesn't require separate executables for different GMM sizes
         - does not check the correctness of the computed jacobians
 
         .NOTES
-        $objectives is an enumerable variable, 
+        $objectives is an enumerable variable,
         where each flag determines whether to run a certain objective:
         GMM, BA, HAND, LSTM
-        
+
         #>
 
         $this.name = $name
@@ -358,16 +358,16 @@ Class Tool {
         This will create a Tool:
         - called "Finite"
         - run from binary executables
-        - runs all four tests 
+        - runs all four tests
         - does not do GMM in separate FULL and SPLIT modes
         - doesn't require separate executables for different GMM sizes
         - does not check the correctness of the computed jacobians
 
         .NOTES
-        $objectives is an enumerable variable, 
+        $objectives is an enumerable variable,
         where each flag determines whether to run a certain objective:
         GMM, BA, HAND, LSTM
-        
+
         #>
 
         $this.name = $name
@@ -477,11 +477,11 @@ Class Tool {
         if ($run_command_status -eq [RunCommandStatus]::OutOfMemory) {
             $status = [RunTestStatus]::OutOfMemory
         }
-        
+
         return $status
     }
 
-    # Get postfix for tool output file name 
+    # Get postfix for tool output file name
     [string] get_out_name_postfix([string]$objective) {
         $postfix = $this.name
 
@@ -697,7 +697,10 @@ $tool_descriptors = @(
     [Tool]::new("ManualEigen", "cpp", [ObjectiveType] "GMM, BA, Hand, LSTM", $true, $default_tolerance)
     [Tool]::new("ManualEigenVector", "cpp", [ObjectiveType] "GMM", $true, $default_tolerance)
     [Tool]::new("DiffSharpModule", "dotnet", [ObjectiveType] "GMM, BA, Hand, LSTM", $true, $default_tolerance)
-    [Tool]::new("Tapenade", "cpp", [ObjectiveType] "BA, LSTM, GMM, Hand", $true, $default_tolerance)
+    [Tool]::new("Tapenade", "cpp", [ObjectiveType] "GMM, BA", $true, $default_tolerance)
+    [Tool]::new("Clad", "cpp", [ObjectiveType] "GMM, BA", $true, $default_tolerance)
+    [Tool]::new("CladLatest", "cpp", [ObjectiveType] "GMM, BA", $true, $default_tolerance)
+    [Tool]::new("CladLatestTBR", "cpp", [ObjectiveType] "GMM, BA", $true, $default_tolerance)
     [Tool]::new("PyTorch", "python", [ObjectiveType] "BA, LSTM, GMM, Hand", $true, 1e-7)
     [Tool]::new("TorchScript", "python", [ObjectiveType] "GMM", $true, 1e-7)
     [Tool]::new("Tensorflow", "python", [ObjectiveType] "BA, LSTM, GMM, Hand", $true, $default_tolerance)

src/cpp/modules/CMakeLists.txt

@@ -1,5 +1,8 @@
 project("ADBench_Modules")
 
+# add_subdirectory ("clad")
+add_subdirectory ("cladLatest")
+add_subdirectory ("cladLatestTBR")
 add_subdirectory ("manual")
 add_subdirectory ("manualEigen")
 add_subdirectory ("manualEigenVector")

src/cpp/modules/cladLatest/CMakeLists.txt

@@ -0,0 +1,15 @@
+project("ADBench_CladLatest" CXX)
+
+add_library("CladLatest" MODULE)
+
+SET(clad_install_root "/home/mad-scientist/clubs/open-source/clad/formal/builds/build-llvm-16/inst")
+
+target_sources("CladLatest" PRIVATE "${CMAKE_SOURCE_DIR}/src/cpp/shared/utils.cpp")
+target_sources("CladLatest" PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/ba/ba_grad.cpp")
+target_sources("CladLatest" PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/CladLatestBA.cpp")
+target_sources("CladLatest" PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/gmm/gmm_grad.cpp")
+target_sources("CladLatest" PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/CladLatestGMM.cpp")
+
+target_include_directories("CladLatest" PRIVATE ${clad_install_root}/include)
+target_include_directories("CladLatest" PRIVATE ${CMAKE_SOURCE_DIR}/src/cpp/shared)
+# add_compile_definitions("CladLatest" CLAD_NO_NUM_DIFF)

src/cpp/modules/cladLatest/CladLatestBA.cpp

@@ -0,0 +1,115 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT license.
+
+#include "CladLatestBA.h"
+#include "ba/ba_grad.h"
+
+// This function must be called before any other function.
+void CladLatestBA::prepare(BAInput &&input) {
+  this->input = input;
+  result = {std::vector<double>(2 * this->input.p),
+            std::vector<double>(this->input.p),
+            BASparseMat(this->input.n, this->input.m, this->input.p)};
+
+  reproj_err_d = std::vector<double>(2 * (BA_NCAMPARAMS + 3 + 1));
+  reproj_err_d_row = std::vector<double>(BA_NCAMPARAMS + 3 + 1);
+}
+
+BAOutput CladLatestBA::output() { return result; }
+
+void CladLatestBA::calculate_objective(int times) {
+  for (int i = 0; i < times; i++) {
+    ba_objective(input.n, input.m, input.p, input.cams.data(), input.X.data(),
+                 input.w.data(), input.obs.data(), input.feats.data(),
+                 result.reproj_err.data(), result.w_err.data());
+  }
+}
+
+void CladLatestBA::calculate_jacobian(int times) {
+  for (int i = 0; i < times; i++) {
+    result.J.clear();
+    calculate_reproj_error_jacobian_part();
+    calculate_weight_error_jacobian_part();
+  }
+}
+
+void CladLatestBA::calculate_reproj_error_jacobian_part() {
+  double
+      errb[2]; // stores dY
+               // (i-th element equals to 1.0 for calculating i-th jacobian row)
+
+  double err[2]; // stores fictive result
+                 // (Tapenade doesn't calculate an original function in
+                 // reverse mode)
+
+  clad::array_ref<double> cam_gradient_part(reproj_err_d_row.data(),
+                                            reproj_err_d_row.size());
+  clad::array_ref<double> x_gradient_part(
+      reproj_err_d_row.data() + BA_NCAMPARAMS,
+      reproj_err_d_row.size() - BA_NCAMPARAMS);
+  clad::array_ref<double> weight_gradient_part(
+      reproj_err_d_row.data() + BA_NCAMPARAMS + 3,
+      reproj_err_d_row.size() - BA_NCAMPARAMS - 3);
+  std::vector<double> feats_gradient_part(2, 0);
+  clad::array_ref<double> feats_gradient_part_ref(feats_gradient_part.data(),
+                                                  feats_gradient_part.size());
+  clad::array_ref<double> errb_ref(errb, 2);
+
+  for (int i = 0; i < input.p; i++) {
+    std::fill(reproj_err_d_row.begin(), reproj_err_d_row.end(), 0);
+    std::fill(feats_gradient_part.begin(), feats_gradient_part.end(), 0);
+
+    int camIdx = input.obs[2 * i + 0];
+    int ptIdx = input.obs[2 * i + 1];
+
+    // calculate first row
+    errb[0] = 1.0;
+    errb[1] = 0.0;
+    computeReprojError_pullback(
+        &input.cams[camIdx * BA_NCAMPARAMS], &input.X[ptIdx * 3], &input.w[i],
+        &input.feats[i * 2], err, cam_gradient_part, x_gradient_part,
+        weight_gradient_part, feats_gradient_part_ref, errb_ref);
+
+    // fill first row elements
+    for (int j = 0; j < BA_NCAMPARAMS + 3 + 1; j++) {
+      reproj_err_d[2 * j] = reproj_err_d_row[j];
+    }
+
+    std::fill(reproj_err_d_row.begin(), reproj_err_d_row.end(), 0);
+    std::fill(feats_gradient_part.begin(), feats_gradient_part.end(), 0);
+    // calculate second row
+    errb[0] = 0.0;
+    errb[1] = 1.0;
+    computeReprojError_pullback(
+        &input.cams[camIdx * BA_NCAMPARAMS], &input.X[ptIdx * 3], &input.w[i],
+        &input.feats[i * 2], err, cam_gradient_part, x_gradient_part,
+        weight_gradient_part, feats_gradient_part_ref, errb);
+
+    // fill second row elements
+    for (int j = 0; j < BA_NCAMPARAMS + 3 + 1; j++) {
+      reproj_err_d[2 * j + 1] = reproj_err_d_row[j];
+    }
+
+    result.J.insert_reproj_err_block(i, camIdx, ptIdx, reproj_err_d.data());
+  }
+}
+
+void CladLatestBA::calculate_weight_error_jacobian_part() {
+  for (int j = 0; j < input.p; j++) {
+    double wb = 0; // stores calculated derivative
+
+    double err = 0.0; // stores fictive result
+                      // (Tapenade doesn't calculate an original function in
+                      // reverse mode)
+
+    double errb = 1.0; // stores dY
+                       // (equals to 1.0 for derivative calculation)
+
+    computeZachWeightError_pullback(&input.w[j], &err, &wb, &errb);
+    result.J.insert_w_err_block(j, wb);
+  }
+}
+
+extern "C" DLL_PUBLIC ITest<BAInput, BAOutput> *get_ba_test() {
+  return new CladLatestBA();
+}

src/cpp/modules/cladLatest/CladLatestBA.h

@@ -0,0 +1,39 @@
+// Copyright (c) Microsoft Corporation.
+// Licensed under the MIT license.
+
+#pragma once
+
+#include "../../shared/ITest.h"
+#include "../../shared/BAData.h"
+#include "../../shared/defs.h"
+
+#include "ba/ba.h"
+
+#include <vector>
+
+class CladLatestBA : public ITest<BAInput, BAOutput>
+{
+private:
+    BAInput input;
+    BAOutput result;
+    std::vector<double> state;
+
+    // buffer for reprojection error jacobian part holding (column-major)
+    std::vector<double> reproj_err_d;
+
+    // buffer for reprojection error jacobian block row holding
+    std::vector<double> reproj_err_d_row;
+public:
+    // This function must be called before any other function.
+    virtual void prepare(BAInput&& input) override;
+
+    virtual void calculate_objective(int times) override;
+    virtual void calculate_jacobian(int times) override;
+    virtual BAOutput output() override;
+
+    ~CladLatestBA() {}
+
+private:
+    void calculate_weight_error_jacobian_part();
+    void calculate_reproj_error_jacobian_part();
+};