program Fst_like;

const
  max_length=100;

type
  sample=array[1..max_length]of integer;
  frequs=array[1..max_length]of real;


var
  eachp,eachl,p,popn,l,locus,numalleles,allele,datum,resolution,
  num_pops,num_loci,index,database:integer;
  counts:sample;
  ps:frequs;
  f,tot,maxf:real;
  fname:string;
  infile,outfile:text;
  legal,continue:boolean;
  response,choice:char;

function like(counts:sample;ps:frequs;f:real):real;
{caculates the likelihood for the multinomial Dirichlet
 using the recursive formulae in Balding & Nichols 1995
 Genetica 96:3-12}

  var index,samplesize,t:integer;
      c:sample;
      l:real;

   begin
   c:=counts;
   samplesize:=0;
   for index:=1 to max_length do samplesize:=samplesize+c[index];
   l:=0;index:=1;
   if samplesize>0 then
      repeat
            while c[index]=0 do index:=index+1;
            c[index]:=c[index]-1;
            l:=l+ln(  (c[index]*f+(1-f)*ps[index])/
                      (1+(samplesize-2)*f)  );
            samplesize:=samplesize-1;
            until samplesize=0;

    like:=l
    end;




procedure name(var fn:string;p,l:integer);
{Generates file names for the output}
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{prog}
writeln('This program will calculate likelihood curves for Fst');
writeln('A curve will be calculated every combination of population');
writeln('and stored in the current directory as P_xxL_yy where x & y');
writeln('are the population and locus numbers');
writeln;
writeln('The information from different loci and populations can then');
writeln('be combined using the companion program called Collate');
writeln;writeln;
writeln('First you must chose between three options.');
writeln('You can use Fst to quantify');
writeln;
writeln('A: the differentiation of a population ');
writeln('   from the average gene frequency in all populations');
writeln;
writeln('B: the differentiation of a population');
writeln('   from the average gene frequency in all populations');
writeln('   apart from the population in question');
writeln;
writeln('C: the differentiation of a population from ');
writeln('   another specific population');
writeln('   which we will call the DATABASE population');
writeln;
writeln('Enter A,B or C to chose your option');
repeat readln(choice);
       legal:=(choice in ['A','a','B','b','C','c']);
       if (not legal) then writeln('Your reply must be one of the letters A,B or C');
       until legal;


writeln('enter name of data file');
readln(fname);
assign(infile,fname);
reset(infile);
readln(infile,num_pops);
readln(infile,num_loci);
writeln;writeln;
writeln('Your input file specified a dataset with');
writeln(num_pops:5,' populations');
writeln(num_loci:5,' loci');
if (choice in ['C','c']) then
   begin
   writeln('Which population is the database population?');
   writeln('Enter its number.');
   repeat readln(database);
       legal:=(database in [1..num_pops]);
       if (not legal) then writeln('Try again. There are ',num_pops:5,' populations');
       until legal
   end;
writeln;writeln;
writeln('Specify the maximum plausible value of Fst');
repeat readln(maxf);
       legal:=((maxf<=1)and(maxf>0));
       if (not legal) then writeln('the value should be between 0 and 1');
       until legal;
writeln;writeln;
writeln('Specify the number of points (10-1000) to be estimated on the likelihood curve');
writeln('eg. 100');
repeat readln(resolution);
       legal:=((resolution<=1000)and(resolution>0));
       if (not legal) then writeln('the value should be between 10 and 1000');
       until legal;

writeln;writeln;
writeln('Do you want to calculate the likelihood curves (y/n)');
repeat readln(response);
       legal:=(response in ['y','Y','n','N']);
       if (not legal) then writeln('reply y or n');
       until legal;
continue:=(response in ['y','Y']);

if continue then for popn:=1 to num_pops do for locus:=1 to num_loci do
          begin
          for index:=1 to max_length do
              begin ps[index]:=0;counts[index]:=0 end;

          name(fname,popn,locus);
          assign(outfile,fname);
          rewrite(outfile);

          reset(infile);
          readln(infile,num_pops);
          readln(infile,num_loci);

          for l:=1 to num_loci do if (l=locus) then
              begin
              readln(infile,numalleles);
              writeln('Presently doing calculations for');
              writeln('Population ',popn:4,' Locus ',locus:4,' with ',numalleles:4,' alleles');
              writeln(numalleles:5);
              if numalleles>max_length then
                 begin
                 writeln('Sorry a maximum of ',max_length:5,' alleles alowed in this version');
                 writeln('Error at locus ',l:5);
                 readln
                 end;
              for p:=1 to num_pops do
                  begin
                  for allele:=1 to numalleles do
                      begin
                      read(infile,datum);write(datum:4);
                      if p=popn then counts[allele]:=counts[allele]+datum;
                      if (choice in ['A','a']) then ps[allele]:=ps[allele]+datum;
                      if ((choice in ['B','b'])and(not (p=popn))) then ps[allele]:=ps[allele]+datum;
                      if ((choice in ['C','c'])and(p=database)) then ps[allele]:=ps[allele]+datum;
                      end;
                  readln(infile);writeln
                  end;
              end
          else begin readln(infile,datum); for p:=1 to num_pops do readln(infile,datum) end;

          tot:=0;
          for index:=1 to numalleles do tot:=tot+ps[index];
          for index:=1 to numalleles do ps[index]:=(ps[index]+1)/(tot+numalleles);
          for index:=numalleles+1 to max_length do counts[index]:=0;
          for index:=0 to resolution-1 do
              begin
              write('.');
              writeln(outfile,index*maxf/resolution:10:5,
                      like(counts,ps,index*maxf/resolution):14:5)
              end;
          writeln;
          close(outfile);
          end;{if continue}


end.{prog}