* Re-enable support for fitting multipole datasets

* Multipole fits now write results in a format suitable for input to a multipole fit * General cleanup of obsolete parts of the code and the BOSS DR9 Lya legacy fit * Update demo/ with a companion wiki page http://darkmatter.ps.uci.edu/wiki/DeepZot/Baofit/Validation * Minor bug fix that affected error on beta (except for cross correlation fits)
igmhub · Jun 14, 2013 · 34320ef · 34320ef
2 parents e80bfd3 + 8151be1
commit 34320ef
Show file tree

Hide file tree

Showing 45 changed files with 9,412 additions and 41,822 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1,3 +1,5 @@
 *~
 # ignore files automatically generated by running autoreconf
 autom4te.cache
+# ignore sublime text editor config files
+*.sublime-*
diff --git a/Makefile.am b/Makefile.am
@@ -31,7 +31,6 @@ libbaofit_la_SOURCES = \
 	baofit/AbsCorrelationData.cc \
 	baofit/QuasarCorrelationData.cc \
 	baofit/ComovingCorrelationData.cc \
-	baofit/MultipoleCorrelationData.cc \
 	baofit/CorrelationFitter.cc \
 	baofit/CorrelationAnalyzer.cc \
 	baofit/boss.cc
@@ -51,7 +50,6 @@ nobase_include_HEADERS = \
 	baofit/AbsCorrelationData.h \
 	baofit/QuasarCorrelationData.h \
 	baofit/ComovingCorrelationData.h \
-	baofit/MultipoleCorrelationData.h \
 	baofit/CorrelationFitter.h \
 	baofit/CorrelationAnalyzer.h \
 	baofit/boss.h

diff --git a/Makefile.in b/Makefile.in
@@ -84,8 +84,7 @@ am_libbaofit_la_OBJECTS = AbsCorrelationModel.lo \
 	BaoCorrelationModel.lo BroadbandModel.lo XiCorrelationModel.lo \
 	PkCorrelationModel.lo AbsCorrelationData.lo \
 	QuasarCorrelationData.lo ComovingCorrelationData.lo \
-	MultipoleCorrelationData.lo CorrelationFitter.lo \
-	CorrelationAnalyzer.lo boss.lo
+	CorrelationFitter.lo CorrelationAnalyzer.lo boss.lo
 libbaofit_la_OBJECTS = $(am_libbaofit_la_OBJECTS)
 PROGRAMS = $(bin_PROGRAMS)
 am_baofit_OBJECTS = baofit.$(OBJEXT)
@@ -280,7 +279,6 @@ libbaofit_la_SOURCES = \
 	baofit/AbsCorrelationData.cc \
 	baofit/QuasarCorrelationData.cc \
 	baofit/ComovingCorrelationData.cc \
-	baofit/MultipoleCorrelationData.cc \
 	baofit/CorrelationFitter.cc \
 	baofit/CorrelationAnalyzer.cc \
 	baofit/boss.cc
@@ -301,7 +299,6 @@ nobase_include_HEADERS = \
 	baofit/AbsCorrelationData.h \
 	baofit/QuasarCorrelationData.h \
 	baofit/ComovingCorrelationData.h \
-	baofit/MultipoleCorrelationData.h \
 	baofit/CorrelationFitter.h \
 	baofit/CorrelationAnalyzer.h \
 	baofit/boss.h
@@ -462,7 +459,6 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ComovingCorrelationData.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CorrelationAnalyzer.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CorrelationFitter.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/MultipoleCorrelationData.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PkCorrelationModel.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/QuasarCorrelationData.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/XiCorrelationModel.Plo@am__quote@
@@ -546,13 +542,6 @@ ComovingCorrelationData.lo: baofit/ComovingCorrelationData.cc
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(LIBTOOL)  --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o ComovingCorrelationData.lo `test -f 'baofit/ComovingCorrelationData.cc' || echo '$(srcdir)/'`baofit/ComovingCorrelationData.cc
 
-MultipoleCorrelationData.lo: baofit/MultipoleCorrelationData.cc
-@am__fastdepCXX_TRUE@	$(LIBTOOL)  --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT MultipoleCorrelationData.lo -MD -MP -MF $(DEPDIR)/MultipoleCorrelationData.Tpo -c -o MultipoleCorrelationData.lo `test -f 'baofit/MultipoleCorrelationData.cc' || echo '$(srcdir)/'`baofit/MultipoleCorrelationData.cc
-@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/MultipoleCorrelationData.Tpo $(DEPDIR)/MultipoleCorrelationData.Plo
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='baofit/MultipoleCorrelationData.cc' object='MultipoleCorrelationData.lo' libtool=yes @AMDEPBACKSLASH@
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
-@am__fastdepCXX_FALSE@	$(LIBTOOL)  --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o MultipoleCorrelationData.lo `test -f 'baofit/MultipoleCorrelationData.cc' || echo '$(srcdir)/'`baofit/MultipoleCorrelationData.cc
-
 CorrelationFitter.lo: baofit/CorrelationFitter.cc
 @am__fastdepCXX_TRUE@	$(LIBTOOL)  --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT CorrelationFitter.lo -MD -MP -MF $(DEPDIR)/CorrelationFitter.Tpo -c -o CorrelationFitter.lo `test -f 'baofit/CorrelationFitter.cc' || echo '$(srcdir)/'`baofit/CorrelationFitter.cc
 @am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/CorrelationFitter.Tpo $(DEPDIR)/CorrelationFitter.Plo

diff --git a/baofit/AbsCorrelationData.cc b/baofit/AbsCorrelationData.cc
@@ -4,9 +4,18 @@
 #include "baofit/RuntimeError.h"
 
 #include "likely/CovarianceMatrix.h"
+#include "likely/AbsBinning.h"
+#include "likely/RuntimeError.h"
 
 #include "boost/lexical_cast.hpp"
+#include "boost/algorithm/string.hpp"
+#include "boost/spirit/include/qi.hpp"
+#include "boost/spirit/include/phoenix_core.hpp"
+#include "boost/spirit/include/phoenix_operator.hpp"
+#include "boost/spirit/include/phoenix_stl.hpp"
+#include "boost/smart_ptr.hpp"
 
+#include <fstream>
 #include <iostream>
 
 namespace local = baofit;
@@ -73,3 +82,153 @@ void local::AbsCorrelationData::_applyFinalCuts(std::set<int> &keep) const {
         keep.insert(index);
     }
 }
+
+likely::BinnedGrid local::createCorrelationGrid(std::string const &axis1Bins, std::string const &axis2Bins,
+    std::string const &axis3Bins, std::string const &axisLabels, bool verbose) {
+    // Extract the comma-separated 3 axis labels
+    std::vector<std::string> labels;
+    boost::split(labels,axisLabels,boost::is_any_of(","));
+    if(labels.size() != 3) {
+        throw RuntimeError("createCorrelationGrid: expected 3 axis labels.");
+    }
+    // Parse the binning from the strings provided for each axis
+    likely::AbsBinningCPtr axis1ptr,axis2ptr,axis3ptr;
+    try {
+        axis1ptr = likely::createBinning(axis1Bins);
+        if(verbose) {
+            std::cout << labels[0] << " bin centers:";
+            int nbins = axis1ptr->getNBins();
+            for(int bin = 0; bin < nbins; ++bin) {
+                std::cout << (bin ? ',':' ') << axis1ptr->getBinCenter(bin);
+            }
+            std::cout << " (n = " << nbins << ")" << std::endl;
+        }
+    }
+    catch(likely::RuntimeError const &e) {
+        throw RuntimeError("createCorrelationGrid: error in axis 1 (" + labels[0] + ") binning.");
+    }
+    try {
+        axis2ptr = likely::createBinning(axis2Bins);
+        if(verbose) {
+            std::cout << labels[1] << " bin centers:";
+            int nbins = axis2ptr->getNBins();
+            for(int bin = 0; bin < nbins; ++bin) {
+                std::cout << (bin ? ',':' ') << axis2ptr->getBinCenter(bin);
+            }
+            std::cout << " (n = " << nbins << ")" << std::endl;
+        }
+    }
+    catch(likely::RuntimeError const &e) {
+        throw RuntimeError("createCorrelationGrid: error in axis 2 (" + labels[1] + ") binning.");
+    }
+    try {
+        axis3ptr = likely::createBinning(axis3Bins);
+        if(verbose) {
+            std::cout << labels[2] << " bin centers:";
+            int nbins = axis3ptr->getNBins();
+            for(int bin = 0; bin < nbins; ++bin) {
+                std::cout << (bin ? ',':' ') << axis3ptr->getBinCenter(bin);
+            }
+            std::cout << " (n = " << nbins << ")" << std::endl;
+        }
+    }
+    catch(likely::RuntimeError const &e) {
+        throw RuntimeError("createCorrelationGrid: error in axis 3 (" + labels[2] + ") binning.");
+    }
+    // Return a BinnedGrid object for these 3 axes.
+    return likely::BinnedGrid(axis1ptr,axis2ptr,axis3ptr);    
+}
+
+namespace qi = boost::spirit::qi;
+namespace ascii = boost::spirit::ascii;
+namespace phoenix = boost::phoenix;
+
+baofit::AbsCorrelationDataPtr local::loadCorrelationData(std::string const &dataName,
+baofit::AbsCorrelationDataCPtr prototype, bool verbose, bool icov, bool weighted) {
+
+    // Create the new AbsCorrelationData that we will fill.
+    baofit::AbsCorrelationDataPtr binnedData(dynamic_cast<AbsCorrelationData*>(prototype->clone(true)));
+
+    // General stuff we will need for reading both files.
+    std::string line;
+    int lines;
+
+    // import boost spirit parser symbols
+    using qi::double_;
+    using qi::int_;
+    using qi::_1;
+    using phoenix::ref;
+    using phoenix::push_back;
+
+    // Loop over lines in the parameter file.
+    std::string paramsName = dataName + (weighted ? ".wdata" : ".data");
+    std::ifstream paramsIn(paramsName.c_str());
+    if(!paramsIn.good()) throw RuntimeError("loadCorrelationData: Unable to open " + paramsName);
+    lines = 0;
+    int index;
+    double data;
+    std::vector<double> bin(3);
+    while(std::getline(paramsIn,line)) {
+        lines++;
+        bin.resize(0);
+        bool ok = qi::phrase_parse(line.begin(),line.end(),
+            (
+                int_[ref(index) = _1] >> double_[ref(data) = _1]
+            ),
+            ascii::space);
+        if(!ok) {
+            throw RuntimeError("loadCorrelationData: error reading line " +
+                boost::lexical_cast<std::string>(lines) + " of " + paramsName);
+        }
+        binnedData->setData(index,data,weighted);
+    }
+    paramsIn.close();
+    int ndata = binnedData->getNBinsWithData();
+    int nbins = binnedData->getGrid().getNBinsTotal();
+    if(verbose) {
+        std::cout << "Read " << ndata << " of " << nbins << " data values from "
+            << paramsName << std::endl;
+    }
+
+    // Loop over lines in the (inverse) covariance file.
+    std::string covName = dataName + (icov ? ".icov" : ".cov");
+    std::ifstream covIn(covName.c_str());
+    if(!covIn.good()) throw RuntimeError("loadCorrelationData: Unable to open " + covName);
+    lines = 0;
+    double value;
+    int index1,index2;
+    while(std::getline(covIn,line)) {
+        lines++;
+        bin.resize(0);
+        bool ok = qi::phrase_parse(line.begin(),line.end(),
+            (
+                int_[ref(index1) = _1] >> int_[ref(index2) = _1] >> double_[ref(value) = _1]
+            ),
+            ascii::space);
+        if(!ok) {
+            throw RuntimeError("loadCorrelationData: error reading line " +
+                boost::lexical_cast<std::string>(lines) + " of " + covName);
+        }
+        // Check for invalid offsets.
+        if(index1 < 0 || index2 < 0 || index1 >= nbins || index2 >= nbins ||
+        !binnedData->hasData(index1) || !binnedData->hasData(index2)) {
+            throw RuntimeError("loadCorrelationData: invalid covariance indices on line " +
+                boost::lexical_cast<std::string>(lines) + " of " + covName);
+        }
+        // Add this covariance to our dataset.
+        if(icov) {
+            binnedData->setInverseCovariance(index1,index2,value);
+        }
+        else {
+            binnedData->setCovariance(index1,index2,value);            
+        }
+    }
+    covIn.close();
+    if(verbose) {
+        int ncov = (ndata*(ndata+1))/2;
+        std::cout << "Read " << lines << " of " << ncov
+            << " covariance values from " << covName << std::endl;
+    }
+
+    return binnedData;
+}
diff --git a/baofit/AbsCorrelationData.h b/baofit/AbsCorrelationData.h
@@ -3,6 +3,8 @@
 #ifndef BAOFIT_ABS_CORRELATION_DATA
 #define BAOFIT_ABS_CORRELATION_DATA
 
+#include "baofit/types.h"
+
 #include "likely/BinnedData.h"
 #include "likely/BinnedGrid.h"
 #include "likely/types.h"
@@ -58,6 +60,20 @@ namespace baofit {
 	inline AbsCorrelationData::TransverseBinningType
 	    AbsCorrelationData::getTransverseBinningType() const { return _type; }
 
+    // Creates and returns a BinnedGrid for a correlation function data set using the
+    // strings provided to specify the binning along each of the three axes. The
+    // axis labels provided as a comma-separated string (e.g., "r,mu,z") will be
+    // used for error messages and verbose printing (if requested).
+    likely::BinnedGrid createCorrelationGrid(std::string const &axis1Bins,
+        std::string const &axis2Bins, std::string const &axis3Bins,
+        std::string const &axisLabels, bool verbose);
+
+    // Loads a binned correlation function using the specified prototype
+    // and returns a BinnedData object. Set icov true to read .icov files instead of .cov.
+    // Set weighted true to read .wdata files instead of .data.
+    AbsCorrelationDataPtr loadCorrelationData(std::string const &dataName,
+        AbsCorrelationDataCPtr prototype, bool verbose, bool icov, bool weighted);
+
 } // baofit
 
 #endif // BAOFIT_ABS_CORRELATION_DATA
diff --git a/baofit/AbsCorrelationModel.cc b/baofit/AbsCorrelationModel.cc
@@ -20,7 +20,7 @@ local::AbsCorrelationModel::~AbsCorrelationModel() { }
 double local::AbsCorrelationModel::evaluate(double r, double mu, double z,
 likely::Parameters const &params) {
     bool anyChanged = updateParameterValues(params);
-    _applyVelocityShift(r,mu,z);
+    if(_crossCorrelation) _applyVelocityShift(r,mu,z);
     double result = _evaluate(r,mu,z,anyChanged);
     resetParameterValuesChanged();
     return result;

diff --git a/baofit/AbsCorrelationModel.h b/baofit/AbsCorrelationModel.h
@@ -45,6 +45,8 @@ namespace baofit {
         // of bias and beta parameters if crossCorrelation is true. Returns the index
         // of the last parameter defined.
         int _defineLinearBiasParameters(double zref, bool crossCorrelation = false);
+        // Returns the reference redshift passed to _defineLinearBiasParameters
+        double _getZRef() const;
         // Applies a shift dpi in the parallel direction to the separation (r,mu) using a fiducial
         // cosmology to convert from dv in km/s to dpi in Mpc/h at the specified redshift z. The
         // value of dv is obtained from the "delta-v" parameter defined by _defineLinearBiasParameters.
@@ -62,6 +64,9 @@ namespace baofit {
         };
         double _zref;
 	}; // AbsCorrelationModel
+
+    inline double AbsCorrelationModel::_getZRef() const { return _zref; }
+
 } // baofit
 
 #endif // BAOFIT_ABS_CORRELATION_MODEL
diff --git a/baofit/ComovingCorrelationData.cc b/baofit/ComovingCorrelationData.cc
@@ -8,8 +8,9 @@
 namespace local = baofit;
 
 local::ComovingCorrelationData::ComovingCorrelationData(likely::BinnedGrid grid,
-CoordinateSystem coordinateSystem) 
-: AbsCorrelationData(grid,Coordinate), _coordinateSystem(coordinateSystem), _lastIndex(-1)
+CoordinateSystem coordinateSystem)
+: AbsCorrelationData(grid,coordinateSystem==MultipoleCoordinates ? Multipole : Coordinate),
+_coordinateSystem(coordinateSystem), _lastIndex(-1)
 {
 }
 
@@ -38,24 +39,37 @@ void local::ComovingCorrelationData::_setIndex(int index) const {
 
 double local::ComovingCorrelationData::getRadius(int index) const {
     _setIndex(index);
-    if(_coordinateSystem == PolarCoordinates) {
-        return _binCenter[0];
-    }
-    else {
+    if(_coordinateSystem == CartesianCoordinates) {
         double rpar = _binCenter[0], rperp = _binCenter[1];
         return std::sqrt(rpar*rpar+rperp*rperp);
     }
+    else {
+        return _binCenter[0];
+    }
 }
 
 double local::ComovingCorrelationData::getCosAngle(int index) const {
     _setIndex(index);
     if(_coordinateSystem == PolarCoordinates) {
         return _binCenter[1];
     }
-    else {
+    else if(_coordinateSystem == CartesianCoordinates) {
         double rpar = _binCenter[0], rperp = _binCenter[1];
         return rpar/std::sqrt(rpar*rpar+rperp*rperp);
     }
+    else {
+        throw RuntimeError("ComovingCorrelationData::getMultipole: invalid coordinate system.");
+    }    
+}
+
+cosmo::Multipole local::ComovingCorrelationData::getMultipole(int index) const {
+    _setIndex(index);
+    if(_coordinateSystem == MultipoleCoordinates) {
+        return static_cast<cosmo::Multipole>(std::floor(_binCenter[1]+0.5));
+    }
+    else {
+        throw RuntimeError("ComovingCorrelationData::getMultipole: invalid coordinate system.");
+    }
 }
 
 double local::ComovingCorrelationData::getRedshift(int index) const {

diff --git a/baofit/ComovingCorrelationData.h b/baofit/ComovingCorrelationData.h
@@ -11,7 +11,7 @@ namespace baofit {
 	class ComovingCorrelationData : public AbsCorrelationData {
 	// Represents 3D correlation data in comoving coordinates (r,mu,z) or (rpar,rperp,z).
 	public:
-        enum CoordinateSystem { PolarCoordinates, CartesianCoordinates };
+        enum CoordinateSystem { PolarCoordinates, CartesianCoordinates, MultipoleCoordinates };
 		ComovingCorrelationData(likely::BinnedGrid grid,
 		    CoordinateSystem coordinateSystem = PolarCoordinates);
 		virtual ~ComovingCorrelationData();
@@ -22,6 +22,8 @@ namespace baofit {
         // Returns the cosine of the angle between the separation vector and
         // the line of sight (aka mu) associated with the specified global index.
         virtual double getCosAngle(int index) const;
+        // Returns the multipole (0,2,4) associated with the specified global index.
+        virtual cosmo::Multipole getMultipole(int index) const;
         // Returns the redshift associated with the specified global index.
         virtual double getRedshift(int index) const;
         // Finalize a comoving dataset by pruning to the limits specified in our constructor.