mx.h

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