AurynVector.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 AURYNVECTOR_H_
#define AURYNVECTOR_H_

#include <stdio.h>
#include <string.h>
#include <ctime>
#include <assert.h>
#include "auryn_definitions.h"

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

#include <boost/random/mersenne_twister.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/random/normal_distribution.hpp>

namespace auryn {

        template <typename T, typename IndexType = NeuronID > 
        class AurynVector { 
                private: 

                        void * mem;

                        friend class boost::serialization::access;
                        template<class Archive>
                        void serialize(Archive & ar, const unsigned int version)
                        {
                                ar & size;
                                for ( IndexType i = 0 ; i < size ; ++i ) 
                                        ar & data[i];
                        }

                protected:
                        void check_size(IndexType x)
                        {
#ifndef NDEBUG
                                if ( x >= size ) {
                                        throw AurynVectorDimensionalityException();
                                }
#endif 
                        }

                        void check_size(AurynVector * v)
                        {
#ifndef NDEBUG
                                if ( v->size != size ) {
                                        throw AurynVectorDimensionalityException();
                                }
#endif 
                        }

                        void allocate(const NeuronID n) {
#if defined(CODE_USE_SIMD_INSTRUCTIONS_EXPLICITLY) && defined(CODE_ALIGNED_SIMD_INSTRUCTIONS)
                                std::size_t mem_alignment = sizeof(T)*SIMD_NUM_OF_PARALLEL_FLOAT_OPERATIONS;
                                std::size_t mem_size = sizeof(T)*n;
                                mem = malloc(mem_size+mem_alignment-1); // adds padding to allocated memory
                                T * ptr = (T*)mem; 
                                if ( (unsigned long)mem%mem_alignment ) ptr = (T*)(((unsigned long)mem/mem_alignment+1)*mem_alignment);
                                // T * ptr = (T*)boost::alignment::aligned_alloc(sizeof(T)*SIMD_NUM_OF_PARALLEL_FLOAT_OPERATIONS,sizeof(T)*n);
                                if ( mem == NULL ) {
                                        // TODO implement proper exception handling
                                        throw AurynMemoryAlignmentException(); 
                                }
                                assert(((unsigned long)ptr % mem_alignment) == 0);
#else
                                T * ptr = new T[n];
#endif
                                data = ptr;
                                size = n;
                                set_zero();
                        }

                        void freebuf() {
#if defined(CODE_USE_SIMD_INSTRUCTIONS_EXPLICITLY) && defined(CODE_ALIGNED_SIMD_INSTRUCTIONS)
                                free(mem);
                                // boost::alignment::aligned_free(data);
#else
                                delete [] data;
#endif
                        }

            T fast_exp256(T x)
            {
                    x = 1.0 + x / 256.0;
                    x *= x; x *= x; x *= x; x *= x;
                    x *= x; x *= x; x *= x; x *= x;
                    return x;
            }

                protected:

                public:
                        IndexType size;

                        T * data;

                        AurynVector(IndexType n) 
                        {
                                allocate(n);
                        }

                        AurynVector(AurynVector * vec) 
                        {
                                allocate(vec->size);
                                copy(vec);
                        }


                        virtual ~AurynVector() 
                        {
                                freebuf();
                        }

                        virtual void resize(IndexType new_size) 
                        {
                                if ( size != new_size ) {
                                        T * old_data = data;
                                        IndexType old_size = size;
                                        allocate(new_size);
                                        // copy old data
                                        const size_t copy_size = std::min(old_size,new_size) * sizeof(T);
                                        memcpy(data, old_data, copy_size);
                                        free(old_data);
                                }
                        }

                        void copy(AurynVector * v) 
                        {
                                check_size(v);
                                std::copy(v->data, v->data+v->size, data);
                        }

                        T get(IndexType i)
                        {
                                check_size(i);
                                return data[i];
                        }

                        T * ptr(IndexType i = 0)
                        {
                                check_size(i);
                                return data+i;
                        }

                        void set(IndexType i, T value)
                        {
                                check_size(i);
                                data[i] = value;
                        }


                        void set_all(const T v) 
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] = v;
                                }
                        }

                        void set_zero() 
                        {
                                set_all(0.0);
                        }

                        void scale(const T a) 
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] *= a;
                                }
                        }


                        void add(const T c) 
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] += c;
                                }
                        }

                        void add_specific(const IndexType i, const T c) 
                        {
                                check_size(i);
                                data[i] += c;
                        }

                        void mul_specific(const IndexType i, const T c) 
                        {
                                check_size(i);
                                data[i] *= c;
                        }

                        void add(AurynVector * v) 
                        {
                                check_size(v);
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] += v->data[i];
                                }
                        }

                        void sub(const T c) 
                        {
                                add(-c);
                        }

                        void sub(AurynVector * v) 
                        {
                                check_size(v);
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] -= v->data[i];
                                }
                        }

                        void mul(const T a) 
                        {
                                scale(a);
                        }

                        void mul(AurynVector * v) 
                        {
                                check_size(v);
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] *= v->data[i];
                                }
                        }

                        void div(const T a) 
                        {
                                scale(1.0/a);
                        }

                        void div(AurynVector * v) 
                        {
                                check_size(v);
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] /= v->data[i];
                                }
                        }

                        void div(AurynVector * a, AurynVector * b) 
                        {
                                check_size(a);
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] = a->data[i]/b->data[i];
                                }
                        }

                        void saxpy(const T a, AurynVector * x) 
                        {
                                check_size(x);
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] += a * x->data[i];
                                }
                        }

                        void follow(AurynVector<T,IndexType> * v, const T rate) 
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] += rate*(v->data[i]-data[i]);
                                }
                        }

                        void follow_scalar(const T a, const T rate) 
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] += rate*(a-data[i]);
                                }
                        }

                        void elementwise_max(AurynVector * v1, AurynVector * v2) 
                        {
                                check_size(v1);
                                check_size(v2);
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] = std::max(v1->data[i],v2->data[i]);
                                }
                        }

                        void elementwise_max(AurynVector * v1) 
                        {
                                elementwise_max(this,v1);
                        }

                        void pow(const unsigned int n) 
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] = std::pow(data[i],n);
                                }
                        }

                        void fast_exp()
                        {
                                // mul(0.00390625); // ie. 1.0/256.0
                                // add(1.0);
                                // for ( int i = 0 ; i < 8 ; ++i ) {
                                //      sqr();
                                // }
                                // seems as if the naive version is faster
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] = fast_exp256(data[i]);
                                }
                        }

                        void exp()
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] = std::exp(data[i]);
                                }
                        }


                        void sigmoid(AurynVector * x, const T beta, const T thr)
                        {
                                diff(x,thr);
                                mul(-beta);
                                exp();
                                add(1.0);
                                inv();
                        }


                        void sqrt() 
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] = std::sqrt(data[i]);
                                }
                        }

                        void neg()
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] = -data[i];
                                }
                        }

                        void inv()
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] = 1.0/data[i];
                                }
                        }

                        void sum(AurynVector * a, AurynVector * b) 
                        {
                                check_size(a);
                                check_size(b);
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] = a->data[i]+b->data[i];
                                }
                        }

                        void sum(AurynVector * a, const T b) 
                        {
                                check_size(a);
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] = a->data[i]+b;
                                }
                        }

                        void diff(AurynVector * a, AurynVector * b) 
                        {
                                check_size(a);
                                check_size(b);
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        data[i] = a->data[i]-b->data[i];
                                }
                        }

                        void diff(AurynVector * a, const T b) 
                        {
                                sum(a,-b);
                        }

                        void diff(const T a, AurynVector * b) 
                        {
                                check_size(b);
                                sum(b,-a);
                                neg();
                        }


                        void sqr() 
                        {       
                                this->mul(this);
                        }

                        void abs() 
                        {       
                                sqr();
                                sqrt();
                        }

                        void rect()
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        if ( data[i] < 0.0 ) {
                                                 data[i] = 0.0;
                                        } 
                                }
                        }

                        void neg_rect()
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        if ( data[i] > 0.0 ) {
                                                 data[i] = 0.0;
                                        } 
                                }
                        }

                        void clip(T min, T max)
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        if ( data[i] < min ) {
                                                 data[i] = min;
                                        } else 
                                                if ( data[i] > max ) 
                                                         data[i] = max;
                                }
                        }

                        double var()
                        {
                                double sum = 0.0;
                                double sum2 = 0.0;
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        double elem = get(i);
                                        sum  += elem;
                                        sum2 += std::pow(elem,2);
                                }
                                double var  = (sum2-(sum*sum)/size)/(size-1);
                                return var;
                        }


                        double std()
                        {
                                return std::sqrt(var());
                        }

                        double mean()
                        {
                                return element_sum()/size;
                        }

                        double element_sum()
                        {
                                double sum = 0.0;
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        sum += get(i);
                                }
                                return sum;
                        }

                        double l1norm()
                        {
                                double sum = 0.0;
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        double e = get(i);
                                        sum += std::abs(e);
                                }
                                return sum;
                        }


                        double l2norm()
                        {
                                double sum = 0.0;
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        double e = get(i);
                                        sum += e*e;
                                }
                                return std::sqrt(sum);
                        }

                        double max()
                        {
                                double max = -1e64; // TODO use numeric limits for this
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        double el = get(i);
                                        if ( el > max ) max = el;
                                }
                                return max;
                        }

                        double min()
                        {
                                double min = 1e64; // TODO use numeric limits for this
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        double el = get(i);
                                        if ( el < min ) min = el;
                                }
                                return min;
                        }

                        IndexType nonzero()
                        {
                                IndexType sum = 0;
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        if ( get(i) != 0 ) ++sum;
                                }
                                return sum;
                        }

                        void zero_effective_zeros( const T epsilon = 1e-3)
                        {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        if ( std::abs(get(i)) < epsilon ) set(i, 0.0);
                                }
                        }

                        void add_random_normal(AurynState mean=0.0, AurynState sigma=1.0, unsigned int seed=8721)
                        {
                                if ( seed == 0 )
                                        seed = static_cast<unsigned int>(std::time(0));
                                boost::mt19937 randgen(seed); 
                                boost::normal_distribution<> dist((double)mean, (double)sigma);
                                boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > die(randgen, dist);
                                AurynState rv;
                                for ( IndexType i = 0 ; i<size ; ++i ) {
                                        rv = die();
                                        data[i] += rv;
                                }
                        }

                        void set_random_normal(AurynState mean=0.0, AurynState sigma=1.0, unsigned int seed=8721)
                        {
                                set_all(0.0);
                                add_random_normal(mean, sigma);
                        }

                        void set_random(unsigned int seed = 0) 
                        {
                                set_random_normal(0.0,1.0,seed);
                        }

                        bool any() {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        if ( get(i) ) return true;
                                }
                                return false;
                        }

                        bool any(T eps) {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        if ( std::abs(get(i))>eps ) return true;
                                }
                                return false;
                        }

                        void print() {
                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        std::cout << get(i) << " ";
                                }
                                std::cout << std::endl;
                        }

                        void write_to_file(std::string filename) {

                                std::ofstream outfile;
                                outfile.open(filename.c_str(),std::ios::out);
                                if (!outfile) {
                                        std::stringstream oss;
                                        throw AurynOpenFileException();
                                }

                                outfile << std::setprecision(7);

                                for ( IndexType i = 0 ; i < size ; ++i ) {
                                        outfile << get(i) << "\n";
                                }

                                outfile.close();
                        }
        };




        class AurynVectorFloat : public AurynVector<float,NeuronID> 
        {
                private:
                        typedef AurynVector<float,NeuronID> super;

                public:
                        AurynVectorFloat(NeuronID n);

                        ~AurynVectorFloat() 
                        {
                        };


                        virtual void resize(NeuronID new_size);
                        void scale(const float a);
                        void saxpy(const float a, AurynVectorFloat * x);
                        void clip(const float min, const float max);
                        void add(const float c);
                        void add(AurynVectorFloat * v);
                        void sum(AurynVectorFloat * a, AurynVectorFloat * b);
                        void sum(AurynVectorFloat * a, const float b);
                        void mul(const float a) { scale(a); };
                        void mul(AurynVectorFloat * v);
                        void diff(AurynVectorFloat * a, AurynVectorFloat * b);
                        void diff(AurynVectorFloat * a, const float b);
                        void diff(const float a, AurynVectorFloat * b );
                        void follow(AurynVectorFloat * v, const float rate);

                        // TODO add pow function with intrinsics _mm_pow_ps

        };

}


#endif /*AURYNVECTOR_H_*/