#!/usr/bin/env python
# coding: utf-8

# # Performing Scenario Discovery in Python
# 
# The purpose of this blog post is to demonstrate how one can do scenario discovery in python. This blogpost will use the exploratory modeling workbench available on github. I will demonstrate how we can perform both PRIM in an interactive way, as well as briefly show how to use CART, which is also available in the exploratory modeling workbench. There is ample literature on both CART and PRIM and their relative merits for use in scenario discovery. So I won't be discussing that here in any detail.
# 
# The workbench is mend as a one stop shop for doing exploratory modeling, scenario discovery, and (multi-objective) robust decision making. To support this, the workbench is split into several packages. The most important packages are the expWorkbench that contains the support for setting up and executing computational experiments or (multi-objective) optimization with models; The connectors package, which contains connectors to vensim (system dynamics modeling package), netlogo (agent based modeling environment), and excel; and the analysis package that contains a wide range of techniques for visualization and analysis of the results from series of computational experiments. Here, we will focus on the analysis package. It some future blog post, I plan to demonstrate the use of the workbench for performing computational experimentation and multi-objective (robust) optimization. 
# 
# The workbench can be found on github and downloaded from there. At present, the workbench is only available for python 2.7. There is a seperate branch where I am working on making a version of the workbench that works under both python 2.7 and 3. The workbench is depended on various scientific python libraries. If you have a standard scientific python distribution, like anaconda, installed, the main dependencies will be met. In addition to the standard scientific python libraries, the workbench is also dependend on `deap` for genetic algorithms. There are also some optional dependencies. These include `seaborn` and `mpld3` for nicer and interactive visualizations, and `jpype` for controlling models implemented in Java, like netlogo, from within the workbench. 
# 
# In order to demonstrate the use of the exploratory modeling workbench for scenario discovery, I am using a published example. I am using the data used in the original article by Ben Bryant and Rob Lempert where they first introduced 2010. Ben Bryant kindly made this data available for my use. The data comes as a csv file. We can import the data easily using pandas. columns 2 up to and including 10 contain the experimental design, while the classification is presented in column 15
# 

# In[1]:


import pandas as pd

data = pd.DataFrame.from_csv('./data/bryant et al 2010 data.csv', index_col=False)
x = data.ix[:, 2:11]
y = data.ix[:, 15]


# The exploratory modeling workbench is built on top of numpy rather than pandas. This is partly a path dependecy issue. The earliest version of prim in the workbench is from 2012, when pandas was still under heavy development. Another problem is that the pandas does not contain explicit information on the datatypes of the columns. The implementation of prim in the exploratory workbench is however datatype aware, in contrast to the scenario discovery toolkit in R. That is, it will handle categorical data differently than continuous data. Internally, prim uses a numpy structured array for x, and a numpy arrea for y. We can easily transform the pandas dataframe to either. 

# In[2]:


x = x.to_records()
y = y.values


# the exploratory modeling workbench comes with a seperate analysis package. This analysis package contains prim. So let's import prim. The workbench also has its own logging functionality. We can turn this on to get some more insight into prim while it is running.

# In[3]:


from analysis import prim
from expWorkbench import ema_logging
ema_logging.log_to_stderr(ema_logging.INFO);


# Next, we need to instantiate the prim algorithm. To mimic the original work of Ben Bryant and Rob Lempert, we set the peeling alpha to 0.1. The peeling alpha determines how much data is peeled off in each iteration of the algorithm. The lower the value, the less data is removed in each iteration. The minimium coverage threshold that a box should meet is set to 0.8. Next, we can use the instantiated algorithm to find a first box. 

# In[4]:


prim_alg = prim.Prim(x, y, threshold=0.8, peel_alpha=0.1)
box1 = prim_alg.find_box()


# Let's investigate this first box is some detail. A first thing to look at is the trade off between coverage and density. The box has a convenience function for this called `show_tradeoff`. To support working in the ipython notebook, this method returns a matplotlib figure with some additional information than can be used by mpld3. 

# In[6]:


import mpld3

box1.show_tradeoff()
mpld3.display()


# As you can see, mpdl3 takes a regular matplotlib figure and transforms it into a d3 figure. In this case we use that to add pop up boxes to the tradeoff scatter plot with additional information on each point. We can also inspect each individual point in more detail. So let's look at point 21, just as in the original paper. For this, we can use the `inspect` method. By default this will display two tables, but we can also make a nice graph instead that contains the same information. 

# In[35]:


box1.inspect(21)
box1.inspect(21, style='graph')
plt.show()


# If one where to do a detailed comparison with the results reported in the original article, one would see small numerical differences. These differences arise out of subtle differences in implementation. The most important difference is that the exploratory modeling workbench uses a custom objective function inside prim which is different from the one used in the scenario discovery toolkit. Other differences have to do with details about the hill climbing optimization that is used in prim, and in particular how ties are handled in selected the next step. The differences between the two implementations are only numerical, and don't affect the overarching conclusions drawn from the analysis. 
# 
# Let's select this 21 box, and get a more detailed view of what the box looks like. Following Bryant et al., we can use scatter plots for this. 

# In[36]:


box1.select(21)
fig = box1.show_pairs_scatter()
fig.set_size_inches((12,12))
plt.show()


# We have now found a first box that explains close to 80% of the cases of interest. Let's see if we can find a second box that explains the remainder of the cases.

# In[37]:


box2 = prim_alg.find_box()


# As we can see, we are unable to find a second box. The best coverage we can achieve is 0.35, which is well below the specified 0.8 threshold. Let's look at the final overal results from interactively fitting PRIM to the data. For this, we can use to convenience functions that transform the stats and boxes to pandas data frames.

# In[38]:


print prim_alg.stats_to_dataframe()
print prim_alg.boxes_to_dataframe()


# # CART

# The way of interacting with CART is quite similar to how we setup the prim analysis. We import cart from the analysis package. We instantiate the algorithm, and next fit CART to the data. This is done via the `build_tree` method.

# In[11]:


from analysis import cart
cart_alg = cart.CART(x,y, 0.05)
cart_alg.build_tree()


# Now that we have trained CART on the data, we can investigate its results. Just like PRIM, we can use `stats_to_dataframe` and `boxes_to_dataframe` to get an overview. 

# In[13]:


print cart_alg.stats_to_dataframe()
print cart_alg.boxes_to_dataframe()


# Alternatively, we might want to look at the classification tree directly. For this, we can use the `show_tree` method. 

# In[14]:


from IPython.display import Image
Image(cart_alg.show_tree()) 


# In[ ]: