admutils.c

#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <math.h>  
// #include <openssl/md2.h>
#include "admutils.h" 

#include  "packit.h"   
static int fastdupnum = 10 ;
static double fastdupthresh = 0.75 ;
static double fastdupkill = 0.75 ;
static int snptab = NO ;

extern int numchrom ;

#define MAXSTR 1024
#define MAXFF   128

unsigned int hashit (char *str)  ;

int countcol(char *fname) {
 FILE *fp ;
 int t ;

 openit(fname, &fp, "r") ;
 t = countcolumns(fp) ;
 fclose(fp) ;
 return t ;
}


int countcolumns(FILE *fp)
{  /* count number of text columns separated by whitespace */
  int i=0,c;
  fpos_t ptr;

  if (fgetpos(fp,&ptr)) {
    printf("error counting columns\n");
    return 0;
  }

  while ( (c = getc(fp)) != '\n' ) {
    if (isgraph(c)) {
      i++;
      while (isgraph(c = getc(fp))) {}
      ungetc(c,fp);
    }
  }
  fsetpos(fp,&ptr);
  return i;
}

void sett1(double *tt, double theta, int numstates) 
{
    if (numstates==2)  { 
     tt[0] = 1.0-theta ;
     tt[1] = theta ;
     tt[2] = 0.0 ;
    }
    if (numstates==3)  { 
     tt[0] = (1.0-theta)*(1.0-theta) ;
     tt[1] = 2.0*theta*(1.0-theta) ;
     tt[2] = theta*theta ;
    }
}

void sett1r(double *tt, double theta, int numstates, double risk)
{
  double y ;
  sett1(tt, theta, numstates) ;
  tt[1] *= risk ;
  tt[2] *= risk*risk ;  
  y = asum(tt, numstates) ;
  vst(tt, tt, 1.0/y, numstates) ;
}

void gettln(SNP *cupt, Indiv *indx, 
  double *ptheta, double *plambda, int *pnumstates, int *pignore) 
/* set theta, lambda numstates */
{
   double theta, lambda ; 
   int numstates, chrom, ignore ;

   theta  = indx->theta_mode ;
   lambda = indx->lambda_mode ;
   ignore = indx->ignore ;
   numstates = 3 ;

   chrom = cupt -> chrom ;

   if (chrom == 23) { 
    theta  = indx->Xtheta_mode ;
    lambda = indx->Xlambda_mode ;
    if (indx -> gender == 'M') numstates = 2;
    if (indx -> gender == 'U') ignore = YES ;
   }
   if (ptheta !=NULL) *ptheta = theta ;
   if (plambda != NULL) *plambda = lambda ;
   if (pnumstates != NULL) *pnumstates = numstates ;
   if (pignore != NULL) *pignore = ignore ;
}


void puttln(SNP *cupt, Indiv *indx, 
  double ptheta, double plambda) 
/* put theta, lambda */
{

   int chrom ;


   chrom = cupt -> chrom ;

   if (chrom == 23) { 
    indx->Xtheta_mode = ptheta;
    indx->Xlambda_mode  = plambda ;
    return ;
   }
   indx->theta_mode = ptheta;
   indx->lambda_mode  = plambda ;
   return ;
}


double hwcheck (SNP *cupt, double *cc) 
{

      int i, n, g ;


      vzero(cc, 3) ;

      n = cupt -> ngtypes ;
      for (i=0; i<n; i++)  {
        g = getgtypes(cupt, i) ;
        if (g<0) continue ;
        ++cc[g] ;
     }
     return (hwstat(cc)) ;

}

double hwcheckx (SNP *cupt, Indiv **indm, double *cc) 
// deals with X  
{

      int i, n, g, ignore, numstates ;
      Indiv *indx ;
      double t, l ; 


      vzero(cc, 3) ;

      n = cupt -> ngtypes ;
      for (i=0; i<n; i++)  {
        indx = indm[i] ;
        gettln(cupt, indx, &t, &l, &numstates, &ignore) ;
        if (ignore) continue ; 
        if (numstates != 3) continue ;
        g = getgtypes(cupt, i) ;
        if (g<0) continue ;
        ++cc[g] ;
     }
     return (hwstat(cc)) ;

}
void hap2dip(SNP *cupt) 
// duplicate chromosomes           
{
   int i, n, g ;

   n = cupt -> ngtypes ;
   for (i=0; i<n; i++) {  
    g = getgtypes(cupt, i) ; 
    if (g<0) continue ;
    if (g==2) g = -1 ; 
    if (g==1) g = 2 ; 
    putgtypes(cupt, i, g) ;
   }
}

void flipalleles(SNP *cupt) 
// flip reference, variant counts
{
   int i, n, g ;
   n = cupt -> ngtypes ;
   for (i=0; i<n; i++) {  
    g = getgtypes(cupt, i) ; 
    if (g<0) continue ;
    putgtypes(cupt, i, 2-g) ;
   }
}

void flipalleles_phased(SNP *cupt) 
// flip reference, variant counts
{
   int i, n, g ;
   n = cupt -> ngtypes ;
   for (i=0; i<n; i++) {  
    g = getgtypes(cupt, i) ; 
    if (g<0) continue ;
    if (g>1) fatalx("genotype > 1 in phased file\n") ;
    putgtypes(cupt, i, 1-g) ;
   }
}


void cntit(double *xc, SNP *c1, SNP *c2) 
{
 int n, i, e, f ;

 n = MIN(c1->ngtypes, c2->ngtypes) ;
 vzero(xc, 9) ;
 for (i=0; i<n; i++) {  
  e = getgtypes(c1, i) ;
  f = getgtypes(c2, i) ;
  if (e<0) continue ;
  if (f<0) continue ;
  ++xc[3*e+f] ;
 }
}

/******** UTILITY FUNCTIONS **********/
void fataly(const char *name)
{
  printf("%s",name);
  exit(1);
}

int compare_doubles (const void *a, const void *b)
{
  const double *da = (const double *) a;
  const double *db = (const double *) b;

  return (*da > *db) - (*da < *db);
}

void pcheck (char *name, char x) 
{
 
  if (name != NULL) return ;
  if (x != CNULL) 
   fatalx ("parameter %c compulsory\n",x) ;
  else fatalx("missing argument\n") ;
}

void printm(double **M, int numstates)  
{
     int i,j ;
     printf("M:\n") ;
     for (i=0; i<numstates; i++)  {  
      for (j=0; j<numstates; j++)  {  
       printf("%9.3f ", M[j][i]) ;
      }
      printf("\n") ;
     }
}


int numvalidgtypes(SNP *cupt) 
{
  int nvalid, n, i, k ;
  if (cupt -> isfake) return 0 ;
  n = cupt -> ngtypes ;
  nvalid = 0 ;
  for (i=0; i<n; i++)  {
   k = getgtypes(cupt, i) ;
   if (k >= 0) ++nvalid ;  
  }
  return nvalid ;
}  

int numvalids(Indiv *indx, SNP **snpmarkers, int fc, int lc)  
{ 
  SNP *cupt ; 
  int idnum, numstates, ignore ;  
  int k, nval= 0 ; 
 
  if (fc>lc) return 0 ; 
  if (lc<0) return 0 ; 
  idnum = indx -> idnum ; 
  for (k=fc; k<=lc; ++k) {   
    cupt = snpmarkers[k] ; 
    if (cupt -> isfake) continue ; 
    gettln(cupt, indx, NULL, NULL, &numstates, &ignore) ; 
    if (ignore) continue ; 
    if (cupt -> ngtypes == 0) continue ; 
    if (getgtypes(cupt, idnum) >= 0) ++nval ; 
  } 
  return nval ; 
} 

double malefreq(Indiv **indmarkers, int numindivs) 

/* pop freq of males in sample */
{
   int i ;
   Indiv *indx ; 
   double cmale, cfemale ;

   cmale = 0 ;
   for (i=0; i<numindivs; ++i) { 
    indx = indmarkers[i] ;
    if (indx -> gender == 'M') ++cmale ;
   }

   cmale /= (double) numindivs ;

   return cmale ;
}

int isimatch(int a, int b)
{
   if (a < 0)  return YES ;
   if (b < 0)  return YES ;
   if (a==b)   return YES ;
   return NO;
}

void gethpos(int *fc, int *lc, SNP **snpm, int numsnps,  
 int xchrom, int lo, int hi)  
{ 
   int k, xfc, xlc, pos ; 
   SNP *cupt ; 
 
   xfc = 9999999 ; 
   xlc = -9999999 ; 
   for (k=0; k<numsnps; k++) {  
    cupt = snpm[k] ; 
    if (cupt -> chrom != xchrom) continue ; 
    pos = cupt -> physpos ; 
    if (pos < lo) continue ; 
    if (pos > hi) continue ; 
    xfc = MIN(xfc, k) ; 
    xlc = MAX(xlc, k) ; 
   } 
   *fc = xfc ;  
   *lc = xlc ; 
} 

void makedir (char * dirname) 
// AT wrote original 
// sets up directory.  Will fail hard if directory does not 
// exist and can't be made 
{
 int fdir ; 
   fdir  = open(dirname,O_RDONLY,0);
   if (fdir >= 0) {  
    close (fdir) ;
    return ;
   }
   fdir = mkdir(dirname,0775);
   if (fdir < 0) {  
    perror("makedir") ;
    fatalx("(makedir) directory %s not created\n") ;
   }
   printf("Created a new directory %s\n",dirname);
}


int
indxindex(char **namelist, int len, char *strid)  
// look for string in list
{
     int k ; 
     for (k=0; k< len; k++) {  
      if (strcmp(namelist[k], strid) == 0) return k ;
     }
     return -1 ;
}

int indindex(Indiv **indivmarkers, int numindivs, char *indid) 
/* hash table would be good here */
{
     int k ; 
     for (k=0; k< numindivs; k++) {  
      if (strcmp(indivmarkers[k] -> ID, indid) == 0) return k ;
     }
     return -1 ;
}

void inddupcheck (Indiv ** indivmarkers, int numindivs)
{
  // fail hard if duplicate
  int t, k, n;
  char **ss;

  n = numindivs ;
  if (n==0) fatalx("numdivs zerp!\n") ;

  ZALLOC (ss, n, char *);
  for (k = 0; k < n; k++) {
    ss[k] = strdup (indivmarkers[k]->ID);
  }

  if (numindivs < 100*1000) { 
   t = finddupalt (ss, n);
  }
  else { 
   freesnpindex ();
   t = finddup (ss, n);
  }
  if (t >= 0)
    fatalx ("duplicate sample: %s\n", ss[t]);

  freeup (ss, n);
  free (ss);
}


int snpindex(SNP **snpmarkers, int numsnps, char *snpid) 
{
     int k ; 
     char **ss ;  

     if (snptab==NO) {  
// set up hash table
      snptab = YES ;
      ZALLOC(ss, numsnps, char *) ;
      for (k=0; k< numsnps; k++) {  
       ss[k] = strdup(snpmarkers[k] -> ID) ;
      }
      xloadsearch(ss, numsnps) ;
      freeup(ss, numsnps) ;
      free(ss) ;
     }
     
     k = xfindit(snpid) ;
     return k ;
}
void freesnpindex() 
{
 if (snptab == NO) return ;
 snptab = NO ; 
 xdestroy() ;
}

int ignoresnp(SNP *cupt)  
{
 if (cupt -> ignore) return YES ;
 if (cupt -> isfake) return YES ;
 if (cupt -> ngtypes == 0) return YES ;
 if (cupt -> isrfake) return NO ; 
 return NO ;
}

double entrop(double *a, int n) 
{
  int i ;
  double ysum, t, ans ;

  ans = 0.0 ; 
  ysum = asum(a,n) ;
  for (i=0; i<n ; i++) {
   t = a[i]/ysum ;
   ans += xxlog2(t) ;
  }
  return -ans ;
}
double xxlog2(double t) 
{
  if (t<=0.0) return t ;
  return t * log(t) / log(2.0) ;
}

void 
testnan(double *a, int n) 
{
  int i ;

  for (i=0; i<n; i++) {
   if (!isfinite(a[i])) fatalx("(testnan) fails:  index %d\n",i) ;
  }
}
void getgall(SNP *cupt, int *x, int n) 
{
  int k, t, a ;
  unsigned char b, w ;

  ivclear(x, -1, n) ;
  if (cupt -> gtypes == NULL) {          
   return ;
  }

  if (!packmode) {  
   copyiarr(cupt-> gtypes, x, n) ;
  }

  if (packmode) {
   k = 0 ;
   for (a=0; 4*a<n; ++a) {  
    w = cupt -> pbuff[a] ;
    for (t = 0; t < 4 ; t++)  { 
     b = w >> 2*(3-t) ;
     x[k] = b & 3 ;
     ++k ; 
     if (k>=n) break ;
    }
   }
  }
  for (k=0; k<n; ++k) {
   if (x[k] == 3) x[k] = -1 ;
  }
  return ;
}

int getgtypes(SNP *cupt, int k) 
{
  char *buff ; 
  int  g ;

  if (cupt -> gtypes == NULL) return -1 ;

  if (packmode)  {  
   buff = cupt -> pbuff ;          
   g = rbuff((unsigned char *)buff, k) ;  
   if (g==3) g = -1 ; 
   return g ;
  }

  return cupt -> gtypes[k] ;
}

void putgtypes(SNP *cupt, int k, int val) 
{
  char *buff ; 
  int   g ;

  if (k>= cupt -> ngtypes) fatalx("(putgtypes)\n") ;
  if (packmode)  {  
   buff = cupt -> pbuff ;          
   g = val ; 
   if (g <0) g=3 ;
   wbuff((unsigned char *)buff, k, g) ;  
   return ;
  }
  cupt ->gtypes[k] = val ;
}

int getep(SNP *cupt, int k) 
{
  char *buff ; 
  int   g ;

  if (cupt -> gtypes == NULL) return -1 ;
  if (k>= cupt -> ngtypes) return -1 ; 
   buff = cupt -> ebuff ;          
   g = rbuff((unsigned char *)buff, k) ;  
   if (g==3) g = -1 ; 
   return g ;

}

void putep(SNP *cupt, int k, int val) 
{
  char *buff ; 
  int   g ;

  buff = cupt -> ebuff ;          
  if (buff == NULL) fatalx("(putep) no buffer\n") ;
  if (k>= cupt -> ngtypes) fatalx("(putep)\n") ;
   g = val ; 
   if (g <0) g=3 ;
   wbuff((unsigned char *)buff, k, g) ;  
   return ;

}

unsigned int hasharr(char **xarr, int nxarr)  
// in application ordering  matters so we hash order dependent 
// code modified at suggestion of Chris Chang
{
  
   int i, n ;  
   unsigned int hash, thash ;

   hash = 0 ;  

   for (i=0; i< nxarr; i++)  { 
    thash = hashit(xarr[i]) ;
    hash *= 17 ;
    hash ^= thash ;
   }
   return hash ;
}
unsigned int hashit (char *str) 
{
/* simple and unimpressive hash function NJP */
   int j, len ; 
   unsigned int hash ;

   hash = 0 ;
   len = strlen(str) ;

   for (j=0; j<len ; j++) {  
     hash *= 23 ;
     hash += (unsigned int) (unsigned char) str[j] ;
   }
   return  hash ;
}
void
wbuff(unsigned char *buff, int num, int g) 
// low level routine writes 2 bits to buffer 
// g should be 0 1 2 or 3   (3 = missing) 
{

   int wnum, wplace ;  
   unsigned char mm, msk, ones ;
   static int ncall = 0 ;
 
   if ((g<0) || (g>3)) fatalx("(wbuff) invalid g value\n", g) ;

   ++ncall ;

   msk = 3 << 6 ;  
   mm = g << 6 ;
   ones = 0XFF  ;

   wnum = num/4 ;  
   wplace = num%4 ;

   mm = mm >> (wplace * 2) ;  
   msk = (msk >> (wplace *2)) ^ ones ;  
   buff[wnum] &= msk ;
   buff[wnum] |= mm ;

/**
     printf("zz %d  %d %d %d %02x\n", num, wnum, wplace, g, buff[wnum]) ;
     printf("yyy %d %d\n", g, rbuff(buff,num)) ;
*/
}
int
rbuff(unsigned char *buff, int num) 
{
   int wnum, wplace, rshft ;  
   unsigned char  b ;
   static int ncall = 0 ;
 
// ++ncall ;


   wnum = num >> 2 ;
   wplace = num & 3 ; 

   rshft = (3-wplace) << 1 ;
   b = buff[wnum] >> rshft ;

   b = b & 3 ;
   return b ;
}
// fast dup code
void setfastdupnum(int  num) 
{
 fastdupnum = num ;
}
void setfastdupthresh(double thresh, double kill) 
{
 fastdupthresh = thresh ;
 fastdupkill = kill ;
}

void
killxhets(SNP **snpmarkers, Indiv **indivmarkers, int numsnps, int numindivs)
{
 SNP *cupt ; 
 Indiv *indx ;
 int i, k, g ;

 for (i=0; i<numsnps; i++) {  
  cupt = snpmarkers[i]  ;
  if (cupt -> ignore) continue ;  
  if (cupt -> isfake) continue ; 
  if (cupt -> chrom != 23) continue ;
  for (k=0; k<numindivs; k++) { 
    indx = indivmarkers[k]  ;
    if (indx -> gender != 'M') continue ;
    g = getgtypes(cupt, k) ; 
    if (g != 1) continue ;
    putgtypes(cupt, k, -1) ;
  }
 }
}
  

void
fastdupcheck(SNP **snpmarkers, Indiv **indivmarkers, int numsnps, int numindivs)
{
 SNP *cupt ; 
 Indiv *indx ;
 int *gtypes ;
 int i, j, k, n ;
 int *snphets, *indsnp, tab[15], ww[15], **codeit, *cc, g, *cbuff ; 
 int *buff, val, vv, lbuff, itry, ilo, ihi ;

 ZALLOC(gtypes, numindivs, int) ;
 ZALLOC(cbuff, 2*numindivs, int) ;
 ZALLOC(codeit, numindivs, int *) ;
 ZALLOC(snphets, numsnps, int) ;
 ZALLOC(indsnp, numsnps, int) ;
 for (i=0; i<numsnps; i++) {  
  cupt = snpmarkers[i]  ;
  if (cupt -> ignore) continue ;  
  if (cupt -> isfake) continue ; 
  if (cupt -> chrom > 22) continue ;
  grabgtypes(gtypes, cupt, numindivs) ;
  for (k=0; k<numindivs; k++) { 
    if (gtypes[k] == 1) ++snphets[i] ;
  }
 }
 ivst(snphets, snphets, -1, numsnps) ;
 isortit(snphets, indsnp, numsnps) ;
// make fastdupnum shots at exact match on 15 snps */
 for (itry = 1; itry < fastdupnum; itry++) {
  ilo = 15*itry ;  
  if ((ilo+15)>=numsnps) break ;
  for (i=0; i<15; i++) { 
   j = indsnp[i+ilo] ; 
   tab[i] = j ;
  }
  n = 0 ;
  for (k=0; k<numindivs; ++k) {  
   indx = indivmarkers[k] ; 
   if (indx -> ignore) continue ;
   for (i=0; i<15; i++) { 
    j = tab[i] ;
    g  = getgtypes(snpmarkers[j], k) ;
    if (g<0) break ;
    ww[i] = g ; 
   }
   if (g < 0 ) continue ;
   cc = codeit[n] = cbuff+2*n ;
   cc[0] = kcode(ww, 15, 4) ;
   cc[1] = k ; 
   ++n  ;
  }


  if (n==0) continue ;
  ipsortit(codeit, NULL, n, 2) ;  

  buff = gtypes ; lbuff = 0; val = -1 ;  

  for (i=0; i<n; i++)  {  
   cc=codeit[i] ;  
   vv = cc[0] ;
   if (vv != val) {  
    cdup(snpmarkers, indivmarkers, numsnps, buff, lbuff) ;
    lbuff = 0 ; 
    val = vv ; 
   }
   buff[lbuff] = cc[1] ;
   ++lbuff ;
  }
  cdup(snpmarkers, indivmarkers, numsnps, buff, lbuff) ;
 }  // itry

 free(snphets) ;
 free(indsnp) ;
 free(gtypes) ;
 free(codeit) ;
 free(cbuff) ;
}

void slowdupcheck(SNP **snpm, Indiv **indm, int nsnp, int k1, int k2) 
{ 
 static long ncall = 0 ;
 SNP * cupt ;
 Indiv *inda, *indb ;
 double ytot, yhit ;
 int g1, g2,  match, nomatch ;
 int i1, i2, j ;

 ++ncall ;

 if (ncall<=1) {
  printf("fastdupthresh, kill: %9.3f %9.3f\n", fastdupthresh, fastdupkill) ;
 }
   match = nomatch = 0 ;
   for (j=0; j<nsnp; ++j) {  
    cupt = snpm[j] ;
    if (cupt -> ignore) continue ;
    if (cupt -> isfake) continue ;
    g1 = getgtypes(cupt, k1) ;
    g2 = getgtypes(cupt, k2) ;
    if ( (g1<0) || (g2<0) ) continue ;
    if (g1==g2) ++match ;
    if (g1!=g2) ++nomatch ;
   }

   inda = indm[k1] ;
   indb = indm[k2] ;
   ytot = (double) (match + nomatch) ;
   yhit = ((double) match) / ytot  ; 

   if (yhit>fastdupthresh) { 
    printdup(snpm, nsnp, inda, indb, match, nomatch)  ;
    if (yhit>fastdupkill) killdup(inda, indb, snpm, nsnp) ;
  }

}

void cdup(SNP **snpm, Indiv **indm, int nsnp, int *buff, int lbuff) 
{ 
 static int ncall = 0 ;
 SNP * cupt ;
 Indiv *inda, *indb ;
 double ytot, yhit ;
 int g1, g2, k1, k2, match, nomatch ;
 int i1, i2, j ;
#define MINMARK 200 


 if (lbuff <= 1) return ;
 ++ncall ;

 if (ncall<=1) {
//printf ("cdup: %d\n", ncall) ;
  printf("fastdupthresh, kill: %9.3f %9.3f\n", fastdupthresh, fastdupkill) ;
//printimat(buff, lbuff, 1) ;
 }
// if (ncall==1) printf(" cdup %9.3f %9.3f\n", fastdupthresh, fastdupkill) ;
 for (i1=0; i1<lbuff; ++i1) {  
  for (i2=i1+1; i2<lbuff; ++i2) {
   k1 = buff[i1] ; 
   k2 = buff[i2] ;  
   match = nomatch = 0 ;
   for (j=0; j<nsnp; ++j) {  
    cupt = snpm[j] ;
    if (cupt -> ignore) continue ;
    if (cupt -> isfake) continue ;
    g1 = getgtypes(cupt, k1) ;
    g2 = getgtypes(cupt, k2) ;
    if ( (g1<0) || (g2<0) ) continue ;
    if (g1==g2) ++match ;
    if (g1!=g2) ++nomatch ;
   }

   inda = indm[k1] ;
   indb = indm[k2] ;
   ytot = (double) (match + nomatch) ;
   if (ytot< MINMARK) continue ;
   yhit = ((double) match) / ytot  ; 

   if (yhit>fastdupthresh) { 
    printdup(snpm, nsnp, inda, indb, match, nomatch)  ;
    if (yhit>fastdupkill) killdup(inda, indb, snpm, nsnp) ;
   }
  }
 }
}

void killdup(Indiv *inda, Indiv *indb, SNP **snpm, int nsnp) 
{
       int t1, t2 ;
       Indiv *indx ;

       t1 = numvalids(inda, snpm, 0, nsnp-1) ;
       t2 = numvalids(indb, snpm, 0, nsnp-1) ;
       indx = inda ; 
       if (t1>t2) indx = indb ;
       indx -> ignore = YES ; 
       printf("dup.  %s ignored\n", indx -> ID) ;
}

void printdup(SNP **snpm, int numsnp, Indiv *inda, Indiv *indb, int nmatch, int nnomatch) 
{
       int t1, t2 ;
       double y ;

       if (nmatch<=0) return ;
       if (inda -> ignore) return ;
       if (indb -> ignore) return ;

       t1 = numvalids(inda, snpm, 0, numsnp-1) ;
       t2 = numvalids(indb, snpm, 0, numsnp-1) ;
       printf("dup? %s %s  match: %d mismatch: %d   %d %d ", 
        inda->ID, indb -> ID, nmatch, nnomatch, t1, t2) ;
       printf("%20s ", inda->egroup) ;
       printf("%20s", indb->egroup) ;
       y = nnomatch / (double) (nnomatch+nmatch) ;
       printf(" %9.3f", y) ;
       printnl() ; 
}

int kcode(int *w, int len, int base) 
{
 int i, t ;
 t = 0; 
 for (i=0; i<len; i++) {  
  t *= base ;  
  t += w[i] ;
 }
 return t ;
}

void
 grabgtypes(int *gtypes, SNP *cupt, int numindivs) 
{

 int k ;

 for (k=0; k<numindivs; k++) { 
  gtypes[k] = getgtypes(cupt, k) ;
 }

}

double kurtosis(double *a, int n) 
{

  double  y1, y2, y4 ;  
  double *w  ;

  ZALLOC(w, n, double) ;

  y1 = asum(a,n) / (double) n ;  
  vsp(w, a, -y1, n) ; 
  
  y2 = asum2(w, n) / (double) n ;
  vst(w, w, 1.0/sqrt(y2), n) ;

  vvt(w, w, w, n) ;  

  y4 = asum2(w, n) / (double) n ;


  free(w) ;
  return y4 ; 

}

void printvers(char *progname, char *vers) 
{
 printf("## %s  version: %s", progname, vers) ; 
 printnl() ;
}
int numvalidind(Indiv **indivmarkers, int  numind)  
{ 
  int i ;              
  int nvalids = 0 ;
  Indiv *indx ;

  for (i=0; i< numind;  i++) {  
   indx = indivmarkers[i] ; 
   if (indx -> ignore) continue ;
   ++nvalids ;
  }
  return nvalids ; 
} 

void numvalidgtallind(int *x, SNP **snpm, int numsnps, int numind) { 
// count valids for all 
  int n = 0 ;
  int k, t, j ; 
  SNP *cupt ;
  int *z ;  
  
  ZALLOC(z, numind, int) ;
  ivzero(x, numind) ;
  for (k=0; k< numsnps; ++k) {  
   cupt = snpm[k]  ;
   if (cupt -> ignore) continue ;
   getgall(cupt, z, numind) ;
   for (j=0; j<numind; ++j) { 
    if (z[j] >=0) ++x[j] ;
   }
  }
  free(z) ;
}


int numvalidgtind(SNP **snpm, int numsnps, int ind) { 
  int n = 0 ;
  int k, t ; 
  SNP *cupt ;
  
  for (k=0; k< numsnps; ++k) {  
   cupt = snpm[k]  ;
   if (cupt -> ignore) continue ;
   t = getgtypes(cupt, ind) ;
   if (t >= 0)   ++n ;
  }
  return n ;
}

int numvalidgt(Indiv **indivmarkers, SNP *cupt)  
     // like numvalidgtypes but tests ignore   
{ 
  int n, k, nvalids ; 
  if (cupt -> gtypes == NULL) return 0 ; 
  nvalids = 0 ; 
  n = cupt -> ngtypes ; 
  for (k=0; k<n; k++) {   
    if (indivmarkers[k] -> ignore)  continue ; 
    if (getgtypes(cupt, k) >=0) ++nvalids ; 
  } 
  return nvalids ; 
} 
int numvalidgtx(Indiv **indivmarkers, SNP *cupt, int affst)  
// like numvalidgtypes but tests ignore   counts only when status=affst
{ 
  int n, k, nvalids ; 
  if (cupt -> gtypes == NULL) return 0 ; 
  nvalids = 0 ; 
  n = cupt -> ngtypes ; 
  for (k=0; k<n; k++) {   
    if (indivmarkers[k] -> ignore)  continue ; 
    if (indivmarkers[k] -> affstatus != affst) continue ;
    if (getgtypes(cupt, k) >=0) ++nvalids ; 
  } 
  return nvalids ; 
} 

int isxmale(SNP *cupt, Indiv *indx) 
{
  if (cupt -> chrom != 23) return NO ; 
  if (indx -> gender != 'M') return NO ;
  return YES ;
}

void
printmatz(double *ww, char **eglist, int n) 
{
  int i, j , x ;
  printf(" %4s", "   ") ;
  for (i=0; i<n; i++)  {  
   printf(" %4s", get3(eglist[i])) ; 
  }
  printnl() ;
  for (i=0; i<n; i++)  {  
    printf("%4s", get3(eglist[i])) ;
    for (j=0; j<n; j++) { 
     x = nnint(1000*ww[i*n+j]) ; 
     printf(" %4d", x)  ;              
    }
    printnl() ;
  }
}
void
printmatz5(double *ww, char **eglist, int n) 
{
  int i, j , x ;
  printf(" %5s", "   ") ;
  for (i=0; i<n; i++)  {  
   printf(" %5s", get3(eglist[i])) ; 
  }
  printnl() ;
  for (i=0; i<n; i++)  {  
    printf("%5s", get3(eglist[i])) ;
    for (j=0; j<n; j++) { 
     x = nnint(1000*ww[i*n+j]) ; 
     printf(" %5d", x)  ;              
    }
    printnl() ;
  }
}
void
printmatz10(double *ww, char **eglist, int n) 
{
  int i, j , x ;
  printf(" %5s", "   ") ;
  for (i=0; i<n; i++)  {  
   printf(" %10s", get3(eglist[i])) ; 
  }
  printnl() ;
  for (i=0; i<n; i++)  {  
    printf("%5s", get3(eglist[i])) ;
    for (j=0; j<n; j++) { 
     x = nnint(1000000*ww[i*n+j]) ; 
     printf(" %10d", x)  ;              
    }
    printnl() ;
  }
}



char *get3(char *ss) 
{
  return getshort(ss, 3) ;
}
char *getshort(char *ss, int n) 
{
  static char xxx[MAXSTR] ;
  strcpy(xxx, ss)  ;
  xxx[n-2] = ' ' ;  
  xxx[n-1] = CNULL ;
  return xxx ;
}




int setid2pops(char *idpopstring, Indiv **indmarkers, int numindivs) 
// replace pop by ID for certain samples
{
#define MAXS 1000 
  char *spt[MAXS] ;
  char *sx ;
  int nsplit, k, t ;
  Indiv *indx ;
 
  nsplit = splitupx(idpopstring, spt, MAXS, ':') ;
  for (k=0; k<nsplit; ++k) { 
   sx = spt[k] ;
   t = indindex(indmarkers, numindivs, sx) ; 
   if (t<0) { 
    printf("(setid2pops): %s not found\n", sx) ;
    continue ;
   }
   indx = indmarkers[t] ;
   freestring(&indx -> egroup) ;
   indx -> egroup = strdup(indx -> ID) ;
  }
  freeup(spt, nsplit) ;
  return nsplit ;
}



int gvalm(char cc, char cm, char c1, char c2, int minfilterval) 
{
   int x=0, t ;
   char cb[2] ; 

   if (fvalid(cm, minfilterval) == NO) return 9 ;      
   if (isiub2(cc) == NO) return 9  ;  
   t = iubcbases(cb, cc) ; 

   if (cb[0] == c1) ++x ;
   if (cb[1] == c1) ++x ;


   return x ;

}

int fvalid(char cm, int minfilterval) 
// is cm indicating valid? 
{
  int t ; 

  if (minfilterval < 0) return YES ;
  t = (int) cm - (int) '0' ;

  if (t<0) return NO;
  if (t>=10) return NO ;
  if (t<minfilterval) return NO ;

  return YES ; 
}


int setfalist(char **poplist, int npops, char *dbfile, char **iublist, char *table_path) 
{
	int t;
	for (t = 0; t < npops; ++t) {
		iublist[t] = strdup(table_path);
		iublist[t] = (char*) realloc(iublist[t], 64);
		iublist[t] = strcat(iublist[t], poplist[t]);
		if ((!strcmp (poplist[t], "Chimp") || !strcmp (poplist[t], "Href")) && strcmp (dbfile, ".fa")) {
			free (iublist[t]);
			iublist[t] = "NULL";
		} else 
			iublist[t] = strcat(iublist[t], dbfile);
	}
	return npops;
}
  

int mkmaplist(char *fname, char **kword, char **val) 
// return list of (kword, values) sorted in descending kword lenth
// return # of keywords
{
   FILE *fff ;  
  char line[MAXSTR + 1];
  char ww[MAXSTR];
  char rest[MAXSTR];
  int len, plen, nmap, n ;
  char **tkword, **tval ; 
  int *slenm, *indx, k, x, l ;

   
   
   openit(fname, &fff, "r") ;

   n = numlines(fname) ; 
   ZALLOC(tkword, n, char *) ;
   ZALLOC(tval, n, char *) ;
 
   nmap = 0 ; 

  line[MAXSTR] = '\0';          // defensive programmming
  while (fgets (line, MAXSTR, fff) != NULL) {

    len = strlen (line);

    if (isspace (line[len - 1]))
      line[len - 1] = CNULL ;

    if (first_word (line, ww, rest) > 0) {

      if (ww[0] == '#')
        continue;
      plen = strlen (ww);
      if (ww[plen - 1] != ':') continue ; 
      ww[plen-1] = CNULL ;
      if (upstring(ww) == NO) continue ; 
  
      tkword[nmap] = strdup(ww) ; 
      tval[nmap] = strdup(rest) ; 
      ++nmap ; 

    }

  }

  fclose(fff) ;

  if (nmap == 0) { 
   freeup(tkword, n) ;
   freeup(val, n) ;
   return 0 ;
  }
 
// sort in descending order of length
   ZALLOC(slenm, nmap, int) ; 
   ZALLOC(indx, nmap, int) ; 
   for (k=0; k<nmap; ++k) { 
    x = 10000 * strlen(tkword[k]) + (nmap - k) ; 
    slenm[k] = -x ;
   }
   isortit(slenm, indx, nmap) ; 
   for (k=0; k<nmap; ++k)  { 
    x = indx[k] ;
    kword[k] = strdup(tkword[x]) ;
    val[k] = strdup(tval[x]) ;
   }


  freeup(tkword, n) ;
  freeup(tval, n) ;
  free(slenm) ;
  free(indx) ;

  return nmap ;

}
int getfalist(char **poplist, int npops, char *dbfile, char **iublist)  

{
 char line[MAXSTR+1] ;
 char *spt[MAXFF], *sx ;
 char c ;
 int nsplit ;
 int  tt, t, k, s, nx = 0  ;
 FILE *fff ;
 char *scolon ; 
 char **kword, **val ;
 int nn, nu ; 
 int isempty = -1  ; 
  

  for (k=0; k<npops; ++k) { 
   if (iublist[k] == NULL) isempty = k ;
  }
  if (isempty<0) return 2 ;
  if (dbfile == NULL) {
   fatalx("no database and file %s not set\n", poplist[isempty]) ;
  }

  openit(dbfile, &fff, "r") ;

  nn = numlines(dbfile) ;

  ZALLOC(kword, nn, char *) ;
  ZALLOC(val, nn, char *) ;

  nu =  mkmaplist(dbfile, kword, val)  ; 

  while (fgets(line, MAXSTR, fff) != NULL)  {
   nsplit = splitup(line, spt, MAXFF) ; 
   if (nsplit<1) continue ;
   sx = spt[0] ;
   if (sx[0] == '#') { 
    freeup(spt, nsplit) ;
    continue ;
   }
   t = indxstring(poplist, npops, spt[0]) ; 
   if (t<0) { 
    freeup(spt, nsplit) ; 
    continue ;
   }
   if (iublist[t] != NULL) continue ;
    sx = spt[2] ;
    tt = strcmp(sx, "NULL") ; 
    if (tt == 0) {  
     iublist[t] = NULL ;
    }
    else {
     mapstrings(&sx, kword, val, nu) ;  
     iublist[t] = strdup(sx) ; 
    }
    freeup(spt, nsplit) ;
    ++nx ;
  }

   fclose(fff) ;
   return nx ;

}


int getdbname(char *dbase, char *name, char **pfqname) 
{
 int n, k, t, i ; 
 char **tpoplist, **tfalist ; 
 int numfalist ; 

// fprintf(stderr, "call numlines in [getdbname]\n");

 ZALLOC(tpoplist, 1, char *) ;
 ZALLOC(tfalist, 1, char *) ;
 tpoplist[0] = strdup(name) ; 
 numfalist = getfalist(tpoplist, 1, dbase, tfalist) ;	
 if (numfalist == 0) fatalx("name: %s not found in %s\n", name, dbase) ;

 *pfqname = strdup(tfalist[0]) ;
 
 freeup(tpoplist, 1) ;
 freeup(tfalist, 1) ;
 
 return 1 ; 
}


int cmap(SNP **snpmarkers, int numsnps) 
{
   int t, k ; 
   double y1, y2 ; 
   SNP *cupt ;  
   for (k=1 ; k<=10; ++k) { 
    t = ranmod(numsnps) ;
    cupt = snpmarkers[t] ; 
    y1 = cupt -> genpos ; 
    y2 = cupt -> physpos / 1.0e8  ; 
    if (fabs(y1-y2) > .001) return YES ;
   }
   return NO ;
}


void setoutfiles(char **pind, char **psnp, char **pgeno, char *stem) 
{
  char ss[MAXSTR] ; 

  snprintf(ss, MAXSTR, "%s.ind", stem) ; 
  checkwrite(ss) ; 
  *pind = strdup(ss) ; 

  snprintf(ss, MAXSTR, "%s.snp", stem) ; 
  checkwrite(ss) ; 
  *psnp = strdup(ss) ; 


  snprintf(ss, MAXSTR, "%s.geno", stem) ; 
  checkwrite(ss) ; 
  *pgeno = strdup(ss) ; 

  printf("output files set from %s\n", stem) ; 

}

void setinfiles(char **pind, char **psnp, char **pgeno, char *stem) 
{
  char ss[MAXSTR] ; 

  snprintf(ss, MAXSTR, "%s.ind", stem) ; 
  fcheckr(ss) ; 
  *pind = strdup(ss) ; 

  snprintf(ss, MAXSTR, "%s.snp", stem) ; 
  fcheckr(ss) ; 
  *psnp = strdup(ss) ; 


  snprintf(ss, MAXSTR, "%s.geno", stem) ; 
  fcheckr(ss) ; 
  *pgeno = strdup(ss) ; 

  printf("input files set from %s\n", stem) ; 

}


int loaddiplike (double *dip, unsigned char *sp) 
{
  int rl2 = 4, sval ;  
  unsigned short bb[2] ; 
  
  sval = (1 << 16) - 1 ;
  memcpy(bb, sp, rl2) ;

  vclear(dip, -1, 3) ;
  if ((bb[0] == sval) && (bb[1] == sval)) return -1 ;

  dip[0] = bb[0] ; 
  dip[2] = bb[1] ; 
  dip[1] = sval - (bb[0] + bb[1]) ; 
  if (dip[1] < 0) fatalx("(loaddiplike) bad bb\n") ;
  bal1(dip, 3) ;
  return 1 ;

}


int geteg(char *name, char **list)
// modified from getlist 
{
  FILE *fff ;
  char line[MAXSTR] ;
  char *spt[MAXFF] ;
  char *sx ;
  int nsplit, num=0, t, k ;

  num = 0;
  if (name == NULL) fatalx("(geteg)  no name")  ;
  openit(name, &fff, "r") ;
  while (fgets(line, MAXSTR, fff) != NULL)  {
   nsplit = splitup(line, spt, MAXFF) ;
   if (nsplit==0) continue ;
   sx = spt[0] ;
   if (sx[0] == '#') {
    freeup(spt, nsplit) ;
    continue ;
   }
   for (k=0; k<nsplit ; ++k) { 
    sx = spt[k] ;
    t = indxstring(list, num, sx) ;
    if (t>=0) continue ;
    list[num] = strdup(sx) ;
    ++num ;
   }
   freeup(spt, nsplit) ;
  }
  fclose(fff) ;
  return num ;
}

int setgklist(Indiv **indm, int numindivs, char **eglist, int numeg)  
{

 int k, t ; 
 Indiv *indx ; 

 for (k=0 ; k<numindivs; ++k) { 
  indx = indm[k] ;
  indx -> gkode = t ; 
  t = indxstring(eglist, numeg, indx -> egroup) ;
  indx -> gkode = t ; // correct when t < 0  
 }

 return 1 ; 

}

int setgk(Indiv **indm, int numindivs, char *pname, char ***peglist, int *pnumeg) 
{
 int k, numeg, t ; 

 char **eglist ;
 static char **egl ; 
 Indiv *indx ; 
 int ncol ; 

 for (k=0 ; k<numindivs; ++k) { 
  indm[k] -> gkode = -1 ;
 }
 if (pname == NULL) { 
   if (pnumeg == NULL) return 0 ;
   *pnumeg = 0 ; 
   *peglist = NULL  ; 
   return 0 ; 
  }
 t = numlines(pname) ; 
 ncol = numcols(pname) ; 
 t *= MAX(ncol, 4) ;
 ZALLOC(eglist, t, char *) ; 
 numeg = geteg(pname, eglist) ;
 for (k=0 ; k<numindivs; ++k) { 
  indx = indm[k] ;
  t = indxstring(eglist, numeg, indx -> egroup) ;
  indx -> gkode = t ; // correct when t < 0  
 }

 if (pnumeg != NULL) { 
  *pnumeg = numeg ; 
  ZALLOC(egl, numeg, char *) ; 
  copystrings(eglist, egl, numeg) ; 
  *peglist = egl ;

 }
 
 return 1 ;
}
int gkignore(Indiv **indivmarkers, int numindivs)
{
 int i, n=0 ; 

 Indiv *indx ; 
  for (i = 0; i < numindivs; i++) {
    indx = indivmarkers[i];

    if (indx -> gkode < 0) { 
     indx -> ignore = YES ;
     ++n ;
    }
  }
  return n ;
}

void setpopsize(int *popsize, char **poplist, int numeg, Indiv **indivmarkers, int numindivs) 
{

   int k, d ; 
   char *sx ;
   Indiv *indx ; 

 
   if (indivmarkers == NULL) { 
    ivclear(popsize, -1, numeg) ; 
    return ;
   }
   ivzero(popsize, numeg) ; 

   for (k=0; k<numindivs; ++k) { 
    indx = indivmarkers[k] ;
    if (indx -> ignore) continue ;
    sx = indx -> egroup ;
    d = indxstring(poplist, numeg, sx) ; 
    if (d >= 0) ++popsize[d] ;  
  }

}


void atopchrom (char *ss, int chrom)
{

  // ancestry chromosome -> map convention  

/**
  if ( chrom == numchrom+1 )  {
    strcpy(ss, "X") ;
    return ;
  }
  else if ( chrom == numchrom+2 )  {
    strcpy(ss, "Y") ;
    return ;
  }
*/
  sprintf (ss, "%d", chrom);
}


int ptoachrom (char *ss)
{
  // map -> ancestry  
  char c;
  c = ss[0];

  if (c == 'X')
    return (numchrom + 1);
  if (c == 'Y')
    return (numchrom + 2);
  return atoi (ss);
}

int strcmpreg(char *a, char *b) 
// deal with X 23, Y 24 ambiguity
{
  int sa, sb, t ; 
    
  t = strcmp(a, b) ; 
  if (t==0) return t ; 

  sa = ptoachrom(a) ; 
  sb = ptoachrom(b) ; 
  if (sa==sb) return 0 ; 

  return t ; 

}

void checkwrite(char *fname) 
{
  if (fname == NULL) return ; 
  if (canwrite(fname) == YES) return ;

  fatalx("file %s can not be written.  check file permissions!\n", fname) ;  


}
#define BIGL 40

void testindlen(Indiv **indm, int n) 
{
 int k, maxlen = 0 ; 
 Indiv *indx, *indmax ; 

 if (n==0) return ;
 for (k=0; k<n; ++k) {
  indx = indm[k] ;
  if (strlen(indx -> ID) > maxlen) { 
   maxlen = strlen(indx -> ID) ;
   indmax = indx  ;
  }
 }
 if (maxlen > BIGL) printf("deprecated long ID: %s\n", indmax -> ID) ;
}

void testsnplen(SNP  **snpm, int n) 
{
 int k, maxlen = 0 ; 
 SNP *cupt, *cuptmax ; 

 if (n==0) return ;
 for (k=0; k<n; ++k) {
  cupt = snpm[k] ;
  if (strlen(cupt -> ID) > maxlen) { 
   maxlen = strlen(cupt -> ID) ;
   cuptmax = cupt  ;
  }
 }
 if (maxlen > BIGL) printf("deprecated long ID: %s\n", cuptmax -> ID) ;
}


int getsnpvals(double *vals,  SNP **snpmarkers, int numsnps, char *vname)
{
  char **names, *sx;
  double **xx ;
  int n, t, k ;
  SNP *cupt ;
  int nvalid = 0 ; 
  int *iwork ;

  n = numlines(vname) ;
  ZALLOC(names, n, char *) ;
  xx = initarray_2Ddouble(1, n, 0) ;
  n = getxxnames(&names, xx, n, 1, vname) ;

  vzero(vals, numsnps) ;
  ZALLOC(iwork, numsnps, int) ;
  
  for (t=0; t<n; ++t) {
   sx = names[t]  ;
   k = snpindex(snpmarkers, numsnps, sx) ;
// if (t<10) printf("zzv %s %d %d %9.3f\n", sx, t, k, xx[0][t])  ;  

   if (k<0) { 
    continue ;
   }
   vals[k] = xx[0][t] ;
   iwork[k] = 1 ;
  }
  for (k=0; k<numsnps; ++k) { 
   cupt = snpmarkers[k] ;
   if (iwork[k] == 0) cupt -> ignore = YES ;
   if (cupt -> ignore == NO) ++nvalid ; 
   if (cupt -> ignore) vals[k] = 0 ; 
  }

  freeup(names, n) ;
  free(iwork) ;
  free2D(&xx, 1) ;
  return nvalid ;
}

void setref(char **preffasta, char *iubfile, char *tag) 
{     

 char *reffasta, **rrr  ;
 int t ; 
 reffasta = *preffasta ; 
 if (reffasta != NULL) return ;
 ZALLOC(rrr, 1, char *) ; 
 if (iubfile == NO) fatalx("no database\n") ;
 t = getfalist(&tag, 1, iubfile, rrr) ; 
 if (t < 0) fatalx("no %s tag in database %s\n", tag, iubfile) ;
 *preffasta = strdup(rrr[0]) ;
 freeup(rrr, 1) ; 
}





void setplimit (Indiv ** indivmarkers, int numindivs,
           char **eglist, int numeg, int plimit)
{
  int *indnums;
  int *psize;
  int i, k, kk;
  Indiv *indx;

  ZALLOC (indnums, numindivs, int);
  ZALLOC (psize, numeg, int);


  idperm (indnums, numindivs);
  ranperm (indnums, numindivs);

  for (i = 0; i < numindivs; i++) {
    k = indnums[i];
    indx = indivmarkers[k];
    if (indx->ignore)
      continue;
    kk = indxindex (eglist, numeg, indx->egroup);
    if (kk < 0)
      continue;
    ++psize[kk];
    if (psize[kk] > plimit)
      indx->ignore = YES;
  }



  free (psize);
  free (indnums);

}

int getrb(int pos, char *string, int len) 
// dangerous bend string is 1-based
{
 char  cx ;
 int pa = pos -1  ; 

 if (string == NULL) return -4 ;
 if (pa<0) return -3 ; 
 if (pa>=len) return -3 ; 

 cx = string[pa] ; 
 if (isdigit(cx)) return 1000 + cx -'0' ;
 
 return base2num(cx) ;
}