LOOS  v2.3.2
water-survival.cpp
1 // (c) 2014 Tod D. Romo, Grossfield Lab, URMC
2 
3 #include <loos.hpp>
4 
5 using namespace loos;
6 using namespace std;
7 
9 
10 
11 
12 int main(int argc, char *argv[]) {
13  string hdr = invocationHeader(argc, argv);
14 
15  if (argc == 1) {
16  cerr << "Usage- " << argv[0] << " water_matrix [max-t] >output.asc\n";
17  exit(-1);
18  }
19 
20  int k = 1;
21  string matname(argv[k++]);
22  uint max_t = 0;
23  if (k != argc)
24  max_t = strtoul(argv[k++], 0, 10);
25 
27  cerr << "Reading matrix...\n";
28  readAsciiMatrix(matname, M);
29  uint m = M.rows();
30  uint n = M.cols();
31 
32  if (max_t == 0)
33  max_t = n/10;
34 
35  cerr << boost::format("Water matrix is %d x %d\n") % m % n;
36 
37  cout << "# " << hdr << endl;
38  cout << "# tau\tavg\tstdev\tsterr\n";
39 
40  cerr << "Processing- ";
41 
42 
43  for (uint tau=0; tau<max_t; ++tau) {
44  vector<double> survivals;
45  if (tau % 100 == 0)
46  cerr << '.';
47 
48  for (uint j=0; j<m; ++j) {
49  uint inside = 0;
50  uint pairs = 0;
51 
52  for (uint t=0; t<n-tau-1; ++t)
53  if (M(j, t))
54  {
55  ++pairs;
56  if (M(j, t+tau))
57  ++inside;
58  }
59  if (pairs)
60  survivals.push_back(static_cast<double>(inside) / (pairs));
61  }
62  dTimeSeries ts(survivals);
63  cout << tau << '\t' << ts.average() << '\t' << ts.stdev() << '\t' << ts.sterr() << endl;
64  }
65 
66  cerr << " Done\n";
67 
68 }
69 
70 
71 
72 
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.