text:
/***************************************************************************
*   Copyright 2012 - 2013 Advanced Micro Devices, Inc.
*
*   Licensed under the Apache License, Version 2.0 (the "License");
*   you may not use this file except in compliance with the License.
*   You may obtain a copy of the License at
*
*       http://www.apache.org/licenses/LICENSE-2.0
*
*   Unless required by applicable law or agreed to in writing, software
*   distributed under the License is distributed on an "AS IS" BASIS,
*   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
*   See the License for the specific language governing permissions and
*   limitations under the License.

***************************************************************************/

#if !defined( AMP_SORT_INL )
#define AMP_SORT_INL
#pragma once

#include <algorithm>
#include <type_traits>

//#include <boost/bind.hpp>
//#include <boost/thread/once.hpp>
//#include <boost/shared_array.hpp>

#include "bolt/amp/bolt.h"
#include "bolt/amp/scan.h"
#include "bolt/amp/functional.h"
#include "bolt/amp/device_vector.h"
#include <amp.h>

#ifdef ENABLE_TBB
#include "tbb/parallel_sort.h"
#include "tbb/task_scheduler_init.h"
#endif

#define BOLT_UINT_MAX 0xFFFFFFFFU
#define BOLT_UINT_MIN 0x0U
#define BOLT_INT_MAX 0x7FFFFFFFU
#define BOLT_INT_MIN 0x80000000U

#define BITONIC_SORT_WGSIZE 64
/* \brief - SORT_CPU_THRESHOLD should be atleast 2 times the BITONIC_SORT_WGSIZE*/
#define SORT_CPU_THRESHOLD 128

namespace bolt {
namespace amp {
template<typename RandomAccessIterator>
void sort(RandomAccessIterator first,
          RandomAccessIterator last)
{
    typedef std::iterator_traits< RandomAccessIterator >::value_type T;

    detail::sort_detect_random_access( bolt::amp::control::getDefault( ),
                                       first, last, less< T >( ),
                                       std::iterator_traits< RandomAccessIterator >::iterator_category( ) );
    return;
}

template<typename RandomAccessIterator, typename StrictWeakOrdering>
void sort(RandomAccessIterator first,
            RandomAccessIterator last,
            StrictWeakOrdering comp)
{
    detail::sort_detect_random_access( bolt::amp::control::getDefault( ),
                                       first, last, comp,
                                       std::iterator_traits< RandomAccessIterator >::iterator_category( ) );
    return;
}

template<typename RandomAccessIterator>
void sort(bolt::amp::control &ctl,
            RandomAccessIterator first,
            RandomAccessIterator last)
{
    typedef std::iterator_traits< RandomAccessIterator >::value_type T;
    detail::sort_detect_random_access(ctl,
                                        first, last, less< T >( ),
                                        std::iterator_traits< RandomAccessIterator >::iterator_category( ) );
    return;
}

template<typename RandomAccessIterator, typename StrictWeakOrdering>
void sort(bolt::amp::control &ctl,
            RandomAccessIterator first,
            RandomAccessIterator last,
            StrictWeakOrdering comp)
{
    detail::sort_detect_random_access(ctl,
                                        first, last, comp,
                                        std::iterator_traits< RandomAccessIterator >::iterator_category( ) );
    return;
}
}// End of amp namespace
};// End of bolt namespace

namespace bolt {
namespace amp {
namespace detail {

// Wrapper that uses default control class, iterator interface
template<typename RandomAccessIterator, typename StrictWeakOrdering>
void sort_detect_random_access( bolt::amp::control &ctl,
                                const RandomAccessIterator& first, const RandomAccessIterator& last,
                                const StrictWeakOrdering& comp, std::input_iterator_tag )
{
    //  \TODO:  It should be possible to support non-random_access_iterator_tag iterators, if we copied the data
    //  to a temporary buffer.  Should we?
    static_assert( false, "Bolt only supports random access iterator types" );
};

template<typename RandomAccessIterator, typename StrictWeakOrdering>
void sort_detect_random_access( bolt::amp::control &ctl,
                                const RandomAccessIterator& first, const RandomAccessIterator& last,
                                const StrictWeakOrdering& comp, std::random_access_iterator_tag )
{
    return sort_pick_iterator(ctl, first, last, comp,
                              std::iterator_traits< RandomAccessIterator >::iterator_category( ) );
};

//Device Vector specialization
template<typename DVRandomAccessIterator, typename StrictWeakOrdering>
void sort_pick_iterator( bolt::amp::control &ctl,
                         const DVRandomAccessIterator& first, const DVRandomAccessIterator& last,
                         const StrictWeakOrdering& comp, bolt::amp::device_vector_tag )
{
    // User defined Data types are not supported with device_vector. Hence we have a static assert here.
    // The code here should be in compliant with the routine following this routine.
    typedef typename std::iterator_traits<DVRandomAccessIterator>::value_type T;
    unsigned int szElements = static_cast< unsigned int >( std::distance( first, last ) );
    if(szElements > (1<<31))
  {
        throw std::exception( "AMP device_vector shall support only upto 2 power 31 elements" );
    return;
  }
    if (szElements == 0 )
        return;
    const bolt::amp::control::e_RunMode runMode = ctl.getForceRunMode();  // could be dynamic choice some day.
    if ((runMode == bolt::amp::control::SerialCpu) || (szElements < SORT_CPU_THRESHOLD)) {
        std::vector<T> localBuffer(szElements);
        //Copy the device_vector buffer to a CPU buffer
        for(unsigned int index=0; index<szElements; index++)
            localBuffer[index] = first.getBuffer()[index];
        //Compute sort using std::sort
        std::sort(localBuffer.begin(), localBuffer.end(), comp);
        //Copy the CPU buffer back to device_vector
        for(unsigned int index=0; index<szElements; index++)
            first.getBuffer()[index] = localBuffer[index];
        return;
    } else if (runMode == bolt::amp::control::MultiCoreCpu) {
#ifdef ENABLE_TBB
        std::vector<T> localBuffer(szElements);
        //Copy the device_vector buffer to a CPU buffer
        for(unsigned int index=0; index<szElements; index++)
            localBuffer[index] = first.getBuffer()[index];
        tbb::task_scheduler_init initialize(tbb::task_scheduler_init::automatic);
        tbb::parallel_sort(localBuffer.begin(), localBuffer.end(), comp);

        std::sort(localBuffer.begin(), localBuffer.end(), comp);

        //Copy the CPU buffer back to device_vector
        for(unsigned int index=0; index<szElements; index++)
            first.getBuffer()[index] = localBuffer[index];
        return;

#else
 //       std::cout << "The MultiCoreCpu version of sort is not enabled. " << std ::endl;
        throw std::exception( "The MultiCoreCpu version of sort is not enabled to be built." );
        return;
#endif
    } else {
        sort_enqueue(ctl,first,last,comp);
    }
    return;
}

#if 0
template<typename DVRandomAccessIterator, typename StrictWeakOrdering>
void sort_pick_iterator( control &ctl,
                         const DVRandomAccessIterator& first, const DVRandomAccessIterator& last,
                         const StrictWeakOrdering& comp, bolt::amp::fancy_iterator_tag )
{
    static_assert( false, "It is not possible to sort fancy iterators. They are not mutable" );
}
#endif

//Non Device Vector specialization.
//This implementation creates a cl::Buffer and passes the cl buffer to the sort specialization whichtakes
//the cl buffer as a parameter. In the future, Each input buffer should be mapped to the device_vector
//and the specialization specific to device_vector should be called.
template<typename RandomAccessIterator, typename StrictWeakOrdering>
void sort_pick_iterator( bolt::amp::control &ctl,
                         const RandomAccessIterator& first, const RandomAccessIterator& last,
                         const StrictWeakOrdering& comp, std::random_access_iterator_tag )
{
    typedef typename std::iterator_traits<RandomAccessIterator>::value_type T;
    unsigned int szElements = static_cast< unsigned int >( std::distance( first, last ) );
    if(szElements > (1<<31))
        throw std::exception( "AMP device vectors shall support only upto 2 power 31 elements" );
    if (szElements == 0)
        return;

    const bolt::amp::control::e_RunMode runMode = ctl.getForceRunMode();  // could be dynamic choice some day.
    if ((runMode == bolt::amp::control::SerialCpu) || (szElements < SORT_CPU_THRESHOLD)) {
        std::sort(first, last, comp);
        return;
    } else if (runMode == bolt::amp::control::MultiCoreCpu) {
#ifdef ENABLE_TBB
        tbb::task_scheduler_init initialize(tbb::task_scheduler_init::automatic);
        tbb::parallel_sort(first,last, comp);
#else
//        std::cout << "The MultiCoreCpu version of sort is not enabled. " << std ::endl;
        throw std::exception( "The MultiCoreCpu version of sort is not enabled to be built." );
        return;
#endif
    } else {
        device_vector< T, concurrency::array_view > dvInputOutput( first, last, false, ctl );
        //Now call the actual amp algorithm
        sort_enqueue(ctl,dvInputOutput.begin(),dvInputOutput.end(),comp);
        //Map the buffer back to the host
        dvInputOutput.data( );
        return;
    }
}

/****** sort_enqueue specailization for unsigned int data types. ******
 * THE FOLLOWING CODE IMPLEMENTS THE RADIX SORT ALGORITHM FOR sort()
 *********************************************************************/

/*AMP Kernels for unsigned integer sorting*/
template <typename T, int N, typename Container>
void AMP_RadixSortHistogramAscendingKernel(bolt::amp::control &ctl,
                                           Container &unsortedData, /*this can be either array_view or array*/
                                           Container &buckets,
                                           Container &histScanBuckets,
                                           unsigned int shiftCount,
                                           unsigned int szElements,
                                           const unsigned int groupSize)
{
    const int RADICES = 1 << N;
    concurrency::extent< 1 > globalSizeK0( szElements/RADICES );

    concurrency::tiled_extent< RADICES > tileK0 = globalSizeK0.tile< RADICES >();
    concurrency::accelerator_view av = ctl.getAccelerator().default_view;
    printf("globalSizeK0 = %d   tileK0.tile_dim0 = %d    tileK0[0]=%d\n",globalSizeK0, tileK0.tile_dim0, tileK0[0]);
    concurrency::parallel_for_each( av, tileK0,
    [
        unsortedData,
        buckets,
        histScanBuckets,
        shiftCount,
        tileK0
    ]
    ( concurrency::tiled_index< RADICES > t_idx ) restrict(amp)
    {
        const int RADIX_T     = N;
        const int RADICES_T   = (1 << RADIX_T);
        const int NUM_OF_ELEMENTS_PER_WORK_ITEM_T = RADICES_T;
        const int MASK_T      = (1<<RADIX_T)  - 1;

        int localId     = t_idx.local[ 0 ];
        int globalId    = t_idx.global[ 0 ];
        int groupId     = t_idx.tile[ 0 ];
        int groupSize   = tileK0.tile_dim0;
        //tileK0[ 0 ] gives the global size
        //tileK0.tile_dim0 gives the number of threads in a tile
        //\TODO - find a API to get the total number of numOfGroups
        int numOfGroups = tileK0[ 0 ]/tileK0.tile_dim0;//get_num_groups(0);
        unsigned int bucketPos   = groupId * RADICES_T * groupSize;
        for(int i = 0; i < RADICES_T; ++i)
        {
            buckets[bucketPos + localId * RADICES_T + i] = 0;
        }
        //barrier(CLK_GLOBAL_MEM_FENCE)
        t_idx.barrier.wait_with_all_memory_fence();

        /* Calculate thread-histograms */
        for(int i = 0; i < NUM_OF_ELEMENTS_PER_WORK_ITEM_T; ++i)
        {
            unsigned int value = unsortedData[globalId * NUM_OF_ELEMENTS_PER_WORK_ITEM_T + i];
            value = (value >> shiftCount) & MASK_T;
            buckets[bucketPos + localId * RADICES_T + value]++;
        }

        //barrier(CLK_GLOBAL_MEM_FENCE)
        t_idx.barrier.wait_with_all_memory_fence();

        //Start First step to scan
        int sum =0;
        for(int i = 0; i < groupSize; i++)
        {
            sum = sum + buckets[bucketPos + localId + groupSize*i];
        }
        histScanBuckets[localId*numOfGroups + groupId + 1] = sum;

    } );// end of concurrency::parallel_for_each
}

template <typename T, int N, typename Container>
void AMP_RadixSortHistogramDescendingKernel(bolt::amp::control &ctl,
                                           Container &unsortedData, /*this can be either array_view or array*/
                                           Container &buckets,
                                           Container &histScanBuckets,
                                           unsigned int shiftCount,
                                           unsigned int szElements,
                                           const unsigned int groupSize)
{
    const int RADICES = 1 << N;
    concurrency::extent< 1 > globalSizeK0( szElements/RADICES );
    concurrency::tiled_extent< RADICES > tileK0 = globalSizeK0.tile< RADICES >();
    concurrency::accelerator_view av = ctl.getAccelerator().default_view;
    concurrency::parallel_for_each( av, tileK0,
    [
        unsortedData,
        buckets,
        histScanBuckets,
        shiftCount,
        tileK0
    ]
    ( concurrency::tiled_index< RADICES > t_idx ) restrict(amp)
    {
        const int RADIX_T     = N;
        const int RADICES_T   = (1 << RADIX_T);
        const int NUM_OF_ELEMENTS_PER_WORK_ITEM_T = RADICES_T;
        const int MASK_T      = (1<<RADIX_T)  -1;
        int localId     = t_idx.local[ 0 ];
        int globalId    = t_idx.global[ 0 ];
        int groupId     = t_idx.tile[ 0 ];
        int groupSize   = tileK0.tile_dim0;
        //tileK0[ 0 ] gives the global size
        //tileK0.tile_dim0 gives the number of threads in a tile
        //\TODO - find a API to get the total number of numOfGroups
        int numOfGroups = tileK0[ 0 ]/tileK0.tile_dim0;
        /*size_t localId     = get_local_id(0);
        size_t globalId    = get_global_id(0);
        size_t groupId     = get_group_id(0);
        size_t groupSize   = get_local_size(0);
        size_t numOfGroups = get_num_groups(0);*/
        unsigned int bucketPos   = groupId * RADICES_T * groupSize;

        for(int i = 0; i < RADICES_T; ++i)
        {
            buckets[bucketPos + localId * RADICES_T + i] = 0;
        }
        //barrier(CLK_GLOBAL_MEM_FENCE)
        t_idx.barrier.wait_with_all_memory_fence();

        /* Calculate thread-histograms */
        for(int i = 0; i < NUM_OF_ELEMENTS_PER_WORK_ITEM_T; ++i)
        {
            unsigned int value = unsortedData[globalId * NUM_OF_ELEMENTS_PER_WORK_ITEM_T + i];
            value = (value >> shiftCount) & MASK_T;
            buckets[bucketPos + localId * RADICES_T + (RADICES_T - value -1)]++;
        }

        //barrier(CLK_GLOBAL_MEM_FENCE)
        t_idx.barrier.wait_with_all_memory_fence();
        //Start First step to scan
        int sum = 0;
        for(int i = 0; i < groupSize; i++)
            sum = sum + buckets[bucketPos + localId + groupSize*i];
        histScanBuckets[localId*numOfGroups + groupId + 1] = sum;
    } );// end of concurrency::parallel_for_each
}

template <typename T, int N, typename Container>
void AMP_scanLocalTemplate(bolt::amp::control &ctl,
                           Container &buckets,
                           Container &histScanBuckets,
                           unsigned int szElements)
{
    const int RADICES = 1<<N;
    concurrency::extent< 1 > globalSizeK0( szElements/RADICES );
    concurrency::tiled_extent< RADICES > tileK0 = globalSizeK0.tile< RADICES >();
    concurrency::accelerator_view av = ctl.getAccelerator().default_view;
    concurrency::parallel_for_each( av, tileK0,
    [
        buckets,
        histScanBuckets,
        tileK0
    ]
    ( concurrency::tiled_index< RADICES > t_idx ) restrict(amp)
    {
        const int RADIX_T     = N;
        const int RADICES_T   = (1 << RADIX_T);
        //create tiled_static for localScanArray
        tile_static T localScanArray[2*RADICES_T];
        /*size_t localId     = get_local_id(0);
        size_t numOfGroups = get_num_groups(0);
        size_t groupId     = get_group_id(0);
        size_t groupSize   = get_local_size(0);*/
        int localId     = t_idx.local[ 0 ];
        int groupId     = t_idx.tile[ 0 ];
        int groupSize   = tileK0.tile_dim0;
        //\TODO - find a API to get the total number of numOfGroups
        int numOfGroups = tileK0[ 0 ]/tileK0.tile_dim0;

        localScanArray[localId] = histScanBuckets[localId*numOfGroups + groupId];
        localScanArray[RADICES_T+localId] = 0;

        //barrier(CLK_LOCAL_MEM_FENCE);
        t_idx.barrier.wait();
        for(int i = 0; i < RADICES_T; ++i)
        {
            unsigned int bucketPos = groupId * RADICES_T * groupSize + i * RADICES_T + localId;
            unsigned int temp = buckets[bucketPos];
            buckets[bucketPos] = localScanArray[RADICES_T+localId] + localScanArray[localId];
            localScanArray[RADICES_T+localId] += temp;
        }
    });
}

template <typename T, int N, typename Container>
void AMP_permuteAscendingRadixNTemplate(bolt::amp::control &ctl,
                                    Container &unsortedData,
                                    Container &scanedBuckets,
                                    unsigned int shiftCount,
                                    Container &sortedData,
                                    unsigned int szElements)
{
    const int RADICES = 1 << N;
    concurrency::extent< 1 > globalSizeK0( szElements/RADICES );
    concurrency::tiled_extent< RADICES > tileK0 = globalSizeK0.tile< RADICES >();
    concurrency::accelerator_view av = ctl.getAccelerator().default_view;
    concurrency::parallel_for_each( av, tileK0,
    [
        unsortedData,
        scanedBuckets,
        shiftCount,
        sortedData,
        tileK0
    ]
    ( concurrency::tiled_index< RADICES > t_idx ) restrict(amp)
    {
        const int RADIX_T     = N;
        const int RADICES_T   = (1 << RADIX_T);
        //const int NUM_OF_ELEMENTS_PER_WORK_ITEM_T = RADICES_T;
        const int MASK_T      = (1<<RADIX_T)  -1;
        /*size_t groupId   = get_group_id(0);
        size_t localId   = get_local_id(0);
        size_t globalId  = get_global_id(0);
        size_t groupSize = get_local_size(0);*/
        int localId     = t_idx.local[ 0 ];
        int globalId    = t_idx.global[ 0 ];
        int groupId     = t_idx.tile[ 0 ];
        int groupSize   = tileK0.tile_dim0;

        unsigned int bucketPos   = groupId * RADICES_T * groupSize;

        /* Premute elements to appropriate location */
        for(int i = 0; i < RADICES_T; ++i)
        {
            unsigned int value = unsortedData[globalId * RADICES_T + i];
            value = (value >> shiftCount) & MASK_T;
            unsigned int index = scanedBuckets[bucketPos+localId * RADICES_T + value];
            sortedData[index] = unsortedData[globalId * RADICES_T + i];
            scanedBuckets[bucketPos+localId * RADICES_T + value] = index + 1;
            //barrier(CLK_LOCAL_MEM_FENCE);
            t_idx.barrier.wait();
        }
    } );
}

template <typename T, int N, typename Container>
void AMP_permuteDescendingRadixNTemplate(bolt::amp::control &ctl,
                                    Container &unsortedData,
                                    Container &scanedBuckets,
                                    unsigned int shiftCount,
                                    Container &sortedData,
                                    unsigned int szElements)
{
    const int RADICES = 1 << N;
    concurrency::extent< 1 > globalSizeK0( szElements/RADICES );
    concurrency::tiled_extent< RADICES > tileK0 = globalSizeK0.tile< RADICES >();
    concurrency::accelerator_view av = ctl.getAccelerator().default_view;
    concurrency::parallel_for_each( av, tileK0,
    [
        unsortedData,
        scanedBuckets,
        shiftCount,
        sortedData,
        tileK0
    ]
    ( concurrency::tiled_index< RADICES > t_idx ) restrict(amp)
    {
        const int RADIX_T     = N;
        const int RADICES_T   = (1 << RADIX_T);
        const int MASK_T      = (1<<RADIX_T)  -1;
        int localId     = t_idx.local[ 0 ];
        int globalId    = t_idx.global[ 0 ];
        int groupId     = t_idx.tile[ 0 ];
        int groupSize   = tileK0.tile_dim0;
        /*size_t groupId   = get_group_id(0);
        size_t localId   = get_local_id(0);
        size_t globalId  = get_global_id(0);
        size_t groupSize = get_local_size(0);*/
        unsigned int bucketPos   = groupId * RADICES_T * groupSize;

        /* Premute elements to appropriate location */
        for(int i = 0; i < RADICES_T; ++i)
        {
            unsigned int value = unsortedData[globalId * RADICES_T + i];
            value = (value >> shiftCount) & MASK_T;
            unsigned int index = scanedBuckets[bucketPos+localId * RADICES_T + (RADICES_T -value-1)];
            sortedData[index] = unsortedData[globalId * RADICES_T + i];
            scanedBuckets[bucketPos+localId * RADICES_T + (RADICES_T -value-1)] = index + 1;
            //barrier(CLK_LOCAL_MEM_FENCE);
            t_idx.barrier.wait();
        }
    } );
}

/*AMP Kernels for Signed integer sorting*/
/*!
* \brief This template function is a permutation Algorithm for one work item.
*        This is called by the the sort_enqueue for input buffers with siigned numbers routine for other than int's and unsigned int's
*
* \details
*         The idea behind sorting signed numbers is that when we sort the MSB which contains the sign bit. We should
*         sort in the signed bit in the descending order. For example in a radix 4 sort of signed numbers if we are
*         doing ascending sort. Then the bits 28:31 should be sorted in an descending order, after doing the following
*         transformation of the input bits.
*         Shift bits 31:28 ->  3:0
*               value = (value >> shiftCount);
*         Retain the sign bit in the bit location 3.
*               signBit = value & (1<<(RADIX_T-1));
*         XOR bits 2:1  with 111
*               value = ( ( ( value & MASK_T ) ^ MASK_T ) & MASK_T )
*         Or the sign bit with the xor'ed bit. This forms your index.
*               value = value | signBitl
*
*       As an Example for Ascending order sort of 4 bit signed numbers.
*       ________________________________
*       -Flip the three bits other than the sign bit
*       -sort the resultant array considering them as unsigned numbers
*       -map the corresponding flipped numbers with the corresponding values.
*       1st Example
*       1111     1000      0000     0111
*       0111     0000      0101     0010
*       1000 ->  1111 ->   0111 ->  0000
*       1010     1101      1000     1111
*       0000     0111      1101     1010
*       0010     0101      1111     1000
*
*       2nd Example
*       1111     1000     1111      1000
*       0111     0000     1101      1010
*       1101     1010     1010      1101
*       1000     1111     1000      1111
*       1010 ->  1101 ->  0111 ->   0000
*       0000     0111     0101      0010
*       0010     0101     0010      0101
*       0101     0010     0001      0110
*       0110     0001     0000      0111
*
*/
/*Descending*/
template <typename T, int N, typename Container>
void AMP_RadixSortHistogramSignedAscendingKernel(bolt::amp::control &ctl,
                                           Container &unsortedData, /*this can be either array_view or array*/
                                           Container &buckets,
                                           Container &histScanBuckets,
                                           unsigned int shiftCount,
                                           unsigned int szElements,
                                           const unsigned int groupSize)
{
    const int RADICES = 1 << N;
    concurrency::extent< 1 > globalSizeK0( szElements/RADICES );

    concurrency::tiled_extent< RADICES > tileK0 = globalSizeK0.tile< RADICES >();
    concurrency::accelerator_view av = ctl.getAccelerator().default_view;
    printf("globalSizeK0 = %d   tileK0.tile_dim0 = %d    tileK0[0]=%d\n",globalSizeK0, tileK0.tile_dim0, tileK0[0]);
    concurrency::parallel_for_each( av, tileK0,
    [
        unsortedData,
        buckets,
        histScanBuckets,
        shiftCount,
        tileK0
    ]
    ( concurrency::tiled_index< RADICES > t_idx ) restrict(amp)
    {
        const int RADIX_T     = N;
        const int RADICES_T   = (1 << RADIX_T);
        const int NUM_OF_ELEMENTS_PER_WORK_ITEM_T = RADICES_T;
        const int MASK_T      = ( 1 << ( RADIX_T - 1 ) ) - 1;

        int localId     = t_idx.local[ 0 ];
        int globalId    = t_idx.global[ 0 ];
        int groupId     = t_idx.tile[ 0 ];
        int groupSize   = tileK0.tile_dim0;
        //tileK0[ 0 ] gives the global size
        //tileK0.tile_dim0 gives the number of threads in a tile
        //\TODO - find a API to get the total number of numOfGroups
        int numOfGroups = tileK0[ 0 ]/tileK0.tile_dim0;//get_num_groups(0);
        unsigned int bucketPos   = groupId * RADICES_T * groupSize;
        for(int i = 0; i < RADICES_T; ++i)
        {
            buckets[bucketPos + localId * RADICES_T + i] = 0;
        }
        //barrier(CLK_GLOBAL_MEM_FENCE)
        t_idx.barrier.wait_with_all_memory_fence();

        /* Calculate thread-histograms */
        for(int i = 0; i < NUM_OF_ELEMENTS_PER_WORK_ITEM_T; ++i)
        {
            unsigned int value = unsortedData[globalId * NUM_OF_ELEMENTS_PER_WORK_ITEM_T + i];
            value = (value >> shiftCount);
            unsigned int signBit = value & (1<<(RADIX_T-1));
            value = ( ( ( value & MASK_T ) ^ MASK_T ) & MASK_T ) | signBit;
            buckets[bucketPos + localId * RADICES_T + value]++;
        }

        //barrier(CLK_GLOBAL_MEM_FENCE)
        t_idx.barrier.wait_with_all_memory_fence();

        //Start First step to scan
        int sum =0;
        for(int i = 0; i < groupSize; i++)
        {
            sum = sum + buckets[bucketPos + localId + groupSize*i];
        }
        histScanBuckets[localId*numOfGroups + groupId + 1] = sum;

    } );// end of concurrency::parallel_for_each
}

template <typename T, int N, typename Container>
void AMP_RadixSortHistogramSignedDescendingKernel(bolt::amp::control &ctl,
                                           Container &unsortedData, /*this can be either array_view or array*/
                                           Container &buckets,
                                           Container &histScanBuckets,
                                           unsigned int shiftCount,
                                           unsigned int szElements,
                                           const unsigned int groupSize)
{
    const int RADICES = 1 << N;
    concurrency::extent< 1 > globalSizeK0( szElements/RADICES );
    concurrency::tiled_extent< RADICES > tileK0 = globalSizeK0.tile< RADICES >();
    concurrency::accelerator_view av = ctl.getAccelerator().default_view;
    concurrency::parallel_for_each( av, tileK0,
    [
        unsortedData,
        buckets,
        histScanBuckets,
        shiftCount,
        tileK0
    ]
    ( concurrency::tiled_index< RADICES > t_idx ) restrict(amp)
    {
        const int RADIX_T     = N;
        const int RADICES_T   = (1 << RADIX_T);
        const int NUM_OF_ELEMENTS_PER_WORK_ITEM_T = RADICES_T;
        const int MASK_T      = ( 1 << ( RADIX_T - 1 ) ) - 1;
        int localId     = t_idx.local[ 0 ];
        int globalId    = t_idx.global[ 0 ];
        int groupId     = t_idx.tile[ 0 ];
        int groupSize   = tileK0.tile_dim0;
        //tileK0[ 0 ] gives the global size
        //tileK0.tile_dim0 gives the number of threads in a tile
        //\TODO - find a API to get the total number of numOfGroups
        int numOfGroups = tileK0[ 0 ]/tileK0.tile_dim0;
        /*size_t localId     = get_local_id(0);
        size_t globalId    = get_global_id(0);
        size_t groupId     = get_group_id(0);
        size_t groupSize   = get_local_size(0);
        size_t numOfGroups = get_num_groups(0);*/
        unsigned int bucketPos   = groupId * RADICES_T * groupSize;

        for(int i = 0; i < RADICES_T; ++i)
        {
            buckets[bucketPos + localId * RADICES_T + i] = 0;
        }
        //barrier(CLK_GLOBAL_MEM_FENCE)
        t_idx.barrier.wait_with_all_memory_fence();

        /* Calculate thread-histograms */
        for(int i = 0; i < NUM_OF_ELEMENTS_PER_WORK_ITEM_T; ++i)
        {
            unsigned int value = unsortedData[globalId * NUM_OF_ELEMENTS_PER_WORK_ITEM_T + i];
            value = (value >> shiftCount);
            unsigned int signBit = value & (1<<(RADIX_T-1));
            value = ( ( ( ( value & MASK_T ) ^ MASK_T ) & MASK_T ) | signBit );
            buckets[bucketPos + localId * RADICES_T + (RADICES_T - value -1)]++;
        }

        //barrier(CLK_GLOBAL_MEM_FENCE)
        t_idx.barrier.wait_with_all_memory_fence();
        //Start First step to scan
        int sum = 0;
        for(int i = 0; i < groupSize; i++)
            sum = sum + buckets[bucketPos + localId + groupSize*i];
        histScanBuckets[localId*numOfGroups + groupId + 1] = sum;
    } );// end of concurrency::parallel_for_each
}

template <typename T, int N, typename Container>
void AMP_permuteSignedAscendingRadixNTemplate(bolt::amp::control &ctl,
                                    Container &unsortedData,
                                    Container &scanedBuckets,
                                    unsigned int shiftCount,
                                    Container &sortedData,
                                    unsigned int szElements)
{
    const int RADICES = 1 << N;
    concurrency::extent< 1 > globalSizeK0( szElements/RADICES );
    concurrency::tiled_extent< RADICES > tileK0 = globalSizeK0.tile< RADICES >();
    concurrency::accelerator_view av = ctl.getAccelerator().default_view;
    concurrency::parallel_for_each( av, tileK0,
    [
        unsortedData,
        scanedBuckets,
        shiftCount,
        sortedData,
        tileK0
    ]
    ( concurrency::tiled_index< RADICES > t_idx ) restrict(amp)
    {
        const int RADIX_T     = N;
        const int RADICES_T   = (1 << RADIX_T);
        //const int NUM_OF_ELEMENTS_PER_WORK_ITEM_T = RADICES_T;
        const int MASK_T      = ( 1 << ( RADIX_T - 1 ) )  -1;
        /*size_t groupId   = get_group_id(0);
        size_t localId   = get_local_id(0);
        size_t globalId  = get_global_id(0);
        size_t groupSize = get_local_size(0);*/
        int localId     = t_idx.local[ 0 ];
        int globalId    = t_idx.global[ 0 ];
        int groupId     = t_idx.tile[ 0 ];
        int groupSize   = tileK0.tile_dim0;

        unsigned int bucketPos   = groupId * RADICES_T * groupSize;

        /* Premute elements to appropriate location */
        for(int i = 0; i < RADICES_T; ++i)
        {
            unsigned int value = unsortedData[globalId * RADICES_T + i];
            unsigned int resultValue = value;
            value = (value >> shiftCount);
            unsigned int signBit = value & (1<<(RADIX_T-1));
            value = ( ( ( ( value & MASK_T ) ^ MASK_T ) & MASK_T ) | signBit );

            unsigned int index = scanedBuckets[bucketPos+localId * RADICES_T + value];
            sortedData[index] = unsortedData[globalId * RADICES_T + i];
            scanedBuckets[bucketPos+localId * RADICES_T + value] = index + 1;
            //barrier(CLK_LOCAL_MEM_FENCE);
            t_idx.barrier.wait();
        }
    } );
}

template <typename T, int N, typename Container>
void AMP_permuteSignedDescendingRadixNTemplate(bolt::amp::control &ctl,
                                    Container &unsortedData,
                                    Container &scanedBuckets,
                                    unsigned int shiftCount,
                                    Container &sortedData,
                                    unsigned int szElements)
{
    const int RADICES = 1 << N;
    concurrency::extent< 1 > globalSizeK0( szElements/RADICES );
    concurrency::tiled_extent< RADICES > tileK0 = globalSizeK0.tile< RADICES >();
    concurrency::accelerator_view av = ctl.getAccelerator().default_view;
    concurrency::parallel_for_each( av, tileK0,
    [
        unsortedData,
        scanedBuckets,
        shiftCount,
        sortedData,
        tileK0
    ]
    ( concurrency::tiled_index< RADICES > t_idx ) restrict(amp)
    {
        const int RADIX_T     = N;
        const int RADICES_T   = (1 << RADIX_T);
        const int MASK_T      = ( 1 << ( RADIX_T - 1 ) )  -1;
        int localId     = t_idx.local[ 0 ];
        int globalId    = t_idx.global[ 0 ];
        int groupId     = t_idx.tile[ 0 ];
        int groupSize   = tileK0.tile_dim0;
        /*size_t groupId   = get_group_id(0);
        size_t localId   = get_local_id(0);
        size_t globalId  = get_global_id(0);
        size_t groupSize = get_local_size(0);*/
        unsigned int bucketPos   = groupId * RADICES_T * groupSize;

        /* Premute elements to appropriate location */
        for(int i = 0; i < RADICES_T; ++i)
        {
            unsigned int value = unsortedData[globalId * RADICES_T + i];
            //value = (value >> shiftCount) & MASK_T;
            value = (value >> shiftCount);
            unsigned int signBit = value & (1<<(RADIX_T-1));
            value = ( ( ( ( value & MASK_T ) ^ MASK_T ) & MASK_T ) | signBit );
            unsigned int index = scanedBuckets[bucketPos+localId * RADICES_T + (RADICES_T -value-1)];
            sortedData[index] = unsortedData[globalId * RADICES_T + i];
            scanedBuckets[bucketPos+localId * RADICES_T + (RADICES_T -value-1)] = index + 1;
            //barrier(CLK_LOCAL_MEM_FENCE);
            t_idx.barrier.wait();
        }
    } );
}

#define BOLT_SORT_INL_DEBUG 0
template<typename DVRandomAccessIterator, typename StrictWeakOrdering>
typename std::enable_if< std::is_same< typename std::iterator_traits<DVRandomAccessIterator >::value_type, unsigned int >::value >::type
sort_enqueue(bolt::amp::control &ctl,
             DVRandomAccessIterator &first, DVRandomAccessIterator &last,
             StrictWeakOrdering comp)
{
    typedef typename std::iterator_traits< DVRandomAccessIterator >::value_type T;
    const int RADIX = 8;
    const int RADICES = (1 << RADIX);
    unsigned int orig_szElements = static_cast<unsigned int>(std::distance(first, last));
    unsigned int szElements = orig_szElements;
    bool  newBuffer = false;
    concurrency::array_view<T> *pLocalArrayView = NULL;
    concurrency::array<T>      *pLocalArray = NULL;
    unsigned int groupSize = RADICES;
    unsigned int mulFactor = groupSize * RADICES;
    concurrency::extent<1> ext( static_cast< int >( orig_szElements ) );

    if(orig_szElements%mulFactor != 0)
    {
        szElements  = ((orig_szElements + mulFactor) /mulFactor) * mulFactor;
        concurrency::extent<1> modified_ext( static_cast< int >( szElements ) );
        pLocalArray     = new concurrency::array<T>( modified_ext );
        pLocalArrayView = new concurrency::array_view<T>(pLocalArray->view_as(modified_ext));
        concurrency::array_view<T> dest = pLocalArrayView->section( ext );
        first.getBuffer( ).copy_to( dest );
        dest.synchronize( );
        newBuffer = true;
    }
    else
    {
        pLocalArrayView = new concurrency::array_view<T>( first.getBuffer( ) );
    }

    unsigned int numGroups = szElements / mulFactor;
    concurrency::extent<1> modified_ext( static_cast< int >( szElements ) );

    device_vector< T, concurrency::array > dvSwapInputData(static_cast<size_t>(szElements), 0);
    device_vector< T, concurrency::array > dvHistogramBins(static_cast<size_t>(numGroups* groupSize * RADICES), 0);
    device_vector< T, concurrency::array > dvHistogramScanBuffer(static_cast<size_t>(numGroups* RADICES + 10), 0 );

    auto& clInputData = *pLocalArrayView;
    auto& clSwapData = dvSwapInputData.begin( ).getBuffer( );
    auto& clHistData = dvHistogramBins.begin( ).getBuffer( );
    auto& clHistScanData = dvHistogramScanBuffer.begin( ).getBuffer( );
    int swap = 0;
    if(comp(2,3))
    {
        /*If the buffer is a local buffer which is more than the usual buffer size then */
        if(newBuffer == true)
        {
            concurrency::index<1> origin(orig_szElements);
            concurrency::array_view<T> dest = pLocalArrayView->section( origin, modified_ext - ext);
            //arrayview_type m_devMemoryAV( *m_devMemory );
            Concurrency::parallel_for_each( dest.extent, [dest]
                (Concurrency::index<1> idx) restrict(amp)
                {
                    dest[idx] = BOLT_UINT_MAX;
                }
            );
        }
        /*Ascending Sort*/
        for(int bits = 0; bits < (sizeof(T) * 8)/*Bits per Byte*/; bits += RADIX)
        {
            if (swap == 0)
                AMP_RadixSortHistogramAscendingKernel<T, RADIX>(ctl,
                                                      clInputData, /*this can be either array_view or array*/
                                                      clHistData,
                                                      clHistScanData,
                                                      bits,
                                                      szElements,
                                                      groupSize); // \TODO - i don't need gropu size
            else
                AMP_RadixSortHistogramAscendingKernel<T, RADIX>(ctl,
                                                      clSwapData, /*this can be either array_view or array*/
                                                      clHistData,
                                                      clHistScanData,
                                                      bits,
                                                      szElements,
                                                      groupSize); // \TODO - i don't need gropu size
#if BOLT_SORT_INL_DEBUG

printf("\n\n\n\n\nBITS = %d\nAfter Histogram", bits);
for (unsigned int ng=0; ng<numGroups; ng++)
{ printf ("\nGroup-Block =%d",ng);
    for(unsigned int gS=0;gS<groupSize; gS++)
    { printf ("\nGroup =%d\n",gS);
        for(int i=0; i<RADICES;i++)
        {
            int index = ng * groupSize * RADICES + gS * RADICES + i;
            int value = clHistData[ index ];
            printf("%2x %2x, ",index, value);
        }
    }
}
int temp = 0;
printf("\n Printing Histogram scan SUM\n");
for(int i=0; i<RADICES;i++)
{
    printf ("\nRadix = %d\n",i);
    for (unsigned int ng=0; ng<numGroups; ng++)
    {
        printf ("%4x, ",clHistScanData[i*numGroups + ng]);
    }
}
#endif
            detail::scan_enqueue( ctl, dvHistogramScanBuffer.begin( ), dvHistogramScanBuffer.end( ),dvHistogramScanBuffer.begin( ), 0, plus< T >( ) );
#if BOLT_SORT_INL_DEBUG
printf("\nprinting scan_enqueue SUM\n");
        for(int i=0; i<RADICES;i++)
        {
            printf ("\nRadix = %d\n",i);
            for (unsigned int ng=0; ng<numGroups; ng++)
            {
                printf ("%4x, ",clHistScanData[i*numGroups + ng]);
            }
        }
#endif
            AMP_scanLocalTemplate<T, RADIX>(ctl,
                                        clHistData,
                                        clHistScanData,
                                        szElements);
#if BOLT_SORT_INL_DEBUG
        printf("\n\nAfter Scan bits = %d", bits);
        for (unsigned int ng=0; ng<numGroups; ng++)
        { printf ("\nGroup-Block =%d",ng);
            for(unsigned int gS=0;gS<groupSize; gS++)
            { printf ("\nGroup =%d\n",gS);
                for(int i=0; i<RADICES;i++)
                {
                    int index = ng * groupSize * RADICES + gS * RADICES + i;
                    int value = clHistData[ index ];
                    printf("%4x %4x, ",index, value);
                }
            }
        }
#endif
            if (swap == 0)
                AMP_permuteAscendingRadixNTemplate<T, RADIX>( ctl,
                                                        clInputData,
                                                        clHistData,
                                                        bits,
                                                        clSwapData,
                                                        szElements);
            else
                AMP_permuteAscendingRadixNTemplate<T, RADIX>( ctl,
                                                        clSwapData,
                                                        clHistData,
                                                        bits,
                                                        clInputData,
                                                        szElements);
#if BOLT_SORT_INL_DEBUG
if (swap == 0)
{
        printf("\n Printing swap data\n");
        for(unsigned int i=0; i<szElements;i+= RADICES)
        {
            for(int j =0;j< RADICES;j++)
                printf("%8x %8x, ",i+j,clSwapData[i+j]);
            printf("\n");
        }
}
else
{
        printf("\n Printing swap data\n");
        for(unsigned int i=0; i<szElements;i+= RADICES)
        {
            for(int j =0;j< RADICES;j++)
                printf("%8x %8x, ",i+j,clInputData[i+j]);
            printf("\n");
        }
}
#endif
            if(swap==0)
                swap = 1;
            else
                swap = 0;
        }
    }
    else
    {
        if(newBuffer == true)
        {
            concurrency::index<1> origin(orig_szElements);
            concurrency::array_view<T> dest = pLocalArrayView->section( origin, modified_ext - ext);
            //arrayview_type m_devMemoryAV( *m_devMemory );
            Concurrency::parallel_for_each( dest.extent, [dest]
                (Concurrency::index<1> idx) restrict(amp)
                {
                    dest[idx] = BOLT_UINT_MIN;
                }
            );
        }
        /*Ascending Sort*/
        for(int bits = 0; bits < (sizeof(T) * 8)/*Bits per Byte*/; bits += RADIX)
        {
            if (swap == 0)
                AMP_RadixSortHistogramDescendingKernel<T, RADIX>(ctl,
                                                      clInputData, /*this can be either array_view or array*/
                                                      clHistData,
                                                      clHistScanData,
                                                      bits,
                                                      szElements,
                                                      groupSize); // \TODO - i don't need gropu size
            else
                AMP_RadixSortHistogramDescendingKernel<T, RADIX>(ctl,
                                                      clSwapData, /*this can be either array_view or array*/
                                                      clHistData,
                                                      clHistScanData,
                                                      bits,
                                                      szElements,
                                                      groupSize); // \TODO - i don't need gropu size

            detail::scan_enqueue( ctl, dvHistogramScanBuffer.begin( ), dvHistogramScanBuffer.end( ),dvHistogramScanBuffer.begin( ), 0, plus< T >( ) );

            AMP_scanLocalTemplate<T, RADIX>(ctl,
                                        clHistData,
                                        clHistScanData,
                                        szElements);

            if (swap == 0)
                AMP_permuteDescendingRadixNTemplate<T, RADIX>( ctl,
                                                        clInputData,
                                                        clHistData,
                                                        bits,
                                                        clSwapData,
                                                        szElements);
            else
                AMP_permuteDescendingRadixNTemplate<T, RADIX>( ctl,
                                                        clSwapData,
                                                        clHistData,
                                                        bits,
                                                        clInputData,
                                                        szElements);
            /*For swapping the buffers*/
            swap = swap? 0: 1;
        }//End of For loop
    }
    if(newBuffer == true)
    {
        //std::cout << "New buffer was allocated So copying back the buffer\n";
        //dest = clInputData.section( ext );
        clInputData.section( ext ).copy_to( first.getBuffer( ) );
        first.getBuffer( ).synchronize( );
        delete pLocalArray;
    }
    delete pLocalArrayView;
    return;
}

/*
 *
 */
template<typename DVRandomAccessIterator, typename StrictWeakOrdering>
typename std::enable_if< std::is_same< typename std::iterator_traits<DVRandomAccessIterator >::value_type,          int >::value >::type
sort_enqueue(bolt::amp::control &ctl,
             DVRandomAccessIterator &first, DVRandomAccessIterator &last,
             StrictWeakOrdering comp)
{
    typedef typename std::iterator_traits< DVRandomAccessIterator >::value_type T;
    const int RADIX = 4;
    const int RADICES = (1 << RADIX);
    unsigned int orig_szElements = static_cast<unsigned int>(std::distance(first, last));
    unsigned int szElements = orig_szElements;
    bool  newBuffer = false;
    concurrency::array_view<T> *pLocalArrayView = NULL;
    concurrency::array<T>      *pLocalArray = NULL;
    unsigned int groupSize = RADICES;
    unsigned int mulFactor = groupSize * RADICES;
    concurrency::extent<1> ext( static_cast< int >( orig_szElements ) );

    if(orig_szElements%mulFactor != 0)
    {
        szElements  = ((orig_szElements + mulFactor) /mulFactor) * mulFactor;
        concurrency::extent<1> modified_ext( static_cast< int >( szElements ) );
        pLocalArray     = new concurrency::array<T>( modified_ext );
        pLocalArrayView = new concurrency::array_view<T>(pLocalArray->view_as( modified_ext ) );
        concurrency::array_view<T> dest = pLocalArrayView->section( ext );
        first.getBuffer( ).copy_to( dest );
        dest.synchronize( );
        newBuffer = true;
    }
    else
    {
        pLocalArrayView = new concurrency::array_view<T>( first.getBuffer( ) );
    }

    unsigned int numGroups = szElements / mulFactor;
    concurrency::extent<1> modified_ext( static_cast< int >( szElements ) );

    device_vector< T, concurrency::array > dvSwapInputData(static_cast<size_t>(szElements), 0);
    device_vector< T, concurrency::array > dvHistogramBins(static_cast<size_t>(numGroups* groupSize * RADICES), 0);
    device_vector< T, concurrency::array > dvHistogramScanBuffer(static_cast<size_t>(numGroups* RADICES + 10), 0 );

    auto& clInputData = *pLocalArrayView;
    auto& clSwapData = dvSwapInputData.begin( ).getBuffer( );
    auto& clHistData = dvHistogramBins.begin( ).getBuffer( );
    auto& clHistScanData = dvHistogramScanBuffer.begin( ).getBuffer( );
    int swap = 0;
    if(comp(2,3))
    {
        /*If the buffer is a local buffer which is more than the usual buffer size then */
        if(newBuffer == true)
        {
            concurrency::index<1> origin(orig_szElements);
            concurrency::array_view<T> dest = pLocalArrayView->section( origin, modified_ext - ext);
            //arrayview_type m_devMemoryAV( *m_devMemory );
            Concurrency::parallel_for_each( dest.extent, [dest]
                (Concurrency::index<1> idx) restrict(amp)
                {
                    dest[idx] = BOLT_INT_MAX;
                }
            );
        }
        /*Ascending Sort*/
        int bits = 0;
        for(bits = 0; bits < ((sizeof(T) * 8) - RADIX); bits += RADIX)
        {
            if (swap == 0)
                AMP_RadixSortHistogramAscendingKernel<T, RADIX>(ctl,
                                                      clInputData, /*this can be either array_view or array*/
                                                      clHistData,
                                                      clHistScanData,
                                                      bits,
                                                      szElements,
                                                      groupSize); // \TODO - i don't need group size
            else
                AMP_RadixSortHistogramAscendingKernel<T, RADIX>(ctl,
                                                      clSwapData, /*this can be either array_view or array*/
                                                      clHistData,
                                                      clHistScanData,
                                                      bits,
                                                      szElements,
                                                      groupSize); // \TODO - i don't need group size
#if BOLT_SORT_INL_DEBUG

printf("\n\n\n\n\nBITS = %d\nAfter Histogram", bits);
for (unsigned int ng=0; ng<numGroups; ng++)
{ printf ("\nGroup-Block =%d",ng);
    for(unsigned int gS=0;gS<groupSize; gS++)
    { printf ("\nGroup =%d\n",gS);
        for(int i=0; i<RADICES;i++)
        {
            int index = ng * groupSize * RADICES + gS * RADICES + i;
            int value = clHistData[ index ];
            printf("%2x %2x, ",index, value);
        }
    }
}
int temp = 0;
printf("\n Printing Histogram scan SUM\n");
for(int i=0; i<RADICES;i++)
{
    printf ("\nRadix = %d\n",i);
    for (unsigned int ng=0; ng<numGroups; ng++)
    {
        printf ("%4x, ",clHistScanData[i*numGroups + ng]);
    }
}
#endif
            detail::scan_enqueue( ctl, dvHistogramScanBuffer.begin( ), dvHistogramScanBuffer.end( ),dvHistogramScanBuffer.begin( ), 0, plus< T >( ) );
#if BOLT_SORT_INL_DEBUG
printf("\nprinting scan_enqueue SUM\n");
        for(int i=0; i<RADICES;i++)
        {
            printf ("\nRadix = %d\n",i);
            for (unsigned int ng=0; ng<numGroups; ng++)
            {
                printf ("%4x, ",clHistScanData[i*numGroups + ng]);
            }
        }
#endif
            AMP_scanLocalTemplate<T, RADIX>(ctl,
                                        clHistData,
                                        clHistScanData,
                                        szElements);
#if BOLT_SORT_INL_DEBUG
        printf("\n\nAfter Scan bits = %d", bits);
        for (unsigned int ng=0; ng<numGroups; ng++)
        { printf ("\nGroup-Block =%d",ng);
            for(unsigned int gS=0;gS<groupSize; gS++)
            { printf ("\nGroup =%d\n",gS);
                for(int i=0; i<RADICES;i++)
                {
                    int index = ng * groupSize * RADICES + gS * RADICES + i;
                    int value = clHistData[ index ];
                    printf("%4x %4x, ",index, value);
                }
            }
        }
#endif
            if (swap == 0)
                AMP_permuteAscendingRadixNTemplate<T, RADIX>( ctl,
                                                        clInputData,
                                                        clHistData,
                                                        bits,
                                                        clSwapData,
                                                        szElements);
            else
                AMP_permuteAscendingRadixNTemplate<T, RADIX>( ctl,
                                                        clSwapData,
                                                        clHistData,
                                                        bits,
                                                        clInputData,
                                                        szElements);
#if BOLT_SORT_INL_DEBUG
if (swap == 0)
{
        printf("\n Printing swap data\n");
        for(unsigned int i=0; i<szElements;i+= RADICES)
        {
            for(int j =0;j< RADICES;j++)
                printf("%8x %8x, ",i+j,clSwapData[i+j]);
            printf("\n");
        }
}
else
{
        printf("\n Printing swap data\n");
        for(unsigned int i=0; i<szElements;i+= RADICES)
        {
            for(int j =0;j< RADICES;j++)
                printf("%8x %8x, ",i+j,clInputData[i+j]);
            printf("\n");
        }
}
#endif
            /*For swapping the buffers*/
            swap = swap? 0: 1;
        }
            /* Do descending for the signed bit
             * IN the case of radix 4 bits = 28 and in the case of radix 8 bits = 24
             */
             AMP_RadixSortHistogramSignedDescendingKernel<T, RADIX>(ctl,
                                                      clSwapData, /*this can be either array_view or array*/
                                                      clHistData,
                                                      clHistScanData,
                                                      bits,
                                                      szElements,
                                                      groupSize); // \TODO - i don't need gropu size

            detail::scan_enqueue( ctl, dvHistogramScanBuffer.begin( ), dvHistogramScanBuffer.end( ),dvHistogramScanBuffer.begin( ), 0, plus< T >( ) );

            AMP_scanLocalTemplate<T, RADIX>(ctl,
                                        clHistData,
                                        clHistScanData,
                                        szElements);

            AMP_permuteSignedDescendingRadixNTemplate<T, RADIX>( ctl,
                                                        clSwapData,
                                                        clHistData,
                                                        bits,
                                                        clInputData,
                                                        szElements);
            /*End of Ascending for the signed bit */
    }
    else
    {
        if(newBuffer == true)
        {
            concurrency::index<1> origin(orig_szElements);
            concurrency::array_view<T> dest = pLocalArrayView->section( origin, modified_ext - ext);
            //arrayview_type m_devMemoryAV( *m_devMemory );
            Concurrency::parallel_for_each( dest.extent, [dest]
                (Concurrency::index<1> idx) restrict(amp)
                {
                    dest[idx] = BOLT_INT_MIN;
                }
            );
        }
        /* Descending Sort */
        int bits=0;
        for(bits = 0; bits < ((sizeof(T) * 8) - RADIX); bits += RADIX)
        {
            if (swap == 0)
                AMP_RadixSortHistogramDescendingKernel<T, RADIX>(ctl,
                                                      clInputData, /*this can be either array_view or array*/
                                                      clHistData,
                                                      clHistScanData,
                                                      bits,
                                                      szElements,
                                                      groupSize); // \TODO - i don't need gropu size
            else
                AMP_RadixSortHistogramDescendingKernel<T, RADIX>(ctl,
                                                      clSwapData, /*this can be either array_view or array*/
                                                      clHistData,
                                                      clHistScanData,
                                                      bits,
                                                      szElements,
                                                      groupSize); // \TODO - i don't need gropu size

            detail::scan_enqueue( ctl, dvHistogramScanBuffer.begin( ),
                                  dvHistogramScanBuffer.end( ),dvHistogramScanBuffer.begin( ),
                                  0, plus< T >( ) );

            AMP_scanLocalTemplate<T, RADIX>(ctl,
                                        clHistData,
                                        clHistScanData,
                                        szElements);

            if (swap == 0)
                AMP_permuteDescendingRadixNTemplate<T, RADIX>( ctl,
                                                        clInputData,
                                                        clHistData,
                                                        bits,
                                                        clSwapData,
                                                        szElements);
            else
                AMP_permuteDescendingRadixNTemplate<T, RADIX>( ctl,
                                                        clSwapData,
                                                        clHistData,
                                                        bits,
                                                        clInputData,
                                                        szElements);
            /*For swapping the buffers*/
            swap = swap? 0: 1;
        }//End of For loop
            /* Do ascending for the signed bit
             * IN the case of radix 4 bits = 28 and in the case of radix 8 bits = 24
             */
             AMP_RadixSortHistogramSignedAscendingKernel<T, RADIX>(ctl,
                                                      clSwapData, /*this can be either array_view or array*/
                                                      clHistData,
                                                      clHistScanData,
                                                      bits,
                                                      szElements,
                                                      groupSize); // \TODO - i don't need gropu size

            detail::scan_enqueue( ctl, dvHistogramScanBuffer.begin( ),
                                  dvHistogramScanBuffer.end( ),dvHistogramScanBuffer.begin( ),
                                  0, plus< T >( ) );

            AMP_scanLocalTemplate<T, RADIX>(ctl,
                                        clHistData,
                                        clHistScanData,
                                        szElements);

            AMP_permuteSignedAscendingRadixNTemplate<T, RADIX>( ctl,
                                                        clSwapData,
                                                        clHistData,
                                                        bits,
                                                        clInputData,
                                                        szElements);
            /*End of Ascending for the signed bit */

    }
    if(newBuffer == true)
    {
        //std::cout << "New buffer was allocated So copying back the buffer\n";
        //dest = clInputData.section( ext );
        clInputData.section( ext ).copy_to( first.getBuffer( ) );
        first.getBuffer( ).synchronize( );
        delete pLocalArray;
    }
    delete pLocalArrayView;
    return;
}

/*!
* \brief This template function is a Bitonic Algorithm for one work item.
*        This is called by the the sort_enqueue routine for other than int's and unsigned int's
* \details Container - can be either array or array_view object and is dependant on the
*                      device_vector of the calling routine
*/
template <typename T, typename Container, typename StrictWeakOrdering>
void AMP_BitonicSortKernel(bolt::amp::control &ctl,
                           Container &A, /*this can be either array_view or array*/
                           int szElements,
                           unsigned int stage,
                           unsigned int passOfStage,
                           StrictWeakOrdering comp)
{
    concurrency::extent< 1 > globalSizeK0( szElements/2 );
    concurrency::tiled_extent< BITONIC_SORT_WGSIZE > tileK0 = globalSizeK0.tile< BITONIC_SORT_WGSIZE >();
    concurrency::accelerator_view av = ctl.getAccelerator().default_view;
    concurrency::parallel_for_each( av, tileK0,
    [
        A,
        passOfStage,
        stage,
        comp
    ]
    ( concurrency::tiled_index< BITONIC_SORT_WGSIZE > t_idx ) restrict(amp)
    {
        unsigned int  threadId = t_idx.global[ 0 ];
        unsigned int  pairDistance = 1 << (stage - passOfStage);
        unsigned int  blockWidth   = 2 * pairDistance;
        unsigned int  temp;
        unsigned int  leftId = (threadId % pairDistance)
                            + (threadId / pairDistance) * blockWidth;
        bool compareResult;

        unsigned int  rightId = leftId + pairDistance;

        T greater, lesser;
        T leftElement = A[leftId];
        T rightElement = A[rightId];

        unsigned int sameDirectionBlockWidth = 1 << stage;

        if((threadId/sameDirectionBlockWidth) % 2 == 1)
        {
            temp = rightId;
            rightId = leftId;
            leftId = temp;
        }

        compareResult = comp(leftElement, rightElement);

        if(compareResult)
        {
            greater = rightElement;
            lesser  = leftElement;
        }
        else
        {
            greater = leftElement;
            lesser  = rightElement;
        }
        A[leftId]  = lesser;
        A[rightId] = greater;
    } );// end of concurrency::parallel_for_each
}

template<typename DVRandomAccessIterator, typename StrictWeakOrdering>
typename std::enable_if<
    !(std::is_same< typename std::iterator_traits<DVRandomAccessIterator >::value_type, unsigned int >::value ||
      std::is_same< typename std::iterator_traits<DVRandomAccessIterator >::value_type,          int >::value)
                       >::type
sort_enqueue(bolt::amp::control &ctl, const DVRandomAccessIterator& first, const DVRandomAccessIterator& last,
const StrictWeakOrdering& comp)
{
    unsigned int numStages,stage,passOfStage;
    size_t  temp;
    typedef typename std::iterator_traits< DVRandomAccessIterator >::value_type T;

    int szElements = static_cast<int>( std::distance(first, last) );

    if(((szElements-1) & (szElements)) != 0)
    {
        sort_enqueue_non_powerOf2(ctl,first,last,comp);
        return;
    }
    /*if((szElements/2) < BITONIC_SORT_WGSIZE)
    {
        wgSize = (int)szElements/2;
    }*/
    auto&  A = first.getBuffer(); //( numElements, av );

    numStages = 0;
    for(temp = szElements; temp > 1; temp >>= 1)
        ++numStages;

    for(stage = 0; stage < numStages; ++stage)
    {
        for(passOfStage = 0; passOfStage < stage + 1; ++passOfStage)
        {
            //\TODO - can we remove this <T> specialization by getting the data type type from the container.
            AMP_BitonicSortKernel<T>(  ctl,
                                    A,
                                    szElements,
                                    stage,
                                    passOfStage,
                                    comp );
        }//end of for passStage = 0:stage-1
    }//end of for stage = 0:numStage-1

    return;
}// END of sort_enqueue

template<typename DVRandomAccessIterator, typename StrictWeakOrdering>
void sort_enqueue_non_powerOf2(control &ctl, const DVRandomAccessIterator& first, const DVRandomAccessIterator& last,
const StrictWeakOrdering& comp)
{
#if 0
    //std::cout << "The BOLT sort routine does not support non power of 2 buffer size. Falling back to CPU std::sort" << std ::endl;
    typedef typename std::iterator_traits< DVRandomAccessIterator >::value_type T;
    static boost::once_flag initOnlyOnce;
    size_t szElements = (size_t)(last - first);

    //Power of 2 buffer size
    // For user-defined types, the user must create a TypeName trait which returns the name of the class - note use of TypeName<>::get to retreive the name here.
    static std::vector< ::cl::Kernel > sortKernels;

    kernelParamsSort args( TypeName< T >::get( ), TypeName< DVRandomAccessIterator >::get( ), TypeName< T >::get( ),
    TypeName< StrictWeakOrdering >::get( ) );

    std::string typeDefinitions = cl_code + ClCode< T >::get( ) + ClCode< DVRandomAccessIterator >::get( );
    if( !boost::is_same< T, StrictWeakOrdering >::value)
    {
        typeDefinitions += ClCode< StrictWeakOrdering >::get( );
    }

    //Power of 2 buffer size
// For user-defined types, the user must create a TypeName trait which returns the name of the class
//  - note use of TypeName<>::get to retreive the name here.
    static std::vector< ::cl::Kernel > sortKernels;
        boost::call_once( initOnlyOnce, boost::bind( CallCompiler_Sort::constructAndCompileSelectionSort, &sortKernels,
            typeDefinitions, &args, &ctl) );

    size_t wgSize  = sortKernels[0].getWorkGroupInfo< CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE >( ctl.device( ), &l_Error );

    size_t totalWorkGroups = (szElements + wgSize)/wgSize;
    size_t globalSize = totalWorkGroups * wgSize;
    V_OPENCL( l_Error, "Error querying kernel for CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE" );

const ::cl::Buffer& in = first.getBuffer( );
    // ::cl::Buffer out(ctl.context(), CL_MEM_READ_WRITE, sizeof(T)*szElements);
    control::buffPointer out = ctl.acquireBuffer( sizeof(T)*szElements );

    ALIGNED( 256 ) StrictWeakOrdering aligned_comp( comp );
// ::cl::Buffer userFunctor(ctl.context(), CL_MEM_USE_HOST_PTR, sizeof( aligned_comp ), &aligned_comp );
control::buffPointer userFunctor = ctl.acquireBuffer( sizeof( aligned_comp ),
                                    CL_MEM_USE_HOST_PTR|CL_MEM_READ_ONLY, &aligned_comp );

    ::cl::LocalSpaceArg loc;
    loc.size_ = wgSize*sizeof(T);

    V_OPENCL( sortKernels[0].setArg(0, in), "Error setting a kernel argument in" );
V_OPENCL( sortKernels[0].setArg(1, first.gpuPayloadSize( ), &first.gpuPayload( ) ), "Error setting a kernel argument" );
V_OPENCL( sortKernels[0].setArg(2, *out), "Error setting a kernel argument out" );
V_OPENCL( sortKernels[0].setArg(3, *userFunctor), "Error setting a kernel argument userFunctor" );
V_OPENCL( sortKernels[0].setArg(4, loc), "Error setting kernel argument loc" );
V_OPENCL( sortKernels[0].setArg(5, static_cast<cl_uint> (szElements)), "Error setting kernel argument szElements" );
    {
            l_Error = ctl.commandQueue().enqueueNDRangeKernel(
                    sortKernels[0],
                    ::cl::NullRange,
                    ::cl::NDRange(globalSize),
                    ::cl::NDRange(wgSize),
                    NULL,
                    NULL);
            V_OPENCL( l_Error, "enqueueNDRangeKernel() failed for sort() kernel" );
            V_OPENCL( ctl.commandQueue().finish(), "Error calling finish on the command queue" );
    }

    wgSize  = sortKernels[1].getWorkGroupInfo< CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE >( ctl.device( ), &l_Error );

    V_OPENCL( sortKernels[1].setArg(0, *out), "Error setting a kernel argument in" );
    V_OPENCL( sortKernels[1].setArg(1, in), "Error setting a kernel argument out" );
V_OPENCL( sortKernels[1].setArg(2, first.gpuPayloadSize( ), &first.gpuPayload( ) ), "Error setting a kernel argument" );
V_OPENCL( sortKernels[1].setArg(3, *userFunctor), "Error setting a kernel argument userFunctor" );
V_OPENCL( sortKernels[1].setArg(4, loc), "Error setting kernel argument loc" );
V_OPENCL( sortKernels[1].setArg(5, static_cast<cl_uint> (szElements)), "Error setting kernel argument szElements" );
    {
            l_Error = ctl.commandQueue().enqueueNDRangeKernel(
                    sortKernels[1],
                    ::cl::NullRange,
                    ::cl::NDRange(globalSize),
                    ::cl::NDRange(wgSize),
                    NULL,
                    NULL);
            V_OPENCL( l_Error, "enqueueNDRangeKernel() failed for sort() kernel" );
            V_OPENCL( ctl.commandQueue().finish(), "Error calling finish on the command queue" );
    }
    // Map the buffer back to the host
    ctl.commandQueue().enqueueMapBuffer(in, true, CL_MAP_READ | CL_MAP_WRITE, 0/*offset*/, sizeof(T) * szElements, NULL, NULL, &l_Error );
    V_OPENCL( l_Error, "Error calling map on the result buffer" );
#endif
    return;
}// END of sort_enqueue_non_powerOf2

}//namespace bolt::cl::detail
}//namespace bolt::cl
}//namespace bolt

#endif

ext: inl

max_stars_repo_name: MattPD/Bolt

max_issues_repo_name: MattPD/Bolt

max_forks_repo_name: MattPD/Bolt


text:
//
// Copyright (c) 2016-2019 Vinnie Falco (vinnie dot falco at gmail dot com)
//
// Distributed under the Boost Software License, Version 1.0. (See accompanying
// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
//
// Official repository: https://github.com/boostorg/beast
//

#ifndef BOOST_BEAST_WEBSOCKET_DETAIL_PRNG_HPP
#define BOOST_BEAST_WEBSOCKET_DETAIL_PRNG_HPP

#include <boost/beast/core/detail/config.hpp>
#include <cstdint>
#include <limits>
#include <random>

namespace boost {
namespace beast {
namespace websocket {
namespace detail {

// Type-erased UniformRandomBitGenerator
// with 32-bit unsigned integer value_type
//
class prng
{
protected:
    virtual ~prng() = default;
    virtual prng& acquire() noexcept = 0;
    virtual void release() noexcept = 0;
    virtual std::uint32_t operator()() noexcept = 0;

public:
    using value_type = std::uint32_t;

    class ref
    {
        prng& p_;

    public:
        ref& operator=(ref const&) = delete;

        ~ref()
        {
            p_.release();
        }

        explicit
        ref(prng& p) noexcept
            : p_(p.acquire())
        {
        }

        ref(ref const& other) noexcept
            : p_(other.p_.acquire())
        {
        }

        value_type
        operator()() noexcept
        {
            return p_();
        }

        static
        value_type constexpr
        min() noexcept
        {
            return (std::numeric_limits<
                value_type>::min)();
        }

        static
        value_type constexpr
        max() noexcept
        {
            return (std::numeric_limits<
                value_type>::max)();
        }
    };
};

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

// Manually seed the prngs, must be called
// before acquiring a prng for the first time.
//
BOOST_BEAST_DECL
std::uint32_t const*
prng_seed(std::seed_seq* ss = nullptr);

// Acquire a PRNG using the no-TLS implementation
//
BOOST_BEAST_DECL
prng::ref
make_prng_no_tls(bool secure);

// Acquire a PRNG using the TLS implementation
//
#ifndef BOOST_NO_CXX11_THREAD_LOCAL
BOOST_BEAST_DECL
prng::ref
make_prng_tls(bool secure);
#endif

// Acquire a PRNG using the TLS implementation if it
// is available, otherwise using the no-TLS implementation.
//
BOOST_BEAST_DECL
prng::ref
make_prng(bool secure);

} // detail
} // websocket
} // beast
} // boost

#ifdef BOOST_BEAST_HEADER_ONLY
#include <boost/beast/websocket/detail/prng.ipp>
#endif

#endif

ext: hpp

max_stars_repo_name: Netis/packet-agent

max_issues_repo_name: Netis/packet-agent

max_forks_repo_name: Netis/packet-agent


text:
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include "main.h"
#include <Eigen/Geometry>

template<typename Scalar,int Size> void homogeneous(void)
{
  /* this test covers the following files:
     Homogeneous.h
  */

  typedef Matrix<Scalar,Size,Size> MatrixType;
  typedef Matrix<Scalar,Size,1, ColMajor> VectorType;

  typedef Matrix<Scalar,Size+1,Size> HMatrixType;
  typedef Matrix<Scalar,Size+1,1> HVectorType;

  typedef Matrix<Scalar,Size,Size+1>   T1MatrixType;
  typedef Matrix<Scalar,Size+1,Size+1> T2MatrixType;
  typedef Matrix<Scalar,Size+1,Size> T3MatrixType;

  VectorType v0 = VectorType::Random(),
             ones = VectorType::Ones();

  HVectorType hv0 = HVectorType::Random();

  MatrixType m0 = MatrixType::Random();

  HMatrixType hm0 = HMatrixType::Random();

  hv0 << v0, 1;
  VERIFY_IS_APPROX(v0.homogeneous(), hv0);
  VERIFY_IS_APPROX(v0, hv0.hnormalized());

  hm0 << m0, ones.transpose();
  VERIFY_IS_APPROX(m0.colwise().homogeneous(), hm0);
  VERIFY_IS_APPROX(m0, hm0.colwise().hnormalized());
  hm0.row(Size-1).setRandom();
  for(int j=0; j<Size; ++j)
    m0.col(j) = hm0.col(j).head(Size) / hm0(Size,j);
  VERIFY_IS_APPROX(m0, hm0.colwise().hnormalized());

  T1MatrixType t1 = T1MatrixType::Random();
  VERIFY_IS_APPROX(t1 * (v0.homogeneous().eval()), t1 * v0.homogeneous());
  VERIFY_IS_APPROX(t1 * (m0.colwise().homogeneous().eval()), t1 * m0.colwise().homogeneous());

  T2MatrixType t2 = T2MatrixType::Random();
  VERIFY_IS_APPROX(t2 * (v0.homogeneous().eval()), t2 * v0.homogeneous());
  VERIFY_IS_APPROX(t2 * (m0.colwise().homogeneous().eval()), t2 * m0.colwise().homogeneous());

  VERIFY_IS_APPROX((v0.transpose().rowwise().homogeneous().eval()) * t2,
                    v0.transpose().rowwise().homogeneous() * t2);
                    m0.transpose().rowwise().homogeneous().eval();
  VERIFY_IS_APPROX((m0.transpose().rowwise().homogeneous().eval()) * t2,
                    m0.transpose().rowwise().homogeneous() * t2);

  T3MatrixType t3 = T3MatrixType::Random();
  VERIFY_IS_APPROX((v0.transpose().rowwise().homogeneous().eval()) * t3,
                    v0.transpose().rowwise().homogeneous() * t3);
  VERIFY_IS_APPROX((m0.transpose().rowwise().homogeneous().eval()) * t3,
                    m0.transpose().rowwise().homogeneous() * t3);

  // test product with a Transform object
  Transform<Scalar, Size, Affine> aff;
  Transform<Scalar, Size, AffineCompact> caff;
  Transform<Scalar, Size, Projective> proj;
  Matrix<Scalar, Size, Dynamic>   pts;
  Matrix<Scalar, Size+1, Dynamic> pts1, pts2;

  aff.affine().setRandom();
  proj = caff = aff;
  pts.setRandom(Size,internal::random<int>(1,20));

  pts1 = pts.colwise().homogeneous();
  VERIFY_IS_APPROX(aff  * pts.colwise().homogeneous(), (aff  * pts1).colwise().hnormalized());
  VERIFY_IS_APPROX(caff * pts.colwise().homogeneous(), (caff * pts1).colwise().hnormalized());
  VERIFY_IS_APPROX(proj * pts.colwise().homogeneous(), (proj * pts1));

  VERIFY_IS_APPROX((aff  * pts1).colwise().hnormalized(),  aff  * pts);
  VERIFY_IS_APPROX((caff * pts1).colwise().hnormalized(), caff * pts);

  pts2 = pts1;
  pts2.row(Size).setRandom();
  VERIFY_IS_APPROX((aff  * pts2).colwise().hnormalized(), aff  * pts2.colwise().hnormalized());
  VERIFY_IS_APPROX((caff * pts2).colwise().hnormalized(), caff * pts2.colwise().hnormalized());
  VERIFY_IS_APPROX((proj * pts2).colwise().hnormalized(), (proj * pts2.colwise().hnormalized().colwise().homogeneous()).colwise().hnormalized());
}

void test_geo_homogeneous()
{
  for(int i = 0; i < g_repeat; i++) {
    CALL_SUBTEST_1(( homogeneous<float,1>() ));
    CALL_SUBTEST_2(( homogeneous<double,3>() ));
    CALL_SUBTEST_3(( homogeneous<double,8>() ));
  }
}

ext: cpp

max_stars_repo_name: subond/tools

max_issues_repo_name: subond/tools

max_forks_repo_name: subond/tools


text:
#ifndef PLOT_WINDOW_HPP
#define PLOT_WINDOW_HPP

// Qt headers
#include <QMainWindow>
#include <QList>
#include <QLayout>
#include <QSettings>
#include <QDockWidget>
#include <QMessageBox>
#include <QFileDialog>
#include <QtDebug>
//#include <QtWidgets>

// C headers
#include <cstdlib>

// C++ headers
#include <iostream>
#include <iomanip>

// Boost headers
//#include <boost/property_tree/ptree.hpp>
//#include <boost/property_tree/ini_parser.hpp>
//#include <boost/tokenizer.hpp>

// Local headers
#include "ode_pathtracer.hpp"
#include "runge_kutta_stepper.hpp"
#include "render_view.hpp"
#include "simulation_loop.hpp"
#include "label_dock_widget.hpp"

// OUT and IN can be redefined as a filestream
// to enable direct file input/output
// In MinGW already defined in windef.h
#ifdef OUT
#undef OUT
#endif // OUT
#define OUT qDebug()
#ifdef IN
#undef IN
#endif
#define IN  std::cin
#define ERROUT qDebug()
#define LF  std::endl
#define ABS std::abs

// Ini file sections
#define SECTION_NAMES     "names"
#define SECTION_PARAM_EQ  "parameter equations"
#define SECTION_VAR_DERIV "variable derivations"
#define SECTION_VAR_INIT  "variable initial"
#define SECTION_TIME      "time"
#define SECTION_PLOT      "plot"

namespace Ui {
class PlotWindow;
}

class PlotWindow : public QMainWindow
{
   Q_OBJECT

public:
   explicit PlotWindow( QWidget *parent = 0 );
   ~PlotWindow();

signals:
   void updateParamLabels( PointValues values );

public slots:
   void updateView( PointValues newPoint );

private:
   Ui::PlotWindow *ui;

   // simulation objects
   RungeKuttaStepper *stepper = NULL;
   SimulationLoop    *simulation = NULL;
   QList<PointValues> pointPath;

   // actions
   QAction *runAction;
   QAction *exitAction;
   QAction *openProblemAction;
   QAction *closeProblemAction;

   // gui elements
   QGridLayout *mainLayout = NULL;
   QDockWidget *labelDock = NULL;
   QList<RenderView *> views;

   // Labels
   QStringList labelNames;
   LabelDockWidget *dockWidget = NULL;

   // Names
   QStringList paramNames;
   QStringList varNames;

   // Equations
   EquationVector   paramRules;
   DerivationVector varRules;

   // Initial values
   PointValues initialValues;

   // Time parameters
   double dt;

   // Plot parameters
   QString plotTransformX;
   QString plotTransformY;
   QRect   plotViewport;
   int plotMaxFPS;
   int plotSkip;
   int plotMaxPathSegments;

   void ThrowError( QString msg );
   void inputData( const QString filename );
   QStringList tokenizeString( QString &str );
   void toggleSimulationRun( bool toggled );
   void exitProgram( bool checked = false );
   void openProblem( bool checked = false );
   void loadProblem( const QString filename );
   void closeProblem( bool checked = false );

   template <class T> T readEntry(
      QSettings *inputFile
    , QString section
    , QString name
    , T defaultValue
   ){
      T value;

      // show warning if key is missing
      if( !(inputFile->contains( section+"/"+name )) ){
         qDebug() << "WARNING: Key \"" << name << "\" in section [" << section << "] not found.\n"
                  << "         Assuming " << name << " = " << defaultValue << "\n";
      }

      // get value
      QVariant var = inputFile->value( section+"/"+name, defaultValue );
      value = var.value<T>();

      //qDebug() << "read value = " << value;

      return value;
   }
};

#endif // PLOT_WINDOW_HPP


text:
#pragma once

#include <ros/ros.h>
#include <tf/transform_listener.h>
#include <tf_conversions/tf_eigen.h>
#include <Eigen/Dense>

#include <sensor_msgs/Image.h>
#include <sensor_msgs/PointCloud2.h>
#include <sensor_msgs/PointCloud.h>
#include <sensor_msgs/point_cloud_conversion.h>
#include <sensor_msgs/CameraInfo.h>

#include <message_filters/subscriber.h>
#include <message_filters/synchronizer.h>
#include <message_filters/time_synchronizer.h>
#include <message_filters/sync_policies/approximate_time.h>

#include <image_transport/image_transport.h>
#include <opencv2/highgui/highgui.hpp>
#include <cv_bridge/cv_bridge.h>

#include "SegmentationMapping/ImageLabelDistribution.h"

#include <chrono>
#include <iostream>
#include <fstream>
#include <opencv2/core/matx.hpp>
#include <unordered_map>
#include <tuple>
#include <math.h>

using namespace std;
using namespace Eigen;
using namespace std::chrono;

namespace SegmentationMapping {

  class LidarProjection {
  public:
    LidarProjection()
    {
      ros::NodeHandle nh("~");
      get_params(nh);

      pcl_pub = nh.advertise<sensor_msgs::PointCloud>(Pub_Topic_, 5);
      pcl_sub.subscribe(nh, cloud_topic_, 5);
      img_sub.subscribe(nh, image_topic_, 10);
      cam_sub.subscribe(nh, cam_info_, 10);
      dist_sub.subscribe(nh, dist_info_, 1);
      sync_.reset(new ApproximateSync(ApproximateSyncPolicy(300), pcl_sub, img_sub,
                           cam_sub, dist_sub));
      //sync_.reset(new ExactSync(ExactSyncPolicy(1000), pcl_sub, img_sub, cam_sub,
      //                          dist_sub ));
      sync_->registerCallback(boost::bind(&LidarProjection::callback,
                                          this, _1, _2, _3, _4));

      label2color[2]  =std::make_tuple(250, 250, 250 ); // road
      //label2color[3]  =std::make_tuple(128, 64,  128 ); // sidewalk
      label2color[3]  =std::make_tuple(250, 250,  250 ); // sidewalk
      label2color[5]  =std::make_tuple(250, 128, 0   ); // building
      label2color[10] =std::make_tuple(192, 192, 192 ); // pole
      label2color[12] =std::make_tuple(250, 250, 0   ); // sign
      label2color[6]  =std::make_tuple(0  , 100, 0   ); // vegetation
      label2color[4]  =std::make_tuple(128, 128, 0   ); // terrain
      label2color[13] =std::make_tuple(135, 206, 235 ); // sky
      label2color[1]  =std::make_tuple( 30, 144, 250 ); // water
      label2color[8]  =std::make_tuple(220, 20,  60  ); // person
      label2color[7]  =std::make_tuple( 0, 0,142     ); // car
      label2color[9]  =std::make_tuple(119, 11, 32   ); // bike
      label2color[11] =std::make_tuple(123, 104, 238 ); // stair
      label2color[0]  =std::make_tuple(255, 255, 255 ); // background

    }

    void get_params(ros::NodeHandle &nh_);
    // Getting rostopics and extrinsic matrix from launch file.

    Matrix4f vector2matrix(vector<float> input);
    // Converting from [x y z qx qy qz qw] into transformation matrix.

    MatrixXf removePoints(const MatrixXf& matrix, vector<int> true_rows);

    MatrixXf dist_info(const MatrixXf& matrix, vector<int> true_rows);
    // vector true_rows contains the points idx we want to keep.
    // Discarding points from matrix.

    vector<int> depth_color(MatrixXf& depth);

    MatrixXf point_in_frame(const MatrixXf& point, const MatrixXf& dists);

    MatrixXf velo_point_filter(MatrixXf& velo_data);
    // Getting points in camera frame.

    void set_pointcloud_matrix(const sensor_msgs::PointCloud2ConstPtr& msg);

    void set_intrinsic(const sensor_msgs::CameraInfoConstPtr& cam);

    void visualize_results(const sensor_msgs::ImageConstPtr& img);

    void set_label_pcl(                    const sensor_msgs::ImageConstPtr& img,
                                           const ImageLabelDistribution::ConstPtr& info);

    void callback(const sensor_msgs::PointCloud2ConstPtr& msg,
                  const sensor_msgs::ImageConstPtr& img,
                  const sensor_msgs::CameraInfoConstPtr& cam,
                  const ImageLabelDistribution::ConstPtr& info);

  private:
    int Lidar_max_range;

    string Pub_Topic_;
    string cloud_topic_;
    string image_topic_;
    string cam_info_;
    string dist_info_;

    Matrix4f Extrinsic;
    sensor_msgs::PointCloud Cloud;
    sensor_msgs::PointCloud DistriCloud;
    MatrixXf CloudMat;
    MatrixXf Intrinsic;
    vector<int> color_info;
    vector<int> pcl_info;
    std::unordered_map<int, std::tuple<uint8_t, uint8_t, uint8_t>> label2color;

    message_filters::Subscriber<sensor_msgs::PointCloud2> pcl_sub;
    message_filters::Subscriber<sensor_msgs::Image> img_sub;
    message_filters::Subscriber<sensor_msgs::CameraInfo> cam_sub;
    message_filters::Subscriber<sensor_msgs::Image> label_sub;
    message_filters::Subscriber<ImageLabelDistribution> dist_sub;

    typedef message_filters::sync_policies::ExactTime<
      sensor_msgs::PointCloud2,
      sensor_msgs::Image,
      sensor_msgs::CameraInfo,
      ImageLabelDistribution> ExactSyncPolicy;

    typedef message_filters::sync_policies::ApproximateTime<
      sensor_msgs::PointCloud2,
      sensor_msgs::Image,
      sensor_msgs::CameraInfo,
      ImageLabelDistribution> ApproximateSyncPolicy;
    typedef message_filters::Synchronizer<ApproximateSyncPolicy> ApproximateSync;
    typedef message_filters::Synchronizer<ExactSyncPolicy> ExactSync;
    boost::shared_ptr<ApproximateSync> sync_;
    //boost::shared_ptr<ExactSync> sync_;
    ros::Publisher pcl_pub;

    const float MIN_DIS = 0;
    const float MAX_DIS = 120;
    int WIDTH, HEIGHT;
    // ofstream file;
  public:
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
  };

  void LidarProjection::get_params(ros::NodeHandle &nh_){
    nh_.getParam("lidar_max_range", Lidar_max_range);
    nh_.getParam("publish_to", Pub_Topic_);
    nh_.getParam("cloud_topic", cloud_topic_);
    nh_.getParam("image_topic", image_topic_);
    nh_.getParam("camera_info", cam_info_);
    nh_.getParam("distribution", dist_info_);
    nh_.getParam("image_width", WIDTH);
    nh_.getParam("image_height", HEIGHT);

    //vector<float> ImuLidar;
    //vector<float> ImuCamera;
    //nh_.getParam("imu2lidar", ImuLidar);
    //nh_.getParam("imu2camera", ImuCamera);
    vector<float> Extrinsic_param;
    nh_.getParam("extrinsic_mat", Extrinsic_param);

    // MatrixXf ImuLidarMat = vector2matrix(ImuLidar);
    // MatrixXf ImuCameraMat = vector2matrix(ImuCamera);
    //Extrinsic = ImuLidarMat * ImuCameraMat.inverse();
    Extrinsic.setIdentity();
    for (int i = 0; i < Extrinsic_param.size(); i++)
      Extrinsic(i / 4, i % 4) = Extrinsic_param[i];
    std::cout<<"Extrinsic calibration mat is \n"<<Extrinsic;
  }

  Matrix4f LidarProjection::vector2matrix(vector<float> input){
    // input is x y z qx qy qz qw;
    Matrix4f TransFull;
    Vector3f translation(input[0], input[1], input[2]);
    Quaternionf q;
    q.x() = input[3];
    q.y() = input[4];
    q.z() = input[5];
    q.w() = input[6];

    TransFull.setIdentity();
    TransFull.block<3, 3>(0, 0) = q.toRotationMatrix();
    TransFull.block<3, 1>(0, 3) = translation;

    return TransFull;
  }

  MatrixXf LidarProjection::removePoints(const MatrixXf& matrix,
                                         vector<int> true_rows){
    MatrixXf new_matrix;
    new_matrix.resize(2, true_rows.size());

    int flag_insert = 0;
    size_t flag_true = 0;
    for (size_t i = 0; i < matrix.cols(); i++){
      // cout << i << " ";
      if (i == true_rows[flag_true]){
        assert(flag_true <= matrix.cols());
        assert(flag_insert <= true_rows.size());
        if ((flag_insert >= true_rows.size())||
            (flag_true >= new_matrix.cols())){
          break;
        }
        new_matrix(0, flag_insert) = matrix(0, i)/matrix(2, i);
        new_matrix(1, flag_insert) = matrix(1, i)/matrix(2, i);
        flag_insert += 1;
        flag_true += 1;
      }
      else{
        continue;
      }
    }
    return new_matrix;
  }

  MatrixXf LidarProjection::dist_info(const MatrixXf& matrix,
                                      vector<int> true_rows){
    MatrixXf new_matrix;
    new_matrix.resize(1, true_rows.size());

    int flag_insert = 0;
    size_t flag_true = 0;
    for (size_t i = 0; i < matrix.rows(); i++){

      if (i == true_rows[flag_true]){
        assert(flag_true <= matrix.rows());
        assert(flag_insert <= true_rows.size());
        if ((flag_insert >= true_rows.size()) ||
            (flag_true >= new_matrix.cols())){
          break;
        }
        new_matrix(0, flag_insert) = matrix(i, 0);
        flag_insert += 1;
        flag_true += 1;
      }
      else{
        continue;
      }
    }
    return new_matrix;
  }

  vector<int> LidarProjection::depth_color(MatrixXf& depth){
    // cout << depth.rows() << " " << depth.cols() << endl;
    vector<int> color;
    for (int i = 0; i < depth.cols(); i++){
      if (depth(0, i) < MIN_DIS)
        color.push_back(MIN_DIS);
      else if (depth(0, i) > MAX_DIS)
        color.push_back(round(MAX_DIS * MAX_DIS/70));
      else
        color.push_back(round(depth(0, i) * MAX_DIS/70));
    }

    return color;
  }

  MatrixXf LidarProjection::point_in_frame(const MatrixXf& point,
                                           const MatrixXf& dists){

    MatrixXf pointmat;
    pointmat.resize(4, point.cols());
    pointmat.block(0, 0, 3, point.cols()) = point;
    pointmat.block(3, 0, 1, point.cols()) = MatrixXf::Constant(
                                                               1, point.cols(), 1);

    // pointmat.resize(point.rows(), 4);
    // pointmat.block(0, 0, point.rows(), 3) = point;
    // pointmat.block(0, 3, point.rows(), 1) = MatrixXf::Constant(
    //                                          point.rows(), 1, 1);
    // pointmat.transposeInPlace();
    pointmat = Intrinsic * Extrinsic * pointmat;
    // cout << "pointmat: " << pointmat.rows() << " " << pointmat.cols() << endl;
    vector<int> v;
    for (size_t i = 0; i < point.cols(); i++){
      if (pointmat(2, i) > 0){
        if (pointmat(0, i)/pointmat(2, i) >= 0
            && pointmat(0, i)/pointmat(2, i) < WIDTH-1
            && pointmat(1, i)/pointmat(2, i) >= 0
            && pointmat(1, i)/pointmat(2, i) < HEIGHT-1){
          v.push_back(i);
        }
      }
    }
    // cout << "vector: " << v.size() << endl;
    MatrixXf inFrame = removePoints(pointmat, v);
    // cout << inFrame.rows() << " " << inFrame.cols() << endl;

    // MatrixXf new_dist = dist_info(dists, v);
    // cout << new_dist.rows() << " " << new_dist.cols() << endl;

    // color_info = depth_color(new_dist);
    pcl_info = v;
    // cout << inFrame.cols() << endl;
    return inFrame;
  }

  MatrixXf LidarProjection::velo_point_filter(MatrixXf& velo_data){
    // cout << velo_data.rows() << " " << velo_data.cols() << endl;
    MatrixXf x = velo_data.block(0, 0, 1, velo_data.cols());
    MatrixXf y = velo_data.block(1, 0, 1, velo_data.cols());
    MatrixXf z = velo_data.block(2, 0, 1, velo_data.cols());

    MatrixXf dist = (x.array().pow(2) +
                     y.array().pow(2) + z.array().pow(2)).array().pow(0.5);
    // cout << dist.cols() << endl;
    // cout << velo_data.rows() << " " << velo_data.cols() << endl;
    MatrixXf points = point_in_frame(velo_data, dist);

    return points;
  }

  void LidarProjection::set_pointcloud_matrix(
                                              const sensor_msgs::PointCloud2ConstPtr& msg){
    sensor_msgs::convertPointCloud2ToPointCloud(*msg, Cloud);
    // CloudMat.resize(3, Cloud.points.size());

    vector<int> in_range;
    int count = 0;
    for (int i = 0; i < Cloud.points.size(); i++){
      if (sqrt(pow(Cloud.points[i].x, 2) +
               pow(Cloud.points[i].y, 2) + pow(Cloud.points[i].z, 2)) <= Lidar_max_range){
        in_range.push_back(i);
        count += 1;
      }
    }
    CloudMat.resize(3, count);
    // for (int i = 0; i < Cloud.points.size(); i++){
    //   CloudMat(0, i) = Cloud.points[i].x;
    //   CloudMat(1, i) = Cloud.points[i].y;
    //   CloudMat(2, i) = Cloud.points[i].z;
    // }
    for (auto i = 0; i < in_range.size(); i++){
      CloudMat(0, i) = Cloud.points[in_range[i]].x;
      CloudMat(1, i) = Cloud.points[in_range[i]].y;
      CloudMat(2, i) = Cloud.points[in_range[i]].z;
    }
  }

  void LidarProjection::set_intrinsic(
                                      const sensor_msgs::CameraInfoConstPtr& cam){
    Intrinsic.resize(12, 1);
    for (size_t j = 0; j < 12; j++){
      Intrinsic(j, 0) = cam->P[j];
    }
    Intrinsic.resize(4, 3);
    Intrinsic.transposeInPlace();
  }

  void LidarProjection::visualize_results(
                                          const sensor_msgs::ImageConstPtr& img){
    cv::Mat image = cv_bridge::toCvCopy(img, "bgr8")->image;
    cv::cvtColor(image, image, cv::COLOR_BGR2HSV);
    for (size_t j = 0; j < color_info.size(); j++){
      cv::circle(image, cv::Point(round(CloudMat(0, j)),
                                  round(CloudMat(1, j))), 2,
                 cv::Scalar(color_info[j], 255, 255), -1);
    }
    cv::cvtColor(image, image, cv::COLOR_HSV2BGR);
    // cout << image.rows << " " << image.cols << endl;
    cv::imshow("image", image);
    cv::waitKey(5);
  }

  void LidarProjection::set_label_pcl(
                                      const sensor_msgs::ImageConstPtr& img,
                                      const ImageLabelDistribution::ConstPtr& info)
  {
    int stride1 = info->distribution.layout.dim[1].stride;
    int stride2 = info->distribution.layout.dim[2].stride;
    int channel = info->distribution.layout.dim[2].size;
    DistriCloud.points.resize(CloudMat.cols());
    DistriCloud.channels.resize(channel+7);

    cv::Mat image = cv_bridge::toCvCopy(img, "rgb8")->image;
    cv::cvtColor(image, image, cv::COLOR_BGR2RGB);

    for (int i = 0; i < channel+7; i++){
      DistriCloud.channels[i].values.resize(CloudMat.cols());
    }

    for (size_t j = 0; j < pcl_info.size(); j++){
      // std::cout<<j<<", \nnew "<<          info->distribution.data[4300799]<<std::endl;
      DistriCloud.points[j].x = Cloud.points[pcl_info[j]].x;
      DistriCloud.points[j].y = Cloud.points[pcl_info[j]].y;
      DistriCloud.points[j].z = Cloud.points[pcl_info[j]].z;
      int col = round(CloudMat(0, j));
      int row = round(CloudMat(1, j));
      float max_prob = 0;
      int max_l = -1;

      for (int k = 0; k < channel; k++){
        // std::cout<<"row "<<row<<", col "<<col<<", index "<<stride1 * row + stride2 * col + k<<std::endl;
        // std::cout<<info->distribution.data[stride1 * row + stride2 * col + k]<<"\n";
        DistriCloud.channels[k+7].values[j] =
          info->distribution.data[stride1 * row + stride2 * col + k];

        if (info->distribution.data[stride1 * row + stride2 * col + k] > max_prob) {
          max_prob = info->distribution.data[stride1 * row + stride2 * col + k];

          max_l = k;
        }
      }
      // std::cout<<"max_l: "<<max_l<<std::endl;
      DistriCloud.channels[0].values[j] = max_l;
      DistriCloud.channels[1].values[j] = std::get<0>(label2color[max_l]);
      DistriCloud.channels[2].values[j] = std::get<1>(label2color[max_l]);
      DistriCloud.channels[3].values[j] = std::get<2>(label2color[max_l]);
      cv::Vec3b intensity = image.at<cv::Vec3b>(row, col);
      uchar blue = intensity.val[0];
      uchar green = intensity.val[1];
      uchar red = intensity.val[2];
      DistriCloud.channels[4].values[j] = red;
      DistriCloud.channels[5].values[j] = green;
      DistriCloud.channels[6].values[j] = blue;
      // std::cout<<"j"<<j<<std::endl;

      // debuging with visualization.
      // cv::circle(image, cv::Point(round(CloudMat(0, j)),
      //            round(CloudMat(1, j))), 2,
      //            cv::Scalar(std::get<0>(label2color[max_l]),
      //                       std::get<1>(label2color[max_l]),
      //                       std::get<2>(label2color[max_l])), -1);
    }
    // std::cout<<"Finish generate distribution cloud\n";
    // cv::cvtColor(image, image, cv::COLOR_RGB2BGR);
    // cv::imshow("image", image);
    // cv::waitKey(5);
  }

  void LidarProjection::callback(const sensor_msgs::PointCloud2ConstPtr& msg,
                                 const sensor_msgs::ImageConstPtr& img,
                                 const sensor_msgs::CameraInfoConstPtr& cam,
                                 const ImageLabelDistribution::ConstPtr& info)
  {
    ROS_INFO("lidar callback!");
    high_resolution_clock::time_point t1 = high_resolution_clock::now();
    ros::Time curr_t = ros::Time::now();
    ROS_INFO_STREAM("Callback starts at time "<<uint32_t(curr_t.toSec())<<". "<<(uint32_t)curr_t.toNSec() <<", the img time is "<<(uint32_t)img->header.stamp.toSec()<< " "<<(uint32_t)img->header.stamp.toNSec()<<", the lidar time is "<< (uint32_t)msg->header.stamp.toSec()<< " "<<(uint32_t)msg->header.stamp.toNSec()<<", the segmentation time is "<< (uint32_t)info->header.stamp.toSec()<< " "<<(uint32_t)info->header.stamp.toNSec());

    set_pointcloud_matrix(msg);

    set_intrinsic(cam);

    CloudMat = velo_point_filter(CloudMat);
    // cout << "CloudMat:" << CloudMat.rows() << " " << CloudMat.cols() << endl;
    // Visulization
    //visualize_results(img);

    set_label_pcl(img, info);
    DistriCloud.header = msg->header;
    DistriCloud.header.frame_id = "/velodyne_actual";
    pcl_pub.publish(DistriCloud);

    high_resolution_clock::time_point t2 = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>( t2 - t1 ).count();
    cout << "1 Frame time: " << duration << endl;
    cout << DistriCloud.header << endl;
    ros::Time curr_t2 = ros::Time::now();
    ROS_DEBUG_STREAM("Callback starts at time "<<uint32_t(curr_t.toSec())<<". "<<(uint32_t)curr_t2.toNSec() );

    ROS_INFO("Published new cloud.");

    // file << duration << "  " << Cloud.points.size() << "  " << pcl_info.size() << endl;
  }

}

max_stars_repo_name: UMich-BipedLab/SegmentationMapping

max_issues_repo_name: crankler/SegmentationMapping

max_forks_repo_name: crankler/SegmentationMapping


text:
#define BOOST_TEST_DYN_LINK

#include <boost/test/unit_test.hpp>
#include <algorithm>
#include "main.cpp"

BOOST_AUTO_TEST_SUITE(GeneralTests)

BOOST_AUTO_TEST_CASE(GetInfo) {

	MyFixture fx;
	getinfo_t info;

	NO_THROW(info = fx.btc.getinfo());
	BOOST_REQUIRE(info.protocolversion >= 70002);

	#ifdef VERBOSE
	std::cout <<  "=== getinfo ===" << std::endl;
	std::cout << "Version: " << info.version << std::endl;
	std::cout << "Protocol Version: " << info.protocolversion << std::endl;
	std::cout << "Wallet Version: " << info.walletversion << std::endl;
	std::cout << "Balance: " << info.balance << std::endl;
	std::cout << "Blocks: " << info.blocks << std::endl;
	std::cout << "Timeoffset: " << info.timeoffset << std::endl;
	std::cout << "Connections: " << info.connections << std::endl;
	std::cout << "Proxy: " << info.proxy << std::endl;
	std::cout << "Difficulty: " << info.difficulty << std::endl;
	std::cout << "Testnet: " << info.testnet << std::endl;
	std::cout << "Keypoololdest: " << info.keypoololdest << std::endl;
	std::cout << "Keypool size: " << info.keypoolsize << std::endl;
	std::cout << "Paytxfee: " << info.paytxfee << std::endl;
	std::cout << "Unlocked_until: " << info.unlocked_until << std::endl;
	std::cout << "Errors: " << info.errors << std::endl << std::endl;
	#endif
}

BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: BankLedger/safe-rpc-api

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: BankLedger/safe-rpc-api

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: BankLedger/safe-rpc-api

max_forks_repo_licenses: MIT Unlicense