ccl_emu17.h File Reference

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Macros
#define	A_MIN_EMU   1./3.
#define	K_MAX_EMU   5.0
#define	K_MIN_EMU   1.0000000474974513E-003
#define	NK_EMU   351

Functions
CCL_BEGIN_DECLS void	ccl_pkemu (double *xstarin, double **Pkemu, int *status, ccl_cosmology *cosmo)

Macro Definition Documentation

#define A_MIN_EMU   1./3.

Definition at line 4 of file ccl_emu17.h.

Referenced by ccl_cosmology_compute_power_emu(), ccl_linear_matter_power(), and ccl_nonlin_matter_power().

#define K_MAX_EMU   5.0

Definition at line 5 of file ccl_emu17.h.

Referenced by ccl_cosmology_compute_power_emu().

#define K_MIN_EMU   1.0000000474974513E-003

Definition at line 6 of file ccl_emu17.h.

Referenced by ccl_cosmology_compute_power_emu().

#define NK_EMU   351

Definition at line 7 of file ccl_emu17.h.

Referenced by ccl_cosmology_compute_power_emu(), and ccl_pkemu().

Function Documentation

CCL_BEGIN_DECLS void ccl_pkemu	(	double *	xstarin,
		double **	Pkemu,
		int *	status,
		ccl_cosmology *	cosmo
	)

Emulator power spectrum Obtain P(k,z) [Mpc^3] for a given set of input parameters.

Parameters

xstarin	vector of input parameters for the emulator, including redshift.
Pkemu	output P(k,z) power spectrum
status	Status flag. 0 if there are no errors, nonzero otherwise.
cosmo	Cosmology parameters and configurations (only relevant for storing status)

Definition at line 141 of file ccl_emu17.c.

References beta, ccl_cosmology_set_status_message(), CCL_ERROR_EMULATOR_BOUND, ccl_raise_exception(), emuInit(), K, KrigBasis, lamz, m, M_PI, mean, mode, neta, NK_EMU, p, peta, pow(), rs, sd, ccl_cosmology::status_message, x, xmax, xmin, xrange, and z.

Referenced by ccl_cosmology_compute_power_emu().

                                                                                 {
    
    static double inited=0;
    int ee, i, j, k;
    double wstar[peta[0]+peta[1]];
    double Sigmastar[2][peta[1]][m[0]];
    double ystaremu[neta];
    *ystar=(double *)malloc(sizeof(double)*NK_EMU);
    double ybyz[rs];
    double logc;
    double xstarstd[p];
    int zmatch;
    gsl_spline *zinterp = gsl_spline_alloc(gsl_interp_cspline, rs);
    gsl_interp_accel *accel = gsl_interp_accel_alloc();
    
    
    // Initialize if necessary
    if(inited==0) {
        emuInit();
        inited = 1;
    }
    
    // Transform w_a into (-w_0-w_a)^(1/4)
    xstar[6] = pow(-xstar[5]-xstar[6], 0.25);
    // Check the inputs to make sure we're interpolating.
    for(i=0; i<p; i++) {
        if((xstar[i] < xmin[i]) || (xstar[i] > xmax[i])) {
            switch(i) {
                case 0:
                    ccl_cosmology_set_status_message(cosmo, 
                            "ccl_pkemu(): omega_m must be between %f and %f.\n", 
                            xmin[i], xmax[i]);
                    break;
                case 1:
                    ccl_cosmology_set_status_message(cosmo, 
                            "ccl_pkemu(): omega_b must be between %f and %f.\n", 
                            xmin[i], xmax[i]);
                    break;
                case 2:
                    ccl_cosmology_set_status_message(cosmo, 
                            "ccl_pkemu(): sigma8 must be between %f and %f.\n", 
                            xmin[i], xmax[i]);
                    break;
                case 3:
                    ccl_cosmology_set_status_message(cosmo, 
                            "ccl_pkemu(): h must be between %f and %f.\n", 
                            xmin[i], xmax[i]);
                    break;
                case 4:
                    ccl_cosmology_set_status_message(cosmo, 
                            "ccl_pkemu(): n_s must be between %f and %f.\n", 
                            xmin[i], xmax[i]);
                    break;
                case 5:
                    ccl_cosmology_set_status_message(cosmo, 
                            "ccl_pkemu(): w_0 must be between %f and %f.\n", 
                            xmin[i], xmax[i]);
                    break;
                case 6:
                    ccl_cosmology_set_status_message(cosmo, 
                            "ccl_pkemu(): (-w_0-w_a)^(1/4) must be between %f and %f.\n", 
                            xmin[i], xmax[i]);
                    break;
                case 7:
                    ccl_cosmology_set_status_message(cosmo, 
                            "ccl_pkemu(): omega_nu must be between %f and %f.\n", 
                            xmin[i], xmax[i]);
                    break;
            }
            *status = CCL_ERROR_EMULATOR_BOUND;
            ccl_raise_exception(*status, cosmo->status_message);
        gsl_spline_free(zinterp);
        gsl_interp_accel_free(accel);
            return;
        }
    } // for(i=0; i<p; i++)
    if((xstar[p] < z[0]) || (xstar[p] > z[rs-1])) {
        ccl_cosmology_set_status_message(cosmo, 
                "ccl_pkemu(): z must be between %f and %f.\n", 
                z[0], z[rs-1]);
        *status = CCL_ERROR_EMULATOR_BOUND;
        ccl_raise_exception(*status, cosmo->status_message);
    gsl_spline_free(zinterp);
    gsl_interp_accel_free(accel);
        return;
    }
    
    // Standardize the inputs
    for(i=0; i<p; i++) {
        xstarstd[i] = (xstar[i] - xmin[i]) / xrange[i];
    }
    
    // compute the covariances between the new input and sims for all the PCs.
    for(ee=0; ee<2; ee++) {
        for(i=0; i<peta[ee]; i++) {
            for(j=0; j<m[ee]; j++) {
                logc = 0.0;
                for(k=0; k<p; k++) {
                    logc -= beta[ee][i][k]*pow(x[j][k] - xstarstd[k], 2.0);
                }
                Sigmastar[ee][i][j] = exp(logc) / lamz[ee][i];
            }
        }
    }
    
    // Compute wstar for the first chunk.
    for(i=0; i<peta[0]; i++) {
        wstar[i]=0.0;
        for(j=0; j<m[0]; j++) {
            wstar[i] += Sigmastar[0][i][j] * KrigBasis[0][i][j];
        }
    }
    // Compute wstar for the second chunk.
    for(i=0; i<peta[1]; i++) {
        wstar[peta[0]+i]=0.0;
        for(j=0; j<m[1]; j++) {
            wstar[peta[0]+i] += Sigmastar[1][i][j] * KrigBasis[1][i][j];
        }
    }
    
    /*
    for(i=0; i<peta[0]+peta[1]; i++) {
        printf("%f\n", wstar[i]);
    }
     */
    
    // Compute ystar, the new output
    for(i=0; i<neta; i++) {
        ystaremu[i] = 0.0;
        for(j=0; j<peta[0]+peta[1]; j++) {
            ystaremu[i] += K[i][j]*wstar[j];
        }
        ystaremu[i] = ystaremu[i]*sd + mean[i];
                
        // Convert to P(k)
        //ystaremu[i] = ystaremu[i] - 1.5*log10(mode[i % NK_EMU]);
        //ystaremu[i] = 2*M_PI*M_PI*pow(10, ystaremu[i]);
    }
    
    
    
    // Interpolate to the desired redshift
    // Natural cubic spline interpolation over z.
    
    // First check to see if the requested z is one of the training z.
    zmatch = -1;
    for(i=0; i<rs; i++) {
        if(xstar[p] == z[i]) {
            zmatch = rs-i-1;
        }
    }
    
    // z doesn't match a training z, interpolate
    if(zmatch == -1) {
        for(i=0; i<NK_EMU; i++) {
            for(j=0; j<rs; j++) {
                ybyz[rs-j-1] = ystaremu[j*NK_EMU+i];
            }
            gsl_spline_init(zinterp, z, ybyz, rs);
            (*ystar)[i] = gsl_spline_eval(zinterp, xstar[p], accel);
            gsl_interp_accel_reset(accel);
        }
        
    } else { //otherwise, copy in the emulated z without interpolating
        for(i=0; i<NK_EMU; i++) {
      (*ystar)[i] = ystaremu[zmatch*NK_EMU + i];
        }
    }
    
    gsl_spline_free(zinterp);
    gsl_interp_accel_free(accel);
    
    // Convert to P(k)
    for(i=0; i<NK_EMU; i++) {
      (*ystar)[i] = (*ystar)[i] - 1.5*log10(mode[i]) + log10(2) + 2*log10(M_PI);
      (*ystar)[i] = pow(10, (*ystar)[i]);
    }
}