CoordinateAxisArray.hpp

#ifndef __FIBER_COORDINATEAXISARRAY_HPP
#define __FIBER_COORDINATEAXISARRAY_HPP
 
#ifndef DONT_USE_PRAGMA_ONCE
#pragma once
#endif
 
#include "MemArray.hpp"
#include <map>
 
namespace Fiber
{
 
struct  LinearCoordinateAxisArray {};
struct  PolynomialCoordinateAxisArray {};
struct  ExplicitCoordinateAxisArray {};
 
template <class T, class Category>
        class   CoordinateAxisArray;
 
template <class T>
class   CoordinateAxisArray<T, PolynomialCoordinateAxisArray>  : public MemCore::Chunk<T>
{
public:
        index_t length;
 
        CoordinateAxisArray(index_t max_idx)
        : MemCore::Chunk<T>(max_idx)
        {}
 
        T operator()(index_t i) const
        {
        index_t p = this->elements()-1;
        T       t = (*this)[p];
                while(p>0)
                {
                        t *= i;
                        p--;
                        t += (*this)[p];
                }
                return t;
        }
 
        index_t Length() const
        {
                return length;
        }
 
#if 0
        bool locate(T&value) const
        {
                if (Length()==2)
                {
                        value = (value - (*this)[0] ) / (*this)[1];
                        return true;
                }
 
                return false;
        }
#endif
};
 
 
template <class T>
class   CoordinateAxisArray<T, LinearCoordinateAxisArray>  : public MemCore::Chunk<T>
{
public:
        index_t length;
 
        CoordinateAxisArray() = delete;
 
        CoordinateAxisArray(index_t max_idx)
        : MemCore::Chunk<T>(2)
        , length(max_idx)
        {}
 
        CoordinateAxisArray(CoordinateAxisArray&Src)
        : MemCore::Chunk<T>( Src.std_vector() )
        , length( Src.length )
        {}
 
        CoordinateAxisArray(const T&C0, const T&C1, index_t max_idx)
        : MemCore::Chunk<T>(2)
        , length(max_idx)
        {
                (*this)[0] = C0;
                (*this)[1] = C1;
        }
 
        T operator()(index_t i) const
        {
                return (*this)[0] + i*(*this)[1];
        }
 
        double  reverse(const T&value) const
        {
        const T&delta = (*this)[1];
 
//              printf("CoordinateAxisArray::reverse(value=%lg) delta=%lg\n", double(value), double(delta) );
                // if there is no delta, then always yield the 0th coordinate index
                if ( delta == 0.0)
                {
//                      printf("CoordinateAxisArray::reverse() delta=%lg ZERO delta! \n", double(delta) );
                        return 0.0;
                }
                return (value - (*this)[0] ) / delta;
        }
 
        index_t Length() const
        {
                return length;
        }
 
 
#if     0
        bool    locate(T&value) const
        {
        double  FloatIndex = reverse(value);
 
                return FloatIndex>=0.0 && FloatIndex<=1.0;
        }
#endif
};
 
 
template <class T>
class   CoordinateAxisArray<T, ExplicitCoordinateAxisArray>  : public MemCore::Chunk<T>
{
public:
        typedef std::map<T, index_t> CoordinateLookupMap_t;
        CoordinateLookupMap_t   CoordinateLookupMap;
 
        void    setupCoordinateLookupMap()
        {
        const std::vector<T>&data = this->std_vector();
 
                CoordinateLookupMap.clear();
                for(index_t i=0; i<data.size(); i++)
                {
                        CoordinateLookupMap[ data[i] ] = i;
                }
        }
 
 
 
        CoordinateAxisArray(index_t max_idx)
        : MemCore::Chunk<T>(max_idx)
        {}
 
        CoordinateAxisArray(const CoordinateAxisArray&Src)
        : MemCore::Chunk<T>( Src.std_vector() )
        {}
 
        T operator()(index_t i) const
        {
                return (*this)[i];
        }
 
        index_t Length() const
        {
                return this->elements();
        }
 
        double  reverse(const T&value)
        {
                if (this->elements()<1) return 0.0;
#if 0
                MEMCORE_PROFILE_THIS("CoordinateAxisArray::reverse(): Non STL Bin Search", 1);
                const std::vector<T>&data = this->std_vector();
 
 
                size_t elements = CoordinateLookupMap.size() - 1;
                size_t low,mid,high;
                /*
                if(value < 0  ||  value > CoordinateLookupMap[elements])
                {
                        //printf("x[82]=%17.14f     px = %17.14f is outside the box\n",x_values[elements],px);
                        //scanf("%d",&n);
                        return false;
                }
                */
                low=0;
                high=elements;
                printf("\nReverse lookup, Explicit Coordinates: Elements= %d, Low= %d, high= %d\n", elements, low, high);
                if(value<CoordinateLookupMap[low])
                {
                        return CoordinateLookupMap[low];
                }
                else if(value>CoordinateLookupMap[high])
                {
                        return CoordinateLookupMap[high];
                }
                printf("\nCoordinate in range");
                while(high>low+1)
                {
                        mid=(low+high)/2;
                        if(value<CoordinateLookupMap[mid])
                        high=mid;
                        else
                        low=mid;
                }
                printf("\nReverse lookup, Explicit Coordinates: Elements= %d, Low= %d, high= %d\n", elements, low, high);
                T       LeftValue  = CoordinateLookupMap[low];
                T       RightValue = CoordinateLookupMap[high];
                index_t IndexValue = low;
 
#else
                if (CoordinateLookupMap.empty() )
                {
                        setupCoordinateLookupMap();
                }
        typename CoordinateLookupMap_t::const_iterator
                LocationJustLargerThanValue = CoordinateLookupMap.lower_bound( value );
 
                if (LocationJustLargerThanValue == CoordinateLookupMap.end() )
                        return double( CoordinateLookupMap.rbegin()->second ); // nothing found, max boundary
 
                if (LocationJustLargerThanValue == CoordinateLookupMap.begin() )
                        return double( CoordinateLookupMap.begin()->second ); // nothing found, min boundary
 
        typename CoordinateLookupMap_t::const_iterator
                LocationJustSmallerThanValue = LocationJustLargerThanValue;
                LocationJustSmallerThanValue--;
 
        T       LeftValue  = LocationJustSmallerThanValue->first;
        T       RightValue = LocationJustLargerThanValue ->first;
 
                if(value < LeftValue)
                {
                        printf("\n%lg, %lg, %lg\n", value, LeftValue, (value - LeftValue));
                        printf("%lg, %d \n", CoordinateLookupMap.begin()->first, int(CoordinateLookupMap.begin()->second) );
 
                }
                //assert( !(LeftValue < key) );
                assert( value >= LeftValue );
                assert( value <= RightValue );
 
        index_t IndexValue = LocationJustSmallerThanValue->second;
#endif
        double  Weight = (value - LeftValue) / ( RightValue - LeftValue);
 
                return IndexValue +  Weight;
 
        }
 
#if 0
        bool locate(T&value) const
        {
                if (Length()==2)
                {
                        value = (value - (*this)[0] ) / (*this)[1];
                        return true;
                }
 
                return false;
        }
#endif
};
 
 
} /* namespace Fiber */
 
#endif /* __FIBER_COORDINATEAXISARRAY_HPP */