<commit_before>#include "qca.hpp"
#include "utilities.hpp"
#include <map>

/*
 * TODO
 * ----
 * 
 * * Dead Cell Systems
 * * functions to measure / output energy spectra (over various parameters)
 * * function to output polarisation (over Vext or for dead cell's P)
 * * polarisation-polarisation response function
 * * set system, parameters from command line or configuration file
 * * nice / meaningful output functions
 * * grand canonical system
 */

class EpsilonLess
{
public:
    EpsilonLess (double e_ = 1E-10)
    : e(e_)
    {}

    bool operator() (double a, double b)
    {
        return b-a > e;
    }
private:
    double e;
};

template<class System>
void printEnergySpectrum (const System& s) 
{
    typedef std::map<double, int, EpsilonLess> myMap;
    typedef myMap::const_iterator mapIt;
    myMap evs;
    for (int i=0; i<s.H.eigenvalues.size(); i++)
        evs[s.H.eigenvalues(i)]++;
    for (mapIt i=evs.begin(); i!=evs.end(); i++)
        std::cout << i->first << "\t" << i->second << std::endl;
}

template<class System>
void printPolarisationPolarisation (System& s, size_t i, size_t j, 
                                    double beta, const std::vector<double>& Vexts)
{
    for (size_t k=0; k<Vexts.size(); k++)
    {
        s.H.Vext = Vexts[k];
        s.H.construct();
        s.H.diagonalise();
        std::cout
            << s.ensembleAverage(beta, s.P(i)) << "\t"
            << s.ensembleAverage(beta, s.P(j)) << std::endl;
    }
}

CommandLineOptions setupCLOptions ()
{
    CommandLineOptions o;
    /*
    o.add("help", "h", "Print this help message.")
     .add("model", "m", "Which QCA model to use.")
     .add("N_p", "p", "Number of plaquets.")
     .add("hopping", "t", "Hopping parameter.")
     .add("hopping-diagonal", "t_d", "Diagonal hopping parameter.")
     .add("on-site", "V_0", "On-site coulomb repulsion (Hubbard U).")
     .add("external", "V_ext", "External potential.")
     .add("intra-plaquet-spacing", "a", "Intra-plaquet spacing.")
     .add("inter-plaquet-spacing", "b", "Inter-plaquet spacing.");
     */
    o.add("help", "h", "Print this help message.")
     .add("model", "m", "Which QCA model to use.")
     .add("N_p", "p", "Number of plaquets.")
     .add("t", "t", "Hopping parameter.")
     .add("td", "td", "Diagonal hopping parameter.")
     .add("V0", "V0", "On-site coulomb repulsion (Hubbard U).")
     .add("Vext", "Vext", "External potential.")
     .add("a", "a", "Intra-plaquet spacing.")
     .add("b", "b", "Inter-plaquet spacing.")
     .add("beta", "beta", "Inverse temperature.")
     .add("test-list", "Test list parsing.")
     .add("energy-spectrum", "E", "Calculate the energy spectrum.")
     .add("polarisation", "P", "Calculate the polarisation for specified plaquet(s).");
    
    o["N_p"] = 1;
    o["beta"] = 1;

    return o;
}

void printUsage (CommandLineOptions& o)
{
    std::cerr << "Usage: " << std::endl;
    std::cerr << o.optionsDescription();
}

template<class System>
void run (CommandLineOptions& opts)
{
    std::vector<double> Vexts = opts["Vext"];
    //TODO: opts["Vext"] = 0 does not work
    if (Vexts.size() > 1) opts["Vext"] = "0"; 

    System qca(opts);

    for (size_t j=0; j<Vexts.size(); j++)
    {
        qca.H.Vext = Vexts[j];
    
        qca.H.construct();
        qca.H.diagonalise();

        if (opts["energy-spectrum"])
        {
            printEnergySpectrum(qca);
            std::cout << std::endl;
        }

        /*
         * TODO: -P 0 is a valid specification for the polarisation, but does not
         * work (yields false below) -> fix this in DescriptionItem
         */
        if (opts["polarisation"])
        {
            std::vector<size_t> ps = opts["polarisation"];
            for (size_t i=0; i<ps.size(); i++)
            {
                if (ps[i] >= qca.N_p) 
                {
                    std::cerr << std::endl << "Polarisation: There is no plaquet " 
                        << ps[i] << " in this system." << std::endl;
                    std::exit(EXIT_FAILURE);
                }
                if (i>0) std::cout << "\t";
                std::cout << qca.ensembleAverage(opts["beta"], qca.P(ps[i]));
            }
            std::cout << std::endl;
        }
    }
            
    /*<commit_msg>Improved energy spectrum and polarisation output.<commit_after>#include "qca.hpp"
#include "utilities.hpp"
#include <map>

/*
 * TODO
 * ----
 * 
 * * Dead Cell Systems
 * * functions to measure / output energy spectra (over various parameters)
 * * function to output polarisation (over Vext or for dead cell's P)
 * * polarisation-polarisation response function
 * * set system, parameters from command line or configuration file
 * * nice / meaningful output functions
 * * grand canonical system
 */

class EpsilonLess
{
public:
    EpsilonLess (double e_ = 1E-10)
    : e(e_)
    {}

    bool operator() (double a, double b)
    {
        return b-a > e;
    }
private:
    double e;
};

template<class System>
void printEnergySpectrum (const System& s) 
{
    typedef std::map<double, int, EpsilonLess> myMap;
    typedef myMap::const_iterator mapIt;
    myMap evs;
    for (int i=0; i<s.H.eigenvalues.size(); i++)
        evs[s.H.eigenvalues(i)]++;
    for (mapIt i=evs.begin(); i!=evs.end(); i++)
        std::cout << i->first << "\t" << i->first - s.H.Emin << "\t" << i->second << std::endl;
}

template<class System>
void printPolarisationPolarisation (System& s, size_t i, size_t j, 
                                    double beta, const std::vector<double>& Vexts)
{
    for (size_t k=0; k<Vexts.size(); k++)
    {
        s.H.Vext = Vexts[k];
        s.H.construct();
        s.H.diagonalise();
        std::cout
            << s.ensembleAverage(beta, s.P(i)) << "\t"
            << s.ensembleAverage(beta, s.P(j)) << std::endl;
    }
}

CommandLineOptions setupCLOptions ()
{
    CommandLineOptions o;
    /*
    o.add("help", "h", "Print this help message.")
     .add("model", "m", "Which QCA model to use.")
     .add("N_p", "p", "Number of plaquets.")
     .add("hopping", "t", "Hopping parameter.")
     .add("hopping-diagonal", "t_d", "Diagonal hopping parameter.")
     .add("on-site", "V_0", "On-site coulomb repulsion (Hubbard U).")
     .add("external", "V_ext", "External potential.")
     .add("intra-plaquet-spacing", "a", "Intra-plaquet spacing.")
     .add("inter-plaquet-spacing", "b", "Inter-plaquet spacing.");
     */
    o.add("help", "h", "Print this help message.")
     .add("model", "m", "Which QCA model to use.")
     .add("N_p", "p", "Number of plaquets.")
     .add("t", "t", "Hopping parameter.")
     .add("td", "td", "Diagonal hopping parameter.")
     .add("V0", "V0", "On-site coulomb repulsion (Hubbard U).")
     .add("Vext", "Vext", "External potential.")
     .add("a", "a", "Intra-plaquet spacing.")
     .add("b", "b", "Inter-plaquet spacing.")
     .add("beta", "beta", "Inverse temperature.")
     .add("test-list", "Test list parsing.")
     .add("energy-spectrum", "E", "Calculate the energy spectrum.")
     .add("polarisation", "P", "Calculate the polarisation for specified plaquet(s).");
    
    o["N_p"] = 1;
    o["beta"] = 1;

    return o;
}

void printUsage (CommandLineOptions& o)
{
    std::cerr << "Usage: " << std::endl;
    std::cerr << o.optionsDescription();
}

template<class System>
void run (CommandLineOptions& opts)
{
    std::vector<double> Vexts = opts["Vext"];
    //TODO: opts["Vext"] = 0 does not work
    if (Vexts.size() > 1) opts["Vext"] = "0"; 

    System qca(opts);

    for (size_t j=0; j<Vexts.size(); j++)
    {
        qca.H.Vext = Vexts[j];
    
        qca.H.construct();
        qca.H.diagonalise();

        if (opts["energy-spectrum"])
        {
            printEnergySpectrum(qca);
            std::cout << std::endl;
        }

        /*
         * TODO: -P 0 is a valid specification for the polarisation, but does not
         * work (yields false below) -> fix this in DescriptionItem
         */
        if (opts["polarisation"])
        {
            std::cout << qca.H.Vext;
            std::vector<size_t> ps = opts["polarisation"];
            for (size_t i=0; i<ps.size(); i++)
            {
                if (ps[i] >= qca.N_p) 
                {
                    std::cerr << std::endl << "Polarisation: There is no plaquet " 
                        << ps[i] << " in this system." << std::endl;
                    std::exit(EXIT_FAILURE);
                }
                std::cout << "\t" << qca.ensembleAverage(opts["beta"], qca.P(ps[i]));
            }
            std::cout << std::endl;
        }
    }
            
    /*
{'fraction_non_alphanumeric': 0.17844727694090382, 'fraction_numerical': 0.003823870220162225, 'mean_word_length': 2.286747905559787, 'pattern_counts': {'":': 0, '<': 113, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1291369', 'n_tokens_mistral': 3010, 'n_tokens_neox': 2885, 'n_words': 723}