1 //uchime by Robert C. Edgar http://drive5.com/uchime This code is donated to the public domain.
\r
4 #include <float.h> // for FLT_MAX
\r
6 #include "alnparams.h"
\r
12 void SetNucSubstMx(double Match, double Mismatch);
\r
13 void ReadSubstMx(const string &FileName, Mx<float> &Mxf);
\r
15 extern Mx<float> g_SubstMxf;
\r
16 extern float **g_SubstMx;
\r
18 void AlnParams::Clear()
\r
21 LocalOpen = OBVIOUSLY_WRONG_PENALTY;
\r
22 LocalExt = OBVIOUSLY_WRONG_PENALTY;
\r
23 OpenA = OBVIOUSLY_WRONG_PENALTY;
\r
24 OpenB = OBVIOUSLY_WRONG_PENALTY;
\r
25 ExtA = OBVIOUSLY_WRONG_PENALTY;
\r
26 ExtB = OBVIOUSLY_WRONG_PENALTY;
\r
27 LOpenA = OBVIOUSLY_WRONG_PENALTY;
\r
28 LOpenB = OBVIOUSLY_WRONG_PENALTY;
\r
29 ROpenA = OBVIOUSLY_WRONG_PENALTY;
\r
30 ROpenB = OBVIOUSLY_WRONG_PENALTY;
\r
31 LExtA = OBVIOUSLY_WRONG_PENALTY;
\r
32 LExtB = OBVIOUSLY_WRONG_PENALTY;
\r
33 RExtA = OBVIOUSLY_WRONG_PENALTY;
\r
34 RExtB = OBVIOUSLY_WRONG_PENALTY;
\r
39 bool AlnParams::Is2() const
\r
43 if (OpenB != g || LOpenA != g || LOpenB != g || ROpenA != g || ROpenB != g)
\r
45 if (ExtB != e || LExtA != e || LExtB != e || RExtA != e || RExtB != e)
\r
50 bool AlnParams::Is4() const
\r
56 if (OpenB != g || LOpenA != tg || LOpenB != tg || ROpenA != tg || ROpenB != tg)
\r
58 if (ExtB != e || LExtA != te || LExtB != te || RExtA != te || RExtB != te)
\r
63 const char *AlnParams::GetType() const
\r
72 void AlnParams::Init2(const float * const *Mx, float Open, float Ext)
\r
75 OpenA = OpenB = LOpenA = LOpenB = ROpenA = ROpenB = Open;
\r
76 ExtA = ExtB = LExtA = LExtB = RExtA = RExtB = Ext;
\r
79 void AlnParams::SetLocal(float Open, float Ext)
\r
85 void AlnParams::Init4(const float * const *Mx, float Open, float Ext,
\r
86 float TermOpen, float TermExt)
\r
89 OpenA = OpenB = Open;
\r
90 LOpenA = LOpenB = ROpenA = ROpenB = TermOpen;
\r
92 LExtA = LExtB = RExtA = RExtB = TermExt;
\r
95 void AlnParams::Init(const AlnParams &AP, const HSPData &HSP,
\r
96 unsigned LA, unsigned LB)
\r
98 SubstMx = AP.SubstMx;
\r
106 LOpenA = AP.LOpenA;
\r
117 LOpenB = AP.LOpenB;
\r
126 if (HSP.RightA(LA))
\r
128 ROpenA = AP.ROpenA;
\r
137 if (HSP.RightB(LB))
\r
139 ROpenB = AP.ROpenB;
\r
149 void AlnParams::LogMe() const
\r
151 Log("AlnParams(%s)", GetType());
\r
153 Log(" g=%.1f e=%.1f", -OpenA, -ExtA);
\r
155 Log(" g=%.1f tg=%.1f e=%.1f te=%.1f", -OpenA, -ExtA, -LOpenA, -LExtA);
\r
158 " 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
159 OpenA, OpenB, LOpenA, LOpenB, ROpenA, ROpenB, ExtA, ExtB, LExtA, LExtB, RExtA, RExtB);
\r
164 Open/Ext format string is one or more:
\r
165 [<flag><flag>...]<value>
\r
167 Value is (positive) penalty or * (disabled).
\r
171 I Internal gaps (defafault internal and terminal).
\r
172 E End gaps (default internal and terminal).
\r
177 static void ParseGapStr(const string &s,
\r
178 float &QI, float &QL, float &QR,
\r
179 float &TI, float &TL, float &TR)
\r
191 const unsigned K = SIZE(s);
\r
193 float Value = FLT_MAX;
\r
194 for (unsigned i = 0; i <= K; ++i)
\r
196 char c = s.c_str()[i];
\r
197 if (c == 0 || c == '/')
\r
199 if (Value == FLT_MAX)
\r
200 Die("Invalid gap penalty string, missing penalty '%s'", s.c_str());
\r
201 if (!Q && !T && !I && !E && !L && !R)
\r
210 if (!E && !I && !L && !R)
\r
221 Die("Invalid gap penalty string (E and L or R) '%s'", s.c_str());
\r
256 if (Value != FLT_MAX)
\r
257 Die("Invalid gap penalty (* in floating point number) '%s'", s.c_str());
\r
258 Value = -MINUS_INFINITY;
\r
260 else if (isdigit(c))
\r
262 if (Value == -MINUS_INFINITY)
\r
263 Die("Invalid gap penalty (* in floating point number) '%s'", s.c_str());
\r
264 if (Value == FLT_MAX)
\r
269 Value += float(c - '0')/Dec;
\r
272 Value = Value*10 + (c - '0');
\r
277 Die("Invalid gap penalty (two decimal points) '%s'", s.c_str());
\r
303 Die("Invalid char '%c' in gap penalty string '%s'", c, s.c_str());
\r
309 void AlnParams::SetPenalties(const string &OpenStr, const string &ExtStr)
\r
311 ParseGapStr(OpenStr, OpenA, LOpenA, ROpenA, OpenB, LOpenB, ROpenB);
\r
312 ParseGapStr(ExtStr, ExtA, LExtA, RExtA, ExtB, LExtB, RExtB);
\r
315 void AlnParams::SetMxFromCmdLine(bool IsNucleo)
\r
318 SetNucSubstMx(opt_match, opt_mismatch);
\r
321 if (opt_matrix == "")
\r
323 SubstMxName = "BLOSUM62";
\r
328 ReadSubstMx(opt_matrix, g_SubstMxf);
\r
329 g_SubstMx = g_SubstMxf.GetData();
\r
330 g_SubstMxf.LogMe();
\r
331 SubstMxName = opt_matrix.c_str();
\r
334 SubstMx = g_SubstMx;
\r
335 asserta(SubstMx != 0);
\r
338 void AlnParams::InitFromCmdLine(bool IsNucleo)
\r
344 SetMxFromCmdLine(IsNucleo);
\r
347 if (optset_lopen || optset_lext)
\r
349 if (!optset_lopen || !optset_lext)
\r
350 Die("Must set both --lopen and --lext");
\r
351 if (opt_lopen < 0.0 || opt_lext < 0.0)
\r
352 Die("Invalid --lopen/--lext, gap penalties must be >= 0");
\r
353 SetLocal(float(-opt_lopen), float(-opt_lext));
\r
357 // Same penalties, if-statement to note could differ.
\r
359 SetLocal(-10.0f, -1.0f);
\r
361 SetLocal(-10.0f, -1.0f);
\r
366 Init4(g_SubstMx, -10.0, -1.0, -0.5, -0.5);
\r
368 Init4(g_SubstMx, -17.0, -1.0, -0.5, -0.5);
\r
369 SetPenalties(opt_gapopen, opt_gapext);
\r
372 float AlnParams::GetLocalOpen() const
\r
377 float AlnParams::GetLocalExt() const
\r
382 bool AlnParams::GetIsNucleo() const
\r
384 asserta(NucleoSet);
\r
388 unsigned GetWindexWordLength(bool Nucleo)
\r
400 static void Test1(const string &os, const string &es)
\r
404 Log("OpenStr %s\n", os.c_str());
\r
405 Log(" ExtStr %s\n", es.c_str());
\r
406 AP.SetPenalties(os, es);
\r
412 Test1("17I/0.5E", "1I/0.5E");
\r
413 Test1("17I/0.5L/0.4R", "1Q/2T");
\r
414 Test1("1QL/2QR/3QI/4TL/5TR/6TI", ".1QL/.2QR/.3QI/.4TL/.5TR/.6TI");
\r