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

#!/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;}


  }
  }
}

}