Annotation of /trunk/gadgetMerluza2015/model/sens.pl

Parent Directory | Revision Log

---

Revision 10 - (view) (download)

1 :	ulcessvp	10	#!/usr/bin/perl --
2 :
3 :			#Usage: sens.pl or sens.pl -i configfile
4 :			#if configfile not given the program will generate one
5 :			#
6 :			#A perl script to automate sensitivity testing of the gadget model
7 :			#Allows for user defined parameters, all controlled by switches
8 :			#in the file sens.cfg, or by a series of questions and answers
9 :			#Now allows program to continue where it left off after being stopped
10 :			#by ^C
11 :			#
12 :			#Version for gadget 2.0 (and hopefully later)
13 :			#
14 :			#Written by Daniel Howell
15 :			#24/2/2003
16 :
17 :			$version=1.0;
18 :
19 :			if(!@ARGV)
20 :			{
21 :			#if we have no commands given then we'll need to generate a config file interactively
22 :			print("No config file specified, so we'll generate one interactively\n");
23 :			generate();
24 :			readfile("$optfilename");
25 :			run();
26 :			}
27 :
28 :			elsif($ARGV[0] eq "-i")
29 :			{
30 :			#config file given, just read it in and away we go
31 :			$optfilename=$ARGV[1];
32 :			print("input file: $optfilename \n");
33 :			readfile("$optfilename");
34 :			run();
35 :			}
36 :
37 :			elsif($ARGV[0] eq "help" or $ARGV[0] eq "-help" or $ARGV[0] eq "--help")
38 :			{
39 :			#They want help, and help they shall have
40 :			print("\nThis program conducts an automated sensitivity test for a gadget model.\n");
41 :			print("To use this program you will need a config file. One was provided with the \ndistribution of this program or the program will help you generate your own.\n");
42 :			print("If you already have a config file ready then type 'sens.pl -i <configfile>'.\n");
43 :			print("Otherwise simply type 'sens.pl' and the select your options interactively.\n\n");
44 :			}
45 :			else
46 :			{
47 :			#Error! Error! Out of Cheese! Redo from Start!
48 :			print("I'm sorry I don't understand @ARGV \n");
49 :			print("Usage: sens.pl or sens.pl -i inputfile \n");
50 :			exit;
51 :			}
52 :
53 :			#Below here are functions called to generate and/or read a config file
54 :
55 :			sub generate
56 :			{
57 :			#This will generate the config file
58 :			#But we need to ask lots of questions
59 :			print("\nWe will now create a config file needed to run the program\n");
60 :			print("I am going to ask a series of questions. \nFor each one a default value will be given in [sqaure brackets]\n");
61 :			print("Pressing return will accept this default value.\n");
62 :			print("At the end of the process the information you select will be written to a configuration file, and the program run with that file\n");
63 :
64 :
65 :			print("\n\nFirst, what would you like this config file to be saved as? [sens.cfg]\n>");
66 :			$optfilename="sens.cfg";
67 :			#Get some input from the keyboard, and loose the newline
68 :			$in=getword();
69 :			#If there was some input then overwrite the default
70 :			if($in) {$optfilename=$in;}
71 :
72 :			print("What is the gadget executable called? If you haven't renamed it then it will still be gadget [gadget]\n>");
73 :			$executable="gadget";
74 :			$in=getword();
75 :			if($in) {$executable=$in;}
76 :
77 :			print("What is the name of the input file for the initial point around which the sensitivity test will be run? [refinputfile]\n>");
78 :			$inputfile="refinputfile";
79 :			$in=getword();
80 :			if($in) {$inputfile=$in;}
81 :
82 :
83 :			print("And what would you like the output file to be called? [sens.txt]\n>");
84 :			$outputfile="sens.txt";
85 :			$in=getword();
86 :			if($in) {$outputfile=$in;}
87 :
88 :
89 :			print("Next we have to specify some parameters that the program will use\n");
90 :
91 :			print("The test will be run by changing each of the variables in your model by up to +/- some percentage of the initial value\n");
92 :			print("What would you like this range to be? [50%]\n>");
93 :			$sensmax=50;
94 :			$in=getnum();
95 :			if($in)
96 :			{
97 :			$sensmax=$in;
98 :			}
99 :
100 :			print("And what stepsize would you like me to use? [5%]\n");
101 :			print("N.B. It is best to use a step size that divides exactly into your specified range.\n>");
102 :			$stepsize=5;
103 :			$in=getnum();
104 :			if($in)
105 :			{
106 :			$stepsize=$in;
107 :			}
108 :
109 :			print("Often we would like a higher resolution near the optimum than is required elsewhwere.\n");
110 :			print("You can specify an inner range which will use smaller steps than the outer range.\n");
111 :			print("Would you like to use such an inner range Y/N? [N]\n>");
112 :			$inner="N";
113 :			$in=getword();
114 :			if($in) {$inner=$in;}
115 :
116 :			if($inner eq "Y" || $inner eq "y" || $inner eq "yes" || $inner eq "YES" || $inner eq "Yes")
117 :			{
118 :			print("What size would you like this inner range to be? [5%]\n>");
119 :			$innermax=5;
120 :			$in=getnum();
121 :			if($in)
122 :			{
123 :			$innermax=$in;
124 :			}
125 :
126 :			print("And what stepsize would you like? [1%]\n");
127 :			print("N.B. It is best to use a step size that divides exactly into your specified range.\n>");
128 :			$innerstepsize=1;
129 :			$in=getnum();
130 :			if($in)
131 :			{
132 :			$innerstepsize=$in;
133 :			}
134 :			}
135 :			else
136 :			{
137 :			$innermax=$sensmax;
138 :			$innerstepsize=$stepsize;
139 :			}
140 :
141 :			print("Do you want me to only examine those variables in $inputfile that are optimized [Y/N]? [Y]\n");
142 :			print("Otherwise I will conduct a sensitivity run on the whole file, including fixed parameters.\n>");
143 :			$flag=0;
144 :			while($flag == 0)
145 :			{
146 :			$opt_only="yes";
147 :			$in=getword();
148 :			if($in)
149 :			{
150 :			$opt_only=$in;
151 :			}
152 :			if($opt_only eq "Y" || $opt_only eq "y" || $opt_only eq "yes" || $opt_only eq "YES" || $only eq "Yes")
153 :			{
154 :			$opt_only="yes";
155 :			$flag=1;
156 :			}
157 :			elsif($opt_only eq "N" || $opt_only eq "n" || $opt_only eq "no" || $opt_only eq "NO" || $opt_only eq "No")
158 :			{
159 :			$opt_only="no";
160 :			$flag=1;
161 :			}
162 :			else
163 :			{
164 :			print("Valid responses are Y,y,yes,N,n,no.\n>");
165 :			}
166 :			}
167 :
168 :			print("And of these, would you like me to examine all variables, or just a selection of them? [all]\n");
169 :			$flag=0;
170 :			while($flag == 0)
171 :			{
172 :			print("Valid options are 'all' and 'few'\n>");
173 :			$howmany="all";
174 :			$in=getword();
175 :			if($in) {$howmany=$in;}
176 :			if($in eq "all" || $in eq "ALL" || $in eq "All" || $in eq "few" | $in eq "FEW" || $in eq "Few" || $in eq "")
177 :			{
178 :			$flag=1;
179 :			}
180 :			}
181 :			if($howmany eq "FEW" || $howmany eq "Few")
182 :			{$howmany="few";}
183 :			if($howmany eq "ALL" || $howmany eq "All")
184 :			{$howmany="all";}
185 :
186 :
187 :			if($howmany eq "few")
188 :			{
189 :			print("Type the names of the variables you would like me to use, seperated by spaces.\n");
190 :			print("Press return when you are done. Sorry no wild cards allowed\n>");
191 :			#N.B. We do not use the getword function here, because in this case we can have several different
192 :			#variable names on a single line
193 :			chomp($in=<STDIN>);
194 :			@vars=$in;
195 :			}
196 :			else
197 :			{
198 :			@vars="";
199 :			}
200 :
201 :
202 :			print("\n\nOK, we're finished. The config file will be written to $optfilename.\n");
203 :
204 :			open(CONFIG, ">$optfilename") || die("couldn't create config file");
205 :			print(CONFIG ";A configuration file for the sens.pl sensitivity testing program.\n");
206 :			print(CONFIG ";This file contains all the set-up information required to perform a sensitivity test on a gadget model\n");
207 :			print(CONFIG ";Use this file by typing 'sens.pl -i <this_filename>' \n");
208 :			chomp($date=`date`);
209 :			chomp($who=`whoami`);
210 :			print(CONFIG ";This file generated automatically at $date by $who using the script sens.pl v.$version\n");
211 :			print(CONFIG "\n;sens.pl v.$version written by Daniel Howell, May 2002\n\n\n\n");
212 :
213 :			print(CONFIG ";name of the gadget executable\nexecutable = $executable \n\n");
214 :			print(CONFIG ";name of the input file with the initial point\ninputfile = $inputfile \n\n");
215 :			print(CONFIG ";name of the output file where results will be stored\noutputfile = $outputfile\n\n");
216 :			print(CONFIG ";What range will we perform our sensitivity tests across\n");
217 :			print(CONFIG "sensmax = $sensmax \nstepsize = $stepsize \ninnermax = $innermax \ninnerstep = $innerstepsize \n\n");
218 :			print(CONFIG ";Will we be looking at only the optimized variables, or all of them?\nopt_only = $opt_only\n\n");
219 :			print(CONFIG ";Will we be looking at all variables or just some? Options are 'few' and 'all'\nvars = $howmany \n\n");
220 :			print(CONFIG ";If only a few, which ones will they be? Can be blank if we are using 'all' variables\nvarnames = @vars \n\n");
221 :			print(CONFIG "\n\n;End of Config file $optfilename\n\n");
222 :
223 :			close(CONFIG);
224 :
225 :			print("\nConfig file written. Would you like to go ahead and run the sensitivity tests now? [Y]\n>");
226 :			$dummy="Y";
227 :			$in=getword();
228 :			if($in) {$dummy=$in;}
229 :			if($dummy eq "N" || $dummy eq "n" || $dummy eq "no" || $dummy eq "NO" || $dummy eq "No")
230 :			{
231 :			print("Program will abort now.\nYour Config file has been saved. To use it restart this program by typing 'sens.p -i $optfilename' . \n");
232 :			exit;
233 :			}
234 :			else
235 :			{
236 :			print("OK, program will now run.\n\n");
237 :			}
238 :			}
239 :
240 :
241 :			sub readfile
242 :			{
243 :			#will read in the config file
244 :			open(CONFIG, "$optfilename") || die "couldn't open config file $optfilename for reading\n";
245 :			#first lets stip off any spaces at the start of the line
246 :
247 :			while(<CONFIG>)
248 :			{
249 :			#Strip leading spaces
250 :			s/^\s*//;
251 :			#Make sure there are spaces between the entries on a line, ready to...
252 :			s/=/ = /;
253 :			#...split them into individual 'words'
254 :			@fields=split(/\s+/,$_);
255 :
256 :			#print("+++ $fields[0]\n");
257 :			if($fields[0])
258 :			{
259 :			#now read in the relevant parts of the file
260 :			if($fields[0] eq "executable")
261 :			{
262 :			$executable=$fields[2];
263 :			}
264 :			if($fields[0] eq "inputfile")
265 :			{
266 :			$inputfile=$fields[2];
267 :			}
268 :			if($fields[0] eq "outputfile")
269 :			{
270 :			$outputfile=$fields[2];
271 :			}
272 :			if($fields[0] eq "sensmax")
273 :			{
274 :			$sensmax=$fields[2];
275 :			$innermax=$sensmax;
276 :			}
277 :			if($fields[0] eq "stepsize")
278 :			{
279 :			$stepsize=$fields[2];
280 :			$innerstepsize=$stepsize;
281 :			}
282 :			if($fields[0] eq "innermax")
283 :			{
284 :			$innermax=$fields[2];
285 :			}
286 :			if($fields[0] eq "innerstep")
287 :			{
288 :			$innerstepsize=$fields[2];
289 :			}
290 :			if($fields[0] eq "opt_only")
291 :			{
292 :			$opt_only=$fields[2];
293 :			print("****opt_only = $opt_only\n");
294 :			}
295 :			if($fields[0] eq "vars")
296 :			{
297 :			$vars=$fields[2];
298 :			$howmany=$vars;
299 :			}
300 :			if($fields[0] eq "varnames" && $howmany eq "few")
301 :			{
302 :			@varnames=@fields;
303 :			shift(@varnames);
304 :			shift(@varnames);
305 :			#print("**** @varnames ****\n");
306 :			}
307 :			}
308 :			}
309 :
310 :			close(CONFIG);
311 :			}
312 :
313 :
314 :			#This is a little function to get a number from the standard input (keyboard)
315 :			#It will only accept a number or a return, not a string
316 :			#Because some of the numbers may be entered as percentages we will strip off any percent signs
317 :			sub getnum
318 :			{
319 :			chomp($in=<STDIN>);
320 :			$in=~ s/\%+//;
321 :			if($in && $in ne $in*1)
322 :			{
323 :			print("Please enter a non-zero number\n");
324 :			$in=getnum();
325 :			}
326 :			return $in;
327 :			}
328 :
329 :
330 :			#This function gets a word from the standard input
331 :			#It will not accept two or more words
332 :			sub getword
333 :			{
334 :			chomp($in=<STDIN>);
335 :			if($in && $in=~ /\w\s+\w/)
336 :			{
337 :			print("Please enter a single word, with no spaces\n");
338 :			$in=getword();
339 :			}
340 :			return $in;
341 :			}
342 :
343 :
344 :
345 :			#This function actually runs the sensitivity analysis once the
346 :			sub run
347 :			{
348 :
349 :			print("\n\nRunning a sensitivity test.\nYou might want to go get a coffee/go for a walk/solve world poverty while you're waiting, this could take a while\n\n");
350 :
351 :
352 :
353 :			#Set up stuff to allow program to catch a ^C
354 :			#set a flag, we haven't been zapped yet
355 :			$zapped=0;
356 :			#set up a sub routine to see if we've been zapped
357 :			sub did_he_get_us {$zapped=1;}
358 :
359 :
360 :			#First check to see if we are starting a new run, or continuing an old one
361 :
362 :			#Does the file exist that tells us our previous run was zapped?
363 :			$test="dummy";
364 :			if(-f "sens.zap")
365 :			{
366 :			print("It looks like the program stopped in the middle of a previous run\n");
367 :			print("Do you want me to do a new run, or continue the old one (n/o)?\n");
368 :			print("Typing 'n' deletes the old output file and starts a new run\n");
369 :			print("Typing 'o' continues from where the old run left off\n");
370 :			while($test ne "o" && $test ne "n")
371 :			{
372 :			chomp($test=<STDIN>);
373 :			if($test eq "n") {print("Starting a brand new run\n")}
374 :			elsif($test eq "o") {print("Continuing old run\n")}
375 :			else {print("I only understand 'o' or 'n'\n")};
376 :			}
377 :			print("\n");
378 :			}
379 :			close(STDIN);
380 :
381 :			#If there's an old file, but we are staring a new run then zap the old file
382 :			if ( -f "sens.txt" && $test ne "o") {
383 :			unlink "sens.txt";
384 :			unlink "sens.zap";
385 :			}
386 :
387 :			#If we are continuing an old run then no need to start from scratch
388 :			if($test eq "o") {
389 :			open(OLD, "sens.zap") || die "couldn't retrieve continuation data";
390 :			while(<OLD>)
391 :			{
392 :			s/^\s*//;
393 :			if(/^VAR/)
394 :			{
395 :			@varfields=split(/ +/,$_);
396 :			$minvar=$varfields[1]+1;
397 :			$oldmod=$varfields[3];
398 :			}
399 :			}
400 :			close(OLD);
401 :			unlink "sens.zap" ;
402 :			}
403 :			#open the refinput file for reading, this should contain the
404 :			#optimized parameters for the model
405 :			open(INPUT, "$inputfile") || die "couldn't open $inputfile to get optimized values";
406 :
407 :			#Now we need to know how many variables and likelihood components there are
408 :
409 :			$varcount=0;
410 :			while(<INPUT>)
411 :			{
412 :			if(/^switch/ || /^;/)
413 :			{
414 :			#these aren't variable names, so we'll do nothing here
415 :			}
416 :			else
417 :			{
418 :			@varfields=split(/\s+/,$_);
419 :			# if($opt_only eq "no" || $varfields[4] == 1)
420 :			# {
421 :			$var_ID[$varcount]=$varfields[0];
422 :			$var_value[$varcount]=$varfields[1];
423 :			$var_lower[$varcount]=$varfields[2];
424 :			$var_upper[$varcount]=$varfields[3];
425 :			$var_opt[$varcount]=$varfields[4];
426 :			@varfields=split(/;+/,$_);
427 :			$var_label[$varcount]=$varfields[1];
428 :			#print("$varcount $var_ID[$varcount]\n");
429 :			$varcount++;
430 :			# }
431 :			}
432 :			}
433 :			close(INPUT);
434 :
435 :
436 :			#if we are using all variables then set the loop bounds accordingly
437 :			if($vars eq "all" && $test ne "o")
438 :			{
439 :			$minvar=1;
440 :			$maxvar=$varcount;
441 :			}
442 :
443 :			#from here on a ^C doesn't automatically kill us
444 :			$SIG{'INT'} = 'did_he_get_us';
445 :
446 :			print("steps $mod: $stepsize, $innerstepsize, $innermax, $sensmax\n");
447 :			print("vars: $minvar, $maxvar\n");
448 :
449 :			#Now have all input values stored, we need to step through all variables
450 :			#in turn, running the sensitivity tests.
451 :			for($i=$minvar-1; $i<$maxvar; $i++)
452 :			{
453 :			if($test eq "o")
454 :			{
455 :			$test="done";
456 :			$mod=$oldmod;
457 :			}
458 :			else{$mod=-$sensmax;}
459 :
460 :			if($opt_only eq "no" || $var_opt[$i] == 1)
461 :			{
462 :			while($mod<=$sensmax)
463 :			{
464 :			$modvalue=(100+$mod)/100 ;
465 :			#print("$var_ID[$i] $mod $modvalue \n");
466 :			open(TEMPINPUT, ">tempinput") || die "couldn't create tempinput file";
467 :			print(TEMPINPUT " switch value lower upper optimize\n");
468 :			print("Running variable $i $var_ID[$i] with mod $mod \n");
469 :
470 :			#create an inputfile for this run
471 :			for($j=0;$j<$varcount;$j++)
472 :			{
473 :			$tmp_value[$j]=$var_value[$j];
474 :			if($i == $j) {$tmp_value[$j]=$var_value[$j]*$modvalue}
475 :			printf(TEMPINPUT "%s %12f %s %s %i \n", $var_ID[$j],$tmp_value[$j],$var_lower[$j],$var_upper[$j],$var_opt[$j]);
476 :			}
477 :			close(TEMPINPUT);
478 :
479 :			#run the model
480 :			{`$executable -s -i tempinput -o lik.out `;}
481 :
482 :			open(LIK, "lik.out") || die "couldn't find lik.out file";
483 :			#pull out the liklihood score
484 :			#if we have an optimizing run there will be many likelihood scores
485 :			#we want the last one
486 :			while(<LIK>)
487 :			{
488 :			if(/Lik/) {@comp=split(/\s+/,$_)};
489 :			# if(/\t/) {chomp($lik=$_);}
490 :			chomp($lik=$_);
491 :			}
492 :
493 :
494 :			#Check to see if we've been zapped
495 :			if($zapped)
496 :			{
497 :			print "\n\nArgh, I've been shot\n";
498 :			print("Data from this run has been saved in sens.txt and sens.zap\n");
499 :			print("You will be able to continue from where this run left off\n");
500 :
501 :			open(ZAPPED, ">sens.zap") || die "couldn't save temporary information";
502 :			#create a temporary file to tell the program where to start next time
503 :			print(ZAPPED "The program sens.pl was stopped in the middle of it's run\n");
504 :			print(ZAPPED "When it starts again it can carry on from where it left off\n");
505 :			print(ZAPPED "When stopped it was just about to run:\n");
506 :			print(ZAPPED "VAR: $i MOD: $mod\n");
507 :
508 :			close(ZAPPED);
509 :			#and exit since we've disabled the normal INTerrupt response
510 :			exit;
511 :			}
512 :
513 :
514 :			open(RESULTS, ">>$outputfile") || die "couldn't create $outputfile";
515 :			#print (RESULTS "$i $var_ID[$i] $modvalue @tmp_value $lik @comp[0...14] \n");
516 :			print (RESULTS "$i $var_ID[$i] $modvalue @tmp_value $lik \n");
517 :			close(RESULTS);
518 :			close(LIK);
519 :
520 :			#and increase the value by the step size
521 :			if($mod>=-$innermax && $mod<$innermax) {$mod+=$innerstepsize;}
522 :			else {$mod+=$stepsize;}
523 :
524 :
525 :			}
526 :			}
527 :			}
528 :
529 :			}

---

root@forge.cesga.es	ViewVC Help
Powered by ViewVC 1.0.0