/*
 * $Id: methods.cc,v 1.2 2002/09/12 07:02:06 eric Exp $
 *
 * $Log: methods.cc,v $
 * Revision 1.2  2002/09/12 07:02:06  eric
 * major rewrite of weave.
 *
 * 0.
 * The underlying library code is significantly re-factored and simpler. There used to be a xxx_spec.py and xxx_info.py file for every group of type conversion classes.  The spec file held the python code that handled the conversion and the info file had most of the C code templates that were generated.  This proved pretty confusing in practice, so the two files have mostly been merged into the spec file.
 *
 * Also, there was quite a bit of code duplication running around.  The re-factoring was able to trim the standard conversion code base (excluding blitz and accelerate stuff) by about 40%.  This should be a huge maintainability and extensibility win.
 *
 * 1.
 * With multiple months of using Numeric arrays, I've found some of weave's "magic variable" names unwieldy and want to change them.  The following are the old declarations for an array x of Float32 type:
 *
 *         PyArrayObject* x = convert_to_numpy(...);
 *         float* x_data = (float*) x->data;
 *         int*   _Nx = x->dimensions;
 *         int*   _Sx = x->strides;
 *         int    _Dx = x->nd;
 *
 * The new declaration looks like this:
 *
 *         PyArrayObject* x_array = convert_to_numpy(...);
 *         float* x = (float*) x->data;
 *         int*   Nx = x->dimensions;
 *         int*   Sx = x->strides;
 *         int    Dx = x->nd;
 *
 * This is obviously not backward compatible, and will break some code (including a lot of mine).  It also makes inline() code more readable and natural to write.
 *
 * 2.
 * I've switched from CXX to Gordon McMillan's SCXX for list, tuples, and dictionaries.  I like CXX pretty well, but its use of advanced C++ (templates, etc.) caused some portability problems.  The SCXX library is similar to CXX but doesn't use templates at all.  This, like (1) is not an
 * API compatible change and requires repairing existing code.
 *
 * I have also thought about boost python, but it also makes heavy use of templates.  Moving to SCXX gets rid of almost all template usage for the standard type converters which should help portability.  std::complex and std::string from the STL are the only templates left.  Of course blitz still uses templates in a major way so weave.blitz will continue to be hard on compilers.
 *
 * I've actually considered scrapping the C++ classes for list, tuples, and
 * dictionaries, and just fall back to the standard Python C API because the classes are waaay slower than the raw API in many cases.  They are also more convenient and less error prone in many cases, so I've decided to stick with them.  The PyObject variable will always be made available for variable "x" under the name "py_x" for more speedy operations.  You'll definitely want to use these for anything that needs to be speedy.
 *
 * 3.
 * strings are converted to std::string now.  I found this to be the most useful type in for strings in my code.  Py::String was used previously.
 *
 * 4.
 * There are a number of reference count "errors" in some of the less tested conversion codes such as instance, module, etc.  I've cleaned most of these up.  I put errors in quotes here because I'm actually not positive that objects passed into "inline" really need reference counting applied to them.  The dictionaries passed in by inline() hold references to these objects so it doesn't seem that they could ever be garbage collected inadvertently.  Variables used by ext_tools, though, definitely need the reference counting done.  I don't think this is a major cost in speed, so it probably isn't worth getting rid of the ref count code.
 *
 * 5.
 * Unicode objects are now supported.  This was necessary to support rendering Unicode strings in the freetype wrappers for Chaco.
 *
 * 6.
 * blitz++ was upgraded to the latest CVS.  It compiles about twice as fast as the old blitz and looks like it supports a large number of compilers (though only gcc 2.95.3 is tested).  Compile times now take about 9 seconds on my 850 MHz PIII laptop.
 *
 * Revision 1.7  2002/08/30 22:14:23  jcumming
 * Added definition of setStorage() method, which lets user set Array storage
 * format after construction.  We check that Array is not allocated first.
 *
 * Revision 1.6  2002/06/27 00:15:35  jcumming
 * Changed T_numtype to P_numtype when used outside the argument list or body
 * of a member function definition (i.e., outside the class scope).  Inside
 * the class scope, we can use the typedef T_numtype.  The IBM xlC compiler
 * gets confused if P_numtype is used as a template parameter name in a member
 * function declaration and then T_numtype is used as the parameter name in
 * the member function definition.  Fixed usage to be more consistent.
 *
 * Revision 1.5  2002/05/27 19:47:22  jcumming
 * Removed use of this->.  Members of templated base class are now declared
 * in derived class Array.
 *
 * Revision 1.4  2002/03/06 16:16:18  patricg
 *
 * data_ replaced by this->data_ everywhere
 * numReferences() replaced by this->numReferences()
 *
 * Revision 1.3  2001/01/26 18:30:50  tveldhui
 * More source code reorganization to reduce compile times.
 *
 * Revision 1.2  2001/01/24 23:41:53  tveldhui
 * Widespread changes to reduce compile time.  For backwards
 * compatibility, #include <blitz/array.h> enables BZ_GANG_INCLUDE
 * mode which includes all array and vector functionality (about
 * 120000 lines of code).  #include <blitz/Array.h> includes
 * a minimal subset of Array funcitonality; other features must
 * be included explicitly.
 *
 * Revision 1.1.1.1  2000/06/19 12:26:13  tveldhui
 * Imported sources
 *
 * Revision 1.4  1998/03/14 00:04:47  tveldhui
 * 0.2-alpha-05
 *
 * Revision 1.3  1997/08/18 19:13:08  tveldhui
 * Just prior to implementing fastRead() optimization for array
 * expression evaluation.
 *
 * Revision 1.2  1997/08/15 21:14:10  tveldhui
 * Just prior to loop-collapse change
 *
 */

#ifndef BZ_ARRAYMETHODS_CC
#define BZ_ARRAYMETHODS_CC

#ifndef BZ_ARRAY_H
 #error <blitz/array/methods.cc> must be included via <blitz/array.h>
#endif

BZ_NAMESPACE(blitz)

template<class P_numtype, int N_rank> template<class T_expr>
Array<P_numtype,N_rank>::Array(_bz_ArrayExpr<T_expr> expr)
{
    // Determine extent of the array expression

    TinyVector<int,N_rank> lbound, extent, ordering;
    TinyVector<bool,N_rank> ascendingFlag;

    for (int i=0; i < N_rank; ++i)
    {
        lbound(i) = expr.lbound(i);
        int ubound = expr.ubound(i);
        extent(i) = ubound - lbound(i) + 1;
        ordering(i) = expr.ordering(i);
        int ascending = expr.ascending(i);
        ascendingFlag(i) = (ascending == 1);

#ifdef BZ_DEBUG
        if ((lbound(i) == INT_MIN) || (ubound == INT_MAX) 
          || (ordering(i) == INT_MAX) || (ascending == INT_MAX))
        {
          BZPRECHECK(0,
           "Attempted to construct an array from an expression " << endl
           << "which does not have a shape.  To use this constructor, "
           << endl 
           << "the expression must contain at least one array operand.");
          return;
        }
#endif
    }

    Array<T_numtype,N_rank> A(lbound,extent,
        GeneralArrayStorage<N_rank>(ordering,ascendingFlag));
    A = expr;
    reference(A);
}

template<class P_numtype, int N_rank>
Array<P_numtype,N_rank>::Array(const TinyVector<int, N_rank>& lbounds,
    const TinyVector<int, N_rank>& extent,
    const GeneralArrayStorage<N_rank>& storage)
    : storage_(storage)
{
    length_ = extent;
    storage_.setBase(lbounds);
    setupStorage(N_rank - 1);
}


/*
 * This routine takes the storage information for the array
 * (ascendingFlag_[], base_[], and ordering_[]) and the size
 * of the array (length_[]) and computes the stride vector
 * (stride_[]) and the zero offset (see explanation in array.h).
 */
template<class P_numtype, int N_rank>
_bz_inline2 void Array<P_numtype, N_rank>::computeStrides()
{
    if (N_rank > 1)
    {
      int stride = 1;

      // This flag simplifies the code in the loop, encouraging
      // compile-time computation of strides through constant folding.
      _bz_bool allAscending = storage_.allRanksStoredAscending();

      // BZ_OLD_FOR_SCOPING
      int n;
      for (n=0; n < N_rank; ++n)
      {
          int strideSign = +1;

          // If this rank is stored in descending order, then the stride
          // will be negative.
          if (!allAscending)
          {
            if (!isRankStoredAscending(ordering(n)))
                strideSign = -1;
          }

          // The stride for this rank is the product of the lengths of
          // the ranks minor to it.
          stride_[ordering(n)] = stride * strideSign;

          stride *= length_[ordering(n)];
      }
    }
    else {
        // Specialization for N_rank == 1
        // This simpler calculation makes it easier for the compiler
        // to propagate stride values.

        if (isRankStoredAscending(0))
            stride_[0] = 1;
        else
            stride_[0] = -1;
    }

    calculateZeroOffset();
}

template<class P_numtype, int N_rank>
void Array<P_numtype, N_rank>::calculateZeroOffset()
{
    // Calculate the offset of (0,0,...,0)
    zeroOffset_ = 0;

    // zeroOffset_ = - sum(where(ascendingFlag_, stride_ * base_,
    //     (length_ - 1 + base_) * stride_))
    for (int n=0; n < N_rank; ++n)
    {
        if (!isRankStoredAscending(n))
            zeroOffset_ -= (length_[n] - 1 + base(n)) * stride_[n];
        else
            zeroOffset_ -= stride_[n] * base(n);
    }
}

template<class P_numtype, int N_rank>
_bz_bool Array<P_numtype, N_rank>::isStorageContiguous() const
{
    // The storage is contiguous if for the set
    // { | stride[i] * extent[i] | }, i = 0..N_rank-1,
    // there is only one value which is not in the set
    // of strides; and if there is one stride which is 1.

    // This algorithm is quadratic in the rank.  It is hard
    // to imagine this being a serious problem.

    int numStridesMissing = 0;
    bool haveUnitStride = _bz_false;

    for (int i=0; i < N_rank; ++i)
    {
        int stride = BZ_MATHFN_SCOPE(abs)(stride_[i]);
        if (stride == 1)
            haveUnitStride = _bz_true;

        int vi = stride * length_[i];

        int j = 0;
        for (j=0; j < N_rank; ++j)
            if (BZ_MATHFN_SCOPE(abs)(stride_[j]) == vi)
                break;

        if (j == N_rank)
        {
            ++numStridesMissing;
            if (numStridesMissing == 2)
                return _bz_false;
        }
    }

    return haveUnitStride;
}

template<class P_numtype, int N_rank>
void Array<P_numtype, N_rank>::dumpStructureInformation(ostream& os) const
{
    os << "Dump of Array<" << BZ_DEBUG_TEMPLATE_AS_STRING_LITERAL(P_numtype) 
       << ", " << N_rank << ">:" << endl
       << "ordering_      = " << storage_.ordering() << endl
       << "ascendingFlag_ = " << storage_.ascendingFlag() << endl
       << "base_          = " << storage_.base() << endl
       << "length_        = " << length_ << endl
       << "stride_        = " << stride_ << endl
       << "zeroOffset_    = " << zeroOffset_ << endl
       << "numElements()  = " << numElements() << endl
       << "isStorageContiguous() = " << isStorageContiguous() << endl;
}

/*
 * Make this array a view of another array's data.
 */
template<class P_numtype, int N_rank>
void Array<P_numtype, N_rank>::reference(const Array<P_numtype, N_rank>& array)
{
    storage_ = array.storage_;
    length_ = array.length_;
    stride_ = array.stride_;
    zeroOffset_ = array.zeroOffset_;

    MemoryBlockReference<P_numtype>::changeBlock(array.noConst(),
        array.zeroOffset_);

    data_ = const_cast<P_numtype*>(array.data_);
}

/*
 * Modify the Array storage.  Array must be unallocated.
 */
template<class P_numtype, int N_rank>
void Array<P_numtype, N_rank>::setStorage(GeneralArrayStorage<N_rank> x)
{
#ifdef BZ_DEBUG
    if (size() != 0) {
        BZPRECHECK(0,"Cannot modify storage format of an Array that has already been allocated!" << endl);
        return;
    }
#endif
    storage_ = x;
    return;
}

/*
 * This method is called to allocate memory for a new array.  
 */
template<class P_numtype, int N_rank>
_bz_inline2 void Array<P_numtype, N_rank>::setupStorage(int lastRankInitialized)
{
    TAU_TYPE_STRING(p1, "Array<T,N>::setupStorage() [T="
        + CT(P_numtype) + ",N=" + CT(N_rank) + "]");
    TAU_PROFILE(" ", p1, TAU_BLITZ);

    /*
     * If the length of some of the ranks was unspecified, fill these
     * in using the last specified value.
     *
     * e.g. Array<int,3> A(40) results in a 40x40x40 array.
     */
    for (int i=lastRankInitialized + 1; i < N_rank; ++i)
    {
        storage_.setBase(i, storage_.base(lastRankInitialized));
        length_[i] = length_[lastRankInitialized];
    }

    // Compute strides
    computeStrides();

    // Allocate a block of memory
    MemoryBlockReference<P_numtype>::newBlock(numElements());

    // Adjust the base of the array to account for non-zero base
    // indices and reversals
    data_ += zeroOffset_;
}

template<class P_numtype, int N_rank>
Array<P_numtype, N_rank> Array<P_numtype, N_rank>::copy() const
{
    if (numElements())
    {
        Array<T_numtype, N_rank> z(length_, storage_);
        z = *this;
        return z;
    }
    else {
        // Null array-- don't bother allocating an empty block.
        return *this;
    }
}

template<class P_numtype, int N_rank>
void Array<P_numtype, N_rank>::makeUnique()
{
    if (numReferences() > 1)
    {
        T_array tmp = copy();
        reference(tmp);
    }
}

template<class P_numtype, int N_rank>
Array<P_numtype, N_rank> Array<P_numtype, N_rank>::transpose(int r0, int r1, 
    int r2, int r3, int r4, int r5, int r6, int r7, int r8, int r9, int r10)
{
    T_array B(*this);
    B.transposeSelf(r0,r1,r2,r3,r4,r5,r6,r7,r8,r9,r10);
    return B;
}

template<class P_numtype, int N_rank>
void Array<P_numtype, N_rank>::transposeSelf(int r0, int r1, int r2, int r3,
    int r4, int r5, int r6, int r7, int r8, int r9, int r10)
{
    BZPRECHECK(r0+r1+r2+r3+r4+r5+r6+r7+r8+r9+r10 == N_rank * (N_rank-1) / 2,
        "Invalid array transpose() arguments." << endl
        << "Arguments must be a permutation of the numerals (0,...,"
        << (N_rank - 1) << ")");

    // Create a temporary reference copy of this array
    Array<T_numtype, N_rank> x(*this);

    // Now reorder the dimensions using the supplied permutation
    doTranspose(0, r0, x);
    doTranspose(1, r1, x);
    doTranspose(2, r2, x);
    doTranspose(3, r3, x);
    doTranspose(4, r4, x);
    doTranspose(5, r5, x);
    doTranspose(6, r6, x);
    doTranspose(7, r7, x);
    doTranspose(8, r8, x);
    doTranspose(9, r9, x);
    doTranspose(10, r10, x);
}

template<class P_numtype, int N_rank>
void Array<P_numtype, N_rank>::doTranspose(int destRank, int sourceRank,
    Array<T_numtype, N_rank>& array)
{
    // BZ_NEEDS_WORK: precondition check

    if (destRank >= N_rank)
        return;

    length_[destRank] = array.length_[sourceRank];
    stride_[destRank] = array.stride_[sourceRank];
    storage_.setAscendingFlag(destRank, 
        array.isRankStoredAscending(sourceRank));
    storage_.setBase(destRank, array.base(sourceRank));

    // BZ_NEEDS_WORK: Handling the storage ordering is currently O(N^2)
    // but it can be done fairly easily in linear time by constructing
    // the appropriate permutation.

    // Find sourceRank in array.storage_.ordering_
    int i=0;
    for (; i < N_rank; ++i)
        if (array.storage_.ordering(i) == sourceRank)
            break;

    storage_.setOrdering(i, destRank);
}

template<class P_numtype, int N_rank>
void Array<P_numtype, N_rank>::reverseSelf(int rank)
{
    BZPRECONDITION(rank < N_rank);

    storage_.setAscendingFlag(rank, !isRankStoredAscending(rank));

    int adjustment = stride_[rank] * (lbound(rank) + ubound(rank));
    zeroOffset_ += adjustment;
    data_ += adjustment;
    stride_[rank] *= -1;
}

template<class P_numtype, int N_rank>
Array<P_numtype, N_rank> Array<P_numtype,N_rank>::reverse(int rank)
{
    T_array B(*this);
    B.reverseSelf(rank);
    return B;
}

template<class P_numtype, int N_rank> template<class P_numtype2>
Array<P_numtype2,N_rank> Array<P_numtype,N_rank>::extractComponent(P_numtype2, 
    int componentNumber, int numComponents) const
{
    BZPRECONDITION((componentNumber >= 0) 
        && (componentNumber < numComponents));

    TinyVector<int,N_rank> stride2;
    for (int i=0; i < N_rank; ++i)
      stride2(i) = stride_(i) * numComponents;
    const P_numtype2* dataFirst2 = 
        ((const P_numtype2*)dataFirst()) + componentNumber;
    return Array<P_numtype2,N_rank>(const_cast<P_numtype2*>(dataFirst2), 
        length_, stride2, storage_);
}

/* 
 * These routines reindex the current array to use a new base vector.
 * The first reindexes the array, the second just returns a reindex view
 * of the current array, leaving the current array unmodified.
 * (Contributed by Derrick Bass)
 */
template<class P_numtype, int N_rank>
_bz_inline2 void Array<P_numtype, N_rank>::reindexSelf(const 
    TinyVector<int, N_rank>& newBase) 
{
    int delta = 0;
    for (int i=0; i < N_rank; ++i)
      delta += (base(i) - newBase(i)) * stride_(i);

    data_ += delta;

    // WAS: dot(base() - newBase, stride_);

    storage_.setBase(newBase);
    calculateZeroOffset();
}

template<class P_numtype, int N_rank>
_bz_inline2 Array<P_numtype, N_rank> 
Array<P_numtype, N_rank>::reindex(const TinyVector<int, N_rank>& newBase) 
{
    T_array B(*this);
    B.reindexSelf(newBase);
    return B;
}

BZ_NAMESPACE_END

#endif // BZ_ARRAY_CC

