LOOS  v2.3.2
gridstat.cpp
1 /*
2  gridstat.cpp
3 
4  Simple grid statistics...
5 */
6 
7 /*
8  This file is part of LOOS.
9 
10  LOOS (Lightweight Object-Oriented Structure library)
11  Copyright (c) 2008, Tod D. Romo, Alan Grossfield
12  Department of Biochemistry and Biophysics
13  School of Medicine & Dentistry, University of Rochester
14 
15  This package (LOOS) is free software: you can redistribute it and/or modify
16  it under the terms of the GNU General Public License as published by
17  the Free Software Foundation under version 3 of the License.
18 
19  This package is distributed in the hope that it will be useful,
20  but WITHOUT ANY WARRANTY; without even the implied warranty of
21  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  GNU General Public License for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with this program. If not, see <http://www.gnu.org/licenses/>.
26 */
27 
28 
29 #include <loos.hpp>
30 #include <boost/format.hpp>
31 #include <DensityGrid.hpp>
32 
33 using namespace std;
34 using namespace loos;
35 using namespace loos::DensityTools;
36 
37 
38 
39 
40 double avgDens(DensityGrid<double>& grid) {
41  long i, n = grid.maxGridIndex();
42  double avg = 0.0;
43 
44  for (i=0; i<n; i++)
45  avg += grid(i);
46 
47  return(avg/n);
48 }
49 
50 
51 double zavgDens(DensityGrid<double>& grid) {
52  long n = grid.maxGridIndex();
53  long m = 0;
54  double avg = 0.0;
55 
56  for (long i=0; i<n; ++i) {
57  double d = grid(i);
58  if (d > 0.0) {
59  avg += d;
60  ++m;
61  }
62  }
63 
64  avg /= m;
65  return(avg);
66 }
67 
68 
69 double stdDens(DensityGrid<double>& grid, const double avg) {
70  long i, n = grid.maxGridIndex();
71  double std = 0.0;
72 
73  for (i=0; i<n; i++)
74  std += (grid(i) - avg) * (grid(i) - avg);
75 
76  return(sqrt(std/(n-1.0)));
77 }
78 
79 
80 double zstdDens(DensityGrid<double>& grid, const double avg) {
81  long i, n = grid.maxGridIndex();
82  double std = 0.0;
83 
84  long m = 0;
85  for (i=0; i<n; i++)
86  if (grid(i) > 0.0) {
87  std += (grid(i) - avg) * (grid(i) - avg);
88  ++m;
89  }
90 
91  return(sqrt(std/(m-1.0)));
92 }
93 
94 
95 double maxDens(DensityGrid<double>& grid) {
96  long i, n = grid.maxGridIndex();
97  double max = 0.0;
98 
99  for (i=0; i<n; i++)
100  if (grid(i) > max)
101  max = grid(i);
102 
103  return(max);
104 }
105 
106 
107 
108 void quickHist(DensityGrid<double>& grid, const double x, const int nbins) {
109  long *bins = new long[nbins];
110  double delta = x / nbins;
111 
112  for (int i=0; i<nbins; i++)
113  bins[i] = 0;
114 
115  long n = grid.maxGridIndex();
116  for (long i=0; i<n; i++) {
117  int k = static_cast<int>(grid(i) / delta);
118  assert(k <= nbins && k >= 0);
119  if (k == nbins)
120  k = nbins - 1;
121  ++bins[k];
122  }
123 
124  cout << "Quick histogram\n";
125  cout << "---------------\n";
126  for (int i=0; i<nbins; i++) {
127  cout << setprecision(6) << setw(10) << i*delta << "\t" << setprecision(4) << static_cast<double>(i)/nbins << "\t";
128  cout << setw(10) << bins[i] << "\t" << setprecision(4) << static_cast<double>(bins[i]) / n << endl;
129  }
130 
131  delete[] bins;
132 }
133 
134 
135 void zAverage(DensityGrid<double>& grid, const int nbins) {
136  DensityGridpoint dims = grid.gridDims();
137 
138  int chunk_size = dims[2] / nbins;
139  long volume = chunk_size * dims[1] * dims[0];
140 
141  cout << endl;
142  cout << "Z-slice averages\n";
143  cout << "----------------\n";
144 
145  int kk = 0;
146  for (int k = 0; k<nbins; k++) {
147  // Calculate z-range...
148  DensityGridpoint bottom(0,0,k*chunk_size);
149  DensityGridpoint top(0,0,chunk_size*(k+1));
150 
151  GCoord wbottom = grid.gridToWorld(bottom);
152  GCoord wtop = grid.gridToWorld(top);
153 
154  double avg = 0.0;
155 
156  for (int sk = 0; sk < chunk_size && sk+kk < dims[2]; sk++, kk++)
157  for (int j=0; j<dims[1]; j++)
158  for (int i=0; i<dims[0]; i++)
159  avg += grid(kk, j, i);
160 
161  avg /= volume;
162  cout << kk << "\t" << wbottom.z() << "\t" << wtop.z() << "\t" << avg << endl;
163  }
164 
165  if (kk < dims[2]) {
166  DensityGridpoint bottom(0,0,kk);
167  GCoord wbottom = grid.gridToWorld(bottom);
168  double avg = 0.0;
169 
170  volume = 0;
171  for (; kk < dims[2]; kk++)
172  for (int j=0; j<dims[1]; j++)
173  for (int i=0; i<dims[0]; i++, volume++)
174  avg += grid(kk, j, i);
175 
176  DensityGridpoint top(0,0,kk);
177  GCoord wtop = grid.gridToWorld(top);
178 
179  avg /= volume;
180  cout << kk << "\t" << wbottom.z() << "\t" << wtop.z() << "\t" << avg << endl;
181  cout << "Warning- last row adjusted\n";
182  }
183 }
184 
185 
186 
187 double rmsdDens(DensityGrid<double>& grid, const double avg) {
188  double rms = 0.0;
189 
190  for (long i=0; i<grid.maxGridIndex(); ++i) {
191  double d = grid(i) - avg;
192  rms += d*d;
193  }
194 
195  rms /= grid.maxGridIndex();
196  return(sqrt(rms));
197 }
198 
199 
200 int main(int argc, char *argv[]) {
201  if (argc != 3) {
202  cerr <<
203  "Usage- gridstat bins zbins <file.grid\n"
204  "\n"
205  "Displays some basic statistics about the density in a grid.\n"
206  "Bins is the number of bins for histogramming the density values.\n"
207  "Zbins is the number of bins in Z (really, K) to calculate density\n"
208  "statistics (useful for membrane systems).\n"
209  "Requires a double-precision floating point grid\n";
210  exit(-1);
211  }
212 
213  double nbins = strtod(argv[1], 0);
214  double zbins = strtod(argv[2], 0);
215 
216  DensityGrid<double> grid;
217  cin >> grid;
218 
219  cout << "Read in grid of size " << grid.gridDims() << endl;
220  cout << "Range is " << grid.minCoord() << " to " << grid.maxCoord() << endl;
221 
222  double gavg = avgDens(grid);
223  double grmsd = rmsdDens(grid, gavg);
224  double gzavg = zavgDens(grid);
225  double gstd = stdDens(grid, gavg);
226  double gzstd = zstdDens(grid, gzavg);
227  double gmax = maxDens(grid);
228 
229  cout << "\n\n* Grid Density Statistics *\n";
230  cout << "Grid density is " << gavg << " (" << gstd << ")\n";
231  cout << "Grid rmsd is " << grmsd << endl;
232  cout << "Grid non-zero avg is " << gzavg << " (" << gzstd << ")\n";
233  cout << "Max density is " << gmax << endl << endl;
234  quickHist(grid, gmax, nbins);
235  zAverage(grid, zbins);
236 
237 
238 }
STL namespace.
Namespace for Density package.
Definition: DensityGrid.hpp:48
loos::GCoord gridToWorld(const DensityGridpoint &v) const
Converts grid coords to real-space (world) coords.
Namespace for most things not already encapsulated within a class.