How to take second derivatives with MetaPhysicL

.gitignore

*~*
*.out
*dSYM

dynamic-sparse-second-derivatives.cpp

#include "metaphysicl/dualdynamicsparsenumbervector.h"
#include <iostream>

#define VectorUnitVector DynamicSparseNumberVectorUnitVector
#define NDIM 2

using namespace MetaPhysicL;

typedef DualNumber<
    DualNumber<double, DynamicSparseNumberVector<double, unsigned int>>,
    DynamicSparseNumberVector<
        DualNumber<double, DynamicSparseNumberVector<double, unsigned int>>,
        unsigned int>>
    SecondDerivType;
typedef double RawScalar;

int main() {

  typedef VectorUnitVector<NDIM, 0, RawScalar>::type XVector;
  XVector xvec = VectorUnitVector<NDIM, 0, RawScalar>::value();

  typedef VectorUnitVector<NDIM, 1, RawScalar>::type YVector;
  YVector yvec = VectorUnitVector<NDIM, 1, RawScalar>::value();

  // Construct "x" variable
  SecondDerivType xn(4);
  xn.value().derivatives().insert<0>() = 1;
  xn.derivatives().insert<0>() = 1;

  // Construct "y" variable
  SecondDerivType yn(5);
  yn.value().derivatives().insert<1>() = 1;
  yn.derivatives().insert<1>() = 1;

  auto x2y3 = xn * xn * yn * yn * yn;
}

second-derivatives.cpp

#include "metaphysicl/dualnumberarray.h"
#include <iostream>

using namespace MetaPhysicL;

int main() {

  // DualNumber takes two typename template arguments. The first template
  // argument specifies the type yielded by DualNumber<>::value(). The second
  // template argument specifies the type yielded by
  // DualNumber<>::derivatives(). If the second template argument is not
  // provided, then it defaults to the type of the first, e.g we assume that the
  // derivative and value type are the same. As a first demonstration of second
  // derivative capability we will consider a single variable system. To do
  // second derivatives, the only real clever bit is that we nest DualNumbers:

  DualNumber<DualNumber<double>> x = 5;

  // Unfortunately, calculation of first derivatives is currently redundant so
  // we also have to redundantly seed them
  x.derivatives().value() = 1;
  x.value().derivatives() = 1;

  auto x2 = x * x;
  std::cout << "second derivative is " << x2.derivatives().derivatives()
            << std::endl;
  std::cout << "first derivative is " << x2.derivatives().value() << std::endl;
  std::cout << "first derivative again is " << x2.value().derivatives()
            << std::endl;
  std::cout << "function value is " << x2.value().value() << std::endl;

  // Now we'll consider a multi-variate system. The typing here probably looks
  // terrifying. We want foo.value().value() to be a single scalar, while we
  // want foo.value().derivatives(), foo.derivatives().value(), and
  // foo.derivatives().derivatives() to all yield NumberArrays. With the below
  // typing, foo.value() yields a type of DualNumber<double, NumberArray<2,
  // double>>, and so consequently, foo.value().value() will yield a double and
  // foo.value().derivatives() will yield a NumberArray<2, double> so that's
  // what we need! foo.derivatives() yields a type of DualNumber<NumberArray<2,
  // double>, NumberArray<2, double>> so foo.derivatives().value() and
  // foo.derivatives().derivatives() will both yield NumberArrays like we want

  // Construct "x" variable
  DualNumber<DualNumber<double, NumberArray<2, double>>,
             DualNumber<NumberArray<2, double>, NumberArray<2, double>>>
      xn = 2;

  // Here we seed the zeroth component of the first derivatives vector for our
  // "x" variable
  xn.derivatives().value()[0] = 1.;
  xn.value().derivatives()[0] = 1.;

  // Construct "y" variable
  DualNumber<DualNumber<double, NumberArray<2, double>>,
             DualNumber<NumberArray<2, double>, NumberArray<2, double>>>
      yn = 3;

  // Here we seed the first component of the first derivatives vector for our
  // "y" variable
  yn.derivatives().value()[1] = 1.;
  yn.value().derivatives()[1] = 1.;

  auto x2y3 = xn * xn * yn * yn * yn;
  std::cout << "second derivatives are " << x2y3.derivatives().derivatives()
            << std::endl;
  std::cout << "first derivatives are " << x2y3.derivatives().value()
            << std::endl;
  std::cout << "first derivatives again are " << x2y3.value().derivatives()
            << std::endl;
  std::cout << "function value is " << x2y3.value().value() << std::endl;

  // Finally let's look at how to obtain cross-derivatives! Courtesy @roystgnr for the typing

  // Construct "x" variable
  DualNumber<DualNumber<double, NumberArray<2, double>>,
             NumberArray<2, DualNumber<double, NumberArray<2, double>>>>
      xc = 4;

  // Here we seed the zeroth component of the first derivatives vector for our
  // "x" variable
  xc.derivatives()[0].value() = 1.;
  xc.value().derivatives()[0] = 1.;

  // Construct "y" variable
  DualNumber<DualNumber<double, NumberArray<2, double>>,
             NumberArray<2, DualNumber<double, NumberArray<2, double>>>>
      yc = 5;

  // Here we seed the first component of the first derivatives vector for our
  // "y" variable
  yc.derivatives()[1].value() = 1.;
  yc.value().derivatives()[1] = 1.;

  auto x2y3c = xc * xc * yc * yc * yc;
  std::cout << "xx derivative is " << x2y3c.derivatives()[0].derivatives()[0]
            << std::endl;
  std::cout << "xy derivative is " << x2y3c.derivatives()[0].derivatives()[1]
            << std::endl;
  std::cout << "yx derivative is " << x2y3c.derivatives()[1].derivatives()[0]
            << std::endl;
  std::cout << "yy derivative is " << x2y3c.derivatives()[1].derivatives()[1]
            << std::endl;

  std::cout << "first derivatives are {" << x2y3c.derivatives()[0].value()
            << "," << x2y3c.derivatives()[1].value() << "}" << std::endl;
  std::cout << "first derivatives again are " << x2y3c.value().derivatives()
            << std::endl;
  std::cout << "function value is " << x2y3c.value().value() << std::endl;
}

Output of above code:

second derivative is 2
first derivative is 10
first derivative again is 10
function value is 25
second derivatives are {54,72}
first derivatives are {108,108}
first derivatives again are {108,108}
function value is 108
xx derivative is 250
xy derivative is 600
yx derivative is 600
yy derivative is 480
first derivatives are {1000,1200}
first derivatives again are {1000,1200}
function value is 2000