pushing the Real Space Interface

RealSpaceInterface/Calc2D/CalculationClass.py

+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

+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

+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

+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
\ No newline at end of file

RealSpaceInterface/Calc2D/TransferFunction.py

+import os.path
+import pickle
+import uuid
+import numpy as np
+from scipy.interpolate import InterpolatedUnivariateSpline, RectBivariateSpline
+import sys
+import logging
+
+from classy import Class
+
+import Calc2D.Database as Database
+import config
+
+TRANSFER_QUANTITIES = ["d_g", "d_ur", "d_cdm", "d_b", "d_g/4 + psi"]
+
+def ComputeTransferData(settings, redshift):
+    database_key = settings.copy()
+    database_key.update({'redshift': tuple(redshift)})
+
+    database = Database.Database(config.DATABASE_DIR)
+    if database_key in database:
+        return database[database_key], redshift
+    else:
+        cosmo = Class()
+        cosmo.set(settings)
+        cosmo.compute()
+
+        outputData = [cosmo.get_transfer(z) for z in redshift]
+        # Calculate d_g/4+psi
+        for transfer_function_dict in outputData:
+            transfer_function_dict["d_g/4 + psi"] = transfer_function_dict["d_g"]/4 + transfer_function_dict["psi"]
+        # Now filter the relevant fields
+        fields = TRANSFER_QUANTITIES + ["k (h/Mpc)"]
+        outputData = [{field: outputData[i][field] for field in fields} for i in range(len(redshift))]
+
+        database[database_key] = outputData
+        return outputData, redshift
+
+
+def ComputeTransferFunctionList(cosmologicalParameters, redshift, kperdecade=200, P_k_max=100):
+    class_settings = cosmologicalParameters.copy()
+    class_settings.update({
+        "output": "mTk",
+        "gauge": "newtonian",
+        "evolver": "1",
+        "P_k_max_h/Mpc": P_k_max,
+        "k_per_decade_for_pk": kperdecade,
+        "z_max_pk": str(max(redshift)),
+    })
+
+    data_dict, redshift = ComputeTransferData(class_settings, redshift)
+    transfer_functions = {field: [] for field in TRANSFER_QUANTITIES}
+
+
+    for i in range(len(redshift)):
+        k_data = data_dict[0]["k (h/Mpc)"] * cosmologicalParameters["h"]  #in order to get k [1/Mpc]
+        k_data_zero = np.concatenate(([0.0], k_data))
+        for field in TRANSFER_QUANTITIES:
+            data = data_dict[i][field] / data_dict[i][field][0]
+            data_zero = np.concatenate(([1.0], data))
+            interpolated_func = InterpolatedUnivariateSpline(k_data_zero, data_zero)
+            transfer_functions[field].append(interpolated_func)
+
+    return transfer_functions

RealSpaceInterface/Calc2D/rFourier.py

+import numpy as np
+import numpy.fft as fft
+
+def realFourier(step, Value):
+    FValue = np.fft.fftshift(
+        np.fft.rfft2(Value), axes=(0))  #shifting only the x axes
+
+    kx = np.fft.fftshift(np.fft.fftfreq(Value.shape[0], d=step)) * 2 * np.pi
+    ky = np.fft.rfftfreq(Value.shape[0], d=step) * 2 * np.pi
+
+    return kx, ky, FValue
+
+def realInverseFourier(FValue):
+    return np.fft.irfft2(np.fft.ifftshift(
+        FValue, axes=(0)))  #shifting only on the x axes
+
+
+def realInverseAllFourier(allFValue):
+    return np.fft.irfftn(
+        np.fft.ifftshift(allFValue, axes=(1)),
+        axes=(1, 2))  #shifting only on the x axes

RealSpaceInterface/README

+For installation of python packages, run
+
+    pip install -r requirements.txt
+
+Launch the application with
+
+    python tornadoserver.py
+
+Then in any browser open the URL
+
+    http://localhost:7777
+
+------------------------------------------------------------
+
+Cache files are located in cache/, so to clear the cache, run
+
+    rm cache/*
+

RealSpaceInterface/colormap_converter.py

+import matplotlib.cm as cm
+import matplotlib.pyplot as plt
+import numpy as np
+from PIL import Image
+import os
+
+OUTPUT_DIR = os.path.join("static", "images", "colormaps")
+WIDTH = 512
+
+def create_image(cmap, width):
+    values = np.linspace(0, 1, width)
+    colors = cmap(values).reshape((1, width, 4))
+    image = Image.fromarray(np.uint8(255 * colors))
+    return image
+
+cmap_names = {}
+cmap_names['Uniform'] = [
+            'viridis', 'plasma', 'inferno', 'magma']
+cmap_names['Diverging'] = [
+            'seismic', 'RdYlBu', 'Spectral'
+            ]
+cmap_names['Miscellaneous'] = ['jet']
+
+if __name__ == "__main__":
+    for category in cmap_names:
+        category_dir = os.path.join(OUTPUT_DIR, category)
+        if not os.path.exists(category_dir):
+            os.mkdir(category_dir)
+        for name in cmap_names[category]:
+            result = create_image(plt.get_cmap(name), width=WIDTH)
+            output_path = os.path.join(category_dir, "{}.png".format(name))
+            print(output_path)
+            result.save(output_path)
+

RealSpaceInterface/config.py

+import os
+
+# Default port number to listen on. Can be overriden by passing a port number
+# as the first command line argument, e.g. `python tornadoserver.py 1234`
+PORT = 7777
+
+# Directory to store previously computed transfer functions, spectra etc. in
+DATABASE_DIR = "cache"
+
+# Maximum number of thread pool workers (only required for multi-user usage)
+MAX_THREADPOOL_WORKERS = 8
+
+# Path of colormap directory relative to the static directory from which
+# tornado serves static files
+COLORMAP_PATH = os.path.join("images", "colormaps")
+
+# number of sample points for the transfer function that is displayed
+# in the client
+TRANSFER_FUNCTION_CLIENT_SAMPLES = 400
+
+# number of sample points for the matter spectrum that is displayed
+# in the client per decade
+MATTER_SPECTRUM_CLIENT_SAMPLES_PER_DECADE = 40

RealSpaceInterface/requirements.txt

+futures>=3.0.2
+numpy>=1.8.2
+scipy>=0.14.0
+tornado>=4.0.2
+opencv-python
+# Building the Python bindings of CLASS requires Cython
+Cython