website/archive/binaries/mac/src/fasta34/res_stats.c

   1 /* calculate stats from results file using scalesws.c */
   2
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5 #include <string.h>
   6
   7 #include <limits.h>
   8 #include <math.h>
   9
  10 #define MAX_LLEN 200
  11
  12 #define LN_FACT 10.0
  13
  14 #include "defs.h"
  15 #include "structs.h"
  16 #include "param.h"
  17
  18 struct beststr {
  19   int score;    /* smith-waterman score */
  20   int sscore;   /* duplicate for compatibility with fasta */
  21   double comp;
  22   double H;
  23   double zscore;
  24   double escore;
  25   int n1;
  26 #ifndef USE_FTELLO
  27   long lseek;   /* position in library file */
  28 #else
  29   off_t lseek;
  30 #endif
  31   int cont;     /* offset into sequence */
  32   int frame;
  33   int lib;
  34   char libstr[13];
  35 } *bbp, *bestptr, **bptr, *best;
  36
  37 struct stat_str {
  38   int score;
  39   int n1;
  40   double comp;
  41   double H;
  42 };
  43
  44 static struct db_str qtt = {0l, 0l, 0};
  45
  46 char gstring2[MAX_STR];                  /* string for label */
  47 char gstring3[MAX_STR];
  48 char hstring1[MAX_STR];
  49
  50 FILE *outfd;
  51
  52 int nbest;      /* number of sequences better than bestcut in best */
  53 int bestcut=1;  /* cut off for getting into MAXBEST */
  54 int bestfull;
  55
  56 int dohist = 0;
  57 int zsflag = 1;
  58 int outtty=1;
  59 int llen=40;
  60
  61 /* statistics functions */
  62 extern void
  63 process_hist(struct stat_str *sptr, int nstat, struct pstruct pst,
  64              struct hist_str *hist, void **);
  65 extern void addhistz(double, struct hist_str *); /* scaleswn.c */
  66 void selectbestz(struct beststr **, int, int );
  67
  68 extern double zs_to_E(double, int, int, long, struct db_str);
  69 extern double zs_to_Ec(double zs, long entries);
  70
  71 extern double (*find_zp)(int score, int length, double comp, void *);
  72
  73 void prhist(FILE *, struct mngmsg, struct pstruct, struct hist_str,
  74             int, struct db_str, char *);
  75
  76 int nshow=20, mshow=50, ashow= -1;
  77 double e_cut=10.0;
  78
  79 main(argc, argv)
  80      int argc; char **argv;
  81 {
  82   FILE *fin;
  83   char line[512];
  84   int max, icol, iarg, i, qsfnum, lsfnum, n0, n1, s[3], frame;
  85   double comp, H;
  86   int idup, ndup, max_s;
  87   char libstr[20], *bp;
  88   char bin_file[80];
  89   FILE *bout=NULL;
  90   struct mngmsg m_msg;          /* Message from host to manager */
  91   struct pstruct pst;
  92   struct stat_str *stats;
  93   int nstats;
  94   double zscor, mu, var;
  95
  96 #if defined(UNIX)
  97   outtty = isatty(1);
  98 #else
  99   outtty = 1;
 100 #endif
 101
 102   if (argc < 2 ) {
 103     fprintf(stderr," useage - res_stats -c col -r bin_file file\n");
 104     exit(1);
 105   }
 106
 107   m_msg.db.length = qtt.length = 0l;
 108   m_msg.db.entries = m_msg.db.carry = qtt.entries = qtt.carry = 0;
 109   m_msg.pstat_void = NULL;
 110   m_msg.hist.hist_a = NULL;
 111   m_msg.nohist = 0;
 112   m_msg.markx = 0;
 113
 114   pst.n0 = 200;         /* sensible dummy value */
 115   pst.zsflag = 1;
 116   pst.dnaseq = 0;
 117   pst.histint = 2;
 118
 119   bin_file[0]='\0';
 120   icol = 1;
 121   iarg = 1;
 122   ndup = 1;
 123   while (1) {
 124     if (argv[iarg][0]=='-' && argv[iarg][1]=='c') {
 125       sscanf(argv[iarg+1],"%d",&icol);
 126       iarg += 2;
 127     }
 128     else if (argv[iarg][0]=='-' && argv[iarg][1]=='r') {
 129       strncpy(bin_file,argv[iarg+1],sizeof(bin_file));
 130       iarg += 2;
 131     }
 132     else if (argv[iarg][0]=='-' && argv[iarg][1]=='z') {
 133       sscanf(argv[iarg+1],"%d",&pst.zsflag);
 134       iarg += 2;
 135     }
 136     else if (argv[iarg][0]=='-' && argv[iarg][1]=='n') {
 137       pst.dnaseq = 1;
 138       iarg += 1;
 139     }
 140     else if (argv[iarg][0]=='-' && argv[iarg][1]=='s') {
 141       sscanf(argv[iarg+1],"%d",&ndup);
 142       iarg += 2;
 143     }
 144     else if (argv[iarg][0]=='-' && argv[iarg][1]=='q') {
 145       outtty = 0;
 146       iarg += 1;
 147     }
 148     else break;
 149   }
 150
 151   icol--;
 152
 153   if ((fin=fopen(argv[iarg],"r"))==NULL) {
 154     fprintf(stderr," cannot open %s\n",argv[1]);
 155     exit(1);
 156   }
 157
 158   if (bin_file[0]!='\0' && ((bout=fopen(bin_file,"w"))==NULL)) {
 159     fprintf(stderr,"cannot open %s for output\n",bin_file);
 160   }
 161
 162   if ((stats =
 163        (struct stat_str *)malloc((MAXSTATS)*sizeof(struct stat_str)))==NULL)
 164     s_abort ("Cannot allocate stats struct","");
 165   nstats = 0;
 166
 167   initbest(MAXBEST+1);  /* +1 required for select() */
 168
 169   for (nbest=0; nbest<MAXBEST+1; nbest++)
 170     bptr[nbest] = &best[nbest];
 171   bptr++; best++;
 172   best[-1].score= BIGNUM;
 173
 174   nbest = 0;
 175
 176   pst.Lambda=0.232;
 177   pst.K = 0.11;
 178   pst.H = 0.34;
 179
 180   /* read the best scores from the results file */
 181
 182   max_s = -1;
 183   idup = 0;
 184
 185   /* get first line with sequence length */
 186   fgets(line,sizeof(line),fin);
 187   sscanf(line,"%d",&n0);
 188   if (n0 > 0) pst.n0 = n0;
 189
 190   while (fgets(line,sizeof(line),fin)!=NULL) {
 191     if (line[0]=='/' && line[1]=='*') {
 192       fputs(line,stdout);
 193       strncpy(gstring2,line,sizeof(gstring2));
 194       if ((bp=strchr(gstring2,'\n'))!=NULL) *bp = '\0';
 195       break;
 196     }
 197     if (line[0]==';') {
 198       if ((bp=strchr(line,'|'))!=NULL) qsfnum = atoi(bp+1);
 199       else continue;
 200       if ((bp=strchr(line,'('))!=NULL) {
 201         n0 = atoi(bp+1);
 202         pst.n0 = n0;
 203       }
 204       else {
 205         fprintf(stderr, "cannot find n0:\n %s\n",line);
 206         continue;
 207       }
 208     }
 209     else {
 210         sscanf(line,"%s %d %d %d %lf %lf %d %d %d",
 211                libstr,&lsfnum,&n1,&frame,&comp, &H, &s[0],&s[1],&s[2]);
 212         if (lsfnum==0 && n1==0) {
 213           fputs(line,stderr);
 214           continue;
 215         }
 216         if (n1 < 10 || s[icol]<=0) fputs(line,stderr);
 217         idup++;
 218
 219         if (s[icol] > max_s) max_s = s[icol];
 220         if (idup < ndup) continue;
 221
 222         m_msg.db.entries++;
 223         m_msg.db.length += n1;
 224
 225         if (dohist) addhistz(zscor=(*find_zp)(max_s,n1,comp,m_msg.pstat_void),
 226                              &m_msg.hist);
 227         else zscor = (double)max_s;
 228
 229         if (nstats < MAXSTATS) {
 230           stats[nstats].n1 = n1;
 231           stats[nstats].comp = comp;
 232           stats[nstats].H = H;
 233           stats[nstats++].score = max_s;
 234         }
 235
 236         else if (!dohist) {
 237           /*      do_bout(bout,stats,nstats); */
 238           process_hist(stats,nstats,pst,&m_msg.hist, &m_msg.pstat_void);
 239           for (i=0; i<nbest; i++)
 240             bptr[i]->zscore =
 241               (*find_zp)(bptr[i]->score,bptr[i]->n1,bptr[i]->comp,
 242                          m_msg.pstat_void);
 243           dohist = 1;
 244         }
 245
 246         if (dohist) {
 247           zscor =(*find_zp)(max_s,n1,comp,m_msg.pstat_void);
 248           addhistz(zscor,&m_msg.hist);
 249         }
 250         else zscor = (double)max_s;
 251
 252         if (nbest >= MAXBEST) {
 253           bestfull = nbest-MAXBEST/4;
 254           selectz(bestfull-1,nbest);
 255           bestcut = (int)(bptr[bestfull-1]->zscore+0.5);
 256           nbest = bestfull;
 257         }
 258         bestptr = bptr[nbest];
 259         bestptr->score = max_s;
 260         bestptr->sscore = max_s;
 261         bestptr->n1 = n1;
 262         bestptr->comp = comp;
 263         bestptr->H = H;
 264         bestptr->lib = lsfnum;
 265         bestptr->zscore = zscor;
 266         strncpy(bestptr->libstr,libstr,12);
 267         bestptr->libstr[12]='\0';
 268         nbest++;
 269
 270         max_s = -1;
 271         idup = 0;
 272     }
 273   }     /* done with reading results */
 274
 275   if (!dohist) {
 276     if (nbest < 20) {
 277       zsflag = 0;
 278     }
 279     else {
 280       /*      do_bout(bout,stats,nstats); */
 281       process_hist(stats,nstats,pst,&m_msg.hist,&m_msg.pstat_void);
 282       for (i=0; i<nbest; i++)
 283         bptr[i]->zscore =
 284           (*find_zp)(bptr[i]->score,bptr[i]->n1,bptr[i]->comp,m_msg.pstat_void);
 285       dohist = 1;
 286     }
 287   }
 288
 289   printf(" using n0: %d\n",pst.n0);
 290
 291   /* print histogram, statistics */
 292
 293   m_msg.nbr_seq = m_msg.db.entries;
 294   pst.zdb_size = m_msg.db.entries;
 295   /* get_param(&pst, gstring2,gstring3); */
 296
 297   prhist(stdout,m_msg,pst,m_msg.hist,nstats,m_msg.db,gstring2);
 298
 299   if (!zsflag) sortbest();
 300   else {
 301     sortbestz(bptr,nbest);
 302     for (i=0; i<nbest; i++)
 303       bptr[i]->escore = zs_to_E(bptr[i]->zscore,bptr[i]->n1,pst.dnaseq,
 304                                 pst.zdb_size, m_msg.db);
 305   }
 306
 307   outfd = stdout;
 308   showbest(m_msg.db);   /* display best matches */
 309 }
 310
 311 initbest(nbest)         /* allocate arrays for best sort */
 312      int nbest;
 313 {
 314
 315   if ((best=(struct beststr *)calloc((size_t)nbest,sizeof(struct beststr)))
 316       == NULL) {fprintf(stderr,"cannot allocate best struct\n"); exit(1);}
 317   if ((bptr=(struct beststr **)calloc((size_t)nbest,sizeof(struct beststr *)))
 318       == NULL) {fprintf(stderr,"cannot allocate bptr\n"); exit(1);}
 319 }
 320
 321 void
 322 prhist(FILE *fd, struct mngmsg m_msg,
 323        struct pstruct pst,
 324        struct hist_str hist,
 325        int nstats,
 326        struct db_str ntt,
 327        char *gstring2)
 328 {
 329   int i,j,hl,hll, el, ell, ev;
 330   char hline[80], pch, *bp;
 331   int mh1, mht;
 332   int maxval, maxvalt, dotsiz, ddotsiz,doinset;
 333   double cur_e, prev_e, f_int;
 334   double max_dev, x_tmp;
 335   double db_tt;
 336   int n_chi_sq, cum_hl, max_i;
 337
 338
 339   fprintf(fd,"\n");
 340
 341   if (pst.zsflag < 0 || nstats <= 10) {
 342     fprintf(fd, "%7ld residues in %5ld sequences\n", ntt.length,ntt.entries);
 343     fprintf(fd,"\n%s\n",gstring2);
 344     return;
 345   }
 346
 347   max_dev = 0.0;
 348   mh1 = hist.maxh-1;
 349   mht = (3*hist.maxh-3)/4 - 1;
 350
 351   if (!m_msg.nohist && mh1 > 0) {
 352     for (i=0,maxval=0,maxvalt=0; i<hist.maxh; i++) {
 353       if (hist.hist_a[i] > maxval) maxval = hist.hist_a[i];
 354       if (i >= mht &&  hist.hist_a[i]>maxvalt) maxvalt = hist.hist_a[i];
 355     }
 356     n_chi_sq = 0;
 357     cum_hl = -hist.hist_a[0];
 358     dotsiz = (maxval-1)/60+1;
 359     ddotsiz = (maxvalt-1)/50+1;
 360     doinset = (ddotsiz < dotsiz && dotsiz > 2);
 361
 362     if (pst.zsflag>=0)
 363       fprintf(fd,"       opt      E()\n");
 364     else
 365       fprintf(fd,"     opt\n");
 366
 367     prev_e =  zs_to_Ec((double)(hist.min_hist-hist.histint/2),hist.entries);
 368     for (i=0; i<=mh1; i++) {
 369       pch = (i==mh1) ? '>' : ' ';
 370       pch = (i==0) ? '<' : pch;
 371       hll = hl = hist.hist_a[i];
 372       if (pst.zsflag>=0) {
 373         cum_hl += hl;
 374         f_int = (double)(i*hist.histint+hist.min_hist)+(double)hist.histint/2.0;
 375         cur_e = (double)zs_to_Ec(f_int,hist.entries);
 376         ev = el = ell = (int)(cur_e - prev_e + 0.5);
 377         if (hl > 0  && i > 5 && i < (90-hist.min_hist)/hist.histint) {
 378           x_tmp  = fabs(cum_hl - cur_e);
 379           if ( x_tmp > max_dev) {
 380             max_dev = x_tmp;
 381             max_i = i;
 382           }
 383           n_chi_sq++;
 384         }
 385         if ((el=(el+dotsiz-1)/dotsiz) > 60) el = 60;
 386         if ((ell=(ell+ddotsiz-1)/ddotsiz) > 40) ell = 40;
 387         fprintf(fd,"%c%3d %5d %5d:",
 388                 pch,(i<mh1)?(i)*hist.histint+hist.min_hist :
 389                 mh1*hist.histint+hist.min_hist,hl,ev);
 390       }
 391       else fprintf(fd,"%c%3d %5d :",
 392                    pch,(i<mh1)?(i)*hist.histint+hist.min_hist :
 393                    mh1*hist.histint+hist.min_hist,hl);
 394
 395       if ((hl=(hl+dotsiz-1)/dotsiz) > 60) hl = 60;
 396       if ((hll=(hll+ddotsiz-1)/ddotsiz) > 40) hll = 40;
 397       for (j=0; j<hl; j++) hline[j]='=';
 398       if (pst.zsflag>=0) {
 399         if (el <= hl ) {
 400           if (el > 0) hline[el-1]='*';
 401           hline[hl]='\0';
 402         }
 403         else {
 404           for (j = hl; j < el; j++) hline[j]=' ';
 405           hline[el-1]='*';
 406           hline[hl=el]='\0';
 407         }
 408       }
 409       else hline[hl] = 0;
 410       if (i==1) {
 411         for (j=hl; j<10; j++) hline[j]=' ';
 412         sprintf(&hline[10]," one = represents %d library sequences",dotsiz);
 413       }
 414       if (doinset && i == mht-2) {
 415         for (j = hl; j < 10; j++) hline[j]=' ';
 416         sprintf(&hline[10]," inset = represents %d library sequences",ddotsiz);
 417       }
 418       if (i >= mht&& doinset ) {
 419         for (j = hl; j < 10; j++) hline[j]=' ';
 420         hline[10]=':';
 421         for (j = 11; j<11+hll; j++) hline[j]='=';
 422         hline[11+hll]='\0';
 423         if (pst.zsflag>=0) {
 424           if (ell <= hll) hline[10+ell]='*';
 425           else {
 426             for (j = 11+hll; j < 10+ell; j++) hline[j]=' ';
 427             hline[10+ell] = '*';
 428             hline[11+ell] = '\0';
 429           }
 430         }
 431       }
 432
 433       fprintf(fd,"%s\n",hline);
 434       prev_e = cur_e;
 435     }
 436   }
 437
 438   if (ntt.carry==0) {
 439     fprintf(fd, "%7ld residues in %5ld sequences\n", ntt.length, ntt.entries);
 440   }
 441   else {
 442     db_tt = (double)ntt.carry*(double)LONG_MAX + (double)ntt.length;
 443     fprintf(fd, "%.0f residues in %5ld library sequences\n", db_tt, ntt.entries);
 444   }
 445
 446   if (pst.zsflag>=0) {
 447     if (MAXSTATS < hist.entries)
 448       fprintf(fd," statistics extrapolated from %d to %ld sequences\n",
 449               MAXSTATS,hist.entries);
 450     /*    summ_stats(stat_info); */
 451     fprintf(fd," %s\n",hist.stat_info);
 452     if (!m_msg.nohist && cum_hl > 0)
 453       fprintf(fd," Kolmogorov-Smirnov  statistic: %6.4f (N=%d) at %3d\n",
 454               max_dev/(double)cum_hl, n_chi_sq,max_i*hist.histint+hist.min_hist);
 455     if (m_msg.markx & MX_M10FORM) {
 456       while ((bp=strchr(hist.stat_info,'\n'))!=NULL) *bp=' ';
 457       if (cum_hl <= 0) cum_hl = -1;
 458       sprintf(hstring1,"; mp_extrap: %d %ld\n; mp_stats: %s\n; mp_KS: %6.4f (N=%d) at %3d\n",
 459               MAXSTATS,hist.entries,hist.stat_info,max_dev/(double)cum_hl, n_chi_sq,max_i*hist.histint+hist.min_hist);
 460     }
 461   }
 462   fprintf(fd,"\n%s\n",gstring2);
 463   fflush(fd);
 464 }
 465
 466 showbest(struct db_str ntt)
 467   {
 468     int ib, istart, istop;
 469     char bline[200], fmt[40], pad[200];
 470     char rline[20];
 471     int ntmp;
 472     int lcont, ccont, loff;
 473     int hcutoff;
 474
 475     sprintf(fmt,"%%-%ds (%%3d)",llen-10);
 476
 477     nshow = min(20,nbest);
 478     mshow = min(20,nbest);
 479
 480     if (outtty) {
 481       printf(" How many scores would you like to see? [%d] ",nshow);
 482       fflush(stdout);
 483       if (fgets(rline,sizeof(rline),stdin)==NULL) exit(0);
 484       if (rline[0]!='\n' && rline[0]!=0) sscanf(rline,"%d",&nshow);
 485       if (nshow<=0) nshow = min(20,nbest);
 486     }
 487     else nshow=mshow;
 488
 489     memset(pad,' ',llen-10);
 490     pad[llen-31]='\0';
 491     if (zsflag)
 492       fprintf(outfd,"The best scores are:%s s-w Z-score E(%ld)\n",pad,ntt.entries);
 493     else
 494       fprintf(outfd,"The best scores are:%s s-w\n",pad);
 495
 496     if (outfd != stdout)
 497       if (zsflag)
 498         fprintf(stdout,"The best scores are:%s s-w Z-score E(%ld)\n",pad,ntt.entries);
 499       else
 500         fprintf(stdout,"The best scores are:%s s-w\n",pad);
 501
 502     istart = 0;
 503   l1:   istop = min(nbest,nshow);
 504   for (ib=istart; ib<istop; ib++) {
 505     bbp = bptr[ib];
 506
 507     if (!outtty && zsflag && bbp->escore > e_cut) {
 508       nshow = ib;
 509       goto done;
 510     }
 511
 512     sprintf(bline,"%-12s %d",bbp->libstr,bbp->lib);
 513     bline[13]='\0';
 514
 515     fprintf(outfd,fmt,bline,bbp->n1);
 516
 517     if (zsflag)
 518       fprintf(outfd,"%4d %4.1f %6.2g\n",
 519               bbp->score,bbp->zscore,
 520               bbp->escore);
 521     else
 522       fprintf(outfd,"%4d\n",bbp->score);
 523
 524     if (outfd!=stdout) {
 525       fprintf(stdout,fmt,bline,bbp->n1);
 526       if (zsflag)
 527         printf("%4d %4.1f %6.2g\n",
 528                bbp->score,bbp->zscore,
 529                bbp->escore);
 530       else
 531         printf("%4d\n",bbp->score);
 532     }
 533   }
 534
 535   fflush(outfd); if (outfd!=stdout) fflush(stdout);
 536
 537   if (outtty) {
 538     printf(" More scores? [0] ");
 539     fflush(stdout);
 540     if (fgets(rline,sizeof(rline),stdin)==NULL) exit(0);
 541     ntmp = 0;
 542     if (rline[0]!='\n' && rline[0]!=0) sscanf(rline,"%d",&ntmp);
 543     if (ntmp<=0) ntmp = 0;
 544     if (ntmp>0) {
 545       istart = istop;
 546       nshow += ntmp;
 547       mshow += ntmp;
 548       goto l1;
 549     }
 550   }
 551   else if (zsflag && bbp->escore < e_cut) {
 552     istart=istop;
 553     nshow += 10;
 554     goto l1;
 555   }
 556
 557   done:
 558   if (outfd!=stdout) fprintf(outfd,"\n");
 559 }
 560
 561 selectz(k,n)    /* k is rank in array */
 562      int k,n;
 563 {
 564   int t, i, j, l, r;
 565   double v;
 566   struct beststr *tmptr;
 567
 568   l=0; r=n-1;
 569
 570   while ( r > l ) {
 571     i = l-1;
 572     j = r;
 573     v = bptr[r]->zscore;
 574     do {
 575       while (bptr[++i]->zscore > v ) ;
 576       while (bptr[--j]->zscore < v ) ;
 577       tmptr = bptr[i]; bptr[i]=bptr[j]; bptr[j]=tmptr;
 578     } while (j > i);
 579     bptr[j]=bptr[i]; bptr[i]=bptr[r]; bptr[r]=tmptr;
 580     if (i>=k) r = i-1;
 581     if (i<=k) l = i+1;
 582   }
 583 }
 584
 585 sortbest()
 586 {
 587   int cmps(), cmp1(), cmpa(), cmpz();
 588   ksort(bptr,nbest,cmps);
 589 }
 590
 591 sortbeste()
 592 {
 593   int cmpe();
 594   ksort(bptr,nbest,cmpe);
 595 }
 596
 597 sortbestz()
 598 {
 599   int cmpz();
 600   ksort(bptr,nbest,cmpz);
 601 }
 602
 603 cmps(ptr1,ptr2)
 604      struct beststr *ptr1, *ptr2;
 605 {
 606   if (ptr1->score < ptr2->score) return (1);
 607   else if (ptr1->score > ptr2->score) return (-1);
 608   else return (0);
 609 }
 610
 611 cmpe(ptr1,ptr2)
 612      struct beststr *ptr1, *ptr2;
 613 {
 614   if (ptr1->escore < ptr2->escore) return (-1);
 615   else if (ptr1->escore > ptr2->escore) return (1);
 616   else return (0);
 617 }
 618
 619 cmpz(ptr1,ptr2)
 620      struct beststr *ptr1, *ptr2;
 621 {
 622   if (ptr1->zscore < ptr2->zscore) return (1);
 623   else if (ptr1->zscore > ptr2->zscore) return (-1);
 624   else return (0);
 625 }
 626
 627 ksort(v,n,comp)
 628      char *v[]; int n, (*comp)();
 629 {
 630   int gap, i, j;
 631   char *tmp;
 632
 633   for (gap=n/2; gap>0; gap/=2)
 634     for (i=gap; i<n; i++)
 635       for (j=i-gap; j>=0; j -= gap) {
 636         if ((*comp)(v[j],v[j+gap]) <=0)
 637           break;
 638         tmp = v[j]; v[j]=v[j+gap]; v[j+gap]=tmp;
 639       }
 640 }
 641
 642 /*
 643 do_bout(FILE *bout,struct stat_str **bptr, int nbest)
 644 {
 645   int i, min_hist, max_hist;
 646   double mu, var;
 647
 648   if (bout==NULL) return;
 649
 650   inithist();
 651   for (i = 0; i<nbest; i++)
 652     addhist(bptr[i]->score,bptr[i]->n1);
 653
 654   for (i=0; i<MAX_LLEN; i++)
 655     if (llen_hist[i]>0) {
 656       min_hist=i;
 657       break;
 658     }
 659
 660   for (i=MAX_LLEN-1; i>=0; i--)
 661     if (llen_hist[i]>0) {
 662       max_hist=i;
 663       break;
 664     }
 665
 666   for (i=min_hist; i<=max_hist; i++) {
 667     mu=(double)score_sums[i]/(double)llen_hist[i];
 668     if (llen_hist[i]>1) {
 669       var = ((double)score2_sums[i]-(double)llen_hist[i]*mu*mu)/
 670         (double)(llen_hist[i]-1);
 671
 672       fprintf(bout,"%d\t%d\t%.1f\t%.1f\t%.1f\t%.4f\t%.4f\n",
 673               i,llen_hist[i],exp(((double)(i))/LN_FACT),
 674               score_sums[i],score2_sums[i],mu,var);
 675     }
 676   }
 677   free_hist();
 678   fclose(bout);
 679 }
 680 */
 681
 682 s_abort()
 683 {
 684   exit(1);
 685 }