<commit_before>#include "ldamcmc.h"


/**
 *  Computes Log Posterior Predictive Value: 
 *  
 *  This is developed by <NAME> (2015) for his dissertation  
 *
 * 	Arguments:
 * 		num_topics_              - number of topics
 * 		vocab_size_              - vocabulary size
 * 		word_ids_                - vocabulary ids of each word in each corpus document
 * 		doc_lengths_             - number of words in each document as a vector
 * 		topic_assignments_       - initial topic assignment of each word in each corpus document
 * 		alpha_v_                 - hyperparameter vector for the document topic Dirichlet
 * 		eta_                     - hyperparameter for the topic Dirichlet
 * 		max_iter_                - max number of Gibbs iterations to be perfomed
 * 		burn_in_                 - burn in period
 * 		spacing_                 - spacing between samples that are stored
 *
 * 	Returns:
 * 		Log posterior predictive value of the corpus 
 *
 */
RcppExport SEXP lda_fgs_lppv(SEXP num_topics_, SEXP vocab_size_, SEXP word_ids_, 
                             SEXP doc_lengths_, SEXP topic_assignments_, 
                             SEXP alpha_v_, SEXP eta_, SEXP max_iter_, 
                             SEXP burn_in_, SEXP spacing_) {
  
  cout << "lda_fgs (c++): init process..." << endl;
  
  // Variable from the R interface  
  
  uvec doc_lengths = as<uvec>(doc_lengths_);
  uvec word_ids = as<uvec>(word_ids_);
  uvec z = as<uvec>(topic_assignments_);
  vec alpha_v = as<vec>(alpha_v_);

  double eta = as<double>(eta_);
  unsigned int num_topics = as<unsigned int>(num_topics_);
  unsigned int vocab_size = as<unsigned int>(vocab_size_);
  unsigned int max_iter = as<unsigned int>(max_iter_);
  unsigned int burn_in = as<unsigned int>(burn_in_);
  unsigned int spacing = as<unsigned int>(spacing_);

  // Function variables 
  unsigned int saved_samples = ceil((max_iter - burn_in) / (double) spacing);
  unsigned int num_docs = doc_lengths.n_elem;

  mat beta_samples = zeros<mat>(num_topics, vocab_size);
  mat beta_counts = zeros<mat>(num_topics, vocab_size);
  mat lppv = zeros<mat>(num_docs, saved_samples);
  double cp, llw_hod, shift_hod, avg_ppv, lppc = 0.0; 
  double smoother = 1e-5;
  
  vector < vector < unsigned int > > doc_word_indices;
  unsigned int d, i, k, iter, ss_idx, instances = 0, hod; 
  unsigned int msg_interval = (100 >= max_iter)?1:100;
  

  // Calculates the indices for each word in a document as a vector
  for (d = 0; d < num_docs; d++){
    vector < unsigned int > word_idx;
    for (i = 0; i < doc_lengths(d); i++){
      word_idx.push_back(instances);
      instances++;
    }
    doc_word_indices.push_back(word_idx);
  }
  
  cout << "lda_fgs (c++): total number of words - " << instances << endl;
  cout << "lda_fgs (c++): completed initialization." << endl;  
  
  
  
  for (hod = 0; hod < num_docs; hod++){ // for each held-out document d 
  
    // double avg_llw_hod = 0.0;  
  
    cout << "lda_fgs (c++): Gibbs sampling [doc #" << hod << "]... " << endl;
    
    ss_idx = 0; 
    
    // Initilizes beta
    
    beta_counts.fill(eta); // initializes with the smoothing parameter
    for (d = 0; d < num_docs; d++){ // for each document
      if (d == hod) continue; // ignores the held-out document    
      vector < unsigned int > word_idx = doc_word_indices[d];
      for (i = 0; i < doc_lengths(d); i++)
        beta_counts(z(i), word_ids(word_idx[i])) += 1;
    }

    // The Gibbs sampling loop
    
    for (iter = 0; iter < max_iter; iter++){ 
      
      if (iter % msg_interval == 0) { 
        cout << "lda_fgs (c++): gibbs iter# " << iter + 1;
      }
      
      // samples \beta
      for(k = 0; k < num_topics; k++)
        beta_samples.row(k) = sample_dirichlet_row_vec(vocab_size, beta_counts.row(k));
      
      
      for (d = 0; d < num_docs; d++){ // for each document
        
        if (d == hod) continue; // ignores the held-out document    
        
        vector < unsigned int > word_idx = doc_word_indices[d];
        
        // samples \theta
        vec partition_counts = alpha_v; // initializes with the smoothing parameter
        for (i = 0; i < doc_lengths(d); i++)
          partition_counts(z(word_idx[i])) += 1;
        vec theta_d = sample_dirichlet(num_topics, partition_counts);

        // samples z and updates \beta counts
        for(i = 0; i < doc_lengths(d); i++)
          beta_counts(z(word_idx[i]), word_ids(word_idx[i])) -= 1; // excludes document d's word-topic counts
        for (i = 0; i < doc_lengths(d); i++)
          z(word_idx[i]) = sample_multinomial(theta_d % beta_samples.col(word_ids(word_idx[i])));
        for(i = 0; i < doc_lengths(d); i++)
          beta_counts(z(word_idx[i]), word_ids(word_idx[i])) += 1; // includes document d's word-topic counts
        
      }
      
      if ((iter >= burn_in) && (iter % spacing == 0)){ // Handles burn in period
        
        // samples theta for held-out document hod  
        vec theta_hod = sample_dirichlet(num_topics, alpha_v); 
        
        // computes the log posterior predictive value of the held-out document 
        // using the current sample/iteration  
        vector < unsigned int > word_idx = doc_word_indices[hod];
        llw_hod = 0.0;  
        for (i = 0; i < doc_lengths(hod); i++){
          cp = accu(theta_hod % beta_samples.col(word_ids(word_idx[i])));
          // cp = accu(theta_hod.t() * beta_samples.col(word_ids(word_idx[i])));
          llw_hod += log(is_finite(cp)?(cp+smoother):smoother);
        }
        lppv(hod, ss_idx) = llw_hod;
        if (iter % msg_interval == 0) { cout << " log(PPV) = " << llw_hod; }

        // Note online average is disabled because of some computational issues 
        // avg_llw_hod = ((ss_idx * avg_llw_hod + exp(llw_hod)) / (ss_idx + 1.)); 

        ss_idx++; // updates the counter  
        
      } // Handles burn in period
      
      
      if (iter % msg_interval == 0) { cout << endl; }
      
    } // The end of the Gibbs loop
    
    cout << "lda_fgs (c++): completed Gibbs sampling." << endl;
    cout << "lda_fgs (c++): number of saved samples - " << ss_idx << endl;
    
    // lppv += log(avg_llw_hod) / (double)num_docs; 
    
    shift_hod = max(lppv.row(hod)) - 20; 
    avg_ppv = log(mean(exp(lppv.row(hod) - shift_hod))) + shift_hod;
    cout << "lda_fgs (c++): log (PPV[avg](" << hod+1 << ")) = " << avg_ppv << endl;
    
    lppc += avg_ppv;
  
  }
  
  
  return List::create(
    Named("lppv") = wrap(lppv), 
    Named("lppc") = wrap(lppc/(double)num_docs)
  );<commit_msg>Code clean-up and documentation<commit_after>#include "ldamcmc.h"


/**
 *  Computes Log Posterior Predictive Value: 
 *  
 *  This is developed by <NAME> (2015) for his dissertation  
 *
 * 	Arguments:
 * 		num_topics_              - number of topics
 * 		vocab_size_              - vocabulary size
 * 		word_ids_                - vocabulary ids of each word in each corpus document
 * 		doc_lengths_             - number of words in each document as a vector
 * 		topic_assignments_       - initial topic assignment of each word in each corpus document
 * 		alpha_v_                 - hyperparameter vector for the document topic Dirichlet
 * 		eta_                     - hyperparameter for the topic Dirichlet
 * 		max_iter_                - max number of Gibbs iterations to be perfomed
 * 		burn_in_                 - burn in period
 * 		spacing_                 - spacing between samples that are stored
 *
 * 	Returns:
 * 		Log posterior predictive value of the corpus 
 *
 */
RcppExport SEXP lda_fgs_lppv(SEXP num_topics_, SEXP vocab_size_, SEXP word_ids_, 
                             SEXP doc_lengths_, SEXP topic_assignments_, 
                             SEXP alpha_v_, SEXP eta_, SEXP max_iter_, 
                             SEXP burn_in_, SEXP spacing_) {
  
  cout << "lda_fgs (c++): init process..." << endl;
  
  // Variable from the R interface  
  
  uvec doc_lengths = as<uvec>(doc_lengths_);
  uvec word_ids = as<uvec>(word_ids_);
  uvec z = as<uvec>(topic_assignments_);
  vec alpha_v = as<vec>(alpha_v_);

  double eta = as<double>(eta_);
  unsigned int num_topics = as<unsigned int>(num_topics_);
  unsigned int vocab_size = as<unsigned int>(vocab_size_);
  unsigned int max_iter = as<unsigned int>(max_iter_);
  unsigned int burn_in = as<unsigned int>(burn_in_);
  unsigned int spacing = as<unsigned int>(spacing_);

  // Function variables 
  unsigned int saved_samples = ceil((max_iter - burn_in) / (double) spacing);
  unsigned int num_docs = doc_lengths.n_elem;

  mat beta_samples = zeros<mat>(num_topics, vocab_size);
  mat beta_counts = zeros<mat>(num_topics, vocab_size);
  mat lppv = zeros<mat>(num_docs, saved_samples);
  double cp, llw_hod, shift_hod, avg_ppv, lppv_sum = 0.0; 
  double smoother = 1e-5;
  
  vector < vector < unsigned int > > doc_word_indices;
  unsigned int d, i, k, iter, ss_idx, instances = 0, hod; 
  unsigned int msg_interval = (100 >= max_iter)?1:100;
  

  // Calculates the indices for each word in a document as a vector
  for (d = 0; d < num_docs; d++){
    vector < unsigned int > word_idx;
    for (i = 0; i < doc_lengths(d); i++){
      word_idx.push_back(instances);
      instances++;
    }
    doc_word_indices.push_back(word_idx);
  }
  
  cout << "lda_fgs (c++): total number of words - " << instances << endl;
  cout << "lda_fgs (c++): completed initialization." << endl;  
  
  
  
  for (hod = 0; hod < num_docs; hod++){ // for each held-out document d 

    cout << "lda_fgs (c++): Gibbs sampling [doc #" << hod << "]... " << endl;
    
    ss_idx = 0; 
    
    // Initilizes beta
    
    beta_counts.fill(eta); // initializes with the smoothing parameter
    for (d = 0; d < num_docs; d++){ // for each document
      if (d == hod) continue; // ignores the held-out document    
      vector < unsigned int > word_idx = doc_word_indices[d];
      for (i = 0; i < doc_lengths(d); i++)
        beta_counts(z(i), word_ids(word_idx[i])) += 1;
    }

    // The Gibbs sampling loop
    
    for (iter = 0; iter < max_iter; iter++){ 
      
      if (iter % msg_interval == 0) { 
        cout << "lda_fgs (c++): gibbs iter# " << iter + 1;
      }
      
      // samples \beta
      
      for(k = 0; k < num_topics; k++)
        beta_samples.row(k) = sample_dirichlet_row_vec(vocab_size, beta_counts.row(k));
      
      
      for (d = 0; d < num_docs; d++){ // for each document
        
        if (d == hod) continue; // ignores the held-out document    
        
        vector < unsigned int > word_idx = doc_word_indices[d];
        
        // samples \theta
        
        vec partition_counts = alpha_v; // initializes with the smoothing parameter
        for (i = 0; i < doc_lengths(d); i++)
          partition_counts(z(word_idx[i])) += 1;
        vec theta_d = sample_dirichlet(num_topics, partition_counts);

        // samples z and updates \beta counts
        
        for(i = 0; i < doc_lengths(d); i++)
          beta_counts(z(word_idx[i]), word_ids(word_idx[i])) -= 1; // excludes document d's word-topic counts
        for (i = 0; i < doc_lengths(d); i++)
          z(word_idx[i]) = sample_multinomial(theta_d % beta_samples.col(word_ids(word_idx[i])));
        for(i = 0; i < doc_lengths(d); i++)
          beta_counts(z(word_idx[i]), word_ids(word_idx[i])) += 1; // includes document d's word-topic counts
        
      }
      
      if ((iter >= burn_in) && (iter % spacing == 0)){ // Handles burn in period
        
        // samples theta for held-out document hod  
        
        vec theta_hod = sample_dirichlet(num_topics, alpha_v); 
        
        // computes the log posterior predictive value of the held-out document 
        // using the current sample/iteration  
        
        vector < unsigned int > word_idx = doc_word_indices[hod];
        llw_hod = 0.0;  
        for (i = 0; i < doc_lengths(hod); i++){
          cp = accu(theta_hod % beta_samples.col(word_ids(word_idx[i])));
          llw_hod += log(is_finite(cp)?(cp+smoother):smoother);
        }
        lppv(hod, ss_idx) = llw_hod;
        if (iter % msg_interval == 0) { cout << " log(PPV) = " << llw_hod; }

        ss_idx++; // updates the counter  
        
      } // Handles burn in period
      
      
      if (iter % msg_interval == 0) { cout << endl; }
      
    } // The end of the Gibbs loop
    
    cout << "lda_fgs (c++): completed Gibbs sampling." << endl;
    cout << "lda_fgs (c++): number of saved samples - " << ss_idx << endl;
    
    // Computes the estimate of the predictive likelihood of the held-out 
    // document hod  
    
    shift_hod = max(lppv.row(hod)) - 20; 
    avg_ppv = log(mean(exp(lppv.row(hod) - shift_hod))) + shift_hod;
    
    cout << "lda_fgs (c++): log (PPV[avg](" << hod+1 << ")) = " << avg_ppv << endl;
    
    lppv_sum += avg_ppv;
  
  }
  
  
  return List::create(
    Named("lppv") = wrap(lppv), 
    Named("lppc") = wrap(lppv_sum/(double)num_docs)
  );
{'fraction_non_alphanumeric': 0.12503856834310398, 'fraction_numerical': 0.007019438444924406, 'mean_word_length': 2.2822784810126584, 'pattern_counts': {'":': 0, '<': 173, '<?xml version=': 0, '>': 47, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6910875', 'n_tokens_mistral': 4647, 'n_tokens_neox': 4376, 'n_words': 1176}