#include <iostream>
#include <fstream>
#include <sstream>
using namespace std;
#include "CausalStateAnalyzer.h"

#include <stdlib.h>

// For past and future cones. Cell states at t, then at t+-1, t+-2, etc...
struct LightCone : public vector<int> {

    // map cones to id so distributions only manipulate the id
    // distributions are possible with the cones directly too, thanks to templatized
    // analyzer, but this is faster
    typedef map<LightCone,int> LightConeMap;

    // Implement a lesser_than operator to get O(log(n)) insert on maps
    bool operator<(const LightCone& cone) {
        // use simple lexicographic ordering
        const int minsize = min(size(), cone.size());
        for (int i=0; i<minsize; ++i) {
            if ((*this)[i] < cone[i]) return true;
            if ((*this)[i] > cone[i]) return false;
        }
        // elements are equal so far, strict < iff second cone has more elements
        return size() < cone.size();
    }

    static LightConeMap conesId;

    // returns the id of this light cone, make a new one if necessary the first time
    int getId() {
        LightConeMap::iterator myPlace = conesId.find(*this);
        if (myPlace != conesId.end()) return myPlace->second;
        return conesId[*this] = ++nextId; // start from 1
    }

    static int nextId;
};
int LightCone::nextId = 0;
LightCone::LightConeMap LightCone::conesId;

struct Rule {

    int nextState[8];

    Rule(unsigned int number) {
        // convert CA rule notation to state array
        for (int i=0; i<8; ++i) {
            nextState[i] = number & 1;
            number = number >> 1;
        }
    }

    // apply the rule to the 3 cells that are either 0 or 1
    int apply(int left, int center, int right) {
        return nextState[ (left << 2) | (center << 1) | right ];
    }
};

// to generate the initial conditions
struct RandFunctor {
    int operator()() {
        return (rand() >> 17) & 1; // use any bit
    }
};


int main(int argc, char** argv) {

    // User-defined variables
    int past_depth = 5;
    int future_depth = 3;
    int rule_number = -1;
    int seed = 1;
    float siglevel = 0.04321;

    bool quietMode = false;
    bool columnMode = false;
    opterr = 0; int c;
    while ((c=getopt(argc,argv,"chqp:f:s:r:g:"))!=-1) switch(c) {
        case 'q': quietMode = true; break;
        case 'c': columnMode = true; break;
        case 'g': {
            float arg = (float)atof(optarg);
            if (arg!=0.0f) siglevel = arg;
            break;
        }
        case 'p': {
            int arg = atoi(optarg);
            if (arg!=0) past_depth=arg;
            break;
        }
        case 'f': {
            int arg = atoi(optarg);
            if (arg!=0) future_depth=arg;
            break;
        }
        case 's': {
            int arg = atoi(optarg);
            if (arg!=0) seed=arg;
            break;
        }
        case 'r': {
            rule_number=atoi(optarg); // 0 is valid
            break;
        }
        case '?':
        case 'h':
        default:
            cout << "Usage:\n-p<integer>: past depth\n-f<integer>: future depth\n-r<integer>: rule number\n-g<float>: chi-square significance level for clustering distributions\n-q: quiet mode\n-s<integer>: random seed\n-c: Use single column instead of whole cones." << endl;
            return 0;
    }

    // derived variables
    int wsize = future_depth*2-1 + (future_depth+past_depth-2)*2;

    if (rule_number==-1) rule_number = 110; // default

    if (!quietMode) {
        std::ios::sync_with_stdio(false); // GNU specific?
        cout << "past depth = " << past_depth << endl;
        cout << "future depth = " << future_depth << endl;
        cout << "world size = " << wsize << endl;
        cout << "rule number = " << rule_number << endl;
        cout << "random seed = " << seed << endl;
        cout << "significance for clustering = " << siglevel << endl;
        if (columnMode) cout << "warning: column mode" << endl;
    }

    srand(seed);

    CausalStateAnalyzer<int,int> analyzer(siglevel);
    Rule rule(rule_number);

    // double-buffering
    int state1[wsize], state2[wsize];
    int* state = state1;
    int* nextState = state2;

    for (int initial_config = 0; initial_config < (1 << wsize); ++initial_config) {
        int config = initial_config;
        for (int i=0; i<wsize; ++i) state[wsize-1-i] = (config >> i) & 1;

        LightCone plc;
        int center = wsize/2;      // wsize necessarily odd, OK
        int side = past_depth - 1; // spacial extent of light cone on each side
        if (columnMode) side = 0;

        // init state in plc
        for (int i=center-side; i<=center+side; ++i) plc.push_back(state[i]);

        // other states in plc
        int t;
        for (t=1; t<past_depth; ++t) {

            // update state where possible
            for (int i=t; i<wsize-t; ++i) nextState[i] = rule.apply(state[i-1], state[i], state[i+1]);
            swap(state, nextState);

            // update plc
            if (!columnMode) --side;
            for (int i=center-side; i<=center+side; ++i) plc.push_back(state[i]);
        }

        assert(side==0);

        // create flc
        LightCone flc;
        // present is in flc as well
        flc.push_back(center);

        for (; t<past_depth+future_depth-1; ++t) {
            // update state where possible
            for (int i=t; i<wsize-t; ++i) nextState[i] = rule.apply(state[i-1], state[i], state[i+1]);
            swap(state, nextState);

            // update flc
            if (!columnMode) ++side;
            for (int i=center-side; i<=center+side; ++i) flc.push_back(state[i]);
        }

        analyzer.addObservation(plc.getId(), flc.getId());

//         for (int i=0; i<wsize; ++i) cout << state[i];
//         cout << endl;

    }

    analyzer.commitObservations();

    if (!quietMode) {
        cout << "Global complexity: " << analyzer.getGlobalComplexity() << endl;
        cout << "Number of states: " << analyzer.clusters.size() << endl;
        cout << "Number of observed pasts: " << analyzer.observedDistributions.size() << endl;
#ifdef MESINCOM_DEBUG
        for (unsigned int i=0; i<analyzer.clusters.size(); ++i) {
            int nf = analyzer.clusters[i]->distribution.size();
            cout << "nf="<<nf<<"("<<analyzer.clusters[i]->distribution.total<<"): ";
            for (int f=0; f<nf;++f) {
                int flcid = analyzer.clusters[i]->distribution[f].first;
                int cf = analyzer.clusters[i]->distribution[f].second;
                cout <<flcid;
                cout << "("<<cf<<") ";
            }
            cout << endl;
        }
        analyzer.assertNoClusterMatch();
#endif
    }

    else cout << "Rule=" << rule_number << ", complexity=" <<  analyzer.getGlobalComplexity() << endl;

    return 0;
}