From 77305f09753a0414ae0090ed1ad966b7414b92bb Mon Sep 17 00:00:00 2001
From: Carl Philipp Klemm <philipp@uvos.xyz>
Date: Tue, 17 Jan 2023 14:22:24 +0100
Subject: [PATCH] eistype: add operators to do math operations on DataPoints as
complex numbers, add eisNyquistDistance
---
eisgenerator/eistype.h | 42 ++++++++++++++++++++++++++++++++--
eistype.cpp | 52 ++++++++++++++++++++++++++++++++++--------
2 files changed, 83 insertions(+), 11 deletions(-)
diff --git a/eisgenerator/eistype.h b/eisgenerator/eistype.h
index 6224047..c1c0466 100644
--- a/eisgenerator/eistype.h
+++ b/eisgenerator/eistype.h
@@ -2,6 +2,7 @@
#include <complex>
#include <vector>
#include <cassert>
+#include <cmath>
typedef double fvalue;
@@ -12,10 +13,45 @@ class DataPoint
public:
std::complex<fvalue> im;
fvalue omega;
- bool operator==(const DataPoint& in)
+ DataPoint() = default;
+ DataPoint(std::complex<fvalue> imIn, fvalue omegaIn = 0): im(imIn), omega(omegaIn){}
+ bool operator==(const DataPoint& in) const
{
return im == in.im;
}
+ DataPoint operator-(const DataPoint& in) const
+ {
+ DataPoint out(*this);
+ out.im = out.im - in.im;
+ return out;
+ }
+ DataPoint operator+(const DataPoint& in) const
+ {
+ DataPoint out(*this);
+ out.im = out.im + in.im;
+ return out;
+ }
+ DataPoint operator/(fvalue in) const
+ {
+ DataPoint out(*this);
+ out.im = out.im / in;
+ return out;
+ }
+ DataPoint operator*(fvalue in) const
+ {
+ DataPoint out(*this);
+ out.im = out.im * in;
+ return out;
+ }
+ DataPoint operator/=(fvalue in)
+ {
+ im = im / in;
+ return *this;
+ }
+ fvalue complexVectorLength() const
+ {
+ return std::sqrt(std::pow(im.real(), 2) + std::pow(im.imag(), 2));
+ }
};
class Range
@@ -75,6 +111,8 @@ public:
bool saveToDisk(const std::vector<DataPoint>& data, const std::string& fileName, std::string headStr = "");
-double eisDistance(const std::vector<eis::DataPoint>& a, const std::vector<eis::DataPoint>& b);
+fvalue eisDistance(const std::vector<eis::DataPoint>& a, const std::vector<eis::DataPoint>& b);
+
+fvalue eisNyquistDistance(const std::vector<eis::DataPoint>& a, const std::vector<eis::DataPoint>& b);
}
diff --git a/eistype.cpp b/eistype.cpp
index 6276fdd..5c43785 100644
--- a/eistype.cpp
+++ b/eistype.cpp
@@ -1,6 +1,7 @@
#include "eistype.h"
#include <fstream>
#include <sstream>
+#include <algorithm>
#include "strops.h"
#include "log.h"
@@ -22,9 +23,7 @@ bool eis::saveToDisk(const std::vector<DataPoint>& data, const std::string& file
file<<"omega,real,im\n";
for(const eis::DataPoint& point : data)
- {
file<<point.omega<<','<<point.im.real()<<','<<point.im.imag()<<'\n';
- }
file.close();
return true;
}
@@ -61,13 +60,13 @@ std::vector<Range> eis::Range::rangesFromParamString(const std::string& paramStr
{
ranges[i] = Range(std::stod(subTokens[0]), std::stod(subTokens[1]), count, log);
if(subTokens.size() > 2)
- Log(Log::WARN)<<"invalid parameter string "<<paramStr<<" more that two arguments at range "<<i;
+ throw std::invalid_argument("");
}
}
catch(const std::invalid_argument& ia)
{
- Log(Log::WARN)<<"invalid parameter string "<<paramStr;
+ throw std::invalid_argument("invalid parameter string \"{"+ paramStr + "}\"");
}
}
return ranges;
@@ -79,9 +78,11 @@ std::string eis::Range::getString() const
ss<<start;
if(count > 1)
+ {
ss<<'~'<<end;
- if(log)
- ss<<'L';
+ if(log)
+ ss<<'L';
+ }
return ss.str();
}
@@ -95,7 +96,8 @@ bool eis::Range::isSane() const
return true;
}
-double eis::eisDistance(const std::vector<eis::DataPoint>& a, const std::vector<eis::DataPoint>& b)
+//Compute simmuliarity on a bode plot
+fvalue eis::eisDistance(const std::vector<eis::DataPoint>& a, const std::vector<eis::DataPoint>& b)
{
assert(a.size() == b.size());
@@ -104,7 +106,39 @@ double eis::eisDistance(const std::vector<eis::DataPoint>& a, const std::vector<
{
double diffRe = std::pow(b[i].im.real() - a[i].im.real(), 2);
double diffIm = std::pow(b[i].im.imag() - a[i].im.imag(), 2);
- accum += std::sqrt(diffRe+diffIm);
+ accum += diffRe+diffIm;
+ }
+ return sqrt(accum/a.size());
+}
+
+//Compute simmuliarity on a nyquist plot
+fvalue eis::eisNyquistDistance(const std::vector<eis::DataPoint>& a, const std::vector<eis::DataPoint>& b)
+{
+ assert(a.size() > 2 && b.size() > 3);
+ double accum = 0;
+ for(size_t i = 0; i < a.size(); ++i)
+ {
+ std::vector<std::pair<double, const eis::DataPoint*>> distances;
+ for(size_t j = 0; j < b.size(); ++j)
+ {
+ double diffRe = std::pow(b[j].im.real() - a[i].im.real(), 2);
+ double diffIm = std::pow(b[j].im.imag() - a[i].im.imag(), 2);
+ std::pair<double, const eis::DataPoint*> dp;
+ dp.first = sqrt(diffRe+diffIm);
+ dp.second = &b[j];
+ distances.push_back(dp);
+ }
+ std::sort(distances.begin(), distances.end(),
+ [](const std::pair<double, const eis::DataPoint*>& a, const std::pair<double, const eis::DataPoint*>& b) -> bool
+ {return a.first < b.first;});
+
+ eis::DataPoint base = (*distances[0].second)-(*distances[1].second);
+ base = base/base.complexVectorLength();
+ eis::DataPoint diff = (*distances[0].second)-a[i];
+ diff = diff/diff.complexVectorLength();
+ fvalue dprod = base.im.real()*diff.im.real() + base.im.imag()*diff.im.imag();
+ fvalue dist = diff.complexVectorLength()*(1-dprod);
+ accum += std::pow(dist, 2);
}
- return accum/a.size();
+ return std::sqrt(accum/a.size());
}
--
GitLab