Simbody
3.5
|
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_