--- /dev/null
+/* copyright (c) 1996, 1997, 1998, 1999, 2002 William R. Pearson and the
+ U. of Virginia */
+
+/* $Name: fa_34_26_5 $ - $Id: comp_lib.c,v 1.100 2007/04/26 18:36:36 wrp Exp $ */
+
+/*
+ * Concurrent read version
+ *
+ * Feb 20, 1998 modifications for prss3
+ *
+ * December, 1998 - DNA searches are now down with forward and reverse
+ * strands
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+#include <time.h>
+
+#include <limits.h>
+#include <float.h>
+#include <math.h>
+
+#ifdef UNIX
+#include <unistd.h>
+#include <sys/types.h>
+#include <signal.h>
+#endif
+
+#include "defs.h"
+#include "mm_file.h"
+
+#include "mw.h" /* defines beststr */
+#include "structs.h" /* mngmsg, libstruct */
+#include "param.h" /* pstruct, thr_str, buf_head, rstruct */
+
+#define XTERNAL
+#include "uascii.h"
+
+char *mp_verstr="34.26";
+
+/********************************/
+/* global variable declarations */
+/********************************/
+char gstring2[MAX_STR]; /* string for label */
+char gstring3[MAX_STR];
+char hstring1[MAX_STR];
+
+extern int max_workers;
+
+#ifdef SUPERFAMNUM
+int nsfnum;
+int sfnum[10];
+extern int sfn_cmp(int *q, int *s);
+int nsfnum_n;
+int sfnum_n[10];
+#endif
+
+/********************************/
+/* extern variable declarations */
+/********************************/
+extern char *prog_func; /* function label */
+extern char *verstr, *iprompt0, *iprompt1, *iprompt2, *refstr;
+
+/********************************/
+/*extern function declarations */
+/********************************/
+/* open sequence file (getseq.c) */
+extern int getseq(char *filen, int *sascii,
+ unsigned char *seq, int maxs,
+ char *libstr, int n_libstr,
+ long *sq0ff);
+
+struct lmf_str *openlib(char *, int, int *, int, struct lmf_str *);
+
+void set_shuffle(struct mngmsg m_msg);
+void closelib(struct lmf_str *m_fptr);
+
+void irand(int);
+int nrand(int);
+
+extern int ann_scan(unsigned char *, int, struct mngmsg *, int );
+extern int scanseq(unsigned char *seq, int n, char *str);
+extern void re_ascii(int *qascii, int *sascii);
+extern int recode(unsigned char *seq, int n, int *qascii, int nsq);
+extern void revcomp(unsigned char *seq, int n, int *c_nt);
+
+extern void init_ascii(int is_ext, int *sascii, int is_dna);
+extern void qshuffle(unsigned char *aa0, int n0, int nm0);
+extern void free_pam2p(int **);
+
+/* initialize environment (doinit.c) */
+extern void initenv (int argc, char **argv, struct mngmsg *m_msg,
+ struct pstruct *ppst, unsigned char **aa0);
+
+/* print timing information */
+extern void ptime (FILE *, time_t);
+
+#ifdef COMP_MLIB
+#define QGETLIB (q_file_p->getlib)
+#endif
+
+#define GETLIB (m_file_p->getlib)
+
+/* calculation functions */
+extern void
+init_work(unsigned char *aa0, int n0,
+ struct pstruct *ppst, void **f_arg );
+#ifndef COMP_THR
+extern void
+do_work(unsigned char *aa0, int n0, unsigned char *aa1, int n1, int frame,
+ struct pstruct *ppst, void *f_str, int qr_flg, struct rstruct *rst);
+#endif
+
+extern void
+close_work(unsigned char *aa0, int n0, struct pstruct *ppst, void **f_arg);
+extern void
+get_param (struct pstruct *pstr, char *pstring1, char *pstring2);
+
+#ifdef COMP_THR
+#ifndef PRSS
+void
+save_best(struct buf_head *cur_buf, struct mngmsg, struct pstruct pst,
+ FILE *fdata, int *, struct hist_str *, void **);
+#else
+void
+save_best(struct buf_head *cur_buf, struct mngmsg, struct pstruct pst,
+ FILE *fdata, int *, struct hist_str *, void **, int *, int *);
+#endif
+#endif
+
+/* statistics functions */
+extern int
+process_hist(struct stat_str *sptr, int nstat,
+ struct mngmsg m_msg,
+ struct pstruct pst,
+ struct hist_str *hist, void **, int);
+extern void addhistz(double, struct hist_str *); /* scaleswn.c */
+void selectbestz(struct beststr **, int, int );
+extern double (*find_zp)(int score, double escore, int length, double comp,void *);
+
+void last_stats(const unsigned char *, int,
+ struct stat_str *sptr, int nstats,
+ struct beststr **bestp_arr, int nbest,
+ struct mngmsg m_msg, struct pstruct pst,
+ struct hist_str *histp, void *);
+
+int last_calc( unsigned char **a0, unsigned char *a1, int maxn,
+ struct beststr **bestp_arr, int nbest,
+ struct mngmsg m_msg, struct pstruct *ppst,
+ void **f_str, void *rs_str);
+
+void scale_scores(struct beststr **bestp_arr, int nbest,
+ struct db_str,struct pstruct pst, void *);
+
+#ifndef COMP_THR
+extern int shuffle(unsigned char *, unsigned char *, int);
+extern int wshuffle(unsigned char *, unsigned char *, int, int, int *);
+#endif
+
+extern void set_db_size(int, struct db_str *, struct hist_str *);
+
+/* display functions */
+extern void
+showbest (FILE *fp, unsigned char **aa0, unsigned char *aa1,
+ int maxn, struct beststr **bestp_arr, int nbest,
+ int qlib, struct mngmsg *m_msg,struct pstruct pst,
+ struct db_str db, char *gstring2, void **f_str);
+
+extern void
+showalign (FILE *fp, unsigned char **aa0, unsigned char *aa1,
+ int maxn, struct beststr **bestp_arr, int nbest,
+ int qlib, struct mngmsg m_msg,struct pstruct pst,
+ char *gstring2, void **f_str);
+
+/* misc functions */
+void h_init(struct pstruct *, struct mngmsg *, char *); /* doinit.c */
+void last_init(struct mngmsg *, struct pstruct *); /* initfa/sw.c */
+void last_params(unsigned char *, int, struct mngmsg *, struct pstruct *);
+
+void s_abort(char *, char *); /* compacc.c */
+
+/* initfa/sw.c */
+void resetp(struct mngmsg *, struct pstruct *);
+
+void gettitle(char *, char *, int); /* nxgetaa.c */
+void libchoice(char *lname, int, struct mngmsg *); /* lib_sel.c */
+void libselect(char *lname, struct mngmsg *); /* lib_sel.c */
+void query_parm(struct mngmsg *, struct pstruct *); /* initfa/sw.c */
+void selectbestz(struct beststr **, int, int);
+
+/* compacc.c */
+void prhist(FILE *, struct mngmsg, struct pstruct,
+ struct hist_str hist, int nstats, struct db_str, char *);
+void printsum(FILE *, struct db_str db);
+int reset_maxn(struct mngmsg *, int); /* set m_msg.maxt, maxn from maxl */
+
+FILE *outfd; /* Output file */
+
+/* this information is global for fsigint() */
+extern time_t s_time(); /* fetches time */
+time_t tstart, tscan, tprev, tdone; /* Timing */
+#ifdef COMP_MLIB
+time_t ttscan, ttdisp;
+#endif
+time_t tdstart, tddone;
+
+static struct db_str qtt = {0l, 0l, 0};
+
+#ifdef COMP_THR
+/***************************************/
+/* thread global variable declarations */
+/***************************************/
+
+/* functions for getting/sending buffers to threads (thr_sub.c) */
+extern void init_thr(int , struct thr_str *, struct mngmsg, struct pstruct *,
+ unsigned char *, int);
+extern void start_thr(void);
+extern void get_rbuf(struct buf_head **cur_buf, int max_wor_buf);
+extern void put_rbuf(struct buf_head *cur_buf, int max_work_buf);
+extern void put_rbuf_done(int nthreads, struct buf_head *cur_buf,
+ int max_work_buf);
+#undef XTERNAL
+#include "thr.h"
+struct buf_head buf_list[NUM_WORK_BUF];
+#endif
+
+/* these variables must be global for comp_thr.c so that savebest()
+ can use them */
+
+static struct beststr
+ *best, /* array of best scores */
+ *bestp,
+ **bestp_arr; /* array of pointers */
+static int nbest; /* number of best scores */
+
+static struct stat_str *stats, *qstats; /* array of scores for statistics */
+
+/* these variables are global so they can be set both by the main()
+ program and savebest() in threaded mode.
+*/
+static int nstats, nqstats, kstats;
+static double zbestcut; /* cut off for best z-score */
+static int bestfull; /* index for selectbest() */
+static int stats_done=0; /* flag for z-value processing */
+void fsigint();
+
+int
+main (int argc, char *argv[])
+{
+ unsigned char *aa0[6], *aa0s, *aa1, *aa1ptr, *aa1s;
+ int n1, n1s; /* n1s needed for PRSS so that when getlib() returns -1 (because no more
+ library sequences, we have a valid n1 for shuffling */
+
+ int *n1tot_ptr=NULL, *n1tot_cur;
+ int n1tot_cnt=0;
+ int n1tot_v, aa1_loff;
+
+ long qoffset; /* qoffset is the equivalent of loffset */
+ /* m_msg.sq0off is the l_off equivalent */
+
+ long loffset, l_off; /* loffset is the coordinate of first residue
+ when lcont > 0; l_off is not used in the
+ main loop, only in showbest and showalign */
+ char lib_label[MAX_FN];
+ char pgm_abbr[MAX_SSTR];
+ char qlabel[MAX_FN];
+#ifdef COMP_MLIB
+ char q_bline[MAX_STR];
+ fseek_t qseek;
+ int qlib;
+ struct lmf_str *q_file_p;
+ int sstart, sstop, is;
+#endif
+ int id;
+ struct lmf_str *m_file_p;
+
+ int t_best, t_rbest, t_qrbest; /* best score of two/six frames */
+ double t_escore, t_rescore, t_qrescore; /* best evalues of two/six frames */
+ int i_score;
+#ifdef PRSS
+ int s_score[3];
+ int s_n1;
+#endif
+
+ struct pstruct pst;
+ void *f_str[6], *qf_str; /* different f_str[]'s for different
+ translation frames, or forward,reverse */
+ int have_f_str=0;
+
+#ifdef COMP_THR
+ long ntbuff;
+ int max_buf_cnt, ave_seq_len, buf_siz;
+ int max_work_buf;
+ struct buf_head *cur_buf;
+ struct buf_str *cur_buf_p;
+ int nseq;
+ struct thr_str *work_info;
+#endif
+
+ struct mngmsg m_msg; /* Message from host to manager */
+ int iln, itt; /* index into library names */
+ char rline[MAX_FN];
+ char argv_line[MAX_STR];
+ int t_quiet;
+
+ struct rstruct rst; /* results structure */
+ struct rstruct rrst; /* results structure for shuffle*/
+ int i;
+
+ FILE *fdata=NULL; /* file for full results */
+ char libstr[MAX_UID]; /* string for labeling full results */
+ char *libstr_p; /* choose between libstr and ltitle */
+ int n_libstr; /* length of libstr */
+ int jstats;
+ int leng; /* leng is length of the descriptive line */
+ int maxn; /* size of the library sequence examined */
+ int maxl; /* size of library buffer */
+ fseek_t lmark; /* seek into library of current sequence */
+ int qlcont; /* continued query sequence */
+ int lcont, ocont, maxt; /* continued sequence */
+ int igncnt=0; /* count for ignoring sequences warning */
+ int ieven=0; /* tmp for wshuffle */
+ double zscore; /* tmp value */
+ char *bp; /* general purpose string ptr */
+
+ /* Initialization */
+
+#if defined(UNIX)
+ m_msg.quiet= !isatty(1);
+#else
+ m_msg.quiet = 0;
+#endif
+
+#ifdef PGM_DOC
+ argv_line[0]='#'; argv_line[1]='\0';
+ for (i=0; i<argc; i++) {
+ strncat(argv_line," ",sizeof(argv_line)-strlen(argv_line)-1);
+ if (strchr(argv[i],' ')) {
+ strncat(argv_line,"\"",sizeof(argv_line)-strlen(argv_line)-1);
+ strncat(argv_line,argv[i],sizeof(argv_line)-strlen(argv_line)-1);
+ strncat(argv_line,"\"",sizeof(argv_line)-strlen(argv_line)-1);
+ }
+ else {
+ strncat(argv_line,argv[i],sizeof(argv_line)-strlen(argv_line)-1);
+ }
+ }
+ argv_line[sizeof(argv_line)-1]='\0';
+#endif
+
+ /* first initialization routine - nothing is known */
+ h_init(&pst, &m_msg, pgm_abbr);
+
+ m_msg.db.length = qtt.length = 0l;
+ m_msg.db.entries = m_msg.db.carry = qtt.entries = qtt.carry = 0;
+ m_msg.pstat_void = NULL;
+ m_msg.hist.entries = 0;
+
+ for (iln=0; iln<MAX_LF; iln++) m_msg.lb_mfd[iln]=NULL;
+
+ f_str[0] = f_str[1] = NULL;
+
+ aa0[0] = NULL;
+ /* second initialiation - get commmand line arguments */
+ initenv (argc, argv, &m_msg, &pst,&aa0[0]);
+
+#ifdef COMP_THR
+ /* now have max_workers - allocate work_info[] */
+ if (max_workers >= MAX_WORKERS) max_workers = MAX_WORKERS;
+ if ((work_info=
+ (struct thr_str *)calloc(max_workers,sizeof(struct thr_str)))==NULL) {
+ fprintf(stderr, " cannot allocate work_info[%d]\n",max_workers);
+ exit(1);
+ }
+#else
+ max_workers = 1;
+#endif
+
+#ifndef PRSS
+ /* label library size limits */
+ if (m_msg.n1_low > 0 && m_msg.n1_high < BIGNUM)
+ sprintf(lib_label,"library (range: %d-%d)",m_msg.n1_low,m_msg.n1_high);
+ else if (m_msg.n1_low > 0)
+ sprintf(lib_label,"library (range: >%d)",m_msg.n1_low);
+ else if (m_msg.n1_high < BIGNUM)
+ sprintf(lib_label,"library (range: <%d)",m_msg.n1_high);
+ else
+ strncpy(lib_label,"library",sizeof(lib_label));
+#else
+ sprintf(lib_label,"shuffled sequence");
+#endif
+ lib_label[sizeof(lib_label)-1]='\0';
+
+ tstart = tscan = s_time();
+ tdstart = time(NULL);
+
+ /* Allocate space for the query and library sequences */
+ /* pad aa0[] with an extra 32 chars for ALTIVEC padding */
+ if (aa0[0]==NULL) {
+ if ((aa0[0] = (unsigned char *)malloc((m_msg.max_tot+1+32)*sizeof(unsigned char)))
+ == NULL)
+ s_abort ("Unable to allocate query sequence", "");
+ *aa0[0]=0;
+ aa0[0]++;
+ }
+ aa0[5]=aa0[4]=aa0[3]=aa0[2]=aa0[1]=aa0[0];
+
+ /* make room for random sequence -
+ also used as storage for COMP_THR library overlaps
+ */
+ if ((aa1s = (unsigned char *)malloc((m_msg.max_tot+1+32)*sizeof (char))) == NULL) {
+ s_abort ("Unable to allocate shuffled library sequence", "");
+ }
+ *aa1s=0;
+ aa1s++;
+
+ irand(0);
+
+ if (m_msg.markx & MX_HTML) {
+#ifdef HTML_HEAD
+ fprintf(stdout,"<html>\n<head>\n<title>%s Results</title>\n</head>\n<body>\n",prog_func);
+#endif
+ fprintf(stdout,"<pre>\n");
+ }
+
+#ifdef PGM_DOC
+ fputs(argv_line,stdout);
+ fputc('\n',stdout);
+#endif
+
+ fprintf(stdout,"%s\n",iprompt0);
+ fprintf(stdout," %s%s\n",verstr,refstr);
+ if (m_msg.markx & MX_HTML) fputs("</pre>\n",stdout);
+
+ /* Query library */
+ if (m_msg.tname[0] == '\0') {
+ if (m_msg.quiet == 1)
+ s_abort("Query sequence undefined","");
+ l1: fputs (iprompt1, stdout);
+ fflush (stdout);
+ if (fgets (m_msg.tname, MAX_FN, stdin) == NULL)
+ s_abort ("Unable to read query library name","");
+ m_msg.tname[MAX_FN-1]='\0';
+ if ((bp=strchr(m_msg.tname,'\n'))!=NULL) *bp='\0';
+ if (m_msg.tname[0] == '\0') goto l1;
+ }
+
+ /* Fetch first sequence */
+ qoffset = 0l;
+ qlcont = 0;
+#ifdef COMP_MLIB
+ /* Open query library */
+ if ((q_file_p= openlib(m_msg.tname, m_msg.qdnaseq,qascii,!m_msg.quiet,NULL))==NULL) {
+ s_abort(" cannot open library ",m_msg.tname);
+ }
+ qlib = 0;
+ m_msg.n0 =
+ QGETLIB (aa0[0], MAXTST, m_msg.qtitle, sizeof(m_msg.qtitle),
+ &qseek, &qlcont,q_file_p,&m_msg.sq0off);
+ if ((bp=strchr(m_msg.qtitle,' '))!=NULL) *bp='\0';
+ strncpy(qlabel,m_msg.qtitle,sizeof(qlabel));
+ if (bp != NULL) *bp = ' ';
+ qlabel[sizeof(qlabel)-1]='\0';
+
+ /* if annotations are included in sequence, remove them */
+ if (m_msg.ann_flg) {
+ m_msg.n0 = ann_scan(aa0[0],m_msg.n0,&m_msg,m_msg.qdnaseq);
+ }
+
+ if (m_msg.term_code && !(m_msg.qdnaseq==SEQT_DNA || m_msg.qdnaseq==SEQT_RNA) &&
+ aa0[0][m_msg.n0-1]!='*') {
+ aa0[0][m_msg.n0++]='*';
+ aa0[0][m_msg.n0]=0;
+ }
+
+ /* check for subset */
+ if (q_file_p->opt_text[0]!='\0') {
+ if (q_file_p->opt_text[0]=='-') {
+ sstart=0; sscanf(&q_file_p->opt_text[1],"%d",&sstop);
+ }
+ else {
+ sscanf(&q_file_p->opt_text[0],"%d-%d",&sstart,&sstop);
+ sstart--;
+ if (sstop <= 0 ) sstop = BIGNUM;
+ }
+ for (id=0,is=sstart; is<min(m_msg.n0,sstop); ) aa0[0][id++]=aa0[0][is++];
+ aa0[0][id]=0;
+ m_msg.n0 = min(m_msg.n0,sstop)-sstart;
+ if (m_msg.sq0off==1) m_msg.sq0off = sstart+1;
+ }
+
+#if defined(SW_ALTIVEC) || defined(SW_SSE2)
+ /* for ALTIVEC, must pad with 15 NULL's */
+ for (id=0; id<SEQ_PAD; id++) {aa0[0][m_msg.n0+id]=0;}
+#endif
+
+ if (qlcont) {
+ qoffset += m_msg.n0 - m_msg.sq0off;
+ }
+ else {
+ qoffset = 0l;
+ }
+
+#else
+ m_msg.n0 = getseq (m_msg.tname, qascii, aa0[0], m_msg.max_tot,
+ m_msg.qtitle, sizeof(m_msg.qtitle),
+ &m_msg.sq0off);
+ strncpy(qlabel,m_msg.tname,sizeof(qlabel));
+ qlabel[sizeof(qlabel)-1]='\0';
+
+ /* if annotations are included in sequence, remove them */
+ if (m_msg.ann_flg) {
+ m_msg.n0 = ann_scan(aa0[0],m_msg.n0,&m_msg,m_msg.qdnaseq);
+ }
+#endif
+
+ if (m_msg.n0 > MAXTST) {
+ fprintf(stderr," sequence truncated to %d\n %s\n",MAXTST,m_msg.sqnam);
+ fprintf(stdout," sequence truncated to %d\n %s\n",MAXTST,m_msg.sqnam);
+ aa0[0][MAXTST]='\0';
+ m_msg.n0=MAXTST;
+ }
+
+ if (m_msg.qdnaseq == SEQT_UNK) {
+
+ /* do automatic sequence recognition,but only for sequences > 20 residues */
+ if (m_msg.n0 > 20 &&
+ (float)scanseq(aa0[0],m_msg.n0,"ACGTUNacgtun")/(float)m_msg.n0 >0.85) {
+ pascii = nascii;
+ m_msg.qdnaseq = SEQT_DNA;
+ }
+ else { /* its protein */
+ pascii = aascii;
+ m_msg.qdnaseq = SEQT_PROT;
+ }
+ /* modify qascii to use encoded version
+ cannot use memcpy() because it loses annotations
+ */
+ re_ascii(qascii,pascii);
+ init_ascii(pst.ext_sq_set,qascii,m_msg.qdnaseq);
+ m_msg.n0 = recode(aa0[0],m_msg.n0,qascii, pst.nsqx);
+ }
+
+ if (m_msg.n0 <= 0)
+ s_abort ("Query sequence length <= 0: ", m_msg.tname);
+
+#ifdef SUPERFAMNUM
+ m_msg.nqsfnum = nsfnum;
+ for (i=0; i <= nsfnum & i<10; i++) m_msg.qsfnum[i] = sfnum[i];
+ m_msg.nqsfnum_n = nsfnum_n;
+ for (i=0; i <= nsfnum_n & i<10; i++) m_msg.qsfnum_n[i] = sfnum_n[i];
+#endif
+
+ resetp (&m_msg, &pst);
+
+#ifndef COMP_MLIB
+ gettitle(m_msg.tname,m_msg.qtitle,sizeof(m_msg.qtitle));
+ if (m_msg.tname[0]=='-' || m_msg.tname[0]=='@') {
+ strncmp(m_msg.tname,m_msg.qtitle,sizeof(m_msg.tname));
+ if ((bp=strchr(m_msg.tname,' '))!=NULL) *bp='\0';
+ }
+#endif
+
+ /* get library file names */
+
+#ifndef PRSS
+ if (strlen (m_msg.lname) == 0) {
+ if (m_msg.quiet == 1) s_abort("Library name undefined","");
+ libchoice(m_msg.lname,sizeof(m_msg.lname),&m_msg);
+ }
+
+ libselect(m_msg.lname, &m_msg);
+#else
+ if (strlen (m_msg.lname) == 0) {
+ if (m_msg.quiet == 1) s_abort("Shuffle sequence undefined","");
+l2: fputs(iprompt2,stdout);
+ fflush(stdout);
+ if (fgets (m_msg.lname, MAX_FN, stdin) == NULL)
+ s_abort ("Unable to read shuffle file name","");
+ m_msg.lname[MAX_FN-1]='\0';
+ if ((bp=strchr(m_msg.lname,'\n'))!=NULL) *bp='\0';
+ if (m_msg.lname[0] == '\0') goto l2;
+ }
+ m_msg.lbnames[0]= m_msg.lname;
+ m_msg.nln = 1;
+ m_msg.nshow = 0;
+#endif
+
+ /* Get additional parameters here */
+ if (!m_msg.quiet) query_parm (&m_msg, &pst);
+
+ last_init(&m_msg, &pst);
+
+ /* Allocate space for saved scores */
+ if ((best =
+ (struct beststr *)calloc((MAXBEST+1),sizeof(struct beststr)))==NULL)
+ s_abort ("Cannot allocate best struct","");
+ if ((bestp_arr =
+ (struct beststr **)malloc((MAXBEST+1)*sizeof(struct beststr *)))==NULL)
+ s_abort ("Cannot allocate bestp_arr","");
+
+ /* Initialize bestp_arr */
+ for (nbest = 0; nbest < MAXBEST+1; nbest++)
+ bestp_arr[nbest] = &best[nbest];
+ best++; bestp_arr++;
+ best[-1].score[0]=best[-1].score[1]=best[-1].score[2]= INT_MAX;
+ best[-1].zscore=FLT_MAX; /* for Z-scores, bigger is best */
+ best[-1].escore=FLT_MIN; /* for E()-values, lower is best */
+
+ if ((stats =
+ (struct stat_str *)calloc(MAXSTATS,sizeof(struct stat_str)))==NULL)
+ s_abort ("Cannot allocate stats struct","");
+
+#ifdef UNIX
+ /* set up signals now that input is done */
+ signal(SIGHUP,SIG_IGN);
+#endif
+
+#ifdef COMP_THR
+ /* Set up buffers for reading the library:
+
+ We will start by using a 2 Mbyte buffer for each worker. For
+ proteins, that means 5,000 sequences of length 400 (average).
+ For DNA, that means 2,000 sequences of length 1000. At the
+ moment, those are good averages.
+ */
+
+ if (m_msg.ldnaseq== SEQT_DNA) {
+ max_buf_cnt = MAX_NT_BUF;
+ ave_seq_len = AVE_NT_LEN;
+ }
+ else {
+ max_buf_cnt = MAX_AA_BUF;
+ ave_seq_len = AVE_AA_LEN;
+ }
+
+ /* however - buffer sizes should be a function of the number of
+ workers so that all the workers are kept busy. Assuming a 10,000
+ entry library is the smallest we want to schedule, then
+ */
+
+ if (max_buf_cnt > 10000/max_workers)
+ max_buf_cnt = 10000/(2*max_workers);
+
+ max_buf_cnt /= m_msg.thr_fact;
+
+ /* finally, max_work_buf should be mod 6 for tfasta */
+ max_buf_cnt -= (max_buf_cnt % 6);
+
+ max_work_buf = 2*max_workers;
+
+ /* allocate space for library buffers and results */
+
+ buf_siz=max_buf_cnt*ave_seq_len;
+ if (buf_siz < m_msg.max_tot) buf_siz = m_msg.max_tot;
+ for (i=0; i<max_work_buf; i++) {
+ if ((buf_list[i].buf =(struct buf_str *)calloc((size_t)(max_buf_cnt+1),
+ sizeof(struct buf_str)))
+ ==NULL) {
+ fprintf(stderr," cannot allocate buffer struct %d %d\n",i,max_buf_cnt+1);
+ exit(1);
+ }
+ buf_list[i].buf_cnt=0;
+ buf_list[i].have_results=0;
+ if ((buf_list[i].start =
+ (unsigned char *)calloc((size_t)(buf_siz),sizeof(unsigned char)))
+ ==NULL) {
+ fprintf(stderr," cannot allocate buffer %d\n",i);
+ exit(1);
+ }
+
+ /* make certain there is a '\0' at the beginning */
+ buf_list[i].start++;
+
+ reader_buf[i] = &buf_list[i];
+ }
+
+ /* initialization of global variables for threads/buffers */
+
+ num_worker_bufs = 0;
+ num_reader_bufs = max_work_buf;
+ reader_done = 0;
+ worker_buf_workp = 0;
+ worker_buf_readp = 0;
+ reader_buf_workp = 0;
+ reader_buf_readp = 0;
+
+ start_thread = 1; /* keeps threads from starting */
+#endif
+
+ /* Label the output */
+ if ((bp = (char *) strchr (m_msg.lname, ' ')) != NULL) *bp = '\0';
+ if (m_msg.ltitle[0] == '\0') {
+ strncpy(m_msg.ltitle,m_msg.lname,sizeof(m_msg.ltitle));
+ m_msg.ltitle[sizeof(m_msg.ltitle)-1]='\0';
+ }
+
+#ifdef COMP_MLIB
+ printf("Query library %s vs %s library\n", m_msg.tname,m_msg.lname);
+ if (m_msg.nln > 0) printf("searching %s library\n\n",m_msg.lbnames[0]);
+#endif
+
+#ifdef COMP_MLIB
+ while(1) {
+ m_msg.db.length = 0l;
+ m_msg.db.entries = m_msg.db.carry = 0;
+ qlib++;
+ stats_done = 0;
+#endif
+
+ maxl = m_msg.max_tot - m_msg.n0 -2; /* maxn = max library sequence space */
+
+ maxn = reset_maxn(&m_msg,maxl);
+ pst.maxlen = maxn;
+
+ outfd = stdout;
+ nbest = 0;
+ zbestcut = -FLT_MAX;
+ nstats = 0;
+
+ /* get the last parameters */
+ last_params(aa0[0],m_msg.n0, &m_msg, &pst);
+
+ /*
+ if our function returns approximate E()-scores, we do not need to
+ work with raw scores and later calculate z-scores. When
+ approx. E()-scores are calculated, we still need various
+ statistics structures, but we can get them immediately. In this
+ case, find_zp() must produce a z_score (large positive is good)
+ from an e_score.
+ */
+
+ if (m_msg.escore_flg) {
+ pst.zsflag_f = process_hist(stats,nstats,m_msg,pst,
+ &m_msg.hist,&m_msg.pstat_void,0);
+ stats_done=1;
+ }
+
+#ifndef COMP_THR
+ if (m_msg.qshuffle) {
+ if ((aa0s=(unsigned char *)calloc(m_msg.n0+2,sizeof(char)))==NULL) {
+ fprintf(stderr,"cannot allocate aa0s[%d]\n",m_msg.n0+2);
+ exit(1);
+ }
+ *aa0s='\0';
+ aa0s++;
+ memcpy(aa0s,aa0[0],m_msg.n0);
+ qshuffle(aa0s,m_msg.n0,m_msg.nm0);
+ }
+
+ /* previous versions of FASTA have stored the reverse complement in
+ the same array as the forward query sequence. This version
+ changes that, by allocating separate space for the reverse complement,
+ and thus reducing the demand for a large MAXLIB/MAXTRN for long queries
+ */
+ if (m_msg.qframe == 2) {
+ if ((aa0[1]=(unsigned char *)calloc(m_msg.n0+2,sizeof(char)))==NULL) {
+ fprintf(stderr,"cannot allocate aa0[1][%d]\n",m_msg.n0+2);
+ exit(1);
+ }
+ *aa0[1] = '\0';
+ aa0[1]++;
+ memcpy(aa0[1],aa0[0],m_msg.n0+1);
+ revcomp(aa0[1],m_msg.n0,&pst.c_nt[0]);
+ }
+ /* set aa1 for serial - threaded points aa1 to buffer */
+
+ aa1 = aa0[0] + m_msg.n0+1; /* modified now that aa0[1] is done separately */
+ *aa1++='\0';
+#else
+ init_thr(max_workers, work_info, m_msg, &pst, aa0[0], max_work_buf);
+#endif
+
+ if (m_msg.qshuffle && qstats==NULL) {
+ if ((qstats =
+ (struct stat_str *)calloc(m_msg.shuff_max+1,sizeof(struct stat_str)))==NULL)
+ s_abort ("Cannot allocate qstats struct","");
+ }
+ nqstats = 0;
+
+ if (m_msg.markx & MX_HTML) fputs("<pre>\n",stdout);
+#ifndef PRSS
+ /* rline[] is a tmp string */
+ if (m_msg.qdnaseq == SEQT_DNA || m_msg.qdnaseq == SEQT_RNA) {
+ strncpy(rline,(m_msg.qframe==1)? " (forward-only)" : "\0",sizeof(rline));
+ rline[sizeof(rline)-1]='\0';
+ }
+ else rline[0]='\0';
+
+ leng = (int)strlen(m_msg.qtitle);
+ if (leng > 50) leng -= 10;
+
+ sprintf (&m_msg.qtitle[leng], " %d %s", m_msg.n0, m_msg.sqnam);
+ m_msg.seqnm = 0;
+
+
+#ifdef COMP_MLIB
+ printf("%3d>>>%s - %d %s%s\n vs %.60s %s\n", qlib,
+ m_msg.qtitle, m_msg.n0, m_msg.sqnam,
+ (m_msg.revcomp ? " (reverse complement)" : rline),
+ m_msg.ltitle,lib_label);
+#else
+ printf("%.50s: %d %s%s\n %s\n vs %.60s %s\n",
+ qlabel, m_msg.n0, m_msg.sqnam,
+ (m_msg.revcomp ? " (reverse complement)" : rline),
+ m_msg.qtitle,m_msg.ltitle,lib_label);
+#endif
+ libstr_p = &libstr[0];
+ n_libstr=sizeof(libstr);
+#else /* PRSS */
+ libstr_p = &m_msg.ltitle[0];
+ n_libstr= sizeof(m_msg.ltitle);
+ set_shuffle(m_msg); /* set count/width parameters in llgetaa.c */
+#endif
+
+ fflush (outfd);
+
+ tprev = s_time();
+
+ if (m_msg.dfile[0] && (fdata=fopen(m_msg.dfile,"w"))!=NULL)
+ fprintf(fdata,"%3d\t%-50s\n",m_msg.n0,m_msg.qtitle);
+
+ qtt.length += m_msg.n0;
+ qtt.entries++;
+
+#ifdef COMP_THR
+ start_thr();
+
+ /* now open the library and start reading */
+ /* get a buffer and fill it up */
+ get_rbuf(&cur_buf,max_work_buf);
+
+ cur_buf->buf_cnt = 0;
+ cur_buf->have_results = 0;
+ cur_buf->buf[0].aa1b = cur_buf->start;
+ ntbuff = 0;
+ nseq = 0;
+#else /* ! COMP_THR */
+ /* initialize the comparison function, returning f_str */
+ init_work (aa0[0], m_msg.n0, &pst, &f_str[0]);
+ have_f_str=1;
+
+ f_str[5] = f_str[4] = f_str[3] = f_str[2] = f_str[1] = f_str[0];
+ if (m_msg.qframe == 2) {
+ init_work ( aa0[1], m_msg.n0, &pst, &f_str[1]);
+ }
+ if (m_msg.qshuffle) {
+ init_work ( aa0s, m_msg.n0, &pst, &qf_str);
+ }
+#endif /* COMP_THR */
+
+ /* open the library - start the search */
+
+ for (iln = 0; iln < m_msg.nln; iln++) {
+ if ((m_msg.lb_mfd[iln] = m_file_p=
+ openlib(m_msg.lbnames[iln], m_msg.ldnaseq, lascii, !m_msg.quiet, m_msg.lb_mfd[iln]))
+ ==NULL) {
+ fprintf(stderr," cannot open library %s\n",m_msg.lbnames[iln]);
+ continue;
+ }
+#if !defined(PRSS) && !defined(COMP_MLIB)
+ else
+ printf ("searching %s %s\n",m_msg.lbnames[iln],lib_label);
+#endif
+
+ loffset = 0l;
+ lcont = 0;
+ ocont = 0;
+ n1tot_v = n1tot_cnt = 0;
+ n1tot_cur = n1tot_ptr = NULL;
+
+ /* get next buffer to read into */
+ maxt = maxn;
+
+#ifndef COMP_THR
+ aa1ptr = aa1;
+#else
+ /* read sequence directly into buffer */
+ aa1ptr = aa1 = cur_buf->buf[nseq].aa1b;
+#endif
+
+ while ((n1=GETLIB(aa1ptr,maxt,libstr_p,n_libstr,&lmark,&lcont,m_file_p,&l_off))>=0) {
+
+ if (n_libstr <= MAX_UID) {
+ if ((bp=strchr(libstr_p,' '))!=NULL) *bp='\0';
+ }
+
+ if (m_msg.term_code && !lcont &&
+ m_msg.ldnaseq==SEQT_PROT && aa1ptr[n1-1]!=m_msg.term_code) {
+ aa1ptr[n1++]=m_msg.term_code;
+ aa1ptr[n1]=0;
+ }
+
+#if defined(SW_ALTIVEC) || defined(SW_SSE2)
+ /* for ALTIVEC, must pad with 15 NULL's */
+ for (id=0; id<SEQ_PAD; id++) {aa1ptr[n1+id]=0;}
+#endif
+
+#ifdef DEBUG
+ if (aa1[-1]!='\0' || aa1ptr[n1]!='\0') {
+ fprintf(stderr,"%s: aa1[%d] missing NULL boundaries: %d %d\n",libstr_p,n1,aa1[-1],aa1ptr[n1]);
+ }
+#endif
+
+ /* check for a continued sequence and provide a pointer to
+ the n1_tot array if lcont || ocont */
+ n1tot_v += n1;
+ if (lcont && !ocont) { /* get a new pointer */
+ if (n1tot_cnt <= 0) {
+ if ((n1tot_ptr=calloc(1000,sizeof(int)))==NULL) {
+ fprintf(stderr," cannot allocate n1tot_ptr\n");
+ exit(1);
+ }
+ else {n1tot_cnt=1000;}
+ }
+ n1tot_cnt--;
+ n1tot_cur = n1tot_ptr++;
+ }
+
+ if (n1tot_v < m_msg.n1_low || n1tot_v > m_msg.n1_high) {
+ goto loop2;
+ }
+
+ m_msg.db.entries++;
+ m_msg.db.length += n1;
+ if (m_msg.db.length > LONG_MAX) {
+ m_msg.db.length -= LONG_MAX; m_msg.db.carry++;
+ }
+
+#ifdef DEBUG
+ /* This finds most reasons for core dumps */
+ if (pst.debug_lib)
+ for (i=0; i<n1; i++)
+ if (aa1[i]>=pst.nsqx)
+ {fprintf(stderr,
+ "%s residue[%d/%d] %d range (%d) lcont/ocont: %d/%d\n%s\n",
+ libstr,i,n1,aa1[i],pst.nsqx,lcont,ocont,aa1ptr+i);
+ aa1[i]=0;
+ n1=i-1;
+ break;
+ }
+#endif
+
+ /* don't count long sequences more than once */
+ if (aa1!=aa1ptr) {n1 += m_msg.loff; m_msg.db.entries--;}
+
+#ifdef PROGRESS
+ if (!m_msg.quiet)
+ if (m_msg.db.entries % 200 == 199) {
+ fputc('.',stderr);
+ if (m_msg.db.entries % 10000 == 9999) fputc('\n',stderr);
+ else if (m_msg.db.entries % 1000 == 999) fputc(' ',stderr);
+
+ }
+#endif
+
+ if (n1<=1) {
+ /* if (igncnt++ <10)
+ fprintf(stderr,"Ignoring: %s\n",libstr);
+ */
+ goto loop2;
+ }
+
+#ifdef PRSS
+ if (lmark==0) {
+ n1s = n1;
+ memcpy(aa1s,aa1,n1s);
+ m_msg.db.entries=0;
+ m_msg.db.length=0;
+ }
+#endif
+
+ /* if COMP_THR - fill and empty buffers */
+#ifdef COMP_THR
+ ntbuff += n1+1;
+
+ for (itt=m_msg.revcomp; itt<=m_msg.nitt1; itt++) {
+
+ cur_buf->buf_cnt++;
+ cur_buf_p = &(cur_buf->buf[nseq++]);
+ cur_buf_p->n1 = n1;
+ cur_buf_p->n1tot_p = n1tot_cur;
+ cur_buf_p->lseek = lmark;
+ cur_buf_p->cont = ocont+1;
+ cur_buf_p->m_file_p = (void *)m_file_p;
+ cur_buf_p->frame = itt;
+ memcpy(cur_buf_p->libstr,libstr,MAX_UID);
+#ifdef SUPERFAMNUM
+ cur_buf_p->nsfnum = nsfnum;
+ if ((cur_buf_p->sfnum[0]=sfnum[0])>0 &&
+ (cur_buf_p->sfnum[1]=sfnum[1])>0 &&
+ (cur_buf_p->sfnum[2]=sfnum[2])>0 &&
+ (cur_buf_p->sfnum[3]=sfnum[3])>0 &&
+ (cur_buf_p->sfnum[4]=sfnum[4])>0 &&
+ (cur_buf_p->sfnum[5]=sfnum[5])>0 &&
+ (cur_buf_p->sfnum[6]=sfnum[6])>0 &&
+ (cur_buf_p->sfnum[7]=sfnum[7])>0 &&
+ (cur_buf_p->sfnum[8]=sfnum[8])>0 &&
+ (cur_buf_p->sfnum[9]=sfnum[9])>0) ;
+#endif
+
+ /* this assumes that max_buf_cnt is guaranteed %6=0 so that
+ additional pointers to the same buffer can be used
+ nseq now points to next buffer
+ */
+
+ cur_buf->buf[nseq].aa1b = cur_buf->buf[nseq-1].aa1b;
+ } /* for (itt .. */
+
+ /* make a copy of the overlap (threaded only) */
+ if (lcont) {
+ memcpy(aa1s,&aa1[n1-m_msg.loff],m_msg.loff);
+ }
+
+ /* if the buffer is filled */
+ if (nseq >= max_buf_cnt || ntbuff >= buf_siz - maxn) {
+
+ /* provide filled buffer to workers */
+ put_rbuf(cur_buf,max_work_buf);
+
+ /* get an empty buffer to fill */
+ get_rbuf(&cur_buf,max_work_buf);
+
+ /* "empty" buffers have results that must be processed */
+ if (cur_buf->buf_cnt && cur_buf->have_results) {
+ save_best(cur_buf,m_msg,pst,fdata,m_msg.qsfnum,&m_msg.hist,
+ &m_msg.pstat_void
+#ifdef PRSS
+ ,s_score,&s_n1
+#endif
+ );
+
+ }
+
+ /* now the buffer is truly empty, fill it up */
+ cur_buf->buf_cnt = 0;
+ cur_buf->have_results = 0;
+ /* point the first aa1 ptr to the buffer start */
+ aa1=cur_buf->buf[0].aa1b = cur_buf->start;
+ ntbuff = 0;
+ nseq=0;
+ }
+ else { /* room left in current buffer, increment ptrs */
+ aa1=cur_buf->buf[nseq].aa1b = cur_buf->buf[nseq-1].aa1b+n1+1;
+ }
+#else /* if !COMP_THR - do a bunch of searches */
+
+ /* t_best and t_rbest are used to save the best score or shuffled
+ score from all the frames */
+
+ t_best = t_rbest = t_qrbest = -1;
+ t_escore = t_rescore = t_qrescore = FLT_MAX;
+ for (itt=m_msg.revcomp; itt<=m_msg.nitt1; itt++) {
+
+ rst.score[0] = rst.score[1] = rst.score[2] = 0;
+ do_work (aa0[itt], m_msg.n0,aa1,n1,itt,&pst,f_str[itt],0,&rst);
+
+ if (rst.score[pst.score_ix] > t_best) {
+ t_best = rst.score[pst.score_ix];
+ }
+
+ if (fdata) {
+ fprintf(fdata,
+ "%-12s %5d %6d %d %.5f %.5f %4d %4d %4d %g %d %d %8lld\n",
+ libstr,
+#ifdef SUPERFAMNUM
+ sfn_cmp(m_msg.qsfnum,sfnum),
+#else
+ 0,
+#endif
+ n1,itt,
+ rst.comp,rst.H,
+ rst.score[0],rst.score[1],rst.score[2],
+ rst.escore, rst.segnum, rst.seglen, lmark);
+ fflush(fdata);
+ }
+
+#ifdef PRSS
+ if (lmark==0) {
+ s_score[0] = rst.score[0];
+ s_score[1] = rst.score[1];
+ s_score[2] = rst.score[2];
+
+ s_n1 = n1;
+ aa1_loff = l_off;
+ }
+ t_best = t_rbest = rst.score[pst.score_ix];
+ t_escore = t_rescore = rst.escore;
+#else
+ if (m_msg.qshuffle) {
+ do_work (aa0s, m_msg.n0,aa1,n1,itt,&pst,qf_str,1,&rrst);
+
+ if (rrst.score[pst.score_ix] > t_qrbest)
+ t_qrbest = rrst.score[pst.score_ix];
+ if (rrst.escore < t_qrescore)
+ t_qrescore = rrst.escore;
+
+ if (itt==m_msg.nitt1 && nqstats < m_msg.shuff_max) {
+ qstats[nqstats].n1 = n1; /* save the best score */
+ qstats[nqstats].comp = rst.comp;
+ qstats[nqstats].H = rst.H;
+ qstats[nqstats].escore = t_qrescore;
+ qstats[nqstats++].score = t_qrbest;
+ t_qrbest = -1; /* reset t_qrbest, t_qrescore */
+ t_qrescore = FLT_MAX;
+ }
+ }
+
+ if (pst.zsflag >= 10) {
+ if (pst.zs_win > 0) wshuffle(aa1,aa1s,n1,pst.zs_win,&ieven);
+ else shuffle(aa1,aa1s,n1);
+ do_work (aa0[itt], m_msg.n0, aa1s, n1,itt,&pst,f_str[itt],0,&rrst);
+ if (rrst.score[pst.score_ix] > t_rbest) {
+ t_rbest = rrst.score[pst.score_ix];
+ t_rescore = rrst.escore;
+ }
+ }
+#endif
+ i_score = rst.score[pst.score_ix];
+
+/* this section saves scores for statistics calculations. For
+ comparisons that can be from one of 2 or 6 frames, it should only
+ be run once, for the best of the 2 or 6 scores. t_rbest,t_rescore
+ have the best of the 2 or 6 scores from the frames. For proteins,
+ this is run for every score.
+
+*/
+#ifdef PRSS /* don't save the first score (unshuffled) with PRSS */
+ if (lmark > 0) {
+#endif
+
+ if (itt == m_msg.nitt1) {
+ if (nstats < MAXSTATS) {
+ stats[nstats].n1 = n1; /* save the best score */
+ stats[nstats].comp = rst.comp;
+ stats[nstats].H = rst.H;
+ if (pst.zsflag >=10) {
+ t_best = t_rbest;
+ t_escore = t_rescore;
+ }
+ stats[nstats].escore = t_escore;
+ stats[nstats++].score = t_best;
+ t_best = t_rbest = -1; /* reset t_rbest, t_best */
+ t_escore = t_rescore = FLT_MAX;
+ }
+ else if (pst.zsflag >= 0) {
+ if (!stats_done) {
+ pst.zsflag_f = process_hist(stats,nstats,m_msg,pst,
+ &m_msg.hist,&m_msg.pstat_void,0);
+ stats_done = 1;
+ kstats = nstats;
+ for (i=0; i<MAXBEST; i++) {
+ bestp_arr[i]->zscore =
+ (*find_zp)(bestp_arr[i]->score[pst.score_ix],
+ bestp_arr[i]->escore, bestp_arr[i]->n1,
+ bestp_arr[i]->comp, m_msg.pstat_void);
+ }
+ zbestcut = bestp_arr[nbest-1]->zscore;
+ }
+
+#ifdef SAMP_STATS
+/* older versions saved the first MAXSTATS scores, and ignored the
+ rest in the statistics. With SAMP_STATS, scores after MAX_STATS
+ are sampled at random, and included in the sample set and the
+ statistics parameters are re-derived at the end of the run using
+ the sampled scores.
+
+ It would be faster not to do the nrand(); if(jstats < MAXSTATS)
+ less often.
+*/
+ if (!m_msg.escore_flg) { /* only for zscores */
+ jstats = nrand(++kstats); /* no mod % 0 */
+ if (jstats < MAXSTATS) {
+ stats[jstats].n1 = n1; /* save the best score */
+ stats[jstats].comp = rst.comp;
+ stats[jstats].H = rst.H;
+ if (pst.zsflag >=10) t_best = t_rbest;
+ stats[jstats].score = t_best;
+ }
+ }
+#endif
+ } /* ( nstats >= MAXSTATS) && zsflag >= 0 */
+ } /* itt1 == nitt1 */
+#ifdef PRSS
+ }
+#endif
+
+ /* this section completes work on the current score */
+ if (stats_done) { /* stats_done > 0 => nstats >= MAXSTATS */
+ zscore=(*find_zp)(i_score, rst.escore, n1, rst.comp,
+ m_msg.pstat_void);
+
+ if (itt == m_msg.nitt1) {
+ if (pst.zsflag >= 10) t_best = t_rbest;
+
+ addhistz((*find_zp)(t_best, t_escore, n1, rst.comp,
+ m_msg.pstat_void),
+ &m_msg.hist);
+ t_best = t_rbest = -1;
+ }
+ }
+ else zscore = (double) i_score;
+
+#ifndef PRSS
+ if (zscore > zbestcut ) {
+ if (nbest >= MAXBEST) {
+ bestfull = nbest-MAXBEST/4;
+ selectbestz(bestp_arr,bestfull-1,nbest);
+ zbestcut = bestp_arr[bestfull-1]->zscore;
+ nbest = bestfull;
+ }
+
+ bestp = bestp_arr[nbest++];
+ bestp->score[0] = rst.score[0];
+ bestp->score[1] = rst.score[1];
+ bestp->score[2] = rst.score[2];
+ bestp->comp = rst.comp;
+ bestp->H = rst.H;
+ bestp->zscore = zscore;
+ bestp->escore = rst.escore;
+ bestp->segnum = rst.segnum;
+ bestp->seglen = rst.seglen;
+ bestp->lseek = lmark;
+ bestp->cont = ocont+1;
+ bestp->m_file_p = m_file_p;
+ bestp->n1 = n1;
+ bestp->n1tot_p=n1tot_cur;
+ bestp->frame = itt;
+ memcpy(bestp->libstr,libstr,MAX_UID);
+#ifdef SUPERFAMNUM
+ bestp->nsfnum = nsfnum;
+ if ((bestp->sfnum[0]=sfnum[0])>0 &&
+ (bestp->sfnum[1]=sfnum[1])>0 &&
+ (bestp->sfnum[2]=sfnum[2])>0 &&
+ (bestp->sfnum[3]=sfnum[3])>0 &&
+ (bestp->sfnum[4]=sfnum[4])>0 &&
+ (bestp->sfnum[5]=sfnum[5])>0 &&
+ (bestp->sfnum[6]=sfnum[6])>0 &&
+ (bestp->sfnum[7]=sfnum[7])>0 &&
+ (bestp->sfnum[8]=sfnum[8])>0 &&
+ (bestp->sfnum[9]=sfnum[9])>0) ;
+#endif
+ }
+#else /* PRSS */
+ if (lmark == 0) {
+ bestp = bestp_arr[nbest++];
+ bestp->score[0] = rst.score[0];
+ bestp->score[1] = rst.score[1];
+ bestp->score[2] = rst.score[2];
+ bestp->comp = rst.comp;
+ bestp->H = rst.H;
+ bestp->zscore = zscore;
+ bestp->escore = rst.escore;
+ bestp->segnum = rst.segnum;
+ bestp->seglen = rst.seglen;
+ bestp->lseek = lmark;
+ bestp->cont = 0;
+ bestp->m_file_p = m_file_p;
+ bestp->n1 = n1;
+ bestp->n1tot_p=n1tot_cur;
+ bestp->frame = itt;
+ memcpy(bestp->libstr,libstr,MAX_UID);
+ bestp->nsfnum = 0;
+ }
+#endif
+ }
+#endif
+
+ loop2:
+ if (lcont) {
+ maxt = m_msg.maxt3;
+#ifndef COMP_THR
+ memcpy(aa1,&aa1[n1-m_msg.loff],m_msg.loff);
+#else
+ memcpy(aa1,aa1s,m_msg.loff);
+#endif
+ aa1ptr= &aa1[m_msg.loff];
+ loffset += n1 - m_msg.loff;
+ ocont = lcont;
+ }
+ else {
+ maxt = maxn;
+ aa1ptr=aa1;
+ if (ocont) *n1tot_cur = n1tot_v;
+ ocont = 0;
+ loffset = 0l;
+ n1tot_v = 0;
+ n1tot_cur = NULL;
+ }
+ } /* end while((n1=getlib())) */
+ } /* end iln=1..nln */
+
+ /* all done */
+
+#ifdef COMP_THR
+ /* check last buffers for any results */
+ put_rbuf_done(max_workers,cur_buf,max_work_buf);
+
+ for (i=0; i < num_reader_bufs; i++) {
+ reader_buf_readp = (reader_buf_readp+1)%(max_work_buf);
+ if (reader_buf[reader_buf_readp]->buf_cnt > 0 &&
+ reader_buf[reader_buf_readp]->have_results) {
+ save_best(reader_buf[reader_buf_readp],m_msg,pst,fdata,m_msg.qsfnum,
+ &m_msg.hist, &m_msg.pstat_void
+#ifdef PRSS
+ ,s_score,&s_n1
+#endif
+ );
+ }
+ }
+#endif
+
+#ifdef PROGRESS
+ if (!m_msg.quiet)
+ if (m_msg.db.entries >= 200) {fprintf(stderr," Done!\n");}
+#endif
+
+ m_msg.nbr_seq = m_msg.db.entries;
+ get_param(&pst, gstring2,gstring3);
+
+/* *************************** */
+/* analyze the last results */
+/* *************************** */
+
+#ifndef PRSS
+#ifndef SAMP_STATS
+ if (!stats_done && nstats > 0) {
+#endif
+ pst.zsflag_f = process_hist(stats,nstats,m_msg,pst,&m_msg.hist,
+ &m_msg.pstat_void,stats_done);
+ if (m_msg.pstat_void != NULL) {
+ stats_done = 1;
+ for (i = 0; i < nbest; i++) {
+ bestp_arr[i]->zscore =
+ (*find_zp)(bestp_arr[i]->score[pst.score_ix],
+ bestp_arr[i]->escore, bestp_arr[i]->n1,
+ bestp_arr[i]->comp, m_msg.pstat_void);
+ }
+#ifndef SAMP_STATS
+ }
+ else pst.zsflag = -1;
+#endif
+ }
+#else /* PRSS */
+ if (pst.zsflag < 10) pst.zsflag += 10;
+ pst.zsflag_f = process_hist(stats,nstats,m_msg,pst,
+ &m_msg.hist, &m_msg.pstat_void,0);
+ stats_done = 1;
+ for (i = 0; i < nbest; i++) {
+ bestp_arr[i]->zscore = (*find_zp)(bestp_arr[i]->score[pst.score_ix],
+ bestp_arr[i]->escore, bestp_arr[i]->n1,
+ bestp_arr[i]->comp, m_msg.pstat_void);
+ }
+#endif
+
+ if (pst.zdb_size <= 1) pst.zdb_size = m_msg.db.entries;
+
+#ifdef COMP_THR
+ /* before I call last_calc/showbest/showalign, I need init_work() to
+ get an f_str. This duplicates some code above, which is used in
+ the non-threaded version
+ */
+
+ if (!have_f_str) {
+ init_work(aa0[0],m_msg.n0,&pst,&f_str[0]);
+ have_f_str = 1;
+ f_str[5] = f_str[4] = f_str[3] = f_str[2] = f_str[1] = f_str[0];
+
+ if (m_msg.qframe == 2) {
+ if ((aa0[1]=(unsigned char *)calloc((size_t)m_msg.n0+2,
+ sizeof(unsigned char)))==NULL) {
+ fprintf(stderr," cannot allocate aa0[1][%d] for alignments\n",
+ m_msg.n0+2);
+ }
+ *aa0[1]='\0';
+ aa0[1]++;
+ memcpy(aa0[1],aa0[0],m_msg.n0+1);
+ revcomp(aa0[1],m_msg.n0,&pst.c_nt[0]);
+ init_work(aa0[1],m_msg.n0,&pst,&f_str[1]);
+ }
+
+ /* I also need a "real" aa1 */
+ aa1 = buf_list[0].start;
+#ifdef PRSS
+ /* for PRSS - I need the original second (non-shuffled) sequence */
+ memcpy(aa1,aa1s,n1s+1);
+#endif
+ }
+#endif
+
+/* now we have one set of scaled scores for in bestp_arr -
+ for FASTS/F, we need to do some additional processing */
+
+ if (!m_msg.qshuffle) {
+ last_stats(aa0[0], m_msg.n0, stats,nstats, bestp_arr,nbest,
+ m_msg, pst, &m_msg.hist, &m_msg.pstat_void);
+ }
+ else {
+ last_stats(aa0[0], m_msg.n0,
+ qstats,nqstats, bestp_arr,nbest, m_msg, pst,
+ &m_msg.hist, &m_msg.pstat_void);
+ }
+
+ /* here is a contradiction: if pst.zsflag < 0, then m_msg.pstat_void
+ should be NULL; if it is not, then process_hist() has been called */
+ if (pst.zsflag < 0 && m_msg.pstat_void != NULL) pst.zsflag = 1;
+
+ if (m_msg.last_calc_flg) {
+ /* last_calc may need coefficients from last_stats() */
+ nbest = last_calc(aa0, aa1, maxn, bestp_arr, nbest, m_msg, &pst,
+ f_str, m_msg.pstat_void);
+ }
+
+ scale_scores(bestp_arr,nbest,m_msg.db,pst,m_msg.pstat_void);
+
+ get_param(&pst, gstring2,gstring3);
+
+#ifdef PRSS
+ /* gettitle(m_msg.lname,m_msg.ltitle,sizeof(m_msg.ltitle)); */
+ printf("%.50s - %s %d %s%s\n vs %.60s - %s shuffled sequence\n",
+ m_msg.tname, m_msg.qtitle,m_msg.n0, m_msg.sqnam,
+ (m_msg.revcomp ? " (reverse complement)" : "\0"),
+ m_msg.lname,m_msg.ltitle);
+#endif
+
+ prhist (stdout, m_msg, pst, m_msg.hist, nstats, m_msg.db, gstring2);
+
+ tscan = s_time();
+ printf (" Scan time: ");
+ ptime(stdout,tscan-tprev);
+ printf ("\n");
+#ifdef COMP_MLIB
+ ttscan += tscan-tprev;
+#endif
+
+ l3:
+ if (!m_msg.quiet) {
+ printf("Enter filename for results [%s]: ", m_msg.outfile);
+ fflush(stdout);
+ }
+
+ rline[0]='\0';
+ if (!m_msg.quiet && fgets(rline,sizeof(rline),stdin)==NULL) goto end_l;
+ if ((bp=strchr(rline,'\n'))!=NULL) *bp = '\0';
+ if (rline[0]!='\0') strncpy(m_msg.outfile,rline,sizeof(m_msg.outfile));
+ if (m_msg.outfile[0]!='\0') {
+ if ((outfd=fopen(m_msg.outfile,"w"))==NULL) {
+ fprintf(stderr," could not open %s\n",m_msg.outfile);
+ if (!m_msg.quiet) goto l3;
+ else goto l4;
+ }
+
+#ifdef PGM_DOC
+ fputs(argv_line,outfd);
+ fputc('\n',outfd);
+#endif
+ fputs(iprompt0,outfd);
+ fprintf(outfd," %s%s\n",verstr,refstr);
+
+ fprintf(outfd," %s%s, %d %s\n vs %s %s\n",
+ qlabel, (m_msg.revcomp ? "-" : "\0"), m_msg.n0,
+ m_msg.sqnam, m_msg.ltitle, lib_label);
+
+ prhist(outfd,m_msg,pst,m_msg.hist, nstats, m_msg.db, gstring2);
+ }
+
+ l4:
+ if (m_msg.markx & MX_HTML) {
+ fputs("</pre>\n<p>\n<hr>\n<p>\n",outfd);
+ }
+
+ /* code from p2_complib.c to pre-calculate -m 9 alignment info -
+ requires -q with -m 9 */
+
+ if (m_msg.quiet || m_msg.markx & MX_M9SUMM) {
+
+ /* to determine how many sequences to re-align (either for
+ do_opt() or calc_id() we need to modify m_msg.mshow to get
+ the correct number of alignments */
+
+ if (m_msg.mshow_flg != 1 && pst.zsflag >= 0) {
+ for (i=0; i<nbest && bestp_arr[i]->escore< m_msg.e_cut; i++) {}
+ m_msg.mshow = i;
+ }
+
+#ifndef PRSS
+ if (m_msg.mshow <= 0) { /* no results to display */
+ fprintf(outfd,"!! No sequences with E() < %f\n",m_msg.e_cut);
+ m_msg.nshow = 0;
+ goto end_l;
+ }
+#endif
+ }
+
+#ifdef PRSS
+ memcpy(aa1,aa1s,n1s);
+ maxn = n1s;
+ nbest = 1;
+#endif
+
+ showbest (stdout, aa0, aa1, maxn, bestp_arr, nbest, qtt.entries, &m_msg, pst,
+ m_msg.db, gstring2, f_str);
+
+ if (outfd != stdout) {
+ t_quiet = m_msg.quiet;
+ m_msg.quiet = -1; /* should guarantee 1..nbest shown */
+ showbest (outfd, aa0, aa1, maxn, bestp_arr, nbest, qtt.entries, &m_msg, pst,
+ m_msg.db, gstring2, f_str);
+ m_msg.quiet = t_quiet;
+ }
+
+ if (m_msg.nshow > 0) {
+ rline[0]='N';
+ if (!m_msg.quiet){
+ printf(" Display alignments also? (y/n) [n] "); fflush(stdout);
+ if (fgets(rline,sizeof(rline),stdin)==NULL) goto end_l;
+ }
+ else rline[0]='Y';
+
+ if (toupper((int)rline[0])=='Y') {
+ if (!m_msg.quiet) {
+ printf(" number of alignments [%d]? ",m_msg.nshow);
+ fflush(stdout);
+ if (fgets(rline,sizeof(rline),stdin)==NULL) goto end_l;
+ if (rline[0]!=0) sscanf(rline,"%d",&m_msg.nshow);
+ m_msg.ashow=m_msg.nshow;
+ }
+
+ if (m_msg.markx & (MX_AMAP+ MX_HTML + MX_M9SUMM)) {
+ fprintf(outfd,"\n>>>%s%s, %d %s vs %s library\n",
+ qlabel,(m_msg.revcomp ? "_rev":"\0"), m_msg.n0,
+ m_msg.sqnam,m_msg.lname);
+ }
+
+ if (m_msg.markx & MX_M10FORM) {
+ fprintf(outfd,"\n>>>%s%s, %d %s vs %s library\n",
+ qlabel,(m_msg.revcomp ? "-":"\0"), m_msg.n0, m_msg.sqnam,
+ m_msg.lname);
+ fprintf(outfd,"; pg_name: %s\n",argv[0]);
+ fprintf(outfd,"; pg_ver: %s\n",mp_verstr);
+ fprintf(outfd,"; pg_argv:");
+ for (i=0; i<argc; i++)
+ fprintf(outfd," %s",argv[i]);
+ fputc('\n',outfd);
+ fputs(gstring3,outfd);
+ fputs(hstring1,outfd);
+ }
+
+#ifndef PRSS
+ showalign (outfd, aa0, aa1, maxn, bestp_arr, nbest, qtt.entries,
+ m_msg, pst, gstring2, f_str);
+#else
+ if (pst.sw_flag > 0 || (!m_msg.quiet && m_msg.nshow>0)) {
+ showalign (outfd, aa0, aa1, maxn, bestp_arr, nbest, qtt.entries,
+ m_msg, pst, gstring2, f_str);
+ }
+#endif
+
+ fflush(outfd);
+ }
+ }
+
+ end_l:
+#if defined(COMP_THR) && defined(COMP_MLIB)
+ for (i=0; i<max_work_buf; i++) {
+ buf_list[i].buf_cnt=0;
+ buf_list[i].have_results=0;
+ }
+
+ num_worker_bufs = 0;
+ num_reader_bufs = max_work_buf;
+ reader_done = 0;
+ worker_buf_workp = 0;
+ worker_buf_readp = 0;
+ reader_buf_workp = 0;
+ reader_buf_readp = 0;
+
+ start_thread = 1; /* stop thread from starting again */
+#endif
+
+ /* clean up alignment encodings */
+ for (i=0; i < m_msg.nshow; i++) {
+ if (bestp_arr[i]->have_ares) {
+ free(bestp_arr[i]->a_res.res);
+ bestp_arr[i]->a_res.res = NULL;
+ bestp_arr[i]->have_ares = 0;
+ }
+ }
+
+ if (m_msg.qframe == 2) free(aa0[1]-1);
+
+ if (have_f_str) {
+ if (f_str[1]!=f_str[0]) {
+ close_work (aa0[1], m_msg.n0, &pst, &f_str[1]);
+ }
+ close_work (aa0[0], m_msg.n0, &pst, &f_str[0]);
+ have_f_str = 0;
+#ifndef COMP_THR
+ if (m_msg.qshuffle) close_work (aa0s, m_msg.n0, &pst, &qf_str);
+#endif
+ if (pst.pam_pssm) {
+ free_pam2p(pst.pam2p[0]);
+ free_pam2p(pst.pam2p[1]);
+ }
+ }
+
+ for (iln=0; iln < m_msg.nln; iln++) {
+ if (m_msg.lb_mfd[iln]!=NULL) closelib(m_msg.lb_mfd[iln]);
+ }
+
+ tddone = time(NULL);
+ tdone = s_time();
+ fflush(outfd);
+
+ if (fdata) {
+ fprintf(fdata,"/** %s **/\n",gstring2);
+ fprintf(fdata,"%3ld%-50s\n",qtt.entries-1,m_msg.qtitle);
+ fflush(fdata);
+ }
+
+#ifdef COMP_MLIB
+ ttdisp += tdone-tscan;
+
+ maxn = m_msg.max_tot;
+ m_msg.n0 =
+ QGETLIB (aa0[0], MAXTST, m_msg.qtitle, sizeof(m_msg.qtitle),
+ &qseek, &qlcont,q_file_p,&m_msg.sq0off);
+ if (m_msg.n0 <= 0) break;
+ if ((bp=strchr(m_msg.qtitle,' '))!=NULL) *bp='\0';
+ strncpy(qlabel, m_msg.qtitle,sizeof(qlabel));
+ if (bp != NULL) *bp=' ';
+ qlabel[sizeof(qlabel)-1]='\0';
+
+ if (m_msg.ann_flg) {
+ m_msg.n0 = ann_scan(aa0[0],m_msg.n0,&m_msg,m_msg.qdnaseq);
+ }
+
+ if (m_msg.term_code && m_msg.qdnaseq==SEQT_PROT &&
+ aa0[0][m_msg.n0-1]!=m_msg.term_code) {
+ aa0[0][m_msg.n0++]=m_msg.term_code;
+ aa0[0][m_msg.n0]=0;
+ }
+
+#if defined(SW_ALTIVEC) || defined(SW_SSE2)
+ /* for ALTIVEC, must pad with 15 NULL's */
+ for (id=0; id<SEQ_PAD; id++) {aa0[0][m_msg.n0+id]=0;}
+#endif
+
+#ifdef SUPERFAMNUM
+ m_msg.nqsfnum = nsfnum;
+ for (i=0; i <= nsfnum & i<10; i++) m_msg.qsfnum[i] = sfnum[i];
+ m_msg.nqsfnum_n = nsfnum_n;
+ for (i=0; i <= nsfnum_n & i<10; i++) m_msg.qsfnum_n[i] = sfnum_n[i];
+#endif
+ }
+#endif
+ if (m_msg.markx & MX_M10FORM)
+ fprintf(outfd,">>><<<\n");
+
+ tdone = s_time();
+ if ( m_msg.markx & MX_HTML) fputs("<p><pre>\n",outfd);
+ printsum(outfd, m_msg.db);
+ if ( m_msg.markx & MX_HTML) fputs("</pre>\n",outfd);
+#ifdef HTML_HEAD
+ if (m_msg.markx & MX_HTML) fprintf(outfd,"</body>\n</html>\n");
+#endif
+ if (outfd!=stdout) printsum(stdout,m_msg.db);
+
+ exit(0);
+} /* End of main program */
+
+void
+printsum(FILE *fd, struct db_str ntt)
+{
+ double db_tt;
+ char tstr1[26], tstr2[26];
+
+ strncpy(tstr1,ctime(&tdstart),sizeof(tstr1));
+ strncpy(tstr2,ctime(&tddone),sizeof(tstr1));
+ tstr1[24]=tstr2[24]='\0';
+
+ /* Print timing to output file as well */
+ fprintf(fd, "\n\n%ld residues in %ld query sequences\n", qtt.length, qtt.entries);
+ if (ntt.carry == 0)
+ fprintf(fd, "%ld residues in %ld library sequences\n", ntt.length, ntt.entries);
+ else {
+ db_tt = (double)ntt.carry*(double)LONG_MAX + (double)ntt.length;
+ fprintf(fd, "%.0f residues in %ld library sequences\n", db_tt, ntt.entries);
+ }
+
+#ifndef COMP_THR
+ fprintf(fd," Scomplib [%s]\n start: %s done: %s\n",mp_verstr,tstr1,tstr2);
+#else
+ fprintf(fd," Tcomplib [%s] (%d proc)\n start: %s done: %s\n", mp_verstr,
+ max_workers,tstr1,tstr2);
+#endif
+#ifndef COMP_MLIB
+ fprintf(fd," Scan time: ");
+ ptime(fd, tscan - tprev);
+ fprintf (fd," Display time: ");
+ ptime (fd, tdone - tscan);
+#else
+ fprintf(fd," Total Scan time: ");
+ ptime(fd, ttscan);
+ fprintf (fd," Total Display time: ");
+ ptime (fd, ttdisp);
+#endif
+ fprintf (fd,"\n");
+ fprintf (fd, "\nFunction used was %s [%s]\n", prog_func,verstr);
+}
+
+void fsigint()
+{
+ struct db_str db;
+
+ db.entries = db.length = db.carry = 0;
+ tdone = s_time();
+ tddone = time(NULL);
+
+ printf(" /*** interrupted ***/\n");
+ if (outfd!=stdout) fprintf(outfd,"/*** interrupted ***/\n");
+ fprintf(stderr,"/*** interrupted ***/\n");
+
+ printsum(stdout,db);
+ if (outfd!=stdout) printsum(outfd,db);
+
+ exit(1);
+}
+
+#ifdef COMP_THR
+void save_best(struct buf_head *cur_buf, struct mngmsg m_msg, struct pstruct pst,
+ FILE *fdata, int *qsfnum, struct hist_str *histp,
+ void **pstat_voidp
+#ifdef PRSS
+ , int *s_score, int *s_n1
+
+#endif
+ )
+{
+ double zscore;
+ int i_score;
+ struct buf_str *p_rbuf, *cur_buf_p;
+ int i, t_best, t_rbest, t_qrbest, tm_best, t_n1, sc_ix;
+ double e_score, tm_escore, t_rescore, t_qrescore;
+ int jstats;
+
+ sc_ix = pst.score_ix;
+
+ cur_buf_p = cur_buf->buf;
+
+ t_best = t_rbest = t_qrbest = -1;
+ tm_escore = t_rescore = t_qrescore = FLT_MAX;
+
+ while (cur_buf->buf_cnt--) { /* count down the number of results */
+ p_rbuf = cur_buf_p++; /* step through the results buffer */
+
+ i_score = p_rbuf->rst.score[sc_ix];
+ e_score = p_rbuf->rst.escore;
+
+ /* need to look for frame 0 if TFASTA, then save stats at frame 6 */
+ if (fdata) {
+ fprintf(fdata,
+ "%-12s %5d %6d %d %.5f %.5f %4d %4d %4d %g %d %d %8ld\n",
+ p_rbuf->libstr,
+#ifdef SUPERFAMNUM
+ sfn_cmp(qsfnum,p_rbuf->sfnum),
+#else
+ 0,
+#endif
+ p_rbuf->n1,p_rbuf->frame,p_rbuf->rst.comp,p_rbuf->rst.H,
+ p_rbuf->rst.score[0],p_rbuf->rst.score[1],p_rbuf->rst.score[2],
+ p_rbuf->rst.escore, p_rbuf->rst.segnum, p_rbuf->rst.seglen, p_rbuf->lseek);
+ }
+
+#ifdef PRSS
+ if (p_rbuf->lseek==0) {
+ s_score[0] = p_rbuf->rst.score[0];
+ s_score[1] = p_rbuf->rst.score[1];
+ s_score[2] = p_rbuf->rst.score[2];
+ *s_n1 = p_rbuf->n1;
+
+ bestp = bestp_arr[nbest++];
+ bestp->score[0] = s_score[0];
+ bestp->score[1] = s_score[1];
+ bestp->score[2] = s_score[2];
+ bestp->n1 = *s_n1;
+ bestp->escore = p_rbuf->rst.escore;
+ bestp->segnum = p_rbuf->rst.segnum;
+ bestp->seglen = p_rbuf->rst.seglen;
+ bestp->zscore = zscore;
+ bestp->lseek = p_rbuf->lseek;
+ bestp->m_file_p = p_rbuf->m_file_p;
+ memcpy(bestp->libstr,p_rbuf->libstr,MAX_UID);
+ bestp->n1tot_p = p_rbuf->n1tot_p;
+ bestp->frame = p_rbuf->frame;
+
+ continue;
+ }
+#endif
+
+ t_n1 = p_rbuf->n1;
+ if (i_score > t_best) tm_best = t_best = i_score;
+ if (e_score < tm_escore) tm_escore = e_score;
+
+ if (m_msg.qshuffle) {
+ if (p_rbuf->qr_score > t_qrbest)
+ t_qrbest = p_rbuf->qr_score;
+ if (p_rbuf->qr_escore < t_qrescore)
+ t_qrescore = p_rbuf->qr_escore;
+
+ if (p_rbuf->frame == m_msg.nitt1 && nqstats < m_msg.shuff_max) {
+ qstats[nqstats].n1 = p_rbuf->n1; /* save the best score */
+ qstats[nqstats].comp = p_rbuf->rst.comp;
+ qstats[nqstats].H = p_rbuf->rst.H;
+ qstats[nqstats].escore = t_qrescore;
+ qstats[nqstats++].score = t_qrbest;
+ t_qrbest = -1; /* reset t_qrbest, t_qrescore */
+ t_qrescore = FLT_MAX;
+ }
+ }
+
+ if (pst.zsflag >= 10 && p_rbuf->r_score > t_rbest) {
+ t_rbest = p_rbuf->r_score;
+ t_rescore = p_rbuf->r_escore;
+ }
+
+ /* statistics done for best score of set */
+
+
+ if (p_rbuf->frame == m_msg.nitt1) {
+ if (nstats < MAXSTATS ) {
+ stats[nstats].n1 = t_n1;
+ stats[nstats].comp = p_rbuf->rst.comp;
+ stats[nstats].H = p_rbuf->rst.H;
+ if (pst.zsflag >= 10) {
+ tm_best = t_rbest;
+ tm_escore = t_rescore;
+ t_rbest = -1;
+ t_rescore = FLT_MAX;
+ }
+ stats[nstats].escore = tm_escore;
+ stats[nstats++].score = tm_best;
+ t_best = -1;
+ tm_escore = FLT_MAX;
+ }
+ else if (pst.zsflag > 0) {
+ if (!stats_done) {
+ pst.zsflag_f = process_hist(stats,nstats,m_msg,pst,
+ histp, pstat_voidp,0);
+ kstats = nstats;
+ stats_done = 1;
+ for (i=0; i<MAXBEST; i++) {
+ bestp_arr[i]->zscore =
+ (*find_zp)(bestp_arr[i]->score[pst.score_ix],
+ bestp_arr[i]->escore, bestp_arr[i]->n1,
+ bestp_arr[i]->comp, *pstat_voidp);
+ }
+ }
+#ifdef SAMP_STATS
+ else {
+ if (!m_msg.escore_flg) {
+ jstats = nrand(++kstats);
+ if (jstats < MAXSTATS) {
+ stats[jstats].n1 = t_n1;
+ stats[jstats].comp = p_rbuf->rst.comp;
+ stats[jstats].H = p_rbuf->rst.H;
+ if (pst.zsflag >= 10) {
+ tm_best = t_rbest;
+ }
+ stats[jstats].score = tm_best;
+ }
+ }
+ }
+#endif
+ }
+ }
+
+ /* best saved for every score */
+ if (stats_done) {
+
+ zscore=(*find_zp)(i_score, e_score, p_rbuf->n1,(double)p_rbuf->rst.comp,
+ *pstat_voidp);
+
+ if (p_rbuf->frame == m_msg.nitt1) {
+ addhistz((*find_zp)(t_best, tm_escore, p_rbuf->n1, (double) p_rbuf->rst.comp,
+ *pstat_voidp), histp);
+ t_best = t_rbest = -1;
+ tm_escore = t_rescore = FLT_MAX;
+ }
+ }
+ else zscore = (double) i_score;
+
+#ifndef PRSS
+ if (zscore > zbestcut) {
+ if (nbest >= MAXBEST) {
+ bestfull = nbest-MAXBEST/4;
+ selectbestz(bestp_arr,bestfull-1,nbest);
+ zbestcut = bestp_arr[bestfull-1]->zscore;
+ nbest = bestfull;
+ }
+ bestp = bestp_arr[nbest++];
+ bestp->score[0] = p_rbuf->rst.score[0];
+ bestp->score[1] = p_rbuf->rst.score[1];
+ bestp->score[2] = p_rbuf->rst.score[2];
+ bestp->comp = (double) p_rbuf->rst.comp;
+ bestp->H = (double) p_rbuf->rst.H;
+ bestp->escore = p_rbuf->rst.escore;
+ bestp->segnum = p_rbuf->rst.segnum;
+ bestp->seglen = p_rbuf->rst.seglen;
+ bestp->zscore = zscore;
+ bestp->lseek = p_rbuf->lseek;
+ memcpy(bestp->libstr,p_rbuf->libstr,MAX_UID);
+ bestp->cont = p_rbuf->cont; /* not cont+1 because incremented already */
+ bestp->m_file_p = p_rbuf->m_file_p;
+ bestp->n1 = p_rbuf->n1;
+ bestp->n1tot_p = p_rbuf->n1tot_p;
+ bestp->frame = p_rbuf->frame;
+ bestp->nsfnum = p_rbuf->nsfnum;
+#ifdef SUPERFAMNUM
+ if ((bestp->sfnum[0] = p_rbuf->sfnum[0])>0 &&
+ (bestp->sfnum[1] = p_rbuf->sfnum[1])>0 &&
+ (bestp->sfnum[2] = p_rbuf->sfnum[2])>0 &&
+ (bestp->sfnum[3] = p_rbuf->sfnum[3])>0 &&
+ (bestp->sfnum[4] = p_rbuf->sfnum[4])>0 &&
+ (bestp->sfnum[5] = p_rbuf->sfnum[5])>0 &&
+ (bestp->sfnum[6] = p_rbuf->sfnum[6])>0 &&
+ (bestp->sfnum[7] = p_rbuf->sfnum[7])>0 &&
+ (bestp->sfnum[8] = p_rbuf->sfnum[8])>0 &&
+ (bestp->sfnum[9] = p_rbuf->sfnum[9])>0) ;
+#endif
+ }
+#endif
+ }
+}
+#endif
git://source.jalview.org / jabaws.git / blobdiff