Main new features are P_cb, P_k_eq, RealSpaceInterface, new notebooks…

… and few bugfixes

* devel: (54 commits)
  updated doc
  added a few comments in spectra related to P_cb
  updated version  number to 2.7.0
  fixed minor issue with dcdm in input
  restored the empty directory RealSpaceInterface/cache/
  added 11 notebooks and equivalent scripts (cambridge version)
  pushing the Real Space Interface
  further cosmetic changes in P_cb
  few extra comkments to Pcb part
  Fixed a bug leading to double free when combining the new quadrature strategies with shooting
  Fixed typo in ncdm number density, cf. issue #225
  final polishing of pk_eq method
  taken care of the case thatan x-array with 2 entries only needs to be splined by using the natural_spline and not the estimated derivative spline method by default in this particular case
  sigma8 and sigma8_cb now computed in two separate functions to avoid confusion
  mem leak fix
  k_nl
  Pcb changes: M and CB share the same tau_min_nl
  Nils pk_type trick
  transfer.c dNdz according to Euclid IST specification
  removed the enum halofit_found_k_max
  ...

RealSpaceInterface/Calc2D/CalculationClass.py

@@ -0,0 +1,187 @@
+import os
+import logging
+
+import cv2
+import numpy as np
+
+from classy import Class
+
+from Calc2D.TransferFunction import ComputeTransferFunctionList
+from Calc2D.DataGeneration import GenerateGaussianData, GenerateSIData
+from Calc2D.DataPropagation import PropagateDatawithList
+from Calc2D.rFourier import *
+from Calc2D.Database import Database
+from collections import namedtuple
+
+import config
+
+ClSpectrum = namedtuple("Cls", ["l", "tCl"])
+PkSpectrum = namedtuple("Pkh", ["kh", "Pkh"])
+
+def normalize(real):
+    """
+    Given the `real` data, i.e. either a 2d array or a flattened 1d array
+    of the relative density perturbations, normalize its values as follows:
+
+    The client expects the values to be in the interval of [-1, 1].
+    Take a symmetric interval around a `real` value of 0 and linearly
+    map it to the required interval [-1, 1].
+    """
+    minimum, maximum = real.min(), real.max()
+    bound = max(abs(maximum), abs(minimum))
+    result = real / bound
+    return result
+
+class Calculation(object):
+    def __init__(self,
+                 kbins,
+                 resolution = 200,
+                 gauge="newtonian",
+                 kperdecade=200,
+                 P_k_max=100,
+                 evolver=1):
+        # also sets `endshape` through setter
+        self.resolution = resolution
+        self.gauge = gauge
+        self.evolver = evolver
+        self.P_k_max = P_k_max
+        self.kperdecade = kperdecade
+
+        self.redshift = None # to be set later
+        self.size = None # to be set later
+
+        self.krange = np.logspace(-4, 1, kbins)
+
+    @property
+    def resolution(self):
+        return self._resolution
+
+    @resolution.setter
+    def resolution(self, resolution):
+        self._resolution = resolution
+        self.endshape = (resolution, resolution)
+
+    def getData(self, redshiftindex):
+        FValuenew = PropagateDatawithList(
+            k=self.k,
+            FValue=self.FValue,
+            zredindex=redshiftindex,
+            transferFunctionlist=self.TransferFunctionList)
+
+        Valuenew = dict()
+        FValue_abs = np.abs(self.FValue)
+        _min, _max = FValue_abs.min(), FValue_abs.max()
+        dimensions = (self.endshape[0] / 2, self.endshape[1])
+        for quantity, FT in FValuenew.items():
+            FT_abs = np.abs(FT)
+            FT_normalized = cv2.resize(FT_abs, dimensions).ravel()
+            FT_normalized = (FT_normalized - _min) / (_max - _min)
+            real = realInverseFourier(FT.reshape(self.FValue.shape))
+            # real = cv2.resize(real, self.endshape).ravel()
+            real = real.ravel()
+
+            minimum, maximum = real.min(), real.max()
+            Valuenew[quantity] = normalize(real)
+
+        return Valuenew, FValuenew, (minimum, maximum)
+
+
+    def getInitialData(self):
+        # for odd values of self._resolution, this is necessary
+        Value = cv2.resize(realInverseFourier(self.FValue), self.endshape)
+
+        minimum, maximum = Value.min(), Value.max()
+        Value = normalize(Value)
+
+        assert Value.size == self.resolution ** 2
+
+        return Value.ravel(), cv2.resize(
+            (np.abs(self.FValue) - np.abs(self.FValue).min()) /
+            (np.abs(self.FValue).max() - np.abs(self.FValue).min()),
+            (self.endshape[0] / 2, self.endshape[1])).ravel(), (minimum,
+                                                                maximum)
+
+    def getTransferData(self, redshiftindex):
+        return {field: transfer_function[redshiftindex](self.krange) for field, transfer_function in self.TransferFunctionList.items()}, self.krange
+
+    def setCosmologialParameters(self, cosmologicalParameters):
+        self.cosmologicalParameters = cosmologicalParameters
+
+        # Calculate transfer functions
+        self.TransferFunctionList = ComputeTransferFunctionList(self.cosmologicalParameters, self.redshift)
+        # Calculate Cl's
+        self.tCl, self.mPk = self.calculate_spectra(self.cosmologicalParameters)
+
+    def calculate_spectra(self, cosmo_params, force_recalc=False):
+        settings = cosmo_params.copy()
+        settings.update({
+            "output": "tCl,mPk",
+            "evolver": "1",
+            "gauge": "newtonian",
+            "P_k_max_1/Mpc": 10,
+            })
+
+        database = Database(config.DATABASE_DIR, "spectra.dat")
+
+        if settings in database and not force_recalc:
+            data = database[settings]
+            ell = data["ell"]
+            tt = data["tt"]
+            kh = data["kh"]
+            Pkh = data["Pkh"]
+            self.z_rec = data["z_rec"]
+        else:
+            cosmo = Class()
+            cosmo.set(settings)
+            cosmo.compute()
+            # Cl's
+            data = cosmo.raw_cl()
+            ell = data["ell"]
+            tt = data["tt"]
+            # Matter spectrum
+            k = np.logspace(-3, 1, config.MATTER_SPECTRUM_CLIENT_SAMPLES_PER_DECADE * 4)
+            Pk = np.vectorize(cosmo.pk)(k, 0)
+            kh = k * cosmo.h()
+            Pkh = Pk / cosmo.h()**3
+            # Get redshift of decoupling
+            z_rec = cosmo.get_current_derived_parameters(['z_rec'])['z_rec']
+            self.z_rec = z_rec
+            # Store to database
+            database[settings] = {
+            "ell": data["ell"],
+            "tt": data["tt"],
+
+            "kh": k,
+            "Pkh": Pk,
+
+            "z_rec": z_rec,
+            }
+
+        return ClSpectrum(ell[2:], tt[2:]), PkSpectrum(kh, Pkh)
+
+    @property
+    def z_dec(self):
+        if self.z_rec is None:
+            raise ValueError("z_rec hasn't been computed yet")
+        return self.z_rec
+
+    def setInitialConditions(self,
+                             A=1,
+                             sigma=2,
+                             initialDataType="SI",
+                             SIlimit=None,
+                             SI_ns=0.96):
+        logging.info("Generating Initial Condition")
+
+        if initialDataType == "Gaussian":
+            self.ValueE, self.FValue, self.k, self.kxE, self.kyE = GenerateGaussianData(
+                sigma, self.size, self.resolution)
+        elif initialDataType == "SI":
+            self.ValueE, self.FValue, self.k, self.kxE, self.kyE = GenerateSIData(
+                A,
+                self.size,
+                self.resolution,
+                limit=SIlimit,
+                ns=SI_ns)
+        else:
+            logging.warn("initialDataType " + str(initialDataType) + " not found")

RealSpaceInterface/Calc2D/DataGeneration.py

@@ -0,0 +1,91 @@
+import logging
+
+import numpy as np
+import cv2
+
+from Calc2D.rFourier import realFourier, realInverseFourier
+
+def GenerateGaussianData(sigma, size, points, A=1):
+    xr = np.linspace(-size / 2.0, size / 2.0, points)
+    yr = np.linspace(-size / 2.0, size / 2.0, points)
+    step = xr[1] - xr[0]
+    x, y = np.meshgrid(
+        xr, yr, indexing='ij', sparse=True)  # indexing is important
+    del xr, yr
+
+    #use the more easy formula
+    Value = A * np.exp(-(x**2 + y**2) / (2 * sigma**2))
+
+    kx, ky, FValue = realFourier(step, Value)
+    kxr, kyr = np.meshgrid(kx, ky, indexing='ij', sparse=True)
+
+    k = np.sqrt(kxr**2 + kyr**2)
+    del kxr, kyr
+
+    kx = (min(kx), max(kx))  #just return the extremal values to save memory
+    ky = (min(ky), max(ky))
+
+    ValueE = (Value.min(), Value.max())
+
+    return ValueE, FValue, k, kx, ky
+
+def GenerateSIData(A, size, points, limit=None, ns=0.96):
+    xr = np.linspace(-size / 2.0, size / 2.0, points)
+    yr = np.linspace(-size / 2.0, size / 2.0, points)
+    step = xr[1] - xr[0]
+
+    x, y = np.meshgrid(
+        xr, yr, indexing='ij', sparse=True)  # indexing is important
+    del xr, yr
+    Value = 0 * x + 0 * y
+
+    kx, ky, FValue = realFourier(step, Value)  #FValue==0
+
+    kxr, kyr = np.meshgrid(kx, ky, indexing='ij', sparse=True)
+
+    k = np.sqrt(kxr**2 + kyr**2)
+    del kxr, kyr
+
+    if limit == None:
+
+        ktilde = k.flatten()
+        ktilde[np.argmin(k)] = 10**9  #just let the background be arbitrary low
+        ktilde = ktilde.reshape(k.shape)
+
+        FValue = np.random.normal(
+            loc=0,
+            scale=np.sqrt(A / ktilde**(
+                2 - (ns - 1) * 2. / 3.)) / np.sqrt(2)) + np.random.normal(
+                    loc=0,
+                    scale=np.sqrt(A / ktilde**
+                                  (2 - (ns - 1) * 2. / 3.)) / np.sqrt(2)) * 1j
+
+    elif type(limit) == list or type(limit) == tuple:
+
+        iunder, junder = np.where(k < limit[1])
+
+        for t in range(len(iunder)):
+
+            if k[iunder[t]][junder[t]] > limit[0] and k[iunder[t]][junder[t]] > 0:
+
+                FValue[iunder[t]][junder[t]] = np.random.normal(
+                    loc=0,
+                    scale=np.sqrt(A / k[iunder[t]][junder[t]]**
+                                  (2 - (ns - 1) * 2. / 3.)) /
+                    np.sqrt(2)) + np.random.normal(
+                        loc=0,
+                        scale=np.sqrt(A / k[iunder[t]][junder[t]]**
+                                      (2 -
+                                       (ns - 1) * 2. / 3.)) / np.sqrt(2)) * 1j
+
+    else:
+        raise ValueError("limit must be None or tuple or list")
+
+    Value = realInverseFourier(FValue)
+
+    kx = (min(kx), max(kx))
+    ky = (min(ky), max(ky))
+
+    ValueE = (Value.min(), Value.max())
+
+    return ValueE, FValue, k, kx, ky

RealSpaceInterface/Calc2D/DataPropagation.py

@@ -0,0 +1,37 @@
+import numpy as np 
+#uses one dimensional interpolation
+def PropagateDatawithListOld(k,FValue,zredindex,transferFunctionlist):
+	return (transferFunctionlist[zredindex](k.ravel()) * FValue.ravel()).reshape(FValue.shape)
+
+def PropagateDatawithList(k, FValue, zredindex, transferFunctionlist):
+    result = {}
+    for field, transfer_function in transferFunctionlist.items():
+        result[field] = (transfer_function[zredindex](k.ravel()) * FValue.ravel()).reshape(FValue.shape)
+    return result
+
+#module with uses two dimensional interpolation and propagates all data at once (fastest but high memory consumption)
+def PropagateAllData(k,FValue,allzred,transferFunction):
+
+	allFValue = np.ones((len(allzred),FValue.shape[0],FValue.shape[1]),dtype=complex)
+
+	for kxindex in range(FValue.shape[0]):
+		allFValue[:,kxindex,:] = transferFunction(allzred,k[kxindex])*FValue[kxindex]
+
+
+	return allFValue
+
+
+#module with uses 2 dimensional interpolation (slowest but can be useful if the set of redshift changes very often)
+def PropagateData(k,FValue,zred,transferFunction):
+
+	FValuenew = np.ones(FValue.shape,dtype=complex)
+
+	for kxindex in range(FValue.shape[0]):
+		allFValue[kxindex,:] = transferFunction(zred,k[kxindex])*FValue[kxindex]
+
+
+	return allFValue
+
+
+
+

RealSpaceInterface/Calc2D/Database.py

@@ -0,0 +1,58 @@
+import pickle
+import os
+import logging
+import uuid
+
+class Database:
+    def __init__(self, directory, db_file="database.dat"):
+        self.directory = directory
+        self.db_file = db_file
+
+        if not os.path.isdir(directory):
+            raise ValueError("'{}' is not a directory!".format(directory))
+
+        self.db_path = os.path.join(directory, db_file)
+        if not os.path.exists(self.db_path):
+            logging.info("No database found; Creating one at {}.".format(self.db_path))
+            with open(self.db_path, "w") as f:
+                pickle.dump(dict(), f)
+
+        self.db = self.__read_database()
+
+    def __read_database(self):
+        with open(self.db_path) as f:
+            return pickle.load(f)
+
+    def __write_database(self):
+        with open(self.db_path, "w") as f:
+            pickle.dump(self.db, f)
+
+    def __create_file(self, data):
+        filename = str(uuid.uuid4())
+        with open(os.path.join(self.directory, filename), "w") as f:
+            pickle.dump(data, f)
+        return filename
+
+    def __get_frozen_key(self, key):
+        return frozenset(key.items())
+
+    def __getitem__(self, key):
+        frozen_key = self.__get_frozen_key(key)
+        if frozen_key in self.db:
+            filename = self.db[frozen_key]
+            with open(os.path.join(self.directory, filename)) as f:
+                return pickle.load(f)
+        else:
+            raise KeyError("No data for key: {}".format(key))
+
+    def __setitem__(self, key, data):
+        frozen_key = self.__get_frozen_key(key)
+        self.db[frozen_key] = self.__create_file(data)
+        self.__write_database()
+
+    def __contains__(self, key):
+        """
+        Return whether `self` contains a record
+        for the given `key`.
+        """
+        return self.__get_frozen_key(key) in self.db