add bounds to optimization

CMakeLists.txt

 cmake_minimum_required(VERSION 3.19)
 
-project(drtwrap LANGUAGES CXX)
+project(eisdrt LANGUAGES CXX)
 
 link_directories(${CMAKE_CURRENT_BINARY_DIR})
 
@@ -17,10 +17,15 @@ function(dump_variables)
 	endforeach()
 endfunction()
 
+if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
+	set(CMAKE_INSTALL_PREFIX "/usr" CACHE PATH "..." FORCE)
+endif(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
+
 add_library(${PROJECT_NAME} SHARED eistotorch.cpp drt.cpp)
 target_link_libraries(${PROJECT_NAME} ${TORCH_LIBRARIES} ${EIGEN3_LIBRARIES} eisgenerator)
 target_include_directories(${PROJECT_NAME} PUBLIC ${TORCH_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ./LBFG)
 set_target_properties(${PROJECT_NAME} PROPERTIES COMPILE_FLAGS "-Wall -O2 -march=native -g" LINK_FLAGS "-flto")
+install(TARGETS ${PROJECT_NAME} DESTINATION lib)
 
 link_directories(${CMAKE_CURRENT_BINARY_DIR})
 add_executable(${PROJECT_NAME}_test main.cpp)
@@ -28,5 +33,10 @@ add_dependencies(${PROJECT_NAME}_test ${PROJECT_NAME})
 target_link_libraries(${PROJECT_NAME}_test -l${PROJECT_NAME} ${TORCH_LIBRARIES} eisgenerator)
 target_include_directories(${PROJECT_NAME}_test PRIVATE . ${TORCH_INCLUDE_DIRS})
 set_target_properties(${PROJECT_NAME} PROPERTIES COMPILE_FLAGS "-Wall -O2 -march=native -g" LINK_FLAGS "-flto")
-
 install(TARGETS ${PROJECT_NAME} RUNTIME DESTINATION bin)
+
+set(API_HEADERS_DIR drt/)
+set(API_HEADERS
+	${API_HEADERS_DIR}/drt.h
+)
+install(FILES ${API_HEADERS} DESTINATION include/${PROJECT_NAME})

drt.cpp

-#include "drt/drt.h"
+#include "eisdrt/drt.h"
 
+#include <ATen/ops/ones.h>
+#include <ATen/ops/zeros.h>
 #include <Eigen/Core>
 #include <eisgenerator/eistype.h>
 
 #include "tensoroptions.h"
 #include "eigentorchconversions.h"
 #include "eistotorch.h"
-#include "LBFG/LBFGS.h"
+#include "LBFG/LBFGSB.h"
 
 static torch::Tensor guesStartingPoint(torch::Tensor& omega, torch::Tensor& impedanceSpectra)
 {
@@ -85,15 +87,13 @@ public:
 
 	fvalue function(const torch::Tensor& x)
 	{
-		assert(checkTorchType<fvalue>(x));
-		assert(x.sizes().size() == 1);
 		auto xAccessor = x.accessor<fvalue, 1>();
 		int64_t size = x.numel();
 		torch::Tensor xLeft = x.narrow(0, 0, x.numel()-1);
 
-		torch::Tensor MSE_re = torch::sum(torch::pow(xAccessor[size-1] + torch::matmul(aMatrixReal, xLeft) - torch::real(impedanceSpectra), 2));
-		torch::Tensor MSE_im = torch::sum(torch::pow(torch::matmul(aMatrixImag, xLeft) - torch::imag(impedanceSpectra), 2));
-		torch::Tensor reg_term = el/2*torch::sum(torch::pow(xLeft, 2));
+		torch::Tensor MSE_re = torch::sum(torch::pow(xAccessor[size-1] + torch::matmul(aMatrixReal, xLeft) - torch::real(impedanceSpectra), 2), torch::typeMetaToScalarType(x.dtype()));
+		torch::Tensor MSE_im = torch::sum(torch::pow(torch::matmul(aMatrixImag, xLeft) - torch::imag(impedanceSpectra), 2), torch::typeMetaToScalarType(x.dtype()));
+		torch::Tensor reg_term = el/2*torch::sum(torch::pow(xLeft, 2), torch::typeMetaToScalarType(x.dtype()));
 		torch::Tensor obj = MSE_re + MSE_im + reg_term;
 		return obj.item().to<fvalue>();
 	}
@@ -127,23 +127,36 @@ public:
 	}
 };
 
+static torch::Tensor calcBounds(torch::Tensor& impedanceSpectra, torch::Tensor startTensor)
+{
+	torch::Tensor lowerBounds = torch::zeros({1, startTensor.numel()}, tensorOptCpu<fvalue>());
+	torch::Tensor upperBounds = torch::ones({1, startTensor.numel()}, tensorOptCpu<fvalue>())*torch::max(torch::abs(impedanceSpectra));
+	return torch::cat({lowerBounds, upperBounds}, 0);
+}
+
 torch::Tensor calcDrt(torch::Tensor& impedanceSpectra, torch::Tensor& omegaTensor, FitMetics& fm, const FitParameters& fp)
 {
 	torch::Tensor aMatrixImag = aImag(omegaTensor);
 	torch::Tensor aMatrixReal = aReal(omegaTensor);
 
-	LBFGSpp::LBFGSParam<fvalue> fitParam;
+	LBFGSpp::LBFGSBParam<fvalue> fitParam;
 	fitParam.epsilon = fp.epsilon;
 	fitParam.max_iterations = fp.maxIter;
 	fitParam.max_linesearch = fp.maxIter*10;
 
-	LBFGSpp::LBFGSSolver<fvalue> solver(fitParam);
+	LBFGSpp::LBFGSBSolver<fvalue> solver(fitParam);
 	RtFunct funct(impedanceSpectra, aMatrixImag, aMatrixReal, 0.01, fp.step);
 
 	torch::Tensor startTensor = guesStartingPoint(omegaTensor, impedanceSpectra);
+	torch::Tensor bounds = calcBounds(impedanceSpectra, startTensor);
+	torch::Tensor lowerBoundTensor = bounds.select(0, 0);
+	torch::Tensor upperBoundTensor = bounds.select(0, 1);
+	Eigen::VectorX<fvalue> lowerbound = libtorch2eigenVector<fvalue>(lowerBoundTensor);
+	Eigen::VectorX<fvalue> upperbound = libtorch2eigenVector<fvalue>(upperBoundTensor);
+
 	Eigen::VectorX<fvalue> x = libtorch2eigenVector<fvalue>(startTensor);
 
-	fm.iterations = solver.minimize(funct, x, fm.fx);
+	fm.iterations = solver.minimize(funct, x, fm.fx, lowerbound, upperbound);
 
 	torch::Tensor xT = eigenVector2libtorch<fvalue>(x);
 	return xT;

eigentorchconversions.h

@@ -67,6 +67,7 @@ Eigen::Matrix<V, Eigen::Dynamic, Eigen::Dynamic> libtorch2eigenMaxtrix(torch::Te
 	MatrixXrm uses Eigen::RowMajor for compatibility.
 	*/
 	assert(checkTorchType<V>(Tin));
+	Tin = Tin.contiguous();
 	auto T = Tin.to(torch::kCPU);
 	Eigen::Map<MatrixXrm<V>> E(T.data_ptr<V>(), T.size(0), T.size(1));
 	return E;
@@ -77,6 +78,7 @@ Eigen::Vector<V, Eigen::Dynamic> libtorch2eigenVector(torch::Tensor &Tin)
 {
 	assert(Tin.sizes().size() == 1);
 	assert(checkTorchType<V>(Tin));
+	Tin = Tin.contiguous();
 	auto T = Tin.to(torch::kCPU);
 	Eigen::Map<Eigen::Vector<V, Eigen::Dynamic>> E(T.data_ptr<V>(), T.numel());
 	return E;

main.cpp

@@ -2,7 +2,7 @@
 #include <eisgenerator/model.h>
 #include <eisgenerator/eistype.h>
 
-#include "drt/drt.h"
+#include "eisdrt/drt.h"
 
 void printImpedance(const std::vector<eis::DataPoint>& data)
 {
@@ -25,13 +25,12 @@ int main(int argc, char** argv)
 {
 	std::cout<<std::scientific;
 
-	for(size_t i = 0; i < 10; ++i)
-	{
-
-	eis::Range omega(1, 1e6, 50, true);
+	eis::Range omega(1, 1e6, 25, true);
 	std::vector<fvalue> omegaVector = omega.getRangeVector();
 	eis::Model model("r{10}-r{50}p{0.02, 0.8}");
 
+	for(size_t i = 0; i < 2; ++i)
+	{
 		std::vector<eis::DataPoint> data = model.executeSweep(omega);
 		printImpedance(data);