mx.cpp

   1 //uchime by Robert C. Edgar http://drive5.com/uchime This code is donated to the public domain.\r
   2 \r
   3 #include "myutils.h"\r
   4 #include "mx.h"\r
   5 #include "seqdb.h"\r
   6 #include "seq.h"\r
   7 \r
   8 char ProbToChar(float p);\r
   9 \r
  10 list<MxBase *> *MxBase::m_Matrices = 0;\r
  11 unsigned MxBase::m_AllocCount;\r
  12 unsigned MxBase::m_ZeroAllocCount;\r
  13 unsigned MxBase::m_GrowAllocCount;\r
  14 double MxBase::m_TotalBytes;\r
  15 double MxBase::m_MaxBytes;\r
  16 \r
  17 static const char *LogizeStr(const char *s)\r
  18         {\r
  19         double d = atof(s);\r
  20         d = log(d);\r
  21         return TypeToStr<float>(float(d));\r
  22         }\r
  23 \r
  24 static const char *ExpizeStr(const char *s)\r
  25         {\r
  26         double d = atof(s);\r
  27         d = exp(d);\r
  28         return TypeToStr<float>(float(d));\r
  29         }\r
  30 \r
  31 void MxBase::OnCtor(MxBase *Mx)\r
  32         {\r
  33         if (m_Matrices == 0)\r
  34                 m_Matrices = new list<MxBase *>;\r
  35         asserta(m_Matrices != 0);\r
  36         m_Matrices->push_front(Mx);\r
  37         }\r
  38 \r
  39 void MxBase::OnDtor(MxBase *Mx)\r
  40         {\r
  41         if (m_Matrices == 0)\r
  42                 {\r
  43                 Warning("MxBase::OnDtor, m_Matrices = 0");\r
  44                 return;\r
  45                 }\r
  46         for (list<MxBase*>::iterator p = m_Matrices->begin();\r
  47           p != m_Matrices->end(); ++p)\r
  48                 {\r
  49                 if (*p == Mx)\r
  50                         {\r
  51                         m_Matrices->erase(p);\r
  52                         if (m_Matrices->empty())\r
  53                                 delete m_Matrices;\r
  54                         return;\r
  55                         }\r
  56                 }\r
  57         Warning("MxBase::OnDtor, not found");\r
  58         }\r
  59 \r
  60 //float **MxBase::Getf(const string &Name)\r
  61 //      {\r
  62 //      Mx<float> *m = (Mx<float> *) Get(Name);\r
  63 //      asserta(m->GetTypeSize() == sizeof(float));\r
  64 //      return m->GetData();\r
  65 //      }\r
  66 //\r
  67 //double **MxBase::Getd(const string &Name)\r
  68 //      {\r
  69 //      Mx<double> *m = (Mx<double> *) Get(Name);\r
  70 //      asserta(m->GetTypeSize() == sizeof(double));\r
  71 //      return m->GetData();\r
  72 //      }\r
  73 //\r
  74 //char **MxBase::Getc(const string &Name)\r
  75 //      {\r
  76 //      Mx<char> *m = (Mx<char> *) Get(Name);\r
  77 //      asserta(m->GetTypeSize() == sizeof(char));\r
  78 //      return m->GetData();\r
  79 //      }\r
  80 \r
  81 void MxBase::Alloc(const char *Name, unsigned RowCount, unsigned ColCount,\r
  82   const SeqDB *DB, unsigned IdA, unsigned IdB)\r
  83         {\r
  84         Alloc(Name, RowCount, ColCount, DB, IdA, IdB, 0, 0);\r
  85         }\r
  86 \r
  87 void MxBase::Alloc(const char *Name, unsigned RowCount, unsigned ColCount,\r
  88   const SeqData *SA, const SeqData *SB)\r
  89         {\r
  90         Alloc(Name, RowCount, ColCount, 0, UINT_MAX, UINT_MAX, SA, SB);\r
  91         }\r
  92 \r
  93 void MxBase::Alloc(const char *Name, unsigned RowCount, unsigned ColCount,\r
  94   const SeqDB *DB, unsigned IdA, unsigned IdB, const SeqData *SA, const SeqData *SB)\r
  95         {\r
  96         StartTimer(MxBase_Alloc);\r
  97 \r
  98         ++m_AllocCount;\r
  99         if (m_AllocatedRowCount == 0)\r
 100                 ++m_ZeroAllocCount;\r
 101 \r
 102         if (DB != 0)\r
 103                 {\r
 104                 asserta(IdA != UINT_MAX);\r
 105                 asserta(IdB != UINT_MAX);\r
 106                 asserta(RowCount >= DB->GetSeqLength(IdA) + 1);\r
 107                 asserta(ColCount >= DB->GetSeqLength(IdB) + 1);\r
 108                 }\r
 109         if (RowCount > m_AllocatedRowCount || ColCount > m_AllocatedColCount)\r
 110                 {\r
 111                 if (m_AllocatedRowCount > 0)\r
 112                         {\r
 113                         if (opt_logmemgrows)\r
 114                                 Log("MxBase::Alloc grow %s %u x %u -> %u x %u, %s bytes\n",\r
 115                                   Name, m_AllocatedRowCount, m_AllocatedColCount,\r
 116                                   RowCount, ColCount,\r
 117                                   IntToStr(GetBytes()));\r
 118                         ++m_GrowAllocCount;\r
 119                         }\r
 120 \r
 121                 m_TotalBytes -= GetBytes();\r
 122 \r
 123                 PauseTimer(MxBase_Alloc);\r
 124                 StartTimer(MxBase_FreeData);\r
 125                 FreeData();\r
 126                 EndTimer(MxBase_FreeData);\r
 127                 StartTimer(MxBase_Alloc);\r
 128 \r
 129                 unsigned N = max(RowCount + 16, m_AllocatedRowCount);\r
 130                 unsigned M = max(ColCount + 16, m_AllocatedColCount);\r
 131                 N = max(N, M);\r
 132 \r
 133                 PauseTimer(MxBase_Alloc);\r
 134                 StartTimer(MxBase_AllocData);\r
 135                 AllocData(N, N);\r
 136                 EndTimer(MxBase_AllocData);\r
 137                 StartTimer(MxBase_Alloc);\r
 138 \r
 139                 m_TotalBytes += GetBytes();\r
 140                 if (m_TotalBytes > m_MaxBytes)\r
 141                         m_MaxBytes = m_TotalBytes;\r
 142                 }\r
 143         \r
 144         unsigned n = sizeof(m_Name)-1;\r
 145         strncpy(m_Name, Name, n);\r
 146         m_Name[n] = 0;\r
 147         m_RowCount = RowCount;\r
 148         m_ColCount = ColCount;\r
 149         m_SeqDB = DB;\r
 150         m_IdA = IdA;\r
 151         m_IdB = IdB;\r
 152         m_SA = SA;\r
 153         m_SB = SB;\r
 154 \r
 155         EndTimer(MxBase_Alloc);\r
 156         }\r
 157 \r
 158 void MxBase::LogMe(bool WithData, int Opts) const\r
 159         {\r
 160         Log("\n");\r
 161         if (Opts & OPT_EXP)\r
 162                 Log("Exp ");\r
 163         else if (Opts & OPT_LOG)\r
 164                 Log("Log ");\r
 165         bool ZeroBased = ((Opts & OPT_ZERO_BASED) != 0);\r
 166         Log("%s(%p) Rows %u/%u, Cols %u/%u",\r
 167           m_Name, this,\r
 168           m_RowCount, m_AllocatedRowCount,\r
 169           m_ColCount, m_AllocatedColCount);\r
 170         if (m_SeqDB != 0 && m_IdA != UINT_MAX)\r
 171                 Log(", A=%s", m_SeqDB->GetLabel(m_IdA));\r
 172         else if (m_SA != 0)\r
 173                 Log(", A=%s", m_SA->Label);\r
 174         if (m_SeqDB != 0 && m_IdB != UINT_MAX)\r
 175                 Log(", B=%s", m_SeqDB->GetLabel(m_IdB));\r
 176         else if (m_SB != 0)\r
 177                 Log(", B=%s", m_SB->Label);\r
 178         Log("\n");\r
 179         if (!WithData || m_RowCount == 0 || m_ColCount == 0)\r
 180                 return;\r
 181 \r
 182         const char *z = GetAsStr(0, 0);\r
 183         unsigned Width = strlen(z);\r
 184         unsigned Mod = 1;\r
 185         for (unsigned i = 0; i < Width; ++i)\r
 186                 Mod *= 10;\r
 187 \r
 188         if (m_Alpha[0] != 0)\r
 189                 {\r
 190                 Log("// Alphabet=%s\n", m_Alpha);\r
 191                 Log("//      ");\r
 192                 unsigned n = strlen(m_Alpha);\r
 193                 for (unsigned j = 0; j < n; ++j)\r
 194                         Log(" %*c", Width, m_Alpha[j]);\r
 195                 Log("\n");\r
 196                 for (unsigned i = 0; i < n; ++i)\r
 197                         {\r
 198                         Log("/* %c */ {", m_Alpha[i]);\r
 199                         unsigned ci = m_Alpha[i];\r
 200                         for (unsigned j = 0; j < n; ++j)\r
 201                                 {\r
 202                                 unsigned cj = m_Alpha[j];\r
 203                                 Log("%s,", GetAsStr(ci, cj));\r
 204                                 }\r
 205                         Log("},  // %c\n", m_Alpha[i]);\r
 206                         }\r
 207                 return;\r
 208                 }\r
 209         else if (m_Alpha2[0] != 0)\r
 210                 {\r
 211                 unsigned n = strlen(m_Alpha2);\r
 212                 Log("// Alphabet=%s\n", m_Alpha2);\r
 213                 Log("//      ");\r
 214                 for (unsigned j = 0; j < n; ++j)\r
 215                         Log(" %*c", Width, m_Alpha2[j]);\r
 216                 Log("\n");\r
 217                 for (unsigned i = 0; i < n; ++i)\r
 218                         {\r
 219                         Log("/* %c */ {", m_Alpha2[i]);\r
 220                         unsigned ci = m_Alpha2[i];\r
 221                         for (unsigned j = 0; j < n; ++j)\r
 222                                 Log("%s,", GetAsStr(i, j));\r
 223                         Log("},  // %c\n", m_Alpha2[i]);\r
 224                         }\r
 225                 return;\r
 226                 }\r
 227 \r
 228         const byte *A = 0;\r
 229         const byte *B = 0;\r
 230         if (m_SeqDB != 0 && m_IdA != UINT_MAX)\r
 231                 A = m_SeqDB->GetSeq(m_IdA);\r
 232         else if (m_SA != 0)\r
 233                 A = m_SA->Seq;\r
 234         if (m_SeqDB != 0 && m_IdB != UINT_MAX)\r
 235                 B = m_SeqDB->GetSeq(m_IdB);\r
 236         else if (m_SB != 0)\r
 237                 B = m_SB->Seq;\r
 238 \r
 239         if (B != 0)\r
 240                 {\r
 241                 if (A != 0)\r
 242                         Log("  ");\r
 243                 Log("%5.5s", "");\r
 244                 if (ZeroBased)\r
 245                         for (unsigned j = 0; j < m_ColCount; ++j)\r
 246                                 Log("%*c", Width, B[j]);\r
 247                 else\r
 248                         for (unsigned j = 0; j < m_ColCount; ++j)\r
 249                                 Log("%*c", Width, j == 0 ? ' ' : B[j-1]);\r
 250                 Log("\n");\r
 251                 }\r
 252 \r
 253         if (A != 0)\r
 254                 Log("  ");\r
 255         Log("%5.5s", "");\r
 256         for (unsigned j = 0; j < m_ColCount; ++j)\r
 257                 Log("%*u", Width, j%Mod);\r
 258         Log("\n");\r
 259 \r
 260         for (unsigned i = 0; i < m_RowCount; ++i)\r
 261                 {\r
 262                 if (A != 0)\r
 263                         {\r
 264                         if (ZeroBased)\r
 265                                 Log("%c ", A[i]);\r
 266                         else\r
 267                                 Log("%c ", i == 0 ? ' ' : A[i-1]);\r
 268                         }\r
 269                 Log("%4u ", i);\r
 270                 \r
 271                 for (unsigned j = 0; j < m_ColCount; ++j)\r
 272                         {\r
 273                         const char *s = GetAsStr(i, j);\r
 274                         if (Opts & OPT_LOG)\r
 275                                 s = LogizeStr(s);\r
 276                         else if (Opts & OPT_EXP)\r
 277                                 s = ExpizeStr(s);\r
 278                         Log("%s", s);\r
 279                         }\r
 280                 Log("\n");\r
 281                 }\r
 282         }\r
 283 static unsigned g_MatrixFileCount;\r
 284 \r
 285 void MxBase::LogCounts()\r
 286         {\r
 287         Log("\n");\r
 288         Log("MxBase::LogCounts()\n");\r
 289         Log("      What           N\n");\r
 290         Log("----------  ----------\n");\r
 291         Log("    Allocs  %10u\n", m_AllocCount);\r
 292         Log("ZeroAllocs  %10u\n", m_ZeroAllocCount);\r
 293         Log("     Grows  %10u\n", m_GrowAllocCount);\r
 294         Log("     Bytes  %10.10s\n", MemBytesToStr(m_TotalBytes));\r
 295         Log(" Max bytes  %10.10s\n", MemBytesToStr(m_MaxBytes));\r
 296         }\r