diff --git a/examples/exRand.cpp b/examples/exRand.cpp index 974f44d..c88e79d 100644 --- a/examples/exRand.cpp +++ b/examples/exRand.cpp @@ -1,19 +1,25 @@ #include #include +#include +#include #include using namespace std; using namespace Latan; -const int seqLength = 25; -const int saveStep = 9; -const string stateFileName = "exRand.seed"; +constexpr int seqLength = 25; +constexpr int saveStep = 9; +constexpr Index nDraw = 20000; +const string stateFileName = "exRand.seed"; int main(void) { RandGenState state; - RandGen gen[2]; - AsciiFile stateFile(stateFileName, File::Mode::write|File::Mode::read); + RandGen gen[2]; + AsciiFile stateFile(stateFileName, File::Mode::write|File::Mode::read); + DVec gauss(nDraw); + Plot p; + Histogram h; cout << "- GENERATOR STATE I/O TESTS" << endl; cout << "-- generating a " << seqLength << " steps random sequence..." @@ -38,5 +44,16 @@ int main(void) { cout << "step " << i << "\t: " << gen[1].uniform() <. + */ + +#include +#include + +using namespace std; +using namespace Latan; + +#define DECL_GSL_HIST(h) \ +gsl_histogram h{static_cast(bin_.size()), x_.data(), bin_.data()} +#define DECL_CONST_GSL_HIST(h) \ +const gsl_histogram h{static_cast(bin_.size()),\ + const_cast(x_.data()),\ + const_cast(bin_.data())} + +/****************************************************************************** + * Histogram implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +Histogram::Histogram(const DVec &data, const double xMin, const double xMax, + const Index nBin) +: Histogram() +{ + setFromData(data, xMin, xMax, nBin); +} + +Histogram::Histogram(const DVec &data, const DVec &w, const double xMin, + const double xMax, const Index nBin) +: Histogram() +{ + setFromData(data, w, xMin, xMax, nBin); +} + +// resize ////////////////////////////////////////////////////////////////////// +void Histogram::resize(const Index nBin) +{ + x_.resize(nBin + 1); + bin_.resize(nBin); +} + +// generate from data ////////////////////////////////////////////////////////// +void Histogram::setFromData(const DVec &data, const DVec &w, const double xMin, + const double xMax, const Index nBin) +{ + if (data.size() != w.size()) + { + LATAN_ERROR(Size, "data vector and weight vector size mismatch"); + } + resize(nBin); + + DECL_GSL_HIST(h); + + xMax_ = xMax; + xMin_ = xMin; + gsl_histogram_set_ranges_uniform(&h, xMin, xMax); + FOR_VEC(data, i) + { + gsl_histogram_accumulate(&h, data(i), w(i)); + } + total_ = w.sum(); + computeNorm(); +} + +void Histogram::setFromData(const DVec &data, const double xMin, + const double xMax, const Index nBin) +{ + resize(nBin); + + DECL_GSL_HIST(h); + + xMax_ = xMax; + xMin_ = xMin; + gsl_histogram_set_ranges_uniform(&h, xMin, xMax); + FOR_VEC(data, i) + { + gsl_histogram_increment(&h, data(i)); + } + total_ = static_cast(data.size()); + computeNorm(); +} + +// compute normalization factor //////////////////////////////////////////////// +void Histogram::computeNorm(void) +{ + norm_ = static_cast(bin_.size())/(total_*(xMax_ - xMin_)); +} + +// normalize as a probablility ///////////////////////////////////////////////// +void Histogram::normalize(const bool n) +{ + normalize_ = n; +} + +bool Histogram::isNormalized(void) const +{ + return normalize_; +} + +// access ////////////////////////////////////////////////////////////////////// +Index Histogram::size(void) const +{ + return bin_.size(); +} + +double Histogram::getX(const Index i) const +{ + return x_(i); +} + +double Histogram::operator[](const Index i) const +{ + return bin_(i)*(isNormalized() ? norm_ : 1.); +} + +double Histogram::operator()(const double x) const +{ + size_t i; + + DECL_CONST_GSL_HIST(h); + + gsl_histogram_find(&h, x, &i); + + return (*this)[static_cast(i)]; +} diff --git a/lib/Histogram.hpp b/lib/Histogram.hpp new file mode 100644 index 0000000..d5b1269 --- /dev/null +++ b/lib/Histogram.hpp @@ -0,0 +1,69 @@ +/* + * Histogram.hpp, part of LatAnalyze 3 + * + * Copyright (C) 2013 - 2015 Antonin Portelli + * + * LatAnalyze 3 is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * LatAnalyze 3 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with LatAnalyze 3. If not, see . + */ + +#ifndef Latan_Histogram_hpp_ +#define Latan_Histogram_hpp_ + +#include +#include +#include + +BEGIN_LATAN_NAMESPACE + +/****************************************************************************** + * Histogram class * + ******************************************************************************/ +class Histogram +{ +public: + // constructor + Histogram(void) = default; + Histogram(const DVec &data, const double xMin, const double xMax, + const Index nBin); + Histogram(const DVec &data, const DVec &w, const double xMin, + const double xMax, const Index nBin); + // destructor + virtual ~Histogram(void) = default; + // generate from data + void setFromData(const DVec &data, const double xMin, const double xMax, + const Index nBin); + void setFromData(const DVec &data, const DVec &w, const double xMin, + const double xMax, const Index nBin); + // normalize as a probablility + void normalize(const bool n = true); + bool isNormalized(void) const; + // access + Index size(void) const; + double getX(const Index i) const; + double operator[](const Index i) const; + double operator()(const double x) const; +private: + // resize + void resize(const Index nBin); + // compute normalization factor + void computeNorm(void); +private: + DVec x_, bin_; + double total_, norm_, xMax_, xMin_; + bool normalize_{false}; +}; + +END_LATAN_NAMESPACE + +#endif // Latan_Histogram_hpp_ diff --git a/lib/Makefile.am b/lib/Makefile.am index 7454a3d..0df21a3 100644 --- a/lib/Makefile.am +++ b/lib/Makefile.am @@ -29,6 +29,7 @@ libLatAnalyze_la_SOURCES = \ Global.cpp \ GslHybridRootFinder.cpp\ GslQagsIntegrator.cpp \ + Histogram.cpp \ includes.hpp \ Mat.cpp \ Math.cpp \ @@ -61,6 +62,7 @@ libLatAnalyze_la_HEADERS = \ Global.hpp \ GslHybridRootFinder.hpp\ GslQagsIntegrator.hpp \ + Histogram.hpp \ Integrator.hpp \ IoObject.hpp \ Mat.hpp \ diff --git a/lib/Plot.cpp b/lib/Plot.cpp index d65d6c9..eb4f79a 100644 --- a/lib/Plot.cpp +++ b/lib/Plot.cpp @@ -236,6 +236,22 @@ PlotPredBand::PlotPredBand(const DoubleFunctionSample &function, " fs solid " + strFrom(opacity) + " noborder"); } +// PlotHistogram constructor /////////////////////////////////////////////////// +PlotHistogram::PlotHistogram(const Histogram &h) +{ + DMat d(h.size(), 2); + string tmpFileName; + + for (Index i = 0; i < h.size(); ++i) + { + d(i, 0) = h.getX(i); + d(i, 1) = h[i]; + } + tmpFileName = dumpToTmpFile(d); + pushTmpFile(tmpFileName); + setCommand("'" + tmpFileName + "' u 1:2 w steps"); +} + /****************************************************************************** * Plot modifiers * ******************************************************************************/ diff --git a/lib/Plot.hpp b/lib/Plot.hpp index 096ca0a..57ffe65 100644 --- a/lib/Plot.hpp +++ b/lib/Plot.hpp @@ -22,6 +22,7 @@ #include #include +#include #include #include #include @@ -126,6 +127,15 @@ public: virtual ~PlotPredBand(void) = default; }; +class PlotHistogram: public PlotObject +{ +public: + // constructor + PlotHistogram(const Histogram &h); + // destructor + virtual ~PlotHistogram(void) = default; +}; + /****************************************************************************** * Plot modifiers * ******************************************************************************/