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