Commit a2cb8a32 authored by Michael Tesch's avatar Michael Tesch
Browse files

add some rough finite-diff tests, minor cleanups to 'dual' 

also re-enable stream teting functions for msvc.
parent efe41a97

.appveyor.yml

0 → 100644
+37 −0
Original line number Diff line number Diff line 
version: 0.2.{build}

clone_folder: c:\projects\cppduals



image:

#- Visual Studio 2013

#- Visual Studio 2015

- Visual Studio 2017



configuration:

- Release

#- Debug # someone with a debugger please investigate this :)



# skip unsupported combinations

init:

- echo %APPVEYOR_BUILD_WORKER_IMAGE%

- if "%APPVEYOR_BUILD_WORKER_IMAGE%"=="Visual Studio 2013" ( set generator="Visual Studio 12 2013" )

- if "%APPVEYOR_BUILD_WORKER_IMAGE%"=="Visual Studio 2015" ( set generator="Visual Studio 14 2015" )

- if "%APPVEYOR_BUILD_WORKER_IMAGE%"=="Visual Studio 2017" ( set generator="Visual Studio 15 2017" )

- echo %generator%



before_build:

- cmd: |-

    mkdir build

    cd build

    cmake --version

    cmake .. -G %generator% -DCPPDUALS_TESTING=ON

    pwd

    ls



build:

  project: c:\projects\cppduals\build\cppduals.sln

  verbosity: minimal

#  parallel: true



test_script:

- pwd

- ctest -C Debug -VV

.gitlab-ci.yml

+4 −1
Original line number Diff line number Diff line
@@ -40,6 +40,8 @@ test: 
  - build



cover:

  variables:

    CTEST_OUTPUT_ON_FAILURE: y

  script:

  - cmake -Bbuild-cov -H. -DCODE_COVERAGE=ON -DCPPDUALS_TESTING=ON

  - cmake --build build-cov --target cov
@@ -53,10 +55,11 @@ cover: 
  - merge_requests



pages:

  variables:

    CTEST_OUTPUT_ON_FAILURE: y

  script:

  - cmake -Bbuild -H. -DCODE_COVERAGE=ON -DCPPDUALS_TESTING=ON

  - cmake --build build --target docs

  - cmake --build build --target cov

  - cmake --build build --target cov-html

  - mv build/docs public/

  - mv build/coverage public/

CMakeLists.txt

+5 −2
Original line number Diff line number Diff line
@@ -6,9 +6,8 @@ project (cppduals C CXX) 


#set (CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})

set (CPPDUALS_VERSION 0.1.2)

set (CMAKE_CXX_STANDARD 11)

set (CMAKE_CXX_STANDARD 11 CACHE STRING "Which C++ standard to test against.")

set (CMAKE_CXX_STANDARD_REQUIRED ON)

set (CMAKE_CXX_EXTENSIONS OFF)

set (CMAKE_DISABLE_IN_SOURCE_BUILD ON)

if (NOT CMAKE_CONFIGURATION_TYPES AND

    NOT CMAKE_NO_BUILD_TYPE AND
@@ -19,6 +18,8 @@ if (NOT CMAKE_CONFIGURATION_TYPES AND 
    STRINGS "Debug" "Release" "MinSizeRel" "RelWithDebInfo")

endif()



set_property (CACHE CMAKE_CXX_STANDARD PROPERTY STRINGS 11 14 17 20)



option (CPPDUALS_TESTING "Enable testing" OFF)

option (CPPDUALS_BENCHMARK "Enable benchmarking" OFF)

option (CPPDUALS_EIGEN_LATEST "Eigen latest" OFF)
@@ -164,3 +165,5 @@ if (ETAGS) 
    )

  add_custom_target (etags DEPENDS tags)

endif (ETAGS)



message ("CMAKE_CXX_FLAGS: ${CMAKE_CXX_FLAGS}")

README.md

+8 −0
Original line number Diff line number Diff line
@@ -6,6 +6,9 @@ used for automatic (forward) differentiation. It can be used in 
conjunction with Eigen to produced very fast vectorized computations

of real and complex matrix functions and their derivatives.



There is a small paper on cppduals here:

[![DOI](https://joss.theoj.org/papers/10.21105/joss.01487/status.svg)](https://doi.org/10.21105/joss.01487)



Documentation

=============
@@ -199,6 +202,11 @@ Questions, bug reports, bug fixes, and contributions are welcome. 
Simply submit an [Issue](https://gitlab.com/tesch1/cppduals/issues)

or [Merge Request](https://gitlab.com/tesch1/cppduals/merge_requests).



Contributors

------------



- [Nestor Demeure](https://gitlab.com/nestordemeure)



Compiler notes

==============

duals/dual

+86 −33
Original line number Diff line number Diff line
@@ -133,6 +133,9 @@ template<class T> class dual { 
    void rpart(T);

    void dpart(T);



    dual<T> operator-() const;

    dual<T> operator+() const;



    dual<T> & operator= (const T &);

    dual<T> & operator+=(const T &);

    dual<T> & operator-=(const T &);
@@ -178,6 +181,7 @@ dual<T> pow(dual<T>, U) 
dual<T> pow(U, dual<T>)

dual<T> pow(dual<T>, dual<T>)

dual<T> sqrt(dual<T>)

dual<T> cbrt(dual<T>)

dual<T> sin(dual<T>)

dual<T> cos(dual<T>)

dual<T> tan(dual<T>)
@@ -520,9 +524,12 @@ public: 
  /// Get the dual part.

  void dpart(value_type du) { _dual = du; }



  /// Unitary negation

  /// Unary negation

  dual<T> operator-() const { return dual<T>(-_real, -_dual); }



  /// Unary nothing

  dual<T> operator+() const { return *this; }



  /// Assignment of `value_type` assigns the real part and zeros the dual part.

  dual<T> & operator= (const T & x) { _real =  x; _dual = value_type(); return *this; }
@@ -782,8 +789,12 @@ template<class T> dual<T> exp(const dual<T> & x) { 
/// Natural log ln(x)

template<class T> dual<T> log(const dual<T> & x) {

  using std::log;

  return dual<T>(log(x.rpart()),

  T v = log(x.rpart());

  if (x.dpart())

    return dual<T>(v,

                   x.dpart() / x.rpart());

  else

    return v;

}



template<class T> dual<T> log10(const dual<T> & x) {
@@ -830,20 +841,20 @@ template <typename T> int sgn(T val) { 
template<class T> dual<T> abs(const dual<T> & x) {

  using std::abs;

  return dual<T>(abs(x.rpart()),

                 utils::sgn(x.rpart()) * x.dpart());

                 x.dpart() * utils::sgn(x.rpart()));

}



template<class T> dual<T> fabs(const dual<T> & x) {

  using std::fabs;

  return dual<T>(fabs(x.rpart()),

                 utils::sgn(x.rpart()) * x.dpart());

                 x.dpart() * utils::sgn(x.rpart()));

}



#if 0

template<class T> dual<T> abs2(const dual<T> & x) {

  using std::abs;

  return dual<T>(x.rpart() * x.rpart(),

                 xxx utils::sgn(x.rpart()) * x.dpart());

                 xxx x.dpart() * utils::sgn(x.rpart()));

}

#endif
@@ -883,15 +894,26 @@ template<class T> bool (signbit)(const duals::dual<T> & d) { using std::signbit; 
template<class T> dual<T> sqrt(const dual<T> & x) {

  using std::sqrt;

  T v = sqrt(x.rpart());

  return dual<T>(v,

                 T(1) / (T(2) * v) * x.dpart());

  if (x.dpart())

    return dual<T>(v, x.dpart() / (T(2) * v) );

  else

    return v;

}



template<class T> dual<T> cbrt(const dual<T> & x) {

  using std::cbrt;

  T v = cbrt(x.rpart());

  if (x.dpart())

    return dual<T>(v, x.dpart() / (T(3) * v * v) );

  else

    return v;

}



template<class T> dual<T> sin(const dual<T> & x) {

  using std::sin;

  using std::cos;

  return dual<T>(sin(x.rpart()),

                 cos(x.rpart()) * x.dpart());

                 x.dpart() * cos(x.rpart()));

}



template<class T> dual<T> cos(const dual<T> & x) {
@@ -903,35 +925,48 @@ template<class T> dual<T> cos(const dual<T> & x) { 


template<class T> dual<T> tan(const dual<T> & x) {

  using std::tan;

  T y = tan(x.rpart());

  return dual<T>(y,

                 (y*y + 1) * x.dpart());

  T v = tan(x.rpart());

  return dual<T>(v, x.dpart() * (v*v + 1));

}



template<class T> dual<T> asin(const dual<T> & x) {

  using std::asin;

  using std::sqrt;

  return dual<T>(asin(x.rpart()),

                 1 / sqrt(1 - x.rpart()*x.rpart()) * x.dpart());

  T v = asin(x.rpart());

  if (x.dpart())

    return dual<T>(v, x.dpart() / sqrt(1 - x.rpart()*x.rpart()));

  else

    return v;

}



template<class T> dual<T> acos(const dual<T> & x) {

  using std::acos;

  using std::sqrt;

  return dual<T>(acos(x.rpart()),

                 -1 / sqrt(1 - x.rpart()*x.rpart()) * x.dpart());

  T v = acos(x.rpart());

  if (x.dpart())

    return dual<T>(v, -x.dpart() / sqrt(1 - x.rpart()*x.rpart()));

  else

    return v;

}



template<class T> dual<T> atan(const dual<T> & x) {

  using std::atan;

  return dual<T>(atan(x.rpart()),

                 1 / (1 + x.rpart()*x.rpart()) * x.dpart());

  T v = atan(x.rpart());

  if (x.dpart())

    return dual<T>(v,

                   x.dpart() / (1 + x.rpart()*x.rpart()));

  else

    return v;

}



template<class T> dual<T> atan2(const dual<T> & x, const dual<T> & y) {

  using std::atan2;

  return dual<T>(atan2(x.rpart(), y.rpart()),

                 1 / (1 + x.rpart()*x.rpart()) * x.dpart());

  T v = atan2(x.rpart(), y.rpart());

  if (x.dpart())

    return dual<T>(v,

                   x.dpart() / (1 + x.rpart()*x.rpart()));

  else

    return v;

}



// TODO
@@ -951,8 +986,12 @@ template<class T> dual<T> erf(const dual<T> & x) { 
  using std::sqrt;

  using std::pow;

  using std::exp;

  return dual<T>(erf(x.rpart()),

                 T(2)/sqrt(T(MY_PI)) * exp(-pow(x.rpart(),T(2))) * x.dpart());

  T v = erf(x.rpart());

  if (x.dpart())

    return dual<T>(v,

                   x.dpart() * T(2)/sqrt(T(MY_PI)) * exp(-pow(x.rpart(),T(2))));

  else

    return v;

}



/// Error function complement (1 - erf()).
@@ -961,25 +1000,38 @@ template<class T> dual<T> erfc(const dual<T> & x) { 
  using std::sqrt;

  using std::pow;

  using std::exp;

  return dual<T>(erfc(x.rpart()),

                 -T(2)/sqrt(T(MY_PI)) * exp(-pow(x.rpart(),T(2))) * x.dpart());

  T v = erfc(x.rpart());

  if (x.dpart())

    return dual<T>(v,

                   x.dpart() * -T(2)/sqrt(T(MY_PI)) * exp(-pow(x.rpart(),T(2))));

  else

    return v;

}



/// Gamma function.  Approximation of the dual part.

// TODO specialize for integers

template<class T> dual<T> tgamma(const dual<T> & x) {

  using std::tgamma;

  T v = tgamma(x.rpart());

  if (x.dpart()) {

    T h(T(1) / (1ull << (std::numeric_limits<T>::digits / 3)));

  return dual<T>(tgamma(x.rpart),

                 (tgamma(x.rpart()+h) - tgamma(x.rpart()))/h * x.dpart());

    return dual<T>(v,

                   x.dpart() * (tgamma(x.rpart()+h) - tgamma(x.rpart()-h))/(2*h));

  }

  else

    return v;

}



/// Log of absolute value of gamma function.  Approximation of the dual part.

template<class T> dual<T> lgamma(const dual<T> & x) {

  using std::lgamma;

  T h(T(1) / (1ull << (std::numeric_limits<T>::digits / 3)));

  return dual<T>(lgamma(x.rpart()),

                 (lgamma(x.rpart()) - lgamma(x.rpart()+h))/h * x.dpart());

  T v = lgamma(x.rpart());

  if (x.dpart())

    return dual<T>(v,

                   x.dpart() * (lgamma(x.rpart()+h) - lgamma(x.rpart()-h))/(2*h));

  else

    return v;

}



/// Putto operator
@@ -994,8 +1046,9 @@ operator<<(std::basic_ostream<_CharT, _Traits> & os, const dual<T> & x) 
  s << '(' << x.rpart()

    << (x.dpart() < 0 ? "" : "+")

    << x.dpart()

    << "_e" << (sizeof(T) == sizeof(float) ? "f" :

                sizeof(T) == sizeof(long double) ? "l" : "")

    << "_e" << (std::is_same<typename std::decay<T>::type,float>::value ? "f" :

                std::is_same<typename std::decay<T>::type,double>::value ? "" :

                std::is_same<typename std::decay<T>::type,long double>::value ? "l" : "")

    << ")";

  return os << s.str();

}
@@ -1028,8 +1081,8 @@ operator>>(std::basic_istream<CharT, Traits> & is, dual<T> & x) 
              if (c == CharT('e')) {

                is.get();

                c = is.peek();

                if ((c == 'f' && sizeof(T) != sizeof(float)) ||

                    (c == 'l' && sizeof(T) != sizeof(long double)))

                if ((c == 'f' && !std::is_same<typename std::decay<T>::type,float>::value) ||

                    (c == 'l' && !std::is_same<typename std::decay<T>::type,long double>::value))

                  is.setstate(std::ios_base::failbit);

                else {

                  if (c == 'f' || c == 'l')