forester/archive/RIO/others/hmmer/src/plan9.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 /* plan9.c
  12  * SRE, Wed Apr  8 07:35:30 1998
  13  *
  14  * alloc, free, and initialization of old Plan9 (HMMER 1.x) functions.
  15  * Rescued from the wreckage of HMMER 1.9m code.
  16  */
  17
  18 #include <stdio.h>
  19 #include <stdlib.h>
  20 #include <string.h>
  21 #include <math.h>
  22 #include "squid.h"
  23 #include "config.h"
  24 #include "structs.h"
  25 #include "funcs.h"
  26
  27 #ifdef MEMDEBUG
  28 #include "dbmalloc.h"
  29 #endif
  30
  31
  32 struct plan9_s *
  33 P9AllocHMM(int M)                               /* length of model to make */
  34 {
  35   struct plan9_s *hmm;        /* RETURN: blank HMM */
  36
  37   hmm        = (struct plan9_s *)     MallocOrDie (sizeof(struct plan9_s));
  38   hmm->ins   = (struct basic_state *) MallocOrDie (sizeof(struct basic_state) * (M+2));
  39   hmm->del   = (struct basic_state *) MallocOrDie (sizeof(struct basic_state) * (M+2));
  40   hmm->mat   = (struct basic_state *) MallocOrDie (sizeof(struct basic_state) * (M+2));
  41   hmm->ref   = (char *)  MallocOrDie ((M+2) * sizeof(char));
  42   hmm->cs    = (char *)  MallocOrDie ((M+2) * sizeof(char));
  43   hmm->xray  = (float *) MallocOrDie ((M+2) * sizeof(float) * NINPUTS);
  44   hmm->M     = M;
  45   hmm->name  = Strdup("unnamed"); /* name is not optional. */
  46
  47   hmm->flags = 0;
  48   P9ZeroHMM(hmm);
  49   return hmm;
  50 }
  51 int
  52 P9FreeHMM(struct plan9_s *hmm)
  53 {
  54   if (hmm == NULL) return 0;
  55   free(hmm->ref);
  56   free(hmm->cs);
  57   free(hmm->xray);
  58   free(hmm->name);
  59   if (hmm->mat != NULL)  free (hmm->mat);
  60   if (hmm->ins != NULL)  free (hmm->ins);
  61   if (hmm->del != NULL)  free (hmm->del);
  62   free(hmm);
  63   return 1;
  64 }
  65
  66
  67 /* Function: P9ZeroHMM()
  68  *
  69  * Purpose:  Zero emission and transition counts in an HMM.
  70  */
  71 void
  72 P9ZeroHMM(struct plan9_s *hmm)
  73 {
  74   int k, ts, idx;
  75
  76   for (k = 0; k <= hmm->M+1; k++)
  77     {
  78       for (ts = 0; ts < 3; ts++)
  79         {
  80           hmm->mat[k].t[ts] = 0.0;
  81           hmm->ins[k].t[ts] = 0.0;
  82           hmm->del[k].t[ts] = 0.0;
  83         }
  84       for (idx = 0; idx < Alphabet_size; idx++)
  85         {
  86           hmm->mat[k].p[idx]   = 0.0;
  87           hmm->ins[k].p[idx]   = 0.0;
  88           hmm->del[k].p[idx]   = 0.0;
  89         }
  90     }
  91 }
  92
  93
  94
  95
  96
  97 /* Function: P9Renormalize()
  98  *
  99  * Normalize all P distributions so they sum to 1.
 100  * P distributions that are all 0, or contain negative
 101  * probabilities, are left untouched.
 102  *
 103  * Returns 1 on success, or 0 on failure.
 104  */
 105 void
 106 P9Renormalize(struct plan9_s *hmm)
 107 {
 108   int    k;                     /* counter for states                  */
 109
 110   for (k = 0; k <= hmm->M ; k++)
 111     {
 112                                 /* match state transition frequencies */
 113       FNorm(hmm->mat[k].t, 3);
 114       FNorm(hmm->ins[k].t, 3);
 115       if (k > 0) FNorm(hmm->del[k].t, 3);
 116
 117       if (k > 0) FNorm(hmm->mat[k].p, Alphabet_size);
 118       FNorm(hmm->ins[k].p, Alphabet_size);
 119     }
 120 }
 121
 122 /* Function: P9DefaultNullModel()
 123  *
 124  * Purpose:  Set up a default random sequence model, using
 125  *           global aafq[]'s for protein or 0.25 for nucleic
 126  *           acid. randomseq is alloc'ed in caller. Alphabet information
 127  *           must already be known.
 128  */
 129 void
 130 P9DefaultNullModel(float *null)
 131 {
 132   int x;
 133   if (Alphabet_type == hmmAMINO)
 134     for (x = 0; x < Alphabet_size; x++)
 135       null[x] = aafq[x];
 136   else if (Alphabet_type == hmmNUCLEIC)
 137     for (x = 0; x < Alphabet_size; x++)
 138       null[x] = 0.25;
 139   else
 140     Die("No support for non-protein, non-nucleic acid alphabets.");
 141 }