Toast++  1.0.2 (r.539)
Forward and inverse modelling in optical tomography
vector.h
1 // -*-C++-*-
2 // ==========================================================================
3 // Module mathlib
4 // File vector.h
5 // Declaration of template class TVector ('template vector')
6 //
7 // The following typedefs for specific template types have been defined
8 // for convenience:
9 // RVector = TVector<double> ('real')
10 // FVector = TVector<float> ('float')
11 // CVector = TVector<complex> ('complex')
12 // IVector = TVector<int> ('integer')
13 // Vector = TVector<double> for backward compatibility
14 // ==========================================================================
15 
16 #ifndef __VECTOR_H
17 #define __VECTOR_H
18 
19 #include <string.h>
20 #include <malloc.h>
21 #include <iostream>
22 #include <sstream>
23 #include "mathdef.h"
24 #include "error.h"
25 //#include "complex.h"
26 //#include "scomplex.h"
27 #include <complex>
28 #include "fblas.h"
29 
30 #ifdef USE_CBLAS
31 #include <cblas++.h>
32 #endif // USE_CBLAS
33 
34 const int END = -1; // "end" index flag
35 
36 // ==========================================================================
37 // Nonmember declarations
38 
39 template<class VT> class TVector;
40 
41 template<class VT>
42 TVector<VT> inv (const TVector<VT> &v); // 1/v
43 
44 template<class VT>
45 TVector<VT> sqr (const TVector<VT> &v);
46 
47 template<class VT>
48 TVector<VT> sqrt (const TVector<VT> &v);
49 
50 template<class VT>
51 TVector<VT> log (const TVector<VT> &v);
52 
53 template<class VT>
54 TVector<VT> exp (const TVector<VT> &v);
55 
56 template<class VT>
57 TVector<VT> pow (const TVector<VT> &v, const VT &s);
58 
59 template<class VT>
60 TVector<VT> conj (const TVector<VT> &v);
61 
62 template<class VT>
63 TVector<VT> sin (const TVector<VT> &v);
64 
65 template<class VT>
66 TVector<VT> cos (const TVector<VT> &v);
67 
68 template<class VT>
69 TVector<VT> tan (const TVector<VT> &v);
70 
71 template<class VT>
72 TVector<VT> asin (const TVector<VT> &v);
73 
74 template<class VT>
75 TVector<VT> acos (const TVector<VT> &v);
76 
77 template<class VT>
78 TVector<VT> atan (const TVector<VT> &v);
79 
80 template<class VT>
81 TVector<VT> sinh (const TVector<VT> &v);
82 
83 template<class VT>
84 TVector<VT> cosh (const TVector<VT> &v);
85 
86 template<class VT>
87 TVector<VT> tanh (const TVector<VT> &v);
88 
89 template<class VT>
90 TVector<VT> asinh (const TVector<VT> &v);
91 
92 template<class VT>
93 TVector<VT> acosh (const TVector<VT> &v);
94 
95 template<class VT>
96 TVector<VT> atanh (const TVector<VT> &v);
97 
98 template<class VT>
99 VT dot (const TVector<VT> &v1, const TVector<VT> &v2);
100 
101 template<class VT>
102 VT doth (const TVector<VT> &v1, const TVector<VT> &v2);
103 
104 template<class VT>
105 VT vmin (const TVector<VT> &v);
106 
107 template<class VT>
108 VT vmax (const TVector<VT> &v);
109 
110 template<class VT>
111 TVector<VT> vsort (const TVector<VT> &v);
112 
113 template<class VT>
114 void vsort (const TVector<VT> &v, TVector<VT> &sorted_v, TVector<int> &sort_order);
115 
116 template<class VT>
117 VT sum (const TVector<VT> &v);
118 
119 template<class VT>
120 VT mean (const TVector<VT> &v);
121 
122 template<class VT>
123 VT median (const TVector<VT> &v);
124 
125 template<class VT>
126 VT variance (const TVector<VT> &v);
127 
128 template<class VT>
129 VT stdv (const TVector<VT> &v);
130 
131 template<class VT>
132 double l1norm (const TVector<VT> &v);
133 
134 template<class VT>
135 double l2norm (const TVector<VT> &v);
136 
137 template<class VT>
138 double linfnorm (const TVector<VT> &v);
139 
140 template<class VT>
141 double l2normsq (const TVector<VT> &v);
142 
143 template<class VT>
144 TVector<VT> &append (TVector<VT> &v1, const TVector<VT> &v2);
145 
146 template<class VT>
147 TVector<VT> cat (const TVector<VT> &v1, const TVector<VT> &v2);
148 
149 template<class VT>
150 TVector<VT> operator+ (const VT &s, const TVector<VT> &v);
151 
152 template<class VT>
153 TVector<VT> operator- (const VT &s, const TVector<VT> &v);
154 
155 template<class VT>
156 TVector<VT> operator* (const VT &s, const TVector<VT> &v);
157 
158 template<class VT>
159 TVector<VT> operator/ (const VT &s, const TVector<VT> &v);
160 
161 template<class VT>
162 bool operator== (const TVector<VT> &v1, const TVector<VT> &v2);
163 
164 template<class VT>
165 bool operator!= (const TVector<VT> &v1, const TVector<VT> &v2);
166 
167 template<class VT>
168 bool visnan (const TVector<VT> &v);
169 
170 template<class VT>
171 std::ostream &operator<< (std::ostream &os, const TVector<VT> &v);
172 
173 template<class VT>
174 std::istream &operator>> (std::istream &is, TVector<VT> &v);
175 
176 template<class VT>
177 VT SparseDotp (const TVector<VT> &v1, idxtype *idx1, int nidx1,
178  const TVector<VT> &v2, idxtype *idx2, int nidx2,
179  int from, int to);
180 
181 template<class VT>
183 
184 /* add specific function for CVector */
185 MATHLIB TVector<float> Re(const TVector<std::complex<float> > &vec);
186 MATHLIB TVector<double> Re(const TVector<std::complex<double> > &vec);
187 MATHLIB TVector<float> Im(const TVector<std::complex<float> > &vec);
188 MATHLIB TVector<double> Im(const TVector<std::complex<double> > &vec);
189 MATHLIB TVector<double> Mod(const TVector<std::complex<double> > &vec);
190 MATHLIB TVector<double> LogMod(const TVector<std::complex<double> > &vec);
191 MATHLIB TVector<double> Arg (const TVector<std::complex<double> > &vec);
192 MATHLIB TVector<std::complex<double> > Conj(const TVector<std::complex<double> > &vec); // obsolete
193 MATHLIB void SelfConj(const TVector<std::complex<double> > &vec);
194 MATHLIB TVector<std::complex<double> > Hadamard(const TVector<std::complex<double> > &a, const TVector<std::complex<double> > &b);
195 MATHLIB void SetReal (TVector<std::complex<double> > &z, const TVector<double> &zre);
196 MATHLIB void SetImag (TVector<std::complex<double> > &z, const TVector<double> &zim);
197 
198 MATHLIB TVector<std::complex<double> > MakeCVector (const TVector<std::complex<float> > &v);
199 
200 // ==========================================================================
201 // class TVector
202 
217 template<class VT> class TVector {
218 public:
222  TVector ();
223 
228  TVector (int dim);
229 
235  TVector (int dim, const VT s);
236 
249  TVector (int dim, VT *values, CopyMode cmode=DEEP_COPY);
250 
261  TVector (int dim, const char *init);
262 
267  TVector (const TVector<VT> &v);
268 
284  TVector (const TVector<VT> &v, int ofs, int dim);
285 
289  ~TVector () { Unlink (); }
290 
295  int Dim () const { return size; }
296 
302  VT &operator[] (int i) const;
303 
309  TVector &operator= (const TVector &v) { Copy (v); return *this; }
310 
312  TVector &operator= (VT s);
313 
318  void Copy (const TVector &v);
319 
331  void Copy (const TVector &v, int tofs, int sofs, int n);
332 
338  void New (int dim) { Unlink (); Allocate (dim); }
339 
345  void Clear ();
346 
353  void Relink (const TVector &v);
354 
365  void Relink (const TVector &v, int ofs, int n);
366 
375  void Relink (VT *values, int dim);
376 
383  VT *data_buffer () { return data; }
384 
391  const VT *data_buffer () const { return data; }
392 
399  bool Clip (VT vmin, VT vmax);
400 
407  void ShiftLeft (int n);
408 
415  void ShiftRight (int n);
416 
417  // arithmetic operators
418 
426  TVector operator+ (const TVector &v) const;
427 
435  TVector operator- (const TVector &v) const;
436 
442  TVector operator* (const TVector &v) const;
443 
449  TVector operator/ (const TVector &v) const;
450 
455  TVector operator+ (const VT &s) const;
456 
461  TVector operator- (const VT &s) const;
462 
467  TVector operator* (const VT &s) const;
468 
473  TVector operator/ (const VT &s) const;
474 
479  TVector operator- () const;
480 
488  friend bool operator== <> (const TVector<VT> &v1,
489  const TVector<VT> &v2);
490 
498  friend bool operator!=<> (const TVector<VT> &v1,
499  const TVector<VT> &v2);
500 
512  VT operator& (const TVector& v) const { return dot (*this, v); }
513 
514  // **** arithmetic/assignment operators ****
515 
521  TVector &operator+= (const TVector &v);
522 
528  TVector &operator-= (const TVector &v);
529 
535  TVector &operator*= (const TVector &v);
536 
542  TVector &operator/= (const TVector &v);
543 
549  TVector &operator+= (const VT &s);
550 
556  TVector &operator-= (const VT &s);
557 
563  TVector &operator*= (const VT &s);
564 
570  TVector &operator/= (const VT &s);
571 
572  // IO
573 
582  void Read (std::istream &is, int start = 0, int n = -1);
583 
595  void ReadIndexed (std::istream &is, int n = -1);
596 
597  // friends
599  friend TVector<VT> inv<> (const TVector<VT> &v);
600 
602  friend TVector<VT> sqr<> (const TVector<VT> &v);
603 
605  friend TVector<VT> sqrt<> (const TVector<VT> &v);
606 
608  friend TVector<VT> pow<> (const TVector<VT> &v, const VT &s);
609 
611  friend TVector<VT> log<> (const TVector<VT> &v);
612 
614  friend TVector<VT> exp<> (const TVector<VT> &v);
615 
618  friend TVector<VT> conj<> (const TVector<VT> &v);
619 
621  friend TVector<VT> sin<> (const TVector<VT> &v);
622 
624  friend TVector<VT> cos<> (const TVector<VT> &v);
625 
627  friend TVector<VT> tan<> (const TVector<VT> &v);
628 
630  friend TVector<VT> asin<> (const TVector<VT> &v);
631 
633  friend TVector<VT> acos<> (const TVector<VT> &v);
634 
636  friend TVector<VT> atan<> (const TVector<VT> &v);
637 
639  friend TVector<VT> sinh<> (const TVector<VT> &v);
640 
642  friend TVector<VT> cosh<> (const TVector<VT> &v);
643 
645  friend TVector<VT> tanh<> (const TVector<VT> &v);
646 
649  friend TVector<VT> asinh<> (const TVector<VT> &v);
650 
653  friend TVector<VT> acosh<> (const TVector<VT> &v);
654 
664  friend TVector<VT> atanh<> (const TVector<VT> &v);
665 
678  friend VT dot<> (const TVector<VT> &v1, const TVector<VT> &v2);
679 
692  friend VT doth<> (const TVector<VT> &v1, const TVector<VT> &v2);
693 
700  friend VT vmin<> (const TVector<VT> &v);
701 
708  friend VT vmax<> (const TVector<VT> &v);
709 
718  friend TVector<VT> vsort<> (const TVector<VT> &v);
719 
730  friend void vsort<> (const TVector<VT> &v, TVector<VT> &sorted_v, TVector<int> &sort_order);
731 
736  friend VT sum<> (const TVector<VT> &v);
737 
745  friend VT mean<> (const TVector<VT> &v);
746 
757  friend VT median<> (const TVector<VT> &v);
758 
764  friend VT variance<> (const TVector<VT> &v);
765 
771  friend VT stdv<> (const TVector<VT> &v);
772 
778  friend double l1norm<> (const TVector<VT> &v);
779 
785  friend double l2norm<> (const TVector<VT> &v);
786 
792  friend double linfnorm<> (const TVector<VT> &v);
793 
799  friend double l2normsq<> (const TVector<VT> &v);
800 
806  friend double length (const TVector<VT> &vec)
807  { return l2norm (vec); }
808 
817  friend TVector<VT> &append<> (TVector<VT> &v1, const TVector<VT> &v2);
818 
826  friend TVector<VT> cat<> (const TVector<VT> &v1, const TVector<VT> &v2);
827 
840 #ifndef USE_CUDA_FLOAT
841 // none of the operator friend definitions seem to work with nvcc
842 #if defined(_WIN32)
843  template <class VT> friend TVector<VT> operator+ (const VT &s,
844  const TVector<VT> &v);
845 #elif (GCC_VERSION < 30404) // old-style definitions
846  friend TVector<VT> operator+ <> (const VT &s, const TVector<VT> &v);
847 #else
848  friend TVector<VT> (::operator+ <>) (const VT &s, const TVector<VT> &v);
849 #endif
850 #endif // !USE_CUDA_FLOAT
851 
864 #if defined(_WIN32)
865  template<class VT> friend TVector<VT> operator- (const VT &s,
866  const TVector<VT> &v);
867 #elif (GCC_VERSION < 30404) // old-style definitions
868  friend TVector<VT> operator- <> (const VT &s, const TVector<VT> &v);
869 #else
870  friend TVector<VT> (::operator- <>) (const VT &s, const TVector<VT> &v);
871 #endif
872 
885 #if defined(_WIN32)
886  template<class VT> friend TVector<VT> operator* (const VT &s,
887  const TVector<VT> &v);
888 #elif (GCC_VERSION < 30404) // old-style definitions
889  friend TVector<VT> operator* <> (const VT &s, const TVector<VT> &v);
890 #else
891  friend TVector<VT> (::operator* <>) (const VT &s, const TVector<VT> &v);
892 #endif
893 
906 #if defined(_WIN32)
907  template<class VT> friend TVector<VT> operator/ (const VT &s,
908  const TVector<VT> &v);
909  // JK also works (TVector::operator/ )
910 #elif (GCC_VERSION < 30404) // old-style definitions
911  friend TVector<VT> operator/ <> (const VT &s, const TVector<VT> &v);
912 #else
913  friend TVector<VT> (::operator/ <>) (const VT &s, const TVector<VT> &v);
914 #endif
915 
916  friend bool visnan<> (const TVector<VT> &v);
917 
935  friend std::ostream &operator<< <> (std::ostream &os,
936  const TVector<VT> &v);
937 
956  friend std::istream &operator>> <> (std::istream &is, TVector<VT> &v);
957 
975  friend VT SparseDotp<> (const TVector<VT> &v1, idxtype *idx1, int nidx1,
976  const TVector<VT> &v2, idxtype *idx2, int nidx2,
977  int from, int to);
978 
988  friend TVector<double> UnfoldComplex<> (const TVector<VT> &v);
989 
997  void Scan (const char *cbuf, int nmax = 0);
998 
999 protected:
1008  void Allocate (int dim);
1009 
1019  void Link (const TVector<VT> &vec);
1020 
1033  void Link (const TVector<VT> &vec, int ofs, int dim);
1034 
1042  void Link (VT *values, int dim);
1043 
1047  void Unlink ();
1048 
1049  int size;
1050  int *base_size;
1051  int *base_nref;
1052  VT *data;
1054 
1055 private:
1056  bool ext_data;
1057 };
1058 
1059 // ==========================================================================
1060 // template typedefs
1061 
1062 typedef TVector<double> RVector;
1063 typedef TVector<float> FVector;
1066 typedef TVector<int> IVector;
1067 
1068 #endif // !__VECTOR_H
void Copy(const TVector &v)
Vector copy. Replaces the vector with a copy of 'v'.
int * base_size
pointer to length of linked data block
Definition: vector.h:1050
const VT * data_buffer() const
Direct access to vector data array.
Definition: vector.h:391
void ShiftRight(int n)
Move vector elements right.
VT * data_buffer()
Direct access to vector data array.
Definition: vector.h:383
friend TVector< VT > acosh(const TVector< VT > &v)
Returns element-wise inverse hyperbolic cosine vector acosh((*this)[i])
VT * base_data
pointer to linked data block
Definition: vector.h:1053
bool Clip(VT vmin, VT vmax)
Truncate vector elements to specified range.
void New(int dim)
Resize the vector.
Definition: vector.h:338
friend TVector< VT > cat(const TVector< VT > &v1, const TVector< VT > &v2)
Concatenate two vectors.
Templated vector class.
Definition: vector.h:39
TVector operator+(const TVector &v) const
Element-wise addition of vectors.
TVector & operator+=(const TVector &v)
Element-wise addition/assignment of vectors.
friend double l1norm(const TVector< VT > &v)
L1 norm of vector elements.
friend TVector< VT > inv(const TVector< VT > &v)
Returns element-wise inverse vector 1/(*this)[i].
friend TVector< VT > pow(const TVector< VT > &v, const VT &s)
Returns element-wise power vector (*this)[i]^s.
friend std::istream & operator>>(std::istream &is, TVector< VT > &v)
Read vector from input stream.
friend double l2normsq(const TVector< VT > &v)
Square of L2 norm of vector elements.
void ReadIndexed(std::istream &is, int n=-1)
Read indexed vector elements from stream.
TVector()
Constructor. Creates a vector of length 0.
friend TVector< VT > & append(TVector< VT > &v1, const TVector< VT > &v2)
Concatenate two vectors.
TVector operator*(const TVector &v) const
Element-wise multiplication of vectors.
friend VT median(const TVector< VT > &v)
Median value of elements.
friend TVector< VT > exp(const TVector< VT > &v)
Returns element-wise base-e exponential vector exp((*this)[i])
friend double length(const TVector< VT > &vec)
Length of a vector (L2 norm of elements)
Definition: vector.h:806
TVector operator-() const
Unary minus.
friend TVector< VT > cos(const TVector< VT > &v)
Returns element-wise cosine vector cos((*this)[i])
friend TVector< double > UnfoldComplex(const TVector< VT > &v)
Unfold real and imaginary parts of a complex vector.
void Relink(const TVector &v)
Link the vector to the data block of another vector.
friend double linfnorm(const TVector< VT > &v)
L-infinity norm of vector elements.
int * base_nref
pointer to data block reference counter
Definition: vector.h:1051
TVector & operator*=(const TVector &v)
Element-wise multiplication/assignment of vectors.
friend bool operator!=(const TVector< VT > &v1, const TVector< VT > &v2)
Vector comparison (relational operator)
int size
vector length
Definition: vector.h:1049
friend TVector< VT > acos(const TVector< VT > &v)
Returns element-wise arc cosine vector acos((*this)[i])
VT operator&(const TVector &v) const
Dot product of two vectors.
Definition: vector.h:512
friend TVector< VT > tan(const TVector< VT > &v)
Returns element-wise tangent vector tan((*this)[i])
TVector & operator-=(const TVector &v)
Element-wise subtraction/assignment of vectors.
void Scan(const char *cbuf, int nmax=0)
Initialise vector elements from a string.
friend TVector< VT > sinh(const TVector< VT > &v)
Returns element-wise hyperbolic sine vector sinh((*this)[i])
friend TVector< VT > tanh(const TVector< VT > &v)
Returns element-wise hyperbolic tangent vector tanh((*this)[i])
int Dim() const
Returns the size of the vector.
Definition: vector.h:295
friend VT SparseDotp(const TVector< VT > &v1, idxtype *idx1, int nidx1, const TVector< VT > &v2, idxtype *idx2, int nidx2, int from, int to)
Sparse dot product of two vectors.
~TVector()
Destructor. Delete vector and deallocate data block.
Definition: vector.h:289
friend TVector< VT > log(const TVector< VT > &v)
Returns element-wise natural logarithm vector ln((*this)[i])
friend VT mean(const TVector< VT > &v)
Mean value of elements.
void Clear()
Zeroes all elements.
void Unlink()
Unlink vector from its data block.
TVector & operator=(const TVector &v)
Assignment operator.
Definition: vector.h:309
friend TVector< VT > sin(const TVector< VT > &v)
Returns element-wise sine vector sin((*this)[i])
friend TVector< VT > vsort(const TVector< VT > &v)
Sort vector.
VT * data
pointer to first vector element
Definition: vector.h:1052
friend TVector< VT > asinh(const TVector< VT > &v)
Returns element-wise inverse hyperbolic sine vector asinh((*this)[i])
void Allocate(int dim)
Allocate data buffer.
void Read(std::istream &is, int start=0, int n=-1)
Read vector from stream.
friend VT stdv(const TVector< VT > &v)
Standard deviation of element values.
friend TVector< VT > asin(const TVector< VT > &v)
Returns element-wise arc sine vector asin((*this)[i])
friend TVector< VT > atan(const TVector< VT > &v)
Returns element-wise arc tangent vector atan((*this)[i])
TVector operator/(const TVector &v) const
Element-wise division of vectors.
void ShiftLeft(int n)
Move vector elements left.
friend VT dot(const TVector< VT > &v1, const TVector< VT > &v2)
Dot product of two vectors.
friend TVector< VT > sqr(const TVector< VT > &v)
Returns element-wise square vector (*this)[i]^2.
friend bool operator==(const TVector< VT > &v1, const TVector< VT > &v2)
Vector comparison (relational operator)
friend TVector< VT > sqrt(const TVector< VT > &v)
Returns element-wise square-root vector (*this)[i]^(1/2)
friend VT doth(const TVector< VT > &v1, const TVector< VT > &v2)
Hermitian product of two vectors.
friend TVector< VT > cosh(const TVector< VT > &v)
Returns element-wise hyperbolic cosine vector cosh((*this)[i])
void Link(const TVector< VT > &vec)
Link the vector to the data block of another vector.
friend double l2norm(const TVector< VT > &v)
L2 norm of vector elements.
friend TVector< VT > atanh(const TVector< VT > &v)
Returns element-wise inverse hyperbolic tangent vector atanh(v[i])
friend VT sum(const TVector< VT > &v)
Sum of elements.
TVector & operator/=(const TVector &v)
Element-wise division/assignment of vectors.
friend VT variance(const TVector< VT > &v)
Variance of element values.
friend VT vmax(const TVector< VT > &v)
Extract largest element in vector.
VT & operator[](int i) const
Vector element access operator (read and write)
friend TVector< VT > conj(const TVector< VT > &v)
Returns element-wise complex conjugate vector ((*this([i].re, -(*this)[i].im)
friend VT vmin(const TVector< VT > &v)
Extract smallest element in vector.