program collate;

const max_res=1000;

var locus,popn,l,p,numfiles,resolution,index,i,num_pops,num_loci:integer;
    outfile,infile:text;
    fname:string;
    totals:array[0..max_res]of real;
    datum,max,dummy:real;
    maxf,cum:real;
    legal:boolean;
    sloci,spops:set of 0..99;
    response:char;

procedure name(var fn:string;p,l:integer);
{Generates file names for the input}
var pc1,pc2,lc1,lc2:char;
begin
pc1:=chr(ord('0')+(p div 10));
lc1:=chr(ord('0')+(l div 10));
pc2:=chr(ord('0')+(p mod 10));
lc2:=chr(ord('0')+(l mod 10));
fn:=('P_'+pc1+pc2+'L_'+lc1+lc2+'.out')
end;



begin
writeln('This program caculates likelihood curves from the output');
writeln('of the program Fst_like');
writeln('The curves are standardised to have area 1, so they can be');
writeln('Treated as probability densities (assuming a uniform prior)');
writeln;
writeln('enter number of populations and number of loci to be combined');
writeln('eg 20 4');
writeln('means 20 population for 4 loci.');writeln;
writeln('Values should be in the range (1-99)');
writeln;
repeat readln(num_pops,num_loci);
       legal:=((num_pops in [1..99])and(num_loci in [1..99]));
       if (not legal) then writeln('Value out of range');
       until legal;

writeln('enter the maximum value of Fst (range: 0-1), ');
writeln('and the number of points (range:10-1000) on the curve');
writeln('which should which should be the values that were used in');
writeln('running Fst_like, and should be common to those all the files');
writeln;
writeln('eg 0.05 100');
writeln('would mean values of Fst up to 0.05 estimated at 100 points');
repeat readln(maxf, resolution);
       legal:=(((maxf>0)and(maxf<1))and((resolution>10)and(resolution<=1000)));
       if (not legal) then writeln('Error values must be in the range (0,1) and (10..1000)');
       until legal;
writeln;
writeln('You can specify which data you want to include, either by');
writeln('Telling the program which to include or which to leave out');
writeln('Which method you you want to use?');
writeln('Reply L or I for Include or Leave out ');
repeat  read(response);
        legal:=(response in ['L','l','I','i']);
        if (not legal) then writeln('Response must be L or I');
        until legal;
If (response in ['L','l']) then
   begin
   spops:=[1..num_pops];sloci:=[1..num_loci];
   repeat
         writeln('At present the following populations will be included');
         writeln('Enter any one you want left out, or 0 if you are satisfied');
         for popn:=1 to num_pops do if (popn in spops) then write(popn:3);
         writeln;
         readln(popn);
         spops:=spops-[popn];
         until popn=0;
   writeln;
   repeat
         writeln('At present the following loci will be included');
         writeln('Enter any one you want left out, or 0 if you are satisfied');
         for locus:=1 to num_loci do if (locus in sloci) then write(locus:3);
         writeln;
         readln(locus);
         sloci:=sloci-[locus];
         until locus=0
   end;

   If (response in ['I','i']) then
   begin
   spops:=[];sloci:=[];
   repeat
         writeln('At present the following populations will be included');
         writeln('Enter any one you want included, or 0 if you are satisfied');
         for popn:=1 to num_pops do if (popn in spops) then write(popn:3);
         writeln('.');
         readln(popn);
         spops:=spops+[popn];
         until popn=0;
   writeln;
   repeat
         writeln('At present the following loci will be included');
         writeln('Enter any one you want included, or 0 if you are satisfied');
         for locus:=1 to num_loci do if (locus in sloci) then write(locus:3);
         writeln;
         readln(locus);
         sloci:=sloci+[locus];
         until locus=0
   end;


for i:=0 to resolution do totals[i]:=0;

for popn:=1 to num_pops do for locus:=1 to num_loci do
  if ((locus in sloci)and(popn in spops)) then
    begin
    name(fname,popn,locus);
    assign(infile,fname);reset(infile);
    for i:=0 to resolution-1 do
        begin
        readln(infile,dummy,datum);
        totals[i]:=totals[i]+datum
        end;
    close(infile)
    end;

max:=-maxint;
for i:=0 to resolution-1 do if totals[i]> max then max:=totals[i];
for i:=0 to resolution-1 do
    if (totals[i]-max)<-80     then totals[i]:=0
                               else totals[i]:=exp(totals[i]-max);
max:=0;
for i:=0 to resolution-1 do max:=max+totals[i];

writeln('enter name of file to score the final curve');
readln(fname);assign(outfile,fname);
rewrite(outfile);
writeln(outfile,'    f        density   cumulative probability');
cum:=0;
for i:=0 to resolution-1 do
    begin
    cum:=cum+totals[i]/max;
    writeln(outfile,i*maxf/resolution:10:5,resolution/maxf*totals[i]/max:10:5,cum:10:5)
    end;
close(outfile)
end.