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