Code Sample: Perl wrapper to run Borg with the Variable Infiltration Capacity (VIC) model

In order to perform a multi-objective calibration of VIC using Borg (on the CU Boulder supercomputer), I needed to write a wrapper, which I chose to do in Perl. Below is the wrapper. Hopefully it will help you run Borg with VIC or another hydrological model. For this scenario, I used eight variables (VIC parameters) and four objective functions (peak, length, nse, rmse). In order to run this code, you will need the following saved in your project folder: borg.exe, global.txt (e.g. global.0.txt), data_lat_long (e.g. data_39.1900_-120.2600), site_swe.txt (e.g. copper_swe.txt). Don’t forget to change all file paths in the Perl script before using it for your own scenario. If you have issues running this script, ensure the file paths in your global.0.txt file are accurate. Enjoy.

#!/usr/bin/perl
$|++; #autoflush

while (1){

##### READ IN BORG PARAMETER SETS ######

    $line = <>;
    my @array = split(/\s/,$line);
    $p1 = $array[0];
    $p2 = $array[1];
    $p3 = $array[2];
    $p4 = $array[3];
    $p5 = $array[4];
    $p6 = $array[5];
    $p7 = $array[6];
    $p8 = $array[7];

#Define filenames and paths to be used throughout the process

#output file directory
$outdir= "/lustre/janus_scratch/elho9743/VIC.borg/"; 

#soil input file template where columns with params to be changed contain "XPX1", "XPX2", etc.
$soilfiletemplate = "soil.copper.txt"; 

#output soil file used as an input for VIC (called in the global.0.txt file)
$outsoil="$outdir/soil.0.txt"; 

##### CREATE SOIL.TXT VIC INPUT FILE ######

#replace each column in the soil template file with Borg parameter values (e.g. "XP1X" in the soil file will be replaced by the $p1 value). For each model run, the soil.0.txt file will be overwritten, so only one soil.txt file will ever exist

@newpattern = ("$p1", "$p2", "$p3","$p4","$p5","$p6","$p7","$p8"); 
@pattern = ("XP1X","XP2X","XP3X","XP4X","XP5X","XP6X","XP7X","XP8X"); 

open(TEMPLATE,$soilfiletemplate) or die "cannot open soil file template $soifiletemplate";
open(SOIL,">$outsoil") or die "cannot open output soil file  $outsoil";
  foreach $line (&lt;TEMPLATE&gt;) {
    for ($j = 0; $j &lt;= $#pattern; $j++) {
      $line =~ s/$pattern[$j]/$newpattern[$j]/g; #where the magic happens (thanks Ben Livneh!)
    }
    printf SOIL "%s", $line;
  }
close(TEMPLATE);
close(SOIL);

###### RUN VIC ######

#/dev/null 2>&1 suppresses the VIC output to the command line to prevent Borg form reading it
@args = ("/lustre/janus_scratch/elho9743/VIC.4/src/vicNl -g /lustre/janus_scratch/elho9743/VIC.borg/global.0.txt > /dev/null 2>&1"); 

#implements command in terminal
system(@args); 

##### CALCULATE OBJECTIVE FUNCTIONS ######

# Open and read observed snow water equivalent (SWE) file

$filename = 'copper_swe.txt';
open (OBS,$filename) || die "Cannot open '$filename': $!";
my $file_as_string1 = join '', <OBS>;
@obs = (split /\n/, $file_as_string1);

#Separate observed SWE water years

@obs1 = @obs[0..363];
@obs2 = @obs[364..729];
@obs3 = @obs[730..1095];
@obs4 = @obs[1096..1460];
@obs5 = @obs[1461..1825];

#Open and read modeled SWE file (output from your VIC run)

$filename2 = "0_output_39.4900_-106.1700";
open (my $fh, $filename2) || die "Cannot open output";
my @file_array;
while (my $line = <$fh>) {
chomp $line;
my @line_array = split(/\s+/, $line);
push (@file_array, \@line_array);
}
@mod = map $_->[14], @file_array;

#Separate modeled SWE water years

@mod1 = @mod[6..369];
@mod2 = @mod[370..735];
@mod3 = @mod[736..1101];
@mod4 = @mod[1102..1466];
@mod5 = @mod[1467..1831];

#Calculate Objective Function 1 (peak - maximum difference between the obs. and mod. SWE peaks)

$max1 = abs( max @obs1 - max @mod1 );
$max2 = abs( max @obs2 - max @mod2 );
$max3 = abs( max @obs3 - max @mod3 );
$max4 = abs( max @obs4 - max @mod4 );
$max5 = abs( max @obs5 - max @mod5 );
$peak = max ($max1, $max2, $max3, $max4, $max5);

#Calculate Objective Function 2 (length - maximum difference between the obs. and mod. number of days with snow on the ground (SWE < 1mm))

@obs_swe1 = grep { $_ &gt; 1 } @obs1;
@obs_swe2 = grep { $_ &gt; 1 } @obs2;
@obs_swe3 = grep { $_ &gt; 1 } @obs3;
@obs_swe4 = grep { $_ &gt; 1 } @obs4;
@obs_swe5 = grep { $_ &gt; 1 } @obs5;

@mod_swe1 = grep { $_ &gt; 1 } @mod1; 
@mod_swe2 = grep { $_ &gt; 1 } @mod2; 
@mod_swe3 = grep { $_ &gt; 1 } @mod3; 
@mod_swe4 = grep { $_ &gt; 1 } @mod4; 
@mod_swe5 = grep { $_ &gt; 1 } @mod5; 

$y1 = scalar @obs_swe1; 
$y2 = scalar @obs_swe2; 
$y3 = scalar @obs_swe3; 
$y4 = scalar @obs_swe4; 
$y5 = scalar @obs_swe5; 

$x1 = scalar @mod_swe1; 
$x2 = scalar @mod_swe2; 
$x3 = scalar @mod_swe3; 
$x4 = scalar @mod_swe4; 
$x5 = scalar @mod_swe5; 

$length1 = abs( $x1 - $y1 );
$length2 = abs( $x2 - $y2 );
$length3 = abs( $x3 - $y3 );
$length4 = abs( $x4 - $y4 );
$length5 = abs( $x5 - $y5 );
$melt = max ($length1, $length2, $length3, $length4, $length5);

#Calculate Objective Functions 3 & 4 (NSE, RMSE)

@mod6 = @mod[6..1831];
$mean = sum( @obs )/@obs;  
@num1 = %obs - $mod5;
@num1 = map { $obs[$_] - $mod6[$_] } 0 .. $#obs;
@num2 = map{ $_ * $_ } @num1;
@denom1 = map { $_ - $mean } @obs;
@denom2 = map{ $_ * $_ } @denom1;
$sum1 = sum( @num2 );
$sum2 = sum( @denom2 );

$nse = -(1 - ( $sum1 / $sum2 ));
$rmse = sqrt( ($sum1) / (scalar @obs) );

##### PRINT OBJECTIVE FUNCTION VALUES TO COMMAND LINE FOR BORG TO READ #####

$line = "$peak $melt $nse $rmse\n";
print STDOUT $line

}

Below is the Linux command to run the Borg wrapper. Ensure you are in the appropriate folder (containing all files created above) when executing this command. See the Borg LRGV blog post (https://waterprogramming.wordpress.com/2014/01/23/lrgv-example-1/) for explanations on the flags in this command.

./borg.exe -v 8 -o 4 -R runtime.txt -F 100 -f solutions.txt -l 0.001,0.02,25,0.8,0.3,0.3,0.1,0.1 -u 0.03,0.15,150,0.9,0.99,0.99,0.99,0.99 -e 0.01,0.1,0.001,0.01 -n 100000 -- perl ./test.pl

Thank you to Dr. Ben Livneh and Dr. Joe Kasprzyk for assistance with writing this wrapper.

Advertisements

Compiling shared libraries on Windows (32 bit and 64 bit systems)

While playing with the Borg R-wrapper, I ran into difficulties loading the shared borg.dll library into R as a result of mingw compiling for a 32-bit system when I was using the 64 bit version of R.  After running the following command from the Windows command line,

 


gcc -shared -O3 -o borg.dll borg.c mt19937ar.c -lm

I compiled a .dll file that worked fine for the python wrapper because I have a 32 bit version of Python and the 32 bit version of the R GUI but resulted in error messages when I tried to optimize a problem using Borg in RStudio.  I was able to resolve this issue by installing MinGW-w64 and selecting the x86_64 architecture option as shown in the following screenshot.

MinGW-64-installation

Then, I added the following to my Windows path:

C:\Program Files\mingw-w64\x86_64-4.9.2-posix-seh-rt_v3-rev1\mingw64\bin

and ran the following command.  After this, everything ran smoothly with the Borg R wrapper.

x86_64-w64-mingw32-gcc -shared -O3 -o borg.dll borg.c mt19937ar.c -lm 

I found the following site useful for appropriate commands to compile using MinGW-w64.

http://fedoraproject.org/wiki/MinGW/Tutorial