uchime_src/mx.h

   1 #ifndef mx_h\r
   2 #define mx_h\r
   3 \r
   4 #include <list>\r
   5 #include <limits.h>\r
   6 #include <math.h>\r
   7 #include "timing.h"\r
   8 #include "myutils.h"\r
   9 \r
  10 const int OPT_LOG = 0x01;\r
  11 const int OPT_EXP = 0x02;\r
  12 const int OPT_ZERO_BASED = 0x04;\r
  13 const float MINUS_INFINITY = -9e9f;\r
  14 const float UNINIT = -8e8f;\r
  15 \r
  16 struct SeqData;\r
  17 \r
  18 template<class T> const char *TypeToStr(T t)\r
  19         {\r
  20         Die("Unspecialised TypeToStr() called");\r
  21         ureturn(0);\r
  22         }\r
  23 \r
  24 template<> inline const char *TypeToStr<unsigned short>(unsigned short f)\r
  25         {\r
  26         static char s[16];\r
  27 \r
  28         sprintf(s, "%12u", f);\r
  29         return s;\r
  30         }\r
  31 \r
  32 template<> inline const char *TypeToStr<short>(short f)\r
  33         {\r
  34         static char s[16];\r
  35 \r
  36         sprintf(s, "%12d", f);\r
  37         return s;\r
  38         }\r
  39 \r
  40 template<> inline const char *TypeToStr<int>(int f)\r
  41         {\r
  42         static char s[16];\r
  43 \r
  44         sprintf(s, "%5d", f);\r
  45         return s;\r
  46         }\r
  47 \r
  48 template<> inline const char *TypeToStr<float>(float f)\r
  49         {\r
  50         static char s[16];\r
  51 \r
  52         if (f == UNINIT)\r
  53                 sprintf(s, "%12.12s", "?");\r
  54         else if (f < MINUS_INFINITY/2)\r
  55                 sprintf(s, "%12.12s", "*");\r
  56         else if (f == 0.0f)\r
  57                 sprintf(s, "%12.12s", ".");\r
  58         else if (f >= -1e5 && f <= 1e5)\r
  59                 sprintf(s, "%12.5f", f);\r
  60         else\r
  61                 sprintf(s, "%12.4g", f);\r
  62         return s;\r
  63         }\r
  64 \r
  65 template<> inline const char *TypeToStr<double>(double f)\r
  66         {\r
  67         static char s[16];\r
  68 \r
  69         if (f < -1e9)\r
  70                 sprintf(s, "%12.12s", "*");\r
  71         else if (f == 0.0f)\r
  72                 sprintf(s, "%12.12s", ".");\r
  73         else if (f >= -1e-5 && f <= 1e5)\r
  74                 sprintf(s, "%12.5f", f);\r
  75         else\r
  76                 sprintf(s, "%12.4g", f);\r
  77         return s;\r
  78         }\r
  79 \r
  80 static inline const char *FloatToStr(float f, string &s)\r
  81         {\r
  82         s = TypeToStr<float>(f);\r
  83         return s.c_str();\r
  84         }\r
  85 \r
  86 template<> inline const char *TypeToStr<char>(char c)\r
  87         {\r
  88         static char s[2];\r
  89         s[0] = c;\r
  90         return s;\r
  91         }\r
  92 \r
  93 template<> inline const char *TypeToStr<byte>(byte c)\r
  94         {\r
  95         static char s[2];\r
  96         s[0] = c;\r
  97         return s;\r
  98         }\r
  99 \r
 100 template<> inline const char *TypeToStr<bool>(bool tof)\r
 101         {\r
 102         static char s[2];\r
 103         s[0] = tof ? 'T' : 'F';\r
 104         return s;\r
 105         }\r
 106 \r
 107 struct SeqDB;\r
 108 \r
 109 struct MxBase\r
 110         {\r
 111 private:\r
 112         MxBase(const MxBase &rhs);\r
 113         MxBase &operator=(const MxBase &rhs);\r
 114 \r
 115 public:\r
 116         char m_Name[32];\r
 117         char m_Alpha[32];\r
 118         char m_Alpha2[32];\r
 119         unsigned m_RowCount;\r
 120         unsigned m_ColCount;\r
 121         unsigned m_AllocatedRowCount;\r
 122         unsigned m_AllocatedColCount;\r
 123         const SeqDB *m_SeqDB;\r
 124         unsigned m_IdA;\r
 125         unsigned m_IdB;\r
 126         const SeqData *m_SA;\r
 127         const SeqData *m_SB;\r
 128 \r
 129         static list<MxBase *> *m_Matrices;\r
 130         //static MxBase *Get(const string &Name);\r
 131         //static float **Getf(const string &Name);\r
 132         //static double **Getd(const string &Name);\r
 133         //static char **Getc(const string &Name);\r
 134 \r
 135         static unsigned m_AllocCount;\r
 136         static unsigned m_ZeroAllocCount;\r
 137         static unsigned m_GrowAllocCount;\r
 138         static double m_TotalBytes;\r
 139         static double m_MaxBytes;\r
 140 \r
 141         static void OnCtor(MxBase *Mx);\r
 142         static void OnDtor(MxBase *Mx);\r
 143 \r
 144         MxBase()\r
 145                 {\r
 146                 m_AllocatedRowCount = 0;\r
 147                 m_AllocatedColCount = 0;\r
 148                 m_RowCount = 0;\r
 149                 m_ColCount = 0;\r
 150                 m_IdA = UINT_MAX;\r
 151                 m_IdB = UINT_MAX;\r
 152                 m_SeqDB = 0;\r
 153                 OnCtor(this);\r
 154                 }\r
 155         virtual ~MxBase()\r
 156                 {\r
 157                 OnDtor(this);\r
 158                 }\r
 159 \r
 160         virtual unsigned GetTypeSize() const = 0;\r
 161         virtual unsigned GetBytes() const = 0;\r
 162 \r
 163         void Clear()\r
 164                 {\r
 165                 FreeData();\r
 166                 m_AllocatedRowCount = 0;\r
 167                 m_AllocatedColCount = 0;\r
 168                 m_RowCount = 0;\r
 169                 m_ColCount = 0;\r
 170                 m_IdA = UINT_MAX;\r
 171                 m_IdB = UINT_MAX;\r
 172                 m_SA = 0;\r
 173                 m_SB = 0;\r
 174                 }\r
 175 \r
 176         bool Empty() const\r
 177                 {\r
 178                 return m_RowCount == 0;\r
 179                 }\r
 180 \r
 181         virtual void AllocData(unsigned RowCount, unsigned ColCount) = 0;\r
 182         virtual void FreeData() = 0;\r
 183         virtual const char *GetAsStr(unsigned i, unsigned j) const = 0;\r
 184 \r
 185         void SetAlpha(const char *Alpha)\r
 186                 {\r
 187                 unsigned n = sizeof(m_Alpha);\r
 188                 strncpy(m_Alpha, Alpha, n);\r
 189                 m_Alpha[n] = 0;\r
 190                 }\r
 191 \r
 192         void Alloc(const char *Name, unsigned RowCount, unsigned ColCount,\r
 193           const SeqDB *DB, unsigned IdA, unsigned IdB,\r
 194           const SeqData *SA, const SeqData *SB);\r
 195 \r
 196         void Alloc(const char *Name, unsigned RowCount, unsigned ColCount,\r
 197           const SeqDB *DB = 0, unsigned IdA = UINT_MAX, unsigned IdB = UINT_MAX);\r
 198 \r
 199         void Alloc(const char *Name, unsigned RowCount, unsigned ColCount,\r
 200           const SeqData *SA, const SeqData *SB);\r
 201 \r
 202         static void LogAll()\r
 203                 {\r
 204                 Log("\n");\r
 205                 if (m_Matrices == 0)\r
 206                         {\r
 207                         Log("MxBase::m_Matrices=0\n");\r
 208                         return;\r
 209                         }\r
 210                 Log("\n");\r
 211                 Log("AllRows  AllCols    Sz        MB  Name\n");\r
 212                 Log("-------  -------  ----  --------  ----\n");\r
 213                 double TotalMB = 0;\r
 214                 for (list<MxBase *>::const_iterator p = m_Matrices->begin();\r
 215                   p != m_Matrices->end(); ++p)\r
 216                         {\r
 217                         const MxBase *Mx = *p;\r
 218                         if (Mx == 0)\r
 219                                 continue;\r
 220                         //if (Mx->m_RowCount != 0 || ShowEmpty)\r
 221                         //      Mx->LogMe(WithData);\r
 222                         unsigned ar = Mx->m_AllocatedRowCount;\r
 223                         if (ar == 0)\r
 224                                 continue;\r
 225                         unsigned ac = Mx->m_AllocatedColCount;\r
 226                         unsigned sz = Mx->GetTypeSize();\r
 227                         double MB = (double) ar*(double) ac*(double) sz/1e6;\r
 228                         TotalMB += MB;\r
 229                         Log("%7u  %7u  %4u  %8.2f  %s\n", ar, ac, sz, MB, Mx->m_Name);\r
 230                         }\r
 231                 Log("                        --------\n");\r
 232                 Log("%7.7s  %7.7s  %4.4s  %8.2f\n", "", "", "", TotalMB);\r
 233                 }\r
 234 \r
 235         void LogMe(bool WithData = true, int Opts = 0) const;\r
 236         static void LogCounts();\r
 237         };\r
 238 \r
 239 template<class T> struct Mx : public MxBase\r
 240         {\r
 241 // Disable unimplemented stuff\r
 242 private:\r
 243         Mx(Mx &rhs);\r
 244         Mx &operator=(Mx &rhs);\r
 245         // const Mx &operator=(const Mx &rhs) const;\r
 246 \r
 247 public:\r
 248         T **m_Data;\r
 249 \r
 250         Mx()\r
 251                 {\r
 252                 m_Data = 0;\r
 253                 }\r
 254         \r
 255         ~Mx()\r
 256                 {\r
 257                 FreeData();\r
 258                 }\r
 259 \r
 260         virtual void AllocData(unsigned RowCount, unsigned ColCount)\r
 261                 {\r
 262                 if (opt_logmemgrows)\r
 263                         Log("MxBase::AllocData(%u,%u) %s bytes, Name=%s\n",\r
 264                           RowCount, ColCount, IntToStr(GetBytes()), m_Name);\r
 265                 // m_Data = myalloc<T *>(RowCount);\r
 266                 m_Data = MYALLOC(T *, RowCount, Mx);\r
 267                 for (unsigned i = 0; i < RowCount; ++i)\r
 268                         // m_Data[i] = myalloc<T>(ColCount);\r
 269                         m_Data[i] = MYALLOC(T, ColCount, Mx);\r
 270                 AddBytes("Mx_AllocData", RowCount*sizeof(T *) + RowCount*ColCount*sizeof(T));\r
 271 \r
 272                 m_AllocatedRowCount = RowCount;\r
 273                 m_AllocatedColCount = ColCount;\r
 274                 }\r
 275 \r
 276         virtual void FreeData()\r
 277                 {\r
 278                 for (unsigned i = 0; i < m_AllocatedRowCount; ++i)\r
 279                         MYFREE(m_Data[i], m_AllocatedColCount, Mx);\r
 280                 MYFREE(m_Data, m_AllocatedRowCount, Mx);\r
 281                 SubBytes("Mx_AllocData",\r
 282                   m_AllocatedRowCount*sizeof(T *) + m_AllocatedRowCount*m_AllocatedColCount*sizeof(T));\r
 283 \r
 284                 m_Data = 0;\r
 285                 m_RowCount = 0;\r
 286                 m_ColCount = 0;\r
 287                 m_AllocatedRowCount = 0;\r
 288                 m_AllocatedColCount = 0;\r
 289                 }\r
 290 \r
 291         T **GetData()\r
 292                 {\r
 293                 return (T **) m_Data;\r
 294                 }\r
 295 \r
 296         T Get(unsigned i, unsigned j) const\r
 297                 {\r
 298                 assert(i < m_RowCount);\r
 299                 assert(j < m_ColCount);\r
 300                 return m_Data[i][j];\r
 301                 }\r
 302 \r
 303         void Put(unsigned i, unsigned j, T x) const\r
 304                 {\r
 305                 assert(i < m_RowCount);\r
 306                 assert(j < m_ColCount);\r
 307                 m_Data[i][j] = x;\r
 308                 }\r
 309 \r
 310         T GetOffDiagAvgs(vector<T> &Avgs) const\r
 311                 {\r
 312                 if (m_RowCount != m_ColCount)\r
 313                         Die("GetOffDiagAvgs, not symmetrical");\r
 314                 Avgs.clear();\r
 315                 T Total = T(0);\r
 316                 for (unsigned i = 0; i < m_RowCount; ++i)\r
 317                         {\r
 318                         T Sum = T(0);\r
 319                         for (unsigned j = 0; j < m_ColCount; ++j)\r
 320                                 {\r
 321                                 if (j == i)\r
 322                                         continue;\r
 323                                 Sum += m_Data[i][j];\r
 324                                 }\r
 325                         T Avg = Sum/(m_RowCount-1);\r
 326                         Total += Avg;\r
 327                         Avgs.push_back(Avg);\r
 328                         }\r
 329                 return m_RowCount == 0 ? T(0) : Total/m_RowCount;\r
 330                 }\r
 331 \r
 332         unsigned GetTypeSize() const\r
 333                 {\r
 334                 return sizeof(T);\r
 335                 }\r
 336 \r
 337         virtual unsigned GetBytes() const\r
 338                 {\r
 339                 return m_AllocatedRowCount*m_AllocatedColCount*GetTypeSize() +\r
 340                   m_AllocatedRowCount*sizeof(T *);\r
 341                 }\r
 342 \r
 343         const char *GetAsStr(unsigned i, unsigned j) const\r
 344                 {\r
 345                 return TypeToStr<T>(Get(i, j));\r
 346                 }\r
 347 \r
 348         const T *const *const GetData() const\r
 349                 {\r
 350                 return (const T *const *) m_Data;\r
 351                 }\r
 352 \r
 353         void Copy(const Mx<T> &rhs)\r
 354                 {\r
 355                 Alloc("Copy", rhs.m_RowCount, rhs.m_ColCount, rhs.m_SeqDB, rhs.m_IdA, rhs.m_IdB);\r
 356                 const T * const *Data = rhs.GetData();\r
 357                 for (unsigned i = 0; i < m_RowCount; ++i)\r
 358                         for (unsigned j = 0; j < m_ColCount; ++j)\r
 359                                 m_Data[i][j] = Data[i][j];\r
 360                 }\r
 361 \r
 362         void Assign(T v)\r
 363                 {\r
 364                 for (unsigned i = 0; i < m_RowCount; ++i)\r
 365                         for (unsigned j = 0; j < m_ColCount; ++j)\r
 366                                 m_Data[i][j] = v;\r
 367                 }\r
 368 \r
 369         bool Eq(const Mx &rhs, bool Bwd = false) const\r
 370                 {\r
 371                 if (rhs.m_ColCount != m_ColCount)\r
 372                         return false;\r
 373                 if (rhs.m_RowCount != m_RowCount)\r
 374                         return false;\r
 375                 const T * const*d = rhs.GetData();\r
 376                 int i1 = Bwd ? m_RowCount : 0;\r
 377                 int j1 = Bwd ? m_ColCount : 0;\r
 378                 int i2 = Bwd ? -1 : m_RowCount;\r
 379                 int j2 = Bwd ? -1 : m_ColCount;\r
 380                 for (int i = i1; i != i2; Bwd ? --i : ++i)\r
 381                         for (int j = j1; j != j2; Bwd ? --j : ++j)\r
 382                                 {\r
 383                                 float x = m_Data[i][j];\r
 384                                 float y = d[i][j];\r
 385                                 if (x < -1e10 && y < -1e10)\r
 386                                         continue;\r
 387                                 if (!feq(x, y))\r
 388                                         {\r
 389                                         Warning("%s[%d][%d] = %g, %s = %g",\r
 390                                           m_Name, i, j, x, rhs.m_Name, y);\r
 391                                         return false;\r
 392                                         }\r
 393                                 }\r
 394                 return true;\r
 395                 }\r
 396 \r
 397         bool EqMask(const Mx &rhs, const Mx<bool> &Mask) const\r
 398                 {\r
 399                 if (rhs.m_ColCount != m_ColCount)\r
 400                         return false;\r
 401                 if (rhs.m_RowCount != m_RowCount)\r
 402                         return false;\r
 403 \r
 404                 if (Mask.m_ColCount != m_ColCount)\r
 405                         return false;\r
 406                 if (Mask.m_RowCount != m_RowCount)\r
 407                         return false;\r
 408 \r
 409                 const T * const*d = rhs.GetData();\r
 410                 bool Bwd = false;\r
 411                 int i1 = Bwd ? m_RowCount : 0;\r
 412                 int j1 = Bwd ? m_ColCount : 0;\r
 413                 int i2 = Bwd ? -1 : m_RowCount;\r
 414                 int j2 = Bwd ? -1 : m_ColCount;\r
 415                 for (int i = i1; i != i2; Bwd ? --i : ++i)\r
 416                         for (int j = j1; j != j2; Bwd ? --j : ++j)\r
 417                                 {\r
 418                                 if (!Mask.m_Data[i][j])\r
 419                                         continue;\r
 420                                 float x = m_Data[i][j];\r
 421                                 float y = d[i][j];\r
 422                                 if (x < -1e10 && y < -1e10)\r
 423                                         continue;\r
 424                                 if (!feq(x, y))\r
 425                                         {\r
 426                                         Warning("%s[%d][%d] = %g, %s = %g",\r
 427                                           m_Name, i, j, x, rhs.m_Name, y);\r
 428                                         return false;\r
 429                                         }\r
 430                                 }\r
 431                 return true;\r
 432                 }\r
 433 \r
 434         void Init(T v)\r
 435                 {\r
 436                 for (unsigned i = 0; i < m_RowCount; ++i)\r
 437                         for (unsigned j = 0; j < m_ColCount; ++j)\r
 438                                 m_Data[i][j] = v;\r
 439                 }\r
 440         };\r
 441 \r
 442 void WriteMx(const string &Name, Mx<float> &Mxf);\r
 443 \r
 444 template<class T> void ReserveMx(Mx<T> &Mxf, unsigned N = UINT_MAX)\r
 445         {\r
 446         if (Mxf.m_AllocatedRowCount > 0)\r
 447                 return;\r
 448         extern unsigned g_MaxInputSeqLength;\r
 449         if (N == UINT_MAX)\r
 450                 N = g_MaxInputSeqLength+1;\r
 451         Mxf.Alloc("(Reserved)", N, N);\r
 452         }\r
 453 \r
 454 #endif // mx_h\r