Range.hpp

#ifndef __FIBEROPERATIONS_RANGE_HPP
#define __FIBEROPERATIONS_RANGE_HPP

#include <field/MemBase.hpp>
#include <grid/Representation.hpp>
#include "fiberopDllApi.h"

#ifdef  USE_OPENMP
#include <omp.h>
#endif

#include <algorithm>

namespace Fiber
{

struct  fiberop_API DataRangeBase : public Interface<DataRangeBase>, public MemCore::Ageable
{
        struct fiberop_API NoDataRange : std::exception
        {
                ~NoDataRange() throw (); 

                override const char*what() const throw();
        }; 

        DataRangeBase(); 

        DataRangeBase(const MemCore::Ageable&Age);

        virtual ~DataRangeBase() = 0; 

        virtual void expandRange(const RefPtr<DataRangeBase>&DRB) = 0;
};

template <class Type>
struct  fiberop_API  DataRange : DataRangeBase
{
        Type    Minimum, Maximum;

public: 
        DataRange()
        {} 

        const Type&Min() const
        {
                return Minimum;
        } 

        const Type&Max() const
        {
                return Maximum;
        } 

        template <class OtherType>
        DataRange(const DataRange<OtherType>&DR)
        : Minimum(DR.Min() )
        , Maximum(DR.Max() )
        {}

        template <class OtherType>
        void    expand(const DataRange<OtherType>&DR)
        {
                if (Minimum > DR.Min() ) Minimum = DR.Min(); 
                if (Maximum < DR.Max() ) Maximum = DR.Max(); 
                update( DR );
        }

        override void expandRange(const RefPtr<DataRangeBase>&DRB)
        {
                if (RefPtr<DataRange<Type> > DR = DRB)
                {
                        expand(*DR);
                }
        } 

        template <class CompareType>
        bool    contains(const CompareType&Value) const
        {
                if (Value<Minimum) return false;
                if (Value>Maximum) return false;
                return true;
        } 

static  RefPtr<DataRange> get(CreativeArrayBase&CAB)
        {
                if (RefPtr<DataRange> DR = interface_cast<DataRange>(CAB) )
                {
                        if ( !DR->isOlderThan( CAB ) )
                                return DR;
                }

                return NullPtr();
        }

static  RefPtr<DataRange> get(Field&F)
        {
                if (RefPtr<DataRange> DR = interface_cast<DataRange>(F) )
                {
                        if ( !DR->isOlderThan( F ) )
                                return DR;
                }

                return NullPtr();
        } 

static  RefPtr<DataRange> compute(CreativeArrayBase&CAB)
        {
        RefPtr<MemBase> ArrayBase = CAB.create();
        RefPtr<TypedArray<Type> > DataArray = ArrayBase;
                if (!DataArray)
                        return NullPtr();

        CreativeIterator<Type>*values = DataArray->creativeIterator();
                if (!values)
                        return NullPtr();

                if (!*values)
                {
                        return NullPtr();
                }
        const Iterator<Type>&Data = *values;

        index_t count = Data.count();
                if (count<1)
                        return NullPtr(); 

//              printf("MINMAX ELEMENTS: %ld\n", count);

        RefPtr<DataRange> DR = new DataRange();

        Type    &MinValue = DR->Minimum,
                &MaxValue = DR->Maximum;
                MinValue = MaxValue = Data[ 0 ]; 

#ifdef  USE_OPENMP
        // parallel version 

        const int max_threads_possible =  omp_get_max_threads(); 

        // stl container (write access) is not thread safe, so every thread gets its own element
        std::vector<Type> minList(max_threads_possible); 
        std::vector<Type> maxList(max_threads_possible); 

        std::vector<int> valueFound(max_threads_possible);

                for(int i=0; i<max_threads_possible; i++)
                {
                        valueFound[i] = 0;      
                }


                #pragma omp parallel
                {
                Type localMin = Data[0];
                Type localMax = Data[0]; 

                        //#pragma omp for shared( Data, minList, maxList, minFound, maxFound)
                        #pragma omp for
                        for(index_t i = 1; i<count ; i++)
                        {
                                // printf("Executing index i=%i in thread %i of %i\n", i, omp_get_thread_num(), omp_get_num_threads()); 
                        
                        const Type&value = Data[i]; 
                                if (value < localMin) localMin = value; 
                                if (value > localMax) localMax = value;
                        } 

                        const int my_thread_id = omp_get_thread_num(); 

                                minList[my_thread_id] = localMin; 
                                maxList[my_thread_id] = localMax; 

                                valueFound[my_thread_id] = 1;
                } 
                // parallel execution ends 

                // ... variables declared locally within a structured block or a routine 
                // that is invoked from within a parallel region are private by default. 

                        for( int i=0; i<max_threads_possible; i++)
                        {
                                if( valueFound[i] == 1 )
                                {       
                                        if( minList[i] < MinValue ) MinValue = minList[i]; 
                                        if( maxList[i] > MaxValue ) MaxValue = maxList[i];
                                } 
                        }          
#else
                // straight forward version
                for(index_t i=1; i<count ; i++)
                {
                const Type&value = Data[i];
                        if (value < MinValue) MinValue = value;
                        if (value > MaxValue) MaxValue = value;
                } 
#endif 
//              std::cout << " Min = " << MinValue << std::endl; 
//              std::cout << " Max = " << MaxValue << std::endl;

                CAB.addInterface( DR ); 
                DR->update( CAB ); 

                return DR;
        }

static  RefPtr<DataRange> retrieve(CreativeArrayBase&CAB)
        {
        RefPtr<DataRange> DR = get(CAB);
                if (DR) return DR;

                return compute(CAB);
        }
};


template <class Type>
struct  FieldRange : Field::Iterator
{
        RefPtr<DataRange<Type> > theRange;

        override bool apply(const RefPtr<FragmentID>&f, const RefPtr<CreativeArrayBase>&F)
        {
                if (!F)
                        return true; 

                if (RefPtr<DataRange<Type> > FragmentRange = DataRange<Type>::retrieve( *F) )
                {
                        if (!theRange)
                                theRange = new DataRange<Type>( *FragmentRange ); 
                        else
                                theRange->expand( *FragmentRange ); 
                }

                return true;
        }

        FieldRange(Field&F, const MemCore::WeakPtr<FragmentSelector>&FS = NullPtr() )
        {
                if (!FS)
                {
                        theRange = DataRange<Type>::get(F);
#if     0
                        if (!theRange)
                                puts("FieldRange: No datarange yet"); 
                        else
                        {
                                printf("FieldRange: Returning old datarange from %ld, Field is %ld (field %p)\n", theRange->time_value(), F.time_value(), &F );
                        }
#endif
                }
                if (!theRange)
                {
                        F.iterate( *this, FS );

                        if (theRange && !FS)
                        {
                                F.addInterface( theRange );
                                theRange->update( F );
//                              puts("Added range to field");
                        }
                }
        }
};

extern fiberop_API template struct DataRange<float>;
extern fiberop_API template struct DataRange<double>;

extern fiberop_API template struct FieldRange<float>;
extern fiberop_API template struct FieldRange<double>;

template <class Type>
inline fiberop_API
bool    Contains(CreativeArrayBase&CAB, const Type&Value, bool ComputeEventually)
{
        if (ComputeEventually)
        {
                if (RefPtr<DataRange<Type> > DR = DataRange<Type>::retrieve(CAB) )
                {
                        return DR->contains( Value );
                }
        }
        else
        {
                if (RefPtr<DataRange<Type> > DR = DataRange<Type>::get(CAB) )
                {
                        return DR->contains( Value );
                }
        }
        throw DataRangeBase::NoDataRange();
}

extern fiberop_API RefPtr<DataRangeBase> getRange(Field&);


} // namespace Fiber


#endif // __FIBEROPERATIONS_RANGE_HPP