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