diff --git a/CMakeLists.txt b/CMakeLists.txt index 12f3e86410133cef5fd5400bf34dcb89f8ed7a6e..eb4f9d1cf3d45d83a90ea9574034bcb9cf518f08 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,7 +6,6 @@ link_directories(${CMAKE_CURRENT_BINARY_DIR}) set (CMAKE_CXX_STANDARD 20) -find_package(Torch REQUIRED) find_package(Eigen3 REQUIRED) function(dump_variables) @@ -21,8 +20,8 @@ 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) +add_library(${PROJECT_NAME} SHARED drt.cpp) +target_link_libraries(${PROJECT_NAME} ${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) @@ -30,10 +29,10 @@ install(TARGETS ${PROJECT_NAME} DESTINATION lib) link_directories(${CMAKE_CURRENT_BINARY_DIR}) add_executable(${PROJECT_NAME}_test main.cpp) 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) +target_link_libraries(${PROJECT_NAME}_test -l${PROJECT_NAME} ${EIGEN3_LIBRARIES} eisgenerator) +target_include_directories(${PROJECT_NAME}_test PRIVATE . ${EIGEN3_INCLUDE_DIRS} ) +set_target_properties(${PROJECT_NAME}_test PROPERTIES COMPILE_FLAGS "-Wall -O2 -march=native -g" LINK_FLAGS "-flto") +install(TARGETS ${PROJECT_NAME}_test RUNTIME DESTINATION bin) set(API_HEADERS_DIR eisdrt/) set(API_HEADERS diff --git a/drt.cpp b/drt.cpp index c997d96c9b0bc07a84d39fe1e7246c1ab6aad91d..f2491dc43819cec094ed0e49b5779a215dcef7fd 100644 --- a/drt.cpp +++ b/drt.cpp @@ -1,65 +1,60 @@ #include "eisdrt/drt.h" -#include <ATen/ops/ones.h> -#include <ATen/ops/zeros.h> #include <Eigen/Core> +#include <Eigen/StdVector> #include <eisgenerator/eistype.h> +#include <iostream> -#include "tensoroptions.h" -#include "eigentorchconversions.h" -#include "eistotorch.h" +#include "Eigen/src/Core/Matrix.h" +#include "eistoeigen.h" #include "LBFG/LBFGSB.h" -static torch::Tensor guesStartingPoint(torch::Tensor& omega, torch::Tensor& impedanceSpectra) +static Eigen::Vector<fvalue, Eigen::Dynamic> guesStartingPoint(Eigen::Vector<fvalue, Eigen::Dynamic>& omega, Eigen::Vector<std::complex<fvalue>, Eigen::Dynamic>& impedanceSpectra) { - std::vector<int64_t> size = omega.sizes().vec(); - ++size[0]; - torch::Tensor startingPoint = torch::zeros(size, tensorOptCpu<fvalue>(false)); - startingPoint[-1] = torch::abs(impedanceSpectra[-1]); + Eigen::Vector<fvalue, Eigen::Dynamic> startingPoint = Eigen::Vector<fvalue, Eigen::Dynamic>::Zero(omega.size()+1); + startingPoint[startingPoint.size()-1] = std::abs(impedanceSpectra[impedanceSpectra.size()-1]); return startingPoint; } -static torch::Tensor aImag(torch::Tensor& omega) +static Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> aImag(Eigen::Vector<fvalue, Eigen::Dynamic>& omega) { - torch::Tensor tau = 1.0/(omega/(2*M_PI)); - torch::Tensor out = torch::zeros({omega.numel(), omega.numel()}, tensorOptCpu<fvalue>()); - auto outAccessor = out.accessor<float, 2>(); - auto omegaAccessor = omega.accessor<float, 1>(); - auto tauAccessor = tau.accessor<float, 1>(); - for(int32_t i = 0; i < out.size(0); ++i) + Eigen::Vector<fvalue, Eigen::Dynamic> tau = (omega * 1/static_cast<fvalue>(2*M_PI)).cwiseInverse(); + Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> out = + Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic>::Zero(omega.size(), omega.size()); + + for(int32_t i = 0; i < out.cols(); ++i) { - for(int32_t j = 0; j < out.size(1); ++j) + for(int32_t j = 0; j < out.rows(); ++j) { - outAccessor[i][j] = 0.5*(omegaAccessor[i]*tauAccessor[j])/(1+std::pow(omegaAccessor[i]*tauAccessor[j], 2)); + out(i,j) = 0.5*(omega[i]*tau[j])/(1+std::pow(omega[i]*tau[j], 2)); if(j == 0) - outAccessor[i][j] = outAccessor[i][j]*std::log(tauAccessor[j+1]/tauAccessor[j]); - else if(j == out.size(1)-1) - outAccessor[i][j] = outAccessor[i][j]*std::log(tauAccessor[j]/tauAccessor[j-1]); + out(i,j) = out(i,j)*std::log(tau[j+1]/tau[j]); + else if(j == out.rows()-1) + out(i,j) = out(i,j)*std::log(tau[j]/tau[j-1]); else - outAccessor[i][j] = outAccessor[i][j]*std::log(tauAccessor[j+1]/tauAccessor[j-1]); + out(i,j) = out(i,j)*std::log(tau[j+1]/tau[j-1]); } } return out; } -static torch::Tensor aReal(torch::Tensor& omega) +static Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> aReal(Eigen::Vector<fvalue, Eigen::Dynamic>& omega) { - torch::Tensor tau = 1.0/(omega/(2*M_PI)); - torch::Tensor out = torch::zeros({omega.numel(), omega.numel()}, torch::TensorOptions().dtype(torch::kFloat32)); - auto outAccessor = out.accessor<float, 2>(); - auto omegaAccessor = omega.accessor<float, 1>(); - auto tauAccessor = tau.accessor<float, 1>(); - for(int32_t i = 0; i < out.size(0); ++i) + Eigen::Vector<fvalue, Eigen::Dynamic> tau = (omega * 1/static_cast<fvalue>(2*M_PI)).cwiseInverse(); + Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> out = + Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic>::Zero(omega.size(), omega.size()); + + for(int32_t i = 0; i < out.cols(); ++i) { - for(int32_t j = 0; j < out.size(1); ++j) + for(int32_t j = 0; j < out.rows(); ++j) { - outAccessor[i][j] = -0.5/(1+std::pow(omegaAccessor[i]*tauAccessor[j], 2)); + out(i, j) = -0.5/(1+std::pow(omega[i]*tau[j], 2)); if(j == 0) - outAccessor[i][j] = outAccessor[i][j]*std::log(tauAccessor[j+1]/tauAccessor[j]); - else if(j == out.size(1)-1) - outAccessor[i][j] = outAccessor[i][j]*std::log(tauAccessor[j]/tauAccessor[j-1]); + out(i, j) = out(i, j)*std::log(tau[j+1]/tau[j]); + else if(j == out.rows()-1) + out(i, j) = out(i, j)*std::log(tau[j]/tau[j-1]); else - outAccessor[i][j] = outAccessor[i][j]*std::log(tauAccessor[j+1]/tauAccessor[j-1]); + out(i, j) = out(i, j)*std::log(tau[j+1]/tau[j-1]); } } return out; @@ -68,14 +63,17 @@ static torch::Tensor aReal(torch::Tensor& omega) class RtFunct { private: - torch::Tensor impedanceSpectra; - torch::Tensor aMatrixImag; - torch::Tensor aMatrixReal; + Eigen::Vector<std::complex<fvalue>, Eigen::Dynamic> impedanceSpectra; + Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> aMatrixImag; + Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> aMatrixReal; fvalue el; fvalue epsilon; public: - RtFunct(torch::Tensor impedanceSpectraI, torch::Tensor aMatrixImagI, torch::Tensor aMatrixRealI, fvalue elI, fvalue epsilonI): + RtFunct(Eigen::Vector<std::complex<fvalue>, Eigen::Dynamic> impedanceSpectraI, + Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> aMatrixImagI, + Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> aMatrixRealI, + fvalue elI, fvalue epsilonI): impedanceSpectra(impedanceSpectraI), aMatrixImag(aMatrixImagI), aMatrixReal(aMatrixRealI), @@ -85,59 +83,68 @@ public: } - fvalue function(const torch::Tensor& x) + fvalue function(const Eigen::Vector<fvalue, Eigen::Dynamic>& x) { - 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::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>(); + int64_t size = x.size(); + Eigen::Vector<fvalue, Eigen::Dynamic> xLeft = x.head(x.size()-1); + + std::cout<<"aMatrixReal:\n"<<aMatrixReal<<"\nxLeft:\n"<<xLeft<<"\nx:\n"<<x<<std::endl; + Eigen::Vector<fvalue, Eigen::Dynamic> t = aMatrixReal*xLeft; + std::cout<<"T1:\n"<<t<<std::endl; + t = t - impedanceSpectra.real(); + std::cout<<"T2:\n"<<t<<std::endl; + t = t.array() + x[size-1]; + std::cout<<"T3:\n"<<t<<std::endl; + t = t.array().pow(2); + std::cout<<"T4:\n"<<t<<std::endl; + fvalue MSE_re = t.sum(); + std::cout<<"T5:\n"<<MSE_re<<std::endl; + + t = (aMatrixImag*xLeft - impedanceSpectra.imag()).array().pow(2); + fvalue MSE_im = t.sum(); + fvalue reg_term = el/2*xLeft.array().pow(2).sum(); + fvalue obj = MSE_re + MSE_im + reg_term; + return obj; } - static torch::Tensor getGrad(std::function<fvalue(const torch::Tensor& x)> fn, const torch::Tensor& xTensor, fvalue epsilon) + static Eigen::Vector<fvalue, Eigen::Dynamic> getGrad(std::function<fvalue(const Eigen::Vector<fvalue, Eigen::Dynamic>& x)> fn, + Eigen::Vector<fvalue, Eigen::Dynamic>& x, fvalue epsilon) { - torch::Tensor out = torch::zeros(xTensor.sizes(), tensorOptCpu<fvalue>(false)); - auto outAccessor = out.accessor<fvalue, 1>(); - assert(checkTorchType<fvalue>(xTensor)); - auto xAccessor = xTensor.accessor<fvalue, 1>(); - for(int64_t i = 0; i < out.size(0); ++i) + Eigen::Vector<fvalue, Eigen::Dynamic> out = Eigen::Vector<fvalue, Eigen::Dynamic>::Zero(x.size()); + for(int64_t i = 0; i < out.size(); ++i) { - xAccessor[i] -= epsilon; - fvalue left = fn(xTensor); - xAccessor[i] += 2*epsilon; - fvalue right = fn(xTensor); - xAccessor[i] -= epsilon; - outAccessor[i] = (right-left)/(2*epsilon); + x[i] -= epsilon; + fvalue left = fn(x); + x[i] += 2*epsilon; + fvalue right = fn(x); + x[i] -= epsilon; + x[i] = (right-left)/(2*epsilon); } return out; } - fvalue operator()(const Eigen::VectorX<fvalue>& x, Eigen::VectorX<fvalue>& grad) + fvalue operator()(Eigen::VectorX<fvalue>& x, Eigen::VectorX<fvalue>& grad) { - Eigen::MatrixX<fvalue> xMatrix = x; - torch::Tensor xTensor = eigen2libtorch(xMatrix); - xTensor = xTensor.reshape({xTensor.numel()}); - torch::Tensor gradTensor = getGrad(std::bind(&RtFunct::function, this, std::placeholders::_1), xTensor, epsilon); - grad = libtorch2eigenVector<fvalue>(gradTensor); - return function(xTensor); + grad = getGrad(std::bind(&RtFunct::function, this, std::placeholders::_1), x, epsilon); + return function(x); } }; -static torch::Tensor calcBounds(torch::Tensor& impedanceSpectra, torch::Tensor startTensor) +static Eigen::Matrix<fvalue, Eigen::Dynamic, 2> calcBounds(Eigen::VectorX<std::complex<fvalue>>& impedanceSpectra, Eigen::VectorX<fvalue> 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); + Eigen::VectorX<fvalue> lowerBounds = Eigen::VectorX<fvalue>::Zero(startTensor.size()); + Eigen::VectorX<fvalue> upperBounds = Eigen::VectorX<fvalue>::Ones(startTensor.size())*impedanceSpectra.cwiseAbs().maxCoeff(); + + Eigen::Matrix<fvalue, Eigen::Dynamic, 2> out(lowerBounds.size(), 2); + out.col(0) = lowerBounds; + out.col(1) = upperBounds; + return out; } -torch::Tensor calcDrt(torch::Tensor& impedanceSpectra, torch::Tensor& omegaTensor, FitMetics& fm, const FitParameters& fp) +Eigen::VectorX<fvalue> calcDrt(Eigen::VectorX<std::complex<fvalue>>& impedanceSpectra, Eigen::VectorX<fvalue>& omegaTensor, FitMetics& fm, const FitParameters& fp) { - torch::Tensor aMatrixImag = aImag(omegaTensor); - torch::Tensor aMatrixReal = aReal(omegaTensor); + Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> aMatrixImag = aImag(omegaTensor); + Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> aMatrixReal = aReal(omegaTensor); LBFGSpp::LBFGSBParam<fvalue> fitParam; fitParam.epsilon = fp.epsilon; @@ -147,33 +154,29 @@ torch::Tensor calcDrt(torch::Tensor& impedanceSpectra, torch::Tensor& omegaTenso 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); + Eigen::VectorX<fvalue> x = guesStartingPoint(omegaTensor, impedanceSpectra); + std::cout<<"StartingPoint\n"<<x<<std::endl; + Eigen::Matrix<fvalue, Eigen::Dynamic, 2> bounds = calcBounds(impedanceSpectra, x); + Eigen::VectorX<fvalue> lowerBounds = bounds.col(0); + Eigen::VectorX<fvalue> upperBounds = bounds.col(1); - fm.iterations = solver.minimize(funct, x, fm.fx, lowerbound, upperbound); + fm.iterations = solver.minimize(funct, x, fm.fx, lowerBounds, upperBounds); - torch::Tensor xT = eigenVector2libtorch<fvalue>(x); - return xT; + return x; } -torch::Tensor calcDrt(const std::vector<eis::DataPoint>& data, const std::vector<fvalue>& omegaVector, FitMetics& fm, const FitParameters& fp) +Eigen::VectorX<fvalue> calcDrt(const std::vector<eis::DataPoint>& data, const std::vector<fvalue>& omegaVector, FitMetics& fm, const FitParameters& fp) { - torch::Tensor impedanceSpectra = eisToComplexTensor(data, nullptr); - torch::Tensor omegaTensor = fvalueVectorToTensor(const_cast<std::vector<fvalue>&>(omegaVector)).clone(); - return calcDrt(impedanceSpectra, omegaTensor, fm, fp); + Eigen::VectorX<std::complex<fvalue>> impedanceSpectra = eistoeigen(data); + Eigen::VectorX<fvalue> omega = Eigen::VectorX<fvalue>::Map(omegaVector.data(), omegaVector.size()); + return calcDrt(impedanceSpectra, omega, fm, fp); } -torch::Tensor calcDrt(const std::vector<eis::DataPoint>& data, FitMetics& fm, const FitParameters& fp) +Eigen::VectorX<fvalue> calcDrt(const std::vector<eis::DataPoint>& data, FitMetics& fm, const FitParameters& fp) { - torch::Tensor omegaTensor; - torch::Tensor impedanceSpectra = eisToComplexTensor(data, &omegaTensor); - return calcDrt(impedanceSpectra, omegaTensor, fm, fp); + Eigen::VectorX<fvalue> omega; + Eigen::VectorX<std::complex<fvalue>> impedanceSpectra = eistoeigen(data, &omega); + return calcDrt(impedanceSpectra, omega, fm, fp); } diff --git a/eigentorchconversions.h b/eigentorchconversions.h deleted file mode 100644 index 79b78fe58cc63044e1a82e33c66ff7c848db536e..0000000000000000000000000000000000000000 --- a/eigentorchconversions.h +++ /dev/null @@ -1,85 +0,0 @@ -#include <climits> -#include <sys/types.h> -#include <torch/torch.h> -#include <Eigen/Dense> -#include <torch/types.h> -#include <vector> - -#include "tensoroptions.h" - -template <typename V> -bool checkTorchType(const torch::Tensor& tensor) -{ - static_assert(std::is_same<V, float>::value || std::is_same<V, double>::value || - std::is_same<V, int64_t>::value || std::is_same<V, int32_t>::value || std::is_same<V, int8_t>::value, - "This function dose not work with this type"); - if constexpr(std::is_same<V, float>::value) - return tensor.dtype() == torch::kFloat32; - else if constexpr(std::is_same<V, double>::value) - return tensor.dtype() == torch::kFloat64; - else if constexpr(std::is_same<V, int64_t>::value) - return tensor.dtype() == torch::kInt64; - else if constexpr(std::is_same<V, int32_t>::value) - return tensor.dtype() == torch::kInt32; - else if constexpr(std::is_same<V, int8_t>::value) - return tensor.dtype() == torch::kInt8; -} - -template <typename V> -using MatrixXrm = typename Eigen::Matrix<V, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>; - -template <typename V> -torch::Tensor eigen2libtorch(Eigen::MatrixX<V> &M) -{ - Eigen::Matrix<V, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> E(M); - std::vector<int64_t> dims = {E.rows(), E.cols()}; - auto T = torch::from_blob(E.data(), dims, tensorOptCpu<V>(false)).clone(); - return T; -} - -template <typename V> -torch::Tensor eigen2libtorch(MatrixXrm<V> &E, bool copydata = true) -{ - std::vector<int64_t> dims = {E.rows(), E.cols()}; - auto T = torch::from_blob(E.data(), dims, tensorOptCpu<V>(false)); - if (copydata) - return T.clone(); - else - return T; -} - -template <typename V> -torch::Tensor eigenVector2libtorch(Eigen::Vector<V, Eigen::Dynamic> &E, bool copydata = true) -{ - std::vector<int64_t> dims = {E.rows()}; - auto T = torch::from_blob(E.data(), dims, tensorOptCpu<V>(false)); - if (copydata) - return T.clone(); - else - return T; -} - -template<typename V> -Eigen::Matrix<V, Eigen::Dynamic, Eigen::Dynamic> libtorch2eigenMaxtrix(torch::Tensor &Tin) -{ - /* - LibTorch is Row-major order and Eigen is Column-major order. - 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; -} - -template<typename V> -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/eisdrt/drt.h b/eisdrt/drt.h index 79f154b173a4b941bc6cb58782a279f4431b9a90..94a78f89d92b80abd4e806f5d9f5381fde92824e 100644 --- a/eisdrt/drt.h +++ b/eisdrt/drt.h @@ -1,6 +1,6 @@ #pragma once #include <eisgenerator/eistype.h> -#include <torch/torch.h> +#include <Eigen/Core> struct FitMetics { @@ -17,8 +17,8 @@ struct FitParameters FitParameters(int maxIterI, double epsilonI = 1e-2, double stepI = 0.001): maxIter(maxIterI), epsilon(epsilonI), step(stepI){} }; -torch::Tensor calcDrt(torch::Tensor& impedanceSpectra, torch::Tensor& omegaTensor, FitMetics& fm, const FitParameters& fp); +Eigen::VectorX<fvalue> calcDrt(Eigen::VectorX<std::complex<fvalue>>& impedanceSpectra, Eigen::VectorX<fvalue>& omegaTensor, FitMetics& fm, const FitParameters& fp); -torch::Tensor calcDrt(const std::vector<eis::DataPoint>& data, const std::vector<fvalue>& omegaVector, FitMetics& fm, const FitParameters& fp); +Eigen::VectorX<fvalue> calcDrt(const std::vector<eis::DataPoint>& data, const std::vector<fvalue>& omegaVector, FitMetics& fm, const FitParameters& fp); -torch::Tensor calcDrt(const std::vector<eis::DataPoint>& data, FitMetics& fm, const FitParameters& fp); +Eigen::VectorX<fvalue> calcDrt(const std::vector<eis::DataPoint>& data, FitMetics& fm, const FitParameters& fp); diff --git a/eistotorch.cpp b/eistotorch.cpp deleted file mode 100644 index 1d5b6dad4b411ae1ba2ce08d40df967a2c25412e..0000000000000000000000000000000000000000 --- a/eistotorch.cpp +++ /dev/null @@ -1,84 +0,0 @@ -#include "eistotorch.h" -#include <cassert> -#include <cmath> -#include <cstdint> - -#include "tensoroptions.h" - -torch::Tensor eisToComplexTensor(const std::vector<eis::DataPoint>& data, torch::Tensor* freqs) -{ - torch::Tensor output = torch::empty({static_cast<long int>(data.size())}, tensorOptCplxCpu<fvalue>()); - if(freqs) - *freqs = torch::empty({static_cast<long int>(data.size())}, tensorOptCpu<fvalue>()); - - torch::Tensor real = torch::real(output); - torch::Tensor imag = torch::imag(output); - - auto realAccessor = real.accessor<fvalue, 1>(); - auto imagAccessor = imag.accessor<fvalue, 1>(); - float* tensorFreqDataPtr = freqs ? freqs->contiguous().data_ptr<float>() : nullptr; - - for(size_t i = 0; i < data.size(); ++i) - { - fvalue real = data[i].im.real(); - fvalue imag = data[i].im.imag(); - if(std::isnan(real) || std::isinf(real)) - real = 0; - if(std::isnan(imag) || std::isinf(imag)) - real = 0; - - realAccessor[i] = real; - imagAccessor[i] = imag; - if(tensorFreqDataPtr) - tensorFreqDataPtr[i] = data[i % data.size()].omega; - } - - return output; -} - -torch::Tensor eisToTorch(const std::vector<eis::DataPoint>& data, torch::Tensor* freqs) -{ - torch::Tensor input = torch::empty({static_cast<long int>(data.size()*2)}, tensorOptCpu<fvalue>()); - if(freqs) - *freqs = torch::empty({static_cast<long int>(data.size()*2)}, tensorOptCpu<fvalue>()); - - float* tensorDataPtr = input.contiguous().data_ptr<float>(); - float* tensorFreqDataPtr = freqs ? freqs->contiguous().data_ptr<float>() : nullptr; - - for(size_t i = 0; i < data.size()*2; ++i) - { - float datapoint = i < data.size() ? data[i].im.real() : data[i - data.size()].im.imag(); - if(std::isnan(datapoint) || std::isinf(datapoint)) - datapoint = 0; - tensorDataPtr[i] = datapoint; - if(tensorFreqDataPtr) - tensorFreqDataPtr[i] = data[i % data.size()].omega; - } - - return input; -} - -std::vector<eis::DataPoint> torchToEis(const torch::Tensor& input) -{ - assert(input.numel() % 2 == 0); - input.reshape({1, input.numel()}); - - std::vector<eis::DataPoint> output(input.numel()/2); - - float* tensorDataPtr = input.contiguous().data_ptr<float>(); - - for(int64_t i = 0; i < input.numel()/2; ++i) - { - output[i].omega = i; - output[i].im.real(tensorDataPtr[i]); - output[i].im.imag(tensorDataPtr[i+input.numel()/2]); - } - - return output; -} - - -torch::Tensor fvalueVectorToTensor(std::vector<fvalue>& vect) -{ - return torch::from_blob(vect.data(), {static_cast<int64_t>(vect.size())}, tensorOptCpu<fvalue>()); -} diff --git a/eistotorch.h b/eistotorch.h deleted file mode 100644 index 2474d74f2ef873b98dc831f9c4e549da7ce48390..0000000000000000000000000000000000000000 --- a/eistotorch.h +++ /dev/null @@ -1,11 +0,0 @@ -#pragma once - -#include <eisgenerator/eistype.h> -#include <torch/torch.h> - -#include <vector> - -torch::Tensor eisToComplexTensor(const std::vector<eis::DataPoint>& data, torch::Tensor* freqs); -torch::Tensor eisToTorch(const std::vector<eis::DataPoint>& data, torch::Tensor* freqs = nullptr); -std::vector<eis::DataPoint> torchToEis(const torch::Tensor& input); -torch::Tensor fvalueVectorToTensor(std::vector<fvalue>& vect); diff --git a/main.cpp b/main.cpp index 8638717f99eaebb529530a632e0110715e386bef..4e2624532b06a7b9b6ed5f4ed91904d1ca688c76 100644 --- a/main.cpp +++ b/main.cpp @@ -3,6 +3,7 @@ #include <eisgenerator/eistype.h> #include "eisdrt/drt.h" +#include "eistoeigen.h" void printImpedance(const std::vector<eis::DataPoint>& data) { @@ -25,7 +26,7 @@ int main(int argc, char** argv) { std::cout<<std::scientific; - eis::Range omega(1, 1e6, 25, true); + eis::Range omega(1, 1e6, 3, true); std::vector<fvalue> omegaVector = omega.getRangeVector(); eis::Model model("r{10}-r{50}p{0.02, 0.8}"); @@ -34,8 +35,10 @@ int main(int argc, char** argv) std::vector<eis::DataPoint> data = model.executeSweep(omega); printImpedance(data); - FitMetics fm; - torch::Tensor x = calcDrt(data, omegaVector, fm, FitParameters(1000)); + FitMetics fm = {}; + Eigen::VectorX<fvalue> omega; + Eigen::VectorX<std::complex<fvalue>> impedanceSpectra = eistoeigen(data, &omega); + Eigen::VectorX<fvalue> x = calcDrt(impedanceSpectra, omega, fm, FitParameters(1000)); std::cout<<"Iterations: "<<fm.iterations<<'\n'; std::cout<<"fx "<<fm.fx<<'\n';