Array.h

Go to the documentation of this file.

#ifndef SimTK_SimTKCOMMON_ARRAY_H_
#define SimTK_SimTKCOMMON_ARRAY_H_

/* -------------------------------------------------------------------------- *
 *                       Simbody(tm): SimTKcommon                             *
 * -------------------------------------------------------------------------- *
 * This is part of the SimTK biosimulation toolkit originating from           *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody.  *
 *                                                                            *
 * Portions copyright (c) 2010-13 Stanford University and the Authors.        *
 * Authors: Michael Sherman                                                   *
 * Contributors:                                                              *
 *                                                                            *
 * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
 * not use this file except in compliance with the License. You may obtain a  *
 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
 *                                                                            *
 * Unless required by applicable law or agreed to in writing, software        *
 * distributed under the License is distributed on an "AS IS" BASIS,          *
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
 * See the License for the specific language governing permissions and        *
 * limitations under the License.                                             *
 * -------------------------------------------------------------------------- */

#include "SimTKcommon/internal/common.h"
#include "SimTKcommon/internal/ExceptionMacros.h"
#include "SimTKcommon/internal/Serialize.h"

#include <algorithm>
#include <iterator>
#include <vector>
#include <ostream>
#include <climits>
#include <typeinfo>

namespace SimTK {

// These are the classes defined in this header.
template <class X>                   struct ArrayIndexTraits;
template <class T, class X=unsigned> class  ArrayViewConst_;
template <class T, class X=unsigned> class  ArrayView_;
template <class T, class X=unsigned> class  Array_;



//==============================================================================
//                           CLASS ArrayIndexTraits
//==============================================================================

template <class X> struct ArrayIndexTraits {
    typedef typename X::size_type       size_type;
    typedef typename X::difference_type difference_type;
    static size_type max_size() {return X::max_size();}
};

template <> struct ArrayIndexTraits<unsigned> {
    typedef unsigned        size_type;
    typedef int             difference_type;
    static size_type        max_size() {return (unsigned)INT_MAX;}
};

template <> struct ArrayIndexTraits<int> {
    typedef int             size_type;
    typedef int             difference_type;
    static size_type        max_size() {return INT_MAX;}
};

template <> struct ArrayIndexTraits<unsigned long> {
    typedef unsigned long       size_type;
    typedef long                difference_type;
    static size_type            max_size() {return (unsigned long)LONG_MAX;}
};

template <> struct ArrayIndexTraits<long> {
    typedef long                size_type;
    typedef long                difference_type;
    static size_type            max_size() {return LONG_MAX;}
};

template <> struct ArrayIndexTraits<unsigned short> {
    typedef unsigned short      size_type;
    typedef int                 difference_type;
    static size_type            max_size() {return USHRT_MAX;}
};

template <> struct ArrayIndexTraits<short> {
    typedef short               size_type;
    typedef short               difference_type;
    static size_type            max_size() {return SHRT_MAX;}
}; 


template <> struct ArrayIndexTraits<unsigned char> {
    typedef unsigned char       size_type;
    typedef short               difference_type;
    static size_type            max_size() {return UCHAR_MAX;} // not CHAR_MAX
};

template <> struct ArrayIndexTraits<signed char> {
    typedef signed char         size_type;
    typedef signed char         difference_type;
    static size_type            max_size() {return SCHAR_MAX;}
};

template <> struct ArrayIndexTraits<char> {
    typedef char                size_type;
    typedef signed char         difference_type;
    static size_type            max_size() {return (char)SCHAR_MAX;}
};

template <> struct ArrayIndexTraits<bool> {
    typedef unsigned char       size_type;
    typedef signed char         difference_type;
    static size_type            max_size() {return 2;}
};

template <> struct ArrayIndexTraits<unsigned long long> {
    typedef unsigned long long  size_type;
    typedef long long           difference_type;
    static size_type            max_size() 
                                    {return (unsigned long long)LLONG_MAX;}
};

template <> struct ArrayIndexTraits<long long> {
    typedef long long           size_type;
    typedef long long           difference_type;
    static size_type            max_size() {return LLONG_MAX;}
};

// Don't show this in Doxygen.
// This helper class decides what integral type we should use to best pack
// the index type's size_type representation. The idea is to pack the whole
// Array_ structure into 8 bytes on a 32 bit machine, 16 bytes on a 64 bit
// machine, using the largest integral type that will work, giving a layout
// like this:          |       data pointer     |
//                     |   nUsed   | nAllocated |

// The default implementation just uses the integral type itself.
template <class Integral, class is64Bit> struct ArrayIndexPackTypeHelper 
{   typedef Integral packed_size_type;};

// On 32 bit machine, pack anything smaller than a short into a short.
template<> struct ArrayIndexPackTypeHelper<bool,FalseType> 
{   typedef unsigned short packed_size_type;};
template<> struct ArrayIndexPackTypeHelper<char,FalseType> 
{   typedef unsigned short packed_size_type;};
template<> struct ArrayIndexPackTypeHelper<unsigned char,FalseType> 
{   typedef unsigned short packed_size_type;};
template<> struct ArrayIndexPackTypeHelper<signed char,FalseType> 
{   typedef short packed_size_type;};

// On 64 bit machine, pack anything smaller than an int into an int.
template<> struct ArrayIndexPackTypeHelper<bool,TrueType> 
{   typedef unsigned int packed_size_type;};
template<> struct ArrayIndexPackTypeHelper<char,TrueType> 
{   typedef unsigned int packed_size_type;};
template<> struct ArrayIndexPackTypeHelper<unsigned char,TrueType> 
{   typedef unsigned int packed_size_type;};
template<> struct ArrayIndexPackTypeHelper<signed char,TrueType> 
{   typedef int packed_size_type;};
template<> struct ArrayIndexPackTypeHelper<unsigned short,TrueType> 
{   typedef unsigned int packed_size_type;};
template<> struct ArrayIndexPackTypeHelper<short,TrueType> 
{   typedef int packed_size_type;};

template <class Integral> struct ArrayIndexPackType
{   typedef typename ArrayIndexPackTypeHelper<Integral,Is64BitPlatformType>
                        ::packed_size_type  packed_size_type;};
//==============================================================================
//                            CLASS ArrayViewConst_
//==============================================================================
template <class T, class X> class ArrayViewConst_ {
public:


//------------------------------------------------------------------------------
typedef T           value_type;
typedef X           index_type;
typedef T*          pointer;
typedef const T*    const_pointer;
typedef T&          reference;
typedef const T&    const_reference;
typedef T*          iterator;
typedef const T*    const_iterator;
typedef std::reverse_iterator<iterator>                 reverse_iterator;
typedef std::reverse_iterator<const_iterator>           const_reverse_iterator;
typedef typename ArrayIndexTraits<X>::size_type         size_type;
typedef typename ArrayIndexTraits<X>::difference_type   difference_type;
typedef typename ArrayIndexPackType<size_type>::packed_size_type 
                                                        packed_size_type;
//------------------------------------------------------------------------------

ArrayViewConst_() : pData(0), nUsed(0), nAllocated(0) {}

ArrayViewConst_(const ArrayViewConst_& src) 
:   pData(0), nUsed(src.nUsed), nAllocated(0) {
    if (nUsed) pData = const_cast<T*>(src.pData);
} 

// Copy assignment is suppressed.

ArrayViewConst_(const T* first, const T* last1) 
:   pData(0),nUsed(0),nAllocated(0) { 
    if (last1==first) return; // empty

    SimTK_ERRCHK((first&&last1)||(first==last1), 
        "ArrayViewConst_<T>(first,last1)", 
        "One of the source pointers was null (0); either both must be"
        " non-null or both must be null.");

    SimTK_ERRCHK3(this->isSizeOK(last1-first), 
        "ArrayViewConst_<T>(first,last1)",
        "The source data's size %llu is too big for this array which"
        " is limited to %llu elements by its index type %s.",
        this->ull(last1-first), ullMaxSize(), indexName());

    pData = const_cast<T*>(first); 
    nUsed = packed_size_type(last1-first); 
    // nAllocated is already zero
}

template <class A>
ArrayViewConst_(const std::vector<T,A>& src) 
:   pData(0),nUsed(0),nAllocated(0) { 
    if (src.empty()) return;

    SimTK_ERRCHK3(this->isSizeOK(src.size()),
        "ArrayViewConst_<T>::ctor(std::vector<T>)",
        "The source std::vector's size %llu is too big for this array which"
        " is limited to %llu elements by its index type %s.",
        this->ull(src.size()), ullMaxSize(), indexName());

    pData = const_cast<T*>(&src.front()); 
    nUsed = packed_size_type(src.size()); 
    // nAllocated is already zero
}
operator const ArrayView_<T,X>&() const
{   return *reinterpret_cast<const ArrayView_<T,X>*>(this); }
operator const Array_<T,X>&() const
{   return *reinterpret_cast<const Array_<T,X>*>(this); }

void disconnect() {
    SimTK_ASSERT(nAllocated==0,
        "ArrayViewConst_::deallocate(): called on an owner Array_");
    nUsed = 0;
    pData = 0;
}

~ArrayViewConst_() {
    disconnect();
}
//------------------------------------------------------------------------------

size_type size() const {return size_type(nUsed);}
size_type max_size() const {return ArrayIndexTraits<X>::max_size();}
bool empty() const {return nUsed==0;}
size_type capacity() const {return size_type(nAllocated?nAllocated:nUsed);}
size_type allocated() const {return size_type(nAllocated);}
bool isOwner() const {return nAllocated || pData==0;}
/*}*/


//------------------------------------------------------------------------------

const T& operator[](index_type i) const {
    SimTK_INDEXCHECK(size_type(i),size(),"ArrayViewConst_<T>::operator[]()");
    return pData[i];
}
const T& at(index_type i) const {
    SimTK_INDEXCHECK_ALWAYS(size_type(i),size(),"ArrayViewConst_<T>::at()");
    return pData[i];
}
const T& getElt(index_type i) const {
    SimTK_INDEXCHECK(size_type(i),size(),"ArrayViewConst_<T>::getElt()");
    return pData[i];
}
const T& front() const 
{   SimTK_ERRCHK(!empty(), "ArrayViewConst_<T>::front()", "Array was empty.");
    return pData[0]; }
const T& back() const 
{   SimTK_ERRCHK(!empty(), "ArrayViewConst_<T>::back()", "Array was empty.");
    return pData[nUsed-1]; }

ArrayViewConst_ operator()(index_type index, size_type length) const {
    const size_type ix(index);
    SimTK_ERRCHK2(isSizeInRange(ix, size()), "ArrayViewConst_<T>(index,length)",
        "For this operator, we must have 0 <= index <= size(), but"
        " index==%llu and size==%llu.", this->ull(ix), ullSize());
    SimTK_ERRCHK2(isSizeInRange(length, size_type(size()-ix)), 
        "ArrayViewConst_<T>(index,length)", 
        "This operator requires 0 <= length <= size()-index, but"
        " length==%llu and size()-index==%llu.",this->ull(length),this->ull(size()-ix));

    return ArrayViewConst_(pData+ix, pData+ix+length);
}
ArrayViewConst_ getSubArray(index_type index, size_type length) const
{   return (*this)(index,length); }

//------------------------------------------------------------------------------

const T* cbegin() const {return pData;}
const T* cend() const {return pData + nUsed;}
const T* begin() const {return pData;}
const T* end() const {return pData + nUsed;}

const_reverse_iterator crbegin() const {return const_reverse_iterator(cend());}
const_reverse_iterator crend() const {return const_reverse_iterator(cbegin());}
const_reverse_iterator rbegin() const {return crbegin();} 
const_reverse_iterator rend() const {return crend();}

const T* cdata() const {return pData;}
const T* data() const {return pData;}
//------------------------------------------------------------------------------
                                 protected:
//------------------------------------------------------------------------------
// The remainder of this class is for the use of the ArrayView_<T,X> and
// Array_<T,X> derived classes only and should not be documented for users to 
// see. 
                                     
// Don't let doxygen see any of this.
packed_size_type psize() const {return nUsed;}
packed_size_type pallocated() const {return nAllocated;}

// These provide direct access to the data members for our trusted friends.
void setData(const T* p)        {pData = const_cast<T*>(p);}
void setSize(size_type n)       {nUsed = packed_size_type(n);}
void incrSize()                 {++nUsed;}
void decrSize()                 {--nUsed;}
void setAllocated(size_type n)  {nAllocated = packed_size_type(n);}

// Check whether a given size is the same as the current size of this array,
// avoiding any compiler warnings due to mismatched integral types.
template <class S> 
bool isSameSize(S sz) const
{   return ull(sz) == ullSize(); }

// Check that a source object's size will fit in the array being careful to
// avoid overflow and warnings in the comparison.
template <class S> 
bool isSizeOK(S srcSz) const
{   return ull(srcSz) <= ullMaxSize(); }

// This is identical in function to std::distance() (reports how many 
// elements lie between two iterators) but avoids any slow 
// Release-build bugcatchers that Microsoft may have felt compelled to add.
// The implementation is specialized for random access iterators because
// they can measure distance very fast.
template<class Iter> static
typename std::iterator_traits<Iter>::difference_type
iterDistance(const Iter& first, const Iter& last1) {
    return iterDistanceImpl(first,last1,
                typename std::iterator_traits<Iter>::iterator_category());
}

// Generic slow implementation for non-random access iterators. This is fine
// for forward and bidirectional iterators, but it will *consume* input
// iterators so is useless for them.
template<class Iter> static
typename std::iterator_traits<Iter>::difference_type
iterDistanceImpl(const Iter& first, const Iter& last1, std::input_iterator_tag) {
    typename std::iterator_traits<Iter>::difference_type d = 0;
    for (Iter src=first; src != last1; ++src, ++d)
        ;
    return d;
}

// Fast specialization for random access iterators (including ordinary
// pointers) -- just subtract.
template<class Iter> static
typename std::iterator_traits<Iter>::difference_type
iterDistanceImpl(const Iter& first, const Iter& last1, 
                 std::random_access_iterator_tag) {
    return last1 - first;
}

// This method attempts to determine whether any elements in the iterator range
// [first,last1) overlap with the elements stored in this array. This is used 
// for error checks for operations where source is not permitted to overlap the
// destination. For random access iterators (including ordinary pointers), we 
// can answer this question definitively because we expect the data to be 
// consecutive in memory. For other kinds of iterators, we will just assume
// there is no overlap. Note that null ranges do not overlap even if the
// pair of equal iterators points within the other range -- what matters is
// the number of overlapping elements.
template<class Iter> bool
overlapsWithData(const Iter& first, const Iter& last1) {
    return overlapsWithDataImpl(first,last1,
                typename std::iterator_traits<Iter>::iterator_category());
}

// This is a partial specialization of the above where the data is given
// with ordinary pointers.
template <class T2> bool
overlapsWithData(const T2* first, const T2* last1) {
    // Find the start and end+1 of the alleged overlap region. There is
    // overlap iff end+1 > start. Note that this works if either range 
    // is [0,0) or [p,p), or if last1 is illegally less than first (we just
    // want to report no overlap in that case -- it is someone else's business
    // to complain).
    const T* obegin = std::max(cbegin(), (const T*)first);
    const T* oend1  = std::min(cend(),   (const T*)last1);

    return obegin < oend1;
}

// This is the generic implementation for any type of input iterator other than
// random access (i.e., bidirectional, forward, or input) -- assume no overlap.
template<class Iter> bool
overlapsWithDataImpl(const Iter&, const Iter&, std::input_iterator_tag) 
{   return false; }

// Here we can actually test for overlap since we have random access iterators.
// We convert them to pointers and then look for memory overlap.
template<class Iter> bool
overlapsWithDataImpl(const Iter& first, const Iter& last1, 
                     std::random_access_iterator_tag) {
    // We must check that the input iterators span a non-zero range before
    // assuming we can dereference them.
    if (last1 <= first)
        return false; // zero or malformed source range: no overlap

    // We now know we can dereference first and last1-1 (can't safely 
    // dereference last1 but we can use pointer arithmetic to point past
    // the (last-1)th element in memory). We then take the dereferenced
    // object's address to get ordinary pointers that we can use to 
    // watch for illegal overlap.
    return overlapsWithData(&*first, &*(last1-1)); // use pointer overload
}

// Cast an integral type to maximal-width unsigned long long to avoid accidental
// overflows that might otherwise occur due to wraparound that can happen 
// with small index types.
template <class S>
static unsigned long long ull(S sz)
{   return (unsigned long long)sz; }

// Return size(), capacity(), and max_size() cast to unsigned long long.
unsigned long long ullSize()     const {return ull(size());}
unsigned long long ullCapacity() const {return ull(capacity());}
unsigned long long ullMaxSize()  const {return ull(max_size());}

// Useful in error messages for explaining why something was too big.
const char* indexName() const {return NiceTypeName<X>::name();}

private:
//------------------------------------------------------------------------------
//                               DATA MEMBERS
// These are the only data members and this layout is guaranteed not to change
// from release to release. If data is null, then nUsed==nAllocated==0.

T*                  pData;      // ptr to data referenced here, or 0 if none
packed_size_type    nUsed;      // number of elements currently present (size)
packed_size_type    nAllocated; // heap allocation; 0 if pData is not owned

ArrayViewConst_& operator=(const ArrayViewConst_& src); // suppressed
};






//==============================================================================
//                            CLASS ArrayView_
//==============================================================================
template <class T, class X> class ArrayView_ : public ArrayViewConst_<T,X> {
typedef ArrayViewConst_<T,X> CBase;
public:
//------------------------------------------------------------------------------
/*{*/
typedef T                                               value_type;
typedef X                                               index_type;
typedef T*                                              pointer;
typedef const T*                                        const_pointer;
typedef T&                                              reference;
typedef const T&                                        const_reference;
typedef T*                                              iterator;
typedef const T*                                        const_iterator;
typedef std::reverse_iterator<iterator>                 reverse_iterator;
typedef std::reverse_iterator<const_iterator>           const_reverse_iterator;
typedef typename ArrayIndexTraits<X>::size_type         size_type;
typedef typename ArrayIndexTraits<X>::difference_type   difference_type;
typedef typename ArrayIndexPackType<size_type>::packed_size_type 
                                                        packed_size_type;
/*}*/


//------------------------------------------------------------------------------

ArrayView_() : CBase() {}

ArrayView_(const ArrayView_& src) : CBase(src) {}

ArrayView_(T* first, const T* last1) : CBase(first,last1) {} 

template <class A>
ArrayView_(std::vector<T,A>& v) : CBase(v) {}

operator const Array_<T,X>&() const 
{   return *reinterpret_cast<const Array_<T,X>*>(this); }  

operator Array_<T,X>&() 
{   return *reinterpret_cast<Array_<T,X>*>(this); } 

void disconnect() {this->CBase::disconnect();}

~ArrayView_() {this->CBase::disconnect();}
//------------------------------------------------------------------------------

ArrayView_& operator=(const ArrayView_& src) {
    if (&src != this)
        avAssignIteratorDispatch(src.cbegin(), src.cend(),
                                 std::random_access_iterator_tag(),
                                 "ArrayView_<T>::operator=(ArrayView_<T>)");
    return *this;
}


template <class T2, class X2>
ArrayView_& operator=(const ArrayViewConst_<T2,X2>& src) {
    if ((const void*)&src != (void*)this)
        avAssignIteratorDispatch(src.cbegin(), src.cend(),
                                 std::random_access_iterator_tag(),
                                 "ArrayView_<T>::operator=(Array_<T2>)");
    return *this;
}

// Help out dumb compilers struggling to match the template arguments and
// promote the Array_ or ArrayView_ to ArrayConstView_ at the same time.

template <class T2, class X2>
ArrayView_& operator=(const ArrayView_<T2,X2>& src)
{   return *this = static_cast<const ArrayViewConst_<T2,X2>&>(src); }
template <class T2, class X2>
ArrayView_& operator=(const Array_<T2,X2>& src)
{   return *this = static_cast<const ArrayViewConst_<T2,X2>&>(src); }

template <class T2, class A2>
ArrayView_& operator=(const std::vector<T2,A2>& src) {
    avAssignIteratorDispatch(src.begin(), src.end(),
                             std::random_access_iterator_tag(),
                             "ArrayView_<T>::operator=(std::vector<T2>)");
    return *this;
}

ArrayView_& operator=(const T& fillValue) 
{   fill(fillValue); return *this; }

ArrayView_& fill(const T& fillValue) {
    for (T* d = begin(); d != end(); ++d)
        *d = fillValue; // using T::operator=(T)
    return *this;
}

void assign(size_type n, const T& fillValue) {
    SimTK_ERRCHK2(n == size(), "ArrayView_<T>::assign(n,value)",
        "Assignment to an ArrayView is permitted only if the source"
        " is the same size. Here n==%llu element(s) but the"
        " ArrayView has a fixed size of %llu.", 
        this->ull(n), this->ull(size()));

    fill(fillValue);
}

template <class T2>
void assign(const T2* first, const T2* last1) {
    const char* methodName = "ArrayView_<T>::assign(T2* first, T2* last1)";
    SimTK_ERRCHK((first&&last1)||(first==last1), methodName, 
        "One of the source pointers was null (0); either both must be"
        " non-null or both must be null.");
    // Valid pointers are random access iterators.
    avAssignIteratorDispatch(first, last1, std::random_access_iterator_tag(),
                             methodName);
}

// Watch out for integral types matching this signature -- they must be
// forwarded to the assign(n, fillValue) signature instead.
template <class Iter>
void assign(const Iter& first, const Iter& last1)
{   avAssignDispatch(first,last1,typename IsIntegralType<Iter>::Result(),
                     "ArrayView_<T>::assign(Iter first, Iter last1)"); }
//------------------------------------------------------------------------------

const T& operator[](index_type i) const {return this->CBase::operator[](i);}

T& operator[](index_type i) {return const_cast<T&>(this->CBase::operator[](i));}

const T& at(index_type i) const {return this->CBase::at(i);}

T& at(index_type i) {return const_cast<T&>(this->CBase::at(i));}

const T& getElt(index_type i) const {return this->CBase::getElt(i);}
T& updElt(index_type i) {return const_cast<T&>(this->CBase::getElt(i));}

const T& front() const {return this->CBase::front();} 

T& front() {return const_cast<T&>(this->CBase::front());}

const T& back() const {return this->CBase::back();}

T& back() {return const_cast<T&>(this->CBase::back());}

ArrayView_ operator()(index_type index, size_type length) {
    const size_type ix(index);
    SimTK_ERRCHK2(isSizeInRange(ix, size()), "ArrayView_<T>(index,length)",
        "For this operator, we must have 0 <= index <= size(), but"
        " index==%llu and size==%llu.", this->ull(ix), ullSize());
    SimTK_ERRCHK2(isSizeInRange(length, size_type(size()-ix)), 
        "ArrayView_<T>(index,length)", 
        "This operator requires 0 <= length <= size()-index, but"
        " length==%llu and size()-index==%llu.",this->ull(length),this->ull(size()-ix));

    return ArrayView_(data()+ix, data()+ix+length);
}
ArrayView_ updSubArray(index_type index, size_type length)
{   return (*this)(index,length); }
//------------------------------------------------------------------------------

const T* cbegin() const {return this->CBase::cbegin();}
const T* begin() const {return this->CBase::cbegin();}
T* begin() {return const_cast<T*>(this->CBase::cbegin());}

const T* cend() const {return this->CBase::cend();}
const T* end() const {return this->CBase::cend();}
T* end() {return const_cast<T*>(this->CBase::cend());}

const_reverse_iterator crbegin() const {return this->CBase::crbegin();}
const_reverse_iterator rbegin() const {return this->CBase::crbegin();} 
reverse_iterator rbegin() {return reverse_iterator(end());}

const_reverse_iterator crend() const {return this->CBase::crend();}
const_reverse_iterator rend() const {return this->CBase::crend();}
reverse_iterator rend() {return reverse_iterator(begin());}

const T* cdata() const {return this->CBase::cdata();}
const T* data() const {return this->CBase::cdata();}
T* data() {return const_cast<T*>(this->CBase::cdata());}
//------------------------------------------------------------------------------

// Note: these have to be explicitly forwarded to the base class methods
// in order to keep gcc from complaining. Note that the "this->" is 
// apparently necessary in order to permit delayed definition of templatized 
// methods. Doxygen picks up the comments from the base class.

size_type size()      const {return this->CBase::size();}
size_type max_size()  const {return this->CBase::max_size();}
bool      empty()     const {return this->CBase::empty();}
size_type capacity()  const {return this->CBase::capacity();}
size_type allocated() const {return this->CBase::allocated();}
bool      isOwner()   const {return this->CBase::isOwner();}
//------------------------------------------------------------------------------
                                   private:
//------------------------------------------------------------------------------
// no data members are allowed

//------------------------------------------------------------------------------
//                       ARRAY VIEW ASSIGN DISPATCH
// This is the assign() implementation for ArrayView_ when the class that 
// matched the alleged InputIterator template argument turned out to be one of 
// the integral types in which case this should match the assign(n, fillValue) 
// signature.
template <class IntegralType>
void avAssignDispatch(IntegralType n, IntegralType v, TrueType isIntegralType,
                      const char*) 
{   assign(size_type(n), value_type(v)); }

// This is the assign() implementation for ArrayView_ when the class that 
// matched the alleged InputIterator template argument is NOT an integral type 
// and may very well be an iterator. 
template <class InputIterator> 
void avAssignDispatch(const InputIterator& first, const InputIterator& last1, 
                      FalseType isIntegralType, const char* methodName) 
{   avAssignIteratorDispatch(first, last1, 
        typename std::iterator_traits<InputIterator>::iterator_category(),
        methodName); }

// This is the assign() implementation for a plain input_iterator
// (i.e., not a forward, bidirectional, or random access iterator). These
// have the unfortunate property that we can't count the elements in advance.
// Here we're going to complain if there aren't enough; but will simply stop
// when we get size() elements and not insist that the input stream reached
// the supplied last1 iterator. Semantics is elementwise assignment.
template <class InputIterator>
void avAssignIteratorDispatch(const InputIterator& first, 
                              const InputIterator& last1, 
                              std::input_iterator_tag, 
                              const char* methodName) 
{
    T* p = begin();
    InputIterator src = first;
    while (src != last1 && p != end())
        *p++ = *src++; // call T's assignment operator

    // p now points just after the last element that was copied.
    const size_type nCopied = size_type(p - begin());
    SimTK_ERRCHK2_ALWAYS(nCopied == size(), methodName,
        "The supplied input_iterator provided only %llu elements but this"
        " ArrayView has a fixed size of %llu elements.",
        this->ull(nCopied), ullSize());

    // We don't care if there are still more input elements available.
}

// This is the assign() implementation that works for forward and bidirectional
// iterators, but is not used for random_access_iterators. Here we'll count
// the elements as we copy them and complain at the end if there were too
// few or too many.
template <class ForwardIterator>
void avAssignIteratorDispatch(const ForwardIterator& first, 
                              const ForwardIterator& last1,
                              std::forward_iterator_tag, 
                              const char* methodName) 
{
    T* p = begin();
    ForwardIterator src = first;
    while (src != last1 && p != end())
        *p++ = *src++; // call T's assignment operator

    // p now points just after the last element that was copied.
    const size_type nCopied = size_type(p - begin());
    SimTK_ERRCHK2_ALWAYS(nCopied == size(), methodName,
        "The supplied forward_ or bidirectional_iterator source range provided"
        " only %llu elements but this ArrayView has a fixed size of"
        " %llu elements.", this->ull(nCopied), ullSize());

    // Make sure we ran out of source elements.
    SimTK_ERRCHK1_ALWAYS(src == last1, methodName,
        "The supplied forward_ or bidirectional_iterator source range"
        " contained too many elements; this ArrayView has a fixed size of"
        " %llu elements.", ullSize());
}

// This is the assign() implementation that works for random_access_iterators
// including ordinary pointers. Here we check the number of elements in advance
// and complain if the source and destination aren't the same size. The 
// copying loop can be done faster in this case.
template <class RandomAccessIterator>
void avAssignIteratorDispatch(const RandomAccessIterator& first, 
                              const RandomAccessIterator& last1,
                              std::random_access_iterator_tag, 
                              const char* methodName) 
{
    SimTK_ERRCHK2_ALWAYS(this->isSameSize(last1-first), methodName,
        "Assignment to an ArrayView is permitted only if the source"
        " is the same size. Here the source had %llu element(s) but the"
        " ArrayView has a fixed size of %llu.", 
        this->ull(last1-first), this->ull(size()));

    SimTK_ERRCHK_ALWAYS(!this->overlapsWithData(first,last1), methodName,
        "Source range can't overlap with the destination data.");

    T* p = begin();
    RandomAccessIterator src = first;
    while (p != end())
        *p++ = *src++; // call T's assignment operator
}


//------------------------------------------------------------------------------
// The following private methods are protected methods in the ArrayViewConst_ 
// base class, so they should not need repeating here. However, we explicitly 
// forward to the base methods to avoid gcc errors. The gcc complaint
// is due to their not depending on any template parameters; the "this->"
// apparently fixes that problem.

packed_size_type psize()      const {return this->CBase::psize();}
packed_size_type pallocated() const {return this->CBase::pallocated();}

// This just cast sizes to unsigned long long so that we can do comparisons
// without getting warnings.
unsigned long long ullSize()     const {return this->CBase::ullSize();}
unsigned long long ullCapacity() const {return this->CBase::ullCapacity();}
unsigned long long ullMaxSize()  const {return this->CBase::ullMaxSize();}
// This is the index type name and is handy for error messages to explain
// why some size was too big.
const char* indexName() const   {return this->CBase::indexName();}
};





//==============================================================================
//                               CLASS Array_
//==============================================================================
template <class T, class X> class Array_ : public ArrayView_<T,X> {
    typedef ArrayView_<T,X>      Base;
    typedef ArrayViewConst_<T,X> CBase;
public:


//------------------------------------------------------------------------------
// Doxygen picks up individual descriptions from the base class.
/*{*/
typedef T                                               value_type;
typedef X                                               index_type;
typedef T*                                              pointer;
typedef const T*                                        const_pointer;
typedef T&                                              reference;
typedef const T&                                        const_reference;
typedef T*                                              iterator;
typedef const T*                                        const_iterator;
typedef std::reverse_iterator<iterator>                 reverse_iterator;
typedef std::reverse_iterator<const_iterator>           const_reverse_iterator;
typedef typename ArrayIndexTraits<X>::size_type         size_type;
typedef typename ArrayIndexTraits<X>::difference_type   difference_type;
typedef typename ArrayIndexPackType<size_type>::packed_size_type 
                                                        packed_size_type;
/*}*/

//------------------------------------------------------------------------------
/*{*/

Array_() : Base() {}

explicit Array_(size_type n) : Base() {
    SimTK_SIZECHECK(n, max_size(), "Array_<T>::ctor(n)");
    allocateNoConstruct(n);
    defaultConstruct(data(), data()+n);
    setSize(n);
}

Array_(size_type n, const T& initVal) : Base() {
    SimTK_SIZECHECK(n, max_size(), "Array_<T>::ctor(n,T)");
    setSize(n);
    allocateNoConstruct(size());
    fillConstruct(begin(), cend(), initVal);
}
template <class InputIter>
Array_(const InputIter& first, const InputIter& last1) : Base() {
    ctorDispatch(first,last1,typename IsIntegralType<InputIter>::Result());
}

template <class T2>
Array_(const T2* first, const T2* last1) : Base() {
    SimTK_ERRCHK((first&&last1)||(first==last1), "Array_<T>(first,last1)", 
        "Pointers must be non-null unless they are both null.");
    SimTK_ERRCHK3(this->isSizeOK(last1-first), "Array_<T>(first,last1)",
        "Source has %llu elements but this array is limited to %llu"
        " elements by its index type %s.",
        this->ull(last1-first), ullMaxSize(), indexName());

    setSize(size_type(last1-first));
    allocateNoConstruct(size());
    copyConstruct(begin(), cend(), first);
}

template <class T2>
explicit Array_(const std::vector<T2>& v) : Base() { 
    if (v.empty()) return;

    SimTK_ERRCHK3(this->isSizeOK(v.size()), "Array_<T>::ctor(std::vector<T2>)",
        "The source std::vector's size %llu is too big for this array which"
        " is limited to %llu elements by its index type %s.",
        this->ull(v.size()), ullMaxSize(), indexName());

    // Call the above constructor, making sure to use pointers into the
    // vector's data rather than the iterators begin() and end() in case
    // they are different types.
    new (this) Array_(&v.front(), (&v.back())+1);
}

Array_(const Array_& src) : Base() {
    setSize(src.size());
    allocateNoConstruct(size());
    copyConstruct(begin(), cend(), src.data());
}

template <class T2, class X2>
Array_(const Array_<T2,X2>& src) : Base() {
    new (this) Array_(src.begin(), src.cend()); // see above
}

Array_(T* first, const T* last1, const DontCopy&) : Base(first,last1) {}

template <class A>
Array_(std::vector<T,A>& v, const DontCopy&) : Base(v) {}

~Array_() {
    deallocate();
}

Array_& deallocate() {
    if (allocated()) { // owner with non-zero allocation
        clear(); // each element is destructed; size()=0; allocated() unchanged
        deallocateNoDestruct(); // free data(); allocated()=0
    }
    this->Base::disconnect(); // clear the handle
    return *this;
}

//------------------------------------------------------------------------------

void assign(size_type n, const T& fillValue) {
    SimTK_ERRCHK3(this->isSizeOK(n), "Array_<T>::assign(n,value)",
        "Requested size %llu is too big for this array which is limited"
        " to %llu elements by its index type %s.",
        this->ull(n), ullMaxSize(), indexName());

    SimTK_ERRCHK2(isOwner() || n==size(), "Array_<T>::assign(n,value)",
        "Requested size %llu is not allowed because this is a non-owner"
        " array of fixed size %llu.", this->ull(n), this->ull(size()));

    if (!isOwner())
        this->Base::fill(fillValue);
    else {
        clear(); // all elements destructed; allocation unchanged
        reallocateIfAdvisable(n); // change size if too small or too big
        fillConstruct(data(), cdata()+n, fillValue);
        setSize(n);
    }
}

void fill(const T& fillValue) {this->Base::fill(fillValue);}


template <class T2>
void assign(const T2* first, const T2* last1) {
    const char* methodName = "Array_<T>::assign(T2* first, T2* last1)";
    SimTK_ERRCHK((first&&last1)||(first==last1), methodName, 
        "Pointers must be non-null unless they are both null.");
    SimTK_ERRCHK(!this->overlapsWithData(first,last1), methodName,
        "Source range can't overlap the current array contents.");
    // Pointers are random access iterators.
    assignIteratorDispatch(first,last1,std::random_access_iterator_tag(),
                           methodName);
}


template <class Iter>
void assign(const Iter& first, const Iter& last1) {
    assignDispatch(first,last1,typename IsIntegralType<Iter>::Result(),
                   "Array_<T>::assign(Iter first, Iter last1)");
}

Array_& operator=(const Array_& src) {
    if (this != &src)
        assignIteratorDispatch(src.begin(), src.end(), 
                               std::random_access_iterator_tag(),
                               "Array_<T>::operator=(Array_<T>)");
    return *this;
}

template <class T2, class X2>
Array_& operator=(const Array_<T2,X2>& src) {
    assignIteratorDispatch(src.begin(), src.end(), 
                           std::random_access_iterator_tag(),
                           "Array_<T>::operator=(Array_<T2,X2>)");
    return *this;
}


template <class T2, class A>
Array_& operator=(const std::vector<T2,A>& src) {
    assignIteratorDispatch(src.begin(), src.end(), 
                           std::random_access_iterator_tag(),
                           "Array_<T>::operator=(std::vector)");
    return *this;
}

void swap(Array_& other) {
    T* const pTmp=data(); setData(other.data()); other.setData(pTmp);
    size_type nTmp=size(); setSize(other.size()); other.setSize(nTmp);
    nTmp=allocated(); setAllocated(other.allocated()); other.setAllocated(nTmp);
}

Array_& adoptData(T* newData, size_type dataSize, 
                  size_type dataCapacity) 
{
    SimTK_SIZECHECK(dataCapacity, max_size(), "Array_<T>::adoptData()");
    SimTK_ERRCHK2(dataSize <= dataCapacity, "Array_<T>::adoptData()", 
        "Specified data size %llu was greater than the specified data"
        " capacity of %llu.", this->ull(dataSize), this->ull(dataCapacity));
    SimTK_ERRCHK(newData || dataCapacity==0, "Array_<T>::adoptData()",
        "A null data pointer is allowed only if the size and capacity are"
        " specified as zero.");
    SimTK_ERRCHK(!this->overlapsWithData(newData, newData+dataSize), 
        "Array_<T>::adoptData()",
        "The new data can't overlap with the old data.");

    deallocate();
    setData(newData);
    setSize(dataSize);
    setAllocated(dataCapacity);
    return *this;
}
Array_& adoptData(T* newData, size_type dataSize) 
{   return adoptData(newData, dataSize, dataSize); }


Array_& shareData(T* newData, size_type dataSize) {
    SimTK_SIZECHECK(dataSize, max_size(), "Array_<T>::shareData()");
    SimTK_ERRCHK(newData || dataSize==0, "Array_<T>::shareData()",
        "A null data pointer is allowed only if the size is zero.");
    SimTK_ERRCHK(!this->overlapsWithData(newData, newData+dataSize), 
        "Array_<T>::shareData()",
        "The new data can't overlap with the old data.");

    deallocate();
    setData(newData);
    setSize(dataSize);
    setAllocated(0); // indicates shared data
    return *this;
}

Array_& shareData(T* first, const T* last1) {
    SimTK_ERRCHK3(this->isSizeOK(last1-first), "Array_<T>::shareData(first,last1)",
        "Requested size %llu is too big for this array which is limited"
        " to %llu elements by its index type %s.",
        this->ull(last1-first), ullMaxSize(), indexName());
    return shareData(first, size_type(last1-first));
}

//------------------------------------------------------------------------------

// Note: these have to be explicitly forwarded to the base class methods
// in order to keep gcc from complaining. Note that the "this->" is 
// apparently necessary in order to permit delayed definition of templatized 
// methods.

size_type size() const {return this->CBase::size();}
size_type max_size() const {return this->CBase::max_size();}
bool empty() const {return this->CBase::empty();}
size_type capacity() const {return this->CBase::capacity();}

void resize(size_type n) {
    if (n == size()) return;

    SimTK_ERRCHK2(isOwner(), "Array_<T>::resize(n)",
        "Requested size change to %llu is not allowed because this is a"
        " non-owner array of fixed size %llu.", this->ull(n), this->ull(size()));

    if (n < size()) {
        erase(data()+n, cend());
        return;
    }
    // n > size()
    reserve(n);
    defaultConstruct(data()+size(), cdata()+n); // data() has changed
    setSize(n);
}

void resize(size_type n, const T& initVal) {
    if (n == size()) return;

    SimTK_ERRCHK2(isOwner(), "Array_<T>::resize(n,value)",
        "Requested size change to %llu is not allowed because this is a"
        " non-owner array of fixed size %llu.", this->ull(n), this->ull(size()));

    if (n < size()) {
        erase(data()+n, cend());
        return;
    }
    // n > size()
    reserve(n);
    fillConstruct(data()+size(), cdata()+n, initVal);
    setSize(n);
}

void reserve(size_type n) {
    if (capacity() >= n)
        return;

    SimTK_ERRCHK2(isOwner(), "Array_<T>::reserve()",
        "Requested capacity change to %llu is not allowed because this is a"
        " non-owner array of fixed size %llu.", this->ull(n), this->ull(size()));

    T* newData = allocN(n); // no construction yet
    copyConstructThenDestructSource(newData, newData+size(), data());
    freeN(data());
    setData(newData);
    setAllocated(n);
}

void shrink_to_fit() {
    // Allow 25% slop, but note that if size()==0 this will always reallocate
    // unless capacity is already zero.
    if (capacity() - size()/4 <= size()) // avoid overflow if size() near max
        return;
    T* newData = allocN(size());
    copyConstructThenDestructSource(newData, newData+size(), data());
    deallocateNoDestruct(); // data()=0, allocated()=0, size() unchanged
    setData(newData);
    setAllocated(size());
}

size_type allocated() const {return this->CBase::allocated();}
bool isOwner() const {return this->CBase::isOwner();}
//------------------------------------------------------------------------------

const T* cbegin() const {return this->CBase::cbegin();}
const T* begin() const {return this->CBase::cbegin();}
T* begin() {return this->Base::begin();}

const T* cend() const {return this->CBase::cend();}
const T* end() const {return this->CBase::cend();}
T* end() {return this->Base::end();}

const_reverse_iterator crbegin() const {return this->CBase::crbegin();}
const_reverse_iterator rbegin() const {return this->CBase::crbegin();} 
reverse_iterator rbegin() {return this->Base::rbegin();}

const_reverse_iterator crend() const {return this->CBase::crend();}
const_reverse_iterator rend() const {return this->CBase::crend();}
reverse_iterator rend() {return this->Base::rend();}

const T* cdata() const {return this->CBase::cdata();}
const T* data() const {return this->CBase::cdata();}
T* data() {return this->Base::data();}

const T& operator[](index_type i) const {return this->CBase::operator[](i);}

T& operator[](index_type i) {return this->Base::operator[](i);}

const T& at(index_type i) const {return this->CBase::at(i);}

T& at(index_type i) {return const_cast<T&>(this->Base::at(i));}

const T& getElt(index_type i) const {return this->CBase::getElt(i);}
T& updElt(index_type i) {return this->Base::updElt(i);}

const T& front() const {return this->CBase::front();} 

T& front() {return const_cast<T&>(this->Base::front());}

const T& back() const {return this->CBase::back();}

T& back() {return const_cast<T&>(this->Base::back());}

ArrayViewConst_<T,X> operator()(index_type index, size_type length) const
{   return CBase::operator()(index,length); }
ArrayViewConst_<T,X> getSubArray(index_type index, size_type length) const
{   return CBase::getSubArray(index,length); }

ArrayView_<T,X> operator()(index_type index, size_type length)
{   return Base::operator()(index,length); }
ArrayView_<T,X> updSubArray(index_type index, size_type length)
{   return Base::updSubArray(index,length); }
//------------------------------------------------------------------------------

void push_back(const T& value) {
    if (pallocated() == psize())
        growAtEnd(1,"Array_<T>::push_back(value)");
    copyConstruct(end(), value);
    incrSize();
}

void push_back() {
    if (pallocated() == psize())
        growAtEnd(1,"Array_<T>::push_back()");
    defaultConstruct(end());
    incrSize();
}

T* raw_push_back() {
    if (pallocated() == psize())
        growAtEnd(1,"Array_<T>::raw_push_back()");
    T* const p = end();
    incrSize();
    return p;
}

void pop_back() {
    SimTK_ERRCHK(!empty(), "Array_<T>::pop_back()", "Array was empty.");
    destruct(&back());
    decrSize();
}

T* erase(T* first, const T* last1) {
    SimTK_ERRCHK(begin() <= first && first <= last1 && last1 <= end(),
    "Array<T>::erase(first,last1)", "Pointers out of range or out of order.");

    const size_type nErased = size_type(last1-first);
    SimTK_ERRCHK(isOwner() || nErased==0, "Array<T>::erase(first,last1)",
        "No elements can be erased from a non-owner array.");

    if (nErased) {
        destruct(first, last1); // Destruct the elements we're erasing.
        moveElementsDown(first+nErased, nErased); // Compress followers into the gap.
        setSize(size()-nErased);
    }
    return first;
}

T* erase(T* p) {
    SimTK_ERRCHK(begin() <= p && p < end(),
        "Array<T>::erase(p)", "Pointer must point to a valid element.");
    SimTK_ERRCHK(isOwner(), "Array<T>::erase(p)",
        "No elements can be erased from a non-owner array.");

    destruct(p);              // Destruct the element we're erasing.
    moveElementsDown(p+1, 1); // Compress followers into the gap.
    decrSize();
    return p;
}


T* eraseFast(T* p) {
    SimTK_ERRCHK(begin() <= p && p < end(),
        "Array<T>::eraseFast(p)", "Pointer must point to a valid element.");
    SimTK_ERRCHK(isOwner(), "Array<T>::eraseFast(p)",
        "No elements can be erased from a non-owner array.");

    destruct(p);
    if (p+1 != end()) 
        moveOneElement(p, &back());
    decrSize();
    return p;
}

void clear() {
    SimTK_ERRCHK(isOwner() || empty(), "Array_<T>::clear()", 
        "clear() is not allowed for a non-owner array.");
    destruct(begin(), end());
    setSize(0);
}


T* insert(T* p, size_type n, const T& value) {
    T* const gap = insertGapAt(p, n, "Array<T>::insert(p,n,value)");
    // Copy construct into the inserted elements and note the size change.
    fillConstruct(gap, gap+n, value);
    setSize(size()+n);
    return gap;
}

T* insert(T* p, const T& value)  {
    T* const gap = insertGapAt(p, 1, "Array<T>::insert(p,value)");
    // Copy construct into the inserted element and note the size change.
    copyConstruct(gap, value);
    incrSize();
    return gap;
}

template <class T2>
T* insert(T* p, const T2* first, const T2* last1) {
    const char* methodName = "Array_<T>::insert(T* p, T2* first, T2* last1)";
    SimTK_ERRCHK((first&&last1) || (first==last1), methodName, 
        "One of first or last1 was null; either both or neither must be null.");
    SimTK_ERRCHK(!this->overlapsWithData(first,last1), methodName,
        "Source range can't overlap with the current array contents.");
    // Pointers are random access iterators.
    return insertIteratorDispatch(p, first, last1,
                                  std::random_access_iterator_tag(),
                                  methodName);
}

template <class Iter>
T* insert(T* p, const Iter& first, const Iter& last1) {
    return insertDispatch(p, first, last1,
                          typename IsIntegralType<Iter>::Result(),
                          "Array_<T>::insert(T* p, Iter first, Iter last1)");
}
//------------------------------------------------------------------------------
                                  private:
//------------------------------------------------------------------------------


// This method is used when we have already decided we need to make room for 
// some new elements by reallocation, by creating a gap somewhere within the
// existing data. We'll issue an error message if this would violate the  
// max_size restriction (we can afford to do that even in a Release build since
// we don't expect to grow very often). Otherwise we'll allocate some more space
// and copy construct the existing elements into the new space, leaving a gap 
// where indicated. Note that this method does not change the current size but
// does change the capacity.
//
// The gapPos must point within the existing data with null OK if the array
// itself is null, and end() being OK any time although you should use the
// more specialized growAtEnd() method if you know that's what's happening.
//
// Don't call this with a gapSz of zero.
T* growWithGap(T* gapPos, size_type gapSz, const char* methodName) {
    assert(gapSz > 0); // <= 0 is a bug, not a user error

    // Note that gapPos may be null if begin() and end() are also.
    SimTK_ERRCHK(begin() <= gapPos && gapPos <= end(), methodName, 
        "Given insertion point is not valid for this array.");

    // Get some new space of a reasonable size.
    setAllocated(calcNewCapacityForGrowthBy(gapSz, methodName));
    T* newData   = allocN(allocated());

    // How many elements will be before the gap?
    const size_type nBefore = (size_type)(gapPos-begin());

    // Locate the gap in the new space allocation.
    T* newGap    = newData + nBefore;
    T* newGapEnd = newGap  + gapSz; // one past the last element in the gap

    // Copy elements before insertion point; destruct source as we go.
    copyConstructThenDestructSource(newData,   newGap,        data());
    // Copy/destruct elements at and after insertion pt; leave gapSz gap.
    copyConstructThenDestructSource(newGapEnd, newData+size(), gapPos);

    // Throw away the old data and switch to the new.
    freeN(data()); setData(newData);
    return newGap;
}

// Same as growWithGap(end(), n, methodName); see above.
void growAtEnd(size_type n, const char* methodName) {
    assert(n > 0); // <= 0 is a bug, not a user error
    // Get some new space of a reasonable size.
    setAllocated(calcNewCapacityForGrowthBy(n, methodName));
    T* newData   = allocN(allocated());
    // Copy all the elements; destruct source as we go.
    copyConstructThenDestructSource(newData, newData+size(), data());
    // Throw away the old data and switch to the new.
    freeN(data()); setData(newData);
}

// This method determines how much we should increase the array's capacity
// when asked to insert n elements, due to an insertion or push_back. We will
// generally allocate more new space than requested, in anticipation of
// further insertions. This has to be based on the current size so that only
// log(n) reallocations are performed to insert n elements one at a time. Our
// policy here is to at least double the capacity unless that would exceed 
// max_size(). There is also a minimum amount of allocation we'll do if the 
// current size is zero or very small. 
size_type calcNewCapacityForGrowthBy(size_type n, const char* methodName) const {
    SimTK_ERRCHK3_ALWAYS(isGrowthOK(n), methodName,
        "Can't grow this Array by %llu element(s) because it would"
        " then exceed the max_size of %llu set by its index type %s.",
        (unsigned long long)n, ullMaxSize(), indexName());

    // At this point we know that capacity()+n <= max_size().
    const size_type mustHave = capacity() + n;

    // Be careful not to overflow size_type as you could if you 
    // double capacity() rather than halving max_size().
    const size_type wantToHave = capacity() <= max_size()/2 
                                    ? 2*capacity() 
                                    : max_size();

    const size_type newCapacity = std::max(std::max(mustHave,wantToHave),
                                           minAlloc());
    return newCapacity;
}

// Insert an unconstructed gap of size n beginning at position p. The return
// value is a pointer to the first element in the gap. It will be p if no
// reallocation occurred, otherwise it will be pointing into the new data.
// On return size() will be unchanged although allocated() may be bigger.
T* insertGapAt(T* p, size_type n, const char* methodName) {
    // Note that p may be null if begin() and end() are also.
    SimTK_ERRCHK(begin() <= p && p <= end(), methodName, 
        "Given insertion point is not valid for this Array.");

    if (n==0) return p; // nothing to do

    SimTK_ERRCHK_ALWAYS(isOwner(), methodName,
        "No elements can be inserted into a non-owner array.");

    // Determine the number of elements before the insertion point and
    // the number at or afterwards (those must be moved up by one slot).
    const size_type before = (size_type)(p-begin()), after = (size_type)(end()-p);

    // Grow the container if necessary. Note that if we have to grow we
    // can create the gap at the same time we copy the old elements over
    // to the new space.
    if (capacity() >= size()+n) {
        moveElementsUp(p, n); // leave a gap at p
    } else { // need to grow
        setAllocated(calcNewCapacityForGrowthBy(n, methodName));
        T* newdata = allocN(allocated());
        // Copy the elements before the insertion point, and destroy source.
        copyConstructThenDestructSource(newdata, newdata+before, data());
        // Copy the elements at and after the insertion point, leaving a gap
        // of n elements.
        copyConstructThenDestructSource(newdata+before+n,
                                        newdata+before+n+after,
                                        p); // i.e., pData+before
        p = newdata + before; // points into newdata now
        freeN(data());
        setData(newdata);
    }

    return p;
}

//------------------------------------------------------------------------------
//                           CTOR DISPATCH
// This is the constructor implementation for when the class that matches
// the alleged InputIterator type turns out to be one of the integral types
// in which case this should be the ctor(n, initValue) constructor.
template <class IntegralType> void
ctorDispatch(IntegralType n, IntegralType v, TrueType isIntegralType) {
    new(this) Array_(size_type(n), value_type(v));
}

// This is the constructor implementation for when the class that matches
// the alleged InputIterator type is NOT an integral type and may very well
// be an iterator. In that case we split into iterators for which we can
// determine the number of elements in advance (forward, bidirectional,
// random access) and input iterators, for which we can't. Note: iterator
// types are arranged hierarchically random->bi->forward->input with each
// deriving from the one on its right, so the forward iterator tag also
// matches bi and random.
template <class InputIterator> void
ctorDispatch(const InputIterator& first, const InputIterator& last1, 
             FalseType isIntegralType) 
{   ctorIteratorDispatch(first, last1, 
        typename std::iterator_traits<InputIterator>::iterator_category()); }

// This is the slow generic ctor implementation for any iterator that can't
// make it up to forward_iterator capability (that is, an input_iterator).
// The issue here is that we can't advance the iterator to count the number
// of elements before allocating because input_iterators are consumed when
// reference so we can't go back to look. That means we may have to reallocate
// memory log n times as we "push back" these elements onto the array.
template <class InputIterator> void
ctorIteratorDispatch(const InputIterator& first, const InputIterator& last1, 
                     std::input_iterator_tag) 
{
    InputIterator src = first;
    while (src != last1) {
        // We can afford to check this always since we are probably doing I/O.
        // Throwing an exception in a constructor is tricky, though -- this
        // won't go through the Array_ destructor although it will call the
        // Base (ArrayView_) destructor. Since we have already allocated
        // some space, we must call deallocate() manually.
        if (size() == max_size()) {
            deallocate();
            SimTK_ERRCHK2_ALWAYS(!"too many elements",
                "Array_::ctor(InputIterator first, InputIterator last1)",
                "There were still source elements available when the array"
                " reached its maximum size of %llu as determined by its index"
                " type %s.", ullMaxSize(), indexName());
        }
        push_back(*src++);
    }
}

// This is the faster constructor implementation for iterator types for which
// we can calculate the number of elements in advance. This will be optimal
// for a random access iterator since we can count in constant time, but for
// forward or bidirectional we'll have to advance n times to count and then
// go back again to do the copy constructions.
template <class ForwardIterator> void
ctorIteratorDispatch(const ForwardIterator& first, const ForwardIterator& last1, 
                     std::forward_iterator_tag) 
{
    typedef typename std::iterator_traits<ForwardIterator>::difference_type
        difference_type;
    // iterDistance() is constant time for random access iterators, but 
    // O(last1-first) for forward and bidirectional since it has to increment 
    // to count how far apart they are.
    const difference_type nInput = this->iterDistance(first,last1);

    SimTK_ERRCHK(nInput >= 0, 
        "Array_(ForwardIterator first, ForwardIterator last1)", 
        "Iterators were out of order.");

    SimTK_ERRCHK3(this->isSizeOK(nInput), 
        "Array_(ForwardIterator first, ForwardIterator last1)",
        "Source has %llu elements but this array is limited to %llu"
        " elements by its index type %s.",
        this->ull(nInput), ullMaxSize(), indexName());

    const size_type n = size_type(nInput);
    setSize(n);
    allocateNoConstruct(n);
    copyConstruct(data(), data()+n, first);
}

//------------------------------------------------------------------------------
//                           INSERT DISPATCH
// This is the insert() implementation for when the class that matches
// the alleged InputIterator type turns out to be one of the integral types
// in which case this should be the insert(p, n, initValue) constructor.
template <class IntegralType> 
T* insertDispatch(T* p, IntegralType n, IntegralType v, 
                  TrueType isIntegralType, const char*) 
{   return insert(p, size_type(n), value_type(v)); }

// This is the insert() implementation for when the class that matches
// the alleged InputIterator type is NOT an integral type and may very well
// be an iterator. See ctorDispatch() above for more information.
template <class InputIterator> 
T* insertDispatch(T* p, const InputIterator& first, const InputIterator& last1, 
                  FalseType isIntegralType, const char* methodName) 
{   return insertIteratorDispatch(p, first, last1, 
        typename std::iterator_traits<InputIterator>::iterator_category(),
        methodName); }

// This is the slow generic insert implementation for any iterator that can't
// make it up to forward_iterator capability (that is, an input_iterator).
// See ctorIteratorDispatch() above for more information.
template <class InputIterator> 
T* insertIteratorDispatch(T* p, InputIterator first, InputIterator last1, 
                          std::input_iterator_tag, const char* methodName) 
{
    size_type nInserted = 0;
    while (first != last1) {
        // We can afford to check this always since we are probably doing I/O.
        SimTK_ERRCHK2_ALWAYS(size() < max_size(), methodName,
            "There were still source elements available when the array"
            " reached its maximum size of %llu as determined by its index"
            " type %s.", ullMaxSize(), indexName());
        p = insert(p, *first++);  // p may now point to reallocated memory
        ++p; ++nInserted;
    }
    // p now points just after the last inserted element; subtract the
    // number inserted to get a pointer to the first inserted element.
    return p-nInserted;
}

// This is the faster constructor implementation for iterator types for which
// we can calculate the number of elements in advance. This will be optimal
// for a random access iterator since we can count in constant time, but for
// forward or bidirectional we'll have to advance n times to count and then
// go back again to do the copy constructions.
template <class ForwardIterator>
T* insertIteratorDispatch(T* p, const ForwardIterator& first, 
                                const ForwardIterator& last1,
                                std::forward_iterator_tag,
                                const char* methodName) 
{
    typedef typename std::iterator_traits<ForwardIterator>::difference_type
        difference_type;
    // iterDistance() is constant time for random access iterators, but 
    // O(last1-first) for forward and bidirectional since it has to increment 
    // to count how far apart they are.
    const difference_type nInput = this->iterDistance(first,last1);

    SimTK_ERRCHK(nInput >= 0, methodName, "Iterators were out of order.");

    SimTK_ERRCHK3(isGrowthOK(nInput), methodName,
        "Source has %llu elements which would make this array exceed the %llu"
        " elements allowed by its index type %s.",
        this->ull(nInput), ullMaxSize(), indexName());

    const size_type n = size_type(nInput);
    p = insertGapAt(p, n, methodName);
    copyConstruct(p, p+n, first);
    setSize(size()+n);
    return p;
}

//------------------------------------------------------------------------------
//                           ASSIGN DISPATCH
// This is the assign() implementation for when the class that matches
// the alleged InputIterator type turns out to be one of the integral types
// in which case this should be the assign(n, initValue) constructor.
template <class IntegralType>
void assignDispatch(IntegralType n, IntegralType v, TrueType isIntegralType,
                    const char* methodName) 
{   assign(size_type(n), value_type(v)); }

// This is the assign() implementation for when the class that matches
// the alleged InputIterator type is NOT an integral type and may very well
// be an iterator. See ctorDispatch() above for more information.
template <class InputIterator> 
void assignDispatch(const InputIterator& first, const InputIterator& last1, 
                    FalseType isIntegralType, const char* methodName) 
{   assignIteratorDispatch(first, last1, 
        typename std::iterator_traits<InputIterator>::iterator_category(),
        methodName); }

// This is the slow generic implementation for a plain input_iterator
// (i.e., not a forward, bidirectional, or random access iterator). These
// have the unfortunate property that we can't count the elements in advance.
template <class InputIterator>
void assignIteratorDispatch(const InputIterator& first, 
                            const InputIterator& last1, 
                            std::input_iterator_tag, 
                            const char* methodName) 
{
    // TODO: should probably allow this and just blow up when the size()+1st
    // element is seen.
    SimTK_ERRCHK_ALWAYS(isOwner(), methodName,
        "Assignment to a non-owner array can only be done from a source"
        " designated with forward iterators or pointers because we"
        " must be able to verify that the source and destination sizes"
        " are the same.");

    clear(); // TODO: change space allocation here?
    InputIterator src = first;
    while (src != last1) {
        // We can afford to check this always since we are probably doing I/O.
        SimTK_ERRCHK2_ALWAYS(size() < max_size(), methodName,
            "There were still source elements available when the array"
            " reached its maximum size of %llu as determined by its index"
            " type %s.", ullMaxSize(), indexName());

        push_back(*src++);
    }
}

// This is the faster implementation that works for forward, bidirectional,
// and random access iterators including ordinary pointers. We can check here that the 
// iterators are in the right order, and that the source is not too big to
// fit in this array. Null pointer checks should be done prior to calling,
// however, since iterators in general aren't pointers.
template <class ForwardIterator>
void assignIteratorDispatch(const ForwardIterator& first, 
                            const ForwardIterator& last1,
                            std::forward_iterator_tag, 
                            const char* methodName) 
{
    typedef typename std::iterator_traits<ForwardIterator>::difference_type
        IterDiffType;
    // iterDistance() is constant time for random access iterators, but 
    // O(last1-first) for forward and bidirectional since it has to increment 
    // to count how far apart they are.
    const IterDiffType nInput = this->iterDistance(first,last1);

    SimTK_ERRCHK(nInput >= 0, methodName, "Iterators were out of order.");

    SimTK_ERRCHK3(this->isSizeOK(nInput), methodName,
        "Source has %llu elements but this Array is limited to %llu"
        " elements by its index type %s.",
        this->ull(nInput), ullMaxSize(), indexName());

    const size_type n = size_type(nInput);
    if (isOwner()) {
        // This is an owner Array; assignment is considered deallocation
        // followed by copy construction.

        clear(); // all elements destructed; allocation unchanged
        reallocateIfAdvisable(n); // change allocation if too small or too big
        copyConstruct(data(), cdata()+n, first);
        setSize(n);
    } else {
        // This is a non-owner Array. Assignment can occur only if the
        // source is the same size as the array, and the semantics are of
        // repeated assignment using T::operator=() not destruction followed
        // by copy construction.
        SimTK_ERRCHK2(n == size(), methodName,
            "Source has %llu elements which does not match the size %llu"
            " of the non-owner array it is being assigned into.",
            this->ull(n), ullSize());

        T* p = begin();
        ForwardIterator src = first;
        while (src != last1)
            *p++ = *src++; // call T's assignment operator
    }
}

// We are going to put a total of n elements into the Array (probably
// because of an assignment or resize) and we want the space allocation
// to be reasonable. That means of course that the allocation must be 
// *at least* n, but we also don't want it to be too big. Our policy
// here is that if it is currently less than twice what we need we
// won't reallocate, otherwise we'll shrink the space. When changing
// the size to zero or something very small we'll treat the Array as
// though its current size is minAlloc, meaning we won't reallocate
// if the existing space is less than twice minAlloc.
// nAllocated will be set appropriately; size() is not touched here.
// No constructors or destructors are called.
void reallocateIfAdvisable(size_type n) {
    if (allocated() < n || allocated()/2 > std::max(minAlloc(), n)) 
        reallocateNoDestructOrConstruct(n);
}


void allocateNoConstruct(size_type n) 
{   setData(allocN(n)); setAllocated(n); } // size() left unchanged
void deallocateNoDestruct() 
{   freeN(data()); setData(0); setAllocated(0); } // size() left unchanged
void reallocateNoDestructOrConstruct(size_type n)
{   deallocateNoDestruct(); allocateNoConstruct(n); }

// This sets the smallest allocation we'll do when growing.
size_type minAlloc() const 
{   return std::min(max_size(), size_type(4)); }

// Allocate without construction. Returns a null pointer if asked
// to allocate 0 elements. In Debug mode we'll fill memory with 
// all 1's as a bug catcher.
static T* allocN(size_type n) {
    if (n==0) return 0;
    unsigned char* newdata = new unsigned char[n * sizeof(T)];
    #ifndef NDEBUG
        unsigned char* b=newdata; 
        const unsigned char* e=newdata+(n*sizeof(T));
        while (b != e) *b++ = 0xff;
    #endif
    return reinterpret_cast<T*>(newdata);
}

// Free memory without calling T's destructor. Nothing happens if passed
// a null pointer.
static void freeN(T* p) {
    delete[] reinterpret_cast<char*>(p);
}

// default construct one element
static void defaultConstruct(T* p) {new(p) T();}
// default construct range [b,e)
static void defaultConstruct(T* b, const T* e) 
{   while (b!=e) new(b++) T(); }

// copy construct range [b,e) with repeats of a given value
static void fillConstruct(T* b, const T* e, const T& v)
{   while(b!=e) new(b++) T(v); }

// copy construct one element from a given value
static void copyConstruct(T* p, const T& v) {new(p) T(v);}
// copy construct range [b,e) from sequence of source values
static void copyConstruct(T* b, const T* e, T* src)
{   while(b!=e) new(b++) T(*src++); }
// Templatized copy construct will work if the source elements are
// assignment compatible with the destination elements.
template <class InputIterator>
static void copyConstruct(T* b, const T* e, InputIterator src)
{   while(b!=e) new(b++) T(*src++); }

// Copy construct range [b,e] from sequence of source values and
// destruct the source after it is copied. It's better to alternate
// copying and destructing than to do this in two passes since we
// will already have touched the memory.
static void copyConstructThenDestructSource(T* b, const T* e, T* src)
{   while(b!=e) {new(b++) T(*src); src++->~T();} }

// We have an element at from that we would like to move into the currently-
// unconstructed slot at to. Both from and to are expected to be pointing to
// elements within the currently allocated space. From's slot will be left
// unconstructed.
void moveOneElement(T* to, T* from) {
    assert(data() <= to   && to   < data()+allocated());
    assert(data() <= from && from < data()+allocated());
    copyConstruct(to, *from); 
    destruct(from);
}


// Move elements from p to end() down by n>0 places to fill an unconstructed 
// gap beginning at p-n. Any leftover space at the end will be unconstructed.
void moveElementsDown(T* p, size_type n) {
    assert(n > 0);
    for (; p != end(); ++p)
        moveOneElement(p-n,p);
}

// Move elements from p to end() up by n>0 places to make an unconstructed gap
// at [p,p+n). Note that this has to be done backwards so that we don't
// write on any elements until after they've been copied.
void moveElementsUp(T* p, size_type n) {
    assert(n > 0);
    T* src = end(); // points one past source element (begin()-1 not allowed)
    while (src != p) {
        --src; // now points to source
        moveOneElement(src+n, src);;
    }
}

// destruct one element
static void destruct(T* p) {p->~T();}
// destruct range [b,e)
static void destruct(T* b, const T* e)
{   while(b!=e) b++->~T(); }

// Check that growing this array by n elements wouldn't cause it to exceed
// its allowable maximum size.
template <class S>
bool isGrowthOK(S n) const
{   return this->isSizeOK(ullCapacity() + this->ull(n)); }

// The following private methods are protected methods in the ArrayView base 
// class, so they should not need repeating here. Howevr, we explicitly 
// forward to the Base methods to avoid gcc errors. The gcc complaint
// is due to their not depending on any template parameters; the "this->"
// apparently fixes that problem.

// These provide direct access to the data members.
packed_size_type psize()      const {return this->CBase::psize();}
packed_size_type pallocated() const {return this->CBase::pallocated();}

void setData(const T* p)        {this->CBase::setData(p);}
void setSize(size_type n)       {this->CBase::setSize(n);}
void incrSize()                 {this->CBase::incrSize();}
void decrSize()                 {this->CBase::decrSize();}
void setAllocated(size_type n)  {this->CBase::setAllocated(n);}
// This just cast sizes to unsigned long long so that we can do comparisons
// without getting warnings.
unsigned long long ullSize()     const {return this->CBase::ullSize();}
unsigned long long ullCapacity() const {return this->CBase::ullCapacity();}
unsigned long long ullMaxSize()  const {return this->CBase::ullMaxSize();}
// This is the index type name and is handy for error messages to explain
// why some size was too big.
const char* indexName() const   {return this->CBase::indexName();}
};


// This "private" static method is used to implement ArrayView's 
// fillArrayViewFromStream() and Array's readArrayFromStream() namespace-scope
// static methods, which are in turn used to implement ArrayView's and 
// Array's stream extraction operators ">>". This method has to be in the 
// header file so that we don't need to pass streams through the API, but it 
// is not intended for use by users and has no Doxygen presence, unlike 
// fillArrayFromStream() and readArrayFromStream() and (more commonly)
// the extraction operators.
template <class T, class X> static inline 
std::istream& readArrayFromStreamHelper
   (std::istream& in, bool isFixedSize, Array_<T,X>& out)
{
    // If already failed, bad, or eof, set failed bit and return without 
    // touching the Array.
    if (!in.good()) {in.setstate(std::ios::failbit); return in;}

    // If the passed-in Array isn't an owner, then we have to treat it as
    // a fixed size ArrayView regardless of the setting of the isFixedSize
    // argument.
    if (!out.isOwner())
        isFixedSize = true; // might be overriding the argument here

    // numRequired will be ignored unless isFixedSize==true.
    const typename Array_<T,X>::size_type 
        numRequired = isFixedSize ? out.size() : 0;

    if (!isFixedSize)
        out.clear(); // We're going to replace the entire contents of the Array.

    // Skip initial whitespace. If that results in eof this may be a successful
    // read of a 0-length, unbracketed Array. That is OK for either a
    // variable-length Array or a fixed-length ArrayView of length zero.
    std::ws(in); if (in.fail()) return in;
    if (in.eof()) {
        if (isFixedSize && numRequired != 0)
            in.setstate(std::ios_base::failbit); // zero elements not OK
        return in;
    }
    
    // Here the stream is good and the next character is non-white.
    assert(in.good());

    // Use this for raw i/o (peeks and gets).
    typename       std::iostream::int_type ch;
    const typename std::iostream::int_type EOFch = 
        std::iostream::traits_type::eof();

    // Now see if the sequence is bare or surrounded by (), [], or {}.
    bool lookForCloser = true;
    char openBracket, closeBracket;
    ch = in.peek(); if (in.fail()) return in;
    assert(ch != EOFch); // we already checked above

    openBracket = (char)ch;
    if      (openBracket=='(') {in.get(); closeBracket = ')';}
    else if (openBracket=='[') {in.get(); closeBracket = ']';}
    else if (openBracket=='{') {in.get(); closeBracket = '}';}
    else lookForCloser = false;

    // If lookForCloser is true, then closeBracket contains the terminating
    // delimiter, otherwise we're not going to quit until eof.

    // Eat whitespace after the opening bracket to see what's next.
    if (in.good()) std::ws(in);

    // If we're at eof now it must be because the open bracket was the
    // last non-white character in the stream, which is an error.
    if (!in.good()) {
        if (in.eof()) {
            assert(lookForCloser); // or we haven't read anything that could eof
            in.setstate(std::ios::failbit);
        }
        return in;
    }

    // istream is good and next character is non-white; ready to read first
    // value or terminator.

    // We need to figure out whether the elements are space- or comma-
    // separated and then insist on consistency.
    bool commaOK = true, commaRequired = false;
    bool terminatorSeen = false;
    X nextIndex(0);
    while (true) {
        char c;

        // Here at the top of this loop, we have already successfully read 
        // n=nextIndex values of type T. For fixed-size reads, it might be
        // the case that n==numRequired already, but we still may need to
        // look for a closing bracket before we can declare victory.
        // The stream is good() (not at eof) but it might be the case that 
        // there is nothing but white space left; we don't know yet because
        // if we have satisfied the fixed-size count and are not expecting
        // a terminator then we should quit without absorbing the trailing
        // white space.
        assert(in.good());

        // Look for closing bracket before trying to read value.
        if (lookForCloser) {
            // Eat white space to find the closing bracket.
            std::ws(in); if (!in.good()) break; // eof?
            ch = in.peek(); assert(ch != EOFch);
            if (!in.good()) break;
            c = (char)ch;
            if (c == closeBracket) {   
                in.get(); // absorb the closing bracket
                terminatorSeen = true; 
                break; 
            }
            // next char not a closing bracket; fall through
        }

        // We didn't look or didn't find a closing bracket. The istream is good 
        // but we might be looking at white space.

        // If we already got all the elements we want, break for final checks.
        if (isFixedSize && (nextIndex == numRequired))
            break; // that's a full count.

        // Look for comma before value, except the first time.
        if (commaOK && nextIndex != 0) {
            // Eat white space to find the comma.
            std::ws(in); if (!in.good()) break; // eof?
            ch = in.peek(); assert(ch != EOFch);
            if (!in.good()) break;
            c = (char)ch;
            if (c == ',') {
                in.get(); // absorb comma
                commaRequired = true; // all commas from now on
            } else { // next char not a comma
                if (commaRequired) // bad, e.g.: v1, v2, v3 v4 
                {   in.setstate(std::ios::failbit); break; }
                else commaOK = false; // saw: v1 v2 (no commas now)
            }
            if (!in.good()) break; // might be eof
        }

        // No closing bracket yet; don't have enough elements; skipped comma 
        // if any; istream is good; might be looking at white space.
        assert(in.good());

        // Now read in an element of type T.
        // The extractor T::operator>>() will ignore leading white space.
        if (!isFixedSize)
            out.push_back(); // grow by one (default consructed)
        in >> out[nextIndex]; if (in.fail()) break;
        ++nextIndex;

        if (!in.good()) break; // might be eof
    }

    // We will get here under a number of circumstances:
    //  - the fail bit is set in the istream, or
    //  - we reached eof
    //  - we saw a closing brace
    //  - we got all the elements we wanted (for a fixed-size read)
    // Note that it is possible that we consumed everything except some
    // trailing white space (meaning we're not technically at eof), but
    // for consistency with built-in operator>>()'s we won't try to absorb
    // that trailing white space.

    if (!in.fail()) {
        if (lookForCloser && !terminatorSeen)
            in.setstate(std::ios::failbit); // missing terminator

        if (isFixedSize && nextIndex != numRequired)
            in.setstate(std::ios::failbit); // wrong number of values
    }

    return in;
}



//------------------------------------------------------------------------------
//                          RELATED GLOBAL OPERATORS
//------------------------------------------------------------------------------
// These are logically part of the Array_<T,X> class but are not actually 
// class members; that is, they are in the SimTK namespace.

// Some of the serialization methods could have been member functions but 
// then an attempt to explicitly instantiate the whole Array_<T> class for
// some type T would fail if T did not support the requisite I/O operations
// even if those operations were never used. This came up when Chris Bruns was
// trying to wrap Array objects for Python, which requires explicit 
// instantiation.


template <class T, class X> inline void
writeUnformatted(std::ostream& o, const Array_<T,X>& v) {
    for (X i(0); i < v.size(); ++i) {
        if (i != 0) o << " ";
        writeUnformatted(o, v[i]);
    }
}


template <class T, class X> inline void
writeFormatted(std::ostream& o, const Array_<T,X>& v) {
    o << '(';
    for (X i(0); i < v.size(); ++i) {
        if (i != 0) o << ',';
        writeFormatted(o, v[i]);
    }
    o << ')';
}

template <class T, class X> inline 
std::ostream&
operator<<(std::ostream& o, 
           const ArrayViewConst_<T,X>& a) 
{
    o << '(';
    if (!a.empty()) {
        o << a.front();
        for (const T* p = a.begin()+1; p != a.end(); ++p)
            o << ',' << *p;
    }
    return o << ')';
} 


template <class T, class X> inline bool 
readUnformatted(std::istream& in, Array_<T,X>& v) {
    v.clear();
    T element;
    std::ws(in); // Make sure we're at eof if stream is all whitespace.
    while (!in.eof() && readUnformatted(in, element))
        v.push_back(element);
    return !in.fail(); // eof is expected and ok
}

template <class T, class X> inline bool 
readUnformatted(std::istream& in, ArrayView_<T,X>& v) {
    for (X i(0); i < v.size(); ++i)
        if (!readUnformatted(in, v[i])) return false;
    return true;
}


template <class T, class X> inline bool 
readFormatted(std::istream& in, Array_<T,X>& v) {
    return !readArrayFromStream(in,v).fail();
}

template <class T, class X> inline bool 
readFormatted(std::istream& in, ArrayView_<T,X>& v) {
    return !fillArrayViewFromStream(in,v).fail();
}

template <class T, class X> static inline 
std::istream& readArrayFromStream(std::istream& in, Array_<T,X>& out)
{   return readArrayFromStreamHelper<T,X>(in, false /*variable sizez*/, out); }



template <class T, class X> static inline 
std::istream& fillArrayFromStream(std::istream& in, Array_<T,X>& out)
{   return readArrayFromStreamHelper<T,X>(in, true /*fixed size*/, out); }

template <class T, class X> static inline 
std::istream& fillArrayViewFromStream(std::istream& in, ArrayView_<T,X>& out)
{   return readArrayFromStreamHelper<T,X>(in, true /*fixed size*/, out); }




template <class T, class X> inline
std::istream& operator>>(std::istream& in, Array_<T,X>& out) 
{   return readArrayFromStream<T,X>(in, out); }

template <class T, class X> inline
std::istream& operator>>(std::istream& in, ArrayView_<T,X>& out) 
{   return fillArrayViewFromStream<T,X>(in, out); }


template <class T1, class X1, class T2, class X2> inline bool 
operator==(const ArrayViewConst_<T1,X1>& a1, const ArrayViewConst_<T2,X2>& a2) {
    // Avoid warnings in size comparison by using common type.
    const ptrdiff_t sz1 = a1.end()-a1.begin();
    const ptrdiff_t sz2 = a2.end()-a2.begin();
    if (sz1 != sz2) return false;
    const T1* p1 = a1.begin();
    const T2* p2 = a2.begin();
    while (p1 != a1.end())
        if (!(*p1++ == *p2++)) return false;
    return true;
}
template <class T1, class X1, class T2, class X2> inline bool 
operator!=(const ArrayViewConst_<T1,X1>& a1, const ArrayViewConst_<T2,X2>& a2)
{   return !(a1 == a2); }

template <class T1, class X1, class T2, class X2> inline bool 
operator<(const ArrayViewConst_<T1,X1>& a1, const ArrayViewConst_<T2,X2>& a2) {
    const T1* p1 = a1.begin();
    const T2* p2 = a2.begin();
    while (p1 != a1.end() && p2 != a2.end()) {
        if (!(*p1 == *p2))
            return *p1 < *p2; // otherwise p1 > p2
        ++p1; ++p2;
    }
    // All elements were equal until one or both arrays ran out of elements.
    // a1 is less than a2 only if a1 ran out and a2 didn't.
    return p1 == a1.end() && p2 != a2.end();
}
template <class T1, class X1, class T2, class X2> inline bool 
operator>=(const ArrayViewConst_<T1,X1>& a1, const ArrayViewConst_<T2,X2>& a2)
{   return !(a1 < a2); }
template <class T1, class X1, class T2, class X2> inline bool 
operator>(const ArrayViewConst_<T1,X1>& a1, const ArrayViewConst_<T2,X2>& a2)
{   return a2 < a1; }
template <class T1, class X1, class T2, class X2> inline bool 
operator<=(const ArrayViewConst_<T1,X1>& a1, const ArrayViewConst_<T2,X2>& a2)
{   return !(a1 > a2); }

template <class T1, class X1, class T2, class A2> inline bool 
operator==(const ArrayViewConst_<T1,X1>& a1, const std::vector<T2,A2>& v2)
{   return a1 == ArrayViewConst_<T2,size_t>(v2); }

template <class T1, class A1, class T2, class X2> inline bool 
operator==(const std::vector<T1,A1>& v1, const ArrayViewConst_<T2,X2>& a2)
{   return a2 == v1; }

template <class T1, class X1, class T2, class A2> inline bool 
operator!=(const ArrayViewConst_<T1,X1>& a1, const std::vector<T2,A2>& v2)
{   return !(a1 == v2); }
template <class T1, class A1, class T2, class X2> inline bool 
operator!=(const std::vector<T1,A1>& v1, const ArrayViewConst_<T2,X2>& a2)
{   return !(a2 == v1); }

template <class T1, class X1, class T2, class A2> inline bool 
operator<(const ArrayViewConst_<T1,X1>& a1, const std::vector<T2,A2>& v2)
{   return a1 < ArrayViewConst_<T2,size_t>(v2); }

template <class T1, class A1, class T2, class X2> inline bool 
operator<(const std::vector<T1,A1>& v1, const ArrayViewConst_<T2,X2>& a2)
{   return ArrayViewConst_<T1,size_t>(v1) < a2; }

template <class T1, class X1, class T2, class A2> inline bool 
operator>=(const ArrayViewConst_<T1,X1>& a1, const std::vector<T2,A2>& v2)
{   return !(a1 < v2); }
template <class T1, class A1, class T2, class X2> inline bool 
operator>=(const std::vector<T1,A1>& v1, const ArrayViewConst_<T2,X2>& a2)
{   return !(v1 < a2); }

template <class T1, class X1, class T2, class A2> inline bool 
operator>(const ArrayViewConst_<T1,X1>& a1, const std::vector<T2,A2>& v2)
{   return v2 < a1; }
template <class T1, class A1, class T2, class X2> inline bool 
operator>(const std::vector<T1,A1>& v1, const ArrayViewConst_<T2,X2>& a2)
{   return a2 < v1; }

template <class T1, class X1, class T2, class A2> inline bool 
operator<=(const ArrayViewConst_<T1,X1>& a1, const std::vector<T2,A2>& v2)
{   return !(a1 > v2); }
template <class T1, class A1, class T2, class X2> inline bool 
operator<=(const std::vector<T1,A1>& v1, const ArrayViewConst_<T2,X2>& a2)
{   return !(v1 > a2); }

} // namespace SimTK

namespace std {
template <class T, class X> inline void
swap(SimTK::Array_<T,X>& a1, SimTK::Array_<T,X>& a2) {
    a1.swap(a2);
}

} // namespace std
  
#endif // SimTK_SimTKCOMMON_ARRAY_H_