histogram class and associated plotting functions

examples/exRand.cpp

@@ -1,19 +1,25 @@
 #include <iostream>
 #include <LatAnalyze/AsciiFile.hpp>
+#include <LatAnalyze/CompiledFunction.hpp>
+#include <LatAnalyze/Plot.hpp>
 #include <LatAnalyze/RandGen.hpp>
 
 using namespace std;
 using namespace Latan;
 
-const int    seqLength     = 25;
+constexpr int    seqLength     = 25;
-const int    saveStep      = 9;
+constexpr int    saveStep      = 9;
-const string stateFileName = "exRand.seed";
+constexpr Index  nDraw         = 20000;
+const     string stateFileName = "exRand.seed";
 
 int main(void)
 {
     RandGenState state;
-    RandGen gen[2];
+    RandGen      gen[2];
-    AsciiFile stateFile(stateFileName, File::Mode::write|File::Mode::read);
+    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() <<endl;
     }
+    cout << "-- generating " << nDraw << " Gaussian random numbers..." << endl;
+    FOR_VEC(gauss, i)
+    {
+        gauss(i) = gen[0].gaussian();
+    }
+    h.setFromData(gauss, -5., 5., 40);
+    h.normalize();
+    p << PlotHistogram(h);
+    p << PlotFunction(compile("return exp(-x_0^2/2)/sqrt(2*pi);", 1), -5., 5.);
+    p.display();
+
     return EXIT_SUCCESS;
 }

lib/Histogram.cpp

@@ -0,0 +1,141 @@
+/*
+ * Histogram.cpp, 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 <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Histogram.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+#define DECL_GSL_HIST(h) \
+gsl_histogram h{static_cast<size_t>(bin_.size()), x_.data(), bin_.data()}
+#define DECL_CONST_GSL_HIST(h) \
+const gsl_histogram h{static_cast<size_t>(bin_.size()),\
+                      const_cast<double *>(x_.data()),\
+                      const_cast<double *>(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<double>(data.size());
+    computeNorm();
+}
+
+// compute normalization factor ////////////////////////////////////////////////
+void Histogram::computeNorm(void)
+{
+    norm_ = static_cast<double>(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<Index>(i)];
+}

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 <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_Histogram_hpp_
+#define Latan_Histogram_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <gsl/gsl_histogram.h>
+#include <gsl/gsl_cdf.h>
+
+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_

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                \

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                                 *
  ******************************************************************************/

lib/Plot.hpp

@@ -22,6 +22,7 @@
 
 #include <LatAnalyze/Global.hpp>
 #include <LatAnalyze/Mat.hpp>
+#include <LatAnalyze/Histogram.hpp>
 #include <LatAnalyze/XYStatData.hpp>
 #include <sstream>
 #include <vector>
@@ -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                                 *
  ******************************************************************************/