DataContainer.h 32.6 KB
Newer Older
Tobias Lasser's avatar
Tobias Lasser committed
1
2
#pragma once

3
#include "elsaDefines.h"
Tobias Lasser's avatar
Tobias Lasser committed
4
5
#include "DataDescriptor.h"
#include "DataHandler.h"
Jens Petit's avatar
Jens Petit committed
6
7
#include "DataHandlerCPU.h"
#include "DataHandlerMapCPU.h"
8
#include "DataContainerIterator.h"
9
#include "Error.h"
Jens Petit's avatar
Jens Petit committed
10
#include "Expression.h"
11
#include "FormatConfig.h"
12
#include "TypeCasts.hpp"
Tobias Lasser's avatar
Tobias Lasser committed
13

14
15
16
17
18
#ifdef ELSA_CUDA_VECTOR
#include "DataHandlerGPU.h"
#include "DataHandlerMapGPU.h"
#endif

Tobias Lasser's avatar
Tobias Lasser committed
19
#include <memory>
20
#include <type_traits>
Tobias Lasser's avatar
Tobias Lasser committed
21
22
23
24

namespace elsa
{
    /**
25
     * @brief class representing and storing a linearized n-dimensional signal
Tobias Lasser's avatar
Tobias Lasser committed
26
27
28
29
     *
     * This class provides a container for a signal that is stored in memory. This signal can
     * be n-dimensional, and will be stored in memory in a linearized fashion. The information
     * on how this linearization is performed is provided by an associated DataDescriptor.
30
31
32
33
34
35
36
37
38
39
     *
     * @tparam data_t data type that is stored in the DataContainer, defaulting to real_t.
     *
     * @author
     * - Matthias Wieczorek - initial code
     * - Tobias Lasser - rewrite, modularization, modernization
     * - David Frank - added DataHandler concept, iterators
     * - Nikola Dinev - add block support
     * - Jens Petit - expression templates
     * - Jonas Jelten - various enhancements, fft, complex handling, pretty formatting
Tobias Lasser's avatar
Tobias Lasser committed
40
     */
Jens Petit's avatar
Jens Petit committed
41
    template <typename data_t>
Jens Petit's avatar
Jens Petit committed
42
43
    class DataContainer
    {
Tobias Lasser's avatar
Tobias Lasser committed
44
    public:
45
46
47
        /// Scalar alias
        using Scalar = data_t;

Tobias Lasser's avatar
Tobias Lasser committed
48
49
50
51
        /// delete default constructor (without metadata there can be no valid container)
        DataContainer() = delete;

        /**
Jonas Jelten's avatar
Jonas Jelten committed
52
53
         * @brief Constructor for empty DataContainer, no initialisation is performed,
         *        but the underlying space is allocated.
Tobias Lasser's avatar
Tobias Lasser committed
54
         *
55
56
         * @param[in] dataDescriptor containing the associated metadata
         * @param[in] handlerType the data handler (default: CPU)
Tobias Lasser's avatar
Tobias Lasser committed
57
58
         */
        explicit DataContainer(const DataDescriptor& dataDescriptor,
59
                               DataHandlerType handlerType = defaultHandlerType);
Jens Petit's avatar
Jens Petit committed
60

Tobias Lasser's avatar
Tobias Lasser committed
61
        /**
62
         * @brief Constructor for DataContainer, initializing it with a DataVector
Jens Petit's avatar
Jens Petit committed
63
         *
64
65
66
         * @param[in] dataDescriptor containing the associated metadata
         * @param[in] data vector containing the initialization data
         * @param[in] handlerType the data handler (default: CPU)
Tobias Lasser's avatar
Tobias Lasser committed
67
         */
Jens Petit's avatar
Jens Petit committed
68
69
        DataContainer(const DataDescriptor& dataDescriptor,
                      const Eigen::Matrix<data_t, Eigen::Dynamic, 1>& data,
70
                      DataHandlerType handlerType = defaultHandlerType);
Tobias Lasser's avatar
Tobias Lasser committed
71
72

        /**
73
         * @brief Copy constructor for DataContainer
Tobias Lasser's avatar
Tobias Lasser committed
74
         *
75
         * @param[in] other DataContainer to copy
Tobias Lasser's avatar
Tobias Lasser committed
76
77
78
79
         */
        DataContainer(const DataContainer<data_t>& other);

        /**
80
         * @brief copy assignment for DataContainer
Tobias Lasser's avatar
Tobias Lasser committed
81
         *
82
         * @param[in] other DataContainer to copy
83
84
85
86
         *
         * Note that a copy assignment with a DataContainer on a different device (CPU vs GPU) will
         * result in an "infectious" copy which means that afterwards the current container will use
         * the same device as "other".
Tobias Lasser's avatar
Tobias Lasser committed
87
88
89
90
         */
        DataContainer<data_t>& operator=(const DataContainer<data_t>& other);

        /**
91
         * @brief Move constructor for DataContainer
Tobias Lasser's avatar
Tobias Lasser committed
92
         *
93
         * @param[in] other DataContainer to move from
94
95
96
97
         *
         * The moved-from objects remains in a valid state. However, as preconditions are not
         * fulfilled for any member functions, the object should not be used. After move- or copy-
         * assignment, this is possible again.
Tobias Lasser's avatar
Tobias Lasser committed
98
         */
99
        DataContainer(DataContainer<data_t>&& other) noexcept;
Tobias Lasser's avatar
Tobias Lasser committed
100
101

        /**
102
         * @brief Move assignment for DataContainer
Tobias Lasser's avatar
Tobias Lasser committed
103
         *
104
         * @param[in] other DataContainer to move from
105
106
107
108
         *
         * The moved-from objects remains in a valid state. However, as preconditions are not
         * fulfilled for any member functions, the object should not be used. After move- or copy-
         * assignment, this is possible again.
109
110
111
112
         *
         * Note that a copy assignment with a DataContainer on a different device (CPU vs GPU) will
         * result in an "infectious" copy which means that afterwards the current container will use
         * the same device as "other".
Tobias Lasser's avatar
Tobias Lasser committed
113
         */
114
        DataContainer<data_t>& operator=(DataContainer<data_t>&& other);
Tobias Lasser's avatar
Tobias Lasser committed
115

Jens Petit's avatar
Jens Petit committed
116
        /**
117
         * @brief Expression evaluation assignment for DataContainer
Jens Petit's avatar
Jens Petit committed
118
         *
119
         * @param[in] source expression to evaluate
Jens Petit's avatar
Jens Petit committed
120
121
122
123
124
125
126
         *
         * This evaluates an expression template term into the underlying data member of
         * the DataHandler in use.
         */
        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
        DataContainer<data_t>& operator=(Source const& source)
        {
127
            if (auto handler = downcast_safe<DataHandlerCPU<data_t>>(_dataHandler.get())) {
128
129
                handler->accessData() = source.template eval<false>();
            } else if (auto handler =
130
                           downcast_safe<DataHandlerMapCPU<data_t>>(_dataHandler.get())) {
131
132
                handler->accessData() = source.template eval<false>();
#ifdef ELSA_CUDA_VECTOR
133
            } else if (auto handler = downcast_safe<DataHandlerGPU<data_t>>(_dataHandler.get())) {
134
135
                handler->accessData().eval(source.template eval<true>());
            } else if (auto handler =
136
                           downcast_safe<DataHandlerMapGPU<data_t>>(_dataHandler.get())) {
137
138
139
                handler->accessData().eval(source.template eval<true>());
#endif
            } else {
140
                throw LogicError("Unknown handler type");
141
            }
Jens Petit's avatar
Jens Petit committed
142
143
144
145
146

            return *this;
        }

        /**
147
         * @brief Expression constructor
Jens Petit's avatar
Jens Petit committed
148
         *
149
         * @param[in] source expression to evaluate
Jens Petit's avatar
Jens Petit committed
150
151
152
153
154
155
         *
         * It creates a new DataContainer out of an expression. For this the meta information which
         * is saved in the expression is used.
         */
        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
        DataContainer<data_t>(Source const& source)
156
            : DataContainer<data_t>(source.getDataMetaInfo().first, source.getDataMetaInfo().second)
Jens Petit's avatar
Jens Petit committed
157
        {
158
            this->operator=(source);
Jens Petit's avatar
Jens Petit committed
159
160
        }

Tobias Lasser's avatar
Tobias Lasser committed
161
162
163
        /// return the current DataDescriptor
        const DataDescriptor& getDataDescriptor() const;

Jens Petit's avatar
Jens Petit committed
164
165
        /// return the size of the stored data (i.e. the number of elements in the linearized
        /// signal)
Tobias Lasser's avatar
Tobias Lasser committed
166
167
168
169
170
171
172
173
174
        index_t getSize() const;

        /// return the index-th element of linearized signal (not bounds-checked!)
        data_t& operator[](index_t index);

        /// return the index-th element of the linearized signal as read-only (not bounds-checked!)
        const data_t& operator[](index_t index) const;

        /// return an element by n-dimensional coordinate (not bounds-checked!)
175
        data_t& operator()(const IndexVector_t& coordinate);
Tobias Lasser's avatar
Tobias Lasser committed
176
177

        /// return an element by n-dimensional coordinate as read-only (not bounds-checked!)
178
        const data_t& operator()(const IndexVector_t& coordinate) const;
Tobias Lasser's avatar
Tobias Lasser committed
179

180
        data_t at(const IndexVector_t& coordinate) const;
David Frank's avatar
David Frank committed
181

182
        /// return an element by its coordinates (not bounds-checked!)
Jens Petit's avatar
Jens Petit committed
183
184
185
186
187
188
189
        template <typename idx0_t, typename... idx_t,
                  typename = std::enable_if_t<
                      std::is_integral_v<idx0_t> && (... && std::is_integral_v<idx_t>)>>
        data_t& operator()(idx0_t idx0, idx_t... indices)
        {
            IndexVector_t coordinate(sizeof...(indices) + 1);
            ((coordinate << idx0), ..., indices);
190
191
192
            return operator()(coordinate);
        }

193
        /// return an element by its coordinates as read-only (not bounds-checked!)
Jens Petit's avatar
Jens Petit committed
194
195
196
197
198
199
200
        template <typename idx0_t, typename... idx_t,
                  typename = std::enable_if_t<
                      std::is_integral_v<idx0_t> && (... && std::is_integral_v<idx_t>)>>
        const data_t& operator()(idx0_t idx0, idx_t... indices) const
        {
            IndexVector_t coordinate(sizeof...(indices) + 1);
            ((coordinate << idx0), ..., indices);
201
202
            return operator()(coordinate);
        }
Tobias Lasser's avatar
Tobias Lasser committed
203
204
205
206

        /// return the dot product of this signal with the one from container other
        data_t dot(const DataContainer<data_t>& other) const;

Jens Petit's avatar
Jens Petit committed
207
208
209
210
        /// return the dot product of this signal with the one from an expression
        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
        data_t dot(const Source& source) const
        {
211
            if (auto handler = downcast_safe<DataHandlerCPU<data_t>>(_dataHandler.get())) {
212
213
                return (*this * source).template eval<false>().sum();
            } else if (auto handler =
214
                           downcast_safe<DataHandlerMapCPU<data_t>>(_dataHandler.get())) {
215
216
                return (*this * source).template eval<false>().sum();
#ifdef ELSA_CUDA_VECTOR
217
            } else if (auto handler = downcast_safe<DataHandlerGPU<data_t>>(_dataHandler.get())) {
218
219
220
                DataContainer temp = (*this * source);
                return temp.sum();
            } else if (auto handler =
221
                           downcast_safe<DataHandlerMapGPU<data_t>>(_dataHandler.get())) {
222
223
224
225
                DataContainer temp = (*this * source);
                return temp.sum();
#endif
            } else {
226
                throw LogicError("Unknown handler type");
227
            }
Jens Petit's avatar
Jens Petit committed
228
229
        }

Tobias Lasser's avatar
Tobias Lasser committed
230
        /// return the squared l2 norm of this signal (dot product with itself)
231
        GetFloatingPointType_t<data_t> squaredL2Norm() const;
Tobias Lasser's avatar
Tobias Lasser committed
232

David Frank's avatar
David Frank committed
233
234
235
        /// return the l2 norm of this signal (square root of dot product with itself)
        GetFloatingPointType_t<data_t> l2Norm() const;

Andi Braimllari's avatar
Andi Braimllari committed
236
237
238
        /// return the l0 pseudo-norm of this signal (number of non-zero values)
        index_t l0PseudoNorm() const;

Tobias Lasser's avatar
Tobias Lasser committed
239
        /// return the l1 norm of this signal (sum of absolute values)
240
        GetFloatingPointType_t<data_t> l1Norm() const;
Tobias Lasser's avatar
Tobias Lasser committed
241
242

        /// return the linf norm of this signal (maximum of absolute values)
243
        GetFloatingPointType_t<data_t> lInfNorm() const;
Tobias Lasser's avatar
Tobias Lasser committed
244
245
246
247

        /// return the sum of all elements of this signal
        data_t sum() const;

248
249
250
251
252
253
        /// return the min of all elements of this signal
        data_t minElement() const;

        /// return the max of all elements of this signal
        data_t maxElement() const;

Jonas Jelten's avatar
Jonas Jelten committed
254
        /// convert to the fourier transformed signal
255
        void fft(FFTNorm norm) const;
Jonas Jelten's avatar
Jonas Jelten committed
256
257

        /// convert to the inverse fourier transformed signal
258
        void ifft(FFTNorm norm) const;
Jonas Jelten's avatar
Jonas Jelten committed
259

260
261
262
263
264
265
266
267
268
269
        /// if the datacontainer is already complex, return itself.
        template <typename _data_t = data_t>
        typename std::enable_if_t<isComplex<_data_t>, DataContainer<_data_t>> asComplex() const
        {
            return *this;
        }

        /// if the datacontainer is not complex,
        /// return a copy and fill in 0 as imaginary values
        template <typename _data_t = data_t>
270
        typename std::enable_if_t<not isComplex<_data_t>, DataContainer<complex<_data_t>>>
271
272
            asComplex() const
        {
273
            DataContainer<complex<data_t>> ret{
274
275
276
277
278
279
                *this->_dataDescriptor,
                this->_dataHandlerType,
            };

            // extend with complex zero value
            for (index_t idx = 0; idx < this->getSize(); ++idx) {
280
                ret[idx] = complex<data_t>{(*this)[idx], 0};
281
282
283
284
285
            }

            return ret;
        }

Tobias Lasser's avatar
Tobias Lasser committed
286
287
288
        /// compute in-place element-wise addition of another container
        DataContainer<data_t>& operator+=(const DataContainer<data_t>& dc);

Jens Petit's avatar
Jens Petit committed
289
290
291
292
293
294
295
296
        /// compute in-place element-wise addition with another expression
        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
        DataContainer<data_t>& operator+=(Source const& source)
        {
            *this = *this + source;
            return *this;
        }

Tobias Lasser's avatar
Tobias Lasser committed
297
298
299
        /// compute in-place element-wise subtraction of another container
        DataContainer<data_t>& operator-=(const DataContainer<data_t>& dc);

Jens Petit's avatar
Jens Petit committed
300
301
302
303
304
305
306
307
        /// compute in-place element-wise subtraction with another expression
        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
        DataContainer<data_t>& operator-=(Source const& source)
        {
            *this = *this - source;
            return *this;
        }

Tobias Lasser's avatar
Tobias Lasser committed
308
309
310
        /// compute in-place element-wise multiplication with another container
        DataContainer<data_t>& operator*=(const DataContainer<data_t>& dc);

Jens Petit's avatar
Jens Petit committed
311
312
313
314
315
316
317
318
        /// compute in-place element-wise multiplication with another expression
        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
        DataContainer<data_t>& operator*=(Source const& source)
        {
            *this = *this * source;
            return *this;
        }

Tobias Lasser's avatar
Tobias Lasser committed
319
320
321
        /// compute in-place element-wise division by another container
        DataContainer<data_t>& operator/=(const DataContainer<data_t>& dc);

Jens Petit's avatar
Jens Petit committed
322
323
324
325
326
327
328
329
        /// compute in-place element-wise division with another expression
        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
        DataContainer<data_t>& operator/=(Source const& source)
        {
            *this = *this / source;
            return *this;
        }

Tobias Lasser's avatar
Tobias Lasser committed
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
        /// compute in-place addition of a scalar
        DataContainer<data_t>& operator+=(data_t scalar);

        /// compute in-place subtraction of a scalar
        DataContainer<data_t>& operator-=(data_t scalar);

        /// compute in-place multiplication with a scalar
        DataContainer<data_t>& operator*=(data_t scalar);

        /// compute in-place division by a scalar
        DataContainer<data_t>& operator/=(data_t scalar);

        /// assign a scalar to the DataContainer
        DataContainer<data_t>& operator=(data_t scalar);

        /// comparison with another DataContainer
        bool operator==(const DataContainer<data_t>& other) const;

        /// comparison with another DataContainer
        bool operator!=(const DataContainer<data_t>& other) const;

351
352
353
354
        /// returns a reference to the i-th block, wrapped in a DataContainer
        DataContainer<data_t> getBlock(index_t i);

        /// returns a const reference to the i-th block, wrapped in a DataContainer
355
        const DataContainer<data_t> getBlock(index_t i) const;
356
357
358
359
360

        /// return a view of this DataContainer with a different descriptor
        DataContainer<data_t> viewAs(const DataDescriptor& dataDescriptor);

        /// return a const view of this DataContainer with a different descriptor
361
        const DataContainer<data_t> viewAs(const DataDescriptor& dataDescriptor) const;
Tobias Lasser's avatar
Tobias Lasser committed
362

363
364
365
366
367
368
369
370
371
372
        /// @brief Slice the container in the last dimension
        ///
        /// Access a portion of the container via a slice. The slicing always is in the last
        /// dimension. So for a 3D volume, the slice would be an sliced in the z direction and would
        /// be a part of the x-y plane.
        ///
        /// A slice is always the same dimension as the original DataContainer, but with a thickness
        /// of 1 in the last dimension (i.e. the coefficient of the last dimension is 1)
        const DataContainer<data_t> slice(index_t i) const;

373
374
375
376
        /// @brief Slice the container in the last dimension, non-const overload
        ///
        /// @overload
        /// @see slice(index_t) const
377
378
        DataContainer<data_t> slice(index_t i);

379
380
381
382
383
384
385
        /// iterator for DataContainer (random access and continuous)
        using iterator = DataContainerIterator<DataContainer<data_t>>;

        /// const iterator for DataContainer (random access and continuous)
        using const_iterator = ConstDataContainerIterator<DataContainer<data_t>>;

        /// alias for reverse iterator
Jens Petit's avatar
Jens Petit committed
386
        using reverse_iterator = std::reverse_iterator<iterator>;
387
        /// alias for const reverse iterator
Jens Petit's avatar
Jens Petit committed
388
        using const_reverse_iterator = std::reverse_iterator<const_iterator>;
389
390
391
392
393
394
395
396
397
398
399
400
401

        /// returns iterator to the first element of the container
        iterator begin();

        /// returns const iterator to the first element of the container (cannot mutate data)
        const_iterator begin() const;

        /// returns const iterator to the first element of the container (cannot mutate data)
        const_iterator cbegin() const;

        /// returns iterator to one past the last element of the container
        iterator end();

Jens Petit's avatar
Jens Petit committed
402
403
        /// returns const iterator to one past the last element of the container (cannot mutate
        /// data)
404
405
        const_iterator end() const;

Jens Petit's avatar
Jens Petit committed
406
407
        /// returns const iterator to one past the last element of the container (cannot mutate
        /// data)
408
409
410
411
412
        const_iterator cend() const;

        /// returns reversed iterator to the last element of the container
        reverse_iterator rbegin();

Jens Petit's avatar
Jens Petit committed
413
414
        /// returns const reversed iterator to the last element of the container (cannot mutate
        /// data)
415
416
        const_reverse_iterator rbegin() const;

Jens Petit's avatar
Jens Petit committed
417
418
        /// returns const reversed iterator to the last element of the container (cannot mutate
        /// data)
419
420
421
422
423
        const_reverse_iterator crbegin() const;

        /// returns reversed iterator to one past the first element of container
        reverse_iterator rend();

Jens Petit's avatar
Jens Petit committed
424
425
        /// returns const reversed iterator to one past the first element of container (cannot
        /// mutate data)
426
427
        const_reverse_iterator rend() const;

Jens Petit's avatar
Jens Petit committed
428
429
        /// returns const reversed iterator to one past the first element of container (cannot
        /// mutate data)
430
431
432
        const_reverse_iterator crend() const;

        /// value_type of the DataContainer elements for iterators
Jens Petit's avatar
Jens Petit committed
433
        using value_type = data_t;
434
        /// pointer type of DataContainer elements for iterators
Jens Petit's avatar
Jens Petit committed
435
        using pointer = data_t*;
436
        /// const pointer type of DataContainer elements for iterators
Jens Petit's avatar
Jens Petit committed
437
        using const_pointer = const data_t*;
438
        /// reference type of DataContainer elements for iterators
Jens Petit's avatar
Jens Petit committed
439
        using reference = data_t&;
440
441
442
443
444
        /// const reference type of DataContainer elements for iterators
        using const_reference = const data_t&;
        /// difference type for iterators
        using difference_type = std::ptrdiff_t;

Jens Petit's avatar
Jens Petit committed
445
446
447
448
        /// returns the type of the DataHandler in use
        DataHandlerType getDataHandlerType() const;

        /// friend constexpr function to implement expression templates
449
        template <bool GPU, class Operand, std::enable_if_t<isDataContainer<Operand>, int>>
Jens Petit's avatar
Jens Petit committed
450
451
        friend constexpr auto evaluateOrReturn(Operand const& operand);

452
        /// write a pretty-formatted string representation to stream
453
        void format(std::ostream& os, format_config cfg = {}) const;
454

455
        /**
456
         * @brief Factory function which returns GPU based DataContainers
457
         *
458
         * @return the GPU based DataContainer
459
         *
460
461
         * Note that if this function is called on a container which is already GPU based, it
         * will throw an exception.
462
463
464
465
         */
        DataContainer loadToGPU();

        /**
466
         * @brief Factory function which returns CPU based DataContainers
467
         *
468
         * @return the CPU based DataContainer
469
470
471
472
473
474
         *
         * Note that if this function is called on a container which is already CPU based, it will
         * throw an exception.
         */
        DataContainer loadToCPU();

Tobias Lasser's avatar
Tobias Lasser committed
475
476
477
    private:
        /// the current DataDescriptor
        std::unique_ptr<DataDescriptor> _dataDescriptor;
Jens Petit's avatar
Jens Petit committed
478

Tobias Lasser's avatar
Tobias Lasser committed
479
        /// the current DataHandler
480
        std::unique_ptr<DataHandler<data_t>> _dataHandler;
Tobias Lasser's avatar
Tobias Lasser committed
481

Jens Petit's avatar
Jens Petit committed
482
483
484
        /// the current DataHandlerType
        DataHandlerType _dataHandlerType;

Jens Petit's avatar
Jens Petit committed
485
486
487
488
489
        /// factory method to create DataHandlers based on handlerType with perfect forwarding of
        /// constructor arguments
        template <typename... Args>
        std::unique_ptr<DataHandler<data_t>> createDataHandler(DataHandlerType handlerType,
                                                               Args&&... args);
Tobias Lasser's avatar
Tobias Lasser committed
490
491

        /// private constructor accepting a DataDescriptor and a DataHandler
Jens Petit's avatar
Jens Petit committed
492
        explicit DataContainer(const DataDescriptor& dataDescriptor,
Jens Petit's avatar
Jens Petit committed
493
                               std::unique_ptr<DataHandler<data_t>> dataHandler,
494
                               DataHandlerType dataType = defaultHandlerType);
495
496

        /**
497
         * @brief Helper function to indicate if a regular assignment or a clone should be performed
498
         *
499
         * @param[in] handlerType the member variable of the other container in
500
501
         * copy-/move-assignment
         *
502
         * @return true if a regular assignment of the pointed to DataHandlers should be done
503
504
505
506
507
508
         *
         * An assignment operation with a DataContainer which does not use the same device (CPU /
         * GPU) has to be handled differently. This helper function indicates if a regular
         * assignment should be performed or not.
         */
        bool canAssign(DataHandlerType handlerType);
Tobias Lasser's avatar
Tobias Lasser committed
509
510
    };

511
512
513
514
515
516
517
518
519
    /// pretty output formatting.
    /// for configurable output, use `DataContainerFormatter` directly.
    template <typename T>
    std::ostream& operator<<(std::ostream& os, const elsa::DataContainer<T>& dc)
    {
        dc.format(os);
        return os;
    }

andibraimllari's avatar
andibraimllari committed
520
521
522
523
    /// clip the container values outside of the interval, to the interval edges
    template <typename data_t>
    DataContainer<data_t> clip(DataContainer<data_t> dc, data_t min, data_t max);

524
525
526
527
528
    /// Concatenate two DataContainers to one (requires copying of both)
    template <typename data_t>
    DataContainer<data_t> concatenate(const DataContainer<data_t>& dc1,
                                      const DataContainer<data_t>& dc2);

529
530
531
532
533
534
535
536
537
538
539
540
    /// Perform the FFT shift operation to the provided signal. Refer to
    /// https://numpy.org/doc/stable/reference/generated/numpy.fft.fftshift.html for further
    /// details.
    template <typename data_t>
    DataContainer<data_t> fftShift2D(DataContainer<data_t> dc);

    /// Perform the IFFT shift operation to the provided signal. Refer to
    /// https://numpy.org/doc/stable/reference/generated/numpy.fft.ifftshift.html for further
    /// details.
    template <typename data_t>
    DataContainer<data_t> ifftShift2D(DataContainer<data_t> dc);

541
542
    /// User-defined template argument deduction guide for the expression based constructor
    template <typename Source>
543
    DataContainer(Source const& source) -> DataContainer<typename Source::data_t>;
544

545
546
547
548
549
550
551
552
    /// Collects callable lambdas for later dispatch
    template <typename... Ts>
    struct Callables : Ts... {
        using Ts::operator()...;
    };

    /// Class template deduction guide
    template <typename... Ts>
553
    Callables(Ts...) -> Callables<Ts...>;
554

Jens Petit's avatar
Jens Petit committed
555
556
557
    /// Multiplying two operands (including scalars)
    template <typename LHS, typename RHS, typename = std::enable_if_t<isBinaryOpOk<LHS, RHS>>>
    auto operator*(LHS const& lhs, RHS const& rhs)
Tobias Lasser's avatar
Tobias Lasser committed
558
    {
559
560
561
562
        auto multiplicationGPU = [](auto const& left, auto const& right, bool /**/) {
            return left * right;
        };

Jens Petit's avatar
Jens Petit committed
563
564
565
566
        if constexpr (isDcOrExpr<LHS> && isDcOrExpr<RHS>) {
            auto multiplication = [](auto const& left, auto const& right) {
                return (left.array() * right.array()).matrix();
            };
567
            return Expression{Callables{multiplication, multiplicationGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
568
569
570
571
        } else if constexpr (isArithmetic<LHS>) {
            auto multiplication = [](auto const& left, auto const& right) {
                return (left * right.array()).matrix();
            };
572
            return Expression{Callables{multiplication, multiplicationGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
573
574
575
576
        } else if constexpr (isArithmetic<RHS>) {
            auto multiplication = [](auto const& left, auto const& right) {
                return (left.array() * right).matrix();
            };
577
            return Expression{Callables{multiplication, multiplicationGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
578
579
        } else {
            auto multiplication = [](auto const& left, auto const& right) { return left * right; };
580
            return Expression{Callables{multiplication, multiplicationGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
581
        }
Tobias Lasser's avatar
Tobias Lasser committed
582
583
    }

Jens Petit's avatar
Jens Petit committed
584
585
586
    /// Adding two operands (including scalars)
    template <typename LHS, typename RHS, typename = std::enable_if_t<isBinaryOpOk<LHS, RHS>>>
    auto operator+(LHS const& lhs, RHS const& rhs)
Tobias Lasser's avatar
Tobias Lasser committed
587
    {
588
589
590
591
        auto additionGPU = [](auto const& left, auto const& right, bool /**/) {
            return left + right;
        };

Jens Petit's avatar
Jens Petit committed
592
593
        if constexpr (isDcOrExpr<LHS> && isDcOrExpr<RHS>) {
            auto addition = [](auto const& left, auto const& right) { return left + right; };
594
            return Expression{Callables{addition, additionGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
595
596
597
598
        } else if constexpr (isArithmetic<LHS>) {
            auto addition = [](auto const& left, auto const& right) {
                return (left + right.array()).matrix();
            };
599
            return Expression{Callables{addition, additionGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
600
601
602
603
        } else if constexpr (isArithmetic<RHS>) {
            auto addition = [](auto const& left, auto const& right) {
                return (left.array() + right).matrix();
            };
604
            return Expression{Callables{addition, additionGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
605
606
        } else {
            auto addition = [](auto const& left, auto const& right) { return left + right; };
607
            return Expression{Callables{addition, additionGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
608
        }
Tobias Lasser's avatar
Tobias Lasser committed
609
610
    }

Jens Petit's avatar
Jens Petit committed
611
612
613
    /// Subtracting two operands (including scalars)
    template <typename LHS, typename RHS, typename = std::enable_if_t<isBinaryOpOk<LHS, RHS>>>
    auto operator-(LHS const& lhs, RHS const& rhs)
Tobias Lasser's avatar
Tobias Lasser committed
614
    {
615
616
617
618
        auto subtractionGPU = [](auto const& left, auto const& right, bool /**/) {
            return left - right;
        };

Jens Petit's avatar
Jens Petit committed
619
620
        if constexpr (isDcOrExpr<LHS> && isDcOrExpr<RHS>) {
            auto subtraction = [](auto const& left, auto const& right) { return left - right; };
621
            return Expression{Callables{subtraction, subtractionGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
622
623
624
625
        } else if constexpr (isArithmetic<LHS>) {
            auto subtraction = [](auto const& left, auto const& right) {
                return (left - right.array()).matrix();
            };
626
            return Expression{Callables{subtraction, subtractionGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
627
628
629
630
        } else if constexpr (isArithmetic<RHS>) {
            auto subtraction = [](auto const& left, auto const& right) {
                return (left.array() - right).matrix();
            };
631
            return Expression{Callables{subtraction, subtractionGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
632
633
        } else {
            auto subtraction = [](auto const& left, auto const& right) { return left - right; };
634
            return Expression{Callables{subtraction, subtractionGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
635
        }
Tobias Lasser's avatar
Tobias Lasser committed
636
637
    }

Jens Petit's avatar
Jens Petit committed
638
639
640
    /// Dividing two operands (including scalars)
    template <typename LHS, typename RHS, typename = std::enable_if_t<isBinaryOpOk<LHS, RHS>>>
    auto operator/(LHS const& lhs, RHS const& rhs)
Tobias Lasser's avatar
Tobias Lasser committed
641
    {
642
643
644
645
        auto divisionGPU = [](auto const& left, auto const& right, bool /**/) {
            return left / right;
        };

Jens Petit's avatar
Jens Petit committed
646
647
648
649
        if constexpr (isDcOrExpr<LHS> && isDcOrExpr<RHS>) {
            auto division = [](auto const& left, auto const& right) {
                return (left.array() / right.array()).matrix();
            };
650
            return Expression{Callables{division, divisionGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
651
652
653
654
        } else if constexpr (isArithmetic<LHS>) {
            auto division = [](auto const& left, auto const& right) {
                return (left / right.array()).matrix();
            };
655
            return Expression{Callables{division, divisionGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
656
657
658
659
        } else if constexpr (isArithmetic<RHS>) {
            auto division = [](auto const& left, auto const& right) {
                return (left.array() / right).matrix();
            };
660
            return Expression{Callables{division, divisionGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
661
        } else {
662
663
            auto division = [](auto const& left, auto const& right) { return left / right; };
            return Expression{Callables{division, divisionGPU}, lhs, rhs};
Jens Petit's avatar
Jens Petit committed
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
    /// Element-wise maximum value operation between two operands
    template <typename LHS, typename RHS, typename = std::enable_if_t<isBinaryOpOk<LHS, RHS>>>
    auto cwiseMax(LHS const& lhs, RHS const& rhs)
    {
        constexpr bool isLHSComplex = isComplex<GetOperandDataType_t<LHS>>;
        constexpr bool isRHSComplex = isComplex<GetOperandDataType_t<RHS>>;

#ifdef ELSA_CUDA_VECTOR
        auto cwiseMaxGPU = [](auto const& lhs, auto const& rhs, bool) {
            return quickvec::cwiseMax(lhs, rhs);
        };
#endif
        auto cwiseMax = [] {
            if constexpr (isLHSComplex && isRHSComplex) {
                return [](auto const& left, auto const& right) {
                    return (left.array().abs().max(right.array().abs())).matrix();
                };
            } else if constexpr (isLHSComplex) {
                return [](auto const& left, auto const& right) {
                    return (left.array().abs().max(right.array())).matrix();
                };
            } else if constexpr (isRHSComplex) {
                return [](auto const& left, auto const& right) {
                    return (left.array().max(right.array().abs())).matrix();
                };
            } else {
                return [](auto const& left, auto const& right) {
                    return (left.array().max(right.array())).matrix();
                };
            }
        }();

#ifdef ELSA_CUDA_VECTOR
        return Expression{Callables{cwiseMax, cwiseMaxGPU}, lhs, rhs};
#else
        return Expression{cwiseMax, lhs, rhs};
#endif
    }

706
707
708
709
710
711
712
713
714
715
716
    /// Element-wise absolute value operation
    template <typename Operand, typename = std::enable_if_t<isDcOrExpr<Operand>>>
    auto cwiseAbs(Operand const& operand)
    {
        auto abs = [](auto const& operand) { return (operand.array().abs()).matrix(); };
#ifdef ELSA_CUDA_VECTOR
        auto absGPU = [](auto const& operand, bool) { return quickvec::cwiseAbs(operand); };
        return Expression{Callables{abs, absGPU}, operand};
#else
        return Expression{abs, operand};
#endif
Tobias Lasser's avatar
Tobias Lasser committed
717
718
    }

Jens Petit's avatar
Jens Petit committed
719
720
721
    /// Element-wise square operation
    template <typename Operand, typename = std::enable_if_t<isDcOrExpr<Operand>>>
    auto square(Operand const& operand)
Tobias Lasser's avatar
Tobias Lasser committed
722
    {
Jens Petit's avatar
Jens Petit committed
723
        auto square = [](auto const& operand) { return (operand.array().square()).matrix(); };
724
725
726
727
#ifdef ELSA_CUDA_VECTOR
        auto squareGPU = [](auto const& operand, bool /**/) { return quickvec::square(operand); };
        return Expression{Callables{square, squareGPU}, operand};
#else
Jens Petit's avatar
Jens Petit committed
728
        return Expression{square, operand};
729
#endif
Tobias Lasser's avatar
Tobias Lasser committed
730
731
    }

Jens Petit's avatar
Jens Petit committed
732
733
734
    /// Element-wise square-root operation
    template <typename Operand, typename = std::enable_if_t<isDcOrExpr<Operand>>>
    auto sqrt(Operand const& operand)
Tobias Lasser's avatar
Tobias Lasser committed
735
    {
Jens Petit's avatar
Jens Petit committed
736
        auto sqrt = [](auto const& operand) { return (operand.array().sqrt()).matrix(); };
737
738
739
740
#ifdef ELSA_CUDA_VECTOR
        auto sqrtGPU = [](auto const& operand, bool /**/) { return quickvec::sqrt(operand); };
        return Expression{Callables{sqrt, sqrtGPU}, operand};
#else
Jens Petit's avatar
Jens Petit committed
741
        return Expression{sqrt, operand};
742
#endif
Tobias Lasser's avatar
Tobias Lasser committed
743
744
    }

Jens Petit's avatar
Jens Petit committed
745
746
747
    /// Element-wise exponenation operation with euler base
    template <typename Operand, typename = std::enable_if_t<isDcOrExpr<Operand>>>
    auto exp(Operand const& operand)
Tobias Lasser's avatar
Tobias Lasser committed
748
    {
Jens Petit's avatar
Jens Petit committed
749
        auto exp = [](auto const& operand) { return (operand.array().exp()).matrix(); };
750
751
752
753
#ifdef ELSA_CUDA_VECTOR
        auto expGPU = [](auto const& operand, bool /**/) { return quickvec::exp(operand); };
        return Expression{Callables{exp, expGPU}, operand};
#else
Jens Petit's avatar
Jens Petit committed
754
        return Expression{exp, operand};
755
#endif
Tobias Lasser's avatar
Tobias Lasser committed
756
757
    }

Jens Petit's avatar
Jens Petit committed
758
759
760
    /// Element-wise natural logarithm
    template <typename Operand, typename = std::enable_if_t<isDcOrExpr<Operand>>>
    auto log(Operand const& operand)
Tobias Lasser's avatar
Tobias Lasser committed
761
    {
Jens Petit's avatar
Jens Petit committed
762
        auto log = [](auto const& operand) { return (operand.array().log()).matrix(); };
763
764
765
766
#ifdef ELSA_CUDA_VECTOR
        auto logGPU = [](auto const& operand, bool /**/) { return quickvec::log(operand); };
        return Expression{Callables{log, logGPU}, operand};
#else
Jens Petit's avatar
Jens Petit committed
767
        return Expression{log, operand};
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
#endif
    }

    /// Element-wise real parts of the Operand
    template <typename Operand, typename = std::enable_if_t<isDcOrExpr<Operand>>>
    auto real(Operand const& operand)
    {
        auto real = [](auto const& operand) { return (operand.array().real()).matrix(); };
#ifdef ELSA_CUDA_VECTOR
        auto realGPU = [](auto const& operand, bool) { return quickvec::real(operand); };
        return Expression{Callables{real, realGPU}, operand};
#else
        return Expression{real, operand};
#endif
    }

    /// Element-wise imaginary parts of the Operand
    template <typename Operand, typename = std::enable_if_t<isDcOrExpr<Operand>>>
    auto imag(Operand const& operand)
    {
        auto imag = [](auto const& operand) { return (operand.array().imag()).matrix(); };
#ifdef ELSA_CUDA_VECTOR
        auto imagGPU = [](auto const& operand, bool) { return quickvec::imag(operand); };
        return Expression{Callables{imag, imagGPU}, operand};
#else
        return Expression{imag, operand};
794
#endif
Tobias Lasser's avatar
Tobias Lasser committed
795
796
    }
} // namespace elsa