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, 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 <cstdio>
  43 #include <cstring>
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/fileio/gmxfio.h"
  48 #include "gromacs/legacyheaders/copyrite.h"
  49 #include "gromacs/math/utilities.h"
  50 #include "gromacs/utility/cstringutil.h"
  51 #include "gromacs/utility/exceptions.h"
  52 #include "gromacs/utility/fatalerror.h"
  53 #include "gromacs/utility/futil.h"
  54 #include "gromacs/utility/programcontext.h"
  55 #include "gromacs/utility/smalloc.h"
  56
  57 #define round(a) (int)(a+0.5)
  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 = NULL;
 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) == NULL)
 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) == NULL)
 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 != NULL)
 231     {
 232         while (((*line)[0] != '\"' ) && ( (*line)[0] != '\0' ))
 233         {
 234             (*line)++;
 235         }
 236
 237         if ((*line)[0] != '\"')
 238         {
 239             return NULL;
 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 = NULL;
 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 = NULL, *line = NULL, *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     = NULL;
 291     mm->axis_x     = NULL;
 292     mm->axis_y     = NULL;
 293     mm->bDiscrete  = FALSE;
 294
 295     llmax = STRLEN;
 296
 297     while ((NULL != 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 && (NULL != 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) && (NULL != 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 == NULL)
 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 == NULL)
 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] != '\"') && (NULL != 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) && (NULL != 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 = NULL;
 556     int   nmat;
 557     int   llalloc = 0;
 558
 559     in = gmx_fio_fopen(fnm, "r");
 560
 561     nmat = 0;
 562     while (NULL != 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 == NULL) || (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 NULL;
 603         }
 604         rmap[i] = tmp;
 605     }
 606
 607     if (out == NULL)
 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     try
 638     {
 639         gmx::BinaryInformationSettings settings;
 640         settings.generatedByHeader(true);
 641         settings.linePrefix("/* ");
 642         settings.lineSuffix(" */");
 643         gmx::printBinaryInformation(out, gmx::getProgramContext(),
 644                                     settings);
 645     }
 646     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 647     fprintf(out,  "/* This file can be converted to EPS by the GROMACS program xpm2ps */\n");
 648     fprintf(out,  "/* title:   \"%s\" */\n", title);
 649     fprintf(out,  "/* legend:  \"%s\" */\n", legend);
 650     fprintf(out,  "/* x-label: \"%s\" */\n", label_x);
 651     fprintf(out,  "/* y-label: \"%s\" */\n", label_y);
 652     if (bDiscrete)
 653     {
 654         fprintf(out, "/* type:    \"Discrete\" */\n");
 655     }
 656     else
 657     {
 658         fprintf(out, "/* type:    \"Continuous\" */\n");
 659     }
 660 }
 661
 662 static int calc_nmid(int nlevels, real lo, real mid, real hi)
 663 {
 664     /* Take care that we have at least 1 entry in the mid to hi range
 665      */
 666     return std::min(std::max(0, static_cast<int>(((mid-lo)/(hi-lo))*(nlevels-1))),
 667                     nlevels-1);
 668 }
 669
 670 static void write_xpm_map3(FILE *out, int n_x, int n_y, int *nlevels,
 671                            real lo, real mid, real hi,
 672                            t_rgb rlo, t_rgb rmid, t_rgb rhi)
 673 {
 674     int    i, nmid;
 675     real   r, g, b, clev_lo, clev_hi;
 676
 677     if (*nlevels > NMAP*NMAP)
 678     {
 679         fprintf(stderr, "Warning, too many levels (%d) in matrix, using %d only\n",
 680                 *nlevels, (int)(NMAP*NMAP));
 681         *nlevels = NMAP*NMAP;
 682     }
 683     else if (*nlevels < 2)
 684     {
 685         fprintf(stderr, "Warning, too few levels (%d) in matrix, using 2 instead\n",
 686                 *nlevels);
 687         *nlevels = 2;
 688     }
 689     if (!((mid >= lo) && (mid < hi)))
 690     {
 691         gmx_fatal(FARGS, "Lo: %f, Mid: %f, Hi: %f\n", lo, mid, hi);
 692     }
 693
 694     fprintf(out, "static char *gromacs_xpm[] = {\n");
 695     fprintf(out, "\"%d %d   %d %d\",\n",
 696             n_x, n_y, *nlevels, (*nlevels <= NMAP) ? 1 : 2);
 697
 698     nmid    = calc_nmid(*nlevels, lo, mid, hi);
 699     clev_lo = nmid;
 700     clev_hi = (*nlevels - 1 - nmid);
 701     for (i = 0; (i < nmid); i++)
 702     {
 703         r   = rlo.r+(i*(rmid.r-rlo.r)/clev_lo);
 704         g   = rlo.g+(i*(rmid.g-rlo.g)/clev_lo);
 705         b   = rlo.b+(i*(rmid.b-rlo.b)/clev_lo);
 706         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 707                 mapper[i % NMAP],
 708                 (*nlevels <= NMAP) ? ' ' : mapper[i/NMAP],
 709                 static_cast<unsigned int>(round(255*r)),
 710                 static_cast<unsigned int>(round(255*g)),
 711                 static_cast<unsigned int>(round(255*b)),
 712                 ((nmid - i)*lo + i*mid)/clev_lo);
 713     }
 714     for (i = 0; (i < (*nlevels-nmid)); i++)
 715     {
 716         r   = rmid.r+(i*(rhi.r-rmid.r)/clev_hi);
 717         g   = rmid.g+(i*(rhi.g-rmid.g)/clev_hi);
 718         b   = rmid.b+(i*(rhi.b-rmid.b)/clev_hi);
 719         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 720                 mapper[(i+nmid) % NMAP],
 721                 (*nlevels <= NMAP) ? ' ' : mapper[(i+nmid)/NMAP],
 722                 static_cast<unsigned int>(round(255*r)),
 723                 static_cast<unsigned int>(round(255*g)),
 724                 static_cast<unsigned int>(round(255*b)),
 725                 ((*nlevels - 1 - nmid - i)*mid + i*hi)/clev_hi);
 726     }
 727 }
 728
 729 static void pr_simple_cmap(FILE *out, real lo, real hi, int nlevel, t_rgb rlo,
 730                            t_rgb rhi, int i0)
 731 {
 732     int    i;
 733     real   r, g, b, fac;
 734
 735     for (i = 0; (i < nlevel); i++)
 736     {
 737         fac = (i+1.0)/(nlevel);
 738         r   = rlo.r+fac*(rhi.r-rlo.r);
 739         g   = rlo.g+fac*(rhi.g-rlo.g);
 740         b   = rlo.b+fac*(rhi.b-rlo.b);
 741         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 742                 mapper[(i+i0) % NMAP],
 743                 (nlevel <= NMAP) ? ' ' : mapper[(i+i0)/NMAP],
 744                 static_cast<unsigned int>(round(255*r)),
 745                 static_cast<unsigned int>(round(255*g)),
 746                 static_cast<unsigned int>(round(255*b)),
 747                 lo+fac*(hi-lo));
 748     }
 749 }
 750
 751 static void pr_discrete_cmap(FILE *out, int *nlevel, int i0)
 752 {
 753     t_rgb  rgbd[16] = {
 754         { 1.0, 1.0, 1.0 }, /* white */
 755         { 1.0, 0.0, 0.0 }, /* red */
 756         { 1.0, 1.0, 0.0 }, /* yellow */
 757         { 0.0, 0.0, 1.0 }, /* blue */
 758         { 0.0, 1.0, 0.0 }, /* green */
 759         { 1.0, 0.0, 1.0 }, /* purple */
 760         { 1.0, 0.4, 0.0 }, /* orange */
 761         { 0.0, 1.0, 1.0 }, /* cyan */
 762         { 1.0, 0.4, 0.4 }, /* pink */
 763         { 1.0, 1.0, 0.0 }, /* yellow */
 764         { 0.4, 0.4, 1.0 }, /* lightblue */
 765         { 0.4, 1.0, 0.4 }, /* lightgreen */
 766         { 1.0, 0.4, 1.0 }, /* lightpurple */
 767         { 1.0, 0.7, 0.4 }, /* lightorange */
 768         { 0.4, 1.0, 1.0 }, /* lightcyan */
 769         { 0.0, 0.0, 0.0 }  /* black */
 770     };
 771
 772     int    i, n;
 773
 774     *nlevel = std::min(16, *nlevel);
 775     n       = *nlevel;
 776     for (i = 0; (i < n); i++)
 777     {
 778         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%3d\" */,\n",
 779                 mapper[(i+i0) % NMAP],
 780                 (n <= NMAP) ? ' ' : mapper[(i+i0)/NMAP],
 781                 static_cast<unsigned int>(round(255*rgbd[i].r)),
 782                 static_cast<unsigned int>(round(255*rgbd[i].g)),
 783                 static_cast<unsigned int>(round(255*rgbd[i].b)),
 784                 i);
 785     }
 786 }
 787
 788
 789
 790 static void write_xpm_map_split(FILE *out, int n_x, int n_y,
 791                                 int *nlevel_top, real lo_top, real hi_top,
 792                                 t_rgb rlo_top, t_rgb rhi_top,
 793                                 gmx_bool bDiscreteColor,
 794                                 int *nlevel_bot, real lo_bot, real hi_bot,
 795                                 t_rgb rlo_bot, t_rgb rhi_bot)
 796 {
 797     int    ntot;
 798
 799     ntot = *nlevel_top + *nlevel_bot;
 800     if (ntot > NMAP)
 801     {
 802         gmx_fatal(FARGS, "Warning, too many levels (%d) in matrix", ntot);
 803     }
 804
 805     fprintf(out, "static char *gromacs_xpm[] = {\n");
 806     fprintf(out, "\"%d %d   %d %d\",\n", n_x, n_y, ntot, 1);
 807
 808     if (bDiscreteColor)
 809     {
 810         pr_discrete_cmap(out, nlevel_bot, 0);
 811     }
 812     else
 813     {
 814         pr_simple_cmap(out, lo_bot, hi_bot, *nlevel_bot, rlo_bot, rhi_bot, 0);
 815     }
 816
 817     pr_simple_cmap(out, lo_top, hi_top, *nlevel_top, rlo_top, rhi_top, *nlevel_bot);
 818 }
 819
 820
 821 static void write_xpm_map(FILE *out, int n_x, int n_y, int *nlevels,
 822                           real lo, real hi, t_rgb rlo, t_rgb rhi)
 823 {
 824     int    i, nlo;
 825     real   invlevel, r, g, b;
 826
 827     if (*nlevels > NMAP*NMAP)
 828     {
 829         fprintf(stderr, "Warning, too many levels (%d) in matrix, using %d only\n",
 830                 *nlevels, (int)(NMAP*NMAP));
 831         *nlevels = NMAP*NMAP;
 832     }
 833     else if (*nlevels < 2)
 834     {
 835         fprintf(stderr, "Warning, too few levels (%d) in matrix, using 2 instead\n", *nlevels);
 836         *nlevels = 2;
 837     }
 838
 839     fprintf(out, "static char *gromacs_xpm[] = {\n");
 840     fprintf(out, "\"%d %d   %d %d\",\n",
 841             n_x, n_y, *nlevels, (*nlevels <= NMAP) ? 1 : 2);
 842
 843     invlevel = 1.0/(*nlevels-1);
 844     for (i = 0; (i < *nlevels); i++)
 845     {
 846         nlo = *nlevels-1-i;
 847         r   = (nlo*rlo.r+i*rhi.r)*invlevel;
 848         g   = (nlo*rlo.g+i*rhi.g)*invlevel;
 849         b   = (nlo*rlo.b+i*rhi.b)*invlevel;
 850         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 851                 mapper[i % NMAP], (*nlevels <= NMAP) ? ' ' : mapper[i/NMAP],
 852                 static_cast<unsigned int>(round(255*r)),
 853                 static_cast<unsigned int>(round(255*g)),
 854                 static_cast<unsigned int>(round(255*b)),
 855                 (nlo*lo+i*hi)*invlevel);
 856     }
 857 }
 858
 859 static void write_xpm_axis(FILE *out, const char *axis, gmx_bool bSpatial,
 860                            int n, real *label)
 861 {
 862     int i;
 863
 864     if (label)
 865     {
 866         for (i = 0; i < (bSpatial ? n+1 : n); i++)
 867         {
 868             if (i % 80 == 0)
 869             {
 870                 if (i)
 871                 {
 872                     fprintf(out, "*/\n");
 873                 }
 874                 fprintf(out, "/* %s-axis:  ", axis);
 875             }
 876             fprintf(out, "%g ", label[i]);
 877         }
 878         fprintf(out, "*/\n");
 879     }
 880 }
 881
 882 static void write_xpm_data(FILE *out, int n_x, int n_y, real **mat,
 883                            real lo, real hi, int nlevels)
 884 {
 885     int  i, j, c;
 886     real invlevel;
 887
 888     invlevel = (nlevels-1)/(hi-lo);
 889     for (j = n_y-1; (j >= 0); j--)
 890     {
 891         if (j%(1+n_y/100) == 0)
 892         {
 893             fprintf(stderr, "%3d%%\b\b\b\b", (100*(n_y-j))/n_y);
 894         }
 895         fprintf(out, "\"");
 896         for (i = 0; (i < n_x); i++)
 897         {
 898             c = gmx_nint((mat[i][j]-lo)*invlevel);
 899             if (c < 0)
 900             {
 901                 c = 0;
 902             }
 903             if (c >= nlevels)
 904             {
 905                 c = nlevels-1;
 906             }
 907             if (nlevels <= NMAP)
 908             {
 909                 fprintf(out, "%c", mapper[c]);
 910             }
 911             else
 912             {
 913                 fprintf(out, "%c%c", mapper[c % NMAP], mapper[c / NMAP]);
 914             }
 915         }
 916         if (j > 0)
 917         {
 918             fprintf(out, "\",\n");
 919         }
 920         else
 921         {
 922             fprintf(out, "\"\n");
 923         }
 924     }
 925 }
 926
 927 static void write_xpm_data3(FILE *out, int n_x, int n_y, real **mat,
 928                             real lo, real mid, real hi, int nlevels)
 929 {
 930     int  i, j, c = 0, nmid;
 931     real invlev_lo, invlev_hi;
 932
 933     nmid      = calc_nmid(nlevels, lo, mid, hi);
 934     invlev_hi = (nlevels-1-nmid)/(hi-mid);
 935     invlev_lo = (nmid)/(mid-lo);
 936
 937     for (j = n_y-1; (j >= 0); j--)
 938     {
 939         if (j%(1+n_y/100) == 0)
 940         {
 941             fprintf(stderr, "%3d%%\b\b\b\b", (100*(n_y-j))/n_y);
 942         }
 943         fprintf(out, "\"");
 944         for (i = 0; (i < n_x); i++)
 945         {
 946             if (mat[i][j] >= mid)
 947             {
 948                 c = nmid+gmx_nint((mat[i][j]-mid)*invlev_hi);
 949             }
 950             else if (mat[i][j] >= lo)
 951             {
 952                 c = gmx_nint((mat[i][j]-lo)*invlev_lo);
 953             }
 954             else
 955             {
 956                 c = 0;
 957             }
 958
 959             if (c < 0)
 960             {
 961                 c = 0;
 962             }
 963             if (c >= nlevels)
 964             {
 965                 c = nlevels-1;
 966             }
 967             if (nlevels <= NMAP)
 968             {
 969                 fprintf(out, "%c", mapper[c]);
 970             }
 971             else
 972             {
 973                 fprintf(out, "%c%c", mapper[c % NMAP], mapper[c / NMAP]);
 974             }
 975         }
 976         if (j > 0)
 977         {
 978             fprintf(out, "\",\n");
 979         }
 980         else
 981         {
 982             fprintf(out, "\"\n");
 983         }
 984     }
 985 }
 986
 987 static void write_xpm_data_split(FILE *out, int n_x, int n_y, real **mat,
 988                                  real lo_top, real hi_top, int nlevel_top,
 989                                  real lo_bot, real hi_bot, int nlevel_bot)
 990 {
 991     int  i, j, c;
 992     real invlev_top, invlev_bot;
 993
 994     invlev_top = (nlevel_top-1)/(hi_top-lo_top);
 995     invlev_bot = (nlevel_bot-1)/(hi_bot-lo_bot);
 996
 997     for (j = n_y-1; (j >= 0); j--)
 998     {
 999         if (j % (1+n_y/100) == 0)
1000         {
1001             fprintf(stderr, "%3d%%\b\b\b\b", (100*(n_y-j))/n_y);
1002         }
1003         fprintf(out, "\"");
1004         for (i = 0; (i < n_x); i++)
1005         {
1006             if (i < j)
1007             {
1008                 c = nlevel_bot+round((mat[i][j]-lo_top)*invlev_top);
1009                 if ((c < nlevel_bot) || (c >= nlevel_bot+nlevel_top))
1010                 {
1011                     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]);
1012                 }
1013             }
1014             else if (i > j)
1015             {
1016                 c = round((mat[i][j]-lo_bot)*invlev_bot);
1017                 if ((c < 0) || (c >= nlevel_bot+nlevel_bot))
1018                 {
1019                     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]);
1020                 }
1021             }
1022             else
1023             {
1024                 c = nlevel_bot;
1025             }
1026
1027             fprintf(out, "%c", mapper[c]);
1028         }
1029         if (j > 0)
1030         {
1031             fprintf(out, "\",\n");
1032         }
1033         else
1034         {
1035             fprintf(out, "\"\n");
1036         }
1037     }
1038 }
1039
1040 void write_xpm_m(FILE *out, t_matrix m)
1041 {
1042     /* Writes a t_matrix struct to .xpm file */
1043
1044     int           i, j;
1045     gmx_bool      bOneChar;
1046     t_xpmelmt     c;
1047
1048     bOneChar = (m.map[0].code.c2 == 0);
1049     write_xpm_header(out, m.title, m.legend, m.label_x, m.label_y,
1050                      m.bDiscrete);
1051     fprintf(out, "static char *gromacs_xpm[] = {\n");
1052     fprintf(out, "\"%d %d   %d %d\",\n", m.nx, m.ny, m.nmap, bOneChar ? 1 : 2);
1053     for (i = 0; (i < m.nmap); i++)
1054     {
1055         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%s\" */,\n",
1056                 m.map[i].code.c1,
1057                 bOneChar ? ' ' : m.map[i].code.c2,
1058                 static_cast<unsigned int>(round(m.map[i].rgb.r*255)),
1059                 static_cast<unsigned int>(round(m.map[i].rgb.g*255)),
1060                 static_cast<unsigned int>(round(m.map[i].rgb.b*255)), m.map[i].desc);
1061     }
1062     write_xpm_axis(out, "x", m.flags & MAT_SPATIAL_X, m.nx, m.axis_x);
1063     write_xpm_axis(out, "y", m.flags & MAT_SPATIAL_Y, m.ny, m.axis_y);
1064     for (j = m.ny-1; (j >= 0); j--)
1065     {
1066         if (j%(1+m.ny/100) == 0)
1067         {
1068             fprintf(stderr, "%3d%%\b\b\b\b", (100*(m.ny-j))/m.ny);
1069         }
1070         fprintf(out, "\"");
1071         if (bOneChar)
1072         {
1073             for (i = 0; (i < m.nx); i++)
1074             {
1075                 fprintf(out, "%c", m.map[m.matrix[i][j]].code.c1);
1076             }
1077         }
1078         else
1079         {
1080             for (i = 0; (i < m.nx); i++)
1081             {
1082                 c = m.map[m.matrix[i][j]].code;
1083                 fprintf(out, "%c%c", c.c1, c.c2);
1084             }
1085         }
1086         if (j > 0)
1087         {
1088             fprintf(out, "\",\n");
1089         }
1090         else
1091         {
1092             fprintf(out, "\"\n");
1093         }
1094     }
1095 }
1096
1097 void write_xpm3(FILE *out, unsigned int flags,
1098                 const char *title, const char *legend,
1099                 const char *label_x, const char *label_y,
1100                 int n_x, int n_y, real axis_x[], real axis_y[],
1101                 real *mat[], real lo, real mid, real hi,
1102                 t_rgb rlo, t_rgb rmid, t_rgb rhi, int *nlevels)
1103 {
1104     /* See write_xpm.
1105      * Writes a colormap varying as rlo -> rmid -> rhi.
1106      */
1107
1108     if (hi <= lo)
1109     {
1110         gmx_fatal(FARGS, "hi (%g) <= lo (%g)", hi, lo);
1111     }
1112
1113     write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1114     write_xpm_map3(out, n_x, n_y, nlevels, lo, mid, hi, rlo, rmid, rhi);
1115     write_xpm_axis(out, "x", flags & MAT_SPATIAL_X, n_x, axis_x);
1116     write_xpm_axis(out, "y", flags & MAT_SPATIAL_Y, n_y, axis_y);
1117     write_xpm_data3(out, n_x, n_y, mat, lo, mid, hi, *nlevels);
1118 }
1119
1120 void write_xpm_split(FILE *out, unsigned int flags,
1121                      const char *title, const char *legend,
1122                      const char *label_x, const char *label_y,
1123                      int n_x, int n_y, real axis_x[], real axis_y[],
1124                      real *mat[],
1125                      real lo_top, real hi_top, int *nlevel_top,
1126                      t_rgb rlo_top, t_rgb rhi_top,
1127                      real lo_bot, real hi_bot, int *nlevel_bot,
1128                      gmx_bool bDiscreteColor,
1129                      t_rgb rlo_bot, t_rgb rhi_bot)
1130 {
1131     /* See write_xpm.
1132      * Writes a colormap varying as rlo -> rmid -> rhi.
1133      */
1134
1135     if (hi_top <= lo_top)
1136     {
1137         gmx_fatal(FARGS, "hi_top (%g) <= lo_top (%g)", hi_top, lo_top);
1138     }
1139     if (hi_bot <= lo_bot)
1140     {
1141         gmx_fatal(FARGS, "hi_bot (%g) <= lo_bot (%g)", hi_bot, lo_bot);
1142     }
1143     if (bDiscreteColor && (*nlevel_bot >= 16))
1144     {
1145         gmx_impl("Can not plot more than 16 discrete colors");
1146     }
1147
1148     write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1149     write_xpm_map_split(out, n_x, n_y, nlevel_top, lo_top, hi_top, rlo_top, rhi_top,
1150                         bDiscreteColor, nlevel_bot, lo_bot, hi_bot, rlo_bot, rhi_bot);
1151     write_xpm_axis(out, "x", flags & MAT_SPATIAL_X, n_x, axis_x);
1152     write_xpm_axis(out, "y", flags & MAT_SPATIAL_Y, n_y, axis_y);
1153     write_xpm_data_split(out, n_x, n_y, mat, lo_top, hi_top, *nlevel_top,
1154                          lo_bot, hi_bot, *nlevel_bot);
1155 }
1156
1157 void write_xpm(FILE *out, unsigned int flags,
1158                const char *title, const char *legend,
1159                const char *label_x, const char *label_y,
1160                int n_x, int n_y, real axis_x[], real axis_y[],
1161                real *mat[], real lo, real hi,
1162                t_rgb rlo, t_rgb rhi, int *nlevels)
1163 {
1164     /* out        xpm file
1165      * title      matrix title
1166      * legend     label for the continuous legend
1167      * label_x    label for the x-axis
1168      * label_y    label for the y-axis
1169      * n_x, n_y   size of the matrix
1170      * axis_x[]   the x-ticklabels
1171      * axis_y[]   the y-ticklables
1172      * *matrix[]  element x,y is matrix[x][y]
1173      * lo         output lower than lo is set to lo
1174      * hi         output higher than hi is set to hi
1175      * rlo        rgb value for level lo
1176      * rhi        rgb value for level hi
1177      * nlevels    number of color levels for the output
1178      */
1179
1180     if (hi <= lo)
1181     {
1182         gmx_fatal(FARGS, "hi (%f) <= lo (%f)", hi, lo);
1183     }
1184
1185     write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1186     write_xpm_map(out, n_x, n_y, nlevels, lo, hi, rlo, rhi);
1187     write_xpm_axis(out, "x", flags & MAT_SPATIAL_X, n_x, axis_x);
1188     write_xpm_axis(out, "y", flags & MAT_SPATIAL_Y, n_y, axis_y);
1189     write_xpm_data(out, n_x, n_y, mat, lo, hi, *nlevels);
1190 }