import pandas as pd

data = pd.read_csv("./data/bryant et al 2010 data.csv", index_col=False)
x = data.iloc[:, 2:11]
y = data.iloc[:, 15].values

"""
Explanation: Performing Scenario Discovery in Python
The purpose of example is to demonstrate how one can do scenario discovery in python. 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.
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 and allowed me to share it. 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
This example is a slightly updated version of a blog post on  https://waterprogramming.wordpress.com/2015/08/05/scenario-discovery-in-python/
End of explanation
"""


from ema_workbench.analysis import prim
from ema_workbench.util import ema_logging

ema_logging.log_to_stderr(ema_logging.INFO);

"""
Explanation: 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.
End of explanation
"""


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

"""
Explanation: 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.
End of explanation
"""


box1.show_tradeoff()
plt.show()

"""
Explanation: 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.
End of explanation
"""


box1.inspect_tradeoff()

"""
Explanation: Since we are doing this analysis in a notebook, we can take advantage of the interactivity that the browser offers. A relatively recent addition to the python ecosystem is the library altair. Altair can be used to create interactive plots for use in a browser. Altair is an optional dependency for the workbench. If available, we can create the following visual.
End of explanation
"""


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

"""
Explanation: Here we can interactively explore the boxes associated with each point in the density coverage trade-off. It also offers mouse overs for the various points on the trade off curve. Given the id of each point, we can also use the workbench to manually inpect the peeling trajectory. Following Bryant & Lempert, we inspect box 21.
End of explanation
"""


box1.select(21)
fig = box1.show_pairs_scatter(21)
plt.show()

"""
Explanation: 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.
End of explanation
"""


box1.drop_restriction("Cellulosic cost")
box1.inspect(style="graph")
plt.show()

"""
Explanation: Because the last restriction is not significant, we can choose to drop this restriction from the box.
End of explanation
"""


box2 = prim_alg.find_box()

"""
Explanation: We have now found a first box that explains over 75% of the cases of interest. Let's see if we can find a second box that explains the remainder of the cases.
End of explanation
"""


prim_alg.stats_to_dataframe()

prim_alg.boxes_to_dataframe()

"""
Explanation: 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.
End of explanation
"""


from ema_workbench.analysis import cart

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

"""
Explanation: 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.
End of explanation
"""


cart_alg.stats_to_dataframe()

cart_alg.boxes_to_dataframe()

"""
Explanation: 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.
End of explanation
"""


fig = cart_alg.show_tree()
fig.set_size_inches((18, 12))
plt.show()

"""
Explanation: Alternatively, we might want to look at the classification tree directly. For this, we can use the show_tree method.
End of explanation
"""


