2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
50 #include "gromacs/commandline/pargs.h"
51 #include "gromacs/commandline/viewit.h"
52 #include "gromacs/fileio/matio.h"
53 #include "gromacs/fileio/readinp.h"
54 #include "gromacs/fileio/trxio.h"
55 #include "gromacs/fileio/warninp.h"
56 #include "gromacs/fileio/writeps.h"
57 #include "gromacs/gmxana/gmx_ana.h"
58 #include "gromacs/utility/arrayref.h"
59 #include "gromacs/utility/arraysize.h"
60 #include "gromacs/utility/cstringutil.h"
61 #include "gromacs/utility/exceptions.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/filestream.h"
64 #include "gromacs/utility/futil.h"
65 #include "gromacs/utility/gmxassert.h"
66 #include "gromacs/utility/smalloc.h"
67 #include "gromacs/utility/stringutil.h"
85 char tickfont[STRLEN];
100 char legfont[STRLEN];
101 char leglabel[STRLEN];
102 char leg2label[STRLEN];
112 /* MUST correspond to char *legend[] in main() */
123 /* MUST correspond to char *combine[] in main() */
135 static void get_params(const char* mpin, const char* mpout, t_psrec* psr)
137 static const char* gmx_bools[BOOL_NR + 1] = { "no", "yes", nullptr };
138 /* this must correspond to t_rgb *linecolors[] below */
139 static const char* colors[] = { "none", "black", "white", nullptr };
141 std::vector<t_inpfile> inp;
143 wi = init_warning(FALSE, 0);
147 std::string libmpin = gmx::findLibraryFile(mpin);
148 gmx::TextInputFile stream(libmpin);
149 inp = read_inpfile(&stream, libmpin.c_str(), wi);
156 psr->bw = get_eenum(&inp, "black&white", gmx_bools);
157 psr->linewidth = get_ereal(&inp, "linewidth", 1.0, wi);
158 setStringEntry(&inp, "titlefont", psr->titfont, "Helvetica");
159 psr->titfontsize = get_ereal(&inp, "titlefontsize", 20.0, wi);
160 psr->legend = (get_eenum(&inp, "legend", gmx_bools) != 0);
161 setStringEntry(&inp, "legendfont", psr->legfont, psr->titfont);
162 setStringEntry(&inp, "legendlabel", psr->leglabel, "");
163 setStringEntry(&inp, "legend2label", psr->leg2label, psr->leglabel);
164 psr->legfontsize = get_ereal(&inp, "legendfontsize", 14.0, wi);
165 psr->xboxsize = get_ereal(&inp, "xbox", 0.0, wi);
166 psr->yboxsize = get_ereal(&inp, "ybox", 0.0, wi);
167 psr->boxspacing = get_ereal(&inp, "matrixspacing", 20.0, wi);
168 psr->xoffs = get_ereal(&inp, "xoffset", 0.0, wi);
169 psr->yoffs = get_ereal(&inp, "yoffset", psr->xoffs, wi);
170 psr->boxlinewidth = get_ereal(&inp, "boxlinewidth", psr->linewidth, wi);
171 psr->ticklinewidth = get_ereal(&inp, "ticklinewidth", psr->linewidth, wi);
172 psr->zerolinewidth = get_ereal(&inp, "zerolinewidth", psr->ticklinewidth, wi);
173 psr->X.lineatzero = get_eenum(&inp, "x-lineat0value", colors);
174 psr->X.major = get_ereal(&inp, "x-major", -1, wi);
175 psr->X.minor = get_ereal(&inp, "x-minor", -1, wi);
176 psr->X.offset = get_ereal(&inp, "x-firstmajor", 0.0, wi);
177 psr->X.first = (get_eenum(&inp, "x-majorat0", gmx_bools) != 0);
178 psr->X.majorticklen = get_ereal(&inp, "x-majorticklen", 8.0, wi);
179 psr->X.minorticklen = get_ereal(&inp, "x-minorticklen", 4.0, wi);
180 setStringEntry(&inp, "x-label", psr->X.label, "");
181 psr->X.fontsize = get_ereal(&inp, "x-fontsize", 16.0, wi);
182 setStringEntry(&inp, "x-font", psr->X.font, psr->titfont);
183 psr->X.tickfontsize = get_ereal(&inp, "x-tickfontsize", 10.0, wi);
184 setStringEntry(&inp, "x-tickfont", psr->X.tickfont, psr->X.font);
185 psr->Y.lineatzero = get_eenum(&inp, "y-lineat0value", colors);
186 psr->Y.major = get_ereal(&inp, "y-major", psr->X.major, wi);
187 psr->Y.minor = get_ereal(&inp, "y-minor", psr->X.minor, wi);
188 psr->Y.offset = get_ereal(&inp, "y-firstmajor", psr->X.offset, wi);
189 psr->Y.first = (get_eenum(&inp, "y-majorat0", gmx_bools) != 0);
190 psr->Y.majorticklen = get_ereal(&inp, "y-majorticklen", psr->X.majorticklen, wi);
191 psr->Y.minorticklen = get_ereal(&inp, "y-minorticklen", psr->X.minorticklen, wi);
192 setStringEntry(&inp, "y-label", psr->Y.label, psr->X.label);
193 psr->Y.fontsize = get_ereal(&inp, "y-fontsize", psr->X.fontsize, wi);
194 setStringEntry(&inp, "y-font", psr->Y.font, psr->X.font);
195 psr->Y.tickfontsize = get_ereal(&inp, "y-tickfontsize", psr->X.tickfontsize, wi);
196 setStringEntry(&inp, "y-tickfont", psr->Y.tickfont, psr->Y.font);
198 check_warning_error(wi, FARGS);
200 if (mpout != nullptr)
202 gmx::TextOutputFile stream(mpout);
203 write_inpfile(&stream, mpout, &inp, TRUE, WriteMdpHeader::yes, wi);
207 done_warning(wi, FARGS);
210 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
211 static t_rgb black = { 0, 0, 0 };
212 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
213 static t_rgb white = { 1, 1, 1 };
214 #define BLACK (&black)
215 /* this must correspond to *colors[] in get_params */
216 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
217 static t_rgb* linecolors[] = { nullptr, &black, &white, nullptr };
219 static void leg_discrete(t_psdata* ps,
222 const std::string& label,
225 gmx::ArrayRef<const t_mapping> map)
230 boxhh = fontsize + DDD;
233 ps_strfont(ps, font, fontsize);
234 yhh = y0 + fontsize + 3 * DDD;
237 ps_ctext(ps, x0, yhh, label, eXLeft);
239 ps_moveto(ps, x0, y0);
240 for (const auto& m : map)
244 ps_fillbox(ps, DDD, DDD, DDD + fontsize, boxhh - DDD);
246 ps_box(ps, DDD, DDD, DDD + fontsize, boxhh - DDD);
247 ps_ctext(ps, boxhh + 2 * DDD, fontsize / 3, m.desc, eXLeft);
249 ps_moverel(ps, DDD, -fontsize / 3);
253 static void leg_continuous(t_psdata* ps,
257 const std::string& label,
260 gmx::ArrayRef<const t_mapping> map,
264 real yhh, boxxh, boxyh;
265 gmx::index mapIndex = gmx::ssize(map) - mapoffset;
268 if (x < 8 * fontsize)
272 boxxh = x / mapIndex;
273 if (boxxh > fontsize)
278 GMX_RELEASE_ASSERT(!map.empty(), "NULL map array provided to leg_continuous()");
281 xx0 = x0 - (mapIndex * boxxh) / 2.0;
283 for (gmx::index i = 0; (i < mapIndex); i++)
285 ps_rgb(ps, &(map[i + mapoffset].rgb));
286 ps_fillbox(ps, xx0 + i * boxxh, y0, xx0 + (i + 1) * boxxh, y0 + boxyh);
288 ps_strfont(ps, font, fontsize);
290 ps_box(ps, xx0, y0, xx0 + mapIndex * boxxh, y0 + boxyh);
292 yhh = y0 + boxyh + 3 * DDD;
293 ps_ctext(ps, xx0 + boxxh / 2, yhh, map[0].desc, eXCenter);
296 ps_ctext(ps, x0, yhh, label, eXCenter);
298 ps_ctext(ps, xx0 + (mapIndex * boxxh) - boxxh / 2, yhh, map[map.size() - 1].desc, eXCenter);
301 static void leg_bicontinuous(t_psdata* ps,
305 const std::string& label1,
306 const std::string& label2,
309 gmx::ArrayRef<const t_mapping> map1,
310 gmx::ArrayRef<const t_mapping> map2)
312 real xx1, xx2, x1, x2;
314 x1 = x / (map1.size() + map2.size()) * map1.size(); /* width of legend 1 */
315 x2 = x / (map1.size() + map2.size()) * map2.size(); /* width of legend 2 */
316 xx1 = x0 - (x2 / 2.0) - fontsize; /* center of legend 1 */
317 xx2 = x0 + (x1 / 2.0) + fontsize; /* center of legend 2 */
318 x1 -= fontsize / 2; /* adjust width */
319 x2 -= fontsize / 2; /* adjust width */
320 leg_continuous(ps, xx1, x1, y0, label1, fontsize, font, map1, 0);
321 leg_continuous(ps, xx2, x2, y0, label2, fontsize, font, map2, 0);
324 static real box_height(const t_matrix& mat, t_psrec* psr)
326 return mat.ny * psr->yboxsize;
329 static real box_dh(t_psrec* psr)
331 return psr->boxspacing;
334 static real box_dh_top(gmx_bool bOnce, t_psrec* psr)
338 if (psr->bTitle || (psr->bTitleOnce && bOnce))
340 dh = 2 * psr->titfontsize;
350 static gmx_bool box_do_all_x_maj_ticks(t_psrec* psr)
352 return (psr->boxspacing > (1.5 * psr->X.majorticklen));
355 static gmx_bool box_do_all_x_min_ticks(t_psrec* psr)
357 return (psr->boxspacing > (1.5 * psr->X.minorticklen));
360 static void draw_boxes(t_psdata* ps, real x0, real y0, real w, gmx::ArrayRef<t_matrix> mat, t_psrec* psr)
364 char **xtick, **ytick;
365 real xx, yy, dy, xx00, yy00, offset_x, offset_y;
369 /* Only necessary when there will be no y-labels */
374 ps_linewidth(ps, static_cast<int>(psr->boxlinewidth));
376 for (auto m = mat.begin(); m != mat.end(); ++m)
378 dy = box_height(*m, psr);
379 ps_box(ps, x0 - 1, yy00 - 1, x0 + w + 1, yy00 + dy + 1);
380 yy00 += dy + box_dh(psr) + box_dh_top(m + 1 == mat.end(), psr);
383 /* Draw the ticks on the axes */
384 ps_linewidth(ps, static_cast<int>(psr->ticklinewidth));
387 auto halfway = mat.begin() + (mat.size() / 2);
388 for (auto m = mat.begin(); m != mat.end(); ++m)
390 if (m->flags & MAT_SPATIAL_X)
400 if (m->flags & MAT_SPATIAL_Y)
411 for (int j = 0; (j < ntx); j++)
413 sprintf(buf, "%g", m->axis_x[j]);
414 xtick[j] = gmx_strdup(buf);
416 ps_strfont(ps, psr->X.tickfont, psr->X.tickfontsize);
417 for (x = 0; (x < ntx); x++)
419 xx = xx00 + (x + offset_x) * psr->xboxsize;
420 if ((bRmod(m->axis_x[x], psr->X.offset, psr->X.major) || (psr->X.first && (x == 0)))
421 && (m == mat.begin() || box_do_all_x_maj_ticks(psr)))
423 /* Longer tick marks */
424 ps_line(ps, xx, yy00, xx, yy00 - psr->X.majorticklen);
425 /* Plot label on lowest graph only */
426 if (m == mat.begin())
428 ps_ctext(ps, xx, yy00 - DDD - psr->X.majorticklen - psr->X.tickfontsize * 0.8, xtick[x], eXCenter);
431 else if (bRmod(m->axis_x[x], psr->X.offset, psr->X.minor)
432 && ((m == mat.begin()) || box_do_all_x_min_ticks(psr)))
434 /* Shorter tick marks */
435 ps_line(ps, xx, yy00, xx, yy00 - psr->X.minorticklen);
437 else if (bRmod(m->axis_x[x], psr->X.offset, psr->X.major))
439 /* Even shorter marks, only each X.major */
440 ps_line(ps, xx, yy00, xx, yy00 - (psr->boxspacing / 2));
443 ps_strfont(ps, psr->Y.tickfont, psr->Y.tickfontsize);
445 for (int j = 0; (j < nty); j++)
447 sprintf(buf, "%g", m->axis_y[j]);
448 ytick[j] = gmx_strdup(buf);
451 for (int y = 0; (y < nty); y++)
453 yy = yy00 + (y + offset_y) * psr->yboxsize;
454 if (bRmod(m->axis_y[y], psr->Y.offset, psr->Y.major) || (psr->Y.first && (y == 0)))
457 strlength = std::max(strlength, std::strlen(ytick[y]));
458 ps_line(ps, xx00, yy, xx00 - psr->Y.majorticklen, yy);
459 ps_ctext(ps, xx00 - psr->Y.majorticklen - DDD, yy - psr->Y.tickfontsize / 3.0, ytick[y], eXRight);
461 else if (bRmod(m->axis_y[y], psr->Y.offset, psr->Y.minor))
464 ps_line(ps, xx00, yy, xx00 - psr->Y.minorticklen, yy);
470 /* Label on Y-axis */
471 if (!psr->bYonce || m == halfway)
474 if (strlen(psr->Y.label) > 0)
476 mylab = psr->Y.label;
484 ps_strfont(ps, psr->Y.font, psr->Y.fontsize);
486 xxx = x0 - psr->X.majorticklen - psr->X.tickfontsize * strlength - DDD;
487 ps_ctext(ps, yy00 + box_height(*m, psr) / 2.0, 612.5 - xxx, mylab, eXCenter);
492 yy00 += box_height(*m, psr) + box_dh(psr) + box_dh_top(m + 1 == mat.end(), psr);
494 /* Label on X-axis */
496 if (strlen(psr->X.label) > 0)
498 mylab = psr->X.label;
502 mylab = mat[0].label_x;
506 ps_strfont(ps, psr->X.font, psr->X.fontsize);
509 y0 - DDD - psr->X.majorticklen - psr->X.tickfontsize * FUDGE - psr->X.fontsize,
515 static void draw_zerolines(t_psdata* out, real x0, real y0, real w, gmx::ArrayRef<t_matrix> mat, t_psrec* psr)
517 real xx, yy, dy, xx00, yy00;
522 ps_linewidth(out, static_cast<int>(psr->zerolinewidth));
523 for (auto m = mat.begin(); m != mat.end(); ++m)
525 dy = box_height(*m, psr);
526 /* m->axis_x and _y were already set by draw_boxes */
527 if (psr->X.lineatzero)
529 ps_rgb(out, linecolors[psr->X.lineatzero]);
530 for (x = 0; (x < m->nx); x++)
532 xx = xx00 + (x + 0.7) * psr->xboxsize;
533 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
534 if (x != 0 && x < m->nx - 1
535 && std::abs(m->axis_x[x]) < 0.1 * std::abs(m->axis_x[x + 1] - m->axis_x[x]))
537 ps_line(out, xx, yy00, xx, yy00 + dy + 2);
541 if (psr->Y.lineatzero)
543 ps_rgb(out, linecolors[psr->Y.lineatzero]);
544 for (y = 0; (y < m->ny); y++)
546 yy = yy00 + (y + 0.7) * psr->yboxsize;
547 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
548 if (y != 0 && y < m->ny - 1
549 && std::abs(m->axis_y[y]) < 0.1 * std::abs(m->axis_y[y + 1] - m->axis_y[y]))
551 ps_line(out, xx00, yy, xx00 + w + 2, yy);
555 yy00 += box_height(*m, psr) + box_dh(psr) + box_dh_top(m + 1 == mat.end(), psr);
559 static void box_dim(gmx::ArrayRef<t_matrix> mat,
560 gmx::ArrayRef<t_matrix> mat2,
570 real ww, hh, dww, dhh;
576 for (const auto& m : mat)
578 ww = std::max(ww, m.nx * psr->xboxsize);
579 hh += box_height(m, psr);
580 maxytick = std::max(maxytick, m.nx);
584 if (mat[0].label_y[0])
586 dww += 2.0 * (psr->Y.fontsize + DDD);
588 if (psr->Y.major > 0)
590 dww += psr->Y.majorticklen + DDD
591 + psr->Y.tickfontsize * (std::log(static_cast<real>(maxytick)) / std::log(10.0));
593 else if (psr->Y.minor > 0)
595 dww += psr->Y.minorticklen;
598 if (mat[0].label_x[0])
600 dhh += psr->X.fontsize + 2 * DDD;
602 if (/* fool emacs auto-indent */
603 (elegend == elBoth && (mat[0].legend[0] || (!mat2.empty() && mat2[0].legend[0])))
604 || (elegend == elFirst && mat[0].legend[0])
605 || (elegend == elSecond && (!mat2.empty() && mat2[0].legend[0])))
607 dhh += 2 * (psr->legfontsize * FUDGE + 2 * DDD);
611 dhh += psr->legfontsize * FUDGE + 2 * DDD;
613 if (psr->X.major > 0)
615 dhh += psr->X.tickfontsize * FUDGE + 2 * DDD + psr->X.majorticklen;
617 else if (psr->X.minor > 0)
619 dhh += psr->X.minorticklen;
622 hh += (mat.size() - 1) * box_dh(psr);
623 hh += box_dh_top(TRUE, psr);
626 hh += (mat.size() - 1) * box_dh_top(FALSE, psr);
635 static std::vector<t_mapping> add_maps(gmx::ArrayRef<t_mapping> map1, gmx::ArrayRef<t_mapping> map2)
637 static char mapper[] =
638 "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<"
640 std::vector<t_mapping> map(map1.size() + map2.size());
642 size_t nsymbols = std::strlen(mapper);
643 if (map.size() > nsymbols * nsymbols)
645 gmx_fatal(FARGS, "Not enough symbols to merge the two colormaps\n");
647 printf("Combining colormaps of %zu and %zu elements into one of %zu elements\n",
652 for (gmx::index j = 0; j < gmx::ssize(map1) && k < gmx::ssize(map); ++j, ++k)
654 map[k].code.c1 = mapper[k % nsymbols];
655 if (map.size() > nsymbols)
657 map[k].code.c2 = mapper[k / nsymbols];
659 map[k].rgb.r = map1[j].rgb.r;
660 map[k].rgb.g = map1[j].rgb.g;
661 map[k].rgb.b = map1[j].rgb.b;
662 map[k].desc = map1[j].desc;
664 for (gmx::index j = 0; j < gmx::ssize(map2) && k < gmx::ssize(map); ++j, ++k)
666 map[k].code.c1 = mapper[k % nsymbols];
667 if (map.size() > nsymbols)
669 map[k].code.c2 = mapper[k / nsymbols];
671 map[k].rgb.r = map2[j].rgb.r;
672 map[k].rgb.g = map2[j].rgb.g;
673 map[k].rgb.b = map2[j].rgb.b;
674 map[k].desc = map2[j].desc;
680 static bool operator==(const t_mapping& lhs, const t_mapping& rhs)
682 return (lhs.rgb.r == rhs.rgb.r && lhs.rgb.g == rhs.rgb.g && lhs.rgb.b == rhs.rgb.b);
685 static void xpm_mat(const char* outf,
686 gmx::ArrayRef<t_matrix> mat,
687 gmx::ArrayRef<t_matrix> mat2,
694 out = gmx_ffopen(outf, "w");
696 GMX_RELEASE_ASSERT(mat.size() == mat2.size(),
697 "Combined matrix write requires matrices of the same size");
698 for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
700 // Color maps that differ only in RGB value are considered different
701 if (mat2.empty() || std::equal(mat[i].map.begin(), mat[i].map.end(), mat2[i].map.begin()))
703 write_xpm_m(out, mat[0]);
707 auto map = add_maps(mat[i].map, mat2[i].map);
708 for (x = 0; (x < mat[i].nx); x++)
710 for (y = 0; (y < mat[i].nx); y++)
712 if ((x < y) || ((x == y) && bFirstDiag)) /* upper left -> map1 */
714 col = mat[i].matrix(x, y);
716 else /* lower right -> map2 */
718 col = mat[i].map.size() + mat[i].matrix(x, y);
720 if ((bDiag) || (x != y))
722 mat[i].matrix(x, y) = col;
726 mat[i].matrix(x, y) = 0;
731 if (mat[i].title != mat2[i].title)
733 mat[i].title += " / " + mat2[i].title;
735 if (mat[i].legend != mat2[i].legend)
737 mat[i].legend += " / " + mat2[i].legend;
739 write_xpm_m(out, mat[i]);
745 static void tick_spacing(int n, real axis[], real offset, char axisnm, real* major, real* minor)
749 int i, j, t, f = 0, ten;
751 real major_fact[NFACT] = { 5, 4, 2, 1 };
752 real minor_fact[NFACT] = { 5, 4, 4, 5 };
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);
759 for (t = ten + 2; t > ten - 3 && bTryAgain; t--)
761 for (f = 0; f < NFACT && bTryAgain; f++)
763 space = std::pow(10.0_real, static_cast<real>(t)) * major_fact[f];
764 /* count how many ticks we would get: */
766 for (j = 0; j < n; j++)
768 if (bRmod(axis[j], offset, space))
773 /* do we have a reasonable number of ticks ? */
774 bTryAgain = (i > std::min(10, n - 1)) || (i < 5);
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", axisnm, space);
783 *minor = space / minor_fact[(f > 0) ? f - 1 : 0];
784 fprintf(stderr, "Auto tick spacing for %c-axis: major %g, minor %g\n", axisnm, *major, *minor);
787 static void ps_mat(const char* outf,
788 gmx::ArrayRef<t_matrix> mat,
789 gmx::ArrayRef<t_matrix> mat2,
809 gmx_bool bMap1, bNextMap1, bDiscrete;
813 get_params(m2p, m2pout, &psrec);
815 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
819 if (psr->X.major <= 0)
821 tick_spacing((mat[0].flags & MAT_SPATIAL_X) ? mat[0].nx + 1 : mat[0].nx,
822 mat[0].axis_x.data(),
828 if (psr->X.minor <= 0)
830 psr->X.minor = psr->X.major / 2;
832 if (psr->Y.major <= 0)
834 tick_spacing((mat[0].flags & MAT_SPATIAL_Y) ? mat[0].ny + 1 : mat[0].ny,
835 mat[0].axis_y.data(),
841 if (psr->Y.minor <= 0)
843 psr->Y.minor = psr->Y.major / 2;
848 psr->xboxsize = boxx;
849 psr->yboxsize = boxx;
853 psr->yboxsize = boxy;
856 if (psr->xboxsize == 0)
858 psr->xboxsize = size / mat[0].nx;
859 printf("Set the x-size of the box to %.3f\n", psr->xboxsize);
861 if (psr->yboxsize == 0)
863 psr->yboxsize = size / mat[0].nx;
864 printf("Set the y-size of the box to %.3f\n", psr->yboxsize);
867 gmx::ArrayRef<const t_mapping> map1;
869 for (const auto& m : mat)
871 if (m.map.size() > map1.size())
875 printf("Selected legend of matrix # %d for display\n", legendIndex);
881 gmx::ArrayRef<const t_mapping> map2;
884 for (const auto& m : mat2)
886 if (m.map.size() > map2.size())
890 printf("Selected legend of matrix # %d for second display\n", legendIndex);
897 if (mat[0].legend.empty() && psr->legend)
899 mat[0].legend = psr->leglabel;
902 bTitle = bTitle && !mat.back().title.empty();
903 bTitleOnce = bTitleOnce && !mat.back().title.empty();
904 psr->bTitle = bTitle;
905 psr->bTitleOnce = bTitleOnce;
906 psr->bYonce = bYonce;
908 /* Set up size of box for nice colors */
909 box_dim(mat, mat2, psr, elegend, bFrame, &w, &h, &dw, &dh);
911 /* Set up bounding box */
912 W = static_cast<int>(w + dw);
913 H = static_cast<int>(h + dh);
918 x = static_cast<int>(W + psr->xoffs);
919 y = static_cast<int>(H + psr->yoffs);
925 t_psdata out = ps_open(outf, 0, 0, x, y);
926 ps_linewidth(&out, static_cast<int>(psr->linewidth));
927 ps_init_rgb_box(&out, psr->xboxsize, psr->yboxsize);
928 ps_init_rgb_nbox(&out, psr->xboxsize, psr->yboxsize);
929 ps_translate(&out, psr->xoffs, psr->yoffs);
933 ps_comment(&out, "Here starts the BOX drawing");
934 draw_boxes(&out, x0, y0, w, mat, psr);
937 for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
939 if (bTitle || (bTitleOnce && i == gmx::ssize(mat) - 1))
941 /* Print title, if any */
943 ps_strfont(&out, psr->titfont, psr->titfontsize);
945 if (mat2.empty() || mat[i].title == mat2[i].title)
951 buf = mat[i].title + " / " + mat2[i].title;
953 ps_ctext(&out, x0 + w / 2, y0 + box_height(mat[i], psr) + psr->titfontsize, buf, eXCenter);
955 ps_comment(&out, gmx::formatString("Here starts the filling of box #%zd", i).c_str());
956 for (x = 0; (x < mat[i].nx); x++)
961 xx = x0 + x * psr->xboxsize;
962 ps_moveto(&out, xx, y0);
964 bMap1 = (mat2.empty() || (x < y || (x == y && bFirstDiag)));
965 if ((bDiag) || (x != y))
967 col = mat[i].matrix(x, y);
973 for (nexty = 1; (nexty <= mat[i].ny); nexty++)
975 bNextMap1 = (mat2.empty() || (x < nexty || (x == nexty && bFirstDiag)));
976 /* TRUE: upper left -> map1 */
977 /* FALSE: lower right -> map2 */
978 if ((nexty == mat[i].ny) || (!bDiag && (x == nexty)))
984 nextcol = mat[i].matrix(x, nexty);
986 if ((nexty == mat[i].ny) || (col != nextcol) || (bMap1 != bNextMap1))
992 ps_rgb_nbox(&out, &(mat[i].map[col].rgb), nexty - y);
996 assert(!mat2.empty());
997 ps_rgb_nbox(&out, &(mat2[i].map[col].rgb), nexty - y);
1002 ps_moverel(&out, 0, psr->yboxsize);
1010 y0 += box_height(mat[i], psr) + box_dh(psr) + box_dh_top(i + 1 == gmx::ssize(mat), psr);
1013 if (psr->X.lineatzero || psr->Y.lineatzero)
1015 /* reset y0 for first box */
1017 ps_comment(&out, "Here starts the zero lines drawing");
1018 draw_zerolines(&out, x0, y0, w, mat, psr);
1021 if (elegend != elNone)
1024 gmx::ArrayRef<const t_mapping> leg_map;
1025 ps_comment(&out, "Now it's legend time!");
1026 ps_linewidth(&out, static_cast<int>(psr->linewidth));
1027 if (mat2.empty() || elegend != elSecond)
1029 bDiscrete = mat[0].bDiscrete;
1030 legend = mat[0].legend;
1035 bDiscrete = mat2[0].bDiscrete;
1036 legend = mat2[0].legend;
1041 leg_discrete(&out, psr->legfontsize, DDD, legend, psr->legfontsize, psr->legfont, leg_map);
1045 if (elegend != elBoth)
1048 &out, x0 + w / 2, w / 2, DDD, legend, psr->legfontsize, psr->legfont, leg_map, mapoffset);
1052 assert(!mat2.empty());
1054 &out, x0 + w / 2, w, DDD, mat[0].legend, mat2[0].legend, psr->legfontsize, psr->legfont, map1, map2);
1057 ps_comment(&out, "Done processing");
1063 static void make_axis_labels(gmx::ArrayRef<t_matrix> mat1)
1065 for (auto& m : mat1)
1067 /* Make labels for x axis */
1068 if (m.axis_x.empty())
1070 m.axis_x.resize(m.nx);
1071 std::iota(m.axis_x.begin(), m.axis_x.end(), 0);
1073 /* Make labels for y axis */
1074 if (m.axis_y.empty())
1076 m.axis_y.resize(m.ny);
1077 std::iota(m.axis_y.begin(), m.axis_y.end(), 0);
1082 static void prune_mat(gmx::ArrayRef<t_matrix> mat, gmx::ArrayRef<t_matrix> mat2, int skip)
1084 GMX_RELEASE_ASSERT(mat.size() == mat2.size() || mat2.empty(),
1085 "Matrix pruning requires matrices of the same size");
1086 for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
1089 "converting %dx%d matrix to %dx%d\n",
1092 (mat[i].nx + skip - 1) / skip,
1093 (mat[i].ny + skip - 1) / skip);
1094 /* walk through matrix */
1096 for (int x = 0; (x < mat[i].nx); x++)
1100 mat[i].axis_x[xs] = mat[i].axis_x[x];
1103 mat2[i].axis_x[xs] = mat2[i].axis_x[x];
1106 for (int y = 0; (y < mat[i].ny); y++)
1110 mat[i].axis_y[ys] = mat[i].axis_y[y];
1113 mat2[i].axis_y[ys] = mat2[i].axis_y[y];
1118 mat[i].matrix(xs, ys) = mat[i].matrix(x, y);
1121 mat2[i].matrix(xs, ys) = mat2[i].matrix(x, y);
1129 /* adjust parameters */
1130 mat[i].nx = (mat[i].nx + skip - 1) / skip;
1131 mat[i].ny = (mat[i].ny + skip - 1) / skip;
1134 mat2[i].nx = (mat2[i].nx + skip - 1) / skip;
1135 mat2[i].ny = (mat2[i].ny + skip - 1) / skip;
1140 static void zero_lines(gmx::ArrayRef<t_matrix> mat, gmx::ArrayRef<t_matrix> mat2)
1142 GMX_RELEASE_ASSERT(mat.size() == mat2.size(), "zero_lines requires matrices of the same size");
1143 for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
1145 for (int m = 0; m < (!mat2.empty() ? 2 : 1); m++)
1147 gmx::ArrayRef<t_matrix> mats;
1156 for (int x = 0; x < mats[i].nx - 1; x++)
1158 if (std::abs(mats[i].axis_x[x + 1]) < 1e-5)
1160 for (int y = 0; y < mats[i].ny; y++)
1162 mats[i].matrix(x, y) = 0;
1166 for (int y = 0; y < mats[i].ny - 1; y++)
1168 if (std::abs(mats[i].axis_y[y + 1]) < 1e-5)
1170 for (int x = 0; x < mats[i].nx; x++)
1172 mats[i].matrix(x, y) = 0;
1180 static void write_combined_matrix(int ecombine,
1182 gmx::ArrayRef<t_matrix> mat1,
1183 gmx::ArrayRef<t_matrix> mat2,
1188 real **rmat1, **rmat2;
1191 out = gmx_ffopen(fn, "w");
1192 GMX_RELEASE_ASSERT(mat1.size() == mat2.size(),
1193 "Combined matrix write requires matrices of the same size");
1194 for (gmx::index k = 0; k != gmx::ssize(mat1); k++)
1196 if (mat2[k].nx != mat1[k].nx || mat2[k].ny != mat1[k].ny)
1199 "Size of frame %zd in 1st (%dx%d) and 2nd matrix (%dx%d) do"
1207 printf("Combining two %dx%d matrices\n", mat1[k].nx, mat1[k].ny);
1208 rmat1 = matrix2real(&mat1[k], nullptr);
1209 rmat2 = matrix2real(&mat2[k], nullptr);
1210 if (nullptr == rmat1 || nullptr == rmat2)
1213 "Could not extract real data from %s xpm matrices. Note that, e.g.,\n"
1214 "g_rms and g_mdmat provide such data, but not do_dssp.\n",
1215 (nullptr == rmat1 && nullptr == rmat2) ? "both" : "one of the");
1219 for (int j = 0; j < mat1[k].ny; j++)
1221 for (int i = 0; i < mat1[k].nx; i++)
1225 case ecAdd: rmat1[i][j] += rmat2[i][j]; break;
1226 case ecSub: rmat1[i][j] -= rmat2[i][j]; break;
1227 case ecMult: rmat1[i][j] *= rmat2[i][j]; break;
1228 case ecDiv: rmat1[i][j] /= rmat2[i][j]; break;
1229 default: gmx_fatal(FARGS, "No such combination rule %d for matrices", ecombine);
1231 rlo = std::min(rlo, rmat1[i][j]);
1232 rhi = std::max(rhi, rmat1[i][j]);
1243 int nlevels = static_cast<int>(std::max(mat1[k].map.size(), mat2[k].map.size()));
1246 fprintf(stderr, "combination results in uniform matrix (%g), no output\n", rhi);
1258 mat1[k].axis_x.data(),
1259 mat1[k].axis_y.data(),
1271 static void do_mat(gmx::ArrayRef<t_matrix> mat,
1272 gmx::ArrayRef<t_matrix> mat2,
1276 gmx_bool bFirstDiag,
1278 gmx_bool bTitleOnce,
1284 const char* epsfile,
1285 const char* xpmfile,
1291 GMX_RELEASE_ASSERT(mat.size() == mat2.size() || mat2.empty(),
1292 "Combined matrix write requires matrices of the same size");
1295 for (gmx::index k = 0; k != gmx::ssize(mat); k++)
1297 if ((mat2[k].nx != mat[k].nx) || (mat2[k].ny != mat[k].ny))
1300 "WAKE UP!! Size of frame %zd in 2nd matrix file (%dx%d) does not match "
1301 "size of 1st matrix (%dx%d) or the other way around.\n",
1308 for (int j = 0; (j < mat[k].ny); j++)
1310 for (int i = bFirstDiag ? j + 1 : j; (i < mat[k].nx); i++)
1312 mat[k].matrix(i, j) = mat2[k].matrix(i, j);
1317 for (gmx::index i = 0; i != gmx::ssize(mat); i++)
1319 fprintf(stderr, "Matrix %zd is %d x %d\n", i, mat[i].nx, mat[i].ny);
1322 make_axis_labels(mat);
1326 prune_mat(mat, mat2, skip);
1331 zero_lines(mat, mat);
1334 if (epsfile != nullptr)
1336 ps_mat(epsfile, mat, mat2, bFrame, bDiag, bFirstDiag, bTitle, bTitleOnce, bYonce, elegend, size, boxx, boxy, m2p, m2pout, mapoffset);
1338 if (xpmfile != nullptr)
1340 xpm_mat(xpmfile, mat, mat2, bDiag, bFirstDiag);
1344 static void gradient_map(const rvec grad, gmx::ArrayRef<t_mapping> map)
1347 real fraction = 1.0 / (map.size() - 1.0);
1350 real c = i * fraction;
1351 m.rgb.r = 1 - c * (1 - grad[XX]);
1352 m.rgb.g = 1 - c * (1 - grad[YY]);
1353 m.rgb.b = 1 - c * (1 - grad[ZZ]);
1358 static void gradient_mat(rvec grad, gmx::ArrayRef<t_matrix> mat)
1362 gradient_map(grad, m.map);
1366 static void rainbow_map(gmx_bool bBlue, gmx::ArrayRef<t_mapping> map)
1370 real c = (m.rgb.r + m.rgb.g + m.rgb.b) / 3;
1380 if (c <= 0.25) /* 0-0.25 */
1383 g = std::pow(4.0 * c, 2.0 / 3.0);
1386 else if (c <= 0.5) /* 0.25-0.5 */
1390 b = std::pow(2.0 - 4.0 * c, 2.0 / 3.0);
1392 else if (c <= 0.75) /* 0.5-0.75 */
1394 r = std::pow(4.0 * c - 2.0, 2.0 / 3.0);
1401 g = std::pow(4.0 - 4.0 * c, 2.0 / 3.0);
1410 static void rainbow_mat(gmx_bool bBlue, gmx::ArrayRef<t_matrix> mat)
1414 rainbow_map(bBlue, m.map);
1418 int gmx_xpm2ps(int argc, char* argv[])
1420 const char* desc[] = {
1421 "[THISMODULE] makes a beautiful color plot of an XPixelMap file.",
1422 "Labels and axis can be displayed, when they are supplied",
1423 "in the correct matrix format.",
1424 "Matrix data may be generated by programs such as [gmx-do_dssp], [gmx-rms] or",
1425 "[gmx-mdmat].[PAR]",
1426 "Parameters are set in the [TT].m2p[tt] file optionally supplied with",
1427 "[TT]-di[tt]. Reasonable defaults are provided. Settings for the [IT]y[it]-axis",
1428 "default to those for the [IT]x[it]-axis. Font names have a defaulting hierarchy:",
1429 "titlefont -> legendfont; titlefont -> (xfont -> yfont -> ytickfont)",
1430 "-> xtickfont, e.g. setting titlefont sets all fonts, setting xfont",
1431 "sets yfont, ytickfont and xtickfont.[PAR]",
1432 "When no [TT].m2p[tt] file is supplied, many settings are taken from",
1433 "command line options. The most important option is [TT]-size[tt],",
1434 "which sets the size of the whole matrix in postscript units.",
1435 "This option can be overridden with the [TT]-bx[tt] and [TT]-by[tt]",
1436 "options (and the corresponding parameters in the [TT].m2p[tt] file),",
1437 "which set the size of a single matrix element.[PAR]",
1438 "With [TT]-f2[tt] a second matrix file can be supplied. Both matrix",
1439 "files will be read simultaneously and the upper left half of the",
1440 "first one ([TT]-f[tt]) is plotted together with the lower right",
1441 "half of the second one ([TT]-f2[tt]). The diagonal will contain",
1442 "values from the matrix file selected with [TT]-diag[tt].",
1443 "Plotting of the diagonal values can be suppressed altogether by",
1444 "setting [TT]-diag[tt] to [TT]none[tt].",
1445 "In this case, a new color map will be generated with",
1446 "a red gradient for negative numbers and a blue for positive.",
1447 "If the color coding and legend labels of both matrices are identical,",
1448 "only one legend will be displayed, else two separate legends are",
1450 "With [TT]-combine[tt], an alternative operation can be selected",
1451 "to combine the matrices. The output range is automatically set",
1452 "to the actual range of the combined matrix. This can be overridden",
1453 "with [TT]-cmin[tt] and [TT]-cmax[tt].[PAR]",
1454 "[TT]-title[tt] can be set to [TT]none[tt] to suppress the title, or to",
1455 "[TT]ylabel[tt] to show the title in the Y-label position (alongside",
1456 "the [IT]y[it]-axis).[PAR]",
1457 "With the [TT]-rainbow[tt] option, dull grayscale matrices can be turned",
1458 "into attractive color pictures.[PAR]",
1459 "Merged or rainbowed matrices can be written to an XPixelMap file with",
1460 "the [TT]-xpm[tt] option."
1463 gmx_output_env_t* oenv;
1464 const char * fn, *epsfile = nullptr, *xpmfile = nullptr;
1465 int i, etitle, elegend, ediag, erainbow, ecombine;
1466 gmx_bool bTitle, bTitleOnce, bDiag, bFirstDiag, bGrad;
1467 static gmx_bool bFrame = TRUE, bZeroLine = FALSE, bYonce = FALSE;
1468 static real size = 400, boxx = 0, boxy = 0, cmin = 0, cmax = 0;
1469 static rvec grad = { 0, 0, 0 };
1479 const char* title[] = { nullptr, "top", "once", "ylabel", "none", nullptr };
1480 /* MUST correspond to enum elXxx as defined at top of file */
1481 const char* legend[] = { nullptr, "both", "first", "second", "none", nullptr };
1490 const char* diag[] = { nullptr, "first", "second", "none", nullptr };
1499 const char* rainbow[] = { nullptr, "no", "blue", "red", nullptr };
1500 /* MUST correspond to enum ecXxx as defined at top of file */
1501 const char* combine[] = { nullptr, "halves", "add", "sub", "mult", "div", nullptr };
1502 static int skip = 1, mapoffset = 0;
1504 { "-frame", FALSE, etBOOL, { &bFrame }, "Display frame, ticks, labels, title and legend" },
1505 { "-title", FALSE, etENUM, { title }, "Show title at" },
1506 { "-yonce", FALSE, etBOOL, { &bYonce }, "Show y-label only once" },
1507 { "-legend", FALSE, etENUM, { legend }, "Show legend" },
1508 { "-diag", FALSE, etENUM, { diag }, "Diagonal" },
1509 { "-size", FALSE, etREAL, { &size }, "Horizontal size of the matrix in ps units" },
1514 "Element x-size, overrides [TT]-size[tt] (also y-size when [TT]-by[tt] is not set)" },
1515 { "-by", FALSE, etREAL, { &boxy }, "Element y-size" },
1516 { "-rainbow", FALSE, etENUM, { rainbow }, "Rainbow colors, convert white to" },
1521 "Re-scale colormap to a smooth gradient from white {1,1,1} to {r,g,b}" },
1522 { "-skip", FALSE, etINT, { &skip }, "only write out every nr-th row and column" },
1527 "insert line in [REF].xpm[ref] matrix where axis label is zero" },
1532 "Skip first N colors from [REF].xpm[ref] file for the legend" },
1533 { "-combine", FALSE, etENUM, { combine }, "Combine two matrices" },
1534 { "-cmin", FALSE, etREAL, { &cmin }, "Minimum for combination output" },
1535 { "-cmax", FALSE, etREAL, { &cmax }, "Maximum for combination output" }
1537 #define NPA asize(pa)
1538 t_filenm fnm[] = { { efXPM, "-f", nullptr, ffREAD }, { efXPM, "-f2", "root2", ffOPTRD },
1539 { efM2P, "-di", nullptr, ffLIBOPTRD }, { efM2P, "-do", "out", ffOPTWR },
1540 { efEPS, "-o", nullptr, ffOPTWR }, { efXPM, "-xpm", nullptr, ffOPTWR } };
1541 #define NFILE asize(fnm)
1543 if (!parse_common_args(
1544 &argc, argv, PCA_CAN_VIEW, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
1549 etitle = nenum(title);
1550 elegend = nenum(legend);
1551 ediag = nenum(diag);
1552 erainbow = nenum(rainbow);
1553 ecombine = nenum(combine);
1554 bGrad = opt2parg_bSet("-gradient", NPA, pa);
1555 for (i = 0; i < DIM; i++)
1557 if (grad[i] < 0 || grad[i] > 1)
1559 gmx_fatal(FARGS, "RGB value %g out of range (0.0-1.0)", grad[i]);
1568 epsfile = ftp2fn_null(efEPS, NFILE, fnm);
1569 xpmfile = opt2fn_null("-xpm", NFILE, fnm);
1570 if (epsfile == nullptr && xpmfile == nullptr)
1572 if (ecombine != ecHalves)
1574 xpmfile = opt2fn("-xpm", NFILE, fnm);
1578 epsfile = ftp2fn(efEPS, NFILE, fnm);
1581 if (ecombine != ecHalves && epsfile)
1584 "WARNING: can only write result of arithmetic combination "
1585 "of two matrices to .xpm file\n"
1586 " file %s will not be written\n",
1591 bDiag = ediag != edNone;
1592 bFirstDiag = ediag != edSecond;
1594 fn = opt2fn("-f", NFILE, fnm);
1595 std::vector<t_matrix> mat, mat2;
1596 mat = read_xpm_matrix(fn);
1598 "There %s %zu matri%s in %s\n",
1599 (mat.size() > 1) ? "are" : "is",
1601 (mat.size() > 1) ? "ces" : "x",
1603 fn = opt2fn_null("-f2", NFILE, fnm);
1606 mat2 = read_xpm_matrix(fn);
1608 "There %s %zu matri%s in %s\n",
1609 (mat2.size() > 1) ? "are" : "is",
1611 (mat2.size() > 1) ? "ces" : "x",
1613 if (mat.size() != mat2.size())
1615 fprintf(stderr, "Different number of matrices, using the smallest number.\n");
1616 if (mat.size() > mat2.size())
1618 mat.resize(mat2.size());
1622 mat2.resize(mat.size());
1628 if (ecombine != ecHalves)
1631 "WARNING: arithmetic matrix combination selected (-combine), "
1632 "but no second matrix (-f2) supplied\n"
1633 " no matrix combination will be performed\n");
1637 bTitle = etitle == etTop;
1638 bTitleOnce = etitle == etOnce;
1639 if (etitle == etYlabel)
1643 m.label_y = m.title;
1645 for (auto& m : mat2)
1647 m.label_y = m.title;
1652 gradient_mat(grad, mat);
1655 gradient_mat(grad, mat2);
1658 if (erainbow != erNo)
1660 rainbow_mat(erainbow == erBlue, mat);
1663 rainbow_mat(erainbow == erBlue, mat2);
1667 if (mat2.empty() && (elegend != elNone))
1672 if (ecombine && ecombine != ecHalves)
1674 write_combined_matrix(ecombine,
1678 opt2parg_bSet("-cmin", NPA, pa) ? &cmin : nullptr,
1679 opt2parg_bSet("-cmax", NPA, pa) ? &cmax : nullptr);
1698 opt2fn_null("-di", NFILE, fnm),
1699 opt2fn_null("-do", NFILE, fnm),
1704 view_all(oenv, NFILE, fnm);