/*
 * #%L
 * gitools-core
 * %%
 * Copyright (C) 2013 Universitat Pompeu Fabra - Biomedical Genomics group
 * %%
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as
 * published by the Free Software Foundation, either version 3 of the 
 * License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public 
 * License along with this program.  If not, see
 * <http://www.gnu.org/licenses/gpl-3.0.html>.
 * #L%
 */
package org.gitools.analysis.stats;

/**
 * Fisher's Exact Test
 * <p/>
 * Fisher's exact test is a statistical significance test used in
 * the analysis of categorical data where sample sizes are small.
 * <p/>
 * The test is used to examine the significance of the association
 * between two variables in a 2 x 2 contingency table.
 * <p/>
 * More info in http://en.wikipedia.org/wiki/Fisher's_exact_test
 */
public class FisherExactTest {

    private int a;
    private int b;
    private int c;
    private int d;

    private double leftPValue;
    private double rightPValue;
    private double twoTailPValue;


    private final HyperGeometric hyper = new HyperGeometric();

    /**
     * Constructor from an array contingency table:
     *
     * @param ctable ctable[0] = a, ctable[1] = b, ctable[2] = c, ctable[3] = d
     */
    public FisherExactTest(int[] ctable) {
        this.a = ctable[0];
        this.b = ctable[1];
        this.c = ctable[2];
        this.d = ctable[3];
    }

    /**
     * Constructor from the contingency table:
     * <p/>
     * | x | !x|
     * ---+---+---+
     * y | a | b |
     * ---+---+---+
     * !y | c | d |
     * ---+---+---+
     *
     * @param a frequency of x ^ y
     * @param b frequency of !x ^ y
     * @param c frequency of x ^ !y
     * @param d frequency of !x ^ !y
     */
    public FisherExactTest(int a, int b, int c, int d) {
        this.a = a;
        this.b = b;
        this.c = c;
        this.d = d;
    }

    /**
     * Test from a given contingency table like:
     * <p/>
     * | x | !x|
     * ---+---+---+
     * y | a | b |
     * ---+---+---+
     * !y | c | d |
     * ---+---+---+
     *
     * @param a frequency of x ^ y
     * @param b frequency of !x ^ y
     * @param c frequency of x ^ !y
     * @param d frequency of !x ^ !y
     * @return probability = C(a+b, a) * C(c+d, c) / C(n, a+c)
     */
    public double testContingencyTable(int a, int b, int c, int d) {
        this.a = a;
        this.b = b;
        this.c = c;
        this.d = d;

        return testContingencyTable();
    }

    /**
     * Run Fisher's Exact Test using internal contingency table.
     *
     * @return probability = C(a+b, a) * C(c+d, c) / C(n, a+c)
     */
    public double testContingencyTable() {
        int n11 = a;
        int n1_ = a + b;
        int n_1 = a + c;
        int n = a + b + c + d;

        return internalTest(n11, n1_, n_1, n);
    }

    public int getA() {
        return a;
    }

    public int getB() {
        return b;
    }

    public int getC() {
        return c;
    }

    public int getD() {
        return d;
    }

    public double getLeftPValue() {
        return leftPValue;
    }

    public double getRightPValue() {
        return rightPValue;
    }

    public double getTwoTailPValue() {
        return twoTailPValue;
    }

    /**
     * Adapted from http://www.langsrud.com/fisher.htm
     * Created by Oyvind Langsrud
     *
     * @param n11 = a
     * @param n1_ = a + b
     * @param n_1 = a + c
     * @param n   = a + b + c + d
     * @return probability = C(a+b, a) * C(c+d, c) / C(n, a+c)
     */
    private double internalTest(int n11, int n1_, int n_1, int n) {
        double sleft, sright, sless, slarg;
        double p, prob;
        int i, j;

        int max = n1_;
        if (n_1 < max) {
            max = n_1;
        }

        int min = n1_ + n_1 - n;
        if (min < 0) {
            min = 0;
        }

        if (min == max) {
            calcPValues(1.0, 1.0, 1.0, 1.0);
            return 1.0;
        }

        prob = hyper.hyper0(n11, n1_, n_1, n);

        sleft = 0.0;
        p = hyper.hyper(min);
        for (i = min + 1; p < 0.99999999 * prob; i++) {
            sleft += p;
            p = hyper.hyper(i);
        }
        i--;
        if (p < 1.00000001 * prob) {
            sleft += p;
        } else {
            i--;
        }

        sright = 0.0;
        p = hyper.hyper(max);
        for (j = max - 1; p < 0.99999999 * prob; j--) {
            sright += p;
            p = hyper.hyper(j);
        }
        j++;
        if (p < 1.00000001 * prob) {
            sright += p;
        } else {
            j++;
        }

        if (Math.abs(i - n11) < Math.abs(j - n11)) {
            sless = sleft;
            slarg = 1 - sleft + prob;
        } else {
            sless = 1 - sright + prob;
            slarg = sright;
        }

        calcPValues(sleft, sright, sless, slarg);

        return prob;
    }

    private void calcPValues(double sleft, double sright, double sless, double slarg) {

        leftPValue = sless < 1.0 ? sless : 1.0;
        rightPValue = slarg < 1.0 ? slarg : 1.0;
        twoTailPValue = sleft + sright;
        if (twoTailPValue > 1.0) {
            twoTailPValue = 1.0;
        }
    }
}
