SGMatrix.h

Go to the documentation of this file.

/*
 * 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 3 of the License, or
 * (at your option) any later version.
 *
 * Written (W) 2012 Fernando José Iglesias García
 * Written (W) 2010,2012 Soeren Sonnenburg
 * Copyright (C) 2010 Berlin Institute of Technology
 * Copyright (C) 2012 Soeren Sonnenburg
 */
#ifndef __SGMATRIX_H__
#define __SGMATRIX_H__

#include <shogun/lib/config.h>
#include <shogun/lib/DataType.h>
#include <shogun/lib/SGReferencedData.h>

namespace shogun
{
    template<class T> class SGVector;
    template<class T> class SGMatrixList;
template<class T> class SGMatrix : public SGReferencedData
{
    public:
        SGMatrix() : SGReferencedData()
        {
            init_data();
        }

        SGMatrix(T* m, index_t nrows, index_t ncols, bool ref_counting=true)
            : SGReferencedData(ref_counting), matrix(m),
            num_rows(nrows), num_cols(ncols) { }

        SGMatrix(index_t nrows, index_t ncols, bool ref_counting=true)
            : SGReferencedData(ref_counting), num_rows(nrows), num_cols(ncols)
        {
            matrix=SG_MALLOC(T, ((int64_t) nrows)*ncols);
        }

        SGMatrix(const SGMatrix &orig) : SGReferencedData(orig)
        {
            copy_data(orig);
        }

        virtual ~SGMatrix()
        {
            unref();
        }

        T* get_column_vector(index_t col) const
        {
            return &matrix[col*num_rows];
        }

        inline const T& operator()(index_t i_row, index_t i_col) const
        {
            return matrix[i_col*num_rows + i_row];
        }

        inline const T& operator[](index_t index) const
        {
            return matrix[index];
        }

        inline T& operator()(index_t i_row, index_t i_col)
        {
            return matrix[i_col*num_rows + i_row];
        }

        inline T& operator[](index_t index)
        {
            return matrix[index];
        }

        inline bool operator==(SGMatrix<T>& other)
        {
            if (num_rows!=other.num_rows || num_cols!=other.num_cols)
                return false;

            if (matrix!=other.matrix)
                return false;

            return true;
        }

        inline bool equals(SGMatrix<T>& other)
        {
            if (num_rows!=other.num_rows || num_cols!=other.num_cols)
                return false;

            for (index_t i=0; i<num_rows*num_cols; ++i)
            {
                if (matrix[i]!=other.matrix[i])
                    return false;
            }

            return true;
        }

        void set_const(T const_elem)
        {
            for (index_t i=0; i<num_rows*num_cols; i++)
                matrix[i]=const_elem ;
        }

        void zero()
        {
            if (matrix && (num_rows*num_cols))
                set_const(0);
        }

        SGMatrix<T> clone()
        {
            return SGMatrix<T>(clone_matrix(matrix, num_rows, num_cols),
                    num_rows, num_cols);
        }

        static T* clone_matrix(const T* matrix, int32_t nrows, int32_t ncols)
        {
            T* result = SG_MALLOC(T, int64_t(nrows)*ncols);
            for (int64_t i=0; i<int64_t(nrows)*ncols; i++)
                result[i]=matrix[i];

            return result;
        }

        static void transpose_matrix(
            T*& matrix, int32_t& num_feat, int32_t& num_vec);

        static void create_diagonal_matrix(T* matrix, T* v,int32_t size)
        {
            for(int32_t i=0;i<size;i++)
            {
                for(int32_t j=0;j<size;j++)
                {
                    if(i==j)
                        matrix[j*size+i]=v[i];
                    else
                        matrix[j*size+i]=0;
                }
            }
        }

        static SGMatrix<T> create_identity_matrix(index_t size, T scale);

        static SGMatrix<float64_t> create_centering_matrix(index_t size);

#ifdef HAVE_LAPACK

        static SGVector<float64_t> compute_eigenvectors(
                SGMatrix<float64_t> matrix);

        static double* compute_eigenvectors(double* matrix, int n, int m);

        void compute_few_eigenvectors(double* matrix_, double*& eigenvalues, double*& eigenvectors,
                                      int n, int il, int iu);
#endif

        static SGMatrix<float64_t> matrix_multiply(
                SGMatrix<float64_t> A, SGMatrix<float64_t> B,
                bool transpose_A=false, bool transpose_B=false,
                float64_t scale=1.0);
#ifdef HAVE_LAPACK

        static void inverse(SGMatrix<float64_t> matrix);

        static float64_t* pinv(
            float64_t* matrix, int32_t rows, int32_t cols,
            float64_t* target=NULL);

#endif

        static inline float64_t trace(
            float64_t* mat, int32_t cols, int32_t rows)
        {
            float64_t trace=0;
            for (int32_t i=0; i<rows; i++)
                trace+=mat[i*cols+i];
            return trace;
        }

        static T* get_row_sum(T* matrix, int32_t m, int32_t n)
        {
            T* rowsums=SG_CALLOC(T, n);

            for (int32_t i=0; i<n; i++)
            {
                for (int32_t j=0; j<m; j++)
                    rowsums[i]+=matrix[j+int64_t(i)*m];
            }
            return rowsums;
        }

        static T* get_column_sum(T* matrix, int32_t m, int32_t n)
        {
            T* colsums=SG_CALLOC(T, m);

            for (int32_t i=0; i<n; i++)
            {
                for (int32_t j=0; j<m; j++)
                    colsums[j]+=matrix[j+int64_t(i)*m];
            }
            return colsums;
        }

        void center()
        {
            center_matrix(matrix, num_rows, num_cols);
        }

        static void center_matrix(T* matrix, int32_t m, int32_t n);

        void remove_column_mean();

        void display_matrix(const char* name="matrix") const;

        static void display_matrix(
            const T* matrix, int32_t rows, int32_t cols,
            const char* name="matrix", const char* prefix="");

        static void display_matrix(
            const SGMatrix<T> matrix, const char* name="matrix",
            const char* prefix="");

        static SGMatrix<T> get_allocated_matrix(index_t num_rows,
                index_t num_cols, SGMatrix<T> pre_allocated=SGMatrix<T>());

    protected:
        virtual void copy_data(const SGReferencedData &orig)
        {
            matrix=((SGMatrix*)(&orig))->matrix;
            num_rows=((SGMatrix*)(&orig))->num_rows;
            num_cols=((SGMatrix*)(&orig))->num_cols;
        }

        virtual void init_data()
        {
            matrix=NULL;
            num_rows=0;
            num_cols=0;
        }

        virtual void free_data()
        {
            SG_FREE(matrix);
            matrix=NULL;
            num_rows=0;
            num_cols=0;
        }

    public:
        T* matrix;
        index_t num_rows;
        index_t num_cols;
};
}
#endif // __SGMATRIX_H__