odindata/data.h
1 /***************************************************************************
2  data.h - description
3  -------------------
4  begin : Fri Apr 6 2001
5  copyright : (C) 2000-2015 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 #include <tjutils/tjfeedback.h>
27 
28 #ifdef ODIN_DEBUG
29 #define BZ_DEBUG
30 #endif
31 
32 
33 
34 
35 #include <blitz/blitz.h> // include 1st to overwrite mutex defines
36 
37 // undef mutex defines
38 #ifdef BZ_MUTEX_DECLARE
39 #undef BZ_MUTEX_DECLARE
40 #endif
41 #ifdef BZ_MUTEX_INIT
42 #undef BZ_MUTEX_INIT
43 #endif
44 #ifdef BZ_MUTEX_LOCK
45 #undef BZ_MUTEX_LOCK
46 #endif
47 #ifdef BZ_MUTEX_UNLOCK
48 #undef BZ_MUTEX_UNLOCK
49 #endif
50 #ifdef BZ_MUTEX_DESTROY
51 #undef BZ_MUTEX_DESTROY
52 #endif
53 
54 // do not use mutex in Blitz, even if pthread is available
55 // to get consistent behaviour across all platforms
56 #define BZ_MUTEX_DECLARE(name)
57 #define BZ_MUTEX_INIT(name)
58 #define BZ_MUTEX_LOCK(name)
59 #define BZ_MUTEX_UNLOCK(name)
60 #define BZ_MUTEX_DESTROY(name)
61 
62 
63 
64 #include <blitz/array.h>
65 #ifdef HAVE_BLITZ_TINYVEC_ET
66 #include <blitz/tinyvec-et.h> // to use arithmetics of TinyVectors
67 #endif
68 using namespace blitz;
69 
70 
71 #include <tjutils/tjtypes.h>
72 #include <tjutils/tjprofiler.h>
73 
74 #include <odinpara/protocol.h>
75 
76 #include <odindata/converter.h>
77 #include <odindata/fileio_opts.h>
78 
79 
80 // some additional macros for complex numbers
81 BZ_DECLARE_FUNCTION_RET (float2real, STD_complex);
82 BZ_DECLARE_FUNCTION_RET (float2imag, STD_complex);
83 BZ_DECLARE_FUNCTION_RET (creal, float)
84 BZ_DECLARE_FUNCTION_RET (cimag, float)
85 BZ_DECLARE_FUNCTION_RET (cabs, float)
86 BZ_DECLARE_FUNCTION_RET (phase, float)
87 BZ_DECLARE_FUNCTION (logc)
88 BZ_DECLARE_FUNCTION (expc)
89 BZ_DECLARE_FUNCTION (conjc)
90 
91 // Other useful macros
92 BZ_DECLARE_FUNCTION2(secureDivision)
93 BZ_DECLARE_FUNCTION (secureInv)
94 
95 
96 #ifdef STL_REPLACEMENT
97 BZ_DECLARE_ARRAY_ET_SCALAR_OPS(tjstd_complex)
98 #endif
99 
100 
106 
108 struct FileMapHandle {
109  FileMapHandle() : fd(-1), offset(0), refcount(1) {}
110  int fd;
111  LONGEST_INT offset; // in bytes
112  int refcount;
113  Mutex mutex; // protect refcount
114 };
115 
117 
128 template <typename T, int N_rank>
129 class Data : public Array<T,N_rank> {
130 
131 public:
132 
133 
137  Data(const TinyVector<int,N_rank>& dimvec, const T& val=0) : Array<T,N_rank>(dimvec), fmap(0) {(*this)=val;}
138 
139 
140  // Variable extent constructors
141  Data() : fmap(0) {}
142  Data(int extent1) : Array<T,N_rank>(extent1), fmap(0) {}
143  Data(int extent1, int extent2) : Array<T,N_rank>(extent1,extent2), fmap(0) {}
144  Data(int extent1, int extent2,int extent3) : Array<T,N_rank>(extent1,extent2,extent3), fmap(0) {}
145  Data(int extent1, int extent2,int extent3,int extent4) : Array<T,N_rank>(extent1,extent2,extent3,extent4), fmap(0) {}
146  Data(int extent1, int extent2,int extent3,int extent4,int extent5) : Array<T,N_rank>(extent1,extent2,extent3,extent4,extent5), fmap(0) {}
147  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) {}
148  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) {}
149  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) {}
150  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) {}
151 
152 
161  Data(const STD_string& filename, bool readonly, const TinyVector<int,N_rank>& shape, LONGEST_INT offset=0 );
162 
167 
171  Data(const Array<T,N_rank>& a) : Array<T,N_rank>(a), fmap(0) {}
172 
176  Data(const tjarray<tjvector<T>,T>& a) : fmap(0) { (*this)=a;}
177 
181  template<class T_expr>
182  Data(BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr) : Array<T,N_rank>(expr), fmap(0) {}
183 
184 
188  ~Data();
189 
190 
194  Data<T,N_rank>& operator = (const Data<T,N_rank>& d) {Array<T,N_rank>::operator=(d); return *this;}
195 
199  Data<T,N_rank>& operator = (const Array<T,N_rank>& a) {Array<T,N_rank>::operator=(a); return *this;}
200 
204  Data<T,N_rank>& operator = (const T& val) {Array<T,N_rank>::operator=(val); return *this;}
205 
209  template<class T_expr>
210  inline Data<T,N_rank>& operator = (BZ_ETPARM(_bz_ArrayExpr<T_expr>) expr) {
211  Data<T,N_rank>::reference(Data<T,N_rank>(expr)); // forwarding to constructor
212  return *this;
213  }
214 
218  Data<T,N_rank>& operator = (const tjarray<tjvector<T>,T>& a);
219 
220 
221  // resolve ambiguities, appears to be obsolete with blitz-0.9 and later versions
222 // 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));}
223 // 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));}
224 // 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));}
225 // 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));}
226 
227 
228  // resolve ambiguities, appears to be obsolete with blitz-0.9 and later versions
229 // Array<T,N_rank> operator + (const T& v) const {return Array<T,N_rank>(Array<T,N_rank>(*this)+v);}
230 // Array<T,N_rank> operator - (const T& v) const {return Array<T,N_rank>(Array<T,N_rank>(*this)-v);}
231 // Array<T,N_rank> operator * (const T& v) const {return Array<T,N_rank>(Array<T,N_rank>(*this)*v);}
232 // Array<T,N_rank> operator / (const T& v) const {return Array<T,N_rank>(Array<T,N_rank>(*this)/v);}
233 
234 
235 
242  int read(const STD_string& format,const STD_string& filename, LONGEST_INT offset=0);
243 
244 
250  template<typename T2>
251  int read(const STD_string& filename, LONGEST_INT offset=0);
252 
253 
260  int write(const STD_string& format,const STD_string& filename, bool autoscale=true) const;
261 
262 
268  template<typename T2>
269  int write(const STD_string& filename, bool autoscale=true) const;
270 
271 
272 
278  int write(const STD_string& filename, fopenMode mode=overwriteMode) const;
279 
289  int write_asc_file(const STD_string& filename, const Array<T,N_rank>& pre=defaultArray, const Array<T,N_rank>& post=defaultArray) const;
290 
295  int read_asc_file(const STD_string& filename);
296 
297 
303  int autowrite(const STD_string& filename, const FileWriteOpts& opts=FileWriteOpts(), const Protocol* prot=0) const;
304 
311  int autoread(const STD_string& filename, const FileReadOpts& opts=FileReadOpts(), Protocol* prot=0, ProgressMeter* progmeter=0);
312 
316  bool is_filemapped() const {return fmap;}
317 
321  operator tjarray<tjvector<T>,T> () const;
322 
323 
328  TinyVector<int, N_rank> create_index(unsigned long index) const;
329 
333  unsigned long create_linear_index(const TinyVector<int, N_rank>& indexvec) const;
334 
340  void shift(unsigned int shift_dim, int shift);
341 
351  void congrid(const TinyVector<int,N_rank>& newshape, const TinyVector<float,N_rank>* subpixel_shift=0, bool left_to_right=false);
352 
353 
360  template<typename T2, int N_rank2> Data<T2,N_rank2>& convert_to(Data<T2,N_rank2>& dst, bool autoscale=true) const;
361 
362  // specialization which uses reference if type/rank are same and scale/offset==1.0/0.0 to reduce memory usage
363  Data<T,N_rank>& convert_to(Data<T,N_rank>& dst, bool autoscale=true) const;
364 
365 
369  T* c_array();
370 
371 
372  // Reimplemented Blitz::Array functions
373  void reference(const Data<T,N_rank>& d);
374 
375 
376  // dummy objects used for default arguments
377  static Array<T,N_rank> defaultArray;
378 
379 
380  private:
381 
382  void interpolate1dim(unsigned int dim, int newsize, float subpixel_shift);
383 
384  void detach_fmap();
385 
386  FileMapHandle* fmap;
387 
388 };
389 
390 template <typename T, int N_rank>
391 Array<T,N_rank> Data<T,N_rank>::defaultArray;
392 
393 
397 
399 #ifdef STREAM_REPLACEMENT
400 
401 template<typename T, int N_rank>
402 inline STD_ostream& operator << (STD_ostream& s,const Array<T,N_rank>&) {return s;}
403 
404 template<typename T, int N_rank>
405 inline STD_ostream& operator << (STD_ostream& s,const TinyVector<T,N_rank>&) { return s;}
406 
407 inline STD_ostream& operator << (STD_ostream& s,const Range&) { return s;}
408 
409 #endif
410 
412 
413 
414 /*
415  * Create n-dim index from linear 'index' for a given 'shape'
416  */
417 template<int N_rank>
418 TinyVector<int,N_rank> index2extent(const TinyVector<int,N_rank>& shape, unsigned int index) {
419  TinyVector<int,N_rank> result;
420  unsigned int temp=index;
421  for(int i=N_rank-1;i>=0;i--){
422  result(i)=temp%shape(i);
423  temp=temp/shape(i);
424  }
425  return result;
426 }
427 
428 
430 
431 
432 // leave this function outside of class Data, otherwise it screws up
433 // when using optimization
434 template<typename T, int N_rank, typename T2>
435 void convert_from_ptr(Data<T,N_rank>& dst, const T2* src, const TinyVector<int, N_rank>& shape, bool autoscale=true) {
436  Log<OdinData> odinlog("Data","convert_from_ptr");
437 
438  int dstsize=product(shape);
439  int srcsize=dstsize*Converter::get_elements((T)0)/Converter::get_elements((T2)0);
440  ODINLOG(odinlog,normalDebug) << "dstsize/srcsize=" << dstsize << "/" << srcsize << STD_endl;
441 
442  dst.resize(shape);
443 
444  Converter::convert_array(src, dst.c_array(), srcsize, dstsize, autoscale);
445 }
446 
447 // specialization in case the type is the same
448 template<typename T, int N_rank>
449 void convert_from_ptr(Data<T,N_rank>& dst, const T* src, const TinyVector<int, N_rank>& shape) {
450  dst.reference(Array<T,N_rank>((T*)src, shape, duplicateData));
451 }
452 
453 
455 
456 template <typename T, int N_rank>
457 Data<T,N_rank>::Data(const STD_string& filename, bool readonly, const TinyVector<int,N_rank>& shape, LONGEST_INT offset ) : fmap(new FileMapHandle) {
458  T* ptr=(T*)filemap(filename, (LONGEST_INT)product(shape)*sizeof(T), offset, readonly, fmap->fd);
459  if(ptr && (fmap->fd)>=0) { // only negative file descriptor indicates filemap failure
460  Array<T,N_rank>::reference(Array<T,N_rank>(ptr, shape, neverDeleteData));
461  fmap->offset=offset;
462  } else {
463  delete fmap;
464  fmap=0;
465  }
466 }
467 
469 
470 template <typename T, int N_rank>
472  Log<OdinData> odinlog("Data","detach_fmap");
473  if(fmap) {
474  fmap->mutex.lock();
475  (fmap->refcount)--;
476  ODINLOG(odinlog,normalDebug) << "fd/offset/refcount=" << fmap->fd << "/" << fmap->offset << "/" << fmap->refcount << STD_endl;
477  if(!(fmap->refcount)) {
478  fileunmap(fmap->fd, Array<T,N_rank>::data(), (LONGEST_INT)Array<T,N_rank>::size()*sizeof(T), fmap->offset);
479  fmap->mutex.unlock();
480  delete fmap;
481  fmap=0;
482  }
483  if(fmap) fmap->mutex.unlock();
484  }
485 }
486 
488 
489 template <typename T, int N_rank>
491  detach_fmap();
492 }
493 
495 
496 template <typename T, int N_rank>
498  Log<OdinData> odinlog("Data","reference");
499  detach_fmap();
500  fmap=d.fmap;
501  if(fmap) {
502  MutexLock lock(fmap->mutex);
503  (fmap->refcount)++;
504  ODINLOG(odinlog,normalDebug) << "fmap->refcount=" << fmap->refcount << STD_endl;
505  }
506  Array<T,N_rank>::reference(d);
507 }
508 
510 
511 template <typename T, int N_rank>
513  Log<OdinData> odinlog("Data","=");
514  if( a.dim() <= N_rank ) {
515  ndim nn=a.get_extent();
516  int npad=N_rank-nn.dim();
517  ODINLOG(odinlog,normalDebug) << "npad=" << npad << STD_endl;
518  for(int ipad=0; ipad<npad; ipad++) nn.add_dim(1, true); // pad at front
519  ODINLOG(odinlog,normalDebug) << "nn(pad)" << nn << STD_endl;
520  TinyVector<int,N_rank> tv;
521  for(unsigned int i=0; i<N_rank; i++) tv(i)=nn[i];
522  Data<T,N_rank>::resize(tv); // scope reuired by GCC 4.7
523  for(unsigned int i=0;i<a.total();i++) (*this)(create_index(i))=a[i];
524  } else ODINLOG(odinlog,errorLog) << "dimension mismatch: this=" << N_rank<< " < tjarray=" << a.dim() << STD_endl;
525  return *this;
526 }
527 
528 
530 
531 template<typename T,int N_rank>
532 int Data<T,N_rank>::read(const STD_string& format,const STD_string& filename, LONGEST_INT offset) {
533  Log<OdinData> odinlog("Data","read");
534  if(format==TypeTraits::type2label((u8bit)0)) return read<u8bit> (filename);
535  if(format==TypeTraits::type2label((s8bit)0)) return read<s8bit> (filename);
536  if(format==TypeTraits::type2label((u16bit)0)) return read<u16bit>(filename);
537  if(format==TypeTraits::type2label((s16bit)0)) return read<s16bit>(filename);
538  if(format==TypeTraits::type2label((u32bit)0)) return read<u32bit>(filename);
539  if(format==TypeTraits::type2label((s32bit)0)) return read<s32bit>(filename);
540  if(format==TypeTraits::type2label((float)0)) return read<float> (filename);
541  if(format==TypeTraits::type2label((double)0)) return read<double>(filename);
542  ODINLOG(odinlog,errorLog) << "Unable to read file " << filename << " with data type " << format << STD_endl;
543  return -1;
544 }
545 
546 
548 
549 template<typename T,int N_rank>
550 template<typename T2>
551 int Data<T,N_rank>::read(const STD_string& filename, LONGEST_INT offset) {
552  Log<OdinData> odinlog("Data","read");
553  LONGEST_INT fsize=filesize(filename.c_str())-offset;
554  LONGEST_INT nelements_file=fsize/sizeof(T2);
555  LONGEST_INT length=product(Array<T,N_rank>::shape());
556  ODINLOG(odinlog,normalDebug) << "fsize / nelements_file / length = " << fsize << " / " << nelements_file << " / " << length << STD_endl;
557 
558  if(!length) return 0;
559 
560  if(nelements_file<length) {
561  ODINLOG(odinlog,errorLog) << "Size of file " << filename << " to small for reading" << STD_endl;
562  return -1;
563  }
564 
565  STD_string srctype=TypeTraits::type2label((T2)0);
566  STD_string dsttype=TypeTraits::type2label((T)0);
567  ODINLOG(odinlog,normalDebug) << "srctype/dsttype=" << srctype << "/" << dsttype << STD_endl;
568 
569  TinyVector<int,N_rank> fileshape(Array<T,N_rank>::shape());
570  fileshape(N_rank-1) *= (Converter::get_elements((T)0)/Converter::get_elements((T2)0)); // Adjust extent of last dim for copmlex data
571  Data<T2,N_rank> filedata(filename, true, fileshape, offset);
572  filedata.convert_to(*this);
573 
574  return 0;
575 }
576 
577 
579 
580 template<typename T,int N_rank>
581 int Data<T,N_rank>::write(const STD_string& format, const STD_string& filename, bool autoscale) const {
582  Log<OdinData> odinlog("Data","write");
583  if(format==TypeTraits::type2label((u8bit)0)) return write<u8bit> (filename,autoscale);
584  if(format==TypeTraits::type2label((s8bit)0)) return write<s8bit> (filename,autoscale);
585  if(format==TypeTraits::type2label((u16bit)0)) return write<u16bit>(filename,autoscale);
586  if(format==TypeTraits::type2label((s16bit)0)) return write<s16bit>(filename,autoscale);
587  if(format==TypeTraits::type2label((u32bit)0)) return write<u32bit>(filename,autoscale);
588  if(format==TypeTraits::type2label((s32bit)0)) return write<s32bit>(filename,autoscale);
589  if(format==TypeTraits::type2label((float)0)) return write<float> (filename,autoscale);
590  if(format==TypeTraits::type2label((double)0)) return write<double>(filename,autoscale);
591  ODINLOG(odinlog,errorLog) << "Unable to write file " << filename << " with data type " << format << STD_endl;
592  return -1;
593 }
594 
596 
597 template<typename T,int N_rank>
598 int Data<T,N_rank>::write(const STD_string& filename, fopenMode mode) const {
599  Log<OdinData> odinlog("Data","write");
600  if(filename=="") return 0;
601 
602  FILE* file_ptr=ODIN_FOPEN(filename.c_str(),modestring(mode));
603 
604  if(file_ptr==NULL) {
605  ODINLOG(odinlog,errorLog) << "unable to create/open file >" << filename << "< - " << lasterr() << STD_endl;
606  return -1;
607  }
608 
609  Data<T,N_rank> data_copy(*this); // for contig memory
610 
611  LONGEST_INT nmemb=Array<T,N_rank>::numElements();
612  LONGEST_INT count=fwrite(data_copy.c_array(),sizeof(T),nmemb,file_ptr);
613  if(count!=nmemb) {
614  ODINLOG(odinlog,errorLog) << "unable to fwrite to file >" << filename << "< - " << lasterr() << STD_endl;
615  return -1;
616  }
617 
618  if(file_ptr!=NULL) fclose(file_ptr);
619  return 0;
620 }
621 
623 
624 template<typename T,int N_rank>
625 template<typename T2>
626  int Data<T,N_rank>::write(const STD_string& filename, bool autoscale) const {
627  Log<OdinData> odinlog("Data","write");
628 
629  rmfile(filename.c_str()); // remove old file to get new file size with mmap
630 
631  ODINLOG(odinlog,normalDebug) << "mem(pre): " << Profiler::get_memory_usage() << STD_endl;
632  Data<T2,N_rank> converted_data; convert_to(converted_data,autoscale);
633  ODINLOG(odinlog,normalDebug) << "mem(cnv): " << Profiler::get_memory_usage() << STD_endl;
634  Data<T2,N_rank> filedata(filename,false,converted_data.shape());
635  ODINLOG(odinlog,normalDebug) << "mem(map): " << Profiler::get_memory_usage() << STD_endl;
636  filedata=converted_data;
637 
638  return 0;
639 }
640 
641 
643 
644 template<typename T,int N_rank>
645 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 {
646  bool have_pre=false;
647  bool have_post=false;
648 
649  Data<T,N_rank> pre_data(pre);
650  Data<T,N_rank> post_data(post);
651  unsigned int n=Array<T,N_rank>::numElements();
652 
653  if(pre_data.numElements()==n) have_pre=true;
654  if(post_data.numElements()==n) have_post=true;
655 
656  STD_ofstream ofs(filename.c_str());
657  if(ofs.bad()) return -1;
658 
659  T val;
660 
661  for(unsigned int i=0; i<n; i++) {
662  if(have_pre) {
663  val=pre_data(pre_data.create_index(i));
664  ofs << val << " ";
665  }
666  val=(*this)(create_index(i));
667  ofs << val;
668  if(have_post) {
669  val=post_data(post_data.create_index(i));
670  ofs << " " << val;
671  }
672  ofs << "\n";
673  }
674 
675  ofs.close();
676  return 0;
677 }
678 
680 
681 template<typename T,int N_rank>
682 int Data<T,N_rank>::read_asc_file(const STD_string& filename) {
683  STD_ifstream ifs(filename.c_str());
684  if(ifs.bad()) return -1;
685 
686  STD_string valstr;
687  for(unsigned int i=0; i<Array<T,N_rank>::numElements(); i++) {
688  if(ifs.bad()) return -1;
689  ifs >> valstr;
690  TypeTraits::string2type(valstr,(*this)(create_index(i)));
691  }
692 
693  ifs.close();
694  return 0;
695 }
696 
697 
699 
700 // Interface to FileIO functionality
701 int fileio_autowrite(const Data<float,4>& data, const STD_string& filename, const FileWriteOpts& opts, const Protocol* prot);
702 int fileio_autoread(Data<float,4>& data, const STD_string& filename, const FileReadOpts& opts, Protocol* prot, ProgressMeter* progmeter);
703 
704 
705 template<typename T,int N_rank>
706 int Data<T,N_rank>::autowrite(const STD_string& filename, const FileWriteOpts& opts, const Protocol* prot) const {
707  Data<float,4> filedata; convert_to(filedata);
708  return fileio_autowrite(filedata, filename, opts, prot);
709 }
710 
711 template<typename T,int N_rank>
712 int Data<T,N_rank>::autoread(const STD_string& filename, const FileReadOpts& opts, Protocol* prot, ProgressMeter* progmeter) {
713  Data<float,4> filedata;
714  int result=fileio_autoread(filedata, filename, opts, prot, progmeter);
715  if(result>0) filedata.convert_to(*this);
716  return result;
717 }
718 
719 
721 
722 template<typename T,int N_rank>
724  tjarray<tjvector<T>,T> a;
725  ndim nn(N_rank);
726  for(unsigned int i=0; i<N_rank; i++) nn[i]=Array<T,N_rank>::extent(i);
727  a.redim(nn);
728  for(unsigned int i=0;i<a.total();i++) a[i]=(*this)(create_index(i));
729  return a;
730 }
731 
732 
733 
735 
736 template<typename T,int N_rank>
737 template<typename T2, int N_rank2>
738 Data<T2,N_rank2>& Data<T,N_rank>::convert_to(Data<T2,N_rank2>& dst, bool autoscale) const {
739  Log<OdinData> odinlog("Data","convert_to");
740 
741  TinyVector<int,N_rank2> newshape; newshape=1;
742  for(int i=0; i<(N_rank-N_rank2+1); i++) { // collapse extra dims into first dim of new shape, including first dim
743  int srcindex=i;
744  if(srcindex>=0 && srcindex<N_rank) newshape(0)*=Array<T,N_rank>::extent(srcindex);
745  }
746  for(int i=0; i<(N_rank2-1); i++) { // Fill new shape with extents of original dim, except for the first dim
747  int srcindex=N_rank-N_rank2+1+i;
748  if(srcindex>=0 && srcindex<N_rank) newshape(1+i)=Array<T,N_rank>::extent(srcindex);
749  }
750 
751  // modify last dimension to account for different element sizes
752  newshape(N_rank2-1)=newshape(N_rank2-1)*Converter::get_elements((T)0)/Converter::get_elements((T2)0);
753 
754  dst.resize(newshape); // dst will be unique now
755 
756  Data<T,N_rank> src_copy(*this); // make read/write copy of this
757 
758  Converter::convert_array(src_copy.c_array(), dst.c_array(), src_copy.numElements(), dst.numElements(), autoscale);
759 
760  return dst;
761 }
762 
763 
764 // specialization which uses reference if type/rank are same and no scaling is necessary to reduce memory usage
765 template<typename T,int N_rank>
766 Data<T,N_rank>& Data<T,N_rank>::convert_to(Data<T,N_rank>& dst, bool autoscale) const {
767  Log<OdinData> odinlog("Data","convert_to");
768  if(!autoscale || !std::numeric_limits<T>::is_integer) {
769  ODINLOG(odinlog,normalDebug) << "Using reference" << STD_endl;
770  dst.reference(*this);
771  } else {
772  Data<T,N_rank> src_copy(*this); // make read/write copy of this
773  dst.resize(src_copy.shape()); // dst will be unique now
774  Converter::convert_array(src_copy.c_array(), dst.c_array(), src_copy.numElements(), dst.numElements(), autoscale);
775  }
776  return dst;
777 }
778 
780 
781 template<typename T,int N_rank>
782 TinyVector<int, N_rank> Data<T,N_rank>::create_index(unsigned long index) const {
783  return index2extent<N_rank>(Array<T,N_rank>::shape(), index);
784 }
785 
787 
788 template<typename T,int N_rank>
789 unsigned long Data<T,N_rank>::create_linear_index(const TinyVector<int, N_rank>& indexvec) const {
790  unsigned long totalIndex=0,subsize;
791  TinyVector<int, N_rank> nn(Array<T,N_rank>::extent());
792  for(unsigned long i=0;i<N_rank;i++) {
793  nn(i)=1;
794  subsize=product(nn);
795  if(!subsize) subsize=1;
796  totalIndex+=subsize*indexvec(i);
797  }
798  return totalIndex;
799 }
800 
802 
803 template<typename T,int N_rank>
804 void Data<T,N_rank>::shift(unsigned int shift_dim, int shift) {
805  Log<OdinData> odinlog("Data","shift");
806 
807  if(!shift) return;
808 
809  if( shift_dim >= N_rank ){
810  ODINLOG(odinlog,errorLog) << "shift dimension(" << shift_dim << ") >= rank of data (" << N_rank << ") !\n";
811  return;
812  }
813 
814  int shift_extent=Array<T,N_rank>::extent(shift_dim);
815  int abs_shift=abs(shift);
816  if(shift_extent < abs_shift){
817  ODINLOG(odinlog,errorLog) << "extent(" << shift_extent << ") less than shift(" << abs_shift << ") !\n";
818  return;
819  }
820 
821  Data<T,N_rank> data_copy(Array<T,N_rank>::copy());
822 
823  TinyVector<int,N_rank> index;
824  for(unsigned int i=0; i<Array<T,N_rank>::numElements(); i++) {
825  index=create_index(i);
826  T val=data_copy(index);
827  int shiftindex=index(shift_dim)+shift;
828  if(shiftindex>=shift_extent) shiftindex-=shift_extent;
829  if(shiftindex<0) shiftindex+=shift_extent;
830  index(shift_dim)=shiftindex;
831  (*this)(index)=val;
832  }
833 
834 }
835 
836 
837 
839 
840 template<typename T,int N_rank>
841 void Data<T,N_rank>::congrid(const TinyVector<int,N_rank>& newshape, const TinyVector<float,N_rank>* subpixel_shift, bool left_to_right) {
842  Log<OdinData> odinlog("Data","congrid");
843  for(int i=0;i<N_rank;i++) {
844  ODINLOG(odinlog,normalDebug) << "oldshape/newshape(" << i << ")=" << Array<T,N_rank>::shape()(i) << "/" << newshape(i) << STD_endl;
845  if(subpixel_shift) {
846  ODINLOG(odinlog,normalDebug) << "subpixel_shift(" << i << ")=" << (*subpixel_shift)(i) << STD_endl;
847  }
848  }
849 
850  for(int irank=0; irank<N_rank; irank++) {
851  int dim=irank;
852  if(!left_to_right) dim=N_rank-1-irank;
853  float shift=0.0;
854  if(subpixel_shift) shift=(*subpixel_shift)(dim);
855  interpolate1dim(dim,newshape(dim),shift);
856  }
857 }
858 
860 
861 
862 template<typename T,int N_rank>
863 void Data<T,N_rank>::interpolate1dim(unsigned int dim, int newsize, float subpixel_shift) {
864  Log<OdinData> odinlog("Data","interpolate1dim");
865  ODINLOG(odinlog,normalDebug) << "dim/newsize=" << dim << "/" << newsize << STD_endl;
866 
867  if(newsize==Array<T,N_rank>::shape()(dim) && subpixel_shift==0.0) {
868  ODINLOG(odinlog,normalDebug) << "interpolation not required" << STD_endl;
869  return;
870  }
871 
872  if(dim>=N_rank) {
873  ODINLOG(odinlog,errorLog) << "dim is larger than N_rank" << STD_endl;
874  return;
875  }
876  if(newsize<0) {
877  ODINLOG(odinlog,errorLog) << "newsize is negative" << STD_endl;
878  return;
879  }
880 
881  Array<T,N_rank> olddata(*this);
882  olddata.makeUnique();
883 
884  TinyVector<int,N_rank> oldshape(olddata.shape());
885  int oldsize=oldshape(dim);
886 
887  TinyVector<int,N_rank> newshape(oldshape);
888  newshape(dim)=newsize;
889  Data<T,N_rank>::resize(newshape); // scope reuired by GCC 4.7
890 
891  // This holds the the shape of the subspace which is orthogonal to direction 'dim'
892  TinyVector<int,N_rank> ortho_shape(oldshape);
893  ortho_shape(dim)=1;
894 
895  // The total number of elements in the orthogonal subspace
896  unsigned long n_ortho=product(ortho_shape);
897 
898  ODINLOG(odinlog,normalDebug) << "n_ortho=" << n_ortho << STD_endl;
899 
900  TinyVector<int,N_rank> indexvec;
901  T* oldoneline=new T[oldsize];
902 
903  for(unsigned long iortho=0; iortho<n_ortho; iortho++) {
904 
905  indexvec=index2extent<N_rank>(ortho_shape,iortho);
906 
907  for(int j=0; j<oldsize; j++) {
908  indexvec(dim)=j;
909  oldoneline[j]=olddata(indexvec);
910  }
911 
912  T* newoneline=interpolate1D(oldoneline,oldsize,newsize,subpixel_shift);
913 
914  for(int j=0; j<newsize; j++) {
915  indexvec(dim)=j;
916  (*this)(indexvec)=newoneline[j];
917  }
918 
919  delete[] newoneline;
920  }
921 
922  delete[] oldoneline;
923 }
924 
925 
927 
928 template<typename T,int N_rank>
930  Log<OdinData> odinlog("Data","c_array");
931 
932  bool need_copying=false;
933 
934  // check storage order
935  TinyVector<int,N_rank>ord=Array<T,N_rank>::ordering();
936  for(int i=0; i<N_rank-1; i++) {
937  if(ord(i)<ord(i+1)) need_copying=true;
938  }
939 
940  // check storage direction
941  for(int i=0; i<N_rank; i++) {
942  if(!Array<T,N_rank>::isRankStoredAscending(i)) need_copying=true;
943  }
944 
945  // check for slicing
946  if(!Array<T,N_rank>::isStorageContiguous()) need_copying=true;
947 
948  if(need_copying) {
949  ODINLOG(odinlog,normalDebug) << "need_copying" << STD_endl;
950  Data<T,N_rank> tmp(Array<T,N_rank>::shape());
951  tmp=(*this); // reference(copy()) would not be enough since it would preserve ordering
952  reference(tmp);
953  }
954 
955  return Array<T,N_rank>::data();
956 }
957 
958 #endif
TinyVector< int, N_rank > create_index(unsigned long index) const
void shift(unsigned int shift_dim, int shift)
LONGEST_INT filesize(const char *filename)
int read(const STD_string &format, const STD_string &filename, LONGEST_INT offset=0)
const char * lasterr()
Data(const Array< T, N_rank > &a)
void * filemap(const STD_string &filename, LONGEST_INT nbytes, LONGEST_INT offset, bool readonly, int &fd)
static void convert_array(const Src *src, Dst *dst, unsigned int srcsize, unsigned int dstsize, bool autoscale=true)
Definition: converter.h:79
Definition: tjlog.h:218
tjarray< V, T > & redim(const ndim &nn)
Definition: tjthread.h:46
ndim & add_dim(unsigned long e, bool first=false)
double secureInv(double denominator)
Definition: tjtools.h:271
int write_asc_file(const STD_string &filename, const Array< T, N_rank > &pre=defaultArray, const Array< T, N_rank > &post=defaultArray) const
Data(BZ_ETPARM(_bz_ArrayExpr< T_expr >) expr)
void congrid(const TinyVector< int, N_rank > &newshape, const TinyVector< float, N_rank > *subpixel_shift=0, bool left_to_right=false)
int write(const STD_string &str, const STD_string &filename, fopenMode mode=overwriteMode)
unsigned long dim() const
Definition: tjarray.h:68
void fileunmap(int fd, void *start, LONGEST_INT nbytes, LONGEST_INT offset)
static STD_string get_memory_usage()
Data< T2, N_rank2 > & convert_to(Data< T2, N_rank2 > &dst, bool autoscale=true) const
Data(const Data< T, N_rank > &d)
int rmfile(const char *fname)
fopenMode
Definition: tjtools.h:107
bool is_filemapped() const
Protcol proxy.
Definition: protocol.h:33
T * c_array()
int autoread(const STD_string &filename, const FileReadOpts &opts=FileReadOpts(), Protocol *prot=0, ProgressMeter *progmeter=0)
unsigned long create_linear_index(const TinyVector< int, N_rank > &indexvec) const
Data(const tjarray< tjvector< T >, T > &a)
fvector phase(const cvector &cv)
double secureDivision(double numerator, double denominator)
Definition: tjarray.h:41
Data< T, N_rank > & operator=(const Data< T, N_rank > &d)
int read_asc_file(const STD_string &filename)
const char * modestring(fopenMode mode)
unsigned long total() const
int autowrite(const STD_string &filename, const FileWriteOpts &opts=FileWriteOpts(), const Protocol *prot=0) const
int write(const STD_string &format, const STD_string &filename, bool autoscale=true) const
Data(const TinyVector< int, N_rank > &dimvec, const T &val=0)