Extrap_Pk< T, N > Class Template Reference

: linear interpolation of data [k, P(k)] within 'useful' range fit to primordial P_i(k) below the 'useful' range fit to Padé approximant R [0/3] above the 'useful' range More...

#include <core_power.h>

Inheritance diagram for Extrap_Pk< T, N >:

Collaboration diagram for Extrap_Pk< T, N >:

Public Member Functions
	Extrap_Pk (const Data_Vec< T, N > &data, const Sim_Param &sim)
	Extrap_Pk (const Data_Vec< T, N > &data, const Sim_Param &sim, const size_t m_l, const size_t n_u)
	Extrap_Pk (const Data_Vec< T, N > &data, const Sim_Param &sim, const size_t m_l, const size_t n_l, const size_t m_u, const size_t n_u)
double	operator() (double k) const
void	fit_lin (const Data_Vec< T, N > &data, const size_t m, const size_t n, double &A)
void	fit_power_law (const Data_Vec< T, N > &data, const size_t m, const size_t n, double &A, double &n_s)
Public Member Functions inherited from Interp_obj
	Interp_obj ()
	~Interp_obj ()
double	operator() (double x) const
template<typename T , size_t N>
void	init (const Data_Vec< T, N > &data)

Public Attributes
double	A_low
	amplitude of linear power in lower range More...
const Cosmo_Param &	cosmo
double	A_up
double	n_s
	scale-free power spectrum in upper range More...
T	k_min
T	k_max
	interpolation range More...
Public Attributes inherited from Interp_obj
double	x_min
double	x_max

Detailed Description

template<typename T, size_t N>
class Extrap_Pk< T, N >

: linear interpolation of data [k, P(k)] within 'useful' range fit to primordial P_i(k) below the 'useful' range fit to Padé approximant R [0/3] above the 'useful' range

Steffen interpolation of data [k, P(k)] within range k_min, k_max fit to primordial P_i(k) below this range, fit A*k^ns above

Definition at line 122 of file core_power.h.

Constructor & Destructor Documentation

template<typename T , size_t N>

Extrap_Pk< T, N >::Extrap_Pk	(	const Data_Vec< T, N > &	data,
		const Sim_Param &	sim
	)

Definition at line 548 of file core_power.cpp.

                                                                          :
    Extrap_Pk(data, sim,
        // trust simulation up to HALF k_nq
        0, get_nearest(sim.other_par.nyquist.at("particle")/2, data[0]) + 1
    ) {}

template<typename T , size_t N>

Extrap_Pk< T, N >::Extrap_Pk	(	const Data_Vec< T, N > &	data,
		const Sim_Param &	sim,
		const size_t	m_l,
		const size_t	n_u
	)

Definition at line 540 of file core_power.cpp.

                                                                                                              :
    Extrap_Pk(data, sim,
        // fit over first and last half of a decade
        m_l, m_l + sim.out_opt.bins_per_decade/2,
        n_u - sim.out_opt.bins_per_decade/2, n_u
    ) {}

template<typename T , size_t N>

Extrap_Pk< T, N >::Extrap_Pk	(	const Data_Vec< T, N > &	data,
		const Sim_Param &	sim,
		const size_t	m_l,
		const size_t	n_l,
		const size_t	m_u,
		const size_t	n_u
	)

Definition at line 523 of file core_power.cpp.

References Extrap_Pk< T, N >::A_low, Extrap_Pk< T, N >::A_up, Extrap_Pk< T, N >::fit_lin(), Extrap_Pk< T, N >::fit_power_law(), Interp_obj::init(), Extrap_Pk< T, N >::k_max, Extrap_Pk< T, N >::k_min, and Extrap_Pk< T, N >::n_s.

                                       :
    cosmo(sim.cosmo)
{
    this->init(data);//< initialize Interp_obj
    // LOWER RANGE -- fit linear power spectrum to data[m_l:n_l)
    BOOST_LOG_TRIVIAL(debug) << "Fitting amplitude of P(k) in lower range.";
    k_min = data[0][m_l]; // first k in data
    fit_lin(data, m_l, n_l, A_low);

    // UPPER RANGE -- fit Ak^ns to data[m_u,n_u)
    BOOST_LOG_TRIVIAL(debug) << "Fitting amplitude of P(k) in upper range.";
    k_max =  data[0][n_u-1]; // last k in data
    fit_power_law(data, m_u, n_u, A_up, n_s);
}

Member Function Documentation

template<typename T , size_t N>

void Extrap_Pk< T, N >::fit_lin	(	const Data_Vec< T, N > &	data,
		const size_t	m,
		const size_t	n,
		double &	A
	)

Definition at line 560 of file core_power.cpp.

References Extrap_Pk< T, N >::cosmo, Cosmo_Param::k2_G, lin_pow_spec(), m, pow2(), sqrt(), Cosmo_Param::truncated_pk, truncation_fce(), and w.

Referenced by Extrap_Pk< T, N >::Extrap_Pk().

{
    // FIT linear power spectrum to data[m:n)
    // fit 'P(k) = A * P_lin(k)' via A
    // for TZA correct for truncation
    // for N = 2 perform non-weighted least-square fitting
    // for N = 3 use data[2] as sigma, w = 1/sigma^2
    double Pk;
    std::vector<double> Pk_res, w;
    Pk_res.reserve(n-m);
    if (N == 3){
        w.reserve(n-m);
    }
    std::vector<double> A_vec(n-m, 1);

    if (cosmo.truncated_pk)
    {
        const FTYPE_t k2_G = cosmo.k2_G;
        for(size_t i = m; i < n; i++){
            FTYPE_t k = data[0][i];
            Pk = lin_pow_spec(1, k, cosmo) * truncation_fce(k, k2_G);;
            Pk_res.push_back(data[1][i] / Pk);
            if (N == 3) w.push_back(pow2(Pk/data[2][i]));
        }
    }
    else
    {
        for(size_t i = m; i < n; i++){
            Pk = lin_pow_spec(1, data[0][i], cosmo);
            Pk_res.push_back(data[1][i] / Pk);
            if (N == 3) w.push_back(pow2(Pk/data[2][i]));
        }
    }
    
    
    double A_sigma2, sumsq;

    int gsl_errno;
    if (N == 3) gsl_errno = gsl_fit_wmul(A_vec.data(), 1, w.data(), 1, Pk_res.data(), 1, n-m, &A, &A_sigma2, &sumsq);
    else gsl_errno = gsl_fit_mul(A_vec.data(), 1, Pk_res.data(), 1, n-m, &A, &A_sigma2, &sumsq);
    if (gsl_errno) throw std::runtime_error("GSL integration error: " + std::string(gsl_strerror(gsl_errno)));

    BOOST_LOG_TRIVIAL(debug) << "\t[" << (N == 3 ? "weighted-" : "") << "fit A = " << A << ", err = " << 100*sqrt(A_sigma2)/A << "%%]";
}

template<typename T , size_t N>

void Extrap_Pk< T, N >::fit_power_law	(	const Data_Vec< T, N > &	data,
		const size_t	m,
		const size_t	n,
		double &	A,
		double &	n_s
	)

Definition at line 606 of file core_power.cpp.

References m, Extrap_Pk< T, N >::n_s, pow2(), sqrt(), and w.

Referenced by Extrap_Pk< T, N >::Extrap_Pk().

{
    // FIT scale-free power spectrum to data[m,n)
    // fit 'log P(k) = log A + n_s * log k' via A, n_s
    // for N = 2 perform non-weighted least-square fitting
    // for N = 3 use data[2] as sigma, w = 1/sigma^2
    std::vector<double> k, Pk, w; //< need double for GSL
    k.reserve(n-m);
    Pk.reserve(n-m);
    if (N == 3){
        w.reserve(n-m);
    }
    for(size_t i = m; i < n; i++){
        k.push_back(log(data[0][i]));
        Pk.push_back(log(data[1][i]));
        if (N == 3)w.push_back(pow2(data[1][i]/data[2][i])); // weight = 1/sigma^2, approx log(1+x) = x for x rel. error
    }
    double cov00, cov01, cov11, sumsq;
    int gsl_errno;
    if (N == 3) gsl_errno = gsl_fit_wlinear(k.data(), 1, w.data(), 1, Pk.data(), 1, n-m, &A, &n_s, &cov00, &cov01, &cov11, &sumsq);
    else gsl_errno = gsl_fit_linear(k.data(), 1, Pk.data(), 1, n-m, &A, &n_s, &cov00, &cov01, &cov11, &sumsq);
    if (gsl_errno) throw std::runtime_error("GSL integration error: " + std::string(gsl_strerror(gsl_errno)));

    A = exp(A); // log A => A
    BOOST_LOG_TRIVIAL(debug) << "\t[" << (N == 3 ? "weighted-" : "") << "fit A = " << A << ", err = " << 100*sqrt(cov00) << "%%"
                             << "n_s = " << n_s << ", err = " << 100*sqrt(cov11)/std::abs(n_s) << "%%, corr = " << 100*cov01/sqrt(cov00*cov11) << "%%]";
}

template<typename T , size_t N>

double Extrap_Pk< T, N >::operator()

(

double

k

)

const

Definition at line 635 of file core_power.cpp.

References Extrap_Pk< T, N >::A_low, Extrap_Pk< T, N >::A_up, Extrap_Pk< T, N >::cosmo, Extrap_Pk< T, N >::k_max, Extrap_Pk< T, N >::k_min, lin_pow_spec(), Extrap_Pk< T, N >::n_s, Interp_obj::operator()(), and pow().

Referenced by Extrap_Pk_Nl< T, N >::operator()().

{
    if (k < k_min)
    {
        return A_low*lin_pow_spec(1, k, cosmo);
    }
    else if (k <= k_max) return Interp_obj::operator()(k);
    else return A_up*pow(k, n_s);
}

Member Data Documentation

template<typename T, size_t N>

double Extrap_Pk< T, N >::A_low

amplitude of linear power in lower range

Definition at line 134 of file core_power.h.

Referenced by Extrap_Pk< T, N >::Extrap_Pk(), and Extrap_Pk< T, N >::operator()().

template<typename T, size_t N>

double Extrap_Pk< T, N >::A_up

Definition at line 136 of file core_power.h.

Referenced by Extrap_Pk< T, N >::Extrap_Pk(), and Extrap_Pk< T, N >::operator()().

template<typename T, size_t N>

const Cosmo_Param& Extrap_Pk< T, N >::cosmo

Definition at line 135 of file core_power.h.

Referenced by Extrap_Pk< T, N >::fit_lin(), Extrap_Pk< T, N >::operator()(), and Extrap_Pk_Nl< T, N >::operator()().

template<typename T, size_t N>

T Extrap_Pk< T, N >::k_max

interpolation range

Definition at line 137 of file core_power.h.

Referenced by Extrap_Pk< T, N >::Extrap_Pk(), and Extrap_Pk< T, N >::operator()().

template<typename T, size_t N>

T Extrap_Pk< T, N >::k_min

Definition at line 137 of file core_power.h.

Referenced by Extrap_Pk< T, N >::Extrap_Pk(), and Extrap_Pk< T, N >::operator()().

template<typename T, size_t N>

double Extrap_Pk< T, N >::n_s

scale-free power spectrum in upper range

Definition at line 136 of file core_power.h.

Referenced by Extrap_Pk< T, N >::Extrap_Pk(), Extrap_Pk< T, N >::fit_power_law(), and Extrap_Pk< T, N >::operator()().

---

The documentation for this class was generated from the following files:

• src/core/include/core_power.h
• src/core/core_power.cpp