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