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.
47 #include "gmx_fatal.h"
53 #include "gromacs/commandline/pargs.h"
54 #include "gromacs/fileio/futil.h"
55 #include "gromacs/fileio/matio.h"
56 #include "gromacs/fileio/trxio.h"
57 #include "gromacs/fileio/writeps.h"
74 char tickfont[STRLEN];
89 char leglabel[STRLEN];
90 char leg2label[STRLEN];
100 /* MUST correspond to char *legend[] in main() */
102 elSel, elBoth, elFirst, elSecond, elNone, elNR
105 /* MUST correspond to char *combine[] in main() */
107 ecSel, ecHalves, ecAdd, ecSub, ecMult, ecDiv, ecNR
110 void get_params(const char *mpin, const char *mpout, t_psrec *psr)
112 static const char *gmx_bools[BOOL_NR+1] = { "no", "yes", NULL };
113 /* this must correspond to t_rgb *linecolors[] below */
114 static const char *colors[] = { "none", "black", "white", NULL };
120 wi = init_warning(FALSE, 0);
124 inp = read_inpfile(mpin, &ninp, wi);
130 ETYPE("black&white", psr->bw, gmx_bools);
131 RTYPE("linewidth", psr->linewidth, 1.0);
132 STYPE("titlefont", psr->titfont, "Helvetica");
133 RTYPE("titlefontsize", psr->titfontsize, 20.0);
134 ETYPE("legend", psr->legend, gmx_bools);
135 STYPE("legendfont", psr->legfont, psr->titfont);
136 STYPE("legendlabel", psr->leglabel, "");
137 STYPE("legend2label", psr->leg2label, psr->leglabel);
138 RTYPE("legendfontsize", psr->legfontsize, 14.0);
139 RTYPE("xbox", psr->xboxsize, 0.0);
140 RTYPE("ybox", psr->yboxsize, 0.0);
141 RTYPE("matrixspacing", psr->boxspacing, 20.0);
142 RTYPE("xoffset", psr->xoffs, 0.0);
143 RTYPE("yoffset", psr->yoffs, psr->xoffs);
144 RTYPE("boxlinewidth", psr->boxlinewidth, psr->linewidth);
145 RTYPE("ticklinewidth", psr->ticklinewidth, psr->linewidth);
146 RTYPE("zerolinewidth", psr->zerolinewidth, psr->ticklinewidth);
147 ETYPE("x-lineat0value", psr->X.lineatzero, colors);
148 RTYPE("x-major", psr->X.major, NOTSET);
149 RTYPE("x-minor", psr->X.minor, NOTSET);
150 RTYPE("x-firstmajor", psr->X.offset, 0.0);
151 ETYPE("x-majorat0", psr->X.first, gmx_bools);
152 RTYPE("x-majorticklen", psr->X.majorticklen, 8.0);
153 RTYPE("x-minorticklen", psr->X.minorticklen, 4.0);
154 STYPE("x-label", psr->X.label, "");
155 RTYPE("x-fontsize", psr->X.fontsize, 16.0);
156 STYPE("x-font", psr->X.font, psr->titfont);
157 RTYPE("x-tickfontsize", psr->X.tickfontsize, 10.0);
158 STYPE("x-tickfont", psr->X.tickfont, psr->X.font);
159 ETYPE("y-lineat0value", psr->Y.lineatzero, colors);
160 RTYPE("y-major", psr->Y.major, psr->X.major);
161 RTYPE("y-minor", psr->Y.minor, psr->X.minor);
162 RTYPE("y-firstmajor", psr->Y.offset, psr->X.offset);
163 ETYPE("y-majorat0", psr->Y.first, gmx_bools);
164 RTYPE("y-majorticklen", psr->Y.majorticklen, psr->X.majorticklen);
165 RTYPE("y-minorticklen", psr->Y.minorticklen, psr->X.minorticklen);
166 STYPE("y-label", psr->Y.label, psr->X.label);
167 RTYPE("y-fontsize", psr->Y.fontsize, psr->X.fontsize);
168 STYPE("y-font", psr->Y.font, psr->X.font);
169 RTYPE("y-tickfontsize", psr->Y.tickfontsize, psr->X.tickfontsize);
170 STYPE("y-tickfont", psr->Y.tickfont, psr->Y.font);
174 write_inpfile(mpout, ninp, inp, TRUE, wi);
177 done_warning(wi, FARGS);
180 t_rgb black = { 0, 0, 0 };
181 t_rgb white = { 1, 1, 1 };
182 t_rgb red = { 1, 0, 0 };
183 t_rgb blue = { 0, 0, 1 };
184 #define BLACK (&black)
185 /* this must correspond to *colors[] in get_params */
186 t_rgb *linecolors[] = { NULL, &black, &white, NULL };
188 gmx_bool diff_maps(int nmap1, t_mapping *map1, int nmap2, t_mapping *map2)
191 gmx_bool bDiff, bColDiff = FALSE;
200 for (i = 0; i < nmap1; i++)
202 if (!matelmt_cmp(map1[i].code, map2[i].code))
206 if (strcmp(map1[i].desc, map2[i].desc) != 0)
210 if ((map1[i].rgb.r != map2[i].rgb.r) ||
211 (map1[i].rgb.g != map2[i].rgb.g) ||
212 (map1[i].rgb.b != map2[i].rgb.b))
217 if (!bDiff && bColDiff)
219 fprintf(stderr, "Warning: two colormaps differ only in RGB value, using one colormap.\n");
226 void leg_discrete(t_psdata ps, real x0, real y0, char *label,
227 real fontsize, char *font, int nmap, t_mapping map[])
233 boxhh = fontsize+DDD;
236 ps_strfont(ps, font, fontsize);
237 yhh = y0+fontsize+3*DDD;
238 if ((int)strlen(label) > 0)
240 ps_ctext(ps, x0, yhh, label, eXLeft);
242 ps_moveto(ps, x0, y0);
243 for (i = 0; (i < nmap); i++)
246 ps_rgb(ps, &(map[i].rgb));
247 ps_fillbox(ps, DDD, DDD, DDD+fontsize, boxhh-DDD);
249 ps_box(ps, DDD, DDD, DDD+fontsize, boxhh-DDD);
250 ps_ctext(ps, boxhh+2*DDD, fontsize/3, map[i].desc, eXLeft);
252 ps_moverel(ps, DDD, -fontsize/3);
256 void leg_continuous(t_psdata ps, real x0, real x, real y0, char *label,
257 real fontsize, char *font,
258 int nmap, t_mapping map[],
263 real yhh, boxxh, boxyh;
270 boxxh = (real)x/(real)(nmap-mapoffset);
271 if (boxxh > fontsize)
277 xx0 = x0-((nmap-mapoffset)*boxxh)/2.0;
279 for (i = 0; (i < nmap-mapoffset); i++)
281 ps_rgb(ps, &(map[i+mapoffset].rgb));
282 ps_fillbox(ps, xx0+i*boxxh, y0, xx0+(i+1)*boxxh, y0+boxyh);
284 ps_strfont(ps, font, fontsize);
286 ps_box(ps, xx0, y0, xx0+(nmap-mapoffset)*boxxh, y0+boxyh);
288 yhh = y0+boxyh+3*DDD;
289 ps_ctext(ps, xx0+boxxh/2, yhh, map[0].desc, eXCenter);
290 if ((int)strlen(label) > 0)
292 ps_ctext(ps, x0, yhh, label, eXCenter);
294 ps_ctext(ps, xx0+((nmap-mapoffset)*boxxh)
295 - boxxh/2, yhh, map[nmap-1].desc, eXCenter);
298 void leg_bicontinuous(t_psdata ps, real x0, real x, real y0, char *label1,
299 char *label2, real fontsize, char *font,
300 int nmap1, t_mapping map1[], int nmap2, t_mapping map2[])
302 real xx1, xx2, x1, x2;
304 x1 = x/(nmap1+nmap2)*nmap1; /* width of legend 1 */
305 x2 = x/(nmap1+nmap2)*nmap2; /* width of legend 2 */
306 xx1 = x0-(x2/2.0)-fontsize; /* center of legend 1 */
307 xx2 = x0+(x1/2.0)+fontsize; /* center of legend 2 */
308 x1 -= fontsize/2; /* adjust width */
309 x2 -= fontsize/2; /* adjust width */
310 leg_continuous(ps, xx1, x1, y0, label1, fontsize, font, nmap1, map1, 0);
311 leg_continuous(ps, xx2, x2, y0, label2, fontsize, font, nmap2, map2, 0);
314 static real box_height(t_matrix *mat, t_psrec *psr)
316 return mat->ny*psr->yboxsize;
319 static real box_dh(t_psrec *psr)
321 return psr->boxspacing;
324 #define IS_ONCE (i == nmat-1)
325 static real box_dh_top(gmx_bool bOnce, t_psrec *psr)
329 if (psr->bTitle || (psr->bTitleOnce && bOnce) )
331 dh = 2*psr->titfontsize;
341 static gmx_bool box_do_all_x_maj_ticks(t_psrec *psr)
343 return (psr->boxspacing > (1.5*psr->X.majorticklen));
346 static gmx_bool box_do_all_x_min_ticks(t_psrec *psr)
348 return (psr->boxspacing > (1.5*psr->X.minorticklen));
351 static void draw_boxes(t_psdata ps, real x0, real y0, real w,
352 int nmat, t_matrix mat[], t_psrec *psr)
357 char **xtick, **ytick;
358 real xx, yy, dy, xx00, yy00, offset_x, offset_y;
359 int i, j, x, y, ntx, nty, strlength;
361 /* Only necessary when there will be no y-labels */
366 ps_linewidth(ps, psr->boxlinewidth);
368 for (i = 0; (i < nmat); i++)
370 dy = box_height(&(mat[i]), psr);
371 ps_box(ps, x0-1, yy00-1, x0+w+1, yy00+dy+1);
372 yy00 += dy+box_dh(psr)+box_dh_top(IS_ONCE, psr);
375 /* Draw the ticks on the axes */
376 ps_linewidth(ps, psr->ticklinewidth);
379 for (i = 0; (i < nmat); i++)
381 if (mat[i].flags & MAT_SPATIAL_X)
391 if (mat[i].flags & MAT_SPATIAL_Y)
402 for (j = 0; (j < ntx); j++)
404 sprintf(buf, "%g", mat[i].axis_x[j]);
405 xtick[j] = strdup(buf);
407 ps_strfont(ps, psr->X.tickfont, psr->X.tickfontsize);
408 for (x = 0; (x < ntx); x++)
410 xx = xx00 + (x + offset_x)*psr->xboxsize;
411 if ( ( bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.major) ||
412 (psr->X.first && (x == 0))) &&
413 ( (i == 0) || box_do_all_x_maj_ticks(psr) ) )
415 /* Longer tick marks */
416 ps_line (ps, xx, yy00, xx, yy00-psr->X.majorticklen);
417 /* Plot label on lowest graph only */
421 yy00-DDD-psr->X.majorticklen-psr->X.tickfontsize*0.8,
425 else if (bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.minor) &&
426 ( (i == 0) || box_do_all_x_min_ticks(psr) ) )
428 /* Shorter tick marks */
429 ps_line(ps, xx, yy00, xx, yy00-psr->X.minorticklen);
431 else if (bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.major) )
433 /* Even shorter marks, only each X.major */
434 ps_line(ps, xx, yy00, xx, yy00-(psr->boxspacing/2));
437 ps_strfont(ps, psr->Y.tickfont, psr->Y.tickfontsize);
439 for (j = 0; (j < nty); j++)
441 sprintf(buf, "%g", mat[i].axis_y[j]);
442 ytick[j] = strdup(buf);
445 for (y = 0; (y < nty); y++)
447 yy = yy00 + (y + offset_y)*psr->yboxsize;
448 if (bRmod(mat[i].axis_y[y], psr->Y.offset, psr->Y.major) ||
449 (psr->Y.first && (y == 0)))
452 strlength = max(strlength, (int)strlen(ytick[y]));
453 ps_line (ps, xx00, yy, xx00-psr->Y.majorticklen, yy);
454 ps_ctext(ps, xx00-psr->Y.majorticklen-DDD,
455 yy-psr->Y.tickfontsize/3.0, ytick[y], eXRight);
457 else if (bRmod(mat[i].axis_y[y], psr->Y.offset, psr->Y.minor) )
460 ps_line(ps, xx00, yy, xx00-psr->Y.minorticklen, yy);
466 /* Label on Y-axis */
467 if (!psr->bYonce || i == nmat/2)
469 if (strlen(psr->Y.label) > 0)
471 mylab = psr->Y.label;
475 mylab = mat[i].label_y;
477 if (strlen(mylab) > 0)
479 ps_strfont(ps, psr->Y.font, psr->Y.fontsize);
481 xxx = x0-psr->X.majorticklen-psr->X.tickfontsize*strlength-DDD;
482 ps_ctext(ps, yy00+box_height(&mat[i], psr)/2.0, 612.5-xxx,
488 yy00 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
490 /* Label on X-axis */
491 if (strlen(psr->X.label) > 0)
493 mylab = psr->X.label;
497 mylab = mat[0].label_x;
499 if (strlen(mylab) > 0)
501 ps_strfont(ps, psr->X.font, psr->X.fontsize);
502 ps_ctext(ps, x0+w/2, y0-DDD-psr->X.majorticklen-psr->X.tickfontsize*FUDGE-
503 psr->X.fontsize, mylab, eXCenter);
507 static void draw_zerolines(t_psdata out, real x0, real y0, real w,
508 int nmat, t_matrix mat[], t_psrec *psr)
510 real xx, yy, dy, xx00, yy00;
515 ps_linewidth(out, psr->zerolinewidth);
516 for (i = 0; (i < nmat); i++)
518 dy = box_height(&(mat[i]), psr);
519 /* mat[i].axis_x and _y were already set by draw_boxes */
520 if (psr->X.lineatzero)
522 ps_rgb(out, linecolors[psr->X.lineatzero]);
523 for (x = 0; (x < mat[i].nx); x++)
525 xx = xx00+(x+0.7)*psr->xboxsize;
526 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
527 if (x != 0 && x < mat[i].nx-1 &&
528 abs(mat[i].axis_x[x]) <
529 0.1*abs(mat[i].axis_x[x+1]-mat[i].axis_x[x]) )
531 ps_line (out, xx, yy00, xx, yy00+dy+2);
535 if (psr->Y.lineatzero)
537 ps_rgb(out, linecolors[psr->Y.lineatzero]);
538 for (y = 0; (y < mat[i].ny); y++)
540 yy = yy00+(y+0.7)*psr->yboxsize;
541 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
542 if (y != 0 && y < mat[i].ny-1 &&
543 abs(mat[i].axis_y[y]) <
544 0.1*abs(mat[i].axis_y[y+1]-mat[i].axis_y[y]) )
546 ps_line (out, xx00, yy, xx00+w+2, yy);
550 yy00 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
554 static void box_dim(int nmat, t_matrix mat[], t_matrix *mat2, t_psrec *psr,
555 int elegend, gmx_bool bFrame,
556 real *w, real *h, real *dw, real *dh)
559 real ww, hh, dww, dhh;
565 for (i = 0; (i < nmat); i++)
567 ww = max(ww, mat[i].nx*psr->xboxsize);
568 hh += box_height(&(mat[i]), psr);
569 maxytick = max(maxytick, mat[i].nx);
573 if (mat[0].label_y[0])
575 dww += 2.0*(psr->Y.fontsize+DDD);
577 if (psr->Y.major > 0)
579 dww += psr->Y.majorticklen + DDD +
580 psr->Y.tickfontsize*(log(maxytick)/log(10.0));
582 else if (psr->Y.minor > 0)
584 dww += psr->Y.minorticklen;
587 if (mat[0].label_x[0])
589 dhh += psr->X.fontsize+2*DDD;
591 if ( /* fool emacs auto-indent */
592 (elegend == elBoth && (mat[0].legend[0] || mat2[0].legend[0])) ||
593 (elegend == elFirst && mat[0].legend[0]) ||
594 (elegend == elSecond && mat2[0].legend[0]) )
596 dhh += 2*(psr->legfontsize*FUDGE+2*DDD);
600 dhh += psr->legfontsize*FUDGE+2*DDD;
602 if (psr->X.major > 0)
604 dhh += psr->X.tickfontsize*FUDGE+2*DDD+psr->X.majorticklen;
606 else if (psr->X.minor > 0)
608 dhh += psr->X.minorticklen;
611 hh += (nmat-1)*box_dh(psr);
612 hh += box_dh_top(TRUE, psr);
615 hh += (nmat-1)*box_dh_top(FALSE, psr);
624 int add_maps(t_mapping **newmap,
625 int nmap1, t_mapping map1[], int nmap2, t_mapping map2[])
627 static char mapper[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<.>/?";
632 nsymbols = strlen(mapper);
634 if (nmap > nsymbols*nsymbols)
636 gmx_fatal(FARGS, "Not enough symbols to merge the two colormaps\n");
638 printf("Combining colormaps of %d and %d elements into one of %d elements\n",
641 for (j = 0; j < nmap1; j++)
643 map[j].code.c1 = mapper[j % nsymbols];
646 map[j].code.c2 = mapper[j/nsymbols];
648 map[j].rgb.r = map1[j].rgb.r;
649 map[j].rgb.g = map1[j].rgb.g;
650 map[j].rgb.b = map1[j].rgb.b;
651 map[j].desc = map1[j].desc;
653 for (j = 0; j < nmap2; j++)
656 map[k].code.c1 = mapper[k % nsymbols];
659 map[k].code.c2 = mapper[k/nsymbols];
661 map[k].rgb.r = map2[j].rgb.r;
662 map[k].rgb.g = map2[j].rgb.g;
663 map[k].rgb.b = map2[j].rgb.b;
664 map[k].desc = map2[j].desc;
671 void xpm_mat(const char *outf, int nmat, t_matrix *mat, t_matrix *mat2,
672 gmx_bool bDiag, gmx_bool bFirstDiag)
676 int i, j, k, x, y, col;
678 t_mapping *map = NULL;
680 out = gmx_ffopen(outf, "w");
682 for (i = 0; i < nmat; i++)
684 if (!mat2 || !diff_maps(mat[i].nmap, mat[i].map, mat2[i].nmap, mat2[i].map))
686 write_xpm_m(out, mat[0]);
690 nmap = add_maps(&map, mat[i].nmap, mat[i].map, mat2[i].nmap, mat2[i].map);
691 for (x = 0; (x < mat[i].nx); x++)
693 for (y = 0; (y < mat[i].nx); y++)
695 if ((x < y) || ((x == y) && bFirstDiag)) /* upper left -> map1 */
697 col = mat[i].matrix[x][y];
699 else /* lower right -> map2 */
701 col = mat[i].nmap+mat[i].matrix[x][y];
703 if ((bDiag) || (x != y))
705 mat[i].matrix[x][y] = col;
709 mat[i].matrix[x][y] = 0;
716 if (mat2 && (strcmp(mat[i].title, mat2[i].title) != 0))
718 sprintf(mat[i].title+strlen(mat[i].title), " / %s", mat2[i].title);
720 if (mat2 && (strcmp(mat[i].legend, mat2[i].legend) != 0))
722 sprintf(mat[i].legend+strlen(mat[i].legend), " / %s", mat2[i].legend);
724 write_xpm_m(out, mat[i]);
730 static void tick_spacing(int n, real axis[], real offset, char axisnm,
731 real *major, real *minor)
734 gmx_bool bTryAgain, bFive;
735 int i, j, t, f = 0, ten;
737 real major_fact[NFACT] = {5, 4, 2, 1};
738 real minor_fact[NFACT] = {5, 4, 4, 5};
740 /* start with interval between 10 matrix points: */
741 space = max(10*axis[1]-axis[0], axis[min(10, n-1)]-axis[0]);
742 /* get power of 10 */
743 ten = (int)ceil(log(space)/log(10))-1;
745 for (t = ten+2; t > ten-3 && bTryAgain; t--)
747 for (f = 0; f < NFACT && bTryAgain; f++)
749 space = pow(10, t) * major_fact[f];
750 /* count how many ticks we would get: */
752 for (j = 0; j < n; j++)
754 if (bRmod(axis[j], offset, space) )
759 /* do we have a reasonable number of ticks ? */
760 bTryAgain = (i > min(10, n-1)) || (i < 5);
765 space = max(10*axis[1]-axis[0], axis[min(10, n-1)]-axis[0]);
766 fprintf(stderr, "Auto tick spacing failed for %c-axis, guessing %g\n",
770 *minor = space / minor_fact[f-1];
771 fprintf(stderr, "Auto tick spacing for %c-axis: major %g, minor %g\n",
772 axisnm, *major, *minor);
775 void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
776 gmx_bool bFrame, gmx_bool bDiag, gmx_bool bFirstDiag,
777 gmx_bool bTitle, gmx_bool bTitleOnce, gmx_bool bYonce, int elegend,
778 real size, real boxx, real boxy, const char *m2p, const char *m2pout,
782 char buf[256], *legend;
786 int i, j, x, y, col, leg = 0;
789 int nmap1 = 0, nmap2 = 0, leg_nmap;
790 t_mapping *map1 = NULL, *map2 = NULL, *leg_map;
791 gmx_bool bMap1, bNextMap1, bDiscrete;
794 libm2p = m2p ? gmxlibfn(m2p) : m2p;
795 get_params(libm2p, m2pout, &psrec);
799 if (psr->X.major <= 0)
801 tick_spacing((mat[0].flags & MAT_SPATIAL_X) ? mat[0].nx + 1 : mat[0].nx,
802 mat[0].axis_x, psr->X.offset, 'X',
803 &(psr->X.major), &(psr->X.minor) );
805 if (psr->X.minor <= 0)
807 psr->X.minor = psr->X.major / 2;
809 if (psr->Y.major <= 0)
811 tick_spacing((mat[0].flags & MAT_SPATIAL_Y) ? mat[0].ny + 1 : mat[0].ny,
812 mat[0].axis_y, psr->Y.offset, 'Y',
813 &(psr->Y.major), &(psr->Y.minor) );
815 if (psr->Y.minor <= 0)
817 psr->Y.minor = psr->Y.major / 2;
822 psr->xboxsize = boxx;
823 psr->yboxsize = boxx;
827 psr->yboxsize = boxy;
830 if (psr->xboxsize == 0)
832 psr->xboxsize = size/mat[0].nx;
833 printf("Set the x-size of the box to %.3f\n", psr->xboxsize);
835 if (psr->yboxsize == 0)
837 psr->yboxsize = size/mat[0].nx;
838 printf("Set the y-size of the box to %.3f\n", psr->yboxsize);
842 for (i = 0; (i < nmat); i++)
844 if (mat[i].nmap > nmap1)
853 printf("Selected legend of matrix # %d for display\n", leg);
858 for (i = 0; (i < nmat); i++)
860 if (mat2[i].nmap > nmap2)
862 nmap2 = mat2[i].nmap;
869 printf("Selected legend of matrix # %d for second display\n", leg);
872 if ( (mat[0].legend[0] == 0) && psr->legend)
874 strcpy(mat[0].legend, psr->leglabel);
877 bTitle = bTitle && mat[nmat-1].title[0];
878 bTitleOnce = bTitleOnce && mat[nmat-1].title[0];
879 psr->bTitle = bTitle;
880 psr->bTitleOnce = bTitleOnce;
881 psr->bYonce = bYonce;
883 /* Set up size of box for nice colors */
884 box_dim(nmat, mat, mat2, psr, elegend, bFrame, &w, &h, &dw, &dh);
886 /* Set up bounding box */
900 out = ps_open(outf, 0, 0, x, y);
901 ps_linewidth(out, psr->linewidth);
902 ps_init_rgb_box(out, psr->xboxsize, psr->yboxsize);
903 ps_init_rgb_nbox(out, psr->xboxsize, psr->yboxsize);
904 ps_translate(out, psr->xoffs, psr->yoffs);
908 ps_comment(out, "Here starts the BOX drawing");
909 draw_boxes(out, x0, y0, w, nmat, mat, psr);
912 for (i = 0; (i < nmat); i++)
914 if (bTitle || (bTitleOnce && i == nmat-1) )
916 /* Print title, if any */
918 ps_strfont(out, psr->titfont, psr->titfontsize);
919 if (!mat2 || (strcmp(mat[i].title, mat2[i].title) == 0))
921 strcpy(buf, mat[i].title);
925 sprintf(buf, "%s / %s", mat[i].title, mat2[i].title);
927 ps_ctext(out, x0+w/2, y0+box_height(&(mat[i]), psr)+psr->titfontsize,
930 sprintf(buf, "Here starts the filling of box #%d", i);
931 ps_comment(out, buf);
932 for (x = 0; (x < mat[i].nx); x++)
937 xx = x0+x*psr->xboxsize;
938 ps_moveto(out, xx, y0);
940 bMap1 = (!mat2 || (x < y || (x == y && bFirstDiag)));
941 if ((bDiag) || (x != y))
943 col = mat[i].matrix[x][y];
949 for (nexty = 1; (nexty <= mat[i].ny); nexty++)
951 bNextMap1 = (!mat2 || (x < nexty || (x == nexty && bFirstDiag)));
952 /* TRUE: upper left -> map1 */
953 /* FALSE: lower right -> map2 */
954 if ((nexty == mat[i].ny) || (!bDiag && (x == nexty)))
960 nextcol = mat[i].matrix[x][nexty];
962 if ( (nexty == mat[i].ny) || (col != nextcol) || (bMap1 != bNextMap1) )
968 ps_rgb_nbox(out, &(mat[i].map[col].rgb), nexty-y);
972 ps_rgb_nbox(out, &(mat2[i].map[col].rgb), nexty-y);
977 ps_moverel(out, 0, psr->yboxsize);
985 y0 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
988 if (psr->X.lineatzero || psr->Y.lineatzero)
990 /* reset y0 for first box */
992 ps_comment(out, "Here starts the zero lines drawing");
993 draw_zerolines(out, x0, y0, w, nmat, mat, psr);
996 if (elegend != elNone)
998 ps_comment(out, "Now it's legend time!");
999 ps_linewidth(out, psr->linewidth);
1000 if (mat2 == NULL || elegend != elSecond)
1002 bDiscrete = mat[0].bDiscrete;
1003 legend = mat[0].legend;
1009 bDiscrete = mat2[0].bDiscrete;
1010 legend = mat2[0].legend;
1016 leg_discrete(out, psr->legfontsize, DDD, legend,
1017 psr->legfontsize, psr->legfont, leg_nmap, leg_map);
1021 if (elegend != elBoth)
1023 leg_continuous(out, x0+w/2, w/2, DDD, legend,
1024 psr->legfontsize, psr->legfont, leg_nmap, leg_map,
1029 leg_bicontinuous(out, x0+w/2, w, DDD, mat[0].legend, mat2[0].legend,
1030 psr->legfontsize, psr->legfont, nmap1, map1, nmap2, map2);
1033 ps_comment(out, "Were there, dude");
1039 void make_axis_labels(int nmat, t_matrix *mat)
1043 for (i = 0; (i < nmat); i++)
1045 /* Make labels for x axis */
1046 if (mat[i].axis_x == NULL)
1048 snew(mat[i].axis_x, mat[i].nx);
1049 for (j = 0; (j < mat[i].nx); j++)
1051 mat[i].axis_x[j] = j;
1054 /* Make labels for y axis */
1055 if (mat[i].axis_y == NULL)
1057 snew(mat[i].axis_y, mat[i].ny);
1058 for (j = 0; (j < mat[i].ny); j++)
1060 mat[i].axis_y[j] = j;
1066 void prune_mat(int nmat, t_matrix *mat, t_matrix *mat2, int skip)
1068 int i, x, y, xs, ys;
1070 for (i = 0; i < nmat; i++)
1072 fprintf(stderr, "converting %dx%d matrix to %dx%d\n",
1073 mat[i].nx, mat[i].ny,
1074 (mat[i].nx+skip-1)/skip, (mat[i].ny+skip-1)/skip);
1075 /* walk through matrix */
1077 for (x = 0; (x < mat[i].nx); x++)
1081 mat[i].axis_x[xs] = mat[i].axis_x[x];
1084 mat2[i].axis_x[xs] = mat2[i].axis_x[x];
1087 for (y = 0; (y < mat[i].ny); y++)
1091 mat[i].axis_y[ys] = mat[i].axis_y[y];
1094 mat2[i].axis_y[ys] = mat2[i].axis_y[y];
1099 mat[i].matrix[xs][ys] = mat[i].matrix[x][y];
1102 mat2[i].matrix[xs][ys] = mat2[i].matrix[x][y];
1110 /* adjust parameters */
1111 mat[i].nx = (mat[i].nx+skip-1)/skip;
1112 mat[i].ny = (mat[i].ny+skip-1)/skip;
1115 mat2[i].nx = (mat2[i].nx+skip-1)/skip;
1116 mat2[i].ny = (mat2[i].ny+skip-1)/skip;
1121 void zero_lines(int nmat, t_matrix *mat, t_matrix *mat2)
1126 for (i = 0; i < nmat; i++)
1128 for (m = 0; m < (mat2 ? 2 : 1); m++)
1138 for (x = 0; x < mats[i].nx-1; x++)
1140 if (abs(mats[i].axis_x[x+1]) < 1e-5)
1142 for (y = 0; y < mats[i].ny; y++)
1144 mats[i].matrix[x][y] = 0;
1148 for (y = 0; y < mats[i].ny-1; y++)
1150 if (abs(mats[i].axis_y[y+1]) < 1e-5)
1152 for (x = 0; x < mats[i].nx; x++)
1154 mats[i].matrix[x][y] = 0;
1162 void write_combined_matrix(int ecombine, const char *fn,
1163 int nmat, t_matrix *mat1, t_matrix *mat2,
1164 real *cmin, real *cmax)
1166 int i, j, k, nlevels;
1167 t_mapping *map = NULL;
1169 real **rmat1, **rmat2;
1172 out = gmx_ffopen(fn, "w");
1173 for (k = 0; k < nmat; k++)
1175 if (mat2[k].nx != mat1[k].nx || mat2[k].ny != mat1[k].ny)
1177 gmx_fatal(FARGS, "Size of frame %d in 1st (%dx%d) and 2nd matrix (%dx%d) do"
1178 " not match.\n", k, mat1[k].nx, mat1[k].ny, mat2[k].nx, mat2[k].ny);
1180 printf("Combining two %dx%d matrices\n", mat1[k].nx, mat1[k].ny);
1181 rmat1 = matrix2real(&mat1[k], NULL);
1182 rmat2 = matrix2real(&mat2[k], NULL);
1183 if (NULL == rmat1 || NULL == rmat2)
1185 gmx_fatal(FARGS, "Could not extract real data from %s xpm matrices. Note that, e.g.,\n"
1186 "g_rms and g_mdmat provide such data, but not do_dssp.\n",
1187 (NULL == rmat1 && NULL == rmat2) ? "both" : "one of the" );
1191 for (j = 0; j < mat1[k].ny; j++)
1193 for (i = 0; i < mat1[k].nx; i++)
1197 case ecAdd: rmat1[i][j] += rmat2[i][j]; break;
1198 case ecSub: rmat1[i][j] -= rmat2[i][j]; break;
1199 case ecMult: rmat1[i][j] *= rmat2[i][j]; break;
1200 case ecDiv: rmat1[i][j] /= rmat2[i][j]; break;
1202 gmx_fatal(FARGS, "No such combination rule %d for matrices", ecombine);
1204 rlo = min(rlo, rmat1[i][j]);
1205 rhi = max(rhi, rmat1[i][j]);
1216 nlevels = max(mat1[k].nmap, mat2[k].nmap);
1220 "combination results in uniform matrix (%g), no output\n", rhi);
1223 else if (rlo>=0 || rhi<=0)
1224 write_xpm(out, mat1[k].flags, mat1[k].title, mat1[k].legend,
1225 mat1[k].label_x, mat1[k].label_y,
1226 mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
1227 rmat1, rlo, rhi, rhi<=0?red:white, rhi<=0?white:blue,
1230 write_xpm3(out, mat2[k].flags, mat1[k].title, mat1[k].legend,
1231 mat1[k].label_x, mat1[k].label_y,
1232 mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
1233 rmat1, rlo, 0, rhi, red, white, blue, &nlevels);
1237 write_xpm(out, mat1[k].flags, mat1[k].title, mat1[k].legend,
1238 mat1[k].label_x, mat1[k].label_y,
1239 mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
1240 rmat1, rlo, rhi, white, black, &nlevels);
1246 void do_mat(int nmat, t_matrix *mat, t_matrix *mat2,
1247 gmx_bool bFrame, gmx_bool bZeroLine, gmx_bool bDiag, gmx_bool bFirstDiag, gmx_bool bTitle,
1248 gmx_bool bTitleOnce, gmx_bool bYonce, int elegend,
1249 real size, real boxx, real boxy,
1250 const char *epsfile, const char *xpmfile, const char *m2p,
1251 const char *m2pout, int skip, int mapoffset)
1257 for (k = 0; (k < nmat); k++)
1259 if ((mat2[k].nx != mat[k].nx) || (mat2[k].ny != mat[k].ny))
1261 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",
1262 k, mat2[k].nx, mat2[k].ny, mat[k].nx, mat[k].ny);
1264 for (j = 0; (j < mat[k].ny); j++)
1266 for (i = bFirstDiag ? j+1 : j; (i < mat[k].nx); i++)
1268 mat[k].matrix[i][j] = mat2[k].matrix[i][j];
1273 for (i = 0; (i < nmat); i++)
1275 fprintf(stderr, "Matrix %d is %d x %d\n", i, mat[i].nx, mat[i].ny);
1278 make_axis_labels(nmat, mat);
1282 prune_mat(nmat, mat, mat2, skip);
1287 zero_lines(nmat, mat, mat);
1290 if (epsfile != NULL)
1292 ps_mat(epsfile, nmat, mat, mat2, bFrame, bDiag, bFirstDiag,
1293 bTitle, bTitleOnce, bYonce, elegend,
1294 size, boxx, boxy, m2p, m2pout, mapoffset);
1296 if (xpmfile != NULL)
1298 xpm_mat(xpmfile, nmat, mat, mat2, bDiag, bFirstDiag);
1302 void gradient_map(rvec grad, int nmap, t_mapping map[])
1307 for (i = 0; i < nmap; i++)
1310 map[i].rgb.r = 1-c*(1-grad[XX]);
1311 map[i].rgb.g = 1-c*(1-grad[YY]);
1312 map[i].rgb.b = 1-c*(1-grad[ZZ]);
1316 void gradient_mat(rvec grad, int nmat, t_matrix mat[])
1320 for (m = 0; m < nmat; m++)
1322 gradient_map(grad, mat[m].nmap, mat[m].map);
1326 void rainbow_map(gmx_bool bBlue, int nmap, t_mapping map[])
1331 for (i = 0; i < nmap; i++)
1333 c = (map[i].rgb.r + map[i].rgb.g + map[i].rgb.b)/3;
1342 if (c <= 0.25) /* 0-0.25 */
1345 g = pow(4*c, 0.666);
1348 else if (c <= 0.5) /* 0.25-0.5 */
1352 b = pow(2-4*c, 0.666);
1354 else if (c <= 0.75) /* 0.5-0.75 */
1356 r = pow(4*c-2, 0.666);
1363 g = pow(4-4*c, 0.666);
1372 void rainbow_mat(gmx_bool bBlue, int nmat, t_matrix mat[])
1376 for (m = 0; m < nmat; m++)
1378 rainbow_map(bBlue, mat[m].nmap, mat[m].map);
1382 int gmx_xpm2ps(int argc, char *argv[])
1384 const char *desc[] = {
1385 "[THISMODULE] makes a beautiful color plot of an XPixelMap file.",
1386 "Labels and axis can be displayed, when they are supplied",
1387 "in the correct matrix format.",
1388 "Matrix data may be generated by programs such as [gmx-do_dssp], [gmx-rms] or",
1389 "[gmx-mdmat].[PAR]",
1390 "Parameters are set in the [TT].m2p[tt] file optionally supplied with",
1391 "[TT]-di[tt]. Reasonable defaults are provided. Settings for the [IT]y[it]-axis",
1392 "default to those for the [IT]x[it]-axis. Font names have a defaulting hierarchy:",
1393 "titlefont -> legendfont; titlefont -> (xfont -> yfont -> ytickfont)",
1394 "-> xtickfont, e.g. setting titlefont sets all fonts, setting xfont",
1395 "sets yfont, ytickfont and xtickfont.[PAR]",
1396 "When no [TT].m2p[tt] file is supplied, many settings are taken from",
1397 "command line options. The most important option is [TT]-size[tt],",
1398 "which sets the size of the whole matrix in postscript units.",
1399 "This option can be overridden with the [TT]-bx[tt] and [TT]-by[tt]",
1400 "options (and the corresponding parameters in the [TT].m2p[tt] file),",
1401 "which set the size of a single matrix element.[PAR]",
1402 "With [TT]-f2[tt] a second matrix file can be supplied. Both matrix",
1403 "files will be read simultaneously and the upper left half of the",
1404 "first one ([TT]-f[tt]) is plotted together with the lower right",
1405 "half of the second one ([TT]-f2[tt]). The diagonal will contain",
1406 "values from the matrix file selected with [TT]-diag[tt].",
1407 "Plotting of the diagonal values can be suppressed altogether by",
1408 "setting [TT]-diag[tt] to [TT]none[tt].",
1409 "In this case, a new color map will be generated with",
1410 "a red gradient for negative numbers and a blue for positive.",
1411 "If the color coding and legend labels of both matrices are identical,",
1412 "only one legend will be displayed, else two separate legends are",
1414 "With [TT]-combine[tt], an alternative operation can be selected",
1415 "to combine the matrices. The output range is automatically set",
1416 "to the actual range of the combined matrix. This can be overridden",
1417 "with [TT]-cmin[tt] and [TT]-cmax[tt].[PAR]",
1418 "[TT]-title[tt] can be set to [TT]none[tt] to suppress the title, or to",
1419 "[TT]ylabel[tt] to show the title in the Y-label position (alongside",
1420 "the [IT]y[it]-axis).[PAR]",
1421 "With the [TT]-rainbow[tt] option, dull grayscale matrices can be turned",
1422 "into attractive color pictures.[PAR]",
1423 "Merged or rainbowed matrices can be written to an XPixelMap file with",
1424 "the [TT]-xpm[tt] option."
1428 const char *fn, *epsfile = NULL, *xpmfile = NULL;
1429 int i, nmat, nmat2, etitle, elegend, ediag, erainbow, ecombine;
1430 t_matrix *mat = NULL, *mat2 = NULL;
1431 gmx_bool bTitle, bTitleOnce, bDiag, bFirstDiag, bGrad;
1432 static gmx_bool bFrame = TRUE, bZeroLine = FALSE, bYonce = FALSE, bAdd = FALSE;
1433 static real size = 400, boxx = 0, boxy = 0, cmin = 0, cmax = 0;
1434 static rvec grad = {0, 0, 0};
1436 etSel, etTop, etOnce, etYlabel, etNone, etNR
1438 const char *title[] = { NULL, "top", "once", "ylabel", "none", NULL };
1439 /* MUST correspond to enum elXxx as defined at top of file */
1440 const char *legend[] = { NULL, "both", "first", "second", "none", NULL };
1442 edSel, edFirst, edSecond, edNone, edNR
1444 const char *diag[] = { NULL, "first", "second", "none", NULL };
1446 erSel, erNo, erBlue, erRed, erNR
1448 const char *rainbow[] = { NULL, "no", "blue", "red", NULL };
1449 /* MUST correspond to enum ecXxx as defined at top of file */
1450 const char *combine[] = {
1451 NULL, "halves", "add", "sub", "mult", "div", NULL
1453 static int skip = 1, mapoffset = 0;
1455 { "-frame", FALSE, etBOOL, {&bFrame},
1456 "Display frame, ticks, labels, title and legend" },
1457 { "-title", FALSE, etENUM, {title}, "Show title at" },
1458 { "-yonce", FALSE, etBOOL, {&bYonce}, "Show y-label only once" },
1459 { "-legend", FALSE, etENUM, {legend}, "Show legend" },
1460 { "-diag", FALSE, etENUM, {diag}, "Diagonal" },
1461 { "-size", FALSE, etREAL, {&size},
1462 "Horizontal size of the matrix in ps units" },
1463 { "-bx", FALSE, etREAL, {&boxx},
1464 "Element x-size, overrides [TT]-size[tt] (also y-size when [TT]-by[tt] is not set)" },
1465 { "-by", FALSE, etREAL, {&boxy}, "Element y-size" },
1466 { "-rainbow", FALSE, etENUM, {rainbow},
1467 "Rainbow colors, convert white to" },
1468 { "-gradient", FALSE, etRVEC, {grad},
1469 "Re-scale colormap to a smooth gradient from white {1,1,1} to {r,g,b}" },
1470 { "-skip", FALSE, etINT, {&skip},
1471 "only write out every nr-th row and column" },
1472 { "-zeroline", FALSE, etBOOL, {&bZeroLine},
1473 "insert line in [TT].xpm[tt] matrix where axis label is zero"},
1474 { "-legoffset", FALSE, etINT, {&mapoffset},
1475 "Skip first N colors from [TT].xpm[tt] file for the legend" },
1476 { "-combine", FALSE, etENUM, {combine}, "Combine two matrices" },
1477 { "-cmin", FALSE, etREAL, {&cmin}, "Minimum for combination output" },
1478 { "-cmax", FALSE, etREAL, {&cmax}, "Maximum for combination output" }
1480 #define NPA asize(pa)
1482 { efXPM, "-f", NULL, ffREAD },
1483 { efXPM, "-f2", "root2", ffOPTRD },
1484 { efM2P, "-di", NULL, ffLIBOPTRD },
1485 { efM2P, "-do", "out", ffOPTWR },
1486 { efEPS, "-o", NULL, ffOPTWR },
1487 { efXPM, "-xpm", NULL, ffOPTWR }
1489 #define NFILE asize(fnm)
1491 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW,
1492 NFILE, fnm, NPA, pa,
1493 asize(desc), desc, 0, NULL, &oenv))
1498 etitle = nenum(title);
1499 elegend = nenum(legend);
1500 ediag = nenum(diag);
1501 erainbow = nenum(rainbow);
1502 ecombine = nenum(combine);
1503 bGrad = opt2parg_bSet("-gradient", NPA, pa);
1504 for (i = 0; i < DIM; i++)
1506 if (grad[i] < 0 || grad[i] > 1)
1508 gmx_fatal(FARGS, "RGB value %g out of range (0.0-1.0)", grad[i]);
1517 epsfile = ftp2fn_null(efEPS, NFILE, fnm);
1518 xpmfile = opt2fn_null("-xpm", NFILE, fnm);
1519 if (epsfile == NULL && xpmfile == NULL)
1521 if (ecombine != ecHalves)
1523 xpmfile = opt2fn("-xpm", NFILE, fnm);
1527 epsfile = ftp2fn(efEPS, NFILE, fnm);
1530 if (ecombine != ecHalves && epsfile)
1533 "WARNING: can only write result of arithmetic combination "
1534 "of two matrices to .xpm file\n"
1535 " file %s will not be written\n", epsfile);
1539 bDiag = ediag != edNone;
1540 bFirstDiag = ediag != edSecond;
1542 fn = opt2fn("-f", NFILE, fnm);
1543 nmat = read_xpm_matrix(fn, &mat);
1544 fprintf(stderr, "There %s %d matri%s in %s\n", (nmat > 1) ? "are" : "is", nmat, (nmat > 1) ? "ces" : "x", fn);
1545 fn = opt2fn_null("-f2", NFILE, fnm);
1548 nmat2 = read_xpm_matrix(fn, &mat2);
1549 fprintf(stderr, "There %s %d matri%s in %s\n", (nmat2 > 1) ? "are" : "is", nmat2, (nmat2 > 1) ? "ces" : "x", fn);
1552 fprintf(stderr, "Different number of matrices, using the smallest number.\n");
1553 nmat = nmat2 = min(nmat, nmat2);
1558 if (ecombine != ecHalves)
1561 "WARNING: arithmetic matrix combination selected (-combine), "
1562 "but no second matrix (-f2) supplied\n"
1563 " no matrix combination will be performed\n");
1568 bTitle = etitle == etTop;
1569 bTitleOnce = etitle == etOnce;
1570 if (etitle == etYlabel)
1572 for (i = 0; (i < nmat); i++)
1574 strcpy(mat[i].label_y, mat[i].title);
1577 strcpy(mat2[i].label_y, mat2[i].title);
1583 gradient_mat(grad, nmat, mat);
1586 gradient_mat(grad, nmat2, mat2);
1589 if (erainbow != erNo)
1591 rainbow_mat(erainbow == erBlue, nmat, mat);
1594 rainbow_mat(erainbow == erBlue, nmat2, mat2);
1598 if ((mat2 == NULL) && (elegend != elNone))
1603 if (ecombine && ecombine != ecHalves)
1605 write_combined_matrix(ecombine, xpmfile, nmat, mat, mat2,
1606 opt2parg_bSet("-cmin", NPA, pa) ? &cmin : NULL,
1607 opt2parg_bSet("-cmax", NPA, pa) ? &cmax : NULL);
1611 do_mat(nmat, mat, mat2, bFrame, bZeroLine, bDiag, bFirstDiag,
1612 bTitle, bTitleOnce, bYonce,
1613 elegend, size, boxx, boxy, epsfile, xpmfile,
1614 opt2fn_null("-di", NFILE, fnm), opt2fn_null("-do", NFILE, fnm), skip,
1618 view_all(oenv, NFILE, fnm);