function adeq = Adequacy2(resultsfile,model,data,pfailure,nsim)
%model = a structure containing the model you are testing
%data = array containing the data the model is intended to fit
%pfailure = percentage of EPSP failures in the data (assessed from EPSP traces)
%nsim = number of Monte Carlo simulations to use, typically 5000.
% adeq = a structure containing details of the testing

fid=fopen(resultsfile,'a');

ndata=length(data);

% Tests:
onetests = { 'C' 'D' 'chisq20' 'chisq30' 'chisq50' 'chisq75' 'chisq100'};
nonetests = length(onetests);
if isnumeric(pfailure) & ~isnan(pfailure)  & pfailure>=0 & pfailure<= 1
    twotests = {'negloglikelihood'  'skew' 'pfailure' };
else
    twotests = {'negloglikelihood'  'skew' };
end
ntwotests = length(twotests);


disp('Calculating statistics for experimental data...')
% See how well data fits model
exptGOF.negloglikelihood = NegLogLikelihood(data,model);
% Try chisq with several different binsizes:
[tmp,bincentre20] = hist(data,20);
propexpected20 = PropExpected(model,bincentre20); % works out how many are expected in each bin.
exptGOF.chisq20 = ChiSquared(data,propexpected20,bincentre20);

[tmp,bincentre30] = hist(data,30);
propexpected30 = PropExpected(model,bincentre30);
exptGOF.chisq30 = ChiSquared(data,propexpected30,bincentre30);

[tmp,bincentre50] = hist(data,50);
propexpected50 = PropExpected(model,bincentre50);
exptGOF.chisq50 = ChiSquared(data,propexpected50,bincentre50);

[tmp,bincentre75] = hist(data,75);
propexpected75 = PropExpected(model,bincentre75);
exptGOF.chisq75 = ChiSquared(data,propexpected75,bincentre75);

[tmp,bincentre100] = hist(data,100);
propexpected100 = PropExpected(model,bincentre100);
exptGOF.chisq100 = ChiSquared(data,propexpected100,bincentre100);

exptGOF.C = SumSqDiffCDF(data,model);
exptGOF.D = KolmogorovSmirnov(data,model);
exptGOF.skew = skewness(data);
exptGOF.pfailure = pfailure;

% Generate MC data and see how well it fits the model it was generated from
simGOF.negloglikelihood=zeros(1,nsim);
simGOF.chisq20=zeros(1,nsim);
simGOF.chisq30=zeros(1,nsim);
simGOF.chisq50=zeros(1,nsim);
simGOF.chisq75=zeros(1,nsim);
simGOF.chisq100=zeros(1,nsim);
simGOF.C=zeros(1,nsim);
simGOF.D=zeros(1,nsim);

disp('Generating Monte Carlo data ...');
tmp=0;
for sim = 1:nsim
	tmp=tmp+1;
	[simdata{sim},simGOF.pfailure(sim)] = Simulate(model,ndata);
	% Do some goodness of fit tests
	simGOF.skew(sim) = skewness(simdata{sim});   
	simGOF.negloglikelihood(sim) = NegLogLikelihood(simdata{sim},model);
	simGOF.chisq20(sim) = ChiSquared(simdata{sim},propexpected20,bincentre20);
	simGOF.chisq30(sim)  = ChiSquared(simdata{sim},propexpected30,bincentre30);
	simGOF.chisq50(sim)  = ChiSquared(simdata{sim},propexpected50,bincentre50);
	simGOF.chisq75(sim)  = ChiSquared(simdata{sim},propexpected75,bincentre75);
	simGOF.chisq100(sim) = ChiSquared(simdata{sim},propexpected100,bincentre100);
	simGOF.C(sim) = SumSqDiffCDF(simdata{sim},model);
	simGOF.D(sim) = KolmogorovSmirnov(simdata{sim},model);
	disp(['Calculating .... ' num2str(sim/nsim*100) '% '])
end


% Do the One-tailed tests
fig1=figure('position',[1 90 1025 588],'name','one-tailed adequacy tests','paperpositionmode','auto');
fprintf(fid,'%s\n','--------- one-tailed tests');
for test=1:nonetests
	simulatedresults = getfield(simGOF,onetests{test});
	exptldataresults = getfield(exptGOF,onetests{test});
	
	% What percentage of simulated results have a worse GOF stat than the experimental?
	signif(test) = 100*sum(simulatedresults >= exptldataresults) / length(simulatedresults);
	quantile95 = prctile( simulatedresults , 95 );
    if signif(test)>5 % assuming 5% signif level - change this if desired
        outcome = 'PASS';
    else
        outcome = 'FAIL';
    end
	fprintf(fid,'%s: %10.2f %s %s %s %s %f %s %f %s %f %s\n',...
        outcome,signif(test),' % of the simulated data-sets had ',onetests{test},...
        ' worse than or equal to the fitted  model',' (fit had ',exptldataresults,...
        '; mean Monte Carlo was ',mean(simulatedresults),' and 5% quantile was',quantile95,' )');
	
    % Write results to output structure
    adeq.onetests(test).name = onetests{test};
    adeq.onetests(test).signif = signif(test);
    adeq.onetests(test).outcome = outcome;
    
	txt=strrep(onetests{test},'chisq','');
	if ~(strcmp(txt,onetests{test}))
		nbins = str2num(txt);
		fprintf(fid,'%10.2f %s\n',100-100*chi2cdf(exptldataresults,nbins-1),' is the prediction from the chi-squared formula');
	end
	fprintf(fid,'\n');
	% Plots
	if signif(test)<5
		col='r';
	else
		col='b';
	end
	nx=ceil(nonetests/4); ny = 4;
	subplot(nx,ny,test); 
	plot(sort(simulatedresults),[1:nsim]/nsim,col);
	hold on
	plot(exptldataresults*[1 1],[0 1],col);
	title([onetests{test} ' test'])
	ylabel('CDF')
	xlabel(onetests{test})
	plot(quantile95*[1 1],[0 1],'k:');
	plot(get(gca,'xlim'),[1 1]*95/100,'k:');  
    drawnow
end
xlabel(['vertical lines: solid = data, dotted = 95% quantile'])
print(fig1,'-djpeg',[resultsfile '.adeqgraph1.jpg'])

% Do the two-tailed tests.
fprintf(fid,'%s\n','--------- two-tailed tests');
fig2=figure('position',[1 90 1025 300],'name','two-tailed adequacy tests','paperpositionmode','auto');

for test=1:ntwotests
	simulatedresults = getfield(simGOF,twotests{test});
	exptldataresults = getfield(exptGOF,twotests{test});
	
	% What percentage of simulated results have a higher/lower stat than the experimental?
	prophieroreq = 100*sum(simulatedresults >= exptldataresults) / length(simulatedresults);
	proplower = 100*sum(simulatedresults < exptldataresults) / length(simulatedresults);
	loquantile = prctile( simulatedresults , 2.5 );
	hiquantile = prctile( simulatedresults , 97.5 );
	
	signif(test+nonetests) = 2 * min([prophieroreq proplower]);
	if signif(test+nonetests)>5 % assuming 5% signif level - change this if desired
        outcome = 'PASS';
    else
        outcome = 'FAIL';
    end

    fprintf(fid,'%s:\n',outcome)
	fprintf(fid,'%10.2f %s %s %s %f %s %f %s %f %s\n',prophieroreq,' % of the simulated data-sets had ',twotests{test},' higher than or equal to the fitted  model (fit had ',exptldataresults,'; mean Monte Carlo was ',mean(simulatedresults),' and 2.5% quantile was',loquantile,' )');
	fprintf(fid,'%10.2f %s %s %s %f %s %f %s %f %s\n',proplower,' % of the simulated data-sets had ',twotests{test},' lower than the fitted  model (fit had ',exptldataresults,'; mean Monte Carlo was ',mean(simulatedresults),' and 97.5% quantile was',hiquantile,' )');
    fprintf(fid,'%10.5f %s\n\n',signif(test+nonetests),' = significance.');
	
    adeq.twotests(test).name = twotests{test};
    adeq.twotests(test).signif = signif(test);
    adeq.twotests(test).outcome = outcome;

    
    % Plots
	nx=ceil(ntwotests/4); ny = 4;
	if signif(test+nonetests)<5
		col='r';
	else
		col='b';
	end
	subplot(nx,ny,test)
	plot(sort(simulatedresults),[1:nsim]/nsim,col);
	hold on
	plot(exptldataresults*[1 1],[0 1],col);
	title([twotests{test} ' test'])
	ylabel('CDF')
	plot(loquantile*[1 1],[0 1],'k:');
	plot(hiquantile*[1 1],[0 1],'k:');
	plot(get(gca,'xlim'),[1 1]*97.5/100,'k:');
	plot(get(gca,'xlim'),[1 1]*2.5/100,'k:');
	
end
xlabel(['vertical lines: solid = data, dotted = 2.5 and 97.5% quantiles'])

print(fig2,'-djpeg',[resultsfile '.adeqgraph2.jpg'])
fclose('all')

% Return summary:
adeq.worstfailsignif = min(signif);
