Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Decide to drop ODEs

  • Loading branch information...
commit 0a82a5a6356c76c9289fe577f99b42ceb1f1ae25 1 parent 925f28a
@dagss dagss authored
Showing with 0 additions and 2,440 deletions.
  1. +0 −432 einsteinboltzmann/cl_mod.pyx
  2. BIN  einsteinboltzmann/equations.pdf
  3. +0 −81 einsteinboltzmann/equations.tex
  4. +0 −407 einsteinboltzmann/evolution_mod.py
  5. 0  einsteinboltzmann/gsl/__init__.py
  6. 0  einsteinboltzmann/gsl/__init__.pyx
  7. +0 −12 einsteinboltzmann/gsl/error.pxd
  8. +0 −20 einsteinboltzmann/gsl/error.pyx
  9. +0 −189 einsteinboltzmann/gsl/ode.pxd
  10. +0 −3  einsteinboltzmann/gsl/ode.pyx
  11. +0 −21 einsteinboltzmann/gsl/special.pyx
  12. +0 −127 einsteinboltzmann/gsl/spline.pxd
  13. +0 −258 einsteinboltzmann/gsl/spline.pyx
  14. +0 −7 einsteinboltzmann/gsl/wscript
  15. +0 −29 einsteinboltzmann/params.py
  16. +0 −3  einsteinboltzmann/rec_mod.pxd
  17. +0 −230 einsteinboltzmann/rec_mod.pyx
  18. +0 −55 einsteinboltzmann/run.py
  19. +0 −4 einsteinboltzmann/time_mod.pxd
  20. +0 −132 einsteinboltzmann/time_mod.pyx
  21. +0 −43 einsteinboltzmann/tools/cython.py
  22. +0 −181 einsteinboltzmann/tools/fwrapktp.py
  23. +0 −64 einsteinboltzmann/tools/fwraptool.py
  24. +0 −58 einsteinboltzmann/tools/inplace.py
  25. +0 −29 einsteinboltzmann/tools/numpy.py
  26. +0 −10 einsteinboltzmann/utils.py
  27. BIN  einsteinboltzmann/waf
  28. +0 −45 einsteinboltzmann/wscript
View
432 einsteinboltzmann/cl_mod.pyx
@@ -1,432 +0,0 @@
-from __future__ import division
-
-cdef double pi, c, x_start_rec, x_end_rec
-
-import numpy as np
-from numpy import pi, exp
-from time_mod import get_H_p, get_dH_p, get_eta, get_H, x_start_rec, x_end_rec
-from rec_mod import get_g, get_dg, get_ddg, get_tau, get_dtau, get_ddtau
-from evolution_mod import compute_einstein_boltzmann_grid, k_min, k_max, x_init
-from time_mod cimport get_H_p_fast
-import os
-from time import clock
-import cPickle as pickle
-from gsl.spline cimport Spline, MultiSpline
-cimport gsl.spline
-import gsl.spline
-cimport numpy as np
-from params import Omega_b, Omega_m, Omega_r, H_0, c
-from libc.math cimport fabs
-
-
-ctypedef int (*spline_eval_func_type)(void *, double* xa, double* ya, size_t size, double x,
- gsl.spline.gsl_interp_accel *, double * y)
-
-cimport cython
-
-#
-# More zeroth order stuff related to HH'. Logically belongs to
-# time_mod perhaps.
-#
-
-def get_H_dH(x):
- return - H_0**2 / 2 * (3 * (Omega_m + Omega_b) * exp(-3*x)
- + 4 * Omega_r * exp(-4*x))
-
-def get_H_dH_derivative(x):
- return H_0**2 / 2 * (9 * (Omega_m + Omega_b) * exp(-3*x)
- + 16 * Omega_r * exp(-4*x))
-
-def get_H_p_dH_p(x):
- return exp(2*x) * (get_H(x)**2 + get_H_dH(x))
-
-def get_H_p_dH_p_derivative(x):
- return exp(2*x) * (2 * get_H(x)**2
- + 4 * get_H_dH(x)
- + get_H_dH_derivative(x))
-
-#
-# Grid evaluation and resampling
-#
-
-def resample(values, oldgrid, newgrid, out=None):
- """
- Resamples a 1D function evaluated in oldgrid onto newgrid, by
- creating a spline and evaluating it.
- """
- spline = Spline(oldgrid, values)
- return spline.evaluate(newgrid, out=out)
-
-def resample_S(S_values, old_k_grid, new_k_grid, old_x_grid, x_nums, new_x_grids):
- """
- Resamples any function in (k,x), such as S, to a new (presumably
- higher) resolution in k, and then resamples along the x axis
- according to provided grids. x_grid[ik, :x_nums[ik]] is the
- x-grid for k value ik.
- """
- assert S_values.shape == (old_k_grid.shape[0], old_x_grid.shape[0])
-
- cdef int ik, ix
- # Resample k
- tmp = np.zeros(shape=(new_k_grid.shape[0], old_x_grid.shape[0]))
- for ix in range(old_x_grid.shape[0]):
- tmp[:, ix] = resample(S_values[:, ix].copy(), old_k_grid, new_k_grid)
- # Resample x
- out = np.zeros(shape=(new_k_grid.shape[0], np.max(x_nums)))
- for ik in range(new_k_grid.shape[0]):
- x_grid = new_x_grids[ik, :x_nums[ik]]
- resample(tmp[ik, :], old_x_grid, x_grid, out=out[ik, :x_nums[ik]])
- return out
-
-def make_x_grids(np.ndarray[double] k_grid, double N_per_wave=10):
- """
- Creates grids in x for each value of k. The start point is
- x_start_rec, and endpoint is today.
-
- Each k will yield have a different number of grid points.
- The grid is created so that N_per_wave evaluation points
- are made per oscillation of the spherical bessel function.
- Still, there's a minimum of 300 points during recombination
- and 200 points after, which will capture any features of the
- source function.
-
- Returns (x_nums, x_grids), where x_grid[ik, :x_nums[ik]]
- contains the grid for ik.
- """
- import evolution_mod
-
- cdef double delta_x, x, k, min_delta_x_rec, min_delta_x_after
- cdef int ix, ik, Nx
- cdef np.ndarray[double, ndim=2] x_grid
- cdef np.ndarray[int] x_nums = np.zeros(k_grid.shape[0], dtype=np.intc)
-
- min_delta_x_rec = (x_end_rec - x_start_rec) / 200
- min_delta_x_after = - x_end_rec / 300
-
- # Overallocate x_grid
- delta_x = 2 * pi * get_H_p_fast(0) / c / np.max(k_grid) / N_per_wave
- x_grid = np.zeros((k_grid.shape[0], np.max(
- [np.round(-x_start_rec / delta_x), 510])), dtype=np.float64)
-
- for ik in range(k_grid.shape[0]):
- k = k_grid[ik]
- # Start from end, then dynamically increase delta_x.
- # Finally reverse the resulting grid (we don't know
- # how big it will be up front).
- x_grid[ik, 0] = x = 0
- ix = 1
- while x > x_start_rec:
- delta_x = 2 * pi * get_H_p_fast(x) / c / k / N_per_wave
- if x < x_end_rec:
- delta_x = delta_x if delta_x < min_delta_x_rec else min_delta_x_rec
- else:
- delta_x = delta_x if delta_x < min_delta_x_after else min_delta_x_after
- x -= delta_x
- x_grid[ik, ix] = x
- ix += 1
- # Store grid size
- x_nums[ik] = Nx = ix
- # Reverse grid
- for ix in range(0, Nx // 2):
- x_grid[ik, ix], x_grid[ik, Nx - ix - 1] = x_grid[ik, Nx - ix - 1], x_grid[ik, ix]
-
- return x_nums, x_grid[:, :np.max(x_nums)].copy()
-
-def compute_S(eb):
- """
- Computes the source function S (=Stilde in Callin).
- eb should be an EinsteinBoltzmannSolution object
- from evolution_mod.
- """
- x_grid = eb.x_grid
- k_grid = eb.k_grid
-
- # All variables here will be arrays with k along the first
- # axis and x along the second. Variables such as g which only
- # depend on one variable is stored in a 1-by-n or n-by-1 array
- # which then automatically "broadcast" (repeats itself to match).
- # arr[None, :] reshapes an n-array to a 1-by-n-array.
- H_p_dH_p = get_H_p_dH_p(x_grid)[None, :]
- H_p_dH_p_derivative = get_H_p_dH_p_derivative(x_grid)[None, :]
- H_p = get_H_p(x_grid)[None, :]
- dH_p = get_dH_p(x_grid)[None, :]
- tau = get_tau(x_grid)[None, :]
- dtau = get_dtau(x_grid)[None, :]
- ddtau = get_ddtau(x_grid)[None, :]
- np.seterr(all='ignore')
- g = get_g(x_grid)[None, :]
- dg = get_dg(x_grid)[None, :]
- ddg = get_ddg(x_grid)[None, :]
- k = k_grid[:, None]
- Theta = eb.Theta
- dTheta = eb.dTheta
-
- ddTheta2 = (+ (2/5) * c * k / H_p * (-(dH_p / H_p) * Theta[1, :, :] + dTheta[1, :, :])
- - (3/5) * c * k / H_p * (-(dH_p / H_p) * Theta[3, :, :] + dTheta[3, :, :])
- + (9/10) * (ddtau * Theta[2, :, :] + dtau * dTheta[2, :, :]))
-
- # First term
- A = g * (eb.Theta[0, :, :] + eb.Psi + (1/4)*eb.Theta[2, :, :])
-
- # Second term
- B = np.exp(-tau) * (eb.dPsi - eb.dPhi)
-
- # Third term
- C = -1/(c*k) * (dH_p * g * eb.v_b + H_p * dg * eb.v_b + H_p * g * eb.dv_b) #* -10
-
- # Fourth term
- D = 3/4/(c*k)**2 * (
- H_p_dH_p_derivative * g * Theta[2, :, :]
- + 3 * H_p_dH_p * (dg * Theta[2, :, :] + g * dTheta[2, :, :])
- + H_p**2 * (ddg * Theta[2, :, :] + 2 * dg * dTheta[2, :, :]
- + g * ddTheta2))
-
- S = A + B + C + D
- return S
-
-_spherical_bessel_cache = None # in-session cache
-def get_spherical_bessels(zmax, lmax=1200, cache_filename='sphbessel.pickle',
- double atol=1e-30, N_per_wave=5):
- """
- Function to compute splines for the spherical bessel functions.
-
- Returns (splines, zmins), each of length lmax+1, the first is a list
- of spline objects and the second an array of zmin values (the functions
- are 0 before this point).
-
- atol is used to determine whether the function is 0, for determining
- zmin, only.
-
- The resulting splines are cached to the file given by cache_filename.
- zmax should be the largest parameter needed. Changing either zmax
- or lmax to something larger than what is present in cache, the
- cache is recomputed.
- """
- global _spherical_bessel_cache
- import cPickle as pickle
- from sphbess import sphbes_sj
-
- cdef np.ndarray[double, ndim=2, mode='c'] j_values
- cdef Py_ssize_t l, zi
-
- zmax += 10 # cut some slack
-
- cache = _spherical_bessel_cache
- if cache is None and os.path.exists(cache_filename):
- print 'Loading spherical bessel functions from cache %s...' % cache_filename
- with file(cache_filename, 'r') as f:
- cache = pickle.load(f)
-
- if cache is not None:
- if cache['lmax'] < lmax:
- print 'Cached for lmax=%d, need for lmax=%d, recomputing' % (cache['lmax'], lmax)
- elif cache['zmax'] < zmax:
- print 'Cached for zmax=%e, need for zmax=%e, recomputing' % (cache['zmax'], zmax)
- else:
- # We're done, return the cached splines
- return cache['splines'], cache['zmins']
-
- print 'Computing spherical bessel functions for lmax=%d, zmax=%e' % (lmax, zmax)
- # Need to (re)compute.
- # Choose our delta_z by the argument from Callin. We select the number of
- # samples dynamically in response to zmax, so only grid step size matter
- delta_z = 2 * pi / N_per_wave
-
- # Sample function at every combination of l and z
- z_grid = np.arange(0, zmax, delta_z) # from 0 to zmax by delta_z
- ls = np.arange(lmax + 1)
- t0 = clock()
- j_values = sphbes_sj(ls[:, None], z_grid[None, :])
- print 'Done in %.3e secs' % (clock() - t0)
-
- # Now, find the zmins arrays
- cdef np.ndarray[double, mode='c'] zmins = np.zeros(lmax + 1, np.double)
- zmins[0] = 0 # first bessel starts at 1
- for l in range(1, lmax + 1):
- for zi in range(1, z_grid.shape[0]):
- if fabs(j_values[l, zi]) > atol:
- zmins[l] = z_grid[zi - 1]
- break
-
- # Make splines
- splines = [Spline(z_grid, j_values[l, :]) for l in range(lmax + 1)]
-
- # Finally, cache to file and return result
- cache = dict(splines=splines, zmins=zmins, lmax=lmax, zmax=zmax)
- _spherical_bessel_cache = cache
- with file(cache_filename, 'w') as f:
- pickle.dump(cache, f, protocol=2)
- return splines, zmins
-
-class PowerSpectrumIntegrator:
- """
- Class for computing the power spectrum. We use a class to store
- the various dynamic grids etc. and their derived quantities (which
- are precomputed to speed things up and therefore adds some complexity,
- which we handle through using a class).
-
- eb - If provided, contains a prevously computed solution to the
- Einstein-Boltzmann equations (otherwise they'll be computed
- when constructing the object)
-
- n_k - Number of k values to use
-
- lmax - Maximum l value queried
- """
-
- def __init__(self, eb=None, n_k=1700, k_grid=None, lmax=1200):
- if eb is None:
- # Must compute E-B solutions. Use coarse grid described in Callin.
- import time_mod
- print 'Computing E-B solutions...'
- t0 = clock()
- eb = compute_einstein_boltzmann_grid(x_grid=time_mod.x_grid)
- print 'done in %.2f secs!' % (clock() - t0)
- if k_grid is None:
- k_grid = np.linspace(eb.k_grid[0], eb.k_grid[-1], n_k)
- self.eb = eb
- self.k_grid = k_grid
- self.lmax = lmax
-
- # Create dynamic grid
- self.x_nums, self.x_grids = make_x_grids(k_grid, N_per_wave=10)
-
- # Compute source function
- S_lowres = compute_S(eb)
- self.S = resample_S(S_lowres, eb.k_grid, self.k_grid, eb.x_grid,
- self.x_nums, self.x_grids)
-
- # Compute dynamic-grid-derived quantities: horizon_delta and dx.
- # These correspond to x_grids w.r.t. indexing
- self.horizon_delta = get_eta(0) - get_eta(self.x_grids)
- self.dxs = self.x_grids[:, 1:] - self.x_grids[:, :-1]
-
- # Get spherical Bessel functions
- self.j_splines, self.j_zmins = get_spherical_bessels(lmax=lmax,
- zmax=get_eta(0) * k_max)
-
- # Done, everything set up for computing Theta_k and C_l
-
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- def compute_Theta_k(self, int l, np.ndarray[double, mode='c'] out=None):
- """
- Grid-evaluate Theta_l(k) for all values in our k_grid, for the l given.
-
- This is done through line-of-sight integration.
- """
- cdef Py_ssize_t ik, ix, nx
- cdef double k, I, jval, z, prev_dx, dx, prev_x, x
-
- cdef Spline j_spline = self.j_splines[l]
- cdef double zmin = self.j_zmins[l]
- cdef np.ndarray[double, mode='c'] k_grid = self.k_grid
-
- cdef np.ndarray[double, ndim=2, mode='c'] S = self.S
- cdef np.ndarray[double, ndim=2, mode='c'] horizon_delta = self.horizon_delta
- cdef np.ndarray[double, ndim=2, mode='c'] dxs = self.dxs
- cdef np.ndarray[int, ndim=1, mode='c'] x_nums = self.x_nums
-
- # Here follows some pretty low-level stuff which speeds up
- # spline evaluation in this particular case. j_eval is the
- # spline evaluation routine pointer.
- cdef gsl.spline.gsl_interp_accel* spline_acc = gsl.spline.gsl_interp_accel_alloc()
- cdef spline_eval_func_type j_eval = j_spline.interp.type.eval
- cdef void* j_state = j_spline.interp.state
- cdef size_t j_len = j_spline.interp.size
- cdef double* xa = j_spline.xa
- cdef double* ya = j_spline.ya
-
- if out is None:
- out = np.zeros_like(k_grid)
-
- # Evaluate the integral as we go, by using the trapezoidal rule
- # with *non-uniform* grid cell sizes
- for ik in range(k_grid.shape[0]):
- k = k_grid[ik]
- nx = x_nums[ik]
- I = 0
-
- if k * horizon_delta[ik, 0] < zmin:
- # Wohoo, bessel is all zero
- out[ik] = 0
- continue # In Python/C/Java/everything but Fortran,
- # continue means "next loop iteration"
-
- # Evaluate first grid point
- z = k * horizon_delta[ik, 0]
- j_eval(j_state, xa, ya, j_len, z, spline_acc, &jval)
- I += (S[ik, 0] * jval) * dxs[ik, 0]
-
- # Inner grid points
- for ix in range(1, nx - 1):
- z = k * horizon_delta[ik, ix]
- if z < zmin:
- continue
- j_eval(j_state, xa, ya, j_len, z, spline_acc, &jval)
- I += (S[ik, ix] * jval) * (dxs[ik, ix - 1] + dxs[ik, ix])
-
- # Last grid point
- z = k * horizon_delta[ik, nx - 1]
- j_eval(j_state, xa, ya, j_len, z, spline_acc, &jval)
- I += (S[ik, nx - 1] * jval) * dxs[ik, nx - 1 - 1]
-
- out[ik] = I * 0.5
-
- gsl.spline.gsl_interp_accel_free(spline_acc)
- return out
-
-
-
- def compute_power_spectrum(self, n_s=1, ls=None):
- """
- Computes the power spectrum. n_s is the tilt parameter from inflation.
-
- High-level driver function; to get and plot intermediate results
- one can call the lower-level functions.
-
- ls is a grid at which to do the actual computation. This will
- be splined before it is returned.
-
- The resulting spectrum is unscaled. You may want to multiply it with
- l*(l+1) (and rescale amplitude) prior to plotting.
-
- Return value: (l, Cl), each going from 2..max(ls)
- """
- from scipy.integrate import simps
-
- if ls is None:
- ls = [ 2, 3, 4, 6, 8, 10, 12, 15, 20, 30, 40, 50, 60, 70, 80, 90, 100,
- 120, 140, 160, 180, 200, 225, 250, 275, 300, 350, 400, 450, 500, 550,
- 600, 650, 700, 750, 800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200]
- lmax = max(ls)
- if lmax > self.lmax:
- raise ValueError('lmax too high (Bessel functions not computed)')
- Theta_k_lores = np.zeros(self.k_grid.shape[0])
- Theta_k_hires = np.zeros(5000)
- k_grid_hires = np.linspace(self.k_grid[0], self.k_grid[-1], 5000)
- integrand_factor = (c * k_grid_hires / H_0)**(n_s - 1) / k_grid_hires
-
- # Sample Cls
- Cl_samples = np.zeros(len(ls))
- t0 = clock()
- print 'Computing power spectrum'
- for idx, l in enumerate(ls):
- self.compute_Theta_k(l, out=Theta_k_lores)
- resample(Theta_k_lores, self.k_grid, k_grid_hires, out=Theta_k_hires)
- integrand = Theta_k_hires **2 * integrand_factor
-
- # Use SciPy's simpson's method
- Cl_samples[idx] = simps(integrand, x=k_grid_hires)
-
- print 'Time taken: %.3f' % (clock() - t0)
-
- # Resample Cls
- ls_hires = np.arange(2, lmax+1)
- Cl = resample(Cl_samples, np.array(ls, np.double), np.array(ls_hires, np.double))
- return ls_hires, Cl
-
-
-
-
View
BIN  einsteinboltzmann/equations.pdf
Binary file not shown
View
81 einsteinboltzmann/equations.tex
@@ -1,81 +0,0 @@
-\documentclass[handout]{beamer}
-
-%\setbeameroption{show notes}
-
-% Vary the color applet (try out your own if you like)
-%\colorlet{structure}{red!65!black}
-
-%\beamertemplateshadingbackground{yellow!100}{white}
-
-%\usepackage{beamerthemesplit}
-\usepackage{tikz}
-\usepackage{graphics}
-\usepackage{graphicx}
-\usepackage{hyperref}
-\usepackage{setspace}
-\usepackage{colortbl}
-\usepackage{listings}
-\usepackage{amsfonts}
-\usepackage{amsmath}
-\usepackage{multicol}
-
-%\usetikzlibrary{shapes}
-
-
-\usetheme{Madrid}
-\useinnertheme{rounded}
-\usecolortheme{crane}
-
-\definecolor{dredcolor}{rgb}{.5,0.0,0.0}
-\definecolor{blackcolor}{rgb}{0.0,0.0,0.0}
-\definecolor{dgreencolor}{rgb}{0.0,0.4,0.0}
-\definecolor{dbluecolor}{rgb}{.02,.02,.808}
-\newcommand{\dblue}{\color{dbluecolor}\bf}
-\newcommand{\dgreen}{\color{dgreencolor}\bf}
-
-
-\title[Cython]{Cython Tutorial}
-\author[Seljebotn]{Dag Sverre Seljebotn \\ \texttt{dagss@student.matnat.uio.no}}
-\date{Open Cython Day, Munich}
-\institute{University of Oslo}
-
-%\pgfdeclaremask{fsu}{fsu_logo_ybkgrd}
-%\pgfdeclareimage[mask=fsu,width=1cm]{fsu-logo}{fsu_logo_ybkgrd}
-
-%\logo{\vbox{\vskip0.1cm\hbox{\pgfuseimage{fsu-logo}}}}
-
-\begin{document}
-\section[Examples]{Practical examples}
-\frame{
-\includegraphics[width=\textwidth]{V1masked.png} \\
-{\tiny Perturbations to the Cosmic Microwave Background
-(http://lambda.gsfc.nasa.gov)}
-
-}
-\frame{
-\def\H{\mathcal{H}}
-
-\small
-\begin{align*}
-\Theta'_0 &= -\frac{ck}{\H} \Theta_1 - \Phi' \\
-\Theta'_1 &= \frac{ck}{3\H} \Theta_0 - \frac{2ck}{3\H}\Theta_2 +
-\frac{ck}{3\H}\Psi + \tau'\left[\Theta_1 + \frac{1}{3}v_b\right], \\
-\Theta'_\ell &= \frac{\ell ck}{(2\ell+1)\H}\Theta_{\ell-1} - \frac{(\ell+1)ck}{(2\ell+1)\H}
-\Theta_{\ell+1} + \tau'\left[\Theta_\ell - \frac{1}{10}\Theta_\ell
- \delta_{\ell,2}\right], \quad \ell \ge 2 \\
-\delta' &= \frac{ck}{\H} v - 3\Phi'
- \quad\quad
-\\ v' &= -v -\frac{ck}{\H} \Psi \\
-\delta_b' &= \frac{ck}{\H}v_b -3\Phi'
- \quad\quad \\
- v_b' &= -v_b - \frac{ck}{\H}\Psi + \tau'R(3\Theta_1 + v_b) \\
-\Phi' &= \Psi - \frac{c^2k^2}{3\H^2} \Phi + \frac{H_0^2}{2\H^2}
-\left[\Omega_m a^{-1} \delta + \Omega_b a^{-1} \delta_b + 4\Omega_r
- a^{-2}\Theta_0\right] \\
-\Psi &= -\Phi - \frac{12H_0^2}{c^2k^2a^2}\Omega_r\Theta_2, \quad\quad R = \frac{4\Omega_r}{3\Omega_b a}
-\end{align*}
-
-
-}
-
-\end{document}
View
407 einsteinboltzmann/evolution_mod.py
@@ -1,407 +0,0 @@
-from __future__ import division
-
-#
-# evolution_mod -- AST5220 Milestone 3
-#
-
-from params import c, H_0, Omega_b, Omega_r, Omega_m
-from time_mod import x_start_rec, get_eta, get_H, get_H_p, get_dH_p
-from rec_mod import get_dtau, get_ddtau
-import numpy as np
-from scipy.integrate import ode
-
-from numpy import inf, arange, array, abs, nan, pi
-from numpy import exp, log, sqrt
-np.seterr(all='raise') # fail hard on NaN and over/underflow
-
-#
-# Accuracy parameters
-#
-a_init = 1e-8
-k_min = 0.1 * H_0 / c
-k_max = 1e3 * H_0 / c
-n_k = 100
-lmax = 6
-# Our n_par here includes Psi, for program simplicity. This may
-# be changed depending on requirements of milestone 4. (The ODE
-# solving doesn't use n_par, only the final parameter computation).
-n_par = 6 + lmax + 1
-
-x_init = log(a_init)
-
-def get_tight_coupling_switch(x_grid, k):
- """
- Find the index in x_grid at which one should switch from tight
- coupling to full equations. Simply do linear searches.
- """
- dtau_values = get_dtau(x_grid)
- idx_cond1 = np.nonzero(-dtau_values < 10)[0][0]
- idx_cond2 = np.nonzero(x_grid >= x_start_rec)[0][0]
- idx_cond3 = np.nonzero(get_H_p(x_grid) * dtau_values > -10 * c * k)[0][0]
- return min([idx_cond1, idx_cond2, idx_cond3])
-
-#
-# ODE system
-#
-
-i_delta = 0
-i_delta_b = 1
-i_v = 2
-i_v_b = 3
-i_Phi = 4
-i_Theta = 5
-i_Theta0 = 5
-i_Theta1 = 6
-i_Theta2 = 7
-
-class EinsteinBoltzmannEquations:
- """
- Evaluates the Einstein-Boltzmann equations for given k, lmax,
- both in tight coupling regime and using full equation set.
-
- The Jacobian is only defined for the full equation set, since
- there are so few oscillations during tight coupling.
-
- Statistics about how often the function and the Jacobian is
- evaluated is available through f_count and jacobian_count
- (and is the reason for making this a class).
- """
-
- def __init__(self, tight_coupling, k, lmax):
- if tight_coupling:
- n = 7 # Only include up to Theta_1
- else:
- n = 6 + lmax
- self.dy = np.zeros(n)
- self.jacobian_array = np.zeros((n,n))
- self.tight_coupling = tight_coupling
- self.k = k
- self._lmax = lmax
- self.f_count = self.jacobian_count = 0
-
- def f(self, x, y):
- """
- Evaluate the derivative vector/right hand side of the E-B ODE.
- """
- dy = self.dy
- k = self.k
- tight_coupling = self.tight_coupling
- lmax = self._lmax
-
- self.f_count += 1
-
- # Fetch time-dependant variables
- a = exp(x)
- H = get_H(x)
- H_p = a * H
- dtau = get_dtau(x)
- ddtau = get_ddtau(x)
- R = 4 / 3 * Omega_r / Omega_b / a
- if not tight_coupling:
- eta = get_eta(x)
-
- # Parse the y array
- delta = y[i_delta]
- delta_b = y[i_delta_b]
- v = y[i_v]
- v_b = y[i_v_b]
- Phi = y[i_Phi]
- Theta0 = y[i_Theta0]
- Theta1 = y[i_Theta1]
- if tight_coupling:
- Theta2 = - 20 / 45 * c * k / H_p / dtau * Theta1
- else:
- Theta2 = y[i_Theta2]
-
- # Evaluate components
- Psi = -Phi - 12 * H_0**2 / c**2 / k**2 / a**2 * Omega_r * Theta2
- dy[i_Phi] = dPhi = (Psi
- - (c * k / H_p)**2 * Phi / 3
- + H_0**2 / H_p**2 / 2
- * (Omega_m / a * delta
- + Omega_b / a * delta_b
- + 4 * Omega_r / a**2 * Theta0))
- dy[i_Theta0] = dTheta0 = -(c * k / H_p) * y[i_Theta + 1] - dPhi
-
- if tight_coupling:
- q = ( (-((1 - 2*R) * dtau + (1 + R) * ddtau) * (3*Theta1 + v_b)
- - (c * k / H_p) * Psi
- + (1 - (H_p / H)) * (c * k / H_p) * (-Theta0 + 2*Theta2)
- - (c * k / H_p) * dTheta0)
- / ((1 + R) * dtau + (H_p / H) - 1))
- dy[i_v_b] = dv_b = ((-v_b - (c * k / H_p) * Psi
- + R * (q + (c * k / H_p) * (-Theta0 + 2*Theta2)
- - (c * k / H_p) * Psi)) / (1 + R))
- else:
- dy[i_v_b] = dv_b = -v_b - (c * k / H_p) * Psi + dtau * R * (3*Theta1 + v_b)
-
- dy[i_delta_b] = (c * k / H_p) * v_b - 3 * dPhi
- dy[i_v] = -v - (c * k / H_p) * Psi
- dy[i_delta] = (c * k / H_p) * v - 3 * dPhi
-
- if tight_coupling:
- dy[i_Theta + 1] = 1/3 * (q - dv_b)
- else:
- dy[i_Theta + 1] = ((c * k / H_p) / 3 * (Theta0 - 2*Theta2 + Psi)
- + dtau * (Theta1 + v_b / 3))
- for l in range(2, lmax): # Python: Endpoint is exclusive!
- dy[i_Theta + l] = (
- l / (2*l + 1) * (c * k / H_p) * y[i_Theta + l - 1]
- - (l+1)/(2*l + 1) * (c * k / H_p) * y[i_Theta + l + 1]
- + dtau * (y[i_Theta + l] - 0.1 * y[i_Theta + l] * (l == 2)))
- dy[i_Theta + lmax] = ((c * k / H_p) * y[i_Theta + lmax - 1]
- - (lmax+1) * c / H_p / eta * y[i_Theta + lmax]
- + dtau * y[i_Theta + lmax])
- return dy
-
- def jacobian(self, x, y):
- """
- Evaluate the Jacobian of the E-B ODE.
- """
- jac = self.jacobian_array
- k = self.k
- lmax = self._lmax
-
- self.jacobian_count += 1
-
- n = jac.shape[0]
- if self.tight_coupling:
- raise NotImplementedError("Only bothered to define Jacobian after tight coupling")
-
- # Fetch time-dependant variables
- a = exp(x)
- H = get_H(x)
- H_p = a * H
- dtau = get_dtau(x)
- ddtau = get_ddtau(x)
- R = 4 / 3 * Omega_r / Omega_b / a
- eta = get_eta(x)
-
- # Evaluate Jacobian and store it in
- # jac[derivative_of_param, with_respect_to_param]
- # The zero elements will always stay zero from the first initialization in
- # __init__.
-
- partial_dPsi_Phi = -1
- partial_dPsi_Theta2 = -12*(H_0 / c / k / a)**2 * Omega_r
-
- jac[i_Phi, i_Phi] = -1 - (c*k/H_p)**2 / 3
- tmp = (H_0/H_p)**2 / 2
- jac[i_Phi, i_delta] = tmp * Omega_m / a
- jac[i_Phi, i_delta_b] = tmp * Omega_b / a
- jac[i_Phi, i_Theta0] = tmp * 4 * Omega_r / a**2
- jac[i_Phi, i_Theta0] = (H_0/H_p)**2 * 2 * Omega_r / a**2
- jac[i_Phi, i_Theta2] = partial_dPsi_Theta2
-
- # delta, delta_b
- for i in range(n):
- jac[i_delta, i] = -3 * jac[i_Phi, i]
- jac[i_delta_b, i] = -3 * jac[i_Phi, i]
- jac[i_delta, i_v] = jac[i_delta_b, i_v_b] = c * k / H_p
-
- # v
- jac[i_v, i_v] = -1
- jac[i_v, i_Phi] = - c * k / H_p * partial_dPsi_Phi
- jac[i_v, i_Theta2] = - c * k / H_p * partial_dPsi_Theta2
-
- # v_b
- jac[i_v_b, i_v_b] = -1 + dtau * R
- jac[i_v_b, i_Theta1] = dtau * R * 3
- jac[i_v, i_Theta2] = - c * k / H_p * partial_dPsi_Theta2
- jac[i_v, i_Phi] = - c * k / H_p * partial_dPsi_Phi
-
- # Thetas
- for i in range(n):
- jac[i_Theta0, i] = -jac[i_Phi, i]
- jac[i_Theta0, i_Theta1] = -c * k / H_p
-
- jac[i_Theta1, i_Theta0] = c * k / H_p / 3
- jac[i_Theta1, i_Theta1] = dtau
- jac[i_Theta1, i_Theta2] = c * k / H_p * (- (2/3) + partial_dPsi_Theta2)
- jac[i_Theta1, i_v_b] = dtau / 3
-
- for l in range(2, lmax):
- jac[i_Theta + l, i_Theta + l - 1] = l / (2*l + 1) * c * k / H_p
- jac[i_Theta + l, i_Theta + l + 0] = - (l + 1) / (2*l + 1) * c * k / H_p
- jac[i_Theta + l, i_Theta + l + 1] = dtau * (1 - (l == 2) / 10)
-
- jac[i_Theta + lmax, i_Theta + lmax - 1] = c * k / H_p
- jac[i_Theta + lmax, i_Theta + lmax] = (l + 1) / H_p / eta + dtau
-
- return jac
-
-def solve_einstein_boltzmann(x_grid,
- k, rtol=1e-15,
- nsteps=10**10, max_step=None,
- min_step=None,
- use_jacobian=True,
- out_y=None, out_dy=None):
- """
- Solve the Einstein-Boltzmann equations for a given k.
-
- On return, out_y and out_dy contains the quantities and their
- derivatives. The first axis is the parameters and the second the
- x grid. The order of the parameters is:
-
- delta, delta_b, v, v_b, Phi, Theta_0, Theta_1, ..., Theta_lmax, Psi
-
- Return value: (out_y, out_dy, nf, nj), where nf and nj are the
- number of times the rhs and the Jacobian were evaluated.
-
- """
- i_Psi = i_Theta + lmax + 1
-
- assert np.isscalar(k)
- if out_y is None:
- out_y = np.zeros((n_par, x_grid.shape[0]))
- if out_dy is None:
- out_dy = np.zeros_like(out_y)
-
- # Options for the ODE solver
- ode_opts = dict(name='vode',
- method='bdf',
- rtol=rtol,
- max_step=max_step,
- # The default maximum steps is 500 per grid point, which is
- # way too low. It doesn't hurt to let this be really big
- # (set to small to get partial results when debugging)
- nsteps=nsteps)
-
- # Find tight coupling regime
- idx_switch = get_tight_coupling_switch(x_grid, k)
-
- # Initial conditions for tight coupling
- x_init_tight = x_grid[0]
- H_p = get_H_p(x_init_tight)
- Phi = 1
- delta = delta_b = Phi * 3 / 2
- v = v_b = 1/2 * Phi * (c * k / H_p)
- Theta0 = Phi / 2
- Theta1 = - 1/6 * Phi * (c * k / H_p)
- y_init_tight = np.array([delta, delta_b, v, v_b, Phi, Theta0, Theta1])
-
- # Integrate in tight coupling regime -- Theta[2:] is not included here
- rhs_tight = EinsteinBoltzmannEquations(tight_coupling=True, k=k, lmax=lmax)
- integrator = ode(rhs_tight.f)
- integrator.set_integrator(**ode_opts)
- integrator.set_initial_value(y_init_tight, x_init_tight)
- out_y[:i_Theta1 + 1, 0] = y_init_tight # Python: Slice endpoint not included
- out_dy[:i_Theta1 + 1, 0] = rhs_tight.f(x_init_tight, y_init_tight)
- for i in range(1, idx_switch): # idx_switch is not included!
- integrator.integrate(x_grid[i])
- if not integrator.successful():
- raise Exception()
- out_y[:i_Theta1 + 1, i] = integrator.y
- out_dy[:i_Theta1 + 1, i] = rhs_tight.f(x_grid[i], integrator.y)
-
- # Assign values to Theta[2:] and dTheta[2:] in the tight coupling
- # regime. This also gives initial conditions for the full equation
- # set. This is all array arithmetic, working on all values until
- # the switchpoint.
- x_grid_tight = x_grid[:idx_switch]
- H_p = get_H_p(x_grid_tight)
- dH_p = get_dH_p(x_grid_tight)
- dtau = get_dtau(x_grid_tight)
- ddtau = get_ddtau(x_grid_tight)
- # Python comment: Changing Theta below will also change out_y
- # accordingly, because Theta is a slice into out_y (NOT
- # a copy). Just like pointers in Fortran.
- Theta = out_y[i_Theta:, :idx_switch]
- dTheta = out_dy[i_Theta:, :idx_switch]
-
- Theta[2,:] = -20/45 * (c * k / H_p) / dtau * Theta[1,:]
- dTheta[2,:] = -20/45 * c * k * (
- dTheta[1,:] / (H_p * dtau)
- - Theta[1,:] * (dH_p * dtau + H_p * ddtau) / (H_p * dtau)**2)
- for l in range(2, lmax + 1):
- Theta[l,:] = - l/(2*l + 1) * (c * k / H_p) / dtau * Theta[l - 1,:]
- dTheta[l,:] = - l/(2*l + 1) * c * k * (
- dTheta[l-1,:] / (H_p * dtau)
- - Theta[l-1,:] * (dH_p * dtau + H_p * ddtau) / (H_p * dtau)**2)
-
- # Integrate the full equation set to the end. Remember to not include Psi
- # in the ODE.
- rhs_full = EinsteinBoltzmannEquations(tight_coupling=False, k=k, lmax=lmax)
- integrator = ode(rhs_full.f, None if not use_jacobian else rhs_full.jacobian)
- integrator.set_integrator(**ode_opts)
- integrator.set_initial_value(out_y[:i_Psi, idx_switch - 1], x_grid[idx_switch - 1])
- for i in range(idx_switch, x_grid.shape[0]):
- x = x_grid[i]
- integrator.integrate(x)
- if not integrator.successful():
- raise Exception()
- out_y[:i_Psi, i] = integrator.y
- out_dy[:i_Psi, i] = rhs_full.f(x_grid[i], integrator.y)
-
- # Finally, (re)compute Psi and its derivative. This may be removed
- # depending on needs of milestone 4.
- a = np.exp(x_grid)
-
- # Set up references ("pointers")
- Theta2 = out_y[i_Theta2, :]
- Phi = out_y[i_Phi, :]
- Psi = out_y[i_Psi, :]
- dTheta2 = out_dy[i_Theta2, :]
- dPhi = out_dy[i_Phi, :]
- dPsi = out_dy[i_Psi, :]
-
- # Then, compute Psi. Psi[:] assigns values rather than reassigning
- # reference.
- Psi[:] = -Phi - 12 * H_0**2 / c**2 / k**2 / a**2 * Omega_r * Theta2
- dPsi[:] = (-dPhi - 12 * H_0**2 / c**2 / k**2 * Omega_r / a**2
- * (dTheta2 - 2*Theta2))
-
- f_count = rhs_tight.f_count + rhs_full.f_count
- jacobian_count = rhs_tight.jacobian_count + rhs_full.jacobian_count
- return out_y, out_dy, f_count, jacobian_count
-
-#
-# Make a class for the solution, so that we can easily pickle it to file
-# (see part3.py)
-#
-class EinsteinBoltzmannSolution:
- def __init__(self, x_grid, k_grid, y, dy):
- self.x_grid = x_grid
- self.k_grid = k_grid
-
- # y has indices (parameter, k, x); unpack into seperate
- # variables for each parameter
- self.delta, self.delta_b, self.v, self.v_b, self.Phi = y[:i_Theta, :, :]
- self.Theta = y[i_Theta:i_Theta + lmax + 1, :, :]
- self.Psi = y[i_Theta + lmax + 1, :, :]
-
- self.ddelta, self.ddelta_b, self.dv, self.dv_b, self.dPhi = dy[:i_Theta, :, :]
- self.dTheta = dy[i_Theta:i_Theta + lmax + 1, :, :]
- self.dPsi = dy[i_Theta + lmax + 1, :, :]
-
-#
-# Set up grids in k and x and solve ODE for each k. This is a callable
-# function which is *not* called automatically due to the computing
-# time, so that plotting scripts etc. can use a pickled version.
-#
-def compute_einstein_boltzmann_grid(x_grid=None, k_grid=None, n_x=2500,
- use_jacobian=True):
- """
- Does all the computation for solving the Einstein-Boltzmann
- equations over a default grid over k. The result is an
- EinsteinBoltzmannSolution object.
- """
- if k_grid is None:
- # Grid in k is quadratic!
- k_grid = k_min + (k_max - k_min) * (np.arange(n_k) / (n_k-1))**2
- if x_grid is None:
- x_grid = np.linspace(x_init, 0, n_x)
-
- y = np.zeros((n_par, n_k, x_grid.shape[0]), np.double)
- dy = np.zeros_like(y)
- fctot = 0; jctot = 0
- for idx_k, k in enumerate(k_grid):
- a, b, fc, jc = solve_einstein_boltzmann(x_grid, k,
- out_y=y[:, idx_k, :],
- out_dy=dy[:, idx_k, :],
- use_jacobian=use_jacobian)
- fctot += fc
- jctot += jc
- return EinsteinBoltzmannSolution(x_grid, k_grid, y, dy)
-
-
View
0  einsteinboltzmann/gsl/__init__.py
No changes.
View
0  einsteinboltzmann/gsl/__init__.pyx
No changes.
View
12 einsteinboltzmann/gsl/error.pxd
@@ -1,12 +0,0 @@
-from stdio cimport FILE
-
-cdef extern from "gsl/gsl_errno.h":
- ctypedef void gsl_error_handler_t (char * reason, char * file,int line, int gsl_errno)
- ctypedef void gsl_stream_handler_t (char * label, char * file,int line, char * reason)
- void gsl_error (char * reason, char * file, int line, int gsl_errno)
- void gsl_stream_printf (char *label, char *file, int line, char *reason)
- char * gsl_strerror (int gsl_errno)
- gsl_error_handler_t * gsl_set_error_handler(gsl_error_handler_t * new_handler)
- gsl_error_handler_t * gsl_set_error_handler_off()
- gsl_stream_handler_t * gsl_set_stream_handler(gsl_stream_handler_t * new_handler)
- FILE * gsl_set_stream(FILE * new_stream)
View
20 einsteinboltzmann/gsl/error.pyx
@@ -1,20 +0,0 @@
-errors = []
-
-cdef extern from *:
- ctypedef char const_char "const char"
-
-cdef void _collecting_error_handler(const_char* reason,
- const_char* file,
- int line,
- int gsl_errno) with gil:
- print 'ERROR HANDLER CALLED'
- errors.append('%s:%d: GSL error %d: %s' % (
- (<bytes><char*>file).decode('ASCII'),
- line,
- gsl_errno,
- (<bytes><char*>reason).decode('ASCII')))
-
-
-def collect_errors():
- print 'starting to collect errors'
- gsl_set_error_handler(&_collecting_error_handler)
View
189 einsteinboltzmann/gsl/ode.pxd
@@ -1,189 +0,0 @@
-#
-# Cython pxd file for Gnu Scientific Library ODE solvers
-#
-
-
-cdef extern from *:
- cdef double const_double "const double"
-
-ctypedef (*gsl_odeiv_system_function)(double t, const_double *y, double *dydt, void *params)
-int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);
-
-
-cdef extern from "gsl/gsl_odeiv.h":
-
-
- ctypedef struct gsl_odeiv_system:
- gsl_odeiv_system_function function
- gsl_odeiv_system_jacobian jacobian
-
-typedef struct
-{
- int (* function) (double t, const double y[], double dydt[], void * params);
- int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);
- size_t dimension;
- void * params;
-}
-gsl_odeiv_system;
-
-#define GSL_ODEIV_FN_EVAL(S,t,y,f) (*((S)->function))(t,y,f,(S)->params)
-#define GSL_ODEIV_JA_EVAL(S,t,y,dfdy,dfdt) (*((S)->jacobian))(t,y,dfdy,dfdt,(S)->params)
-
-
-/* General stepper object.
- *
- * Opaque object for stepping an ODE system from t to t+h.
- * In general the object has some state which facilitates
- * iterating the stepping operation.
- */
-
-typedef struct
-{
- const char * name;
- int can_use_dydt_in;
- int gives_exact_dydt_out;
- void * (*alloc) (size_t dim);
- int (*apply) (void * state, size_t dim, double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv_system * dydt);
- int (*reset) (void * state, size_t dim);
- unsigned int (*order) (void * state);
- void (*free) (void * state);
-}
-gsl_odeiv_step_type;
-
-typedef struct {
- const gsl_odeiv_step_type * type;
- size_t dimension;
- void * state;
-}
-gsl_odeiv_step;
-
-
-/* Available stepper types.
- *
- * rk2 : embedded 2nd(3rd) Runge-Kutta
- * rk4 : 4th order (classical) Runge-Kutta
- * rkck : embedded 4th(5th) Runge-Kutta, Cash-Karp
- * rk8pd : embedded 8th(9th) Runge-Kutta, Prince-Dormand
- * rk2imp : implicit 2nd order Runge-Kutta at Gaussian points
- * rk4imp : implicit 4th order Runge-Kutta at Gaussian points
- * gear1 : M=1 implicit Gear method
- * gear2 : M=2 implicit Gear method
- */
-
-GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk2;
-GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk4;
-GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rkf45;
-GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rkck;
-GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk8pd;
-GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk2imp;
-GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk2simp;
-GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_rk4imp;
-GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_bsimp;
-GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_gear1;
-GSL_VAR const gsl_odeiv_step_type *gsl_odeiv_step_gear2;
-
-
-/* Constructor for specialized stepper objects.
- */
-gsl_odeiv_step * gsl_odeiv_step_alloc(const gsl_odeiv_step_type * T, size_t dim);
-int gsl_odeiv_step_reset(gsl_odeiv_step * s);
-void gsl_odeiv_step_free(gsl_odeiv_step * s);
-
-/* General stepper object methods.
- */
-const char * gsl_odeiv_step_name(const gsl_odeiv_step *);
-unsigned int gsl_odeiv_step_order(const gsl_odeiv_step * s);
-
-int gsl_odeiv_step_apply(gsl_odeiv_step *, double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv_system * dydt);
-
-/* General step size control object.
- *
- * The hadjust() method controls the adjustment of
- * step size given the result of a step and the error.
- * Valid hadjust() methods must return one of the codes below.
- *
- * The general data can be used by specializations
- * to store state and control their heuristics.
- */
-
-typedef struct
-{
- const char * name;
- void * (*alloc) (void);
- int (*init) (void * state, double eps_abs, double eps_rel, double a_y, double a_dydt);
- int (*hadjust) (void * state, size_t dim, unsigned int ord, const double y[], const double yerr[], const double yp[], double * h);
- void (*free) (void * state);
-}
-gsl_odeiv_control_type;
-
-typedef struct
-{
- const gsl_odeiv_control_type * type;
- void * state;
-}
-gsl_odeiv_control;
-
-/* Possible return values for an hadjust() evolution method.
- */
-#define GSL_ODEIV_HADJ_INC 1 /* step was increased */
-#define GSL_ODEIV_HADJ_NIL 0 /* step unchanged */
-#define GSL_ODEIV_HADJ_DEC (-1) /* step decreased */
-
-gsl_odeiv_control * gsl_odeiv_control_alloc(const gsl_odeiv_control_type * T);
-int gsl_odeiv_control_init(gsl_odeiv_control * c, double eps_abs, double eps_rel, double a_y, double a_dydt);
-void gsl_odeiv_control_free(gsl_odeiv_control * c);
-int gsl_odeiv_control_hadjust (gsl_odeiv_control * c, gsl_odeiv_step * s, const double y[], const double yerr[], const double dydt[], double * h);
-const char * gsl_odeiv_control_name(const gsl_odeiv_control * c);
-
-/* Available control object constructors.
- *
- * The standard control object is a four parameter heuristic
- * defined as follows:
- * D0 = eps_abs + eps_rel * (a_y |y| + a_dydt h |y'|)
- * D1 = |yerr|
- * q = consistency order of method (q=4 for 4(5) embedded RK)
- * S = safety factor (0.9 say)
- *
- * / (D0/D1)^(1/(q+1)) D0 >= D1
- * h_NEW = S h_OLD * |
- * \ (D0/D1)^(1/q) D0 < D1
- *
- * This encompasses all the standard error scaling methods.
- *
- * The y method is the standard method with a_y=1, a_dydt=0.
- * The yp method is the standard method with a_y=0, a_dydt=1.
- */
-
-gsl_odeiv_control * gsl_odeiv_control_standard_new(double eps_abs, double eps_rel, double a_y, double a_dydt);
-gsl_odeiv_control * gsl_odeiv_control_y_new(double eps_abs, double eps_rel);
-gsl_odeiv_control * gsl_odeiv_control_yp_new(double eps_abs, double eps_rel);
-
-/* This controller computes errors using different absolute errors for
- * each component
- *
- * D0 = eps_abs * scale_abs[i] + eps_rel * (a_y |y| + a_dydt h |y'|)
- */
-gsl_odeiv_control * gsl_odeiv_control_scaled_new(double eps_abs, double eps_rel, double a_y, double a_dydt, const double scale_abs[], size_t dim);
-
-/* General evolution object.
- */
-typedef struct {
- size_t dimension;
- double * y0;
- double * yerr;
- double * dydt_in;
- double * dydt_out;
- double last_step;
- unsigned long int count;
- unsigned long int failed_steps;
-}
-gsl_odeiv_evolve;
-
-/* Evolution object methods.
- */
-gsl_odeiv_evolve * gsl_odeiv_evolve_alloc(size_t dim);
-int gsl_odeiv_evolve_apply(gsl_odeiv_evolve *, gsl_odeiv_control * con, gsl_odeiv_step * step, const gsl_odeiv_system * dydt, double * t, double t1, double * h, double y[]);
-int gsl_odeiv_evolve_reset(gsl_odeiv_evolve *);
-void gsl_odeiv_evolve_free(gsl_odeiv_evolve *);
-
-
View
3  einsteinboltzmann/gsl/ode.pyx
@@ -1,3 +0,0 @@
-
-
-cdef
View
21 einsteinboltzmann/gsl/special.pyx
@@ -1,21 +0,0 @@
-cimport numpy as np
-import numpy as np
-
-cdef extern from "gsl/gsl_sf_bessel.h":
- int gsl_sf_bessel_jl_array (int lmax, double x, double* result_array)
- int gsl_sf_bessel_jl_steed_array (int lmax, double x, double* result_array)
-
-def bessel_jl_array(int lmax, double x, np.ndarray[double, mode='c'] out=None,
- algorithm='default'):
- if out is None:
- out = np.zeros((lmax + 1), np.double)
- cdef int ret
- if algorithm == 'steed':
- ret = gsl_sf_bessel_jl_steed_array(lmax, x, <double*>out.data)
- elif algorithm == 'default':
- ret = gsl_sf_bessel_jl_array(lmax, x, <double*>out.data)
- else:
- raise ValueError()
- if ret != 0:
- raise ValueError("GSL Error: %d" % ret)
- return out
View
127 einsteinboltzmann/gsl/spline.pxd
@@ -1,127 +0,0 @@
-cimport numpy as np
-
-cdef extern from *:
- cdef double const_double "const double"
-
-cdef extern from "gsl/gsl_interp.h":
-
- ctypedef struct gsl_interp_accel:
- size_t cache
- size_t miss_count
- size_t hit_count
-
- ctypedef struct gsl_interp_type:
- char* name
- unsigned int min_size
- void * (*alloc) (size_t size)
-# int (*init) (void *, const double xa[], const double ya[], size_t size);
- int (*eval) (void *, double* xa, double* ya, size_t size, double x,
- gsl_interp_accel *, double * y)
-## int (*eval_deriv) (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y_p);
-## int (*eval_deriv2) (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y_pp);
-## int (*eval_integ) (const void *, const double xa[], const double ya[], size_t size, gsl_interp_accel *, double a, double b, double * result);
-## void (*free) (void *);
-
-
- gsl_interp_type * gsl_interp_linear
- gsl_interp_type * gsl_interp_polynomial
- gsl_interp_type * gsl_interp_cspline
- gsl_interp_type * gsl_interp_cspline_periodic
- gsl_interp_type * gsl_interp_akima
- gsl_interp_type * gsl_interp_akima_periodic
-
- ctypedef struct gsl_interp:
- gsl_interp_type * type
- double xmin
- double xmax
- size_t size
- void * state
-
- gsl_interp_accel * gsl_interp_accel_alloc()
- int gsl_interp_accel_reset (gsl_interp_accel * a)
- void gsl_interp_accel_free(gsl_interp_accel * a)
-
- gsl_interp * gsl_interp_alloc(gsl_interp_type* T, int n)
- void gsl_interp_free(gsl_interp * interp)
-
- double gsl_interp_eval(gsl_interp * obj, double* xa, double* ya, double x,
- gsl_interp_accel * a)
-
- double gsl_interp_eval_deriv(gsl_interp * obj,
- double* xa, double* ya, double x,
- gsl_interp_accel * a)
-
- double gsl_interp_eval_deriv2(gsl_interp * obj,
- double* xa, double *ya, double x,
- gsl_interp_accel * a)
-
-
- int gsl_interp_init(gsl_interp * obj, double* xa, double* ya, size_t size)
- char * gsl_interp_name(gsl_interp * interp)
- unsigned int gsl_interp_min_size(gsl_interp * interp)
-
-
-cdef extern from "gsl/gsl_spline.h":
- ctypedef struct gsl_spline:
- gsl_interp * interp
- double * x
- double * y
- size_t size
-
- gsl_spline *gsl_spline_alloc(gsl_interp_type* T, size_t size)
-
- int gsl_spline_init(gsl_spline *spline, double *xa, double *ya, size_t size)
- char *gsl_spline_name(gsl_spline * spline)
- unsigned int gsl_spline_min_size(gsl_spline spline)
-
- int gsl_spline_eval_e(gsl_spline * spline, double x, gsl_interp_accel * a, double * y)
-
- double gsl_spline_eval(gsl_spline * spline, double x, gsl_interp_accel * a)
-
- int gsl_spline_eval_deriv_e(gsl_spline * spline,
- double x,
- gsl_interp_accel * a,
- double * y)
-
- double gsl_spline_eval_deriv(gsl_spline * spline,
- double x,
- gsl_interp_accel * a)
-
- int gsl_spline_eval_deriv2_e(gsl_spline * spline,
- double x,
- gsl_interp_accel * a,
- double * y)
-
- double gsl_spline_eval_deriv2(gsl_spline * spline,
- double x,
- gsl_interp_accel * a)
-
- int gsl_spline_eval_integ_e(gsl_spline * spline,
- double a, double b,
- gsl_interp_accel * acc,
- double * y)
-
- double gsl_spline_eval_integ(gsl_spline * spline,
- double a, double b,
- gsl_interp_accel * acc)
-
- void gsl_spline_free(gsl_spline * spline)
-
-cdef class Spline:
- cdef gsl_interp_accel *acc
- cdef gsl_interp *interp
- cdef size_t npoints
- cdef readonly object x, y
- cdef double *xa
- cdef double *ya
- cdef gsl_interp_type *interp_type
-
- cdef double evaluate_fast(self, double x)
- cdef double derivative_fast(self, double x)
- cdef double second_derivative_fast(self, double x)
-
-cdef class MultiSpline:
- cdef readonly np.ndarray x, y
- cdef gsl_interp_accel *acc
- cdef gsl_interp **interps
- cdef Py_ssize_t nsplines
View
258 einsteinboltzmann/gsl/spline.pyx
@@ -1,258 +0,0 @@
-
-cimport numpy as np
-import numpy as np
-cimport cython
-from libc.stdlib cimport malloc, free
-
-np.import_array()
-
-cdef class Spline:
-
- def __cinit__(self,
- np.ndarray[double, mode='c'] x,
- np.ndarray[double, mode='c'] y,
- algorithm='cspline'):
- self.acc = NULL
- self.interp = NULL
- if y.shape[0] != x.shape[0]:
- raise ValueError('x and y shapes do not match')
- self.npoints = x.shape[0]
- self.x = x
- self.y = y
- self.xa = <double*>x.data
- self.ya = <double*>y.data
- if algorithm == 'cspline':
- self.interp_type = gsl_interp_cspline
- else:
- raise NotImplementedError()
- self.interp = gsl_interp_alloc(self.interp_type, self.npoints)
- gsl_interp_init(self.interp, self.xa, self.ya, self.npoints)
- self.acc = gsl_interp_accel_alloc()
-
- def __dealloc__(self):
- if self.interp != NULL:
- gsl_interp_free(self.interp)
- if self.acc != NULL:
- gsl_interp_accel_free(self.acc)
-
- cdef double evaluate_fast(self, double x):
- return gsl_interp_eval(self.interp, self.xa, self.ya, x, self.acc)
-
- cdef double derivative_fast(self, double x):
- return gsl_interp_eval_deriv(self.interp, self.xa, self.ya, x, self.acc)
-
- cdef double second_derivative_fast(self, double x):
- return gsl_interp_eval_deriv2(self.interp, self.xa, self.ya, x, self.acc)
-
- def evaluate(self, x, out=None, hack=False):
- cdef Py_ssize_t i,j
- cdef gsl_interp_accel *acc = self.acc
- cdef gsl_interp *interp = self.interp
- cdef double* xa = self.xa
- cdef double* ya = self.ya
- cdef bint was_scalar = np.isscalar(x)
- cdef double* xbuf, *outbuf
- if was_scalar:
- x = [x]
- x = np.asarray(x, dtype=np.double)
- if out is None:
- out = np.empty_like(x)
- elif out.shape != x.shape:
- raise ValueError()
- cdef np.broadcast mit = np.broadcast(x, out)
- cdef int innerdim = np.PyArray_RemoveSmallest(mit)
- cdef Py_ssize_t xstride = x.strides[innerdim], outstride = out.strides[innerdim]
- cdef Py_ssize_t length = x.shape[innerdim]
- if xstride == sizeof(double) and outstride == sizeof(double):
- # Contiguous case
- while np.PyArray_MultiIter_NOTDONE(mit):
- xbuf = <double*>np.PyArray_MultiIter_DATA(mit, 0)
- outbuf = <double*>np.PyArray_MultiIter_DATA(mit, 1)
- for i in range(length):
- outbuf[i] = gsl_interp_eval(interp, xa, ya, xbuf[i], acc)
- np.PyArray_MultiIter_NEXT(mit)
- else:
- # Non-contiguous case
- raise NotImplementedError()
- if was_scalar:
- return out[0]
- else:
- return out
-
- def derivative(self, x, out=None):
- cdef Py_ssize_t i
- was_scalar = np.isscalar(x)
- if was_scalar:
- x = [x]
- cdef np.ndarray[double] xbuf = np.asarray(x, dtype=np.double)
- cdef np.ndarray[double] outbuf = out
- cdef gsl_interp_accel *acc = self.acc
- cdef gsl_interp *interp = self.interp
- cdef double* xa = self.xa
- cdef double* ya = self.ya
- if xbuf is None:
- raise ValueError()
- if out is None:
- outbuf = np.empty_like(xbuf)
- for i in range(xbuf.shape[0]):
- outbuf[i] = gsl_interp_eval_deriv(interp, xa, ya, xbuf[i], acc)
- if was_scalar:
- return outbuf[0]
- else:
- return outbuf
-
- def second_derivative(self, x, out=None):
- cdef Py_ssize_t i
- was_scalar = np.isscalar(x)
- if was_scalar:
- x = [x]
- cdef np.ndarray[double] xbuf = np.asarray(x, dtype=np.double)
- cdef np.ndarray[double] outbuf = out
- cdef gsl_interp_accel *acc = self.acc
- cdef gsl_interp *interp = self.interp
- cdef double* xa = self.xa
- cdef double* ya = self.ya
- if xbuf is None:
- raise ValueError()
- if out is None:
- outbuf = np.empty_like(xbuf)
- for i in range(xbuf.shape[0]):
- outbuf[i] = gsl_interp_eval_deriv2(interp, xa, ya, xbuf[i], acc)
- if was_scalar:
- return outbuf[0]
- else:
- return outbuf
-
- def __reduce__(self):
- version = 0
- return (_unpickle, (version, self.x, self.y, 'cspline'))
-
-cdef class MultiSpline:
-
- def __cinit__(self,
- np.ndarray[double] x,
- np.ndarray[double, ndim=2] y,
- algorithm='cspline'):
- self.acc = NULL
- self.interps = NULL
- cdef Py_ssize_t i, size
- size = x.shape[0]
- if y.shape[1] != x.shape[0]:
- raise ValueError()
- if y.strides[1] != sizeof(double):
- y = y.copy()
- if x.strides[0] != sizeof(double):
- x = x.copy()
- self.nsplines = y.shape[0]
- self.x = x
- self.y = y
-
- self.interps = <gsl_interp**>malloc(sizeof(gsl_interp*) * y.shape[0])
- self.acc = gsl_interp_accel_alloc()
- if algorithm == 'cspline':
- for i in range(y.shape[0]):
- self.interps[i] = gsl_interp_alloc(gsl_interp_cspline, size)
- gsl_interp_init(self.interps[i], <double*>x.data,
- <double*>(y.data + y.strides[0] * i), size)
- else:
- raise NotImplementedError()
-
- def __dealloc__(self):
- if self.acc == NULL:
- return
-
- cdef Py_ssize_t i
- for i in range(self.y.shape[0]):
- gsl_interp_free(self.interps[i])
- gsl_interp_accel_free(self.acc)
- free(self.interps)
-
- def evaluate(self, Py_ssize_t splineidx, x, out=None):
- cdef Py_ssize_t i,j
- cdef gsl_interp_accel *acc = self.acc
- if splineidx >= self.nsplines:
- raise ValueError('spline index out of range')
- cdef gsl_interp *interp = self.interps[splineidx]
- cdef bint was_scalar = np.isscalar(x)
- if was_scalar:
- x = [x]
- cdef np.ndarray[double] xbuf = np.asarray(x, dtype=np.double)
- if out is None:
- out = np.empty_like(x)
- elif out.shape != x.shape:
- raise ValueError()
- cdef np.ndarray[double] outbuf = out
- cdef double* xa = <double*>self.x.data
- cdef double* ya = <double*>(self.y.data + self.y.strides[0] * splineidx)
- for i in range(xbuf.shape[0]):
- outbuf[i] = gsl_interp_eval(interp, xa, ya, xbuf[i], acc)
- if was_scalar:
- return out[0]
- else:
- return out
-
-
-def _unpickle(version, x, y, algorithm):
- if version == 0:
- return Spline(x, y, algorithm)
- else:
- raise ValueError()
-
-def example1():
- """
- The example from the GSL documentation
- """
- cdef int i
- cdef double xi, yi
- cdef np.ndarray[double] x, y
- cdef np.ndarray[double] xhigh, yhigh
-
- ir = np.arange(10)
- x = (ir + 0.5 * np.sin(ir)).astype(np.double)
- y = (ir + np.cos(ir*ir)).astype(np.double)
-
- xhigh = np.linspace(x[0], x[9], 100)
- yhigh = np.zeros(100, dtype=np.double)
-
- cdef gsl_interp_accel *acc = gsl_interp_accel_alloc()
- cdef gsl_spline *spline = gsl_spline_alloc(gsl_interp_cspline, 10)
-
- try:
- gsl_spline_init(spline, <double*>x.data, <double*>y.data, 10)
-
- for i in range(xhigh.shape[0]):
- yhigh[i] = gsl_spline_eval(spline, xhigh[i], acc)
-
- finally:
- gsl_spline_free (spline);
- gsl_interp_accel_free (acc);
-
- import matplotlib.pyplot as plt
- plt.plot(x, y, 'r', xhigh, yhigh, 'g')
- plt.show()
-
-def example2():
- """
- The example from the GSL documentation using the higher-level interface
- """
- cdef int i
- cdef double xi, yi
- cdef np.ndarray[double] x, y
- cdef np.ndarray[double] xhigh, yhigh
-
- ir = np.arange(10)
- x = (ir + 0.5 * np.sin(ir)).astype(np.double)
- y = (ir + np.cos(ir*ir)).astype(np.double)
-
- xhigh = np.linspace(x[0], x[9], 100)
-
- spline = Spline(x, y)
- yhigh = spline.evaluate(xhigh)
-
-
- import matplotlib.pyplot as plt
- plt.plot(x, y, 'r', xhigh, yhigh, 'g')
- plt.show()
-
-
-
View
7 einsteinboltzmann/gsl/wscript
@@ -1,7 +0,0 @@
-def build(bld):
- bld.env['LIB'] = ['mkl_intel_lp64', 'mkl_sequential', 'mkl_core', 'pthread', 'gsl']
- for x in 'special spline error'.split():
- bld(target=x, source='%s.pyx' % x,
- features='c fc pyext cshlib',
- use='NUMPY'
- )
View
29 einsteinboltzmann/params.py
@@ -1,29 +0,0 @@
-from __future__ import division
-from numpy import pi
-
-# Units
-eV = 1.60217646e-19
-Mpc = 3.08568025e22
-
-# Cosmological parameters
-Omega_b = 0.046
-Omega_m = 0.224
-Omega_r = 8.3e-5 # 5.042e-5
-Omega_lambda = 1. - Omega_m - Omega_b - Omega_r
-T_0 = 2.725
-n_s = 1.
-A_s = 1.
-h0 = 0.7
-H_0 = h0 * 100. * 1.e3 / Mpc
-
-# General constants
-c = 2.99792458e8
-epsilon_0 = 13.605698 * eV
-m_e = 9.10938188e-31
-m_H = 1.673534e-27
-sigma_T = 6.652462e-29
-G_grav = 6.67258e-11
-rho_c = 3.*H_0**2 / (8.*pi*G_grav)
-alpha = 7.29735308e-3
-hbar = 1.05457148e-34
-k_b = 1.3806503e-23
View
3  einsteinboltzmann/rec_mod.pxd
@@ -1,3 +0,0 @@
-cdef double get_tau_fast(double x)
-cdef double get_dtau_fast(double x)
-cdef double get_ddtau_fast(double x)
View
230 einsteinboltzmann/rec_mod.pyx
@@ -1,230 +0,0 @@
-from __future__ import division
-import numpy as np
-from numpy import log, exp, sqrt, inf, pi
-from gsl.spline cimport Spline
-from scipy.integrate import ode
-import pprint
-
-from params import *
-from time_mod import get_H, get_H_p
-
-cimport libc.math as cm
-
-pp = pprint.PrettyPrinter(indent=4)
-np.seterr(all='raise') # fail hard on NaN and over/underflow
-
-#
-# Code for finding n_e
-#
-def solve_X_e_saha(x_grid, out, treshold=.99):
- """
- Solves for X_e using the Saha equation. Solving stops
- when the value falls below the given treshold.
-
- "out" is filled with the resulting values for each
- grid point.
-
- Returns the last grid index for which a value was computed.
- """
- a_grid = exp(x_grid)
- n_b = Omega_b * rho_c / m_H / a_grid**3
- T_b = T_0 / a_grid
-
- assert a_grid.shape == out.shape
- for i in range(a_grid.shape[0]):
- a = a_grid[i]
- # rhs of Saha equation
- rhs = (sqrt(m_e * T_b[i] * k_b / hbar**2 / 2 / pi)**3
- / n_b[i]
- * exp(-epsilon_0 / k_b / T_b[i]))
- # solve 2nd degree equation
- out[i] = (-rhs + sqrt(rhs**2 + 4*rhs)) / 2
- if out[i] < treshold:
- return i
- raise Exception("Never reached the %f treshold" % treshold)
-
-
-def peebles_X_e_derivative(x, y):
- """
- Right-hand side of the Peebles' ODE.
- """
- X_e = y[0]
-
- # Time-dependant "input" quantities
- a = exp(x)
- T_b = T_0 / a
- n_H = Omega_b * rho_c / m_H / a**3
- H = get_H(x)
-
- # Compute quantities in reverse order of listed in assignment
- phi_2 = 0.448 * log(epsilon_0 / T_b / k_b) # dimensionless
- alpha_2 = (64 * pi / sqrt(27 * pi)
- * (alpha * hbar / m_e)**2 / c
- * sqrt(epsilon_0 / T_b / k_b)
- * phi_2) # units of m^3/s
-
- # The following beta will become 0 numerically when T_b becomes
- # too small, and is only used in the final expression for the
- # derivative.
- beta = (alpha_2
- * (m_e * T_b * k_b / 2 / pi / hbar**2)**(3/2)
- * exp(-epsilon_0 / T_b / k_b)) # units of 1/s
-
- # In beta_2 we instead insert beta analytically,
- # to get a valid result even when beta underflows above
- beta_2 = (alpha_2
- * (m_e * T_b * k_b / 2 / pi / hbar**2)**(3/2)
- * exp(-(1/4) * epsilon_0 / T_b / k_b)) # units of 1/s
- n_1s = (1 - X_e) * n_H # units of 1/m^3
- Lambda_alpha = (H
- * (3 * epsilon_0 / c / hbar)**3
- / (8*pi)**2
- / n_1s) # units of 1/s
- Lambda_2s_to_1s = 8.227 # units of 1/s
-
- C_r = ((Lambda_2s_to_1s + Lambda_alpha)
- / (Lambda_2s_to_1s + Lambda_alpha + beta_2)) # dim.less
-
- # Finally, the main differential equation.
- diff = (beta * (1 - X_e) - n_H * alpha_2 * X_e**2)
- deriv = C_r / H * diff
- return [deriv]
-
-def solve_X_e_peebles(x_grid, X_e_init, out):
- """
- Solves for X_e using Peebles' equation. X_e_init is taken
- to be the initial condition for X_e at the point x_grid[0].
- """
- r = ode(peebles_X_e_derivative)
- r.set_integrator('vode', method='adams', with_jacobian=False,
- max_step=x_grid[1] - x_grid[0])
- r.set_initial_value(X_e_init, x_grid[0])
- out[0] = X_e_init
- for i in range(1, x_grid.shape[0]):
- r.integrate(x_grid[i])
- if not r.successful():
- raise Exception()
- out[i] = r.y
-
-def compute_X_e(x_grid):
- X_e = np.zeros_like(x_grid)
-
- # Solve using Saha to the treshold (idx gets the grid point this
- # happens at), then solve using Peebles' on the remainder of the
- # grid/values.
- switch_idx = solve_X_e_saha(x_grid, X_e)
- solve_X_e_peebles(x_grid[switch_idx:], X_e[switch_idx], X_e[switch_idx:])
-
- # Return gridded X_e
- return X_e
-
-def compute_n_e(x_grid, X_e_values):
- a_grid = exp(x_grid)
- n_H_values = Omega_b * rho_c / m_H / a_grid ** 3
- return X_e_values * n_H_values
-
-xstart = log(1e-10) # Start grid at a = 10^-10
-xstop = 0 # Stop grid at a = 1
-n = 1000 # Number of grid points between xstart and xstop
-x_grid = np.linspace(xstart, xstop, n)
-
-X_e_grid = x_grid
-X_e_values = compute_X_e(x_grid)
-
-n_e_values = compute_n_e(x_grid, X_e_values)
-log_n_e_spline = Spline(x_grid, log(compute_n_e(x_grid, X_e_values)))
-
-def get_n_e(x):
- return exp(log_n_e_spline.evaluate(x))
-
-#
-# Code for finding tau.
-#
-# tau is truncated to 0 between the second-to-last and last grid point.
-#
-def tau_derivative(x, y):
- tau = y[0]
- n_e = get_n_e(x)
- H_p = get_H_p(x)
- return -n_e * sigma_T * exp(x) * c / H_p
-
-def solve_tau(x_grid):
- # Integrate from the end of x_grid, i.e. reverse it
- reverse_x_grid = x_grid[::-1]
- out = np.zeros_like(x_grid)
- r = ode(tau_derivative)
- r.set_integrator('vode', method='adams', with_jacobian=False,
- max_step=x_grid[1] - x_grid[0])
- # There's no probability of scattering during a duration of zero length:
- r.set_initial_value(0, 0)
- out[0] = 0
- for i in range(1, reverse_x_grid.shape[0]):
- r.integrate(reverse_x_grid[i])
- if not r.successful():
- raise Exception()
- out[i] = r.y
- return out[::-1] # remeber to reverse values
-
-tau_values = solve_tau(x_grid)
-# When splining it, don't include the last (zero) point
-cdef double tau_spline_xstop = x_grid[-2]
-cdef Spline log_tau_spline = Spline(x_grid[:-1], log(tau_values[:-1]))
-
-# Python comment: get_tau is supposed to take both array and scalar
-# arguments (for easy plotting etc.), which means we must do without
-# if-tests (unless we want to code a loop). Just multiply instead...
-def get_tau(x):
- return exp(log_tau_spline.evaluate(x)) * (x <= tau_spline_xstop)
-
-def get_dtau(x):
- log_tau = log_tau_spline.evaluate(x)
- d_log_tau = log_tau_spline.derivative(x)
- return exp(log_tau) * d_log_tau * (x <= tau_spline_xstop)
-
-def get_ddtau(x):
- log_tau = log_tau_spline.evaluate(x)
- d_log_tau = log_tau_spline.derivative(x)
- dd_log_tau = log_tau_spline.second_derivative(x)
- return exp(log_tau) * (d_log_tau**2 + dd_log_tau) * (x <= tau_spline_xstop)
-
-# Fast, non-array version
-cdef double get_tau_fast(double x):
- return cm.exp(log_tau_spline.evaluate_single(x)) * (x <= tau_spline_xstop)
-
-cdef double get_dtau_fast(double x):
- cdef double log_tau, d_log_tau
- log_tau = log_tau_spline.evaluate_fast(x)
- d_log_tau = log_tau_spline.derivative_fast(x)
- return cm.exp(log_tau) * d_log_tau * (x <= tau_spline_xstop)
-
-cdef double get_ddtau_fast(double x):
- cdef double log_tau, d_log_tau, dd_log_tau
- log_tau = log_tau_spline.evaluate_fast(x)
- d_log_tau = log_tau_spline.derivative_fast(x)
- dd_log_tau = log_tau_spline.second_derivative_fast(x)
- return cm.exp(log_tau) * (d_log_tau**2 + dd_log_tau) * (x <= tau_spline_xstop)
-
-
-#
-# Code for finding g = g_tilde.
-#
-# g is set to 0 between the second-to-last and last grid point.
-#
-log_g_values = log(-get_dtau(x_grid[:-1])) - get_tau(x_grid[:-1])
-log_g_spline = Spline(x_grid[:-1], log_g_values)
-
-def get_g(x):
- return exp(log_g_spline.evaluate(x)) * (x <= tau_spline_xstop)
-
-def get_dg(x):
- log_g = log_g_spline.evaluate(x)
- d_log_g = log_g_spline.derivative(x)
- return exp(log_g) * d_log_g * (x <= tau_spline_xstop)
-
-def get_ddg(x):
- log_g = log_g_spline.evaluate(x)
- d_log_g = log_g_spline.derivative(x)
- dd_log_g = log_g_spline.second_derivative(x)
- return exp(log_g) * (d_log_g**2 + dd_log_g) * (x <= tau_spline_xstop)
-
-
View
55 einsteinboltzmann/run.py
@@ -1,55 +0,0 @@
-from __future__ import division
-
-from utils import print_time_taken
-import matplotlib.pyplot as plt
-import os
-import cPickle as pickle
-import numpy as np
-
-from evolution_mod import compute_einstein_boltzmann_grid, k_min, k_max
-
-# k_min = 7.56e-28
-# k_max = 7.56e-24
-
-k = 7.56e-25
-
-with print_time_taken('Solving...'):
- sol = compute_einstein_boltzmann_grid(n_x=1000, k_grid=np.array([k]))
-
-x_grid = sol.x_grid
-
-fig = plt.gcf()
-fig.clear()
-axs = [fig.add_subplot(2, 2, idx + 1) for idx in range(4)]
-axsiter = iter(axs)
-
-# delta plot
-idx_k = 0
-ax = axsiter.next()
-ax.semilogy(x_grid, abs(sol.delta[idx_k,:]), ':k', label=r'$|\delta|$')
-ax.semilogy(x_grid, abs(sol.delta_b[idx_k,:]), '-k', label=r'$|\delta_b|$')
-ax.set_ylim(1e-3, 1e5)
-ax.legend(loc='upper left')
-
-# v plot
-ax = axsiter.next()
-ax.semilogy(x_grid, abs(sol.v[idx_k,:]), ':k', label=r'$|v|$')
-ax.semilogy(x_grid, abs(sol.v[idx_k,:]), '-k', label=r'$|v_b|$')
-ax.set_ylim(1e-4, 1e2)
-ax.legend(loc='upper left')
-
-# Phi, Psi plot
-ax = axsiter.next()
-ax.plot(x_grid, sol.Phi[idx_k,:], ':k', label=r'$\Phi$')
-ax.plot(x_grid, -sol.Psi[idx_k,:], '-k', label=r'$-\Psi$')
-ax.set_ylim(-.05, 1.05)
-ax.legend(loc='lower left')
-
-# Theta plot
-ax = axsiter.next()
-ax.semilogy(x_grid, abs(sol.Theta[0, idx_k, :]), ':k', label=r'$|\Theta_0|$')
-ax.semilogy(x_grid, abs(sol.Theta[1, idx_k, :]), '-k', label=r'$|\Theta_1|$')
-ax.legend(loc='lower left')
-ax.set_ylim(1e-4, 10)
-
-plt.show()
View
4 einsteinboltzmann/time_mod.pxd
@@ -1,4 +0,0 @@
-cdef double get_eta_fast(double x)
-cdef double get_H_fast(double x)
-cdef double get_H_p_fast(double x)
-
View
132 einsteinboltzmann/time_mod.pyx
@@ -1,132 +0,0 @@
-from __future__ import division
-
-cdef double Omega_b, Omega_m, Omega_r, Omega_lambda
-cdef double H_0, c
-
-import numpy as np
-from numpy import log, inf
-from scipy.integrate import odeint
-from gsl.spline cimport Spline
-
-cimport libc.math as cm# cimport exp, sqrt
-
-from params import Omega_b, Omega_m, Omega_r, Omega_lambda, H_0, c
-
-def debug(x):
- print x
- import sys
- sys.stdout.flush()
-
-#
-# Define the grids
-#
-
-# Define two epochs, 1) during and 2) after recombination.
-
-n1 = 200 # Number of grid points during recombination
-n2 = 300 # Number of grid points after recombination
-z_start_rec = 1630.4 # Redshift of start of recombination
-z_end_rec = 614.2 # Redshift of end of recombination
-z_0 = 0 # Redshift today
-x_start_rec = -log(1 + z_start_rec) # x of start of recombination
-x_end_rec = -log(1 + z_end_rec) # x of end of recombination
-x_0 = 0 # x today
-
-n_eta = 1000 # Number of eta grid points (for spline)
-a_init = 1e-10 # Start value of a for eta evaluation
-x_eta1 = log(a_init) # Start value of x for eta evaluation
-x_eta2 = 0 # End value of x for eta evaluation
-
-x_grid = np.concatenate((
- np.linspace(np.log(1e-8), x_start_rec, 20, endpoint=False),
- np.linspace(x_start_rec, x_end_rec, n1, endpoint=False),
- np.linspace(x_end_rec, x_0, n2, endpoint=True)))
-
-x_eta_grid = np.linspace(x_eta1, x_eta2, n_eta)
-
-#
-# Utility
-#
-
-def redshift_to_x(rs):
- return -log(1 + rs)
-
-def x_to_redshift(x):
- return np.exp(-x) - 1
-
-
-#
-# Computation of H. Note: These also work elementwise on whole
-# arrays of "x" passed as argument.
-#
-
-def get_H(x):
- return H_0 * np.sqrt(
- (Omega_b + Omega_m) * np.exp(-3 * x)
- + Omega_r * np.exp(-4 * x)
- + Omega_lambda)
-
-def get_H_p(x):
- return np.exp(x) * get_H(x)
-
-def get_dH_p(x):
- dHdx = (-H_0 / 2
- * (3 * (Omega_b + Omega_m) * np.exp(-3 * x)
- + 4 * Omega_r * np.exp(-4 * x))
- / np.sqrt((Omega_b + Omega_m) * np.exp(-3 * x) +
- Omega_r * np.exp(-4*x) + Omega_lambda))
- return np.exp(x) * (get_H(x) + dHdx)
-
-# Fast, compiled versions which don't work with arrays
-cdef double get_H_fast(double x):
- return H_0 * cm.sqrt(
- (Omega_b + Omega_m) * cm.exp(-3 * x)
- + Omega_r * cm.exp(-4 * x)
- + Omega_lambda)
-
-cdef double get_H_p_fast(double x):
- return cm.exp(x) * get_H_fast(x)
-
-#
-# Computation and retrieval of eta. get_eta is the front-end
-# for user code.
-#
-
-def eta_derivative(eta, a):
- # Simplified the expression up-front in order to avoid
- # problems when a is close to 0
- # (since exp(log(0)) != 0...)
- return c / H_0 / np.sqrt(
- (Omega_b + Omega_m) * a
- + Omega_r
- + Omega_lambda * a**4)
-
-def compute_eta(x_grid):
- # Compute it on a grid on a = exp(x).
- # Add an initial 0 to the grid for the boundary condition
- a_grid = np.concatenate(([0], np.exp(x_grid)))
- eta_at_0 = 0 # Light has not travelled before any time has passed
-
- # Solve ODE
- eta, info = odeint(func=eta_derivative,
- y0=eta_at_0,
- t=a_grid,
- rtol=1e-12,
- atol=1e-20,
- full_output=True)
- # trim off first point (t=0) on return (1:)
- # also, odeint returns an n-by-1 array, so index to return a 1-dim array
- return eta[1:,0]
-
-def compute_eta_spline(x_grid):
- eta_vals = compute_eta(x_grid)
- result = Spline(x_grid, eta_vals, algorithm='cspline')
- return result
-
-cdef Spline eta_spline = compute_eta_spline(x_eta_grid)
-def get_eta(x):
- r = eta_spline.evaluate(x)
- return r
-
-cdef double get_eta_fast(double x):
- return eta_spline.evaluate_fast(x)
View
43 einsteinboltzmann/tools/cython.py
@@ -1,43 +0,0 @@
-# Cython tool for waf
-
-from waflib.Configure import conf
-from waflib import TaskGen
-
-TaskGen.declare_chain(
- name = "cython",
- rule = "${CYTHON} ${CYTHONFLAGS} ${CPPPATH_ST:INCPATHS} ${SRC} -o ${TGT}",
- ext_in = ['.pyx'],
- ext_out = ['.c'],
- reentrant = True,
- )
-
-@conf
-def check_cython_version(conf, minver):
- conf.start_msg("Checking cython version")
- minver = tuple(minver)
- import re
- version_re = re.compile(r'cython\s*version\s*(?P<major>\d*)\.(?P<minor>\d*)(?:\.(?P<micro>\d*))?', re.I).search
- cmd = conf.cmd_to_list(conf.env['CYTHON'])
- cmd = cmd + ['--version']
- from waflib.Tools import fc_config
- stdout, stderr = fc_config.getoutput(conf, cmd)
- if stdout:
- match = version_re(stdout)
- else:
- match = version_re(stderr)
- if not match:
- conf.fatal("cannot determine the Cython version")
- cy_ver = [match.group('major'), match.group('minor')]
- if match.group('micro'):
- cy_ver.append(match.group('micro'))
- else:
- cy_ver.append('0')
- cy_ver = tuple([int(x) for x in cy_ver])
- if cy_ver < minver:
- conf.end_msg(False)
- conf.fatal("cython version %s < %s" % (cy_ver, minver))
- conf.end_msg(str(cy_ver))
-
-def configure(conf):
- conf.find_program('cython', var='CYTHON')
-
View
181 einsteinboltzmann/tools/fwrapktp.py
@@ -1,181 +0,0 @@
-# fwrap kind-type probing tool for waf
-#
-# TODO: Get rid of fwrap dependency; bundle necessarry parts of
-# gen_config.py as a tool.
-
-from waflib import Logs, Build, Utils
-#from waflib.Task import Task
-from waflib import TaskGen, Task
-from waflib.TaskGen import after_method, before_method, feature, taskgen_method, extension
-import os
-
-gc = None
-
-@TaskGen.feature('fwrapktp')
-@TaskGen.after('process_source')
-@TaskGen.before('apply_link')
-def generate_fwrap_ktp_headers(self):
- node = self.path.find_or_declare(getattr(self, 'typemap', generate_fwrap_ktp.typemap_in))
-
- typemap_f90 = self.path.find_or_declare(generate_fwrap_ktp.typemap_f90)
- typemap_h = self.path.find_or_declare(generate_fwrap_ktp.typemap_h)
- typemap_pxd = self.path.find_or_declare(generate_fwrap_ktp.typemap_pxd)
- typemap_pxi = self.path.find_or_declare(generate_fwrap_ktp.typemap_pxi)
-
- inputs = [node]
- outputs = [typemap_f90, typemap_h, typemap_pxd, typemap_pxi]
-
- gen_task = self.create_task("generate_fwrap_ktp", inputs, outputs)
-
- fc_task = self.create_compiled_task('fc', typemap_f90)
- self.includes_nodes = [] # *shrug*
- fc_task.nomod = True # the fortran files won't compile unless all the .mod files are set, ick
-
-
-class generate_fwrap_ktp(Task.Task):
- before = ['c', 'fc', 'cython']
-
- typemap_in = 'fwrap_type_specs.in'
- typemap_f90 = 'fwrap_ktp_mod.f90'
- typemap_h = 'fwrap_ktp_header.h'
- typemap_pxd = 'fwrap_ktp.pxd'
- typemap_pxi = 'fwrap_ktp.pxi'
-
- def run(self):
- obld = self.generator.bld
- bld = Build.BuildContext(top_dir=obld.srcnode.abspath(), out_dir=obld.bldnode.abspath())
- bld.init_dirs()
- bld.in_msg = 1 # disable all that comes from bld.msg(..), bld.start_msg(..) and bld.end_msg(...)
- bld.env = self.env.derive()
- ktp_in, = self.inputs
- bld.logger = Logs.make_logger(bld.bldnode.abspath() + os.sep + ktp_in.name + '.log', 'build')
-
- ctps = gc.read_type_spec(ktp_in.abspath())
- find_types(bld, ctps)
-
- typemap_f90, typemap_h, typemap_pxd, typemap_pxi = self.outputs
-
- def genfile(generator, node, *args):
- generator(ctps, node, *args)
-
- genfile(gc.write_f_mod, typemap_f90)
- genfile(gc.write_header, typemap_h)
- genfile(gc.write_pxd, typemap_pxd, self.typemap_h)
- genfile(gc.write_pxi, typemap_pxi)
-
-
-def configure(conf):
- # TODO: Get rid of this dependency
- conf.check_python_module('fwrap')
- global gc
- from fwrap import gen_config as gc
-
-def process_typemap(self, **kw):
- return self(rule=_process_typemap_builder, **kw)
-Build.BuildContext.process_typemap = process_typemap
-
-
-def find_types(bld, ctps):
- for ctp in ctps:
- fc_type = None
- if ctp.lang == 'fortran':
- fc_type = find_fc_type(bld, ctp.basetype,
- ctp.odecl, ctp.possible_modules)
- elif ctp.lang == 'c':
- fc_type = find_c_type(bld, ctp)
- if not fc_type:
- raise bld.errors.WafError(
- "unable to find C type for type %s" % ctp.odecl)
- ctp.fc_type = fc_type
-
-
-fc_type_memo = {}
-def find_fc_type(bld, basetype, decl, possible_modules):
- res = fc_type_memo.get((basetype, decl), None)
- if res is not None:
- return res
-
- MODULES = '\n'.join('use %s' % x for x in possible_modules)
-
- if basetype == 'logical':
- basetype = 'integer'
- decl = decl.replace('logical', 'integer')
-
- fsrc_tmpl = '''\
-subroutine outer(a)
- use, intrinsic :: iso_c_binding
- %(MODULES)s
- implicit none
- %(TEST_DECL)s, intent(inout) :: a
- interface
- subroutine inner(a)
- use, intrinsic :: iso_c_binding
- %(MODULES)s
- implicit none
- %(TYPE_DECL)s, intent(inout) :: a
- end subroutine inner
- end interface
- call inner(a)
-end subroutine outer
-'''
- for ctype in gc.type_dict[basetype]:
- test_decl = '%s(kind=%s)' % (basetype, ctype)
- fsrc = fsrc_tmpl % {'TYPE_DECL' : decl,
- 'TEST_DECL' : test_decl,
- 'MODULES' : MODULES}
- try:
- bld.check_cc(
- fragment=fsrc,
- compile_filename='test.f90',
- features='fc',
- includes = bld.bldnode.abspath())
- except bld.errors.ConfigurationError:
- pass
- else:
- res = ctype
- break
- else:
- res = ''
- fc_type_memo[basetype, decl] = res
- return res
-
-def find_c_type(bld, ctp):
- if ctp.lang != 'c':
- raise ValueError("wrong language, given %s, expected 'c'" % ctp.lang)
- if ctp.basetype != 'integer':
- raise ValueError(
- "only integer basetype supported for C type discovery.")
-
- tmpl = r'''
-#include "Python.h"
-#include "numpy/arrayobject.h"
-
-typedef %(type)s npy_check_sizeof_type;
-int main(int argc, char **argv)
-{
- static int test_array [1 - 2 * !(((long) (sizeof (npy_check_sizeof_type))) == sizeof(%(ctype)s))];
- test_array [0] = 0
-
- ;
- return 0;
-}
-'''
- ctypes = ('signed char', 'short int',
- 'int', 'long int', 'long long int')
- for ctype in ctypes:
- cfrag = tmpl % {'type' : ctp.odecl, 'ctype' : ctype}
- try:
- bld.check_cc(
- fragment=cfrag,
- features = 'c',
- compile_filename='test.c',
- use='NUMPY PYEXT')
- except bld.errors.ConfigurationError:
- pass
- else:
- res = ctype
- break
- else:
- res = ''
- return gc.c2f[res]
-
View
64 einsteinboltzmann/tools/fwraptool.py
@@ -1,64 +0,0 @@
-from waflib.Configure import conf
-from waflib import TaskGen
-from waflib import Task
-
-@TaskGen.feature('fwrap')
-@TaskGen.before('process_source')
-def fwrap_fortran_sources(self):
- input_sources = list(self.source)
- f90_wrapped = False
- for node in input_sources:
- node = self.bld.srcnode.find_or_declare(node)
- if node.suffix() in ('.f90', '.f'):
- if f90_wrapped:
- raise ValueError('fwrap feature only works with a single Fortran source')
- f90_wrapped = True
-
- if node.suffix() == '.f':
- f77binding = True
- self.env['FWRAP_F77BINDING'] = '--f77binding'
- suffix = '.f'
- else:
- f77binding = False
- suffix = '.f90'
-
- # Automatically set target name if not provided
- name = node.name[:-len(node.suffix())]
- if not getattr(self, 'target', None):
- self.target = name
-
- pyx = node.change_ext('.pyx')
- typemap_h = node.parent.find_or_declare('fwrap_ktp_header.h')
- typemap_pxd = node.parent.find_or_declare('fwrap_ktp.pxd')
- typemap_pxi = node.parent.find_or_declare('fwrap_ktp.pxi')
-
- if not f77binding:
- fc_f = node.change_ext('_fc' + suffix)
- ktp = node.parent.find_or_declare('fwrap_type_specs.in')
- typemap_f = node.parent.find_or_declare('fwrap_ktp_mod' + suffix)
-
- self.create_task('fwrap', src=[node], tgt=[pyx, fc_f, ktp])
- self.create_task('generate_fwrap_ktp',
- src=[ktp],