Logo Search packages:      
Sourcecode: magics++ version File versions  Download package

MatrixHandler.h

Go to the documentation of this file.
/******************************** LICENSE ********************************

 Copyright 2007 European Centre for Medium-Range Weather Forecasts (ECMWF)

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at 

    http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.

 ******************************** LICENSE ********************************/

/*! \file MatrixHandler.h
    \brief Definition of the Template class MatrixHandler.
    
    Magics Team - ECMWF 2004
    
    Started: Wed 18-Feb-2004 
    
    Changes:
    
*/ 

#ifndef MatrixHandler_H
#define MatrixHandler_H
 
#include "magics.h"
#include "Matrix.h"
#include "BasePointsHandler.h"
#include "Transformation.h"
#include "VectorOfPointers.h"
#include "Polyline.h"


namespace magics {
      
struct CoordIndex
{
      CoordIndex(double coordinate, int index) : coordinate_(coordinate), index_(index) {}
      ~CoordIndex() {}
      bool  operator<(const CoordIndex& other) const
    {
        return coordinate_ < other.coordinate_;
    }
    
    bool  operator<(double& val) const
    {
        return coordinate_ < val;
    }
    
    bool  operator == (double val) const
    {
        return coordinate_ == val;
    }
 
    bool  operator == (CoordIndex& other) const
    {
        return coordinate_ == other.coordinate_;
    }
      
      double coordinate_;
      int   index_;
       
};

struct FinderHelper 
{
    FinderHelper(double ref) : reference_(ref) {}
    ~FinderHelper() {}
    bool operator()(const CoordIndex& left, const CoordIndex& right) const
    {
        return (left.coordinate_ <= reference_ && reference_ <= right.coordinate_);
    }
    
    double reference_;
    
    
};

class CoordinatesFinder : public vector<CoordIndex>
{
public :
    void find(double val, double& lowerval, int& lower, double& upperval, int& upper) const {
        const_iterator coord = std::adjacent_find(begin(), end(), FinderHelper(val));
        
        if ( coord == end()) {
            upper = -1;
            lower = -1;        
            return;     
        }
        upper = (coord+1)->index_;
        upperval = (coord+1)->coordinate_;
        lower = coord->index_;
        lowerval = coord->coordinate_;
        return;
        
    }
    
    int find(double val) const {          
        const_iterator coord = std::find(begin(), end(), val);        
        return  ( coord == end() ) ? -1 : coord->index_;;        
    }
    
    void add(float val, int index) {
        push_back(CoordIndex(val, index));
    }
    
    void print() {
       
        for (const_iterator coord = begin(); coord != end(); ++coord) {
              Log::dev()<< "[" << coord->coordinate_ << ", " << coord->index_ << "]" << endl;
        }
    }


};


template <class P>
class MatrixHandler : public AbstractMatrix, public AbstractPoints<P>
{
public : 
    MatrixHandler(const AbstractMatrix& matrix) : AbstractMatrix(),
        AbstractPoints<P>(),
        matrix_(matrix) {}
    MatrixHandler(const MatrixHandler<P>& matrix) :        
        AbstractMatrix(),
        AbstractPoints<P>(),
        matrix_(matrix) {}
    virtual ~MatrixHandler() {}
   
    virtual double operator()(int  i, int  j) const { return matrix_(i, j); } 
    
    virtual int rowIndex(double r) const     { return matrix_.rowIndex(r); } 
    virtual int columnIndex(double c) const  { return matrix_.columnIndex(c); }  
    virtual bool  akimaEnable() const  { return matrix_.akimaEnable(); }
    
    
    virtual void boundRow(double r, 
      double& row1, int& index1, double& row2, int& index2) const 
            { return matrix_.boundRow(r, row1, index1, row2, index2); }
    virtual void boundColumn(double r, 
      double& column1, int& index1, double& column2, int& index2) const 
            { return matrix_.boundColumn(r, column1, index1, column2, index2); }          
    
    double  left() const { return matrix_.left(); }
         double bottom() const { return matrix_.bottom(); } 
         double  right() const { return matrix_.right(); }
         double  top() const { return matrix_.top(); }
         
         double x(double x, double y) const  { return matrix_.x(x, y); }
         double y(double x, double y) const { return matrix_.y(x, y); }
         
    
   
   
 
      
    virtual double interpolate(double  i, double  j) const 
    {
      int ii = rowIndex(i);
      if (ii == -1) {
            // interpolate between 2 rows.
            double v1, v2;
            int i1, i2;
            boundRow(i, v1, i1, v2, i2);
            
            if (i1 == -1) return missing(); 
            
          double a = (*this).interpolate(v1, j);
            double b = (*this).interpolate(v2, j);
          
            if ( same(a, missing()) || same(b, missing()) ) return missing();
            
            double da = (v2-i)/(v2-v1);
            double db = (i-v1)/(v2-v1);
            double val = (a*da) + (b*db);
            return val;
        }
        int jj = columnIndex(j);
        if (jj == -1) {
            double v1, v2;
            int i1, i2;
            boundColumn(j, v1, i1, v2, i2);
            if (i1 == -1) return missing(); 
            
            
            double a = (*this)(ii, i1);
            double b = (*this)(ii, i2);
            
            if ( same(a, missing()) || same(b, missing()) ) return missing();
            
            double da = (v2-j)/(v2-v1);
            double db = (j-v1)/(v2-v1);
            double val = (a*da) + (b*db);
            return val;
            
        }
      
      return (*this)(ii, jj);
                  
      
    }  
    
  
    
    virtual int    rows() const { return matrix_.rows(); }
    virtual int    columns() const { return matrix_.columns(); }
    virtual int    lowerRow(double v) const { return matrix_.lowerRow(v); }
    virtual int    lowerColumn(double v) const { return matrix_.lowerColumn(v); } 
    virtual double  XResolution() const { return matrix_.XResolution(); } 
    virtual double  YResolution() const { return matrix_.YResolution(); } 
    virtual double  width() const { return matrix_.width(); } 
    virtual double  height() const { return matrix_.height(); } 
    
    virtual const AbstractMatrix& original() const { return matrix_.original(); }
    virtual int firstRow() const    { return matrix_.firstRow(); }
    virtual int nextRow(int i, int f) const   { return matrix_.nextRow(i, f); }
    virtual int firstColumn() const { return matrix_.firstColumn(); }
    virtual int nextColumn(int j, int f) const  { return matrix_.nextColumn(j, f); } 
     
 
    virtual double  minX() const { return matrix_.minX(); }
    virtual double  maxX() const { return matrix_.maxX(); }
    virtual double  minY() const { return matrix_.minY(); }
    virtual double  maxY() const { return matrix_.maxY(); }

    // Implements the AbstractPoints interface
    virtual void setToFirst() const  {
        if (points_.empty()) {
            int nb_rows = rows();
            int nb_columns = columns();
            
            points_.reserve(nb_rows * nb_columns);
            
            for (int r = 0; r < nb_rows; r++) {
                for (int c = 0; c < nb_columns; c++) {
                    if ( matrix_.accept(column(r, c), row(r, c)) ) 
                        if ( !same((*this)(r, c), matrix_.missing() ) ) 
                              points_.push_back(new P(column(r,c), row(r,c), (*this)(r, c)));
                }
            }            
        }
        current_ = points_.begin();
    }
    
    //! Method to test the end of collection.
    virtual bool more() const  {
        return current_ != points_.end();
    }
   
    virtual bool accept(double x, double y) const { return matrix_.accept(x, y); }
  
    virtual const P& current()  const 
    {
        return **current_;
       
    }
    
    virtual void advance() const {
        current_++;
    }
    
    
    virtual vector<double>&  rowsAxis()  const { return const_cast<MatrixHandler*>(this)->matrix_.rowsAxis(); }
    virtual vector<double>&  columnsAxis() const  { return const_cast<MatrixHandler*>(this)->matrix_.columnsAxis(); }
    
    virtual double  row(int i, int j) const { 
        return matrix_.row(i, j); 
    }
    virtual double  column(int i, int j) const { 
        return matrix_.column(i, j); }
    
    virtual double  real_row(double i, double j) const { 
        return matrix_.real_row(i, j); 
    }
    virtual double  real_column(double i, double j) const { 
        return matrix_.real_column(i, j); }
    virtual double  regular_row(int i) const { 
        return matrix_.regular_row(i); 
    }
    virtual double  regular_column(int i) const { 
        return matrix_.regular_column(i); 
    }
    
    
    virtual double  missing() const  { return matrix_.missing(); }
    virtual bool  hasMissingValues() const  { 
       
        for (int r = 0; r < rows(); r++) {
                for (int c = 0; c < columns(); c++) {                   
                        if ( this->operator()(r, c) == matrix_.missing() ) 
                              return true;
                }
            }
        return false;
    }    
    
    
    
 
    
protected:    
    const AbstractMatrix& matrix_;
    mutable VectorOfPointers<vector<P*> > points_;
    mutable typename VectorOfPointers<vector<P*> >::const_iterator current_;
    vector<CoordIndex> rowsIndex_;
    vector<CoordIndex> columnsIndex_;
};



   









template <class P> 
class TransformMatrixHandler : public MatrixHandler<P>
{
public :
      TransformMatrixHandler(const AbstractMatrix& matrix) : MatrixHandler<P>(matrix)
      {}


    double  operator()(int  i, int  j) const
    {
      int column = columnrevert_ ? maxcolumn_ - j : j + mincolumn_; 
      int row = rowrevert_ ? maxrow_ - i : i + minrow_; 
      return this->matrix_(row, column ); 
       
    }  
    
   

    int  rows() const { return maxrow_ - minrow_ +1; } 
    int  columns() const { return maxcolumn_ - mincolumn_+1; } 
    double regular_row(int i) const {
      int index = rowrevert_ ? maxrow_ - i : i + minrow_;   
       
      return this->matrix_.regular_row(index); 
    }
    double regular_column(int i) const { 
      int index = columnrevert_ ? maxcolumn_ - i : i + mincolumn_; 
      return this->matrix_.regular_column(index); 
    }
    double real_row(int i) const  { 
      int index = rowrevert_ ? maxrow_ - i : i + minrow_; 
      return this->matrix_.real_row( index); 
    }
    double real_column(int i) const { 
      int index = columnrevert_ ? maxcolumn_ - i : i + mincolumn_; 
      return this->matrix_.real_column(index); 
    }
    inline double column(int i, int j) const {  
      int column = columnrevert_ ? maxcolumn_ - j : j + mincolumn_; 
      int row = rowrevert_ ? maxrow_ - i : i + minrow_; 
      return this->matrix_.column(row, column ); 
    }
    virtual double  real_row(double row, double) const { 
         return row; 
     }
     virtual double  real_column(double, double column) const { 
         return column; 
    }
    inline double row(int i, int j) const {  
      int column = columnrevert_ ? maxcolumn_ - j : j + mincolumn_; 
        int row = rowrevert_ ? maxrow_ - i : i + minrow_; 
        return this->matrix_.row(row, column ); 
    }
    virtual bool  hasMissingValues() const {  return this->matrix_.hasMissingValues(); }
    double interpolate(double  i, double  j) const { return this->matrix_.interpolate(i, j);}
    double  missing() const { return this->matrix_.missing(); }
      int    lowerRow(double r) const {
            if ( rowsMap_.empty() ) {
                  for ( int i = 0; i < this->rows(); i++ ) 
                        rowsMap_.insert(make_pair(this->regular_row(i), i));              
            }
            
            int last = -1;
            for ( map<double, int>::const_iterator i = rowsMap_.begin(); i != rowsMap_.end(); ++i) {  
                  if ( i->first >  r  ) {                   
                        return last;
                  }                       
                  last = i->second;
            }           
            return -1;  
      }
      int    lowerColumn(double c) const { 
            if ( columnsMap_.empty() ) {
                  for ( int i = 0; i <= this->columns(); i++ ) 
                        columnsMap_.insert(make_pair(this->regular_column(i), i));              
            }
            int last = -1;
            for ( map<double, int>::const_iterator i = columnsMap_.begin(); i != columnsMap_.end(); ++i) {
                        if ( i->first > c  ) 
                              return last;
                        last = i->second;
            }
            return -1;
    } 
  
    

protected :
   int minrow_;
   int maxrow_;
   int mincolumn_;
   int maxcolumn_;
   mutable map<double, int> rowsMap_;
   mutable map<double, int> columnsMap_;
   bool rowrevert_;
   bool columnrevert_;
   
   
};
    
struct LonHelper {
      LonHelper(double min , double max, double regular, int index) : 
            min_(min), max_(max), regular_(regular), index_(index) { 
                  if ( min_ == max_ ) regular_ = min;
      }
      bool  operator<(const LonHelper& other) const
      {
              return min_ < other.min_;
      }
      double x(double xpos) {
        
        double xx = xpos;
        
            if ( same(xx, min_) ) return xx;
            if ( same(xx, max_) ) return xx;
            
        // make sure that x is in the range 
            //Log::dev() << xx << "->" << xpos << endl;
        while ( xx > min_ ) xx-=360;
            while ( true ) {
            if ( same(xx, min_) ) return xx;
                if ( same(xx, max_) ) return xx;
                  if ( min_ <= xx && xx <= max_ ) 
                return xx;
                  if ( xx > max_ ) {                        
                        return xx-360; 
                  }
                xx += 360;
            }
      }
      
      double regular() {
                  return regular_;
      }
      
      double min_;
      double max_;
      double regular_;
      int index_;
};

template <class P> 
class BoxMatrixHandler : public TransformMatrixHandler<P>
{
public:
    BoxMatrixHandler(const AbstractMatrix& matrix, const Transformation& transformation) : 
        TransformMatrixHandler<P>(matrix), 
      transformation_(transformation),
      original_(0) 
   { 
 
        
        double minx = std::min(transformation.getMinX(), transformation.getMaxX());
        double maxx = std::max(transformation.getMinX(), transformation.getMaxX());
        double miny = std::min(transformation.getMinY(), transformation.getMaxY());
        double maxy = std::max(transformation.getMinY(), transformation.getMaxY());
        
       
        
        int rows = this->matrix_.rows();
        int columns = this->matrix_.columns();
        
       
        
        this->mincolumn_ = columns-1;
        this->maxcolumn_ = 0;
        this->minrow_ = rows-1;
        this->maxrow_ = 0;
        
        for ( int row = 0; row < rows;  row++) {
            for ( int column = 0; column < columns;  column++) {
                  double x = this->matrix_.column(row, column);
                  double y = this->matrix_.row(row, column);
                  if ( minx <= x && x < maxx && miny <= y && y <= maxy) {
                        this->mincolumn_ = std::min(this->mincolumn_, column);
                        this->maxcolumn_ = std::max(this->maxcolumn_, column);
                        this->minrow_ = std::min(this->minrow_, row);
                        this->maxrow_ = std::max(this->maxrow_, row);
                  }                 
            }
            }
        
        
        if ( this->mincolumn_ > this->maxcolumn_ ) {
            this->mincolumn_ = this->maxcolumn_;
            Log::warning() << "No data to plot in the requested area" << endl;
        } 
        if ( this->minrow_ > this->maxrow_ ) {
            this->minrow_ = this->maxrow_;
            Log::warning() << "No data to plot in the requested area" << endl;
        }
        //Log::broadcast();
        
        this->mincolumn_ = std::max(this->mincolumn_-1, 0);
        this->maxcolumn_ = std::min(this->maxcolumn_+1, columns-1);  
        
        this->columnrevert_ =   this->matrix_.column(0,  this->maxcolumn_ ) < this->matrix_.column(0,  this->mincolumn_) ;
        
        this->minrow_ = std::max(this->minrow_-1, 0);
        this->maxrow_ = std::min(this->maxrow_+1, rows-1); 
       
        this->rowrevert_ =   this->matrix_.row(this->maxrow_, 0 ) < this->matrix_.row(this->minrow_, 0 ) ;
               
    }
    
    virtual const AbstractMatrix& original() const { 
      if ( !original_) 
            original_ = new BoxMatrixHandler<P>(this->matrix_.original(), this->transformation_);
      return *original_;
    }
    
    
    
  
    
    virtual ~BoxMatrixHandler() { delete original_; }
     
    // Implements the AbstractPoints interface
    virtual bool accept(double x, double y) const {
        return this->transformation_.in(x, y); 
    }
   
    double  minX() const {return std::min(this->transformation_.getMinX(), transformation_.getMaxX());  }
    double  maxX() const { return std::max(this->transformation_.getMinX(), transformation_.getMaxX()); } 
    double  minY() const { return std::min(transformation_.getMinY(), transformation_.getMaxY()); }
    double  maxY() const { return std::max(transformation_.getMinY(), transformation_.getMaxY());}

protected :
    const Transformation& transformation_;
    mutable BoxMatrixHandler<P>*   original_;
};

template <>
class BoxMatrixHandler <GeoPoint>: public TransformMatrixHandler<GeoPoint>
{
public:
    BoxMatrixHandler(const AbstractMatrix& matrix, const Transformation& transformation) : 
        TransformMatrixHandler<GeoPoint>(matrix), 
      transformation_(transformation),
      original_(0) 
   { 
 
        
        double minx = std::min(transformation.getMinX(), transformation.getMaxX()) -5 ;
        double maxx = std::max(transformation.getMinX(), transformation.getMaxX()) +5;
        double miny = std::min(transformation.getMinY(), transformation.getMaxY())-5;
        double maxy = std::max(transformation.getMinY(), transformation.getMaxY())+5;

        int rows = this->matrix_.rows();
        int columns = this->matrix_.columns();
        
        this->mincolumn_ = columns;
        this->maxcolumn_ = 0;
        this->minrow_ = rows;
        this->maxrow_ = 0;
        
        vector<pair<double, int> > all_lon;
                
        vector<pair< double, double> > longitudes_range;
        vector<pair< double, double> > latitudes_range;
        
        vector<LonHelper> longitudes;
        
        for ( int column = 0; column < columns;  column++) {
            vector<double> helper;
            for ( int row = 0; row < rows;  row++) {
                  helper.push_back(this->matrix_.column(row, column));
            }
            double min = *std::min_element(helper.begin(), helper.end());
            double max = *std::max_element(helper.begin(), helper.end());
            longitudes_range.push_back(make_pair(min, max));
        }
        
        for ( int row = 0; row < rows;  row++) {
                  vector<double> helper;
                  for ( int column = 0; column < columns;  column++)  {
                        helper.push_back(this->matrix_.row(row, column));
                  }
                  double min = *std::min_element(helper.begin(), helper.end());
                  double max = *std::max_element(helper.begin(), helper.end());
                    latitudes_range.push_back(make_pair(min, max));
          }
        
        
           //Log::dev()<< "minx=" << minx << " maxx=" << maxx << endl;
            for ( int column = 0; column < columns;  column++) {
                  
                  double min = longitudes_range[column].first;
                  double max = longitudes_range[column].second;
                  
                  //Log::dev()<< "min=" << min << " max=" << max << endl;
                  // first we make sure 
                  
                        while ( max > minx ) {
                              min -= 360;
                              max -= 360;                         
                              
            
                        }
                  
                  while ( true ) {
                        if (minx <= max && min <= maxx ) {    
                          //Log::dev()<< "new min=" << min << " new max=" << max << "[" << column << "]->" << this->matrix_.regular_column(column) << endl;                         
                              longitudes.push_back(LonHelper(min, max, this->matrix_.regular_column(column), column));  
                              
                        }
                        min += 360;
                        max += 360;
                        
                        
                        if (min > maxx) {                               
                              break;
                        }
                  }
            }
            
        
        
        std::sort(longitudes.begin(), longitudes.end());
        
        columns_ = longitudes.size();
        // Make sure that it the matrix is always  monotonic increasing 
        // This is needed by the Isoline algorithm!
        map<double, int> lats;
        for ( int row = 1; row < rows-1 ;  row++) {
            double min = latitudes_range[row].first;
            double max = latitudes_range[row].second;
                                                      
            if (miny <= max && min <= maxy ) {
                       lats.insert(make_pair(this->matrix_.regular_row(row-1), row-1));
                       lats.insert(make_pair(this->matrix_.regular_row(row), row));
                       lats.insert(make_pair(this->matrix_.regular_row(row+1), row+1));
            }
        }
     
        
        for ( map<double, int>::iterator lat = lats.begin(); lat != lats.end(); ++lat) {
             regular_latitudes_.push_back(lat->first);
             for (vector<LonHelper>::iterator lon = longitudes.begin(); lon != longitudes.end(); ++lon) {
                  
                   latitudes_.push_back(this->matrix_.row(lat->second,  lon->index_ ));
                  
                   longitudes_.push_back(lon->x(this->matrix_.column(lat->second, lon->index_)));
                 
                   values_.push_back(this->matrix_(lat->second,  lon->index_ ));
             }
        }
        for (vector<LonHelper>::iterator lon = longitudes.begin(); lon != longitudes.end(); ++lon) {
                  regular_longitudes_.push_back(lon->regular());
                
        }
                     
                  
        
     
        columns_ = regular_longitudes_.size();
        
    }
    
    virtual const AbstractMatrix& original() const { 
      if ( !original_) 
            original_ = new BoxMatrixHandler<GeoPoint>(this->matrix_.original(), this->transformation_);
      return *original_;
    }
    
    int columns() const { return regular_longitudes_.size(); }
    int rows() const { return regular_latitudes_.size(); }
    
    inline double column(int row, int column) const {
      
            
            return this->longitudes_[row*this->columns_ + column]; 
   }
    inline double row(int row, int column) const {
            
      return this->latitudes_[row*this->columns_ + column]; 
    }
     double  operator()(int  row, int  column) const {
        return this->values_[row*this->columns_ + column]; 
     }
     
     int  lowerRow(double r) const {   
       assert( !regular_latitudes_.empty());
        assert ( r >= regular_latitudes_.front() && r <= regular_latitudes_.back());
         for (unsigned int i = 0; i < regular_latitudes_.size(); i++)
                  if ( regular_latitudes_[i] >= r  ) return i-1;
        return -1;
     }  
     
     
        
     int  lowerColumn(double c) const {         
       assert( !regular_longitudes_.empty());
        assert ( c >= regular_longitudes_.front() && c <= regular_longitudes_.back());

        for (unsigned int i = 0; i < regular_longitudes_.size(); i++)
            if ( regular_longitudes_[i] >= c  ) return i-1;
        return -1;
     }
    double regular_row(int i) const { return this->regular_latitudes_[i]; }
    double regular_column(int i) const { return this->regular_longitudes_[i]; }
    
  
    
    virtual ~BoxMatrixHandler() { delete original_; }
     
    // Implements the AbstractPoints interface
    virtual bool accept(double x, double y) const {
        return this->transformation_.in(x, y); 
    }
   
    double  minX() const {return std::min(this->transformation_.getMinX(), transformation_.getMaxX());  }
    double  maxX() const { return std::max(this->transformation_.getMinX(), transformation_.getMaxX()); } 
    double  minY() const { return std::min(transformation_.getMinY(), transformation_.getMaxY()); }
    double  maxY() const { return std::max(transformation_.getMinY(), transformation_.getMaxY());}
   
    double  left() const { return regular_longitudes_.front(); }
    double bottom() const { return regular_latitudes_.front();; } 
    double  right() const { return regular_longitudes_.back(); }
    double  top() const { return regular_latitudes_.back(); }
           
protected :
    const Transformation& transformation_;
    mutable BoxMatrixHandler<GeoPoint>*   original_;
    vector<double> longitudes_;
    vector<double> latitudes_;
    vector<double> regular_latitudes_;
    vector<double> regular_longitudes_;
    vector<double> values_;
    int columns_;
};



template <class P>
class MonotonicIncreasingMatrixHandler : public MatrixHandler<P>
{
public:
    MonotonicIncreasingMatrixHandler(const AbstractMatrix& matrix) : 
        MatrixHandler<P>(matrix) { 
      // Check RowAxis...
      int row = this->matrix_.rows();
      if (this->matrix_.regular_row(1) - this->matrix_.regular_row(0) >= 0) // Increasing Axis...
           for (int i = 0; i < row; i++) {
             rows_[i] = i;
             newRowsMap_[matrix.regular_row(i)] = i;
           }
      else // Decreasing axis...
           for (int i = 0; i < row; i++) {
             rows_[i] = ( row - 1) - i;
             newRowsMap_[matrix.regular_row(( row - 1) - i)] = i;
           }
      // Check ColumnAxis
      int column = this->matrix_.columns();
      if (this->matrix_.regular_column(1) - this->matrix_.regular_column(0) >= 0) // Increasing Axis...
           for (int j = 0; j < column; j++) {
             columns_[j] = j;
             newColumnsMap_[matrix.regular_column(j)] = j;
           }
      else // Decreasing axis...
           for (int j = 0; j < column; j++) {
             columns_[j] = (column - 1) - j;
             newColumnsMap_[matrix.regular_column((column - 1) - j)] = j;
           }
      
      }
    virtual ~MonotonicIncreasingMatrixHandler() {}
    
    double  operator()(int  i, int  j) const
    {
        int x = const_cast<MonotonicIncreasingMatrixHandler*>(this)->rows_[i];
        int y = const_cast<MonotonicIncreasingMatrixHandler*>(this)->columns_[j];
        
        return this->matrix_(x, y);
    }  

    int rows() const { return this->matrix_.rows(); } 
    virtual int columns() const { return this->matrix_.columns(); } 
    virtual double regular_column(int i) const { return this->matrix_.regular_column(const_cast<MonotonicIncreasingMatrixHandler*>(this)->columns_[i]); }
    virtual double regular_row(int j) const { return this->matrix_.regular_row( const_cast<MonotonicIncreasingMatrixHandler*>(this)->rows_[j]); }
    virtual double interpolate(double  i, double  j) const {return this->matrix_.interpolate(i, j);}
    virtual double missing() const { return this->matrix_.missing(); }
     void print() 
     {
          Log::debug() << "MonotonicIncreasingMatrixHandler--->" << "\n";
          for (int j = 0; j < rows() ; j++) {
                  
            for (int i = 0; i < columns(); i++) {
                Log::dev()<< (*this)(j,i) << " ";
            }
            Log::dev()<<  "\n";
            
        }
        Log::debug() << "<-------------------" << "\n";
     }
    
    
     int  lowerRow(double r) const {    
        map<double, int>::const_iterator bound = newRowsMap_.find(r);
        if ( bound != newRowsMap_.end() ) return (*bound).second;
      
        bound = newRowsMap_.lower_bound(r);
        if ( bound == newRowsMap_.end() ) return -1; 
        return (*bound).second - 1;
   
    }  
    
    int  lowerColumn(double c) const {          
        map<double, int>::const_iterator bound = newColumnsMap_.find(c);
        if ( bound != newColumnsMap_.end() ) return (*bound).second;
            
        bound = newColumnsMap_.lower_bound(c);
        if ( bound == newColumnsMap_.end() ) return -1;
        return (*bound).second - 1;
    }
    
    protected :
    map<int, int>     rows_;
    map<int, int>     columns_;
    map<double, int>   newRowsMap_;
    map<double, int>   newColumnsMap_;
};
 

template <class P>
class OriginalMatrixHandler : public MatrixHandler<P>
{
public:
      OriginalMatrixHandler(AbstractMatrix& matrix) : 
            MatrixHandler<P>(matrix.original())
      {
      }
};

template <class P>
class ThinningMatrixHandler : public MatrixHandler<P>
{
public:
      ThinningMatrixHandler(const AbstractMatrix& matrix, int  fr, int  fc) : 
            MatrixHandler<P>(matrix), frequencyRow_(fr), frequencyColumn_(fc)
      {       

        int rows = this->matrix_.rows();
        int columns = this->matrix_.columns();
        
     
        int row = 0;
        for (int i = 0; i < rows; i+=this->frequencyRow_) {
            rowIndex_.insert(make_pair(row, i));
            row++;           
        }   
        int column=0;
        for (int i = 0; i < columns; i+=this->frequencyColumn_) {
            //Log::dev()<< "Sample --> " << column << "=" << i << endl;
                    columnIndex_.insert(make_pair(column, i));
                    //Log::dev()<< "Sample --> " << column << "=" << i << "[" << this->regular_column(column) << "]" << endl;
                    
                    column++;           
        }   
        
        columnIndex_.insert(make_pair(column, columns-1));
        //Log::dev()<< "Sample --> " << column << "=" << columns-1 << "[" << this->regular_column(column) << "]"<< endl;
      }
      int rows() const { return rowIndex_.size(); }
      int columns() const { return columnIndex_.size(); }
      
      double operator()(int row, int column) const {
      
            return this->matrix_(rowIndex(row), columnIndex(column));
      }
      
      double column(int row, int column) const {
                  return this->matrix_.column(rowIndex(row), columnIndex(column));
      }
      double row(int row, int column) const {
                  return this->matrix_.row(rowIndex(row), columnIndex(column));
      }
      double regular_row(int row) const {
                        return this->matrix_.regular_row(rowIndex(row));
            }
      double regular_column(int column) const {
                        return this->matrix_.regular_column(columnIndex(column));
            }
      double real_row(int row, int column) const {
                              return this->matrix_.real_row(rowIndex(row), columnIndex(column));
                  }
      double real_column(int row, int column) const {
                              return this->matrix_.real_column(rowIndex(row), columnIndex(column));
                  }
      

protected :
      int columnIndex(int column) const {
            map<int, int>::const_iterator index = columnIndex_.find(column);
            assert( index != columnIndex_.end() );
            return index->second;
      }
      int rowIndex(int row) const {
                  map<int, int>::const_iterator index = rowIndex_.find(row);
                  assert( index != rowIndex_.end() );
                  return index->second;
      }
      int  frequencyRow_;
      int  frequencyColumn_;
    map<int, int> rowIndex_;
    map<int, int> columnIndex_;
};

template <class P>
class CacheMatrixHandler : public MatrixHandler<P>
{    
public:
    CacheMatrixHandler(const AbstractMatrix& matrix) : 
        MatrixHandler<P>(matrix), 
        cache_(matrix.rows(), matrix.columns(), int_MIN) {
        
        cache_.missing(matrix.missing());
        rows_ = matrix.rows();
        columns_ = matrix.columns();
        for (int r = 0; r < rows_; r++) {
            cache_.rowsAxis()[r] = matrix.regular_row(r);
            rowsFinder_.add(matrix.regular_row(r), r);
        }
      for (int c = 0; c < columns_; c++) {
            cache_.columnsAxis()[c] = matrix.regular_column(c);
            columnsFinder_.add(matrix.regular_column(c), c);
      }
    std::sort(this->columnsFinder_.begin(), this->columnsFinder_.end());
    std::sort(this->rowsFinder_.begin(), this->rowsFinder_.end());
      
      cache_.setMapsAxis();
        
       
      fast_.reserve(rows_*columns_);
      
        for (int r = 0; r < rows_; r++)               
            for (int c = 0; c < columns_; c++)
                   fast_.push_back(matrix(r, c));
        
      }
       virtual void boundRow(double r, 
      double& row1, int& index1, double& row2, int& index2) const 
            { rowsFinder_.find(r, row1, index1, row2, index2); }
    virtual void boundColumn(double r, 
      double& column1, int& index1, double& column2, int& index2) const 
            { columnsFinder_.find(r, column1, index1, column2, index2); }           
    
      
    virtual inline double operator()(int  row, int  column) const {
      return fast_[(row * columns_) + column];
        
    }  
    int    lowerRow(double v) const { return cache_.lowerRow(v); }
    int    lowerColumn(double v) const { return cache_.lowerColumn(v); } 
    
    double  regular_row(int i) const { 
        return cache_.regular_row(i); 
    }
    double  regular_column(int j) const { 
        return cache_.regular_column(j); } 
    
   
    int rowIndex(double row)   const    { return  rowsFinder_.find(row); }
    int columnIndex(double column) const { return  columnsFinder_.find(column); }
    
    mutable Matrix cache_;   
    mutable vector<double> fast_; 
    int rows_;
    int columns_;
    CoordinatesFinder rowsFinder_;
    CoordinatesFinder columnsFinder_;
};

} // namespace magics
#endif

Generated by  Doxygen 1.6.0   Back to index