2 #include <float.h> // for FLT_MAX
\r
4 #include "alnparams.h"
\r
10 void SetNucSubstMx(double Match, double Mismatch);
\r
11 void ReadSubstMx(const string &FileName, Mx<float> &Mxf);
\r
13 extern Mx<float> g_SubstMxf;
14 extern float **g_SubstMx;
16 void AlnParams::Clear()
\r
19 LocalOpen = OBVIOUSLY_WRONG_PENALTY;
\r
20 LocalExt = OBVIOUSLY_WRONG_PENALTY;
\r
21 OpenA = OBVIOUSLY_WRONG_PENALTY;
\r
22 OpenB = OBVIOUSLY_WRONG_PENALTY;
\r
23 ExtA = OBVIOUSLY_WRONG_PENALTY;
\r
24 ExtB = OBVIOUSLY_WRONG_PENALTY;
\r
25 LOpenA = OBVIOUSLY_WRONG_PENALTY;
\r
26 LOpenB = OBVIOUSLY_WRONG_PENALTY;
\r
27 ROpenA = OBVIOUSLY_WRONG_PENALTY;
\r
28 ROpenB = OBVIOUSLY_WRONG_PENALTY;
\r
29 LExtA = OBVIOUSLY_WRONG_PENALTY;
\r
30 LExtB = OBVIOUSLY_WRONG_PENALTY;
\r
31 RExtA = OBVIOUSLY_WRONG_PENALTY;
\r
32 RExtB = OBVIOUSLY_WRONG_PENALTY;
\r
37 bool AlnParams::Is2() const
\r
41 if (OpenB != g || LOpenA != g || LOpenB != g || ROpenA != g || ROpenB != g)
\r
43 if (ExtB != e || LExtA != e || LExtB != e || RExtA != e || RExtB != e)
\r
48 bool AlnParams::Is4() const
\r
54 if (OpenB != g || LOpenA != tg || LOpenB != tg || ROpenA != tg || ROpenB != tg)
\r
56 if (ExtB != e || LExtA != te || LExtB != te || RExtA != te || RExtB != te)
\r
61 const char *AlnParams::GetType() const
\r
70 void AlnParams::Init2(const float * const *Mx, float Open, float Ext)
\r
73 OpenA = OpenB = LOpenA = LOpenB = ROpenA = ROpenB = Open;
\r
74 ExtA = ExtB = LExtA = LExtB = RExtA = RExtB = Ext;
\r
77 void AlnParams::SetLocal(float Open, float Ext)
\r
83 void AlnParams::Init4(const float * const *Mx, float Open, float Ext,
\r
84 float TermOpen, float TermExt)
\r
87 OpenA = OpenB = Open;
\r
88 LOpenA = LOpenB = ROpenA = ROpenB = TermOpen;
\r
90 LExtA = LExtB = RExtA = RExtB = TermExt;
\r
93 void AlnParams::Init(const AlnParams &AP, const HSPData &HSP,
\r
94 unsigned LA, unsigned LB)
\r
96 SubstMx = AP.SubstMx;
\r
104 LOpenA = AP.LOpenA;
\r
115 LOpenB = AP.LOpenB;
\r
124 if (HSP.RightA(LA))
\r
126 ROpenA = AP.ROpenA;
\r
135 if (HSP.RightB(LB))
\r
137 ROpenB = AP.ROpenB;
\r
147 void AlnParams::LogMe() const
\r
149 Log("AlnParams(%s)", GetType());
\r
151 Log(" g=%.1f e=%.1f", -OpenA, -ExtA);
\r
153 Log(" g=%.1f tg=%.1f e=%.1f te=%.1f", -OpenA, -ExtA, -LOpenA, -LExtA);
\r
156 " gA=%.1f gB=%.1f gAL=%.1f gBL=%.1f gAR=%.1f gBR=%.1f eA=%.1f eB=%.1f eAL=%.1f eBL=%.1f eAR=%.1f eBR=%.1f",
\r
157 OpenA, OpenB, LOpenA, LOpenB, ROpenA, ROpenB, ExtA, ExtB, LExtA, LExtB, RExtA, RExtB);
\r
162 Open/Ext format string is one or more:
\r
163 [<flag><flag>...]<value>
\r
165 Value is (positive) penalty or * (disabled).
\r
169 I Internal gaps (defafault internal and terminal).
\r
170 E End gaps (default internal and terminal).
\r
175 static void ParseGapStr(const string &s,
\r
176 float &QI, float &QL, float &QR,
\r
177 float &TI, float &TL, float &TR)
\r
189 const unsigned K = SIZE(s);
\r
191 float Value = FLT_MAX;
\r
192 for (unsigned i = 0; i <= K; ++i)
\r
194 char c = s.c_str()[i];
\r
195 if (c == 0 || c == '/')
\r
197 if (Value == FLT_MAX)
\r
198 Die("Invalid gap penalty string, missing penalty '%s'", s.c_str());
\r
199 if (!Q && !T && !I && !E && !L && !R)
\r
208 if (!E && !I && !L && !R)
\r
219 Die("Invalid gap penalty string (E and L or R) '%s'", s.c_str());
\r
254 if (Value != FLT_MAX)
\r
255 Die("Invalid gap penalty (* in floating point number) '%s'", s.c_str());
\r
256 Value = -MINUS_INFINITY;
\r
258 else if (isdigit(c))
\r
260 if (Value == -MINUS_INFINITY)
\r
261 Die("Invalid gap penalty (* in floating point number) '%s'", s.c_str());
\r
262 if (Value == FLT_MAX)
\r
267 Value += float(c - '0')/Dec;
\r
270 Value = Value*10 + (c - '0');
\r
275 Die("Invalid gap penalty (two decimal points) '%s'", s.c_str());
\r
301 Die("Invalid char '%c' in gap penalty string '%s'", c, s.c_str());
\r
307 void AlnParams::SetPenalties(const string &OpenStr, const string &ExtStr)
\r
309 ParseGapStr(OpenStr, OpenA, LOpenA, ROpenA, OpenB, LOpenB, ROpenB);
\r
310 ParseGapStr(ExtStr, ExtA, LExtA, RExtA, ExtB, LExtB, RExtB);
\r
313 void AlnParams::SetMxFromCmdLine(bool IsNucleo)
\r
316 SetNucSubstMx(opt_match, opt_mismatch);
319 if (opt_matrix == "")
\r
321 SubstMxName = "BLOSUM62";
\r
326 ReadSubstMx(opt_matrix, g_SubstMxf);
\r
327 g_SubstMx = g_SubstMxf.GetData();
\r
328 g_SubstMxf.LogMe();
\r
329 SubstMxName = opt_matrix.c_str();
\r
332 SubstMx = g_SubstMx;
\r
333 asserta(SubstMx != 0);
\r
336 void AlnParams::InitFromCmdLine(bool IsNucleo)
\r
342 SetMxFromCmdLine(IsNucleo);
\r
345 if (optset_lopen || optset_lext)
\r
347 if (!optset_lopen || !optset_lext)
\r
348 Die("Must set both --lopen and --lext");
\r
349 if (opt_lopen < 0.0 || opt_lext < 0.0)
\r
350 Die("Invalid --lopen/--lext, gap penalties must be >= 0");
\r
351 SetLocal(float(-opt_lopen), float(-opt_lext));
\r
355 // Same penalties, if-statement to note could differ.
\r
357 SetLocal(-10.0f, -1.0f);
\r
359 SetLocal(-10.0f, -1.0f);
\r
364 Init4(g_SubstMx, -10.0, -1.0, -0.5, -0.5);
366 Init4(g_SubstMx, -17.0, -1.0, -0.5, -0.5);
367 SetPenalties(opt_gapopen, opt_gapext);
\r
370 float AlnParams::GetLocalOpen() const
\r
375 float AlnParams::GetLocalExt() const
\r
380 bool AlnParams::GetIsNucleo() const
\r
382 asserta(NucleoSet);
\r
386 unsigned GetWindexWordLength(bool Nucleo)
\r
398 static void Test1(const string &os, const string &es)
\r
402 Log("OpenStr %s\n", os.c_str());
\r
403 Log(" ExtStr %s\n", es.c_str());
\r
404 AP.SetPenalties(os, es);
\r
410 Test1("17I/0.5E", "1I/0.5E");
\r
411 Test1("17I/0.5L/0.4R", "1Q/2T");
\r
412 Test1("1QL/2QR/3QI/4TL/5TR/6TI", ".1QL/.2QR/.3QI/.4TL/.5TR/.6TI");
\r