LOOS  v2.3.2
TimeSeries.hpp
1 /*
2  This file is part of LOOS.
3 
4  LOOS (Lightweight Object-Oriented Structure library)
5  Copyright (c) 2008, Alan Grossfield
6  Department of Biochemistry and Biophysics
7  School of Medicine & Dentistry, University of Rochester
8 
9  This package (LOOS) 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 under version 3 of the License.
12 
13  This package is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 
23 
24 #if !defined(LOOS_TIMESERIES_HPP)
25 #define LOOS_TIMESERIES_HPP
26 
27 #include <vector>
28 #include <stdexcept>
29 #include <cmath>
30 #include <cstdlib>
31 #include <iostream>
32 #include <fstream>
33 #include <sstream>
34 
35 #include <loos_defs.hpp>
36 
37 namespace loos {
38 
40 
50 template<class T>
51 class TimeSeries {
52 public:
53  typedef typename std::vector<T>::iterator iterator;
54  typedef typename std::vector<T>::const_iterator const_iterator;
55  typedef const T& const_reference;
56  typedef T& reference;
57  typedef T value_type;
58 
59 
60  TimeSeries() {
61  init();
62  }
63 
64  TimeSeries(const std::vector<T> &inp) {
65  _data = inp;
66  }
67 
68  TimeSeries(const uint size, const T *array) {
69  _data.reserve(size);
70  for (unsigned int i=0; i<size; i++) {
71  _data.push_back(array[i]);
72  }
73  }
74 
75  TimeSeries(const TimeSeries<T> &inp) {
76  _data = inp._data;
77  }
78 
79  TimeSeries(const uint n) {
80  _data = std::vector<T>(n, 0);
81  }
82 
83  TimeSeries(const uint n, const T val) {
84  _data.assign(n, (T) val);
85  }
86 
88  void resize(const uint n, const T val= (T) 0.0) {
89  _data.resize(n, val);
90  }
91 
92 
96  TimeSeries (const std::string &filename, const int col=2) {
97  std::ifstream ifs(filename.c_str());
98  if (!ifs) {
99  throw(std::runtime_error("Cannot open timeseries file "
100  + filename));
101  }
102 
103  std::string line;
104  while (ifs.good()) {
105  getline(ifs, line);
106  if ( (line.substr(0,1) == "#") || (line.empty()) ){
107  // comment -- do nothing
108  }
109  else {
110  std::stringstream s(line);
111  double val;
112  int i=0;
113  while (i < col) {
114  if (!s.good()) {
115  throw(std::runtime_error("Problem reading timeseries file "
116  + filename));
117  }
118  s >> val;
119  i++;
120  }
121  _data.push_back(val);
122  }
123  }
124  }
125 
126 
127 
128  void init(void) {
129  _data.clear();
130  }
131 
132  void zero(void) {
133  _data.assign(_data.size(), (T)0.0);
134  }
135 
136 
137 #if !defined(SWIG)
138  T& operator[](const unsigned int i) {
139  return (_data.at(i)); // this way we get range checking
140  }
141 
142  const T& operator[](const unsigned int i) const {
143  return (_data.at(i)); // this way we get range checking
144  }
145 #endif // !defined(SWIG)
146 
147 
148  unsigned int size(void) const {
149  return (_data.size());
150  }
151 
152 
153  TimeSeries<T> operator+=(const T val) {
154  for (unsigned int i=0; i<_data.size(); i++) {
155  _data[i] += val;
156  }
157  return(*this);
158  }
159 
160  TimeSeries<T> operator+=(const TimeSeries<T> &rhs) {
161  if (_data.size() != rhs.size())
162  throw(std::runtime_error("mismatched timeseries sizes in +="));
163  for (unsigned int i=0; i<_data.size(); i++) {
164  _data[i] += rhs[i];
165  }
166  return(*this);
167  }
168 
169  TimeSeries<T> operator+(const T val) const {
170  TimeSeries<T> res(*this);
171  for (unsigned int i=0; i<_data.size(); i++) {
172  res[i] += val;
173  }
174  return(res);
175  }
176 
177 #if !defined(SWIG)
178  friend TimeSeries<T> operator+(const T lhs, const TimeSeries<T> &rhs) {
179  TimeSeries<T> res( rhs.size(), (T) 0.0 );
180  for (unsigned int i=0; i<rhs.size(); i++) {
181  res[i] = lhs + rhs[i];
182  }
183  return(res);
184  }
185 #endif
186 
187  TimeSeries<T> operator+(const TimeSeries<T> &rhs) const {
188  if (_data.size() != rhs.size())
189  throw(std::runtime_error("mismatched timeseries sizes in +="));
190 
191  TimeSeries<T> res(*this);
192  for (unsigned int i=0; i<res.size(); i++) {
193  res[i] += rhs[i];
194  }
195  return(res);
196  }
197 
198  TimeSeries<T> operator-=(const T val) {
199  for (unsigned int i=0; i<_data.size(); i++) {
200  _data[i] -= val;
201  }
202  return(*this);
203  }
204 
205  TimeSeries<T> operator-=(const TimeSeries<T> &rhs) {
206  if (_data.size() != rhs.size())
207  throw(std::runtime_error("mismatched sizes of time series"));
208  for (unsigned int i=0; i<_data.size(); i++) {
209  _data[i] -= rhs[i];
210  }
211  return(*this);
212  }
213 
214  TimeSeries<T> operator-(const T val) const {
215  TimeSeries<T> res(*this);
216  for (unsigned int i=0; i<_data.size(); i++) {
217  res[i] -= val;
218  }
219  return(res);
220  }
221 
222  TimeSeries<T> operator-(const TimeSeries<T> &rhs) const {
223  if (_data.size() != rhs.size())
224  throw(std::runtime_error("mismatched timeseries sizes in +="));
225 
226  TimeSeries<T> res(*this);
227  for (unsigned int i=0; i<res.size(); i++) {
228  res[i] -= rhs[i];
229  }
230  return(res);
231  }
232 
233 
234 
235  TimeSeries<T> operator-() const {
236  TimeSeries<T> res(*this);
237  for (unsigned int i=0; i<_data.size(); i++) {
238  res[i] = -res[i];
239  }
240  return(res);
241  }
242 
243 #if !defined(SWIG)
244  friend TimeSeries<T> operator-(const T lhs, const TimeSeries<T> &rhs) {
245  TimeSeries<T> res( rhs.size() );
246  for (unsigned int i=0; i<rhs.size(); i++) {
247  res[i] = lhs - rhs[i];
248  }
249  return(res);
250  }
251 #endif
252 
253  TimeSeries<T> operator*=(const T val) {
254  for (unsigned int i=0; i<_data.size(); i++) {
255  _data[i] *= val;
256  }
257  return(*this);
258  }
259 
260  TimeSeries<T> operator*(const T val) const {
261  TimeSeries<T> res(*this);
262  for (unsigned int i=0; i<_data.size(); i++) {
263  res[i] *= val;
264  }
265  return(res);
266  }
267 
268  TimeSeries<T> operator*=(const TimeSeries<T> &rhs) {
269  if ( _data.size() != rhs.size() )
270  throw(std::runtime_error("mismatched timeseries sizes in *="));
271  for (unsigned int i=0; i<_data.size(); i++) {
272  _data[i] *= rhs[i];
273  }
274  return(*this);
275  }
276 
277  TimeSeries<T> operator*(const TimeSeries<T> &rhs) const {
278  if ( _data.size() != rhs.size() )
279  throw(std::runtime_error("mismatched timeseries sizes in *="));
280 
281  TimeSeries<T> res(*this);
282  for (unsigned int i=0; i<_data.size(); i++) {
283  res[i] *= rhs[i];
284  }
285  return(res);
286  }
287 
288 #if !defined(SWIG)
289  friend TimeSeries<T> operator*(const T lhs, const TimeSeries<T> &rhs) {
290  TimeSeries<T> res( rhs.size(), (T) 0.0 );
291  for (unsigned int i=0; i<rhs.size(); i++) {
292  res[i] = lhs * rhs[i];
293  }
294  return(res);
295  }
296 #endif
297 
298  TimeSeries<T> operator/=(const T val) {
299  for (unsigned int i=0; i<_data.size(); i++) {
300  _data[i] /= val;
301  }
302  return(*this);
303  }
304 
305  TimeSeries<T> operator/(const T val) const {
306  TimeSeries<T> res(*this);
307  for (unsigned int i=0; i<_data.size(); i++) {
308  res[i] /= val;
309  }
310  return(res);
311  }
312 
313  TimeSeries<T> operator/=(const TimeSeries<T> &rhs) {
314  if ( _data.size() != rhs.size() )
315  throw(std::runtime_error("mismatched timeseries sizes in *="));
316  for (unsigned int i=0; i<_data.size(); i++) {
317  _data[i] /= rhs[i];
318  }
319  return(*this);
320  }
321 
322  TimeSeries<T> operator/(const TimeSeries<T> &rhs) const {
323  if ( _data.size() != rhs.size() )
324  throw(std::runtime_error("mismatched timeseries sizes in *="));
325 
326  TimeSeries<T> res(*this);
327  for (unsigned int i=0; i<_data.size(); i++) {
328  res[i] /= rhs[i];
329  }
330  return(res);
331  }
332 
333 #if !defined(SWIG)
334  friend TimeSeries<T> operator/(const T lhs, const TimeSeries<T> &rhs) {
335  TimeSeries<T> res( rhs.size() );
336  for (unsigned int i=0; i<rhs.size(); i++) {
337  res[i] = lhs / rhs[i];
338  }
339  return(res);
340  }
341 #endif
342 
343  TimeSeries<T> copy(void) const {
344  TimeSeries<T> result(_data.size(), 0.0);
345  for (unsigned int i = 0; i < _data.size(); i++) {
346  result[i] = _data[i];
347  }
348  return(result);
349  }
350 
353  void set_skip(unsigned int num_points) {
354  if (num_points < _data.size()) {
355  _data.erase(_data.begin(), _data.begin()+num_points);
356  }
357  else {
358  throw(std::out_of_range("num_points has invalid value in set_skip"));
359  }
360  }
361 
363  T average(void) const {
364  T ave = 0.0;
365  for (unsigned int i=0; i < _data.size(); i++) {
366  ave += _data[i];
367  }
368  ave /= _data.size();
369  return (ave);
370  }
371 
373  T variance(void) const {
374  T ave = 0.0;
375  T ave2 = 0.0;
376  for (unsigned int i=0; i < _data.size(); i++) {
377  ave += _data[i];
378  ave2 += _data[i] * _data[i];
379  }
380 
381  ave /= _data.size();
382  ave2 /= _data.size();
383  return (ave2 - ave*ave);
384  }
385 
387  T stdev(void) const {
388  return (sqrt(variance() ));
389  }
390 
395  T sterr(void) const {
396  return (stdev() / sqrt(_data.size()));
397  }
398 
402  TimeSeries<T> result(_data.size());
403  T sum = 0.0;
404  for (unsigned int i=0; i<_data.size(); i++) {
405  sum += _data[i];
406  result[i] = sum / (i+1);
407  }
408  return(result);
409 
410  }
411 
416  TimeSeries<T> windowed_average(const uint window) const {
417 
418  if (window > _data.size() )
419  throw(std::out_of_range("Error in windowed_average: window too large"));
420 
421  TimeSeries<T> result(_data.size() - window);
422  T sum = 0;
423  for (uint i=0; i<window; i++) {
424  sum += _data[i];
425  }
426  result[0] = sum / window;
427 
428 
429  for (unsigned int i=1; i < result.size(); i++) {
430  sum = sum - _data[i-1] + _data[i+window-1];
431  result[i] = sum / window;
432  }
433 
434  return(result);
435  }
436 
443  //
444  T block_var(const int num_blocks) const {
445  int points_per_block = size() / num_blocks;
446  T block_ave = 0.0;
447  T block_ave2 = 0.0;
448  for (int i=0; i<num_blocks; i++) {
449  T block_sum = 0.0;
450  int offset = i*points_per_block;
451  for (int j=0; j< points_per_block; j++) {
452  block_sum += _data[offset + j];
453  }
454  T ave = block_sum / points_per_block;
455  block_ave += ave;
456  block_ave2 += ave*ave;
457  }
458  block_ave /= num_blocks;
459  block_ave2 /= num_blocks;
460 
461  // The variance must be computed with N-1, not N, because
462  // we determine the mean from the data (as opposed to independently
463  // specifying it)
464  T ratio = num_blocks/(num_blocks-1.0);
465  return (block_ave2 - block_ave*block_ave)*ratio;
466  }
467 
468  TimeSeries<T> correl(const int max_time,
469  const int interval=1,
470  const bool normalize=true,
471  T tol=1.0e-8) const {
472 
473  TimeSeries<T> data = copy();
474  uint n = abs(max_time);
475  if (n > data.size()) {
476  throw(std::runtime_error("Can't take correlation time longer than time series"));
477  }
478 
479  n /= interval;
480  TimeSeries<T> c(n, 0.0);
481 
482  // normalize the data
483  if (normalize) {
484  data -= data.average();
485  T dev = data.stdev();
486 
487  // drop through if this is a constant array
488  if (dev < tol) {
489  c._data.assign(n, 1.0);
490  return(c);
491  }
492 
493  data /= dev;
494  }
495 
496  std::vector<int> num_pairs(n);
497  num_pairs.assign(n, 0);
498  // TODO: This is the O(N^2) way -- there are much faster
499  // algorithms for long time series
500  for (int i = 0; i < max_time; i+=interval) {
501  int index = i / interval;
502  for (unsigned int j = 0; j < data.size() - i; j++) {
503  c[index] += data[j] * data[j+i];
504  num_pairs[index]++;
505  }
506 
507  }
508 
509  // Divide each value by the number of pairs used to generated it
510  for (uint i = 0; i < n; i++) {
511  c[i] /= num_pairs[i];
512  }
513 
514  return(c);
515  }
516 
517  // Vector interface...
518  void push_back(const T& x) { _data.push_back(x); }
519 
520  iterator begin() { return(_data.begin()); }
521  iterator end() { return(_data.end()); }
522 
523 #if !defined(SWIG)
524  const_iterator begin() const { return(_data.begin()); }
525  const_iterator end() const { return(_data.end()); }
526 #endif
527 
528 
529 private:
530  std::vector<T> _data;
531 };
532 
533 
534 
535 typedef TimeSeries<double> dTimeSeries;
536 typedef TimeSeries<float> fTimeSeries;
537 
538 
539 
540 }
541 
542 #endif
T block_var(const int num_blocks) const
Definition: TimeSeries.hpp:444
void resize(const uint n, const T val=(T) 0.0)
Resize the TimeSeries by calling the underlying vector's resize.
Definition: TimeSeries.hpp:88
T stdev(void) const
Return standard deviation of time series.
Definition: TimeSeries.hpp:387
TimeSeries(const std::string &filename, const int col=2)
Definition: TimeSeries.hpp:96
void set_skip(unsigned int num_points)
Definition: TimeSeries.hpp:353
TimeSeries< T > running_average(void) const
Definition: TimeSeries.hpp:401
TimeSeries< T > windowed_average(const uint window) const
Definition: TimeSeries.hpp:416
Time Series Class.
Definition: TimeSeries.hpp:51
T average(void) const
Return average of time series.
Definition: TimeSeries.hpp:363
Namespace for most things not already encapsulated within a class.
T sterr(void) const
Definition: TimeSeries.hpp:395
T variance(void) const
Return variance of time series.
Definition: TimeSeries.hpp:373