wordfish_mt.cpp at a5c4412 ⋅ quanteda/quanteda.textmodels

#include "lib.h"
using namespace quanteda;
using namespace RcppParallel;
using namespace Rcpp;
using namespace arma;

//# define RESIDUALS_LIM 0.5
# define NUMSVD 1
# define OUTERITER 100
# define INNERITER 10
# define LASTLP -2000000000000.0
# define SA 0.99 //simulated annealing algorithm

//find the principle elements for the sparse residual matrix
void create_residual(const std::size_t row_num, const arma::sp_mat& wfm, const arma::colvec &rsum, const arma::rowvec &csum, const double &asum,
                     const double residual_floor, const std::size_t K, Triplets &residual_tri) {
    for (std::size_t k = 0; k < K; k++) {
        double residual = (wfm(row_num, k) - rsum(row_num) * csum(k) / asum) / sqrt(rsum(row_num) * csum(k) / asum);
        if (fabs(residual) > residual_floor) {
            Triplet mat_triplet = std::make_tuple(row_num, k, residual);
            residual_tri.push_back(mat_triplet);
        }
    }
}
//Create the residual matrix
struct Residual : public Worker {
    const arma::sp_mat& wfm;
    const arma::colvec &rsum;
    const arma::rowvec &csum;
    const double asum;
    const double residual_floor;
    const std::size_t K;
    
    // output: residual[index0, index1, residual_value]
    Triplets &residual_tri;
    
    //constructor
    Residual(const arma::sp_mat& wfm, const arma::colvec &rsum, const arma::rowvec &csum, const double asum,
             const double residual_floor, const std::size_t K, Triplets &residual_tri):
             wfm(wfm), rsum(rsum), csum(csum), asum(asum), residual_floor(residual_floor), 
             K(K), residual_tri(residual_tri) {}

    
    void operator() (std::size_t begin, std::size_t end) {
        for (std::size_t i = begin; i < end; i++) {
            create_residual(i, wfm, rsum, csum, asum, residual_floor, K, residual_tri);
        }
    }
};

//Compute LOG-POSTERIOR 
struct LogPos : public Worker {
    const arma::colvec& alpha; 
    const arma::rowvec& psi;
    const arma::rowvec& beta;
    const arma::colvec& theta;
    const arma::sp_mat& wfm;
    const std::size_t K;
    
    // accumulated value
    double lp;
    
    // constructors
    LogPos(const arma::colvec& alpha, const arma::rowvec& psi, const arma::rowvec& beta, 
           const arma::colvec& theta, const arma::sp_mat& wfm, const std::size_t K) :
           alpha(alpha), psi(psi), beta(beta), theta(theta), wfm(wfm), K(K), lp(0) {}
    
    LogPos(const LogPos& sum, Split) : 
           alpha(sum.alpha), psi(sum.psi), beta(sum.beta),
           theta(sum.theta), wfm(sum.wfm), K(sum.K), lp(0) {}
    
    
    void operator() (std::size_t begin, std::size_t end) {
        for (std::size_t i = begin; i < end; i++) {
            double temp = lp;
            for (std::size_t k = 0; k < K; k++){
                double loglambdaik = alpha(i) + psi(k) + beta(k) * theta(i);
                temp += loglambdaik * wfm(i,k) - exp(loglambdaik);
            }
            lp = temp;
        }
    }
    
    // join lp with that of another LogPos
    void join(const LogPos& rhs) { 
        lp += rhs.lp; 
    }
};

// UPDATE WORD PARAMETERS
struct WordPar : public Worker {
    const arma::colvec& alpha; 
    RVector<double> psi; 
    RVector<double> beta; 
    const arma::colvec& theta; 
    const arma::rowvec& phi;
    const arma::sp_mat& wfm;
    const NumericVector& tolvec;
    const int outeriter;
    const double priorprecpsi;
    const double priorprecbeta;
    double stepsize;
    
    WordPar(const arma::colvec& alpha,  NumericVector& psi,  NumericVector& beta, 
            const arma::colvec& theta, const arma::rowvec& phi, const arma::sp_mat& wfm,
            const NumericVector& tolvec, const int outeriter, 
            const double priorprecpsi, const double priorprecbeta, double stepsize) :
            alpha(alpha), psi(psi), beta(beta), theta(theta), phi(phi), wfm(wfm), tolvec(tolvec), 
            outeriter(outeriter), priorprecpsi(priorprecpsi), priorprecbeta(priorprecbeta), stepsize(stepsize) {}
    
    void operator() (std::size_t begin, std::size_t end) {
        for (std::size_t k = begin; k < end; k++) {
            double cc = 1;
            int inneriter = 0;
            if (outeriter == 1) stepsize = 0.5;
            arma::colvec col_k(wfm.col(k));
            arma::colvec lambdak(alpha.size());
            arma::mat pars(2,1);
            arma::mat newpars(2,1);
            arma::mat G(2,1);
            arma::mat H(2,2);
            while ((cc > tolvec(1)) && inneriter < INNERITER){
                inneriter++;
                lambdak = exp(alpha + psi[k] + beta[k] * theta);
                G(0,0) = accu(col_k - lambdak) / phi(k) - psi[k] * priorprecpsi;
                G(1,0) = accu(theta % (col_k - lambdak)) / phi(k) - beta[k] * priorprecbeta;
                H(0,0) = -accu(lambdak) / phi(k) - priorprecpsi;
                H(1,0) = -accu(theta % lambdak) / phi(k);
                H(0,1) = H(1,0);
                H(1,1) = -accu((theta % theta) % lambdak) / phi(k) - priorprecbeta;
                pars(0,0) = psi[k];
                pars(1,0) = beta[k];
                newpars(0,0) = pars(0,0) - stepsize * (H(1,1) * G(0,0) - H(0,1) * G(1,0)) / (H(0,0) * H(1,1) - H(0,1) * H(1,0));
                newpars(1,0) = pars(1,0) - stepsize * (H(0,0) * G(1,0) - H(1,0) * G(0,0)) / (H(0,0) * H(1,1) - H(0,1) * H(1,0));
                newpars(0,0) *= SA;  //simulated annealing algorithm
                newpars(1,0) *= SA;
                psi[k] = newpars(0,0);
                beta[k] = newpars(1,0);
                cc = abs(newpars - pars).max();
                stepsize = 1.0;
            }
        }
    } //End of operator
};

// UPDATE DOCUMENT PARAMETERS
struct DocPar : public Worker {
    RVector<double> alpha; 
    const arma::rowvec& psi;
    const arma::rowvec& beta; 
    RVector<double> theta; 
    const arma::rowvec& phi;
    const arma::sp_mat& wfm;
    const NumericVector& tolvec;
    const int outeriter;
    const double priorprecalpha;
    const double priorprectheta;
    double stepsize;
    
    DocPar(NumericVector& alpha,  const arma::rowvec& psi,  const arma::rowvec& beta, 
           NumericVector& theta, const arma::rowvec& phi, const arma::sp_mat& wfm,
           const NumericVector& tolvec, const int outeriter, 
           const double priorprecalpha, const double priorprectheta, double stepsize) :
           alpha(alpha), psi(psi), beta(beta), theta(theta), phi(phi), wfm(wfm), tolvec(tolvec), 
           outeriter(outeriter), priorprecalpha(priorprecalpha), priorprectheta(priorprectheta), stepsize(stepsize) {}
    
    void operator() (std::size_t begin, std::size_t end) {
        for (std::size_t i = begin; i < end; i++) {
            double cc = 1;
            int inneriter = 0;
            if (outeriter == 1) stepsize = 0.5;
            arma::rowvec row_i(wfm.row(i));
            arma::rowvec lambdai(alpha.size());
            arma::mat pars(2,1);
            arma::mat newpars(2,1);
            arma::mat G(2,1);
            arma::mat H(2,2);
            while ((cc > tolvec(1)) && inneriter < INNERITER){
                inneriter++;
                lambdai = exp(alpha[i] + psi + beta * theta[i]);
                G(0,0) = accu((row_i - lambdai) / phi) - alpha[i] * priorprecalpha;
                G(1,0) = accu((beta % (row_i - lambdai)) / phi) - theta[i] * priorprectheta;		
                H(0,0) = -accu(lambdai / phi) - priorprecalpha;
                H(1,0) = -accu((beta % lambdai) / phi);
                H(0,1) = H(1,0);
                H(1,1) = -accu(((beta % beta) % lambdai) / phi) - priorprectheta;
                pars(0,0) = alpha[i];
                pars(1,0) = theta[i];
                newpars(0,0) = pars(0,0) - stepsize*(H(1,1) * G(0,0) - H(0,1) * G(1,0)) / (H(0,0) * H(1,1) - H(0,1) * H(1,0));
                newpars(1,0) = pars(1,0) - stepsize*(H(0,0) * G(1,0) - H(1,0) * G(0,0)) / (H(0,0) * H(1,1) - H(0,1) * H(1,0));
                newpars(0,0) *= SA;
                newpars(0,0) *= SA;
                alpha[i] = newpars(0,0);
                theta[i] = newpars(1,0);
                cc = abs(newpars - pars).max();	
                stepsize = 1.0;
            }
        }
    } //End of operator
};

// Update dispersion parameters -- single dispersion parameter for all words
struct DispPar : public Worker {
    const arma::colvec& alpha; 
    const arma::rowvec& psi;
    const arma::rowvec& beta;
    const arma::colvec& theta;
    const arma::sp_mat& wfm;
    const std::size_t& N;
    
    // accumulated value
    double phitmp;
    
    // constructors
    DispPar(const arma::colvec& alpha, const arma::rowvec& psi, const arma::rowvec& beta, 
            const arma::colvec& theta, const arma::sp_mat& wfm, const std::size_t& N) :
            alpha(alpha), psi(psi), beta(beta), theta(theta), wfm(wfm), N(N), phitmp(0) {}
    
    DispPar(const DispPar& sum, Split) : 
            alpha(sum.alpha), psi(sum.psi), beta(sum.beta),
            theta(sum.theta), wfm(sum.wfm), N(sum.N), phitmp(0) {}
    
    
    void operator() (std::size_t begin, std::size_t end) {
        for (std::size_t k = begin; k < end; k++) {
            float temp = phitmp;
            for (std::size_t i = 0; i < N; i++){
                double mutmp = exp(alpha(i) + psi(k) + beta(k) * theta(i));
                temp += (wfm(i,k) - mutmp) * (wfm(i,k) - mutmp) / mutmp;
            }
            phitmp = temp;
        }
    }
    
    // join phitmp with that of another DispPar
    void join(const DispPar& rhs) { 
        phitmp += rhs.phitmp; 
    }
};

// Update dispersion parameters -- individual dispersion parameter for each word
struct DispPar2 : public Worker {
    const arma::colvec& alpha; 
    const arma::rowvec& psi;
    const arma::rowvec& beta;
    const arma::colvec& theta;
    const arma::sp_mat& wfm;
    const IntegerVector& disptype;
    const NumericVector& dispmin;
    const std::size_t& N;
    const std::size_t& K;
    // output vector
    RVector<double> phi;
    
    // constructors
    DispPar2(const arma::colvec& alpha, const arma::rowvec& psi, const arma::rowvec& beta, 
             const arma::colvec& theta, const arma::sp_mat& wfm, const IntegerVector& disptype,
             const NumericVector& dispmin, const std::size_t& N, const std::size_t& K, NumericVector& phi) :
             alpha(alpha), psi(psi), beta(beta), theta(theta), wfm(wfm), disptype(disptype), 
             dispmin(dispmin), N(N), K(K), phi(phi){}
    
    void operator() (std::size_t begin, std::size_t end) {
        for (std::size_t k = begin; k < end; k++) {
            double phitmp = 0.0;
            for (std::size_t i = 0; i < N; i++){
                double mutmp = exp(alpha(i) + psi(k) + beta(k) * theta(i));
                phitmp += (wfm(i,k) - mutmp) * (wfm(i,k) - mutmp) / mutmp;
            }
            phitmp = ((K)*phitmp)/(N * K - 2 * N - 2 * K);
            if (disptype(0) == 4) {
                phi[k] = fmax(dispmin(0), phitmp);
            }else{
                phi[k] = phitmp;
            }
        }
    }
};

// Calculate Document Standard Errors
struct DocErr : public Worker {
    const arma::colvec& alpha; 
    const arma::rowvec& psi;
    const arma::rowvec& beta;
    const arma::colvec& theta;
    const arma::rowvec& phi;
    const double& priorprecalpha;
    const double& priorprectheta;
    // output vector
    RVector<double> thetaSE;
    
    // constructors
    DocErr(const arma::colvec& alpha, const arma::rowvec& psi, const arma::rowvec& beta, 
           const arma::colvec& theta, const arma::rowvec& phi, 
           const double& priorprecalpha, const double& priorprectheta, NumericVector& thetaSE) :
           alpha(alpha), psi(psi), beta(beta), theta(theta), phi(phi), 
           priorprecalpha(priorprecalpha), priorprectheta(priorprectheta), thetaSE(thetaSE){}
    
    void operator() (std::size_t begin, std::size_t end) {
        arma::rowvec lambdai(alpha.size());
        for (std::size_t i = begin; i < end; i++) {
            lambdai = exp(alpha(i) + psi + beta * theta(i));
            arma::mat H(2,2);
            H(0,0) = -accu(lambdai / phi) - priorprecalpha;
            H(1,0) = -accu((beta % lambdai) / phi);
            H(0,1) = H(1,0);
            H(1,1) = -accu(((beta % beta) % lambdai) / phi) - priorprectheta;
            thetaSE[i] = sqrt(-1.0 * H(0,0) / (H(0,0) * H(1,1) - H(1,0) * H(0,1)));
        }
    }
};
// [[Rcpp::export]]

Rcpp::List qatd_cpp_wordfish(arma::sp_mat &wfm, 
                             IntegerVector& dirvec, 
                             NumericVector& priorvec, 
                             NumericVector& tolvec, 
                             IntegerVector& disptype, 
                             NumericVector& dispmin, 
                             bool ABS,bool svd_sparse, 
                             double residual_floor){
    
    // DEFINE INPUTS
    double priorprecalpha = priorvec(0);
    double priorprecpsi = priorvec(1);
    double priorprecbeta = priorvec(2);
    double priorprectheta = priorvec(3);		
    
    // std::random_device rd; // random engine
    // std::mt19937 mt(rd());
    std::mt19937 mt(time(0)); // issue #1063
    std::uniform_real_distribution<double> dist(0.0, 1.0);
    
    std::size_t N = wfm.n_rows;
    std::size_t K = wfm.n_cols;
    
    // Set initial values
    arma::colvec alpha(N); 
    arma::rowvec psi(K); 
    arma::rowvec beta(K); 
    arma::colvec theta(N); 
    arma::colvec thetaSE(N); // document position standard errors
    arma::rowvec phi(K); // word-level dispersion parameters
    phi.fill(1.0);
 
    // Construct Chi-Sq Residuals	
    arma::colvec rsum(sum(wfm,1));
    arma::rowvec csum(sum(wfm,0));
    double asum = accu(wfm);
    if (svd_sparse == true){
     
        //create the residual matrix
        Triplets residual_tri;

#if QUANTEDA_USE_TBB
        Residual residual(wfm, rsum, csum, asum, residual_floor, K, residual_tri);
        parallelFor(0, N, residual);
        //Rcout<<"used TBB"<<std::endl;
#else        
        for (std::size_t i = 0; i < N; i++) {
            create_residual(i, wfm, rsum, csum, asum, residual_floor, K, residual_tri);
        }
        //Rcout<<"not use TBB"<<std::endl;
#endif
        // Convert to Rcpp objects
        std::size_t mat_size = residual_tri.size();
        arma::umat index_mat(2, mat_size, arma::fill::zeros);
        arma::vec w_values(mat_size, arma::fill::zeros);
        for (std::size_t k = 0; k < residual_tri.size(); k++) {
            index_mat(0,k) = std::get<0>(residual_tri[k]);
            index_mat(1,k) = std::get<1>(residual_tri[k]);
            w_values(k) = std::get<2>(residual_tri[k]);
        }
        
        // constract the sparse matrix
        arma::sp_mat C(index_mat, w_values, N, K);

        const int svdk = NUMSVD;
        arma::mat U(N, svdk);
        arma::vec s(svdk);
        arma::mat V(K, svdk);
        arma::svds(U, s, V, C, svdk);
        for (std::size_t i = 0; i < N; i++) theta(i) = pow(rsum(i) / asum, -0.5) * U(i, 0);
    } else {
        // Load initial values
        for (std::size_t i=0; i < N; i++) theta(i) = pow(rsum(i) / asum, -0.5) - dist(mt);//* U(i, 0);
    }
    for (std::size_t k = 0; k < K; k++) beta(k) = 0;//  pow(csum(k)/asum,-0.5) * V(k,0);
    //beta.fill(0.0);
    alpha = log(rsum);
    psi = log(csum/N);
    
    alpha = alpha - log(mean(rsum));
    theta = (theta - mean(theta)) / stddev(theta);  
    
    // Create temporary variables
    arma::mat pars(2,1);
    arma::mat newpars(2,1);		
    arma::mat G(2,1);
    arma::mat H(2,2);
    arma::rowvec lambdai(K);
    arma::colvec lambdak(N);
    double stepsize = 1.0;
    int outeriter = 0;

    //Initialize LOG-POSTERIOR 
    double lastlp = LASTLP;
    double lp = -1.0 * (accu(0.5 * ((alpha % alpha) * priorprecalpha)) + accu(0.5 * ((psi  % psi) * priorprecpsi)) 
                          + accu(0.5 * ((beta  % beta) * priorprecbeta)) + accu(0.5 * ((theta % theta) * priorprectheta)));
    
    // compute LOG posterior
    LogPos logPos(alpha, psi, beta, theta, wfm, K);
    parallelReduce(0, N, logPos);
    lp += logPos.lp;
    
    
    // BEGIN WHILE LOOP
    while (((lp - lastlp) > tolvec(0)) && outeriter < OUTERITER) {	
        outeriter++;
        //Rcout<<"alphs_b="<<mean(alpha)<<" theta_B="<<mean(theta)<<std::endl;
        
        // UPDATE WORD PARAMETERS
        NumericVector psi_N(psi.begin(), psi.end());
        NumericVector beta_N(beta.begin(), beta.end());
        WordPar wordPar(alpha, psi_N, beta_N, theta, phi, wfm, tolvec,
                        outeriter, priorprecpsi, priorprecbeta, stepsize);
        parallelFor(0, K, wordPar);
        
        psi = as<arma::rowvec>(psi_N);
        beta = as<arma::rowvec>(beta_N);

        // UPDATE DOCUMENT PARAMETERS
        NumericVector alpha_N(alpha.begin(), alpha.end());
        NumericVector theta_N(theta.begin(), theta.end());
        DocPar docPar(alpha_N, psi, beta, theta_N, phi, wfm, tolvec,
                        outeriter, priorprecalpha, priorprectheta, stepsize);
        parallelFor(0, N, docPar);
        
        alpha = as<arma::colvec>(alpha_N);
        theta = as<arma::colvec>(theta_N);
        
        
        // UPDATE DISPERSION PARAMETERS	  
        
        if (disptype(0) == 2) { // single dispersion parameter for all words
            DispPar dispPar(alpha, psi, beta, theta, wfm, N);
            parallelReduce(0, K, dispPar);
            double phitmp = dispPar.phitmp;
            phitmp /= N * K - 2 * N - 2 * K;
            phi.fill(phitmp);
        }
        
        if (disptype(0) >= 3) { // individual dispersion parameter for each word
            NumericVector phi_N(phi.begin(), phi.end());
            DispPar2 dispPar2(alpha, psi, beta, theta, wfm, disptype, dispmin, N, K, phi_N);
            parallelFor(0, K, dispPar2);
            phi = as<arma::rowvec>(phi_N);
        }	  			
        
        alpha = alpha - mean(alpha);
        theta = (theta - mean(theta))/stddev(theta);		

        // CHECK LOG-POSTERIOR FOR CONVERGENCE
        lastlp = lp;
        lp = -1.0 * (accu(0.5 * ((alpha % alpha) * priorprecalpha)) + accu(0.5 * ((psi % psi) * priorprecpsi)) +
                     accu(0.5 * ((beta % beta) * priorprecbeta)) + accu(0.5 * ((theta % theta) * priorprectheta)));

        LogPos logPos2(alpha, psi, beta, theta, wfm, K);
        parallelReduce(0, N, logPos2);
        lp += logPos2.lp;
        // END WHILE LOOP		
    } 
    
    // Fix Global Polarity  
    //Rcout<<"outeriter="<<outeriter<<"  lp - lastlp= "<<lp - lastlp<<std::endl;
    // added the -1 because C counts from ZERO...  -- KB
    if (theta(dirvec(0) - 1) > theta(dirvec(1) - 1)) {
        beta = -beta;
        theta = -theta;
    }
    
    // COMPUTE DOCUMENT STANDARD ERRORS
    NumericVector thetaSE_N(thetaSE.begin(), thetaSE.end());
    DocErr docErr(alpha, psi, beta, theta, phi, priorprecalpha, priorprectheta, thetaSE_N);
    parallelFor(0, N, docErr);
    thetaSE = as<arma::colvec>(thetaSE_N);
    
    // DEFINE OUTPUT	
    return Rcpp::List::create(Rcpp::Named("theta") = theta,
                              Rcpp::Named("alpha") = alpha,
                              Rcpp::Named("psi") = psi,
                              Rcpp::Named("beta") = beta,
                              Rcpp::Named("phi") = phi,
                              Rcpp::Named("thetaSE") = thetaSE);
    
}

Read our documentation on viewing source code .

Navigation	Overlay
`t` Navigate files	`h` Toggle hits
`y` Change url to tip of branch	`m` Toggle misses
`b / v` Jump to prev/next hit line	`p` Toggle partial
`z / x` Jump to prev/next missed or partial line	`1..9` Toggle flags
`shift + o` Open current page in GitHub	`a` Toggle all on
`/ or ?` Show keyboard shortcuts dialog	`c` Toggle context lines or commits

1		#include "lib.h"
2		using namespace quanteda;
3		using namespace RcppParallel;
4		using namespace Rcpp;
5		using namespace arma;
6
7		//# define RESIDUALS_LIM 0.5
8		# define NUMSVD 1
9		# define OUTERITER 100
10		# define INNERITER 10
11		# define LASTLP -2000000000000.0
12		# define SA 0.99 //simulated annealing algorithm
13
14		//find the principle elements for the sparse residual matrix
15	1	void create_residual(const std::size_t row_num, const arma::sp_mat& wfm, const arma::colvec &rsum, const arma::rowvec &csum, const double &asum,
16		const double residual_floor, const std::size_t K, Triplets &residual_tri) {
17	1	for (std::size_t k = 0; k < K; k++) {
18	1	double residual = (wfm(row_num, k) - rsum(row_num) * csum(k) / asum) / sqrt(rsum(row_num) * csum(k) / asum);
19	1	if (fabs(residual) > residual_floor) {
20	1	Triplet mat_triplet = std::make_tuple(row_num, k, residual);
21	1	residual_tri.push_back(mat_triplet);
22		}
23		}
24		}
25		//Create the residual matrix
26	1	struct Residual : public Worker {
27		const arma::sp_mat& wfm;
28		const arma::colvec &rsum;
29		const arma::rowvec &csum;
30		const double asum;
31		const double residual_floor;
32		const std::size_t K;
33
34		// output: residual[index0, index1, residual_value]
35		Triplets &residual_tri;
36
37		//constructor
38	1	Residual(const arma::sp_mat& wfm, const arma::colvec &rsum, const arma::rowvec &csum, const double asum,
39		const double residual_floor, const std::size_t K, Triplets &residual_tri):
40	1	wfm(wfm), rsum(rsum), csum(csum), asum(asum), residual_floor(residual_floor),
41	1	K(K), residual_tri(residual_tri) {}
42
43
44	1	void operator() (std::size_t begin, std::size_t end) {
45	1	for (std::size_t i = begin; i < end; i++) {
46	1	create_residual(i, wfm, rsum, csum, asum, residual_floor, K, residual_tri);
47		}
48		}
49		};
50
51		//Compute LOG-POSTERIOR
52	1	struct LogPos : public Worker {
53		const arma::colvec& alpha;
54		const arma::rowvec& psi;
55		const arma::rowvec& beta;
56		const arma::colvec& theta;
57		const arma::sp_mat& wfm;
58		const std::size_t K;
59
60		// accumulated value
61		double lp;
62
63		// constructors
64	1	LogPos(const arma::colvec& alpha, const arma::rowvec& psi, const arma::rowvec& beta,
65		const arma::colvec& theta, const arma::sp_mat& wfm, const std::size_t K) :
66	1	alpha(alpha), psi(psi), beta(beta), theta(theta), wfm(wfm), K(K), lp(0) {}
67
68	1	LogPos(const LogPos& sum, Split) :
69	1	alpha(sum.alpha), psi(sum.psi), beta(sum.beta),
70	1	theta(sum.theta), wfm(sum.wfm), K(sum.K), lp(0) {}
71
72
73	1	void operator() (std::size_t begin, std::size_t end) {
74	1	for (std::size_t i = begin; i < end; i++) {
75	1	double temp = lp;
76	1	for (std::size_t k = 0; k < K; k++){
77	1	double loglambdaik = alpha(i) + psi(k) + beta(k) * theta(i);
78	1	temp += loglambdaik * wfm(i,k) - exp(loglambdaik);
79		}
80	1	lp = temp;
81		}
82		}
83
84		// join lp with that of another LogPos
85	1	void join(const LogPos& rhs) {
86	1	lp += rhs.lp;
87		}
88		};
89
90		// UPDATE WORD PARAMETERS
91	1	struct WordPar : public Worker {
92		const arma::colvec& alpha;
93		RVector<double> psi;
94		RVector<double> beta;
95		const arma::colvec& theta;
96		const arma::rowvec& phi;
97		const arma::sp_mat& wfm;
98		const NumericVector& tolvec;
99		const int outeriter;
100		const double priorprecpsi;
101		const double priorprecbeta;
102		double stepsize;
103
104	1	WordPar(const arma::colvec& alpha, NumericVector& psi, NumericVector& beta,
105		const arma::colvec& theta, const arma::rowvec& phi, const arma::sp_mat& wfm,
106		const NumericVector& tolvec, const int outeriter,
107		const double priorprecpsi, const double priorprecbeta, double stepsize) :
108	1	alpha(alpha), psi(psi), beta(beta), theta(theta), phi(phi), wfm(wfm), tolvec(tolvec),
109	1	outeriter(outeriter), priorprecpsi(priorprecpsi), priorprecbeta(priorprecbeta), stepsize(stepsize) {}
110
111	1	void operator() (std::size_t begin, std::size_t end) {
112	1	for (std::size_t k = begin; k < end; k++) {
113	1	double cc = 1;
114	1	int inneriter = 0;
115	1	if (outeriter == 1) stepsize = 0.5;
116	1	arma::colvec col_k(wfm.col(k));
117	1	arma::colvec lambdak(alpha.size());
118	1	arma::mat pars(2,1);
119	1	arma::mat newpars(2,1);
120	1	arma::mat G(2,1);
121	1	arma::mat H(2,2);
122	1	while ((cc > tolvec(1)) && inneriter < INNERITER){
123	1	inneriter++;
124	1	lambdak = exp(alpha + psi[k] + beta[k] * theta);
125	1	G(0,0) = accu(col_k - lambdak) / phi(k) - psi[k] * priorprecpsi;
126	1	G(1,0) = accu(theta % (col_k - lambdak)) / phi(k) - beta[k] * priorprecbeta;
127	1	H(0,0) = -accu(lambdak) / phi(k) - priorprecpsi;
128	1	H(1,0) = -accu(theta % lambdak) / phi(k);
129	1	H(0,1) = H(1,0);
130	1	H(1,1) = -accu((theta % theta) % lambdak) / phi(k) - priorprecbeta;
131	1	pars(0,0) = psi[k];
132	1	pars(1,0) = beta[k];
133	1	newpars(0,0) = pars(0,0) - stepsize * (H(1,1) * G(0,0) - H(0,1) * G(1,0)) / (H(0,0) * H(1,1) - H(0,1) * H(1,0));
134	1	newpars(1,0) = pars(1,0) - stepsize * (H(0,0) * G(1,0) - H(1,0) * G(0,0)) / (H(0,0) * H(1,1) - H(0,1) * H(1,0));
135	1	newpars(0,0) *= SA; //simulated annealing algorithm
136	1	newpars(1,0) *= SA;
137	1	psi[k] = newpars(0,0);
138	1	beta[k] = newpars(1,0);
139	1	cc = abs(newpars - pars).max();
140	1	stepsize = 1.0;
141		}
142		}
143	1	} //End of operator
144		};
145
146		// UPDATE DOCUMENT PARAMETERS
147	1	struct DocPar : public Worker {
148		RVector<double> alpha;
149		const arma::rowvec& psi;
150		const arma::rowvec& beta;
151		RVector<double> theta;
152		const arma::rowvec& phi;
153		const arma::sp_mat& wfm;
154		const NumericVector& tolvec;
155		const int outeriter;
156		const double priorprecalpha;
157		const double priorprectheta;
158		double stepsize;
159
160	1	DocPar(NumericVector& alpha, const arma::rowvec& psi, const arma::rowvec& beta,
161		NumericVector& theta, const arma::rowvec& phi, const arma::sp_mat& wfm,
162		const NumericVector& tolvec, const int outeriter,
163		const double priorprecalpha, const double priorprectheta, double stepsize) :
164	1	alpha(alpha), psi(psi), beta(beta), theta(theta), phi(phi), wfm(wfm), tolvec(tolvec),
165	1	outeriter(outeriter), priorprecalpha(priorprecalpha), priorprectheta(priorprectheta), stepsize(stepsize) {}
166
167	1	void operator() (std::size_t begin, std::size_t end) {
168	1	for (std::size_t i = begin; i < end; i++) {
169	1	double cc = 1;
170	1	int inneriter = 0;
171	1	if (outeriter == 1) stepsize = 0.5;
172	1	arma::rowvec row_i(wfm.row(i));
173	1	arma::rowvec lambdai(alpha.size());
174	1	arma::mat pars(2,1);
175	1	arma::mat newpars(2,1);
176	1	arma::mat G(2,1);
177	1	arma::mat H(2,2);
178	1	while ((cc > tolvec(1)) && inneriter < INNERITER){
179	1	inneriter++;
180	1	lambdai = exp(alpha[i] + psi + beta * theta[i]);
181	1	G(0,0) = accu((row_i - lambdai) / phi) - alpha[i] * priorprecalpha;
182	1	G(1,0) = accu((beta % (row_i - lambdai)) / phi) - theta[i] * priorprectheta;
183	1	H(0,0) = -accu(lambdai / phi) - priorprecalpha;
184	1	H(1,0) = -accu((beta % lambdai) / phi);
185	1	H(0,1) = H(1,0);
186	1	H(1,1) = -accu(((beta % beta) % lambdai) / phi) - priorprectheta;
187	1	pars(0,0) = alpha[i];
188	1	pars(1,0) = theta[i];
189	1	newpars(0,0) = pars(0,0) - stepsize(H(1,1) G(0,0) - H(0,1) * G(1,0)) / (H(0,0) * H(1,1) - H(0,1) * H(1,0));
190	1	newpars(1,0) = pars(1,0) - stepsize(H(0,0) G(1,0) - H(1,0) * G(0,0)) / (H(0,0) * H(1,1) - H(0,1) * H(1,0));
191	1	newpars(0,0) *= SA;
192	1	newpars(0,0) *= SA;
193	1	alpha[i] = newpars(0,0);
194	1	theta[i] = newpars(1,0);
195	1	cc = abs(newpars - pars).max();
196	1	stepsize = 1.0;
197		}
198		}
199	1	} //End of operator
200		};
201
202		// Update dispersion parameters -- single dispersion parameter for all words
203	0	struct DispPar : public Worker {
204		const arma::colvec& alpha;
205		const arma::rowvec& psi;
206		const arma::rowvec& beta;
207		const arma::colvec& theta;
208		const arma::sp_mat& wfm;
209		const std::size_t& N;
210
211		// accumulated value
212		double phitmp;
213
214		// constructors
215	0	DispPar(const arma::colvec& alpha, const arma::rowvec& psi, const arma::rowvec& beta,
216		const arma::colvec& theta, const arma::sp_mat& wfm, const std::size_t& N) :
217	0	alpha(alpha), psi(psi), beta(beta), theta(theta), wfm(wfm), N(N), phitmp(0) {}
218
219	0	DispPar(const DispPar& sum, Split) :
220	0	alpha(sum.alpha), psi(sum.psi), beta(sum.beta),
221	0	theta(sum.theta), wfm(sum.wfm), N(sum.N), phitmp(0) {}
222
223
224	0	void operator() (std::size_t begin, std::size_t end) {
225	0	for (std::size_t k = begin; k < end; k++) {
226	0	float temp = phitmp;
227	0	for (std::size_t i = 0; i < N; i++){
228	0	double mutmp = exp(alpha(i) + psi(k) + beta(k) * theta(i));
229	0	temp += (wfm(i,k) - mutmp) * (wfm(i,k) - mutmp) / mutmp;
230		}
231	0	phitmp = temp;
232		}
233		}
234
235		// join phitmp with that of another DispPar
236	0	void join(const DispPar& rhs) {
237	0	phitmp += rhs.phitmp;
238		}
239		};
240
241		// Update dispersion parameters -- individual dispersion parameter for each word
242	1	struct DispPar2 : public Worker {
243		const arma::colvec& alpha;
244		const arma::rowvec& psi;
245		const arma::rowvec& beta;
246		const arma::colvec& theta;
247		const arma::sp_mat& wfm;
248		const IntegerVector& disptype;
249		const NumericVector& dispmin;
250		const std::size_t& N;
251		const std::size_t& K;
252		// output vector
253		RVector<double> phi;
254
255		// constructors
256	1	DispPar2(const arma::colvec& alpha, const arma::rowvec& psi, const arma::rowvec& beta,
257		const arma::colvec& theta, const arma::sp_mat& wfm, const IntegerVector& disptype,
258		const NumericVector& dispmin, const std::size_t& N, const std::size_t& K, NumericVector& phi) :
259	1	alpha(alpha), psi(psi), beta(beta), theta(theta), wfm(wfm), disptype(disptype),
260	1	dispmin(dispmin), N(N), K(K), phi(phi){}
261
262	1	void operator() (std::size_t begin, std::size_t end) {
263	1	for (std::size_t k = begin; k < end; k++) {
264	1	double phitmp = 0.0;
265	1	for (std::size_t i = 0; i < N; i++){
266	1	double mutmp = exp(alpha(i) + psi(k) + beta(k) * theta(i));
267	1	phitmp += (wfm(i,k) - mutmp) * (wfm(i,k) - mutmp) / mutmp;
268		}
269	1	phitmp = ((K)phitmp)/(N K - 2 * N - 2 * K);
270	1	if (disptype(0) == 4) {
271	0	phi[k] = fmax(dispmin(0), phitmp);
272		}else{
273	1	phi[k] = phitmp;
274		}
275		}
276		}
277		};
278
279		// Calculate Document Standard Errors
280	1	struct DocErr : public Worker {
281		const arma::colvec& alpha;
282		const arma::rowvec& psi;
283		const arma::rowvec& beta;
284		const arma::colvec& theta;
285		const arma::rowvec& phi;
286		const double& priorprecalpha;
287		const double& priorprectheta;
288		// output vector
289		RVector<double> thetaSE;
290
291		// constructors
292	1	DocErr(const arma::colvec& alpha, const arma::rowvec& psi, const arma::rowvec& beta,
293		const arma::colvec& theta, const arma::rowvec& phi,
294		const double& priorprecalpha, const double& priorprectheta, NumericVector& thetaSE) :
295	1	alpha(alpha), psi(psi), beta(beta), theta(theta), phi(phi),
296	1	priorprecalpha(priorprecalpha), priorprectheta(priorprectheta), thetaSE(thetaSE){}
297
298	1	void operator() (std::size_t begin, std::size_t end) {
299	1	arma::rowvec lambdai(alpha.size());
300	1	for (std::size_t i = begin; i < end; i++) {
301	1	lambdai = exp(alpha(i) + psi + beta * theta(i));
302	1	arma::mat H(2,2);
303	1	H(0,0) = -accu(lambdai / phi) - priorprecalpha;
304	1	H(1,0) = -accu((beta % lambdai) / phi);
305	1	H(0,1) = H(1,0);
306	1	H(1,1) = -accu(((beta % beta) % lambdai) / phi) - priorprectheta;
307	1	thetaSE[i] = sqrt(-1.0 * H(0,0) / (H(0,0) * H(1,1) - H(1,0) * H(0,1)));
308		}
309		}
310		};
311		// [[Rcpp::export]]
312
313	1	Rcpp::List qatd_cpp_wordfish(arma::sp_mat &wfm,
314		IntegerVector& dirvec,
315		NumericVector& priorvec,
316		NumericVector& tolvec,
317		IntegerVector& disptype,
318		NumericVector& dispmin,
319		bool ABS,bool svd_sparse,
320		double residual_floor){
321
322		// DEFINE INPUTS
323	1	double priorprecalpha = priorvec(0);
324	1	double priorprecpsi = priorvec(1);
325	1	double priorprecbeta = priorvec(2);
326	1	double priorprectheta = priorvec(3);
327
328		// std::random_device rd; // random engine
329		// std::mt19937 mt(rd());
330	1	std::mt19937 mt(time(0)); // issue #1063
331	1	std::uniform_real_distribution<double> dist(0.0, 1.0);
332
333	1	std::size_t N = wfm.n_rows;
334	1	std::size_t K = wfm.n_cols;
335
336		// Set initial values
337	1	arma::colvec alpha(N);
338	1	arma::rowvec psi(K);
339	1	arma::rowvec beta(K);
340	1	arma::colvec theta(N);
341	1	arma::colvec thetaSE(N); // document position standard errors
342	1	arma::rowvec phi(K); // word-level dispersion parameters
343	1	phi.fill(1.0);
344
345		// Construct Chi-Sq Residuals
346	1	arma::colvec rsum(sum(wfm,1));
347	1	arma::rowvec csum(sum(wfm,0));
348	1	double asum = accu(wfm);
349	1	if (svd_sparse == true){
350
351		//create the residual matrix
352	1	Triplets residual_tri;
353
354		#if QUANTEDA_USE_TBB
355	1	Residual residual(wfm, rsum, csum, asum, residual_floor, K, residual_tri);
356	1	parallelFor(0, N, residual);
357		//Rcout<<"used TBB"<<std::endl;
358		#else
359		for (std::size_t i = 0; i < N; i++) {
360		create_residual(i, wfm, rsum, csum, asum, residual_floor, K, residual_tri);
361		}
362		//Rcout<<"not use TBB"<<std::endl;
363		#endif
364		// Convert to Rcpp objects
365	1	std::size_t mat_size = residual_tri.size();
366	1	arma::umat index_mat(2, mat_size, arma::fill::zeros);
367	1	arma::vec w_values(mat_size, arma::fill::zeros);
368	1	for (std::size_t k = 0; k < residual_tri.size(); k++) {
369	1	index_mat(0,k) = std::get<0>(residual_tri[k]);
370	1	index_mat(1,k) = std::get<1>(residual_tri[k]);
371	1	w_values(k) = std::get<2>(residual_tri[k]);
372		}
373
374		// constract the sparse matrix
375	1	arma::sp_mat C(index_mat, w_values, N, K);
376
377	1	const int svdk = NUMSVD;
378	1	arma::mat U(N, svdk);
379	1	arma::vec s(svdk);
380	1	arma::mat V(K, svdk);
381	1	arma::svds(U, s, V, C, svdk);
382	1	for (std::size_t i = 0; i < N; i++) theta(i) = pow(rsum(i) / asum, -0.5) * U(i, 0);
383		} else {
384		// Load initial values
385	0	for (std::size_t i=0; i < N; i++) theta(i) = pow(rsum(i) / asum, -0.5) - dist(mt);//* U(i, 0);
386		}
387	1	for (std::size_t k = 0; k < K; k++) beta(k) = 0;// pow(csum(k)/asum,-0.5) * V(k,0);
388		//beta.fill(0.0);
389	1	alpha = log(rsum);
390	1	psi = log(csum/N);
391
392	1	alpha = alpha - log(mean(rsum));
393	1	theta = (theta - mean(theta)) / stddev(theta);
394
395		// Create temporary variables
396	1	arma::mat pars(2,1);
397	1	arma::mat newpars(2,1);
398	1	arma::mat G(2,1);
399	1	arma::mat H(2,2);
400	1	arma::rowvec lambdai(K);
401	1	arma::colvec lambdak(N);
402	1	double stepsize = 1.0;
403	1	int outeriter = 0;
404
405		//Initialize LOG-POSTERIOR
406	1	double lastlp = LASTLP;
407	1	double lp = -1.0 * (accu(0.5 * ((alpha % alpha) * priorprecalpha)) + accu(0.5 * ((psi % psi) * priorprecpsi))
408	1	+ accu(0.5 * ((beta % beta) * priorprecbeta)) + accu(0.5 * ((theta % theta) * priorprectheta)));
409
410		// compute LOG posterior
411	1	LogPos logPos(alpha, psi, beta, theta, wfm, K);
412	1	parallelReduce(0, N, logPos);
413	1	lp += logPos.lp;
414
415
416		// BEGIN WHILE LOOP
417	1	while (((lp - lastlp) > tolvec(0)) && outeriter < OUTERITER) {
418	1	outeriter++;
419		//Rcout<<"alphs_b="<<mean(alpha)<<" theta_B="<<mean(theta)<<std::endl;
420
421		// UPDATE WORD PARAMETERS
422	1	NumericVector psi_N(psi.begin(), psi.end());
423	1	NumericVector beta_N(beta.begin(), beta.end());
424	1	WordPar wordPar(alpha, psi_N, beta_N, theta, phi, wfm, tolvec,
425	1	outeriter, priorprecpsi, priorprecbeta, stepsize);
426	1	parallelFor(0, K, wordPar);
427
428	1	psi = as<arma::rowvec>(psi_N);
429	1	beta = as<arma::rowvec>(beta_N);
430
431		// UPDATE DOCUMENT PARAMETERS
432	1	NumericVector alpha_N(alpha.begin(), alpha.end());
433	1	NumericVector theta_N(theta.begin(), theta.end());
434	1	DocPar docPar(alpha_N, psi, beta, theta_N, phi, wfm, tolvec,
435	1	outeriter, priorprecalpha, priorprectheta, stepsize);
436	1	parallelFor(0, N, docPar);
437
438	1	alpha = as<arma::colvec>(alpha_N);
439	1	theta = as<arma::colvec>(theta_N);
440
441
442		// UPDATE DISPERSION PARAMETERS
443
444	1	if (disptype(0) == 2) { // single dispersion parameter for all words
445	0	DispPar dispPar(alpha, psi, beta, theta, wfm, N);
446	0	parallelReduce(0, K, dispPar);
447	0	double phitmp = dispPar.phitmp;
448	0	phitmp /= N * K - 2 * N - 2 * K;
449	0	phi.fill(phitmp);
450		}
451
452	1	if (disptype(0) >= 3) { // individual dispersion parameter for each word
453	1	NumericVector phi_N(phi.begin(), phi.end());
454	1	DispPar2 dispPar2(alpha, psi, beta, theta, wfm, disptype, dispmin, N, K, phi_N);
455	1	parallelFor(0, K, dispPar2);
456	1	phi = as<arma::rowvec>(phi_N);
457		}
458
459	1	alpha = alpha - mean(alpha);
460	1	theta = (theta - mean(theta))/stddev(theta);
461
462		// CHECK LOG-POSTERIOR FOR CONVERGENCE
463	1	lastlp = lp;
464	1	lp = -1.0 * (accu(0.5 * ((alpha % alpha) * priorprecalpha)) + accu(0.5 * ((psi % psi) * priorprecpsi)) +
465	1	accu(0.5 * ((beta % beta) * priorprecbeta)) + accu(0.5 * ((theta % theta) * priorprectheta)));
466
467	1	LogPos logPos2(alpha, psi, beta, theta, wfm, K);
468	1	parallelReduce(0, N, logPos2);
469	1	lp += logPos2.lp;
470		// END WHILE LOOP
471		}
472
473		// Fix Global Polarity
474		//Rcout<<"outeriter="<<outeriter<<" lp - lastlp= "<<lp - lastlp<<std::endl;
475		// added the -1 because C counts from ZERO... -- KB
476	1	if (theta(dirvec(0) - 1) > theta(dirvec(1) - 1)) {
477	1	beta = -beta;
478	1	theta = -theta;
479		}
480
481		// COMPUTE DOCUMENT STANDARD ERRORS
482	1	NumericVector thetaSE_N(thetaSE.begin(), thetaSE.end());
483	1	DocErr docErr(alpha, psi, beta, theta, phi, priorprecalpha, priorprectheta, thetaSE_N);
484	1	parallelFor(0, N, docErr);
485	1	thetaSE = as<arma::colvec>(thetaSE_N);
486
487		// DEFINE OUTPUT
488	1	return Rcpp::List::create(Rcpp::Named("theta") = theta,
489	1	Rcpp::Named("alpha") = alpha,
490	1	Rcpp::Named("psi") = psi,
491	1	Rcpp::Named("beta") = beta,
492	1	Rcpp::Named("phi") = phi,
493	1	Rcpp::Named("thetaSE") = thetaSE);
494
495		}