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