From fd064cd63af38b771a92cd5f7a8156c65d3a4c4b Mon Sep 17 00:00:00 2001
From: Carl Philipp Klemm <philipp@uvos.xyz>
Date: Tue, 30 Apr 2024 16:40:22 +0200
Subject: [PATCH] allow limiting of noise intenisty

---
 eisnoise.cpp        | 36 +++++++++++++++++++++++---
 eisnoise/eisnoise.h |  5 ++--
 main.cpp            | 61 ++++++++++++++++++++++++++++++++-------------
 3 files changed, 80 insertions(+), 22 deletions(-)

diff --git a/eisnoise.cpp b/eisnoise.cpp
index a2e28d9..f9cfd29 100644
--- a/eisnoise.cpp
+++ b/eisnoise.cpp
@@ -1,13 +1,16 @@
 #include "eisnoise/eisnoise.h"
 
+#include <limits>
 #include <numeric>
 #include <sstream>
 #include <random>
+#include <iostream>
 #include <list>
 #define INCBIN_PREFIX r
 #include "incbin.h"
 #include "microtar.h"
-#include "eistype.h"
+#include <eisgenerator/eistype.h>
+#include <eisgenerator/basicmath.h>
 
 INCBIN(Tar, DATA_DIR "/data.tar");
 
@@ -40,7 +43,7 @@ void EisNoise::addImpl(std::vector<eis::DataPoint>& spectra, const Noise& noise,
 		if(std::abs(spectra[i].im) > 0.001)
 			spectra[i].im = (noiseValue * spectra[i].im) * intensity + spectra[i].im;
 		else
-			(noiseValue * static_cast<fvalue>(0.001)) * intensity + spectra[i].im;
+			spectra[i].im = (noiseValue * static_cast<fvalue>(0.001)) * intensity + spectra[i].im;
 	}
 }
 
@@ -96,7 +99,25 @@ std::vector<std::string> EisNoise::getAvailable()
 	return out;
 }
 
-EisNoise::EisNoise()
+fvalue EisNoise::calcMaxIntensity(std::vector<eis::DataPoint>& spectra)
+{
+	fvalue reMax = std::numeric_limits<fvalue>::min();
+	fvalue imMax = std::numeric_limits<fvalue>::min();
+
+	for(const eis::DataPoint& point : spectra)
+	{
+		fvalue re = std::abs(point.im.real());
+		fvalue im = std::abs(point.im.imag());
+		if(im > imMax)
+			imMax = im;
+		if(re > reMax)
+			reMax = re;
+	}
+
+	return std::max(reMax, imMax);
+}
+
+EisNoise::EisNoise(fvalue maxIntensity)
 {
 	mtar_t tar;
 
@@ -115,6 +136,15 @@ EisNoise::EisNoise()
 		noise.omegaMax = noise.spectra.data.back().omega;
 		noise.omegaMin = noise.spectra.data.front().omega;
 		assert(noise.omegaMax > noise.omegaMin);
+
+
+		if(maxIntensity > 0)
+		{
+			fvalue maxIntensityForSpectra = calcMaxIntensity(noise.spectra.data);
+			if(maxIntensityForSpectra > maxIntensity)
+				eis::mulitplyAdd(noise.spectra.data, maxIntensity/maxIntensityForSpectra, 0);
+		}
+
 		noises.push_back(noise);
 		mtar_next(&tar);
 	}
diff --git a/eisnoise/eisnoise.h b/eisnoise/eisnoise.h
index 8408db4..6310699 100644
--- a/eisnoise/eisnoise.h
+++ b/eisnoise/eisnoise.h
@@ -20,6 +20,8 @@ private:
 
 	static void addImpl(std::vector<eis::DataPoint>& spectra, const Noise& noise, fvalue intensity, bool ignFreq);
 
+	fvalue calcMaxIntensity(std::vector<eis::DataPoint>& spectra);
+
 public:
 	class apply_err: public std::exception
 	{
@@ -33,11 +35,10 @@ public:
 		}
 	};
 
-	EisNoise();
+	EisNoise(fvalue maxIntensity = -1);
 	~EisNoise() = default;
 
 	void add(std::vector<eis::DataPoint>& spectra, fvalue intensity = 1.0, bool ignFreq = false);
 	void addSpecific(std::vector<eis::DataPoint>& spectra, const std::string& name , fvalue intensity = 1.0,  bool ignFreq = false);
 	std::vector<std::string> getAvailable();
-
 };
diff --git a/main.cpp b/main.cpp
index ea795a0..dd5e199 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,10 +1,13 @@
 #include "eisnoise/eisnoise.h"
 #include "eistype.h"
+#include "log.h"
 #include <eisgenerator/model.h>
 #include <iostream>
 #include <valarray>
 #include <vector>
 #include <sciplot/sciplot.hpp>
+#include <random>
+#include <eisgenerator/basicmath.h>
 
 
 int main(int argc, char** argv)
@@ -18,38 +21,62 @@ int main(int argc, char** argv)
 
 	eis::Model model(argv[1], 100, true);
 	eis::Range omega(1, 1e6, 50, true);
+	eis::Log::level = eis::Log::DEBUG;
 
 	std::vector<size_t> recommendedIndecies = model.getRecommendedParamIndices(omega, 0.01, true);
 
-	std::vector<eis::DataPoint> data = model.executeSweep(omega);
-	EisNoise noise;
+	EisNoise noise(0.2);
+	std::random_device rd;
+	std::mt19937_64 gen(rd());
+	std::uniform_int_distribution<size_t> distrib(0, recommendedIndecies.size()-1);
 
-	std::vector<std::string> names = noise.getAvailable();
+	/*std::vector<std::string> names = noise.getAvailable();
 	for(const std::string& name : names)
-		std::cout<<name<<'\n';
+		std::cout<<name<<'\n';*/
 
 	sciplot::Plot2D plot;
 	plot.xlabel("Re");
 	plot.ylabel("Im");
 
+	double xmin = std::numeric_limits<double>::max();
+	double xmax = std::numeric_limits<double>::min();
+	double ymin = std::numeric_limits<double>::max();
+	double ymax = std::numeric_limits<double>::min();
+
 	for(size_t i = 0; i < 100; ++i)
 	{
-		std::vector<eis::DataPoint> nosiyData = data;
-		noise.add(nosiyData);
-
-		std::pair<std::valarray<fvalue>, std::valarray<fvalue>> arrays = eis::eisToValarrays(nosiyData);
-
-		double xMin = arrays.first.min();
-		double xMax = arrays.first.max();
-		double yMin = arrays.second.min();
-		double yMax = arrays.second.max();
-
-		plot.xrange(xMin, xMax);
-		plot.yrange(yMin, yMax);
+		size_t index = distrib(gen);
+		std::vector<eis::DataPoint> data = model.executeSweep(omega, recommendedIndecies[index]);
+		noise.add(data);
+
+		std::pair<std::valarray<fvalue>, std::valarray<fvalue>> arrays = eis::eisToValarrays(data);
+
+		double cxmin = arrays.first.min();
+		cxmin = cxmin - std::abs(cxmin*0.05);
+		double cymin = arrays.second.min();
+		cymin = cymin - std::abs(cymin*0.05);
+
+		double cxmax = arrays.first.max();
+		cxmax = cxmax - std::abs(cxmax*0.05);
+		double cymax = arrays.second.max();
+		cymax = cymax - std::abs(cymax*0.05);
+
+		if(xmax < cxmax)
+			xmax = cxmax;
+		if(ymax < cymax)
+			ymax = cymax;
+
+		if(xmin > cxmin)
+			xmin = cxmin;
+		if(ymin > cymin)
+			ymin = cymin;
+		plot.xrange(xmin, xmax);
+		plot.yrange(ymin, ymax);
 
 		plot.drawCurve(arrays.first, arrays.second);
 
-		eis::EisSpectra(nosiyData, model.getModelStrWithParam(), "").saveToStream(std::cout);
+		eis::EisSpectra(data, model.getModelStrWithParam(index), "").saveToStream(std::cout);
+		std::cout<<eis::nonConstantScore(data)<<'\n';
 
 	}
 
-- 
GitLab