fvalueEq: handle eadge case at 0 where epsilon is zero

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)

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.

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.
 	*

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" : "");

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