#ifndef CSC_H__

#define CSC_H__

#include "./shared.h"

#include "./timer.h"

#include "./blaze_adaptor.h"

#include "./merge.h"

#include "./io.h"

#include "./exception.h"

#include "thirdparty/mio.hpp"

#include <fstream>

namespace minicore {

namespace util {

template<typename T>

INLINE T abs_diff(T x, T y) {

if constexpr(std::is_unsigned_v<T>) {

return std::max(x, y) - std::min(x, y);

} else {

return std::abs(x - y);

}

}

template<typename T, typename T2>

INLINE auto abs_diff(T x, T2 y) {

using CT = std::common_type_t<T, T2>;

return abs_diff(CT(x), CT(y));

}

static constexpr size_t MINICORE_UTIL_ALN =

#ifdef __AVX512F__

sizeof(__m512) / sizeof(char);

#elif __AVX2__ || __AVX__

sizeof(__m256) / sizeof(char);

#elif __SSE4_1__ || __SSE2__

sizeof(__m128) / sizeof(char);

#else

1;

#endif

static inline bool is_file(std::string path) noexcept {

return ::access(path.data(), F_OK) != -1;

}

template<typename DataType>

struct ConstSViewMul: public std::pair<const DataType*, size_t> {

// Same as ConstSView, but multiplies the value

// by a constant during dereferencing

const DataType mul_;

ConstSViewMul(DataType mul): mul_(mul) {}

size_t index() const {return this->second;}

size_t &index() {return this->second;}

const DataType value() const {return *this->first * mul_;}

};

template<typename IndPtrType=uint64_t, typename IndicesType=uint64_t, typename DataType=uint32_t>

struct CSCMatrixView {

using ElementType = DataType;

IndPtrType *const indptr_;

IndicesType *const indices_;

DataType *const data_;

const uint64_t nnz_;

const uint32_t nf_, n_;

static_assert(std::is_integral_v<IndPtrType>, "IndPtr must be integral");

static_assert(std::is_integral_v<IndicesType>, "Indices must be integral");

static_assert(std::is_arithmetic_v<IndicesType>, "Data must be arithmetic");

CSCMatrixView(IndPtrType *indptr, IndicesType *indices, DataType *data,

uint64_t nnz, uint32_t nfeat, uint32_t nitems):

indptr_(indptr),

indices_(indices),

data_(data),

nnz_(nnz),

nf_(nfeat), n_(nitems)

{

}

struct Column {

const CSCMatrixView &mat_;

size_t start_;

size_t stop_;

Column(const CSCMatrixView &mat, size_t start, size_t stop)

: mat_(mat), start_(start), stop_(stop)

{

}

size_t nnz() const {return stop_ - start_;}

size_t size() const {return mat_.columns();}

template<bool is_const>

struct ColumnIteratorBase {

struct ViewType {

using VT = DataType;

ViewType(const ColumnIteratorBase &it): it_(it) {}

const ColumnIteratorBase &it_;

INLINE size_t index() const {return it_.col_.indices_[it_.index_];}

INLINE std::conditional_t<is_const, std::add_const_t<VT>, VT> &value() {return it_.col_.data_[it_.index_];}

INLINE std::add_const_t<VT> &value() const {return it_.col_.data_[it_.index_];}

};

using ColType = std::conditional_t<is_const, const Column, Column>;

using ViewedType = std::conditional_t<is_const, const DataType, DataType>;

using difference_type = std::ptrdiff_t;

using value_type = ViewedType;

using reference = ViewedType &;

using pointer = ViewedType *;

using iterator_category = std::random_access_iterator_tag;

ColType &col_;

size_t index_;

private:

mutable ViewType data_;

public:

template<bool oconst>

bool operator==(const ColumnIteratorBase<oconst> &o) const {

return index_ == o.index_;

}

template<bool oconst>

bool operator!=(const ColumnIteratorBase<oconst> &o) const {

return index_ != o.index_;

}

template<bool oconst>

bool operator<(const ColumnIteratorBase<oconst> &o) const {

return index_ < o.index_;

}

template<bool oconst>

bool operator>(const ColumnIteratorBase<oconst> &o) const {

return index_ > o.index_;

}

template<bool oconst>

bool operator<=(const ColumnIteratorBase<oconst> &o) const {

return index_ <= o.index_;

}

template<bool oconst>

bool operator>=(const ColumnIteratorBase<oconst> &o) const {

return index_ >= o.index_;

}

template<bool oconst>

difference_type operator-(const ColumnIteratorBase<oconst> &o) const {

return this->index_ - o.index_;

}

ColumnIteratorBase<is_const> &operator++() {

++index_;

return *this;

}

ColumnIteratorBase<is_const> operator++(int) {

ColumnIteratorBase ret(col_, index_);

++index_;

return ret;

}

const auto &operator*() const {

return data_;

}

auto &operator*() {

return data_;

}

ViewType *operator->() {

return &data_;

}

const ViewType *operator->() const {

return &data_;

}

ColumnIteratorBase(ColType &col, size_t ind): col_(col), index_(ind) {

}

};

using ColumnIterator = ColumnIteratorBase<false>;

using ConstColumnIterator = ColumnIteratorBase<true>;

ColumnIterator begin() {return ColumnIterator(*this, start_);}

ColumnIterator end() {return ColumnIterator(*this, stop_);}

ConstColumnIterator begin() const {return ConstColumnIterator(*this, start_);}

ConstColumnIterator end() const {return ConstColumnIterator(*this, stop_);}

};

auto column(size_t i) const {

return Column(*this, indptr_[i], indptr_[i + 1]);

}

size_t nnz() const {return nnz_;}

size_t rows() const {return n_;}

size_t columns() const {return nf_;}

};

template<typename T>

struct is_csc_view: public std::false_type {};

template<typename IndPtrType, typename IndicesType, typename DataType>

struct is_csc_view<CSCMatrixView<IndPtrType, IndicesType, DataType>>: public std::true_type {};

template<typename T>

static constexpr bool is_csc_view_v = is_csc_view<T>::value;

template<typename VT, typename IT>

struct CSparseVector {

using ElementType = VT;

static constexpr const bool transposeFlag = false;

VT *data_;

IT *indices_;

size_t n_, dim_;

//std::unique_ptr<VT> sum_;

CSparseVector(VT *data, IT *indices, size_t n, size_t dim): data_(data), indices_(indices), n_(n), dim_(dim)

{

}

size_t nnz() const {return n_;}

size_t size() const {return dim_;}

using NCVT = std::remove_const_t<VT>;

NCVT sum() const {

#if 0

std::remove_const_t<VT> ret;

auto di = reinterpret_cast<uint64_t>(data_);

if(di % MINICORE_UTIL_ALN && n_ > (MINICORE_UTIL_ALN / sizeof(VT)) && di % sizeof(VT) == 0) {

// Break into short unaligned + long aligned sum

const auto offset = (MINICORE_UTIL_ALN - (di % MINICORE_UTIL_ALN));

const auto offset_n = offset / sizeof(VT);

assert(reinterpret_cast<uint64_t>((NCVT *)data_ + offset_n) % MINICORE_UTIL_ALN == 0);

ret = blz::sum(blz::make_cv((NCVT *)data_, offset_n))

+ blz::sum(blz::make_cv<blz::aligned>((NCVT *)data_ + offset_n, n_ - offset_n));

} else ret = blz::sum(blz::make_cv<blz::aligned>((NCVT *)data_, n_));

return ret ;

#else

//std::fprintf(stderr, "Calling sum. value: %g\n", double(blz::sum(blz::make_cv((NCVT *)data_, n_))));

return blz::sum(blz::make_cv((NCVT *)data_, n_));

#endif

}

using DataType = VT;

template<bool is_const>

struct CSparseVectorIteratorBase {

using ColType = std::conditional_t<is_const, std::add_const_t<CSparseVector>, CSparseVector>;

using ViewedType = std::conditional_t<is_const, std::add_const_t<DataType>, DataType>;

using difference_type = std::ptrdiff_t;

using value_type = ViewedType;

using reference = ViewedType &;

using pointer = ViewedType *;

using iterator_category = std::random_access_iterator_tag;

ColType &col_;

size_t index_;

INLINE size_t index() const {return col_.indices_[index_];}

INLINE std::conditional_t<is_const, std::add_const_t<VT>, VT> &value() {return col_.data_[index_];}

INLINE std::add_const_t<VT> &value() const {return col_.data_[index_];}

using ViewType = CSparseVectorIteratorBase<is_const>;

template<bool oconst>

bool operator==(const CSparseVectorIteratorBase<oconst> &o) const {

return index_ == o.index_;

}

template<bool oconst>

bool operator!=(const CSparseVectorIteratorBase<oconst> &o) const {

return index_ != o.index_;

}

template<bool oconst>

bool operator<(const CSparseVectorIteratorBase<oconst> &o) const {

return index_ < o.index_;

}

template<bool oconst>

bool operator>(const CSparseVectorIteratorBase<oconst> &o) const {

return index_ > o.index_;

}

template<bool oconst>

bool operator<=(const CSparseVectorIteratorBase<oconst> &o) const {

return index_ <= o.index_;

}

template<bool oconst>

bool operator>=(const CSparseVectorIteratorBase<oconst> &o) const {

return index_ >= o.index_;

}

template<bool oconst>

difference_type operator-(const CSparseVectorIteratorBase<oconst> &o) const {

return this->index_ - o.index_;

}

CSparseVectorIteratorBase<is_const> &operator++() {

//std::fprintf(stderr, "before incrementing: indptr: %zu. index: %zu. value: %g\n", index_, size_t(col_.indices_[index_]), col_.data_[index_]);

++index_;

//std::fprintf(stderr, "after incrementing: indptr: %zu. index: %zu. value: %g\n", index_, size_t(col_.indices_[index_]), col_.data_[index_]);

return *this;

}

#if 0

CSparseVectorIteratorBase<is_const> operator++(int) {

CSparseVectorIteratorBase ret(col_, index_);

++index_;

return ret;

}

#endif

const ViewType &operator*() const {

return *this;

}

ViewType &operator*() {

return *this;

}

const ViewType *operator->() const {

//std::fprintf(stderr, "Calling const -> operator\n");

return this;

}

CSparseVectorIteratorBase(ColType &col, size_t ind): col_(col), index_(ind) {

}

};

double l2Norm() const {

double ret;

auto di = reinterpret_cast<uint64_t>(data_);

if(di % MINICORE_UTIL_ALN) {

if(n_ > (MINICORE_UTIL_ALN / sizeof(VT)) && di % sizeof(VT) == 0) {

// Break into short unaligned + long aligned sum

const auto offset = (MINICORE_UTIL_ALN - (di % MINICORE_UTIL_ALN));

const auto offset_n = offset / sizeof(VT);

assert(reinterpret_cast<uint64_t>((NCVT *)data_ + offset_n) % MINICORE_UTIL_ALN == 0);

ret = sqrNorm(blz::make_cv((NCVT *)data_, offset_n))

+ sqrNorm(blz::make_cv<blz::aligned>((NCVT *)data_ + offset_n, n_ - offset_n));

} else {

ret = sqrNorm(blz::make_cv((NCVT *)data_, n_));

}

} else ret = sqrNorm(blz::make_cv<blz::aligned>((NCVT *)data_, n_));

return std::sqrt(ret);

}

using ConstCSparseIterator = CSparseVectorIteratorBase<true>;

using CSparseIterator = CSparseVectorIteratorBase<false>;

CSparseIterator begin() {return CSparseIterator(*this, 0);}

CSparseIterator end() {return CSparseIterator(*this, n_);}

ConstCSparseIterator begin() const {return ConstCSparseIterator(*this, 0);}

ConstCSparseIterator end() const {return ConstCSparseIterator(*this, n_);}

};

template<typename VT, typename IT, typename IPtrT>

struct CSparseMatrix;

template<typename VT, typename IT>

std::ostream& operator<< (std::ostream& out, const CSparseVector<VT, IT> & item);

template<typename VT, typename IT, typename IPtrT>

std::ostream& operator<< (std::ostream& out, const CSparseMatrix<VT, IT, IPtrT> & item);

template<typename VT, typename IT>

struct ProdCSparseVector {

VT *data_;

IT *indices_;

const size_t n_;

const size_t dim_;

const double prod_;

ProdCSparseVector(const CSparseVector<VT, IT> &ovec, double prod): data_(ovec.data_), indices_(ovec.indices_), n_(ovec.n_), dim_(ovec.dim_), prod_(prod) {

}

ProdCSparseVector(const ProdCSparseVector<VT, IT> &ovec, double prod): data_(ovec.data_), indices_(ovec.indices_), n_(ovec.n_), dim_(ovec.dim_), prod_(prod * ovec.prod_) {

}

size_t nnz() const {return n_;}

size_t size() const {return dim_;}

using NCVT = std::remove_const_t<VT>;

double sum() const {

#if 0

double ret;

auto di = reinterpret_cast<uint64_t>(data_);

if(di % MINICORE_UTIL_ALN && n_ > (MINICORE_UTIL_ALN / sizeof(VT)) && di % sizeof(VT) == 0) {

// Break into short unaligned + long aligned sum

const auto offset = (MINICORE_UTIL_ALN - (di % MINICORE_UTIL_ALN));

const auto offset_n = offset / sizeof(VT);

assert(reinterpret_cast<uint64_t>((NCVT *)data_ + offset_n) % MINICORE_UTIL_ALN == 0);

ret = blz::sum(blz::make_cv((NCVT *)data_, offset_n))

+ blz::sum(blz::make_cv<blz::aligned>((NCVT *)data_ + offset_n, n_ - offset_n));

} else ret = blz::sum(blz::make_cv<blz::aligned>((NCVT *)data_, n_));

return ret * prod_;

#else

std::fprintf(stderr, "Calling sum. value: %g * %g = %g\n", double(blz::sum(blz::make_cv((NCVT *)data_, n_))), prod_, blz::sum(blz::make_cv((NCVT *)data_, n_)) * prod_);

return blz::sum(blz::make_cv((NCVT *)data_, n_)) * prod_;

#endif

}

using ConstCView = ConstSViewMul<VT>;

using DataType = VT;

struct ProdCSparseVectorIteratorBase {

using ColType = std::add_const_t<ProdCSparseVector>;

using ViewedType = std::add_const_t<DataType>;

using difference_type = std::ptrdiff_t;

using value_type = ViewedType;

using reference = ViewedType &;

using pointer = ViewedType *;

using iterator_category = std::random_access_iterator_tag;

using ViewType = ProdCSparseVectorIteratorBase;

ColType &col_;

size_t index_;

size_t index() const {return col_.indices_[index_];}

double value() const {

return col_.data_[index_] * col_.prod_;

}

public:

bool operator==(const ProdCSparseVectorIteratorBase &o) const {

return index_ == o.index_;

}

bool operator!=(const ProdCSparseVectorIteratorBase &o) const {

return index_ != o.index_;

}

bool operator<(const ProdCSparseVectorIteratorBase &o) const {

return index_ < o.index_;

}

bool operator>(const ProdCSparseVectorIteratorBase &o) const {

return index_ > o.index_;

}

bool operator<=(const ProdCSparseVectorIteratorBase &o) const {

return index_ <= o.index_;

}

bool operator>=(const ProdCSparseVectorIteratorBase &o) const {

return index_ >= o.index_;

}

difference_type operator-(const ProdCSparseVectorIteratorBase &o) const {

return this->index_ - o.index_;

}

ProdCSparseVectorIteratorBase &operator++() {

++index_;

return *this;

}

ProdCSparseVectorIteratorBase operator++(int) {

ProdCSparseVectorIteratorBase ret(col_, index_);

++index_;

return ret;

}

const ViewType &operator*() const {

return *this;

}

ViewType &operator*() {

return *this;

}

ViewType *operator->() {

return this;

}

const ViewType *operator->() const {

return this;

}

ProdCSparseVectorIteratorBase(ColType &col, size_t ind): col_(col), index_(ind) {

}

};

double l2Norm() const {

#if 0

double ret;

auto di = reinterpret_cast<uint64_t>(data_);

if(di % MINICORE_UTIL_ALN) {

if(n_ > (MINICORE_UTIL_ALN / sizeof(VT)) && di % sizeof(VT) == 0) {

// Break into short unaligned + long aligned sum

const auto offset = (MINICORE_UTIL_ALN - (di % MINICORE_UTIL_ALN));

const auto offset_n = offset / sizeof(VT);

assert(reinterpret_cast<uint64_t>((NCVT *)data_ + offset_n) % MINICORE_UTIL_ALN == 0);

ret = sqrNorm(blz::make_cv((NCVT *)data_, offset_n))

+ sqrNorm(blz::make_cv<blz::aligned>((NCVT *)data_ + offset_n, n_ - offset_n));

} else {

ret = sqrNorm(blz::make_cv((NCVT *)data_, n_));

}

} else ret = sqrNorm(blz::make_cv<blz::aligned>((NCVT *)data_, n_));

#else

double ret = blz::sqrNorm(blz::make_cv((NCVT *)data_, n_));

#endif

return prod_ * std::sqrt(ret);

}

using CSparseIterator = ProdCSparseVectorIteratorBase;

CSparseIterator begin() {return CSparseIterator(*this, 0);}

CSparseIterator end() {return CSparseIterator(*this, n_);}

CSparseIterator begin() const {return CSparseIterator(*this, 0);}

CSparseIterator end() const {return CSparseIterator(*this, n_);}

};

template<typename VT, typename IT>

inline double l2Norm(const CSparseVector<VT, IT> &x) {

return x.l2Norm();

}

template<typename VT, typename IT>

inline double l2Norm(const ProdCSparseVector<VT, IT> &x) {

return x.l2Norm();

}

template<typename VT, typename IT, typename OVT>

ProdCSparseVector<VT, IT> operator*(const CSparseVector<VT, IT> &lhs, OVT rhs) {

return ProdCSparseVector<VT, IT>(lhs, rhs);

}

template<typename VT, typename IT, typename OVT>

ProdCSparseVector<VT, IT> operator/(const CSparseVector<VT, IT> &lhs, OVT rhs) {

return lhs * double(1. / rhs);

}

template<typename VT, typename IT, typename OVT>

ProdCSparseVector<VT, IT> operator/(const ProdCSparseVector<VT, IT> &lhs, OVT rhs) {

return ProdCSparseVector<VT, IT>(lhs, 1. / rhs);

}

template<typename VT1, typename IT1, typename VT2, bool TF>

auto l2Dist(const ProdCSparseVector<VT1, IT1> &lhs, const blaze::DenseVector<VT2, TF> &rhs) {

if(lhs.size() != (*rhs).size()) throw std::invalid_argument("lhs and rhs have mismatched sizes");

auto &rr = *rhs;

using CT = std::common_type_t<VT1, blaze::ElementType_t<VT2>>;

CT ret = 0;

size_t si = 0, di = 0;

for(;si != lhs.n_ || di != (*rhs).size(); ++di) {

if(si == lhs.n_) {

ret += sum(subvector(*rhs, di, (*rhs).size() - di));

break;

} else if(di == (*rhs).size()) {

ret += sum(blz::make_cv(&lhs.data_[si], lhs.n_ - si)) * lhs.prod_;

break;

} else {

ret += sum(abs(subvector(*rhs, di, si - di)));

di = si;

ret += abs_diff((*rhs)[di], lhs.data_[si]);

++si;

}

}

return ret;

}

template<typename VT1, typename IT1, typename VT2, bool TF>

auto l2Dist(const CSparseVector<VT1, IT1> &lhs, const blaze::DenseVector<VT2, TF> &rhs) {

if(lhs.size() != (*rhs).size()) throw std::invalid_argument("lhs and rhs have mismatched sizes");

using CT = std::common_type_t<VT1, blaze::ElementType_t<VT2>>;

CT ret = 0;

size_t si = 0, di = 0;

for(;si != lhs.n_ || di != (*rhs).size(); ++di) {

if(si == lhs.n_) {

ret += sum(subvector(*rhs, di, (*rhs).size() - di));

break;

} else if(di == (*rhs).size()) {

ret += sum(blz::make_cv(&lhs.data_[si], lhs.n_ - si));

break;

} else {

while(di < lhs.indices_[si]) {

ret += std::abs((*rhs)[di++]);

}

ret += abs_diff((*rhs)[di], lhs.data_[si]);

++si;

}

}

return ret;

}

template<typename VT1, typename IT1, typename VT2, bool TF>

auto l2Dist(const blaze::DenseVector<VT2, TF> &rhs,

const CSparseVector<VT1, IT1> &lhs)

{

return l2Dist(lhs, rhs);

}

template<typename VT1, typename IT1, typename VT2, bool TF>

auto l2Dist(const blaze::DenseVector<VT2, TF> &rhs,

const ProdCSparseVector<VT1, IT1> &lhs)

{

return l2Dist(lhs, rhs);

}

template<typename VT1, typename IT1, typename VT2, bool TF>

auto l2Dist(const CSparseVector<VT1, IT1> &lhs, const blaze::SparseVector<VT2, TF> &rhs) {

if(lhs.size() != (*rhs).size()) throw std::invalid_argument("lhs and rhs have mismatched sizes");

auto &rr = *rhs;

using CT = std::common_type_t<VT1, blaze::ElementType_t<VT2>>;

CT ret = 0;

merge::for_each_by_case(lhs.size(), lhs.begin(), lhs.end(), rr.begin(), rr.end(),

[&ret](auto, auto lhv, auto rhv) {

auto v = abs_diff(lhv, rhv); ret += v * v;},

[&ret](auto, auto rhv) {ret += rhv * rhv;},

[&ret](auto, auto lhv) {ret += lhv * lhv;});

return ret;

}

template<typename VT1, typename IT1, typename VT2, bool TF>

auto l2Dist(const blaze::SparseVector<VT2, TF> &rhs, const CSparseVector<VT1, IT1> &lhs) {

return l2Dist(lhs, rhs);

}

template<typename VT1, typename IT1, typename VT2, bool TF>

auto l2Dist(const ProdCSparseVector<VT1, IT1> &lhs, const blaze::SparseVector<VT2, TF> &rhs) {

if(lhs.size() != (*rhs).size()) throw std::invalid_argument("lhs and rhs have mismatched sizes");

std::common_type_t<VT1, blz::ElementType_t<VT2>> ret = 0;

merge::for_each_by_case(lhs.size(), lhs.begin(), lhs.end(), (*rhs).begin(), (*rhs).end(),

[&ret](auto, auto lhv, auto rhv) {

auto v = abs_diff(lhv, rhv);

ret += v * v;},

[&ret](auto, auto rhv) {ret += rhv * rhv;},

[&ret](auto, auto lhv) {ret += lhv * lhv;});

return ret;

}

template<typename VT1, typename IT1, typename VT2, typename IT2>

auto l2Dist(const CSparseVector<VT1, IT1> &lhs, const CSparseVector<VT2, IT2> &rhs) {

if(lhs.size() != rhs.size()) throw std::invalid_argument("lhs and rhs have mismatched sizes");

std::common_type_t<VT1, VT2> ret = 0;

merge::for_each_by_case(lhs.size(), lhs.begin(), lhs.end(), rhs.begin(), rhs.end(),

[&ret](auto, auto lhv, auto rhv) {

auto v = abs_diff(lhv, rhv); ret += v * v;},

[&ret](auto, auto rhv) {ret += rhv * rhv;},

[&ret](auto, auto lhv) {ret += lhv * lhv;});

return ret;

}

template<typename VT1, typename IT1, typename VT2, typename IT2>

auto l2Dist(const ProdCSparseVector<VT1, IT1> &lhs, const CSparseVector<VT2, IT2> &rhs) {

if(lhs.size() != rhs.size()) throw std::invalid_argument("lhs and rhs have mismatched sizes");

std::common_type_t<VT1, VT2> ret = 0;

merge::for_each_by_case(lhs.size(), lhs.begin(), lhs.end(), rhs.begin(), rhs.end(),

[&ret](auto, auto lhv, auto rhv) {

auto v = abs_diff(lhv, rhv); ret += v * v;},

[&ret](auto, auto rhv) {ret += rhv * rhv;},

[&ret](auto, auto lhv) {ret += lhv * lhv;});

return ret;

}

template<typename VT1, typename IT1, typename VT2, typename IT2>

auto l2Dist(const CSparseVector<VT1, IT1> &lhs, const ProdCSparseVector<VT2, IT2> &rhs) {

return l2Dist(rhs, lhs);

}

template<typename VT1, typename IT1, typename VT2, typename IT2>

std::common_type_t<VT1, VT2> l2Dist(const ProdCSparseVector<VT1, IT1> &lhs, const ProdCSparseVector<VT2, IT2> &rhs) {

if(lhs.size() != rhs.size()) throw std::invalid_argument("lhs and rhs have mismatched sizes");

std::common_type_t<VT1, VT2> ret = 0;

merge::for_each_by_case(lhs.size(), lhs.begin(), lhs.end(), rhs.begin(), rhs.end(),

[&ret](auto, auto lhv, auto rhv) {

auto v = abs_diff(lhv, rhv); ret += v * v;},

[&ret](auto, auto rhv) {ret += rhv * rhv;},

[&ret](auto, auto lhv) {ret += lhv * lhv;});

return ret;

}

template<typename VT1, typename IT1, typename VT2, bool TF>

std::common_type_t<blz::ElementType_t<VT2>, VT1> l2Dist(const blaze::SparseVector<VT2, TF> &rhs, const ProdCSparseVector<VT1, IT1> &lhs) {

return l2Dist(lhs, rhs);

}

template<typename T1, typename T2>

auto sqrDist(const T1 &lhs, const T2 &rhs) {

auto ret = l2Dist(lhs, rhs);

return ret * ret;

}

template<typename T1, typename T2>

auto sqrl2Dist(const T1 &lhs, const T2 &rhs) {

return sqrDist(lhs, rhs);

}

template<typename T>

INLINE auto abs(T x) {

if constexpr(std::is_unsigned_v<T>) {

return x;

} else {

return std::abs(x);

}

}

template<typename VT1, typename IT1, typename VT2, bool TF>

auto l1Dist(const CSparseVector<VT1, IT1> &lhs, const blaze::SparseVector<VT2, TF> &rhs) {

//for(const auto &x: lhs) std::fprintf(stderr, "lhs %zu/%g\n", x.index(), x.value());

//for(const auto &x: *rhs) std::fprintf(stderr, "rhs %zu/%g\n", x.index(), x.value());

//std::fprintf(stderr, "%s l1dist with %zu/%zu sizes\n", __PRETTY_FUNCTION__, lhs.size(), (*rhs).size());

if(lhs.size() != (*rhs).size()) throw std::invalid_argument("lhs and rhs have mismatched sizes");

std::common_type_t<VT1, blaze::ElementType_t<VT2>, float> ret = 0;

merge::for_each_by_case(lhs.size(), lhs.begin(), lhs.end(), (*rhs).begin(), (*rhs).end(),

[&ret](auto, auto lhv, auto rhv) {ret += abs_diff(lhv, rhv);},

[&ret](auto, auto rhv) {ret += abs(rhv);},

[&ret](auto, auto lhv) {ret += abs(lhv);});

return ret;

}

template<typename VT1, typename IT1, typename VT2, bool TF>

auto l1Dist(const blaze::SparseVector<VT2, TF> &rhs, const CSparseVector<VT1, IT1> &lhs) {

return l1Dist(lhs, rhs);

}

template<typename VT1, typename IT1, typename VT2, bool TF>

auto l1Dist(const ProdCSparseVector<VT1, IT1> &lhs, const blaze::SparseVector<VT2, TF> &rhs) {

if(lhs.size() != (*rhs).size()) throw std::invalid_argument("lhs and rhs have mismatched sizes");

std::common_type_t<VT1, blz::ElementType_t<VT2>, float> ret = 0;

merge::for_each_by_case(lhs.size(), lhs.begin(), lhs.end(), (*rhs).begin(), (*rhs).end(),

[&ret](auto, auto lhv, auto rhv) {ret += abs_diff(lhv, rhv);},

[&ret](auto, auto rhv) {ret += abs(rhv);},

[&ret](auto, auto lhv) {ret += abs(lhv);});

return ret;

}

template<typename VT1, typename IT1, typename VT2, typename IT2>

auto l1Dist(const CSparseVector<VT1, IT1> &lhs, const CSparseVector<VT2, IT2> &rhs) {

if(lhs.size() != rhs.size()) throw std::invalid_argument("lhs and rhs have mismatched sizes");

std::common_type_t<VT1, VT2, float> ret = 0;

merge::for_each_by_case(lhs.size(), lhs.begin(), lhs.end(), rhs.begin(), rhs.end(),

[&ret](auto, auto lhv, auto rhv) {ret += abs_diff(lhv, rhv);},

[&ret](auto, auto rhv) {ret += abs(rhv);},

[&ret](auto, auto lhv) {ret += abs(lhv);});

return ret;

}

template<typename VT1, typename IT1, typename VT2, typename IT2>

auto l1Dist(const ProdCSparseVector<VT1, IT1> &lhs, const CSparseVector<VT2, IT2> &rhs) {

if(lhs.size() != rhs.size()) throw std::invalid_argument("lhs and rhs have mismatched sizes");

std::common_type_t<VT1, VT2, float> ret = 0;

merge::for_each_by_case(lhs.size(), lhs.begin(), lhs.end(), rhs.begin(), rhs.end(),

[&ret](auto, auto lhv, auto rhv) {ret += abs_diff(lhv, rhv);},

[&ret](auto, auto rhv) {ret += abs(rhv);},

[&ret](auto, auto lhv) {ret += abs(lhv);});

return ret;

}

template<typename VT1, typename IT1, typename VT2, typename IT2>

auto l1Dist(const CSparseVector<VT1, IT1> &lhs, const ProdCSparseVector<VT2, IT2> &rhs) {

return l1Dist(rhs, lhs);

}

template<typename VT1, typename IT1, typename VT2, typename IT2>

std::common_type_t<VT1, VT2> l1Dist(const ProdCSparseVector<VT1, IT1> &lhs, const ProdCSparseVector<VT2, IT2> &rhs) {

if(lhs.size() != rhs.size()) throw std::invalid_argument("lhs and rhs have mismatched sizes");

std::common_type_t<VT1, VT2, float> ret = 0;

merge::for_each_by_case(lhs.size(), lhs.begin(), lhs.end(), rhs.begin(), rhs.end(),

[&ret](auto, auto lhv, auto rhv) {ret += abs_diff(lhv, rhv);},

[&ret](auto, auto rhv) {ret += abs(rhv);},

[&ret](auto, auto lhv) {ret += abs(lhv);});

std::fprintf(stderr, "ret for l1Dist: %g\n", ret);

return ret;

}

template<typename VT1, typename IT1, typename VT2, bool TF>

std::common_type_t<blz::ElementType_t<VT2>, VT1> l1Dist(const blaze::SparseVector<VT2, TF> &rhs, const ProdCSparseVector<VT1, IT1> &lhs) {

return l1Dist(lhs, rhs);

}

template<typename T1, typename T2>

double dot_by_case(const T1 &lhs, const T2 &rhs) {

double ret = 0.;

merge::for_each_by_case(lhs.size(), lhs.begin(), lhs.end(), rhs.begin(), rhs.end(),

[&ret](auto, auto lhv, auto rhv) {ret += lhv * rhv;},

[&ret](auto, auto) {},

[&ret](auto, auto) {});

return ret;

}

template<typename VT, typename IT, typename VT2, typename IT2>

double dot(const CSparseVector<VT, IT> &lhs, const ProdCSparseVector<VT2, IT2> &rhs) {

return dot_by_case(lhs, rhs);

}

template<typename VT, typename IT, typename VT2, typename IT2>

double dot(const ProdCSparseVector<VT, IT> &lhs, const CSparseVector<VT2, IT2> &rhs) {

return dot(rhs, lhs);

}

template<typename VT, typename IT, typename VT2, typename IT2>

double dot(const CSparseVector<VT, IT> &lhs, const CSparseVector<VT2, IT2> &rhs) {

return dot_by_case(lhs, rhs);

}

template<typename VT, typename IT, typename VT2, typename IT2>

double dot(const ProdCSparseVector<VT, IT> &lhs, const ProdCSparseVector<VT2, IT2> &rhs) {

return dot_by_case(lhs, rhs);

}

template<typename VT, typename IT, typename VT2, bool TF>

double dot(const CSparseVector<VT, IT> &lhs, const blaze::SparseVector<VT2, TF> &rhs) {

return dot_by_case(lhs, *rhs);

}

template<typename VT, typename IT, typename VT2, bool TF>

double dot(const ProdCSparseVector<VT, IT> &lhs, const blaze::SparseVector<VT2, TF> &rhs) {

return dot_by_case(lhs, *rhs);

}

template<typename VT, typename IT, typename VT2, bool TF>

double dot(const blz::SparseVector<VT2, TF> &lhs, const CSparseVector<VT, IT> &rhs) {

return dot(rhs, *lhs);

}

template<typename VT, typename IT, typename VT2, bool TF>

double dot(const blz::SparseVector<VT2, TF> &lhs, const ProdCSparseVector<VT, IT> &rhs) {

return dot(rhs, *lhs);

}

template<typename T>

struct IsCSparseVector {

static constexpr bool value = false;

};

template<typename VT, typename IT>

struct IsCSparseVector<CSparseVector<VT, IT>>: public std::true_type {};

template<typename VT, typename IT>

struct IsCSparseVector<ProdCSparseVector<VT, IT>>: public std::true_type {};

template<typename T>

static constexpr const bool IsCSparseVector_v = IsCSparseVector<T>::value;

template<typename VT, typename IT>

auto make_csparse_view(VT *data, IT *idx, size_t n, size_t dim=-1) {

return CSparseVector<VT, IT>(data, idx, n, dim);

}

template<typename DataType, typename IndPtrType, typename IndicesType>

size_t nonZeros(const typename CSCMatrixView<IndPtrType, IndicesType, DataType>::Column &col) {

return col.nnz();

}

template<typename DataType, typename IndPtrType, typename IndicesType>

size_t nonZeros(const CSCMatrixView<IndPtrType, IndicesType, DataType> &mat) {

return mat.nnz();

}

template<typename VT, typename IT>

struct COOMatrixView {

IT *x, *y;

VT *data;

size_t nr_, nc_, nnz_;

size_t nnz() const {return nnz_;}

size_t rows() const {return nr_;}

size_t columns() const {return nc_;}

};

template<typename VT, typename IT>

struct COOMatrix {

std::vector<IT> x_, y_;

std::vector<VT> data_;

operator COOMatrixView<VT, IT> &() {

return COOMatrixView<VT, IT> {x_.data(), y_.data(), data_.data()};

}

operator COOMatrixView<const VT, const IT> &() const {

return COOMatrixView<const VT, const IT> {x_.data(), y_.data(), data_.data()};

}

size_t nr_, nc_;

size_t nnz() const {return x_.size();}

size_t rows() const {return nr_;}

size_t columns() const {return nc_;}

void add(IT x, IT y, VT data) {

x_.push_back(x); y_.push_back(y); data_.push_back(data_);

}

};

template<typename VT, typename IT, typename IPtrT>

struct CSparseMatrix {

VT *__restrict__ data_;

IT *__restrict__ indices_;

IPtrT *__restrict__ indptr_;

size_t nr_, nc_, nnz_;

constexpr CSparseMatrix(VT *__restrict__ data, IT *__restrict__ indices, IPtrT *__restrict__ indptr, size_t nr, size_t nc, size_t nnz):

data_(data), indices_(indices), indptr_(indptr), nr_(nr), nc_(nc), nnz_(nnz)

{

}

size_t nnz() const {return nnz_;}

size_t rows() const {return nr_;}

size_t columns() const {return nc_;}

auto row(size_t i) {

return CSparseVector<VT, IT>(data_ + indptr_[i], indices_ + indptr_[i], indptr_[i + 1] - indptr_[i], nc_);

}

auto row(size_t i) const {

return CSparseVector<const VT, IT>(data_ + indptr_[i], indices_ + indptr_[i], indptr_[i + 1] - indptr_[i], nc_);

}

auto sum() const {

return blz::sum(blz::CustomVector<VT, blaze::unaligned, blaze::unpadded>(data_, nnz_));

}

auto &operator~() {return *this;}

const auto &operator~() const {return *this;}

auto &operator*() {return *this;}

const auto &operator*() const {return *this;}

using ElementType = VT;

};

using blaze::unchecked;

template<typename VT, typename ORVT, typename IT, typename IPtr, bool TF, typename OIT=IT, typename WeightT=blz::DV<VT>, bool rowwise=true>

void l1_median(const CSparseMatrix<VT, IT, IPtr> &mat, blaze::Vector<ORVT, TF> &ret, OIT *asnptr=static_cast<OIT *>(nullptr), size_t asnsz=0, WeightT *weightc=static_cast<WeightT *>(nullptr)) {

if constexpr(rowwise) {

using RVT = blaze::ElementType_t<ORVT>;

shared::flat_hash_map<IT, std::vector<RVT>> nzfeatures;

nzfeatures.reserve(mat.columns());

const size_t nrows = asnsz ? asnsz: mat.rows();

if(weightc) throw NotImplementedError("Not implemented: weighted L1 median for CSparseMatrix");

for(size_t i = 0; i < nrows; ++i) {

auto rownum = asnsz ? OIT(asnptr[i]): OIT(i);

auto r = row(mat, rownum, unchecked);

auto nnz = nonZeros(r);

if(!nnz) continue;

for(size_t i = 0; i < r.n_; ++i) {

auto idx = r.indices_[i];

auto val = r.data_[i];

auto it = nzfeatures.find(idx);

if(it == nzfeatures.end()) {

it = nzfeatures.emplace(idx, std::vector<RVT>{RVT(val)}).first;

} else {

it->second.push_back(val);

}

}

}

const size_t nfeat = nzfeatures.size();

(*ret).reset();

(*ret).resize(mat.nc_);

(*ret).reserve(nfeat);

const bool nr_is_odd = nrows & 1;

for(auto &pair: nzfeatures) {

auto it = &pair.second;

auto i = pair.first;

auto ibeg = it->begin(), iend = it->end();

shared::sort(ibeg, iend);

size_t nel = it->size();

if(nel > nrows / 2) {

if(it->front() > 0.) {

if(nr_is_odd) {

auto mid = ibeg + nrows / 2 - (nrows - nel);

(*ret)[i] = *mid;

} else {

auto midm1 = ibeg + nrows / 2 - (nrows - nel), mid = midm1 + 1;

(*ret)[i] = .5 * (*mid + *midm1);

}

} else if(it->back() < 0.) {

if(nr_is_odd) {

(*ret)[i] = *(iend - nrows / 2 + (nrows - nel));

} else {

auto midm1 = iend - nrows / 2 + (nrows - nel), mid = midm1 + 1;

(*ret)[i] = .5 * (*mid + *midm1);

}

} else {

it->resize(nrows, RVT(0));

shared::sort(it->begin(), it->end());

if(nr_is_odd) (*ret)[i] = (*it)[nrows / 2];

else (*ret)[i] = .5 * ((*it)[nrows / 2] + (*it)[nrows / 2 + 1]);

}

} else if((nrows & 1) && nel == nrows / 2) {

if(it->front() > 0.) (*ret)[i] = .5 * it->front();

else if(it->back() < 0.) (*ret)[i] = .5 * it->back();

} // else the value is 0

}

} else {

throw std::runtime_error("Not implemented");

}

}

template<typename VT, typename ORVT, typename IT, typename IPtr, bool TF, typename OIT, typename WeightT=blz::DV<VT>, bool rowwise=true, typename RSums>

void tvd_median(const CSparseMatrix<VT, IT, IPtr> &mat, blaze::SparseVector<ORVT, TF> &ret, OIT *asnptr=static_cast<OIT *>(nullptr), size_t asnsz=0, WeightT *weightc=static_cast<WeightT *>(nullptr), const RSums &rsums=RSums()) {

if constexpr(rowwise) {

using RVT = blaze::ElementType_t<ORVT>;

shared::flat_hash_map<IT, std::vector<RVT>> nzfeatures;

nzfeatures.reserve(mat.columns());

const size_t nrows = asnsz ? asnsz: mat.rows();

if(weightc) throw NotImplementedError("Not implemented");

for(size_t i = 0; i < nrows; ++i) {

auto rownum = asnsz ? OIT(asnptr[i]): OIT(i);

auto r = row(mat, rownum, unchecked);

auto nnz = nonZeros(r);

if(!nnz) continue;

const typename std::conditional_t<(sizeof(VT) <= 4), float, double> invmul = 1. / rsums[rownum];

for(size_t i = 0; i < r.n_; ++i) {

auto idx = r.indices_[i];

const RVT val = RVT(r.data_[i]) * invmul;

auto it = nzfeatures.find(idx);

if(it == nzfeatures.end()) {

it = nzfeatures.emplace(idx, std::vector<RVT>{val}).first;

} else {

it->second.push_back(val);

}

}

}

const size_t nfeat = nzfeatures.size();

(*ret).reset();

(*ret).resize(mat.nc_);

(*ret).reserve(nfeat);

const bool nr_is_odd = nrows & 1;

for(auto &pair: nzfeatures) {

auto it = &pair.second;

auto i = pair.first;

auto ibeg = it->begin(), iend = it->end();

shared::sort(ibeg, iend);

size_t nel = it->size();

if(nel > nrows / 2) {

if(it->front() > 0.) {

#define append insert

if(nr_is_odd) {

auto mid = ibeg + nrows / 2 - (nrows - nel);

(*ret).append(i, *mid);

} else {

auto midm1 = ibeg + nrows / 2 - (nrows - nel), mid = midm1 + 1;

(*ret).append(i, .5 * (*mid + *midm1));

}

} else if(it->back() < 0.) {

if(nr_is_odd) {

auto mid = iend - nrows / 2 + (nrows - nel);

(*ret).append(i, *mid);

} else {

auto midm1 = ibeg - nrows / 2 + (nrows - nel), mid = midm1 + 1;

(*ret).append(i, (.5) * (*mid + *midm1));

}

} else {

it->resize(nrows, RVT(0));

shared::sort(it->begin(), it->end());

if(nr_is_odd) (*ret).append(i, (*it)[nrows / 2]);

else (*ret).append(i, .5 * ((*it)[nrows / 2 - 1] + (*it)[nrows / 2]));

}

} else if((nrows & 1) && nel == nrows / 2) {

if(it->front() > 0.) (*ret).append(i, .5 * it->front());

else if(it->back() < 0.) (*ret).append(i, .5 * it->back());

} // else it's zero.

#undef append

}

} else {

throw std::runtime_error("Not implemented");

}

}

template<typename VT, typename IT, typename IPtrT>

inline CSparseMatrix<VT, IT, IPtrT> make_csparse_matrix(VT *__restrict__ data, IT *__restrict__ indices, IPtrT *__restrict__ indptr, size_t nr, size_t nc, size_t nnz) {

return CSparseMatrix<VT, IT, IPtrT>(data, indices, indptr, nr, nc, nnz);

}

template<typename VT, typename IT, typename IPtrT, bool checked>

inline auto row(const CSparseMatrix<VT, IT, IPtrT> &mat, size_t i, blaze::Check<checked>) {

if constexpr(checked) {

if(unlikely(i > mat.rows())) throw std::out_of_range(std::string("Out of range: row ") + std::to_string(i) + " out of " + std::to_string(mat.rows()));

}

return mat.row(i);

}

template<typename VT, typename IT, typename IPtrT>

inline auto row(const CSparseMatrix<VT, IT, IPtrT> &mat, size_t i) {

return row(mat, i, blaze::Check<true>());

}

template<typename VT, typename IT, typename IPtrT, bool checked>

inline auto row(CSparseMatrix<VT, IT, IPtrT> &mat, size_t i, blaze::Check<checked>) {

if constexpr(checked) {

if(unlikely(i > mat.rows())) throw std::out_of_range(std::string("Out of range: row ") + std::to_string(i) + " out of " + std::to_string(mat.rows()));

}

return mat.row(i);

}

template<typename VT, typename IT, typename IPtrT>

inline auto row(CSparseMatrix<VT, IT, IPtrT> &mat, size_t i) {

return row(mat, i, blaze::Check<true>());

}

template<typename VT, typename IT, typename IPtrT>

INLINE auto sum(const CSparseMatrix<VT, IT, IPtrT> &sm) {

return sm.sum();

}

template<typename VT, typename IT>

INLINE auto sum(const ProdCSparseVector<VT, IT> &sm) {