message.c

来自「EM算法的改进」· C语言 代码 · 共 531 行 · 第 1/2 页

/*
 * $Id: message.c 1339 2006-09-21 19:46:28Z tbailey $
 * 
 * $Log$
 * Revision 1.1  2005/07/29 17:23:03  nadya
 * Initial revision
 *
 */

/***********************************************************************
*								       *
*	MEME++							       *
*	Copyright 1996, The Regents of the University of California    *
*	Author: Tim Bailey & Bill Grundy			        *
*								       *
***********************************************************************/

#ifdef PARALLEL
 
#include <macros.h>
#include <mpi.h>
#include <mp.h>
#include <message.h>

/**********************************************************************
 *
 * void balance_loop1
 *
 * This function does load-balancing on the search for starting
 * points. It divides the dataset by the number of processors at the
 * character level. The work assigned to each processor is stored in a
 * global struct.
 *
 * IN: SAMPLE **samples             the dataset
 *     int n_samples                size of the dataset
 * OUT: SEQ_PARAMS *start_n_end     (global) start and end points for this node
 *
 * This function is called in 'init.c'.
 **********************************************************************/
void balance_loop1(SAMPLE **samples, int n_samples)
{
  int n_residues;
  float residues_per_node, target;
  int iseq, inode, seq_start, seq_end;

  /* Allocate storage for the sequence parameters. */
  start_n_end = (SEQ_PARAMS *)mymalloc(sizeof(SEQ_PARAMS));

  /* First count how many residues we have. */
  for (iseq = 0, n_residues = 0; iseq < n_samples; iseq++)
    n_residues += samples[iseq]->length;
  
  /* Calculate the number of residues per node. */
  residues_per_node = MAX(1, n_residues / mpNodes());
  
  /*
  fprintf(stderr, "n_residues=%d residues_per_node=%g\n", n_residues,
    residues_per_node);
  */
  
  /* Set the starting points for node 0. */
  if (mpMyID() == 0) {
    start_n_end->start_seq = 0;
    start_n_end->start_off = 0;
  }
  
  /* Calculate per-node starting and ending points. */
  for (iseq = 0, seq_start = 0, seq_end = samples[0]->length - 1,
       inode = 0, target = residues_per_node;
       (target < n_residues) || (iseq < n_samples);
       /* Incrementing is done in the loop */) {
    
    /* See whether we've moved past the target. */      
    if (seq_end >= target) {
      if (mpMyID() == inode) {
	start_n_end->end_seq = iseq;
	start_n_end->end_off = (int)target - seq_start - 1;
      }
      if (mpMyID() == inode + 1) {
	start_n_end->start_seq = iseq;
	start_n_end->start_off = (int)target - seq_start;
      }

      /*
      fprintf(stderr, "end %d %d\nnode %d start %d %d ",
	iseq, (int)target-seq_start-1, inode+1, iseq, (int)target-seq_start);
      */

      /* Move to a new node. */
      inode++;
      /* Calculate the target for this node. */
      target = residues_per_node * (inode+1);

    } else {
      
      /* Move to a new sequence. */
      if (++iseq < n_samples) {
	/* Calculate the starting and ending points of this sequence. */
	seq_start = seq_end;
	seq_end += samples[iseq]->length; /* 1 more than the index. */
      }
    }
  }
  
  /* Set the ending points for the last node. */
  if (mpMyID() == mpNodes() - 1) {
    start_n_end->end_seq = n_samples - 1;
    start_n_end->end_off = samples[n_samples - 1]->length - 1;
  }
  
  /* Search for (and correct) fencepost problems. */
  if (start_n_end->start_off >= samples[start_n_end->start_seq]->length) {
    start_n_end->start_seq++;
    start_n_end->start_off = 0;
  } else if (start_n_end->start_off <= -1) {
    start_n_end->start_seq--;
    start_n_end->start_off = samples[start_n_end->start_seq]->length - 1;
  }
  
  if (start_n_end->end_off >= samples[start_n_end->end_seq]->length) {
    start_n_end->end_seq++;
    start_n_end->end_off = 0;
  } else if (start_n_end->end_off <= -1) {
    start_n_end->end_seq--;
    start_n_end->end_off = samples[start_n_end->end_seq]->length - 1;
  }
  
#ifdef DEBUG
  fprintf(stderr, "%d: (%d,%d) -> (%d,%d)\n", mpMyID(), 
	  start_n_end->start_seq, start_n_end->start_off,
	  start_n_end->end_seq, start_n_end->end_off);
#endif /* DEBUG */
  
  /* Do some error checking. */
  /*
     if (inode != mpNodes()) {
       fprintf(stderr, "Warning! inode(%d) != mpNodes(%d)\n",
	     inode, mpNodes());
     }
  */
  if (iseq != n_samples)
    fprintf(stderr, "Warning! iseq(%d) != n_samples(%d)\n",
	    iseq, n_samples);
}

/**********************************************************************
 *
 * void store_consensus
 *
 * Because the reduction across models does not include the
 * candidates, we must transfer the human-readable consensus string
 * into the model so that it can be printed later.
 *
 * This function is called by 'meme.c'.
 **********************************************************************/
void store_consensus(MODEL *model, CANDIDATE *candidates)
{
    int w = model->w;				/* width of motif */
    S_POINT *s_point = candidates[w].s_point; 	/* starting point for model */

    /* Make sure we have at least one candidate. */
    if (s_point == NULL) {
      /* Guarantee that this model won't get chosen. */
      model->logev = log(BIG);
    } else {
      /*fprintf(stderr, "%d: BEFORE consensus=%s\n", mpMyID(),
	s_point->cons0);*/
      strcpy(model->cons0, s_point->cons0);
    }
} /* store_consensus */

/**********************************************************************
 *
 * void save_theta_ptrs
 *
 * When we broadcast the winning model, the array of pointers in that
 * model gets sent to everyone. Since those pointers are nonsense to
 * other nodes, those nodes have to store their own pointers and then
 * restore them after the broadcast.
 *
 * The maxima array is also saved and restored here.
 *
 **********************************************************************/

void save_theta_ptrs(MODEL *model, int save_or_restore)
{
  static THETA saved_theta;
  static THETA saved_logtheta;
  static THETA saved_obs;
  static P_PROB saved_maxima;

  if (save_or_restore == 1) { 		/* Save the theta matrix pointer */
    saved_theta = model->theta;
    saved_logtheta = model->logtheta;
    saved_obs = model->obs;
    saved_maxima = model->maxima;
  } else {
    model->theta = saved_theta;
    model->logtheta = saved_logtheta;
    model->obs = saved_obs;
    model->maxima = saved_maxima;
    Resize(model->maxima, model->nsites_dis, p_prob);
  }
} /* save_theta_ptrs */

/**********************************************************************
 *
 * void theta_packer
 *
 * This function packs (or unpacks) the 2-D theta matrix into a 1-D
 * array so that it can be sent in a single broadcast.
 * 
 **********************************************************************/
void theta_packer(
   MODEL *model,
   double *theta_matrix,
   double *obs_matrix,
   int pack_or_unpack,
   int alength
)
{
  THETA theta, obs;

  int width, i, j;

  /* An index for stepping through the linear theta matrix. */
  int i_theta = 0;

  /* Get the theta matrix for this component. */
  theta = model->theta;
  obs = model->obs;

  /* Find the width of motif. */
  width = model->w;
  
  /*
  fprintf(stderr, "%d: alength=%d width=%d\n", mpMyID(), alength, width);
  */

  /* Iterate across positions. */
  for (i = 0; i < width; i++) {
    /* Iterate through the letters in the dataset. */
    for (j = 0; j < alength; j++) {

      /* Copy from theta into the theta matrix or vice versa. */
      if (pack_or_unpack == 1) {
	theta_matrix[i_theta] = theta(i, j);
	obs_matrix[i_theta++] = obs[i][j];
      } else {
	theta(i, j) = theta_matrix[i_theta];
	obs[i][j] = obs_matrix[i_theta++];
      }
    }
  }

} /* void theta_packer */

/**********************************************************************
 *
 * void max_packets
 *
 * Given two packets, select the best logev.
 * Ties are resolved in favor of:
 *   1) lowest start width
 *   2) lowest start nsites
 * This function is used in the reduction across models.