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