wrappers.hpp

Go to the documentation of this file.

#ifndef IRLBA_WRAPPERS_HPP
#define IRLBA_WRAPPERS_HPP
 
#include "utils.hpp"
#include "Eigen/Dense"
#include <type_traits>
 
namespace irlba {
 
// This painful setup for the workspaces is because std::conditional requires
// both options to be valid types, so we can't just make a conditional to
// extract Matrix_::(Adjoint)Workspace for non-Eigen matrices (as this won't be
// valid for Eigen matrices that lack the typedef).
 
template<class Matrix_, typename = int>
struct get_workspace {
    typedef typename Matrix_::Workspace type;
};
 
template<class Matrix_>
struct get_workspace<Matrix_, decltype((void) std::declval<typename Matrix_::Index>(), 0)> {
    typedef bool type;
};
 
template<class Matrix_>
using WrappedWorkspace = typename get_workspace<Matrix_>::type;
 
template<class Matrix_>
WrappedWorkspace<Matrix_> wrapped_workspace(const Matrix_& matrix) {
    if constexpr(std::is_same<WrappedWorkspace<Matrix_>, bool>::value) {
        return false;
    } else {
        return matrix.workspace();
    }
}
 
template<class Matrix_, typename = int>
struct get_adjoint_workspace {
    typedef typename Matrix_::AdjointWorkspace type;
};
 
template<class Matrix_>
struct get_adjoint_workspace<Matrix_, decltype((void) std::declval<typename Matrix_::Index>(), 0)> {
    typedef bool type;
};
 
template<class Matrix_>
using WrappedAdjointWorkspace = typename get_adjoint_workspace<Matrix_>::type;
 
template<class Matrix_>
WrappedAdjointWorkspace<Matrix_> wrapped_adjoint_workspace(const Matrix_& matrix) {
    if constexpr(std::is_same<WrappedAdjointWorkspace<Matrix_>, bool>::value) {
        return false;
    } else {
        return matrix.adjoint_workspace();
    }
}
 
namespace internal {
 
template<class Matrix_, typename = int>
struct is_eigen {
    static constexpr bool value = false;
};
 
template<class Matrix_>
struct is_eigen<Matrix_, decltype((void) std::declval<typename Matrix_::Index>(), 0)> {
    static constexpr bool value = true;
};
 
}
template<class Matrix_, class Right_, class EigenVector_>
void wrapped_multiply(const Matrix_& matrix, const Right_& rhs, WrappedWorkspace<Matrix_>& work, EigenVector_& out) {
    if constexpr(internal::is_eigen<Matrix_>::value) {
        out.noalias() = matrix * rhs;
    } else {
        matrix.multiply(rhs, work, out);
    }
}
 
template<class Matrix_, class Right_, class EigenVector_>
void wrapped_adjoint_multiply(const Matrix_& matrix, const Right_& rhs, WrappedAdjointWorkspace<Matrix_>& work, EigenVector_& out) {
    if constexpr(internal::is_eigen<Matrix_>::value) {
        out.noalias() = matrix.adjoint() * rhs;
    } else {
        matrix.adjoint_multiply(rhs, work, out);
    }
}
 
template<class EigenMatrix_, class Matrix_>
EigenMatrix_ wrapped_realize(const Matrix_& matrix) {
    if constexpr(internal::is_eigen<Matrix_>::value) {
        return EigenMatrix_(matrix);
    } else {
        return matrix.template realize<EigenMatrix_>();
    }
}
 
template<class Matrix_, class EigenVector_>
class Centered {
public:
    Centered(const Matrix_& matrix, const EigenVector_& center) : my_matrix(matrix), my_center(center) {}
 
public:
    Eigen::Index rows() const { return my_matrix.rows(); }
 
    Eigen::Index cols() const { return my_matrix.cols(); }
 
public:
    struct Workspace {
        Workspace(WrappedWorkspace<Matrix_> i) : inner(std::move(i)) {}
        WrappedWorkspace<Matrix_> inner;
        EigenVector_ buffer;
    };
 
    Workspace workspace() const {
        return Workspace(wrapped_workspace(my_matrix));
    }
 
    struct AdjointWorkspace {
        AdjointWorkspace(WrappedAdjointWorkspace<Matrix_> i) : inner(std::move(i)) {}
        WrappedAdjointWorkspace<Matrix_> inner;
        EigenVector_ buffer;
    };
 
    AdjointWorkspace adjoint_workspace() const {
        return AdjointWorkspace(wrapped_adjoint_workspace(my_matrix));
    }
 
public:
    template<class Right_>
    void multiply(const Right_& rhs, Workspace& work, EigenVector_& out) const {
        const auto& realized_rhs = internal::realize_rhs(rhs, work.buffer);
        wrapped_multiply(my_matrix, realized_rhs, work.inner, out);
        auto beta = realized_rhs.dot(my_center);
        for (auto& o : out) {
            o -= beta;
        }
        return;
    }
 
    template<class Right_>
    void adjoint_multiply(const Right_& rhs, AdjointWorkspace& work, EigenVector_& out) const {
        const auto& realized_rhs = internal::realize_rhs(rhs, work.buffer);
        wrapped_adjoint_multiply(my_matrix, realized_rhs, work.inner, out);
        auto beta = realized_rhs.sum();
        out -= beta * (my_center);
        return;
    }
 
    template<class EigenMatrix_>
    Eigen::MatrixXd realize() const {
        auto output = wrapped_realize<EigenMatrix_>(my_matrix);
        output.array().rowwise() -= my_center.adjoint().array();
        return output;
    }
private:
    const Matrix_& my_matrix;
    const EigenVector_& my_center;
};
 
template<bool column_, class Matrix_, class EigenVector_>
class Scaled {
public:
    Scaled(const Matrix_& matrix, const EigenVector_& scale, bool divide) : 
        my_matrix(matrix), my_scale(scale), my_divide(divide) {}
 
public:
    Eigen::Index rows() const { return my_matrix.rows(); }
 
    Eigen::Index cols() const { return my_matrix.cols(); }
 
public:
    template<template<class> class Wrapper_>
    struct BufferedWorkspace {
        BufferedWorkspace(size_t n, Wrapper_<Matrix_> c) : buffer(n), child(std::move(c)) {}
        EigenVector_ buffer;
        Wrapper_<Matrix_> child;
    };
 
    typedef typename std::conditional<column_, BufferedWorkspace<WrappedWorkspace>, WrappedWorkspace<Matrix_> >::type Workspace;
 
    Workspace workspace() const {
        if constexpr(column_) {
            return BufferedWorkspace<WrappedWorkspace>(my_matrix.cols(), wrapped_workspace(my_matrix));
        } else {
            return wrapped_workspace(my_matrix);
        }
    }
 
    typedef typename std::conditional<column_, WrappedAdjointWorkspace<Matrix_>, BufferedWorkspace<WrappedAdjointWorkspace> >::type AdjointWorkspace;
 
    AdjointWorkspace adjoint_workspace() const {
        if constexpr(column_) {
            return wrapped_adjoint_workspace(my_matrix);
        } else {
            return BufferedWorkspace<WrappedAdjointWorkspace>(my_matrix.rows(), wrapped_adjoint_workspace(my_matrix));
        }
    }
 
public:
    template<class Right_>
    void multiply(const Right_& rhs, Workspace& work, EigenVector_& out) const {
        if constexpr(column_) {
            if (my_divide) {
                // We store the result here, because the underlying matrix's multiply()
                // might need to access rhs/scale multiple times, especially if it's
                // parallelized. Better to pay the cost of accessing a separate memory
                // space than computing the quotient repeatedly.
                work.buffer = rhs.cwiseQuotient(my_scale);
            } else {
                work.buffer = rhs.cwiseProduct(my_scale);
            }
            wrapped_multiply(my_matrix, work.buffer, work.child, out);
 
        } else {
            wrapped_multiply(my_matrix, rhs, work, out);
            if (my_divide) {
                out.array() /= my_scale.array();
            } else {
                out.array() *= my_scale.array();
            }
        }
    }
 
    template<class Right_>
    void adjoint_multiply(const Right_& rhs, AdjointWorkspace& work, EigenVector_& out) const {
        if constexpr(column_) {
            wrapped_adjoint_multiply(my_matrix, rhs, work, out);
            if (my_divide) {
                out.array() /= my_scale.array();
            } else {
                out.array() *= my_scale.array();
            }
 
        } else {
            if (my_divide) {
                work.buffer = rhs.cwiseQuotient(my_scale);
            } else {
                work.buffer = rhs.cwiseProduct(my_scale);
            }
            wrapped_adjoint_multiply(my_matrix, work.buffer, work.child, out);
        }
    }
 
    template<class EigenMatrix_>
    EigenMatrix_ realize() const {
        auto output = wrapped_realize<EigenMatrix_>(my_matrix);
 
        if constexpr(column_) {
            if (my_divide) {
                output.array().rowwise() /= my_scale.adjoint().array();
            } else {
                output.array().rowwise() *= my_scale.adjoint().array();
            }
        } else {
            if (my_divide) {
                output.array().colwise() /= my_scale.array();
            } else {
                output.array().colwise() *= my_scale.array();
            }
        }
 
        return output;
    }
private:
    const Matrix_& my_matrix;
    const EigenVector_& my_scale;
    bool my_divide;
};
 
template<bool column_, class Matrix_, class EigenVector_>
Scaled<column_, Matrix_, EigenVector_> make_Scaled(const Matrix_& matrix, const EigenVector_& scale, bool divide) {
    return Scaled<column_, Matrix_, EigenVector_>(matrix, scale, divide);
}
 
}
 
#endif