text:
#ifndef _PS_
#define _PS_

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <ctype.h>
#include <math.h>
#include <unistd.h>
#include "../Parameter_files/COSMOLOGY.H"
#include "../Parameter_files/INIT_PARAMS.H"
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include "cosmo_progs.c"
#include "misc.c"

/* New in v1.1 */
#define ERFC_NPTS (int) 75
#define ERFC_PARAM_DELTA (float) 0.1
static double log_erfc_table[ERFC_NPTS], erfc_params[ERFC_NPTS];
static gsl_interp_accel *erfc_acc;
static gsl_spline *erfc_spline;

double sigma_norm, R, theta_cmb, omhh, z_equality, y_d, sound_horizon, alpha_nu, f_nu, f_baryon, beta_c, d2fact, R_CUTOFF, DEL_CURR, SIG_CURR;

/*****     FUNCTION PROTOTYPES     *****/
double init_ps(); /* initialize global variables, MUST CALL THIS FIRST!!! returns R_CUTOFF */
void free_ps(); /* deallocates the gsl structures from init_ps */
double splined_erfc(double); /* returns erfc for x>=0, using cubic spline in logy-x space */
double deltolindel(float del, float z); /* converts a non-linear overdensity, del, at z to a linear overdensity at z=0 */
double lindeltodel(float lindel, float z); /* converts a linear overdensity, del, at z=0 to a non-linear overdensity at redshift z */
double power_in_k(double k); /* Returns the value of the linear power spectrum density (i.e. <|delta_k|^2>/V) at a given k mode at z=0 */
double RtoM(double); /* R in Mpc, M in Msun */
double MtoR(double); /* R in Mpc, M in Msun */
double M_J_WDM(); /* returns the "effective Jeans mass" corresponding to the gas analog of WDM ; eq. 10 in BHO 2001 */
double sheth_delc(double del, double sig);
double dNdM_st(double z, double M);
double dNdM(double z, double M);
double dnbiasdM(double M, float z, double M_o, float del_o); /* dnbiasdM */
double FgtrM(double z, double M);  //calculates the fraction of mass contained in haloes with mass > M at redshift z
double FgtrM_st(double z, double M);  //calculates the fraction of mass contained in haloes with mass > M at redshift z, with Sheth-Tormen correction
double FgtrM_bias(double z, double M, double del_bias, double sig_bias);  //calculates the fraction of mass contained in haloes with mass > M at redshift z, in regions with a linear overdensity of del_bias, and standard deviation sig_bias
double sigmaparam_FgtrM_bias(float z, float sigsmallR, float del_bias, float sig_bias);/* Uses sigma parameters instead of Mass for scale */

double FgtrM_bias_BL08(double z, double M, double del_bias, double sig_bias); // as above, but this version uses the hybrid perscription of Barkana & Loeb 2004 (specifically the separate integral version of eq. 2 in Barkana & Loeb 2008)
double dicke(double z); //calculates the dicke growth function at redshift z
double ddickedz(double z); /* Redshift derivative of the growth function at z */
double ddickedt(double z); /* Time derivative of the growth function at z */
double sigma_z0(double M); //calculates sigma at z=0 (no dicke)
double dsigmasqdm_z0(double M); //calculates d(sigma^2)/dm at z=0 (i.e. does not include dicke growth)
double TFmdm(double k); //Eisenstien & Hu power spectrum transfer function
void TFset_parameters();
float get_R_c();  // returns R_CUTOFF
double get_M_min_ion(float z);
/***************************************/

/* Returns the minimum source mass for ionizing sources, according to user specifications */
double get_M_min_ion(float z){
  double MMIN;
  if (ION_M_MIN < 0){ // use the virial temperature for Mmin
    if (ION_Tvir_MIN < 9.99999e3) // neutral IGM
      MMIN = TtoM(z, ION_Tvir_MIN, 1.22);
    else // ionized IGM
      MMIN = TtoM(z, ION_Tvir_MIN, 0.6);
  }
  else if (ION_Tvir_MIN < 0){ // use the mass
    MMIN = ION_M_MIN;
  }
  else{
    fprintf(stderr, "You have to \"turn-off\" either the ION_M_MIN or \
                     the ION_Tvir_MIN option in ANAL_PARAMS.H\nAborting...\n");
    return -1;
  }
  // check for WDM
  if (P_CUTOFF && ( MMIN < M_J_WDM()))
    MMIN = M_J_WDM();
  //  printf("Mmin is %e\n", MMIN);
  return MMIN;
}

/* Returns the minimum source mass for x-ray sources, according to user specifications */
double get_M_min_xray(float z){
  double MMIN;
  if (X_RAY_Tvir_MIN < 9.99999e3) //neutral IGM
    MMIN = TtoM(z, X_RAY_Tvir_MIN, 1.22);
  else // ionized IGM
    MMIN = TtoM(z, X_RAY_Tvir_MIN, 0.6);

  // check for WDM
  if (P_CUTOFF && ( MMIN < M_J_WDM()))
    MMIN = M_J_WDM();
  //  printf("Mmin is %e\n", MMIN);
  return MMIN;
}

/* returns the "effective Jeans mass" in Msun
   corresponding to the gas analog of WDM ; eq. 10 in Barkana+ 2001 */
double M_J_WDM(){
  double z_eq, fudge=60;
  if (!P_CUTOFF)
    return 0;
  z_eq = 3600*(OMm-OMb)*hlittle*hlittle/0.15;
  return fudge*3.06e8 * (1.5/g_x) * sqrt((OMm-OMb)*hlittle*hlittle/0.15) * pow(M_WDM, -4) * pow(z_eq/3000.0, 1.5);
}

/* converts a non-linear overdensity, del, at z to a linear overdensity at z=0 */
double deltolindel(float del, float z){
  float onepdel = 1.0+del;
  return ( 1.68647 - 1.35*pow(onepdel,-2/3.0) + 0.78785*pow(onepdel,-0.58661) - 1.12431*pow(onepdel,-0.5) )/dicke(z);
}

/* converts a linear overdensity, del, at z=0 to a non-linear overdensity at redshift z */
double lindeltodel(float lindel, float z){
  float prev_lindelguess, delcrit, delguess;
  float lindelguess, delmin, delmax, epsilon = 1.0e-7;

  // set the critical density corresponding to virialization
  // this will be maximum allowed del
  delcrit = Deltac_nonlinear(z)*rho_critz(z)/(OMm*RHOcrit*pow(1+z, 3)) - 1;

  delmin = -1;
  delmax = 500;
  prev_lindelguess = -1e10;
  while (1){
    delguess = 0.5*(delmax+delmin);
    lindelguess = deltolindel(delguess, z);
    //fprintf(stderr, "%e\t%e\n", delmin, delmax);
	      //    fprintf(stderr, "%e\t%e\t%e\n\n", delguess, lindelguess, lindel);
    if ((fabs((lindelguess-lindel)/lindel) < epsilon ) ||
	(fabs(lindelguess-lindel) < epsilon ) ||
	(fabs(prev_lindelguess - lindelguess) < TINY ))// close enough, or resolution loop
      return delguess;

    if (lindelguess > lindel)
      delmax = delguess;
    else
      delmin = delguess;

    // check if we are above delcrit (see above)
    if (delmin > delcrit){
      //      printf("exced max at lindel=%e\n", lindel);
      return delcrit;
    }

    prev_lindelguess = lindelguess;
  }
}

/* R in Mpc, M in Msun */
double RtoM(double R){
  // set M according to M<->R conversion defined by the filter type in ../Parameter_files/COSMOLOGY.H
  if (FILTER == 0) //top hat M = (4/3) PI <rho> R^3
    return (4.0/3.0)*PI*pow(R,3)*(OMm*RHOcrit);
  else if (FILTER == 1) //gaussian: M = (2PI)^1.5 <rho> R^3
    return pow(2*PI, 1.5) * OMm*RHOcrit * pow(R, 3);
  else // filter not defined
    fprintf(stderr, "No such filter = %i.\nResults are bogus.\n", FILTER);
  return -1;
}

/* R in Mpc, M in Msun */
double MtoR(double M){
  // set R according to M<->R conversion defined by the filter type in ../Parameter_files/COSMOLOGY.H
  if (FILTER == 0) //top hat M = (4/3) PI <rho> R^3
    return pow(3*M/(4*PI*OMm*RHOcrit), 1.0/3.0);
  else if (FILTER == 1) //gaussian: M = (2PI)^1.5 <rho> R^3
    return pow( M/(pow(2*PI, 1.5) * OMm * RHOcrit), 1.0/3.0 );
  else // filter not defined
    fprintf(stderr, "No such filter = %i.\nResults are bogus.\n", FILTER);
  return -1;
}

/* equation (5) from jenkis et al. (2001) */
double f_jenkins(float del, double sigsq){
  if (del < 0){  fprintf(stderr, "ERROR:  In function f_jenkins del_o must be less than del_1 = del_crit/dicke(z)!\nAborting...\n"); return 0; }

  //  fprintf(stderr, "%f\t%f\n", del, sqrt(sigsq));
  return sqrt(2/PI) * del/sqrt(sigsq) * pow(E, -0.5*del*del/sigsq);
}

float get_R_c(){
  return R_CUTOFF;
}

/* sheth correction to delta crit */
double sheth_delc(double del, double sig){
  return sqrt(SHETH_a)*del*(1 + SHETH_b*pow(sig*sig/(SHETH_a*del*del), SHETH_c));
}

/* dnbiasdM */
double dnbiasdM(double M, float z, double M_o, float del_o){
  double sigsq, del, sig_one, sig_o;

  if ((M_o-M) < TINY){
    fprintf(stderr, "WARNING:  In function dnbiasdM: M must be less than M_o!\nAborting...\n");
    return -1;
  }
  del = Deltac/dicke(z) - del_o;
  if (del < 0){  fprintf(stderr, "ERROR:  In function dnbiasdM: del_o must be less than del_1 = del_crit/dicke(z)!\nAborting...\n"); return 0; }
  sig_o = sigma_z0(M_o);
  sig_one = sigma_z0(M);
  sigsq = sig_one*sig_one - sig_o*sig_o;
  return -(RHOcrit*OMm)/M /sqrt(2*PI) *del*pow(sigsq,-1.5)*pow(E, -0.5*del*del/sigsq)*dsigmasqdm_z0(M);
}

/*
  FUNCTION dNdM(z, M)
  Computes the Press_schechter mass function with Sheth-Torman correction for ellipsoidal collapse at
  redshift z, and dark matter halo mass M (in solar masses).

  The return value is the number density per unit mass of halos in the mass range M to M+dM in units of:
  comoving Mpc^-3 Msun^-1

  Reference: Sheth, Mo, Torman 2001
*/
double dNdM_st(double z, double M){
  double sigma, dsigmadm, nuhat, dicke_growth;

  dicke_growth = dicke(z);
  sigma = sigma_z0(M) * dicke_growth;
  dsigmadm = dsigmasqdm_z0(M) * dicke_growth*dicke_growth/(2.0*sigma);
  nuhat = sqrt(SHETH_a) * Deltac / sigma;

  return (-OMm*RHOcrit/M) * (dsigmadm/sigma) * sqrt(2/PI)*SHETH_A * (1+ pow(nuhat, -2*SHETH_p)) * nuhat * pow(E, -nuhat*nuhat/2.0);
}

/*
  FUNCTION dNdM(z, M)
  Computes the Press_schechter mass function at
  redshift z, and dark matter halo mass M (in solar masses).

  The return value is the number density per unit mass of halos in the mass range M to M+dM in units of:
  comoving Mpc^-3 Msun^-1

  Reference: Padmanabhan, pg. 214
*/
double dNdM(double z, double M){
  double sigma, dsigmadm, dicke_growth;

  dicke_growth = dicke(z);
  sigma = sigma_z0(M) * dicke_growth;
  dsigmadm = dsigmasqdm_z0(M) * (dicke_growth*dicke_growth/(2*sigma));

  return (-OMm*RHOcrit/M) * sqrt(2/PI) * (Deltac/(sigma*sigma)) * dsigmadm * pow(E, -(Deltac*Deltac)/(2*sigma*sigma));
}

/*
  FUNCTION FgtrM_st(z, M)
  Computes the fraction of mass contained in haloes with mass > M at redshift z
  Uses Sheth-Torman correction
*/
double dFdlnM_st (double lnM, void *params){
  double z = *(double *)params;
  double M = exp(lnM);
  return dNdM_st(z, M) * M * M;
}
double FgtrM_st(double z, double M){
  double result, error, lower_limit, upper_limit;
  gsl_function F;
  double rel_tol  = 0.001; //<- relative tolerance
  gsl_integration_workspace * w
    = gsl_integration_workspace_alloc (1000);

  F.function = &dFdlnM_st;
  F.params = &z;
  lower_limit = log(M);
  upper_limit = log(FMAX(1e16, M*100));

 // printf("%f %f\n",lower_limit, upper_limit);
  gsl_integration_qag (&F, lower_limit, upper_limit, 0, rel_tol,
		       1000, GSL_INTEG_GAUSS61, w, &result, &error);
  gsl_integration_workspace_free (w);

  return result / (OMm*RHOcrit);
}

/*
  FUNCTION FgtrM(z, M)
  Computes the fraction of mass contained in haloes with mass > M at redshift z
*/
double FgtrM(double z, double M){
  double del, sig;

  del = Deltac/dicke(z); //regular spherical collapse delta
  sig = sigma_z0(M);

  return splined_erfc(del / (sqrt(2)*sig));
}

/*
  calculates the fraction of mass contained in haloes with mass > M at redshift z, in regions with a linear overdensity of del_bias, and standard deviation sig_bias
 */
double FgtrM_bias(double z, double M, double del_bias, double sig_bias){
  double del, sig, sigsmallR;

  sigsmallR = sigma_z0(M);

  if (!(sig_bias < sigsmallR)){ // biased region is smaller that halo!
    fprintf(stderr, "FgtrM_bias: Biased region is smaller than halo!\nResult is bogus.\n");
    return 0;
  }

  del = Deltac/dicke(z) - del_bias;
  sig = sqrt(sigsmallR*sigsmallR - sig_bias*sig_bias);

   return splined_erfc(del / (sqrt(2)*sig));
}

/* Uses sigma parameters instead of Mass for scale */
double sigmaparam_FgtrM_bias(float z, float sigsmallR, float del_bias, float sig_bias){
  double del, sig;

  if (!(sig_bias < sigsmallR)){ // biased region is smaller that halo!
    fprintf(stderr, "local_FgtrM_bias: Biased region is smaller than halo!\nResult is bogus.\n");
    return 0;
  }

  del = Deltac/dicke(z) - del_bias;
  sig = sqrt(sigsmallR*sigsmallR - sig_bias*sig_bias);

  return splined_erfc(del / (sqrt(2)*sig));
}

/*
  Calculates the fraction of mass contained in haloes with mass > M at redshift z, in regions with a linear overdensity of del_bias, and standard deviation sig_bias.
  This version uses the hybrid perscription of Barkana & Loeb 2004 (specifically the separate
  integral version of eq. 2 in Barkana & Loeb 2008)
 */
double FgtrM_bias_BL08(double z, double M, double del_bias, double sig_bias){
  return FgtrM_st(z, M) / FgtrM(z, M) * FgtrM_bias(z, M, del_bias, sig_bias);
}

/*
  FUNCTION dicke(z)
  Computes the dicke growth function at redshift z, i.e. the z dependance part of sigma

  References: Peebles, "Large-Scale...", pg.53 (eq. 11.16). Includes omega<=1
  Nonzero Lambda case from Liddle et al, astro-ph/9512102, eqs. 6-8.
  and quintessence case from Wang et al, astro-ph/9804015

  Normalized to dicke(z=0)=1
*/
double dicke(double z){
  double omegaM_z, dick_z, dick_0, x, x_0;
  double tiny = 1e-4;

  if (fabs(OMm-1.0) < tiny){ //OMm = 1 (Einstein de-Sitter)
    return 1.0/(1.0+z);
  }
  else if ( (OMl > (-tiny)) && (fabs(OMl+OMm+OMr-1.0) < 0.01) && (fabs(wl+1.0) < tiny) ){
    //this is a flat, cosmological CONSTANT universe, with only lambda, matter and radiation
    //it is taken from liddle et al.
    omegaM_z = OMm*pow(1+z,3) / ( OMl + OMm*pow(1+z,3) + OMr*pow(1+z,4) );
    dick_z = 2.5*omegaM_z / ( 1.0/70.0 + omegaM_z*(209-omegaM_z)/140.0 + pow(omegaM_z, 4.0/7.0) );
    dick_0 = 2.5*OMm / ( 1.0/70.0 + OMm*(209-OMm)/140.0 + pow(OMm, 4.0/7.0) );
    return dick_z / (dick_0 * (1.0+z));
  }
  else if ( (OMtot < (1+tiny)) && (fabs(OMl) < tiny) ){ //open, zero lambda case (peebles, pg. 53)
    x_0 = 1.0/(OMm+0.0) - 1.0;
    dick_0 = 1 + 3.0/x_0 + 3*log(sqrt(1+x_0)-sqrt(x_0))*sqrt(1+x_0)/pow(x_0,1.5);
    x = fabs(1.0/(OMm+0.0) - 1.0) / (1+z);
    dick_z = 1 + 3.0/x + 3*log(sqrt(1+x)-sqrt(x))*sqrt(1+x)/pow(x,1.5);
    return dick_z/dick_0;
  }
  else if ( (OMl > (-tiny)) && (fabs(OMtot-1.0) < tiny) && (fabs(wl+1) > tiny) ){
    fprintf(stderr, "IN WANG\n");
    return -1;
  }

  fprintf(stderr, "No growth function!!! Output will be fucked up.");
  return -1;
}

/* redshift derivative of the growth function at z */
double ddicke_dz(double z){
  float dz = 1e-10;
  double omegaM_z, ddickdz, dick_0, x, x_0, domegaMdz;

   return (dicke(z+dz)-dicke(z))/dz;
}

/* Time derivative of the growth function at z */
double ddickedt(double z){
  float dz = 1e-10;
  double omegaM_z, ddickdz, dick_0, x, x_0, domegaMdz;
  double tiny = 1e-4;

   return (dicke(z+dz)-dicke(z))/dz/dtdz(z); // lazy non-analytic form getting

  if (fabs(OMm-1.0) < tiny){ //OMm = 1 (Einstein de-Sitter)
    return -pow(1+z,-2)/dtdz(z);
  }
  else if ( (OMl > (-tiny)) && (fabs(OMl+OMm+OMr-1.0) < 0.01) && (fabs(wl+1.0) < tiny) ){
    //this is a flat, cosmological CONSTANT universe, with only lambda, matter and radiation
    //it is taken from liddle et al.
    omegaM_z = OMm*pow(1+z,3) / ( OMl + OMm*pow(1+z,3) + OMr*pow(1+z,4) );
    domegaMdz = omegaM_z*3/(1+z) - OMm*pow(1+z,3)*pow(OMl + OMm*pow(1+z,3) + OMr*pow(1+z,4), -2) * (3*OMm*(1+z)*(1+z) + 4*OMr*pow(1+z,3));
    dick_0 = OMm / ( 1.0/70.0 + OMm*(209-OMm)/140.0 + pow(OMm, 4.0/7.0) );

    ddickdz = (domegaMdz/(1+z)) * (1.0/70.0*pow(omegaM_z,-2) + 1.0/140.0 + 3.0/7.0*pow(omegaM_z, -10.0/3.0)) * pow(1.0/70.0/omegaM_z + (209.0-omegaM_z)/140.0 + pow(omegaM_z, -3.0/7.0) , -2);
    ddickdz -= pow(1+z,-2)/(1.0/70.0/omegaM_z + (209.0-omegaM_z)/140.0 + pow(omegaM_z, -3.0/7.0));

    return ddickdz / dick_0 / dtdz(z);
  }

  fprintf(stderr, "No growth function!!! Output will be fucked up.");
  return -1;
}

/*
  FUNCTION sigma_z0(M)
  Returns the standard deviation of the normalized, density excess (delta(x)) field,
  smoothed on the comoving scale of M (see filter definitions for M<->R conversion).
  The sigma is evaluated at z=0, with the time evolution contained in the dicke(z) factor,
  i.e. sigma(M,z) = sigma_z0(m) * dicke(z)

  normalized so that sigma_z0(M->8/h Mpc) = SIGMA8 in ../Parameter_files/COSMOLOGY.H

  NOTE: volume is normalized to = 1, so this is equvalent to the mass standard deviation

  M is in solar masses

  References: Padmanabhan, pg. 210, eq. 5.107
*/
double dsigma_dk(double k, void *params){
  double p, w, T, gamma, q, aa, bb, cc, kR;

  // get the power spectrum.. choice of 5:
  if (POWER_SPECTRUM == 0){ // Eisenstein & Hu
    T = TFmdm(k);
    // check if we should cuttoff power spectrum according to Bode et al. 2000 transfer function
    if (P_CUTOFF) T *= pow(1 + pow(BODE_e*k*R_CUTOFF, 2*BODE_v), -BODE_n/BODE_v);
    p = pow(k, POWER_INDEX) * T * T;
  }
  else if (POWER_SPECTRUM == 1){ // BBKS
    gamma = OMm * hlittle * pow(E, -OMb - OMb/OMm);
    q = k / (hlittle*gamma);
    T = (log(1.0+2.34*q)/(2.34*q)) *
      pow( 1.0+3.89*q + pow(16.1*q, 2) + pow( 5.46*q, 3) + pow(6.71*q, 4), -0.25);
     p = pow(k, POWER_INDEX) * T * T;
  }
  else if (POWER_SPECTRUM == 2){ // Efstathiou,G., Bond,J.R., and White,S.D.M., MNRAS,258,1P (1992)
    gamma = 0.25;
    aa = 6.4/(hlittle*gamma);
    bb = 3.0/(hlittle*gamma);
    cc = 1.7/(hlittle*gamma);
    p = pow(k, POWER_INDEX) / pow( 1+pow( aa*k + pow(bb*k, 1.5) + pow(cc*k,2), 1.13), 2.0/1.13 );
  }
  else if (POWER_SPECTRUM == 3){ // Peebles, pg. 626
    gamma = OMm * hlittle * pow(E, -OMb - OMb/OMm);
    aa = 8.0 / (hlittle*gamma);
    bb = 4.7 / pow(hlittle*gamma, 2);
    p = pow(k, POWER_INDEX) / pow(1 + aa*k + bb*k*k, 2);
  }
  else if (POWER_SPECTRUM == 4){ // White, SDM and Frenk, CS, 1991, 379, 52
    gamma = OMm * hlittle * pow(E, -OMb - OMb/OMm);
    aa = 1.7/(hlittle*gamma);
    bb = 9.0/pow(hlittle*gamma, 1.5);
    cc = 1.0/pow(hlittle*gamma, 2);
    p = pow(k, POWER_INDEX) * 19400.0 / pow(1 + aa*k + bb*pow(k, 1.5) + cc*k*k, 2);
  }
  else{
    fprintf(stderr, "No such power spectrum defined: %i\nOutput is bogus.\n", POWER_SPECTRUM);
    p = 0;
  }

  // now get the value of the window function
  // NOTE: only use top hat for SIGMA8 normalization
  kR = k*R;
  if ( (FILTER == 0) || (sigma_norm < 0) ){ // top hat
    if ( (kR) < 1.0e-4 ){ w = 1.0;} // w converges to 1 as (kR) -> 0
    else { w = 3.0 * (sin(kR)/pow(kR, 3) - cos(kR)/pow(kR, 2));}
  }
  else if (FILTER == 1){ // gaussian of width 1/R
    w = pow(E, -kR*kR/2.0);
  }
  else {
    fprintf(stderr, "No such filter: %i\nOutput is bogus.\n", FILTER);
    w=0;
  }

  return k*k*p*w*w;
}
double sigma_z0(double M){
  double result, error, lower_limit, upper_limit;
  gsl_function F;
  double rel_tol  = FRACT_FLOAT_ERR*10; //<- relative tolerance
  gsl_integration_workspace * w
    = gsl_integration_workspace_alloc (1000);
  double kstart, kend;

  R = MtoR(M);

  // now lets do the integral for sigma and scale it with sigma_norm
  kstart = 1.0e-99/R;
  kend = 350.0/R;
  lower_limit = kstart;//log(kstart);
  upper_limit = kend;//log(kend);

  F.function = &dsigma_dk;
  gsl_integration_qag (&F, lower_limit, upper_limit, 0, rel_tol,
		       1000, GSL_INTEG_GAUSS61, w, &result, &error);
  gsl_integration_workspace_free (w);

  return sigma_norm * sqrt(result);
}

/*
  Returns the value of the linear power spectrum DENSITY (i.e. <|delta_k|^2>/V)
  at a given k mode linearly extrapolated to z=0
*/
double power_in_k(double k){
  double p, T, gamma, q, aa, bb, cc;

  // get the power spectrum.. choice of 5:
  if (POWER_SPECTRUM == 0){ // Eisenstein & Hu
    T = TFmdm(k);
    // check if we should cuttoff power spectrum according to Bode et al. 2000 transfer function
    if (P_CUTOFF) T *= pow(1 + pow(BODE_e*k*R_CUTOFF, 2*BODE_v), -BODE_n/BODE_v);
    p = pow(k, POWER_INDEX) * T * T;
    //p = pow(k, POWER_INDEX - 0.05*log(k/0.05)) * T * T; //running, alpha=0.05
  }
  else if (POWER_SPECTRUM == 1){ // BBKS
    gamma = OMm * hlittle * pow(E, -OMb - OMb/OMm);
    q = k / (hlittle*gamma);
    T = (log(1.0+2.34*q)/(2.34*q)) *
      pow( 1.0+3.89*q + pow(16.1*q, 2) + pow( 5.46*q, 3) + pow(6.71*q, 4), -0.25);
     p = pow(k, POWER_INDEX) * T * T;
  }
  else if (POWER_SPECTRUM == 2){ // Efstathiou,G., Bond,J.R., and White,S.D.M., MNRAS,258,1P (1992)
    gamma = 0.25;
    aa = 6.4/(hlittle*gamma);
    bb = 3.0/(hlittle*gamma);
    cc = 1.7/(hlittle*gamma);
    p = pow(k, POWER_INDEX) / pow( 1+pow( aa*k + pow(bb*k, 1.5) + pow(cc*k,2), 1.13), 2.0/1.13 );
  }
  else if (POWER_SPECTRUM == 3){ // Peebles, pg. 626
    gamma = OMm * hlittle * pow(E, -OMb - OMb/OMm);
    aa = 8.0 / (hlittle*gamma);
    bb = 4.7 / pow(hlittle*gamma, 2);
    p = pow(k, POWER_INDEX) / pow(1 + aa*k + bb*k*k, 2);
  }
  else if (POWER_SPECTRUM == 4){ // White, SDM and Frenk, CS, 1991, 379, 52
    gamma = OMm * hlittle * pow(E, -OMb - OMb/OMm);
    aa = 1.7/(hlittle*gamma);
    bb = 9.0/pow(hlittle*gamma, 1.5);
    cc = 1.0/pow(hlittle*gamma, 2);
    p = pow(k, POWER_INDEX) * 19400.0 / pow(1 + aa*k + bb*pow(k, 1.5) + cc*k*k, 2);
  }
  else{
    fprintf(stderr, "No such power spectrum defined: %i\nOutput is bogus.\n", POWER_SPECTRUM);
    p = 0;
  }

  return p*TWOPI*PI*sigma_norm*sigma_norm;
}

/*
  FUNCTION dsigmasqdm_z0(M)
  returns  d/dm (sigma^2) (see function sigma), in units of Msun^-1
*/
double dsigmasq_dm(double k, void *params){
  double p, w, T, gamma, q, aa, bb, cc, dwdr, drdm, kR;

  // get the power spectrum.. choice of 5:
  if (POWER_SPECTRUM == 0){ // Eisenstein & Hu ApJ, 1999, 511, 5
    T = TFmdm(k);
    // check if we should cuttoff power spectrum according to Bode et al. 2000 transfer function
    if (P_CUTOFF) T *= pow(1 + pow(BODE_e*k*R_CUTOFF, 2*BODE_v), -BODE_n/BODE_v);
    p = pow(k, POWER_INDEX) * T * T;
    //p = pow(k, POWER_INDEX - 0.05*log(k/0.05)) * T * T; //running, alpha=0.05
  }
  else if (POWER_SPECTRUM == 1){ // BBKS
    gamma = OMm * hlittle * pow(E, -OMb - OMb/OMm);
    q = k / (hlittle*gamma);
    T = (log(1.0+2.34*q)/(2.34*q)) *
      pow( 1.0+3.89*q + pow(16.1*q, 2) + pow( 5.46*q, 3) + pow(6.71*q, 4), -0.25);
    p = pow(k, POWER_INDEX) * T * T;
  }
  else if (POWER_SPECTRUM == 2){ // Efstathiou,G., Bond,J.R., and White,S.D.M., MNRAS,258,1P (1992)
    gamma = 0.25;
    aa = 6.4/(hlittle*gamma);
    bb = 3.0/(hlittle*gamma);
    cc = 1.7/(hlittle*gamma);
    p = pow(k, POWER_INDEX) / pow( 1+pow( aa*k + pow(bb*k, 1.5) + pow(cc*k,2), 1.13), 2.0/1.13 );
  }
  else if (POWER_SPECTRUM == 3){ // Peebles, pg. 626
    gamma = OMm * hlittle * pow(E, -OMb - OMb/OMm);
    aa = 8.0 / (hlittle*gamma);
    bb = 4.7 / (hlittle*gamma);
    p = pow(k, POWER_INDEX) / pow(1 + aa*k + bb*k*k, 2);
  }
  else if (POWER_SPECTRUM == 4){ // White, SDM and Frenk, CS, 1991, 379, 52
    gamma = OMm * hlittle * pow(E, -OMb - OMb/OMm);
    aa = 1.7/(hlittle*gamma);
    bb = 9.0/pow(hlittle*gamma, 1.5);
    cc = 1.0/pow(hlittle*gamma, 2);
    p = pow(k, POWER_INDEX) * 19400.0 / pow(1 + aa*k + pow(bb*k, 1.5) + cc*k*k, 2);
  }
  else{
    fprintf(stderr, "No such power spectrum defined: %i\nOutput is bogus.\n", POWER_SPECTRUM);
    p = 0;
  }

  // now get the value of the window function
  kR = k * R;
  if (FILTER == 0){ // top hat
    if ( (kR) < 1.0e-4 ){ w = 1.0; }// w converges to 1 as (kR) -> 0
    else { w = 3.0 * (sin(kR)/pow(kR, 3) - cos(kR)/pow(kR, 2));}

    // now do d(w^2)/dm = 2 w dw/dr dr/dm
    if ( (kR) < 1.0e-10 ){  dwdr = 0;}
    else{ dwdr = 9*cos(kR)*k/pow(kR,3) + 3*sin(kR)*(1 - 3/(kR*kR))/(kR*R);}
	    //3*k*( 3*cos(kR)/pow(kR,3) + sin(kR)*(-3*pow(kR, -4) + 1/(kR*kR)) );}
      //     dwdr = -1e8 * k / (R*1e3);
    drdm = 1.0 / (4.0*PI * OMm*RHOcrit * R*R);
  }
  else if (FILTER == 1){ // gaussian of width 1/R
    w = pow(E, -kR*kR/2.0);
    dwdr = - k*kR * w;
    drdm = 1.0 / (pow(2*PI, 1.5) * OMm*RHOcrit * 3*R*R);
  }
  else {
    fprintf(stderr, "No such filter: %i\nOutput is bogus.\n", FILTER);
    w=0;
  }

  //  printf("%e\t%e\t%e\t%e\t%e\t%e\t%e\n", k, R, p, w, dwdr, drdm, dsigmadk[1]);
  return k*k*p*2*w*dwdr*drdm * d2fact;
}
double dsigmasqdm_z0(double M){
  double result, error, lower_limit, upper_limit;
  gsl_function F;
  double rel_tol  = FRACT_FLOAT_ERR*10; //<- relative tolerance
  gsl_integration_workspace * w
    = gsl_integration_workspace_alloc (1000);
  double kstart, kend;

  R = MtoR(M);

  // now lets do the integral for sigma and scale it with sigma_norm
  kstart = 1.0e-99/R;
  kend = 350.0/R;
  lower_limit = kstart;//log(kstart);
  upper_limit = kend;//log(kend);
  d2fact = M*10000/sigma_z0(M);

  F.function = &dsigmasq_dm;
  gsl_integration_qag (&F, lower_limit, upper_limit, 0, rel_tol,
		       1000, GSL_INTEG_GAUSS61, w, &result, &error);
  gsl_integration_workspace_free (w);

  return sigma_norm * sigma_norm * result /d2fact;
}

/*
  FUNCTION TFmdm is the power spectrum transfer function from Eisenstein & Hu ApJ, 1999, 511, 5
*/
double TFmdm(double k){
  double q, gamma_eff, q_eff, TF_m, q_nu;

  q = k*pow(theta_cmb,2)/omhh;
  gamma_eff=sqrt(alpha_nu) + (1.0-sqrt(alpha_nu))/(1.0+pow(0.43*k*sound_horizon, 4));
  q_eff = q/gamma_eff;
  TF_m= log(E+1.84*beta_c*sqrt(alpha_nu)*q_eff);
  TF_m /= TF_m + pow(q_eff,2) * (14.4 + 325.0/(1.0+60.5*pow(q_eff,1.11)));
  q_nu = 3.92*q/sqrt(f_nu/N_nu);
  TF_m *= 1.0 + (1.2*pow(f_nu,0.64)*pow(N_nu,0.3+0.6*f_nu)) /
    (pow(q_nu,-1.6)+pow(q_nu,0.8));

  //   printf("%f  %e  %f  %f  %f  %f\n",omhh,f_nu,f_baryon,N_nu,y_d,alpha_nu);
  // printf("%f  %f  %f  %f\n", beta_c,sound_horizon,theta_cmb,z_equality);
  //printf("%f  %e  %f  %f  %f\n\n",q, k, gamma_eff, q_nu, TF_m);

  return TF_m;
}

void TFset_parameters(){
  double z_drag, R_drag, R_equality, p_c, p_cb, f_c, f_cb, f_nub, k_equality;

  z_equality = 25000*omhh*pow(theta_cmb, -4) - 1.0;
  k_equality = 0.0746*omhh/(theta_cmb*theta_cmb);

  z_drag = 0.313*pow(omhh,-0.419) * (1 + 0.607*pow(omhh, 0.674));
  z_drag = 1 + z_drag*pow(OMb*hlittle*hlittle, 0.238*pow(omhh, 0.223));
  z_drag *= 1291 * pow(omhh, 0.251) / (1 + 0.659*pow(omhh, 0.828));

  y_d = (1 + z_equality) / (1.0 + z_drag);

  R_drag = 31.5 * OMb*hlittle*hlittle * pow(theta_cmb, -4) * 1000 / (1.0 + z_drag);
  R_equality = 31.5 * OMb*hlittle*hlittle * pow(theta_cmb, -4) * 1000 / (1.0 + z_equality);

  sound_horizon = 2.0/3.0/k_equality * sqrt(6.0/R_equality) *
    log( (sqrt(1+R_drag) + sqrt(R_drag+R_equality)) / (1.0 + sqrt(R_equality)) );

  p_c = -(5 - sqrt(1 + 24*(1 - f_nu-f_baryon)))/4.0;
  p_cb = -(5 - sqrt(1 + 24*(1 - f_nu)))/4.0;
  f_c = 1 - f_nu - f_baryon;
  f_cb = 1 - f_nu;
  f_nub = f_nu+f_baryon;

  alpha_nu = (f_c/f_cb) * (2*(p_c+p_cb)+5)/(4*p_cb+5.0);
  alpha_nu *= 1 - 0.553*f_nub+0.126*pow(f_nub,3);
  alpha_nu /= 1-0.193*sqrt(f_nu)+0.169*f_nu;
  alpha_nu *= pow(1+y_d, p_c-p_cb);
  alpha_nu *= 1+ (p_cb-p_c)/2.0 * (1.0+1.0/(4.0*p_c+3.0)/(4.0*p_cb+7.0))/(1.0+y_d);
  beta_c = 1.0/(1.0-0.949*f_nub);
}

double init_ps(){
  double result, error, lower_limit, upper_limit;
  gsl_function F;
  double rel_tol  = FRACT_FLOAT_ERR*10; //<- relative tolerance
  gsl_integration_workspace * w
    = gsl_integration_workspace_alloc (1000);
  double kstart, kend;
  int i;
  double x;

  // Set cuttoff scale for WDM (eq. 4 in Barkana et al. 2001) in comoving Mpc
  R_CUTOFF = 0.201*pow((OMm-OMb)*hlittle*hlittle/0.15, 0.15)*pow(g_x/1.5, -0.29)*pow(M_WDM, -1.15);

//  fprintf(stderr, "For M_DM = %.2e keV, R_CUTOFF is: %.2e comoving Mpc\n", M_WDM, R_CUTOFF);
  if (!P_CUTOFF)
//    fprintf(stderr, "But you have selected CDM, so this is ignored\n");

  omhh = OMm*hlittle*hlittle;
  theta_cmb = T_cmb / 2.7;

  // Translate Parameters into forms GLOBALVARIABLES form
  f_nu = OMn/OMm;
  f_baryon = OMb/OMm;
  if (f_nu < TINY) f_nu = 1e-10;
  if (f_baryon < TINY) f_baryon = 1e-10;

  TFset_parameters();

  sigma_norm = -1;

  R = 8.0/hlittle;
  kstart = 1.0e-99/R;
  kend = 350.0/R;
  lower_limit = kstart;//log(kstart);
  upper_limit = kend;//log(kend);

  F.function = &dsigma_dk;

  gsl_integration_qag (&F, lower_limit, upper_limit, 0, rel_tol,
		       1000, GSL_INTEG_GAUSS61, w, &result, &error);
  gsl_integration_workspace_free (w);

  sigma_norm = SIGMA8/sqrt(result); //takes care of volume factor

  /* initialize the lookup table for erfc */
  /*
  for (i=0; i<=ERFC_NPTS; i++){
    erfc_params[i] = i*ERFC_PARAM_DELTA;
    log_erfc_table[i] = log(erfcc(erfc_params[i]));
  }
  // Set up spline table
  erfc_acc   = gsl_interp_accel_alloc ();
  erfc_spline  = gsl_spline_alloc (gsl_interp_cspline, ERFC_NPTS);
  gsl_spline_init(erfc_spline, erfc_params, log_erfc_table, ERFC_NPTS);
  */

  return R_CUTOFF;
}

void free_ps(){
  /*    gsl_spline_free (erfc_spline);
    gsl_interp_accel_free(erfc_acc);
  */
  return;
}

double splined_erfc(double x){
  if (x < 0){
    //printf("WARNING: Negative value %e passed to splined_erfc. Returning 1\n", x);
    return 1;
  }
  return erfcc(x); // the interpolation below doesn't seem to be stable in Ts.c
  if (x > ERFC_PARAM_DELTA*(ERFC_NPTS-1))
    return erfcc(x);
  else
    return exp(gsl_spline_eval(erfc_spline, x, erfc_acc));
}

#endif


text:
//////////////////////////////////////////////////////////////////////////////////////
// This file is distributed under the University of Illinois/NCSA Open Source License.
// See LICENSE file in top directory for details.
//
// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
//
// File developed by: Paul R. C. Kent, kentpr@ornl.gov, Oak Ridge National Laboratory
//                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Paul R. C. Kent, kentpr@ornl.gov, Oak Ridge National Laboratory
//////////////////////////////////////////////////////////////////////////////////////

//           http://pathintegrals.info                     //
/////////////////////////////////////////////////////////////

#ifndef VECTOR_OPS_H
#define VECTOR_OPS_H

// extern "C"{
// #ifdef USE_MKL
//   #include <mkl_cblas.h>
// #else
//   #include <cblas.h>
// #endif
// }
#include "Blitz.h"
#include "../config.h"

typedef TinyVector<int, 3> Int3;
typedef Array<std::complex<double>, 1> zVec;
typedef Array<cVec3, 1> zVecVec;
typedef Array<cMat3, 1> zMatVec;

#ifdef USE_CBLAS
extern "C"
{
#ifdef USE_MKL_CBLAS
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
}
#else

#define F77_DZNRM2 F77_FUNC(dznrm2, DZNRM2)
#define F77_ZDSCAL F77_FUNC(zdscal, ZDSCAL)
#define F77_ZDOTC F77_FUNC(zdotc, ZDOTC)
#define F77_ZGEMV F77_FUNC(zgemv, ZGEMV)
#define F77_ZAXPY F77_FUNC(zaxpy, ZAXPY)

extern "C" double F77_DZNRM2(const int* N, const void* X, const int* INC);
extern "C" void F77_ZDSCAL(const int* N, double* ALPHA, const void* X, const int* INC);
// extern "C" void F77_ZDOTC (std::complex<double> *z, const int *N,
// 			   const void *X, const int *INCX,
// 			   const void *Y, const int *INCY);
extern "C" std::complex<double> F77_ZDOTC(const int* N, const void* X, const int* INCX, const void* Y, const int* INCY);
extern "C" void F77_ZGEMV(char* TRANS,
                          const int* M,
                          const int* N,
                          std::complex<double>* alpha,
                          const void* A,
                          const int* LDA,
                          const void* X,
                          const int* INCX,
                          std::complex<double>* beta,
                          const void* Y,
                          const int* INCY);

extern "C" void F77_ZAXPY(const int* N, std::complex<double>* ALPHA, void* X, int* INCX, void* Y, int* INCY);
#endif

#ifdef USE_CBLAS
inline void Normalize(zVec& c)
{
  double norm = cblas_dznrm2(c.size(), c.data(), 1);
  norm        = 1.0 / norm;
  cblas_zdscal(c.size(), norm, c.data(), 1);
}

inline double norm(const zVec& c) { return cblas_dznrm2(c.size(), c.data(), 1); }

inline std::complex<double> conjdot(zVec& cA, zVec& cB)
{
  std::complex<double> z;
  cblas_zdotc_sub(cA.size(), cA.data(), 1, cB.data(), 1, &z);
  return z;
}

inline void Orthogonalize(const Array<std::complex<double>, 2>& A, zVec& x)
{
  int m = A.rows();
  int n = A.cols();
  assert(n == x.size());
  std::complex<double> zero(0.0, 0.0);
  std::complex<double> one(1.0, 0.0);
  std::complex<double> minusone(-1.0, 0.0);
  Array<std::complex<double>, 1> S(m);

  cblas_zgemv(CblasColMajor, CblasConjTrans, n, m, &one, A.data(), n, x.data(), 1, &zero, S.data(), 1);
  cblas_zgemv(CblasRowMajor, CblasTrans, m, n, &minusone, A.data(), n, S.data(), 1, &one, x.data(), 1);
}

#else

inline void Normalize(zVec& c)
{
  const int inc = 1;
  int n         = c.size();
  double norm   = F77_DZNRM2(&n, c.data(), &inc);
  norm          = 1.0 / norm;
  F77_ZDSCAL(&n, &norm, c.data(), &inc);
}

inline double norm(const zVec& c)
{
  double norm;
  int n         = c.size();
  const int inc = 1;
  return F77_DZNRM2(&n, c.data(), &inc);
}

inline std::complex<double> conjdot(zVec& cA, zVec& cB)
{
  const int n    = cA.size();
  const int incA = 1;
  const int incB = 1;
  return F77_ZDOTC(&n, cA.data(), &incA, cB.data(), &incB);

  //   std::complex<double> z;
  //   F77_ZDOTC (&z, &n, cA.data(), &incA, cB.data(), &incB);
  //   return z;
}

inline void Orthogonalize(const Array<std::complex<double>, 2>& A, zVec& x)
{
  int m = A.rows();
  int n = A.cols();
  assert(n == x.size());
  std::complex<double> zero(0.0, 0.0);
  std::complex<double> one(1.0, 0.0);
  std::complex<double> minusone(-1.0, 0.0);
  Array<std::complex<double>, 1> S(m);

  // Calculate overlaps
  // Calling with column major and ConjTrans is equivalent to
  // conjugate of untransposed row major
  char trans    = 'C';
  const int inc = 1;

  F77_ZGEMV(&trans, &n, &m, &one, A.data(), &n, x.data(), &inc, &zero, S.data(), &inc);

  //   cblas_zgemv(CblasColMajor, CblasConjTrans, n, m, &one,
  // 	      A.data(), n, x.data(), 1, &zero, S.data(), 1);

  //   for (int i=0; i<m; i++) {
  //     fprintf (stderr, "S[%d] = %18.14f + %18.14fi\n",
  // 	     real(S[i]), imag(S[i]));

  // Now, subtract off components * overlaps
  trans = 'T';
  F77_ZGEMV(&trans, &m, &n, &minusone, A.data(), &n, S.data(), &inc, &one, x.data(), &inc);
  //   cblas_zgemv(CblasRowMajor, CblasTrans, m, n, &minusone,
  //  	      A.data(), n, S.data(), 1, &one, x.data(), 1);
}

#endif

inline double realconjdot(zVec& cA, zVec& cB) { return conjdot(cA, cB).real(); }

inline double mag(std::complex<double> x) { return (x.real() * x.real() + x.imag() * x.imag()); }

inline void Orthogonalize2(Array<std::complex<double>, 2>& A, zVec& x, int lastBand)
{
  int m = A.rows();
  int n = A.cols();
  assert(n == x.size());
  zVec Ar;
  Array<std::complex<double>, 1> S(m);

  for (int row = 0; row <= lastBand; row++)
  {
    Ar.reference(A(row, Range::all()));
    S(row) = conjdot(Ar, x);
  }
  for (int row = 0; row <= lastBand; row++)
    x -= S(row) * A(row, Range::all());
  //   for (int row=0; row<=lastBand; row++) {
  //     Ar.reference (A(row,Range::all()));
  //     S(row) = conjdot (Ar, x);
  //     if (mag(S(row)) > 1.0e-14) {
  //       std::cerr << "row = " << row << " lastband = " << lastBand << std::endl;
  //       std::cerr << "Error in Orthogonalize2!, s = " << S(row) << std::endl;
  //       double norm = realconjdot (Ar, Ar);
  //       std::cerr << "norm = " << norm << std::endl;
  //     }
  //   }
}

inline void OrthogExcluding(const Array<std::complex<double>, 2>& A, zVec& x, int excluding)
{
  int m = A.rows();
  int n = A.cols();
  assert(n == x.size());
  zVec Ar;

#ifndef __INTEL_COMPILER
  std::complex<double> S[m];
  for (int row = 0; row < m; row++)
  {
    Ar.reference(A(row, Range::all()));
    S[row] = conjdot(Ar, x);
  }
  for (int row = 0; row < m; row++)
    if (row != excluding)
      x -= S[row] * A(row, Range::all());
#else

  double Sre[m], Sim[m];

  for (int row = 0; row < m; row++)
  {
    Ar.reference(A(row, Range::all()));
    std::complex<double> S = conjdot(Ar, x);
    Sre[row]               = real(S);
    Sim[row]               = imag(S);
  }
  for (int row = 0; row < m; row++)
    if (row != excluding)
      x -= std::complex<double>(Sre[row], Sim[row]) * A(row, Range::all());
#endif
}

inline void OrthogLower(const Array<std::complex<double>, 2>& A, zVec& x, int currBand)
{
  int m = currBand;
  int n = A.cols();
  assert(n == x.size());
  zVec Ar;
  std::complex<double> S;

  for (int row = 0; row < m; row++)
  {
    Ar.reference(A(row, Range::all()));
    S = conjdot(Ar, x);
    x -= S * A(row, Range::all());
  }
}

inline void GramSchmidt(Array<std::complex<double>, 2>& A)
{
  zVec a, b;
  for (int i = 0; i < A.rows(); i++)
  {
    a.reference(A(i, Range::all()));
    Normalize(a);
    for (int j = i + 1; j < A.rows(); j++)
    {
      b.reference(A(j, Range::all()));
      b = b - (conjdot(a, b) * a);
      Normalize(b);
    }
  }
}

inline void Orthogonalize(Array<std::complex<double>, 2>& A)
{
  zVec x, y;
  Array<std::complex<double>, 1> S(A.rows());
  for (int iter = 0; iter < 40; iter++)
  {
    for (int i = 0; i < A.rows(); i++)
    {
      x.reference(A(i, Range::all()));
      for (int j = i + 1; j < A.rows(); j++)
      {
        y.reference(A(j, Range::all()));
        S(j) = conjdot(y, x);
      }
      for (int j = i + 1; j < A.rows(); j++)
      {
        y.reference(A(j, Range::all()));
        x -= S(j) * y;
      }
      Normalize(x);
    }
  }
}

inline void CheckOrthog(const Array<std::complex<double>, 2>& A, zVec& x)
{
  zVec Ai;
  double normInv = 1.0 / norm(x);
  for (int i = 0; i < A.rows(); i++)
  {
    Ai.reference(A(i, Range::all()));
    if (normInv * mag(conjdot(Ai, x)) > 1.0e-13)
    {
      std::cerr << "CheckOrthog failed for i=" << i << ".\n";
      exit(1);
    }
  }
}

inline void zaxpy(std::complex<double> alpha,
                  const Array<std::complex<double>, 1>& x,
                  const std::complex<double> y,
                  Array<std::complex<double>, 1>& axpy)
{}

// inline Array<std::complex<double>,3>&
// operator*= (Array<std::complex<double>,3> &A, const Array<std::complex<double>,3> &B)
// {

// }

#endif

ext: h

max_stars_repo_name: djstaros/qmcpack

max_issues_repo_name: djstaros/qmcpack

max_forks_repo_name: djstaros/qmcpack


text:
#ifndef __ESTIMATE_INC_
#define __ESTIMATE__INC_

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

double getLogLikelyhood(gsl_matrix *cinverse, double det_cmatrix,  gsl_matrix *xmodel,
												gsl_vector *trainingvector, gsl_vector *thetas, gsl_matrix *h_matrix, int nmodel_points,
												int nthetas, int nparams, int nregression_fns,
												void (*makeHVector)(gsl_vector *, gsl_vector*, int));

#endif

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: MADAI/MADAIEmulator

max_forks_repo_name: MADAI/MADAIEmulator


text:
/* vector/gsl_vector_long.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, 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.
 */

#ifndef __GSL_VECTOR_LONG_H__
#define __GSL_VECTOR_LONG_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_long.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  long *data;
  gsl_block_long *block;
  int owner;
}
gsl_vector_long;

typedef struct
{
  gsl_vector_long vector;
} _gsl_vector_long_view;

typedef _gsl_vector_long_view gsl_vector_long_view;

typedef struct
{
  gsl_vector_long vector;
} _gsl_vector_long_const_view;

typedef const _gsl_vector_long_const_view gsl_vector_long_const_view;

/* Allocation */

gsl_vector_long *gsl_vector_long_alloc (const size_t n);
gsl_vector_long *gsl_vector_long_calloc (const size_t n);

gsl_vector_long *gsl_vector_long_alloc_from_block (gsl_block_long * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

gsl_vector_long *gsl_vector_long_alloc_from_vector (gsl_vector_long * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

void gsl_vector_long_free (gsl_vector_long * v);

/* Views */

_gsl_vector_long_view
gsl_vector_long_view_array (long *v, size_t n);

_gsl_vector_long_view
gsl_vector_long_view_array_with_stride (long *base,
                                         size_t stride,
                                         size_t n);

_gsl_vector_long_const_view
gsl_vector_long_const_view_array (const long *v, size_t n);

_gsl_vector_long_const_view
gsl_vector_long_const_view_array_with_stride (const long *base,
                                               size_t stride,
                                               size_t n);

_gsl_vector_long_view
gsl_vector_long_subvector (gsl_vector_long *v,
                            size_t i,
                            size_t n);

_gsl_vector_long_view
gsl_vector_long_subvector_with_stride (gsl_vector_long *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

_gsl_vector_long_const_view
gsl_vector_long_const_subvector (const gsl_vector_long *v,
                                  size_t i,
                                  size_t n);

_gsl_vector_long_const_view
gsl_vector_long_const_subvector_with_stride (const gsl_vector_long *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

void gsl_vector_long_set_zero (gsl_vector_long * v);
void gsl_vector_long_set_all (gsl_vector_long * v, long x);
int gsl_vector_long_set_basis (gsl_vector_long * v, size_t i);

int gsl_vector_long_fread (FILE * stream, gsl_vector_long * v);
int gsl_vector_long_fwrite (FILE * stream, const gsl_vector_long * v);
int gsl_vector_long_fscanf (FILE * stream, gsl_vector_long * v);
int gsl_vector_long_fprintf (FILE * stream, const gsl_vector_long * v,
                              const char *format);

int gsl_vector_long_memcpy (gsl_vector_long * dest, const gsl_vector_long * src);

int gsl_vector_long_reverse (gsl_vector_long * v);

int gsl_vector_long_swap (gsl_vector_long * v, gsl_vector_long * w);
int gsl_vector_long_swap_elements (gsl_vector_long * v, const size_t i, const size_t j);

long gsl_vector_long_max (const gsl_vector_long * v);
long gsl_vector_long_min (const gsl_vector_long * v);
void gsl_vector_long_minmax (const gsl_vector_long * v, long * min_out, long * max_out);

size_t gsl_vector_long_max_index (const gsl_vector_long * v);
size_t gsl_vector_long_min_index (const gsl_vector_long * v);
void gsl_vector_long_minmax_index (const gsl_vector_long * v, size_t * imin, size_t * imax);

int gsl_vector_long_add (gsl_vector_long * a, const gsl_vector_long * b);
int gsl_vector_long_sub (gsl_vector_long * a, const gsl_vector_long * b);
int gsl_vector_long_mul (gsl_vector_long * a, const gsl_vector_long * b);
int gsl_vector_long_div (gsl_vector_long * a, const gsl_vector_long * b);
int gsl_vector_long_scale (gsl_vector_long * a, const double x);
int gsl_vector_long_add_constant (gsl_vector_long * a, const double x);

int gsl_vector_long_equal (const gsl_vector_long * u,
                            const gsl_vector_long * v);

int gsl_vector_long_isnull (const gsl_vector_long * v);
int gsl_vector_long_ispos (const gsl_vector_long * v);
int gsl_vector_long_isneg (const gsl_vector_long * v);
int gsl_vector_long_isnonneg (const gsl_vector_long * v);

INLINE_DECL long gsl_vector_long_get (const gsl_vector_long * v, const size_t i);
INLINE_DECL void gsl_vector_long_set (gsl_vector_long * v, const size_t i, long x);
INLINE_DECL long * gsl_vector_long_ptr (gsl_vector_long * v, const size_t i);
INLINE_DECL const long * gsl_vector_long_const_ptr (const gsl_vector_long * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
long
gsl_vector_long_get (const gsl_vector_long * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

INLINE_FUN
void
gsl_vector_long_set (gsl_vector_long * v, const size_t i, long x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

INLINE_FUN
long *
gsl_vector_long_ptr (gsl_vector_long * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (long *) (v->data + i * v->stride);
}

INLINE_FUN
const long *
gsl_vector_long_const_ptr (const gsl_vector_long * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const long *) (v->data + i * v->stride);
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_LONG_H__ */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: ielomariala/Hex-Game


text:
#ifndef _RICO_
#define _RICO_

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <ctype.h>
#include <math.h>
#include <unistd.h>
#include "../Parameter_files/COSMOLOGY.H"
#include "../Parameter_files/INIT_PARAMS.H"
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include "cosmo_progs.c"
#include "misc.c"
#include "ps.c"

#define A_NPTS (int) (60) /*Warning: the calculation of the MHR model parameters is valid only from redshift 2 to A_NPTS+2*/
static double A_table[A_NPTS], A_params[A_NPTS];
static gsl_interp_accel *A_acc;
static gsl_spline *A_spline;

#define C_NPTS (int) (12)
static double C_table[C_NPTS], C_params[C_NPTS];
static gsl_interp_accel *C_acc;
static gsl_spline *C_spline;

#define beta_NPTS (int) (5)
static double beta_table[beta_NPTS], beta_params[beta_NPTS];
static gsl_interp_accel *beta_acc;
static gsl_spline *beta_spline;

#define RR_Z_NPTS (int) (300) // number of points in redshift axis;  we will only interpolate over gamma, and just index sample in redshift
#define RR_DEL_Z (float) (0.2)
#define RR_lnGamma_NPTS (int) (150) // number of samples of gamma for the interpolation tables
#define RR_lnGamma_min (double) (-10) // min ln gamma12 used
#define RR_DEL_lnGamma (float) (0.1)
static double RR_table[RR_Z_NPTS][RR_lnGamma_NPTS], lnGamma_values[RR_lnGamma_NPTS];
static gsl_interp_accel *RR_acc[RR_Z_NPTS];
static gsl_spline *RR_spline[RR_Z_NPTS];

/***  FUNCTION PROTOTYPES ***/
double splined_recombination_rate(double z_eff, double gamma12_bg); // assumes T=1e4 and case B

double recombination_rate(double z_eff, double gamma12_bg, double T4, int usecaseB);
void init_MHR(); /*initializes the lookup table for the PDF density integral in MHR00 model at redshift z*/
void free_MHR(); /* deallocates the gsl structures from init_MHR */
double Gamma_SS(double Gamma_bg, double Delta, double T_4, double z);//ionization rate w. self shielding
double MHR_rr (double D, void *params);
double A_MHR(double z); /*returns the A parameter in MHR00model*/
double C_MHR(double z); /*returns the C parameter in MHR00model*/
double beta_MHR(double z); /*returns the beta parameter in MHR00model*/
double splined_A_MHR(double z); /*returns the splined A parameter in MHR00model*/
double splined_C_MHR(double z); /*returns the splined C parameter in MHR00model*/
double splined_beta_MHR(double z);/*returns the splined beta parameter in MHR00*/
void free_A_MHR(); /* deallocates the gsl structures from init_A */
void free_C_MHR(); /* deallocates the gsl structures from init_C */
void free_beta_MHR(); /* deallocates the gsl structures from init_beta */
void init_A_MHR(); /*initializes the lookup table for the A paremeter in MHR00 model*/
void init_C_MHR(); /*initializes the lookup table for the C paremeter in MHR00 model*/
void init_beta_MHR(); /*initializes the lookup table for the beta paremeter in MHR00 model*/

double splined_recombination_rate(double z_eff, double gamma12_bg){
  int z_ct = (int) (z_eff / RR_DEL_Z + 0.5); // round to nearest int
  double lnGamma = log(gamma12_bg);

  // check out of bounds
  if ( z_ct < 0 ){ // out of array bounds
//    fprintf(stderr, "WARNING: splined_recombination_rate: effective redshift %g is outside of array bouds\n", z_eff);
    z_ct = 0;
  }
  else if (z_ct  >= RR_Z_NPTS){
//    fprintf(stderr, "WARNING: splined_recombination_rate: effective redshift %g is outside of array bouds\n", z_eff);
    z_ct = RR_Z_NPTS-1;
  }

  if (lnGamma < RR_lnGamma_min){
    return 0;
  }
  else if (lnGamma >= (RR_lnGamma_min + RR_DEL_lnGamma * RR_lnGamma_NPTS) ){
//    fprintf(stderr, "WARNING: splined_recombination_rate: Gamma12 of %g is outside of interpolation array\n", gamma12_bg);
    lnGamma =  RR_lnGamma_min + RR_DEL_lnGamma * RR_lnGamma_NPTS - FRACT_FLOAT_ERR;
  }

  return gsl_spline_eval(RR_spline[z_ct], lnGamma, RR_acc[z_ct]);
}

void init_MHR(){
  int z_ct, gamma_ct;
  float z, gamma;

  // first initialize the MHR parameter look up tables
  init_C_MHR(); /*initializes the lookup table for the C paremeter in MHR00 model*/
  init_beta_MHR(); /*initializes the lookup table for the beta paremeter in MHR00 model*/
  init_A_MHR(); /*initializes the lookup table for the A paremeter in MHR00 model*/

  // now the recombination rate look up tables
  for (z_ct=0; z_ct < RR_Z_NPTS; z_ct++){

    z = z_ct * RR_DEL_Z; // redshift corresponding to index z_ct of the array

    // Intialize the Gamma values
    for (gamma_ct=0; gamma_ct < RR_lnGamma_NPTS; gamma_ct++){
      lnGamma_values[gamma_ct] = RR_lnGamma_min  + gamma_ct*RR_DEL_lnGamma;  // ln of Gamma12
      gamma = exp(lnGamma_values[gamma_ct]);
      RR_table[z_ct][gamma_ct] = recombination_rate(z, gamma, 1, 1); // CHANGE THIS TO INCLUDE TEMPERATURE
    }

    // set up the spline in gamma
    RR_acc[z_ct] = gsl_interp_accel_alloc();
    RR_spline[z_ct] = gsl_spline_alloc (gsl_interp_cspline, RR_lnGamma_NPTS);
    gsl_spline_init(RR_spline[z_ct], lnGamma_values, RR_table[z_ct], RR_lnGamma_NPTS);

  } // go to next redshift

  return;
}

void free_MHR(){
  int z_ct;

  free_A_MHR();
  free_C_MHR();
  free_beta_MHR();

  // now the recombination rate look up tables
  for (z_ct=0; z_ct < RR_Z_NPTS; z_ct++){
    gsl_spline_free (RR_spline[z_ct]);
    gsl_interp_accel_free(RR_acc[z_ct]);
  }

  return;
}

//calculates the attenuated photoionization rate due to self-shielding (in units of 1e-12 s^-1)
// input parameters are the background ionization rate, overdensity, temperature (in 10^4k), redshift, respectively
//  Uses the fitting formula from Rahmati et al, assuming a UVB power law index of alpha=5
double Gamma_SS(double Gamma_bg, double D, double T_4, double z){
  double D_ss = 26.7*pow(T_4, 0.17) * pow( (1+z)/10.0, -3) * pow(Gamma_bg, 2.0/3.0);
  return Gamma_bg * (0.98 * pow( (1.0+pow(D/D_ss, 1.64)), -2.28) + 0.02*pow( 1.0+D/D_ss, -0.84));
}

typedef struct {double z, gamma12_bg, T4, A, C_0, beta, avenH; int usecaseB;} RR_par;

double MHR_rr (double lnD, void *params){
  double D=exp(lnD);
  double alpha;
  RR_par p = *(RR_par *) params;
  double z = p.z;
  double gamma = Gamma_SS(p.gamma12_bg, D, p.T4, z);
  double n_H = p.avenH*D;
  double x_e = 1.0 - neutral_fraction(n_H, p.T4, gamma, p.usecaseB);
  double PDelta;

  PDelta = p.A * exp( - 0.5*pow((pow(D,-2.0/3.0)- p.C_0 ) / ((2.0*7.61/(3.0*(1.0+z)))), 2)) * pow(D, p.beta);

  if (p.usecaseB)
    alpha = alpha_B(p.T4*1e4);
  else
    alpha = alpha_A(p.T4*1e4);

  //  fprintf(stderr, "%g\t%g\t%g\t%g\t%g\n", n_H, PDelta, alpha, x_e, D);

  return n_H * PDelta * alpha * x_e * x_e * D * D;//note extra D since we are integrating over lnD
}

// returns the recombination rate per baryon (1/s), integrated over the MHR density PDF,
// given an ionizing background of gamma12_bg
// temeperature T4 (in 1e4 K), and usecaseB rate coefficient
// Assumes self-shielding according to Rahmati+ 2013
double recombination_rate(double z, double gamma12_bg, double T4, int usecaseB){
  double result, error, lower_limit, upper_limit, A, C_0, beta, avenH;
  gsl_function F;
  double rel_tol  = 0.01; //<- relative tolerance
  gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
  RR_par p = {z, gamma12_bg, T4, A_MHR(z), C_MHR(z), beta_MHR(z), No*pow( 1+z, 3), usecaseB};

  F.function = &MHR_rr;
  F.params=&p;
  lower_limit = log(0.01);
  upper_limit = log(200);

  gsl_integration_qag (&F, lower_limit, upper_limit, 0, rel_tol,
		       1000, GSL_INTEG_GAUSS61, w, &result, &error);
  gsl_integration_workspace_free (w);

  return result;
}

double aux_function(double D, void *params){
  double result;
  double z = *(double *) params;

  result = exp(-(pow(D,-2.0/3.0)-C_MHR(z))*(pow(D,-2.0/3.0)-C_MHR(z))/(2.0*(2.0*7.61/(3.0*(1.0+z)))*(2.0*7.61/(3.0*(1.0+z)))))*pow(D, beta_MHR(z));

  return result;
}

double A_aux_integral(double z){
  double result, error, lower_limit, upper_limit;
  gsl_function F;
  double rel_tol  = 0.001; //<- relative tolerance
  gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

  F.function = &aux_function;
  F.params = &z;
  lower_limit = 1e-25;
  upper_limit = 1e25;

  gsl_integration_qag (&F, lower_limit, upper_limit, 0, rel_tol,
		       1000, GSL_INTEG_GAUSS61, w, &result, &error);
  gsl_integration_workspace_free (w);

  return result;
}

double A_MHR(double z){
  double result;
  if(z>=2.0+(float)A_NPTS)
    result = splined_A_MHR(2.0+(float)A_NPTS);
  else
    if(z<=2.0)
      result = splined_A_MHR(2.0);
    else
      result = splined_A_MHR(z);
  return result;
}

void init_A_MHR(){
/* initialize the lookup table for the parameter A in the MHR00 model */
   int i;

   for (i=0; i<A_NPTS; i++){
     A_params[i] = 2.0+(float)i;
     A_table[i] = 1.0/A_aux_integral(2.0+(float)i);
   }

  // Set up spline table
  A_acc   = gsl_interp_accel_alloc();
  A_spline  = gsl_spline_alloc (gsl_interp_cspline, A_NPTS);
  gsl_spline_init(A_spline, A_params, A_table, A_NPTS);

  return;
 }

double splined_A_MHR(double x){
  return gsl_spline_eval(A_spline, x, A_acc);
}

void free_A_MHR(){

  gsl_spline_free (A_spline);
  gsl_interp_accel_free(A_acc);

  return;
}

double C_MHR(double z){
  double result;
  if(z>=13.0)
    result = 1.0;
  else
    if(z<=2.0)
      result = 0.558;
    else
      result = splined_C_MHR(z);
  return result;
}

void init_C_MHR(){
/* initialize the lookup table for the parameter C in the MHR00 model */
   int i;

  for (i=0; i<C_NPTS; i++)
    C_params[i] = (float)i+2.0;

  C_table[0] = 0.558;
  C_table[1] = 0.599;
  C_table[2] = 0.611;
  C_table[3] = 0.769;
  C_table[4] = 0.868;
  C_table[5] = 0.930;
  C_table[6] = 0.964;
  C_table[7] = 0.983;
  C_table[8] = 0.993;
  C_table[9] = 0.998;
  C_table[10] = 0.999;
  C_table[11] = 1.00;

  // Set up spline table
  C_acc   = gsl_interp_accel_alloc ();
  C_spline  = gsl_spline_alloc (gsl_interp_cspline, C_NPTS);
  gsl_spline_init(C_spline, C_params, C_table, C_NPTS);

  return;
 }

double splined_C_MHR(double x){
  return gsl_spline_eval(C_spline, x, C_acc);
}

void free_C_MHR(){

  gsl_spline_free (C_spline);
  gsl_interp_accel_free(C_acc);

  return;
}

double beta_MHR(double z){
  double result;
  if(z>=6.0)
    result = -2.50;
  else
    if(z<=2.0)
      result = -2.23;
    else
      result = splined_beta_MHR(z);
  return result;
}

void init_beta_MHR(){
/* initialize the lookup table for the parameter C in the MHR00 model */
   int i;

  for (i=0; i<beta_NPTS; i++)
    beta_params[i] = (float)i+2.0;

  beta_table[0] = -2.23;
  beta_table[1] = -2.35;
  beta_table[2] = -2.48;
  beta_table[3] = -2.49;
  beta_table[4] = -2.50;

  // Set up spline table
  beta_acc   = gsl_interp_accel_alloc ();
  beta_spline  = gsl_spline_alloc (gsl_interp_cspline, beta_NPTS);
  gsl_spline_init(beta_spline, beta_params, beta_table, beta_NPTS);

  return;
 }

double splined_beta_MHR(double x){
  return gsl_spline_eval(beta_spline, x, beta_acc);
}

void free_beta_MHR(){

  gsl_spline_free(beta_spline);
  gsl_interp_accel_free(beta_acc);

  return;
}

#endif

ext: c


text:
/* multilarge_nlinear/test.c
 *
 * Copyright (C) 2015, 2016 Patrick Alken
 *
 * 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.
 */

/* These tests are based on the NIST Statistical Reference Datasets
   See http://www.nist.gov/itl/div898/strd/index.html for more
   information. */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_multilarge_nlinear.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_ieee_utils.h>

#include "test_fdf.c"

static void
test_proc(const gsl_multilarge_nlinear_trs *trs,
          const gsl_multilarge_nlinear_scale *scale,
          const int fdtype)
{
  gsl_multilarge_nlinear_parameters fdf_params =
    gsl_multilarge_nlinear_default_parameters();

  fdf_params.trs = trs;
  fdf_params.scale = scale;
  fdf_params.fdtype = fdtype;

  if (trs == gsl_multilarge_nlinear_trs_lm ||
      trs == gsl_multilarge_nlinear_trs_lmaccel)
    fdf_params.solver = gsl_multilarge_nlinear_solver_cholesky;
  else
    fdf_params.solver = gsl_multilarge_nlinear_solver_mcholesky;

  test_fdf_main(&fdf_params);
}

int
main (void)
{
  const gsl_multilarge_nlinear_trs **nlinear_trs[7];
  const gsl_multilarge_nlinear_scale **nlinear_scales[3];
  const gsl_multilarge_nlinear_trs **trs;
  const gsl_multilarge_nlinear_scale **scale;
  int fdtype;
  size_t i = 0;

  gsl_ieee_env_setup();

  /* initialize arrays */

  nlinear_trs[0] = &gsl_multilarge_nlinear_trs_lm;
  nlinear_trs[1] = &gsl_multilarge_nlinear_trs_lmaccel;
  nlinear_trs[2] = &gsl_multilarge_nlinear_trs_dogleg;
  nlinear_trs[3] = &gsl_multilarge_nlinear_trs_ddogleg;
  nlinear_trs[4] = &gsl_multilarge_nlinear_trs_subspace2D;
  nlinear_trs[5] = &gsl_multilarge_nlinear_trs_cgst;
  nlinear_trs[6] = NULL;

  nlinear_scales[0] = &gsl_multilarge_nlinear_scale_levenberg;
  nlinear_scales[1] = &gsl_multilarge_nlinear_scale_more;
  nlinear_scales[2] = NULL;

  /* run testsuite over all parameter combinations */

  for (trs = nlinear_trs[i]; trs != NULL; trs = nlinear_trs[++i])
    {
      size_t j = 0;

      fprintf(stderr, "trs = %s\n", (*trs)->name);

      for (scale = nlinear_scales[j]; scale != NULL; scale = nlinear_scales[++j])
        {
          for (fdtype = GSL_MULTILARGE_NLINEAR_FWDIFF;
               fdtype <= GSL_MULTILARGE_NLINEAR_CTRDIFF; ++fdtype)
            {
              test_proc(*trs, *scale, fdtype);
            }
        }
    }

  exit (gsl_test_summary ());
}

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for functions windowing and computing FFT/IFFT of time/frequency series.
 *
 */

#ifndef _FFT_H
#define _FFT_H

#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 <fftw3.h> /* Note: when included AFTER complex.h, fftw_complex type defaults to the native double complex */

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

/* Window functions */
double WindowFunction(double x, double xi, double xf, double deltaxi, double deltaxf);
double WindowFunctionLeft(double x, double xf, double deltaxf);
double WindowFunctionRight(double x, double xi, double deltaxi);

/* FFT of real time series */
/* Note: FFT uses flipped convention (i.e. h(f) = int e^(+2ipift)h(t)) */
int FFTRealTimeSeries(
  ReImFrequencySeries** freqseries,   /* Output: frequency series */
  RealTimeSeries* timeseries,         /* Input: real time series */
  double twindowbeg,                  /* Extent of the window at beginning (starts at the first point) */
  double twindowend,                  /* Extent of the window at the end (end at the last point) */
  int nzeropad);                      /* For 0-padding: length will be (upper power of 2)*2^nzeropad */

/* FFT of Re/Im time series */
/* Note: FFT uses flipped convention (i.e. h(f) = int e^(+2ipift)h(t)) */
int FFTTimeSeries(
  ReImFrequencySeries** freqseries,   /* Output: frequency series */
  ReImTimeSeries* timeseries,         /* Input: Re/Im time series */
  double twindowbeg,                  /* Extent of the window at beginning (starts at the first point) */
  double twindowend,                  /* Extent of the window at the end (end at the last point) */
  int nzeropad);                      /* For 0-padding: length will be (upper power of 2)*2^nzeropad */

/* IFFT of frequency series */
/* Note: assumes frequency series is FT of real data */
/* Note: FFT uses flipped convention (i.e. h(f) = int e^(+2ipift)h(t)) */
int IFFTFrequencySeriesReal(
  RealTimeSeries** timeseries,       /* Output: real time series*/
  ReImFrequencySeries* freqseries,   /* Input: complex frequency series, assumed to be the FT of a real time series */
  double f1windowbeg,                /* Start of window at the beginning */
  double f2windowbeg,                /* End of window at the beginning */
  double f1windowend,                /* Start of window at the end */
  double f2windowend,                /* End of window at the end */
  int nzeropad);                     /* For 0-padding: length will be (upper power of 2)*2^nzeropad */

/* IFFT of frequency series */
/* Note: assumes frequency series is FT of complex data - produces complex output */
/* Note: FFT uses flipped convention (i.e. h(f) = int e^(+2ipift)h(t)) */
int IFFTFrequencySeries(
  ReImTimeSeries** timeseries,       /* Output: complex time series */
  ReImFrequencySeries* freqseries,   /* Input: complex frequency series, assumed to be the FT of a real time series */
  double f1windowbeg,                /* Start of window at the beginning */
  double f2windowbeg,                /* End of window at the beginning */
  double f1windowend,                /* Start of window at the end */
  double f2windowend,                /* End of window at the end */
  int nzeropad);                     /* For 0-padding: length will be (upper power of 2)*2^nzeropad */

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _FFT_H */

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
#ifndef AMGIF_H
#define AMGIF_H

#include <memory>
#include <utility>
#include <vector>

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>

#include "constants.h"
#include "operators_buffer.h"

using std::pair;
using std::shared_ptr;
using std::tuple;
using std::vector;

// Smart pointer to matrix

// Shorthand for template types
using VecTriplet = tuple<gsl_vector *, gsl_vector *, gsl_vector *>;
using OpBuf = OperatorBuffer;

// Compute the AM-GIF estimates
VecTriplet computeAMGIF(
        OpBuf& Cbuf,
        gsl_vector *me,
        gsl_vector *pe,
        const gsl_matrix *L,
        const double alpha,
        const double beta,
        const double tau,
        const double eps
);

#endif // AMGIF_H

max_stars_repo_name: ichi-rika/glottal-inverse

max_issues_repo_name: ichi-rika/glottal-inverse

max_forks_repo_name: ichi-rika/glottal-inverse


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

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex64_t
 *
 *  CORE_zhegst - reduces a complex Hermitian-definite generalized
 *  eigenproblem to standard form.
 *  If PlasmaItype == 1, the problem is A*x = lambda*B*x, and A is
 *  overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
 *  If PlasmaItype == 2 or 3, the problem is A*B*x = lambda*x or B*A*x
 *  = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L.  B must
 *  have been previously factorized as U**H*U or L*L**H by
 *  CORE_zpotrf.
 *
 *******************************************************************************
 *
 * @param[in] itype
 *          Intended usage:
 *          = 1: A*x=(lambda)*B*x
 *          = 2: A*Bx=(lambda)*x
 *          = 3: B*A*x=(lambda)*x
 *
 * @param[in] uplo
 *          Specifies whether the matrix A is upper triangular or
 *          lower triangular:
 *          = PlasmaUpper: Upper triangle of A is stored;
 *          = PlasmaLower: Lower triangle of A is stored.
 *
 * @param[in] N
 *          The order of the matrices A and B. N >= 0.
 *
 * @param[in,out] A
 *          On entry, the symmetric (or Hermitian) matrix A.
 *          If uplo = PlasmaUpper, the leading N-by-N upper triangular
 *          part of A contains the upper triangular part of the matrix
 *          A, and the strictly lower triangular part of A is not
 *          referenced.
 *          If uplo = PlasmaLower, the leading N-by-N lower triangular
 *          part of A contains the lower triangular part of the matrix
 *          A, and the strictly upper triangular part of A is not
 *          referenced.
 *          On exit, if return value == 0, the transformed matrix,
 *          stored in the same format as A.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,N).
 *
 * @param[in,out] B
 *          On entry, the triangular factor from the Cholesky
 *          factorization of B, as returned by PLASMA_ZPOTRF.
 *
 * @param[in] LDB
 *          The leading dimension of the array B. LDB >= max(1,N).
 *
 * @param[out] INFO
 *          - 0 on successful exit
 *          - <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zhegst = PCORE_zhegst
#define CORE_zhegst PCORE_zhegst
#endif
void CORE_zhegst(int itype, PLASMA_enum uplo, int N,
                 PLASMA_Complex64_t *A, int LDA,
                 PLASMA_Complex64_t *B, int LDB, int *INFO)
{
    *INFO = LAPACKE_zhegst_work(
        LAPACK_COL_MAJOR,
        itype,
        lapack_const(uplo),
        N, A, LDA, B, LDB );
}


text:
/* ieee-utils/fp-solaris.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 <math.h>
#include <ieeefp.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>

int
gsl_ieee_set_mode (int precision, int rounding, int exception_mask)
{
  fp_except mode = 0 ;
  fp_rnd    rnd  = 0 ;

  switch (precision)
    {
    case GSL_IEEE_SINGLE_PRECISION:
      GSL_ERROR ("solaris only supports default precision rounding",
                 GSL_EUNSUP) ;
      break ;
    case GSL_IEEE_DOUBLE_PRECISION:
      GSL_ERROR ("solaris only supports default precision rounding",
                 GSL_EUNSUP) ;
      break ;
    case GSL_IEEE_EXTENDED_PRECISION:
      GSL_ERROR ("solaris only supports default precision rounding",
                 GSL_EUNSUP) ;
      break ;
    }

  switch (rounding)
    {
    case GSL_IEEE_ROUND_TO_NEAREST:
      rnd = FP_RN ;
      fpsetround (rnd) ;
      break ;
    case GSL_IEEE_ROUND_DOWN:
      rnd = FP_RM ;
      fpsetround (rnd) ;
      break ;
    case GSL_IEEE_ROUND_UP:
      rnd = FP_RP ;
      fpsetround (rnd) ;
      break ;
    case GSL_IEEE_ROUND_TO_ZERO:
      rnd = FP_RZ ;
      fpsetround (rnd) ;
      break ;
    default:
      rnd = FP_RN ;
      fpsetround (rnd) ;
    }

  /* Turn on all the exceptions apart from 'inexact' */

  mode = FP_X_INV | FP_X_DZ | FP_X_OFL | FP_X_UFL ;

  if (exception_mask & GSL_IEEE_MASK_INVALID)
    mode &= ~ FP_X_INV ;

  if (exception_mask & GSL_IEEE_MASK_DENORMALIZED)
    {
      /* do nothing */
    }
  else
    {
      GSL_ERROR ("solaris does not support the denormalized operand exception. "
                 "Use 'mask-denormalized' to work around this.",
                 GSL_EUNSUP) ;
    }

  if (exception_mask & GSL_IEEE_MASK_DIVISION_BY_ZERO)
    mode &= ~ FP_X_DZ ;

  if (exception_mask & GSL_IEEE_MASK_OVERFLOW)
    mode &= ~ FP_X_OFL ;

  if (exception_mask & GSL_IEEE_MASK_UNDERFLOW)
    mode &=  ~ FP_X_UFL ;

  if (exception_mask & GSL_IEEE_TRAP_INEXACT)
    {
      mode |= FP_X_IMP ;
    }
  else
    {
      mode &= ~ FP_X_IMP ;
    }

  fpsetmask (mode) ;

  return GSL_SUCCESS ;

}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* specfunc/gsl_sf_ellint.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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.
 */

/* Author: G. Jungman */

#ifndef __GSL_SF_ELLINT_H__
#define __GSL_SF_ELLINT_H__

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

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

__BEGIN_DECLS

/* Legendre form of complete elliptic integrals
 *
 * K(k) = Integral[1/Sqrt[1 - k^2 Sin[t]^2], {t, 0, Pi/2}]
 * E(k) = Integral[  Sqrt[1 - k^2 Sin[t]^2], {t, 0, Pi/2}]
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_ellint_Kcomp_e(double k, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_ellint_Kcomp(double k, gsl_mode_t mode);

int gsl_sf_ellint_Ecomp_e(double k, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_ellint_Ecomp(double k, gsl_mode_t mode);

int gsl_sf_ellint_Pcomp_e(double k, double n, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_ellint_Pcomp(double k, double n, gsl_mode_t mode);

int gsl_sf_ellint_Dcomp_e(double k, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_ellint_Dcomp(double k, gsl_mode_t mode);

/* Legendre form of incomplete elliptic integrals
 *
 * F(phi,k)   = Integral[1/Sqrt[1 - k^2 Sin[t]^2], {t, 0, phi}]
 * E(phi,k)   = Integral[  Sqrt[1 - k^2 Sin[t]^2], {t, 0, phi}]
 * P(phi,k,n) = Integral[(1 + n Sin[t]^2)^(-1)/Sqrt[1 - k^2 Sin[t]^2], {t, 0, phi}]
 * D(phi,k,n) = R_D(1-Sin[phi]^2, 1-k^2 Sin[phi]^2, 1.0)
 *
 * F: [Carlson, Numerische Mathematik 33 (1979) 1, (4.1)]
 * E: [Carlson, ", (4.2)]
 * P: [Carlson, ", (4.3)]
 * D: [Carlson, ", (4.4)]
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_ellint_F_e(double phi, double k, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_ellint_F(double phi, double k, gsl_mode_t mode);

int gsl_sf_ellint_E_e(double phi, double k, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_ellint_E(double phi, double k, gsl_mode_t mode);

int gsl_sf_ellint_P_e(double phi, double k, double n, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_ellint_P(double phi, double k, double n, gsl_mode_t mode);

int gsl_sf_ellint_D_e(double phi, double k, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_ellint_D(double phi, double k, gsl_mode_t mode);

/* Carlson's symmetric basis of functions
 *
 * RC(x,y)   = 1/2 Integral[(t+x)^(-1/2) (t+y)^(-1)], {t,0,Inf}]
 * RD(x,y,z) = 3/2 Integral[(t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-3/2), {t,0,Inf}]
 * RF(x,y,z) = 1/2 Integral[(t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2), {t,0,Inf}]
 * RJ(x,y,z,p) = 3/2 Integral[(t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2) (t+p)^(-1), {t,0,Inf}]
 *
 * exceptions: GSL_EDOM
 */
int gsl_sf_ellint_RC_e(double x, double y, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_ellint_RC(double x, double y, gsl_mode_t mode);

int gsl_sf_ellint_RD_e(double x, double y, double z, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_ellint_RD(double x, double y, double z, gsl_mode_t mode);

int gsl_sf_ellint_RF_e(double x, double y, double z, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_ellint_RF(double x, double y, double z, gsl_mode_t mode);

int gsl_sf_ellint_RJ_e(double x, double y, double z, double p, gsl_mode_t mode, gsl_sf_result * result);
double gsl_sf_ellint_RJ(double x, double y, double z, double p, gsl_mode_t mode);

__END_DECLS

#endif /* __GSL_SF_ELLINT_H__ */

max_stars_repo_name: shi-bash-cmd/qtTest


text:
/* examples/C/ssmfe/precond_core.f90 */
/* Laplacian on a square grid (using SPRAL_SSMFE_CORE routines) */
#include "spral.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <cblas.h>

/* Header that implements Laplacian and preconditioners */
#include "laplace2d.h"

int main(void) {
   const int ngrid = 20;         /* grid points along each side */
   const int n = ngrid*ngrid;    /* problem size */
   const int nep = 5;            /* eigenpairs wanted */
   const int m = 3;              /* dimension of the iterated subspace */
   const double tol = 1.e-6;     /* eigenvector tolerance */

   int state = SPRAL_RANDOM_INITIAL_SEED; /* PRNG state */

   int ind[m];                   /* permutation index */
   double lambda[n];             /* eigenvalues */
   double X[n][n];               /* eigenvectors */
   /* Work arrays */
   double lmd[m];
   double rr[3][2*m][2*m];
   double W[7][m][n];
   double U[m][n];

   /* Derived types */
   struct spral_ssmfe_rcid rci;              /* reverse communication data */
   struct spral_ssmfe_core_options options;  /* options */
   void *keep;                               /* private data */
   struct spral_ssmfe_inform inform;         /* information */

   /* Initialize options to default values */
   spral_ssmfe_core_default_options(&options);

   /* Initialize W to lin indep vectors by randomizing */
   for(int i=0; i<n; i++)
   for(int j=0; j<m; j++)
      W[0][j][i] = spral_random_real(&state, true);

   int ncon = 0;        /* number of converged eigenpairs */
   rci.job = 0; keep = NULL;
   while(true) { /* reverse communication loop */
      spral_ssmfe_double(&rci, 0, nep, 0, m, lmd, &rr[0][0][0], ind,
         &keep, &options, &inform);
      switch ( rci.job ) {
      case 1:
         apply_laplacian(
            ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0], &W[rci.ky][rci.jy][0]
            );
         break;
      case 2:
         apply_gauss_seidel_step (
            ngrid, ngrid, rci.nx, &W[rci.kx][rci.jx][0], &W[rci.ky][rci.jy][0]
            );
         break;
      case 4:
         for(int j=0; j<m; j++) {
            if ( inform.converged[j] != 0 ) continue;
            if ( inform.err_X[j] > 0 && inform.err_X[j] < tol )
               inform.converged[j] = 1;
         }
         break;
      case 5:
         if ( rci.i < 0 ) continue;
         for(int k=0; k<rci.nx; k++) {
           int j = ncon + k;
           lambda[j] = lmd[rci.jx + k];
           cblas_dcopy( n, &W[0][rci.jx+k][0], 1, &X[j][0], 1 );
         }
         ncon += rci.nx;
         if ( ncon >= nep || inform.iteration > 300 ) goto finished;
         break;
      case 11:
         if ( rci.i == 0 ) {
            if ( rci.kx != rci.ky || rci.jx > rci.jy ) {
               cblas_dcopy(n*rci.nx, &W[rci.kx][rci.jx][0], 1, &W[rci.ky][rci.jy][0], 1);
            } else if ( rci.jx < rci.jy ) {
               for(int j=rci.nx-1; j>=0; j--)
                  cblas_dcopy(n, &W[rci.kx][rci.jx+j][0], 1, &W[rci.ky][rci.jy+j][0], 1);
            }
         } else {
            for(int i=0; i<n; i++) {
               for(int j=0; j<rci.nx; j++)
                  U[j][i] = W[rci.kx][ind[j]][i];
               for(int j=0; j<rci.nx; j++)
                  W[rci.kx][j][i] = U[j][i];
               if(rci.ky != rci.kx) {
                  for(int j=0; j<rci.nx; j++)
                    U[j][i] = W[rci.ky][ind[j]][i];
                  for(int j=0; j<rci.nx; j++)
                    W[rci.ky][j][i] = U[j][i];
               }
            }
         }
         break;
      case 12:
         for(int i=0; i<rci.nx; i++)
           rr[rci.k][rci.j+i][rci.i+i] =
             cblas_ddot(n, &W[rci.kx][rci.jx+i][0], 1, &W[rci.ky][rci.jy+i][0], 1);
         break;
      case 13:
         for(int i=0; i<rci.nx; i++) {
            if( rci.kx == rci.ky ) {
               double s = cblas_dnrm2(n, &W[rci.kx][rci.jx+i][0], 1);
               if( s > 0 )
                  cblas_dscal(n, 1/s, &W[rci.kx][rci.jx+i][0], 1);
            } else {
               double s = sqrt(fabs(cblas_ddot(
                  n, &W[rci.kx][rci.jx+i][0], 1, &W[rci.ky][rci.jy+i][0], 1)
                  ));
               if ( s > 0 ) {
                  cblas_dscal(n, 1/s, &W[rci.kx][rci.jx+i][0], 1);
                  cblas_dscal(n, 1/s, &W[rci.ky][rci.jy+i][0], 1);
               } else {
                  for(int j=0; j<n; j++)
                     W[rci.ky][rci.jy+i][j] = 0.0;
               }
            }
         }
         break;
      case 14:
         for(int i=0; i<rci.nx; i++) {
           double s = -rr[rci.k][rci.j+i][rci.i+i];
           cblas_daxpy(n, s, &W[rci.kx][rci.jx+i][0], 1, &W[rci.ky][rci.jy+i][0], 1);
         }
         break;
      case 15:
         if ( rci.nx > 0 && rci.ny > 0 )
            cblas_dgemm(
               CblasColMajor, CblasTrans, CblasNoTrans, rci.nx, rci.ny, n,
               rci.alpha, &W[rci.kx][rci.jx][0], n, &W[rci.ky][rci.jy][0], n,
               rci.beta, &rr[rci.k][rci.j][rci.i], 2*m
               );
         break;
      case 16: // Fall through to 17
      case 17:
         if( rci.ny < 1 ) continue;
         if( rci.nx < 1 ) {
            if( rci.job == 17 ) continue;
            if( rci.beta == 1.0 ) continue;
            for(int j=rci.jy; j<rci.jy+rci.ny; j++)
               cblas_dscal(n, rci.beta, &W[rci.ky][j][0], 1);
            continue;
         }
         if( rci.job == 17 ) {
            cblas_dgemm(
               CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.ny, rci.nx,
               1.0, &W[rci.kx][rci.jx][0], n, &rr[rci.k][rci.j][rci.i], 2*m,
               0.0, &W[rci.ky][rci.jy][0], n
               );
            cblas_dcopy(n*rci.ny, &W[rci.ky][rci.jy][0], 1, &W[rci.kx][rci.jx][0], 1);
         } else {
            cblas_dgemm(
               CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.ny, rci.nx,
               rci.alpha, &W[rci.kx][rci.jx][0], n, &rr[rci.k][rci.j][rci.i],
               2*m, rci.beta, &W[rci.ky][rci.jy][0], n
               );
         }
         break;
      case 21: // Fall through to 22
      case 22:
         if( ncon > 0 ) {
            cblas_dgemm(
               CblasColMajor, CblasTrans, CblasNoTrans, ncon, rci.nx, n,
               1.0, &X[0][0], n, &W[rci.ky][rci.jy][0], n, 0.0, &U[0][0], n
               );
            cblas_dgemm(
               CblasColMajor, CblasNoTrans, CblasNoTrans, n, rci.nx, ncon,
               -1.0, &X[0][0], n, &U[0][0], n, 1.0, &W[rci.kx][rci.jx][0], n
               );
         }
         break;
      default:
         goto finished;
      }
   }
finished:
   if(inform.flag != 0) printf("inform.flag = %d\n", inform.flag);
   printf("%3d eigenpairs converged in %d iterations\n", ncon, inform.iteration);
   for(int i=0; i<ncon; i++)
      printf(" lambda[%1d] = %13.7e\n", i, lambda[i]);
   spral_ssmfe_core_free(&keep, &inform);

   /* Success */
   return 0;
}

max_stars_repo_name: mjacobse/spral

max_issues_repo_name: mjacobse/spral

max_forks_repo_name: mjacobse/spral


text:
/* Copyright (C) 2010-2019, The Regents of The University of Michigan.
 All rights reserved.

 This software was developed as part of the The Vulcan project in the Intelligent Robotics Lab
 under the direction of Benjamin Kuipers, kuipers@umich.edu. Use of this code is governed by an
 MIT-style License that can be found at "https://github.com/h2ssh/Vulcan".
*/

/**
* \file     local_trajectory_regression.h
* \author   Jong Jin Park and Collin Johnson
*
* Declaration of LocalTrajectoryRegression.
*/

#ifndef TRACKER_MOTIONS_LOCAL_TRAJECTORY_REGRESSION_H
#define TRACKER_MOTIONS_LOCAL_TRAJECTORY_REGRESSION_H

#include <tracker/object_state.h>
#include <core/point.h>
#include <boost/optional.hpp>
#include <nlopt.h>
#include <deque>
#include <vector>
#include <memory>

namespace vulcan
{
namespace tracker
{

struct regression_data_t
{
    int64_t timestamp;
    float   x;
    float   y;

    explicit regression_data_t(int64_t timestamp = 0, float x = 0.0f, float y = 0.0f)
    : timestamp(timestamp)
    , x(x)
    , y(y)
    {
    }
};

/**
* LocalTrajectoryRegression
*/
class LocalTrajectoryRegression
{
public:

    static const int kRegressionDimension = 4;
    static const int kIdxX  = 0;
    static const int kIdxY  = 1;
    static const int kIdxVX = 2;
    static const int kIdxVY = 3;

    // constructor
    LocalTrajectoryRegression(void);

    // Copy constructor
    LocalTrajectoryRegression(const LocalTrajectoryRegression& rhs);

    // destructor
    ~LocalTrajectoryRegression(void);

    // fit trajectory. Estimated velocity and trajectory are the outputs
    boost::optional<velocity_t> fitTrajectory(const std::deque<regression_data_t>& recentTrajectory,
                                              std::vector<Position>*               estimatedTrajectory = nullptr);

    // cost funtion for fitting the trajectory
    double evaluate(const double x[kRegressionDimension], double* grad);

private:

    void initializeData (const std::deque<regression_data_t>& recentTrajectory);

    // construct estimated trajectory using series of timestamps from data
    void predictTrajectory(double x0, double y0, double vx, double vy, std::vector<Position>& estimatedTrajectory);

    // computing energy and gradient
    double computeCost(void);
    void   computeGradient(double* grad);

    // setting up the optimizer
    void setupNLOpt(void);

    // optimizer
    nlopt_opt optimizer_;

    // data members
    std::vector<regression_data_t> data_;
    std::vector<Position>          prediction_;
    Position                       prior_;
    Point<double>            priorWeights_;

    int  numIter_;
    bool haveEnoughData_;

    // parameters
//     regression_params_t params_; // NOTE: hard coded for now
};

} // namespace mpepc
} // namespace vulcan

#endif // TRACKER_MOTIONS_LOCAL_TRAJECTORY_REGRESSION_H


text:
/*************************************************************
  Generalized-ICP Copyright (c) 2009 Aleksandr Segal.
  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.
* The names of the contributors may not be used to endorse
  or promote products derived from this software
  without specific prior written permission.

  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS 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 OWNER 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 OPTIMIZE_H_
#define OPTIMIZE_H_

#include <ANN.h>
#include "gicp.h"
#include <vector>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_multifit_nlin.h>

namespace dgc {
namespace gicp {
inline void print_gsl_matrix(gsl_matrix *mat, const char * name) {
    std::cout << name << "= [";
    for(unsigned int i = 0; i < mat->size1; i++) {
        for(unsigned int j = 0; j < mat->size2; j++) {
            std::cout << gsl_matrix_get(mat, i, j) << " ";
        }
        std::cout << ";" << std::endl;
    }
    std::cout << "]" << std::endl;
}

struct GICPOptData {
    GICPPointSet *p1;
    GICPPointSet *p2;
    ANNidx *nn_indecies; // nearest point indecies
    gicp_mat_t *M;      // mahalanobis matrices for each pair
    dgc_transform_t base_t;
    int num_matches;
    bool solve_rotation;
};

class GICPOptimizer {
public:
    GICPOptimizer();
    ~GICPOptimizer();

    int Iterations() { return iter; }
    const char* Status() { return gsl_strerror(status); }

    bool Optimize(dgc_transform_t t, GICPOptData &opt_data);
    bool OptimizeLM(dgc_transform_t t, GICPOptData &opt_data);

    void SetDebug(bool d) { debug = d; }

    void SetMaxIterations(int iter) { max_iter = iter; }

    void PlotError(dgc_transform_t t, GICPOptData &opt_data, const char* filename);

private:
    static double f(const gsl_vector * x, void * params);
    static void df(const gsl_vector * x, void * params, gsl_vector * g);
    static void fdf(const gsl_vector * x, void * params, double * f, gsl_vector * g);

    static void compute_dr(gsl_vector const* x, gsl_matrix const* gsl_temp_mat_r, gsl_vector *g);
    static double mat_inner_prod(gsl_matrix const* mat1, gsl_matrix const* mat2);
    static void apply_state(dgc_transform_t t, gsl_vector const* x);

    gsl_multimin_fdfminimizer *gsl_minimizer;
    gsl_vector *x;
    int max_iter;
    int iter;
    int status;
    bool debug;
    const static int N = 6;
    const gsl_multimin_fdfminimizer_type *T_min;
};

}
}

#endif

max_stars_repo_name: raandoom/gicp

max_issues_repo_name: raandoom/gicp

max_forks_repo_name: raandoom/gicp


text:
#include <gsl/gsl_interp.h>
#include "gsl_wrapper.h"

typedef struct tag_IntWrapper
{
    gsl_interp* pInterpolator;
    gsl_interp_accel *acc;

    double* pX;
    double* pY;
    int nPts;
} IntWrapper;

void* makeIntWrapper( double* x, double* y, int nPts)
{
    IntWrapper* pNewWrapper = (IntWrapper*) malloc( sizeof(IntWrapper));

    pNewWrapper->pInterpolator = gsl_interp_alloc (gsl_interp_linear, nPts );
    gsl_interp_init(pNewWrapper->pInterpolator, x,y, nPts );
    pNewWrapper->acc = gsl_interp_accel_alloc ();

    //Copy the data:
    pNewWrapper->pX = (double*) malloc( sizeof(double) * nPts);
    pNewWrapper->pY = (double*) malloc( sizeof(double) * nPts);
    int i;
    for(i=0;i< nPts;i++)
    {
        pNewWrapper->pX[i] = x[i];
        pNewWrapper->pY[i] = y[i];
    }
    pNewWrapper->nPts = nPts;

    return pNewWrapper;
}

double interpolate2( double x0, void * pInterpolator)
{
     //printf("%lf\n",x0);
    IntWrapper* pInt = (IntWrapper*) pInterpolator;

    if(x0 <= pInt->pX[0])
        return pInt->pY[0];
    if(x0 >= pInt->pX[pInt->nPts-1])
        return pInt->pY[pInt->nPts-1];
    return gsl_interp_eval ( pInt->pInterpolator, pInt->pX, pInt->pY, x0, pInt->acc);
}

int hello()
{
    return 0;
}

float MyInf()
{
    return 0.5;
}

max_stars_repo_name: mikehulluk/morphforge

max_issues_repo_name: mikehulluk/morphforge

max_forks_repo_name: mikehulluk/morphforge