TBTK
Smooth.h
Go to the documentation of this file.
1 /* Copyright 2017 Kristofer Bj√∂rnson
2  *
3  * Licensed under the Apache License, Version 2.0 (the "License");
4  * you may not use this file except in compliance with the License.
5  * You may obtain a copy of the License at
6  *
7  * http://www.apache.org/licenses/LICENSE-2.0
8  *
9  * Unless required by applicable law or agreed to in writing, software
10  * distributed under the License is distributed on an "AS IS" BASIS,
11  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12  * See the License for the specific language governing permissions and
13  * limitations under the License.
14  */
15 
23 #ifndef COM_DAFER45_TBTK_SMOOTH
24 #define COM_DAFER45_TBTK_SMOOTH
25 
26 #include "TBTK/Array.h"
27 #include "TBTK/Property/DOS.h"
28 #include "TBTK/TBTKMacros.h"
29 
30 #include <cmath>
31 #include <vector>
32 
33 namespace TBTK{
34 
35 class Smooth{
36 public:
38  static Array<double> gaussian(
39  const Array<double> &data,
40  double sigma,
41  int windowSize
42  );
43 
45  static std::vector<double> gaussian(
46  const std::vector<double> &data,
47  double sigma,
48  int windowSize
49  );
50 
52  static std::vector<double> gaussian(
53  const double *data,
54  unsigned int size,
55  double sigma,
56  int windowSize
57  );
58 
60  static Property::DOS gaussian(
61  const Property::DOS &dos,
62  double sigma,
63  int windowSize
64  );
65 private:
66 };
67 
69  const Array<double> &data,
70  double sigma,
71  int windowSize
72 ){
73  TBTKAssert(
74  windowSize > 0,
75  "Smooth::gaussian()",
76  "'windowSize' must be larger than zero.",
77  ""
78  );
79  TBTKAssert(
80  windowSize%2 == 1,
81  "Smooth::gaussian()",
82  "'windowSize' must be odd.",
83  ""
84  );
85  TBTKAssert(
86  data.getRanges().size() == 1,
87  "Smooth::gaussian()",
88  "Array must have rank 1, but the rank is "
89  << data.getRanges().size() << ".",
90  ""
91  );
92 
93  double normalization = 0;
94  for(int n = -windowSize/2; n <= windowSize/2; n++){
95  normalization += exp(-n*n/(2*sigma*sigma));
96  }
97  normalization = 1/normalization;
98 
99  Array<double> result({data.getRanges()[0]}, 0);
100  for(unsigned int n = 0; n < data.getRanges()[0]; n++){
101  for(
102  int c = std::max(0, (int)n - (int)windowSize/2);
103  c < std::min(
104  (int)n + (int)windowSize/2 + 1,
105  (int)data.getRanges()[0]
106  );
107  c++
108  ){
109  result[{(unsigned int)n}]
110  += data[
111  {(unsigned int)c}
112  ]*exp(-(c-n)*(c-n)/(2*sigma*sigma));
113  }
114  result[{(unsigned int)n}] *= normalization;
115  }
116 
117  return result;
118 }
119 
120 inline std::vector<double> Smooth::gaussian(
121  const std::vector<double> &data,
122  double sigma,
123  int windowSize
124 ){
125  TBTKAssert(
126  windowSize > 0,
127  "Smooth::gaussian()",
128  "'windowSize' must be larger than zero.",
129  ""
130  );
131  TBTKAssert(
132  windowSize%2 == 1,
133  "Smooth::gaussian()",
134  "'windowSize' must be odd.",
135  ""
136  );
137 
138  double normalization = 0;
139  for(int n = -windowSize/2; n <= windowSize/2; n++){
140  normalization += exp(-n*n/(2*sigma*sigma));
141  }
142  normalization = 1/normalization;
143 
144  std::vector<double> result;
145  for(unsigned int n = 0; n < data.size(); n++){
146  result.push_back(0);
147  for(
148  int c = std::max(0, (int)n - (int)windowSize/2);
149  c < std::min(
150  (int)n + (int)windowSize/2 + 1,
151  (int)data.size()
152  );
153  c++
154  ){
155  result.at(n) += data.at(c)*exp(-(c-n)*(c-n)/(2*sigma*sigma));
156  }
157  result.at(n) *= normalization;
158  }
159 
160  return result;
161 }
162 
163 inline std::vector<double> Smooth::gaussian(
164  const double *data,
165  unsigned int size,
166  double sigma,
167  int windowSize
168 ){
169  std::vector<double> dataVector;
170  for(unsigned int n = 0; n < size; n++)
171  dataVector.push_back(data[n]);
172 
173  return gaussian(
174  dataVector,
175  sigma,
176  windowSize
177  );
178 }
179 
181  const Property::DOS &dos,
182  double sigma,
183  int windowSize
184 ){
185  const double *data = dos.getData();
186  std::vector<double> dataVector;
187  for(unsigned int n = 0; n < dos.getSize(); n++)
188  dataVector.push_back(data[n]);
189 
190  double lowerBound = dos.getLowerBound();
191  double upperBound = dos.getUpperBound();
192  int resolution = dos.getResolution();
193  double scaledSigma = sigma/(upperBound - lowerBound)*resolution;
194 
195  std::vector<double> smoothedData = gaussian(
196  dataVector,
197  scaledSigma,
198  windowSize
199  );
200 
201  return Property::DOS(
202  lowerBound,
203  upperBound,
204  resolution,
205  smoothedData.data()
206  );
207 }
208 
209 }; //End of namespace TBTK
210 
211 #endif
const DataType * getData() const
Definition: AbstractProperty.h:226
Precompiler macros.
Multi-dimensional array.
Definition: Array.h:34
double getUpperBound() const
Definition: DOS.h:89
unsigned int getSize() const
Definition: AbstractProperty.h:217
int getResolution() const
Definition: DOS.h:93
Definition: ModelFactory.h:35
Property container for density of states (DOS).
static Array< double > gaussian(const Array< double > &data, double sigma, int windowSize)
Definition: Smooth.h:68
Property container for density of states (DOS).
Definition: DOS.h:32
const std::vector< unsigned int > & getRanges() const
Definition: Array.h:411
Definition: Smooth.h:35
double getLowerBound() const
Definition: DOS.h:85