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