#!/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 '.\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=); @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 ' \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() { #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=); $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=); 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=); 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() { 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() { 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() { 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;} } } } }