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