VectorBase.h

Go to the documentation of this file.

#ifndef SimTK_SIMMATRIX_VECTORBASE_H_
#define SimTK_SIMMATRIX_VECTORBASE_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) 2005-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 <initializer_list>
 
namespace SimTK {
 
//==============================================================================
//                                VECTOR BASE
//==============================================================================
template <class ELT> class VectorBase : public MatrixBase<ELT> {
    typedef MatrixBase<ELT>                             Base;
    typedef typename Base::ScalarNormSq                 ScalarNormSq;
    typedef typename Base::EAbs                         EAbs;
    typedef typename CNT<ELT>::Scalar                   Scalar;
    typedef typename CNT<ELT>::Number                   Number;
    typedef typename CNT<ELT>::StdNumber                StdNumber;
    typedef VectorBase<ELT>                             T;
    typedef VectorBase<typename CNT<ELT>::TAbs>         TAbs;
    typedef VectorBase<typename CNT<ELT>::TNeg>         TNeg;
    typedef RowVectorView_<typename CNT<ELT>::THerm>    THerm;
public:  
    //  ------------------------------------------------------------------------
 
    explicit VectorBase(int m=0) : Base(MatrixCommitment::Vector(), m, 1) {}
 
    VectorBase(const VectorBase& source) : Base(source) {}
 
    VectorBase(const TNeg& source) : Base(source) {}
 
    VectorBase(int m, const ELT& initialValue)
    :   Base(MatrixCommitment::Vector(),m,1,initialValue) {}  
 
    VectorBase(int m, const ELT* cppInitialValues)
    :   Base(MatrixCommitment::Vector(),m,1,cppInitialValues) {}
 
    VectorBase(std::initializer_list<ELT> init)
    :   VectorBase(static_cast<int>(init.size()), init.begin()) {}
 
    //  ------------------------------------------------------------------------
 
    VectorBase(int m, int stride, const Scalar* s)
    :   Base(MatrixCommitment::Vector(m), MatrixCharacter::Vector(m),stride,s) { }
    VectorBase(int m, int stride, Scalar* s)
    :   Base(MatrixCommitment::Vector(m), MatrixCharacter::Vector(m),stride,s) { }
        
    //  ------------------------------------------------------------------------
 
    VectorBase(MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) 
    :   Base(MatrixCommitment::Vector(), h,s) { }
    VectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) 
    :   Base(MatrixCommitment::Vector(), h,s) { }
    VectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::DeepCopy& d)    
    :   Base(MatrixCommitment::Vector(), h,d) { }
 
    // This gives the resulting vector type when (v[i] op P) is applied to each element.
    // It will have element types which are the regular composite result of ELT op P.
    template <class P> struct EltResult { 
        typedef VectorBase<typename CNT<ELT>::template Result<P>::Mul> Mul;
        typedef VectorBase<typename CNT<ELT>::template Result<P>::Dvd> Dvd;
        typedef VectorBase<typename CNT<ELT>::template Result<P>::Add> Add;
        typedef VectorBase<typename CNT<ELT>::template Result<P>::Sub> Sub;
    };
 
    VectorBase& operator=(const VectorBase& b) {
        Base::operator=(b); return *this;
    }
 
    // default destructor
 
 
    VectorBase& operator*=(const StdNumber& t)  { Base::operator*=(t); return *this; }
    VectorBase& operator/=(const StdNumber& t)  { Base::operator/=(t); return *this; }
    VectorBase& operator+=(const VectorBase& r) { Base::operator+=(r); return *this; }
    VectorBase& operator-=(const VectorBase& r) { Base::operator-=(r); return *this; }  
 
 
    template <class EE> VectorBase& operator=(const VectorBase<EE>& b) 
      { Base::operator=(b);  return *this; } 
    template <class EE> VectorBase& operator+=(const VectorBase<EE>& b) 
      { Base::operator+=(b); return *this; } 
    template <class EE> VectorBase& operator-=(const VectorBase<EE>& b) 
      { Base::operator-=(b); return *this; } 
 
 
    VectorBase& operator=(const ELT& t) { Base::setTo(t); return *this; }  
 
    template <class EE> VectorBase& rowScaleInPlace(const VectorBase<EE>& v)
    { Base::template rowScaleInPlace<EE>(v); return *this; }
    template <class EE> inline void rowScale(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
    { Base::rowScale(v,out); }
    template <class EE> SimTK_NODISCARD inline typename EltResult<EE>::Mul rowScale(const VectorBase<EE>& v) const
    { typename EltResult<EE>::Mul out(nrow()); Base::rowScale(v,out); return out; }
 
    typename CNT<ScalarNormSq>::TSqrt 
    normRMS(int* worstOne=0) const {
        if (!CNT<ELT>::IsScalar)
            SimTK_THROW1(Exception::Cant, 
                "Vector::normRMS() only defined for scalar elements.");
        const int n = nrow();
        if (n == 0) {
            if (worstOne) *worstOne = -1;
            return typename CNT<ScalarNormSq>::TSqrt(0);
        }
 
        ScalarNormSq sumsq = 0;
        if (worstOne) {
            *worstOne = 0;
            ScalarNormSq maxsq = 0; 
            for (int i=0; i<n; ++i) {
                const ScalarNormSq v2 = square((*this)[i]);
                if (v2 > maxsq) maxsq=v2, *worstOne=i;
                sumsq += v2;
            }
        } else { // don't track the worst element
            for (int i=0; i<n; ++i) {
                const ScalarNormSq v2 = square((*this)[i]);
                sumsq += v2;
            }
        }
 
        return CNT<ScalarNormSq>::sqrt(sumsq/n);
    }
 
    template <class EE>
    typename CNT<ScalarNormSq>::TSqrt 
    weightedNormRMS(const VectorBase<EE>& w, int* worstOne=0) const {
        if (!CNT<ELT>::IsScalar || !CNT<EE>::IsScalar)
            SimTK_THROW1(Exception::Cant, 
            "Vector::weightedNormRMS() only defined for scalar elements"
            " and weights.");
        const int n = nrow();
        assert(w.nrow()==n);
        if (n == 0) {
            if (worstOne) *worstOne = -1;
            return typename CNT<ScalarNormSq>::TSqrt(0);
        }
 
        ScalarNormSq sumsq = 0;
        if (worstOne) {
            *worstOne = 0;
            ScalarNormSq maxsq = 0; 
            for (int i=0; i<n; ++i) {
                const ScalarNormSq wv2 = square(w[i]*(*this)[i]);
                if (wv2 > maxsq) maxsq=wv2, *worstOne=i;
                sumsq += wv2;
            }
        } else { // don't track the worst element
            for (int i=0; i<n; ++i) {
                const ScalarNormSq wv2 = square(w[i]*(*this)[i]);
                sumsq += wv2;
            }
        }
 
        return CNT<ScalarNormSq>::sqrt(sumsq/n);
    }
 
    EAbs normInf(int* worstOne=0) const {
        if (!CNT<ELT>::IsScalar)
            SimTK_THROW1(Exception::Cant, 
                "Vector::normInf() only defined for scalar elements.");
        const int n = nrow();
        if (n == 0) {
            if (worstOne) *worstOne = -1;
            return EAbs(0);
        }
 
        EAbs maxabs = 0;
        if (worstOne) {
            *worstOne = 0;
            for (int i=0; i<n; ++i) {
                const EAbs a = std::abs((*this)[i]);
                if (a > maxabs) maxabs=a, *worstOne=i;
            }
        } else { // don't track the worst element
            for (int i=0; i<n; ++i) {
                const EAbs a = std::abs((*this)[i]);
                if (a > maxabs) maxabs=a;
            }
        }
 
        return maxabs;
    }
 
    template <class EE>
    EAbs weightedNormInf(const VectorBase<EE>& w, int* worstOne=0) const {
        if (!CNT<ELT>::IsScalar || !CNT<EE>::IsScalar)
            SimTK_THROW1(Exception::Cant, 
            "Vector::weightedNormInf() only defined for scalar elements"
            " and weights.");
        const int n = nrow();
        assert(w.nrow()==n);
        if (n == 0) {
            if (worstOne) *worstOne = -1;
            return EAbs(0);
        }
 
        EAbs maxabs = 0;
        if (worstOne) {
            *worstOne = 0;
            for (int i=0; i<n; ++i) {
                const EAbs wv = std::abs(w[i]*(*this)[i]);
                if (wv > maxabs) maxabs=wv, *worstOne=i;
            }
        } else { // don't track the worst element
            for (int i=0; i<n; ++i) {
                const EAbs wv = std::abs(w[i]*(*this)[i]);
                if (wv > maxabs) maxabs=wv;
            }
        }
 
        return maxabs;
    }
 
    VectorBase& elementwiseInvertInPlace() {
        Base::elementwiseInvertInPlace();
        return *this;
    }
 
    void elementwiseInvert(VectorBase<typename CNT<ELT>::TInvert>& out) const {
        Base::elementwiseInvert(out);
    }
 
    SimTK_NODISCARD VectorBase<typename CNT<ELT>::TInvert> elementwiseInvert() const {
        VectorBase<typename CNT<ELT>::TInvert> out(nrow());
        Base::elementwiseInvert(out);
        return out;
    }
 
    // elementwise multiply
    template <class EE> VectorBase& elementwiseMultiplyInPlace(const VectorBase<EE>& r)
    { Base::template elementwiseMultiplyInPlace<EE>(r); return *this; }
    template <class EE> inline void elementwiseMultiply(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
    { Base::template elementwiseMultiply<EE>(v,out); }
    template <class EE> inline typename EltResult<EE>::Mul elementwiseMultiply(const VectorBase<EE>& v) const
    { typename EltResult<EE>::Mul out(nrow()); Base::template elementwiseMultiply<EE>(v,out); return out; }
 
    // elementwise multiply from left
    template <class EE> VectorBase& elementwiseMultiplyFromLeftInPlace(const VectorBase<EE>& r)
    { Base::template elementwiseMultiplyFromLeftInPlace<EE>(r); return *this; }
    template <class EE> inline void 
    elementwiseMultiplyFromLeft(
        const VectorBase<EE>& v, 
        typename VectorBase<EE>::template EltResult<ELT>::Mul& out) const
    { 
        Base::template elementwiseMultiplyFromLeft<EE>(v,out);
    }
    template <class EE> inline typename VectorBase<EE>::template EltResult<ELT>::Mul 
    elementwiseMultiplyFromLeft(const VectorBase<EE>& v) const
    { 
        typename VectorBase<EE>::template EltResult<ELT>::Mul out(nrow()); 
        Base::template elementwiseMultiplyFromLeft<EE>(v,out); 
        return out;
    }
 
    // elementwise divide
    template <class EE> VectorBase& elementwiseDivideInPlace(const VectorBase<EE>& r)
    { Base::template elementwiseDivideInPlace<EE>(r); return *this; }
    template <class EE> inline void elementwiseDivide(const VectorBase<EE>& v, typename EltResult<EE>::Dvd& out) const
    { Base::template elementwiseDivide<EE>(v,out); }
    template <class EE> inline typename EltResult<EE>::Dvd elementwiseDivide(const VectorBase<EE>& v) const
    { typename EltResult<EE>::Dvd out(nrow()); Base::template elementwiseDivide<EE>(v,out); return out; }
 
    // elementwise divide from left
    template <class EE> VectorBase& elementwiseDivideFromLeftInPlace(const VectorBase<EE>& r)
    { Base::template elementwiseDivideFromLeftInPlace<EE>(r); return *this; }
    template <class EE> inline void 
    elementwiseDivideFromLeft(
        const VectorBase<EE>& v, 
        typename VectorBase<EE>::template EltResult<ELT>::Dvd& out) const
    { 
        Base::template elementwiseDivideFromLeft<EE>(v,out);
    }
    template <class EE> inline typename VectorBase<EE>::template EltResult<ELT>::Dvd 
    elementwiseDivideFromLeft(const VectorBase<EE>& v) const
    { 
        typename VectorBase<EE>::template EltResult<ELT>::Dvd out(nrow()); 
        Base::template elementwiseDivideFromLeft<EE>(v,out); 
        return out;
    }
 
    // Implicit conversions are allowed to Vector or Matrix, but not to RowVector.   
    operator const Vector_<ELT>&()     const { return *reinterpret_cast<const Vector_<ELT>*>(this); }
    operator       Vector_<ELT>&()           { return *reinterpret_cast<      Vector_<ELT>*>(this); }
    operator const VectorView_<ELT>&() const { return *reinterpret_cast<const VectorView_<ELT>*>(this); }
    operator       VectorView_<ELT>&()       { return *reinterpret_cast<      VectorView_<ELT>*>(this); }
    
    operator const Matrix_<ELT>&()     const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
    operator       Matrix_<ELT>&()           { return *reinterpret_cast<      Matrix_<ELT>*>(this); } 
    operator const MatrixView_<ELT>&() const { return *reinterpret_cast<const MatrixView_<ELT>*>(this); }
    operator       MatrixView_<ELT>&()       { return *reinterpret_cast<      MatrixView_<ELT>*>(this); } 
 
 
    // size() for Vectors is Base::nelt() but returns int instead of ptrdiff_t.
    int size() const { 
    assert(Base::nelt() <= (ptrdiff_t)std::numeric_limits<int>::max()); 
    assert(Base::ncol()==1);
    return (int)Base::nelt();
    }
    int       nrow() const {assert(Base::ncol()==1); return Base::nrow();}
    int       ncol() const {assert(Base::ncol()==1); return Base::ncol();}
    ptrdiff_t nelt() const {assert(Base::ncol()==1); return Base::nelt();}
 
    // Override MatrixBase operators to return the right shape
    SimTK_NODISCARD TAbs abs() const {TAbs result; Base::abs(result); return result;}
    
    // Override MatrixBase indexing operators          
    const ELT& operator[](int i) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(i));}
    ELT&       operator[](int i)       {return *reinterpret_cast<ELT*>      (Base::updHelper().updElt(i));}
    const ELT& operator()(int i) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(i));}
    ELT&       operator()(int i)       {return *reinterpret_cast<ELT*>      (Base::updHelper().updElt(i));}
         
    // Block (contiguous subvector) view creation      
    VectorView_<ELT> operator()(int i, int m) const {return Base::operator()(i,0,m,1).getAsVectorView();}
    VectorView_<ELT> operator()(int i, int m)       {return Base::operator()(i,0,m,1).updAsVectorView();}
 
    // Indexed view creation (arbitrary subvector). Indices must be 
    // monotonically increasing.
    VectorView_<ELT> index(const Array_<int>& indices) const {
        MatrixHelper<Scalar> h(Base::getHelper().getCharacterCommitment(), 
                               Base::getHelper(), indices);
        return VectorView_<ELT>(h);
    }
    VectorView_<ELT> updIndex(const Array_<int>& indices) {
        MatrixHelper<Scalar> h(Base::getHelper().getCharacterCommitment(), 
                               Base::updHelper(), indices);
        return VectorView_<ELT>(h);
    }
 
    VectorView_<ELT> operator()(const Array_<int>& indices) const {return index(indices);}
    VectorView_<ELT> operator()(const Array_<int>& indices)       {return updIndex(indices);}
 
    // Hermitian transpose.
    SimTK_NODISCARD THerm transpose() const {return Base::transpose().getAsRowVectorView();}
    THerm updTranspose()    {return Base::updTranspose().updAsRowVectorView();}
 
    THerm operator~() const {return transpose();}
    THerm operator~()       {return updTranspose();}
 
    const VectorBase& operator+() const {return *this; }
 
    // Negation
 
    SimTK_NODISCARD const TNeg& negate()    const {return *reinterpret_cast<const TNeg*>(this); }
    TNeg&       updNegate()       {return *reinterpret_cast<TNeg*>(this); }
 
    const TNeg& operator-() const {return negate();}
    TNeg&       operator-()       {return updNegate();}
 
    VectorBase& resize(int m)     {Base::resize(m,1); return *this;}
    VectorBase& resizeKeep(int m) {Base::resizeKeep(m,1); return *this;}
 
    //TODO: this is not re-locking the number of columns at 1.
    void clear() {Base::clear(); Base::resize(0,1);}
 
    ELT sum() const {ELT s; Base::getHelper().sum(reinterpret_cast<Scalar*>(&s)); return s; } // add all the elements        
    VectorIterator<ELT, VectorBase<ELT> > begin() {
        return VectorIterator<ELT, VectorBase<ELT> >(*this, 0);
    }
    VectorIterator<ELT, VectorBase<ELT> > end() {
        return VectorIterator<ELT, VectorBase<ELT> >(*this, size());
    }
 
protected:
    // Create a VectorBase handle using a given helper rep. 
    explicit VectorBase(MatrixHelperRep<Scalar>* hrep) : Base(hrep) {}
 
private:
    // NO DATA MEMBERS ALLOWED
};
 
} //namespace SimTK
 
#endif // SimTK_SIMMATRIX_VECTORBASE_H_