diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9fec7abfcf24439f97b3e8a2747fb3e2f7a404a8..df8a4d769477f175ceaedf9369418b24851f07ac 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -27,6 +27,7 @@ set(SRC_FILES
translators.cpp
randomgen.cpp
compcache.cpp
+ linearregession.cpp
)
set(API_HEADERS_CPP_DIR eisgenerator/)
diff --git a/basicmath.cpp b/basicmath.cpp
index 84db4b396f9b9343d8d73ee5bc13dfdf2c3b0d0e..dda06ece53ba619863f11998215a12907a5b17e7 100644
--- a/basicmath.cpp
+++ b/basicmath.cpp
@@ -2,9 +2,11 @@
#include <algorithm>
#include <random>
#include <limits>
+#include <stdexcept>
#include "eistype.h"
#include "log.h"
+#include "linearregession.h"
static size_t gradIndex(size_t dataSize, size_t inputIndex)
{
@@ -259,6 +261,25 @@ getLrClosest(const eis::DataPoint& dp, std::vector<eis::DataPoint>::const_iterat
return {left, right};
}
+static bool omegaCompeare(std::pair<fvalue, std::vector<eis::DataPoint>::const_iterator> a,
+ std::pair<fvalue, std::vector<eis::DataPoint>::const_iterator> b)
+{
+ return a.first < b.first;
+}
+
+static std::vector<std::pair<fvalue, std::vector<eis::DataPoint>::const_iterator>>
+getSortedOmegaDistances(fvalue omega,
+ std::vector<eis::DataPoint>::const_iterator start,
+ std::vector<eis::DataPoint>::const_iterator end)
+{
+ std::vector<std::pair<fvalue, std::vector<eis::DataPoint>::const_iterator>> out;
+
+ for(std::vector<eis::DataPoint>::const_iterator it = start; it != end; it++)
+ out.push_back({std::abs(it->omega-omega), it});
+ std::sort(out.begin(), out.end(), omegaCompeare);
+ return out;
+}
+
static eis::DataPoint linearInterpolatePoint(fvalue omega, const eis::DataPoint& left, const eis::DataPoint& right)
{
assert(left.omega <= omega);
@@ -269,7 +290,42 @@ static eis::DataPoint linearInterpolatePoint(fvalue omega, const eis::DataPoint&
return eis::DataPoint(left.im+sloap*(omega - left.omega), omega);
}
-std::vector<eis::DataPoint> eis::fitToFrequencies(std::vector<fvalue> omegas, const std::vector<eis::DataPoint>& data)
+static eis::DataPoint linearExtrapoloatePoint(fvalue omega, const std::vector<eis::DataPoint>& data)
+{
+ if(data.size() < 3)
+ throw std::invalid_argument("extrapolation requires at least 3 points");
+
+ std::vector<std::pair<fvalue, std::vector<eis::DataPoint>::const_iterator>> dist;
+ dist = getSortedOmegaDistances(omega, data.begin(), data.end());
+
+ std::vector<fvalue> omegas;
+ std::vector<fvalue> im;
+ std::vector<fvalue> re;
+ omegas.reserve(dist.size());
+ im.reserve(dist.size());
+ re.reserve(dist.size());
+ for(size_t i = 0; i < 3; ++i)
+ {
+ omegas.push_back(dist[i].second->omega);
+ im.push_back(dist[i].second->im.imag());
+ re.push_back(dist[i].second->im.real());
+ }
+
+ LinearRegessionCalculator realReg(omegas, re);
+ LinearRegessionCalculator imagReg(omegas, im);
+
+ eis::Log(eis::Log::DEBUG)<<"Real regression:\n\toffset: "<<realReg.offset
+ <<"\n\tsloap: "<<realReg.slope<<"\n\tstderror: "<<realReg.stdError;
+ eis::Log(eis::Log::DEBUG)<<"Imag regression:\n\toffset: "<<imagReg.offset
+ <<"\n\tsloap: "<<imagReg.slope<<"\n\tstderror: "<<imagReg.stdError;
+
+ if(realReg.stdError > 1 || imagReg.stdError > 1)
+ throw std::invalid_argument("input data must be sufficantly linear");
+
+ return eis::DataPoint({realReg.slope*omega+realReg.offset, imagReg.slope*omega+imagReg.offset}, omega);
+}
+
+std::vector<eis::DataPoint> eis::fitToFrequencies(std::vector<fvalue> omegas, const std::vector<eis::DataPoint>& data, bool linearExtrapolation)
{
std::vector<eis::DataPoint> out;
out.reserve(omegas.size());
@@ -289,15 +345,45 @@ std::vector<eis::DataPoint> eis::fitToFrequencies(std::vector<fvalue> omegas, co
eis::Log(eis::Log::DEBUG)<<"\tRight "<<lr.second->omega<<','<<lr.second->im;
if(lr.first == lr.second)
+ {
dp.im = lr.first->im;
+ }
else if(lr.first != data.end() && lr.second != data.end())
+ {
dp = linearInterpolatePoint(dp.omega, *lr.first, *lr.second);
+ }
else if(lr.first != data.end() && lr.second == data.end())
- dp.im = lr.first->im;
+ {
+ try
+ {
+ if(!linearExtrapolation)
+ dp.im = lr.first->im;
+ else
+ dp = linearExtrapoloatePoint(dp.omega, data);
+ }
+ catch (const std::invalid_argument& ex)
+ {
+ dp.im = lr.first->im;
+ }
+ }
else if(lr.first == data.end() && lr.second != data.end())
- dp.im = lr.second->im;
+ {
+ try
+ {
+ if(!linearExtrapolation)
+ dp.im = lr.second->im;
+ else
+ dp = linearExtrapoloatePoint(dp.omega, data);
+ }
+ catch (const std::invalid_argument& ex)
+ {
+ dp.im = lr.second->im;
+ }
+ }
else
+ {
assert(false);
+ }
}
return out;
diff --git a/eisgenerator/basicmath.h b/eisgenerator/basicmath.h
index f1c73813c03ebb40b65043e9a673f9471c66daf6..1daab85c16628f8648497f53caa38730a36446b8 100644
--- a/eisgenerator/basicmath.h
+++ b/eisgenerator/basicmath.h
@@ -18,6 +18,8 @@ namespace eis
void noise(std::vector<eis::DataPoint>& data, double amplitude, bool relative);
void removeDuplicates(std::vector<eis::DataPoint>& data);
bool fvalueEq(fvalue a, fvalue b, unsigned int ulp = 4);
- std::vector<eis::DataPoint> fitToFrequencies(std::vector<fvalue> omegas, const std::vector<eis::DataPoint>& data);
+ std::vector<eis::DataPoint> fitToFrequencies(std::vector<fvalue> omegas,
+ const std::vector<eis::DataPoint>& data,
+ bool linearExtrapolation = false);
}
diff --git a/linearregession.cpp b/linearregession.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9d555670d1e1c822ac281a1afbb63ac92cea2c69
--- /dev/null
+++ b/linearregession.cpp
@@ -0,0 +1,38 @@
+#include "linearregession.h"
+#include <stdexcept>
+#include <cmath>
+
+LinearRegessionCalculator::LinearRegessionCalculator(const std::vector < fvalue > &xValues, const std::vector < fvalue > &yValues)
+{
+ if(xValues.size() == yValues.size())
+ {
+ this->xValues = xValues;
+ this->yValues = yValues;
+
+ fvalue sumY = 0;
+ fvalue sumX = 0;
+ fvalue sumXTimesY = 0;
+ fvalue sumSquaredX = 0;
+ fvalue sumSquaredY = 0;
+ for(unsigned i = 0; i < xValues.size(); i++)
+ {
+ sumY += yValues[i];
+ sumX += xValues[i];
+ sumSquaredX += xValues[i] * xValues[i];
+ sumSquaredY += yValues[i] * yValues[i];
+ sumXTimesY += xValues[i] * yValues[i];
+ }
+
+ slope = (sumXTimesY - (sumX * sumY) / xValues.size()) / (sumSquaredX - (sumX * sumX) / xValues.size());
+ offset = sumY / xValues.size() - slope * (sumX / xValues.size());
+
+ fvalue error =
+ (xValues.size() * sumSquaredY - sumY * sumY -
+ slope * slope * (xValues.size() * sumSquaredX - sumX * sumX)) / (xValues.size() * (xValues.size() - 2));
+ stdError = sqrt((error * error * xValues.size()) / (xValues.size() * sumSquaredX - sumX * sumX));
+ }
+ else
+ {
+ throw std::invalid_argument("xValues and yValues need to be the same size");
+ }
+}
diff --git a/linearregession.h b/linearregession.h
new file mode 100644
index 0000000000000000000000000000000000000000..5061e02c5413b67742238ebf4db9d24bfd3b1abe
--- /dev/null
+++ b/linearregession.h
@@ -0,0 +1,17 @@
+#pragma once
+
+#include <vector>
+#include "eisgenerator/eistype.h"
+
+class LinearRegessionCalculator
+{
+public:
+ std::vector < fvalue > xValues;
+ std::vector < fvalue > yValues;
+
+ fvalue slope;
+ fvalue offset;
+ fvalue stdError;
+
+ LinearRegessionCalculator(const std::vector < fvalue > &xValues, const std::vector < fvalue > &yValues);
+};