src/gromacs/fileio/matio.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team.
   7  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "matio.h"
  41
  42 #include <cctype>
  43 #include <cmath>
  44 #include <cstdio>
  45 #include <cstring>
  46
  47 #include <algorithm>
  48 #include <regex>
  49 #include <string>
  50
  51 #include "gromacs/fileio/gmxfio.h"
  52 #include "gromacs/math/functions.h"
  53 #include "gromacs/math/utilities.h"
  54 #include "gromacs/utility/binaryinformation.h"
  55 #include "gromacs/utility/cstringutil.h"
  56 #include "gromacs/utility/exceptions.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/futil.h"
  59 #include "gromacs/utility/programcontext.h"
  60 #include "gromacs/utility/smalloc.h"
  61
  62 using namespace gmx;
  63
  64 static const char mapper[] =
  65         "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<.>/"
  66         "?";
  67 #define NMAP static_cast<long int>(sizeof(mapper) / sizeof(mapper[0]))
  68
  69 #define MAX_XPM_LINELENGTH 4096
  70
  71 real** mk_matrix(int nx, int ny, gmx_bool b1D)
  72 {
  73     int    i;
  74     real** m;
  75
  76     snew(m, nx);
  77     if (b1D)
  78     {
  79         snew(m[0], nx * ny);
  80     }
  81
  82     for (i = 0; (i < nx); i++)
  83     {
  84         if (b1D)
  85         {
  86             m[i] = &(m[0][i * ny]);
  87         }
  88         else
  89         {
  90             snew(m[i], ny);
  91         }
  92     }
  93
  94     return m;
  95 }
  96
  97 void done_matrix(int nx, real*** m)
  98 {
  99     int i;
 100
 101     for (i = 0; (i < nx); i++)
 102     {
 103         sfree((*m)[i]);
 104     }
 105     sfree(*m);
 106     *m = nullptr;
 107 }
 108
 109 static bool operator==(t_xpmelmt e1, t_xpmelmt e2)
 110 {
 111     return (e1.c1 == e2.c1) && (e1.c2 == e2.c2);
 112 }
 113
 114 //! Return the index of the first element that matches \c c, or -1 if not found.
 115 t_matelmt searchcmap(ArrayRef<const t_mapping> map, t_xpmelmt c)
 116 {
 117     auto findIt = std::find_if(map.begin(), map.end(), [&c](const auto& m) { return (m.code == c); });
 118     return (findIt == map.end()) ? -1 : (findIt - map.begin());
 119 }
 120
 121 //! Read the mapping table from in, return number of entries
 122 static std::vector<t_mapping> getcmap(FILE* in, const char* fn)
 123 {
 124     int                    i, n;
 125     char                   line[STRLEN];
 126     char                   code[STRLEN], desc[STRLEN];
 127     double                 r, g, b;
 128     std::vector<t_mapping> m;
 129
 130     if (fgets2(line, STRLEN - 1, in) == nullptr)
 131     {
 132         gmx_fatal(FARGS,
 133                   "Not enough lines in colormap file %s"
 134                   "(just wanted to read number of entries)",
 135                   fn);
 136     }
 137     sscanf(line, "%d", &n);
 138     m.resize(n);
 139     for (i = 0; (i < n); i++)
 140     {
 141         if (fgets2(line, STRLEN - 1, in) == nullptr)
 142         {
 143             gmx_fatal(FARGS,
 144                       "Not enough lines in colormap file %s"
 145                       "(should be %d, found only %d)",
 146                       fn, n + 1, i);
 147         }
 148         sscanf(line, "%s%s%lf%lf%lf", code, desc, &r, &g, &b);
 149         m[i].code.c1 = code[0];
 150         m[i].code.c2 = 0;
 151         m[i].desc    = gmx_strdup(desc);
 152         m[i].rgb.r   = r;
 153         m[i].rgb.g   = g;
 154         m[i].rgb.b   = b;
 155     }
 156
 157     return m;
 158 }
 159
 160 std::vector<t_mapping> readcmap(const char* fn)
 161 {
 162     FilePtr in = openLibraryFile(fn);
 163     return getcmap(in.get(), fn);
 164 }
 165
 166 void printcmap(FILE* out, int n, t_mapping map[])
 167 {
 168     int i;
 169
 170     fprintf(out, "%d\n", n);
 171     for (i = 0; (i < n); i++)
 172     {
 173         fprintf(out, "%c%c  %20s  %10g  %10g  %10g\n", map[i].code.c1 ? map[i].code.c1 : ' ',
 174                 map[i].code.c2 ? map[i].code.c2 : ' ', map[i].desc, map[i].rgb.r, map[i].rgb.g,
 175                 map[i].rgb.b);
 176     }
 177 }
 178
 179 void writecmap(const char* fn, int n, t_mapping map[])
 180 {
 181     FILE* out;
 182
 183     out = gmx_fio_fopen(fn, "w");
 184     printcmap(out, n, map);
 185     gmx_fio_fclose(out);
 186 }
 187
 188 static char* fgetline(char** line, int llmax, int* llalloc, FILE* in)
 189 {
 190     char* fg;
 191
 192     if (llmax > *llalloc)
 193     {
 194         srenew(*line, llmax + 1);
 195         *llalloc = llmax;
 196     }
 197     fg = fgets(*line, llmax, in);
 198     trim(*line);
 199
 200     return fg;
 201 }
 202
 203 static void skipstr(char* line)
 204 {
 205     int i, c;
 206
 207     ltrim(line);
 208     c = 0;
 209     while ((line[c] != ' ') && (line[c] != '\0'))
 210     {
 211         c++;
 212     }
 213     i = c;
 214     while (line[c] != '\0')
 215     {
 216         line[c - i] = line[c];
 217         c++;
 218     }
 219     line[c - i] = '\0';
 220 }
 221
 222 static char* line2string(char** line)
 223 {
 224     int i;
 225
 226     if (*line != nullptr)
 227     {
 228         while (((*line)[0] != '\"') && ((*line)[0] != '\0'))
 229         {
 230             (*line)++;
 231         }
 232
 233         if ((*line)[0] != '\"')
 234         {
 235             return nullptr;
 236         }
 237         (*line)++;
 238
 239         i = 0;
 240         while (((*line)[i] != '\"') && ((*line)[i] != '\0'))
 241         {
 242             i++;
 243         }
 244
 245         if ((*line)[i] != '\"')
 246         {
 247             *line = nullptr;
 248         }
 249         else
 250         {
 251             (*line)[i] = 0;
 252         }
 253     }
 254
 255     return *line;
 256 }
 257
 258 //! If a label named \c label is found in \c line, return it. Otherwise return empty string.
 259 static std::string findLabelInLine(const std::string& line, const std::string& label)
 260 {
 261     std::regex  re(".*" + label + "\"(.*)\"");
 262     std::smatch match;
 263     if (std::regex_search(line, match, re) && match.size() > 1)
 264     {
 265         return match.str(1);
 266     }
 267     return std::string();
 268 }
 269
 270 //! Read and return a matrix from \c in
 271 static t_matrix read_xpm_entry(FILE* in)
 272 {
 273     char *       line_buf = nullptr, *line = nullptr, *str, buf[256] = { 0 };
 274     int          i, m, col_len, nch                                  = 0, llmax;
 275     int          llalloc = 0;
 276     unsigned int r, g, b;
 277     double       u;
 278     gmx_bool     bGetOnWithIt, bSetLine;
 279     t_xpmelmt    c;
 280
 281     t_matrix mm;
 282
 283     llmax = STRLEN;
 284
 285     while ((nullptr != fgetline(&line_buf, llmax, &llalloc, in))
 286            && (std::strncmp(line_buf, "static", 6) != 0))
 287     {
 288         std::string lineString = line_buf;
 289         mm.title               = findLabelInLine(lineString, "title");
 290         mm.legend              = findLabelInLine(lineString, "legend");
 291         mm.label_x             = findLabelInLine(lineString, "x-label");
 292         mm.label_y             = findLabelInLine(lineString, "y-label");
 293         findLabelInLine(lineString, "type"); // discard the returned string
 294     }
 295
 296     if (!line_buf || strncmp(line_buf, "static", 6) != 0)
 297     {
 298         gmx_input("Invalid XPixMap");
 299     }
 300
 301     if (buf[0] && (gmx_strcasecmp(buf, "Discrete") == 0))
 302     {
 303         mm.bDiscrete = TRUE;
 304     }
 305
 306     if (debug)
 307     {
 308         fprintf(debug, "%s %s %s %s\n", mm.title.c_str(), mm.legend.c_str(), mm.label_x.c_str(),
 309                 mm.label_y.c_str());
 310     }
 311
 312     /* Read sizes */
 313     int nmap     = 0;
 314     bGetOnWithIt = FALSE;
 315     while (!bGetOnWithIt && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)))
 316     {
 317         line = line_buf;
 318         while ((line[0] != '\"') && (line[0] != '\0'))
 319         {
 320             line++;
 321         }
 322
 323         if (line[0] == '\"')
 324         {
 325             line2string(&line);
 326             sscanf(line, "%d %d %d %d", &(mm.nx), &(mm.ny), &nmap, &nch);
 327             if (nch > 2)
 328             {
 329                 gmx_fatal(FARGS, "Sorry can only read xpm's with at most 2 characters per pixel\n");
 330             }
 331             if (mm.nx <= 0 || mm.ny <= 0)
 332             {
 333                 gmx_fatal(FARGS, "Dimensions of xpm-file have to be larger than 0\n");
 334             }
 335             llmax        = std::max(STRLEN, mm.nx + 10);
 336             bGetOnWithIt = TRUE;
 337         }
 338     }
 339     if (debug)
 340     {
 341         fprintf(debug, "mm.nx %d mm.ny %d nmap %d nch %d\n", mm.nx, mm.ny, nmap, nch);
 342     }
 343     if (nch == 0)
 344     {
 345         gmx_fatal(FARGS, "Number of characters per pixel not found in xpm\n");
 346     }
 347
 348     /* Read color map */
 349     mm.map.resize(nmap);
 350     m = 0;
 351     while ((m < nmap) && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)))
 352     {
 353         line = std::strchr(line_buf, '\"');
 354         if (line)
 355         {
 356             line++;
 357             /* Read xpm color map entry */
 358             mm.map[m].code.c1 = line[0];
 359             if (nch == 1)
 360             {
 361                 mm.map[m].code.c2 = 0;
 362             }
 363             else
 364             {
 365                 mm.map[m].code.c2 = line[1];
 366             }
 367             line += nch;
 368             str = std::strchr(line, '#');
 369             if (str)
 370             {
 371                 str++;
 372                 col_len = 0;
 373                 while (std::isxdigit(str[col_len]))
 374                 {
 375                     col_len++;
 376                 }
 377                 if (col_len == 6)
 378                 {
 379                     sscanf(line, "%*s #%2x%2x%2x", &r, &g, &b);
 380                     mm.map[m].rgb.r = r / 255.0;
 381                     mm.map[m].rgb.g = g / 255.0;
 382                     mm.map[m].rgb.b = b / 255.0;
 383                 }
 384                 else if (col_len == 12)
 385                 {
 386                     sscanf(line, "%*s #%4x%4x%4x", &r, &g, &b);
 387                     mm.map[m].rgb.r = r / 65535.0;
 388                     mm.map[m].rgb.g = g / 65535.0;
 389                     mm.map[m].rgb.b = b / 65535.0;
 390                 }
 391                 else
 392                 {
 393                     gmx_file("Unsupported or invalid colormap in X PixMap");
 394                 }
 395             }
 396             else
 397             {
 398                 str = std::strchr(line, 'c');
 399                 if (str)
 400                 {
 401                     str += 2;
 402                 }
 403                 else
 404                 {
 405                     gmx_file("Unsupported or invalid colormap in X PixMap");
 406                 }
 407                 fprintf(stderr, "Using white for color \"%s", str);
 408                 mm.map[m].rgb.r = 1;
 409                 mm.map[m].rgb.g = 1;
 410                 mm.map[m].rgb.b = 1;
 411             }
 412             line = std::strchr(line, '\"');
 413             line++;
 414             line2string(&line);
 415             mm.map[m].desc = gmx_strdup(line);
 416             m++;
 417         }
 418     }
 419     if (m != nmap)
 420     {
 421         gmx_fatal(FARGS,
 422                   "Number of read colors map entries (%d) does not match the number in the header "
 423                   "(%d)",
 424                   m, nmap);
 425     }
 426
 427     /* Read axes, if there are any */
 428     bSetLine = FALSE;
 429     do
 430     {
 431         if (bSetLine)
 432         {
 433             line = line_buf;
 434         }
 435         bSetLine = TRUE;
 436         GMX_RELEASE_ASSERT(line, "Need to have valid line to parse");
 437         if (strstr(line, "x-axis"))
 438         {
 439             line = std::strstr(line, "x-axis");
 440             skipstr(line);
 441             mm.axis_x.resize(0);
 442             mm.axis_x.reserve(mm.nx + 1);
 443             while (sscanf(line, "%lf", &u) == 1)
 444             {
 445                 if (ssize(mm.axis_x) > mm.nx)
 446                 {
 447                     gmx_fatal(FARGS, "Too many x-axis labels in xpm (max %d)", mm.nx);
 448                 }
 449                 else if (ssize(mm.axis_x) == mm.nx)
 450                 {
 451                     mm.flags |= MAT_SPATIAL_X;
 452                 }
 453                 mm.axis_x.push_back(u);
 454                 skipstr(line);
 455             }
 456         }
 457         else if (std::strstr(line, "y-axis"))
 458         {
 459             line = std::strstr(line, "y-axis");
 460             skipstr(line);
 461             mm.axis_y.resize(0);
 462             mm.axis_y.reserve(mm.ny + 1);
 463             while (sscanf(line, "%lf", &u) == 1)
 464             {
 465                 if (ssize(mm.axis_y) > mm.ny)
 466                 {
 467                     gmx_fatal(FARGS, "Too many y-axis labels in xpm (max %d)", mm.ny);
 468                 }
 469                 else if (ssize(mm.axis_y) == mm.ny)
 470                 {
 471                     mm.flags |= MAT_SPATIAL_Y;
 472                 }
 473                 mm.axis_y.push_back(u);
 474                 skipstr(line);
 475             }
 476         }
 477     } while ((line[0] != '\"') && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)));
 478
 479     /* Read matrix */
 480     mm.matrix.resize(mm.nx, mm.ny);
 481     int rowIndex = mm.ny - 1;
 482     bSetLine     = FALSE;
 483     do
 484     {
 485         if (bSetLine)
 486         {
 487             line = line_buf;
 488         }
 489         bSetLine = TRUE;
 490         if (rowIndex % (1 + mm.ny / 100) == 0)
 491         {
 492             fprintf(stderr, "%3d%%\b\b\b\b", (100 * (mm.ny - rowIndex)) / mm.ny);
 493         }
 494         while ((line[0] != '\"') && (line[0] != '\0'))
 495         {
 496             line++;
 497         }
 498         if (line[0] != '\"')
 499         {
 500             gmx_fatal(FARGS, "Not enough characters in row %d of the matrix\n", rowIndex + 1);
 501         }
 502         else
 503         {
 504             line++;
 505             for (i = 0; i < mm.nx; i++)
 506             {
 507                 c.c1 = line[nch * i];
 508                 if (nch == 1)
 509                 {
 510                     c.c2 = 0;
 511                 }
 512                 else
 513                 {
 514                     c.c2 = line[nch * i + 1];
 515                 }
 516                 mm.matrix(i, rowIndex) = searchcmap(mm.map, c);
 517             }
 518             rowIndex--;
 519         }
 520     } while ((rowIndex >= 0) && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)));
 521     if (rowIndex >= 0)
 522     {
 523         gmx_incons("Not enough rows in the matrix");
 524     }
 525
 526     sfree(line_buf);
 527     return mm;
 528 }
 529
 530 std::vector<t_matrix> read_xpm_matrix(const char* fnm)
 531 {
 532     FILE* in;
 533     char* line    = nullptr;
 534     int   llalloc = 0;
 535
 536     in = gmx_fio_fopen(fnm, "r");
 537
 538     std::vector<t_matrix> mat;
 539     while (nullptr != fgetline(&line, STRLEN, &llalloc, in))
 540     {
 541         if (std::strstr(line, "/* XPM */"))
 542         {
 543             mat.emplace_back(read_xpm_entry(in));
 544         }
 545     }
 546     gmx_fio_fclose(in);
 547
 548     if (mat.empty())
 549     {
 550         gmx_file("Invalid XPixMap");
 551     }
 552
 553     sfree(line);
 554
 555     return mat;
 556 }
 557
 558 real** matrix2real(t_matrix* in, real** out)
 559 {
 560     double tmp;
 561
 562     std::vector<real> rmap(in->map.size());
 563
 564     for (gmx::index i = 0; i != ssize(in->map); ++i)
 565     {
 566         if ((in->map[i].desc == nullptr) || (sscanf(in->map[i].desc, "%lf", &tmp) != 1))
 567         {
 568             fprintf(stderr,
 569                     "Could not convert matrix to reals,\n"
 570                     "color map entry %zd has a non-real description: \"%s\"\n",
 571                     i, in->map[i].desc);
 572             return nullptr;
 573         }
 574         rmap[i] = tmp;
 575     }
 576
 577     if (out == nullptr)
 578     {
 579         snew(out, in->nx);
 580         for (int i = 0; i < in->nx; i++)
 581         {
 582             snew(out[i], in->ny);
 583         }
 584     }
 585     for (int i = 0; i < in->nx; i++)
 586     {
 587         for (int j = 0; j < in->ny; j++)
 588         {
 589             out[i][j] = rmap[in->matrix(i, j)];
 590         }
 591     }
 592
 593     fprintf(stderr, "Converted a %dx%d matrix with %zu levels to reals\n", in->nx, in->ny, in->map.size());
 594
 595     return out;
 596 }
 597
 598 static void write_xpm_header(FILE*              out,
 599                              const std::string& title,
 600                              const std::string& legend,
 601                              const std::string& label_x,
 602                              const std::string& label_y,
 603                              gmx_bool           bDiscrete)
 604 {
 605     fprintf(out, "/* XPM */\n");
 606     fprintf(out, "/* This file can be converted to EPS by the GROMACS program xpm2ps */\n");
 607     fprintf(out, "/* title:   \"%s\" */\n", title.c_str());
 608     fprintf(out, "/* legend:  \"%s\" */\n", legend.c_str());
 609     fprintf(out, "/* x-label: \"%s\" */\n", label_x.c_str());
 610     fprintf(out, "/* y-label: \"%s\" */\n", label_y.c_str());
 611     if (bDiscrete)
 612     {
 613         fprintf(out, "/* type:    \"Discrete\" */\n");
 614     }
 615     else
 616     {
 617         fprintf(out, "/* type:    \"Continuous\" */\n");
 618     }
 619 }
 620
 621 static int calc_nmid(int nlevels, real lo, real mid, real hi)
 622 {
 623     /* Take care that we have at least 1 entry in the mid to hi range
 624      */
 625     return std::min(std::max(0, static_cast<int>(((mid - lo) / (hi - lo)) * (nlevels - 1))), nlevels - 1);
 626 }
 627
 628 static void
 629 write_xpm_map3(FILE* out, int n_x, int n_y, int* nlevels, real lo, real mid, real hi, t_rgb rlo, t_rgb rmid, t_rgb rhi)
 630 {
 631     int    i, nmid;
 632     double r, g, b, clev_lo, clev_hi;
 633
 634     if (*nlevels > NMAP * NMAP)
 635     {
 636         fprintf(stderr, "Warning, too many levels (%d) in matrix, using %d only\n", *nlevels,
 637                 static_cast<int>(NMAP * NMAP));
 638         *nlevels = NMAP * NMAP;
 639     }
 640     else if (*nlevels < 2)
 641     {
 642         fprintf(stderr, "Warning, too few levels (%d) in matrix, using 2 instead\n", *nlevels);
 643         *nlevels = 2;
 644     }
 645     if (!((mid >= lo) && (mid < hi)))
 646     {
 647         gmx_fatal(FARGS, "Lo: %f, Mid: %f, Hi: %f\n", lo, mid, hi);
 648     }
 649
 650     fprintf(out, "static char *gromacs_xpm[] = {\n");
 651     fprintf(out, "\"%d %d   %d %d\",\n", n_x, n_y, *nlevels, (*nlevels <= NMAP) ? 1 : 2);
 652
 653     nmid    = calc_nmid(*nlevels, lo, mid, hi);
 654     clev_lo = nmid;
 655     clev_hi = (*nlevels - 1 - nmid);
 656     for (i = 0; (i < nmid); i++)
 657     {
 658         r = rlo.r + (i * (rmid.r - rlo.r) / clev_lo);
 659         g = rlo.g + (i * (rmid.g - rlo.g) / clev_lo);
 660         b = rlo.b + (i * (rmid.b - rlo.b) / clev_lo);
 661         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n", mapper[i % NMAP],
 662                 (*nlevels <= NMAP) ? ' ' : mapper[i / NMAP],
 663                 static_cast<unsigned int>(std::round(255 * r)),
 664                 static_cast<unsigned int>(std::round(255 * g)),
 665                 static_cast<unsigned int>(std::round(255 * b)), ((nmid - i) * lo + i * mid) / clev_lo);
 666     }
 667     for (i = 0; (i < (*nlevels - nmid)); i++)
 668     {
 669         r = rmid.r + (i * (rhi.r - rmid.r) / clev_hi);
 670         g = rmid.g + (i * (rhi.g - rmid.g) / clev_hi);
 671         b = rmid.b + (i * (rhi.b - rmid.b) / clev_hi);
 672         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n", mapper[(i + nmid) % NMAP],
 673                 (*nlevels <= NMAP) ? ' ' : mapper[(i + nmid) / NMAP],
 674                 static_cast<unsigned int>(std::round(255 * r)),
 675                 static_cast<unsigned int>(std::round(255 * g)),
 676                 static_cast<unsigned int>(std::round(255 * b)),
 677                 ((*nlevels - 1 - nmid - i) * mid + i * hi) / clev_hi);
 678     }
 679 }
 680
 681 static void pr_simple_cmap(FILE* out, real lo, real hi, int nlevel, t_rgb rlo, t_rgb rhi, int i0)
 682 {
 683     int  i;
 684     real r, g, b, fac;
 685
 686     for (i = 0; (i < nlevel); i++)
 687     {
 688         fac = (i + 1.0) / (nlevel);
 689         r   = rlo.r + fac * (rhi.r - rlo.r);
 690         g   = rlo.g + fac * (rhi.g - rlo.g);
 691         b   = rlo.b + fac * (rhi.b - rlo.b);
 692         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n", mapper[(i + i0) % NMAP],
 693                 (nlevel <= NMAP) ? ' ' : mapper[(i + i0) / NMAP],
 694                 static_cast<unsigned int>(std::round(255 * r)),
 695                 static_cast<unsigned int>(std::round(255 * g)),
 696                 static_cast<unsigned int>(std::round(255 * b)), lo + fac * (hi - lo));
 697     }
 698 }
 699
 700 static void pr_discrete_cmap(FILE* out, int* nlevel, int i0)
 701 {
 702     t_rgb rgbd[16] = {
 703         { 1.0, 1.0, 1.0 }, /* white */
 704         { 1.0, 0.0, 0.0 }, /* red */
 705         { 1.0, 1.0, 0.0 }, /* yellow */
 706         { 0.0, 0.0, 1.0 }, /* blue */
 707         { 0.0, 1.0, 0.0 }, /* green */
 708         { 1.0, 0.0, 1.0 }, /* purple */
 709         { 1.0, 0.4, 0.0 }, /* orange */
 710         { 0.0, 1.0, 1.0 }, /* cyan */
 711         { 1.0, 0.4, 0.4 }, /* pink */
 712         { 1.0, 1.0, 0.0 }, /* yellow */
 713         { 0.4, 0.4, 1.0 }, /* lightblue */
 714         { 0.4, 1.0, 0.4 }, /* lightgreen */
 715         { 1.0, 0.4, 1.0 }, /* lightpurple */
 716         { 1.0, 0.7, 0.4 }, /* lightorange */
 717         { 0.4, 1.0, 1.0 }, /* lightcyan */
 718         { 0.0, 0.0, 0.0 }  /* black */
 719     };
 720
 721     int i, n;
 722
 723     *nlevel = std::min(16, *nlevel);
 724     n       = *nlevel;
 725     for (i = 0; (i < n); i++)
 726     {
 727         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%3d\" */,\n", mapper[(i + i0) % NMAP],
 728                 (n <= NMAP) ? ' ' : mapper[(i + i0) / NMAP],
 729                 static_cast<unsigned int>(round(255 * rgbd[i].r)),
 730                 static_cast<unsigned int>(round(255 * rgbd[i].g)),
 731                 static_cast<unsigned int>(round(255 * rgbd[i].b)), i);
 732     }
 733 }
 734
 735
 736 static void write_xpm_map_split(FILE*      out,
 737                                 int        n_x,
 738                                 int        n_y,
 739                                 const int* nlevel_top,
 740                                 real       lo_top,
 741                                 real       hi_top,
 742                                 t_rgb      rlo_top,
 743                                 t_rgb      rhi_top,
 744                                 gmx_bool   bDiscreteColor,
 745                                 int*       nlevel_bot,
 746                                 real       lo_bot,
 747                                 real       hi_bot,
 748                                 t_rgb      rlo_bot,
 749                                 t_rgb      rhi_bot)
 750 {
 751     int ntot;
 752
 753     ntot = *nlevel_top + *nlevel_bot;
 754     if (ntot > NMAP)
 755     {
 756         gmx_fatal(FARGS, "Warning, too many levels (%d) in matrix", ntot);
 757     }
 758
 759     fprintf(out, "static char *gromacs_xpm[] = {\n");
 760     fprintf(out, "\"%d %d   %d %d\",\n", n_x, n_y, ntot, 1);
 761
 762     if (bDiscreteColor)
 763     {
 764         pr_discrete_cmap(out, nlevel_bot, 0);
 765     }
 766     else
 767     {
 768         pr_simple_cmap(out, lo_bot, hi_bot, *nlevel_bot, rlo_bot, rhi_bot, 0);
 769     }
 770
 771     pr_simple_cmap(out, lo_top, hi_top, *nlevel_top, rlo_top, rhi_top, *nlevel_bot);
 772 }
 773
 774
 775 static void write_xpm_map(FILE* out, int n_x, int n_y, int* nlevels, real lo, real hi, t_rgb rlo, t_rgb rhi)
 776 {
 777     int  i, nlo;
 778     real invlevel, r, g, b;
 779
 780     if (*nlevels > NMAP * NMAP)
 781     {
 782         fprintf(stderr, "Warning, too many levels (%d) in matrix, using %d only\n", *nlevels,
 783                 static_cast<int>(NMAP * NMAP));
 784         *nlevels = NMAP * NMAP;
 785     }
 786     else if (*nlevels < 2)
 787     {
 788         fprintf(stderr, "Warning, too few levels (%d) in matrix, using 2 instead\n", *nlevels);
 789         *nlevels = 2;
 790     }
 791
 792     fprintf(out, "static char *gromacs_xpm[] = {\n");
 793     fprintf(out, "\"%d %d   %d %d\",\n", n_x, n_y, *nlevels, (*nlevels <= NMAP) ? 1 : 2);
 794
 795     invlevel = 1.0 / (*nlevels - 1);
 796     for (i = 0; (i < *nlevels); i++)
 797     {
 798         nlo = *nlevels - 1 - i;
 799         r   = (nlo * rlo.r + i * rhi.r) * invlevel;
 800         g   = (nlo * rlo.g + i * rhi.g) * invlevel;
 801         b   = (nlo * rlo.b + i * rhi.b) * invlevel;
 802         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n", mapper[i % NMAP],
 803                 (*nlevels <= NMAP) ? ' ' : mapper[i / NMAP],
 804                 static_cast<unsigned int>(std::round(255 * r)),
 805                 static_cast<unsigned int>(std::round(255 * g)),
 806                 static_cast<unsigned int>(std::round(255 * b)), (nlo * lo + i * hi) * invlevel);
 807     }
 808 }
 809
 810 static void writeXpmAxis(FILE* out, const char* axis, ArrayRef<const real> label)
 811 {
 812     if (label.empty())
 813     {
 814         return;
 815     }
 816     for (gmx::index i = 0; i != ssize(label); ++i)
 817     {
 818         if (i % 80 == 0)
 819         {
 820             if (i)
 821             {
 822                 fprintf(out, "*/\n");
 823             }
 824             fprintf(out, "/* %s-axis:  ", axis);
 825         }
 826         fprintf(out, "%g ", label[i]);
 827     }
 828     fprintf(out, "*/\n");
 829 }
 830
 831 static void write_xpm_data(FILE* out, int n_x, int n_y, real** mat, real lo, real hi, int nlevels)
 832 {
 833     int  i, j, c;
 834     real invlevel;
 835
 836     invlevel = (nlevels - 1) / (hi - lo);
 837     for (j = n_y - 1; (j >= 0); j--)
 838     {
 839         if (j % (1 + n_y / 100) == 0)
 840         {
 841             fprintf(stderr, "%3d%%\b\b\b\b", (100 * (n_y - j)) / n_y);
 842         }
 843         fprintf(out, "\"");
 844         for (i = 0; (i < n_x); i++)
 845         {
 846             c = roundToInt((mat[i][j] - lo) * invlevel);
 847             if (c < 0)
 848             {
 849                 c = 0;
 850             }
 851             if (c >= nlevels)
 852             {
 853                 c = nlevels - 1;
 854             }
 855             if (nlevels <= NMAP)
 856             {
 857                 fprintf(out, "%c", mapper[c]);
 858             }
 859             else
 860             {
 861                 fprintf(out, "%c%c", mapper[c % NMAP], mapper[c / NMAP]);
 862             }
 863         }
 864         if (j > 0)
 865         {
 866             fprintf(out, "\",\n");
 867         }
 868         else
 869         {
 870             fprintf(out, "\"\n");
 871         }
 872     }
 873 }
 874
 875 static void write_xpm_data3(FILE* out, int n_x, int n_y, real** mat, real lo, real mid, real hi, int nlevels)
 876 {
 877     int  i, j, c = 0, nmid;
 878     real invlev_lo, invlev_hi;
 879
 880     nmid      = calc_nmid(nlevels, lo, mid, hi);
 881     invlev_hi = (nlevels - 1 - nmid) / (hi - mid);
 882     invlev_lo = (nmid) / (mid - lo);
 883
 884     for (j = n_y - 1; (j >= 0); j--)
 885     {
 886         if (j % (1 + n_y / 100) == 0)
 887         {
 888             fprintf(stderr, "%3d%%\b\b\b\b", (100 * (n_y - j)) / n_y);
 889         }
 890         fprintf(out, "\"");
 891         for (i = 0; (i < n_x); i++)
 892         {
 893             if (mat[i][j] >= mid)
 894             {
 895                 c = nmid + roundToInt((mat[i][j] - mid) * invlev_hi);
 896             }
 897             else if (mat[i][j] >= lo)
 898             {
 899                 c = roundToInt((mat[i][j] - lo) * invlev_lo);
 900             }
 901             else
 902             {
 903                 c = 0;
 904             }
 905
 906             if (c < 0)
 907             {
 908                 c = 0;
 909             }
 910             if (c >= nlevels)
 911             {
 912                 c = nlevels - 1;
 913             }
 914             if (nlevels <= NMAP)
 915             {
 916                 fprintf(out, "%c", mapper[c]);
 917             }
 918             else
 919             {
 920                 fprintf(out, "%c%c", mapper[c % NMAP], mapper[c / NMAP]);
 921             }
 922         }
 923         if (j > 0)
 924         {
 925             fprintf(out, "\",\n");
 926         }
 927         else
 928         {
 929             fprintf(out, "\"\n");
 930         }
 931     }
 932 }
 933
 934 static void write_xpm_data_split(FILE*  out,
 935                                  int    n_x,
 936                                  int    n_y,
 937                                  real** mat,
 938                                  real   lo_top,
 939                                  real   hi_top,
 940                                  int    nlevel_top,
 941                                  real   lo_bot,
 942                                  real   hi_bot,
 943                                  int    nlevel_bot)
 944 {
 945     int  i, j, c;
 946     real invlev_top, invlev_bot;
 947
 948     invlev_top = (nlevel_top - 1) / (hi_top - lo_top);
 949     invlev_bot = (nlevel_bot - 1) / (hi_bot - lo_bot);
 950
 951     for (j = n_y - 1; (j >= 0); j--)
 952     {
 953         if (j % (1 + n_y / 100) == 0)
 954         {
 955             fprintf(stderr, "%3d%%\b\b\b\b", (100 * (n_y - j)) / n_y);
 956         }
 957         fprintf(out, "\"");
 958         for (i = 0; (i < n_x); i++)
 959         {
 960             if (i < j)
 961             {
 962                 c = nlevel_bot + roundToInt((mat[i][j] - lo_top) * invlev_top);
 963                 if ((c < nlevel_bot) || (c >= nlevel_bot + nlevel_top))
 964                 {
 965                     gmx_fatal(FARGS,
 966                               "Range checking i = %d, j = %d, c = %d, bot = %d, top = %d "
 967                               "matrix[i,j] = %f",
 968                               i, j, c, nlevel_bot, nlevel_top, mat[i][j]);
 969                 }
 970             }
 971             else if (i > j)
 972             {
 973                 c = roundToInt((mat[i][j] - lo_bot) * invlev_bot);
 974                 if ((c < 0) || (c >= nlevel_bot + nlevel_bot))
 975                 {
 976                     gmx_fatal(FARGS,
 977                               "Range checking i = %d, j = %d, c = %d, bot = %d, top = %d "
 978                               "matrix[i,j] = %f",
 979                               i, j, c, nlevel_bot, nlevel_top, mat[i][j]);
 980                 }
 981             }
 982             else
 983             {
 984                 c = nlevel_bot;
 985             }
 986
 987             fprintf(out, "%c", mapper[c]);
 988         }
 989         if (j > 0)
 990         {
 991             fprintf(out, "\",\n");
 992         }
 993         else
 994         {
 995             fprintf(out, "\"\n");
 996         }
 997     }
 998 }
 999
1000 void write_xpm_m(FILE* out, t_matrix m)
1001 {
1002     gmx_bool  bOneChar;
1003     t_xpmelmt c;
1004
1005     bOneChar = (m.map[0].code.c2 == 0);
1006     write_xpm_header(out, m.title, m.legend, m.label_x, m.label_y, m.bDiscrete);
1007     fprintf(out, "static char *gromacs_xpm[] = {\n");
1008     fprintf(out, "\"%d %d   %zu %d\",\n", m.nx, m.ny, m.map.size(), bOneChar ? 1 : 2);
1009     for (const auto& map : m.map)
1010     {
1011         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%s\" */,\n", map.code.c1,
1012                 bOneChar ? ' ' : map.code.c2, static_cast<unsigned int>(round(map.rgb.r * 255)),
1013                 static_cast<unsigned int>(round(map.rgb.g * 255)),
1014                 static_cast<unsigned int>(round(map.rgb.b * 255)), map.desc);
1015     }
1016     writeXpmAxis(out, "x", m.axis_x);
1017     writeXpmAxis(out, "y", m.axis_y);
1018     for (int j = m.ny - 1; (j >= 0); j--)
1019     {
1020         if (j % (1 + m.ny / 100) == 0)
1021         {
1022             fprintf(stderr, "%3d%%\b\b\b\b", (100 * (m.ny - j)) / m.ny);
1023         }
1024         fprintf(out, "\"");
1025         if (bOneChar)
1026         {
1027             for (int i = 0; (i < m.nx); i++)
1028             {
1029                 fprintf(out, "%c", m.map[m.matrix(i, j)].code.c1);
1030             }
1031         }
1032         else
1033         {
1034             for (int i = 0; (i < m.nx); i++)
1035             {
1036                 c = m.map[m.matrix(i, j)].code;
1037                 fprintf(out, "%c%c", c.c1, c.c2);
1038             }
1039         }
1040         if (j > 0)
1041         {
1042             fprintf(out, "\",\n");
1043         }
1044         else
1045         {
1046             fprintf(out, "\"\n");
1047         }
1048     }
1049 }
1050
1051 void write_xpm3(FILE*              out,
1052                 unsigned int       flags,
1053                 const std::string& title,
1054                 const std::string& legend,
1055                 const std::string& label_x,
1056                 const std::string& label_y,
1057                 int                n_x,
1058                 int                n_y,
1059                 real               axis_x[],
1060                 real               axis_y[],
1061                 real*              mat[],
1062                 real               lo,
1063                 real               mid,
1064                 real               hi,
1065                 t_rgb              rlo,
1066                 t_rgb              rmid,
1067                 t_rgb              rhi,
1068                 int*               nlevels)
1069 {
1070     /* See write_xpm.
1071      * Writes a colormap varying as rlo -> rmid -> rhi.
1072      */
1073
1074     if (hi <= lo)
1075     {
1076         gmx_fatal(FARGS, "hi (%g) <= lo (%g)", hi, lo);
1077     }
1078
1079     write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1080     write_xpm_map3(out, n_x, n_y, nlevels, lo, mid, hi, rlo, rmid, rhi);
1081     writeXpmAxis(out, "x", ArrayRef<real>(axis_x, axis_x + n_x + ((flags & MAT_SPATIAL_X) != 0U ? 1 : 0)));
1082     writeXpmAxis(out, "y", ArrayRef<real>(axis_y, axis_y + n_y + ((flags & MAT_SPATIAL_Y) != 0U ? 1 : 0)));
1083     write_xpm_data3(out, n_x, n_y, mat, lo, mid, hi, *nlevels);
1084 }
1085
1086 void write_xpm_split(FILE*              out,
1087                      unsigned int       flags,
1088                      const std::string& title,
1089                      const std::string& legend,
1090                      const std::string& label_x,
1091                      const std::string& label_y,
1092                      int                n_x,
1093                      int                n_y,
1094                      real               axis_x[],
1095                      real               axis_y[],
1096                      real*              mat[],
1097                      real               lo_top,
1098                      real               hi_top,
1099                      int*               nlevel_top,
1100                      t_rgb              rlo_top,
1101                      t_rgb              rhi_top,
1102                      real               lo_bot,
1103                      real               hi_bot,
1104                      int*               nlevel_bot,
1105                      gmx_bool           bDiscreteColor,
1106                      t_rgb              rlo_bot,
1107                      t_rgb              rhi_bot)
1108 {
1109     /* See write_xpm.
1110      * Writes a colormap varying as rlo -> rmid -> rhi.
1111      */
1112
1113     if (hi_top <= lo_top)
1114     {
1115         gmx_fatal(FARGS, "hi_top (%g) <= lo_top (%g)", hi_top, lo_top);
1116     }
1117     if (hi_bot <= lo_bot)
1118     {
1119         gmx_fatal(FARGS, "hi_bot (%g) <= lo_bot (%g)", hi_bot, lo_bot);
1120     }
1121     if (bDiscreteColor && (*nlevel_bot >= 16))
1122     {
1123         gmx_impl("Can not plot more than 16 discrete colors");
1124     }
1125
1126     write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1127     write_xpm_map_split(out, n_x, n_y, nlevel_top, lo_top, hi_top, rlo_top, rhi_top, bDiscreteColor,
1128                         nlevel_bot, lo_bot, hi_bot, rlo_bot, rhi_bot);
1129     writeXpmAxis(out, "x", ArrayRef<real>(axis_x, axis_x + n_x + ((flags & MAT_SPATIAL_X) != 0U ? 1 : 0)));
1130     writeXpmAxis(out, "y", ArrayRef<real>(axis_y, axis_y + n_y + ((flags & MAT_SPATIAL_Y) != 0U ? 1 : 0)));
1131     write_xpm_data_split(out, n_x, n_y, mat, lo_top, hi_top, *nlevel_top, lo_bot, hi_bot, *nlevel_bot);
1132 }
1133
1134 void write_xpm(FILE*              out,
1135                unsigned int       flags,
1136                const std::string& title,
1137                const std::string& legend,
1138                const std::string& label_x,
1139                const std::string& label_y,
1140                int                n_x,
1141                int                n_y,
1142                real               axis_x[],
1143                real               axis_y[],
1144                real*              mat[],
1145                real               lo,
1146                real               hi,
1147                t_rgb              rlo,
1148                t_rgb              rhi,
1149                int*               nlevels)
1150 {
1151     /* out        xpm file
1152      * title      matrix title
1153      * legend     label for the continuous legend
1154      * label_x    label for the x-axis
1155      * label_y    label for the y-axis
1156      * n_x, n_y   size of the matrix
1157      * axis_x[]   the x-ticklabels
1158      * axis_y[]   the y-ticklables
1159      * *matrix[]  element x,y is matrix[x][y]
1160      * lo         output lower than lo is set to lo
1161      * hi         output higher than hi is set to hi
1162      * rlo        rgb value for level lo
1163      * rhi        rgb value for level hi
1164      * nlevels    number of color levels for the output
1165      */
1166
1167     if (hi <= lo)
1168     {
1169         gmx_fatal(FARGS, "hi (%f) <= lo (%f)", hi, lo);
1170     }
1171
1172     write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1173     write_xpm_map(out, n_x, n_y, nlevels, lo, hi, rlo, rhi);
1174     writeXpmAxis(out, "x", ArrayRef<real>(axis_x, axis_x + n_x + ((flags & MAT_SPATIAL_X) != 0U ? 1 : 0)));
1175     writeXpmAxis(out, "y", ArrayRef<real>(axis_y, axis_y + n_y + ((flags & MAT_SPATIAL_Y) != 0U ? 1 : 0)));
1176     write_xpm_data(out, n_x, n_y, mat, lo, hi, *nlevels);
1177 }