/**
 * Created with IntelliJ IDEA.
 * User: rebeccamancy
 * Date: 16/12/2013
 * Time: 14:25
 * To change this template use File | Settings | File Templates.
 */
public class SimCluster {

    // Attributes are an array of simulators
    private Simulator[] clust;     // one per simulator run
    private Ecology eco;

    private Sequences[] seqs;      // one Sequences object per patch (array of these of length numPatches)
    private int numDraws, timeWindowLength; // for convergence ...

    private double[] extinctionTimes;
    private int[][] extPaths;

    /**
     * Constructor that is used for SimCluster used to generate QSD until convergence (thus needs Sequences)
     * @param numSims
     * @param numSeqsPerSimulator
     * @param eco
     */
    public SimCluster(int numSims, int numSeqsPerSimulator, Ecology eco, int numDraws, int timeWindowLength, int maxIndexed){
        // Set up array of simulators
        clust = new Simulator[numSims];
        this.eco = eco;
        for (int i=0; i<numSims; i++) {
            clust[i] = new Simulator(eco, numDraws, timeWindowLength, maxIndexed);
        }

        // Sequences - start with zeros everywhere
        double[] prevSum = new double[numSims];  // One sum per simulator
        double[] prevSumSq = new double[numSims];
        seqs = new Sequences[numSeqsPerSimulator]; // i.e. one seqs object per patch
        for (int i=0; i<numSeqsPerSimulator; i++) {
            seqs[i] = new Sequences(0, prevSum, prevSumSq);
        }
        // Set numDraws and timeWindowLength
        this.numDraws = numDraws;
        this.timeWindowLength = timeWindowLength;
    }


    /**
     * Constructor used when we have a QSD to draw from for initial states for the simulation cluster
     * @param numSims  number of stochastic repetitions
     * @param eco      ecology (combination of habitat and model, with parameter values of e,c,alpha)
     * @param numDraws reporting points in the case where we're meant to be storing these for a trace
     * @param maxTime  time at which to bail out simulating
     * @param qsd      qsd from which to draw starting state
     */
    public SimCluster(int numSims, int numPaths, Ecology eco, int numDraws, int maxTime, QSD qsd){
        // Set up array of simulators drawing initial state from QSD
        clust = new Simulator[numSims];
        this.eco = eco;
        for (int i=0; i<numSims; i++) {
            clust[i] = new Simulator(eco, numDraws, maxTime, qsd.drawFromQSD());
        }
        // Initialise array storing extinction times
        extinctionTimes = new double[numSims];
        extPaths = new int[numPaths][numDraws];
    }

    /**
     * Sets QSD in relevant Ecology by running cluster of simulators until the QSD using the rhat conv test and deOliveira
     * @param maxIndexed         maximum number of states we're allowed to store
     * @param rhatThreshold
     * @param trace              whether to output to the console
     */
    public void run2QSD(int maxIndexed, double rhatThreshold, boolean trace) {

        int numSims = clust.length;
        if (numSims < 2) {
            System.err.println("Can't run to qsd on only one sequence as can't test for convergence! Going to exit now. See properties file.");
            System.exit(1);
        }

        int numPatches = eco.getHab().getNumPatches();
        double[] rhat = new double[numPatches];    // assume zero initialisation far enough away from one

        // Storage for patch status for each patch at each of the drawpoints
        int[][] curSlice = new int[numSims][this.numDraws]; // holds state at each drawpoint for one patch

        // Simulate until convergence of chains - i.e. until the rhat value for each *patch* has converged
        while (Util.largestAbsVal(Util.minus(rhat,1)) > rhatThreshold) {
            // Simulate the next time window, collecting states at sample points into an array (rows=sample points; cols=seqs)
            for (Simulator sim:clust) {
                sim.runDeOliveiraDiscrete(this.timeWindowLength, trace, 0.1, 5); // TODO
            }
            // For each patch, collect current excerpts per patch from sim into array, update sequences, and calculate rhat
            for (int patch = 0; patch<numPatches; patch++) {
                for (int s = 0; s < numSims; s++) {
                    curSlice[s] = clust[s].getExcerptPatch(patch);
                    //System.out.println("numSims:" + numSims);
                    //System.out.println("s:" + s + " patch:" + patch);
                    //Util.print(curSlice[s]);
                }
                rhat[patch] = seqs[patch].testConvergence(Util.transposeMatrix(Util.int2dbl(curSlice)));
                seqs[patch].setrHat(rhat[patch]);
            }
            System.out.print("rhat values: "); Util.print(rhat);

        }
    }

    public void run2extinction(boolean trace) {
        int numSims = clust.length;
        for (int simID = 0; simID < numSims; simID++) {
            // Simulate
            clust[simID].run2Extinction(trace);
            // Store extinction time into array
            extinctionTimes[simID] = clust[simID].getExtinctionTime();
            if (simID < extPaths.length) {
                extPaths[simID] = clust[simID].getExtPath();
            }
        }
    }

    /**
     * Returns the first simulator of the cluster
     * @return
     */
    public Simulator getFirstSim() {
        return clust[0];
    }

    public double[] getExtinctionTimes() {
        return extinctionTimes;
    }

    public int[][] getExtPaths() { return extPaths; }

    /*public int[][] getExtPaths(int startRow, int endRow) {
        System.out.println(extPaths.length + " x " + extPaths[0].length);
        return Util.subArray(extPaths, startRow, endRow, 0, extPaths[0].length-1);
    }*/
}