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