/* ----------------------------------------------------------------------- *//** * * @file matrix_agg.cpp * * @brief Build a matrix using the given column vectors * *//* ----------------------------------------------------------------------- */ #include #include #include #include "matrix_agg.hpp" namespace madlib { // Use Eigen using namespace dbal::eigen_integration; namespace modules { namespace linalg { /** * @brief Transition state for building a matrix * * Note: We assume that the DOUBLE PRECISION array is initialized by the * database with length 3, and all elemenets are 0. Handle::operator[] will * perform bounds checking. */ template class MatrixAggState { template friend class MatrixAggState; public: MatrixAggState(const AnyType &inArray) : mStorage(inArray.getAs()) { rebind(static_cast(mStorage[0]), static_cast(mStorage[1])); } operator AnyType() const { return mStorage; } void initialize(const Allocator& inAllocator, uint64_t inNumRows) { // Allocate the storage for 1 column mStorage = inAllocator.allocateArray(arraySize(inNumRows, 1)); rebind(inNumRows, 1); numRows = inNumRows; numCols = 0; } typename HandleTraits::MatrixTransparentHandleMap::ColXpr newColumn(const Allocator& inAllocator) { uint64_t numColsReserved = utils::nextPowerOfTwo( static_cast(numCols)); if (numColsReserved <= numCols) { if (numColsReserved == 0) numColsReserved = 1; else numColsReserved *= 2; // Shallow copy to save the references MatrixAggState oldSelf(*this); // Allocate new storage mStorage = inAllocator.allocateArray( arraySize(numRows, numColsReserved)); rebind(numRows, numColsReserved); // Copy from old storage to new storage numRows = oldSelf.numRows; numCols = oldSelf.numCols; matrix.leftCols(static_cast(numCols)) = oldSelf.matrix.leftCols(static_cast(numCols)); } // Not necessary for numCols = 0 rebind(numRows, numCols + 1); return matrix.col(static_cast(numCols++)); } private: static inline size_t arraySize(uint64_t inNumRows, uint64_t inNumCols) { return static_cast(2 + inNumRows * inNumCols); } /** * @brief Rebind to a new storage array * * @param inNumRows The number of rows * @param inNumCols The number of columns * * Array layout (iteration refers to one aggregate-function call): * Inter-iteration components (updated in final function): * - 0: numRows (number of rows) * - 1: numCols (number of columns) * - 2: sumOfPoints (matrix with \c numRows rows and \c numCols columns) */ void rebind(uint64_t inNumRows, uint64_t inNumCols) { numRows.rebind(&mStorage[0]); numCols.rebind(&mStorage[1]); matrix.rebind(&mStorage[2], static_cast(inNumRows), static_cast(inNumCols)); madlib_assert(mStorage.size() >= arraySize(inNumRows, inNumCols), std::runtime_error("Out-of-bounds array access detected.")); } Handle mStorage; public: typename HandleTraits::ReferenceToUInt64 numRows; typename HandleTraits::ReferenceToUInt64 numCols; typename HandleTraits::MatrixTransparentHandleMap matrix; }; AnyType matrix_agg_transition::run(AnyType& args) { MatrixAggState > state = args[0]; MappedColumnVector x = args[1].getAs(); if (state.numCols == 0) state.initialize(*this, x.size()); else if (x.size() != state.matrix.rows() || state.numRows != static_cast(state.matrix.rows()) || state.numCols != static_cast(state.matrix.cols())) throw std::invalid_argument("Invalid arguments: Dimensions of vectors " "not consistent."); state.newColumn(*this) = x; return state; } AnyType matrix_agg_final::run(AnyType& args) { MatrixAggState > state = args[0]; return MappedMatrix(state.matrix.leftCols( static_cast(state.numCols))); } AnyType matrix_column::run(AnyType& args) { MappedMatrix M = args[0].getAs(); uint32_t column = args[1].getAs(); return MappedColumnVector(M.col(column)); } } // namespace linalg } // namespace modules } // namespace regress