text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include <gsl/gsl_integration.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>

#include "ccl.h"

/*----- ROUTINE: dc_NakamuraSuto -----
INPUT: cosmology, scale factor
TASK: Computes the peak threshold: delta_c(z) assuming LCDM.
Cosmology dependence of the critical linear density according to the spherical-collapse model.
Fitting function from Nakamura & Suto (1997; arXiv:astro-ph/9710107).
*/
double dc_NakamuraSuto(ccl_cosmology *cosmo, double a, int *status){

  double Om_mz = ccl_omega_x(cosmo, a, ccl_species_m_label, status);
  double dc0 = (3./20.)*pow(12.*M_PI,2./3.);
  double dc = dc0*(1.+0.012299*log10(Om_mz));

  return dc;

}

/*----- ROUTINE: Dv_BryanNorman -----
INPUT: cosmology, scale factor
TASK: Computes the virial collapse density contrast with respect to the matter density assuming LCDM.
Cosmology dependence of the virial collapse density according to the spherical-collapse model
Fitting function from Bryan & Norman (1998; arXiv:astro-ph/9710107)
*/
double Dv_BryanNorman(ccl_cosmology *cosmo, double a, int *status){

  double Om_mz = ccl_omega_x(cosmo, a, ccl_species_m_label, status);
  double x = Om_mz-1.;
  double Dv0 = 18.*pow(M_PI,2);
  double Dv = (Dv0+82.*x-39.*pow(x,2))/Om_mz;

  return Dv;
}

/*----- ROUTINE: r_delta -----
INPUT: cosmology, halo mass, scale factor, halo overdensity
TASK: Computes comoving halo radius assuming the overdensity criteria
*/
double r_delta(ccl_cosmology *cosmo, double halomass, double a, double odelta, int *status){

  double rho_matter = ccl_rho_x(cosmo, 1., 1, 1, status);

  return pow(halomass*3.0/(4.0*M_PI*rho_matter*odelta),1.0/3.0);

}

// This checks to make sure all necessary halo mass function parameters have been set-up,
// as well as associated splines.
void ccl_cosmology_compute_hmfparams(ccl_cosmology *cosmo, int *status) {
  if(cosmo->computed_hmfparams)
    return;

  gsl_spline* alphahmf = NULL;
  gsl_spline* betahmf = NULL;
  gsl_spline* gammahmf = NULL;
  gsl_spline* phihmf = NULL;
  gsl_spline* etahmf = NULL;

  // declare parameter splines on case-by-case basis
  switch(cosmo->config.mass_function_method) {
  case ccl_tinker10:{
    double delta[9] = {200.0, 300.0, 400.0, 600.0, 800.0, 1200.0, 1600.0, 2400.0, 3200.0};
    double lgdelta[9];
    double alpha[9] = {0.368, 0.363, 0.385, 0.389, 0.393, 0.365, 0.379, 0.355, 0.327};
    double beta[9] = {0.589, 0.585, 0.544, 0.543, 0.564, 0.623, 0.637, 0.673, 0.702};
    double gamma[9] ={0.864, 0.922, 0.987, 1.09, 1.20, 1.34, 1.50, 1.68, 1.81};
    double phi[9] = {-0.729, -0.789, -0.910, -1.05, -1.20, -1.26, -1.45, -1.50, -1.49};
    double eta[9] = {-0.243, -0.261, -0.261, -0.273, -0.278, -0.301, -0.301, -0.319, -0.336};
    int nd = 9;
    int i;

    for(i=0; i<nd; i++) {
      lgdelta[i] = log10(delta[i]);
    }

    if (*status == 0) {
      alphahmf = gsl_spline_alloc(cosmo->spline_params.D_SPLINE_TYPE, nd);
      if (alphahmf == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error allocating alpha(D) spline\n");
      }
    }
    if (*status == 0) {
      *status = gsl_spline_init(alphahmf, lgdelta, alpha, nd);
      if (*status) {
        *status = CCL_ERROR_SPLINE ;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating alpha(D) spline\n");
      }
    }

    if (*status == 0) {
      betahmf = gsl_spline_alloc(cosmo->spline_params.D_SPLINE_TYPE, nd);
      if (betahmf == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error allocating beta(D) spline\n");
      }
    }
    if (*status == 0) {
      *status = gsl_spline_init(betahmf, lgdelta, beta, nd);
      if (*status) {
        *status = CCL_ERROR_SPLINE;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating beta(D) spline\n");
      }
    }

    if (*status == 0) {
      gammahmf = gsl_spline_alloc(cosmo->spline_params.D_SPLINE_TYPE, nd);
      if (gammahmf == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error allocating gamma(D) spline\n");
      }
    }
    if (*status == 0) {
      *status = gsl_spline_init(gammahmf, lgdelta, gamma, nd);
      if (*status) {
        *status = CCL_ERROR_SPLINE;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating gamma(D) spline\n");
      }
    }

    if (*status == 0) {
      phihmf = gsl_spline_alloc(cosmo->spline_params.D_SPLINE_TYPE, nd);
      if (phihmf == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error allocating phi(D) spline\n");
      }
    }
    if (*status == 0) {
      *status = gsl_spline_init(phihmf, lgdelta, phi, nd);
      if (*status) {
        *status = CCL_ERROR_SPLINE;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating phi(D) spline\n");
      }
    }

    if (*status == 0) {
      etahmf = gsl_spline_alloc(cosmo->spline_params.D_SPLINE_TYPE, nd);
      if (etahmf == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error allocating eta(D) spline\n");
      }
    }
    if (*status == 0) {
      *status = gsl_spline_init(etahmf, lgdelta, eta, nd);
      if (*status) {
        *status = CCL_ERROR_SPLINE;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating eta(D) spline\n");
      }
    }

    if (*status == 0) {
      cosmo->data.alphahmf = alphahmf;
      cosmo->data.betahmf = betahmf;
      cosmo->data.gammahmf = gammahmf;
      cosmo->data.phihmf = phihmf;
      cosmo->data.etahmf = etahmf;
      cosmo->computed_hmfparams = true;
    }
    break;
  }
  case ccl_tinker:{
    double delta[9] = {200.0, 300.0, 400.0, 600.0, 800.0, 1200.0, 1600.0, 2400.0, 3200.0};
    double lgdelta[9];
    double alpha[9] = {0.186, 0.200, 0.212, 0.218, 0.248, 0.255, 0.260, 0.260, 0.260};
    double beta[9] = {1.47, 1.52, 1.56, 1.61, 1.87, 2.13, 2.30, 2.53, 2.66};
    double gamma[9] ={2.57, 2.25, 2.05, 1.87, 1.59, 1.51, 1.46, 1.44, 1.41};
    double phi[9] = {1.19, 1.27, 1.34, 1.45, 1.58, 1.80, 1.97, 2.24, 2.44};
    int nd = 9;
    int i;

    for(i=0; i<nd; i++) {
      lgdelta[i] = log10(delta[i]);
    }

    if (*status == 0) {
      alphahmf = gsl_spline_alloc(cosmo->spline_params.D_SPLINE_TYPE, nd);
      if (alphahmf == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error allocating alpha(D) spline\n");
      }
    }
    if (*status == 0) {
      *status = gsl_spline_init(alphahmf, lgdelta, alpha, nd);
      if (*status) {
        *status = CCL_ERROR_SPLINE;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating alpha(D) spline\n");
      }
    }

    if (*status == 0) {
      betahmf = gsl_spline_alloc(cosmo->spline_params.D_SPLINE_TYPE, nd);
      if (betahmf == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error allocating beta(D) spline\n");
      }
    }
    if (*status == 0) {
      *status = gsl_spline_init(betahmf, lgdelta, beta, nd);
      if (*status) {
        *status = CCL_ERROR_SPLINE;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating beta(D) spline\n");
      }
    }

    if (*status == 0) {
      gammahmf = gsl_spline_alloc(cosmo->spline_params.D_SPLINE_TYPE, nd);
      if (gammahmf == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error allocating gamma(D) spline\n");
      }
    }
    if (*status == 0) {
      *status = gsl_spline_init(gammahmf, lgdelta, gamma, nd);
      if (*status) {
        *status = CCL_ERROR_SPLINE;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating gamma(D) spline\n");
      }
    }

    if (*status == 0) {
      phihmf = gsl_spline_alloc(cosmo->spline_params.D_SPLINE_TYPE, nd);
      if (phihmf == NULL) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error allocating phi(D) spline\n");
      }
    }
    if (*status == 0) {
      *status = gsl_spline_init(phihmf, lgdelta, phi, nd);
      if (*status) {
        *status = CCL_ERROR_SPLINE;
        ccl_cosmology_set_status_message(
          cosmo,
          "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating phi(D) spline\n");
      }
    }

    if (*status == 0) {
      cosmo->data.alphahmf = alphahmf;
      cosmo->data.betahmf = betahmf;
      cosmo->data.gammahmf = gammahmf;
      cosmo->data.phihmf = phihmf;
      cosmo->computed_hmfparams = true;
    }

    break;
  }
  default:
    // Error message could go here if we decide to make this public facing.
    // Currently not accessible from the API though.
    break;
  }

  if (*status) {
    gsl_spline_free(alphahmf);
    gsl_spline_free(betahmf);
    gsl_spline_free(gammahmf);
    gsl_spline_free(phihmf);
    gsl_spline_free(etahmf);
  }
}

//TODO: some of these are unused, many are included in ccl.h

/*----- ROUTINE: ccl_massfunc_f -----
INPUT: cosmology+parameters, a halo mass, and scale factor
TASK: Outputs fitting function for use in halo mass function calculation;
  currently only supports:
    ccl_tinker (arxiv 0803.2706 )
    ccl_tinker10 (arxiv 1001.3162 )
    ccl_angulo (arxiv 1203.3216 )
    ccl_watson (arxiv 1212.0095 )
    ccl_shethtormen (arxiv 9901122)
-*/
static double massfunc_f(ccl_cosmology *cosmo, double halomass, double a, double odelta, int *status)
{
  double fit_A, fit_a, fit_b, fit_c, fit_d, fit_p, overdensity_delta;
  double Omega_m_a;
  double delta_c_Tinker, nu;

  double sigma=ccl_sigmaM(cosmo, log10(halomass), a, status);
  int gslstatus;

  switch(cosmo->config.mass_function_method) {

  // Equation (10) in arxiv: 9901122
  // Note that Sheth & Tormen (1999) use nu=(dc/sigma)^2 whereas we use nu=dc/sigma
  case ccl_shethtormen:

    // Check if odelta is outside the interpolated range
    if (odelta != Dv_BryanNorman(cosmo, a, status)) {
      *status = CCL_ERROR_HMF_DV;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_massfunc.c: massfunc_f(): Sheth-Tormen called with not virial Delta_v\n");
      return NAN;
    }

    // ST mass function fitting parameters
    fit_A = 0.21616;
    fit_p = 0.3;
    fit_a = 0.707;

    // nu = delta_c(z) / sigma(M)
    nu = dc_NakamuraSuto(cosmo, a, status)/ccl_sigmaM(cosmo, log10(halomass), a, status);

    return nu*fit_A*(1.+pow(fit_a*pow(nu,2),-fit_p))*exp(-fit_a*pow(nu,2)/2.);

  case ccl_tinker:

    // Check if odelta is outside the interpolated range
    if ((odelta < 200) || (odelta > 3200)) {
      *status = CCL_ERROR_HMF_INTERP;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_massfunc.c: massfunc_f(): Tinker 2008 only "
        "supported in range of Delta = 200 to Delta = 3200.\n");
      return NAN;
    }

    if (!cosmo->computed_hmfparams) {
      *status = CCL_ERROR_HMF_INIT;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_massfunc.c: massfunc_f(): mass function parameters splines have not been computed!");
      return NAN;
    }

    // Compute HMF parameter (alpha, beta, gamma, phi) splines if they haven't
    // been computed already
    gslstatus = gsl_spline_eval_e(cosmo->data.alphahmf, log10(odelta), NULL, &fit_A);
    gslstatus |= gsl_spline_eval_e(cosmo->data.betahmf, log10(odelta), NULL, &fit_a);
    gslstatus |= gsl_spline_eval_e(cosmo->data.gammahmf, log10(odelta), NULL, &fit_b);
    gslstatus |= gsl_spline_eval_e(cosmo->data.phihmf, log10(odelta), NULL, &fit_c);
    fit_d = pow(10, -1.0*pow(0.75 / log10(odelta / 75.0), 1.2));

    fit_A = fit_A*pow(a, 0.14);
    fit_a = fit_a*pow(a, 0.06);
    fit_b = fit_b*pow(a, fit_d);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_massfunc.c: massfunc_f():");
      *status |= gslstatus;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_massfunc.c: massfunc_f(): interpolation error for Tinker MF\n");
      return NAN;
    }
    return fit_A*(pow(sigma/fit_b,-fit_a)+1.0)*exp(-fit_c/sigma/sigma);
    break;
  case ccl_tinker10:
    // this version uses f(nu) parameterization from Eq. 8 in Tinker et al. 2010
    // use this for consistency with Tinker et al. 2010 fitting function for halo bias

    // Check if odelta is outside the interpolated range
    if ((odelta < 200) || (odelta > 3200)) {
      *status = CCL_ERROR_HMF_INTERP;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_massfunc.c: massfunc_f(): Tinker 2010 only "
        "supported in range of Delta = 200 to Delta = 3200.\n");
      return 0;
    }

    //critical collapse overdensity assumed in this model
    delta_c_Tinker = 1.686;
    nu = delta_c_Tinker/(sigma);

    if (!cosmo->computed_hmfparams) {
      *status = CCL_ERROR_HMF_INIT;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_massfunc.c: massfunc_f(): mass function parameters splines have not been computed!");
      return NAN;
    }

    gslstatus = gsl_spline_eval_e(cosmo->data.alphahmf, log10(odelta), NULL, &fit_A); //alpha in Eq. 8
    gslstatus |= gsl_spline_eval_e(cosmo->data.etahmf, log10(odelta), NULL, &fit_a); //eta in Eq. 8
    gslstatus |= gsl_spline_eval_e(cosmo->data.betahmf, log10(odelta), NULL, &fit_b); //beta in Eq. 8
    gslstatus |= gsl_spline_eval_e(cosmo->data.gammahmf, log10(odelta), NULL, &fit_c); //gamma in Eq. 8
    gslstatus |= gsl_spline_eval_e(cosmo->data.phihmf, log10(odelta), NULL, &fit_d); //phi in Eq. 8;

    fit_a *=pow(a, -0.27);
    fit_b *=pow(a, -0.20);
    fit_c *=pow(a, 0.01);
    fit_d *=pow(a, 0.08);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_massfunc.c: massfunc_f():");
      *status |= gslstatus;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_massfunc.c: massfunc_f(): interpolation error for Tinker 2010 MF\n");
      return NAN;
    }
    return nu*fit_A*(1.+pow(fit_b*nu,-2.*fit_d))*pow(nu, 2.*fit_a)*exp(-0.5*fit_c*nu*nu);
    break;

  case ccl_watson:
    if(odelta!=200.) {
      *status = CCL_ERROR_HMF_INTERP;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_massfunc.c: massfunc_f(): Watson HMF only supported for Delta = 200.\n");
      return NAN;
    }
    // these parameters from: Angulo et al 2012 (arxiv 1203.3216 )
    Omega_m_a = ccl_omega_x(cosmo, a, ccl_species_m_label,status);
    fit_A = Omega_m_a*(0.990*pow(a,3.216)+0.074);
    fit_a = Omega_m_a*(5.907*pow(a,3.599)+2.344);
    fit_b = Omega_m_a*(3.136*pow(a,3.058)+2.349);
    fit_c = 1.318;

    return fit_A*(pow(sigma/fit_b,-fit_a)+1.0)*exp(-fit_c/sigma/sigma);

  case ccl_angulo:
    if(odelta!=200.) {
      *status = CCL_ERROR_HMF_INTERP;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_massfunc.c: massfunc_f(): Angulo HMF only supported for Delta = 200.\n");
      return NAN;
    }
    // these parameters from: Watson et al 2012 (arxiv 1212.0095 )
    fit_A = 0.201;
    fit_a = 2.08;
    fit_b = 1.7;
    fit_c = 1.172;

    return fit_A*pow( (fit_a/sigma)+1.0, fit_b)*exp(-fit_c/sigma/sigma);

  default:
    *status = CCL_ERROR_MF;
    ccl_cosmology_set_status_message(cosmo ,
	    "ccl_massfunc.c: massfunc_f(): Unknown or non-implemented mass function method: %d \n",
	    cosmo->config.mass_function_method);
    return NAN;
  }
}

static double ccl_halo_b1(ccl_cosmology *cosmo, double halomass, double a, double odelta, int *status)
{
  double fit_A, fit_B, fit_C, fit_a, fit_b, fit_c, fit_p, overdensity_delta, y;
  double delta_c_Tinker, nu;
  double sigma=ccl_sigmaM(cosmo,log10(halomass),a, status);
  switch(cosmo->config.mass_function_method) {

  // Equation (12) in  arXiv: 9901122
  // Derived using the peak-background split applied to the mass function in the same paper
  // Note that Sheth & Tormen (1999) use nu=(dc/sigma)^2 whereas we use nu=dc/sigma
  case ccl_shethtormen:

    // Check if Delta_v is the virial Delta_v
    if (odelta != Dv_BryanNorman(cosmo, a, status)) {
      *status = CCL_ERROR_HMF_DV;
      ccl_cosmology_set_status_message(
        cosmo, "ccl_massfunc.c: halo_b1(): Sheth-Tormen called with not virial Delta_v\n");
      return NAN;
    }

    // ST bias fitting parameters (which are the same as for the mass function)
    fit_p = 0.3;
    fit_a = 0.707;

    // Cosmology dependent delta_c and nu
    double delta_c = dc_NakamuraSuto(cosmo, a, status);
    nu = delta_c/ccl_sigmaM(cosmo, log10(halomass), a, status);

    return 1.+(fit_a*pow(nu,2)-1.+2.*fit_p/(1.+pow(fit_a*pow(nu,2),fit_p)))/delta_c;

    //this version uses b(nu) parameterization, Eq. 6 in Tinker et al. 2010
    // use this for consistency with Tinker et al. 2010 fitting function for halo bias
  case ccl_tinker10:
    y = log10(odelta);
    //critical collapse overdensity assumed in this model
    delta_c_Tinker = 1.686;
    //peak height - note that this factorization is incorrect for e.g. massive neutrino cosmologies
    nu = delta_c_Tinker/(sigma);
    // Table 2 in https://arxiv.org/pdf/1001.3162.pdf
    fit_A = 1.0 + 0.24*y*exp(-pow(4./y,4.));
    fit_a = 0.44*y-0.88;
    fit_B = 0.183;
    fit_b = 1.5;
    fit_C = 0.019+0.107*y+0.19*exp(-pow(4./y,4.));
    fit_c = 2.4;

    return 1.-fit_A*pow(nu,fit_a)/(pow(nu,fit_a)+pow(delta_c_Tinker,fit_a))+fit_B*pow(nu,fit_b)+fit_C*pow(nu,fit_c);
    break;

  default:
    *status = CCL_ERROR_MF;
    ccl_cosmology_set_status_message(cosmo ,
	    "ccl_massfunc.c: ccl_halo_b1(): No b(M) fitting function implemented for mass_function_method: %d \n",
      cosmo->config.mass_function_method);
    return 0;
  }
}

void ccl_cosmology_compute_sigma(ccl_cosmology *cosmo, int *status)
{
  if(cosmo->computed_sigma)
    return;

  // create linearly-spaced values of the mass.
  int nm = cosmo->spline_params.LOGM_SPLINE_NM;
  double *m = NULL;
  double *y = NULL;
  double smooth_radius;
  double na, nb;

  m = ccl_linear_spacing(cosmo->spline_params.LOGM_SPLINE_MIN, cosmo->spline_params.LOGM_SPLINE_MAX, nm);
  if (m == NULL ||
      (fabs(m[0]-cosmo->spline_params.LOGM_SPLINE_MIN)>1e-5) ||
      (fabs(m[nm-1]-cosmo->spline_params.LOGM_SPLINE_MAX)>1e-5) ||
      (m[nm-1]>10E17)) {
    *status = CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo,"ccl_cosmology_compute_sigmas(): Error creating linear spacing in m\n");
  }

  if (*status == 0) {
    // create space for y, to be filled with sigma and dlnsigma_dlogm
    y = malloc(sizeof(double)*nm);
    if (y == NULL) {
      *status = CCL_ERROR_MEMORY;
    }
  }

  // start up of GSL pointers
  int gslstatus = 0;
  gsl_spline *logsigma = NULL;
  gsl_spline *dlnsigma_dlogm = NULL;

  // fill in sigma, if no errors have been triggered at this time.
  if (*status == 0) {
    for (int i=0; i<nm; i++) {
      smooth_radius = ccl_massfunc_m2r(cosmo, pow(10,m[i]), status);
      y[i] = log(ccl_sigmaR(cosmo, smooth_radius, 1., status));
    }
    logsigma = gsl_spline_alloc(cosmo->spline_params.M_SPLINE_TYPE, nm);
    if (logsigma == NULL) {
      *status = CCL_ERROR_MEMORY;
    }
  }

  if (*status == 0) {
    gslstatus = gsl_spline_init(logsigma, m, y, nm);
    if (gslstatus != GSL_SUCCESS) {
      *status = CCL_ERROR_SPLINE ;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_massfunc.c: ccl_cosmology_compute_sigma(): Error creating sigma(M) spline\n");
    }
  }

  // again, making splines assuming nothing bad has happened to this point
  if (*status == 0 ) {
    for (int i=0; i<nm; i++) {
      if(i==0) {
        gslstatus |= gsl_spline_eval_e(logsigma, m[i], NULL,&na);
        gslstatus |= gsl_spline_eval_e(logsigma, m[i]+cosmo->spline_params.LOGM_SPLINE_DELTA/2., NULL,&nb);
        y[i] = 2.*(na-nb)*y[i] / cosmo->spline_params.LOGM_SPLINE_DELTA;
      }
      else if (i==nm-1) {
        gslstatus |= gsl_spline_eval_e(logsigma, m[i]-cosmo->spline_params.LOGM_SPLINE_DELTA/2., NULL,&na);
        gslstatus |= gsl_spline_eval_e(logsigma, m[i], NULL,&nb);
        y[i] = 2.*(na-nb)*y[i] / cosmo->spline_params.LOGM_SPLINE_DELTA;
      }
      else {
        gslstatus |= gsl_spline_eval_e(logsigma, m[i]-cosmo->spline_params.LOGM_SPLINE_DELTA/2., NULL,&na);
        gslstatus |= gsl_spline_eval_e(logsigma, m[i]+cosmo->spline_params.LOGM_SPLINE_DELTA/2., NULL,&nb);
        y[i] = (na-nb) / cosmo->spline_params.LOGM_SPLINE_DELTA;
      }
    }

    if(gslstatus != GSL_SUCCESS ) {
      ccl_raise_gsl_warning(
        gslstatus, "ccl_massfunc.c: ccl_cosmology_compute_sigma():");
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_massfunc.c: ccl_cosmology_compute_sigma(): "
        "Error evaluating grid points for dlnsigma/dlogM spline\n");
    }
  }

  if (*status == 0) {
    dlnsigma_dlogm = gsl_spline_alloc(cosmo->spline_params.M_SPLINE_TYPE, nm);
    if (dlnsigma_dlogm == NULL) {
      *status = CCL_ERROR_MEMORY;
    }
  }

  if (*status == 0) {
    gslstatus = gsl_spline_init(dlnsigma_dlogm, m, y, nm);
    if (gslstatus != GSL_SUCCESS) {
      *status = CCL_ERROR_SPLINE ;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_massfunc.c: ccl_cosmology_compute_sigma(): Error creating dlnsigma/dlogM spline\n");
    }
  }

  if (*status == 0) {
    cosmo->data.logsigma = logsigma;
    cosmo->data.dlnsigma_dlogm = dlnsigma_dlogm;
    cosmo->computed_sigma = true;
  } else {
    gsl_spline_free(logsigma);
    gsl_spline_free(dlnsigma_dlogm);
  }

  free(m);
  free(y);
}

/*----- ROUTINE: ccl_massfunc -----
INPUT: ccl_cosmology * cosmo, double halo mass in units of Msun, double scale factor
TASK: returns halo mass function as dn/dlog10(m) in comoving Msun^-1 Mpc^-3 (haloes per mass interval per volume)
*/
double ccl_massfunc(ccl_cosmology *cosmo, double halomass, double a, double odelta, int *status)
{
  if (cosmo->params.N_nu_mass>0){
	  *status = CCL_ERROR_NOT_IMPLEMENTED;
	  ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Support for the halo mass function in cosmologies with massive neutrinos is not yet implemented.\n");
	  return NAN;
  }

  double logmass;
  logmass = log10(halomass);
  if (logmass > cosmo->spline_params.LOGM_SPLINE_MAX || logmass < cosmo->spline_params.LOGM_SPLINE_MIN){
    *status = CCL_ERROR_HMF_INTERP;
    ccl_cosmology_set_status_message(cosmo, "ccl_massfunc(): The specified halo mass is outside of the range.");
    return NAN;
  }

  if (fabs(cosmo->params.mu_0)>1e-14 || fabs(cosmo->params.sigma_0)>1e-14){
	  *status = CCL_ERROR_NOT_IMPLEMENTED;
	  strcpy(cosmo->status_message,"ccl_massfunc.c: ccl_massfunc(): The halo mass funcion is not implemented the mu / Sigma modified gravity parameterisation.\n");
	  return NAN;
  }

  double f, rho_m;

  rho_m = ccl_constants.RHO_CRITICAL*cosmo->params.Omega_m*cosmo->params.h*cosmo->params.h;
  f=massfunc_f(cosmo,halomass,a,odelta,status);

  return f*rho_m*ccl_dlnsigM_dlogM(cosmo,logmass,status)/halomass;
}

/*----- ROUTINE: ccl_halob1 -----
INPUT: ccl_cosmology * cosmo, double halo mass in units of Msun, double scale factor
TASK: returns dimensionless linear halo bias
*/
double ccl_halo_bias(ccl_cosmology *cosmo, double halomass, double a, double odelta, int *status)
{
  if (cosmo->params.N_nu_mass>0){
	  *status = CCL_ERROR_NOT_IMPLEMENTED;
	  ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Support for the halo bias in cosmologies with massive neutrinos is not yet implemented.\n");
	  return NAN;
  }

  if (fabs(cosmo->params.mu_0)>1e-14 || fabs(cosmo->params.sigma_0)>1e-14){
	  *status = CCL_ERROR_NOT_IMPLEMENTED;
	  strcpy(cosmo->status_message,"ccl_massfunc.c: ccl_halobias(): The halo bias is not implemented the mu / Sigma modified gravity parameterisation.\n");
	  return NAN;
  }

  double f;
  f = ccl_halo_b1(cosmo,halomass,a,odelta, status);
  return f;
}
/*---- ROUTINE: ccl_massfunc_m2r -----
INPUT: ccl_cosmology * cosmo, halomass in units of Msun
TASK: takes halo mass and converts to halo radius
  in units of Mpc.
*/
double ccl_massfunc_m2r(ccl_cosmology *cosmo, double halomass, int *status)
{
  double rho_m, smooth_radius;

  // Comoving matter density
  //rho_m = ccl_constants.RHO_CRITICAL*cosmo->params.Omega_m*cosmo->params.h*cosmo->params.h;
  rho_m = ccl_rho_x(cosmo, 1., ccl_species_m_label, 1, status);

  smooth_radius = pow((3.0*halomass) / (4*M_PI*rho_m), (1.0/3.0));

  return smooth_radius;
}

/*----- ROUTINE: ccl_sigma_M -----
INPUT: ccl_cosmology * cosmo, double halo mass in units of Msun, double scale factor
TASK: returns sigma from the sigmaM interpolation. Also computes the sigma interpolation if
necessary.
*/
double ccl_sigmaM(ccl_cosmology *cosmo, double log_halomass, double a, int *status)
{
  double sigmaM;
  // Check if sigma has already been calculated
  if (!cosmo->computed_sigma) {
    *status = CCL_ERROR_SIGMA_INIT;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_massfunc.c: ccl_sigmaM(): linear power spctrum has not been computed!");
    return NAN;
  }

  double lgsigmaM;

  int gslstatus = gsl_spline_eval_e(cosmo->data.logsigma, log_halomass, NULL, &lgsigmaM);

  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_massfunc.c: ccl_sigmaM():");
    *status |= gslstatus;
  }

  // Interpolate to get sigma
  sigmaM = exp(lgsigmaM)*ccl_growth_factor(cosmo, a, status);
  return sigmaM;
}

/*----- ROUTINE: ccl_dlnsigM_dlogM -----
INPUT: ccl_cosmology *cosmo, double halo mass in units of Msun
TASK: returns the value of the derivative of ln(sigma^-1) with respect to log10 in halo mass.
*/
double ccl_dlnsigM_dlogM(ccl_cosmology *cosmo, double log_halomass, int *status)
{
  // Check if sigma has already been calculated
  if (!cosmo->computed_sigma) {
    *status = CCL_ERROR_SIGMA_INIT;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_massfunc.c: ccl_sigmaM(): linear power spctrum has not been computed!");
    return NAN;
  }

  double dlsdlgm;
  int gslstatus = gsl_spline_eval_e(cosmo->data.dlnsigma_dlogm,
				    log_halomass, NULL, &dlsdlgm);
  if(gslstatus) {
    ccl_raise_gsl_warning(gslstatus, "ccl_massfunc.c: ccl_dlnsigM_dlogM():");
    *status |= gslstatus;
  }
  return dlsdlgm;
}

max_stars_repo_name: benediktdiemer/CCL

max_issues_repo_name: benediktdiemer/CCL

max_forks_repo_name: benediktdiemer/CCL


text:
/* -*- c++ -*- */
/*
 * Copyright 2015 Free Software Foundation, Inc.
 *
 * This is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published
 * by the Free Software Foundation; either version 3, or (at your
 * option) any later version.
 *
 * This software is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this software; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

#ifndef INCLUDED_fec_mtrx_impl_H
#define INCLUDED_fec_mtrx_impl_H

#include <string>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gnuradio/fec/fec_mtrx.h>

namespace gr {
  namespace fec {
    namespace code {

      class fec_mtrx_impl : public fec_mtrx
      {
      protected:
        //! Constructor
        fec_mtrx_impl();

        //! Codeword length n
        unsigned int d_n;

        //! Information word length k
        unsigned int d_k;

        //! Number of rows in the matrix read in from alist file
        unsigned int d_num_rows;

        //! Number of columns in the matrix read in from alist file
        unsigned int d_num_cols;

        //! GSL matrix structure for the parity check matrix
        matrix_sptr d_H_sptr;

        //! Flag for whether or not the parity bits come first or last
        bool d_par_bits_last;

      public:
        //! Returns the parity check matrix H (needed by decoder)
        const gsl_matrix *H() const;

        //!Get the codeword length n
        unsigned int n() const;

        //! Get the information word length k
        unsigned int k() const;

        //! Subtract matrices using mod2 operations
        void add_matrices_mod2(gsl_matrix *result,
                               const gsl_matrix *,
                               const gsl_matrix *) const;

        //! Multiply matrices using mod2 operations
        void mult_matrices_mod2(gsl_matrix *result,
                                const gsl_matrix *,
                                const gsl_matrix *) const;

        //! Invert a square matrix using mod2 operations
        gsl_matrix *calc_inverse_mod2(const gsl_matrix *) const;

        /*!
         * \brief Get Boolean for whether or not parity bits come first or last
         * \details
         * The decoder will need to know if the parity bits are
         * coming first or last
         */
        bool parity_bits_come_last() const;

        virtual ~fec_mtrx_impl();
      };
    }
  }
}

#endif /* INCLUDED_fec_mtrx_impl_H */

ext: h


text:
// demo3d.c

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <time.h>
#include <string.h>
#include <assert.h>
#include <stdbool.h>
#include <math.h>
#include <gsl/gsl_math.h>

#define SDL_MAIN_HANDLED
#include <SDL2/SDL.h>

#include "libvecdisp.h"

int main(int argc, char * argv[]) {
	VECDISP_T ret;

	// Init
	ret = vecdisp_init();
	assert(ret == VECDISP_SUCCESS);
	ret = vecdisp_out_init();
	assert(ret == VECDISP_SUCCESS);

	SDL_GameController * controller_p0 = NULL, * controller_p1 = NULL;
	int controller_cnt = 0;

	SDL_Init( SDL_INIT_GAMECONTROLLER );
	SDL_GameControllerEventState(SDL_ENABLE);

	vecdisp_shape_t * shape_cube3d = vecdisp_shape_create(VECDISP_SHAPE_LINES, NULL, 0 );

	// format: {x, y, z}, x = horizontal, y = vertical, z = depth
	double cube3d[][3] = {
		{ -1, -1, -1 },
		{ 1, -1, -1 },
		{ 1, -1, -1 },
		{ 1, 1, -1 },
		{ 1, 1, -1 },
		{ -1, 1, -1 },
		{ -1, 1, -1 },
		{ -1, -1, -1 },
		{ -1, -1, -1 },
		{ -1, -1, 1 },
		{ -1, 1, -1 },
		{ -1, 1, 1 },
		{ 1, 1, -1 },
		{ 1, 1, 1 },
		{ 1, -1, -1 },
		{ 1, -1, 1 },
		{ -1, -1, 1 },
		{ 1, -1, 1 },
		{ 1, -1, 1 },
		{ 1, 1, 1 },
		{ 1, 1, 1 },
		{ -1, 1, 1 },
		{ -1, 1, 1 },
		{ -1, -1, 1 },
	 };
	for(int i = 0; i < 24; i++) {
		uint16_t data[2] = {0,0};
		vecdisp_shape_data_add(shape_cube3d, &data, 1 );

	}
	printf("%i\n", shape_cube3d->data_len);

	double rotx = 0, roty = 0, rotz = 0;
	double sqrt_3 = sqrt(3);
	double sqrt_6 = sqrt(6);

	// main loop
	SDL_Event event;
	bool quit = false;
	while( quit == false ) {
		while( SDL_PollEvent(&event) != 0 ) {
			switch (event.type) {
			case SDL_QUIT:
				quit = true;
				break;
			case SDL_CONTROLLERDEVICEADDED:
				if(controller_p0 == NULL) {
					controller_p0 = SDL_GameControllerOpen(event.cdevice.which);
					printf("adding Controller p0\n");
					controller_cnt++;
				}
				else if(controller_p1 == NULL) {
					controller_p1 = SDL_GameControllerOpen(event.cdevice.which);
					printf("adding Controller p1\n");
					controller_cnt++;
				}
				else {
					printf("Only two controllers allowed.\n");
				}
				break;
			case SDL_CONTROLLERDEVICEREMOVED:
				if (event.cdevice.which == SDL_JoystickInstanceID(SDL_GameControllerGetJoystick(controller_p0) ) )  {
					SDL_GameControllerClose(controller_p0);
					printf("freeing controller_p0\n");
					controller_p0 = NULL;
					controller_cnt--;
				}
				else if (event.cdevice.which == SDL_JoystickInstanceID(SDL_GameControllerGetJoystick(controller_p1) ) )  {
					SDL_GameControllerClose(controller_p1);
					printf("freeing controller_p1\n");
					controller_p1 = NULL;
					controller_cnt--;
				}
				break;
			case SDL_CONTROLLERBUTTONDOWN:
				switch (event.cbutton.button) {
				case SDL_CONTROLLER_BUTTON_BACK:
					if(event.cbutton.state == 1) {
						SDL_Event user_quit_event;
						user_quit_event.type = SDL_QUIT;
						SDL_PushEvent( &user_quit_event );
					}
					break;
				}
				break;
			case SDL_KEYDOWN:
				switch (event.key.keysym.sym) {
				case SDLK_q: {
					SDL_Event user_quit_event;
					user_quit_event.type = SDL_QUIT;
					SDL_PushEvent( &user_quit_event );
					break;
					}
				}
				break;
			}
		}

		rotx += 0.00004;
		roty += 0.00002;
		rotz += 0.00001;
		if( rotx >= M_PI ) rotx -= M_PI;
		if( roty >= M_PI ) roty -= M_PI;
		if( rotz >= M_PI ) rotz -= M_PI;

		for(int i = 0; i < 24; i++) {
			double c[3];
			cube3d[i][0] = cos(rotx) * cube3d[i][0] - sin(rotx) * cube3d[i][1]; // Rotation around X Axis
			cube3d[i][1] = sin(rotx) * cube3d[i][0] + cos(rotx) * cube3d[i][1];
			cube3d[i][2] = 1 * cube3d[i][2];

			cube3d[i][0] = cos(roty) * cube3d[i][0] - sin(roty) * cube3d[i][2]; // Rotation around Y Axis
			cube3d[i][1] = 1 * cube3d[i][1];
			cube3d[i][2] = cos(roty) * cube3d[i][2] + sin(roty) * cube3d[i][0];

			cube3d[i][0] = 1 * cube3d[i][0]; // Rotation around Z Axis
			cube3d[i][1] = cos(rotz) * cube3d[i][1] - sin(rotz) * cube3d[i][2]; // Rotation around Z Axis
			cube3d[i][2] = cos(rotz) * cube3d[i][2] + sin(rotz) * cube3d[i][1];

			c[0] = ( sqrt(3) * cube3d[i][0] + (-1) * sqrt_3 * cube3d[i][2] ) / sqrt_6;
			c[1] = ( 1 * cube3d[i][0] + 2 * cube3d[i][1] + 1 * cube3d[i][2] ) / sqrt_6;

			shape_cube3d->data[i][0] =  lround((DRAW_RES / 4) * c[0] + (DRAW_RES / 2));
			shape_cube3d->data[i][1] = lround((DRAW_RES / 4) * c[1] + (DRAW_RES / 2));
		}

		vecdisp_draw_shape( shape_cube3d, 0, 0, DRAW_RES - 1, DRAW_RES - 1, DRAW_BRTNS_BRIGHT);
		vecdisp_dbg_showfps();
		vecdisp_draw_update();
	}

	// cleaning up
	vecdisp_shape_destroy(shape_cube3d);
	SDL_QuitSubSystem( SDL_INIT_EVERYTHING );
	vecdisp_out_end();
	vecdisp_end();
	return 0;
}

ext: c

max_stars_repo_name: felixzwettler/vecdisp

max_issues_repo_name: fpunktz/vecdisp

max_forks_repo_name: fpunktz/vecdisp


text:
// Author:  Adam J. Robinson
// Summary: This file contains definition for the integral that
//          characterizes the interaction between nuclei and electrons,
//          as defined in the mathematical documentation.
// Performs the Electron Nuclei integral for the given
// Gaussian wavefunction.

#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_errno.h>

// Can't just include "gaussian.h", it will create a circular
// reference.
class GaussianWFN;

#ifndef INTEGRAL_EN
#define INTEGRAL_EN

class ENIntegrator {

public:

	gsl_monte_vegas_state * state;

	double * lowerBounds;
	double * upperBounds;
	GaussianWFN * wavefn;

	int    maxCalls;
	double width;
	double epsilon;

	// Initializes an integral calculator.
	//     fn     = The wavefunction that the computation is being carried
	//              out for.
	//     width  = The number of half width half maxes to
	//              limit the integral range to. Monte Carlo
	//              integration does not work for bounds at infinity.
	//              typically, 3 < width < 6
	//
	//     ncalls = The maximum number of calls to the integrand function
	//              that the integrator is permitted to make in order to
	//              calculate the integral. Higher values are slower and
	//              more accurate.
	ENIntegrator(GaussianWFN * fn, double width, int ncalls);

	~ENIntegrator();

	// These are set temporarily before an integral is performed.
	double A1;
	double A2;
	double A3;

	double s1;
	double s2;
	double s3;

	double R1;
	double R2;
	double R3;

	// Used to set the number of calls to the integrand that can be made.
	void setMaxCalls(int ncalls);

	// Given the three indices that define the integral being calculated,
	// sets the lowerBounds and upperBounds members of the class to
	// values appropriate for the given width, nuclei location, gaussian
	// center and gaussian width.
	void setBounds(int iu, int w, int l);

	double Integrate(int iu, int w, int l, double * error);

private:

	static double integrand(double * x, size_t dim, void * params);

};

#endif


text:
#pragma once

#include <memory>
#include <halley/maths/colour.h>
#include <halley/maths/vector2.h>
#include "halley/maths/rect.h"
#include "halley/data_structures/maybe.h"
#include <gsl/span>
#include <map>
#include "halley/core/graphics/sprite/sprite.h"
#include "halley/core/graphics/text/font.h"

namespace Halley
{
	class LocalisedString;
	class Painter;
	class Material;

	using ColourOverride = std::pair<size_t, std::optional<Colour4f>>;

	class TextRenderer
	{
	public:
		using SpriteFilter = std::function<void(gsl::span<Sprite>)>;

		TextRenderer();
		explicit TextRenderer(std::shared_ptr<const Font> font, String text = "", float size = 20, Colour colour = Colour(1, 1, 1, 1), float outline = 0, Colour outlineColour = Colour(0, 0, 0, 1));

		TextRenderer& setPosition(Vector2f pos);
		TextRenderer& setFont(std::shared_ptr<const Font> font);
		TextRenderer& setText(const String& text);
		TextRenderer& setText(const StringUTF32& text);
		TextRenderer& setText(const LocalisedString& text);
		TextRenderer& setSize(float size);
		TextRenderer& setColour(Colour colour);
		TextRenderer& setOutlineColour(Colour colour);
		TextRenderer& setOutline(float width);
		TextRenderer& setAlignment(float align);
		TextRenderer& setOffset(Vector2f align);
		TextRenderer& setClip(Rect4f clip);
		TextRenderer& setClip();
		TextRenderer& setSmoothness(float smoothness);
		TextRenderer& setPixelOffset(Vector2f offset);
		TextRenderer& setColourOverride(const std::vector<ColourOverride>& colOverride);
		TextRenderer& setLineSpacing(float spacing);

		TextRenderer clone() const;

		void generateSprites(std::vector<Sprite>& sprites) const;
		void draw(Painter& painter, const std::optional<Rect4f>& extClip = {}) const;

		void setSpriteFilter(SpriteFilter f);

		Vector2f getExtents() const;
		Vector2f getExtents(const StringUTF32& str) const;
		Vector2f getCharacterPosition(size_t character) const;
		Vector2f getCharacterPosition(size_t character, const StringUTF32& str) const;
		size_t getCharacterAt(const Vector2f& position) const;
		size_t getCharacterAt(const Vector2f& position, const StringUTF32& str) const;

		StringUTF32 split(const String& str, float width) const;
		StringUTF32 split(const StringUTF32& str, float width, std::function<bool(int32_t)> filter = {}) const;
		StringUTF32 split(float width) const;

		Vector2f getPosition() const;
		String getText() const;
		const StringUTF32& getTextUTF32() const;
		Colour getColour() const;
		float getOutline() const;
		Colour getOutlineColour() const;
		float getSmoothness() const;
		std::optional<Rect4f> getClip() const;
		float getLineHeight() const;
		float getAlignment() const;

		bool empty() const;

	private:
		std::shared_ptr<const Font> font;
		mutable std::map<const Font*, std::shared_ptr<Material>> materials;
		StringUTF32 text;
		SpriteFilter spriteFilter;

		float size = 20;
		float outline = 0;
		float align = 0;
		float smoothness = 1.0f;
		float lineSpacing = 0.0f;

		Vector2f position;
		Vector2f offset;
		Vector2f pixelOffset;
		Colour colour;
		Colour outlineColour;
		std::optional<Rect4f> clip;

		std::vector<ColourOverride> colourOverrides;

		mutable Vector<Sprite> spritesCache;
		mutable bool materialDirty = true;
		mutable bool glyphsDirty = true;
		mutable bool positionDirty = true;

		std::shared_ptr<Material> getMaterial(const Font& font) const;
		void updateMaterial(Material& material, const Font& font) const;
		void updateMaterialForFont(const Font& font) const;
		void updateMaterials() const;
		float getScale(const Font& font) const;
	};

	class ColourStringBuilder {
	public:
		void append(std::string_view text, std::optional<Colour4f> col = {});

		std::pair<String, std::vector<ColourOverride>> moveResults();

	private:
		std::vector<String> strings;
		std::vector<ColourOverride> colours;
		size_t len = 0;
	};
}

max_stars_repo_name: akien-mga/halley

max_issues_repo_name: akien-mga/halley

max_forks_repo_name: akien-mga/halley


text:
/* integration/qawf.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <math.h>
#include <float.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>

#include "initialise.c"
#include "append.c"
#include "qelg.c"

int
gsl_integration_qawf (gsl_function * f,
		      const double a,
		      const double epsabs,
		      const size_t limit,
		      gsl_integration_workspace * workspace,
		      gsl_integration_workspace * cycle_workspace,
		      gsl_integration_qawo_table * wf,
		      double *result, double *abserr)
{
  double area, errsum;
  double res_ext, err_ext;
  double correc, total_error = 0.0, truncation_error;

  size_t ktmin = 0;
  size_t iteration = 0;

  struct extrapolation_table table;

  double cycle;
  double omega = wf->omega;

  const double p = 0.9;
  double factor = 1;
  double initial_eps, eps;
  int error_type = 0;

  /* Initialize results */

  initialise (workspace, a, a);

  *result = 0;
  *abserr = 0;

  if (limit > workspace->limit)
    {
      GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
    }

  /* Test on accuracy */

  if (epsabs <= 0)
    {
      GSL_ERROR ("absolute tolerance epsabs must be positive", GSL_EBADTOL) ;
    }

  if (omega == 0.0)
    {
      if (wf->sine == GSL_INTEG_SINE)
	{
	  /* The function sin(w x) f(x) is always zero for w = 0 */

	  *result = 0;
	  *abserr = 0;

	  return GSL_SUCCESS;
	}
      else
	{
	  /* The function cos(w x) f(x) is always f(x) for w = 0 */

	  int status = gsl_integration_qagiu (f, a, epsabs, 0.0,
					      cycle_workspace->limit,
					      cycle_workspace,
					      result, abserr);
	  return status;
	}
    }

  if (epsabs > GSL_DBL_MIN / (1 - p))
    {
      eps = epsabs * (1 - p);
    }
  else
    {
      eps = epsabs;
    }

  initial_eps = eps;

  area = 0;
  errsum = 0;

  res_ext = 0;
  err_ext = GSL_DBL_MAX;
  correc = 0;

  cycle = (2 * floor (fabs (omega)) + 1) * M_PI / fabs (omega);

  gsl_integration_qawo_table_set_length (wf, cycle);

  initialise_table (&table);

  for (iteration = 0; iteration < limit; iteration++)
    {
      double area1, error1, reseps, erreps;

      double a1 = a + iteration * cycle;
      double b1 = a1 + cycle;

      double epsabs1 = eps * factor;

      int status = gsl_integration_qawo (f, a1, epsabs1, 0.0, limit,
					 cycle_workspace, wf,
					 &area1, &error1);

      append_interval (workspace, a1, b1, area1, error1);

      factor *= p;

      area = area + area1;
      errsum = errsum + error1;

      /* estimate the truncation error as 50 times the final term */

      truncation_error = 50 * fabs (area1);

      total_error = errsum + truncation_error;

      if (total_error < epsabs && iteration > 4)
	{
	  goto compute_result;
	}

      if (error1 > correc)
	{
	  correc = error1;
	}

      if (status)
	{
	  eps = GSL_MAX_DBL (initial_eps, correc * (1.0 - p));
	}

      if (status && total_error < 10 * correc && iteration > 3)
	{
	  goto compute_result;
	}

      append_table (&table, area);

      if (table.n < 2)
	{
	  continue;
	}

      qelg (&table, &reseps, &erreps);

      ktmin++;

      if (ktmin >= 15 && err_ext < 0.001 * total_error)
	{
	  error_type = 4;
	}

      if (erreps < err_ext)
	{
	  ktmin = 0;
	  err_ext = erreps;
	  res_ext = reseps;

	  if (err_ext + 10 * correc <= epsabs)
	    break;
	  if (err_ext <= epsabs && 10 * correc >= epsabs)
	    break;
	}

    }

  if (iteration == limit)
    error_type = 1;

  if (err_ext == GSL_DBL_MAX)
    goto compute_result;

  err_ext = err_ext + 10 * correc;

  *result = res_ext;
  *abserr = err_ext;

  if (error_type == 0)
    {
      return GSL_SUCCESS ;
    }

  if (res_ext != 0.0 && area != 0.0)
    {
      if (err_ext / fabs (res_ext) > errsum / fabs (area))
	goto compute_result;
    }
  else if (err_ext > errsum)
    {
      goto compute_result;
    }
  else if (area == 0.0)
    {
      goto return_error;
    }

  if (error_type == 4)
    {
      err_ext = err_ext + truncation_error;
    }

  goto return_error;

compute_result:

  *result = area;
  *abserr = total_error;

return_error:

  if (error_type > 2)
    error_type--;

  if (error_type == 0)
    {
      return GSL_SUCCESS;
    }
  else if (error_type == 1)
    {
      GSL_ERROR ("number of iterations was insufficient", GSL_EMAXITER);
    }
  else if (error_type == 2)
    {
      GSL_ERROR ("cannot reach tolerance because of roundoff error",
		 GSL_EROUND);
    }
  else if (error_type == 3)
    {
      GSL_ERROR ("bad integrand behavior found in the integration interval",
		 GSL_ESING);
    }
  else if (error_type == 4)
    {
      GSL_ERROR ("roundoff error detected in the extrapolation table",
		 GSL_EROUND);
    }
  else if (error_type == 5)
    {
      GSL_ERROR ("integral is divergent, or slowly convergent",
		 GSL_EDIVERGE);
    }
  else
    {
      GSL_ERROR ("could not integrate function", GSL_EFAILED);
    }

}


text:
/*
 *  Define digamma() and polygamma() here.
 *
 *  We use the GSL special functions, GSL is available
 *  on most platforms including MAC and Windows.
 *
 *  Another option is the Cephes library with psi() and zeta()
 *           http://www.netlib.org/cephes/
 *  which have nice, small self contained functions without all the
 *  GSL cruft.
 */

#ifndef __DIGAMMA_H
#define __DIGAMMA_H

/*
 *  Define this to switch off use of polygamma.
 *  Polygamma (trigamma, etc.) is used in some
 *  functions to speed things up.  So disabling
 *  it wont stop the functions working:
 *            gammadiff(), psidiff(), S_approx()
 *  But it also disables digammaInv() which
 *  uses a Newton-Raphson step.
 */
#define LS_NOPOLYGAMMA
/*
 *  Leaving this undefined means we use the library functions grabbed
 *  from the Mathlib and under GPL.  This makes the code self
 *  contained since no other libraries will be needed.
 *  Define this if you have GSL available to use instead.
 *  In which case, you will also need to modify the Makefile
 *  to include GSL.
 */
// #define GSL_POLYGAMMA

#ifdef LS_NOPOLYGAMMA
/*
 *   Radford Neal's implementation in digamma.c
 */
double digammaRN(double x);
#define digamma(x) digammaRN(x)

#else
double digammaInv(double x);

#ifndef GSL_POLYGAMMA
/*
 *   Mathlib implementation
 */
double MLpsigamma(double x, double deriv);
double MLdigamma(double x);
double MLtrigamma(double x);
double MLtetragamma(double x);
double MLpentagamma(double x);

#define digamma(x) MLdigamma(x)
#define trigamma(x) MLtrigamma(x)
#define tetragamma(x) MLtetragamma(x)
#define pentagamma(x) MLpentagamma(x)
#else
/*
 *   GSL library
 */
#include <gsl/gsl_sf.h>

#define digamma(x) gsl_sf_psi(x)
#define trigamma(x) gsl_sf_psi_n(1,x)
#define tetragamma(x) gsl_sf_psi_n(2,x)
#define pentagamma(x) gsl_sf_psi_n(3,x)

#endif
#endif

#endif

max_stars_repo_name: awbennett/HCA-WSI

max_issues_repo_name: awbennett/HCA-WSI

max_forks_repo_name: awbennett/HCA-WSI


text:
/*
 * resampling.c
 *
 *  Created on: 5.5.2017
 *      Author: heine
 */
#define _XOPEN_SOURCE 600

#include <math.h>
#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <mpi.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
//#include "ziggurat.h"

const gsl_rng_type *T; // Generator type
const gsl_rng *rgen;   // Generator

void set_resamplingseed( int id ) {

  struct timespec spec;

  clock_gettime( CLOCK_REALTIME, &spec );

  gsl_rng_env_setup();
  T = gsl_rng_default;
  rgen = gsl_rng_alloc( T );

  gsl_rng_set( rgen, ( unsigned long int ) ( round( spec.tv_nsec / 1.0e6 ) + 223 * id ) );

}

double depth_first_sum(double *x, long len) {

  long n = len;
  double *x_tmp, out;
  x_tmp = malloc(sizeof(double) * len);
  for (long i = 0; i < len; i++)
    x_tmp[i] = x[i];

  while (n > 1) {
    for (long i = 0; i < n / 2; i++) {
      x_tmp[i] = x_tmp[2 * i] + x_tmp[2 * i + 1];
    }
    if (n % 2 != 0) {
      x_tmp[n / 2] = x_tmp[n - 1];
    }
    n = n / 2 + n % 2;
  }
  out = x_tmp[0];
  free(x_tmp);
  return out;
}

double *depth_first_cumsum(double *x, long len) {

  double *W;
  W = malloc(sizeof(double) * len);

  for (long i = 0; i < len; i++) {
    W[i] = depth_first_sum(x, i);
  }
  return W;
}

double *exclusive_prefix_sum(double *w, long len) {
  double *W;
  W = malloc(sizeof(double) * len);

  W[0] = 0;
  for (long i = 1; i < len; i++) {
    W[i] = W[i - 1] + w[i - 1];
  }

  return W;
}

double *combining_exclusive_prefix_sum(double *w, double *w_recv, long n) {

  double *W;
  W = malloc(sizeof(double) * 2 * n);

  W[0] = 0;

  for (long i = 1; i <= n; i++) {
    W[i] = W[i - 1] + w[i - 1];
  }

  for (long i = 0; i < n-1; i++) {
    W[n + i + 1] = W[n + i] + w_recv[i];
  }

  return W;
}

double serial_multinomial_resample(long n, double *w, double *w_recv, long *inds) {

  double *W;
  double total_weight, lnmax = 0, u;
  long j = 2 * n - 1;

  W = combining_exclusive_prefix_sum(w, w_recv, n);

  total_weight = log(W[2 * n - 1] + w_recv[n - 1]);

  for (long i = n; i >= 1; i--) {

    u = gsl_rng_uniform_pos( rgen );

    lnmax = lnmax + log(u) / (double) i;

    u = exp(total_weight + lnmax);

    while (u < W[j])
      j--;

    inds[i - 1] = j;
  }

  free(W);

  return exp(total_weight);
}

double serial_multinomial_resample_sng(long n, double *w, long *inds) {

  double *W;
  double total_weight, lnmax = 0, u;
  long j = n - 1;

  W = exclusive_prefix_sum(w, n);

  total_weight = log(W[n - 1] + w[n - 1]);

  for (long i = n; i >= 1; i--) {

    u = gsl_rng_uniform_pos( rgen );

    lnmax = lnmax + log(u) / (double) i;

    u = exp(total_weight + lnmax);

    while (u < W[j])
      j--;

    inds[i - 1] = j;
  }

  free(W);

  return exp(total_weight);
}

double serial_multinomial_resample_sng_noutk(long n, double *w, long *inds, long k) {

  double *W;
  double total_weight, lnmax = 0, u;
  long j = k - 1;

  W = exclusive_prefix_sum(w, k);

  total_weight = log( W[ k - 1 ] + w[ k - 1 ] );

  for (long i = n; i >= 1; i--) {

    u = gsl_rng_uniform_pos( rgen );

    lnmax = lnmax + log(u) / (double) i;

    u = exp(total_weight + lnmax);

    while (u < W[j])
      j--;

    inds[i - 1] = j;
  }

  free(W);

  return exp(total_weight);
}

double serial_multinomial_resample_noutk(long n, double *w, double* W, long *inds, long k) {

  double total_weight, lnmax = 0, u;
  long j = k - 1;

  total_weight = log(W[k - 1] + w[k - 1]);

  for (long i = n; i >= 1; i--) {

    u = gsl_rng_uniform_pos( rgen );

    lnmax = lnmax + log(u) / (double) i;

    u = exp(total_weight + lnmax);

    while (u < W[j])
      j--;

    inds[i - 1] = j;
  }

  return -1;
}

void binary_resampler(long n, double w, long *counts)  {

  counts[0] = 0;
  counts[1] = 0;

  for ( long i = 0; i < n; i++ )
    counts[  gsl_rng_uniform_pos( rgen ) < w ? 0 : 1 ]++;

}

void
simplified_resample(
		    long n,
		    long *inds,
		    long k)
{

  for ( long i = 0; i < n; i++ )
    inds[ i ] = ( long ) ( gsl_rng_uniform_int( rgen, ( unsigned long int ) k ) );
}

void serial_multinomial_resample_counts(long n, long k, double *w, long *cnts) {

  double *W;
  double total_weight, lnmax = 0, u;
  long j = k - 1;

  for( long i = 0; i < k ; i++ )
    cnts[ i ] = 0;

  W = exclusive_prefix_sum(w, k);

  total_weight = log(W[k - 1] + w[k - 1]);

  for (long i = n; i >= 1; i--) {

    u = gsl_rng_uniform_pos( rgen );

    lnmax = lnmax + log(u) / (double) i;

    u = exp(total_weight + lnmax);

    while (u < W[j])
      j--;

    cnts[j]++;
  }

  free(W);

}

max_stars_repo_name: heinekmp/AIRPF

max_issues_repo_name: heinekmp/AIRPF

max_forks_repo_name: heinekmp/AIRPF


text:
/*
 *  vbHmmGaussDiffusion.c
 *  Model-specific core functions for VB-HMM-GAUSS-DIFFUSION.
 *
 *  Created by OKAMOTO Kenji, SAKO Yasushi and RIKEN
 *  Copyright 2011-2016
 *  Cellular Informatics Laboratory, Advance Science Institute, RIKEN, Japan.
 *  All rights reserved.
 *
 *  Ver. 1.0.0
 *  Last modified on 2016.02.08
 */

#include "vbHmmGaussDiffusion.h"
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_psi.h>
#include <string.h>
#include "rand.h"

#ifdef _OPENMP
#include "omp.h"
#endif

#define  MAX(a,b)  ((a)>(b)?(a):(b))
#define  MIN(a,b)  ((a)<(b)?(a):(b))

static int isGlobalAnalysis = 0;

void setFunctions_gaussDiff(){
    commonFunctions funcs;
    funcs.newModelParameters  = newModelParameters_gaussDiff;
    funcs.freeModelParameters = freeModelParameters_gaussDiff;
    funcs.newModelStats       = newModelStats_gaussDiff;
    funcs.freeModelStats      = freeModelStats_gaussDiff;
    funcs.initializeVbHmm     = initializeVbHmm_gaussDiff;
    funcs.pTilde_z1           = pTilde_z1_gaussDiff;
    funcs.pTilde_zn_zn1       = pTilde_zn_zn1_gaussDiff;
    funcs.pTilde_xn_zn        = pTilde_xn_zn_gaussDiff;
    funcs.calcStatsVars       = calcStatsVars_gaussDiff;
    funcs.maximization        = maximization_gaussDiff;
    funcs.varLowerBound       = varLowerBound_gaussDiff;
    funcs.reorderParameters   = reorderParameters_gaussDiff;
    funcs.outputResults       = outputResults_gaussDiff;
    setFunctions( funcs );
}

void setGFunctions_gaussDiff(){
    gCommonFunctions funcs;
    funcs.newModelParameters    = newModelParameters_gaussDiff;
    funcs.freeModelParameters   = freeModelParameters_gaussDiff;
    funcs.newModelStats         = newModelStats_gaussDiff;
    funcs.freeModelStats        = freeModelStats_gaussDiff;
    funcs.newModelStatsG        = newModelStatsG_gaussDiff;
    funcs.freeModelStatsG       = freeModelStatsG_gaussDiff;
    funcs.initializeVbHmmG      = initializeVbHmmG_gaussDiff;
    funcs.pTilde_z1             = pTilde_z1_gaussDiff;
    funcs.pTilde_zn_zn1         = pTilde_zn_zn1_gaussDiff;
    funcs.pTilde_xn_zn          = pTilde_xn_zn_gaussDiff;
    funcs.calcStatsVarsG        = calcStatsVarsG_gaussDiff;
    funcs.maximizationG         = maximizationG_gaussDiff;
    funcs.varLowerBoundG        = varLowerBoundG_gaussDiff;
    funcs.reorderParametersG    = reorderParametersG_gaussDiff;
    funcs.outputResultsG        = outputResultsG_gaussDiff;
    setGFunctions( funcs );
    isGlobalAnalysis = 1;
}

void outputResults_gaussDiff( xn, gv, iv, logFP )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
FILE *logFP;
{
    outputGaussDiffResults( xn, gv, iv, logFP );
}

void outputResultsG_gaussDiff( xns, gv, ivs, logFP )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
FILE *logFP;
{
    outputGaussDiffResultsG( xns, gv, ivs, logFP );
}

void *newModelParameters_gaussDiff( xn, sNo )
xnDataSet *xn;
int sNo;
{
    int i;
    gaussDiffParameters *p = (gaussDiffParameters*)malloc( sizeof(gaussDiffParameters) );

    p->uPiArr = (double*)malloc( sNo * sizeof(double) );
    p->sumUPi = 0.0;
    p->uAMat = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ ){
        p->uAMat[i] = (double*)malloc( sNo * sizeof(double) );
    }
    p->sumUAArr = (double*)malloc( sNo * sizeof(double) );

    p->avgPi = (double *)malloc( sNo * sizeof(double) );
    p->avgLnPi = (double *)malloc( sNo * sizeof(double) );
    p->avgA = (double **)malloc( sNo * sizeof(double*) );
    p->avgLnA = (double **)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ ){
        p->avgA[i] = (double *)malloc( sNo * sizeof(double) );
        p->avgLnA[i] = (double *)malloc( sNo * sizeof(double) );
    }
    p->avgDlt = (double *)malloc( sNo * sizeof(double) );
    p->avgLnDlt = (double *)malloc( sNo * sizeof(double) );

    p->uAArr = (double *)malloc( sNo * sizeof(double) );
    p->uBArr = (double *)malloc( sNo * sizeof(double) );
    p->aDlt = (double *)malloc( sNo * sizeof(double) );
    p->bDlt = (double *)malloc( sNo * sizeof(double) );

    return p;
}

void freeModelParameters_gaussDiff( p, xn, sNo )
void **p;
xnDataSet *xn;
int sNo;
{
    gaussDiffParameters *gp = *p;
    int i;

    free( gp->uPiArr );
    for( i = 0 ; i < sNo ; i++ ){
        free( gp->uAMat[i] );
    }
    free( gp->uAMat );
    free( gp->sumUAArr );

    free( gp->avgPi );
    free( gp->avgLnPi );
    for( i = 0 ; i < sNo ; i++ ){
        free( gp->avgA[i] );
        free( gp->avgLnA[i] );
    }
    free( gp->avgA );
    free( gp->avgLnA );
    free( gp->avgDlt );
    free( gp->avgLnDlt );

    free( gp->uAArr );
    free( gp->uBArr );
    free( gp->aDlt );
    free( gp->bDlt );

    free( *p );
    *p = NULL;
}

void *newModelStats_gaussDiff( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    if( isGlobalAnalysis == 0 ){
        int sNo = gv->sNo;
        gaussDiffStats *s = (gaussDiffStats*)malloc( sizeof(gaussDiffStats) );

        int i;
        s->Ni = (double *)malloc( sNo * sizeof(double) );
        s->Ri = (double *)malloc( sNo * sizeof(double) );
        s->Nij = (double **)malloc( sNo * sizeof(double*) );
        for( i = 0 ; i < sNo ; i++ )
        {   s->Nij[i] = (double *)malloc( sNo * sizeof(double) );   }
        s->Nii = (double *)malloc( sNo * sizeof(double) );

        return s;

    } else {

        return NULL;

    }
}

void freeModelStats_gaussDiff( s, xn, gv, iv )
void **s;
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    if( isGlobalAnalysis == 0 ){
        int sNo = gv->sNo;
        gaussDiffStats *gs = *s;
        int i;

        free( gs->Ni );
        free( gs->Ri );
        for( i = 0 ; i < sNo ; i++ )
        {   free( gs->Nij[i] );   }
        free( gs->Nij );
        free( gs->Nii );

        free( gs );
        *s = NULL;
    }
}

void *newModelStatsG_gaussDiff( xns, gv, ivs)
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    int sNo = gv->sNo;
    gaussDiffGlobalStats *gs = (gaussDiffGlobalStats*)malloc( sizeof(gaussDiffGlobalStats) );

    int i;
    gs->NiR = (double *)malloc( sNo * sizeof(double) );
    gs->RiR = (double *)malloc( sNo * sizeof(double) );
    gs->NijR = (double **)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ )
    {   gs->NijR[i] = (double *)malloc( sNo * sizeof(double) );   }
    gs->NiiR = (double *)malloc( sNo * sizeof(double) );
    gs->z1iR = (double *)malloc( sNo * sizeof(double) );

    return gs;
}

void freeModelStatsG_gaussDiff( gs, xns, gv, ivs )
void **gs;
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    int sNo = gv->sNo;
    gaussDiffGlobalStats *ggs = *gs;
    int i;
    free( ggs->NiR );
    for( i = 0 ; i < sNo ; i++ )
    {   free( ggs->NijR[i] );   }
    free( ggs->NijR );
    free( ggs->NiiR );
    free( ggs->RiR );
    free( ggs->z1iR );

    free( *gs );
    *gs = NULL;
}

void initializeVbHmm_gaussDiff( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    int sNo = gv->sNo;
    gaussDiffParameters *p = gv->params;
    int i, j;

    // hyper parameter for p( pi(i) )
    p->sumUPi = 0.0;
    for( i = 0 ; i < sNo ; i++ ){
        p->uPiArr[i] = 1.0;
        p->sumUPi += p->uPiArr[i];
    }

    // hyper parameter for p( A(i,j) )
    for( i = 0 ; i < sNo ; i++ ){
        p->sumUAArr[i] = 0.0;
        for( j = 0 ; j < sNo ; j++ ){
            if( j == i ){
                p->uAMat[i][j] = 5.0;
            } else {
                p->uAMat[i][j] = 1.0;
            }
            p->sumUAArr[i] += p->uAMat[i][j];
        }
    }

    // hyper parameter for p( delta(k) )
    for( i = 0 ; i < sNo ; i++ ){
        p->uAArr[i] = 1.0;
        p->uBArr[i] = 0.0005;
    }

    initialize_indVars_gaussDiff( xn, gv, iv );

    calcStatsVars_gaussDiff( xn, gv, iv );
    maximization_gaussDiff( xn, gv, iv );
}

void initializeVbHmmG_gaussDiff( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    int sNo = gv->sNo, rNo = xns->R;
    gaussDiffParameters *p = gv->params;
    int i, j, r;

    p->sumUPi = 0.0;
    for( i = 0 ; i < sNo ; i++ ){
        p->uPiArr[i] = 1.0;
        p->sumUPi += p->uPiArr[i];
    }

    for( i = 0 ; i < sNo ; i++ ){
        p->sumUAArr[i] = 0.0;
        for( j = 0 ; j < sNo ; j++ ){
            if( j == i ){
                p->uAMat[i][j] = 5.0;
            } else {
                p->uAMat[i][j] = 1.0;
            }
            p->sumUAArr[i] += p->uAMat[i][j];
        }
    }

    // hyper parameter for p( delta(k) )
    for( i = 0 ; i < sNo ; i++ ){
        p->uAArr[i] = 1.0;
        p->uBArr[i] = 0.0005;
    }

    for( r = 0 ; r < rNo ; r++ ){
        initialize_indVars_gaussDiff( xns->xn[r], gv, ivs->indVars[r] );
    }

    calcStatsVarsG_gaussDiff( xns, gv, ivs );
    maximizationG_gaussDiff( xns, gv, ivs );
}

void initialize_indVars_gaussDiff( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat;

    int i;
    size_t n;
    double sumPar;
    for( n = 0 ; n < dLen ; n++ ){
        sumPar = 0.0;
        for( i = 0 ; i < sNo ; i++ ){
            gmMat[n][i] = enoise(1.0) + 1.0;
            sumPar += gmMat[n][i];
        }
        for( i = 0 ; i < sNo ; i++ ){
            gmMat[n][i] /= sumPar;
        }
    }
}

xnDataSet *newXnDataSet_gaussDiff( filename )
const char *filename;
{
    xnDataSet *xn = (xnDataSet*)malloc( sizeof(xnDataSet) );
    xn->name = (char*)malloc( strlen(filename) + 2 );
    strncpy( xn->name, filename, strlen(filename)+1 );
    xn->data = (gaussDiffData*)malloc( sizeof(gaussDiffData) );
    gaussDiffData *d = (gaussDiffData*)xn->data;
    d->v = NULL;
    return xn;
}

void freeXnDataSet_gaussDiff( xn )
xnDataSet **xn;
{
    gaussDiffData *d = (gaussDiffData*)(*xn)->data;
    free( d->v );
    free( (*xn)->data );
    free( (*xn)->name );
    free( *xn );
    *xn = NULL;
}

double pTilde_z1_gaussDiff( i, params )
int i;
void *params;
{
    gaussDiffParameters *p = (gaussDiffParameters*)params;
    return exp( p->avgLnPi[i] );
}

double pTilde_zn_zn1_gaussDiff( i, j, params )
int i, j;
void *params;
{
    gaussDiffParameters *p = (gaussDiffParameters*)params;
    return exp( p->avgLnA[i][j] );
}

double pTilde_xn_zn_gaussDiff( xnWv, n, i, params )
xnDataSet *xnWv;
size_t n;
int i;
void *params;
{
    gaussDiffParameters *p = (gaussDiffParameters*)params;
    gaussDiffData *xn = (gaussDiffData*)xnWv->data;
    double val;
    val  = p->avgLnDlt[i] - log(2.0);
    val -= log(xn->v[n]) + p->avgDlt[i] * pow( xn->v[n], 2.0) / 4.0;
    return exp(val);
}

void calcStatsVars_gaussDiff( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    gaussDiffData *d = (gaussDiffData*)xn->data;
    gaussDiffStats *s = (gaussDiffStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, ***xiMat = iv->xiMat;
    double *Ni = s->Ni, *Ri = s->Ri, *Nii = s->Nii, **Nij = s->Nij;
    size_t n;
    int i, j;

    for( i = 0 ; i < sNo ; i++ ){
        Ni[i]   = 1e-10;
        Ri[i]   = 1e-10;
        Nii[i]  = 1e-10;
        for( j = 0 ; j < sNo ; j++ ){
            Nij[i][j] = 1e-10;
        }
        for( n = 0 ; n < dLen ; n++ ){
            Ni[i] += gmMat[n][i];
            Ri[i] += gmMat[n][i] * pow( d->v[n], 2.0 );
            for( j = 0 ; j < sNo ; j++ ){
                Nii[i]    += xiMat[n][i][j];
                Nij[i][j] += xiMat[n][i][j];
            }
        }
    }
}

void calcStatsVarsG_gaussDiff( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    gaussDiffData *d;
    gaussDiffGlobalStats *gs = (gaussDiffGlobalStats*)ivs->stats;
    int sNo = gv->sNo, rNo = xns->R;
    double **gmMat, ***xiMat;
    double *NiR = gs->NiR, *RiR = gs->RiR, *NiiR = gs->NiiR;
    double **NijR = gs->NijR, *z1iR = gs->z1iR;
    size_t dLen, n;
    int i, j, r;

    for( i = 0 ; i < sNo ; i++ ){
        NiR[i]   = 1e-10;
        RiR[i]  = 1e-10;
        NiiR[i]  = 1e-10;
        for( j = 0 ; j < sNo ; j++ ){
            NijR[i][j] = 1e-10;
        }
        z1iR[i] = 1e-10;
    }
    for( r = 0 ; r < rNo ; r++ ){
        d = (gaussDiffData*)xns->xn[r]->data;
        dLen = xns->xn[r]->N;
        gmMat = ivs->indVars[r]->gmMat;
        xiMat = ivs->indVars[r]->xiMat;
        for( i = 0 ; i < sNo ; i++ ){
            z1iR[i] += gmMat[0][i];
            for( n = 0 ; n < dLen ; n++ ){
                NiR[i] += gmMat[n][i];
                RiR[i] += gmMat[n][i] * pow( d->v[n], 2.0 );
                for( j = 0 ; j < sNo ; j++ ){
                    NiiR[i]    += xiMat[n][i][j];
                    NijR[i][j] += xiMat[n][i][j];
                }
            }
        }
    }
}

void maximization_gaussDiff( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    gaussDiffParameters *p = (gaussDiffParameters*)gv->params;
    gaussDiffStats *s = (gaussDiffStats*)iv->stats;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double *uAArr = p->uAArr, *uBArr = p->uBArr, *aDlt = p->aDlt, *bDlt = p->bDlt;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    double *avgDlt = p->avgDlt, *avgLnDlt = p->avgLnDlt;
    double *Ni = s->Ni, *Ri = s->Ri, *Nii = s->Nii, **Nij = s->Nij;
    int i, j;

    for( i = 0 ; i < sNo ; i++ ){
        avgPi[i] = ( uPiArr[i] + gmMat[0][i] ) / ( sumUPi + 1.0 );
        avgLnPi[i] = gsl_sf_psi( uPiArr[i] + gmMat[0][i] ) - gsl_sf_psi( sumUPi + 1.0 );

        for( j = 0 ; j < sNo ; j++ ){
            avgA[i][j] = ( uAMat[i][j] + Nij[i][j] ) / ( sumUAArr[i] + Nii[i] );
            avgLnA[i][j] = gsl_sf_psi( uAMat[i][j] + Nij[i][j] ) - gsl_sf_psi( sumUAArr[i] + Nii[i] );
        }

        aDlt[i] = uAArr[i] + Ni[i];
        bDlt[i] = uBArr[i] + Ri[i] / 4.0;

        avgDlt[i]   = aDlt[i] / bDlt[i];
        avgLnDlt[i] = gsl_sf_psi( aDlt[i] ) - log( bDlt[i] );
    }
}

void maximizationG_gaussDiff( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    gaussDiffParameters *p = (gaussDiffParameters*)gv->params;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double *uAArr = p->uAArr, *uBArr = p->uBArr, *aDlt = p->aDlt, *bDlt = p->bDlt;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    double *avgDlt = p->avgDlt, *avgLnDlt = p->avgLnDlt;
    gaussDiffGlobalStats *gs = (gaussDiffGlobalStats*)ivs->stats;
    double *NiR = gs->NiR, *NiiR = gs->NiiR, **NijR = gs->NijR, *RiR = gs->RiR;
    double *z1iR = gs->z1iR, dR = (double)(xns->R);
    int sNo = gv->sNo;
    int i, j;

    for( i = 0 ; i < sNo ; i++ ){
        avgPi[i] = ( uPiArr[i] + z1iR[i] ) / ( sumUPi + dR );
        avgLnPi[i] = gsl_sf_psi( uPiArr[i] + z1iR[i] ) - gsl_sf_psi( sumUPi + dR );

        for( j = 0 ; j < sNo ; j++ ){
            avgA[i][j] = ( uAMat[i][j] + NijR[i][j] ) / ( sumUAArr[i] + NiiR[i] );
            avgLnA[i][j] = gsl_sf_psi( uAMat[i][j] + NijR[i][j] ) - gsl_sf_psi( sumUAArr[i] + NiiR[i] );
        }

        aDlt[i] = uAArr[i] + NiR[i];
        bDlt[i] = uBArr[i] + RiR[i] / 4.0;

        avgDlt[i]   = aDlt[i] / bDlt[i];
        avgLnDlt[i] = gsl_sf_psi( aDlt[i] ) - log( bDlt[i] );
    }
}

double varLowerBound_gaussDiff( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    gaussDiffParameters *p = (gaussDiffParameters*)gv->params;
    gaussDiffStats *s = (gaussDiffStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, *cn = iv->cn;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double *uAArr = p->uAArr, *uBArr = p->uBArr;
    double *avgLnPi = p->avgLnPi, **avgLnA = p->avgLnA;
    double *avgDlt = p->avgDlt, *avgLnDlt = p->avgLnDlt;
    double *Nii = s->Nii, **Nij = s->Nij;
    double *aDlt = p->aDlt, *bDlt = p->bDlt;
    size_t n;
    int i, j;

    double lnpPi = gsl_sf_lngamma(sumUPi);
    double lnpA = 0.0;
    double lnpDlt = 0.0;
    double lnqPi = gsl_sf_lngamma(sumUPi + 1.0);
    double lnqA = 0.0;
    double lnqDlt = - sNo / 2.0;
    for( i = 0 ; i < sNo ; i++ ){
        lnpPi += (uPiArr[i]-1.0) * avgLnPi[i] - gsl_sf_lngamma(uPiArr[i]);

        lnpDlt += - gsl_sf_lngamma(uAArr[i]) + uAArr[i] * log(uBArr[i]);
        lnpDlt += (uAArr[i] - 1.0) * avgLnDlt[i] - uBArr[i] * avgDlt[i];

        lnqPi += (uPiArr[i]+gmMat[0][i]-1.0) * (gsl_sf_psi(uPiArr[i]+gmMat[0][i]) - gsl_sf_psi(sumUPi+1.0));
        lnqPi -= gsl_sf_lngamma(uPiArr[i] + gmMat[0][i]);

        lnpA += gsl_sf_lngamma(sumUAArr[i]);
        lnqA += gsl_sf_lngamma(sumUAArr[i] + Nii[i]);
        for( j = 0 ; j < sNo ; j++ ){
            lnpA += (uAMat[i][j]-1.0)*avgLnA[i][j] - gsl_sf_lngamma(uAMat[i][j]);

            lnqA += (uAMat[i][j] + Nij[i][j] - 1.0) * (gsl_sf_psi(uAMat[i][j]+Nij[i][j]) - gsl_sf_psi(sumUAArr[i]+Nii[i]));
            lnqA -= gsl_sf_lngamma( uAMat[i][j] + Nij[i][j] );
        }

        lnqDlt += - gsl_sf_lngamma(aDlt[i]) + aDlt[i] * log(bDlt[i]);
        lnqDlt += (aDlt[i] - 1.0) * avgLnDlt[i] - aDlt[i];
    }

    double lnpX = 0.0;
    for( n = 0 ; n < dLen ; n++ ){
        lnpX += log( cn[n] );
    }

    double val;
    val  = lnpPi + lnpA + lnpDlt;
    val -= lnqPi + lnqA + lnqDlt;
    val += lnpX;
    val += log(gsl_sf_fact(sNo));

    return val;
}

double varLowerBoundG_gaussDiff( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    gaussDiffParameters *p = (gaussDiffParameters*)gv->params;
    int sNo = gv->sNo, rNo = xns->R;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi;
    double **uAMat = p->uAMat, *sumUAArr = p->sumUAArr;
    double *uAArr = p->uAArr, *uBArr = p->uBArr;
    double *avgLnPi = p->avgLnPi, **avgLnA = p->avgLnA;
    double *avgDlt = p->avgDlt, *avgLnDlt = p->avgLnDlt;
    double *aDlt = p->aDlt, *bDlt = p->bDlt;
    gaussDiffGlobalStats *gs = (gaussDiffGlobalStats*)ivs->stats;
    double *NiiR = gs->NiiR, **NijR = gs->NijR, *z1iR = gs->z1iR, dR = (double)xns->R;
    size_t n;
    int i, j, r;

    double lnpPi = gsl_sf_lngamma(sumUPi);
    double lnpA = 0.0;
    double lnpDlt = 0.0;
    double lnqPi = gsl_sf_lngamma(sumUPi + dR);
    double lnqA = 0.0;
    double lnqDlt = - sNo / 2.0;
    for( i = 0 ; i < sNo ; i++ ){
        lnpPi += (uPiArr[i]-1.0) * avgLnPi[i] - gsl_sf_lngamma(uPiArr[i]);

        lnpDlt += - gsl_sf_lngamma(uAArr[i]) + uAArr[i] * log(uBArr[i]);
        lnpDlt += (uAArr[i] - 1.0) * avgLnDlt[i] - uBArr[i] * avgDlt[i];

        lnqPi += (uPiArr[i]+z1iR[i]-1.0) * (gsl_sf_psi(uPiArr[i]+z1iR[i]) - gsl_sf_psi(sumUPi+dR));
        lnqPi -= gsl_sf_lngamma(uPiArr[i] + z1iR[i]);

        lnpA += gsl_sf_lngamma(sumUAArr[i]);
        lnqA += gsl_sf_lngamma(sumUAArr[i] + NiiR[i]);
        for( j = 0 ; j < sNo ; j++ ){
            lnpA += (uAMat[i][j]-1.0)*avgLnA[i][j] - gsl_sf_lngamma(uAMat[i][j]);

            lnqA += (uAMat[i][j] + NijR[i][j] - 1.0) * (gsl_sf_psi(uAMat[i][j]+NijR[i][j]) - gsl_sf_psi(sumUAArr[i]+NiiR[i]));
            lnqA -= gsl_sf_lngamma( uAMat[i][j] + NijR[i][j] );
        }

        lnqDlt += - gsl_sf_lngamma(aDlt[i]) + aDlt[i] * log(bDlt[i]);
        lnqDlt += (aDlt[i] - 1.0) * avgLnDlt[i] - aDlt[i];
    }

    double lnpX = 0.0;
    for( r = 0 ; r < rNo ; r++ ){
        size_t dLen = xns->xn[r]->N;
        for( n = 0 ; n < dLen ; n++ ){
            lnpX += log( ivs->indVars[r]->cn[n] );
        }
    }

    double val;
    val  = lnpPi + lnpA + lnpDlt;
    val -= lnqPi + lnqA + lnqDlt;
    val += lnpX;
    val += log(gsl_sf_fact(sNo));

    return val;
}

void reorderParameters_gaussDiff( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    gaussDiffParameters *p = (gaussDiffParameters*)gv->params;
    gaussDiffStats *s = (gaussDiffStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, ***xiMat = iv->xiMat;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    double *avgDlt = p->avgDlt, *avgLnDlt = p->avgLnDlt;
    double *Ni = s->Ni;
    size_t n;
    int i, j;

    int *index = (int*)malloc( sNo * sizeof(int) );
    double *store = (double*)malloc( sNo * sizeof(double) );
    double **s2D = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ )
    {   s2D[i] = (double*)malloc( MAX(sNo,2) * sizeof(double) );   }

    // index indicates order of avgDlt values (0=biggest avgDlt -- sNo=smallest avgDlt).
    for( i = 0 ; i < sNo ; i++ ){
        index[i] = sNo - 1;
        for( j = 0 ; j < sNo ; j++ ){
            if( j != i ){
                if( avgDlt[i] < avgDlt[j] ){
                    index[i]--;
                } else if( avgDlt[i] == avgDlt[j] ){
                    if( j > i )
                    {   index[i]--;   }
                }
            }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgDlt[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgDlt[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnDlt[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnDlt[i] = store[i];   }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgA[i][j] = s2D[i][j];   }
    }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgLnA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgLnA[i][j] = s2D[i][j];   }
    }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = Ni[i];   }
    for( i = 0 ; i < sNo ; i++ ){   Ni[i] = store[i];   }

    for( n = 0 ; n < dLen ; n++ ){
        for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = gmMat[n][i];   }
        for( i = 0 ; i < sNo ; i++ ){   gmMat[n][i] = store[i];   }
    }

    for( n = 0 ; n < dLen ; n++ ){
        for( j = 0 ; j < sNo ; j++ ){
            for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = xiMat[n][i][j];   }
        }
        for( j = 0 ; j < sNo ; j++ ){
            for( i = 0 ; i < sNo ; i++ ){   xiMat[n][i][j] = s2D[i][j];   }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   free( s2D[i] );   }
    free( s2D );
    free( store );
    free( index );
}

void reorderParametersG_gaussDiff( xns, gv, ivs )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
{
    gaussDiffParameters *p = (gaussDiffParameters*)gv->params;
    size_t dLen;
    int sNo = gv->sNo, rNo = xns->R;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgA = p->avgA, **avgLnA = p->avgLnA;
    double *avgDlt = p->avgDlt, *avgLnDlt = p->avgLnDlt;
    size_t n;
    int i, j, r;

    int *index = (int*)malloc( sNo * sizeof(int) );
    double *store = (double*)malloc( sNo * sizeof(double) );
    double **s2D = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ )
    {   s2D[i] = (double*)malloc( MAX(sNo,2) * sizeof(double) );   }

    // index indicates order of avgMu values (0=biggest avgMu -- sNo=smallest avgMu).
    for( i = 0 ; i < sNo ; i++ ){
        index[i] = sNo - 1;
        for( j = 0 ; j < sNo ; j++ ){
            if( j != i ){
                if( avgDlt[i] < avgDlt[j] ){
                    index[i]--;
                } else if( avgDlt[i] == avgDlt[j] ){
                    if( j > i )
                    {   index[i]--;   }
                }
            }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgDlt[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgDlt[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnDlt[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnDlt[i] = store[i];   }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgA[i][j] = s2D[i][j];   }
    }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgLnA[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgLnA[i][j] = s2D[i][j];   }
    }

    double *NiR = ((gaussDiffGlobalStats*)ivs->stats)->NiR;
    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = NiR[i];   }
    for( i = 0 ; i < sNo ; i++ ){   NiR[i] = store[i];   }

    for( r = 0 ; r < rNo ; r++ ){
        double **gmMat = ivs->indVars[r]->gmMat;
        dLen = xns->xn[r]->N;
        for( n = 0 ; n < dLen ; n++ ){
            for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = gmMat[n][i];   }
            for( i = 0 ; i < sNo ; i++ ){   gmMat[n][i] = store[i];   }
        }
    }

    for( r = 0 ; r < rNo ; r++ ){
        double ***xiMat = ivs->indVars[r]->xiMat;
        dLen = xns->xn[r]->N;
        for( n = 0 ; n < dLen ; n++ ){
            for( j = 0 ; j < sNo ; j++ ){
                for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = xiMat[n][i][j];   }
            }
            for( j = 0 ; j < sNo ; j++ ){
                for( i = 0 ; i < sNo ; i++ ){   xiMat[n][i][j] = s2D[i][j];   }
            }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   free( s2D[i] );   }
    free( s2D );
    free( store );
    free( index );
}

void outputGaussDiffResults( xn, gv, iv, logFP )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
FILE *logFP;
{
    gaussDiffParameters *p = (gaussDiffParameters*)gv->params;
    int sNo = gv->sNo;
    int i, j;
    fprintf(logFP, "  results: K = %d \n", sNo);

    fprintf(logFP, "   delta: ( %g", p->avgDlt[0]);
    for( i = 1 ; i < sNo ; i++ ){
        fprintf(logFP, ", %g", p->avgDlt[i]);
    }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   pi: ( %g", p->avgPi[0]);
    for( i = 1 ; i < sNo ; i++ ){
        fprintf(logFP, ", %g", p->avgPi[i]);
    }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   A_matrix: [");
    for( i = 0 ; i < sNo ; i++ ){
            fprintf(logFP, " ( %g", p->avgA[i][0]);
            for( j = 1 ; j < sNo ; j++ )
            {   fprintf(logFP, ", %g", p->avgA[i][j]);   }
            fprintf(logFP, ")");
    }
    fprintf(logFP, " ] \n\n");

    char fn[256];
    FILE *fp;
    size_t n;

    sprintf( fn, "%s.param%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        fprintf(fp, "delta, pi");
        for( i = 0 ; i < sNo ; i++ )
        {   fprintf(fp, ", A%dx", i);   }
        fprintf(fp, "\n");

        for( i = 0 ; i < sNo ; i++ ){
            fprintf(fp, "%g, %g", p->avgDlt[i], p->avgPi[i]);
            for( j = 0 ; j < sNo ; j++ )
            {   fprintf(fp, ", %g", p->avgA[j][i]);   }
            fprintf(fp, "\n");
        }
        fclose(fp);
    }

    sprintf( fn, "%s.Lq%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < gv->iteration ; n++ ){
            fprintf( fp, "%24.20e\n", gv->LqArr[n] );
        }
        fclose(fp);
    }

    sprintf( fn, "%s.maxS%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < xn->N ; n++ ){
            fprintf( fp, "%d\n", iv->stateTraj[n] );
        }
        fclose(fp);
    }

}

void outputGaussDiffResultsG( xns, gv, ivs, logFP )
xnDataBundle *xns;
globalVars *gv;
indVarBundle *ivs;
FILE *logFP;
{
    int sNo = gv->sNo, rNo = xns->R;
    gaussDiffParameters *p = (gaussDiffParameters*)gv->params;
    int i, j, r;

    fprintf(logFP, "  results: K = %d \n", sNo);

    fprintf(logFP, "   delta: ( %g", p->avgDlt[0]);
    for( i = 1 ; i < sNo ; i++ )
    {   fprintf(logFP, ", %g", p->avgDlt[i]);   }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   pi: ( %g", p->avgPi[0]);
    for( i = 1 ; i < sNo ; i++ ){
        fprintf(logFP, ", %g", p->avgPi[i]);
    }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   A_matrix: [");
    for( i = 0 ; i < sNo ; i++ ){
        fprintf(logFP, " ( %g", p->avgA[i][0]);
        for( j = 1 ; j < sNo ; j++ )
        {   fprintf(logFP, ", %g", p->avgA[i][j]);   }
        fprintf(logFP, ")");
    }
    fprintf(logFP, " ] \n\n");

    char fn[256];
    FILE *fp;
    size_t n;

    sprintf( fn, "%s.param%03d", xns->xn[0]->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        fprintf(fp, "delta, pi");
        for( i = 0 ; i < sNo ; i++ )
        {   fprintf(fp, ", A%dx", i);   }
        fprintf(fp, "\n");

        for( i = 0 ; i < sNo ; i++ ){
            fprintf(fp, "%g, %g", p->avgDlt[i], p->avgPi[i]);
            for( j = 0 ; j < sNo ; j++ )
            {   fprintf(fp, ", %g", p->avgA[j][i]);   }
            fprintf(fp, "\n");
        }
        fclose(fp);
    }

    sprintf( fn, "%s.Lq%03d", xns->xn[0]->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < gv->iteration ; n++ ){
            fprintf( fp, "%24.20e\n", gv->LqArr[n] );
        }
        fclose(fp);
    }

    sprintf( fn, "%s.maxS%03d", xns->xn[0]->name, sNo );
    int flag = 0;
    if( (fp = fopen( fn, "w")) != NULL ){
        n = 0;
        do{
            flag = 1;
            for( r = 0 ; r < rNo ; r++ ){
                xnDataSet *xn = xns->xn[r];
                indVars *iv = ivs->indVars[r];

                if( r > 0 ){
                    fprintf( fp, "," );
                }
                if( n < xn->N ){
                    fprintf( fp, "%d", iv->stateTraj[n] );
                }
                flag &= (n >= (xn->N - 1));
            }
            fprintf( fp, "\n" );
            n++;
        }while( !flag );
        fclose(fp);
    }
}

//

max_stars_repo_name: okamoto-kenji/varBayes-HMM

max_issues_repo_name: okamoto-kenji/varBayes-HMM

max_forks_repo_name: okamoto-kenji/varBayes-HMM


text:
/* randist/demo.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int
main ()
{
  gsl_rng * r ;

  int i, n = 10;
  double mu = 3.0;

  /* create a generator chosen by the environment variable GSL_RNG_TYPE */

  gsl_rng_env_setup();

  r = gsl_rng_alloc (gsl_rng_default);

  /* print n random variates chosen from the poisson distribution with
     mean parameter mu */

  for (i = 0; i < n; i++)
    {
      unsigned int k = gsl_ran_poisson (r, mu);

      printf(" %u", k);
    }

  printf("\n");
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#pragma once

#define OEMRESOURCE

#include <atomic>
#include <concepts>

#include <ppl.h>

#include <gsl/gsl>

// Included before windows.h, because pfc.h includes winsock2.h
#include "../pfc/pfc.h"

#include <windows.h>
#include <SHLWAPI.H>

#include "../foobar2000/SDK/foobar2000.h"

#include "../mmh/stdafx.h"
#include "../ui_helpers/stdafx.h"

#include "config_var.h"
#include "config_object.h"
#include "console.h"
#include "fcl.h"
#include "info_box.h"
#include "initquit.h"
#include "library.h"
#include "low_level_hook.h"
#include "main_thread_callback.h"
#include "sort.h"
#include "stream.h"

max_stars_repo_name: reupen/fbh

max_issues_repo_name: reupen/fbh

max_forks_repo_name: reupen/fbh


text:
/**
 *
 * @file core_zsyssq.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.8.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <math.h>
#include <lapacke.h>
#include "common.h"

#define COMPLEX

#define UPDATE( __nb, __value )                                         \
    if (__value != 0. ){                                                \
        if ( *scale < __value ) {                                       \
            *sumsq = __nb + (*sumsq) * ( *scale / __value ) * ( *scale / __value ); \
            *scale = __value;                                           \
        } else {                                                        \
            *sumsq = *sumsq + __nb * ( __value / *scale ) *  ( __value / *scale ); \
        }                                                               \
    }

/*****************************************************************************
 *
 * @ingroup dplasma_cores_complex64
 *
 *  CORE_zsyssq returns the values scl and ssq such that
 *
 *    ( scl**2 )*ssq = sum( A( i, j )**2 ) + ( scale**2 )*sumsq,
 *                     i,j
 *
 * with i from 0 to N-1 and j form 0 to N-1. The value of sumsq is
 * assumed to be at least unity and the value of ssq will then satisfy
 *
 *    1.0 .le. ssq .le. ( sumsq + 2*n*n ).
 *
 * scale is assumed to be non-negative and scl returns the value
 *
 *    scl = max( scale, abs( real( A( i, j ) ) ), abs( aimag( A( i, j ) ) ) ),
 *          i,j
 *
 * scale and sumsq must be supplied in SCALE and SUMSQ respectively.
 * SCALE and SUMSQ are overwritten by scl and ssq respectively.
 *
 * The routine makes only one pass through the tile triangular part of the
 * symmetric tile A defined by uplo.
 * See also LAPACK zlassq.f
 *
 *******************************************************************************
 *
 *  @param[in] uplo
 *          Specifies whether the upper or lower triangular part of
 *          the symmetric matrix A is to be referenced as follows:
 *          = PlasmaLower:     Only the lower triangular part of the
 *                             symmetric matrix A is to be referenced.
 *          = PlasmaUpper:     Only the upper triangular part of the
 *                             symmetric matrix A is to be referenced.
 *
 *  @param[in] N
 *          The number of columns and rows in the tile A.
 *
 *  @param[in] A
 *          The N-by-N matrix on which to compute the norm.
 *
 *  @param[in] LDA
 *          The leading dimension of the tile A. LDA >= max(1,N).
 *
 *  @param[in,out] scale
 *          On entry, the value  scale  in the equation above.
 *          On exit, scale is overwritten with the value scl.
 *
 *  @param[in,out] sumsq
 *          On entry, the value  sumsq  in the equation above.
 *          On exit, SUMSQ is overwritten with the value ssq.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval -k, the k-th argument had an illegal value
 *
 */
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zsyssq = PCORE_zsyssq
#define CORE_zsyssq PCORE_zsyssq
#endif
int CORE_zsyssq(PLASMA_enum uplo, int N,
                const PLASMA_Complex64_t *A, int LDA,
                double *scale, double *sumsq)
{
    int i, j;
    double tmp;
    double *ptr;

    if ( uplo == PlasmaUpper ) {
        for(j=0; j<N; j++) {
            ptr = (double*) ( A + j * LDA );

            for(i=0; i<j; i++, ptr++) {

                tmp = fabs(*ptr);
                UPDATE( 2., tmp );

#ifdef COMPLEX
                ptr++;
                tmp = fabs(*ptr);
                UPDATE( 2., tmp );
#endif
            }

            /* Diagonal */
            tmp = fabs(*ptr);
            UPDATE( 1., tmp );

#ifdef COMPLEX
            ptr++;
            tmp = fabs(*ptr);
            UPDATE( 1., tmp );
#endif
        }
    } else {

        for(j=0; j<N; j++) {
            ptr = (double*) ( A + j * LDA + j);

            /* Diagonal */
            tmp = fabs(*ptr);
            UPDATE( 1., tmp );
            ptr++;

#ifdef COMPLEX
            tmp = fabs(*ptr);
            UPDATE( 1., tmp );
            ptr++;
#endif

            for(i=j+1; i<N; i++, ptr++) {

                tmp = fabs(*ptr);
                UPDATE( 2., tmp );

#ifdef COMPLEX
                ptr++;
                tmp = fabs(*ptr);
                UPDATE( 2., tmp );
#endif
            }
        }
    }
    return PLASMA_SUCCESS;
}

max_stars_repo_name: therault/dplasma

max_issues_repo_name: therault/dplasma

max_forks_repo_name: therault/dplasma


text:
/* -*- linux-c -*- */
/* fewbody_nonks.c

   Copyright (C) 2002-2004 John M. Fregeau

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include "fewbody.h"

/* the derivatives function for the GSL ODE integrator */
#define FB_FM(i, j, k) fm[nstar*3*i + 3*j + k]
int fb_nonks_func(double t, const double *y, double *f, void *params)
{
	int i, j, k, nstar;
	double r[3], *m, *fm, val;

	nstar = (*(fb_nonks_params_t *) params).nstar;
	m = (*(fb_nonks_params_t *) params).m;

	fm = fb_malloc_vector(nstar * nstar * 3);

	/* calculate the matrix */
	for (i=0; i<nstar; i++) {
		for (j=0; j<i; j++) {
			for (k=0; k<3; k++) {
				FB_FM(i, j, k) = -FB_FM(j, i, k);
			}
		}

		for (j=i+1; j<nstar; j++) {
			for (k=0; k<3; k++) {
				r[k] = y[j*6+k] - y[i*6+k];
			}

			val = 1.0 / fb_cub(fb_mod(r));

			for (k=0; k<3; k++) {
				FB_FM(i, j, k) = val * r[k];
			}
		}
	}

	/* calculate derivatives */
	for (i=0; i<nstar; i++) {
		for (k=0; k<3; k++) {
			f[i*6+k] = y[i*6+k+3];
			f[i*6+k+3] = 0.0;
		}

		for (j=0; j<i; j++) {
			for (k=0; k<3; k++) {
				f[i*6+k+3] += m[j] * FB_FM(i, j, k);
			}
		}

		for (j=i+1; j<nstar; j++) {
			for (k=0; k<3; k++) {
				f[i*6+k+3] += m[j] * FB_FM(i, j, k);
			}
		}
	}

	fb_free_vector(fm);

	return(GSL_SUCCESS);
}
#undef FB_FM

/* the Jacobian for the GSL ODE integrator */
int fb_nonks_jac(double t, const double *y, double *dfdy, double *dfdt, void *params)
{
	unsigned int i, j, k, a, b, kk;
	int nstar;
	double r[3], *m, val;
	gsl_matrix_view dfdy_mat;
	gsl_matrix *matrix;

	nstar = (*(fb_nonks_params_t *) params).nstar;
	m = (*(fb_nonks_params_t *) params).m;

	/* allocate matrices */
	dfdy_mat = gsl_matrix_view_array(dfdy, nstar*6, nstar*6);
	matrix = &dfdy_mat.matrix;

	/* set dfdt to zero */
	for (j=0; j< (unsigned int) nstar*6; j++) {
		dfdt[j] = 0.0;
	}

	/* set the matrix to zero to begin with */
	gsl_matrix_set_zero(matrix);

	/* then set the actual values */
	for (i=0; i< (unsigned int) nstar; i++) {
		for (a=0; a<3; a++) {
			gsl_matrix_set(matrix, i+a, i+a+3, 1.0);

			for (k=0; k< (unsigned int) nstar; k++) {
				for (b=0; b<3; b++) {
					val = 0.0;
					for (j=0; j< (unsigned int) nstar; j++) {
						if (j != i) {
							for (kk=0; kk<3; kk++) {
								r[kk] = y[i*6+kk] - y[j*6+kk];
							}
							val -= m[j] * ((double) (FB_DELTA(i, k) - FB_DELTA(j, k))) / fb_cub(fb_mod(r)) * \
								(((double) FB_DELTA(a, b)) - 3.0 * (y[i*6+a] - y[j*6+a])*(y[i*6+b] - y[j*6+b])/fb_dot(r, r));
						}
					}
					gsl_matrix_set(matrix, i+a+3, k+b, val);
				}
			}
		}
	}

	return(GSL_SUCCESS);
}

void fb_euclidean_to_nonks(fb_obj_t **star, double *y, int nstar)
{
	int i, j;

	for (i=0; i<nstar; i++) {
		for (j=0; j<3; j++) {
			y[i*6+j] = star[i]->x[j];
			y[i*6+j+3] = star[i]->v[j];
		}
	}
}

void fb_nonks_to_euclidean(double *y, fb_obj_t **star, int nstar)
{
	int i, j;

	for (i=0; i<nstar; i++) {
		for (j=0; j<3; j++) {
			star[i]->x[j] = y[i*6+j];
			star[i]->v[j] = y[i*6+j+3];
		}
	}
}

max_stars_repo_name: gnodvi/cosmos

max_issues_repo_name: gnodvi/cosmos

max_forks_repo_name: gnodvi/cosmos


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef uuid_93a5858f_7040_4f67_9767_808d2a4f7ff9_h
#define uuid_93a5858f_7040_4f67_9767_808d2a4f7ff9_h

#include <gslib/type.h>
#include <gslib/string.h>

__gslib_begin__

struct uuid_timestamp
{
    byte        tm_sec;
    byte        tm_min;
    byte        tm_hour;
    byte        tm_mday;    /* day of month */
    byte        tm_mon;
    byte        tm_wday;    /* day of week */
    int16       tm_year;
    int16       tm_yday;    /* day of year */
    int32       tm_fraction;
};

enum uuid_version
{
    uuid_v1,                /* date-time & MAC address */
    uuid_v2,                /* DCE security */
    uuid_v3,                /* MD5 hash & namespace */
    uuid_v4,                /* random */
    uuid_v5,                /* SHA-1 hash & namespace */
    uuid_ver_default = uuid_v4,
    uuid_ver_invalid = uuid_v2,
};

struct uuid_raw
{
    uint        data1;
    uint16      data2;
    uint16      data3;
    byte        data4[8];
};

class uuid:
    protected uuid_raw
{
public:
    uuid();
    uuid(uuid_version ver) { generate(ver); }
    uuid(uuid_version ver, const gchar* name, int len) { generate(ver, name, len); }
    uuid(const string& str) { from_string(str.c_str(), str.length()); }
    uuid(const gchar* str, int len) { from_string(str, len); }
    bool operator==(const uuid& that) const;
    bool is_valid() const;
    void generate(uuid_version ver);
    void generate(uuid_version ver, const gchar* name, int len);
    void from_string(const gchar* str, int len);
    const gchar* to_string(string& str) const;
    void get_timestamp(uuid_timestamp& ts) const;
};

__gslib_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib