Fancyvoigtfit

.gitignore

@@ -13,3 +13,9 @@ fake_spectra.egg-info/SOURCES.txt
 fake_spectra.egg-info/PKG-INFO
 fake_spectra.egg-info/dependency_links.txt
 fake_spectra/*.so
+.*.cache
+*.pyc
+build.log
+dist/
+.pytest_cache/
+cextract/*.o

Makefile

@@ -54,7 +54,7 @@ cextract/build/read_snapshot.o: cextract/read_snapshot.cpp cextract/global_vars.
 btest/%.o: fake_spectra/%.cpp fake_spectra/%.h
 	$(CXX) $(CFLAGS) -c $< -o $@
 
-btest/absorption.o: fake_spectra/absorption.cpp fake_spectra/absorption.h fake_spectra/singleabs.h
+btest/absorption.o: fake_spectra/absorption.cpp fake_spectra/absorption.h fake_spectra/singleabs.h fake_spectra/fastvoigt.h
 btest/part_int.o: fake_spectra/part_int.cpp fake_spectra/part_int.h fake_spectra/absorption.h fake_spectra/index_table.h
 
 clean:

README.md

@@ -19,11 +19,9 @@ how to use it.
 
 Required Python libraries:
 - numpy (core functionality)
+- scipy (core functionality)
 - h5py (for saving)
 
-Required C libraries:
-- GSL
-
 Optional libraries:
 - matplotlib (if you want to plot)
 - bigfile (to install, do 'pip install --user bigfile') for reading BigFile snapshot outputs from Yu Feng's MP-Gadget.

cextract/global_vars.h

@@ -9,13 +9,13 @@
 
 struct particle_data
 {
-  float *Pos;
-  float *Vel;
-  float *Mass;
-  float *U, *Ne;
-  float *temp;
-  float *h;
-  float *fraction;
+  double *Pos;
+  double *Vel;
+  double *Mass;
+  double *U, *Ne;
+  double *temp;
+  double *h;
+  double *fraction;
 };
 typedef struct particle_data pdata;
 

cextract/init.c

@@ -56,14 +56,14 @@ void populate_los_table(double * cofm, int * aaxis, int NumLos, char * ext_table
 /*Note this assumes only one species*/
 int alloc_parts(pdata* P, int np)
 {
-    return ((*P).Vel=(float *)malloc(np*3*sizeof(float))) &&
-    ((*P).Pos=(float *)malloc(np*3*sizeof(float))) &&
-     ((*P).Mass=(float *) malloc(np*sizeof(float))) &&
-    ((*P).U=(float *)malloc(np*sizeof(float))) &&
-    ((*P).fraction=(float *)malloc(np*sizeof(float))) &&
-    ((*P).temp=(float *)malloc(np*sizeof(float))) &&
-    ((*P).Ne=(float *)malloc(np*sizeof(float))) &&
-    ((*P).h=(float *)malloc(np*sizeof(float)));
+    return ((*P).Vel=(double *)malloc(np*3*sizeof(double))) &&
+    ((*P).Pos=(double *)malloc(np*3*sizeof(double))) &&
+     ((*P).Mass=(double *) malloc(np*sizeof(double))) &&
+    ((*P).U=(double *)malloc(np*sizeof(double))) &&
+    ((*P).fraction=(double *)malloc(np*sizeof(double))) &&
+    ((*P).temp=(double *)malloc(np*sizeof(double))) &&
+    ((*P).Ne=(double *)malloc(np*sizeof(double))) &&
+    ((*P).h=(double *)malloc(np*sizeof(double)));
 }
 
 void free_parts(pdata* P)

cextract/read_hdf_snapshot.c

@@ -23,7 +23,9 @@ hsize_t get_single_dataset(const char *name, void * data_ptr,  hsize_t data_leng
              return 0;
           }
           H5LTget_dataset_info(*hdf_group, name, &vlength, &class_id, &type_size);
-          if(type_size != 4 || class_id != H5T_FLOAT  || vlength > data_length || H5LTread_dataset_float(*hdf_group, name, data_ptr) < 0 ){
+          /*The interpolation wants double: HDF5 widens a single precision
+           * dataset on read, so both sizes are fine here.*/
+          if((type_size != 4 && type_size != 8) || class_id != H5T_FLOAT  || vlength > data_length || H5LTread_dataset_double(*hdf_group, name, data_ptr) < 0 ){
               fprintf(stderr, "File %d: Failed reading %s (%lu)\n",fileno,name, (uint64_t)vlength);
               return 0;
           }
@@ -41,7 +43,7 @@ hsize_t get_triple_dataset(const char *name, void * data_ptr, hsize_t data_lengt
              return 0;
           }
           H5LTget_dataset_info(*hdf_group, name, &vlength[0], &class_id, &type_size);
-          if(type_size != 4 || class_id != H5T_FLOAT || vlength[1] != 3 || vlength[0] > data_length || H5LTread_dataset_float(*hdf_group, name, data_ptr) < 0 ){
+          if((type_size != 4 && type_size != 8) || class_id != H5T_FLOAT || vlength[1] != 3 || vlength[0] > data_length || H5LTread_dataset_double(*hdf_group, name, data_ptr) < 0 ){
               fprintf(stderr, "File %d: Failed reading %s (%lu)\n",fileno,name, (uint64_t)vlength[0]);
               return 0;
           }
@@ -213,7 +215,9 @@ int load_hdf5_snapshot(const char *ffname, pdata *P, int fileno)
 exit:
   H5Gclose(hdf_group);
   H5Fclose(hdf_file);
-  if(fileno < 1){
+  /*Note length is unsigned, so length-1 is not a valid index when nothing was
+   * read: we get here with length == 0 whenever reading the file failed.*/
+  if(fileno < 1 && length > 0){
         printf("\nP[%d].Pos = [%g %g %g]\n", 0, (*P).Pos[0], (*P).Pos[1],(*P).Pos[2]);
         printf("P[%d].Vel = [%g %g %g]\n", 0, (*P).Vel[0], (*P).Vel[1],(*P).Vel[2]);
         printf("P[-1].Density = %e\n", (*P).Mass[0]);

cextract/read_snapshot.cpp

@@ -16,13 +16,25 @@
 #ifndef NOGREAD
 
 #include <math.h>
+#include <vector>
 #include "global_vars.h"
 #include "gadgetreader.hpp"
 
 
 /* Snapshot information */
 #define PARTTYPE 0/* Particle type required */
 
+/* GadgetReader hands back the single precision the file holds, but the
+ * interpolation wants double, so read into a scratch buffer and widen.
+ * ncomp is how many values a particle has in this block: 3 for a position. */
+static void GetBlockDouble(GadgetReader::GSnap& snap, const char * name, double * out, int64_t npart, int ncomp, int64_t start, int mask)
+{
+  std::vector<float> buf((size_t)npart*ncomp);
+  snap.GetBlock(name, &buf[0], npart, start, mask);
+  for(int64_t i = 0; i < npart*ncomp; i++)
+      out[i] = buf[i];
+}
+
 int load_header(const char *fname,double  *atime, double *redshift, double * Hz, double *box100, double *h100)
 {
 #ifndef SPLIT_NE
@@ -70,10 +82,10 @@ int64_t load_snapshot(const char *fname,int64_t StartPart, pdata *P)
     exit(1);
   }
   printf("Reading from %ld to %ld\n",StartPart,StartPart+NumPart);
-  snap.GetBlock("POS ",(*P).Pos,NumPart,StartPart, (1<<N_TYPE)-1-(1<<PARTTYPE));
-  snap.GetBlock("VEL ",(*P).Vel,NumPart,StartPart, (1<<N_TYPE)-1-(1<<PARTTYPE));
+  GetBlockDouble(snap, "POS ",(*P).Pos,NumPart,3,StartPart, (1<<N_TYPE)-1-(1<<PARTTYPE));
+  GetBlockDouble(snap, "VEL ",(*P).Vel,NumPart,3,StartPart, (1<<N_TYPE)-1-(1<<PARTTYPE));
   /* Particle densities */
-  snap.GetBlock("RHO ",(*P).Mass,NumPart,StartPart, (1<<N_TYPE)-1-(1<<PARTTYPE));
+  GetBlockDouble(snap, "RHO ",(*P).Mass,NumPart,1,StartPart, (1<<N_TYPE)-1-(1<<PARTTYPE));
   for(int i=0; i< NumPart;i++)
   if ((*P).Mass[i] != (*P).Mass[0]){
         fprintf(stderr, "i=%d N = %ld Mass change: %g\n",i,NumPart, (*P).Mass[i]);
@@ -83,7 +95,7 @@ int64_t load_snapshot(const char *fname,int64_t StartPart, pdata *P)
   if(PARTTYPE == 0)
     { 
       /*The internal energy of all the Sph particles is read in */
-      snap.GetBlock("U   ",(*P).U,NumPart,StartPart,0);
+      GetBlockDouble(snap, "U   ",(*P).U,NumPart,1,StartPart,0);
       /* The free electron fraction */
       if(snap.GetHeader().flag_cooling)
         {
@@ -92,29 +104,29 @@ int64_t load_snapshot(const char *fname,int64_t StartPart, pdata *P)
           /* Use that the universe is neutral, so 
            * NE = NHP + NHEP +2 NHEPP*/
       #ifdef SPLIT_NE
-          snap.GetBlock("NHP ",(*P).Ne,NumPart,StartPart,0);
+          GetBlockDouble(snap, "NHP ",(*P).Ne,NumPart,1,StartPart,0);
           /*Use the space for HSML as temp space*/
-          snap.GetBlock("NHEP",(*P).h,NumPart,StartPart,0);
+          GetBlockDouble(snap, "NHEP",(*P).h,NumPart,1,StartPart,0);
           for(int k=0;k<NumPart;k++){
                   (*P).Ne[k]+=(*P).h[k];
           }
-          snap.GetBlock("NHEQ",(*P).h,NumPart,StartPart,0);
+          GetBlockDouble(snap, "NHEQ",(*P).h,NumPart,1,StartPart,0);
           for(int k=0;k<NumPart;k++){
                   (*P).Ne[k]+=2*(*P).h[k];
           }
       #else
-          snap.GetBlock("NE  ",(*P).Ne,NumPart,StartPart,0);
+          GetBlockDouble(snap, "NE  ",(*P).Ne,NumPart,1,StartPart,0);
       #endif
       #ifdef HELIUM
-          snap.GetBlock("NHE ",(*P).NHep,NumPart,StartPart,0);
+          GetBlockDouble(snap, "NHE ",(*P).NHep,NumPart,1,StartPart,0);
       #endif 
           /* The HI fraction, nHI/nH */
-          snap.GetBlock("NH  ",(*P).fraction,NumPart,StartPart,0);
+          GetBlockDouble(snap, "NH  ",(*P).fraction,NumPart,1,StartPart,0);
           /*An NHE block*/
         }
      /* The smoothing length */
-     snap.GetBlock("HSML",(*P).h,NumPart,StartPart,0);
-     for(i=0;i<NumPart;i++)
+     GetBlockDouble(snap, "HSML",(*P).h,NumPart,1,StartPart,0);
+     for(int64_t i=0;i<NumPart;i++)
              (*P).h[i] = (*P).h[i]/2;
     }
 

fake_spectra/absorption.cpp

@@ -3,7 +3,7 @@
  * Permission to use, copy, modify, and/or distribute this software for any
  * purpose with or without fee is hereby granted, provided that the above
  * copyright notice and this permission notice appear in all copies.
- * 
+ *
  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
@@ -32,6 +32,9 @@
 /*Conversion factor between internal energy and mu and temperature in K */
 #define TSCALE ((GAMMA-1.0) * PROTONMASS * ESCALE / BOLTZMANN)
 
+/*The tabulated Voigt function, built once before main().*/
+const VoigtTable voigt_table;
+
 /* Find the integral of the particle density in this pixel by integrating an SPH kernel
  * over the z direction.
  * Arguments:
@@ -164,7 +167,7 @@ kernel(kernel_i)
  * tau, and the density from the particle to the array colden
  * The slightly C-style interface is so we can easily use the data in python
  */
-void LineAbsorption::add_colden_particle(double * colden, const int nbins, const double dr2, const float dens, const float pos, const float smooth)
+void LineAbsorption::add_colden_particle(double * colden, const int nbins, const double dr2, const double dens, const double pos, const double smooth)
 {
   double pos1 = pos;
   if(kernel == VORONOI_MESH)
@@ -209,7 +212,7 @@ void LineAbsorption::add_colden_particle(double * colden, const int nbins, const
   }
 }
 
-void LineAbsorption::add_tau_particle(double * tau, const int nbins, const double dr2, const float dens, const float ppos, const float pvel, const float temp, const float smooth)
+void LineAbsorption::add_tau_particle(double * tau, const int nbins, const double dr2, const double dens, const double ppos, const double pvel, const double temp, const double smooth)
 {
   /* Velocity of particle parallel to los: pos in kpc/h comoving
      to vel in km/s physical. Note that gadget velocities come comoving,

fake_spectra/absorption.h

@@ -38,15 +38,15 @@ class LineAbsorption
          *
          * Input:
          * dr2: transverse distance to spectra from particle (comoving kpc/h)
-         * dens: value of density field for absorbing species (1e10 M_sun / h / (comoving kpc/h)^3 from C extractor 
+         * dens: value of density field for absorbing species (1e10 M_sun / h / (comoving kpc/h)^3 from C extractor
          * or amu/cm^2 kpc/h from the python module)
          * ppos: particle distance from box edge parallel to spectrum (comoving kpc/h)
          * pvel: particle velocity parallel to spectrum (physical km/s) (not needed for colden)
          * temp: particle temperature (K) (not needed for colden)
          * smooth: particle smoothing length (comoving kpc/h)
          */
-        void add_colden_particle(double * colden, const int nbins, const double dr2, const float dens, const float ppos, const float smooth);
-        void add_tau_particle(double * tau, const int nbins, const double dr2, const float dens, const float ppos, const float pvel, const float temp, const float smooth);
+        void add_colden_particle(double * colden, const int nbins, const double dr2, const double dens, const double ppos, const double smooth);
+        void add_tau_particle(double * tau, const int nbins, const double dr2, const double dens, const double ppos, const double pvel, const double temp, const double smooth);
 
     private:
         /* Threshold of tau below which we stop computing profiles.

fake_spectra/emission.py

@@ -14,6 +14,7 @@
 import numpy as np
 
 from . import spectra as ss
+from .near_lines import near_lines
 
 def maginJy(mag, band):
     """Convert a magnitude to flux in Jansky, according to wikipedia's table"""
@@ -51,7 +52,7 @@ def _read_stellar_data(self,fn, band, hhmult=10.):
         #Set each stellar radius to the pixel size
         hh = hhmult*np.ones(np.shape(pos)[0], dtype=np.float32)
         #Find particles we care about
-        ind = self.particles_near_lines(pos, hh,self.axis,self.cofm)
+        ind = near_lines(self.box, pos, hh, self.axis, self.cofm, pool=self.pool)
         #print np.size(ind)
         #Do nothing if there aren't any, and return a suitably shaped zero array
         if np.size(ind) == 0:
@@ -60,7 +61,7 @@ def _read_stellar_data(self,fn, band, hhmult=10.):
         hh = hh[ind]
         #Find the magnitude of stars in this band
         emflux = maginJy(self.snapshot_set.get_data(4,"GFM_StellarPhotometrics", segment = fn).astype(np.float32)[ind][:,nband],band)
-        fluxx = np.array([ np.sum(emflux[self.particles_near_lines(pos, hh,np.array([ax,]),np.array([cofm,]))]) for (ax, cofm) in zip(self.axis, self.cofm)])
+        fluxx = np.array([ np.sum(emflux[near_lines(self.box, pos, hh, np.array([ax,]), np.array([cofm,]), pool=self.pool)]) for (ax, cofm) in zip(self.axis, self.cofm)])
         #print np.sum(emflux)
         return fluxx
         #return (pos, emflux, hh)
@@ -139,7 +140,7 @@ def calculator(H0, Omega_M, zz):
     a = np.logspace(np.log10(az), 0, n)
     a2H = a*a*np.sqrt(Omega_M/a**3+WR/(a**4)+WV)
     #Comoving distance
-    DCMR = np.trapz(1./a2H, a)
+    DCMR = np.trapezoid(1./a2H, a)
     #In Mpc
     DC_Mpc = (light/H0) * DCMR
     # angular size distance In Mpc

fake_spectra/fastvoigt.h

@@ -0,0 +1,82 @@
+#ifndef FAST_VOIGT_H
+#define FAST_VOIGT_H
+
+#include <cmath>
+#include "Faddeeva.h"
+
+/* A fast evaluation of the Voigt function H(a,u) = Re[w(u + i a)], which is
+ * where almost all of the time in compute_tau goes.
+ *
+ * The damping parameter a = gamma lambda / (4 pi b) is small: 5e-4 for
+ * Lyman alpha at 10^4 K, and below 3e-3 for any line and temperature we
+ * encounter. In that regime the expansion in a (Harris 1948) is essentially
+ * exact:
+ *
+ *      H(a,u) = exp(-u^2) + a H1(u) + O(a^2)
+ *      H1(u)  = (2/sqrt(pi)) (2 u Daw(u) - 1) = 2 u Im[w(u)] - 2/sqrt(pi)
+ *
+ * The neglected term is a^2 (1 - 2u^2) exp(-u^2), so for a = 5e-4 the error
+ * is 2.5e-7 absolute and 6e-6 relative. Above VOIGT_AMAX we call Faddeeva
+ * rather than extrapolate the expansion into a regime it is not good for.
+ *
+ * H1 is smooth, O(1), and multiplied by a small number, so tabulating it and
+ * interpolating linearly leaves an error well below the a^2 truncation.
+ * The Gaussian is tabulated on the same grid: exp() is otherwise the most
+ * expensive thing left once Faddeeva is gone.
+ *
+ * Past VOIGT_UMAX the Gaussian has underflowed and only the asymptotic series
+ * for H1 is needed. That is where the damping wings are, which is where most
+ * of the bins are for a strong absorber, so those bins cost a few flops. */
+
+#define VOIGT_UMAX 6.0
+#define VOIGT_NTAB 1024
+/* Above this a the expansion is not accurate enough and we call Faddeeva. */
+#define VOIGT_AMAX 0.001
+
+struct VoigtTable
+{
+    /*The Gaussian, the first order damping correction, and the difference
+     * between neighbouring entries of the latter (so interpolating is a
+     * multiply-add rather than two loads and a subtraction).*/
+    double gau[VOIGT_NTAB+2];
+    double h1[VOIGT_NTAB+2];
+    double dh1[VOIGT_NTAB+2];
+    VoigtTable()
+    {
+        for(int i = 0; i < VOIGT_NTAB+2; ++i) {
+            const double u = i*(VOIGT_UMAX/VOIGT_NTAB);
+            gau[i] = exp(-u*u);
+            h1[i] = 2*u*Faddeeva::w_im(u) - 2/sqrt(M_PI);
+        }
+        for(int i = 0; i < VOIGT_NTAB+1; ++i)
+            dh1[i] = h1[i+1] - h1[i];
+        dh1[VOIGT_NTAB+1] = 0;
+    }
+};
+
+/*Built once, at startup: defined in absorption.cpp.*/
+extern const VoigtTable voigt_table;
+
+/* H(a,u) to first order in a. Only call this for a < VOIGT_AMAX;
+ * profile() in singleabs.h does the check. */
+inline double fast_voigt(const double uu, const double aa)
+{
+    const double u = fabs(uu);
+    if (u >= VOIGT_UMAX) {
+        /*Asymptotic series for H1: the Gaussian is < 3e-16 here.*/
+        const double u2 = 1./(u*u);
+        return aa*M_2_SQRTPI/2*u2*(1 + u2*(1.5 + u2*(3.75 + u2*(13.125 + 59.0625*u2))));
+    }
+    const double x = u*(VOIGT_NTAB/VOIGT_UMAX);
+    const int i = (int) x;
+    const double f = x - i;
+    /*exp(-u^2) = exp(-u_i^2) exp(-(2 u_i d + d^2)) where d is the distance to
+     * the tabulated point below. The second factor is within 0.08 of one, so
+     * a short Taylor series for it is good to 1e-8.*/
+    const double d = f*(VOIGT_UMAX/VOIGT_NTAB);
+    const double t = -d*(2*i*(VOIGT_UMAX/VOIGT_NTAB) + d);
+    const double ex = voigt_table.gau[i]*(1+t*(1+t*(0.5+t*(1./6+t*(1./24)))));
+    return ex + aa*(voigt_table.h1[i] + f*voigt_table.dh1[i]);
+}
+
+#endif