Commit 15d1fbb1 authored by Tobias Lasser's avatar Tobias Lasser

Add DataHandlerMapCPU, move copy-on-write functionality to DataHandler, add...

Add DataHandlerMapCPU, move copy-on-write functionality to DataHandler, add block functionality to DataContainer
parent bf27c029
Pipeline #186717 passed with stages
in 18 minutes and 42 seconds
......@@ -33,12 +33,12 @@ namespace elsa
index_t BlockDescriptor::getNumberOfBlocks() const { return _blockDescriptors.size(); }
const DataDescriptor& BlockDescriptor::getIthDescriptor(index_t i) const
const DataDescriptor& BlockDescriptor::getDescriptorOfBlock(index_t i) const
{
return *_blockDescriptors.at(i);
}
index_t BlockDescriptor::getIthBlockOffset(elsa::index_t i) const
index_t BlockDescriptor::getOffsetOfBlock(elsa::index_t i) const
{
if (i < 0 || i >= _blockOffsets.size())
throw std::invalid_argument("BlockDescriptor: index i is out of bounds");
......
......@@ -50,10 +50,10 @@ namespace elsa
index_t getNumberOfBlocks() const;
/// return the DataDescriptor of the i-th block
const DataDescriptor& getIthDescriptor(index_t i) const;
const DataDescriptor& getDescriptorOfBlock(index_t i) const;
/// return the offset to access the data of the i-th block
index_t getIthBlockOffset(index_t i) const;
index_t getOffsetOfBlock(index_t i) const;
protected:
/// vector of DataDescriptors describing the individual blocks
......
......@@ -16,6 +16,7 @@ set(MODULE_HEADERS
DataContainerIterator.h
DataHandler.h
DataHandlerCPU.h
DataHandlerMapCPU.h
LinearOperator.h)
# list all the code files of the module
......@@ -24,6 +25,7 @@ set(MODULE_SOURCES
BlockDescriptor.cpp
DataContainer.cpp
DataHandlerCPU.cpp
DataHandlerMapCPU.cpp
LinearOperator.cpp)
......
......@@ -27,10 +27,7 @@ namespace elsa
std::unique_ptr<Derived> clone() const { return std::unique_ptr<Derived>(cloneImpl()); }
/// comparison operators
bool operator==(const Derived& other) const
{
return typeid(*this) == typeid(other) && isEqual(other);
}
bool operator==(const Derived& other) const { return isEqual(other); }
bool operator!=(const Derived& other) const { return !(*this == other); }
......
#include "DataContainer.h"
#include "DataHandlerCPU.h"
#include "DataHandlerMapCPU.h"
#include "BlockDescriptor.h"
#include <stdexcept>
#include <utility>
......@@ -31,7 +33,7 @@ namespace elsa
template <typename data_t>
DataContainer<data_t>::DataContainer(const DataContainer<data_t>& other)
: _dataDescriptor{other._dataDescriptor->clone()}, _dataHandler{other._dataHandler}
: _dataDescriptor{other._dataDescriptor->clone()}, _dataHandler{other._dataHandler->clone()}
{
}
......@@ -40,7 +42,12 @@ namespace elsa
{
if (this != &other) {
_dataDescriptor = other._dataDescriptor->clone();
_dataHandler = other._dataHandler;
if (_dataHandler) {
*_dataHandler = *other._dataHandler;
} else {
_dataHandler = other._dataHandler->clone();
}
}
return *this;
......@@ -57,10 +64,14 @@ namespace elsa
}
template <typename data_t>
DataContainer<data_t>& DataContainer<data_t>::operator=(DataContainer<data_t>&& other) noexcept
DataContainer<data_t>& DataContainer<data_t>::operator=(DataContainer<data_t>&& other)
{
_dataDescriptor = std::move(other._dataDescriptor);
_dataHandler = std::move(other._dataHandler);
if (_dataHandler) {
*_dataHandler = std::move(*other._dataHandler);
} else {
_dataHandler = std::move(other._dataHandler);
}
// leave other in a valid state
other._dataDescriptor = nullptr;
......@@ -84,27 +95,26 @@ namespace elsa
template <typename data_t>
data_t& DataContainer<data_t>::operator[](index_t index)
{
detach();
return (*_dataHandler)[index];
}
template <typename data_t>
const data_t& DataContainer<data_t>::operator[](index_t index) const
{
return (*_dataHandler)[index];
return static_cast<const DataHandler<data_t>&>(*_dataHandler)[index];
}
template <typename data_t>
data_t& DataContainer<data_t>::operator()(IndexVector_t coordinate)
{
detach();
return (*_dataHandler)[_dataDescriptor->getIndexFromCoordinate(coordinate)];
}
template <typename data_t>
const data_t& DataContainer<data_t>::operator()(IndexVector_t coordinate) const
{
return (*_dataHandler)[_dataDescriptor->getIndexFromCoordinate(coordinate)];
return static_cast<const DataHandler<data_t>&>(
*_dataHandler)[_dataDescriptor->getIndexFromCoordinate(coordinate)];
}
template <typename data_t>
......@@ -164,7 +174,6 @@ namespace elsa
template <typename data_t>
DataContainer<data_t>& DataContainer<data_t>::operator+=(const DataContainer<data_t>& dc)
{
detach();
*_dataHandler += *dc._dataHandler;
return *this;
}
......@@ -172,7 +181,6 @@ namespace elsa
template <typename data_t>
DataContainer<data_t>& DataContainer<data_t>::operator-=(const DataContainer<data_t>& dc)
{
detach();
*_dataHandler -= *dc._dataHandler;
return *this;
}
......@@ -180,7 +188,6 @@ namespace elsa
template <typename data_t>
DataContainer<data_t>& DataContainer<data_t>::operator*=(const DataContainer<data_t>& dc)
{
detach();
*_dataHandler *= *dc._dataHandler;
return *this;
}
......@@ -188,7 +195,6 @@ namespace elsa
template <typename data_t>
DataContainer<data_t>& DataContainer<data_t>::operator/=(const DataContainer<data_t>& dc)
{
detach();
*_dataHandler /= *dc._dataHandler;
return *this;
}
......@@ -196,7 +202,6 @@ namespace elsa
template <typename data_t>
DataContainer<data_t>& DataContainer<data_t>::operator+=(data_t scalar)
{
detach();
*_dataHandler += scalar;
return *this;
}
......@@ -204,7 +209,6 @@ namespace elsa
template <typename data_t>
DataContainer<data_t>& DataContainer<data_t>::operator-=(data_t scalar)
{
detach();
*_dataHandler -= scalar;
return *this;
}
......@@ -212,7 +216,6 @@ namespace elsa
template <typename data_t>
DataContainer<data_t>& DataContainer<data_t>::operator*=(data_t scalar)
{
detach();
*_dataHandler *= scalar;
return *this;
}
......@@ -220,7 +223,6 @@ namespace elsa
template <typename data_t>
DataContainer<data_t>& DataContainer<data_t>::operator/=(data_t scalar)
{
detach();
*_dataHandler /= scalar;
return *this;
}
......@@ -228,7 +230,6 @@ namespace elsa
template <typename data_t>
DataContainer<data_t>& DataContainer<data_t>::operator=(data_t scalar)
{
detach();
*_dataHandler = scalar;
return *this;
}
......@@ -272,19 +273,67 @@ namespace elsa
}
template <typename data_t>
void DataContainer<data_t>::detach()
DataContainer<data_t> DataContainer<data_t>::getBlock(index_t i)
{
if (_dataHandler.use_count() != 1) {
#pragma omp barrier
#pragma omp single
_dataHandler = _dataHandler->clone();
}
const auto blockDesc = dynamic_cast<const BlockDescriptor*>(_dataDescriptor.get());
if (!blockDesc)
throw std::logic_error("DataContainer: cannot get block from not-blocked container");
if (i >= blockDesc->getNumberOfBlocks() || i < 0)
throw std::invalid_argument("DataContainer: block index out of bounds");
index_t startIndex = blockDesc->getOffsetOfBlock(i);
const auto& ithDesc = blockDesc->getDescriptorOfBlock(i);
index_t blockSize = ithDesc.getNumberOfCoefficients();
return DataContainer<data_t>{ithDesc, _dataHandler->getBlock(startIndex, blockSize)};
}
template <typename data_t>
const DataContainer<data_t> DataContainer<data_t>::getBlock(index_t i) const
{
const auto blockDesc = dynamic_cast<const BlockDescriptor*>(_dataDescriptor.get());
if (!blockDesc)
throw std::logic_error("DataContainer: cannot get block from not-blocked container");
if (i >= blockDesc->getNumberOfBlocks() || i < 0)
throw std::invalid_argument("DataContainer: block index out of bounds");
index_t startIndex = blockDesc->getOffsetOfBlock(i);
const auto& ithDesc = blockDesc->getDescriptorOfBlock(i);
index_t blockSize = ithDesc.getNumberOfCoefficients();
// getBlock() returns a pointer to non-const DH, but that's fine as it gets wrapped in a
// constant container
return DataContainer<data_t>{ithDesc, _dataHandler->getBlock(startIndex, blockSize)};
}
template <typename data_t>
DataContainer<data_t> DataContainer<data_t>::viewAs(const DataDescriptor& dataDescriptor)
{
if (dataDescriptor.getNumberOfCoefficients() != getSize())
throw std::invalid_argument("DataContainer: view must have same size as container");
return DataContainer<data_t>{dataDescriptor, _dataHandler->getBlock(0, getSize())};
}
template <typename data_t>
const DataContainer<data_t>
DataContainer<data_t>::viewAs(const DataDescriptor& dataDescriptor) const
{
if (dataDescriptor.getNumberOfCoefficients() != getSize())
throw std::invalid_argument("DataContainer: view must have same size as container");
// getBlock() returns a pointer to non-const DH, but that's fine as it gets wrapped in a
// constant container
return DataContainer<data_t>{dataDescriptor, _dataHandler->getBlock(0, getSize())};
}
template <typename data_t>
typename DataContainer<data_t>::iterator DataContainer<data_t>::begin()
{
detach();
return iterator(&(*this)[0]);
}
......@@ -303,7 +352,6 @@ namespace elsa
template <typename data_t>
typename DataContainer<data_t>::iterator DataContainer<data_t>::end()
{
detach();
return iterator(&(*this)[0] + getSize());
}
......@@ -322,7 +370,6 @@ namespace elsa
template <typename data_t>
typename DataContainer<data_t>::reverse_iterator DataContainer<data_t>::rbegin()
{
detach();
return reverse_iterator(end());
}
......@@ -341,7 +388,6 @@ namespace elsa
template <typename data_t>
typename DataContainer<data_t>::reverse_iterator DataContainer<data_t>::rend()
{
detach();
return reverse_iterator(begin());
}
......
......@@ -10,12 +10,6 @@
namespace elsa
{
// forward declaration for friend test function
template <typename data_t = real_t>
class DataContainer;
// used for testing and defined in test file (declared as friend)
template <typename data_t>
int useCount(const DataContainer<data_t>& /*dc*/);
/**
* \brief class representing and storing a linearized n-dimensional signal
......@@ -23,17 +17,15 @@ namespace elsa
* \author Matthias Wieczorek - initial code
* \author Tobias Lasser - rewrite, modularization, modernization
* \author David Frank - added DataHandler concept, iterators
* \author Nikola Dinev - add block support
*
* \tparam data_t - data type that is stored in the DataContainer, defaulting to real_t.
*
* 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.
*
* The class implements copy-on-write. Therefore any non-const functions should call the
* detach() function first to trigger the copy-on-write mechanism.
*/
template <typename data_t>
template <typename data_t = real_t>
class DataContainer
{
public:
......@@ -42,8 +34,6 @@ namespace elsa
*
* The following handler types are currently supported:
* - CPU: data is stored as an Eigen::Matrix in CPU main memory
* - MAP: data is not explicitly stored, but using an Eigen::Map to refer to other
* storage
*/
enum class DataHandlerType { CPU };
......@@ -104,7 +94,7 @@ namespace elsa
* fulfilled for any member functions, the object should not be used. After move- or copy-
* assignment, this is possible again.
*/
DataContainer<data_t>& operator=(DataContainer<data_t>&& other) noexcept;
DataContainer<data_t>& operator=(DataContainer<data_t>&& other);
/// return the current DataDescriptor
const DataDescriptor& getDataDescriptor() const;
......@@ -207,8 +197,17 @@ namespace elsa
/// comparison with another DataContainer
bool operator!=(const DataContainer<data_t>& other) const;
/// used for testing only and defined in test file
friend int useCount<>(const DataContainer<data_t>& dc);
/// 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
const DataContainer<data_t> getBlock(index_t i) const;
/// 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
const DataContainer<data_t> viewAs(const DataDescriptor& dataDescriptor) const;
/// iterator for DataContainer (random access and continuous)
using iterator = DataContainerIterator<DataContainer<data_t>>;
......@@ -280,7 +279,7 @@ namespace elsa
/// the current DataDescriptor
std::unique_ptr<DataDescriptor> _dataDescriptor;
/// the current DataHandler
std::shared_ptr<DataHandler<data_t>> _dataHandler;
std::unique_ptr<DataHandler<data_t>> _dataHandler;
/// factory method to create DataHandlers based on handlerType with perfect forwarding of
/// constructor arguments
......@@ -291,9 +290,6 @@ namespace elsa
/// private constructor accepting a DataDescriptor and a DataHandler
explicit DataContainer(const DataDescriptor& dataDescriptor,
std::unique_ptr<DataHandler<data_t>> dataHandler);
/// creates the deep copy for the copy-on-write mechanism
void detach();
};
/// element-wise addition of two DataContainers
......
......@@ -12,6 +12,7 @@ namespace elsa
*
* \author David Frank - initial code
* \author Tobias Lasser - modularization, modernization
* \author Nikola Dinev - add block support
*
* This abstract base class serves as an interface for data handlers, which encapsulate the
* actual data being stored e.g. in main memory of the CPU or in various memory types of GPUs.
......@@ -89,6 +90,36 @@ namespace elsa
/// assign a scalar to all elements of the data vector
virtual DataHandler<data_t>& operator=(data_t scalar) = 0;
/// copy assignment operator
DataHandler<data_t>& operator=(const DataHandler<data_t>& other)
{
if (other.getSize() != getSize())
throw std::invalid_argument("DataHandler: assignment argument has wrong size");
assign(other);
return *this;
}
/// move assignment operator
DataHandler<data_t>& operator=(DataHandler<data_t>&& other)
{
if (other.getSize() != getSize())
throw std::invalid_argument("DataHandler: assignment argument has wrong size");
assign(std::move(other));
return *this;
}
/// return a reference to the sequential block starting at startIndex and containing
/// numberOfElements elements
virtual std::unique_ptr<DataHandler<data_t>> getBlock(index_t startIndex,
index_t numberOfElements) = 0;
/// return a const reference to the sequential block starting at startIndex and containing
/// numberOfElements elements
virtual std::unique_ptr<const DataHandler<data_t>>
getBlock(index_t startIndex, index_t numberOfElements) const = 0;
protected:
/// slow element-wise dot product fall-back for when DataHandler types do not match
data_t slowDotProduct(const DataHandler<data_t>& v) const
......@@ -126,6 +157,19 @@ namespace elsa
for (index_t i = 0; i < getSize(); ++i)
(*this)[i] /= v[i];
}
/// slow element-wise assignment fall-back for when DataHandler types do not match
void slowAssign(const DataHandler<data_t>& other)
{
for (index_t i = 0; i < getSize(); ++i)
(*this)[i] = other[i];
}
/// derived classes should override this method to implement copy assignment
virtual void assign(const DataHandler<data_t>& other) = 0;
/// derived classes should override this method to implement move assignment
virtual void assign(DataHandler<data_t>&& other) = 0;
};
/// element-wise addition of two DataHandlers
......
#include "DataHandlerCPU.h"
#include "DataHandlerMapCPU.h"
namespace elsa
{
template <typename data_t>
DataHandlerCPU<data_t>::DataHandlerCPU(index_t size, bool initialize) : _data(size)
DataHandlerCPU<data_t>::DataHandlerCPU(index_t size, bool initialize)
: _data(std::make_shared<DataVector_t>(size))
{
if (initialize)
_data.setZero();
_data->setZero();
}
template <typename data_t>
DataHandlerCPU<data_t>::DataHandlerCPU(DataVector_t vector) : _data{vector}
DataHandlerCPU<data_t>::DataHandlerCPU(DataVector_t vector)
: _data{std::make_shared<DataVector_t>(vector)}
{
}
template <typename data_t>
DataHandlerCPU<data_t>::DataHandlerCPU(const DataHandlerCPU<data_t>& other)
: _data{other._data}, _associatedMaps{}
{
}
template <typename data_t>
DataHandlerCPU<data_t>::DataHandlerCPU(DataHandlerCPU<data_t>&& other)
: _data{std::move(other._data)}, _associatedMaps{std::move(other._associatedMaps)}
{
for (auto& map : _associatedMaps)
map->_dataOwner = this;
}
template <typename data_t>
DataHandlerCPU<data_t>::~DataHandlerCPU()
{
for (auto& map : _associatedMaps)
map->_dataOwner = nullptr;
}
template <typename data_t>
index_t DataHandlerCPU<data_t>::getSize() const
{
return static_cast<index_t>(_data.size());
return static_cast<index_t>(_data->size());
}
template <typename data_t>
data_t& DataHandlerCPU<data_t>::operator[](index_t index)
{
return _data[index];
detach();
return (*_data)[index];
}
template <typename data_t>
const data_t& DataHandlerCPU<data_t>::operator[](index_t index) const
{
return _data[index];
return (*_data)[index];
}
template <typename data_t>
......@@ -39,44 +64,45 @@ namespace elsa
if (v.getSize() != getSize())
throw std::invalid_argument("DataHandlerCPU: dot product argument has wrong size");
// if the other handler is not CPU, use the slow element-wise fallback version of dot
// product
auto otherHandler = dynamic_cast<const DataHandlerCPU*>(&v);
if (!otherHandler)
// use Eigen if the other handler is CPU or Map, otherwise use the slow fallback version
if (auto otherHandler = dynamic_cast<const DataHandlerCPU*>(&v)) {
return _data->dot(*otherHandler->_data);
} else if (auto otherHandler = dynamic_cast<const DataHandlerMapCPU<data_t>*>(&v)) {
return _data->dot(otherHandler->_map);
} else {
return this->slowDotProduct(v);
return _data.dot(otherHandler->_data);
}
}
template <typename data_t>
data_t DataHandlerCPU<data_t>::squaredL2Norm() const
{
return _data.squaredNorm();
return _data->squaredNorm();
}
template <typename data_t>
data_t DataHandlerCPU<data_t>::l1Norm() const
{
return _data.array().abs().sum();
return _data->array().abs().sum();
}
template <typename data_t>
data_t DataHandlerCPU<data_t>::lInfNorm() const
{
return _data.array().abs().maxCoeff();
return _data->array().abs().maxCoeff();
}
template <typename data_t>
data_t DataHandlerCPU<data_t>::sum() const
{
return _data.sum();
return _data->sum();
}
template <typename data_t>
std::unique_ptr<DataHandler<data_t>> DataHandlerCPU<data_t>::square() const
{
auto result = std::make_unique<DataHandlerCPU<data_t>>(getSize(), false);
result->_data = _data.array().square();
*result->_data = _data->array().square();
return result;
}
......@@ -85,7 +111,7 @@ namespace elsa
std::unique_ptr<DataHandler<data_t>> DataHandlerCPU<data_t>::sqrt() const
{
auto result = std::make_unique<DataHandlerCPU<data_t>>(getSize(), false);
result->_data = _data.array().sqrt();
*result->_data = _data->array().sqrt();
return result;
}
......@@ -94,7 +120,7 @@ namespace elsa
std::unique_ptr<DataHandler<data_t>> DataHandlerCPU<data_t>::exp() const
{
auto result = std::make_unique<DataHandlerCPU<data_t>>(getSize(), false);
result->_data = _data.array().exp();
*result->_data = _data->array().exp();
return result;
}
......@@ -103,7 +129,7 @@ namespace elsa
std::unique_ptr<DataHandler<data_t>> DataHandlerCPU<data_t>::log() const
{
auto result = std::make_unique<DataHandlerCPU<data_t>>(getSize(), false);
result->_data = _data.array().log();
*result->_data = _data->array().log();
return result;
}
......@@ -114,14 +140,17 @@ namespace elsa
if (v.getSize() != getSize())
throw std::invalid_argument("DataHandler: addition argument has wrong size");
// if the other handler is not CPU, use the slow element-wise fallback version of addition
auto otherHandler = dynamic_cast<const DataHandlerCPU*>(&v);
if (!otherHandler) {
detach();
// use Eigen if the other handler is CPU or Map, otherwise use the slow fallback version
if (auto otherHandler = dynamic_cast<const DataHandlerCPU*>(&v)) {
*_data += *otherHandler->_data;