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

Matrix.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 Matrix.h
    \brief Definition of the Template class Matrix.
    
    Magics Team - ECMWF 2004
    
    Started: Wed 18-Feb-2004
    
    Changes:
    
*/

#ifndef Matrix_H
#define Matrix_H

#include "magics.h"

namespace magics {

class XmlNode;

class AbstractMatrix 
{
public :
    virtual ~AbstractMatrix(){}
    virtual double  operator()(int  i, int  j) const = 0;   
    virtual int    rows() const = 0;
    virtual int    columns() const = 0;
    virtual double  regular_row(int) const = 0;
    virtual double  regular_column(int) const = 0;
    virtual double  row(int, int) const = 0;
    virtual double  column(int, int) const = 0;
     virtual double  real_row(double, double) const = 0;
    virtual double  real_column(double, double) const = 0;
    virtual int    lowerRow(double) const = 0;
    virtual int    lowerColumn(double) const = 0;    
    virtual double  interpolate(double  i, double  j) const = 0;   
    virtual double  missing() const = 0;
    virtual double  XResolution() const = 0;
    virtual double  YResolution() const = 0;
    virtual double  width() const = 0;
    virtual double   height() const = 0;
    virtual bool akimaEnable() const { return false; }
  

   
   

    
     virtual const AbstractMatrix&  original() const { return *this; }
     
     virtual int firstRow() const = 0;
     virtual int nextRow(int, int) const = 0;
     
     virtual int firstColumn() const = 0;
     virtual int nextColumn(int, int) const = 0;
     
     template <class O> 
     void for_each(int xf, int yf, const O& object)
     {
         int i,j; // STEPHAN: I assume that is right ... ???
       for ( i = firstRow(); i > 0; i = nextRow(i, xf) )
            for ( j = firstColumn(); j > 0; j = nextColumn(j, yf) )
                  object(row(i,j), column(i,j), (*this)(i, j));
     }

    virtual double  minX() const = 0;
    virtual double  minY() const = 0;
    virtual double  maxX() const = 0;
    virtual double  maxY() const = 0;
    
    virtual double  left() const = 0;
       virtual double  top() const = 0;
       virtual double  right() const = 0;
       virtual double  bottom() const = 0;
       
       virtual double  x(double, double) const = 0;
       virtual double  y(double, double) const = 0;      
       
    virtual int rowIndex(double r) const  = 0;
    virtual int columnIndex(double c) const = 0;
    
    virtual void boundRow(double r, 
      double& row1, int& index1, double& row2, int& index2) const = 0;
    
    virtual void boundColumn(double r, 
      double& column1, int& index1, double& column2, int& index2) const = 0;
    
    virtual bool accept(double, double) const { return true; }
   
    virtual vector<double>&  rowsAxis() const = 0;
    
    virtual bool  hasMissingValues() const  { return false; }
  
    
    virtual vector<double>&  columnsAxis() const = 0;
    virtual void print(ostream& out) const 
        { out << "No Print implemented for this MatrixHandler" << "\n"; }
    //! Overloaded << operator to call print().
      friend ostream& operator<<(ostream& s,const AbstractMatrix& p)
            { p.print(s); return s; }
};

struct Plus 
{
    Plus(double offset, double missing) : offset_(offset), missing_(missing) {}
    double operator()(double x) const { return ( x == missing_ ) ? missing_ : x + offset_; }
    double offset_;
    double missing_;
};

struct Multiply 
{
    Multiply(double factor, double missing) : factor_(factor), missing_(missing) {}
    double operator()(double x) const { return ( x == missing_ ) ? missing_ : x * factor_; }
    double factor_;
    double missing_;
};


class OutOfRange : public MagicsException
{
public:
    
    OutOfRange(double r, double c) 
    { 
        ostringstream s;
        s << "Out of Range: Cannot access [" << r << ", " << c << "]" << ends;
        what_ = s.str();
    }
    OutOfRange(double x) 
    { 
        ostringstream s;
        s << "Out of Range: Cannot access [" << x << "]" << ends;
        what_ = s.str();
    }
};



class Matrix: public AbstractMatrix, public magvector<double> {

public:
      Matrix(int rows, int columns): 
            rows_(rows), 
            columns_(columns),
            missing_(int_MIN),
            akima_(false)
            
      {  
         reserve(rows_*columns_); 
         
      }
      
      Matrix* clone() { return new Matrix(); }
      void set(const XmlNode&) { }
      
    Matrix(int rows, int columns, double val):     
        rows_(rows), 
        columns_(columns), missing_(int_MIN), akima_(false)
     {       
         
         resize(rows_ * columns_, val); 
         rowsAxis_.resize(rows_, val);
         columnsAxis_.resize(columns_, val); 
      }
    
    Matrix(): missing_(int_MIN), akima_(false) {}
    
    void set(int rows, int columns) 
    {
         rows_ = rows;
         columns_ = columns;
         reserve(rows_*columns_); 
         rowsAxis_.reserve(rows); 
         columnsAxis_.reserve(columns);  
    }
    
    
    
    virtual ~Matrix() {}
    
    double width() const { return regular_column(columns_ - 1) - regular_column(0); }
    double height() const { return regular_row(rows_ - 1) - regular_row(0); }
    
    int rows() const { return rows_; }
    int columns() const { return columns_; }
     
    double regular_row(int i) const { return rowsAxis_[i]; } 
    double row(int i, int) const { return regular_row(i); }
    double real_row(double r, double) const { return r; }
    double real_column(double, double c) const { return c; }
    
    int  lowerRow(double r) const {    
            map<double, int>::const_iterator bound = rowsMap_.find(r);
            if ( bound != rowsMap_.end() ) return (*bound).second;
      
            bound = rowsMap_.lower_bound(r);
        if ( bound == rowsMap_.end() ) return -1;
            return (*bound).second - 1;
   
    }  
      
  
   
    int  lowerColumn(double c) const {    
            map<double, int>::const_iterator bound = columnsMap_.find(c);
            if ( bound != columnsMap_.end() ) return (*bound).second;
            
            bound = columnsMap_.lower_bound(c);
              if ( bound == columnsMap_.end() ) return -1;  
            return (*bound).second - 1;
    }
    
    double regular_column(int j) const { return columnsAxis_[j];  }
    double column(int, int j) const { return columnsAxis_[j];  }
     
    void missing(double missing) { missing_ = missing; }
    double missing() const { return missing_; }
    
    void setRowsAxis(const vector<double>& axis) 
    {
        int ind = 0;
        rowsAxis_.reserve(axis.size());
        for (vector<double>::const_iterator val = axis.begin(); val != axis.end(); ++val) {
            rowsAxis_.push_back(*val);
            rowsMap_[*val] = ind++;
        }
        rows_ = axis.size();          
    }
    void setColumnsAxis(const vector<double>& axis) 
    {
        int ind = 0;
        columnsAxis_.reserve(axis.size());
        for ( vector<double>::const_iterator val = axis.begin(); val != axis.end(); ++val) {
            columnsAxis_.push_back(*val);
            columnsMap_[*val] = ind++;
        }
        columns_ = axis.size();
            
    }
    
    virtual void setMapsAxis() 
    {
        int ind = 0;
        for (vector<double>::const_iterator val = rowsAxis_.begin(); val != rowsAxis_.end(); ++val) {
            rowsMap_[*val] = ind++;
        }
        rows_ = ind;
        
        
        ind = 0;
        for (vector<double>::const_iterator val = columnsAxis_.begin(); val != columnsAxis_.end(); ++val) {
            columnsMap_[*val] = ind++;
        }
        columns_ = ind;
    }
    
  
    double interpolate(double r, double c) const 
    {   
        int ii, jj;
        //Log::debug() << " interpolate---> " << r << " , " << c << "\n";
       
            ii = row_ind(r);
        
       if (ii == -1) {
           
            map<double, int>::const_iterator bound = rowsMap_.lower_bound(r);
            
            if ( bound == rowsMap_.end() ) return missing_;
            double xr = (*bound).first;
            if (bound == rowsMap_.begin()) 
                ii = (*bound).second;
                
            else {
                double xl = (*--bound).first;
          
           //Log::debug() << xl << " < " << r << " < " << xr << "\n";
            
            double a = (*this).interpolate(xl, c);
            double b = (*this).interpolate(xr, c);
          
            if ( same(a, missing_) || same(b, missing_) ) return missing_;
            
            double da = (xr-r)/(xr-xl);
            double db = (r-xl)/(xr-xl);
            double val = (a*da) + (b*db);
            //Log::debug() << a << "*" << da << "+" << b << "*" << db << "=" << val << "\n";
            return val;
            }
        }    
       
        jj = column_ind(c);
        
        if ( jj == -1 ) {
             map<double, int>::const_iterator bound = columnsMap_.lower_bound(c);
            if ( bound == columnsMap_.end() ) return missing_;
            
            double yr = (*bound).first;
          
            
            if (bound == columnsMap_.begin()) 
                jj = (*bound).second;
            else {
            double yl = (*--bound).first;
            //Log::debug() << yl << " < " << c << " < " << yr << "\n";
            double a = (*this)(ii, column_ind(yl));           
            double b = (*this)(ii, column_ind(yr));
            
            if ( same(a, missing_) || same(b, missing_) ) return missing_;
            
            double da = ((yr-c)/(yr-yl));
            double db = (c-yl)/(yr-yl);
            double val = (a*da) + (b*db);
            //Log::debug() << a << "*" << da << "+" << b << "*" << db << "=" << val << "\n";
            return val;
            }
            
        }    
      
        return (*this)(ii, jj);
    }
    
   void multiply(double factor) 
    {
        if (factor == 1 ) return;
        std::transform(begin(), end(), begin(), Multiply(factor, missing_));
    }
    
    void plus(double offset) 
    {
        if (offset == 0 ) return;
        std::transform(begin(), end(), begin(), Plus(offset, missing_));
       
    }
    
     virtual int firstRow() const { return 0; }
     virtual int nextRow(int i, int f) const   
     { 
      i += f; 
      return ( i < rows_ ) ? i : -1;
     }
     
     virtual int firstColumn() const { return 0; }
     virtual int nextColumn(int j, int f) const   
     { 
      j += f; 
      return ( j < rows_ ) ? j : -1;
     }
     
    
    double operator()(int row, int column) const
    { 
           // Here for perfrormance we are trusting the user we do not catch exception 
                        return (*this)[ row * columns_ + column];
          
    }
    
    double YResolution() const {
           magvector<double> diff;
           diff.reserve(rowsAxis_.size());
           std::adjacent_difference(rowsAxis_.begin(), rowsAxis_.end(), back_inserter(diff));
           double resol = std::accumulate(diff.begin()+1, diff.end(), 0.)/(diff.size()-1);
           //Log::dev() << "Matrix::YResolution()--->" << resol << "\n";
           return resol;
    }
     double XResolution() const {
           magvector<double> diff;
           diff.reserve(columnsAxis_.size());
           std::adjacent_difference(columnsAxis_.begin(), columnsAxis_.end(), back_inserter(diff));
           double resol = std::accumulate(diff.begin()+1, diff.end(), 0.)/(diff.size()-1);
           //Log::dev() << "Matrix::XResolution()--->" << resol << "\n";
           return resol;
    }
   
    vector<double>& rowsAxis() const { return rowsAxis_; }
    vector<double>& columnsAxis() const  { return columnsAxis_; }
    
    
    double  minX() const { return columnsAxis_.front(); }
    double  minY() const { return rowsAxis_.front(); } 
    double  maxX() const { return columnsAxis_.back(); }
    double  maxY() const { return rowsAxis_.back(); }
    
    double  left() const { return columnsAxis_.front(); }
      double bottom() const { return rowsAxis_.front(); } 
      double  right() const { return columnsAxis_.back(); }
      double  top() const { return rowsAxis_.back(); }
      
      double x(double x, double) const  { return x; }
      double y(double, double y) const { return y; }
    
    virtual int rowIndex(double r) const     { return row_ind(r); } 
    virtual int columnIndex(double c) const  { return column_ind(c); } 
    virtual bool akimaEnable() const  { return akima_; }  
    void akimaEnabled()  { akima_ = true; }  
    void akimaDisabled()  { akima_ = false; } 
    virtual void boundRow(double r, 
      double& row1, int& index1, double& row2, int& index2) const { 
      
      map<double, int>::const_iterator bound = rowsMap_.lower_bound(r);
      row1 = bound->first;
      index1 = bound->second;
      bound++;
      row2 = bound->first;
      index2 = bound->second;
    } 
    
    virtual void boundColumn(double r, 
      double& column1, int& index1, double& column2, int& index2) const { 
      
      map<double, int>::const_iterator bound = columnsMap_.lower_bound(r);
      column1 = bound->first;
      index1 = bound->second;
      bound++;
      column2 = bound->first;
      index2 = bound->second;
    } 
    
    
    
protected:
     //! Method to print string about this class on to a stream of type ostream (virtual).
       virtual void print(ostream& out) const {
            out << "Matrix<P>[";
      out << "rowsAxis=" << rowsAxis_;
      out << ", columnsAxis=" << columnsAxis_;
            out << ", values=";
      magvector<double>::print(out);
      out << "]"; 
       }
       
     map<double, int>   rowsMap_;
     mutable magvector<double>     rowsAxis_;
     
     map<double, int>    columnsMap_;
     mutable magvector<double>      columnsAxis_;
     
     int rows_;
     int columns_;
     double missing_;
     bool akima_;

    int row_ind(double row) const {
        map<double, int>::const_iterator i = rowsMap_.find(row);
        if ( i != rowsMap_.end() ) return (*i).second;
        return -1;
    }
    int column_ind(double column) const {
        map<double, int>::const_iterator j = columnsMap_.find(column);
        if ( j != columnsMap_.end() ) return (*j).second;
        return -1;
    }    
private:
   
    
// -- Friends
    //! Overloaded << operator to call print().
      friend ostream& operator<<(ostream& s,const Matrix& p)
            { p.print(s); return s; }

};


class ProjectedMatrix: public Matrix
{
public:
          
            ProjectedMatrix(int rows, int columns): Matrix( rows, columns) {
                   rows_ = rows;
                           columns_ = columns;
                           reserve(rows_*columns_); 
                           rowsAxis_.reserve(rows);  
                           rowsArray_.reserve(rows_*columns_);  
                           columnsArray_.reserve(rows_*columns_);  
                           columnsAxis_.reserve(columns);  
            }
      
            double row(int row, int column) const {
                   return rowsArray_[ row * columns_ + column];
            }
            
            double column(int row, int column) const {
                               return columnsArray_[ row * columns_ + column];
                        }
              double real_column(int row, int column) const { 
                    return columnsArray_[ row * columns_ + column];
                  
              }
              double real_row(int row,  int column) const { 
                    return rowsArray_[ row * columns_ + column];
               }
              vector<double>&  rowsArray() const { return rowsArray_; }
              vector<double>&  columnsArray() const { return columnsArray_; }
              void setMapsAxis() 
                {
                    
                    for (int i = 0; i < rows_; i++) {
                        rowsAxis_.push_back(i);
                        rowsMap_[i] = i;
                    }
                    for (int i = 0; i < columns_; i++) {
                        columnsAxis_.push_back(i);
                        columnsMap_[i] = i;
                   }
                    
                    
                   
                }
protected:
              mutable vector<double> rowsArray_;
              mutable vector<double> columnsArray_; 
             
                  
};

} // namespace magics

#endif

Generated by  Doxygen 1.6.0   Back to index