Skip to content

Commit

Permalink
Removed likelihood computation from pll_compute_likelihood_derivatives
Browse files Browse the repository at this point in the history 
Although per-site likelihoods need to be computed in pll_compute_likelihood_derivatives,
applying the logarithm to each site can produce a significant overhead during N-R
optimization process, while it is not needed.
  • Loading branch information
ddarriba committed Jul 23, 2016
1 parent 45e7aea commit c0f1cb0
Show file tree
Hide file tree
Showing 5 changed files with 183 additions and 94 deletions.
53 changes: 20 additions & 33 deletions src/core_likelihood.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -832,28 +832,27 @@ static void core_site_likelihood_derivatives(unsigned int states,
}
}

PLL_EXPORT double pll_core_likelihood_derivatives(unsigned int states,
unsigned int sites,
unsigned int rate_cats,
const double * rate_weights,
const unsigned int * parent_scaler,
const unsigned int * child_scaler,
const int * invariant,
const unsigned int * pattern_weights,
double branch_length,
const double * prop_invar,
double ** freqs,
const const double * rates,
double ** eigenvals,
const double * sumtable,
double * d_f,
double * dd_f,
unsigned int attrib)
PLL_EXPORT int pll_core_likelihood_derivatives(unsigned int states,
unsigned int sites,
unsigned int rate_cats,
const double * rate_weights,
const unsigned int * parent_scaler,
const unsigned int * child_scaler,
const int * invariant,
const unsigned int * pattern_weights,
double branch_length,
const double * prop_invar,
double ** freqs,
const const double * rates,
double ** eigenvals,
const double * sumtable,
double * d_f,
double * dd_f,
unsigned int attrib)
{
unsigned int n, i, j;
unsigned int ef_sites;
double site_lk[3];
double logLK = 0.0;

const double * sum;
double deriv1, deriv2;
Expand Down Expand Up @@ -897,7 +896,7 @@ PLL_EXPORT double pll_core_likelihood_derivatives(unsigned int states,
{
pll_errno = PLL_ERROR_MEM_ALLOC;
snprintf (pll_errmsg, 200, "Cannot allocate memory for diagptable");
return -INFINITY;
return PLL_FAILURE;
}

/* pre-compute the derivatives of the P matrix for all discrete GAMMA rates */
Expand Down Expand Up @@ -936,12 +935,6 @@ PLL_EXPORT double pll_core_likelihood_derivatives(unsigned int states,
scale_factors = (parent_scaler) ? parent_scaler[n] : 0;
scale_factors += (child_scaler) ? child_scaler[n] : 0;

logLK += log (site_lk[0]) * pattern_weights[n];
if (scale_factors)
{
logLK += scale_factors * log (PLL_SCALE_THRESHOLD);
}

/* build derivatives */
deriv1 = (-site_lk[1] / site_lk[0]);
deriv2 = (deriv1 * deriv1 - (site_lk[2] / site_lk[0]));
Expand Down Expand Up @@ -994,19 +987,13 @@ PLL_EXPORT double pll_core_likelihood_derivatives(unsigned int states,
switch(asc_bias_type)
{
case PLL_ATTRIB_ASC_BIAS_LEWIS:
/* correct log-likelihood */
logLK -= pattern_weight_sum * log(1.0 - asc_Lk[0]);

/* derivatives of log(1.0 - (sum Li(s) over states 's')) */
*d_f -= pattern_weight_sum * (asc_Lk[1] / (asc_Lk[0] - 1.0));
*dd_f -= pattern_weight_sum *
(((asc_Lk[0] - 1.0) * asc_Lk[2] - asc_Lk[1] * asc_Lk[1]) /
((asc_Lk[0] - 1.0) * (asc_Lk[0] - 1.0)));
break;
case PLL_ATTRIB_ASC_BIAS_FELSENSTEIN:
/* correct log-likelihood */
logLK += sum_w_inv * log(asc_Lk[0]);

/* derivatives of log(sum Li(s) over states 's') */
*d_f += sum_w_inv * (asc_Lk[1] / asc_Lk[0]);
*dd_f += sum_w_inv *
Expand All @@ -1016,11 +1003,11 @@ PLL_EXPORT double pll_core_likelihood_derivatives(unsigned int states,
default:
pll_errno = PLL_ERROR_ASC_BIAS;
snprintf(pll_errmsg, 200, "Illegal ascertainment bias algorithm");
return -INFINITY;
return PLL_FAILURE;
}
}
}

free (diagptable);
return logLK;
return PLL_SUCCESS;
}
52 changes: 26 additions & 26 deletions src/likelihood.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -930,14 +930,14 @@ PLL_EXPORT int pll_update_sumtable(pll_partition_t * partition,
* d_f: [output] first derivative
* dd_f: [output] second derivative
*/
PLL_EXPORT double pll_compute_likelihood_derivatives(pll_partition_t * partition,
int parent_scaler_index,
int child_scaler_index,
double branch_length,
const unsigned int * params_indices,
const double * sumtable,
double * d_f,
double * dd_f)
PLL_EXPORT int pll_compute_likelihood_derivatives(pll_partition_t * partition,
int parent_scaler_index,
int child_scaler_index,
double branch_length,
const unsigned int * params_indices,
const double * sumtable,
double * d_f,
double * dd_f)
{
unsigned int * parent_scaler;
unsigned int * child_scaler;
Expand Down Expand Up @@ -966,27 +966,27 @@ PLL_EXPORT double pll_compute_likelihood_derivatives(pll_partition_t * partition
else
child_scaler = partition->scale_buffer[child_scaler_index];

double logLK = pll_core_likelihood_derivatives(partition->states,
partition->sites,
partition->rate_cats,
partition->rate_weights,
parent_scaler,
child_scaler,
partition->invariant,
partition->pattern_weights,
branch_length,
prop_invar,
freqs,
partition->rates,
eigenvals,
sumtable,
d_f,
dd_f,
partition->attributes);
int retval = pll_core_likelihood_derivatives(partition->states,
partition->sites,
partition->rate_cats,
partition->rate_weights,
parent_scaler,
child_scaler,
partition->invariant,
partition->pattern_weights,
branch_length,
prop_invar,
freqs,
partition->rates,
eigenvals,
sumtable,
d_f,
dd_f,
partition->attributes);

free (freqs);
free (prop_invar);
free (eigenvals);

return logLK;
return retval;
}
36 changes: 18 additions & 18 deletions src/pll.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -440,7 +440,7 @@ PLL_EXPORT int pll_update_sumtable(pll_partition_t * partition,
const unsigned int * params_indices,
double *sumtable);

PLL_EXPORT double pll_compute_likelihood_derivatives(pll_partition_t * partition,
PLL_EXPORT int pll_compute_likelihood_derivatives(pll_partition_t * partition,
int parent_scaler_index,
int child_scaler_index,
double branch_length,
Expand Down Expand Up @@ -671,23 +671,23 @@ PLL_EXPORT int pll_core_update_sumtable_ti(unsigned int states,
double *sumtable,
unsigned int attrib);

PLL_EXPORT double pll_core_likelihood_derivatives(unsigned int states,
unsigned int sites,
unsigned int rate_cats,
const double * rate_weights,
const unsigned int * parent_scaler,
const unsigned int * child_scaler,
const int * invariant,
const unsigned int * pattern_weights,
double branch_length,
const double * prop_invar,
double ** freqs,
const double * rates,
double ** eigenvals,
const double * sumtable,
double * d_f,
double * dd_f,
unsigned int attrib);
PLL_EXPORT int pll_core_likelihood_derivatives(unsigned int states,
unsigned int sites,
unsigned int rate_cats,
const double * rate_weights,
const unsigned int * parent_scaler,
const unsigned int * child_scaler,
const int * invariant,
const unsigned int * pattern_weights,
double branch_length,
const double * prop_invar,
double ** freqs,
const const double * rates,
double ** eigenvals,
const double * sumtable,
double * d_f,
double * dd_f,
unsigned int attrib);

/* functions in core_likelihood_avx.c */

Expand Down
35 changes: 27 additions & 8 deletions test/src/asc-bias.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,14 +118,33 @@ static double eval(pll_partition_t * partition,
for (i=0; i<NUM_BRANCH_LENGTHS; ++i)
{
double branch_length = test_branch_lengths[i];
logl = pll_compute_likelihood_derivatives(partition,
tree->scaler_index,
tree->back->scaler_index,
branch_length,
params_indices,
sumtable,
&d_f,
&dd_f);
if (!pll_compute_likelihood_derivatives(partition,
tree->scaler_index,
tree->back->scaler_index,
branch_length,
params_indices,
sumtable,
&d_f,
&dd_f))
{
printf("Error computing likelihood derivatives\n");
exit(1);
}

/* update logLikelihood */
pll_update_prob_matrices(partition,
params_indices,
&(tree->pmatrix_index),
&branch_length,
1);
logl = pll_compute_edge_loglikelihood(partition,
tree->clv_index,
tree->scaler_index,
tree->back->clv_index,
tree->back->scaler_index,
tree->pmatrix_index,
params_indices,
NULL);

printf("%8.4f %18.6f %15.8e %15.8e ", branch_length, logl, d_f, dd_f);
if (logl > max_logl)
Expand Down
Loading

0 comments on commit c0f1cb0

Please sign in to comment.