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