/**
 * Created with IntelliJ IDEA.
 * User: rebeccamancy
 * Date: 25/10/2013
 * Time: 16:16
 * To change this template use File | Settings | File Templates.
 */
public class Sequences {
    // The sequence class is used to test for convergence using the Rhat statistic, which approaches 1 as n->infty

    // Properties
    private int prevN;
    private double[] prevSum, prevSumSq;
    private double rHat;

    /**
     * Constructor - before rhat is calculated - sets rhat to zero initially.
     * @param prevN
     * @param prevSum
     * @param prevSumSq
     */
    public Sequences (int prevN, double[] prevSum, double[] prevSumSq) {
        // Set up a sequence from the information read in ...
        this.prevN = prevN;
        this.prevSum = prevSum;
        this.prevSumSq = prevSumSq;
        this.rHat = 0; // Some arbitrary value to start
    }

    /**
     * Tests convergence given the existing information on the sequences and an array of new data, returning current rhat
     * @param curSlice
     * @return
     */
    public double testConvergence( double[][] curSlice) {
        // Calculate rHat on a Sequences object (i.e. given an array of doubles)
        int n, curN, m;
        double[] sumij, sumsq;
        double[] curSum, curSumSq;
        m = this.prevSum.length;

        // Calculate n, sum and sumsq for current window
        curN = curSlice.length;
        curSum = Util.sumArray(curSlice);
        curSumSq = Util.sumSqArray(curSlice);

        // Calculate total n, sum and sumsq
        n = this.prevN + curN;
        sumij = Util.add(this.prevSum, curSum);
        sumsq = Util.add(this.prevSumSq, curSumSq);

        // Intermediate calculations
        double[] psidotjbar = Util.div(sumij,n);

        double psidotdotbar = Util.sumArray(psidotjbar)/m;
        double[] psidotdotbararr = new double[m];
        for (int i=0; i<m; i++) {
            psidotdotbararr[i] = psidotdotbar;
        }

        double mm = m, nn=n;     // convert to doubles to avoid maths problems with / not giving remainders
        double[] sjsq = Util.times( ( Util.minus(sumsq, (Util.div((Util.expArr(sumij,2)),nn) ) ) ), (1/(nn-1))); // row vector containing variance within run for each sequence
        double W = (1/mm) * Util.sumArray(sjsq); // single value - mean within-chain variance

        double B = Util.sumArray(Util.expArr(Util.minus(psidotjbar, psidotdotbararr),2)) * (nn/(mm-1));     // single value ...

        double varhat = ( ((nn-1)/nn) * W )    +    ( B/nn );
        double rhat = Math.sqrt(varhat / W);

        // Update sequences object
        this.prevN = n;
        this.prevSum = sumij;
        this.prevSumSq = sumsq;
        this.rHat = rhat;

        // Return rhat for convenience
        return rhat;
    }

    /**
     * Converts the object to a string for output and debugging
     * @return
     */
    public String toString() {
        String mystring = "Seqs object, prevN:" + this.prevN + " prevSum: " + Util.arr2str(this.prevSum) + " prevSumSq: " + Util.arr2str(this.prevSumSq) + " rhat: " + this.rHat;
        return mystring;
    }

    /**
     * Returns the rhat that was calculated previously
     * @return
     */
    public double getrHat() {
        return this.rHat;
    }

    /**
     * Returns the total number of sample points
     * @return
     */
    public int getPrevN() {
        return this.prevN;
    }

    public void setrHat(double rHat) {
        this.rHat = rHat;
    }

}