Commit b5c4ec5c authored by Plaszczynski Stephane's avatar Plaszczynski Stephane
Browse files

modify hasPk_NL to match class structure

parent 753f34a2
......@@ -97,6 +97,7 @@ public:
virtual double get_PkNL(double k, double z) { return klass->get_PkNL(k,z);}
inline bool get_Pkvec(const std::vector<double>& knodes, const std::vector<double>& zvec, std::vector<double>& pks) {return klass->get_Pkvec(knodes,zvec,pks);}
inline bool get_Pkvec(const std::vector<double>& knodes, const double& z, std::vector<double>& pks) {return klass->get_Pkvec(knodes,z,pks);}
void writeCls(std::ostream &o) {klass->writeCls(o);}
......
......@@ -40,7 +40,7 @@ public:
virtual bool has_CMB_Lensing() const {return false;}
virtual bool has_Pklin() const {return false;}
virtual bool has_PkNL() const {return _nonlin;} //modify _nonlin in derived class
virtual bool has_PkNL() const {return false;}
virtual std::string name() const=0;
......@@ -87,7 +87,7 @@ public:
//P(K): convention if NL is defined return it it, oterwise lin
inline virtual double get_Pk(double k,double z=0)
{
return (_nonlin ? get_PkNL(k,z) : get_Pklin(k,z));
return (has_PkNL() ? get_PkNL(k,z) : get_Pklin(k,z));
}
//rather prefer explicit calls
virtual double get_Pklin(double k, double z) { return undef("get_Pklin");}
......
#include "Parser.hh"
#include "Timer.hh"
#include "Class/MnClassEngine.hh"
#include"cxxsupport/cxxutils.h"
#include "Timer.hh"
//#include"cxxsupport/arr.h"
#include"cxxsupport/fitshandle.h"
#include <iostream>
#include<fstream>
#include<cassert>
#include<numeric>
#include<string>
#include<cmath>
using namespace std;
......@@ -19,9 +16,13 @@ Utlities to check the CLASS Pk spectrum: fixed nodes location + inetrpolated val
author : S. Plaszczynski, Aug.2017
*/
static bool abs_compare(double a, double b)
{
return (std::abs(a) < std::abs(b));
};
void logspace(double kmin,double kmax, int npks,vector<double>&v){
static void logspace(double kmin,double kmax, int npks,vector<double>&v){
double lmax=std::log(kmax);
double lmin=std::log(kmin);
......@@ -35,7 +36,7 @@ void logspace(double kmin,double kmax, int npks,vector<double>&v){
}
void linspace(double kmin,double kmax, int npks,vector<double>&v){
static void linspace(double kmin,double kmax, int npks,vector<double>&v){
double step=(kmax-kmin)/(npks-1);
v.resize(npks);
......@@ -129,8 +130,6 @@ int main(int argc,char** argv){
}
//precision file
string pre="";
if (parser.params.param_present("precisionFile")){
......@@ -147,7 +146,34 @@ int main(int argc,char** argv){
cout << "Computing P(k) at z=" << z0 << endl;
ClassEngine* klass=e->getClass();
//output
cout << "sigma8(0)=" << e->get_sigma8(0) << "\tsigma8(" << z0 << ")=" << e->get_sigma8(z0) << endl;
//test fast Pk
Timer timer;
const size_t nloop=100;
cout << "hasPK? " << (int)e->has_PkNL() << "\t" << klass->has_PkNL() << endl;
cout <<"Pk timer test..." << endl;
vector<double> pkslow(kval.size());
for (size_t i=0;i<nloop;i++){
for (size_t k =0; k<kval.size() ; k++){
pkslow[k]=e->get_Pk(kval[k], z0);
}
}
cout << "std=" << timer.partial() <<endl;
vector<double> pkfast(kval.size());
for (size_t i=0;i<nloop;i++) e->get_Pkvec(kval,z0,pkfast);
cout << "fast=" << timer.partial() <<endl;
//check ouptut
vector<double> residual(kval.size());
std::transform(pkfast.begin(),pkfast.end(),pkslow.begin(),
residual.begin(),std::minus<double>());
double diff=*max_element(residual.begin(),residual.end(),abs_compare);
cout << "max diff=" << diff << endl;
//fits output
const std::string fileName(argv[2]);
cout << "writing FITS file =" << fileName << endl;
if (file_present(fileName)) {
......@@ -162,7 +188,7 @@ int main(int argc,char** argv){
{
std::vector<fitscolumn> pkcols;
pkcols.push_back(fitscolumn("knode", "1/Mpc",1,TDOUBLE));
pkcols.push_back(fitscolumn("pknode", "(Mpc)^3" ,1,TDOUBLE));
pkcols.push_back(fitscolumn("pklinnode", "(Mpc)^3" ,1,TDOUBLE));
if (do_nl) pkcols.push_back(fitscolumn("pknodenl", "(Mpc)^3" ,1,TDOUBLE));
fout.insert_bintab(pkcols);
......@@ -181,25 +207,6 @@ int main(int argc,char** argv){
}
//interp values
//test fast Pk
/*
Timer timer;
cout <<"Pk timer test..." << endl;
vector<double> pkvec(kval.size());
for (size_t i=0;i<100;i++)
e->get_Pkvec(kval,z0,pkvec);
cout << "fast=" << timer.partial() <<endl;
for (size_t i=0;i<100;i++){
for (int k =0; k<kval.size() ; k++){
double pk=e->get_Pklin(kval[k], 0);
}
}
cout << "std=" << timer.partial() <<endl;
*/
{
std::vector<fitscolumn> pkcols;
pkcols.push_back(fitscolumn("k", "1/Mpc",1,TDOUBLE));
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment