Skip to content

Commit

Permalink
Fix higher order derivatives (#48)
Browse files Browse the repository at this point in the history 
* Add pip as dep

* Fix aliasing bug when computing higher order derivatives.

* Increment version
  • Loading branch information
@allanleal
allanleal committed Jul 29, 2019
1 parent 0d4d79d commit c62ddc2
Show file tree
Hide file tree
Showing 4 changed files with 36 additions and 16 deletions.
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Expand Up @@ -8,7 +8,7 @@ set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
include(CCache)

# Name and details of the project
project(autodiff VERSION 0.5.5 LANGUAGES CXX)
project(autodiff VERSION 0.5.7 LANGUAGES CXX)

# Include the cmake variables with values for installation directories
include(GNUInstallDirs)
Expand Down
26 changes: 13 additions & 13 deletions autodiff/forward/forward.hpp
View file
Expand Up @@ -1152,9 +1152,9 @@ constexpr void assignMul(Dual<T, G>& self, U&& other)
}
// ASSIGN-MULTIPLY A DUAL NUMBER: self *= dual
else if constexpr (isDual<U>) {
const auto aux = self.val * other.grad; // to avoid aliasing when self === other
const G aux = other.grad; // to avoid aliasing when self === other
self.grad *= other.val;
self.grad += aux;
self.grad += self.val * aux;
self.val *= other.val;
}
// ASSIGN-MULTIPLY A NEGATIVE EXPRESSION: self *= (-expr)
Expand Down Expand Up @@ -1219,7 +1219,7 @@ constexpr void assignDiv(Dual<T, G>& self, U&& other)
}
// ASSIGN-DIVIDE A DUAL NUMBER: self /= dual
else if constexpr (isDual<U>) {
const auto aux = One<T> / other.val; // to avoid aliasing when self === other
const T aux = One<T> / other.val; // to avoid aliasing when self === other
self.val *= aux;
self.grad -= self.val * other.grad;
self.grad *= aux;
Expand Down Expand Up @@ -1287,14 +1287,14 @@ constexpr void assignPow(Dual<T, G>& self, U&& other)
{
// ASSIGN-POW A NUMBER: self = pow(self, number)
if constexpr (isNumber<U>) {
const auto aux = std::pow(self.val, other - 1);
const T aux = std::pow(self.val, other - 1);
self.grad *= other * aux;
self.val = aux * self.val;
}
// ASSIGN-POW A DUAL NUMBER: self = pow(self, dual)
else if constexpr (isDual<U>) {
const auto aux1 = std::pow(self.val, other.val);
const auto aux2 = std::log(self.val);
const T aux1 = std::pow(self.val, other.val);
const T aux2 = std::log(self.val);
self.grad *= other.val/self.val;
self.grad += aux2 * other.grad;
self.grad *= aux1;
Expand Down Expand Up @@ -1350,31 +1350,31 @@ constexpr void apply(Dual<T, G>& self, CosOp)
template<typename T, typename G>
constexpr void apply(Dual<T, G>& self, TanOp)
{
const auto aux = One<T> / std::cos(self.val);
const T aux = One<T> / std::cos(self.val);
self.val = tan(self.val);
self.grad *= aux * aux;
}

template<typename T, typename G>
constexpr void apply(Dual<T, G>& self, ArcSinOp)
{
const auto aux = One<T> / sqrt(1.0 - self.val * self.val);
const T aux = One<T> / sqrt(1.0 - self.val * self.val);
self.val = asin(self.val);
self.grad *= aux;
}

template<typename T, typename G>
constexpr void apply(Dual<T, G>& self, ArcCosOp)
{
const auto aux = -One<T> / sqrt(1.0 - self.val * self.val);
const T aux = -One<T> / sqrt(1.0 - self.val * self.val);
self.val = acos(self.val);
self.grad *= aux;
}

template<typename T, typename G>
constexpr void apply(Dual<T, G>& self, ArcTanOp)
{
const auto aux = One<T> / (1.0 + self.val * self.val);
const T aux = One<T> / (1.0 + self.val * self.val);
self.val = atan(self.val);
self.grad *= aux;
}
Expand All @@ -1389,7 +1389,7 @@ constexpr void apply(Dual<T, G>& self, ExpOp)
template<typename T, typename G>
constexpr void apply(Dual<T, G>& self, LogOp)
{
const auto aux = One<T> / self.val;
const T aux = One<T> / self.val;
self.val = log(self.val);
self.grad *= aux;
}
Expand All @@ -1398,7 +1398,7 @@ template<typename T, typename G>
constexpr void apply(Dual<T, G>& self, Log10Op)
{
constexpr T ln10 = 2.3025850929940456840179914546843;
const auto aux = One<T> / (ln10 * self.val);
const T aux = One<T> / (ln10 * self.val);
self.val = log10(self.val);
self.grad *= aux;
}
Expand All @@ -1413,7 +1413,7 @@ constexpr void apply(Dual<T, G>& self, SqrtOp)
template<typename T, typename G>
constexpr void apply(Dual<T, G>& self, AbsOp)
{
const auto aux = self.val;
const T aux = self.val;
self.val = abs(self.val);
self.grad *= aux / self.val;
}
Expand Down
1 change: 1 addition & 0 deletions environment.devenv.yml
View file
Expand Up @@ -6,6 +6,7 @@ dependencies:
- eigen
- gxx_linux-64=7.3.0 # [linux]
- ccache # [unix]
- pip
- pip:
- mkdocs
- mkdocs-material
Expand Down
23 changes: 21 additions & 2 deletions test/forward.test.cpp
View file
Expand Up @@ -620,11 +620,13 @@ TEST_CASE("autodiff::dual tests", "[dual]")
{
using dual2nd = HigherOrderDual<2>;

std::function<dual2nd(dual2nd, dual2nd)> f;

dual2nd x = 5;
dual2nd y = 7;

// Testing complex function involving sin and cos
auto f = [](dual2nd x, dual2nd y) -> dual2nd
// Testing simpler functions
f = [](dual2nd x, dual2nd y) -> dual2nd
{
return x*x + x*y + y*y;
};
Expand All @@ -634,6 +636,23 @@ TEST_CASE("autodiff::dual tests", "[dual]")
REQUIRE( derivative(f, wrt(x, y), at(x, y)) == Approx(1.0) );
REQUIRE( derivative(f, wrt(y, x), at(x, y)) == Approx(1.0) );
REQUIRE( derivative(f, wrt(y, y), at(x, y)) == Approx(2.0) );

// Testing complex function involving log
x = 2.0;
y = 1.0;

f = [](dual2nd x, dual2nd y) -> dual2nd
{
return 1 + x + y + x * y + y / x + log(x / y);
};

REQUIRE( val(f(x, y)) == Approx( 1 + val(x) + val(y) + val(x) * val(y) + val(y) / val(x) + log(val(x) / val(y)) ) );
REQUIRE( val(derivative(f, wrt(x), at(x, y))) == Approx( 1 + val(y) - val(y) / (val(x) * val(x)) + 1.0/val(x) - log(val(y)) ) );
REQUIRE( val(derivative(f, wrt(y), at(x, y))) == Approx( 1 + val(x) + 1.0 / val(x) - 1.0/val(y) ) );
REQUIRE( val(derivative(f, wrt(x, x), at(x, y))) == Approx( 2 * val(y) / (val(x) * val(x) * val(x)) + -1.0/(val(x) * val(x)) ) );
REQUIRE( val(derivative(f, wrt(y, y), at(x, y))) == Approx( 1.0/(val(y)*val(y)) ) );
REQUIRE( val(derivative(f, wrt(y, x), at(x, y))) == Approx( 1 - 1.0 / (val(x) * val(x)) ) );
REQUIRE( val(derivative(f, wrt(x, y), at(x, y))) == Approx( 1 - 1.0 / (val(x) * val(x)) ) );
}

SECTION("testing higher order derivatives")
Expand Down

0 comments on commit c62ddc2

Please sign in to comment.