Fixes and unit tests

Fix hessian calculation for sin and cos operations.

Added few unit tests for automatic differentiation

Source/Headers/inmost_expression.h

@@ -1045,7 +1045,10 @@ namespace INMOST
     }
     __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
     {
-        arg.GetHessian(multJ*dmult,J,-multH*value,H);
+		Sparse::HessianRow htmp;
+        arg.GetHessian(multJ,J,multH,htmp);
+		Sparse::HessianRow::MergeJacobianHessian(-value,J,J,dmult,htmp,H);
+		for(Sparse::Row::iterator it = J.Begin(); it != J.End(); ++it) it->second*=dmult;
     }
 	};
 
@@ -1073,7 +1076,11 @@ namespace INMOST
     }
     __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
     {
-        arg.GetHessian(multJ*dmult,J,-multH*value,H);
+        //arg.GetHessian(multJ*dmult,J,-multH*value,H);
+		Sparse::HessianRow htmp;
+        arg.GetHessian(multJ,J,multH,htmp);
+		Sparse::HessianRow::MergeJacobianHessian(-value,J,J,dmult,htmp,H);
+		for(Sparse::Row::iterator it = J.Begin(); it != J.End(); ++it) it->second*=dmult;
     }
 	};
 

Source/Headers/inmost_sparse.h

@@ -506,17 +506,10 @@ namespace INMOST
 		/// This class may be used to sum multiple sparse rows.
         /// \warning
         /// In parallel column indices of the matrix may span wider then
-        /// local row indices, to prevent any problem you are currently
-        /// advised to set total size of the matrix as interval of the
-        /// RowMerger. In future this may change, see todo 2 below.
-        /// \todo
-        /// 1. (testing!) Add iterators over entries.
-        /// 2. Implement multiple intervals for distributed computation,
-        ///    then in parallel the user may specify additional range of indexes
-        ///    for elements that lay on the borders between each pair of processors.
-        ///    Or even better implement mapping that will remap nonlocal entries to the
-        ///    end of current linked list when added and put them back in correct places
-        ///    when retrived. May use algorithm from class OrderInfo.
+        /// local row indices, to prevent any problem you can safely
+        /// set total size of the matrix as interval of the RowMerger.
+        /// For efficiency you should use RowMerger::SetNonlocal function
+		/// to provide information about non-local elements.
         class RowMerger
         {
         public:

Tests/CMakeLists.txt

+if(USE_AUTODIFF)
+add_subdirectory(autodiff_test000)
+endif(USE_AUTODIFF)
+
 if(USE_SOLVER)
 add_subdirectory(solver_test000)
 #add_subdirectory(solver_test001)

Tests/autodiff_test000/CMakeLists.txt

+project(autodiff_test000)
+set(SOURCE main.cpp)
+
+add_executable(autodiff_test000 ${SOURCE})
+target_link_libraries(autodiff_test000 inmost)
+
+if(USE_MPI)
+  message("linking autodiff_test000 with MPI")
+  target_link_libraries(autodiff_test000 ${MPI_LIBRARIES}) 
+  if(MPI_LINK_FLAGS)
+    set_target_properties(autodiff_test000 PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
+  endif() 
+endif(USE_MPI)
+
+
+
+add_test(NAME autodiff_test000_hessian_sin                 COMMAND $<TARGET_FILE:autodiff_test000> 0)
+add_test(NAME autodiff_test000_hessian_sin_mixed           COMMAND $<TARGET_FILE:autodiff_test000> 1)
+add_test(NAME autodiff_test000_hessian_cos_mixed           COMMAND $<TARGET_FILE:autodiff_test000> 2)

Tests/autodiff_test000/main.cpp

+#include <cstdio>
+#include <cmath>
+#include "inmost.h"
+using namespace INMOST;
+
+int main(int argc,char ** argv)
+{
+	int test = 0;
+	if (argc > 1)  test = atoi(argv[1]);
+
+
+	if( test == 0 ) //check derivative and hessian of sin(x*x+y*y)
+	{
+		// (sin(x*x+y*y))'' = (2cos(x*x+y*y)(xdx+ydy))' =
+		// 2cos(x*x+y*y)(dxdx+dydy) - 4sin(x*x+y*y)(x*x dxdx+y*y dydy + 2x*y dxdy)(xdx+ydy)
+		//
+		// Expected Jacobian
+		// 0: 2cos(x*x+y*y)x dx
+		// 1: 2cos(x*x+y*y)y dy
+		//
+		// Expected Hessian
+		// 0,0: 2*cos(x*x+y*y)-4*sin(x*x+y*y) x*x dxdx
+		// 0,1: -8*sin(x*x+y*y) x*y dx dy
+		// 1,1: 2*cos(x*x+y*y)-4*sin(x*x+y*y) y*y dydy
+		double _x = 0.5;
+		double _y = 0.25;
+		double _dx = 2*cos(_x*_x+_y*_y)*_x;
+		double _dy = 2*cos(_x*_x+_y*_y)*_y;
+		double _dxdx = 2*cos(_x*_x+_y*_y)-4*sin(_x*_x+_y*_y)*_x*_x;
+		double _dydy = 2*cos(_x*_x+_y*_y)-4*sin(_x*_x+_y*_y)*_y*_y;
+		double _dxdy = -8*sin(_x*_x+_y*_y)*_x*_y;
+		unknown x(_x,0), y(_y,1);
+		hessian_variable f = sin(x*x+y*y);
+		double dx = f.GetRow()[0];
+		double dy = f.GetRow()[1];
+		double dxdx = f.GetHessianRow()[Sparse::HessianRow::make_index(0,0)];
+		double dxdy = f.GetHessianRow()[Sparse::HessianRow::make_index(0,1)];
+		double dydy = f.GetHessianRow()[Sparse::HessianRow::make_index(1,1)];
+
+		bool error = false;
+		std::cout << std::setw(10) << "derivative " << std::setw(10) << "original " << std::setw(10) << "computed" << std::endl;
+		std::cout << std::setw(10) << "dx " << std::setw(10) << _dx << std::setw(10) << dx << std::endl;
+		std::cout << std::setw(10) << "dy " << std::setw(10) << _dy << std::setw(10) << dy << std::endl;
+		std::cout << std::setw(10) << "dxdx " << std::setw(10) << _dxdx << std::setw(10) << dxdx << std::endl;
+		std::cout << std::setw(10) << "dxdy " << std::setw(10) << _dxdy << std::setw(10) << dxdy << std::endl;
+		std::cout << std::setw(10) << "dydy " << std::setw(10) << _dydy << std::setw(10) << dydy << std::endl;
+		if( std::abs(dx-_dx) > 1.0e-9 ) error = true, std::cout << "Error in dx: " << std::abs(dx-_dx) << " original " << _dx << " computed " << dx << std::endl;
+		if( std::abs(dy-_dy) > 1.0e-9 ) error = true, std::cout << "Error in dy: " << std::abs(dy-_dy) << " original " << _dy << " computed " << dy << std::endl;
+		if( std::abs(dxdx-_dxdx) > 1.0e-9 ) error = true, std::cout << "Error in dxdx: " << std::abs(dxdx-_dxdx) << " original " << _dxdx << " computed " << dxdx << std::endl;
+		if( std::abs(dxdy-_dxdy) > 1.0e-9 ) error = true, std::cout << "Error in dxdy: " << std::abs(dxdy-_dxdy) << " original " << _dxdy << " computed " << dxdy << std::endl;
+		if( std::abs(dydy-_dydy) > 1.0e-9 ) error = true, std::cout << "Error in dydy: " << std::abs(dydy-_dydy) << " original " << _dydy << " computed " << dydy << std::endl;
+
+		if( error ) return -1;
+	}
+	else if( test == 1 )
+	{
+		// (sin(x*x+y*y+x*y))'' = (cos(x*x+y*y+x*y)(2xdx+2ydy+xdy+ydx))' =
+		// 2cos(x*x+y*y+x*y)(dxdx+dydy+dxdy) - sin(x*x+y*y+x*y)((2x+y)dx+(2y+x)dy)*((2x+y)dx+(2y+x)dy) =
+		// 2cos(x*x+y*y+x*y)(dxdx+dydy+dxdy) - sin(x*x+y*y+x*y)((2x+y)(2x+y)dxdx + 2(2x+y)(2y+x)dxdy + (2y+x)(2y+x)dydy)
+		//
+		// Expected Jacobian
+		// 0: cos(x*x+y*y+x*y)(2x+y) dx
+		// 1: cos(x*x+y*y+x*y)(2y+x) dy
+		//
+		// Expected Hessian
+		// 0,0: 2cos(x*x+y*y+x*y) -  sin(x*x+y*y+x*y) (2x+y)(2x+y) dxdx
+		// 0,1: 2cos(x*x+y*y+x*y) - 2sin(x*x+y*y+x*y) (2x+y)(2y+x) dxdy
+		// 1,1: 2cos(x*x+y*y+x*y) -  sin(x*x+y*y+x*y) (x+2y)(x+2y) dydy
+		double _x = 0.5;
+		double _y = 0.25;
+		double _dx = cos(_x*_x+_y*_y+_x*_y)*(2*_x+_y);
+		double _dy = cos(_x*_x+_y*_y+_x*_y)*(_x+2*_y);
+		double _dxdx = 2*cos(_x*_x+_y*_y+_x*_y)-sin(_x*_x+_y*_y+_x*_y)*(2*_x+_y)*(2*_x+_y);
+		double _dydy = 2*cos(_x*_x+_y*_y+_x*_y)-sin(_x*_x+_y*_y+_x*_y)*(_x+2*_y)*(_x+2*_y);
+		double _dxdy = 2*cos(_x*_x+_y*_y+_x*_y)-2*sin(_x*_x+_y*_y+_x*_y)*(2*_x+_y)*(_x+2*_y);
+		unknown x(_x,0), y(_y,1);
+		hessian_variable f = sin(x*x+y*y+x*y);
+		double dx = f.GetRow()[0];
+		double dy = f.GetRow()[1];
+		double dxdx = f.GetHessianRow()[Sparse::HessianRow::make_index(0,0)];
+		double dxdy = f.GetHessianRow()[Sparse::HessianRow::make_index(0,1)];
+		double dydy = f.GetHessianRow()[Sparse::HessianRow::make_index(1,1)];
+
+		bool error = false;
+		std::cout << std::setw(10) << "derivative " << std::setw(10) << "original " << std::setw(10) << "computed" << std::endl;
+		std::cout << std::setw(10) << "dx " << std::setw(10) << _dx << std::setw(10) << dx << std::endl;
+		std::cout << std::setw(10) << "dy " << std::setw(10) << _dy << std::setw(10) << dy << std::endl;
+		std::cout << std::setw(10) << "dxdx " << std::setw(10) << _dxdx << std::setw(10) << dxdx << std::endl;
+		std::cout << std::setw(10) << "dxdy " << std::setw(10) << _dxdy << std::setw(10) << dxdy << std::endl;
+		std::cout << std::setw(10) << "dydy " << std::setw(10) << _dydy << std::setw(10) << dydy << std::endl;
+		if( std::abs(dx-_dx) > 1.0e-9 ) error = true, std::cout << "Error in dx: " << std::abs(dx-_dx) << " original " << _dx << " computed " << dx << std::endl;
+		if( std::abs(dy-_dy) > 1.0e-9 ) error = true, std::cout << "Error in dy: " << std::abs(dy-_dy) << " original " << _dy << " computed " << dy << std::endl;
+		if( std::abs(dxdx-_dxdx) > 1.0e-9 ) error = true, std::cout << "Error in dxdx: " << std::abs(dxdx-_dxdx) << " original " << _dxdx << " computed " << dxdx << std::endl;
+		if( std::abs(dxdy-_dxdy) > 1.0e-9 ) error = true, std::cout << "Error in dxdy: " << std::abs(dxdy-_dxdy) << " original " << _dxdy << " computed " << dxdy << std::endl;
+		if( std::abs(dydy-_dydy) > 1.0e-9 ) error = true, std::cout << "Error in dydy: " << std::abs(dydy-_dydy) << " original " << _dydy << " computed " << dydy << std::endl;
+
+		if( error ) return -1;
+	}
+	else if( test == 2 )
+	{
+		// (cos(x*x+y*y+x*y))'' = (-sin(x*x+y*y+x*y)(2xdx+2ydy+xdy+ydx))' =
+		// -2sin(x*x+y*y+x*y)(dxdx+dydy+dxdy) - cos(x*x+y*y+x*y)((2x+y)dx+(2y+x)dy)*((2x+y)dx+(2y+x)dy) =
+		// -2sin(x*x+y*y+x*y)(dxdx+dydy+dxdy) - cos(x*x+y*y+x*y)((2x+y)(2x+y)dxdx + 2(2x+y)(2y+x)dxdy + (2y+x)(2y+x)dydy)
+		//
+		// Expected Jacobian
+		// 0: -sin(x*x+y*y+x*y)(2x+y) dx
+		// 1: -sin(x*x+y*y+x*y)(2y+x) dy
+		//
+		// Expected Hessian
+		// 0,0: -2sin(x*x+y*y+x*y) -  cos(x*x+y*y+x*y) (2x+y)(2x+y) dxdx
+		// 0,1: -2sin(x*x+y*y+x*y) - 2cos(x*x+y*y+x*y) (2x+y)(2y+x) dxdy
+		// 1,1: -2sin(x*x+y*y+x*y) -  cos(x*x+y*y+x*y) (x+2y)(x+2y) dydy
+		double _x = 0.5;
+		double _y = 0.25;
+		double _dx = -sin(_x*_x+_y*_y+_x*_y)*(2*_x+_y);
+		double _dy = -sin(_x*_x+_y*_y+_x*_y)*(_x+2*_y);
+		double _dxdx = -2*sin(_x*_x+_y*_y+_x*_y)-cos(_x*_x+_y*_y+_x*_y)*(2*_x+_y)*(2*_x+_y);
+		double _dydy = -2*sin(_x*_x+_y*_y+_x*_y)-cos(_x*_x+_y*_y+_x*_y)*(_x+2*_y)*(_x+2*_y);
+		double _dxdy = -2*sin(_x*_x+_y*_y+_x*_y)-2*cos(_x*_x+_y*_y+_x*_y)*(2*_x+_y)*(_x+2*_y);
+		unknown x(_x,0), y(_y,1);
+		hessian_variable f = cos(x*x+y*y+x*y);
+		double dx = f.GetRow()[0];
+		double dy = f.GetRow()[1];
+		double dxdx = f.GetHessianRow()[Sparse::HessianRow::make_index(0,0)];
+		double dxdy = f.GetHessianRow()[Sparse::HessianRow::make_index(0,1)];
+		double dydy = f.GetHessianRow()[Sparse::HessianRow::make_index(1,1)];
+
+		bool error = false;
+		std::cout << std::setw(10) << "derivative " << std::setw(10) << "original " << std::setw(10) << "computed" << std::endl;
+		std::cout << std::setw(10) << "dx " << std::setw(10) << _dx << std::setw(10) << dx << std::endl;
+		std::cout << std::setw(10) << "dy " << std::setw(10) << _dy << std::setw(10) << dy << std::endl;
+		std::cout << std::setw(10) << "dxdx " << std::setw(10) << _dxdx << std::setw(10) << dxdx << std::endl;
+		std::cout << std::setw(10) << "dxdy " << std::setw(10) << _dxdy << std::setw(10) << dxdy << std::endl;
+		std::cout << std::setw(10) << "dydy " << std::setw(10) << _dydy << std::setw(10) << dydy << std::endl;
+		if( std::abs(dx-_dx) > 1.0e-9 ) error = true, std::cout << "Error in dx: " << std::abs(dx-_dx) << " original " << _dx << " computed " << dx << std::endl;
+		if( std::abs(dy-_dy) > 1.0e-9 ) error = true, std::cout << "Error in dy: " << std::abs(dy-_dy) << " original " << _dy << " computed " << dy << std::endl;
+		if( std::abs(dxdx-_dxdx) > 1.0e-9 ) error = true, std::cout << "Error in dxdx: " << std::abs(dxdx-_dxdx) << " original " << _dxdx << " computed " << dxdx << std::endl;
+		if( std::abs(dxdy-_dxdy) > 1.0e-9 ) error = true, std::cout << "Error in dxdy: " << std::abs(dxdy-_dxdy) << " original " << _dxdy << " computed " << dxdy << std::endl;
+		if( std::abs(dydy-_dydy) > 1.0e-9 ) error = true, std::cout << "Error in dydy: " << std::abs(dydy-_dydy) << " original " << _dydy << " computed " << dydy << std::endl;
+
+		if( error ) return -1;
+	}
+	else if( test == 3 )
+	{
+		// (sqrt(x*x+y*y+x*y))'
+	}
+
+	return 0;
+}