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