C++ Branch::set_len方法代码示例

/**
 * Calculates the log likelihood (log probability of model given the
 * data) and gradient. If the calc_ll flag is FALSE, the
 * log-likelihood is not calculated and 0.0 is returned. If the
 * calc_grad flag is FALSE the gradient is not calculated and the
 * model partial derivatives are not updated.
 */
double bkgd_evo_mdl_calc_ll(Model *mdl, int calc_ll, int calc_grad) {
  BkgdEvoMdlParam param, ll_deriv;
  Branch *br;
  ColType *cltype;
  BkgdEvoMdlConfig *conf;
  BkgdEvoMdlData *data;
  double ll;
  ColTypeStat *clstat, *cons_clstat;
  long i;
  int j;
   
  conf = mdl->config;
  data = mdl->data;

  if(data == NULL) {
    g_error("bkgd_evo_mdl_calc_ll: data must be set before "
	    "likelihood can be calculated\n");
  }

  param_set_zero(&param);
  param_set_zero(&ll_deriv);
  set_bkgd_evo_param_from_model(mdl, conf, &param);

  ll = 0.0;

  clstat  = g_new(ColTypeStat, conf->n_cltype);
  cons_clstat = NULL;
  for(i = 0; i < conf->n_cltype; i++) {
    clstat[i].n = 0;
    clstat[i].ttl_prob = 0.0;
    clstat[i].ttl_prob_double = 0.0;
    clstat[i].ttl_prob_single = 0.0;
    clstat[i].ttl_ll = 0.0;
  }

  for(i = 0; i < data->n_bin; i++) {
    /* Combine non-exonic and exonic B vals and rescale by new
     * deleterious rate estimates. Use precomputed log B values so
     * that rescaling is more efficient (do not need to take logs
     * first and then re-exponentiate)
     */ 
    param.B = exp(data->bin[i].lB_ex  * param.u_ex_scale + 
		  data->bin[i].lB_nex * param.u_nex_scale);

    if(param.B < 0.0 || param.B > 1.0 || isnan(param.B) || isinf(param.B)) {
      g_error("bkgd_evo_mdl_calc_ll: invalid B (%g)\n"
	      "B_ex=%g, B_nex=%g, u_ex_scale=%g, u_nex_scale=%g\n",
	      param.B, data->bin[i].B_ex, data->bin[i].B_nex,
	      param.u_ex_scale, param.u_nex_scale);
    }

    /* calculate mutation rate for current bin */
    bkgd_evo_mdl_calc_mu(&param, mdl, &data->bin[i]);

    /* calculate probability HC coalescent predates HC/G speciation */
    if(param.T_hcg == 0.0) {
      param.k_hcg = 1.0;
    } else {
      param.k_hcg = exp(-0.5 * param.T_hcg / (param.B * param.N_hc));

      if(param.k_hcg > 1.0 || param.k_hcg < 0.0) {
	g_error("bkgd_evo_mdl_calc_ll: invalid probability of HC "
		"coalescent predating gorilla speciation: k_hcg=%g\n",
		param.k_hcg);
      }
    }

    /* calculate branch lengths and substitution probs for current bin */
    for(j = 0; j < conf->n_branch; j++) {
      br = &conf->branches[j];

      br->set_len(br, &param);

      if(br->len < 0.0 || isnan(br->len) || isinf(br->len)) {
	g_warning("bkgd_evo_mdl_calc_ll: invalid length (%g) for branch %s",
		br->len, br->name);
	br->len = BRANCH_LEN_SMALL;
      }

      branch_set_prob(br, &param, conf);
    }

    /* calc prob of observing each observed column type for bin */
    for(j = 0; j < conf->n_cltype; j++) {
      cltype = conf->cltypes[j];

      /* do only non-conserved columns first */
      if(cltype->subst_type == SUBST_TYPE_CONSERVED) {
	cons_clstat = &clstat[j];
	continue;
      }

      cltype->set_n(cltype, &data->bin[i]);
//.........这里部分代码省略.........