binaries/src/fasta34/workacc.c

   1
   2 /* copyright (c) 1996, 1997, 1998, 1999 William R. Pearson and the
   3    U. of Virginia */
   4
   5 /* $Name: fa_34_26_5 $ - $Id: workacc.c,v 1.19 2006/02/07 17:58:19 wrp Exp $ */
   6
   7 /* Concurrent read version */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <string.h>
  12
  13 #include "param.h"
  14
  15 #define XTERNAL
  16 #include "uascii.h"
  17 #include "upam.h"
  18 #undef XTERNAL
  19
  20 char err_str[128];
  21
  22 /* Initialization - set up defaults - assume protein sequence */
  23 void w_init ()
  24 {
  25   pascii=aascii;
  26 }
  27
  28 #ifndef MPI_SRC
  29 /* allocate memory for pam matrix - identical to version in initfa/sw.c */
  30 alloc_pam (int d1, int d2, struct pstruct *ppst)
  31 {
  32   int     i, *d2p;
  33   char err_str[128];
  34
  35   if ((ppst->pam2[0] = (int **) malloc (d1 * sizeof (int *))) == NULL) {
  36      sprintf(err_str,"Cannot allocate 2D pam matrix: %d",d1);
  37      return -1;
  38   }
  39
  40   if ((ppst->pam2[1] = (int **) malloc (d1 * sizeof (int *))) == NULL) {
  41      sprintf(err_str,"Cannot allocate 2D pam matrix: %d",d1);
  42      return -1;
  43   }
  44
  45   if ((d2p = pam12 = (int *) malloc (d1 * d2 * sizeof (int))) == NULL) {
  46      sprintf(err_str,"Cannot allocate 2D pam matrix: %d",d1);
  47      return -1;
  48    }
  49    for (i = 0; i < d1; i++, d2p += d2) ppst->pam2[0][i] = d2p;
  50
  51    if ((d2p=pam12x= (int *) malloc (d1 * d2 * sizeof (int))) == NULL) {
  52      sprintf(err_str,"Cannot allocate 2d pam matrix: %d",d2);
  53      return -1;
  54    }
  55    for (i = 0;  i < d1; i++, d2p += d2) ppst->pam2[1][i] = d2p;
  56
  57    return 1;
  58 }
  59
  60 int **
  61 alloc_pam2p(int len, int nsq) {
  62   int i;
  63   int **pam2p;
  64
  65   if ((pam2p = (int **)calloc(len,sizeof(int *)))==NULL) {
  66     fprintf(stderr," Cannot allocate pam2p: %d\n",len);
  67     return NULL;
  68   }
  69
  70   if((pam2p[0] = (int *)calloc((nsq+1)*len,sizeof(int)))==NULL) {
  71     fprintf(stderr, "Cannot allocate pam2p[0]: %d\n", (nsq+1)*len);
  72     free(pam2p);
  73     return NULL;
  74   }
  75
  76   for (i=1; i<len; i++) {
  77     pam2p[i] = pam2p[0] + (i*(nsq+1));
  78   }
  79
  80   return pam2p;
  81 }
  82
  83 void free_pam2p(int **pam2p) {
  84   if (pam2p) {
  85     free(pam2p[0]);
  86     free(pam2p);
  87   }
  88 }
  89
  90 void
  91 aancpy(char *to, char *from, int count, struct pstruct pst)
  92 {
  93   char *tp, *sq;
  94   int nsq;
  95
  96   tp=to;
  97
  98   if (pst.ext_sq_set) {
  99     nsq = pst.nsqx;
 100     sq = pst.sqx;
 101   }
 102   else {
 103     nsq = pst.nsq;
 104     sq = pst.sq;
 105   }
 106
 107   while (count-- && *from) {
 108     if (*from <= nsq) *tp++ = sq[*(from++)];
 109     else *tp++ = *from++;
 110   }
 111   *tp='\0';
 112 }
 113 #endif
 114
 115 /* copies from from to to shuffling */
 116
 117 void
 118 shuffle(unsigned char *from, unsigned char *to, int n)
 119 {
 120   int i,j; unsigned char tmp;
 121
 122   if (from != to) memcpy((void *)to,(void *)from,(size_t)n);
 123
 124   for (i=n; i>0; i--) {
 125     j = nrand(i);
 126     tmp = to[j];
 127     to[j] = to[i-1];
 128     to[i-1] = tmp;
 129   }
 130   to[n] = 0;
 131 }
 132
 133 /* this shuffle is for FASTS */
 134 /* convert ',' -> '\0', shuffle each of the substrings */
 135 qshuffle(unsigned char *aa0, int n0, int nm0)
 136 {
 137   unsigned char **aa0start, *aap, tmp;
 138   int i,j,k, ns;
 139
 140   if ((aa0start=(unsigned char **)calloc(nm0+1,
 141                                          sizeof(unsigned char *)))==NULL) {
 142     fprintf(stderr,"cannot calloc for qshuffle %d\n",nm0);
 143     exit(1);
 144   }
 145   aa0start[0]=aa0;
 146   for (k=1,i=0; i<n0; i++) {
 147     if (aa0[i]==EOSEQ || aa0[i]==ESS) {
 148       aa0[i]='\0';
 149       aa0start[k++] = &aa0[i+1];
 150     }
 151   }
 152
 153   /* aa0start has the beginning of each substring */
 154   for (k=0; k<nm0; k++) {
 155     aap=aa0start[k];
 156     ns = strlen((char *)aap);
 157     for (i=ns; i>1; i--) {
 158       j = nrand(i);
 159       tmp = aap[j];
 160       aap[j] = aap[i-1];
 161       aap[i-1] = tmp;
 162     }
 163     aap[ns] = 0;
 164   }
 165
 166   for (k=1; k<nm0; k++) {
 167 /*  aap = aa0start[k];
 168     while (*aap) fputc(pst.sq[*aap++],stderr);
 169     fputc('\n',stderr);
 170 */
 171     aa0start[k][-1]=ESS;
 172   }
 173
 174   free(aa0start);
 175 }
 176
 177 /* copies from from to from shuffling */
 178 void
 179 wshuffle(unsigned char *from, unsigned char *to, int n, int wsiz, int *ieven)
 180 {
 181   int i,j, k, mm;
 182   unsigned char tmp, *top;
 183
 184   memcpy((void *)to,(void *)from,n);
 185
 186   mm = n%wsiz;
 187
 188   if (*ieven) {
 189     for (k=0; k<(n-wsiz); k += wsiz) {
 190       top = &to[k];
 191       for (i=wsiz; i>0; i--) {
 192         j = nrand(i);
 193         tmp = top[j];
 194         top[j] = top[i-1];
 195         top[i-1] = tmp;
 196       }
 197     }
 198     top = &to[n-mm];
 199     for (i=mm; i>0; i--) {
 200       j = nrand(i);
 201       tmp = top[j];
 202       top[j] = top[i-1];
 203       top[i-1] = tmp;
 204     }
 205     *ieven = 0;
 206   }
 207   else {
 208     for (k=n; k>=wsiz; k -= wsiz) {
 209       top = &to[k-wsiz];
 210       for (i=wsiz; i>0; i--) {
 211         j = nrand(i);
 212         tmp = top[j];
 213         top[j] = top[i-1];
 214         top[i-1] = tmp;
 215       }
 216     }
 217     top = &to[0];
 218     for (i=mm; i>0; i--) {
 219       j = nrand(i);
 220       tmp = top[j];
 221       top[j] = top[i-1];
 222       top[i-1] = tmp;
 223     }
 224     *ieven = 1;
 225   }
 226   to[n] = 0;
 227 }
 228
 229 void initseq(char **seqc0, char **seqc0a, char **seqc1, char **seqca, int seqsiz)       /* initialize arrays */
 230 {
 231   *seqc0=(char *)calloc((size_t)(seqsiz+1)*4,sizeof(char));
 232   *seqc0a= *seqc0+seqsiz+1;
 233   *seqc1= *seqc0a+seqsiz+1;
 234   *seqca= *seqc1+seqsiz+1;
 235   if (*seqc0==NULL)
 236     {fprintf(stderr,"cannot allocate consensus arrays %d\n",seqsiz);
 237      exit(1);}
 238 }
 239
 240 void freeseq(char **seqc0, char **seqc1, char **seqca)
 241 {
 242   free(*seqc0);
 243 }
 244
 245 #define ESS 49
 246
 247 void
 248 revcomp(unsigned char *seq, int n, int *c_nt)
 249 {
 250   unsigned char tmp;
 251   int i, ni;
 252
 253   for (i=0, ni = n-1; i< n/2; i++,ni--) {
 254     tmp = c_nt[seq[i]];
 255     seq[i] = c_nt[seq[ni]];
 256     seq[ni] = tmp;
 257   }
 258   if ((n%2)==1) {
 259     i = n/2;
 260     seq[i] = c_nt[seq[i]];
 261   }
 262   seq[n]=0;
 263 }