c_index_map_imp.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 20, 2013
//     Author: Eugen Wintersberger <eugen.wintersberger@desy.de>
//

#pragma once

#include <iostream>
#include <sstream>
#include <array>
#include <numeric>
#include <algorithm>
#include <functional>

#include "../../error/exceptions.hpp"
#include "../../error/exception_utils.hpp"

namespace pni{
namespace core{

    class c_index_map_imp
    {
        private:
            template<typename IITERT,
                     typename SITERT> 
            static size_t offset(IITERT &&index_start,
                                 IITERT &&index_stop, 
                                 SITERT &&shape_start)
            {
                //initialize the offset and the stride variable
                size_t offset = *index_start++,stride=1;

                //loop over all indices - remember that the iterators are
                //reverse iterators
                while(index_start!=index_stop)
                {
                    //compute the actuall stride 
                    stride *= *shape_start++;
                    //compute the offset contribution
                    offset += stride*(*index_start++);
                }

                return offset;
            }

            //------------------------------------------------------------------
            template<
                     typename SELITER,
                     typename SITER,
                     typename IITER
                    > 
            static size_t offset(SELITER &&sel_shape_start,
                                 SELITER &&sel_shape_end,
                                 SELITER &&sel_offset, 
                                 SELITER &&sel_stride,
                                 SITER   &&shape_start,
                                 IITER   &&sel_index)
            {
                size_t index = *sel_offset++;

                if(*sel_shape_start++ != 1) index += (*sel_index++)*(*sel_stride);
                ++sel_stride;

                //initialize the offset and the stride variable
                size_t offset = index ,dim_stride=1;

                //loop over all indices - remember that the iterators are
                //reverse iterators
                while(sel_shape_start!=sel_shape_end)
                {
                    //compute the index from the selection
                    index = *sel_offset++;
                    if(*sel_shape_start++ != 1) 
                        index += (*sel_index++)*(*sel_stride);

                    ++sel_stride; //increment the selection stride in any case

                    //compute the actuall stride 
                    dim_stride *= *shape_start++;
                    //compute the offset contribution
                    offset += dim_stride*index;
                }

                return offset;
            }


            //------------------------------------------------------------------
            template<
                     typename IITERT,
                     typename SITERT
                    >
            static void index(SITERT &&shape_start,
                              SITERT &&shape_stop,
                              IITERT &&index_start,
                              size_t offset)
            {
                size_t t;
                size_t stride = std::accumulate(++shape_start,shape_stop,1,
                                                std::multiplies<size_t>());
                t = offset%stride;
                *(index_start++) = (offset-t)/stride;
                offset = t;
                while(shape_start != shape_stop)
                {
                    //increment here the shape_start iterator - we already start
                    //with start+1 with the stride computation
                    stride /= *shape_start++;
                    t = offset%stride;
                    *(index_start++) = (offset-t)/stride;
                    offset = t;
                }
            }
        public:

            //-----------------------------------------------------------------
            template<
                     typename CSHAPE,
                     typename CINDEX
                    >
            static size_t offset(const CSHAPE &shape,const CINDEX &index)
            {
                //use here reverse iterators as required by the c-ordering
                return offset(index.rbegin(),index.rend(),shape.rbegin());
            }

            //-----------------------------------------------------------------
            template<
                     typename SELTYPE,
                     typename CSHAPE,
                     typename SINDEX
                    >
            static size_t offset(const SELTYPE &sel,const CSHAPE &shape,
                                 const SINDEX &index)
            {
                return offset(sel.full_shape().rbegin(), //original selection shape
                              sel.full_shape().rend(),   //
                              sel.offset().rbegin(), //selection offset
                              sel.stride().rbegin(), //selection stride
                              shape.rbegin(), //original shape begin
                              index.rbegin()); //selection index start
            }

            //-----------------------------------------------------------------
            template<
                     typename CINDEX,
                     typename CSHAPE
                    >
            static void index(const CSHAPE &shape,CINDEX &idx,size_t offset)
            {
                //for index computation we can use the forward iterators
                index(shape.begin(),shape.end(),idx.begin(),offset); 
            }

            //-----------------------------------------------------------------
            template<
                     typename CTYPE,
                     typename SELTYPE
                    >
            static bool is_contiguous(const CTYPE &shape,const SELTYPE &sel)
            {
                auto oshape = sel.template full_shape<CTYPE>();
                
                auto oiter = oshape.begin();
                auto siter = shape.begin();

                for(;siter != shape.end(); ++oiter,++siter)
                {
                    if((*oiter!=1)&&(*oiter!=*siter)) return false;
                }

                return true;

            }

    };
//end of namespace
}
}