add test for autodiff about forwarding func

.gitignore

@@ -0,0 +1,2 @@
+.cache
+build

.vscode/launch.json

@@ -0,0 +1,21 @@
+{
+    // Use IntelliSense to learn about possible attributes.
+    // Hover to view descriptions of existing attributes.
+    // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
+    "version": "0.2.0",
+    "configurations": [
+        {
+            "type": "vgdb",
+            "request": "launch",
+            "name": "C/C++ Debug",
+            "program": "${command:cmake.launchTargetPath}",
+            "args": [],
+            // "cwd": "${command:cmake.getLaunchTargetDirectory}",
+            "cwd": "${workspaceFolder}",
+            "env": {
+                "name": "PATH",
+                "value": "${env:PATH}:${command:cmake.getLaunchTargetDirectory}"
+            }
+        }
+    ]
+}

CMakeLists.txt

@@ -0,0 +1,24 @@
+cmake_minimum_required(VERSION 3.12)
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
+set(CMAKE_BUILD_TYPE Debug)
+
+set(conda_env_path "/share/software/anaconda3/envs/unisolver")
+if(DEFINED CMAKE_PREFIX_PATH)
+    list(APPEND CMAKE_PREFIX_PATH ${conda_path})
+else()
+    set(CMAKE_PREFIX_PATH ${conda_env_path})
+endif()
+message(STATUS "cmake_prefix_path: ${CMAKE_PREFIX_PATH}")
+
+project(test_autodiff LANGUAGES CXX)
+
+find_package(autodiff)
+find_package(Eigen3 3.4 REQUIRED)
+include_directories(${EIGEN3_INCLUDE_DIRS})
+
+
+add_subdirectory(test)

test/CMakeLists.txt

@@ -0,0 +1,7 @@
+file(GLOB V_GLOB LIST_DIRECTORIES true "*")
+# list(FILTER V_GLOB EXCLUDE REGEX ".*ThirdParty")
+foreach(item ${V_GLOB})
+	if(IS_DIRECTORY ${item})
+		add_subdirectory(${item})
+	endif()
+endforeach()

test/forward/10jacobian_matrix_of_vector_func_using_memory_maps.cpp

@@ -0,0 +1,37 @@
+// NOTE: 使用内存映射的向量函数的雅可比矩阵 
+
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/real.hpp>
+#include <autodiff/forward/real/eigen.hpp>
+using namespace autodiff;
+
+// The vector function for which the Jacobian is needed
+VectorXreal f(const VectorXreal& x)
+{
+    return x * x.sum();
+}
+
+int main()
+{
+    using Eigen::Map;
+    using Eigen::MatrixXd;
+
+    VectorXreal x(5);                           // the input vector x with 5 variables
+    x << 1, 2, 3, 4, 5;                         // x = [1, 2, 3, 4, 5]
+    double y[25];                               // the output Jacobian as a flat array
+
+    VectorXreal F;                              // the output vector F = f(x) evaluated together with Jacobian matrix below
+    Map<MatrixXd> J(y, 5, 5);                   // the output Jacobian dF/dx mapped onto the flat array
+
+    jacobian(f, wrt(x), at(x), F, J);           // evaluate the output vector F and the Jacobian matrix dF/dx
+
+    std::cout << "F = \n" << F << std::endl;    // print the evaluated output vector F
+    std::cout << "J = \n" << J << std::endl;    // print the evaluated Jacobian matrix dF/dx
+    std::cout << "y = ";                        // print the flat array
+    for(int i = 0 ; i < 25 ; ++i)
+        std::cout << y[i] << " ";
+    std::cout << std::endl;
+}

test/forward/11higher_order_der_of_scalar_func.cpp

@@ -0,0 +1,40 @@
+// NOTE：标量函数的高阶交叉导
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/dual.hpp>
+using namespace autodiff;
+
+// The multi-variable function for which higher-order derivatives are needed (up to 4th order)
+dual4th f(dual4th x, dual4th y, dual4th z)
+{
+    return 1 + x + y + z + x*y + y*z + x*z + x*y*z + exp(x/y + y/z);
+}
+
+int main()
+{
+    dual4th x = 1.0;
+    dual4th y = 2.0;
+    dual4th z = 3.0;
+
+    auto [u0, ux, uxy, uxyx, uxyxz] = derivatives(f, wrt(x, y, x, z), at(x, y, z));
+
+    std::cout << "u0 = " << u0 << std::endl;       // print the evaluated value of u = f(x, y, z)
+    std::cout << "ux = " << ux << std::endl;       // print the evaluated derivative du/dx
+    std::cout << "uxy = " << uxy << std::endl;     // print the evaluated derivative d²u/dxdy
+    std::cout << "uxyx = " << uxyx << std::endl;   // print the evaluated derivative d³u/dx²dy
+    std::cout << "uxyxz = " << uxyxz << std::endl; // print the evaluated derivative d⁴u/dx²dydz
+}
+
+/*-------------------------------------------------------------------------------------------------
+=== Note ===
+---------------------------------------------------------------------------------------------------
+In most cases, dual can be replaced by real, as commented in other examples.
+However, computing higher-order cross derivatives has definitely to be done
+using higher-order dual types (e.g., dual3rd, dual4th)! This is because real
+types (e.g., real2nd, real3rd, real4th) are optimally designed for computing
+higher-order directional derivatives.
+对于计算高阶交叉导数,必须使用更高阶的 dual 类型(如 dual3rd、dual4th 等)。
+否则可以使用real
+-------------------------------------------------------------------------------------------------*/

test/forward/12higher_order_directional_der_of_scalar_func.cpp

@@ -0,0 +1,36 @@
+// NOTE: 标量函数的高阶方向导数
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/real.hpp>
+using namespace autodiff;
+
+// The multi-variable function for which higher-order derivatives are needed (up to 4th order)
+real4th f(real4th x, real4th y, real4th z)
+{
+    return sin(x) * cos(y) * exp(z);
+}
+
+int main()
+{
+    real4th x = 1.0;
+    real4th y = 2.0;
+    real4th z = 3.0;
+
+    auto dfdv = derivatives(f, along(1.0, 1.0, 2.0), at(x, y, z)); // the directional derivatives of f along direction v = (1, 1, 2) at (x, y, z) = (1, 2, 3)
+
+    std::cout << "dfdv[0] = " << dfdv[0] << std::endl; // print the evaluated 0th order directional derivative of f along v (equivalent to f(x, y, z))
+    std::cout << "dfdv[1] = " << dfdv[1] << std::endl; // print the evaluated 1st order directional derivative of f along v
+    std::cout << "dfdv[2] = " << dfdv[2] << std::endl; // print the evaluated 2nd order directional derivative of f along v
+    std::cout << "dfdv[3] = " << dfdv[3] << std::endl; // print the evaluated 3rd order directional derivative of f along v
+    std::cout << "dfdv[4] = " << dfdv[4] << std::endl; // print the evaluated 4th order directional derivative of f along v
+}
+
+/*-------------------------------------------------------------------------------------------------
+=== Note ===
+---------------------------------------------------------------------------------------------------
+This example would also work if dual was used instead of real. However, real
+types are your best option for directional derivatives, as they were optimally
+designed for this kind of derivatives.
+-------------------------------------------------------------------------------------------------*/

test/forward/13higher_order_directional_der_of_vector_func.cpp

@@ -0,0 +1,73 @@
+// NOTE: 向量的高阶方向导数
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/real.hpp>
+#include <autodiff/forward/real/eigen.hpp>
+using namespace autodiff;
+
+// The vector function for which higher-order directional derivatives are needed (up to 4th order).
+ArrayXreal4th f(const ArrayXreal4th& x, real4th p)
+{
+    return p * x.log();
+}
+
+int main()
+{
+    using Eigen::ArrayXd;
+
+    ArrayXreal4th x(5); // the input vector x
+    x << 1.0, 2.0, 3.0, 4.0, 5.0;
+
+    real4th p = 7.0; // the input parameter p = 1
+
+    ArrayXd vx(5); // the direction vx in the direction vector v = (vx, vp)
+    vx << 1.0, 1.0, 1.0, 1.0, 1.0;
+
+    double vp = 1.0; // the direction vp in the direction vector v = (vx, vp)
+
+    auto dfdv = derivatives(f, along(vx, vp), at(x, p)); // the directional derivatives of f along direction v = (vx, vp) at (x, p)
+
+    std::cout << std::scientific << std::showpos;
+    std::cout << "Directional derivatives of f along v = (vx, vp) up to 4th order:" << std::endl;
+    std::cout << "dfdv[0] = " << dfdv[0].transpose() << std::endl; // print the evaluated 0th order directional derivative of f along v (equivalent to f(x, p))
+    std::cout << "dfdv[1] = " << dfdv[1].transpose() << std::endl; // print the evaluated 1st order directional derivative of f along v
+    std::cout << "dfdv[2] = " << dfdv[2].transpose() << std::endl; // print the evaluated 2nd order directional derivative of f along v
+    std::cout << "dfdv[3] = " << dfdv[3].transpose() << std::endl; // print the evaluated 3rd order directional derivative of f along v
+    std::cout << "dfdv[4] = " << dfdv[4].transpose() << std::endl; // print the evaluated 4th order directional derivative of f along v
+    std::cout << std::endl;
+
+    double t = 0.1; // the step length along direction v = (vx, vp) used to compute 4th order Taylor estimate of f
+
+    ArrayXreal4th u = f(x + t * vx, p + t * vp); // start from (x, p), walk a step length t = 0.1 along direction v = (vx, vp) and evaluate f at this current point
+
+    ArrayXd utaylor = dfdv[0] + t*dfdv[1] + (t*t/2.0)*dfdv[2] + (t*t*t/6.0)*dfdv[3] + (t*t*t*t/24.0)*dfdv[4]; // evaluate the 4th order Taylor estimate of f along direction v = (vx, vp) at a step length of t = 0.1
+
+    std::cout << "Comparison between exact evaluation and 4th order Taylor estimate:" << std::endl;
+    std::cout << "u(exact)  = " << u.transpose() << std::endl;
+    std::cout << "u(taylor) = " << utaylor.transpose() << std::endl;
+}
+
+/*-------------------------------------------------------------------------------------------------
+=== Output ===
+---------------------------------------------------------------------------------------------------
+Directional derivatives of f along v = (vx, vp) up to 4th order:
+dfdv[0] = +0.000000e+00 +4.852030e+00 +7.690286e+00 +9.704061e+00 +1.126607e+01
+dfdv[1] = +7.000000e+00 +4.193147e+00 +3.431946e+00 +3.136294e+00 +3.009438e+00
+dfdv[2] = -5.000000e+00 -7.500000e-01 -1.111111e-01 +6.250000e-02 +1.200000e-01
+dfdv[3] = +1.100000e+01 +1.000000e+00 +1.851852e-01 +3.125000e-02 -8.000000e-03
+dfdv[4] = -3.400000e+01 -1.625000e+00 -2.222222e-01 -3.906250e-02 -3.200000e-03
+
+Comparison between exact evaluation and 4th order Taylor estimate:
+u(exact)  = +6.767023e-01 +5.267755e+00 +8.032955e+00 +1.001801e+01 +1.156761e+01
+u(taylor) = +6.766917e-01 +5.267755e+00 +8.032955e+00 +1.001801e+01 +1.156761e+01
+-------------------------------------------------------------------------------------------------*/
+
+/*-------------------------------------------------------------------------------------------------
+=== Note ===
+---------------------------------------------------------------------------------------------------
+This example would also work if dual was used instead of real. However, real
+types are your best option for directional derivatives, as they were optimally
+designed for this kind of derivatives.
+-------------------------------------------------------------------------------------------------*/

test/forward/14taylor_series_of_scalar_func_along_a_direction.cpp

@@ -0,0 +1,42 @@
+// NOTE: 标量函数沿某个方向的泰勒级数
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/real.hpp>
+#include <autodiff/forward/real/eigen.hpp>
+using namespace autodiff;
+
+// The scalar function for which a 4th order directional Taylor series will be computed.
+real4th f(const real4th& x, const real4th& y, const real4th& z)
+{
+    return sin(x * y) * cos(x * z) * exp(z);
+}
+
+int main()
+{
+    real4th x = 1.0;                                       // the input vector x
+    real4th y = 2.0;                                       // the input vector y
+    real4th z = 3.0;                                       // the input vector z
+
+    auto g = taylorseries(f, along(1, 1, 2), at(x, y, z)); // the function g(t) as a 4th order Taylor approximation of f(x + t, y + t, z + 2t)
+
+    double t = 0.1;                                        // the step length used to evaluate g(t), the Taylor approximation of f(x + t, y + t, z + 2t)
+
+    real4th u = f(x + t, y + t, z + 2*t);                  // the exact value of f(x + t, y + t, z + 2t)
+
+    double utaylor = g(t);                                 // the 4th order Taylor estimate of f(x + t, y + t, z + 2t)
+
+    std::cout << std::fixed;
+    std::cout << "Comparison between exact evaluation and 4th order Taylor estimate of f(x + t, y + t, z + 2t):" << std::endl;
+    std::cout << "u(exact)  = " << u << std::endl;
+    std::cout << "u(taylor) = " << utaylor << std::endl;
+}
+
+/*-------------------------------------------------------------------------------------------------
+=== Output ===
+---------------------------------------------------------------------------------------------------
+Comparison between exact evaluation and 4th order Taylor estimate of f(x + t, y + t, z + 2t):
+u(exact)  = -16.847071
+u(taylor) = -16.793986
+-------------------------------------------------------------------------------------------------*/

test/forward/15taylor_series_of_vector_func_along_a_direction.cpp

@@ -0,0 +1,46 @@
+// NOTE: 向量函数沿某个方向的泰勒级数
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/real.hpp>
+#include <autodiff/forward/real/eigen.hpp>
+using namespace autodiff;
+
+// The vector function for which a 4th order directional Taylor series will be computed.
+ArrayXreal4th f(const ArrayXreal4th& x)
+{
+    return x.sin() / x;
+}
+
+int main()
+{
+    using Eigen::ArrayXd;
+
+    ArrayXreal4th x(5);                        // the input vector x
+    x << 1.0, 2.0, 3.0, 4.0, 5.0;
+
+    ArrayXd v(5);                              // the direction vector v
+    v << 1.0, 1.0, 1.0, 1.0, 1.0;
+
+    auto g = taylorseries(f, along(v), at(x)); // the function g(t) as a 4th order Taylor approximation of f(x + t·v)
+
+    double t = 0.1;                            // the step length used to evaluate g(t), the Taylor approximation of f(x + t·v)
+
+    ArrayXreal4th u = f(x + t * v);            // the exact value of f(x + t·v)
+
+    ArrayXd utaylor = g(t);                    // the 4th order Taylor estimate of f(x + t·v)
+
+    std::cout << std::fixed;
+    std::cout << "Comparison between exact evaluation and 4th order Taylor estimate of f(x + t·v):" << std::endl;
+    std::cout << "u(exact)  = " << u.transpose() << std::endl;
+    std::cout << "u(taylor) = " << utaylor.transpose() << std::endl;
+}
+
+/*-------------------------------------------------------------------------------------------------
+=== Output ===
+---------------------------------------------------------------------------------------------------
+Comparison between exact evaluation and 4th order Taylor estimate of f(x + t·v):
+u(exact)  =  0.810189  0.411052  0.013413 -0.199580 -0.181532
+u(taylor) =  0.810189  0.411052  0.013413 -0.199580 -0.181532
+-------------------------------------------------------------------------------------------------*/

test/forward/1single_variable_function.cpp

@@ -0,0 +1,25 @@
+// NOTE: 单变量函数的导数
+
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/dual.hpp>
+using namespace autodiff;
+
+// The single-variable function for which derivatives are needed
+dual f(dual x)
+{
+    return 1 + x + x*x + 1/x + log(x);
+}
+
+int main()
+{
+    dual x = 2.0;                                 // the input variable x
+    dual u = f(x);                                // the output variable u
+
+    double dudx = derivative(f, wrt(x), at(x));   // evaluate the derivative du/dx
+
+    std::cout << "u = " << u << std::endl;        // print the evaluated output u
+    std::cout << "du/dx = " << dudx << std::endl; // print the evaluated derivative du/dx
+}

test/forward/2single_variable_function_with_custom_scalar.cpp

@@ -0,0 +1,32 @@
+// NOTE: 使用自定义标量（复数）的单变量函数的导数
+
+// C++ includes
+#include <iostream>
+#include <complex>
+using namespace std;
+
+// autodiff include
+#include <autodiff/forward/dual.hpp>
+using namespace autodiff;
+
+// Specialize isArithmetic for complex to make it compatible with dual
+namespace autodiff::detail {
+    template<typename T>
+    struct ArithmeticTraits<complex<T>>: ArithmeticTraits<T> {};
+} // autodiff::detail
+
+using cxdual = Dual<complex<double>, complex<double>>;
+
+cxdual f(cxdual x){
+    return 1+x+x*x+1/x+log(x);
+}
+
+int main(){
+    cxdual x = 2.0; // the input variable x
+    cxdual u = f(x); // the output variable u
+    cxdual dudx = derivative(f, wrt(x), at(x)); // evaluate the derivative du/dx
+
+    cout << "u = " << u << endl;
+    cout << "du/dx = " << dudx << endl;
+    return 0;
+}

test/forward/3multi_varialbe_function.cpp

@@ -0,0 +1,32 @@
+// NOTE: 多变量函数的导数
+
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/dual.hpp>
+using namespace autodiff;
+
+// The multi-variable function for which derivatives are needed
+dual f(dual x, dual y, dual z)
+{
+    return 1 + x + y + z + x*y + y*z + x*z + x*y*z + exp(x/y + y/z);
+}
+
+int main()
+{
+    dual x = 1.0;
+    dual y = 2.0;
+    dual z = 3.0;
+
+    dual u = f(x, y, z);
+
+    double dudx = derivative(f, wrt(x), at(x, y, z));
+    double dudy = derivative(f, wrt(y), at(x, y, z));
+    double dudz = derivative(f, wrt(z), at(x, y, z));
+
+    std::cout << "u = " << u << std::endl;         // print the evaluated output u = f(x, y, z)
+    std::cout << "du/dx = " << dudx << std::endl;  // print the evaluated derivative du/dx
+    std::cout << "du/dy = " << dudy << std::endl;  // print the evaluated derivative du/dy
+    std::cout << "du/dz = " << dudz << std::endl;  // print the evaluated derivative du/dz
+}

test/forward/4multi_variable_function_with_params.cpp

@@ -0,0 +1,59 @@
+// NOTE: 带参数的多变量函数的导数
+
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/dual.hpp>
+using namespace autodiff;
+
+// A type defining parameters for a function of interest
+struct Params
+{
+    dual a;
+    dual b;
+    dual c;
+};
+
+// The function that depends on parameters for which derivatives are needed
+dual f(dual x, const Params& params)
+{
+    return params.a * sin(x) + params.b * cos(x) + params.c * sin(x)*cos(x);
+}
+
+int main()
+{
+    Params params;   // initialize the parameter variables
+    params.a = 1.0;  // the parameter a of type dual, not double!
+    params.b = 2.0;  // the parameter b of type dual, not double!
+    params.c = 3.0;  // the parameter c of type dual, not double!
+
+    dual x = 0.5;  // the input variable x
+
+    dual u = f(x, params);  // the output variable u
+
+    double dudx = derivative(f, wrt(x), at(x, params));        // evaluate the derivative du/dx
+    double duda = derivative(f, wrt(params.a), at(x, params)); // evaluate the derivative du/da
+    double dudb = derivative(f, wrt(params.b), at(x, params)); // evaluate the derivative du/db
+    double dudc = derivative(f, wrt(params.c), at(x, params)); // evaluate the derivative du/dc
+
+    std::cout << "u = " << u << std::endl;         // print the evaluated output u
+    std::cout << "du/dx = " << dudx << std::endl;  // print the evaluated derivative du/dx
+    std::cout << "du/da = " << duda << std::endl;  // print the evaluated derivative du/da
+    std::cout << "du/db = " << dudb << std::endl;  // print the evaluated derivative du/db
+    std::cout << "du/dc = " << dudc << std::endl;  // print the evaluated derivative du/dc
+}
+
+/* NOTE: 
+
+This example would also work if real was used instead of dual. Should you
+need higher-order cross derivatives, however, e.g.,:
+
+    double d2udxda = derivative(f, wrt(x, params.a), at(x, params));
+
+then higher-order dual types are the right choicesince real types are
+optimally designed for higher-order directional derivatives.
+
+如果需要计算更高阶的交叉倒数（例如d2udxda）时，应该使用更高阶的dual类型，而不是real类型
+
+*/

test/forward/5multi_varialbe_function_rely_on_analytical_der.cpp

@@ -0,0 +1,95 @@
+// NOTE: 依赖于解析导数的多变量函数的导数
+
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/dual.hpp>
+using namespace autodiff;
+
+// Define functions A, Ax, Ay using double; analytical derivatives are available.
+double  A(double x, double y) { return x*y; }
+double Ax(double x, double y) { return x; }
+double Ay(double x, double y) { return y; }
+
+// Define functions B, Bx, By using double; analytical derivatives are available.
+double  B(double x, double y) { return x + y; }
+double Bx(double x, double y) { return 1.0; }
+double By(double x, double y) { return 1.0; }
+
+// Wrap A into Adual function so that it can be used within autodiff-enabled codes.
+// A是一个double相关的函数，无法与autodiff混合使用；这里将A包装成Adual
+dual Adual(dual const& x, dual const& y)
+{
+    dual res = A(x.val, y.val);
+
+    if(x.grad != 0.0)
+        res.grad += x.grad * Ax(x.val, y.val);
+
+    if(y.grad != 0.0)
+        res.grad += y.grad * Ay(x.val, y.val);
+
+    return res;
+}
+
+// Wrap B into Bdual function so that it can be used within autodiff-enabled codes.
+dual Bdual(dual const& x, dual const& y)
+{
+    dual res = B(x.val, y.val);
+
+    if(x.grad != 0.0)
+        res.grad += x.grad * Bx(x.val, y.val);
+
+    if(y.grad != 0.0)
+        res.grad += y.grad * By(x.val, y.val);
+
+    return res;
+}
+
+// Define autodiff-enabled C function that relies on Adual and Bdual
+dual C(dual const& x, dual const& y)
+{
+    const auto A = Adual(x, y);
+    const auto B = Bdual(x, y);
+    return A*A + B;
+}
+
+int main()
+{
+    // --- begin of test ---
+    double tmp=0.0;
+    dual t1 = tmp;
+    dual t2 = 0.0;
+    std::cout << t1 << " " << t2 << std::endl;
+    // --- end of test --- 
+
+    dual x = 1.0;
+    dual y = 2.0;
+
+    auto C0 = C(x, y);
+
+    // Compute derivatives of C with respect to x and y using autodiff!
+    auto Cx = derivative(C, wrt(x), at(x, y));
+    auto Cy = derivative(C, wrt(y), at(x, y));
+
+    // Compute expected analytical derivatives of C with respect to x and y
+    auto x0 = x.val;
+    auto y0 = y.val;
+    auto expectedCx = 2.0*A(x0, y0)*Ax(x0, y0) + Bx(x0, y0);
+    auto expectedCy = 2.0*A(x0, y0)*Ay(x0, y0) + By(x0, y0);
+
+    std::cout << "C0 = " << C0 << "\n";
+
+    std::cout << "Cx(computed) = " << Cx << "\n";
+    std::cout << "Cx(expected) = " << expectedCx << "\n";
+
+    std::cout << "Cy(computed) = " << Cy << "\n";
+    std::cout << "Cy(expected) = " << expectedCy << "\n";
+}
+
+// Output:
+// C0 = 7
+// Cx(computed) = 5
+// Cx(expected) = 5
+// Cy(computed) = 9
+// Cy(expected) = 9

test/forward/6gradient_vector_of_scalar_func.cpp

@@ -0,0 +1,29 @@
+// NOTE: 标量函数的梯度向量
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/real.hpp>
+#include <autodiff/forward/real/eigen.hpp>
+using namespace autodiff;
+
+// The scalar function for which the gradient is needed
+real f(const ArrayXreal& x)
+{
+    return (x * x.exp()).sum(); // sum([xi * exp(xi) for i = 1:5])
+}
+
+int main()
+{
+    using Eigen::VectorXd;
+
+    ArrayXreal x(5);                            // the input array x with 5 variables
+    x << 1, 2, 3, 4, 5;                         // x = [1, 2, 3, 4, 5]
+
+    real u;                                     // the output scalar u = f(x) evaluated together with gradient below
+
+    VectorXd g = gradient(f, wrt(x), at(x), u); // evaluate the function value u and its gradient vector g = du/dx
+
+    std::cout << "u = " << u << std::endl;      // print the evaluated output u
+    std::cout << "g = \n" << g << std::endl;    // print the evaluated gradient vector g = du/dx
+}

test/forward/7gradient_vector_of_scalar_func_with_params.cpp

@@ -0,0 +1,51 @@
+// NOTE: 带参数的标量函数的梯度向量 
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/real.hpp>
+#include <autodiff/forward/real/eigen.hpp>
+using namespace autodiff;
+
+// The scalar function for which the gradient is needed
+real f(const ArrayXreal& x, const ArrayXreal& p, const real& q)
+{
+    return (x * x).sum() * p.sum() * exp(q); // sum([xi * xi for i = 1:5]) * sum(p) * exp(q)
+}
+
+int main()
+{
+    using Eigen::VectorXd;
+
+    ArrayXreal x(5);    // the input vector x with 5 variables
+    x << 1, 2, 3, 4, 5; // x = [1, 2, 3, 4, 5]
+
+    ArrayXreal p(3);    // the input parameter vector p with 3 variables
+    p << 1, 2, 3;       // p = [1, 2, 3]
+
+    real q = -2;        // the input parameter q as a single variable
+
+    real u;             // the output scalar u = f(x, p, q) evaluated together with gradient below
+
+    VectorXd gx   = gradient(f, wrt(x), at(x, p, q), u);       // evaluate the function value u and its gradient vector gx = du/dx
+    VectorXd gp   = gradient(f, wrt(p), at(x, p, q), u);       // evaluate the function value u and its gradient vector gp = du/dp
+    VectorXd gq   = gradient(f, wrt(q), at(x, p, q), u);       // evaluate the function value u and its gradient vector gq = du/dq
+    VectorXd gqpx = gradient(f, wrt(q, p, x), at(x, p, q), u); // evaluate the function value u and its gradient vector gqpx = [du/dq, du/dp, du/dx]
+
+    std::cout << "u = " << u << std::endl;       // print the evaluated output u
+    std::cout << "gx = \n" << gx << std::endl;   // print the evaluated gradient vector gx = du/dx
+    std::cout << "gp = \n" << gp << std::endl;   // print the evaluated gradient vector gp = du/dp
+    std::cout << "gq = \n" << gq << std::endl;   // print the evaluated gradient vector gq = du/dq
+    std::cout << "gqpx = \n" << gqpx << std::endl; // print the evaluated gradient vector gqpx = [du/dq, du/dp, du/dx]
+}
+
+/*-------------------------------------------------------------------------------------------------
+=== Note ===
+---------------------------------------------------------------------------------------------------
+This example would also work if dual was used instead. However, if gradient,
+Jacobian, and directional derivatives are all you need, then real types are your
+best option. You want to use dual types when evaluating higher-order cross
+derivatives, which is not supported for real types.
+如果只需要计算梯度、雅可比矩阵和方向导数，使用real类型
+是最佳的。如果需要计算更高阶的交叉倒数，则应该使用dual类型。
+-------------------------------------------------------------------------------------------------*/

test/forward/8jacobian_matrix_of_vector_func.cpp

@@ -0,0 +1,32 @@
+// NOTE: 向量函数的雅可比矩阵 
+
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/real.hpp>
+#include <autodiff/forward/real/eigen.hpp>
+using namespace autodiff;
+
+// The vector function for which the Jacobian is needed
+VectorXreal f(const VectorXreal& x)
+{
+    return x * x.sum();
+}
+
+int main()
+{
+    using Eigen::MatrixXd;
+
+    VectorXreal x(5);                           // the input vector x with 5 variables
+    x << 1, 2, 3, 4, 5;                         // x = [1, 2, 3, 4, 5]
+
+    VectorXreal F;                              // the output vector F = f(x) evaluated together with Jacobian matrix below
+    MatrixXd J1 = jacobian(f, wrt(x), at(x), F); // evaluate the output vector F and the Jacobian matrix dF/dx
+
+    std::cout << "F = \n" << F << std::endl;    // print the evaluated output vector F
+    std::cout << "J1 = \n" << J1 << std::endl;    // print the evaluated Jacobian matrix dF/dx
+
+    MatrixXd J2 = jacobian(f, wrt(x), at(x));
+    std::cout << "J2 = \n" << J2 << std::endl;    // print the evaluated Jacobian matrix dF/dx
+}

test/forward/9jacobian_matrix_of_vector_func_with_params.cpp

@@ -0,0 +1,40 @@
+// NOTE: 带参数的向量函数的雅可比矩阵
+// C++ includes
+#include <iostream>
+
+// autodiff include
+#include <autodiff/forward/real.hpp>
+#include <autodiff/forward/real/eigen.hpp>
+using namespace autodiff;
+
+// The vector function with parameters for which the Jacobian is needed
+VectorXreal f(const VectorXreal& x, const VectorXreal& p, const real& q)
+{
+    return x * p.sum() * exp(q);
+}
+
+int main()
+{
+    using Eigen::MatrixXd;
+
+    VectorXreal x(5);    // the input vector x with 5 variables
+    x << 1, 2, 3, 4, 5;  // x = [1, 2, 3, 4, 5]
+
+    VectorXreal p(3);    // the input parameter vector p with 3 variables
+    p << 1, 2, 3;        // p = [1, 2, 3]
+
+    real q = -2;         // the input parameter q as a single variable
+
+    VectorXreal F;       // the output vector F = f(x, p, q) evaluated together with Jacobian below
+
+    MatrixXd Jx   = jacobian(f, wrt(x), at(x, p, q), F);       // evaluate the function and the Jacobian matrix Jx = dF/dx
+    MatrixXd Jp   = jacobian(f, wrt(p), at(x, p, q), F);       // evaluate the function and the Jacobian matrix Jp = dF/dp
+    MatrixXd Jq   = jacobian(f, wrt(q), at(x, p, q), F);       // evaluate the function and the Jacobian matrix Jq = dF/dq
+    MatrixXd Jqpx = jacobian(f, wrt(q, p, x), at(x, p, q), F); // evaluate the function and the Jacobian matrix Jqpx = [dF/dq, dF/dp, dF/dx]
+
+    std::cout << "F = \n" << F << std::endl;     // print the evaluated output vector F
+    std::cout << "Jx = \n" << Jx << std::endl;   // print the evaluated Jacobian matrix dF/dx
+    std::cout << "Jp = \n" << Jp << std::endl;   // print the evaluated Jacobian matrix dF/dp
+    std::cout << "Jq = \n" << Jq << std::endl;   // print the evaluated Jacobian matrix dF/dq
+    std::cout << "Jqpx = \n" << Jqpx << std::endl; // print the evaluated Jacobian matrix [dF/dq, dF/dp, dF/dx]
+}

test/forward/CMakeLists.txt

@@ -0,0 +1,9 @@
+
+file(GLOB CPPFILES "${CMAKE_CURRENT_SOURCE_DIR}/*.cpp")
+foreach(item ${CPPFILES})
+    get_filename_component(itemname ${item} NAME_WE)
+    add_executable(${itemname} ${item})
+    message(STATUS "add_executable(${itemname} ${item})")
+    # target_link_libraries(${itemname} Eigen3::Eigen)
+endforeach()
+