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