src/gromacs/gmxana/gmx_xpm2ps.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,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,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 <cassert>
  41 #include <cmath>
  42 #include <cstdio>
  43 #include <cstdlib>
  44 #include <cstring>
  45
  46 #include <algorithm>
  47 #include <numeric>
  48 #include <string>
  49
  50 #include "gromacs/commandline/pargs.h"
  51 #include "gromacs/commandline/viewit.h"
  52 #include "gromacs/fileio/matio.h"
  53 #include "gromacs/fileio/readinp.h"
  54 #include "gromacs/fileio/trxio.h"
  55 #include "gromacs/fileio/warninp.h"
  56 #include "gromacs/fileio/writeps.h"
  57 #include "gromacs/gmxana/gmx_ana.h"
  58 #include "gromacs/utility/arrayref.h"
  59 #include "gromacs/utility/arraysize.h"
  60 #include "gromacs/utility/cstringutil.h"
  61 #include "gromacs/utility/exceptions.h"
  62 #include "gromacs/utility/fatalerror.h"
  63 #include "gromacs/utility/filestream.h"
  64 #include "gromacs/utility/futil.h"
  65 #include "gromacs/utility/gmxassert.h"
  66 #include "gromacs/utility/smalloc.h"
  67 #include "gromacs/utility/stringutil.h"
  68
  69 #define FUDGE 1.2
  70 #define DDD 2
  71
  72 typedef struct
  73 {
  74     real     major;
  75     real     minor;
  76     real     offset;
  77     gmx_bool first;
  78     int      lineatzero;
  79     real     majorticklen;
  80     real     minorticklen;
  81     char     label[STRLEN];
  82     real     fontsize;
  83     char     font[STRLEN];
  84     real     tickfontsize;
  85     char     tickfont[STRLEN];
  86 } t_axisdef;
  87
  88 typedef struct
  89 {
  90     int       bw;
  91     real      linewidth;
  92     real      xoffs, yoffs;
  93     gmx_bool  bTitle;
  94     gmx_bool  bTitleOnce;
  95     gmx_bool  bYonce;
  96     real      titfontsize;
  97     char      titfont[STRLEN];
  98     gmx_bool  legend;
  99     real      legfontsize;
 100     char      legfont[STRLEN];
 101     char      leglabel[STRLEN];
 102     char      leg2label[STRLEN];
 103     real      xboxsize;
 104     real      yboxsize;
 105     real      boxspacing;
 106     real      boxlinewidth;
 107     real      ticklinewidth;
 108     real      zerolinewidth;
 109     t_axisdef X, Y;
 110 } t_psrec;
 111
 112 /* MUST correspond to char *legend[] in main() */
 113 enum
 114 {
 115     elSel,
 116     elBoth,
 117     elFirst,
 118     elSecond,
 119     elNone,
 120     elNR
 121 };
 122
 123 /* MUST correspond to char *combine[] in main() */
 124 enum
 125 {
 126     ecSel,
 127     ecHalves,
 128     ecAdd,
 129     ecSub,
 130     ecMult,
 131     ecDiv,
 132     ecNR
 133 };
 134
 135 static void get_params(const char* mpin, const char* mpout, t_psrec* psr)
 136 {
 137     static const char* gmx_bools[BOOL_NR + 1] = { "no", "yes", nullptr };
 138     /* this must correspond to t_rgb *linecolors[] below */
 139     static const char*     colors[] = { "none", "black", "white", nullptr };
 140     warninp_t              wi;
 141     std::vector<t_inpfile> inp;
 142
 143     wi = init_warning(FALSE, 0);
 144
 145     if (mpin != nullptr)
 146     {
 147         std::string        libmpin = gmx::findLibraryFile(mpin);
 148         gmx::TextInputFile stream(libmpin);
 149         inp = read_inpfile(&stream, libmpin.c_str(), wi);
 150     }
 151     else
 152     {
 153         inp.clear();
 154     }
 155
 156     psr->bw        = get_eenum(&inp, "black&white", gmx_bools);
 157     psr->linewidth = get_ereal(&inp, "linewidth", 1.0, wi);
 158     setStringEntry(&inp, "titlefont", psr->titfont, "Helvetica");
 159     psr->titfontsize = get_ereal(&inp, "titlefontsize", 20.0, wi);
 160     psr->legend      = (get_eenum(&inp, "legend", gmx_bools) != 0);
 161     setStringEntry(&inp, "legendfont", psr->legfont, psr->titfont);
 162     setStringEntry(&inp, "legendlabel", psr->leglabel, "");
 163     setStringEntry(&inp, "legend2label", psr->leg2label, psr->leglabel);
 164     psr->legfontsize    = get_ereal(&inp, "legendfontsize", 14.0, wi);
 165     psr->xboxsize       = get_ereal(&inp, "xbox", 0.0, wi);
 166     psr->yboxsize       = get_ereal(&inp, "ybox", 0.0, wi);
 167     psr->boxspacing     = get_ereal(&inp, "matrixspacing", 20.0, wi);
 168     psr->xoffs          = get_ereal(&inp, "xoffset", 0.0, wi);
 169     psr->yoffs          = get_ereal(&inp, "yoffset", psr->xoffs, wi);
 170     psr->boxlinewidth   = get_ereal(&inp, "boxlinewidth", psr->linewidth, wi);
 171     psr->ticklinewidth  = get_ereal(&inp, "ticklinewidth", psr->linewidth, wi);
 172     psr->zerolinewidth  = get_ereal(&inp, "zerolinewidth", psr->ticklinewidth, wi);
 173     psr->X.lineatzero   = get_eenum(&inp, "x-lineat0value", colors);
 174     psr->X.major        = get_ereal(&inp, "x-major", -1, wi);
 175     psr->X.minor        = get_ereal(&inp, "x-minor", -1, wi);
 176     psr->X.offset       = get_ereal(&inp, "x-firstmajor", 0.0, wi);
 177     psr->X.first        = (get_eenum(&inp, "x-majorat0", gmx_bools) != 0);
 178     psr->X.majorticklen = get_ereal(&inp, "x-majorticklen", 8.0, wi);
 179     psr->X.minorticklen = get_ereal(&inp, "x-minorticklen", 4.0, wi);
 180     setStringEntry(&inp, "x-label", psr->X.label, "");
 181     psr->X.fontsize = get_ereal(&inp, "x-fontsize", 16.0, wi);
 182     setStringEntry(&inp, "x-font", psr->X.font, psr->titfont);
 183     psr->X.tickfontsize = get_ereal(&inp, "x-tickfontsize", 10.0, wi);
 184     setStringEntry(&inp, "x-tickfont", psr->X.tickfont, psr->X.font);
 185     psr->Y.lineatzero   = get_eenum(&inp, "y-lineat0value", colors);
 186     psr->Y.major        = get_ereal(&inp, "y-major", psr->X.major, wi);
 187     psr->Y.minor        = get_ereal(&inp, "y-minor", psr->X.minor, wi);
 188     psr->Y.offset       = get_ereal(&inp, "y-firstmajor", psr->X.offset, wi);
 189     psr->Y.first        = (get_eenum(&inp, "y-majorat0", gmx_bools) != 0);
 190     psr->Y.majorticklen = get_ereal(&inp, "y-majorticklen", psr->X.majorticklen, wi);
 191     psr->Y.minorticklen = get_ereal(&inp, "y-minorticklen", psr->X.minorticklen, wi);
 192     setStringEntry(&inp, "y-label", psr->Y.label, psr->X.label);
 193     psr->Y.fontsize = get_ereal(&inp, "y-fontsize", psr->X.fontsize, wi);
 194     setStringEntry(&inp, "y-font", psr->Y.font, psr->X.font);
 195     psr->Y.tickfontsize = get_ereal(&inp, "y-tickfontsize", psr->X.tickfontsize, wi);
 196     setStringEntry(&inp, "y-tickfont", psr->Y.tickfont, psr->Y.font);
 197
 198     check_warning_error(wi, FARGS);
 199
 200     if (mpout != nullptr)
 201     {
 202         gmx::TextOutputFile stream(mpout);
 203         write_inpfile(&stream, mpout, &inp, TRUE, WriteMdpHeader::yes, wi);
 204         stream.close();
 205     }
 206
 207     done_warning(wi, FARGS);
 208 }
 209
 210 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
 211 static t_rgb black = { 0, 0, 0 };
 212 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
 213 static t_rgb white = { 1, 1, 1 };
 214 #define BLACK (&black)
 215 /* this must correspond to *colors[] in get_params */
 216 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
 217 static t_rgb* linecolors[] = { nullptr, &black, &white, nullptr };
 218
 219 static void leg_discrete(t_psdata*                      ps,
 220                          real                           x0,
 221                          real                           y0,
 222                          const std::string&             label,
 223                          real                           fontsize,
 224                          char*                          font,
 225                          gmx::ArrayRef<const t_mapping> map)
 226 {
 227     real yhh;
 228     real boxhh;
 229
 230     boxhh = fontsize + DDD;
 231     /* LANDSCAPE */
 232     ps_rgb(ps, BLACK);
 233     ps_strfont(ps, font, fontsize);
 234     yhh = y0 + fontsize + 3 * DDD;
 235     if (!label.empty())
 236     {
 237         ps_ctext(ps, x0, yhh, label, eXLeft);
 238     }
 239     ps_moveto(ps, x0, y0);
 240     for (const auto& m : map)
 241     {
 242         ps_setorigin(ps);
 243         ps_rgb(ps, &m.rgb);
 244         ps_fillbox(ps, DDD, DDD, DDD + fontsize, boxhh - DDD);
 245         ps_rgb(ps, BLACK);
 246         ps_box(ps, DDD, DDD, DDD + fontsize, boxhh - DDD);
 247         ps_ctext(ps, boxhh + 2 * DDD, fontsize / 3, m.desc, eXLeft);
 248         ps_unsetorigin(ps);
 249         ps_moverel(ps, DDD, -fontsize / 3);
 250     }
 251 }
 252
 253 static void leg_continuous(t_psdata*                      ps,
 254                            real                           x0,
 255                            real                           x,
 256                            real                           y0,
 257                            const std::string&             label,
 258                            real                           fontsize,
 259                            char*                          font,
 260                            gmx::ArrayRef<const t_mapping> map,
 261                            int                            mapoffset)
 262 {
 263     real       xx0;
 264     real       yhh, boxxh, boxyh;
 265     gmx::index mapIndex = gmx::ssize(map) - mapoffset;
 266
 267     boxyh = fontsize;
 268     if (x < 8 * fontsize)
 269     {
 270         x = 8 * fontsize;
 271     }
 272     boxxh = x / mapIndex;
 273     if (boxxh > fontsize)
 274     {
 275         boxxh = fontsize;
 276     }
 277
 278     GMX_RELEASE_ASSERT(!map.empty(), "NULL map array provided to leg_continuous()");
 279
 280     /* LANDSCAPE */
 281     xx0 = x0 - (mapIndex * boxxh) / 2.0;
 282
 283     for (gmx::index i = 0; (i < mapIndex); i++)
 284     {
 285         ps_rgb(ps, &(map[i + mapoffset].rgb));
 286         ps_fillbox(ps, xx0 + i * boxxh, y0, xx0 + (i + 1) * boxxh, y0 + boxyh);
 287     }
 288     ps_strfont(ps, font, fontsize);
 289     ps_rgb(ps, BLACK);
 290     ps_box(ps, xx0, y0, xx0 + mapIndex * boxxh, y0 + boxyh);
 291
 292     yhh = y0 + boxyh + 3 * DDD;
 293     ps_ctext(ps, xx0 + boxxh / 2, yhh, map[0].desc, eXCenter);
 294     if (!label.empty())
 295     {
 296         ps_ctext(ps, x0, yhh, label, eXCenter);
 297     }
 298     ps_ctext(ps, xx0 + (mapIndex * boxxh) - boxxh / 2, yhh, map[map.size() - 1].desc, eXCenter);
 299 }
 300
 301 static void leg_bicontinuous(t_psdata*                      ps,
 302                              real                           x0,
 303                              real                           x,
 304                              real                           y0,
 305                              const std::string&             label1,
 306                              const std::string&             label2,
 307                              real                           fontsize,
 308                              char*                          font,
 309                              gmx::ArrayRef<const t_mapping> map1,
 310                              gmx::ArrayRef<const t_mapping> map2)
 311 {
 312     real xx1, xx2, x1, x2;
 313
 314     x1  = x / (map1.size() + map2.size()) * map1.size(); /* width of legend 1 */
 315     x2  = x / (map1.size() + map2.size()) * map2.size(); /* width of legend 2 */
 316     xx1 = x0 - (x2 / 2.0) - fontsize;                    /* center of legend 1 */
 317     xx2 = x0 + (x1 / 2.0) + fontsize;                    /* center of legend 2 */
 318     x1 -= fontsize / 2;                                  /* adjust width */
 319     x2 -= fontsize / 2;                                  /* adjust width */
 320     leg_continuous(ps, xx1, x1, y0, label1, fontsize, font, map1, 0);
 321     leg_continuous(ps, xx2, x2, y0, label2, fontsize, font, map2, 0);
 322 }
 323
 324 static real box_height(const t_matrix& mat, t_psrec* psr)
 325 {
 326     return mat.ny * psr->yboxsize;
 327 }
 328
 329 static real box_dh(t_psrec* psr)
 330 {
 331     return psr->boxspacing;
 332 }
 333
 334 static real box_dh_top(gmx_bool bOnce, t_psrec* psr)
 335 {
 336     real dh;
 337
 338     if (psr->bTitle || (psr->bTitleOnce && bOnce))
 339     {
 340         dh = 2 * psr->titfontsize;
 341     }
 342     else
 343     {
 344         dh = 0;
 345     }
 346
 347     return dh;
 348 }
 349
 350 static gmx_bool box_do_all_x_maj_ticks(t_psrec* psr)
 351 {
 352     return (psr->boxspacing > (1.5 * psr->X.majorticklen));
 353 }
 354
 355 static gmx_bool box_do_all_x_min_ticks(t_psrec* psr)
 356 {
 357     return (psr->boxspacing > (1.5 * psr->X.minorticklen));
 358 }
 359
 360 static void draw_boxes(t_psdata* ps, real x0, real y0, real w, gmx::ArrayRef<t_matrix> mat, t_psrec* psr)
 361 {
 362     char   buf[128];
 363     real   xxx;
 364     char **xtick, **ytick;
 365     real   xx, yy, dy, xx00, yy00, offset_x, offset_y;
 366     int    x, ntx, nty;
 367     size_t strlength;
 368
 369     /* Only necessary when there will be no y-labels */
 370     strlength = 0;
 371
 372     /* Draw the box */
 373     ps_rgb(ps, BLACK);
 374     ps_linewidth(ps, static_cast<int>(psr->boxlinewidth));
 375     yy00 = y0;
 376     for (auto m = mat.begin(); m != mat.end(); ++m)
 377     {
 378         dy = box_height(*m, psr);
 379         ps_box(ps, x0 - 1, yy00 - 1, x0 + w + 1, yy00 + dy + 1);
 380         yy00 += dy + box_dh(psr) + box_dh_top(m + 1 == mat.end(), psr);
 381     }
 382
 383     /* Draw the ticks on the axes */
 384     ps_linewidth(ps, static_cast<int>(psr->ticklinewidth));
 385     xx00         = x0 - 1;
 386     yy00         = y0 - 1;
 387     auto halfway = mat.begin() + (mat.size() / 2);
 388     for (auto m = mat.begin(); m != mat.end(); ++m)
 389     {
 390         if (m->flags & MAT_SPATIAL_X)
 391         {
 392             ntx      = m->nx + 1;
 393             offset_x = 0.1;
 394         }
 395         else
 396         {
 397             ntx      = m->nx;
 398             offset_x = 0.6;
 399         }
 400         if (m->flags & MAT_SPATIAL_Y)
 401         {
 402             nty      = m->ny + 1;
 403             offset_y = 0.1;
 404         }
 405         else
 406         {
 407             nty      = m->ny;
 408             offset_y = 0.6;
 409         }
 410         snew(xtick, ntx);
 411         for (int j = 0; (j < ntx); j++)
 412         {
 413             sprintf(buf, "%g", m->axis_x[j]);
 414             xtick[j] = gmx_strdup(buf);
 415         }
 416         ps_strfont(ps, psr->X.tickfont, psr->X.tickfontsize);
 417         for (x = 0; (x < ntx); x++)
 418         {
 419             xx = xx00 + (x + offset_x) * psr->xboxsize;
 420             if ((bRmod(m->axis_x[x], psr->X.offset, psr->X.major) || (psr->X.first && (x == 0)))
 421                 && (m == mat.begin() || box_do_all_x_maj_ticks(psr)))
 422             {
 423                 /* Longer tick marks */
 424                 ps_line(ps, xx, yy00, xx, yy00 - psr->X.majorticklen);
 425                 /* Plot label on lowest graph only */
 426                 if (m == mat.begin())
 427                 {
 428                     ps_ctext(ps, xx, yy00 - DDD - psr->X.majorticklen - psr->X.tickfontsize * 0.8, xtick[x], eXCenter);
 429                 }
 430             }
 431             else if (bRmod(m->axis_x[x], psr->X.offset, psr->X.minor)
 432                      && ((m == mat.begin()) || box_do_all_x_min_ticks(psr)))
 433             {
 434                 /* Shorter tick marks */
 435                 ps_line(ps, xx, yy00, xx, yy00 - psr->X.minorticklen);
 436             }
 437             else if (bRmod(m->axis_x[x], psr->X.offset, psr->X.major))
 438             {
 439                 /* Even shorter marks, only each X.major */
 440                 ps_line(ps, xx, yy00, xx, yy00 - (psr->boxspacing / 2));
 441             }
 442         }
 443         ps_strfont(ps, psr->Y.tickfont, psr->Y.tickfontsize);
 444         snew(ytick, nty);
 445         for (int j = 0; (j < nty); j++)
 446         {
 447             sprintf(buf, "%g", m->axis_y[j]);
 448             ytick[j] = gmx_strdup(buf);
 449         }
 450
 451         for (int y = 0; (y < nty); y++)
 452         {
 453             yy = yy00 + (y + offset_y) * psr->yboxsize;
 454             if (bRmod(m->axis_y[y], psr->Y.offset, psr->Y.major) || (psr->Y.first && (y == 0)))
 455             {
 456                 /* Major ticks */
 457                 strlength = std::max(strlength, std::strlen(ytick[y]));
 458                 ps_line(ps, xx00, yy, xx00 - psr->Y.majorticklen, yy);
 459                 ps_ctext(ps, xx00 - psr->Y.majorticklen - DDD, yy - psr->Y.tickfontsize / 3.0, ytick[y], eXRight);
 460             }
 461             else if (bRmod(m->axis_y[y], psr->Y.offset, psr->Y.minor))
 462             {
 463                 /* Minor ticks */
 464                 ps_line(ps, xx00, yy, xx00 - psr->Y.minorticklen, yy);
 465             }
 466         }
 467         sfree(xtick);
 468         sfree(ytick);
 469
 470         /* Label on Y-axis */
 471         if (!psr->bYonce || m == halfway)
 472         {
 473             std::string mylab;
 474             if (strlen(psr->Y.label) > 0)
 475             {
 476                 mylab = psr->Y.label;
 477             }
 478             else
 479             {
 480                 mylab = m->label_y;
 481             }
 482             if (!mylab.empty())
 483             {
 484                 ps_strfont(ps, psr->Y.font, psr->Y.fontsize);
 485                 ps_flip(ps, TRUE);
 486                 xxx = x0 - psr->X.majorticklen - psr->X.tickfontsize * strlength - DDD;
 487                 ps_ctext(ps, yy00 + box_height(*m, psr) / 2.0, 612.5 - xxx, mylab, eXCenter);
 488                 ps_flip(ps, FALSE);
 489             }
 490         }
 491
 492         yy00 += box_height(*m, psr) + box_dh(psr) + box_dh_top(m + 1 == mat.end(), psr);
 493     }
 494     /* Label on X-axis */
 495     std::string mylab;
 496     if (strlen(psr->X.label) > 0)
 497     {
 498         mylab = psr->X.label;
 499     }
 500     else
 501     {
 502         mylab = mat[0].label_x;
 503     }
 504     if (!mylab.empty())
 505     {
 506         ps_strfont(ps, psr->X.font, psr->X.fontsize);
 507         ps_ctext(ps,
 508                  x0 + w / 2,
 509                  y0 - DDD - psr->X.majorticklen - psr->X.tickfontsize * FUDGE - psr->X.fontsize,
 510                  mylab,
 511                  eXCenter);
 512     }
 513 }
 514
 515 static void draw_zerolines(t_psdata* out, real x0, real y0, real w, gmx::ArrayRef<t_matrix> mat, t_psrec* psr)
 516 {
 517     real xx, yy, dy, xx00, yy00;
 518     int  x, y;
 519
 520     xx00 = x0 - 1.5;
 521     yy00 = y0 - 1.5;
 522     ps_linewidth(out, static_cast<int>(psr->zerolinewidth));
 523     for (auto m = mat.begin(); m != mat.end(); ++m)
 524     {
 525         dy = box_height(*m, psr);
 526         /* m->axis_x and _y were already set by draw_boxes */
 527         if (psr->X.lineatzero)
 528         {
 529             ps_rgb(out, linecolors[psr->X.lineatzero]);
 530             for (x = 0; (x < m->nx); x++)
 531             {
 532                 xx = xx00 + (x + 0.7) * psr->xboxsize;
 533                 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
 534                 if (x != 0 && x < m->nx - 1
 535                     && std::abs(m->axis_x[x]) < 0.1 * std::abs(m->axis_x[x + 1] - m->axis_x[x]))
 536                 {
 537                     ps_line(out, xx, yy00, xx, yy00 + dy + 2);
 538                 }
 539             }
 540         }
 541         if (psr->Y.lineatzero)
 542         {
 543             ps_rgb(out, linecolors[psr->Y.lineatzero]);
 544             for (y = 0; (y < m->ny); y++)
 545             {
 546                 yy = yy00 + (y + 0.7) * psr->yboxsize;
 547                 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
 548                 if (y != 0 && y < m->ny - 1
 549                     && std::abs(m->axis_y[y]) < 0.1 * std::abs(m->axis_y[y + 1] - m->axis_y[y]))
 550                 {
 551                     ps_line(out, xx00, yy, xx00 + w + 2, yy);
 552                 }
 553             }
 554         }
 555         yy00 += box_height(*m, psr) + box_dh(psr) + box_dh_top(m + 1 == mat.end(), psr);
 556     }
 557 }
 558
 559 static void box_dim(gmx::ArrayRef<t_matrix> mat,
 560                     gmx::ArrayRef<t_matrix> mat2,
 561                     t_psrec*                psr,
 562                     int                     elegend,
 563                     gmx_bool                bFrame,
 564                     real*                   w,
 565                     real*                   h,
 566                     real*                   dw,
 567                     real*                   dh)
 568 {
 569     int  maxytick;
 570     real ww, hh, dww, dhh;
 571
 572     hh = dww = dhh = 0;
 573     maxytick       = 0;
 574
 575     ww = 0;
 576     for (const auto& m : mat)
 577     {
 578         ww = std::max(ww, m.nx * psr->xboxsize);
 579         hh += box_height(m, psr);
 580         maxytick = std::max(maxytick, m.nx);
 581     }
 582     if (bFrame)
 583     {
 584         if (mat[0].label_y[0])
 585         {
 586             dww += 2.0 * (psr->Y.fontsize + DDD);
 587         }
 588         if (psr->Y.major > 0)
 589         {
 590             dww += psr->Y.majorticklen + DDD
 591                    + psr->Y.tickfontsize * (std::log(static_cast<real>(maxytick)) / std::log(10.0));
 592         }
 593         else if (psr->Y.minor > 0)
 594         {
 595             dww += psr->Y.minorticklen;
 596         }
 597
 598         if (mat[0].label_x[0])
 599         {
 600             dhh += psr->X.fontsize + 2 * DDD;
 601         }
 602         if (/* fool emacs auto-indent */
 603             (elegend == elBoth && (mat[0].legend[0] || (!mat2.empty() && mat2[0].legend[0])))
 604             || (elegend == elFirst && mat[0].legend[0])
 605             || (elegend == elSecond && (!mat2.empty() && mat2[0].legend[0])))
 606         {
 607             dhh += 2 * (psr->legfontsize * FUDGE + 2 * DDD);
 608         }
 609         else
 610         {
 611             dhh += psr->legfontsize * FUDGE + 2 * DDD;
 612         }
 613         if (psr->X.major > 0)
 614         {
 615             dhh += psr->X.tickfontsize * FUDGE + 2 * DDD + psr->X.majorticklen;
 616         }
 617         else if (psr->X.minor > 0)
 618         {
 619             dhh += psr->X.minorticklen;
 620         }
 621
 622         hh += (mat.size() - 1) * box_dh(psr);
 623         hh += box_dh_top(TRUE, psr);
 624         if (mat.size() > 1)
 625         {
 626             hh += (mat.size() - 1) * box_dh_top(FALSE, psr);
 627         }
 628     }
 629     *w  = ww;
 630     *h  = hh;
 631     *dw = dww;
 632     *dh = dhh;
 633 }
 634
 635 static std::vector<t_mapping> add_maps(gmx::ArrayRef<t_mapping> map1, gmx::ArrayRef<t_mapping> map2)
 636 {
 637     static char mapper[] =
 638             "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<"
 639             ".>/?";
 640     std::vector<t_mapping> map(map1.size() + map2.size());
 641
 642     size_t nsymbols = std::strlen(mapper);
 643     if (map.size() > nsymbols * nsymbols)
 644     {
 645         gmx_fatal(FARGS, "Not enough symbols to merge the two colormaps\n");
 646     }
 647     printf("Combining colormaps of %zu and %zu elements into one of %zu elements\n",
 648            map1.size(),
 649            map2.size(),
 650            map.size());
 651     gmx::index k = 0;
 652     for (gmx::index j = 0; j < gmx::ssize(map1) && k < gmx::ssize(map); ++j, ++k)
 653     {
 654         map[k].code.c1 = mapper[k % nsymbols];
 655         if (map.size() > nsymbols)
 656         {
 657             map[k].code.c2 = mapper[k / nsymbols];
 658         }
 659         map[k].rgb.r = map1[j].rgb.r;
 660         map[k].rgb.g = map1[j].rgb.g;
 661         map[k].rgb.b = map1[j].rgb.b;
 662         map[k].desc  = map1[j].desc;
 663     }
 664     for (gmx::index j = 0; j < gmx::ssize(map2) && k < gmx::ssize(map); ++j, ++k)
 665     {
 666         map[k].code.c1 = mapper[k % nsymbols];
 667         if (map.size() > nsymbols)
 668         {
 669             map[k].code.c2 = mapper[k / nsymbols];
 670         }
 671         map[k].rgb.r = map2[j].rgb.r;
 672         map[k].rgb.g = map2[j].rgb.g;
 673         map[k].rgb.b = map2[j].rgb.b;
 674         map[k].desc  = map2[j].desc;
 675     }
 676
 677     return map;
 678 }
 679
 680 static bool operator==(const t_mapping& lhs, const t_mapping& rhs)
 681 {
 682     return (lhs.rgb.r == rhs.rgb.r && lhs.rgb.g == rhs.rgb.g && lhs.rgb.b == rhs.rgb.b);
 683 }
 684
 685 static void xpm_mat(const char*             outf,
 686                     gmx::ArrayRef<t_matrix> mat,
 687                     gmx::ArrayRef<t_matrix> mat2,
 688                     gmx_bool                bDiag,
 689                     gmx_bool                bFirstDiag)
 690 {
 691     FILE* out;
 692     int   x, y, col;
 693
 694     out = gmx_ffopen(outf, "w");
 695
 696     GMX_RELEASE_ASSERT(mat.size() == mat2.size(),
 697                        "Combined matrix write requires matrices of the same size");
 698     for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
 699     {
 700         // Color maps that differ only in RGB value are considered different
 701         if (mat2.empty() || std::equal(mat[i].map.begin(), mat[i].map.end(), mat2[i].map.begin()))
 702         {
 703             write_xpm_m(out, mat[0]);
 704         }
 705         else
 706         {
 707             auto map = add_maps(mat[i].map, mat2[i].map);
 708             for (x = 0; (x < mat[i].nx); x++)
 709             {
 710                 for (y = 0; (y < mat[i].nx); y++)
 711                 {
 712                     if ((x < y) || ((x == y) && bFirstDiag)) /* upper left  -> map1 */
 713                     {
 714                         col = mat[i].matrix(x, y);
 715                     }
 716                     else /* lower right -> map2 */
 717                     {
 718                         col = mat[i].map.size() + mat[i].matrix(x, y);
 719                     }
 720                     if ((bDiag) || (x != y))
 721                     {
 722                         mat[i].matrix(x, y) = col;
 723                     }
 724                     else
 725                     {
 726                         mat[i].matrix(x, y) = 0;
 727                     }
 728                 }
 729             }
 730             mat[i].map = map;
 731             if (mat[i].title != mat2[i].title)
 732             {
 733                 mat[i].title += " / " + mat2[i].title;
 734             }
 735             if (mat[i].legend != mat2[i].legend)
 736             {
 737                 mat[i].legend += " / " + mat2[i].legend;
 738             }
 739             write_xpm_m(out, mat[i]);
 740         }
 741     }
 742     gmx_ffclose(out);
 743 }
 744
 745 static void tick_spacing(int n, real axis[], real offset, char axisnm, real* major, real* minor)
 746 {
 747     real     space;
 748     gmx_bool bTryAgain;
 749     int      i, j, t, f = 0, ten;
 750 #define NFACT 4
 751     real major_fact[NFACT] = { 5, 4, 2, 1 };
 752     real minor_fact[NFACT] = { 5, 4, 4, 5 };
 753
 754     /* start with interval between 10 matrix points: */
 755     space = std::max(10 * axis[1] - axis[0], axis[std::min(10, n - 1)] - axis[0]);
 756     /* get power of 10 */
 757     ten       = static_cast<int>(std::ceil(std::log(space) / std::log(10.0)) - 1);
 758     bTryAgain = TRUE;
 759     for (t = ten + 2; t > ten - 3 && bTryAgain; t--)
 760     {
 761         for (f = 0; f < NFACT && bTryAgain; f++)
 762         {
 763             space = std::pow(10.0_real, static_cast<real>(t)) * major_fact[f];
 764             /* count how many ticks we would get: */
 765             i = 0;
 766             for (j = 0; j < n; j++)
 767             {
 768                 if (bRmod(axis[j], offset, space))
 769                 {
 770                     i++;
 771                 }
 772             }
 773             /* do we have a reasonable number of ticks ? */
 774             bTryAgain = (i > std::min(10, n - 1)) || (i < 5);
 775         }
 776     }
 777     if (bTryAgain)
 778     {
 779         space = std::max(10 * axis[1] - axis[0], axis[std::min(10, n - 1)] - axis[0]);
 780         fprintf(stderr, "Auto tick spacing failed for %c-axis, guessing %g\n", axisnm, space);
 781     }
 782     *major = space;
 783     *minor = space / minor_fact[(f > 0) ? f - 1 : 0];
 784     fprintf(stderr, "Auto tick spacing for %c-axis: major %g, minor %g\n", axisnm, *major, *minor);
 785 }
 786
 787 static void ps_mat(const char*             outf,
 788                    gmx::ArrayRef<t_matrix> mat,
 789                    gmx::ArrayRef<t_matrix> mat2,
 790                    gmx_bool                bFrame,
 791                    gmx_bool                bDiag,
 792                    gmx_bool                bFirstDiag,
 793                    gmx_bool                bTitle,
 794                    gmx_bool                bTitleOnce,
 795                    gmx_bool                bYonce,
 796                    int                     elegend,
 797                    real                    size,
 798                    real                    boxx,
 799                    real                    boxy,
 800                    const char*             m2p,
 801                    const char*             m2pout,
 802                    int                     mapoffset)
 803 {
 804     t_psrec  psrec, *psr;
 805     int      W, H;
 806     int      x, y, col;
 807     real     x0, y0, xx;
 808     real     w, h, dw, dh;
 809     gmx_bool bMap1, bNextMap1, bDiscrete;
 810
 811     try
 812     {
 813         get_params(m2p, m2pout, &psrec);
 814     }
 815     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 816
 817     psr = &psrec;
 818
 819     if (psr->X.major <= 0)
 820     {
 821         tick_spacing((mat[0].flags & MAT_SPATIAL_X) ? mat[0].nx + 1 : mat[0].nx,
 822                      mat[0].axis_x.data(),
 823                      psr->X.offset,
 824                      'X',
 825                      &(psr->X.major),
 826                      &(psr->X.minor));
 827     }
 828     if (psr->X.minor <= 0)
 829     {
 830         psr->X.minor = psr->X.major / 2;
 831     }
 832     if (psr->Y.major <= 0)
 833     {
 834         tick_spacing((mat[0].flags & MAT_SPATIAL_Y) ? mat[0].ny + 1 : mat[0].ny,
 835                      mat[0].axis_y.data(),
 836                      psr->Y.offset,
 837                      'Y',
 838                      &(psr->Y.major),
 839                      &(psr->Y.minor));
 840     }
 841     if (psr->Y.minor <= 0)
 842     {
 843         psr->Y.minor = psr->Y.major / 2;
 844     }
 845
 846     if (boxx > 0)
 847     {
 848         psr->xboxsize = boxx;
 849         psr->yboxsize = boxx;
 850     }
 851     if (boxy > 0)
 852     {
 853         psr->yboxsize = boxy;
 854     }
 855
 856     if (psr->xboxsize == 0)
 857     {
 858         psr->xboxsize = size / mat[0].nx;
 859         printf("Set the x-size of the box to %.3f\n", psr->xboxsize);
 860     }
 861     if (psr->yboxsize == 0)
 862     {
 863         psr->yboxsize = size / mat[0].nx;
 864         printf("Set the y-size of the box to %.3f\n", psr->yboxsize);
 865     }
 866
 867     gmx::ArrayRef<const t_mapping> map1;
 868     int                            legendIndex = 0;
 869     for (const auto& m : mat)
 870     {
 871         if (m.map.size() > map1.size())
 872         {
 873             if (map1.empty())
 874             {
 875                 printf("Selected legend of matrix # %d for display\n", legendIndex);
 876             }
 877             map1 = m.map;
 878         }
 879         ++legendIndex;
 880     }
 881     gmx::ArrayRef<const t_mapping> map2;
 882     if (!mat2.empty())
 883     {
 884         for (const auto& m : mat2)
 885         {
 886             if (m.map.size() > map2.size())
 887             {
 888                 if (map2.empty())
 889                 {
 890                     printf("Selected legend of matrix # %d for second display\n", legendIndex);
 891                 }
 892                 map2 = m.map;
 893             }
 894             ++legendIndex;
 895         }
 896     }
 897     if (mat[0].legend.empty() && psr->legend)
 898     {
 899         mat[0].legend = psr->leglabel;
 900     }
 901
 902     bTitle          = bTitle && !mat.back().title.empty();
 903     bTitleOnce      = bTitleOnce && !mat.back().title.empty();
 904     psr->bTitle     = bTitle;
 905     psr->bTitleOnce = bTitleOnce;
 906     psr->bYonce     = bYonce;
 907
 908     /* Set up size of box for nice colors */
 909     box_dim(mat, mat2, psr, elegend, bFrame, &w, &h, &dw, &dh);
 910
 911     /* Set up bounding box */
 912     W = static_cast<int>(w + dw);
 913     H = static_cast<int>(h + dh);
 914
 915     /* Start box at */
 916     x0 = dw;
 917     y0 = dh;
 918     x  = static_cast<int>(W + psr->xoffs);
 919     y  = static_cast<int>(H + psr->yoffs);
 920     if (bFrame)
 921     {
 922         x += 5 * DDD;
 923         y += 4 * DDD;
 924     }
 925     t_psdata out = ps_open(outf, 0, 0, x, y);
 926     ps_linewidth(&out, static_cast<int>(psr->linewidth));
 927     ps_init_rgb_box(&out, psr->xboxsize, psr->yboxsize);
 928     ps_init_rgb_nbox(&out, psr->xboxsize, psr->yboxsize);
 929     ps_translate(&out, psr->xoffs, psr->yoffs);
 930
 931     if (bFrame)
 932     {
 933         ps_comment(&out, "Here starts the BOX drawing");
 934         draw_boxes(&out, x0, y0, w, mat, psr);
 935     }
 936
 937     for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
 938     {
 939         if (bTitle || (bTitleOnce && i == gmx::ssize(mat) - 1))
 940         {
 941             /* Print title, if any */
 942             ps_rgb(&out, BLACK);
 943             ps_strfont(&out, psr->titfont, psr->titfontsize);
 944             std::string buf;
 945             if (mat2.empty() || mat[i].title == mat2[i].title)
 946             {
 947                 buf = mat[i].title;
 948             }
 949             else
 950             {
 951                 buf = mat[i].title + " / " + mat2[i].title;
 952             }
 953             ps_ctext(&out, x0 + w / 2, y0 + box_height(mat[i], psr) + psr->titfontsize, buf, eXCenter);
 954         }
 955         ps_comment(&out, gmx::formatString("Here starts the filling of box #%zd", i).c_str());
 956         for (x = 0; (x < mat[i].nx); x++)
 957         {
 958             int nexty;
 959             int nextcol;
 960
 961             xx = x0 + x * psr->xboxsize;
 962             ps_moveto(&out, xx, y0);
 963             y     = 0;
 964             bMap1 = (mat2.empty() || (x < y || (x == y && bFirstDiag)));
 965             if ((bDiag) || (x != y))
 966             {
 967                 col = mat[i].matrix(x, y);
 968             }
 969             else
 970             {
 971                 col = -1;
 972             }
 973             for (nexty = 1; (nexty <= mat[i].ny); nexty++)
 974             {
 975                 bNextMap1 = (mat2.empty() || (x < nexty || (x == nexty && bFirstDiag)));
 976                 /* TRUE:  upper left  -> map1 */
 977                 /* FALSE: lower right -> map2 */
 978                 if ((nexty == mat[i].ny) || (!bDiag && (x == nexty)))
 979                 {
 980                     nextcol = -1;
 981                 }
 982                 else
 983                 {
 984                     nextcol = mat[i].matrix(x, nexty);
 985                 }
 986                 if ((nexty == mat[i].ny) || (col != nextcol) || (bMap1 != bNextMap1))
 987                 {
 988                     if (col >= 0)
 989                     {
 990                         if (bMap1)
 991                         {
 992                             ps_rgb_nbox(&out, &(mat[i].map[col].rgb), nexty - y);
 993                         }
 994                         else
 995                         {
 996                             assert(!mat2.empty());
 997                             ps_rgb_nbox(&out, &(mat2[i].map[col].rgb), nexty - y);
 998                         }
 999                     }
1000                     else
1001                     {
1002                         ps_moverel(&out, 0, psr->yboxsize);
1003                     }
1004                     y     = nexty;
1005                     bMap1 = bNextMap1;
1006                     col   = nextcol;
1007                 }
1008             }
1009         }
1010         y0 += box_height(mat[i], psr) + box_dh(psr) + box_dh_top(i + 1 == gmx::ssize(mat), psr);
1011     }
1012
1013     if (psr->X.lineatzero || psr->Y.lineatzero)
1014     {
1015         /* reset y0 for first box */
1016         y0 = dh;
1017         ps_comment(&out, "Here starts the zero lines drawing");
1018         draw_zerolines(&out, x0, y0, w, mat, psr);
1019     }
1020
1021     if (elegend != elNone)
1022     {
1023         std::string                    legend;
1024         gmx::ArrayRef<const t_mapping> leg_map;
1025         ps_comment(&out, "Now it's legend time!");
1026         ps_linewidth(&out, static_cast<int>(psr->linewidth));
1027         if (mat2.empty() || elegend != elSecond)
1028         {
1029             bDiscrete = mat[0].bDiscrete;
1030             legend    = mat[0].legend;
1031             leg_map   = map1;
1032         }
1033         else
1034         {
1035             bDiscrete = mat2[0].bDiscrete;
1036             legend    = mat2[0].legend;
1037             leg_map   = map2;
1038         }
1039         if (bDiscrete)
1040         {
1041             leg_discrete(&out, psr->legfontsize, DDD, legend, psr->legfontsize, psr->legfont, leg_map);
1042         }
1043         else
1044         {
1045             if (elegend != elBoth)
1046             {
1047                 leg_continuous(
1048                         &out, x0 + w / 2, w / 2, DDD, legend, psr->legfontsize, psr->legfont, leg_map, mapoffset);
1049             }
1050             else
1051             {
1052                 assert(!mat2.empty());
1053                 leg_bicontinuous(
1054                         &out, x0 + w / 2, w, DDD, mat[0].legend, mat2[0].legend, psr->legfontsize, psr->legfont, map1, map2);
1055             }
1056         }
1057         ps_comment(&out, "Done processing");
1058     }
1059
1060     ps_close(&out);
1061 }
1062
1063 static void make_axis_labels(gmx::ArrayRef<t_matrix> mat1)
1064 {
1065     for (auto& m : mat1)
1066     {
1067         /* Make labels for x axis */
1068         if (m.axis_x.empty())
1069         {
1070             m.axis_x.resize(m.nx);
1071             std::iota(m.axis_x.begin(), m.axis_x.end(), 0);
1072         }
1073         /* Make labels for y axis */
1074         if (m.axis_y.empty())
1075         {
1076             m.axis_y.resize(m.ny);
1077             std::iota(m.axis_y.begin(), m.axis_y.end(), 0);
1078         }
1079     }
1080 }
1081
1082 static void prune_mat(gmx::ArrayRef<t_matrix> mat, gmx::ArrayRef<t_matrix> mat2, int skip)
1083 {
1084     GMX_RELEASE_ASSERT(mat.size() == mat2.size() || mat2.empty(),
1085                        "Matrix pruning requires matrices of the same size");
1086     for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
1087     {
1088         fprintf(stderr,
1089                 "converting %dx%d matrix to %dx%d\n",
1090                 mat[i].nx,
1091                 mat[i].ny,
1092                 (mat[i].nx + skip - 1) / skip,
1093                 (mat[i].ny + skip - 1) / skip);
1094         /* walk through matrix */
1095         int xs = 0;
1096         for (int x = 0; (x < mat[i].nx); x++)
1097         {
1098             if (x % skip == 0)
1099             {
1100                 mat[i].axis_x[xs] = mat[i].axis_x[x];
1101                 if (!mat2.empty())
1102                 {
1103                     mat2[i].axis_x[xs] = mat2[i].axis_x[x];
1104                 }
1105                 int ys = 0;
1106                 for (int y = 0; (y < mat[i].ny); y++)
1107                 {
1108                     if (x == 0)
1109                     {
1110                         mat[i].axis_y[ys] = mat[i].axis_y[y];
1111                         if (!mat2.empty())
1112                         {
1113                             mat2[i].axis_y[ys] = mat2[i].axis_y[y];
1114                         }
1115                     }
1116                     if (y % skip == 0)
1117                     {
1118                         mat[i].matrix(xs, ys) = mat[i].matrix(x, y);
1119                         if (!mat2.empty())
1120                         {
1121                             mat2[i].matrix(xs, ys) = mat2[i].matrix(x, y);
1122                         }
1123                         ys++;
1124                     }
1125                 }
1126                 xs++;
1127             }
1128         }
1129         /* adjust parameters */
1130         mat[i].nx = (mat[i].nx + skip - 1) / skip;
1131         mat[i].ny = (mat[i].ny + skip - 1) / skip;
1132         if (!mat2.empty())
1133         {
1134             mat2[i].nx = (mat2[i].nx + skip - 1) / skip;
1135             mat2[i].ny = (mat2[i].ny + skip - 1) / skip;
1136         }
1137     }
1138 }
1139
1140 static void zero_lines(gmx::ArrayRef<t_matrix> mat, gmx::ArrayRef<t_matrix> mat2)
1141 {
1142     GMX_RELEASE_ASSERT(mat.size() == mat2.size(), "zero_lines requires matrices of the same size");
1143     for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
1144     {
1145         for (int m = 0; m < (!mat2.empty() ? 2 : 1); m++)
1146         {
1147             gmx::ArrayRef<t_matrix> mats;
1148             if (m == 0)
1149             {
1150                 mats = mat;
1151             }
1152             else
1153             {
1154                 mats = mat2;
1155             }
1156             for (int x = 0; x < mats[i].nx - 1; x++)
1157             {
1158                 if (std::abs(mats[i].axis_x[x + 1]) < 1e-5)
1159                 {
1160                     for (int y = 0; y < mats[i].ny; y++)
1161                     {
1162                         mats[i].matrix(x, y) = 0;
1163                     }
1164                 }
1165             }
1166             for (int y = 0; y < mats[i].ny - 1; y++)
1167             {
1168                 if (std::abs(mats[i].axis_y[y + 1]) < 1e-5)
1169                 {
1170                     for (int x = 0; x < mats[i].nx; x++)
1171                     {
1172                         mats[i].matrix(x, y) = 0;
1173                     }
1174                 }
1175             }
1176         }
1177     }
1178 }
1179
1180 static void write_combined_matrix(int                     ecombine,
1181                                   const char*             fn,
1182                                   gmx::ArrayRef<t_matrix> mat1,
1183                                   gmx::ArrayRef<t_matrix> mat2,
1184                                   const real*             cmin,
1185                                   const real*             cmax)
1186 {
1187     FILE*  out;
1188     real **rmat1, **rmat2;
1189     real   rhi, rlo;
1190
1191     out = gmx_ffopen(fn, "w");
1192     GMX_RELEASE_ASSERT(mat1.size() == mat2.size(),
1193                        "Combined matrix write requires matrices of the same size");
1194     for (gmx::index k = 0; k != gmx::ssize(mat1); k++)
1195     {
1196         if (mat2[k].nx != mat1[k].nx || mat2[k].ny != mat1[k].ny)
1197         {
1198             gmx_fatal(FARGS,
1199                       "Size of frame %zd in 1st (%dx%d) and 2nd matrix (%dx%d) do"
1200                       " not match.\n",
1201                       k,
1202                       mat1[k].nx,
1203                       mat1[k].ny,
1204                       mat2[k].nx,
1205                       mat2[k].ny);
1206         }
1207         printf("Combining two %dx%d matrices\n", mat1[k].nx, mat1[k].ny);
1208         rmat1 = matrix2real(&mat1[k], nullptr);
1209         rmat2 = matrix2real(&mat2[k], nullptr);
1210         if (nullptr == rmat1 || nullptr == rmat2)
1211         {
1212             gmx_fatal(FARGS,
1213                       "Could not extract real data from %s xpm matrices. Note that, e.g.,\n"
1214                       "g_rms and g_mdmat provide such data, but not do_dssp.\n",
1215                       (nullptr == rmat1 && nullptr == rmat2) ? "both" : "one of the");
1216         }
1217         rlo = 1e38;
1218         rhi = -1e38;
1219         for (int j = 0; j < mat1[k].ny; j++)
1220         {
1221             for (int i = 0; i < mat1[k].nx; i++)
1222             {
1223                 switch (ecombine)
1224                 {
1225                     case ecAdd: rmat1[i][j] += rmat2[i][j]; break;
1226                     case ecSub: rmat1[i][j] -= rmat2[i][j]; break;
1227                     case ecMult: rmat1[i][j] *= rmat2[i][j]; break;
1228                     case ecDiv: rmat1[i][j] /= rmat2[i][j]; break;
1229                     default: gmx_fatal(FARGS, "No such combination rule %d for matrices", ecombine);
1230                 }
1231                 rlo = std::min(rlo, rmat1[i][j]);
1232                 rhi = std::max(rhi, rmat1[i][j]);
1233             }
1234         }
1235         if (cmin)
1236         {
1237             rlo = *cmin;
1238         }
1239         if (cmax)
1240         {
1241             rhi = *cmax;
1242         }
1243         int nlevels = static_cast<int>(std::max(mat1[k].map.size(), mat2[k].map.size()));
1244         if (rhi == rlo)
1245         {
1246             fprintf(stderr, "combination results in uniform matrix (%g), no output\n", rhi);
1247         }
1248         else
1249         {
1250             write_xpm(out,
1251                       mat1[k].flags,
1252                       mat1[k].title,
1253                       mat1[k].legend,
1254                       mat1[k].label_x,
1255                       mat1[k].label_y,
1256                       mat1[k].nx,
1257                       mat1[k].ny,
1258                       mat1[k].axis_x.data(),
1259                       mat1[k].axis_y.data(),
1260                       rmat1,
1261                       rlo,
1262                       rhi,
1263                       white,
1264                       black,
1265                       &nlevels);
1266         }
1267     }
1268     gmx_ffclose(out);
1269 }
1270
1271 static void do_mat(gmx::ArrayRef<t_matrix> mat,
1272                    gmx::ArrayRef<t_matrix> mat2,
1273                    gmx_bool                bFrame,
1274                    gmx_bool                bZeroLine,
1275                    gmx_bool                bDiag,
1276                    gmx_bool                bFirstDiag,
1277                    gmx_bool                bTitle,
1278                    gmx_bool                bTitleOnce,
1279                    gmx_bool                bYonce,
1280                    int                     elegend,
1281                    real                    size,
1282                    real                    boxx,
1283                    real                    boxy,
1284                    const char*             epsfile,
1285                    const char*             xpmfile,
1286                    const char*             m2p,
1287                    const char*             m2pout,
1288                    int                     skip,
1289                    int                     mapoffset)
1290 {
1291     GMX_RELEASE_ASSERT(mat.size() == mat2.size() || mat2.empty(),
1292                        "Combined matrix write requires matrices of the same size");
1293     if (!mat2.empty())
1294     {
1295         for (gmx::index k = 0; k != gmx::ssize(mat); k++)
1296         {
1297             if ((mat2[k].nx != mat[k].nx) || (mat2[k].ny != mat[k].ny))
1298             {
1299                 gmx_fatal(FARGS,
1300                           "WAKE UP!! Size of frame %zd in 2nd matrix file (%dx%d) does not match "
1301                           "size of 1st matrix (%dx%d) or the other way around.\n",
1302                           k,
1303                           mat2[k].nx,
1304                           mat2[k].ny,
1305                           mat[k].nx,
1306                           mat[k].ny);
1307             }
1308             for (int j = 0; (j < mat[k].ny); j++)
1309             {
1310                 for (int i = bFirstDiag ? j + 1 : j; (i < mat[k].nx); i++)
1311                 {
1312                     mat[k].matrix(i, j) = mat2[k].matrix(i, j);
1313                 }
1314             }
1315         }
1316     }
1317     for (gmx::index i = 0; i != gmx::ssize(mat); i++)
1318     {
1319         fprintf(stderr, "Matrix %zd is %d x %d\n", i, mat[i].nx, mat[i].ny);
1320     }
1321
1322     make_axis_labels(mat);
1323
1324     if (skip > 1)
1325     {
1326         prune_mat(mat, mat2, skip);
1327     }
1328
1329     if (bZeroLine)
1330     {
1331         zero_lines(mat, mat);
1332     }
1333
1334     if (epsfile != nullptr)
1335     {
1336         ps_mat(epsfile, mat, mat2, bFrame, bDiag, bFirstDiag, bTitle, bTitleOnce, bYonce, elegend, size, boxx, boxy, m2p, m2pout, mapoffset);
1337     }
1338     if (xpmfile != nullptr)
1339     {
1340         xpm_mat(xpmfile, mat, mat2, bDiag, bFirstDiag);
1341     }
1342 }
1343
1344 static void gradient_map(const rvec grad, gmx::ArrayRef<t_mapping> map)
1345 {
1346     int  i        = 0;
1347     real fraction = 1.0 / (map.size() - 1.0);
1348     for (auto& m : map)
1349     {
1350         real c  = i * fraction;
1351         m.rgb.r = 1 - c * (1 - grad[XX]);
1352         m.rgb.g = 1 - c * (1 - grad[YY]);
1353         m.rgb.b = 1 - c * (1 - grad[ZZ]);
1354         ++i;
1355     }
1356 }
1357
1358 static void gradient_mat(rvec grad, gmx::ArrayRef<t_matrix> mat)
1359 {
1360     for (auto& m : mat)
1361     {
1362         gradient_map(grad, m.map);
1363     }
1364 }
1365
1366 static void rainbow_map(gmx_bool bBlue, gmx::ArrayRef<t_mapping> map)
1367 {
1368     for (auto& m : map)
1369     {
1370         real c = (m.rgb.r + m.rgb.g + m.rgb.b) / 3;
1371         real r, g, b;
1372         if (c > 1)
1373         {
1374             c = 1;
1375         }
1376         if (bBlue)
1377         {
1378             c = 1 - c;
1379         }
1380         if (c <= 0.25) /* 0-0.25 */
1381         {
1382             r = 0;
1383             g = std::pow(4.0 * c, 2.0 / 3.0);
1384             b = 1;
1385         }
1386         else if (c <= 0.5) /* 0.25-0.5 */
1387         {
1388             r = 0;
1389             g = 1;
1390             b = std::pow(2.0 - 4.0 * c, 2.0 / 3.0);
1391         }
1392         else if (c <= 0.75) /* 0.5-0.75 */
1393         {
1394             r = std::pow(4.0 * c - 2.0, 2.0 / 3.0);
1395             g = 1;
1396             b = 0;
1397         }
1398         else /* 0.75-1 */
1399         {
1400             r = 1;
1401             g = std::pow(4.0 - 4.0 * c, 2.0 / 3.0);
1402             b = 0;
1403         }
1404         m.rgb.r = r;
1405         m.rgb.g = g;
1406         m.rgb.b = b;
1407     }
1408 }
1409
1410 static void rainbow_mat(gmx_bool bBlue, gmx::ArrayRef<t_matrix> mat)
1411 {
1412     for (auto& m : mat)
1413     {
1414         rainbow_map(bBlue, m.map);
1415     }
1416 }
1417
1418 int gmx_xpm2ps(int argc, char* argv[])
1419 {
1420     const char* desc[] = {
1421         "[THISMODULE] makes a beautiful color plot of an XPixelMap file.",
1422         "Labels and axis can be displayed, when they are supplied",
1423         "in the correct matrix format.",
1424         "Matrix data may be generated by programs such as [gmx-do_dssp], [gmx-rms] or",
1425         "[gmx-mdmat].[PAR]",
1426         "Parameters are set in the [TT].m2p[tt] file optionally supplied with",
1427         "[TT]-di[tt]. Reasonable defaults are provided. Settings for the [IT]y[it]-axis",
1428         "default to those for the [IT]x[it]-axis. Font names have a defaulting hierarchy:",
1429         "titlefont -> legendfont; titlefont -> (xfont -> yfont -> ytickfont)",
1430         "-> xtickfont, e.g. setting titlefont sets all fonts, setting xfont",
1431         "sets yfont, ytickfont and xtickfont.[PAR]",
1432         "When no [TT].m2p[tt] file is supplied, many settings are taken from",
1433         "command line options. The most important option is [TT]-size[tt],",
1434         "which sets the size of the whole matrix in postscript units.",
1435         "This option can be overridden with the [TT]-bx[tt] and [TT]-by[tt]",
1436         "options (and the corresponding parameters in the [TT].m2p[tt] file),",
1437         "which set the size of a single matrix element.[PAR]",
1438         "With [TT]-f2[tt] a second matrix file can be supplied. Both matrix",
1439         "files will be read simultaneously and the upper left half of the",
1440         "first one ([TT]-f[tt]) is plotted together with the lower right",
1441         "half of the second one ([TT]-f2[tt]). The diagonal will contain",
1442         "values from the matrix file selected with [TT]-diag[tt].",
1443         "Plotting of the diagonal values can be suppressed altogether by",
1444         "setting [TT]-diag[tt] to [TT]none[tt].",
1445         "In this case, a new color map will be generated with",
1446         "a red gradient for negative numbers and a blue for positive.",
1447         "If the color coding and legend labels of both matrices are identical,",
1448         "only one legend will be displayed, else two separate legends are",
1449         "displayed.",
1450         "With [TT]-combine[tt], an alternative operation can be selected",
1451         "to combine the matrices. The output range is automatically set",
1452         "to the actual range of the combined matrix. This can be overridden",
1453         "with [TT]-cmin[tt] and [TT]-cmax[tt].[PAR]",
1454         "[TT]-title[tt] can be set to [TT]none[tt] to suppress the title, or to",
1455         "[TT]ylabel[tt] to show the title in the Y-label position (alongside",
1456         "the [IT]y[it]-axis).[PAR]",
1457         "With the [TT]-rainbow[tt] option, dull grayscale matrices can be turned",
1458         "into attractive color pictures.[PAR]",
1459         "Merged or rainbowed matrices can be written to an XPixelMap file with",
1460         "the [TT]-xpm[tt] option."
1461     };
1462
1463     gmx_output_env_t* oenv;
1464     const char *      fn, *epsfile = nullptr, *xpmfile = nullptr;
1465     int               i, etitle, elegend, ediag, erainbow, ecombine;
1466     gmx_bool          bTitle, bTitleOnce, bDiag, bFirstDiag, bGrad;
1467     static gmx_bool   bFrame = TRUE, bZeroLine = FALSE, bYonce = FALSE;
1468     static real       size = 400, boxx = 0, boxy = 0, cmin = 0, cmax = 0;
1469     static rvec       grad = { 0, 0, 0 };
1470     enum
1471     {
1472         etSel,
1473         etTop,
1474         etOnce,
1475         etYlabel,
1476         etNone,
1477         etNR
1478     };
1479     const char* title[] = { nullptr, "top", "once", "ylabel", "none", nullptr };
1480     /* MUST correspond to enum elXxx as defined at top of file */
1481     const char* legend[] = { nullptr, "both", "first", "second", "none", nullptr };
1482     enum
1483     {
1484         edSel,
1485         edFirst,
1486         edSecond,
1487         edNone,
1488         edNR
1489     };
1490     const char* diag[] = { nullptr, "first", "second", "none", nullptr };
1491     enum
1492     {
1493         erSel,
1494         erNo,
1495         erBlue,
1496         erRed,
1497         erNR
1498     };
1499     const char* rainbow[] = { nullptr, "no", "blue", "red", nullptr };
1500     /* MUST correspond to enum ecXxx as defined at top of file */
1501     const char* combine[] = { nullptr, "halves", "add", "sub", "mult", "div", nullptr };
1502     static int  skip = 1, mapoffset = 0;
1503     t_pargs     pa[] = {
1504         { "-frame", FALSE, etBOOL, { &bFrame }, "Display frame, ticks, labels, title and legend" },
1505         { "-title", FALSE, etENUM, { title }, "Show title at" },
1506         { "-yonce", FALSE, etBOOL, { &bYonce }, "Show y-label only once" },
1507         { "-legend", FALSE, etENUM, { legend }, "Show legend" },
1508         { "-diag", FALSE, etENUM, { diag }, "Diagonal" },
1509         { "-size", FALSE, etREAL, { &size }, "Horizontal size of the matrix in ps units" },
1510         { "-bx",
1511           FALSE,
1512           etREAL,
1513           { &boxx },
1514           "Element x-size, overrides [TT]-size[tt] (also y-size when [TT]-by[tt] is not set)" },
1515         { "-by", FALSE, etREAL, { &boxy }, "Element y-size" },
1516         { "-rainbow", FALSE, etENUM, { rainbow }, "Rainbow colors, convert white to" },
1517         { "-gradient",
1518           FALSE,
1519           etRVEC,
1520           { grad },
1521           "Re-scale colormap to a smooth gradient from white {1,1,1} to {r,g,b}" },
1522         { "-skip", FALSE, etINT, { &skip }, "only write out every nr-th row and column" },
1523         { "-zeroline",
1524           FALSE,
1525           etBOOL,
1526           { &bZeroLine },
1527           "insert line in [REF].xpm[ref] matrix where axis label is zero" },
1528         { "-legoffset",
1529           FALSE,
1530           etINT,
1531           { &mapoffset },
1532           "Skip first N colors from [REF].xpm[ref] file for the legend" },
1533         { "-combine", FALSE, etENUM, { combine }, "Combine two matrices" },
1534         { "-cmin", FALSE, etREAL, { &cmin }, "Minimum for combination output" },
1535         { "-cmax", FALSE, etREAL, { &cmax }, "Maximum for combination output" }
1536     };
1537 #define NPA asize(pa)
1538     t_filenm fnm[] = { { efXPM, "-f", nullptr, ffREAD },      { efXPM, "-f2", "root2", ffOPTRD },
1539                        { efM2P, "-di", nullptr, ffLIBOPTRD }, { efM2P, "-do", "out", ffOPTWR },
1540                        { efEPS, "-o", nullptr, ffOPTWR },     { efXPM, "-xpm", nullptr, ffOPTWR } };
1541 #define NFILE asize(fnm)
1542
1543     if (!parse_common_args(
1544                 &argc, argv, PCA_CAN_VIEW, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
1545     {
1546         return 0;
1547     }
1548
1549     etitle   = nenum(title);
1550     elegend  = nenum(legend);
1551     ediag    = nenum(diag);
1552     erainbow = nenum(rainbow);
1553     ecombine = nenum(combine);
1554     bGrad    = opt2parg_bSet("-gradient", NPA, pa);
1555     for (i = 0; i < DIM; i++)
1556     {
1557         if (grad[i] < 0 || grad[i] > 1)
1558         {
1559             gmx_fatal(FARGS, "RGB value %g out of range (0.0-1.0)", grad[i]);
1560         }
1561     }
1562     if (!bFrame)
1563     {
1564         etitle  = etNone;
1565         elegend = elNone;
1566     }
1567
1568     epsfile = ftp2fn_null(efEPS, NFILE, fnm);
1569     xpmfile = opt2fn_null("-xpm", NFILE, fnm);
1570     if (epsfile == nullptr && xpmfile == nullptr)
1571     {
1572         if (ecombine != ecHalves)
1573         {
1574             xpmfile = opt2fn("-xpm", NFILE, fnm);
1575         }
1576         else
1577         {
1578             epsfile = ftp2fn(efEPS, NFILE, fnm);
1579         }
1580     }
1581     if (ecombine != ecHalves && epsfile)
1582     {
1583         fprintf(stderr,
1584                 "WARNING: can only write result of arithmetic combination "
1585                 "of two matrices to .xpm file\n"
1586                 "         file %s will not be written\n",
1587                 epsfile);
1588         epsfile = nullptr;
1589     }
1590
1591     bDiag      = ediag != edNone;
1592     bFirstDiag = ediag != edSecond;
1593
1594     fn = opt2fn("-f", NFILE, fnm);
1595     std::vector<t_matrix> mat, mat2;
1596     mat = read_xpm_matrix(fn);
1597     fprintf(stderr,
1598             "There %s %zu matri%s in %s\n",
1599             (mat.size() > 1) ? "are" : "is",
1600             mat.size(),
1601             (mat.size() > 1) ? "ces" : "x",
1602             fn);
1603     fn = opt2fn_null("-f2", NFILE, fnm);
1604     if (fn)
1605     {
1606         mat2 = read_xpm_matrix(fn);
1607         fprintf(stderr,
1608                 "There %s %zu matri%s in %s\n",
1609                 (mat2.size() > 1) ? "are" : "is",
1610                 mat2.size(),
1611                 (mat2.size() > 1) ? "ces" : "x",
1612                 fn);
1613         if (mat.size() != mat2.size())
1614         {
1615             fprintf(stderr, "Different number of matrices, using the smallest number.\n");
1616             if (mat.size() > mat2.size())
1617             {
1618                 mat.resize(mat2.size());
1619             }
1620             else
1621             {
1622                 mat2.resize(mat.size());
1623             }
1624         }
1625     }
1626     else
1627     {
1628         if (ecombine != ecHalves)
1629         {
1630             fprintf(stderr,
1631                     "WARNING: arithmetic matrix combination selected (-combine), "
1632                     "but no second matrix (-f2) supplied\n"
1633                     "         no matrix combination will be performed\n");
1634         }
1635         ecombine = 0;
1636     }
1637     bTitle     = etitle == etTop;
1638     bTitleOnce = etitle == etOnce;
1639     if (etitle == etYlabel)
1640     {
1641         for (auto& m : mat)
1642         {
1643             m.label_y = m.title;
1644         }
1645         for (auto& m : mat2)
1646         {
1647             m.label_y = m.title;
1648         }
1649     }
1650     if (bGrad)
1651     {
1652         gradient_mat(grad, mat);
1653         if (!mat2.empty())
1654         {
1655             gradient_mat(grad, mat2);
1656         }
1657     }
1658     if (erainbow != erNo)
1659     {
1660         rainbow_mat(erainbow == erBlue, mat);
1661         if (!mat2.empty())
1662         {
1663             rainbow_mat(erainbow == erBlue, mat2);
1664         }
1665     }
1666
1667     if (mat2.empty() && (elegend != elNone))
1668     {
1669         elegend = elFirst;
1670     }
1671
1672     if (ecombine && ecombine != ecHalves)
1673     {
1674         write_combined_matrix(ecombine,
1675                               xpmfile,
1676                               mat,
1677                               mat2,
1678                               opt2parg_bSet("-cmin", NPA, pa) ? &cmin : nullptr,
1679                               opt2parg_bSet("-cmax", NPA, pa) ? &cmax : nullptr);
1680     }
1681     else
1682     {
1683         do_mat(mat,
1684                mat2,
1685                bFrame,
1686                bZeroLine,
1687                bDiag,
1688                bFirstDiag,
1689                bTitle,
1690                bTitleOnce,
1691                bYonce,
1692                elegend,
1693                size,
1694                boxx,
1695                boxy,
1696                epsfile,
1697                xpmfile,
1698                opt2fn_null("-di", NFILE, fnm),
1699                opt2fn_null("-do", NFILE, fnm),
1700                skip,
1701                mapoffset);
1702     }
1703
1704     view_all(oenv, NFILE, fnm);
1705
1706     return 0;
1707 }