From 9629b244874caba9caa71c208a34390c0fc0f5ca Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 24 Feb 2015 17:41:30 +0100 Subject: [PATCH] Fixed bug in CosmoPower/ Eisenstein & Hu --- python/_cosmo_cic.pyx | 24 +++++++++++++++++------- python/_cosmo_power.pyx | 9 ++++++--- src/cic.cpp | 5 +++++ src/cosmopower.cpp | 23 +++++++++++++++-------- src/cosmopower.hpp | 2 +- 5 files changed, 44 insertions(+), 19 deletions(-) diff --git a/python/_cosmo_cic.pyx b/python/_cosmo_cic.pyx index 76a540a..cdb6a23 100644 --- a/python/_cosmo_cic.pyx +++ b/python/_cosmo_cic.pyx @@ -1,5 +1,6 @@ from libcpp cimport bool from libcpp cimport string as cppstring +from libcpp cimport vector as cppvector import numpy as np cimport numpy as np from cpython cimport PyObject, Py_INCREF @@ -29,7 +30,7 @@ cdef extern from "cic.hpp" namespace "CosmoTool": @cython.cdivision(True) @cython.wraparound(False) def leanCic(float[:,:] particles, float L, int Resolution): - cdef CICParticles p + cdef cppvector.vector[CICParticles] *p cdef CICFilter *cic cdef np.uint64_t i cdef CICType *field @@ -40,25 +41,34 @@ def leanCic(float[:,:] particles, float L, int Resolution): cdef np.uint64_t j cic = new CICFilter(Resolution, L) + print("Reset mesh") cic.resetMesh() if particles.shape[1] != 3: raise ValueError("Particles must be Nx3 array") - p.mass = 1 - for i in xrange(particles.shape[0]): - p.coords[0] = particles[i,0] - p.coords[1] = particles[i,1] - p.coords[2] = particles[i,2] - cic.putParticles(&p, 1) + print("Inserting particles") +# p = new cppvector.vector[CICParticles](particles.shape[0]) +# for i in xrange(particles.shape[0]): +# *p[i].mass = 1 +# *p[i].coords[0] = particles[i,0] +# *p[i].coords[1] = particles[i,1] +# *p[i].coords[2] = particles[i,2] + +# cic.putParticles(&p[0], particles.shape[0]) + del p + print("Done") field = 0 dummyRes = 0 cic.getDensityField(field, dummyRes) + + print("Got to allocate a numpy %dx%dx%d" % (dummyRes, dummyRes,dummyRes)) out_field = np.empty((dummyRes, dummyRes, dummyRes), dtype=np.float64) out_field0 = out_field.reshape(out_field.size) out_field_buf = out_field + print("Copy") for j in xrange(out_field_buf.size): out_field_buf[j] = field[j] diff --git a/python/_cosmo_power.pyx b/python/_cosmo_power.pyx index 2e87c5c..c23da1f 100644 --- a/python/_cosmo_power.pyx +++ b/python/_cosmo_power.pyx @@ -44,7 +44,7 @@ cdef extern from "cosmopower.hpp" namespace "CosmoTool": void setFunction(CosmoFunction) void updateCosmology() void updatePhysicalCosmology() - void normalize(double) + void normalize(double,double) void setNormalization(double) double power(double) @@ -74,8 +74,11 @@ cdef class CosmologyPower: assert self.power.OMEGA_C > 0 self.power.updateCosmology() + + def setNormalization(self,A): + self.power.setNormalization(A) - def normalize(self,s8,k_max=-1): + def normalize(self,s8,k_min=-1,k_max=-1): """normalize(self, sigma8) Compute the normalization of the power spectrum using sigma8. @@ -84,7 +87,7 @@ cdef class CosmologyPower: sigma8 (float): standard deviation of density field smoothed at 8 Mpc/h """ self.power.SIGMA8 = s8 - self.power.normalize(k_max) + self.power.normalize(k_min, k_max) def setFunction(self,funcname): diff --git a/src/cic.cpp b/src/cic.cpp index 31596ca..497bb49 100644 --- a/src/cic.cpp +++ b/src/cic.cpp @@ -186,6 +186,11 @@ void CICFilter::putParticles(CICParticles *particles, uint32_t N) for (long p = tUsedMin[t]; p < tUsedMax[t]; p++) densityGrid[p] += threadedDensity[t][p]; } + + for (int t = 0; t < threadUsed; t++) + { + delete[] threadedDensity; + } } void CICFilter::getDensityField(CICType*& field, uint32_t& res) diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp index 5104148..a7a90fb 100644 --- a/src/cosmopower.cpp +++ b/src/cosmopower.cpp @@ -90,7 +90,8 @@ static double powC(double q, double alpha_c) static double T_tilde_0(double q, double alpha_c, double beta_c) { - double a = log(M_E + 1.8 * beta_c * q); + static const double c_E = 2.718282; //M_E; + double a = log(c_E + 1.8 * beta_c * q); return a / ( a + powC(q, alpha_c) * q * q); } @@ -123,6 +124,7 @@ double CosmoPower::powerEfstathiou(double k) void CosmoPower::updateHuWigglesConsts() { + double k_silk = 1.6 * pow(OMEGA_B * h * h, 0.52) * pow(OmegaEff, 0.73) * (1 + pow(10.4 * OmegaEff, -0.95)); double z_eq = 2.50e4 * OmegaEff * pow(Theta_27, -4); //double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75))); @@ -130,7 +132,7 @@ void CosmoPower::updateHuWigglesConsts() double b1_zd = 0.313 * pow(OmegaEff, -0.419) * (1 + 0.607 * pow(OmegaEff, 0.674)); double b2_zd = 0.238 * pow(OmegaEff, 0.223); - double z_d = 1291 * pow(OmegaEff, 0.251) / (1 + 0.659 * pow(OmegaEff, 0.828)) * (1 + b1_zd * pow(OmegaEff, b2_zd)); + double z_d = 1291 * pow(OmegaEff, 0.251) / (1 + 0.659 * pow(OmegaEff, 0.828)) * (1 + b1_zd * pow(OMEGA_B*h*h, b2_zd)); double R_d = 31.5 * OMEGA_B * h * h * pow(Theta_27, -4) * 1e3 / z_d; double Req = 31.5 * OMEGA_B * h * h * pow(Theta_27, -4) * 1e3 / z_eq; @@ -180,7 +182,10 @@ double CosmoPower::powerHuWiggles(double k) double q = k / (13.41 * k_eq); double T_c = f * T_tilde_0(q, 1, beta_c) + (1 - f) * T_tilde_0(q, alpha_c, beta_c); - double T_b = (T_tilde_0(q, 1, 1) / (1 + pow(k * s / 5.2, 2)) + alpha_b / (1 + pow(beta_b / (k * s), 3)) * exp(-pow(k/k_silk, 1.4))) * j_0(k * s_tilde); + double T_b = ( + T_tilde_0(q, 1, 1) / (1 + pow(xx / 5.2, 2)) + + alpha_b / (1 + pow(beta_b / xx, 3)) * exp(-pow(k/k_silk, 1.4)) + ) * j_0(k * s_tilde); double T_k = OMEGA_B/OMEGA_0 * T_b + OMEGA_C/OMEGA_0 * T_c; @@ -262,16 +267,18 @@ double CosmoPower::integrandNormalize(double x) return power(k)*k*k*f*f/(x*x); } -void CosmoPower::normalize(double k_max) +void CosmoPower::normalize(double k_min, double k_max) { double normVal = 0; double abserr; gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION); gsl_function f; - double x_min = 0; + double x_min = 0, x_max = 1; if (k_max > 0) x_min = 1/(1+k_max); + if (k_min > 0) + x_max = 1/(1+k_min); f.function = gslPowSpecNorm; f.params = this; @@ -281,12 +288,12 @@ void CosmoPower::normalize(double k_max) ofstream ff("PP_k.txt"); for (int i = 0; i < 100; i++) { - double k = pow(10.0, 4.0*i/100.-2); + double k = pow(10.0, 8.0*i/100.-4); ff << k << " " << power(k) << endl; } // gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr); - gsl_integration_qag(&f, x_min, 1, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr); + gsl_integration_qag(&f, x_min, x_max, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr); gsl_integration_workspace_free(w); normVal /= (2*M_PI*M_PI); @@ -371,5 +378,5 @@ void CosmoPower::setFunction(CosmoFunction f) void CosmoPower::setNormalization(double A_K) { - normPower = A_K/power(0.002); + normPower = A_K;///power(0.002); } diff --git a/src/cosmopower.hpp b/src/cosmopower.hpp index f1ba45d..b32c6a7 100644 --- a/src/cosmopower.hpp +++ b/src/cosmopower.hpp @@ -92,7 +92,7 @@ namespace CosmoTool { void updateCosmology(); void updatePhysicalCosmology(); - void normalize(double k_max = -1); + void normalize(double k_min = -1, double k_max = -1); void setNormalization(double A_K); void updateHuWigglesConsts();