#!/usr/bin/perl --
#Usage: sens.pl or sens.pl -i configfile
#if configfile not given the program will generate one
#
#A perl script to automate sensitivity testing of the gadget model
#Allows for user defined parameters, all controlled by switches
#in the file sens.cfg, or by a series of questions and answers
#Now allows program to continue where it left off after being stopped
#by ^C
#
#Version for gadget 2.0 (and hopefully later)
#
#Written by Daniel Howell
#24/2/2003
$version=1.0;
if(!@ARGV)
{
#if we have no commands given then we'll need to generate a config file interactively
print("No config file specified, so we'll generate one interactively\n");
generate();
readfile("$optfilename");
run();
}
elsif($ARGV[0] eq "-i")
{
#config file given, just read it in and away we go
$optfilename=$ARGV[1];
print("input file: $optfilename \n");
readfile("$optfilename");
run();
}
elsif($ARGV[0] eq "help" or $ARGV[0] eq "-help" or $ARGV[0] eq "--help")
{
#They want help, and help they shall have
print("\nThis program conducts an automated sensitivity test for a gadget model.\n");
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");
print("If you already have a config file ready then type 'sens.pl -i <configfile>'.\n");
print("Otherwise simply type 'sens.pl' and the select your options interactively.\n\n");
}
else
{
#Error! Error! Out of Cheese! Redo from Start!
print("I'm sorry I don't understand @ARGV \n");
print("Usage: sens.pl or sens.pl -i inputfile \n");
exit;
}
#Below here are functions called to generate and/or read a config file
sub generate
{
#This will generate the config file
#But we need to ask lots of questions
print("\nWe will now create a config file needed to run the program\n");
print("I am going to ask a series of questions. \nFor each one a default value will be given in [sqaure brackets]\n");
print("Pressing return will accept this default value.\n");
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");
print("\n\nFirst, what would you like this config file to be saved as? [sens.cfg]\n>");
$optfilename="sens.cfg";
#Get some input from the keyboard, and loose the newline
$in=getword();
#If there was some input then overwrite the default
if($in) {$optfilename=$in;}
print("What is the gadget executable called? If you haven't renamed it then it will still be gadget [gadget]\n>");
$executable="gadget";
$in=getword();
if($in) {$executable=$in;}
print("What is the name of the input file for the initial point around which the sensitivity test will be run? [refinputfile]\n>");
$inputfile="refinputfile";
$in=getword();
if($in) {$inputfile=$in;}
print("And what would you like the output file to be called? [sens.txt]\n>");
$outputfile="sens.txt";
$in=getword();
if($in) {$outputfile=$in;}
print("Next we have to specify some parameters that the program will use\n");
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");
print("What would you like this range to be? [50%]\n>");
$sensmax=50;
$in=getnum();
if($in)
{
$sensmax=$in;
}
print("And what stepsize would you like me to use? [5%]\n");
print("N.B. It is best to use a step size that divides exactly into your specified range.\n>");
$stepsize=5;
$in=getnum();
if($in)
{
$stepsize=$in;
}
print("Often we would like a higher resolution near the optimum than is required elsewhwere.\n");
print("You can specify an inner range which will use smaller steps than the outer range.\n");
print("Would you like to use such an inner range Y/N? [N]\n>");
$inner="N";
$in=getword();
if($in) {$inner=$in;}
if($inner eq "Y" || $inner eq "y" || $inner eq "yes" || $inner eq "YES" || $inner eq "Yes")
{
print("What size would you like this inner range to be? [5%]\n>");
$innermax=5;
$in=getnum();
if($in)
{
$innermax=$in;
}
print("And what stepsize would you like? [1%]\n");
print("N.B. It is best to use a step size that divides exactly into your specified range.\n>");
$innerstepsize=1;
$in=getnum();
if($in)
{
$innerstepsize=$in;
}
}
else
{
$innermax=$sensmax;
$innerstepsize=$stepsize;
}
print("Do you want me to only examine those variables in $inputfile that are optimized [Y/N]? [Y]\n");
print("Otherwise I will conduct a sensitivity run on the whole file, including fixed parameters.\n>");
$flag=0;
while($flag == 0)
{
$opt_only="yes";
$in=getword();
if($in)
{
$opt_only=$in;
}
if($opt_only eq "Y" || $opt_only eq "y" || $opt_only eq "yes" || $opt_only eq "YES" || $only eq "Yes")
{
$opt_only="yes";
$flag=1;
}
elsif($opt_only eq "N" || $opt_only eq "n" || $opt_only eq "no" || $opt_only eq "NO" || $opt_only eq "No")
{
$opt_only="no";
$flag=1;
}
else
{
print("Valid responses are Y,y,yes,N,n,no.\n>");
}
}
print("And of these, would you like me to examine all variables, or just a selection of them? [all]\n");
$flag=0;
while($flag == 0)
{
print("Valid options are 'all' and 'few'\n>");
$howmany="all";
$in=getword();
if($in) {$howmany=$in;}
if($in eq "all" || $in eq "ALL" || $in eq "All" || $in eq "few" | $in eq "FEW" || $in eq "Few" || $in eq "")
{
$flag=1;
}
}
if($howmany eq "FEW" || $howmany eq "Few")
{$howmany="few";}
if($howmany eq "ALL" || $howmany eq "All")
{$howmany="all";}
if($howmany eq "few")
{
print("Type the names of the variables you would like me to use, seperated by spaces.\n");
print("Press return when you are done. Sorry no wild cards allowed\n>");
#N.B. We do not use the getword function here, because in this case we can have several different
#variable names on a single line
chomp($in=<STDIN>);
@vars=$in;
}
else
{
@vars="";
}
print("\n\nOK, we're finished. The config file will be written to $optfilename.\n");
open(CONFIG, ">$optfilename") || die("couldn't create config file");
print(CONFIG ";A configuration file for the sens.pl sensitivity testing program.\n");
print(CONFIG ";This file contains all the set-up information required to perform a sensitivity test on a gadget model\n");
print(CONFIG ";Use this file by typing 'sens.pl -i <this_filename>' \n");
chomp($date=`date`);
chomp($who=`whoami`);
print(CONFIG ";This file generated automatically at $date by $who using the script sens.pl v.$version\n");
print(CONFIG "\n;sens.pl v.$version written by Daniel Howell, May 2002\n\n\n\n");
print(CONFIG ";name of the gadget executable\nexecutable = $executable \n\n");
print(CONFIG ";name of the input file with the initial point\ninputfile = $inputfile \n\n");
print(CONFIG ";name of the output file where results will be stored\noutputfile = $outputfile\n\n");
print(CONFIG ";What range will we perform our sensitivity tests across\n");
print(CONFIG "sensmax = $sensmax \nstepsize = $stepsize \ninnermax = $innermax \ninnerstep = $innerstepsize \n\n");
print(CONFIG ";Will we be looking at only the optimized variables, or all of them?\nopt_only = $opt_only\n\n");
print(CONFIG ";Will we be looking at all variables or just some? Options are 'few' and 'all'\nvars = $howmany \n\n");
print(CONFIG ";If only a few, which ones will they be? Can be blank if we are using 'all' variables\nvarnames = @vars \n\n");
print(CONFIG "\n\n;End of Config file $optfilename\n\n");
close(CONFIG);
print("\nConfig file written. Would you like to go ahead and run the sensitivity tests now? [Y]\n>");
$dummy="Y";
$in=getword();
if($in) {$dummy=$in;}
if($dummy eq "N" || $dummy eq "n" || $dummy eq "no" || $dummy eq "NO" || $dummy eq "No")
{
print("Program will abort now.\nYour Config file has been saved. To use it restart this program by typing 'sens.p -i $optfilename' . \n");
exit;
}
else
{
print("OK, program will now run.\n\n");
}
}
sub readfile
{
#will read in the config file
open(CONFIG, "$optfilename") || die "couldn't open config file $optfilename for reading\n";
#first lets stip off any spaces at the start of the line
while(<CONFIG>)
{
#Strip leading spaces
s/^\s*//;
#Make sure there are spaces between the entries on a line, ready to...
s/=/ = /;
#...split them into individual 'words'
@fields=split(/\s+/,$_);
#print("+++ $fields[0]\n");
if($fields[0])
{
#now read in the relevant parts of the file
if($fields[0] eq "executable")
{
$executable=$fields[2];
}
if($fields[0] eq "inputfile")
{
$inputfile=$fields[2];
}
if($fields[0] eq "outputfile")
{
$outputfile=$fields[2];
}
if($fields[0] eq "sensmax")
{
$sensmax=$fields[2];
$innermax=$sensmax;
}
if($fields[0] eq "stepsize")
{
$stepsize=$fields[2];
$innerstepsize=$stepsize;
}
if($fields[0] eq "innermax")
{
$innermax=$fields[2];
}
if($fields[0] eq "innerstep")
{
$innerstepsize=$fields[2];
}
if($fields[0] eq "opt_only")
{
$opt_only=$fields[2];
print("****opt_only = $opt_only\n");
}
if($fields[0] eq "vars")
{
$vars=$fields[2];
$howmany=$vars;
}
if($fields[0] eq "varnames" && $howmany eq "few")
{
@varnames=@fields;
shift(@varnames);
shift(@varnames);
#print("**** @varnames ****\n");
}
}
}
close(CONFIG);
}
#This is a little function to get a number from the standard input (keyboard)
#It will only accept a number or a return, not a string
#Because some of the numbers may be entered as percentages we will strip off any percent signs
sub getnum
{
chomp($in=<STDIN>);
$in=~ s/\%+//;
if($in && $in ne $in*1)
{
print("Please enter a non-zero number\n");
$in=getnum();
}
return $in;
}
#This function gets a word from the standard input
#It will not accept two or more words
sub getword
{
chomp($in=<STDIN>);
if($in && $in=~ /\w\s+\w/)
{
print("Please enter a single word, with no spaces\n");
$in=getword();
}
return $in;
}
#This function actually runs the sensitivity analysis once the
sub run
{
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");
#Set up stuff to allow program to catch a ^C
#set a flag, we haven't been zapped yet
$zapped=0;
#set up a sub routine to see if we've been zapped
sub did_he_get_us {$zapped=1;}
#First check to see if we are starting a new run, or continuing an old one
#Does the file exist that tells us our previous run was zapped?
$test="dummy";
if(-f "sens.zap")
{
print("It looks like the program stopped in the middle of a previous run\n");
print("Do you want me to do a new run, or continue the old one (n/o)?\n");
print("Typing 'n' deletes the old output file and starts a new run\n");
print("Typing 'o' continues from where the old run left off\n");
while($test ne "o" && $test ne "n")
{
chomp($test=<STDIN>);
if($test eq "n") {print("Starting a brand new run\n")}
elsif($test eq "o") {print("Continuing old run\n")}
else {print("I only understand 'o' or 'n'\n")};
}
print("\n");
}
close(STDIN);
#If there's an old file, but we are staring a new run then zap the old file
if ( -f "sens.txt" && $test ne "o") {
unlink "sens.txt";
unlink "sens.zap";
}
#If we are continuing an old run then no need to start from scratch
if($test eq "o") {
open(OLD, "sens.zap") || die "couldn't retrieve continuation data";
while(<OLD>)
{
s/^\s*//;
if(/^VAR/)
{
@varfields=split(/ +/,$_);
$minvar=$varfields[1]+1;
$oldmod=$varfields[3];
}
}
close(OLD);
unlink "sens.zap" ;
}
#open the refinput file for reading, this should contain the
#optimized parameters for the model
open(INPUT, "$inputfile") || die "couldn't open $inputfile to get optimized values";
#Now we need to know how many variables and likelihood components there are
$varcount=0;
while(<INPUT>)
{
if(/^switch/ || /^;/)
{
#these aren't variable names, so we'll do nothing here
}
else
{
@varfields=split(/\s+/,$_);
# if($opt_only eq "no" || $varfields[4] == 1)
# {
$var_ID[$varcount]=$varfields[0];
$var_value[$varcount]=$varfields[1];
$var_lower[$varcount]=$varfields[2];
$var_upper[$varcount]=$varfields[3];
$var_opt[$varcount]=$varfields[4];
@varfields=split(/;+/,$_);
$var_label[$varcount]=$varfields[1];
#print("$varcount $var_ID[$varcount]\n");
$varcount++;
# }
}
}
close(INPUT);
#if we are using all variables then set the loop bounds accordingly
if($vars eq "all" && $test ne "o")
{
$minvar=1;
$maxvar=$varcount;
}
#from here on a ^C doesn't automatically kill us
$SIG{'INT'} = 'did_he_get_us';
print("steps $mod: $stepsize, $innerstepsize, $innermax, $sensmax\n");
print("vars: $minvar, $maxvar\n");
#Now have all input values stored, we need to step through all variables
#in turn, running the sensitivity tests.
for($i=$minvar-1; $i<$maxvar; $i++)
{
if($test eq "o")
{
$test="done";
$mod=$oldmod;
}
else{$mod=-$sensmax;}
if($opt_only eq "no" || $var_opt[$i] == 1)
{
while($mod<=$sensmax)
{
$modvalue=(100+$mod)/100 ;
#print("$var_ID[$i] $mod $modvalue \n");
open(TEMPINPUT, ">tempinput") || die "couldn't create tempinput file";
print(TEMPINPUT " switch value lower upper optimize\n");
print("Running variable $i $var_ID[$i] with mod $mod \n");
#create an inputfile for this run
for($j=0;$j<$varcount;$j++)
{
$tmp_value[$j]=$var_value[$j];
if($i == $j) {$tmp_value[$j]=$var_value[$j]*$modvalue}
printf(TEMPINPUT "%s %12f %s %s %i \n", $var_ID[$j],$tmp_value[$j],$var_lower[$j],$var_upper[$j],$var_opt[$j]);
}
close(TEMPINPUT);
#run the model
{`$executable -s -i tempinput -o lik.out `;}
open(LIK, "lik.out") || die "couldn't find lik.out file";
#pull out the liklihood score
#if we have an optimizing run there will be many likelihood scores
#we want the last one
while(<LIK>)
{
if(/Lik/) {@comp=split(/\s+/,$_)};
# if(/\t/) {chomp($lik=$_);}
chomp($lik=$_);
}
#Check to see if we've been zapped
if($zapped)
{
print "\n\nArgh, I've been shot\n";
print("Data from this run has been saved in sens.txt and sens.zap\n");
print("You will be able to continue from where this run left off\n");
open(ZAPPED, ">sens.zap") || die "couldn't save temporary information";
#create a temporary file to tell the program where to start next time
print(ZAPPED "The program sens.pl was stopped in the middle of it's run\n");
print(ZAPPED "When it starts again it can carry on from where it left off\n");
print(ZAPPED "When stopped it was just about to run:\n");
print(ZAPPED "VAR: $i MOD: $mod\n");
close(ZAPPED);
#and exit since we've disabled the normal INTerrupt response
exit;
}
open(RESULTS, ">>$outputfile") || die "couldn't create $outputfile";
#print (RESULTS "$i $var_ID[$i] $modvalue @tmp_value $lik @comp[0...14] \n");
print (RESULTS "$i $var_ID[$i] $modvalue @tmp_value $lik \n");
close(RESULTS);
close(LIK);
#and increase the value by the step size
if($mod>=-$innermax && $mod<$innermax) {$mod+=$innerstepsize;}
else {$mod+=$stepsize;}
}
}
}
}