forester/archive/RIO/others/hmmer/src/hmmpfam-pvm.c

   1 /************************************************************
   2  * HMMER - Biological sequence analysis with profile HMMs
   3  * Copyright (C) 1992-1999 Washington University School of Medicine
   4  * All Rights Reserved
   5  *
   6  *     This source code is distributed under the terms of the
   7  *     GNU General Public License. See the files COPYING and LICENSE
   8  *     for details.
   9  ************************************************************/
  10
  11 #ifdef HMMER_PVM
  12
  13 /* hmmslave-pvm.c
  14  * SRE, Sun Jul 12 17:15:36 1998
  15  *
  16  * PVM slave for hmmpfam-pvm and hmmsearch-pvm.
  17  * RCS $Id: hmmpfam-pvm.c,v 1.1.1.1 2005/03/22 08:34:15 cmzmasek Exp $
  18  */
  19
  20 #include <stdio.h>
  21 #include <stdlib.h>
  22 #include <pvm3.h>
  23
  24 #include "version.h"
  25 #include "structs.h"            /* data structures, macros, #define's   */
  26 #include "config.h"             /* compile-time configuration constants */
  27 #include "funcs.h"              /* function declarations                */
  28 #include "globals.h"            /* alphabet global variables            */
  29 #include "squid.h"              /* general sequence analysis library    */
  30
  31 static void leave_pvm(void);
  32
  33 int
  34 main(void)
  35 {
  36   struct p7trace_s *tr;         /* traceback of an alignment               */
  37   int      master_tid;          /* PVM TID of our master */
  38   char    *hmmfile;             /* file to read HMM(s) from                */
  39   HMMFILE *hmmfp;               /* opened hmmfile for reading              */
  40   struct plan7_s *hmm;
  41   char    *seq;
  42   char    *dsq;
  43   int      len;
  44   int      nhmm;                /* number of HMM to work on                */
  45   float    sc;
  46   int      my_idx = -1;         /* my index, 0..nslaves-1 */
  47   double   pvalue;              /* Z*pvalue = Evalue                        */
  48   double   evalue;              /* upper bound on evalue                    */
  49   struct threshold_s thresh;    /* threshold settings                       */
  50   int      send_trace;          /* TRUE if score is significant             */
  51   int      do_xnu;              /* TRUE to do XNU filter on seq             */
  52   int      do_forward;          /* TRUE to use Forward() scores not Viterbi */
  53   int      do_null2;            /* TRUE to correct scores w/ ad hoc null2   */
  54   int      alphatype;           /* alphabet type, hmmAMINO or hmmNUCLEIC    */
  55   int      code;                /* return code after initialization         */
  56
  57
  58   SQD_DPRINTF1(("a slave reporting for duty!\n"));
  59
  60   /* Register leave_pvm() cleanup function so any exit() call
  61    * first calls pvm_exit().
  62    */
  63   if (atexit(leave_pvm) != 0) { pvm_exit(); Die("slave couldn't register leave_pvm()"); }
  64
  65   /*****************************************************************
  66    * initialization.
  67    * Master broadcasts to us:
  68    *     1) len of HMM file name        (int)
  69    *     2) name of HMM file            (string)
  70    *     3) length of sequence string   (int)
  71    *     4) sequence                    (string)
  72    *     5) globT threshold
  73    *     6) globE threshold
  74    *     7) Z
  75    *     8) autocut setting
  76    *     9) do_xnu flag
  77    *    10) do_forward flag
  78    *    11) do_null2 flag
  79    *    12) alphabet type
  80    * We receive the broadcast and open the files.
  81    ******************************************************************/
  82
  83   master_tid = pvm_parent();    /* who's our master? */
  84   SQD_DPRINTF1(("I know my master is %d\n", master_tid));
  85
  86   pvm_recv(master_tid, HMMPVM_INIT);
  87   pvm_upkint(&len, 1, 1);
  88   hmmfile = MallocOrDie(sizeof(char *) * (len+1));
  89   pvm_upkstr(hmmfile);
  90   pvm_upkint(&len, 1, 1);
  91   seq = MallocOrDie(sizeof(char *) * (len+1));
  92   pvm_upkstr(seq);
  93   pvm_upkfloat(&(thresh.globT), 1, 1);
  94   pvm_upkdouble(&(thresh.globE), 1, 1);
  95   pvm_upkint(&(thresh.Z), 1, 1);
  96   pvm_upkint((int *) &(thresh.autocut), 1, 1);
  97   pvm_upkint(&do_xnu, 1, 1);
  98   pvm_upkint(&do_forward, 1, 1);
  99   pvm_upkint(&do_null2, 1, 1);
 100   pvm_upkint(&alphatype, 1, 1);
 101   SQD_DPRINTF1(("My master has told me how to initialize, and I am happy.\n"));
 102
 103   SetAlphabet(alphatype);
 104                                 /* Open HMM file (maybe in HMMERDB) */
 105   code = HMMPVM_OK;
 106   if ((hmmfp = HMMFileOpen(hmmfile, "HMMERDB")) == NULL)
 107     code = HMMPVM_NO_HMMFILE;
 108   else if (hmmfp->ssi == NULL)
 109     code = HMMPVM_NO_INDEX;
 110
 111   /* report our status.
 112    */
 113   pvm_initsend(PvmDataDefault);
 114   pvm_pkint(&code, 1, 1);
 115   PVMPackString(RELEASE);       /* proofing against bug#1 */
 116   pvm_send(master_tid, HMMPVM_RESULTS);
 117   SQD_DPRINTF1(("I have told my master my initialization status and I await his command.\n"));
 118
 119   dsq = DigitizeSequence(seq, len);
 120   if (do_xnu) XNU(dsq, len);
 121
 122   /*****************************************************************
 123    * Main loop.
 124    * Receive an integer 0..nhmm-1 for which HMM to search against.
 125    * If we receive a -1, we shut down.
 126    *****************************************************************/
 127
 128   for (;;)
 129     {
 130       pvm_recv(master_tid, HMMPVM_WORK);
 131       pvm_upkint(&nhmm, 1, 1);
 132       if (my_idx < 0) my_idx = nhmm; /* first time thru, remember what index we are. */
 133
 134       if (nhmm == -1) { /* shutdown signal */
 135         SQD_DPRINTF1(("I've been told to shut down."));
 136         break;
 137       }
 138
 139       /* move to our assigned HMM in the HMM file, and read it
 140        */
 141       SQD_DPRINTF1(("The master says to do HMM #%d - I hear and obey\n", nhmm));
 142       if (! HMMFilePositionByIndex(hmmfp, nhmm)) Die("didn't position the HMM file");
 143       if (! HMMFileRead(hmmfp, &hmm))            Die("unexpected end of HMM file");
 144       if (hmm == NULL)                           Die("unexpected failure to parse HMM file");
 145       P7Logoddsify(hmm, TRUE);
 146
 147                         /* set Pfam specific score thresholds if needed */
 148       if (! SetAutocuts(&thresh, hmm))
 149         Die("HMM %s doesn't have the score cutoffs you wanted", hmm->name);
 150
 151       /* Score sequence, do alignment (Viterbi), recover trace
 152        */
 153       if (P7ViterbiSize(len, hmm->M) <= RAMLIMIT)
 154         {
 155           SQD_DPRINTF1(("P7Viterbi(): Estimated size %d Mb\n", P7ViterbiSize(len, hmm->M)));
 156           sc = P7Viterbi(dsq, len, hmm, &tr);
 157         }
 158       else
 159         {
 160           SQD_DPRINTF1(("P7SmallViterbi() called; %d Mb > %d\n", P7ViterbiSize(len, hmm->M), RAMLIMIT));
 161           sc = P7SmallViterbi(dsq, len, hmm, &tr);
 162         }
 163
 164       /* The Forward score override.
 165        * See comments in hmmpfam.c in serial version.
 166        */
 167       if (do_forward) {
 168         sc  = P7Forward(dsq, len, hmm, NULL);
 169         if (do_null2)   sc -= TraceScoreCorrection(hmm, tr, dsq);
 170       }
 171
 172       pvalue = PValue(hmm, sc);
 173       evalue = thresh.Z ? (double) thresh.Z * pvalue : (double) nhmm * pvalue;
 174       send_trace = (sc >= thresh.globT && evalue <= thresh.globE) ? 1 : 0;
 175
 176       /* return output
 177        */
 178       pvm_initsend(PvmDataDefault);
 179       pvm_pkint(&my_idx, 1, 1); /* tell master who we are */
 180       pvm_pkstr(hmm->name);     /* double check that we did the right thing */
 181       pvm_pkfloat(&sc, 1, 1);
 182       pvm_pkdouble(&pvalue, 1, 1);
 183       pvm_pkint(&send_trace, 1, 1); /* flag for whether a trace structure is coming */
 184       if (send_trace) PVMPackTrace(tr);
 185       pvm_send(master_tid, HMMPVM_RESULTS);
 186
 187       /* cleanup
 188        */
 189       FreePlan7(hmm);
 190       P7FreeTrace(tr);
 191     }
 192
 193   /***********************************************
 194    * Cleanup, return.
 195    ***********************************************/
 196
 197   HMMFileClose(hmmfp);
 198   free(seq);
 199   free(dsq);
 200   free(hmmfile);
 201   return 0;
 202 }
 203
 204
 205 /* Function: leave_pvm()
 206  *
 207  * Purpose:  Cleanup function, to deal with crashes. We register
 208  *           this function using atexit() so it gets called before
 209  *           the slave dies.
 210  */
 211 static void leave_pvm(void)
 212 {
 213   SQD_DPRINTF1(("slave leaving PVM.\n"));
 214   pvm_exit();
 215 }
 216
 217
 218
 219 #else /* if HMMER_PVM not defined: include a dummy */
 220
 221 #include <stdio.h>
 222 int main(void)
 223 {
 224   printf("hmmpfam-slave is disabled. PVM support was not compiled into HMMER.\n");
 225   exit(0);
 226 }
 227
 228 #endif
 229