adding in paolos script for calculating fs_strength

#include <iostream>

#include <cstdint>

#include <string>

#include <sstream>

#include <fstream>

#include <vector>

#include <cmath>

using namespace std;

struct OutData {

    double I;

    double V;

    double p;

    OutData(double ii, double vv, double pp) : I(ii), V(vv), p(pp) {}

};

class OutIV { 

    using OutArray = std::vector<OutData>;

    int const realizations;

    OutArray accArray; //avg envelope

    int ppv = 100; //points per volt interval

    OutArray::iterator thit;

public:

    OutIV(int const realizations) : realizations(realizations) {}

    static bool compareV(const OutData& a, const OutData& b) {return (a.V < b.V);}

    double interpolate_y(double x, double x0, double x1, double y0, double y1)

    {

        return y0 + (x-x0)*(y1-y0)/(x1-x0);

    }

    OutArray exctract_envelope(const OutArray& s)

    {

        OutArray enve;

        OutData prev(0.,0.,0.); //last analyzed

        OutData last(0.,0.,0.); //last stored

        enve.push_back(prev);

        bool previous_was_discarded = false;

        for(auto& d: s) {

            if(compareV(last,d)) {

                if(previous_was_discarded) {

                    double new_V = last.V;

                    //double new_I = prev.I+((last.V-prev.V)/(d.V-prev.V))*(d.I-prev.I);

                    double new_I = interpolate_y(last.V,prev.V,d.V,prev.I,d.I);

                    //double new_p = interpolate_y(last.V,prev.V,d.V,prev.p,d.p);

                    double new_p = prev.p; //not interpolated for the single realization! int value

                    enve.emplace_back(new_I,new_V,new_p);               

                }

                enve.push_back(d);

                last = d;

                previous_was_discarded = false;

            }

            else previous_was_discarded = true;

            prev = d;

        }

        enve.emplace_back(0.,last.V,last.p); //last point

        return enve;    

    }

    void export_single_envelope(const OutArray& ar, const char* filename)

    {

        std::ofstream f;

        f.open(filename);

        for(auto& d: ar)  {

            f << d.V << " " << d.I << " "

            << d.p << std::endl;

        }

        f.close();

    }

    OutArray resample_envelope(const OutArray& enve)

    {

        double largest_V = enve.back().V; 

        int sgn_V = static_cast<int>(fabs(largest_V)/largest_V);

        int hat_V = 1+static_cast<int>(ceil(sgn_V*largest_V));

        //hat_V is the largest possible V as integer in abs value

        //ceil guarantees that it is in fact much lager than that in most cases

        //This is because one wants to make sure that "a" has enough elements 

        //to write on, if floating point comparisons lead to unexpected results

        //Yet there are some cases in which during the last step *it->V

        //overcomes hat_V upon an increment of 1./ppv:

        //the +1 becomes essential in these cases, in which otherwise

        //*it would end up out of bounds. As a result, the output array

        //is usually padded by a >= 1-volt interval of zero currents.

        //This last (clearly exaggerate) measure exemplifies 

        //my absolute discomfort with floating point arithmetics

        int target_size = hat_V * ppv;

        OutArray a;

        a.reserve(target_size);

        for(int i=0; i<target_size; ++i) { 

            double v = i / static_cast<double>(ppv);

            a.emplace_back(0.,v,0);

        }

        auto it = a.begin();

        auto fr = enve.begin();

        auto to = enve.begin()+1;

        while(to != enve.end()) {

            while(compareV(*it,*fr)) ++it; //should not be necessary, mostly a guard

            while(compareV(*it,*to)) {

                it->I = interpolate_y(it->V,fr->V,to->V,fr->I,to->I); 

                it->p = interpolate_y(it->V,fr->V,to->V,fr->p,to->p);

                ++it;

            }

            ++fr; ++to;

        }

        //Next point: I=0, p same as last in unsampled data

        it->p = enve.back().p;

        //Next point to end: remove everything

        a.erase(it+1,a.end());

        return a; //Make sure RVO or move semantics are ON!

    }

    void accumulate_envelope(const OutArray& s, bool optional) {

        OutArray enve = exctract_envelope(s);

        if(optional) export_single_envelope(enve, "envelope"); 

        OutArray resampled = resample_envelope(enve);

        const int old_size = accArray.size();

        const int resampled_size = resampled.size();

        //If necessary, resize the target array:

        if(old_size<resampled_size) {

            accArray.reserve(resampled_size);

            const double last_p = accArray.back().p;

            for(int i=0; i<resampled_size-old_size; ++i) { 

                const double v = (old_size+i) / ((double)ppv);

                accArray.emplace_back(0.,v,last_p);

            }

        }

        //Accumulate resampled data into target array:

        const int new_size = accArray.size();//it is at least resampled_size

        for(int i=0; i<resampled_size; ++i) {

            accArray[i].I += resampled[i].I;

            accArray[i].p += resampled[i].p;

        }

        //In case resampled data have smaller support,

        //accumulate equivalent data for the missing support interval:

        //(current=0, broken link density=constant)

        if(resampled_size<new_size) {

            const double last_p = resampled.back().p;

            for(int i=resampled_size; i<new_size; ++i) {

                //accArray[i].I += 0; left here for clarity

                accArray[i].p += last_p;

            }

        }   

    }

    void export_average_envelope()

    {

        std::ofstream f;

        const char* filename = "average";

        f.open(filename);

        double norm = 1.0/realizations;

        for(auto d: accArray)  {

            f << d.V << " " << (d.I)*norm << " " << (d.p)*norm << "\n";

        }

        f.close();  

    }

};

struct IVstorage {

    uint64_t ncurves = 0;

    std::vector<float> i;

    std::vector<float> v;

    std::vector<uint64_t> ptrs;

    IVstorage() {

        i.reserve(150000000);

        v.reserve(150000000);

        ptrs.reserve(200000); 

    }

    void import_iv_from_ivcat_to_storage(std::string filename) {

        uint64_t ptr = 0;

        std::ifstream f;

        f.open(filename);

        std::string line;

        while(std::getline(f, line)) {

            if(line.front() == '#') {

                ++ncurves;

                ptrs.push_back(ptr);

            }

            else {

                std::istringstream iss(line);

                double ii;

                double vv;

                if(!(iss >> vv >> ii)) break;

                i.push_back(ii);

                v.push_back(vv);

                ++ptr;

            }

        }

        f.close();

        ptrs.push_back(ptr);

    }

};

class ToughnessCounter {

    std::string filename;

    uint64_t ncurves;

    std::vector<double> ip; // peak i

    //std::vector<double> vp; // v for peak i

    std::vector<double> tp; // the actual toughness

    IVstorage& ivs;

    uint32_t Nx;

public:   

    ToughnessCounter(IVstorage& ivs, uint32_t Nx) : ivs(ivs), Nx(Nx)

        {ncurves = ivs.ncurves;} 

    void accumulate_true_toughness_from_range(uint64_t fr, uint64_t to) {

        std::vector<OutData> od;

        for(uint64_t point=fr; point<to; ++point) {

            od.emplace_back(ivs.i[point],ivs.v[point],point-fr);

        }

        OutIV out(1);

        std::vector<OutData> enve = out.exctract_envelope(od);

        //out.export_single_envelope(enve,"envelope");

        double toughness = 0;

        OutData prev(0.,0.,0.);

        for(uint32_t pp=1; pp<enve.size(); ++pp) {

            OutData curr = enve[pp];

            double deltaV = curr.V-prev.V;

            toughness += 0.5*(curr.I+prev.I)*deltaV;

            prev = curr; 

        }

        tp.push_back(toughness);

    }

    void accumulate_strength_from_range(uint64_t fr, uint64_t to) {

        std::vector<OutData> od;

        for(uint64_t point=fr; point<to; ++point) {

            od.emplace_back(ivs.i[point],ivs.v[point],point-fr);

        }

        OutIV out(1);

        std::vector<OutData> enve = out.exctract_envelope(od);

        //out.export_single_envelope(enve,"envelope");

        double toughness = 0;

        double candidate_ip = 0;

        OutData prev(0.,0.,0.);

        for(uint32_t pp=1; pp<enve.size(); ++pp) {

            OutData curr = enve[pp];

            double deltaV = curr.V-prev.V;

            toughness += 0.5*(curr.I+prev.I)*deltaV;

            prev = curr; 

            //now the peak current

            if(curr.I>candidate_ip) {

                candidate_ip = curr.I; 

            }

        }

        tp.push_back(toughness);

        ip.push_back(candidate_ip);

    }

    void accumulate_strength() {

        for(uint64_t n=0; n<ncurves; ++n) {

            //std::cout << n << std::endl;

            uint64_t fr = ivs.ptrs[n];

            uint64_t to = ivs.ptrs[n+1];

            accumulate_strength_from_range(fr, to);

        }

    }

    void accumulate_true_toughness() {

        for(uint64_t n=0; n<ncurves; ++n) {

            //std::cout << n << std::endl;

            uint64_t fr = ivs.ptrs[n];

            uint64_t to = ivs.ptrs[n+1];

            accumulate_true_toughness_from_range(fr, to);

        }

    }

    void export_strength(const char* filename) 

    {

        double tpavg = 0;

        double tp2avg = 0;

        std::size_t npoints = tp.size();

        for(uint64_t i=0; i<npoints; ++i) {

            tpavg += tp[i];

            tp2avg += tp[i]*tp[i];

        }

        tpavg = tpavg/npoints;

        tp2avg = tp2avg/npoints;

        double tpsigma = sqrt( tp2avg-tpavg*tpavg );

        double ipavg = 0;

        double ip2avg = 0;

        npoints = ip.size();

        for(uint64_t i=0; i<npoints; ++i) {

            ipavg += ip[i];

            ip2avg += ip[i]*ip[i];

        }

        ipavg = ipavg/npoints;

        ip2avg = ip2avg/npoints;

        double ipsigma = sqrt( ip2avg-ipavg*ipavg );

        std::ofstream f;

        f.open(filename);

        f << ipavg << " " << ipsigma << "\n";

        f << tpavg << " " << tpsigma << "\n";

        f.close();

    }

    void export_true_toughness(const char* filename) 

    {

        double tpavg = 0;

        std::size_t npoints = tp.size();

        for(uint64_t i=0; i<npoints; ++i) {

            tpavg += tp[i];

        }

        tpavg = tpavg/npoints;

        std::ofstream f;

        f.open(filename);

        f << tpavg << "\n";

        f.close();

    }

};

template<typename T>

void get_from_command_line(int argc, char **argv, std::string name, T &a) {

    std::string value;

    for(int i=1; i<argc; ++i) {

        std::string argument(argv[i]);

        size_t found = argument.find(name);

        if(found!=std::string::npos) {

            size_t eqsgn = argument.find("=");

            value = std::string(argument,eqsgn+1,argument.size());

            break;

        }

    }

    std::stringstream vss(value);

    vss >> a;

}

int main(int argc, char **argv) 

{    

    uint32_t Nx = 512;

    get_from_command_line<>(argc,argv,"Nx",Nx);

    std::cout << "Nx=" << Nx << std::endl;

    IVstorage ivs;

    ivs.import_iv_from_ivcat_to_storage("ivcat");

    //ivs.import_iv_from_ivcat_to_storage("ivcat_steps9");

    ToughnessCounter ac(ivs, Nx);

    ac.accumulate_strength();

    ac.export_strength("strength");

}