IGBops.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <https://www.gnu.org/licenses/>.
// ----------------------------------------------------------------------------


#include<stdlib.h>
#include<stdio.h>
#include<unistd.h>
#include<iostream>
#include "IGBheader.h"
#include <assert.h>
#include "ops_cmdline.h"
#include <vector>
#include "muParser.h"
#ifdef _OPENMP
#include <omp.h>
#endif

const int MAX_NC=4;

using namespace std;

enum Data_file_t { NONE=0, SPACE_TIME=5, SPACE=4, PT_TIME=3, PT=2, TIME=1 };
int  Tbit=1;  // time bit

class Interp {
    public:
        Interp( char *ind, char *weight );
        int size() { return ind.size(); }
        vector<float*>w;
        vector<int*>ind;
        int n=0;
};

Interp::Interp( char *indf, char *wf )
{
  char buf[8192];
  FILE *iin = fopen( indf, "r" );
  while( fgets( buf, 8192, iin )!=NULL ) {
    if( !n ) {
      char *l1 = strdup(buf);
      char *ptr = strtok( l1, " " );
      while((ptr=strtok( NULL, " " ))) n++;
      free(l1);
    }
    int *ip = new int[n];
    istringstream oss( buf );
    for( int i=0; i<n; i++ )
      oss>>ip[i];
    ind.push_back(ip);
  }
  FILE *win = fopen( wf, "r" );
  for( int i=0; i<ind.size(); i++ ) {
    fgets( buf, 8192, win );
    float *wp = new float[n];
    istringstream oss( buf );
    for( int j=0; j<n; j++ )
      oss >> wp[j];
    w.push_back( wp );
  }
}



#define data_t  float
int IGB_DATA_T[5] = { 0, IGB_FLOAT, IGB_COMPLEX, IGB_VEC3_f, IGB_VEC4_f };

void output_scratch( data_t* s, FILE *out, int n, bool bin );
void output_scratch( vector<data_t>& s, FILE *out, bool bin );


bool overlap(IGBheader *first, IGBheader *second )
{
  return fabs(second->org_t() - (first->org_t()+(first->t()-1)*first->inc_t()) ) < first->inc_t()/2.;
}


// concatenate X and Y in time
void joinXY( IGBheader *h1, IGBheader *h2, IGBheader *hout )
{
  if( h2==NULL ) {
    fprintf( stderr, "Y must be an IGB file\n\n" );
    exit(1);
  }

  gzFile in1 = (gzFile)(h1->fileptr());
  gzFile in2 = (gzFile)(h2->fileptr());
  FILE*  out = (FILE *)(hout->fileptr());

  std::vector<data_t> dat1(hout->slice_sz() * hout->num_components());
  for( int i=0; i<h1->t(); i++ ) {
    h1->read_data( dat1.data() );
    fwrite( dat1.data(), sizeof(data_t), hout->slice_sz() * hout->num_components(), out );
  }

  bool skip_second_torg = overlap( h1, h2 );
  for( int i=0; i<h2->t(); i++ ) {
    h2->read_data( dat1.data() );
    if(skip_second_torg)
      skip_second_torg = false;
    else
      fwrite( dat1.data(), sizeof(data_t), hout->slice_sz() * hout->num_components(), out );
  }
  fclose(out);
}

// put Y beside X
void join_space( data_t *x, int nx, data_t *y, int ny, FILE *out )
{
  fwrite( x, sizeof(data_t), nx, out );
  fwrite( y, sizeof(data_t), ny, out );
}


// evaluate arbitrary expression for one file
void eval_expr( data_t *x, int n, char *expr, FILE *out, data_t *result )
{
  using namespace mu;

  static Parser hParser;
  static bool init=false;

  if( !init ) {
    hParser.DefineVar( "X", x );
    hParser.SetExpr(expr);
    init = true;
  }
  try {
    hParser.Eval(result, n);
  }
  catch (Parser::exception_type &e)
  {
    std::cerr << e.GetMsg() << std::endl;
    exit(1);
  }
  output_scratch( result, out, n, true );
}


// evaluate arbitrary expression as function of two files
// nc>0 means that x1 is a point
template<class T>
void eval_expr( data_t *x0, T* x1, int n, char *expr, FILE *out, int nc=0 )
{
  using namespace mu;

  Parser hParser;
  data_t a, b;
  hParser.DefineVar( "X", &a );
  hParser.DefineVar( "Y", &b );
  hParser.SetExpr(expr);

  data_t *y = new data_t[n];

  for( int i=0; i<n; i++ ) {
    try {
      a = x0[i];
      b = x1[nc>0 ? i%nc : i];
      y[i] = hParser.Eval();
    }
    catch (Parser::exception_type &e)
    {
      std::cerr << e.GetMsg() << std::endl;
      exit(1);
    }
  }
  output_scratch( y, out, n, true );
  delete[] y;
}


// central difference
void center_diff( data_t* x, int slsz, vector<data_t>& odat, int t, float dt, FILE *out)
{
  int older= (t%2) ? 0 : slsz;
  int old  = older ? 0 : slsz;

  if( t==1 ) {  // forward difference for first frame
    for(int i=0;i<slsz;i++){
      data_t d = (x[i]-odat[old+i])/dt;
      fwrite( &d, sizeof(data_t), 1, out );
    }
  } else if( t>1 ) {
    for(int i=0;i<slsz;i++){
      data_t d = (x[i]-odat[older+i])/dt/2.;
      fwrite( &d, sizeof(data_t), 1, out );
    }
  }
  for(int i=0;i<slsz;i++)
    odat[older+i]   = x[i];
}


// forward difference
void forward_diff( data_t* x, int slsz, vector<data_t>& old, int t, float dt, FILE *out)
{
  for(int i=0;i<slsz;i++){
    if( t ) {
      data_t d = (x[i]-old[i])/dt;
      fwrite( &d, sizeof(data_t), 1, out );
    }
    old[i] = x[i];
  }
}


// aX[i] + bY[i%nc]
void aXPLUS_bY( data_t* x, data_t *y, data_t a, data_t b, int slsz, int nc, FILE *out)
{
  for(int i=0;i<slsz;i++){
    data_t datum = x[i]*a + b*y[i%nc];
    fwrite( &datum, sizeof(data_t), 1, out );
  }
}

// aX[i] + bY[i]
template<class T>
void aXPLUS_bY( data_t* x, T* y, data_t a, data_t b, int slsz, FILE *out)
{
  for(int i=0;i<slsz;i++){
    data_t datum = x[i]*a + b*y[i];
    fwrite( &datum, sizeof(data_t), 1, out );
  }
}

// avgX[*,t]
void avgXt( data_t* x, int slsz, int nc, float t, FILE *out)
{
  data_t sumf[MAX_NC];
  for( int c=0; c<nc; c++ ) {
    double s=0;
#pragma omp parallel for reduction( +:s )
    for(int i=c;i<slsz;i+=nc)
      s += x[i];
    sumf[c] = s/double(slsz);
  }
  fwrite( &sumf, sizeof(sumf[0]), nc, out );
}


// avgX[i,*]
void avgXi( data_t* x, vector<data_t>& sum, const data_t fac)
{
#pragma omp parallel for
  for(int i=0; i<sum.size(); i++){
    sum[i] += fac * x[i];
  }
}

void varXi( data_t* x, data_t* mean, vector<data_t>& sum, const data_t fac)
{
#pragma omp parallel for
  for(int i=0; i<sum.size(); i++){
    sum[i] += fac * (x[i] - mean[i]) * (x[i] - mean[i]);
  }
}

// aX[i,t]+b
void aXPLUS_b( data_t* x, data_t  a, data_t b, int slsz, FILE *out)
{
  for(int i=0;i<slsz;i++){
    data_t datum = x[i]*a+b;
    fwrite( &datum, sizeof(data_t), 1, out );
  }
}

// X[i]/Y[i]
void XdivY( data_t* x, data_t* y, int slsz, FILE *out)
{
  bool NAN_FLAG = 0;
  bool INF_FLAG = 0;

  for(int i=0;i<slsz;i++){
    data_t datum = x[i]/y[i];
    if (std::isinf(datum)) {
      datum = 0.;  // handle infinite numbers
      INF_FLAG = 1;
    }
    if (std::isnan(datum)){
      datum = 0.;  // handle not a number
      NAN_FLAG = 1;
    }
    fwrite( &datum, sizeof(data_t), 1, out );
  }
  if (INF_FLAG) std::cout << "Infinite number(s) detected (and set to zero)" << endl;
  if (NAN_FLAG) std::cout << "NaN(s) detected (and set to zero)" << endl;
}


// X[i]/Y[0]
void XdivY( data_t* x, data_t* y, int slsz, int nc, FILE *out)
{
  bool NAN_FLAG = 0;
  bool INF_FLAG = 0;

  for(int i=0;i<slsz;i++){
    data_t datum = x[i]/y[i%nc];
    if (std::isinf(datum)) {
      datum = 0.;  // handle infinite numbers
      INF_FLAG = 1;
    }
    if (std::isnan(datum)){
      datum = 0.;  // handle not a number
      NAN_FLAG = 1;
    }
    fwrite( &datum, sizeof(datum), 1, out );
  }
  if (INF_FLAG) std::cout << "Infinite number(s) detected (and set to zero)" << endl;
  if (NAN_FLAG) std::cout << "NaN(s) detected (and set to zero)" << endl;
}

// X[i]*Y[i]
void XmulY( data_t* x, data_t* y, int slsz, FILE *out)
{
  for(int i=0;i<slsz;i++){
    data_t datum = x[i]*y[i];
    fwrite( &datum, sizeof(data_t), 1, out );
  }
}

// X[i]*Y[0]
void XmulY( data_t* x, data_t* y, int slsz, int nc, FILE *out)
{
  for(int i=0;i<slsz;i++){
    data_t datum = x[i]*y[i%nc];
    fwrite( &datum, sizeof(datum), 1, out );
  }
}

// MIN_i( X[i,t] )
void minXt(  data_t* x, int slsz, int nc, float t, FILE *out )
{
  data_t min[MAX_NC];
  for( int i=0; i<nc; i++ ) min[i]=x[i];
#pragma omp parallel for
  for(int i=nc;i<slsz*nc;i++){
    if( x[i]<min[i%nc] )
      min[i%nc] = x[i];
  }
  fwrite( min, sizeof(data_t), nc, out );
}


// t(MAX_t( X[i,t] ))
void t_maxX(  data_t* x, data_t* mx, data_t* tm, size_t n, float t )
{
#pragma omp parallel for
  for( int i=0; i<n; i++ ){
    if( x[i]>mx[i] ){
      mx[i] = x[i];
      tm[i] = t;
    }
  }
}

// MAX_t( X[i,t] )
void maxXi(  data_t* x, data_t* mx, size_t n )
{
#pragma omp parallel for
  for( int i=0; i<n; i++ ){
    if( x[i]>mx[i] )
      mx[i] = x[i];
  }
}


// MIN_t( X[i,t] )
void minXi( data_t* x, data_t* mn, size_t n )
{
#pragma omp parallel for
  for( int i=0; i<n; i++ ){
    if( x[i]<mn[i] )
      mn[i] = x[i];
  }
}


// MAX_i( X[i,t] )
void maxXt(  data_t* x, int slsz, int nc, float t, FILE *out )
{
  data_t max[MAX_NC];
  for( int i=0; i<nc; i++ ) max[i]=x[i];
#pragma omp parallel for
  for(int i=nc;i<slsz*nc;i++){
    if( x[i]>max[i%nc] )
      max[i%nc] = x[i];
  }
  fwrite( max, sizeof(data_t), nc, out );
}

// X[i,t]+a
void XPLUS_a(   data_t* x, data_t  a, int slsz, FILE *out )
{
  for(int i=0;i<slsz;i++){
    data_t datum = x[i]+a;
    fwrite( &datum, sizeof(data_t), 1, out );
  }
}


// apply a simple box filter to the xXy grid where the center pixel
// is the average of the enclosing dXd square, do not filter edges
void Xbox( data_t *dat, int x, int y, FILE *out, int d, int nc )
{
  assert(nc==1);

  int   fm   = d/2,
        fp   = d - 1 - fm;
  float area = d*d;
  for( int j=0; j<y; j++ )
    for( int i=0; i<x; i++ ) {
      data_t datum = 0;
      if( j<fm || j>y-fp-1 || i<fm || i>x-fp-1 )
        datum = dat[j*x+i];
      else {
        for( int m=j-fm; m<=j+fp; m++ )
          for( int n=i-fm; n<=i+fp; n++ )
            datum += dat[m*x+n];
        datum /= area;
      }
      fwrite( &datum, sizeof(data_t), 1, out );
    }
}


// soft clipping: remove offending point and interpolate
void sclip_ab( data_t *dat, int sz, data_t upper, data_t lower, vector<data_t>& scr, int t,
                                   int tfinal, FILE *out )
{
  unsigned int curr  = (t%3)*sz;
  unsigned int old   = (curr+2*sz)%(3*sz);
  unsigned int older = (curr+sz)%(3*sz);

  memcpy( scr.data()+curr, dat, sizeof(data_t)*sz );

  if( !t )
    fwrite( dat, sizeof(data_t), sz, out );

  if( t>1 ) {
    for( int i=0; i<sz; i++ )
      if( scr[old+i]>upper || scr[old+i]<lower ) {
        if(  scr[curr+i]>upper || scr[curr+i]<lower )
          scr[old+i] = scr[older+i];
        else
          scr[old+i] = (scr[older+i]+scr[curr+i])/2.;
      }
    fwrite( scr.data()+old, sizeof(data_t), sz, out );
  }

  if( t==tfinal-1 )
    fwrite( dat, sizeof(data_t), sz, out );
}


void 
magnitude( data_t* x, int n, int nc, FILE *out)
{
  for(int i=0; i<n; i+=nc){
    data_t mag = 0;
    for( int j=0; j<nc; j++ )
      mag += x[i+j]*x[i+j];
    mag = sqrt( mag );
    fwrite( &mag, sizeof(mag), 1, out );
  }
}

void 
dot( data_t* x, data_t *y, int n, int nc, FILE *out, bool pt)
{
  int k = !pt;
  for(int i=0; i<n; i+=nc){
    data_t dot = 0;
    for( int j=0; j<nc; j++ )
      dot += x[i+j]*y[k*i+j];
    fwrite( &dot, sizeof(dot), 1, out );
  }
}

template<class T>
void
simp38( T *x, int n, int t, data_t *runInt, float dt )
{
  int s=t%3;
  T scale = dt;

  if( t<0 )
    scale *= 3./8.;
  else if( s==0 )
    scale *= 3./4.;
  else if( s==1 || s==2 )
    scale *= 9./8.;

#pragma omp parallel for
  for( int i=0; i<n; i++ )
    runInt[i] += x[i]*scale;
}


template<class T>
void
midpt( T *x, int n, bool end, data_t *reduce, float dt )
{
#pragma omp parallel for
  for( int i=0; i<n; i++ )
    if( end )
      reduce[i] += dt*x[i]/2.;
    else 
      reduce[i] += dt*x[i];
}


void
interpolate( data_t *dat, vector<int*>&idx, vector<float*>& w, int n, FILE *out )
{
  for( int i=0; i<idx.size(); i++ ) {
    data_t datum=0.;
    for( int j=0; j<n; j++ )
        datum += dat[idx[i][j]]*w[i][j];
    fwrite( &datum, sizeof(data_t), 1, out );
  }
}


template<class T>
void
normalize( IGBheader *h, FILE *out, T *min, T* max )
{
  size_t slice_mem = h->slice_sz()*h->num_components();
  h->go2slice(0);
  vector<data_t> slice( slice_mem );

  for( int t=0; t< h->t(); t++ ) {
    h->read_data( slice.data() );
    for( int j=0; j<slice_mem; j++ ){
      data_t datum = (slice[j]-min[j])/(max[j]-min[j]);
      fwrite( &datum, sizeof(data_t), 1, out );
    }
  }
}


void
maxP2P( data_t* dat, data_t* big, data_t* maxpp, int sz, int t )
{
  data_t *old    = big;
  data_t *change = big+sz;
  data_t *extrem = big+2*sz;

  if( !t ) memcpy( old, dat, sizeof(data_t)*sz );

#pragma omp parallel for
  for( int i=0; i<sz; i++ ) {
    data_t dV = dat[i] - old[i];
    if( dV*change[i] < 0 ) {
      data_t pp = fabs( extrem[i]-old[i] );
      if( pp > maxpp[i] ) maxpp[i] = pp;
      change[i] = copysign( 1., dV );
      extrem[i] = old[i];
    } else if( change[i]==0. && dV != 0 ) {
      change[i] = copysign( 1., dV );
    }
  }
  memcpy( old, dat, sizeof(data_t)*sz );
}


//output scratch vector
void
output_scratch( data_t* s, FILE *out, int n, bool bin )
{
  for( int i=0; i<n; i++ )
    if( bin ) {
      fwrite( s+i, sizeof(data_t), 1, out );
    }else
      fprintf( out, "%g\n", s[i] );
}

//output scratch vector
void
output_scratch( vector<data_t>& s, FILE *out, bool bin=false )
{
  for( auto i : s )
    if( bin ) {
      data_t datum = i;
      fwrite( &datum, sizeof(data_t), 1, out );
    }else
      fprintf( out, "%g\n", i );
}



// temporal running average
void ravgX_ta( data_t *dat1, int sz, int t, int a,  vector<data_t>& tmp, FILE *out )
{
  vector <data_t> sum( sz, 0 );

#pragma omp parallel for
  for( int i=0; i<sz; i++ ) {

    tmp[i+(t%a)*sz] = dat1[i];

    for( int j=0; j<a; j++ ) {
        sum[i] += tmp[i+j*sz];
    }
    sum[i] /= (float)a;
  }

  if( t<a-1 )
    return;

  output_scratch( sum, out, true );
}


//read in an ascii vector file
void
read_vector( vector<data_t> &y, gzFile in, int nc )
{
  assert( nc<=MAX_NC );
  char   buf[1024];
  float  datum[MAX_NC];
  while( gzgets( in, buf, 1024 ) != Z_NULL &&
           sscanf( buf, "%f %f %f %f", datum, datum+1, datum+2, datum+3 )==nc )
    for( int i=0; i<nc; i++ ) y.push_back( datum[i] );
}


std::string 
cmdline_str( int n, char *w[] )
{
  std::string cl;
  for( int i=0; i<n; i++ ) {
    cl += w[i];
    cl += " ";
  }
  return cl;
}


int
main( int argc, char *argv[] )
{
#ifdef _OPENMP
  if( !getenv("OMP_NUM_THREADS") )
    omp_set_num_threads(1);
#endif

  gengetopt_args_info opts;

  std::string cmdline = cmdline_str( argc, argv );

  // let's call our cmdline parser
  if (cmdline_parser (argc, argv, &opts) != 0)
    exit(1);

  // the last arg must be a file name
  gzFile in1 = (!opts.inputs_num || !strcmp(opts.inputs[0],"-")) ?
               gzdopen( dup(0), "r" ):       // duplicate stdin
               gzopen( opts.inputs[0], "r" );
  if( in1 == NULL ) {
    cerr << "File not found: " << opts.inputs[0] << endl;
    exit(1);
  }
  // make sure we can open X which must be an IGB file
  IGBheader* h1= new IGBheader( in1 );
  if( h1->read(true) & 0xff ) {
    cerr << "File 1 is not a proper IGB file: "<< opts.inputs[0] << endl;
    exit(1);
  }

  Interp *interp=NULL;
  if( opts.interpolate_flag )
    interp = new Interp( opts.inputs[1], opts.inputs[2] );

  // is a second file necessary?
  int numf=1;
  if( opts.op_given )
    switch( opts.op_arg ) {
        case op_arg_aXPLUS_bY:
        case op_arg_XdivY:
        case op_arg_XmulY:
        case op_arg_joinXY_t:
        case op_arg_joinXY_i:
        case op_arg_dotXY:
            numf = 2;
            break;
    }
  else if( opts.expression_given )
    if( strchr( opts.expression_arg, 'Y' ) )
      numf =2;

  // read header and determine size of second file if necessary
  gzFile        in2;
  IGBheader    *h2 = NULL;
  vector<float> Y;
  Data_file_t   Ytype = NONE;
  bool          Yigb=false;

  if( numf==2 ) {
    if( opts.inputs_num<2 ) {
      cerr << "A 2nd file is needed : " << endl;
      exit(1);
    }
    if( (in2=gzopen( opts.inputs[1], "r" ))==NULL ) {
      cerr << "file 2 not found : " << opts.inputs[1] << endl;
      exit(1);
    }
    // check file type
    h2 = new IGBheader( in2 );
    if( h2->read(true) & 0xff ) {
      gzrewind( in2 );
      read_vector( Y, in2, h1->num_components() );
      if( Y.size()==h1->slice_sz()*h1->num_components() )
        Ytype = SPACE;
      else if( Y.size()==h1->t()*h1->num_components() )
        Ytype = PT_TIME;
      else if( Y.size()==h1->num_components() )
        Ytype = PT;
      else if( opts.op_arg != op_arg_joinXY_i ) {
        cerr << "incompatile file sizes! Y must have " << h1->t()*h1->num_components() 
             << " or " << h1->slice_sz()*h1->num_components() << " or " << h1->num_components()
             << " entries" << endl;
        exit(1);
      }
    } else { 
      Yigb = true;
      if( h2->num_components() != h1->num_components() ) {
        cerr << "Number of data components do not match" << endl;
        exit(1);
      }
      if( h2->slice_sz()==h1->slice_sz() && h2->t()==h1->t() )
        Ytype = SPACE_TIME;
      else if( h2->slice_sz()==h1->slice_sz() ) {
        Ytype = SPACE;
        if( h2->t()!=1 && opts.op_arg != op_arg_joinXY_t ) {
          cerr << "Second IGB file time does not match" << endl;
          exit(1);
        }
        Y.resize(h1->slice_sz()*h1->num_components());
      } else if(  h2->slice_sz()==1 && h2->t()==h1->t() )
        Ytype = PT_TIME;
      else if ( h2->t()==h1->t() )
        Ytype = TIME;
      else if ( h2->slice_sz()==1 && h2->t()==1 )
        Ytype = PT;
      else if( opts.op_arg != op_arg_joinXY_i ) {
        cerr << "incompatile IGB file sizes!" << endl;
        exit(1);
      }
    }
  }

  // determine the output file size
  IGBheader *hout;
  FILE *out;
  string fname =  opts.output_file_arg;
  if( fname=="-" )
    out= stdout;
  else{
    if( fname.length()<5 || fname.substr(fname.length()-4,4)!=".igb" )
      fname += opts.extension_arg;
    out = fopen( fname.c_str(), "w" );
  }
  hout = new IGBheader;
  *hout = *h1;
  hout->type(IGB_DATA_T[h1->num_components()]);
  switch( opts.op_arg ) {
      case op_arg_ravgX_ta:
          hout->t(h1->t()-opts.a_arg+1);
          break;
      case op_arg_joinXY_i:
          if( Yigb )
            hout->x(h1->slice_sz()+h2->slice_sz());
          else
            hout->x(h1->slice_sz()+Y.size());
          hout->y(1);
          hout->z(1);
          break;
      case op_arg_joinXY_t:
          hout->type(h1->type());
          hout->t(h1->t()+h2->t()-overlap(h1,h2));
          break;
      case op_arg_diff_f:
      case op_arg_diff_c:
          hout->t(h1->t()-1);
          break;
      case op_arg_minX_i:
      case op_arg_maxX_i:
      case op_arg_t_maxX:
      case op_arg_avgX_i:
      case op_arg_stdD_i:
      case op_arg_intX_i:
      case op_arg_maxP2P:
          hout->t(1);
          hout->type(IGB_DATA_T[1]);
          break;
      case op_arg_minX_t:
      case op_arg_maxX_t:
      case op_arg_avgX_t:
      case op_arg_stdD_t:
          hout->x(1);
          hout->y(1);
          hout->z(1);
          break;
      case op_arg_magX:
      case op_arg_dotXY:
          hout->type(IGB_DATA_T[1]);
          break;
  }
  if( opts.interpolate_flag ) {
    hout->x( interp->size() );
    hout->y(1);
    hout->z(1);
  }
  //hout->comment( cmdline.c_str() );
  hout->fileptr( out );
  hout->write();

  if( opts.op_arg == op_arg_joinXY_t ) {
    joinXY( h1, h2, hout );
    exit(0);
  }

  // allocate some input memory
  std::vector<float> dat1(h1->slice_sz() * h1->num_components());
  std::vector<std::vector<float>> dat3;
  if( Yigb )
    Y.assign(h2->slice_sz()*h2->num_components(), 0.0f);

  // initialize if we need scratch data
  int h1_slice_comps = h1->slice_sz()*h1->num_components();
  std::vector<data_t> t_reduce; // automatically output at the end for a temporal reduction
  std::vector<data_t> big_scratch;
  switch( opts.op_arg ) {
    case op_arg_minX_i:
    case op_arg_maxX_i:
    case op_arg_avgX_i:
    case op_arg_stdD_i:
    case op_arg_intX_i :
    case op_arg_intX_t :
        t_reduce.assign( h1_slice_comps, 0 );
        break;
    case op_arg_t_maxX:
        t_reduce.assign( h1_slice_comps, 0 );
        big_scratch.assign( h1_slice_comps, 0 );
        break;
    case op_arg_ravgX_ta :
        big_scratch.assign( h1_slice_comps*opts.a_arg, 0 );
        break;
    case op_arg_sclip_ab :
        big_scratch.assign( h1_slice_comps*3, 0 );
        break;
    case op_arg_diff_f :
        big_scratch.assign( h1_slice_comps, 0 );
        break;
    case op_arg_diff_c :
    case op_arg_norm_i :
        big_scratch.assign( h1_slice_comps*3, 0 );
        break;
    case op_arg_maxP2P :
        big_scratch.assign( h1_slice_comps*3, 0 );
        t_reduce.assign( h1_slice_comps, 0 );
        break;
  }
  if( opts.expression_given )
    big_scratch.assign( h1_slice_comps*opts.a_arg, 0 );

  // loop through the data
  for( int t=0; t<h1->t(); t++ ) {

    float tm = h1->org_t() + t*h1->inc_t();

    h1->read_data( dat1.data() );

    if(opts.op_arg == op_arg_stdD_i)
      dat3.push_back(dat1);
    
    if( Yigb && (Ytype&Tbit || !t) )
      h2->read_data( Y.data() );

    float *Yptr = Y.data();
    if( !Yigb && Ytype==PT_TIME )
      Yptr += h1->num_components()*t;

    // first data read initialization
    if( !t ) { 
      if( opts.op_arg==op_arg_minX_i || opts.op_arg==op_arg_maxX_i ) {
        t_reduce = dat1;
      } else if( opts.op_arg==op_arg_norm_i ) {
        copy( dat1.begin(), dat1.end(), big_scratch.begin() );
        copy( dat1.begin(), dat1.end(), big_scratch.begin()+h1_slice_comps );
      } else if( opts.op_arg==op_arg_t_maxX)  {
        big_scratch = dat1;
      }
    }

    if( opts.expression_given ){
      if( Ytype==NONE )
        eval_expr( dat1.data(), h1_slice_comps, opts.expression_arg, out, big_scratch.data() );
      else {
        bool space = Ytype==SPACE_TIME || Ytype==SPACE;
        eval_expr( dat1.data(), Yptr, h1_slice_comps, opts.expression_arg, 
                                out, space? 0 : h1->num_components() );
      }
    }

    if( opts.interpolate_given )
        interpolate( dat1.data(), interp->ind, interp->w, interp->n, out );

    switch( opts.op_arg ) {
      case op_arg_aXPLUS_bY:
          if( Ytype==SPACE_TIME || Ytype==SPACE )
            aXPLUS_bY(dat1.data(), Yptr, opts.a_arg, opts.b_arg,
                                                   h1_slice_comps,out);
          else 
            aXPLUS_bY(dat1.data(), Yptr, opts.a_arg, opts.b_arg,
                                         h1_slice_comps,h1->num_components(), out);
          break;
      case op_arg_aXPLUS_b:
          aXPLUS_b(dat1.data(),opts.a_arg,opts.b_arg,h1_slice_comps,out);
          break;
      case op_arg_XdivY:
          if( Ytype==SPACE_TIME || Ytype==SPACE )
            XdivY(dat1.data(), Yptr, h1_slice_comps, out);
          else
            XdivY(dat1.data(), Yptr, h1_slice_comps, h1->num_components(), out);
          break;
      case op_arg_XmulY:
          if( Ytype==SPACE_TIME || Ytype==SPACE )
            XmulY( dat1.data(), Yptr, h1_slice_comps, out);
          else 
            XmulY( dat1.data(), Yptr, h1_slice_comps, h1->num_components(), out);
          break;
      case op_arg_minX_t:
          minXt( dat1.data(), h1->slice_sz(), h1->num_components(), tm, out );
          break;
      case op_arg_maxX_t:
          maxXt( dat1.data(), h1->slice_sz(), h1->num_components(), tm, out );
          break;
      case op_arg_avgX_t:
          avgXt( dat1.data(), h1->slice_sz(), h1->num_components(), tm, out );
          break;
      case op_arg_minX_i:
          minXi( dat1.data(), t_reduce.data(), t_reduce.size() );
          break;
      case op_arg_maxX_i:
          maxXi( dat1.data(), t_reduce.data(), t_reduce.size() );
          break;
      case op_arg_t_maxX:
          t_maxX( dat1.data(), big_scratch.data(), t_reduce.data(), t_reduce.size(), t );
          break;
      case op_arg_avgX_i: 
      case op_arg_stdD_i:
          avgXi( dat1.data(), t_reduce, 1.0f);
          break;
      case op_arg_ravgX_ta :
          ravgX_ta( dat1.data(), h1_slice_comps, t, opts.a_arg, big_scratch, out );
          break;
      case op_arg_Xbox :
          Xbox( dat1.data(), h1->x(), h1->y(), out, opts.a_arg, h1->num_components() );
          break;
      case op_arg_sclip_ab :
          sclip_ab( dat1.data(), h1_slice_comps, opts.a_arg, opts.b_arg, big_scratch, t, h1->t(), out );
          break;
      case op_arg_diff_f :
          forward_diff( dat1.data(), h1_slice_comps, big_scratch, t, h1->inc_t(), out );
          break;
      case op_arg_diff_c :
          center_diff( dat1.data(), h1_slice_comps, big_scratch, t, h1->inc_t(), out );
          break;
      case op_arg_joinXY_i:
          if( Yigb )
            join_space( dat1.data(), h1_slice_comps, Yptr, 
                                                 h2->slice_sz()*h2->num_components(), out);
          else
            join_space( dat1.data(), h1_slice_comps, Y.data(), Y.size(), out);
          break;
      case op_arg_magX:
          magnitude( dat1.data(), h1_slice_comps, h1->num_components(), out ); 
          break;
      case op_arg_dotXY:
          dot( dat1.data(), Yptr, h1_slice_comps, h1->num_components(), out, 
                                                                     Ytype==PT_TIME|| Ytype==PT ); 
          break;
      case op_arg_intX_i :
      case op_arg_intX_t :
          {
          int lastT    = h1->t()-1;
          int lastSimp = lastT - (lastT%3);
          if( t<=lastSimp  && lastSimp!=0 )
            simp38( dat1.data(), h1_slice_comps, t==0||t==lastSimp?-1:t, t_reduce.data(),
                                          h1->inc_t() );
          if( t>=lastSimp && lastSimp!=lastT )
            midpt( dat1.data(), h1_slice_comps, t==lastSimp||t==lastT, t_reduce.data(),
                                          h1->inc_t() );
          }
          if( opts.op_arg == op_arg_intX_t )
            fwrite( t_reduce.data(), sizeof(t_reduce[0]), t_reduce.size(), out );
          break;
      case op_arg_norm_i:
          minXi( dat1.data(), big_scratch.data(), h1_slice_comps );
          maxXi( dat1.data(), big_scratch.data()+h1_slice_comps, h1_slice_comps );
          break;
      case op_arg_maxP2P:
          maxP2P( dat1.data(), big_scratch.data(), t_reduce.data(), h1_slice_comps, t );
          break;
      case op_arg_NOP:
          break;
    }
  }

  // post process
  switch( opts.op_arg ) {
    case op_arg_avgX_i:
        for( int i=0; i<h1->slice_sz()*h1->num_components(); i++ )
          t_reduce[i] /= (data_t)(h1->t());
        break;
    case op_arg_stdD_i:
    {
      for( int i=0; i<h1->slice_sz()*h1->num_components(); i++ )
        t_reduce[i] /= (data_t)(h1->t());
      std::vector<data_t> mean;
      mean.assign(t_reduce.begin(), t_reduce.end());
      std::fill(t_reduce.begin(), t_reduce.end(), 0);
      for( int t=0; t<h1->t(); t++ ) {
        float tm = h1->org_t() + t*h1->inc_t();
        varXi( dat3[t].data(), mean.data(), t_reduce, 1.0f);
      }
      for( int i=0; i<h1->slice_sz()*h1->num_components(); i++ )
        t_reduce[i] = std::sqrt(t_reduce[i] / ((data_t)(h1->t() - 1)));
      break;
    }
    case op_arg_norm_i:
      normalize( h1, out, big_scratch.data(), big_scratch.data()+h1_slice_comps );
      break;
    case op_arg_t_maxX:
      for( int i=0; i<t_reduce.size(); i++ )
        t_reduce[i] = h1->org_t() + t_reduce[i]*h1->inc_t();
      break;
  }

  if( t_reduce.size() )
    output_scratch( t_reduce, out, true );

  return 0;
}