<commit_before>#ifndef ALEPH_TOPOLOGY_IO_VTK_HH__
#define ALEPH_TOPOLOGY_IO_VTK_HH__

#include <cstddef>

#include <algorithm>
#include <fstream>
#include <regex>
#include <stdexcept>
#include <sstream>
#include <string>
#include <vector>

#include "utilities/String.hh"

namespace aleph
{

namespace topology
{

namespace io
{

/**
  @class VTKStructuredGridReader
  @brief Simple reader class for VTK structured grids

  This class is a simple parser for VTK files in 'legacy format'. It is
  capable of parsing a structured grid and converting it to a simplicial
  complex. Data and weights of the simplicial complex will be taken from
  the VTK file.
*/

class VTKStructuredGridReader
{
public:
  template <class SimplicialComplex> void operator()( const std::string& filename, SimplicialComplex& K )
  {
    using Simplex    = typename SimplicialComplex::ValueType;
    using DataType   = typename Simplex::DataType;

    this->operator()( filename, K, [] ( DataType a, DataType b ) { return std::max(a,b); } );
  }

  template <class SimplicialComplex, class Functor> void operator()( const std::string& filename, SimplicialComplex& K, Functor f )
  {
    std::ifstream in( filename );
    if( !in )
      throw std::runtime_error( "Unable to read input file" );

    this->operator()( in, K, f );
  }

  template <class SimplicialComplex> void operator()( std::ifstream& in, SimplicialComplex& K )
  {
    using Simplex    = typename SimplicialComplex::ValueType;
    using DataType   = typename Simplex::DataType;

    this->operator()( in, K, [] ( DataType a, DataType b ) { return std::max(a,b); } );
  }

  template <class SimplicialComplex, class Functor> void operator()( std::ifstream& in, SimplicialComplex& K, Functor f )
  {
    using namespace aleph::utilities;

    using Simplex    = typename SimplicialComplex::ValueType;
    using DataType   = typename Simplex::DataType;
    using VertexType = typename Simplex::VertexType;

    std::string line;

    // Parse header first ----------------------------------------------

    std::size_t nx, ny, nz,
                n,
                s;

    bool parsedHeader = this->parseHeader( in, nx, ny, nz, n, s );
    if( !parsedHeader )
      return;

    // TODO: Check data type size against 's' and report if there are
    // issues such as insufficient storage space

    // Parse body ------------------------------------------------------
    //
    // The body contains coordinates for each of the points, which are
    // dutifully ignored for now, and attributes. For now, point-based
    // attributes are supported.

    std::regex rePointData( "POINT_DATA[[:space:]]+([[:digit:]]+)" );
    std::regex reScalars( "SCALARS[[:space:]]+([[:alnum:]]+)[[:space:]]+([[:alnum:]]+)[[:space:]]*([[:digit:]]*)" );
    std::regex reLookupTable( "LOOKUP_TABLE[[:space:]]+([[:alnum:]]+)" );

    std::vector<DataType> coordinates;
    coordinates.reserve( n*3 );

    while( std::getline( in, line ) )
    {
      line         = trim( line );
      auto strings = split( line );

      for( auto&& coordinate : strings )
        coordinates.push_back( convert<DataType>( coordinate ) );

      if( coordinates.size() == n*3 )
        break;
    }

    std::vector<DataType> values;
    values.reserve( n );

    {
      std::smatch matches;

      while( std::getline( in, line ) )
      {
        if( std::regex_match( line, matches, rePointData ) )
        {
          auto sn = matches[1];
          if( static_cast<std::size_t>( std::stoull( sn ) ) != n )
            throw std::runtime_error( "Format error: number of point data attributes does not match number of points" );
        }
        else if( std::regex_match( line, matches, reScalars ) )
        {
          auto name       = matches[1];
          auto type       = matches[2];
          auto components = matches[3];

          // TODO:
          //  - Use name
          //  - Check type
          //  - Check number of components (if present)

          (void) name;
          (void) type;
          (void) components;
        }
        else if( std::regex_match( line, matches, reLookupTable ) )
        {
          auto name = matches[1];
          if( name != "default" )
            throw std::runtime_error( "Handling non-default lookup tables is not yet implemented" );
        }
        else
        {
          line         = trim( line );
          auto strings = split( line );

          for( auto&& value : strings )
            values.push_back( convert<DataType>( value ) );
        }
      }
    }

    // Create topology -------------------------------------------------
    //
    // Notice that this class only adds 0-simplices and 1-simplices to
    // the simplicial complex for now. While it is possible to include
    // triangles (i.e. 2-simplices), their creation order is not clear
    // and may subtly influence calculations.

    std::vector<Simplex> simplices;

    // Create 0-simplices ----------------------------------------------

    {
      VertexType v = VertexType();

      for( auto&& value : values )
        simplices.push_back( Simplex(v, value) );
    }

    // Create 1-simplices ----------------------------------------------

    for( std::size_t z = 0; z < nz; z++ )
    {
      for( std::size_t y = 0; y < ny; y++ )
      {
        for( std::size_t x = 0; x < nx; x++ )
        {
          auto i = coordinatesToIndex(nx,ny,x,y,z);
          auto N = neighbours(nx,ny,nz,x,y,z);

          for( auto&& j : N )
          {
            auto wi = simplices.at(i).data();
            auto wj = simplices.at(j).data();

            // Use the functor specified by the client in order to
            // assign a weight for the new simplex.
            auto w  = f(wi, wj);

            simplices.push_back( Simplex( {i,j}, w ) );
          }
        }
      }<commit_msg>Refraining from creating duplicate neighbours<commit_after>#ifndef ALEPH_TOPOLOGY_IO_VTK_HH__
#define ALEPH_TOPOLOGY_IO_VTK_HH__

#include <cstddef>

#include <algorithm>
#include <fstream>
#include <regex>
#include <stdexcept>
#include <sstream>
#include <string>
#include <vector>

#include "utilities/String.hh"

namespace aleph
{

namespace topology
{

namespace io
{

/**
  @class VTKStructuredGridReader
  @brief Simple reader class for VTK structured grids

  This class is a simple parser for VTK files in 'legacy format'. It is
  capable of parsing a structured grid and converting it to a simplicial
  complex. Data and weights of the simplicial complex will be taken from
  the VTK file.
*/

class VTKStructuredGridReader
{
public:
  template <class SimplicialComplex> void operator()( const std::string& filename, SimplicialComplex& K )
  {
    using Simplex    = typename SimplicialComplex::ValueType;
    using DataType   = typename Simplex::DataType;

    this->operator()( filename, K, [] ( DataType a, DataType b ) { return std::max(a,b); } );
  }

  template <class SimplicialComplex, class Functor> void operator()( const std::string& filename, SimplicialComplex& K, Functor f )
  {
    std::ifstream in( filename );
    if( !in )
      throw std::runtime_error( "Unable to read input file" );

    this->operator()( in, K, f );
  }

  template <class SimplicialComplex> void operator()( std::ifstream& in, SimplicialComplex& K )
  {
    using Simplex    = typename SimplicialComplex::ValueType;
    using DataType   = typename Simplex::DataType;

    this->operator()( in, K, [] ( DataType a, DataType b ) { return std::max(a,b); } );
  }

  template <class SimplicialComplex, class Functor> void operator()( std::ifstream& in, SimplicialComplex& K, Functor f )
  {
    using namespace aleph::utilities;

    using Simplex    = typename SimplicialComplex::ValueType;
    using DataType   = typename Simplex::DataType;
    using VertexType = typename Simplex::VertexType;

    std::string line;

    // Parse header first ----------------------------------------------

    std::size_t nx, ny, nz,
                n,
                s;

    bool parsedHeader = this->parseHeader( in, nx, ny, nz, n, s );
    if( !parsedHeader )
      return;

    // TODO: Check data type size against 's' and report if there are
    // issues such as insufficient storage space

    // Parse body ------------------------------------------------------
    //
    // The body contains coordinates for each of the points, which are
    // dutifully ignored for now, and attributes. For now, point-based
    // attributes are supported.

    std::regex rePointData( "POINT_DATA[[:space:]]+([[:digit:]]+)" );
    std::regex reScalars( "SCALARS[[:space:]]+([[:alnum:]]+)[[:space:]]+([[:alnum:]]+)[[:space:]]*([[:digit:]]*)" );
    std::regex reLookupTable( "LOOKUP_TABLE[[:space:]]+([[:alnum:]]+)" );

    std::vector<DataType> coordinates;
    coordinates.reserve( n*3 );

    while( std::getline( in, line ) )
    {
      line         = trim( line );
      auto strings = split( line );

      for( auto&& coordinate : strings )
        coordinates.push_back( convert<DataType>( coordinate ) );

      if( coordinates.size() == n*3 )
        break;
    }

    std::vector<DataType> values;
    values.reserve( n );

    {
      std::smatch matches;

      while( std::getline( in, line ) )
      {
        if( std::regex_match( line, matches, rePointData ) )
        {
          auto sn = matches[1];
          if( static_cast<std::size_t>( std::stoull( sn ) ) != n )
            throw std::runtime_error( "Format error: number of point data attributes does not match number of points" );
        }
        else if( std::regex_match( line, matches, reScalars ) )
        {
          auto name       = matches[1];
          auto type       = matches[2];
          auto components = matches[3];

          // TODO:
          //  - Use name
          //  - Check type
          //  - Check number of components (if present)

          (void) name;
          (void) type;
          (void) components;
        }
        else if( std::regex_match( line, matches, reLookupTable ) )
        {
          auto name = matches[1];
          if( name != "default" )
            throw std::runtime_error( "Handling non-default lookup tables is not yet implemented" );
        }
        else
        {
          line         = trim( line );
          auto strings = split( line );

          for( auto&& value : strings )
            values.push_back( convert<DataType>( value ) );
        }
      }
    }

    // Create topology -------------------------------------------------
    //
    // Notice that this class only adds 0-simplices and 1-simplices to
    // the simplicial complex for now. While it is possible to include
    // triangles (i.e. 2-simplices), their creation order is not clear
    // and may subtly influence calculations.

    std::vector<Simplex> simplices;

    // Create 0-simplices ----------------------------------------------

    {
      VertexType v = VertexType();

      for( auto&& value : values )
        simplices.push_back( Simplex(v, value) );
    }

    // Create 1-simplices ----------------------------------------------

    for( std::size_t z = 0; z < nz; z++ )
    {
      for( std::size_t y = 0; y < ny; y++ )
      {
        for( std::size_t x = 0; x < nx; x++ )
        {
          auto i = coordinatesToIndex(nx,ny,x,y,z);
          auto N = neighbours(nx,ny,nz,x,y,z);

          for( auto&& j : N )
          {
            if( j > i )
              continue;

            auto wi = simplices.at(i).data();
            auto wj = simplices.at(j).data();

            // Use the functor specified by the client in order to
            // assign a weight for the new simplex.
            auto w  = f(wi, wj);

            simplices.push_back( Simplex( {i,j}, w ) );
          }
        }
      }
{'fraction_non_alphanumeric': 0.1611615861424811, 'fraction_numerical': 0.0025473380317568143, 'mean_word_length': 2.0355670103092782, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '68309', 'n_tokens_mistral': 3569, 'n_tokens_neox': 3400, 'n_words': 1040}