hmmcalibrate.c

hmmer源程序

       * The 9999. is an arbitrary high number that means we won't trim
       * outliers on the right.
       */
      if (! ExtremeValueFitHistogram(hist, TRUE, 9999.))
	Die("fit failed; -n may be set too small?\n");
      
      mu[nhmm]     = hist->param[EVD_MU];
      lambda[nhmm] = hist->param[EVD_LAMBDA];
      nhmm++;
      if (nhmm % 100 == 0) {
	mu     = ReallocOrDie(mu,     sizeof(float) * (nhmm+mu_lumpsize));
	lambda = ReallocOrDie(lambda, sizeof(float) * (nhmm+mu_lumpsize));
      }      

      /* Output
       */
      printf("HMM    : %s\n",   hmm->name);
      printf("mu     : %12f\n", hist->param[EVD_MU]);
      printf("lambda : %12f\n", hist->param[EVD_LAMBDA]);
      printf("max    : %12f\n", max);
      printf("//\n");

      if (hfp != NULL) 
	{
	  fprintf(hfp, "HMM: %s\n", hmm->name);
	  PrintASCIIHistogram(hfp, hist);
	  fprintf(hfp, "//\n");
	}

      FreeHistogram(hist);
    }
  SQD_DPRINTF1(("Main body believes it has calibrations for %d HMMs\n", nhmm));

  /*****************************************************************
   * Rewind the HMM file for a second pass.
   * Write a temporary HMM file with new mu, lambda values in it
   *****************************************************************/

  HMMFileRewind(hmmfp);
  if (FileExists(tmpfile))
    Die("Ouch. Temporary file %s appeared during the run.", tmpfile);
  if ((outfp = fopen(tmpfile, mode)) == NULL)
    Die("Ouch. Temporary file %s couldn't be opened for writing.", tmpfile); 
  
  for (idx = 0; idx < nhmm; idx++)
    {
      /* Sanity checks 
       */
      if (!HMMFileRead(hmmfp, &hmm))
	Die("Ran out of HMMs too early in pass 2");
      if (hmm == NULL) 
	Die("HMM file %s was corrupted? Parse failed in pass 2", hmmfile);

      /* Put results in HMM
       */
      hmm->mu     = mu[idx];
      hmm->lambda = lambda[idx];
      hmm->flags |= PLAN7_STATS;
      Plan7ComlogAppend(hmm, argc, argv);

      /* Save HMM to tmpfile
       */
      if (hmmfp->is_binary) WriteBinHMM(outfp, hmm);
      else                  WriteAscHMM(outfp, hmm); 

      FreePlan7(hmm);
    }
  
  /*****************************************************************
   * Now, carefully remove original file and replace it
   * with the tmpfile. Note the protection from signals;
   * we wouldn't want a user to ctrl-C just as we've deleted
   * their HMM file but before the new one is moved.
   *****************************************************************/

  HMMFileClose(hmmfp);
  if (fclose(outfp)   != 0) PANIC;

  if (sigemptyset(&blocksigs) != 0) PANIC;
  if (sigaddset(&blocksigs, SIGINT) != 0) PANIC;
  if (sigprocmask(SIG_BLOCK, &blocksigs, NULL) != 0)   PANIC;
  if (remove(hmmfile) != 0)                            PANIC;
  if (rename(tmpfile, hmmfile) != 0)                   PANIC;
  if (sigprocmask(SIG_UNBLOCK, &blocksigs, NULL) != 0) PANIC;

  /***********************************************
   * Exit
   ***********************************************/

  StopwatchStop(&stopwatch);
  if (do_pvm > 0) {
    printf("PVM processors used: %d\n", pvm_nslaves);
    StopwatchInclude(&stopwatch, &extrawatch);
  }
#ifdef PTHREAD_TIMES_HACK
  else if (num_threads > 0) StopwatchInclude(&stopwatch, &extrawatch);
#endif

  /*  StopwatchDisplay(stdout, "CPU Time: ", &stopwatch); */

  free(mu);
  free(lambda);
  free(tmpfile);
  if (hfp != NULL) fclose(hfp);
  SqdClean();
  return 0;
}

/* Function: main_loop_serial()
 * Date:     SRE, Tue Aug 18 16:18:28 1998 [St. Louis]
 *
 * Purpose:  Given an HMM and parameters for synthesizing random
 *           sequences; return a histogram of scores.
 *           (Serial version)  
 *
 * Args:     hmm      - an HMM to calibrate.
 *           seed     - random number seed
 *           nsample  - number of seqs to synthesize
 *           lenmean  - mean length of random sequence
 *           lensd    - std dev of random seq length
 *           fixedlen - if nonzero, override lenmean, always this len
 *           ret_hist - RETURN: the score histogram 
 *           ret_max  - RETURN: highest score seen in simulation
 *
 * Returns:  (void)
 *           hist is alloc'ed here, and must be free'd by caller.
 */
static void
main_loop_serial(struct plan7_s *hmm, int seed, int nsample, 
		 float lenmean, float lensd, int fixedlen,
		 struct histogram_s **ret_hist, float *ret_max)
{
  struct histogram_s *hist;
  float  randomseq[MAXABET];
  float  p1;
  float  max;
  char  *seq;
  char  *dsq;
  float  score;
  int    sqlen;
  int    idx;
  
  /* Initialize.
   * We assume we've already set the alphabet (safe, because
   * HMM input sets the alphabet).
   */
  sre_srandom(seed);
  P7Logoddsify(hmm, TRUE);
  P7DefaultNullModel(randomseq, &p1);
  hist = AllocHistogram(-200, 200, 100);
  max = -FLT_MAX;

  for (idx = 0; idx < nsample; idx++)
    {
				/* choose length of random sequence */
      if (fixedlen) sqlen = fixedlen;
      else do sqlen = (int) Gaussrandom(lenmean, lensd); while (sqlen < 1);
				/* generate it */
      seq = RandomSequence(Alphabet, randomseq, Alphabet_size, sqlen);
      dsq = DigitizeSequence(seq, sqlen);

      if (P7ViterbiSize(sqlen, hmm->M) <= RAMLIMIT)
	score = P7Viterbi(dsq, sqlen, hmm, NULL);
      else
	score = P7SmallViterbi(dsq, sqlen, hmm, NULL);

      AddToHistogram(hist, score);
      if (score > max) max = score;

      free(dsq); 
      free(seq);
    }

  *ret_hist   = hist;
  *ret_max    = max;
  return;
}


#ifdef HMMER_THREADS
/* Function: main_loop_threaded()
 * Date:     SRE, Wed Dec  1 12:43:09 1999 [St. Louis]
 *
 * Purpose:  Given an HMM and parameters for synthesizing random
 *           sequences; return a histogram of scores.
 *           (Threaded version.)  
 *
 * Args:     hmm      - an HMM to calibrate.
 *           seed     - random number seed
 *           nsample  - number of seqs to synthesize
 *           lenmean  - mean length of random sequence
 *           lensd    - std dev of random seq length
 *           fixedlen - if nonzero, override lenmean, always this len
 *           nthreads - number of threads to start
 *           ret_hist - RETURN: the score histogram 
 *           ret_max  - RETURN: highest score seen in simulation
 *           twatch   - RETURN: accumulation of thread times
 *
 * Returns:  (void)
 *           hist is alloc'ed here, and must be free'd by caller.
 */
static void
main_loop_threaded(struct plan7_s *hmm, int seed, int nsample, 
		   float lenmean, float lensd, int fixedlen,
		   int nthreads,
		   struct histogram_s **ret_hist, float *ret_max,
		   Stopwatch_t *twatch)
{
  struct histogram_s *hist;
  float  randomseq[MAXABET];
  float  p1;
  struct workpool_s *wpool;     /* pool of worker threads  */
  
  /* Initialize.
   * We assume we've already set the alphabet (safe, because
   * HMM input sets the alphabet).
   */
  sre_srandom(seed);
  P7Logoddsify(hmm, TRUE);
  P7DefaultNullModel(randomseq, &p1);
  hist = AllocHistogram(-200, 200, 100);

  wpool = workpool_start(hmm, lenmean, lensd, fixedlen, randomseq, nsample,
			 hist, nthreads);
  workpool_stop(wpool);

  *ret_hist = hist;
  *ret_max  = wpool->max_score;
  StopwatchInclude(twatch, &(wpool->watch));

  workpool_free(wpool);
  return;
}

/*****************************************************************
 * POSIX threads implementation.
 * API:
 *      workpool_start()   (makes a workpool_s structure. Starts calculations.)
 *      workpool_stop()    (waits for threads to finish.)
 *      [process histogram]
 *      workpool_free()    (destroys the structure)
 *      
 * Threads:
 *      worker_thread()    (the actual parallelized worker thread).
 *****************************************************************/

/* Function: workpool_start()
 * Date:     SRE, Thu Jul 16 11:09:05 1998 [St. Louis]
 *
 * Purpose:  Initialize a workpool_s structure, and return it.
 *
 * Args:     hmm      - the HMM to calibrate
 *           fixedlen - 0, or a fixed length for seqs (bypass of Gaussian)
 *           lenmean  - mean sequence length 
 *           lensd    - std. dev. for sequence length
 *           randomseq- i.i.d. frequencies for residues, 0..Alphabet_size-1
 *           nsample  - how many seqs to calibrate on
 *           hist     - histogram structure for storing results
 *           num_threads - how many processors to run on
 *
 * Returns:  ptr to struct workpool_s.
 *           Caller must wait for threads to finish with workpool_stop(),
 *           then free the structure with workpool_free().
 */
static struct workpool_s *
workpool_start(struct plan7_s *hmm, float lenmean, float lensd, int fixedlen,
	       float *randomseq, int nsample, struct histogram_s *hist, 
	       int num_threads)
{
  struct workpool_s *wpool;
  pthread_attr_t    attr;
  int i;
  int rtn;

  wpool         = MallocOrDie(sizeof(struct workpool_s));
  wpool->thread = MallocOrDie(num_threads * sizeof(pthread_t));
  wpool->hmm        = hmm;
  wpool->fixedlen   = fixedlen;
  wpool->lenmean    = lenmean;
  wpool->lensd      = lensd;
  wpool->randomseq  = randomseq;
  wpool->nsample    = nsample;
  
  wpool->nseq       = 0;
  wpool->hist       = hist;
  wpool->max_score  = -FLT_MAX;
  wpool->num_threads= num_threads;

  StopwatchZero(&(wpool->watch));
  
  if ((rtn = pthread_mutex_init(&(wpool->input_lock), NULL)) != 0)
    Die("pthread_mutex_init FAILED; %s\n", strerror(rtn));
  if ((rtn = pthread_mutex_init(&(wpool->output_lock), NULL)) != 0)
    Die("pthread_mutex_init FAILED; %s\n", strerror(rtn));

  /* Create slave threads.
   * Note the crazy machinations we have to go through to achieve concurrency.
   * You'd think that POSIX threads were portable... ha.
   * On IRIX 6.5, system scope threads are only available to root, or if
   *   /etc/capability has been configured specially, so to avoid strange
   *   permissions errors we can't set PTHREAD_SCOPE_SYSTEM for IRIX.
   * On IRIX pre-6.5, we can't get good concurrency, period. As of 6.5,
   *   SGI provides the nonportable pthread_setconcurrency() call.
   * On FreeBSD (3.0 snapshots), the pthread_attr_setscope() call isn't
   *   even provided, apparently on grounds of "if it doesn't do anything,
   *   why provide it?" Hello? POSIX compliance, perhaps?
   * On Sun Solaris, we need to set system scope to achieve concurrency.
   * Linux and DEC Digital UNIX seem to work fine in either process scope
   *   or system scope, without a pthread_setconcurrency call.
   */
  pthread_attr_init(&attr);
#ifndef __sgi
#ifdef HAVE_PTHREAD_ATTR_SETSCOPE
  pthread_attr_setscope(&attr, PTHREAD_SCOPE_SYSTEM);
#endif
#endif
#ifdef HAVE_PTHREAD_SETCONCURRENCY
  pthread_setconcurrency(num_threads+1);
#endif
  for (i = 0; i < num_threads; i++)