src/gromacs/fileio/xvgr.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,2016,2017,2018,2019, 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 "xvgr.h"
  40
  41 #include <cassert>
  42 #include <cctype>
  43 #include <cstring>
  44
  45 #include <string>
  46
  47 #include "gromacs/fileio/gmxfio.h"
  48 #include "gromacs/fileio/oenv.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/utility/binaryinformation.h"
  51 #include "gromacs/utility/coolstuff.h"
  52 #include "gromacs/utility/cstringutil.h"
  53 #include "gromacs/utility/exceptions.h"
  54 #include "gromacs/utility/fatalerror.h"
  55 #include "gromacs/utility/futil.h"
  56 #include "gromacs/utility/smalloc.h"
  57 #include "gromacs/utility/sysinfo.h"
  58
  59 gmx_bool output_env_get_print_xvgr_codes(const gmx_output_env_t* oenv)
  60 {
  61     int xvg_format;
  62
  63     xvg_format = output_env_get_xvg_format(oenv);
  64
  65     return (xvg_format == exvgXMGRACE || xvg_format == exvgXMGR);
  66 }
  67
  68 static char* xvgrstr(const std::string& gmx, const gmx_output_env_t* oenv, char* buf, int buflen)
  69 {
  70     /* Supported greek letter names and corresponding xmgrace/xmgr symbols */
  71     const char* sym[]  = { "beta",  "chi", "delta", "eta", "lambda", "mu",
  72                           "omega", "phi", "psi",   "rho", "theta",  nullptr };
  73     const char  symc[] = { 'b', 'c', 'd', 'h', 'l', 'm', 'w', 'f', 'y', 'r', 'q', '\0' };
  74     int         xvgf;
  75     gmx_bool    bXVGR;
  76     int         g, b, i;
  77     char        c;
  78
  79     xvgf  = output_env_get_xvg_format(oenv);
  80     bXVGR = (xvgf == exvgXMGRACE || xvgf == exvgXMGR);
  81
  82     g = 0;
  83     b = 0;
  84     while (gmx[g] != '\0')
  85     {
  86         /* Check with the largest string we have ("lambda"), add one for \0 */
  87         if (b + 6 + 1 >= buflen)
  88         {
  89             gmx_fatal(FARGS,
  90                       "Output buffer length in xvgstr (%d) too small to process xvg input string "
  91                       "'%s'",
  92                       buflen, gmx.c_str());
  93         }
  94         if (gmx[g] == '\\')
  95         {
  96             g++;
  97             if (gmx[g] == 's')
  98             {
  99                 /* Subscript */
 100                 if (bXVGR)
 101                 {
 102                     buf[b++] = '\\';
 103                     buf[b++] = 's';
 104                 }
 105                 else
 106                 {
 107                     buf[b++] = '_';
 108                 }
 109                 g++;
 110             }
 111             else if (gmx[g] == 'S')
 112             {
 113                 /* Superscript */
 114                 if (bXVGR)
 115                 {
 116                     buf[b++] = '\\';
 117                     buf[b++] = 'S';
 118                 }
 119                 else
 120                 {
 121                     buf[b++] = '^';
 122                 }
 123                 g++;
 124             }
 125             else if (gmx[g] == 'N')
 126             {
 127                 /* End sub/superscript */
 128                 if (bXVGR)
 129                 {
 130                     buf[b++] = '\\';
 131                     buf[b++] = 'N';
 132                 }
 133                 else
 134                 {
 135                     if (gmx[g + 1] != ' ')
 136                     {
 137                         buf[b++] = ' ';
 138                     }
 139                 }
 140                 g++;
 141             }
 142             else if (gmx[g] == '4')
 143             {
 144                 /* Backward compatibility for xmgr normal font "\4" */
 145                 switch (xvgf)
 146                 {
 147                     case exvgXMGRACE: sprintf(buf + b, "%s", "\\f{}"); break;
 148                     case exvgXMGR: sprintf(buf + b, "%s", "\\4"); break;
 149                     default: buf[b] = '\0'; break;
 150                 }
 151                 g++;
 152                 b = strlen(buf);
 153             }
 154             else if (gmx[g] == '8')
 155             {
 156                 /* Backward compatibility for xmgr symbol font "\8" */
 157                 switch (xvgf)
 158                 {
 159                     case exvgXMGRACE: sprintf(buf + b, "%s", "\\x"); break;
 160                     case exvgXMGR: sprintf(buf + b, "%s", "\\8"); break;
 161                     default: buf[b] = '\0'; break;
 162                 }
 163                 g++;
 164                 b = std::strlen(buf);
 165             }
 166             else
 167             {
 168                 /* Check for special symbol */
 169                 i = 0;
 170                 while (sym[i] != nullptr && gmx_strncasecmp(sym[i], &gmx[g], std::strlen(sym[i])) != 0)
 171                 {
 172                     i++;
 173                 }
 174                 if (sym[i] != nullptr)
 175                 {
 176                     c = symc[i];
 177                     if (std::isupper(gmx[g]))
 178                     {
 179                         c = std::toupper(c);
 180                     }
 181                     switch (xvgf)
 182                     {
 183                         case exvgXMGRACE: sprintf(buf + b, "%s%c%s", "\\x", c, "\\f{}"); break;
 184                         case exvgXMGR: sprintf(buf + b, "%s%c%s", "\\8", c, "\\4"); break;
 185                         default:
 186                             std::strncat(buf + b, &gmx[g], std::strlen(sym[i]));
 187                             b += std::strlen(sym[i]);
 188                             if (gmx[g + std::strlen(sym[i])] != ' ')
 189                             {
 190                                 buf[b++] = ' ';
 191                             }
 192                             buf[b] = '\0';
 193                             break;
 194                     }
 195                     g += std::strlen(sym[i]);
 196                     b = std::strlen(buf);
 197                 }
 198                 else
 199                 {
 200                     /* Unknown escape sequence, this should not happen.
 201                      * We do not generate a fatal error, since that might
 202                      * stop otherwise functioning code from working.
 203                      * Copy the backslash to the output and continue processing.
 204                      */
 205                     buf[b++] = '\\';
 206                 }
 207             }
 208         }
 209         else
 210         {
 211             buf[b++] = gmx[g++];
 212         }
 213     }
 214
 215     buf[b++] = '\0';
 216
 217     return buf;
 218 }
 219
 220 void xvgr_header(FILE*                   fp,
 221                  const char*             title,
 222                  const std::string&      xaxis,
 223                  const std::string&      yaxis,
 224                  int                     exvg_graph_type,
 225                  const gmx_output_env_t* oenv)
 226 {
 227     char buf[STRLEN];
 228
 229     if (output_env_get_print_xvgr_codes(oenv))
 230     {
 231         fprintf(fp, "# This file was created %s", gmx_format_current_time().c_str());
 232         try
 233         {
 234             gmx::BinaryInformationSettings settings;
 235             settings.generatedByHeader(true);
 236             settings.linePrefix("# ");
 237             gmx::printBinaryInformation(fp, output_env_get_program_context(oenv), settings);
 238         }
 239         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 240         fprintf(fp, "# %s is part of G R O M A C S:\n#\n", output_env_get_program_display_name(oenv));
 241         fprintf(fp, "# %s\n#\n", gmx::bromacs().c_str());
 242         fprintf(fp, "@    title \"%s\"\n", xvgrstr(title, oenv, buf, STRLEN));
 243         fprintf(fp, "@    xaxis  label \"%s\"\n", xvgrstr(xaxis, oenv, buf, STRLEN));
 244         fprintf(fp, "@    yaxis  label \"%s\"\n", xvgrstr(yaxis, oenv, buf, STRLEN));
 245         switch (exvg_graph_type)
 246         {
 247             case exvggtXNY:
 248                 if (output_env_get_xvg_format(oenv) == exvgXMGR)
 249                 {
 250                     fprintf(fp, "@TYPE nxy\n");
 251                 }
 252                 else
 253                 {
 254                     fprintf(fp, "@TYPE xy\n");
 255                 }
 256                 break;
 257             case exvggtXYDY: fprintf(fp, "@TYPE xydy\n"); break;
 258             case exvggtXYDYDY: fprintf(fp, "@TYPE xydydy\n"); break;
 259         }
 260     }
 261 }
 262
 263 FILE* xvgropen_type(const char*             fn,
 264                     const char*             title,
 265                     const std::string&      xaxis,
 266                     const std::string&      yaxis,
 267                     int                     exvg_graph_type,
 268                     const gmx_output_env_t* oenv)
 269 {
 270     FILE* fp;
 271
 272     fp = gmx_fio_fopen(fn, "w");
 273
 274     xvgr_header(fp, title, xaxis, yaxis, exvg_graph_type, oenv);
 275
 276     return fp;
 277 }
 278
 279 FILE* xvgropen(const char*             fn,
 280                const char*             title,
 281                const std::string&      xaxis,
 282                const std::string&      yaxis,
 283                const gmx_output_env_t* oenv)
 284 {
 285     return xvgropen_type(fn, title, xaxis, yaxis, exvggtXNY, oenv);
 286 }
 287
 288 void xvgrclose(FILE* fp)
 289 {
 290     gmx_fio_fclose(fp);
 291 }
 292
 293 void xvgr_subtitle(FILE* out, const char* subtitle, const gmx_output_env_t* oenv)
 294 {
 295     char buf[STRLEN];
 296
 297     if (output_env_get_print_xvgr_codes(oenv))
 298     {
 299         fprintf(out, "@ subtitle \"%s\"\n", xvgrstr(subtitle, oenv, buf, STRLEN));
 300     }
 301 }
 302
 303 void xvgr_view(FILE* out, real xmin, real ymin, real xmax, real ymax, const gmx_output_env_t* oenv)
 304 {
 305     if (output_env_get_print_xvgr_codes(oenv))
 306     {
 307         fprintf(out, "@ view %g, %g, %g, %g\n", xmin, ymin, xmax, ymax);
 308     }
 309 }
 310
 311 void xvgr_world(FILE* out, real xmin, real ymin, real xmax, real ymax, const gmx_output_env_t* oenv)
 312 {
 313     if (output_env_get_print_xvgr_codes(oenv))
 314     {
 315         fprintf(out,
 316                 "@ world xmin %g\n"
 317                 "@ world ymin %g\n"
 318                 "@ world xmax %g\n"
 319                 "@ world ymax %g\n",
 320                 xmin, ymin, xmax, ymax);
 321     }
 322 }
 323
 324 static bool stringIsEmpty(const std::string& s)
 325 {
 326     return s.empty();
 327 }
 328
 329 static bool stringIsEmpty(const char* s)
 330 {
 331     return (s == nullptr || s[0] == '\0');
 332 }
 333
 334 template<typename T>
 335 static void xvgr_legend(FILE* out, int nsets, const T* setname, const gmx_output_env_t* oenv)
 336 {
 337     int  i;
 338     char buf[STRLEN];
 339
 340     if (output_env_get_print_xvgr_codes(oenv))
 341     {
 342         xvgr_view(out, 0.15, 0.15, 0.75, 0.85, oenv);
 343         fprintf(out, "@ legend on\n");
 344         fprintf(out, "@ legend box on\n");
 345         fprintf(out, "@ legend loctype view\n");
 346         fprintf(out, "@ legend %g, %g\n", 0.78, 0.8);
 347         fprintf(out, "@ legend length %d\n", 2);
 348         for (i = 0; (i < nsets); i++)
 349         {
 350             if (!stringIsEmpty(setname[i]))
 351             {
 352                 if (output_env_get_xvg_format(oenv) == exvgXMGR)
 353                 {
 354                     fprintf(out, "@ legend string %d \"%s\"\n", i, xvgrstr(setname[i], oenv, buf, STRLEN));
 355                 }
 356                 else
 357                 {
 358                     fprintf(out, "@ s%d legend \"%s\"\n", i, xvgrstr(setname[i], oenv, buf, STRLEN));
 359                 }
 360             }
 361         }
 362     }
 363 }
 364
 365 void xvgrLegend(FILE* out, const std::vector<std::string>& setNames, const struct gmx_output_env_t* oenv)
 366 {
 367     xvgr_legend(out, setNames.size(), setNames.data(), oenv);
 368 }
 369 void xvgr_legend(FILE* out, int nsets, const char* const* setnames, const struct gmx_output_env_t* oenv)
 370 {
 371     xvgr_legend<const char*>(out, nsets, setnames, oenv);
 372 }
 373
 374 void xvgr_new_dataset(FILE* out, int nr_first, int nsets, const char** setname, const gmx_output_env_t* oenv)
 375 {
 376     int  i;
 377     char buf[STRLEN];
 378
 379     if (output_env_get_print_xvgr_codes(oenv))
 380     {
 381         fprintf(out, "@\n");
 382         for (i = 0; (i < nsets); i++)
 383         {
 384             if (setname[i])
 385             {
 386                 if (output_env_get_xvg_format(oenv) == exvgXMGR)
 387                 {
 388                     fprintf(out, "@ legend string %d \"%s\"\n", i + nr_first,
 389                             xvgrstr(setname[i], oenv, buf, STRLEN));
 390                 }
 391                 else
 392                 {
 393                     fprintf(out, "@ s%d legend \"%s\"\n", i + nr_first,
 394                             xvgrstr(setname[i], oenv, buf, STRLEN));
 395                 }
 396             }
 397         }
 398     }
 399     else
 400     {
 401         fprintf(out, "\n");
 402     }
 403 }
 404
 405 void xvgr_line_props(FILE* out, int NrSet, int LineStyle, int LineColor, const gmx_output_env_t* oenv)
 406 {
 407     if (output_env_get_print_xvgr_codes(oenv))
 408     {
 409         fprintf(out, "@    with g0\n");
 410         fprintf(out, "@    s%d linestyle %d\n", NrSet, LineStyle);
 411         fprintf(out, "@    s%d color %d\n", NrSet, LineColor);
 412     }
 413 }
 414
 415 static const char* LocTypeStr[] = { "view", "world" };
 416 static const char* BoxFillStr[] = { "none", "color", "pattern" };
 417
 418 void xvgr_box(FILE*                   out,
 419               int                     LocType,
 420               real                    xmin,
 421               real                    ymin,
 422               real                    xmax,
 423               real                    ymax,
 424               int                     LineStyle,
 425               int                     LineWidth,
 426               int                     LineColor,
 427               int                     BoxFill,
 428               int                     BoxColor,
 429               int                     BoxPattern,
 430               const gmx_output_env_t* oenv)
 431 {
 432     if (output_env_get_print_xvgr_codes(oenv))
 433     {
 434         fprintf(out, "@with box\n");
 435         fprintf(out, "@    box on\n");
 436         fprintf(out, "@    box loctype %s\n", LocTypeStr[LocType]);
 437         fprintf(out, "@    box %g, %g, %g, %g\n", xmin, ymin, xmax, ymax);
 438         fprintf(out, "@    box linestyle %d\n", LineStyle);
 439         fprintf(out, "@    box linewidth %d\n", LineWidth);
 440         fprintf(out, "@    box color %d\n", LineColor);
 441         fprintf(out, "@    box fill %s\n", BoxFillStr[BoxFill]);
 442         fprintf(out, "@    box fill color %d\n", BoxColor);
 443         fprintf(out, "@    box fill pattern %d\n", BoxPattern);
 444         fprintf(out, "@box def\n");
 445     }
 446 }
 447
 448 /* reads a line into ptr, adjusting len and renewing ptr if neccesary */
 449 static char* fgets3(FILE* fp, char** ptr, int* len, int maxlen)
 450 {
 451     int len_remaining = *len; /* remaining amount of allocated bytes in buf */
 452     int curp          = 0;    /* current position in buf to read into */
 453
 454     do
 455     {
 456         if (len_remaining < 2)
 457         {
 458             if (*len + STRLEN < maxlen)
 459             {
 460                 /* This line is longer than len characters, let's increase len! */
 461                 *len += STRLEN;
 462                 len_remaining += STRLEN;
 463                 srenew(*ptr, *len);
 464             }
 465             else
 466             {
 467                 /*something is wrong, we'll just keep reading and return NULL*/
 468                 len_remaining = STRLEN;
 469                 curp          = 0;
 470             }
 471         }
 472         if (fgets(*ptr + curp, len_remaining, fp) == nullptr)
 473         {
 474             /* if last line, skip */
 475             return nullptr;
 476         }
 477         curp += len_remaining - 1; /* overwrite the nul char in next iteration */
 478         len_remaining = 1;
 479     } while ((std::strchr(*ptr, '\n') == nullptr) && (feof(fp) == 0));
 480
 481     if (*len + STRLEN >= maxlen)
 482     {
 483         return nullptr; /* this line was too long */
 484     }
 485
 486     if (feof(fp))
 487     {
 488         /* We reached EOF before '\n', skip this last line. */
 489         return nullptr;
 490     }
 491     {
 492         /* now remove newline */
 493         int slen = std::strlen(*ptr);
 494         if ((*ptr)[slen - 1] == '\n')
 495         {
 496             (*ptr)[slen - 1] = '\0';
 497         }
 498     }
 499
 500     return *ptr;
 501 }
 502
 503 static int wordcount(char* ptr)
 504 {
 505     int i, n = 0, is[2];
 506     int cur = 0;
 507 #define prev (1 - cur)
 508
 509     if (nullptr != ptr)
 510     {
 511         for (i = 0; (ptr[i] != '\0'); i++)
 512         {
 513             is[cur] = std::isspace(ptr[i]);
 514             if ((0 == i) && !is[cur])
 515             {
 516                 n++;
 517             }
 518             else if ((i > 0) && (!is[cur] && is[prev]))
 519             {
 520                 n++;
 521             }
 522             cur = prev;
 523         }
 524     }
 525     return n;
 526 }
 527
 528 static char* read_xvgr_string(const char* line)
 529 {
 530     const char *ptr0, *ptr1;
 531     char*       str;
 532
 533     ptr0 = std::strchr(line, '"');
 534     if (ptr0 != nullptr)
 535     {
 536         ptr0++;
 537         ptr1 = std::strchr(ptr0, '"');
 538         if (ptr1 != nullptr)
 539         {
 540             str              = gmx_strdup(ptr0);
 541             str[ptr1 - ptr0] = '\0';
 542         }
 543         else
 544         {
 545             str = gmx_strdup("");
 546         }
 547     }
 548     else
 549     {
 550         str = gmx_strdup("");
 551     }
 552
 553     return str;
 554 }
 555
 556 int read_xvg_legend(const char* fn, double*** y, int* ny, char** subtitle, char*** legend)
 557 {
 558     FILE*    fp;
 559     char*    ptr;
 560     char*    base = nullptr;
 561     char*    fmt  = nullptr;
 562     int      k, line = 0, nny, nx, maxx, rval, legend_nalloc, set, nchar;
 563     double   lf;
 564     double** yy = nullptr;
 565     char*    tmpbuf;
 566     int      len = STRLEN;
 567     *ny          = 0;
 568     nny          = 0;
 569     nx           = 0;
 570     maxx         = 0;
 571     fp           = gmx_fio_fopen(fn, "r");
 572
 573     snew(tmpbuf, len);
 574     if (subtitle != nullptr)
 575     {
 576         *subtitle = nullptr;
 577     }
 578     legend_nalloc = 0;
 579     if (legend != nullptr)
 580     {
 581         *legend = nullptr;
 582     }
 583
 584     while ((ptr = fgets3(fp, &tmpbuf, &len, 10 * STRLEN)) != nullptr && ptr[0] != '&')
 585     {
 586         line++;
 587         trim(ptr);
 588         if (ptr[0] == '@')
 589         {
 590             if (legend != nullptr)
 591             {
 592                 ptr++;
 593                 trim(ptr);
 594                 set = -1;
 595                 if (std::strncmp(ptr, "subtitle", 8) == 0)
 596                 {
 597                     ptr += 8;
 598                     if (subtitle != nullptr)
 599                     {
 600                         *subtitle = read_xvgr_string(ptr);
 601                     }
 602                 }
 603                 else if (std::strncmp(ptr, "legend string", 13) == 0)
 604                 {
 605                     ptr += 13;
 606                     sscanf(ptr, "%d%n", &set, &nchar);
 607                     ptr += nchar;
 608                 }
 609                 else if (ptr[0] == 's')
 610                 {
 611                     ptr++;
 612                     sscanf(ptr, "%d%n", &set, &nchar);
 613                     ptr += nchar;
 614                     trim(ptr);
 615                     if (std::strncmp(ptr, "legend", 6) == 0)
 616                     {
 617                         ptr += 6;
 618                     }
 619                     else
 620                     {
 621                         set = -1;
 622                     }
 623                 }
 624                 if (set >= 0)
 625                 {
 626                     if (set >= legend_nalloc)
 627                     {
 628                         legend_nalloc = set + 1;
 629                         srenew(*legend, legend_nalloc);
 630                         (*legend)[set] = read_xvgr_string(ptr);
 631                     }
 632                 }
 633             }
 634         }
 635         else if (ptr[0] != '#')
 636         {
 637             if (nny == 0)
 638             {
 639                 (*ny) = nny = wordcount(ptr);
 640                 /* fprintf(stderr,"There are %d columns in your file\n",nny);*/
 641                 if (nny == 0)
 642                 {
 643                     return 0;
 644                 }
 645                 snew(yy, nny);
 646                 snew(fmt, 3 * nny + 1);
 647                 snew(base, 3 * nny + 1);
 648             }
 649             /* Allocate column space */
 650             if (nx >= maxx)
 651             {
 652                 maxx += 1024;
 653                 for (k = 0; (k < nny); k++)
 654                 {
 655                     srenew(yy[k], maxx);
 656                 }
 657             }
 658             /* Initiate format string */
 659             fmt[0]  = '\0';
 660             base[0] = '\0';
 661
 662             /* fprintf(stderr,"ptr='%s'\n",ptr);*/
 663             for (k = 0; (k < nny); k++)
 664             {
 665                 std::strcpy(fmt, base);
 666                 std::strcat(fmt, "%lf");
 667                 rval = sscanf(ptr, fmt, &lf);
 668                 /* fprintf(stderr,"rval = %d\n",rval);*/
 669                 if ((rval == EOF) || (rval == 0))
 670                 {
 671                     break;
 672                 }
 673                 yy[k][nx] = lf;
 674                 srenew(fmt, 3 * (nny + 1) + 1);
 675                 srenew(base, 3 * nny + 1);
 676                 std::strcat(base, "%*s");
 677             }
 678             if (k != nny)
 679             {
 680                 fprintf(stderr, "Only %d columns on line %d in file %s\n", k, line, fn);
 681                 for (; (k < nny); k++)
 682                 {
 683                     yy[k][nx] = 0.0;
 684                 }
 685             }
 686             nx++;
 687         }
 688     }
 689     gmx_fio_fclose(fp);
 690
 691     *y = yy;
 692     sfree(tmpbuf);
 693     sfree(base);
 694     sfree(fmt);
 695
 696     if (legend_nalloc > 0)
 697     {
 698         if (*ny - 1 > legend_nalloc)
 699         {
 700             assert(legend);
 701             srenew(*legend, *ny - 1);
 702             for (set = legend_nalloc; set < *ny - 1; set++)
 703             {
 704                 (*legend)[set] = nullptr;
 705             }
 706         }
 707     }
 708
 709     return nx;
 710 }
 711
 712 int read_xvg(const char* fn, double*** y, int* ny)
 713 {
 714     return read_xvg_legend(fn, y, ny, nullptr, nullptr);
 715 }
 716
 717 void write_xvg(const char* fn, const char* title, int nx, int ny, real** y, const char** leg, const gmx_output_env_t* oenv)
 718 {
 719     FILE* fp;
 720     int   i, j;
 721
 722     fp = xvgropen(fn, title, "X", "Y", oenv);
 723     if (leg)
 724     {
 725         xvgr_legend(fp, ny - 1, leg, oenv);
 726     }
 727     for (i = 0; (i < nx); i++)
 728     {
 729         for (j = 0; (j < ny); j++)
 730         {
 731             fprintf(fp, "  %12.5e", y[j][i]);
 732         }
 733         fprintf(fp, "\n");
 734     }
 735     xvgrclose(fp);
 736 }
 737
 738 real** read_xvg_time(const char* fn,
 739                      gmx_bool    bHaveT,
 740                      gmx_bool    bTB,
 741                      real        tb,
 742                      gmx_bool    bTE,
 743                      real        te,
 744                      int         nsets_in,
 745                      int*        nset,
 746                      int*        nval,
 747                      real*       dt,
 748                      real**      t)
 749 {
 750     FILE* fp;
 751 #define MAXLINELEN 16384
 752     char     line0[MAXLINELEN];
 753     char*    line;
 754     int      t_nalloc, *val_nalloc, a, narg, n, sin, set, nchar;
 755     double   dbl;
 756     gmx_bool bEndOfSet, bTimeInRange, bFirstLine = TRUE;
 757     real**   val;
 758
 759     t_nalloc   = 0;
 760     *t         = nullptr;
 761     val        = nullptr;
 762     val_nalloc = nullptr;
 763     *nset      = 0;
 764     *dt        = 0;
 765     fp         = gmx_fio_fopen(fn, "r");
 766     for (sin = 0; sin < nsets_in; sin++)
 767     {
 768         if (nsets_in == 1)
 769         {
 770             narg = 0;
 771         }
 772         else
 773         {
 774             narg = bHaveT ? 2 : 1;
 775         }
 776         n         = 0;
 777         bEndOfSet = FALSE;
 778         while (!bEndOfSet && fgets(line0, MAXLINELEN, fp))
 779         {
 780             line = line0;
 781             /* Remove whitespace */
 782             while (line[0] == ' ' || line[0] == '\t')
 783             {
 784                 line++;
 785             }
 786             bEndOfSet = (line[0] == '&');
 787             if (line[0] != '#' && line[0] != '@' && line[0] != '\n' && !bEndOfSet)
 788             {
 789                 if (bFirstLine && bHaveT)
 790                 {
 791                     /* Check the first line that should contain data */
 792                     a = sscanf(line, "%lf%lf", &dbl, &dbl);
 793                     if (a == 0)
 794                     {
 795                         gmx_fatal(FARGS, "Expected a number in %s on line:\n%s", fn, line0);
 796                     }
 797                     else if (a == 1)
 798                     {
 799                         fprintf(stderr,
 800                                 "Found only 1 number on line, "
 801                                 "assuming no time is present.\n");
 802                         bHaveT = FALSE;
 803                         if (nsets_in > 1)
 804                         {
 805                             narg = 1;
 806                         }
 807                     }
 808                 }
 809
 810                 a            = 0;
 811                 bTimeInRange = TRUE;
 812                 while ((a < narg || (nsets_in == 1 && n == 0))
 813                        && sscanf(line, "%lf%n", &dbl, &nchar) == 1 && bTimeInRange)
 814                 {
 815                     /* Use set=-1 as the time "set" */
 816                     if (sin)
 817                     {
 818                         if (!bHaveT || (a > 0))
 819                         {
 820                             set = sin;
 821                         }
 822                         else
 823                         {
 824                             set = -1;
 825                         }
 826                     }
 827                     else
 828                     {
 829                         if (!bHaveT)
 830                         {
 831                             set = a;
 832                         }
 833                         else
 834                         {
 835                             set = a - 1;
 836                         }
 837                     }
 838                     if (set == -1 && ((bTB && dbl < tb) || (bTE && dbl > te)))
 839                     {
 840                         bTimeInRange = FALSE;
 841                     }
 842
 843                     if (bTimeInRange)
 844                     {
 845                         if (n == 0)
 846                         {
 847                             if (nsets_in == 1)
 848                             {
 849                                 narg++;
 850                             }
 851                             if (set >= 0)
 852                             {
 853                                 *nset = set + 1;
 854                                 srenew(val, *nset);
 855                                 srenew(val_nalloc, *nset);
 856                                 val_nalloc[set] = 0;
 857                                 val[set]        = nullptr;
 858                             }
 859                         }
 860                         if (set == -1)
 861                         {
 862                             if (sin == 0)
 863                             {
 864                                 if (n >= t_nalloc)
 865                                 {
 866                                     t_nalloc = over_alloc_small(n);
 867                                     srenew(*t, t_nalloc);
 868                                 }
 869                                 (*t)[n] = dbl;
 870                             }
 871                             /* else we should check the time of the next sets with set 0 */
 872                         }
 873                         else
 874                         {
 875                             if (n >= val_nalloc[set])
 876                             {
 877                                 val_nalloc[set] = over_alloc_small(n);
 878                                 srenew(val[set], val_nalloc[set]);
 879                             }
 880                             val[set][n] = static_cast<real>(dbl);
 881                         }
 882                     }
 883                     a++;
 884                     line += nchar;
 885                 }
 886                 if (line0[strlen(line0) - 1] != '\n')
 887                 {
 888                     fprintf(stderr, "File %s does not end with a newline, ignoring the last line\n", fn);
 889                 }
 890                 else if (bTimeInRange)
 891                 {
 892                     if (a == 0)
 893                     {
 894                         fprintf(stderr, "Ignoring invalid line in %s:\n%s", fn, line0);
 895                     }
 896                     else
 897                     {
 898                         if (a != narg)
 899                         {
 900                             fprintf(stderr,
 901                                     "Invalid line in %s:\n%s"
 902                                     "Using zeros for the last %d sets\n",
 903                                     fn, line0, narg - a);
 904                         }
 905                         n++;
 906                     }
 907                 }
 908                 if (a > 0)
 909                 {
 910                     bFirstLine = FALSE;
 911                 }
 912             }
 913         }
 914         if (sin == 0)
 915         {
 916             *nval = n;
 917             if (!bHaveT)
 918             {
 919                 snew(*t, n);
 920                 for (a = 0; a < n; a++)
 921                 {
 922                     (*t)[a] = a;
 923                 }
 924             }
 925             if (n > 1)
 926             {
 927                 *dt = static_cast<real>((*t)[n - 1] - (*t)[0]) / (n - 1.0);
 928             }
 929             else
 930             {
 931                 *dt = 1;
 932             }
 933         }
 934         else
 935         {
 936             if (n < *nval)
 937             {
 938                 fprintf(stderr, "Set %d is shorter (%d) than the previous set (%d)\n", sin + 1, n, *nval);
 939                 *nval = n;
 940                 fprintf(stderr, "Will use only the first %d points of every set\n", *nval);
 941             }
 942         }
 943     }
 944     gmx_fio_fclose(fp);
 945
 946     sfree(val_nalloc);
 947
 948     return val;
 949 }