diff --git a/CMakeLists.txt b/CMakeLists.txt index 006c316b7ef7412d6cedc7ef3fa3532bf769c2f6..5932d7129bd2cf4df1370363e047edd0791c5460 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,6 +11,8 @@ set(SRC_FILES model.cpp tokenize.cpp paralellseriel.cpp + log.cpp + normalize.cpp ) set(API_HEADERS @@ -20,7 +22,7 @@ set(API_HEADERS eisgenerator/constantphase.h eisgenerator/warburg.h eisgenerator/model.h - eisgenerator/tokenize.h + eisgenerator/log.h eisgenerator/paralellseriel.h ) @@ -32,14 +34,28 @@ set_target_properties(${PROJECT_NAME} PROPERTIES COMPILE_FLAGS "-std=c++20 -Wall set(CMAKE_INSTALL_PREFIX "/usr") install(TARGETS ${PROJECT_NAME} DESTINATION lib) install(FILES ${API_HEADERS} DESTINATION include/${PROJECT_NAME}) +install(FILES eisgenerator_plot DESTINATION bin PERMISSIONS WORLD_EXECUTE WORLD_READ +OWNER_READ OWNER_WRITE OWNER_EXECUTE +GROUP_READ GROUP_EXECUTE ) link_directories(${CMAKE_CURRENT_BINARY_DIR}) -set(SRC_FILES_TEST_APP main.cpp) -set(LIBS_TEST -l${PROJECT_NAME}) +set(SRC_FILES_TEST_APP test.cpp) +set(LIBS_TEST -L. -l${PROJECT_NAME}) add_executable(${PROJECT_NAME}_test ${SRC_FILES_TEST_APP}) add_dependencies(${PROJECT_NAME}_test ${PROJECT_NAME}) target_link_libraries(${PROJECT_NAME}_test ${LIBS_TEST}) target_include_directories(${PROJECT_NAME}_test PUBLIC eisgenerator) set_target_properties(${PROJECT_NAME}_test PROPERTIES COMPILE_FLAGS "-std=c++20 -Wall -O3 -march=native -g" LINK_FLAGS "-flto") - install(TARGETS ${PROJECT_NAME}_test DESTINATION bin) + +link_directories(${CMAKE_CURRENT_BINARY_DIR}) +set(SRC_FILES_TEST_APP main.cpp) +set(LIBS_TEST -L. -l${PROJECT_NAME}) +add_executable(${PROJECT_NAME}_export ${SRC_FILES_TEST_APP}) +add_dependencies(${PROJECT_NAME}_export ${PROJECT_NAME}) +target_link_libraries(${PROJECT_NAME}_export ${LIBS_TEST}) +target_include_directories(${PROJECT_NAME}_export PUBLIC eisgenerator) +set_target_properties(${PROJECT_NAME}_export PROPERTIES COMPILE_FLAGS "-std=c++20 -Wall -O3 -march=native -g" LINK_FLAGS "-flto") +install(TARGETS ${PROJECT_NAME}_export DESTINATION bin) + + diff --git a/cap.cpp b/cap.cpp index 0bed2004a2a9a707be626f2ba0d5facb04af49cf..2eee9f3ec047ccec9f534cc48dcc7087d11c2361 100644 --- a/cap.cpp +++ b/cap.cpp @@ -3,6 +3,10 @@ #include <cstdlib> #include <math.h> +#include "log.h" + +using namespace eis; + Cap::Cap(double c): _C(c) { @@ -13,7 +17,7 @@ Cap::Cap(std::string paramStr) std::vector<std::string> tokens = tokenize(paramStr, ','); if(tokens.empty()) { - std::cout<<"Warning: to few parameters in "<<__func__<<" parameter string: "<<paramStr<<'\n'; + Log(Log::WARN)<<"to few parameters in "<<__func__<<" parameter string: "<<paramStr<<'\n'; _C = 1e-6; return; } @@ -25,7 +29,7 @@ Cap::Cap(std::string paramStr) } catch(const std::invalid_argument& ia) { - std::cout<<"Warning: cant parse parameter in "<<__func__<<" parameter: "<<tokens[0]<<'\n'; + Log(Log::WARN)<<"Warning: cant parse parameter in "<<__func__<<" parameter: "<<tokens[0]<<'\n'; _C = 1e3; } } @@ -45,7 +49,7 @@ void Cap::setParam(const std::vector<double>& param) { if(param.empty()) { - std::cout<<"Warning: invalid parameter list sent to "<<__func__<<'\n'; + Log(Log::WARN)<<"invalid parameter list sent to "<<__func__<<'\n'; return; } diff --git a/componant.cpp b/componant.cpp index 395f12ac4551d725758aa614ab29a2068925c398..c0542bad0d53cf79f3bd395e19720af5c2219dca 100644 --- a/componant.cpp +++ b/componant.cpp @@ -5,6 +5,9 @@ #include "cap.h" #include "constantphase.h" #include "warburg.h" +#include "log.h" + +using namespace eis; Componant* Componant::copy(Componant* componant) { @@ -23,7 +26,7 @@ Componant* Componant::copy(Componant* componant) case 's': return new Serial(*dynamic_cast<Serial*>(componant)); default: - std::cout<<"Error: unimplmented type copy for "<<getComponantChar(componant)<<'\n'; + Log(Log::ERROR)<<"unimplmented type copy for "<<getComponantChar(componant)<<'\n'; assert(0); break; } diff --git a/constantphase.cpp b/constantphase.cpp index fa86aa0b129f4e24fdf3c13bc916c63c87054b47..a3098eadaecde80c36435e2e84af5967b96b365f 100644 --- a/constantphase.cpp +++ b/constantphase.cpp @@ -3,6 +3,10 @@ #include <cstdlib> #include <math.h> +#include "log.h" + +using namespace eis; + Cpe::Cpe(double q, double alpha): _Q(q), _alpha(alpha) { diff --git a/eisgenerator/cap.h b/eisgenerator/cap.h index e8c72bd28e37dc3377d626046ecb85c68818710d..dc5001eece98ddbf480dea6fe6ee0c48db84f616 100644 --- a/eisgenerator/cap.h +++ b/eisgenerator/cap.h @@ -4,6 +4,9 @@ #include <vector> #include "componant.h" +namespace eis +{ + class Cap: public Componant { private: @@ -17,3 +20,5 @@ public: virtual size_t paramCount() override; virtual ~Cap() = default; }; + +} diff --git a/eisgenerator/componant.h b/eisgenerator/componant.h index ebf924ab96524ccfa5b114ab175c163dd9a84664..8d2ca10387b36c8ec85738c687e8a111f3974b3a 100644 --- a/eisgenerator/componant.h +++ b/eisgenerator/componant.h @@ -3,6 +3,9 @@ #include <iostream> #include <vector> +namespace eis +{ + class Componant { public: @@ -23,3 +26,5 @@ class Componant static Componant* copy(Componant* componant); static char getComponantChar(Componant* componant); }; + +} diff --git a/eisgenerator/constantphase.h b/eisgenerator/constantphase.h index 8092832e32ba719964e7863c23866effb990ac6c..c9e35a14d5b45b7071138838e5eb1d024909fcfd 100644 --- a/eisgenerator/constantphase.h +++ b/eisgenerator/constantphase.h @@ -4,6 +4,9 @@ #include <vector> #include "componant.h" +namespace eis +{ + class Cpe: public Componant { private: @@ -18,3 +21,5 @@ public: virtual size_t paramCount() override; virtual ~Cpe() = default; }; + +} diff --git a/eisgenerator/model.h b/eisgenerator/model.h index 7e53901e53f48581b8cbc32692312741705deeff..45f7a4569cbff11416e04ac534645db33de9d5f5 100644 --- a/eisgenerator/model.h +++ b/eisgenerator/model.h @@ -7,28 +7,32 @@ #include "componant.h" -class Model +namespace eis { -public: - struct DataPoint - { - std::complex<double> im; - double omega; - }; - struct Range + +struct DataPoint +{ + std::complex<double> im; + double omega; +}; + +struct Range +{ + double start; + double end; + size_t count; + bool log = false; + + double stepSize() const { - double start; - double end; - size_t count; - - double stepSize() const - { - return (end-start)/count; - } - Range(double startI, double endI, size_t countI): start(startI), end(endI), count(countI){} - Range() = default; - }; + return (end-start)/count; + } + Range(double startI, double endI, size_t countI, bool logI = false): start(startI), end(endI), count(countI), log(logI){} + Range() = default; +}; +class Model +{ private: size_t opposingBraket(const std::string& str, size_t index, char bracketChar = ')'); size_t deepestBraket(const std::string& str); @@ -55,9 +59,12 @@ public: std::string getModelStr(); std::vector<Componant*> getFlatComponants(); + std::vector<double> getFlatParameters(); size_t getFlatParametersCount(); bool setFlatParameters(const std::vector<double>& parameters); static std::vector<double> getSweepParamByIndex(const std::vector<Range>& componantRanges, size_t index); static size_t getRequiredStepsForSweeps(const std::vector<Range>& componantRanges); }; + +} diff --git a/eisgenerator/normalize.h b/eisgenerator/normalize.h new file mode 100644 index 0000000000000000000000000000000000000000..035b5844a5682db19276d3592f53f7d8598a370c --- /dev/null +++ b/eisgenerator/normalize.h @@ -0,0 +1,16 @@ +#pragma once +#include <vector> +#include <string> + +#include "eistype.h" +#include "model.h" + +namespace eis +{ + +void normalize(std::vector<eis::DataPoint>& data); +std::vector<eis::DataPoint> reduceRegion(const std::vector<eis::DataPoint>& data, fvalue gradThreshFactor = 0.01); +std::complex<fvalue> absGrad(const std::vector<eis::DataPoint>& data, size_t index); +void eraseSingularites(std::vector<eis::DataPoint>& data); + +} diff --git a/eisgenerator/paralellseriel.h b/eisgenerator/paralellseriel.h index 896370d2432a6bb3721fb00b1e05ac996947fc4e..2e027cb9df0f4e86115c3cd91c672fe0bf27bb47 100644 --- a/eisgenerator/paralellseriel.h +++ b/eisgenerator/paralellseriel.h @@ -3,6 +3,9 @@ #include <complex> #include "componant.h" +namespace eis +{ + class Parallel: public Componant { public: @@ -26,3 +29,5 @@ public: ~Serial(); virtual std::complex<double> execute(double omaga) override; }; + +} diff --git a/eisgenerator/resistor.h b/eisgenerator/resistor.h index d21a20e7326989787578362e819bd97436f99d25..7462b673492e428e401fdfa1b0a3a5a0406e6133 100644 --- a/eisgenerator/resistor.h +++ b/eisgenerator/resistor.h @@ -2,6 +2,9 @@ #include "componant.h" #include <string> +namespace eis +{ + class Resistor: public Componant { private: @@ -16,3 +19,5 @@ public: virtual size_t paramCount() override; virtual ~Resistor() = default; }; + +} diff --git a/eisgenerator/tokenize.h b/eisgenerator/tokenize.h deleted file mode 100644 index 5f36a46bfc73ef851a72abb638c739f8c27839cf..0000000000000000000000000000000000000000 --- a/eisgenerator/tokenize.h +++ /dev/null @@ -1,6 +0,0 @@ -#pragma once -#include <string> -#include <vector> -#include <sstream> - -std::vector<std::string> tokenize(const std::string& str, const char delim = ' ', const char ignBracketStart = '\0', const char ignBracketEnd = '\0'); diff --git a/eisgenerator/warburg.h b/eisgenerator/warburg.h index f9f4011fb3fc711fd79598c0c23f00faa6d56af3..635e49b266535ae8c12a2d5555747f38095f6e9e 100644 --- a/eisgenerator/warburg.h +++ b/eisgenerator/warburg.h @@ -4,6 +4,9 @@ #include <vector> #include "componant.h" +namespace eis +{ + class Warburg: public Componant { private: @@ -17,3 +20,5 @@ public: virtual size_t paramCount() override; virtual ~Warburg() = default; }; + +} diff --git a/main.cpp b/main.cpp index 417fa70860ee94d32a62cbe55775e3796e7faf82..dbb50f97865ce81620f328a500084c432b55939b 100644 --- a/main.cpp +++ b/main.cpp @@ -3,80 +3,107 @@ #include <chrono> #include "model.h" +#include "log.h" +#include "options.h" +#include "normalize.h" -void runSingle() +static void printComponants(eis::Model& model) { - std::string modelStr("w{20e3}p{1e-7, 0.9}"); - - std::vector<Model::DataPoint> results; - - Model model(modelStr); - std::cout<<"Compnants: \n"; - for(Componant* componant : model.getFlatComponants()) + eis::Log(eis::Log::DEBUG)<<"Compnants:"; + for(eis::Componant* componant : model.getFlatComponants()) { - std::cout<<Componant::getComponantChar(componant)<<"{"; + eis::Log(eis::Log::DEBUG)<<eis::Componant::getComponantChar(componant)<<"{"; for(size_t i = 0; i < componant->paramCount(); ++i) { - std::cout<<componant->getParam()[i]; + eis::Log(eis::Log::DEBUG)<<componant->getParam()[i]; if(i != componant->paramCount()-1) - std::cout<<", "; + eis::Log(eis::Log::DEBUG)<<", "; } - std::cout<<"}\n"; + eis::Log(eis::Log::DEBUG)<<"}"; } +} + +static void paramSweepCb(std::vector<eis::DataPoint>& data, const std::vector<double>& parameters) +{ + static size_t i = 0; + ++i; + if((i & 0x3FF) == 0) + std::cout<<'.'<<std::flush; +} + +static void runSweep(const std::string& modelString, eis::Range omega, bool normalize = false, bool reduce = false) +{ + std::vector<eis::DataPoint> results; + + eis::Model model(modelString); - Model::Range omega(0, 1e6, 50); + printComponants(model); auto start = std::chrono::high_resolution_clock::now(); results = model.executeSweep(omega); auto end = std::chrono::high_resolution_clock::now(); auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start); - for(const Model::DataPoint& res : results) - std::cout<<"omega: "<<res.omega<<" real = "<<res.im.real()<<" im = "<<res.im.imag()<<'\n'; - - Model modelCopy(model); - results = modelCopy.executeSweep(omega); - - for(const Model::DataPoint& res : results) - std::cout<<"omega: "<<res.omega<<" real = "<<res.im.real()<<" im = "<<res.im.imag()<<'\n'; + if(reduce) + { + eis::Log(eis::Log::INFO)<<"reduced normalized results:"; + results = eis::reduceRegion(results); + } + else if(normalize) + { + eis::Log(eis::Log::INFO)<<"normalized results:"; + eis::normalize(results); + } + else + { + eis::Log(eis::Log::INFO)<<"results:"; + } + eis::Log(eis::Log::INFO)<<"omega,real,im"; - std::cout<<"time taken: "<<duration.count()<<" us"<<'\n'; -} + for(const eis::DataPoint& res : results) + std::cout<<res.omega<<','<<res.im.real()<<','<<res.im.imag()<<'\n'; -void sweepCb(std::vector<Model::DataPoint>& data, const std::vector<double>& parameters) -{ - static size_t i = 0; - ++i; - if((i & 0x3FF) == 0) - std::cout<<'.'<<std::flush; + eis::Log(eis::Log::INFO)<<"time taken: "<<duration.count()<<" us"; } -void runSweep() +static void runParamSweep() { + eis::Log(eis::Log::INFO)<<__func__; std::string modelStr("w{20e3}p{1e-7, 0.9}"); - std::vector<Model::DataPoint> results; - Model model(modelStr); + eis::Model model(modelStr); - std::vector<Model::Range> parameters; - parameters.push_back(Model::Range(1e3, 50e3, 100)); - parameters.push_back(Model::Range(1e-7, 20e-7, 100)); - parameters.push_back(Model::Range(0.7, 1.2, 100)); + std::vector<eis::Range> parameters; + parameters.push_back(eis::Range(1e3, 50e3, 100)); + parameters.push_back(eis::Range(1e-7, 20e-7, 100)); + parameters.push_back(eis::Range(0.7, 1.2, 100)); - Model::Range omega(0, 1e6, 25); + eis::Range omega(0, 1e6, 25); auto start = std::chrono::high_resolution_clock::now(); - model.executeParamSweep(parameters, omega, &sweepCb); + model.executeParamSweep(parameters, omega, ¶mSweepCb); auto end = std::chrono::high_resolution_clock::now(); auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start); + std::cout<<std::endl; - std::cout<<"\ntime taken: "<<duration.count()<<" ms"<<'\n'; + eis::Log(eis::Log::INFO)<<"time taken: "<<duration.count()<<" ms"; } int main(int argc, char** argv) { - runSingle(); - runSweep(); + eis::Log::level = eis::Log::INFO; + Config config; + argp_parse(&argp, argc, argv, 0, 0, &config); + + switch(config.mode) + { + case MODE_SWEEP: + runSweep(config.modelStr, config.omegaRange, config.normalize, config.reduce); + break; + case MODE_PARAM_SWEEP: + runParamSweep(); + break; + } return 0; } diff --git a/model.cpp b/model.cpp index da2bf5cf90e2c6d2efb038e9daa538a11d5b8f1e..b099c84a45e9fa2ddd7e8f7930700ba53fe4ed9f 100644 --- a/model.cpp +++ b/model.cpp @@ -7,6 +7,9 @@ #include "constantphase.h" #include "warburg.h" #include "paralellseriel.h" +#include "log.h" + +using namespace eis; size_t Model::opposingBraket(const std::string& str, size_t index, char bracketChar) { @@ -41,15 +44,13 @@ size_t Model::deepestBraket(const std::string& str) Componant *Model::processBrackets(std::string& str, size_t& bracketCounter) { size_t bracketStart = deepestBraket(str); - std::cout<<str<<" bracket start "<<(bracketStart == std::string::npos ? std::string("npos") : std::to_string(bracketStart))<<'\n'; + Log(Log::DEBUG)<<str<<" bracket start "<<(bracketStart == std::string::npos ? std::string("npos") : std::to_string(bracketStart)); if(bracketStart == std::string::npos) { Componant* componant = processBracket(str); if(!componant) - { - std::cout<<"Warning: can not create componant type B for "<<str<<'\n'; - } + Log(Log::DEBUG)<<"Warning: can not create componant type B for "<<str; return componant; } @@ -65,7 +66,7 @@ Componant *Model::processBrackets(std::string& str, size_t& bracketCounter) Componant* componant = processBracket(bracket); if(!componant) { - std::cout<<"Warning: can not create componant type A for "<<bracket<<'\n'; + Log(Log::DEBUG)<<"can not create componant type A for "<<bracket; } _bracketComponants.push_back(componant); @@ -79,30 +80,30 @@ std::string Model::getParamStr(const std::string& str, size_t index) { if(str.size()-index < 3 || str[index+1] != '{') { - std::cout<<"Warning: missing parameter for "<<str[index]<<'\n'; + Log(Log::DEBUG)<<"missing parameter for "<<str[index]; return 0; } size_t end = opposingBraket(str, index, '}'); std::string parameterStr = str.substr(index+2, end-index-2); - std::cout<<"param for "<<str[index]<<' '<<parameterStr<<'\n'; + Log(Log::DEBUG)<<"param for "<<str[index]<<' '<<parameterStr; return parameterStr; } Componant *Model::processBracket(std::string& str) { - std::cout<<__func__<<'('<<str<<")\n"; + Log(Log::DEBUG)<<__func__<<'('<<str<<')'; std::vector<std::string> tokens = tokenize(str, '-', '{', '}'); std::vector<Componant*> nodes; for(const std::string& nodeStr : tokens) { - std::cout<<__func__<<" full node str: "<<nodeStr<<'\n'; + Log(Log::DEBUG)<<__func__<<" full node str: "<<nodeStr; std::vector<Componant*> componants; for(size_t i = 0; i < nodeStr.size(); ++i) { - std::cout<<__func__<<" arg: "<<nodeStr[i]<<'\n'; + Log(Log::DEBUG)<<__func__<<" arg: "<<nodeStr[i]; switch(nodeStr[i]) { case 'c': @@ -128,7 +129,7 @@ Componant *Model::processBracket(std::string& str) case '{': i = opposingBraket(nodeStr, i, '}'); case '}': - std::cout<<"Warning: stray "<<nodeStr[i]<<" in model string\n"; + Log(Log::WARN)<<"stray "<<nodeStr[i]<<" in model string"; break; case '0' ... '9': { @@ -146,7 +147,7 @@ Componant *Model::processBracket(std::string& str) else if(componants.size() == 1) nodes.push_back(componants[0]); else - std::cout<<"Warning: empty node for "<<nodeStr<<'\n'; + Log(Log::WARN)<<"empty node for "<<nodeStr; } if(nodes.size() > 1) @@ -184,7 +185,7 @@ Model::~Model() delete _model; } -Model::DataPoint Model::execute(double omega) +DataPoint Model::execute(double omega) { if(_model) { @@ -195,7 +196,7 @@ Model::DataPoint Model::execute(double omega) } else { - std::cout<<"Warning: model not ready\n"; + Log(Log::WARN)<<"model not ready"; } return DataPoint({std::complex<double>(0,0), 0}); } @@ -230,24 +231,42 @@ std::vector<Componant*> Model::getFlatComponants() return getFlatComponants(); } -std::vector<Model::DataPoint> Model::executeSweep(const Range& omega) +std::vector<DataPoint> Model::executeSweep(const Range& omega) { - double step = (omega.end - omega.start)/omega.count; - double currOmega = omega.start; - std::vector<DataPoint> results; results.reserve(omega.count); - for(size_t i = 0; i < omega.count; ++i) + + if(!omega.log) + { + double currOmega = omega.start; + double step = (omega.end - omega.start)/(omega.count-1); + + for(size_t i = 0; i < omega.count; ++i) + { + results.push_back(execute(currOmega)); + currOmega+=step; + } + } + else { - results.push_back(execute(currOmega)); - currOmega+=step; + double start = log10(omega.start); + double end = log10(omega.end); + double step = (end-start)/(omega.count-1); + double currOmegaL = start; + + for(size_t i = 0; i < omega.count; ++i) + { + results.push_back(execute(pow(10, currOmegaL))); + currOmegaL+=step; + } } return results; } -std::vector<Model::DataPoint> Model::executeParamByIndex(const std::vector<Range>& componantRanges, const Range& omega, size_t index) +std::vector<DataPoint> Model::executeParamByIndex(const std::vector<Range>& componantRanges, const Range& omega, size_t index) { - assert(setFlatParameters(getSweepParamByIndex(componantRanges, index))); + std::vector<double> parameters = getSweepParamByIndex(componantRanges, index); + assert(setFlatParameters(parameters)); return executeSweep(omega); } @@ -273,7 +292,7 @@ std::vector<double> Model::getSweepParamByIndex(const std::vector<Range>& compon std::vector<double> parameters(parametersCount, 0); for(size_t i = 0; i < parametersCount; ++i) - parameters[i] = parameterIndexies[i]*componantRanges[i].stepSize(); + parameters[i] = parameterIndexies[i]*componantRanges[i].stepSize()+componantRanges[i].start; return parameters; } @@ -283,7 +302,7 @@ bool Model::executeParamSweep(const std::vector<Range>& componantRanges, const R size_t parametersCount = getFlatParametersCount(); if(componantRanges.size() != parametersCount) { - std::cout<<"Error: a parameter range must be provided for eatch componant parameter\n"; + Log(Log::ERROR)<<"a parameter range must be provided for eatch componant parameter"; return false; } @@ -291,12 +310,12 @@ bool Model::executeParamSweep(const std::vector<Range>& componantRanges, const R { if(componantRanges[i].count == 0 || (componantRanges[i].count < 2 && componantRanges[i].start != componantRanges[i].end)) { - std::cout<<"Error: paramter range must specify at least one paramter point if only one paramer point is specified start and end must be the same\n"; + Log(Log::ERROR)<<"paramter range must specify at least one paramter point if only one paramer point is specified start and end must be the same"; return false; } else if(componantRanges[i].start > componantRanges[i].end) { - std::cout<<"Error: paramter range end-start must be positive\n"; + Log(Log::ERROR)<<"paramter range end-start must be positive"; return false; } } @@ -307,7 +326,7 @@ bool Model::executeParamSweep(const std::vector<Range>& componantRanges, const R for(size_t i = 0; i < parametersCount; ++i) currentParam[i] = componantRanges[i].start; - std::cout<<"Executing sweep. Steps requried: "<<stepsRequired<<std::endl; + Log(Log::INFO)<<"Executing sweep. Steps requried: "<<stepsRequired; for(size_t i = 0; i < stepsRequired; ++i) { @@ -336,6 +355,18 @@ bool Model::setFlatParameters(const std::vector<double>& parameters) return true; } +std::vector<double> Model::getFlatParameters() +{ + std::vector<double> params; + std::vector<Componant*> componants = getFlatComponants(); + for(Componant* componant : componants) + { + for(double param : componant->getParam()) + params.push_back(param); + } + return params; +} + size_t Model::getFlatParametersCount() { size_t count = 0; @@ -343,4 +374,3 @@ size_t Model::getFlatParametersCount() count += componant->paramCount(); return count; } - diff --git a/normalize.cpp b/normalize.cpp new file mode 100644 index 0000000000000000000000000000000000000000..e06fc4ec3bf857c1d8ebb13f22a4f73dc8b5d13f --- /dev/null +++ b/normalize.cpp @@ -0,0 +1,108 @@ +#include "normalize.h" +#include <cmath> +#include <complex> +#include <limits> +#include "eistype.h" +#include "log.h" + +std::complex<fvalue> eis::absGrad(const std::vector<eis::DataPoint>& data, size_t index) +{ + if(data.size() < 3) + return std::complex<fvalue>(1,1); + + if(index == 0) + index = 1; + else if(index > data.size()-2) + index = data.size()-2; + + return std::complex<fvalue>(std::abs((data[index+1].im.real()-data[index-1].im.real())/(data[index+1].omega-data[index-1].omega)), + std::abs((data[index+1].im.imag()-data[index-1].im.imag())/(data[index+1].omega-data[index-1].omega))); +} + + + +void eis::eraseSingularites(std::vector<eis::DataPoint>& data) +{ + for(size_t i = 0; i < data.size(); ++i) + { + if(std::isnan(data[i].im.real()) || std::isnan(data[i].im.imag())) + { + size_t left = i-1; + size_t right = i+1; + + while(left > 1 && std::abs(eis::absGrad(data, left)) > 10) + --left; + while(right < data.size()-2 && std::abs(eis::absGrad(data, right)) > 10) + ++right; + + std::complex<fvalue> mean = (data[left].im + data[right].im)/2.0; + for(size_t j = left; j < right; ++j) + data[j].im = mean; + } + } +} + +void eis::normalize(std::vector<eis::DataPoint>& data) +{ + fvalue max = std::numeric_limits<fvalue>::min(); + for(const DataPoint& dataPoint : data) + { + fvalue norm = std::abs(dataPoint.im); + if(norm > max) + max = norm; + } + + for(DataPoint& dataPoint : data) + dataPoint.im = dataPoint.im / max; +} + +std::vector<eis::DataPoint> eis::reduceRegion(const std::vector<eis::DataPoint>& inData, fvalue gradThreshFactor) +{ + if(inData.size() < 3) + return inData; + + std::vector<eis::DataPoint> data = inData; + eis::eraseSingularites(data); + eis::normalize(data); + + std::vector<fvalue> grads; + grads.reserve(data.size()); + fvalue meanGrad = 0; + eis::Log(eis::Log::DEBUG)<<"Grads:"; + for(size_t i = 0; i < data.size(); ++i) + { + grads.push_back(std::abs(eis::absGrad(data, i))); + meanGrad += grads.back(); + eis::Log(eis::Log::DEBUG)<<i<<": "<<inData[i].omega<<','<<grads.back(); + } + + meanGrad = meanGrad / grads.size(); + fvalue gradThresh = meanGrad*gradThreshFactor; + + eis::Log(eis::Log::DEBUG)<<"Grad thresh is:"<<','<<gradThresh; + + size_t start = 0; + for(size_t i = 1; i < data.size()-1; ++i) + { + if(grads[i] < gradThresh) + start = i; + else + break; + } + + size_t end = data.size()-1; + for(size_t i = data.size()-1; i > 1; --i) + { + if(grads[i] < gradThresh) + end = i; + else + break; + } + + eis::Log(eis::Log::DEBUG)<<"reduced range "<<start<<'-'<<end; + + data.erase(data.begin(), data.begin()+start); + data.erase(data.begin()+end, data.end()); + + return data; +} diff --git a/paralellseriel.cpp b/paralellseriel.cpp index 041d0f6afb8fed835194af31440bb9274f835134..ffae512b3bafab4e692c0d35badaa0225b4267a4 100644 --- a/paralellseriel.cpp +++ b/paralellseriel.cpp @@ -1,5 +1,7 @@ #include "paralellseriel.h" +using namespace eis; + Parallel::Parallel(std::vector<Componant*> componantsIn): componants(componantsIn) { } diff --git a/resistor.cpp b/resistor.cpp index d81c17f5cdd45d1f82373e4e40d5caf5d1dca19f..58b524a5872db40f4bddb9bfe66a58167ee07788 100644 --- a/resistor.cpp +++ b/resistor.cpp @@ -3,6 +3,10 @@ #include <vector> #include <math.h> +#include "log.h" + +using namespace eis; + Resistor::Resistor(double r): _R(r) {} @@ -11,7 +15,7 @@ Resistor::Resistor(std::string paramStr) std::vector<std::string> tokens = tokenize(paramStr, ','); if(tokens.empty()) { - std::cout<<"Warning: to few parameters in "<<__func__<<" parameter string: "<<paramStr<<'\n'; + Log(Log::WARN)<<"to few parameters in "<<__func__<<" parameter string: "<<paramStr; _R = 1e3; return; } @@ -23,7 +27,7 @@ Resistor::Resistor(std::string paramStr) } catch(const std::invalid_argument& ia) { - std::cout<<"Warning: cant parse parameter in "<<__func__<<" parameter: "<<tokens[0]<<'\n'; + Log(Log::WARN)<<"can't parse parameter in "<<__func__<<" parameter: "<<tokens[0]; _R = 1e3; } } @@ -44,7 +48,7 @@ void Resistor::setParam(const std::vector<double>& param) { if(param.empty()) { - std::cout<<"Warning: invalid parameter list sent to "<<__func__<<'\n'; + Log(Log::WARN)<<"invalid parameter list sent to "<<__func__; return; } diff --git a/warburg.cpp b/warburg.cpp index 96bb1e5b21e98b4b0e04a96f0d62a5cb9e7c7725..6700f9c8e373be60c26011954a54ab97a2588eaa 100644 --- a/warburg.cpp +++ b/warburg.cpp @@ -3,6 +3,10 @@ #include <cstdlib> #include <math.h> +#include "log.h" + +using namespace eis; + Warburg::Warburg(double a): _A(a) { @@ -13,7 +17,7 @@ Warburg::Warburg(std::string paramStr) std::vector<std::string> tokens = tokenize(paramStr, ','); if(tokens.size() < 1) { - std::cout<<"Warning: to few parameters in "<<__func__<<" parameter string: "<<paramStr<<'\n'; + Log(Log::WARN)<<"to few parameters in "<<__func__<<" parameter string: "<<paramStr<<'\n'; return; } else @@ -24,7 +28,7 @@ Warburg::Warburg(std::string paramStr) } catch(const std::invalid_argument& ia) { - std::cout<<"Warning: cant parse parameter in "<<__func__<<" parameter: "<<tokens[0]<<'\n'; + Log(Log::WARN)<<"can't parse parameter in "<<__func__<<" parameter: "<<tokens[0]<<'\n'; } } } @@ -44,7 +48,7 @@ void Warburg::setParam(const std::vector<double>& param) { if(param.size() != paramCount()) { - std::cout<<"Warning: invalid parameter list sent to "<<__func__<<'\n'; + Log(Log::WARN)<<"Warning: invalid parameter list sent to "<<__func__<<'\n'; return; }