LOOS  v2.3.2
gridautoscale.cpp
1 /*
2  gridautoscale.cpp
3 
4 */
5 
6 /*
7  This file is part of LOOS.
8 
9  LOOS (Lightweight Object-Oriented Structure library)
10  Copyright (c) 2012, Tod D. Romo, Alan Grossfield
11  Department of Biochemistry and Biophysics
12  School of Medicine & Dentistry, University of Rochester
13 
14  This package (LOOS) is free software: you can redistribute it and/or modify
15  it under the terms of the GNU General Public License as published by
16  the Free Software Foundation under version 3 of the License.
17 
18  This package is distributed in the hope that it will be useful,
19  but WITHOUT ANY WARRANTY; without even the implied warranty of
20  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  GNU General Public License for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with this program. If not, see <http://www.gnu.org/licenses/>.
25 */
26 
27 
28 #include <loos.hpp>
29 #include <DensityGrid.hpp>
30 
31 using namespace std;
32 using namespace loos;
33 using namespace loos::DensityTools;
34 
36 
37 
38 double findPeakDensitySlice(const Grid& grid, const int nbins) {
39  DensityGridpoint dims = grid.gridDims();
40 
41  int chunk_size = dims[2] / nbins;
42  long volume = chunk_size * dims[1] * dims[0];
43 
44  double max_avg_density = 0.0;
45 
46  int kk = 0;
47  for (int k = 0; k<nbins; k++) {
48  // Calculate z-range...
49  DensityGridpoint bottom(0,0,k*chunk_size);
50  DensityGridpoint top(0,0,chunk_size*(k+1));
51 
52 
53  double avg = 0.0;
54 
55  for (int sk = 0; sk < chunk_size && sk+kk < dims[2]; sk++, kk++)
56  for (int j=0; j<dims[1]; j++)
57  for (int i=0; i<dims[0]; i++)
58  avg += grid(kk, j, i);
59 
60  avg /= volume;
61  if (avg >= max_avg_density)
62  max_avg_density = avg;
63  }
64 
65  if (kk < dims[2]) {
66  double avg = 0.0;
67 
68  volume = 0;
69  for (; kk < dims[2]; kk++)
70  for (int j=0; j<dims[1]; j++)
71  for (int i=0; i<dims[0]; i++, volume++)
72  avg += grid(kk, j, i);
73 
74  DensityGridpoint top(0,0,kk);
75 
76  avg /= volume;
77  if (avg >= max_avg_density)
78  max_avg_density = avg;
79  }
80 
81  return(max_avg_density);
82 }
83 
84 
85 
86 int main(int argc, char *argv[]) {
87  string hdr = invocationHeader(argc, argv);
88  if (argc != 1) {
89  cerr << "Usage- gridautoscale <input.grid >output.grid\n";
90  cerr <<
91  "DESCRIPTION\n\tgridautoscale is used to normalize the density\n"
92  "values in a grid for solvated membrane systems. It divides the system\n"
93  "into bins in Z (normal to the membrane) and looks for the bulk water peak.\n"
94  "The entire grid is then scaled so that the average density in the bulk\n"
95  "regions is 1.\n";
96  exit(-1);
97  }
98 
99 
100  Grid grid;
101  cin >> grid;
102 
103  DensityGridpoint dims = grid.gridDims();
104  double best_avg = 0.0;
105  uint best_bin = 0;
106  cerr << "Autoscaling- ";
107  for (int nbins = 5; nbins <= dims[2]; ++nbins) {
108  if (nbins % 10 == 0)
109  cerr << '.';
110  double avg = findPeakDensitySlice(grid, nbins);
111  if (avg > best_avg) {
112  best_avg = avg;
113  best_bin = nbins;
114  }
115  }
116 
117  cerr << "\nScaling to 1/" << best_avg << " based on " << best_bin << " bins" << endl;
118  double konst = 1.0 / best_avg;
119  grid.scale(konst);
120  grid.addMetadata(hdr);
121  ostringstream oss;
122  oss << "Auto scaling = " << konst << ", bins = " << best_bin;
123  grid.addMetadata(oss.str());
124  cout << grid;
125 }
STL namespace.
std::string invocationHeader(int argc, char *argv[])
Create an invocation header.
Definition: utils.cpp:124
Namespace for Density package.
Definition: DensityGrid.hpp:48
Namespace for most things not already encapsulated within a class.