1 | /*
|
---|
2 | * analysis.hpp
|
---|
3 | *
|
---|
4 | * Created on: Oct 13, 2009
|
---|
5 | * Author: heber
|
---|
6 | */
|
---|
7 |
|
---|
8 | #ifndef ANALYSIS_HPP_
|
---|
9 | #define ANALYSIS_HPP_
|
---|
10 |
|
---|
11 | using namespace std;
|
---|
12 |
|
---|
13 | /*********************************************** includes ***********************************/
|
---|
14 |
|
---|
15 | // include config.h
|
---|
16 | #ifdef HAVE_CONFIG_H
|
---|
17 | #include <config.h>
|
---|
18 | #endif
|
---|
19 |
|
---|
20 | #include <cmath>
|
---|
21 | #include <fstream>
|
---|
22 |
|
---|
23 | // STL headers
|
---|
24 | #include <map>
|
---|
25 | #include <vector>
|
---|
26 |
|
---|
27 | #include "defs.hpp"
|
---|
28 |
|
---|
29 | #include "atom.hpp"
|
---|
30 | #include "Helpers/Verbose.hpp"
|
---|
31 |
|
---|
32 | /****************************************** forward declarations *****************************/
|
---|
33 |
|
---|
34 | class BoundaryTriangleSet;
|
---|
35 | class element;
|
---|
36 | class LinkedCell;
|
---|
37 | class Tesselation;
|
---|
38 | class Vector;
|
---|
39 |
|
---|
40 | /********************************************** definitions *********************************/
|
---|
41 |
|
---|
42 | typedef multimap<double, pair<atom *, atom *> > PairCorrelationMap;
|
---|
43 | typedef multimap<double, pair<atom *, const Vector *> > CorrelationToPointMap;
|
---|
44 | typedef multimap<double, pair<atom *, BoundaryTriangleSet *> > CorrelationToSurfaceMap;
|
---|
45 | typedef map<double, int> BinPairMap;
|
---|
46 |
|
---|
47 | /********************************************** declarations *******************************/
|
---|
48 |
|
---|
49 | PairCorrelationMap *PairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements);
|
---|
50 | CorrelationToPointMap *CorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point );
|
---|
51 | CorrelationToSurfaceMap *CorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Tesselation * const Surface, const LinkedCell *LC );
|
---|
52 | PairCorrelationMap *PeriodicPairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const int ranges[NDIM] );
|
---|
53 | CorrelationToPointMap *PeriodicCorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point, const int ranges[NDIM] );
|
---|
54 | CorrelationToSurfaceMap *PeriodicCorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] );
|
---|
55 | int GetBin ( const double value, const double BinWidth, const double BinStart );
|
---|
56 | void OutputCorrelation( ofstream * const file, const BinPairMap * const map );
|
---|
57 | void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map );
|
---|
58 | void OutputCorrelationToPoint( ofstream * const file, const CorrelationToPointMap * const map );
|
---|
59 | void OutputCorrelationToSurface( ofstream * const file, const CorrelationToSurfaceMap * const map );
|
---|
60 |
|
---|
61 |
|
---|
62 | /** Searches for lowest and highest value in a given map of doubles.
|
---|
63 | * \param *map map of doubles to scan
|
---|
64 | * \param &min minimum on return
|
---|
65 | * \param &max maximum on return
|
---|
66 | */
|
---|
67 | template <typename T> void GetMinMax ( T *map, double &min, double &max)
|
---|
68 | {
|
---|
69 | min = 0.;
|
---|
70 | max = 0.;
|
---|
71 | bool FirstMinFound = false;
|
---|
72 | bool FirstMaxFound = false;
|
---|
73 |
|
---|
74 | if (map == NULL) {
|
---|
75 | DoeLog(0) && (eLog()<< Verbose(0) << "Nothing to min/max, map is NULL!" << endl);
|
---|
76 | performCriticalExit();
|
---|
77 | return;
|
---|
78 | }
|
---|
79 |
|
---|
80 | for (typename T::iterator runner = map->begin(); runner != map->end(); ++runner) {
|
---|
81 | if ((min > runner->first) || (!FirstMinFound)) {
|
---|
82 | min = runner->first;
|
---|
83 | FirstMinFound = true;
|
---|
84 | }
|
---|
85 | if ((max < runner->first) || (!FirstMaxFound)) {
|
---|
86 | max = runner->first;
|
---|
87 | FirstMaxFound = true;
|
---|
88 | }
|
---|
89 | }
|
---|
90 | };
|
---|
91 |
|
---|
92 | /** Puts given correlation data into bins of given size (histogramming).
|
---|
93 | * Note that BinStart and BinEnd are desired to fill the complete range, even where Bins are zero. If this is
|
---|
94 | * not desired, give equal BinStart and BinEnd and the map will contain only Bins where the count is one or greater. If a
|
---|
95 | * certain start value is desired, give BinStart and a BinEnd that is smaller than BinStart, then only BinEnd will be
|
---|
96 | * calculated automatically, i.e. BinStart = 0. and BinEnd = -1. .
|
---|
97 | * Also note that the range is given inclusive, i.e. [ BinStart, BinEnd ].
|
---|
98 | * \param *map map of doubles to count
|
---|
99 | * \param BinWidth desired width of the binds
|
---|
100 | * \param BinStart first bin
|
---|
101 | * \param BinEnd last bin
|
---|
102 | * \return Map of doubles (the bins) with counts as values
|
---|
103 | */
|
---|
104 | template <typename T> BinPairMap *BinData ( T *map, const double BinWidth, const double BinStart, const double BinEnd )
|
---|
105 | {
|
---|
106 | BinPairMap *outmap = new BinPairMap;
|
---|
107 | int bin = 0;
|
---|
108 | double start = 0.;
|
---|
109 | double end = 0.;
|
---|
110 | pair <BinPairMap::iterator, bool > BinPairMapInserter;
|
---|
111 |
|
---|
112 | if (map == NULL) {
|
---|
113 | DoeLog(0) && (eLog()<< Verbose(0) << "Nothing to bin, is NULL!" << endl);
|
---|
114 | performCriticalExit();
|
---|
115 | return outmap;
|
---|
116 | }
|
---|
117 |
|
---|
118 | if (BinStart == BinEnd) { // if same, find range ourselves
|
---|
119 | GetMinMax( map, start, end);
|
---|
120 | } else if (BinEnd < BinStart) { // if BinEnd smaller, just look for new max
|
---|
121 | GetMinMax( map, start, end);
|
---|
122 | start = BinStart;
|
---|
123 | } else { // else take both values
|
---|
124 | start = BinStart;
|
---|
125 | end = BinEnd;
|
---|
126 | }
|
---|
127 | for (int runner = 0; runner <= ceil((end-start)/BinWidth); runner++)
|
---|
128 | outmap->insert( pair<double, int> ((double)runner*BinWidth+start, 0) );
|
---|
129 |
|
---|
130 | for (typename T::iterator runner = map->begin(); runner != map->end(); ++runner) {
|
---|
131 | bin = GetBin (runner->first, BinWidth, start);
|
---|
132 | BinPairMapInserter = outmap->insert ( pair<double, int> ((double)bin*BinWidth+start, 1) );
|
---|
133 | if (!BinPairMapInserter.second) { // bin already present, increase
|
---|
134 | BinPairMapInserter.first->second += 1;
|
---|
135 | }
|
---|
136 | }
|
---|
137 |
|
---|
138 | return outmap;
|
---|
139 | };
|
---|
140 |
|
---|
141 | #endif /* ANALYSIS_HPP_ */
|
---|