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