mdarray.hpp

//
// (c) Copyright 2013 DESY, Eugen Wintersberger <eugen.wintersberger@desy.de>
//
// This file is part of libpnicore.
//
// libpnicore 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.
//
// libpnicore 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 libpnicore.  If not, see <http://www.gnu.org/licenses/>.
//
// ============================================================================
//
// Created on: Oct 28, 2013
//     Author: Eugen Wintersberger <eugen.wintersberger@desy.de>
//

#pragma once

//#include<memory>
#include<iostream>
#include<sstream>
#include<stdexcept>
#include<utility>
#include<complex>
#include<cstdarg>
#include<cstdio>
#include<memory>


#include "../error/exception_utils.hpp"
#include "../types/types.hpp"
#include "../utilities.hpp"
#include "slice.hpp"
#include "array_view.hpp"
#include "array_factory.hpp"
#include "array_view_utils.hpp"
#include "index_map/index_maps.hpp"
#include "../algorithms.hpp"


namespace pni {
namespace core {
    
    
    template<
             typename STORAGE,
             typename IMAP=dynamic_cindex_map,
             typename IPA =inplace_arithmetics
            > 
    class mdarray
    {
        public:
            //================public types=====================================
            typedef STORAGE storage_type;
            typedef typename storage_type::value_type value_type;  
            typedef IMAP map_type;
            typedef mdarray<storage_type,map_type,IPA> array_type;
            typedef typename storage_type::iterator iterator;
            typedef typename storage_type::const_iterator const_iterator; 
            typedef typename storage_type::reverse_iterator reverse_iterator;
            typedef typename storage_type::const_reverse_iterator
                const_reverse_iterator;
            typedef IPA inplace_arithmetic;
            typedef array_view<array_type> view_type;
            typedef array_view<const array_type> const_view_type;
            typedef size_t size_type;

            
            //==================public members=================================
            static const type_id_t type_id; 
        private:
            IMAP _imap;  
            STORAGE _data;  
        public:

            //=================constructors and destructor=====================
            mdarray():_imap(),_data() {}

            //-----------------------------------------------------------------
            explicit mdarray(const map_type &map,const storage_type &s):
                _imap(map),
                _data(s)
            {}

            //-----------------------------------------------------------------
            explicit mdarray(map_type &&map,storage_type &&s):
                _imap(std::move(map)),
                _data(std::move(s))
            {}

            //-----------------------------------------------------------------
            template<typename ATYPE>
            explicit mdarray(const array_view<ATYPE> &view):
                _imap(map_utils<map_type>::create(view.template shape<shape_t>())),
                _data(container_utils<storage_type>::create(view.size()))
            {
                std::copy(view.begin(),view.end(),_data.begin());
            }


            //-----------------------------------------------------------------
            template<typename ...MDARGS>
            explicit mdarray(const mdarray<MDARGS...> &array):
                _imap(map_utils<map_type>::create(array.template shape<shape_t>())),
                _data(container_utils<storage_type>::create(array.size()))
            {
                //copy data
                for(size_t i=0;i<array.size();++i) (*this)[i] = array[i];
                //std::copy(array.begin(),array.end(),this->begin());
            }

            //====================static methods to create arrays==============
            template<typename ...ARGS>
            static array_type create(ARGS... arguments)
            {
                typedef array_factory<array_type> factory_type;

                return factory_type::create(std::forward<ARGS>(arguments)...);
            }

            //-----------------------------------------------------------------
            template<
                     typename T1,
                     typename T2
                    >
            static array_type create(const std::initializer_list<T1> &l1,
                                     const std::initializer_list<T2> &l2)
            {
                typedef array_factory<array_type> factory_type;

                return factory_type::create(l1,l2);
            }



            //====================assignment operations========================
            template<typename ...MDARGS>
            array_type &operator=(const mdarray<MDARGS...> &array)
            {
                if((void*)this == (void*)&array) return *this;
    
                size_t s = array.size();
                for(size_t i=0;i<s;++i) (*this)[i] = array[i];

                return *this;
            }

            //-----------------------------------------------------------------
            array_type &operator=(const std::initializer_list<value_type> &l)
            {
                if(l.size() != this->size())
                    throw size_mismatch_error(EXCEPTION_RECORD,
                            "Size of array and initializer list do not match");

                std::copy(l.begin(),l.end(),this->begin());
                return *this;
            }
            

            //================public methods===================================
            const map_type &map() const { return this->_imap; }

            //----------------------------------------------------------------
            template<typename CTYPE> CTYPE shape() const
            {
                auto c = container_utils<CTYPE>::create(_imap.rank());
                std::copy(_imap.begin(),_imap.end(),c.begin());
                return c;
            }

            //-----------------------------------------------------------------
            size_t size() const 
            { 
                return _data.size(); 
            }

            //-----------------------------------------------------------------
            size_t rank() const 
            { 
                return _imap.rank(); 
            }

            //=============operators and methods to access array data==========
            value_type& operator[](size_t i) 
            { 
#ifdef DEBUG
                return at(i);
#else
                return _data[i]; 
#endif
            }

            //-----------------------------------------------------------------
            value_type operator[](size_t i) const 
            { 
#ifdef DEBUG
                return at(i);
#else
                return _data[i]; 
#endif
            }

            //-----------------------------------------------------------------
            value_type &at(size_t i) 
            { 
                try
                {
                    return _data.at(i); 
                }
                catch(std::out_of_range &error)
                {
                    std::stringstream ss;
                    ss<<"Index "<<i<<" is out of range ("<<size()<<")!";
                    throw index_error(EXCEPTION_RECORD,ss.str());
                }
            
            } 

            //-----------------------------------------------------------------
            value_type at(size_t i) const 
            { 
                try
                {
                    return _data.at(i); 
                }
                catch(std::out_of_range &error)
                {
                    std::stringstream ss;
                    ss<<"Index "<<i<<" is out of range ("<<size()<<")!";
                    throw index_error(EXCEPTION_RECORD,ss.str());
                }
            
            } 

            //-----------------------------------------------------------------
            void insert(size_t i,const value_type &value)
            {
                at(i)=value;
            }
            //-----------------------------------------------------------------
            template<
                     typename CTYPE,
                     typename = typename enable_element_cont<CTYPE>::type
                    >
            value_type &operator()(const CTYPE &index)
            {
               return _data[_imap.offset(index)]; 
            }
            
            //-----------------------------------------------------------------
            template<
                     typename CTYPE,
                     typename = typename enable_element_cont<CTYPE>::type
                    >
            value_type operator()(const CTYPE &index) const
            {
                return _data[_imap.offset(index)];
            }
            
            //-----------------------------------------------------------------
            template<
                     typename CTYPE,
                     typename = typename enable_view_cont<CTYPE>::type
                    >
            array_view<const array_type> 
            operator()(const CTYPE &slices) const
            {
                typedef array_view<const array_type> view_type;

                return view_type(*this,array_selection::create(slices));
                
            }

            //-----------------------------------------------------------------
            template<
                     typename CTYPE,
                     typename = typename enable_view_cont<CTYPE>::type
                    >
            array_view<array_type> 
            operator()(const CTYPE &slices)
            {
                typedef array_view<array_type> view_type;

                return view_type(*this, array_selection::create(slices));
                
            }

            //-----------------------------------------------------------------
            template<
                     typename... ITYPES,
                     typename = typename enable_valid_index<ITYPES...>::type
                    >
            typename view_type_trait<array_type,ITYPES...>::type 
            operator()(ITYPES... indexes)
            {
                typedef view_provider<array_type,is_view_index<ITYPES...>::value>
                    provider_type;
                return provider_type::get_reference(*this,_imap,indexes...);
            }


            //-----------------------------------------------------------------
            template<
                     typename... ITYPES,
                     typename = typename enable_valid_index<ITYPES...>::type
                    >
            typename view_type_trait<const array_type,ITYPES...>::const_type
            operator()(ITYPES ...indexes) const
            {
                typedef view_provider<array_type,is_view_index<ITYPES...>::value>
                    provider_type;
                return provider_type::get_value(*this,_imap,indexes...);
            }


            
            //-----------------------------------------------------------------
            const value_type *data() const
            {
                return _data.data();
            }

            //-----------------------------------------------------------------
            value_type *data()
            {
                return _data.data();
            }

            //-----------------------------------------------------------------
            value_type &front() { return _data.front(); }

            //-----------------------------------------------------------------
            value_type front() const { return _data.front(); }

            //-----------------------------------------------------------------
            value_type &back() { return _data.back(); }

            //-----------------------------------------------------------------
            value_type back() const { return _data.back(); }

            //-----------------------------------------------------------------
            iterator begin() { return _data.begin(); }

            //-----------------------------------------------------------------
            iterator end() { return _data.end(); }

            //-----------------------------------------------------------------
            const_iterator begin() const { return _data.begin(); }

            //-----------------------------------------------------------------
            const_iterator end() const { return _data.end(); }

            //-----------------------------------------------------------------
            reverse_iterator rbegin() { return _data.rbegin(); }

            //-----------------------------------------------------------------
            const_reverse_iterator rbegin() const { return _data.rbegin(); }

            //-----------------------------------------------------------------
            reverse_iterator rend() { return this->_data.rend(); }

            //-----------------------------------------------------------------
            const_reverse_iterator rend() const { return _data.rend(); }

            //==========implementation of unary arithmetic operators===========

            array_type &operator+=(value_type s) 
            { 
                IPA::add(*this,s); 
                return *this;
            }

            //-----------------------------------------------------------------

            template<typename ATYPE> 
            array_type &operator+=(const ATYPE &v) 
            { 
                IPA::add(*this,v); 
                return *this;
            }

            //-----------------------------------------------------------------

            array_type &operator-=(value_type s) 
            { 
                IPA::sub(*this,s); 
                return *this;
            }

            //-----------------------------------------------------------------

            template<typename ATYPE> 
            array_type &operator-=(const ATYPE &v) 
            { 
                IPA::sub(*this,v); 
                return *this; 
            }

            //-----------------------------------------------------------------

            array_type &operator*=(value_type s) 
            { 
                IPA::mult(*this,s); 
                return *this; 
            }

            //-----------------------------------------------------------------

            template<typename ATYPE>
            array_type &operator*=(const ATYPE &v) 
            { 
                IPA::mult(*this,v); 
                return *this;
            }

            //-----------------------------------------------------------------

            array_type &operator/=(value_type s) 
            { 
                IPA::div(*this,s); 
                return *this;
            }

            //-----------------------------------------------------------------

            template<typename ATYPE>
            array_type &operator/=(const ATYPE &v) 
            { 
                IPA::div(*this,v); 
                return *this;
            }


    };

    //-------------------------------------------------------------------------
    template<
             typename STORAGE,
             typename IMAP,
             typename IPA
            >
    struct container_trait<mdarray<STORAGE,IMAP,IPA>>
    {
        static const bool is_random_access = true;
        static const bool is_iterable = true;
        static const bool is_contiguous =
            container_trait<STORAGE>::is_contiguous;
        static const bool is_multidim = true;
    };

    //set data for static member attribute
    template<typename STORAGE,typename IMAP,typename IPA>
    const type_id_t mdarray<STORAGE,IMAP,IPA>::type_id = 
    type_id_map<typename STORAGE::value_type>::type_id;
    //=====================non-member operators================================

    template<
             typename STORAGE,
             typename IMAP,
             typename IPA
            >
    std::ostream &operator<<(std::ostream &o,const mdarray<STORAGE,IMAP,IPA> &a)
    {
        for(auto iter = a.begin();iter!=a.end();++iter)
            o<<*iter<<" ";

        return o;
    }

    //-------------------------------------------------------------------------
    template<
             typename STORAGE,
             typename IMAP,
             typename IPA
            >
    std::istream &operator>>(std::istream &is,mdarray<STORAGE,IMAP,IPA> &a)
    {
        for(auto iter=a.begin();iter!=a.end();++iter)
            is>>*iter;

        return is;
    }
   
    //-------------------------------------------------------------------------
    template<
             typename STORAGE,
             typename IMAP,
             typename IPA
            >
    bool operator==(const mdarray<STORAGE,IMAP,IPA> &b1, 
                    const mdarray<STORAGE,IMAP,IPA> &b2) 
    {
        return std::equal(b1.begin(),b1.end(),b2.begin());
    }

    //-------------------------------------------------------------------------
    template<
             typename STORAGE,
             typename IMAP,
             typename IPA
            >
    bool operator!=(const mdarray<STORAGE,IMAP,IPA> &b1, 
                    const mdarray<STORAGE,IMAP,IPA> &b2) 
    {
        return !(b1==b2);
    }

    
}
}