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

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

#include "ccl.h"

/*--------ROUTINE: taper_cl ------
TASK:n Apply cosine tapering to Cls to reduce aliasing
INPUT: number of ell bins for Cl, ell vector, C_ell vector, limits for tapering
       e.g., ell_limits=[low_ell_limit_lower,low_ell_limit_upper,high_ell_limit_lower,high_ell_limit_upper]
*/
static int taper_cl(int n_ell,double *ell,double *cl, double *ell_limits)
{

  for(int i=0;i<n_ell;i++) {
    if(ell[i]<ell_limits[0] || ell[i]>ell_limits[3]) {
      cl[i]=0;//ell outside desirable range
      continue;
    }
    if(ell[i]>=ell_limits[1] && ell[i]<=ell_limits[2])
      continue;//ell within good ell range

    if(ell[i]<ell_limits[1])//tapering low ell
      cl[i]*=cos((ell[i]-ell_limits[1])/(ell_limits[1]-ell_limits[0])*M_PI/2.);

    if(ell[i]>ell_limits[2])//tapering high ell
      cl[i]*=cos((ell[i]-ell_limits[2])/(ell_limits[3]-ell_limits[2])*M_PI/2.);
  }

  return 0;
}

/*--------ROUTINE: ccl_tracer_corr_fftlog ------
TASK: For a given tracer, get the correlation function
      Following function takes a function to calculate angular cl as well.
      By default above function will call it using ccl_angular_cl
INPUT: type of tracer, number of theta values to evaluate = NL, theta vector
 */
static void ccl_tracer_corr_fftlog(ccl_cosmology *cosmo,
                                   int n_ell,double *ell,double *cls,
                                   int n_theta,double *theta,double *wtheta,
                                   int corr_type,int do_taper_cl,double *taper_cl_limits,
                                   int *status) {
  int i;
  double *l_arr,*cl_arr,*th_arr,*wth_arr;

  l_arr=ccl_log_spacing(cosmo->spline_params.ELL_MIN_CORR,cosmo->spline_params.ELL_MAX_CORR,cosmo->spline_params.N_ELL_CORR);
  if(l_arr==NULL) {
    *status=CCL_ERROR_LINSPACE;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }
  cl_arr=malloc(cosmo->spline_params.N_ELL_CORR*sizeof(double));
  if(cl_arr==NULL) {
    free(l_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }

  //Interpolate input Cl into array needed for FFTLog
  ccl_f1d_t *cl_spl=ccl_f1d_t_new(n_ell,ell,cls,cls[0],0,
				  ccl_f1d_extrap_const,
				  ccl_f1d_extrap_logx_logy,
				  status);
  if(cl_spl==NULL) {
    free(l_arr);
    free(cl_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo,
				     "ccl_correlation.c: ccl_tracer_corr_fftlog "
				     "ran out of memory\n");
    return;
  }

  for(i=0;i<cosmo->spline_params.N_ELL_CORR;i++)
    cl_arr[i]=ccl_f1d_t_eval(cl_spl,l_arr[i]);
  ccl_f1d_t_free(cl_spl);

  if (do_taper_cl)
    taper_cl(cosmo->spline_params.N_ELL_CORR,l_arr,cl_arr,taper_cl_limits);

  th_arr=malloc(sizeof(double)*cosmo->spline_params.N_ELL_CORR);
  if(th_arr==NULL) {
    free(l_arr);
    free(cl_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }
  wth_arr=(double *)malloc(sizeof(double)*cosmo->spline_params.N_ELL_CORR);
  if(wth_arr==NULL) {
    free(l_arr);
    free(cl_arr);
    free(th_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }

  for(i=0;i<cosmo->spline_params.N_ELL_CORR;i++)
    th_arr[i]=0;
  //Although set here to 0, theta is modified by FFTlog to obtain the correlation at ~1/l

  int i_bessel=0;
  if(corr_type==CCL_CORR_GG) i_bessel=0;
  if(corr_type==CCL_CORR_GL) i_bessel=2;
  if(corr_type==CCL_CORR_LP) i_bessel=0;
  if(corr_type==CCL_CORR_LM) i_bessel=4;
  ccl_fftlog_ComputeXi2D(i_bessel,0,
			 1, cosmo->spline_params.N_ELL_CORR,l_arr,&cl_arr,
			 th_arr,&wth_arr, status);

  // Interpolate to output values of theta
  ccl_f1d_t *wth_spl=ccl_f1d_t_new(cosmo->spline_params.N_ELL_CORR,th_arr,
				   wth_arr,wth_arr[0],0,
				   ccl_f1d_extrap_const,
				   ccl_f1d_extrap_const, status);
  if (wth_spl == NULL) {
    free(l_arr);
    free(cl_arr);
    free(th_arr);
    free(wth_arr);
    *status = CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }
  for(i=0;i<n_theta;i++)
    wtheta[i]=ccl_f1d_t_eval(wth_spl,theta[i]*M_PI/180.);
  ccl_f1d_t_free(wth_spl);

  free(l_arr);
  free(cl_arr);
  free(th_arr);
  free(wth_arr);

  return;
}

typedef struct {
  ccl_f1d_t *cl_spl;
  int i_bessel;
  double th;
} corr_int_par;

static double corr_bessel_integrand(double l,void *params)
{
  double cl,jbes;
  corr_int_par *p=(corr_int_par *)params;
  double x=l*p->th;

  cl=ccl_f1d_t_eval(p->cl_spl,l);

  jbes=gsl_sf_bessel_Jn(p->i_bessel,x);

  return l*jbes*cl;
}

static void ccl_tracer_corr_bessel(ccl_cosmology *cosmo,
                                   int n_ell,double *ell,double *cls,
                                   int n_theta,double *theta,double *wtheta,
                                   int corr_type,int *status) {
  corr_int_par cp;
  ccl_f1d_t *cl_spl = NULL;
  cl_spl = ccl_f1d_t_new(n_ell, ell, cls, cls[0], 0,
			 ccl_f1d_extrap_const,
			 ccl_f1d_extrap_logx_logy, status);
  if(cl_spl == NULL) {
    *status = CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_correlation.c: ccl_tracer_corr_bessel ran out of memory\n");
    return;
  }

  int ith, gslstatus;
  double result,eresult;
  gsl_function F;
  gsl_integration_workspace *w = NULL;
  int local_status;

#pragma omp parallel default(none) \
                     shared(cosmo, status, wtheta, n_ell, ell, cls, \
			    corr_type, cl_spl, theta, n_theta)	\
                     private(w, F, result, eresult, local_status, ith, \
			     gslstatus, cp)
  {
    local_status = *status;

    switch(corr_type) {
      case CCL_CORR_GG:
        cp.i_bessel = 0;
        break;
      case CCL_CORR_GL:
        cp.i_bessel = 2;
        break;
      case CCL_CORR_LP:
        cp.i_bessel = 0;
        break;
      case CCL_CORR_LM:
        cp.i_bessel = 4;
        break;
    }

    cp.cl_spl = cl_spl;

    w = gsl_integration_workspace_alloc(cosmo->gsl_params.N_ITERATION);

    if (w == NULL) {
      local_status = CCL_ERROR_MEMORY;
    }
    F.function = &corr_bessel_integrand;
    F.params = &cp;

    #pragma omp for schedule(dynamic)
    for(ith=0; ith < n_theta; ith++) {
      if (local_status == 0) {
        cp.th = theta[ith]*M_PI/180;
        //TODO: Split into intervals between first bessel zeros before integrating
        //This will help both speed and accuracy of the integral.
        gslstatus = gsl_integration_qag(&F, 0, cosmo->spline_params.ELL_MAX_CORR, 0,
                                        cosmo->gsl_params.INTEGRATION_EPSREL, cosmo->gsl_params.N_ITERATION,
                                        cosmo->gsl_params.INTEGRATION_GAUSS_KRONROD_POINTS,
                                        w, &result, &eresult);
        if(gslstatus != GSL_SUCCESS) {
          ccl_raise_gsl_warning(gslstatus, "ccl_correlation.c: ccl_tracer_corr_bessel():");
          local_status |= gslstatus;
        }
        wtheta[ith] = result/(2*M_PI);
      }
    }

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

    gsl_integration_workspace_free(w);
  }
  ccl_f1d_t_free(cl_spl);
}

/*--------ROUTINE: ccl_compute_legendre_polynomial ------
TASK: Compute input factor for ccl_tracer_corr_legendre
INPUT: tracer 1, tracer 2, i_bessel, theta array, n_theta, L_max, output Pl_theta
 */
static void ccl_compute_legendre_polynomial(int corr_type,double theta,int ell_max,double *Pl_theta)
{
  int i,j;
  double k=0;
  double cth=cos(theta*M_PI/180);

  //Initialize Pl_theta
  for (j=0;j<=ell_max;j++)
      Pl_theta[j]=0.;

  if(corr_type==CCL_CORR_GG) {
    gsl_sf_legendre_Pl_array(ell_max,cth,Pl_theta);
    for (j=0;j<=ell_max;j++)
      Pl_theta[j]*=(2*j+1);
  }
  else if(corr_type==CCL_CORR_GL) {
    for (j=2;j<=ell_max;j++) {//https://arxiv.org/pdf/1007.4809.pdf
      Pl_theta[j]=gsl_sf_legendre_Plm(j,2,cth);
      Pl_theta[j]*=(2*j+1.)/((j+0.)*(j+1.));
    }
  }
}

/*--------ROUTINE: ccl_tracer_corr_legendre ------
TASK: Compute correlation function via Legendre polynomials
INPUT: cosmology, number of theta bins, theta array, tracer 1, tracer 2, i_bessel, boolean
       for tapering, vector of tapering limits, correlation vector, angular_cl function.
 */
static void ccl_tracer_corr_legendre(ccl_cosmology *cosmo,
                                     int n_ell,double *ell,double *cls,
                                     int n_theta,double *theta,double *wtheta,
                                     int corr_type,int do_taper_cl,double *taper_cl_limits,
                                     int *status) {
  int i;
  double *l_arr = NULL, *cl_arr = NULL, *Pl_theta = NULL;
  ccl_f1d_t *cl_spl;

  if(corr_type==CCL_CORR_LM || corr_type==CCL_CORR_LP){
    *status=CCL_ERROR_NOT_IMPLEMENTED;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: CCL does not support full-sky xi+- calcuations.\nhttps://arxiv.org/abs/1702.05301 indicates flat-sky to be sufficient.\n");
  }

  if(*status==0) {
    l_arr=malloc(((int)(cosmo->spline_params.ELL_MAX_CORR)+1)*sizeof(double));
    if(l_arr==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_legendre ran out of memory\n");
    }
  }

  if(*status==0) {
    cl_arr=malloc(((int)(cosmo->spline_params.ELL_MAX_CORR)+1)*sizeof(double));
    if(cl_arr==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_legendre ran out of memory\n");
    }
  }

  if(*status==0) {
    //Interpolate input Cl into
    cl_spl=ccl_f1d_t_new(n_ell,ell,cls,cls[0],0,
			 ccl_f1d_extrap_const,
			 ccl_f1d_extrap_logx_logy, status);
    if(cl_spl==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_legendre ran out of memory\n");
    }
  }

  if(*status==0) {
    for(i=0;i<=(int)(cosmo->spline_params.ELL_MAX_CORR);i++) {
      double l=(double)i;
      l_arr[i]=l;
      cl_arr[i]=ccl_f1d_t_eval(cl_spl,l);
    }
    ccl_f1d_t_free(cl_spl);

    if (do_taper_cl)
      *status=taper_cl((int)(cosmo->spline_params.ELL_MAX_CORR)+1,l_arr,cl_arr,taper_cl_limits);
  }

  int local_status, i_L;
#pragma omp parallel default(none) \
                     shared(cosmo, theta, cl_arr, wtheta, n_theta, status, corr_type) \
                     private(Pl_theta, i, i_L, local_status)
  {
    Pl_theta = NULL;
    local_status = *status;

    if (local_status == 0) {
      Pl_theta = malloc(sizeof(double)*((int)(cosmo->spline_params.ELL_MAX_CORR)+1));
      if (Pl_theta == NULL) {
        local_status = CCL_ERROR_MEMORY;
      }
    }

    #pragma omp for schedule(dynamic)
    for (int i=0; i < n_theta; i++) {
      if (local_status == 0) {
        wtheta[i] = 0;
        ccl_compute_legendre_polynomial(corr_type, theta[i], (int)(cosmo->spline_params.ELL_MAX_CORR), Pl_theta);
        for (i_L=1; i_L < (int)(cosmo->spline_params.ELL_MAX_CORR); i_L+=1)
          wtheta[i] += cl_arr[i_L]*Pl_theta[i_L];
        wtheta[i] /= (M_PI*4);
      }
    }

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

    free(Pl_theta);
  }
  free(l_arr);
  free(cl_arr);
}

/*--------ROUTINE: ccl_tracer_corr ------
TASK: For a given tracer, get the correlation function. Do so by running
      ccl_angular_cls. If you already have Cls calculated, go to the next
      function to pass them directly.
INPUT: cosmology, number of theta values to evaluate = NL, theta vector,
       tracer 1, tracer 2, i_bessel, key for tapering, limits of tapering
       correlation function.
 */
void ccl_correlation(ccl_cosmology *cosmo,
                     int n_ell,double *ell,double *cls,
                     int n_theta,double *theta,double *wtheta,
                     int corr_type,int do_taper_cl,double *taper_cl_limits,int flag_method,
                     int *status) {
  switch(flag_method) {
  case CCL_CORR_FFTLOG :
    ccl_tracer_corr_fftlog(cosmo,n_ell,ell,cls,n_theta,theta,wtheta,corr_type,
                           do_taper_cl,taper_cl_limits,status);
    break;
  case CCL_CORR_LGNDRE :
    ccl_tracer_corr_legendre(cosmo,n_ell,ell,cls,n_theta,theta,wtheta,corr_type,
                             do_taper_cl,taper_cl_limits,status);
    break;
  case CCL_CORR_BESSEL :
    ccl_tracer_corr_bessel(cosmo,n_ell,ell,cls,n_theta,theta,wtheta,corr_type,status);
    break;
  default :
    *status=CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation. Unknown algorithm\n");
  }

}

/*--------ROUTINE: ccl_correlation_3d ------
TASK: Calculate the 3d-correlation function. Do so by using FFTLog.

INPUT: cosmology, scale factor a,
       number of r values, r values,
       key for tapering, limits of tapering

Correlation function result will be in array xi
 */

void ccl_correlation_3d(ccl_cosmology *cosmo, double a,
                        int n_r,double *r,double *xi,
                        int do_taper_pk,double *taper_pk_limits,
                        int *status) {
  int i,N_ARR;
  double *k_arr,*pk_arr,*r_arr,*xi_arr;

  if (!cosmo->computed_nonlin_power) {
    *status = CCL_ERROR_NONLIN_POWER_INIT;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_correlation.c: ccl_correlation_3d(): non-linear power spctrum has not been computed!");
    return;
  }

  //number of data points for k and pk array
  N_ARR=(int)(cosmo->spline_params.N_K_3DCOR*log10(cosmo->spline_params.K_MAX/cosmo->spline_params.K_MIN));

  k_arr=ccl_log_spacing(cosmo->spline_params.K_MIN,cosmo->spline_params.K_MAX,N_ARR);
  if(k_arr==NULL) {
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_3d ran out of memory\n");
    return;
  }

  pk_arr=malloc(N_ARR*sizeof(double));
  if(pk_arr==NULL) {
    free(k_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_3d ran out of memory\n");
    return;
  }

  for (i=0; i<N_ARR; i++){
    pk_arr[i] = ccl_nonlin_matter_power(cosmo, k_arr[i], a, status);
  }
  if (do_taper_pk)
    taper_cl(N_ARR,k_arr,pk_arr,taper_pk_limits);

  r_arr=malloc(sizeof(double)*N_ARR);
  if(r_arr==NULL) {
    free(k_arr);
    free(pk_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_3d ran out of memory\n");
    return;
  }
  xi_arr=malloc(sizeof(double)*N_ARR);
  if(xi_arr==NULL) {
    free(k_arr);
    free(pk_arr);
    free(r_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_3d ran out of memory\n");
    return;
  }

  for(i=0;i<N_ARR;i++)
    r_arr[i]=0;

  ccl_fftlog_ComputeXi3D(0, 0, 1, N_ARR, k_arr, &pk_arr, r_arr, &xi_arr, status);

  // Interpolate to output values of r
  ccl_f1d_t *xi_spl=ccl_f1d_t_new(N_ARR,r_arr,xi_arr,xi_arr[0],0,
				  ccl_f1d_extrap_const,
				  ccl_f1d_extrap_const, status);
  if (xi_spl == NULL) {
    free(k_arr);
    free(pk_arr);
    free(r_arr);
    free(xi_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_3d ran out of memory\n");
    return;
  }
  for(i=0;i<n_r;i++)
    xi[i]=ccl_f1d_t_eval(xi_spl,r[i]);
  ccl_f1d_t_free(xi_spl);

  free(k_arr);
  free(pk_arr);
  free(r_arr);
  free(xi_arr);

  return;
}

/*--------ROUTINE: ccl_correlation_multipole ------
TASK: Calculate multipole of the redshift space correlation function. Do so using FFTLog.

INPUT:  cosmology, scale factor a, beta (= growth rate / bias),
        multipole order l = 0, 2, or 4, number of s values, s values

Multipole function result will be in array xi
 */

void ccl_correlation_multipole(ccl_cosmology *cosmo, double a, double beta,
                               int l, int n_s, double *s, double *xi,
                               int *status) {
  int i, N_ARR;
  double *k_arr, *pk_arr, *s_arr, *xi_arr, *xi_arr0;

  if (!cosmo->computed_nonlin_power) {
    *status = CCL_ERROR_NONLIN_POWER_INIT;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_correlation.c: ccl_correlation_multipole(): non-linear power spctrum has not been computed!");
    return;
  }

  N_ARR = (int)(cosmo->spline_params.N_K_3DCOR * log10(cosmo->spline_params.K_MAX / cosmo->spline_params.K_MIN));

  k_arr = ccl_log_spacing(cosmo->spline_params.K_MIN, cosmo->spline_params.K_MAX, N_ARR);
  if (k_arr == NULL) {
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole ran out of memory\n");
    return;
  }

  pk_arr = malloc(N_ARR * sizeof(double));
  if (pk_arr == NULL) {
    free(k_arr);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole ran out of memory\n");
    return;
  }

  for (i = 0; i < N_ARR; i++)
    pk_arr[i] = ccl_nonlin_matter_power(cosmo, k_arr[i], a, status);

  s_arr = malloc(sizeof(double) * N_ARR);
  if (s_arr == NULL) {
    free(k_arr);
    free(pk_arr);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole ran out of memory\n");
    return;
  }
  xi_arr = malloc(sizeof(double) * N_ARR);
  if (xi_arr == NULL) {
    free(k_arr);
    free(pk_arr);
    free(s_arr);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole ran out of memory\n");
    return;
  }
  xi_arr0 = malloc(sizeof(double) * N_ARR);
  if (xi_arr0 == NULL) {
    free(k_arr);
    free(pk_arr);
    free(s_arr);
    free(xi_arr);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole ran out of memory\n");
    return;
  }

  for (i = 0; i < N_ARR; i++) s_arr[i] = 0;

  // Calculate multipoles

  if (l == 0) {
    ccl_fftlog_ComputeXi3D(0, 0, 1, N_ARR, k_arr, &pk_arr, s_arr, &xi_arr0, status);
    for (i = 0; i < N_ARR; i++)
      xi_arr[i] = (1. + 2. / 3 * beta + 1. / 5 * beta * beta) * xi_arr0[i];
  } else if (l == 2) {
    ccl_fftlog_ComputeXi3D(2, 0, 1, N_ARR, k_arr, &pk_arr, s_arr, &xi_arr0, status);
    for (i = 0; i < N_ARR; i++)
      xi_arr[i] = -(4. / 3 * beta + 4. / 7 * beta * beta) * xi_arr0[i];
  } else if (l == 4) {
    ccl_fftlog_ComputeXi3D(4, 0, 1, N_ARR, k_arr, &pk_arr, s_arr, &xi_arr0, status);
    for (i = 0; i < N_ARR; i++) xi_arr[i] = 8. / 35 * beta * beta * xi_arr0[i];
  } else {
    strcpy(cosmo->status_message, "unavailable value of l\n");
    return;
  }

  // Interpolate to output values of s
  ccl_f1d_t *xi_spl = ccl_f1d_t_new(N_ARR, s_arr, xi_arr, xi_arr[0], 0,
				    ccl_f1d_extrap_const,
				    ccl_f1d_extrap_const, status);
  if (xi_spl == NULL) {
    free(k_arr);
    free(pk_arr);
    free(s_arr);
    free(xi_arr);
    free(xi_arr0);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole ran out of memory\n");
  }
  for (i = 0; i < n_s; i++) xi[i] = ccl_f1d_t_eval(xi_spl,s[i]);
  ccl_f1d_t_free(xi_spl);

  free(k_arr);
  free(pk_arr);
  free(s_arr);
  free(xi_arr);
  free(xi_arr0);

  return;
}

/*--------ROUTINE: ccl_correlation_multipole_spline ------
TASK: Store multipoles of the redshift-space correlation in global splines

INPUT:  cosmology, scale factor a

Result is stored in cosmo->data.rsd_splines[]
 */

void ccl_correlation_multipole_spline(ccl_cosmology *cosmo, double a,
                                      int *status) {
  int i, N_ARR;
  double *k_arr, *pk_arr, *s_arr, *xi_arr, *xi_arr0, *xi_arr2, *xi_arr4;

  if (!cosmo->computed_nonlin_power) {
    *status = CCL_ERROR_NONLIN_POWER_INIT;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_correlation.c: ccl_correlation_multipole_spline(): non-linear power spctrum has not been computed!");
    return;
  }

  N_ARR = (int)(cosmo->spline_params.N_K_3DCOR * log10(cosmo->spline_params.K_MAX / cosmo->spline_params.K_MIN));

  k_arr = ccl_log_spacing(cosmo->spline_params.K_MIN, cosmo->spline_params.K_MAX, N_ARR);
  if (k_arr == NULL) {
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole_spline ran out of "
           "memory\n");
    return;
  }

  pk_arr = malloc(N_ARR * sizeof(double));
  if (pk_arr == NULL) {
    free(k_arr);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole_spline ran out of "
           "memory\n");
    return;
  }

  for (i = 0; i < N_ARR; i++)
    pk_arr[i] = ccl_nonlin_matter_power(cosmo, k_arr[i], a, status);

  s_arr = malloc(sizeof(double) * N_ARR);
  if (s_arr == NULL) {
    free(k_arr);
    free(pk_arr);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole_spline ran out of "
           "memory\n");
    return;
  }
  xi_arr = malloc(sizeof(double) * N_ARR);
  if (xi_arr == NULL) {
    free(k_arr);
    free(pk_arr);
    free(s_arr);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole_spline ran out of "
           "memory\n");
    return;
  }
  xi_arr0 = malloc(sizeof(double) * N_ARR);
  if (xi_arr0 == NULL) {
    free(k_arr);
    free(pk_arr);
    free(s_arr);
    free(xi_arr);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole_spline ran out of "
           "memory\n");
    return;
  }
  xi_arr2 = malloc(sizeof(double) * N_ARR);
  if (xi_arr2 == NULL) {
    free(k_arr);
    free(pk_arr);
    free(s_arr);
    free(xi_arr);
    free(xi_arr0);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole_spline ran out of "
           "memory\n");
    return;
  }
  xi_arr4 = malloc(sizeof(double) * N_ARR);
  if (xi_arr4 == NULL) {
    free(k_arr);
    free(pk_arr);
    free(s_arr);
    free(xi_arr);
    free(xi_arr0);
    free(xi_arr2);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole_spline ran out of "
           "memory\n");
    return;
  }

  for (i = 0; i < N_ARR; i++) s_arr[i] = 0;

  // Calculate multipoles
  ccl_fftlog_ComputeXi3D(0, 0, 1, N_ARR, k_arr, &pk_arr, s_arr, &xi_arr0, status);
  ccl_fftlog_ComputeXi3D(2, 0, 1, N_ARR, k_arr, &pk_arr, s_arr, &xi_arr2, status);
  ccl_fftlog_ComputeXi3D(4, 0, 1, N_ARR, k_arr, &pk_arr, s_arr, &xi_arr4, status);

  // free any memory that may have been allocated
  ccl_f1d_t_free(cosmo->data.rsd_splines[0]);
  ccl_f1d_t_free(cosmo->data.rsd_splines[1]);
  ccl_f1d_t_free(cosmo->data.rsd_splines[2]);
  cosmo->data.rsd_splines[0] = NULL;
  cosmo->data.rsd_splines[1] = NULL;
  cosmo->data.rsd_splines[1] = NULL;

  // Interpolate to output values of s
  cosmo->data.rsd_splines[0] = ccl_f1d_t_new(N_ARR, s_arr, xi_arr0, xi_arr0[0], 0,
					     ccl_f1d_extrap_const,
					     ccl_f1d_extrap_const, status);
  if (cosmo->data.rsd_splines[0] == NULL) {
    free(k_arr);
    free(pk_arr);
    free(s_arr);
    free(xi_arr);
    free(xi_arr0);
    free(xi_arr2);
    free(xi_arr4);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole_spline ran out of "
           "memory\n");
    return;
  }

  cosmo->data.rsd_splines[1] = ccl_f1d_t_new(N_ARR, s_arr, xi_arr2, xi_arr2[0], 0,
					     ccl_f1d_extrap_const,
					     ccl_f1d_extrap_const, status);
  if (cosmo->data.rsd_splines[1] == NULL) {
    free(k_arr);
    free(pk_arr);
    free(s_arr);
    free(xi_arr);
    free(xi_arr0);
    free(xi_arr2);
    free(xi_arr4);
    ccl_f1d_t_free(cosmo->data.rsd_splines[0]);
    cosmo->data.rsd_splines[0] = NULL;
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole_spline ran out of "
           "memory\n");
    return;
  }

  cosmo->data.rsd_splines[2] = ccl_f1d_t_new(N_ARR, s_arr, xi_arr4, xi_arr4[0], 0,
					     ccl_f1d_extrap_const,
					     ccl_f1d_extrap_const, status);
  if (cosmo->data.rsd_splines[2] == NULL) {
    free(k_arr);
    free(pk_arr);
    free(s_arr);
    free(xi_arr);
    free(xi_arr0);
    free(xi_arr2);
    free(xi_arr4);
    ccl_f1d_t_free(cosmo->data.rsd_splines[0]);
    cosmo->data.rsd_splines[0] = NULL;
    ccl_f1d_t_free(cosmo->data.rsd_splines[1]);
    cosmo->data.rsd_splines[1] = NULL;
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_multipole_spline ran out of "
           "memory\n");
    return;
  }

  // set the scale factor
  cosmo->data.rsd_splines_scalefactor = a;

  free(k_arr);
  free(pk_arr);
  free(s_arr);
  free(xi_arr);
  free(xi_arr0);
  free(xi_arr2);
  free(xi_arr4);

  return;
}

/*--------ROUTINE: ccl_correlation_3dRsd ------
TASK: Calculate the redshift-space correlation function.

INPUT:  cosmology, scale factor a, number of s values, s values,
        mu = cosine of galaxy separation angle w.r.t. line of sight,
        beta (= growth rate / bias), key for using spline

Correlation function result will be in array xi
 */

void ccl_correlation_3dRsd(ccl_cosmology *cosmo, double a, int n_s, double *s,
                           double mu, double beta, double *xi, int use_spline,
                           int *status) {
  int i;
  double *xi_arr0, *xi_arr2, *xi_arr4;

  if (!cosmo->computed_nonlin_power) {
    *status = CCL_ERROR_NONLIN_POWER_INIT;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_correlation.c: ccl_correlation_3dRsd(): non-linear power spctrum has not been computed!");
    return;
  }

  if (use_spline == 0) {
    xi_arr0 = malloc(sizeof(double) * n_s);
    if (xi_arr0 == NULL) {
      *status = CCL_ERROR_MEMORY;
      strcpy(cosmo->status_message,
             "ccl_correlation.c: ccl_correlation_3dRsd ran out of memory\n");
      return;
    }
    xi_arr2 = malloc(sizeof(double) * n_s);
    if (xi_arr2 == NULL) {
      free(xi_arr0);
      *status = CCL_ERROR_MEMORY;
      strcpy(cosmo->status_message,
             "ccl_correlation.c: ccl_correlation_3dRsd ran out of memory\n");
      return;
    }
    xi_arr4 = malloc(sizeof(double) * n_s);
    if (xi_arr4 == NULL) {
      free(xi_arr0);
      free(xi_arr2);
      *status = CCL_ERROR_MEMORY;
      strcpy(cosmo->status_message,
             "ccl_correlation.c: ccl_correlation_3dRsd ran out of memory\n");
      return;
    }

    ccl_correlation_multipole(cosmo, a, beta, 0, n_s, s, xi_arr0, status);
    ccl_correlation_multipole(cosmo, a, beta, 2, n_s, s, xi_arr2, status);
    ccl_correlation_multipole(cosmo, a, beta, 4, n_s, s, xi_arr4, status);
    for (i = 0; i < n_s; i++)
      xi[i] = xi_arr0[i] + xi_arr2[i] * gsl_sf_legendre_Pl(2, mu) +
              xi_arr4[i] * gsl_sf_legendre_Pl(4, mu);
    free(xi_arr0);
    free(xi_arr2);
    free(xi_arr4);

  } else {
    if ((cosmo->data.rsd_splines[0] == NULL) ||
        (cosmo->data.rsd_splines[1] == NULL) ||
        (cosmo->data.rsd_splines[2] == NULL) ||
        (cosmo->data.rsd_splines_scalefactor != a))
      ccl_correlation_multipole_spline(cosmo, a, status);

    for (i = 0; i < n_s; i++)
      xi[i] = (1. + 2. / 3 * beta + 1. / 5 * beta * beta) *
        ccl_f1d_t_eval(cosmo->data.rsd_splines[0],s[i]) -
        (4. / 3 * beta + 4. / 7 * beta * beta) *
        ccl_f1d_t_eval(cosmo->data.rsd_splines[1],s[i]) *
        gsl_sf_legendre_Pl(2, mu) +
        8. / 35 * beta * beta * ccl_f1d_t_eval(cosmo->data.rsd_splines[2],s[i]) *
        gsl_sf_legendre_Pl(4, mu);
  }

  return;
}

/*--------ROUTINE: ccl_correlation_3dRsd_avgmu ------
TASK: Calculate the average of redshift-space correlation function xi(s,mu) over mu at constant s

INPUT:  cosmology, scale factor a, number of s values, s values, beta (= growth rate / bias)

The result will be in array xi
*/

void ccl_correlation_3dRsd_avgmu(ccl_cosmology *cosmo, double a, int n_s, double *s,
                                 double beta, double *xi,
                                 int *status) {
// The average is just the l=0 multipole - the higher multiples inetegrate to zero.
  ccl_correlation_multipole(cosmo, a, beta, 0, n_s, s, xi, status);

  return;
}

/*--------ROUTINE: ccl_correlation_pi_sigma ------
TASK: Calculate the redshift-space correlation function using longitudinal and
      transverse coordinates pi and sigma.

INPUT:  cosmology, scale factor a, beta (= growth rate / bias),
        pi, number of sigma values, sigma values,
        key for using spline

Correlation function result will be in array xi
*/

void ccl_correlation_pi_sigma(ccl_cosmology *cosmo, double a, double beta,
                              double pi, int n_sig, double *sig, double *xi,
                              int use_spline, int *status) {
  int i;
  double *mu_arr, *s_arr, *xi_arr;

  if (!cosmo->computed_nonlin_power) {
    *status = CCL_ERROR_NONLIN_POWER_INIT;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_correlation.c: ccl_correlation_pi_sigma(): non-linear power spctrum has not been computed!");
    return;
  }

  mu_arr = malloc(sizeof(double) * n_sig);
  if (mu_arr == NULL) {
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_pi_sigma ran out of memory\n");
    return;
  }

  s_arr = malloc(sizeof(double) * n_sig);
  if (s_arr == NULL) {
    free(mu_arr);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_pi_sigma ran out of memory\n");
    return;
  }

  xi_arr = malloc(sizeof(double) * n_sig);
  if (xi_arr == NULL) {
    free(mu_arr);
    free(s_arr);
    *status = CCL_ERROR_MEMORY;
    strcpy(cosmo->status_message,
           "ccl_correlation.c: ccl_correlation_pi_sigma ran out of memory\n");
    return;
  }

  for (i = 0; i < n_sig; i++) {
    s_arr[i] = sqrt(pi * pi + sig[i] * sig[i]);
    mu_arr[i] = pi / s_arr[i];
  }

  for (i = 0; i < n_sig; i++) {
    ccl_correlation_3dRsd(cosmo, a, n_sig, s_arr, mu_arr[i], beta, xi_arr,
                          use_spline, status);
    xi[i] = xi_arr[i];
  }

  free(mu_arr);
  free(xi_arr);
  free(s_arr);

  return;
}

max_stars_repo_name: damonge/CCL

max_issues_repo_name: damonge/CCL

max_forks_repo_name: damonge/CCL


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

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_erf.h>

double
gsl_ran_gaussian_tail (const gsl_rng * r, const double a, const double sigma)
{
  /* Returns a gaussian random variable larger than a
   * This implementation does one-sided upper-tailed deviates.
   */

  double s = a / sigma;

  if (s < 1)
    {
      /* For small s, use a direct rejection method. The limit s < 1
         can be adjusted to optimise the overall efficiency */

      double x;

      do
        {
          x = gsl_ran_gaussian (r, 1.0);
        }
      while (x < s);
      return x * sigma;
    }
  else
    {
      /* Use the "supertail" deviates from the last two steps
       * of Marsaglia's rectangle-wedge-tail method, as described
       * in Knuth, v2, 3rd ed, pp 123-128.  (See also exercise 11, p139,
       * and the solution, p586.)
       */

      double u, v, x;

      do
        {
          u = gsl_rng_uniform (r);
          do
            {
              v = gsl_rng_uniform (r);
            }
          while (v == 0.0);
          x = sqrt (s * s - 2 * log (v));
        }
      while (x * u > s);
      return x * sigma;
    }
}

double
gsl_ran_gaussian_tail_pdf (const double x, const double a, const double sigma)
{
  if (x < a)
    {
      return 0;
    }
  else
    {
      double N, p;
      double u = x / sigma ;

      double f = gsl_sf_erfc (a / (sqrt (2.0) * sigma));

      N = 0.5 * f;

      p = (1 / (N * sqrt (2 * M_PI) * sigma)) * exp (-u * u / 2);

      return p;
    }
}

double
gsl_ran_ugaussian_tail (const gsl_rng * r, const double a)
{
  return gsl_ran_gaussian_tail (r, a, 1.0) ;
}

double
gsl_ran_ugaussian_tail_pdf (const double x, const double a)
{
  return gsl_ran_gaussian_tail_pdf (x, a, 1.0) ;
}

ext: c

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 "schedule.h"
#include <gsl/gsl-lite.hpp>
#include <range/v3/view/span.hpp>
#include <vector>

namespace angonoka::stun {
using ranges::span;
struct ScheduleParams;

/**
    Makespan estimator.

    Further reading:
    https://en.wikipedia.org/wiki/Makespan
*/
class Makespan {
public:
    /**
        Constructor.

        @param params       An instance of ScheduleParams
        @param tasks_count  Total number of tasks
        @param agents_count Total number of agents
    */
    Makespan(const ScheduleParams& params);

    Makespan(const Makespan& other);
    Makespan& operator=(const Makespan& other) noexcept;
    Makespan(Makespan&& other) noexcept;
    Makespan& operator=(Makespan&& other) noexcept;
    ~Makespan() noexcept;

    /**
        Calculate the makespan of a given scheduling configuration.

        @param schedule Scheduling configuration

        @return Makespan in seconds
    */
    float operator()(Schedule schedule) noexcept;

    /**
        Get the current ScheduleParams object.

        @return Schedule parameters.
    */
    [[nodiscard]] const ScheduleParams& params() const;

    /**
        Set the ScheduleParams object.

        @param params ScheduleParams object
    */
    void params(const ScheduleParams& params);

private:
    struct Impl;
    gsl::not_null<const ScheduleParams*> params_;
    std::vector<float> sum_buffer;
    span<float> task_done;
    span<float> work_done;
};

class RandomUtils;

/**
    Shuffle tasks and agents in-place.

    Randomly swaps two adjacent tasks within the schedule and
    reassigns an agent of a random task.
*/
class Mutator {
public:
    /**
        Mutator options.

        @var params Schedule parameters
        @var random Random utils.
    */
    struct Options {
        gsl::not_null<const ScheduleParams*> params;
        gsl::not_null<RandomUtils*> random;
    };

    /**
        Constructor.

        @param params An instance of ScheduleParams
        @param random An instance of RandomUtils
    */
    Mutator(const ScheduleParams& params, RandomUtils& random);
    Mutator(const Options& options);

    /**
        Mutates the scheduling configuration in-place.

        @param schedule Scheduling configuration
    */
    void operator()(MutSchedule schedule) const noexcept;

    /**
        Get current options.

        @return Options.
    */
    [[nodiscard]] Options options() const;

    /**
        Set options.

        @param options Options.
    */
    void options(const Options& options);

private:
    struct Impl;
    gsl::not_null<const ScheduleParams*> params;
    gsl::not_null<RandomUtils*> random;
};
} // namespace angonoka::stun

ext: h

max_stars_repo_name: coffee-lord/angonoka

max_issues_repo_name: coffee-lord/angonoka

max_forks_repo_name: coffee-lord/angonoka


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "parmt_invert.h"
#ifdef PARMT_USE_INTEL
#include <mkl_lapacke.h>
#include <mkl_cblas.h>
#else
#include <lapacke.h>
#include <cblas.h>
#endif
#include "parmt_utils.h"
#include "iscl/array/array.h"
#include "iscl/linalg/linalg.h"
#include "iscl/memory/memory.h"

/*!
 * @brief Inverts for a moment tensor at each location in the grid-search.
 *
 * @param[in] data     parMT data structure with observations and
 *                     corresponding green's functions at all locations
 *                     for all waveforms.
 * @param[in] ldm      Leading dimension of mts.  This must be at least 6.
 * @param[in] lcDev    If true then apply the deviatoric constraint.
 * @param[in] lcDev    Otherwise invert for all moment tensors terms.
 *
 * @param[out] mts     Moment tensors tabulated at each location.  This has
 *                     dimension [data.nlocs x ldm].
 * @param[out] phi     L2 objective function for each moment tensor.  This has
 *                     dimension [data.nlocs].
 *
 * @result 0 indicates success.
 *
 */
int parmt_invertLocationSearch(struct parmtData_struct data,
                               const int ldm,
                               const bool lcDev,
                               double *__restrict__ mts,
                               double *__restrict__ phi)
{
    double *G, *gxxAll, *gyyAll, *gzzAll, *gxyAll, *gxzAll, *gyzAll,
           *obsAll, *wts;
    //double m6[6], phiLoc;
    int *dataPtr, i1, ierr, ierr1, ierr2, iloc, iobs, k, mrows;
    ierr = 0;
    if (ldm < 6 || mts == NULL || phi == NULL)
    {
        if (ldm < 6){fprintf(stderr, "%s: ldm must be at least 6\n", __func__);}
        if (mts == NULL){fprintf(stderr, "%s: mts is NULL\n", __func__);}
        if (phi == NULL){fprintf(stderr, "%s: phi is NULL\n", __func__);}
        return -1;
    }
    G = NULL;
    wts = NULL;
    gxxAll = NULL;
    gyyAll = NULL;
    gzzAll = NULL;
    gxyAll = NULL;
    gxzAll = NULL;
    gyzAll = NULL;
    // 6 observations required for full MT, 5 for deviatoric constraint
    if (data.nobs < 6 - (int) lcDev)
    {
        fprintf(stderr, "%s: Insufficient number of observations \n",
                __func__);
        return -1;
    }
    // Figure out the workspace
    dataPtr = memory_calloc32i(data.nobs+1);
    for (iobs=0; iobs<data.nobs; iobs++)
    {
        dataPtr[iobs+1] = dataPtr[iobs] + data.data[iobs].npts;
    }
    if (dataPtr[data.nobs] < 6 - (int) lcDev)
    {
        fprintf(stderr, "%s: Insufficient number of points to invert with\n",
                __func__);
        goto ERROR;
    }
    obsAll = memory_calloc64f(dataPtr[data.nobs]);
    gxxAll = memory_calloc64f(dataPtr[data.nobs]);
    gyyAll = memory_calloc64f(dataPtr[data.nobs]);
    gzzAll = memory_calloc64f(dataPtr[data.nobs]);
    gxyAll = memory_calloc64f(dataPtr[data.nobs]);
    gxzAll = memory_calloc64f(dataPtr[data.nobs]);
    gyzAll = memory_calloc64f(dataPtr[data.nobs]);
    // Loop on locations
    for (iloc=0; iloc<data.nlocs; iloc++)
    {
        // Loop on observations
        ierr1 = 0;
        ierr2 = 0;
        for (iobs=0; iobs<data.nobs; iobs++)
        {
            k = iobs*data.nlocs + iloc;;
            G = parmt_utils_sacGrns2GrnsMat(data.sacGxx[k], data.sacGyy[k],
                                            data.sacGzz[k], data.sacGxy[k],
                                            data.sacGxz[k], data.sacGyz[k],
                                            &mrows, &ierr1);
            memory_free64f(&G);
            if (ierr1 != 0)
            {
                fprintf(stderr, "%s: Failed to make greens functions matrix\n",
                        __func__);
                ierr2 = ierr2 + 1;
            }
            if (mrows != dataPtr[iobs+1] - dataPtr[iobs])
            {
                fprintf(stderr, "%s: Internal error on mrows and npts\n",
                        __func__);
                ierr2 = ierr2 + 1;
                goto NEXT_OBS;
            }
            // Copy the observations and Green's functions
            i1 = dataPtr[iobs];
            cblas_dcopy(mrows, &G[0*mrows], 1, &gxxAll[i1], 1);
            cblas_dcopy(mrows, &G[1*mrows], 1, &gyyAll[i1], 1);
            cblas_dcopy(mrows, &G[2*mrows], 1, &gzzAll[i1], 1);
            cblas_dcopy(mrows, &G[3*mrows], 1, &gxyAll[i1], 1);
            cblas_dcopy(mrows, &G[4*mrows], 1, &gxzAll[i1], 1);
            cblas_dcopy(mrows, &G[5*mrows], 1, &gyzAll[i1], 1);
            cblas_dcopy(mrows, data.data[iobs].data, 1, &obsAll[i1], 1);
NEXT_OBS:;
            memory_free64f(&G);

        } // loop on observations
        // Invert it
        if (ierr2 != 0)
        {
            ierr1 = parmt_invertMT64f(data.nobs,
                                      lcDev, dataPtr,
                                      wts,
                                      gxxAll, gyyAll, gzzAll,
                                      gxyAll, gxzAll, gyzAll,
                                      obsAll, &phi[iloc], &mts[iloc*ldm]);
            //&phiLoc, m6);
            if (ierr1 != 0)
            {
                fprintf(stderr, "%s: Error inverting location %d\n",
                        __func__, iloc);
                ierr = ierr + 1;
            }
        }
        else
        {
            ierr = ierr + 1;
        }
    }
ERROR:;
    memory_free32i(&dataPtr);
    memory_free64f(&gxxAll);
    memory_free64f(&gyyAll);
    memory_free64f(&gzzAll);
    memory_free64f(&gxyAll);
    memory_free64f(&gxzAll);
    memory_free64f(&gyzAll);
    return 0;
}
//============================================================================//
/*!
 * @brief
 *
 * @param[in] nobs     number of waveforms (observations)
 * @param[in] lcDev    if true then force the solution to be purely deviatoric.
 *                     otherwise all moment tensor terms will be inverted for.
 * @param[in] dataPtr  points from iobs'th observation to start index of
 *                     observation [nobs+1]
 * @param[in] wts      if not NULL then these are the data weights for each
 *                     observations [nobs].
 * @param[in] gxxAll   all Gxx greens functions [dataPtr[nobs]]
 * @param[in] gyyAll   all Gyy greens functions [dataPtr[nobs]]
 * @param[in] gzzAll   all Gzz greens functions [dataPtr[nobs]]
 * @param[in] gxyAll   all Gxy greens functions [dataPtr[nobs]]
 * @param[in] gxzAll   all Gxz greens functions [dataPtr[nobs]]
 * @param[in] gyzAll   all Gyz greens functions [dataPtr[nobs]]
 * @param[in] obsAll   observations [dataPtr[nobs]]
 *
 * @param[out] phi     L2 difference between observed and estimate seismograms.
 * @param[out] m       on successful contains the moment tensor elements packed
 *                     \f$\{ m_{xx}, m_{yy}, m_{zz},
 *                           m_{xy}, m_{xz}, m_{yz} \}\f$.
 *
 * @result 0 indicates success
 *
 * @author Ben Baker
 *
 * @copyright ISTI distributed under Apache 2
 *
 */
int parmt_invertMT64f(const int nobs,
                      const bool lcDev,
                      const int *dataPtr,
                      const double *__restrict__ wts,
                      const double *__restrict__ gxxAll,
                      const double *__restrict__ gyyAll,
                      const double *__restrict__ gzzAll,
                      const double *__restrict__ gxyAll,
                      const double *__restrict__ gxzAll,
                      const double *__restrict__ gyzAll,
                      const double *__restrict__ obsAll,
                      double *phi, double *__restrict__ m)
{
    double *G, *b, *est, m6[6];
    int i1, i2, ierr, iobs, ldg, mrows, ncols, nploc;
    *phi = 0.0;
    ierr = 0;
    // Create the forward modeling matrix
    if (nobs < 5 || dataPtr == NULL ||
        gxxAll == NULL || gyyAll == NULL || gzzAll == NULL ||
        gxyAll == NULL || gxzAll == NULL || gyzAll == NULL ||
        m == NULL)
    {
        if (nobs < 5)
        {
            fprintf(stderr, "%s: At least 6 observations required\n", __func__);
        }
        if (dataPtr == NULL)
        {
            fprintf(stderr, "%s: Error dataPtr is NULL\n", __func__);
        }
        if (gxxAll == NULL)
        {
            fprintf(stderr, "%s: Error gxxAll is NULL\n", __func__);
        }
        if (gyyAll == NULL)
        {
            fprintf(stderr, "%s: Error gyyAll is NULL\n", __func__);
        }
        if (gzzAll == NULL)
        {
            fprintf(stderr, "%s: Error gzzAll is NULL\n", __func__);
        }
        if (gxyAll == NULL)
        {
            fprintf(stderr, "%s: Error gxyAll is NULL\n", __func__);
        }
        if (gxzAll == NULL)
        {
            fprintf(stderr, "%s: Error gxzAll is NULL\n", __func__);
        }
        if (gyzAll == NULL)
        {
            fprintf(stderr, "%s: Error gyzAll is NULL\n", __func__);
        }
        if (m == NULL){fprintf(stderr, "%s: Error m is NULL\n", __func__);}
        return -1;
    }
    memset(m, 0, 6*sizeof(double));
    ncols = 6;
    if (lcDev){ncols = 5;}
    // Get the data size
    mrows = 0;
    for (iobs=0; iobs<nobs; iobs++)
    {
        mrows = mrows + dataPtr[iobs-1] - dataPtr[iobs];
    }
    if (mrows < ncols)
    {
        fprintf(stderr, "%s: Error problem is undetermined\n", __func__);
        return -1;
    }
    ldg = mrows;
    G = memory_calloc64f(ldg*ncols);
    b = memory_calloc64f(mrows);
    est = memory_calloc64f(mrows);
    // Set the forward modeling matrix
    if (wts == NULL)
    {
        if (lcDev)
        {
            cblas_dcopy(mrows, gxxAll, 1, &G[0*ldg], 1);
            cblas_daxpy(mrows, -1.0, gzzAll, 1, &G[0*ldg], 1); // mxx - mzz
            cblas_dcopy(mrows, gyyAll, 1, &G[1*ldg], 1);
            cblas_daxpy(mrows, -1.0, gzzAll, 1, &G[1*ldg], 1); // myy - mzz
            cblas_dcopy(mrows, gxyAll, 1, &G[2*ldg], 1);
            cblas_dcopy(mrows, gxzAll, 1, &G[3*ldg], 1);
            cblas_dcopy(mrows, gyzAll, 1, &G[4*ldg], 1);
        }
        else
        {
            cblas_dcopy(mrows, gxxAll, 1, &G[0*ldg], 1);
            cblas_dcopy(mrows, gyyAll, 1, &G[1*ldg], 1);
            cblas_dcopy(mrows, gyyAll, 1, &G[2*ldg], 1);
            cblas_dcopy(mrows, gxyAll, 1, &G[3*ldg], 1);
            cblas_dcopy(mrows, gxzAll, 1, &G[4*ldg], 1);
            cblas_dcopy(mrows, gyzAll, 1, &G[5*ldg], 1);
        }
        cblas_dcopy(mrows, obsAll, 1, b, 1);
    }
    else
    {
        if (lcDev)
        {
            for (iobs=0; iobs<nobs; iobs++)
            {
                i1 = dataPtr[iobs];
                i2 = dataPtr[iobs];
                nploc = i2 - i1;
                cblas_daxpy(nploc, wts[iobs], &gxxAll[i1], 1, &G[0*ldg+i1], 1);
                cblas_daxpy(nploc,-wts[iobs], &gzzAll[i1], 1, &G[0*ldg+i1], 1);
                cblas_daxpy(nploc, wts[iobs], &gyyAll[i1], 1, &G[1*ldg+i1], 1);
                cblas_daxpy(nploc,-wts[iobs], &gzzAll[i1], 1, &G[1*ldg+i1], 1);
                cblas_daxpy(nploc, wts[iobs], &gxyAll[i1], 1, &G[2*ldg+i1], 1);
                cblas_daxpy(nploc, wts[iobs], &gxzAll[i1], 1, &G[3*ldg+i1], 1);
                cblas_daxpy(nploc, wts[iobs], &gyzAll[i1], 1, &G[4*ldg+i1], 1);
                cblas_daxpy(nploc, wts[iobs], &obsAll[i1], 1, &b[i1], 1);
            }
        }
        else
        {
            for (iobs=0; iobs<nobs; iobs++)
            {
                i1 = dataPtr[iobs];
                i2 = dataPtr[iobs];
                nploc = i2 - i1;
                cblas_daxpy(nploc, wts[iobs], &gxxAll[i1], 1, &G[0*ldg+i1], 1);
                cblas_daxpy(nploc, wts[iobs], &gyyAll[i1], 1, &G[1*ldg+i1], 1);
                cblas_daxpy(nploc, wts[iobs], &gzzAll[i1], 1, &G[2*ldg+i1], 1);
                cblas_daxpy(nploc, wts[iobs], &gxyAll[i1], 1, &G[3*ldg+i1], 1);
                cblas_daxpy(nploc, wts[iobs], &gxzAll[i1], 1, &G[4*ldg+i1], 1);
                cblas_daxpy(nploc, wts[iobs], &gyzAll[i1], 1, &G[5*ldg+i1], 1);
                cblas_daxpy(nploc, wts[iobs], &obsAll[i1], 1, &b[i1], 1);
            }
        }
    }
    // Solve \tilde{G} m = \tilde{b}
    ierr = linalg_lstsq_qr64f_work(LAPACK_COL_MAJOR,
                                   mrows, ncols, 1,
                                   false, G, b,
                                   m6, NULL);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error solving Gm=u\n", __func__);
        ierr = 1;
        goto END;
    }
    if (lcDev)
    {
        m[0] = m6[0];
        m[1] = m6[1];
        m[2] =-(m6[0] + m6[1]);
        m[3] = m6[2];
        m[4] = m6[3];
        m[5] = m6[4];
    }
    else
    {
        cblas_dcopy(6, m6, 1, m, 1);
    }
    // Apply u = \tilde{G} m
    cblas_dgemv(CblasColMajor, CblasNoTrans, mrows, ncols,
                1.0, G, ncols, m6, 1, 0.0, est, 1);
    // Compute the least-squares residual
    *phi = array_normDiff64f(mrows, est, b, TWO_NORM, 2.0, &ierr);
END:;
    // Release memory
    memory_free64f(&G);
    memory_free64f(&b);
    memory_free64f(&est);
    return ierr;
}


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2018 Couchbase, Inc.
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 */
#pragma once

#include <array>
#include <gsl/gsl>
#include <optional>
#include <type_traits>

namespace couchbase::protocol
{

/**
 * Helper code for encode and decode of LEB128 values.
 * - mcbp encodes collection-ID as an unsigned LEB128
 * - see https://en.wikipedia.org/wiki/LEB128
 */

struct Leb128NoThrow {
};

/**
 * decode_unsigned_leb128 returns the decoded T and a const_byte_buffer
 * initialised with the data following the leb128 data. This form of the decode
 * does not throw for invalid input and the caller should always check
 * second.data() for success or error (see returns info).
 *
 * @param buf buffer containing a leb128 encoded value (of size T)
 * @returns On error a std::pair where first is set to 0 and second is nullptr/0
 *          const_byte_buffer. On success a std::pair where first is the decoded
 *          value and second is a buffer initialised with the data following the
 *          leb128 data.
 */
template<class T>
typename std::enable_if_t<std::is_unsigned_v<T>, std::pair<T, std::string_view>>
decode_unsigned_leb128(std::string_view buf, struct Leb128NoThrow /* unused */)
{
    T rv = static_cast<uint8_t>(buf[0]) & 0x7fULL;
    size_t end = 0;
    if ((static_cast<uint8_t>(buf[0]) & 0x80ULL) == 0x80ULL) {
        T shift = 7;
        // shift in the remaining data
        for (end = 1; end < buf.size(); end++) {
            rv |= (static_cast<uint8_t>(buf[end]) & 0x7fULL) << shift;
            if ((static_cast<uint8_t>(buf[end]) & 0x80ULL) == 0) {
                break; // no more
            }
            shift += 7;
        }

        // We should of stopped for a stop byte, not the end of the buffer
        if (end == buf.size()) {
            return { 0, std::string_view{} };
        }
    }
    // Return the decoded value and a buffer for any remaining data
    return { rv, std::string_view{ buf.data() + end + 1, buf.size() - (end + 1) } };
}

/**
 * decode_unsigned_leb128 returns the decoded T and a const_byte_buffer
 * initialised with the data following the leb128 data. This form of the decode
 * throws for invalid input.
 *
 * @param buf buffer containing a leb128 encoded value (of size T)
 * @returns std::pair first is the decoded value and second a buffer for the
 *          remaining data (size will be 0 for no more data)
 * @throws std::invalid_argument if buf[0] does not encode a leb128 value with
 *         a stop byte.
 */
template<class T>
typename std::enable_if_t<std::is_unsigned_v<T>, std::pair<T, std::string_view>>
decode_unsigned_leb128(std::string_view buf)
{
    if (buf.size() > 0) {
        auto rv = decode_unsigned_leb128<T>(buf, Leb128NoThrow());
        if (rv.second.data()) {
            return rv;
        }
    }
    throw std::invalid_argument("decode_unsigned_leb128: invalid buf size:" + std::to_string(buf.size()));
}

/**
 * @return a buffer to the data after the leb128 prefix
 */
template<class T>
typename std::enable_if_t<std::is_unsigned_v<T>, std::string_view>
skip_unsigned_leb128(std::string_view buf)
{
    return decode_unsigned_leb128<T>(buf).second;
}

// Empty, non specialised version of the decoder class
template<class T, class Enable = void>
class unsigned_leb128
{
};

/**
 * For encoding a unsigned T leb128, class constructs from a T value and
 * provides a const_byte_buffer for access to the encoded
 */
template<class T>
class unsigned_leb128<T, typename std::enable_if_t<std::is_unsigned_v<T>>>
{
  public:
    explicit unsigned_leb128(T in)
    {
        while (in > 0) {
            auto byte = gsl::narrow_cast<uint8_t>(in & 0x7fULL);
            in >>= 7;

            // In has more data?
            if (in > 0) {
                byte |= 0x80;
                encodedData[encodedSize - 1U] = byte;
                // Increase the size
                encodedSize++;
            } else {
                encodedData[encodedSize - 1U] = byte;
            }
        }
    }

    [[nodiscard]] std::string get() const
    {
        return { begin(), end() };
    }

    [[nodiscard]] const uint8_t* begin() const
    {
        return encodedData.data();
    }

    [[nodiscard]] const uint8_t* end() const
    {
        return encodedData.data() + encodedSize;
    }

    [[nodiscard]] const uint8_t* data() const
    {
        return encodedData.data();
    }

    [[nodiscard]] size_t size() const
    {
        return encodedSize;
    }

    constexpr static size_t getMaxSize()
    {
        return maxSize;
    }

  private:
    // Larger T may need a larger array
    static_assert(sizeof(T) <= 8, "Class is only valid for uint 8/16/64");

    // value is large enough to store ~0 as leb128
    static constexpr size_t maxSize = sizeof(T) + (((sizeof(T) + 1) / 8) + 1);
    std::array<uint8_t, maxSize> encodedData{};
    uint8_t encodedSize{ 1 };
};

} // namespace couchbase::protocol

max_stars_repo_name: avsej/couchbase-ruby-client

max_issues_repo_name: avsej/couchbase-ruby-client

max_forks_repo_name: avsej/couchbase-ruby-client


text:
#pragma once
/****************************************************************************
 *
 *  nc.h
 *      ($\nativecoin-cpp\src)
 *
 *  by icedac
 *
 ***/
#ifndef ___NATIVECOIN__NATIVECOIN__NC_H_
#define ___NATIVECOIN__NATIVECOIN__NC_H_

#include <string>
#include <vector>
#include <unordered_map>
#include <set>
#include <iostream>
#include <sstream>

#define GSL_THROW_ON_CONTRACT_VIOLATION
#include <gsl/gsl>

#include <utility>
#include <stdio.h>
#include <tuple>
#include <iomanip>

#include "platform.h"

namespace nc {
    typedef unsigned char byte;
    typedef std::int32_t int32;
    typedef std::int64_t int64;
    typedef std::uint32_t uint32;
    typedef std::uint64_t uint64;
    using string = std::string;
    template < typename T >
    using vector = std::vector<T>;
    template < typename T >
    using set = std::set<T>;
    using stringstream = std::stringstream;
    using ostringstream = std::ostringstream;

    template < typename K, typename T >
    using unordered_map = std::unordered_map< K, T>;

    typedef uint32 difficulty_t;
    typedef uint64 nonce_t;
    typedef uint64 timestamp_t;
    typedef uint32 blockindex_t;
}

namespace nc {

    class coin {
    public:
        static const auto SATOSHI_PER_COIN = 1000000;
        explicit inline coin() noexcept : satoshi(0) {}
        explicit inline coin(int32 amount) noexcept : satoshi(amount*SATOSHI_PER_COIN) {}

        inline coin(const coin& c) = default;
        inline coin& operator = (const coin& rhs) = default;
        inline coin(coin&& c) = default;
        inline coin& operator = (coin&& rhs) = default;

        inline operator int32() const {
            return get();
        }

        inline coin& operator += (const coin& rhs) {
            satoshi += rhs.get_as_satoshi();
            return *this;
        }

        inline coin& operator -= (const coin& rhs) {
            satoshi -= rhs.get_as_satoshi();
            return *this;
        }

        inline coin operator + (const coin& rhs) {
            return coin::from_satoshi(this->get_as_satoshi() + rhs.get_as_satoshi());
        }
        inline coin operator - (const coin& rhs) {
            return coin::from_satoshi(this->get_as_satoshi() - rhs.get_as_satoshi());
        }

        inline bool operator == (const coin& rhs) const {
            return this->satoshi == rhs.satoshi;
        }
        inline bool operator != (const coin& rhs) const {
            return !(*this == rhs);
        }

        inline int32 get() const {
            return (int32)(satoshi / SATOSHI_PER_COIN);
        }

        inline int64 get_as_satoshi() const { return satoshi; }

        inline int64 get_satoshi() const {
            return satoshi % SATOSHI_PER_COIN;
        }

        inline string as_string() const {

            int64 satoshi = get_satoshi();

            ostringstream ss;
            ss << get() << ".";
            if (satoshi)
                ss << std::setfill('0') << std::setw(6) << satoshi;
            else
                ss << "00";

            return ss.str();
        }

        static inline coin from_satoshi(int64 satoshi) {
            coin c{ 0 };
            c.satoshi = satoshi;
            return c;
        }

    private:
        int64		satoshi;
    };
}

#endif// ___NATIVECOIN__NATIVECOIN__NC_H_

max_stars_repo_name: icedac/nativecoin-cpp

max_issues_repo_name: icedac/nativecoin-cpp

max_forks_repo_name: icedac/nativecoin-cpp


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C code for structures representing a waveform as a list of modes in amplitude/phase form.
 *
 */

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"

/* Call a function that references an object */
double ObjectFunctionCall(const ObjectFunction* this,double f){return this->function(this->object,f);};

/**************************************************************/
/* Functions computing the max and min between two int */
int max (int a, int b) { return a > b ? a : b; }
int min (int a, int b) { return a < b ? a : b; }

/************** GSL error handling and I/O ********************/

/* GSL error handler */
void Err_Handler(const char *reason, const char *file, int line, int gsl_errno) {
  printf("gsl: %s:%d: %s - %d\n", file, line, reason, gsl_errno);
  exit(1);
}

/* Functions to read binary data from files */
int Read_Vector(const char dir[], const char fname[], gsl_vector *v) {
  char *path=malloc(strlen(dir)+64);
  sprintf(path,"%s/%s", dir, fname);
  FILE *f = fopen(path, "rb");
  if (!f) {
    fprintf(stderr, "Error reading data from %s\n", path);
    free(path);
    return(FAILURE);
  }
  int ret = gsl_vector_fread(f, v);
  if (ret != 0) {
    fprintf(stderr, "Error reading data from %s.\n",path);
    free(path);
    return(FAILURE);
  }
  fclose(f);
  free(path);
  return(SUCCESS);
}
int Read_Matrix(const char dir[], const char fname[], gsl_matrix *m) {
  char *path=malloc(strlen(dir)+64);
  sprintf(path,"%s/%s", dir, fname);
  FILE *f = fopen(path, "rb");
  if (!f) {
    fprintf(stderr, "Error reading data from %s\n", path);
    free(path);
    return(FAILURE);
  }
  int ret = gsl_matrix_fread(f, m);
  if (ret != 0) {
    fprintf(stderr, "Error reading data from %s\n", path);
    free(path);
    return(FAILURE);
  }
  fclose(f);
  free(path);
  return(SUCCESS);
}
/* Functions to read text data from files */
int Read_Text_Vector(const char dir[], const char fname[], gsl_vector *v) {
  char *path=malloc(strlen(dir)+64);
  sprintf(path,"%s/%s", dir, fname);
  FILE *f = fopen(path, "rb");
  if (!f) {
    fprintf(stderr, "Error reading data from %s\n", path);
    free(path);
    return(FAILURE);
  }
  int ret = gsl_vector_fscanf(f, v);
  if (ret != 0) {
    fprintf(stderr, "Error reading data from %s.\n",path);
    free(path);
    return(FAILURE);
  }
  fclose(f);
  free(path);
  return(SUCCESS);
}
int Read_Text_Matrix(const char dir[], const char fname[], gsl_matrix *m) {
  char *path=malloc(strlen(dir)+64);
  sprintf(path,"%s/%s", dir, fname);
  FILE *f = fopen(path, "rb");
  if (!f) {
    fprintf(stderr, "Error reading data from %s\n", path);
    free(path);
    return(FAILURE);
  }
  int ret = gsl_matrix_fscanf(f, m);
  if (ret != 0) {
    fprintf(stderr, "Error reading data from %s.\n",path);
    free(path);
    return(FAILURE);
  }
  fclose(f);
  free(path);
  return(SUCCESS);
}
/* Functions to write data to files */
int Write_Vector(const char dir[], const char fname[], gsl_vector *v) {
  char *path=malloc(strlen(dir)+64);
  sprintf(path,"%s/%s", dir, fname);
  FILE *f = fopen(path, "w");
  if (!f) {
    fprintf(stderr, "Error writing data to %s\n", path);
    free(path);
    return(FAILURE);
  }
  int ret = gsl_vector_fwrite(f, v);
  if (ret != 0) {
    fprintf(stderr, "Error writing data to %s\n",path);
    free(path);
    return(FAILURE);
  }
  fclose(f);
  free(path);
  return(SUCCESS);
}
int Write_Matrix(const char dir[], const char fname[], gsl_matrix *m) {
  char *path=malloc(strlen(dir)+64);

  sprintf(path,"%s/%s", dir, fname);
  FILE *f = fopen(path, "w");
  if (!f) {
    fprintf(stderr, "Error writing data to %s\n", path);
    free(path);
    return(FAILURE);
  }
  int ret = gsl_matrix_fwrite(f, m);
  if (ret != 0) {
    fprintf(stderr, "Error writing data to %s\n", path);
    free(path);
    return(FAILURE);
  }
  fclose(f);
  free(path);
  return(SUCCESS);
}
/* Functions to write text data to files */
int Write_Text_Vector(const char dir[], const char fname[], gsl_vector *v) {
  char *path=malloc(strlen(dir)+64);
  sprintf(path,"%s/%s", dir, fname);
  FILE *f = fopen(path, "w");
  if (!f) {
    fprintf(stderr, "Error writing data to %s\n", path);
    free(path);
    return(FAILURE);
  }
  int ret = gsl_vector_fprintf(f, v, "%.16e");
  if (ret != 0) {
    fprintf(stderr, "Error writing data to %s\n",path);
    free(path);
    return(FAILURE);
  }
  fclose(f);
  free(path);
  return(SUCCESS);
}
int Write_Text_Matrix(const char dir[], const char fname[], gsl_matrix *m) {
  char *path=malloc(strlen(dir)+64);
  int ret = 0;

  sprintf(path,"%s/%s", dir, fname);
  FILE *f = fopen(path, "w");
  if (!f) {
    fprintf(stderr, "Error writing data to %s\n",path);
    free(path);
    return(FAILURE);
  }
  int N = (int) m->size1;
  int M = (int) m->size2;
  for(int i=0; i<N; i++){
    for(int j=0; j<M; j++){
      ret |= (fprintf(f, "%.16e ", gsl_matrix_get(m, i, j)) < 0);
    }
    if(i < N-1) ret |= (fprintf(f, "\n") < 0);
  }
  if (ret != 0) {
    fprintf(stderr, "Error writing data to %s\n",path);
    free(path);
    return(FAILURE);
  }
  fclose(f);
  free(path);
  return(SUCCESS);
}

/******** Functions to initialize and clean up CAmpPhaseFrequencySeries structure ********/
void CAmpPhaseFrequencySeries_Init(CAmpPhaseFrequencySeries **freqseries, const int n) {
  if(!freqseries) exit(1);
  /* Create storage for structures */
  if(!*freqseries) *freqseries=malloc(sizeof(CAmpPhaseFrequencySeries));
  else
  {
    CAmpPhaseFrequencySeries_Cleanup(*freqseries);
    *freqseries=malloc(sizeof(CAmpPhaseFrequencySeries));
  }
  gsl_set_error_handler(&Err_Handler);
  (*freqseries)->freq = gsl_vector_alloc(n);
  (*freqseries)->amp_real = gsl_vector_alloc(n);
  (*freqseries)->amp_imag = gsl_vector_alloc(n);
  (*freqseries)->phase = gsl_vector_alloc(n);
}
void CAmpPhaseFrequencySeries_Cleanup(CAmpPhaseFrequencySeries *freqseries) {
  if(freqseries->freq) gsl_vector_free(freqseries->freq);
  if(freqseries->amp_real) gsl_vector_free(freqseries->amp_real);
  if(freqseries->amp_imag) gsl_vector_free(freqseries->amp_imag);
  if(freqseries->phase) gsl_vector_free(freqseries->phase);
  free(freqseries);
}

/******** Functions to initialize and clean up CAmpPhaseSpline structure ********/
void CAmpPhaseSpline_Init(CAmpPhaseSpline **splines, const int n) {
  if(!splines) exit(1);
  /* Create storage for structures */
  if(!*splines) *splines=malloc(sizeof(CAmpPhaseSpline));
  else
  {
    CAmpPhaseSpline_Cleanup(*splines);
  }
  gsl_set_error_handler(&Err_Handler);
  (*splines)->spline_amp_real = gsl_matrix_alloc(n, 5);
  (*splines)->spline_amp_imag = gsl_matrix_alloc(n, 5);
  (*splines)->quadspline_phase = gsl_matrix_alloc(n, 4);
}
void CAmpPhaseSpline_Cleanup(CAmpPhaseSpline *splines) {
  if(splines->spline_amp_real) gsl_matrix_free(splines->spline_amp_real);
  if(splines->spline_amp_imag) gsl_matrix_free(splines->spline_amp_imag);
  if(splines->quadspline_phase) gsl_matrix_free(splines->quadspline_phase);
  free(splines);
}

/******** Functions to initialize and clean up CAmpPhaseGSLSpline structure ********/
void CAmpPhaseGSLSpline_Init(CAmpPhaseGSLSpline **splines, const int n) {
  if(!splines) exit(1);
  /* Create storage for structures */
  if(!*splines) *splines=malloc(sizeof(CAmpPhaseGSLSpline));
  else
  {
    CAmpPhaseGSLSpline_Cleanup(*splines);
  }
  gsl_set_error_handler(&Err_Handler);
  /* Note: for freq we won't copy the vector but simply copy the pointer to the existing one, so we don't allocate anything here */
  (*splines)->spline_amp_real = gsl_spline_alloc(gsl_interp_cspline, n);
  (*splines)->spline_amp_imag = gsl_spline_alloc(gsl_interp_cspline, n);
  (*splines)->spline_phase = gsl_spline_alloc(gsl_interp_cspline, n);
  (*splines)->accel_amp_real = gsl_interp_accel_alloc();
  (*splines)->accel_amp_imag = gsl_interp_accel_alloc();
  (*splines)->accel_phase = gsl_interp_accel_alloc();
}
void CAmpPhaseGSLSpline_Cleanup(CAmpPhaseGSLSpline *splines) {
  if(splines->spline_amp_real) gsl_spline_free(splines->spline_amp_real);
  if(splines->spline_amp_imag) gsl_spline_free(splines->spline_amp_imag);
  if(splines->spline_phase) gsl_spline_free(splines->spline_phase);
  if(splines->accel_amp_real) gsl_interp_accel_free(splines->accel_amp_real);
  if(splines->accel_amp_imag) gsl_interp_accel_free(splines->accel_amp_imag);
  if(splines->accel_phase) gsl_interp_accel_free(splines->accel_phase);
  free(splines);
}

/******** Functions to initialize and clean up ReImFrequencySeries structure ********/
void ReImFrequencySeries_Init(ReImFrequencySeries **freqseries, const int n) {
  if(!freqseries) exit(1);
  /* Create storage for structures */
  if(!*freqseries) *freqseries=malloc(sizeof(ReImFrequencySeries));
  else
  {
    ReImFrequencySeries_Cleanup(*freqseries);
  }
  gsl_set_error_handler(&Err_Handler);
  (*freqseries)->freq = gsl_vector_alloc(n);
  (*freqseries)->h_real = gsl_vector_alloc(n);
  (*freqseries)->h_imag = gsl_vector_alloc(n);
}
void ReImFrequencySeries_Cleanup(ReImFrequencySeries *freqseries) {
  if(freqseries->freq) gsl_vector_free(freqseries->freq);
  if(freqseries->h_real) gsl_vector_free(freqseries->h_real);
  if(freqseries->h_imag) gsl_vector_free(freqseries->h_imag);
  free(freqseries);
}

/******** Functions to initialize and clean up ReImTimeSeries structure ********/
void ReImTimeSeries_Init(ReImTimeSeries **timeseries, const int n) {
  if(!timeseries) exit(1);
  /* Create storage for structures */
  if(!*timeseries) *timeseries=malloc(sizeof(ReImTimeSeries));
  else
  {
    ReImTimeSeries_Cleanup(*timeseries);
  }
  gsl_set_error_handler(&Err_Handler);
  (*timeseries)->times = gsl_vector_alloc(n);
  (*timeseries)->h_real = gsl_vector_alloc(n);
  (*timeseries)->h_imag = gsl_vector_alloc(n);
}
void ReImTimeSeries_Cleanup(ReImTimeSeries *timeseries) {
  if(timeseries->times) gsl_vector_free(timeseries->times);
  if(timeseries->h_real) gsl_vector_free(timeseries->h_real);
  if(timeseries->h_imag) gsl_vector_free(timeseries->h_imag);
  free(timeseries);
}
/******** Functions to initialize and clean up AmpPhaseTimeSeries structure ********/
void AmpPhaseTimeSeries_Init(AmpPhaseTimeSeries **timeseries, const int n) {
  if(!timeseries) exit(1);
  /* Create storage for structures */
  if(!*timeseries) *timeseries=malloc(sizeof(AmpPhaseTimeSeries));
  else
  {
    AmpPhaseTimeSeries_Cleanup(*timeseries);
  }
  gsl_set_error_handler(&Err_Handler);
  (*timeseries)->times = gsl_vector_alloc(n);
  (*timeseries)->h_amp = gsl_vector_alloc(n);
  (*timeseries)->h_phase = gsl_vector_alloc(n);
}
void AmpPhaseTimeSeries_Cleanup(AmpPhaseTimeSeries *timeseries) {
  if(timeseries->times) gsl_vector_free(timeseries->times);
  if(timeseries->h_amp) gsl_vector_free(timeseries->h_amp);
  if(timeseries->h_phase) gsl_vector_free(timeseries->h_phase);
  free(timeseries);
}

/******** Functions to initialize and clean up RealTimeSeries structure ********/
void RealTimeSeries_Init(RealTimeSeries **timeseries, const int n) {
  if(!timeseries) exit(1);
  /* Create storage for structures */
  if(!*timeseries) *timeseries=malloc(sizeof(RealTimeSeries));
  else
  {
    RealTimeSeries_Cleanup(*timeseries);
  }
  gsl_set_error_handler(&Err_Handler);
  (*timeseries)->times = gsl_vector_alloc(n);
  (*timeseries)->h = gsl_vector_alloc(n);
}
void RealTimeSeries_Cleanup(RealTimeSeries *timeseries) {
  if(timeseries->times) gsl_vector_free(timeseries->times);
  if(timeseries->h) gsl_vector_free(timeseries->h);
  free(timeseries);
}

/***************** Functions for the ListmodesCAmpPhaseFrequencySeries structure ****************/
ListmodesCAmpPhaseFrequencySeries* ListmodesCAmpPhaseFrequencySeries_AddModeNoCopy(
	   ListmodesCAmpPhaseFrequencySeries* appended,  /* List structure to prepend to */
	   CAmpPhaseFrequencySeries* freqseries,  /* data to contain */
	   int l, /* major mode number */
	   int m  /* minor mode number */)
{
    ListmodesCAmpPhaseFrequencySeries* list;
    /* Check if the node with this mode already exists */
    list = appended;
    while( list ){
      if( l == list->l && m == list->m ){
	break;
      }
      list = list->next;
    }
    if( list ){ /* We don't allow for the case where the mode already exists in the list*/
      printf("Error: Tried to add an already existing mode to a ListmodesCAmpPhaseFrequencySeries ");
      return(NULL);
    } else { /* In that case, we do NOT COPY the input interpolated data, which therefore can't be
		used anywhere else; this will be acceptable as these operations will only be done
		when interpolating the initialization data */
      list = malloc( sizeof(ListmodesCAmpPhaseFrequencySeries) );
    }
    list->l = l;
    list->m = m;
    if( freqseries ){
      list->freqseries = freqseries;
    } else {
      list->freqseries = NULL;
    }
    if( appended ){
      list->next = appended;
    } else {
        list->next = NULL;
    }
    return list;
}
/* Get the element of a ListmodesCAmpPhaseFrequencySeries with a given index */
ListmodesCAmpPhaseFrequencySeries* ListmodesCAmpPhaseFrequencySeries_GetMode(
	   ListmodesCAmpPhaseFrequencySeries* const list,  /* List structure to get a particular mode from */
	   int l, /*< major mode number */
	   int m  /*< minor mode number */ )
{
    if( !list ) return NULL;

    ListmodesCAmpPhaseFrequencySeries *itr = list;
    while( itr->l != l || itr->m != m ){
        itr = itr->next;
        if( !itr ) return NULL;
    }
    return itr; /* The element returned is itself a pointer to a ListmodesCAmpPhaseFrequencySeries */
}
void ListmodesCAmpPhaseFrequencySeries_Destroy(
	   ListmodesCAmpPhaseFrequencySeries* list  /* List structure to destroy; notice that the data is destroyed too */
)
{
  ListmodesCAmpPhaseFrequencySeries* pop;
  while( (pop = list) ){
    if( pop->freqseries ){ /* Destroying the CAmpPhaseFrequencySeries data */
      CAmpPhaseFrequencySeries_Cleanup( pop->freqseries );
    }
    /* Notice that the mode indices l and m are not freed, like in SphHarmTimeSeries struct indices l and m */
    list = pop->next;
    free( pop );
  }
}

/***************** Functions for the ListmodesCAmpPhaseSpline structure ****************/
ListmodesCAmpPhaseSpline* ListmodesCAmpPhaseSpline_AddModeNoCopy(
	   ListmodesCAmpPhaseSpline* appended,  /* List structure to prepend to */
	   CAmpPhaseSpline* splines,  /* data to contain */
	   int l, /* major mode number */
	   int m  /* minor mode number */)
{
    ListmodesCAmpPhaseSpline* list;
    /* Check if the node with this mode already exists */
    list = appended;
    while( list ){
      if( l == list->l && m == list->m ){
	break;
      }
      list = list->next;
    }
    if( list ){ /* We don't allow for the case where the mode already exists in the list*/
      printf("Error: Tried to add an already existing mode to a ListmodesCAmpPhaseSpline ");
      return(NULL);
    } else { /* In that case, we do NOT COPY the input interpolated data, which therefore can't be
		used anywhere else; this will be acceptable as these operations will only be done
		when interpolating the initialization data */
      list = malloc( sizeof(ListmodesCAmpPhaseSpline) );
    }
    list->l = l;
    list->m = m;
    if( splines ){
      list->splines = splines;
    } else {
      list->splines = NULL;
    }
    if( appended ){
      list->next = appended;
    } else {
        list->next = NULL;
    }
    return list;
}
/* Get the element of a ListmodesCAmpPhaseSpline with a given index */
ListmodesCAmpPhaseSpline* ListmodesCAmpPhaseSpline_GetMode(
	   ListmodesCAmpPhaseSpline* const list,  /* List structure to get a particular mode from */
	   int l, /*< major mode number */
	   int m  /*< minor mode number */ )
{
    if( !list ) return NULL;

    ListmodesCAmpPhaseSpline *itr = list;
    while( itr->l != l || itr->m != m ){
        itr = itr->next;
        if( !itr ) return NULL;
    }
    return itr; /* The element returned is itself a pointer to a ListmodesCAmpPhaseSpline */
}
void ListmodesCAmpPhaseSpline_Destroy(
	   ListmodesCAmpPhaseSpline* list  /* List structure to destroy; notice that the data is destroyed too */
)
{
  ListmodesCAmpPhaseSpline* pop;
  while( (pop = list) ){
    if( pop->splines ){ /* Destroying the CAmpPhaseSpline data */
      CAmpPhaseSpline_Cleanup( pop->splines );
    }
    /* Notice that the mode indices l and m are not freed, like in SphHarmTimeSeries struct indices l and m */
    list = pop->next;
    free( pop );
  }
}

/***********************************************************************/
/**************** I/O functions for internal structures ****************/

/* Read waveform Real time series */
int Read_RealTimeSeries(RealTimeSeries** timeseries, const char dir[], const char file[], const int nblines, const int binary)
{
  /* Initalize and read input */
  int ret;
  gsl_matrix* inmatrix =  gsl_matrix_alloc(nblines, 2);
  if(!binary) ret = Read_Text_Matrix(dir, file, inmatrix);
  else ret = Read_Matrix(dir, file, inmatrix);

  /* Initialize structures */
  RealTimeSeries_Init(timeseries, nblines);

  /* Set values */
  gsl_vector_view timesview = gsl_matrix_column(inmatrix, 0);
  gsl_vector_view hview = gsl_matrix_column(inmatrix, 1);
  gsl_vector_memcpy((*timeseries)->times, &timesview.vector);
  gsl_vector_memcpy((*timeseries)->h, &hview.vector);

  /* Clean up */
  gsl_matrix_free(inmatrix);

  return ret;
}

/* Read waveform Amp/Phase time series */
int Read_AmpPhaseTimeSeries(AmpPhaseTimeSeries** timeseries, const char dir[], const char file[], const int nblines, const int binary)
{
  /* Initalize and read input */
  int ret;
  gsl_matrix* inmatrix =  gsl_matrix_alloc(nblines, 3);
  if(!binary) ret = Read_Text_Matrix(dir, file, inmatrix);
  else ret = Read_Matrix(dir, file, inmatrix);

  /* Initialize structures */
  AmpPhaseTimeSeries_Init(timeseries, nblines);

  /* Set values */
  gsl_vector_view timesview = gsl_matrix_column(inmatrix, 0);
  gsl_vector_view hampview = gsl_matrix_column(inmatrix, 1);
  gsl_vector_view hphaseview = gsl_matrix_column(inmatrix, 2);
  gsl_vector_memcpy((*timeseries)->times, &timesview.vector);
  gsl_vector_memcpy((*timeseries)->h_amp, &hampview.vector);
  gsl_vector_memcpy((*timeseries)->h_phase, &hphaseview.vector);

  /* Clean up */
  gsl_matrix_free(inmatrix);

  return ret;
}

/* Read waveform Re/Im time series */
int Read_ReImTimeSeries(ReImTimeSeries** timeseries, const char dir[], const char file[], const int nblines, const int binary)
{
  /* Initalize and read input */
  int ret;
  gsl_matrix* inmatrix =  gsl_matrix_alloc(nblines, 3);
  if(!binary) ret = Read_Text_Matrix(dir, file, inmatrix);
  else ret = Read_Matrix(dir, file, inmatrix);

  /* Initialize structures */
  ReImTimeSeries_Init(timeseries, nblines);

  /* Set values */
  gsl_vector_view timesview = gsl_matrix_column(inmatrix, 0);
  gsl_vector_view hrealview = gsl_matrix_column(inmatrix, 1);
  gsl_vector_view himagview = gsl_matrix_column(inmatrix, 2);
  gsl_vector_memcpy((*timeseries)->times, &timesview.vector);
  gsl_vector_memcpy((*timeseries)->h_real, &hrealview.vector);
  gsl_vector_memcpy((*timeseries)->h_imag, &himagview.vector);

  /* Clean up */
  gsl_matrix_free(inmatrix);

  return ret;
}

/* Output Re/Im frequency series */
int Write_ReImFrequencySeries(const char dir[], const char file[], ReImFrequencySeries* freqseries, const int binary)
{
  /* Initialize output */
  /* Note: assumes hplus, hcross have same length as expected */
  int nbfreq = freqseries->freq->size;
  gsl_matrix* outmatrix = gsl_matrix_alloc(nbfreq, 3);

  /* Set output matrix */
  gsl_matrix_set_col(outmatrix, 0, freqseries->freq);
  gsl_matrix_set_col(outmatrix, 1, freqseries->h_real);
  gsl_matrix_set_col(outmatrix, 2, freqseries->h_imag);

  /* Output */
  int ret;
  if (!binary) ret = Write_Text_Matrix(dir, file, outmatrix);
  else ret = Write_Matrix(dir, file, outmatrix);

  return ret;
}

/* Output real time series */
int Write_RealTimeSeries(const char dir[], const char file[], RealTimeSeries* timeseries, int binary)
{
  /* Initialize output */
  int nbtimes = timeseries->times->size;
  gsl_matrix* outmatrix = gsl_matrix_alloc(nbtimes, 2);

  /* Set data */
  gsl_matrix_set_col(outmatrix, 0, timeseries->times);
  gsl_matrix_set_col(outmatrix, 1, timeseries->h);

  /* Output */
  int ret;
  if(!binary) ret = Write_Text_Matrix(dir, file, outmatrix);
  else ret = Write_Matrix(dir, file, outmatrix);

  return ret;
}

/* Output Amp/Phase time series */
int Write_AmpPhaseTimeSeries(const char dir[], const char file[], AmpPhaseTimeSeries* timeseries, int binary)
{
  /* Initialize output */
  int nbtimes = timeseries->times->size;
  gsl_matrix* outmatrix = gsl_matrix_alloc(nbtimes, 3);

  /* Set data */
  gsl_matrix_set_col(outmatrix, 0, timeseries->times);
  gsl_matrix_set_col(outmatrix, 1, timeseries->h_amp);
  gsl_matrix_set_col(outmatrix, 2, timeseries->h_phase);

  /* Output */
  int ret;
  if(!binary) ret = Write_Text_Matrix(dir, file, outmatrix);
  else ret = Write_Matrix(dir, file, outmatrix);

  return ret;
}

/* Output Re/Im time series */
int Write_ReImTimeSeries(const char dir[], const char file[], ReImTimeSeries* timeseries, int binary)
{
  /* Initialize output */
  int nbtimes = timeseries->times->size;
  gsl_matrix* outmatrix = gsl_matrix_alloc(nbtimes, 3);

  /* Set data */
  gsl_matrix_set_col(outmatrix, 0, timeseries->times);
  gsl_matrix_set_col(outmatrix, 1, timeseries->h_real);
  gsl_matrix_set_col(outmatrix, 2, timeseries->h_imag);

  /* Output */
  int ret;
  if(!binary) ret = Write_Text_Matrix(dir, file, outmatrix);
  else ret = Write_Matrix(dir, file, outmatrix);

  return ret;
}

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
/*
 * 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.
 *
 * Written (W) 1999-2009 Soeren Sonnenburg
 * Written (W) 1999-2008 Gunnar Raetsch
 * Written (W) 2006-2007 Mikio L. Braun
 * Written (W) 2008 Jochen Garcke
 * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
 */

#ifndef _LAPACK_H__
#define _LAPACK_H__

#include <shogun/lib/config.h>
#include <shogun/lib/common.h>

#ifdef HAVE_LAPACK

extern "C" {

#ifdef HAVE_MKL
#include <mkl_cblas.h>
#include <mkl_lapack.h>
#elif defined(HAVE_MVEC)
//FIXME: Accelerate framework's vForce.h forward declares
// std::complex<> classes that causes major errors
// in c++11 mode and Eigen3
// this define basically disables the include of vForce.h
#ifdef HAVE_CXX11
#define __VFORCE_H 1
#endif
#include <Accelerate/Accelerate.h>
#else
#include <cblas.h>
#endif

#ifdef HAVE_ACML
#include <acml.h>
#endif

#ifdef HAVE_ATLAS
#include <clapack.h>
#else
// ACML and MKL do not provide clapack_* routines
// double precision
int clapack_dpotrf(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
		const int N, double *A, const int lda);
int clapack_dposv(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
		const int N, const int NRHS, double *A, const int lda,
		double *B, const int ldb);
int clapack_dpotri(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
		const int N, double *A, const int LDA);
int clapack_dpotrs(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
                   const int N, const int NRHS, double *A, const int lda,
                   double *B, const int ldb);
int clapack_dgetrf(const CBLAS_ORDER Order, const int M, const int N,
                   double *A, const int lda, int *ipiv);
int clapack_dgetri(const CBLAS_ORDER Order, const int N, double *A,
                   const int lda, int *ipiv);
int clapack_dgetrs(const CBLAS_ORDER Order, const CBLAS_TRANSPOSE Transpose,
                   const int N, const int NRHS, double *A, const int lda,
                   int *ipiv, double *B, const int ldb);
#endif

namespace shogun
{
#ifndef DOXYGEN_SHOULD_SKIP_THIS
// double precision
void wrap_dsyev(char jobz, char uplo, int n, double *a, int lda,
		double *w, int *info);
void wrap_dgesvd(char jobu, char jobvt, int m, int n, double *a, int lda,
		double *sing, double *u, int ldu, double *vt, int ldvt,
		int *info);
void wrap_dgeqrf(int m, int n, double *a, int lda, double *tau, int *info);
void wrap_dorgqr(int m, int n, int k, double *a, int lda, double *tau, int *info);
void wrap_dsyevr(char jobz, char uplo, int n, double *a, int lda, int il, int iu,
                 double *eigenvalues, double *eigenvectors, int *info);
void wrap_dsygvx(int itype, char jobz, char uplo, int n, double *a, int lda, double *b,
                 int ldb, int il, int iu, double *eigenvalues, double *eigenvectors, int *info);
void wrap_dstemr(char jobz, char range, int n, double* d__, double *e, double vl, double vu,
		int il,	int iu, int* m, double* w, double* z__, int ldz, int nzc, int *isuppz,
		int tryrac, int *info);
#endif
}

// only MKL, ACML and Mac OS vector library provide a header file for the lapack routines
#if !defined(HAVE_ACML) && !defined(HAVE_MKL) && !defined(HAVE_MVEC)
// double precision
int dsyev_(char*, char*, int*, double*, int*, double*, double*, int*, int*);
int dgesvd_(char* jobu, char* jobvt, int* m, int* n, double* a, int* lda,
		double* s, double* u, int* ldu, double* vt, int* ldvt, double* work,
		int* lwork, int* info);
int dposv_(const char *uplo, const int *n, const int *nrhs, double *a, const int *lda, double *b, const int *ldb, int *info);
int dpotrf_(const char *uplo, int *n, double *a, int * lda, int *info);
int dpotri_(const char *uplo, int *n, double *a, int * lda, int *info);
int dgetri_(int *n, double *a, int *lda, int *ipiv, double *work, int *lwork, int *info);
int dgetrf_(int *m, int *n, double *a, int *lda, int *ipiv, int *info);
int dgeqrf_(int*, int*, double*, int*, double*, double*, int*, int*);
int dorgqr_(int*, int*, int*, double*, int*, double*, double*, int*, int*);
int dsyevr_(const char*, const char*, const char*, int*, double*, int*,
            double*, double*, int*, int*, double*, int*, double*, double*,
            int*, int*, double*, int*, int*, int*, int*);
int dgetrs_(const char*, int*, int*, double*, int*, int*, double*, int*, int*);
int dpotrs_(const char*, int*, int*, double*, int*, double*, int*, int*);
int dsygvx_(int*, const char*, const char*, const char*, int*, double*, int*,
            double*, int*, double* , double*, int*, int*, double*,
            int*, double*, double*, int*, double*, int*, int*, int*, int*);
int dstemr_(char*, char*, int*, double*, double*, double*, double*, int*,
	int*, int*, double*, double*, int*, int*, int*, int*, double*,
	int*, int*, int*, int*);
#endif
}

#endif //HAVE_LAPACK
#endif //_LAPACK_H__

max_stars_repo_name: srgnuclear/shogun

max_stars_repo_licenses: Ruby MIT

max_issues_repo_name: waderly/shogun

max_issues_repo_licenses: Ruby MIT

max_forks_repo_name: waderly/shogun

max_forks_repo_licenses: Ruby MIT


text:
/* Copyright (c) 2011-2012, Jérémy Fix. All rights reserved. */

/* Redistribution and use in source and binary forms, with or without */
/* modification, are permitted provided that the following conditions are met: */

/* * Redistributions of source code must retain the above copyright notice, */
/* this list of conditions and the following disclaimer. */
/* * Redistributions in binary form must reproduce the above copyright notice, */
/* this list of conditions and the following disclaimer in the documentation */
/* and/or other materials provided with the distribution. */
/* * None of the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. */

/* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND */
/* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED */
/* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE */
/* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE */
/* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL */
/* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR */
/* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */
/* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, */
/* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE */
/* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */

#ifndef EKF_H
#define EKF_H

#include "ekf_types.h"
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

namespace ekf
{

  void ekf_init(ekf_param &p, ekf_state &s)
  {

    // Matrices for the state
    s.xk = gsl_vector_alloc(p.n);
    gsl_vector_set_zero(s.xk);

    s.xkm = gsl_vector_alloc(p.n);
    gsl_vector_set_zero(s.xkm);

    s.Pxk = gsl_matrix_alloc(p.n, p.n);
    gsl_matrix_set_zero(s.Pxk);

    s.Fxk = gsl_matrix_alloc(p.n, p.n);
    gsl_matrix_set_zero(s.Fxk);

    s.Rv = gsl_matrix_alloc(p.n, p.n);
    gsl_matrix_set_zero(s.Rv);

    // Matrices for the observations
    s.yk = gsl_vector_alloc(p.no);
    gsl_vector_set_zero(s.yk);

    s.ino_yk = gsl_vector_alloc(p.no);
    gsl_vector_set_zero(s.ino_yk);

    s.Hyk = gsl_matrix_alloc(p.no, p.n);
    gsl_matrix_set_zero(s.Hyk);

    s.Rn = gsl_matrix_alloc(p.no, p.no);
    gsl_matrix_set_zero(s.Rn);

    // Matrices for the kalman gain and the updates
    s.Kk = gsl_matrix_alloc(p.n, p.no);
    gsl_matrix_set_zero(s.Kk);

    // Temporary matrices
    s.temp_n_n = gsl_matrix_alloc(p.n, p.n);
    gsl_matrix_set_zero(s.temp_n_n);

    s.temp_n_1 = gsl_matrix_alloc(p.n, 1);
    gsl_matrix_set_zero(s.temp_n_1);

    s.temp_no_no = gsl_matrix_alloc(p.no, p.no);
    gsl_matrix_set_zero(s.temp_no_no);

    s.temp_n_no = gsl_matrix_alloc(p.n, p.no);
    gsl_matrix_set_zero(s.temp_n_no);

    s.temp_2_n_n = gsl_matrix_alloc(p.n, p.n);
    gsl_matrix_set_zero(s.temp_2_n_n);

    s.temp_no = gsl_vector_alloc(p.no);
    gsl_vector_set_zero(s.temp_no);

    // Initialize the noises
    //gsl_matrix_set_identity(s.Rv);
    //gsl_matrix_scale(s.Rv, p.evolution_noise);
    p.evolution_noise->init(p, s);
    //printf("Evolution noise : max = %e , min = %e \n", gsl_matrix_max(s.Rv), gsl_matrix_min(s.Rv));

    gsl_matrix_set_identity(s.Rn);
    gsl_matrix_scale(s.Rn, p.observation_noise);

    // Initialize the covariance of the parameters
    gsl_matrix_set_identity(s.Pxk);
    gsl_matrix_scale(s.Pxk, p.prior_pk);

  }

  void ekf_free(ekf_param &p, ekf_state &s)
  {
    gsl_vector_free(s.xk);
    gsl_vector_free(s.xkm);
    gsl_matrix_free(s.Pxk);
    gsl_matrix_free(s.Fxk);
    gsl_matrix_free(s.Rv);

    gsl_vector_free(s.yk);
    gsl_vector_free(s.ino_yk);
    gsl_matrix_free(s.Hyk);
    gsl_matrix_free(s.Rn);

    gsl_matrix_free(s.Kk);

    gsl_matrix_free(s.temp_n_n);
    gsl_matrix_free(s.temp_n_1);
    gsl_matrix_free(s.temp_no_no);
    gsl_matrix_free(s.temp_n_no);
    gsl_matrix_free(s.temp_2_n_n);
    gsl_vector_free(s.temp_no);
  }
  template<typename FFunc,
    typename DFFunc,
    typename HFunc,
    typename DHFunc>
    void ekf_iterate(ekf_param &p, ekf_state &s,
		     FFunc f, DFFunc df,
		     HFunc h, DHFunc dh,
		     gsl_vector* yk)
  {
    /****************************/
    /***** Prediction step  *****/
    /****************************/

    // Compute the Jacobian of the evolution
    // Eq. 2.34
    df(s.params, s.xk, s.Fxk);

    // Compute the predicted state mean and covariance
    // Eq. 2.36
    // s.xk will now hold the predicted state !
    f(s.params, s.xk, s.xkm);

    // Eq. 2.37
    // s.Pxk will now hold Pxk^-
    gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, s.Pxk, s.Fxk, 0.0, s.temp_n_n);
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, s.Fxk, s.temp_n_n, 0.0, s.Pxk);
    gsl_matrix_add(s.Pxk, s.Rv);

    /****************************/
    /***** Correction step  *****/
    /****************************/

    // Compute the Jacobian of the observation model
    // Eq. 2.38
    dh(s.params, s.xkm, s.Hyk);

    if(!p.observation_gradient_is_diagonal)
      {
	// Update the estimates
	// Eq 2.40
	// 1 - Compute H.P^-.H^T + Rn
	gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, s.Pxk, s.Hyk, 0.0, s.temp_n_no);
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, s.Hyk, s.temp_n_no, 0.0, s.temp_no_no);
	gsl_matrix_add(s.temp_no_no, s.Rn);
	// 2 - Compute its inverse
	gsl_linalg_cholesky_decomp(s.temp_no_no);
	gsl_linalg_cholesky_invert(s.temp_no_no);

	// 3 - Compute P^-.H ^T.( H P H^T + R)^-1
	gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, s.Hyk, s.temp_no_no, 0.0, s.temp_n_no);
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, s.Pxk, s.temp_n_no, 0.0, s.Kk);
      }
    else
      {
	// We can make some simplifications when computing H P H^T = H P H
	// and also when computing P . H . (H P H^T + R) ^-1

	// Update the estimates
	// Eq 2.40
	// 1 - Compute H.P^-.H^T + Rn
	for(int i = 0 ; i < p.no ; ++i)
	  {
	    for(int j = 0 ; j < p.no ; ++j)
	      {
		gsl_matrix_set(s.temp_no_no, i, j, gsl_matrix_get(s.Hyk, i, i) * gsl_matrix_get(s.Hyk, j, j) * gsl_matrix_get(s.Pxk, i, j) + gsl_matrix_get(s.Rn, i, j));
	      }
	  }
	// Compute its inverse : U = (H P H^T + R) ^-1
	gsl_linalg_cholesky_decomp(s.temp_no_no);
	gsl_linalg_cholesky_invert(s.temp_no_no);

	// 3 - Compute P^- H^T . U

	// Compute H^T . U
	for(int i = 0 ; i < p.no; ++i)
	  for(int j = 0 ; j < p.no ; ++j)
	    gsl_matrix_set(s.temp_n_no, i, j, gsl_matrix_get(s.Hyk, i,i) * gsl_matrix_get(s.temp_no_no, i, j));
	for(int i = p.no ; i < p.n; ++i)
	  for(int j = 0 ; j < p.no ; ++j)
	    gsl_matrix_set(s.temp_n_no, i, j, 0.0);
	// Compute P^- . (H^T . U)
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, s.Pxk, s.temp_n_no, 0.0, s.Kk);

      }
    // Update the current estimate
    // Eq 2.41
    // 1 - We need the observations
    h(s.params, s.xkm, s.yk);
    gsl_vector_memcpy(s.ino_yk, yk);
    gsl_vector_sub(s.ino_yk, s.yk);
    gsl_vector_memcpy(s.xk, s.xkm);
    gsl_blas_dgemv(CblasNoTrans, 1.0, s.Kk, s.ino_yk, 1.0,s.xk);

    // Update the variance/covariance matrix
    // Compute  -Kk * Hk
    if(!p.observation_gradient_is_diagonal)
      {
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1.0, s.Kk, s.Hyk, 0.0, s.temp_n_n);

	// Add identity : I - Kk * Hk
	for(int i = 0 ; i < p.n ; ++i)
	  gsl_matrix_set(s.temp_n_n, i, i,gsl_matrix_get(s.temp_n_n,i,i) + 1.0);
      }
    else
      {
	// We can make some simplifications when H is diagonal
	for(int i = 0 ; i < p.n ; ++i)
	  {
	    for(int j = 0 ; j < p.no ; ++j)
	      {
		gsl_matrix_set(s.temp_n_n, i, j, (i==j?1.0:0.0) - gsl_matrix_get(s.Kk, i,j) * gsl_matrix_get(s.Hyk, j,j) );
	      }
	    for(int j = p.no ; j < p.n ; ++j)
	      {
		gsl_matrix_set(s.temp_n_n, i, j, (i == j ? 1.0 : 0.0));
	      }
	  }
      }

    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, s.temp_n_n, s.Pxk, 0.0, s.temp_2_n_n);
    gsl_matrix_memcpy(s.Pxk, s.temp_2_n_n);

    /***********************************/
    /***** Evolution noise update  *****/
    /***********************************/

    p.evolution_noise->updateEvolutionNoise(p, s);
  }
}

#endif // EKF_H


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2017 Couchbase, Inc.
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 */

#pragma once

#include "config.h"

#include <cJSON.h>
#include <platform/socket.h>
#include <gsl/gsl>

class NetworkInterface {
public:
    /// Should a protocol be enabled, and if so how?
    enum class Protocol {
        Off, /// Do not enable protocol.
        Optional, /// Protocol should be enabled, failure is non-fatal.
        Required, /// Protocol must be enabled; failure is fatal.
    };

    NetworkInterface() = default;
    explicit NetworkInterface(gsl::not_null<const cJSON*> json);

    std::string host;
    struct {
        std::string key;
        std::string cert;
    } ssl;
    int maxconn = 1000;
    int backlog = 1024;
    in_port_t port = 11211;
    Protocol ipv6 = Protocol::Optional;
    Protocol ipv4 = Protocol::Optional;
    bool tcp_nodelay = true;
    bool management = false;
};

std::string to_string(const NetworkInterface::Protocol& proto);