TBTK
Need a break? Support the development by playing Polarity Puzzles
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/Property/LDOS.h"
30 #include "TBTK/TBTKMacros.h"
31 
32 #include <cmath>
33 #include <vector>
34 
35 namespace TBTK{
36 
37 class Smooth{
38 public:
40  template<typename DataType>
42  const Array<DataType> &data,
43  double sigma,
44  int windowSize
45  );
46 
48  template<typename DataType>
49  static std::vector<DataType> gaussian(
50  const std::vector<DataType> &data,
51  double sigma,
52  int windowSize
53  );
54 
56  template<typename DataType>
57  static std::vector<DataType> gaussian(
58  const DataType *data,
59  unsigned int size,
60  double sigma,
61  int windowSize
62  );
63 
65  static Property::DOS gaussian(
66  const Property::DOS &dos,
67  double sigma,
68  int windowSize
69  );
70 
72  static Property::LDOS gaussian(
73  const Property::LDOS &ldos,
74  double sigma,
75  int windowSize
76  );
77 
80  const Property::SpinPolarizedLDOS &ldos,
81  double sigma,
82  int windowSize
83  );
84 private:
85 };
86 
87 template<typename DataType>
89  const Array<DataType> &data,
90  double sigma,
91  int windowSize
92 ){
93  TBTKAssert(
94  windowSize > 0,
95  "Smooth::gaussian()",
96  "'windowSize' must be larger than zero.",
97  ""
98  );
99  TBTKAssert(
100  windowSize%2 == 1,
101  "Smooth::gaussian()",
102  "'windowSize' must be odd.",
103  ""
104  );
105  TBTKAssert(
106  data.getRanges().size() == 1,
107  "Smooth::gaussian()",
108  "Array must have rank 1, but the rank is "
109  << data.getRanges().size() << ".",
110  ""
111  );
112 
113  DataType normalization = 0;
114  for(int n = -windowSize/2; n <= windowSize/2; n++){
115  normalization += exp(-n*n/(2*sigma*sigma));
116  }
117  normalization = DataType(1)/normalization;
118 
119  Array<DataType> result({data.getRanges()[0]}, 0);
120  for(int n = 0; n < (int)data.getRanges()[0]; n++){
121  for(
122  int c = std::max(0, (int)n - (int)windowSize/2);
123  c < std::min(
124  (int)n + (int)windowSize/2 + 1,
125  (int)data.getRanges()[0]
126  );
127  c++
128  ){
129  result[{(unsigned int)n}]
130  += data[
131  {(unsigned int)c}
132  ]*exp(-(c-n)*(c-n)/(2*sigma*sigma));
133  }
134  result[{(unsigned int)n}] *= normalization;
135  }
136 
137  return result;
138 }
139 
140 template<typename DataType>
141 inline std::vector<DataType> Smooth::gaussian(
142  const std::vector<DataType> &data,
143  double sigma,
144  int windowSize
145 ){
146  TBTKAssert(
147  windowSize > 0,
148  "Smooth::gaussian()",
149  "'windowSize' must be larger than zero.",
150  ""
151  );
152  TBTKAssert(
153  windowSize%2 == 1,
154  "Smooth::gaussian()",
155  "'windowSize' must be odd.",
156  ""
157  );
158 
159  DataType normalization = 0;
160  for(int n = -windowSize/2; n <= windowSize/2; n++){
161  normalization += exp(-n*n/(2*sigma*sigma));
162  }
163  normalization = DataType(1)/normalization;
164 
165  std::vector<DataType> result;
166  for(int n = 0; n < (int)data.size(); n++){
167  result.push_back(0);
168  for(
169  int c = std::max(0, (int)n - (int)windowSize/2);
170  c < std::min(
171  (int)n + (int)windowSize/2 + 1,
172  (int)data.size()
173  );
174  c++
175  ){
176  result.at(n) += data.at(c)*exp(-(c-n)*(c-n)/(2*sigma*sigma));
177  }
178  result.at(n) *= normalization;
179  }
180 
181  return result;
182 }
183 
184 template<typename DataType>
185 inline std::vector<DataType> Smooth::gaussian(
186  const DataType *data,
187  unsigned int size,
188  double sigma,
189  int windowSize
190 ){
191  std::vector<DataType> dataVector;
192  for(unsigned int n = 0; n < size; n++)
193  dataVector.push_back(data[n]);
194 
195  return gaussian(
196  dataVector,
197  sigma,
198  windowSize
199  );
200 }
201 
203  const Property::DOS &dos,
204  double sigma,
205  int windowSize
206 ){
207  const std::vector<double> &data = dos.getData();
208  std::vector<double> dataVector;
209  for(unsigned int n = 0; n < data.size(); n++)
210  dataVector.push_back(data[n]);
211 
212  double lowerBound = dos.getLowerBound();
213  double upperBound = dos.getUpperBound();
214  int resolution = dos.getResolution();
215  double scaledSigma = sigma/(upperBound - lowerBound)*resolution;
216 
217  std::vector<double> smoothedData = gaussian(
218  dataVector,
219  scaledSigma,
220  windowSize
221  );
222 
223  return Property::DOS(
224  lowerBound,
225  upperBound,
226  resolution,
227  smoothedData.data()
228  );
229 }
230 
232  const Property::LDOS &ldos,
233  double sigma,
234  int windowSize
235 ){
236  Property::LDOS newLdos = ldos;
237  std::vector<double> &newData = newLdos.getDataRW();
238 
239  const std::vector<double> &data = ldos.getData();
240  unsigned int blockSize = ldos.getBlockSize();
241  unsigned int numBlocks = ldos.getSize()/blockSize;
242  double lowerBound = ldos.getLowerBound();
243  double upperBound = ldos.getUpperBound();
244  int resolution = ldos.getResolution();
245  double scaledSigma = sigma/(upperBound - lowerBound)*resolution;
246  for(unsigned int block = 0; block < numBlocks; block++){
247  std::vector<double> blockData(blockSize);
248  for(unsigned int n = 0; n < blockSize; n++)
249  blockData[n] = data[block*blockSize + n];
250 
251 
252  std::vector<double> smoothedData = gaussian(
253  blockData,
254  scaledSigma,
255  windowSize
256  );
257 
258  for(unsigned int n = 0; n < blockSize; n++)
259  newData[block*blockSize + n] = smoothedData[n];
260  }
261 
262  return newLdos;
263 }
264 
266  const Property::SpinPolarizedLDOS &spinPolarizedLDOS,
267  double sigma,
268  int windowSize
269 ){
270  Property::SpinPolarizedLDOS newSpinPolarizedLDOS = spinPolarizedLDOS;
271  std::vector<SpinMatrix> &newData = newSpinPolarizedLDOS.getDataRW();
272 
273  const std::vector<SpinMatrix> &data = spinPolarizedLDOS.getData();
274  unsigned int blockSize = spinPolarizedLDOS.getBlockSize();
275  unsigned int numBlocks = spinPolarizedLDOS.getSize()/blockSize;
276  double lowerBound = spinPolarizedLDOS.getLowerBound();
277  double upperBound = spinPolarizedLDOS.getUpperBound();
278  int resolution = spinPolarizedLDOS.getResolution();
279  double scaledSigma = sigma/(upperBound - lowerBound)*resolution;
280  for(unsigned int block = 0; block < numBlocks; block++){
281  std::vector<SpinMatrix> blockData(blockSize);
282  for(unsigned int n = 0; n < blockSize; n++)
283  blockData[n] = data[block*blockSize + n];
284 
285  std::vector<
286  std::vector<std::complex<double>>
287  > spinMatrixComponents(
288  4,
289  std::vector<std::complex<double>>(blockSize)
290  );
291  for(unsigned int n = 0; n < blockSize; n++){
292  spinMatrixComponents[0][n] = blockData[n].at(0, 0);
293  spinMatrixComponents[1][n] = blockData[n].at(0, 1);
294  spinMatrixComponents[2][n] = blockData[n].at(1, 0);
295  spinMatrixComponents[3][n] = blockData[n].at(1, 1);
296  }
297 
298  for(unsigned int n = 0; n < 4; n++){
299  spinMatrixComponents[n] = gaussian(
300  spinMatrixComponents[n],
301  scaledSigma,
302  windowSize
303  );
304  }
305 
306  for(unsigned int n = 0; n < blockSize; n++){
307  newData[block*blockSize + n].at(0, 0)
308  = spinMatrixComponents[0][n];
309  newData[block*blockSize + n].at(0, 1)
310  = spinMatrixComponents[1][n];
311  newData[block*blockSize + n].at(1, 0)
312  = spinMatrixComponents[2][n];
313  newData[block*blockSize + n].at(1, 1)
314  = spinMatrixComponents[3][n];
315  }
316  }
317 
318  return newSpinPolarizedLDOS;
319 }
320 
321 }; //End of namespace TBTK
322 
323 #endif
Property container for the local density of states (LDOS).
Definition: LDOS.h:48
unsigned int getResolution() const
Definition: EnergyResolvedProperty.h:881
static Array< DataType > gaussian(const Array< DataType > &data, double sigma, int windowSize)
Definition: Smooth.h:88
Precompiler macros.
Property container for local density of states (LDOS).
const std::vector< DataType > & getData() const
Definition: AbstractProperty.h:567
double getLowerBound() const
Definition: EnergyResolvedProperty.h:857
const std::vector< unsigned int > & getRanges() const
Definition: Array.h:1228
Multi-dimensional array.
Property container for spin-polarized local density of states (spin-polarized LDOS).
double getUpperBound() const
Definition: EnergyResolvedProperty.h:869
Multi-dimensional array.
Definition: Array.h:98
unsigned int getBlockSize() const
Definition: AbstractProperty.h:555
std::vector< DataType > & getDataRW()
Definition: AbstractProperty.h:572
Definition: Boolean.h:32
Property container for spin-polarized local density of states (spin-polarized LDOS).
Definition: SpinPolarizedLDOS.h:46
Property container for density of states (DOS).
unsigned int getSize() const
Definition: AbstractProperty.h:560
Property container for density of states (DOS).
Definition: DOS.h:48
Definition: Smooth.h:37