LOOS  v2.3.2
water-autocorrel.cpp
1 /*
2  This file is part of LOOS.
3 
4  LOOS (Lightweight Object-Oriented Structure library)
5  Copyright (c) 2014, Tod D. Romo, 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 #include <loos.hpp>
24 
25 
26 using namespace std;
27 using namespace loos;
28 
29 
30 int main(int argc, char *argv[]) {
31  string hdr = invocationHeader(argc, argv);
32 
33  if (argc == 1) {
34  cerr << "Usage- " << argv[0] << " water_matrix [max-t] >output.asc\n";
35  exit(-1);
36  }
37 
38  int k = 1;
39  string matname(argv[k++]);
40  uint max_t = 0;
41  if (k != argc)
42  max_t = strtoul(argv[k++], 0, 10);
43 
45  cerr << "Reading matrix...\n";
46  readAsciiMatrix(matname, M);
47  uint m = M.rows();
48  uint n = M.cols();
49 
50  if (max_t == 0)
51  max_t = n/10;
52 
53  cerr << boost::format("Water matrix is %d x %d\n") % m % n;
54  cerr << "Processing- ";
55  vector< TimeSeries<double> > waters;
56  for (uint j=0; j<m; ++j) {
57  if (j % 250 == 0)
58  cerr << '.';
59 
60  vector<double> tmp(n, 0);
61  bool flag = false;
62  for (uint i=0; i<n; ++i) {
63  tmp[i] = M(j, i);
64  if (tmp[i])
65  flag = true;
66  }
67  if (flag) {
68  TimeSeries<double> wtmp(tmp);
69  waters.push_back(wtmp.correl(max_t));
70  }
71  }
72  cerr << boost::format(" done\nFound %d unique waters inside\n") % waters.size();
73  cout << "# " << hdr << endl;
74  for (uint j=0; j<max_t; ++j) {
75  vector<double> tmp(m, 0);
76  for (uint i=0; i<m; ++i)
77  tmp[i] = waters[i][j];
78  dTimeSeries dtmp(tmp);
79  cout << j << '\t' << dtmp.average() << '\t' << dtmp.stdev() << '\t' << dtmp.sterr() << endl;
80  }
81 }
STL namespace.
std::string invocationHeader(int argc, char *argv[])
Create an invocation header.
Definition: utils.cpp:124
Math::Matrix< T, P, S > readAsciiMatrix(std::istream &is)
Read in a matrix from a stream returning a newly created matrix.
Definition: MatrixRead.hpp:72
Time Series Class.
Definition: TimeSeries.hpp:51
Namespace for most things not already encapsulated within a class.