VectorD.h

Go to the documentation of this file.

/**
 *  \file IMP/algebra/VectorD.h   \brief Simple D vector class.
 *
 *  Copyright 2007-2013 IMP Inventors. All rights reserved.
 *
 */

#ifndef IMPALGEBRA_VECTOR_D_H
#define IMPALGEBRA_VECTOR_D_H

#include <IMP/algebra/algebra_config.h>
#include "GeometricPrimitiveD.h"
#include <IMP/base/types.h>
#include <IMP/base/check_macros.h>
#include <IMP/base/exception.h>
#include <IMP/base/utility.h>
#include <IMP/base/InputAdaptor.h>
#include <IMP/base/random.h>
#include "internal/vector.h"

#include <limits>
#include <cmath>
#include <boost/random/normal_distribution.hpp>
#include <boost/static_assert.hpp>

#if IMP_HAS_CHECKS >= IMP_USAGE
#define IMP_VECTOR_CHECK check_vector()
#define IMP_VECTOR_CHECK_INDEX(i) check_index(i)
#define IMP_VECTOR_CHECK_COMPATIBLE(o) \
  check_compatible_vector(o);          \
  o.check_vector()
#else
#define IMP_VECTOR_CHECK
#define IMP_VECTOR_CHECK_INDEX(i)
#define IMP_VECTOR_CHECK_COMPATIBLE(o)
#endif

IMPALGEBRA_BEGIN_NAMESPACE
//! A Cartesian vector in D-dimensions.
/** Store a vector of Cartesian coordinates. It supports all expected
    mathematical operators, including using * for the dot product.
    \see Vector3D
    \see Vector2D

    \geometry
 */
template <int D>
class VectorD : public GeometricPrimitiveD<D> {
  void check_vector() const {
    IMP_USAGE_CHECK(!data_.get_is_null(),
                    "Attempt to use uninitialized vector.");
  }
  template <int OD>
  void check_compatible_vector(const VectorD<OD> &o) const {
    IMP_USAGE_CHECK_VARIABLE(o);
    IMP_USAGE_CHECK(o.get_dimension() == get_dimension(),
                    "Dimensions don't match: " << get_dimension() << " vs "
                                               << o.get_dimension());
  }
  void check_index(unsigned int i) const {
#if IMP_HAS_CHECKS < IMP_INTERNAL
    IMP_UNUSED(i);
#endif
    IMP_INTERNAL_CHECK(i < data_.get_dimension(),
                       "Invalid component of vector requested: "
                           << i << " of " << get_dimension());
  }

 public:
#if !defined(SWIG) && !defined(IMP_DOXYGEN)
  // to shut off deprecation warning
  VectorD(const Floats &f, bool) {
    if (D != -1 && static_cast<int>(f.size()) != D) {
      IMP_THROW("Expected " << D << " but got " << f.size(),
                base::ValueException);
    }
    data_.set_coordinates(f.begin(), f.end());
  }
  template <int OD>
  VectorD(const VectorD<OD> &o) {
    BOOST_STATIC_ASSERT(D == -1 || OD == -1 || D == OD);
    IMP_USAGE_CHECK(
        D == -1 || o.get_dimension() == static_cast<unsigned int>(D),
        "Dimensions don't match in conversion");
    data_.set_coordinates(o.coordinates_begin(), o.coordinates_end());
  }
  template <int OD>
  VectorD &operator=(const VectorD<OD> &o) {
    BOOST_STATIC_ASSERT(D == -1 || OD == -1 || D == OD);
    IMP_USAGE_CHECK(D == -1 || o.get_dimension() == D,
                    "Dimensions don't match in conversion");
    data_.set_coordinates(o.coordinates_begin(), o.coordinates_end());
  }
#endif

/** \throw base::ValueException if f.size() is not appropriate.
    \note Only use this from python. */
#ifndef SWIG
  IMP_DEPRECATED_ATTRIBUTE
#endif
  VectorD(const Floats &f) {
    if (D != -1 && static_cast<int>(f.size()) != D) {
      IMP_THROW("Expected " << D << " but got " << f.size(),
                base::ValueException);
    }
    data_.set_coordinates(f.begin(), f.end());
  }

  /** The distance between b and e must be equal to D.
   */
  template <class It>
  VectorD(It b, It e) {
    data_.set_coordinates(b, e);
  }

  //! Initialize the 1-vector from its value.
  explicit VectorD(double x) {
/* Note that MSVC gets confused with static asserts if we try to subclass
   VectorD, as we do for example in the various IMP::display Geometry
   subclasses, so replace with runtime checks. */
#if defined(IMP_SWIG_WRAPPER) || defined(_MSC_VER)
    IMP_USAGE_CHECK(D == 1 || D == -1, "Need " << D << " to construct a " << D
                                               << "-vector.");
#else
    BOOST_STATIC_ASSERT(D == 1);
#endif
    data_.set_coordinates(&x, &x + 1);
  }

  //! Initialize a 2-vector from separate x,y values.
  VectorD(double x, double y) {
#if defined(IMP_SWIG_WRAPPER) || defined(_MSC_VER)
    IMP_USAGE_CHECK(D == 2 || D == -1, "Need " << D << " to construct a " << D
                                               << "-vector.");
#else
    BOOST_STATIC_ASSERT(D == 2);
#endif
    double d[] = {x, y};
    data_.set_coordinates(d, d + 2);
  }

  //! Initialize a 3-vector from separate x,y,z values.
  VectorD(double x, double y, double z) {
#ifdef IMP_SWIG_WRAPPER
    IMP_USAGE_CHECK(D == 3 || D == -1, "Need " << D << " to construct a " << D
                                               << "-vector.");
#else
    BOOST_STATIC_ASSERT(D == 3);
#endif
    double d[] = {x, y, z};
    data_.set_coordinates(d, d + 3);
  }

  //! Initialize a 4-vector from separate w,x,y,z values.
  VectorD(double x0, double x1, double x2, double x3) {
#if defined(IMP_SWIG_WRAPPER) || defined(_MSC_VER)
    IMP_USAGE_CHECK(D == 4 || D == -1, "Need " << D << " to construct a " << D
                                               << "-vector.");
#else
    BOOST_STATIC_ASSERT(D == 4);
#endif
    double d[] = {x0, x1, x2, x3};
    data_.set_coordinates(d, d + 4);
  }

  //! Default constructor
  VectorD() {}

  /** Return the ith Cartesian coordinate. In 3D use [0] to get
      the x coordinate etc.*/
  inline double operator[](unsigned int i) const {
    IMP_VECTOR_CHECK_INDEX(i);
    IMP_VECTOR_CHECK;
    return data_.get_data()[i];
  }
  /** Return the ith Cartesian coordinate. In 3D use [0] to get
      the x coordinate etc. */
  inline double &operator[](unsigned int i) {
    IMP_VECTOR_CHECK_INDEX(i);
    return data_.get_data()[i];
  }

  double get_scalar_product(const VectorD<D> &o) const {
    IMP_VECTOR_CHECK_COMPATIBLE(o);
    IMP_VECTOR_CHECK;
    double ret = 0;
    for (unsigned int i = 0; i < get_dimension(); ++i) {
      ret += operator[](i) * o.operator[](i);
    }
    return ret;
  }

  double get_squared_magnitude() const { return get_scalar_product(*this); }

  double get_magnitude() const { return std::sqrt(get_squared_magnitude()); }

  /**
     Returns a unit vector pointing at the same direction as this vector.

     @note If the magnitude of this vector is smaller than 1e-12
           (an arbitrarily selected small number), returns a unit
           vector pointing at a random direction
   */
  VectorD get_unit_vector() const {
    const double tiny_double = 1e-12;
    double mag = get_magnitude();
    if (mag > tiny_double) {
      return operator/(mag);
    } else {
      // avoid division by zero - return random unit v
      // NOTE: (1) avoids vector_generators / SphereD to prevent recursiveness
      //       (2) D might be -1, so use get_dimension()
      VectorD<D> ret(*this);
      boost::variate_generator<boost::rand48, boost::normal_distribution<> >
          generator(IMP::base::random_number_generator,
                    ::boost::normal_distribution<>(0, 1.0));
      for (unsigned int i = 0; i < get_dimension(); ++i) {
        ret[i] = generator();
      }
      return ret.get_unit_vector();
    }
  }

#ifndef IMP_DOXYGEN
  double operator*(const VectorD<D> &o) const {
    IMP_VECTOR_CHECK_COMPATIBLE(o);
    return get_scalar_product(o);
  }

  VectorD operator*(double s) const {
    IMP_VECTOR_CHECK;
    VectorD ret = *this;
    ret *= s;
    return ret;
  }

  VectorD operator/(double s) const {
    IMP_VECTOR_CHECK;
    VectorD ret = *this;
    ret /= s;
    return ret;
  }

  VectorD operator-() const {
    IMP_VECTOR_CHECK;
    VectorD ret = *this;
    for (unsigned int i = 0; i < get_dimension(); ++i) {
      ret[i] = -ret[i];
    }
    return ret;
  }

  VectorD operator-(const VectorD &o) const {
    IMP_VECTOR_CHECK_COMPATIBLE(o);
    IMP_VECTOR_CHECK;
    VectorD ret = *this;
    ret -= o;
    return ret;
  }

  VectorD operator+(const VectorD &o) const {
    IMP_VECTOR_CHECK_COMPATIBLE(o);
    IMP_VECTOR_CHECK;
    VectorD ret = *this;
    ret += o;
    return ret;
  }

  VectorD &operator+=(const VectorD &o) {
    IMP_VECTOR_CHECK_COMPATIBLE(o);
    IMP_VECTOR_CHECK;
    for (unsigned int i = 0; i < get_dimension(); ++i) {
      operator[](i) += o[i];
    }
    return *this;
  }

  VectorD &operator-=(const VectorD &o) {
    IMP_VECTOR_CHECK_COMPATIBLE(o);
    IMP_VECTOR_CHECK;
    for (unsigned int i = 0; i < get_dimension(); ++i) {
      operator[](i) -= o[i];
    }
    return *this;
  }

  VectorD &operator/=(double f) {
    IMP_VECTOR_CHECK;
    for (unsigned int i = 0; i < get_dimension(); ++i) {
      operator[](i) /= f;
    }
    return *this;
  }

  VectorD &operator*=(double f) {
    IMP_VECTOR_CHECK;
    for (unsigned int i = 0; i < get_dimension(); ++i) {
      operator[](i) *= f;
    }
    return *this;
  }

  void show(std::ostream &out, std::string delim, bool parens = true) const {
    IMP_VECTOR_CHECK;
    if (parens) out << "(";
    for (unsigned int i = 0; i < get_dimension(); ++i) {
      out << operator[](i);
      if (i != get_dimension() - 1) {
        out << delim;
      }
    }
    if (parens) out << ")";
  }
  IMP_SHOWABLE_INLINE(VectorD, show(out, ", "););
#endif

#ifndef SWIG
  typedef double *CoordinateIterator;
  CoordinateIterator coordinates_begin() { return data_.get_data(); }
  CoordinateIterator coordinates_end() {
    return data_.get_data() + get_dimension();
  }
  typedef const double *CoordinateConstIterator;
  CoordinateConstIterator coordinates_begin() const { return data_.get_data(); }
  CoordinateConstIterator coordinates_end() const {
    return data_.get_data() + get_dimension();
  }
#endif

#ifndef SWIG
  // For some reason, this method breaks IMP::atom::get_rmsd() in Python, so
  // hide it from SWIG
  Floats get_coordinates() const {
    return Floats(coordinates_begin(), coordinates_end());
  }
#endif

#if !defined(IMP_DOXYGEN) && !defined(SWIG)
  const double *get_data() const { return data_.get_data(); }
#endif
  unsigned int get_dimension() const { return data_.get_dimension(); }

 private:

  internal::VectorData<double, D, false> data_;
};

#ifndef IMP_DOXYGEN

template <int D>
inline std::ostream &operator<<(std::ostream &out, const VectorD<D> &v) {
  v.show(out);
  return out;
}

template <int D>
inline std::istream &operator>>(std::istream &in, VectorD<D> &v) {
  for (unsigned int i = 0; i < D; ++i) {
    in >> v[i];
  }
  return in;
}

#endif

//! lexicographic comparison of two vectors
/** Note that this is not very reliable and probably should not be used.
    See VectorD
 */
template <int D>
inline int compare(const VectorD<D> &a, const VectorD<D> &b) {
  IMP_USAGE_CHECK(a.get_dimension() == b.get_dimension(),
                  "Dimensions don't match.");
  for (unsigned int i = 0; i < a.get_dimension(); ++i) {
    if (a[i] < b[i])
      return -1;
    else if (a[i] > b[i])
      return 1;
  }
  return 0;
}

/** See VectorD */
template <int D>
inline VectorD<D> operator*(double s, const VectorD<D> &o) {
  return o * s;
}

//! compute the squared distance between two vectors
/** See VectorD
 */
template <int D>
inline double get_squared_distance(const VectorD<D> &v1, const VectorD<D> &v2) {
  return (v1 - v2).get_squared_magnitude();
}

//! compute the distance between two vectors
/** See VectorD
 */
template <int D>
inline double get_distance(const VectorD<D> &v1, const VectorD<D> &v2) {
  return std::sqrt(get_squared_distance(v1, v2));
}

//! Return the basis vector for the given coordinate
/** Return the unit vector pointing in the direction of the requested
    coordinate. That is
    \code
    get_basis_vector_d<3>(2)== Vector3D(0,0,1);
    \endcode
    See VectorD
 */
template <int D>
inline VectorD<D> get_basis_vector_d(unsigned int coordinate) {
  IMP_USAGE_CHECK(coordinate < D, "There are only " << D << " basis vectors");
  double vs[D];
  for (unsigned int i = 0; i < D; ++i) {
    if (i == coordinate)
      vs[i] = 1;
    else
      vs[i] = 0;
  }
  return VectorD<D>(vs, vs + D);
}

//! Return a dynamically sized basis vector
inline VectorD<-1> get_basis_vector_kd(int D, unsigned int coordinate) {
  IMP_USAGE_CHECK(D > 0, "D must be positive");
  IMP_USAGE_CHECK(coordinate < static_cast<unsigned int>(D),
                  "There are only " << D << " basis vectors");
  boost::scoped_array<double> vs(new double[D]);
  for (int i = 0; i < D; ++i) {
    if (i == static_cast<int>(coordinate))
      vs[i] = 1;
    else
      vs[i] = 0;
  }
  return VectorD<-1>(vs.get(), vs.get() + D);
}

//! Return a vector of zeros
template <int D>
inline VectorD<D> get_zero_vector_d() {
  IMP_USAGE_CHECK(D > 0, "D must be positive");
  VectorD<D> ret;
  for (int i = 0; i < D; ++i) {
    ret[i] = 0;
  }
  return ret;
}

//! Return a dynamically sized vector of zeros
template <int D>
inline VectorD<D> get_zero_vector_kd(int Di) {
  IMP_USAGE_CHECK(D == Di, "D must be positive");
  IMP_UNUSED(Di);
  return get_zero_vector_d<D>();
}

//! Return a dynamically sized vector of zeros
inline VectorD<-1> get_zero_vector_kd(int D) {
  IMP_USAGE_CHECK(D > 0, "D must be positive");
  Floats vs(D, 0);
  return VectorD<-1>(vs.begin(), vs.end());
}

//! Return a vector of ones (or another constant)
template <int D>
inline VectorD<D> get_ones_vector_d(double v = 1) {
  IMP_USAGE_CHECK(D > 0, "D must be positive");
  VectorD<D> ret;
  for (unsigned int i = 0; i < D; ++i) {
    ret[i] = v;
  }
  return ret;
}

//! Return a vector of ones (or another constant)
/** Di must equal D. */
template <int D>
inline VectorD<D> get_ones_vector_kd(unsigned int Di, double v = 1) {
  IMP_USAGE_CHECK(D == Di, "D must be equal");
  IMP_UNUSED(Di);
  return get_ones_vector_d<D>(v);
}

//! Return a vector of ones (or another constant)
inline VectorD<-1> get_ones_vector_kd(unsigned int D, double v = 1) {
  IMP_USAGE_CHECK(D > 0, "D must be positive");
  boost::scoped_array<double> vv(new double[D]);
  for (unsigned int i = 0; i < D; ++i) {
    vv[i] = v;
  }
  return VectorD<-1>(vv.get(), vv.get() + D);
}

#ifndef SWIG

/** \name Norms
    We define a number of standard, \f$L^p\f$, norms on VectorD.
    - \f$L^1\f$ is the Manhattan distance, the sum of the components
    - \f$L^2\f$ is the standard Euclidean length
    - \f$L^{\inf}\f$ is the maximum of the components
    @{
*/

template <int D>
inline double get_l2_norm(const VectorD<D> &v) {
  return v.get_magnitude();
}

template <int D>
inline double get_l1_norm(const VectorD<D> &v) {
  double n = std::abs(v[0]);
  for (unsigned int i = 1; i < v.get_dimension(); ++i) {
    n += std::abs(v[i]);
  }
  return n;
}

template <int D>
inline double get_linf_norm(const VectorD<D> &v) {
  double n = std::abs(v[0]);
  for (unsigned int i = 1; i < v.get_dimension(); ++i) {
    n = std::max(n, std::abs(v[i]));
  }
  return n;
}

/** @} */

#ifndef IMP_DOXYGEN

template <int D>
struct SpacesIO {
  const VectorD<D> &v_;
  SpacesIO(const VectorD<D> &v) : v_(v) {}
};

template <int D>
struct CommasIO {
  const VectorD<D> &v_;
  CommasIO(const VectorD<D> &v) : v_(v) {}
};
template <int D>
inline std::ostream &operator<<(std::ostream &out, const SpacesIO<D> &s) {
  s.v_.show(out, " ", false);
  return out;
}
template <int D>
inline std::ostream &operator<<(std::ostream &out, const CommasIO<D> &s) {
  s.v_.show(out, ", ", false);
  return out;
}

//! Use this before outputing to delimited vector entries with a space
/** std::cout << spaces_io(v);
    produces "1.0 2.0 3.0"
    See VectorD
 */
template <int D>
inline SpacesIO<D> spaces_io(const VectorD<D> &v) {
  return SpacesIO<D>(v);
}

//! Use this before outputing to delimited vector entries with a comma
/** std::cout << commas_io(v);
    produces "1.0, 2.0, 3.0"
    See VectorD
 */
template <int D>
inline CommasIO<D> commas_io(const VectorD<D> &v) {
  return CommasIO<D>(v);
}
#endif  // doxygen

#endif  // swig
/** 1D vector typedef for swig */
typedef VectorD<1> Vector1D;
/** 1D vectors typedef for swig */
typedef base::Vector<VectorD<1> > Vector1Ds;
/** 2D vector typedef for swig */
typedef VectorD<2> Vector2D;
/** 2D vectors typedef for swig */
typedef base::Vector<VectorD<2> > Vector2Ds;
/** 3D vector typedef for swig */
typedef VectorD<3> Vector3D;
/** 3D vectors typedef for swig */
typedef base::Vector<VectorD<3> > Vector3Ds;
/** 4D vector typedef for swig */
typedef VectorD<4> Vector4D;
/** 4D vectors typedef for swig */
typedef base::Vector<VectorD<4> > Vector4Ds;
/** 5D vector typedef for swig */
typedef VectorD<5> Vector5D;
/** 5D vectors typedef for swig */
typedef base::Vector<VectorD<5> > Vector5Ds;
/** 6D vector typedef for swig */
typedef VectorD<6> Vector6D;
/** 6D vector typedef for swig */
typedef base::Vector<VectorD<6> > Vector6Ds;
/** KD vector typedef for swig */
typedef VectorD<-1> VectorKD;
/** KD vectors typedef for swig */
typedef base::Vector<VectorD<-1> > VectorKDs;

/** See VectorD */
template <int D>
inline const VectorD<D> &get_vector_d_geometry(const VectorD<D> &g) {
  return g;
}
/** See VectorD */
template <int D>
inline void set_vector_d_geometry(VectorD<D> &g, const VectorD<D> &v) {
  g = v;
}

/** See VectorD
    Return the vector that is the elementwise product of the two.
*/
template <int D>
inline VectorD<D> get_elementwise_product(const algebra::VectorD<D> &a,
                                          const algebra::VectorD<D> &b) {
  VectorD<D> ret(a);
  for (unsigned int i = 0; i < ret.get_dimension(); ++i) {
    ret[i] *= b[i];
  }
  return ret;
}

/** See VectorD
    Return the vector that is the elementwise product of the two.
*/
template <int D>
inline VectorD<D> get_elementwise_product(const Ints &a,
                                          const algebra::VectorD<D> &b) {
  IMP_USAGE_CHECK(a.size() == b.get_dimension(), "Dimensions don't match,");
  VectorD<D> ret(b);
  for (unsigned int i = 0; i < ret.get_dimension(); ++i) {
    ret[i] *= a[i];
  }
  return ret;
}

/** A class to flexibly accept vectors as inputs to functions.
    See VectorD
 */
template <int D>
class VectorInputD : public VectorD<D>, public base::InputAdaptor {
 public:
  VectorInputD(const VectorD<D> &v) : VectorD<D>(v) {}
  VectorInputD(const Floats &v) : VectorD<D>(v, true) {}
};

/** Also accept floating point values for Vector1Ds

    See VectorD
 */
template <>
class VectorInputD<1> : public VectorD<1>, public base::InputAdaptor {
 public:
  VectorInputD(const VectorD<1> &v) : VectorD<1>(v) {}
  VectorInputD(const Floats &v) : VectorD<1>(v, true) {}
  VectorInputD(double v) : VectorD<1>(v) {}
};

/** Typedef for python. */
typedef VectorInputD<1> VectorInput1D;
/** Typedef for python. */
typedef base::Vector<VectorInputD<1> > VectorInput1Ds;
/** Typedef for python. */
typedef VectorInputD<2> VectorInput2D;
/** Typedef for python.*/
typedef base::Vector<VectorInputD<2> > VectorInput2Ds;
/** Typedef for python. */
typedef VectorInputD<3> VectorInput3D;
/** Typedef for python. */
typedef base::Vector<VectorInputD<3> > VectorInput3Ds;
/** Typedef for python. */
typedef VectorInputD<4> VectorInput4D;
/** Typedef for python. */
typedef base::Vector<VectorInputD<4> > VectorInput4Ds;
/** Typedef for python. */
typedef VectorInputD<5> VectorInput5D;
/** Typedef for python. */
typedef base::Vector<VectorInputD<5> > VectorInput5Ds;
/** Typedef for python. */
typedef VectorInputD<6> VectorInput6D;
/** Typedef for python. */
typedef base::Vector<VectorInputD<6> > VectorInput6Ds;
/** Typedef for python. */
typedef VectorInputD<-1> VectorInputKD;
/** Typedef for python. */
typedef base::Vector<VectorInputD<-1> > VectorInputKDs;

IMPALGEBRA_END_NAMESPACE

#endif /* IMPALGEBRA_VECTOR_D_H */