src/gromacs/gmxana/gmx_xpm2ps.cpp

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