/*
Copyright (C) 2012, University of Cambridge, Department of Psychiatry.
Created by Rudolf Cardinal (rnc1001@cam.ac.uk).
This file is part of CamCOPS.
CamCOPS is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
CamCOPS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with CamCOPS. If not, see <https://www.gnu.org/licenses/>.
*/
#include "eigenfunc.h"
#include <cmath>
#include <QDebug>
#include <QtMath>
#include "maths/logisticregression.h"
using namespace Eigen;
namespace eigenfunc
{
// ============================================================================
// Making Eigen containers from std::vector
// ============================================================================
IndexArray makeIndexArray(const std::vector<Index>& v)
{
return eigenIndexVectorFromStdVector(v).array();
}
IndexArray makeIndexArray(std::initializer_list<Index> vlist)
{
return makeIndexArray(std::vector<Index>(vlist));
}
ArrayXb makeBoolArray(const std::vector<bool>& v)
{
return eigenColumnVectorFromStdVector<bool>(v).array();
}
ArrayXb makeBoolArray(std::initializer_list<bool> vlist)
{
return makeBoolArray(std::vector<bool>(vlist));
}
// ============================================================================
// Miscellaneous helpers
// ============================================================================
IndexArray indexSeq(const Index first, const Index last, const Index step)
{
const Index n = (1 + last - first) / step;
IndexArray indices(n);
Index idx = 0;
if (step > 0) {
for (Index i = first; i <= last; i += step) {
indices(idx++) = i;
}
} else if (step < 0) {
for (Index i = first; i >= last; i -= step) {
indices(idx++) = i;
}
}
return indices;
}
ArrayXb selectBoolFromIndices(const IndexArray& indices, const Index size)
{
ArrayXb select_bool(size);
select_bool.setConstant(false);
const Index n_indices = indices.size();
for (Index i = 0; i < n_indices; ++i) {
Eigen::Index idx = normalizeIndex(indices(i), size);
select_bool(idx) = true;
}
return select_bool;
}
ArrayXXb selectBoolFromIndices(const IndexArray& indices,
const Index n_rows,
const Index n_cols)
{
ArrayXXb select_bool(n_rows, n_cols);
select_bool.setConstant(false);
const Index n_indices = indices.size();
const Index size = select_bool.size();
Index row, col;
for (Index i = 0; i < n_indices; ++i) {
calcRowColFromIndex(indices(i), row, col, n_rows, size);
select_bool(row, col) = true;
}
return select_bool;
}
MatrixXd addOnesAsFirstColumn(const MatrixXd& m)
{
MatrixXd X_design(m.rows(), m.cols() + 1);
X_design << MatrixXd::Ones(m.rows(), 1.0), m; // first column is 1.0
return X_design;
}
// ============================================================================
// Subsetting Eigen arrays and matrices
// ============================================================================
Eigen::Index normalizeIndex(Eigen::Index idx, Eigen::Index size)
{
// Takes an index idx, makes sure it fits within size (cycling if need be),
// and returns it.
while (idx < 0) {
idx += size;
} // because (negative % positive) gives negative
return idx % size;
}
void calcRowColFromIndex(Eigen::Index idx,
Eigen::Index& row,
Eigen::Index& col,
Eigen::Index n_rows,
Eigen::Index size)
{
// Take an index idx; treat it using R's matrix approach, where the index
// increases down columns, then across rows, i.e. COLUMN-MAJOR ORDER:
//
// 0 3 6
// 1 4 7
// 2 5 8
//
// (which is like R except zero-based). Calculate the row and column
// indices (also zero-based), using "size" and "n_rows". Return them in
// "row" and "col".
idx = normalizeIndex(idx, size);
#ifdef __ANDROID__
std::ldiv_t division = std::ldiv(idx, n_rows);
#else
std::lldiv_t division = std::lldiv(idx, n_rows);
#endif
col = division.quot;
row = division.rem;
}
// ============================================================================
// Other R functions
// ============================================================================
MatrixXd scale(const MatrixXd& x,
const bool centre_on_column_mean,
const bool scale_divide_by_column_rms,
const ArrayXd& centre_values,
const ArrayXd& scale_values)
{
// To see R code:
// scale
// methods(scale)
// getAnywhere(scale.default)
Index nc = x.cols();
ArrayXXd xa = x.array();
// 1. Centre
if (centre_values.size() != nc && centre_values.size() != 0) {
qCritical() << Q_FUNC_INFO << "centre_values.size() is"
<< centre_values.size()
<< "which doesn't match number of columns" << nc;
} else if (centre_on_column_mean) {
// Centre each column on its mean
for (Index i = 0; i < nc; ++i) {
xa.col(i) -= xa.col(i).mean();
}
} else if (centre_values.size() == nc) {
// Centre each column on the specified value
for (Index i = 0; i < nc; ++i) {
xa.col(i) -= centre_values(i);
}
}
// 2. Scale
if (scale_values.size() != nc && scale_values.size() != 0) {
qCritical() << Q_FUNC_INFO << "scale_values.size() is"
<< scale_values.size()
<< "which doesn't match number of columns" << nc;
} else if (scale_divide_by_column_rms) {
// Scale each column on its mean
for (Index i = 0; i < nc; ++i) {
double divisor = std::sqrt(xa.col(i).square().mean());
// If "centre" is true, this will be the standard deviation.
// If not, it's just the root mean square.
// (This is because the SD is the RMS of deviations from the mean.)
xa.col(i) /= divisor;
}
} else if (scale_values.size() == nc) {
// Scale each column by dividing by the specified value
for (Index i = 0; i < nc; ++i) {
xa.col(i) /= scale_values(i);
}
}
// Note that even though xa is made with an ArrayWrapper() from x, it's not
// a "view" on x in the sense of a C++ reference.
// At this point, xa.matrix() is entirely different from x.
return xa.matrix();
}
Eigen::MatrixXd chol(const Eigen::MatrixXd& x, const bool pivot)
{
if (x.rows() != x.cols()) {
qFatal("Cholesky decomposition requires a SQUARE matrix");
}
// R's chol() returns the upper triangular matrix; see ?chol
// - https://stats.stackexchange.com/questions/117661/why-use-upper-triangular-cholesky
// - https://eigen.tuxfamily.org/dox/group__Cholesky__Module.html
if (pivot) {
// Eigen's LDLT: "Robust Cholesky decomposition of a matrix with pivoting"
LDLT<MatrixXd> c(x);
return c.matrixU(); // upper triangular
}
// Eigen's LLT: "Standard Cholesky decomposition (LL^T) of a matrix..."
LLT<MatrixXd> c(x);
return c.matrixU();
}
Eigen::MatrixXd backsolve(const Eigen::MatrixXd& r,
const Eigen::MatrixXd& x,
const Eigen::Index k,
const bool transpose,
const bool upper_tri)
{
return forwardOrBackSolve(r, x, k, transpose, upper_tri);
}
Eigen::MatrixXd forwardsolve(const Eigen::MatrixXd& l,
const Eigen::MatrixXd& x,
const Eigen::Index k,
const bool transpose,
const bool upper_tri)
{
return forwardOrBackSolve(l, x, k, transpose, upper_tri);
}
Eigen::MatrixXd forwardOrBackSolve(Eigen::MatrixXd lr,
Eigen::MatrixXd x,
Eigen::Index k,
const bool transpose,
bool upper_tri)
{
// - http://lists.r-forge.r-project.org/pipermail/rcpp-devel/2014-June/007781.html
// - http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20140627/2034608b/attachment.cpp
// - http://skip.ucsc.edu/leslie_MOUSE/r/R-1.8.1/src/library/base/R/backsolve.R
// PROPER SOURCE:
// - https://github.com/krlmlr/cxxr/blob/master/src/appl/bakslv.c
if (k < 0) {
k = lr.cols();
}
if (lr.rows() < k || lr.cols() < k) {
qCritical() << Q_FUNC_INFO << "lr too small (#rows or #cols < k)";
return MatrixXd();
}
if (x.rows() < k) {
qCritical() << Q_FUNC_INFO << "x too small (#rows < k)";
return MatrixXd();
}
if (transpose) {
// lr = lr.transpose() crashes:
// "aliasing detected during transposition, use transposeInPlace() or
// evaluate the rhs into a temporary using .eval()"
lr.transposeInPlace();
upper_tri = !upper_tri; // a hunch, but a correct one
// ... See
// https://github.com/krlmlr/cxxr/blob/master/src/appl/bakslv.c
// The phrasing is:
// if job is
// 00 solve t * x = b, t lower triangular,
// 01 solve t * x = b, t upper triangular,
// 10 solve t' * x = b, t lower triangular,
// 11 solve t' * x = b, t upper triangular.
}
if (k < lr.cols() || k < x.rows()) {
lr = lr.block(0, 0, lr.rows(), k);
x = x.block(0, 0, k, x.cols());
}
if (lr.cols() != x.rows()) {
qCritical() << Q_FUNC_INFO << "Size mismatch: lr.cols() != x.rows()";
return MatrixXd();
}
// qDebug() << Q_FUNC_INFO << "lr:" << qStringFromEigenMatrixOrArray(lr);
if (upper_tri) {
return lr.triangularView<Upper>().solve(x);
}
return lr.triangularView<Lower>().solve(x);
}
// ============================================================================
// Testing
// ============================================================================
const QString LINE("===============================================================================");
#define ASSERT_ARRAYS_SAME(a, b) Q_ASSERT((a).matrix() == (b).matrix());
#define OUTPUT_LINE() lines.append(LINE);
#define STATE(x) lines.append(x);
#define REPORT(x) lines.append(QString("%1: %2").arg(#x, qStringFromEigenMatrixOrArray(x)));
QStringList testEigenFunctions()
{
QStringList lines;
STATE("Testing eigenfunc...");
QVector<int> qv1{-1, 0, 1, 2};
OUTPUT_LINE();
VectorXi ev1_a = eigenColumnVectorFromQVector<int>(qv1);
VectorXi ev1_b(4); // must pre-size vector for <<
ev1_b << -1, 0, 1, 2;
Q_ASSERT(ev1_a == ev1_b);
Q_ASSERT(qVectorFromEigenVector<int>(ev1_a) == qv1);
STATE("Example column vector:");
REPORT(ev1_a);
RowVectorXi ev2_a = eigenColumnVectorFromQVector<int>(qv1);
RowVectorXi ev2_b(4); // must pre-size vector for <<
ev2_b << -1, 0, 1, 2;
Q_ASSERT(ev2_a == ev2_b);
Q_ASSERT(qVectorFromEigenVector<int>(ev2_a) == qv1);
STATE("Example row vector:");
REPORT(ev2_a);
IndexArray idxarr1_a = (IndexArray(3) << 3, 4, 5).finished();
IndexArray idxarr1_b = indexSeq(3, 5);
ASSERT_ARRAYS_SAME(idxarr1_a, idxarr1_b);
OUTPUT_LINE();
MatrixXi m1(4, 3);
m1 << 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12;
// m1: m1 "indices", cf. R:
// 1 2 3 0 4 8
// 4 5 6 1 5 9
// 7 8 9 2 6 10
// 10 11 12 3 7 11
IndexArray m1_which_a = which(m1.array() > 7);
IndexArray m1_which_b(5);
m1_which_b << 3, 6, 7, 10, 11;
REPORT(m1);
REPORT(m1_which_a);
REPORT(m1_which_b);
ASSERT_ARRAYS_SAME(m1_which_a, m1_which_b);
OUTPUT_LINE();
IndexArray m1_subset_cols = makeIndexArray({1, 2});
MatrixXi m2_a = subsetByColumnIndex(m1, m1_subset_cols);
MatrixXi m2_b(4, 2);
m2_b << 2, 3, 5, 6, 8, 9, 11, 12;
REPORT(m2_a);
Q_ASSERT(m2_a == m2_b);
OUTPUT_LINE();
IndexArray m1_subset_rows = makeIndexArray({1, 2});
MatrixXi m3_a = subsetByRowIndex(m1, m1_subset_rows);
MatrixXi m3_b(2, 3);
m3_b << 4, 5, 6, 7, 8, 9;
REPORT(m3_a);
Q_ASSERT(m3_a == m3_b);
OUTPUT_LINE();
IndexArray m1_subset_elements = makeIndexArray({1, 2, 10, 11});
ColumnArray<int> m4_a = subsetByElementIndex(m1, m1_subset_elements);
ColumnArray<int> m4_b(m1_subset_elements.size());
m4_b << 4, 7, 9, 12;
REPORT(m4_a);
ASSERT_ARRAYS_SAME(m4_a, m4_b);
OUTPUT_LINE();
ArrayXb m1_subset_cols_bool = makeBoolArray({false, true, true});
MatrixXi m5_a = subsetByColumnBoolean(m1, m1_subset_cols_bool);
REPORT(m5_a);
Q_ASSERT(m5_a == m2_b); // re-use
OUTPUT_LINE();
ArrayXb m1_subset_rows_bool = makeBoolArray({false, true, true, false});
MatrixXi m6_a = subsetByRowBoolean(m1, m1_subset_rows_bool);
REPORT(m6_a);
Q_ASSERT(m6_a == m3_b); // re-use
OUTPUT_LINE();
ArrayXXb m1_subset_elements_bool(4, 3);
m1_subset_elements_bool << // fill rows then cols
false, false, false,
true, false, false,
true, false, true,
false, false, true; // total of 4 true values
ColumnArray<int> m7_a = subsetByElementBoolean(m1, m1_subset_elements_bool);
REPORT(m7_a);
ASSERT_ARRAYS_SAME(m7_a, m4_b); // re-use
OUTPUT_LINE();
MatrixXi m8 = m1;
ArrayXi m9(4);
m9 << 100, 101, 102, 103;
MatrixXi m10(4, 3);
m10 << 1, 2, 3, 100, 5, 6, 101, 8, 102, 10, 11, 103;
assignByBooleanSequentially(m8, m1_subset_elements_bool, m9);
REPORT(m8);
REPORT(m9);
REPORT(m10);
Q_ASSERT(m8 == m10);
OUTPUT_LINE();
MatrixXi m11 = m1;
assignByIndexSequentially(m11, m1_subset_elements, m9);
REPORT(m11);
Q_ASSERT(m11 == m10); // re-use
OUTPUT_LINE();
MatrixXi m12 = sorted(m11);
MatrixXi m13(4, 3);
m13 << 1, 6, 100,
2, 8, 101,
3, 10, 102,
5, 11, 103;
REPORT(m12);
Q_ASSERT(m12 == m13);
OUTPUT_LINE();
MatrixXd m14 = m13.cast<double>(); // move into the double arena, not int
MatrixXd m15 = scale(m14);
STATE("[NOT TESTED BY AN ASSERT] Testing scale():");
REPORT(m15);
/*
R code:
m16 <- matrix(c(4, 12, -16, 12, 37, -43, -16, -43, 98), nrow=3)
chol(m16, pivot=FALSE)
chol(m16, pivot=TRUE)
*/
OUTPUT_LINE();
Matrix3d m16;
// https://en.wikipedia.org/wiki/Cholesky_decomposition#Example
m16 << 4, 12, -16,
12, 37, -43,
-16, -43, 98;
STATE("Matrix to undergo Cholesky decomposition:");
REPORT(m16);
Matrix3d m16_llt_u;
m16_llt_u << 2, 6, -8, // the U = Upper triangular version
0, 1, 5,
0, 0, 3;
// Matrix3d m16_ldlt_u;
// m16_ldlt_u << 1, 3, -4, // the U = Upper triangular version
// 0, 1, 5,
// 0, 0, 1;
STATE("Testing chol(pivot=false):");
Matrix3d m17 = chol(m16, false);
REPORT(m17);
Q_ASSERT(m17 == m16_llt_u);
STATE("[NOT TESTED BY AN ASSERT] Testing chol(pivot=true):");
Matrix3d m18 = chol(m16, true);
REPORT(m18);
// Q_ASSERT(m18 == m16_ldlt_u);
// Doesn't match Wikipedia; nor does it match R; nor does R match
// Wikipedia; nor are we using the pivot=false version...
OUTPUT_LINE();
STATE("Testing backsolve:");
// ?backsolve
Matrix3d r;
r << 1, 2, 3,
0, 1, 1,
0, 0, 2;
Vector3d x;
x << 8, 4, 2;
Vector3d backsolve_solution_a;
backsolve_solution_a << -1, 3, 1;
Vector3d backsolve_solution_b = backsolve(r, x);
REPORT(r);
REPORT(x);
REPORT(backsolve_solution_b);
Q_ASSERT(backsolve_solution_a == backsolve_solution_b);
Vector3d backsolve_tr_solution_a;
backsolve_tr_solution_a << 8, -12, -5;
Vector3d backsolve_tr_solution_b = backsolve(r, x, -1, true);
REPORT(backsolve_tr_solution_b);
Q_ASSERT(backsolve_tr_solution_a == backsolve_tr_solution_b);
Vector3d other_backsolve_solution;
other_backsolve_solution << 8, 4, 1;
// transpose, upper_tri:
Q_ASSERT(backsolve(r, x, -1, false, true) == backsolve_solution_a);
Q_ASSERT(backsolve(r, x, -1, true, true) == backsolve_tr_solution_a);
Q_ASSERT(backsolve(r, x, -1, true, false) == other_backsolve_solution);
Q_ASSERT(backsolve(r, x, -1, false, false) == other_backsolve_solution);
Matrix3d tr = r.transpose();
Q_ASSERT(backsolve(tr, x, -1, false, true) == other_backsolve_solution);
Q_ASSERT(backsolve(tr, x, -1, true, true) == other_backsolve_solution);
Q_ASSERT(backsolve(tr, x, -1, true, false) == backsolve_solution_a);
Q_ASSERT(backsolve(tr, x, -1, false, false) == backsolve_tr_solution_a);
OUTPUT_LINE();
STATE("Testing multiply:");
ArrayXXi m19(4, 4);
m19 << 1, 5, 9, 13,
2, 6, 10, 14,
3, 7, 11, 15,
4, 8, 12, 16;
ArrayXi m20(3);
m20 << 1, 10, 100;
ArrayXXi m21_a = multiply(m19, m20);
ArrayXXi m21_b = multiply(m20, m19);
ArrayXXi m21_c(4, 4);
m21_c << 1, 50, 900, 13,
20, 600, 10, 140,
300, 7, 110, 1500,
4, 80, 1200, 16;
REPORT(m19);
REPORT(m20);
REPORT(m21_a);
REPORT(m21_b);
ASSERT_ARRAYS_SAME(m21_a, m21_c);
ASSERT_ARRAYS_SAME(m21_b, m21_c);
OUTPUT_LINE();
STATE("... all eigenfunc tests completed correctly.");
return lines;
}
} // namespace eigenfunc