diff --git a/LBFG/LBFGSpp/LineSearchNocedalWright.h b/LBFG/LBFGSpp/LineSearchNocedalWright.h
index 3248ff7bc5ab7002abb6f953214272529a86b1a1..b41b99cff5a33fea67fa3154ea806772203c434b 100644
--- a/LBFG/LBFGSpp/LineSearchNocedalWright.h
+++ b/LBFG/LBFGSpp/LineSearchNocedalWright.h
@@ -84,7 +84,7 @@ public:
test_curv = -param.wolfe * dg_init; // Curvature
// Ends of the line search range (step_lo > step_hi is allowed)
- Scalar step_hi, fx_hi, dg_hi;
+ Scalar step_hi, fx_hi;//, dg_hi;
Scalar step_lo = Scalar(0), fx_lo = fx_init, dg_lo = dg_init;
// STEP 1: Bracketing Phase
@@ -110,7 +110,7 @@ public:
{
step_hi = step;
fx_hi = fx;
- dg_hi = dg;
+ //dg_hi = dg;
break;
}
// If reaching here, then the sufficient decrease condition is satisfied
@@ -121,7 +121,7 @@ public:
step_hi = step_lo;
fx_hi = fx_lo;
- dg_hi = dg_lo;
+ //dg_hi = dg_lo;
step_lo = step;
fx_lo = fx;
dg_lo = dg;
@@ -188,7 +188,7 @@ public:
step_hi = step;
fx_hi = fx;
- dg_hi = dg;
+ //dg_hi = dg;
}
else
{
@@ -200,7 +200,7 @@ public:
{
step_hi = step_lo;
fx_hi = fx_lo;
- dg_hi = dg_lo;
+ //dg_hi = dg_lo;
}
if (step == step_lo)
diff --git a/drt.cpp b/drt.cpp
index 023b91e19091a2a4679aab6725bcd5247a759338..b69fb4731df7d8806fa8cc5e7ecbf383b3b71b60 100644
--- a/drt.cpp
+++ b/drt.cpp
@@ -3,14 +3,11 @@
#include <Eigen/Core>
#include <eisgenerator/eistype.h>
-
#include "tensoroptions.h"
#include "eigentorchconversions.h"
#include "eistotorch.h"
#include "LBFG/LBFGS.h"
-
-
static torch::Tensor guesStartingPoint(torch::Tensor& omega, torch::Tensor& impedanceSpectra)
{
std::vector<int64_t> size = omega.sizes().vec();
@@ -94,14 +91,6 @@ public:
int64_t size = x.numel();
torch::Tensor xLeft = x.narrow(0, 0, x.numel()-1);
- std::cout<<"x:\n"<<x<<'\n';
- std::cout<<"xLeft:\n"<<xLeft<<'\n';
- std::cout<<"real:\n"<<torch::real(impedanceSpectra)<<'\n';
- std::cout<<"imag:\n"<<torch::imag(impedanceSpectra)<<'\n';
- std::cout<<"aMatrixReal:\n"<<aMatrixReal<<'\n';
- std::cout<<"aMatrixImag:\n"<<aMatrixImag<<'\n';
-
-
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));
@@ -138,30 +127,31 @@ public:
}
};
-torch::Tensor calcDrt(torch::Tensor& impedanceSpectra, torch::Tensor& omegaTensor, FitMetics& fm)
+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;
- fitParam.epsilon = 1e-2;
- fitParam.max_iterations = 100;
- fitParam.max_linesearch = 1000;
+ fitParam.epsilon = fp.epsilon;
+ fitParam.max_iterations = fp.maxIter;
+ fitParam.max_linesearch = fp.maxIter*10;
LBFGSpp::LBFGSSolver<fvalue> solver(fitParam);
- RtFunct funct(impedanceSpectra, aMatrixImag, aMatrixReal, 0.01, 0.001);
+ RtFunct funct(impedanceSpectra, aMatrixImag, aMatrixReal, 0.01, fp.step);
torch::Tensor startTensor = guesStartingPoint(omegaTensor, impedanceSpectra);
Eigen::VectorX<fvalue> x = libtorch2eigenVector<fvalue>(startTensor);
fm.iterations = solver.minimize(funct, x, fm.fx);
- return eigenVector2libtorch<fvalue>(x);
+ torch::Tensor xT = eigenVector2libtorch<fvalue>(x);
+ return xT;
}
-torch::Tensor calcDrt(const std::vector<eis::DataPoint>& data, const std::vector<fvalue>& omegaVector, FitMetics& fm)
+torch::Tensor 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);
+ return calcDrt(impedanceSpectra, omegaTensor, fm, fp);
}
diff --git a/drt/drt.h b/drt/drt.h
index 192ee3e2cf750772b93406296b6ffee505f99d35..e3b9c6a22648eb8c0a7b961ad287d46e7b36467d 100644
--- a/drt/drt.h
+++ b/drt/drt.h
@@ -6,8 +6,17 @@ struct FitMetics
{
int iterations;
fvalue fx;
+ bool compleated;
};
-torch::Tensor calcDrt(torch::Tensor& impedanceSpectra, torch::Tensor& omegaTensor, FitMetics& fm);
+struct FitParameters
+{
+ int maxIter;
+ double epsilon;
+ double step;
+ 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);
-torch::Tensor calcDrt(const std::vector<eis::DataPoint>& data, const std::vector<fvalue>& omegaVector, FitMetics& fm);
+torch::Tensor calcDrt(const std::vector<eis::DataPoint>& data, const std::vector<fvalue>& omegaVector, FitMetics& fm, const FitParameters& fp);
diff --git a/eigentorchconversions.h b/eigentorchconversions.h
index efd7a488b17305fe88cb2048616d8563cd9a255b..2a3355c01825c630f3e3722730ce4321935e7fa7 100644
--- a/eigentorchconversions.h
+++ b/eigentorchconversions.h
@@ -51,7 +51,7 @@ torch::Tensor eigen2libtorch(MatrixXrm<V> &E, bool copydata = true)
template <typename V>
torch::Tensor eigenVector2libtorch(Eigen::Vector<V, Eigen::Dynamic> &E, bool copydata = true)
{
- std::vector<int64_t> dims = {E.cols()};
+ std::vector<int64_t> dims = {E.rows()};
auto T = torch::from_blob(E.data(), dims, tensorOptCpu<V>(false));
if (copydata)
return T.clone();
diff --git a/main.cpp b/main.cpp
index a17f0d3d7bd304cf509d3271a87c9da0d9a7e35a..d859e4a0c62a4f65e6f756e4cbe799aac3afa0a4 100644
--- a/main.cpp
+++ b/main.cpp
@@ -25,7 +25,10 @@ int main(int argc, char** argv)
{
std::cout<<std::scientific;
- eis::Range omega(1, 1e6, 4, true);
+ for(size_t i = 0; i < 10; ++i)
+ {
+
+ eis::Range omega(1, 1e6, 50, true);
std::vector<fvalue> omegaVector = omega.getRangeVector();
eis::Model model("r{10}-r{50}p{0.02, 0.8}");
@@ -33,11 +36,12 @@ int main(int argc, char** argv)
printImpedance(data);
FitMetics fm;
- torch::Tensor x = calcDrt(data, omegaVector, 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';
+ }
return 0;
}