15.1.239. tablet_qt/maths/eigenfunc.h

/*
    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/>.
*/

#pragma once

#include <algorithm>
#include <cstdlib>
#include <functional>
#include <QDebug>
#include <QStringList>
#include <QVector>
#include <vector>
#include "maths/include_eigen_core.h"  // IWYU pragma: keep
#include "maths/include_eigen_dense.h"  // IWYU pragma: keep

namespace eigenfunc
{

/*

Reminders about Eigen types:

    Matrix<ContentsT, 3, 4>
        fixed rows=3 x cols=3 matrix of type ContentsT
    Matrix<ContentsT, Dynamic, 1>
        column vector, size unknown at compile-time, of type ContentsT
    MatrixXd
        shorthand for Matrix<double, Dynamic, Dynamic>

    Array<ContentsT, Dynamic, Dynamic>
        arbitrary 2x2 matrix of type ContentsT
    ArrayXXd
        shorthand for Array<double, Dynamic, Dynamic>
    ArrayXd
        shorthand for Array<double, Dynamic, 1>

The common base class of Matrix and Array is DenseBase.

SVD:

    with JacobiSVD<...> S:

        R       Eigen
        ----------------------------------------------
        S$d     S.singularValues()  // column vector
        S$u     S.matrixU()         // matrix
        S$v     S.matrixV()         // matrix

Inheritance and templates

    - ArrayBase<> and MatrixBase<> are both derived from DenseBase<> and
      ultimately from EigenBase<>

    - If you just want to use any old Eigen container based on its behaviour,
      you can use
            template<typename EigenContainerT>
            void doSomething(const EigenContainerT& m) { ... }

    - You can access the contents type from DenseBase<Derived> like this:
            template<typename Derived>
            void doSomething(const Eigen::DenseBase<Derived>& m)
            {
                using Scalar = typename Derived::Scalar;
                // ...
            }

    - But can you define a return value like this, where the thing you're
      returning uses the Scalar part, but is of a different type to the
      parameters? For example:

            Array<Scalar, Dynamic, 1> returnSomething(
                const ArrayOrMatrix<Scalar, ...>& m) { ... }

      Turns out you can:

            template<typename Derived>
            ColumnArray<typename Derived::Scalar> something(
                    const Eigen::DenseBase<Derived>& m)
            {
                Eigen::Array<typename Derived::Scalar, Eigen::Dynamic, 1> v;
                v << m(0, 0);
                return v;
            }

Conditional assignment

    No need for special functions for many cases; use

        X = boolean_array.select(Y, X)

    to achieve X = boolean_array ? Y : X, but elementwise.
    Note that either parameter of select() can be a scalar.
    HOWEVER, when it's a matrix, it must match by size.

    For size mismatch, we have assignByBoolean and assignByIndex.

Index

    Eigen::Index is of type std::ptrdiff_t
        http://eigen.tuxfamily.org/index.php?title=3.3#Index_typedef
        ... which is of type long (on my machine, anyway).
    So when dividing them, use
        std::ldiv_t division = std::ldiv(idx, nr);

Comma initializer

    - You have to pre-size a vector before using <<.
    - You can't insert elements one by one using <<; you can do

            v << 1, 2, 3, 4, 5;

      but not
            v << 1;
            v << 2;
            // ...

      ("Too few coefficients passed to comma initializer (operator<<)")

Comparison

    - It seems you can do
            vector1 == vector2

      but not
            array1 == array2

      ... which doesn't seem to match
      https://eigen.tuxfamily.org/dox/group__QuickRefPage.html

      ... or is it that you can compare, but only if the sizes are
      identical?

*/

// ============================================================================
// Type shorthands
// ============================================================================

template<typename ContentsT>
using ColumnVector = Eigen::Matrix<ContentsT, Eigen::Dynamic, 1>;

template<typename ContentsT>
using RowVector = Eigen::Matrix<ContentsT, 1, Eigen::Dynamic>;

template<typename ContentsT>
using ColumnArray = Eigen::Array<ContentsT, Eigen::Dynamic, 1>;

template<typename ContentsT>
using RowArray = Eigen::Array<ContentsT, 1, Eigen::Dynamic>;

template<typename ContentsT>
using GenericMatrix = Eigen::Matrix<ContentsT, Eigen::Dynamic, Eigen::Dynamic>;

template<typename ContentsT>
using GenericArray = Eigen::Matrix<ContentsT, Eigen::Dynamic, Eigen::Dynamic>;

using IndexArray = Eigen::Array<Eigen::Index, Eigen::Dynamic, 1>;  // 1-dimensional (column) array of indices
// Default storage is column-major, i.e. column vectors should be faster
// (though you can change this on a per-object basis);
// https://eigen.tuxfamily.org/dox/group__TopicStorageOrders.html
using IndexVector = Eigen::Matrix<Eigen::Index, Eigen::Dynamic, 1>;  // 1-dimensional (column) vector of indices

using ArrayXb = Eigen::Array<bool, Eigen::Dynamic, 1>;
// Eigen doesn't define ArrayXb, but it's helpful as a column vector of bool

using ArrayXXb = Eigen::Array<bool, Eigen::Dynamic, Eigen::Dynamic>;
// Eigen doesn't define ArrayXXb, but it's helpful as a n x n vector of bool


// ============================================================================
// Conversion between Qt and Eigen types
// ============================================================================

template<typename DestContainerT, typename SourceContentsT>
DestContainerT eigenVectorFromQVector(const QVector<SourceContentsT>& qv)
{
    // Takes a Qt QVector and returns an Eigen vector with the same contents.
    // The precise Eigen vector type (e.g. column vector, row vector) is
    // determined by the template.
    // This function is usually called indirectly, e.g. from
    //      eigenColumnVectorFromQVector()
    //      eigenRowVectorFromQVector()
    int n = qv.size();
    DestContainerT ev(n);
    for (int i = 0; i < n; ++i) {
        ev(i) = qv.at(i);
    }
    return ev;
}


template<typename DestContentsT, typename SourceContentsT>
ColumnVector<DestContentsT> eigenColumnVectorFromQVector(
        const QVector<SourceContentsT>& qv)
{
    // Takes a Qt QVector and returns an Eigen Vector (i.e. an Eigen Matrix
    // of dimensions nx1).
    return eigenVectorFromQVector<ColumnVector<DestContentsT>>(qv);
}


template<typename DestContentsT, typename SourceContentsT>
RowVector<DestContentsT> eigenRowVectorFromQVector(
        const QVector<SourceContentsT>& v)
{
    // Takes a Qt QVector and returns an Eigen RowVector (i.e. an Eigen Matrix
    // of dimensions 1xn).
    return eigenVectorFromQVector<RowVector<DestContentsT>>(v);
}


template<typename DestContentsT, typename SourceContainerT>
QVector<DestContentsT> qVectorFromEigenVector(const SourceContainerT& ev)
{
    // Takes an Eigen vector (row vector or column vector) and returns a Qt
    // QVector.
    int n = ev.size();
    QVector<DestContentsT> qv(n);
    for (int i = 0; i < n; ++i) {
        qv[i] = ev(i);
    }
    return qv;
}


template<typename Derived>
QString qStringFromEigenMatrixOrArray(const Eigen::DenseBase<Derived>& m,
                                      const QString type_name = "DenseBase")
{
    // Formats an Eigen Matrix or Array for display using a Qt QString.
    // - https://eigen.tuxfamily.org/dox/structEigen_1_1IOFormat.html
    // - C++ std::type_traits doesn't allow detection of template
    //   specializations, I don't think! Hence the type_name system, via which
    //   the detection is done at compile-time.
    Eigen::IOFormat heavy_fmt(
                Eigen::FullPrecision,  // precision
                0,  // flags
                ", ",  // coeffSeparator
                ";\n",  // rowSeparator
                "[",  // rowPrefix
                "]",  // rowSuffix
                "[",  // matPrefix
                "]");  // matSuffix
    std::stringstream ss;
    ss << m.format(heavy_fmt);
    QString description = QString("%1 (%2 rows x %3 cols)")
            .arg(type_name).arg(m.rows()).arg(m.cols());
    return description + "\n" + QString::fromStdString(ss.str());
}


// Specialization of qStringFromEigenMatrixOrArray for a Matrix.
template<typename Derived>
QString qStringFromEigenMatrixOrArray(const Eigen::MatrixBase<Derived>& m)
{
    return qStringFromEigenMatrixOrArray(m, "Matrix");
}


// Specialization of qStringFromEigenMatrixOrArray for an Array.
template<typename Derived>
QString qStringFromEigenMatrixOrArray(const Eigen::ArrayBase<Derived>& m)
{
    return qStringFromEigenMatrixOrArray(m, "Array");
}


// ============================================================================
// Making Eigen containers from std::vector
// ============================================================================

template<typename DestContainerT, typename SourceContentsT>
DestContainerT eigenVectorFromStdVector(const std::vector<SourceContentsT>& sv)
{
    // Takes a C++ std::vector and returns an Eigen vector with the same
    // contents. The precise Eigen vector type (e.g. column vector, row vector)
    // is determined by the template.
    // This function is usually called indirectly, e.g. from
    //      eigenColumnVectorFromStdVector()
    //      eigenRowVectorFromStdVector()
    //
    // Note that std::vector::size() returns size_t, but Eigen containers use
    // int for their dimensions. We need to make sure that everything is
    // compatible without any compiler warnings.
    // Remember: size_t is unsigned (but its size is implementation-specific);
    // int is of course signed (and at least 16 bits, but it could be larger).
    // More generally, re integer overflow in C/C++:
    //      http://www.cs.utah.edu/~regehr/papers/overflow12.pdf
    // and signed/unsigned comparisons:
    //      https://peter.bourgon.org/blog/2009/05/01/comparison-between-signed-and-unsigned-integer-expressions.html
    const size_t n_as_size_t = sv.size();
    const int n = static_cast<int>(n_as_size_t);
    // ... cannot be negative; if it is, there's overflow
    if (n < 0 || static_cast<size_t>(n) != n_as_size_t) {
        qWarning()
                << "Unable to represent std::vector of size_t" << n_as_size_t
                << "in Eigen container of int size" << n
                << "(likely large container overflowing int dimension). "
                   "RESULTS WILL BE WRONG.";
    }
    DestContainerT ev(n);
    for (int i = 0; i < n; ++i) {
        ev(i) = sv.at(i);
    }
    return ev;
}


template<typename DestContentsT, typename SourceContentsT>
ColumnVector<DestContentsT> eigenColumnVectorFromStdVector(
        const std::vector<SourceContentsT>& sv)
{
    // Takes a C++ std::vector and returns an Eigen Vector (i.e. an Eigen
    // Matrix of dimensions nx1).
    return eigenVectorFromStdVector<ColumnVector<DestContentsT>>(sv);
}


template<typename Scalar>
ColumnVector<Scalar> eigenColumnVectorFromInitList(
        std::initializer_list<Scalar> vlist)
{
    // Equivalent, but taking an initializer list
    return eigenColumnVectorFromStdVector<Scalar>(std::vector<Scalar>(vlist));
}


template<typename DestContentsT, typename SourceContentsT>
RowVector<DestContentsT> eigenRowVectorFromStdVector(
        const std::vector<SourceContentsT>& sv)
{
    // Takes a C++ std::vector and returns an Eigen RowVector (i.e. an Eigen
    // Matrix of dimensions 1xn).
    return eigenVectorFromStdVector<RowVector<DestContentsT>>(sv);
}


template<typename Scalar>
ColumnVector<Scalar> eigenRowVectorFromInitList(
        std::initializer_list<Scalar> vlist)
{
    // Equivalent, but taking an initializer list
    return eigenRowVectorFromStdVector<Scalar>(std::vector<Scalar>(vlist));
}


// A quick shorthand:
template<typename SourceContentsT>
ColumnVector<Eigen::Index> eigenIndexVectorFromStdVector(
        const std::vector<SourceContentsT>& sv)
{
    // Takes a C++ std::vector and returns an Eigen column Vector (Matrix)
    // containing Eigen::Index (= long int, typically).
    return eigenVectorFromStdVector<ColumnVector<Eigen::Index>>(sv);
}


template<typename DestContentsT, typename SourceContentsT>
ColumnArray<DestContentsT> eigenColumnArrayFromStdVector(
        const std::vector<SourceContentsT>& sv)
{
    // Takes a C++ std::vector and returns an Eigen column Array containing
    // the same contents.
    return eigenVectorFromStdVector<ColumnArray<DestContentsT>>(sv);
}


// Quick functions to make Eigen arrays of Eigen::Index or bool, from
// std::vector objects or initializer lists:
IndexArray makeIndexArray(const std::vector<Eigen::Index>& v);
IndexArray makeIndexArray(std::initializer_list<Eigen::Index> vlist);
ArrayXb makeBoolArray(const std::vector<bool>& v);
ArrayXb makeBoolArray(std::initializer_list<bool> vlist);


// ============================================================================
// Miscellaneous helpers
// ============================================================================

// Make a sequence of Eigen::Index, and return it in a column array:
IndexArray indexSeq(Eigen::Index first, Eigen::Index last,
                    Eigen::Index step = 1);

// Take a sequence of Eigen::Index, and return a sequence of bool for use as
// the condition in an Eigen "condition.select(then, _else)" statement. There's
// a version for vectors and another for 2d arrays:
ArrayXb selectBoolFromIndices(const IndexArray& indices, Eigen::Index size);
ArrayXXb selectBoolFromIndices(const IndexArray& indices, Eigen::Index n_rows,
                               Eigen::Index n_cols);

// Add a column of ones as the first column, for creating design matrices in
// which an intercept term is required.
Eigen::MatrixXd addOnesAsFirstColumn(const Eigen::MatrixXd& m);


// ============================================================================
// Subsetting Eigen arrays and matrices
// ============================================================================

// Takes an index idx, makes sure it fits within size (cycling if need be),
// and returns it.
Eigen::Index normalizeIndex(Eigen::Index idx, 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".
void calcRowColFromIndex(Eigen::Index idx,
                         Eigen::Index& row,
                         Eigen::Index& col,
                         Eigen::Index n_rows,
                         Eigen::Index size);


template<typename Derived>
void getRowColFromIndex(const Eigen::DenseBase<Derived>& m,
                        Eigen::Index idx,
                        Eigen::Index& row,
                        Eigen::Index& col)
{
    // Calculates row/col as per calcRowColFromIndex(), but taking an Eigen
    // Matrix/Array to work out the size from.

    const Eigen::Index nr = m.rows();
    const Eigen::Index size = m.size();
    calcRowColFromIndex(idx, row, col, nr, size);
}


template<typename Derived>
IndexArray which(const Eigen::DenseBase<Derived>& m)
{
    // As per R: produce a single vector of "indices"; go down columns before
    // going across rows (see calcRowColFromIndex, the converse). Tests the
    // truth value of every element in m, and put the index into the array if
    // its element is true.

    std::vector<Eigen::Index> indices;
    const Eigen::Index nr = m.rows();
    const Eigen::Index nc = m.cols();
    Eigen::Index i = 0;
    for (Eigen::Index c = 0; c < nc; ++c) {
        for (Eigen::Index r = 0; r < nr; ++r) {
            if (m(r, c)) {
                indices.push_back(i);
            }
            ++i;
        }
    }
    ColumnVector<Eigen::Index> v = eigenIndexVectorFromStdVector(indices);
    return v.array();
}


template<typename EigenContainerT>
EigenContainerT subsetByColumnIndex(const EigenContainerT& m,
                                    const IndexArray& column_indices)
{
    // Takes an Eigen Matrix/Array, and creates a corresponding object with
    // columns specified by "column_indices".
    // The number of output rows will be m.rows().
    // The number of output columns will be column_indices.size().

    const Eigen::Index m_ncols = m.cols();
    const Eigen::Index n_indices = column_indices.size();
    EigenContainerT subset(m.rows(), n_indices);
    for (Eigen::Index i = 0; i < n_indices; ++i) {
        Eigen::Index col_idx = normalizeIndex(column_indices(i), m_ncols);
        subset.col(i) = m.col(col_idx);
    }
    return subset;
}


template<typename EigenContainerT>
EigenContainerT subsetByRowIndex(const EigenContainerT& m,
                                 const IndexArray& row_indices)
{
    // Takes an Eigen Matrix/Array, and creates a corresponding object with
    // rows specified by "row_indices".
    // The number of output rows will be row_indices.size().
    // The number of output columns will be m.cols().

    const Eigen::Index m_nrows = m.rows();
    const Eigen::Index n_indices = row_indices.size();
    EigenContainerT subset(n_indices, m.cols());
    for (Eigen::Index i = 0; i < n_indices; ++i) {
        Eigen::Index row_idx = normalizeIndex(row_indices(i), m_nrows);
        subset.row(i) = m.row(row_idx);
    }
    return subset;
}


template<typename Derived>
ColumnArray<typename Derived::Scalar> subsetByElementIndex(
        const Eigen::DenseBase<Derived>& m,
        const IndexArray& indices)
{
    // Fetches elements of m by their index.
    // Treats "indices" as going down columns before going across rows.
    // Returns a column array, of size indices.size().

    const Eigen::Index nr = m.rows();
    const Eigen::Index oldlength = m.size();
    const Eigen::Index newlength = indices.size();
    ColumnArray<typename Derived::Scalar> result(newlength);
    for (Eigen::Index i = 0; i < newlength; ++i) {
        Eigen::Index idx = indices(i);
        Eigen::Index col, row;
        calcRowColFromIndex(idx, row, col, nr, oldlength);
        result(i) = m(row, col);
    }
    return result;
}


template<typename Derived>
GenericArray<typename Derived::Scalar> subsetByColumnBoolean(
        const Eigen::DenseBase<Derived>& m,
        const ArrayXb& use_column)
{
    // Takes an Eigen Matrix/Array, and creates an array whose columns are
    // specified by "use_column". Each element of "use_column" is true/false,
    // and columns are included or not depending on this.
    // CYCLES THROUGH use_column TO MAX LENGTH OF m.cols().
    // The number of output rows will be m.rows().
    // The number of output columns will be the number of "true" values
    // encounted while cycling through "use_column" up to m.cols().

    const Eigen::Index m_nc = m.cols();
    const Eigen::Index m_nr = m.rows();
    const Eigen::Index use_column_size = use_column.size();
    // Phase 1: calculate # columns
    Eigen::Index dest_nc = 0;
    for (Eigen::Index i = 0; i < m_nc; ++i) {
        Eigen::Index uc_idx = normalizeIndex(i, use_column_size);
        if (use_column(uc_idx)) {
            ++dest_nc;
        }
    }
    // Phase 2: assign
    GenericArray<typename Derived::Scalar> subset(m_nr, dest_nc);
    Eigen::Index dest_idx = 0;
    for (Eigen::Index i = 0; i < m_nc; ++i) {
        Eigen::Index uc_idx = normalizeIndex(i, use_column_size);
        if (use_column(uc_idx)) {
            subset.col(dest_idx++) = m.col(i);
        }
    }
    return subset;
}


template<typename Derived>
GenericArray<typename Derived::Scalar> subsetByRowBoolean(
        const Eigen::DenseBase<Derived>& m,
        const ArrayXb& use_row)
{
    // Takes an Eigen Matrix/Array, and creates an array whose rows are
    // specified by "use_row". Each element of "use_row" is true/false,
    // and columns are included or not depending on this.
    // CYCLES THROUGH use_row TO MAX LENGTH OF m.rows().
    // The number of output rows will be the number of "true" values
    // encounted while cycling through "use_row" up to m.rows().
    // The number of output columns will be m.cols().

    const Eigen::Index m_nc = m.cols();
    const Eigen::Index m_nr = m.rows();
    const Eigen::Index use_row_size = use_row.size();
    // Phase 1: calculate # rows
    Eigen::Index dest_nr = 0;
    for (Eigen::Index i = 0; i < m_nr; ++i) {
        Eigen::Index ur_idx = normalizeIndex(i, use_row_size);
        if (use_row(ur_idx)) {
            ++dest_nr;
        }
    }
    // Phase 2: assign
    GenericArray<typename Derived::Scalar> subset(dest_nr, m_nc);
    Eigen::Index dest_idx = 0;
    for (Eigen::Index i = 0; i < m_nr; ++i) {
        Eigen::Index ur_idx = normalizeIndex(i, use_row_size);
        if (use_row(ur_idx)) {
            subset.row(dest_idx++) = m.row(i);
        }
    }
    return subset;
}


template<typename Derived>
ColumnArray<typename Derived::Scalar> subsetByElementBoolean(
        const Eigen::DenseBase<Derived>& m,
        const ArrayXXb& which)
{
    // As per R (approximately): reads values of m for which "which" is true,
    // and spits them out into a vector, reading down columns first, then
    // across rows (i.e. column-wise, then row-wise).
    // Also following R: we don't care about the dimensionality of "which",
    // and will cycle through it. That is:
    // CYCLES THROUGH which TO SIZE OF m.

    using Scalar = typename Derived::Scalar;
    std::vector<Scalar> v;
    const Eigen::Index m_size = m.size();
    const Eigen::Index which_size = which.size();
    const Eigen::Index m_nr = m.rows();
    const Eigen::Index which_nr = which.rows();
    Eigen::Index m_row, m_col, which_row, which_col;
    for (Eigen::Index i = 0; i < m_size; ++i) {
        calcRowColFromIndex(i, m_row, m_col, m_nr, m_size);
        calcRowColFromIndex(i, which_row, which_col, which_nr, which_size);
        if (which(which_row, which_col)) {
            v.push_back(m(m_row, m_col));
        }
    }
    return eigenColumnArrayFromStdVector<Scalar>(v);
}


// ============================================================================
// Assigning to parts of Eigen object from source objects matching the "change"
// size, not the "recipient" size (for which select() works fine).
// Note that select() also works fine for assignment of scalars.
// ============================================================================

template <typename DerivedTo, typename DerivedFrom>
void assignByBooleanSequentially(Eigen::DenseBase<DerivedTo>& to,
                                 const ArrayXXb& which,
                                 const Eigen::DenseBase<DerivedFrom>& from)
{
    // Assigns values to "to", from "from", according to "which".
    // As a simple example:
    //      to    = [ 1  2  3  4  5  6  7  8  9 10]
    //      which = [ f  T  f  f  T  f  f  f  T  f]
    //      from  = [97 98 99]
    //
    //     result = [ 1 97  3  4 98  6  7  8 99 10]
    //
    // CYCLES THROUGH which TO MAX LENGTH OF to.
    // Also cycles through "from" as required.
    // (If "from" is two-dimensional, cycles through "from" in column-major
    // order.)

    const Eigen::Index num_to_replace = which.cast<int>().sum();
    if (num_to_replace == 0) {
        return;
    }
    const Eigen::Index from_size = from.size();
    if (from_size == 0) {
        qCritical() << Q_FUNC_INFO << "Empty 'from'";
        return;
    }
    if (from_size > num_to_replace || from_size % num_to_replace != 0) {
        qCritical() << Q_FUNC_INFO << "Number of items to replace is not a "
                                      "multiple of replacement length";
        return;
    }
    const Eigen::Index to_nr = to.rows();
    const Eigen::Index to_size = to.size();
    const Eigen::Index which_nr = which.rows();
    const Eigen::Index which_size = which.size();
    const Eigen::Index from_nr = from.rows();
    Eigen::Index to_row, to_col, which_row, which_col, from_row, from_col;
    Eigen::Index from_idx = 0;
    for (Eigen::Index i = 0; i < to_size; ++i) {
        calcRowColFromIndex(i, which_row, which_col, which_nr, which_size);
        if (!which(which_row, which_col)) {
            continue;
        }
        calcRowColFromIndex(i, to_row, to_col, to_nr, to_size);
        calcRowColFromIndex(from_idx, from_row, from_col, from_nr, from_size);
        to(to_row, to_col) = from(from_row, from_col);
        from_idx = normalizeIndex(++from_idx, from_size);
    }
}


template <typename DerivedTo, typename DerivedFrom>
void assignByIndexSequentially(Eigen::DenseBase<DerivedTo>& to,
                               const IndexArray& indices,
                               const Eigen::DenseBase<DerivedFrom>& from)
{
    // Assigns values to "to", according to element indices in "indices".
    // (Those indices are treated as down-columns-before-across-rows, i.e.
    // column-major order, zero-based.)
    // As a simple example:
    //      to      = [ 1  2  3  4  5  6  7  8  9 10]
    //      indices = [ 1 4 8 ]
    //      from    = [97 98 99]
    //
    //     result   = [ 1 97  3  4 98  6  7  8 99 10]
    //
    // Also cycles through "from" as required.

    const Eigen::Index n_indices = indices.size();
    const Eigen::Index from_size = from.size();
    if (n_indices == 0) {
        return;
    }

    // To mimic R behaviour:
    if (from_size > n_indices || from_size % n_indices != 0) {
        qCritical() << Q_FUNC_INFO << "Number of items to replace is not a "
                                      "multiple of replacement length";
        return;
    }

    const Eigen::Index to_nr = to.rows();
    const Eigen::Index to_size = to.size();
    Eigen::Index from_idx = 0;
    for (Eigen::Index i = 0; i < n_indices; ++i) {
        const Eigen::Index idx = indices(i);
        Eigen::Index row, col;
        calcRowColFromIndex(idx, row, col, to_nr, to_size);
        to(row, col) = from(from_idx++);
        if (from_idx >= from_size) {
            // Cycle round
            from_idx = 0;
        }
    }
}


// ============================================================================
// Array-by-vector elementwise operations, following R
// ============================================================================

template <typename Derived1, typename Derived2>
Eigen::Array<typename Derived1::Scalar,
             Eigen::Dynamic, Eigen::Dynamic> multiply(
        const Eigen::ArrayBase<Derived1>& a,
        const Eigen::ArrayBase<Derived2>& b)
{
    // Multiplies either (a) an array by an array of the same shape, or
    // (b) an array by a vector (not necessarily of the "right" length).
    //
    // # R example:
    // a = matrix(1:16, nrow=4)
    // b = c(1, 10, 100)
    // a * b
    //
    // # ... gives:
    // #     [,1] [,2] [,3] [,4]
    // # [1,]    1   50  900   13
    // # [2,]   20  600   10  140
    // # [3,]  300    7  110 1500
    // # [4,]    4   80 1200   16
    // # Warning message:
    // # In a * b : longer object length is not a multiple of shorter object length
    //
    // b * a  # same as a * b

    using ArrayType = Eigen::Array<typename Derived1::Scalar,
                                   Eigen::Dynamic, Eigen::Dynamic>;

    const Eigen::Index a_nr = a.rows();
    const Eigen::Index a_nc = a.cols();
    const Eigen::Index b_nr = b.rows();
    const Eigen::Index b_nc = b.cols();
    const Eigen::Index b_size = b.size();

    if (a_nr == b_nr && a_nc == b_nc) {
        // Arrays have the same dimensions
        return a * b;
    }

    bool a_is_vec = a_nr == 1 || a_nc == 1;
    bool b_is_vec = b_nr == 1 || b_nc == 1;
    if (!a_is_vec && !b_is_vec) {
        qFatal("multiply: non-conformable arrays");
    }
    if (!b_is_vec) {
        // It's very hard, in a template, to do some sort of neat swap.
        // For example, you can't create a reference or pointer and point to
        // either type (since the template isn't aware they are *of* the same
        // or a compatible type).
        return multiply(b, a);
    }

    ArrayType dest = a;
    Eigen::Index b_idx = 0;
    for (Eigen::Index c = 0; c < a_nc; ++c) {
        for (Eigen::Index r = 0; r < a_nr; ++r) {
            dest(r, c) *= b(b_idx++);
            b_idx %= b_size;  // cycle
        }
    }
    return dest;
}


#if 0
// ============================================================================
// Cyclic select()
// ============================================================================
// PROBLEM: mimicking R's   y = x[boolvec]  # achieved
//                          x[boolvec] = x  # achieved only indirectly:
// The problems are:
// (*) Eigen's select() requires conformable arrays.
// (*) R's x[boolvec], in contrast, cycles through boolvec up to the
//     length of x.
// (*) If we subset, e.g. with
//          const ArrayXd t_good = subsetByElementBoolean(t, good);
//     we have lots of short vectors floating around and it becomes hard
//     to be sure we're replacing them all correctly.
// Perhaps the answer is an extended select() function.
// See DenseBase.h, Select.h for the original.
// Let's write cyclicSelect();
//
// The original is condition.select(then, else).
// It can handle matrix-matrix, scalar-matrix, matrix-scalar

template <typename DerivedCondition, typename DerivedThenElse>
Eigen::Array<typename Derived::Scalar, Eigen::Dynamic, 1> cyclicSelect(
        const Eigen::DenseBase<DerivedCondition>& condition,
        const Eigen::DenseBase<DerivedThenElse>& then,
        const Eigen::DenseBase<DerivedThenElse>& _else)  // else is a C++ keyword
{
    // Produces a column array whose length is the maximum of the lengths of
    // (condition, then, else), cycling as necessary.

    using ArrayType = Eigen::Array<typename DerivedThenElse::Scalar,
                                   Eigen::Dynamic, Eigen::Dynamic>;

    // ABANDONED FOR NOW.
}
#endif


// ============================================================================
// Sorting Eigen things
// ============================================================================

template<typename Derived>
void sort(Eigen::PlainObjectBase<Derived>& m, bool decreasing = false)
{
    // Sorts an Eigen Matrix/Array in place.

    // An Array can be a view on a Matrix, I think, judging by
    //      array() functions at lines 326/329 of MatrixBase.h
    //          ... which return an ArrayWrapper to the derived() matrix object
    // The data() function is defined in PlainObjectBase.h and DenseStorage.h
    // Let's try using PlainObjectBase.
    using Scalar = typename Derived::Scalar;
    Scalar* p_data = m.data();
    Eigen::Index size = m.size();
    if (decreasing) {
        std::sort(p_data, p_data + size, std::greater<Scalar>());
    } else {
        std::sort(p_data, p_data + size);
    }
}


template<typename EigenContainerT>
EigenContainerT sorted(const EigenContainerT& vec)
{
    // Returns a sorted copy of an Eigen Matrix/Array.

    EigenContainerT newvec = vec;
    sort(newvec);
    return newvec;
}


// ============================================================================
// Other R functions
// ============================================================================
// See the identically named functions in R.

Eigen::MatrixXd scale(const Eigen::MatrixXd& x,
                      bool centre_on_column_mean = true,
                      bool scale_divide_by_column_rms = true,
                      const Eigen::ArrayXd& centre_values = Eigen::ArrayXd(),
                      const Eigen::ArrayXd& scale_values = Eigen::ArrayXd());

Eigen::MatrixXd chol(const Eigen::MatrixXd& x, bool pivot = false);

Eigen::MatrixXd backsolve(const Eigen::MatrixXd& r,
                          const Eigen::MatrixXd& x,
                          Eigen::Index k = -1,  // default to ncol(r)
                          bool transpose = false,
                          bool upper_tri = true);  // upper_tri = true MAKES IT A BACKSOLVE
Eigen::MatrixXd forwardsolve(const Eigen::MatrixXd& l,
                             const Eigen::MatrixXd& x,
                             Eigen::Index k = -1,  // default to ncol(l)
                             bool transpose = false,
                             bool upper_tri = false);  // upper_tri = false MAKES IT A FORWARDSOLVE
// I think that forwardsolve and backsolve are basically the same thing, except
// for the default relating to upper_tri. That's what R's documentation
// suggests:
//           Solves a system of linear equations where the coefficient matrix
//           is upper (or ‘right’, ‘R’) or lower (‘left’, ‘L’) triangular.
//           ‘x <- backsolve(R, b)’ solves R x = b [RNC: for x], and
//           ‘x <- forwardsolve(L, b)’ solves L x = b [RNC: for x], respectively.
// And this:
//      http://www.stat.berkeley.edu/~paciorek/research/techVignettes/techVignette6.pdf
//  "A quick note on backsolves and forwardsolves. If you have a
//  lower-triangular L, a forwardsolve is the calculation L^−1 b. A backsolve
//  is the calculation U^−1 b for upper-triangular U.
//  In R, we can do either type of solve with either L or U without an explicit
//  transpose:
//      U^−1 b: backsolve(U, b)
//      L^−1 b = U[T]^−1 b: backsolve(U, transpose = TRUE)
//      L^−1 b: forwardsolve(L, b)
//      U^−1 b = L[T]^−1 b: forwardsolve(L, b, transpose = TRUE)
// The "transpose" option:
// - In R help, it normally solves R x = b (giving X),
//   and the transpose option solves t(R) y = x (for y).
// - I wish it would be consistent within a single help page in its
//   use of variable names...

// My internal function to implement both forwardsolve() and backsolve():
Eigen::MatrixXd forwardOrBackSolve(Eigen::MatrixXd lr,
                                   Eigen::MatrixXd x,
                                   Eigen::Index k,
                                   bool transpose,
                                   bool upper_tri);


// ============================================================================
// Testing
// ============================================================================

// Test our Eigen functions.
QStringList testEigenFunctions();


}  // namespace eigenfunc