Hex  1.0
Hydrogen-electron collision solver
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
arrays.h
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
2  * *
3  * / / / / __ \ \ / / *
4  * / /__ / / / _ \ \ \/ / *
5  * / ___ / | |/_/ / /\ \ *
6  * / / / / \_\ / / \ \ *
7  * *
8  * Jakub Benda (c) 2014 *
9  * Charles University in Prague *
10  * *
11 \* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
12 
13 #ifndef HEX_ARRAYS
14 #define HEX_ARRAYS
15 
16 #include <cassert>
17 #include <complex>
18 #include <cstdlib>
19 #include <cstring>
20 #include <fstream>
21 #include <iostream>
22 #include <iomanip>
23 #include <map>
24 #include <numeric>
25 #include <typeinfo>
26 #include <type_traits>
27 #include <vector>
28 
29 #ifndef NO_HDF
30 #include "hdffile.h"
31 #endif
32 
33 #include "complex.h"
34 #include "misc.h"
35 
43 template <class T> class PlainAllocator
44 {
45  public:
46 
54  static T * alloc (size_t n)
55  {
56  return new T[n]();
57  }
58 
68  static void free (T * ptr)
69  {
70  if (ptr != nullptr)
71  delete [] ptr;
72  }
73 };
74 
75 #ifndef NO_ALIGN
76 
89 template <class T, size_t alignment = std::alignment_of<T>::value> class AlignedAllocator
90 {
91  public:
92 
102  static T * alloc (size_t n)
103  {
104  // is there anything to allocate?
105  if (n < 1)
106  return nullptr;
107 
108  // allocate the aligned memory; make sure there will be even number of elements
109  // so that we can always use pairs
110  size_t bytes = (n + (n % 2)) * sizeof(T);
111  size_t align = std::max(alignment, sizeof(void*));
112 
113  // use standard function
114  void* aligned_ptr = nullptr;
115  int err = posix_memalign(&aligned_ptr, align, bytes);
116 
117  // check the return value
118  if (err != 0 or aligned_ptr == nullptr)
119  throw exception ("[AlignedAllocator<T>::alloc] Aligned memory allocation error. Probably out of memory.");
120 
121  // get the number pointer
122  T* ptr = reinterpret_cast<T*>(aligned_ptr);
123 
124  // clear the last element so that we may disregard it during multiplication
125  *(ptr + n + (n % 2) - 1) = 0;
126 
127  // return the pointer
128  return ptr;
129  }
130 
139  static void free (T * ptr)
140  {
141  if (ptr != nullptr)
142  std::free (ptr);
143  }
144 };
145 #endif
146 
147 // Forward declaration of Array (unaligned array of items).
148 template <
149  class T,
150  class Alloc = PlainAllocator<T>
151 > class Array;
152 
153 #ifndef NO_ALIGN
154 // Forward declaration of NumberArray (aligned array of numbers).
155 // - Align at least at multiples of sizeof(void*), but preferably on multiples
156 // of 2*sizeof(Complex). This should enable vectorization of Complex
157 // operations using AVX instruction, bacause two Complex numbers occupy
158 // 2*2*64 = 256 bits, which is the size of AVX register.
159 template <
160  class T,
161  class Alloc = AlignedAllocator <
162  T,
163  larger_of ( sizeof(void*), 2*sizeof(T) )
164  >
165 > class NumberArray;
166 #else
167 // Forward declaration of NumberArray (unaligned array of numbers).
168 template <
169  class T,
170  class Alloc = PlainAllocator<T>
171 > class NumberArray;
172 #endif
173 
186 template <class T> class ArrayView
187 {
188  public:
189 
190  // aliases
191  typedef T DataType;
192  typedef T * iterator;
193  typedef T const * const_iterator;
194 
195  protected:
196 
198  size_t N_;
199 
201  T * array_;
202 
203  public:
204 
205  // constructor
207  : N_(0), array_(nullptr) {}
208 
209  // construct from pointer and size
210  ArrayView (size_t n, T * ptr)
211  : N_(n), array_(ptr) {}
212 
213  // copy constructor
214  ArrayView (ArrayView<T> const & a, size_t i = 0, size_t n = 0)
215  : N_((n > 0) ? n : a.size()), array_(const_cast<T*>(a.data()) + i) { assert(i + N_ <= a.size()); }
216 
217  // construct view from Array const lvalue reference
218  ArrayView (Array<T> const & a, size_t i = 0, size_t n = 0)
219  : N_((n > 0) ? n : a.size()), array_(const_cast<T*>(a.data()) + i) { assert(i + N_ <= a.size()); }
220 
221  // construct view from NumberArray const lvalue reference
222  ArrayView (NumberArray<T> const & a, size_t i = 0, size_t n = 0)
223  : N_((n > 0) ? n : a.size()), array_(const_cast<T*>(a.data()) + i) { assert(i + N_ <= a.size()); }
224 
225  // construct from consecutive memory segment
227  : N_(j - i), array_(const_cast<T*>(&(*i))) {}
228 
229  // construct from right-value reference
231  : ArrayView() { std::swap(N_, r.N_); std::swap(array_, r.array_); }
232 
233  // destructor
234  virtual ~ArrayView () {}
235 
238  {
239  if (v.size() != size())
240  throw exception ("[ArrayView::operator=] Cannot copy %ld elements to %ld fields!", v.size(), N_);
241 
242  for (size_t i = 0; i < size(); i++)
243  array_[i] = v[i];
244 
245  return *this;
246  }
247 
249  T & operator[] (size_t i)
250  {
251 #ifdef NDEBUG
252  return array_[i];
253 #else
254  // bounds check
255  if (i < size())
256  return array_[i];
257  else
258  throw exception ("[ArrayView::operator[]] Index %ld out of bounds (size = %ld) !", i, N_);
259 #endif
260  }
261 
263  T const & operator[] (size_t i) const
264  {
265 #ifdef NDEBUG
266  return array_[i];
267 #else
268  if (i < size())
269  return array_[i];
270  else
271  throw exception ("[ArrayView::operator[]] Index %ld out of bounds (size = %ld) !", i, N_);
272 #endif
273  }
274 
276  size_t size () const { return N_; }
277 
279 
280  virtual T * data () { return array_; }
281  virtual T const * data () const { return array_; }
283 
284  //
285  // STL-like iterator interface
286  //
287 
288  iterator begin () { return data(); }
289  const_iterator begin () const { return data(); }
290  iterator end () { return data() + size(); }
291  const_iterator end () const { return data() + size(); }
292  T & front (int i = 0) { return *(data() + i); }
293  T const & front (int i = 0) const { return *(data() + i); }
294  T & back (int i = 0) { return *(data() + size() - 1 - i); }
295  T const & back (int i = 0) const { return *(data() + size() - 1 - i); }
296 
298  void fill (T x) { for (T & y : *this) y = x; }
299 
301  bool empty () const { return size() == 0; }
302 
304  template <class = typename std::enable_if<is_scalar<T>::value>> double norm () const
305  {
306  double sqrnorm = 0.;
307  for (T const & x : *this)
308  sqrnorm += sqrabs(x);
309  return sqrt(sqrnorm);
310  }
311 };
312 
326 template <class T, class Alloc> class Array : public ArrayView<T>
327 {
328  public:
329 
330  // aliases
331  typedef T DataType;
332  typedef T * iterator;
333  typedef T const * const_iterator;
334 
335  public:
336 
337  // constructors, creates an empty array
338  Array ()
339  : ArrayView<T>() {}
340 
341  // constructor, creates a length-n "x"-filled array
342  Array (size_t n, T x = T(0))
343  : ArrayView<T>(n, Alloc::alloc(n)) { for (size_t i = 0; i < size(); i++) (*this)[i] = x; }
344 
345  // constructor, copies a length-n "array
346  Array (size_t n, T const * const x)
347  : ArrayView<T>(n, Alloc::alloc(n)) { for (size_t i = 0; i < size(); i++) (*this)[i] = x[i]; }
348 
349  // copy constructor from Array const lvalue reference
350  Array (Array<T> const & a)
351  : ArrayView<T>(a.size(), Alloc::alloc(a.size())) { for (size_t i = 0; i < size(); i++) (*this)[i] = a[i]; }
352 
353  // copy constructor from const ArrayView
354  Array (const ArrayView<T> a)
355  : ArrayView<T>(a.size(), Alloc::alloc(a.size())) { for (size_t i = 0; i < size(); i++) (*this)[i] = a[i]; }
356 
357  // copy constructor from Array rvalue reference
359  : ArrayView<T>() { std::swap (ArrayView<T>::N_, a.ArrayView<T>::N_); std::swap (ArrayView<T>::array_, a.ArrayView<T>::array_); }
360 
361  // copy constructor from std::vector
362  Array (std::vector<T> const & a)
363  : ArrayView<T>(a.size(), Alloc::alloc(a.size())) { for (size_t i = 0; i < size(); i++) (*this)[i] = a[i]; }
364 
365  // copy constructor from initializer list
366  Array (std::initializer_list<T> a)
367  : ArrayView<T>(a.end()-a.begin(), Alloc::alloc(a.end()-a.begin())) { size_t i = 0; for (auto it = a.begin(); it != a.end(); it++) (*this)[i++] = *it; }
368 
369  // copy constructor from two forward iterators
370  template <typename ForwardIterator> Array (ForwardIterator i, ForwardIterator j)
371  : ArrayView<T>(std::distance(i,j), Alloc::alloc(std::distance(i,j)))
372  {
373  size_t n = 0;
374  for (ForwardIterator k = i; k != j; k++)
375  (*this)[n++] = *k;
376  }
377 
378  // destructor
380  {
381  if (ArrayView<T>::array_ != nullptr)
382  {
383  Alloc::free (ArrayView<T>::array_);
384  ArrayView<T>::array_ = nullptr;
385  }
386  }
387 
389  size_t size () const { return ArrayView<T>::size(); }
390 
400  virtual size_t resize (size_t n)
401  {
402  if (n == 0)
403  {
404  if (ArrayView<T>::array_ != nullptr)
405  {
406  Alloc::free (ArrayView<T>::array_);
407  ArrayView<T>::array_ = nullptr;
408  ArrayView<T>::N_ = 0;
409  }
410  return 0;
411  }
412 
413  T * new_array = Alloc::alloc(n);
414  for (size_t i = 0; i < n; i++)
415  new_array[i] = (i < size()) ? (*this)[i] : T(0);
416 
417  Alloc::free (ArrayView<T>::array_);
418 
419  ArrayView<T>::N_ = n;
420  ArrayView<T>::array_ = new_array;
421 
422  return size();
423  }
424 
426  inline T & operator[] (size_t i) { return ArrayView<T>::operator[](i); }
427 
429  inline T const & operator[] (size_t i) const { return ArrayView<T>::operator[](i); }
430 
432 
433  virtual T* data () { return ArrayView<T>::data(); }
434  virtual T const * data () const { return ArrayView<T>::data(); }
436 
437  //
438  // STL-like iterator interface
439  //
440 
442  const_iterator begin () const { return ArrayView<T>::begin(); }
443  iterator end () { return ArrayView<T>::end(); }
444  const_iterator end () const { return ArrayView<T>::end(); }
445  T & front (int i = 0) { return ArrayView<T>::front(i); }
446  T const & front (int i = 0) const { return ArrayView<T>::front(i); }
447  T & back (int i = 0) { return ArrayView<T>::back(i); }
448  T const & back (int i = 0) const { return ArrayView<T>::back(i); }
449 
459  virtual void push_back (T const & a)
460  {
461  T* new_array = Alloc::alloc (size() + 1);
462  for (size_t i = 0; i < size(); i++)
463  new_array[i] = std::move(ArrayView<T>::array_[i]);
464  new_array[size()] = a;
466  Alloc::free (ArrayView<T>::array_);
467  ArrayView<T>::array_ = new_array;
468  }
469 
480  virtual T pop_back ()
481  {
482  if (size() > 0)
483  return *(data() + (--ArrayView<T>::N_));
484  else
485  throw exception ("Array has no element to pop!");
486  }
487 
495  template <class InputIterator> void append (
496  InputIterator first, InputIterator last
497  ) {
498  T* new_array = Alloc::alloc(size() + last - first);
499  for (size_t i = 0; i < size(); i++)
500  new_array[i] = (*this)[i];
501  for (InputIterator it = first; it != last; it++)
502  new_array[size() + it - first] = *it;
503  ArrayView<T>::N_ += last - first;
504  Alloc::free (ArrayView<T>::array_);
505  ArrayView<T>::array_ = new_array;
506  }
507 
518  void insert (iterator it, T x)
519  {
520  // create new array (one element longer)
521  T* new_array = Alloc::alloc(size() + 1);
522 
523  // copy everything to the new location
524  for (int i = 0; i < it - begin(); i++)
525  new_array[i] = std::move((*this)[i]);
526 
527  // insert new element
528  *(new_array + (it - begin())) = std::move(x);
529 
530  // copy the rest
531  for (int i = it - begin(); i < (int)size(); i++)
532  new_array[i+1] = std::move((*this)[i]);
533 
534  // change pointers
535  Alloc::free (ArrayView<T>::array_);
537  ArrayView<T>::array_ = new_array;
538  }
539 
541  bool empty () const { return size() == 0; }
542 
543  //
544  // assignment operators
545  //
546 
548  {
549  // if we already have some allocated space, check its size,
550  // so that we do not free it uselessly
551  if (data() != nullptr and size() != b.size())
552  {
553  Alloc::free (ArrayView<T>::array_);
554  ArrayView<T>::array_ = nullptr;
555  }
556 
557  // set the new dimension
558  ArrayView<T>::N_ = b.size();
559 
560  // if necessary, reserve space
561  if (data() == nullptr)
562  ArrayView<T>::array_ = Alloc::alloc(size());
563 
564  // run over the elements
565  for (size_t i = 0; i < size(); i++)
566  (*this)[i] = b[i];
567 
568  return *this;
569  }
570 
572  {
573  // if we already have some allocated space, check its size,
574  // so that we do not free it uselessly
575  if (data() != nullptr)
576  {
577  Alloc::free (ArrayView<T>::array_);
578  ArrayView<T>::array_ = nullptr;
579  ArrayView<T>::N_ = 0;
580  }
581 
582  // swap content
583  std::swap (ArrayView<T>::N_, b.ArrayView<T>::N_);
584  std::swap (ArrayView<T>::array_, b.ArrayView<T>::array_);
585 
586  return *this;
587  }
588 };
589 
604 template <class T, class Alloc> class NumberArray : public Array<T, Alloc>
605 {
606  public:
607 
608  // aliases
609  typedef T DataType;
610  typedef T * iterator;
611  typedef T const * const_iterator;
612 
613  protected:
614 
616  size_t Nres_;
617 
619  std::string name_;
620 
621  public:
622 
623  // default constructor, creates an empty array
625  : Array<T,Alloc>(), Nres_(size()), name_() {}
626 
627  // constructor, creates a length-n "x"-filled array
628  NumberArray (size_t n, T x = T(0))
629  : Array<T,Alloc>(n, x), Nres_(size()), name_() {}
630 
631  // constructor, copies a length-n "array
632  NumberArray (size_t n, T const * x)
633  : Array<T,Alloc>(n, x), Nres_(size()), name_() {}
634 
635  // copy constructor from ArrayView const lvalue reference
637  : Array<T,Alloc>(a), Nres_(size()), name_() {}
638 
639  // copy constructor from Array const lvalue reference
641  : Array<T,Alloc>((ArrayView<T> const &)a), Nres_(size()), name_() {}
642 
643  // copy constructor from std::vector
644  NumberArray (std::vector<T> const & a)
645  : Array<T,Alloc>(a), Nres_(size()), name_() {}
646 
647  // copy constructor from initializer list
648  NumberArray (std::initializer_list<T> a)
649  : Array<T,Alloc>(a), Nres_(size()), name_() {}
650 
651  // copy constructor from two forward iterators
652  template <typename ForwardIterator> NumberArray (ForwardIterator i, ForwardIterator j)
653  : Array<T,Alloc>(i,j), Nres_(size()), name_() {}
654 
655  // copy constructor from Array rvalue reference
657  : Array<T,Alloc>(), name_()
658  {
659  std::swap(ArrayView<T>::N_, a.ArrayView<T>::N_);
660  std::swap(ArrayView<T>::array_, a.ArrayView<T>::array_);
661  std::swap(Nres_,a.Nres_);
662  std::swap(name_, a.name_);
663  }
664 
665  // destructor
667  {
668  // Array<T> will do the destruction
669  // ... do nothing here
670  }
671 
673  size_t size () const { return ArrayView<T>::size(); }
674 
685  virtual size_t resize (size_t n)
686  {
687  if (n <= Nres_)
688  {
689  ArrayView<T>::N_ = n;
690  return n;
691  }
692 
693  Nres_ = n;
694  T * new_array = Alloc::alloc(Nres_);
695 
696  for (size_t i = 0; i < n; i++)
697  new_array[i] = (i < size()) ? (*this)[i] : T(0);
698 
699  Alloc::free (ArrayView<T>::array_);
700  ArrayView<T>::N_ = n;
701  ArrayView<T>::array_ = new_array;
702 
703  return size();
704  }
705 
715  size_t reserve (size_t n)
716  {
717  if (n > Nres_)
718  {
719  Nres_ = n;
720  T* new_array = Alloc::alloc(Nres_);
721 
722  if (size() > 0)
723  {
724  memcpy (new_array, data(), size() * sizeof(T));
725  Alloc::free (ArrayView<T>::array_);
726  }
727 
728  ArrayView<T>::array_ = new_array;
729  }
730 
731  return Nres_;
732  }
733 
735  inline T & operator[] (size_t i) { return ArrayView<T>::operator[](i); }
736 
738  inline T const & operator[] (size_t i) const { return ArrayView<T>::operator[](i); }
739 
757  virtual T * data ()
758  {
759  return ArrayView<T>::array_;
760  }
761  virtual T const * data () const
762  {
763  return ArrayView<T>::array_;
764  }
766 
767  //
768  // STL-like iterator interface
769  //
770 
772  const_iterator begin () const { return ArrayView<T>::begin(); }
774  iterator end () { return ArrayView<T>::end(); }
775  const_iterator end () const { return ArrayView<T>::end(); }
776  const_iterator cend () const { return ArrayView<T>::end(); }
777  T & front (int i = 0) { return ArrayView<T>::front(i); }
778  T const & front (int i = 0) const { return ArrayView<T>::front(i); }
779  T & back (int i = 0) { return ArrayView<T>::back(i); }
780  T const & back (int i = 0) const { return ArrayView<T>::back(i); }
781 
789  virtual void push_front (T const & a)
790  {
791  // can we just whift the whole data?
792  if (size() + 1 <= Nres_)
793  {
794  // shift data by one element
795  for (size_t i = size(); i > 0; i--)
796  (*this)[i] = (*this)[i-1];
797 
798  // set the new front element
799  (*this)[0] = a;
800 
801  // update size of the array
803  }
804  else
805  {
806  // reset storage size
807  Nres_ = 2 * std::max<size_t>(1, Nres_);
808 
809  // reallocate
810  T* new_array = Alloc::alloc(Nres_);
811 
812  // copy original data
813  memcpy(new_array + 1, data(), size() * sizeof(T));
814 
815  // destroy original array
816  if (data() != nullptr)
817  Alloc::free (ArrayView<T>::array_);
818 
819  // use new array
820  ArrayView<T>::array_ = new_array;
821 
822  // set the new front element
823  (*this)[0] = a;
824  }
825  }
826 
833  virtual void push_back (T const & a)
834  {
835  if (size() + 1 > Nres_)
836  {
837  // double the capacity
838  Nres_ = 2 * Nres_ + 1;
839 
840  // allocate space
841  T* new_array = Alloc::alloc(Nres_);
842 
843  // copy original data
844  memcpy(new_array, data(), size() * sizeof(T));
845 
846  // destroy original array
847  if (data() != nullptr)
848  Alloc::free (ArrayView<T>::array_);
849 
850  // use new array
851  ArrayView<T>::array_ = new_array;
852  }
853 
854  // copy new element
855  (*this)[ArrayView<T>::N_++] = a;
856  }
857 
865  {
866  return Array<T,Alloc>::pop_back();
867  }
868 
876  template <class InputIterator> void append
877  (
878  InputIterator first, InputIterator last
879  )
880  {
881  if (size() + last - first > (int)Nres_)
882  {
883  // raise the capacity
884  Nres_ += last - first;
885 
886  // allocate space
887  T* new_array = Alloc::alloc(Nres_);
888 
889  // copy original data
890  memcpy(new_array, data(), size() * sizeof(T));
891 
892  // destroy original array
893  if (data() != nullptr)
894  Alloc::free (ArrayView<T>::array_);
895 
896  // use new array
897  ArrayView<T>::array_ = new_array;
898  }
899 
900  // copy new elements
901  for (InputIterator it = first; it != last; it++)
902  (*this)[ArrayView<T>::N_++] = *it;
903  }
904 
906  bool empty () const
907  {
908  return size() == 0;
909  }
910 
912  void clear ()
913  {
914  memset(ArrayView<T>::array_, 0, size() * sizeof(T));
915  }
916 
918  void drop ()
919  {
920  Nres_ = ArrayView<T>::N_ = 0;
921  Alloc::free (ArrayView<T>::array_);
922  ArrayView<T>::array_ = nullptr;
923  }
924 
925  //
926  // assignment operators
927  //
928 
930  {
931  if (data() == nullptr or Nres_ < b.size())
932  {
933  // delete insufficient storage
934  Alloc::free (ArrayView<T>::array_);
935 
936  // allocate new storage
937  Nres_ = b.size();
938  ArrayView<T>::array_ = Alloc::alloc(Nres_);
939  }
940 
941  // set the new dimension
942  ArrayView<T>::N_ = b.size();
943 
944  // run over the elements
945  for (size_t i = 0; i < size(); i++)
946  (*this)[i] = b[i];
947 
948  return *this;
949  }
950 
952  {
953  std::swap(ArrayView<T>::N_, b.ArrayView<T>::N_);
954  std::swap(ArrayView<T>::array_, b.ArrayView<T>::array_);
955  std::swap(Nres_, b.Nres_);
956  return *this;
957  }
958 
961  {
962  NumberArray<T> c = *this;
963  for (size_t i = 0; i < size(); i++)
964  {
965  Complex z = c[i];
966  c[i] = Complex(z.real(), -z.imag());
967  }
968  return c;
969  }
970 
972  double norm () const
973  {
974  double ret = 0.;
975  for (size_t i = 0; i < size(); i++)
976  {
977  Complex z = (*this)[i];
978  ret += z.real() * z.real() + z.imag() * z.imag();
979  }
980  return sqrt(ret);
981  }
982 
990  template <class Functor> auto transform (Functor f) const -> NumberArray<decltype(f(T(0)))>
991  {
992  size_t n = size();
994 
995  for (size_t i = 0; i < n; i++)
996  c[i] = f((*this)[i]);
997 
998  return c;
999  }
1000 
1002  ArrayView<T> slice (size_t left, size_t right) const
1003  {
1004  assert (right >= left);
1005 
1006  return ArrayView<T>(begin() + left, begin() + right);
1007  }
1008 
1016  std::string toBlob () const
1017  {
1018  // get byte pointer
1019  unsigned char const * dataptr = reinterpret_cast<unsigned char const*>(data());
1020 
1021  // get byte count
1022  size_t count = size() * sizeof(T);
1023 
1024  // resulting string
1025  std::ostringstream hexa;
1026  hexa << "x'" << std::hex << std::setfill('0');
1027 
1028  // for all bytes
1029  for (size_t i = 0; i < count; i++)
1030  hexa << std::setw(2) << static_cast<unsigned>(dataptr[i]);
1031 
1032  hexa << "'";
1033  return hexa.str();
1034  }
1035 
1043  void fromBlob (std::string const & s)
1044  {
1045  if (data() != nullptr and size() != 0)
1046  Alloc::free (ArrayView<T>::array_);
1047 
1048  // the first character outght to be "x" or "X"
1049  // the second character outght to be "'"
1050  // the last character ought to be "'" as well
1051  if ((s[0] != 'x' and s[0] != 'X') or s[1] != '\'' or s.back() != '\'')
1052  throw exception ("[NumberArray::fromBlob] Blob has wrong format, %s.", s.c_str());
1053 
1054  // create substring
1055  std::string ss (s.begin() + 2, s.end() - 1);
1056 
1057  // compute size
1058  size_t bytes = ss.size() / 2;
1059  ArrayView<T>::N_ = bytes / sizeof(T);
1060 
1061  // allocate space
1062  ArrayView<T>::array_ = Alloc::alloc(size());
1063 
1064  // for all bytes
1065  for (size_t i = 0; i < bytes; i++)
1066  {
1067  unsigned byte;
1068  std::stringstream sst;
1069 
1070  // put two hexa digits from "ss" to "sst"
1071  sst << std::hex << ss.substr(2*i, 2);
1072 
1073  // read those two digits as a single byte
1074  sst >> byte;
1075 
1076  // store this byte
1077  reinterpret_cast<char*>(data())[i] = byte;
1078  }
1079  }
1080 
1088  void link (std::string name)
1089  {
1090  name_ = name;
1091  }
1092 
1108  bool hdfsave (std::string name, bool docompress = false, int consec = 10) const
1109  {
1110 #ifndef NO_HDF
1111  // save to HDF file
1112  HDFFile hdf(name, HDFFile::overwrite);
1113 
1114  if (not hdf.valid())
1115  return false;
1116 
1117  if (docompress)
1118  {
1119  NumberArray<int> zero_blocks;
1120  NumberArray<T> elements;
1121  std::tie(zero_blocks,elements) = compress(consec);
1122 
1123  if (not zero_blocks.empty())
1124  {
1125  if (not hdf.write (
1126  "zero_blocks",
1127  &(zero_blocks[0]),
1128  zero_blocks.size()
1129  )) return false;
1130  }
1131  if (not elements.empty())
1132  {
1133  if (not hdf.write (
1134  "array",
1135  &(elements[0]),
1136  elements.size()
1137  )) return false;
1138  }
1139  }
1140  else
1141  {
1142  if (not hdf.write("array", data(), size()))
1143  return false;
1144  }
1145 
1146  return true;
1147 #else
1148  return false;
1149 #endif
1150  }
1151  bool hdfsave () const
1152  {
1153  return hdfsave (name_);
1154  }
1156 
1164  bool hdfload (std::string name)
1165  {
1166 #ifndef NO_HDF
1167  // open the file
1168  HDFFile hdf(name, HDFFile::readonly);
1169 
1170  if (not hdf.valid())
1171  return false;
1172 
1173  NumberArray<T> elements;
1174  NumberArray<int> zero_blocks;
1175 
1176  // read zero blocks
1177  if (zero_blocks.resize(hdf.size("zero_blocks")))
1178  if (not hdf.read("zero_blocks", &(zero_blocks[0]), zero_blocks.size()))
1179  return false;
1180 
1181  // get data size
1182  size_t size = hdf.size("array");
1183 
1184  // scale size if this is a complex array
1185  if (typeid(T) == typeid(Complex))
1186  size /= 2;
1187 
1188  // read packed elements
1189  if (elements.resize(size))
1190  {
1191  if (not hdf.read("array", &(elements[0]), elements.size()))
1192  return false;
1193  }
1194 
1195  // remove previous data
1196  if (data() != nullptr)
1197  {
1198  Alloc::free (ArrayView<T>::array_);
1199  ArrayView<T>::array_ = nullptr;
1200  ArrayView<T>::N_ = Nres_ = 0;
1201  }
1202 
1203  // unpack
1204  *this = elements.decompress(zero_blocks);
1205 
1206  return true;
1207 #else
1208  return false;
1209 #endif
1210  }
1211  bool hdfload ()
1212  {
1213  return hdfload (name_);
1214  }
1216 
1226  std::tuple<NumberArray<int>,NumberArray<T>> compress (int consec) const
1227  {
1228  // compressed array
1229  NumberArray<T> carray;
1230  carray.reserve(size());
1231 
1232  // zero blocks
1233  NumberArray<int> zero_blocks;
1234 
1235  // consecutive zeros counter
1236  int zero_counter = 0;
1237 
1238  // analyze: find compressible segments
1239  for (size_t i = 0; i < size(); i++)
1240  {
1241  if ((*this)[i] == 0.)
1242  {
1243  // another consecutive zero
1244  zero_counter++;
1245  }
1246  else if (zero_counter >= consec)
1247  {
1248  // end of large zero block -> compress
1249  zero_blocks.push_back(i-zero_counter);
1250  zero_blocks.push_back(i);
1251  zero_counter = 0;
1252  }
1253  else
1254  {
1255  // end of tiny zero block -> do not bother with compression
1256  zero_counter = 0;
1257  }
1258 
1259  }
1260  if (zero_counter >= 10)
1261  {
1262  zero_blocks.push_back(size()-zero_counter);
1263  zero_blocks.push_back(size());
1264  }
1265 
1266  // invert selection: get non-zero blocks
1267  NumberArray<int> nonzero_blocks;
1268  nonzero_blocks.push_back(0);
1269  nonzero_blocks.append(zero_blocks.begin(), zero_blocks.end());
1270  nonzero_blocks.push_back(size());
1271 
1272  // compress: copy only nonzero elements
1273  for (size_t iblock = 0; iblock < nonzero_blocks.size()/2; iblock++)
1274  {
1275  int start = nonzero_blocks[2*iblock];
1276  int end = nonzero_blocks[2*iblock+1];
1277  carray.append (
1278  begin() + start,
1279  begin() + end
1280  );
1281  }
1282 
1283  return std::make_tuple(zero_blocks,carray);
1284  }
1285 
1296  NumberArray<T> decompress (NumberArray<int> const & zero_blocks) const
1297  {
1298  if (zero_blocks.empty())
1299  return *this;
1300 
1301  // compute final size
1302  size_t final_size = size();
1303  for (size_t i = 0; i < zero_blocks.size()/2; i++)
1304  final_size += zero_blocks[2*i+1] - zero_blocks[2*i];
1305 
1306  // resize and clean internal storage
1307  NumberArray<DataType> unpack(final_size);
1308  memset(&(unpack[0]), 0, final_size * sizeof(T));
1309 
1310  // copy nonzero chunks
1311  int this_end = 0; // index of last updated element in "this"
1312  int load_end = 0; // index of last used element in "nnz_array"
1313  for (size_t i = 0; i < zero_blocks.size()/2; i++)
1314  {
1315  int zero_start = zero_blocks[2*i];
1316  int zero_end = zero_blocks[2*i+1];
1317 
1318  // append nonzero data before this zero block
1319  memcpy (
1320  &(unpack[0]) + this_end,
1321  begin() + load_end,
1322  (zero_start - this_end) * sizeof(T)
1323  );
1324 
1325  // move cursors
1326  load_end += zero_start - this_end;
1327  this_end = zero_end;
1328  }
1329 
1330  // append remaining data
1331  memcpy (
1332  &(unpack[0]) + this_end,
1333  begin() + load_end,
1334  (final_size - this_end) * sizeof(T)
1335  );
1336 
1337  return unpack;
1338  }
1339 };
1340 
1341 // load array arithmetic operators
1342 #include "arrithm.h"
1343 
1345 template <class T> T operator | (const ArrayView<T> a, const ArrayView<T> b)
1346 {
1347  // get size; check if sizes match
1348  size_t N = a.size();
1349  assert(N == b.size());
1350 
1351  // the scalar product
1352  T result = 0;
1353 
1354  // iterators
1355  T const * const restrict pa = a.data();
1356  T const * const restrict pb = b.data();
1357 
1358  // sum the products
1359  for (size_t i = 0; i < N; i++)
1360  result += pa[i] * pb[i];
1361 
1362  return result;
1363 }
1364 
1373 template <class T1, class T2> auto outer_product (
1374  const ArrayView<T1> a,
1375  const ArrayView<T2> b
1377 {
1378  NumberArray<decltype(T1(0) * T2(0))> c (a.size() * b.size());
1379 
1380  auto ic = c.begin();
1381 
1382  for (auto a_ : a)
1383  for (auto b_ : b)
1384  *(ic++) = a_ * b_;
1385 
1386  return c;
1387 }
1388 
1390 template <typename T> std::ostream & operator << (std::ostream & out, ArrayView<T> const & a)
1391 {
1392  out << "[";
1393  for (size_t i = 0; i < a.size(); i++)
1394  {
1395  if (i == 0)
1396  out << a[i];
1397  else
1398  out << "," << a[i];
1399  }
1400  out << "]";
1401 
1402  return out;
1403 }
1404 
1422 template <typename T> NumberArray<T> linspace (T start, T end, unsigned samples)
1423 {
1424  NumberArray<T> space(samples);
1425 
1426  if (samples == 0)
1427  {
1428  return space;
1429  }
1430 
1431  if (samples == 1)
1432  {
1433  space[0] = start;
1434  return space;
1435  }
1436 
1437  for (unsigned i = 0; i < samples; i++)
1438  {
1439  space[i] = start + ((end - start) * T(i)) / T(samples - 1);
1440  }
1441 
1442  return space;
1443 }
1444 
1462 template <typename T> NumberArray<T> logspace (T x0, T x1, size_t samples)
1463 {
1464  if (x0 <= 0 or x1 <= 0 or x1 < x0)
1465  throw exception ("[logspace] It must be 0 < x1 <= x2 !");
1466 
1467  NumberArray<T> grid(samples);
1468 
1469  if (samples == 0)
1470  {
1471  return grid;
1472  }
1473 
1474  if (samples == 1)
1475  {
1476  grid[0] = x0;
1477  return grid;
1478  }
1479 
1480  for (unsigned i = 0; i < samples; i++)
1481  {
1482  grid[i] = x0 * pow(x1 / x0, i / T(samples - 1));
1483  }
1484 
1485  return grid;
1486 }
1487 
1505 template <typename T> NumberArray<T> geomspace (T x0, T x1, size_t samples, double q)
1506 {
1507  NumberArray<T> grid(samples);
1508 
1509  if (samples == 0)
1510  {
1511  return grid;
1512  }
1513 
1514  if (samples == 1)
1515  {
1516  grid[0] = x0;
1517  return grid;
1518  }
1519 
1520  for (unsigned i = 0; i < samples; i++)
1521  {
1522  grid[i] = x0 + (x1 - x0) * (1. - std::pow(q,i)) / (1. - std::pow(q,samples));
1523  }
1524 
1525  return grid;
1526 }
1527 
1534 template <class T> void write_array
1535 (
1536  const ArrayView<T> array,
1537  const char* filename
1538 );
1539 
1547 template <class T1, class T2> void write_array
1548 (
1549  const ArrayView<T1> grid,
1550  const ArrayView<T2> array,
1551  const char* filename
1552 );
1553 
1560 template <typename Fetcher> bool write_1D_data (size_t m, const char* filename, Fetcher fetch)
1561 {
1562  std::ofstream f(filename);
1563 
1564  if (f.bad())
1565  return false;
1566 
1567  for (size_t i = 0; i < m; i++)
1568  f << fetch(i) << "\n";
1569 
1570  return true;
1571 }
1572 
1586 template <class Fetcher> bool write_2D_data (size_t m, size_t n, const char* filename, Fetcher fetch)
1587 {
1588  std::ofstream f(filename);
1589 
1590  if (f.bad())
1591  return false;
1592 
1593  for (size_t i = 0; i < m; i++)
1594  {
1595  for (size_t j = 0; j < n; j++)
1596  f << fetch(i,j) << " ";
1597  f << "\n";
1598  }
1599 
1600  return true;
1601 }
1602 
1603 //
1604 // aliases
1605 //
1606 
1611 
1616 
1621 
1626 inline rArray concatenate ()
1627 {
1628  return rArray();
1629 }
1630 
1640 template <typename ...Params> rArray concatenate (rArray const & v1, Params ...p)
1641 {
1642  if (sizeof...(p) == 0)
1643  {
1644  return v1;
1645  }
1646  else
1647  {
1648  rArray v2 = concatenate (p...);
1649  rArray v (v1.size() + v2.size());
1650  for (size_t i = 0; i < v1.size(); i++)
1651  v[i] = v1[i];
1652  for (size_t i = 0; i < v2.size(); i++)
1653  v[i + v1.size()] = v2[i];
1654  return v;
1655  }
1656 }
1657 
1658 // return absolute values
1659 rArray abs (const cArrayView u);
1660 rArrays abs (cArrays const &u);
1661 
1663 template <typename T> T min (const ArrayView<T> a)
1664 {
1665  T z = a.front();
1666  for (T const * it = a.begin(); it != a.end(); it++)
1667  if (*it < z)
1668  z = *it;
1669  return z;
1670 }
1671 
1673 template <typename T> T max (const ArrayView<T> a)
1674 {
1675  T z = a.front();
1676  for (T const * it = a.begin(); it != a.end(); it++)
1677  if (*it > z)
1678  z = *it;
1679  return z;
1680 }
1681 
1683 template <typename T> NumberArray<T> pow (NumberArray<T> const & u, double e)
1684 {
1685  NumberArray<T> v(u.size());
1686 
1687  auto iu = u.begin();
1688  auto iv = v.begin();
1689 
1690  while (iu != u.end())
1691  *(iv++) = std::pow(*(iu++), e);
1692 
1693  return v;
1694 }
1695 
1697 template <typename T> Array<T> pow (Array<T> const & u, double e)
1698 {
1699  Array<T> v(u.size());
1700 
1701  auto iu = u.begin();
1702  auto iv = v.begin();
1703 
1704  while (iu != u.end())
1705  *(iv++) = std::pow(*(iu++), e);
1706 
1707  return v;
1708 }
1709 
1711 template <class T> NumberArray<T> sqrt (NumberArray<T> const & A)
1712 {
1713  size_t N = A.size();
1714  NumberArray<T> B (N);
1715 
1716  for (size_t i = 0; i < N; i++)
1717  B[i] = std::sqrt(A[i]);
1718 
1719  return B;
1720 }
1721 
1723 template <class T> NumberArray<T> sin (NumberArray<T> const & A)
1724 {
1725  size_t N = A.size();
1726  NumberArray<T> B (N);
1727 
1728  for (size_t i = 0; i < N; i++)
1729  B[i] = std::sin(A[i]);
1730 
1731  return B;
1732 }
1733 
1735 template <class T> NumberArray<T> asin (NumberArray<T> const & A)
1736 {
1737  size_t N = A.size();
1738  NumberArray<T> B (N);
1739 
1740  for (size_t i = 0; i < N; i++)
1741  B[i] = std::asin(A[i]);
1742 
1743  return B;
1744 }
1745 
1747 template <class T> NumberArray<T> cos (NumberArray<T> const & A)
1748 {
1749  size_t N = A.size();
1750  NumberArray<T> B (N);
1751 
1752  for (size_t i = 0; i < N; i++)
1753  B[i] = std::cos(A[i]);
1754 
1755  return B;
1756 }
1757 
1768 
1770 template <typename T> T sum (const ArrayView<T> v)
1771 {
1772  return std::accumulate(v.begin(), v.end(), T(0));
1773 }
1774 
1776 template <typename T> NumberArray<T> sums (const ArrayView<NumberArray<T>> v)
1777 {
1778  if (v.size() == 0)
1779  return NumberArray<T>(); // empty array
1780 
1781  return std::accumulate (
1782  v.begin(),
1783  v.end(),
1784  NumberArray<T> (v[0].size()),
1786  return a + b;
1787  }
1788  );
1789 }
1790 
1797 template <typename T>
1799 {
1800  Array<bool> v(u.size());
1801  for (size_t i = 0; i < u.size(); i++)
1802  v[i] = (u[i] == x);
1803  return v;
1804 }
1805 
1807 template <typename T>
1809 {
1810  assert(u.size() == v.size());
1811 
1812  Array<bool> w(u.size());
1813  for (size_t i = 0; i < u.size(); i++)
1814  w[i] = (u[i] == v[i]);
1815  return w;
1816 }
1817 
1819 inline bool all (const ArrayView<bool> B)
1820 {
1821  bool ok = true;
1822  for (bool b : B)
1823  ok = ok and b;
1824  return ok;
1825 }
1826 
1828 inline bool any(const ArrayView<bool> B)
1829 {
1830  bool ok = false;
1831  for (bool b : B)
1832  ok = ok or b;
1833  return ok;
1834 }
1835 
1843 template <typename TFunctor, typename TArray>
1844 void eval (TFunctor f, TArray grid, TArray& vals)
1845 {
1846  size_t N = grid.size();
1847  assert(N == vals.size());
1848 
1849  for (size_t i = 0; i < N; i++)
1850  vals[i] = f(grid[i]);
1851 }
1852 
1863 template <typename Tidx, typename Tval> void merge (
1864  NumberArray<Tidx> & idx1, NumberArray<Tval> & arr1,
1865  NumberArray<Tidx> const & idx2, NumberArray<Tval> const & arr2
1866 ){
1867  // positions in arrays
1868  size_t i1 = 0;
1869  size_t i2 = 0;
1870 
1871  // output arrays
1872  NumberArray<Tidx> idx;
1873  NumberArray<Tval> arr;
1874 
1875  // while there is anything to merge
1876  while (i1 < arr1.size() and i2 < arr2.size())
1877  {
1878  if (idx1[i1] == idx2[i2])
1879  {
1880  idx.push_back(idx1[i1]);
1881  arr.push_back(arr1[i1] + arr2[i2]);
1882  i1++;
1883  i2++;
1884  }
1885  else if (idx1[i1] < idx2[i2])
1886  {
1887  idx.push_back(idx1[i1]);
1888  arr.push_back(arr1[i1]);
1889  i1++;
1890  }
1891  else /* idx1[i2] > idx2[i2] */
1892  {
1893  idx.push_back(idx2[i2]);
1894  arr.push_back(arr2[i2]);
1895  i2++;
1896  }
1897  }
1898 
1899  // the rest will be done by a single copy
1900  if (i1 == arr1.size() and i2 < arr2.size())
1901  {
1902  idx.append(idx2.begin() + i2, idx2.end());
1903  arr.append(arr2.begin() + i2, arr2.end());
1904  }
1905 
1906  // copy to the first pair
1907  idx1 = idx;
1908  arr1 = arr;
1909 }
1910 
1912 template <typename T> NumberArray<T> join (const ArrayView<NumberArray<T>> arrays)
1913 {
1914  NumberArray<size_t> partial_sizes(arrays.size() + 1);
1915 
1916  // get partial sizes
1917  partial_sizes[0] = 0;
1918  for (size_t i = 0; i < arrays.size(); i++)
1919  partial_sizes[i+1] = partial_sizes[i] + arrays[i].size();
1920 
1921  // result array
1922  NumberArray<T> res(partial_sizes.back());
1923 
1924  // concatenate arrays
1925  for (size_t i = 0; i < arrays.size(); i++)
1926  {
1927  if (arrays[i].size() > 0)
1928  ArrayView<T>(res, partial_sizes[i], arrays[i].size()) = arrays[i];
1929  }
1930 
1931  return res;
1932 }
1933 
1935 template <class T> NumberArray<T> sorted_unique (const ArrayView<T> v, int n = 1)
1936 {
1937  // create output array
1938  NumberArray<T> w(v.size());
1939 
1940  // iterators
1941  T * iw = w.begin();
1942  T const * iv = v.begin();
1943  int repeat = 0;
1944 
1945  // copy all elements, drop repetitions above "repeat"
1946  while (iv != v.end())
1947  {
1948  // use this element definitely if
1949  // - it is the first element
1950  // - it differs from the preceding element
1951  if (iw == w.begin() or *(iw - 1) != *iv)
1952  repeat = 0;
1953 
1954  // use the conforming element
1955  if (++repeat <= n)
1956  *(iw++) = *iv;
1957 
1958  // move on to the next element
1959  iv++;
1960  }
1961 
1962  // resize and return output array
1963  w.resize(iw - w.begin());
1964 
1965  return w;
1966 }
1967 
1974 template <class T> void smoothen (ArrayView<T> v)
1975 {
1976  if (v.size() == 0)
1977  return;
1978 
1979  T prev = v[0], newv;
1980  for (int i = 1; i < v.size(); i++)
1981  {
1982  newv = 0.5 * (prev + v[i]);
1983  prev = v[i];
1984  v[i] = newv;
1985  }
1986 }
1987 
1989 template <class T> std::string to_string (const ArrayView<T> v, char sep = ' ')
1990 {
1991  std::ostringstream ss;
1992  for (size_t i = 0; i < v.size(); i++)
1993  {
1994  if (i == 0)
1995  ss << v[i];
1996  else
1997  ss << sep << v[i];
1998  }
1999  return ss.str();
2000 }
2001 
2003 rArray threshold (const rArrayView a, double eps);
2004 
2005 #endif
auto outer_product(const ArrayView< T1 > a, const ArrayView< T2 > b) -> NumberArray< decltype(T1(0)*T2(0))>
Outer product of two arrays.
Definition: arrays.h:1373
virtual T * data()
Data pointer.
Definition: arrays.h:757
Array view.
Definition: arrays.h:186
virtual void push_front(T const &a)
Add element to beginning.
Definition: arrays.h:789
NumberArray()
Definition: arrays.h:624
virtual ~ArrayView()
Definition: arrays.h:234
Array(const ArrayView< T > a)
Definition: arrays.h:354
bool hdfsave(std::string name, bool docompress=false, int consec=10) const
Save array to HDF file.
Definition: arrays.h:1108
NumberArray< double > sqrabs(NumberArray< Complex > const &A)
Return per-element square of absolute value.
Definition: arrays.cpp:71
bool empty() const
Check that size equals to zero.
Definition: arrays.h:906
NumberArray(ArrayView< T > const &a)
Definition: arrays.h:636
NumberArray(std::vector< T > const &a)
Definition: arrays.h:644
T & back(int i=0)
Definition: arrays.h:447
Basic memory allocator.
Definition: arrays.h:43
Array(ForwardIterator i, ForwardIterator j)
Definition: arrays.h:370
void fill(T x)
Fill the array with a value.
Definition: arrays.h:298
bool all(const ArrayView< bool > B)
Check that all values are &quot;true&quot;.
Definition: arrays.h:1819
T const * const_iterator
Definition: arrays.h:611
rArray threshold(const rArrayView a, double eps)
Drop small elements of array (replace by zero).
Definition: arrays.cpp:136
T const * const_iterator
Definition: arrays.h:193
T const & front(int i=0) const
Definition: arrays.h:446
const_iterator cbegin() const
Definition: arrays.h:773
T & back(int i=0)
Definition: arrays.h:294
void merge(NumberArray< Tidx > &idx1, NumberArray< Tval > &arr1, NumberArray< Tidx > const &idx2, NumberArray< Tval > const &arr2)
Definition: arrays.h:1863
T * iterator
Definition: arrays.h:610
bool empty() const
Check that size equals to zero.
Definition: arrays.h:541
T & front(int i=0)
Definition: arrays.h:777
T pop_back()
Remove the last element.
Definition: arrays.h:864
const_iterator begin() const
Definition: arrays.h:289
NumberArray< T > sorted_unique(const ArrayView< T > v, int n=1)
Drop all redundant repetitions from sorted array.
Definition: arrays.h:1935
~Array()
Definition: arrays.h:379
std::string to_string(const ArrayView< T > v, char sep= ' ')
Convert ArrayView&lt;T&gt; to a string.
Definition: arrays.h:1989
NumberArray< T > cos(NumberArray< T > const &A)
Return per-element cosine.
Definition: arrays.h:1747
virtual size_t resize(size_t n)
Resize array.
Definition: arrays.h:685
NumberArray< double > atan2(NumberArray< double > const &A, NumberArray< double > const &B)
Return per-element atan2.
Definition: arrays.cpp:58
NumberArray< Complex > cArray
Definition: arrays.h:1610
T & front(int i=0)
Definition: arrays.h:445
A comfortable number array class.
Definition: arrays.h:171
T min(const ArrayView< T > a)
Minimal element of array.
Definition: arrays.h:1663
T DataType
Definition: arrays.h:609
NumberArray< T > & operator=(NumberArray< T > const &b)
Definition: arrays.h:929
virtual void push_back(T const &a)
Append an element to the end.
Definition: arrays.h:459
const_iterator begin() const
Definition: arrays.h:772
NumberArray< double > hypot(NumberArray< double > const &A, NumberArray< double > const &B)
Return per-element hypot.
Definition: arrays.cpp:45
NumberArray(NumberArray< T > const &a)
Definition: arrays.h:640
int size
Definition: misc.h:311
#define restrict
Definition: misc.h:88
class RadialFunction exception
Array< rArray > rArrays
Definition: arrays.h:1619
std::tuple< NumberArray< int >, NumberArray< T > > compress(int consec) const
Get compressed array.
Definition: arrays.h:1226
virtual T const * data() const
Definition: arrays.h:281
NumberArray< double > imagpart(NumberArray< Complex > const &A)
Return per-element imag part.
Definition: arrays.cpp:93
Array< T > & operator=(Array< T > const &b)
Definition: arrays.h:547
const_iterator end() const
Definition: arrays.h:775
constexpr T const & larger_of(T const &a, T const &b)
Constant-expression max.
Definition: misc.h:218
T operator|(const ArrayView< T > a, const ArrayView< T > b)
Scalar product of two arrays.
Definition: arrays.h:1345
NumberArray(size_t n, T const *x)
Definition: arrays.h:632
T const * const_iterator
Definition: arrays.h:333
void smoothen(ArrayView< T > v)
Running average.
Definition: arrays.h:1974
double norm() const
Compute usual 2-norm.
Definition: arrays.h:972
T * array_
Pointer to the array.
Definition: arrays.h:201
bool any(const ArrayView< bool > B)
Check if any value is &quot;true&quot;.
Definition: arrays.h:1828
T & back(int i=0)
Definition: arrays.h:779
NumberArray< T > conj() const
Return complex conjugated array.
Definition: arrays.h:960
iterator begin()
Definition: arrays.h:288
const double e
Definition: special.h:44
Array(std::vector< T > const &a)
Definition: arrays.h:362
T const & front(int i=0) const
Definition: arrays.h:293
iterator end()
Definition: arrays.h:774
void eval(TFunctor f, TArray grid, TArray &vals)
Definition: arrays.h:1844
ArrayView< double > rArrayView
Definition: arrays.h:1614
NumberArray< T > decompress(NumberArray< int > const &zero_blocks) const
Decompress array.
Definition: arrays.h:1296
NumberArray< T > join(const ArrayView< NumberArray< T >> arrays)
Join elements from all subarrays.
Definition: arrays.h:1912
void append(InputIterator first, InputIterator last)
Append more items.
Definition: arrays.h:495
size_t Nres_
Allocated memory.
Definition: arrays.h:616
ArrayView(Array< T > const &a, size_t i=0, size_t n=0)
Definition: arrays.h:218
T & operator[](size_t i)
Element-wise access (non-const).
Definition: arrays.h:735
T * iterator
Definition: arrays.h:192
NumberArray< T > linspace(T start, T end, unsigned samples)
Generate uniform grid.
Definition: arrays.h:1422
ArrayView(size_t n, T *ptr)
Definition: arrays.h:210
bool empty() const
Check whether the size is equal to zero.
Definition: arrays.h:301
T DataType
Definition: arrays.h:191
Array< lArray > lArrays
Definition: arrays.h:1618
double norm() const
Two-norm (defined only for scalar data type).
Definition: arrays.h:304
NumberArray< T > sums(const ArrayView< NumberArray< T >> v)
Sum arrays.
Definition: arrays.h:1776
rArray concatenate()
Definition: arrays.h:1626
Array(Array< T > &&a)
Definition: arrays.h:358
Array(size_t n, T const *const x)
Definition: arrays.h:346
Array< iArray > iArrays
Definition: arrays.h:1617
T & operator[](size_t i)
Element-wise access (non-const).
Definition: arrays.h:249
void link(std::string name)
Link to HDF file.
Definition: arrays.h:1088
NumberArray< double > realpart(NumberArray< Complex > const &A)
Return per-element real part.
Definition: arrays.cpp:82
void insert(iterator it, T x)
Insert item.
Definition: arrays.h:518
ArrayView(ArrayView< T > const &a, size_t i=0, size_t n=0)
Definition: arrays.h:214
const_iterator begin() const
Definition: arrays.h:442
auto transform(Functor f) const -> NumberArray< decltype(f(T(0)))>
Apply a user transformation.
Definition: arrays.h:990
size_t size() const
Length of the array (number of elements).
Definition: arrays.h:276
T & front(int i=0)
Definition: arrays.h:292
T * iterator
Definition: arrays.h:332
const_iterator cend() const
Definition: arrays.h:776
T sum(const ArrayView< T > v)
Sum elements in array.
Definition: arrays.h:1770
std::string name_
HDF file to store to / load from.
Definition: arrays.h:619
iterator begin()
Definition: arrays.h:771
static T * alloc(size_t n)
Allocate array.
Definition: arrays.h:54
ArrayView< Complex > cArrayView
Definition: arrays.h:1615
T const & front(int i=0) const
Definition: arrays.h:778
iterator end()
Definition: arrays.h:443
NumberArray< long > lArray
Definition: arrays.h:1608
NumberArray< T > sin(NumberArray< T > const &A)
Return per-element sine.
Definition: arrays.h:1723
NumberArray(NumberArray< T > &&a)
Definition: arrays.h:656
size_t size() const
Return number of elements.
Definition: arrays.h:389
size_t N_
Number of elements in the array.
Definition: arrays.h:198
ArrayView(ArrayView< T > &&r)
Definition: arrays.h:230
size_t reserve(size_t n)
Reserve memory.
Definition: arrays.h:715
ArrayView< int > iArrayView
Definition: arrays.h:1612
bool hdfload()
Definition: arrays.h:1211
ArrayView(NumberArray< T > const &a, size_t i=0, size_t n=0)
Definition: arrays.h:222
bool hdfsave() const
Definition: arrays.h:1151
NumberArray< T > pow(NumberArray< T > const &u, double e)
Return per-element power.
Definition: arrays.h:1683
void append(InputIterator first, InputIterator last)
Append a range of values at end.
Definition: arrays.h:877
A comfortable data array class.
Definition: arrays.h:151
virtual T * data()
Data pointer.
Definition: arrays.h:433
static void free(T *ptr)
Free array.
Definition: arrays.h:68
ArrayView(const_iterator i, const_iterator j)
Definition: arrays.h:226
NumberArray< T > asin(NumberArray< T > const &A)
Return per-element arcsine.
Definition: arrays.h:1735
virtual size_t resize(size_t n)
Resize array.
Definition: arrays.h:400
NumberArray< double > rArray
Definition: arrays.h:1609
Array(std::initializer_list< T > a)
Definition: arrays.h:366
ArrayView()
Definition: arrays.h:206
name_()
Definition: arrays.h:629
virtual T * data()
Pointer to the data.
Definition: arrays.h:280
T DataType
Definition: arrays.h:331
virtual T pop_back()
Remove the last element from the array.
Definition: arrays.h:480
const_iterator end() const
Definition: arrays.h:291
Exception class.
Definition: misc.h:39
T const & back(int i=0) const
Definition: arrays.h:295
T const & back(int i=0) const
Definition: arrays.h:448
NumberArray< T > geomspace(T x0, T x1, size_t samples, double q)
Generate geometric grid.
Definition: arrays.h:1505
void fromBlob(std::string const &s)
Convert from SQL BLOB.
Definition: arrays.h:1043
ArrayView< long > lArrayView
Definition: arrays.h:1613
virtual T const * data() const
Definition: arrays.h:434
bool hdfload(std::string name)
Load array from HDF file.
Definition: arrays.h:1164
T const & back(int i=0) const
Definition: arrays.h:780
NumberArray(std::initializer_list< T > a)
Definition: arrays.h:648
T & operator[](size_t i)
Element-wise access (non-const).
Definition: arrays.h:426
ArrayView< T > slice(size_t left, size_t right) const
Return a subarray using ArrayView.
Definition: arrays.h:1002
NumberArray< T > sqrt(NumberArray< T > const &A)
Return per-element square root.
Definition: arrays.h:1711
virtual void push_back(T const &a)
Add element to end.
Definition: arrays.h:833
Array(Array< T > const &a)
Definition: arrays.h:350
rArray abs(const cArrayView u)
Definition: arrays.cpp:19
std::complex< double > Complex
Definition: complex.h:20
iterator end()
Definition: arrays.h:290
NumberArray(ForwardIterator i, ForwardIterator j)
Definition: arrays.h:652
bool write_2D_data(size_t m, size_t n, const char *filename, Fetcher fetch)
Definition: arrays.h:1586
Array()
Definition: arrays.h:338
T max(const ArrayView< T > a)
Maximal element of array.
Definition: arrays.h:1673
~NumberArray()
Definition: arrays.h:666
ArrayView< T > & operator=(const ArrayView< T > v)
Assignment operator.
Definition: arrays.h:237
virtual T const * data() const
Definition: arrays.h:761
std::string toBlob() const
Convert to SQL BLOB.
Definition: arrays.h:1016
const_iterator end() const
Definition: arrays.h:444
void write_array(const ArrayView< double > array, const char *filename)
Definition: arrays.cpp:104
bool write_1D_data(size_t m, const char *filename, Fetcher fetch)
Write elements.
Definition: arrays.h:1560
void drop()
Reset array: deallocate everything, resize to zero.
Definition: arrays.h:918
Array< cArray > cArrays
Definition: arrays.h:1620
Array< bool > operator==(const ArrayView< T > u, T x)
Definition: arrays.h:1798
size_t size() const
Item count.
Definition: arrays.h:673
iterator begin()
Definition: arrays.h:441
void clear()
Fill array with zeros.
Definition: arrays.h:912
NumberArray< T > logspace(T x0, T x1, size_t samples)
Generate logarithmic grid.
Definition: arrays.h:1462
NumberArray< int > iArray
Definition: arrays.h:1607