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,
 426                              xtick[x], eXCenter);
 427                 }
 428             }
 429             else if (bRmod(m->axis_x[x], psr->X.offset, psr->X.minor)
 430                      && ((m == mat.begin()) || box_do_all_x_min_ticks(psr)))
 431             {
 432                 /* Shorter tick marks */
 433                 ps_line(ps, xx, yy00, xx, yy00 - psr->X.minorticklen);
 434             }
 435             else if (bRmod(m->axis_x[x], psr->X.offset, psr->X.major))
 436             {
 437                 /* Even shorter marks, only each X.major */
 438                 ps_line(ps, xx, yy00, xx, yy00 - (psr->boxspacing / 2));
 439             }
 440         }
 441         ps_strfont(ps, psr->Y.tickfont, psr->Y.tickfontsize);
 442         snew(ytick, nty);
 443         for (int j = 0; (j < nty); j++)
 444         {
 445             sprintf(buf, "%g", m->axis_y[j]);
 446             ytick[j] = gmx_strdup(buf);
 447         }
 448
 449         for (int y = 0; (y < nty); y++)
 450         {
 451             yy = yy00 + (y + offset_y) * psr->yboxsize;
 452             if (bRmod(m->axis_y[y], psr->Y.offset, psr->Y.major) || (psr->Y.first && (y == 0)))
 453             {
 454                 /* Major ticks */
 455                 strlength = std::max(strlength, std::strlen(ytick[y]));
 456                 ps_line(ps, xx00, yy, xx00 - psr->Y.majorticklen, yy);
 457                 ps_ctext(ps, xx00 - psr->Y.majorticklen - DDD, yy - psr->Y.tickfontsize / 3.0,
 458                          ytick[y], eXRight);
 459             }
 460             else if (bRmod(m->axis_y[y], psr->Y.offset, psr->Y.minor))
 461             {
 462                 /* Minor ticks */
 463                 ps_line(ps, xx00, yy, xx00 - psr->Y.minorticklen, yy);
 464             }
 465         }
 466         sfree(xtick);
 467         sfree(ytick);
 468
 469         /* Label on Y-axis */
 470         if (!psr->bYonce || m == halfway)
 471         {
 472             std::string mylab;
 473             if (strlen(psr->Y.label) > 0)
 474             {
 475                 mylab = psr->Y.label;
 476             }
 477             else
 478             {
 479                 mylab = m->label_y;
 480             }
 481             if (!mylab.empty())
 482             {
 483                 ps_strfont(ps, psr->Y.font, psr->Y.fontsize);
 484                 ps_flip(ps, TRUE);
 485                 xxx = x0 - psr->X.majorticklen - psr->X.tickfontsize * strlength - DDD;
 486                 ps_ctext(ps, yy00 + box_height(*m, psr) / 2.0, 612.5 - xxx, mylab, eXCenter);
 487                 ps_flip(ps, FALSE);
 488             }
 489         }
 490
 491         yy00 += box_height(*m, psr) + box_dh(psr) + box_dh_top(m + 1 == mat.end(), psr);
 492     }
 493     /* Label on X-axis */
 494     std::string mylab;
 495     if (strlen(psr->X.label) > 0)
 496     {
 497         mylab = psr->X.label;
 498     }
 499     else
 500     {
 501         mylab = mat[0].label_x;
 502     }
 503     if (!mylab.empty())
 504     {
 505         ps_strfont(ps, psr->X.font, psr->X.fontsize);
 506         ps_ctext(ps, x0 + w / 2,
 507                  y0 - DDD - psr->X.majorticklen - psr->X.tickfontsize * FUDGE - psr->X.fontsize,
 508                  mylab, 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", map1.size(),
 645            map2.size(), map.size());
 646     gmx::index k = 0;
 647     for (gmx::index j = 0; j < gmx::ssize(map1) && k < gmx::ssize(map); ++j, ++k)
 648     {
 649         map[k].code.c1 = mapper[k % nsymbols];
 650         if (map.size() > nsymbols)
 651         {
 652             map[k].code.c2 = mapper[k / nsymbols];
 653         }
 654         map[k].rgb.r = map1[j].rgb.r;
 655         map[k].rgb.g = map1[j].rgb.g;
 656         map[k].rgb.b = map1[j].rgb.b;
 657         map[k].desc  = map1[j].desc;
 658     }
 659     for (gmx::index j = 0; j < gmx::ssize(map2) && k < gmx::ssize(map); ++j, ++k)
 660     {
 661         map[k].code.c1 = mapper[k % nsymbols];
 662         if (map.size() > nsymbols)
 663         {
 664             map[k].code.c2 = mapper[k / nsymbols];
 665         }
 666         map[k].rgb.r = map2[j].rgb.r;
 667         map[k].rgb.g = map2[j].rgb.g;
 668         map[k].rgb.b = map2[j].rgb.b;
 669         map[k].desc  = map2[j].desc;
 670     }
 671
 672     return map;
 673 }
 674
 675 static bool operator==(const t_mapping& lhs, const t_mapping& rhs)
 676 {
 677     return (lhs.rgb.r == rhs.rgb.r && lhs.rgb.g == rhs.rgb.g && lhs.rgb.b == rhs.rgb.b);
 678 }
 679
 680 static void xpm_mat(const char*             outf,
 681                     gmx::ArrayRef<t_matrix> mat,
 682                     gmx::ArrayRef<t_matrix> mat2,
 683                     gmx_bool                bDiag,
 684                     gmx_bool                bFirstDiag)
 685 {
 686     FILE* out;
 687     int   x, y, col;
 688
 689     out = gmx_ffopen(outf, "w");
 690
 691     GMX_RELEASE_ASSERT(mat.size() == mat2.size(),
 692                        "Combined matrix write requires matrices of the same size");
 693     for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
 694     {
 695         // Color maps that differ only in RGB value are considered different
 696         if (mat2.empty() || std::equal(mat[i].map.begin(), mat[i].map.end(), mat2[i].map.begin()))
 697         {
 698             write_xpm_m(out, mat[0]);
 699         }
 700         else
 701         {
 702             auto map = add_maps(mat[i].map, mat2[i].map);
 703             for (x = 0; (x < mat[i].nx); x++)
 704             {
 705                 for (y = 0; (y < mat[i].nx); y++)
 706                 {
 707                     if ((x < y) || ((x == y) && bFirstDiag)) /* upper left  -> map1 */
 708                     {
 709                         col = mat[i].matrix(x, y);
 710                     }
 711                     else /* lower right -> map2 */
 712                     {
 713                         col = mat[i].map.size() + mat[i].matrix(x, y);
 714                     }
 715                     if ((bDiag) || (x != y))
 716                     {
 717                         mat[i].matrix(x, y) = col;
 718                     }
 719                     else
 720                     {
 721                         mat[i].matrix(x, y) = 0;
 722                     }
 723                 }
 724             }
 725             mat[i].map = map;
 726             if (mat[i].title != mat2[i].title)
 727             {
 728                 mat[i].title += " / " + mat2[i].title;
 729             }
 730             if (mat[i].legend != mat2[i].legend)
 731             {
 732                 mat[i].legend += " / " + mat2[i].legend;
 733             }
 734             write_xpm_m(out, mat[i]);
 735         }
 736     }
 737     gmx_ffclose(out);
 738 }
 739
 740 static void tick_spacing(int n, real axis[], real offset, char axisnm, real* major, real* minor)
 741 {
 742     real     space;
 743     gmx_bool bTryAgain;
 744     int      i, j, t, f = 0, ten;
 745 #define NFACT 4
 746     real major_fact[NFACT] = { 5, 4, 2, 1 };
 747     real minor_fact[NFACT] = { 5, 4, 4, 5 };
 748
 749     /* start with interval between 10 matrix points: */
 750     space = std::max(10 * axis[1] - axis[0], axis[std::min(10, n - 1)] - axis[0]);
 751     /* get power of 10 */
 752     ten       = static_cast<int>(std::ceil(std::log(space) / std::log(10.0)) - 1);
 753     bTryAgain = TRUE;
 754     for (t = ten + 2; t > ten - 3 && bTryAgain; t--)
 755     {
 756         for (f = 0; f < NFACT && bTryAgain; f++)
 757         {
 758             space = std::pow(10.0_real, static_cast<real>(t)) * major_fact[f];
 759             /* count how many ticks we would get: */
 760             i = 0;
 761             for (j = 0; j < n; j++)
 762             {
 763                 if (bRmod(axis[j], offset, space))
 764                 {
 765                     i++;
 766                 }
 767             }
 768             /* do we have a reasonable number of ticks ? */
 769             bTryAgain = (i > std::min(10, n - 1)) || (i < 5);
 770         }
 771     }
 772     if (bTryAgain)
 773     {
 774         space = std::max(10 * axis[1] - axis[0], axis[std::min(10, n - 1)] - axis[0]);
 775         fprintf(stderr, "Auto tick spacing failed for %c-axis, guessing %g\n", axisnm, space);
 776     }
 777     *major = space;
 778     *minor = space / minor_fact[(f > 0) ? f - 1 : 0];
 779     fprintf(stderr, "Auto tick spacing for %c-axis: major %g, minor %g\n", axisnm, *major, *minor);
 780 }
 781
 782 static void ps_mat(const char*             outf,
 783                    gmx::ArrayRef<t_matrix> mat,
 784                    gmx::ArrayRef<t_matrix> mat2,
 785                    gmx_bool                bFrame,
 786                    gmx_bool                bDiag,
 787                    gmx_bool                bFirstDiag,
 788                    gmx_bool                bTitle,
 789                    gmx_bool                bTitleOnce,
 790                    gmx_bool                bYonce,
 791                    int                     elegend,
 792                    real                    size,
 793                    real                    boxx,
 794                    real                    boxy,
 795                    const char*             m2p,
 796                    const char*             m2pout,
 797                    int                     mapoffset)
 798 {
 799     t_psrec  psrec, *psr;
 800     int      W, H;
 801     int      x, y, col;
 802     real     x0, y0, xx;
 803     real     w, h, dw, dh;
 804     gmx_bool bMap1, bNextMap1, bDiscrete;
 805
 806     try
 807     {
 808         get_params(m2p, m2pout, &psrec);
 809     }
 810     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 811
 812     psr = &psrec;
 813
 814     if (psr->X.major <= 0)
 815     {
 816         tick_spacing((mat[0].flags & MAT_SPATIAL_X) ? mat[0].nx + 1 : mat[0].nx,
 817                      mat[0].axis_x.data(), psr->X.offset, 'X', &(psr->X.major), &(psr->X.minor));
 818     }
 819     if (psr->X.minor <= 0)
 820     {
 821         psr->X.minor = psr->X.major / 2;
 822     }
 823     if (psr->Y.major <= 0)
 824     {
 825         tick_spacing((mat[0].flags & MAT_SPATIAL_Y) ? mat[0].ny + 1 : mat[0].ny,
 826                      mat[0].axis_y.data(), psr->Y.offset, 'Y', &(psr->Y.major), &(psr->Y.minor));
 827     }
 828     if (psr->Y.minor <= 0)
 829     {
 830         psr->Y.minor = psr->Y.major / 2;
 831     }
 832
 833     if (boxx > 0)
 834     {
 835         psr->xboxsize = boxx;
 836         psr->yboxsize = boxx;
 837     }
 838     if (boxy > 0)
 839     {
 840         psr->yboxsize = boxy;
 841     }
 842
 843     if (psr->xboxsize == 0)
 844     {
 845         psr->xboxsize = size / mat[0].nx;
 846         printf("Set the x-size of the box to %.3f\n", psr->xboxsize);
 847     }
 848     if (psr->yboxsize == 0)
 849     {
 850         psr->yboxsize = size / mat[0].nx;
 851         printf("Set the y-size of the box to %.3f\n", psr->yboxsize);
 852     }
 853
 854     gmx::ArrayRef<const t_mapping> map1;
 855     int                            legendIndex = 0;
 856     for (const auto& m : mat)
 857     {
 858         if (m.map.size() > map1.size())
 859         {
 860             if (map1.empty())
 861             {
 862                 printf("Selected legend of matrix # %d for display\n", legendIndex);
 863             }
 864             map1 = m.map;
 865         }
 866         ++legendIndex;
 867     }
 868     gmx::ArrayRef<const t_mapping> map2;
 869     if (!mat2.empty())
 870     {
 871         for (const auto& m : mat2)
 872         {
 873             if (m.map.size() > map2.size())
 874             {
 875                 if (map2.empty())
 876                 {
 877                     printf("Selected legend of matrix # %d for second display\n", legendIndex);
 878                 }
 879                 map2 = m.map;
 880             }
 881             ++legendIndex;
 882         }
 883     }
 884     if (mat[0].legend.empty() && psr->legend)
 885     {
 886         mat[0].legend = psr->leglabel;
 887     }
 888
 889     bTitle          = bTitle && !mat.back().title.empty();
 890     bTitleOnce      = bTitleOnce && !mat.back().title.empty();
 891     psr->bTitle     = bTitle;
 892     psr->bTitleOnce = bTitleOnce;
 893     psr->bYonce     = bYonce;
 894
 895     /* Set up size of box for nice colors */
 896     box_dim(mat, mat2, psr, elegend, bFrame, &w, &h, &dw, &dh);
 897
 898     /* Set up bounding box */
 899     W = static_cast<int>(w + dw);
 900     H = static_cast<int>(h + dh);
 901
 902     /* Start box at */
 903     x0 = dw;
 904     y0 = dh;
 905     x  = static_cast<int>(W + psr->xoffs);
 906     y  = static_cast<int>(H + psr->yoffs);
 907     if (bFrame)
 908     {
 909         x += 5 * DDD;
 910         y += 4 * DDD;
 911     }
 912     t_psdata out = ps_open(outf, 0, 0, x, y);
 913     ps_linewidth(&out, static_cast<int>(psr->linewidth));
 914     ps_init_rgb_box(&out, psr->xboxsize, psr->yboxsize);
 915     ps_init_rgb_nbox(&out, psr->xboxsize, psr->yboxsize);
 916     ps_translate(&out, psr->xoffs, psr->yoffs);
 917
 918     if (bFrame)
 919     {
 920         ps_comment(&out, "Here starts the BOX drawing");
 921         draw_boxes(&out, x0, y0, w, mat, psr);
 922     }
 923
 924     for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
 925     {
 926         if (bTitle || (bTitleOnce && i == gmx::ssize(mat) - 1))
 927         {
 928             /* Print title, if any */
 929             ps_rgb(&out, BLACK);
 930             ps_strfont(&out, psr->titfont, psr->titfontsize);
 931             std::string buf;
 932             if (!mat2.empty() || mat[i].title == mat2[i].title)
 933             {
 934                 buf = mat[i].title;
 935             }
 936             else
 937             {
 938                 buf = mat[i].title + " / " + mat2[i].title;
 939             }
 940             ps_ctext(&out, x0 + w / 2, y0 + box_height(mat[i], psr) + psr->titfontsize, buf, eXCenter);
 941         }
 942         ps_comment(&out, gmx::formatString("Here starts the filling of box #%zd", i).c_str());
 943         for (x = 0; (x < mat[i].nx); x++)
 944         {
 945             int nexty;
 946             int nextcol;
 947
 948             xx = x0 + x * psr->xboxsize;
 949             ps_moveto(&out, xx, y0);
 950             y     = 0;
 951             bMap1 = (mat2.empty() || (x < y || (x == y && bFirstDiag)));
 952             if ((bDiag) || (x != y))
 953             {
 954                 col = mat[i].matrix(x, y);
 955             }
 956             else
 957             {
 958                 col = -1;
 959             }
 960             for (nexty = 1; (nexty <= mat[i].ny); nexty++)
 961             {
 962                 bNextMap1 = (mat2.empty() || (x < nexty || (x == nexty && bFirstDiag)));
 963                 /* TRUE:  upper left  -> map1 */
 964                 /* FALSE: lower right -> map2 */
 965                 if ((nexty == mat[i].ny) || (!bDiag && (x == nexty)))
 966                 {
 967                     nextcol = -1;
 968                 }
 969                 else
 970                 {
 971                     nextcol = mat[i].matrix(x, nexty);
 972                 }
 973                 if ((nexty == mat[i].ny) || (col != nextcol) || (bMap1 != bNextMap1))
 974                 {
 975                     if (col >= 0)
 976                     {
 977                         if (bMap1)
 978                         {
 979                             ps_rgb_nbox(&out, &(mat[i].map[col].rgb), nexty - y);
 980                         }
 981                         else
 982                         {
 983                             assert(!mat2.empty());
 984                             ps_rgb_nbox(&out, &(mat2[i].map[col].rgb), nexty - y);
 985                         }
 986                     }
 987                     else
 988                     {
 989                         ps_moverel(&out, 0, psr->yboxsize);
 990                     }
 991                     y     = nexty;
 992                     bMap1 = bNextMap1;
 993                     col   = nextcol;
 994                 }
 995             }
 996         }
 997         y0 += box_height(mat[i], psr) + box_dh(psr) + box_dh_top(i + 1 == gmx::ssize(mat), psr);
 998     }
 999
1000     if (psr->X.lineatzero || psr->Y.lineatzero)
1001     {
1002         /* reset y0 for first box */
1003         y0 = dh;
1004         ps_comment(&out, "Here starts the zero lines drawing");
1005         draw_zerolines(&out, x0, y0, w, mat, psr);
1006     }
1007
1008     if (elegend != elNone)
1009     {
1010         std::string                    legend;
1011         gmx::ArrayRef<const t_mapping> leg_map;
1012         ps_comment(&out, "Now it's legend time!");
1013         ps_linewidth(&out, static_cast<int>(psr->linewidth));
1014         if (mat2.empty() || elegend != elSecond)
1015         {
1016             bDiscrete = mat[0].bDiscrete;
1017             legend    = mat[0].legend;
1018             leg_map   = map1;
1019         }
1020         else
1021         {
1022             bDiscrete = mat2[0].bDiscrete;
1023             legend    = mat2[0].legend;
1024             leg_map   = map2;
1025         }
1026         if (bDiscrete)
1027         {
1028             leg_discrete(&out, psr->legfontsize, DDD, legend, psr->legfontsize, psr->legfont, leg_map);
1029         }
1030         else
1031         {
1032             if (elegend != elBoth)
1033             {
1034                 leg_continuous(&out, x0 + w / 2, w / 2, DDD, legend, psr->legfontsize, psr->legfont,
1035                                leg_map, mapoffset);
1036             }
1037             else
1038             {
1039                 assert(!mat2.empty());
1040                 leg_bicontinuous(&out, x0 + w / 2, w, DDD, mat[0].legend, mat2[0].legend,
1041                                  psr->legfontsize, psr->legfont, map1, map2);
1042             }
1043         }
1044         ps_comment(&out, "Done processing");
1045     }
1046
1047     ps_close(&out);
1048 }
1049
1050 static void make_axis_labels(gmx::ArrayRef<t_matrix> mat1)
1051 {
1052     for (auto& m : mat1)
1053     {
1054         /* Make labels for x axis */
1055         if (m.axis_x.empty())
1056         {
1057             m.axis_x.resize(m.nx);
1058             std::iota(m.axis_x.begin(), m.axis_x.end(), 0);
1059         }
1060         /* Make labels for y axis */
1061         if (m.axis_y.empty())
1062         {
1063             m.axis_y.resize(m.ny);
1064             std::iota(m.axis_y.begin(), m.axis_y.end(), 0);
1065         }
1066     }
1067 }
1068
1069 static void prune_mat(gmx::ArrayRef<t_matrix> mat, gmx::ArrayRef<t_matrix> mat2, int skip)
1070 {
1071     GMX_RELEASE_ASSERT(mat.size() == mat2.size() || mat2.empty(),
1072                        "Matrix pruning requires matrices of the same size");
1073     for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
1074     {
1075         fprintf(stderr, "converting %dx%d matrix to %dx%d\n", mat[i].nx, mat[i].ny,
1076                 (mat[i].nx + skip - 1) / skip, (mat[i].ny + skip - 1) / skip);
1077         /* walk through matrix */
1078         int xs = 0;
1079         for (int x = 0; (x < mat[i].nx); x++)
1080         {
1081             if (x % skip == 0)
1082             {
1083                 mat[i].axis_x[xs] = mat[i].axis_x[x];
1084                 if (!mat2.empty())
1085                 {
1086                     mat2[i].axis_x[xs] = mat2[i].axis_x[x];
1087                 }
1088                 int ys = 0;
1089                 for (int y = 0; (y < mat[i].ny); y++)
1090                 {
1091                     if (x == 0)
1092                     {
1093                         mat[i].axis_y[ys] = mat[i].axis_y[y];
1094                         if (!mat2.empty())
1095                         {
1096                             mat2[i].axis_y[ys] = mat2[i].axis_y[y];
1097                         }
1098                     }
1099                     if (y % skip == 0)
1100                     {
1101                         mat[i].matrix(xs, ys) = mat[i].matrix(x, y);
1102                         if (!mat2.empty())
1103                         {
1104                             mat2[i].matrix(xs, ys) = mat2[i].matrix(x, y);
1105                         }
1106                         ys++;
1107                     }
1108                 }
1109                 xs++;
1110             }
1111         }
1112         /* adjust parameters */
1113         mat[i].nx = (mat[i].nx + skip - 1) / skip;
1114         mat[i].ny = (mat[i].ny + skip - 1) / skip;
1115         if (!mat2.empty())
1116         {
1117             mat2[i].nx = (mat2[i].nx + skip - 1) / skip;
1118             mat2[i].ny = (mat2[i].ny + skip - 1) / skip;
1119         }
1120     }
1121 }
1122
1123 static void zero_lines(gmx::ArrayRef<t_matrix> mat, gmx::ArrayRef<t_matrix> mat2)
1124 {
1125     GMX_RELEASE_ASSERT(mat.size() == mat2.size(), "zero_lines requires matrices of the same size");
1126     for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
1127     {
1128         for (int m = 0; m < (!mat2.empty() ? 2 : 1); m++)
1129         {
1130             gmx::ArrayRef<t_matrix> mats;
1131             if (m == 0)
1132             {
1133                 mats = mat;
1134             }
1135             else
1136             {
1137                 mats = mat2;
1138             }
1139             for (int x = 0; x < mats[i].nx - 1; x++)
1140             {
1141                 if (std::abs(mats[i].axis_x[x + 1]) < 1e-5)
1142                 {
1143                     for (int y = 0; y < mats[i].ny; y++)
1144                     {
1145                         mats[i].matrix(x, y) = 0;
1146                     }
1147                 }
1148             }
1149             for (int y = 0; y < mats[i].ny - 1; y++)
1150             {
1151                 if (std::abs(mats[i].axis_y[y + 1]) < 1e-5)
1152                 {
1153                     for (int x = 0; x < mats[i].nx; x++)
1154                     {
1155                         mats[i].matrix(x, y) = 0;
1156                     }
1157                 }
1158             }
1159         }
1160     }
1161 }
1162
1163 static void write_combined_matrix(int                     ecombine,
1164                                   const char*             fn,
1165                                   gmx::ArrayRef<t_matrix> mat1,
1166                                   gmx::ArrayRef<t_matrix> mat2,
1167                                   const real*             cmin,
1168                                   const real*             cmax)
1169 {
1170     FILE*  out;
1171     real **rmat1, **rmat2;
1172     real   rhi, rlo;
1173
1174     out = gmx_ffopen(fn, "w");
1175     GMX_RELEASE_ASSERT(mat1.size() == mat2.size(),
1176                        "Combined matrix write requires matrices of the same size");
1177     for (gmx::index k = 0; k != gmx::ssize(mat1); k++)
1178     {
1179         if (mat2[k].nx != mat1[k].nx || mat2[k].ny != mat1[k].ny)
1180         {
1181             gmx_fatal(FARGS,
1182                       "Size of frame %zd in 1st (%dx%d) and 2nd matrix (%dx%d) do"
1183                       " not match.\n",
1184                       k, mat1[k].nx, mat1[k].ny, mat2[k].nx, mat2[k].ny);
1185         }
1186         printf("Combining two %dx%d matrices\n", mat1[k].nx, mat1[k].ny);
1187         rmat1 = matrix2real(&mat1[k], nullptr);
1188         rmat2 = matrix2real(&mat2[k], nullptr);
1189         if (nullptr == rmat1 || nullptr == rmat2)
1190         {
1191             gmx_fatal(FARGS,
1192                       "Could not extract real data from %s xpm matrices. Note that, e.g.,\n"
1193                       "g_rms and g_mdmat provide such data, but not do_dssp.\n",
1194                       (nullptr == rmat1 && nullptr == rmat2) ? "both" : "one of the");
1195         }
1196         rlo = 1e38;
1197         rhi = -1e38;
1198         for (int j = 0; j < mat1[k].ny; j++)
1199         {
1200             for (int i = 0; i < mat1[k].nx; i++)
1201             {
1202                 switch (ecombine)
1203                 {
1204                     case ecAdd: rmat1[i][j] += rmat2[i][j]; break;
1205                     case ecSub: rmat1[i][j] -= rmat2[i][j]; break;
1206                     case ecMult: rmat1[i][j] *= rmat2[i][j]; break;
1207                     case ecDiv: rmat1[i][j] /= rmat2[i][j]; break;
1208                     default: gmx_fatal(FARGS, "No such combination rule %d for matrices", ecombine);
1209                 }
1210                 rlo = std::min(rlo, rmat1[i][j]);
1211                 rhi = std::max(rhi, rmat1[i][j]);
1212             }
1213         }
1214         if (cmin)
1215         {
1216             rlo = *cmin;
1217         }
1218         if (cmax)
1219         {
1220             rhi = *cmax;
1221         }
1222         int nlevels = static_cast<int>(std::max(mat1[k].map.size(), mat2[k].map.size()));
1223         if (rhi == rlo)
1224         {
1225             fprintf(stderr, "combination results in uniform matrix (%g), no output\n", rhi);
1226         }
1227         else
1228         {
1229             write_xpm(out, mat1[k].flags, mat1[k].title, mat1[k].legend, mat1[k].label_x,
1230                       mat1[k].label_y, mat1[k].nx, mat1[k].ny, mat1[k].axis_x.data(),
1231                       mat1[k].axis_y.data(), rmat1, rlo, rhi, white, black, &nlevels);
1232         }
1233     }
1234     gmx_ffclose(out);
1235 }
1236
1237 static void do_mat(gmx::ArrayRef<t_matrix> mat,
1238                    gmx::ArrayRef<t_matrix> mat2,
1239                    gmx_bool                bFrame,
1240                    gmx_bool                bZeroLine,
1241                    gmx_bool                bDiag,
1242                    gmx_bool                bFirstDiag,
1243                    gmx_bool                bTitle,
1244                    gmx_bool                bTitleOnce,
1245                    gmx_bool                bYonce,
1246                    int                     elegend,
1247                    real                    size,
1248                    real                    boxx,
1249                    real                    boxy,
1250                    const char*             epsfile,
1251                    const char*             xpmfile,
1252                    const char*             m2p,
1253                    const char*             m2pout,
1254                    int                     skip,
1255                    int                     mapoffset)
1256 {
1257     GMX_RELEASE_ASSERT(mat.size() == mat2.size() || mat2.empty(),
1258                        "Combined matrix write requires matrices of the same size");
1259     if (!mat2.empty())
1260     {
1261         for (gmx::index k = 0; k != gmx::ssize(mat); k++)
1262         {
1263             if ((mat2[k].nx != mat[k].nx) || (mat2[k].ny != mat[k].ny))
1264             {
1265                 gmx_fatal(FARGS,
1266                           "WAKE UP!! Size of frame %zd in 2nd matrix file (%dx%d) does not match "
1267                           "size of 1st matrix (%dx%d) or the other way around.\n",
1268                           k, mat2[k].nx, mat2[k].ny, mat[k].nx, mat[k].ny);
1269             }
1270             for (int j = 0; (j < mat[k].ny); j++)
1271             {
1272                 for (int i = bFirstDiag ? j + 1 : j; (i < mat[k].nx); i++)
1273                 {
1274                     mat[k].matrix(i, j) = mat2[k].matrix(i, j);
1275                 }
1276             }
1277         }
1278     }
1279     for (gmx::index i = 0; i != gmx::ssize(mat); i++)
1280     {
1281         fprintf(stderr, "Matrix %zd is %d x %d\n", i, mat[i].nx, mat[i].ny);
1282     }
1283
1284     make_axis_labels(mat);
1285
1286     if (skip > 1)
1287     {
1288         prune_mat(mat, mat2, skip);
1289     }
1290
1291     if (bZeroLine)
1292     {
1293         zero_lines(mat, mat);
1294     }
1295
1296     if (epsfile != nullptr)
1297     {
1298         ps_mat(epsfile, mat, mat2, bFrame, bDiag, bFirstDiag, bTitle, bTitleOnce, bYonce, elegend,
1299                size, boxx, boxy, m2p, m2pout, mapoffset);
1300     }
1301     if (xpmfile != nullptr)
1302     {
1303         xpm_mat(xpmfile, mat, mat2, bDiag, bFirstDiag);
1304     }
1305 }
1306
1307 static void gradient_map(const rvec grad, gmx::ArrayRef<t_mapping> map)
1308 {
1309     int  i        = 0;
1310     real fraction = 1.0 / (map.size() - 1.0);
1311     for (auto& m : map)
1312     {
1313         real c  = i * fraction;
1314         m.rgb.r = 1 - c * (1 - grad[XX]);
1315         m.rgb.g = 1 - c * (1 - grad[YY]);
1316         m.rgb.b = 1 - c * (1 - grad[ZZ]);
1317         ++i;
1318     }
1319 }
1320
1321 static void gradient_mat(rvec grad, gmx::ArrayRef<t_matrix> mat)
1322 {
1323     for (auto& m : mat)
1324     {
1325         gradient_map(grad, m.map);
1326     }
1327 }
1328
1329 static void rainbow_map(gmx_bool bBlue, gmx::ArrayRef<t_mapping> map)
1330 {
1331     for (auto& m : map)
1332     {
1333         real c = (m.rgb.r + m.rgb.g + m.rgb.b) / 3;
1334         real r, g, b;
1335         if (c > 1)
1336         {
1337             c = 1;
1338         }
1339         if (bBlue)
1340         {
1341             c = 1 - c;
1342         }
1343         if (c <= 0.25) /* 0-0.25 */
1344         {
1345             r = 0;
1346             g = std::pow(4.0 * c, 2.0 / 3.0);
1347             b = 1;
1348         }
1349         else if (c <= 0.5) /* 0.25-0.5 */
1350         {
1351             r = 0;
1352             g = 1;
1353             b = std::pow(2.0 - 4.0 * c, 2.0 / 3.0);
1354         }
1355         else if (c <= 0.75) /* 0.5-0.75 */
1356         {
1357             r = std::pow(4.0 * c - 2.0, 2.0 / 3.0);
1358             g = 1;
1359             b = 0;
1360         }
1361         else /* 0.75-1 */
1362         {
1363             r = 1;
1364             g = std::pow(4.0 - 4.0 * c, 2.0 / 3.0);
1365             b = 0;
1366         }
1367         m.rgb.r = r;
1368         m.rgb.g = g;
1369         m.rgb.b = b;
1370     }
1371 }
1372
1373 static void rainbow_mat(gmx_bool bBlue, gmx::ArrayRef<t_matrix> mat)
1374 {
1375     for (auto& m : mat)
1376     {
1377         rainbow_map(bBlue, m.map);
1378     }
1379 }
1380
1381 int gmx_xpm2ps(int argc, char* argv[])
1382 {
1383     const char* desc[] = {
1384         "[THISMODULE] makes a beautiful color plot of an XPixelMap file.",
1385         "Labels and axis can be displayed, when they are supplied",
1386         "in the correct matrix format.",
1387         "Matrix data may be generated by programs such as [gmx-do_dssp], [gmx-rms] or",
1388         "[gmx-mdmat].[PAR]",
1389         "Parameters are set in the [TT].m2p[tt] file optionally supplied with",
1390         "[TT]-di[tt]. Reasonable defaults are provided. Settings for the [IT]y[it]-axis",
1391         "default to those for the [IT]x[it]-axis. Font names have a defaulting hierarchy:",
1392         "titlefont -> legendfont; titlefont -> (xfont -> yfont -> ytickfont)",
1393         "-> xtickfont, e.g. setting titlefont sets all fonts, setting xfont",
1394         "sets yfont, ytickfont and xtickfont.[PAR]",
1395         "When no [TT].m2p[tt] file is supplied, many settings are taken from",
1396         "command line options. The most important option is [TT]-size[tt],",
1397         "which sets the size of the whole matrix in postscript units.",
1398         "This option can be overridden with the [TT]-bx[tt] and [TT]-by[tt]",
1399         "options (and the corresponding parameters in the [TT].m2p[tt] file),",
1400         "which set the size of a single matrix element.[PAR]",
1401         "With [TT]-f2[tt] a second matrix file can be supplied. Both matrix",
1402         "files will be read simultaneously and the upper left half of the",
1403         "first one ([TT]-f[tt]) is plotted together with the lower right",
1404         "half of the second one ([TT]-f2[tt]). The diagonal will contain",
1405         "values from the matrix file selected with [TT]-diag[tt].",
1406         "Plotting of the diagonal values can be suppressed altogether by",
1407         "setting [TT]-diag[tt] to [TT]none[tt].",
1408         "In this case, a new color map will be generated with",
1409         "a red gradient for negative numbers and a blue for positive.",
1410         "If the color coding and legend labels of both matrices are identical,",
1411         "only one legend will be displayed, else two separate legends are",
1412         "displayed.",
1413         "With [TT]-combine[tt], an alternative operation can be selected",
1414         "to combine the matrices. The output range is automatically set",
1415         "to the actual range of the combined matrix. This can be overridden",
1416         "with [TT]-cmin[tt] and [TT]-cmax[tt].[PAR]",
1417         "[TT]-title[tt] can be set to [TT]none[tt] to suppress the title, or to",
1418         "[TT]ylabel[tt] to show the title in the Y-label position (alongside",
1419         "the [IT]y[it]-axis).[PAR]",
1420         "With the [TT]-rainbow[tt] option, dull grayscale matrices can be turned",
1421         "into attractive color pictures.[PAR]",
1422         "Merged or rainbowed matrices can be written to an XPixelMap file with",
1423         "the [TT]-xpm[tt] option."
1424     };
1425
1426     gmx_output_env_t* oenv;
1427     const char *      fn, *epsfile = nullptr, *xpmfile = nullptr;
1428     int               i, etitle, elegend, ediag, erainbow, ecombine;
1429     gmx_bool          bTitle, bTitleOnce, bDiag, bFirstDiag, bGrad;
1430     static gmx_bool   bFrame = TRUE, bZeroLine = FALSE, bYonce = FALSE;
1431     static real       size = 400, boxx = 0, boxy = 0, cmin = 0, cmax = 0;
1432     static rvec       grad = { 0, 0, 0 };
1433     enum
1434     {
1435         etSel,
1436         etTop,
1437         etOnce,
1438         etYlabel,
1439         etNone,
1440         etNR
1441     };
1442     const char* title[] = { nullptr, "top", "once", "ylabel", "none", nullptr };
1443     /* MUST correspond to enum elXxx as defined at top of file */
1444     const char* legend[] = { nullptr, "both", "first", "second", "none", nullptr };
1445     enum
1446     {
1447         edSel,
1448         edFirst,
1449         edSecond,
1450         edNone,
1451         edNR
1452     };
1453     const char* diag[] = { nullptr, "first", "second", "none", nullptr };
1454     enum
1455     {
1456         erSel,
1457         erNo,
1458         erBlue,
1459         erRed,
1460         erNR
1461     };
1462     const char* rainbow[] = { nullptr, "no", "blue", "red", nullptr };
1463     /* MUST correspond to enum ecXxx as defined at top of file */
1464     const char* combine[] = { nullptr, "halves", "add", "sub", "mult", "div", nullptr };
1465     static int  skip = 1, mapoffset = 0;
1466     t_pargs     pa[] = {
1467         { "-frame", FALSE, etBOOL, { &bFrame }, "Display frame, ticks, labels, title and legend" },
1468         { "-title", FALSE, etENUM, { title }, "Show title at" },
1469         { "-yonce", FALSE, etBOOL, { &bYonce }, "Show y-label only once" },
1470         { "-legend", FALSE, etENUM, { legend }, "Show legend" },
1471         { "-diag", FALSE, etENUM, { diag }, "Diagonal" },
1472         { "-size", FALSE, etREAL, { &size }, "Horizontal size of the matrix in ps units" },
1473         { "-bx",
1474           FALSE,
1475           etREAL,
1476           { &boxx },
1477           "Element x-size, overrides [TT]-size[tt] (also y-size when [TT]-by[tt] is not set)" },
1478         { "-by", FALSE, etREAL, { &boxy }, "Element y-size" },
1479         { "-rainbow", FALSE, etENUM, { rainbow }, "Rainbow colors, convert white to" },
1480         { "-gradient",
1481           FALSE,
1482           etRVEC,
1483           { grad },
1484           "Re-scale colormap to a smooth gradient from white {1,1,1} to {r,g,b}" },
1485         { "-skip", FALSE, etINT, { &skip }, "only write out every nr-th row and column" },
1486         { "-zeroline",
1487           FALSE,
1488           etBOOL,
1489           { &bZeroLine },
1490           "insert line in [REF].xpm[ref] matrix where axis label is zero" },
1491         { "-legoffset",
1492           FALSE,
1493           etINT,
1494           { &mapoffset },
1495           "Skip first N colors from [REF].xpm[ref] file for the legend" },
1496         { "-combine", FALSE, etENUM, { combine }, "Combine two matrices" },
1497         { "-cmin", FALSE, etREAL, { &cmin }, "Minimum for combination output" },
1498         { "-cmax", FALSE, etREAL, { &cmax }, "Maximum for combination output" }
1499     };
1500 #define NPA asize(pa)
1501     t_filenm fnm[] = { { efXPM, "-f", nullptr, ffREAD },      { efXPM, "-f2", "root2", ffOPTRD },
1502                        { efM2P, "-di", nullptr, ffLIBOPTRD }, { efM2P, "-do", "out", ffOPTWR },
1503                        { efEPS, "-o", nullptr, ffOPTWR },     { efXPM, "-xpm", nullptr, ffOPTWR } };
1504 #define NFILE asize(fnm)
1505
1506     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW, NFILE, fnm, NPA, pa, asize(desc), desc, 0,
1507                            nullptr, &oenv))
1508     {
1509         return 0;
1510     }
1511
1512     etitle   = nenum(title);
1513     elegend  = nenum(legend);
1514     ediag    = nenum(diag);
1515     erainbow = nenum(rainbow);
1516     ecombine = nenum(combine);
1517     bGrad    = opt2parg_bSet("-gradient", NPA, pa);
1518     for (i = 0; i < DIM; i++)
1519     {
1520         if (grad[i] < 0 || grad[i] > 1)
1521         {
1522             gmx_fatal(FARGS, "RGB value %g out of range (0.0-1.0)", grad[i]);
1523         }
1524     }
1525     if (!bFrame)
1526     {
1527         etitle  = etNone;
1528         elegend = elNone;
1529     }
1530
1531     epsfile = ftp2fn_null(efEPS, NFILE, fnm);
1532     xpmfile = opt2fn_null("-xpm", NFILE, fnm);
1533     if (epsfile == nullptr && xpmfile == nullptr)
1534     {
1535         if (ecombine != ecHalves)
1536         {
1537             xpmfile = opt2fn("-xpm", NFILE, fnm);
1538         }
1539         else
1540         {
1541             epsfile = ftp2fn(efEPS, NFILE, fnm);
1542         }
1543     }
1544     if (ecombine != ecHalves && epsfile)
1545     {
1546         fprintf(stderr,
1547                 "WARNING: can only write result of arithmetic combination "
1548                 "of two matrices to .xpm file\n"
1549                 "         file %s will not be written\n",
1550                 epsfile);
1551         epsfile = nullptr;
1552     }
1553
1554     bDiag      = ediag != edNone;
1555     bFirstDiag = ediag != edSecond;
1556
1557     fn = opt2fn("-f", NFILE, fnm);
1558     std::vector<t_matrix> mat, mat2;
1559     mat = read_xpm_matrix(fn);
1560     fprintf(stderr, "There %s %zu matri%s in %s\n", (mat.size() > 1) ? "are" : "is", mat.size(),
1561             (mat.size() > 1) ? "ces" : "x", fn);
1562     fn = opt2fn_null("-f2", NFILE, fnm);
1563     if (fn)
1564     {
1565         mat2 = read_xpm_matrix(fn);
1566         fprintf(stderr, "There %s %zu matri%s in %s\n", (mat2.size() > 1) ? "are" : "is",
1567                 mat2.size(), (mat2.size() > 1) ? "ces" : "x", fn);
1568         if (mat.size() != mat2.size())
1569         {
1570             fprintf(stderr, "Different number of matrices, using the smallest number.\n");
1571             if (mat.size() > mat2.size())
1572             {
1573                 mat.resize(mat2.size());
1574             }
1575             else
1576             {
1577                 mat2.resize(mat.size());
1578             }
1579         }
1580     }
1581     else
1582     {
1583         if (ecombine != ecHalves)
1584         {
1585             fprintf(stderr,
1586                     "WARNING: arithmetic matrix combination selected (-combine), "
1587                     "but no second matrix (-f2) supplied\n"
1588                     "         no matrix combination will be performed\n");
1589         }
1590         ecombine = 0;
1591     }
1592     bTitle     = etitle == etTop;
1593     bTitleOnce = etitle == etOnce;
1594     if (etitle == etYlabel)
1595     {
1596         for (auto& m : mat)
1597         {
1598             m.label_y = m.title;
1599         }
1600         for (auto& m : mat2)
1601         {
1602             m.label_y = m.title;
1603         }
1604     }
1605     if (bGrad)
1606     {
1607         gradient_mat(grad, mat);
1608         if (!mat2.empty())
1609         {
1610             gradient_mat(grad, mat2);
1611         }
1612     }
1613     if (erainbow != erNo)
1614     {
1615         rainbow_mat(erainbow == erBlue, mat);
1616         if (!mat2.empty())
1617         {
1618             rainbow_mat(erainbow == erBlue, mat2);
1619         }
1620     }
1621
1622     if (mat2.empty() && (elegend != elNone))
1623     {
1624         elegend = elFirst;
1625     }
1626
1627     if (ecombine && ecombine != ecHalves)
1628     {
1629         write_combined_matrix(ecombine, xpmfile, mat, mat2, opt2parg_bSet("-cmin", NPA, pa) ? &cmin : nullptr,
1630                               opt2parg_bSet("-cmax", NPA, pa) ? &cmax : nullptr);
1631     }
1632     else
1633     {
1634         do_mat(mat, mat2, bFrame, bZeroLine, bDiag, bFirstDiag, bTitle, bTitleOnce, bYonce, elegend,
1635                size, boxx, boxy, epsfile, xpmfile, opt2fn_null("-di", NFILE, fnm),
1636                opt2fn_null("-do", NFILE, fnm), skip, mapoffset);
1637     }
1638
1639     view_all(oenv, NFILE, fnm);
1640
1641     return 0;
1642 }