Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Newer
Older
100644 357 lines (284 sloc) 11.881 kb
72ac9cd @schwehr cleanup and LGPL license
authored
1 // $Revision$ $Author$ $Date$
2 /*
3 Copyright (C) 2004 Kurt Schwehr
4
5 This library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public
7 License as published by the Free Software Foundation; either
8 version 2.1 of the License, or (at your option) any later version.
9
10 This library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public
16 License along with this library; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18
19 */
b155d68 @schwehr Initial revision
authored
20
21 /*
22
23 Random: double gsl_ran_gaussian (const gsl_rng * r, double sigma)
24
25 This function returns a Gaussian random variate, with mean zero and
26 standard deviation sigma. The probability distribution for Gaussian
27 random variates is,
28
29 p(x) dx = {1 \over \sqrt{2 \pi \sigma^2}} \exp (-x^2 / 2\sigma^2) dx
30
31 for x in the range -\infty to +\infty. Use the transformation z = \mu
32 + x on the numbers returned by gsl_ran_gaussian to obtain a Gaussian
33 distribution with mean \mu. This function uses the Box-Mueller
34 algorithm which requires two calls to the random number generator r.
35
36 */
37
38 #include <gsl/gsl_rng.h>
39 #include <gsl/gsl_randist.h>
40
41 #include <iostream>
20d9616 @schwehr major cleanup. Moved code to BootStrap
authored
42 #include <iomanip>
d290f0c @schwehr Starting to get Test1() going. Only does sitesigma so far.
authored
43 #include <fstream>
20d9616 @schwehr major cleanup. Moved code to BootStrap
authored
44 #include <string>
d290f0c @schwehr Starting to get Test1() going. Only does sitesigma so far.
authored
45
46 #include "kdsPmagL.H" // L is for local
47 #include "SiteSigma.H"
e6fae9e @schwehr Pulled a lot of code out to Bootstrap.[CH]
authored
48 #include "Bootstrap.H"
b155d68 @schwehr Initial revision
authored
49
50 using namespace std;
51
815f89e @schwehr Version 0.4. Now has verbosity and debug.H
authored
52
53 /***************************************************************************
54 * LOCAL MACROS and DEFINES
55 ***************************************************************************/
56
57 #include "debug.H" // provides FAILED_HERE, UNUSED, DebugPrintf
58 #ifndef NDEBUG
59 int debug_level;
60 #endif
61
62 /***************************************************************************
63 * GLOBALS
64 ***************************************************************************/
65
66 /// Let the debugger find out which version is being used.
67 static const UNUSED char* RCSid ="@(#) $Id$";
68
69 /***************************************************************************
70 * LOCAL FUNCTIONS
71 ***************************************************************************/
72
73
d290f0c @schwehr Starting to get Test1() going. Only does sitesigma so far.
authored
74 // return true if all went well.
75 // false if trouble of any kind
76
77 // Unlike Lisa's code, this one does NOT alter the sigmas on loading
78 // which is what her adread subroutine did.
e6fae9e @schwehr Pulled a lot of code out to Bootstrap.[CH]
authored
79 // Have to call SiteSigma if doing a Site based Parametric Bootstrap
d290f0c @schwehr Starting to get Test1() going. Only does sitesigma so far.
authored
80 bool
81 LoadS(const string filename,vector <SVec> &s,vector<float> &sigmas) {
82 ifstream in(filename.c_str(),ios::in);
83 if (!bool(in)) {cerr << "failed to open file: " << filename << endl; return false; }
84
85 // FIX: detect formats - {s[6]}, {s[6],sigma}, {name, sigma, s[6]}
86 // FIX: only do {s[6],sigma} for now
e6fae9e @schwehr Pulled a lot of code out to Bootstrap.[CH]
authored
87 SVec tmp(6,0.); float tmpSigma;
d290f0c @schwehr Starting to get Test1() going. Only does sitesigma so far.
authored
88 while (in >> tmp[0] >> tmp[1] >> tmp[2] >> tmp[3] >> tmp[4] >> tmp[5] >> tmpSigma) {
e6fae9e @schwehr Pulled a lot of code out to Bootstrap.[CH]
authored
89 s.push_back(tmp); sigmas.push_back(tmpSigma);
d290f0c @schwehr Starting to get Test1() going. Only does sitesigma so far.
authored
90 }
e6fae9e @schwehr Pulled a lot of code out to Bootstrap.[CH]
authored
91 // FIX: do we need to normalize so that the trace is 1?
92 // Not all data will have a trace==1??
d290f0c @schwehr Starting to get Test1() going. Only does sitesigma so far.
authored
93 return (true);
94 }
95
e6fae9e @schwehr Pulled a lot of code out to Bootstrap.[CH]
authored
96 //////////////////////////////////////////////////////////////////////
97 // MAIN
98 //////////////////////////////////////////////////////////////////////
99 #ifndef REGRESSION_TEST // NOT testing
d290f0c @schwehr Starting to get Test1() going. Only does sitesigma so far.
authored
100
c62545f @schwehr Works with the new UI.
authored
101 enum BootTypeEnum { BAD_PARAMETRIC, SITE_PARAMETRIC, SAMPLE_PARAMETRIC };
102 enum FormatEnum {BAD_FORMAT, XYZ_FORMAT,TPR_FORMAT,S_FORMAT};
103 #include "s_bootstrap_cmd.h" // Command line args
104 #include "Eigs.H" // Let's us convert to other coords
105
106 bool GetFiles(char **in_arg,vector<string> &inFiles) {
107 assert(in_arg);
108 assert(*in_arg); // Must be at least 1 file in!
109 // can not have more than 64000 input files
110 for (size_t i=0;i<64000 && in_arg[i];i++) {
111 inFiles.push_back(string(in_arg[i]));
112 }
113 return (true);
114 }
115
116 BootTypeEnum GetParametricType(const int site_given, const int sample_given) {
117 if (site_given && sample_given) {
118 cerr << "ERROR: can not specify both site and sample paramteric bootstrap" <<endl;
119 exit (EXIT_FAILURE);
120 } else if (site_given) return(SITE_PARAMETRIC);
121 return(SAMPLE_PARAMETRIC); // default
122 }
123
124
125 FormatEnum GetFormat(const char *format_arg) {
126 assert(format_arg);
127 string format(format_arg);
128 if (0==string("s").compare(format)) return(S_FORMAT);
129 if (0==string("tpr").compare(format)) return(TPR_FORMAT);
130 if (0==string("xyz").compare(format)) return(XYZ_FORMAT);
131 cerr << "ERROR: format must be one of s, tpr, or xyz" << endl;
132 exit(EXIT_FAILURE);
133 //return(BAD_FORMAT);
134 }
135
136
137 /// \param oneFile true if out1, out2, and out3 are all the same file
138 bool DoS_Bootstrap(const vector<string> &inFiles,
139 ofstream &out1, ofstream &out2, ofstream &out3,
140 const int numout_arg, const FormatEnum format, const BootTypeEnum type,
141 const int draw)
142 {
143 bool ok=true;
144 assert(1==numout_arg || 3==numout_arg);
145 assert(0<draw);
146
147 // Setup random number generator
148 gsl_rng * r; /* global generator */
149 const gsl_rng_type *T;
150 gsl_rng_env_setup();
151 T = gsl_rng_default;
152 r = gsl_rng_alloc (T);
153 { unsigned long int s; getDevRandom(s); gsl_rng_set (r, s); } // Set the Seed
154
155 S_Engine sengine; // For converting to xyz or ptr
156
157 vector<SVec> s;
158 vector<float> sigmas;
159
160 vector<float> Vmin(3,0), Vint(3,0), Vmax(3,0); // for xyz or tpr
161
162 for(size_t i=0;i<inFiles.size();i++) {
163 cout << "reading file " << inFiles[i] << endl;
164 s.clear(); sigmas.clear();
165 if (!LoadS(inFiles[i],s,sigmas)) {
166 cerr << "ERROR - can't load datafile, skipping: " << inFiles[i] << endl;
167 ok=false; continue;
168 }
169
170 const float siteSigma = (SITE_PARAMETRIC==type)?SiteSigma(s):-666.;
171
172 SVec newSample(6,0.);
173 out1 << setiosflags(ios::fixed) << setprecision(10);
174 out2 << setiosflags(ios::fixed) << setprecision(10);
175 out3 << setiosflags(ios::fixed) << setprecision(10);
176
177 // Draw out and bootstrap 'draw' number of samples
178 for (size_t i=0;i<size_t(draw); i++) {
179 switch (type) {
180 case SITE_PARAMETRIC: BootstrapParametricSite (s,siteSigma,newSample, r); break;
181 case SAMPLE_PARAMETRIC: BootstrapParametricSample (s,sigmas ,newSample, r); break;
182 default: assert(false);
183 }
184 // Now what do we do with the new sample?
185 switch(format) {
186 case S_FORMAT:
187 for (size_t i=0;i<newSample.size();i++) out1 << newSample[i] << " ";
188 break;
189 //case PTR_FORMAT: break; NOT SUPPORTED YET
190 case XYZ_FORMAT:
191 {
192 sengine.setS(newSample);
193 sengine.getXYZ(KMIN, Vmin);
194 sengine.getXYZ(KINT, Vint);
195 sengine.getXYZ(KMAX, Vmax);
196 if (1==numout_arg) {
197 for (size_t i=0;i<3;i++) out3 << Vmin[i] << " "; // V3
198 for (size_t i=0;i<3;i++) out2 << Vint[i] << " "; // V2
199 for (size_t i=0;i<3;i++) out1 << Vmax[i] << " "; // V1
200 } else {
201 for (size_t i=0;i<3;i++){out1<<Vmin[i]<<" "; out2<<Vint[i]<<" "; out3<<Vmax[i]<<" ";}
202 }
203 } // case XYZ
204 break;
205 default: assert(false && "Hell in a hand basket");
206 }
207
208 switch (numout_arg) {
209 case 1: out1 << endl; break;
210 case 3: out1 << endl; out2 << endl; out3 << endl; break;
211 default: assert(false && "What are we gonna do now, man?!?!");
212 }
213 } // for draws
214
215
216
217 } // for inFiles
218
219
220 return (ok);
221 }
222
223 //////////////////////////////////////////////////////////////////////
815f89e @schwehr Version 0.4. Now has verbosity and debug.H
authored
224 // MAIN
c62545f @schwehr Works with the new UI.
authored
225 //////////////////////////////////////////////////////////////////////
226
227 int main (const int argc, char *argv[]) {
228 bool ok=true;
b155d68 @schwehr Initial revision
authored
229
c62545f @schwehr Works with the new UI.
authored
230 gengetopt_args_info a;
231 if (0!=cmdline_parser(argc,argv,&a)) {
232 cerr << "FIX: should never get here" << endl;
233 cerr << "Early exit" << endl;
234 return (EXIT_FAILURE);
235 }
236
815f89e @schwehr Version 0.4. Now has verbosity and debug.H
authored
237 #ifdef NDEBUG
238 if (a.verbosity_given) {
239 cerr << "Verbosity is totally ignored for optimized code. Continuing in silent mode" << endl;
240 }
241 #else // debugging
242 debug_level = a.verbosity_arg;
243 DebugPrintf(TRACE,("Debug level = %d",debug_level));
244 #endif
245
c62545f @schwehr Works with the new UI.
authored
246 vector<string> inFiles;
247 if (!GetFiles(a.in_arg,inFiles)) {cerr << "Doh! What happened?" << endl; return(EXIT_FAILURE);}
248
249 if (1!=a.numout_arg && 3!=a.numout_arg) {
250 cerr << "ERROR: numout must be either 1 or 3" << endl
251 << " found: " << a.numout_arg << endl;
252 return(EXIT_FAILURE);
253 }
254
255 const FormatEnum format = GetFormat(a.format_arg);
815f89e @schwehr Version 0.4. Now has verbosity and debug.H
authored
256 DebugPrintf(VERBOSE,("format: %s (%d)",a.format_arg,int(format)));
c62545f @schwehr Works with the new UI.
authored
257
258 // Can't have S_FORMAT and 3 outfiles. Does not make sense!!!
259 if (3==a.numout_arg && S_FORMAT==format) {
260 cerr << "ERROR: you can't have 3 output files with the S_FORMAT." << endl
261 << " That would be ludicrous" << endl;
262 return (EXIT_FAILURE);
263 }
264
265
266 const BootTypeEnum type = GetParametricType(a.site_given,a.sample_given);
267 #ifndef NDEBUG
268 cerr << "Bootstrap type: " << (SITE_PARAMETRIC==type?"site":"sample") << "parametric"<<endl;
269 #endif
270
271 if (1==a.numout_arg) {
272 // just one file
273 ofstream out(a.out_arg,ios::out);
274 if (out.is_open()) {
275 if (!DoS_Bootstrap(inFiles, out,out,out, a.numout_arg, format, type, a.draw_arg)) {
276 ok=false; cerr << "ERROR: " << argv[0] << " failed in bootstrap routine." << endl;
277 }
278 } else {ok=false; cerr << "Failed to open output file" << endl;}
279 } else {
280 // 3 different files.
281 const string o1Name(string(a.out_arg)+string("1"));
282 const string o2Name(string(a.out_arg)+string("2"));
283 const string o3Name(string(a.out_arg)+string("3"));
284
285 ofstream o1(o1Name.c_str(),ios::out);
286 ofstream o2(o2Name.c_str(),ios::out);
287 ofstream o3(o3Name.c_str(),ios::out);
288
289 if (!o1.is_open() || !o2.is_open() || !o3.is_open() ) ok=false;
290 if (ok && !DoS_Bootstrap(inFiles, o1,o2,o3, a.numout_arg, format, type, a.draw_arg)) {
291 ok=false; cerr << "ERROR: " << argv[0] << " failed in bootstrap routine." << endl;
292 }
293 }
294
295 return (ok?EXIT_SUCCESS:EXIT_FAILURE);
a6ea48a @schwehr removed old UI.
authored
296 } // main
20d9616 @schwehr major cleanup. Moved code to BootStrap
authored
297
d290f0c @schwehr Starting to get Test1() going. Only does sitesigma so far.
authored
298 #endif // !REGRESSION_TEST
299
300 //////////////////////////////////////////////////////////////////////
301 // REGRESSION TESTS
302 //////////////////////////////////////////////////////////////////////
303
304 #ifdef REGRESSION_TEST
305 static bool
306 isEqual (const float a, const float b, const float del) {
307 return ( ( a<b+del && a > b-del) ? true : false );
308 }
309
310 bool Test1 (void) {
311 vector<SVec> s;
312 vector<float> sigmas;
313 if (!LoadS(string("as1-crypt.s"),s,sigmas)) {FAILED_HERE; return false;};
e6fae9e @schwehr Pulled a lot of code out to Bootstrap.[CH]
authored
314
d290f0c @schwehr Starting to get Test1() going. Only does sitesigma so far.
authored
315 const float siteSigma = SiteSigma(s);
316 const float expectedSiteSigma = 0.00133387523;
317 if (!isEqual(siteSigma, expectedSiteSigma, 0.000001)) {FAILED_HERE; return false;}
e6fae9e @schwehr Pulled a lot of code out to Bootstrap.[CH]
authored
318
319 gsl_rng *r;
320 const gsl_rng_type *T;
321 gsl_rng_env_setup();
322 T = gsl_rng_default;
323 r = gsl_rng_alloc (T);
324 // Let it always start with the same value
325 gsl_rng_set(r,0);
326 //{ unsigned long int seed; getDevRandom(seed); gsl_rng_set(r,s); }
327
328 SVec newSample;
329 const size_t siteNum = BootstrapParametricSite (s,siteSigma,newSample, r);
330 cout << "Site boot picked sample " << siteNum << endl;
331 const size_t sampleNum = BootstrapParametricSample(s,sigmas ,newSample, r);
332 cout << "Sample boot picked sample " << sampleNum << endl;
333
d290f0c @schwehr Starting to get Test1() going. Only does sitesigma so far.
authored
334 return (true);
335 }
336
e6fae9e @schwehr Pulled a lot of code out to Bootstrap.[CH]
authored
337 bool Test2 () {
338 { char a; cout << " devRandom char : " << short(getDevRandom (a)) << endl; }
339 { short a; cout << " devRandom short : " << getDevRandom (a) << endl; }
340 { int a; cout << " devRandom int : " << getDevRandom (a) << endl; }
341 { long a; cout << " devRandom long : " << getDevRandom (a) << endl; }
342 { float a; cout << " devRandom float : " << getDevRandom (a) << endl; }
343 { double a; cout << " devRandom double: " << getDevRandom (a) << endl; }
344 return (true);
345 }
d290f0c @schwehr Starting to get Test1() going. Only does sitesigma so far.
authored
346
347 int main(UNUSED int argc, char *argv[]) {
348 bool ok=true;
349
350 if (!Test1()) {FAILED_HERE;ok=false;};
e6fae9e @schwehr Pulled a lot of code out to Bootstrap.[CH]
authored
351 if (!Test2()) {FAILED_HERE;ok=false;};
d290f0c @schwehr Starting to get Test1() going. Only does sitesigma so far.
authored
352
353 cout << argv[0] << " :" << (ok?"ok":"FAILED") << endl;
354 return (ok?EXIT_SUCCESS:EXIT_FAILURE);
355 }
356 #endif // REGRESSION_TEST
Something went wrong with that request. Please try again.