% Function: This is the program that calculates the p-value based on the
% consistency and number of experiments performed. The program will
% automatically calculate all of the p-value cutoffs for every possible
% consistency.

% Inputs:
% BinomCoef - a matrix that lists out all the possible permutations of 
%    significant hits given the number of experimental replicates performed
% MH_p - the p-value desired as a final false positive rate using the data
%    across all the experiments, given a particular consistency (usually a
%    multiple hypothesis corrected value of 0.05)

% Outputs:
% POutputList - a list of the p-value cutoffs for every possible
% consistency level given the number of experiments performed 

function [POutputList] = BinomPValue(BinomCoef, MH_p)

[BinomLength,~] =  size(BinomCoef);

ProbFunc = cell(BinomLength,1);
N = max(BinomCoef(:));
syms p


for i = 1:BinomLength
    
    ProbFunc{i} = (p/2)^(BinomCoef(i,1))*...
        (1-p)^(BinomCoef(i,2))*(p/2)^(BinomCoef(i,3))*...
        factorial(N)/(factorial(BinomCoef(i,1))*...
        factorial(BinomCoef(i,2))*factorial(BinomCoef(i,3)));

end


CoefList = cell(1,N);
PValueList = cell(1,N);
POutputList = cell(1,N);

for i = 1:N
    CoefList{i} = cell(1,i);
    PValueList{i} = zeros(1,i);
    POutputList{i} = zeros(1,i);
    for j = 1:i
        F_j = find(BinomCoef(:,1)==N+1-j | BinomCoef(:,3)==N+1-j);
        CoefList{i}{j} = setdiff(F_j,vertcat(CoefList{i}{1:j}));
        PValueList{i}(j) = MH_p/i;
    end
end


for alpha = 1: length(CoefList)
    
    for beta = 1:length(CoefList{alpha})
        
        C_i = CoefList{alpha}{beta};

        for gamma = 1:length(C_i)

             if gamma == 1
                P_Fun_a = ProbFunc{C_i(gamma)};
             else
                P_Fun_a = P_Fun_a + ProbFunc{C_i(gamma)};
             end

        end

        P_FunSolve_a = double(vpasolve(P_Fun_a == PValueList{alpha}(beta),...
            p,[0,Inf]));

        P_Out_a = min(P_FunSolve_a(P_FunSolve_a >= 0));

        POutputList{alpha}(beta) = P_Out_a; 
        
    end
    
end