From 0e818d315d0f72152fb55dc690603b75dda47e7d Mon Sep 17 00:00:00 2001
From: Carl Philipp Klemm <philipp@uvos.xyz>
Date: Fri, 17 Nov 2023 14:30:09 +0100
Subject: [PATCH] basicmath: add support for extrapolateing spectra
---
CMakeLists.txt | 3 --
basicmath.cpp | 76 ++++++++++++++++++++++++++++++++++++++++
eisgenerator/basicmath.h | 2 +-
eisgenerator/eistype.h | 1 +
eistype.cpp | 55 ++++++++++++++++++-----------
main.cpp | 7 +++-
options.h | 32 ++++++++---------
7 files changed, 135 insertions(+), 41 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2d3c8de..9fec7ab 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -112,9 +112,6 @@ install(TARGETS ${PROJECT_NAME} DESTINATION lib)
install(TARGETS ${PROJECT_NAME}_static DESTINATION lib)
install(TARGETS ${PROJECT_NAME}_obj FILE_SET header_set DESTINATION include/${PROJECT_NAME})
install(FILES ${CMAKE_CURRENT_BINARY_DIR}/pkgconfig/libeisgenerator.pc DESTINATION lib/pkgconfig)
-install(FILES eisgenerator_plot DESTINATION bin PERMISSIONS WORLD_EXECUTE WORLD_READ
-OWNER_READ OWNER_WRITE OWNER_EXECUTE
-GROUP_READ GROUP_EXECUTE)
install(FILES kissplotcsv DESTINATION bin PERMISSIONS WORLD_EXECUTE WORLD_READ
OWNER_READ OWNER_WRITE OWNER_EXECUTE
GROUP_READ GROUP_EXECUTE )
diff --git a/basicmath.cpp b/basicmath.cpp
index d451c8b..84db4b3 100644
--- a/basicmath.cpp
+++ b/basicmath.cpp
@@ -4,6 +4,7 @@
#include <limits>
#include "eistype.h"
+#include "log.h"
static size_t gradIndex(size_t dataSize, size_t inputIndex)
{
@@ -226,3 +227,78 @@ 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;
}
+
+static std::pair<std::vector<eis::DataPoint>::const_iterator, std::vector<eis::DataPoint>::const_iterator>
+getLrClosest(const eis::DataPoint& dp, std::vector<eis::DataPoint>::const_iterator start, std::vector<eis::DataPoint>::const_iterator end)
+{
+
+ std::vector<eis::DataPoint>::const_iterator left = end;
+ fvalue distLeft = std::numeric_limits<fvalue>::max();
+ std::vector<eis::DataPoint>::const_iterator right = end;
+ fvalue distRight = std::numeric_limits<fvalue>::max();
+
+ for(std::vector<eis::DataPoint>::const_iterator it = start; it != end; it++)
+ {
+ if(eis::fvalueEq(it->omega, dp.omega))
+ return {it, it};
+ fvalue dist = it->omega-dp.omega;
+ bool sign = std::signbit(dist);
+ dist = std::abs(dist);
+
+ if(sign && (left == end || dist < distLeft))
+ {
+ distLeft = dist;
+ left = it;
+ }
+ else if(!sign && (right == end || dist < distRight))
+ {
+ distRight = dist;
+ right = it;
+ }
+ }
+ return {left, right};
+}
+
+static eis::DataPoint linearInterpolatePoint(fvalue omega, const eis::DataPoint& left, const eis::DataPoint& right)
+{
+ assert(left.omega <= omega);
+ assert(right.omega >= omega);
+
+ fvalue omegaDif = right.omega - left.omega;
+ std::complex<fvalue> sloap = (right.im - left.im)/omegaDif;
+ return eis::DataPoint(left.im+sloap*(omega - left.omega), omega);
+}
+
+std::vector<eis::DataPoint> eis::fitToFrequencies(std::vector<fvalue> omegas, const std::vector<eis::DataPoint>& data)
+{
+ std::vector<eis::DataPoint> out;
+ out.reserve(omegas.size());
+ for(fvalue omega : omegas)
+ out.push_back(eis::DataPoint({0,0}, omega));
+
+ eis::Log(eis::Log::DEBUG)<<__func__<<':';
+
+ for(eis::DataPoint& dp : out)
+ {
+ auto lr = getLrClosest(dp, data.begin(), data.end());
+
+ eis::Log(eis::Log::DEBUG)<<"\tValue for "<<dp.omega;
+ if(lr.first != data.end())
+ eis::Log(eis::Log::DEBUG)<<"\tLeft "<<lr.first->omega<<','<<lr.first->im;
+ if(lr.second != data.end())
+ eis::Log(eis::Log::DEBUG)<<"\tRight "<<lr.second->omega<<','<<lr.second->im;
+
+ if(lr.first == lr.second)
+ dp.im = lr.first->im;
+ else if(lr.first != data.end() && lr.second != data.end())
+ dp = linearInterpolatePoint(dp.omega, *lr.first, *lr.second);
+ else if(lr.first != data.end() && lr.second == data.end())
+ dp.im = lr.first->im;
+ else if(lr.first == data.end() && lr.second != data.end())
+ dp.im = lr.second->im;
+ else
+ assert(false);
+ }
+
+ return out;
+}
diff --git a/eisgenerator/basicmath.h b/eisgenerator/basicmath.h
index a14b4ab..f1c7381 100644
--- a/eisgenerator/basicmath.h
+++ b/eisgenerator/basicmath.h
@@ -1,6 +1,5 @@
#pragma once
#include <vector>
-#include "model.h"
#include "eistype.h"
namespace eis
@@ -19,5 +18,6 @@ namespace eis
void noise(std::vector<eis::DataPoint>& data, double amplitude, bool relative);
void removeDuplicates(std::vector<eis::DataPoint>& data);
bool fvalueEq(fvalue a, fvalue b, unsigned int ulp = 4);
+ std::vector<eis::DataPoint> fitToFrequencies(std::vector<fvalue> omegas, const std::vector<eis::DataPoint>& data);
}
diff --git a/eisgenerator/eistype.h b/eisgenerator/eistype.h
index 09aa1ae..78ed81b 100644
--- a/eisgenerator/eistype.h
+++ b/eisgenerator/eistype.h
@@ -122,6 +122,7 @@ public:
bool isSane() const;
std::vector<fvalue> getRangeVector() const;
+ static Range fromString(std::string str, size_t count);
static std::vector<Range> rangesFromParamString(const std::string& paramStr, size_t count);
};
diff --git a/eistype.cpp b/eistype.cpp
index e1c97dd..a17f4c3 100644
--- a/eistype.cpp
+++ b/eistype.cpp
@@ -36,6 +36,40 @@ std::vector<fvalue> eis::Range::getRangeVector() const
return out;
}
+eis::Range eis::Range::fromString(std::string str, size_t count)
+{
+ bool log = false;
+ std::vector<std::string> tokens = tokenize(str, '~');
+ eis::Range out;
+
+ if(str.back() == 'l' || str.back() == 'L')
+ {
+ log = true;
+ str.pop_back();
+ }
+
+ try
+ {
+ if(tokens.size() == 1)
+ {
+ out = Range(std::stod(tokens[0]), std::stod(tokens[0]), 1, log);
+ }
+ else
+ {
+ out = Range(std::stod(tokens[0]), std::stod(tokens[1]), count, log);
+ if(tokens.size() > 2)
+ throw std::invalid_argument("");
+ }
+
+ }
+ catch(const std::invalid_argument& ia)
+ {
+ throw std::invalid_argument("invalid parameter \""+ str + "\"");
+ }
+
+ return out;
+}
+
std::vector<Range> eis::Range::rangesFromParamString(const std::string& paramStr, size_t count)
{
std::vector<std::string> tokens = tokenize(paramStr, ',');
@@ -43,29 +77,10 @@ std::vector<Range> eis::Range::rangesFromParamString(const std::string& paramStr
std::vector<Range> ranges(tokens.size());
for(size_t i = 0; i < tokens.size(); ++i)
{
- bool log = false;
std::string& token = tokens[i];
- std::vector<std::string> subTokens = tokenize(token, '~');
-
- if(token.back() == 'l' || token.back() == 'L')
- {
- log = true;
- token.pop_back();
- }
-
try
{
- if(subTokens.size() == 1)
- {
- ranges[i] = Range(std::stod(subTokens[0]), std::stod(subTokens[0]), 1, log);
- }
- else
- {
- ranges[i] = Range(std::stod(subTokens[0]), std::stod(subTokens[1]), count, log);
- if(subTokens.size() > 2)
- throw std::invalid_argument("");
- }
-
+ ranges[i] = fromString(token, count);
}
catch(const std::invalid_argument& ia)
{
diff --git a/main.cpp b/main.cpp
index 8667dca..6c058da 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,4 +1,3 @@
-#include <eisgenerator/eistype.h>
#include <iostream>
#include <complex>
#include <chrono>
@@ -60,6 +59,12 @@ static void runSweep(const Config& config, eis::Model& model)
eis::Log(eis::Log::INFO)<<"normalized results:";
eis::normalize(results);
}
+ else if(config.extrapolate)
+ {
+ std::vector<fvalue> exrapolateOmegas = config.extrapolateRange.getRangeVector();
+ results = eis::fitToFrequencies(exrapolateOmegas, results);
+ eis::Log(eis::Log::INFO)<<"extrapolated results:";
+ }
else
{
eis::Log(eis::Log::INFO)<<"results:";
diff --git a/options.h b/options.h
index f866547..4db314e 100644
--- a/options.h
+++ b/options.h
@@ -17,7 +17,8 @@ static struct argp_option options[] =
{"quiet", 'q', 0, 0, "only output data" },
{"param", 's', "[COUNT]", 0, "parameter sweep steps" },
{"model", 'm', "[STRING]", 0, "set model string" },
- {"omega", 'o', "[START-END]", 0, "set omega range" },
+ {"omega", 'o', "[START~END]", 0, "set omega range" },
+ {"extrapolate", 'a', "[START~END]", 0, "extrapolate a spectra simulated on the omega range to the one given here" },
{"omegasteps", 'c', "[COUNT]", 0, "set omega range steps" },
{"linear", 'l', 0, 0, "use linear instead of logarithmic steps" },
{"normalize", 'n', 0, 0, "normalize values" },
@@ -62,6 +63,8 @@ struct Config
int mode = MODE_NORMAL;
size_t paramSteps = 10;
eis::Range omegaRange;
+ bool extrapolate = false;
+ eis::Range extrapolateRange;
bool normalize = false;
bool reduce = false;
bool hertz = false;
@@ -152,26 +155,23 @@ parse_opt (int key, char *arg, struct argp_state *state)
break;
case 'o':
{
- std::vector<std::string> tokens = tokenize(std::string(arg), '-');
- if(tokens.size() != 2)
+ try
+ {
+ config->omegaRange = eis::Range::fromString(std::string(arg), config->omegaRange.count);
+ }
+ catch(const std::invalid_argument& ia)
{
argp_usage(state);
- break;
}
+ break;
+ }
+ case 'a':
+ {
try
{
- double start = std::stod(tokens[0]);
- double end = std::stod(tokens[1]);
- if(start < end)
- {
- config->omegaRange.start = start;
- config->omegaRange.end = end;
- }
- else
- {
- config->omegaRange.start = end;
- config->omegaRange.end = start;
- }
+ config->extrapolateRange = eis::Range::fromString(std::string(arg), config->omegaRange.count);
+ config->extrapolateRange.log = config->omegaRange.log;
+ config->extrapolate = true;
}
catch(const std::invalid_argument& ia)
{
--
GitLab