pk_interp: allows for fewer z' (down to 1)

[JesusTorrado]

Hi Antony,

This patch allows to define the P(k) interpolator with fewer z's (it failed before for len(z)<4 because of the cubic degree of the spline):

• If len(z)<4: reduces the degree of the 2DInterpolator
• If len(z)=1: creates an scipy.interp1d object, that returns equivalent arrays in terms of dimensionality, so the new behaviour is transparent the user (except it produces an error if one requests a z different from the single one computed)

Can be tested with the following code (notice the difference in the 2-z's case (intermediate z is linearly interpolated):

import sys, platform, os
import matplotlib
from matplotlib import pyplot as plt
import numpy as np

# Path to the patched version!!!
sys.path.insert(0,"/home/jesus/projects/cobaya/CAMB_pk/CAMB")
import camb
from camb import model

params = {'As': 2.204090400490399e-09, 'ns': 0.9647522, 'cosmomc_theta': 0.01040778,
         'ombh2': 0.02225203, 'omch2': 0.1198657, 'mnu': 0.06, 'tau': 0.07888604,
         'num_massive_neutrinos': 1, 'nnu': 3.046, 'theta_H0_range': [40, 100],
         'halofit_version': 'mead', 'H0': None, 'kmax': 15}
zs_1 = [3]
zs_2 = [3,1]
zs_4 = [3,2,1]

from copy import deepcopy

pars_1 = deepcopy(params)
pars_1.update({"redshifts": zs_1})
cambpars_1 = camb.set_params(**pars_1)
results_1 = camb.get_results(cambpars_1)
PK_1 = results_1.get_matter_power_interpolator(nonlinear=True,
    hubble_units=False, k_hunit=False, var1="delta_tot", var2="delta_tot")

pars_2 = deepcopy(params)
pars_2.update({"redshifts": zs_2})
cambpars_2 = camb.set_params(**pars_2)
results_2 = camb.get_results(cambpars_2)
PK_2 = results_2.get_matter_power_interpolator(nonlinear=True,
    hubble_units=False, k_hunit=False, var1="delta_tot", var2="delta_tot")

pars_4 = deepcopy(params)
pars_4.update({"redshifts": zs_4})
cambpars_4 = camb.set_params(**pars_4)
results_4 = camb.get_results(cambpars_4)
PK_4 = results_4.get_matter_power_interpolator(nonlinear=True,
    hubble_units=False, k_hunit=False, var1="delta_tot", var2="delta_tot")

# Plots!
plt.figure(figsize=(8,5))
k=np.exp(np.log(10)*np.linspace(-4,1,200))
zplot = zs_4
pks = [PK_1, PK_2, PK_4]
names = ["one", "two", "four"]
styles = [":", "--", "-"]
colors = ["r", "g", "b"]
sizes = ["3", "2", "1"]
for z in zplot:
    for i, pk in enumerate(pks):
        try:  # catches the single-z case
            plt.loglog(k, pk.P(z,k), c=colors[i], ls=styles[i],
                       lw=sizes[i], label="%s, z=%g"%(names[i], z))
        except:
            pass
plt.xlim([1e-4,params["kmax"]])
plt.xlabel('k Mpc')
plt.ylabel('$P_\Psi\, Mpc^{-3}$')
plt.legend();

# Error if wrong z requested in single-z case
PK_1.P(1,k)