Commit f987fea0 authored by Benjamin Audren's avatar Benjamin Audren
Browse files

Fixed minor bugs and improved functionalities of CPU.py

Warning: this is a cherry picked version of the devel branch (see the cherry
picked commit).

* 2.3.2:
  Cherry picked version of the devel branch (2.3.2)
  added the derived param r_0002
  CPU.py now computes the ratios
parents 9c6673f4 4d2e4f05
......@@ -21,11 +21,18 @@ python CPU.py output/lcdm_z2_pk.dat output/lncdm_z2_pk.dat -r
"""
import numpy as np
from numpy import ma
MaskedArray = ma.MaskedArray
from scipy.interpolate import splrep, splev
import os
import matplotlib.pyplot as plt
import sys
import argparse
import itertools
from math import floor
from matplotlib import scale as mscale
from matplotlib.transforms import Transform
from matplotlib.ticker import FixedLocator
START_LINE = {}
START_LINE['error'] = [r' /|\ ',
......@@ -57,22 +64,151 @@ def create_parser():
parser.add_argument('-s', '--selection', dest='selection',
nargs='+',
help='specify the fields you want to plot.')
parser.add_argument('--scale', choices=['lin', 'loglog', 'loglin'],
type=str,
parser.add_argument('--scale', type=str,
choices=['lin', 'loglog', 'loglin', 'george'],
help='Specify the scale to use for the plot')
parser.add_argument('--xlim', dest='xlim', nargs='+', type=float,
default=[], help='Specify the x range')
parser.add_argument(
'-p, --print',
dest='printfile', action='store_true', default=False,
help='print the graph directly in a .png file')
help='print the graph directly in a .pdf file')
parser.add_argument(
'-r, --repeat',
dest='repeat', action='store_true', default=False,
help='repeat the step for all redshifts with same base name')
return parser
# Helper code from cosmo_mini_toolbox, by Jesus Torrado, available fully at
# https://github.com/JesusTorrado/cosmo_mini_toolbox, to use the log then
# linear scale for the multipole axis when plotting Cl.
nonpos = "mask"
change = 50.0
factor = 500.
def plot_CLASS_output(files, selection, ratio=False, output_name='',
extension='', x_variable='', scale='lin'):
def _mask_nonpos(a):
"""
Return a Numpy masked array where all non-positive 1 are
masked. If there are no non-positive, the original array
is returned.
"""
mask = a <= 0.0
if mask.any():
return ma.MaskedArray(a, mask=mask)
return a
def _clip_smaller_than_one(a):
a[a <= 0.0] = 1e-300
return a
class PlanckScale(mscale.ScaleBase):
"""
Scale used by the Planck collaboration to plot Temperature power spectra:
base-10 logarithmic up to l=50, and linear from there on.
Care is taken so non-positive values are not plotted.
"""
name = 'planck'
def __init__(self, axis, **kwargs):
pass
def set_default_locators_and_formatters(self, axis):
axis.set_major_locator(
FixedLocator(
np.concatenate((np.array([2, 10, change]),
np.arange(500, 2500, 500)))))
axis.set_minor_locator(
FixedLocator(
np.concatenate((np.arange(2, 10),
np.arange(10, 50, 10),
np.arange(floor(change/100), 2500, 100)))))
def get_transform(self):
"""
Return a :class:`~matplotlib.transforms.Transform` instance
appropriate for the given logarithm base.
"""
return self.PlanckTransform(nonpos)
def limit_range_for_scale(self, vmin, vmax, minpos):
"""
Limit the domain to positive values.
"""
return (vmin <= 0.0 and minpos or vmin,
vmax <= 0.0 and minpos or vmax)
class PlanckTransform(Transform):
input_dims = 1
output_dims = 1
is_separable = True
has_inverse = True
def __init__(self, nonpos):
Transform.__init__(self)
if nonpos == 'mask':
self._handle_nonpos = _mask_nonpos
else:
self._handle_nonpos = _clip_nonpos
def transform_non_affine(self, a):
lower = a[np.where(a<=change)]
greater = a[np.where(a> change)]
if lower.size:
lower = self._handle_nonpos(lower * 10.0)/10.0
if isinstance(lower, MaskedArray):
lower = ma.log10(lower)
else:
lower = np.log10(lower)
lower = factor*lower
if greater.size:
greater = (factor*np.log10(change) + (greater-change))
# Only low
if not(greater.size):
return lower
# Only high
if not(lower.size):
return greater
return np.concatenate((lower, greater))
def inverted(self):
return PlanckScale.InvertedPlanckTransform()
class InvertedPlanckTransform(Transform):
input_dims = 1
output_dims = 1
is_separable = True
has_inverse = True
def transform_non_affine(self, a):
lower = a[np.where(a<=factor*np.log10(change))]
greater = a[np.where(a> factor*np.log10(change))]
if lower.size:
if isinstance(lower, MaskedArray):
lower = ma.power(10.0, lower/float(factor))
else:
lower = np.power(10.0, lower/float(factor))
if greater.size:
greater = (greater + change - factor*np.log10(change))
# Only low
if not(greater.size):
return lower
# Only high
if not(lower.size):
return greater
return np.concatenate((lower, greater))
def inverted(self):
return PlanckTransform()
# Finished. Register the scale!
mscale.register_scale(PlanckScale)
def plot_CLASS_output(files, selection, ratio=False, printing=False,
output_name='', extension='', x_variable='',
scale='lin', xlim=[]):
"""
Load the data to numpy arrays, write a Python script and plot them.
......@@ -105,6 +241,7 @@ def plot_CLASS_output(files, selection, ratio=False, output_name='',
# Create the python script, and initialise it
python_script_path = files[0]+'.py'
pdf_path = files[0]+'.pdf'
text = '''
import matplotlib.pyplot as plt
import numpy as np\n'''
......@@ -149,6 +286,7 @@ for data_file in files:
# if ratio is not set, then simply plot them all
if not ratio:
loc = ''
text += 'for curve in data:\n'
for idx, curve in enumerate(data):
_, curve_names, _ = extract_headers(files[idx])
......@@ -161,15 +299,52 @@ for data_file in files:
elif scale == 'loglog':
text += 'loglog(curve[:, 0], curve[:, %i])\n' % index
ax.loglog(curve[:, 0], curve[:, index])
elif scale == 'loglin':
text += 'semilogx(curve[:, 0], curve[:, %i])\n' % index
ax.semilogx(curve[:, 0], curve[:, index])
elif scale == 'george':
text += 'plot(curve[:, 0], curve[:, %i])\n' % index
ax.plot(curve[:, 0], curve[:, index])
ax.set_xscale('planck')
loc = 'upper right'
if not loc:
loc = 'lower right'
ax.legend([root+': '+elem for (root, elem) in
itertools.product(roots, selection)], loc='lower right')
#ax.legend([
itertools.product(roots, selection)], loc=loc)
else:
ref = data[0]
#for index in range(1, len(data)):
#current = data[index]
_, ref_curve_names, _ = extract_headers(files[0])
for idx in range(1, len(data)):
current = data[idx]
_, curve_names, _ = extract_headers(files[idx])
for selec in selection:
# Do the interpolation
axis = ref[:, 0]
reference = ref[:, ref_curve_names.index(selec)]
interpolated = splrep(current[:, 0],
current[:, curve_names.index(selec)])
if scale == 'lin':
ax.plot(axis, splev(ref[:, 0], interpolated)/reference-1)
elif scale == 'loglin':
ax.semilogx(axis, splev(ref[:, 0],
interpolated)/reference-1)
elif scale == 'loglog':
raise InputError(
"loglog plot is not available for ratios")
#if np.allclose(current[0], ref[0]):
#ax.plot(current[0], current[colnum]/ref[colnum])
if 'TT' in curve_names:
ax.set_xlabel('$\ell$', fontsize=20)
elif 'P' in curve_names:
ax.set_xlabel('$k$ [$h$/Mpc]', fontsize=20)
if xlim:
if len(xlim) > 1:
ax.set_xlim(xlim)
else:
ax.set_xlim(xlim[0])
ax.set_ylim()
text += 'plt.show()\n'
plt.show()
......@@ -178,6 +353,10 @@ for data_file in files:
with open(python_script_path, 'w') as python_script:
python_script.write(text)
# If the use wants to print the figure to a file
if printing:
fig.savefig(pdf_path)
class FormatError(Exception):
"""Format not recognised"""
......@@ -247,6 +426,8 @@ def process_long_names(long_names):
elif name.find('[') != -1:
names[index] = name[:name.index('[')]
# Finally, remove any extra spacing
names = [''.join(elem.split()) for elem in names]
return names, tex_names
......@@ -282,18 +463,14 @@ def main():
parser.print_usage()
return
# Ratio is not implemented yet, catch it
if args.ratio:
raise InputError(
"Sorry, this is not working yet")
# if the first file name contains cl or pk, infer the type of desired
# spectrum
if not args.selection:
if args.files[0].rfind('cl') != -1:
selection = 'TT'
selection = ['TT']
scale = 'loglog'
elif args.files[0].rfind('pk') != -1:
selection = 'P'
selection = ['P']
scale = 'loglog'
else:
raise TypeError(
......@@ -307,6 +484,8 @@ def main():
else:
args.scale = 'lin'
# Remove extra spacing in the selection list
args.selection = [''.join(elem.split()) for elem in args.selection]
# If ratio is asked, but only one file was passed in argument, politely
# complain
if args.ratio:
......@@ -317,8 +496,9 @@ def main():
# actual plotting. By default, a simple superposition of the graph is
# performed. If asked to be divided, the ratio is shown - whether a need
# for interpolation arises or not.
plot_CLASS_output(args.files, args.selection,
ratio=args.ratio, scale=args.scale)
plot_CLASS_output(args.files, args.selection, ratio=args.ratio,
printing=args.printfile, scale=args.scale,
xlim=args.xlim)
if __name__ == '__main__':
sys.exit(main())
/** @file bessel.h Documented includes for Bessel module */
#ifndef __BESSEL__
#define __BESSEL__
#include "common.h"
#include "arrays.h"
/**
* Structure containing everything about spherical bessel functions
* that other modules need to know.
*
* Once initialized by bessel_init(), contains table of
* spherical Bessel functions \f$ j_l(x) \f$).
*/
struct bessels {
/** @name - input parameters initialized by user in input module
(all other quantitites are computed in this module, given these
parameters and the content of the 'precision' structure) */
//@{
int l_max; /**< last l value */
double x_max; /**< maximum value of x (always multiple of x-step) */
double x_step; /**< step dx for sampling Bessel functions */
short bessel_always_recompute; /**< if set to true, Bessels are never read from / written in files */
short use_pbs; /**< if _TRUE_, try to get bessel function from this module, if _FALSE_, from new hypershperical module. For non-flat models this parameter is forced by input module to be _FALSE_ */
int get_HIS_from_shared_memory; /**< flag specifying if class should try to get HIS from shared memory */
//@}
/** @name - parameters defining uniquely the exact content of the Bessel table
(hence when reading a file, will compare these values with needed value
in order to take the decision to recompute or not) */
//@{
int l_size; /**< number of multipole values */
int * l; /**< list of multipole values, l[index_l] */
double j_cut; /**< value of \f$ j_l \f$ below which it is approximated by zero (in the region x << l) */
int has_dj; /**< set to true means j_l'(x) also need to be stored */
//@}
/** @name - Bessel table, and arrays necessary for reading it */
//@{
int * x_size; /**< x_size[index_l] is the number of x values for l[index_l]; hence *x_min[index_l]+x_step*(x_size[index_l]-1) = x_max */
int x_size_max;
double ** buffer; /**< buffer[index_l] is a pointer towards a memory zone containing x_min, all j_l(x) and all j_l''(x) for each value of l */
double ** x_min; /**< x_min[index_l] is a pointer towards the minimum value of x for l[index_l], given j_cut; always a multiple of x_step */
double ** j; /* (j[index_l])[index_x] is \f$ j_l(x) \f$ for l[index_l] and x=x_min[index_l]+x_step*index_x */
double ** ddj; /* (ddj[index_l])[index_x] \f$ j_l''(x) \f$ for l[index_l] and x=x_min[index_l]+x_step*index_x (in view of spline interpolation) */
double ** dj; /* (dj[index_l])[index_x] is \f$ j_l'(x) \f$ for l[index_l] and x=x_min[index_l]+x_step*index_x */
double ** dddj; /* (dddj[index_l])[index_x] \f$ j_l'''(x) \f$ for l[index_l] and x=x_min[index_l]+x_step*index_x (in view of spline interpolation) */
//@}
/** @name - technical parameters */
//@{
short bessels_verbose; /**< flag regulating the amount of information sent to standard output (none if set to zero) */
ErrorMsg error_message; /**< zone for writing error messages */
//@}
};
/*************************************************************************************************************/
/*
* Boilerplate for C++
*/
#ifdef __cplusplus
extern "C" {
#endif
int bessel_at_x(
struct bessels * pbs,
double x,
int l,
double * j
);
int bessel_init(
struct precision * ppr,
struct bessels * pbs
);
int bessel_free(
struct bessels * pbs
);
int bessel_get_l_list(
struct precision * ppr,
struct bessels * pbs
);
int bessel_j_for_l(
struct precision * ppr,
struct bessels * pbs,
int index_l
);
int bessel_j(
struct bessels * pbs,
int l,
double x,
double * jl
);
#ifdef __cplusplus
}
#endif
/**
* @name Constants used for computing Bessel function
*/
//@{
#define _GAMMA1_ 2.6789385347077476336556
#define _GAMMA2_ 1.3541179394264004169452
//@}
#endif
......@@ -15,7 +15,7 @@
#ifndef __COMMON__
#define __COMMON__
#define _VERSION_ "v2.3.1"
#define _VERSION_ "v2.3.2"
#define _TRUE_ 1 /**< integer associated to true statement */
#define _FALSE_ 0 /**< integer associated to false statement */
......
/**
* definitions for module thermodynamics.c
*/
#ifndef __TOOLS_ARRAYS__
#define __TOOLS_ARRAYS__
#include "common.h"
#define _SPLINE_NATURAL_ 0 /**< natural spline: ddy0=ddyn=0 */
#define _SPLINE_EST_DERIV_ 1 /**< spline with estimation of first derivative on both edges */
/**
* Boilerplate for C++
*/
#ifdef __cplusplus
extern "C" {
#endif
int array_derive(
double * array,
int n_columns,
int n_lines,
int index_x, /** from 0 to (n_columns-1) */
int index_y,
int index_dydx,
ErrorMsg errmsg);
int array_derive_spline(
double * x_array,
int n_lines,
double * array,
double * array_splined,
int n_columns,
int index_y,
int index_dydx,
ErrorMsg errmsg);
int array_derive_spline_table_line_to_line(
double * x_array,
int n_lines,
double * array,
int n_columns,
int index_y,
int index_ddy,
int index_dy,
ErrorMsg errmsg);
int array_derive1_order2_table_line_to_line(
double * x_array,
int n_lines,
double * array,
int n_columns,
int index_y,
int index_dy,
ErrorMsg errmsg);
int array_derive2_order2_table_line_to_line(
double * x_array,
int n_lines,
double * array,
int n_columns,
int index_y,
int index_dy,
int index_ddy,
ErrorMsg errmsg);
int array_derive_two(
double * array,
int n_columns,
int n_lines,
int index_x, /** from 0 to (n_columns-1) */
int index_y,
int index_dydx,
int index_ddydxdx,
ErrorMsg errmsg);
int array_spline(
double * array,
int n_columns,
int n_lines,
int index_x, /** from 0 to (n_columns-1) */
int index_y,
int index_ddydx2,
short spline_mode,
ErrorMsg errmsg);
int array_spline_table_line_to_line(
double * x, /* vector of size x_size */
int x_size,
double * array,
int n_columns,
int index_y,
int index_ddydx2,
short spline_mode,
ErrorMsg errmsg);
int array_spline_table_columns(
double * x,
int x_size,
double * y_array,
int y_size,
double * ddy_array,
short spline_mode,
ErrorMsg errmsg);
int array_spline_table_lines(
double * x,
int x_size,
double * y_array,
int y_size,
double * ddy_array,
short spline_mode,
ErrorMsg errmsg
);
int array_logspline_table_lines(
double * x,
int x_size,
double * y_array,
int y_size,
double * ddlny_array,
short spline_mode,
ErrorMsg errmsg
);
int array_spline_table_one_column(
double * x, /* vector of size x_size */
int x_size,
double * y_array, /* array of size x_size*y_size with elements
y_array[index_y*x_size+index_x] */
int y_size,
int index_y,
double * ddy_array, /* array of size x_size*y_size */
short spline_mode,
ErrorMsg errmsg
);
int array_logspline_table_one_column(
double * x, /* vector of size x_size */
int x_size,
int x_stop,
double * y_array, /* array of size x_size*y_size with elements
y_array[index_y*x_size+index_x] */
int y_size,
int index_y,
double * ddlogy_array, /* array of size x_size*y_size */
short spline_mode,
ErrorMsg errmsg
);
int array_integrate_all_spline(
double * array,
int n_columns,
int n_lines,
int index_x,
int index_y,
int index_ddy,
double * result,
ErrorMsg errmsg
);
int array_integrate(
double * array,
int n_columns,
int n_lines,
int index_x, /** from 0 to (n_columns-1) */
int index_y,
int index_int_y_dx,
ErrorMsg errmsg);