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