TMB Documentation  v1.9.1
array.hpp
1 // Copyright (C) 2013-2015 Kasper Kristensen
2 // License: GPL-2
3 
21 template<class Type>
22 struct array:Map< Array<Type,Dynamic,1> >{
23 
24  typedef Array<Type,Dynamic,1> Base;
25  typedef Map< Base > MapBase;
26 
27  vector<int> dim;
28 
29  /* def: inner product sum(mult*tuple) gives vector index
30  example: array a(4,5,6);
31  dim=[4,5,6] first dimension runs fastest, last dimension slowest:
32  mult=[1,1*4,1*4*5]
33  */
34  vector<int> mult;
35 
36  Base vectorcopy; /* Array data */
37 
39  void setdim(vector<int> dim_){
40  dim=dim_;
41  mult.resize(dim.size());
42  mult[0]=1;
43  for(int k=1;k<dim.size();k++){
44  mult[k]=mult[k-1]*dim[k-1];
45  }
46  }
47  void initZeroArray(vector<int> dim_){
48  vectorcopy.resize(dim_.prod());
49  vectorcopy.setZero();
50  if (vectorcopy.size() > 0) {
51  new (this) MapBase(&vectorcopy[0],vectorcopy.size()); /* Eigen manual: Despite appearances, this does not invoke the memory allocator... */
52  }
53  setdim(dim_);
54  }
55 
56  /* Utility function */
57  vector<int> c(int n1){
58  vector<int> ans(1);
59  ans << n1;
60  return ans;
61  }
62  vector<int> c(int n1, int n2){
63  vector<int> ans(2);
64  ans << n1,n2;
65  return ans;
66  }
67  vector<int> c(int n1, int n2, int n3){
68  vector<int> ans(3);
69  ans << n1,n2,n3;
70  return ans;
71  }
72  vector<int> c(int n1, int n2, int n3, int n4){
73  vector<int> ans(4);
74  ans << n1,n2,n3,n4;
75  return ans;
76  }
77  vector<int> c(int n1, int n2, int n3, int n4, int n5){
78  vector<int> ans(5);
79  ans << n1,n2,n3,n4,n5;
80  return ans;
81  }
82  vector<int> c(int n1, int n2, int n3, int n4, int n5, int n6){
83  vector<int> ans(6);
84  ans << n1,n2,n3,n4,n5,n6;
85  return ans;
86  }
87  vector<int> c(int n1, int n2, int n3, int n4, int n5, int n6, int n7){
88  vector<int> ans(7);
89  ans << n1,n2,n3,n4,n5,n6,n7;
90  return ans;
91  }
92 
93 
94  /* Default constructor: e.g. array<double> a; */
95  array():MapBase(NULL,0){};
96 
98  array(vector<int> dim_):MapBase(NULL,0){
99  initZeroArray(dim_);
100  }
101  array(int n1):MapBase(NULL,0){initZeroArray(c(n1));}
102  array(int n1, int n2):MapBase(NULL,0){initZeroArray(c(n1,n2));}
103  array(int n1, int n2, int n3):MapBase(NULL,0){initZeroArray(c(n1,n2,n3));}
104  array(int n1, int n2, int n3, int n4):MapBase(NULL,0){initZeroArray(c(n1,n2,n3,n4));}
105  array(int n1, int n2, int n3, int n4, int n5):MapBase(NULL,0){initZeroArray(c(n1,n2,n3,n4,n5));}
106  array(int n1, int n2, int n3, int n4, int n5, int n6):MapBase(NULL,0){initZeroArray(c(n1,n2,n3,n4,n5,n6));}
107  array(int n1, int n2, int n3, int n4, int n5, int n6, int n7):MapBase(NULL,0){initZeroArray(c(n1,n2,n3,n4,n5,n6,n7));}
108 
109  /* Default construction: always deep copy ! */
110  template<class T>
111  array(T &x, vector<int> dim_):MapBase(NULL,0),vectorcopy(x){
112  if(x.size()>0){
113  new (this) MapBase(&vectorcopy[0],x.size()); /* Eigen manual: Despite appearances, this does not invoke the memory allocator... */
114  }
115  setdim(dim_);
116  }
117  /* Deep copy existing array - as above with dim_=x.dim */
118  template<class T>
119  array(array<T> &x):MapBase(NULL,0),vectorcopy(x){
120  if(x.size()>0){
121  new (this) MapBase(&vectorcopy[0],x.size()); /* Eigen manual: Despite appearances, this does not invoke the memory allocator... */
122  }
123  setdim(x.dim);
124  }
125  /* Templated CTOR above may be ignored by the
126  compiler. E.g. clang-13 used a compiler generated default copy
127  CTOR which did not synchronize the MapBase::data() pointer (!) */
128  array(const array &x) : MapBase(NULL,0), vectorcopy(x) {
129  if ( x.size() > 0 ) {
130  // sync pointer
131  new (this) MapBase(&vectorcopy[0], x.size());
132  }
133  setdim(x.dim);
134  }
135 
136  /* Sometimes we want a reference only... See col() */
137  array(Type *p, vector<int> dim_):MapBase(p,dim_.prod()){
138  setdim(dim_);
139  }
140 
141  /*
142  Example:
143  array<double> a(3,4);
144  array<double> c=a; // Assign
145  a.col(1)=c.col(1); // Assign
146  array<double> b=a.col(1); // Initialize (not this method)
147  */
148  array & operator= (const array & other)
149  {
150  if(this->dim.size() == 0){ // Not initialized (default constructed)
151  this->initZeroArray(other.dim);
152  }
153  this->MapBase::operator=(other);
154  this->setdim(other.dim);
155  return *this;
156  }
157 
158  /* Assign array from any other object. First convert other object to
159  1D object (because we assign to underlying 1D array data).
160 
161  Example:
162  a = m (assign from matrix)
163  a.col(0) = m; (assign from matrix)
164  a = m*m (assign from unevaluated expression template)
165  a = v (assign from vector)
166  */
167  template <class T>
168  array<Type> operator=(T y){
169  Array<Type, Dynamic, Dynamic> a = y;
170  a.resize(a.size(),1);
171  return array(MapBase::operator=(a), dim);
172  }
173 
174  void print(){
175  std::cout << "Array dim: ";
176  for(int i=0;i<dim.size();i++)std::cout << dim[i] << " ";
177  std::cout << "\n";
178  std::cout << "Array val: ";
179  for(int i=0;i<this->MapBase::size();i++)std::cout << this->MapBase::operator[](i) << " ";
180  std::cout << "\n";
181  };
182 
187  int cols(){
188  return dim[dim.size()-1];
189  }
190 
194  int rows(){
195  return dim[0];
196  }
197 
201  array<Type> col(int i){
202  int nslice=this->MapBase::size()/this->cols();
203  Type* p=&(this->MapBase::operator()(i*nslice));
204  vector<int> newdim;
205  if(dim.size()>1){
206  newdim=dim.segment(0,dim.size()-1);
207  } else {
208  newdim.resize(1);
209  newdim << 1;
210  }
211  return array(p,newdim);
212  }
213 
216  Type& operator()(int i1){
217  return this->MapBase::operator[](i1);
218  }
220  Type& operator()(int i1, int i2){
221  return this->MapBase::operator[](index(c(i1,i2)));
222  }
224  Type& operator()(int i1, int i2, int i3){
225  return this->MapBase::operator[](index(c(i1,i2,i3)));
226  }
228  Type& operator()(int i1, int i2, int i3, int i4){
229  return this->MapBase::operator[](index(c(i1,i2,i3,i4)));
230  }
232  Type& operator()(int i1, int i2, int i3, int i4, int i5){
233  return this->MapBase::operator[](index(c(i1,i2,i3,i4,i5)));
234  }
236  Type& operator()(int i1, int i2, int i3, int i4, int i5, int i6){
237  return this->MapBase::operator[](index(c(i1,i2,i3,i4,i5,i6)));
238  }
240  Type& operator()(int i1, int i2, int i3, int i4, int i5, int i6, int i7){
241  return this->MapBase::operator[](index(c(i1,i2,i3,i4,i5,i6,i7)));
242  }
243 
244 
245  int index(vector<int> tup){
246  eigen_assert( tup.size() == dim.size() );
247  eigen_assert( ( (dim*0 <= tup) && (tup < dim) ).all() );
248  return (tup*mult).sum();
249  }
250  vector<int> tuple(int i){
251  vector<int> revtup(dim.size());
252  vector<int> revmult=mult.reverse();
253  for(int k=0;k<dim.size();k++){
254  revtup[k]=i/revmult[k];
255  i=i-revtup[k]*revmult[k];
256  }
257  return revtup.reverse();
258  }
259 
264  vector<Type> x(this->size());
265  array<Type> ans(x,dim(p)); /* Create new array with permuted dimension */
266  for(int i=0;i<this->size();i++){ /* Loop through values of old array */
267  ans[ans.index(tuple(i)(p))]=this->operator[](i);
268  }
269  return ans;
270  }
277  vector<int> p(dim.size());
278  for(int i=0;i<p.size();i++)p[i]=i;
279  return this->perm(p.reverse());
280  }
281 
282  int mod(int i,int n){return ((i%n)+n)%n;}
291  vector<int> p(dim.size());
292  for(int i=0;i<p.size();i++)p[i]=mod(i-n,p.size());
293  return this->perm(p);
294  }
295 
296 #define INHERIT(OP) \
297  template <class T> \
298  array<Type> OP(T y){return array(MapBase::OP(y),dim);}
299  INHERIT(operator+)
300  INHERIT(operator-)
301  INHERIT(operator*)
302  INHERIT(operator/)
303 #undef INHERIT
304 
311  tmbutils::matrix<Type> ans = this->MapBase::matrix();
312  ans.resize(this->rows(), ans.size() / this->rows() );
313  return ans;
314  }
315 
317  tmbutils::vector<Type> vec() { return *this; }
318 
319  /* Methods this class should *not* inherit (generate compile time error if used) */
320  private:
321  using MapBase::row;
322 };
323 
Vector class used by TMB.
Definition: vector.hpp:17
void setdim(vector< int > dim_)
Sets dimension attribute and updates internal stride. Must be used when e.g. collapsing array dimensi...
Definition: array.hpp:39
array< Type > perm(vector< int > p)
Array permutation. Permutes array dimensions corresponding to permutation vector p.
Definition: array.hpp:263
tmbutils::matrix< Type > matrix()
Convert TMB array to matrix by keeping the first dimension and collapsing remaining dimensions...
Definition: array.hpp:310
int cols()
Number of outer-most dimensions.
Definition: array.hpp:187
Type & operator()(int i1)
Elementwise subsetting 1D array. Also allowed in general to access the underlying vector of n-dim arr...
Definition: array.hpp:216
array< Type > rotate(int n)
Array rotate (Special case of array permutation)
Definition: array.hpp:290
Array class used by TMB.
Definition: array.hpp:22
int rows()
Number of inner-most dimensions.
Definition: array.hpp:194
tmbutils::vector< Type > vec()
Convert TMB array to vector.
Definition: array.hpp:317
Matrix class used by TMB.
Definition: tmbutils.hpp:102
array< Type > transpose()
Array transpose (Special case of array permutation)
Definition: array.hpp:276
Vector class used by TMB.
Definition: tmbutils.hpp:18
array< Type > col(int i)
Extract sub-array with write access Index i refers to the outer-most (i.e. final) dimension...
Definition: array.hpp:201
array(vector< int > dim_)
Construct array from dimension vector and fill with zeros.
Definition: array.hpp:98
License: GPL v2