1 /************************************************************
2  * HMMER - Biological sequence analysis with profile-HMMs
3  * Copyright (C) 1992-1998 Washington University School of Medicine
4  *
5  *   This source code is distributed under the terms of the
6  *   GNU General Public License. See the files COPYING and
7  *   GNULICENSE for details.
8  *
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 *
P9AllocHMM(int M)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
P9FreeHMM(struct plan9_s * hmm)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
P9ZeroHMM(struct plan9_s * hmm)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
P9Renormalize(struct plan9_s * hmm)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
P9DefaultNullModel(float * null)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 }
142