Simbody  3.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VectorBase.h
Go to the documentation of this file.
1 #ifndef SimTK_SIMMATRIX_VECTORBASE_H_
2 #define SimTK_SIMMATRIX_VECTORBASE_H_
3 
4 /* -------------------------------------------------------------------------- *
5  * Simbody(tm): SimTKcommon *
6  * -------------------------------------------------------------------------- *
7  * This is part of the SimTK biosimulation toolkit originating from *
8  * Simbios, the NIH National Center for Physics-Based Simulation of *
9  * Biological Structures at Stanford, funded under the NIH Roadmap for *
10  * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. *
11  * *
12  * Portions copyright (c) 2005-13 Stanford University and the Authors. *
13  * Authors: Michael Sherman *
14  * Contributors: *
15  * *
16  * Licensed under the Apache License, Version 2.0 (the "License"); you may *
17  * not use this file except in compliance with the License. You may obtain a *
18  * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
19  * *
20  * Unless required by applicable law or agreed to in writing, software *
21  * distributed under the License is distributed on an "AS IS" BASIS, *
22  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
23  * See the License for the specific language governing permissions and *
24  * limitations under the License. *
25  * -------------------------------------------------------------------------- */
26 
31 namespace SimTK {
32 
33 //==============================================================================
34 // VECTOR BASE
35 //==============================================================================
42 template <class ELT> class VectorBase : public MatrixBase<ELT> {
43  typedef MatrixBase<ELT> Base;
44  typedef typename Base::ScalarNormSq ScalarNormSq;
45  typedef typename Base::EAbs EAbs;
46  typedef typename CNT<ELT>::Scalar Scalar;
47  typedef typename CNT<ELT>::Number Number;
48  typedef typename CNT<ELT>::StdNumber StdNumber;
49  typedef VectorBase<ELT> T;
50  typedef VectorBase<typename CNT<ELT>::TAbs> TAbs;
51  typedef VectorBase<typename CNT<ELT>::TNeg> TNeg;
52  typedef RowVectorView_<typename CNT<ELT>::THerm> THerm;
53 public:
54  // ------------------------------------------------------------------------
64 
67  explicit VectorBase(int m=0) : Base(MatrixCommitment::Vector(), m, 1) {}
68 
72  VectorBase(const VectorBase& source) : Base(source) {}
73 
75  VectorBase(const TNeg& source) : Base(source) {}
76 
79  VectorBase(int m, const ELT& initialValue)
80  : Base(MatrixCommitment::Vector(),m,1,initialValue) {}
81 
86  VectorBase(int m, const ELT* cppInitialValues)
87  : Base(MatrixCommitment::Vector(),m,1,cppInitialValues) {}
89 
90  // ------------------------------------------------------------------------
99 
101  VectorBase(int m, int stride, const Scalar* s)
102  : Base(MatrixCommitment::Vector(m), MatrixCharacter::Vector(m),stride,s) { }
104  VectorBase(int m, int stride, Scalar* s)
105  : Base(MatrixCommitment::Vector(m), MatrixCharacter::Vector(m),stride,s) { }
107 
108  // ------------------------------------------------------------------------
115 
117  VectorBase(MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s)
118  : Base(MatrixCommitment::Vector(), h,s) { }
120  VectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s)
121  : Base(MatrixCommitment::Vector(), h,s) { }
123  VectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::DeepCopy& d)
124  : Base(MatrixCommitment::Vector(), h,d) { }
126 
127  // This gives the resulting vector type when (v[i] op P) is applied to each element.
128  // It will have element types which are the regular composite result of ELT op P.
129  template <class P> struct EltResult {
130  typedef VectorBase<typename CNT<ELT>::template Result<P>::Mul> Mul;
131  typedef VectorBase<typename CNT<ELT>::template Result<P>::Dvd> Dvd;
132  typedef VectorBase<typename CNT<ELT>::template Result<P>::Add> Add;
133  typedef VectorBase<typename CNT<ELT>::template Result<P>::Sub> Sub;
134  };
135 
138  VectorBase& operator=(const VectorBase& b) {
139  Base::operator=(b); return *this;
140  }
141 
142  // default destructor
143 
144 
145  VectorBase& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
146  VectorBase& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
147  VectorBase& operator+=(const VectorBase& r) { Base::operator+=(r); return *this; }
148  VectorBase& operator-=(const VectorBase& r) { Base::operator-=(r); return *this; }
149 
150 
151  template <class EE> VectorBase& operator=(const VectorBase<EE>& b)
152  { Base::operator=(b); return *this; }
153  template <class EE> VectorBase& operator+=(const VectorBase<EE>& b)
154  { Base::operator+=(b); return *this; }
155  template <class EE> VectorBase& operator-=(const VectorBase<EE>& b)
156  { Base::operator-=(b); return *this; }
157 
158 
162  VectorBase& operator=(const ELT& t) { Base::setTo(t); return *this; }
163 
168  template <class EE> VectorBase& rowScaleInPlace(const VectorBase<EE>& v)
169  { Base::template rowScaleInPlace<EE>(v); return *this; }
170  template <class EE> inline void rowScale(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
171  { Base::rowScale(v,out); }
172  template <class EE> inline typename EltResult<EE>::Mul rowScale(const VectorBase<EE>& v) const
173  { typename EltResult<EE>::Mul out(nrow()); Base::rowScale(v,out); return out; }
174 
179  typename CNT<ScalarNormSq>::TSqrt
180  normRMS(int* worstOne=0) const {
181  if (!CNT<ELT>::IsScalar)
182  SimTK_THROW1(Exception::Cant,
183  "Vector::normRMS() only defined for scalar elements.");
184  const int n = nrow();
185  if (n == 0) {
186  if (worstOne) *worstOne = -1;
187  return typename CNT<ScalarNormSq>::TSqrt(0);
188  }
189 
190  ScalarNormSq sumsq = 0;
191  if (worstOne) {
192  *worstOne = 0;
193  ScalarNormSq maxsq = 0;
194  for (int i=0; i<n; ++i) {
195  const ScalarNormSq v2 = square((*this)[i]);
196  if (v2 > maxsq) maxsq=v2, *worstOne=i;
197  sumsq += v2;
198  }
199  } else { // don't track the worst element
200  for (int i=0; i<n; ++i) {
201  const ScalarNormSq v2 = square((*this)[i]);
202  sumsq += v2;
203  }
204  }
205 
206  return CNT<ScalarNormSq>::sqrt(sumsq/n);
207  }
208 
214  template <class EE>
215  typename CNT<ScalarNormSq>::TSqrt
216  weightedNormRMS(const VectorBase<EE>& w, int* worstOne=0) const {
217  if (!CNT<ELT>::IsScalar || !CNT<EE>::IsScalar)
218  SimTK_THROW1(Exception::Cant,
219  "Vector::weightedNormRMS() only defined for scalar elements"
220  " and weights.");
221  const int n = nrow();
222  assert(w.nrow()==n);
223  if (n == 0) {
224  if (worstOne) *worstOne = -1;
225  return typename CNT<ScalarNormSq>::TSqrt(0);
226  }
227 
228  ScalarNormSq sumsq = 0;
229  if (worstOne) {
230  *worstOne = 0;
231  ScalarNormSq maxsq = 0;
232  for (int i=0; i<n; ++i) {
233  const ScalarNormSq wv2 = square(w[i]*(*this)[i]);
234  if (wv2 > maxsq) maxsq=wv2, *worstOne=i;
235  sumsq += wv2;
236  }
237  } else { // don't track the worst element
238  for (int i=0; i<n; ++i) {
239  const ScalarNormSq wv2 = square(w[i]*(*this)[i]);
240  sumsq += wv2;
241  }
242  }
243 
244  return CNT<ScalarNormSq>::sqrt(sumsq/n);
245  }
246 
251  EAbs normInf(int* worstOne=0) const {
252  if (!CNT<ELT>::IsScalar)
253  SimTK_THROW1(Exception::Cant,
254  "Vector::normInf() only defined for scalar elements.");
255  const int n = nrow();
256  if (n == 0) {
257  if (worstOne) *worstOne = -1;
258  return EAbs(0);
259  }
260 
261  EAbs maxabs = 0;
262  if (worstOne) {
263  *worstOne = 0;
264  for (int i=0; i<n; ++i) {
265  const EAbs a = std::abs((*this)[i]);
266  if (a > maxabs) maxabs=a, *worstOne=i;
267  }
268  } else { // don't track the worst element
269  for (int i=0; i<n; ++i) {
270  const EAbs a = std::abs((*this)[i]);
271  if (a > maxabs) maxabs=a;
272  }
273  }
274 
275  return maxabs;
276  }
277 
283  template <class EE>
284  EAbs weightedNormInf(const VectorBase<EE>& w, int* worstOne=0) const {
285  if (!CNT<ELT>::IsScalar || !CNT<EE>::IsScalar)
286  SimTK_THROW1(Exception::Cant,
287  "Vector::weightedNormInf() only defined for scalar elements"
288  " and weights.");
289  const int n = nrow();
290  assert(w.nrow()==n);
291  if (n == 0) {
292  if (worstOne) *worstOne = -1;
293  return EAbs(0);
294  }
295 
296  EAbs maxabs = 0;
297  if (worstOne) {
298  *worstOne = 0;
299  for (int i=0; i<n; ++i) {
300  const EAbs wv = std::abs(w[i]*(*this)[i]);
301  if (wv > maxabs) maxabs=wv, *worstOne=i;
302  }
303  } else { // don't track the worst element
304  for (int i=0; i<n; ++i) {
305  const EAbs wv = std::abs(w[i]*(*this)[i]);
306  if (wv > maxabs) maxabs=wv;
307  }
308  }
309 
310  return maxabs;
311  }
312 
314  VectorBase& elementwiseInvertInPlace() {
315  Base::elementwiseInvertInPlace();
316  return *this;
317  }
318 
320  void elementwiseInvert(VectorBase<typename CNT<ELT>::TInvert>& out) const {
321  Base::elementwiseInvert(out);
322  }
323 
325  VectorBase<typename CNT<ELT>::TInvert> elementwiseInvert() const {
326  VectorBase<typename CNT<ELT>::TInvert> out(nrow());
327  Base::elementwiseInvert(out);
328  return out;
329  }
330 
331  // elementwise multiply
332  template <class EE> VectorBase& elementwiseMultiplyInPlace(const VectorBase<EE>& r)
333  { Base::template elementwiseMultiplyInPlace<EE>(r); return *this; }
334  template <class EE> inline void elementwiseMultiply(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
335  { Base::template elementwiseMultiply<EE>(v,out); }
336  template <class EE> inline typename EltResult<EE>::Mul elementwiseMultiply(const VectorBase<EE>& v) const
337  { typename EltResult<EE>::Mul out(nrow()); Base::template elementwiseMultiply<EE>(v,out); return out; }
338 
339  // elementwise multiply from left
340  template <class EE> VectorBase& elementwiseMultiplyFromLeftInPlace(const VectorBase<EE>& r)
341  { Base::template elementwiseMultiplyFromLeftInPlace<EE>(r); return *this; }
342  template <class EE> inline void
343  elementwiseMultiplyFromLeft(
344  const VectorBase<EE>& v,
345  typename VectorBase<EE>::template EltResult<ELT>::Mul& out) const
346  {
347  Base::template elementwiseMultiplyFromLeft<EE>(v,out);
348  }
349  template <class EE> inline typename VectorBase<EE>::template EltResult<ELT>::Mul
350  elementwiseMultiplyFromLeft(const VectorBase<EE>& v) const
351  {
352  typename VectorBase<EE>::template EltResult<ELT>::Mul out(nrow());
353  Base::template elementwiseMultiplyFromLeft<EE>(v,out);
354  return out;
355  }
356 
357  // elementwise divide
358  template <class EE> VectorBase& elementwiseDivideInPlace(const VectorBase<EE>& r)
359  { Base::template elementwiseDivideInPlace<EE>(r); return *this; }
360  template <class EE> inline void elementwiseDivide(const VectorBase<EE>& v, typename EltResult<EE>::Dvd& out) const
361  { Base::template elementwiseDivide<EE>(v,out); }
362  template <class EE> inline typename EltResult<EE>::Dvd elementwiseDivide(const VectorBase<EE>& v) const
363  { typename EltResult<EE>::Dvd out(nrow()); Base::template elementwiseDivide<EE>(v,out); return out; }
364 
365  // elementwise divide from left
366  template <class EE> VectorBase& elementwiseDivideFromLeftInPlace(const VectorBase<EE>& r)
367  { Base::template elementwiseDivideFromLeftInPlace<EE>(r); return *this; }
368  template <class EE> inline void
369  elementwiseDivideFromLeft(
370  const VectorBase<EE>& v,
371  typename VectorBase<EE>::template EltResult<ELT>::Dvd& out) const
372  {
373  Base::template elementwiseDivideFromLeft<EE>(v,out);
374  }
375  template <class EE> inline typename VectorBase<EE>::template EltResult<ELT>::Dvd
376  elementwiseDivideFromLeft(const VectorBase<EE>& v) const
377  {
378  typename VectorBase<EE>::template EltResult<ELT>::Dvd out(nrow());
379  Base::template elementwiseDivideFromLeft<EE>(v,out);
380  return out;
381  }
382 
383  // Implicit conversions are allowed to Vector or Matrix, but not to RowVector.
384  operator const Vector_<ELT>&() const { return *reinterpret_cast<const Vector_<ELT>*>(this); }
385  operator Vector_<ELT>&() { return *reinterpret_cast< Vector_<ELT>*>(this); }
386  operator const VectorView_<ELT>&() const { return *reinterpret_cast<const VectorView_<ELT>*>(this); }
387  operator VectorView_<ELT>&() { return *reinterpret_cast< VectorView_<ELT>*>(this); }
388 
389  operator const Matrix_<ELT>&() const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
390  operator Matrix_<ELT>&() { return *reinterpret_cast< Matrix_<ELT>*>(this); }
391  operator const MatrixView_<ELT>&() const { return *reinterpret_cast<const MatrixView_<ELT>*>(this); }
392  operator MatrixView_<ELT>&() { return *reinterpret_cast< MatrixView_<ELT>*>(this); }
393 
394 
395  // size() for Vectors is Base::nelt() but returns int instead of ptrdiff_t.
396  int size() const {
397  assert(Base::nelt() <= (ptrdiff_t)std::numeric_limits<int>::max());
398  assert(Base::ncol()==1);
399  return (int)Base::nelt();
400  }
401  int nrow() const {assert(Base::ncol()==1); return Base::nrow();}
402  int ncol() const {assert(Base::ncol()==1); return Base::ncol();}
403  ptrdiff_t nelt() const {assert(Base::ncol()==1); return Base::nelt();}
404 
405  // Override MatrixBase operators to return the right shape
406  TAbs abs() const {TAbs result; Base::abs(result); return result;}
407 
408  // Override MatrixBase indexing operators
409  const ELT& operator[](int i) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(i));}
410  ELT& operator[](int i) {return *reinterpret_cast<ELT*> (Base::updHelper().updElt(i));}
411  const ELT& operator()(int i) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(i));}
412  ELT& operator()(int i) {return *reinterpret_cast<ELT*> (Base::updHelper().updElt(i));}
413 
414  // Block (contiguous subvector) view creation
415  VectorView_<ELT> operator()(int i, int m) const {return Base::operator()(i,0,m,1).getAsVectorView();}
416  VectorView_<ELT> operator()(int i, int m) {return Base::operator()(i,0,m,1).updAsVectorView();}
417 
418  // Indexed view creation (arbitrary subvector). Indices must be
419  // monotonically increasing.
420  VectorView_<ELT> index(const Array_<int>& indices) const {
421  MatrixHelper<Scalar> h(Base::getHelper().getCharacterCommitment(),
422  Base::getHelper(), indices);
423  return VectorView_<ELT>(h);
424  }
425  VectorView_<ELT> updIndex(const Array_<int>& indices) {
426  MatrixHelper<Scalar> h(Base::getHelper().getCharacterCommitment(),
427  Base::updHelper(), indices);
428  return VectorView_<ELT>(h);
429  }
430 
431  VectorView_<ELT> operator()(const Array_<int>& indices) const {return index(indices);}
432  VectorView_<ELT> operator()(const Array_<int>& indices) {return updIndex(indices);}
433 
434  // Hermitian transpose.
435  THerm transpose() const {return Base::transpose().getAsRowVectorView();}
436  THerm updTranspose() {return Base::updTranspose().updAsRowVectorView();}
437 
438  THerm operator~() const {return transpose();}
439  THerm operator~() {return updTranspose();}
440 
441  const VectorBase& operator+() const {return *this; }
442 
443  // Negation
444 
445  const TNeg& negate() const {return *reinterpret_cast<const TNeg*>(this); }
446  TNeg& updNegate() {return *reinterpret_cast<TNeg*>(this); }
447 
448  const TNeg& operator-() const {return negate();}
449  TNeg& operator-() {return updNegate();}
450 
451  VectorBase& resize(int m) {Base::resize(m,1); return *this;}
452  VectorBase& resizeKeep(int m) {Base::resizeKeep(m,1); return *this;}
453 
454  //TODO: this is not re-locking the number of columns at 1.
455  void clear() {Base::clear(); Base::resize(0,1);}
456 
457  ELT sum() const {ELT s; Base::getHelper().sum(reinterpret_cast<Scalar*>(&s)); return s; } // add all the elements
458  VectorIterator<ELT, VectorBase<ELT> > begin() {
459  return VectorIterator<ELT, VectorBase<ELT> >(*this, 0);
460  }
461  VectorIterator<ELT, VectorBase<ELT> > end() {
462  return VectorIterator<ELT, VectorBase<ELT> >(*this, size());
463  }
464 
465 protected:
466  // Create a VectorBase handle using a given helper rep.
467  explicit VectorBase(MatrixHelperRep<Scalar>* hrep) : Base(hrep) {}
468 
469 private:
470  // NO DATA MEMBERS ALLOWED
471 };
472 
473 } //namespace SimTK
474 
475 #endif // SimTK_SIMMATRIX_VECTORBASE_H_