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 static t_rgb black = { 0, 0, 0 };
211 static t_rgb white = { 1, 1, 1 };
212 #define BLACK (&black)
213 /* this must correspond to *colors[] in get_params */
214 static t_rgb* linecolors[] = { nullptr, &black, &white, nullptr };
216 static void leg_discrete(t_psdata* ps,
219 const std::string& label,
222 gmx::ArrayRef<const t_mapping> map)
227 boxhh = fontsize + DDD;
230 ps_strfont(ps, font, fontsize);
231 yhh = y0 + fontsize + 3 * DDD;
234 ps_ctext(ps, x0, yhh, label, eXLeft);
236 ps_moveto(ps, x0, y0);
237 for (const auto& m : map)
241 ps_fillbox(ps, DDD, DDD, DDD + fontsize, boxhh - DDD);
243 ps_box(ps, DDD, DDD, DDD + fontsize, boxhh - DDD);
244 ps_ctext(ps, boxhh + 2 * DDD, fontsize / 3, m.desc, eXLeft);
246 ps_moverel(ps, DDD, -fontsize / 3);
250 static void leg_continuous(t_psdata* ps,
254 const std::string& label,
257 gmx::ArrayRef<const t_mapping> map,
261 real yhh, boxxh, boxyh;
262 gmx::index mapIndex = gmx::ssize(map) - mapoffset;
265 if (x < 8 * fontsize)
269 boxxh = x / mapIndex;
270 if (boxxh > fontsize)
275 GMX_RELEASE_ASSERT(!map.empty(), "NULL map array provided to leg_continuous()");
278 xx0 = x0 - (mapIndex * boxxh) / 2.0;
280 for (gmx::index i = 0; (i < mapIndex); i++)
282 ps_rgb(ps, &(map[i + mapoffset].rgb));
283 ps_fillbox(ps, xx0 + i * boxxh, y0, xx0 + (i + 1) * boxxh, y0 + boxyh);
285 ps_strfont(ps, font, fontsize);
287 ps_box(ps, xx0, y0, xx0 + mapIndex * boxxh, y0 + boxyh);
289 yhh = y0 + boxyh + 3 * DDD;
290 ps_ctext(ps, xx0 + boxxh / 2, yhh, map[0].desc, eXCenter);
293 ps_ctext(ps, x0, yhh, label, eXCenter);
295 ps_ctext(ps, xx0 + (mapIndex * boxxh) - boxxh / 2, yhh, map[map.size() - 1].desc, eXCenter);
298 static void leg_bicontinuous(t_psdata* ps,
302 const std::string& label1,
303 const std::string& label2,
306 gmx::ArrayRef<const t_mapping> map1,
307 gmx::ArrayRef<const t_mapping> map2)
309 real xx1, xx2, x1, x2;
311 x1 = x / (map1.size() + map2.size()) * map1.size(); /* width of legend 1 */
312 x2 = x / (map1.size() + map2.size()) * map2.size(); /* width of legend 2 */
313 xx1 = x0 - (x2 / 2.0) - fontsize; /* center of legend 1 */
314 xx2 = x0 + (x1 / 2.0) + fontsize; /* center of legend 2 */
315 x1 -= fontsize / 2; /* adjust width */
316 x2 -= fontsize / 2; /* adjust width */
317 leg_continuous(ps, xx1, x1, y0, label1, fontsize, font, map1, 0);
318 leg_continuous(ps, xx2, x2, y0, label2, fontsize, font, map2, 0);
321 static real box_height(const t_matrix& mat, t_psrec* psr)
323 return mat.ny * psr->yboxsize;
326 static real box_dh(t_psrec* psr)
328 return psr->boxspacing;
331 static real box_dh_top(gmx_bool bOnce, t_psrec* psr)
335 if (psr->bTitle || (psr->bTitleOnce && bOnce))
337 dh = 2 * psr->titfontsize;
347 static gmx_bool box_do_all_x_maj_ticks(t_psrec* psr)
349 return (psr->boxspacing > (1.5 * psr->X.majorticklen));
352 static gmx_bool box_do_all_x_min_ticks(t_psrec* psr)
354 return (psr->boxspacing > (1.5 * psr->X.minorticklen));
357 static void draw_boxes(t_psdata* ps, real x0, real y0, real w, gmx::ArrayRef<t_matrix> mat, t_psrec* psr)
361 char **xtick, **ytick;
362 real xx, yy, dy, xx00, yy00, offset_x, offset_y;
366 /* Only necessary when there will be no y-labels */
371 ps_linewidth(ps, static_cast<int>(psr->boxlinewidth));
373 for (auto m = mat.begin(); m != mat.end(); ++m)
375 dy = box_height(*m, psr);
376 ps_box(ps, x0 - 1, yy00 - 1, x0 + w + 1, yy00 + dy + 1);
377 yy00 += dy + box_dh(psr) + box_dh_top(m + 1 == mat.end(), psr);
380 /* Draw the ticks on the axes */
381 ps_linewidth(ps, static_cast<int>(psr->ticklinewidth));
384 auto halfway = mat.begin() + (mat.size() / 2);
385 for (auto m = mat.begin(); m != mat.end(); ++m)
387 if (m->flags & MAT_SPATIAL_X)
397 if (m->flags & MAT_SPATIAL_Y)
408 for (int j = 0; (j < ntx); j++)
410 sprintf(buf, "%g", m->axis_x[j]);
411 xtick[j] = gmx_strdup(buf);
413 ps_strfont(ps, psr->X.tickfont, psr->X.tickfontsize);
414 for (x = 0; (x < ntx); x++)
416 xx = xx00 + (x + offset_x) * psr->xboxsize;
417 if ((bRmod(m->axis_x[x], psr->X.offset, psr->X.major) || (psr->X.first && (x == 0)))
418 && (m == mat.begin() || box_do_all_x_maj_ticks(psr)))
420 /* Longer tick marks */
421 ps_line(ps, xx, yy00, xx, yy00 - psr->X.majorticklen);
422 /* Plot label on lowest graph only */
423 if (m == mat.begin())
425 ps_ctext(ps, xx, yy00 - DDD - psr->X.majorticklen - psr->X.tickfontsize * 0.8, xtick[x], eXCenter);
428 else if (bRmod(m->axis_x[x], psr->X.offset, psr->X.minor)
429 && ((m == mat.begin()) || box_do_all_x_min_ticks(psr)))
431 /* Shorter tick marks */
432 ps_line(ps, xx, yy00, xx, yy00 - psr->X.minorticklen);
434 else if (bRmod(m->axis_x[x], psr->X.offset, psr->X.major))
436 /* Even shorter marks, only each X.major */
437 ps_line(ps, xx, yy00, xx, yy00 - (psr->boxspacing / 2));
440 ps_strfont(ps, psr->Y.tickfont, psr->Y.tickfontsize);
442 for (int j = 0; (j < nty); j++)
444 sprintf(buf, "%g", m->axis_y[j]);
445 ytick[j] = gmx_strdup(buf);
448 for (int y = 0; (y < nty); y++)
450 yy = yy00 + (y + offset_y) * psr->yboxsize;
451 if (bRmod(m->axis_y[y], psr->Y.offset, psr->Y.major) || (psr->Y.first && (y == 0)))
454 strlength = std::max(strlength, std::strlen(ytick[y]));
455 ps_line(ps, xx00, yy, xx00 - psr->Y.majorticklen, yy);
456 ps_ctext(ps, xx00 - psr->Y.majorticklen - DDD, yy - psr->Y.tickfontsize / 3.0, ytick[y], eXRight);
458 else if (bRmod(m->axis_y[y], psr->Y.offset, psr->Y.minor))
461 ps_line(ps, xx00, yy, xx00 - psr->Y.minorticklen, yy);
467 /* Label on Y-axis */
468 if (!psr->bYonce || m == halfway)
471 if (strlen(psr->Y.label) > 0)
473 mylab = psr->Y.label;
481 ps_strfont(ps, psr->Y.font, psr->Y.fontsize);
483 xxx = x0 - psr->X.majorticklen - psr->X.tickfontsize * strlength - DDD;
484 ps_ctext(ps, yy00 + box_height(*m, psr) / 2.0, 612.5 - xxx, mylab, eXCenter);
489 yy00 += box_height(*m, psr) + box_dh(psr) + box_dh_top(m + 1 == mat.end(), psr);
491 /* Label on X-axis */
493 if (strlen(psr->X.label) > 0)
495 mylab = psr->X.label;
499 mylab = mat[0].label_x;
503 ps_strfont(ps, psr->X.font, psr->X.fontsize);
506 y0 - DDD - psr->X.majorticklen - psr->X.tickfontsize * FUDGE - psr->X.fontsize,
512 static void draw_zerolines(t_psdata* out, real x0, real y0, real w, gmx::ArrayRef<t_matrix> mat, t_psrec* psr)
514 real xx, yy, dy, xx00, yy00;
519 ps_linewidth(out, static_cast<int>(psr->zerolinewidth));
520 for (auto m = mat.begin(); m != mat.end(); ++m)
522 dy = box_height(*m, psr);
523 /* m->axis_x and _y were already set by draw_boxes */
524 if (psr->X.lineatzero)
526 ps_rgb(out, linecolors[psr->X.lineatzero]);
527 for (x = 0; (x < m->nx); x++)
529 xx = xx00 + (x + 0.7) * psr->xboxsize;
530 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
531 if (x != 0 && x < m->nx - 1
532 && std::abs(m->axis_x[x]) < 0.1 * std::abs(m->axis_x[x + 1] - m->axis_x[x]))
534 ps_line(out, xx, yy00, xx, yy00 + dy + 2);
538 if (psr->Y.lineatzero)
540 ps_rgb(out, linecolors[psr->Y.lineatzero]);
541 for (y = 0; (y < m->ny); y++)
543 yy = yy00 + (y + 0.7) * psr->yboxsize;
544 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
545 if (y != 0 && y < m->ny - 1
546 && std::abs(m->axis_y[y]) < 0.1 * std::abs(m->axis_y[y + 1] - m->axis_y[y]))
548 ps_line(out, xx00, yy, xx00 + w + 2, yy);
552 yy00 += box_height(*m, psr) + box_dh(psr) + box_dh_top(m + 1 == mat.end(), psr);
556 static void box_dim(gmx::ArrayRef<t_matrix> mat,
557 gmx::ArrayRef<t_matrix> mat2,
567 real ww, hh, dww, dhh;
573 for (const auto& m : mat)
575 ww = std::max(ww, m.nx * psr->xboxsize);
576 hh += box_height(m, psr);
577 maxytick = std::max(maxytick, m.nx);
581 if (mat[0].label_y[0])
583 dww += 2.0 * (psr->Y.fontsize + DDD);
585 if (psr->Y.major > 0)
587 dww += psr->Y.majorticklen + DDD
588 + psr->Y.tickfontsize * (std::log(static_cast<real>(maxytick)) / std::log(10.0));
590 else if (psr->Y.minor > 0)
592 dww += psr->Y.minorticklen;
595 if (mat[0].label_x[0])
597 dhh += psr->X.fontsize + 2 * DDD;
599 if (/* fool emacs auto-indent */
600 (elegend == elBoth && (mat[0].legend[0] || (!mat2.empty() && mat2[0].legend[0])))
601 || (elegend == elFirst && mat[0].legend[0])
602 || (elegend == elSecond && (!mat2.empty() && mat2[0].legend[0])))
604 dhh += 2 * (psr->legfontsize * FUDGE + 2 * DDD);
608 dhh += psr->legfontsize * FUDGE + 2 * DDD;
610 if (psr->X.major > 0)
612 dhh += psr->X.tickfontsize * FUDGE + 2 * DDD + psr->X.majorticklen;
614 else if (psr->X.minor > 0)
616 dhh += psr->X.minorticklen;
619 hh += (mat.size() - 1) * box_dh(psr);
620 hh += box_dh_top(TRUE, psr);
623 hh += (mat.size() - 1) * box_dh_top(FALSE, psr);
632 static std::vector<t_mapping> add_maps(gmx::ArrayRef<t_mapping> map1, gmx::ArrayRef<t_mapping> map2)
634 static char mapper[] =
635 "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<"
637 std::vector<t_mapping> map(map1.size() + map2.size());
639 size_t nsymbols = std::strlen(mapper);
640 if (map.size() > nsymbols * nsymbols)
642 gmx_fatal(FARGS, "Not enough symbols to merge the two colormaps\n");
644 printf("Combining colormaps of %zu and %zu elements into one of %zu elements\n",
649 for (gmx::index j = 0; j < gmx::ssize(map1) && k < gmx::ssize(map); ++j, ++k)
651 map[k].code.c1 = mapper[k % nsymbols];
652 if (map.size() > nsymbols)
654 map[k].code.c2 = mapper[k / nsymbols];
656 map[k].rgb.r = map1[j].rgb.r;
657 map[k].rgb.g = map1[j].rgb.g;
658 map[k].rgb.b = map1[j].rgb.b;
659 map[k].desc = map1[j].desc;
661 for (gmx::index j = 0; j < gmx::ssize(map2) && k < gmx::ssize(map); ++j, ++k)
663 map[k].code.c1 = mapper[k % nsymbols];
664 if (map.size() > nsymbols)
666 map[k].code.c2 = mapper[k / nsymbols];
668 map[k].rgb.r = map2[j].rgb.r;
669 map[k].rgb.g = map2[j].rgb.g;
670 map[k].rgb.b = map2[j].rgb.b;
671 map[k].desc = map2[j].desc;
677 static bool operator==(const t_mapping& lhs, const t_mapping& rhs)
679 return (lhs.rgb.r == rhs.rgb.r && lhs.rgb.g == rhs.rgb.g && lhs.rgb.b == rhs.rgb.b);
682 static void xpm_mat(const char* outf,
683 gmx::ArrayRef<t_matrix> mat,
684 gmx::ArrayRef<t_matrix> mat2,
691 out = gmx_ffopen(outf, "w");
693 GMX_RELEASE_ASSERT(mat.size() == mat2.size(),
694 "Combined matrix write requires matrices of the same size");
695 for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
697 // Color maps that differ only in RGB value are considered different
698 if (mat2.empty() || std::equal(mat[i].map.begin(), mat[i].map.end(), mat2[i].map.begin()))
700 write_xpm_m(out, mat[0]);
704 auto map = add_maps(mat[i].map, mat2[i].map);
705 for (x = 0; (x < mat[i].nx); x++)
707 for (y = 0; (y < mat[i].nx); y++)
709 if ((x < y) || ((x == y) && bFirstDiag)) /* upper left -> map1 */
711 col = mat[i].matrix(x, y);
713 else /* lower right -> map2 */
715 col = mat[i].map.size() + mat[i].matrix(x, y);
717 if ((bDiag) || (x != y))
719 mat[i].matrix(x, y) = col;
723 mat[i].matrix(x, y) = 0;
728 if (mat[i].title != mat2[i].title)
730 mat[i].title += " / " + mat2[i].title;
732 if (mat[i].legend != mat2[i].legend)
734 mat[i].legend += " / " + mat2[i].legend;
736 write_xpm_m(out, mat[i]);
742 static void tick_spacing(int n, real axis[], real offset, char axisnm, real* major, real* minor)
746 int i, j, t, f = 0, ten;
748 real major_fact[NFACT] = { 5, 4, 2, 1 };
749 real minor_fact[NFACT] = { 5, 4, 4, 5 };
751 /* start with interval between 10 matrix points: */
752 space = std::max(10 * axis[1] - axis[0], axis[std::min(10, n - 1)] - axis[0]);
753 /* get power of 10 */
754 ten = static_cast<int>(std::ceil(std::log(space) / std::log(10.0)) - 1);
756 for (t = ten + 2; t > ten - 3 && bTryAgain; t--)
758 for (f = 0; f < NFACT && bTryAgain; f++)
760 space = std::pow(10.0_real, static_cast<real>(t)) * major_fact[f];
761 /* count how many ticks we would get: */
763 for (j = 0; j < n; j++)
765 if (bRmod(axis[j], offset, space))
770 /* do we have a reasonable number of ticks ? */
771 bTryAgain = (i > std::min(10, n - 1)) || (i < 5);
776 space = std::max(10 * axis[1] - axis[0], axis[std::min(10, n - 1)] - axis[0]);
777 fprintf(stderr, "Auto tick spacing failed for %c-axis, guessing %g\n", axisnm, space);
780 *minor = space / minor_fact[(f > 0) ? f - 1 : 0];
781 fprintf(stderr, "Auto tick spacing for %c-axis: major %g, minor %g\n", axisnm, *major, *minor);
784 static void ps_mat(const char* outf,
785 gmx::ArrayRef<t_matrix> mat,
786 gmx::ArrayRef<t_matrix> mat2,
806 gmx_bool bMap1, bNextMap1, bDiscrete;
810 get_params(m2p, m2pout, &psrec);
812 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
816 if (psr->X.major <= 0)
818 tick_spacing((mat[0].flags & MAT_SPATIAL_X) ? mat[0].nx + 1 : mat[0].nx,
819 mat[0].axis_x.data(),
825 if (psr->X.minor <= 0)
827 psr->X.minor = psr->X.major / 2;
829 if (psr->Y.major <= 0)
831 tick_spacing((mat[0].flags & MAT_SPATIAL_Y) ? mat[0].ny + 1 : mat[0].ny,
832 mat[0].axis_y.data(),
838 if (psr->Y.minor <= 0)
840 psr->Y.minor = psr->Y.major / 2;
845 psr->xboxsize = boxx;
846 psr->yboxsize = boxx;
850 psr->yboxsize = boxy;
853 if (psr->xboxsize == 0)
855 psr->xboxsize = size / mat[0].nx;
856 printf("Set the x-size of the box to %.3f\n", psr->xboxsize);
858 if (psr->yboxsize == 0)
860 psr->yboxsize = size / mat[0].nx;
861 printf("Set the y-size of the box to %.3f\n", psr->yboxsize);
864 gmx::ArrayRef<const t_mapping> map1;
866 for (const auto& m : mat)
868 if (m.map.size() > map1.size())
872 printf("Selected legend of matrix # %d for display\n", legendIndex);
878 gmx::ArrayRef<const t_mapping> map2;
881 for (const auto& m : mat2)
883 if (m.map.size() > map2.size())
887 printf("Selected legend of matrix # %d for second display\n", legendIndex);
894 if (mat[0].legend.empty() && psr->legend)
896 mat[0].legend = psr->leglabel;
899 bTitle = bTitle && !mat.back().title.empty();
900 bTitleOnce = bTitleOnce && !mat.back().title.empty();
901 psr->bTitle = bTitle;
902 psr->bTitleOnce = bTitleOnce;
903 psr->bYonce = bYonce;
905 /* Set up size of box for nice colors */
906 box_dim(mat, mat2, psr, elegend, bFrame, &w, &h, &dw, &dh);
908 /* Set up bounding box */
909 W = static_cast<int>(w + dw);
910 H = static_cast<int>(h + dh);
915 x = static_cast<int>(W + psr->xoffs);
916 y = static_cast<int>(H + psr->yoffs);
922 t_psdata out = ps_open(outf, 0, 0, x, y);
923 ps_linewidth(&out, static_cast<int>(psr->linewidth));
924 ps_init_rgb_box(&out, psr->xboxsize, psr->yboxsize);
925 ps_init_rgb_nbox(&out, psr->xboxsize, psr->yboxsize);
926 ps_translate(&out, psr->xoffs, psr->yoffs);
930 ps_comment(&out, "Here starts the BOX drawing");
931 draw_boxes(&out, x0, y0, w, mat, psr);
934 for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
936 if (bTitle || (bTitleOnce && i == gmx::ssize(mat) - 1))
938 /* Print title, if any */
940 ps_strfont(&out, psr->titfont, psr->titfontsize);
942 if (!mat2.empty() || mat[i].title == mat2[i].title)
948 buf = mat[i].title + " / " + mat2[i].title;
950 ps_ctext(&out, x0 + w / 2, y0 + box_height(mat[i], psr) + psr->titfontsize, buf, eXCenter);
952 ps_comment(&out, gmx::formatString("Here starts the filling of box #%zd", i).c_str());
953 for (x = 0; (x < mat[i].nx); x++)
958 xx = x0 + x * psr->xboxsize;
959 ps_moveto(&out, xx, y0);
961 bMap1 = (mat2.empty() || (x < y || (x == y && bFirstDiag)));
962 if ((bDiag) || (x != y))
964 col = mat[i].matrix(x, y);
970 for (nexty = 1; (nexty <= mat[i].ny); nexty++)
972 bNextMap1 = (mat2.empty() || (x < nexty || (x == nexty && bFirstDiag)));
973 /* TRUE: upper left -> map1 */
974 /* FALSE: lower right -> map2 */
975 if ((nexty == mat[i].ny) || (!bDiag && (x == nexty)))
981 nextcol = mat[i].matrix(x, nexty);
983 if ((nexty == mat[i].ny) || (col != nextcol) || (bMap1 != bNextMap1))
989 ps_rgb_nbox(&out, &(mat[i].map[col].rgb), nexty - y);
993 assert(!mat2.empty());
994 ps_rgb_nbox(&out, &(mat2[i].map[col].rgb), nexty - y);
999 ps_moverel(&out, 0, psr->yboxsize);
1007 y0 += box_height(mat[i], psr) + box_dh(psr) + box_dh_top(i + 1 == gmx::ssize(mat), psr);
1010 if (psr->X.lineatzero || psr->Y.lineatzero)
1012 /* reset y0 for first box */
1014 ps_comment(&out, "Here starts the zero lines drawing");
1015 draw_zerolines(&out, x0, y0, w, mat, psr);
1018 if (elegend != elNone)
1021 gmx::ArrayRef<const t_mapping> leg_map;
1022 ps_comment(&out, "Now it's legend time!");
1023 ps_linewidth(&out, static_cast<int>(psr->linewidth));
1024 if (mat2.empty() || elegend != elSecond)
1026 bDiscrete = mat[0].bDiscrete;
1027 legend = mat[0].legend;
1032 bDiscrete = mat2[0].bDiscrete;
1033 legend = mat2[0].legend;
1038 leg_discrete(&out, psr->legfontsize, DDD, legend, psr->legfontsize, psr->legfont, leg_map);
1042 if (elegend != elBoth)
1045 &out, x0 + w / 2, w / 2, DDD, legend, psr->legfontsize, psr->legfont, leg_map, mapoffset);
1049 assert(!mat2.empty());
1051 &out, x0 + w / 2, w, DDD, mat[0].legend, mat2[0].legend, psr->legfontsize, psr->legfont, map1, map2);
1054 ps_comment(&out, "Done processing");
1060 static void make_axis_labels(gmx::ArrayRef<t_matrix> mat1)
1062 for (auto& m : mat1)
1064 /* Make labels for x axis */
1065 if (m.axis_x.empty())
1067 m.axis_x.resize(m.nx);
1068 std::iota(m.axis_x.begin(), m.axis_x.end(), 0);
1070 /* Make labels for y axis */
1071 if (m.axis_y.empty())
1073 m.axis_y.resize(m.ny);
1074 std::iota(m.axis_y.begin(), m.axis_y.end(), 0);
1079 static void prune_mat(gmx::ArrayRef<t_matrix> mat, gmx::ArrayRef<t_matrix> mat2, int skip)
1081 GMX_RELEASE_ASSERT(mat.size() == mat2.size() || mat2.empty(),
1082 "Matrix pruning requires matrices of the same size");
1083 for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
1086 "converting %dx%d matrix to %dx%d\n",
1089 (mat[i].nx + skip - 1) / skip,
1090 (mat[i].ny + skip - 1) / skip);
1091 /* walk through matrix */
1093 for (int x = 0; (x < mat[i].nx); x++)
1097 mat[i].axis_x[xs] = mat[i].axis_x[x];
1100 mat2[i].axis_x[xs] = mat2[i].axis_x[x];
1103 for (int y = 0; (y < mat[i].ny); y++)
1107 mat[i].axis_y[ys] = mat[i].axis_y[y];
1110 mat2[i].axis_y[ys] = mat2[i].axis_y[y];
1115 mat[i].matrix(xs, ys) = mat[i].matrix(x, y);
1118 mat2[i].matrix(xs, ys) = mat2[i].matrix(x, y);
1126 /* adjust parameters */
1127 mat[i].nx = (mat[i].nx + skip - 1) / skip;
1128 mat[i].ny = (mat[i].ny + skip - 1) / skip;
1131 mat2[i].nx = (mat2[i].nx + skip - 1) / skip;
1132 mat2[i].ny = (mat2[i].ny + skip - 1) / skip;
1137 static void zero_lines(gmx::ArrayRef<t_matrix> mat, gmx::ArrayRef<t_matrix> mat2)
1139 GMX_RELEASE_ASSERT(mat.size() == mat2.size(), "zero_lines requires matrices of the same size");
1140 for (gmx::index i = 0; i != gmx::ssize(mat); ++i)
1142 for (int m = 0; m < (!mat2.empty() ? 2 : 1); m++)
1144 gmx::ArrayRef<t_matrix> mats;
1153 for (int x = 0; x < mats[i].nx - 1; x++)
1155 if (std::abs(mats[i].axis_x[x + 1]) < 1e-5)
1157 for (int y = 0; y < mats[i].ny; y++)
1159 mats[i].matrix(x, y) = 0;
1163 for (int y = 0; y < mats[i].ny - 1; y++)
1165 if (std::abs(mats[i].axis_y[y + 1]) < 1e-5)
1167 for (int x = 0; x < mats[i].nx; x++)
1169 mats[i].matrix(x, y) = 0;
1177 static void write_combined_matrix(int ecombine,
1179 gmx::ArrayRef<t_matrix> mat1,
1180 gmx::ArrayRef<t_matrix> mat2,
1185 real **rmat1, **rmat2;
1188 out = gmx_ffopen(fn, "w");
1189 GMX_RELEASE_ASSERT(mat1.size() == mat2.size(),
1190 "Combined matrix write requires matrices of the same size");
1191 for (gmx::index k = 0; k != gmx::ssize(mat1); k++)
1193 if (mat2[k].nx != mat1[k].nx || mat2[k].ny != mat1[k].ny)
1196 "Size of frame %zd in 1st (%dx%d) and 2nd matrix (%dx%d) do"
1204 printf("Combining two %dx%d matrices\n", mat1[k].nx, mat1[k].ny);
1205 rmat1 = matrix2real(&mat1[k], nullptr);
1206 rmat2 = matrix2real(&mat2[k], nullptr);
1207 if (nullptr == rmat1 || nullptr == rmat2)
1210 "Could not extract real data from %s xpm matrices. Note that, e.g.,\n"
1211 "g_rms and g_mdmat provide such data, but not do_dssp.\n",
1212 (nullptr == rmat1 && nullptr == rmat2) ? "both" : "one of the");
1216 for (int j = 0; j < mat1[k].ny; j++)
1218 for (int i = 0; i < mat1[k].nx; i++)
1222 case ecAdd: rmat1[i][j] += rmat2[i][j]; break;
1223 case ecSub: rmat1[i][j] -= rmat2[i][j]; break;
1224 case ecMult: rmat1[i][j] *= rmat2[i][j]; break;
1225 case ecDiv: rmat1[i][j] /= rmat2[i][j]; break;
1226 default: gmx_fatal(FARGS, "No such combination rule %d for matrices", ecombine);
1228 rlo = std::min(rlo, rmat1[i][j]);
1229 rhi = std::max(rhi, rmat1[i][j]);
1240 int nlevels = static_cast<int>(std::max(mat1[k].map.size(), mat2[k].map.size()));
1243 fprintf(stderr, "combination results in uniform matrix (%g), no output\n", rhi);
1255 mat1[k].axis_x.data(),
1256 mat1[k].axis_y.data(),
1268 static void do_mat(gmx::ArrayRef<t_matrix> mat,
1269 gmx::ArrayRef<t_matrix> mat2,
1273 gmx_bool bFirstDiag,
1275 gmx_bool bTitleOnce,
1281 const char* epsfile,
1282 const char* xpmfile,
1288 GMX_RELEASE_ASSERT(mat.size() == mat2.size() || mat2.empty(),
1289 "Combined matrix write requires matrices of the same size");
1292 for (gmx::index k = 0; k != gmx::ssize(mat); k++)
1294 if ((mat2[k].nx != mat[k].nx) || (mat2[k].ny != mat[k].ny))
1297 "WAKE UP!! Size of frame %zd in 2nd matrix file (%dx%d) does not match "
1298 "size of 1st matrix (%dx%d) or the other way around.\n",
1305 for (int j = 0; (j < mat[k].ny); j++)
1307 for (int i = bFirstDiag ? j + 1 : j; (i < mat[k].nx); i++)
1309 mat[k].matrix(i, j) = mat2[k].matrix(i, j);
1314 for (gmx::index i = 0; i != gmx::ssize(mat); i++)
1316 fprintf(stderr, "Matrix %zd is %d x %d\n", i, mat[i].nx, mat[i].ny);
1319 make_axis_labels(mat);
1323 prune_mat(mat, mat2, skip);
1328 zero_lines(mat, mat);
1331 if (epsfile != nullptr)
1333 ps_mat(epsfile, mat, mat2, bFrame, bDiag, bFirstDiag, bTitle, bTitleOnce, bYonce, elegend, size, boxx, boxy, m2p, m2pout, mapoffset);
1335 if (xpmfile != nullptr)
1337 xpm_mat(xpmfile, mat, mat2, bDiag, bFirstDiag);
1341 static void gradient_map(const rvec grad, gmx::ArrayRef<t_mapping> map)
1344 real fraction = 1.0 / (map.size() - 1.0);
1347 real c = i * fraction;
1348 m.rgb.r = 1 - c * (1 - grad[XX]);
1349 m.rgb.g = 1 - c * (1 - grad[YY]);
1350 m.rgb.b = 1 - c * (1 - grad[ZZ]);
1355 static void gradient_mat(rvec grad, gmx::ArrayRef<t_matrix> mat)
1359 gradient_map(grad, m.map);
1363 static void rainbow_map(gmx_bool bBlue, gmx::ArrayRef<t_mapping> map)
1367 real c = (m.rgb.r + m.rgb.g + m.rgb.b) / 3;
1377 if (c <= 0.25) /* 0-0.25 */
1380 g = std::pow(4.0 * c, 2.0 / 3.0);
1383 else if (c <= 0.5) /* 0.25-0.5 */
1387 b = std::pow(2.0 - 4.0 * c, 2.0 / 3.0);
1389 else if (c <= 0.75) /* 0.5-0.75 */
1391 r = std::pow(4.0 * c - 2.0, 2.0 / 3.0);
1398 g = std::pow(4.0 - 4.0 * c, 2.0 / 3.0);
1407 static void rainbow_mat(gmx_bool bBlue, gmx::ArrayRef<t_matrix> mat)
1411 rainbow_map(bBlue, m.map);
1415 int gmx_xpm2ps(int argc, char* argv[])
1417 const char* desc[] = {
1418 "[THISMODULE] makes a beautiful color plot of an XPixelMap file.",
1419 "Labels and axis can be displayed, when they are supplied",
1420 "in the correct matrix format.",
1421 "Matrix data may be generated by programs such as [gmx-do_dssp], [gmx-rms] or",
1422 "[gmx-mdmat].[PAR]",
1423 "Parameters are set in the [TT].m2p[tt] file optionally supplied with",
1424 "[TT]-di[tt]. Reasonable defaults are provided. Settings for the [IT]y[it]-axis",
1425 "default to those for the [IT]x[it]-axis. Font names have a defaulting hierarchy:",
1426 "titlefont -> legendfont; titlefont -> (xfont -> yfont -> ytickfont)",
1427 "-> xtickfont, e.g. setting titlefont sets all fonts, setting xfont",
1428 "sets yfont, ytickfont and xtickfont.[PAR]",
1429 "When no [TT].m2p[tt] file is supplied, many settings are taken from",
1430 "command line options. The most important option is [TT]-size[tt],",
1431 "which sets the size of the whole matrix in postscript units.",
1432 "This option can be overridden with the [TT]-bx[tt] and [TT]-by[tt]",
1433 "options (and the corresponding parameters in the [TT].m2p[tt] file),",
1434 "which set the size of a single matrix element.[PAR]",
1435 "With [TT]-f2[tt] a second matrix file can be supplied. Both matrix",
1436 "files will be read simultaneously and the upper left half of the",
1437 "first one ([TT]-f[tt]) is plotted together with the lower right",
1438 "half of the second one ([TT]-f2[tt]). The diagonal will contain",
1439 "values from the matrix file selected with [TT]-diag[tt].",
1440 "Plotting of the diagonal values can be suppressed altogether by",
1441 "setting [TT]-diag[tt] to [TT]none[tt].",
1442 "In this case, a new color map will be generated with",
1443 "a red gradient for negative numbers and a blue for positive.",
1444 "If the color coding and legend labels of both matrices are identical,",
1445 "only one legend will be displayed, else two separate legends are",
1447 "With [TT]-combine[tt], an alternative operation can be selected",
1448 "to combine the matrices. The output range is automatically set",
1449 "to the actual range of the combined matrix. This can be overridden",
1450 "with [TT]-cmin[tt] and [TT]-cmax[tt].[PAR]",
1451 "[TT]-title[tt] can be set to [TT]none[tt] to suppress the title, or to",
1452 "[TT]ylabel[tt] to show the title in the Y-label position (alongside",
1453 "the [IT]y[it]-axis).[PAR]",
1454 "With the [TT]-rainbow[tt] option, dull grayscale matrices can be turned",
1455 "into attractive color pictures.[PAR]",
1456 "Merged or rainbowed matrices can be written to an XPixelMap file with",
1457 "the [TT]-xpm[tt] option."
1460 gmx_output_env_t* oenv;
1461 const char * fn, *epsfile = nullptr, *xpmfile = nullptr;
1462 int i, etitle, elegend, ediag, erainbow, ecombine;
1463 gmx_bool bTitle, bTitleOnce, bDiag, bFirstDiag, bGrad;
1464 static gmx_bool bFrame = TRUE, bZeroLine = FALSE, bYonce = FALSE;
1465 static real size = 400, boxx = 0, boxy = 0, cmin = 0, cmax = 0;
1466 static rvec grad = { 0, 0, 0 };
1476 const char* title[] = { nullptr, "top", "once", "ylabel", "none", nullptr };
1477 /* MUST correspond to enum elXxx as defined at top of file */
1478 const char* legend[] = { nullptr, "both", "first", "second", "none", nullptr };
1487 const char* diag[] = { nullptr, "first", "second", "none", nullptr };
1496 const char* rainbow[] = { nullptr, "no", "blue", "red", nullptr };
1497 /* MUST correspond to enum ecXxx as defined at top of file */
1498 const char* combine[] = { nullptr, "halves", "add", "sub", "mult", "div", nullptr };
1499 static int skip = 1, mapoffset = 0;
1501 { "-frame", FALSE, etBOOL, { &bFrame }, "Display frame, ticks, labels, title and legend" },
1502 { "-title", FALSE, etENUM, { title }, "Show title at" },
1503 { "-yonce", FALSE, etBOOL, { &bYonce }, "Show y-label only once" },
1504 { "-legend", FALSE, etENUM, { legend }, "Show legend" },
1505 { "-diag", FALSE, etENUM, { diag }, "Diagonal" },
1506 { "-size", FALSE, etREAL, { &size }, "Horizontal size of the matrix in ps units" },
1511 "Element x-size, overrides [TT]-size[tt] (also y-size when [TT]-by[tt] is not set)" },
1512 { "-by", FALSE, etREAL, { &boxy }, "Element y-size" },
1513 { "-rainbow", FALSE, etENUM, { rainbow }, "Rainbow colors, convert white to" },
1518 "Re-scale colormap to a smooth gradient from white {1,1,1} to {r,g,b}" },
1519 { "-skip", FALSE, etINT, { &skip }, "only write out every nr-th row and column" },
1524 "insert line in [REF].xpm[ref] matrix where axis label is zero" },
1529 "Skip first N colors from [REF].xpm[ref] file for the legend" },
1530 { "-combine", FALSE, etENUM, { combine }, "Combine two matrices" },
1531 { "-cmin", FALSE, etREAL, { &cmin }, "Minimum for combination output" },
1532 { "-cmax", FALSE, etREAL, { &cmax }, "Maximum for combination output" }
1534 #define NPA asize(pa)
1535 t_filenm fnm[] = { { efXPM, "-f", nullptr, ffREAD }, { efXPM, "-f2", "root2", ffOPTRD },
1536 { efM2P, "-di", nullptr, ffLIBOPTRD }, { efM2P, "-do", "out", ffOPTWR },
1537 { efEPS, "-o", nullptr, ffOPTWR }, { efXPM, "-xpm", nullptr, ffOPTWR } };
1538 #define NFILE asize(fnm)
1540 if (!parse_common_args(
1541 &argc, argv, PCA_CAN_VIEW, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
1546 etitle = nenum(title);
1547 elegend = nenum(legend);
1548 ediag = nenum(diag);
1549 erainbow = nenum(rainbow);
1550 ecombine = nenum(combine);
1551 bGrad = opt2parg_bSet("-gradient", NPA, pa);
1552 for (i = 0; i < DIM; i++)
1554 if (grad[i] < 0 || grad[i] > 1)
1556 gmx_fatal(FARGS, "RGB value %g out of range (0.0-1.0)", grad[i]);
1565 epsfile = ftp2fn_null(efEPS, NFILE, fnm);
1566 xpmfile = opt2fn_null("-xpm", NFILE, fnm);
1567 if (epsfile == nullptr && xpmfile == nullptr)
1569 if (ecombine != ecHalves)
1571 xpmfile = opt2fn("-xpm", NFILE, fnm);
1575 epsfile = ftp2fn(efEPS, NFILE, fnm);
1578 if (ecombine != ecHalves && epsfile)
1581 "WARNING: can only write result of arithmetic combination "
1582 "of two matrices to .xpm file\n"
1583 " file %s will not be written\n",
1588 bDiag = ediag != edNone;
1589 bFirstDiag = ediag != edSecond;
1591 fn = opt2fn("-f", NFILE, fnm);
1592 std::vector<t_matrix> mat, mat2;
1593 mat = read_xpm_matrix(fn);
1595 "There %s %zu matri%s in %s\n",
1596 (mat.size() > 1) ? "are" : "is",
1598 (mat.size() > 1) ? "ces" : "x",
1600 fn = opt2fn_null("-f2", NFILE, fnm);
1603 mat2 = read_xpm_matrix(fn);
1605 "There %s %zu matri%s in %s\n",
1606 (mat2.size() > 1) ? "are" : "is",
1608 (mat2.size() > 1) ? "ces" : "x",
1610 if (mat.size() != mat2.size())
1612 fprintf(stderr, "Different number of matrices, using the smallest number.\n");
1613 if (mat.size() > mat2.size())
1615 mat.resize(mat2.size());
1619 mat2.resize(mat.size());
1625 if (ecombine != ecHalves)
1628 "WARNING: arithmetic matrix combination selected (-combine), "
1629 "but no second matrix (-f2) supplied\n"
1630 " no matrix combination will be performed\n");
1634 bTitle = etitle == etTop;
1635 bTitleOnce = etitle == etOnce;
1636 if (etitle == etYlabel)
1640 m.label_y = m.title;
1642 for (auto& m : mat2)
1644 m.label_y = m.title;
1649 gradient_mat(grad, mat);
1652 gradient_mat(grad, mat2);
1655 if (erainbow != erNo)
1657 rainbow_mat(erainbow == erBlue, mat);
1660 rainbow_mat(erainbow == erBlue, mat2);
1664 if (mat2.empty() && (elegend != elNone))
1669 if (ecombine && ecombine != ecHalves)
1671 write_combined_matrix(ecombine,
1675 opt2parg_bSet("-cmin", NPA, pa) ? &cmin : nullptr,
1676 opt2parg_bSet("-cmax", NPA, pa) ? &cmax : nullptr);
1695 opt2fn_null("-di", NFILE, fnm),
1696 opt2fn_null("-do", NFILE, fnm),
1701 view_all(oenv, NFILE, fnm);