autodiff / autodiff Goto Github PK

automatic differentiation made easier for C++

License: MIT License

CMake 2.68% C++ 82.75% Starlark 0.26% Python 0.25% Cuda 14.06%

automatic-differentiation auto-differentiation autodiff autodifferentiation differentiation derivatives numerical-derivation

autodiff's Introduction

Overview

autodiff is a C++17 library that uses modern and advanced programming techniques to enable automatic computation of derivatives in an efficient, easy, and intuitive way.

Demonstration

Consider the following function f(x, y, z):

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

which we use use to evaluate the variable u = f(x, y, z):

double x = 1.0;
double y = 2.0;
double z = 3.0;
double u = f(x, y, z);

How can we minimally transform this code so that not only u, but also its derivatives ∂u/∂x, ∂u/∂y, and ∂u/∂z, can be computed?

The next two sections present how this can be achieved using two automatic differentiation algorithms implemented in autodiff: forward mode and reverse mode.

Forward mode

In a forward mode automatic differentiation algorithm, both output variables and one or more of their derivatives are computed together. For example, the function evaluation f(x, y, z) can be transformed in a way that it will not only produce the value of u, the output variable, but also one or more of its derivatives (∂u/∂x, ∂u/∂y, ∂u/∂z) with respect to the input variables (x, y, z).

Enabling forward automatic differentiation for the calculation of derivatives using autodiff is relatively simple. For our previous function f, we only need to replace the floating-point type double with autodiff::dual for both input and output variables:

dual f(const dual& x, const dual& y, const dual& z)
{
    return (x + y + z) * exp(x * y * z);
}

We can now compute the derivatives ∂u/∂x, ∂u/∂y, and ∂u/∂z as follows:

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));

The auxiliary function autodiff::wrt, an acronym for with respect to, is used to indicate which input variable (x, y, z) is the selected one to compute the partial derivative of f. The auxiliary function autodiff::at is used to indicate where (at which values of its parameters) the derivative of f is evaluated.

Reverse mode

In a reverse mode automatic differentiation algorithm, the output variable of a function is evaluated first. During this function evaluation, all mathematical operations between the input variables are "recorded" in an expression tree. By traversing this tree from top-level (output variable as the root node) to bottom-level (input variables as the leaf nodes), it is possible to compute the contribution of each branch on the derivatives of the output variable with respect to input variables.

Thus, a single pass in a reverse mode calculation computes all derivatives, in contrast with forward mode, which requires one pass for each input variable. Note, however, that it is possible to change the behavior of a forward pass so that many (perhaps even all) derivatives of an output variable are computed simultaneously (e.g., in a single forward pass, ∂u/∂x, ∂u/∂y, and ∂u/∂z are evaluated together with u, in contrast with three forward passes, each one computing the individual derivatives).

Similar as before, we can use autodiff to enable reverse automatic differentiation for our function f by simply replacing type double with autodiff::var as follows:

var f(var x, var y, var z)
{
    return (x + y + z) * exp(x * y * z);
}

The code below demonstrates how the derivatives ∂u/∂x, ∂u/∂y, and ∂u/∂z can be calculated:

var x = 1.0;
var y = 2.0;
var z = 3.0;
var u = f(x, y, z);

Derivatives dud = derivatives(u);

double dudx = dud(x);
double dudy = dud(y);
double dudz = dud(z);

The function autodiff::derivatives will traverse the expression tree stored in variable u and compute all its derivatives with respect to the input variables (x, y, z), which are then stored in the object dud. The derivative of u with respect to input variable x (i.e., ∂u/∂x) can then be extracted from dud using dud(x). The operations dud(x), dud(y), dud(z) involve no computations! Just extraction of derivatives previously computed with a call to function autodiff::derivatives.

Documentation

Check the documentation website for more details:

License

MIT License

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

autodiff's People

Contributors

Stargazers

Watchers

Forkers

antoniouf allanleal pehlivanian yochju keris sprinterzzj haitaozhao qiaohaijun zxhsama volpatto jyhyun00 ludkinm supersega backyes feixh pariterre cyy1991 yichencityu jinyaozhu nemoleiyang jyakaranda ibell ram-nad jianantian zoq aaditya3015 dclcs chrisdrozdowski stevemtse danilomendesdias igordzreyev gjaeger1 xyyeh dev-fennek daniellayeghi stephenberry yatli s-telescope xumaoxuan capri2014 hasantahir xtsxisaxns hendy3271 zhuxinyan forkbackups maze516 pettni flying-tiger simipro vermosen dmitrylyakh chengwei920412 miquelramirez andrewolal ogsdave mihirk284 hossen-code xrrain maksimpiriyev gayatri-a-b shejingbingxiansheng yamaha-bps yeongseok94 externalrepositories saikishor gokhanyu heronimus renatomatz mceowen vguerra mythlee gearslogy pipigenius mfkiwl varasinguluri jgornet rath3t wbthomason geo-lianliu yuanweizhang jan-grimo c-renton adam-ce jjaraalm artivis markisus darlanjurak baruuum perfcv org-mars mattarroz gregvernon louyu2015 ajunlonglive laplacekorea gatsby23 hs5530hs zhengjun-du dungducphan ddisqq

autodiff's Issues

Unable to nest functions

This one puzzles me, and I suspect it's a bug. If I expand the calls to pow2(x), the code compiles and works, but if I don't, I get the error below:

// C++ includes
#include <iostream>

// autodiff include
#include <autodiff/forward.hpp>
#include "Eigen/Dense"
#include <autodiff/forward/eigen.hpp>

template <typename T>
T pow2(const T x0) {
    return x0 * x0;
}

template <typename T>
T Rosenbrock(const T x0, const T x1) {
    return 100.0 * pow2(pow2(x0) - x1) + pow2(1 - x0);
}

int main()
{
    using namespace autodiff;
    Eigen::VectorXdual x(2);  // the input vector x
    x << -0.5, 0.5;
    auto Rosenbrockdual = [](const Eigen::VectorXdual& x) {
        return Rosenbrock(x[0], x[1]);
    };
    dual F;
    auto grad = gradient(Rosenbrockdual, wrt(x), at(x), F);
    std::cout << grad;
}

yields

1>------ Build started: Project: bugdualvec, Configuration: Debug x64 ------
1>bugdualvec.cpp
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:  'return': cannot convert from 'autodiff::forward::BinaryExpr<autodiff::forward::MulOp,L,R>' to 'autodiff::forward::BinaryExpr<autodiff::forward::AddOp,L,R>'
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:         with
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:         [
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:             L=const autodiff::forward::BinaryExpr<autodiff::forward::AddOp,autodiff::forward::Dual<double,double>,autodiff::forward::UnaryExpr<autodiff::forward::NegOp,const autodiff::forward::Dual<double,double> &>> &,
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:             R=const autodiff::forward::BinaryExpr<autodiff::forward::AddOp,autodiff::forward::Dual<double,double>,autodiff::forward::UnaryExpr<autodiff::forward::NegOp,const autodiff::forward::Dual<double,double> &>> &
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:         ]
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:         and
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:         [
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:             L=autodiff::forward::Dual<double,double>,
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:             R=autodiff::forward::UnaryExpr<autodiff::forward::NegOp,const autodiff::forward::Dual<double,double> &>
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:         ]
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): message :  No constructor could take the source type, or constructor overload resolution was ambiguous
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :  see reference to function template instantiation 'T pow2<autodiff::forward::BinaryExpr<autodiff::forward::AddOp,L,R>>(const T)' being compiled
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :         with
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :         [
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :             T=autodiff::forward::BinaryExpr<autodiff::forward::AddOp,autodiff::forward::Dual<double,double>,autodiff::forward::UnaryExpr<autodiff::forward::NegOp,const autodiff::forward::Dual<double,double> &>>,
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :             L=autodiff::forward::Dual<double,double>,
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :             R=autodiff::forward::UnaryExpr<autodiff::forward::NegOp,const autodiff::forward::Dual<double,double> &>
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :         ]
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(26): message :  see reference to function template instantiation 'T Rosenbrock<autodiff::forward::Dual<double,double>>(const T,const T)' being compiled
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(26): message :         with
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(26): message :         [
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(26): message :             T=autodiff::forward::Dual<double,double>
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(26): message :         ]
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:  'return': cannot convert from 'autodiff::forward::BinaryExpr<autodiff::forward::MulOp,L,R>' to 'autodiff::forward::BinaryExpr<autodiff::forward::AddOp,L,R>'
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:         with
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:         [
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:             L=const autodiff::forward::BinaryExpr<autodiff::forward::AddOp,int,autodiff::forward::UnaryExpr<autodiff::forward::NegOp,const autodiff::forward::Dual<double,double> &>> &,
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:             R=const autodiff::forward::BinaryExpr<autodiff::forward::AddOp,int,autodiff::forward::UnaryExpr<autodiff::forward::NegOp,const autodiff::forward::Dual<double,double> &>> &
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:         ]
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:         and
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:         [
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:             L=int,
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:             R=autodiff::forward::UnaryExpr<autodiff::forward::NegOp,const autodiff::forward::Dual<double,double> &>
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): error C2440:         ]
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(12,1): message :  No constructor could take the source type, or constructor overload resolution was ambiguous
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :  see reference to function template instantiation 'T pow2<autodiff::forward::BinaryExpr<autodiff::forward::AddOp,L,R>>(const T)' being compiled
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :         with
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :         [
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :             T=autodiff::forward::BinaryExpr<autodiff::forward::AddOp,int,autodiff::forward::UnaryExpr<autodiff::forward::NegOp,const autodiff::forward::Dual<double,double> &>>,
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :             L=int,
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :             R=autodiff::forward::UnaryExpr<autodiff::forward::NegOp,const autodiff::forward::Dual<double,double> &>
1>D:\Code\helmderiv\snippets\bugdualvec.cpp(17): message :         ]
1>Done building project "bugdualvec.vcxproj" -- FAILED.

Determine function dependency once, and speed up future gradient calculations

Given a function f of n variables, implement a function dependency such that:

auto ideps = dependency(f, n);

returns ideps, a vector with the indices of all variables x_i for which ∂f/∂n_i is not identical to zero.

Create then overloads of grad functions that accept this indices to speed up gradient evaluations (by skiping those variables x_j for which f is independent of -- ∂f/∂n_j ≡ 0!):

auto dfdx = grad(f, x, ideps);

Consider that:

if f is a scalar function, ideps is a vector.
if f is a vector function, ideps is a matrix (or vector of vectors).

This feature will possibly require some conditions about the function, such as they should be continuous on the variables.

Potential difficulties: what to do in case of conditional branches (e.g., if/else)?

This feature might require function f to be templated on the numerical type, so that a special type is used to evaluate f and obtain some introspection about its variable dependency.

Forward mode: Higher-order derivatives of a multi-variable function

Hi Allan:

Thank you very much for making your library available, and in particular the thought you have invested into making the API so intuittive.

I plan to use your library extensively to compute sensitivities of complicated financial instruments.

At the moment, I am unit testing, with the Catch2 C++ unit test framework, the examples in your tutorial section, and testing the output from autodiff against the answer computed by symbolic differentiation. My development environment is VS2017 (Professional).

Everything looks fine, until I reach the example "Higher-order derivatives of a multi-variable function" in the forward mode. The answer I obtain from autodiff for ux (first derivative of u with respect to x) is 3.192695, but the correct answer should be 2.25; I have not proceeded beyond this point, as I wish to resolve the issue.

Interestingly, I progressively added terms into the definition of the function f, and found that the term log(x/y) seems to be a problem - if I omit this term, then ux is fine.

In the attached zip folder you will find my unit test.

Appreciate any assitance you could provide to resolve my issue.

autodiff.zip

Really nice documentation!

I know I have opened rather a lot of issues in the last 24 hours (and have a few more to go), but I just wanted to comment that I think you have some of the nicest documentation I have ever read for any C++ project. It's well put together and builds on itself very nicely so it is easy(ish) to see how to build up from simple hello world examples to something a bit more involved.

Speeding up Compile Time

On my laptop it takes about 3 seconds to compile the reverse single variable example. Any ideas how to speed this up?

INSTALL and CONTRIBUTING files

I think that it would be nice if autodiff provides INSTALL and CONTRIBUTING files with instructions for those that want to use and help developing autodiff.

Reverse mode run-time

Hi,

I have an optimization problem where my cost function f depends on n (usually >100) input variables (vector FA) and also a handful of vector parameters with the same size:

var f(const VectorXd& TR, const VectorXd& TE, VectorXvar FA, const VectorXd& RF)
{  
       //cost function calculation;
}

void cpp17_EPG_GRE_autodiff (const int& n, const VectorXd TR0, VectorXd TE0, VectorXd FA0, VectorXd RF0, double& F, VectorXd& grad)
{
        VectorXvar FA(n);
        for(int i=0; i<n; i++){
            FA(i) = FA0(i);
        }
        
     	var fval = f(TR0, TE0, FA, RF0);
        std::cout << "I'm here 1 " << val(fval) << std::endl
       
	VectorXd dydx = gradient(fval, FA);
        std::cout << "I'm here 2 " << val(fval) << std::endl;

        F = val(fval);

        for(int i=0; i<n; i++){
            grad(i) = dydx(i);
        }
}

For testing, I set n=6 variables and performed a simple function evaluation which takes <0.001 seconds, while calculating the cost function gradient takes 15 seconds. With the std::cout commands I can see that it is taking a lot of time at VectorXd dydx = gradient(fval, FA);. Is there anything wrong in my implementation from what I have shown?

Thanks for your help,
David

Derivatives w.r.t. `var` and references

Can you tell me if the following code makes sense?

var w = energy(m, D);
VectorXvar u(2*m.nn);
        
int i; for (i=0; i<m.nn; i++)
{
    u(2*i) = m.one[i];
    u(2*i+1) = m.two[i];
}
        
VectorXd dwdu = gradient(w,u);

So I'm building up a VectorXvar from two lists I have (my degrees of freedom) and then differentiation the energy with this new list. Will this work? I'm getting all zeros and I wonder if it's because I'm not using the variables directly but rather building a new list....

Multi Vector-variable Jacobian

Previous discussion:

Regarding the multi variable jacobian, I don't know how to add a new item to your tutorial.
https://autodiff.github.io/tutorials/#jacobian-of-a-vector-function
Specifically, I'd like to add a section after:

Jacobian of a vector function

And call it

Jacobian of a multi-vector function

Then attached with a section of code. Could you point me to which file I should change? I can then go ahead change it and do a PR.

As to the tutorial, I propose we discuss this in a new issue. Would you please open a new issue and explain (briefly) about the tutorial (with explanation about this multi-vector function)?

The tutorial would essentially be just an example of how to compute a Jacobian of a multi vector-variable function, with the help of lambda.

VectorXdual g(const VectorXdual& x, const VectorXdual& v)
{
  Vector3dual vec_p;        vec_p << x(0), x(1), x(2);
  Vector3dual vec_o;        vec_o << x(3), x(4), x(5);
  Vector3dual vec_pv;       vec_pv<< x(6), x(7), x(8);
  Vector3dual measured_vel = R(vec_o).transpose() * vec_pv;
  VectorXdual g_value(9);   g_value << vec_p, vec_o, measured_vel;
  g_value = g_value + v;
  return g_value;
}

    VectorXdual x = eMU_t;
    VectorXdual v = VectorXd::Zero(Zt.size());
    auto gCt = [=](const VectorXdual& x) {return g(x, v); };
    auto gWt = [=](const VectorXdual& v) {return g(x, v); };
    VectorXdual g_value = g(x,v);

    MatrixXd Ct = autodiff::jacobian(gCt,g_value,x);
    MatrixXd Wt = autodiff::jacobian(gWt,g_value,v);

This is part of the code for implementing EKF. Essentially g(x,v) is the function, where both x and v are vector, and I want their corresponding jacobian.

In the tutorial I'll clean up the code and make a easier one with complete variable definitions and stuff.

Also, what is the type of gCt and gWt? auto is the only way I can get it to compile.

Higher-order derivatives using forward mode?

Hi thanks for this work.

I was just wondering is it possible to do higher-order derivatives using forward mode?

Problems Installing and linking with Eigen

None of the CMakeList files contain a command to find the Eigen package so the final installation command fails unless the files are altered. The example files that use the Eigen directories have an issue in the include command #include <eigen3/Eigen/Core> should just be #include <Eigen/Core>.

Products of MatrixXd and MatrixXdual not working

In my code the products of type double and type dual inside Eigen Vectors and Matrices seems to be not working. Do i overlook something? My MWE looks like the following

#include <iostream>
#include <Eigen/Core>
#include <autodiff/forward.hpp>
#include <autodiff/forward/eigen.hpp>

using namespace Eigen;
using namespace autodiff;
using namespace std;

dual testEnergy(const Vector3dual& DDual) {

    MatrixX3dual   C(2,3);
    Matrix2Xd A(2,4);
    MatrixX3dual B(4 ,3);
    A << 1, 2, 3, 4, 5, 6, 7, 8;

    B.setIdentity(4,3);
    B = DDual.norm()*B;
    C = A*B; //Not working
//    C = A.cast<dual>() * B; //Working but useless temporary constructed of type dual type?
return C.norm();
}

int main() {
    Vector3dual DDual;
    DDual<< 1,5,13;
    dual dualenergy;
    VectorXd g = gradient(testEnergy, wrt(DDual), at(DDual), dualenergy);
    cout<<g<<endl;
    return 0;
}

which produces the following compilation error:

Scanning dependencies of target Test_AutoDiff
[ 50%] Building CXX object CMakeFiles/Test_AutoDiff.dir/main.cpp.obj
In file included from PATH/src/eigen3/Eigen/Core:313:0,
                 from C:\Users\Alex\CLionProjects\Test_AutoDiff\main.cpp:2:
PATH/src/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h: In instantiation of 'void Eigen::internal::gebp_traits<_LhsScalar, _RhsScalar, _ConjLhs, _ConjRhs, Arch, _PacketSize>::madd(const LhsPacketType&, const RhsPacketType&, AccPacketType&, RhsPacketType&, const LaneIdType&) const [with LhsPacketType = autodiff::forward::Dual<double, double>; RhsPacketType = double; AccPacketType = autodiff::forward::Dual<double, double>; LaneIdType = Eigen::internal::FixedInt<0>; _LhsScalar = autodiff::forward::Dual<double, double>; _RhsScalar = double; bool _ConjLhs = false; bool _ConjRhs = false; int Arch = 0; int _PacketSize = 0]':
PATH/src/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:2094:15:   required from 'void Eigen::internal::gebp_kernel<LhsScalar, RhsScalar, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>::operator()(const DataMapper&, const LhsScalar*, const RhsScalar*, Index, Index, Index, Eigen::internal::gebp_kernel<LhsScalar, RhsScalar, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>::ResScalar, Index, Index, Index, Index) [with LhsScalar = double; RhsScalar = autodiff::forward::Dual<double, double>; Index = int; DataMapper = Eigen::internal::blas_data_mapper<autodiff::forward::Dual<double, double>, int, 0, 0, 1>; int mr = 1; int nr = 4; bool ConjugateLhs = false; bool ConjugateRhs = false; Eigen::internal::gebp_kernel<LhsScalar, RhsScalar, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs>::ResScalar = autodiff::forward::Dual<double, double>]'
PATH/src/eigen3/Eigen/src/Core/products/GeneralMatrixMatrix.h:198:15:   required from 'static void Eigen::internal::general_matrix_matrix_product<Index, LhsScalar, LhsStorageOrder, ConjugateLhs, RhsScalar, RhsStorageOrder, ConjugateRhs, 0, ResInnerStride>::run(Index, Index, Index, const LhsScalar*, Index, const RhsScalar*, Index, Eigen::internal::general_matrix_matrix_product<Index, LhsScalar, LhsStorageOrder, ConjugateLhs, RhsScalar, RhsStorageOrder, ConjugateRhs, 0, ResInnerStride>::ResScalar*, Index, Index, Eigen::internal::general_matrix_matrix_product<Index, LhsScalar, LhsStorageOrder, ConjugateLhs, RhsScalar, RhsStorageOrder, ConjugateRhs, 0, ResInnerStride>::ResScalar, Eigen::internal::level3_blocking<LhsScalar, RhsScalar>&, Eigen::internal::GemmParallelInfo<Index>*) [with Index = int; LhsScalar = double; int LhsStorageOrder = 0; bool ConjugateLhs = false; RhsScalar = autodiff::forward::Dual<double, double>; int RhsStorageOrder = 0; bool ConjugateRhs = false; int ResInnerStride = 1; Eigen::internal::general_matrix_matrix_product<Index, LhsScalar, LhsStorageOrder, ConjugateLhs, RhsScalar, RhsStorageOrder, ConjugateRhs, 0, ResInnerStride>::ResScalar = autodiff::forward::Dual<double, double>]'
PATH/src/eigen3/Eigen/src/Core/products/GeneralMatrixMatrix.h:230:14:   required from 'void Eigen::internal::gemm_functor<Scalar, Index, Gemm, Lhs, Rhs, Dest, BlockingType>::operator()(Index, Index, Index, Index, Eigen::internal::GemmParallelInfo<Index>*) const [with Scalar = autodiff::forward::Dual<double, double>; Index = int; Gemm = Eigen::internal::general_matrix_matrix_product<int, double, 0, false, autodiff::forward::Dual<double, double>, 0, false, 0, 1>; Lhs = Eigen::Matrix<double, 2, -1>; Rhs = Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>; Dest = Eigen::Matrix<autodiff::forward::Dual<double, double>, 2, 3, 0, 2, 3>; BlockingType = Eigen::internal::gemm_blocking_space<0, double, autodiff::forward::Dual<double, double>, 2, 3, -1, 1, false>]'
PATH/src/eigen3/Eigen/src/Core/products/Parallelizer.h:114:7:   required from 'void Eigen::internal::parallelize_gemm(const Functor&, Index, Index, Index, bool) [with bool Condition = false; Functor = Eigen::internal::gemm_functor<autodiff::forward::Dual<double, double>, int, Eigen::internal::general_matrix_matrix_product<int, double, 0, false, autodiff::forward::Dual<double, double>, 0, false, 0, 1>, Eigen::Matrix<double, 2, -1>, Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>, Eigen::Matrix<autodiff::forward::Dual<double, double>, 2, 3, 0, 2, 3>, Eigen::internal::gemm_blocking_space<0, double, autodiff::forward::Dual<double, double>, 2, 3, -1, 1, false> >; Index = int]'
PATH/src/eigen3/Eigen/src/Core/products/GeneralMatrixMatrix.h:509:9:   required from 'static void Eigen::internal::generic_product_impl<Lhs, Rhs, Eigen::DenseShape, Eigen::DenseShape, 8>::scaleAndAddTo(Dest&, const Lhs&, const Rhs&, const Scalar&) [with Dest = Eigen::Matrix<autodiff::forward::Dual<double, double>, 2, 3, 0, 2, 3>; Lhs = Eigen::Matrix<double, 2, -1>; Rhs = Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>; Eigen::internal::generic_product_impl<Lhs, Rhs, Eigen::DenseShape, Eigen::DenseShape, 8>::Scalar = autodiff::forward::Dual<double, double>]'
PATH/src/eigen3/Eigen/src/Core/products/GeneralMatrixMatrix.h:445:20:   [ skipping 5 instantiation contexts, use -ftemplate-backtrace-limit=0 to disable ]
PATH/src/eigen3/Eigen/src/Core/Matrix.h:332:31:   required from 'Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::Matrix(const T&) [with T = Eigen::Product<Eigen::Matrix<double, 2, -1>, Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>, 0>; _Scalar = autodiff::forward::Dual<double, double>; int _Rows = 2; int _Cols = 3; int _Options = 0; int _MaxRows = 2; int _MaxCols = 3]'
PATH/src/eigen3/Eigen/src/Core/AssignEvaluator.h:822:41:   required from 'void Eigen::internal::call_assignment(Dst&, const Src&, const Func&, typename Eigen::internal::enable_if<Eigen::internal::evaluator_assume_aliasing<Src>::value, void*>::type) [with Dst = Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>; Src = Eigen::Product<Eigen::Matrix<double, 2, -1>, Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>, 0>; Func = Eigen::internal::assign_op<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >; typename Eigen::internal::enable_if<Eigen::internal::evaluator_assume_aliasing<Src>::value, void*>::type = void*]'
PATH/src/eigen3/Eigen/src/Core/AssignEvaluator.h:808:18:   required from 'void Eigen::internal::call_assignment(Dst&, const Src&) [with Dst = Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>; Src = Eigen::Product<Eigen::Matrix<double, 2, -1>, Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>, 0>]'
PATH/src/eigen3/Eigen/src/Core/PlainObjectBase.h:779:32:   required from 'Derived& Eigen::PlainObjectBase<Derived>::_set(const Eigen::DenseBase<OtherDerived>&) [with OtherDerived = Eigen::Product<Eigen::Matrix<double, 2, -1>, Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>, 0>; Derived = Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>]'
PATH/src/eigen3/Eigen/src/Core/Matrix.h:225:24:   required from 'Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::operator=(const Eigen::DenseBase<OtherDerived>&) [with OtherDerived = Eigen::Product<Eigen::Matrix<double, 2, -1>, Eigen::Matrix<autodiff::forward::Dual<double, double>, -1, 3, 0, -1, 3>, 0>; _Scalar = autodiff::forward::Dual<double, double>; int _Rows = -1; int _Cols = 3; int _Options = 0; int _MaxRows = -1; int _MaxCols = 3]'
C:\Users\Alex\CLionProjects\Test_AutoDiff\main.cpp:20:11:   required from here
PATH/src/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:527:18: error: cannot convert 'Eigen::internal::conj_helper<autodiff::forward::Dual<double, double>, double, false, false>::Scalar {aka autodiff::forward::Dual<double, double>}' to 'double' in assignment
     tmp = b; tmp = cj.pmul(a,tmp); c = padd(c,tmp);
PATH/src/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:527:44: error: no matching function for call to 'padd(autodiff::forward::Dual<double, double>&, double&)'
     tmp = b; tmp = cj.pmul(a,tmp); c = padd(c,tmp);
                                        ~~~~^~~~~~~
In file included from PATH/src/eigen3/Eigen/Core:160:0,
                 from C:\Users\Alex\CLionProjects\Test_AutoDiff\main.cpp:2:
PATH/src/eigen3/Eigen/src/Core/GenericPacketMath.h:163:1: note: candidate: template<class Packet> Packet Eigen::internal::padd(const Packet&, const Packet&)
 padd(const Packet& a, const Packet& b) { return a+b; }
 ^~~~
PATH/src/eigen3/Eigen/src/Core/GenericPacketMath.h:163:1: note:   template argument deduction/substitution failed:
In file included from PATH/src/eigen3/Eigen/Core:313:0,
                 from C:\Users\Alex\CLionProjects\Test_AutoDiff\main.cpp:2:
PATH/src/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:527:44: note:   deduced conflicting types for parameter 'const Packet' ('autodiff::forward::Dual<double, double>' and 'double')
     tmp = b; tmp = cj.pmul(a,tmp); c = padd(c,tmp);
                                        ~~~~^~~~~~~
CMakeFiles\Test_AutoDiff.dir\build.make:62: recipe for target 'CMakeFiles/Test_AutoDiff.dir/main.cpp.obj' failed
CMakeFiles\Makefile2:74: recipe for target 'CMakeFiles/Test_AutoDiff.dir/all' failed

N-th derivative

Working through the examples in the tutorial (very nice, by the way), I can see how to hard-code the construction of the wrt argument for first- and second-order derivatives of a one-dimensional (and multivariate) function, and that works fine. But when it comes to higher derivatives, I've had no luck. I've tried every combination and permutation of std::apply, perfect forwarding, and so on that I can think of, but I have to think that this might just be a current limitation of how argument forwarding is handled in wrt()? If not, I would be very happy to be wrong. I feel like I'm very, very close to getting this working. Of course I can hard-code the higher derivatives (as I will do for now), but I don't want to have to do that.

// C++ includes
#include <iostream>
#include <chrono>
#include <iomanip>
using namespace std;

// autodiff include
#include <autodiff/forward.hpp>
using namespace autodiff;

/*
Calculate the n-th derivative of the (one-dimensional) function sin(x)
*/
template <std::size_t n> double deriv(double x)
{
    using vartype = HigherOrderDual<n>;
    vartype _x = x; // explicit copy construction (makes the shared_ptr)
    std::array<vartype, n> args; args.fill(_x); // fill copy constructor
    return derivative(
        [](vartype x)->vartype{ return sin(x); }, 
        std::apply(wrt, args), 
        at(_x)
    );
}

int main()
{
    std::cout << deriv<1>(0) << std::endl;
}

yields

[ 50%] Building CXX object CMakeFiles/Nderiv.dir/snippets/Nderiv.cpp.o
/home/ian/Code/helmderiv/snippets/Nderiv.cpp: In instantiation of ‘double deriv(double) [with long unsigned int n = 1]’:
/home/ian/Code/helmderiv/snippets/Nderiv.cpp:28:28:   required from here
/home/ian/Code/helmderiv/snippets/Nderiv.cpp:21:19: error: no matching function for call to ‘apply(<unresolved overloaded function type>, std::array<autodiff::forward::Dual<double, double>, 1>&)’
         std::apply(wrt, args),
         ~~~~~~~~~~^~~~~~~~~~~
In file included from /usr/include/c++/7/bits/unique_ptr.h:37:0,
                 from /usr/include/c++/7/bits/locale_conv.h:41,
                 from /usr/include/c++/7/locale:43,
                 from /usr/include/c++/7/iomanip:43,
                 from /home/ian/Code/helmderiv/snippets/Nderiv.cpp:4:
/usr/include/c++/7/tuple:1668:5: note: candidate: template<class _Fn, class _Tuple> constexpr decltype(auto) std::apply(_Fn&&, _Tuple&&)
     apply(_Fn&& __f, _Tuple&& __t)
     ^~~~~
/usr/include/c++/7/tuple:1668:5: note:   template argument deduction/substitution failed:
/home/ian/Code/helmderiv/snippets/Nderiv.cpp:21:19: note:   couldn't deduce template parameter ‘_Fn’
         std::apply(wrt, args),
         ~~~~~~~~~~^~~~~~~~~~~
/home/ian/Code/helmderiv/snippets/Nderiv.cpp:22:11: error: use of ‘auto autodiff::forward::at(Args&& ...) [with Args = {autodiff::forward::Dual<double, double>&}]’ before deduction of ‘auto’
         at(_x)
         ~~^~~~
CMakeFiles/Nderiv.dir/build.make:62: recipe for target 'CMakeFiles/Nderiv.dir/snippets/Nderiv.cpp.o' failed
make[3]: *** [CMakeFiles/Nderiv.dir/snippets/Nderiv.cpp.o] Error 1
CMakeFiles/Makefile2:141: recipe for target 'CMakeFiles/Nderiv.dir/all' failed
make[2]: *** [CMakeFiles/Nderiv.dir/all] Error 2
CMakeFiles/Makefile2:153: recipe for target 'CMakeFiles/Nderiv.dir/rule' failed
make[1]: *** [CMakeFiles/Nderiv.dir/rule] Error 2
Makefile:144: recipe for target 'Nderiv' failed
make: *** [Nderiv] Error 2

How to convert a MatrixXdual to MatrixXd?

I'm expecting something like:

MatrixXdual mat_dual;
MatrixXd mat_d = val(mat_dual);

However apparently val() can only be applied to dual, rather than any container of dual. Is there an existing way to do that?

Also, it took me quite a while to figure out how to do jacobian on multi-vector-variable function and conversion between dual and double, the readme should be updated. I can help if you want.

Implement class TaylorSeries

This is a proposal for the implementation of a class TaylorSeries and a few other convenient methods to deal with Taylor series of multi-variable functions along a direction vector.

Let f denote a multi-variable function whose input vector x and output vector y have dimensions n and m respectively, with n, m > 1. Let v denote a direction vector and xo be a reference input vector, both given. We are interested in the Taylor series of the single-variable function g(t) = f(xo + tv):

g(t) = a0 + a1*t + a2*t**2 + a3*t**3 + ...

Some usage example (say, the Taylor series up to order 4) of function f:

auto g = taylor(f, at(x0), along(v), order(4));

Alternatively:

TaylorSeries g(f, x0, v, 4); // or
TaylorSeries g(f, at(x0), along(v), order(4));

The following members in class TaylorSeries make sense:

TaylorSeries::direction() returns the direction vector v along which the Taylor series were generated.
TaylorSeries::order() returns the truncation order of the Taylor series.
TaylorSeries::at() returns the reference vector x0 at which the Taylor series were generated.
TaylorSeries::coefficients() returns a matrix whose columns are the coefficient vectors a0, ..., am.
TaylorSeries::derivatives() returns a matrix whose columns are the directional derivatives at order 0, 1, ..., order
TaylorSeries::eval(t) evaluates the Taylor series at parameter t.

Implement gradient function for forward mode so that other arguments can be provided

Inspired from discussion in issue #20, there is a need to have the gradient function (in forward.hpp) defined so that the following usage becomes possible:

VectorXd dudx = gradient(f, wrt(x), x, a, b);

Currently, it is assumed that f is defined with a single argument of type VectorXdual. Thus, we currently need to wrap function f in a lambda function to use gradient:

auto fab = [=](const VectorXdual& x) { return f(x, a, b); };
VectorXd dudx = gradient(fab, x);

Location of VectorXdual

In my opinion, the VectorXdual typedef doesn't belong in the Eigen namespace, I think it would be more sensible to put it in the autodiff namespace. It's hard to see this from the examples because the using namespace ... paradigm is used which obfuscates where things are coming from.

Side effects

Are there side effects in autodiff? I ask because when I put my solver into a loop I get strange results, even though I use fresh variables each time.

int n; for (n=1; n<20; n*=2)
{
    Mesh m = get_mesh(n);
    VectorXd U(nn), V(nn);

    solve(mesh,U,V);
}

I'm not sure how the variable tracking is done. Is there a way to reset autodiff and say "I'm solving a new problem now" ?

Derivatives of General Functions

I think I might have made a mistake in my understanding of autodiff. Taking the reverse single variable example as a starting point:

var f(var x)
{
    return 1 + x + x*x + 1/x + log(x);
}

int main()
{
    var x = 2.0;                         // the input variable x
    var u = f(x);                        // the output variable u

    Derivatives dud = derivatives(u);    // evaluate all derivatives of u

    var dudx = dud(x);                   // extract the derivative du/dx

    cout << "u = " << u << endl;         // print the evaluated output u
    cout << "du/dx = " << dudx << endl;  // print the evaluated derivative du/dx
}

This outputs

u = 8.19315
du/dx = 5.25

Now look what happens when I replace this with the following:

var get_x()
{
    var x = 0.0;
    int i; for (i=0; i<2; i++)
        x += 1;
        
    return x;
}

and then

    var x = get_x();                         // the input variable x
    var u = f(x);                        // the output variable u

So I'm getting x with a function which is a loop. This is now the output:

u = 8.19315
du/dx = 0

What happened? Even though x is still equal 2.0 I can't get a derivative. How did my var x suddenly not be usable for the target derivative variable?

Optimize reverse mode

I did fast look too the code. And my first idea to add allocator, but I'm not sure. @allanleal If you have any ideas, please tell =)

Can't extract a value as double from an VectorXdual

I'm trying to calculate the jacobian of a given vector function - for that I kept pretty close to the tutorial Jacobian of a vector function. Everything seems to work fine except for the conversion of values extracted from a VectorXdual - I've tried access via (), [] or coeff() I've received dual or scalar objects from these methods but never double values. I've tried casting the duals or scalars to double but unfortunately that did not work for me - What is the ideal process of extraction/conversion here?

Compiling errors

I got a ton of errors when trying to compile the Jacobian example

Karls-Mac-mini:scrap karl$ g++ jacobian.cpp -o jacobian.o && ./jacobian.o
In file included from jacobian.cpp:10:
In file included from /usr/local/include/autodiff/forward.hpp:33:
/usr/local/include/autodiff/forward/forward.hpp:213:17: error: unknown type name 'constexpr'
struct isExpr { constexpr static bool value = false; };
^
/usr/local/include/autodiff/forward/forward.hpp:213:27: error: expected member name or ';' after
declaration specifiers
struct isExpr { constexpr static bool value = false; };

(This is just the first few lines)

Am I using the wrong compiler?

Cast to double fails in reverse mode

In the examples, in reverse mode, piping derivatives into std::cout works. But I am unable to automatically cast them to double precision, which should be a supported cast I think. Printing the value (via std::cout) without the double cast does work.

// C++ includes
#include <iostream>
using namespace std;

// autodiff include
#include <autodiff/reverse.hpp>
using namespace autodiff;

// The function that depends on parameters for which derivatives are needed
template <typename T>
T alpha(T tau, T delta)
{
    return tau + delta*exp(-tau);
}

int main()
{
    using vartype = var;
    vartype tau = 0.5, delta = 0.7;  // the input variables
    vartype a = alpha<vartype>(tau, delta);  // the output variable u
    DerivativesX dud = derivativesx(a);
    double dudtau = dud(tau);
    std::cout << dudtau << std::endl;
}

with

Scanning dependencies of target double_cast
[ 16%] Building CXX object CMakeFiles/double_cast.dir/double_cast.cpp.o
/home/ian/Code/helmderiv/double_cast.cpp: In function ‘int main()’:
/home/ian/Code/helmderiv/double_cast.cpp:22:28: error: cannot convert ‘autodiff::var’ to ‘double’ in initialization
     double dudtau = dud(tau);
                            ^
CMakeFiles/double_cast.dir/build.make:62: recipe for target 'CMakeFiles/double_cast.dir/double_cast.cpp.o' failed
make[2]: *** [CMakeFiles/double_cast.dir/double_cast.cpp.o] Error 1
CMakeFiles/Makefile2:67: recipe for target 'CMakeFiles/double_cast.dir/all' failed
make[1]: *** [CMakeFiles/double_cast.dir/all] Error 2
Makefile:83: recipe for target 'all' failed
make: *** [all] Error 2

Implement output operator '<<' for dual expressions

We need to support output of dual expressions, for example:

dual a = 10.0;
std::cout << sin(a) << std::endl;

Error on computing the Jacobian of a vector-valued function when output and input vectors have different dimensions

I have a vector valued function that returns a 3 vector but I get an error when I take the gradient with respect to the input vector which has 12 variables. Can autodiff do a non-square jacobian of a vector valued function?

VectorXdual f(const VectorXdual& x)
{
    VectorXdual res = VectorXdual::Zero(3);
    dual theta = 0.0;

    for(int i(0); i < x.rows()/3; i++)
    {
        Eigen::Matrix3dual R;
        R << cos(theta), -sin(theta), 0, sin(theta), cos(theta), 0, 0, 0, 1;
        res.head(3) += R * x.segment(3*i, 3);
        res(2) = wrap(res(2));
        theta = res(2);
    }

    return res;
}

int main()
{
    //Vector function. Doesn't like it not being square
    VectorXdual u = f(x);
    Eigen::MatrixXd J = autodiff::jacobian(f, u, x);

    std::cout << "dux_dx = \n" << dux_dx << "\n\n";
    std::cout << "duy_dx = \n" << duy_dx << "\n\n";
    std::cout << "duphi_dx = \n" << duphi_dx << "\n\n";

    std::cout << "J = \n" << J << "\n\n";

    return 0;
}

Cannot compile with g++ 6.3.0 or clang 3.8.1

Hello,

I am looking to use autodiff in my project.
Standard cmake install fails on my device.

mkdir build; cd build; cmake ../autodiff; make install

output.txt

I have tried using with g++ 6.3.0 or clang 3.8.1
I have tried adding -std=c++17 to CMAKE_CXX_FLAGS but it does not improve (and C++17 is already requested in CMakeLists.txt).

By the way, it would be nice if a cmake option was added to disable building examples.

Overload method fabs and abs

Overload methods fabs and abs so that abs(a - 0.1) works without doing abs(val(a) - 0.1).

This comes from issue #13 .

Multiple derivatives in forward mode in one pass

In the documentation for reverse mode is written: Note, however, that it is possible to change the behavior of a forward pass so that many (even all) derivatives of an output variable are computed simultaneously (e.g., in a single forward pass, ∂u/∂x, ∂u/∂y, and ∂u/∂z are evaluated together with u, in contrast with three forward passes, each one computing the individual derivatives).

Do you have an example of this somewhere?

Can't add to VectorXdual

On MSVC 2019, this code

// C++ includes
#include <iostream>

// autodiff include
#include <autodiff/forward.hpp>
#include "Eigen/Dense"
#include <autodiff/forward/eigen.hpp>
using namespace autodiff;

int main()
{
    Eigen::VectorXdual x(2);  // the input vector x
    x << -0.5, 0.5;
    Eigen::VectorXdual x2 = x + 0.5 * 0.001;
    std::cout << x2 << "\n";
}

yields the compile error

1>D:\Code\helmderiv\snippets\bugdualvec.cpp(14,31): error C2676:  binary '+': 'Eigen::VectorXdual' does not define this operator or a conversion to a type acceptable to the predefined operator

Readme is a bit out-dated

For example, it feels like that the eigen stuff should be used with include instead of define?

Example: Gradient of a scalar function with parameters

It was my understanding that the example "Gradient of a scalar function with parameters" should be computing the derivative of f wrt x and wrt p
In maths, I expected gx = df/dx = 2*x*exp(sum(p)) and gp=df/dp = p * sum(x*x) * exp(sum(p))
Which on the example x=[1,2,3,4,5] and p=[1,2,3] would yield to one decimal place:
expect_gx = [806.9, 1613.7, 2420.6, 3227.4, 4034.3]
and
expect_gp = [22188.6, 44377.2, 66565.8]

However, when I run the example I get:
gp = [22188.6, 22188.6, 22188.6]

Is this because it is a parameter or is my maths bad?

Let Derivatives (in reverse.hpp) store derivatives of ALL nodes in the expression tree

From the discussion held in #17, it makes sense that Derivatives in reverse.hpp not only stores the derivatives of var variables with respect to leaf var variables (those created out of doubles such as var a = 2.0;) but also intermediate var variables (e.g., var z = x * y; with both x and y being var).

Show Propogation

Is there any way to see what is inside a var or perhaps it's derivative? I keep getting zero for my dud(x) and I want to see what's inside. Any suggestions?

Derivative of pow fails

I tried to do higher-order derivatives of the pow(x,y) function, and that failed with an rather large amount of output to the shell. Other functions (e.g., sin(x)) work just fine with higher-order duals.

Failing example, based on file from tutorial:

// C++ includes
#include <iostream>
using namespace std;

// autodiff include
#include <autodiff/forward.hpp>
using namespace autodiff;

// Define a Nth order dual type using HigherOrderDual<N> construct.
using dualN = HigherOrderDual<4>;
using dualNm1 = HigherOrderDual<3>;

// The single-variable function for which derivatives are needed
dualN f(dualN x, dualN y)
{
    // return sin(y); // OK!
    return pow(x, y);
}

int main()
{
    dualN x = 2.0;      // the input variable x
    dualN y = 1.0;      // the input variable y
    dualN u = f(x, y);  // the output variable u

    dualNm1 ux = derivative(f, wrt(x), at(x, y));  // evaluate the derivative du/dx
    dualNm1 uy = derivative(f, wrt(y), at(x, y));  // evaluate the derivative du/dy
}

with shell output:

Scanning dependencies of target Nthdual
[ 50%] Building CXX object CMakeFiles/Nthdual.dir/Nthdual.cpp.o
In file included from /home/ian/Code/helmderiv/externals/autodiff/autodiff/forward.hpp:33:0,
                 from /home/ian/Code/helmderiv/Nthdual.cpp:6:
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp: In instantiation of ‘constexpr void autodiff::forward::assignPow(autodiff::forward::Dual<T, G>&, U&&) [with T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; U = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&]’:
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:976:18:   required from ‘constexpr void autodiff::forward::assign(autodiff::forward::Dual<T, G>&, U&&) [with T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; U = autodiff::forward::BinaryExpr<autodiff::forward::PowOp, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&>]’
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:471:15:   required from ‘autodiff::forward::Dual<T, G>::Dual(U&&) [with U = autodiff::forward::BinaryExpr<autodiff::forward::PowOp, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&>; typename std::enable_if<(isExpr<U> && (! isDual<U>)), void>::type ...<anonymous> = {}; T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >]’
/home/ian/Code/helmderiv/Nthdual.cpp:17:20:   required from here
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:1296:32: error: no matching function for call to ‘pow(autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >&, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >&)’
         const T aux1 = std::pow(self.val, other.val);
                        ~~~~~~~~^~~~~~~~~~~~~~~~~~~~~
In file included from /home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:33:0,
                 from /home/ian/Code/helmderiv/externals/autodiff/autodiff/forward.hpp:33,
                 from /home/ian/Code/helmderiv/Nthdual.cpp:6:
/usr/include/c++/7/cmath:415:5: note: candidate: template<class _Tp, class _Up> constexpr typename __gnu_cxx::__promote_2<_Tp, _Up>::__type std::pow(_Tp, _Up)
     pow(_Tp __x, _Up __y)
     ^~~
/usr/include/c++/7/cmath:415:5: note:   template argument deduction/substitution failed:
/usr/include/c++/7/cmath: In substitution of ‘template<class _Tp, class _Up> constexpr typename __gnu_cxx::__promote_2<_Tp, _Up>::__type std::pow(_Tp, _Up) [with _Tp = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; _Up = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >]’:
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:1296:32:   required from ‘constexpr void autodiff::forward::assignPow(autodiff::forward::Dual<T, G>&, U&&) [with T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; U = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&]’
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:976:18:   required from ‘constexpr void autodiff::forward::assign(autodiff::forward::Dual<T, G>&, U&&) [with T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; U = autodiff::forward::BinaryExpr<autodiff::forward::PowOp, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&>]’
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:471:15:   required from ‘autodiff::forward::Dual<T, G>::Dual(U&&) [with U = autodiff::forward::BinaryExpr<autodiff::forward::PowOp, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&>; typename std::enable_if<(isExpr<U> && (! isDual<U>)), void>::type ...<anonymous> = {}; T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >]’
/home/ian/Code/helmderiv/Nthdual.cpp:17:20:   required from here
/usr/include/c++/7/cmath:415:5: error: no type named ‘__type’ in ‘struct __gnu_cxx::__promote<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, false>’
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp: In instantiation of ‘constexpr void autodiff::forward::assignPow(autodiff::forward::Dual<T, G>&, U&&) [with T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; U = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&]’:
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:976:18:   required from ‘constexpr void autodiff::forward::assign(autodiff::forward::Dual<T, G>&, U&&) [with T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; U = autodiff::forward::BinaryExpr<autodiff::forward::PowOp, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&>]’
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:471:15:   required from ‘autodiff::forward::Dual<T, G>::Dual(U&&) [with U = autodiff::forward::BinaryExpr<autodiff::forward::PowOp, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&>; typename std::enable_if<(isExpr<U> && (! isDual<U>)), void>::type ...<anonymous> = {}; T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >]’
/home/ian/Code/helmderiv/Nthdual.cpp:17:20:   required from here
/usr/include/c++/7/cmath:392:3: note: candidate: constexpr long double std::pow(long double, long double)
   pow(long double __x, long double __y)
   ^~~
/usr/include/c++/7/cmath:392:3: note:   no known conversion for argument 1 from ‘autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >’ to ‘long double’
/usr/include/c++/7/cmath:388:3: note: candidate: constexpr float std::pow(float, float)
   pow(float __x, float __y)
   ^~~
/usr/include/c++/7/cmath:388:3: note:   no known conversion for argument 1 from ‘autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >’ to ‘float’
In file included from /usr/include/features.h:424:0,
                 from /usr/include/x86_64-linux-gnu/c++/7/bits/os_defines.h:39,
                 from /usr/include/x86_64-linux-gnu/c++/7/bits/c++config.h:533,
                 from /usr/include/c++/7/iostream:38,
                 from /home/ian/Code/helmderiv/Nthdual.cpp:2:
/usr/include/x86_64-linux-gnu/bits/mathcalls.h:140:1: note: candidate: double pow(double, double)
 __MATHCALL_VEC (pow,, (_Mdouble_ __x, _Mdouble_ __y));
 ^
/usr/include/x86_64-linux-gnu/bits/mathcalls.h:140:1: note:   no known conversion for argument 1 from ‘autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >’ to ‘double’
In file included from /home/ian/Code/helmderiv/externals/autodiff/autodiff/forward.hpp:33:0,
                 from /home/ian/Code/helmderiv/Nthdual.cpp:6:
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:1297:32: error: no matching function for call to ‘log(autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >&)’
         const T aux2 = std::log(self.val);
                        ~~~~~~~~^~~~~~~~~~
In file included from /home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:33:0,
                 from /home/ian/Code/helmderiv/externals/autodiff/autodiff/forward.hpp:33,
                 from /home/ian/Code/helmderiv/Nthdual.cpp:6:
/usr/include/c++/7/cmath:350:5: note: candidate: template<class _Tp> constexpr typename __gnu_cxx::__enable_if<std::__is_integer<_Tp>::__value, double>::__type std::log(_Tp)
     log(_Tp __x)
     ^~~
/usr/include/c++/7/cmath:350:5: note:   template argument deduction/substitution failed:
/usr/include/c++/7/cmath: In substitution of ‘template<class _Tp> constexpr typename __gnu_cxx::__enable_if<std::__is_integer<_Tp>::__value, double>::__type std::log(_Tp) [with _Tp = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >]’:
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:1297:32:   required from ‘constexpr void autodiff::forward::assignPow(autodiff::forward::Dual<T, G>&, U&&) [with T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; U = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&]’
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:976:18:   required from ‘constexpr void autodiff::forward::assign(autodiff::forward::Dual<T, G>&, U&&) [with T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; U = autodiff::forward::BinaryExpr<autodiff::forward::PowOp, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&>]’
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:471:15:   required from ‘autodiff::forward::Dual<T, G>::Dual(U&&) [with U = autodiff::forward::BinaryExpr<autodiff::forward::PowOp, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&>; typename std::enable_if<(isExpr<U> && (! isDual<U>)), void>::type ...<anonymous> = {}; T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >]’
/home/ian/Code/helmderiv/Nthdual.cpp:17:20:   required from here
/usr/include/c++/7/cmath:350:5: error: no type named ‘__type’ in ‘struct __gnu_cxx::__enable_if<false, double>’
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp: In instantiation of ‘constexpr void autodiff::forward::assignPow(autodiff::forward::Dual<T, G>&, U&&) [with T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; U = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&]’:
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:976:18:   required from ‘constexpr void autodiff::forward::assign(autodiff::forward::Dual<T, G>&, U&&) [with T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; U = autodiff::forward::BinaryExpr<autodiff::forward::PowOp, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&>]’
/home/ian/Code/helmderiv/externals/autodiff/autodiff/forward/forward.hpp:471:15:   required from ‘autodiff::forward::Dual<T, G>::Dual(U&&) [with U = autodiff::forward::BinaryExpr<autodiff::forward::PowOp, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >, autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > > >&>; typename std::enable_if<(isExpr<U> && (! isDual<U>)), void>::type ...<anonymous> = {}; T = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >; G = autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >]’
/home/ian/Code/helmderiv/Nthdual.cpp:17:20:   required from here
/usr/include/c++/7/cmath:342:3: note: candidate: constexpr long double std::log(long double)
   log(long double __x)
   ^~~
/usr/include/c++/7/cmath:342:3: note:   no known conversion for argument 1 from ‘autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >’ to ‘long double’
/usr/include/c++/7/cmath:338:3: note: candidate: constexpr float std::log(float)
   log(float __x)
   ^~~
/usr/include/c++/7/cmath:338:3: note:   no known conversion for argument 1 from ‘autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >’ to ‘float’
In file included from /usr/include/features.h:424:0,
                 from /usr/include/x86_64-linux-gnu/c++/7/bits/os_defines.h:39,
                 from /usr/include/x86_64-linux-gnu/c++/7/bits/c++config.h:533,
                 from /usr/include/c++/7/iostream:38,
                 from /home/ian/Code/helmderiv/Nthdual.cpp:2:
/usr/include/x86_64-linux-gnu/bits/mathcalls.h:104:1: note: candidate: double log(double)
 __MATHCALL_VEC (log,, (_Mdouble_ __x));
 ^
/usr/include/x86_64-linux-gnu/bits/mathcalls.h:104:1: note:   no known conversion for argument 1 from ‘autodiff::forward::Dual<autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> >, autodiff::forward::Dual<autodiff::forward::Dual<double, double>, autodiff::forward::Dual<double, double> > >’ to ‘double’
CMakeFiles/Nthdual.dir/build.make:62: recipe for target 'CMakeFiles/Nthdual.dir/Nthdual.cpp.o' failed
make[3]: *** [CMakeFiles/Nthdual.dir/Nthdual.cpp.o] Error 1
CMakeFiles/Makefile2:67: recipe for target 'CMakeFiles/Nthdual.dir/all' failed
make[2]: *** [CMakeFiles/Nthdual.dir/all] Error 2
CMakeFiles/Makefile2:79: recipe for target 'CMakeFiles/Nthdual.dir/rule' failed
make[1]: *** [CMakeFiles/Nthdual.dir/rule] Error 2
Makefile:118: recipe for target 'Nthdual' failed
make: *** [Nthdual] Error 2

Implement ConditionExpr for reverse mode algorithm

@ludkinm @supersega

This is a follow up of this discussion.

Consider the following function:

var f(var x) { return (x > 0) ? x**x : x*x*x; }

Consider now what happens in the code below:

var x, u;

x = 1.0; 
u = f(x); // expression tree constructed for x*x

x = -1.0; 
u = f(x); // expression tree constructed for x*x*x

Each function call is re-constructing the expression tree. To be able to avoid this, we need a ConditionExpr and condition function, so that future evaluations of u could be done by just changing x, and then applying an update request on u, for example:

var f(var x) { return condition(x > 0, x**x, x*x*x); }

which would permit us to do:

x = 1.0; 
u = f(x);

x = -1.0; 
u.update(); // this functionality also needs to be implemented

I think with this addition, will truly be able to no longer require expression tree reconstruction, which has happened so far because of potential branching in expressions involving var numbers.

Originally posted by @allanleal in #40 (comment)

Scalar operations on VectorXvar and VectorXd

This results in Eigen compilation errors

#include "Eigen/Core"
#include "autodiff/reverse.hpp"
#include "autodiff/reverse/eigen.hpp"

int main(void){
  Eigen::VectorXvar Xv(10);
  Eigen::VectorXd Y(10);
  autodiff::var z = (Xv -Y).squaredNorm();
  return 0;
}

Originally posted by @ludkinm in #73 (comment)

Add bibtex reference

Hi Allan,
it would be nice to have a bibtex reference in the documentation to just copy paste it and to properly cite this autodiff tool. This could fit in the FAQ or About page.

How do you think about this?

Differentiation of member functions

Hi Allan Leal,
the Autodiff tool is very cool!
I looked into your code and stumbled into the following problem when i wanted to insert Autodiff in my code

When i try to differentiate a member function "foo:energy" using the gradient function. It doesn't work and upcoming warning of my Clion IDE is "reference to non-static member function must be called".

For the full error message see comment below

As a workaround i tried to pass a dummy function to gradient. This only works up to the when i change a variable of another class in the class is changed.

Maybe this issue is by design and the function which i want to differentiate should be Independent of member variables?

In the following i have a MWE appended.

#include <math.h>
#include "Packages/Eigen/Core"
using namespace Eigen;

#include <iostream>
using namespace std;
// autodiff include
#include "Packages/autodiff/forward.hpp"
#include "Packages/autodiff/forward/eigen.hpp"
using namespace autodiff;

class foo2
{
public:
    foo2();
    foo2(double c);
    void increase_d(double f)
    {d += f;};

    double get_d()
    {return d;};


    // Copy constructor
    foo2(const foo2 &p2) {d = p2.d;}

private:
    double d=5;
};

foo2::foo2(double c)
{
    this->d=c;
}

class foo
{
public:
    foo();
    foo(double c,foo2 foo2instance );
    dual energy(VectorXdual);
    VectorXdual calculategradient();

    void increase_a(double b)
    {a += b;};

private:
    double a=1.0;
    foo2*  foo2Infoo = new foo2(2.0);

};
dual foo::energy(VectorXdual x)
{

    return x.dot(x)*a*foo2Infoo->get_d();
}

foo*  foodummy;

foo::foo(double c,foo2 foo2instance)
{

    this->a=c;
    this->foo2Infoo = &foo2instance;

    foodummy=this;
}

dual energydummy(const Eigen::VectorXdual& x)
{

    return foodummy->energy(x);
}

VectorXdual foo::calculategradient()
{
    Eigen::VectorXdual xD(3);
    xD << 1,2,3;
    dual dualenergy;

    //Eigen::VectorXd gx = gradient(energy, wrt(xD), at(xD), dualenergy); // gives "reference to non-static member function must be called" error
    Eigen::VectorXd gx = gradient(energydummy, wrt(xD), at(xD), dualenergy);


    return gx;
}

int main()
{
   foo2 foo2instance(1);
   foo fooinstance(2,foo2instance);
    cout<<"Gradient at the start"<< endl;
    cout<<fooinstance.calculategradient()<< endl;
    fooinstance.increase_a(2.0);


    cout<<"Gradient after increasing value of member double:"<< endl;
    cout<<fooinstance.calculategradient()<< endl;

    foo2instance.increase_d(1.0);
    cout<<"Gradient after increasing double of reference member :"<< endl;
    cout<<fooinstance.calculategradient()<< endl;

    return 0;
}

Thanks for such a nice Package.

Inaccuracies in the processing of a special case for a power function

Hello, thanks for developing a great math library.
Below I will describe the problem that I discovered and which seems to me to be a bug.

#include <iostream>
#include <Eigen/Core>
#include <autodiff/reverse.hpp>

using namespace std;
using namespace autodiff;

var f1(var x)
{
    return x*x;
}

var f2(var x)
{
    return pow(x, 2);
}

int main(){
    var x = 0.0;                           // zero is special case but when power >= 1 derivative should exist
    var u1 = f1(x);                        // f1 = x * x
    var u2 = f2(x);                        // f2 = pow(x, 2)

    Derivatives dud1 = derivatives(u1);
    Derivatives dud2 = derivatives(u2);

    var dudx1 = dud1(x);
    var dudx2 = dud2(x);

    cout << "du/dx1 = " << dudx1 << endl;  // here du/dx = 0, that's OK
    cout << "du/dx2 = " << dudx2 << endl;  // but here du/dx = -nan, that's possible bug
}

This simple program gives different results, although the functions are actually the same. Perhaps the problem is that I need an implementation of the function of raising to the INTEGER power, but I did not find it when I looked at the documentation fluently. However, in any case, in the zero argument with powers of at least 1, the derivative (at least one-sided) still exists and it would be nice to give it out.

Incorrect result when variable self assignment is used

I have noticed the following behaviour when assigning the result of an operation to a variable used in the operation itself:

#include <autodiff/forward.hpp>

int main() {

	using DoubleT = dual;
	DoubleT a = 1;
	a = a - 2*a;
	cout << "it's " << a << " should be -1\n"; // dual returns -4

	a = 1;
	a = a + 2 * a;
	cout << "it's " << a << " should be 3\n";  // dual returns 4

	a = 1;
	a = a - a;
	cout << "it's " << a << " should be 0\n";  // dual returns -2

	a = 1;
	a = a + a;
	cout << "it's " << a << " should be 2\n"; // dual returns 2

	a = 1;
	a = 2*a - a;
	cout << "it's " << a << " should be 1\n";  // dual returns -3

	a = 1;
	a = 2*a + a;
	cout << "it's " << a << " should be 3\n";  // dual returns 3

	a = 1;
	a = DoubleT{ 2 * a - a };
	cout << "it's " << a << " should be 1\n";  // dual correctly returns 1
}

cmake error: install DIRECTORY given no DESTINATION!

cmake .. -DCMAKE_INSTALL_PREFIX=/some/local/dir
CMake Error at CMakeLists.txt:20 (install):
install DIRECTORY given no DESTINATION!

CMake Warning (dev) in CMakeLists.txt:
No cmake_minimum_required command is present. A line of code such as

cmake_minimum_required(VERSION 3.10)

should be added at the top of the file. The version specified may be lower
if you wish to support older CMake versions for this project. For more
information run "cmake --help-policy CMP0000".
This warning is for project developers. Use -Wno-dev to suppress it.

-- Configuring incomplete, errors occurred!
See also "/mnt/d/research/Shuang/auto_diff/autodiff/autodiff/CMakeFiles/CMakeOutput.log".

Support for complex and matrices

Hi!

I was wondering if you plan to support complex duals in the future? I wanted to use your library, but this would be a big requirement for me.

Thanks for this amazing work! Cheers,
David

Multi-level derivatives (and seeding)

I'm getting strange behavior when trying to calculate a gradient. I wonder if it's because in my objective function I use calculate derivatives.

To explain: I see the way one calculates involves seeding autodiff with the target variable:

VectorXdual u(5);
u << 1.0, 1.1, 1.2, 1.3, 1.4;

dual w = energy(u);

VectorXd dwdu(5);
int i; for (i=0; i<5; i++)
{
    seed(wrt(u[i]));
    dwdu[i] = energy(u).grad;
    unseed(wrt(u[i]));
}

So if I understand, seeding is a way to (globally?) track the variable of interest. But what happens if my target function also calculates derivatives?

dual f(VectorXdual u)
{
   // ....
}

dual energy(VectorXdual u)
{
  double dfdx = derivative(f, wrt(x), u);
  // ...
}

I see in forward.hpp that derivative() also uses seeding. Do the two clash with one another? Can a nested derivative calculation like this work?

Potential bug when using a unary minus

I've spotted the following bug when using the dual type in forward.hpp

dual f_bug(dual x)
{
  return -x - 2.0 * log(1.0 + exp(-x));
}

Running f(1.0) should result in -1.62652... however I get -1.0

Adding parentheses fixes the problem:

dual f_ok(dual x)
{
  return -(x + 2.0 * log(1.0 + exp(-x)));
}

Is there an issue with the order of operations? Or is this expected behaviour?

Conversion to Double

I have a function which returns a var. How do I then do a calculation such as fabs(a-0.001)? I'm getting the error

error: cannot convert ‘autodiff::var’ to ‘double’ in assignment

No header guards?

Causes compiler to spit redefinition errors.

In Real<N, T> type, ensure T is a numeric type

In the definition of type Real<N, T>:

template<size_t N, typename T>
class Real
{
...
};

Ensure T is a numeric type (not necessarily standard number types such as float, double, but also permit other custom numbers to be used - maybe even var or dual, or emulated high-precision for floating-point types).

Modify NumberTraits in numbertraits.hpp to allow one to register a number type that autodiff can understand.

Originally posted by @supersega in #69

Gradient of General Function in Forward Mode

In forward mode to take the gradient of a scalar you define some function f which returns a dual and takes in a VectorXdual

dual f(const VectorXdual& x)
{
    // ...
}

The gradient is then taken with respect to another VectorXdual

VectorXdual x = ...;
VectorXd dudx = gradient(f, x);

What happens if my function f needs to take in more than just a VectorXdual?

dual f(const VectorXdual& x, double a, int b)
{
    // ...
}

How do I take the gradient then?

Valgrind Uninitialised Value

I'm getting a strange error on valgrind. Here is the code:

mesh.nx = (var*) malloc(sizeof(var)*mesh.nn);
mesh.ny = (var*) malloc(sizeof(var)*mesh.nn);


    real dx = w/(2*wn);
    real dy = h/(2*hn);

    int xi, yi;
    for (xi=0; xi<=2*wn; xi++)
      for (yi=0; yi<=2*hn; yi++)
      {
        int i = yi + (2*hn+1) * xi;
        mesh.nx[i] = xi * dx;
        mesh.ny[i] = yi * dy;
      }

And here is what valgrind says:

==6502== Conditional jump or move depends on uninitialised value(s)
==6502==    at 0x163634: std::__shared_count<(__gnu_cxx::_Lock_policy)2>::~__shared_count() (shared_ptr_base.h:707)
==6502==    by 0x160FFD: std::__shared_ptr<autodiff::reverse::Expr const, (__gnu_cxx::_Lock_policy)2>::~__shared_ptr() (shared_ptr_base.h:1147)
==6502==    by 0x16587F: std::__shared_ptr<autodiff::reverse::Expr const, (__gnu_cxx::_Lock_policy)2>::operator=(std::__shared_ptr<autodiff::reverse::Expr const, (__gnu_cxx::_Lock_policy)2>&&) (shared_ptr_base.h:1243)
==6502==    by 0x163673: std::shared_ptr<autodiff::reverse::Expr const>::operator=(std::shared_ptr<autodiff::reverse::Expr const>&&) (shared_ptr.h:335)
==6502==    by 0x162FF6: autodiff::var::operator=(autodiff::var&&) (reverse.hpp:649)
==6502==    by 0x1597B7: square_tri(double, double, int, int) (mesh.cpp:46)
==6502==    by 0x160405: main (energy.cpp:98)

Any idea what it means?

Python wrapper

As I am a hybrid Python and C++ developer (primarily), I was thinking about bringing autodiff into Python, and I started playing with pybind11 to write a wrapper, but that wasn't very successful. Wrapping individual objects like dual was fine, but dealing with the expression outputs for operations on duals wasn't super obvious how to handle. I haven't entirely given up on this topic, but it would be really slick to be able to use this library in Python.

Has anybody given any thought to this topic? I know it is not trivial, but if possible, I would be keen to give it a shot.