Branch data Line data Source code
1 : : #include "RDC.hpp"
2 : :
3 : : #include <Eigen/Dense>
4 : : #include <mutable/util/macro.hpp>
5 : : #include <vector>
6 : :
7 : :
8 : : using namespace m;
9 : : using namespace Eigen;
10 : :
11 : :
12 : : /*======================================================================================================================
13 : : * Helper functions
14 : : *====================================================================================================================*/
15 : :
16 : 838 : MatrixXf m::create_CDF_matrix(const MatrixXf &data)
17 : : {
18 : : /* Allocate CDF matrix. */
19 : 838 : MatrixXf CDF_matrix(data.rows(), data.cols() + 1); // +1 for y-intercepts in non-linear transformations in RDC
20 : :
21 : : /* Allocate index vector and CDF vector. */
22 : 838 : std::vector<unsigned> indices;
23 [ + - ]: 838 : indices.resize(data.rows());
24 : :
25 [ + + ]: 1676 : for (unsigned col_id = 0; col_id != data.cols(); ++col_id) {
26 : : /* Initialize index vector. */
27 [ + - ]: 838 : M_insist(indices.size() == (unsigned long)(data.rows()));
28 [ + - ]: 838 : std::iota(indices.begin(), indices.end(), 0);
29 : :
30 : : /* Compute index vector that sorts the column. */
31 [ + - + - ]: 96148 : std::sort(indices.begin(), indices.end(), [col=data.col(col_id)](unsigned first, unsigned second) {
32 : 95310 : return col(first) < col(second);
33 : : });
34 : :
35 : : /* Translate index vector to CDF. */
36 : 838 : unsigned last_pos = data.rows() - 1;
37 [ + - ]: 838 : float last_value = data(indices[last_pos], col_id);
38 : 839 : unsigned num_same = 0;
39 [ + + ]: 24698 : for (unsigned i = data.rows(); i-->0;) {
40 [ + - + + ]: 23860 : if (data(indices[i], col_id) != last_value) {
41 [ + - ]: 21613 : last_value = data(indices[i], col_id);
42 : 21613 : last_pos -= num_same;
43 : 21613 : num_same = 0;
44 : 21613 : }
45 [ + - ]: 23860 : CDF_matrix(indices[i], col_id) = float(last_pos + 1) / data.rows();
46 : 23860 : ++num_same;
47 : : }
48 : 838 : }
49 : :
50 : : /* Add last, all-ones column for Y-intercept of non-liner transformations, used later for RDC computation. */
51 [ + - + - : 838 : CDF_matrix.col(data.cols()) = VectorXf::Ones(data.rows());
+ - ]
52 : :
53 : 838 : return CDF_matrix;
54 [ + - ]: 838 : }
55 : :
56 : : namespace {
57 : 1 :
58 : : template<typename URBG>
59 : 880 : MatrixXf create_w_b(float stddev, unsigned num_rows, unsigned k, URBG &&g)
60 : : {
61 : 880 : std::normal_distribution<float> normal_distribution(0.f, stddev);
62 : 880 : MatrixXf w_b(num_rows, k);
63 [ + + ]: 5280 : for (unsigned col_id = 0; col_id != k; ++col_id) {
64 [ + - ]: 4400 : auto col = w_b.col(col_id);
65 [ + + ]: 13200 : for (unsigned row_id = 0; row_id != num_rows; ++row_id)
66 [ + - + - ]: 8800 : col(row_id) = normal_distribution(g);
67 : 4400 : }
68 : 880 : return w_b;
69 [ + - ]: 880 : }
70 : :
71 : : }
72 : :
73 : :
74 : 1 : /*======================================================================================================================
75 : : * RDC
76 : : *====================================================================================================================*/
77 : :
78 : 0 : float m::rdc(const Eigen::MatrixXf &X, const Eigen::MatrixXf &Y, unsigned k, float stddev)
79 : : {
80 : 0 : M_insist(X.rows() == Y.rows());
81 : :
82 : : /* Estimate copula transformation by CDFs. */
83 : 0 : MatrixXf CDFs_of_X = create_CDF_matrix(X);
84 [ # # ]: 0 : MatrixXf CDFs_of_Y = create_CDF_matrix(Y);
85 [ # # ]: 0 : M_insist(CDFs_of_X.rows() == X.rows());
86 [ # # ]: 0 : M_insist(CDFs_of_Y.rows() == X.rows());
87 [ # # ]: 0 : M_insist(CDFs_of_X.cols() == X.cols() + 1);
88 [ # # ]: 0 : M_insist(CDFs_of_Y.cols() == Y.cols() + 1);
89 : :
90 [ # # # # ]: 0 : return rdc_precomputed_CDF(CDFs_of_X, CDFs_of_Y, k, stddev);
91 : 0 : }
92 : :
93 : : template<typename URBG>
94 : 440 : float m::rdc_precomputed_CDF(MatrixXf &CDFs_of_X, MatrixXf &CDFs_of_Y, unsigned k, float stddev, URBG &&g)
95 : : {
96 : 440 : unsigned num_rows = CDFs_of_X.rows();
97 : 440 : M_insist(num_rows == CDFs_of_Y.rows());
98 : :
99 : : /* Generate random non-linear projections. */
100 : 440 : MatrixXf X_w_b = create_w_b(stddev, CDFs_of_X.cols(), k, g);
101 [ + - ]: 440 : MatrixXf Y_w_b = create_w_b(stddev, CDFs_of_Y.cols(), k, g);
102 [ + - ]: 440 : M_insist(X_w_b.rows() == CDFs_of_X.cols());
103 [ + - ]: 440 : M_insist(Y_w_b.rows() == CDFs_of_Y.cols());
104 [ + - ]: 440 : M_insist(X_w_b.cols() == k);
105 [ + - ]: 440 : M_insist(Y_w_b.cols() == k);
106 : :
107 : : /* Compute canonical correlations. */
108 [ + - ]: 440 : MatrixXf concat_matrix(CDFs_of_X.rows(), 2 * k);
109 [ + - + - : 440 : concat_matrix << (CDFs_of_X * X_w_b), (CDFs_of_Y * Y_w_b);
+ - + - ]
110 [ + - + - : 440 : concat_matrix.array() = concat_matrix.array().sin();
+ - + - ]
111 : :
112 : : /* Center concat matrix. */
113 [ + - + - : 440 : const RowVectorXf means = concat_matrix.colwise().mean();
+ - ]
114 [ + - + - ]: 440 : concat_matrix.rowwise() -= means;
115 : :
116 : : /* Compute covariance matrix. */
117 [ + - + - : 440 : MatrixXf covariance_matrix = concat_matrix.adjoint() * concat_matrix / float(num_rows - 1);
+ - + - ]
118 : :
119 [ + - + - ]: 440 : FullPivLU<MatrixXf> CXX(covariance_matrix.block(0, 0, k, k));
120 [ + - + - ]: 440 : FullPivLU<MatrixXf> CYY(covariance_matrix.block(k, k, k, k));
121 [ + - + - ]: 440 : MatrixXf CXY = covariance_matrix.block(0, k, k, k);
122 [ + - + - ]: 440 : MatrixXf CYX = covariance_matrix.block(k, 0, k, k);
123 : :
124 [ + - + - : 440 : SelfAdjointEigenSolver<MatrixXf> eigensolver((CXX.inverse() * CXY) * (CYY.inverse() * CYX));
+ - + - +
- - + ]
125 : :
126 [ + - + - ]: 440 : return sqrtf(eigensolver.eigenvalues().maxCoeff());
127 : 440 : }
|
