SimpleMatrix.h

Go to the documentation of this file.

/* 
* Copyright 2014-2018 Friedemann Zenke
*
* This file is part of Auryn, a simulation package for plastic
* spiking neural networks.
* 
* Auryn 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.
* 
* Auryn 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 Auryn.  If not, see <http://www.gnu.org/licenses/>.
*
* If you are using Auryn or parts of it for your work please cite:
* Zenke, F. and Gerstner, W., 2014. Limits to high-speed simulations 
* of spiking neural networks using general-purpose computers. 
* Front Neuroinform 8, 76. doi: 10.3389/fninf.2014.00076
*/

#ifndef SIMPLEMATRIX_H_
#define SIMPLEMATRIX_H_

#include "auryn_global.h"
#include <string.h>

#include <boost/serialization/utility.hpp>
#include <boost/serialization/split_member.hpp>

namespace auryn {

template <typename T>
class SimpleMatrix 
{
private:
        friend class boost::serialization::access;
        template<class Archive>
        void save(Archive & ar, const unsigned int version) const
                {
                        ar & m_rows;
                        ar & n_cols;
                        ar & datasize;
                        ar & n_nonzero;
                        for (NeuronID i = 0 ; i < m_rows ; ++i) {
                                for (NeuronID * r = rowptrs[i] ; r < rowptrs[i+1] ; ++r ) {
                                        ar & i;
                                        ar & *r;
                                        ar & coldata[r-colinds];
                                }
                        }
                }
        template<class Archive>
        void load(Archive & ar, const unsigned int version)
                {
                        AurynLong load_datasize;
                        AurynLong nnz;

                        ar & m_rows;
                        ar & n_cols;
                        ar & load_datasize;
                        ar & nnz;

                        resize_buffer_and_clear(load_datasize);

                        for (AurynLong c = 0 ; c < nnz ; ++c) {
                                NeuronID i,j;
                                T val;
                                ar & i;
                                ar & j;
                                ar & val;
                                push_back(i,j,val);
                        }
                }
        BOOST_SERIALIZATION_SPLIT_MEMBER()

        
        NeuronID m_rows,n_cols;
        NeuronID current_row,current_col;
        AurynLong n_nonzero;
        AurynLong datasize;
protected:
        NeuronID ** rowptrs;
        NeuronID * colinds;
        T * coldata; 

        void init(NeuronID m, NeuronID n, AurynLong size=256);
        void free();

public:
        SimpleMatrix();
        SimpleMatrix(SimpleMatrix * mat);
        SimpleMatrix(NeuronID rows, NeuronID cols);
        SimpleMatrix(NeuronID rows, NeuronID cols, AurynLong size);
        virtual ~SimpleMatrix();

        void clear();

        void resize_buffer(AurynLong newsize);

        void resize_buffer_and_clear(AurynLong size);

        void prune();

        void push_back(const NeuronID i, const NeuronID j, const T value);
        void copy(SimpleMatrix * mat);
        void set_data(AurynLong i, T value);
        void scale_data(AurynLong i, T value);
        T * get_data_ptr(const AurynLong i);
        T get_data(const AurynLong i);
        T * get_data_ptr(const NeuronID * ind_ptr);
        T get_data(const NeuronID * ind_ptr);
        void fill_zeros();
        void fill_na();
        AurynDouble get_fill_level();
        T get(NeuronID i, NeuronID j);
        bool exists(NeuronID i, NeuronID j);
        T * get_ptr(NeuronID i, NeuronID j);
        T * get_ptr(AurynLong data_index);
        T get_value(AurynLong data_index);
        void add_value(AurynLong data_index, T value);
        NeuronID get_colind(AurynLong data_index);
        bool set(NeuronID i, NeuronID j, T value);
        void set_all(T value);
        void set_row(NeuronID i, T value);
        void scale_row(NeuronID i, T value);
        void scale_all(T value);
        void set_col(NeuronID j, T value);
        void scale_col(NeuronID j, T value);
        double sum_col(NeuronID j);
        AurynLong get_datasize();
        AurynLong get_nonzero();
        void print();
        double mean();
        NeuronID * get_ind_begin();
        NeuronID * get_row_begin(NeuronID i);
        AurynLong get_row_begin_index(NeuronID i);
        NeuronID * get_row_end(NeuronID i);
        AurynLong get_row_end_index(NeuronID i);
        NeuronID get_m_rows();
        NeuronID get_n_cols();
        NeuronID ** get_rowptrs();
        T * get_data_begin();
        T * get_data_end();
        T get_value(NeuronID i);
        T get_value(NeuronID * r);
        T * get_value_ptr(NeuronID i);
        T * get_value_ptr(NeuronID * r);
        NeuronID get_data_offset(NeuronID * r);
};

template <typename T>
T * SimpleMatrix<T>::get_data_ptr(AurynLong i)
{
        return coldata+i;
}

template <typename T>
T SimpleMatrix<T>::get_data(AurynLong i)
{
        return coldata[i];
}

template <typename T>
T * SimpleMatrix<T>::get_data_ptr(const NeuronID * ind_ptr) 
{
        size_t ptr_offset = ind_ptr-get_ind_begin();
        return coldata+ptr_offset;
}

template <typename T>
T SimpleMatrix<T>::get_data(const NeuronID * ind_ptr) 
{
        return *(get_data_ptr(ind_ptr));
}

template <typename T>
void SimpleMatrix<T>::set_data(AurynLong i, T value)
{
        if (i<datasize)
                coldata[i] = value;
}

template <typename T>
void SimpleMatrix<T>::scale_data(AurynLong i, T value)
{
        if (i<datasize)
                coldata[i] *= value;
}


template <typename T>
void SimpleMatrix<T>::clear()
{
        current_row = 0;
        current_col = 0;
        n_nonzero = 0;
        rowptrs[0] = colinds;
        rowptrs[1] = colinds;
}

template <typename T>
void SimpleMatrix<T>::init(NeuronID m, NeuronID n, AurynLong size)
{
        m_rows = m;
        n_cols = n;
        datasize = size;
        rowptrs = new NeuronID * [m_rows+1];
        colinds = new NeuronID [datasize];
        coldata = new T [datasize];
        clear();
}

template <typename T>
void SimpleMatrix<T>::resize_buffer(AurynLong newsize)
{
        const AurynLong oldsize = datasize;
        if ( oldsize == newsize ) return;
        const AurynLong copysize = std::min(oldsize,newsize);

        // std::cout << " copy colinds " << std::endl;

        NeuronID * new_colinds = new NeuronID [newsize];
        std::copy(colinds, colinds+copysize, new_colinds);

        // std::cout << " update rowptrs " << std::endl;

        // update rowpointers
        ptrdiff_t offset = new_colinds-colinds;
        for ( NeuronID i = 0 ; i < m_rows+1 ; ++i ) {
                rowptrs[i] += offset;
        }
        
        // std::cout << " delete old colinds " << std::endl;
        // FIXME I get regular crashes in this line:
        //  invalid pointer: 0x0000000001d954c0 ***
        delete [] colinds;

        // std::cout << " replace colinds " << std::endl;
        colinds = new_colinds;

        // std::cout << " create new data array " << std::endl;
        T * new_coldata = new T [newsize];
        std::copy(coldata, coldata+copysize, new_coldata);
        delete [] coldata;
        coldata = new_coldata;

        datasize = newsize;
        // std::cout << " done " << std::endl;
}

template <typename T>
void SimpleMatrix<T>::resize_buffer_and_clear(AurynLong size)
{
        // free();
        // init(m_rows, n_cols, size);
        resize_buffer(size);
        clear();
}

template <typename T>
void SimpleMatrix<T>::prune()
{
        if ( get_datasize() > get_nonzero() )
                resize_buffer(get_nonzero());
}

template <typename T>
SimpleMatrix<T>::SimpleMatrix()
{
        init(1,1);
}

template <typename T>
SimpleMatrix<T>::SimpleMatrix(SimpleMatrix * mat)
{
        init(mat->get_m_rows(),mat->get_n_cols(),mat->get_nonzero());
        copy(mat);
}



template <typename T>
SimpleMatrix<T>::SimpleMatrix(NeuronID rows, NeuronID cols)
{
        init(rows,cols);
}

template <typename T>
SimpleMatrix<T>::SimpleMatrix(NeuronID rows, NeuronID cols, AurynLong datasize)
{
        init(rows,cols,datasize);
}

template <typename T>
void SimpleMatrix<T>::free()
{
        delete [] rowptrs;
        delete [] colinds;
        delete [] coldata;
}

template <typename T>
void SimpleMatrix<T>::copy(SimpleMatrix * mat)
{
        if ( get_m_rows() != mat->get_m_rows() || get_n_cols() != mat->get_n_cols() )
                throw AurynMatrixDimensionalityException();

        clear();

        for ( NeuronID i = 0 ; i < mat->get_m_rows() ; ++i ) {
                for ( NeuronID * r = mat->get_row_begin(i) ; r != mat->get_row_end(i) ; ++r ) {
                        push_back(i,*r,mat->get_value(r));
                }
        }
}

template <typename T>
SimpleMatrix<T>::~SimpleMatrix()
{
        free();
}

template <typename T>
void SimpleMatrix<T>::push_back(NeuronID i, NeuronID j, T value)
{   
        if ( i >= m_rows || j >= n_cols ) throw AurynMatrixDimensionalityException();
        while ( i > current_row )
        {
                rowptrs[current_row+2] = rowptrs[current_row+1];  // save value of last element
                current_col = 0;
                current_row++;
        } 
        current_col = j;
        if (i >= current_row && j >= current_col) {
                if ( n_nonzero >= datasize ) throw AurynMatrixBufferException();
                *(rowptrs[i+1]) = j; // write last j to end of index array
                coldata[rowptrs[i+1]-colinds] = value; // write value to end of data array
                ++rowptrs[i+1]; //increment end by one
                rowptrs[m_rows] = rowptrs[i+1]; // last (m_row+1) marks end of last row
                n_nonzero++;
        } else {
                throw AurynMatrixPushBackException();
        }
}

template <typename T>
AurynLong SimpleMatrix<T>::get_datasize()
{
        return datasize;
}


template <typename T>
AurynLong SimpleMatrix<T>::get_nonzero()
{
        return n_nonzero;
}

template <typename T>
void SimpleMatrix<T>::fill_na()
{
        for ( NeuronID i = current_row ; i < m_rows-1 ; ++i )
        {
                rowptrs[i+2] = rowptrs[i+1];  // save value of last element
        }
        current_row = get_m_rows();
}

template <typename T>
void SimpleMatrix<T>::fill_zeros()
{
        fill_na();
}

template <typename T>
T SimpleMatrix<T>::get(NeuronID i, NeuronID j)
{
        return *get_ptr(i,j);
}

template <typename T>
bool SimpleMatrix<T>::exists(NeuronID i, NeuronID j)
{
        if ( get_ptr(i,j) == NULL )
                return false;
        else 
                return true;
}

template <typename T>
T * SimpleMatrix<T>::get_ptr(NeuronID i, NeuronID j)
{
        // check bounds
        if ( !(i < m_rows && j < n_cols) ) return NULL;

#ifdef DEBUG
        std::cout << "enter search" << std::endl;
#endif // DEBUG
        // perform binary search
        NeuronID * lo = rowptrs[i];
        NeuronID * hi = rowptrs[i+1];

        if ( lo >= hi ) // no elements/targets in this row
                return NULL; 
        
        while ( lo < hi ) {
                NeuronID * c = lo + (hi-lo)/2;
                if ( *c < j ) lo = c+1;
                else hi = c;
#ifdef DEBUG
                std::cout << i << ":" << j << "   " << *lo << ":" << *hi << std::endl;
#endif // DEBUG
        }
        
        if ( *lo == j ) {
                return coldata+(lo-colinds);
        }

        return NULL; 
}

template <typename T>
T * SimpleMatrix<T>::get_ptr(AurynLong data_index)
{
        return &coldata[data_index];
}

template <typename T>
T SimpleMatrix<T>::get_value(AurynLong data_index)
{
        return coldata[data_index];
}

template <typename T>
void SimpleMatrix<T>::add_value(AurynLong data_index, T value)
{
        coldata[data_index] += value;
}

template <typename T>
NeuronID SimpleMatrix<T>::get_colind(AurynLong data_index)
{
        return colinds[data_index];
}

template <typename T>
bool SimpleMatrix<T>::set(NeuronID i, NeuronID j, T value)
{
        T * ptr = get_ptr(i,j);
        if ( ptr != NULL) {
                *ptr = value;
                return true;
        }
        else
                return false;
}

template <typename T>
void SimpleMatrix<T>::scale_row(NeuronID i, T value)
{
        NeuronID * rowbegin = rowptrs[i];
        NeuronID * rowend = rowptrs[i+1]--;

        for (NeuronID * c = rowbegin ; c <= rowend ; ++c) 
        {
                coldata[c-colinds] *= value;
        }
}

template <typename T>
void SimpleMatrix<T>::scale_all(T value)
{
        for ( AurynLong i = 0 ; i < n_nonzero ; ++i ) 
                scale_data( i , value );
}

template <typename T>
void SimpleMatrix<T>::set_row(NeuronID i, T value)
{
        NeuronID * rowbegin = rowptrs[i];
        NeuronID * rowend = rowptrs[i+1]--;

        for (NeuronID * c = rowbegin ; c <= rowend ; ++c) 
        {
                coldata[c-colinds] = value;
        }
}

template <typename T>
void SimpleMatrix<T>::set_all(T value)
{
        for ( AurynLong i = 0 ; i < n_nonzero ; ++i ) 
                set_data( i , value );
}

template <typename T>
void SimpleMatrix<T>::scale_col(NeuronID j, T value)
{
        for ( AurynLong i = 0 ; i < n_nonzero ; ++i ) {
                if ( colinds[i] == j ) scale_data(i,value);
        }
}

template <typename T>
double SimpleMatrix<T>::sum_col(NeuronID j)
{
        double sum = 0;
        for ( AurynLong i = 0 ; i < n_nonzero ; ++i ) {
                if ( colinds[i] == j ) 
                        sum += *(get_data_begin()+i);
        }
        return sum;
}

template <typename T>
void SimpleMatrix<T>::set_col(NeuronID j, T value)
{
        for ( AurynLong i = 0 ; i < n_nonzero ; ++i ) {
                if ( colinds[i] == j ) set_data(i,value);
        }
}

template <typename T>
NeuronID * SimpleMatrix<T>::get_ind_begin()
{
        return rowptrs[0];
}

template <typename T>
NeuronID * SimpleMatrix<T>::get_row_begin(NeuronID i)
{
        return rowptrs[i];
}

template <typename T>
AurynLong SimpleMatrix<T>::get_row_begin_index(NeuronID i)
{
        return rowptrs[i]-rowptrs[0];
}

template <typename T>
NeuronID * SimpleMatrix<T>::get_row_end(NeuronID i)
{
        return rowptrs[i+1];
}

template <typename T>
AurynLong SimpleMatrix<T>::get_row_end_index(NeuronID i)
{
        return rowptrs[i+1]-rowptrs[0];
}


template <typename T>
NeuronID ** SimpleMatrix<T>::get_rowptrs()
{
        return rowptrs;
}

template <typename T>
T * SimpleMatrix<T>::get_data_begin()
{
        return coldata;
}

template <typename T>
T * SimpleMatrix<T>::get_data_end()
{
        return coldata+get_nonzero();
}


template <typename T>
AurynDouble SimpleMatrix<T>::get_fill_level()
{
        return 1.*n_nonzero/datasize;
}

template <typename T>
NeuronID SimpleMatrix<T>::get_m_rows()
{
        return m_rows;
}

template <typename T>
NeuronID SimpleMatrix<T>::get_n_cols()
{
        return n_cols;
}

template <typename T>
void SimpleMatrix<T>::print()
{
        for (NeuronID i = 0 ; i < m_rows ; ++i) {
                for (NeuronID * r = get_row_begin(i) ; r != get_row_end(i) ; ++r ) {
                        std::cout << i << " " << *r << " " << coldata[r-colinds] << "\n";
                }
        }
}

template <typename T>
double SimpleMatrix<T>::mean()
{
        double sum = 0;
        for (NeuronID i = 0 ; i < get_nonzero() ; ++i) {
                sum += coldata[i];
        }
        return sum/get_nonzero();
}

template <typename T>
T SimpleMatrix<T>::get_value(NeuronID i)
{
        return get_data_begin()[i];
}

template <typename T>
T SimpleMatrix<T>::get_value(NeuronID * r)
{
        return get_data_begin()[r-get_ind_begin()];
}

template <typename T>
T * SimpleMatrix<T>::get_value_ptr(NeuronID i)
{
        return &get_data_begin()[i];
}

template <typename T>
T * SimpleMatrix<T>::get_value_ptr(NeuronID * r)
{
        return &get_data_begin()[r-get_ind_begin()];
}

template <typename T>
NeuronID SimpleMatrix<T>::get_data_offset(NeuronID * r)
{
        return r-get_ind_begin();
}

}


#endif /*SIMPLEMATRIX_H_*/