//Purpose:
//  Show affect of storage order.
//
#define INDENT_OUTPUT
#ifdef INDENT_OUTPUT
#include <boost/utility/trace_scope.hpp>
#endif
#include <numeric>
#include <boost/multi_array.hpp>
#include <boost/array.hpp>

  template
  < typename T
  >
  std::ostream&
operator<<
  ( std::ostream& sout
  , std::vector<T> const& a_container
  )
{
   sout<<"{ ";
   unsigned const size=a_container.size();
   for(unsigned i=0; i<size; ++i)
   {
      if(0<i) sout<<", ";
      sout<<a_container[i];
   }
   sout<<"} ";
   return sout;
}

  template
  < typename Permut
  , typename Strides
  >
  bool
mono_tonic_print
  ( Permut const& a_permut
  , Strides const& a_strides
  )
{
    std::cout<<":monotonically increasing={ ";
    unsigned const n=a_permut.size();
    bool ok=true;
    for(unsigned i=0; i<n-1; ++i)
    {
         if(0<i)std::cout<<", ";
         bool permut_increase=a_permut[i]<a_permut[i+1];
         bool stride_increase=a_strides[i]<a_strides[i+1];
         bool monotonic=permut_increase==stride_increase;
         std::cout<<monotonic;
         ok=ok && monotonic;
    }
    std::cout<<"}\n";
    return ok;
}
   
  template
  < typename Permut
  >
  Permut
invert_permut
  ( Permut const& a_permut
  )
{
    unsigned const n=a_permut.size();
    Permut result(n);
    for( unsigned i=0; i<n; ++i) result[a_permut[i]]=i;
    return result;
}

  template
  < std::size_t NumDims
  >
  std::vector<unsigned>
extents2lengths
  ( boost::detail::multi_array::extent_gen<NumDims> const& a_extents
  )
{
    unsigned const n=a_extents.ranges_.size();
    std::vector<unsigned> result(n);
    for( unsigned i=0; i<n; ++i) 
      result[i]=a_extents.ranges_[i].size();
    return result;
}    
  template
  < typename Permut
  >
  std::vector<unsigned>
expected_strides
  ( Permut const& a_permut
  , std::vector<unsigned> const& a_lengths
  )
{
    unsigned const n=a_permut.size();
    std::vector<unsigned> strides(n);
  #define NO_INVERT_PERMUT
  #ifdef NO_INVERT_PERMUT
    Permut permut_i=a_permut;
  #else
    Permut permut_i=invert_permut(a_permut);
  #endif
    strides[permut_i[0]]=1;
    for( unsigned i=1; i<n; ++i) 
    {
        unsigned i_prev=permut_i[i-1];
        strides[permut_i[i]]=strides[i_prev]*a_lengths[i_prev];
    }
    return strides;
}

unsigned const rank=3;
typedef boost::multi_array<int,rank> array_t;

void strides_print( array_t const& a_array)
{
    std::cout<<"strides()={ ";
    for(unsigned i=0; i<rank; ++i)
    {
         if(0<i)std::cout<<", ";
         std::cout<<a_array.strides()[i];
    }
    std::cout<<"}\n";
 }

typedef std::vector<unsigned> axes_permutation_t;

  axes_permutation_t
mk_permut( unsigned rank)
{
    axes_permutation_t result(rank);
    std::iota( result.begin(), result.end(), 0);
    return result;
}

  template
  < std::size_t NumDims
  >
  bool
operator==
  ( axes_permutation_t const& a_permut
  , boost::general_storage_order<NumDims>const& a_so
  )
{
    unsigned const n=a_permut.size();
    bool result=n==NumDims;
    for( unsigned i=0; result && i<n; ++i) result=a_permut[i]==a_so.ordering(i);
    return result;
}    

int main()
{
  #ifdef INDENT_OUTPUT
    boost::iostreams::indent_scoped_ostreambuf<char> indent_outbuf(std::cout,2);
  #endif
    std::vector<bool> const directions(rank,true);
    auto const extents=array_t::extent_gen()[2][3][4];
    auto const lengths=extents2lengths(extents);
    std::cout<<":lengths="<<lengths<<"\n";
    {
        std::string scope_name="fortran_order";
      #ifdef INDENT_OUTPUT
        ::boost::trace_scope ts(scope_name);
      #else
        std::cout<<"::"<<scope_name<<"::\n";
      #endif
        array_t ma( extents, boost::fortran_storage_order());
        strides_print(ma);
    }
    {
        std::string scope_name="c_order";
      #ifdef INDENT_OUTPUT
        ::boost::trace_scope ts(scope_name);
      #else
        std::cout<<"::"<<scope_name<<"::\n";
      #endif
        array_t ma( extents, boost::c_storage_order());
        strides_print(ma);
    }
    {
        std::string scope_name="all_permutations_order";
      #ifdef INDENT_OUTPUT
        ::boost::trace_scope ts(scope_name);
      #else
        std::cout<<"::"<<scope_name<<"::\n";
      #endif
        axes_permutation_t permut_v=mk_permut(rank);
        {
            unsigned count=0;
            do
            {
                std::cout<<":count="<<count<<"\n";
              #ifdef INDENT_OUTPUT
                ::boost::indent_scope<> is;
              #endif
                std::cout<<":permut_v="<<permut_v<<"\n";
                axes_permutation_t permut_i=invert_permut(permut_v);
                std::cout<<":permut_i="<<permut_i<<"\n";
                boost::general_storage_order<rank>
              store_order
                ( permut_v.begin()
                , directions.begin()
                );
                array_t
              ma
                ( extents
                , store_order
                );
                std::cout<<":ma.ordering(*)==permut_v="<<(permut_v==ma.storage_order())<<"\n";
                std::cout<<":multi_array_strides:";
                {
                  #ifdef INDENT_OUTPUT
                    ::boost::indent_scope<> is;
                  #endif
                    strides_print(ma);
                    bool ok=mono_tonic_print( permut_v, ma.strides());
                    std::cout<<":is_mono_tonic="<<ok<<"\n";
                }
                std::cout<<":expected_strides:";
                {
                  #ifdef INDENT_OUTPUT
                    ::boost::indent_scope<> is;
                  #endif
                    auto const strides_e=expected_strides( permut_v, lengths);
                    std::cout<<strides_e<<"\n";
                    bool ok=mono_tonic_print( permut_v, strides_e);
                    std::cout<<":is_mono_tonic="<<ok<<"\n";
                }
            }while 
              (  ++count< 10
              && std::next_permutation
                 ( permut_v.begin()
                 , permut_v.end()
                 )
              );
        }
    }
    return 0;
}