1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
|
#include "muscle.h"
#include "dpreglist.h"
#include "diaglist.h"
#include "pwpath.h"
#include "profile.h"
#include "timing.h"
#define TRACE 0
#define TRACE_PATH 0
#define LIST_DIAGS 0
static double g_dDPAreaWithoutDiags = 0.0;
static double g_dDPAreaWithDiags = 0.0;
static void OffsetPath(PWPath &Path, unsigned uOffsetA, unsigned uOffsetB)
{
const unsigned uEdgeCount = Path.GetEdgeCount();
for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)
{
const PWEdge &Edge = Path.GetEdge(uEdgeIndex);
// Nasty hack -- poke new values back into path, circumventing class
PWEdge &NonConstEdge = (PWEdge &) Edge;
NonConstEdge.uPrefixLengthA += uOffsetA;
NonConstEdge.uPrefixLengthB += uOffsetB;
}
}
static void DiagToPath(const Diag &d, PWPath &Path)
{
Path.Clear();
const unsigned uLength = d.m_uLength;
for (unsigned i = 0; i < uLength; ++i)
{
PWEdge Edge;
Edge.cType = 'M';
Edge.uPrefixLengthA = d.m_uStartPosA + i + 1;
Edge.uPrefixLengthB = d.m_uStartPosB + i + 1;
Path.AppendEdge(Edge);
}
}
static void AppendRegPath(PWPath &Path, const PWPath &RegPath)
{
const unsigned uRegEdgeCount = RegPath.GetEdgeCount();
for (unsigned uRegEdgeIndex = 0; uRegEdgeIndex < uRegEdgeCount; ++uRegEdgeIndex)
{
const PWEdge &RegEdge = RegPath.GetEdge(uRegEdgeIndex);
Path.AppendEdge(RegEdge);
}
}
SCORE GlobalAlignDiags(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
unsigned uLengthB, PWPath &Path)
{
#if LIST_DIAGS
TICKS t1 = GetClockTicks();
#endif
DiagList DL;
if (ALPHA_Amino == g_Alpha)
FindDiags(PA, uLengthA, PB, uLengthB, DL);
else if (ALPHA_DNA == g_Alpha || ALPHA_RNA == g_Alpha)
FindDiagsNuc(PA, uLengthA, PB, uLengthB, DL);
else
Quit("GlobalAlignDiags: bad alpha");
#if TRACE
Log("GlobalAlignDiags, diag list:\n");
DL.LogMe();
#endif
DL.Sort();
DL.DeleteIncompatible();
#if TRACE
Log("After DeleteIncompatible:\n");
DL.LogMe();
#endif
MergeDiags(DL);
#if TRACE
Log("After MergeDiags:\n");
DL.LogMe();
#endif
DPRegionList RL;
DiagListToDPRegionList(DL, RL, uLengthA, uLengthB);
#if TRACE
Log("RegionList:\n");
RL.LogMe();
#endif
#if LIST_DIAGS
{
TICKS t2 = GetClockTicks();
unsigned uArea = RL.GetDPArea();
Log("ticks=%ld\n", (long) (t2 - t1));
Log("area=%u\n", uArea);
}
#endif
g_dDPAreaWithoutDiags += uLengthA*uLengthB;
double dDPAreaWithDiags = 0.0;
const unsigned uRegionCount = RL.GetCount();
for (unsigned uRegionIndex = 0; uRegionIndex < uRegionCount; ++uRegionIndex)
{
const DPRegion &r = RL.Get(uRegionIndex);
PWPath RegPath;
if (DPREGIONTYPE_Diag == r.m_Type)
{
DiagToPath(r.m_Diag, RegPath);
#if TRACE_PATH
Log("DiagToPath, path=\n");
RegPath.LogMe();
#endif
}
else if (DPREGIONTYPE_Rect == r.m_Type)
{
const unsigned uRegStartPosA = r.m_Rect.m_uStartPosA;
const unsigned uRegStartPosB = r.m_Rect.m_uStartPosB;
const unsigned uRegLengthA = r.m_Rect.m_uLengthA;
const unsigned uRegLengthB = r.m_Rect.m_uLengthB;
const ProfPos *RegPA = PA + uRegStartPosA;
const ProfPos *RegPB = PB + uRegStartPosB;
dDPAreaWithDiags += uRegLengthA*uRegLengthB;
GlobalAlignNoDiags(RegPA, uRegLengthA, RegPB, uRegLengthB, RegPath);
#if TRACE_PATH
Log("GlobalAlignNoDiags RegPath=\n");
RegPath.LogMe();
#endif
OffsetPath(RegPath, uRegStartPosA, uRegStartPosB);
#if TRACE_PATH
Log("After offset path, RegPath=\n");
RegPath.LogMe();
#endif
}
else
Quit("GlobalAlignDiags, Invalid region type %u", r.m_Type);
AppendRegPath(Path, RegPath);
#if TRACE_PATH
Log("After AppendPath, path=");
Path.LogMe();
#endif
}
#if TRACE
{
double dDPAreaWithoutDiags = uLengthA*uLengthB;
Log("DP area with diags %.3g without %.3g pct saved %.3g %%\n",
dDPAreaWithDiags, dDPAreaWithoutDiags, (1.0 - dDPAreaWithDiags/dDPAreaWithoutDiags)*100.0);
}
#endif
g_dDPAreaWithDiags += dDPAreaWithDiags;
return 0;
}
void ListDiagSavings()
{
if (!g_bVerbose || !g_bDiags)
return;
double dAreaSaved = g_dDPAreaWithoutDiags - g_dDPAreaWithDiags;
double dPct = dAreaSaved*100.0/g_dDPAreaWithoutDiags;
Log("DP area saved by diagonals %-4.1f%%\n", dPct);
}
|
