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