databuffer.cpp 44.1 KB
Newer Older
1
// Copyright 2017-2021, Schlumberger
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//      http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include "databuffer.h"
#include "roundandclip.h"
#include "../exception.h"
#include "environment.h"
19
#include "fancy_timers.h"
20
21
22
23
24
25
26
27
28
29

#include <cstdint>
#include <array>
#include <algorithm>
#include <vector>
#include <memory.h>
#include <typeinfo>
#include <sstream>
#include <cmath>
#include <iostream>
Paal Kvamme's avatar
Paal Kvamme committed
30
#include <iomanip>
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

namespace InternalZGY {
#if 0
}
#endif

namespace {
#if 0
}
#endif

namespace {
  /**
   * \brief Make a fake shared_ptr from a plain unsafe pointer.
   *
   * Used when a method or a data member expects a smart pointer while
   * the application has an unsafe pointer is guaranteed to remain
   * valid long enough. The returned smart pointer does nothing on
   * cleanup.
   *
   * CAVEAT: Using this method indicates a code smell.
   * Note: The method is duplicated in bulk.cpp, databuffer.cpp,
   * and api.cpp. There isn't much point in consolidatimg them.
   */
  template<typename T>
  static inline std::shared_ptr<T>
  fake_shared(T* p)
  {
    return std::shared_ptr<T>(p, [](T* p){});
  }

  /**
   * \brief Instrumenting DataBufferNd for performance measurements.
   *
   * The following lists functions that are public, potentially
   * expensive, and are not called from other members of DataBuffer
   * which might cause us to count the same task twice. Beware that
   * I did not check this list very thoroughly.
   *
   * Declaring all timers in one spot is not really needed, at least
71
   * not when using static instances that report when the app exits.
72
73
   * If I need a report each time a DataBuffer goes out of scope
   * it would be different.
74
   *
75
76
77
   * For slightly more advanced usage, a global cleartimers() can
   * be called to output all collected details at that point.
   *
78
   * \details Thread safety: Safe because SummaryPrintingTimerEx is
79
   * thread safe when used correctly.
80
81
82
83
   */
  class AdHocTimers
  {
  public:
84
85
86
87
88
89
90
91
92
    SummaryPrintingTimerEx ctor;
    SummaryPrintingTimerEx allsame;
    SummaryPrintingTimerEx fill;
    SummaryPrintingTimerEx range;
    SummaryPrintingTimerEx clone;
    SummaryPrintingTimerEx scaletof;
    SummaryPrintingTimerEx scaletos;
    SummaryPrintingTimerEx copyscalar;
    SummaryPrintingTimerEx copysubset;
93
    AdHocTimers()
94
95
96
97
98
99
100
101
102
      : ctor      ("DBNd.ctor")
      , allsame   ("DBNd.allsame")
      , fill      ("DBNd.fill")
      , range     ("DBNd.range")
      , clone     ("DBNd.clone")
      , scaletof  ("DBNd.scaletof")
      , scaletos  ("DBNd.scaletos")
      , copyscalar("DBNd.copyscalar")
      , copysubset("DBNd.copysubset")
103
104
105
106
107
108
    {}
    static AdHocTimers& instance()
    {
      static AdHocTimers instance_;
      return instance_;
    }
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
    void cleartimers(bool show) {
      if (show) {
        copysubset.print();
        copyscalar.print();
        scaletos.print();
        scaletof.print();
        clone.print();
        range.print();
        fill.print();
        allsame.print();
        ctor.print();
      }
      else {
        copysubset.reset();
        copyscalar.reset();
        scaletos.reset();
        scaletof.reset();
        clone.reset();
        range.reset();
        fill.reset();
        allsame.reset();
        ctor.reset();
      }
    }
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
  };
}

/**
 * Copy subset of one array into a subset of another.
 *
 * In typical usage, either source or destination will be a single
 * brick read from or being written to the file. "orig" is the
 * survey-relative position corresponding to the first sample in the
 * memory buffer. The other buffer is typically something the user
 * supplied. With "orig" being the first sample the user wants to read
 * or write.
 *
 * @param ndim      Highest dimension to consider.
 *
 * @param srcorig   Survey coordinates of the first sample in the destination.
 * @param srcsize   Size of the source array.
 * @param srcstride The stride to use, per dimension, for the source array.
 * @param srcbuf    Pointer to buffer holding source array.
 *
 * @param dstorig   Survey coordinates of the first sample in the destination.
 * @param dstsize   Size of the destination array.
 * @param dststride The stride to use, per dimension, for the destination array.
 * @param dstbuf    Pointer to buffer holding destination array.
 *
 * @param cpyorig   if non-null: Don't copy data before this, even if possible.
 * @param cpysize   If non-null: Don't copy past cpyorig + cpysize, even if
 *                  possible. If setting either of cpyorig or copysize then
 *                  both should be set.
 *
 * Based on CopySubset() in ArrayConvert.cpp in the old accessor.
 * This can also replace Permute(), as that one is
 * apparently just a special case with
 *
 *   srcorig = dstorig = cpyorig = (0,0,0);
 *   srcsize = dstsize = cpysize = size;
 *
 * As for the other code in ArrayConvert: RescaleConvert may be too general
 * for my use because OpenZGY doesn't allow conversion between arbitrary
 * representations such as between int8/int16 or int16/uint16. Big-endian
 * architectures are currently not supported so this makes many of the other
 * ArrayConvert methods moot as well.
 *
 * CopySubset needs to be templated on item type only in order to efficiently
 * handle the case where the fastest varying dimension isn't 1 in both source
 * and target. If this was not a requirement I could have used void* and an
 * explicit elementsize instead.
 *
 * TODO-Performance: The code is supposed to handle any buffer strides,
 * but is a LOT more efficient if the last dimension has step 1. This is
 * normally the case in OpenZGY but somebody might want to re-use it in
 * pther contexts.
 */
template <typename T>
static void
CopySubset(std::int32_t ndim,
           const std::int64_t* srcorig,
           const std::int64_t* srcsize,
           const std::int64_t* srcstride,
           const T* srcbuf,
           const std::int64_t* dstorig,
           const std::int64_t* dstsize,
           const std::int64_t* dststride,
           T* dstbuf,
           const std::int64_t* cpyorig = 0,
           const std::int64_t* cpysize = 0)
{
  // Sanity check.
  if (ndim < 1)
    return;

  if (cpyorig == nullptr || cpysize == nullptr) {
    // This basically ignores the "copy" limits.
    cpyorig = srcorig;
    cpysize = srcsize;
  }

  // Calculate overlap range.
  const std::int64_t beg = std::max(std::max(dstorig[0], srcorig[0]), cpyorig[0]);
  const std::int64_t end = std::min(std::min(dstorig[0] + dstsize[0], srcorig[0] + srcsize[0]), cpyorig[0] + cpysize[0]);

  if (beg >= end) {
    // degenerate case: empty overlap. Nothing to do.
  }
  else if (ndim == 1) {
    // Degenerate case: 1-dimensional.
    if (srcstride[0] == 1 && dststride[0] == 1) {
      memcpy(dstbuf + (beg - dstorig[0]),
             srcbuf + (beg - srcorig[0]),
             (end - beg)*sizeof(T));
    }
    else {
      // TODO-Test remember to test this case with a nontrivial region.
      // Copy one element at a time. Using memcpy is not as expensive
      // as it looks. As long as the compiler can see the constant
      // size it is free to replace memcpy with an assignment.
      // memcpy() is more robust wrt. misaligned data. Also, we probably
      // won't get here very often because the stride of the last dim is
      // normally 1.
      T* dst       = &dstbuf[(beg - dstorig[0])*dststride[0]];
      const T* src = &srcbuf[(beg - srcorig[0])*srcstride[0]];
      for (std::int64_t ii = beg; ii < end; ++ii,
             dst += dststride[0], src += srcstride[0])
        memcpy(dst, src, sizeof(T));
    }
  }
  else {
    // general case: recursive
    for (std::int64_t ii = beg; ii < end; ++ii) {
      CopySubset(ndim - 1,
                 srcorig + 1, srcsize + 1, srcstride + 1, srcbuf + (ii - srcorig[0])*srcstride[0],
                 dstorig + 1, dstsize + 1, dststride + 1, dstbuf + (ii - dstorig[0])*dststride[0],
                 cpyorig + 1, cpysize + 1);
    }
  }
}

/**
 * Copy a constant value into a subset of an array.
 * Apart from the source being a constant this is identical to CopySubset.
 */
template <typename T>
static void
CopyScalar(std::int32_t ndim,
           const std::int64_t* srcorig,
           const std::int64_t* srcsize,
           const T scalar,
           const std::int64_t* dstorig,
           const std::int64_t* dstsize,
           const std::int64_t* dststride,
           T* dstbuf,
           const std::int64_t* cpyorig = 0,
           const std::int64_t* cpysize = 0)
{
  // Sanity check.
  if (ndim < 1)
    return;

  if (cpyorig == nullptr || cpysize == nullptr) {
    // This basically ignores the "copy" limits.
    cpyorig = dstorig;
    cpysize = dstsize;
  }

  // Calculate overlap range.
  const std::int64_t beg = std::max(std::max(dstorig[0], srcorig[0]), cpyorig[0]);
  const std::int64_t end = std::min(std::min(dstorig[0] + dstsize[0], srcorig[0] + srcsize[0]), cpyorig[0] + cpysize[0]);

  if (beg >= end) {
    // degenerate case: empty overlap. Nothing to do.
  }
  else if (ndim == 1) {
    std::int64_t step = dststride[0];
    T* it_dst = &dstbuf[(beg - dstorig[0]) * step];
    T* it_end = &dstbuf[(end - dstorig[0]) * step];
    if (step == 1) {
      std::fill(it_dst, it_end, scalar);
    }
    else {
      for (; it_dst < it_end; it_dst += step)
          *it_dst = scalar;
    }
  }
  else {
    // general case: recursive
    for (std::int64_t ii = beg; ii < end; ++ii) {
      CopyScalar(ndim - 1,
                 srcorig + 1, srcsize + 1, scalar,
                 dstorig + 1, dstsize + 1, dststride + 1, dstbuf + (ii - dstorig[0])*dststride[0],
                 cpyorig + 1, cpysize + 1);
    }
  }
}

/**
 * Convert void* + size to a DataBuffer. See makeDataBuffer3d for details.
 *
 * TODO-Worry: Does raw need to be correctly aligned for T, or do
 * all the places where this might matter use memcpy anyway?
 * Should I perhaps automatically copy the buffer if misaligned?
 * Can it ever *be* misaligned? Probably not. makeDataBuffer() with an
 * external buffer is called when delivering a chunk of data that was
 * read or fetched from the cache. This probably came from malloc'd data
 * which was already correctly aligned. Probably. The same applies to
 * read() even though that doesn't use makeDataBuffer(). Because the type
 * is known in that case, so it calls the constructor directly.  Either way
 * the caller of read() ought to have given a correctly aligned buffer.
 */
template<typename T, int NDim>
static std::shared_ptr<DataBuffer>
_makeDataBufferT(void *raw, std::int64_t rawsize,
                   const std::array<std::int64_t,NDim> bricksize)
{
  const std::int64_t nsamples = bricksize[0] * bricksize[1] * bricksize[2];
  if (rawsize == (std::int64_t)sizeof(double) &&
        rawsize == (std::int64_t)sizeof(T) * nsamples &&
        typeid(T) != typeid(double))
    throw OpenZGY::Errors::ZgyInternalError("Ambiguous arguments to makeDataBuffer.");
  T *rawdata = static_cast<T*>(raw);
  if (rawdata == nullptr) {
    if (rawsize != 0)
        throw OpenZGY::Errors::ZgyInternalError("Invalid arguments to makeDataBuffer.");
    return std::shared_ptr<DataBufferNd<T,NDim>>
        (new DataBufferNd<T,NDim>(bricksize));
  }
  else if (rawsize == sizeof(T)) {
    T value;
    memcpy(&value, raw, sizeof(T));
    return std::make_shared<DataBufferNd<T,NDim>>(value, bricksize);
  }
  else if (rawsize == (std::int64_t)sizeof(T) * nsamples) {
    // TODO-Worry: Risk of accessing freed data.
    std::shared_ptr<T> fakeshared = fake_shared(rawdata);
    return std::shared_ptr<DataBufferNd<T,NDim>>
        (new DataBufferNd<T,NDim>(fakeshared, bricksize));
  }
  else if (rawsize == (std::int64_t)sizeof(double)) {
    // TODO-Worry: Technically ambiguous in some cases.
    double dblvalue;
    memcpy(&dblvalue, raw, sizeof(double));
    T value = static_cast<T>(dblvalue);
    return std::make_shared<DataBufferNd<T,NDim>>(value, bricksize);
  }
  else
    throw OpenZGY::Errors::ZgyInternalError("Raw byte count does not match brick size");
}

std::shared_ptr<DataBuffer>
_makeDataBuffer3d
    (void *raw, std::int64_t nbytes,
     const std::array<std::int64_t,3>& size,
     RawDataType dtype)
{
  switch (dtype) {
  case RawDataType::SignedInt8:    return _makeDataBufferT<std::int8_t,3>(raw, nbytes, size);
  case RawDataType::UnsignedInt8:  return _makeDataBufferT<std::uint8_t,3>(raw, nbytes, size);
  case RawDataType::SignedInt16:   return _makeDataBufferT<std::int16_t,3>(raw, nbytes, size);
  case RawDataType::UnsignedInt16: return _makeDataBufferT<std::uint16_t,3>(raw, nbytes, size);
  case RawDataType::SignedInt32:   return _makeDataBufferT<std::int32_t,3>(raw, nbytes, size);
  case RawDataType::UnsignedInt32: return _makeDataBufferT<std::uint32_t,3>(raw, nbytes, size);
  case RawDataType::Float32:       return _makeDataBufferT<float,3>(raw, nbytes, size);
  case RawDataType::IbmFloat32:
  default: throw OpenZGY::Errors::ZgyInternalError("Unrecognized type enum");
  }
}

} // end anonymous namespace

DataBuffer::~DataBuffer()
{
}

/**
 * Convert voidptr + nbytes + size3d to a DataBuffer.
 *
 * To avoid too much copy/paste, this method this able to create several
 * different buffer types by calling different DataBufferNd constructors.
 *
 * If voidptr is null this creates an uninitialized 3d DataBuffer with
 * the given size and type and isScalar() == false. nbytes should be 0.
 *
 * If voidptr is not null then this creates a scalar DataBuffer if
 * nbytes indicates there is room for just one sample. If there are
 * enough bytes it creates a regular DataBuffer pointing to the voidptr
 * that was passed in.
 *
 * The valid values for nbytes are:
 *
 *    0              - If and only if voidptr is null.
 *    sizeof(T)      - (T*)voidptr contains the scalar to use.
 *    sizeof(double) - (double*)voidptr contains the scalar to use,
 *    size[0]*size[1]*size[2]*sizeof(T) - (T*)voidptr is used as-is.
 *
 * Any other value will raise an exception.
 *
 * Technically the scheme is ambiguous if sizeof(T) is 8 bytes and
 * size indicates less than 8 bytes. Cannot know for sure whether the
 * buffer points to a double or to 8/4/2/1 element of T. This should
 * never happen in practice. If worried, remove the "double" feature
 * or introduce a new is_double flag.
 *
 * The scheme also cannot distinguish between a scaler size 1x1x1 and
 * a regular buffer size 1x1x1. But those two are practically equivalent.
 *
 * Caveat: If the data pointer is not correctly aligned to the type of data
 * then the returned DataBuffer will also have a misaligned pointer.
 *
 * Caveat: Since the user can pass in a void* buffer, this will not be
 * reference counted and could go out of scope beore the databuffer
 * does. This might change. To mitigate this somehow, if the code knows
 * the buffer is no longer valid it can call databuffer.clear().
 */
std::shared_ptr<DataBuffer>
DataBuffer::makeDataBuffer3d
    (void *raw, std::int64_t nbytes,
     const std::array<std::int64_t,3>& size,
     RawDataType dtype)
{
  std::shared_ptr<DataBuffer> result =
    _makeDataBuffer3d(raw, nbytes, size, dtype);
  if (result->datatype() != dtype || result->size3d() != size)
    throw OpenZGY::Errors::ZgyInternalError("makeDataBuffer3d made a mistake.");
  return result;
}

template <typename T, int NDim>
DataBufferNd<T,NDim>::DataBufferNd(T scalar, const std::array<std::int64_t,NDim>& size)
  : _size(size)
  , _stride({0})
  , _is_scalar(true)
  , _data(std::shared_ptr<T>(new T[1]{scalar}, std::default_delete<T[]>()))
  , _ownsdata(true)
{
  // Effectively this only checks _size for all-positive.
  check_stride(_size, cstride(_size));
}

template <typename T, int NDim>
DataBufferNd<T,NDim>::DataBufferNd(const std::array<std::int64_t,NDim>& size)
  : _size(size)
  , _stride(cstride(size))
  , _is_scalar(false)
  , _data(std::shared_ptr<T>(new T[make_product(size)], std::default_delete<T[]>()))
  , _ownsdata(true)
{
458
  SimpleTimerEx tt(AdHocTimers::instance().ctor);
459
  check_stride(_size, _stride);
460
461
462
463
464
  // TODO-Worry: Leaving the buffer not initialized might turn an easily
  // reproducible bug into a heissenbug. But zeroing the buffers here
  // is too wasteful. If you must, then do a fill() just after constructing
  // a new buffer,
  //memset(_data.get(), 0, make_product(size) * sizeof(T));
465
466
467
468
469
470
471
472
473
474
475
}

template <typename T, int NDim>
DataBufferNd<T,NDim>::DataBufferNd(const std::array<std::int64_t,NDim>& size,
                                   const std::array<std::int64_t,NDim>& stride)
  : _size(size)
  , _stride(stride)
  , _is_scalar(false)
  , _data(std::shared_ptr<T>(new T[make_allocsize(size, stride)], std::default_delete<T[]>()))
  , _ownsdata(true)
{
476
  SimpleTimerEx tt(AdHocTimers::instance().ctor);
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
  check_stride(_size, _stride);
  memset(_data.get(), 0, make_product(size) * sizeof(T));
}

template <typename T, int NDim>
DataBufferNd<T,NDim>::DataBufferNd(const std::shared_ptr<T>& data, const std::array<std::int64_t,NDim>& size)
  : _size(size)
  , _stride(cstride(size))
  , _is_scalar(false)
  , _data(data)
  , _ownsdata(false)
{
  check_stride(_size, _stride);
  if (!data)
    throw OpenZGY::Errors::ZgyInternalError("DataBuffer storage cannot be null");
}

template <typename T, int NDim>
DataBufferNd<T,NDim>::DataBufferNd(const std::shared_ptr<T>& data, const std::array<std::int64_t,NDim>& size, const std::array<std::int64_t,NDim>& stride)
  : _size(size)
  , _stride(stride)
  , _is_scalar(false)
  , _data(data)
  , _ownsdata(false)
{
  check_stride(_size, _stride);
  if (!data)
    throw OpenZGY::Errors::ZgyInternalError("DataBuffer storage cannot be null");
}

/**
 * \brief Check that size and stride are valid.
 *
 * Check both size and stride, and throw an exception on clearly invalid
 * cases such as overlapping areas. Non-fatal situations such as a
 * layout that is not C-style contiguous are returned in a bitmask.
 *
 * TODO-Low: Do I need to support anything other than C-style contiguous?
 * See discussion in the documentation for the class itself.
 * If no, this function can be written much simpler. If yes, the result
 * of this check should probably be made available to the caller. And/or
 * I might write a method make_contiguous_if_not_already().
 */
template <typename T, int NDim>
std::uint32_t
DataBufferNd<T,NDim>::check_stride(const ndsize_t& size, const ndsize_t& stride)
{
  std::uint32_t flags{ C_ORDERING | CONTIGUOUS | POSITIVE_STRIDE };

  // Order the size- and stride fastest changing first, i.e. by
  // descending stride. All sizes must be positive. Completely ignore
  // dimensions that have size 1 because the corresponding stride is
  // then irrelevant. Thpse entries will be missing from v, as the
  // rules do not apply to them.
  std::vector<std::pair<std::int64_t,std::int64_t>> v;
  for (int ii=0; ii<NDim; ++ii) {
    if (size[ii] == 0)
      throw OpenZGY::Errors::ZgyInternalError("Size cannot be 0");
    else if (size[ii] < 0)
      throw OpenZGY::Errors::ZgyInternalError("Size cannot be negative");
    else if (size[ii] == 1)
      flags |= DEGENERATE;
    else if (stride[ii] == 0)
      throw OpenZGY::Errors::ZgyInternalError("Stride cannot be 0");
    else /* size[ii] > 1 */ {
      if (stride[ii] < 0)
        flags &= (~POSITIVE_STRIDE);
      v.push_back(std::make_pair(std::abs(stride[ii]), size[ii]));
    }
  }

  // The "standard" ordering is first index slowest, last index fastest.
  if (v.empty() && !std::is_sorted(v.begin(), v.end()))
    flags &= (~C_ORDERING);

  // To check the strides, need then in "standard" order.
  std::sort(v.begin(), v.end());
  std::pair<std::int64_t,std::int64_t> prev{1,1};
  for (const auto& it : v) {
    if (it.first < prev.first * prev.second)
      throw OpenZGY::Errors::ZgyInternalError("Invalid stride");
    else if (it.first > prev.first * prev.second)
      flags &= (~CONTIGUOUS);
    prev = it;
  }
  return flags;
}

template <typename T, int NDim>
std::uint32_t
DataBufferNd<T,NDim>::layout() const
{
  // Scalars are equivalent to a plain c-style array.
  if (isScalar())
    return C_ORDERING | CONTIGUOUS | POSITIVE_STRIDE;
  // TODO-Low: Could also store the flags in the instance.
  return check_stride(_size, _stride);
}

template <typename T, int NDim>
bool
DataBufferNd<T,NDim>::is_cstride() const
{
  return ((layout() & (C_ORDERING | CONTIGUOUS | POSITIVE_STRIDE)) ==
          (C_ORDERING | CONTIGUOUS | POSITIVE_STRIDE));
}

template <typename T, int NDim>
void
DataBufferNd<T,NDim>::check_cstride() const
{
  if (!is_cstride())
    throw OpenZGY::Errors::ZgyInternalError("the DataBuffer is not c-contiguous");
}

/**
 * Set up a default stride with the last dimension varying fastest
 * and with no holes in the data.
 */
template <typename T, int NDim>
typename DataBufferNd<T,NDim>::ndsize_t
DataBufferNd<T,NDim>::cstride(const ndsize_t& size)
{
  ndsize_t result;
  result[NDim-1] = 1;
  for (int dim = NDim-2; dim >= 0; --dim)
    result[dim] = result[dim+1] * size[dim+1];
  return result;
}

/**
 * Convert a raw pointer to a fixed-length array. There is an unsafe
 * assumption thay the array passed by the caller contains the correct
 * number of elements.
 */
template <typename T, int NDim>
typename DataBufferNd<T,NDim>::ndsize_t
DataBufferNd<T,NDim>::make_array(const std::int64_t *in)
{
  ndsize_t result;
  for (int ii=0; ii<NDim; ++ii)
    result[ii] = in[ii];
  return result;
}

template <typename T, int NDim>
std::int64_t
DataBufferNd<T,NDim>::make_product(const ndsize_t& size)
{
  std::int64_t result = 1;
  for (int dim = 0; dim < NDim; ++dim)
    result *= size[dim];
  return result;
}

/**
 * Return the size in samples (not in bytes) of the data buffer,
 * including any padding added by adjusting the stride upwards,
 * but not padding added to the slowest varying index because
 * we cannot figure that out based on just size and stride.
 *
 * TODO-Low: this functions makes several assumptions about the stride.
 * The only deviation from the "obvious" stride is that padding can be
 * added by making the strides larger. I don't know yet whether this
 * feature will be needed at all. Or whether I will need support for
 * arbitrary size.
 *
 * Either way, the assumptions need to be explicitly checked.
 *
 * The current code assumes a simple layout possibly with padding
 * in i,j,k. The resulting size will not take i-padding into account
 * (it cannot be deduced from just size and stride) but it does
 * take j- and k-padding into account. So the result will be
 * i_size * padded_j_size * padded_k_size. Which might be slightly
 * more (padded_k_size-k_size) than actually needed. But it might
 * be less surprising to the caller.
 */
template <typename T, int NDim>
std::int64_t
DataBufferNd<T,NDim>::make_allocsize(const ndsize_t& size, const ndsize_t& stride)
{
  std::int64_t slowest = 0;
  for (int dim = 0; dim < NDim; ++dim)
    if (stride[dim] > stride[slowest])
      slowest = dim;
  return size[slowest] * stride[slowest];
}

template <typename T, int NDim>
std::string
DataBufferNd<T,NDim>::toString() const
{
  std::stringstream ss;
  ss << typeid(T).name() << " size (";
  for (int dim = 0; dim < NDim; ++dim)
    ss << _size[dim] << (dim == NDim-1 ? ")" : ", ");
  ss << " stride (";
  for (int dim = 0; dim < NDim; ++dim)
    ss << _stride[dim] << (dim == NDim-1 ? ")" : ", ");
  if (isScalar())
    ss << " constant " << (double)scalarValue();
  else
    ss << " ptr 0x" << std::hex << (intptr_t)voidData().get() << std::dec;
  return ss.str();
}

template <typename T, int NDim>
bool
DataBufferNd<T,NDim>::contiguous() const
{
  return isScalar() ? true : allocsize() == totalsize();
}

template <typename T, int NDim>
std::int64_t
DataBufferNd<T,NDim>::allocsize() const
{
  return isScalar() ? 1 : make_allocsize(_size, _stride);
}

template <typename T, int NDim>
std::int64_t
DataBufferNd<T,NDim>::totalsize() const
{
  return make_product(_size);
}

template <typename T, int NDim>
std::int64_t
DataBufferNd<T,NDim>::itemsize() const
{
  return sizeof(T);
}

/**
 * Return the size of the buffer, as a pointer to the first dimension.
 * This is unsafe since the compiler cannot check the length. This
 * function is meant for internal use in this file. Used by CopySubset
 * and CopySize. TODO-Low consider making those friends so that this method
 * can be declared private or inlined.
 */
template <typename T, int NDim>
const std::int64_t*
DataBufferNd<T,NDim>::sizeptr() const
{
  return _size.data();
}

/**
 * Return the layout of the buffer, assuming it is 3d. See sizeptr().
 */
template <typename T, int NDim>
const std::int64_t*
DataBufferNd<T,NDim>::strideptr() const
{
  return _stride.data();
}

/**
 * Return the size of the buffer, assuming it is 3d. Which it almost always is.
 * As a code smell this suggests that DataBuffer maybe shouldn't have been
 * templated on NDim in the first place. Or have NDim as a regular varable.
 * If fewer than 3 dimensions exist then the first index or indices are
 * treated as size=1, stride=0. If more then 3 dimensions only the last 3
 * will be returned. Which is probably not very useful.
 */
template <typename T, int NDim>
std::array<std::int64_t,3>
DataBufferNd<T,NDim>::size3d() const
{
  std::array<std::int64_t,3> result;
  result[2] = (NDim >= 1) ? _size[NDim-1] : 0;
  result[1] = (NDim >= 2) ? _size[NDim-2] : 1;
  result[0] = (NDim >= 3) ? _size[NDim-3] : 1;
  return result;
}

/**
 * Return the layout of the buffer, assuming it is 3d. See size3d().
 */
template <typename T, int NDim>
std::array<std::int64_t,3>
DataBufferNd<T,NDim>::stride3d() const
{
  std::array<std::int64_t,3> result;
  result[2] = (NDim >= 1) ? _stride[NDim-1] : 0;
  result[1] = (NDim >= 2) ? _stride[NDim-2] : 0;
  result[0] = (NDim >= 3) ? _stride[NDim-3] : 0;
  return result;
}

/**
 * \brief Set all samples to the same value.
 *
 * The value is passed as a float and will be cast to the correct type.
 * There is NO check for overflow.
 *
 * Some optimizing is attempted because gprof claims that (a) calling
 * this function amounts to 50% of the read overhead, and (b) memset,
 * if possible, would reduce that time significantly. I suspect that
 * both of those numbers are false. An ad-hoc Timer reports 10% in both
 * cases. But just in case I am wrong I have added the optimization.
 *
 * Using memset is possible if the value to be set is zero (the common
 * case for float cubes) and if the sample size is 1 (all int8 cubes).
 * Int16 cubes are less likely to benefit.
 */
template <typename T, int NDim>
void
DataBufferNd<T,NDim>::fill(double value)
{
788
  SimpleTimerEx tt(AdHocTimers::instance().fill);
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
  const T val = static_cast<T>(value);
  if (isScalar())
    _data.get()[0] = val;
  else {
    if (true && value == 0)
      memset(_data.get(), 0, allocsize() * sizeof(T));
    else if (false && sizeof(T) == 1)
      memset(_data.get(), static_cast<int>(val), allocsize() * sizeof(T));
    else
      std::fill(_data.get(), _data.get() + allocsize(), val);
  }
}

/**
 * Make the buffer empty. Size and stride are unchanged but the data will be
 * nullptr. This will cause an access violation if trying to read or write
 * data. The reason this is useful is that _data might point to something
 * that is not reference counted by our smart pointer. The pointer's deleter
 * is then a no-op. clear() should be called if we can no longer trust
 * that the pointer is valid. A reproducible null pointer exception is way
 * better than random crashes.
 */
template <typename T, int NDim>
void
DataBufferNd<T,NDim>::clear()
{
  _data.reset();
}

818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
/**
 * \brief Find lowest and highest sample excluding +/- Inf and NaN.
 *
 * \details The OpemMP directive is commented out because it might
 * unfortunately turn out to run slower. The inner loop is very
 * inexpensive which means the code is sensitive to how much overhead
 * OpenMP will add. I have made real-world measurements but only in
 * one hardware / OS / compiler environment and only for float input.
 * I observed that 1 million elements was the point where
 * multi-threading barely started to have a benefit. In a typical
 * scenario I don't expect to see more then 5 times that number. That
 * is a typical size of one brick-column. With 5 million floats the
 * task ran 3.5 times faster using 16 threads. With a less capable
 * OpenMP installation it might easily have run slower. It could
 * potentially run a *lot* slower. In my test, range() accounted for
 * 5% of total elapsed time. So I am, weighing a potential 5%
 * performance improvement against the risk of slowing down. I am
 * chicken. I hope I can instead parallelize this at a higher level.
 * Making this issue moot.
 */
838
839
840
841
template <typename T, int NDim>
std::pair<double, double>
DataBufferNd<T,NDim>::range() const
{
842
  SimpleTimerEx tt(AdHocTimers::instance().range);
843
844
  T min = std::numeric_limits<T>::max();
  T max = std::numeric_limits<T>::lowest();
845
  if (isScalar()) {
846
847
    if (IsFiniteT(scalarValue()))
      min = max = scalarValue();
848
849
850
  }
  else if (contiguous()) {
    // Short cut: No holes. Treat as 1d.
851
852
853
854
855
    const T * const ptr = data();
    const std::int64_t totalsize = allocsize();
//#pragma omp parallel for reduction(min: min) reduction(max: max) if(totalsize >= 1*1024*1024)
    for (std::int64_t ii=0; ii<totalsize; ++ii) {
      const T value = ptr[ii];
856
      if (!IsFiniteT(value)) continue;
857
858
      if (min > value) min = value;
      if (max < value) max = value;
859
860
861
862
    }
  }
  else {
    const T * const ptr = data();
863
//#pragma omp parallel for -- NOT TESTED YET
864
865
866
867
868
869
870
871
872
873
874
875
    for (int ii=0; ii<_size[0]; ++ii) {
      for (int jj=0; jj<_size[1]; ++jj) {
        for (int kk=0; kk<_size[2]; ++kk) {
          const T value = ptr[ii*_stride[0]+jj*_stride[1]+kk*_stride[2]];
          if (!IsFiniteT(value)) continue;
          else if (min > max) min = max = value;
          else if (min > value) min = value;
          else if (max < value) max = value;
        }
      }
    }
  }
876
  return std::make_pair(static_cast<double>(min), static_cast<double>(max));
877
878
879
880
881
882
}

template <typename T, int NDim>
bool
DataBufferNd<T,NDim>::isAllSame(const std::int64_t *used_in) const
{
883
  SimpleTimerEx tt(AdHocTimers::instance().allsame);
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
  // TODO-Low yet another place I wish I hadn't templeted on NDim.
  std::array<std::int64_t,NDim> used;
  for (int dim=0; dim<NDim; ++dim)
    used[dim] = used_in[dim];

  if (used[0]<=0 || used[1] <= 0 || used[2] <= 0 || isScalar()) {
    // Has zero or one sample, so clearly all samples are the same.
    return true;
  }
  else if (used[0] > _size[0] || used[1] > _size[1] || used[2] > _size[2]) {
    // Caller made an error. Should maybe assert.
    return false;
  }
  else if (used == _size && contiguous()) {
    // Short cut: No holes and no unused area. Treat as 1d.
    const T *ptr = data();
    const T *end = ptr + allocsize();
    const T first = *ptr++;
    if (IsNanT(first)) {
      for (; ptr < end; ++ptr)
        if (!IsNanT(*ptr))
          return false;
    }
    else {
      for (; ptr < end; ++ptr)
        if (*ptr != first)
          return false;
    }
  }
  else {
    const T *ptr = data();
    const T first = *ptr;
    if (IsNanT(first)) {
      for (int ii=0; ii<used[0]; ++ii)
        for (int jj=0; jj<used[1]; ++jj)
          for (int kk=0; kk<used[2]; ++kk)
            if (!IsNanT(ptr[ii*_stride[0]+jj*_stride[1]+kk*_stride[2]]))
              return false;
    }
    else{
      for (int ii=0; ii<used[0]; ++ii)
        for (int jj=0; jj<used[1]; ++jj)
          for (int kk=0; kk<used[2]; ++kk)
            if (ptr[ii*_stride[0]+jj*_stride[1]+kk*_stride[2]] != first)
              return false;
    }
  }
  return true;
}

template <typename T, int NDim>
void
DataBufferNd<T,NDim>::copySubset(const ndsize_t& srcorig, const self_type& src,
                                 const ndsize_t& dstorig,       self_type& dst,
                                 const ndsize_t& cpyorig, const ndsize_t& cpysize)
{
  if (src.isScalar())
    CopyScalar(NDim,
               srcorig.data(), src._size.data(), src.scalarValue(),
               dstorig.data(), dst._size.data(), dst._stride.data(), dst.data(),
               cpyorig.data(), cpysize.data());
  else
    CopySubset(NDim,
               srcorig.data(), src._size.data(), src._stride.data(), src.data(),
               dstorig.data(), dst._size.data(), dst._stride.data(), dst.data(),
               cpyorig.data(), cpysize.data());
}

template <typename T, int NDim>
void
DataBufferNd<T,NDim>::copySubset(const ndsize_t& srcorig, const self_type& src,
                                 const ndsize_t& dstorig,       self_type& dst)
{
  if (src.isScalar())
    CopyScalar(NDim,
               srcorig.data(), src._size.data(), src.scalarValue(),
               dstorig.data(), dst._size.data(), dst._stride.data(), dst.data(),
               nullptr, nullptr);
  else
    CopySubset(NDim,
               srcorig.data(), src._size.data(), src._stride.data(), src.data(),
               dstorig.data(), dst._size.data(), dst._stride.data(), dst.data(),
               nullptr, nullptr);
}

template <typename T, int NDim>
void
DataBufferNd<T,NDim>::copyFrom(const DataBuffer* src,
                               const std::int64_t *srcorig,
                               const std::int64_t *dstorig,
                               const std::int64_t *cpyorig,
                               const std::int64_t *cpysize)
{
977
978
979
  SimpleTimerEx tt(src && src->isScalar() ?
                   AdHocTimers::instance().copyscalar :
                   AdHocTimers::instance().copysubset);
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000

  const self_type* source = dynamic_cast<const self_type*>(src);
  if (!source) {
    if (src)
      throw OpenZGY::Errors::ZgyInternalError("Attempted copy from wrong buffer type.");
    else
      throw OpenZGY::Errors::ZgyInternalError("Attempted copy from null buffer.");
  }
  if (cpyorig && cpysize)
    copySubset(make_array(srcorig), *source,
               make_array(dstorig), *this,
               make_array(cpyorig), make_array(cpysize));
  else
    copySubset(make_array(srcorig), *source,
               make_array(dstorig), *this);
}

/**
 * Copy the input DataBuffer into a newly allocated buffer of float.
 * The input is assumed to have been read from a ZGY file,
 * so this function will also apply the storage to float transform.
For faster browsing, not all history is shown. View entire blame