nnfw/include/matrices.h

Go to the documentation of this file.

/********************************************************************************
 *  Neural Network Framework.                                                   *
 *  Copyright (C) 2005-2011 Gianluca Massera <emmegian@yahoo.it>                *
 *                                                                              *
 *  This program is free software; you can redistribute it and/or modify        *
 *  it under the terms of the GNU General Public License as published by        *
 *  the Free Software Foundation; either version 2 of the License, or           *
 *  (at your option) any later version.                                         *
 *                                                                              *
 *  This program is distributed in the hope that it will be useful,             *
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of              *
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               *
 *  GNU General Public License for more details.                                *
 *                                                                              *
 *  You should have received a copy of the GNU General Public License           *
 *  along with this program; if not, write to the Free Software                 *
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA  *
 ********************************************************************************/

#ifndef MATRICES_H
#define MATRICES_H

#include "memutils.h"
#include "vectors.h"
#include <exception>

namespace farsa {

class FARSA_NNFW_API DoubleMatrix {
public:
    DoubleMatrix( unsigned int rows, unsigned int cols, bool isinternal = false ) {
        shData = new sharedData();
        shData->refcounts = 1;
        shData->nrows = rows;
        shData->ncols = cols;
        shData->tsize = rows*cols;
        shData->temporary = false;
        shData->alldata = new double[ shData->tsize ];
        shData->alldataref = new doubleRef[ shData->tsize ];
        shData->rowdata = new DoubleVector[ rows ];
        shData->coldata = new DoubleVector[ cols ];
        for( unsigned int i=0; i<rows; i++ ) {
            // this will also activate isinternal flag on DoubleVectors
            shData->rowdata[i].resizeNoData( cols );
        }
        for( unsigned int i=0; i<cols; i++ ) {
            // this will also activate isinternal flag on DoubleVectors
            shData->coldata[i].resizeNoData( rows );
        }
        for( unsigned int t=0; t<shData->tsize; t++ ) {
            unsigned int r = t/cols; //--- division may be expensive
            unsigned int c = t%cols; //--- module may be expensive
            shData->rowdata[r].shData->dataref[c].setRef( shData->alldata + t );
            shData->coldata[c].shData->dataref[r].setRef( shData->alldata + t );
            shData->alldataref[t].setRef( shData->alldata + t );
        }
        this->isinternal = isinternal;
    };
    DoubleMatrix( const DoubleMatrix& src ) {
        shData = src.shData;
        shData->refcounts++;
        //--- is not a temporary anymore !
        shData->temporary = false;
        isinternal = false;
    };
    DoubleMatrix(  unsigned int rows, unsigned int cols, bool temp, unsigned int dummy ) {
        (void)dummy; // This stops warnings. Assigning dummy a value would have produce a "parameter ‘dummy’ set but not used" warning
        shData = new sharedData();
        shData->refcounts = 1;
        shData->nrows = rows;
        shData->ncols = cols;
        shData->tsize = rows*cols;
        shData->temporary = temp;
        shData->alldata = new double[ shData->tsize ];
        shData->alldataref = new doubleRef[ shData->tsize ];
        shData->rowdata = new DoubleVector[ rows ];
        shData->coldata = new DoubleVector[ cols ];
        for( unsigned int i=0; i<rows; i++ ) {
            // this will also activate isinternal flag on DoubleVectors
            shData->rowdata[i].resizeNoData( cols );
        }
        for( unsigned int i=0; i<cols; i++ ) {
            // this will also activate isinternal flag on DoubleVectors
            shData->coldata[i].resizeNoData( rows );
        }
        for( unsigned int t=0; t<shData->tsize; t++ ) {
            unsigned int r = t/cols; //--- division may be expensive
            unsigned int c = t%cols; //--- module may be expensive
            shData->rowdata[r].shData->dataref[c].setRef( shData->alldata + t );
            shData->coldata[c].shData->dataref[r].setRef( shData->alldata + t );
            shData->alldataref[t].setRef( shData->alldata + t );
        }
        this->isinternal = false;
    };
    ~DoubleMatrix() {
        shData->refcounts -= 1;
        if ( shData->refcounts == 0 ) {
            //--- the last destroy the data
            if ( shData->tsize > 0 ) {
                delete [](shData->rowdata);
                delete [](shData->coldata);
                delete [](shData->alldataref);
                delete [](shData->alldata);
            }
            delete shData;
        }
    };
    unsigned int rows() const {
        return shData->nrows;
    };
    unsigned int cols() const {
        return shData->ncols;
    };
    unsigned int size() const {
        return shData->tsize;
    };
    bool operator==( const DoubleMatrix& right ) const {
        if ( shData == right.shData ) return true;
        if ( shData->tsize != right.shData->tsize ) return false;
        if ( shData->nrows != right.shData->nrows ) return false;
        if ( shData->ncols != right.shData->ncols ) return false;
        for( unsigned int i=0; i<shData->tsize; i++ ) {
            if ( shData->alldata[i] != right.shData->alldata[i] ) return false;
        }
        return true;
    };
    bool operator!=( const DoubleMatrix& right ) const {
        return !( *this == right );
    };
    DoubleMatrix& resize( unsigned int rows, unsigned int cols ) {
        if ( isinternal ) throw MatrixResizeNotAllowed();
        DoubleMatrix newMatrix( rows, cols );
        newMatrix.copyValues( *this );
        //--- the assignment will release the shared data, so there is no need to call detach()
        (*this) = newMatrix;
        return (*this);
    };
    DoubleMatrix& operator=( const DoubleMatrix& src ) {
        if ( isinternal ) throw MatrixAssignmentNotAllowed();
        if ( shData == src.shData ) return (*this);
        //--- eliminate the previous data
        shData->refcounts -= 1;
        if ( shData->refcounts == 0 ) {
            //--- the last destroy the data
            if ( shData->tsize > 0 ) {
                delete [](shData->rowdata);
                delete [](shData->coldata);
                delete [](shData->alldata);
            }
            delete shData;
        }
        //--- set the new data taken from src
        shData = src.shData;
        shData->refcounts++;
        //--- is not a temporary anymore !
        shData->temporary = false;
        return (*this);
    };
    DoubleMatrix& copyValues( const DoubleMatrix& src ) {
        // if shared data is the same, they already share the same values
        if ( shData == src.shData ) return (*this);
        detach();
        unsigned int rowMax = qMin( shData->nrows, src.shData->nrows );
        unsigned int colMax = qMin( shData->ncols, src.shData->ncols );
        for( unsigned int r=0; r<rowMax; r++ ) {
            for( unsigned int c=0; c<colMax; c++ ) {
                shData->rowdata[r][c] = src.shData->rowdata[r][c];
            }
        }
        return (*this);
    };
    DoubleMatrix& copyValues( const DoubleMatrix& src, unsigned int srcRowOffset, unsigned int srcColOffset, unsigned int thisRowOffset, unsigned int thisColOffset ) {
        // if shared data is the same and offsets are zero, then there is nothing to do
        if ( shData == src.shData && srcRowOffset == 0 && srcColOffset == 0 && thisRowOffset == 0 && thisColOffset == 0 ) {
            return (*this);
        }
        detach();
        unsigned int rowMax = qMin( shData->nrows-thisRowOffset, src.shData->nrows-srcRowOffset );
        unsigned int colMax = qMin( shData->ncols-thisColOffset, src.shData->ncols-thisColOffset );
        for( unsigned int r=0; r<rowMax; r++ ) {
            for( unsigned int c=0; c<colMax; c++ ) {
                shData->rowdata[r+thisRowOffset][c+thisColOffset] = src.shData->rowdata[r+srcRowOffset][c+srcColOffset];
            }
        }
        return (*this);
    };
    DoubleMatrix& copyValues( const DoubleMatrix& src, unsigned int srcRowOffset, unsigned int srcColOffset, unsigned int thisRowOffset, unsigned int thisColOffset, unsigned int rowStride, unsigned int colStride ) {
        if ( rowStride == 0 || colStride == 0 ) return (*this);
        // if shared data is the same and offsets are zero, then there is nothing to do
        if ( shData == src.shData && srcRowOffset == 0 && srcColOffset == 0 && thisRowOffset == 0 && thisColOffset == 0 ) {
            return (*this);
        }
        detach();
        unsigned int rowMax = qMin( shData->nrows-thisRowOffset, src.shData->nrows-srcRowOffset );
        unsigned int colMax = qMin( shData->ncols-thisColOffset, src.shData->ncols-thisColOffset );
        for( unsigned int r=0; r<rowMax; r=r+rowStride ) {
            for( unsigned int c=0; c<colMax; c=c+colStride ) {
                shData->rowdata[r+thisRowOffset][c+thisColOffset] = src.shData->rowdata[r+srcRowOffset][c+srcColOffset];
            }
        }
        return (*this);
    };
    DoubleMatrix& detach() {
        if ( shData->refcounts > 1 ) {
            sharedData* tmp = new sharedData();
            tmp->refcounts = 1;
            tmp->nrows = shData->nrows;
            tmp->ncols = shData->ncols;
            tmp->tsize = shData->tsize;
            tmp->temporary = false;
            tmp->alldata = new double[shData->tsize];
            tmp->rowdata = new DoubleVector[ shData->nrows ];
            tmp->coldata = new DoubleVector[ shData->ncols ];
            for( unsigned int i=0; i<shData->nrows; i++ ) {
                tmp->rowdata[i].resizeNoData( shData->ncols );
            }
            for( unsigned int i=0; i<shData->ncols; i++ ) {
                tmp->coldata[i].resizeNoData( shData->nrows );
            }
            for( unsigned int t=0; t<tmp->tsize; t++ ) {
                unsigned int r = t/( tmp->ncols ); //--- division may be expensive
                unsigned int c = t%( tmp->ncols ); //--- module may be expensive
                tmp->rowdata[r].shData->dataref[c].setRef( tmp->alldata + t );
                tmp->coldata[c].shData->dataref[r].setRef( tmp->alldata + t );
            }
            shData->refcounts--;
            shData = tmp;
        }
        return (*this);
    };
    doubleRef& at( unsigned int row, unsigned int col ) {
#ifdef FARSA_DEBUG
        if ( row >= shData->nrows ) throw OutsideMatrixBoundaries("Accessing an element beyond Row boundary of matrix");
        if ( col >= shData->ncols ) throw OutsideMatrixBoundaries("Accessing an element beyond Column boundary of matrix");
#endif
        detach();
        return shData->rowdata[row][col];
    };
    double at( unsigned int row, unsigned int col ) const {
#ifdef FARSA_DEBUG
        if ( row >= shData->nrows ) throw OutsideMatrixBoundaries("Accessing an element beyond Row boundary of matrix");
        if ( col >= shData->ncols ) throw OutsideMatrixBoundaries("Accessing an element beyond Column boundary of matrix");
#endif
        return shData->rowdata[row][col];
    };
    DoubleVector& operator[]( unsigned int row ) {
#ifdef FARSA_DEBUG
        if ( row >= shData->nrows ) throw OutsideMatrixBoundaries("Accessing an element beyond Row boundary of matrix");
#endif
        detach();
        return shData->rowdata[row];
    };
    DoubleVector operator[]( unsigned int row ) const {
#ifdef FARSA_DEBUG
        if ( row >= shData->nrows ) throw OutsideMatrixBoundaries("Accessing an element beyond Row boundary of matrix");
#endif
        return shData->rowdata[row];
    };
    DoubleVector& row( unsigned int r ) {
#ifdef FARSA_DEBUG
        if ( r >= shData->nrows ) throw OutsideMatrixBoundaries("Accessing an element beyond Row boundary of matrix");
#endif
        detach();
        return shData->rowdata[r];
    };
    DoubleVector row( unsigned int r ) const {
#ifdef FARSA_DEBUG
        if ( r >= shData->nrows ) throw OutsideMatrixBoundaries("Accessing an element beyond Row boundary of matrix");
#endif
        return shData->rowdata[r];
    };
    DoubleVector& column( unsigned int c ) {
#ifdef FARSA_DEBUG
        if ( c >= shData->ncols ) throw OutsideMatrixBoundaries("Accessing an element beyond Column boundary of matrix");
#endif
        detach();
        return shData->coldata[c];
    };
    DoubleVector column( unsigned int c ) const {
#ifdef FARSA_DEBUG
        if ( c >= shData->ncols ) throw OutsideMatrixBoundaries("Accessing an element beyond Column boundary of matrix");
#endif
        return shData->coldata[c];
    };
    DoubleMatrix& steady( unsigned int i, unsigned int j ) {
#ifdef FARSA_DEBUG
        if( i >= shData->nrows || j >= shData->ncols ) throw OutsideMatrixBoundaries("Fixating elements outside boundary");
#endif
        detach();
        shData->rowdata[i].steady(j);
        shData->coldata[j].steady(i);
        shData->alldataref[ i*shData->ncols + j ].setSteady();
        return (*this);
    };
    DoubleMatrix& unsteady( unsigned int i, unsigned int j ) {
#ifdef FARSA_DEBUG
        if( i >= shData->nrows || j >= shData->ncols ) throw OutsideMatrixBoundaries("Un-Fixating elements outside boundary");
#endif
        detach();
        shData->rowdata[i].unsteady(j);
        shData->coldata[j].unsteady(i);
        shData->alldataref[ i*shData->ncols + j ].setNoSteady();
        return (*this);
    };
    bool isSteady( unsigned int i, unsigned int j ) const {
#ifdef FARSA_DEBUG
        if( i >= shData->nrows || j >= shData->ncols ) throw OutsideMatrixBoundaries("Accessing elements outside boundary");
#endif
        return shData->rowdata[i].isSteady(j);
    };
    DoubleMatrix& setAll( const double value ) {
        detach();
        for( unsigned int i=0; i<shData->tsize; i++ ) {
            shData->alldataref[i] = value;
        }
        return (*this);
    };
    DoubleMatrix& setIdentity() {
        detach();
        for( unsigned int r=0; r<shData->nrows; r++ ) {
            for( unsigned int c=0; c<shData->ncols; c++ ) {
                shData->rowdata[r][c] = (r==c);
            }
        }
        return (*this);
    };
    DoubleMatrix& zeroing() {
        detach();
        for( unsigned int i=0; i<shData->tsize; i++ ) {
            shData->alldataref[i] = 0;
        }
        return (*this);
    };
    const DoubleMatrix operator+( const DoubleMatrix& right ) const {
#ifdef FARSA_DEBUG
        if ( shData->tsize != right.shData->tsize || shData->nrows != right.shData->nrows || shData->ncols != right.shData->ncols ) {
            throw IncompatibleMatrices("Incompatibles Matrices in operator+ (dimension must be equals)");
        }
#endif
        if ( right.shData->temporary ) {
            for( unsigned int i=0; i<shData->tsize; i++ ) {
                right.shData->alldataref[i] = shData->alldata[i] + right.shData->alldata[i];
            }
            return right;
        } else if ( shData->temporary ) {
            for( unsigned int i=0; i<shData->tsize; i++ ) {
                shData->alldataref[i] = shData->alldata[i] + right.shData->alldata[i];
            }
            return (*this);
        } else {
            DoubleMatrix ret( shData->nrows, shData->ncols, true, 0 );
            for( unsigned int i=0; i<shData->tsize; i++ ) {
                ret.shData->alldataref[i] = shData->alldata[i] + right.shData->alldata[i];
            }
            return ret;
        }
    };
    const DoubleMatrix operator-( const DoubleMatrix& right ) const {
#ifdef FARSA_DEBUG
        if ( shData->tsize != right.shData->tsize || shData->nrows != right.shData->nrows || shData->ncols != right.shData->ncols ) {
            throw IncompatibleMatrices("Incompatibles Matrices in operator- (dimension must be equals)");
        }
#endif
        if ( right.shData->temporary ) {
            for( unsigned int i=0; i<shData->tsize; i++ ) {
                right.shData->alldataref[i] = shData->alldata[i] - right.shData->alldata[i];
            }
            return right;
        } else if ( shData->temporary ) {
            for( unsigned int i=0; i<shData->tsize; i++ ) {
                shData->alldataref[i] = shData->alldata[i] - right.shData->alldata[i];
            }
            return (*this);
        } else {
            DoubleMatrix ret( shData->nrows, shData->ncols, true, 0 );
            for( unsigned int i=0; i<shData->tsize; i++ ) {
                ret.shData->alldataref[i] = shData->alldata[i] - right.shData->alldata[i];
            }
            return ret;
        }
    };
    const DoubleMatrix operator*( const DoubleMatrix& right ) const {
#ifdef FARSA_DEBUG
        if ( shData->ncols != right.shData->nrows ) {
            throw IncompatibleMatrices("Incompatibles Matrices in operator* (left matrix column dimension must be equal to right matrix row dimension)");
        }
#endif
        // FIXME: this is very slow and inefficient way to multiple matrices
        DoubleMatrix ret( shData->nrows, right.shData->ncols, true, 0 );
        ret.zeroing();
        for( unsigned int i=0; i<ret.shData->nrows; i++ ) {
            for( unsigned int j=0; j<ret.shData->ncols; j++ ) {
                for( unsigned int r=0; r<shData->ncols; r++ ) {
                    ret.shData->rowdata[i][j] += shData->rowdata[i][r] * right.shData->rowdata[r][j];
                }
            }
        }
        return ret;
    };
    const DoubleMatrix operator%( const DoubleMatrix& right ) const {
#ifdef FARSA_DEBUG
        if ( shData->tsize != right.shData->tsize || shData->nrows != right.shData->nrows || shData->ncols != right.shData->ncols ) {
            throw IncompatibleMatrices("Incompatibles Matrices in operator% (dimension must be equals)");
        }
#endif
        if ( right.shData->temporary ) {
            for( unsigned int i=0; i<shData->tsize; i++ ) {
                right.shData->alldataref[i] = shData->alldata[i] * right.shData->alldata[i];
            }
            return right;
        } else if ( shData->temporary ) {
            for( unsigned int i=0; i<shData->tsize; i++ ) {
                shData->alldataref[i] = shData->alldata[i] * right.shData->alldata[i];
            }
            return (*this);
        } else {
            DoubleMatrix ret( shData->nrows, shData->ncols, true, 0 );
            for( unsigned int i=0; i<shData->tsize; i++ ) {
                ret.shData->alldataref[i] = shData->alldata[i] * right.shData->alldata[i];
            }
            return ret;
        }
    };
    const DoubleMatrix operator/( const DoubleMatrix& right ) const {
#ifdef FARSA_DEBUG
        if ( shData->tsize != right.shData->tsize || shData->nrows != right.shData->nrows || shData->ncols != right.shData->ncols ) {
            throw IncompatibleMatrices("Incompatibles Matrices in operator/ (dimension must be equals)");
        }
#endif
        if ( right.shData->temporary ) {
            for( unsigned int i=0; i<shData->tsize; i++ ) {
                right.shData->alldataref[i] = shData->alldata[i] / right.shData->alldata[i];
            }
            return right;
        } else if ( shData->temporary ) {
            for( unsigned int i=0; i<shData->tsize; i++ ) {
                shData->alldataref[i] = shData->alldata[i] / right.shData->alldata[i];
            }
            return (*this);
        } else {
            DoubleMatrix ret( shData->nrows, shData->ncols, true, 0 );
            for( unsigned int i=0; i<shData->tsize; i++ ) {
                ret.shData->alldataref[i] = shData->alldata[i] / right.shData->alldata[i];
            }
            return ret;
        }
    };
    DoubleMatrix& operator+=( const DoubleMatrix& right ) {
#ifdef FARSA_DEBUG
        if ( shData->tsize != right.shData->tsize || shData->nrows != right.shData->nrows || shData->ncols != right.shData->ncols ) {
            throw IncompatibleMatrices("Incompatibles Matrices in operator += (dimension must be equals)");
        }
#endif
        detach();
        for( unsigned int i=0; i<shData->tsize; i++ ) {
            shData->alldataref[i] += right.shData->alldata[i];
        }
        return (*this);
    };
    DoubleMatrix& operator-=( const DoubleMatrix& right ) {
#ifdef FARSA_DEBUG
        if ( shData->tsize != right.shData->tsize || shData->nrows != right.shData->nrows || shData->ncols != right.shData->ncols ) {
            throw IncompatibleMatrices("Incompatibles Matrices in operator -= (dimension must be equals)");
        }
#endif
        detach();
        for( unsigned int i=0; i<shData->tsize; i++ ) {
            shData->alldataref[i] -= right.shData->alldata[i];
        }
        return (*this);
    };
    DoubleMatrix& operator%=( const DoubleMatrix& right ) {
#ifdef FARSA_DEBUG
        if ( shData->tsize != right.shData->tsize || shData->nrows != right.shData->nrows || shData->ncols != right.shData->ncols ) {
            throw IncompatibleMatrices("Incompatibles Matrices in operator %= (dimension must be equals)");
        }
#endif
        detach();
        for( unsigned int i=0; i<shData->tsize; i++ ) {
            shData->alldataref[i] *= right.shData->alldata[i];
        }
        return (*this);
    };
    DoubleMatrix& operator/=( const DoubleMatrix& right ) {
#ifdef FARSA_DEBUG
        if ( shData->tsize != right.shData->tsize || shData->nrows != right.shData->nrows || shData->ncols != right.shData->ncols ) {
            throw IncompatibleMatrices("Incompatibles Matrices in operator /= (dimension must be equals)");
        }
#endif
        detach();
        for( unsigned int i=0; i<shData->tsize; i++ ) {
            shData->alldataref[i] /= right.shData->alldata[i];
        }
        return (*this);
    };

private:
    class sharedData {
    public:
        unsigned int nrows;
        unsigned int ncols;
        unsigned int tsize;
        double* alldata;
        doubleRef* alldataref;
        DoubleVector* rowdata;
        DoubleVector* coldata;
        int refcounts;
        bool temporary;
    };
    sharedData* shData;
    bool isinternal;
};

}

#endif