Simbody  3.5
VectorBase.h
Go to the documentation of this file.
00001 #ifndef SimTK_SIMMATRIX_VECTORBASE_H_
00002 #define SimTK_SIMMATRIX_VECTORBASE_H_
00003 
00004 /* -------------------------------------------------------------------------- *
00005  *                       Simbody(tm): SimTKcommon                             *
00006  * -------------------------------------------------------------------------- *
00007  * This is part of the SimTK biosimulation toolkit originating from           *
00008  * Simbios, the NIH National Center for Physics-Based Simulation of           *
00009  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
00010  * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody.  *
00011  *                                                                            *
00012  * Portions copyright (c) 2005-13 Stanford University and the Authors.        *
00013  * Authors: Michael Sherman                                                   *
00014  * Contributors:                                                              *
00015  *                                                                            *
00016  * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
00017  * not use this file except in compliance with the License. You may obtain a  *
00018  * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
00019  *                                                                            *
00020  * Unless required by applicable law or agreed to in writing, software        *
00021  * distributed under the License is distributed on an "AS IS" BASIS,          *
00022  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
00023  * See the License for the specific language governing permissions and        *
00024  * limitations under the License.                                             *
00025  * -------------------------------------------------------------------------- */
00026 
00031 namespace SimTK {
00032 
00033 //==============================================================================
00034 //                                VECTOR BASE
00035 //==============================================================================
00042 template <class ELT> class VectorBase : public MatrixBase<ELT> {
00043     typedef MatrixBase<ELT>                             Base;
00044     typedef typename Base::ScalarNormSq                 ScalarNormSq;
00045     typedef typename Base::EAbs                         EAbs;
00046     typedef typename CNT<ELT>::Scalar                   Scalar;
00047     typedef typename CNT<ELT>::Number                   Number;
00048     typedef typename CNT<ELT>::StdNumber                StdNumber;
00049     typedef VectorBase<ELT>                             T;
00050     typedef VectorBase<typename CNT<ELT>::TAbs>         TAbs;
00051     typedef VectorBase<typename CNT<ELT>::TNeg>         TNeg;
00052     typedef RowVectorView_<typename CNT<ELT>::THerm>    THerm;
00053 public:  
00054     //  ------------------------------------------------------------------------
00064 
00067     explicit VectorBase(int m=0) : Base(MatrixCommitment::Vector(), m, 1) {}
00068 
00072     VectorBase(const VectorBase& source) : Base(source) {}
00073 
00075     VectorBase(const TNeg& source) : Base(source) {}
00076 
00079     VectorBase(int m, const ELT& initialValue)
00080     :   Base(MatrixCommitment::Vector(),m,1,initialValue) {}  
00081 
00086     VectorBase(int m, const ELT* cppInitialValues)
00087     :   Base(MatrixCommitment::Vector(),m,1,cppInitialValues) {}
00089 
00090     //  ------------------------------------------------------------------------
00099 
00101     VectorBase(int m, int stride, const Scalar* s)
00102     :   Base(MatrixCommitment::Vector(m), MatrixCharacter::Vector(m),stride,s) { }
00104     VectorBase(int m, int stride, Scalar* s)
00105     :   Base(MatrixCommitment::Vector(m), MatrixCharacter::Vector(m),stride,s) { }
00107         
00108     //  ------------------------------------------------------------------------
00115 
00117     VectorBase(MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) 
00118     :   Base(MatrixCommitment::Vector(), h,s) { }
00120     VectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) 
00121     :   Base(MatrixCommitment::Vector(), h,s) { }
00123     VectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::DeepCopy& d)    
00124     :   Base(MatrixCommitment::Vector(), h,d) { }
00126 
00127     // This gives the resulting vector type when (v[i] op P) is applied to each element.
00128     // It will have element types which are the regular composite result of ELT op P.
00129     template <class P> struct EltResult { 
00130         typedef VectorBase<typename CNT<ELT>::template Result<P>::Mul> Mul;
00131         typedef VectorBase<typename CNT<ELT>::template Result<P>::Dvd> Dvd;
00132         typedef VectorBase<typename CNT<ELT>::template Result<P>::Add> Add;
00133         typedef VectorBase<typename CNT<ELT>::template Result<P>::Sub> Sub;
00134     };
00135 
00138     VectorBase& operator=(const VectorBase& b) {
00139         Base::operator=(b); return *this;
00140     }
00141 
00142     // default destructor
00143 
00144 
00145     VectorBase& operator*=(const StdNumber& t)  { Base::operator*=(t); return *this; }
00146     VectorBase& operator/=(const StdNumber& t)  { Base::operator/=(t); return *this; }
00147     VectorBase& operator+=(const VectorBase& r) { Base::operator+=(r); return *this; }
00148     VectorBase& operator-=(const VectorBase& r) { Base::operator-=(r); return *this; }  
00149 
00150 
00151     template <class EE> VectorBase& operator=(const VectorBase<EE>& b) 
00152       { Base::operator=(b);  return *this; } 
00153     template <class EE> VectorBase& operator+=(const VectorBase<EE>& b) 
00154       { Base::operator+=(b); return *this; } 
00155     template <class EE> VectorBase& operator-=(const VectorBase<EE>& b) 
00156       { Base::operator-=(b); return *this; } 
00157 
00158 
00162     VectorBase& operator=(const ELT& t) { Base::setTo(t); return *this; }  
00163 
00168     template <class EE> VectorBase& rowScaleInPlace(const VectorBase<EE>& v)
00169     { Base::template rowScaleInPlace<EE>(v); return *this; }
00170     template <class EE> inline void rowScale(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
00171     { Base::rowScale(v,out); }
00172     template <class EE> inline typename EltResult<EE>::Mul rowScale(const VectorBase<EE>& v) const
00173     { typename EltResult<EE>::Mul out(nrow()); Base::rowScale(v,out); return out; }
00174 
00179     typename CNT<ScalarNormSq>::TSqrt 
00180     normRMS(int* worstOne=0) const {
00181         if (!CNT<ELT>::IsScalar)
00182             SimTK_THROW1(Exception::Cant, 
00183                 "Vector::normRMS() only defined for scalar elements.");
00184         const int n = nrow();
00185         if (n == 0) {
00186             if (worstOne) *worstOne = -1;
00187             return typename CNT<ScalarNormSq>::TSqrt(0);
00188         }
00189 
00190         ScalarNormSq sumsq = 0;
00191         if (worstOne) {
00192             *worstOne = 0;
00193             ScalarNormSq maxsq = 0; 
00194             for (int i=0; i<n; ++i) {
00195                 const ScalarNormSq v2 = square((*this)[i]);
00196                 if (v2 > maxsq) maxsq=v2, *worstOne=i;
00197                 sumsq += v2;
00198             }
00199         } else { // don't track the worst element
00200             for (int i=0; i<n; ++i) {
00201                 const ScalarNormSq v2 = square((*this)[i]);
00202                 sumsq += v2;
00203             }
00204         }
00205 
00206         return CNT<ScalarNormSq>::sqrt(sumsq/n);
00207     }
00208 
00214     template <class EE>
00215     typename CNT<ScalarNormSq>::TSqrt 
00216     weightedNormRMS(const VectorBase<EE>& w, int* worstOne=0) const {
00217         if (!CNT<ELT>::IsScalar || !CNT<EE>::IsScalar)
00218             SimTK_THROW1(Exception::Cant, 
00219             "Vector::weightedNormRMS() only defined for scalar elements"
00220             " and weights.");
00221         const int n = nrow();
00222         assert(w.nrow()==n);
00223         if (n == 0) {
00224             if (worstOne) *worstOne = -1;
00225             return typename CNT<ScalarNormSq>::TSqrt(0);
00226         }
00227 
00228         ScalarNormSq sumsq = 0;
00229         if (worstOne) {
00230             *worstOne = 0;
00231             ScalarNormSq maxsq = 0; 
00232             for (int i=0; i<n; ++i) {
00233                 const ScalarNormSq wv2 = square(w[i]*(*this)[i]);
00234                 if (wv2 > maxsq) maxsq=wv2, *worstOne=i;
00235                 sumsq += wv2;
00236             }
00237         } else { // don't track the worst element
00238             for (int i=0; i<n; ++i) {
00239                 const ScalarNormSq wv2 = square(w[i]*(*this)[i]);
00240                 sumsq += wv2;
00241             }
00242         }
00243 
00244         return CNT<ScalarNormSq>::sqrt(sumsq/n);
00245     }
00246 
00251     EAbs normInf(int* worstOne=0) const {
00252         if (!CNT<ELT>::IsScalar)
00253             SimTK_THROW1(Exception::Cant, 
00254                 "Vector::normInf() only defined for scalar elements.");
00255         const int n = nrow();
00256         if (n == 0) {
00257             if (worstOne) *worstOne = -1;
00258             return EAbs(0);
00259         }
00260 
00261         EAbs maxabs = 0;
00262         if (worstOne) {
00263             *worstOne = 0;
00264             for (int i=0; i<n; ++i) {
00265                 const EAbs a = std::abs((*this)[i]);
00266                 if (a > maxabs) maxabs=a, *worstOne=i;
00267             }
00268         } else { // don't track the worst element
00269             for (int i=0; i<n; ++i) {
00270                 const EAbs a = std::abs((*this)[i]);
00271                 if (a > maxabs) maxabs=a;
00272             }
00273         }
00274 
00275         return maxabs;
00276     }
00277 
00283     template <class EE>
00284     EAbs weightedNormInf(const VectorBase<EE>& w, int* worstOne=0) const {
00285         if (!CNT<ELT>::IsScalar || !CNT<EE>::IsScalar)
00286             SimTK_THROW1(Exception::Cant, 
00287             "Vector::weightedNormInf() only defined for scalar elements"
00288             " and weights.");
00289         const int n = nrow();
00290         assert(w.nrow()==n);
00291         if (n == 0) {
00292             if (worstOne) *worstOne = -1;
00293             return EAbs(0);
00294         }
00295 
00296         EAbs maxabs = 0;
00297         if (worstOne) {
00298             *worstOne = 0;
00299             for (int i=0; i<n; ++i) {
00300                 const EAbs wv = std::abs(w[i]*(*this)[i]);
00301                 if (wv > maxabs) maxabs=wv, *worstOne=i;
00302             }
00303         } else { // don't track the worst element
00304             for (int i=0; i<n; ++i) {
00305                 const EAbs wv = std::abs(w[i]*(*this)[i]);
00306                 if (wv > maxabs) maxabs=wv;
00307             }
00308         }
00309 
00310         return maxabs;
00311     }
00312 
00314     VectorBase& elementwiseInvertInPlace() {
00315         Base::elementwiseInvertInPlace();
00316         return *this;
00317     }
00318 
00320     void elementwiseInvert(VectorBase<typename CNT<ELT>::TInvert>& out) const {
00321         Base::elementwiseInvert(out);
00322     }
00323 
00325     VectorBase<typename CNT<ELT>::TInvert> elementwiseInvert() const {
00326         VectorBase<typename CNT<ELT>::TInvert> out(nrow());
00327         Base::elementwiseInvert(out);
00328         return out;
00329     }
00330 
00331     // elementwise multiply
00332     template <class EE> VectorBase& elementwiseMultiplyInPlace(const VectorBase<EE>& r)
00333     { Base::template elementwiseMultiplyInPlace<EE>(r); return *this; }
00334     template <class EE> inline void elementwiseMultiply(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
00335     { Base::template elementwiseMultiply<EE>(v,out); }
00336     template <class EE> inline typename EltResult<EE>::Mul elementwiseMultiply(const VectorBase<EE>& v) const
00337     { typename EltResult<EE>::Mul out(nrow()); Base::template elementwiseMultiply<EE>(v,out); return out; }
00338 
00339     // elementwise multiply from left
00340     template <class EE> VectorBase& elementwiseMultiplyFromLeftInPlace(const VectorBase<EE>& r)
00341     { Base::template elementwiseMultiplyFromLeftInPlace<EE>(r); return *this; }
00342     template <class EE> inline void 
00343     elementwiseMultiplyFromLeft(
00344         const VectorBase<EE>& v, 
00345         typename VectorBase<EE>::template EltResult<ELT>::Mul& out) const
00346     { 
00347         Base::template elementwiseMultiplyFromLeft<EE>(v,out);
00348     }
00349     template <class EE> inline typename VectorBase<EE>::template EltResult<ELT>::Mul 
00350     elementwiseMultiplyFromLeft(const VectorBase<EE>& v) const
00351     { 
00352         typename VectorBase<EE>::template EltResult<ELT>::Mul out(nrow()); 
00353         Base::template elementwiseMultiplyFromLeft<EE>(v,out); 
00354         return out;
00355     }
00356 
00357     // elementwise divide
00358     template <class EE> VectorBase& elementwiseDivideInPlace(const VectorBase<EE>& r)
00359     { Base::template elementwiseDivideInPlace<EE>(r); return *this; }
00360     template <class EE> inline void elementwiseDivide(const VectorBase<EE>& v, typename EltResult<EE>::Dvd& out) const
00361     { Base::template elementwiseDivide<EE>(v,out); }
00362     template <class EE> inline typename EltResult<EE>::Dvd elementwiseDivide(const VectorBase<EE>& v) const
00363     { typename EltResult<EE>::Dvd out(nrow()); Base::template elementwiseDivide<EE>(v,out); return out; }
00364 
00365     // elementwise divide from left
00366     template <class EE> VectorBase& elementwiseDivideFromLeftInPlace(const VectorBase<EE>& r)
00367     { Base::template elementwiseDivideFromLeftInPlace<EE>(r); return *this; }
00368     template <class EE> inline void 
00369     elementwiseDivideFromLeft(
00370         const VectorBase<EE>& v, 
00371         typename VectorBase<EE>::template EltResult<ELT>::Dvd& out) const
00372     { 
00373         Base::template elementwiseDivideFromLeft<EE>(v,out);
00374     }
00375     template <class EE> inline typename VectorBase<EE>::template EltResult<ELT>::Dvd 
00376     elementwiseDivideFromLeft(const VectorBase<EE>& v) const
00377     { 
00378         typename VectorBase<EE>::template EltResult<ELT>::Dvd out(nrow()); 
00379         Base::template elementwiseDivideFromLeft<EE>(v,out); 
00380         return out;
00381     }
00382 
00383     // Implicit conversions are allowed to Vector or Matrix, but not to RowVector.   
00384     operator const Vector_<ELT>&()     const { return *reinterpret_cast<const Vector_<ELT>*>(this); }
00385     operator       Vector_<ELT>&()           { return *reinterpret_cast<      Vector_<ELT>*>(this); }
00386     operator const VectorView_<ELT>&() const { return *reinterpret_cast<const VectorView_<ELT>*>(this); }
00387     operator       VectorView_<ELT>&()       { return *reinterpret_cast<      VectorView_<ELT>*>(this); }
00388     
00389     operator const Matrix_<ELT>&()     const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
00390     operator       Matrix_<ELT>&()           { return *reinterpret_cast<      Matrix_<ELT>*>(this); } 
00391     operator const MatrixView_<ELT>&() const { return *reinterpret_cast<const MatrixView_<ELT>*>(this); }
00392     operator       MatrixView_<ELT>&()       { return *reinterpret_cast<      MatrixView_<ELT>*>(this); } 
00393 
00394 
00395     // size() for Vectors is Base::nelt() but returns int instead of ptrdiff_t.
00396     int size() const { 
00397     assert(Base::nelt() <= (ptrdiff_t)std::numeric_limits<int>::max()); 
00398     assert(Base::ncol()==1);
00399     return (int)Base::nelt();
00400     }
00401     int       nrow() const {assert(Base::ncol()==1); return Base::nrow();}
00402     int       ncol() const {assert(Base::ncol()==1); return Base::ncol();}
00403     ptrdiff_t nelt() const {assert(Base::ncol()==1); return Base::nelt();}
00404 
00405     // Override MatrixBase operators to return the right shape
00406     TAbs abs() const {TAbs result; Base::abs(result); return result;}
00407     
00408     // Override MatrixBase indexing operators          
00409     const ELT& operator[](int i) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(i));}
00410     ELT&       operator[](int i)       {return *reinterpret_cast<ELT*>      (Base::updHelper().updElt(i));}
00411     const ELT& operator()(int i) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(i));}
00412     ELT&       operator()(int i)       {return *reinterpret_cast<ELT*>      (Base::updHelper().updElt(i));}
00413          
00414     // Block (contiguous subvector) view creation      
00415     VectorView_<ELT> operator()(int i, int m) const {return Base::operator()(i,0,m,1).getAsVectorView();}
00416     VectorView_<ELT> operator()(int i, int m)       {return Base::operator()(i,0,m,1).updAsVectorView();}
00417 
00418     // Indexed view creation (arbitrary subvector). Indices must be 
00419     // monotonically increasing.
00420     VectorView_<ELT> index(const Array_<int>& indices) const {
00421         MatrixHelper<Scalar> h(Base::getHelper().getCharacterCommitment(), 
00422                                Base::getHelper(), indices);
00423         return VectorView_<ELT>(h);
00424     }
00425     VectorView_<ELT> updIndex(const Array_<int>& indices) {
00426         MatrixHelper<Scalar> h(Base::getHelper().getCharacterCommitment(), 
00427                                Base::updHelper(), indices);
00428         return VectorView_<ELT>(h);
00429     }
00430 
00431     VectorView_<ELT> operator()(const Array_<int>& indices) const {return index(indices);}
00432     VectorView_<ELT> operator()(const Array_<int>& indices)       {return updIndex(indices);}
00433  
00434     // Hermitian transpose.
00435     THerm transpose() const {return Base::transpose().getAsRowVectorView();}
00436     THerm updTranspose()    {return Base::updTranspose().updAsRowVectorView();}
00437 
00438     THerm operator~() const {return transpose();}
00439     THerm operator~()       {return updTranspose();}
00440 
00441     const VectorBase& operator+() const {return *this; }
00442 
00443     // Negation
00444 
00445     const TNeg& negate()    const {return *reinterpret_cast<const TNeg*>(this); }
00446     TNeg&       updNegate()       {return *reinterpret_cast<TNeg*>(this); }
00447 
00448     const TNeg& operator-() const {return negate();}
00449     TNeg&       operator-()       {return updNegate();}
00450 
00451     VectorBase& resize(int m)     {Base::resize(m,1); return *this;}
00452     VectorBase& resizeKeep(int m) {Base::resizeKeep(m,1); return *this;}
00453 
00454     //TODO: this is not re-locking the number of columns at 1.
00455     void clear() {Base::clear(); Base::resize(0,1);}
00456 
00457     ELT sum() const {ELT s; Base::getHelper().sum(reinterpret_cast<Scalar*>(&s)); return s; } // add all the elements        
00458     VectorIterator<ELT, VectorBase<ELT> > begin() {
00459         return VectorIterator<ELT, VectorBase<ELT> >(*this, 0);
00460     }
00461     VectorIterator<ELT, VectorBase<ELT> > end() {
00462         return VectorIterator<ELT, VectorBase<ELT> >(*this, size());
00463     }
00464 
00465 protected:
00466     // Create a VectorBase handle using a given helper rep. 
00467     explicit VectorBase(MatrixHelperRep<Scalar>* hrep) : Base(hrep) {}
00468 
00469 private:
00470     // NO DATA MEMBERS ALLOWED
00471 };
00472 
00473 } //namespace SimTK
00474 
00475 #endif // SimTK_SIMMATRIX_VECTORBASE_H_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines