Skip to content
Snippets Groups Projects
Commit a8d6a61b authored by Carl Philipp Klemm's avatar Carl Philipp Klemm
Browse files

convert to pure eigen

parent 77d0fd61
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 6
− 7
View file @ a8d6a61b
...@@ -6,7 +6,6 @@ link_directories(${CMAKE_CURRENT_BINARY_DIR}) ...@@ -6,7 +6,6 @@ link_directories(${CMAKE_CURRENT_BINARY_DIR})
set (CMAKE_CXX_STANDARD 20) set (CMAKE_CXX_STANDARD 20)
find_package(Torch REQUIRED)
find_package(Eigen3 REQUIRED) find_package(Eigen3 REQUIRED)
function(dump_variables) function(dump_variables)
...@@ -21,8 +20,8 @@ if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) ...@@ -21,8 +20,8 @@ if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
set(CMAKE_INSTALL_PREFIX "/usr" CACHE PATH "..." FORCE) set(CMAKE_INSTALL_PREFIX "/usr" CACHE PATH "..." FORCE)
endif(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) endif(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
add_library(${PROJECT_NAME} SHARED eistotorch.cpp drt.cpp) add_library(${PROJECT_NAME} SHARED drt.cpp)
target_link_libraries(${PROJECT_NAME} ${TORCH_LIBRARIES} ${EIGEN3_LIBRARIES} eisgenerator) target_link_libraries(${PROJECT_NAME} ${EIGEN3_LIBRARIES} eisgenerator)
target_include_directories(${PROJECT_NAME} PUBLIC ${TORCH_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ./LBFG) 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") set_target_properties(${PROJECT_NAME} PROPERTIES COMPILE_FLAGS "-Wall -O2 -march=native -g" LINK_FLAGS "-flto")
install(TARGETS ${PROJECT_NAME} DESTINATION lib) install(TARGETS ${PROJECT_NAME} DESTINATION lib)
...@@ -30,10 +29,10 @@ install(TARGETS ${PROJECT_NAME} DESTINATION lib) ...@@ -30,10 +29,10 @@ install(TARGETS ${PROJECT_NAME} DESTINATION lib)
link_directories(${CMAKE_CURRENT_BINARY_DIR}) link_directories(${CMAKE_CURRENT_BINARY_DIR})
add_executable(${PROJECT_NAME}_test main.cpp) add_executable(${PROJECT_NAME}_test main.cpp)
add_dependencies(${PROJECT_NAME}_test ${PROJECT_NAME}) add_dependencies(${PROJECT_NAME}_test ${PROJECT_NAME})
target_link_libraries(${PROJECT_NAME}_test -l${PROJECT_NAME} ${TORCH_LIBRARIES} eisgenerator) target_link_libraries(${PROJECT_NAME}_test -l${PROJECT_NAME} ${EIGEN3_LIBRARIES} eisgenerator)
target_include_directories(${PROJECT_NAME}_test PRIVATE . ${TORCH_INCLUDE_DIRS}) target_include_directories(${PROJECT_NAME}_test PRIVATE . ${EIGEN3_INCLUDE_DIRS} )
set_target_properties(${PROJECT_NAME} PROPERTIES COMPILE_FLAGS "-Wall -O2 -march=native -g" LINK_FLAGS "-flto") set_target_properties(${PROJECT_NAME}_test PROPERTIES COMPILE_FLAGS "-Wall -O2 -march=native -g" LINK_FLAGS "-flto")
install(TARGETS ${PROJECT_NAME} RUNTIME DESTINATION bin) install(TARGETS ${PROJECT_NAME}_test RUNTIME DESTINATION bin)
set(API_HEADERS_DIR eisdrt/) set(API_HEADERS_DIR eisdrt/)
set(API_HEADERS set(API_HEADERS
......
drt.cpp
+ 98
− 95
View file @ a8d6a61b
#include "eisdrt/drt.h" #include "eisdrt/drt.h"
#include <ATen/ops/ones.h>
#include <ATen/ops/zeros.h>
#include <Eigen/Core> #include <Eigen/Core>
#include <Eigen/StdVector>
#include <eisgenerator/eistype.h> #include <eisgenerator/eistype.h>
#include <iostream>
#include "tensoroptions.h" #include "Eigen/src/Core/Matrix.h"
#include "eigentorchconversions.h" #include "eistoeigen.h"
#include "eistotorch.h"
#include "LBFG/LBFGSB.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(); Eigen::Vector<fvalue, Eigen::Dynamic> startingPoint = Eigen::Vector<fvalue, Eigen::Dynamic>::Zero(omega.size()+1);
++size[0]; startingPoint[startingPoint.size()-1] = std::abs(impedanceSpectra[impedanceSpectra.size()-1]);
torch::Tensor startingPoint = torch::zeros(size, tensorOptCpu<fvalue>(false));
startingPoint[-1] = torch::abs(impedanceSpectra[-1]);
return startingPoint; 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)); Eigen::Vector<fvalue, Eigen::Dynamic> tau = (omega * 1/static_cast<fvalue>(2*M_PI)).cwiseInverse();
torch::Tensor out = torch::zeros({omega.numel(), omega.numel()}, tensorOptCpu<fvalue>()); Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> out =
auto outAccessor = out.accessor<float, 2>(); Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic>::Zero(omega.size(), omega.size());
auto omegaAccessor = omega.accessor<float, 1>();
auto tauAccessor = tau.accessor<float, 1>(); for(int32_t i = 0; i < out.cols(); ++i)
for(int32_t i = 0; i < out.size(0); ++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) if(j == 0)
outAccessor[i][j] = outAccessor[i][j]*std::log(tauAccessor[j+1]/tauAccessor[j]); out(i,j) = out(i,j)*std::log(tau[j+1]/tau[j]);
else if(j == out.size(1)-1) else if(j == out.rows()-1)
outAccessor[i][j] = outAccessor[i][j]*std::log(tauAccessor[j]/tauAccessor[j-1]); out(i,j) = out(i,j)*std::log(tau[j]/tau[j-1]);
else 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; 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)); Eigen::Vector<fvalue, Eigen::Dynamic> tau = (omega * 1/static_cast<fvalue>(2*M_PI)).cwiseInverse();
torch::Tensor out = torch::zeros({omega.numel(), omega.numel()}, torch::TensorOptions().dtype(torch::kFloat32)); Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> out =
auto outAccessor = out.accessor<float, 2>(); Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic>::Zero(omega.size(), omega.size());
auto omegaAccessor = omega.accessor<float, 1>();
auto tauAccessor = tau.accessor<float, 1>(); for(int32_t i = 0; i < out.cols(); ++i)
for(int32_t i = 0; i < out.size(0); ++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) if(j == 0)
outAccessor[i][j] = outAccessor[i][j]*std::log(tauAccessor[j+1]/tauAccessor[j]); out(i, j) = out(i, j)*std::log(tau[j+1]/tau[j]);
else if(j == out.size(1)-1) else if(j == out.rows()-1)
outAccessor[i][j] = outAccessor[i][j]*std::log(tauAccessor[j]/tauAccessor[j-1]); out(i, j) = out(i, j)*std::log(tau[j]/tau[j-1]);
else 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; return out;
...@@ -68,14 +63,17 @@ static torch::Tensor aReal(torch::Tensor& omega) ...@@ -68,14 +63,17 @@ static torch::Tensor aReal(torch::Tensor& omega)
class RtFunct class RtFunct
{ {
private: private:
torch::Tensor impedanceSpectra; Eigen::Vector<std::complex<fvalue>, Eigen::Dynamic> impedanceSpectra;
torch::Tensor aMatrixImag; Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> aMatrixImag;
torch::Tensor aMatrixReal; Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> aMatrixReal;
fvalue el; fvalue el;
fvalue epsilon; fvalue epsilon;
public: 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), impedanceSpectra(impedanceSpectraI),
aMatrixImag(aMatrixImagI), aMatrixImag(aMatrixImagI),
aMatrixReal(aMatrixRealI), aMatrixReal(aMatrixRealI),
...@@ -85,59 +83,68 @@ public: ...@@ -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.size();
int64_t size = x.numel(); Eigen::Vector<fvalue, Eigen::Dynamic> xLeft = x.head(x.size()-1);
torch::Tensor xLeft = x.narrow(0, 0, x.numel()-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;
torch::Tensor MSE_re = torch::sum(torch::pow(xAccessor[size-1] + torch::matmul(aMatrixReal, xLeft) - torch::real(impedanceSpectra), 2), torch::typeMetaToScalarType(x.dtype())); t = (aMatrixImag*xLeft - impedanceSpectra.imag()).array().pow(2);
torch::Tensor MSE_im = torch::sum(torch::pow(torch::matmul(aMatrixImag, xLeft) - torch::imag(impedanceSpectra), 2), torch::typeMetaToScalarType(x.dtype())); fvalue MSE_im = t.sum();
torch::Tensor reg_term = el/2*torch::sum(torch::pow(xLeft, 2), torch::typeMetaToScalarType(x.dtype())); fvalue reg_term = el/2*xLeft.array().pow(2).sum();
torch::Tensor obj = MSE_re + MSE_im + reg_term; fvalue obj = MSE_re + MSE_im + reg_term;
return obj.item().to<fvalue>(); 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)); Eigen::Vector<fvalue, Eigen::Dynamic> out = Eigen::Vector<fvalue, Eigen::Dynamic>::Zero(x.size());
auto outAccessor = out.accessor<fvalue, 1>(); for(int64_t i = 0; i < out.size(); ++i)
assert(checkTorchType<fvalue>(xTensor));
auto xAccessor = xTensor.accessor<fvalue, 1>();
for(int64_t i = 0; i < out.size(0); ++i)
{ {
xAccessor[i] -= epsilon; x[i] -= epsilon;
fvalue left = fn(xTensor); fvalue left = fn(x);
xAccessor[i] += 2*epsilon; x[i] += 2*epsilon;
fvalue right = fn(xTensor); fvalue right = fn(x);
xAccessor[i] -= epsilon; x[i] -= epsilon;
outAccessor[i] = (right-left)/(2*epsilon); x[i] = (right-left)/(2*epsilon);
} }
return out; 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; grad = getGrad(std::bind(&RtFunct::function, this, std::placeholders::_1), x, epsilon);
torch::Tensor xTensor = eigen2libtorch(xMatrix); return function(x);
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);
} }
}; };
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>()); Eigen::VectorX<fvalue> lowerBounds = Eigen::VectorX<fvalue>::Zero(startTensor.size());
torch::Tensor upperBounds = torch::ones({1, startTensor.numel()}, tensorOptCpu<fvalue>())*torch::max(torch::abs(impedanceSpectra)); Eigen::VectorX<fvalue> upperBounds = Eigen::VectorX<fvalue>::Ones(startTensor.size())*impedanceSpectra.cwiseAbs().maxCoeff();
return torch::cat({lowerBounds, upperBounds}, 0);
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); Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> aMatrixImag = aImag(omegaTensor);
torch::Tensor aMatrixReal = aReal(omegaTensor); Eigen::Matrix<fvalue, Eigen::Dynamic, Eigen::Dynamic> aMatrixReal = aReal(omegaTensor);
LBFGSpp::LBFGSBParam<fvalue> fitParam; LBFGSpp::LBFGSBParam<fvalue> fitParam;
fitParam.epsilon = fp.epsilon; fitParam.epsilon = fp.epsilon;
...@@ -147,33 +154,29 @@ torch::Tensor calcDrt(torch::Tensor& impedanceSpectra, torch::Tensor& omegaTenso ...@@ -147,33 +154,29 @@ torch::Tensor calcDrt(torch::Tensor& impedanceSpectra, torch::Tensor& omegaTenso
LBFGSpp::LBFGSBSolver<fvalue> solver(fitParam); LBFGSpp::LBFGSBSolver<fvalue> solver(fitParam);
RtFunct funct(impedanceSpectra, aMatrixImag, aMatrixReal, 0.01, fp.step); RtFunct funct(impedanceSpectra, aMatrixImag, aMatrixReal, 0.01, fp.step);
torch::Tensor startTensor = guesStartingPoint(omegaTensor, impedanceSpectra); Eigen::VectorX<fvalue> x = guesStartingPoint(omegaTensor, impedanceSpectra);
torch::Tensor bounds = calcBounds(impedanceSpectra, startTensor); std::cout<<"StartingPoint\n"<<x<<std::endl;
torch::Tensor lowerBoundTensor = bounds.select(0, 0); Eigen::Matrix<fvalue, Eigen::Dynamic, 2> bounds = calcBounds(impedanceSpectra, x);
torch::Tensor upperBoundTensor = bounds.select(0, 1); Eigen::VectorX<fvalue> lowerBounds = bounds.col(0);
Eigen::VectorX<fvalue> lowerbound = libtorch2eigenVector<fvalue>(lowerBoundTensor); Eigen::VectorX<fvalue> upperBounds = bounds.col(1);
Eigen::VectorX<fvalue> upperbound = libtorch2eigenVector<fvalue>(upperBoundTensor);
Eigen::VectorX<fvalue> x = libtorch2eigenVector<fvalue>(startTensor);
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 x;
return xT;
} }
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); Eigen::VectorX<std::complex<fvalue>> impedanceSpectra = eistoeigen(data);
torch::Tensor omegaTensor = fvalueVectorToTensor(const_cast<std::vector<fvalue>&>(omegaVector)).clone(); Eigen::VectorX<fvalue> omega = Eigen::VectorX<fvalue>::Map(omegaVector.data(), omegaVector.size());
return calcDrt(impedanceSpectra, omegaTensor, fm, fp); 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; Eigen::VectorX<fvalue> omega;
torch::Tensor impedanceSpectra = eisToComplexTensor(data, &omegaTensor); Eigen::VectorX<std::complex<fvalue>> impedanceSpectra = eistoeigen(data, &omega);
return calcDrt(impedanceSpectra, omegaTensor, fm, fp); return calcDrt(impedanceSpectra, omega, fm, fp);
} }
eigentorchconversions.h deleted 100644 → 0
+ 0
− 85
View file @ 77d0fd61
#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;
}
eisdrt/drt.h
+ 4
− 4
View file @ a8d6a61b
#pragma once #pragma once
#include <eisgenerator/eistype.h> #include <eisgenerator/eistype.h>
#include <torch/torch.h> #include <Eigen/Core>
struct FitMetics struct FitMetics
{ {
...@@ -17,8 +17,8 @@ struct FitParameters ...@@ -17,8 +17,8 @@ struct FitParameters
FitParameters(int maxIterI, double epsilonI = 1e-2, double stepI = 0.001): maxIter(maxIterI), epsilon(epsilonI), step(stepI){} 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);
eistotorch.cpp deleted 100644 → 0
+ 0
− 84
View file @ 77d0fd61
#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>());
}
eistotorch.h deleted 100644 → 0
+ 0
− 11
View file @ 77d0fd61
#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);
main.cpp
+ 6
− 3
View file @ a8d6a61b
...@@ -3,6 +3,7 @@ ...@@ -3,6 +3,7 @@
#include <eisgenerator/eistype.h> #include <eisgenerator/eistype.h>
#include "eisdrt/drt.h" #include "eisdrt/drt.h"
#include "eistoeigen.h"
void printImpedance(const std::vector<eis::DataPoint>& data) void printImpedance(const std::vector<eis::DataPoint>& data)
{ {
...@@ -25,7 +26,7 @@ int main(int argc, char** argv) ...@@ -25,7 +26,7 @@ int main(int argc, char** argv)
{ {
std::cout<<std::scientific; std::cout<<std::scientific;
eis::Range omega(1, 1e6, 25, true); eis::Range omega(1, 1e6, 3, true);
std::vector<fvalue> omegaVector = omega.getRangeVector(); std::vector<fvalue> omegaVector = omega.getRangeVector();
eis::Model model("r{10}-r{50}p{0.02, 0.8}"); eis::Model model("r{10}-r{50}p{0.02, 0.8}");
...@@ -34,8 +35,10 @@ int main(int argc, char** argv) ...@@ -34,8 +35,10 @@ int main(int argc, char** argv)
std::vector<eis::DataPoint> data = model.executeSweep(omega); std::vector<eis::DataPoint> data = model.executeSweep(omega);
printImpedance(data); printImpedance(data);
FitMetics fm; FitMetics fm = {};
torch::Tensor x = calcDrt(data, omegaVector, fm, FitParameters(1000)); 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<<"Iterations: "<<fm.iterations<<'\n';
std::cout<<"fx "<<fm.fx<<'\n'; std::cout<<"fx "<<fm.fx<<'\n';
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment