source: branches/0.7-stable/yat/statistics/Histogram.cc @ 2405

Last change on this file since 2405 was 2405, checked in by Peter, 11 years ago

help compiler choose correct load function. fixes #649

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.0 KB
Line 
1// $Id: Histogram.cc 2405 2011-01-12 14:56:30Z peter $
2
3/*
4  Copyright (C) 2004 Jari Häkkinen
5  Copyright (C) 2005 Peter Johansson
6  Copyright (C) 2006 Jari Häkkinen
7  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
8  Copyright (C) 2009, 2010 Peter Johansson
9
10  This file is part of the yat library, http://dev.thep.lu.se/yat
11
12  The yat library is free software; you can redistribute it and/or
13  modify it under the terms of the GNU General Public License as
14  published by the Free Software Foundation; either version 3 of the
15  License, or (at your option) any later version.
16
17  The yat library is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20  General Public License for more details.
21
22  You should have received a copy of the GNU General Public License
23  along with yat. If not, see <http://www.gnu.org/licenses/>.
24*/
25
26#include "Histogram.h"
27
28#include <yat/utility/utility.h>
29
30#include <algorithm>
31#include <cassert>
32#include <cmath>
33#include <functional>
34#include <istream>
35#include <ostream>
36
37namespace theplu {
38namespace yat {
39namespace statistics {
40
41
42  Histogram::Histogram(void)
43    : xmax_(0), xmin_(0), sum_all_(), sum_histogram_()
44  {
45  }
46
47
48  Histogram::Histogram(std::istream& is)
49  {
50    std::string line;
51    getline(is, line, ':');
52    getline(is, line, '\n');
53    xmin_ = utility::convert<double>(line);
54    getline(is, line, ':');
55    getline(is, line, '\n');
56    xmax_ = utility::convert<double>(line);
57    getline(is, line, ':');
58    getline(is, line, '\n');
59    size_t n = utility::convert<size_t>(line);
60    histogram_.resize(n);
61    getline(is, line, '\n');
62    getline(is, line, ':');
63    getline(is, line, '\n');
64    double total = utility::convert<double>(line);
65    getline(is, line, '\n');
66    getline(is, line, '\n');
67    std::vector<std::vector<double> > data;
68    utility::load(is, data, '\0', '\n');
69    for (size_t i=0; i<histogram_.size(); ++i) {
70      assert(data[i].size()==2);
71      add(data[i][0], data[i][1]);
72    }
73    add(xmax_, total - averager_all().sum_w());
74  }
75
76
77  Histogram::Histogram(const Histogram& b)
78  {
79    *this=b;
80  }
81
82
83  Histogram::Histogram(const double min, const double max, const size_t n)
84    : histogram_(std::vector<double>(n,0.0)),
85      xmax_(max), xmin_(min),
86      sum_all_(), sum_histogram_()
87  {
88  }
89
90
91  Histogram::~Histogram(void)
92  {
93  }
94
95
96  int Histogram::add(const double x, const double w)
97  {
98    sum_all_.add(x,w);
99    if (x<xmin_)
100      return -1;
101    else if (x>=xmax_)
102      return 1;
103
104    sum_histogram_.add(x,w);
105    histogram_[bin(x)] += w;
106    return 0;
107  }
108
109
110  const AveragerWeighted& Histogram::averager_all(void) const
111  {
112    return sum_all_;
113  }
114
115
116  const AveragerWeighted& Histogram::averager_histogram(void) const
117  {
118    return sum_histogram_;
119  }
120
121
122  size_t Histogram::bin(double d)
123  {
124    return (((d<xmin_) || (d>xmax_)) ? 0 :
125            static_cast<size_t>(floor((d-xmin_)/spacing() )));
126  }
127
128
129  size_t Histogram::nof_bins(void) const
130  {
131    return histogram_.size();
132  }
133
134
135  void Histogram::normalize(bool choice)
136  {
137    if (choice) 
138      rescale(1.0/sum_all_.sum_w());
139    else
140      rescale(1.0/sum_all_.sum_w()/spacing());
141  }
142
143
144  double Histogram::observation_value(const size_t k) const
145  {
146    return xmin_+spacing()*(k+0.5);
147  }
148
149
150  void Histogram::rescale(double factor)
151  {
152    for (size_t i=0; i<histogram_.size(); i++)
153      histogram_[i]*=factor;
154  }
155
156
157  void Histogram::reset(void)
158  {
159    for (size_t i=0; i<histogram_.size(); i++)
160      histogram_[i]=0;
161    sum_all_.reset();
162    sum_histogram_.reset();
163  }
164
165
166  double Histogram::spacing(void) const
167  {
168    return (xmax_-xmin_)/nof_bins();
169  }
170
171
172  double Histogram::xmax(void) const
173  {
174    return xmax_;
175  }
176
177
178  double Histogram::xmin(void) const
179  {
180    return xmin_;
181  }
182
183
184  double Histogram::operator[](size_t k) const
185  {
186    return histogram_[k];
187  }
188
189
190  const Histogram& Histogram::operator=(const Histogram& b)
191  {
192    if (this==&b)
193      return *this;
194    histogram_=b.histogram_;
195    xmax_=b.xmax_;
196    xmin_=b.xmin_;
197    sum_all_=b.sum_all_;
198    sum_histogram_=b.sum_histogram_;
199    return *this;
200  }
201
202
203  const Histogram& Histogram::operator+=(const Histogram& rhs)
204  {
205    assert(histogram_.size()==rhs.histogram_.size());
206    std::transform(histogram_.begin(), histogram_.end(),
207                   rhs.histogram_.begin(), histogram_.begin(),
208                   std::plus<double>());
209    assert(xmax_==rhs.xmax_);
210    assert(xmin_==rhs.xmin_);
211    sum_all_ += rhs.sum_all_;
212    sum_histogram_ += rhs.sum_histogram_;
213    return *this;
214  }
215
216
217  std::ostream& operator<<(std::ostream& s,const Histogram& histogram)
218  {
219    s << "# histogram min : " << histogram.xmin() << '\n';
220    s << "# histogram max : " << histogram.xmax() << '\n';
221    s << "# number of bins: " << histogram.nof_bins() << '\n';
222    s << "# nof points in histogram : " 
223      << histogram.averager_histogram().sum_w() << '\n';
224    s << "# nof points in total:      " 
225      << histogram.averager_all().sum_w() << '\n';
226    s << "# column 1: center of observation bin\n"
227      << "# column 2: frequency\n";
228
229    for (size_t i=0; i<histogram.nof_bins(); i++) {
230      s.width(12);
231      s << histogram.observation_value(i);
232      s.width(12);
233      s << histogram[i] << '\n';
234    }
235
236    return s;
237  }
238
239}}} // of namespace statistics, yat, and theplu
Note: See TracBrowser for help on using the repository browser.