hmmcalibrate.c

hmmer源程序

    if ((rtn = pthread_create(&(wpool->thread[i]), &attr,
			      worker_thread , (void *) wpool)) != 0)
      Die("Failed to create thread %d; return code %d\n", i, rtn);

  pthread_attr_destroy(&attr);

  return wpool;
}

/* Function: workpool_stop()
 * Date:     SRE, Thu Jul 16 11:20:16 1998 [St. Louis]
 *
 * Purpose:  Waits for threads in a workpool to finish.
 *
 * Args:     wpool -- ptr to the workpool structure
 *
 * Returns:  (void)
 */
static void
workpool_stop(struct workpool_s *wpool)
{
  int i;
				/* wait for threads to stop */
  for (i = 0; i < wpool->num_threads; i++)
    if (pthread_join(wpool->thread[i],NULL) != 0)
      Die("pthread_join failed");
  return;
}

/* Function: workpool_free()
 * Date:     SRE, Thu Jul 16 11:26:27 1998 [St. Louis]
 *
 * Purpose:  Free a workpool_s structure, after the threads
 *           have finished.
 *
 * Args:     wpool -- ptr to the workpool.
 *
 * Returns:  (void)
 */
static void
workpool_free(struct workpool_s *wpool)
{
  free(wpool->thread);
  free(wpool);
  return;
}

/* Function: worker_thread()
 * Date:     SRE, Thu Jul 16 10:41:02 1998 [St. Louis]
 *
 * Purpose:  The procedure executed by the worker threads.
 *
 * Args:     ptr  - (void *) that is recast to a pointer to
 *                  the workpool.
 *
 * Returns:  (void *)
 */
void *
worker_thread(void *ptr)
{
  struct plan7_s    *hmm;
  struct workpool_s *wpool;
  char       *seq;
  char       *dsq;
  int         len;
  float       sc;
  int         rtn;
  Stopwatch_t thread_watch;

  StopwatchStart(&thread_watch);
  wpool = (struct workpool_s *) ptr;
  hmm   = wpool->hmm;
  for (;;)
    {
      /* 1. Synthesize a random sequence. 
       *    The input sequence number is a shared resource,
       *    and sre_random() isn't thread-safe, so protect
       *    the whole section with mutex.
       */
				/* acquire a lock */
      if ((rtn = pthread_mutex_lock(&(wpool->input_lock))) != 0)
	Die("pthread_mutex_lock failure: %s\n", strerror(rtn));
				/* generate a sequence */
      wpool->nseq++;
      if (wpool->nseq > wpool->nsample) 
	{ /* we're done; release input lock, break loop */
	  if ((rtn = pthread_mutex_unlock(&(wpool->input_lock))) != 0)
	    Die("pthread_mutex_unlock failure: %s\n", strerror(rtn));
	  break;
	}
      if (wpool->fixedlen) len = wpool->fixedlen;
      else do len = (int) Gaussrandom(wpool->lenmean, wpool->lensd); while (len < 1);
      seq = RandomSequence(Alphabet, wpool->randomseq, Alphabet_size, len);

				/* release the lock */
      if ((rtn = pthread_mutex_unlock(&(wpool->input_lock))) != 0)
	Die("pthread_mutex_unlock failure: %s\n", strerror(rtn));

      /* 2. Score the sequence against the model.
       */
      dsq = DigitizeSequence(seq, len);
      
      if (P7ViterbiSize(len, hmm->M) <= RAMLIMIT)
	sc = P7Viterbi(dsq, len, hmm, NULL);
      else
	sc = P7SmallViterbi(dsq, len, hmm, NULL);
      free(dsq); 
      free(seq);
      
      /* 3. Save the output; hist and max_score are shared,
       *    so protect this section with the output mutex.
       */
				/* acquire lock on the output queue */
      if ((rtn = pthread_mutex_lock(&(wpool->output_lock))) != 0)
	Die("pthread_mutex_lock failure: %s\n", strerror(rtn));
				/* save output */
      AddToHistogram(wpool->hist, sc);
      if (sc > wpool->max_score) wpool->max_score = sc;
    				/* release our lock */
      if ((rtn = pthread_mutex_unlock(&(wpool->output_lock))) != 0)
	Die("pthread_mutex_unlock failure: %s\n", strerror(rtn));
    }

  StopwatchStop(&thread_watch);
				/* acquire lock on the output queue */
  if ((rtn = pthread_mutex_lock(&(wpool->output_lock))) != 0)
    Die("pthread_mutex_lock failure: %s\n", strerror(rtn));
				/* accumulate cpu time into main stopwatch */
  StopwatchInclude(&(wpool->watch), &thread_watch);
    				/* release our lock */
  if ((rtn = pthread_mutex_unlock(&(wpool->output_lock))) != 0)
    Die("pthread_mutex_unlock failure: %s\n", strerror(rtn));

  pthread_exit(NULL);
  return NULL; /* solely to silence compiler warnings */
}
#endif /* HMMER_THREADS */



#ifdef HMMER_PVM
/* Function: main_loop_pvm()
 * Date:     SRE, Wed Aug 19 13:59:54 1998 [St. Louis]
 *
 * Purpose:  Given an HMM and parameters for synthesizing random
 *           sequences; return a histogram of scores.
 *           (PVM version)  
 *
 * Args:     hmm     - an HMM to calibrate.
 *           seed    - random number seed
 *           nsample - number of seqs to synthesize
 *           lumpsize- # of seqs per slave exchange; controls granularity
 *           lenmean - mean length of random sequence
 *           lensd   - std dev of random seq length
 *           fixedlen- if nonzero, override lenmean, always this len
 *           hist       - RETURN: the score histogram 
 *           ret_max    - RETURN: highest score seen in simulation
 *           extrawatch - RETURN: total CPU time spend in slaves.
 *           ret_nslaves- RETURN: number of PVM slaves run.
 *
 * Returns:  (void)
 *           hist is alloc'ed here, and must be free'd by caller.
 */
static void
main_loop_pvm(struct plan7_s *hmm, int seed, int nsample, int lumpsize,
	      float lenmean, float lensd, int fixedlen,
	      struct histogram_s **ret_hist, float *ret_max, 
	      Stopwatch_t *extrawatch, int *ret_nslaves)
{
  struct histogram_s *hist;
  int                 master_tid;
  int                *slave_tid;
  int                 nslaves;
  int                 nsent;	/* # of seqs we've asked for so far       */
  int                 ndone;	/* # of seqs we've got results for so far */
  int		      packet;	/* # of seqs to have a slave do           */
  float               max;
  int                 slaveidx;	/* id of a slave */
  float              *sc;        /* scores returned by a slave */
  Stopwatch_t         slavewatch;
  int                 i;
  
  StopwatchZero(extrawatch);
  hist = AllocHistogram(-200, 200, 100);
  max  = -FLT_MAX;

  /* Initialize PVM
   */
  if ((master_tid = pvm_mytid()) < 0)
    Die("pvmd not responding -- do you have PVM running?");
#if DEBUGLEVEL >= 1
  pvm_catchout(stderr);		/* catch output for debugging */
#endif
  PVMSpawnSlaves("hmmcalibrate-pvm", &slave_tid, &nslaves);

  /* Initialize the slaves
   */
  pvm_initsend(PvmDataDefault);
  pvm_pkfloat(&lenmean,       1, 1);
  pvm_pkfloat(&lensd,         1, 1);
  pvm_pkint(  &fixedlen,      1, 1);
  pvm_pkint(  &Alphabet_type, 1, 1);
  pvm_pkint(  &seed,          1, 1);
  if (! PVMPackHMM(hmm)) Die("Failed to pack the HMM");
  pvm_mcast(slave_tid, nslaves, HMMPVM_INIT);
  SQD_DPRINTF1(("Initialized %d slaves\n", nslaves));

  /* Confirm slaves' OK status.
   */
  PVMConfirmSlaves(slave_tid, nslaves);
  SQD_DPRINTF1(("Slaves confirm that they're ok...\n"));
 
  /* Load the slaves
   */
  nsent = ndone = 0;
  for (slaveidx = 0; slaveidx < nslaves; slaveidx++)
    {
      packet    = (nsample - nsent > lumpsize ? lumpsize : nsample - nsent);

      pvm_initsend(PvmDataDefault);
      pvm_pkint(&packet,    1, 1);
      pvm_pkint(&slaveidx,  1, 1);
      pvm_send(slave_tid[slaveidx], HMMPVM_WORK);
      nsent += packet;
    }
  SQD_DPRINTF1(("Loaded %d slaves\n", nslaves));

  /* Receive/send loop
   */
  sc = MallocOrDie(sizeof(float) * lumpsize);
  while (nsent < nsample)
    {
				/* integrity check of slaves */
      PVMCheckSlaves(slave_tid, nslaves);

				/* receive results */
      SQD_DPRINTF2(("Waiting for results...\n"));
      pvm_recv(-1, HMMPVM_RESULTS);
      pvm_upkint(&slaveidx,   1, 1);
      pvm_upkint(&packet,     1, 1);
      pvm_upkfloat(sc,   packet, 1);
      SQD_DPRINTF2(("Got results.\n"));
      ndone += packet;

				/* store results */
      for (i = 0; i < packet; i++) {
	AddToHistogram(hist, sc[i]);
	if (sc[i] > max) max = sc[i];
      }
				/* send new work */
      packet    = (nsample - nsent > lumpsize ? lumpsize : nsample - nsent);

      pvm_initsend(PvmDataDefault);
      pvm_pkint(&packet,    1, 1);
      pvm_pkint(&slaveidx,  1, 1);
      pvm_send(slave_tid[slaveidx], HMMPVM_WORK);
      SQD_DPRINTF2(("Told slave %d to do %d more seqs.\n", slaveidx, packet));
      nsent += packet;
    }
      
  /* Wait for the last output to come in.
   */
  while (ndone < nsample)
    {
				/* integrity check of slaves */
      PVMCheckSlaves(slave_tid, nslaves);

				/* receive results */
      SQD_DPRINTF1(("Waiting for final results...\n"));
      pvm_recv(-1, HMMPVM_RESULTS);
      pvm_upkint(&slaveidx, 1, 1);
      pvm_upkint(&packet,   1, 1);
      pvm_upkfloat(sc, packet, 1);
      SQD_DPRINTF2(("Got some final results.\n"));
      ndone += packet;
				/* store results */
      for (i = 0; i < packet; i++) {
	AddToHistogram(hist, sc[i]);
	if (sc[i] > max) max = sc[i];
      }
    }

  /* Shut down the slaves: send -1,-1,-1.
   */
  pvm_initsend(PvmDataDefault);
  packet = -1;
  pvm_pkint(&packet, 1, 1);
  pvm_pkint(&packet, 1, 1);
  pvm_pkint(&packet, 1, 1);
  pvm_mcast(slave_tid, nslaves, HMMPVM_WORK);

  /* Collect stopwatch results; quit the VM; return.
   */
  for (i = 0; i < nslaves; i++)
    {
      pvm_recv(-1, HMMPVM_RESULTS);
      pvm_upkint(&slaveidx, 1, 1);
      StopwatchPVMUnpack(&slavewatch);

      SQD_DPRINTF1(("Slave %d finished; says it used %.2f cpu, %.2f sys\n",
		    slaveidx, slavewatch.user, slavewatch.sys));

      StopwatchInclude(extrawatch, &slavewatch);
    }

  free(slave_tid);
  free(sc);
  pvm_exit();
  *ret_hist    = hist;
  *ret_max     = max;
  *ret_nslaves = nslaves;
  return;
}
#endif /* HMMER_PVM */