From a13bbd25c444ea9e3f2050604154b5fa97035871 Mon Sep 17 00:00:00 2001 From: Carl Philipp Klemm <philipp@uvos.xyz> Date: Fri, 12 May 2023 11:42:51 +0200 Subject: [PATCH] add bounds to optimization --- CMakeLists.txt | 14 ++++++++++++-- drt.cpp | 33 +++++++++++++++++++++++---------- eigentorchconversions.h | 2 ++ {drt => eisdrt}/drt.h | 0 main.cpp | 23 +++++++++++------------ 5 files changed, 48 insertions(+), 24 deletions(-) rename {drt => eisdrt}/drt.h (100%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 8b108cc..0a133f0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ 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}) diff --git a/drt.cpp b/drt.cpp index b69fb47..f7429ee 100644 --- a/drt.cpp +++ b/drt.cpp @@ -1,12 +1,14 @@ -#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; diff --git a/eigentorchconversions.h b/eigentorchconversions.h index 2a3355c..79b78fe 100644 --- a/eigentorchconversions.h +++ b/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; diff --git a/drt/drt.h b/eisdrt/drt.h similarity index 100% rename from drt/drt.h rename to eisdrt/drt.h diff --git a/main.cpp b/main.cpp index d859e4a..8638717 100644 --- a/main.cpp +++ b/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,22 +25,21 @@ 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}"); - std::vector<eis::DataPoint> data = model.executeSweep(omega); - printImpedance(data); + for(size_t i = 0; i < 2; ++i) + { + std::vector<eis::DataPoint> data = model.executeSweep(omega); + printImpedance(data); - FitMetics fm; - torch::Tensor x = calcDrt(data, omegaVector, fm, FitParameters(1000)); + FitMetics fm; + torch::Tensor x = calcDrt(data, omegaVector, fm, FitParameters(1000)); - std::cout<<"Iterations: "<<fm.iterations<<'\n'; - std::cout<<"fx "<<fm.fx<<'\n'; - std::cout<<"xVect\n"<<x<<'\n'; + std::cout<<"Iterations: "<<fm.iterations<<'\n'; + std::cout<<"fx "<<fm.fx<<'\n'; + std::cout<<"xVect\n"<<x<<'\n'; } return 0; -- GitLab