Page MenuHomeHEPForge
Files F11222235

No OneTemporary
Actions

View Options

diff --git a/ChangeLog b/ChangeLog
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,102 +1,106 @@
2013-01-13 Andy Buckley <andy.buckley@cern.ch>
+ * Using std::auto_ptr to handle Interpolator/Extrapolator
+ assignment to PDFGrid, and cleaning up the interface and
+ implementation a bit.
+
* Using boost::shared_ptr to implement the memory handling in
LHAGLUE. NB. not std::auto_ptr, or boost::scoped_ptr, and although
std::unique_ptr would be ideal it requires C++11.
* Adding alpha_s(Q2) calculation mechanism to the PDF
interface. AlphaS calculation is not yet properly implemented and
tested, though, cf. grid interpolation!
* Adding src/Factories.cc to hide the implementation details /
dependencies for factory users.
* Added Info::metadata<bool> template specialisation to handle
true/false/on/off/yes/no as well as 0/1 strings.
* Add src/Paths.cc to avoid circular header dependencies.
* Add LHAPDF_{INSTALL,DATA}_PREFIX variables and PwdInSearchPath
config flag, and use them in search path determination.
2013-01-11 Andy Buckley <andy.buckley@cern.ch>
* Start of the AlphaS refactoring.
* Metadata methods and flavour list caching added to PDF.
* Renaming PDFInfo.h/.cc to just Info.h/.cc
* Renaming LHAPDFConfig.h -> Config.h
* DESIGN doc updates.
* Factory renaming and adding factory functions for PDF construction.
* PDF and PDFInfo extra constructors from LHAPDF ID code.
* Moving PDFInfo YAML reading into a new PDFInfo.cc to avoid an API header dependency.
* Refactoring PDF filename construction etc. to minimise code duplication.
* Adding the PDFIndex.h header and index lookup functions.
* Adding the start of the LHA Fortran wrapper, based on the Py8
wrapper by Steve and Martin.
2013-01-10 Andy Buckley <andy.buckley@cern.ch>
* More PDFGrid data loading development: improvements to the API
and the data parser now works.
* Adding the error extrapolator to the extrapolator factory function.
2013-01-09 Andy Buckley <andy.buckley@cern.ch>
* A few tweaks to the info/data loader system: PDF loading from
set name + member ID with cascading info/set/config levels WORKS!
* Adding more constructors for PDFGrid, calling the new ones on
PDF to populate the Info class.
* Adding the EXAMPLEPDF directory for testing.
2013-01-08 Andy Buckley <andy.buckley@cern.ch>
* Adding example lhapdf.conf.
* Inlining, clean-up and other tweaks.
* Adding PDF base constructors with loading apparatus for info
from member files + discovery and loading of set-level info.
2013-01-07 Andy Buckley <andy.buckley@cern.ch>
* Adding to_str variants for both single objects and for vectors
of objects which are convertible to string via lexical_cast.
* Adding the YAML parser operation to the Info class (and its subclasses).
* Removing Types.h in favour of split Utils.h and Exceptions.h
* include/LHAPDF/Paths.h: Adding path searching machinery (using Boost.Filesystem).
2013-01-02 Andy Buckley <andy.buckley@cern.ch>
* Adding some of the necessary infrastructure for subgrids, and
rewriting the grid format parser to handle separated subgrid
blocks. Boost.MultiArray being used for the internal storage in
place of the dynamically allocated C arrays. Metadata methods are
currently disabled for porting to the cascading Info system.
2013-01-01 Andy Buckley <andy.buckley@cern.ch>
* Major restructuring of the API: lots of subtle problems with the
first attempt, and room for improvement. Painful, but worth the
reworking to get it right. Lots left to do... see the TODO and
DESIGN files for the ideas and tasks.
2012-10-22 Andy Buckley <andy.buckley@cern.ch>
* Start of autotools build system on top of initial summer student
work by Martin and Steve.
diff --git a/include/LHAPDF/PDFGrid.h b/include/LHAPDF/PDFGrid.h
--- a/include/LHAPDF/PDFGrid.h
+++ b/include/LHAPDF/PDFGrid.h
@@ -1,289 +1,304 @@
#pragma once
#include "LHAPDF/PDF.h"
#include "LHAPDF/Interpolator.h"
#include "LHAPDF/Extrapolator.h"
#include "boost/multi_array.hpp"
namespace LHAPDF {
/// @brief A PDF defined via an interpolation grid
class PDFGrid : public PDF {
public:
/// @name Creation and deletion
//@{
// /// Default constructor.
// PDFGrid()
// : _interpolator(0), _extrapolator(0)
// { }
/// Constructor from a file path.
PDFGrid(const std::string& path)
- : PDF(path), _interpolator(0), _extrapolator(0)
+ : PDF(path)
{
_loadData(_mempath);
_init();
}
/// Constructor from a set name and member ID.
PDFGrid(const std::string& setname, int member)
- : PDF(setname, member), _interpolator(0), _extrapolator(0)
+ : PDF(setname, member)
{
_loadData(_mempath);
_init();
}
/// Constructor from a set name and member ID.
PDFGrid(int lhaid)
- : PDF(lhaid), _interpolator(0), _extrapolator(0)
+ : PDF(lhaid)
{
_loadData(_mempath);
_init();
}
- /// Destructor
- ~PDFGrid() {
- delete _interpolator;
- delete _extrapolator;
- }
+ /// Virtual destructor to allow inheritance.
+ virtual ~PDFGrid() { }
//@}
protected:
+ /// Load the PDF grid data block (not the metadata) from the given PDF member file
void _loadData(const path& mempath);
+ /// Set default inter/extrapolators, etc.
void _init();
public:
/// Metadata
//@{
/// Get the list of available flavours by PDG ID code.
/// @todo Or get the flavour list from the set?
std::vector<int> flavors() const {
std::vector<int> rtn;
for (std::map<int, double*>::const_iterator i = _ptdata.begin(); i != _ptdata.end(); ++i) {
rtn.push_back(i->first);
}
return rtn;
}
/// Check if x is in the grid range
bool inRangeX(double x) const {
if (x < xKnots().front()) return false;
if (x > xKnots().back()) return false;
return true;
}
/// Check if q2 is in the grid range
bool inRangeQ2(double q2) const {
if (q2 < q2Knots().front()) return false;
if (q2 > q2Knots().back()) return false;
return true;
}
//@}
/// @name Interpolators and extrapolators
//@{
+ /// @brief Set the interpolator by pointer
+ ///
+ /// The provided Interpolator must have been new'd, as it will not be copied
+ /// and ownership passes to this PDFGrid: delete will be called on this ptr
+ /// when this PDFGrid goes out of scope or another setInterpolator call is made.
+ void setInterpolator(Interpolator* ipol) {
+ _interpolator.reset(ipol);
+ _interpolator->bind(this);
+ }
+
/// @brief Set the interpolator by value
///
/// The passed value must be a concrete instantiation of the Interpolator
/// interface. It will be copied and heap-assigned for use inside this PDFGrid.
template <typename INTERPOLATOR>
void setInterpolator(INTERPOLATOR ipol) {
- _interpolator = new INTERPOLATOR(ipol);
- _interpolator->bind(this);
- }
-
- /// @brief Set the interpolator by pointer
- ///
- /// This interpolator argument is exactly the one that will be used by this
- /// PDFGrid: it will not be copied and the PDF takes ownership of the
- /// pointer and will delete the interpolator when the PDF goes out of scope.
- ///
- /// @todo Use smart pointers?
- void setInterpolator(Interpolator* ipol) {
- _interpolator = ipol;
- _interpolator->bind(this);
+ setInterpolator(new INTERPOLATOR(ipol));
}
/// @brief Set the interpolator by name
///
/// Use the interpolator specified by the given name, as passed to the
/// createInterpolator factory function.
- void setInterpolator(const std::string& ipolname);
-
- /// Get the current interpolator (ownership remains with the PDFGrid).
- const Interpolator* interpolator() const {
- return _interpolator;
+ void setInterpolator(const std::string& ipolname) {
+ setInterpolator(mkInterpolator(ipolname));
}
+ /// Find whether an extrapolator has been set on this PDF
+ bool hasInterpolator() const {
+ return _interpolator.get() != 0;
+ }
+ /// Get the current interpolator
+ const Interpolator& interpolator() const {
+ if (_interpolator.get() == 0)
+ throw GridError("No interpolator has been set on this PDFGrid");
+ return *_interpolator;
+ }
+
+
+
+ /// @brief Set the extrapolator by pointer
+ ///
+ /// The provided Extrapolator must have been new'd, as it will not be copied
+ /// and ownership passes to this PDFGrid: delete will be called on this ptr
+ /// when this PDFGrid goes out of scope or another setExtrapolator call is made.
+ void setExtrapolator(Extrapolator* xpol) {
+ _extrapolator.reset(xpol);
+ _extrapolator->bind(this);
+ }
/// @brief Set the extrapolator by value
///
/// The passed value must be a concrete instantiation of the Extrapolator
/// interface. It will be copied and heap-assigned for use inside this PDFGrid.
template <typename EXTRAPOLATOR>
void setExtrapolator(EXTRAPOLATOR xpol) {
- _extrapolator = new EXTRAPOLATOR(xpol);
- _extrapolator->bind(this);
- }
-
- /// @brief Set the extrapolator by pointer
- ///
- /// This extrapolator argument is exactly the one that will be used by this
- /// PDFGrid: it will not be copied and the PDF takes ownership of the
- /// pointer and will delete the extrapolator when the PDF goes out of scope.
- ///
- /// @todo Use smart pointers?
- void setExtrapolator(Extrapolator* xpol) {
- _extrapolator = xpol;
- _extrapolator->bind(this);
+ setExtrapolator(new EXTRAPOLATOR(xpol));
}
/// @brief Set the extrapolator by name
///
/// Use the extrapolator specified by the given name, as passed to the
/// createExtrapolator factory function.
- void setExtrapolator(const std::string& xpolname);
+ void setExtrapolator(const std::string& xpolname) {
+ setExtrapolator(mkExtrapolator(xpolname));
+ }
- /// Get the current extrapolator (ownership remains with the PDFGrid).
- const Extrapolator* extrapolator() const {
- return _extrapolator;
+ /// Find whether an extrapolator has been set on this PDF
+ bool hasExtrapolator() const {
+ return _extrapolator.get() != 0;
+ }
+
+ /// Get the current extrapolator
+ const Extrapolator& extrapolator() const {
+ if (_extrapolator.get() == 0)
+ throw GridError("No extrapolator has been set on this PDFGrid");
+ return *_extrapolator;
}
//@}
/// @name Info about the grid, and access to the raw data points
//@{
/// Return knot values in x
const std::vector<double>& xKnots() const {
return _xknots;
}
/// Return knot values in Q2
const std::vector<double>& q2Knots() const {
return _q2knots;
}
/// Get the index of the closest x knot row <= x
size_t xKnotLow(double x) const {
/// @todo Test for x in grid range
size_t i = lower_bound(xKnots().begin(), xKnots().end(), x) - xKnots().begin();
if (i == xKnots().size()-1) --i; // if last row, step back
return i;
}
/// Get the index of the closest Q2 knot column <= q2
size_t q2KnotLow(double q2) const {
size_t i = lower_bound(q2Knots().begin(), q2Knots().end(), q2) - q2Knots().begin();
if (i == q2Knots().size()-1) --i; // if last col, step back
return i;
}
/// Get the raw xf(x,Q2) data points
const double* ptdata(int id) const {
if (!hasFlavor(id)) {
std::stringstream error;
error << "Undefined particle ID requested: " << id;
throw FlavorError(error.str());
}
return _ptdata.find(id)->second;
}
/// @brief Transform a (ix, iQ2) pair into a 1D "raw" index
size_t ptindex(size_t ix, size_t iq2) const {
if (ix >= xKnots().size()) throw GridError("Invalid x index");
if (iq2 >= q2Knots().size()) throw GridError("Invalid Q2 index");
return ix + iq2 * xKnots().size();
}
//@}
protected:
/// @brief Get PDF xf(x,Q2) value
double _xfxQ2(int, double x, double q2) const;
/// @name Internal storage
//@{
/// We use "array" to refer to the "raw" knot grid, while "grid" means a grid-based PDF.
/// The "1F" means that this is a single-flavour array
class KnotArray1F {
public:
// Use the Boost multi_array for efficiency and ease of indexing
typedef boost::multi_array<double, 2> valarray;
KnotArray1F() {} //< for std::map insertability
KnotArray1F(const std::vector<double>& xknots, const std::vector<double>& q2knots, const valarray& xfs)
: _xs(xknots), _q2s(q2knots), _xfs(xfs)
{ assert(_xfs.shape()[0] == xknots.size() && _xfs.shape()[1] == q2knots.size()); }
KnotArray1F(const std::vector<double>& xknots, const std::vector<double>& q2knots)
: _xs(xknots), _q2s(q2knots), _xfs(boost::extents[xknots.size()][q2knots.size()])
{ assert(_xfs.shape()[0] == xknots.size() && _xfs.shape()[1] == q2knots.size()); }
const std::vector<double>& xs() const { return _xs; }
void setxs(const std::vector<double>& xs) { _xs = xs; _xfs.resize(boost::extents[_xs.size()][_q2s.size()]); }
const std::vector<double>& q2s() const { return _q2s; }
void setq2s(const std::vector<double>& q2s) { _q2s = q2s; _xfs.resize(boost::extents[_xs.size()][_q2s.size()]); }
const valarray& xfs() const { return _xfs; }
valarray& xfs() { return _xfs; }
void setxfs(const valarray& xfs) { _xfs = xfs; }
/// @todo Add index converter and finder methods
private:
std::vector<double> _xs, _q2s;
valarray _xfs;
};
/// Typedef for a collection of KnotArray1F accessed by PID code
/// The "NF" means "> 1 flavour", cf. the KnotArray1F name for a single flavour data array.
typedef std::map<int, KnotArray1F> KnotArrayNF;
//@}
private:
/// Interpolation grid anchor point lists in x and Q2
std::vector<double> _xknots, _q2knots;
/// Raw data grids, indexed by flavour
/// @todo Need an intermediate type for the subgrids
std::map<int, double*> _ptdata;
/// Map of multi-flavour KnotArrays "binned" for lookup by low edge in Q2
std::map<double, KnotArrayNF> _knotarrays;
+ /// Typedefs of smart pointer types for ipol/xpol memory handling
+ typedef auto_ptr<Interpolator> InterpolatorPtr;
+ typedef auto_ptr<Extrapolator> ExtrapolatorPtr;
+
/// Associated interpolator
- Interpolator* _interpolator;
+ InterpolatorPtr _interpolator;
/// Associated extrapolator
- Extrapolator* _extrapolator;
+ ExtrapolatorPtr _extrapolator;
};
}
diff --git a/src/PDFGrid.cc b/src/PDFGrid.cc
--- a/src/PDFGrid.cc
+++ b/src/PDFGrid.cc
@@ -1,139 +1,125 @@
#include "LHAPDF/PDFGrid.h"
#include "LHAPDF/Interpolator.h"
#include "LHAPDF/Factories.h"
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <stdexcept>
#include <cstring>
using namespace std;
namespace LHAPDF {
double PDFGrid::_xfxQ2(int id, double x, double q2 ) const {
// Decide whether to use interpolation or extrapolation... the sanity checks
// are done in the public PDF::xfxQ2 function.
if (inRangeXQ2(x, q2)) {
- if (interpolator() == 0) throw GridError("Undefined interpolator");
- return interpolator()->interpolateXQ2(id, x, q2);
+ return interpolator().interpolateXQ2(id, x, q2);
} else {
- if (extrapolator() == 0) throw GridError("Undefined extrapolator");
- return extrapolator()->extrapolateXQ2(id, x, q2);
+ return extrapolator().extrapolateXQ2(id, x, q2);
}
}
void PDFGrid::_loadData(const path& mempath) {
string line;
int iblock(0), iblockline(0), iline(0);
vector<double> xs, q2s;
const size_t npid = flavors().size();
vector< vector<double> > ipid_xfs(npid);
try {
ifstream file(mempath.c_str());
while (getline(file, line)) {
iline += 1;
iblockline += 1;
if (line != "---") { // if we are not on a block separator line...
// Block 0 is the metadata, which we ignore here
if (iblock == 0) continue;
/// @todo Check that there is no effect of leading spaces, etc.: itrim the lines
// Parse the data lines
double token;
istringstream tokens(line);
if (iblockline == 1) { // x knots line
while (tokens >> token) xs.push_back(token);
} if (iblockline == 2) { // Q2 knots line
while (tokens >> token) q2s.push_back(token);
} else {
if (iblockline == 3) { // on the first line of the xf block, resize the arrays
for (size_t ipid = 0; ipid < npid; ++ipid) { ipid_xfs[ipid].reserve(xs.size() * q2s.size()); }
}
size_t ipid = 0;
while (tokens >> token) {
ipid_xfs[ipid].push_back(token);
ipid += 1;
}
// Check that each line has many tokens as there should be flavours
if (ipid != flavors().size())
throw ReadError("PDF grid data error on line " + to_str(iline) + ": " + to_str(ipid) +
" flavor entries seen but " + to_str(flavors().size()) + " were expected");
}
} else { // we *are* on a block separator line
// Check that the expected number of data lines were seen in the last block
if (iblock > 0 && iblockline - 1 != int(xs.size()*q2s.size()) + 2)
throw ReadError("PDF grid data error on line " + to_str(iline) + ": " +
to_str(iblockline-1) + " data lines were seen in block " + to_str(iblock-1) +
" but " + to_str(xs.size()*q2s.size() + 2) + " were expected");
// Increment/reset the block and line counters
iblock += 1;
iblockline = 0;
// Escape here if we've just finished reading the 0th (metadata) block
if (iblock == 1) continue;
// Die with an assert if the block was of zero size
/// @todo Convert to throwing some exception? Is this ever allowable?
assert(xs.size() > 0);
assert(q2s.size() > 0);
assert(ipid_xfs.size() > 0);
/// @todo Define the ordering of the values in x and Q2 indexing
// Register data from the previous (>0th) block into the PDFGrid data structure
// KnotArrayNF arraynf;
KnotArrayNF& arraynf = _knotarrays[q2s.front()]; //< Reference to newly created subgrid object
for (size_t ipid = 0; ipid < npid; ++ipid) {
int pid = flavors()[ipid];
arraynf[pid] = KnotArray1F(xs, q2s); // create the 2D array with the x and Q2 knot positions
arraynf[pid].xfs().assign(ipid_xfs[ipid].begin(), ipid_xfs[ipid].end()); // populate the xf array
}
cout << q2s.size() << endl;
cout << q2s.front() << endl;
// _knotarrays[q2s.front()] = arraynf; //< @todo Prefer getting the reference as above, to avoid re-copying a big data array
xs.clear(); q2s.clear(); ipid_xfs.clear();
}
}
} catch (Exception& e) {
throw;
} catch (std::exception& e) {
throw ReadError("Read error while parsing " + mempath.native() + " as a PDFGrid data file");
}
}
void PDFGrid::_init() {
// Set default inter/extrapolators
const string ipolname = info().metadata("Interpolator");
setInterpolator(ipolname);
const string xpolname = info().metadata("Extrapolator");
setExtrapolator(xpolname);
}
- // Defined here to avoid circular dependencies in the headers
- void PDFGrid::setInterpolator(const std::string& ipolname) {
- setInterpolator(mkInterpolator(ipolname));
- }
-
-
- // Defined here to avoid circular dependencies in the headers
- void PDFGrid::setExtrapolator(const std::string& xpolname) {
- setExtrapolator(mkExtrapolator(xpolname));
- }
-
-
}

File Metadata

Mime Type
text/x-diff
Expires
Wed, May 14, 11:38 AM (14 h, 25 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111471
Default Alt Text
(20 KB)

Event Timeline