Commit 7aee1fe9 authored by Julien Lesgourgues's avatar Julien Lesgourgues
Browse files

includes dcdm, theta_s, new root system, faster tensors, new matlab plotting

parents 8318fa9f 8ed36f27
.. module:: CPU
:synopsis: CPU, a CLASS Plotting Utility
.. moduleauthor:: Benjamin Audren <>
.. version:: 2.0
This is a small python program aimed to gain time when comparing two spectra,
i.e. from CAMB and CLASS, or a non-linear spectrum to a linear one. It is
designed to be used in a command line fashion, not being restricted to your
CLASS directory, though it recognized mainly CLASS output format. Far from
perfect, or complete, it could use any suggestion for enhancing it, just to
avoid losing time on useless matters for others. Be warned that, when
comparing with other format, the following is assumed: there are no empty line
(especially at the end of file). Gnuplot comment lines (starting with a # are
allowed). This issue will cause a non-very descriptive error in CPU, any
suggestion for testing it is welcome. Example of use: To superimpose two
different spectra and see their global shape :
python output/lcdm_z2_pk.dat output/lncdm_z2_pk.dat
To see in details their ratio:
python output/lcdm_z2_pk.dat output/lncdm_z2_pk.dat -r
import numpy as np
import os
import matplotlib.pyplot as plt
import sys
import argparse
import itertools
START_LINE['error'] = [r' /|\ ',
r'/_o_\ ',
r' ']
START_LINE['warning'] = [r' /!\ ',
r' ']
START_LINE['info'] = [r' /!\ ',
r' ']
STANDARD_LENGTH = 80 # standard, increase if you have a big screen
def create_parser():
parser = argparse.ArgumentParser(
'CPU, a CLASS Plotting Utility, specify wether you want\n'
'to superimpose, or plot the ratio of different files.'),
'A standard usage would be, for instance:\n'
'python output/test_pk.dat output/test_pk_nl_density.dat'
' -r\npython output/wmap_cl.dat output/planck_cl.dat'),
'files', type=str, nargs='*', help='Files to plot')
parser.add_argument('-r', '--ratio', dest='ratio', action='store_true',
help='Plot the ratio of the spectra')
parser.add_argument('-s', '--selection', dest='selection',
help='specify the fields you want to plot.')
parser.add_argument('--scale', choices=['lin', 'loglog', 'loglin'],
help='Specify the scale to use for the plot')
'-p, --print',
dest='printfile', action='store_true', default=False,
help='print the graph directly in a .png file')
'-r, --repeat',
dest='repeat', action='store_true', default=False,
help='repeat the step for all redshifts with same base name')
return parser
def plot_CLASS_output(files, selection, ratio=False, output_name='',
extension='', x_variable='', scale='lin'):
Load the data to numpy arrays, write a Python script and plot them.
Inspired heavily by the matlab version by Thomas Tram
files : list
List of files to plot
selection : list, or string
List of items to plot, which should match the way they appear in the
file, for instance: ['TT', 'BB]
Keyword Arguments
ratio : bool
If set to yes, plots the ratio of the files, taking as a reference the
first one
output_name : str
Specify a different name for the produced figure (by default, it takes
the name of the first file, and replace the .dat by .pdf)
extension : str
# Load all the graphs
data = []
for data_file in files:
# Create the python script, and initialise it
python_script_path = files[0]+'.py'
text = '''
import matplotlib.pyplot as plt
import numpy as np\n'''
# Create the full_path_files list, that contains the absolute path, so that
# the future python script can import them directly.
full_path_files = [os.path.abspath(elem) for elem in files]
# Recover the base name of the files, everything before the .
roots = [elem.split(os.path.sep)[-1].split('.')[0] for elem in files]
text += '''files = %s\n''' % full_path_files
text += '''
data = []
for data_file in files:
# Recover the number of columns in the first file, as well as their title.
num_columns, names, tex_names = extract_headers(files[0])
# Check if everything is in order
if num_columns == 2:
selection = [names[1]]
elif num_columns > 2:
# in case selection was only a string, cast it to a list
if isinstance(selection, str):
selection = [selection]
for elem in selection:
if elem not in names:
raise InputError(
"The entry 'selection' must contain names of the fields "
"in the specified files. You asked for %s " % elem +
"where I only found %s." % names)
# Store the selected text and tex_names to the script
text += 'selection = %s\n' % selection
text += 'tex_names = %s\n' % [elem for (elem, name) in
zip(tex_names, names) if name in selection]
# Create the figure and ax objects
fig, ax = plt.subplots()
text += '\nfig, ax = plt.subplots()\n'
# if ratio is not set, then simply plot them all
if not ratio:
text += 'for curve in data:\n'
for idx, curve in enumerate(data):
_, curve_names, _ = extract_headers(files[idx])
for selec in selection:
index = curve_names.index(selec)
text += ' ax.'
if scale == 'lin':
text += 'plot(curve[:, 0], curve[:, %i])\n' % index
ax.plot(curve[:, 0], curve[:, index])
elif scale == 'loglog':
text += 'loglog(curve[:, 0], curve[:, %i])\n' % index
ax.loglog(curve[:, 0], curve[:, index])
ax.legend([root+': '+elem for (root, elem) in
itertools.product(roots, selection)], loc='lower right')
ref = data[0]
#for index in range(1, len(data)):
#current = data[index]
#if np.allclose(current[0], ref[0]):
#ax.plot(current[0], current[colnum]/ref[colnum])
text += '\n'
# Write to the python file all the issued commands. You can then reproduce
# the plot by running "python output/"
with open(python_script_path, 'w') as python_script:
class FormatError(Exception):
"""Format not recognised"""
class TypeError(Exception):
"""Spectrum type not recognised"""
class NumberOfFilesError(Exception):
"""Invalid number of files"""
class InputError(Exception):
"""Incompatible input requirements"""
def replace_scale(string):
This assumes that the string starts with "(.)", which will be replaced by
>>> print replace_scale('(.)toto')
>>> '(8\\pi G/3)toto'
string_list = list(string)
string_list[1:1] = list('8\\pi G/3')
return ''.join(string_list)
def process_long_names(long_names):
Given the names extracted from the header, return two arrays, one with the
short version, and one tex version
>>> names, tex_names = process_long_names(['(.)toto', 'proper time [Gyr]'])
>>> print names
>>> ['toto', 'proper time']
>>> print tex_names
>>> ['(8\\pi G/3)toto, 'proper time [Gyr]']
names = []
tex_names = []
# First pass, to remove the leading scales
for name in long_names:
# This can happen in the background file
if name.startswith('(.)', 0):
temp_name = name[3:]
# Otherwise, we simply
# Second pass, to remove from the short names the indication of scale,
# which should look like something between parenthesis, or square brackets,
# and located at the end of the string
for index, name in enumerate(names):
if name.find('(') != -1:
names[index] = name[:name.index('(')]
elif name.find('[') != -1:
names[index] = name[:name.index('[')]
return names, tex_names
def extract_headers(header_path):
with open(header_path, 'r') as header_file:
header = [line for line in header_file if line[0] == '#']
header = header[-1]
# Count the number of columns in the file, and recover their name. Thanks
# Thomas Tram for the trick
indices = [i+1 for i in range(len(header)) if
header.startswith(':', i)]
num_columns = len(indices)
long_names = [header[indices[i]:indices[(i+1)]-3].strip() if i < num_columns-1
else header[indices[i]:].strip()
for i in range(num_columns)]
# Process long_names further to handle special cases, and extract names,
# which will correspond to the tags specified in "selection".
names, tex_names = process_long_names(long_names)
return num_columns, names, tex_names
def main():
print '~~~ Running CPU, a CLASS Plotting Utility ~~~'
parser = create_parser()
# Parse the command line arguments
args = parser.parse_args()
# if there are no argument in the input, print usage
if len(args.files) == 0:
# 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'
scale = 'loglog'
elif args.files[0].rfind('pk') != -1:
selection = 'P'
scale = 'loglog'
raise TypeError(
"Please specify a field to plot")
args.selection = selection
scale = ''
if not args.scale:
if scale:
args.scale = scale
args.scale = 'lin'
# If ratio is asked, but only one file was passed in argument, politely
# complain
if args.ratio:
if len(args.files) < 2:
raise NumberOfFilesError(
"If you want me to compute a ratio between two files, "
"I strongly encourage you to give me at least two of them.")
# 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)
if __name__ == '__main__':
......@@ -24,6 +24,9 @@ CC = gcc -Wall #-g -pg -E#-ggdb
# your tool for creating static libraries:
AR = ar rv
# (OPT) your python interpreter
PYTHON = python
# your optimization flag
OPTFLAG = -O4 -ffast-math #-march=native
#OPTFLAG = -Ofast -ffast-math #-march=native
......@@ -122,12 +125,12 @@ H_ALL = $(addprefix include/, common.h svnversion.h $(addsuffix .h, $(basename $
PRE_ALL = cl_ref.pre clt_permille.pre
INI_ALL = explanatory.ini lcdm.ini
MISC_FILES = Makefile CPU psd_FD_single.dat myselection.dat myevolution.dat README bbn/sBBN.dat external_Pk/* cpp
PYTHON = python/classy.pyx python/ python/cclassy.pxd python/
PYTHON_FILES = python/classy.pyx python/ python/cclassy.pxd python/
all: class libclass.a
all: class libclass.a classy
libclass.a: $(TOOLS) $(SOURCE) $(EXTERNAL)
$(AR) $@ $(addprefix build/, $(TOOLS) $(SOURCE) $(EXTERNAL))
......@@ -162,9 +165,14 @@ test_thermodynamics: $(TOOLS) $(INPUT) $(BACKGROUND) $(THERMO) $(EXTERNAL) $(TES
$(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm
tar czvf class.tar.gz $(C_ALL) $(H_ALL) $(PRE_ALL) $(INI_ALL) $(MISC_FILES) $(HYREC) $(PYTHON)
tar czvf class.tar.gz $(C_ALL) $(H_ALL) $(PRE_ALL) $(INI_ALL) $(MISC_FILES) $(HYREC) $(PYTHON_FILES)
classy: libclass.a python/classy.pyx python/cclassy.pxd
cd python; $(PYTHON) install --user
clean: .base
rm -rf $(WRKDIR);
rm -f libclass.a
rm -f $(MDIR)/python/classy.c
rm -rf $(MDIR)/python/build
......@@ -20,6 +20,9 @@
#include <stdexcept>
//#define DBUG
using namespace std;
......@@ -67,16 +70,26 @@ ClassEngine::ClassEngine(const ClassParams& pars): cl(0),dofree(true){
//prepare fp structure
size_t n=pars.size();
for (size_t i=0;i<pars.size();i++){
//identify lmax
cout << pars.key(i) << "\t" << pars.value(i) <<endl;
if (pars.key(i)=="l_max_scalars") {
istringstream strstrm(pars.value(i));
strstrm >> _lmax;
cout << __FILE__ << " : using lmax=" << _lmax <<endl;
if (input_init(&fc,&pr,&ba,&th,&pt,&tr,&pm,&sp,&nl,&le,&op,_errmsg) == _FAILURE_)
if (input_init(&fc,&pr,&ba,&th,&pt,&tr,&pm,&sp,&nl,&le,&op,_errmsg) == _FAILURE_)
throw invalid_argument(_errmsg);
//proetction parametres mal defini
......@@ -86,7 +99,7 @@ ClassEngine::ClassEngine(const ClassParams& pars): cl(0),dofree(true){
//calcul class
//cout <<"creating " << sp.ct_size << " arrays" <<endl;
cl=new double[sp.ct_size];
......@@ -107,24 +120,32 @@ ClassEngine::ClassEngine(const ClassParams& pars,const string & precision_file):
struct file_content fc_input;
fc_input.size = 0;
fc_input.filename=new char[1];
//prepare fc par structure
size_t n=pars.size();
for (size_t i=0;i<pars.size();i++){
if (pars.key(i)=="l_max_scalars") {
istringstream strstrm(pars.value(i));
strstrm >> _lmax;
cout << __FILE__ << " : using lmax=" << _lmax <<endl;
//concatenate both
if (parser_cat(&fc_input,&fc_precision,&fc,_errmsg) == _FAILURE_) throw invalid_argument(_errmsg);
if (input_init(&fc,&pr,&ba,&th,&pt,&tr,&pm,&sp,&nl,&le,&op,_errmsg) == _FAILURE_)
if (input_init(&fc,&pr,&ba,&th,&pt,&tr,&pm,&sp,&nl,&le,&op,_errmsg) == _FAILURE_)
throw invalid_argument(_errmsg);
//proetction parametres mal defini
......@@ -134,7 +155,7 @@ ClassEngine::ClassEngine(const ClassParams& pars,const string & precision_file):
//calcul class
//cout <<"creating " << sp.ct_size << " arrays" <<endl;
cl=new double[sp.ct_size];
......@@ -142,7 +163,7 @@ ClassEngine::ClassEngine(const ClassParams& pars,const string & precision_file):
......@@ -161,13 +182,22 @@ ClassEngine::~ClassEngine()
// Member functions --
int ClassEngine::updateParValues(const std::vector<double>& par){
bool ClassEngine::updateParValues(const std::vector<double>& par){
dofree && freeStructs();
for (size_t i=0;i<par.size();i++) {
double val=par[i];
#ifdef DBUG
cout << "update par values #" << i << "\t" << val << "\t" << str(val).c_str() << endl;
return computeCls();
int status=computeCls();
#ifdef DBUG
cout << "update par status=" << status << " succes=" << _SUCCESS_ << endl;
return (status==_SUCCESS_);
//print content of file_content
......@@ -177,21 +207,20 @@ void ClassEngine::printFC() {
int ClassEngine::class(
struct file_content *pfc,
struct precision * ppr,
struct background * pba,
struct thermo * pth,
struct perturbs * ppt,
struct primordial * ppm,
struct nonlinear * pnl,
struct transfers * ptr,
struct spectra * psp,
struct lensing * ple,
struct output * pop,
ErrorMsg errmsg) {
int ClassEngine::class_main(
struct file_content *pfc,
struct precision * ppr,
struct background * pba,
struct thermo * pth,
struct perturbs * ppt,
struct transfers * ptr,
struct primordial * ppm,
struct spectra * psp,
struct nonlinear * pnl,
struct lensing * ple,
struct output * pop,
ErrorMsg errmsg) {
if (input_init(pfc,ppr,pba,pth,ppt,ptr,ppm,psp,pnl,ple,pop,errmsg) == _FAILURE_) {
printf("\n\nError running input_init_from_arguments \n=>%s\n",errmsg);
......@@ -275,7 +304,6 @@ int ClassEngine::class(
return _FAILURE_;
//fprintf(stderr,"%d %e %e %e\n",l,cl[l][0],cl[l][1],cl[l][2]);
return _SUCCESS_;
......@@ -284,35 +312,43 @@ int ClassEngine::class(
int ClassEngine::computeCls(){
#ifdef DBUG
cout <<"call computecls" << endl;
//new call
return this->class_main(&fc,&pr,&ba,&th,&pt,&tr,&pm,&sp,&nl,&le,&op,_errmsg);
int status=this->class_main(&fc,&pr,&ba,&th,&pt,&tr,&pm,&sp,&nl,&le,&op,_errmsg);
#ifdef DBUG
cout <<"status=" << status << endl;
return status;
if (lensing_free(&le) == _FAILURE_) {
printf("\n\nError in spectra_free \n=>%s\n",le.error_message);
return _FAILURE_;
if (nonlinear_free(&nl) == _FAILURE_) {
printf("\n\nError in nonlinear_free \n=>%s\n",nl.error_message);
return _FAILURE_;
if (spectra_free(&sp) == _FAILURE_) {
printf("\n\nError in spectra_free \n=>%s\n",sp.error_message);
return _FAILURE_;
if (primordial_free(&pm) == _FAILURE_) {
printf("\n\nError in primordial_free \n=>%s\n",pm.error_message);
return _FAILURE_;
if (transfer_free(&tr) == _FAILURE_) {
printf("\n\nError in transfer_free \n=>%s\n",tr.error_message);
return _FAILURE_;
......@@ -336,46 +372,13 @@ ClassEngine::freeStructs(){
return _SUCCESS_;
// int
// ClassEngine::l_size(Engine::cltype t){
// int lmax(-1);
// switch(t)
// {
// case TT:
// if (sp.has_tt==_TRUE_) lmax=sp.l_size[sp.index_ct_tt];
// break;
// case TE:
// if (sp.has_te==_TRUE_) lmax=sp.l_size[sp.index_ct_te] ;
// break;
// case EE:
// if (sp.has_ee==_TRUE_) lmax=sp.l_size[sp.index_ct_ee] ;
// break;
// case BB:
// if (sp.has_bb==_TRUE_) lmax=sp.l_size[sp.index_ct_bb] ;
// break;
// case PP:
// if (sp.has_pp==_TRUE_) lmax=sp.l_size[sp.index_ct_pp] ;
// break;
// case TP:
// if (sp.has_tp==_TRUE_) lmax=sp.l_size[sp.index_ct_tp] ;
// break;
// case EP:
// if (sp.has_ep==_TRUE_) lmax=sp.l_size[sp.index_ct_ep] ;
// break;
// }
// return lmax;
// }
ClassEngine::getCl(Engine::cltype t,const long &l){
if (!dofree) throw out_of_range("no Cl available because CLASS failed");
if (output_total_cl_at_l(&sp,&le,&op,static_cast<double>(l),cl) == _FAILURE_){
cerr << ">>>fail getting Cl type=" << (int)t << " @l=" << l <<endl;
cerr << ">>>fail getting Cl type=" << (int)t << " @l=" << l <<endl;