src/gromacs/fileio/matio.cpp

   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, 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 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <stdio.h>
  42 #include <ctype.h>
  43
  44 #include <algorithm>
  45
  46 #include "sysstuff.h"
  47 #include "futil.h"
  48 #include "string2.h"
  49 #include "macros.h"
  50 #include "smalloc.h"
  51 #include "gmx_fatal.h"
  52 #include "matio.h"
  53 #include "statutil.h"
  54 #include "gmxfio.h"
  55 #include "maths.h"
  56 #include "copyrite.h"
  57
  58 #include "gromacs/utility/exceptions.h"
  59 #include "gromacs/utility/programinfo.h"
  60
  61 #define round(a) (int)(a+0.5)
  62
  63 static const char mapper[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<.>/?";
  64 #define NMAP (long int)strlen(mapper)
  65
  66 #define MAX_XPM_LINELENGTH 4096
  67
  68 real **mk_matrix(int nx, int ny, gmx_bool b1D)
  69 {
  70     int    i;
  71     real **m;
  72
  73     snew(m, nx);
  74     if (b1D)
  75     {
  76         snew(m[0], nx*ny);
  77     }
  78
  79     for (i = 0; (i < nx); i++)
  80     {
  81         if (b1D)
  82         {
  83             m[i] = &(m[0][i*ny]);
  84         }
  85         else
  86         {
  87             snew(m[i], ny);
  88         }
  89     }
  90
  91     return m;
  92 }
  93
  94 void done_matrix(int nx, real ***m)
  95 {
  96     int i;
  97
  98     for (i = 0; (i < nx); i++)
  99     {
 100         sfree((*m)[i]);
 101     }
 102     sfree(*m);
 103     *m = NULL;
 104 }
 105
 106 void clear_matrix(int nx, int ny, real **m)
 107 {
 108     int x, y;
 109
 110     for (x = 0; x < nx; x++)
 111     {
 112         for (y = 0; y < ny; y++)
 113         {
 114             m[x][y] = 0;
 115         }
 116     }
 117 }
 118
 119 gmx_bool matelmt_cmp(t_xpmelmt e1, t_xpmelmt e2)
 120 {
 121     return (e1.c1 == e2.c1) && (e1.c2 == e2.c2);
 122 }
 123
 124 t_matelmt searchcmap(int n, t_mapping map[], t_xpmelmt c)
 125 {
 126     int i;
 127
 128     for (i = 0; (i < n); i++)
 129     {
 130         if (matelmt_cmp(map[i].code, c))
 131         {
 132             return i;
 133         }
 134     }
 135
 136     return -1;
 137 }
 138
 139 int getcmap(FILE *in, const char *fn, t_mapping **map)
 140 {
 141     int        i, n;
 142     char       line[STRLEN];
 143     char       code[STRLEN], desc[STRLEN];
 144     double     r, g, b;
 145     t_mapping *m;
 146
 147     if (fgets2(line, STRLEN-1, in) == NULL)
 148     {
 149         gmx_fatal(FARGS, "Not enough lines in colormap file %s"
 150                   "(just wanted to read number of entries)", fn);
 151     }
 152     sscanf(line, "%d", &n);
 153     snew(m, n);
 154     for (i = 0; (i < n); i++)
 155     {
 156         if (fgets2(line, STRLEN-1, in) == NULL)
 157         {
 158             gmx_fatal(FARGS, "Not enough lines in colormap file %s"
 159                       "(should be %d, found only %d)", fn, n+1, i);
 160         }
 161         sscanf(line, "%s%s%lf%lf%lf", code, desc, &r, &g, &b);
 162         m[i].code.c1 = code[0];
 163         m[i].code.c2 = 0;
 164         m[i].desc    = strdup(desc);
 165         m[i].rgb.r   = r;
 166         m[i].rgb.g   = g;
 167         m[i].rgb.b   = b;
 168     }
 169     *map = m;
 170
 171     return n;
 172 }
 173
 174 int readcmap(const char *fn, t_mapping **map)
 175 {
 176     FILE      *in;
 177     int        n;
 178
 179     in = libopen(fn);
 180     n  = getcmap(in, fn, map);
 181     gmx_fio_fclose(in);
 182
 183     return n;
 184 }
 185
 186 void printcmap(FILE *out, int n, t_mapping map[])
 187 {
 188     int i;
 189
 190     fprintf(out, "%d\n", n);
 191     for (i = 0; (i < n); i++)
 192     {
 193         fprintf(out, "%c%c  %20s  %10g  %10g  %10g\n",
 194                 map[i].code.c1 ? map[i].code.c1 : ' ',
 195                 map[i].code.c2 ? map[i].code.c2 : ' ',
 196                 map[i].desc, map[i].rgb.r, map[i].rgb.g, map[i].rgb.b);
 197     }
 198 }
 199
 200 void writecmap(const char *fn, int n, t_mapping map[])
 201 {
 202     FILE *out;
 203
 204     out = gmx_fio_fopen(fn, "w");
 205     printcmap(out, n, map);
 206     gmx_fio_fclose(out);
 207 }
 208
 209 static char *fgetline(char **line, int llmax, int *llalloc, FILE *in)
 210 {
 211     char *fg;
 212
 213     if (llmax > *llalloc)
 214     {
 215         srenew(*line, llmax+1);
 216         *llalloc = llmax;
 217     }
 218     fg = fgets(*line, llmax, in);
 219     trim(*line);
 220
 221     return fg;
 222 }
 223
 224 void skipstr(char *line)
 225 {
 226     int i, c;
 227
 228     ltrim(line);
 229     c = 0;
 230     while ((line[c] != ' ') && (line[c] != '\0'))
 231     {
 232         c++;
 233     }
 234     i = c;
 235     while (line[c] != '\0')
 236     {
 237         line[c-i] = line[c];
 238         c++;
 239     }
 240     line[c-i] = '\0';
 241 }
 242
 243 char *line2string(char **line)
 244 {
 245     int i;
 246
 247     if (*line != NULL)
 248     {
 249         while (((*line)[0] != '\"' ) && ( (*line)[0] != '\0' ))
 250         {
 251             (*line)++;
 252         }
 253
 254         if ((*line)[0] != '\"')
 255         {
 256             return NULL;
 257         }
 258         (*line)++;
 259
 260         i = 0;
 261         while (( (*line)[i] != '\"' ) && ( (*line)[i] != '\0' ))
 262         {
 263             i++;
 264         }
 265
 266         if ((*line)[i] != '\"')
 267         {
 268             *line = NULL;
 269         }
 270         else
 271         {
 272             (*line)[i] = 0;
 273         }
 274     }
 275
 276     return *line;
 277 }
 278
 279 void parsestring(char *line, const char *label, char *string)
 280 {
 281     if (strstr(line, label))
 282     {
 283         if (strstr(line, label) < strchr(line, '\"'))
 284         {
 285             line2string(&line);
 286             strcpy(string, line);
 287         }
 288     }
 289 }
 290
 291 void read_xpm_entry(FILE *in, t_matrix *mm)
 292 {
 293     t_mapping   *map;
 294     char        *line_buf = NULL, *line = NULL, *str, buf[256] = {0};
 295     int          i, m, col_len, nch = 0, n_axis_x, n_axis_y, llmax;
 296     int          llalloc = 0;
 297     unsigned int r, g, b;
 298     double       u;
 299     gmx_bool     bGetOnWithIt, bSetLine;
 300     t_xpmelmt    c;
 301
 302     mm->flags      = 0;
 303     mm->title[0]   = 0;
 304     mm->legend[0]  = 0;
 305     mm->label_x[0] = 0;
 306     mm->label_y[0] = 0;
 307     mm->matrix     = NULL;
 308     mm->axis_x     = NULL;
 309     mm->axis_y     = NULL;
 310     mm->bDiscrete  = FALSE;
 311
 312     llmax = STRLEN;
 313
 314     while ((NULL != fgetline(&line_buf, llmax, &llalloc, in)) &&
 315            (strncmp(line_buf, "static", 6) != 0))
 316     {
 317         line = line_buf;
 318         parsestring(line, "title", (mm->title));
 319         parsestring(line, "legend", (mm->legend));
 320         parsestring(line, "x-label", (mm->label_x));
 321         parsestring(line, "y-label", (mm->label_y));
 322         parsestring(line, "type", buf);
 323     }
 324
 325     if (!line || strncmp(line, "static", 6) != 0)
 326     {
 327         gmx_input("Invalid XPixMap");
 328     }
 329
 330     if (buf[0] && (gmx_strcasecmp(buf, "Discrete") == 0))
 331     {
 332         mm->bDiscrete = TRUE;
 333     }
 334
 335     if (debug)
 336     {
 337         fprintf(debug, "%s %s %s %s\n",
 338                 mm->title, mm->legend, mm->label_x, mm->label_y);
 339     }
 340
 341     /* Read sizes */
 342     bGetOnWithIt = FALSE;
 343     while (!bGetOnWithIt && (NULL != fgetline(&line_buf, llmax, &llalloc, in)))
 344     {
 345         line = line_buf;
 346         while (( line[0] != '\"' ) && ( line[0] != '\0' ))
 347         {
 348             line++;
 349         }
 350
 351         if  (line[0] == '\"')
 352         {
 353             line2string(&line);
 354             sscanf(line, "%d %d %d %d", &(mm->nx), &(mm->ny), &(mm->nmap), &nch);
 355             if (nch > 2)
 356             {
 357                 gmx_fatal(FARGS, "Sorry can only read xpm's with at most 2 caracters per pixel\n");
 358             }
 359             if (mm->nx <= 0 || mm->ny <= 0)
 360             {
 361                 gmx_fatal(FARGS, "Dimensions of xpm-file have to be larger than 0\n");
 362             }
 363             llmax        = std::max(STRLEN, mm->nx+10);
 364             bGetOnWithIt = TRUE;
 365         }
 366     }
 367     if (debug)
 368     {
 369         fprintf(debug, "mm->nx %d mm->ny %d mm->nmap %d nch %d\n",
 370                 mm->nx, mm->ny, mm->nmap, nch);
 371     }
 372     if (nch == 0)
 373     {
 374         gmx_fatal(FARGS, "Number of characters per pixel not found in xpm\n");
 375     }
 376
 377     /* Read color map */
 378     snew(map, mm->nmap);
 379     m = 0;
 380     while ((m < mm->nmap) && (NULL != fgetline(&line_buf, llmax, &llalloc, in)))
 381     {
 382         line = strchr(line_buf, '\"');
 383         if  (line)
 384         {
 385             line++;
 386             /* Read xpm color map entry */
 387             map[m].code.c1 = line[0];
 388             if (nch == 1)
 389             {
 390                 map[m].code.c2 = 0;
 391             }
 392             else
 393             {
 394                 map[m].code.c2 = line[1];
 395             }
 396             line += nch;
 397             str   = strchr(line, '#');
 398             if (str)
 399             {
 400                 str++;
 401                 col_len = 0;
 402                 while (isxdigit(str[col_len]))
 403                 {
 404                     col_len++;
 405                 }
 406                 if (col_len == 6)
 407                 {
 408                     sscanf(line, "%*s #%2x%2x%2x", &r, &g, &b);
 409                     map[m].rgb.r = r/255.0;
 410                     map[m].rgb.g = g/255.0;
 411                     map[m].rgb.b = b/255.0;
 412                 }
 413                 else if (col_len == 12)
 414                 {
 415                     sscanf(line, "%*s #%4x%4x%4x", &r, &g, &b);
 416                     map[m].rgb.r = r/65535.0;
 417                     map[m].rgb.g = g/65535.0;
 418                     map[m].rgb.b = b/65535.0;
 419                 }
 420                 else
 421                 {
 422                     gmx_file("Unsupported or invalid colormap in X PixMap");
 423                 }
 424             }
 425             else
 426             {
 427                 str = strchr(line, 'c');
 428                 if (str)
 429                 {
 430                     str += 2;
 431                 }
 432                 else
 433                 {
 434                     gmx_file("Unsupported or invalid colormap in X PixMap");
 435                 }
 436                 fprintf(stderr, "Using white for color \"%s", str);
 437                 map[m].rgb.r = 1;
 438                 map[m].rgb.g = 1;
 439                 map[m].rgb.b = 1;
 440             }
 441             line = strchr(line, '\"');
 442             line++;
 443             line2string(&line);
 444             map[m].desc = strdup(line);
 445             m++;
 446         }
 447     }
 448     if  (m != mm->nmap)
 449     {
 450         gmx_fatal(FARGS, "Number of read colors map entries (%d) does not match the number in the header (%d)", m, mm->nmap);
 451     }
 452     mm->map = map;
 453
 454     /* Read axes, if there are any */
 455     n_axis_x     = 0;
 456     n_axis_y     = 0;
 457     bSetLine     = FALSE;
 458     do
 459     {
 460         if (bSetLine)
 461         {
 462             line = line_buf;
 463         }
 464         bSetLine = TRUE;
 465         if (strstr(line, "x-axis"))
 466         {
 467             line = strstr(line, "x-axis");
 468             skipstr(line);
 469             if (mm->axis_x == NULL)
 470             {
 471                 snew(mm->axis_x, mm->nx + 1);
 472             }
 473             while (sscanf(line, "%lf", &u) == 1)
 474             {
 475                 if (n_axis_x > mm->nx)
 476                 {
 477                     gmx_fatal(FARGS, "Too many x-axis labels in xpm (max %d)", mm->nx);
 478                 }
 479                 else if (n_axis_x == mm->nx)
 480                 {
 481                     mm->flags |= MAT_SPATIAL_X;
 482                 }
 483                 mm->axis_x[n_axis_x] = u;
 484                 n_axis_x++;
 485                 skipstr(line);
 486             }
 487         }
 488         else if (strstr(line, "y-axis"))
 489         {
 490             line = strstr(line, "y-axis");
 491             skipstr(line);
 492             if (mm->axis_y == NULL)
 493             {
 494                 snew(mm->axis_y, mm->ny + 1);
 495             }
 496             while (sscanf(line, "%lf", &u) == 1)
 497             {
 498                 if (n_axis_y > mm->ny)
 499                 {
 500                     gmx_file("Too many y-axis labels in xpm");
 501                 }
 502                 else if (n_axis_y == mm->ny)
 503                 {
 504                     mm->flags |= MAT_SPATIAL_Y;
 505                 }
 506                 mm->axis_y[n_axis_y] = u;
 507                 n_axis_y++;
 508                 skipstr(line);
 509             }
 510         }
 511     }
 512     while ((line[0] != '\"') && (NULL != fgetline(&line_buf, llmax, &llalloc, in)));
 513
 514     /* Read matrix */
 515     snew(mm->matrix, mm->nx);
 516     for (i = 0; i < mm->nx; i++)
 517     {
 518         snew(mm->matrix[i], mm->ny);
 519     }
 520     m        = mm->ny-1;
 521     bSetLine = FALSE;
 522     do
 523     {
 524         if (bSetLine)
 525         {
 526             line = line_buf;
 527         }
 528         bSetLine = TRUE;
 529         if (m%(1+mm->ny/100) == 0)
 530         {
 531             fprintf(stderr, "%3d%%\b\b\b\b", (100*(mm->ny-m))/mm->ny);
 532         }
 533         while ((line[0] != '\"') && (line[0] != '\0'))
 534         {
 535             line++;
 536         }
 537         if (line[0] != '\"')
 538         {
 539             gmx_fatal(FARGS, "Not enough caracters in row %d of the matrix\n", m+1);
 540         }
 541         else
 542         {
 543             line++;
 544             for (i = 0; i < mm->nx; i++)
 545             {
 546                 c.c1 = line[nch*i];
 547                 if (nch == 1)
 548                 {
 549                     c.c2 = 0;
 550                 }
 551                 else
 552                 {
 553                     c.c2 = line[nch*i+1];
 554                 }
 555                 mm->matrix[i][m] = searchcmap(mm->nmap, mm->map, c);
 556             }
 557             m--;
 558         }
 559     }
 560     while ((m >= 0) && (NULL != fgetline(&line_buf, llmax, &llalloc, in)));
 561     if (m >= 0)
 562     {
 563         gmx_incons("Not enough rows in the matrix");
 564     }
 565
 566     sfree(line_buf);
 567 }
 568
 569 int read_xpm_matrix(const char *fnm, t_matrix **matrix)
 570 {
 571     FILE *in;
 572     char *line = NULL;
 573     int   nmat;
 574     int   llalloc = 0;
 575
 576     in = gmx_fio_fopen(fnm, "r");
 577
 578     nmat = 0;
 579     while (NULL != fgetline(&line, STRLEN, &llalloc, in))
 580     {
 581         if (strstr(line, "/* XPM */"))
 582         {
 583             srenew(*matrix, nmat+1);
 584             read_xpm_entry(in, &(*matrix)[nmat]);
 585             nmat++;
 586         }
 587     }
 588     gmx_fio_fclose(in);
 589
 590     if (nmat == 0)
 591     {
 592         gmx_file("Invalid XPixMap");
 593     }
 594
 595     sfree(line);
 596
 597     return nmat;
 598 }
 599
 600 real **matrix2real(t_matrix *matrix, real **mat)
 601 {
 602     t_mapping *map;
 603     double     tmp;
 604     real      *rmap;
 605     int        i, j, nmap;
 606
 607     nmap = matrix->nmap;
 608     map  = matrix->map;
 609     snew(rmap, nmap);
 610
 611     for (i = 0; i < nmap; i++)
 612     {
 613         if ((map[i].desc == NULL) || (sscanf(map[i].desc, "%lf", &tmp) != 1))
 614         {
 615             fprintf(stderr, "Could not convert matrix to reals,\n"
 616                     "color map entry %d has a non-real description: \"%s\"\n",
 617                     i, map[i].desc);
 618             sfree(rmap);
 619             return NULL;
 620         }
 621         rmap[i] = tmp;
 622     }
 623
 624     if (mat == NULL)
 625     {
 626         snew(mat, matrix->nx);
 627         for (i = 0; i < matrix->nx; i++)
 628         {
 629             snew(mat[i], matrix->ny);
 630         }
 631     }
 632     for (i = 0; i < matrix->nx; i++)
 633     {
 634         for (j = 0; j < matrix->ny; j++)
 635         {
 636             mat[i][j] = rmap[matrix->matrix[i][j]];
 637         }
 638     }
 639
 640     sfree(rmap);
 641
 642     fprintf(stderr, "Converted a %dx%d matrix with %d levels to reals\n",
 643             matrix->nx, matrix->ny, nmap);
 644
 645     return mat;
 646 }
 647
 648 void write_xpm_header(FILE *out,
 649                       const char *title, const char *legend,
 650                       const char *label_x, const char *label_y,
 651                       gmx_bool bDiscrete)
 652 {
 653     fprintf(out,  "/* XPM */\n");
 654     try
 655     {
 656         gmx::BinaryInformationSettings settings;
 657         settings.generatedByHeader(true);
 658         settings.linePrefix("/* ");
 659         settings.lineSuffix(" */");
 660         gmx::printBinaryInformation(out, gmx::ProgramInfo::getInstance(),
 661                                     settings);
 662     }
 663     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 664     fprintf(out,  "/* This file can be converted to EPS by the GROMACS program xpm2ps */\n");
 665     fprintf(out,  "/* title:   \"%s\" */\n", title);
 666     fprintf(out,  "/* legend:  \"%s\" */\n", legend);
 667     fprintf(out,  "/* x-label: \"%s\" */\n", label_x);
 668     fprintf(out,  "/* y-label: \"%s\" */\n", label_y);
 669     if (bDiscrete)
 670     {
 671         fprintf(out, "/* type:    \"Discrete\" */\n");
 672     }
 673     else
 674     {
 675         fprintf(out, "/* type:    \"Continuous\" */\n");
 676     }
 677 }
 678
 679 static int calc_nmid(int nlevels, real lo, real mid, real hi)
 680 {
 681     /* Take care that we have at least 1 entry in the mid to hi range
 682      */
 683     return std::min(std::max(0, static_cast<int>(((mid-lo)/(hi-lo))*(nlevels-1))),
 684                     nlevels-1);
 685 }
 686
 687 void write_xpm_map3(FILE *out, int n_x, int n_y, int *nlevels,
 688                     real lo, real mid, real hi,
 689                     t_rgb rlo, t_rgb rmid, t_rgb rhi)
 690 {
 691     int    i, nmid;
 692     real   r, g, b, clev_lo, clev_hi;
 693
 694     if (*nlevels > NMAP*NMAP)
 695     {
 696         fprintf(stderr, "Warning, too many levels (%d) in matrix, using %d only\n",
 697                 *nlevels, (int)(NMAP*NMAP));
 698         *nlevels = NMAP*NMAP;
 699     }
 700     else if (*nlevels < 2)
 701     {
 702         fprintf(stderr, "Warning, too few levels (%d) in matrix, using 2 instead\n",
 703                 *nlevels);
 704         *nlevels = 2;
 705     }
 706     if (!((mid >= lo) && (mid < hi)))
 707     {
 708         gmx_fatal(FARGS, "Lo: %f, Mid: %f, Hi: %f\n", lo, mid, hi);
 709     }
 710
 711     fprintf(out, "static char *gromacs_xpm[] = {\n");
 712     fprintf(out, "\"%d %d   %d %d\",\n",
 713             n_x, n_y, *nlevels, (*nlevels <= NMAP) ? 1 : 2);
 714
 715     nmid    = calc_nmid(*nlevels, lo, mid, hi);
 716     clev_lo = nmid;
 717     clev_hi = (*nlevels - 1 - nmid);
 718     for (i = 0; (i < nmid); i++)
 719     {
 720         r   = rlo.r+(i*(rmid.r-rlo.r)/clev_lo);
 721         g   = rlo.g+(i*(rmid.g-rlo.g)/clev_lo);
 722         b   = rlo.b+(i*(rmid.b-rlo.b)/clev_lo);
 723         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 724                 mapper[i % NMAP],
 725                 (*nlevels <= NMAP) ? ' ' : mapper[i/NMAP],
 726                 (unsigned int)round(255*r),
 727                 (unsigned int)round(255*g),
 728                 (unsigned int)round(255*b),
 729                 ((nmid - i)*lo + i*mid)/clev_lo);
 730     }
 731     for (i = 0; (i < (*nlevels-nmid)); i++)
 732     {
 733         r   = rmid.r+(i*(rhi.r-rmid.r)/clev_hi);
 734         g   = rmid.g+(i*(rhi.g-rmid.g)/clev_hi);
 735         b   = rmid.b+(i*(rhi.b-rmid.b)/clev_hi);
 736         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 737                 mapper[(i+nmid) % NMAP],
 738                 (*nlevels <= NMAP) ? ' ' : mapper[(i+nmid)/NMAP],
 739                 (unsigned int)round(255*r),
 740                 (unsigned int)round(255*g),
 741                 (unsigned int)round(255*b),
 742                 ((*nlevels - 1 - nmid - i)*mid + i*hi)/clev_hi);
 743     }
 744 }
 745
 746 static void pr_simple_cmap(FILE *out, real lo, real hi, int nlevel, t_rgb rlo,
 747                            t_rgb rhi, int i0)
 748 {
 749     int    i;
 750     real   r, g, b, fac;
 751
 752     for (i = 0; (i < nlevel); i++)
 753     {
 754         fac = (i+1.0)/(nlevel);
 755         r   = rlo.r+fac*(rhi.r-rlo.r);
 756         g   = rlo.g+fac*(rhi.g-rlo.g);
 757         b   = rlo.b+fac*(rhi.b-rlo.b);
 758         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 759                 mapper[(i+i0) % NMAP],
 760                 (nlevel <= NMAP) ? ' ' : mapper[(i+i0)/NMAP],
 761                 (unsigned int)round(255*r),
 762                 (unsigned int)round(255*g),
 763                 (unsigned int)round(255*b),
 764                 lo+fac*(hi-lo));
 765     }
 766 }
 767
 768 static void pr_discrete_cmap(FILE *out, int *nlevel, int i0)
 769 {
 770     t_rgb  rgbd[16] = {
 771         { 1.0, 1.0, 1.0 }, /* white */
 772         { 1.0, 0.0, 0.0 }, /* red */
 773         { 1.0, 1.0, 0.0 }, /* yellow */
 774         { 0.0, 0.0, 1.0 }, /* blue */
 775         { 0.0, 1.0, 0.0 }, /* green */
 776         { 1.0, 0.0, 1.0 }, /* purple */
 777         { 1.0, 0.4, 0.0 }, /* orange */
 778         { 0.0, 1.0, 1.0 }, /* cyan */
 779         { 1.0, 0.4, 0.4 }, /* pink */
 780         { 1.0, 1.0, 0.0 }, /* yellow */
 781         { 0.4, 0.4, 1.0 }, /* lightblue */
 782         { 0.4, 1.0, 0.4 }, /* lightgreen */
 783         { 1.0, 0.4, 1.0 }, /* lightpurple */
 784         { 1.0, 0.7, 0.4 }, /* lightorange */
 785         { 0.4, 1.0, 1.0 }, /* lightcyan */
 786         { 0.0, 0.0, 0.0 }  /* black */
 787     };
 788
 789     int    i, n;
 790
 791     *nlevel = std::min(16, *nlevel);
 792     n       = *nlevel;
 793     for (i = 0; (i < n); i++)
 794     {
 795         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%3d\" */,\n",
 796                 mapper[(i+i0) % NMAP],
 797                 (n <= NMAP) ? ' ' : mapper[(i+i0)/NMAP],
 798                 (unsigned int)round(255*rgbd[i].r),
 799                 (unsigned int)round(255*rgbd[i].g),
 800                 (unsigned int)round(255*rgbd[i].b),
 801                 i);
 802     }
 803 }
 804
 805
 806
 807 void write_xpm_map_split(FILE *out, int n_x, int n_y,
 808                          int *nlevel_top, real lo_top, real hi_top,
 809                          t_rgb rlo_top, t_rgb rhi_top,
 810                          gmx_bool bDiscreteColor,
 811                          int *nlevel_bot, real lo_bot, real hi_bot,
 812                          t_rgb rlo_bot, t_rgb rhi_bot)
 813 {
 814     int    ntot;
 815
 816     ntot = *nlevel_top + *nlevel_bot;
 817     if (ntot > NMAP)
 818     {
 819         gmx_fatal(FARGS, "Warning, too many levels (%d) in matrix", ntot);
 820     }
 821
 822     fprintf(out, "static char *gromacs_xpm[] = {\n");
 823     fprintf(out, "\"%d %d   %d %d\",\n", n_x, n_y, ntot, 1);
 824
 825     if (bDiscreteColor)
 826     {
 827         pr_discrete_cmap(out, nlevel_bot, 0);
 828     }
 829     else
 830     {
 831         pr_simple_cmap(out, lo_bot, hi_bot, *nlevel_bot, rlo_bot, rhi_bot, 0);
 832     }
 833
 834     pr_simple_cmap(out, lo_top, hi_top, *nlevel_top, rlo_top, rhi_top, *nlevel_bot);
 835 }
 836
 837
 838 void write_xpm_map(FILE *out, int n_x, int n_y, int *nlevels, real lo, real hi,
 839                    t_rgb rlo, t_rgb rhi)
 840 {
 841     int    i, nlo;
 842     real   invlevel, r, g, b;
 843
 844     if (*nlevels > NMAP*NMAP)
 845     {
 846         fprintf(stderr, "Warning, too many levels (%d) in matrix, using %d only\n",
 847                 *nlevels, (int)(NMAP*NMAP));
 848         *nlevels = NMAP*NMAP;
 849     }
 850     else if (*nlevels < 2)
 851     {
 852         fprintf(stderr, "Warning, too few levels (%d) in matrix, using 2 instead\n", *nlevels);
 853         *nlevels = 2;
 854     }
 855
 856     fprintf(out, "static char *gromacs_xpm[] = {\n");
 857     fprintf(out, "\"%d %d   %d %d\",\n",
 858             n_x, n_y, *nlevels, (*nlevels <= NMAP) ? 1 : 2);
 859
 860     invlevel = 1.0/(*nlevels-1);
 861     for (i = 0; (i < *nlevels); i++)
 862     {
 863         nlo = *nlevels-1-i;
 864         r   = (nlo*rlo.r+i*rhi.r)*invlevel;
 865         g   = (nlo*rlo.g+i*rhi.g)*invlevel;
 866         b   = (nlo*rlo.b+i*rhi.b)*invlevel;
 867         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 868                 mapper[i % NMAP], (*nlevels <= NMAP) ? ' ' : mapper[i/NMAP],
 869                 (unsigned int)round(255*r),
 870                 (unsigned int)round(255*g),
 871                 (unsigned int)round(255*b),
 872                 (nlo*lo+i*hi)*invlevel);
 873     }
 874 }
 875
 876 void write_xpm_axis(FILE *out, const char *axis, gmx_bool bSpatial, int n,
 877                     real *label)
 878 {
 879     int i;
 880
 881     if (label)
 882     {
 883         for (i = 0; i < (bSpatial ? n+1 : n); i++)
 884         {
 885             if (i % 80 == 0)
 886             {
 887                 if (i)
 888                 {
 889                     fprintf(out, "*/\n");
 890                 }
 891                 fprintf(out, "/* %s-axis:  ", axis);
 892             }
 893             fprintf(out, "%g ", label[i]);
 894         }
 895         fprintf(out, "*/\n");
 896     }
 897 }
 898
 899 void write_xpm_data(FILE *out, int n_x, int n_y, real **matrix,
 900                     real lo, real hi, int nlevels)
 901 {
 902     int  i, j, c;
 903     real invlevel;
 904
 905     invlevel = (nlevels-1)/(hi-lo);
 906     for (j = n_y-1; (j >= 0); j--)
 907     {
 908         if (j%(1+n_y/100) == 0)
 909         {
 910             fprintf(stderr, "%3d%%\b\b\b\b", (100*(n_y-j))/n_y);
 911         }
 912         fprintf(out, "\"");
 913         for (i = 0; (i < n_x); i++)
 914         {
 915             c = gmx_nint((matrix[i][j]-lo)*invlevel);
 916             if (c < 0)
 917             {
 918                 c = 0;
 919             }
 920             if (c >= nlevels)
 921             {
 922                 c = nlevels-1;
 923             }
 924             if (nlevels <= NMAP)
 925             {
 926                 fprintf(out, "%c", mapper[c]);
 927             }
 928             else
 929             {
 930                 fprintf(out, "%c%c", mapper[c % NMAP], mapper[c / NMAP]);
 931             }
 932         }
 933         if (j > 0)
 934         {
 935             fprintf(out, "\",\n");
 936         }
 937         else
 938         {
 939             fprintf(out, "\"\n");
 940         }
 941     }
 942 }
 943
 944 void write_xpm_data3(FILE *out, int n_x, int n_y, real **matrix,
 945                      real lo, real mid, real hi, int nlevels)
 946 {
 947     int  i, j, c = 0, nmid;
 948     real invlev_lo, invlev_hi;
 949
 950     nmid      = calc_nmid(nlevels, lo, mid, hi);
 951     invlev_hi = (nlevels-1-nmid)/(hi-mid);
 952     invlev_lo = (nmid)/(mid-lo);
 953
 954     for (j = n_y-1; (j >= 0); j--)
 955     {
 956         if (j%(1+n_y/100) == 0)
 957         {
 958             fprintf(stderr, "%3d%%\b\b\b\b", (100*(n_y-j))/n_y);
 959         }
 960         fprintf(out, "\"");
 961         for (i = 0; (i < n_x); i++)
 962         {
 963             if (matrix[i][j] >= mid)
 964             {
 965                 c = nmid+gmx_nint((matrix[i][j]-mid)*invlev_hi);
 966             }
 967             else if (matrix[i][j] >= lo)
 968             {
 969                 c = gmx_nint((matrix[i][j]-lo)*invlev_lo);
 970             }
 971             else
 972             {
 973                 c = 0;
 974             }
 975
 976             if (c < 0)
 977             {
 978                 c = 0;
 979             }
 980             if (c >= nlevels)
 981             {
 982                 c = nlevels-1;
 983             }
 984             if (nlevels <= NMAP)
 985             {
 986                 fprintf(out, "%c", mapper[c]);
 987             }
 988             else
 989             {
 990                 fprintf(out, "%c%c", mapper[c % NMAP], mapper[c / NMAP]);
 991             }
 992         }
 993         if (j > 0)
 994         {
 995             fprintf(out, "\",\n");
 996         }
 997         else
 998         {
 999             fprintf(out, "\"\n");
1000         }
1001     }
1002 }
1003
1004 void write_xpm_data_split(FILE *out, int n_x, int n_y, real **matrix,
1005                           real lo_top, real hi_top, int nlevel_top,
1006                           real lo_bot, real hi_bot, int nlevel_bot)
1007 {
1008     int  i, j, c;
1009     real invlev_top, invlev_bot;
1010
1011     invlev_top = (nlevel_top-1)/(hi_top-lo_top);
1012     invlev_bot = (nlevel_bot-1)/(hi_bot-lo_bot);
1013
1014     for (j = n_y-1; (j >= 0); j--)
1015     {
1016         if (j % (1+n_y/100) == 0)
1017         {
1018             fprintf(stderr, "%3d%%\b\b\b\b", (100*(n_y-j))/n_y);
1019         }
1020         fprintf(out, "\"");
1021         for (i = 0; (i < n_x); i++)
1022         {
1023             if (i < j)
1024             {
1025                 c = nlevel_bot+round((matrix[i][j]-lo_top)*invlev_top);
1026                 if ((c < nlevel_bot) || (c >= nlevel_bot+nlevel_top))
1027                 {
1028                     gmx_fatal(FARGS, "Range checking i = %d, j = %d, c = %d, bot = %d, top = %d matrix[i,j] = %f", i, j, c, nlevel_bot, nlevel_top, matrix[i][j]);
1029                 }
1030             }
1031             else if (i > j)
1032             {
1033                 c = round((matrix[i][j]-lo_bot)*invlev_bot);
1034                 if ((c < 0) || (c >= nlevel_bot+nlevel_bot))
1035                 {
1036                     gmx_fatal(FARGS, "Range checking i = %d, j = %d, c = %d, bot = %d, top = %d matrix[i,j] = %f", i, j, c, nlevel_bot, nlevel_top, matrix[i][j]);
1037                 }
1038             }
1039             else
1040             {
1041                 c = nlevel_bot;
1042             }
1043
1044             fprintf(out, "%c", mapper[c]);
1045         }
1046         if (j > 0)
1047         {
1048             fprintf(out, "\",\n");
1049         }
1050         else
1051         {
1052             fprintf(out, "\"\n");
1053         }
1054     }
1055 }
1056
1057 void write_xpm_m(FILE *out, t_matrix m)
1058 {
1059     /* Writes a t_matrix struct to .xpm file */
1060
1061     int           i, j;
1062     gmx_bool      bOneChar;
1063     t_xpmelmt     c;
1064
1065     bOneChar = (m.map[0].code.c2 == 0);
1066     write_xpm_header(out, m.title, m.legend, m.label_x, m.label_y,
1067                      m.bDiscrete);
1068     fprintf(out, "static char *gromacs_xpm[] = {\n");
1069     fprintf(out, "\"%d %d   %d %d\",\n", m.nx, m.ny, m.nmap, bOneChar ? 1 : 2);
1070     for (i = 0; (i < m.nmap); i++)
1071     {
1072         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%s\" */,\n",
1073                 m.map[i].code.c1,
1074                 bOneChar ? ' ' : m.map[i].code.c2,
1075                 (unsigned int)round(m.map[i].rgb.r*255),
1076                 (unsigned int)round(m.map[i].rgb.g*255),
1077                 (unsigned int)round(m.map[i].rgb.b*255), m.map[i].desc);
1078     }
1079     write_xpm_axis(out, "x", m.flags & MAT_SPATIAL_X, m.nx, m.axis_x);
1080     write_xpm_axis(out, "y", m.flags & MAT_SPATIAL_Y, m.ny, m.axis_y);
1081     for (j = m.ny-1; (j >= 0); j--)
1082     {
1083         if (j%(1+m.ny/100) == 0)
1084         {
1085             fprintf(stderr, "%3d%%\b\b\b\b", (100*(m.ny-j))/m.ny);
1086         }
1087         fprintf(out, "\"");
1088         if (bOneChar)
1089         {
1090             for (i = 0; (i < m.nx); i++)
1091             {
1092                 fprintf(out, "%c", m.map[m.matrix[i][j]].code.c1);
1093             }
1094         }
1095         else
1096         {
1097             for (i = 0; (i < m.nx); i++)
1098             {
1099                 c = m.map[m.matrix[i][j]].code;
1100                 fprintf(out, "%c%c", c.c1, c.c2);
1101             }
1102         }
1103         if (j > 0)
1104         {
1105             fprintf(out, "\",\n");
1106         }
1107         else
1108         {
1109             fprintf(out, "\"\n");
1110         }
1111     }
1112 }
1113
1114 void write_xpm3(FILE *out, unsigned int flags,
1115                 const char *title, const char *legend,
1116                 const char *label_x, const char *label_y,
1117                 int n_x, int n_y, real axis_x[], real axis_y[],
1118                 real *matrix[], real lo, real mid, real hi,
1119                 t_rgb rlo, t_rgb rmid, t_rgb rhi, int *nlevels)
1120 {
1121     /* See write_xpm.
1122      * Writes a colormap varying as rlo -> rmid -> rhi.
1123      */
1124
1125     if (hi <= lo)
1126     {
1127         gmx_fatal(FARGS, "hi (%g) <= lo (%g)", hi, lo);
1128     }
1129
1130     write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1131     write_xpm_map3(out, n_x, n_y, nlevels, lo, mid, hi, rlo, rmid, rhi);
1132     write_xpm_axis(out, "x", flags & MAT_SPATIAL_X, n_x, axis_x);
1133     write_xpm_axis(out, "y", flags & MAT_SPATIAL_Y, n_y, axis_y);
1134     write_xpm_data3(out, n_x, n_y, matrix, lo, mid, hi, *nlevels);
1135 }
1136
1137 void write_xpm_split(FILE *out, unsigned int flags,
1138                      const char *title, const char *legend,
1139                      const char *label_x, const char *label_y,
1140                      int n_x, int n_y, real axis_x[], real axis_y[],
1141                      real *matrix[],
1142                      real lo_top, real hi_top, int *nlevel_top,
1143                      t_rgb rlo_top, t_rgb rhi_top,
1144                      real lo_bot, real hi_bot, int *nlevel_bot,
1145                      gmx_bool bDiscreteColor,
1146                      t_rgb rlo_bot, t_rgb rhi_bot)
1147 {
1148     /* See write_xpm.
1149      * Writes a colormap varying as rlo -> rmid -> rhi.
1150      */
1151
1152     if (hi_top <= lo_top)
1153     {
1154         gmx_fatal(FARGS, "hi_top (%g) <= lo_top (%g)", hi_top, lo_top);
1155     }
1156     if (hi_bot <= lo_bot)
1157     {
1158         gmx_fatal(FARGS, "hi_bot (%g) <= lo_bot (%g)", hi_bot, lo_bot);
1159     }
1160     if (bDiscreteColor && (*nlevel_bot >= 16))
1161     {
1162         gmx_impl("Can not plot more than 16 discrete colors");
1163     }
1164
1165     write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1166     write_xpm_map_split(out, n_x, n_y, nlevel_top, lo_top, hi_top, rlo_top, rhi_top,
1167                         bDiscreteColor, nlevel_bot, lo_bot, hi_bot, rlo_bot, rhi_bot);
1168     write_xpm_axis(out, "x", flags & MAT_SPATIAL_X, n_x, axis_x);
1169     write_xpm_axis(out, "y", flags & MAT_SPATIAL_Y, n_y, axis_y);
1170     write_xpm_data_split(out, n_x, n_y, matrix, lo_top, hi_top, *nlevel_top,
1171                          lo_bot, hi_bot, *nlevel_bot);
1172 }
1173
1174 void write_xpm(FILE *out, unsigned int flags,
1175                const char *title, const char *legend,
1176                const char *label_x, const char *label_y,
1177                int n_x, int n_y, real axis_x[], real axis_y[],
1178                real *matrix[], real lo, real hi,
1179                t_rgb rlo, t_rgb rhi, int *nlevels)
1180 {
1181     /* out        xpm file
1182      * title      matrix title
1183      * legend     label for the continuous legend
1184      * label_x    label for the x-axis
1185      * label_y    label for the y-axis
1186      * n_x, n_y   size of the matrix
1187      * axis_x[]   the x-ticklabels
1188      * axis_y[]   the y-ticklables
1189      * *matrix[]  element x,y is matrix[x][y]
1190      * lo         output lower than lo is set to lo
1191      * hi         output higher than hi is set to hi
1192      * rlo        rgb value for level lo
1193      * rhi        rgb value for level hi
1194      * nlevels    number of color levels for the output
1195      */
1196
1197     if (hi <= lo)
1198     {
1199         gmx_fatal(FARGS, "hi (%f) <= lo (%f)", hi, lo);
1200     }
1201
1202     write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1203     write_xpm_map(out, n_x, n_y, nlevels, lo, hi, rlo, rhi);
1204     write_xpm_axis(out, "x", flags & MAT_SPATIAL_X, n_x, axis_x);
1205     write_xpm_axis(out, "y", flags & MAT_SPATIAL_Y, n_y, axis_y);
1206     write_xpm_data(out, n_x, n_y, matrix, lo, hi, *nlevels);
1207 }