1 #include "p7_config.h"
2 #include "easel.h"
3
4 #include "hmmer.h"
5 #include "p7_gbands.h"
6 #include "p7_gmxb.h"
7
8 P7_GMXB *
p7_gmxb_Create(P7_GBANDS * bnd)9 p7_gmxb_Create(P7_GBANDS *bnd)
10 {
11 P7_GMXB *gxb = NULL;
12 int status;
13
14 ESL_ALLOC(gxb, sizeof(P7_GMXB));
15 gxb->dp = NULL;
16 gxb->xmx = NULL;
17 gxb->bnd = bnd;
18 gxb->dalloc = 0;
19 gxb->xalloc = 0;
20
21 ESL_ALLOC(gxb->dp, sizeof(float) * bnd->ncell * p7G_NSCELLS); /* i.e. *3, for MID (0..2) */
22 ESL_ALLOC(gxb->xmx, sizeof(float) * bnd->nrow * p7G_NXCELLS); /* i.e. *5, for ENJBC (0..4) */
23 gxb->dalloc = bnd->ncell;
24 gxb->xalloc = bnd->nrow;
25 return gxb;
26
27 ERROR:
28 p7_gmxb_Destroy(gxb);
29 return NULL;
30 }
31
32 int
p7_gmxb_Reinit(P7_GMXB * gxb,P7_GBANDS * bnd)33 p7_gmxb_Reinit(P7_GMXB *gxb, P7_GBANDS *bnd)
34 {
35 int status;
36
37 if (bnd->ncell > gxb->dalloc) {
38 ESL_REALLOC(gxb->dp, sizeof(float) * bnd->ncell * p7G_NSCELLS);
39 gxb->dalloc = bnd->ncell;
40 }
41
42 if (bnd->nrow > gxb->xalloc) {
43 ESL_REALLOC(gxb->xmx, sizeof(float) * bnd->nrow * p7G_NXCELLS);
44 gxb->xalloc = bnd->nrow;
45 }
46
47 gxb->bnd = bnd;
48 return eslOK;
49
50 ERROR:
51 return status;
52 }
53
54 int
p7_gmxb_Reuse(P7_GMXB * gxb)55 p7_gmxb_Reuse(P7_GMXB *gxb)
56 {
57 gxb->bnd = NULL;
58 return eslOK;
59 }
60
61
62 void
p7_gmxb_Destroy(P7_GMXB * gxb)63 p7_gmxb_Destroy(P7_GMXB *gxb)
64 {
65 if (gxb)
66 {
67 if (gxb->dp) free(gxb->dp);
68 if (gxb->xmx) free(gxb->xmx);
69 /* gxb->bnd is a reference ptr copy; memory remains caller's responsibility */
70 free(gxb);
71 }
72 }
73
74 static void
print_val(FILE * ofp,float val,int width,int precision,int flags)75 print_val(FILE *ofp, float val, int width, int precision, int flags)
76 {
77 if (flags & p7_SHOW_LOG) val = log(val);
78 fprintf(ofp, "%*.*f ", width, precision, val);
79 }
80
81 int
p7_gmxb_Dump(FILE * ofp,P7_GMXB * gxb,int flags)82 p7_gmxb_Dump(FILE *ofp, P7_GMXB *gxb, int flags)
83 {
84 int g, i, k, x;
85 int *bnd_ip = gxb->bnd->imem;
86 int *bnd_kp = gxb->bnd->kmem;
87 float *dp = gxb->dp;
88 float *xp = gxb->xmx;
89 int M = gxb->bnd->M;
90 int ia, ib;
91 int ka, kb;
92 int width = 9;
93 int precision = 4;
94
95
96 /* Header */
97 fprintf(ofp, " ");
98 for (k = 0; k <= M; k++) fprintf(ofp, "%*d ", width, k);
99 if (! (flags & p7_HIDE_SPECIALS)) fprintf(ofp, "%*s %*s %*s %*s %*s\n", width, "E", width, "N", width, "J", width, "B", width, "C");
100
101 fprintf(ofp, " ");
102 for (k = 0; k <= M; k++) fprintf(ofp, "%*.*s ", width, width, "----------");
103 if (! (flags & p7_HIDE_SPECIALS)) fprintf(ofp, "%*.*s ", width, width, "----------");
104 fprintf(ofp, "\n");
105
106 i = 0;
107 for (g = 0; g < gxb->bnd->nseg; g++)
108 {
109 ia = *bnd_ip; bnd_ip++;
110 ib = *bnd_ip; bnd_ip++;
111
112 if (ia > i+1) fprintf(ofp, "...\n\n");
113
114 for (i = ia; i <= ib; i++)
115 {
116 ka = *bnd_kp; bnd_kp++;
117 kb = *bnd_kp; bnd_kp++;
118
119 /* match cells */
120 fprintf(ofp, "%3d M ", i);
121 for (k = 0; k < ka; k++) fprintf (ofp, "%*s ", width, ".....");
122 for ( ; k <= kb; k++) print_val(ofp, dp[(k-ka)*p7G_NSCELLS + p7G_M], width, precision, flags);
123 for ( ; k <= M; k++) fprintf (ofp, "%*s ", width, ".....");
124
125 /* ENJBC specials */
126 if (! (flags & p7_HIDE_SPECIALS)) {
127 for (x = 0; x < p7G_NXCELLS; x++) print_val(ofp, xp[x], width, precision, flags);
128 }
129 fprintf(ofp, "\n");
130
131 /* insert cells */
132 fprintf(ofp, "%3d I ", i);
133 for (k = 0; k < ka; k++) fprintf (ofp, "%*s ", width, ".....");
134 for ( ; k <= kb; k++) print_val(ofp, dp[(k-ka)*p7G_NSCELLS + p7G_I], width, precision, flags);
135 for ( ; k <= M; k++) fprintf (ofp, "%*s ", width, ".....");
136 fprintf(ofp, "\n");
137
138 /* delete cells */
139 fprintf(ofp, "%3d D ", i);
140 for (k = 0; k < ka; k++) fprintf(ofp, "%*s ", width, ".....");
141 for ( ; k <= kb; k++) print_val(ofp, dp[(k-ka)*p7G_NSCELLS + p7G_D], width, precision, flags);
142 for ( ; k <= M; k++) fprintf(ofp, "%*s ", width, ".....");
143 fprintf(ofp, "\n\n");
144
145 dp += p7G_NSCELLS * (kb-ka+1); /* skip ahead to next dp sparse "row" */
146 xp += p7G_NXCELLS;
147 }
148 }
149 if (i <= gxb->bnd->L) fprintf(ofp, "...\n");
150 return eslOK;
151 }
152