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, 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.
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.
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.
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.
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.
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.
43 #include "gromacs/commandline/pargs.h"
44 #include "gromacs/fileio/matio.h"
45 #include "gromacs/fileio/trxio.h"
46 #include "gromacs/fileio/writeps.h"
47 #include "gromacs/gmxana/gmx_ana.h"
48 #include "gromacs/legacyheaders/macros.h"
49 #include "gromacs/legacyheaders/typedefs.h"
50 #include "gromacs/legacyheaders/viewit.h"
51 #include "gromacs/utility/fatalerror.h"
52 #include "gromacs/utility/futil.h"
53 #include "gromacs/utility/smalloc.h"
70 char tickfont[STRLEN];
85 char leglabel[STRLEN];
86 char leg2label[STRLEN];
96 /* MUST correspond to char *legend[] in main() */
98 elSel, elBoth, elFirst, elSecond, elNone, elNR
101 /* MUST correspond to char *combine[] in main() */
103 ecSel, ecHalves, ecAdd, ecSub, ecMult, ecDiv, ecNR
106 void get_params(const char *mpin, const char *mpout, t_psrec *psr)
108 static const char *gmx_bools[BOOL_NR+1] = { "no", "yes", NULL };
109 /* this must correspond to t_rgb *linecolors[] below */
110 static const char *colors[] = { "none", "black", "white", NULL };
116 wi = init_warning(FALSE, 0);
120 inp = read_inpfile(mpin, &ninp, wi);
126 ETYPE("black&white", psr->bw, gmx_bools);
127 RTYPE("linewidth", psr->linewidth, 1.0);
128 STYPE("titlefont", psr->titfont, "Helvetica");
129 RTYPE("titlefontsize", psr->titfontsize, 20.0);
130 ETYPE("legend", psr->legend, gmx_bools);
131 STYPE("legendfont", psr->legfont, psr->titfont);
132 STYPE("legendlabel", psr->leglabel, "");
133 STYPE("legend2label", psr->leg2label, psr->leglabel);
134 RTYPE("legendfontsize", psr->legfontsize, 14.0);
135 RTYPE("xbox", psr->xboxsize, 0.0);
136 RTYPE("ybox", psr->yboxsize, 0.0);
137 RTYPE("matrixspacing", psr->boxspacing, 20.0);
138 RTYPE("xoffset", psr->xoffs, 0.0);
139 RTYPE("yoffset", psr->yoffs, psr->xoffs);
140 RTYPE("boxlinewidth", psr->boxlinewidth, psr->linewidth);
141 RTYPE("ticklinewidth", psr->ticklinewidth, psr->linewidth);
142 RTYPE("zerolinewidth", psr->zerolinewidth, psr->ticklinewidth);
143 ETYPE("x-lineat0value", psr->X.lineatzero, colors);
144 RTYPE("x-major", psr->X.major, NOTSET);
145 RTYPE("x-minor", psr->X.minor, NOTSET);
146 RTYPE("x-firstmajor", psr->X.offset, 0.0);
147 ETYPE("x-majorat0", psr->X.first, gmx_bools);
148 RTYPE("x-majorticklen", psr->X.majorticklen, 8.0);
149 RTYPE("x-minorticklen", psr->X.minorticklen, 4.0);
150 STYPE("x-label", psr->X.label, "");
151 RTYPE("x-fontsize", psr->X.fontsize, 16.0);
152 STYPE("x-font", psr->X.font, psr->titfont);
153 RTYPE("x-tickfontsize", psr->X.tickfontsize, 10.0);
154 STYPE("x-tickfont", psr->X.tickfont, psr->X.font);
155 ETYPE("y-lineat0value", psr->Y.lineatzero, colors);
156 RTYPE("y-major", psr->Y.major, psr->X.major);
157 RTYPE("y-minor", psr->Y.minor, psr->X.minor);
158 RTYPE("y-firstmajor", psr->Y.offset, psr->X.offset);
159 ETYPE("y-majorat0", psr->Y.first, gmx_bools);
160 RTYPE("y-majorticklen", psr->Y.majorticklen, psr->X.majorticklen);
161 RTYPE("y-minorticklen", psr->Y.minorticklen, psr->X.minorticklen);
162 STYPE("y-label", psr->Y.label, psr->X.label);
163 RTYPE("y-fontsize", psr->Y.fontsize, psr->X.fontsize);
164 STYPE("y-font", psr->Y.font, psr->X.font);
165 RTYPE("y-tickfontsize", psr->Y.tickfontsize, psr->X.tickfontsize);
166 STYPE("y-tickfont", psr->Y.tickfont, psr->Y.font);
170 write_inpfile(mpout, ninp, inp, TRUE, wi);
173 done_warning(wi, FARGS);
176 t_rgb black = { 0, 0, 0 };
177 t_rgb white = { 1, 1, 1 };
178 t_rgb red = { 1, 0, 0 };
179 t_rgb blue = { 0, 0, 1 };
180 #define BLACK (&black)
181 /* this must correspond to *colors[] in get_params */
182 t_rgb *linecolors[] = { NULL, &black, &white, NULL };
184 gmx_bool diff_maps(int nmap1, t_mapping *map1, int nmap2, t_mapping *map2)
187 gmx_bool bDiff, bColDiff = FALSE;
196 for (i = 0; i < nmap1; i++)
198 if (!matelmt_cmp(map1[i].code, map2[i].code))
202 if (strcmp(map1[i].desc, map2[i].desc) != 0)
206 if ((map1[i].rgb.r != map2[i].rgb.r) ||
207 (map1[i].rgb.g != map2[i].rgb.g) ||
208 (map1[i].rgb.b != map2[i].rgb.b))
213 if (!bDiff && bColDiff)
215 fprintf(stderr, "Warning: two colormaps differ only in RGB value, using one colormap.\n");
222 void leg_discrete(t_psdata ps, real x0, real y0, char *label,
223 real fontsize, char *font, int nmap, t_mapping map[])
229 boxhh = fontsize+DDD;
232 ps_strfont(ps, font, fontsize);
233 yhh = y0+fontsize+3*DDD;
234 if ((int)strlen(label) > 0)
236 ps_ctext(ps, x0, yhh, label, eXLeft);
238 ps_moveto(ps, x0, y0);
239 for (i = 0; (i < nmap); i++)
242 ps_rgb(ps, &(map[i].rgb));
243 ps_fillbox(ps, DDD, DDD, DDD+fontsize, boxhh-DDD);
245 ps_box(ps, DDD, DDD, DDD+fontsize, boxhh-DDD);
246 ps_ctext(ps, boxhh+2*DDD, fontsize/3, map[i].desc, eXLeft);
248 ps_moverel(ps, DDD, -fontsize/3);
252 void leg_continuous(t_psdata ps, real x0, real x, real y0, char *label,
253 real fontsize, char *font,
254 int nmap, t_mapping map[],
259 real yhh, boxxh, boxyh;
266 boxxh = (real)x/(real)(nmap-mapoffset);
267 if (boxxh > fontsize)
273 xx0 = x0-((nmap-mapoffset)*boxxh)/2.0;
275 for (i = 0; (i < nmap-mapoffset); i++)
277 ps_rgb(ps, &(map[i+mapoffset].rgb));
278 ps_fillbox(ps, xx0+i*boxxh, y0, xx0+(i+1)*boxxh, y0+boxyh);
280 ps_strfont(ps, font, fontsize);
282 ps_box(ps, xx0, y0, xx0+(nmap-mapoffset)*boxxh, y0+boxyh);
284 yhh = y0+boxyh+3*DDD;
285 ps_ctext(ps, xx0+boxxh/2, yhh, map[0].desc, eXCenter);
286 if ((int)strlen(label) > 0)
288 ps_ctext(ps, x0, yhh, label, eXCenter);
290 ps_ctext(ps, xx0+((nmap-mapoffset)*boxxh)
291 - boxxh/2, yhh, map[nmap-1].desc, eXCenter);
294 void leg_bicontinuous(t_psdata ps, real x0, real x, real y0, char *label1,
295 char *label2, real fontsize, char *font,
296 int nmap1, t_mapping map1[], int nmap2, t_mapping map2[])
298 real xx1, xx2, x1, x2;
300 x1 = x/(nmap1+nmap2)*nmap1; /* width of legend 1 */
301 x2 = x/(nmap1+nmap2)*nmap2; /* width of legend 2 */
302 xx1 = x0-(x2/2.0)-fontsize; /* center of legend 1 */
303 xx2 = x0+(x1/2.0)+fontsize; /* center of legend 2 */
304 x1 -= fontsize/2; /* adjust width */
305 x2 -= fontsize/2; /* adjust width */
306 leg_continuous(ps, xx1, x1, y0, label1, fontsize, font, nmap1, map1, 0);
307 leg_continuous(ps, xx2, x2, y0, label2, fontsize, font, nmap2, map2, 0);
310 static real box_height(t_matrix *mat, t_psrec *psr)
312 return mat->ny*psr->yboxsize;
315 static real box_dh(t_psrec *psr)
317 return psr->boxspacing;
320 #define IS_ONCE (i == nmat-1)
321 static real box_dh_top(gmx_bool bOnce, t_psrec *psr)
325 if (psr->bTitle || (psr->bTitleOnce && bOnce) )
327 dh = 2*psr->titfontsize;
337 static gmx_bool box_do_all_x_maj_ticks(t_psrec *psr)
339 return (psr->boxspacing > (1.5*psr->X.majorticklen));
342 static gmx_bool box_do_all_x_min_ticks(t_psrec *psr)
344 return (psr->boxspacing > (1.5*psr->X.minorticklen));
347 static void draw_boxes(t_psdata ps, real x0, real y0, real w,
348 int nmat, t_matrix mat[], t_psrec *psr)
353 char **xtick, **ytick;
354 real xx, yy, dy, xx00, yy00, offset_x, offset_y;
355 int i, j, x, y, ntx, nty, strlength;
357 /* Only necessary when there will be no y-labels */
362 ps_linewidth(ps, psr->boxlinewidth);
364 for (i = 0; (i < nmat); i++)
366 dy = box_height(&(mat[i]), psr);
367 ps_box(ps, x0-1, yy00-1, x0+w+1, yy00+dy+1);
368 yy00 += dy+box_dh(psr)+box_dh_top(IS_ONCE, psr);
371 /* Draw the ticks on the axes */
372 ps_linewidth(ps, psr->ticklinewidth);
375 for (i = 0; (i < nmat); i++)
377 if (mat[i].flags & MAT_SPATIAL_X)
387 if (mat[i].flags & MAT_SPATIAL_Y)
398 for (j = 0; (j < ntx); j++)
400 sprintf(buf, "%g", mat[i].axis_x[j]);
401 xtick[j] = gmx_strdup(buf);
403 ps_strfont(ps, psr->X.tickfont, psr->X.tickfontsize);
404 for (x = 0; (x < ntx); x++)
406 xx = xx00 + (x + offset_x)*psr->xboxsize;
407 if ( ( bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.major) ||
408 (psr->X.first && (x == 0))) &&
409 ( (i == 0) || box_do_all_x_maj_ticks(psr) ) )
411 /* Longer tick marks */
412 ps_line (ps, xx, yy00, xx, yy00-psr->X.majorticklen);
413 /* Plot label on lowest graph only */
417 yy00-DDD-psr->X.majorticklen-psr->X.tickfontsize*0.8,
421 else if (bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.minor) &&
422 ( (i == 0) || box_do_all_x_min_ticks(psr) ) )
424 /* Shorter tick marks */
425 ps_line(ps, xx, yy00, xx, yy00-psr->X.minorticklen);
427 else if (bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.major) )
429 /* Even shorter marks, only each X.major */
430 ps_line(ps, xx, yy00, xx, yy00-(psr->boxspacing/2));
433 ps_strfont(ps, psr->Y.tickfont, psr->Y.tickfontsize);
435 for (j = 0; (j < nty); j++)
437 sprintf(buf, "%g", mat[i].axis_y[j]);
438 ytick[j] = gmx_strdup(buf);
441 for (y = 0; (y < nty); y++)
443 yy = yy00 + (y + offset_y)*psr->yboxsize;
444 if (bRmod(mat[i].axis_y[y], psr->Y.offset, psr->Y.major) ||
445 (psr->Y.first && (y == 0)))
448 strlength = max(strlength, (int)strlen(ytick[y]));
449 ps_line (ps, xx00, yy, xx00-psr->Y.majorticklen, yy);
450 ps_ctext(ps, xx00-psr->Y.majorticklen-DDD,
451 yy-psr->Y.tickfontsize/3.0, ytick[y], eXRight);
453 else if (bRmod(mat[i].axis_y[y], psr->Y.offset, psr->Y.minor) )
456 ps_line(ps, xx00, yy, xx00-psr->Y.minorticklen, yy);
462 /* Label on Y-axis */
463 if (!psr->bYonce || i == nmat/2)
465 if (strlen(psr->Y.label) > 0)
467 mylab = psr->Y.label;
471 mylab = mat[i].label_y;
473 if (strlen(mylab) > 0)
475 ps_strfont(ps, psr->Y.font, psr->Y.fontsize);
477 xxx = x0-psr->X.majorticklen-psr->X.tickfontsize*strlength-DDD;
478 ps_ctext(ps, yy00+box_height(&mat[i], psr)/2.0, 612.5-xxx,
484 yy00 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
486 /* Label on X-axis */
487 if (strlen(psr->X.label) > 0)
489 mylab = psr->X.label;
493 mylab = mat[0].label_x;
495 if (strlen(mylab) > 0)
497 ps_strfont(ps, psr->X.font, psr->X.fontsize);
498 ps_ctext(ps, x0+w/2, y0-DDD-psr->X.majorticklen-psr->X.tickfontsize*FUDGE-
499 psr->X.fontsize, mylab, eXCenter);
503 static void draw_zerolines(t_psdata out, real x0, real y0, real w,
504 int nmat, t_matrix mat[], t_psrec *psr)
506 real xx, yy, dy, xx00, yy00;
511 ps_linewidth(out, psr->zerolinewidth);
512 for (i = 0; (i < nmat); i++)
514 dy = box_height(&(mat[i]), psr);
515 /* mat[i].axis_x and _y were already set by draw_boxes */
516 if (psr->X.lineatzero)
518 ps_rgb(out, linecolors[psr->X.lineatzero]);
519 for (x = 0; (x < mat[i].nx); x++)
521 xx = xx00+(x+0.7)*psr->xboxsize;
522 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
523 if (x != 0 && x < mat[i].nx-1 &&
524 fabs(mat[i].axis_x[x]) <
525 0.1*fabs(mat[i].axis_x[x+1]-mat[i].axis_x[x]) )
527 ps_line (out, xx, yy00, xx, yy00+dy+2);
531 if (psr->Y.lineatzero)
533 ps_rgb(out, linecolors[psr->Y.lineatzero]);
534 for (y = 0; (y < mat[i].ny); y++)
536 yy = yy00+(y+0.7)*psr->yboxsize;
537 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
538 if (y != 0 && y < mat[i].ny-1 &&
539 fabs(mat[i].axis_y[y]) <
540 0.1*fabs(mat[i].axis_y[y+1]-mat[i].axis_y[y]) )
542 ps_line (out, xx00, yy, xx00+w+2, yy);
546 yy00 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
550 static void box_dim(int nmat, t_matrix mat[], t_matrix *mat2, t_psrec *psr,
551 int elegend, gmx_bool bFrame,
552 real *w, real *h, real *dw, real *dh)
555 real ww, hh, dww, dhh;
561 for (i = 0; (i < nmat); i++)
563 ww = max(ww, mat[i].nx*psr->xboxsize);
564 hh += box_height(&(mat[i]), psr);
565 maxytick = max(maxytick, mat[i].nx);
569 if (mat[0].label_y[0])
571 dww += 2.0*(psr->Y.fontsize+DDD);
573 if (psr->Y.major > 0)
575 dww += psr->Y.majorticklen + DDD +
576 psr->Y.tickfontsize*(log(maxytick)/log(10.0));
578 else if (psr->Y.minor > 0)
580 dww += psr->Y.minorticklen;
583 if (mat[0].label_x[0])
585 dhh += psr->X.fontsize+2*DDD;
587 if ( /* fool emacs auto-indent */
588 (elegend == elBoth && (mat[0].legend[0] || mat2[0].legend[0])) ||
589 (elegend == elFirst && mat[0].legend[0]) ||
590 (elegend == elSecond && mat2[0].legend[0]) )
592 dhh += 2*(psr->legfontsize*FUDGE+2*DDD);
596 dhh += psr->legfontsize*FUDGE+2*DDD;
598 if (psr->X.major > 0)
600 dhh += psr->X.tickfontsize*FUDGE+2*DDD+psr->X.majorticklen;
602 else if (psr->X.minor > 0)
604 dhh += psr->X.minorticklen;
607 hh += (nmat-1)*box_dh(psr);
608 hh += box_dh_top(TRUE, psr);
611 hh += (nmat-1)*box_dh_top(FALSE, psr);
620 int add_maps(t_mapping **newmap,
621 int nmap1, t_mapping map1[], int nmap2, t_mapping map2[])
623 static char mapper[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<.>/?";
628 nsymbols = strlen(mapper);
630 if (nmap > nsymbols*nsymbols)
632 gmx_fatal(FARGS, "Not enough symbols to merge the two colormaps\n");
634 printf("Combining colormaps of %d and %d elements into one of %d elements\n",
637 for (j = 0; j < nmap1; j++)
639 map[j].code.c1 = mapper[j % nsymbols];
642 map[j].code.c2 = mapper[j/nsymbols];
644 map[j].rgb.r = map1[j].rgb.r;
645 map[j].rgb.g = map1[j].rgb.g;
646 map[j].rgb.b = map1[j].rgb.b;
647 map[j].desc = map1[j].desc;
649 for (j = 0; j < nmap2; j++)
652 map[k].code.c1 = mapper[k % nsymbols];
655 map[k].code.c2 = mapper[k/nsymbols];
657 map[k].rgb.r = map2[j].rgb.r;
658 map[k].rgb.g = map2[j].rgb.g;
659 map[k].rgb.b = map2[j].rgb.b;
660 map[k].desc = map2[j].desc;
667 void xpm_mat(const char *outf, int nmat, t_matrix *mat, t_matrix *mat2,
668 gmx_bool bDiag, gmx_bool bFirstDiag)
672 int i, j, k, x, y, col;
674 t_mapping *map = NULL;
676 out = gmx_ffopen(outf, "w");
678 for (i = 0; i < nmat; i++)
680 if (!mat2 || !diff_maps(mat[i].nmap, mat[i].map, mat2[i].nmap, mat2[i].map))
682 write_xpm_m(out, mat[0]);
686 nmap = add_maps(&map, mat[i].nmap, mat[i].map, mat2[i].nmap, mat2[i].map);
687 for (x = 0; (x < mat[i].nx); x++)
689 for (y = 0; (y < mat[i].nx); y++)
691 if ((x < y) || ((x == y) && bFirstDiag)) /* upper left -> map1 */
693 col = mat[i].matrix[x][y];
695 else /* lower right -> map2 */
697 col = mat[i].nmap+mat[i].matrix[x][y];
699 if ((bDiag) || (x != y))
701 mat[i].matrix[x][y] = col;
705 mat[i].matrix[x][y] = 0;
712 if (strcmp(mat[i].title, mat2[i].title) != 0)
714 sprintf(mat[i].title+strlen(mat[i].title), " / %s", mat2[i].title);
716 if (strcmp(mat[i].legend, mat2[i].legend) != 0)
718 sprintf(mat[i].legend+strlen(mat[i].legend), " / %s", mat2[i].legend);
720 write_xpm_m(out, mat[i]);
726 static void tick_spacing(int n, real axis[], real offset, char axisnm,
727 real *major, real *minor)
730 gmx_bool bTryAgain, bFive;
731 int i, j, t, f = 0, ten;
733 real major_fact[NFACT] = {5, 4, 2, 1};
734 real minor_fact[NFACT] = {5, 4, 4, 5};
736 /* start with interval between 10 matrix points: */
737 space = max(10*axis[1]-axis[0], axis[min(10, n-1)]-axis[0]);
738 /* get power of 10 */
739 ten = (int)ceil(log(space)/log(10))-1;
741 for (t = ten+2; t > ten-3 && bTryAgain; t--)
743 for (f = 0; f < NFACT && bTryAgain; f++)
745 space = pow(10, t) * major_fact[f];
746 /* count how many ticks we would get: */
748 for (j = 0; j < n; j++)
750 if (bRmod(axis[j], offset, space) )
755 /* do we have a reasonable number of ticks ? */
756 bTryAgain = (i > min(10, n-1)) || (i < 5);
761 space = max(10*axis[1]-axis[0], axis[min(10, n-1)]-axis[0]);
762 fprintf(stderr, "Auto tick spacing failed for %c-axis, guessing %g\n",
766 *minor = space / minor_fact[f-1];
767 fprintf(stderr, "Auto tick spacing for %c-axis: major %g, minor %g\n",
768 axisnm, *major, *minor);
771 void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
772 gmx_bool bFrame, gmx_bool bDiag, gmx_bool bFirstDiag,
773 gmx_bool bTitle, gmx_bool bTitleOnce, gmx_bool bYonce, int elegend,
774 real size, real boxx, real boxy, const char *m2p, const char *m2pout,
778 char buf[256], *legend;
782 int i, j, x, y, col, leg = 0;
785 int nmap1 = 0, nmap2 = 0, leg_nmap;
786 t_mapping *map1 = NULL, *map2 = NULL, *leg_map;
787 gmx_bool bMap1, bNextMap1, bDiscrete;
790 libm2p = m2p ? gmxlibfn(m2p) : m2p;
791 get_params(libm2p, m2pout, &psrec);
795 if (psr->X.major <= 0)
797 tick_spacing((mat[0].flags & MAT_SPATIAL_X) ? mat[0].nx + 1 : mat[0].nx,
798 mat[0].axis_x, psr->X.offset, 'X',
799 &(psr->X.major), &(psr->X.minor) );
801 if (psr->X.minor <= 0)
803 psr->X.minor = psr->X.major / 2;
805 if (psr->Y.major <= 0)
807 tick_spacing((mat[0].flags & MAT_SPATIAL_Y) ? mat[0].ny + 1 : mat[0].ny,
808 mat[0].axis_y, psr->Y.offset, 'Y',
809 &(psr->Y.major), &(psr->Y.minor) );
811 if (psr->Y.minor <= 0)
813 psr->Y.minor = psr->Y.major / 2;
818 psr->xboxsize = boxx;
819 psr->yboxsize = boxx;
823 psr->yboxsize = boxy;
826 if (psr->xboxsize == 0)
828 psr->xboxsize = size/mat[0].nx;
829 printf("Set the x-size of the box to %.3f\n", psr->xboxsize);
831 if (psr->yboxsize == 0)
833 psr->yboxsize = size/mat[0].nx;
834 printf("Set the y-size of the box to %.3f\n", psr->yboxsize);
838 for (i = 0; (i < nmat); i++)
840 if (mat[i].nmap > nmap1)
849 printf("Selected legend of matrix # %d for display\n", leg);
854 for (i = 0; (i < nmat); i++)
856 if (mat2[i].nmap > nmap2)
858 nmap2 = mat2[i].nmap;
865 printf("Selected legend of matrix # %d for second display\n", leg);
868 if ( (mat[0].legend[0] == 0) && psr->legend)
870 strcpy(mat[0].legend, psr->leglabel);
873 bTitle = bTitle && mat[nmat-1].title[0];
874 bTitleOnce = bTitleOnce && mat[nmat-1].title[0];
875 psr->bTitle = bTitle;
876 psr->bTitleOnce = bTitleOnce;
877 psr->bYonce = bYonce;
879 /* Set up size of box for nice colors */
880 box_dim(nmat, mat, mat2, psr, elegend, bFrame, &w, &h, &dw, &dh);
882 /* Set up bounding box */
896 out = ps_open(outf, 0, 0, x, y);
897 ps_linewidth(out, psr->linewidth);
898 ps_init_rgb_box(out, psr->xboxsize, psr->yboxsize);
899 ps_init_rgb_nbox(out, psr->xboxsize, psr->yboxsize);
900 ps_translate(out, psr->xoffs, psr->yoffs);
904 ps_comment(out, "Here starts the BOX drawing");
905 draw_boxes(out, x0, y0, w, nmat, mat, psr);
908 for (i = 0; (i < nmat); i++)
910 if (bTitle || (bTitleOnce && i == nmat-1) )
912 /* Print title, if any */
914 ps_strfont(out, psr->titfont, psr->titfontsize);
915 if (!mat2 || (strcmp(mat[i].title, mat2[i].title) == 0))
917 strcpy(buf, mat[i].title);
921 sprintf(buf, "%s / %s", mat[i].title, mat2[i].title);
923 ps_ctext(out, x0+w/2, y0+box_height(&(mat[i]), psr)+psr->titfontsize,
926 sprintf(buf, "Here starts the filling of box #%d", i);
927 ps_comment(out, buf);
928 for (x = 0; (x < mat[i].nx); x++)
933 xx = x0+x*psr->xboxsize;
934 ps_moveto(out, xx, y0);
936 bMap1 = (!mat2 || (x < y || (x == y && bFirstDiag)));
937 if ((bDiag) || (x != y))
939 col = mat[i].matrix[x][y];
945 for (nexty = 1; (nexty <= mat[i].ny); nexty++)
947 bNextMap1 = (!mat2 || (x < nexty || (x == nexty && bFirstDiag)));
948 /* TRUE: upper left -> map1 */
949 /* FALSE: lower right -> map2 */
950 if ((nexty == mat[i].ny) || (!bDiag && (x == nexty)))
956 nextcol = mat[i].matrix[x][nexty];
958 if ( (nexty == mat[i].ny) || (col != nextcol) || (bMap1 != bNextMap1) )
964 ps_rgb_nbox(out, &(mat[i].map[col].rgb), nexty-y);
968 ps_rgb_nbox(out, &(mat2[i].map[col].rgb), nexty-y);
973 ps_moverel(out, 0, psr->yboxsize);
981 y0 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
984 if (psr->X.lineatzero || psr->Y.lineatzero)
986 /* reset y0 for first box */
988 ps_comment(out, "Here starts the zero lines drawing");
989 draw_zerolines(out, x0, y0, w, nmat, mat, psr);
992 if (elegend != elNone)
994 ps_comment(out, "Now it's legend time!");
995 ps_linewidth(out, psr->linewidth);
996 if (mat2 == NULL || elegend != elSecond)
998 bDiscrete = mat[0].bDiscrete;
999 legend = mat[0].legend;
1005 bDiscrete = mat2[0].bDiscrete;
1006 legend = mat2[0].legend;
1012 leg_discrete(out, psr->legfontsize, DDD, legend,
1013 psr->legfontsize, psr->legfont, leg_nmap, leg_map);
1017 if (elegend != elBoth)
1019 leg_continuous(out, x0+w/2, w/2, DDD, legend,
1020 psr->legfontsize, psr->legfont, leg_nmap, leg_map,
1025 leg_bicontinuous(out, x0+w/2, w, DDD, mat[0].legend, mat2[0].legend,
1026 psr->legfontsize, psr->legfont, nmap1, map1, nmap2, map2);
1029 ps_comment(out, "Were there, dude");
1035 void make_axis_labels(int nmat, t_matrix *mat)
1039 for (i = 0; (i < nmat); i++)
1041 /* Make labels for x axis */
1042 if (mat[i].axis_x == NULL)
1044 snew(mat[i].axis_x, mat[i].nx);
1045 for (j = 0; (j < mat[i].nx); j++)
1047 mat[i].axis_x[j] = j;
1050 /* Make labels for y axis */
1051 if (mat[i].axis_y == NULL)
1053 snew(mat[i].axis_y, mat[i].ny);
1054 for (j = 0; (j < mat[i].ny); j++)
1056 mat[i].axis_y[j] = j;
1062 void prune_mat(int nmat, t_matrix *mat, t_matrix *mat2, int skip)
1064 int i, x, y, xs, ys;
1066 for (i = 0; i < nmat; i++)
1068 fprintf(stderr, "converting %dx%d matrix to %dx%d\n",
1069 mat[i].nx, mat[i].ny,
1070 (mat[i].nx+skip-1)/skip, (mat[i].ny+skip-1)/skip);
1071 /* walk through matrix */
1073 for (x = 0; (x < mat[i].nx); x++)
1077 mat[i].axis_x[xs] = mat[i].axis_x[x];
1080 mat2[i].axis_x[xs] = mat2[i].axis_x[x];
1083 for (y = 0; (y < mat[i].ny); y++)
1087 mat[i].axis_y[ys] = mat[i].axis_y[y];
1090 mat2[i].axis_y[ys] = mat2[i].axis_y[y];
1095 mat[i].matrix[xs][ys] = mat[i].matrix[x][y];
1098 mat2[i].matrix[xs][ys] = mat2[i].matrix[x][y];
1106 /* adjust parameters */
1107 mat[i].nx = (mat[i].nx+skip-1)/skip;
1108 mat[i].ny = (mat[i].ny+skip-1)/skip;
1111 mat2[i].nx = (mat2[i].nx+skip-1)/skip;
1112 mat2[i].ny = (mat2[i].ny+skip-1)/skip;
1117 void zero_lines(int nmat, t_matrix *mat, t_matrix *mat2)
1122 for (i = 0; i < nmat; i++)
1124 for (m = 0; m < (mat2 ? 2 : 1); m++)
1134 for (x = 0; x < mats[i].nx-1; x++)
1136 if (fabs(mats[i].axis_x[x+1]) < 1e-5)
1138 for (y = 0; y < mats[i].ny; y++)
1140 mats[i].matrix[x][y] = 0;
1144 for (y = 0; y < mats[i].ny-1; y++)
1146 if (fabs(mats[i].axis_y[y+1]) < 1e-5)
1148 for (x = 0; x < mats[i].nx; x++)
1150 mats[i].matrix[x][y] = 0;
1158 void write_combined_matrix(int ecombine, const char *fn,
1159 int nmat, t_matrix *mat1, t_matrix *mat2,
1160 real *cmin, real *cmax)
1162 int i, j, k, nlevels;
1163 t_mapping *map = NULL;
1165 real **rmat1, **rmat2;
1168 out = gmx_ffopen(fn, "w");
1169 for (k = 0; k < nmat; k++)
1171 if (mat2[k].nx != mat1[k].nx || mat2[k].ny != mat1[k].ny)
1173 gmx_fatal(FARGS, "Size of frame %d in 1st (%dx%d) and 2nd matrix (%dx%d) do"
1174 " not match.\n", k, mat1[k].nx, mat1[k].ny, mat2[k].nx, mat2[k].ny);
1176 printf("Combining two %dx%d matrices\n", mat1[k].nx, mat1[k].ny);
1177 rmat1 = matrix2real(&mat1[k], NULL);
1178 rmat2 = matrix2real(&mat2[k], NULL);
1179 if (NULL == rmat1 || NULL == rmat2)
1181 gmx_fatal(FARGS, "Could not extract real data from %s xpm matrices. Note that, e.g.,\n"
1182 "g_rms and g_mdmat provide such data, but not do_dssp.\n",
1183 (NULL == rmat1 && NULL == rmat2) ? "both" : "one of the" );
1187 for (j = 0; j < mat1[k].ny; j++)
1189 for (i = 0; i < mat1[k].nx; i++)
1193 case ecAdd: rmat1[i][j] += rmat2[i][j]; break;
1194 case ecSub: rmat1[i][j] -= rmat2[i][j]; break;
1195 case ecMult: rmat1[i][j] *= rmat2[i][j]; break;
1196 case ecDiv: rmat1[i][j] /= rmat2[i][j]; break;
1198 gmx_fatal(FARGS, "No such combination rule %d for matrices", ecombine);
1200 rlo = min(rlo, rmat1[i][j]);
1201 rhi = max(rhi, rmat1[i][j]);
1212 nlevels = max(mat1[k].nmap, mat2[k].nmap);
1216 "combination results in uniform matrix (%g), no output\n", rhi);
1219 else if (rlo>=0 || rhi<=0)
1220 write_xpm(out, mat1[k].flags, mat1[k].title, mat1[k].legend,
1221 mat1[k].label_x, mat1[k].label_y,
1222 mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
1223 rmat1, rlo, rhi, rhi<=0?red:white, rhi<=0?white:blue,
1226 write_xpm3(out, mat2[k].flags, mat1[k].title, mat1[k].legend,
1227 mat1[k].label_x, mat1[k].label_y,
1228 mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
1229 rmat1, rlo, 0, rhi, red, white, blue, &nlevels);
1233 write_xpm(out, mat1[k].flags, mat1[k].title, mat1[k].legend,
1234 mat1[k].label_x, mat1[k].label_y,
1235 mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
1236 rmat1, rlo, rhi, white, black, &nlevels);
1242 void do_mat(int nmat, t_matrix *mat, t_matrix *mat2,
1243 gmx_bool bFrame, gmx_bool bZeroLine, gmx_bool bDiag, gmx_bool bFirstDiag, gmx_bool bTitle,
1244 gmx_bool bTitleOnce, gmx_bool bYonce, int elegend,
1245 real size, real boxx, real boxy,
1246 const char *epsfile, const char *xpmfile, const char *m2p,
1247 const char *m2pout, int skip, int mapoffset)
1253 for (k = 0; (k < nmat); k++)
1255 if ((mat2[k].nx != mat[k].nx) || (mat2[k].ny != mat[k].ny))
1257 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",
1258 k, mat2[k].nx, mat2[k].ny, mat[k].nx, mat[k].ny);
1260 for (j = 0; (j < mat[k].ny); j++)
1262 for (i = bFirstDiag ? j+1 : j; (i < mat[k].nx); i++)
1264 mat[k].matrix[i][j] = mat2[k].matrix[i][j];
1269 for (i = 0; (i < nmat); i++)
1271 fprintf(stderr, "Matrix %d is %d x %d\n", i, mat[i].nx, mat[i].ny);
1274 make_axis_labels(nmat, mat);
1278 prune_mat(nmat, mat, mat2, skip);
1283 zero_lines(nmat, mat, mat);
1286 if (epsfile != NULL)
1288 ps_mat(epsfile, nmat, mat, mat2, bFrame, bDiag, bFirstDiag,
1289 bTitle, bTitleOnce, bYonce, elegend,
1290 size, boxx, boxy, m2p, m2pout, mapoffset);
1292 if (xpmfile != NULL)
1294 xpm_mat(xpmfile, nmat, mat, mat2, bDiag, bFirstDiag);
1298 void gradient_map(rvec grad, int nmap, t_mapping map[])
1303 for (i = 0; i < nmap; i++)
1306 map[i].rgb.r = 1-c*(1-grad[XX]);
1307 map[i].rgb.g = 1-c*(1-grad[YY]);
1308 map[i].rgb.b = 1-c*(1-grad[ZZ]);
1312 void gradient_mat(rvec grad, int nmat, t_matrix mat[])
1316 for (m = 0; m < nmat; m++)
1318 gradient_map(grad, mat[m].nmap, mat[m].map);
1322 void rainbow_map(gmx_bool bBlue, int nmap, t_mapping map[])
1327 for (i = 0; i < nmap; i++)
1329 c = (map[i].rgb.r + map[i].rgb.g + map[i].rgb.b)/3;
1338 if (c <= 0.25) /* 0-0.25 */
1341 g = pow(4*c, 0.666);
1344 else if (c <= 0.5) /* 0.25-0.5 */
1348 b = pow(2-4*c, 0.666);
1350 else if (c <= 0.75) /* 0.5-0.75 */
1352 r = pow(4*c-2, 0.666);
1359 g = pow(4-4*c, 0.666);
1368 void rainbow_mat(gmx_bool bBlue, int nmat, t_matrix mat[])
1372 for (m = 0; m < nmat; m++)
1374 rainbow_map(bBlue, mat[m].nmap, mat[m].map);
1378 int gmx_xpm2ps(int argc, char *argv[])
1380 const char *desc[] = {
1381 "[THISMODULE] makes a beautiful color plot of an XPixelMap file.",
1382 "Labels and axis can be displayed, when they are supplied",
1383 "in the correct matrix format.",
1384 "Matrix data may be generated by programs such as [gmx-do_dssp], [gmx-rms] or",
1385 "[gmx-mdmat].[PAR]",
1386 "Parameters are set in the [TT].m2p[tt] file optionally supplied with",
1387 "[TT]-di[tt]. Reasonable defaults are provided. Settings for the [IT]y[it]-axis",
1388 "default to those for the [IT]x[it]-axis. Font names have a defaulting hierarchy:",
1389 "titlefont -> legendfont; titlefont -> (xfont -> yfont -> ytickfont)",
1390 "-> xtickfont, e.g. setting titlefont sets all fonts, setting xfont",
1391 "sets yfont, ytickfont and xtickfont.[PAR]",
1392 "When no [TT].m2p[tt] file is supplied, many settings are taken from",
1393 "command line options. The most important option is [TT]-size[tt],",
1394 "which sets the size of the whole matrix in postscript units.",
1395 "This option can be overridden with the [TT]-bx[tt] and [TT]-by[tt]",
1396 "options (and the corresponding parameters in the [TT].m2p[tt] file),",
1397 "which set the size of a single matrix element.[PAR]",
1398 "With [TT]-f2[tt] a second matrix file can be supplied. Both matrix",
1399 "files will be read simultaneously and the upper left half of the",
1400 "first one ([TT]-f[tt]) is plotted together with the lower right",
1401 "half of the second one ([TT]-f2[tt]). The diagonal will contain",
1402 "values from the matrix file selected with [TT]-diag[tt].",
1403 "Plotting of the diagonal values can be suppressed altogether by",
1404 "setting [TT]-diag[tt] to [TT]none[tt].",
1405 "In this case, a new color map will be generated with",
1406 "a red gradient for negative numbers and a blue for positive.",
1407 "If the color coding and legend labels of both matrices are identical,",
1408 "only one legend will be displayed, else two separate legends are",
1410 "With [TT]-combine[tt], an alternative operation can be selected",
1411 "to combine the matrices. The output range is automatically set",
1412 "to the actual range of the combined matrix. This can be overridden",
1413 "with [TT]-cmin[tt] and [TT]-cmax[tt].[PAR]",
1414 "[TT]-title[tt] can be set to [TT]none[tt] to suppress the title, or to",
1415 "[TT]ylabel[tt] to show the title in the Y-label position (alongside",
1416 "the [IT]y[it]-axis).[PAR]",
1417 "With the [TT]-rainbow[tt] option, dull grayscale matrices can be turned",
1418 "into attractive color pictures.[PAR]",
1419 "Merged or rainbowed matrices can be written to an XPixelMap file with",
1420 "the [TT]-xpm[tt] option."
1424 const char *fn, *epsfile = NULL, *xpmfile = NULL;
1425 int i, nmat, nmat2, etitle, elegend, ediag, erainbow, ecombine;
1426 t_matrix *mat = NULL, *mat2 = NULL;
1427 gmx_bool bTitle, bTitleOnce, bDiag, bFirstDiag, bGrad;
1428 static gmx_bool bFrame = TRUE, bZeroLine = FALSE, bYonce = FALSE, bAdd = FALSE;
1429 static real size = 400, boxx = 0, boxy = 0, cmin = 0, cmax = 0;
1430 static rvec grad = {0, 0, 0};
1432 etSel, etTop, etOnce, etYlabel, etNone, etNR
1434 const char *title[] = { NULL, "top", "once", "ylabel", "none", NULL };
1435 /* MUST correspond to enum elXxx as defined at top of file */
1436 const char *legend[] = { NULL, "both", "first", "second", "none", NULL };
1438 edSel, edFirst, edSecond, edNone, edNR
1440 const char *diag[] = { NULL, "first", "second", "none", NULL };
1442 erSel, erNo, erBlue, erRed, erNR
1444 const char *rainbow[] = { NULL, "no", "blue", "red", NULL };
1445 /* MUST correspond to enum ecXxx as defined at top of file */
1446 const char *combine[] = {
1447 NULL, "halves", "add", "sub", "mult", "div", NULL
1449 static int skip = 1, mapoffset = 0;
1451 { "-frame", FALSE, etBOOL, {&bFrame},
1452 "Display frame, ticks, labels, title and legend" },
1453 { "-title", FALSE, etENUM, {title}, "Show title at" },
1454 { "-yonce", FALSE, etBOOL, {&bYonce}, "Show y-label only once" },
1455 { "-legend", FALSE, etENUM, {legend}, "Show legend" },
1456 { "-diag", FALSE, etENUM, {diag}, "Diagonal" },
1457 { "-size", FALSE, etREAL, {&size},
1458 "Horizontal size of the matrix in ps units" },
1459 { "-bx", FALSE, etREAL, {&boxx},
1460 "Element x-size, overrides [TT]-size[tt] (also y-size when [TT]-by[tt] is not set)" },
1461 { "-by", FALSE, etREAL, {&boxy}, "Element y-size" },
1462 { "-rainbow", FALSE, etENUM, {rainbow},
1463 "Rainbow colors, convert white to" },
1464 { "-gradient", FALSE, etRVEC, {grad},
1465 "Re-scale colormap to a smooth gradient from white {1,1,1} to {r,g,b}" },
1466 { "-skip", FALSE, etINT, {&skip},
1467 "only write out every nr-th row and column" },
1468 { "-zeroline", FALSE, etBOOL, {&bZeroLine},
1469 "insert line in [TT].xpm[tt] matrix where axis label is zero"},
1470 { "-legoffset", FALSE, etINT, {&mapoffset},
1471 "Skip first N colors from [TT].xpm[tt] file for the legend" },
1472 { "-combine", FALSE, etENUM, {combine}, "Combine two matrices" },
1473 { "-cmin", FALSE, etREAL, {&cmin}, "Minimum for combination output" },
1474 { "-cmax", FALSE, etREAL, {&cmax}, "Maximum for combination output" }
1476 #define NPA asize(pa)
1478 { efXPM, "-f", NULL, ffREAD },
1479 { efXPM, "-f2", "root2", ffOPTRD },
1480 { efM2P, "-di", NULL, ffLIBOPTRD },
1481 { efM2P, "-do", "out", ffOPTWR },
1482 { efEPS, "-o", NULL, ffOPTWR },
1483 { efXPM, "-xpm", NULL, ffOPTWR }
1485 #define NFILE asize(fnm)
1487 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW,
1488 NFILE, fnm, NPA, pa,
1489 asize(desc), desc, 0, NULL, &oenv))
1494 etitle = nenum(title);
1495 elegend = nenum(legend);
1496 ediag = nenum(diag);
1497 erainbow = nenum(rainbow);
1498 ecombine = nenum(combine);
1499 bGrad = opt2parg_bSet("-gradient", NPA, pa);
1500 for (i = 0; i < DIM; i++)
1502 if (grad[i] < 0 || grad[i] > 1)
1504 gmx_fatal(FARGS, "RGB value %g out of range (0.0-1.0)", grad[i]);
1513 epsfile = ftp2fn_null(efEPS, NFILE, fnm);
1514 xpmfile = opt2fn_null("-xpm", NFILE, fnm);
1515 if (epsfile == NULL && xpmfile == NULL)
1517 if (ecombine != ecHalves)
1519 xpmfile = opt2fn("-xpm", NFILE, fnm);
1523 epsfile = ftp2fn(efEPS, NFILE, fnm);
1526 if (ecombine != ecHalves && epsfile)
1529 "WARNING: can only write result of arithmetic combination "
1530 "of two matrices to .xpm file\n"
1531 " file %s will not be written\n", epsfile);
1535 bDiag = ediag != edNone;
1536 bFirstDiag = ediag != edSecond;
1538 fn = opt2fn("-f", NFILE, fnm);
1539 nmat = read_xpm_matrix(fn, &mat);
1540 fprintf(stderr, "There %s %d matri%s in %s\n", (nmat > 1) ? "are" : "is", nmat, (nmat > 1) ? "ces" : "x", fn);
1541 fn = opt2fn_null("-f2", NFILE, fnm);
1544 nmat2 = read_xpm_matrix(fn, &mat2);
1545 fprintf(stderr, "There %s %d matri%s in %s\n", (nmat2 > 1) ? "are" : "is", nmat2, (nmat2 > 1) ? "ces" : "x", fn);
1548 fprintf(stderr, "Different number of matrices, using the smallest number.\n");
1549 nmat = nmat2 = min(nmat, nmat2);
1554 if (ecombine != ecHalves)
1557 "WARNING: arithmetic matrix combination selected (-combine), "
1558 "but no second matrix (-f2) supplied\n"
1559 " no matrix combination will be performed\n");
1564 bTitle = etitle == etTop;
1565 bTitleOnce = etitle == etOnce;
1566 if (etitle == etYlabel)
1568 for (i = 0; (i < nmat); i++)
1570 strcpy(mat[i].label_y, mat[i].title);
1573 strcpy(mat2[i].label_y, mat2[i].title);
1579 gradient_mat(grad, nmat, mat);
1582 gradient_mat(grad, nmat2, mat2);
1585 if (erainbow != erNo)
1587 rainbow_mat(erainbow == erBlue, nmat, mat);
1590 rainbow_mat(erainbow == erBlue, nmat2, mat2);
1594 if ((mat2 == NULL) && (elegend != elNone))
1599 if (ecombine && ecombine != ecHalves)
1601 write_combined_matrix(ecombine, xpmfile, nmat, mat, mat2,
1602 opt2parg_bSet("-cmin", NPA, pa) ? &cmin : NULL,
1603 opt2parg_bSet("-cmax", NPA, pa) ? &cmax : NULL);
1607 do_mat(nmat, mat, mat2, bFrame, bZeroLine, bDiag, bFirstDiag,
1608 bTitle, bTitleOnce, bYonce,
1609 elegend, size, boxx, boxy, epsfile, xpmfile,
1610 opt2fn_null("-di", NFILE, fnm), opt2fn_null("-do", NFILE, fnm), skip,
1614 view_all(oenv, NFILE, fnm);