ODIN
odindata/data.h
1 /***************************************************************************
2  data.h - description
3  -------------------
4  begin : Fri Apr 6 2001
5  copyright : (C) 2000-2021 by Thies Jochimsen & Michael von Mengershausen
6  email : thies@jochimsen.de mengers@cns.mpg.de
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #ifndef DATA_H
19 #define DATA_H
20 
21 #ifndef TJUTILS_CONFIG_H
22 #define TJUTILS_CONFIG_H
23 #include <tjutils/config.h>
24 #endif
25 
26 
27 #include <tjutils/tjstd.h>
28 #include <tjutils/tjcomplex.h>
29 #include <tjutils/tjfeedback.h>
30 
31 
32 #ifdef ODIN_DEBUG
33 #define BZ_DEBUG
34 #endif
35 
36 
37 #ifdef BLITZ_REPLACEMENT
38 
39 #include <../replacements/blitz_replacement.h>
40 
41 
42 #else
43 
44 
45 #include <blitz/blitz.h> // include 1st to overwrite mutex defines
46 
47 // undef mutex defines
48 #ifdef BZ_MUTEX_DECLARE
49 #undef BZ_MUTEX_DECLARE
50 #endif
51 #ifdef BZ_MUTEX_INIT
52 #undef BZ_MUTEX_INIT
53 #endif
54 #ifdef BZ_MUTEX_LOCK
55 #undef BZ_MUTEX_LOCK
56 #endif
57 #ifdef BZ_MUTEX_UNLOCK
58 #undef BZ_MUTEX_UNLOCK
59 #endif
60 #ifdef BZ_MUTEX_DESTROY
61 #undef BZ_MUTEX_DESTROY
62 #endif
63 
64 // do not use mutex in Blitz, even if pthread is available
65 // to get consistent behaviour across all platforms
66 #define BZ_MUTEX_DECLARE(name)
67 #define BZ_MUTEX_INIT(name)
68 #define BZ_MUTEX_LOCK(name)
69 #define BZ_MUTEX_UNLOCK(name)
70 #define BZ_MUTEX_DESTROY(name)
71 
72 
73 
74 #include <blitz/array.h>
75 #ifdef HAVE_BLITZ_TINYVEC_ET
76 #include <blitz/tinyvec-et.h> // to use arithmetics of TinyVectors
77 #endif
78 using namespace blitz;
79 
80 #endif // BLITZ_REPLACEMENT
81 
82 
83 
84 #include <tjutils/tjtypes.h>
85 #include <tjutils/tjprofiler.h>
86 
87 #include <odinpara/protocol.h>
88 
89 #include <odindata/converter.h>
90 #include <odindata/fileio_opts.h>
91 
92 
93 
94 #ifndef BLITZ_REPLACEMENT
95 
96 // some additional macros for complex numbers
97 BZ_DECLARE_FUNCTION_RET (float2real, STD_complex);
98 BZ_DECLARE_FUNCTION_RET (float2imag, STD_complex);
99 BZ_DECLARE_FUNCTION_RET (creal, float)
100 BZ_DECLARE_FUNCTION_RET (cimag, float)
101 BZ_DECLARE_FUNCTION_RET (cabs, float)
102 BZ_DECLARE_FUNCTION_RET (phase, float)
103 BZ_DECLARE_FUNCTION (logc)
104 BZ_DECLARE_FUNCTION (expc)
105 BZ_DECLARE_FUNCTION (conjc)
106 
107 // Other useful macros
108 BZ_DECLARE_FUNCTION2(secureDivision)
109 BZ_DECLARE_FUNCTION (secureInv)
110 
111 
112 #ifdef STL_REPLACEMENT
113 BZ_DECLARE_ARRAY_ET_SCALAR_OPS(tjstd_complex)
114 #endif
115 
116 #endif // BLITZ_REPLACEMENT
117 
118 
119 
126 
127 struct FileMapHandle {
128  FileMapHandle() : fd(-1), offset(0), refcount(1) {}
129  int fd;
130  LONGEST_INT offset; // in bytes
131  int refcount;
132  Mutex mutex; // protect refcount
133 };
134 
136 
147 template <typename T, int N_rank>
148 class Data : public Array<T,N_rank> {
149 
150 public:
151 
152 
156  Data(const TinyVector<int,N_rank>& dimvec, const T& val=0) : Array<T,N_rank>(dimvec), fmap(0) {(*this)=val;}
157 
158 
159  // Variable extent constructors
160  Data() : fmap(0) {}
161  Data(int extent1) : Array<T,N_rank>(extent1), fmap(0) {}
162  Data(int extent1, int extent2) : Array<T,N_rank>(extent1,extent2), fmap(0) {}
163  Data(int extent1, int extent2,int extent3) : Array<T,N_rank>(extent1,extent2,extent3), fmap(0) {}
164  Data(int extent1, int extent2,int extent3,int extent4) : Array<T,N_rank>(extent1,extent2,extent3,extent4), fmap(0) {}
165  Data(int extent1, int extent2,int extent3,int extent4,int extent5) : Array<T,N_rank>(extent1,extent2,extent3,extent4,extent5), fmap(0) {}
166 // Data(int extent1, int extent2,int extent3,int extent4,int extent5,int extent6) : Array<T,N_rank>(extent1,extent2,extent3,extent4,extent5,extent6), fmap(0) {}
167 // Data(int extent1, int extent2,int extent3,int extent4,int extent5,int extent6,int extent7) : Array<T,N_rank>(extent1,extent2,extent3,extent4,extent5,extent6,extent7), fmap(0) {}
168 // Data(int extent1, int extent2,int extent3,int extent4,int extent5,int extent6,int extent7,int extent8) : Array<T,N_rank>(extent1,extent2,extent3,extent4,extent5,extent6,extent7,extent8), fmap(0) {}
169 // Data(int extent1, int extent2,int extent3,int extent4,int extent5,int extent6,int extent7,int extent8,int extent9) : Array<T,N_rank>(extent1,extent2,extent3,extent4,extent5,extent6,extent7,extent8,extent9), fmap(0) {}
170 
171 
180  Data(const STD_string& filename, bool readonly, const TinyVector<int,N_rank>& shape, LONGEST_INT offset=0 );
181 
185  Data(const Data<T,N_rank>& d) : fmap(0) {Data<T,N_rank>::reference(d);}
186 
190  Data(const Array<T,N_rank>& a) : Array<T,N_rank>(a), fmap(0) {}
191 
195  Data(const tjarray<tjvector<T>,T>& a) : fmap(0) { (*this)=a;}
196 
197 
198 #ifdef BLITZ_REPLACEMENT
202  template<int M_rank>
203  Data(const Array<T,M_rank>& a) : Array<T,N_rank>(a), fmap(0) {}
204 #else
208  template<class T_expr>
209  Data(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr) : Array<T,N_rank>(expr), fmap(0) {}
210 #endif
211 
212 
216  ~Data();
217 
218 
222  Data<T,N_rank>& operator = (const Data<T,N_rank>& d) {Array<T,N_rank>::operator=(d); return *this;}
223 
227  Data<T,N_rank>& operator = (const Array<T,N_rank>& a) {Array<T,N_rank>::operator=(a); return *this;}
228 
232  Data<T,N_rank>& operator = (const T& val) {Array<T,N_rank>::operator=(val); return *this;}
233 
234 
235 #ifdef BLITZ_REPLACEMENT
239  template<int M_rank>
240  inline Data<T,N_rank>& operator = (const RectView<const Array<T,M_rank>&,T,M_rank>& rhs) {
241  Array<T,N_rank>::operator = (rhs);
242  return *this;
243  }
244 #else
248  template<class T_expr>
249  inline Data<T,N_rank>& operator = (BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr) {
250  Data<T,N_rank>::reference(Data<T,N_rank>(expr)); // forwarding to constructor
251  return *this;
252  }
253 #endif
254 
255 
259  Data<T,N_rank>& operator = (const tjarray<tjvector<T>,T>& a);
260 
261 
262  // resolve ambiguities, appears to be obsolete with blitz-0.9 and later versions
263 // Array<T,N_rank> operator + (const Data<T,N_rank>& b) const {return Array<T,N_rank>(Array<T,N_rank>(*this)+Array<T,N_rank>(b));}
264 // Array<T,N_rank> operator - (const Data<T,N_rank>& b) const {return Array<T,N_rank>(Array<T,N_rank>(*this)-Array<T,N_rank>(b));}
265 // Array<T,N_rank> operator * (const Data<T,N_rank>& b) const {return Array<T,N_rank>(Array<T,N_rank>(*this)*Array<T,N_rank>(b));}
266 // Array<T,N_rank> operator / (const Data<T,N_rank>& b) const {return Array<T,N_rank>(Array<T,N_rank>(*this)/Array<T,N_rank>(b));}
267 
268 
269  // resolve ambiguities, appears to be obsolete with blitz-0.9 and later versions
270 // Array<T,N_rank> operator + (const T& v) const {return Array<T,N_rank>(Array<T,N_rank>(*this)+v);}
271 // Array<T,N_rank> operator - (const T& v) const {return Array<T,N_rank>(Array<T,N_rank>(*this)-v);}
272 // Array<T,N_rank> operator * (const T& v) const {return Array<T,N_rank>(Array<T,N_rank>(*this)*v);}
273 // Array<T,N_rank> operator / (const T& v) const {return Array<T,N_rank>(Array<T,N_rank>(*this)/v);}
274 
275 
276 
283  int read(const STD_string& format,const STD_string& filename, LONGEST_INT offset=0);
284 
285 
291  template<typename T2>
292  int read(const STD_string& filename, LONGEST_INT offset=0);
293 
294 
301  int write(const STD_string& format,const STD_string& filename, bool autoscale=true) const;
302 
303 
309  template<typename T2>
310  int write(const STD_string& filename, bool autoscale=true) const;
311 
312 
313 
319  int write(const STD_string& filename, fopenMode mode=overwriteMode) const;
320 
330  int write_asc_file(const STD_string& filename, const Array<T,N_rank>& pre=defaultArray, const Array<T,N_rank>& post=defaultArray) const;
331 
336  int read_asc_file(const STD_string& filename);
337 
338 
344  int autowrite(const STD_string& filename, const FileWriteOpts& opts=FileWriteOpts(), const Protocol* prot=0) const;
345 
352  int autoread(const STD_string& filename, const FileReadOpts& opts=FileReadOpts(), Protocol* prot=0, ProgressMeter* progmeter=0);
353 
357  bool is_filemapped() const {return fmap;}
358 
362  operator tjarray<tjvector<T>,T> () const;
363 
364 
369  TinyVector<int, N_rank> create_index(unsigned long index) const;
370 
374  unsigned long create_linear_index(const TinyVector<int, N_rank>& indexvec) const;
375 
381  void shift(unsigned int shift_dim, int shift);
382 
392  void congrid(const TinyVector<int,N_rank>& newshape, const TinyVector<float,N_rank>* subpixel_shift=0, bool left_to_right=false);
393 
394 
401  template<typename T2, int N_rank2> Data<T2,N_rank2>& convert_to(Data<T2,N_rank2>& dst, bool autoscale=true) const;
402 
403  // specialization which uses reference if type/rank are same and scale/offset==1.0/0.0 to reduce memory usage
404  Data<T,N_rank>& convert_to(Data<T,N_rank>& dst, bool autoscale=true) const;
405 
406 
410  T* c_array();
411 
412 
413  // Reimplemented Blitz::Array functions
414  void reference(const Data<T,N_rank>& d);
415 
416 
417  // dummy objects used for default arguments
418  static Array<T,N_rank> defaultArray;
419 
420 
421 
422  protected:
423  void detach_fmap(); // to be accessed by ComplexData
424 
425 
426  private:
427 
428  void interpolate1dim(unsigned int dim, int newsize, float subpixel_shift);
429 
430 
431  FileMapHandle* fmap;
432 
433 };
434 
435 template <typename T, int N_rank>
436 Array<T,N_rank> Data<T,N_rank>::defaultArray;
437 
438 
443 
444 #ifdef STREAM_REPLACEMENT
445 
446 template<typename T, int N_rank>
447 inline STD_ostream& operator << (STD_ostream& s,const Array<T,N_rank>&) {return s;}
448 
449 template<typename T, int N_rank>
450 inline STD_ostream& operator << (STD_ostream& s,const TinyVector<T,N_rank>&) { return s;}
451 
452 inline STD_ostream& operator << (STD_ostream& s,const Range&) { return s;}
453 
454 #endif
455 
457 
458 
459 /*
460  * Create n-dim index from linear 'index' for a given 'shape'
461  */
462 template<int N_rank>
463 TinyVector<int,N_rank> index2extent(const TinyVector<int,N_rank>& shape, unsigned int index) {
464  TinyVector<int,N_rank> result;
465  unsigned int subindex=index;
466  for(int i=N_rank-1;i>=0;i--){
467  result(i)=subindex%shape(i);
468  subindex/=shape(i);
469  }
470  return result;
471 }
472 
473 
475 
476 
477 // leave this function outside of class Data, otherwise it screws up
478 // when using optimization
479 template<typename T, int N_rank, typename T2>
480 void convert_from_ptr(Data<T,N_rank>& dst, const T2* src, const TinyVector<int, N_rank>& shape, bool autoscale=true) {
481  Log<OdinData> odinlog("Data","convert_from_ptr");
482 
483  int dstsize=product(shape);
484  int srcsize=dstsize*Converter::get_elements((T)0)/Converter::get_elements((T2)0);
485  ODINLOG(odinlog,normalDebug) << "dstsize/srcsize=" << dstsize << "/" << srcsize << STD_endl;
486 
487  dst.resize(shape);
488 
489  Converter::convert_array(src, dst.c_array(), srcsize, dstsize, autoscale);
490 }
491 
492 // specialization in case the type is the same
493 template<typename T, int N_rank>
494 void convert_from_ptr(Data<T,N_rank>& dst, const T* src, const TinyVector<int, N_rank>& shape) {
495  dst.reference(Array<T,N_rank>((T*)src, shape, duplicateData));
496 }
497 
498 
500 
501 template <typename T, int N_rank>
502 Data<T,N_rank>::Data(const STD_string& filename, bool readonly, const TinyVector<int,N_rank>& shape, LONGEST_INT offset ) : fmap(new FileMapHandle) {
503  T* ptr=(T*)filemap(filename, (LONGEST_INT)product(shape)*sizeof(T), offset, readonly, fmap->fd);
504  if(ptr && (fmap->fd)>=0) { // only negative file descriptor indicates filemap failure
505  Array<T,N_rank>::reference(Array<T,N_rank>(ptr, shape, neverDeleteData));
506  fmap->offset=offset;
507  } else {
508  delete fmap;
509  fmap=0;
510  }
511 }
512 
514 
515 template <typename T, int N_rank>
516 void Data<T,N_rank>::detach_fmap() {
517  Log<OdinData> odinlog("Data","detach_fmap");
518  if(fmap) {
519  fmap->mutex.lock();
520  (fmap->refcount)--;
521  ODINLOG(odinlog,normalDebug) << "fd/offset/refcount=" << fmap->fd << "/" << fmap->offset << "/" << fmap->refcount << STD_endl;
522  if(!(fmap->refcount)) {
523  fileunmap(fmap->fd, Array<T,N_rank>::data(), (LONGEST_INT)Array<T,N_rank>::size()*sizeof(T), fmap->offset);
524  fmap->mutex.unlock();
525  delete fmap;
526  fmap=0;
527  }
528  if(fmap) fmap->mutex.unlock();
529  }
530 }
531 
533 
534 template <typename T, int N_rank>
535 Data<T,N_rank>::~Data() {
536  detach_fmap();
537 }
538 
540 
541 template <typename T, int N_rank>
542 void Data<T,N_rank>::reference(const Data<T,N_rank>& d) {
543  Log<OdinData> odinlog("Data","reference");
544  detach_fmap();
545  fmap=d.fmap;
546  if(fmap) {
547  MutexLock lock(fmap->mutex);
548  (fmap->refcount)++;
549  ODINLOG(odinlog,normalDebug) << "fmap->refcount=" << fmap->refcount << STD_endl;
550  }
551  Array<T,N_rank>::reference(d);
552 }
553 
555 
556 template <typename T, int N_rank>
557 Data<T,N_rank>& Data<T,N_rank>::operator = (const tjarray<tjvector<T>,T>& a) {
558  Log<OdinData> odinlog("Data","=");
559  if( a.dim() <= N_rank ) {
560  ndim nn=a.get_extent();
561  int npad=N_rank-nn.dim();
562  ODINLOG(odinlog,normalDebug) << "npad=" << npad << STD_endl;
563  for(int ipad=0; ipad<npad; ipad++) nn.add_dim(1, true); // pad at front
564  ODINLOG(odinlog,normalDebug) << "nn(pad)" << nn << STD_endl;
565  TinyVector<int,N_rank> tv;
566  for(unsigned int i=0; i<N_rank; i++) tv(i)=nn[i];
567  Data<T,N_rank>::resize(tv); // scope reuired by GCC 4.7
568  for(unsigned int i=0;i<a.total();i++) (*this)(create_index(i))=a[i];
569  } else ODINLOG(odinlog,errorLog) << "dimension mismatch: this=" << N_rank<< " < tjarray=" << a.dim() << STD_endl;
570  return *this;
571 }
572 
573 
575 
576 template<typename T,int N_rank>
577 int Data<T,N_rank>::read(const STD_string& format,const STD_string& filename, LONGEST_INT offset) {
578  Log<OdinData> odinlog("Data","read");
579  if(format==TypeTraits::type2label((u8bit)0)) return read<u8bit> (filename);
580  if(format==TypeTraits::type2label((s8bit)0)) return read<s8bit> (filename);
581  if(format==TypeTraits::type2label((u16bit)0)) return read<u16bit>(filename);
582  if(format==TypeTraits::type2label((s16bit)0)) return read<s16bit>(filename);
583  if(format==TypeTraits::type2label((u32bit)0)) return read<u32bit>(filename);
584  if(format==TypeTraits::type2label((s32bit)0)) return read<s32bit>(filename);
585  if(format==TypeTraits::type2label((float)0)) return read<float> (filename);
586  if(format==TypeTraits::type2label((double)0)) return read<double>(filename);
587  ODINLOG(odinlog,errorLog) << "Unable to read file " << filename << " with data type " << format << STD_endl;
588  return -1;
589 }
590 
591 
593 
594 template<typename T,int N_rank>
595 template<typename T2>
596 int Data<T,N_rank>::read(const STD_string& filename, LONGEST_INT offset) {
597  Log<OdinData> odinlog("Data","read");
598  LONGEST_INT fsize=filesize(filename.c_str())-offset;
599  LONGEST_INT nelements_file=fsize/sizeof(T2);
600  LONGEST_INT length=product(Array<T,N_rank>::shape());
601  ODINLOG(odinlog,normalDebug) << "fsize / nelements_file / length = " << fsize << " / " << nelements_file << " / " << length << STD_endl;
602 
603  if(!length) return 0;
604 
605  if(nelements_file<length) {
606  ODINLOG(odinlog,errorLog) << "Size of file " << filename << " to small for reading" << STD_endl;
607  return -1;
608  }
609 
610  STD_string srctype=TypeTraits::type2label((T2)0);
611  STD_string dsttype=TypeTraits::type2label((T)0);
612  ODINLOG(odinlog,normalDebug) << "srctype/dsttype=" << srctype << "/" << dsttype << STD_endl;
613 
614  TinyVector<int,N_rank> fileshape(Array<T,N_rank>::shape());
615  fileshape(N_rank-1) *= (Converter::get_elements((T)0)/Converter::get_elements((T2)0)); // Adjust extent of last dim for copmlex data
616  Data<T2,N_rank> filedata(filename, true, fileshape, offset);
617  filedata.convert_to(*this);
618 
619  return 0;
620 }
621 
622 
624 
625 template<typename T,int N_rank>
626 int Data<T,N_rank>::write(const STD_string& format, const STD_string& filename, bool autoscale) const {
627  Log<OdinData> odinlog("Data","write");
628  if(format==TypeTraits::type2label((u8bit)0)) return write<u8bit> (filename,autoscale);
629  if(format==TypeTraits::type2label((s8bit)0)) return write<s8bit> (filename,autoscale);
630  if(format==TypeTraits::type2label((u16bit)0)) return write<u16bit>(filename,autoscale);
631  if(format==TypeTraits::type2label((s16bit)0)) return write<s16bit>(filename,autoscale);
632  if(format==TypeTraits::type2label((u32bit)0)) return write<u32bit>(filename,autoscale);
633  if(format==TypeTraits::type2label((s32bit)0)) return write<s32bit>(filename,autoscale);
634  if(format==TypeTraits::type2label((float)0)) return write<float> (filename,autoscale);
635  if(format==TypeTraits::type2label((double)0)) return write<double>(filename,autoscale);
636  ODINLOG(odinlog,errorLog) << "Unable to write file " << filename << " with data type " << format << STD_endl;
637  return -1;
638 }
639 
641 
642 template<typename T,int N_rank>
643 int Data<T,N_rank>::write(const STD_string& filename, fopenMode mode) const {
644  Log<OdinData> odinlog("Data","write");
645  if(filename=="") return 0;
646 
647  FILE* file_ptr=ODIN_FOPEN(filename.c_str(),modestring(mode));
648 
649  if(file_ptr==NULL) {
650  ODINLOG(odinlog,errorLog) << "unable to create/open file >" << filename << "< - " << lasterr() << STD_endl;
651  return -1;
652  }
653 
654  Data<T,N_rank> data_copy(*this); // for contig memory
655 
656  LONGEST_INT nmemb=Array<T,N_rank>::numElements();
657  LONGEST_INT count=fwrite(data_copy.c_array(),sizeof(T),nmemb,file_ptr);
658  if(count!=nmemb) {
659  ODINLOG(odinlog,errorLog) << "unable to fwrite to file >" << filename << "< - " << lasterr() << STD_endl;
660  return -1;
661  }
662 
663  if(file_ptr!=NULL) fclose(file_ptr);
664  return 0;
665 }
666 
668 
669 template<typename T,int N_rank>
670 template<typename T2>
671  int Data<T,N_rank>::write(const STD_string& filename, bool autoscale) const {
672  Log<OdinData> odinlog("Data","write");
673 
674  rmfile(filename.c_str()); // remove old file to get new file size with mmap
675 
676  ODINLOG(odinlog,normalDebug) << "mem(pre): " << Profiler::get_memory_usage() << STD_endl;
677  Data<T2,N_rank> converted_data; convert_to(converted_data,autoscale);
678  ODINLOG(odinlog,normalDebug) << "mem(cnv): " << Profiler::get_memory_usage() << STD_endl;
679  Data<T2,N_rank> filedata(filename,false,converted_data.shape());
680  ODINLOG(odinlog,normalDebug) << "mem(map): " << Profiler::get_memory_usage() << STD_endl;
681  filedata=converted_data;
682 
683  return 0;
684 }
685 
686 
688 
689 template<typename T,int N_rank>
690 int Data<T,N_rank>::write_asc_file(const STD_string& filename, const Array<T,N_rank>& pre, const Array<T,N_rank>& post) const {
691  bool have_pre=false;
692  bool have_post=false;
693 
694  Data<T,N_rank> pre_data(pre);
695  Data<T,N_rank> post_data(post);
696  unsigned int n=Array<T,N_rank>::numElements();
697 
698  if(pre_data.numElements()==n) have_pre=true;
699  if(post_data.numElements()==n) have_post=true;
700 
701  STD_ofstream ofs(filename.c_str());
702  if(ofs.bad()) return -1;
703 
704  T val;
705 
706  for(unsigned int i=0; i<n; i++) {
707  if(have_pre) {
708  val=pre_data(pre_data.create_index(i));
709  ofs << val << " ";
710  }
711  val=(*this)(create_index(i));
712  ofs << val;
713  if(have_post) {
714  val=post_data(post_data.create_index(i));
715  ofs << " " << val;
716  }
717  ofs << "\n";
718  }
719 
720  ofs.close();
721  return 0;
722 }
723 
725 
726 template<typename T,int N_rank>
727 int Data<T,N_rank>::read_asc_file(const STD_string& filename) {
728  STD_ifstream ifs(filename.c_str());
729  if(ifs.bad()) return -1;
730 
731  STD_string valstr;
732  for(unsigned int i=0; i<Array<T,N_rank>::numElements(); i++) {
733  if(ifs.bad()) return -1;
734  ifs >> valstr;
735  TypeTraits::string2type(valstr,(*this)(create_index(i)));
736  }
737 
738  ifs.close();
739  return 0;
740 }
741 
742 
744 
745 // Interface to FileIO functionality
746 int fileio_autowrite(const Data<float,4>& data, const STD_string& filename, const FileWriteOpts& opts, const Protocol* prot);
747 int fileio_autoread(Data<float,4>& data, const STD_string& filename, const FileReadOpts& opts, Protocol* prot, ProgressMeter* progmeter);
748 
749 
750 template<typename T,int N_rank>
751 int Data<T,N_rank>::autowrite(const STD_string& filename, const FileWriteOpts& opts, const Protocol* prot) const {
752  Data<float,4> filedata; convert_to(filedata);
753  return fileio_autowrite(filedata, filename, opts, prot);
754 }
755 
756 template<typename T,int N_rank>
757 int Data<T,N_rank>::autoread(const STD_string& filename, const FileReadOpts& opts, Protocol* prot, ProgressMeter* progmeter) {
758  Data<float,4> filedata;
759  int result=fileio_autoread(filedata, filename, opts, prot, progmeter);
760  if(result>0) filedata.convert_to(*this);
761  return result;
762 }
763 
764 
766 
767 template<typename T,int N_rank>
768 Data<T,N_rank>::operator tjarray<tjvector<T>,T> () const {
769  tjarray<tjvector<T>,T> a;
770  ndim nn(N_rank);
771  for(unsigned int i=0; i<N_rank; i++) nn[i]=Array<T,N_rank>::extent(i);
772  a.redim(nn);
773  for(unsigned int i=0;i<a.total();i++) a[i]=(*this)(create_index(i));
774  return a;
775 }
776 
777 
778 
780 
781 template<typename T,int N_rank>
782 template<typename T2, int N_rank2>
783 Data<T2,N_rank2>& Data<T,N_rank>::convert_to(Data<T2,N_rank2>& dst, bool autoscale) const {
784  Log<OdinData> odinlog("Data","convert_to");
785 
786  TinyVector<int,N_rank2> newshape; newshape=1;
787  for(int i=0; i<(N_rank-N_rank2+1); i++) { // collapse extra dims into first dim of new shape, including first dim
788  int srcindex=i;
789  if(srcindex>=0 && srcindex<N_rank) newshape(0)*=Array<T,N_rank>::extent(srcindex);
790  }
791  for(int i=0; i<(N_rank2-1); i++) { // Fill new shape with extents of original dim, except for the first dim
792  int srcindex=N_rank-N_rank2+1+i;
793  if(srcindex>=0 && srcindex<N_rank) newshape(1+i)=Array<T,N_rank>::extent(srcindex);
794  }
795 
796  // modify last dimension to account for different element sizes
797  newshape(N_rank2-1)=newshape(N_rank2-1)*Converter::get_elements((T)0)/Converter::get_elements((T2)0);
798 
799  dst.resize(newshape); // dst will be unique now
800 
801  Data<T,N_rank> src_copy(*this); // make read/write copy of this
802 
803  Converter::convert_array(src_copy.c_array(), dst.c_array(), src_copy.numElements(), dst.numElements(), autoscale);
804 
805  return dst;
806 }
807 
808 
809 // specialization which uses reference if type/rank are same and no scaling is necessary to reduce memory usage
810 template<typename T,int N_rank>
811 Data<T,N_rank>& Data<T,N_rank>::convert_to(Data<T,N_rank>& dst, bool autoscale) const {
812  Log<OdinData> odinlog("Data","convert_to");
813  if(!autoscale || !std::numeric_limits<T>::is_integer) {
814  ODINLOG(odinlog,normalDebug) << "Using reference" << STD_endl;
815  dst.reference(*this);
816  } else {
817  Data<T,N_rank> src_copy(*this); // make read/write copy of this
818  dst.resize(src_copy.shape()); // dst will be unique now
819  Converter::convert_array(src_copy.c_array(), dst.c_array(), src_copy.numElements(), dst.numElements(), autoscale);
820  }
821  return dst;
822 }
823 
825 
826 template<typename T,int N_rank>
827 TinyVector<int, N_rank> Data<T,N_rank>::create_index(unsigned long index) const {
828  return index2extent<N_rank>(Array<T,N_rank>::shape(), index);
829 }
830 
832 
833 template<typename T,int N_rank>
834 unsigned long Data<T,N_rank>::create_linear_index(const TinyVector<int, N_rank>& indexvec) const {
835  unsigned long result=0;
836  unsigned long subsize=1;
837  for(int i=N_rank-1;i>=0;i--) {
838  result+=subsize*indexvec(i);
839  subsize*=Array<T,N_rank>::extent(i);
840  }
841  return result;
842 }
843 
845 
846 template<typename T,int N_rank>
847 void Data<T,N_rank>::shift(unsigned int shift_dim, int shift) {
848  Log<OdinData> odinlog("Data","shift");
849 
850  if(!shift) return;
851 
852  if( shift_dim >= N_rank ){
853  ODINLOG(odinlog,errorLog) << "shift dimension(" << shift_dim << ") >= rank of data (" << N_rank << ") !\n";
854  return;
855  }
856 
857  int shift_extent=Array<T,N_rank>::extent(shift_dim);
858  int abs_shift=abs(shift);
859  if(shift_extent < abs_shift){
860  ODINLOG(odinlog,errorLog) << "extent(" << shift_extent << ") less than shift(" << abs_shift << ") !\n";
861  return;
862  }
863 
864  Data<T,N_rank> data_copy(Array<T,N_rank>::copy());
865 
866  TinyVector<int,N_rank> index;
867  for(unsigned int i=0; i<Array<T,N_rank>::numElements(); i++) {
868  index=create_index(i);
869  T val=data_copy(index);
870  int shiftindex=index(shift_dim)+shift;
871  if(shiftindex>=shift_extent) shiftindex-=shift_extent;
872  if(shiftindex<0) shiftindex+=shift_extent;
873  index(shift_dim)=shiftindex;
874  (*this)(index)=val;
875  }
876 
877 }
878 
879 
880 
882 
883 template<typename T,int N_rank>
884 void Data<T,N_rank>::congrid(const TinyVector<int,N_rank>& newshape, const TinyVector<float,N_rank>* subpixel_shift, bool left_to_right) {
885  Log<OdinData> odinlog("Data","congrid");
886  for(int i=0;i<N_rank;i++) {
887  ODINLOG(odinlog,normalDebug) << "oldshape/newshape(" << i << ")=" << Array<T,N_rank>::shape()(i) << "/" << newshape(i) << STD_endl;
888  if(subpixel_shift) {
889  ODINLOG(odinlog,normalDebug) << "subpixel_shift(" << i << ")=" << (*subpixel_shift)(i) << STD_endl;
890  }
891  }
892 
893  for(int irank=0; irank<N_rank; irank++) {
894  int dim=irank;
895  if(!left_to_right) dim=N_rank-1-irank;
896  float shift=0.0;
897  if(subpixel_shift) shift=(*subpixel_shift)(dim);
898  interpolate1dim(dim,newshape(dim),shift);
899  }
900 }
901 
903 
904 
905 template<typename T,int N_rank>
906 void Data<T,N_rank>::interpolate1dim(unsigned int dim, int newsize, float subpixel_shift) {
907  Log<OdinData> odinlog("Data","interpolate1dim");
908  ODINLOG(odinlog,normalDebug) << "dim/newsize=" << dim << "/" << newsize << STD_endl;
909 
910  if(newsize==Array<T,N_rank>::shape()(dim) && subpixel_shift==0.0) {
911  ODINLOG(odinlog,normalDebug) << "interpolation not required" << STD_endl;
912  return;
913  }
914 
915  if(dim>=N_rank) {
916  ODINLOG(odinlog,errorLog) << "dim is larger than N_rank" << STD_endl;
917  return;
918  }
919  if(newsize<0) {
920  ODINLOG(odinlog,errorLog) << "newsize is negative" << STD_endl;
921  return;
922  }
923 
924  Array<T,N_rank> olddata(*this);
925  olddata.makeUnique();
926 
927  TinyVector<int,N_rank> oldshape(olddata.shape());
928  int oldsize=oldshape(dim);
929 
930  TinyVector<int,N_rank> newshape(oldshape);
931  newshape(dim)=newsize;
932  Data<T,N_rank>::resize(newshape); // scope reuired by GCC 4.7
933 
934  // This holds the the shape of the subspace which is orthogonal to direction 'dim'
935  TinyVector<int,N_rank> ortho_shape(oldshape);
936  ortho_shape(dim)=1;
937 
938  // The total number of elements in the orthogonal subspace
939  unsigned long n_ortho=product(ortho_shape);
940 
941  ODINLOG(odinlog,normalDebug) << "n_ortho=" << n_ortho << STD_endl;
942 
943  TinyVector<int,N_rank> indexvec;
944  T* oldoneline=new T[oldsize];
945 
946  for(unsigned long iortho=0; iortho<n_ortho; iortho++) {
947 
948  indexvec=index2extent<N_rank>(ortho_shape,iortho);
949 
950  for(int j=0; j<oldsize; j++) {
951  indexvec(dim)=j;
952  oldoneline[j]=olddata(indexvec);
953  }
954 
955  T* newoneline=interpolate1D(oldoneline,oldsize,newsize,subpixel_shift);
956 
957  for(int j=0; j<newsize; j++) {
958  indexvec(dim)=j;
959  (*this)(indexvec)=newoneline[j];
960  }
961 
962  delete[] newoneline;
963  }
964 
965  delete[] oldoneline;
966 }
967 
968 
970 
971 template<typename T,int N_rank>
972 T* Data<T,N_rank>::c_array() {
973  Log<OdinData> odinlog("Data","c_array");
974 
975  bool need_copying=false;
976 
977 #ifndef BLITZ_REPLACEMENT
978 
979  // check storage order
980  TinyVector<int,N_rank>ord=Array<T,N_rank>::ordering();
981  for(int i=0; i<N_rank-1; i++) {
982  if(ord(i)<ord(i+1)) need_copying=true;
983  }
984 
985  // check storage direction
986  for(int i=0; i<N_rank; i++) {
987  if(!Array<T,N_rank>::isRankStoredAscending(i)) need_copying=true;
988  }
989 
990  // check for slicing
991  if(!Array<T,N_rank>::isStorageContiguous()) need_copying=true;
992 
993 #endif
994 
995  if(need_copying) {
996  ODINLOG(odinlog,normalDebug) << "need_copying" << STD_endl;
997  Data<T,N_rank> tmp(Array<T,N_rank>::shape());
998  tmp=(*this); // reference(copy()) would not be enough since it would preserve ordering
999  reference(tmp);
1000  }
1001 
1002  return Array<T,N_rank>::data();
1003 }
1004 
1005 #endif
Converter::convert_array
static void convert_array(const Src *src, Dst *dst, unsigned int srcsize, unsigned int dstsize, bool autoscale=true)
Definition: converter.h:79
Data::Data
Data(BZ_ETPARM(_bz_ArrayExpr< T_expr >) expr)
Definition: odindata/data.h:209
Data::write_asc_file
int write_asc_file(const STD_string &filename, const Array< T, N_rank > &pre=defaultArray, const Array< T, N_rank > &post=defaultArray) const
Definition: odindata/data.h:690
Data::autowrite
int autowrite(const STD_string &filename, const FileWriteOpts &opts=FileWriteOpts(), const Protocol *prot=0) const
Definition: odindata/data.h:751
Data::Data
Data(const STD_string &filename, bool readonly, const TinyVector< int, N_rank > &shape, LONGEST_INT offset=0)
Definition: odindata/data.h:502
Data::autoread
int autoread(const STD_string &filename, const FileReadOpts &opts=FileReadOpts(), Protocol *prot=0, ProgressMeter *progmeter=0)
Definition: odindata/data.h:757
Data::Data
Data(const TinyVector< int, N_rank > &dimvec, const T &val=0)
Definition: odindata/data.h:156
Data::convert_to
Data< T2, N_rank2 > & convert_to(Data< T2, N_rank2 > &dst, bool autoscale=true) const
Definition: odindata/data.h:783
Data::operator=
Data< T, N_rank > & operator=(const Data< T, N_rank > &d)
Definition: odindata/data.h:222
Data::create_linear_index
unsigned long create_linear_index(const TinyVector< int, N_rank > &indexvec) const
Definition: odindata/data.h:834
Data::read
int read(const STD_string &format, const STD_string &filename, LONGEST_INT offset=0)
Definition: odindata/data.h:577
Data::congrid
void congrid(const TinyVector< int, N_rank > &newshape, const TinyVector< float, N_rank > *subpixel_shift=0, bool left_to_right=false)
Definition: odindata/data.h:884
Data::Data
Data(const tjarray< tjvector< T >, T > &a)
Definition: odindata/data.h:195
Data::write
int write(const STD_string &filename, fopenMode mode=overwriteMode) const
Definition: odindata/data.h:643
Data::Data
Data(const Data< T, N_rank > &d)
Definition: odindata/data.h:185
Data::is_filemapped
bool is_filemapped() const
Definition: odindata/data.h:357
Data::c_array
T * c_array()
Definition: odindata/data.h:972
Data::shift
void shift(unsigned int shift_dim, int shift)
Definition: odindata/data.h:847
Data::write
int write(const STD_string &filename, bool autoscale=true) const
Definition: odindata/data.h:671
Data::Data
Data(const Array< T, N_rank > &a)
Definition: odindata/data.h:190
Data::read
int read(const STD_string &filename, LONGEST_INT offset=0)
Definition: odindata/data.h:596
Data::write
int write(const STD_string &format, const STD_string &filename, bool autoscale=true) const
Definition: odindata/data.h:626
Data::read_asc_file
int read_asc_file(const STD_string &filename)
Definition: odindata/data.h:727
Data::create_index
TinyVector< int, N_rank > create_index(unsigned long index) const
Definition: odindata/data.h:827
Log
Definition: tjlog.h:218
Mutex
Definition: tjthread.h:46
Profiler::get_memory_usage
static STD_string get_memory_usage()
Protocol
Protocol proxy.
Definition: protocol.h:33
ndim
Definition: tjarray.h:41
ndim::dim
unsigned long dim() const
Definition: tjarray.h:68
ndim::add_dim
ndim & add_dim(unsigned long e, bool first=false)
tjarray::total
unsigned long total() const
tjarray::redim
tjarray< V, T > & redim(const ndim &nn)
secureInv
double secureInv(double denominator)
Definition: tjtools.h:271
modestring
const char * modestring(fopenMode mode)
filesize
LONGEST_INT filesize(const char *filename)
secureDivision
double secureDivision(double numerator, double denominator)
fileunmap
void fileunmap(int fd, void *start, LONGEST_INT nbytes, LONGEST_INT offset)
filemap
void * filemap(const STD_string &filename, LONGEST_INT nbytes, LONGEST_INT offset, bool readonly, int &fd)
lasterr
const char * lasterr()
fopenMode
fopenMode
Definition: tjtools.h:107
phase
fvector phase(const cvector &cv)
rmfile
int rmfile(const char *fname)
Generated on Fri Oct 15 2021 11:36:53 for ODIN by doxygen 1.9.1