1
#include "lib.h"
2
using namespace quanteda;
3
using namespace arma;
4

5

6
//find the principle elements for the sparse residual matrix
7 1
void create_residual_ca(std::size_t row_num, const arma::sp_mat& dfm, const arma::colvec &rsum, const arma::rowvec &csum,
8
                        const double residual_floor, const std::size_t K, Triplets &residual_tri)
9
{
10 1
    for (std::size_t k = 0; k < K; k++){
11 1
        double residual = (dfm(row_num, k) - rsum(row_num) * csum(k)) / sqrt(rsum(row_num) * csum(k) );
12 1
        if (fabs(residual) > residual_floor) {
13 1
            Triplet mat_triplet = std::make_tuple(row_num, k, residual);
14 1
            residual_tri.push_back(mat_triplet);
15
        }
16
    }
17
}
18

19
//Create the residual matrix
20 1
struct Res : public Worker {
21
    const arma::sp_mat &dfm;
22
    const arma::colvec &rsum;
23
    const arma::rowvec &csum;
24
    const double residual_floor;
25
    const std::size_t K;
26
    
27
    // output: residual[index0, index1, residual_value]
28
    Triplets &residual_tri;
29
    
30
    //constructor
31 1
    Res(const arma::sp_mat& dfm, const arma::colvec &rsum, const arma::rowvec &csum, 
32
        const double residual_floor, const std::size_t K, Triplets &residual_tri):
33 1
        dfm(dfm), rsum(rsum), csum(csum), residual_floor(residual_floor), K(K), residual_tri(residual_tri) {}
34
    
35 1
    void operator() (std::size_t begin, std::size_t end) {
36 1
        for (std::size_t i = begin; i < end; i++) {
37 1
            create_residual_ca(i, dfm, rsum, csum, residual_floor, K, residual_tri);
38
        }
39
    }
40
};
41

42
// [[Rcpp::export]]
43 1
S4 qatd_cpp_ca(const arma::sp_mat &dfm, const double residual_floor){
44
    
45 1
    const std::size_t N = dfm.n_rows;
46 1
    const std::size_t K = dfm.n_cols;
47
    
48
    // Construct Chi-Sq Residuals	
49 1
    const arma::colvec rsum(sum(dfm,1));
50 1
    const arma::rowvec csum(sum(dfm,0));
51
    
52
    //create the residual matrix
53 1
    Triplets residual_tri;
54 1
    residual_tri.reserve(N * K / 1000); // assume 99.9% sparsity
55

56
#if QUANTEDA_USE_TBB
57 1
    Res res(dfm, rsum, csum, residual_floor, K, residual_tri);
58 1
    parallelFor(0, N, res);
59
#else
60
    for (std::size_t i = 0; i < N; i++) {
61
        create_residual_ca(i, dfm, rsum, csum, residual_floor, K, residual_tri);
62
    }
63
#endif
64
    
65 1
    return to_matrix( residual_tri, N, K, false );
66
//     std::size_t residual_size = residual_tri.size();
67
//     IntegerVector dim_ = IntegerVector::create(N, K);
68
//     IntegerVector i_(residual_size), j_(residual_size);
69
//     NumericVector x_(residual_size);
70
//
71
//     for (std::size_t k = 0; k < residual_tri.size(); k++) {
72
//         i_[k] = std::get<0>(residual_tri[k]);
73
//         j_[k] = std::get<1>(residual_tri[k]);
74
//         x_[k] = std::get<2>(residual_tri[k]);
75
//     }
76
//
77
//     S4 res_("dgTMatrix");
78
//     res_.slot("i") = i_;
79
//     res_.slot("j") = j_;
80
//     res_.slot("x") = x_;
81
//     res_.slot("Dim") = dim_;
82
//
83
//     return res_;
84
    
85
}
86

87
/***R
88
smoke <- matrix(c(4,2,3,2, 4,5,7,4,25,10,12,4,18,24,33,13,10,6,7,2), nrow = 5, ncol = 4, byrow = T)
89
residual_floor <- 0.1
90
n = 195
91
P <- as.dfm(smoke) / n
92
qatd_cpp_ca(P, residual_floor/sqrt(n))
93
*/

Read our documentation on viewing source code .

Loading