source: src/Fragmentation/SetValues/Histogram.cpp@ beb16e

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since beb16e was beb16e, checked in by Frederik Heber <heber@…>, 12 years ago

Added template function ZeroInstance to have zero value in OrthogonalSummation::Sum() initially.

  • the problem is that e.g. double's default cstor does not initialize to 0 as the underlying representation is not a bunch of 0 bits as given by the OS.
  • Property mode set to 100644
File size: 11.1 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2012 University of Bonn. All rights reserved.
5 *
6 *
7 * This file is part of MoleCuilder.
8 *
9 * MoleCuilder is free software: you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation, either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * MoleCuilder is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
21 */
22
23/*
24 * Histogram.cpp
25 *
26 * Created on: Jul 26, 2012
27 * Author: heber
28 */
29
30
31// include config.h
32#ifdef HAVE_CONFIG_H
33#include <config.h>
34#endif
35
36#include "CodePatterns/MemDebug.hpp"
37
38#include "Histogram.hpp"
39
40#include <algorithm>
41#include <cmath>
42#include <iostream>
43#include <sstream>
44#include <string>
45
46#include <boost/bind.hpp>
47#include <boost/foreach.hpp>
48
49#include "CodePatterns/Assert.hpp"
50#include "CodePatterns/Log.hpp"
51
52/** Tiny helper struct to create equally spaced bins with count of zero.
53 *
54 */
55struct BinCreator_t {
56 BinCreator_t( const double lowerend, const double _width ) :
57 currentstart(lowerend),
58 width(_width)
59 {}
60 Histogram::Bin_t operator()() {
61 Histogram::Bin_t bin( make_pair( currentstart, 0.) );
62 currentstart+=width;
63 return bin;
64 }
65private:
66 double currentstart;
67 const double width;
68};
69
70
71// see http://stackoverflow.com/questions/634087/stdcopy-to-stdcout-for-stdpair
72// printing a pair is not easy, especially so with ostream_iterator as it cannot find
73// overload operator<<() as it is not in namespace std.
74template <class T>
75struct PrintPair : public std::unary_function<T, void>
76{
77 std::ostream& os;
78 PrintPair(std::ostream& strm) : os(strm) {}
79
80 void operator()(const T& elem) const
81 {
82 os << "(" << elem.first << "," << elem.second << ") ";
83 }
84};
85
86std::ostream & operator<<(std::ostream &ost, const Histogram::Bin_t &elem)
87{
88 ost << "(" << elem.first << "," << elem.second << ") ";
89 return ost;
90}
91
92std::ostream & operator<<(std::ostream &ost, const Histogram &histogram)
93{
94 for (Histogram::Bins_t::const_iterator iter = histogram.bins.begin();
95 iter != histogram.bins.end(); ++iter)
96 ost << *iter;
97 return ost;
98}
99
100Histogram::Histogram(const samples_t &samples) :
101 binwidth(0.5),
102 offset(0.)
103{
104 if (!samples.empty()) {
105 // set offset to first value
106 const_cast<BinLowerEnd &>(offset) = *samples.begin();
107 // set binwidth to statistical sensible value
108 MinMax_t MinMax = getMinMaxFromSamples(samples);
109 const_cast<double &>(binwidth) = (*(MinMax.second) - *(MinMax.first))/pow(samples.size(), 1./3.);
110 // and add all samples to these histogram bins
111 addSamples(samples);
112 }
113}
114
115Histogram::Histogram(const samples_t &samples, const BinLowerEnd _offset, const double _binwidth) :
116 binwidth(_binwidth),
117 offset(_offset)
118{
119 addSamples(samples);
120}
121
122Histogram::MinMax_t Histogram::getMinMaxFromSamples(const samples_t &samples) const
123{
124 std::pair<samples_t::const_iterator, samples_t::const_iterator> returnpair;
125 returnpair.first = min_element(samples.begin(), samples.end());
126 returnpair.second = max_element(samples.begin(), samples.end());
127 ASSERT((returnpair.second != samples.end()) || (returnpair.first != samples.end()),
128 "Histogram::Histogram() - cannot find min/max despite non-empty range.");
129 return returnpair;
130}
131
132void Histogram::addSamples(const samples_t &samples)
133{
134 // build histogram from samples
135 if (!samples.empty()) {
136 // 1. get min and max and determine width
137 MinMax_t MinMax = getMinMaxFromSamples(samples);
138// LOG(1, "DEBUG: min is " << *(MinMax.first) << " and max is " << *(MinMax.second) << ".");
139
140 // 2. create each bin
141 {
142 std::vector<Bin_t> vectorbins;
143 const BinLowerEnd HistogramStart = getLowerEnd(*(MinMax.first));
144 BinCreator_t BinCreator( HistogramStart, binwidth );
145 // we need one extra bin for get...Bin()'s find to work properly
146 const int CountBins = ceil((*(MinMax.second) - HistogramStart)/binwidth)+1;
147 vectorbins.resize(CountBins+1, Bin_t( make_pair(0., 0.) ) );
148 std::generate( vectorbins.begin(), vectorbins.end(), BinCreator );
149 bins.insert(vectorbins.begin(), vectorbins.end());
150 }
151
152 // 3. place each sample into bin
153 BOOST_FOREACH( double value, samples) {
154 const Bins_t::iterator biniter = getLowerEndBin(value);
155 ASSERT( biniter != bins.end(),
156 "Histogram::Histogram() - cannot find bin for value from given samples.");
157 // (make bin count normalized, i.e. always equal to surface area of 1)
158 biniter->second += 1./binwidth;
159 }
160 LOG(2, "DEBUG: Bins are " << printBins() << ".");
161 } else {
162// LOG(1, "INFO: Given vector of samples is empty.");
163 }
164}
165
166std::string Histogram::printBins() const
167{
168 std::stringstream output;
169 std::for_each( bins.begin(), bins.end(), PrintPair<Bin_t>(output));
170 return output.str();
171}
172
173void Histogram::extendMissingBins(const BinLowerEnd LowerEnd, const BinLowerEnd NextLowerEnd)
174{
175 Bins_t::const_iterator loweriter = getLowerEndBin(LowerEnd);
176 Bins_t::const_iterator upperiter = getHigherEndBin(NextLowerEnd);
177 if (loweriter == bins.end()) {
178 // we need bins beneath our first, add them
179 for(BinLowerEnd offset = getLowerEnd(LowerEnd);
180 offset < getLowerEnd(NextLowerEnd);
181 offset += binwidth)
182 bins.insert( std::make_pair (offset, 0. ) );
183// LOG(2, "DEBUG: Bins after adding empties before start are " << printBins() << ".");
184 loweriter = getLowerEndBin(LowerEnd);
185 ASSERT( loweriter != bins.end(),
186 "Histogram::extendMissingBins() - still no lower bound after adding empties.");
187 }
188 if (upperiter == bins.end()) {
189 // we need bins after our last, add them
190 for(BinLowerEnd offset = getLowerEnd(LowerEnd);
191 offset <= getLowerEnd(NextLowerEnd)+binwidth;
192 offset += binwidth)
193 bins.insert( std::make_pair (offset, 0. ) );
194// LOG(2, "DEBUG: Bins after adding empties after end are " << printBins() << ".");
195 upperiter = getHigherEndBin(NextLowerEnd);
196 ASSERT( upperiter != bins.end(),
197 "Histogram::extendMissingBins() - still no upper bound after adding empties.");
198 }
199}
200
201void Histogram::superposeOtherHistogram(const Histogram &other, const double prefactor)
202{
203 // go through each of the other histogram's bins
204 Bins_t::const_iterator enditer = --other.bins.end(); // (except internal last one)
205 for (Bins_t::const_iterator biniter = other.bins.begin();
206 biniter != enditer; /* we advance ourselves in loop */) {
207 const Bin_t &bin = *biniter;
208 ++biniter;
209 const Bin_t &nextbin = *biniter;
210
211// LOG(2, "DEBUG: Current bin is " << bin << ", next bin is " << nextbin << ".");
212
213 // Check first whether start or end actually fit into our histogram, if not extend.
214 extendMissingBins(bin.first, nextbin.first);
215
216 // The bin will in general not fit into one bin in this histogram, but overlap.
217 // Hence, we determine the contribution of the bin to each bin in this histogram
218 // its overlaps into and add this weight to the bin.
219 Bins_t::const_iterator loweriter = getLowerEndBin(bin.first);
220 Bins_t::const_iterator upperiter = getHigherEndBin(nextbin.first);
221
222 ASSERT( loweriter->first < upperiter->first,
223 "Histogram::superposeOtherHistogram() - the bin range is invalid.");
224// LOG(2, "DEBUG: bin range here is ["
225// << loweriter->first << "," << upperiter->first << ").");
226
227 // Next, we create a vector of offsets
228 typedef std::vector< BinLowerEnd > offsets_t;
229 offsets_t offsets;
230 {
231 offsets.push_back(bin.first);
232 Bins_t::const_iterator iter = loweriter;
233 for (++iter; iter != upperiter; ++iter)
234 if (offsets.back() != iter->first)
235 offsets.push_back(iter->first);
236 if (offsets.back() != nextbin.first)
237 offsets.push_back(nextbin.first);
238// LOG(2, "DEBUG: Offsets are " << offsets << ".");
239 }
240
241 // then, we go through the offsets but the last one and add the respective area
242 {
243 offsets_t::const_iterator iter = offsets.begin();
244 offsets_t::const_iterator nextiter = ++offsets.begin();
245 for (; iter != --offsets.end(); ++iter, ++nextiter) {
246 const double length = *nextiter - *iter;
247 const double weight = bin.second * (length/binwidth);
248 Bins_t::iterator filliter = getLowerEndBin(*iter);
249 filliter->second += prefactor * weight;
250 }
251 }
252
253 {
254 std::stringstream output;
255 std::for_each( bins.begin(), bins.end(), PrintPair<Bin_t>(output));
256 LOG(2, "DEBUG: Bins are after summation " << output.str() << ".");
257 }
258 }
259}
260
261Histogram& Histogram::operator=(const Histogram &other)
262{
263 // check for self-assigment
264 if (&other != this) {
265 bins.clear();
266 const_cast<BinLowerEnd &>(offset) = other.offset;
267 const_cast<double &>(binwidth) = other.binwidth;
268 bins.insert(other.bins.begin(), other.bins.end());
269 }
270 return *this;
271}
272
273Histogram& Histogram::operator+=(const Histogram &other)
274{
275 superposeOtherHistogram(other, +1.);
276 LOG(2, "DEBUG: Bins after addition are " << printBins() << ".");
277 return *this;
278}
279
280Histogram& Histogram::operator-=(const Histogram &other)
281{
282 superposeOtherHistogram(other, -1.);
283 LOG(2, "DEBUG: Bins after subtraction are " << printBins() << ".");
284 return *this;
285}
286
287bool Histogram::isEmpty() const
288{
289 bool status = true;
290 for (Bins_t::const_iterator iter = bins.begin(); iter != bins.end(); ++iter)
291 status &= iter->second == 0;
292 return status;
293}
294
295Histogram::Bins_t::iterator Histogram::getLowerEndBin(const double _value)
296{
297 // lower bound returns key that is equal or greater
298 Bins_t::iterator iter = bins.lower_bound(_value);
299 if (iter != bins.end()) {
300 // if we are not on the boundary we always have to step back by one
301 if (_value != iter->first) {
302 if (iter != bins.begin()) {
303 --iter;
304 } else {
305 iter = bins.end();
306 }
307 } else if (iter == --bins.end()) {
308 // if we actually are on boundary of "last bin", set to end
309 iter = bins.end();
310 }
311 }
312 return iter;
313}
314
315Histogram::Bins_t::iterator Histogram::getHigherEndBin(const double _value)
316{
317 // upper bound return key that is strictly greater
318 Bins_t::iterator iter = bins.upper_bound(_value);
319 // if we are on the boundary we have to step back by one
320 if (iter != bins.end())
321 if (_value == iter->first)
322 if (iter != bins.begin())
323 --iter;
324
325 return iter;
326}
327
328Histogram::BinLowerEnd Histogram::getLowerEnd(const double _value) const
329{
330 BinLowerEnd start = _value - offset;
331 // then divide by binwidth
332 const int integral = floor(start/binwidth);
333 //const double fraction = fmod(start,binwidth);
334 start = offset + binwidth*integral;
335 return start;
336}
337
338double Histogram::area() const
339{
340 double area = 0.;
341 for(Bins_t::const_iterator iter = bins.begin(); iter != bins.end(); ++iter)
342 area += iter->second*binwidth;
343 return area;
344}
345
346template<> Histogram ZeroInstance<Histogram>()
347{
348 Histogram returnvalue;
349 return returnvalue;
350}
Note: See TracBrowser for help on using the repository browser.