DataOnDemand.h 6.23 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
#ifndef DATAONDEMAND_H
#define DATAONDEMAND_H

#include "IGBheader.h"
#include <map>
#include <string>
#include <exception>
#include <algorithm>
#include "DataClass.h"
#include "AsciiReader.hpp"

using namespace std;

static const size_t IGB_OFFSET = 1024;

/** Class for reading data requested from the file on disk */
template<class T>
class DataOnDemand: public DataClass<T>
{
  public:
vigmond's avatar
vigmond committed
21 22 23 24 25 26 27
    virtual T      max(int t){return extremum(t,true);}	//!< maximum data value at a time
    virtual T      max(){assert(0);}                    //!< maximum data value at a time
    virtual T      min(int t){return extremum(t,false);}//!< minimum data value at a time
    virtual T      min(){assert(0);}                    //!< minimum data value at a time
    virtual T*     slice(int);                          //!< pointer to time slice of data
    virtual void   time_series( int, T* );              //!< time series for a point
    virtual void   increment(int){}                     //!< time slice increment
28 29 30 31 32 33 34 35

    DataOnDemand(const char *fn, int);
    virtual ~DataOnDemand(){free(data);if(databuf)free(databuf);}

  private:
    IGBheader   *hdr;
    AsciiReader *asciiRdr=NULL;
    char        *databuf=NULL;
vigmond's avatar
vigmond committed
36 37 38 39
    const char  *_scanstr;             //!< format code to read ascii file
    const int    bufsz=1024;            
    char         buf[1024];            //!< storage tor read ascii line
    T            extremum(int,bool);
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
    using DataClass<T> :: data;
    using DataClass<T> :: maxtm;
    using DataClass<T> :: last_tm;
    using DataClass<T> :: slice_size;
    using DataClass<T> :: filename;
    using DataClass<T> :: _dt;
    using DataClass<T> :: _t0;
    using DataClass<T> :: _ftype;
};

template<class T>
DataOnDemand<T>::DataOnDemand( const char *fn, int slsz )
{
  bool            IGBdata;
  int             j = 0;
  FILE           *in;
  string          scanner;
  map<int,string> CGfiles;
  map<int,string>::iterator CGp;
#ifdef USE_HDF5
  hid_t                  hin;
  struct ch5s_nodal_grid info;
  unsigned int           index;
#endif

  filename       = fn;
  slice_size     = slsz;
  _ftype         = FileTypeFinder( fn );

  if ( _ftype == FThdf5 )  {
#ifdef USE_HDF5
    if(ch5_open( fname.substr(0, filename.find_last_of(":")).c_str(), &hin ) )
      throw(1);
#endif
  } else {
    if ( (in=fopen(fn, "r")) == NULL )
      throw( 1 );
  }

  // ugly but I don't know what else to do besides specialization which is ugly
vigmond's avatar
vigmond committed
80 81 82 83 84
  if (      typeid(T) == typeid(double) ) _scanstr = "%lf";
  else if ( typeid(T) == typeid(float) )  _scanstr = "%f";
  else if ( typeid(T) == typeid(int) )    _scanstr = "%d";
  else if ( typeid(T) == typeid(short) )  _scanstr = "%hd";
  else if ( typeid(T) == typeid(long) )   _scanstr = "%ld";
85 86 87 88 89 90 91 92 93
  else exit(1);

  // open file
  if ( _ftype == FTIGB ) {
    hdr = new IGBheader( in, true );
    if ( hdr->slice_sz() != slice_size ) {
      maxtm = -1;
      throw PointMismatch(slice_size,hdr->slice_sz());
    }
vigmond's avatar
vigmond committed
94
    maxtm = hdr->t();
95 96 97 98 99 100 101
    _dt = hdr->inc_t();
    _t0 = hdr->org_t();
    databuf = (char *)malloc( slice_size*hdr->data_size() );
  } else if ( _ftype == FTfileSeqCG ) {
    CG_file_list( CGfiles, fn );
    CGp = CGfiles.begin();
    scanner = "%*f %*f";
vigmond's avatar
vigmond committed
102
    scanner += _scanstr;
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133
  } else if( _ftype == FThdf5 ) {
    string gtype;
#ifdef USE_HDF5
    if( parse_HDF5_grid( fn, gtype, index ) || gtype!="nodal" )
      throw(1);
    ch5s_nodal_grid_info( hin, index, &info );
    _dt = info.time_delta;
    maxtm = 0;
#endif
    _t0 = 0;
  } else if( _ftype == FTascii ) {
    asciiRdr = new AsciiReader;
    asciiRdr->read_file( filename );
    if( asciiRdr->num_lines() < slice_size ){
      maxtm = -1;
      throw PointMismatch(slice_size,asciiRdr->num_lines());
    }
    maxtm = asciiRdr->num_lines()/slice_size - 1;
  }

  data    = (T *)malloc( slice_size*sizeof(T) );
  last_tm = 0;
  
  // get one slice
  switch( _ftype ) {
      case FTIGB:
          hdr->read_data( data, 1, databuf );
          break;
      case FTascii:
          for( int i=0; i<slice_size; i++ ) {
            asciiRdr->get_line( i, buf, bufsz );
vigmond's avatar
vigmond committed
134
            if( sscanf(buf, _scanstr,  data+i ) != 1 ) {
135 136 137 138 139 140 141 142 143
              maxtm = -1;
              throw( 1 );
            }
          }
          break;
#ifdef USE_HDF5
      case FThdf5:
          assert(0);
#endif
144 145
      default:
          assert(0);
146 147 148 149
  }
}


vigmond's avatar
vigmond committed
150 151 152 153 154
/** return an extremum 
 *
 * \param tm  frame number
 * \param max determine max?, min o.w.
 */
155
template<class T>
vigmond's avatar
vigmond committed
156
T DataOnDemand<T>::extremum( int tm, bool max )
157
{
vigmond's avatar
vigmond committed
158
  if ( tm>maxtm ) return 0;
159

vigmond's avatar
vigmond committed
160 161 162 163 164
  T* (*exfcn)(T*, T*);
  if( max )
    exfcn = max_element<T*>;
  else
    exfcn = min_element<T*>;
165 166

  if( last_tm == tm )
vigmond's avatar
vigmond committed
167
    return *exfcn( data, data+slice_size );
168 169 170 171 172 173
 
  // analyze a temporary slice
  T* sdata = (T*)malloc(slice_size*sizeof(T));
  if( _ftype == FTIGB ) {
    hdr->data_frame( tm );
    hdr->read_data( sdata, 1, NULL );
174
  } else if( _ftype == FTascii ) {
175 176
    for( int i=0; i<slice_size; i++ ) {
      asciiRdr->get_line( i+slice_size*tm, buf, bufsz );
vigmond's avatar
vigmond committed
177
      sscanf( buf, _scanstr, sdata+i );
178 179 180 181 182
    }
  } else {
    assert(0);
  }

vigmond's avatar
vigmond committed
183
  T maxdat = *exfcn( sdata, sdata+slice_size ) ;
184
  free( sdata );
vigmond's avatar
vigmond committed
185
  return maxdat;
186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202
}


template<class T>
T* DataOnDemand<T>::slice( int tm )
{
  if ( tm>maxtm )
    return NULL;
  if( tm == last_tm )
    return data;

  if( _ftype == FTIGB ) {
    hdr->data_frame( tm );
    hdr->read_data( data, 1, databuf );
  } else if( _ftype == FTascii ) {
    for( int i=0; i<slice_size; i++ ) {
      asciiRdr->get_line( i+slice_size*tm, buf, bufsz );
vigmond's avatar
vigmond committed
203
      sscanf( buf, _scanstr, data+i );
204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229
    }
  } else
    assert(0);

  last_tm = tm;
  return data;
}


template<class T>
void DataOnDemand<T>::time_series( int offset, T* buffer )
{
  FILE *in = (FILE *)hdr->fileptr();
  if( _ftype == FTIGB ) {
    size_t slsz = hdr->slice_sz()*hdr->data_size();
    char datum[hdr->data_size()];
    offset *=hdr->data_size();
    for( size_t i=0; i<=maxtm; i++ ) {
       fseek( in, IGB_OFFSET+offset+i*slsz, SEEK_SET );
       fread( datum, hdr->data_size(), 1, in );
       if ( hdr->endian() != hdr->systeme() ) hdr->swab( datum, 1 );
       buffer[i] = hdr->convert_buffer_datum<T>( datum, 0 );
    }
  } else if( _ftype == FTascii ) {
    for( int i=0; i<=maxtm; i ++ ) {
      asciiRdr->get_line( offset+slice_size, buf, bufsz );
vigmond's avatar
vigmond committed
230
      sscanf( buf, _scanstr, buffer+i );
231 232 233 234
    }
  }
}
#endif