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 : |
|
|
}
|