diff --git a/main.cpp b/main.cpp index b4d9f1d20ae4e157ccbfaace901df6449830d17a..d7a66a06d8d0eb51c27fecfd42096b4c03da3e95 100644 --- a/main.cpp +++ b/main.cpp @@ -16,7 +16,7 @@ #define M_PI 3.14159265358979323846 #endif -static constexpr double DIST_THRESH = 0.01; +static constexpr double STEP_THRESH = 0.30; static constexpr char PARA_SWEEP_OUTPUT_DIR[] = "./sweep"; static void printComponants(eis::Model& model) @@ -52,7 +52,8 @@ static void runSweep(const Config& config, eis::Model& model) if(config.reduce) { eis::Log(eis::Log::INFO)<<"reduced normalized results:"; - results = eis::reduceRegion(results); + results = eis::reduceRegion(results, 0.01, false); + results = eis::reduceRegion(results, 0.01, true); } else if(config.normalize) { @@ -67,7 +68,7 @@ static void runSweep(const Config& config, eis::Model& model) if(config.noise > 0) eis::noise(results, config.noise, false); - std::cout<<(config.hertz ? "freqency" : "omega")<<",real,im\n"; + eis::Log(eis::Log::INFO)<<(config.hertz ? "freqency" : "omega")<<",real,im"; for(const eis::DataPoint& res : results) std::cout<<(config.hertz ? res.omega/(2*M_PI) : res.omega)<<','<< @@ -97,9 +98,32 @@ static void runParamSweep(const Config& config, eis::Model& model) data = allSweeps[i]; else data = model.executeSweep(config.omegaRange, i); - size_t outputSize = data.size(); - data = eis::reduceRegion(data); - data = eis::rescale(data, outputSize); + if(config.normalize) + eis::normalize(data); + if(config.reduce) + { + size_t initalDataSize = data.size(); + data = eis::reduceRegion(data); + if(data.size() < initalDataSize/8) + { + eis::Log(eis::Log::INFO)<<"\nskipping output for step "<<i + <<" as data has no interesting region"; + continue; + } + //data = eis::rescale(data, initalDataSize); + } + + if(config.skipLinear && i > 0) + { + fvalue correlation = std::abs(pearsonCorrelation(data)); + if(correlation > 0.5) + { + eis::Log(eis::Log::INFO)<<"skipping output for step "<<i + <<" as data is too linear: "<<correlation; + continue; + } + } + eis::saveToDisk(data, std::string(PARA_SWEEP_OUTPUT_DIR)+std::string("/")+std::to_string(i)+".csv", model.getModelStrWithParam(i)); eis::Log(eis::Log::INFO, false)<<'.'; } @@ -110,25 +134,18 @@ static void runParamSweep(const Config& config, eis::Model& model) eis::Log(eis::Log::INFO)<<"time taken: "<<duration.count()<<" ms"; } -void findRanges(const Config& config, eis::Model& model) +static std::vector<std::vector<fvalue>> getRangeValuesForModel(const Config& config, eis::Model& model, std::vector<size_t>* indices = nullptr) { std::vector<std::vector<fvalue>> values; std::vector<std::vector<eis::DataPoint>> sweeps; size_t count = model.getRequiredStepsForSweeps(); + eis::Log(eis::Log::INFO)<<"Executeing "<<count<<" steps"; std::vector<eis::Componant*> componants = model.getFlatComponants(); std::vector<std::vector<eis::DataPoint>> allSweeps; if(config.threaded) allSweeps = model.executeAllSweeps(config.omegaRange); - std::vector<char> componantChars; - for(eis::Componant* componant : componants) - { - std::vector<eis::Range>& ranges = componant->getParamRanges(); - for(size_t i = 0; i < ranges.size(); ++i) - componantChars.push_back(componant->getComponantChar()); - } - for(size_t i = 0; i < count; ++i) { std::vector<eis::DataPoint> data; @@ -136,14 +153,29 @@ void findRanges(const Config& config, eis::Model& model) data = allSweeps[i]; else data = model.executeSweep(config.omegaRange, i); - size_t outputSize = data.size(); - data = eis::reduceRegion(data); - data = eis::rescale(data, outputSize); + normalize(data); + + fvalue maxJump = maximumNyquistJump(data); + if(maxJump > STEP_THRESH) + { + eis::Log(eis::Log::DEBUG)<<"skipping output for step "<<i + <<" is not well centered: "<<maxJump; + continue; + } + + fvalue correlation = std::abs(pearsonCorrelation(data)); + if(correlation > 0.8) + { + eis::Log(eis::Log::DEBUG)<<"skipping output for step "<<i + <<" as data is too linear: "<<correlation; + continue; + } bool found = false; - for(ssize_t i = static_cast<ssize_t>(sweeps.size())-1; i >= 0; --i) + for(ssize_t j = static_cast<ssize_t>(sweeps.size())-1; j >= 0; --j) { - if(eisDistance(data, sweeps[i]) < DIST_THRESH) + fvalue dist = eisDistance(data, sweeps[j]); + if(dist < config.rangeDistance) { found = true; break; @@ -151,6 +183,8 @@ void findRanges(const Config& config, eis::Model& model) } if(!found) { + if(indices) + indices->push_back(i); sweeps.push_back(data); values.push_back(std::vector<fvalue>()); if(config.threaded) @@ -170,6 +204,20 @@ void findRanges(const Config& config, eis::Model& model) } eis::Log(eis::Log::INFO, false)<<'\n'; + return values; +} + +static void findRanges(const Config& config, eis::Model& model) +{ + std::vector<std::vector<fvalue>> values = getRangeValuesForModel(config, model); + std::vector<eis::Componant*> componants = model.getFlatComponants(); + + if(values.empty()) + { + eis::Log(eis::Log::INFO)<<"Cant recommend a range"; + return; + } + std::vector<fvalue> maxValues(values[0].size(), std::numeric_limits<fvalue>::min()); std::vector<fvalue> minValues(values[0].size(), std::numeric_limits<fvalue>::max()); for(size_t i = 0; i < values.size(); ++i) @@ -183,11 +231,45 @@ void findRanges(const Config& config, eis::Model& model) } } - eis::Log(eis::Log::INFO)<<"Recommended ranges:"; - for(size_t i = 0; i < maxValues.size(); ++i) + size_t i = 0; + for(eis::Componant* componant : componants) + { + std::vector<eis::Range>& ranges = componant->getParamRanges(); + for(eis::Range& range : ranges) + { + fvalue mean = (maxValues[i]+minValues[i])/2; + fvalue var = std::pow(maxValues[i]-minValues[i], 2)/std::pow(mean, 2); + if(var < 0.1) + { + range.start = mean; + range.end = mean; + range.count = 1; + } + else + { + range.start = minValues[i]; + range.end = maxValues[i]; + } + ++i; + } + } + + eis::Log(eis::Log::INFO)<<"Recommended ranges:\n"<<model.getModelStrWithParam(); +} + +static void outputRanges(const Config& config, eis::Model& model) +{ + std::vector<size_t> indices; + getRangeValuesForModel(config, model, &indices); + + if(indices.empty()) { - eis::Log(eis::Log::INFO)<<componantChars[i]<<": "<<minValues[i]<<'-'<<maxValues[i]; + eis::Log(eis::Log::INFO)<<"Cant recommend a range"; + return; } + + for(size_t index : indices) + std::cout<<model.getModelStrWithParam(index)<<'\n'; } int main(int argc, char** argv) @@ -212,23 +294,51 @@ int main(int argc, char** argv) else if(config.inputType == INPUT_TYPE_UNKOWN) eis::Log(eis::Log::WARN)<<"Invalid input type specified, assumeing eis"; - eis::Log(eis::Log::INFO)<<"Using model string: "<<config.modelStr; + if(config.inputType != INPUT_TYPE_EIS) + eis::Log(eis::Log::INFO)<<"Using translated model string: "<<config.modelStr; - eis::Model model(config.modelStr, config.paramSteps); - printComponants(model); + try + { + eis::Model model(config.modelStr, config.paramSteps); + + eis::Log(eis::Log::INFO)<<"Using parsed model string: "<<model.getModelStrWithParam(); - if(config.findRange && !model.isParamSweep()) + printComponants(model); + + if(config.mode == MODE_INVALID) + { + eis::Log(eis::Log::ERROR)<<"Invalid mode selected"; + return 1; + } + + if(config.mode == MODE_FIND_RANGE && !model.isParamSweep()) + { + eis::Log(eis::Log::ERROR)<<"Cant find range on a non-sweep model"; + return 1; + } + + if(config.mode == MODE_FIND_RANGE) + { + findRanges(config, model); + } + if(config.mode == MODE_OUTPUT_RANGE_DATAPOINTS) + { + outputRanges(config, model); + } + else + { + if(model.isParamSweep()) + runParamSweep(config, model); + else + runSweep(config, model); + } + + } + catch(const std::invalid_argument& ia) { - eis::Log(eis::Log::ERROR)<<"Cant find range on a non-sweep model"; + eis::Log(eis::Log::ERROR)<<"Unable to parse model string, "<<ia.what(); return 1; } - if(model.isParamSweep() && config.findRange) - findRanges(config, model); - else if(model.isParamSweep()) - runParamSweep(config, model); - else - runSweep(config, model); - return 0; } diff --git a/options.h b/options.h index a58bb650da5fdc64f193a6f2a3e5106a77817d50..cb5fcac8795e458678b4aa4552ad628a30304382 100644 --- a/options.h +++ b/options.h @@ -26,8 +26,10 @@ static struct argp_option options[] = {"invert", 'i', 0, 0, "inverts the imaginary axis"}, {"noise", 'x', "[AMPLITUDE]", 0, "add noise to output"}, {"input-type", 't', "[STRING]", 0, "set input string type, possible values: eis, boukamp, relaxis"}, - {"find-range", 'f', 0, 0, "find sensible part of parameter range"}, + {"find-range", 'f', "[STRING]", 0, "mode, possible values: export, find-range, export-ranges"}, + {"range-distance", 'd', "[DISTANCE]", 0, "distance from a previous point where a range is considered \"new\""}, {"parallel", 'p', 0, 0, "run on multiple threads"}, + {"skip-linear", 'e', 0, 0, "dont output param sweeps that create linear nyquist plots"}, { 0 } }; @@ -39,19 +41,29 @@ enum INPUT_TYPE_UNKOWN }; +enum +{ + MODE_NORMAL, + MODE_FIND_RANGE, + MODE_OUTPUT_RANGE_DATAPOINTS, + MODE_INVALID, +}; + struct Config { std::string modelStr = "c{1e-6}r{1e3}-r{1e3}"; int inputType = INPUT_TYPE_EIS; + int mode = MODE_NORMAL; size_t paramSteps = 10; eis::Range omegaRange; bool normalize = false; bool reduce = false; bool hertz = false; bool invert = false; - bool findRange = false; bool threaded = false; + bool skipLinear = false; double noise = 0; + double rangeDistance = 0.35; Config(): omegaRange(1, 1e6, 50, true) {} @@ -68,6 +80,17 @@ static int parseInputType(const std::string& str) return INPUT_TYPE_UNKOWN; } +static int parseMode(const std::string& str) +{ + if(str == "export") + return INPUT_TYPE_EIS; + else if(str == "find-range") + return MODE_FIND_RANGE; + else if(str == "export-ranges") + return MODE_OUTPUT_RANGE_DATAPOINTS; + return MODE_INVALID; +} + static error_t parse_opt (int key, char *arg, struct argp_state *state) { @@ -148,11 +171,17 @@ parse_opt (int key, char *arg, struct argp_state *state) config->inputType = parseInputType(std::string(arg)); break; case 'f': - config->findRange = true; + config->mode = parseMode(std::string(arg)); break; case 'p': config->threaded = true; break; + case 'e': + config->skipLinear = true; + break; + case 'd': + config->rangeDistance = std::stod(std::string(arg)); + break; default: return ARGP_ERR_UNKNOWN; } diff --git a/strops.cpp b/strops.cpp index 3c1bc0fcfb220ad6a679c0c71b71a84ee8f5b820..67c726fb99171745627388d285ed8500863a42c8 100644 --- a/strops.cpp +++ b/strops.cpp @@ -104,6 +104,16 @@ size_t eisRemoveUnneededBrackets(std::string& in, long int bracketStart) { bool bracketNeeded = false; size_t paramBracketCount = 0; + + if(bracketStart == -1 && in.size() > 2) + { + if(in[0] == '(' && in.back() == ')') + { + in.erase(in.begin()); + in.pop_back(); + } + } + for(size_t i = (bracketStart >= 0 ? bracketStart+1 : 0); i < in.size(); ++i) { if(paramBracketCount == 0)