]> git.donarmstrong.com Git - mothur.git/blob - uchime_src/diagbox.h
changes while testing
[mothur.git] / uchime_src / diagbox.h
1 #ifndef diagbox_h\r
2 #define diagbox_h\r
3 \r
4 struct DiagBox;\r
5 \r
6 void GetDiagBox(unsigned LA, unsigned LB, unsigned DiagLo, unsigned DiagHi, DiagBox &Box);\r
7 void GetDiagRange(unsigned LA, unsigned LB, unsigned d,\r
8   unsigned &mini, unsigned &minj, unsigned &maxi, unsigned &maxj);\r
9 void GetDiagLoHi(unsigned LA, unsigned LB, const char *Path,\r
10   unsigned &dlo, unsigned &dhi);\r
11 \r
12 struct DiagBox\r
13         {\r
14         DiagBox()\r
15                 {\r
16                 }\r
17 \r
18         DiagBox(unsigned LA_, unsigned LB_, unsigned DiagLo, unsigned DiagHi)\r
19                 {\r
20                 //GetDiagBox(LA, LB, DiagLo, DiagHi, *this);\r
21                 //Validate();\r
22                 Init(LA_, LB_, DiagLo, DiagHi);\r
23                 }\r
24 \r
25         void Init(unsigned LA_, unsigned LB_, unsigned DiagLo, unsigned DiagHi)\r
26                 {\r
27                 GetDiagBox(LA_, LB_, DiagLo, DiagHi, *this);\r
28                 Validate();\r
29                 }\r
30 \r
31         unsigned LA;\r
32         unsigned LB;\r
33 \r
34         unsigned dlo;\r
35         unsigned dhi;\r
36 \r
37         unsigned dlo_mini;\r
38         unsigned dlo_minj;\r
39 \r
40         unsigned dlo_maxi;\r
41         unsigned dlo_maxj;\r
42 \r
43         unsigned dhi_mini;\r
44         unsigned dhi_minj;\r
45 \r
46         unsigned dhi_maxi;\r
47         unsigned dhi_maxj;\r
48 \r
49         unsigned GetDiag(unsigned i, unsigned j) const\r
50                 {\r
51                 return LA - i + j;\r
52                 }\r
53 \r
54 // i, j are positions 0..LA-1, 0..LB-1.\r
55         bool InBox(unsigned i, unsigned j) const\r
56                 {\r
57                 unsigned d = GetDiag(i, j);\r
58                 return d >= dlo && d <= dhi;\r
59                 }\r
60 \r
61 /***\r
62 i, j are 0-based prefix lengths 0..LA, 0..LB.\r
63 \r
64 A full path is in the box iff all match pairs are in the box.\r
65 \r
66 A partial path that aligns a prefix of A to a prefix of B as\r
67 in D.P.) is in the box iff it is is the prefix of at least\r
68 one full path that is in the box.\r
69 \r
70 A D.P. matrix entry X[i][j] is in the box iff there is at\r
71 least one full path aligning the first i letters of A and\r
72 the first j letters of B ending in a column of type X, i.e.\r
73 if there exists a partial path in the box that ends in X.\r
74 \r
75 Assume terminals appear in all paths, and DI/ID forbidden.\r
76 \r
77 Intuitively seems that by these definitions D is in box iff\r
78 DM or MD is in box, I is in box iff IM or MI is in box.\r
79 Don't have proof..\r
80 ***/\r
81         bool InBoxDPM(unsigned i, unsigned j) const\r
82                 {\r
83         // Special case for M[0][0]\r
84                 if (i == 0 && j == 0)\r
85                         return true;\r
86                 if (i == 0 || j == 0)\r
87                         return false;\r
88                 unsigned d = GetDiag(i-1, j-1);\r
89                 return d >= dlo && d <= dhi;\r
90                 }\r
91 \r
92         bool InBoxDPD(unsigned i, unsigned j) const\r
93                 {\r
94                 bool MD = i == 0 ? false : InBoxDPM(i-1, j);\r
95                 bool DM = (i == LA || j == LB) ? false : InBoxDPM(i+1, j+1);\r
96                 return MD || DM;\r
97                 }\r
98 \r
99         bool InBoxDPI(unsigned i, unsigned j) const\r
100                 {\r
101                 bool MI = j == 0 ? false : InBoxDPM(i, j-1);\r
102                 bool IM = (i == LA || j == LB) ? false : InBoxDPM(i+1, j+1);\r
103                 return MI || IM;\r
104                 }\r
105 \r
106         // d = LA - i + j = 1 .. LA+LB-1\r
107         void Validate() const\r
108                 {\r
109                 asserta(dlo <= dhi);\r
110                 asserta(dlo >= GetDiag(LA-1, 0));\r
111                 asserta(dhi <= GetDiag(0, LB-1));\r
112 \r
113                 asserta(GetDiag(dlo_mini, dlo_minj) == dlo);\r
114                 asserta(GetDiag(dlo_maxi, dlo_maxj) == dlo);\r
115                 asserta(GetDiag(dhi_mini, dhi_minj) == dhi);\r
116                 asserta(GetDiag(dhi_maxi, dhi_maxj) == dhi);\r
117 \r
118                 asserta(dlo_mini >= dhi_mini);\r
119                 asserta(dlo_minj <= dhi_minj);\r
120                 asserta(dlo_maxi >= dhi_maxi);\r
121                 asserta(dlo_maxj <= dhi_maxj);\r
122                 }\r
123 \r
124         unsigned GetMini() const\r
125                 {\r
126                 return dhi_mini;\r
127                 }\r
128 \r
129         unsigned GetMaxi() const\r
130                 {\r
131                 return dlo_maxi;\r
132                 }\r
133 \r
134         unsigned GetMinj() const\r
135                 {\r
136                 return dlo_minj;\r
137                 }\r
138 \r
139         unsigned GetMaxj() const\r
140                 {\r
141                 return dhi_maxj;\r
142                 }\r
143 /***\r
144         i = 0..LA-1\r
145         j = 0..LB-1\r
146         d = LA - i + j = 1 .. LA+LB-1\r
147         j = d - LA + i\r
148         i = LA - d + j\r
149 ***/\r
150         void GetRange_j(unsigned i, unsigned &Startj, unsigned &Endj) const\r
151                 {\r
152         // j = d - LA + i\r
153                 if (dlo + i >= LA)\r
154                         Startj = dlo + i - LA;\r
155                 else\r
156                         Startj = 0;\r
157 \r
158                 if (Startj >= LB)\r
159                         Startj = LB - 1;\r
160 \r
161                 if (dhi + i + 1 >= LA)\r
162                         Endj = dhi + i + 1 - LA;\r
163                 else\r
164                         Endj = 0;\r
165 \r
166                 if (Endj > LB)\r
167                         Endj = LB;\r
168 \r
169                 asserta(Endj >= Startj);\r
170                 }\r
171 \r
172         void LogMe() const\r
173                 {\r
174                 Log("LA=%u LB=%d dlo(%u): (%u,%u)-(%u,%u) dhi(%u): (%u,%u)-(%u,%u) i=[%u-%u] j=[%u-%u]\n",\r
175                   LA, LB,\r
176                   dlo,\r
177                   dlo_mini, dlo_minj,\r
178                   dlo_maxi, dlo_maxj,\r
179                   dhi,\r
180                   dhi_mini, dhi_minj,\r
181                   dhi_maxi, dhi_maxj,\r
182                   GetMini(), GetMaxi(),\r
183                   GetMinj(), GetMaxj());\r
184                 }\r
185         };\r
186 \r
187 typedef const char *(*NWDIAG)(const byte *A, unsigned LA, const byte *B, unsigned LB,
188   unsigned DiagLo, unsigned DiagHi, bool LeftTerm, bool RightTerm);
189
190 const char *NWBandWrap(NWDIAG NW, const byte *A, unsigned LA, const byte *B, unsigned LB,
191   unsigned DiagLo, unsigned DiagHi, bool LeftTerm, bool RightTerm);
192 \r
193 #endif // diagbox_h\r