GridData1.h

Go to the documentation of this file.

// -*- C++ -*-

#ifndef _included_GridData_1_h
#define _included_GridData_1_h

#include "DAGHParams.h"

#ifdef DEBUG_PRINT_GD_MEMORY
#include "DAGHMemoryTrace.h"
#endif

#include "BBox.h"
#include "Coords.h"
#include "PackedGridDataBucket.h"

#ifdef DEBUG_PRINT
#include "CommServer.h"
#endif

#include "IndexGridData1.h"
#include "generic.h"

#ifndef gd_OperateRegion
#define gd_OperateRegion(op)    name2(gd_OperateRegion,op)
#endif

// #include "GDIterator.h"

#include <iosfwd>
#include <iostream>
#include <fstream>
#include <sstream>
#include <cstdlib>
#include <cassert>

template <class Type, int dim> class GridData;
template <class Type>
std::ostream& operator<<(std::ostream&, const GridData<Type,1>&);
template <class Type>
std::ofstream& operator<<(std::ofstream&, const GridData<Type,1>&);
template <class Type>
std::ifstream& operator>>(std::ifstream&, GridData<Type,1>&);
template <class Type>
std::stringstream& operator<<(std::stringstream&, const GridData<Type,1>&);
template <class Type>
std::stringstream& operator>>(std::stringstream&, GridData<Type,1>&);

template <class Type> class GridData<Type,1> { 
  BBox _bbox;
  Coords _extents;
  Coords _step;
  int _size;
  int _bottom;
  Type *_data;

public:
  GridData(void) : _bbox(), _extents(0,0), _step(1,1),
                   _size(0), _bottom(0), _data((Type *) 0) {}

  GridData(const BBox& bb)  : _bbox(bb) {
    _bbox.rank = 1; 
    _bbox.lower().rank = 1;
    _bbox.upper().rank = 1;
    _bbox.stepsize().rank = 1;
    _extents = _bbox.extents();
    _step = _bbox.stepsize();
    _bottom = _bbox.bottom();
    _size = _bbox.size();
    _data = _size ? new Type[_size] : ((Type *) 0);
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size); 
#endif
  }

  GridData(const int i, const int ii)  : 
    _bbox(1,i,ii,1), _extents(_bbox.extents()), _step(_bbox.stepsize()),
    _size(_bbox.size()), _bottom(_bbox.bottom()),
    _data(_size ? new Type[_size] : ((Type *) 0)) {
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size); 
#endif
  }

  GridData(const int i, const int ii,  const int s)  : 
    _bbox(1,i,ii,s), _extents(_bbox.extents()), _step(_bbox.stepsize()),
    _size(_bbox.size()), _bottom(_bbox.bottom()),
    _data(_size ? new Type[_size] : ((Type *) 0)) {
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size); 
#endif
  }

  GridData(const BBox& bb, Type *databuf) : 
    _bbox(bb), _data(databuf) {
    _bbox.rank = 1; 
    _bbox.lower().rank = 1;
    _bbox.upper().rank = 1;
    _bbox.stepsize().rank = 1;
    _extents = _bbox.extents();
    _step = _bbox.stepsize();
    _bottom = _bbox.bottom();
    _size = _bbox.size();
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size); 
#endif
  }

  GridData(GridDataBucket<Type> &gdbkt) : 
    _bbox(gdbkt.bbox()), _data(gdbkt.data()) {
    _bbox.rank = 1; 
    _bbox.lower().rank = 1;
    _bbox.upper().rank = 1;
    _bbox.stepsize().rank = 1;
    _extents = _bbox.extents();
    _step = _bbox.stepsize();
    _bottom = _bbox.bottom();
    _size = _bbox.size();
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size); 
#endif
  }

   GridData(GridDataBucket<Type> &gdbkt, const int n): _bbox(gdbkt.bbox(n)),
    _data(gdbkt.data(n))
    {
      _bbox.rank = 1; 
      _bbox.lower().rank = 1;
      _bbox.upper().rank = 1;
      _bbox.stepsize().rank = 1;
      _extents = _bbox.extents();
      _step = _bbox.stepsize();
      _bottom = _bbox.bottom();
      _size = _bbox.size();
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
                          << _bbox 
                          << ": " << _data[_size-1]
                          << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size); 
#endif
    }

  /* A "pseudo" copy-constructor */
  GridData (const GridData<Type,1> &other) : 
    _bbox(other._bbox), _extents(other._extents), _step(other._step),
    _size(other._size), _bottom(other._bottom),
    _data(_size ? new Type[_size] : ((Type *) 0)) {
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size); 
#endif
  }

   inline ~GridData(void) { 
     if (_data) {
       delete [] _data; 
#ifdef DEBUG_PRINT_GD_MEMORY
       DAGHMemoryTrace::free(sizeof(Type)*_size);
#endif
     }
   }

  void allocate(const BBox &bb) {
    _bbox = bb;
    _bbox.rank = 1; 
    _bbox.lower().rank = 1;
    _bbox.upper().rank = 1;
    _bbox.stepsize().rank = 1;
    _extents = _bbox.extents();
    _step = _bbox.stepsize();
    _bottom = _bbox.bottom();
    _size = _bbox.size();
    _data = _size ? new Type[_size] : (Type *) 0;
#ifdef DEBUG_PRINT_GD
    assert (_data!=0);
    fill(9999);
    ( comm_service::log() << "Data Check"
      << _bbox 
      << ": " << _data[_size-1]
      << std::endl ).flush();
#endif
#ifdef DEBUG_PRINT_GD_MEMORY
    DAGHMemoryTrace::alloc(sizeof(Type)*_size); 
#endif
  }

   void allocate(const BBox &bb, Type *databuf) {
     _bbox = bb;
     _bbox.rank = 1; 
     _bbox.lower().rank = 1;
     _bbox.upper().rank = 1;
     _bbox.stepsize().rank = 1;
     _extents = _bbox.extents();
     _step = _bbox.stepsize();
     _bottom = _bbox.bottom();
     _size = _bbox.size();
#ifdef DEBUG_PRINT_GD
     assert(_data==0);
#endif
     _data = databuf;
#ifdef DEBUG_PRINT_GD
     assert (_data!=0);
     fill(9999);
     ( comm_service::log() << "Data Check"
       << _bbox 
       << ": " << _data[_size-1]
       << std::endl ).flush();
#endif
   }

   inline void allocate(Type *databuf) {
     if (_data) {
       delete [] _data; _data = (Type*)0;
#ifdef DEBUG_PRINT_GD_MEMORY
       DAGHMemoryTrace::free(sizeof(Type)*_size);
#endif
     }
     _data = databuf;
   }
  
  inline void deallocate() {
    if (_data) {
      delete [] _data; _data = (Type*)0;
#ifdef DEBUG_PRINT_GD_MEMORY
      DAGHMemoryTrace::free(sizeof(Type)*_size);
#endif
    }
  }
  inline void deallocate(Type*& databuf)
  { databuf = _data; _data = (Type*) 0; }
  
  void* databuffer() { return ((void*) _data); }

  /* Query funtions */
  /* Define some inline functions for getting at region extents */
  
  inline const Coords& lower() const { return(_bbox.lower()); }
  inline const Coords& upper() const { return(_bbox.upper()); }
  inline const Coords& extents() const { return(_extents); }
  inline const Coords& stepsize() const { return(_step); }
  inline Coords lower() { return(_bbox.lower()); }
  inline Coords upper() { return(_bbox.upper()); }
  inline Coords extents() { return(_extents); }
  inline Coords stepsize() { return(_step); }
  inline int bottom() const { return(_bottom); }
  inline int lower(const int i) const { return(_bbox.lower(i)); }
  inline int upper(const int i) const { return(_bbox.upper(i)); }
  inline int extents(const int i) const { return(_extents(i)); }
  inline int stepsize(const int i) const { return(_step(i)); }
  inline const BBox& bbox() const { return(_bbox); }

  /* Define the inline functions for indexing the GridData object */

  inline int ok_to_index() { return(_data != (Type *) 0); }
  inline int size() const { return(_size); }
  
  inline int idx(const int i) const
  { return( _bottom+(i/_step(0)) ); }
  
  inline const Type& operator () (const int i) const
  { assert(idx(i) >= 0 && idx(i) < _size); return(_data[idx(i)]); }
  inline Type& operator () (const int i)
  { assert(idx(i) >= 0 && idx(i) < _size); return(_data[idx(i)]); }
  
  inline const Type& operator () (const Coords& c) const
  { assert(idx(c(0)) >= 0 && idx(c(0)) < _size); return(_data[idx(c(0))]); }
  inline Type& operator () (const Coords& c)
  { assert(idx(c(0)) >= 0 && idx(c(0)) < _size); return(_data[idx(c(0))]); }
  
  inline const Type *ptr(const int i) const
  { assert(idx(i) >= 0 && idx(i) < _size); return(_data + idx(i)); }
  inline Type *ptr(const int i)
  { assert(idx(i) >= 0 && idx(i) < _size); return(_data + idx(i)); }
  
  inline const Type *ptr(const Coords& c) const
  { assert(idx(c(0)) >= 0 && idx(c(0)) < _size); return(_data + idx(c(0))); }
  inline Type *ptr(const Coords& c)
  { assert(idx(c(0)) >= 0 && idx(c(0)) < _size); return(_data + idx(c(0))); }
  
  inline const Type *data() const { return(_data); }
  inline Type *data() { return(_data); }

   void fill(const Type &val)
  { register int i; for (i = 0; i < _size; i++) _data[i] = val; }

  /* Define the copy functions for communication between objects */
  void copy(const GridData<Type,1> &gd) {
    if (&gd != this) {
      BBox intersection = _bbox * gd._bbox;
      if (!intersection.empty()) { 
        Coords max_step = _step.getmax(gd._step);
        BBox to(intersection); to.setstepsize(_step); 
        BBox from(intersection); from.setstepsize(gd._step);
        gd_CopyRegion(gd, to, from, max_step); 
      }
    }
  }

  void copy(const GridData<Type,1> &gd, const BBox &where) {
    if (&gd != this) {
      BBox intersection = _bbox * gd._bbox * where;
      if (!intersection.empty()) { 
        Coords max_step = _step.getmax(gd._step);
        BBox to(intersection); to.setstepsize(_step); 
        BBox from(intersection); from.setstepsize(gd._step);
         gd_CopyRegion(gd, to, from, max_step); 
      }
    }
  }

  void copy(const GridData<Type,1> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_CopyRegion(gd, newto, newfrom, max_step);
    }
  }

  void copy(const GridDataBucket<Type> &gdbkt)  {
    BBox const &gdbktbb = gdbkt.bbox();
    BBox intersection = _bbox * gdbktbb;
    if (!intersection.empty()) { 
      Coords max_step = _step.getmax(gdbktbb.stepsize());
      BBox to(intersection); to.setstepsize(_step); 
      BBox from(intersection); from.setstepsize(gdbktbb.stepsize());
      gdb_CopyRegion(gdbkt, to, from, max_step); 
    }
  }

  void copy(const GridDataBucket<Type> &gdbkt, const BBox &where) {
    BBox const &gdbktbb = gdbkt.bbox();
    BBox intersection = _bbox * gdbktbb * where;
    if (!intersection.empty()) { 
      Coords max_step = _step.getmax(gdbktbb.stepsize());
      BBox to(intersection); to.setstepsize(_step); 
      BBox from(intersection); from.setstepsize(gdbktbb.stepsize());
      gdb_CopyRegion(gdbkt, to, from, max_step); 
    }
  }

  void copy(const GridDataBucket<Type> &gdbkt, const BBox &to, 
            const BBox &from) {
    BBox const &gdbktbb = gdbkt.bbox();
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gdbktbb * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gdbktbb.stepsize());
      newto.setstepsize(_step);
      newfrom.setstepsize(gdbktbb.stepsize());
      gdb_CopyRegion(gdbkt, newto, newfrom, max_step);
    }
  }
  
  void copy(const GridDataBucket<Type> &gdbkt, const int n) {
    BBox const &gdbktbb = gdbkt.bbox(n);
    BBox intersection = _bbox * gdbktbb;
    if (!intersection.empty()) { 
      Coords max_step = _step.getmax(gdbktbb.stepsize());
      BBox to(intersection); to.setstepsize(_step); 
      BBox from(intersection); from.setstepsize(gdbktbb.stepsize());
      gdb_CopyRegion(gdbkt, n, to, from, max_step); 
    }
  }
  
  void copy(const GridDataBucket<Type> &gdbkt, const int n, const BBox &where) {
    BBox const &gdbktbb = gdbkt.bbox(n);
    BBox intersection = _bbox * gdbktbb * where;
    if (!intersection.empty()) { 
      Coords max_step = _step.getmax(gdbktbb.stepsize());
      BBox to(intersection); to.setstepsize(_step); 
      BBox from(intersection); from.setstepsize(gdbktbb.stepsize());
      gdb_CopyRegion(gdbkt, n, to, from, max_step); 
    }
  }

  void copy(const GridDataBucket<Type> &gdbkt, const int n, const BBox &to, 
            const BBox &from) {
    BBox const &gdbktbb = gdbkt.bbox(n);
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gdbktbb * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gdbktbb.stepsize());
      newto.setstepsize(_step);
      newfrom.setstepsize(gdbktbb.stepsize());
      gdb_CopyRegion(gdbkt, n, newto, newfrom, max_step);
    }
  }

  /* Linear interpolation */
  void lin_interp(const GridData<Type,1> &gd1, const double frac1,
                  const GridData<Type,1> &gd2, const double frac2,
                  const BBox &where) {
    BBox intersection = _bbox * gd1._bbox * gd2._bbox * where;
    if (!intersection.empty()) {

      Coords max_step = _step.getmax(gd1._step);
      max_step.max(gd2._step);

      GridData<Type,1> &dst = *this;

      BeginFastIndex1(src1, gd1._bbox, gd1._data, const Type);
      BeginFastIndex1(src2, gd2._bbox, gd2._data, const Type);
      BeginFastIndex1(dst, dst._bbox, dst._data, Type);

      for_1(i, intersection, max_step)
        FastIndex1(dst,i) = (Type) ((FastIndex1(src1,i) * frac1)
                                    + (FastIndex1(src2,i) * frac2));
      end_for

      EndFastIndex1(dst);
      EndFastIndex1(src1);
      EndFastIndex1(src2);
    }
  }
    
  inline void lin_interp(const GridData<Type,1> &gd1, const double frac1,
                         const GridData<Type,1> &gd2, const double frac2)
  { GridData<Type,1>::lin_interp(gd1,frac1,gd2,frac2,_bbox); }

  /* first moment (i think) sum(q(i)*i) */
  Type moment1(const int axis, const BBox &where) {
    Type m1 = (Type) 0;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);
      
      GridData<Type,1> &dst = *this;

      BeginFastIndex1(dst, dst._bbox, dst._data, Type);

      if (axis == DAGH_X) {
        for_1(i, intersection, max_step)
          m1 += FastIndex1(dst,i) * i;
        end_for
      }

      EndFastIndex1(dst);
    }
    return (m1);
  }

  inline Type moment1(const int axis)
  { return (GridData<Type,1>::moment1(axis,_bbox)); }

  /* sum of squares - for calculating L2-norm */
  double sumsqrd(const BBox &where) {
    double s = 0.0;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);
      
      GridData<Type,1> &dst = *this;
      
      BeginFastIndex1(dst, dst._bbox, dst._data, Type);

      for_1(i, intersection, max_step)
        s += Abs(FastIndex1(dst,i)) * Abs(FastIndex1(dst,i));
      end_for

      EndFastIndex1(dst);

    }
    return s;
  }

  inline double sumsqrd()
  { return (GridData<Type,1>::sumsqrd(_bbox)); }

  /* sum of absolute values - for calculating L1-norm */
  double sumabs(const BBox &where) {
    double s = 0.0;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);
      
      GridData<Type,1> &dst = *this;
      
      BeginFastIndex1(dst, dst._bbox, dst._data, Type);

      for_1(i, intersection, max_step)
        s += Abs(FastIndex1(dst,i));
      end_for

      EndFastIndex1(dst);

    }
    return s;
  }

  inline double sumabs()
  { return (GridData<Type,1>::sumabs(_bbox)); }

  /* Maximum of absolute values - for calculating Loo-norm */
  double maxabs(const BBox &where) {
    double m = 0.0;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);
      
      GridData<Type,1> &dst = *this;

      BeginFastIndex1(dst, dst._bbox, dst._data, Type);

      for_1(i, intersection, max_step)
        double t = Abs(FastIndex1(dst,i));
        m = Max(m,t);
      end_for

      EndFastIndex1(dst);

    }
    return m;
  }

  inline double maxabs()
  { return (GridData<Type,1>::maxabs(_bbox)); }

  /* operator == */
  BBox is_eq(const Type &val, const BBox &where)  {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      BBox bb(1,max_step);
      gd_OperateRegion(is_eq)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(1,1); return (bb); }
  }

  inline BBox is_eq(const Type &val) 
  { return (GridData<Type,1>::is_eq(val, _bbox)); }

  BBox is_eq(const GridData<Type,1> &gd, const BBox &to, const BBox &from)  {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      BBox bb(1,max_step);
      gd_OperateRegion(is_eq)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(1,1); return (bb); }
  }

  inline BBox is_eq(const GridData<Type,1> &gd)
  { return (GridData<Type,1>::is_eq(gd, gd._bbox, gd._bbox)); }
  inline BBox is_eq(const GridData<Type,1> &gd, const BBox &where)
  { return (GridData<Type,1>::is_eq(gd, where, where)); }

  /* operator != */
  BBox is_neq(const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      BBox bb(1,max_step);
      gd_OperateRegion(is_neq)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(1,1); return (bb); }
  }

  inline BBox is_neq(const Type &val) 
  { return (GridData<Type,1>::is_neq(val, _bbox)); }

  BBox is_neq(const GridData<Type,1> &gd, const BBox &to, const BBox &from)  {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      BBox bb(1,max_step);
      gd_OperateRegion(is_neq)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(1,1); return (bb); }
  }

  inline BBox is_neq(const GridData<Type,1> &gd)
  { return (GridData<Type,1>::is_neq(gd, gd._bbox, gd._bbox)); }
  inline BBox is_neq(const GridData<Type,1> &gd, const BBox &where)
  { return (GridData<Type,1>::is_neq(gd, where, where)); }

  /* operator > */
  BBox is_gt(const Type &val, const BBox &where)  {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      BBox bb(1,max_step);
      gd_OperateRegion(is_gt)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(1,1); return (bb); }
  }

  inline BBox is_gt(const Type &val) 
  { return (GridData<Type,1>::is_gt(val, _bbox)); }
  
  BBox is_gt(const GridData<Type,1> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      BBox bb(1,max_step);
      gd_OperateRegion(is_gt)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(1,1); return (bb); }
  }

  inline BBox is_gt(const GridData<Type,1> &gd)
  { return (GridData<Type,1>::is_gt(gd, gd._bbox, gd._bbox)); }
  inline BBox is_gt(const GridData<Type,1> &gd, const BBox &where)
  { return (GridData<Type,1>::is_gt(gd, where, where)); }

  /* operator >= */
  BBox is_ge(const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      BBox bb(1,max_step);
      gd_OperateRegion(is_ge)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(1,1); return (bb); }
  }

  inline BBox is_ge(const Type &val) 
  { return (GridData<Type,1>::is_ge(val, _bbox)); }

  BBox is_ge(const GridData<Type,1> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      BBox bb(1,max_step);
      gd_OperateRegion(is_ge)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
     { BBox bb(1,1); return (bb); }
  }

  inline BBox is_ge(const GridData<Type,1> &gd)
  { return (GridData<Type,1>::is_ge(gd, gd._bbox, gd._bbox)); }
  inline BBox is_ge(const GridData<Type,1> &gd, const BBox &where)
  { return (GridData<Type,1>::is_ge(gd, where, where)); }

  /* operator < */
  BBox is_lt(const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      BBox bb(1,max_step);
      gd_OperateRegion(is_lt)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(1,1); return (bb); }
  }

  inline BBox is_lt(const Type &val) 
  { return (GridData<Type,1>::is_lt(val, _bbox)); }

  BBox is_lt(const GridData<Type,1> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      BBox bb(1,max_step);
      gd_OperateRegion(is_lt)(gd, newto, newfrom, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(1,1); return (bb); }
  }

  inline BBox is_lt(const GridData<Type,1> &gd)
  { return (GridData<Type,1>::is_lt(gd, gd._bbox, gd._bbox)); }
  inline BBox is_lt(const GridData<Type,1> &gd, const BBox &where)
  { return (GridData<Type,1>::is_lt(gd, where, where)); }

  /* operator <= */
  BBox is_le(const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      BBox bb(1,max_step);
      gd_OperateRegion(is_le)(val, to, max_step, bb);
      return (bb);
    }
    else 
      { BBox bb(1,1); return (bb); }
  }

  inline BBox is_le(const Type &val) 
  { return (GridData<Type,1>::is_le(val, _bbox)); }

  BBox is_le(const GridData<Type,1> &gd, const BBox &to, const BBox &from) {
   const Coords  toshift = from.lower() - to.lower();
   BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
   if (!newfrom.empty()) {
     BBox newto = shiftabs(newfrom, -toshift);
     Coords max_step = _step.getmax(gd._step);
     newto.setstepsize(_step);
     newfrom.setstepsize(gd._step);
     BBox bb(1,max_step);
     gd_OperateRegion(is_le)(gd, newto, newfrom, max_step, bb);
     return (bb);
   }
   else 
     { BBox bb(1,1); return (bb); }
  }

  inline BBox is_le(const GridData<Type,1> &gd)
  { return (GridData<Type,1>::is_le(gd, gd._bbox, gd._bbox)); }
  inline BBox is_le(const GridData<Type,1> &gd, const BBox &where)
  { return (GridData<Type,1>::is_le(gd, where, where)); }

  /* operator = */
  void equals(const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(equal)(val, to, max_step);
    }
  }

  inline void equals(const Type &val) 
  { GridData<Type,1>::equals(val, _bbox); }

  void equals(const GridData<Type,1> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(equal)(gd, newto, newfrom, max_step);
    }
  }
    
  inline void equals(const GridData<Type,1> &gd)
  { GridData<Type,1>::equals(gd, gd._bbox, gd._bbox); }
  inline void equals(const GridData<Type,1> &gd, const BBox &where)
  { GridData<Type,1>::equals(gd, where, where); }

  /* operator += */
  void plus (const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(plus)(val, to, max_step);
    }
  }
    
  inline void plus (const Type &val) 
  { GridData<Type,1>::plus(val, _bbox); }

  void plus (const GridData<Type,1> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(plus)(gd, newto, newfrom, max_step);
    }
  }

  inline void plus (const GridData<Type,1> &gd)
  { GridData<Type,1>::plus(gd, gd._bbox, gd._bbox); }
  inline void plus (const GridData<Type,1> &gd, const BBox &where)
  { GridData<Type,1>::plus(gd, where, where); }

  /* operator -= */
  void minus (const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(minus)(val, to, max_step);
    }
  }

  inline void minus (const Type &val) 
  { GridData<Type,1>::minus(val, _bbox); }

  void minus (const GridData<Type,1> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(minus)(gd, newto, newfrom, max_step);
    }
  }

  inline void minus (const GridData<Type,1> &gd)
  { GridData<Type,1>::minus(gd, gd._bbox, gd._bbox); }
  inline void minus (const GridData<Type,1> &gd, const BBox &where)
  { GridData<Type,1>::minus(gd, where, where); }

  /* operator *= */
  void multiply(const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(mult)(val, to, max_step);
    }
  }

  inline void multiply(const Type &val) 
  { GridData<Type,1>::multiply(val, _bbox); }

  void multiply (const GridData<Type,1> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(mult)(gd, newto, newfrom, max_step);
    }
  }

  inline void multiply (const GridData<Type,1> &gd)
  { GridData<Type,1>::multiply(gd, gd._bbox, gd._bbox); }
  inline void multiply (const GridData<Type,1> &gd, const BBox &where)
  { GridData<Type,1>::multiply(gd, where, where); }

  /* operator /= */
  void divide (const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(div)(val, to, max_step);
    }
  }

  inline void divide (const Type &val) 
  { GridData<Type,1>::divide(val, _bbox); }

  void divide (const GridData<Type,1> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(div)(gd, newto, newfrom, max_step);
    }
  }

  inline void divide (const GridData<Type,1> &gd)
  { GridData<Type,1>::divide(gd, gd._bbox, gd._bbox); }
  inline void divide (const GridData<Type,1> &gd, const BBox &where)
  { GridData<Type,1>::divide(gd, where, where); }

  /* point-wise minimum */
  void minimum (const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(min)(val, to, max_step);
    }
  }

  inline void minimum (const Type &val) 
  { GridData<Type,1>::minimum(val, _bbox); }

  void minimum (const GridData<Type,1> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(min)(gd, newto, newfrom, max_step);
    }
  }

  inline void minimum (const GridData<Type,1> &gd)
  { GridData<Type,1>::minimum(gd, gd._bbox, gd._bbox); }
  inline void minimum (const GridData<Type,1> &gd, const BBox &where)
  { GridData<Type,1>::minimum(gd, where, where); }

  /* point-wise maximum */
  void maximum (const Type &val, const BBox &where) {
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      BBox to(intersection); to.setstepsize(_step);
      gd_OperateRegion(max)(val, to, max_step);
    }
  }

  inline void maximum (const Type &val) 
  { GridData<Type,1>::maximum(val, _bbox); }

  void maximum (const GridData<Type,1> &gd, const BBox &to, const BBox &from) {
    const Coords  toshift = from.lower() - to.lower();
    BBox newfrom = gd._bbox * from * shiftabs(_bbox * to, toshift);
    if (!newfrom.empty()) {
      BBox newto = shiftabs(newfrom, -toshift);
      Coords max_step = _step.getmax(gd._step);
      newto.setstepsize(_step);
      newfrom.setstepsize(gd._step);
      gd_OperateRegion(max)(gd, newto, newfrom, max_step);
    }
  }

  inline void maximum (const GridData<Type,1> &gd)
  { GridData<Type,1>::maximum(gd, gd._bbox, gd._bbox); }
  inline void maximum (const GridData<Type,1> &gd, const BBox &where)
  { GridData<Type,1>::maximum(gd, where, where); }

  /**********************************************************************/
  inline BBox operator == (const Type &val)
  { return (GridData<Type,1>::is_eq(val, _bbox)); }
  inline BBox operator == (const GridData<Type,1> &gd)
  { return (GridData<Type,1>::is_eq(gd, gd._bbox, gd._bbox)); }
  
  inline BBox operator != (const Type &val)
  { return (GridData<Type,1>::is_neq(val, _bbox)); }
  inline BBox operator != (const GridData<Type,1> &gd)
  { return (GridData<Type,1>::is_neq(gd, gd._bbox, gd._bbox)); }
  
  inline BBox operator > (const Type &val)
  { return (GridData<Type,1>::is_gt(val, _bbox)); }
  inline BBox operator > (const GridData<Type,1> &gd)
  { return (GridData<Type,1>::is_gt(gd, gd._bbox, gd._bbox)); }

  inline BBox operator >= (const Type &val)
  { return (GridData<Type,1>::is_ge(val, _bbox)); }
  inline BBox operator >= (const GridData<Type,1> &gd)
  { return (GridData<Type,1>::is_ge(gd, gd._bbox, gd._bbox)); }

  inline BBox operator < (const Type &val)
  { return (GridData<Type,1>::is_lt(val, _bbox)); }
  inline BBox operator < (const GridData<Type,1> &gd)
  { return (GridData<Type,1>::is_lt(gd, gd._bbox, gd._bbox)); }
  
  inline BBox operator <= (const Type &val)
  { return (GridData<Type,1>::is_le(val, _bbox)); }
  inline BBox operator <= (const GridData<Type,1> &gd)
  { return (GridData<Type,1>::is_le(gd, gd._bbox, gd._bbox)); }
  
  inline void operator = (const Type &val)
  { GridData<Type,1>::equals(val, _bbox); }
  inline void operator = (const GridData<Type,1> &gd)
  { GridData<Type,1>::equals(gd, gd._bbox, gd._bbox); }
  
  inline void operator += (const Type &val)
  { GridData<Type,1>::plus(val, _bbox); }
  inline void operator += (const GridData<Type,1> &gd)
  { GridData<Type,1>::plus(gd, gd._bbox, gd._bbox); }
  
  inline void operator -= (const Type &val)
  { GridData<Type,1>::minus(val, _bbox); }
  inline void operator -= (const GridData<Type,1> &gd)
  { GridData<Type,1>::minus(gd, gd._bbox, gd._bbox); }
  
  inline void operator *= (const Type &val)
  { GridData<Type,1>::multiply(val, _bbox); }
  inline void operator *= (const GridData<Type,1> &gd)
  { GridData<Type,1>::multiply(gd, gd._bbox, gd._bbox); }
  
  inline void operator /= (const Type &val)
  { GridData<Type,1>::divide(val, _bbox); }
  inline void operator /= (const GridData<Type,1> &gd)
  { GridData<Type,1>::divide(gd, gd._bbox, gd._bbox); }
  /**********************************************************************/

  /**********************************************************************/
  // Reduction Operations
  /**********************************************************************/
  
  /****************************** maxval ***************************************/
  Type maxval (const BBox &where) {
    Type max_val = (Type) 0;  
    short flag = DAGHTrue;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);
      
      GridData<Type,1> &dst = *this;
      
      BeginFastIndex1(dst, dst._bbox, dst._data, Type);
      
      for_1(i, intersection, max_step)
        if (flag == DAGHTrue) 
          { max_val = FastIndex1(dst,i); flag = DAGHFalse; }
        else 
          max_val = Max(max_val,FastIndex1(dst,i));
      end_for

      EndFastIndex1(dst);

    }
    return (max_val);
  }

  inline Type maxval (void)
  { return (GridData<Type,1>::maxval(_bbox)); }
  /************************************************************************/

  /****************************** minval ***************************************/
  Type minval (const BBox &where) {
    Type min_val = (Type) 0;
    short flag = DAGHTrue;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);
      
      GridData<Type,1> &dst = *this;
      
      BeginFastIndex1(dst, dst._bbox, dst._data, Type);

      for_1(i, intersection, max_step)
        if (flag == DAGHTrue) 
          { min_val = FastIndex1(dst,i); flag = DAGHFalse; }
        else 
          min_val = Min(min_val,FastIndex1(dst,i));
      end_for

      EndFastIndex1(dst);

    }
    return (min_val);
  }

  inline Type minval (void)
  { return (GridData<Type,1>::minval(_bbox)); }
  /************************************************************************/

  /****************************** sum ***************************************/
  Type sum (const BBox &where) {
    Type sum_val = (Type) 0;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);
      
      GridData<Type,1> &dst = *this;

      BeginFastIndex1(dst, dst._bbox, dst._data, Type);

      for_1(i, intersection, max_step)
        sum_val += FastIndex1(dst,i);
      end_for

      EndFastIndex1(dst);

    }
    return (sum_val);
  }

  inline Type sum (void)
  { return (GridData<Type,1>::sum(_bbox)); }
  /************************************************************************/
  
  /****************************** product ***************************************/
  Type product (const BBox &where) {
    Type prod_val = (Type) 1;
    BBox intersection = _bbox * where;
    if (!intersection.empty()) {
      Coords max_step = _step.getmax(where.stepsize());
      intersection.setstepsize(_step);

      GridData<Type,1> &dst = *this;

      BeginFastIndex1(dst, dst._bbox, dst._data, Type);

      for_1(i, intersection, max_step)
        prod_val *= FastIndex1(dst,i);
      end_for

      EndFastIndex1(dst);

    }
    return (prod_val);
  }

  inline Type product (void)
  { return (GridData<Type,1>::product(_bbox)); }
  /**********************************************************************/

private:
   void gd_CopyRegion(const GridData<Type,1> &src, const BBox &to,
                      const BBox &from, const Coords &step) {
     GridData<Type,1> &dst = *this;
     
     if (dst._bbox==src._bbox && to==from && to==dst._bbox) 
       std::memcpy((void *) dst._data, (void *) src._data, sizeof(Type)*to.size());
     else {
       const int di = from.lower(0)-to.lower(0);
       
       BeginFastIndex1(src, src._bbox, src._data, const Type);
       BeginFastIndex1(dst, dst._bbox, dst._data, Type);
       
       for_1(i, to, step)
         FastIndex1(dst,i) = FastIndex1(src,i+di);
       end_for

       EndFastIndex1(dst);
       EndFastIndex1(src);
     }
   }

   void gdb_CopyRegion(const GridDataBucket<Type> &gdbkt, const BBox &to,
                       const BBox &from, const Coords &step) {
     GridData<Type,1> &dst = *this;

     if (dst._bbox==gdbkt.bbox() && to==from && to==dst._bbox) 
       std::memcpy((void *) dst._data, (void *) gdbkt.data(), sizeof(Type)*to.size());
     else {
       const int di = from.lower(0)-to.lower(0);

       BeginFastIndex1(src, gdbkt.bbox(), gdbkt.data(), const Type);
       BeginFastIndex1(dst, dst.bbox(), dst.data(), Type);
       for_1(i, to, step)
         FastIndex1(dst,i) = FastIndex1(src,i+di);
       end_for
       EndFastIndex1(dst);
       EndFastIndex1(src);
     }
   }

  void gdb_CopyRegion(const GridDataBucket<Type> &gdbkt, const int n, const BBox &to,
                      const BBox &from, const Coords &step) {
    GridData<Type,1> &dst = *this;

    if (dst._bbox==gdbkt.bbox(n) && to==from && to==dst._bbox) 
      std::memcpy((void *) dst._data, (void *) gdbkt.data(n), sizeof(Type)*to.size());
    else {
      const int di = from.lower(0)-to.lower(0);
      
      BeginFastIndex1(src, gdbkt.bbox(n), gdbkt.data(n), const Type);
      BeginFastIndex1(dst, dst.bbox(), dst.data(), Type);
      for_1(i, to, step)
        FastIndex1(dst,i) = FastIndex1(src,i+di);
      end_for
      EndFastIndex1(dst);
      EndFastIndex1(src);
    }
  }

  /****************************** == ***************************************/
  void gd_OperateRegion(is_eq)(const Type &val, const BBox &to,
                               const Coords &step, BBox &bb) {
    GridData<Type,1> &dst = *this;
    
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);
   
    for_1(i, to, step)
      if (FastIndex1(dst,i) == val) bb += Coords(1,i);
    end_for

    EndFastIndex1(dst);
  }
    
  void gd_OperateRegion(is_eq)(const GridData<Type,1> &src, const BBox &to,
                               const BBox &from, const Coords &step, BBox &bb) {
    GridData<Type,1> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    
    BeginFastIndex1(src, src._bbox, src._data, const Type);
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      if (FastIndex1(dst,i) == FastIndex1(src,i+di))
        bb += Coords(1,i);
    end_for

    EndFastIndex1(dst);
    EndFastIndex1(src);
  }
  /************************************************************************/
  
  /****************************** != ***************************************/
  void gd_OperateRegion(is_neq)(const Type &val, const BBox &to,
                                const Coords &step, BBox &bb) {
    GridData<Type,1> &dst = *this;

    BeginFastIndex1(dst, dst._bbox, dst._data, Type);
    
    for_1(i, to, step)
     if (FastIndex1(dst,i) != val) bb += Coords(1,i);
    end_for

    EndFastIndex1(dst);
  }

  void gd_OperateRegion(is_neq)(const GridData<Type,1> &src, const BBox &to,
                                const BBox &from, const Coords &step, BBox &bb)  {
    GridData<Type,1> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    
    BeginFastIndex1(src, src._bbox, src._data, const Type);
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      if (FastIndex1(dst,i) != FastIndex1(src,i+di))
        bb += Coords(1,i);
    end_for

    EndFastIndex1(dst);
    EndFastIndex1(src);
  }
  /************************************************************************/

  /****************************** > ***************************************/
  void gd_OperateRegion(is_gt)(const Type &val, const BBox &to,
                               const Coords &step, BBox &bb)  {
    GridData<Type,1> &dst = *this;

    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      if (FastIndex1(dst,i) > val) bb += Coords(1,i);
    end_for

    EndFastIndex1(dst);
  }

  void gd_OperateRegion(is_gt)(const GridData<Type,1> &src, const BBox &to,
                               const BBox &from, const Coords &step, BBox &bb)  {
    GridData<Type,1> &dst = *this;
    const int di = from.lower(0)-to.lower(0);

    BeginFastIndex1(src, src._bbox, src._data, const Type);
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      if (FastIndex1(dst,i) > FastIndex1(src,i+di))
        bb += Coords(1,i);
    end_for

    EndFastIndex1(dst);
    EndFastIndex1(src);
  }
  /************************************************************************/

  /****************************** >= ***************************************/
  void gd_OperateRegion(is_ge)(const Type &val, const BBox &to,
                               const Coords &step, BBox &bb) {
    GridData<Type,1> &dst = *this;
    
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      if (FastIndex1(dst,i) >= val) bb += Coords(1,i);
    end_for

    EndFastIndex1(dst);
  }

  void gd_OperateRegion(is_ge)(const GridData<Type,1> &src, const BBox &to,
                               const BBox &from, const Coords &step, BBox &bb)  {
    GridData<Type,1> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    
    BeginFastIndex1(src, src._bbox, src._data, const Type);
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      if (FastIndex1(dst,i) >= FastIndex1(src,i+di))
        bb += Coords(1,i);
    end_for

    EndFastIndex1(dst);
    EndFastIndex1(src);
  }
  /************************************************************************/

  /****************************** < ***************************************/
  void gd_OperateRegion(is_lt)(const Type &val, const BBox &to,
                               const Coords &step, BBox &bb) {
    GridData<Type,1> &dst = *this;

    BeginFastIndex1(dst, dst._bbox, dst._data, Type);
    
    for_1(i, to, step)
      if (FastIndex1(dst,i) < val) bb += Coords(1,i);
    end_for

    EndFastIndex1(dst);
  }

  void gd_OperateRegion(is_lt)(const GridData<Type,1> &src, const BBox &to,
                               const BBox &from, const Coords &step, BBox &bb)  {
    GridData<Type,1> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    
    BeginFastIndex1(src, src._bbox, src._data, const Type);
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);
   
    for_1(i, to, step)
      if (FastIndex1(dst,i) < FastIndex1(src,i+di))
        bb += Coords(1,i);
    end_for

    EndFastIndex1(dst);
    EndFastIndex1(src);
  }
  /************************************************************************/

  /****************************** <= ***************************************/
  void gd_OperateRegion(is_le)(const Type &val, const BBox &to,
                               const Coords &step, BBox &bb)  {
    GridData<Type,1> &dst = *this;

    BeginFastIndex1(dst, dst._bbox, dst._data, Type);
    
    for_1(i, to, step)
      if (FastIndex1(dst,i) <= val) bb += Coords(1,i);
    end_for

    EndFastIndex1(dst);
  }

  void gd_OperateRegion(is_le)(const GridData<Type,1> &src, const BBox &to,
                               const BBox &from, const Coords &step, BBox &bb)  {
    GridData<Type,1> &dst = *this;
    const int di = from.lower(0)-to.lower(0);

    BeginFastIndex1(src, src._bbox, src._data, const Type);
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      if (FastIndex1(dst,i) <= FastIndex1(src,i+di))
        bb += Coords(1,i);
    end_for

    EndFastIndex1(dst);
    EndFastIndex1(src);
  }

  /****************************** = ***************************************/
  void gd_OperateRegion(equal)(const Type &val, const BBox &to, const Coords &step) {
    GridData<Type,1> &dst = *this;

    BeginFastIndex1(dst, dst._bbox, dst._data, Type);
    
    for_1(i, to, step)
      FastIndex1(dst,i) = val;
    end_for
      
    EndFastIndex1(dst);
  }

  void gd_OperateRegion(equal)(const GridData<Type,1> &src, const BBox &to,
                               const BBox &from, const Coords &step) {
    GridData<Type,1> &dst = *this;
    const int di = from.lower(0)-to.lower(0);

    BeginFastIndex1(src, src._bbox, src._data, const Type);
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      FastIndex1(dst,i) = FastIndex1(src,i+di);
    end_for

    EndFastIndex1(dst);
    EndFastIndex1(src);
  }
  /************************************************************************/
    
  /****************************** + ***************************************/
  void gd_OperateRegion(plus)(const Type &val, const BBox &to, const Coords &step) {
    GridData<Type,1> &dst = *this;

    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      FastIndex1(dst,i) += val;
    end_for

    EndFastIndex1(dst);
  }

  void gd_OperateRegion(plus)(const GridData<Type,1> &src, const BBox &to,
                              const BBox &from, const Coords &step) {
    GridData<Type,1> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    
    BeginFastIndex1(src, src._bbox, src._data, const Type);
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      FastIndex1(dst,i) += FastIndex1(src,i+di);
    end_for

    EndFastIndex1(dst);
    EndFastIndex1(src);
  }
  /************************************************************************/

  /****************************** - ***************************************/
  void gd_OperateRegion(minus)(const Type &val, const BBox &to, const Coords &step) {
    GridData<Type,1> &dst = *this;
    
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);
    
    for_1(i, to, step)
      FastIndex1(dst,i) -= val;
    end_for

    EndFastIndex1(dst);
  }

  void gd_OperateRegion(minus)(const GridData<Type,1> &src, const BBox &to,
                               const BBox &from, const Coords &step) {
    GridData<Type,1> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    
    BeginFastIndex1(src, src._bbox, src._data, const Type);
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      FastIndex1(dst,i) -= FastIndex1(src,i+di);
    end_for

    EndFastIndex1(dst);
    EndFastIndex1(src);
  }
  /************************************************************************/

  /****************************** * ***************************************/
  void gd_OperateRegion(mult)(const Type &val, const BBox &to, const Coords &step) {
    GridData<Type,1> &dst = *this;

   BeginFastIndex1(dst, dst._bbox, dst._data, Type);

   for_1(i, to, step)
     FastIndex1(dst,i) *= val;
   end_for
     
   EndFastIndex1(dst);
  }

  void gd_OperateRegion(mult)(const GridData<Type,1> &src, const BBox &to,
                              const BBox &from, const Coords &step) {
    GridData<Type,1> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    
    BeginFastIndex1(src, src._bbox, src._data, const Type);
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      FastIndex1(dst,i) *= FastIndex1(src,i+di);
    end_for

    EndFastIndex1(dst);
    EndFastIndex1(src);
  }
  /************************************************************************/

  /****************************** / ***************************************/
  void gd_OperateRegion(div)(const Type &val, const BBox &to, const Coords &step) {
    assert (val != (Type)0);
    GridData<Type,1> &dst = *this;
    
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);
   
    for_1(i, to, step)
      FastIndex1(dst,i) /= val;
    end_for

    EndFastIndex1(dst);
  }

  void gd_OperateRegion(div)(const GridData<Type,1> &src, const BBox &to,
                             const BBox &from, const Coords &step) {
    GridData<Type,1> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    
    BeginFastIndex1(src, src._bbox, src._data, const Type);
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      FastIndex1(dst,i) /= FastIndex1(src,i+di);
    end_for

    EndFastIndex1(dst);
    EndFastIndex1(src);
  }
  /************************************************************************/
  
  /****************************** min ***************************************/
  void gd_OperateRegion(min)(const Type &val, const BBox &to, Coords const &step) {
    GridData<Type,1> &dst = *this;

    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      FastIndex1(dst,i) = Min(val, FastIndex1(dst,i));
    end_for

    EndFastIndex1(dst);
  }

  void gd_OperateRegion(min)(const GridData<Type,1> &src, const BBox &to,
                             const BBox &from, Coords const &step) {
    GridData<Type,1> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    
    BeginFastIndex1(src, src._bbox, src._data, const Type);
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      FastIndex1(dst,i) = Min(FastIndex1(dst,i), FastIndex1(src,i+di));
    end_for

    EndFastIndex1(dst);
    EndFastIndex1(src);
  }
  /************************************************************************/

  /****************************** max ***************************************/
  void gd_OperateRegion(max)(const Type &val, const BBox &to, Coords const &step) {
    GridData<Type,1> &dst = *this;

    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      FastIndex1(dst,i) = Max(val, FastIndex1(dst,i));
    end_for

    EndFastIndex1(dst);
  }

  void gd_OperateRegion(max)(const GridData<Type,1> &src, const BBox &to,
                             const BBox &from, Coords const &step) {
    GridData<Type,1> &dst = *this;
    const int di = from.lower(0)-to.lower(0);
    
    BeginFastIndex1(src, src._bbox, src._data, const Type);
    BeginFastIndex1(dst, dst._bbox, dst._data, Type);

    for_1(i, to, step)
      FastIndex1(dst,i) = Max(FastIndex1(dst,i), FastIndex1(src,i+di));
    end_for

    EndFastIndex1(dst);
    EndFastIndex1(src);
  }
  /************************************************************************/

public:
  /* Perform efficient buffer packing/unpacking and data movement */
  void PackRegion(Type *sendbuf, const BBox &from) const {
    const GridData<Type,1> &src = *this;

    if (from==src._bbox) 
      std::memcpy((void *) sendbuf, (void *) src._data, sizeof(Type)*from.size());
    else {
      register int p = 0;
      
      BeginFastIndex1(src, src._bbox, src._data, const Type);
      
      for_1(i, from, _step)
        sendbuf[p++] = FastIndex1(src,i);
      end_for
        
      EndFastIndex1(src);
    }
  } 

  void UnPackRegion(const Type *recvbuf, const BBox &to) {
    GridData<Type,1> &dst = *this;

    if (to==dst._bbox) 
      std::memcpy((void *) dst._data, (void *) recvbuf, sizeof(Type)*to.size());
    else {
      register int p = 0;

      BeginFastIndex1(dst, dst._bbox, dst._data, Type);

      for_1(i, to, _step)
        FastIndex1(dst,i) = recvbuf[p++];
      end_for

      EndFastIndex1(dst);
    }
  }

  friend std::ostream& operator<< <>(std::ostream &os, 
                                     const GridData<Type,1> &gd);
  friend std::ofstream& operator<< <>(std::ofstream &ofs, 
                                      const GridData<Type,1> &gd);
  friend std::ifstream& operator>> <>(std::ifstream &ifs, 
                                      GridData<Type,1> &gd);  
  friend std::stringstream& operator<< <>(std::stringstream &ofs, 
                                          const GridData<Type,1> &gd);
  friend std::stringstream& operator>> <>(std::stringstream &ifs, 
                                          GridData<Type,1> &gd);
};

template <class Type>
std::ostream&  operator << (std::ostream& os, const GridData<Type,1> &gd) {
  if (&gd == (const GridData<Type,1> *)0) return os;

  os << "BBox: " << gd.bbox() << " " << "Extents: " << gd.extents() << " ";
  os << "Step: " << gd.stepsize() << " ";
  os << "Size: " << gd.size() << " " << "Bottom: " << gd.bottom() << " ";
  
  os << "\n";
  
  if (gd.bbox().empty()) return os;

  const Coords &l = gd.lower();
  const Coords &u = gd.upper();
  const Coords &step = gd.stepsize();
  for (register int i=l(0);i<=u(0);i+=step(0)) 
    { os << "[" << i << "]=" << gd(i) << " "; }
  os << "\n";

  return os;
}

template <class Type>
std::ofstream&  operator << (std::ofstream& ofs, const GridData<Type,1> &gd) {
  if (&gd == (GridData<Type,1> *)0) return ofs;
  
  ofs.write((char*)&gd.bbox(),sizeof(BBox));
  ofs.write((char*)gd.data(),gd.size()*sizeof(Type));

  return ofs;
}

template <class Type>
std::ifstream&  operator >> (std::ifstream& ifs, GridData<Type,1> &gd) {
  if (&gd == (GridData<Type,1> *)0) return ifs;
  
  BBox bb;
  ifs.read((char*)&bb,sizeof(BBox));
  
  if (!gd.ok_to_index())
    { gd.allocate(bb); }
  
  if (bb == gd.bbox()) {
    ifs.read((char*)gd.data(),gd.size()*sizeof(Type));
  }
  else {
    GridData<Type,1> tmpgd(bb);
    ifs.read((char*)tmpgd.data(),tmpgd.size()*sizeof(Type));
    gd.copy(tmpgd);  
  } 
  
  return ifs;
}

template <class Type>
std::stringstream& operator<<(std::stringstream& ofs, const GridData<Type,1> &gd) {
  if (&gd == (GridData<Type,1> *)0) return ofs;
  
  ofs.write((char*)&gd.bbox(),sizeof(BBox));
  ofs.write((char*)gd.data(),gd.size()*sizeof(Type));
  
  return ofs;
}

template <class Type>
std::stringstream& operator>>(std::stringstream& ifs, GridData<Type,1> &gd) {
  if (&gd == (GridData<Type,1> *)0) return ifs;
  
  BBox bb;
  ifs.read((char*)&bb,sizeof(BBox));
  
  if (!gd.ok_to_index())
    { gd.allocate(bb); }
  
  if (bb == gd.bbox()) {
    ifs.read((char*)gd.data(),gd.size()*sizeof(Type));
  }
  else {
    GridData<Type,1> tmpgd(bb);
    ifs.read((char*)tmpgd.data(),tmpgd.size()*sizeof(Type));
    gd.copy(tmpgd);  
  } 

  return ifs;
}

/*********************************************************************/
template <class Type, int dim>
inline Coords lower(GridData<Type,1> &gd) { return(gd.lower()); }
template <class Type, int dim>
inline Coords upper(GridData<Type,1> &gd) { return(gd.upper()); }
template <class Type, int dim>
inline Coords extents(GridData<Type,1> &gd) { return(gd.extents()); }
template <class Type, int dim>
inline int lower(GridData<Type,1> &gd, const int i) { return(gd.lower(i)); }
template <class Type, int dim>
inline int upper(GridData<Type,1> &gd, const int i) { return(gd.upper(i)); }
template <class Type, int dim>
inline int extents(GridData<Type,1> &gd, const int i) { return(gd.extents(i)); }
template <class Type, int dim>
inline BBox bbox(GridData<Type,1> &gd) { return(gd.bbox()); }
template <class Type, int dim>
inline int size(GridData<Type,1> &gd) { return(gd.size()); }
template <class Type, int dim>
inline Coords stepsize(GridData<Type,1> &gd) { return(gd.stepsize()); }

#endif

---