text:
#ifndef EX_UTIL
#define EX_UTIL

#include <gsl/gsl_vector.h>

typedef double (*db_fxy)(double x, double y);
typedef double (*db_fx) (double);

/* print_table_fxy
 * prints f(x,y) tables given range for x and y
 */
void print_table_fxy(db_fxy fun,
        double x_start, double x_end, double x_step, char xch,
        double y_start, double y_end, double y_step, char ych);

/*
 * vec_linspace
 * returns a pointer to a gsl_vector from xlow (inclusive) to xhigh (exclusive)
 * with steps xstep
 */
gsl_vector *vec_linspace(double xlow, double xhigh, double xstep);

/*
 * vec_fx
 * given a function db_fx and vector of xpoints xdata, returns a vector f(x)
 */
gsl_vector *vec_fx(db_fx f, gsl_vector *xdata);

/*
 * vec_from_arr
 * generates a vector from double array
 */
gsl_vector *vec_from_arr(double *arr, int arr_len);

/*
 * sign
 */
int sign(double x);
#endif

ext: h


text:
//////////////////////////////////////////////////////////////////////////////////
// Stefano Allesina Dec 2009
// sallesina@uchicago.edu
// Given a matrix generated by the niche model and a vector containing
// the positions of the species used to generate the matrix, finds the likelihood.
//////////////////////////////////////////////////////////////////////////////////
// Compile
// gcc NewNicheLikelihood.c -o LNICHE -lgsl -lgslcblas -Wall -DHAVE_INLINE -O3
//////////////////////////////////////////////////////////////////////////////////
// Run
// ./LNICHE 4 Matrix4.txt Vec4.txt
// Where:
// * 4 is the number of species
// * Matrix4.txt is a text file containing a space separated adjacency matrix
//               where the rows are the resources and the cols are the consumers
// * Vec4.txt is a space separated text file containing the positions (in [0,1.0])
//            of the species
//////////////////////////////////////////////////////////////////////////////////
// Output
// NONE
//////////////////////////////////////////////////////////////////////////////////

#include<stdio.h>
#include<math.h>
#include <stdlib.h>
#include<gsl/gsl_math.h>
#include<gsl/gsl_integration.h>
#include<gsl/gsl_vector.h>
#include<gsl/gsl_matrix.h>
#include<gsl/gsl_rng.h>
#include<gsl/gsl_randist.h>
#include<gsl/gsl_cdf.h>
#include<gsl/gsl_blas.h>
#include<gsl/gsl_linalg.h>
#include<gsl/gsl_permutation.h>
#include<gsl/gsl_sort_vector_double.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>

// GLOBAL FILE FOR OUTPUT
FILE * FOUT;

// Structure needed for numerical integration
struct myP { double nlm1; double nl; double nr; double nrp1; double n; double beta;};

// UTILITIES
int PrintAdjMat(gsl_matrix * X, int a){
  int i,j;
  for (i=0;i<a;i++){
      for (j=0;j<a;j++){
	fprintf(FOUT,"%1.0f ",gsl_matrix_get(X,i,j));
      }
      fprintf(FOUT,"\n");
  }
  return 0;
}
int PrintVector(gsl_vector * X, int a){
  int i;
  for (i=0;i<a;i++){
    fprintf(FOUT,"%1.3f ",gsl_vector_get(X,i));
  }
  fprintf(FOUT,"\n");
  return 0;
}
// END UTILITIES

// INTEGRATION FOR CONSUMERS
double Inte(double x, void *p){
  struct myP * params  = (struct myP *)p;
  double nlm1=(params->nlm1);
  double nl=(params->nl);
  double nr=(params->nr);
  double nrp1=(params->nrp1);
  double n=(params->n);
  double beta=(params->beta);
  double res=0;
  res=gsl_ran_beta_pdf(x,1.0,beta);
  //fprintf(stderr,"Pdf %f\n",res);
  // Difference
  double c0,c1;
  double r2=(double )x*n/2.;
  if (x*n<=GSL_MIN_DBL(nr-nlm1,nrp1-nl)){
    // case '
    c0=GSL_MIN_DBL(n,nr-r2);
    c1=GSL_MIN_DBL(n,nl+r2);
  }
  else{
    if (x*n<=GSL_MAX_DBL(nr-nlm1,nrp1-nl)){
      // case ''
      if ((nr-nlm1)<=(nrp1-nl)){
	// a<=b
	c0=GSL_MIN_DBL(n,nlm1+r2);
	c1=GSL_MIN_DBL(n,nl+r2);
      }
      else{
	// b<a
	c0=GSL_MIN_DBL(n,nr-r2);
	c1=GSL_MIN_DBL(n,nrp1-r2);
      }
    }
    else{
      // case '''
      c0=GSL_MIN_DBL(n,nlm1+r2);
      c1=GSL_MIN_DBL(n,nrp1-r2);
    }
  }
  res*=(c1-c0);
  //fprintf(stderr,"Center %f (%f)\n",res,(c1-c0));
  // Normalization
  if ((n+r2)<1.0){
    res/=(n-r2);
  }
  else{
    res/=(1.0-2.0*r2);
  }
  //fprintf(stderr,"Rescale %f\n",res);
  return (res);
}

double Integrate (double nlm1, double nl, double n, double nr, double nrp1, double Beta){
  double error;
  double result;
  double low,upp;
  gsl_function F;
  struct myP params = { 0., 0., 0., 0., 0., 0.};
  // setup the environment
  gsl_integration_workspace * w= gsl_integration_workspace_alloc (10000);
  params.beta=Beta;
  params.nlm1=nlm1;
  params.nl=nl;
  params.n=n;
  params.nr=nr;
  params.nrp1=nrp1;

  F.function = &Inte;
  F.params = &params;
  low=GSL_MIN_DBL((nr-nl)/n,1.0);
  upp=GSL_MIN_DBL((nrp1-nlm1)/n,1.0);
  gsl_integration_qags(&F,low,upp,0,1e-4,10000,w,&result, &error);
  //fprintf(stderr,"Prob %f\n",result);
  if (result==0.0) {
    fprintf(FOUT,"The diet is incompatible with the niche model!!\n");
    return -10000.0; // The diet is impossible!
  }
  return log(result);
}

// INTEGRATION FOR BASAL SPECIES
double InteBasal(double x, void *p){
  struct myP * params  = (struct myP *)p;
  double nlm1=(params->nlm1);
  double nrp1=(params->nrp1);
  double n=(params->n);
  double beta=(params->beta);
  double res=0;
  res=gsl_ran_beta_pdf(x,1.0,beta);
  double r2=(double )x*n/2.;
  double c1=nrp1-r2;
  double c0=nlm1+r2;
  res*=(c1-c0);
  //fprintf(stderr,"Center %f (%f)\n",res,(c1-c0));
  // Normalization
  if ((n+r2)<1.0){
    res/=(n-r2);
  }
  else{
    res/=(1.0-2.0*r2);
  }
  return (res);
}

double IntegrateBasal(double GapLeft, double GapRight, double n, double Beta){
  double error;
  double result;
  double low,upp;
  gsl_function F;
  struct myP params = { 0., 0., 0., 0., 0., 0.};
  // setup the environment
  gsl_integration_workspace * w= gsl_integration_workspace_alloc (10000);
  params.beta=Beta;
  params.nlm1=GapLeft;
  params.nl=0.0;
  params.n=n;
  params.nr=0.0;
  params.nrp1=GapRight;

  F.function = &InteBasal;
  F.params = &params;
  low=GSL_MIN_DBL(0.0,1.0);
  upp=GSL_MIN_DBL((GapRight-GapLeft)/n,1.0);
  gsl_integration_qags(&F,low,upp,0,1e-4,10000,w,&result, &error);
  //fprintf(stderr,"Prob %f\n",result);
  if (result==0.0) {
    fprintf(FOUT,"The diet is incompatible with the niche model!!\n");
    return -10000000.0; // The diet is impossible!
  }
  return (result);
}

// COMPUTE THE LIKELIHOOD OF A MATRIX PERFECTLY COMPATIBLE WITH THE NICHE MODEL
double ComputeLikelihood(gsl_matrix *A, gsl_vector *Pos, int N, double Beta){
  double Prob=0.0;
  int i,j;
  double nlm1, nl, n, nr, nrp1;
  double TmpProb;
  for (i=0;i<N;i++){
    if (gsl_matrix_get(A,i,0)==1){
      fprintf(FOUT,"The first species is not basal!!\n");
      Prob-=100000000.0;
    }
  }
  if (Prob>-1.0) fprintf(FOUT,"First    Sp 0, LogLikelihood Diet 0.0 (1.0)\n");
  for (i=1;i<N;i++){
    n=gsl_vector_get(Pos,i);
    nlm1=nl=0.0;
    nr=nrp1=1.0;
    // Find first prey
    for (j=0;j<N;j++){
      if (gsl_matrix_get(A,j,i)==1){
	nl=gsl_vector_get(Pos,j);
	if (j>0){
	  nlm1=gsl_vector_get(Pos,j-1);
	}
	break;
      }
    }
    if (nl>0.0){
      // i is a predator
      // Find the last prey
      for (j=N-1;j>=0;j--){
	if (gsl_matrix_get(A,j,i)==1){
	  nr=gsl_vector_get(Pos,j);
	  if (j<(N-1)){
	    nrp1=gsl_vector_get(Pos,j+1);
	  }
	  break;
	}
      }
      // Now compute the probability
      TmpProb=0.0;
      //fprintf(stderr,"Computing probability of %d (%f %f %f %f %f) Beta: %f\n", i, nlm1, nl, n, nr, nrp1, Beta);
      TmpProb=Integrate(nlm1, nl, n, nr, nrp1, Beta);
      fprintf(FOUT,"Consumer Sp %d, LogLikelihood Diet %f (%f)\n",i,TmpProb,exp(TmpProb));
      Prob+=TmpProb;
    }
    else{
      TmpProb=0.0;
      // for each species compute the probability that the species is not eating any one in between
      double GapLeft, GapRight;
      for (j=0;j<=i;j++){
	if (j==0){
	  GapLeft=0.0;
	}
	else{
	  GapLeft=gsl_vector_get(Pos,j-1);
	}
	GapRight=gsl_vector_get(Pos,j);
	//fprintf(stderr,"Left: %f Right:%f Cumulative: %f\n",GapLeft,GapRight,TmpProb);
	TmpProb+=IntegrateBasal(GapLeft,GapRight, n,Beta);
      }
      fprintf(FOUT,"Basal    Sp %d, LogLikelihood Diet %f (%f)\n",i,log(TmpProb),TmpProb);
      // now that you have summed across all possibilities take the log and sum it
      Prob+=log(TmpProb);
    }
  }
  return Prob;
}

// READ THE MATRIX AND RUN THE PROBABILITY
// Main
int main(int argc, char *argv[]){
  // First argument: number of species
  int Sp=atoi(argv[1]);
  char * FileMat=(argv[2]);
  char * FilePos=(argv[3]);
  // Read the matrix
  FILE * F;
  // matrix
  gsl_matrix * A =gsl_matrix_calloc(Sp,Sp);
  F=fopen(FileMat,"rb");
  gsl_matrix_fscanf(F,A);
  fclose(F);
  int i,j;
  double C,Beta;
  C=Beta=0.0;
  // compute the C and Beta
  for (i=0;i<Sp;i++){
    for (j=0;j<Sp;j++){
      if(gsl_matrix_get(A,i,j)==1) C+=1.0;
    }
  }
  C=(double)C/(double)Sp/(double)Sp;
  Beta=((double) Sp -1.0)/(2.0*C*(double) Sp) -1.0;
  char OutFile[200];
  sprintf(OutFile,"Output-%s",FileMat);
  if ((FOUT = fopen(OutFile, "w")) == NULL){
    puts("Unable to open the file");
    return 0;
  }
  // read the vector
  gsl_vector * Eta=gsl_vector_calloc(Sp);
  F=fopen(FilePos,"rb");
  gsl_vector_fscanf(F,Eta);
  fclose(F);
  fprintf(FOUT,"Compatible Adjacency Matrix\n");
  PrintAdjMat(A,Sp);
  fprintf(FOUT,"Species Positions\n");
  PrintVector(Eta,Sp);
  fprintf(FOUT,"Beta: %f\n",Beta);
  fprintf(FOUT,"LogLikelihood: %f\n",ComputeLikelihood(A,Eta,Sp,Beta));
  gsl_matrix_free(A);
  gsl_vector_free(Eta);
  fclose(FOUT);
  return 0;
}

ext: c


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

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

#include "ccl.h"

typedef struct{
  double l;
  ccl_cosmology *cosmo;
  ccl_cl_tracer_collection_t *trc1;
  ccl_cl_tracer_collection_t *trc2;
  ccl_f2d_t *psp;
  int *status;
} integ_cl_par;

typedef struct{
  int chipow;
  double l1;
  double l2;
  ccl_cosmology *cosmo;
  ccl_cl_tracer_collection_t *trc1;
  ccl_cl_tracer_collection_t *trc2;
  ccl_cl_tracer_collection_t *trc3;
  ccl_cl_tracer_collection_t *trc4;
  ccl_f3d_t *tsp;
  ccl_f1d_t *ker_extra;
  ccl_a_finder *finda;
  int *status;
} integ_cov_par;

static void update_chi_limits(ccl_cl_tracer_collection_t *trc,
                              double *chimin, double *chimax,
                              int is_union)
{
  int itr;
  double chimin_h=1E15;
  double chimax_h=-1E15;
  for(itr=0; itr < trc->n_tracers; itr++) {
    if (trc->ts[itr]->chi_min < chimin_h)
      chimin_h = trc->ts[itr]->chi_min;
    if (trc->ts[itr]->chi_max > chimax_h)
      chimax_h = trc->ts[itr]->chi_max;
  }

  if(is_union) {
    if(chimin_h < *chimin)
      *chimin = chimin_h;
    if(chimax_h > *chimax)
      *chimax = chimax_h;
  }
  else {
    if(chimin_h > *chimin)
      *chimin = chimin_h;
    if(chimax_h < *chimax)
      *chimax = chimax_h;
  }
}

static void get_k_interval(ccl_cosmology *cosmo,
                           ccl_cl_tracer_collection_t *trc1,
                           ccl_cl_tracer_collection_t *trc2,
                           double l, double *lkmin, double *lkmax) {
  int itr;

  // Loop through all tracers and find distance bounds
  double chi_min1 = 1E15;
  double chi_max1 = -1E15;
  for (itr=0; itr < trc1->n_tracers; itr++) {
    if (trc1->ts[itr]->chi_min < chi_min1)
      chi_min1 = trc1->ts[itr]->chi_min;
    if (trc1->ts[itr]->chi_max > chi_max1)
      chi_max1 = trc1->ts[itr]->chi_max;
  }

  double chi_min2 = 1E15;
  double chi_max2 = -1E15;
  for (itr=0; itr < trc2->n_tracers; itr++) {
    if (trc2->ts[itr]->chi_min < chi_min2)
      chi_min2 = trc2->ts[itr]->chi_min;
    if (trc2->ts[itr]->chi_max > chi_max2)
      chi_max2 = trc2->ts[itr]->chi_max;
  }

  // Find maximum of minima and minimum of maxima
  // (i.e. edges where the product of both kernels will have support).
  double chi_min = fmax(chi_min1, chi_min2);
  double chi_max = fmin(chi_max1, chi_max2);

  if (chi_min <= 0)
    chi_min = 0.5*(l+0.5)/cosmo->spline_params.K_MAX;

  // Don't go beyond kmax
  *lkmax = log(fmin(cosmo->spline_params.K_MAX, 2*(l+0.5)/chi_min));
  *lkmin = log(fmax(cosmo->spline_params.K_MIN, (l+0.5)/chi_max));
}

static double transfer_limber_single(ccl_cl_tracer_t *tr, double l, double lk,
                                     double k, double chi_l, double a_l,
                                     ccl_cosmology *cosmo, ccl_f2d_t *psp,
                                     int ignore_jbes_deriv,
                                     int *status) {
  double dd = 0;

  // Kernel and transfer evaluated at chi_l
  double w = ccl_cl_tracer_t_get_kernel(tr, chi_l, status);
  double t = ccl_cl_tracer_t_get_transfer(tr, lk,a_l, status);
  double fl = ccl_cl_tracer_t_get_f_ell(tr, l, status);

  if (tr->der_bessel < 1) { //We don't need l+1
    dd = w*t;
    if (tr->der_bessel == -1) { //If we divide by (chi*k)^2
      double lp1h = l+0.5;
      dd /= (lp1h*lp1h);
    }
  }
  else { // We will need l+1
    if(ignore_jbes_deriv)
      dd = 0;
    else {
      // Compute chi_{l+1} and a_{l+1}
      double lp1h = l+0.5;
      double lp3h = l+1.5;
      double chi_lp = lp3h/k;
      double a_lp = ccl_scale_factor_of_chi(cosmo, chi_lp, status);

      // Compute power spectrum ratio there
      double pk_ratio = fabs(ccl_f2d_t_eval(psp, lk, a_lp, cosmo, status) /
                             ccl_f2d_t_eval(psp, lk, a_l, cosmo, status));

      // Compute kernel and trasfer at chi_{l+1}
      double w_p = ccl_cl_tracer_t_get_kernel(tr, chi_lp, status);
      double t_p = ccl_cl_tracer_t_get_transfer(tr, lk,a_lp, status);

      // sqrt(2l+1/2l+3)
      double sqell = sqrt(lp1h*pk_ratio/lp3h);
      if (tr->der_bessel == 1)
        dd = l*w*t/lp1h-sqell*w_p*t_p;
      else //we assume der_bessel=2 here to avoid extra if clause
        dd = sqell*2*w_p*t_p/lp3h - (0.25+2*l)*w*t/(lp1h*lp1h);
    }
  }
  return dd*fl;
}

static double transfer_limber_wrap(double l,double lk, double k, double chi,
                                   double a, ccl_cl_tracer_collection_t *trc,
                                   ccl_cosmology *cosmo,ccl_f2d_t *psp,
                                   int ignore_jbes_deriv, int *status) {
  int itr;
  double transfer = 0;

  for (itr=0; itr < trc->n_tracers; itr++) {
    transfer += transfer_limber_single(
      trc->ts[itr], l, lk, k, chi, a, cosmo, psp, ignore_jbes_deriv, status);
    if (*status != 0)
      return -1;
  }
  return transfer;
}

static double cl_integrand(double lk, void *params) {
  double d1, d2;
  integ_cl_par *p = (integ_cl_par *)params;
  double k = exp(lk);
  double chi = (p->l+0.5)/k;
  double a = ccl_scale_factor_of_chi(p->cosmo, chi, p->status);

  d1 = transfer_limber_wrap(p->l, lk, k, chi, a, p->trc1,
                            p->cosmo, p->psp, 0, p->status);
  if (d1 == 0)
    return 0;

  d2 = transfer_limber_wrap(p->l, lk, k, chi, a, p->trc2,
                            p->cosmo, p->psp, 0, p->status);

  if (d2 == 0)
    return 0;

  double pk = ccl_f2d_t_eval(p->psp, lk, a, p->cosmo, p->status);

  return k*pk*d1*d2;
}

static void integ_cls_limber_spline(ccl_cosmology *cosmo,
				    integ_cl_par *ipar,
				    double lkmin, double lkmax,
				    double *result, int *status) {
  int ik;
  int nk = (int)(fmax((lkmax - lkmin) / cosmo->spline_params.DLOGK_INTEGRATION + 0.5,
		      1))+1;
  double *fk_arr = NULL;
  double *lk_arr = NULL;
  lk_arr = ccl_linear_spacing(lkmin, lkmax, nk);
  if(lk_arr == NULL)
    *status = CCL_ERROR_LOGSPACE;

  if(*status == 0) {
    fk_arr = malloc(nk * sizeof(double));
    if(fk_arr == NULL)
      *status = CCL_ERROR_MEMORY;
  }

  if(*status == 0) {
    for(ik=0; ik<nk; ik++) {
      fk_arr[ik] = cl_integrand(lk_arr[ik], ipar);
      if(*(ipar->status)) {
	*status = *(ipar->status);
	break;
      }
    }
  }

  if(*status == 0) {
    ccl_integ_spline(1, nk, lk_arr, &fk_arr,
                     1, -1, result, gsl_interp_akima,
                     status);
  }
  free(fk_arr);
  free(lk_arr);
}

static void integ_cls_limber_qag_quad(ccl_cosmology *cosmo,
				      gsl_function *F,
				      double lkmin, double lkmax,
				      gsl_integration_workspace *w,
				      double *result, double *eresult,
				      int *status) {
  int gslstatus;
  size_t nevals;
  gsl_integration_cquad_workspace *w_cquad = NULL;
  // Integrate
  gslstatus = gsl_integration_qag(F, lkmin, lkmax, 0,
				  cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL,
				  cosmo->gsl_params.N_ITERATION,
				  cosmo->gsl_params.INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS,
				  w, result, eresult);

  // Test if a round-off error occured in the evaluation of the integral
  // If so, try another integration function, more robust but potentially slower
  if (gslstatus == GSL_EROUND) {
    ccl_raise_gsl_warning(gslstatus,
			  "ccl_cls.c: integ_cls_limber_qag_quad(): "
			  "Default GSL integration failure, attempting backup method.");
    w_cquad = gsl_integration_cquad_workspace_alloc(cosmo->gsl_params.N_ITERATION);
    if (w_cquad == NULL)
      *status = CCL_ERROR_MEMORY;

    if (*status == 0) {
      nevals = 0;
      gslstatus = gsl_integration_cquad(F, lkmin, lkmax, 0,
					cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL,
					w_cquad, result, eresult, &nevals);
    }
  }
  gsl_integration_cquad_workspace_free(w_cquad);
  if(*status == 0)
    *status = gslstatus;
}

void ccl_angular_cls_limber(ccl_cosmology *cosmo,
			    ccl_cl_tracer_collection_t *trc1,
			    ccl_cl_tracer_collection_t *trc2,
			    ccl_f2d_t *psp,
			    int nl_out, double *l_out, double *cl_out,
			    ccl_integration_t integration_method,
			    int *status) {

  // make sure to init core things for safety
  if (!cosmo->computed_distances) {
    *status = CCL_ERROR_DISTANCES_INIT;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_cls.c: ccl_angular_cls_limber(): distance splines have not been precomputed!");
    return;
  }

  #pragma omp parallel shared(cosmo, trc1, trc2, l_out, cl_out, \
                              nl_out, status, psp, integration_method) \
                       default(none)
  {
    int clastatus, lind;
    integ_cl_par ipar;
    gsl_integration_workspace *w = NULL;
    int local_status = *status;
    gsl_function F;
    double lkmin, lkmax, l, result, eresult;

    if (local_status == 0) {
      // Set up integrating function parameters
      ipar.cosmo = cosmo;
      ipar.trc1 = trc1;
      ipar.trc2 = trc2;
      ipar.psp = psp;
      ipar.status = &clastatus;
    }

    if(integration_method == ccl_integration_qag_quad) {
      if (local_status == 0) {
	w = gsl_integration_workspace_alloc(cosmo->gsl_params.N_ITERATION);
	if (w == NULL) {
	  local_status = CCL_ERROR_MEMORY;
	}
      }

      if (local_status == 0) {
	// Set up integrating function
	F.function = &cl_integrand;
	F.params = &ipar;
      }
    }

    #pragma omp for schedule(dynamic)
    for (lind=0; lind < nl_out; ++lind) {
      if (local_status == 0) {
        l = l_out[lind];
        clastatus = 0;
        ipar.l = l;

        // Get integration limits
        get_k_interval(cosmo, trc1, trc2, l, &lkmin, &lkmax);

	// Integrate
	if(integration_method == ccl_integration_qag_quad) {
	  integ_cls_limber_qag_quad(cosmo, &F, lkmin, lkmax, w,
				    &result, &eresult, &local_status);
	}
	else if(integration_method == ccl_integration_spline) {
	  integ_cls_limber_spline(cosmo, &ipar, lkmin, lkmax,
				  &result, &local_status);
	}
	else
	  local_status = CCL_ERROR_NOT_IMPLEMENTED;

        if ((*ipar.status == 0) && (local_status == 0)) {
          cl_out[lind] = result / (l+0.5);
        }
        else {
          ccl_raise_gsl_warning(local_status, "ccl_cls.c: ccl_angular_cls_limber():");
          cl_out[lind] = NAN;
          local_status = CCL_ERROR_INTEG;
        }
      }
    }

    gsl_integration_workspace_free(w);

    if (local_status) {
      #pragma omp atomic write
      *status = local_status;
    }
  }

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

void ccl_angular_cls_nonlimber(ccl_cosmology *cosmo,
                               ccl_cl_tracer_collection_t *trc1,
                               ccl_cl_tracer_collection_t *trc2,
                               ccl_f2d_t *psp,
                               int nl_out, int *l_out, double *cl_out,
                               int *status) {
  *status = CCL_ERROR_INCONSISTENT;
  ccl_cosmology_set_status_message(
    cosmo,
    "ccl_cls.c: ccl_angular_cls_nonlimber(); non-Limber integrator not implemented yet\n");
}

static double cov_integrand(double chi, void *params)
{
  double d1, d2, d3, d4, tkk, ker=1;
  integ_cov_par *p = (integ_cov_par *)params;
  double k1=(p->l1+0.5)/chi;
  double k2=(p->l2+0.5)/chi;
  double lk1=log(k1);
  double lk2=log(k2);
  double a = ccl_scale_factor_of_chi(p->cosmo, chi, p->status);

  d1 = transfer_limber_wrap(p->l1, lk1, k1, chi, a, p->trc1,
                            p->cosmo, NULL, 1, p->status);
  if (d1 == 0)
    return 0;
  d2 = transfer_limber_wrap(p->l1, lk1, k1, chi, a, p->trc2,
                            p->cosmo, NULL, 1, p->status);
  if (d2 == 0)
    return 0;
  d3 = transfer_limber_wrap(p->l2, lk2, k2, chi, a, p->trc3,
                            p->cosmo, NULL, 1, p->status);
  if (d3 == 0)
    return 0;
  d4 = transfer_limber_wrap(p->l2, lk2, k2, chi, a, p->trc4,
                            p->cosmo, NULL, 1, p->status);
  if (d4 == 0)
    return 0;

  tkk = ccl_f3d_t_eval(p->tsp, lk1, lk2, a,
                       p->finda, p->cosmo, p->status);
  if(p->ker_extra!=NULL)
    ker = ccl_f1d_t_eval(p->ker_extra, chi);

  return d1*d2*d3*d4*tkk*ker/pow(chi, p->chipow);
}

static void integ_cov_limber_spline(ccl_cosmology *cosmo,
				    integ_cov_par *ipar,
				    double chimin, double chimax,
				    double *result, int *status)
{
  int ichi;
  int nchi = (int)(fmax((chimax - chimin) / cosmo->spline_params.DCHI_INTEGRATION + 0.5,
		      1))+1;
  double *fchi_arr = NULL;
  double *chi_arr = NULL;
  chi_arr = ccl_linear_spacing(chimin, chimax, nchi);
  if(chi_arr == NULL)
    *status = CCL_ERROR_LOGSPACE;

  if(*status == 0) {
    fchi_arr = malloc(nchi * sizeof(double));
    if(fchi_arr == NULL)
      *status = CCL_ERROR_MEMORY;
  }

  if(*status == 0) {
    for(ichi=0; ichi<nchi; ichi++) {
      fchi_arr[ichi] = cov_integrand(chi_arr[ichi], ipar);
      if(*(ipar->status)) {
	*status = *(ipar->status);
	break;
      }
    }
  }

  if(*status == 0) {
    ccl_integ_spline(1, nchi, chi_arr, &fchi_arr,
                     1, -1, result, gsl_interp_akima,
                     status);
  }

  free(fchi_arr);
  free(chi_arr);
}

static void integ_cov_limber_qag_quad(ccl_cosmology *cosmo,
				      gsl_function *F,
				      double chimin, double chimax,
				      gsl_integration_workspace *w,
				      double *result, double *eresult,
				      int *status) {
  int gslstatus;
  size_t nevals;
  gsl_integration_cquad_workspace *w_cquad = NULL;
  // Integrate
  gslstatus = gsl_integration_qag(F, chimin, chimax, 0,
				  cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL,
				  cosmo->gsl_params.N_ITERATION,
				  cosmo->gsl_params.INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS,
				  w, result, eresult);

  // Test if a round-off error occured in the evaluation of the integral
  // If so, try another integration function, more robust but potentially slower
  if (gslstatus == GSL_EROUND) {
    ccl_raise_gsl_warning(gslstatus,
			  "ccl_cls.c: ccl_angular_cov_limber(): "
			  "Default GSL integration failure, attempting backup method.");
    w_cquad = gsl_integration_cquad_workspace_alloc(cosmo->gsl_params.N_ITERATION);
    if (w_cquad == NULL)
      *status = CCL_ERROR_MEMORY;

    if (*status == 0) {
      nevals = 0;
      gslstatus = gsl_integration_cquad(F, chimin, chimax, 0,
					cosmo->gsl_params.INTEGRATION_LIMBER_EPSREL,
					w_cquad, result, eresult, &nevals);
    }
  }
  gsl_integration_cquad_workspace_free(w_cquad);
  if(*status == 0)
    *status = gslstatus;
}

void ccl_angular_cl_covariance(ccl_cosmology *cosmo,
                               ccl_cl_tracer_collection_t *trc1,
                               ccl_cl_tracer_collection_t *trc2,
                               ccl_cl_tracer_collection_t *trc3,
                               ccl_cl_tracer_collection_t *trc4,
                               ccl_f3d_t *tsp,
                               int nl1_out, double *l1_out,
                               int nl2_out, double *l2_out,
                               double *cov_out,
                               ccl_integration_t integration_method,
                               int chi_exponent, ccl_f1d_t *kernel_extra,
                               double prefactor_extra, int *status)
{
  if(!cosmo->computed_distances) {
    *status = CCL_ERROR_DISTANCES_INIT;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_cls.c: ccl_angular_cl_limber(): distance splines have not been precomputed!");
    return;
  }

  #pragma omp parallel shared(cosmo, trc1, trc2, trc3, trc4, tsp, \
                              nl1_out, l1_out, nl2_out, l2_out, cov_out, \
                              integration_method, chi_exponent, \
                              kernel_extra, prefactor_extra, status) \
                       default(none)
  {
    int clastatus, lind1,lind2;
    integ_cov_par ipar;
    gsl_integration_workspace *w = NULL;
    int local_status = *status;
    gsl_function F;
    double chimin, chimax;
    double l1, l2, result, eresult;
    ccl_a_finder *finda = ccl_a_finder_new_from_f3d(tsp);

    // Find integration limits
    chimin = 1E15;
    chimax = -1E15;
    update_chi_limits(trc1, &chimin, &chimax, 1);
    update_chi_limits(trc2, &chimin, &chimax, 0);
    update_chi_limits(trc3, &chimin, &chimax, 0);
    update_chi_limits(trc4, &chimin, &chimax, 0);

    if (local_status == 0) {
      // Set up integrating function parameters
      ipar.cosmo = cosmo;
      ipar.trc1 = trc1;
      ipar.trc2 = trc2;
      ipar.trc3 = trc3;
      ipar.trc4 = trc4;
      ipar.tsp = tsp;
      ipar.ker_extra = kernel_extra;
      ipar.finda = finda;
      ipar.status = &clastatus;
      ipar.chipow = chi_exponent;
    }

    if(integration_method == ccl_integration_qag_quad) {
      if (local_status == 0) {
	w = gsl_integration_workspace_alloc(cosmo->gsl_params.N_ITERATION);
	if (w == NULL) {
	  local_status = CCL_ERROR_MEMORY;
	}
      }

      if (local_status == 0) {
	// Set up integrating function
	F.function = &cov_integrand;
	F.params = &ipar;
      }
    }

    #pragma omp for schedule(dynamic)
    for (lind1=0; lind1 < nl1_out; ++lind1) {
      l1 = l1_out[lind1];
      ipar.l1 = l1;

      for (lind2=0; lind2 < nl2_out; ++lind2) {
        if (local_status == 0) {
          l2 = l2_out[lind2];
          clastatus = 0;
          ipar.l2 = l2;

          // Integrate
          if(integration_method == ccl_integration_qag_quad) {
            integ_cov_limber_qag_quad(cosmo, &F, chimin, chimax, w,
                                      &result, &eresult, &local_status);
          }
          else if(integration_method == ccl_integration_spline) {
            integ_cov_limber_spline(cosmo, &ipar, chimin, chimax,
                                    &result, &local_status);
          }
          else
            local_status = CCL_ERROR_NOT_IMPLEMENTED;

          if ((*ipar.status == 0) && (local_status == 0)) {
            cov_out[lind1+nl1_out*lind2] = result * prefactor_extra;
          }
          else {
            ccl_raise_gsl_warning(local_status, "ccl_cls.c: ccl_angular_cov_limber():");
            cov_out[lind1+nl1_out*lind2] = NAN;
            local_status = CCL_ERROR_INTEG;
          }
        }
      }
    }

    gsl_integration_workspace_free(w);

    if (local_status) {
      #pragma omp atomic write
        *status = local_status;
    }
    ccl_a_finder_free(finda);
  }

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

max_stars_repo_name: carlosggarcia/CCL

max_issues_repo_name: carlosggarcia/CCL

max_forks_repo_name: carlosggarcia/CCL


text:
#ifndef SQ_INCLUDE_GUARD_core_SharedRange_h_
#define SQ_INCLUDE_GUARD_core_SharedRange_h_

#include <gsl/gsl>
#include <memory>
#include <range/v3/range/concepts.hpp>
#include <vector>

namespace sq {

template <typename R> requires ranges::cpp20::range<R> struct SharedRange {
public:
  explicit SharedRange(const std::shared_ptr<R> &base) : base_{base} {}

  explicit SharedRange(std::shared_ptr<R> &&base) : base_{std::move(base)} {}

  constexpr SharedRange() noexcept = default;
  SharedRange(const SharedRange &) noexcept = default;
  SharedRange(SharedRange &&) noexcept = default;
  SharedRange &operator=(const SharedRange &) noexcept = default;
  SharedRange &operator=(SharedRange &&) noexcept = default;
  ~SharedRange() noexcept = default;

  auto begin() {
    Expects(base_ != nullptr);
    return ranges::begin(*base_);
  }
  auto end() {
    Expects(base_ != nullptr);
    return ranges::end(*base_);
  }

private:
  // views must be semiregular, which means that SharedRange must have a
  // default constructor, which, unfortunately, means that we can't use
  // gsl::not_null here.
  std::shared_ptr<R> base_ = nullptr;
};

} // namespace sq

namespace ranges {

template <typename R>
inline constexpr bool enable_view<sq::SharedRange<R>> = true;

} // namespace ranges

namespace sq {

static_assert(ranges::cpp20::range<SharedRange<std::vector<int>>>);
static_assert(ranges::cpp20::view<SharedRange<std::vector<int>>>);
static_assert(ranges::viewable_range<SharedRange<std::vector<int>>>);

} // namespace sq

#endif // SQ_INCLUDE_GUARD_core_SharedRange_h_

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Mon Mar  8 17:45:11 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -notwiddleinv 32 */

/*
 * This function contains 372 FP additions, 84 FP multiplications,
 * (or, 340 additions, 52 multiplications, 32 fused multiply/add),
 * 92 stack variables, and 128 memory accesses
 */
static const fftw_real K195090322 = FFTW_KONST(+0.195090322016128267848284868477022240927691618);
static const fftw_real K980785280 = FFTW_KONST(+0.980785280403230449126182236134239036973933731);
static const fftw_real K831469612 = FFTW_KONST(+0.831469612302545237078788377617905756738560812);
static const fftw_real K555570233 = FFTW_KONST(+0.555570233019602224742830813948532874374937191);
static const fftw_real K382683432 = FFTW_KONST(+0.382683432365089771728459984030398866761344562);
static const fftw_real K923879532 = FFTW_KONST(+0.923879532511286756128183189396788286822416626);
static const fftw_real K707106781 = FFTW_KONST(+0.707106781186547524400844362104849039284835938);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.36 1999/02/19 17:22:11 athena Exp $
 * $Id: fft.ml,v 1.41 1999/02/19 17:22:13 athena Exp $
 * $Id: to_c.ml,v 1.24 1999/02/19 17:22:17 athena Exp $
 */

void fftwi_no_twiddle_32(const fftw_complex *input, fftw_complex *output, int istride, int ostride)
{
     fftw_real tmp7;
     fftw_real tmp339;
     fftw_real tmp70;
     fftw_real tmp313;
     fftw_real tmp97;
     fftw_real tmp215;
     fftw_real tmp179;
     fftw_real tmp241;
     fftw_real tmp14;
     fftw_real tmp314;
     fftw_real tmp77;
     fftw_real tmp340;
     fftw_real tmp182;
     fftw_real tmp216;
     fftw_real tmp104;
     fftw_real tmp242;
     fftw_real tmp153;
     fftw_real tmp236;
     fftw_real tmp53;
     fftw_real tmp60;
     fftw_real tmp287;
     fftw_real tmp336;
     fftw_real tmp360;
     fftw_real tmp290;
     fftw_real tmp293;
     fftw_real tmp294;
     fftw_real tmp170;
     fftw_real tmp233;
     fftw_real tmp333;
     fftw_real tmp359;
     fftw_real tmp164;
     fftw_real tmp234;
     fftw_real tmp173;
     fftw_real tmp237;
     fftw_real tmp22;
     fftw_real tmp318;
     fftw_real tmp343;
     fftw_real tmp85;
     fftw_real tmp112;
     fftw_real tmp185;
     fftw_real tmp220;
     fftw_real tmp245;
     fftw_real tmp29;
     fftw_real tmp321;
     fftw_real tmp342;
     fftw_real tmp92;
     fftw_real tmp119;
     fftw_real tmp184;
     fftw_real tmp223;
     fftw_real tmp244;
     fftw_real tmp126;
     fftw_real tmp229;
     fftw_real tmp38;
     fftw_real tmp45;
     fftw_real tmp278;
     fftw_real tmp329;
     fftw_real tmp357;
     fftw_real tmp281;
     fftw_real tmp284;
     fftw_real tmp285;
     fftw_real tmp143;
     fftw_real tmp226;
     fftw_real tmp326;
     fftw_real tmp356;
     fftw_real tmp137;
     fftw_real tmp227;
     fftw_real tmp146;
     fftw_real tmp230;
     ASSERT_ALIGNED_DOUBLE();
     {
	  fftw_real tmp3;
	  fftw_real tmp177;
	  fftw_real tmp66;
	  fftw_real tmp96;
	  fftw_real tmp6;
	  fftw_real tmp95;
	  fftw_real tmp69;
	  fftw_real tmp178;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp1;
	       fftw_real tmp2;
	       fftw_real tmp64;
	       fftw_real tmp65;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp1 = c_re(input[0]);
	       tmp2 = c_re(input[16 * istride]);
	       tmp3 = tmp1 + tmp2;
	       tmp177 = tmp1 - tmp2;
	       tmp64 = c_im(input[0]);
	       tmp65 = c_im(input[16 * istride]);
	       tmp66 = tmp64 + tmp65;
	       tmp96 = tmp64 - tmp65;
	  }
	  {
	       fftw_real tmp4;
	       fftw_real tmp5;
	       fftw_real tmp67;
	       fftw_real tmp68;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp4 = c_re(input[8 * istride]);
	       tmp5 = c_re(input[24 * istride]);
	       tmp6 = tmp4 + tmp5;
	       tmp95 = tmp4 - tmp5;
	       tmp67 = c_im(input[8 * istride]);
	       tmp68 = c_im(input[24 * istride]);
	       tmp69 = tmp67 + tmp68;
	       tmp178 = tmp67 - tmp68;
	  }
	  tmp7 = tmp3 + tmp6;
	  tmp339 = tmp3 - tmp6;
	  tmp70 = tmp66 + tmp69;
	  tmp313 = tmp66 - tmp69;
	  tmp97 = tmp95 + tmp96;
	  tmp215 = tmp96 - tmp95;
	  tmp179 = tmp177 - tmp178;
	  tmp241 = tmp177 + tmp178;
     }
     {
	  fftw_real tmp10;
	  fftw_real tmp98;
	  fftw_real tmp73;
	  fftw_real tmp99;
	  fftw_real tmp13;
	  fftw_real tmp102;
	  fftw_real tmp76;
	  fftw_real tmp101;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp8;
	       fftw_real tmp9;
	       fftw_real tmp71;
	       fftw_real tmp72;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp8 = c_re(input[4 * istride]);
	       tmp9 = c_re(input[20 * istride]);
	       tmp10 = tmp8 + tmp9;
	       tmp98 = tmp8 - tmp9;
	       tmp71 = c_im(input[4 * istride]);
	       tmp72 = c_im(input[20 * istride]);
	       tmp73 = tmp71 + tmp72;
	       tmp99 = tmp71 - tmp72;
	  }
	  {
	       fftw_real tmp11;
	       fftw_real tmp12;
	       fftw_real tmp74;
	       fftw_real tmp75;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp11 = c_re(input[28 * istride]);
	       tmp12 = c_re(input[12 * istride]);
	       tmp13 = tmp11 + tmp12;
	       tmp102 = tmp11 - tmp12;
	       tmp74 = c_im(input[28 * istride]);
	       tmp75 = c_im(input[12 * istride]);
	       tmp76 = tmp74 + tmp75;
	       tmp101 = tmp74 - tmp75;
	  }
	  tmp14 = tmp10 + tmp13;
	  tmp314 = tmp10 - tmp13;
	  tmp77 = tmp73 + tmp76;
	  tmp340 = tmp76 - tmp73;
	  {
	       fftw_real tmp180;
	       fftw_real tmp181;
	       fftw_real tmp100;
	       fftw_real tmp103;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp180 = tmp98 - tmp99;
	       tmp181 = tmp102 + tmp101;
	       tmp182 = K707106781 * (tmp180 + tmp181);
	       tmp216 = K707106781 * (tmp180 - tmp181);
	       tmp100 = tmp98 + tmp99;
	       tmp103 = tmp101 - tmp102;
	       tmp104 = K707106781 * (tmp100 + tmp103);
	       tmp242 = K707106781 * (tmp103 - tmp100);
	  }
     }
     {
	  fftw_real tmp49;
	  fftw_real tmp149;
	  fftw_real tmp169;
	  fftw_real tmp288;
	  fftw_real tmp52;
	  fftw_real tmp166;
	  fftw_real tmp152;
	  fftw_real tmp289;
	  fftw_real tmp56;
	  fftw_real tmp154;
	  fftw_real tmp157;
	  fftw_real tmp291;
	  fftw_real tmp59;
	  fftw_real tmp159;
	  fftw_real tmp162;
	  fftw_real tmp292;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp47;
	       fftw_real tmp48;
	       fftw_real tmp167;
	       fftw_real tmp168;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp47 = c_re(input[31 * istride]);
	       tmp48 = c_re(input[15 * istride]);
	       tmp49 = tmp47 + tmp48;
	       tmp149 = tmp47 - tmp48;
	       tmp167 = c_im(input[31 * istride]);
	       tmp168 = c_im(input[15 * istride]);
	       tmp169 = tmp167 - tmp168;
	       tmp288 = tmp167 + tmp168;
	  }
	  {
	       fftw_real tmp50;
	       fftw_real tmp51;
	       fftw_real tmp150;
	       fftw_real tmp151;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp50 = c_re(input[7 * istride]);
	       tmp51 = c_re(input[23 * istride]);
	       tmp52 = tmp50 + tmp51;
	       tmp166 = tmp50 - tmp51;
	       tmp150 = c_im(input[7 * istride]);
	       tmp151 = c_im(input[23 * istride]);
	       tmp152 = tmp150 - tmp151;
	       tmp289 = tmp150 + tmp151;
	  }
	  {
	       fftw_real tmp54;
	       fftw_real tmp55;
	       fftw_real tmp155;
	       fftw_real tmp156;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp54 = c_re(input[3 * istride]);
	       tmp55 = c_re(input[19 * istride]);
	       tmp56 = tmp54 + tmp55;
	       tmp154 = tmp54 - tmp55;
	       tmp155 = c_im(input[3 * istride]);
	       tmp156 = c_im(input[19 * istride]);
	       tmp157 = tmp155 - tmp156;
	       tmp291 = tmp155 + tmp156;
	  }
	  {
	       fftw_real tmp57;
	       fftw_real tmp58;
	       fftw_real tmp160;
	       fftw_real tmp161;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp57 = c_re(input[27 * istride]);
	       tmp58 = c_re(input[11 * istride]);
	       tmp59 = tmp57 + tmp58;
	       tmp159 = tmp57 - tmp58;
	       tmp160 = c_im(input[27 * istride]);
	       tmp161 = c_im(input[11 * istride]);
	       tmp162 = tmp160 - tmp161;
	       tmp292 = tmp160 + tmp161;
	  }
	  {
	       fftw_real tmp334;
	       fftw_real tmp335;
	       fftw_real tmp331;
	       fftw_real tmp332;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp153 = tmp149 - tmp152;
	       tmp236 = tmp149 + tmp152;
	       tmp53 = tmp49 + tmp52;
	       tmp60 = tmp56 + tmp59;
	       tmp287 = tmp53 - tmp60;
	       tmp334 = tmp49 - tmp52;
	       tmp335 = tmp292 - tmp291;
	       tmp336 = tmp334 - tmp335;
	       tmp360 = tmp334 + tmp335;
	       tmp290 = tmp288 + tmp289;
	       tmp293 = tmp291 + tmp292;
	       tmp294 = tmp290 - tmp293;
	       tmp170 = tmp166 + tmp169;
	       tmp233 = tmp169 - tmp166;
	       tmp331 = tmp288 - tmp289;
	       tmp332 = tmp56 - tmp59;
	       tmp333 = tmp331 - tmp332;
	       tmp359 = tmp332 + tmp331;
	       {
		    fftw_real tmp158;
		    fftw_real tmp163;
		    fftw_real tmp171;
		    fftw_real tmp172;
		    ASSERT_ALIGNED_DOUBLE();
		    tmp158 = tmp154 - tmp157;
		    tmp163 = tmp159 + tmp162;
		    tmp164 = K707106781 * (tmp158 + tmp163);
		    tmp234 = K707106781 * (tmp158 - tmp163);
		    tmp171 = tmp154 + tmp157;
		    tmp172 = tmp162 - tmp159;
		    tmp173 = K707106781 * (tmp171 + tmp172);
		    tmp237 = K707106781 * (tmp172 - tmp171);
	       }
	  }
     }
     {
	  fftw_real tmp18;
	  fftw_real tmp106;
	  fftw_real tmp81;
	  fftw_real tmp110;
	  fftw_real tmp21;
	  fftw_real tmp109;
	  fftw_real tmp84;
	  fftw_real tmp107;
	  fftw_real tmp316;
	  fftw_real tmp317;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp16;
	       fftw_real tmp17;
	       fftw_real tmp79;
	       fftw_real tmp80;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp16 = c_re(input[2 * istride]);
	       tmp17 = c_re(input[18 * istride]);
	       tmp18 = tmp16 + tmp17;
	       tmp106 = tmp16 - tmp17;
	       tmp79 = c_im(input[2 * istride]);
	       tmp80 = c_im(input[18 * istride]);
	       tmp81 = tmp79 + tmp80;
	       tmp110 = tmp79 - tmp80;
	  }
	  {
	       fftw_real tmp19;
	       fftw_real tmp20;
	       fftw_real tmp82;
	       fftw_real tmp83;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp19 = c_re(input[10 * istride]);
	       tmp20 = c_re(input[26 * istride]);
	       tmp21 = tmp19 + tmp20;
	       tmp109 = tmp19 - tmp20;
	       tmp82 = c_im(input[10 * istride]);
	       tmp83 = c_im(input[26 * istride]);
	       tmp84 = tmp82 + tmp83;
	       tmp107 = tmp82 - tmp83;
	  }
	  tmp22 = tmp18 + tmp21;
	  tmp316 = tmp18 - tmp21;
	  tmp317 = tmp81 - tmp84;
	  tmp318 = tmp316 - tmp317;
	  tmp343 = tmp316 + tmp317;
	  tmp85 = tmp81 + tmp84;
	  {
	       fftw_real tmp108;
	       fftw_real tmp111;
	       fftw_real tmp218;
	       fftw_real tmp219;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp108 = tmp106 - tmp107;
	       tmp111 = tmp109 + tmp110;
	       tmp112 = (K923879532 * tmp108) - (K382683432 * tmp111);
	       tmp185 = (K923879532 * tmp111) + (K382683432 * tmp108);
	       tmp218 = tmp106 + tmp107;
	       tmp219 = tmp110 - tmp109;
	       tmp220 = (K382683432 * tmp218) - (K923879532 * tmp219);
	       tmp245 = (K382683432 * tmp219) + (K923879532 * tmp218);
	  }
     }
     {
	  fftw_real tmp25;
	  fftw_real tmp116;
	  fftw_real tmp88;
	  fftw_real tmp114;
	  fftw_real tmp28;
	  fftw_real tmp113;
	  fftw_real tmp91;
	  fftw_real tmp117;
	  fftw_real tmp319;
	  fftw_real tmp320;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp23;
	       fftw_real tmp24;
	       fftw_real tmp86;
	       fftw_real tmp87;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp23 = c_re(input[30 * istride]);
	       tmp24 = c_re(input[14 * istride]);
	       tmp25 = tmp23 + tmp24;
	       tmp116 = tmp23 - tmp24;
	       tmp86 = c_im(input[30 * istride]);
	       tmp87 = c_im(input[14 * istride]);
	       tmp88 = tmp86 + tmp87;
	       tmp114 = tmp86 - tmp87;
	  }
	  {
	       fftw_real tmp26;
	       fftw_real tmp27;
	       fftw_real tmp89;
	       fftw_real tmp90;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp26 = c_re(input[6 * istride]);
	       tmp27 = c_re(input[22 * istride]);
	       tmp28 = tmp26 + tmp27;
	       tmp113 = tmp26 - tmp27;
	       tmp89 = c_im(input[6 * istride]);
	       tmp90 = c_im(input[22 * istride]);
	       tmp91 = tmp89 + tmp90;
	       tmp117 = tmp89 - tmp90;
	  }
	  tmp29 = tmp25 + tmp28;
	  tmp319 = tmp25 - tmp28;
	  tmp320 = tmp88 - tmp91;
	  tmp321 = tmp319 + tmp320;
	  tmp342 = tmp320 - tmp319;
	  tmp92 = tmp88 + tmp91;
	  {
	       fftw_real tmp115;
	       fftw_real tmp118;
	       fftw_real tmp221;
	       fftw_real tmp222;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp115 = tmp113 + tmp114;
	       tmp118 = tmp116 - tmp117;
	       tmp119 = (K382683432 * tmp115) + (K923879532 * tmp118);
	       tmp184 = (K923879532 * tmp115) - (K382683432 * tmp118);
	       tmp221 = tmp114 - tmp113;
	       tmp222 = tmp116 + tmp117;
	       tmp223 = (K923879532 * tmp221) + (K382683432 * tmp222);
	       tmp244 = (K382683432 * tmp221) - (K923879532 * tmp222);
	  }
     }
     {
	  fftw_real tmp34;
	  fftw_real tmp122;
	  fftw_real tmp142;
	  fftw_real tmp279;
	  fftw_real tmp37;
	  fftw_real tmp139;
	  fftw_real tmp125;
	  fftw_real tmp280;
	  fftw_real tmp41;
	  fftw_real tmp127;
	  fftw_real tmp130;
	  fftw_real tmp282;
	  fftw_real tmp44;
	  fftw_real tmp132;
	  fftw_real tmp135;
	  fftw_real tmp283;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp32;
	       fftw_real tmp33;
	       fftw_real tmp140;
	       fftw_real tmp141;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp32 = c_re(input[istride]);
	       tmp33 = c_re(input[17 * istride]);
	       tmp34 = tmp32 + tmp33;
	       tmp122 = tmp32 - tmp33;
	       tmp140 = c_im(input[istride]);
	       tmp141 = c_im(input[17 * istride]);
	       tmp142 = tmp140 - tmp141;
	       tmp279 = tmp140 + tmp141;
	  }
	  {
	       fftw_real tmp35;
	       fftw_real tmp36;
	       fftw_real tmp123;
	       fftw_real tmp124;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp35 = c_re(input[9 * istride]);
	       tmp36 = c_re(input[25 * istride]);
	       tmp37 = tmp35 + tmp36;
	       tmp139 = tmp35 - tmp36;
	       tmp123 = c_im(input[9 * istride]);
	       tmp124 = c_im(input[25 * istride]);
	       tmp125 = tmp123 - tmp124;
	       tmp280 = tmp123 + tmp124;
	  }
	  {
	       fftw_real tmp39;
	       fftw_real tmp40;
	       fftw_real tmp128;
	       fftw_real tmp129;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp39 = c_re(input[5 * istride]);
	       tmp40 = c_re(input[21 * istride]);
	       tmp41 = tmp39 + tmp40;
	       tmp127 = tmp39 - tmp40;
	       tmp128 = c_im(input[5 * istride]);
	       tmp129 = c_im(input[21 * istride]);
	       tmp130 = tmp128 - tmp129;
	       tmp282 = tmp128 + tmp129;
	  }
	  {
	       fftw_real tmp42;
	       fftw_real tmp43;
	       fftw_real tmp133;
	       fftw_real tmp134;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp42 = c_re(input[29 * istride]);
	       tmp43 = c_re(input[13 * istride]);
	       tmp44 = tmp42 + tmp43;
	       tmp132 = tmp42 - tmp43;
	       tmp133 = c_im(input[29 * istride]);
	       tmp134 = c_im(input[13 * istride]);
	       tmp135 = tmp133 - tmp134;
	       tmp283 = tmp133 + tmp134;
	  }
	  {
	       fftw_real tmp327;
	       fftw_real tmp328;
	       fftw_real tmp324;
	       fftw_real tmp325;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp126 = tmp122 - tmp125;
	       tmp229 = tmp122 + tmp125;
	       tmp38 = tmp34 + tmp37;
	       tmp45 = tmp41 + tmp44;
	       tmp278 = tmp38 - tmp45;
	       tmp327 = tmp34 - tmp37;
	       tmp328 = tmp283 - tmp282;
	       tmp329 = tmp327 - tmp328;
	       tmp357 = tmp327 + tmp328;
	       tmp281 = tmp279 + tmp280;
	       tmp284 = tmp282 + tmp283;
	       tmp285 = tmp281 - tmp284;
	       tmp143 = tmp139 + tmp142;
	       tmp226 = tmp142 - tmp139;
	       tmp324 = tmp279 - tmp280;
	       tmp325 = tmp41 - tmp44;
	       tmp326 = tmp324 - tmp325;
	       tmp356 = tmp325 + tmp324;
	       {
		    fftw_real tmp131;
		    fftw_real tmp136;
		    fftw_real tmp144;
		    fftw_real tmp145;
		    ASSERT_ALIGNED_DOUBLE();
		    tmp131 = tmp127 - tmp130;
		    tmp136 = tmp132 + tmp135;
		    tmp137 = K707106781 * (tmp131 + tmp136);
		    tmp227 = K707106781 * (tmp131 - tmp136);
		    tmp144 = tmp127 + tmp130;
		    tmp145 = tmp135 - tmp132;
		    tmp146 = K707106781 * (tmp144 + tmp145);
		    tmp230 = K707106781 * (tmp145 - tmp144);
	       }
	  }
     }
     {
	  fftw_real tmp277;
	  fftw_real tmp301;
	  fftw_real tmp304;
	  fftw_real tmp306;
	  fftw_real tmp296;
	  fftw_real tmp300;
	  fftw_real tmp299;
	  fftw_real tmp305;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp275;
	       fftw_real tmp276;
	       fftw_real tmp302;
	       fftw_real tmp303;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp275 = tmp70 - tmp77;
	       tmp276 = tmp22 - tmp29;
	       tmp277 = tmp275 - tmp276;
	       tmp301 = tmp276 + tmp275;
	       tmp302 = tmp278 + tmp285;
	       tmp303 = tmp294 - tmp287;
	       tmp304 = K707106781 * (tmp302 + tmp303);
	       tmp306 = K707106781 * (tmp303 - tmp302);
	  }
	  {
	       fftw_real tmp286;
	       fftw_real tmp295;
	       fftw_real tmp297;
	       fftw_real tmp298;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp286 = tmp278 - tmp285;
	       tmp295 = tmp287 + tmp294;
	       tmp296 = K707106781 * (tmp286 - tmp295);
	       tmp300 = K707106781 * (tmp286 + tmp295);
	       tmp297 = tmp7 - tmp14;
	       tmp298 = tmp92 - tmp85;
	       tmp299 = tmp297 + tmp298;
	       tmp305 = tmp297 - tmp298;
	  }
	  c_im(output[28 * ostride]) = tmp277 - tmp296;
	  c_im(output[12 * ostride]) = tmp277 + tmp296;
	  c_re(output[20 * ostride]) = tmp299 - tmp300;
	  c_re(output[4 * ostride]) = tmp299 + tmp300;
	  c_im(output[20 * ostride]) = tmp301 - tmp304;
	  c_im(output[4 * ostride]) = tmp301 + tmp304;
	  c_re(output[28 * ostride]) = tmp305 - tmp306;
	  c_re(output[12 * ostride]) = tmp305 + tmp306;
     }
     {
	  fftw_real tmp31;
	  fftw_real tmp311;
	  fftw_real tmp310;
	  fftw_real tmp312;
	  fftw_real tmp62;
	  fftw_real tmp63;
	  fftw_real tmp94;
	  fftw_real tmp307;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp15;
	       fftw_real tmp30;
	       fftw_real tmp308;
	       fftw_real tmp309;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp15 = tmp7 + tmp14;
	       tmp30 = tmp22 + tmp29;
	       tmp31 = tmp15 + tmp30;
	       tmp311 = tmp15 - tmp30;
	       tmp308 = tmp281 + tmp284;
	       tmp309 = tmp290 + tmp293;
	       tmp310 = tmp308 + tmp309;
	       tmp312 = tmp309 - tmp308;
	  }
	  {
	       fftw_real tmp46;
	       fftw_real tmp61;
	       fftw_real tmp78;
	       fftw_real tmp93;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp46 = tmp38 + tmp45;
	       tmp61 = tmp53 + tmp60;
	       tmp62 = tmp46 + tmp61;
	       tmp63 = tmp46 - tmp61;
	       tmp78 = tmp70 + tmp77;
	       tmp93 = tmp85 + tmp92;
	       tmp94 = tmp78 - tmp93;
	       tmp307 = tmp78 + tmp93;
	  }
	  c_re(output[16 * ostride]) = tmp31 - tmp62;
	  c_re(output[0]) = tmp31 + tmp62;
	  c_im(output[8 * ostride]) = tmp63 + tmp94;
	  c_im(output[24 * ostride]) = tmp94 - tmp63;
	  c_im(output[16 * ostride]) = tmp307 - tmp310;
	  c_im(output[0]) = tmp307 + tmp310;
	  c_re(output[24 * ostride]) = tmp311 - tmp312;
	  c_re(output[8 * ostride]) = tmp311 + tmp312;
     }
     {
	  fftw_real tmp121;
	  fftw_real tmp189;
	  fftw_real tmp187;
	  fftw_real tmp193;
	  fftw_real tmp148;
	  fftw_real tmp190;
	  fftw_real tmp175;
	  fftw_real tmp191;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp105;
	       fftw_real tmp120;
	       fftw_real tmp183;
	       fftw_real tmp186;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp105 = tmp97 - tmp104;
	       tmp120 = tmp112 - tmp119;
	       tmp121 = tmp105 - tmp120;
	       tmp189 = tmp105 + tmp120;
	       tmp183 = tmp179 - tmp182;
	       tmp186 = tmp184 - tmp185;
	       tmp187 = tmp183 + tmp186;
	       tmp193 = tmp183 - tmp186;
	  }
	  {
	       fftw_real tmp138;
	       fftw_real tmp147;
	       fftw_real tmp165;
	       fftw_real tmp174;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp138 = tmp126 - tmp137;
	       tmp147 = tmp143 - tmp146;
	       tmp148 = (K555570233 * tmp138) - (K831469612 * tmp147);
	       tmp190 = (K831469612 * tmp138) + (K555570233 * tmp147);
	       tmp165 = tmp153 - tmp164;
	       tmp174 = tmp170 - tmp173;
	       tmp175 = (K555570233 * tmp165) + (K831469612 * tmp174);
	       tmp191 = (K555570233 * tmp174) - (K831469612 * tmp165);
	  }
	  {
	       fftw_real tmp176;
	       fftw_real tmp188;
	       fftw_real tmp192;
	       fftw_real tmp194;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp176 = tmp148 - tmp175;
	       c_im(output[29 * ostride]) = tmp121 - tmp176;
	       c_im(output[13 * ostride]) = tmp121 + tmp176;
	       tmp188 = tmp148 + tmp175;
	       c_re(output[21 * ostride]) = tmp187 - tmp188;
	       c_re(output[5 * ostride]) = tmp187 + tmp188;
	       tmp192 = tmp190 + tmp191;
	       c_im(output[21 * ostride]) = tmp189 - tmp192;
	       c_im(output[5 * ostride]) = tmp189 + tmp192;
	       tmp194 = tmp191 - tmp190;
	       c_re(output[29 * ostride]) = tmp193 - tmp194;
	       c_re(output[13 * ostride]) = tmp193 + tmp194;
	  }
     }
     {
	  fftw_real tmp197;
	  fftw_real tmp209;
	  fftw_real tmp207;
	  fftw_real tmp213;
	  fftw_real tmp200;
	  fftw_real tmp210;
	  fftw_real tmp203;
	  fftw_real tmp211;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp195;
	       fftw_real tmp196;
	       fftw_real tmp205;
	       fftw_real tmp206;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp195 = tmp97 + tmp104;
	       tmp196 = tmp185 + tmp184;
	       tmp197 = tmp195 - tmp196;
	       tmp209 = tmp195 + tmp196;
	       tmp205 = tmp179 + tmp182;
	       tmp206 = tmp112 + tmp119;
	       tmp207 = tmp205 + tmp206;
	       tmp213 = tmp205 - tmp206;
	  }
	  {
	       fftw_real tmp198;
	       fftw_real tmp199;
	       fftw_real tmp201;
	       fftw_real tmp202;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp198 = tmp126 + tmp137;
	       tmp199 = tmp143 + tmp146;
	       tmp200 = (K980785280 * tmp198) - (K195090322 * tmp199);
	       tmp210 = (K195090322 * tmp198) + (K980785280 * tmp199);
	       tmp201 = tmp153 + tmp164;
	       tmp202 = tmp170 + tmp173;
	       tmp203 = (K980785280 * tmp201) + (K195090322 * tmp202);
	       tmp211 = (K980785280 * tmp202) - (K195090322 * tmp201);
	  }
	  {
	       fftw_real tmp204;
	       fftw_real tmp208;
	       fftw_real tmp212;
	       fftw_real tmp214;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp204 = tmp200 - tmp203;
	       c_im(output[25 * ostride]) = tmp197 - tmp204;
	       c_im(output[9 * ostride]) = tmp197 + tmp204;
	       tmp208 = tmp200 + tmp203;
	       c_re(output[17 * ostride]) = tmp207 - tmp208;
	       c_re(output[ostride]) = tmp207 + tmp208;
	       tmp212 = tmp210 + tmp211;
	       c_im(output[17 * ostride]) = tmp209 - tmp212;
	       c_im(output[ostride]) = tmp209 + tmp212;
	       tmp214 = tmp211 - tmp210;
	       c_re(output[25 * ostride]) = tmp213 - tmp214;
	       c_re(output[9 * ostride]) = tmp213 + tmp214;
	  }
     }
     {
	  fftw_real tmp323;
	  fftw_real tmp347;
	  fftw_real tmp350;
	  fftw_real tmp352;
	  fftw_real tmp338;
	  fftw_real tmp346;
	  fftw_real tmp345;
	  fftw_real tmp351;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp315;
	       fftw_real tmp322;
	       fftw_real tmp348;
	       fftw_real tmp349;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp315 = tmp313 - tmp314;
	       tmp322 = K707106781 * (tmp318 - tmp321);
	       tmp323 = tmp315 + tmp322;
	       tmp347 = tmp315 - tmp322;
	       tmp348 = (K382683432 * tmp329) - (K923879532 * tmp326);
	       tmp349 = (K923879532 * tmp333) + (K382683432 * tmp336);
	       tmp350 = tmp348 - tmp349;
	       tmp352 = tmp348 + tmp349;
	  }
	  {
	       fftw_real tmp330;
	       fftw_real tmp337;
	       fftw_real tmp341;
	       fftw_real tmp344;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp330 = (K382683432 * tmp326) + (K923879532 * tmp329);
	       tmp337 = (K382683432 * tmp333) - (K923879532 * tmp336);
	       tmp338 = tmp330 + tmp337;
	       tmp346 = tmp337 - tmp330;
	       tmp341 = tmp339 - tmp340;
	       tmp344 = K707106781 * (tmp342 - tmp343);
	       tmp345 = tmp341 - tmp344;
	       tmp351 = tmp341 + tmp344;
	  }
	  c_im(output[22 * ostride]) = tmp323 - tmp338;
	  c_im(output[6 * ostride]) = tmp323 + tmp338;
	  c_re(output[30 * ostride]) = tmp345 - tmp346;
	  c_re(output[14 * ostride]) = tmp345 + tmp346;
	  c_im(output[30 * ostride]) = tmp347 - tmp350;
	  c_im(output[14 * ostride]) = tmp347 + tmp350;
	  c_re(output[22 * ostride]) = tmp351 - tmp352;
	  c_re(output[6 * ostride]) = tmp351 + tmp352;
     }
     {
	  fftw_real tmp355;
	  fftw_real tmp367;
	  fftw_real tmp370;
	  fftw_real tmp372;
	  fftw_real tmp362;
	  fftw_real tmp366;
	  fftw_real tmp365;
	  fftw_real tmp371;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp353;
	       fftw_real tmp354;
	       fftw_real tmp368;
	       fftw_real tmp369;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp353 = tmp314 + tmp313;
	       tmp354 = K707106781 * (tmp343 + tmp342);
	       tmp355 = tmp353 + tmp354;
	       tmp367 = tmp353 - tmp354;
	       tmp368 = (K923879532 * tmp357) - (K382683432 * tmp356);
	       tmp369 = (K382683432 * tmp359) + (K923879532 * tmp360);
	       tmp370 = tmp368 - tmp369;
	       tmp372 = tmp368 + tmp369;
	  }
	  {
	       fftw_real tmp358;
	       fftw_real tmp361;
	       fftw_real tmp363;
	       fftw_real tmp364;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp358 = (K923879532 * tmp356) + (K382683432 * tmp357);
	       tmp361 = (K923879532 * tmp359) - (K382683432 * tmp360);
	       tmp362 = tmp358 + tmp361;
	       tmp366 = tmp361 - tmp358;
	       tmp363 = tmp339 + tmp340;
	       tmp364 = K707106781 * (tmp318 + tmp321);
	       tmp365 = tmp363 - tmp364;
	       tmp371 = tmp363 + tmp364;
	  }
	  c_im(output[18 * ostride]) = tmp355 - tmp362;
	  c_im(output[2 * ostride]) = tmp355 + tmp362;
	  c_re(output[26 * ostride]) = tmp365 - tmp366;
	  c_re(output[10 * ostride]) = tmp365 + tmp366;
	  c_im(output[26 * ostride]) = tmp367 - tmp370;
	  c_im(output[10 * ostride]) = tmp367 + tmp370;
	  c_re(output[18 * ostride]) = tmp371 - tmp372;
	  c_re(output[2 * ostride]) = tmp371 + tmp372;
     }
     {
	  fftw_real tmp225;
	  fftw_real tmp249;
	  fftw_real tmp247;
	  fftw_real tmp253;
	  fftw_real tmp232;
	  fftw_real tmp250;
	  fftw_real tmp239;
	  fftw_real tmp251;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp217;
	       fftw_real tmp224;
	       fftw_real tmp243;
	       fftw_real tmp246;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp217 = tmp215 - tmp216;
	       tmp224 = tmp220 - tmp223;
	       tmp225 = tmp217 + tmp224;
	       tmp249 = tmp217 - tmp224;
	       tmp243 = tmp241 - tmp242;
	       tmp246 = tmp244 - tmp245;
	       tmp247 = tmp243 - tmp246;
	       tmp253 = tmp243 + tmp246;
	  }
	  {
	       fftw_real tmp228;
	       fftw_real tmp231;
	       fftw_real tmp235;
	       fftw_real tmp238;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp228 = tmp226 - tmp227;
	       tmp231 = tmp229 - tmp230;
	       tmp232 = (K195090322 * tmp228) + (K980785280 * tmp231);
	       tmp250 = (K195090322 * tmp231) - (K980785280 * tmp228);
	       tmp235 = tmp233 - tmp234;
	       tmp238 = tmp236 - tmp237;
	       tmp239 = (K195090322 * tmp235) - (K980785280 * tmp238);
	       tmp251 = (K980785280 * tmp235) + (K195090322 * tmp238);
	  }
	  {
	       fftw_real tmp240;
	       fftw_real tmp248;
	       fftw_real tmp252;
	       fftw_real tmp254;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp240 = tmp232 + tmp239;
	       c_im(output[23 * ostride]) = tmp225 - tmp240;
	       c_im(output[7 * ostride]) = tmp225 + tmp240;
	       tmp248 = tmp239 - tmp232;
	       c_re(output[31 * ostride]) = tmp247 - tmp248;
	       c_re(output[15 * ostride]) = tmp247 + tmp248;
	       tmp252 = tmp250 - tmp251;
	       c_im(output[31 * ostride]) = tmp249 - tmp252;
	       c_im(output[15 * ostride]) = tmp249 + tmp252;
	       tmp254 = tmp250 + tmp251;
	       c_re(output[23 * ostride]) = tmp253 - tmp254;
	       c_re(output[7 * ostride]) = tmp253 + tmp254;
	  }
     }
     {
	  fftw_real tmp257;
	  fftw_real tmp269;
	  fftw_real tmp267;
	  fftw_real tmp273;
	  fftw_real tmp260;
	  fftw_real tmp270;
	  fftw_real tmp263;
	  fftw_real tmp271;
	  ASSERT_ALIGNED_DOUBLE();
	  {
	       fftw_real tmp255;
	       fftw_real tmp256;
	       fftw_real tmp265;
	       fftw_real tmp266;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp255 = tmp215 + tmp216;
	       tmp256 = tmp245 + tmp244;
	       tmp257 = tmp255 + tmp256;
	       tmp269 = tmp255 - tmp256;
	       tmp265 = tmp241 + tmp242;
	       tmp266 = tmp220 + tmp223;
	       tmp267 = tmp265 - tmp266;
	       tmp273 = tmp265 + tmp266;
	  }
	  {
	       fftw_real tmp258;
	       fftw_real tmp259;
	       fftw_real tmp261;
	       fftw_real tmp262;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp258 = tmp226 + tmp227;
	       tmp259 = tmp229 + tmp230;
	       tmp260 = (K831469612 * tmp258) + (K555570233 * tmp259);
	       tmp270 = (K831469612 * tmp259) - (K555570233 * tmp258);
	       tmp261 = tmp233 + tmp234;
	       tmp262 = tmp236 + tmp237;
	       tmp263 = (K831469612 * tmp261) - (K555570233 * tmp262);
	       tmp271 = (K555570233 * tmp261) + (K831469612 * tmp262);
	  }
	  {
	       fftw_real tmp264;
	       fftw_real tmp268;
	       fftw_real tmp272;
	       fftw_real tmp274;
	       ASSERT_ALIGNED_DOUBLE();
	       tmp264 = tmp260 + tmp263;
	       c_im(output[19 * ostride]) = tmp257 - tmp264;
	       c_im(output[3 * ostride]) = tmp257 + tmp264;
	       tmp268 = tmp263 - tmp260;
	       c_re(output[27 * ostride]) = tmp267 - tmp268;
	       c_re(output[11 * ostride]) = tmp267 + tmp268;
	       tmp272 = tmp270 - tmp271;
	       c_im(output[27 * ostride]) = tmp269 - tmp272;
	       c_im(output[11 * ostride]) = tmp269 + tmp272;
	       tmp274 = tmp270 + tmp271;
	       c_re(output[19 * ostride]) = tmp273 - tmp274;
	       c_re(output[3 * ostride]) = tmp273 + tmp274;
	  }
     }
}

fftw_codelet_desc fftwi_no_twiddle_32_desc =
{
     "fftwi_no_twiddle_32",
     (void (*)()) fftwi_no_twiddle_32,
     32,
     FFTW_BACKWARD,
     FFTW_NOTW,
     521,
     0,
     (const int *) 0,
};

max_stars_repo_name: jlprieur/jlplib

max_issues_repo_name: jlprieur/jlplib

max_forks_repo_name: jlprieur/jlplib


text:
#pragma once
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_cdf.h>

#include <QVector>

#include <algorithm>
#include <tuple>

class Gaussian1DFit
{
    struct data
    {
        double* x;
        double* y;
        size_t n;
    };

    enum
    {
        numOfPara = 4
    };
    /* model function: a * exp( -1/2 * [ (t - b) / c ]^2 ) + d */
private:
    gsl_vector*                      p0; /* initial fitting parameters:amplitude,center,width,offset*/
    gsl_vector*                      p; /* fitting parameters */
    gsl_vector*                      f; /* residual yi-y(xi) */
    gsl_matrix*                      covar; /*covariance matrix, not yet multiplied by sigma, so is not variance-covariance matrix*/

    gsl_multifit_nlinear_fdf         fdf;
    gsl_multifit_nlinear_parameters  fdf_params;
    gsl_multifit_nlinear_workspace*  work;

    data                             fit_data;

    size_t                           max_iter;
    double                           ptol; /* tolerance on fitting parameter p */
    double                           gtol; /* tolerance on gradient */
    double                           ftol;

    int                              info;/*fiting stop reason, for debug*/
    double                           rss; /*mean of sum of residual square, dof corrected*/
    gsl_vector*                      confid95;

public:
    Gaussian1DFit() {};
    Gaussian1DFit(size_t n, double* datax, double* datay,
        double a0, double b0, double c0, double d0,
        size_t max_iter = 200, double ptol = 1.0e-8, double gtol = 1.0e-8, double ftol = 1.0e-8);
    ~Gaussian1DFit();
    void solve_system();
    void set_data(double* datay);
    void set_data(size_t n, double* datax, double* datay, bool free = true);
    void set_initialP(double a0, double b0, double c0, double d0);

    static double gaussian(const double a, const double b, const double c, const double d, const double t);
    QVector<double> calcFittedGaussian();

    const QVector<double> fittedPara() const;
    const QVector<double> confidence95Interval() const;

private:
    static int func_f (const gsl_vector* p, void* datafit, gsl_vector* f);
    static int func_df(const gsl_vector* p, void* datafit, gsl_matrix* J);
    static int func_fvv(const gsl_vector* p, const gsl_vector* v,
        void* datafit, gsl_vector* fvv);

};

max_stars_repo_name: zzpwahaha/VimbaCamJILA

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
/**
 * nn.c
 *
 * Created by mizu-bai
 *
 * This is a implementation of Tensors Warm-up: numpy in LEARNING PYTORCH WITH EXAMPLES
 * (https://pytorch.org/tutorials/beginner/pytorch_with_examples.html#) using C with GSL.
 */

#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>

int main(int argc, const char *argv[]) {
    // create random input and output data
    int size = 2000;
    gsl_vector *x = gsl_vector_alloc(size);
    gsl_vector *y = gsl_vector_alloc(size);
    for (int i = 0; i < size; i++) {
        double x_tmp = 2 * M_PI / (double)(size - 1) * (double)i - M_PI;
        gsl_vector_set(x, i, x_tmp);
        gsl_vector_set(y, i, sin(x_tmp));
    }

    // randomly init weights
    int w_size = 4;
    gsl_vector *w = gsl_vector_alloc(w_size);
    gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
    for (int i = 0; i < w_size; i++) {
        gsl_vector_set(w, i, gsl_ran_gaussian(r, 1.0));
    }
    gsl_rng_free(r);

    double learning_rate = 1e-6;
    for(int t = 0; t < 2000; t++) {
        // forward pass, compute predicted y & compute loss
        gsl_vector *y_pred = gsl_vector_alloc(size);
        double loss = 0;
        for (int i = 0; i < size; i++) {
            double w_tmp = gsl_poly_eval(w->data, w_size, gsl_vector_get(x, i));
            gsl_vector_set(y_pred, i, w_tmp);
            loss += gsl_pow_2(gsl_vector_get(y_pred, i) - gsl_vector_get(y, i));
        }
        if (t % 100 == 99) {
            printf("t: %d, loss: %.4f\n", t, loss);
        }

        // compute gradients of w with respect to loss
        gsl_vector *grad_y_pred = gsl_vector_alloc(size);
        gsl_vector_memcpy(grad_y_pred, y_pred);
        gsl_vector_sub(grad_y_pred, y);
        gsl_vector_scale(grad_y_pred, 2);
        gsl_vector *grad_w = gsl_vector_alloc(w_size);
        gsl_vector_set_all(grad_w, 0.0);
        for(int i = 0; i < size; i++) {
            gsl_vector *grad_y_pred_tmp = gsl_vector_alloc(size);
            gsl_vector *x_tmp = gsl_vector_alloc(size);

            gsl_vector_set(grad_w, 0, gsl_vector_get(grad_w, 0) + gsl_vector_get(grad_y_pred, i));

            gsl_vector_memcpy(grad_y_pred_tmp, grad_y_pred);
            gsl_vector_mul(grad_y_pred_tmp, x);
            gsl_vector_set(grad_w, 1, gsl_vector_get(grad_w, 1) + gsl_vector_get(grad_y_pred_tmp, i));

            gsl_vector_memcpy(grad_y_pred_tmp, grad_y_pred);
            gsl_vector_memcpy(x_tmp, x);
            gsl_vector_mul(x_tmp, x);
            gsl_vector_mul(grad_y_pred_tmp, x_tmp);
            gsl_vector_set(grad_w, 2, gsl_vector_get(grad_w, 2) + gsl_vector_get(grad_y_pred_tmp, i));

            gsl_vector_memcpy(grad_y_pred_tmp, grad_y_pred);
            gsl_vector_memcpy(x_tmp, x);
            gsl_vector_mul(x_tmp, x);
            gsl_vector_mul(x_tmp, x);
            gsl_vector_mul(grad_y_pred_tmp, x_tmp);
            gsl_vector_set(grad_w, 3, gsl_vector_get(grad_w, 3) + gsl_vector_get(grad_y_pred_tmp, i));

            gsl_vector_free(grad_y_pred_tmp);
            gsl_vector_free(x_tmp);
        }

        gsl_vector_scale(grad_w, learning_rate);
        gsl_vector_sub(w, grad_w);

        // free
        gsl_vector_free(y_pred);
        gsl_vector_free(grad_y_pred);
        gsl_vector_free(grad_w);
    }

    printf("Result: y = %.4f + %.4f * x + %.4f * x^2 + %.4f * x^3\n", gsl_vector_get(w, 0), gsl_vector_get(w, 1), gsl_vector_get(w, 2), gsl_vector_get(w, 3));

    // free
    gsl_vector_free(x);
    gsl_vector_free(y);
    gsl_vector_free(w);

    return 0;
}

max_stars_repo_name: mizu-bai/Data-Science-Demos

max_issues_repo_name: mizu-bai/Data-Science-Demos

max_forks_repo_name: mizu-bai/Data-Science-Demos


text:
#ifndef LIB_INCLUDE_TICK_ARRAY_VECTOR_OPERATIONS_H_
#define LIB_INCLUDE_TICK_ARRAY_VECTOR_OPERATIONS_H_

// License: BSD 3 clause

#include <numeric>

#include <algorithm>
#include <atomic>
#include <type_traits>

#include "promote.h"
#include "tick/base/defs.h"

namespace tick {
namespace detail {

template <typename T>
struct DLL_PUBLIC vector_operations_unoptimized {
  template <typename K>
  tick::promote_t<K> sum(const ulong n, const T *x) const;

  template <typename K>
  typename std::enable_if<std::is_same<T, std::atomic<K>>::value, T>::type dot(
      const ulong n, const T *x, const K *y) const;

  template <typename K>
  typename std::enable_if<std::is_same<T, std::atomic<K>>::value, K>::type dot(
      const ulong n, const K *x, const T *y) const;

  template <typename K>
  typename std::enable_if<!std::is_same<T, std::atomic<K>>::value, T>::type dot(
      const ulong n, const T *x, const K *y) const;

  template <typename K>
  typename std::enable_if<std::is_same<T, std::atomic<K>>::value>::type scale(
      const ulong n, const K alpha, T *x) const;

  template <typename K>
  typename std::enable_if<!std::is_same<T, std::atomic<K>>::value>::type scale(
      const ulong n, const K alpha, T *x) const;

  template <typename K>
  typename std::enable_if<std::is_same<T, std::atomic<K>>::value>::type set(
      const ulong n, const K alpha, T *x) const;

  template <typename K>
  typename std::enable_if<!std::is_same<T, std::atomic<K>>::value>::type set(
      const ulong n, const K alpha, T *x) const;

  template <typename K, typename Y>
  typename std::enable_if<std::is_same<T, std::atomic<K>>::value &&
                          !std::is_same<Y, std::atomic<K>>::value>::type
  mult_incr(const uint64_t n, const K alpha, const Y *x, T *y) const;

  template <typename K, typename Y>
  typename std::enable_if<std::is_same<Y, std::atomic<K>>::value &&
                          !std::is_same<T, std::atomic<K>>::value>::type
  mult_incr(const uint64_t n, const K alpha, const Y *x, T *y) const;

  template <typename K, typename Y>
  typename std::enable_if<std::is_same<T, std::atomic<K>>::value &&
                          std::is_same<Y, std::atomic<K>>::value>::type
  mult_incr(const uint64_t n, const K alpha, const Y *x, T *y) const;

  template <typename K, typename Y>
  typename std::enable_if<std::is_same<T, K>::value &&
                          std::is_same<Y, K>::value>::type
  mult_incr(const uint64_t n, const K alpha, const Y *x, T *y) const;
};

}  // namespace detail
}  // namespace tick

#if !defined(TICK_CBLAS_AVAILABLE)

namespace tick {

template <typename T>
using vector_operations = detail::vector_operations_unoptimized<T>;

}  // namespace tick

#else  // if defined(TICK_CBLAS_AVAILABLE)

// Find available blas distribution
#if defined(TICK_USE_MKL)

#include <mkl.h>

#elif defined(__APPLE__)

#include <Accelerate/Accelerate.h>

// TODO(svp) Disabling this feature until we find a good way to determine if
// ATLAS is actually available #define XDATA_CATLAS_AVAILABLE

#else

extern "C" {
#include <cblas.h>
}

#endif  // defined(__APPLE__)

#include "tick/array/vector/cblas.hpp"

#endif  // if !defined(TICK_CBLAS_AVAILABLE)

#include "tick/array/vector/un_optimized.hpp"

#endif  // LIB_INCLUDE_TICK_ARRAY_VECTOR_OPERATIONS_H_


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

#include "../inc/knnring.h"

knnresult kNN(double *X, double *Y, int n, int m, int d, int k){

	knnresult k_nearest;

	//initalize the result's variable
	k_nearest.nidx  = (int *)malloc(m*k*sizeof(int));
	k_nearest.ndist = (double *)malloc(m*k*sizeof(double));
	k_nearest.m = m;
	k_nearest.k = k;

	  //compute the distance Matrix D
	double *D = compute_D(X,Y,n,m,d,k);

	 //select k-nearest neighbors for every point in QUERY
	k_select(D,k_nearest.ndist, k_nearest.nidx,n,m,k);
	return k_nearest;
}

double *compute_D(double *X, double *Y, int n, int m, int d, int k){

	//compute distances using
	//D = (X*X)ee_t - 2XY_t + ee_t(Y*Y)_t
	double *tempD   = (double *)malloc(n*m*sizeof(double));
	double *e_1     = (double *)malloc(d*m*sizeof(double));
	double *XX      = (double *)malloc(n*d*sizeof(double));

	for(int i=0;i<d*m;i++)
		e_1[i]=1;

	for(int i=0;i<n*d;i++)
		XX[i] = X[i]*X[i];

	cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n,m,d,1,XX,d,e_1,m,0,tempD,m);

	free(e_1);
	free(XX);

	double *e_2 = (double *)malloc(n*d*sizeof(double));
	double *YY = (double *)malloc(m*d*sizeof(double));

	for(int i=0;i<n*d;i++)
		e_2[i]=1;

	for(int i=0;i<m*d;i++)
		YY[i] = Y[i]*Y[i];

	cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, n,m,d,-2,X,d,Y,d,1,tempD,m);
	cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, n,m,d,1,e_2,d,YY,d,1,tempD,m);

	for(int i=0;i<n*m;i++)
		tempD[i] = sqrt(fabs(tempD[i]));

	free(e_2);
	free(YY);

	return tempD;
}

void k_select(double *D, double *ndist, int *nidx,int n, int m, int k){
	for(int j=0;j<m;j++){
		for(int i=0;i<k;i++){
			ndist[j*k+i] = D[j+i*m];
			nidx[j*k+i]  = i;
		}
		//sort the initial k neighbors
		quick_sort(ndist,nidx,j*k,k*(j+1)-1);
	}

	//search the rest points to find it there are closer neighbors
	for(int j=0;j<m;j++){
		for(int i=k;i<n;i++){
			if(D[j+i*m] < ndist[k*(j+1)-1]){
				ndist[k*(j+1)-1] = D[j+i*m];
				nidx[k*(j+1)-1] = i;
				quick_sort(ndist,nidx,j*k,(j+1)*k-1);
			}
		}
	}
}

void quick_sort(double *k_dist, int *k_neigh, int low, int high){
	if(low<high){

		//actuall sorting an element
		int pi = partition(k_dist,k_neigh,low,high);

		//recursivelly call quick_sort
		quick_sort(k_dist,k_neigh,low,pi-1);
		quick_sort(k_dist,k_neigh,pi+1,high);
	}
}

int partition(double *k_dist, int *k_neigh,int low, int high){

	double temp_dist;
	int temp_point;
	int checker = 0;

	double pivot = k_dist[high];
	//printf("the pivot is %lf \n",pivot);
	int i = (low-1);
	for(int j=low;j<=high-1;j++){
		if(k_dist[j]<=pivot){
			/*swap the two elements*/

			i=i+1;
			temp_dist  = k_dist[i];
			temp_point = k_neigh[i];

			k_dist[i]  = k_dist[j];
			k_neigh[i] = k_neigh[j];

			k_dist[j]  = temp_dist;
			k_neigh[j] = temp_point;
			checker = 1;
		}
	}

	/* swap the pivot */
	i=i+1;
	temp_dist  = k_dist[i];
	temp_point = k_neigh[i];

	k_dist[i]  = k_dist[high];
	k_neigh[i] = k_neigh[high];

	k_dist[high]  = temp_dist;
	k_neigh[high] = temp_point;

	return i;
}

void create_points(double *X,int size){
	double *temp;
	temp = X;
	srand(time(NULL));
	for(int i=0;i<size;i++)
		temp[i] = (double) 1000*rand()/RAND_MAX;
}

void write_to_file(char *str1, int *nidx, double *ndist, int amount, int k){
  FILE *fp;
  //file for indexes
  fp = fopen(str1,"a+");

  int counter = 1;
  for(int i=0;i<amount;i++){
    fprintf(fp,"%d)%d -- %lf\n",counter,nidx[i],ndist[i]);
    counter ++;
	if((i+1)%k == 0){
      fprintf(fp,"k nearest neighbors for next point \n");
      counter = 1;
    }
  }
	fclose(fp);
}

//######################### V1 FUNCTIONS ############################

knnresult distrAllkNN(double *X, int n, int d, int k){

	//results variables
	 knnresult knn_process;
	knnresult knn_temp;

	int p_rank;                            //PROCESS RANK
	int p_size;                            //TOTAL NUMBER OF PROCESSES
	int next,prev;                         //NEXT AND PREVIOUS PROCESSES
	double *receive_buff;                  //BUFFERS FOR COMMUNICATION
	double *Y;                             //QUERY POINTS

	//variables for MPI COMMUNICATION.
	MPI_Status status;
	MPI_Request send_request,receive_request;
	int tag;

	//time variables
	struct timespec ts_start;
	struct timespec ts_end;
	double time;

	//find the number of processes
	MPI_Comm_size(MPI_COMM_WORLD, &p_size);

	//find the id of the process
	MPI_Comm_rank(MPI_COMM_WORLD, &p_rank);

	//define the communication tag
	tag = 1;

	//define the next and the previous process.
	//and create the circular communication
	next = p_rank+1;
	if(next == p_size) next = 0;

	prev = p_rank-1;
	if(prev== -1) prev = p_size-1;

	//initialize Y block
	Y = (double *)malloc(n*d*sizeof(double));

	//in the first iteration every process computes distances between its own points
	memcpy(Y,X,n*d*sizeof(double));

	//initialize receive buffer
	receive_buff = (double *)malloc(n*d*sizeof(double));

	//MASK COMMUNICATION TIME
	//Send and receive point for the next iteration
	//clock_gettime(CLOCK_MONOTONIC,&ts_start);

	MPI_Isend(Y,n*d,MPI_DOUBLE,next,tag,MPI_COMM_WORLD, &send_request);
	MPI_Irecv(receive_buff,n*d, MPI_DOUBLE, prev, tag, MPI_COMM_WORLD, &receive_request);

	//find k-nearest neighbors in this set of points.
	knn_process = kNN(X,Y,n,n,d,k);

	for(int i=0;i<n*k;i++) knn_process.nidx[i] = knn_process.nidx[i]+n*p_rank;

	MPI_Wait(&send_request,&status);
	//make sure receive_buffer is ready
	MPI_Wait(&receive_request,&status);

	//clock_gettime(CLOCK_MONOTONIC,&ts_end);
	//time = 1000000*(double)(ts_end.tv_sec-ts_start.tv_sec)+(double)(ts_end.tv_nsec-ts_start.tv_nsec)/1000;
	//printf("communication time and computation time%lf\n",time);

  for(int i=0;i<p_size-1;i++){
	  memcpy(Y,receive_buff,n*d*sizeof(double));

	  //MASK COMMUNICATION TIME
	  //send the block to the next node
	  clock_gettime(CLOCK_MONOTONIC,&ts_start);
	  MPI_Isend(Y,n*d,MPI_DOUBLE,next,tag,MPI_COMM_WORLD, &send_request);
	  MPI_Irecv(receive_buff,n*d,MPI_DOUBLE,prev,tag,MPI_COMM_WORLD,&receive_request);

	  knn_temp = kNN(Y,X,n,n,d,k);

	  //compare k-nearest points
	  knn_combine(knn_process, knn_temp,n,k,(p_rank-i-1+p_size)%p_size);

	  //make sure Y buffer has been send in order to change it
	  MPI_Wait(&send_request,&status);
	  //make sure receive_buffer is ready
	  MPI_Wait(&receive_request,&status);

	  //clock_gettime(CLOCK_MONOTONIC,&ts_end);
	  //time = 1000000*(double)(ts_end.tv_sec-ts_start.tv_sec)+(double)(ts_end.tv_nsec-ts_start.tv_nsec)/1000;
	  //printf("communication time and knn compute time%lf\n",time);
  }
  free(Y);
  free(knn_temp.nidx);
  free(knn_temp.ndist);

  for(int i=0;i<knn_process.m;i++) quick_sort(knn_process.ndist,knn_process.nidx,i*k,k*(i+1)-1);
  return knn_process;
}

void knn_combine(knnresult knn_total,knnresult knn_temp, int n, int k, int no_block){
	for(int m=0;m<knn_total.m;m++){
		for(int j=0;j<k;j++){
			if(knn_temp.ndist[m*k+j]>knn_total.ndist[k*(m+1)-1])
				break;
			else{
				knn_total.ndist[k*(m+1)-1] = knn_temp.ndist[m*k+j];
				knn_total.nidx[k*(m+1)-1] = knn_temp.nidx[m*k+j] + no_block*n;

				//sort the new element
				quick_sort(knn_total.ndist,knn_total.nidx,m*k,(m+1)*k-1);
			}
		}
	}
}


text:
/**
 * Created              :   2020.09.02;
 * Last Update          :   2021.06.11;
 * Author               :   Gabriel Marino de Oliveira <ra115114@uem.br>;
 * Supervisor/Advisor   :   Breno Ferraz de Oliveira <>;
 * Notes                :
 * Compilation          :   gcc -O3 src/main.c -static -lgls -lglscblas -lm -o simulationRPS.out
 * Execution            :   time ./simulationRPS.out
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>           //  add mathematical operations;
#include <time.h>           //  add functions to work with time and date, in this case used only to get the seconds counter from 01/01/1970;
#include <complex.h>        //  add complex number operations;
#include <gsl/gsl_rng.h>    //  add pseudo-random number generators;

#define Ni  500     //  Size of the grid in "i", or abscissa axis;
#define Nj  500     //  Size of the grid in "j", or ordinate axis;
#define dt  0.1     //  Time interval;
#define tf  10000   //  Total Time, or number of generations;
#define Pp  0.3     //  Predation probability;
#define Pm  0.4     //  Mobility probability for sthocastic method;
                    /** Is not necessary to sthocastic method define the Reproduction probability, here we assume which "Pp + Pm + Pr = 1" or "Pr = 1 - (Pp + Pm)",
                     * and this is why Reproduction probability is only needed to define for RK method**/
#define Cm  1.0     //  Mobility coeficient for RK method;
#define Pr  0.4     //  Reproduction probability;

#define np  100                                 //  Number of points which will be printed when used a logarithmic scale;
// const double r = pow((tf/dt), (1.0/np));        //  r stands for ratio, will used to generate a arithmetic progression, useful to print equally spaced point


text:
#ifndef __INC_EMULATORSTRUCT__
#define __INC_EMULATORSTRUCT__

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

/**
 * \@file emulator_struct.h
 * \brief defines the emulator_struct which contains the details of a working emulator at runtime
 */

struct modelstruct;

/*********************************************************************
This data structure holds pointers to the information used to execute
the emulator at a new point.  This is especially helpful since I want
to stream in queries and get results streamed back out.
*********************************************************************/
typedef struct emulator_struct {
	int nparams;
	int nmodel_points;
	int nregression_fns;
	int nthetas;
	struct modelstruct * model;
	gsl_matrix * cinverse;
	gsl_vector * beta_vector;
	gsl_matrix * h_matrix;
} emulator_struct;

#include "modelstruct.h"

emulator_struct * alloc_emulator_struct(modelstruct * model);
void free_emulator_struct(emulator_struct * e);

void makeHMatrix_es(gsl_matrix *h_matrix, emulator_struct *e);
void makeCovMatrix_es(gsl_matrix *cov_matrix, emulator_struct *e);
void makeKVector_es(gsl_vector *kvector, gsl_vector *point, emulator_struct *e);
void estimateBeta_es(gsl_vector *beta_vector, emulator_struct *e);
void emulate_point(emulator_struct* e, gsl_vector * point,  double * mean, double * variance);
#endif

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: MADAI/MADAIEmulator

max_forks_repo_name: MADAI/MADAIEmulator