#include "tjvector.h"
#include "tjtypes.h"
#include "tjtest.h"
#include "tjlog_code.h"

// required for interpolation
#ifdef HAVE_LIBGSL
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#else
#include "replacements/gsl_replacement.cpp"
#endif

#ifdef HAVE_DIRENT_H
#include <dirent.h>
#endif

#ifdef USING_WIN32
#include <io.h>
#endif

const char* VectorComp::get_compName() {return "vector";}
LOGGROUNDWORK(VectorComp)


/////////////////////////////////////////////////////////////////////////////////////////////////////////////

fvector real(const cvector& cv) {
  unsigned int n=cv.size();
  fvector result(n);
  for (unsigned int i=0; i<n; i++) result[i]=cv[i].real();
  return result;
}

////////////////////////////////////////

fvector imag(const cvector& cv) {
  unsigned int n=cv.size();
  fvector result(n);
  for (unsigned int i=0; i<n; i++) result[i]=cv[i].imag();
  return result;
}

////////////////////////////////////////

fvector amplitude(const cvector& cv) {
  unsigned int n=cv.size();
  fvector result(n);
  for (unsigned int i=0; i<n; i++) result[i]=STD_abs(cv[i]);
  return result;
}

////////////////////////////////////////

fvector phase(const cvector& cv) {
  unsigned int n=cv.size();
  fvector result(n);
  for (unsigned int i=0; i<n; i++) result[i]=STD_arg(cv[i]);
  return result;
}

////////////////////////////////////////

cvector real2complex(const fvector& fv) {
  unsigned int n=fv.size();
  cvector result(n);
  for (unsigned int i=0; i<n; i++) {
    result[i]=STD_complex(fv[i]);
  }
  return result;
}

////////////////////////////////////////

dvector fvector2dvector(const fvector& fv) {
  unsigned int n=fv.size();
  dvector result(n);
  for (unsigned int i=0; i<n; i++) result[i]=fv[i];
  return result;
}

////////////////////////////////////////

fvector dvector2fvector(const dvector& dv) {
  unsigned int n=dv.size();
  fvector result(n);
  for (unsigned int i=0; i<n; i++) result[i]=dv[i];
  return result;
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////

svector tokens(const STD_string& tokenstring, char custom_separator, char escape_begin,  char escape_end) {
  Log<VectorComp> odinlog("","tokens");

  svector result;

  const int n=tokenstring.length();

  STD_string sepstr(" ");
  if(custom_separator) sepstr[0]=custom_separator;

  // fill the list with the tokens
  int beginpos=0;
  int endpos;
  STD_string onetoken;
  int beginquote=0;
  int endquote=0;
  while(beginpos>=0 && beginpos <n) {
    beginpos=textbegin(tokenstring,beginpos,custom_separator);
    endpos=sepbegin(tokenstring,beginpos,custom_separator);
    if(endpos<beginpos) endpos=n;
    if(beginpos>=0 && endpos>=0) {
      STD_string substring(tokenstring.substr(beginpos,endpos-beginpos));
      beginquote+=noccur(substring,STD_string(1,escape_begin));
      endquote+=noccur(substring,STD_string(1,escape_end));
      onetoken+=substring;
      bool token_complete=true;
      if( (escape_begin==escape_end) && (beginquote%2) ) token_complete=false;
      if( (escape_begin!=escape_end) && (beginquote>endquote) ) token_complete=false;
      if(token_complete) {
        result.push_back(onetoken);
        onetoken="";
        beginquote=endquote=0;
      } else {
        onetoken+=sepstr;
      }
    }
    beginpos=endpos;
  }

  return result;
}


////////////////////////////////////////

STD_string tokenstring(const svector& tokens,unsigned int linewidth) {
  Log<VectorComp> odinlog("","tokenstring");
  unsigned int stringsize=0;
  unsigned int n=tokens.size();
  ODINLOG(odinlog,normalDebug) << "creating string from " << n << " tokens" << STD_endl;
  ODINLOG(odinlog,normalDebug) << "linewidth = " << linewidth << STD_endl;
  unsigned int i;
  for(i=0;i<n;i++) stringsize+=tokens[i].length();
  ODINLOG(odinlog,normalDebug) << "stringsize=" << stringsize << STD_endl;
  char* tmpstring=new char[stringsize+n+100];
  char* stringptr=tmpstring;

  unsigned int width=0;
  for(i=0;i<n;i++) {
    strcpy(stringptr,tokens[i].c_str());
    unsigned int l=tokens[i].length();
    stringptr+=l;

    if(linewidth && (width>linewidth || i==(n-1))) {
      if(width>linewidth) {
        *stringptr='\n';
        width=0;
        stringptr++;
      } else *stringptr='\0';
    } else {
      if(l) {
        width+=l;
        *stringptr=' ';
        width++;
        stringptr++;
      }
    }

  }
  *stringptr='\0';
  STD_string result((const char*)tmpstring);
  delete[] tmpstring;
  return result;
}


////////////////////////////////////////

#ifndef NO_FILEHANDLING

svector browse_dir(const STD_string& dirname,bool only_dirs, bool discard_dotfiles) {
  Log<VectorComp> odinlog("","browse_dir");

  STD_list<STD_string> result;

  bool success=true;


#ifdef USING_WIN32
  // browse directory the Windows way

  _finddata_t fileinfo;

  STD_string files_match=dirname+"\\*";

  long findhandle=_findfirst(files_match.c_str(), &fileinfo );
  if(findhandle==-1) {
    ODINLOG(odinlog,normalDebug) << "(" << dirname << ") " << lasterr() << STD_endl;
    success=false;
  }


  while(findhandle!=-1) {
    STD_string curr_file(fileinfo.name);

    bool store=true;

    if(only_dirs) {
      if(fileinfo.attrib!=_A_SUBDIR) store=false;
    }
    if(discard_dotfiles) {
      if(curr_file[(unsigned int)0]=='.') store=false;
    }

    if(store) result.push_back(curr_file);

    if(_findnext(findhandle, &fileinfo )==-1) {
      _findclose(findhandle);
      findhandle=-1;
    }
  }

#else

#ifdef HAVE_DIRENT_H
  // browse directory the UNIX way

  struct dirent *ep;

  DIR* dp = opendir(dirname.c_str());

  if (dp != NULL) {
    while ( (ep = readdir (dp)) ) {

      STD_string curr_file(ep->d_name);

      bool store=true;

      if(only_dirs && ep->d_type!=DT_DIR) store=false;

      if(discard_dotfiles) {
        if(curr_file[(unsigned int)0]=='.') store=false;
      }

      if(store) result.push_back(curr_file);


    }
    closedir (dp);
  } else success=false;

#endif

#endif

  if(!success) {
    ODINLOG(odinlog,errorLog) << "Couldn't open directory >" << dirname << "< - " << lasterr() << STD_endl;
  }

  result.sort();
  return list2vector(result);
}

#endif // NO_FILEHANDLING

/////////////////////////////////////////////////////////////////////////////////////////////////////////////


STD_string svector::printbody() const {
  return tokenstring(*this,_DEFAULT_LINEWIDTH_);
/*
  STD_string result;
  for(unsigned int i=0; i<size(); i++) result+=(*this)[i]+" ";
  return result;
*/
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////

STD_complex* interpolate1D(const STD_complex* data, unsigned int oldsize, unsigned int newsize, float subpixel_shift) {
  Log<VectorComp> odinlog("tjvector","interpolate1D(complex)");
  unsigned int i;

  STD_complex* newdata=new STD_complex[newsize];
  for(i=0; i<newsize; i++) newdata[i]=0.0;

  if(oldsize==newsize && subpixel_shift==0.0) {
    for(unsigned int i=0; i<newsize; i++) newdata[i]=data[i];
    return newdata;
  }

  if(oldsize==0) {
    for(i=0; i<newsize; i++) newdata[i]=0.0;
    return newdata;
  }

  if(oldsize==1) {
    for(i=0; i<newsize; i++) newdata[i]=data[0];
    return newdata;
  }

  if(newsize==0) {
    return newdata;
  }

  if(newsize<oldsize && (oldsize%newsize)==0 && subpixel_shift==0.0) { // oldsize is a multiple of newsize -> simple adding blocks
    unsigned int blocksize=oldsize/newsize;
    ODINLOG(odinlog,normalDebug) << "blocksize=" << blocksize << STD_endl;
    for(unsigned int inew=0; inew<newsize; inew++) {
      for(unsigned int iblock=0; iblock<blocksize; iblock++) {
        newdata[inew]+=data[inew*blocksize+iblock];
      }
      newdata[inew]/=STD_complex(blocksize);
    }
    return newdata;
  }

  double* x=new double[oldsize];
  double* data_re=new double[oldsize];
  double* data_im=new double[oldsize];

  double oldstep=secureInv(double(oldsize));
  for(i=0; i<oldsize; i++) {
    x[i] = (double(i)+0.5) * oldstep;
    data_re[i]=data[i].real();
    data_im[i]=data[i].imag();
  }
  double minold=x[0]; //using doubles to prevent round-off errors
  double maxold=x[oldsize-1];
  ODINLOG(odinlog,normalDebug) << "oldstep/minold/maxold=" << oldstep << "/" << minold << "/" << maxold << STD_endl;

  gsl_interp_accel *acc_re = gsl_interp_accel_alloc ();
  gsl_interp_accel *acc_im = gsl_interp_accel_alloc ();

  const gsl_interp_type* itype=gsl_interp_linear; // fallback
//  if(oldsize>=3) itype=gsl_interp_cspline; // for few points
  if(oldsize>=5) itype=gsl_interp_akima; // this seems to be the best algorithm for non-smooth data

  gsl_spline *spline_re = gsl_spline_alloc (itype, oldsize);
  gsl_spline *spline_im = gsl_spline_alloc (itype, oldsize);

  gsl_spline_init (spline_re, x, data_re, oldsize);
  gsl_spline_init (spline_im, x, data_im, oldsize);

  for(i=0; i<newsize; i++) {
    double newx=((double(i)-subpixel_shift)+0.5)/double(newsize);
    if(newx<minold) newx=minold;
    if(newx>maxold) newx=maxold;
    float re=gsl_spline_eval (spline_re, newx, acc_re);
    float im=gsl_spline_eval (spline_im, newx, acc_im);
    newdata[i]=STD_complex(re,im);
  }

  gsl_spline_free (spline_re);
  gsl_spline_free (spline_im);
  gsl_interp_accel_free(acc_re);
  gsl_interp_accel_free(acc_im);

  delete[] x;
  delete[] data_re;
  delete[] data_im;


  // Deal with edge effects, i.e. coordinates outside the old range at the edges
  if(subpixel_shift==0.0) {
    int nedge=int(0.5*secureDivision(newsize,oldsize)+0.5);
    if(nedge>0 && nedge<int(newsize-1)) {
      ODINLOG(odinlog,normalDebug) << "newsize/oldsize/nedge=" << newsize << "/" << oldsize << "/" << nedge << STD_endl;
      for(int i=1; i<=nedge; i++) {

        // lower edge
        int iedge=nedge;
        int isrc=iedge+i;
        int idst=iedge-i;
        STD_complex edgeval2=STD_complex(2.0)*newdata[iedge];
        if(isrc>=0 && isrc<int(newsize) && idst>=0 && idst<int(newsize)) {
          newdata[idst]=edgeval2-newdata[isrc];
        }

        // upper edge
        iedge=newsize-1-nedge;
        isrc=iedge-i;
        idst=iedge+i;
        edgeval2=STD_complex(2.0)*newdata[iedge];
        if(isrc>=0 && isrc<int(newsize) && idst>=0 && idst<int(newsize)) {
          newdata[idst]=edgeval2-newdata[isrc];
        }
      }
    }
  }

  return newdata;
}


template<class T, bool integer_rounding>
T* interpolate1D_impl(const T* olddata, unsigned int oldsize, unsigned int newsize, float subpixel_shift) {
  Log<VectorComp> odinlog("tjvector","interpolate1D_impl");
  STD_complex* oldcomplexdata=new STD_complex[oldsize];
  unsigned int i;
  T minval, maxval;
  for(i=0;i<oldsize;i++) {
    if(i) {
      minval=STD_min(minval, olddata[i]);
      maxval=STD_max(maxval, olddata[i]);
    } else {
      minval=maxval=olddata[i];
    }
    oldcomplexdata[i]=olddata[i];
  }
  ODINLOG(odinlog,normalDebug) << "minval/maxval=" << minval << "/" << maxval << STD_endl;

  STD_complex* newcomplexdata=interpolate1D(oldcomplexdata,oldsize,newsize,subpixel_shift);
  T* newdata=new T[newsize];
  for(i=0;i<newsize;i++) {
    if(integer_rounding)  newdata[i]=T(newcomplexdata[i].real()+0.5);
    else                  newdata[i]=newcomplexdata[i].real();

    // keeping values in old range
    newdata[i]=STD_min(maxval, newdata[i]);
    newdata[i]=STD_max(minval, newdata[i]);
  }
  delete[] oldcomplexdata;
  delete[] newcomplexdata;
  return newdata;
}


float* interpolate1D(const float* olddata,unsigned int oldsize, unsigned int newsize, float subpixel_shift) {
  return interpolate1D_impl<float,false>(olddata, oldsize, newsize, subpixel_shift);
}

double* interpolate1D(const double* olddata,unsigned int oldsize,unsigned int newsize, float subpixel_shift) {
  return interpolate1D_impl<double,false>(olddata, oldsize, newsize, subpixel_shift);
}

int* interpolate1D(const int* olddata,unsigned int oldsize,unsigned int newsize, float subpixel_shift) {
  return interpolate1D_impl<int,true>(olddata, oldsize, newsize, subpixel_shift);
}




/**************************************************************/
/******************       Template code     *******************/
/**************************************************************/

template<class T>
tjvector<T>::tjvector(unsigned int n) : STD_vector<T>(n, T(0)) {
  Log<VectorComp> odinlog("tjvector","tjvector(unsigned int)");
  ODINLOG(odinlog,normalDebug) << "n=" << n << STD_endl;
  c_array_cache=0;
}


template<class T>
tjvector<T>::tjvector(const T *array, unsigned int n) :  STD_vector<T>() {
  STD_vector<T>::resize(n);
  set_c_array((unsigned char*)array,n);
  c_array_cache=0;
}


template<class T>
tjvector<T>::tjvector(const STD_vector<T>& v) :  STD_vector<T>(v) {
  c_array_cache=0;
}


template<class T>
tjvector<T>::tjvector(const tjvector<T>& tv) : STD_vector<T>(tv) {
  c_array_cache=0;
}

template<class T>
tjvector<T>::~tjvector() {
  Log<VectorComp> odinlog("tjvector","~tjvector");
  if(c_array_cache) {
     ODINLOG(odinlog,normalDebug) << "deleting c_array_cache=" << (void*)c_array_cache << STD_endl;
     delete[] c_array_cache;
  }
}


template<class T>
  tjvector<T>& tjvector<T>::resize(unsigned int newsize) {
    Log<VectorComp> odinlog("tjvector","resize");
    unsigned int oldsize=STD_vector<T>::size();
    if(newsize==oldsize) return *this;
    STD_vector<T> old(*this);
    STD_vector<T>::resize(newsize);
    for(unsigned int i=0; i<newsize; i++) {
      if(i<oldsize) (*this)[i]=old[i];
      else (*this)[i]=T(0);
    }
    return *this;
  }


template<class T>
unsigned int tjvector<T>::length () const {return STD_vector<T>::size();}



template<class T>
unsigned int tjvector<T>::fill_linear(const T& min,const T& max) {
  if ((length()-1)==0) tjvector<T>::operator = (min);
  else {
    T incr=(max-min)/(T)(length()-1);
    for(unsigned int i=0;i<length();i++) (*this)[i]=min+(T)(incr*(T)i);
  }
  return length();
}



template<class T>
T tjvector<T>::sum() const {
  Log<VectorComp> odinlog("","sum");
  T result=T(0);
  ODINLOG(odinlog,normalDebug) << "initial value=" << result << STD_endl;
  for(unsigned int i=0;i<length();i++) result+=(*this)[i];
  ODINLOG(odinlog,normalDebug) << "final value=" << result << STD_endl;
  return result;
}



template<class T>
tjvector<T> tjvector<T>::range(unsigned int startindex, unsigned int endindex) const {
  tjvector<T> result;
  if(endindex<=startindex) return result;
  if(endindex>length()) return result;
  if(startindex>=length()) return result;
  result.resize(endindex-startindex);
  for(unsigned int i=startindex; i<endindex; i++) result[i-startindex]=(*this)[i];
  return result;
}


template<class T>
tjvector<T>& tjvector<T>::operator = (const T& value) {
  Log<VectorComp> odinlog("tjvector","operator = (const T&)");
#ifdef STL_REPLACEMENT
  for(unsigned int i=0; i<length(); i++) (*this)[i]=value;
#else
  std::fill(STD_vector<T>::begin(),STD_vector<T>::end(),value);
#endif
  return *this;
}


template<class T>
tjvector<T>& tjvector<T>::operator = (const STD_vector<T>& vec) {
  STD_vector<T>::operator = (vec);
  return *this;
}


template<class T>
tjvector<T>& tjvector<T>::operator = (const tjvector<T>& tv) {
  Log<VectorComp> odinlog("tjvector","operator = (const tjvector<T>&)");
  STD_vector<T>::operator = (tv);
  return *this;
}

template<class T>
STD_string tjvector<T>::printbody() const {
  unsigned int n=length();
  svector tokens; tokens.resize(n);
  for(unsigned int i=0;i<n;i++) {
    tokens[i]=TypeTraits::type2string((*this)[i]);
  }
  return tokenstring(tokens,_DEFAULT_LINEWIDTH_);
}





template<class T>
T tjvector<T>::maxvalue() const {
  if(!length()) return T(0);
  T result=(*this)[0];
  for(unsigned int i=1;i<length();i++) if((*this)[i]>result) result=(*this)[i];
  return result;
}

template<class T>
T tjvector<T>::minvalue() const {
  if(!length()) return T(0);
  T result=(*this)[0];
  for(unsigned int i=1;i<length();i++) if((*this)[i]<result) result=(*this)[i];
  return result;
}


template<class T>
const T* tjvector<T>::c_array() const {
  Log<VectorComp> odinlog("tjvector","c_array");
  if(c_array_cache) {
     ODINLOG(odinlog,normalDebug) << "deleting c_array_cache=" << (void*)c_array_cache << STD_endl;
     delete[] c_array_cache;
     c_array_cache=0;
  }
  c_array_cache=new T[length()];
  ODINLOG(odinlog,normalDebug) << "allocated c_array_cache=" << (void*)c_array_cache << STD_endl;
  for(unsigned int i=0;i<length();i++) c_array_cache[i]=(*this)[i];
  return c_array_cache;
}

template<class T>
tjvector<T>& tjvector<T>::set_c_array(const unsigned char* array, unsigned int n) {
  Log<VectorComp> odinlog("tjvector","set_c_array");
  if(length()==n) {
    for(unsigned int i=0;i<n;i++) {
     (*this)[i]=((const T*)array)[i];
    }
  } else {
    ODINLOG(odinlog,errorLog) << "Size mismatch" << STD_endl;
  }
  return *this;
}
  

template<class T>
T tjvector<T>::maxabs() const{
  return T(STD_max(cabs(minvalue()),cabs(maxvalue())));
}


template<class T>
T tjvector<T>::normalize() {
  Log<VectorComp> odinlog("tjvector","normalize");
  T maxnorn=maxabs();
  ODINLOG(odinlog,normalDebug) << "maxnorn=" << maxnorn << STD_endl;
  if(!(maxnorn==T(0))) (*this)/=maxnorn;
  return maxnorn;
}


template<class T>
tjvector<T>& tjvector<T>::interpolate(unsigned int newsize, float subpixel_shift) {
  Log<VectorComp> odinlog("tjvector","interpolate");
  unsigned int oldsize=length();

  ODINLOG(odinlog,normalDebug) << "oldsize/newsize/subpixel_shift=" << oldsize << "/" << newsize << "/" << subpixel_shift << STD_endl;

  T* olddata=new T[oldsize];
  unsigned int i;
  for(i=0;i<oldsize;i++) olddata[i]=(*this)[i];
  T* newdata=interpolate1D(olddata,oldsize,newsize,subpixel_shift);
  resize(newsize);
  for(i=0;i<newsize;i++) (*this)[i]=newdata[i];
  delete[] olddata;
  delete[] newdata;
  return *this;
}




template<class T>
int tjvector<T>::load(const STD_string& fname) {
  Log<VectorComp> odinlog("tjvector","load");
#ifndef NO_FILEHANDLING
  if(fname=="") return 0;

  LONGEST_INT fsize=filesize(fname.c_str());
  LONGEST_INT nelements=fsize/sizeof(T);

  FILE* file_ptr=ODIN_FOPEN(fname.c_str(),modestring(readMode));
  if(file_ptr==NULL) {
    ODINLOG(odinlog,errorLog) << "unable to open file >" << fname << "<, " << lasterr() << STD_endl;
    return(-1);
  }

  if((LONGEST_INT)STD_vector<T>::size() != nelements) resize(nelements);

  T* buff=new T[nelements];

  if((LONGEST_INT)fread(buff,sizeof(T),nelements,file_ptr)!=nelements) {
    ODINLOG(odinlog,errorLog) << "unable to read data from file >" << fname << "<, " << lasterr() << STD_endl;
  } else {
    (*this)=tjvector<T>(buff,nelements);
  }

  if(file_ptr!=NULL) fclose(file_ptr);

  delete[] buff;
#endif
  return 0;
}



template<class T>
int tjvector<T>::write(const STD_string& fname, fopenMode mode, LONGEST_INT nelements) const {
  Log<VectorComp> odinlog("tjvector","write");
#ifndef NO_FILEHANDLING
  if(fname=="") return 0;

  if(nelements>(LONGEST_INT)(length()) || nelements<0) {
    nelements=length();
  }

  FILE* file_ptr=ODIN_FOPEN(fname.c_str(),modestring(mode));

  if(file_ptr==NULL) {
    ODINLOG(odinlog,errorLog) << "unable to create/open file >" << fname << "<, " << lasterr() << STD_endl;
    return(-1);
  }

  if((LONGEST_INT)fwrite(c_array(),sizeof(T),nelements,file_ptr)!=nelements) {
    ODINLOG(odinlog,errorLog) << "unable to write data to file >" << fname << "<, " << lasterr() << STD_endl;
  }
  if(file_ptr!=NULL) fclose(file_ptr);
#endif
  return 0;
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////////


// template instantiations
template class tjvector<float>;
template class tjvector<double>;
template class tjvector<int>;
template class STD_vector<STD_string>;
template class tjvector<STD_complex>;


#ifndef NO_UNIT_TEST
class VectorTest : public UnitTest {

 public:
  VectorTest() : UnitTest(VectorComp::get_compName()) {}

 private:

  bool check() const {
    Log<UnitTest> odinlog(this,"check");

    // Test sum()
    ivector iv1(2);
    iv1[0]=1;
    iv1[1]=3;
    ivector iv2(2);
    iv2[0]=2;
    iv2[1]=4;
    int calculated=(iv1+iv2).sum();
    int expected=10;
    if(calculated!=expected) {
      ODINLOG(odinlog,errorLog) << "sum(): calculated/expected=" << calculated << "/" << expected << STD_endl;
      return false;
    }


    // Test interpolation by converting high-res sine to and from low-res version
    int srcsize=1000;
    fvector srcvec(srcsize);
    for(int i=0; i<srcsize; i++) {
      srcvec[i]=sin(2.0*PII*float(i)/float(srcsize));
    }
//    STD_ofstream oss("srcvec.asc"); oss << srcvec.printbody();
    fvector interpvec(srcvec);
    int interpsize=20;
    interpvec.interpolate(interpsize);
    interpvec.interpolate(srcsize);
//    STD_ofstream osi("interpvec.asc"); osi << interpvec.printbody();
    float absdiff=0.0;
    for(int i=0; i<srcsize; i++) absdiff+=fabs(srcvec[i]-interpvec[i]);
    if(absdiff>0.005*srcsize) { // No more than 0.5% on average
      ODINLOG(odinlog,errorLog) << "interpolate(): absdiff=" << absdiff << STD_endl;
      return false;
    }


    // Test tokenizer
    svector gottoks=tokens("aa_bbb_<c_dd>_eee", '_', '<',  '>');
    svector exptoks; exptoks.resize(4);
    exptoks[0]="aa";
    exptoks[1]="bbb";
    exptoks[2]="<c_dd>";
    exptoks[3]="eee";
    if(gottoks!=exptoks) {
      ODINLOG(odinlog,errorLog) << "tokens: got/expected=" << gottoks.printbody() << "/" << exptoks.printbody() << STD_endl;
      return false;
    }


    return true;
  }

};

void alloc_VectorTest() {new VectorTest();} // create test instance
#endif