regularData1D.h

Go to the documentation of this file.

// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//

#ifndef BALL_DATATYPE_REGULARDATA1D_H
#define BALL_DATATYPE_REGULARDATA1D_H

#ifndef BALL_COMMON_H
# include <BALL/common.h>
#endif

#ifndef BALL_SYSTEM_FILE_H
# include <BALL/SYSTEM/file.h>
#endif

#ifndef BALL_SYSTEM_BINARYFILEADAPTOR_H
# include <BALL/SYSTEM/binaryFileAdaptor.h>
#endif

#include <vector>
#include <iostream>
#include <fstream>
#include <iterator>
#include <algorithm>

namespace BALL
{
  template <typename ValueType>
  class TRegularData1D
  {
    public:
      
    BALL_CREATE(TRegularData1D<ValueType>)

    

    typedef Position IndexType;
    typedef std::vector<ValueType>  VectorType;
    typedef double CoordinateType;
    typedef typename std::vector<ValueType>::iterator Iterator;
    typedef typename std::vector<ValueType>::const_iterator ConstIterator;

    //  STL compatibility types
    //
    typedef ValueType value_type;
    typedef typename std::vector<ValueType>::iterator iterator;
    typedef typename std::vector<ValueType>::const_iterator const_iterator;
    typedef typename std::vector<ValueType>::reference reference;
    typedef typename std::vector<ValueType>::const_reference const_reference;
    typedef typename std::vector<ValueType>::pointer pointer;
    typedef typename std::vector<ValueType>::difference_type difference_type;
    typedef typename std::vector<ValueType>::size_type size_type;

      
    TRegularData1D();

    TRegularData1D(const TRegularData1D& data);
  
    TRegularData1D(const CoordinateType& origin, const CoordinateType& dimension, const CoordinateType& spacing);

    TRegularData1D(const IndexType& size);

    TRegularData1D(const VectorType& data, const CoordinateType& origin = 0.0, const CoordinateType& dimension = 1.0);

    virtual ~TRegularData1D();

    virtual void clear();



    TRegularData1D& operator = (const TRegularData1D<ValueType>& data);

    TRegularData1D& operator = (const VectorType& data);


    bool operator == (const TRegularData1D& data) const;

    BALL_INLINE bool operator != (const TRegularData1D& data) const { return !this->operator == (data); }

    BALL_INLINE bool empty() const { return data_.empty(); }

    bool isInside(const CoordinateType& x) const;

    BALL_INLINE ConstIterator begin() const { return data_.begin(); }
    BALL_INLINE ConstIterator end() const { return data_.end(); }
    BALL_INLINE Iterator begin() { return data_.begin(); }
    BALL_INLINE Iterator end() { return data_.end(); }

    
    // STL compatibility
    BALL_INLINE size_type size() const { return data_.size(); }
    BALL_INLINE size_type max_size() const { return data_.max_size(); }
    BALL_INLINE void swap(TRegularData1D<ValueType>& data) { std::swap(*this, data); }

    const ValueType& getData(const IndexType& index) const;

    ValueType& getData(const IndexType& index);

    const ValueType& operator [] (const IndexType& index) const { return data_[index]; }
      
    ValueType& operator [] (const IndexType& index) { return data_[index]; }
    
    ValueType operator () (const CoordinateType& x) const;

    ValueType getInterpolatedValue(const CoordinateType& x) const;
      
    void getEnclosingIndices(const CoordinateType& x, Position& lower, Position& upper) const;

    void getEnclosingValues(const CoordinateType& x, ValueType& lower, ValueType& upper) const;

    CoordinateType getCoordinates(const IndexType& index) const;

    IndexType getClosestIndex(const CoordinateType& x) const;
      
    IndexType getLowerIndex(const CoordinateType& x) const;
      
    const ValueType& getClosestValue(const CoordinateType& x) const;
      
    ValueType& getClosestValue(const CoordinateType& x);
      
    BALL_INLINE IndexType getSize() const { return (IndexType)data_.size(); }

    BALL_INLINE const CoordinateType& getOrigin() const { return origin_; }
    
    BALL_INLINE const CoordinateType& getSpacing() const {  return spacing_; }
          
    BALL_INLINE void setOrigin(const CoordinateType& origin) { origin_ = origin; }

    BALL_INLINE const CoordinateType& getDimension() const { return dimension_; }

    BALL_INLINE void setDimension(const CoordinateType& dimension) { dimension_ = dimension; }

    void resize(const IndexType& size);

    void rescale(const IndexType& new_size);

    ValueType calculateMean() const;
    
    ValueType calculateSD() const;
    
    void binaryWrite(const String& filename) const;

    void binaryRead(const String& filename);
  
    
    protected:
    CoordinateType  origin_;

    CoordinateType  dimension_;

    CoordinateType  spacing_;

    VectorType      data_;

    typedef struct { ValueType bt[1024]; } BlockValueType;
  };

  typedef TRegularData1D<float> RegularData1D;
  
  template <typename ValueType>
  TRegularData1D<ValueType>::TRegularData1D()
    : origin_(0.0),
      dimension_(0.0),
      spacing_(1.0),
      data_()
  {
  }

  template <typename ValueType>
  TRegularData1D<ValueType>::~TRegularData1D()
  {
  }

  template <typename ValueType>
  TRegularData1D<ValueType>::TRegularData1D(const TRegularData1D<ValueType>& data)
    : origin_(data.origin_),
      dimension_(data.dimension_),
      spacing_(data.spacing_),
      data_()
  {
    // Try to catch allocation errors and rethrow them as OutOfMemory
    try
    {
      data_ = data.data_;
    }
    catch (std::bad_alloc&)
    {
      throw Exception::OutOfMemory(__FILE__, __LINE__, data.size() * sizeof(ValueType));
    }
  }

  template <typename ValueType>
  TRegularData1D<ValueType>::TRegularData1D
    (const typename TRegularData1D<ValueType>::CoordinateType& origin, 
     const typename TRegularData1D<ValueType>::CoordinateType& dimension, 
     const typename TRegularData1D<ValueType>::CoordinateType& spacing)
    : origin_(origin),
      dimension_(dimension),
      spacing_(spacing),
      data_()
  {
    // Determine the size of the vector
    size_type size = (size_type)(dimension_ / spacing_ + 1.0);

    // Try to catch allocation errors and rethrow them as OutOfMemory
    try
    {
      data_.resize(size);
    }
    catch (std::bad_alloc&)
    {
      throw Exception::OutOfMemory(__FILE__, __LINE__, size * sizeof(ValueType));
    }
  }
      
  template <typename ValueType>
  TRegularData1D<ValueType>::TRegularData1D
    (const typename TRegularData1D<ValueType>::VectorType& data, 
     const typename TRegularData1D<ValueType>::CoordinateType& origin, 
     const typename TRegularData1D<ValueType>::CoordinateType& dimension)
    : origin_(origin),
      dimension_(dimension),
      spacing_(dimension / ((double)data.size()-1)),
      data_()
  {
    // Try to catch allocation errors and rethrow them as OutOfMemory
    try
    {
      data_ = data;
    }
    catch (std::bad_alloc&)
    {
      throw Exception::OutOfMemory(__FILE__, __LINE__, data.size() * sizeof(ValueType));
    }
  }


  template <class ValueType>
  TRegularData1D<ValueType>::TRegularData1D
    (const typename TRegularData1D<ValueType>::IndexType& size)
    : origin_(0.0),
      dimension_(1.0),
      data_()
  {
    // Compute the grid spacing
    spacing_ = dimension_ / (double)(size - 1);

    try
    {
      data_.resize(size);
    }
    catch (std::bad_alloc&)
    {
      data_.resize(0);
      throw Exception::OutOfMemory(__FILE__, __LINE__, size * sizeof(ValueType));
    }
  }

  template <typename ValueType>
  void TRegularData1D<ValueType>::clear()
  {
    // iterate over the data and reset all values to their default
    // boundaries and vector size remain unchanged
    static ValueType default_value = ValueType();
    std::fill(data_.begin(), data_.end(), default_value);
  }

  template <typename ValueType>
  TRegularData1D<ValueType>& TRegularData1D<ValueType>::operator = (const TRegularData1D<ValueType>& rhs)
  {
    // copy all members...
    origin_ = rhs.origin_;
    dimension_ = rhs.dimension_;
    spacing_ = rhs.spacing_;
    try 
    {
      data_ = rhs.data_;
    }
    catch (std::bad_alloc&)
    {
      data_.resize(0);
      throw Exception::OutOfMemory(__FILE__, __LINE__, rhs.size() * sizeof(ValueType));
    }

    return *this;
  }

  template <typename ValueType>
  TRegularData1D<ValueType>& TRegularData1D<ValueType>::operator = (const VectorType& rhs)
  {
    // Copy the data. The boundaries remain unchanged.
    try 
    {
      data_ = rhs;
    }
    catch (std::bad_alloc&)
    {
      data_.resize(0);
      throw Exception::OutOfMemory(__FILE__, __LINE__, rhs.size() * sizeof(ValueType));
    }
    
    return *this;
  }

  template <typename ValueType>
  bool TRegularData1D<ValueType>::operator == (const TRegularData1D<ValueType>& data) const
  {
    return (origin_ == data.origin_ 
            && dimension_ == data.dimension_ 
            && data_ == data.data_);
  }

  template <class ValueType>
  BALL_INLINE
  bool TRegularData1D<ValueType>::isInside(const typename TRegularData1D<ValueType>::CoordinateType& r) const
  {
    return ((r >= origin_) && (r <= (origin_ + dimension_)));
  }

  template <typename ValueType>
  BALL_INLINE
  const ValueType& TRegularData1D<ValueType>::getData(const IndexType& index) const
  {
    if (index >= data_.size())
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }
    return data_[index];
  }
    
  template <typename ValueType>
  BALL_INLINE
  ValueType& TRegularData1D<ValueType>::getData(const IndexType& index)
  {
    if (index >= data_.size())
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }
    return data_[index];
  }
    
  template <typename ValueType>
  void TRegularData1D<ValueType>::getEnclosingIndices
    (const typename TRegularData1D<ValueType>::CoordinateType& x,
     Position& lower, Position& upper) const
  {
    if (!isInside(x) || (data_.size() < 2))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }
    lower = (Position)floor((x - origin_) / spacing_);
    if (lower == data_.size() - 1)
    {
      // If we are on the right most data point, we cannot interpolate to the right!
      lower = data_.size() - 2;
    }
    upper = lower + 1;
  }

  template <typename ValueType>
  void TRegularData1D<ValueType>::getEnclosingValues
    (const typename TRegularData1D<ValueType>::CoordinateType& x,
     ValueType& lower, ValueType& upper) const
  {
    Position lower_index;
    Position upper_index;
    getEnclosingIndices(x, lower_index, upper_index);
    lower = data_[lower_index];
    upper = data_[upper_index];
  }

  template <typename ValueType>
  BALL_INLINE
  ValueType TRegularData1D<ValueType>::getInterpolatedValue(const CoordinateType& x) const
  {
    if (!isInside(x))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }
    return operator () (x);
  }
      
  template <typename ValueType>
  BALL_INLINE
  typename TRegularData1D<ValueType>::CoordinateType TRegularData1D<ValueType>::getCoordinates
    (const typename TRegularData1D<ValueType>::IndexType& index) const
  {
    if ((index >= data_.size()) || (data_.size() == 0))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }

    return (CoordinateType)(origin_ + (double)index / ((double)data_.size()-1) * dimension_);
  }

  template <typename ValueType>
  BALL_INLINE
  typename TRegularData1D<ValueType>::IndexType TRegularData1D<ValueType>::getClosestIndex(const CoordinateType& x) const
  {
    if ((x < origin_) || (x > (origin_ + dimension_)))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }

    return (IndexType)(size_type)floor((x - origin_) / spacing_ + 0.5);
  }
      
  template <typename ValueType>
  BALL_INLINE
  typename TRegularData1D<ValueType>::IndexType TRegularData1D<ValueType>::getLowerIndex(const CoordinateType& x) const
  {
    if ((x < origin_) || (x > (origin_ + dimension_)))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }

    return (IndexType)(size_type)floor((x - origin_) / spacing_);
  }
      
  template <typename ValueType>
  BALL_INLINE
  const ValueType& TRegularData1D<ValueType>::getClosestValue(const CoordinateType& x) const
  {
    if ((x < origin_) || (x > (origin_ + dimension_)))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }
    
    // Round to the closest data point.
    size_type index = (size_type)floor((x - origin_) / spacing_ + 0.5);
    return data_[index];
  }
      
  template <typename ValueType>
  BALL_INLINE
  ValueType& TRegularData1D<ValueType>::getClosestValue(const CoordinateType& x)
  {
    if ((x < origin_) || (x > (origin_ + dimension_)))
    {
      throw Exception::OutOfGrid(__FILE__, __LINE__);
    }
    
    // Round to the closest data point.
    size_type index = (size_type)floor((x - origin_) / spacing_ + 0.5);
    return data_[index];
  }
      
  template <typename ValueType>
  BALL_INLINE
  ValueType TRegularData1D<ValueType>::calculateMean() const
  {
    IndexType data_points = this->getSize();
    ValueType mean = 0;
    for (IndexType i = 0; i < data_points; i++)
    {
      mean += data_[i];
    }
    mean /= data_points;
    return mean;
  }
  
  template <typename ValueType>
  BALL_INLINE
  ValueType TRegularData1D<ValueType>::calculateSD() const
  {
    IndexType data_points = this->getSize();
    ValueType stddev = 0;
    ValueType mean = this->calculateMean();
    for (IndexType i = 0; i < data_points; i++)
    {
      stddev += (pow(data_[i]-mean,2));
    }
    stddev /= (data_points-1);
    stddev = sqrt(stddev);
    return stddev;
  }
    
  template <typename ValueType>
  BALL_INLINE
  ValueType TRegularData1D<ValueType>::operator () (const CoordinateType& x) const
  {
    size_type left_index = (size_type)floor((x - origin_) / spacing_);
    if (left_index == data_.size() - 1)
    {
      // If we are on the right most data point, we cannot interpolate to the right!
      return data_[data_.size() - 1];
    }
    
    // Interpolate between the point to the left and the point to the right.
    double d = 1.0 - (((x - origin_) - (double)left_index * spacing_) / spacing_);
    return data_[left_index] * d + (1.0 - d) * data_[left_index + 1];
  }
      
  template <typename ValueType>
  void TRegularData1D<ValueType>::resize
    (const typename TRegularData1D<ValueType>::IndexType& new_size)
  {
    // Rescale dimension to the new size.
    if (data_.size() > 0)
    {
      dimension_ *= (double)new_size / (double)data_.size();
    }

    // Try to resize the vactor and rethrow any bad_allocs.
    try
    {
      data_.resize(new_size);
    }
    catch (std::bad_alloc&)
    {
      // The resulting vector is empty and thus well-defined.
      data_.resize(0);
      throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * sizeof(ValueType));
    }
  }
    
  template <typename ValueType>
  void TRegularData1D<ValueType>::rescale
    (const typename TRegularData1D<ValueType>::IndexType& new_size)
  {
    // if the new and the old size coincide: done.
    if (new_size == (IndexType)data_.size())
    {
      return;
    }

    // Catch any bad_allocs throw by vector::resize
    try
    {
      // if the data set is empty...
      if (data_.size() == 0)
      {
        // ...there's nothing to do: a resize was requested
        data_.resize(new_size);
        return;
      }
      
      // if the data set contains only a single value,
      // we fill everything with this value
      if ((data_.size() == 1) && (new_size > 1))
      {
        ValueType old_value = data_[0];
        data_.resize(new_size);
        for (IndexType i = 1; i < new_size; i++)
        {
          data_[i] = old_value;
        }

        return;
      }

      // that's the default case: use linear interpolation
      // to determine the values at the new positions
      VectorType new_data(new_size);
      CoordinateType factor1 = (CoordinateType)data_.size() / (CoordinateType)new_size;
      CoordinateType factor2 = (CoordinateType)(data_.size() - 1) / (new_size - 1);

      for (Size i = 0; i < new_size; i++)
      {
        // determine the interval of the old data set we are currently in
        // ([old_idx, old_idx + 1])
        IndexType old_idx = (IndexType)((CoordinateType)i * factor1);

        // consider numerical inaccuracies...
        if (old_idx >= (data_.size() - 1))
        {
          old_idx = data_.size() - 2;
        }
        CoordinateType factor3 = (CoordinateType)i * factor2 - (CoordinateType)old_idx;
        new_data[i] = data_[old_idx] * (1 - factor3) + factor3 * data_[old_idx + 1];
      }

      // assign the new data
      data_ = new_data;
    }
    catch (std::bad_alloc&)
    {
      // Make sure we are in a well-defined state.
      data_.resize(0);
      throw Exception::OutOfMemory(__FILE__, __LINE__, new_size * sizeof(ValueType));     
    }
  }

  template <typename ValueType>
  std::ostream& operator << (std::ostream& os, const TRegularData1D<ValueType>& data)
  {
    // Write the grid origin, dimension, and number of grid points
    os << data.getOrigin() << std::endl
       << data.getOrigin() + data.getDimension() << std::endl
       << data.getSize() - 1 << std::endl;

    // Write the array contents.
    std::copy(data.begin(), data.end(), std::ostream_iterator<ValueType>(os, "\n"));
    return os;
  }

  template <typename ValueType>
  std::istream& operator >> (std::istream& is, TRegularData1D<ValueType>& grid)
  {
    typename TRegularData1D<ValueType>::CoordinateType origin;
    typename TRegularData1D<ValueType>::CoordinateType dimension;
    typename TRegularData1D<ValueType>::IndexType size;

    is >> origin;
    is >> dimension;
    is >> size;
    
    dimension -= origin;
    size++;

    grid.resize(size);
    grid.setOrigin(origin);
    grid.setDimension(dimension);

    std::copy(std::istream_iterator<ValueType>(is), 
              std::istream_iterator<ValueType>(), 
              grid.begin());
    //    std::copy_n(std::istream_iterator<ValueType>(is), grid.size(), grid.begin());
    
    return is;
  }

  template <typename ValueType>
  void TRegularData1D<ValueType>::binaryWrite(const String& filename) const
  {
    File outfile(filename, std::ios::out|std::ios::binary);
    if (!outfile.isValid()) 
    {
      throw Exception::FileNotFound(__FILE__, __LINE__, filename);
    }

    BinaryFileAdaptor<BlockValueType> adapt_block;
    BinaryFileAdaptor<ValueType>      adapt_single;
    
    // write all information we need to recreate the grid
    BinaryFileAdaptor<CoordinateType> adapt_coordinate;
    BinaryFileAdaptor<Size>           adapt_size;

    adapt_size.setData(data_.size());
    outfile << adapt_size;
    
    adapt_coordinate.setData(origin_);
    outfile << adapt_coordinate;

    adapt_coordinate.setData(dimension_);
    outfile << adapt_coordinate;

    adapt_coordinate.setData(spacing_);
    outfile << adapt_coordinate;

    // we slide a window of size 1024 over our data
    Index window_pos = 0;
    while (((int)data_.size() - (1024 + window_pos)) >= 0 )
    {
      adapt_block.setData(*(BlockValueType*)&(data_[window_pos]));
      outfile << adapt_block;
      window_pos += 1024;
    }

    // now we have to write the remaining data one by one
    for (Size i = window_pos; i < data_.size(); i++)
    {
      adapt_single.setData(data_[i]);
      outfile << adapt_single;
    }

    // that's it. I hope...
    outfile.close();
  }

  template <typename ValueType>
  void TRegularData1D<ValueType>::binaryRead(const String& filename)
  {
    File infile(filename, std::ios::in|std::ios::binary);
    if (!infile.isValid()) throw Exception::FileNotFound(__FILE__, __LINE__, filename);
    
    BinaryFileAdaptor<BlockValueType> adapt_block;
    BinaryFileAdaptor<ValueType>      adapt_single;
    
    // read all information we need to recreate the grid
    BinaryFileAdaptor<CoordinateType> adapt_coordinate;
    BinaryFileAdaptor<Size>           adapt_size;

    infile >> adapt_size;
    Size new_size = adapt_size.getData();
  
    infile >> adapt_coordinate;
    origin_ = adapt_coordinate.getData();

    infile >> adapt_coordinate;
    dimension_ = adapt_coordinate.getData();

    infile >> adapt_coordinate;
    spacing_ = adapt_coordinate.getData();

    data_.resize(new_size);

    // we slide a window of size 1024 over our data
    Index window_pos = 0;
    
    while ( ((int)data_.size() - (1024 + window_pos)) >= 0 )
    {
      infile >> adapt_block;
      *(BlockValueType*)(&(data_[window_pos])) = adapt_block.getData();
      /*
      for (Size i=0; i<1024; i++)
      {
        data_[i+window_pos] = adapt_block.getData().bt[i];
      }
      */
      window_pos+=1024;
    }

    // now we have to read the remaining data one by one
    for (Size i=window_pos; i<data_.size(); i++)
    {
      infile >> adapt_single;
      data_[i] = adapt_single.getData();
    }

    // that's it. I hope...
    infile.close();
  }
} // namespace BALL

#endif // BALL_DATATYPE_REGULARDATA1D_H