<commit_before>// __BEGIN_LICENSE__
// 
// Copyright (C) 2006 United States Government as represented by the
// Administrator of the National Aeronautics and Space Administration
// (NASA).  All Rights Reserved.
// 
// Copyright 2006 <NAME>. All rights reserved.
// 
// This software is distributed under the NASA Open Source Agreement
// (NOSA), version 1.3.  The NOSA has been approved by the Open Source
// Initiative.  See the file COPYING at the top of the distribution
// directory tree for the complete NOSA document.
// 
// THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY
// KIND, EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT
// LIMITED TO, ANY WARRANTY THAT THE SUBJECT SOFTWARE WILL CONFORM TO
// SPECIFICATIONS, ANY IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR
// A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY WARRANTY THAT
// THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT
// DOCUMENTATION, IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE.
// 
// __END_LICENSE__

/// \file Vector.h
///
/// Provides the core Vector<> mathematical vector template classes.  
///
/// Currently I believe we support:
///   Dynamic (Vector<T>) and fixed (Vector<T,N>) vectors
///   Element access via vec(i) and vec[i]
///   Seamless conversion between compatibile vector types
///   Vector<T,N> construction with between zero and four initializer elements
///   Dynamic and fixed VectorProxy proxy objects
///   Subvector access via subvector()
///   Vector transposition via transpose()
///   Explicit expression evaluation via eval()
///   Printing of vectors to ostreams
///   Equality of vectors with or without an epsilon
///   Vector negation, addition, subtraction
///   Scalar multiplication and division
///   Elementwise vector addition, subtraction, multiplication, and division
///   Elementwise scalar addition, subtraction, multiplication, and division
///   Vector negation, addition, subtraction of transposed vectors
///   Scalar multiplication and division of transposed vectors
///   Elementwise comprison operations.
///   Norms via norm_1(), norm_2(), norm_2_sqr(), and norm_inf()
///   Sum and product of elements via sum() and prod()
///   Vector normalization via normalize()
///   Dot/inner product via dot_prod(v1,v2) or transpose(v1)*v2
///   Cross product of 3-element vectors via cross_prod()
///
#ifndef __VW_MATH_VECTOR_H__
#define __VW_MATH_VECTOR_H__

#include <boost/utility/enable_if.hpp>
#include <boost/type_traits.hpp>
#include <boost/mpl/if.hpp>
#include <boost/static_assert.hpp>
#include <boost/array.hpp>
#include <boost/iterator/iterator_facade.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/utility/result_of.hpp>

#include <vw/Core/Exception.h>
#include <vw/Core/TypeDeduction.h>
#include <vw/Core/Functors.h>

namespace vw {
namespace math {

  namespace bnu = boost::numeric::ublas;

  /// \cond INTERNAL
  // A wrapper around boost::array that provides a constructor taking 
  // a size parameter, which the BOOST uBLAS containers expect to be 
  // available.  With debugging disabled this whole thing should be 
  // optimized out of existence.
  template <class ElemT, int SizeN>
  class FixedArray : public boost::array<ElemT,SizeN> {
  public:
    FixedArray() {}
    FixedArray( unsigned size ) {
      VW_ASSERT( size==SizeN, LogicErr() << "Internal size mismatch in FixedArray." );
    }
  };
  /// \endcond

  /// A type function to compute the dimension of a vector expression 
  /// at compile time (or zero for dynamically-sized vectors).
  template <class VectorT>
  struct VectorSize {
    const static int value = 0;
  };


  // *******************************************************************
  // class VectorBase<VectorT>
  // A CRTP vector base class.
  // *******************************************************************

  /// A CRTP base class for vectors and vector expressions.  
  /// Provides a mechanism for restricting function arguments to 
  /// vectors, and provides the various arithmetic assignment 
  /// operators.
  template <class VectorT>
  struct VectorBase {

    /// Returns the derived implementation type.
    VectorT& impl() { return *static_cast<VectorT*>(this); }

    /// Returns the derived implementation type.
    VectorT const& impl() const { return *static_cast<VectorT const*>(this); }
    
    /// Sum-assignment operator
    template <class T>
    VectorT& operator+=( T const& v ) {
      return impl() = impl() + v;
    }

    /// Difference-assignment operator
    template <class T>
    VectorT& operator-=( T const& v ) {
      return impl() = impl() - v;
    }

    /// Product-assignment operator
    template <class T>
    VectorT& operator*=( T s ) {
      return impl() = impl() * s;
    }

    /// Quotient-assignment operator
    template <class T>
    VectorT& operator/=( T s ) {
      return impl() = impl() / s;
    }

  };


  // *******************************************************************
  // class Vector<ElemT,SizeN>
  // A statically-allocated fixed-dimension vector class.
  // *******************************************************************

  /// A fixed-dimension mathematical vector class.
  template <class ElemT, int SizeN = 0>
  class Vector : public VectorBase<Vector<ElemT,SizeN> >
  {
    typedef bnu::vector<ElemT,FixedArray<ElemT,SizeN> > core_type;
    core_type core_;
  public:
    typedef ElemT value_type;

    typedef ElemT& reference_type;
    typedef ElemT const& const_reference_type;

    typedef ElemT* iterator;
    typedef const ElemT* const_iterator;

    /// Constructs a vector of zeroes.
    Vector() : core_(SizeN) {
      core_.clear();
    }

    /// Constructs a vector whose first element is as given.
    Vector( ElemT e1 ) : core_(SizeN) {
      BOOST_STATIC_ASSERT( SizeN >= 1 );
      (*this)[0] = e1;
      for( unsigned i=1; i<SizeN; ++i ) (*this)[i] = ElemT();
    }

    /// Constructs a vector whose first two elements are as given.
    Vector( ElemT e1, ElemT e2 ) : core_(SizeN) {
      BOOST_STATIC_ASSERT( SizeN >= 2 );
      (*this)[0] = e1; (*this)[1] = e2;
      for( unsigned i=2; i<SizeN; ++i ) (*this)[i] = ElemT();
    }

    /// Constructs a vector whose first three elements are as given.
    Vector( ElemT e1, ElemT e2, ElemT e3 ) : core_(SizeN) {
      BOOST_STATIC_ASSERT( SizeN >= 3 );
      (*this)[0] = e1; (*this)[1] = e2; (*this)[2] = e3;
      for( unsigned i=3; i<SizeN; ++i ) (*this)[i] = ElemT();
    }

    /// Constructs a vector whose first four elements are as given.
    Vector( ElemT e1, ElemT e2, ElemT e3, ElemT e4 ) : core_(SizeN) {
      BOOST_STATIC_ASSERT( SizeN >= 4 );
      (*this)[0] = e1; (*this)[1] = e2; (*this)[2] = e3; (*this)[3] = e4;
      for( unsigned i=4; i<SizeN; ++i ) (*this)[i] = ElemT();
    }

    /// Constructs a vector from given densely-packed data.  This
    /// constructor copies the data.  If you wish to make a shallow
    /// proxy object instead, see vw::VectorProxy.
    Vector( const ElemT data[SizeN] ) : core_(SizeN) {
      std::copy( data, data+SizeN, begin() );
    }

    /// Standard copy constructor.
    Vector( Vector const& v ) : core_( v.core_ ) {}

    /// Generalized copy constructor, from arbitrary VW vector expressions.
    template <class T>
    Vector( VectorBase<T> const& v ) : core_(SizeN) { 
      VW_ASSERT( v.impl().size()==SizeN, ArgumentErr() << "Vector must have dimension " << SizeN << "." );
      std::copy( v.impl().begin(), v.impl().end(), begin() );
    }

    /// Generalized copy constructor, from arbitrary uBLAS vector expressions.
    template <class T>
    Vector( bnu::vector_expression<T> const& v ) : core_(SizeN) {
      VW_ASSERT( v().size()==SizeN, ArgumentErr() << "Vector must have dimension " << SizeN << "." );
      core_.assign(v);
    }

    /// Standard copy assignment operator.
    Vector& operator=( Vector const& v ) {
      core_ = v.core_;
      return *this;
    }

    /// Generalized assignment operator, from arbitrary VW vector expressions.
    template <class T>
    Vector& operator=( VectorBase<T> const& v ) { 
      VW_ASSERT( v.impl().size()==SizeN, ArgumentErr() << "Vector must have dimension " << SizeN << "." );
      std::copy( v.impl().begin(), v.impl().end(), begin() );
      return *this;
    }

    /// Generalized assignment operator, from arbitrary uBLAS vector expressions.
    template <class T>
    Vector& operator=( bnu::vector_expression<T> const& v ) {
      VW_ASSERT( v().size()==SizeN, ArgumentErr() << "Vector must have dimension " << SizeN << "." );
      core_.assign(v);
      return *this;
    }

    /// Returns the size of the vector.
    unsigned size() const { 
      return SizeN;
    }

    /// Change the size of the vector. Elements in memory are preserved when specified.
    void set_size( unsigned new_size, bool preserve = false ) {
      VW_ASSERT( new_size==size(), ArgumentErr() << "Cannot change size of fixed-size Vector." );
    }

    reference_type operator()( unsigned i ) {
      return core_[i];
    }

    const_reference_type operator()( unsigned i ) const {
      return core_[i];
    }

    reference_type operator[]( unsigned i ) {
      return core_[i];
    }

    const_reference_type operator[]( unsigned i ) const {
      return core_[i];
    }

    reference_type x() {
      BOOST_STATIC_ASSERT( SizeN >= 1 );
      return core_[0];
    }

    const_reference_type x() const {
      BOOST_STATIC_ASSERT( SizeN >= 1 );
      return core_[0];
    }

    reference_type y() {
      BOOST_STATIC_ASSERT( SizeN >= 2 );
      return core_[1];
    }

    const_reference_type y() const {
      BOOST_STATIC_ASSERT( SizeN >= 2 );
      return core_[1];
    }

    reference_type z() {
      BOOST_STATIC_ASSERT( SizeN >= 3 );
      return core_[2];
    }

    const_reference_type z() const {
      BOOST_STATIC_ASSERT( SizeN >= 3 );
      return core_[2];
    }

    iterator begin() { 
      return &(core_(0));
    }

    const_iterator begin() const { 
      return &(core_(0));
    }

    iterator end() {
      return &(core_(0)) + SizeN;
    }

    const_iterator end() const {
      return &(core_(0)) + SizeN;
    }

  };

  template <class ElemT, int SizeN>
  struct VectorSize<Vector<ElemT,SizeN> > {
    const static int value = SizeN;
  };
  

  // *******************************************************************
  // class Vector<ElemT>
  // A dynamically-allocated arbitrary-dimension vector class.
  // *******************************************************************

  /// An arbitrary-dimension mathematical vector class.
  template <class ElemT>
  class Vector<ElemT,0> : public VectorBase<Vector<ElemT> > {
    typedef bnu::vector<ElemT> core_type;
    core_type core_;
  public:
    typedef ElemT value_type;

    typedef ElemT& reference_type;
    typedef ElemT const& const_reference_type;

    typedef ElemT* iterator;
    typedef const ElemT* const_iterator;

    /// Constructs a vector with zero size.
    Vector() {}

    /// Constructs a zero vector of the given size.
    Vector( unsigned size ) : core_(size) {
      core_.clear();
    }

    /// Constructs a vector of the given size from given
    /// densely- packed data.  This constructor copies the data.
    /// If you wish to make a shallow proxy object instead, see
    /// vw::VectorProxy.
    Vector( unsigned size, const ElemT *data ) : core_(size) {
      std::copy( data, data+size, begin() );
    }

    /// Standard copy constructor.
    Vector( Vector const& v ) : core_( v.core_ ) {}

    /// Generalized copy constructor, from arbitrary VW vector expressions.
    template <class T>
    Vector( VectorBase<T> const& v ) {
      set_size( v.impl().size() );
      std::copy( v.impl().begin(), v.impl().end(), begin() );
    }

    /// Generalized copy constructor, from arbitrary uBLAS vector expressions.
    template <class T>
    Vector( bnu::vector_expression<T> const& v ) {
      set_size( v().size() );
      core_.assign(v);
    }

    /// Standard copy assignment operator.
    Vector& operator=( Vector const& v ) {
      core_ = v.core_;
      return *this;
    }

    /// Generalized assignment operator, from arbitrary VW vector expressions.
    template <class T>
    Vector& operator=( VectorBase<T> const& v ) {
      set_size( v.impl().size() );
      std::copy( v.impl().begin(), v.impl().end(), begin() );
      return *this;
    }

    /// Generalized assignment operator, from arbitrary uBLAS vector expressions.
    template <class T>
    Vector& operator=( bnu::vector_expression<T> const& v ) {
      set_size( v().size() );
      core_.assign(v);
      return *this;
    }

    /// Returns the size of the vector.
    unsigned size() const { 
      return core_.size();
    }

    /// Change the size of the vector. Elements in memory are preserved when specified.
    void set_size( unsigned new_size, bool preserve = false ) {
      core_.resize(new_size, preserve);
    }

    reference_type operator()( unsigned i ) {
      return core_[i];
    }

    const_reference_type operator()( unsigned i ) const {
      return core_[i];
    }

    reference_type operator[]( unsigned i ) {
      return core_[i];
    }

    const_reference_type operator[]( unsigned i ) const {
      return core_[i];
    }

    iterator begin() { 
      return &(core_(0));
    }

    const_iterator begin() const { 
      return &(core_(0));
    }

    iterator end() {
      return &(core_(0)) + size();
    }

    const_iterator end() const {
      return &(core_(0)) + size();
    }

  };


  // *******************************************************************
  // class VectorProxy<ElemT,SizeN>
  // A fixed-dimension vector proxy class, treating an arbitrary block 
  // of memory as a Vector.
  // *******************************************************************

  /// A fixed-dimension mathematical vector class.
  template <class ElemT, int SizeN = 0>
  class VectorProxy : public VectorBase<VectorProxy<ElemT,SizeN> >
  {
    ElemT *m_ptr;
  public:
    typedef ElemT value_type;

    typedef ElemT& reference_type;
    typedef ElemT const& const_reference_type;

    typedef ElemT* iterator;
    typedef const ElemT* const_iterator;

    /// Constructs a vector proxy.
    VectorProxy( ElemT *ptr ) : m_ptr(ptr) {}

    /// Standard copy assignment operator.
    VectorProxy& operator=( VectorProxy const& v ) {
      VW_ASSERT( v.size()==size(), 
                 ArgumentErr() << "Vector must have dimension " 
                 << size() << " in vector proxy assignment." );
      std::copy( v.begin(), v.end(), begin() );
      return *this;
    }

    /// Generalized assignment operator, from arbitrary VW vector expressions.
    template <class T>
    VectorProxy& operator=( VectorBase<T> const& v ) { 
      VW_ASSERT( v.impl().size()==size(), 
                 ArgumentErr() << "Vector must have dimension " 
                 << size() << " in vector proxy assignment." );
      std::copy( v.impl().begin(), v.impl().end(), begin() );
      return *this;
    }

    /// Generalized assignment operator, from arbitrary uBLAS vector expressions.
    template <class T>
    VectorProxy& operator=( bnu::vector_expression<T> const& v ) {
      VW_ASSERT( v().size()==size(), 
                 ArgumentErr() << "Vector must have dimension " 
                 << size() << " in vector proxy assignment." );
      std::copy( v().begin(), v().end(), begin() );
      return *this;
    }

    /// Returns the size of the vector.
    unsigned size() const { 
      return SizeN;
    }

    /// Change the size of the vector. Elements in memory are preserved when specified.
    void set_size( unsigned new_size, bool preserve = false ) {
      VW_ASSERT( new_size==size(), ArgumentErr() << "Cannot resize a vector proxy." );
    }

    reference_type operator()( unsigned i ) {
      return m_ptr[i];
    }

    const_reference_type operator()( unsigned i ) const {
      return m_ptr[i];
    }

    reference_type operator[]( unsigned i ) {
      return m_ptr[i];
    }

    const_reference_type operator[]( unsigned i ) const {
      return m_ptr[i];
    }

    reference_type x() {
      BOOST_STATIC_ASSERT( SizeN >= 1 );
      return m_ptr[0];
    }

    const_reference_type x() const {
      BOOST_STATIC_ASSERT( SizeN >= 1 );
      return m_ptr[0];
    }

    reference_type y() {
      BOOST_STATIC_ASSERT( SizeN >= 2 );
      return m_ptr[1];
    }

    const_reference_type y() const {
      BOOST_STATIC_ASSERT( SizeN >= 2 );
      return m_ptr[1];
    }

    reference_type z() {
      BOOST_STATIC_ASSERT( SizeN >= 3 );
      return m_ptr[2];
    }

    const_reference_type z() const {
      BOOST_STATIC_ASSERT( SizeN >= 3 );
      return m_ptr[2];
    }

    iterator begin() { 
      return m_ptr;
    }

    const_iterator begin() const { 
      return m_ptr;
    }

    iterator end() {
      return m_ptr + size();
    }

    const_iterator end() const {
      return m_ptr + size();
    }

  };

  template <class ElemT, int SizeN>
  struct VectorSize<VectorProxy<ElemT,SizeN> > {
    const static int value = SizeN;
  };


  // *******************************************************************
  // class VectorProxy<ElemT>
  // An arbitrary-dimension vector proxy class, treating an arbitrary 
  // block of memory as a Vector.
  // *******************************************************************

  /// An arbitrary-dimension vector proxy class.
  template <class ElemT>
  class VectorProxy<ElemT,0> : public VectorBase<VectorProxy<ElemT> > {
    ElemT *m_ptr;
    unsigned m_size;
  public:
    typedef ElemT value_type;

    typedef ElemT& reference_type;
    typedef ElemT const& const_reference_type;

    typedef ElemT* iterator;
    typedef const ElemT* const_iterator;

    /// Constructs a vector with zero size.
    VectorProxy( unsigned size, ElemT *ptr ) : m_ptr(ptr), m_size(size) {}

    VectorProxy& operator=( VectorProxy const& v ) {
      VW_ASSERT( v.size()==size(), 
                 ArgumentErr() << "Vector must have dimension " 
                 << size() << " in vector proxy assignment." );
      std::copy( v.begin(), v.end(), begin() );
      return *this;
    }

    /// Generalized assignment operator, from arbitrary VW vector expressions.
    template <class T>
    VectorProxy& operator=( VectorBase<T> const& v ) { 
      VW_ASSERT( v.impl().size()==size(), 
                 ArgumentErr() << "Vector must have dimension " 
                 << size() << " in vector proxy assignment." );
      std::copy( v.impl().begin(), v.impl().end(), begin() );
      return *this;
    }

    /// Generalized assignment operator, from arbitrary uBLAS vector expressions.
    template <class T>
    VectorProxy& operator=( bnu::vector_expression<T> const& v ) {
      VW_ASSERT( v().size()==size(), 
                 ArgumentErr() << "Vector must have dimension " 
                 << size() << " in vector proxy assignment." );
      std::copy( v().begin(), v().end(), begin() );
      return *this;
    }

    /// Returns the size of the vector.
    unsigned size() const { 
      return m_size;
    }

    /// Change the size of the vector. Elements in memory are preserved when specified.
    void set_size( unsigned new_size, bool preserve = false ) {
      VW_ASSERT( new_size==size(), ArgumentErr() << "Cannot resize a vector proxy." );
    }

    reference_type operator()( unsigned i ) {
      return m_ptr[i];
    }

    const_reference_type operator()( unsigned i ) const {
      return m_ptr[i];
    }

    reference_type operator[]( unsigned i ) {
      return m_ptr[i];
    }

    const_reference_type operator[]( unsigned i ) const {
      return m_ptr[i];
    }

    iterator begin() { 
      return m_ptr;
    }

    const_iterator begin() const { 
      return m_ptr;
    }

    iterator end() {
      return m_ptr + size();
    }

    const_iterator end() const {
      return m_ptr + size();
    }

  };

  /// Shallow proxy view of an block of memory as a vector.
  template <class DataT>
  VectorProxy<DataT> 
  vector_proxy( DataT* data, int size) {
    return VectorProxy<DataT>( data, size );
  }


  // *******************************************************************
  // class VectorTranspose<VectorT>
  // A transposed vector wrapper class.
  // *******************************************************************

  /// A transposed vector class.  This class represents the transposed
  /// version of a vector.  Note that unlike a transposed matrix, a
  /// transposed vector is not simply another vector.  For this reason
  /// VectorTransposed does not derive from VectorBase<>, since that
  /// would permit a variety of mathematically meaninless expressions
  /// to be created with unexpected results.  An unfortunate side
  /// effect of this is that explicit overloads must be provided for
  /// all operations on transposed vectors.  For the moment only a
  /// limited set of operations are supported, including primarily the
  /// mathematical operators.
  template <class VectorT>
  class VectorTranspose {
    // We want to store Vector objects by reference so that we don't copy 
    // them, but we want to store everything else by value so that we can 
    // return transposed versions of various vector expressions.
    template <class T> struct VectorClosure { typedef T type; };
    template <class ElemT, int SizeN> struct VectorClosure<Vector<ElemT,SizeN> > { typedef Vector<ElemT,SizeN>& type; };
    template <class ElemT, int SizeN> struct VectorClosure<const Vector<ElemT,SizeN> > { typedef Vector<ElemT,SizeN> const& type; };
    typename VectorClosure<VectorT>::type m_vector;
  public:
    typedef typename VectorT::value_type value_type;

    typedef typename boost::mpl::if_<boost::is_const<VectorT>,
                                     typename VectorT::const_reference_type,
                                     typename VectorT::reference_type>::type reference_type;
    typedef typename VectorT::const_reference_type const_reference_type;

    typedef typename boost::mpl::if_<boost::is_const<VectorT>,
                                     typename VectorT::const_iterator,
                                     typename VectorT::iterator>::type iterator;
    typedef typename VectorT::const_iterator const_iterator;

    explicit VectorTranspose( VectorT& v ) : m_vector(v) {}
    
    template <class OtherT>
    VectorTranspose& operator=( VectorTranspose<OtherT> const& v ) {
      VW_ASSERT( v.size()==size(), 
                 ArgumentErr() << "Vectors must have same size in transposed vector assignment" );
      m_vector = v.child();
      return *this;
    }

    VectorT& child() {
      return m_vector;
    }

    VectorT const& child() const {
      return m_vector;
    }

    unsigned size() const {
      return m_vector.size();
    }

    reference_type operator()( int i ) {
      return m_vector(i);
    }

    const_reference_type operator()( int i ) const {
      return m_vector(i);
    }

    iterator begin() {
      return m_vector.begin();
    }

    const_iterator begin() const {
      return m_vector.begin();
    }

    iterator end() {
      return m_vector.begin();
    }

    const_iterator end() const {
      return m_vector.begin();
    }

  };

  /// Vector transpose.
  template <class VectorT>
  inline VectorTranspose<VectorT> transpose( VectorBase<VectorT>& vector ) {
    return VectorTranspose<VectorT>( vector.impl() );
  }

  /// Vector transpose (const overload).
  template <class VectorT>
  inline VectorTranspose<const VectorT> transpose( VectorBase<VectorT> const& vector ) {<commit_msg>Fixed a bug in VectorTranspose.<commit_after>// __BEGIN_LICENSE__
// 
// Copyright (C) 2006 United States Government as represented by the
// Administrator of the National Aeronautics and Space Administration
// (NASA).  All Rights Reserved.
// 
// Copyright 2006 <NAME>. All rights reserved.
// 
// This software is distributed under the NASA Open Source Agreement
// (NOSA), version 1.3.  The NOSA has been approved by the Open Source
// Initiative.  See the file COPYING at the top of the distribution
// directory tree for the complete NOSA document.
// 
// THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY
// KIND, EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT
// LIMITED TO, ANY WARRANTY THAT THE SUBJECT SOFTWARE WILL CONFORM TO
// SPECIFICATIONS, ANY IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR
// A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY WARRANTY THAT
// THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT
// DOCUMENTATION, IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE.
// 
// __END_LICENSE__

/// \file Vector.h
///
/// Provides the core Vector<> mathematical vector template classes.  
///
/// Currently I believe we support:
///   Dynamic (Vector<T>) and fixed (Vector<T,N>) vectors
///   Element access via vec(i) and vec[i]
///   Seamless conversion between compatibile vector types
///   Vector<T,N> construction with between zero and four initializer elements
///   Dynamic and fixed VectorProxy proxy objects
///   Subvector access via subvector()
///   Vector transposition via transpose()
///   Explicit expression evaluation via eval()
///   Printing of vectors to ostreams
///   Equality of vectors with or without an epsilon
///   Vector negation, addition, subtraction
///   Scalar multiplication and division
///   Elementwise vector addition, subtraction, multiplication, and division
///   Elementwise scalar addition, subtraction, multiplication, and division
///   Vector negation, addition, subtraction of transposed vectors
///   Scalar multiplication and division of transposed vectors
///   Elementwise comprison operations.
///   Norms via norm_1(), norm_2(), norm_2_sqr(), and norm_inf()
///   Sum and product of elements via sum() and prod()
///   Vector normalization via normalize()
///   Dot/inner product via dot_prod(v1,v2) or transpose(v1)*v2
///   Cross product of 3-element vectors via cross_prod()
///
#ifndef __VW_MATH_VECTOR_H__
#define __VW_MATH_VECTOR_H__

#include <boost/utility/enable_if.hpp>
#include <boost/type_traits.hpp>
#include <boost/mpl/if.hpp>
#include <boost/static_assert.hpp>
#include <boost/array.hpp>
#include <boost/iterator/iterator_facade.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/utility/result_of.hpp>

#include <vw/Core/Exception.h>
#include <vw/Core/TypeDeduction.h>
#include <vw/Core/Functors.h>

namespace vw {
namespace math {

  namespace bnu = boost::numeric::ublas;

  /// \cond INTERNAL
  // A wrapper around boost::array that provides a constructor taking 
  // a size parameter, which the BOOST uBLAS containers expect to be 
  // available.  With debugging disabled this whole thing should be 
  // optimized out of existence.
  template <class ElemT, int SizeN>
  class FixedArray : public boost::array<ElemT,SizeN> {
  public:
    FixedArray() {}
    FixedArray( unsigned size ) {
      VW_ASSERT( size==SizeN, LogicErr() << "Internal size mismatch in FixedArray." );
    }
  };
  /// \endcond

  /// A type function to compute the dimension of a vector expression 
  /// at compile time (or zero for dynamically-sized vectors).
  template <class VectorT>
  struct VectorSize {
    const static int value = 0;
  };


  // *******************************************************************
  // class VectorBase<VectorT>
  // A CRTP vector base class.
  // *******************************************************************

  /// A CRTP base class for vectors and vector expressions.  
  /// Provides a mechanism for restricting function arguments to 
  /// vectors, and provides the various arithmetic assignment 
  /// operators.
  template <class VectorT>
  struct VectorBase {

    /// Returns the derived implementation type.
    VectorT& impl() { return *static_cast<VectorT*>(this); }

    /// Returns the derived implementation type.
    VectorT const& impl() const { return *static_cast<VectorT const*>(this); }
    
    /// Sum-assignment operator
    template <class T>
    VectorT& operator+=( T const& v ) {
      return impl() = impl() + v;
    }

    /// Difference-assignment operator
    template <class T>
    VectorT& operator-=( T const& v ) {
      return impl() = impl() - v;
    }

    /// Product-assignment operator
    template <class T>
    VectorT& operator*=( T s ) {
      return impl() = impl() * s;
    }

    /// Quotient-assignment operator
    template <class T>
    VectorT& operator/=( T s ) {
      return impl() = impl() / s;
    }

  };


  // *******************************************************************
  // class Vector<ElemT,SizeN>
  // A statically-allocated fixed-dimension vector class.
  // *******************************************************************

  /// A fixed-dimension mathematical vector class.
  template <class ElemT, int SizeN = 0>
  class Vector : public VectorBase<Vector<ElemT,SizeN> >
  {
    typedef bnu::vector<ElemT,FixedArray<ElemT,SizeN> > core_type;
    core_type core_;
  public:
    typedef ElemT value_type;

    typedef ElemT& reference_type;
    typedef ElemT const& const_reference_type;

    typedef ElemT* iterator;
    typedef const ElemT* const_iterator;

    /// Constructs a vector of zeroes.
    Vector() : core_(SizeN) {
      core_.clear();
    }

    /// Constructs a vector whose first element is as given.
    Vector( ElemT e1 ) : core_(SizeN) {
      BOOST_STATIC_ASSERT( SizeN >= 1 );
      (*this)[0] = e1;
      for( unsigned i=1; i<SizeN; ++i ) (*this)[i] = ElemT();
    }

    /// Constructs a vector whose first two elements are as given.
    Vector( ElemT e1, ElemT e2 ) : core_(SizeN) {
      BOOST_STATIC_ASSERT( SizeN >= 2 );
      (*this)[0] = e1; (*this)[1] = e2;
      for( unsigned i=2; i<SizeN; ++i ) (*this)[i] = ElemT();
    }

    /// Constructs a vector whose first three elements are as given.
    Vector( ElemT e1, ElemT e2, ElemT e3 ) : core_(SizeN) {
      BOOST_STATIC_ASSERT( SizeN >= 3 );
      (*this)[0] = e1; (*this)[1] = e2; (*this)[2] = e3;
      for( unsigned i=3; i<SizeN; ++i ) (*this)[i] = ElemT();
    }

    /// Constructs a vector whose first four elements are as given.
    Vector( ElemT e1, ElemT e2, ElemT e3, ElemT e4 ) : core_(SizeN) {
      BOOST_STATIC_ASSERT( SizeN >= 4 );
      (*this)[0] = e1; (*this)[1] = e2; (*this)[2] = e3; (*this)[3] = e4;
      for( unsigned i=4; i<SizeN; ++i ) (*this)[i] = ElemT();
    }

    /// Constructs a vector from given densely-packed data.  This
    /// constructor copies the data.  If you wish to make a shallow
    /// proxy object instead, see vw::VectorProxy.
    Vector( const ElemT data[SizeN] ) : core_(SizeN) {
      std::copy( data, data+SizeN, begin() );
    }

    /// Standard copy constructor.
    Vector( Vector const& v ) : core_( v.core_ ) {}

    /// Generalized copy constructor, from arbitrary VW vector expressions.
    template <class T>
    Vector( VectorBase<T> const& v ) : core_(SizeN) { 
      VW_ASSERT( v.impl().size()==SizeN, ArgumentErr() << "Vector must have dimension " << SizeN << "." );
      std::copy( v.impl().begin(), v.impl().end(), begin() );
    }

    /// Generalized copy constructor, from arbitrary uBLAS vector expressions.
    template <class T>
    Vector( bnu::vector_expression<T> const& v ) : core_(SizeN) {
      VW_ASSERT( v().size()==SizeN, ArgumentErr() << "Vector must have dimension " << SizeN << "." );
      core_.assign(v);
    }

    /// Standard copy assignment operator.
    Vector& operator=( Vector const& v ) {
      core_ = v.core_;
      return *this;
    }

    /// Generalized assignment operator, from arbitrary VW vector expressions.
    template <class T>
    Vector& operator=( VectorBase<T> const& v ) { 
      VW_ASSERT( v.impl().size()==SizeN, ArgumentErr() << "Vector must have dimension " << SizeN << "." );
      std::copy( v.impl().begin(), v.impl().end(), begin() );
      return *this;
    }

    /// Generalized assignment operator, from arbitrary uBLAS vector expressions.
    template <class T>
    Vector& operator=( bnu::vector_expression<T> const& v ) {
      VW_ASSERT( v().size()==SizeN, ArgumentErr() << "Vector must have dimension " << SizeN << "." );
      core_.assign(v);
      return *this;
    }

    /// Returns the size of the vector.
    unsigned size() const { 
      return SizeN;
    }

    /// Change the size of the vector. Elements in memory are preserved when specified.
    void set_size( unsigned new_size, bool preserve = false ) {
      VW_ASSERT( new_size==size(), ArgumentErr() << "Cannot change size of fixed-size Vector." );
    }

    reference_type operator()( unsigned i ) {
      return core_[i];
    }

    const_reference_type operator()( unsigned i ) const {
      return core_[i];
    }

    reference_type operator[]( unsigned i ) {
      return core_[i];
    }

    const_reference_type operator[]( unsigned i ) const {
      return core_[i];
    }

    reference_type x() {
      BOOST_STATIC_ASSERT( SizeN >= 1 );
      return core_[0];
    }

    const_reference_type x() const {
      BOOST_STATIC_ASSERT( SizeN >= 1 );
      return core_[0];
    }

    reference_type y() {
      BOOST_STATIC_ASSERT( SizeN >= 2 );
      return core_[1];
    }

    const_reference_type y() const {
      BOOST_STATIC_ASSERT( SizeN >= 2 );
      return core_[1];
    }

    reference_type z() {
      BOOST_STATIC_ASSERT( SizeN >= 3 );
      return core_[2];
    }

    const_reference_type z() const {
      BOOST_STATIC_ASSERT( SizeN >= 3 );
      return core_[2];
    }

    iterator begin() { 
      return &(core_(0));
    }

    const_iterator begin() const { 
      return &(core_(0));
    }

    iterator end() {
      return &(core_(0)) + SizeN;
    }

    const_iterator end() const {
      return &(core_(0)) + SizeN;
    }

  };

  template <class ElemT, int SizeN>
  struct VectorSize<Vector<ElemT,SizeN> > {
    const static int value = SizeN;
  };
  

  // *******************************************************************
  // class Vector<ElemT>
  // A dynamically-allocated arbitrary-dimension vector class.
  // *******************************************************************

  /// An arbitrary-dimension mathematical vector class.
  template <class ElemT>
  class Vector<ElemT,0> : public VectorBase<Vector<ElemT> > {
    typedef bnu::vector<ElemT> core_type;
    core_type core_;
  public:
    typedef ElemT value_type;

    typedef ElemT& reference_type;
    typedef ElemT const& const_reference_type;

    typedef ElemT* iterator;
    typedef const ElemT* const_iterator;

    /// Constructs a vector with zero size.
    Vector() {}

    /// Constructs a zero vector of the given size.
    Vector( unsigned size ) : core_(size) {
      core_.clear();
    }

    /// Constructs a vector of the given size from given
    /// densely- packed data.  This constructor copies the data.
    /// If you wish to make a shallow proxy object instead, see
    /// vw::VectorProxy.
    Vector( unsigned size, const ElemT *data ) : core_(size) {
      std::copy( data, data+size, begin() );
    }

    /// Standard copy constructor.
    Vector( Vector const& v ) : core_( v.core_ ) {}

    /// Generalized copy constructor, from arbitrary VW vector expressions.
    template <class T>
    Vector( VectorBase<T> const& v ) {
      set_size( v.impl().size() );
      std::copy( v.impl().begin(), v.impl().end(), begin() );
    }

    /// Generalized copy constructor, from arbitrary uBLAS vector expressions.
    template <class T>
    Vector( bnu::vector_expression<T> const& v ) {
      set_size( v().size() );
      core_.assign(v);
    }

    /// Standard copy assignment operator.
    Vector& operator=( Vector const& v ) {
      core_ = v.core_;
      return *this;
    }

    /// Generalized assignment operator, from arbitrary VW vector expressions.
    template <class T>
    Vector& operator=( VectorBase<T> const& v ) {
      set_size( v.impl().size() );
      std::copy( v.impl().begin(), v.impl().end(), begin() );
      return *this;
    }

    /// Generalized assignment operator, from arbitrary uBLAS vector expressions.
    template <class T>
    Vector& operator=( bnu::vector_expression<T> const& v ) {
      set_size( v().size() );
      core_.assign(v);
      return *this;
    }

    /// Returns the size of the vector.
    unsigned size() const { 
      return core_.size();
    }

    /// Change the size of the vector. Elements in memory are preserved when specified.
    void set_size( unsigned new_size, bool preserve = false ) {
      core_.resize(new_size, preserve);
    }

    reference_type operator()( unsigned i ) {
      return core_[i];
    }

    const_reference_type operator()( unsigned i ) const {
      return core_[i];
    }

    reference_type operator[]( unsigned i ) {
      return core_[i];
    }

    const_reference_type operator[]( unsigned i ) const {
      return core_[i];
    }

    iterator begin() { 
      return &(core_(0));
    }

    const_iterator begin() const { 
      return &(core_(0));
    }

    iterator end() {
      return &(core_(0)) + size();
    }

    const_iterator end() const {
      return &(core_(0)) + size();
    }

  };


  // *******************************************************************
  // class VectorProxy<ElemT,SizeN>
  // A fixed-dimension vector proxy class, treating an arbitrary block 
  // of memory as a Vector.
  // *******************************************************************

  /// A fixed-dimension mathematical vector class.
  template <class ElemT, int SizeN = 0>
  class VectorProxy : public VectorBase<VectorProxy<ElemT,SizeN> >
  {
    ElemT *m_ptr;
  public:
    typedef ElemT value_type;

    typedef ElemT& reference_type;
    typedef ElemT const& const_reference_type;

    typedef ElemT* iterator;
    typedef const ElemT* const_iterator;

    /// Constructs a vector proxy.
    VectorProxy( ElemT *ptr ) : m_ptr(ptr) {}

    /// Standard copy assignment operator.
    VectorProxy& operator=( VectorProxy const& v ) {
      VW_ASSERT( v.size()==size(), 
                 ArgumentErr() << "Vector must have dimension " 
                 << size() << " in vector proxy assignment." );
      std::copy( v.begin(), v.end(), begin() );
      return *this;
    }

    /// Generalized assignment operator, from arbitrary VW vector expressions.
    template <class T>
    VectorProxy& operator=( VectorBase<T> const& v ) { 
      VW_ASSERT( v.impl().size()==size(), 
                 ArgumentErr() << "Vector must have dimension " 
                 << size() << " in vector proxy assignment." );
      std::copy( v.impl().begin(), v.impl().end(), begin() );
      return *this;
    }

    /// Generalized assignment operator, from arbitrary uBLAS vector expressions.
    template <class T>
    VectorProxy& operator=( bnu::vector_expression<T> const& v ) {
      VW_ASSERT( v().size()==size(), 
                 ArgumentErr() << "Vector must have dimension " 
                 << size() << " in vector proxy assignment." );
      std::copy( v().begin(), v().end(), begin() );
      return *this;
    }

    /// Returns the size of the vector.
    unsigned size() const { 
      return SizeN;
    }

    /// Change the size of the vector. Elements in memory are preserved when specified.
    void set_size( unsigned new_size, bool preserve = false ) {
      VW_ASSERT( new_size==size(), ArgumentErr() << "Cannot resize a vector proxy." );
    }

    reference_type operator()( unsigned i ) {
      return m_ptr[i];
    }

    const_reference_type operator()( unsigned i ) const {
      return m_ptr[i];
    }

    reference_type operator[]( unsigned i ) {
      return m_ptr[i];
    }

    const_reference_type operator[]( unsigned i ) const {
      return m_ptr[i];
    }

    reference_type x() {
      BOOST_STATIC_ASSERT( SizeN >= 1 );
      return m_ptr[0];
    }

    const_reference_type x() const {
      BOOST_STATIC_ASSERT( SizeN >= 1 );
      return m_ptr[0];
    }

    reference_type y() {
      BOOST_STATIC_ASSERT( SizeN >= 2 );
      return m_ptr[1];
    }

    const_reference_type y() const {
      BOOST_STATIC_ASSERT( SizeN >= 2 );
      return m_ptr[1];
    }

    reference_type z() {
      BOOST_STATIC_ASSERT( SizeN >= 3 );
      return m_ptr[2];
    }

    const_reference_type z() const {
      BOOST_STATIC_ASSERT( SizeN >= 3 );
      return m_ptr[2];
    }

    iterator begin() { 
      return m_ptr;
    }

    const_iterator begin() const { 
      return m_ptr;
    }

    iterator end() {
      return m_ptr + size();
    }

    const_iterator end() const {
      return m_ptr + size();
    }

  };

  template <class ElemT, int SizeN>
  struct VectorSize<VectorProxy<ElemT,SizeN> > {
    const static int value = SizeN;
  };


  // *******************************************************************
  // class VectorProxy<ElemT>
  // An arbitrary-dimension vector proxy class, treating an arbitrary 
  // block of memory as a Vector.
  // *******************************************************************

  /// An arbitrary-dimension vector proxy class.
  template <class ElemT>
  class VectorProxy<ElemT,0> : public VectorBase<VectorProxy<ElemT> > {
    ElemT *m_ptr;
    unsigned m_size;
  public:
    typedef ElemT value_type;

    typedef ElemT& reference_type;
    typedef ElemT const& const_reference_type;

    typedef ElemT* iterator;
    typedef const ElemT* const_iterator;

    /// Constructs a vector with zero size.
    VectorProxy( unsigned size, ElemT *ptr ) : m_ptr(ptr), m_size(size) {}

    VectorProxy& operator=( VectorProxy const& v ) {
      VW_ASSERT( v.size()==size(), 
                 ArgumentErr() << "Vector must have dimension " 
                 << size() << " in vector proxy assignment." );
      std::copy( v.begin(), v.end(), begin() );
      return *this;
    }

    /// Generalized assignment operator, from arbitrary VW vector expressions.
    template <class T>
    VectorProxy& operator=( VectorBase<T> const& v ) { 
      VW_ASSERT( v.impl().size()==size(), 
                 ArgumentErr() << "Vector must have dimension " 
                 << size() << " in vector proxy assignment." );
      std::copy( v.impl().begin(), v.impl().end(), begin() );
      return *this;
    }

    /// Generalized assignment operator, from arbitrary uBLAS vector expressions.
    template <class T>
    VectorProxy& operator=( bnu::vector_expression<T> const& v ) {
      VW_ASSERT( v().size()==size(), 
                 ArgumentErr() << "Vector must have dimension " 
                 << size() << " in vector proxy assignment." );
      std::copy( v().begin(), v().end(), begin() );
      return *this;
    }

    /// Returns the size of the vector.
    unsigned size() const { 
      return m_size;
    }

    /// Change the size of the vector. Elements in memory are preserved when specified.
    void set_size( unsigned new_size, bool preserve = false ) {
      VW_ASSERT( new_size==size(), ArgumentErr() << "Cannot resize a vector proxy." );
    }

    reference_type operator()( unsigned i ) {
      return m_ptr[i];
    }

    const_reference_type operator()( unsigned i ) const {
      return m_ptr[i];
    }

    reference_type operator[]( unsigned i ) {
      return m_ptr[i];
    }

    const_reference_type operator[]( unsigned i ) const {
      return m_ptr[i];
    }

    iterator begin() { 
      return m_ptr;
    }

    const_iterator begin() const { 
      return m_ptr;
    }

    iterator end() {
      return m_ptr + size();
    }

    const_iterator end() const {
      return m_ptr + size();
    }

  };

  /// Shallow proxy view of an block of memory as a vector.
  template <class DataT>
  VectorProxy<DataT> 
  vector_proxy( DataT* data, int size) {
    return VectorProxy<DataT>( data, size );
  }


  // *******************************************************************
  // class VectorTranspose<VectorT>
  // A transposed vector wrapper class.
  // *******************************************************************

  /// A transposed vector class.  This class represents the transposed
  /// version of a vector.  Note that unlike a transposed matrix, a
  /// transposed vector is not simply another vector.  For this reason
  /// VectorTransposed does not derive from VectorBase<>, since that
  /// would permit a variety of mathematically meaninless expressions
  /// to be created with unexpected results.  An unfortunate side
  /// effect of this is that explicit overloads must be provided for
  /// all operations on transposed vectors.  For the moment only a
  /// limited set of operations are supported, including primarily the
  /// mathematical operators.
  template <class VectorT>
  class VectorTranspose {
    // We want to store Vector objects by reference so that we don't copy 
    // them, but we want to store everything else by value so that we can 
    // return transposed versions of various vector expressions.
    template <class T> struct VectorClosure { typedef T type; };
    template <class ElemT, int SizeN> struct VectorClosure<Vector<ElemT,SizeN> > { typedef Vector<ElemT,SizeN>& type; };
    template <class ElemT, int SizeN> struct VectorClosure<const Vector<ElemT,SizeN> > { typedef Vector<ElemT,SizeN> const& type; };
    typename VectorClosure<VectorT>::type m_vector;
  public:
    typedef typename VectorT::value_type value_type;

    typedef typename boost::mpl::if_<boost::is_const<VectorT>,
                                     typename VectorT::const_reference_type,
                                     typename VectorT::reference_type>::type reference_type;
    typedef typename VectorT::const_reference_type const_reference_type;

    typedef typename boost::mpl::if_<boost::is_const<VectorT>,
                                     typename VectorT::const_iterator,
                                     typename VectorT::iterator>::type iterator;
    typedef typename VectorT::const_iterator const_iterator;

    explicit VectorTranspose( VectorT& v ) : m_vector(v) {}
    
    template <class OtherT>
    VectorTranspose& operator=( VectorTranspose<OtherT> const& v ) {
      VW_ASSERT( v.size()==size(), 
                 ArgumentErr() << "Vectors must have same size in transposed vector assignment" );
      m_vector = v.child();
      return *this;
    }

    VectorT& child() {
      return m_vector;
    }

    VectorT const& child() const {
      return m_vector;
    }

    unsigned size() const {
      return m_vector.size();
    }

    reference_type operator()( int i ) {
      return m_vector(i);
    }

    const_reference_type operator()( int i ) const {
      return m_vector(i);
    }

    iterator begin() {
      return m_vector.begin();
    }

    const_iterator begin() const {
      return m_vector.begin();
    }

    iterator end() {
      return m_vector.end();
    }

    const_iterator end() const {
      return m_vector.end();
    }

  };

  /// Vector transpose.
  template <class VectorT>
  inline VectorTranspose<VectorT> transpose( VectorBase<VectorT>& vector ) {
    return VectorTranspose<VectorT>( vector.impl() );
  }

  /// Vector transpose (const overload).
  template <class VectorT>
  inline VectorTranspose<const VectorT> transpose( VectorBase<VectorT> const& vector ) {
{'fraction_non_alphanumeric': 0.15168207336119202, 'fraction_numerical': 0.003925541344814486, 'mean_word_length': 2.7563817980022196, 'pattern_counts': {'":': 0, '<': 367, '<?xml version=': 0, '>': 251, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1737726', 'n_tokens_mistral': 14105, 'n_tokens_neox': 13390, 'n_words': 4314}