LOOS  v2.3.2
expfit.cpp
1 /*
2 
3  This file is part of LOOS.
4 
5  LOOS (Lightweight Object-Oriented Structure library)
6  Copyright (c) 2011, Tod D. Romo
7  Department of Biochemistry and Biophysics
8  School of Medicine & Dentistry, University of Rochester
9 
10  This package (LOOS) is free software: you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation under version 3 of the License.
13 
14  This package is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  along with this program. If not, see <http://www.gnu.org/licenses/>.
21 */
22 
23 
24 
25 #include <loos.hpp>
26 #include <Simplex.hpp>
27 #include <limits>
28 
29 using namespace std;
30 using namespace loos;
31 
32 // @cond TOOLS_INTERNAL
33 
34 typedef pair<double, double> LPoint;
35 typedef vector<double> vecDouble;
36 typedef vector<vecDouble> vvecDouble;
37 
38 struct ExponentialFit {
39  ExponentialFit(const vector<LPoint>& datapoints) : _datapoints(datapoints) { }
40 
41  double operator()(const vector<double>& v) {
42  double sum = 0.0;
43 
44  for (vector<double>::const_iterator i = v.begin(); i != v.end(); ++i)
45  if (*i < 0.0)
46  return(numeric_limits<double>::max());
47 
48  for (vector<LPoint>::iterator j = _datapoints.begin(); j != _datapoints.end(); ++j) {
49  double x = j->first;
50  double y = j->second;
51 
52  double val = 1.0;
53  for (vector<double>::const_iterator i = v.begin(); i != v.end();) {
54  double k = *i++;
55  double t = *i++;
56  val += k * exp(-x/t);
57  }
58  double d = y - val;
59  sum += d*d;
60  }
61 
62  return(sum);
63  }
64 
65 
66  vector<LPoint> _datapoints;
67 };
68 
69 // @endcond
70 
71 vvecDouble readData(const string& fname) {
72  ifstream ifs(fname.c_str());
73  if (!ifs) {
74  cerr << "Error- unable to open " << fname << endl;
75  exit(-1);
76  }
77 
78  return(readTable<double>(ifs));
79 }
80 
81 
82 
83 string fullHelpMessage(void) {
84  string msg =
85  "\n"
86  "SYNOPSIS\n"
87  "\tExponential fit for bootstrapped-bcom/bcom output\n"
88  "\n"
89  "DESCRIPTION\n"
90  "\n"
91  "\tThis tool calculates a multi-exponential fit to the bootstrapped-bcom/bcom data.\n"
92  "A Nelder-Mead Simplex is used as the optimizer.\n"
93  "\n"
94  "EXAMPLES\n"
95  "\n"
96  "\texpfit bcom.asc boot_bcom.asc 5 2 0.7 10 0.3 100\n"
97  "This tries to fit bcom.asc and boot_bcom.asc using 5 replicas and using a 2-exponential\n"
98  "with inital coefficients of 0.7 and 0.3 and initial correlation times of 10 and 100\n"
99  "respectively.\n"
100  "\n"
101  "SEE ALSO\n"
102  "\tbcom, boot_bcom, bootstrap_overlap.pl\n";
103 
104  return(msg);
105 }
106 
107 
108 
109 int main(int argc, char *argv[]) {
110  string hdr = invocationHeader(argc, argv);
111  if (argc < 5) {
112  cerr << "Usage- expfit bcom.asc boot_bcom.asc nreps ndims constant-1 time-1 [constant-2 time-2 ...]\n";
113  cerr << fullHelpMessage();
114  exit(-1);
115  }
116 
117 
118  int k = 1;
119  vvecDouble bcom = readData(argv[k++]);
120  vvecDouble bbcom = readData(argv[k++]);
121 
122  if (bcom.size() != bbcom.size()) {
123  cerr << boost::format("Error- bcom has %d datapoints but bbcom has %d\n") % bcom.size() % bbcom.size();
124  exit(-10);
125  }
126 
127  vector<LPoint> datapoints;
128  for (uint i=0; i<bcom.size(); ++i)
129  datapoints.push_back(LPoint(bcom[i][0], bbcom[i][1]/bcom[i][1]));
130 
131  int nreps = atoi(argv[k++]);
132  int ndims = atoi(argv[k++]);
133  vector<double> seeds;
134  vector<double> lens;
135  ndims *= 2;
136  if (argc - 5 != ndims) {
137  cerr << boost::format("Error- only %d seeds were specified, but require %d for %d dimensions\n")
138  % (argc - 5)
139  % (ndims)
140  % (ndims/2);
141  exit(-1);
142  }
143  for (int i=0; i<ndims; ++i) {
144  double d = strtod(argv[k++], 0);
145  seeds.push_back(d);
146  lens.push_back(d/2.0);
147  }
148 
149  ExponentialFit fitFunction(datapoints);
150 
151  while (nreps-- > 0) {
152  Simplex<double> optimizer(ndims);
153  optimizer.tolerance(1e-6);
154  optimizer.maximumIterations(10000);
155  optimizer.seedLengths(lens);
156  vector<double> final = optimizer.optimize(seeds, fitFunction);
157 
158  for (uint i=0; i<final.size(); ++i)
159  cout << boost::format("%12.8g ") % final[i];
160  cout << boost::format("\t\t%.6g\n") % optimizer.finalValue();
161  seeds = final;
162  }
163 
164 }
STL namespace.
std::string invocationHeader(int argc, char *argv[])
Create an invocation header.
Definition: utils.cpp:124
Namespace for most things not already encapsulated within a class.
Nelder-Meade Simplex Optimizer (based loosely on the NRC (1996) implementation)
Definition: Simplex.hpp:58