Note: This is documentation for version 5.4 of Source. For a different version of Source, select the relevant space by using the Spaces menu in the toolbar above

Calculating System Yield using Python and the Source Command Line

/wiki/spaces/SC/pages/51643422Original problem + solution: Fareed Mirza
Python solution: Joel Rahman

This example depends on these files: PythonExample.zip 

Full example is available in the source code under SourceCommunity\PythonExample.

 

Problem: We want to find the demand that the system can meet to a given reliability (eg 95%). We have a source model built (TestYieldCmdLine.rsproj) and it exposes a global expression ($gDemandDailyValue) which we can use set in order to modify the demand in the system. The model also exposes $gSupply as the total water supplied per time step, which is always equal to or less than the corresponding demand. The reliability is simply total supply divided by total demand for the simulation. The yield is the constant demand at which the system reliability crosses our target reliability.


In this exercise, we’ll calculate the yield using the Python language, a binary search and the Source external interface.


Python is a dynamic programming language that’s very popular in the scientific community. Various versions of Python are currently supported. For this exercise, I have used ‘classic’ or ‘C’ Python 3.2, available from www.python.org. Python programs can reside in module files or can be entered interactively at the Python prompt. In this case, we’ll use some helper functions I’ve defined in one module (python_helpers.py), we’ll construct another module to perform the actual binary search. We’ll then invoke the binary search program from the Python interactive prompt.


This isn’t a tutorial on Python. A few things to watch for:

  • It’s dynamically typed – so you don’t declare whether a variable is an integer or a double or a list etc…
  • Whitespace is significant. In particular, consecutive statements need to be indented consistently and the body of compound statements (eg if, while, for)
  • When entering file paths in Python you’ll use string variables (contained in quotes ‘’). You’ll need to ‘escape’ all backslashes with extra backslashes. For example, to access the program stored in C:\src\trunk\Solutions\Output\RiverSystem.Forms.exe you would use ‘C:\\src\\trunk\\Solutions\\Output\\RiverSystem.Forms.exe’

Also, note two conventions used in this tutorial:

  • C:\...> indicates commands to be entered at the windows command prompt. (Type everything to the right of the >)
  • >>> indicates commands that should be entered in the Python interactive interpreter. (Type everything to the right of the third >)

Steps

Copy all files related to the tutorial to somewhere convenient.

Start a Windows command prompt:

  • • In XP: Start|Run then enter cmd
  • • In Windows7: Start then type cmd

In the command prompt change to the directory where you’ve copied the tutorial files. For example, if you’ve copied them to c:\workshop\python

C:\...> cd \workshop\python

Then start up the Python interactive interpreter. If you installed python in C:\Python32, the command will be

C:\...> C:\Python32\Python

At the Python interactive prompt import the python_helpers module. This includes some useful methods and classes for invoking the Source command

line, passing in values of parameters and retrieving results into lists of numbers.

>>> import python_helpers

At this stage it only works for local runs and it restarts the Source command line for each run (thereby needing to load the rsproj each time). This can be improved (smile)

Importantly, the python_helpers module needs to know where to find the Source command line. If you’ve compiled Source from code and your source

code folder is c:\src\trunk, the command is:

>>> python_helpers.configure_source('C:\\src\\trunk\\TIME\\Solutions\\Output')

Open a text editor to edit a new python file (.py). This is where you will build the binary search across the demand parameter. You can use IDLE (which comes with Python) or Visual Studio (or even Notepad!)

First up this file also needs access to the python_helpers module.

import python_helpers
ph=python_helpers # alias the module because I'm lazy :)

Then we’ll write code that computes the supply reliability for a given demand. This will call the Source command line with a given value of demand, retrieve the results and then compute reliability. This function (reliability_for_demand) will be called many times from within the actual binary search, below.

import python_helpers
ph=python_helpers # alias the module :)
proj='C:\\Projects\\YieldCompWithCommandLine\\YieldProject\\TestYieldCmdLine.rsproj'
demand_param = '$gDemandDailyValue'
supply_var = '$gSupply'
output_file = 'C:\\Projects\\YieldCompWithCommandLine\\YieldProject\\py_results.csv'
def reliability_for_demand(demand):
  params=ph.Metaparams()
  params.gDemandDailyValue = demand
  
  results=ph.Results()
  results.gSupply=[]
  results.gDemandDailyValue=[]
  
  ph.run(proj,params,results)
  return compute_reliability(results)
def compute_reliability(results):
    total_demand = 0
    total_supply = 0
    
    for i in range(0,len(results.gSupply)):
        total_demand+=results.gDemandDailyValue[i]
        total_supply+=results.gSupply[i]
    return total_supply/total_demand
 
# alternatively, this function could be implemented as
# return sum(results.gSupply)/sum(results.gDemandDailyValue)

We can test the reliability_for_demand function before moving onto the actual search. First, we need to import the module (python file) you’ve just created. Assuming you’ve called it ‘find_yield.py’ and its stored in the same folder as all the other tutorial files, the statement is:

>>> import find_yield

Importantly: Don’t name the file yield.py! yield is a reserved word in python so import yield will fail!

If there are any syntax errors in your file, you’ll hear about it when you try to import it. If that all works, you can configure your code for a given project file and then run reliability_for_demand for a given demand (in this case 5)

>>> find_yield.proj = 'C:\\workshop\\python\\TestYieldCommandLine.rsproj'

>>> find_yield.reiability_for_demand(5)

If all is well, the last command above should lead to much output, some waiting (for the project to load, run and for the results to be processed) and finally a result:

>>> find_yield.reiability_for_demand(25)

0.2614598394160584

>>> 

In this case, we’re informed that the reliability was 0.26 (or 26%) – certainly less than the reliability we’re looking for!!!

Now we can create the overall yield computation using a binary search. We provide a minimal demand (that we are sure we can meet), a maximal demand (which we’re sure we CANNOT meet) and a target reliability (0..1). We first confirm that we *can* meet the minimum and cannot meet the maximum. Then we conduct a binary search, which involves testing the halfway point between min and max.If we can’t meet the reliability criteria then the yield is between the min and the halfway point, so we bisect that range. If we can meet the reliability criteria for the halfway mark, then the yield is between the halfway point and the max. We keep iterating (dividing the search space, testing midway and then moving down or up) until we reach a point where we are within 1% of meeting or not meeting the reliability criteria.

 

def compute_yield(min,max,target_reliability):
  reliability_min_demand = reliability_for_demand(min)    
  reliability_max_demand = reliability_for_demand(max)
  
  if reliability_min_demand < target_reliability:
    print( "Can't meet reliability with minimum demand")
    return
  if reliability_max_demand >= target_reliability:
    print( "Can meet reliability with maximum demand")
    return
  while True:
    trial = (max+min)/2
    print( 'Trial  = ' + str(trial))
    reliability = reliability_for_demand(trial)
    reliability_diff = reliability-target_reliability
    if abs(reliability_diff) < 0.01:
      return (trial,reliability)
    if reliability_diff < 0:
      max = trial
    else:
      min = trial

Save the file in the editor and then import into python. Then run (try a range of 1 – 20 and a reliability of 95%)

>>> from imp import reload

>>> reload(find_yield)

>>> find_yield.compute_yield(1,10,0.95)

Lots of output… then

(3.8125, 0.9500965154077076)

(In this case the yield is 3.8125 corresponding to a reliability that is just over 95%)

Other adventures with the Source External Interface

If you are familiar with another scripting language (eg bash, Ruby, Perl, R) then you will almost certainly be able to talk to the Source External Interface. The Source External Interface itself supports a wide ranging of invocation methods, including the ability to load the project once and stay alive for multiple simulation runs, step by step execution (changing parameters each timestep) and executing multiple simulations in parallel, either locally or distributed across computers.

You can reuse python_helpers in your own scripts or interactively from the python interpreter. If you want to use the python_helpers module for other problems, they key steps are as follows:

  • Always start with a call to configure_source to point the python_helpers module to your copy of the software.
  • The main function is run, with a function signature run(proj, params, results)
    • proj is a string containing the path to the .rsproj file of interest
    • params is a Metaparams object (see below) with values for each of the global expressions you want to change
    • results is a Results object (see below) initialised with the results you want to retrieve from the run

You can create a new Metaparams object as follows:

params=python_helpers.Metaparams() # metaparam objects are used to setup the values of global expressions

Once created you can set the values of the global expressions you wish to modify by setting attributes of the params object. For example, if you want to set the global expression $X to a random number between 0 and 10, you can use the following:

import random 
params.X = 10 * random.random()

Note that the call to the run function will translate this into the command to set $X to the result of the random number generator.

You can create a new Results object as follows:

results=python_helpers.Results() # Results objects gather the results of the run

Rather than setting values in the results object, you declare your interest in particular outputs by initialising them to empty lists (that are subsequently populated in the call to the run function). For results directly from global expressions you can reference them as attributes of the results object (similar to the treatment of parameters)

results.gSupply=[] # maps to a global expression $gSupply

To retrieve results from regular recorders (eg the Volume ordered from a supply point ) indicate your request by initialising a list within the ‘other’ dictionary of the results object:

results.other['Supply Point\\node #6\\Volume Ordered']=[] # maps to a global expression $gSupply

(The exact text you need for each recorded output is listed on the output when you run the Source Command Line with the project)

Once you have the project filename, a Metaparams object and a Results object, you can call run

python_helpers.run(proj,params,results) # run the model and populate the results

When this function returns, the results will have been populated and you can use a variety of python libraries. NumPy is a good starting point. You can find a list of capabilities in module using the dir() and help() functions

>>> import numpy

>>> dir(numpy)

# lots of output

>>> help(numpy.std)

std(a, axis=None, dtype=None, out=None, ddof=0)

Compute the standard deviation along the specified axis.