KADATH
array.hpp
1 /*
2  Copyright 2017 Philippe Grandclement
3 
4  This file is part of Kadath.
5 
6  Kadath is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  Kadath is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with Kadath. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #ifndef __ARRAY_HPP_
21 #define __ARRAY_HPP_
22 
23 
24 #include "headcpp.hpp"
25 #include "utilities.hpp"
26 #include "dim_array.hpp"
27 #include "index.hpp"
28 #include "array_iterator.hpp"
29 
30 namespace Kadath {
31 
35 enum Array_ordering : bool {
39  first_index,
43  last_index
44 };
45 
46 template <typename T> class Array ;
47 template <typename T> ostream& operator<< (ostream&, const Array<T>& ) ;
48 template <typename T> Array<T> sin(const Array<T>& ) ;
49 template <typename T> Array<T> cos(const Array<T>& ) ;
50 template <typename T> Array<T> sinh(const Array<T>& ) ;
51 template <typename T> Array<T> cosh(const Array<T>& ) ;
52 template <typename T> Array<T> operator+ (const Array<T>&) ;
53 template <typename T> Array<T> operator- (const Array<T>&) ;
54 template <typename T> Array<T> operator+ (const Array<T>&, const Array<T>&) ;
55 template <typename T> Array<T> operator+ (const Array<T>&, T) ;
56 template <typename T> Array<T> operator+ (T, const Array<T>&) ;
57 template <typename T> Array<T> operator- (const Array<T>&, const Array<T>&) ;
58 template <typename T> Array<T> operator- (const Array<T>&, T) ;
59 template <typename T> Array<T> operator- (T, const Array<T>&) ;
60 template <typename T> Array<T> operator* (const Array<T>&, const Array<T>&) ;
61 template <typename T> Array<T> operator* (const Array<T>&, T) ;
62 template <typename T> Array<T> operator* (T, const Array<T>&) ;
63 template <typename T> Array<T> operator/ (const Array<T>&, const Array<T>&) ;
64 template <typename T> Array<T> operator/ (const Array<T>&, T) ;
65 template <typename T> Array<T> operator/ (T, const Array<T>&) ;
66 template <typename T> Array<T> pow (const Array<T>&, int) ;
67 template <typename T> Array<T> pow (const Array<T>&, double) ;
68 template <typename T> Array<T> sqrt (const Array<T>&) ;
69 template <typename T> Array<T> exp (const Array<T>&) ;
70 template <typename T> Array<T> log (const Array<T>&) ;
71 template <typename T> Array<T> atanh (const Array<T>&) ;
72 template <typename T> Array<T> fabs (const Array<T>&) ;
73 template <typename T> T scal (const Array<T>&, const Array<T>&) ;
74 template <typename T> T diffmax (const Array<T>&, const Array<T>&) ;
75 template <typename T> T max (const Array<T>&) ;
76 template <typename T> T min (const Array<T>&) ;
77 template <typename T> T sum (const Array<T>&) ;
78 template <typename T> Array<T> atan (const Array<T>&) ;
79 
86 template <typename T> class Array : public Memory_mapped {
87  static_assert(std::is_arithmetic<T>::value,"Array<T> implementation won't be efficient if T is not a "
88  " either an integral or floating point built-in type.");
89  public:
90  Dim_array dimensions ;
91  int nbr ;
92  Memory_mapped_array<T> data ;
93 
94  public:
96  using reference = T& ;
98  using const_reference = T const &;
100  using pointer = T*;
102  using const_pointer = T const *;
107  explicit Array (const Dim_array& res) : dimensions{res}, nbr{1}, data{nullptr} {
108  for (int i=0 ; i<res.get_ndim() ; i++) nbr *= res(i) ;
109  data.resize(nbr) ;
110  }
116  explicit Array (int i) : dimensions{1}, nbr{i}, data{nbr} {
117  dimensions.set(0) = i ;
118  }
125  explicit Array (int i, int j) : dimensions{2}, nbr{i*j}, data{nbr} {
126  dimensions.set(0) = i ;
127  dimensions.set(1) = j ;
128  }
129 
137  explicit Array (int i, int j, int k) : dimensions{3},nbr{i*j*k},data{nbr} {
138  dimensions.set(0) = i ;
139  dimensions.set(1) = j ;
140  dimensions.set(2) = k ;
141  }
148  explicit Array (FILE* fd,Array_ordering order = last_index) ;
149 
150  public:
152  void swap(Array<T> & so) {dimensions.swap(so.dimensions); std::swap(nbr,so.nbr); data.swap(so.data);}
153 
154  void delete_data() {data.clear();}
155 
160  void resize(Dim_array const & new_dim) {this->operator=(std::move(Array<T>{new_dim}));}
165  void resize(int new_size) {this->operator=(std::move(Array<T>{new_size}));}
166 
176  void save (FILE* fd, Array_ordering order = last_index) const ;
181  Array & operator= (T xx) {for (int i=0 ; i<nbr ; i++) data[i] = xx ; return *this;}
186  reference set(const Index& pos) {
187  assert (pos.sizes == dimensions) ;
188  int index = pos(dimensions.get_ndim()-1) ;
189  for (int i=dimensions.get_ndim()-2 ; i>=0 ; i--) {
190  index *= dimensions(i) ;
191  assert ((pos(i) >=0) && (pos(i)<dimensions(i))) ;
192  index += pos(i) ;
193  }
194  return data[index] ;
195  }
196 // reference set(const Index& pos) {
198 // int index = pos(0) ;
199 // for (int i=1 ; i<dimensions.get_ndim() ; i++) {
200 // index *= dimensions(i) ;
202 // index += pos(i) ;
203 // }
204 // return data[index] ;
205 // }
210  reference set(const Array_iterator &pos) {return data[pos.position];}
215  reference set(int i) {
216  assert (dimensions.get_ndim() == 1) ;
217  assert ((i >=0) && (i<dimensions(0))) ;
218  return data[i] ;
219  }
225  reference set(int i, int j) {
226  assert (dimensions.get_ndim() == 2) ;
227  assert ((i >=0) && (i<dimensions(0))) ;
228  assert ((j >=0) && (j<dimensions(1))) ;
229  return data[i+j*dimensions(0)] ;
230 // return data[i*dimensions(1)+j] ;
231  }
238  reference set(int i, int j, int k) {
239  assert (dimensions.get_ndim() == 3) ;
240  assert ((i >=0) && (i<dimensions(0))) ;
241  assert ((j >=0) && (j<dimensions(1))) ;
242  assert ((k >=0) && (k<dimensions(2))) ;
243  return data[i+dimensions(0)*(j+k*dimensions(1))] ;
244 // return data[(i*dimensions(1)+j)*dimensions(2)+k] ;
245  }
250  T operator() (const Index& pos) const {
251  assert (pos.sizes == dimensions) ;
252  int index = pos(dimensions.get_ndim()-1) ;
253  for (int i=dimensions.get_ndim()-2 ; i>=0 ; i--) {
254  index *= dimensions(i) ;
255  assert ((pos(i) >=0) && (pos(i)<dimensions(i))) ;
256  index += pos(i) ;
257  }
258  return data[index] ;
259  }
260 // T operator() (const Index& pos) const {
262 // int index = pos(0) ;
263 // for (int i=1 ; i<dimensions.get_ndim() ; i++) {
264 // index *= dimensions(i) ;
266 // index += pos(i) ;
267 // }
268 // return data[index] ;
269 // }
274  T operator()(Array_iterator const &pos) const {return data[pos.position];}
279  T operator() (int i) const {
280  assert (dimensions.get_ndim() ==1) ;
281  assert ((i >=0) && (i<dimensions(0))) ;
282  return data[i] ;
283  }
289  T operator() (int i,int j) const {
290  assert (dimensions.get_ndim() ==2) ;
291  assert ((i >=0) && (i<dimensions(0))) ;
292  assert ((j >=0) && (j<dimensions(1))) ;
293  return data[i+j*dimensions(0)] ;
294 // return data[i*dimensions(1)+j] ;
295  }
302  T operator() (int i,int j, int k) const {
303  assert (dimensions.get_ndim() ==3) ;
304  assert ((i >=0) && (i<dimensions(0))) ;
305  assert ((j >=0) && (j<dimensions(1))) ;
306  assert ((k >=0) && (k<dimensions(2))) ;
307  return data[i+dimensions(0)*(j+k*dimensions(1))] ;
308 // return data[i*dimensions(1)*dimensions(2)+j*dimensions(2)+k] ;
309  }
313  CXX_17_ATTRIBUTES(nodiscard) const_pointer get_data() const {return data.get_data() ;} ;
314 
318  CXX_17_ATTRIBUTES(nodiscard) pointer set_data() {return data.set_data() ;} ;
319 
323  CXX_17_ATTRIBUTES(nodiscard) int get_ndim() const {return dimensions.get_ndim() ;} ;
327  CXX_17_ATTRIBUTES(nodiscard) int get_nbr() const {return nbr ;} ;
331  CXX_17_ATTRIBUTES(nodiscard) int get_size(int i) const {return dimensions(i) ;} ;
335  CXX_17_ATTRIBUTES(nodiscard) const Dim_array& get_dimensions() const {return dimensions ;} ;
336 
340  bool is_increasing() const ;
341 
343  Array<T> & operator+= (const Array<T>& so)
344  {assert(nbr == so.nbr); for(int i{0};i<nbr;i++) data[i] += so.data[i]; return *this;}
346  Array<T> & operator-= (const Array<T>& so)
347  {assert(nbr == so.nbr); for(int i{0};i<nbr;i++) data[i] -= so.data[i]; return *this;}
349  Array<T> & operator*= (const Array<T>& so)
350  {assert(nbr == so.nbr); for(int i{0};i<nbr;i++) data[i] *= so.data[i]; return *this;}
352  Array<T> & operator/= (const Array<T>& so)
353  {assert(nbr == so.nbr); for(int i{0};i<nbr;i++) data[i] /= so.data[i]; return *this;}
355  Array<T> & operator+= (const T xx) {for(int i{0};i<nbr;i++) data[i] += xx; return *this;}
357  Array<T> & operator-= (const T xx) {for(int i{0};i<nbr;i++) data[i] -= xx; return *this;}
359  Array<T> & operator*= (const T xx) {for(int i{0};i<nbr;i++) data[i] *= xx; return *this;}
361  Array<T> & operator/= (const T xx) {for(int i{0};i<nbr;i++) data[i] /= xx; return *this;}
362 
368  int to_last_dim_major_index(int i_first_dim_major) const;
374  int to_first_dim_major_index(int i_last_dim_major) const;
375 
376  friend class Matrice ;
377 
379  friend ostream& operator<< <> (ostream&, const Array<T>& ) ;
381  friend Array<T> sin<> (const Array<T>&) ;
383  friend Array<T> cos<> (const Array<T>&) ;
385  friend Array<T> sinh<> (const Array<T>&) ;
387  friend Array<T> cosh<> (const Array<T>&) ;
389  friend Array<T> operator+ <>(const Array<T>&) ;
391  friend Array<T> operator- <> (const Array<T>&) ;
393  friend Array<T> operator+ <> (const Array<T>&, const Array<T>&) ;
395  friend Array<T> operator+ <>(const Array<T>&, T) ;
397  friend Array<T> operator+ <>(T, const Array<T>&) ;
399  friend Array<T> operator- <>(const Array<T>&, const Array<T>&) ;
401  friend Array<T> operator- <>(const Array<T>&, T) ;
403  friend Array<T> operator- <>(T, const Array<T>&) ;
405  friend Array<T> operator* <>(const Array<T>&, const Array<T>&) ;
407  friend Array<T> operator* <>(const Array<T>&, T) ;
409  friend Array<T> operator* <>(T, const Array<T>&) ;
411  friend Array<T> operator/ <>(const Array<T>&, const Array<T>&) ;
413  friend Array<T> operator/ <>(const Array<T>&, T) ;
415  friend Array<T> operator/ <>(T, const Array<T>&) ;
417  friend Array<T> pow <>(const Array<T>&, int) ;
419  friend Array<T> pow <>(const Array<T>&, double) ;
421  friend Array<T> sqrt <>(const Array<T>&) ;
423  friend Array<T> exp <>(const Array<T>&) ;
425  friend Array<T> log <>(const Array<T>&) ;
427  friend Array<T> atanh <>(const Array<T>&) ;
429  friend Array<T> fabs <>(const Array<T>&) ;
431  friend T scal<>(const Array<T>&, const Array<T>&) ;
433  friend T diffmax<>(const Array<T>&, const Array<T>&) ;
435  friend T max<> (const Array<T>&) ;
437  friend T min<> (const Array<T>&) ;
439  friend T sum<> (const Array<T>&) ;
441  friend Array<T> atan <>(const Array<T>&) ;
442 } ;
443 
444 template <typename T> bool Array<T>::is_increasing() const {
445  assert(get_ndim()==1) ;
446  for (int i = 1 ; i < nbr ; i++){
447  if (data[i]-data[i-1] <= 0){
448  return false ;
449  }
450  }
451  return true ;
452  }
453 
454  template<typename T> int Array<T>::to_last_dim_major_index(int i_first_dim_major) const {
455  int const dim {dimensions.get_ndim()};
456  int k{0};
457  Index ret{dimensions};
458  int q=i_first_dim_major;
459  while(q!=0 && k<dim) {
460  auto q_r = div(q,dimensions(k));
461  q = q_r.quot;
462  ret.set(k) = q_r.rem;
463  assert(q==0 || k<dim-1);
464  k++;
465  }
466  k=0;
467  int ip = ret(k);
468  for(;k<dim-1;k++) ip = ret(k+1) + ip * dimensions(k+1);
469  return ip;
470  }
471 
472  template<typename T> int Array<T>::to_first_dim_major_index(int i_last_dim_major) const {
473  int const dim {dimensions.get_ndim()};
474  int k{dim-1};
475  Index ret{dimensions};
476  int q = i_last_dim_major;
477  while(q != 0 && k>=0) {
478  auto q_r = div(q,dimensions(k));
479  q = q_r.quot;
480  ret.set(k) = q_r.rem;
481  assert(q==0 || k>0);
482  k--;
483  }
484  k = dim-1;
485  int ip = ret(k);
486  for(;k>0;k--) ip = ret(k-1) + dimensions(k-1)*ip;
487  return ip;
488  }
489 
490  template <typename T> Array<T>::Array (FILE* fd,Array_ordering order) : dimensions{fd},nbr{0},data{nullptr} {
491  fread_be(&nbr, sizeof(int), 1, fd) ;
492  data.resize(nbr) ;
493  if(order == first_index) fread_be(data.set_data(), sizeof(T), nbr, fd) ;
494  else {
495  Memory_mapped_array<T> lim_data{nbr};
496  fread_be(lim_data.set_data(),sizeof(T),nbr,fd);
497  for(int i=0;i<nbr;i++) data[i] = lim_data[to_last_dim_major_index(i)];
498  }
499  }
500 
501 
502  template <typename T> void Array<T>::save (FILE* fd,Array_ordering order) const {
503  dimensions.save(fd) ;
504  fwrite_be(&nbr, sizeof(int), 1, fd) ;
505  if(order == first_index) fwrite_be(data.get_data(), sizeof(T), nbr, fd) ;
506  else {
507  Memory_mapped_array<T> lim_data {nbr};
508  for(int i=0;i<nbr;i++) lim_data[i] = data[to_first_dim_major_index(i)];
509  fwrite_be(lim_data.get_data(),sizeof(T),nbr,fd);
510  }
511  }
512 
513 // Display
514  template <typename T> ostream& operator<< (ostream& o, const Array<T>& so) {
515  int ndim = so.dimensions.get_ndim() ;
516  o << "Array of " << ndim << " dimension(s)" << endl ;
517  Index xx (so.get_dimensions()) ;
518  switch (ndim) {
519  case 1:
520  for (int i=0 ; i<so.dimensions(0) ; i++) {
521  xx.set(0) = i ;
522  o << so(xx) << " " ;
523  }
524  break ;
525  case 2:
526  for (int i=0 ; i<so.dimensions(1) ; i++) {
527  o << i << " : " ;
528  xx.set(1) = i ;
529  for (int j=0 ; j<so.dimensions(0) ; j++) {
530  xx.set(0) = j ;
531  o << so(xx) << " " ;
532  }
533  o << endl ;
534  }
535  break ;
536  case 3:
537  for (int i=0 ; i<so.dimensions(2) ; i++) {
538  o << " " << i << endl ;
539  xx.set(2) = i ;
540  for (int j=0 ; j<so.dimensions(1) ; j++) {
541  o << j << " : " ;
542  xx.set(1) = j ;
543  for (int k=0 ; k<so.dimensions(0) ; k++) {
544  xx.set(0) = k ;
545  o << so(xx) << " " ;
546  }
547  o << endl ;
548  }
549  o << endl ;
550  }
551  break ;
552  default:
553  for (int i=0 ; i<so.nbr ; i++)
554  o << so.data[i] << " " ;
555  o << endl ;
556  break ;
557  }
558  return o ;
559  }
560 
561 }
562 #endif
Kadath::Array_iterator
Version of the Index class optimized for incremental access to Array components.
Definition: array_iterator.hpp:19
Kadath::Array_iterator::position
int position
Corresponding value for 1D indexing .
Definition: array_iterator.hpp:28
Kadath::Array
Template class for arrays.
Definition: array.hpp:86
Kadath::Array::swap
void swap(Array< T > &so)
Swaps contents between the two arrays (beware when using it with arrays of allocated pointers).
Definition: array.hpp:152
Kadath::Array::pointer
T * pointer
Sylvain's stuff
Definition: array.hpp:100
Kadath::Array::reference
T & reference
Sylvain's stuff.
Definition: array.hpp:96
Kadath::Array::set
reference set(const Index &pos)
Read/write of an element.
Definition: array.hpp:186
Kadath::Array::const_reference
T const & const_reference
Sylvain's stuff.
Definition: array.hpp:98
Kadath::Array::Array
Array(FILE *fd, Array_ordering order=last_index)
Constructor from a file.
Definition: array.hpp:490
Kadath::Array::get_nbr
int get_nbr() const
Returns the total number of elements.
Definition: array.hpp:327
Kadath::Array::delete_data
void delete_data()
Logical destructor (kills the data)
Definition: array.hpp:154
Kadath::Array::operator=
Array & operator=(T xx)
Assigns the same value to all the elements.
Definition: array.hpp:181
Kadath::Array::set
reference set(const Array_iterator &pos)
Read/write of an element.
Definition: array.hpp:210
Kadath::Array::Array
Array(int i, int j)
Constructor for a 2d-array.
Definition: array.hpp:125
Kadath::Array::set_data
pointer set_data()
Direct accessor to the data, read/write version.
Definition: array.hpp:318
Kadath::Array::operator/=
Array< T > & operator/=(const Array< T > &so)
Operator /=.
Definition: array.hpp:352
Kadath::Array::get_ndim
int get_ndim() const
Returns the number of dimensions.
Definition: array.hpp:323
Kadath::Array::set
reference set(int i)
Read/write of an element for a 1d-array.
Definition: array.hpp:215
Kadath::Array::is_increasing
bool is_increasing() const
Checks if a 1D array is increasing.
Definition: array.hpp:444
Kadath::Array::Array
Array(const Dim_array &res)
Constructor from a Dim_array.
Definition: array.hpp:107
Kadath::Array::Array
Array(int i, int j, int k)
Constructor for a 3d-array.
Definition: array.hpp:137
Kadath::Array::get_size
int get_size(int i) const
Returns the size of a given dimension.
Definition: array.hpp:331
Kadath::Array::save
void save(FILE *fd, Array_ordering order=last_index) const
Save in a file.
Definition: array.hpp:502
Kadath::Array::get_dimensions
const Dim_array & get_dimensions() const
Returns the Dim_array of the Array.
Definition: array.hpp:335
Kadath::Array::resize
void resize(Dim_array const &new_dim)
Resize the array by reallocating its ressource.
Definition: array.hpp:160
Kadath::Array::set
reference set(int i, int j)
Read/write of an element for a 2d-array.
Definition: array.hpp:225
Kadath::Array::operator*=
Array< T > & operator*=(const Array< T > &so)
Operator *=.
Definition: array.hpp:349
Kadath::Array::operator()
T operator()(const Index &pos) const
Read only of an element.
Definition: array.hpp:250
Kadath::Array::Array
Array(int i)
Constructor for a 1d-array.
Definition: array.hpp:116
Kadath::Array::dimensions
Dim_array dimensions
Dimensions of the Array.
Definition: array.hpp:88
Kadath::Array::nbr
int nbr
Total number of elements.
Definition: array.hpp:91
Kadath::Array::to_first_dim_major_index
int to_first_dim_major_index(int i_last_dim_major) const
Index formulae to change major dimension indexation from last dimension to first.
Definition: array.hpp:472
Kadath::Array::const_pointer
T const * const_pointer
Sylvain's stuff
Definition: array.hpp:102
Kadath::Array::get_data
const_pointer get_data() const
Direct accessor to the data, read only version.
Definition: array.hpp:313
Kadath::Array::operator-=
Array< T > & operator-=(const Array< T > &so)
Operator -=.
Definition: array.hpp:346
Kadath::Array::operator()
T operator()(Array_iterator const &pos) const
Read only of an element.
Definition: array.hpp:274
Kadath::Array::resize
void resize(int new_size)
Resize overload for the one-dimension array case.
Definition: array.hpp:165
Kadath::Array::to_last_dim_major_index
int to_last_dim_major_index(int i_first_dim_major) const
Index formulae to change major dimension indexation from first dimension to the last.
Definition: array.hpp:454
Kadath::Array::data
Memory_mapped_array< T > data
Elements of the Array.
Definition: array.hpp:92
Kadath::Array::operator+=
Array< T > & operator+=(const Array< T > &so)
Operator +.
Definition: array.hpp:343
Kadath::Array::set
reference set(int i, int j, int k)
Read/write of an element for a 3d-array.
Definition: array.hpp:238
Kadath::Dim_array
Class for storing the dimensions of an array.
Definition: dim_array.hpp:34
Kadath::Dim_array::get_ndim
int get_ndim() const
Returns the number of dimensions.
Definition: dim_array.hpp:63
Kadath::Dim_array::set
int & set(int i)
Read/write of the size of a given dimension.
Definition: dim_array.hpp:54
Kadath::Dim_array::swap
void swap(Dim_array &so)
Sylvain's stuff.
Definition: dim_array.hpp:70
Kadath::Index
Class that gives the position inside a multi-dimensional Array.
Definition: index.hpp:38
Kadath::Index::sizes
Dim_array sizes
Sizes of the associated Array.
Definition: index.hpp:48
Kadath::Index::set
int & set(int i)
Read/write of the position in a given dimension.
Definition: index.hpp:72
Kadath::Matrice
Matrix handling.
Definition: matrice.hpp:38