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

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

#include "ccl.h"
#include "ccl_f2d.h"
#include "ccl_emu17.h"
#include "ccl_emu17_params.h"

// helper functions for BBKS and EH98
static double bbks_power(ccl_parameters *params, void *p, double k) {
  return ccl_bbks_power(params, k);
}

static double eh_power(ccl_parameters *params, void *p, double k) {
  return ccl_eh_power(params, (eh_struct*)p, k);
}

/*------ ROUTINE: ccl_cosmology_compute_power_analytic -----
INPUT: cosmology
TASK: provide spline for an analytic power spectrum with baryonic correction
*/

static ccl_f2d_t *ccl_compute_linpower_analytic(ccl_cosmology* cosmo, void* par,
                                                double (*pk)(ccl_parameters* params,
                                                             void* p, double k),
                                                int* status) {
  ccl_f2d_t *psp_out = NULL;
  double sigma8,log_sigma8;
  //These are the limits of the splining range
  double kmin = cosmo->spline_params.K_MIN;
  double kmax = cosmo->spline_params.K_MAX;
  //Compute nk from number of decades and N_K = # k per decade
  double ndecades = log10(kmax) - log10(kmin);
  int nk = (int)ceil(ndecades*cosmo->spline_params.N_K);
  // Compute na using predefined spline spacing
  double amin = cosmo->spline_params.A_SPLINE_MINLOG_PK;
  double amax = cosmo->spline_params.A_SPLINE_MAX;
  int na = cosmo->spline_params.A_SPLINE_NA_PK+cosmo->spline_params.A_SPLINE_NLOG_PK-1;

  // Exit if sigma8 wasn't specified
  if (isnan(cosmo->params.sigma8)) {
    *status = CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo,
             "ccl_power.c: ccl_compute_linpower_analytic(): "
             "sigma8 not set, required for analytic power spectra\n");
    return NULL;
  }

  // The x array is initially k, but will later
  // be overwritten with log(k)
  double *x=NULL, *y=NULL, *z=NULL, *y2d=NULL;
  x=ccl_log_spacing(kmin, kmax, nk);
  if(x==NULL) {
    *status = CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo,
             "ccl_power.c: ccl_compute_linpower_analytic(): "
             "memory allocation\n");
  }
  if(*status==0) {
    y=malloc(sizeof(double)*nk);
    if(y==NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo,
               "ccl_power.c: ccl_compute_linpower_analytic(): "
               "memory allocation\n");
    }
  }
  if(*status==0) {
    z=ccl_linlog_spacing(amin, cosmo->spline_params.A_SPLINE_MIN_PK,
       amax, cosmo->spline_params.A_SPLINE_NLOG_PK,
       cosmo->spline_params.A_SPLINE_NA_PK);
    if(z==NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo,
               "ccl_power.c: ccl_compute_linpower_analytic(): "
               "memory allocation\n");
    }
  }

  if(*status==0) {
    y2d = malloc(nk * na * sizeof(double));
    if(y2d==NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo,
               "ccl_power.c: ccl_compute_linpower_analytic(): "
               "memory allocation\n");
    }
  }

  if(*status==0) {
    // Calculate P(k) on k grid. After this loop, x will contain log(k) and y
    // will contain log(pk) [which has not yet been normalized]
    // After this loop x will contain log(k)
    for (int i=0; i<nk; i++) {
      y[i] = log((*pk)(&cosmo->params, par, x[i]));
      x[i] = log(x[i]);
    }
  }

  if(*status==0) {
    for (int j = 0; j < na; j++) {
      double gfac = ccl_growth_factor(cosmo,z[j], status);
      double g2 = 2.*log(gfac);
      for (int i=0; i<nk; i++) {
        y2d[j*nk+i] = y[i]+g2;
      }
    }
  }

  if(*status==0) {
    psp_out=ccl_f2d_t_new(na,z,nk,x,y2d,NULL,NULL,0,
                          1,2,ccl_f2d_cclgrowth,1,0,2,
                          ccl_f2d_3,status);
  }
  if(*status==0) {
    sigma8 = ccl_sigma8(cosmo, psp_out, status);
  }

  if(*status==0) {
    // Calculate normalization factor using computed value of sigma8, then
    // recompute P(k, a) using this normalization
    log_sigma8 = 2*(log(cosmo->params.sigma8) - log(sigma8));
    for(int j=0;j<na*nk;j++)
      y2d[j] += log_sigma8;
  }

  if(*status==0) {
    // Free the previous P(k,a) spline, and allocate a new one to store the
    // properly-normalized P(k,a)
    ccl_f2d_t_free(psp_out);
    psp_out = ccl_f2d_t_new(na,z,nk,x,y2d,NULL,NULL,0,
                            1,2,ccl_f2d_cclgrowth,1,0,2,
                            ccl_f2d_3,status);
  }

  free(x);
  free(y);
  free(z);
  free(y2d);
  return psp_out;
}

ccl_f2d_t *ccl_compute_linpower_bbks(ccl_cosmology *cosmo, int *status)
{
  ccl_f2d_t *psp=ccl_compute_linpower_analytic(cosmo, NULL, bbks_power, status);
  return psp;
}

ccl_f2d_t *ccl_compute_linpower_eh(ccl_cosmology *cosmo, int wiggled, int *status)
{
  ccl_f2d_t *psp = NULL;
  eh_struct *eh = NULL;
  eh = ccl_eh_struct_new(&(cosmo->params),wiggled);
  if (eh != NULL) {
    psp=ccl_compute_linpower_analytic(cosmo, eh,
                                      eh_power,
                                      status);
  }
  else
    *status = CCL_ERROR_MEMORY;
  free(eh);
  return psp;
}

/*------ ROUTINE: ccl_compute_power_emu -----
INPUT: cosmology
TASK: provide spline for the emulated power spectrum from Cosmic EMU
*/

ccl_f2d_t *ccl_compute_power_emu(ccl_cosmology * cosmo, int * status)
{
  double Omeganuh2_eq;
  ccl_f2d_t *psp_out=NULL;

  // Check ranges to see if the cosmology is valid
  if(*status==0) {
    if((cosmo->params.h<0.55) || (cosmo->params.h>0.85)){
      *status=CCL_ERROR_INCONSISTENT;
      ccl_cosmology_set_status_message(cosmo,
               "ccl_power.c: ccl_compute_power_emu(): "
               "h is outside allowed range\n");
    }
  }

  if(*status==0) {
   // Check if the cosmology has been set up with equal neutrino masses for the emulator
    // If not, check if the user has forced redistribution of masses and if so do this.
    if(cosmo->params.N_nu_mass>0) {
      if (cosmo->config.emulator_neutrinos_method == ccl_emu_strict){
        if (cosmo->params.N_nu_mass==3){
          if (cosmo->params.m_nu[0] != cosmo->params.m_nu[1] ||
              cosmo->params.m_nu[0] != cosmo->params.m_nu[2] ||
              cosmo->params.m_nu[1] != cosmo->params.m_nu[2]){
            *status = CCL_ERROR_INCONSISTENT;
            ccl_cosmology_set_status_message(cosmo,
                                             "ccl_power.c: ccl_compute_power_emu(): "
                                             "In the default configuration, you must pass a list of 3 "
                                             "equal neutrino masses or pass a sum and set "
                                             "m_nu_type = 'equal'. If you wish to over-ride this, "
                                             "set config->emulator_neutrinos_method = "
                                             "'ccl_emu_equalize'. This will force the neutrinos to "
                                             "be of equal mass but will result in "
                                             "internal inconsistencies.\n");
          }
        }else if (cosmo->params.N_nu_mass!=3){
          *status = CCL_ERROR_INCONSISTENT;
          ccl_cosmology_set_status_message(cosmo,
                                           "ccl_power.c: ccl_compute_power_emu(): "
                                           "In the default configuration, you must pass a list of 3 "
                                           "equal neutrino masses or pass a sum and set "
                                           "m_nu_type = 'equal'. If you wish to over-ride this, "
                                           "set config->emulator_neutrinos_method = "
                                           "'ccl_emu_equalize'. This will force the neutrinos to "
                                           "be of equal mass but will result in "
                                           "internal inconsistencies.\n");
        }
      }else if (cosmo->config.emulator_neutrinos_method == ccl_emu_equalize){
        // Reset the masses to equal
        double mnu_eq[3] = {cosmo->params.sum_nu_masses / 3.,
                            cosmo->params.sum_nu_masses / 3.,
                            cosmo->params.sum_nu_masses / 3.};
        Omeganuh2_eq = ccl_Omeganuh2(1.0, 3, mnu_eq, cosmo->params.T_CMB, status);
      }
    } else {
      if(fabs(cosmo->params.N_nu_rel - 3.04)>1.e-6){
        *status=CCL_ERROR_INCONSISTENT;
        ccl_cosmology_set_status_message(cosmo,
                                         "ccl_power.c: ccl_compute_power_emu(): "
                                         "Set Neff = 3.04 for cosmic emulator predictions in "
                                         "absence of massive neutrinos.\n");
      }
    }
  }

  if(*status==0) {
    double w0wacomb = -cosmo->params.w0 - cosmo->params.wa;
    if(w0wacomb<8.1e-3){ //0.3^4
      *status=CCL_ERROR_INCONSISTENT;
      ccl_cosmology_set_status_message(cosmo,
                                       "ccl_power.c: ccl_compute_power_emu(): "
                                       "w0 and wa do not satisfy the emulator bound\n");
    }
  }

  if(*status==0) {
    if(cosmo->params.Omega_nu_mass*cosmo->params.h*cosmo->params.h>0.01){
      *status=CCL_ERROR_INCONSISTENT;
      ccl_cosmology_set_status_message(cosmo,
                                       "ccl_power.c: ccl_compute_power_emu(): "
                                       "Omega_nu does not satisfy the emulator bound\n");
    }
  }

  if(*status==0) {
    // Check to see if sigma8 was defined
    if(isnan(cosmo->params.sigma8)){
      *status=CCL_ERROR_INCONSISTENT;
      ccl_cosmology_set_status_message(cosmo,
                                       "ccl_power.c: ccl_compute_power_emu(): "
                                       "sigma8 is not defined; specify sigma8 instead of A_s\n");
    }
  }

  int na=cosmo->spline_params.A_SPLINE_NA_PK;
  double *lpk_1a=NULL,*lk=NULL,*aemu=NULL,*lpk_nl=NULL;
  if (*status == 0) {
    //Now start the NL computation with the emulator
    //These are the limits of the splining range
    aemu = ccl_linear_spacing(A_MIN_EMU,cosmo->spline_params.A_SPLINE_MAX, na);
    if(aemu==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo,
                                       "ccl_power.c: ccl_compute_power_emu(): "
                                       "memory allocation error\n");
    }
  }
  if (*status == 0) {
    lk=malloc(NK_EMU*sizeof(double));
    if(lk==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo,
                                       "ccl_power.c: ccl_compute_power_emu(): "
                                       "memory allocation error\n");
    }
  }
  if (*status == 0) {
    //The emulator only computes power spectra at fixed nodes in k,
    //given by the global variable "mode"
    for (int i=0; i<NK_EMU; i++)
      lk[i] = log(mode[i]);
  }
  if (*status == 0) {
    lpk_nl = malloc(NK_EMU * na * sizeof(double));
    if(lpk_nl==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo,
                                       "ccl_power.c: ccl_compute_power_emu(): "
                                       "memory allocation error\n");
    }
  }
  if (*status == 0) {
    lpk_1a=malloc(NK_EMU*sizeof(double));
    if(lpk_1a==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo,
                                       "ccl_power.c: ccl_compute_power_emu(): "
                                       "memory allocation error\n");
    }
  }

  if (*status == 0) {
    double emu_par[9];
    //For each redshift:
    for (int j = 0; j < na; j++){
      //Turn cosmology into emu_par:
      if ((cosmo->params.N_nu_mass>0) &&
          (cosmo->config.emulator_neutrinos_method == ccl_emu_equalize)){
          emu_par[0] = (cosmo->params.Omega_c+cosmo->params.Omega_b)*cosmo->params.h*cosmo->params.h + Omeganuh2_eq;
      }else{
          emu_par[0] = (cosmo->params.Omega_c+cosmo->params.Omega_b+cosmo->params.Omega_nu_mass)*cosmo->params.h*cosmo->params.h;
      }
      emu_par[1] = cosmo->params.Omega_b*cosmo->params.h*cosmo->params.h;
      emu_par[2] = cosmo->params.sigma8;
      emu_par[3] = cosmo->params.h;
      emu_par[4] = cosmo->params.n_s;
      emu_par[5] = cosmo->params.w0;
      emu_par[6] = cosmo->params.wa;
      if ((cosmo->params.N_nu_mass>0) &&
          (cosmo->config.emulator_neutrinos_method == ccl_emu_equalize)){
        emu_par[7] = Omeganuh2_eq;
      }else{
        emu_par[7] = cosmo->params.Omega_nu_mass*cosmo->params.h*cosmo->params.h;
      }
      emu_par[8] = 1./aemu[j]-1;
      //Need to have this here because otherwise overwritten by emu in each loop

      //Call emulator at this redshift
      ccl_pkemu(emu_par,NK_EMU,lpk_1a, status, cosmo);
      if (*status) {
        *status=CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(cosmo,
                                         "ccl_power.c: ccl_compute_power_emu(): "
                                         "memory allocation error\n");
        break;
      }
      for (int i=0; i<NK_EMU; i++)
        lpk_nl[j*NK_EMU+i] = log(lpk_1a[i]);
    }
  }

  if(*status==0) {
    psp_out=ccl_f2d_t_new(na,aemu,NK_EMU,lk,lpk_nl,NULL,NULL,0,
                          1,2,ccl_f2d_no_extrapol,
                          1,0,2,ccl_f2d_3,status);
  }

  free(lpk_1a);
  free(lk);
  free(aemu);
  free(lpk_nl);
  return psp_out;
}

ccl_f2d_t *ccl_apply_halofit(ccl_cosmology* cosmo, ccl_f2d_t *plin, int *status)
{
  ccl_f2d_t *psp_out=NULL;
  size_t nk, na;
  double *x, *z, *y2d=NULL;

  //Halofit structure
  halofit_struct *hf=NULL;
  hf = ccl_halofit_struct_new(cosmo, plin, status);

  if(*status == 0) {
    //Find lk array
    if(plin->fk != NULL) {
      nk = plin->fk->size;
      x = plin->fk->x;
    }
    else {
      nk = plin->fka->interp_object.xsize;
      x = plin->fka->xarr;
    }

    //Find a array
    if(plin->fa != NULL) {
      na = plin->fa->size;
      z = plin->fa->x;
    }
    else {
      na = plin->fka->interp_object.ysize;
      z = plin->fka->yarr;
    }

    //Allocate pka array
    y2d = malloc(nk * na * sizeof(double));
    if (y2d == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo,
        "ccl_power.c: ccl_apply_halofit(): memory allocation\n");
    }
  }

  if (*status == 0) {
    // Calculate P(k) on a, k grid. After this loop, x will contain log(k) and y
    // will contain log(pk) [which has not yet been normalized]
    for (int j = 0; j<na; j++) {
      for (int i=0; i<nk; i++) {
        if (*status == 0) {
          double pk = ccl_halofit_power(cosmo, plin, x[i], z[j], hf, status);
          y2d[j*nk + i] = log(pk);
        }
      }
    }
  }

  if(*status == 0)
    psp_out = ccl_f2d_t_new(na, z, nk, x, y2d, NULL, NULL, 0,
                            1, 2, ccl_f2d_cclgrowth, 1,
                            0, 2, ccl_f2d_3, status);

  free(y2d);
  ccl_halofit_struct_free(hf);
  return psp_out;
}

void ccl_rescale_linpower(ccl_cosmology* cosmo, ccl_f2d_t *psp,
                          int rescale_mg, int rescale_norm,
                          int *status)
{
  if(rescale_mg || rescale_norm)
    ccl_rescale_musigma_s8(cosmo, psp, rescale_mg, status);
}

// Params for sigma(R) integrand
typedef struct {
  ccl_cosmology *cosmo;
  double R;
  double a;
  ccl_f2d_t *psp;
  int* status;
} SigmaR_pars;

typedef struct {
  ccl_cosmology *cosmo;
  double R;
  double a;
  ccl_f2d_t *psp;
  int* status;
} SigmaV_pars;

// Params for k_NL integrand
typedef struct {
  ccl_cosmology *cosmo;
  double a;
  ccl_f2d_t *psp;
  int* status;
} KNL_pars;

/* --------- ROUTINE: w_tophat ---------
INPUT: kR, ususally a wavenumber multiplied by a smoothing radius
TASK: Output W(x)=[sin(x)-x*cos(x)]*(3/x)^3
*/
static double w_tophat_2d(double kR) {
  double w;
  double kR2 = kR*kR;

  // This is the Maclaurin expansion of W(x)=2 J1(x)/x to O(x^10), with x=kR.
  // Necessary numerically because at low x W(x) relies on the fine cancellation of two terms
  if(kR<0.1) {
    w= 1. + kR2*(-1.0/8.0 + kR2*(1.0/192.0 +
      kR2*(-1.0/9216.0 + kR2*(1.0/737280.0 +
      kR2* (-1.0/88473600.0)))));
  }
  else
    w = 2 * gsl_sf_bessel_J1(kR) / kR;
  return w;
}

// Integrand for sigmaB integral (used for the SSC covariance calculation)
static double sigma2B_integrand(double lk,void *params) {
  SigmaR_pars *par=(SigmaR_pars *)params;

  double k=pow(10.,lk);
  double pk=ccl_f2d_t_eval(par->psp, lk * M_LN10, par->a,
                           par->cosmo, par->status);
  double kR=k*par->R;
  double w = w_tophat_2d(kR);

  return pk*k*k*w*w;
}

/* --------- ROUTINE: ccl_sigmaB ---------
INPUT: cosmology, comoving smoothing radius, scale factor
TASK: compute sigmaB, the variance in the projected *linear* density field
smoothed with a 2D tophat filter of comoving size R
*/
double ccl_sigma2B(ccl_cosmology *cosmo,double R,double a,ccl_f2d_t *psp, int *status)
{
  SigmaR_pars par;
  par.status = status;

  par.cosmo=cosmo;
  par.R=R;
  par.a=a;
  par.psp=psp;

  gsl_integration_cquad_workspace *workspace =  NULL;
  gsl_function F;
  F.function=&sigma2B_integrand;
  F.params=∥
  double sigma_B;

  workspace = gsl_integration_cquad_workspace_alloc(cosmo->gsl_params.N_ITERATION);
  if (workspace == NULL) {
    *status = CCL_ERROR_MEMORY;
  }
  if (*status == 0) {
    int gslstatus = gsl_integration_cquad(&F,
                                          log10(cosmo->spline_params.K_MIN),
                                          log10(cosmo->spline_params.K_MAX),
                                          0.0, cosmo->gsl_params.INTEGRATION_SIGMAR_EPSREL,
                                          workspace,&sigma_B,NULL,NULL);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_sigma2B():");
      *status |= gslstatus;
    }
  }
  gsl_integration_cquad_workspace_free(workspace);

  return sigma_B*M_LN10/(2*M_PI);
}

void ccl_sigma2Bs(ccl_cosmology *cosmo,int na, double *a, double *R,
                  double *sigma2B_out, ccl_f2d_t *psp, int *status) {
#pragma omp parallel default(none)                      \
  shared(cosmo, na, a, R, psp, sigma2B_out, status)
  {
    int ia;
    int local_status=*status;

    #pragma omp for
    for(ia=0; ia<na; ia++) {
      if(local_status==0)
        sigma2B_out[ia]=ccl_sigma2B(cosmo,R[ia],a[ia],psp,&local_status);
    } //end omp for
    if(local_status) {
      #pragma omp atomic write
      *status=local_status;
    }
  } //end omp parallel

  if(*status) {
    ccl_cosmology_set_status_message(cosmo,
             "ccl_power.c: ccl_sigma2Bs(): "
             "integration error\n");
  }
}

/* --------- ROUTINE: w_tophat ---------
INPUT: kR, ususally a wavenumber multiplied by a smoothing radius
TASK: Output W(x)=[sin(x)-x*cos(x)]*(3/x)^3
*/
static double w_tophat(double kR) {
  double w;
  double kR2 = kR*kR;

  // This is the Maclaurin expansion of W(x)=[sin(x)-xcos(x)]*3/x**3 to O(x^10), with x=kR.
  // Necessary numerically because at low x W(x) relies on the fine cancellation of two terms
  if(kR<0.1) {
    w= 1. + kR2*(-1.0/10.0 + kR2*(1.0/280.0 +
      kR2*(-1.0/15120.0 + kR2*(1.0/1330560.0 +
      kR2* (-1.0/172972800.0)))));
  }
  else
    w = 3.*(sin(kR) - kR*cos(kR))/(kR2*kR);
  return w;
}

// Integrand for sigmaR integral
static double sigmaR_integrand(double lk,void *params) {
  SigmaR_pars *par=(SigmaR_pars *)params;

  double k=pow(10.,lk);
  double pk=ccl_f2d_t_eval(par->psp, lk * M_LN10, par->a,
                           par->cosmo, par->status);
  double kR=k*par->R;
  double w = w_tophat(kR);

  return pk*k*k*k*w*w;
}

// Integrand for sigmaV integral
static double sigmaV_integrand(double lk,void *params) {
  SigmaV_pars *par=(SigmaV_pars *)params;

  double k=pow(10.,lk);
  double pk=ccl_f2d_t_eval(par->psp, lk * M_LN10, par->a,
                           par->cosmo, par->status);
  double kR=k*par->R;
  double w = w_tophat(kR);

  return pk*k*w*w/3.0;
}

/* --------- ROUTINE: ccl_sigmaR ---------
INPUT: cosmology, comoving smoothing radius, scale factor
TASK: compute sigmaR, the variance in the *linear* density field
smoothed with a tophat filter of comoving size R
*/
double ccl_sigmaR(ccl_cosmology *cosmo,double R,double a,ccl_f2d_t *psp, int *status) {

  SigmaR_pars par;
  par.status = status;

  par.cosmo=cosmo;
  par.R=R;
  par.a=a;
  par.psp=psp;

  gsl_integration_cquad_workspace *workspace =  NULL;
  gsl_function F;
  F.function=&sigmaR_integrand;
  F.params=∥
  double sigma_R;

  workspace = gsl_integration_cquad_workspace_alloc(cosmo->gsl_params.N_ITERATION);
  if (workspace == NULL) {
    *status = CCL_ERROR_MEMORY;
  }
  if (*status == 0) {
    int gslstatus = gsl_integration_cquad(&F,
                                          log10(cosmo->spline_params.K_MIN),
                                          log10(cosmo->spline_params.K_MAX),
                                          0.0, cosmo->gsl_params.INTEGRATION_SIGMAR_EPSREL,
                                          workspace,&sigma_R,NULL,NULL);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_sigmaR():");
      *status |= gslstatus;
    }
  }
  gsl_integration_cquad_workspace_free(workspace);

  return sqrt(sigma_R*M_LN10/(2*M_PI*M_PI));
}

/* --------- ROUTINE: ccl_sigmaV ---------
INPUT: cosmology, comoving smoothing radius, scale factor
TASK: compute sigmaV, the variance in the *linear* displacement field
smoothed with a tophat filter of comoving size R
The linear displacement field is the gradient of the linear density field
*/
double ccl_sigmaV(ccl_cosmology *cosmo,double R,double a,ccl_f2d_t *psp, int *status) {

  SigmaV_pars par;
  par.status = status;

  par.cosmo=cosmo;
  par.R=R;
  par.a=a;
  par.psp=psp;

  gsl_integration_cquad_workspace *workspace = NULL;
  gsl_function F;
  F.function=&sigmaV_integrand;
  F.params=∥
  double sigma_V;

  workspace = gsl_integration_cquad_workspace_alloc(cosmo->gsl_params.N_ITERATION);

  if (workspace == NULL) {
    *status = CCL_ERROR_MEMORY;
  }

  if (*status == 0) {
    int gslstatus = gsl_integration_cquad(&F,
                                          log10(cosmo->spline_params.K_MIN),
                                          log10(cosmo->spline_params.K_MAX),
                                          0.0, cosmo->gsl_params.INTEGRATION_SIGMAR_EPSREL,
                                          workspace,&sigma_V,NULL,NULL);

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

  gsl_integration_cquad_workspace_free(workspace);

  return sqrt(sigma_V*M_LN10/(2*M_PI*M_PI));
}

/* --------- ROUTINE: ccl_sigma8 ---------
INPUT: cosmology
TASK: compute sigma8, the variance in the *linear* density field at a=1
smoothed with a tophat filter of comoving size 8 Mpc/h
*/

double ccl_sigma8(ccl_cosmology *cosmo, ccl_f2d_t *psp, int *status) {
  return ccl_sigmaR(cosmo, 8/cosmo->params.h, 1., psp, status);
}

// Integrand for kNL integral
static double kNL_integrand(double k,void *params) {
  KNL_pars *par=(KNL_pars *)params;

  double pk=ccl_f2d_t_eval(par->psp, log(k), par->a,
                           par->cosmo, par->status);

  return pk;
}

/* --------- ROUTINE: ccl_kNL ---------
INPUT: cosmology, scale factor
TASK: compute kNL, the scale for the non-linear cut
*/
double ccl_kNL(ccl_cosmology *cosmo,double a,ccl_f2d_t *psp, int *status) {

  KNL_pars par;
  par.status = status;
  par.a = a;
  par.psp=psp;

  par.cosmo=cosmo;
  gsl_integration_cquad_workspace *workspace =  NULL;
  gsl_function F;
  F.function=&kNL_integrand;
  F.params=∥
  double PL_integral;

  workspace = gsl_integration_cquad_workspace_alloc(cosmo->gsl_params.N_ITERATION);
  if (workspace == NULL) {
    *status = CCL_ERROR_MEMORY;
  }
  if (*status == 0) {
    int gslstatus = gsl_integration_cquad(&F, cosmo->spline_params.K_MIN, cosmo->spline_params.K_MAX,
                                          0.0, cosmo->gsl_params.INTEGRATION_KNL_EPSREL,
                                          workspace,&PL_integral,NULL,NULL);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_kNL():");
      *status |= gslstatus;
    }
  }
  gsl_integration_cquad_workspace_free(workspace);
  double sigma_eta = sqrt(PL_integral/(6*M_PI*M_PI));
  return pow(sigma_eta, -1);
}

ext: c

max_stars_repo_name: Jappenn/CCL

max_issues_repo_name: Jappenn/CCL

max_forks_repo_name: Jappenn/CCL


text:
#include <math.h>
#include <unistd.h>
#include <stdio.h>
#include <stdbool.h>
#include <ctype.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
//#include <pthread.h>
#include <omp.h>
#include <complex.h>
#include <fftw3.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

#include "21cmFAST.h"
#include "exceptions.h"
#include "logger.h"
#include "Constants.h"
#include "Globals.h"
#include "indexing.c"
#include "UsefulFunctions.c"
#include "ps.c"
#include "dft.c"
#include "PerturbField.c"
#include "bubble_helper_progs.c"
#include "elec_interp.c"
#include "heating_helper_progs.c"
#include "recombinations.c"
#include "IonisationBox.c"
#include "SpinTemperatureBox.c"
#include "BrightnessTemperatureBox.c"
#include "FindHaloes.c"
#include "PerturbHaloField.c"

void adj_complex_conj(fftwf_complex *HIRES_box, struct UserParams *user_params, struct CosmoParams *cosmo_params){
    /*****  Adjust the complex conjugate relations for a real array  *****/

    int i, j, k;

    // corners
    HIRES_box[C_INDEX(0,0,0)] = 0;
    HIRES_box[C_INDEX(0,0,MIDDLE)] = crealf(HIRES_box[C_INDEX(0,0,MIDDLE)]);
    HIRES_box[C_INDEX(0,MIDDLE,0)] = crealf(HIRES_box[C_INDEX(0,MIDDLE,0)]);
    HIRES_box[C_INDEX(0,MIDDLE,MIDDLE)] = crealf(HIRES_box[C_INDEX(0,MIDDLE,MIDDLE)]);
    HIRES_box[C_INDEX(MIDDLE,0,0)] = crealf(HIRES_box[C_INDEX(MIDDLE,0,0)]);
    HIRES_box[C_INDEX(MIDDLE,0,MIDDLE)] = crealf(HIRES_box[C_INDEX(MIDDLE,0,MIDDLE)]);
    HIRES_box[C_INDEX(MIDDLE,MIDDLE,0)] = crealf(HIRES_box[C_INDEX(MIDDLE,MIDDLE,0)]);
    HIRES_box[C_INDEX(MIDDLE,MIDDLE,MIDDLE)] = crealf(HIRES_box[C_INDEX(MIDDLE,MIDDLE,MIDDLE)]);

    // do entire i except corners
#pragma omp parallel shared(HIRES_box) private(i,j,k) num_threads(user_params->N_THREADS)
    {
#pragma omp for
        for (i=1; i<MIDDLE; i++){
            // just j corners
            for (j=0; j<=MIDDLE; j+=MIDDLE){
                for (k=0; k<=MIDDLE; k+=MIDDLE){
                    HIRES_box[C_INDEX(i,j,k)] = conjf(HIRES_box[C_INDEX((user_params->DIM)-i,j,k)]);
                }
            }

            // all of j
            for (j=1; j<MIDDLE; j++){
                for (k=0; k<=MIDDLE; k+=MIDDLE){
                    HIRES_box[C_INDEX(i,j,k)] = conjf(HIRES_box[C_INDEX((user_params->DIM)-i,(user_params->DIM)-j,k)]);
                    HIRES_box[C_INDEX(i,(user_params->DIM)-j,k)] = conjf(HIRES_box[C_INDEX((user_params->DIM)-i,j,k)]);
                }
            }
        } // end loop over i
    }

    // now the i corners
#pragma omp parallel shared(HIRES_box) private(i,j,k) num_threads(user_params->N_THREADS)
    {
#pragma omp for
        for (i=0; i<=MIDDLE; i+=MIDDLE){
            for (j=1; j<MIDDLE; j++){
                for (k=0; k<=MIDDLE; k+=MIDDLE){
                    HIRES_box[C_INDEX(i,j,k)] = conjf(HIRES_box[C_INDEX(i,(user_params->DIM)-j,k)]);
                }
            }
        } // end loop over remaining j
    }
}

// Re-write of init.c for original 21cmFAST

int ComputeInitialConditions(
    unsigned long long random_seed, struct UserParams *user_params,
    struct CosmoParams *cosmo_params, struct InitialConditions *boxes
){

//     Generates the initial conditions: gaussian random density field (user_params->DIM^3) as well as the equal or lower resolution velocity fields, and smoothed density field (user_params->HII_DIM^3).
//
//     Author: Andrei Mesinger
//     Date: 9/29/06

    int status;

    Try{ // This Try wraps the entire function so we don't indent.

    // Makes the parameter structs visible to a variety of functions/macros
    // Do each time to avoid Python garbage collection issues
    Broadcast_struct_global_PS(user_params,cosmo_params);
    Broadcast_struct_global_UF(user_params,cosmo_params);

    unsigned long long ct;
    int n_x, n_y, n_z, i, j, k, ii, thread_num, dimension;
    float k_x, k_y, k_z, k_mag, p, a, b, k_sq;
    double pixel_deltax;
    float p_vcb, vcb_i;

    float f_pixel_factor;

    gsl_rng * r[user_params->N_THREADS];
    gsl_rng * rseed = gsl_rng_alloc(gsl_rng_mt19937); // An RNG for generating seeds for multithreading

    gsl_rng_set(rseed, random_seed);

    omp_set_num_threads(user_params->N_THREADS);

    switch(user_params->PERTURB_ON_HIGH_RES) {
        case 0:
            dimension = user_params->HII_DIM;
            break;
        case 1:
            dimension = user_params->DIM;
            break;
    }

    // ************  INITIALIZATION ********************** //
    unsigned int seeds[user_params->N_THREADS];

    // For multithreading, seeds for the RNGs are generated from an initial RNG (based on the input random_seed) and then shuffled (Author: Fred Davies)
    int num_int = INT_MAX/16;
    unsigned int *many_ints = (unsigned int *)malloc((size_t)(num_int*sizeof(unsigned int))); // Some large number of possible integers
    for (i=0; i<num_int; i++) {
        many_ints[i] = i;
    }

    gsl_ran_choose(rseed, seeds, user_params->N_THREADS, many_ints, num_int, sizeof(unsigned int)); // Populate the seeds array from the large list of integers
    gsl_ran_shuffle(rseed, seeds, user_params->N_THREADS, sizeof(unsigned int)); // Shuffle the randomly selected integers

    int checker;

    checker = 0;
    // seed the random number generators
    for (thread_num = 0; thread_num < user_params->N_THREADS; thread_num++){
        switch (checker){
            case 0:
                r[thread_num] = gsl_rng_alloc(gsl_rng_mt19937);
                gsl_rng_set(r[thread_num], seeds[thread_num]);
                break;
            case 1:
                r[thread_num] = gsl_rng_alloc(gsl_rng_gfsr4);
                gsl_rng_set(r[thread_num], seeds[thread_num]);
                break;
            case 2:
                r[thread_num] = gsl_rng_alloc(gsl_rng_cmrg);
                gsl_rng_set(r[thread_num], seeds[thread_num]);
                break;
            case 3:
                r[thread_num] = gsl_rng_alloc(gsl_rng_mrg);
                gsl_rng_set(r[thread_num], seeds[thread_num]);
                break;
            case 4:
                r[thread_num] = gsl_rng_alloc(gsl_rng_taus2);
                gsl_rng_set(r[thread_num], seeds[thread_num]);
                break;
        } // end switch

        checker += 1;

        if(checker==5) {
            checker = 0;
        }
    }

    free(many_ints);

    // allocate array for the k-space and real-space boxes
    fftwf_complex *HIRES_box = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);
    fftwf_complex *HIRES_box_saved = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);

    // allocate array for the k-space and real-space boxes for vcb
    fftwf_complex *HIRES_box_vcb_saved;
    // HIRES_box_vcb_saved may be needed if FFTW_Wisdom doesn't exist -- currently unused
    // but I am not going to allocate it until I am certain I needed it.

    // find factor of HII pixel size / deltax pixel size
    f_pixel_factor = user_params->DIM/(float)user_params->HII_DIM;

    // ************  END INITIALIZATION ****************** //
    LOG_DEBUG("Finished initialization.");
    // ************ CREATE K-SPACE GAUSSIAN RANDOM FIELD *********** //

    init_ps();

#pragma omp parallel shared(HIRES_box,r) \
                    private(n_x,n_y,n_z,k_x,k_y,k_z,k_mag,p,a,b,p_vcb) num_threads(user_params->N_THREADS)
    {
#pragma omp for
        for (n_x=0; n_x<user_params->DIM; n_x++){
            // convert index to numerical value for this component of the k-mode: k = (2*pi/L) * n
            if (n_x>MIDDLE)
                k_x =(n_x-user_params->DIM) * DELTA_K;  // wrap around for FFT convention
            else
                k_x = n_x * DELTA_K;

            for (n_y=0; n_y<user_params->DIM; n_y++){
                // convert index to numerical value for this component of the k-mode: k = (2*pi/L) * n
                if (n_y>MIDDLE)
                    k_y =(n_y-user_params->DIM) * DELTA_K;
                else
                    k_y = n_y * DELTA_K;

                // since physical space field is real, only half contains independent modes
                for (n_z=0; n_z<=MIDDLE; n_z++){
                    // convert index to numerical value for this component of the k-mode: k = (2*pi/L) * n
                    k_z = n_z * DELTA_K;

                    // now get the power spectrum; remember, only the magnitude of k counts (due to issotropy)
                    // this could be used to speed-up later maybe
                    k_mag = sqrt(k_x*k_x + k_y*k_y + k_z*k_z);
                    p = power_in_k(k_mag);

                    // ok, now we can draw the values of the real and imaginary part
                    // of our k entry from a Gaussian distribution
                    if(user_params->NO_RNG) {
                        a = 1.0;
                        b = -1.0;
                    }
                    else {
                        a = gsl_ran_ugaussian(r[omp_get_thread_num()]);
                        b = gsl_ran_ugaussian(r[omp_get_thread_num()]);
                    }

                    HIRES_box[C_INDEX(n_x, n_y, n_z)] = sqrt(VOLUME*p/2.0) * (a + b*I);

                }
            }
        }
    }
    LOG_DEBUG("Drawn random fields.");

    // *****  Adjust the complex conjugate relations for a real array  ***** //
    adj_complex_conj(HIRES_box,user_params,cosmo_params);

    memcpy(HIRES_box_saved, HIRES_box, sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);

    // FFT back to real space
    int stat = dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->DIM, user_params->N_THREADS, HIRES_box);
    if(stat>0) Throw(stat);
    LOG_DEBUG("FFT'd hires boxes.");

#pragma omp parallel shared(boxes,HIRES_box) private(i,j,k) num_threads(user_params->N_THREADS)
    {
#pragma omp for
        for (i=0; i<user_params->DIM; i++){
            for (j=0; j<user_params->DIM; j++){
                for (k=0; k<user_params->DIM; k++){
                    *((float *)boxes->hires_density + R_INDEX(i,j,k)) = *((float *)HIRES_box + R_FFT_INDEX(i,j,k))/VOLUME;
                }
            }
        }
    }

    // *** If required, let's also create a lower-resolution version of the density field  *** //
    memcpy(HIRES_box, HIRES_box_saved, sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);

    // Only filter if we are perturbing on the low-resolution grid
    if(!user_params->PERTURB_ON_HIGH_RES) {
        if (user_params->DIM != user_params->HII_DIM) {
            filter_box(HIRES_box, 0, 0, L_FACTOR*user_params->BOX_LEN/(user_params->HII_DIM+0.0));
        }

        // FFT back to real space
        dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->DIM, user_params->N_THREADS, HIRES_box);

        // Renormalise the FFT'd box (sample the high-res box if we are perturbing on the low-res grid)
#pragma omp parallel shared(boxes,HIRES_box,f_pixel_factor) private(i,j,k) num_threads(user_params->N_THREADS)
        {
#pragma omp for
            for (i=0; i<user_params->HII_DIM; i++){
                for (j=0; j<user_params->HII_DIM; j++){
                    for (k=0; k<user_params->HII_DIM; k++){
                        boxes->lowres_density[HII_R_INDEX(i,j,k)] =
                        *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                           (unsigned long long)(j*f_pixel_factor+0.5),
                                                           (unsigned long long)(k*f_pixel_factor+0.5)))/VOLUME;
                    }
                }
            }
        }
    }

    // ******* Relative Velocity part ******* //
  if(user_params->USE_RELATIVE_VELOCITIES){
    //JBM: We use the memory allocated to HIRES_box as it's free.

      for(ii=0;ii<3;ii++) {

        memcpy(HIRES_box, HIRES_box_saved, sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);

#pragma omp parallel shared(HIRES_box,ii) private(n_x,n_y,n_z,k_x,k_y,k_z,k_mag,p,p_vcb) num_threads(user_params->N_THREADS)
        {
#pragma omp for
            for (n_x=0; n_x<user_params->DIM; n_x++){
                if (n_x>MIDDLE)
                    k_x =(n_x-user_params->DIM) * DELTA_K;  // wrap around for FFT convention
                else
                    k_x = n_x * DELTA_K;

                for (n_y=0; n_y<user_params->DIM; n_y++){
                    if (n_y>MIDDLE)
                        k_y =(n_y-user_params->DIM) * DELTA_K;
                    else
                        k_y = n_y * DELTA_K;

                    for (n_z=0; n_z<=MIDDLE; n_z++){
                        k_z = n_z * DELTA_K;

                        k_mag = sqrt(k_x*k_x + k_y*k_y + k_z*k_z);
                        p = power_in_k(k_mag);
                        p_vcb = power_in_vcb(k_mag);

                        // now set the velocities
                        if ((n_x==0) && (n_y==0) && (n_z==0)){ // DC mode
                            HIRES_box[0] = 0;
                        }
                        else{
                            if(ii==0) {
                                HIRES_box[C_INDEX(n_x,n_y,n_z)] *= I * k_x/k_mag * sqrt(p_vcb/p) * C_KMS;
                            }
                            if(ii==1) {
                                HIRES_box[C_INDEX(n_x,n_y,n_z)] *= I * k_y/k_mag * sqrt(p_vcb/p) * C_KMS;
                            }
                            if(ii==2) {
                                HIRES_box[C_INDEX(n_x,n_y,n_z)] *= I * k_z/k_mag * sqrt(p_vcb/p) * C_KMS;
                            }
                        }
                    }
                }
            }
        }

//we only care about the lowres vcb box, so we filter it directly.
      if (user_params->DIM != user_params->HII_DIM) {
          filter_box(HIRES_box, 0, 0, L_FACTOR*user_params->BOX_LEN/(user_params->HII_DIM+0.0));
      }

//fft each velocity component back to real space
      dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->DIM, user_params->N_THREADS, HIRES_box);

      #pragma omp parallel shared(boxes,HIRES_box,f_pixel_factor,ii) private(i,j,k,vcb_i) num_threads(user_params->N_THREADS)
              {
      #pragma omp for
                  for (i=0; i<user_params->HII_DIM; i++){
                      for (j=0; j<user_params->HII_DIM; j++){
                          for (k=0; k<user_params->HII_DIM; k++){
                            vcb_i = *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                             (unsigned long long)(j*f_pixel_factor+0.5),
                                                             (unsigned long long)(k*f_pixel_factor+0.5)));
                            boxes->lowres_vcb[HII_R_INDEX(i,j,k)] += vcb_i*vcb_i;
                          }
                      }
                  }
              }

    }

//now we take the sqrt of that and normalize the FFT
    for (i=0; i<user_params->HII_DIM; i++){
        for (j=0; j<user_params->HII_DIM; j++){
            for (k=0; k<user_params->HII_DIM; k++){
              boxes->lowres_vcb[HII_R_INDEX(i,j,k)] = sqrt(boxes->lowres_vcb[HII_R_INDEX(i,j,k)])/VOLUME;
            }
        }
    }

  }
    LOG_DEBUG("Completed Relative velocities.");
    // ******* End of Relative Velocity part ******* //

    // Now look at the velocities

    for(ii=0;ii<3;ii++) {

        memcpy(HIRES_box, HIRES_box_saved, sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);
        // Now let's set the velocity field/dD/dt (in comoving Mpc)

#pragma omp parallel shared(HIRES_box,ii) private(n_x,n_y,n_z,k_x,k_y,k_z,k_sq) num_threads(user_params->N_THREADS)
        {
#pragma omp for
            for (n_x=0; n_x<user_params->DIM; n_x++){
                if (n_x>MIDDLE)
                    k_x =(n_x-user_params->DIM) * DELTA_K;  // wrap around for FFT convention
                else
                    k_x = n_x * DELTA_K;

                for (n_y=0; n_y<user_params->DIM; n_y++){
                    if (n_y>MIDDLE)
                        k_y =(n_y-user_params->DIM) * DELTA_K;
                    else
                        k_y = n_y * DELTA_K;

                    for (n_z=0; n_z<=MIDDLE; n_z++){
                        k_z = n_z * DELTA_K;

                        k_sq = k_x*k_x + k_y*k_y + k_z*k_z;

                        // now set the velocities
                        if ((n_x==0) && (n_y==0) && (n_z==0)){ // DC mode
                            HIRES_box[0] = 0;
                        }
                        else{
                            if(ii==0) {
                                HIRES_box[C_INDEX(n_x,n_y,n_z)] *= k_x*I/k_sq/VOLUME;
                            }
                            if(ii==1) {
                                HIRES_box[C_INDEX(n_x,n_y,n_z)] *= k_y*I/k_sq/VOLUME;
                            }
                            if(ii==2) {
                                HIRES_box[C_INDEX(n_x,n_y,n_z)] *= k_z*I/k_sq/VOLUME;
                            }
                        }
                    }
                }
            }
        }

        // Filter only if we require perturbing on the low-res grid
        if(!user_params->PERTURB_ON_HIGH_RES) {
            if (user_params->DIM != user_params->HII_DIM) {
                filter_box(HIRES_box, 0, 0, L_FACTOR*user_params->BOX_LEN/(user_params->HII_DIM+0.0));
            }
        }

        dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->DIM, user_params->N_THREADS, HIRES_box);

        // now sample to lower res
        // now sample the filtered box
#pragma omp parallel shared(boxes,HIRES_box,f_pixel_factor,ii,dimension) private(i,j,k) num_threads(user_params->N_THREADS)
        {
#pragma omp for
            for (i=0; i<dimension; i++){
                for (j=0; j<dimension; j++){
                    for (k=0; k<dimension; k++){
                        if(user_params->PERTURB_ON_HIGH_RES) {
                            if(ii==0) {
                                boxes->hires_vx[R_INDEX(i,j,k)] =
                                *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i),
                                                                   (unsigned long long)(j),
                                                                   (unsigned long long)(k)));
                            }
                            if(ii==1) {
                                boxes->hires_vy[R_INDEX(i,j,k)] =
                                *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i),
                                                                   (unsigned long long)(j),
                                                                   (unsigned long long)(k)));
                            }
                            if(ii==2) {
                                boxes->hires_vz[R_INDEX(i,j,k)] =
                                *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i),
                                                                   (unsigned long long)(j),
                                                                   (unsigned long long)(k)));
                            }
                        }
                        else {
                            if(ii==0) {
                                boxes->lowres_vx[HII_R_INDEX(i,j,k)] =
                                *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                                   (unsigned long long)(j*f_pixel_factor+0.5),
                                                                   (unsigned long long)(k*f_pixel_factor+0.5)));
                            }
                            if(ii==1) {
                                boxes->lowres_vy[HII_R_INDEX(i,j,k)] =
                                *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                                   (unsigned long long)(j*f_pixel_factor+0.5),
                                                                   (unsigned long long)(k*f_pixel_factor+0.5)));
                            }
                            if(ii==2) {
                                boxes->lowres_vz[HII_R_INDEX(i,j,k)] =
                                *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                                   (unsigned long long)(j*f_pixel_factor+0.5),
                                                                   (unsigned long long)(k*f_pixel_factor+0.5)));
                            }
                        }
                    }
                }
            }
        }
    }

    LOG_DEBUG("Done Inverse FT.");

    // * *************************************************** * //
    // *              BEGIN 2LPT PART                        * //
    // * *************************************************** * //

    // Generation of the second order Lagrangian perturbation theory (2LPT) corrections to the ZA
    // reference: Scoccimarro R., 1998, MNRAS, 299, 1097-1118 Appendix D

    // Parameter set in ANAL_PARAMS.H
    if(user_params->USE_2LPT){

        // use six supplementary boxes to store the gradients of phi_1 (eq. D13b)
        // Allocating the boxes
#define PHI_INDEX(i, j) ((int) ((i) - (j)) + 3*((j)) - ((int)(j))/2  )
        // ij -> INDEX
        // 00 -> 0
        // 11 -> 3
        // 22 -> 5
        // 10 -> 1
        // 20 -> 2
        // 21 -> 4

        fftwf_complex *phi_1 = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);

        // First generate the ii,jj phi_1 boxes

        int phi_component;

        float component_ii,component_jj,component_ij;

        // Indexing for the various phy components
        int phi_directions[3][2] = {{0,1},{0,2},{1,2}};

#pragma omp parallel shared(HIRES_box,phi_1) private(i,j,k) num_threads(user_params->N_THREADS)
        {
#pragma omp for
            for (i=0; i<user_params->DIM; i++){
                for (j=0; j<user_params->DIM; j++){
                    for (k=0; k<user_params->DIM; k++){
                        *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i),
                                                           (unsigned long long)(j),
                                                           (unsigned long long)(k)) ) = 0.;
                    }
                }
            }
        }

        // First iterate over the i = j components to phi
        // We'll also save these temporarily to the hires_vi_2LPT boxes which will get
        // overwritten later with the correct 2LPT velocities
        for(phi_component=0;phi_component<3;phi_component++) {

            i = j = phi_component;

                // generate the phi_1 boxes in Fourier transform
#pragma omp parallel shared(HIRES_box,phi_1,i,j) private(n_x,n_y,n_z,k_x,k_y,k_z,k_sq,k) num_threads(user_params->N_THREADS)
                {
#pragma omp for
                for (n_x=0; n_x<user_params->DIM; n_x++){
                    if (n_x>MIDDLE)
                        k_x =(n_x-user_params->DIM) * DELTA_K;  // wrap around for FFT convention
                    else
                        k_x = n_x * DELTA_K;

                    for (n_y=0; n_y<user_params->DIM; n_y++){
                        if (n_y>MIDDLE)
                            k_y =(n_y-user_params->DIM) * DELTA_K;
                        else
                            k_y = n_y * DELTA_K;

                        for (n_z=0; n_z<=MIDDLE; n_z++){
                            k_z = n_z * DELTA_K;

                            k_sq = k_x*k_x + k_y*k_y + k_z*k_z;

                            float k[] = {k_x, k_y, k_z};
                            // now set the velocities
                            if ((n_x==0) && (n_y==0) && (n_z==0)){ // DC mode
                                phi_1[0] = 0;
                            }
                            else{
                                phi_1[C_INDEX(n_x,n_y,n_z)] = -k[i]*k[j]*HIRES_box_saved[C_INDEX(n_x, n_y, n_z)]/k_sq/VOLUME;
                                // note the last factor of 1/VOLUME accounts for the scaling in real-space, following the FFT
                            }
                        }
                    }
                }
            }

            dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->DIM, user_params->N_THREADS, phi_1);

            // Temporarily store in the allocated hires_vi_2LPT boxes
#pragma omp parallel shared(boxes,phi_1,phi_component) private(i,j,k) num_threads(user_params->N_THREADS)
            {
#pragma omp for
                for (i=0; i<user_params->DIM; i++){
                    for (j=0; j<user_params->DIM; j++){
                        for (k=0; k<user_params->DIM; k++){
                            if(phi_component==0) {
                                boxes->hires_vx_2LPT[R_INDEX(i,j,k)] = *((float *)phi_1 + R_FFT_INDEX((unsigned long long)(i),
                                                                                (unsigned long long)(j),
                                                                                (unsigned long long)(k)));
                            }
                            if(phi_component==1) {
                                boxes->hires_vy_2LPT[R_INDEX(i,j,k)] = *((float *)phi_1 + R_FFT_INDEX((unsigned long long)(i),
                                                                                (unsigned long long)(j),
                                                                                (unsigned long long)(k)));
                            }
                            if(phi_component==2) {
                                boxes->hires_vz_2LPT[R_INDEX(i,j,k)] = *((float *)phi_1 + R_FFT_INDEX((unsigned long long)(i),
                                                                                (unsigned long long)(j),
                                                                                (unsigned long long)(k)));
                            }
                        }
                    }
                }
            }
        }

        for(phi_component=0;phi_component<3;phi_component++) {
            // Now calculate the cross components and start evaluating the 2LPT field
            i = phi_directions[phi_component][0];
            j = phi_directions[phi_component][1];

            // generate the phi_1 boxes in Fourier transform
#pragma omp parallel shared(HIRES_box,phi_1) private(n_x,n_y,n_z,k_x,k_y,k_z,k_sq,k) num_threads(user_params->N_THREADS)
            {
#pragma omp for
                for (n_x=0; n_x<user_params->DIM; n_x++){
                    if (n_x>MIDDLE)
                        k_x =(n_x-user_params->DIM) * DELTA_K;  // wrap around for FFT convention
                    else
                        k_x = n_x * DELTA_K;

                    for (n_y=0; n_y<user_params->DIM; n_y++){
                        if (n_y>MIDDLE)
                            k_y =(n_y-user_params->DIM) * DELTA_K;
                        else
                            k_y = n_y * DELTA_K;

                        for (n_z=0; n_z<=MIDDLE; n_z++){
                            k_z = n_z * DELTA_K;

                            k_sq = k_x*k_x + k_y*k_y + k_z*k_z;

                            float k[] = {k_x, k_y, k_z};
                            // now set the velocities
                            if ((n_x==0) && (n_y==0) && (n_z==0)){ // DC mode
                                phi_1[0] = 0;
                            }
                            else{
                                phi_1[C_INDEX(n_x,n_y,n_z)] = -k[i]*k[j]*HIRES_box_saved[C_INDEX(n_x, n_y, n_z)]/k_sq/VOLUME;
                                // note the last factor of 1/VOLUME accounts for the scaling in real-space, following the FFT
                            }
                        }
                    }
                }
            }

            dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->DIM, user_params->N_THREADS, phi_1);

            // Then we will have the laplacian of phi_2 (eq. D13b)
            // After that we have to return in Fourier space and generate the Fourier transform of phi_2
#pragma omp parallel shared(HIRES_box,phi_1,phi_component) private(i,j,k,component_ii,component_jj,component_ij) num_threads(user_params->N_THREADS)
            {
#pragma omp for
                for (i=0; i<user_params->DIM; i++){
                    for (j=0; j<user_params->DIM; j++){
                        for (k=0; k<user_params->DIM; k++){
                            // Note, I have temporarily stored the components into other arrays to minimise memory usage
                            // phi - {0, 1, 2} -> {hires_vx_2LPT, hires_vy_2LPT, hires_vz_2LPT}
                            // This may be opaque to the user, but this shouldn't need modification
                            if(phi_component==0) {
                                component_ii = boxes->hires_vx_2LPT[R_INDEX(i,j,k)];
                                component_jj = boxes->hires_vy_2LPT[R_INDEX(i,j,k)];
                                component_ij = *((float *)phi_1 + R_FFT_INDEX((unsigned long long)(i),
                                                                              (unsigned long long)(j),
                                                                              (unsigned long long)(k)));
                            }
                            if(phi_component==1) {
                                component_ii = boxes->hires_vx_2LPT[R_INDEX(i,j,k)];
                                component_jj = boxes->hires_vz_2LPT[R_INDEX(i,j,k)];
                                component_ij = *((float *)phi_1 + R_FFT_INDEX((unsigned long long)(i),
                                                                              (unsigned long long)(j),
                                                                              (unsigned long long)(k)));
                            }
                            if(phi_component==2) {
                                component_ii = boxes->hires_vy_2LPT[R_INDEX(i,j,k)];
                                component_jj = boxes->hires_vz_2LPT[R_INDEX(i,j,k)];
                                component_ij = *((float *)phi_1 + R_FFT_INDEX((unsigned long long)(i),
                                                                              (unsigned long long)(j),
                                                                              (unsigned long long)(k)));
                            }

                            // Kept in this form to maintain similar (possible) rounding errors
                            *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i),
                                                               (unsigned long long)(j),
                                                               (unsigned long long)(k)) ) += \
                            ( component_ii * component_jj );

                            *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i),
                                                               (unsigned long long)(j),
                                                               (unsigned long long)(k)) ) -= \
                            ( component_ij * component_ij );
                        }
                    }
                }
            }
        }

#pragma omp parallel shared(HIRES_box,phi_1) private(i,j,k) num_threads(user_params->N_THREADS)
        {
#pragma omp for
            for (i=0; i<user_params->DIM; i++){
                for (j=0; j<user_params->DIM; j++){
                    for (k=0; k<user_params->DIM; k++){
                        *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i),(unsigned long long)(j),(unsigned long long)(k)) ) /= TOT_NUM_PIXELS;
                    }
                }
            }
        }

        // Perform FFTs
        dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->DIM, user_params->N_THREADS, HIRES_box);

        memcpy(HIRES_box_saved, HIRES_box, sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);

        // Now we can store the content of box in a back-up array
        // Then we can generate the gradients of phi_2 (eq. D13b and D9)

        // ***** Store back-up k-box RHS eq. D13b ***** //

        // For each component, we generate the velocity field (same as the ZA part)

        // Now let's set the velocity field/dD/dt (in comoving Mpc)

        // read in the box
        // TODO correct free of phi_1

        for(ii=0;ii<3;ii++) {

            if(ii>0) {
                memcpy(HIRES_box, HIRES_box_saved, sizeof(fftwf_complex)*KSPACE_NUM_PIXELS);
            }

#pragma omp parallel shared(HIRES_box,ii) private(n_x,n_y,n_z,k_x,k_y,k_z,k_sq) num_threads(user_params->N_THREADS)
            {
#pragma omp for
            // set velocities/dD/dt
                for (n_x=0; n_x<user_params->DIM; n_x++){
                    if (n_x>MIDDLE)
                        k_x =(n_x-user_params->DIM) * DELTA_K;  // wrap around for FFT convention
                    else
                        k_x = n_x * DELTA_K;

                    for (n_y=0; n_y<user_params->DIM; n_y++){
                        if (n_y>MIDDLE)
                            k_y =(n_y-user_params->DIM) * DELTA_K;
                        else
                            k_y = n_y * DELTA_K;

                        for (n_z=0; n_z<=MIDDLE; n_z++){
                            k_z = n_z * DELTA_K;

                            k_sq = k_x*k_x + k_y*k_y + k_z*k_z;

                            // now set the velocities
                            if ((n_x==0) && (n_y==0) && (n_z==0)){ // DC mode
                                HIRES_box[0] = 0;
                            }
                            else{
                                if(ii==0) {
                                    HIRES_box[C_INDEX(n_x,n_y,n_z)] *= k_x*I/k_sq;
                                }
                                if(ii==1) {
                                    HIRES_box[C_INDEX(n_x,n_y,n_z)] *= k_y*I/k_sq;
                                }
                                if(ii==2) {
                                    HIRES_box[C_INDEX(n_x,n_y,n_z)] *= k_z*I/k_sq;
                                }
                            }
                        }
                        // note the last factor of 1/VOLUME accounts for the scaling in real-space, following the FFT
                    }
                }
            }

            // Filter only if we require perturbing on the low-res grid
            if(!user_params->PERTURB_ON_HIGH_RES) {
                if (user_params->DIM != user_params->HII_DIM) {
                    filter_box(HIRES_box, 0, 0, L_FACTOR*user_params->BOX_LEN/(user_params->HII_DIM+0.0));
                }
            }

            dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->DIM, user_params->N_THREADS, HIRES_box);

            // now sample to lower res
            // now sample the filtered box
#pragma omp parallel shared(boxes,HIRES_box,f_pixel_factor,ii,dimension) private(i,j,k) num_threads(user_params->N_THREADS)
            {
#pragma omp for
                for (i=0; i<dimension; i++){
                    for (j=0; j<dimension; j++){
                        for (k=0; k<dimension; k++){
                            if(user_params->PERTURB_ON_HIGH_RES) {
                                if(ii==0) {
                                    boxes->hires_vx_2LPT[R_INDEX(i,j,k)] =
                                    *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i),
                                                                       (unsigned long long)(j),
                                                                       (unsigned long long)(k)));
                                }
                                if(ii==1) {
                                    boxes->hires_vy_2LPT[R_INDEX(i,j,k)] =
                                    *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i),
                                                                       (unsigned long long)(j),
                                                                       (unsigned long long)(k)));
                                }
                                if(ii==2) {
                                    boxes->hires_vz_2LPT[R_INDEX(i,j,k)] =
                                    *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i),
                                                                       (unsigned long long)(j),
                                                                       (unsigned long long)(k)));
                                }
                            }
                            else {
                                if(ii==0) {
                                    boxes->lowres_vx_2LPT[HII_R_INDEX(i,j,k)] =
                                    *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                                       (unsigned long long)(j*f_pixel_factor+0.5),
                                                                       (unsigned long long)(k*f_pixel_factor+0.5)));
                                }
                                if(ii==1) {
                                    boxes->lowres_vy_2LPT[HII_R_INDEX(i,j,k)] =
                                    *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                                       (unsigned long long)(j*f_pixel_factor+0.5),
                                                                       (unsigned long long)(k*f_pixel_factor+0.5)));
                                }
                                if(ii==2) {
                                    boxes->lowres_vz_2LPT[HII_R_INDEX(i,j,k)] =
                                    *((float *)HIRES_box + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5),
                                                                       (unsigned long long)(j*f_pixel_factor+0.5),
                                                                       (unsigned long long)(k*f_pixel_factor+0.5)));
                                }
                            }
                        }
                    }
                }
            }
        }

        // deallocate the supplementary boxes
        fftwf_free(phi_1);

    }
    LOG_DEBUG("Done 2LPT.");

    // * *********************************************** * //
    // *               END 2LPT PART                     * //
    // * *********************************************** * //
    fftwf_cleanup_threads();
    fftwf_cleanup();
    fftwf_forget_wisdom();

    // deallocate
    fftwf_free(HIRES_box);
    fftwf_free(HIRES_box_saved);

    free_ps();

    for (i=0; i<user_params->N_THREADS; i++) {
        gsl_rng_free (r[i]);
    }
    gsl_rng_free(rseed);
    LOG_DEBUG("Cleaned Up.");
    } // End of Try{}

    Catch(status){
        return(status);
    }
    return(0);
}


text:
#pragma once

#include <array>

#include <gsl/span>

#include "Joystick.h"
#include "JoystickAxis.h"
#include "JoystickButton.h"
#include "JoystickProvider.h"

namespace qvr
{

class SfmlJoystick : public Joystick
{
public:
	struct ButtonState {
		bool isDown = false;
		bool wasDown = false;
	};

	struct AxisState {
		float position = 0.0f;
		float previousPosition = 0.0f;
	};

private:
	const int m_Index;

	int m_ButtonCount;
	int m_AxisCount;

	bool m_IsConnected  = false;
	bool m_WasConnected = false;

	// Copied from sf::Joystick
	enum
	{
		ButtonCount = 32, ///< Maximum number of supported buttons
		AxisCount = 8     ///< Maximum number of supported axes
	};

	std::array<ButtonState, ButtonCount> m_Buttons;
	std::array<AxisState,   AxisCount>   m_Axes;

	bool JustConnected() const {
		return m_IsConnected && !m_WasConnected;
	}
	bool JustDisconnected() const {
		return m_WasConnected && !m_IsConnected;
	}

	void UpdateButtons();
	void UpdateAxes();

public:
	SfmlJoystick(const int index) : m_Index(index) {}

	void Update();

	bool IsConnected() const {
		return m_IsConnected;
	}

	bool IsDown(const JoystickButton button) const override {
		return m_Buttons[button.GetValue()].isDown;
	}
	bool JustDown(const JoystickButton button) const override {
		return m_Buttons[button.GetValue()].isDown && !m_Buttons[button.GetValue()].wasDown;
	}
	bool JustUp(const JoystickButton button) const override {
		return m_Buttons[button.GetValue()].wasDown && !m_Buttons[button.GetValue()].isDown;
	}

	float GetPosition(const JoystickAxis axis) const override {
		return m_Axes[(int)axis].position;
	}
	float GetPreviousPosition(const JoystickAxis axis) const override {
		return m_Axes[(int)axis].previousPosition;
	}

	const gsl::span<const ButtonState> GetButtons() const {
		return gsl::make_span(m_Buttons.data(), m_ButtonCount);
	}

	const gsl::span<const AxisState> GetAxes() const {
		return gsl::make_span(m_Axes);
	}

	int GetButtonCount() const {
		return m_ButtonCount;
	}

	int GetAxisCount() const {
		return m_AxisCount;
	}
};

class SfmlJoystickSet : public JoystickProvider
{
public:
	// Worst constructor ever.
	SfmlJoystickSet()
		: m_Joysticks({ 0, 1, 2, 3, 4, 5, 6, 7 })
	{}

	void Update();

	const Joystick* GetJoystick(const JoystickIndex index) const override {
		if (m_Joysticks[index.Get()].IsConnected()) {
			return &m_Joysticks[index.Get()];
		}
		return nullptr;
	}

	static const int Count = 8; // The maximum number of supported joysticks.

private:
	std::array<SfmlJoystick, Count> m_Joysticks;
};

void LogSfmlJoystickInfo(const int index);

}

ext: h

max_stars_repo_name: rachelnertia/Quarrel

max_issues_repo_name: rachelnertia/Quarrel

max_forks_repo_name: rachelnertia/Quiver


text:
#ifndef MATCH_H_COIEXYNJ
#define MATCH_H_COIEXYNJ

#include <gsl/gsl>
#include <sens_loc/analysis/recognition_performance.h>
#include <sens_loc/math/pointcloud.h>
#include <utility>
#include <vector>

namespace sens_loc::analysis {

/// Create linear vectors that pair the keypoints from \c kp0 to \c kp1
/// according to \c matches.
///
/// \pre \c query is the query set in \c matches.
/// \pre \c train is the training set in \c matches.
/// \param matches Set of matches in \c query and \c train
/// \param query Set of keypoints that were matched as the query set.
/// \param train Set of keypoints that were matched as the training set.
/// \returns \c pair{QueriedKps, TrainedKps} were each vector contains only
/// matched keypoints. All other keypoints are discarded in the result.
std::pair<math::imagepoints_t, math::imagepoints_t>
gather_correspondences(gsl::span<const keypoint_correspondence> matches,
                       const math::imagepoints_t&               query,
                       const math::imagepoints_t&               train) noexcept;
}  // namespace sens_loc::analysis

#endif /* end of include guard: MATCH_H_COIEXYNJ */

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


text:
#ifndef PLV_DAL_BLOCKCHAIN_STORAGE_H
#define PLV_DAL_BLOCKCHAIN_STORAGE_H

#include <plv/io/file.h>
#include <plv/serializers/netstring.h>
#include <gsl/gsl_assert>
#include <memory>
#include <vector>
#include "chunkiterator.h"
#include "metablock.h"
#include "storage.h"

namespace plv::dal
{
class BlockchainStorage : public Storage
{
public:
    BlockchainStorage(std::shared_ptr<ChunkIterator> chunkIter, std::shared_ptr<io::File> file)
        : file_(file), chunkIterator_(chunkIter)
    {
    }

    auto indexStorage() -> void override
    {
        using namespace serializers::netstring;
        metaBlocks_.clear();
        auto chunkOpt = chunkIterator_->getNext();
        while (chunkOpt.has_value())
        {
            uint32_t objectIdx = 0;
            do
            {
                auto chunk                     = chunkOpt->getData().substr(objectIdx);
                auto [blockStart, blockLength] = decodeOffsetAndLength(chunk);
                metaBlocks_.emplace_back(chunkOpt->getOffset() + objectIdx + blockStart, blockLength);
                objectIdx += calculateOffsetOfNextObject(blockStart, blockLength);
            } while (objectIdx < chunkOpt->getChunkSize() - 1);
            chunkOpt = chunkIterator_->getNext();
        }
    }

    auto getBlock(Id id) -> Buffer override
    {
        Expects(id < metaBlocks_.size());
        auto metaBlock = metaBlocks_[id];
        Buffer buffer(metaBlock.blockSize, 0);
        file_->read(buffer, metaBlock.blockOffset, metaBlock.blockSize);
        return buffer;
    }

    auto addBlock(BufferView block) -> void override
    {
        using namespace serializers::netstring;
        auto blockSerialized           = encode(block);
        auto [blockStart, blockLength] = decodeOffsetAndLength(blockSerialized);
        metaBlocks_.emplace_back(file_->getSize() + blockStart, blockLength);
        file_->write(blockSerialized);
    }
    /*
        auto getBlock(const datamodel::MetaBlock& metaBlock) -> datamodel::Block
        {
            datamodel::Block block{};
            Buffer buffer(metaBlock.blockSize, 0);
            file_->read(buffer, metaBlock.blockOffset, metaBlock.blockSize);
            auto [hash, consumed] = netstring::decode(buffer);
            std::copy(std::begin(hash), std::end(hash), block.prevBlockHash.data());
            auto remaining = metaBlock.blockSize - consumed;
            while (remaining > 0)
            {
                auto nextData                  = buffer.substr(consumed, remaining);
                auto [dataElement, consumedDE] = netstring::decode(nextData);
                block.dataElements.emplace_back(dataElement);
                remaining -= consumedDE;
                consumed += consumedDE;
            }
            return block;
        }

        auto getBlockPrevHash(const datamodel::MetaBlock& metaBlock) -> datamodel::Block
        {
            datamodel::Block block{};
            Buffer buffer(metaBlock.blockSize, 0);
            file_->read(buffer, metaBlock.blockOffset, metaBlock.blockSize);
            auto [hash, consumed] = netstring::decode(buffer);
            std::copy(std::begin(hash), std::end(hash), block.prevBlockHash.data());
            auto remaining = metaBlock.blockSize - consumed;
            while (remaining > 0)
            {
                auto nextData                  = buffer.substr(consumed, remaining);
                auto [dataElement, consumedDE] = netstring::decode(nextData);
                block.dataElements.emplace_back(dataElement);
                remaining -= consumedDE;
                consumed += consumedDE;
            }
            return block;
        }
    */
    auto getChainLength() -> Size override
    {
        return metaBlocks_.size();
    }

    auto getStorageSize() -> Size override
    {
        return file_->getSize();
    }

    auto isIndexed() -> bool override
    {
        return !(metaBlocks_.size() == 0 && file_->getSize() != 0);
    }

private:
    std::shared_ptr<io::File> file_{};
    std::shared_ptr<ChunkIterator> chunkIterator_{};
    std::vector<dal::MetaBlock> metaBlocks_{};
    bool indexed_{false};
};
}  // namespace plv::dal
#endif  // PLV_DAL_BLOCKCHAIN_STORAGE_H

max_stars_repo_name: tymion/blockchain

max_issues_repo_name: tymion/blockchain

max_forks_repo_name: tymion/blockchain


text:
/* randist/dirichlet.c
 *
 * Copyright (C) 2007 Brian Gough
 * Copyright (C) 2002 Gavin E. Crooks <gec@compbio.berkeley.edu>
 *
 * 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 <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>

/* The Dirichlet probability distribution of order K-1 is

     p(\theta_1,...,\theta_K) d\theta_1 ... d\theta_K =
        (1/Z) \prod_i=1,K \theta_i^{alpha_i - 1} \delta(1 -\sum_i=1,K \theta_i)

   The normalization factor Z can be expressed in terms of gamma functions:

      Z = {\prod_i=1,K \Gamma(\alpha_i)} / {\Gamma( \sum_i=1,K \alpha_i)}

   The K constants, \alpha_1,...,\alpha_K, must be positive. The K parameters,
   \theta_1,...,\theta_K are nonnegative and sum to 1.

   The random variates are generated by sampling K values from gamma
   distributions with parameters a=\alpha_i, b=1, and renormalizing.
   See A.M. Law, W.D. Kelton, Simulation Modeling and Analysis (1991).

   Gavin E. Crooks <gec@compbio.berkeley.edu> (2002)
*/

static void ran_dirichlet_small (const gsl_rng * r, const size_t K, const double alpha[], double theta[]);

void
gsl_ran_dirichlet (const gsl_rng * r, const size_t K,
                   const double alpha[], double theta[])
{
  size_t i;
  double norm = 0.0;

  for (i = 0; i < K; i++)
    {
      theta[i] = gsl_ran_gamma (r, alpha[i], 1.0);
    }

  for (i = 0; i < K; i++)
    {
      norm += theta[i];
    }

  if (norm < GSL_SQRT_DBL_MIN)  /* Handle underflow */
    {
      ran_dirichlet_small (r, K, alpha, theta);
      return;
    }

  for (i = 0; i < K; i++)
    {
      theta[i] /= norm;
    }
}

/* When the values of alpha[] are small, scale the variates to avoid
   underflow so that the result is not 0/0.  Note that the Dirichlet
   distribution is defined by a ratio of gamma functions so we can
   take out an arbitrary factor to keep the values in the range of
   double precision. */

static void
ran_dirichlet_small (const gsl_rng * r, const size_t K,
                     const double alpha[], double theta[])
{
  size_t i;
  double norm = 0.0, umax = 0;

  for (i = 0; i < K; i++)
    {
      double u = log(gsl_rng_uniform_pos (r)) / alpha[i];

      theta[i] = u;

      if (u > umax || i == 0) {
        umax = u;
      }
    }

  for (i = 0; i < K; i++)
    {
      theta[i] = exp(theta[i] - umax);
    }

  for (i = 0; i < K; i++)
    {
      theta[i] = theta[i] * gsl_ran_gamma (r, alpha[i] + 1.0, 1.0);
    }

  for (i = 0; i < K; i++)
    {
      norm += theta[i];
    }

  for (i = 0; i < K; i++)
    {
      theta[i] /= norm;
    }
}

double
gsl_ran_dirichlet_pdf (const size_t K,
                       const double alpha[], const double theta[])
{
  return exp (gsl_ran_dirichlet_lnpdf (K, alpha, theta));
}

double
gsl_ran_dirichlet_lnpdf (const size_t K,
                         const double alpha[], const double theta[])
{
  /*We calculate the log of the pdf to minimize the possibility of overflow */
  size_t i;
  double log_p = 0.0;
  double sum_alpha = 0.0;

  for (i = 0; i < K; i++)
    {
      log_p += (alpha[i] - 1.0) * log (theta[i]);
    }

  for (i = 0; i < K; i++)
    {
      sum_alpha += alpha[i];
    }

  log_p += gsl_sf_lngamma (sum_alpha);

  for (i = 0; i < K; i++)
    {
      log_p -= gsl_sf_lngamma (alpha[i]);
    }

  return log_p;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include <map>
#include "DrawableGameComponent.h"
#include "FullScreenRenderTarget.h"
#include "FullScreenQuad.h"
#include "MatrixHelper.h"

namespace Library
{
	class PixelShader;
}

namespace Rendering
{
	class DiffuseLightingDemo;

	enum class ColorFilters
	{
		GrayScale = 0,
		Inverse,
		Sepia,
		Generic,
		End
	};

	class ColorFilteringDemo final : public Library::DrawableGameComponent
	{
	public:
		ColorFilteringDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		ColorFilteringDemo(const ColorFilteringDemo&) = delete;
		ColorFilteringDemo(ColorFilteringDemo&&) = default;
		ColorFilteringDemo& operator=(const ColorFilteringDemo&) = default;
		ColorFilteringDemo& operator=(ColorFilteringDemo&&) = default;
		~ColorFilteringDemo();

		std::shared_ptr<DiffuseLightingDemo> DiffuseLighting() const;

		ColorFilters ActiveColorFilter() const;
		void SetActiveColorFilter(ColorFilters colorFilter);

		float GenericFilterBrightness() const;
		void SetGenericFilterBrightness(float brightness);

		static const std::map<ColorFilters, std::string> ColorFilterNames;

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		struct GenericColorFilterPSConstantBuffer
		{
			DirectX::XMFLOAT4X4 ColorFilter{ Library::MatrixHelper::Identity };
		};

		std::shared_ptr<DiffuseLightingDemo> mDiffuseLightingDemo;
		Library::FullScreenRenderTarget mRenderTarget;
		Library::FullScreenQuad mFullScreenQuad;
		std::map<ColorFilters, std::shared_ptr<Library::PixelShader>> mPixelShadersByColorFilter;
		ColorFilters mActiveColorFilter;
		winrt::com_ptr<ID3D11Buffer> mGenericColorFilterPSConstantBuffer;
		GenericColorFilterPSConstantBuffer mGenericColorFilterPSConstantBufferData;
	};
}


text:
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.

#pragma once

#include <limits>
#include <cassert>
#include <chrono>
#include <vector>
#include <map>
#include <set>
#include <numeric>

#include <wrl/client.h>
#include <wrl/implements.h>

#include <wil/wrl.h>
#include <wil/result.h>

#include <gsl/gsl>

max_stars_repo_name: dennyac/onnxruntime

max_issues_repo_name: dennyac/onnxruntime

max_forks_repo_name: dennyac/onnxruntime


text:
///////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2015 Microsoft Corporation. All rights reserved.
//
// This code is licensed under the MIT License (MIT).
//
// 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.
//
///////////////////////////////////////////////////////////////////////////////

#ifndef GSL_UTIL_H
#define GSL_UTIL_H

#include <gsl/gsl_assert> // for Expects

#include <array>
#include <cstddef>          // for ptrdiff_t, size_t
#include <exception>        // for exception
#include <initializer_list> // for initializer_list
#include <type_traits>      // for is_signed, integral_constant
#include <utility>          // for forward

#if defined(_MSC_VER)

#pragma warning(push)
#pragma warning(disable : 4127) // conditional expression is constant

#if _MSC_VER < 1910
#pragma push_macro("constexpr")
#define constexpr /*constexpr*/
#endif            // _MSC_VER < 1910
#endif            // _MSC_VER

namespace gsl
{
//
// GSL.util: utilities
//

// index type for all container indexes/subscripts/sizes
using index = std::ptrdiff_t;

// final_action allows you to ensure something gets run at the end of a scope
template <class F>
class final_action
{
public:
    explicit final_action(F f) noexcept : f_(std::move(f)) {}

    final_action(final_action&& other) noexcept : f_(std::move(other.f_)), invoke_(other.invoke_)
    {
        other.invoke_ = false;
    }

    final_action(const final_action&) = delete;
    final_action& operator=(const final_action&) = delete;
    final_action& operator=(final_action&&) = delete;

    GSL_SUPPRESS(f.6) // NO-FORMAT: attribute // terminate if throws
    ~final_action() noexcept
    {
        if (invoke_) f_();
    }

private:
    F f_;
    bool invoke_{true};
};

// finally() - convenience function to generate a final_action
template <class F>
final_action<F> finally(const F& f) noexcept
{
    return final_action<F>(f);
}

template <class F>
final_action<F> finally(F&& f) noexcept
{
    return final_action<F>(std::forward<F>(f));
}

// narrow_cast(): a searchable way to do narrowing casts of values
template <class T, class U>
GSL_SUPPRESS(type.1) // NO-FORMAT: attribute
constexpr T narrow_cast(U&& u) noexcept
{
    return static_cast<T>(std::forward<U>(u));
}

struct narrowing_error : public std::exception
{
};

namespace details
{
    template <class T, class U>
    struct is_same_signedness
        : public std::integral_constant<bool, std::is_signed<T>::value == std::is_signed<U>::value>
    {
    };
} // namespace details

// narrow() : a checked version of narrow_cast() that throws if the cast changed the value
template <class T, class U>
GSL_SUPPRESS(type.1) // NO-FORMAT: attribute
GSL_SUPPRESS(f.6) // NO-FORMAT: attribute // TODO: MSVC /analyze does not recognise noexcept(false)
T narrow(U u) noexcept(false)
{
    T t = narrow_cast<T>(u);
    if (static_cast<U>(t) != u) gsl::details::throw_exception(narrowing_error());
    if (!details::is_same_signedness<T, U>::value && ((t < T{}) != (u < U{})))
        gsl::details::throw_exception(narrowing_error());
    return t;
}

//
// at() - Bounds-checked way of accessing builtin arrays, std::array, std::vector
//
template <class T, std::size_t N>
GSL_SUPPRESS(bounds.4) // NO-FORMAT: attribute
GSL_SUPPRESS(bounds.2) // NO-FORMAT: attribute
constexpr T& at(T (&arr)[N], const index i)
{
    Expects(i >= 0 && i < narrow_cast<index>(N));
    return arr[narrow_cast<std::size_t>(i)];
}

template <class Cont>
GSL_SUPPRESS(bounds.4) // NO-FORMAT: attribute
GSL_SUPPRESS(bounds.2) // NO-FORMAT: attribute
constexpr auto at(Cont& cont, const index i) -> decltype(cont[cont.size()])
{
    Expects(i >= 0 && i < narrow_cast<index>(cont.size()));
    using size_type = decltype(cont.size());
    return cont[narrow_cast<size_type>(i)];
}

template <class T>
GSL_SUPPRESS(bounds.1) // NO-FORMAT: attribute
constexpr T at(const std::initializer_list<T> cont, const index i)
{
    Expects(i >= 0 && i < narrow_cast<index>(cont.size()));
    return *(cont.begin() + i);
}

} // namespace gsl

#if defined(_MSC_VER)
#if _MSC_VER < 1910
#undef constexpr
#pragma pop_macro("constexpr")

#endif // _MSC_VER < 1910

#pragma warning(pop)

#endif // _MSC_VER

#endif // GSL_UTIL_H

max_stars_repo_name: Redchards/CVRP

max_issues_repo_name: Redchards/CVRP

max_forks_repo_name: Redchards/CVRP


text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpRNG.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void set_MCMC_covariance(MCMC_info *MCMC, double *V) {
    int i_P, j_P;
    if(V == NULL)
        SID_log("Initializing the covariance matrix...", SID_LOG_OPEN);
    else
        SID_log("Updating the covariance matrix...", SID_LOG_OPEN);
    if(MCMC->V == NULL)
        MCMC->V = (double *)SID_malloc(sizeof(double) * MCMC->n_P * MCMC->n_P);
    if(MCMC->m == NULL)
        MCMC->m = gsl_matrix_calloc(MCMC->n_P, MCMC->n_P);
    if(MCMC->b == NULL)
        MCMC->b = gsl_vector_calloc(MCMC->n_P);
    if(V == NULL) {
        for(i_P = 0; i_P < MCMC->n_P; i_P++) {
            for(j_P = 0; j_P < MCMC->n_P; j_P++) {
                if(i_P == j_P) {
                    // Matrix decomposition fails if we initialize an element to zero here.
                    //   Check for this and start with an different value if true.
                    if(MCMC->P_init[i_P] == 0.)
                        MCMC->V[i_P * MCMC->n_P + j_P] = 1e-3 * MCMC->P_limit_max[i_P];
                    else
                        MCMC->V[i_P * MCMC->n_P + j_P] = pow(MCMC->P_init[i_P], 2.);
                } else
                    MCMC->V[i_P * MCMC->n_P + j_P] = 0.;
            }
        }
    } else
        memcpy(MCMC->V, V, (size_t)(MCMC->n_P * MCMC->n_P) * sizeof(double));
    for(i_P = 0; i_P < MCMC->n_P; i_P++)
        for(j_P = 0; j_P < MCMC->n_P; j_P++)
            gsl_matrix_set(MCMC->m, i_P, j_P, MCMC->V[i_P * MCMC->n_P + j_P]);
    gsl_linalg_cholesky_decomp(MCMC->m);

    // Set the upper-diagonal to zero to get L where covariance=L*L^T
    for(i_P = 0; i_P < MCMC->n_P; i_P++)
        for(j_P = i_P + 1; j_P < MCMC->n_P; j_P++)
            gsl_matrix_set(MCMC->m, i_P, j_P, 0.);
    SID_log("Done.", SID_LOG_CLOSE);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
/* specfunc/gsl_sf_mathieu.h
 *
 * Copyright (C) 2002 Lowell Johnson
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  L. Johnson */

#ifndef __GSL_SF_MATHIEU_H__
#define __GSL_SF_MATHIEU_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_eigen.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

#define GSL_SF_MATHIEU_COEFF 100

typedef struct
{
  size_t size;
  size_t even_order;
  size_t odd_order;
  int extra_values;
  double qa;   /* allow for caching of results: not implemented yet */
  double qb;   /* allow for caching of results: not implemented yet */
  double *aa;
  double *bb;
  double *dd;
  double *ee;
  double *tt;
  double *e2;
  double *zz;
  gsl_vector *eval;
  gsl_matrix *evec;
  gsl_eigen_symmv_workspace *wmat;
} gsl_sf_mathieu_workspace;

/* Compute an array of characteristic (eigen) values from the recurrence
   matrices for the Mathieu equations. */
GSL_FUN int gsl_sf_mathieu_a_array(int order_min, int order_max, double qq, gsl_sf_mathieu_workspace *work, double result_array[]);
GSL_FUN int gsl_sf_mathieu_b_array(int order_min, int order_max, double qq,  gsl_sf_mathieu_workspace *work, double result_array[]);

/* Compute the characteristic value for a Mathieu function of order n and
   type ntype. */
GSL_FUN int gsl_sf_mathieu_a_e(int order, double qq, gsl_sf_result *result);
GSL_FUN double gsl_sf_mathieu_a(int order, double qq);
GSL_FUN int gsl_sf_mathieu_b_e(int order, double qq, gsl_sf_result *result);
GSL_FUN double gsl_sf_mathieu_b(int order, double qq);

/* Compute the Fourier coefficients for a Mathieu function. */
GSL_FUN int gsl_sf_mathieu_a_coeff(int order, double qq, double aa, double coeff[]);
GSL_FUN int gsl_sf_mathieu_b_coeff(int order, double qq, double aa, double coeff[]);

/* Allocate computational storage space for eigenvalue solution. */
GSL_FUN gsl_sf_mathieu_workspace *gsl_sf_mathieu_alloc(const size_t nn,
                                               const double qq);
GSL_FUN void gsl_sf_mathieu_free(gsl_sf_mathieu_workspace *workspace);

/* Compute an angular Mathieu function. */
GSL_FUN int gsl_sf_mathieu_ce_e(int order, double qq, double zz, gsl_sf_result *result);
GSL_FUN double gsl_sf_mathieu_ce(int order, double qq, double zz);
GSL_FUN int gsl_sf_mathieu_se_e(int order, double qq, double zz, gsl_sf_result *result);
GSL_FUN double gsl_sf_mathieu_se(int order, double qq, double zz);
GSL_FUN int gsl_sf_mathieu_ce_array(int nmin, int nmax, double qq, double zz,
                            gsl_sf_mathieu_workspace *work,
                            double result_array[]);
GSL_FUN int gsl_sf_mathieu_se_array(int nmin, int nmax, double qq, double zz,
                            gsl_sf_mathieu_workspace *work,
                            double result_array[]);

/* Compute a radial Mathieu function. */
GSL_FUN int gsl_sf_mathieu_Mc_e(int kind, int order, double qq, double zz,
                      gsl_sf_result *result);
GSL_FUN double gsl_sf_mathieu_Mc(int kind, int order, double qq, double zz);
GSL_FUN int gsl_sf_mathieu_Ms_e(int kind, int order, double qq, double zz,
                      gsl_sf_result *result);
GSL_FUN double gsl_sf_mathieu_Ms(int kind, int order, double qq, double zz);
GSL_FUN int gsl_sf_mathieu_Mc_array(int kind, int nmin, int nmax, double qq,
                            double zz, gsl_sf_mathieu_workspace *work,
                            double result_array[]);
GSL_FUN int gsl_sf_mathieu_Ms_array(int kind, int nmin, int nmax, double qq,
                            double zz, gsl_sf_mathieu_workspace *work,
                            double result_array[]);

__END_DECLS

#endif /* !__GSL_SF_MATHIEU_H__ */

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:
#ifndef DISTREADER_H_
#define DISTREADER_H_

#include <ctime>
#include <string>
#include <iostream>
#include <sstream>
#include <cmath>

//#define WINDOWS_32_Q
#define NOGSL

#include <ctime>
#include <string>

#ifdef WINDOWS_32_Q
#include <windows.h>
#define NOGSL
#endif

#ifndef NOGSL
#include <gsl/gsl_cdf.h>
#endif
// convert from integer to string
std::string itos(int j)
{
	std::stringstream s;
	s << j;
	return s.str();
}

// convert double to string (doesn't work though...)
std::string dtos(double j)
{
    std::stringstream s;
    s << j;

    return s.str();
}

unsigned get_nanoseconds()
{
#ifndef WINDOWS_32_Q

	timespec ts;

	if (0!=clock_gettime(CLOCK_REALTIME,&ts))
	{
		throw "error in getting the real time stamp in nanosecs, quitting!";
	}

	return(ts.tv_nsec);

#endif

#ifdef WINDOWS_32_Q
	SYSTEMTIME SystemTime;

	GetSystemTime(&SystemTime);

	return(SystemTime.wMilliseconds);
#endif
}

std::string create_filename(std::string filename)
{
	time_t timep = time(NULL);

	tm *timestruct = localtime(&timep);

	filename.append("_");
	filename.append(itos(timestruct->tm_mday));
	filename.append("_");
	filename.append(itos(timestruct->tm_mon + 1));
	filename.append("_");
	filename.append(itos(timestruct->tm_year + 1900));
	filename.append("_");

	size_t hour = timestruct->tm_hour;

	if (hour < 10)
		filename.append("0");

	filename.append(itos(hour));

	size_t minutes = timestruct->tm_min;

	if (minutes < 10)
		filename.append("0");

	filename.append(itos(minutes));

	size_t seconds = timestruct->tm_sec;

	if (seconds < 10)
		filename.append("0");

	filename.append(itos(seconds));

	filename.append("_");

	// now we are going to insert microsecs
	// to obtain unique filename
	filename.append(itos(get_nanoseconds()));

	return(filename);
}

int linear_search(const int list[], const int maxsize, int value)
{
    int index = 0;
    int position = -1;
    bool found = false;

    while (index < maxsize && !found)
    {
        if (list[index] == value)
        {
            found = true;
            position = index;
        }

        ++index;
    }

    return(position);
}

struct Stats
{
    double mean;
    double mean_ci;
    double sumsquares;
    double var;
    double sumthirds;
    double skew;
    double sumfourths;
    double kurt;
    double sample;
};

struct JointStats
{
    double cov;
    double sum12;
    double sample;
};

void stat_reset(Stats &results)
{
    results.mean = 0;
    results.mean_ci = 0;
    results.sumsquares = 0;
    results.var = 0;
    results.sumthirds = 0;
    results.skew = 0;
    results.sumfourths = 0;
    results.kurt = 0;
    results.sample = 0;
}

void stat_addval(Stats &results, const double traitvalue)
{
    results.mean += traitvalue;
    results.sumsquares += pow(traitvalue,2);
    results.sumthirds += pow(traitvalue,3);
    results.sumfourths += pow(traitvalue,4);
    ++results.sample;
}

void jstat_reset(JointStats &results)
{
    results.sum12 = 0;
    results.sample = 0;
}

void jstat_addval(JointStats &results, const double trait1, const double trait2)
{
    results.sum12+=trait1*trait2;
    ++results.sample;
}

void jstat_finalize(JointStats &results, const double mean1, const double mean2)
{
    results.cov = results.sum12 / results.sample - mean1 * mean2;
}

void stat_finalize(Stats &results)
{
    if (results.sample == 0)
    {
        results.mean = 0;
        results.var = 0;
        results.mean_ci = 0;
        results.skew = 0;
        results.kurt = 0;
    }
    else
    {
        results.mean /= results.sample;
        results.var = results.sumsquares == 0 ? 0 : results.sumsquares / results.sample - pow(results.mean,2);
#ifdef WINDOWS_32_Q
        results.mean_ci = 0;
#endif

#ifndef NOGSL
        results.mean_ci = gsl_cdf_tdist_Qinv(0.025,results.sample - 1) * results.var / sqrt(results.sample);
#endif
        results.skew = results.sumthirds / results.sample - 3 * results.mean * results.sumsquares / results.sample + 2 * pow(results.mean,3);
        results.kurt = (double) (results.sumfourths / results.sample - 4 * results.sumthirds / results.sample * results.mean + 6 * pow(results.mean,2) * results.sumsquares / results.sample - 3 * pow(results.mean,4));

    }
}

#endif