From a18672760352dcca81fb1c4e74f79d8c02cba2cc Mon Sep 17 00:00:00 2001
From: Carl Philipp Klemm <philipp@uvos.xyz>
Date: Mon, 11 Mar 2024 16:11:40 +0100
Subject: [PATCH] fvalueEq: handle eadge case at 0 where epsilon is zero
---
basicmath.cpp | 4 ++--
eisgenerator/eistype.h | 3 ++-
eisgenerator/model.h | 22 +++++++++++++++++++++-
eistype.cpp | 14 ++++++++++++++
model.cpp | 34 ++++++++++++++++++++++++++++++++--
5 files changed, 71 insertions(+), 6 deletions(-)
diff --git a/basicmath.cpp b/basicmath.cpp
index 74cfcc9..f08bff2 100644
--- a/basicmath.cpp
+++ b/basicmath.cpp
@@ -248,7 +248,7 @@ void eis::removeDuplicates(std::vector<eis::DataPoint>& data)
bool eis::fvalueEq(fvalue a, fvalue b, unsigned int ulp)
{
fvalue epsilon = std::numeric_limits<fvalue>::epsilon()*std::fabs(a+b)*ulp;
- return a - epsilon < b && a + epsilon > b;
+ return a - epsilon <= b && a + epsilon >= b;
}
static std::pair<std::vector<eis::DataPoint>::const_iterator, std::vector<eis::DataPoint>::const_iterator>
@@ -412,7 +412,7 @@ std::vector<eis::DataPoint> eis::fitToFrequencies(std::vector<fvalue> omegas, co
return out;
}
-void difference(std::vector<eis::DataPoint>& a, const std::vector<eis::DataPoint>& b)
+void eis::difference(std::vector<eis::DataPoint>& a, const std::vector<eis::DataPoint>& b)
{
assert(a.size() == b.size());
for(size_t i = 0; i < a.size(); ++i)
diff --git a/eisgenerator/eistype.h b/eisgenerator/eistype.h
index ebb0557..d151162 100644
--- a/eisgenerator/eistype.h
+++ b/eisgenerator/eistype.h
@@ -445,9 +445,10 @@ public:
/**
* @brief Returns the a vector of DataPoints as a pair of valarrays.
*
+* @param data the data to transform to valarrays
* @return A pair of valarrays, first the real part and second the imaginary part.
*/
-std::pair<std::valarray<fvalue>, std::valarray<fvalue>> eisToValarrays(const std::vector<eis::DataPoint>& b);
+std::pair<std::valarray<fvalue>, std::valarray<fvalue>> eisToValarrays(const std::vector<eis::DataPoint>& data);
/**
* @brief Returns the mean l2 element wise distance of the given spectra.
diff --git a/eisgenerator/model.h b/eisgenerator/model.h
index 2f053da..e5f9dce 100644
--- a/eisgenerator/model.h
+++ b/eisgenerator/model.h
@@ -52,7 +52,7 @@ private:
static size_t paramSkipIndex(const std::string& str, size_t index);
static void addComponantToFlat(Componant* componant, std::vector<Componant*>* flatComponants);
- static void sweepThreadFn(std::vector<std::vector<DataPoint>>* data, Model* model, size_t start, size_t stop, const Range& omega);
+ static void sweepThreadFn(std::vector<std::vector<DataPoint>>* data, Model* model, size_t start, size_t stop, const std::vector<fvalue>& omega);
size_t getActiveParameterCount();
@@ -118,6 +118,26 @@ public:
*/
std::vector<DataPoint> executeSweep(const std::vector<fvalue>& omega, size_t index = 0);
+ /**
+ * @brief Executes a frequency and parameter sweep at the given parameter indecies
+ *
+ * @param omega The range along which to execute the frequency sweep.
+ * @param indecies the parameter indecies to include in the sweep
+ * @param parallel if this is set to true, the parameter sweep is executed in parallel
+ * @return A vector of vectors of DataPoint structs containing the impedance at every frequency sweep and parameter index.
+ */
+ std::vector<std::vector<DataPoint>> executeSweeps(const Range& omega, const std::vector<size_t>& indecies, bool parallel = false);
+
+ /**
+ * @brief Executes a frequency and parameter sweep at the given parameter indecies
+ *
+ * @param omega A vector of frequencies in rad/s to calculate the impedance at.
+ * @param indecies the parameter indecies to include in the sweep
+ * @param parallel if this is set to true, the parameter sweep is executed in parallel
+ * @return A vector of vectors of DataPoint structs containing the impedance at every frequency sweep and parameter index.
+ */
+ std::vector<std::vector<DataPoint>> executeSweeps(const std::vector<fvalue>& omega, const std::vector<size_t>& indecies, bool parallel = false);
+
/**
* @brief Executes a frequency sweep with the given omega values for each parameter combination in the applied parameter sweep.
*
diff --git a/eistype.cpp b/eistype.cpp
index 5900750..5d499f8 100644
--- a/eistype.cpp
+++ b/eistype.cpp
@@ -44,6 +44,20 @@ EisSpectra eis::loadFromDisk(const std::filesystem::path& path)
return EisSpectra(path);
}
+std::pair<std::valarray<fvalue>, std::valarray<fvalue>> eis::eisToValarrays(const std::vector<eis::DataPoint>& data)
+{
+ std::valarray<fvalue> real(data.size());
+ std::valarray<fvalue> imag(data.size());
+
+ for(size_t i = 0; i < data.size(); ++i)
+ {
+ real[i] = data[i].im.real();
+ imag[i] = data[i].im.imag();
+ }
+
+ return {real, imag};
+}
+
void eis::Range::print(int level) const
{
Log(static_cast<Log::Level>(level))<<"Range "<<start<<'-'<<end<<' '<<count<<" steps"<<(log ? " Log" : "");
diff --git a/model.cpp b/model.cpp
index 62eab45..c172bff 100644
--- a/model.cpp
+++ b/model.cpp
@@ -372,7 +372,37 @@ std::vector<DataPoint> Model::executeSweep(const std::vector<fvalue>& omega, siz
return results;
}
-void Model::sweepThreadFn(std::vector<std::vector<DataPoint>>* data, Model* model, size_t start, size_t stop, const Range& omega)
+std::vector<std::vector<DataPoint>> Model::executeSweeps(const Range& omega, const std::vector<size_t>& indecies, bool parallel)
+{
+ return executeSweeps(omega.getRangeVector(), indecies, parallel);
+}
+
+std::vector<std::vector<DataPoint>> Model::executeSweeps(const std::vector<fvalue>& omega, const std::vector<size_t>& indecies, bool parallel)
+{
+ unsigned int threadsCount = parallel ? std::thread::hardware_concurrency() : 1;
+
+ if(indecies.size() < threadsCount*10)
+ threadsCount = 1;
+
+ size_t countPerThread = indecies.size()/threadsCount;
+ std::vector<std::thread> threads(threadsCount);
+ std::vector<Model> models(threadsCount, *this);
+
+ std::vector<std::vector<DataPoint>> data(indecies.size());
+
+ for(size_t i = 0; i < threadsCount; ++i)
+ {
+ size_t start = i*countPerThread;
+ size_t stop = i < threadsCount-1 ? (i+1)*countPerThread : indecies.size();
+ threads[i] = std::thread(sweepThreadFn, &data, &models[i], start, stop, omega);
+ }
+ for(size_t i = 0; i < threadsCount; ++i)
+ threads[i].join();
+
+ return data;
+}
+
+void Model::sweepThreadFn(std::vector<std::vector<DataPoint>>* data, Model* model, size_t start, size_t stop, const std::vector<fvalue>& omega)
{
for(size_t i = start; i < stop; ++i)
{
@@ -396,7 +426,7 @@ std::vector<std::vector<DataPoint>> Model::executeAllSweeps(const Range& omega)
{
size_t start = i*countPerThread;
size_t stop = i < threadsCount-1 ? (i+1)*countPerThread : count;
- threads[i] = std::thread(sweepThreadFn, &data, &models[i], start, stop, std::ref(omega));
+ threads[i] = std::thread(sweepThreadFn, &data, &models[i], start, stop, omega.getRangeVector());
}
for(size_t i = 0; i < threadsCount; ++i)
threads[i].join();
--
GitLab