eistype: add operators to do math operations on DataPoints as complex numbers,...

eistype: add operators to do math operations on DataPoints as complex numbers, add eisNyquistDistance

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);
 
 }

eistype.cpp

 #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';
+	}
 
 	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());
 }