Skip to content

Commit

Permalink
Handle errors from GSL in gamma site model
Browse files Browse the repository at this point in the history
  • Loading branch information
@4ment
4ment committed Jan 8, 2024
1 parent ab3524e commit 9c94f57
Show file tree
Hide file tree
Showing 5 changed files with 103 additions and 36 deletions.
10 changes: 5 additions & 5 deletions src/phyc/gamma.c
View file
Expand Up @@ -119,7 +119,7 @@ double gser( double *gamser, double a, double x){
double gln = gammln(a);

if (x <= 0.0) {
if (x < 0.0) error("x less than 0 in routine gser");
if (x < 0.0) return NAN; // error("x less than 0 in routine gser");
*gamser=0.0;
return gln;
}
Expand All @@ -135,8 +135,8 @@ double gser( double *gamser, double a, double x){
return gln;
}
}
error("a too large, ITMAX too small in routine gser");
return -1;
// error("a too large, ITMAX too small in routine gser");
return NAN;
}
}

Expand Down Expand Up @@ -197,8 +197,8 @@ double invgammp( const double p, const double a) {
double afac = 0;
double lna1 = 0;
double eps = 1.e-8; //Accuracy is the square of EPS.
double gln = gammln(a);
if (a <= 0.) error("a must be pos in invgammap");
double gln = gammln(a);
if (a <= 0.) return NAN; //error("a must be pos in invgammap");
if (p >= 1.) return dmax(100.,a + 100.*sqrt(a));
if (p <= 0.) return 0.0;
if (a > 1.) {
Expand Down
5 changes: 4 additions & 1 deletion src/phyc/operator.c
View file
Expand Up @@ -114,7 +114,10 @@ bool operator_interval_scaler(Operator* op, double* logHR){

bool operator_scaler(Operator* op, double* logHR){
if(op->x != NULL){
size_t index = gsl_rng_uniform_int(op->rng, Parameters_count(op->x));
size_t index = 0;
if(Parameters_count(op->x) != 1){
index = gsl_rng_uniform_int(op->rng, Parameters_count(op->x));
}
Parameter* p = Parameters_at(op->x, index);
double scaler_scaleFactor = op->parameters[0];
double v = Parameter_value(p);
Expand Down
101 changes: 76 additions & 25 deletions src/phyc/sitemodel.c
View file
Expand Up @@ -34,7 +34,7 @@
#include <gsl/gsl_cdf.h>
#endif

static void _gamma_approx_quantile( SiteModel *sm );
static bool _gamma_approx_quantile( SiteModel *sm );
static void _calculate_rates_discrete( SiteModel *sm );

static double _get_rate( SiteModel *sm, const int index );
Expand Down Expand Up @@ -264,30 +264,44 @@ double _gamma_shape_derivative( SiteModel *sm, const double* ingrad ){
double shape = Parameters_value(sm->rates, 0);
size_t j = (variableCat == catCount ? 0 : 1); // ignore category 0 with r_0=0
double* temp = dvector(variableCat*2);
for (size_t i = 0; i < variableCat; i++, j++) {
const double eps = pow(DBL_EPSILON, 1.0/3.0);
#ifndef GSL_DISABLED
gsl_error_handler_t* handler = gsl_set_error_handler_off();
#endif
size_t i = 0;
for ( ; i < variableCat; i++, j++) {
double prob = (2.0*i + 1.0)/(2.0*variableCat);
double xp = shape*(1.0 + eps);
double xm = shape*(1.0 - eps);
#ifndef GSL_DISABLED
double plus = gsl_cdf_gamma_Qinv( prob, shape + sm->epsilon, 1.0/(shape + sm->epsilon) );
double minus = gsl_cdf_gamma_Qinv( prob, shape - sm->epsilon, 1.0/(shape - sm->epsilon) );
double plus = gsl_cdf_gamma_Qinv( prob, xp, 1.0/xp );
double minus = gsl_cdf_gamma_Qinv( prob, xm, 1.0/xm );
temp[i+variableCat] = gsl_cdf_gamma_Qinv( prob, shape, 1.0/shape);
#else
double plus = qgamma( prob, shape + sm->epsilon, shape + sm->epsilon );
double minus = qgamma( prob, shape - sm->epsilon, shape - sm->epsilon );
double plus = qgamma( prob, xp, xp );
double minus = qgamma( prob, xm, xm );
temp[i+variableCat] = qgamma( prob, shape, shape);
#endif
temp[i] = (plus - minus)/(2.0*sm->epsilon);

temp[i] = (plus - minus)/(xp-xm);
if(isnan(plus) || isnan(minus) || isnan(temp[i+variableCat])){
break;
}

derivsum += temp[i] * proportions[j];
sum += temp[i+variableCat] * proportions[j];
}

double shape_gradient = 0;
j = (variableCat == catCount ? 0 : 1); // ignore category 0 with r_0=0
for (size_t i = 0; i < variableCat; i++, j++) {
double prob = (2.0*i + 1.0)/(2.0*variableCat);
double deriv_rate = temp[i]/sum - (temp[i+variableCat]*derivsum)/sum/sum;
shape_gradient += ingrad[j] * deriv_rate * proportions[j];
#ifndef GSL_DISABLED
gsl_set_error_handler(handler);
#endif
double shape_gradient = NAN;
if(i == variableCat){
shape_gradient = 0.0;
j = (variableCat == catCount ? 0 : 1); // ignore category 0 with r_0=0
for (size_t i = 0; i < variableCat; i++, j++) {
double prob = (2.0*i + 1.0)/(2.0*variableCat);
double deriv_rate = temp[i]/sum - (temp[i+variableCat]*derivsum)/sum/sum;
shape_gradient += ingrad[j] * deriv_rate * proportions[j];
}
}
free(temp);
return shape_gradient;
Expand Down Expand Up @@ -421,6 +435,17 @@ void _weibull_gradient( SiteModel *sm, const double* ingrad, double* grad ){
}
}

bool _update_gamma_approx_quantile(SiteModel *sm){
if ( sm->need_update ) {
return _gamma_approx_quantile(sm);
}
return true;
}

bool _update_nothing(SiteModel *sm){
return true;
}

// (Gamma or/and Invariant) or one rate
// should not be used directly
SiteModel * new_SiteModel_with_parameters( const Parameters *params, Simplex* proportions, const size_t cat_count, distribution_t distribution, bool invariant, quadrature_t quad){
Expand Down Expand Up @@ -463,7 +488,8 @@ SiteModel * new_SiteModel_with_parameters( const Parameters *params, Simplex* pr
sm->cat_proportions = dvector(sm->cat_count);
sm->gradient = _no_gradient;
sm->derivative = _no_derivative;

sm->update = _update_nothing;

if (distribution == DISTRIBUTION_UNIFORM) {
sm->get_rate = _get_rate;
sm->get_proportion = _get_proportion;
Expand All @@ -484,6 +510,7 @@ SiteModel * new_SiteModel_with_parameters( const Parameters *params, Simplex* pr
sm->get_rate = _get_rate_gamma;
sm->get_proportion = _get_proportion_gamma;
sm->get_proportions = _get_proportions_gamma;
sm->update = _update_gamma_approx_quantile;
}

// Discrete gamma or prop of invariant or both
Expand Down Expand Up @@ -543,7 +570,7 @@ double icdf_weibull_1(double p, double k){
return pow(-log(1.0 - p), 1.0/k);
}

void _gamma_approx_quantile( SiteModel *sm ) {
bool _gamma_approx_quantile( SiteModel *sm ) {

double propVariable = 1.0;
int cat = (sm->invariant ? 1 : 0);
Expand Down Expand Up @@ -622,20 +649,19 @@ void _gamma_approx_quantile( SiteModel *sm ) {
sm->cat_rates[0] = 0;
sm->cat_rates[1] = 1.0 / sm->cat_proportions[1];
sm->need_update = false;
return;
free(quantiles);
return true;
}
// +G+I
else{
cat = 1;
sm->cat_proportions[0] = values[0];
// printf("%f %f\n", values[0], values[1]);
propVariable = values[1];
int gammaCat = sm->cat_count - 1;
for (int i = 0; i < gammaCat; i++) {
quantiles[i+cat] = (2.0 * i + 1.0) / (2.0 * gammaCat);
sm->cat_proportions[i+cat] = propVariable/gammaCat;
}
// exit(2);
}
}
// distribution without invariant and proportions not estimated
Expand All @@ -656,35 +682,56 @@ void _gamma_approx_quantile( SiteModel *sm ) {
if(sm->quadrature == QUADRATURE_QUANTILE_MEDIAN || sm->quadrature == QUADRATURE_DISCRETE ||
sm->quadrature == QUADRATURE_BETA || sm->quadrature == QUADRATURE_KUMARASWAMY){
sm->cat_rates[0] = 0;
size_t i = 0;
if(sm->distribution == DISTRIBUTION_GAMMA){
for (int i = 0; i < sm->cat_count - cat; i++) {
#ifndef GSL_DISABLED
gsl_error_handler_t* handler = gsl_set_error_handler_off();
for ( ; i < sm->cat_count - cat; i++) {
sm->cat_rates[i + cat] = gsl_cdf_gamma_Qinv( quantiles[i+cat], alpha, 1.0/alpha );
if(isnan(sm->cat_rates[i + cat])){
break;
}
}
gsl_set_error_handler(handler);
#else
for ( ; i < sm->cat_count - cat; i++) {
sm->cat_rates[i + cat] = qgamma( quantiles[i+cat], alpha, alpha );
#endif
if(isnan(sm->cat_rates[j + cat])){
break;
}
}
#endif
}
else if(sm->distribution == DISTRIBUTION_WEIBULL){
for (int i = 0; i < sm->cat_count - cat; i++) {
for ( ; i < sm->cat_count - cat; i++) {
// Unit mean Weibull
// sm->cat_rates[i + cat] = gsl_cdf_weibull_Qinv( sm->cat_proportions[i + cat], alpha, 1.0/exp(gammln(1.0 + 1.0/alpha)) );
// Fix lambda:=1
sm->cat_rates[i + cat] = icdf_weibull_1( quantiles[i+cat], alpha);
if(isnan(sm->cat_rates[i + cat])){
break;
}

}
}
#ifndef GSL_DISABLED
else if(sm->distribution == DISTRIBUTION_LOGNORMAL){
for (int i = 0; i < sm->cat_count - cat; i++) {
for ( ; i < sm->cat_count - cat; i++) {
sm->cat_rates[i + cat] = gsl_cdf_lognormal_Qinv( quantiles[i+cat], -alpha*alpha/2, alpha );
}
}
else if(sm->distribution == DISTRIBUTION_BETA){
for (int i = 0; i < sm->cat_count - cat; i++) {
for ( ; i < sm->cat_count - cat; i++) {
sm->cat_rates[i + cat] = gsl_cdf_beta_Qinv( quantiles[i+cat], alpha, Parameters_value(sm->rates, 1) );
}
}
#endif

if(i != sm->cat_count - cat){
free(quantiles);
sm->need_update = false;
return false;
}

if (sm->quadrature == QUADRATURE_BETA || sm->quadrature == QUADRATURE_KUMARASWAMY) {
for (int i = 0; i < sm->cat_count - cat; i++) {
Expand Down Expand Up @@ -728,6 +775,7 @@ void _gamma_approx_quantile( SiteModel *sm ) {
}
free(quantiles);
sm->need_update = false;
return true;
}

// Gamma distribution approximated using Laguerre quadrature
Expand Down Expand Up @@ -886,6 +934,7 @@ SiteModel * new_CATSiteModel_with_parameters( const Parameters *params, const s
sm->get_rate = _get_rate_cat;
sm->get_proportion = _get_proportion;
sm->get_proportions = _get_proportions;
sm->update = _update_nothing;

sm->integrate = false;
sm->need_update = true;
Expand Down Expand Up @@ -932,6 +981,7 @@ SiteModel * clone_SiteModel_with( const SiteModel *sm ){
newsm->get_proportion = sm->get_proportion;
newsm->get_proportions = sm->get_proportions;
newsm->get_site_category = sm->get_site_category;
newsm->update = sm->update;

newsm->distribution = sm->distribution;
newsm->invariant = sm->invariant;
Expand Down Expand Up @@ -985,6 +1035,7 @@ SiteModel * clone_SiteModel_with_parameters( const SiteModel *sm, Simplex* props
newsm->get_proportion = sm->get_proportion;
newsm->get_proportions = sm->get_proportions;
newsm->get_site_category = sm->get_site_category;
newsm->update = sm->update;

newsm->gradient = sm->gradient;
newsm->derivative = sm->derivative;
Expand Down
1 change: 1 addition & 0 deletions src/phyc/sitemodel.h
View file
Expand Up @@ -46,6 +46,7 @@ typedef struct SiteModel{

void (*set_rate)( struct SiteModel *, const int, const double );

bool (*update)( struct SiteModel * );
double (*get_rate)( struct SiteModel *, const int );
double (*get_proportion)( struct SiteModel *, const int );
double * (*get_proportions)( struct SiteModel * );
Expand Down
22 changes: 17 additions & 5 deletions src/phyc/treelikelihood.c
View file
Expand Up @@ -322,11 +322,15 @@ double* TreeLikelihood_gradient(Model *self){
}
double logP = tlk->calculate(tlk); // make sure it is updated
if (isnan(logP) || isinf(logP)) {
// return logP;
for(size_t i = 0; i < tlk->gradient_length; i++){
tlk->gradient[i] = NAN;
}
}
else{
update_upper_partials(tlk, Tree_root(tlk->tree), tlk->include_root_freqs);
TreeLikelihood_calculate_gradient(self, tlk->gradient);
tlk->update_upper = false;
}
update_upper_partials(tlk, Tree_root(tlk->tree), tlk->include_root_freqs);
TreeLikelihood_calculate_gradient(self, tlk->gradient);
tlk->update_upper = false;
}
return tlk->gradient;
}
Expand Down Expand Up @@ -1450,6 +1454,13 @@ double _calculate_simple( SingleTreeLikelihood *tlk ){
if( !tlk->update ){
return tlk->lk;
}

bool success = tlk->sm->update(tlk->sm);
if(!success){
tlk->lk = NAN;
return NAN;
}

if(Tree_is_time_mode(tlk->tree)){
Tree_update_heights(tlk->tree);
}
Expand Down Expand Up @@ -2124,10 +2135,11 @@ void update_upper_partials(SingleTreeLikelihood *tlk, Node* node, bool include_r

if(!Node_isroot(parent)){
int idMatrix = Node_id(parent);

// u_n = (P_p u_p) o (P_s p_s)
tlk->update_partials(tlk, tlk->upper_partial_indexes[idNode], tlk->upper_partial_indexes[Node_id(parent)], idMatrix, idSibling, idSibling);
}
else{
// u_n = P_s u_s
tlk->update_partials(tlk, tlk->upper_partial_indexes[idNode], idSibling, idSibling, -1, -1);
if(include_root_freqs){
const double* freqs = tlk->get_root_frequencies(tlk);
Expand Down

0 comments on commit 9c94f57

Please sign in to comment.