Median.hpp

Go to the documentation of this file.

/*--------------------------------------------------------------------------------------------
Copyright (c) 2019, NDEVR LLC
tyler.parke@ndevr.org
  __    __   ____     _____ __     __  _______
 |  \  |  | |  __ \  |  ___|\ \   / / |   __  \
 |   \ |  | | |  \ \ | |___  \ \ / /  |  |__)  |
 |  . \|  | | |__/ / | |___   \ V /   |   _   /
 |  |\    |_|_____/__|_____|___\_/____|  | \  \
 |__| \__________________________________|  \__\
 
Subject to the terms of the Enterprise+ Agreement, NDEVR hereby grants 
Licensee a limited, non-exclusive, non-transferable, royalty-free license 
(without the right to sublicense) to use the API solely for the purpose of 
Licensee's internal development efforts to develop applications for which 
the API was provided.
 
The above copyright notice and this permission notice shall be included in all 
copies or substantial portions of the Software.
 
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE
FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
 
Library: Base
File: Median
Included in API: True
Author(s): Tyler Parke
 *-----------------------------------------------------------------------------------------**/
#pragma once
#include <NDEVR/Buffer.h>
#include <NDEVR/BinaryHeap.h>
namespace NDEVR
{
    /**--------------------------------------------------------------------------------------------------
    \brief A dummy class for including optimized median functions
    **/
    class Median
    {};
    // return the median value in a vector 
    template<class t_container_type>
    uint04 median(const t_container_type& elements, Buffer<uint04>& indices, const uint04 start, const uint04 end, uint01 dimension)
    {
        uint04* middle = indices.begin(start + ((end - start) / 2));
        std::nth_element(indices.begin(start), middle, indices.begin(end), [&elements, dimension](uint04 i, uint04 j)
        {
            return elements[i].center()[dimension] < elements[j].center()[dimension];
        });
        return start + ((end - start) / 2);
    }
 
    // return the median value in a vector 
    template<class t_node_type, class t_container_type>
    uint04 dynamicHeapMedian(uint04 size, const Buffer<t_node_type>& points, const t_container_type& indices, const uint04 start, const uint04 step, uint01 dimension)
    {
        const sint04 half_heap = cast<sint04>(size / 2);
        const sint04 max_heap = cast<sint04>(size / 2);
        MinHeap<t_node_type, uint04> left;
        MaxHeap<t_node_type, uint04> right;
        sint04 nLeft = 1;
        sint04 nRight = 1;
        sint04 index = start;
        // pick first value as median candidate
        sint04 median = index;
        index += step;
 
        for (;;)
        {
            // get next value
            sint04 val = index;
            sint04 child = nLeft++, parent = (child - 1) / 2;
            index += step;
            // if value is smaller than median, append to left heap
            if (points[indices[val]][dimension] < points[indices[median]][dimension])
            {
                // move biggest value to the heap top
                left.insert(points[indices[val]][dimension], val);
 
                // if left heap is full
                if (left.size() == max_heap)
                {
                    // for each remaining value
                    for (; nRight != max_heap; nRight++)
                    {
                        // get next value
                        val = index;
                        index += step;
                        // if value is to be inserted in the left heap
                        if (points[indices[val]][dimension] < points[indices[median]][dimension])
                        {
 
                            child = points[indices[left[2]]][dimension] > points[indices[left[1]]][dimension] ? 2 : 1;
                            if (points[indices[val]][dimension] >= points[indices[left[child]]][dimension])
                                median = val;
                            else
                            {
                                median = left[child];
                                parent = child;
                                child = parent * 2 + 1;
                                while (child < max_heap)
                                {
                                    if (child < half_heap
                                        && points[indices[left[child + 1]]][dimension] > points[indices[left[child]]][dimension])
                                        ++child;
                                    if (points[indices[val]][dimension] >= points[indices[left[child]]][dimension])
                                        break;
                                    left[parent] = left[child];
                                    parent = child;
                                    child = parent * 2 + 1;
                                }
                                left[parent] = val;
                            }
                        }
                    }
                    return median;
                }
            }
 
            // else append to right heap
            else
            {
                // move smallest value to the heap top
                sint04 child = nRight++, parent = (child - 1) / 2;
                while (parent && points[indices[val]][dimension] < points[indices[right[parent]]][dimension])
                {
                    right[child] = right[parent];
                    child = parent;
                    parent = (parent - 1) / 2;
                }
                right[child] = val;
 
                // if right heap is full
                if (nRight == max_heap)
                {
                    // for each remaining value
                    for (; nLeft != max_heap; nLeft++)
                    {
                        // get next value
                        val = index;
                        index += step;
                        // if value is to be inserted in the right heap
                        if (points[indices[val]][dimension] > points[indices[median]][dimension])
                        {
                            child = points[indices[right[2]]][dimension] < points[indices[right[1]]][dimension] ? 2 : 1;
                            if (points[indices[val]][dimension] <= points[indices[right[child]]][dimension])
                                median = val;
                            else
                            {
                                median = right[child];
                                parent = child;
                                child = parent * 2 + 1;
                                while (child < max_heap)
                                {
                                    if (child < half_heap
                                        && points[indices[right[child + 1]]][dimension] < points[indices[right[child]]][dimension])
                                        ++child;
                                    if (points[indices[val]][dimension] <= points[indices[right[child]]][dimension])
                                        break;
                                    right[parent] = right[child];
                                    parent = child;
                                    child = parent * 2 + 1;
                                }
                                right[parent] = val;
                            }
                        }
                    }
                    return median;
                }
            }
        }
    }
    /**--------------------------------------------------------------------------------------------------
    \brief Return the median value in a Buffer of points based on a given dimension. 
    **/
    template<uint01 t_bucket_size, class t_buffer_type, class t_container_type>
    uint04 StaticHeapMedian(const t_buffer_type& points, const t_container_type& indices, uint04 start, uint04 step, uint01 dimension)
    {
        const uint04 half_heap = t_bucket_size / 2;
        const uint04 max_heap = t_bucket_size / 2 + 1;
        uint04 left[t_bucket_size / 2 + 1];
        uint04 right[t_bucket_size / 2 + 1];
        sint04 nLeft = 1;
        sint04 nRight = 1;
        sint04 index = start;
        // pick first value as median candidate
        sint04 median = index;
        index += step;
 
        for (;;)
        {
            // get next value
            sint04 val = index;
            index += step;
            // if value is smaller than median, append to left heap
            if (points[indices[val]][dimension] < points[indices[median]][dimension])
            {
                // move biggest value to the heap top
                sint04 child = nLeft++, parent = (child - 1) / 2;
                while (parent && points[indices[val]][dimension] > points[indices[left[parent]]][dimension])
                {
                    left[child] = left[parent];
                    child = parent;
                    parent = (parent - 1) / 2;
                }
                left[child] = val;
 
                // if left heap is full
                if (nLeft == max_heap)
                {
                    // for each remaining value
                    for (; nRight != max_heap; nRight++)
                    {
                        // get next value
                        val = index;
                        index += step;
                        // if value is to be inserted in the left heap
                        if (points[indices[val]][dimension] < points[indices[median]][dimension])
                        {
                            child = points[indices[left[2]]][dimension] > points[indices[left[1]]][dimension] ? 2 : 1;
                            if (points[indices[val]][dimension] >= points[indices[left[child]]][dimension])
                                median = val;
                            else
                            {
                                median = left[child];
                                parent = child;
                                child = parent * 2 + 1;
                                while (child < max_heap)
                                {
                                    if (child < half_heap
                                        && points[indices[left[child + 1]]][dimension] > points[indices[left[child]]][dimension])
                                        ++child;
                                    if (points[indices[val]][dimension] >= points[indices[left[child]]][dimension])
                                        break;
                                    left[parent] = left[child];
                                    parent = child;
                                    child = parent * 2 + 1;
                                }
                                left[parent] = val;
                            }
                        }
                    }
                    return median;
                }
            }
 
            // else append to right heap
            else
            {
                // move smallest value to the heap top
                sint04 child = nRight++, parent = (child - 1) / 2;
                while (parent && points[indices[val]][dimension] < points[indices[right[parent]]][dimension])
                {
                    right[child] = right[parent];
                    child = parent;
                    parent = (parent - 1) / 2;
                }
                right[child] = val;
 
                // if right heap is full
                if (nRight == max_heap)
                {
                    // for each remaining value
                    for (; nLeft != max_heap; nLeft++)
                    {
                        // get next value
                        val = index;
                        index += step;
                        // if value is to be inserted in the right heap
                        if (points[indices[val]][dimension] > points[indices[median]][dimension])
                        {
                            child = points[indices[right[2]]][dimension] < points[indices[right[1]]][dimension] ? 2 : 1;
                            if (points[indices[val]][dimension] <= points[indices[right[child]]][dimension])
                                median = val;
                            else
                            {
                                median = right[child];
                                parent = child;
                                child = parent * 2 + 1;
                                while (child < max_heap)
                                {
                                    if (child < half_heap
                                        && points[indices[right[child + 1]]][dimension] < points[indices[right[child]]][dimension])
                                        ++child;
                                    if (points[indices[val]][dimension] <= points[indices[right[child]]][dimension])
                                        break;
                                    right[parent] = right[child];
                                    parent = child;
                                    child = parent * 2 + 1;
                                }
                                right[parent] = val;
                            }
                        }
                    }
                    return median;
                }
            }
        }
    }
 
    /**--------------------------------------------------------------------------------------------------
    \brief Returns the aproximate median, or value that is close to the median, significantly faster than
        it would take to calculate the true median.
    **/
    template<class t_buffer_type>
    uint04 ApxMedian(const t_buffer_type& elements, const Buffer<uint04>& indices, uint04 start, uint04 end, uint01 dimension)
    {
        uint04 skip = (end - start) / 64;
        if (skip == 0)
        {
            if (end - start > 32) return StaticHeapMedian<31>(elements, indices, start, 1, dimension);
            if (end - start > 16) return StaticHeapMedian<15>(elements, indices, start, 1, dimension);
            if (end - start > 8)  return StaticHeapMedian<7>(elements, indices, start, 1, dimension);
            if (end - start > 4)  return StaticHeapMedian<3>(elements, indices, start, 1, dimension);
            else
                return indices[start];
        }
        uint04 number_value;
        if (skip > 64) number_value = 64;
        else if (skip > 32) number_value = 32;
        else if (skip > 16) number_value = 16;
        else if (skip > 8) number_value = 8;
        else return StaticHeapMedian<63>(elements, indices, start, skip, dimension);
        uint04 median_locations[64];
        uint04 median_values[64];
        for (uint04 i = 0; i < number_value; i++)
        {
            uint04 median_loc = StaticHeapMedian<63>(elements, indices, start + i, skip, dimension);
            median_locations[i] = (indices[median_loc]);
            median_values[i] = median_loc;
        }
        switch (number_value)
        {
        case 64: return median_values[StaticHeapMedian<63>(elements, median_locations, 0, 1, dimension)];
        case 32: return median_values[StaticHeapMedian<31>(elements, median_locations, 0, 1, dimension)];
        case 16: return median_values[StaticHeapMedian<15>(elements, median_locations, 0, 1, dimension)];
        case 8:  return median_values[StaticHeapMedian<7>(elements, median_locations, 0, 1, dimension)];
        default: lib_assert(false, "Unknown number_value in roughHeapMedian");
        }
        return Constant<uint04>::Invalid;
    }
 
    /**--------------------------------------------------------------------------------------------------
    \brief Sorts all incoming indices within the given range such that the resulting index matrix will
        have values less than the index first, and greater than the index second.
    **/
    template<class t_type, class t_node_type>
    uint04 SortAboutValue(uint04 value_index, const Buffer<t_node_type>& points, Buffer<uint04>& indices, uint04 start, uint04 end, uint01 dimension)
    {
        const t_type value = points[indices[value_index]][dimension];
 
        std::swap(indices[--end], indices[value_index]);
        uint04 store = start;
        for (uint04 idx = start; idx < end; idx++)
        {
            if (points[indices[idx]][dimension] <= value)
            {
                std::swap(indices[idx], indices[store]);
                store++;
            }
        }
        std::swap(indices[end], indices[store]);
        return store;
    }
 
    /**--------------------------------------------------------------------------------------------------
    \brief Returns the aproximate Nth element, or value that is close to the Nth Element, significantly faster than
        it would take to calculate the true Nth element of the sorted series
    **/
    template<class t_type, class t_node_type>
    uint04 ApxNthElement(const Buffer<t_node_type>& points, Buffer<uint04>& indices, uint04 start, uint04 end, uint01 dimension)
    {
        const uint04 median_loc = ApxMedian(points, indices, start, end, dimension);
        return SortAboutValue<t_type>(median_loc, points, indices, start, end, dimension);
    }
 
    /**--------------------------------------------------------------------------------------------------
    \brief Sorts all incoming indices within the given range such that the resulting index matrix will
        have values less than the index first, and greater than the index second.
    **/
    template<class t_type, class t_buffer_type, class t_bounds_type>
    uint04 SortAboutValue(uint04 value_index, const t_buffer_type& points, Buffer<uint04>& indices, uint04 start, uint04 end, uint01 dimension, t_bounds_type& bounds_left, t_bounds_type& bounds_right)
    {
        const t_type value = points[indices[value_index]].center()[dimension];
        bounds_right.addToBounds(points[indices[value_index]]);
 
        std::swap(indices[--end], indices[value_index]);
        uint04 store = start;
        for (uint04 idx = start; idx < end; idx++)
        {
            auto element = points[indices[idx]];
            if (element.center()[dimension] == value)
            {
                if ((store - start) <= (idx - start) / 2)
                {
                    bounds_left.addToBounds(element);
                    std::swap(indices[idx], indices[store]);
                    store++;
                }
                else
                {
                    bounds_right.addToBounds(element);
                }
            }
            else if (element.center()[dimension] < value)
            {
                bounds_left.addToBounds(element);
                std::swap(indices[idx], indices[store]);
                store++;
            }
            else
            {
                bounds_right.addToBounds(element);
            }
        }
        std::swap(indices[end], indices[store]);
        return store;
    }
};