/*****************************************************************************

Project: pyhlogenic tree model with insetions and/or deletions

File name: indel_std.c

Author: Manami Nishizawa

Contents: The standard functions  for insertion deletion

Date:  May 30, 2000

******************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <conio.h>
#include <string.h>
#include <math.h>
#include "indel_std.h"

/*
	Frees allocated block and returns (void *)NULL
*/

extern void * MemFree(void *ptr)
{
	if(ptr !=NULL)  free(ptr);

   return NULL;
}


/* string copy with memory allocation*/
extern char *StringSave(const char *str)
{
	char *des_str;
   int  length;

   length = strlen(str);
   des_str = (char *)calloc(1,length*sizeof(char));
   strcpy(des_str,str);

   return des_str;
}


static char  BASEs[] = {0,84,67,65,71,85,89,82,77,75,83,87,72,66,86,68,78,63,45};
/*{'\0' T C A G U Y R M K S W H B V D N ? - }*/  /*nucleotides*/

static char  AAs[] = {0,65,82,78,68,67,81,69,71,72,73,76,75,77,70,80,83,84,87,89,86,66,90,88,63,42};
/*{'\0' A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  B  Z  X  - *} */
/*amino acids*/

Uint1 ASCIIToBASEs(char base)
{
	Uint1  A;
   int    i;

  if(base == '\0') {return '\0';}
  if(base>=97){base -= 32;}
  A = 19;

  if(base=='U') A = 1; /*treat as 'T'*/
  for(i=1;BASEs[i]!='-';i++)
  {
   	if(base==BASEs[i]) A = (Uint1)i;
  }
  if(A==19){
   if(base=='-'){return 18;}
   else{	printf("ambiguous residue found,%c\n",base);/*getch();*/return '\0';}
  }

   return A;
}

char BASEsToASCII(Uint1 A)
{
   char   base;

  if(A == 0) {return '\0';}
  if(A>18) {return '*';}
  base = BASEs[A];

   return base;
}

Uint1 ASCIIToAA(char amino, short alphabet_code)
{
	Uint1  A;
   int    i;

  if(amino == '\0') {return 0;}
  A = 26;
  if(alphabet_code==1)
  {
   for(i=1;AAs[i]!='*';i++)
   {
   	if(amino==AAs[i]) A = (Uint1)i;
   }
  }
  if(A==26){printf("ambiguous residue found,%c\n",amino);/*getch();*/ return (Uint1)NULL;}

   return A;
}


char AAToASCII(Uint1 A, short alphabet_code)
{
   char   amino;

  if(A == 0) {return '\0';}
  if(A>25) {return '*';}
  if(alphabet_code==1){
   amino = AAs[A];
  }
   return amino;
}

/*simple function*/
void error (char * message)
{ printf("\nerror: %s.\n", message); exit(-1); }
int xtoy (double x[], double y[], int n)
{ int i; for (i=0; i<n; y[i]=x[i],i++) ;  return(0); }


/*round the real number((A-1).5<=x<A.5-->round(x)=A)*/
int	round(double x)
{

   if(((x-floor(x))<0.5)&&((ceil(x)-x)>0.5)){return (int)floor(x);}
   if(((x-floor(x))>=0.5)&&((ceil(x)-x)<=0.5)){return (int)ceil(x);}
   else{
      if(((x-floor(x))==0.0)&&((ceil(x)-x)==0.0)){return (int)x;}
   	printf("?????x=%lf, x-floor(x)=%lf, ceil(x)-x=%lf int(x)=%d\n",x,(x-floor(x)),(ceil(x)-x),(int)x); getch();
   	return  (int)x;
   }

}


/*random number by Wichmann BA & Hill ID.  1982*/
static int z_rndu=137;
static unsigned w_rndu=13757;

void SetSeed (int seed)
{
   z_rndu = 170*(seed%178) + 137;
   w_rndu=seed;
}


double rndu (void)
{
/* U(0,1): AS 183: Appl. Stat. 31:188-190 
   Wichmann BA & Hill ID.  1982.  An efficient and portable
   pseudo-random number generator.  Appl. Stat. 31:188-190

   x, y, z are any numbers in the range 1-30000.  Integer operation up
   to 30323 required.
*/
   static int x_rndu=11, y_rndu=23;
   double r;

   x_rndu = 171*(x_rndu%177) -  2*(x_rndu/177);
   y_rndu = 172*(y_rndu%176) - 35*(y_rndu/176);
   z_rndu = 170*(z_rndu%178) - 63*(z_rndu/178);
   if (x_rndu<0) x_rndu+=30269;
   if (y_rndu<0) y_rndu+=30307;
   if (z_rndu<0) z_rndu+=30323;
   r = x_rndu/30269.0 + y_rndu/30307.0 + z_rndu/30323.0;
   return (r-(int)r);
}

double  RandomReal(double low,double high)
{
	double d;

   d = rndu();

return(low+d*(high-low));
}


int  RandomInteger(int low, int high)

{
   int    k;
	double d;

   d = rndu();
   k = (int)(d*(high-low+1));

return(low+k);
}



double rndnorm (void)
{
/* standard normal variate, using the -Muller method (1958), improved by 
   Marsaglia and Bray (1964).  The method generates a pair of random
   variates, and only one used.
   See N. L. Johnson et al. (1994), Continuous univariate distributions, 
   vol 1. p.153.
*/
   double u1,u2, sumsquare=0;

   for (; ;) {
      u1=2*rndu()-1;
      u2=2*rndu()-1;
      sumsquare=u1*u1+u2*u2;
      if (sumsquare>=0 && sumsquare<=1) break;
   }
   return (u1*sqrt(-2*log(sumsquare)/sumsquare));
}

#define rndexp(mean) (-mean*log(rndu()))

int rndpoisson (double m)
{
/* m is the rate parameter of the poisson
   Numerical Recipes in C, 2nd ed. pp. 293-295
*/
   static double sq, alm, g, oldm=-1;
   double em, t, y;

/* search from the origin
   if (m<5) { 
      if (m!=oldm) { oldm=m; g=exp(-m); }
      y=rndu();  sq=alm=g;
      for (em=0; ; ) {
         if (y<sq) break;
         sq+= (alm*=m/ ++em);
      }
   }
*/
   if (m<12) { 
      if (m!=oldm) { oldm=m; g=exp(-m); }
      em=-1; t=1;
      for (; ;) {
         em++; t*=rndu();
         if (t<=g) break;
      }
   }
   else {
     if (m!=oldm) {
        oldm=m;  sq=sqrt(2*m);  alm=log(m);
        g=m*alm-LnGamma(m+1);
     }
     do {
        do {
           y=tan(3.141592654*rndu());
           em=sq*y+m;
        } while (em<0);
        em=floor(em);
        t=0.9*(1+y*y)*exp(em*alm-LnGamma(em+1)-g);
     } while (rndu()>t);
   }
   return ((int) em);
}

double rndgamma1 (double s);
double rndgamma2 (double s);

double rndgamma (double s)
{
/* random standard gamma (Mean=Var=s,  with shape par=s, scale par=1)
      r^(s-1)*exp(-r)
   J. Dagpunar (1988) Principles of random variate generation,
   Clarendon Press, Oxford
   calling rndgamma1() if s<1 or
           rndgamma2() if s>1 or
           exponential if s=1
*/
   double r=0;

   if (s<=0)      puts ("jgl gamma..");
   else if (s<1)  r=rndgamma1 (s);
   else if (s>1)  r=rndgamma2 (s);
   else           r=-log(rndu());
   return (r);
}


double rndgamma1 (double s)
{
/* random standard gamma for s<1
   switching method
*/
   double r, x=0,small=1e-37,w;
   static double a,p,uf,ss=10,d;

   if (s!=ss) {
      a=1-s;
      p=a/(a+s*exp(-a));
      uf=p*pow(small/a,s);
      d=a*log(a);
      ss=s;
   }
   for (;;) {
      r=rndu();
      if (r>p)        x=a-log((1-r)/(1-p)), w=a*log(x)-d;
      else if (r>uf)  x=a*pow(r/p,1/s), w=x;
      else            return (0);
      r=rndu ();
      if (1-r<=w && r>0)
         if (r*(w+1)>=1 || -log(r)<=w)  continue;
      break;
   }
   return (x);
}

double rndgamma2 (double s)
{
/* random standard gamma for s>1
   Best's (1978) t distribution method
*/
   double r,d,f,g,x;
   static double b,h,ss=0;
   if (s!=ss) {
      b=s-1;
      h=sqrt(3*s-0.75);
      ss=s;
   }
   for (;;) {
      r=rndu ();
      g=r-r*r;
      f=(r-0.5)*h/sqrt(g);
      x=b+f;
      if (x <= 0) continue;
      r=rndu();
      d=64*r*r*g*g*g;
      if (d*x < x-2*f*f || log(d) < 2*(b*log(x/b)-f))  break;
   }
   return (x);
}


double LnGamma (double x)

{
/* returns ln(gamma(x)) for alpha>0, accurate to 10 decimal places.  
   Stirling's formula is used for the central polynomial part of the procedure.
   Pike MC & Hill ID (1966) Algorithm 291: Logarithm of the gamma function.
   Communications of the Association for Computing Machinery, 9:684
*/
   double f=0, fneg=0, z;

   if(x<=0) {
      error("lnGamma not implemented for x<0");
      if((int)x-x==0) { puts("lnGamma undefined"); return(-1); }
      for (fneg=1; x<0; x++) fneg/=x;
      if(fneg<0) error("strange!! check lngamma");
      fneg=log(fneg);
   }
   if (x<7) {
      f=1;  z=x-1;
      while (++z<7)  f*=z;
      x=z;   f=-log(f);
   }
   z = 1/(x*x);
   return  fneg+ f + (x-0.5)*log(x) - x + .918938533204673 
          + (((-.000595238095238*z+.000793650793651)*z-.002777777777778)*z
               +.083333333333333)/x;  
}