<commit_before>#include "optimize.hpp"

#include <vector>
#include <cmath>
#include <cstdlib>
#include <cassert>
#include <stdexcept>
#include <limits>
#include <algorithm>

#include "pll_util.hpp"
#include "constants.hpp"
#include "Log.hpp"

using namespace std;

static void traverse_update_partials(pll_utree_t * tree, pll_partition_t * partition,
    pll_utree_t ** travbuffer, double * branch_lengths, unsigned int * matrix_indices,
    pll_operation_t * operations)
{
  unsigned int num_matrices, num_ops;
  vector<unsigned int> param_indices(partition->rate_cats, 0);
  /* perform a full traversal*/
  assert(tree->next != nullptr);
  unsigned int traversal_size;
  // TODO this only needs to be done once, outside of this func. pass traversal size also
  // however this is practically nonexistent impact compared to clv comp
  pll_utree_traverse(tree, cb_full_traversal, travbuffer, &traversal_size);

  /* given the computed traversal descriptor, generate the operations
     structure, and the corresponding probability matrix indices that
     may need recomputing */
  pll_utree_create_operations(travbuffer,
                              traversal_size,
                              branch_lengths,
                              matrix_indices,
                              operations,
                              &num_matrices,
                              &num_ops);

  pll_update_prob_matrices(partition,
                           &param_indices[0],
                           matrix_indices,// matrices to update
                           branch_lengths,
                           num_matrices); // how many should be updated

  /* use the operations array to compute all num_ops inner CLVs. Operations
     will be carried out sequentially starting from operation 0 towrds num_ops-1 */
  pll_update_partials(partition, operations, num_ops);

}

static void utree_derivative_func (void * parameters, double proposal,
                                     double *df, double *ddf)
{
  pll_newton_tree_params_t * params = (pll_newton_tree_params_t *) parameters;
  pll_compute_likelihood_derivatives (
      params->partition,
      params->tree->scaler_index,
      params->tree->back->scaler_index,
      proposal,
      params->params_indices,
      params->sumtable, df, ddf);
}

/**
 * Branch length optimization akin to how pplacer does it: slide the pendant branch along
 * insertion branch of the reference tree, optimize pendant length fully.
 * 
 * @param  partition  the partition
 * @param  tree       the tree structure
 * @param  smoothings maximum number of iterations
 * @return            negative log likelihood after optimization
 */
static double opt_branch_lengths_pplacer(pll_partition_t * partition, pll_utree_t * tree, unsigned int smoothings, const double tolerance)
{
  double loglikelihood = 0.0, new_loglikelihood;
  double xmin,    /* min branch length */
         xguess,  /* initial guess */
         xmax,    /* max branch length */
         xtol,    /* tolerance */
         xres;    /* optimal found branch length */

  const auto score_node = tree;
  const auto blo_node = tree->next->back;
  const auto blo_antinode = tree->next->next->back;

  const auto original_length = blo_node->length * 2;
  bool opt_proximal = (blo_node->length > PLLMOD_OPT_MIN_BRANCH_LEN);

  double lengths[3] = {blo_node->length, blo_antinode->length, score_node->length};
  unsigned int p_indices[3] = {blo_node->pmatrix_index, blo_antinode->pmatrix_index, score_node->pmatrix_index};

  if(!opt_proximal && original_length > PLLMOD_OPT_MIN_BRANCH_LEN)
  {
    blo_node->length = PLLMOD_OPT_MIN_BRANCH_LEN;
    blo_antinode->length = original_length - PLLMOD_OPT_MIN_BRANCH_LEN;
    pll_update_prob_matrices(partition, &param_indices[0], p_indices, lengths, 1);
    pll_update_partials(partition, &toward_score, 1);
    opt_proximal = true;
  }

  vector<unsigned int> param_indices(partition->rate_cats, 0);

  /* set parameters for N-R optimization */
  pll_newton_tree_params_t nr_params;
  nr_params.partition         = partition;
  // nr_params.tree              = score_node;
  nr_params.params_indices    = &param_indices[0];
  // nr_params.branch_length_min = PLLMOD_OPT_MIN_BRANCH_LEN;
  // nr_params.branch_length_max = PLLMOD_OPT_MAX_BRANCH_LEN;
  // nr_params.tolerance         = tolerance;
  nr_params.sumtable          = nullptr;

  pll_operation_t toward_score;
  toward_score.parent_clv_index = score_node->clv_index;
  toward_score.parent_scaler_index = score_node->scaler_index;
  toward_score.child1_clv_index = blo_node->clv_index;
  toward_score.child1_scaler_index = blo_node->scaler_index;
  toward_score.child1_matrix_index = blo_node->pmatrix_index;
  toward_score.child2_clv_index = blo_antinode->clv_index;
  toward_score.child2_scaler_index = blo_antinode->scaler_index;
  toward_score.child2_matrix_index = blo_antinode->pmatrix_index;

  pll_operation_t toward_blo_node;
  toward_blo_node.parent_clv_index = blo_node->back->clv_index;
  toward_blo_node.parent_scaler_index = blo_node->back->scaler_index;
  toward_blo_node.child1_clv_index = score_node->back->clv_index;
  toward_blo_node.child1_scaler_index = score_node->back->scaler_index;
  toward_blo_node.child1_matrix_index = score_node->pmatrix_index;
  toward_blo_node.child2_clv_index = blo_antinode->clv_index;
  toward_blo_node.child2_scaler_index = blo_antinode->scaler_index;
  toward_blo_node.child2_matrix_index = blo_antinode->pmatrix_index;

  /* get the initial likelihood score */

  loglikelihood = -pll_compute_edge_loglikelihood (partition,
                                                  score_node->back->clv_index,
                                                  score_node->back->scaler_index,
                                                  score_node->clv_index,
                                                  score_node->scaler_index,
                                                  score_node->pmatrix_index,
                                                  &param_indices[0],
                                                  nullptr);

  /* allocate the sumtable */
  auto sites_alloc = partition->sites;
  if (partition->attributes & PLL_ATTRIB_AB_FLAG)
    sites_alloc += partition->states;

  if ((nr_params.sumtable = (double *) pll_aligned_alloc(
       sites_alloc * partition->rate_cats * partition->states_padded *
       sizeof(double), partition->alignment)) == nullptr)
  {<commit_msg>fixed compile error<commit_after>#include "optimize.hpp"

#include <vector>
#include <cmath>
#include <cstdlib>
#include <cassert>
#include <stdexcept>
#include <limits>
#include <algorithm>

#include "pll_util.hpp"
#include "constants.hpp"
#include "Log.hpp"

using namespace std;

static void traverse_update_partials(pll_utree_t * tree, pll_partition_t * partition,
    pll_utree_t ** travbuffer, double * branch_lengths, unsigned int * matrix_indices,
    pll_operation_t * operations)
{
  unsigned int num_matrices, num_ops;
  vector<unsigned int> param_indices(partition->rate_cats, 0);
  /* perform a full traversal*/
  assert(tree->next != nullptr);
  unsigned int traversal_size;
  // TODO this only needs to be done once, outside of this func. pass traversal size also
  // however this is practically nonexistent impact compared to clv comp
  pll_utree_traverse(tree, cb_full_traversal, travbuffer, &traversal_size);

  /* given the computed traversal descriptor, generate the operations
     structure, and the corresponding probability matrix indices that
     may need recomputing */
  pll_utree_create_operations(travbuffer,
                              traversal_size,
                              branch_lengths,
                              matrix_indices,
                              operations,
                              &num_matrices,
                              &num_ops);

  pll_update_prob_matrices(partition,
                           &param_indices[0],
                           matrix_indices,// matrices to update
                           branch_lengths,
                           num_matrices); // how many should be updated

  /* use the operations array to compute all num_ops inner CLVs. Operations
     will be carried out sequentially starting from operation 0 towrds num_ops-1 */
  pll_update_partials(partition, operations, num_ops);

}

static void utree_derivative_func (void * parameters, double proposal,
                                     double *df, double *ddf)
{
  pll_newton_tree_params_t * params = (pll_newton_tree_params_t *) parameters;
  pll_compute_likelihood_derivatives (
      params->partition,
      params->tree->scaler_index,
      params->tree->back->scaler_index,
      proposal,
      params->params_indices,
      params->sumtable, df, ddf);
}

/**
 * Branch length optimization akin to how pplacer does it: slide the pendant branch along
 * insertion branch of the reference tree, optimize pendant length fully.
 * 
 * @param  partition  the partition
 * @param  tree       the tree structure
 * @param  smoothings maximum number of iterations
 * @return            negative log likelihood after optimization
 */
static double opt_branch_lengths_pplacer(pll_partition_t * partition, pll_utree_t * tree, unsigned int smoothings, const double tolerance)
{
  double loglikelihood = 0.0, new_loglikelihood;
  double xmin,    /* min branch length */
         xguess,  /* initial guess */
         xmax,    /* max branch length */
         xtol,    /* tolerance */
         xres;    /* optimal found branch length */
  vector<unsigned int> param_indices(partition->rate_cats, 0);

  const auto score_node = tree;
  const auto blo_node = tree->next->back;
  const auto blo_antinode = tree->next->next->back;

  pll_operation_t toward_score;
  toward_score.parent_clv_index = score_node->clv_index;
  toward_score.parent_scaler_index = score_node->scaler_index;
  toward_score.child1_clv_index = blo_node->clv_index;
  toward_score.child1_scaler_index = blo_node->scaler_index;
  toward_score.child1_matrix_index = blo_node->pmatrix_index;
  toward_score.child2_clv_index = blo_antinode->clv_index;
  toward_score.child2_scaler_index = blo_antinode->scaler_index;
  toward_score.child2_matrix_index = blo_antinode->pmatrix_index;

  pll_operation_t toward_blo_node;
  toward_blo_node.parent_clv_index = blo_node->back->clv_index;
  toward_blo_node.parent_scaler_index = blo_node->back->scaler_index;
  toward_blo_node.child1_clv_index = score_node->back->clv_index;
  toward_blo_node.child1_scaler_index = score_node->back->scaler_index;
  toward_blo_node.child1_matrix_index = score_node->pmatrix_index;
  toward_blo_node.child2_clv_index = blo_antinode->clv_index;
  toward_blo_node.child2_scaler_index = blo_antinode->scaler_index;
  toward_blo_node.child2_matrix_index = blo_antinode->pmatrix_index;

  const auto original_length = blo_node->length * 2;
  bool opt_proximal = (blo_node->length > PLLMOD_OPT_MIN_BRANCH_LEN);

  double lengths[3] = {blo_node->length, blo_antinode->length, score_node->length};
  unsigned int p_indices[3] = {blo_node->pmatrix_index, blo_antinode->pmatrix_index, score_node->pmatrix_index};

  if(!opt_proximal && original_length > PLLMOD_OPT_MIN_BRANCH_LEN)
  {
    blo_node->length = PLLMOD_OPT_MIN_BRANCH_LEN;
    blo_antinode->length = original_length - PLLMOD_OPT_MIN_BRANCH_LEN;
    pll_update_prob_matrices(partition, &param_indices[0], p_indices, lengths, 1);
    pll_update_partials(partition, &toward_score, 1);
    opt_proximal = true;
  }

  /* set parameters for N-R optimization */
  pll_newton_tree_params_t nr_params;
  nr_params.partition         = partition;
  // nr_params.tree              = score_node;
  nr_params.params_indices    = &param_indices[0];
  // nr_params.branch_length_min = PLLMOD_OPT_MIN_BRANCH_LEN;
  // nr_params.branch_length_max = PLLMOD_OPT_MAX_BRANCH_LEN;
  // nr_params.tolerance         = tolerance;
  nr_params.sumtable          = nullptr;

  /* get the initial likelihood score */
  loglikelihood = -pll_compute_edge_loglikelihood (partition,
                                                  score_node->back->clv_index,
                                                  score_node->back->scaler_index,
                                                  score_node->clv_index,
                                                  score_node->scaler_index,
                                                  score_node->pmatrix_index,
                                                  &param_indices[0],
                                                  nullptr);

  /* allocate the sumtable */
  auto sites_alloc = partition->sites;
  if (partition->attributes & PLL_ATTRIB_AB_FLAG)
    sites_alloc += partition->states;

  if ((nr_params.sumtable = (double *) pll_aligned_alloc(
       sites_alloc * partition->rate_cats * partition->states_padded *
       sizeof(double), partition->alignment)) == nullptr)
  {
{'fraction_non_alphanumeric': 0.08295929422689986, 'fraction_numerical': 0.004178919671877418, 'mean_word_length': 2.8404160475482914, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 143, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 30, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4179215', 'n_tokens_mistral': 3893, 'n_tokens_neox': 3794, 'n_words': 806}