binaries/src/muscle/readmx.cpp

   1 #include "muscle.h"\r
   2 #include "textfile.h"\r
   3 \r
   4 #define TRACE   0\r
   5 \r
   6 const int MAX_LINE = 4096;\r
   7 const int MAX_HEADINGS = 32;\r
   8 static char Heading[MAX_HEADINGS];\r
   9 static unsigned HeadingCount = 0;\r
  10 static float Mx[32][32];\r
  11 \r
  12 static void LogMx()\r
  13         {\r
  14         Log("Matrix\n");\r
  15         Log("     ");\r
  16         for (int i = 0; i < 20; ++i)\r
  17                 Log("    %c", LetterToChar(i));\r
  18         Log("\n");\r
  19 \r
  20         for (int i = 0; i < 20; ++i)\r
  21                 {\r
  22                 Log("%c    ", LetterToChar(i));\r
  23                 for (int j = 0; j < 20; ++j)\r
  24                         Log("%5.1f", Mx[i][j]);\r
  25                 Log("\n");\r
  26                 }\r
  27         Log("\n");\r
  28         }\r
  29 \r
  30 static unsigned MxCharToLetter(char c)\r
  31         {\r
  32         for (unsigned Letter = 0; Letter < HeadingCount; ++Letter)\r
  33                 if (Heading[Letter] == c)\r
  34                         return Letter;\r
  35         Quit("Letter '%c' has no heading", c);\r
  36         return 0;\r
  37         }\r
  38 \r
  39 PTR_SCOREMATRIX ReadMx(TextFile &File)\r
  40         {\r
  41 // Find column headers\r
  42         char Line[MAX_LINE];\r
  43         for (;;)\r
  44                 {\r
  45                 bool EndOfFile = File.GetLine(Line, sizeof(Line));\r
  46                 if (EndOfFile)\r
  47                         Quit("Premature EOF in matrix file");\r
  48 \r
  49                 if (Line[0] == '#')\r
  50                         continue;\r
  51                 else if (Line[0] == ' ')\r
  52                         break;\r
  53                 else\r
  54                         Quit("Invalid line in matrix file: '%s'", Line);\r
  55                 }\r
  56 \r
  57 // Read column headers\r
  58         HeadingCount = 0;\r
  59         for (char *p = Line; *p; ++p)\r
  60                 {\r
  61                 char c = *p;\r
  62                 if (!isspace(c))\r
  63                         Heading[HeadingCount++] = c;\r
  64                 }\r
  65 \r
  66         if (HeadingCount > 0 && Heading[HeadingCount-1] == '*')\r
  67                 --HeadingCount;\r
  68 \r
  69         if (HeadingCount < 20)\r
  70                 Quit("Error in matrix file: < 20 headers, line='%s'", Line);\r
  71 \r
  72 #if TRACE\r
  73         {\r
  74         Log("ReadMx\n");\r
  75         Log("%d headings: ", HeadingCount);\r
  76         for (unsigned i = 0; i < HeadingCount; ++i)\r
  77                 Log("%c", Heading[i]);\r
  78         Log("\n");\r
  79         }\r
  80 #endif\r
  81 \r
  82 // Zero out matrix\r
  83         for (int i = 0; i < MAX_ALPHA; ++i)\r
  84                 for (int j = 0; j < MAX_ALPHA; ++j)\r
  85                         Mx[i][j] = 0.0;\r
  86 \r
  87 // Read data lines\r
  88         for (unsigned RowIndex = 0; RowIndex < HeadingCount; ++RowIndex)\r
  89                 {\r
  90                 bool EndOfFile = File.GetTrimLine(Line, sizeof(Line));\r
  91                 if (EndOfFile)\r
  92                         Quit("Premature EOF in matrix file");\r
  93 #if     TRACE\r
  94                 Log("Line=%s\n", Line);\r
  95 #endif\r
  96                 if (Line[0] == '#')\r
  97                         continue;\r
  98 \r
  99                 char c = Line[0];\r
 100 #if     TRACE\r
 101                 Log("Row char=%c\n", c);\r
 102 #endif\r
 103                 if (!IsResidueChar(c))\r
 104                         continue;\r
 105                 unsigned RowLetter = CharToLetter(c);\r
 106                 if (RowLetter >= 20)\r
 107                         continue;\r
 108 #if     TRACE\r
 109                 Log("Row letter = %u\n", RowLetter);\r
 110 #endif\r
 111 \r
 112                 char *p = Line + 1;\r
 113                 char *maxp = p + strlen(Line);\r
 114                 for (unsigned Col = 0; Col < HeadingCount - 1; ++Col)\r
 115                         {\r
 116                         if (p >= maxp)\r
 117                                 Quit("Too few fields in line of matrix file: '%s'", Line);\r
 118                         while (isspace(*p))\r
 119                                 ++p;\r
 120                         char *Value = p;\r
 121                         while (!isspace(*p))\r
 122                                 ++p;\r
 123                         float v = (float) atof(Value);\r
 124                         char HeaderChar = Heading[Col];\r
 125                         if (IsResidueChar(HeaderChar))\r
 126                                 {\r
 127                                 unsigned ColLetter = CharToLetter(HeaderChar);\r
 128                                 if (ColLetter >= 20)\r
 129                                         continue;\r
 130                                 Mx[RowLetter][ColLetter] = v;\r
 131                                 }\r
 132                         p += 1;\r
 133                         }\r
 134                 }\r
 135 \r
 136 // Sanity check for symmetry\r
 137         for (int i = 0; i < 20; ++i)\r
 138                 for (int j = 0; j < i; ++j)\r
 139                         {\r
 140                         if (Mx[i][j] != Mx[j][i])\r
 141                                 {\r
 142                                 Warning("Matrix is not symmetrical, %c->%c=%g, %c->%c=%g",\r
 143                                   CharToLetter(i),\r
 144                                   CharToLetter(j),\r
 145                                   Mx[i][j],\r
 146                                   CharToLetter(j),\r
 147                                   CharToLetter(i),\r
 148                                   Mx[j][i]);\r
 149                                 goto ExitLoop;\r
 150                                 }\r
 151                         }\r
 152 ExitLoop:;\r
 153 \r
 154         if (g_bVerbose)\r
 155                 LogMx();\r
 156 \r
 157         return &Mx;\r
 158         }\r