<commit_before>#ifndef INSIDE_MONTE_CARLO_FMM_H
#error "This is not a public header file!"
#endif

template<typename TKernelAux>
void MonteCarloFMM<TKernelAux>::BaseCase_
(const Matrix &query_set, const ArrayList<index_t> &query_index_permutation,
 QueryTree *query_node, const ReferenceTree *reference_node,
 Vector &query_kernel_sums, double one_sided_probability,
 Vector &query_extra_probability) const {
 
  // Loop over each query point in the query node.
  for(index_t q = query_node->begin(); q < query_node->end(); q++) {
    
    // Get the pointer to the current query point.
    const double *query_point = query_set.GetColumnPtr(q);
    
    // Loop over each reference point in the reference node.
    for(index_t r = reference_node->begin(); r < reference_node->end(); 
	r++) {
      
      // Get the pointer to the current reference point.
      const double *reference_point = reference_set_.GetColumnPtr(r);
      
      // Compute the pairwise distance and kernel value.
      double squared_distance = 
	la::DistanceSqEuclidean(query_set.n_rows(), query_point,
				reference_point);

      double weighted_kernel_value = reference_weights_[r] *
	ka_.kernel_.EvalUnnormOnSq(squared_distance);
      
      query_kernel_sums[query_index_permutation[q]] += weighted_kernel_value;
      
    } // end of iterating over each reference point.

    query_extra_probability[q] += (1 - one_sided_probability);

  } // end of looping over each query point.
}

template<typename TKernelAux>
bool MonteCarloFMM<TKernelAux>::Prunable_
(const QueryTree *query_node, 
 const ArrayList<index_t> &query_index_permutation,
 const ReferenceTree *reference_node, Vector &query_kernel_sums) {
  
  // Try pruning after bound refinement: first compute distance/kernel
  // value bounds.
  DRange dsqd_range, kernel_value_range;
  dsqd_range.lo = query_node->bound().MinDistanceSq(reference_node->bound());
  dsqd_range.hi = query_node->bound().MaxDistanceSq(reference_node->bound());
  kernel_value_range = ka_.kernel_.RangeUnnormOnSq(dsqd_range);

  double used_error = 0.5 * (kernel_value_range.hi - kernel_value_range.lo) *
    reference_node->count();

  // This can probably done more in a smart way...
  double new_mass_l = DBL_MAX;
  for(index_t q = query_node->begin(); q < query_node->end(); q++) {
    new_mass_l = std::min(new_mass_l, 
			  query_kernel_sums[query_index_permutation[q]]);
  }
  new_mass_l += kernel_value_range.lo * reference_node->count();

  double allowed_error = (relative_error_ * new_mass_l) *
    reference_node->count() / ((double) reference_tree_root_->count());

  if(used_error <= allowed_error) {
    for(index_t q = query_node->begin(); q < query_node->end(); q++) {
      query_kernel_sums[query_index_permutation[q]] +=
	kernel_value_range.lo * reference_node->count();
    } 
    return true;
  }
  else {
    return false;
  }
}

template<typename TKernelAux>
bool MonteCarloFMM<TKernelAux>::MonteCarloPrunable_
(const Matrix &query_set, const ArrayList<index_t> &query_index_permutation,
 QueryTree *query_node, const ReferenceTree *reference_node,
 Vector &query_kernel_sums, Vector &query_kernel_sums_scratch_space,
 Vector &query_squared_kernel_sums_scratch_space,
 double one_sided_probability, Vector &query_extra_probability) {

  if(num_initial_samples_per_query_ > reference_node->count()) {
    return false;
  }

  // For each query point in the query node, take samples and
  // determine how many more samples are needed.
  bool flag = true;
  for(index_t q = query_node->begin(); q < query_node->end() && flag; q++) {

    // Compute the required standard score for the given one-sided
    // probability.
    double standard_score = 
      InverseNormalCDF::Compute
      (std::max(one_sided_probability - query_extra_probability[q], 0.0));
    
    // Get the pointer to the current query point.
    const double *query_point = query_set.GetColumnPtr(q);

    // Reset the current position of the scratch space to zero.
    query_kernel_sums_scratch_space[q] = 0;
    query_squared_kernel_sums_scratch_space[q] = 0;
    
    // The initial number of samples is equal to the default.
    int num_samples = num_initial_samples_per_query_;
    int total_samples = 0;

    do {
      for(index_t s = 0; s < num_samples; s++) {
	index_t random_reference_point_index = 
	  math::RandInt(reference_node->begin(), reference_node->end());
	
	// Get the pointer to the current reference point.
	const double *reference_point = 
	  reference_set_.GetColumnPtr(random_reference_point_index);
	
	// Compute the pairwise distance and kernel value.
	double squared_distance = 
	  la::DistanceSqEuclidean(query_set.n_rows(), query_point,
				  reference_point);

	double weighted_kernel_value = 
	  ka_.kernel_.EvalUnnormOnSq(squared_distance);
	query_kernel_sums_scratch_space[q] += weighted_kernel_value;
	query_squared_kernel_sums_scratch_space[q] +=
	  weighted_kernel_value * weighted_kernel_value;
      }
      total_samples += num_samples;

      // Compute the current estimate of the sample mean and the
      // sample variance.
      double sample_mean = query_kernel_sums_scratch_space[q] / 
	((double) total_samples);
      double sample_variance =
	(query_squared_kernel_sums_scratch_space[q] -
	 total_samples * sample_mean * sample_mean) /
	((double) total_samples - 1);

      // Compute the current threshold for guaranteeing the relative
      // error bound.
      int threshold = (sample_variance > DBL_EPSILON) ?
	(int) ceil(math::Sqr(standard_score * reference_tree_root_->count() /
			     relative_error_ / 
			     (sample_mean * reference_node->count() + 
			      query_kernel_sums[query_index_permutation[q]]))
		   * sample_variance):
	0;
      num_samples = threshold - total_samples;

      // If it will require too many samples, give up.
      if(num_samples > 10 * reference_node->count()) {
	flag = false;
	break;
      }
      
      // If we are done, then move onto the next query.
      else if(num_samples <= 0) {
	query_kernel_sums_scratch_space[q] /= ((double) total_samples);
	break;
      }

    } while(1);
  } // end of looping over each query...

  // If all queries can be pruned, then add the approximations.
  if(flag) {
    for(index_t q = query_node->begin(); q < query_node->end(); q++) {
      query_kernel_sums[query_index_permutation[q]] += 
	query_kernel_sums_scratch_space[q] * reference_node->count();

      // This can potentially be improved...
      query_extra_probability[q] = 0;
    }    
    return true;
  }
  
  return false;<commit_msg>Additional parameter tweaking for performance enhancement.<commit_after>#ifndef INSIDE_MONTE_CARLO_FMM_H
#error "This is not a public header file!"
#endif

template<typename TKernelAux>
void MonteCarloFMM<TKernelAux>::BaseCase_
(const Matrix &query_set, const ArrayList<index_t> &query_index_permutation,
 QueryTree *query_node, const ReferenceTree *reference_node,
 Vector &query_kernel_sums, double one_sided_probability,
 Vector &query_extra_probability) const {
 
  // Loop over each query point in the query node.
  for(index_t q = query_node->begin(); q < query_node->end(); q++) {
    
    // Get the pointer to the current query point.
    const double *query_point = query_set.GetColumnPtr(q);
    
    // Loop over each reference point in the reference node.
    for(index_t r = reference_node->begin(); r < reference_node->end(); 
	r++) {
      
      // Get the pointer to the current reference point.
      const double *reference_point = reference_set_.GetColumnPtr(r);
      
      // Compute the pairwise distance and kernel value.
      double squared_distance = 
	la::DistanceSqEuclidean(query_set.n_rows(), query_point,
				reference_point);

      double weighted_kernel_value = reference_weights_[r] *
	ka_.kernel_.EvalUnnormOnSq(squared_distance);
      
      query_kernel_sums[query_index_permutation[q]] += weighted_kernel_value;
      
    } // end of iterating over each reference point.

    query_extra_probability[q] += (1 - one_sided_probability);

  } // end of looping over each query point.
}

template<typename TKernelAux>
bool MonteCarloFMM<TKernelAux>::Prunable_
(const QueryTree *query_node, 
 const ArrayList<index_t> &query_index_permutation,
 const ReferenceTree *reference_node, Vector &query_kernel_sums) {
  
  // Try pruning after bound refinement: first compute distance/kernel
  // value bounds.
  DRange dsqd_range, kernel_value_range;
  dsqd_range.lo = query_node->bound().MinDistanceSq(reference_node->bound());
  dsqd_range.hi = query_node->bound().MaxDistanceSq(reference_node->bound());
  kernel_value_range = ka_.kernel_.RangeUnnormOnSq(dsqd_range);

  double used_error = 0.5 * (kernel_value_range.hi - kernel_value_range.lo) *
    reference_node->count();

  // This can probably done more in a smart way...
  double new_mass_l = DBL_MAX;
  for(index_t q = query_node->begin(); q < query_node->end(); q++) {
    new_mass_l = std::min(new_mass_l, 
			  query_kernel_sums[query_index_permutation[q]]);
  }
  new_mass_l += kernel_value_range.lo * reference_node->count();

  double allowed_error = (relative_error_ * new_mass_l) *
    reference_node->count() / ((double) reference_tree_root_->count());

  if(used_error <= allowed_error) {
    for(index_t q = query_node->begin(); q < query_node->end(); q++) {
      query_kernel_sums[query_index_permutation[q]] +=
	kernel_value_range.lo * reference_node->count();
    } 
    return true;
  }
  else {
    return false;
  }
}

template<typename TKernelAux>
bool MonteCarloFMM<TKernelAux>::MonteCarloPrunable_
(const Matrix &query_set, const ArrayList<index_t> &query_index_permutation,
 QueryTree *query_node, const ReferenceTree *reference_node,
 Vector &query_kernel_sums, Vector &query_kernel_sums_scratch_space,
 Vector &query_squared_kernel_sums_scratch_space,
 double one_sided_probability, Vector &query_extra_probability) {

  if(num_initial_samples_per_query_ > reference_node->count()) {
    return false;
  }

  // For each query point in the query node, take samples and
  // determine how many more samples are needed.
  bool flag = true;
  for(index_t q = query_node->begin(); q < query_node->end() && flag; q++) {

    // Compute the required standard score for the given one-sided
    // probability.
    double standard_score = 
      InverseNormalCDF::Compute
      (std::max(one_sided_probability - query_extra_probability[q], 0.0));
    
    // Get the pointer to the current query point.
    const double *query_point = query_set.GetColumnPtr(q);

    // Reset the current position of the scratch space to zero.
    query_kernel_sums_scratch_space[q] = 0;
    query_squared_kernel_sums_scratch_space[q] = 0;
    
    // The initial number of samples is equal to the default.
    int num_samples = num_initial_samples_per_query_;
    int total_samples = 0;

    do {
      for(index_t s = 0; s < num_samples; s++) {
	index_t random_reference_point_index = 
	  math::RandInt(reference_node->begin(), reference_node->end());
	
	// Get the pointer to the current reference point.
	const double *reference_point = 
	  reference_set_.GetColumnPtr(random_reference_point_index);
	
	// Compute the pairwise distance and kernel value.
	double squared_distance = 
	  la::DistanceSqEuclidean(query_set.n_rows(), query_point,
				  reference_point);

	double weighted_kernel_value = 
	  ka_.kernel_.EvalUnnormOnSq(squared_distance);
	query_kernel_sums_scratch_space[q] += weighted_kernel_value;
	query_squared_kernel_sums_scratch_space[q] +=
	  weighted_kernel_value * weighted_kernel_value;
      }
      total_samples += num_samples;

      // Compute the current estimate of the sample mean and the
      // sample variance.
      double sample_mean = query_kernel_sums_scratch_space[q] / 
	((double) total_samples);
      double sample_variance =
	(query_squared_kernel_sums_scratch_space[q] -
	 total_samples * sample_mean * sample_mean) /
	((double) total_samples - 1);

      // Compute the current threshold for guaranteeing the relative
      // error bound.
      int threshold = (sample_variance > DBL_EPSILON) ?
	(int) ceil(math::Sqr(standard_score * reference_tree_root_->count() /
			     relative_error_ / 
			     (sample_mean * reference_node->count() + 
			      query_kernel_sums[query_index_permutation[q]]))
		   * sample_variance):
	0;
      num_samples = threshold - total_samples;

      // If it will require too many samples, give up.
      if(num_samples > reference_node->count()) {
	flag = false;
	break;
      }
      
      // If we are done, then move onto the next query.
      else if(num_samples <= 0) {
	query_kernel_sums_scratch_space[q] /= ((double) total_samples);
	break;
      }

    } while(1);
  } // end of looping over each query...

  // If all queries can be pruned, then add the approximations.
  if(flag) {
    for(index_t q = query_node->begin(); q < query_node->end(); q++) {
      query_kernel_sums[query_index_permutation[q]] += 
	query_kernel_sums_scratch_space[q] * reference_node->count();

      // This can potentially be improved...
      query_extra_probability[q] = 0;
    }    
    return true;
  }
  
  return false;
{'fraction_non_alphanumeric': 0.10179549604382228, 'fraction_numerical': 0.002282410225197809, 'mean_word_length': 3.892072943803498, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 83, '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': '1217954', 'n_tokens_mistral': 4390, 'n_tokens_neox': 4147, 'n_words': 985}