function results = testsofuniformity(pdens)
% conducs the Kolmogorov-Smirnov and Cramer-Von-Mises specification
% tests; critical values approximated as in Kroese, Taimre and Botev (2011, Chapter 8, page 339-340)
% INPUTS: pdens - PITS
% OUTPUT: results, which includes results.kspValue  & results.ksStat for
% the KS test p-value and the test statistic; as well as
% results.adpValue & results.adStat for the AD test p-value and test
% statistic
for i = 1:size(pdens,1)
    pit = sort(cleanNaN(pdens(i,:)));
    n = size(pit,2);
    j = 1:n;
    
    dn_up = max(abs(pit-j/n));
    dn_down = max(abs(pit-(j-1)/n));
    dn = max(dn_up, dn_down);
    ksStat = sqrt(n)*dn; 
    
    k = -20:1:20;
    ksProb = (-1).^k.*exp(-2*(k.*ksStat).^2);
    kspValue = 1 - sum(ksProb);
    
    ksStatVec(i,:) = ksStat;
    kspValueVec(i,:) = kspValue;
    
    for j = 1:n
        if (pit(:,j) == 0 || pit(:,n+1-j) == 1)
            atemp(j,1) = NaN;
        else
            atemp(j,1) = (2*j-1)*(log(pit(:,j)*(1-pit(:,n+1-j))));
        end
    end
    
     adStat = -(sum(atemp))/n - n;
     cn = 0.01265 + 0.1757/n;
     adProb = a(adStat);

     
     if adProb < cn
         g1 = g(adProb/cn);
         adProb = adProb + 1/10^5*(370/n^3 + 78/n^2 + 6/n)*g1;
     elseif (adStat >= cn) && (adStat < 0.8)
         [g1 g2] = g((adProb-cn)/(0.8-cn));
         adProb = adProb + 1/10^5*(1365/n^2 + 4213/n)*g2;
     elseif (adProb >= 0.8)
         [g1 g2 g3] = g(adProb);
         adProb = adProb + g3/n;
     end
     adpValue = 1 - adProb;
     adStatVec(i,:) = adStat;
     adpValueVec(i,:) = adpValue;
end

results.kspValue = kspValueVec;
results.ksStat = ksStatVec;

results.adpValue = adpValueVec;
results.adStat = adStatVec;
end

function a = a(x)
b = [1.2337141; 2.00012; 0.247105; 0.0649821; 0.0347962; 0.0116720; 0.00168691]';
c = [1.0776; 2.30695; 0.43424; 0.082433; 0.008056; 0.0003146]';
    if (x > 0) && (x < 2)
        a = x.^(-1.2)*exp(-b(1,1)./x)*(b(1,2) + (b(1,3) - (b(1,4) - (b(1,5) - (b(1,6) - b(1,7)*x)*x)*x)*x)*x);
    else
        a = exp(-exp(c(1,1)-(c(1,2)-(c(1,3) - (c(1,4) - (c(1,5) - c(1,6)*x)*x)*x)*x)*x));
    end
end

function [g1 g2 g3] = g(y)
    d = [0.00022633;    6.54034;    14.6538;    14.458;     8.259;      1.91864]';
    m = [130.2137;      745.2337;   1705.091;   1950.646;   1116.360;   255.7844]';


    g1 = sqrt(y)*(1-y)*(49*y-102);
    g2 = -d(1,1) + (d(1,2) - (d(1,3) - (d(1,4) - (d(1,5) - d(1,6)*y)*y)*y)*y)*y;  
    g3 = -m(1,1) + (m(1,2) - (m(1,3) - (m(1,4) - (m(1,5) - m(1,6)*y)*y)*y)*y)*y;
end