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