import java.util.ArrayList;
import java.util.Collections;

/**
 * Created with IntelliJ IDEA.
 * User: rebeccamancy
 * Date: 08/11/2013
 * Time: 14:46
 * To change this template use File | Settings | File Templates.
 */
public class Simulator {
    private Ecology eco;    // Combination of a habitat and a model
    private double[] rates; // Maintains current transition rates
    private int[] patchStatus;          // list of patches with their occupancy
    private int numAlive;
    private int numPatches; // Copy of this for convenience

    // For simulating to QSD
    private int[][] excerpt;       // rows=patches; cols=time points
    private int[] extPath; // for storing number of patches occupied at any one time (for outputing paths)
    private double time;           // current time
    private double maxTime;
    private int numDraws, timeWindowLength;      // num of draws per time window when recording state for rhat calculations       todo - do I need this?
    private double interDrawTime, nextSampleTime;     // we record system state at various intervals for calculating rhat
    private int numEvents = 0;  // for counting number of events needed before convergence
    private QSD estQSD;

    // For simulating until extinction
    private double extinctionTime = -1;


    /**
     * Constructor for the Simulators used to generate QSD, in which we don't give a starting state
     * (assumes time=0 and random initial state)
     * @param eco
     * @param numDraws
     * @param timeWindowLength
     * @param maxIndexed
     */
    public Simulator(Ecology eco, int numDraws, int timeWindowLength, int maxIndexed) {    // simulator for QSD
        this.eco = eco;
        this.numPatches = eco.getHab().getNumPatches();
        this.numDraws = numDraws;
        this.timeWindowLength = 0; // not sure if this is required ... todo
        this.maxTime = 0;     // initially set to timeWindowLength todo
        initPatchStatus();    // sets an initial patch status drawn at random with equal probability
        initRates();
        time = 0;
        this.interDrawTime = (double)timeWindowLength/(double)numDraws; // sz=size or length of inter-draw period
        nextSampleTime = interDrawTime;
        excerpt = new int[numPatches][numDraws];
        estQSD = new QSD();  // contains arrays of times, states, numAlive
    }

    /**
     * Constructor with starting state (initial Patch Status)
     * @param eco
     * @param numDraws
     * @param maxTime
     * @param initPatchStatus
     */
    public Simulator(Ecology eco, int numDraws, int maxTime, int[] initPatchStatus) {    // Simulator with starting state
        this.eco = eco;
        this.numPatches = eco.getHab().getNumPatches();
        this.numDraws = numDraws;
        this.maxTime = maxTime;
        this.patchStatus = initPatchStatus;
        initRates();
        time = 0;
        this.interDrawTime = (double)this.maxTime/(double)numDraws; // sz=size or length of inter-draw period
        nextSampleTime = 0;
        //excerpt = new int[numPatches][numDraws];
        extPath = new int[numDraws];
        numAlive = Util.sumArray(patchStatus);
    }

    /**
     * Sets up initial patchStatus by choosing a random a non-extinct initial state
     */
    private void initPatchStatus() {
        int seeds = Util.randInt(1,numPatches);        // number of seeds
        numAlive = seeds;

        ArrayList<Integer> listStatus = new ArrayList<Integer>(numPatches);
        // Generate a new array with the right number of entries in the seed state
        for (int i = 0; i < numPatches; i++) {
            if (i < seeds) {
                listStatus.add(eco.getMod().getSeedState());        // set into seed state
            } else {
                listStatus.add(0);
            }
        }
        // Randomise which patches in which state
        Collections.shuffle(listStatus);
        // Convert back to a primitive integer array
        this.patchStatus = Util.convertIntegers(listStatus);
    }

    /**
     * Initialises rates according to the model (using method in Ecology)
     */
    private void initRates() {
        double[] r = new double[numPatches];
        for (int focalPatch = 0; focalPatch < numPatches; focalPatch++) {
            r[focalPatch] = eco.calcRate(focalPatch, patchStatus);
        }
        rates = r;
    }

    public void run(int maxTime, Boolean trace) {
        double tau;
        int focalPatch, nextState;

        // Loop until maxTime
        while (time < maxTime && numAlive > 0) {
            // Get tau (using rates array) and update time
            tau = genTau();

            // Select patch according to rates, update state and update all rates appropriately
            focalPatch = Util.drawPosit(rates);
            if (trace) {
                System.out.println("Updating site " + focalPatch + " from " + patchStatus[focalPatch] + " to " + eco.getMod().getNextState(patchStatus[focalPatch]));
            }
            // Suggested next state ...
            nextState = eco.getMod().getNextState(patchStatus[focalPatch]);

            // Update number alive if going from 0 to something or something to 0
            if (patchStatus[focalPatch] > 0 && nextState == 0) {
                numAlive --;
            } else if (patchStatus[focalPatch] == 0 && nextState > 0) {
                numAlive ++;
            }
            time = time + tau;
            patchStatus[focalPatch] = nextState;
            //rates = eco.getMod().updateRates(focalPatch, patchStatus, rates, eco.getHab());
            initRates();
            if (trace) {
                System.out.println("Have done the update. numAlive:" + numAlive);
            }

            // If extinct, set extinction time
            if (numAlive==0) {
                this.extinctionTime = time;
            }
        }
        System.out.println("[Simulator.run] Extinction time: " + this.extinctionTime);
    }

    /**
     * Runs the simulator following the deOliveira algorithm until estQSD space limits reached or end of window
     * @param windowLength   how long to run before censoring for rhat
     * @param trace          output what's going on to the console
     * @param tickrate
     * @param numRandJumps
     */
    public void runDeOliveiraDiscrete(int windowLength, Boolean trace, double tickrate, int numRandJumps) {

        double tau;
        int focalPatch, nextState;  // next state for one patch
        double deltat, tickpoint = 0;    // deltaTime = time since last saved to QSD
        int recordIndex = 0;
        maxTime = windowLength + this.maxTime;
        excerpt = new int[numPatches][numDraws];
        int numExts = 0, numticks = 0;     // used to send to a random state the first few times

        // Loop until maxTime (until end of time window or stored enough states)
        while (time < maxTime ) {
            // Get tau (using rates array) and update time
            tau = genTau();

            // Select patch according to rates, update state and update all rates appropriately
            focalPatch = Util.drawPosit(rates);
            if (trace) {
                System.out.println("runDeOliveira: Updating site " + focalPatch + " from " + patchStatus[focalPatch] + " to " + eco.getMod().getNextState(patchStatus[focalPatch]));
            }
            // Suggested next state ...
            nextState = eco.getMod().getNextState(patchStatus[focalPatch]);

            // Store information in excerpt as appropriate for calculating rhat values
            // If next event will take us past a sampling time threshold, save current state to windowrecord
            if (time > nextSampleTime && recordIndex < numDraws) {
                // might have gone past several drawpoints so loop until caught up until end of time window
                while (nextSampleTime < time && recordIndex < numDraws) {
                    for (int patch = 0; patch<numPatches; patch++) {
                        excerpt[patch][recordIndex] = patchStatus[patch];
                    }
                    recordIndex ++; // move to next index in array
                    nextSampleTime = nextSampleTime + interDrawTime;
                }
            }

            // Store information on time in previous state to current estimated QSD
            //?? estQSD.addState(this.numAlive, Util.arr2biStr(patchStatus), tau);
            deltat = (time + tau) - tickpoint;
            numticks = (int) Math.floor(deltat / tickrate); // number of ticks crossed since stored
            // System.out.println(" tickrate:" + tickrate + " time: " + time + " tau:" + tau + " tickpoint:" + tickpoint + "numticks:" + numticks);
            estQSD.add(Util.bin2BitSet(patchStatus), numticks);
            tickpoint = tickpoint + (numticks*tickrate); // move tickpoint time along by number crossed

            // If the last patch is due to die
            if (numAlive == 1 && nextState == 0) {
                if (numExts < numRandJumps) {
                    initPatchStatus();
                    numAlive = Util.sumArray(patchStatus); // todo - this won't work for SEIR or SIRS
                    // Calculate new rates
                    initRates();
                } else {
                    // Pick a new next state from the estimated QSD
                    this.patchStatus = estQSD.drawFromQSD();
                    // Calculate numAlive ...
                    numAlive = Util.sumArray(patchStatus); // todo - this won't work for SEIR or SIRS
                    // Calculate new rates
                    initRates();
                }
            } else { // do the original update
                // Update number alive if going from 0 to something or something to 0
                if (patchStatus[focalPatch] > 0 && nextState == 0) {
                    numAlive --;
                } else if (patchStatus[focalPatch] == 0 && nextState > 0) {
                    numAlive ++;
                }
                // Update patch status for focal patch
                patchStatus[focalPatch] = nextState;
                // Update rates
                initRates();
            }
            // Update time, number of events and rates
            time = time + tau;
            numEvents = numEvents + 1;
            if (trace) {
                System.out.println("Have done the update. numAlive:" + numAlive);
                Util.print(getPatchStatus());
                Util.print(getRates());
            }
        }
    }


    /**
     * Runs the simulator following the deOliveira algorithm until estQSD space limits reached or end of window
     * @param windowLength   how long to run before censoring for rhat
     * @param trace          output what's going on to the console
     * @param maxIndexed     max number of states to index
     *//*
    public void runDeOliveira(int windowLength, Boolean trace, int maxIndexed) {

        double tau;
        int focalPatch, nextState;  // next state for one patch
        int recordIndex = 0;
        maxTime = windowLength + this.maxTime;
        excerpt = new int[numPatches][numDraws];

        // Loop until maxTime (until end of time window or stored enough states)
        //?? while (time < maxTime && estQSD.getArrPointer() < maxIndexed ) {
        while (time < maxTime ) {
            // Get tau (using rates array) and update time
            tau = genTau();

            // Select patch according to rates, update state and update all rates appropriately
            focalPatch = Util.drawPosit(rates);
            if (trace) {
                System.out.println("runDeOliveira: Updating site " + focalPatch + " from " + patchStatus[focalPatch] + " to " + eco.getMod().getNextState(patchStatus[focalPatch]));
            }
            // Suggested next state ...
            nextState = eco.getMod().getNextState(patchStatus[focalPatch]);

            // Store information in excerpt as appropriate for calculating rhat values
            // If next event will take us past a sampling time threshold, save current state to windowrecord
            if (time > nextSampleTime && recordIndex < numDraws) {
                // might have gone past several drawpoints so loop until caught up until end of time window
                while (nextSampleTime < time && recordIndex < numDraws) {
                    for (int patch = 0; patch<numPatches; patch++) {
                        excerpt[patch][recordIndex] = patchStatus[patch];
                    }
                    recordIndex ++; // move to next index in array
                    nextSampleTime = nextSampleTime + interDrawTime;
                }
            }

            // Store information on time in previous state to current estimated QSD
            //?? estQSD.addState(this.numAlive, Util.arr2biStr(patchStatus), tau);
            estQSD.add(Util.bin2BitSet(patchStatus), tau);

            // If the last patch is due to die
            if (numAlive == 1 && nextState == 0) {
                // Pick a new next state from the estimated QSD
                this.patchStatus = estQSD.drawFromQSD();
                // Calculate numAlive ...
                numAlive = Util.sumArray(patchStatus); // todo - this won't work for SEIR or SIRS
                // Calculate new rates
                initRates();
            } else { // do the original update
                // Update number alive if going from 0 to something or something to 0
                if (patchStatus[focalPatch] > 0 && nextState == 0) {
                    numAlive --;
                } else if (patchStatus[focalPatch] == 0 && nextState > 0) {
                    numAlive ++;
                }
                // Update patch status for focal patch
                patchStatus[focalPatch] = nextState;
                // Update rates
                initRates();
            }
            // Update time, number of events and rates
            time = time + tau;
            numEvents = numEvents + 1;
            if (trace) {
                System.out.println("Have done the update. numAlive:" + numAlive);
                Util.print(getPatchStatus());
                Util.print(getRates());
            }
        }
    }
*/
    /**
     * Runs to extinction and sets attribute extinctionTime from an initial state set externally in constructor
     * @param trace
     */
    public void run2Extinction(Boolean trace) {

        double tau;
        int focalPatch, nextState;  // next state for one patch
        int recordIndex = 0;
        excerpt = new int[numPatches][numDraws];

        // Loop until maxTime (until end of time window or stored enough states)
        while (time < this.maxTime && numAlive > 0) {
            // Get tau (using rates array) and update time
            tau = genTau();

            // Select patch according to rates, update state and update all rates appropriately
            focalPatch = Util.drawPosit(rates);
            if (trace) {
                System.out.println("run2Extinction: Updating site " + focalPatch + " from " + patchStatus[focalPatch] + " to " + eco.getMod().getNextState(patchStatus[focalPatch]));
            }
            // Suggested next state ...
            nextState = eco.getMod().getNextState(patchStatus[focalPatch]);

            // Store information in excerpt
            // If next event will take us past a sampling time threshold, save current state to windowrecord
            if (time >= nextSampleTime && recordIndex < numDraws) {         // storing the initial state also ...
                // might have gone past several drawpoints so loop until caught up until end of time window
                while (nextSampleTime < time && recordIndex < numDraws) {
                    for (int patch = 0; patch<numPatches; patch++) {
                        //excerpt[patch][recordIndex] = patchStatus[patch];
                        extPath[recordIndex] = this.numAlive;
                    }
                    recordIndex ++; // move to next index in array
                    nextSampleTime = nextSampleTime + interDrawTime;
                }
            }

            // If the last patch is due to die
            if (numAlive == 1 && nextState == 0) {
                // Set numAlive to zero ...  which should break outside loop on next round
                numAlive = 0; // todo - this won't work for SEIR or SIRS
                // Store extinction time
                extinctionTime = time + tau;
            } else { // do the original update
                // Update number alive if going from 0 to something or something to 0
                if (patchStatus[focalPatch] > 0 && nextState == 0) {
                    numAlive --;
                } else if (patchStatus[focalPatch] == 0 && nextState > 0) {
                    numAlive ++;
                }
                // Update patch status for focal patch
                patchStatus[focalPatch] = nextState;
                // Update rates
                initRates();
            }
            // Update time, number of events and rates
            time = time + tau;
            numEvents = numEvents + 1;
            if (trace) {
                System.out.println("Have done the update. numAlive:" + numAlive);
                Util.print(getPatchStatus());
                Util.print(getRates());
            }
        }
    }

    /**
     * Generates a time to next event on basis of current total rate (Gillespie algorithm)
     * @return   time to next event, tau
     */
    private double genTau() {
        double lambda =  Util.sumArray(this.rates);
        return -Math.log(1 - Math.random()) / lambda;
    }

    // ####################################
    // ------- GETTERS AND SETTERS --------
    // ####################################

    public double[] getRates() {
        return this.rates;
    }

    public void setRate (int patchID, double value) {
        this.rates[patchID] = value;
    }

    public int[] getPatchStatus() {
        return this.patchStatus;
    }

    public int[][] getExcerpt() {
        return this.excerpt;
    }

    public int[] getExcerptPatch(int patchID) {
        return this.excerpt[patchID];
    }

    public double getTime() {
        return this.time;
    }

    public int getNumEvents() {
        return this.numEvents;
    }

    public void setTime(double time) {
        this.time = time;
    }

    public QSD getEstQSD() {
        return this.estQSD;
    }

    public double getExtinctionTime() { return this.extinctionTime; }

    public int[] getExtPath() { return this.extPath; }
}