<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Combining Machine Learning and Optimization
With Gurobi and sklearn
Machine Learning topics
Touching the elephant here, but ~~not there~~
Supervised Learning
  * Algorithm selection and hyper-parametric optimization
  * KFold assessment vs overfitting
  * Separating training from prediction 
~~Unsupervised Learning~~
~~Time Series Data~~
~~Deep Learning~~
Optimization topics
Exploratory programming to application deployment
Coping with the combinatorial explosion
Validating optimization with simulation
All Under the Banner of Python!
The Soda Promotion Problem
We have the challenge of designing the upcoming promotion campaign for a Soda Company. The intended objective is to bolster sales while at the same time obeying various business constraints.
The First Challenge
We need to predict impact of different price points on the expected sales for each type of soda.
To do this, we need to train a soda sales predictor from a historical data table.
Examine historical data
Step1:  Convert categorical columns to numeric
Step2:  Picking the right predictor algorithm is of upmost importance
Hence we examine our choices here in great deal.
Step3:  Resist the temptation to overfit!
Instead, split the samples into train, test subsections.
Step4:  By repeatedly splitting, training, and testing, you can create a realistic simulation of prediction accuracy.
Step5:  Now lets visualize our results
Step6:  More visualizations
Step7:  Even more visualizations
Step8:  Visualizations only take us so far
Let's use model_selection.cross_val_score to automate the process of train/test split based assessment. (K-Fold Cross Validation)
Step9:  Now use the complete historical table to create a predictor object with the best algorithm
Step10:  We will make predictions for the "Sales"-less table of current data
Step11:  Note that the current data table might have less categorical range than the historical data.
Step12:  Take care!! sklearn has no concept of columns. We make sure that the df_superbowl columns are ordered consistently with the df_hist independent column sub-matrix.
Step13:  Prediction in hand, we commence optimization!
LaTeX summary of family of equations
$$
\begin{array}{ll}
\max & sales\
s.t. & X_{b,p}\in{0,1}\quad\forall (b,p)\in Prod\
     & \sum\left(X_{b,p}
Step14:  Note that not all estimated discounts yield a boost in sales.
Step15:  Building a MIP model
$$
\begin{array}{l}
X_{b,p}\in{0,1}\quad\forall (b,p)\in Prod\
0 \leq sales\
0 \leq revenue\
0 \leq investment \leq max_{investment}
\end{array}
$$
Step16:  $$ 
sales = \sum\left(f_{b,p} X_{b,p}
Step17:  $$
\sum\left(X_{b,p}
Step18:  Optimize and results
Step19:  Only the paranoid survive
Carefully sanity check the solution.
Step20:  Create a range of predictions to simulate the behavior of our solution under a range of conditions.
<ASSISTANT_TASK:>
Python Code:
import pandas
df_hist = pandas.read_excel("soda_sales_historical_data.xlsx")
df_hist[:5]
df_hist.shape
Explanation: Combining Machine Learning and Optimization
With Gurobi and sklearn
Machine Learning topics
Touching the elephant here, but ~~not there~~
Supervised Learning
  * Algorithm selection and hyper-parametric optimization
  * KFold assessment vs overfitting
  * Separating training from prediction 
~~Unsupervised Learning~~
~~Time Series Data~~
~~Deep Learning~~
Optimization topics
Exploratory programming to application deployment
Coping with the combinatorial explosion
Validating optimization with simulation
All Under the Banner of Python!
The Soda Promotion Problem
We have the challenge of designing the upcoming promotion campaign for a Soda Company. The intended objective is to bolster sales while at the same time obeying various business constraints.
The First Challenge
We need to predict impact of different price points on the expected sales for each type of soda.
To do this, we need to train a soda sales predictor from a historical data table.
Examine historical data
End of explanation
from pandas import DataFrame, get_dummies
categorical_columns = ['Product','Easter Included','Super Bowl Included', 
                       'Christmas Included', 'Other Holiday']
df_hist = get_dummies(df_hist, prefix={k:"dmy_%s"%k for k in categorical_columns},
                      columns = list(categorical_columns))
df_hist[:5]
Explanation: Convert categorical columns to numeric
End of explanation
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import BaggingRegressor
from sklearn import model_selection
experiments = {"Algorithm":["Ordinary Least Squares", "Regression Tree", 
                            "Big Random Forest", "Random Forest", 
                            "Bagging"], 
               "Objects" : [lambda : LinearRegression(), 
                            lambda : DecisionTreeRegressor(), 
                            lambda : RandomForestRegressor(n_estimators=100), 
                            lambda : RandomForestRegressor(), 
                            lambda : BaggingRegressor()], 
               "Predictions":[[] for _ in range(5)]}
actuals = []
Explanation: Picking the right predictor algorithm is of upmost importance
Hence we examine our choices here in great deal.
End of explanation
from sklearn.model_selection import train_test_split
[_.shape for _ in  train_test_split(df_hist.drop("Sales", axis=1), 
                                    df_hist["Sales"], test_size=0.25)]
Explanation: Resist the temptation to overfit!
Instead, split the samples into train, test subsections.
End of explanation
for _ in range (4):
    train_X, test_X, train_y, test_y = (
        train_test_split(df_hist.drop("Sales", axis=1), 
                         df_hist["Sales"], test_size=0.25))
    for i, obj_factory in enumerate(experiments["Objects"]):
        obj = obj_factory()
        obj.fit(y=train_y,X=train_X)
        experiments["Predictions"][i] += list(obj.predict(test_X))
    actuals += list(test_y)
actuals = pandas.Series(actuals)
experiments["Predictions"] = list(map(pandas.Series, experiments["Predictions"]))
len(actuals), map(len, experiments["Predictions"])
Explanation: By repeatedly splitting, training, and testing, you can create a realistic simulation of prediction accuracy.
End of explanation
%matplotlib inline 
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import cm
color=iter(cm.rainbow(np.linspace(0,1,len(experiments)+3)))
plt.figure(figsize=(12,7),dpi=300)
plt.plot(actuals,actuals,c=next(color),markersize=2,label='Data')
for _, row in DataFrame(experiments).iterrows():
    plt.plot(actuals, row["Predictions"],'o',c=next(color),
             markersize=2,label=row['Algorithm'])
plt.title('Scatter Plot Prediction v/s Data')
plt.grid(True)
plt.legend()
plt.show()
Explanation: Now lets visualize our results
End of explanation
color=iter(cm.rainbow(np.linspace(0,1,len(experiments)+3)))
next(color)
plt.figure(figsize=(13,8),dpi=300)
for index, row in DataFrame(experiments).iterrows():
    relative_error = (row["Predictions"] - actuals) / (1 + abs(actuals))
    plt.plot(np.sort(relative_error),'o',c=next(color),
             markersize=2,label=row['Algorithm'])
plt.title('Relative Error Prediction v/s Data')
plt.ylabel('Relative Error')
plt.grid(True)
plt.legend()
plt.axis([0,len(actuals),-1,1])
plt.show()
Explanation: More visualizations
End of explanation
def boxplot(algorithm):
    prediction = (experiments["Predictions"]
                  [experiments["Algorithm"].index(algorithm)])
    plt.title(algorithm)
    plt.boxplot( (prediction - actuals) / (1 + abs(actuals))  )
    plt.show()
boxplot("Bagging")
boxplot("Big Random Forest")
Explanation: Even more visualizations
End of explanation
experiments["Results"] = []
for o in experiments["Objects"]:
    experiments["Results"].append(
        model_selection.cross_val_score(o(), y=df_hist['Sales'], 
                                        X=df_hist.drop("Sales", axis=1),
                                        cv=5).mean())
DataFrame(experiments).drop(["Objects", "Predictions"], 
                            axis=1).set_index("Algorithm")
Explanation: Visualizations only take us so far
Let's use model_selection.cross_val_score to automate the process of train/test split based assessment. (K-Fold Cross Validation)
End of explanation
fitted = (experiments["Objects"]
          [experiments["Algorithm"].index("Big Random Forest")]().
          fit(y=df_hist["Sales"], X=df_hist.drop("Sales", axis=1)))
Explanation: Now use the complete historical table to create a predictor object with the best algorithm
End of explanation
df_superbowl_original = pandas.read_excel("super_bowl_promotion_data.xlsx")
df_superbowl = get_dummies(df_superbowl_original, 
                           prefix={k:"dmy_%s"%k for k in categorical_columns},
                           columns = list(categorical_columns))
assert "Sales" not in df_superbowl.columns 
assert {"Sales"}.union(df_superbowl.columns).issubset(set(df_hist.columns))
len(df_superbowl)
Explanation: We will make predictions for the "Sales"-less table of current data
End of explanation
for fld in set(df_hist.columns).difference(df_superbowl.columns, {"Sales"}):
    assert fld.startswith("dmy_")
    df_superbowl[fld] = 0
Explanation: Note that the current data table might have less categorical range than the historical data.
End of explanation
df_superbowl = df_superbowl[list(df_hist.drop("Sales", axis=1).columns)]
predicted = fitted.predict(df_superbowl)
Explanation: Take care!! sklearn has no concept of columns. We make sure that the df_superbowl columns are ordered consistently with the df_hist independent column sub-matrix.
End of explanation
forecast_sales = df_superbowl_original[["Product", "Cost Per Unit"]].copy()
forecast_sales["Sales"] = predicted
forecast_sales.set_index(['Product','Cost Per Unit'], inplace=True)
soda_family = {'11 Down': 'Clear', 'AB Root Beer': 'Dark', 
               'Alpine Stream': 'Clear', 'Bright': 'Clear', 
               'Crisp Clear': 'Clear', 'DC Kola': 'Dark',
               'Koala Kola': 'Dark', 'Mr. Popper': 'Dark', 
               'Popsi Kola': 'Dark'}
family  = set(soda_family[j] for j in soda_family)
soda    = set(j for j in soda_family)
max_prom = {f:2 for f in family}
max_investment = 750
product_prices = set(forecast_sales.index.values)
normal_price = {b:0 for b in soda}
for b,p in product_prices:
    normal_price[b] = max(normal_price[b],p)
Explanation: Prediction in hand, we commence optimization!
LaTeX summary of family of equations
$$
\begin{array}{ll}
\max & sales\
s.t. & X_{b,p}\in{0,1}\quad\forall (b,p)\in Prod\
     & \sum\left(X_{b,p}:{(b,p)\in Prod} \right)=1\quad\forall b\in Soda\
     & \sum\left(X_{b,p}:{(b,p)\in Prod, p\neq p_o, T(b)=t} \right) \leq max_t\quad\forall t\
     & sales = \sum\left(f_{b,p} X_{b,p}:{(b,p)\in Prod}\right)\
     & revenue = \sum\left(f_{b,p}  p X_{b,p}:{(b,p)\in Prod} \right)\
     & investment = \sum\left(\left(f_{b,p} - f_{b,p_o}\right)+ p_o X{b,p}:{(b,p)\in Prod} \right)\
     & investment <= max_{investment}
\end{array}
$$
Putting the optimization input set together
End of explanation
meaningful_discounts = 0
for b,p in product_prices:
    if forecast_sales.Sales[b,p] > forecast_sales.Sales[b,normal_price[b]]:
        meaningful_discounts += 1
meaningful_discounts, len(forecast_sales) - len(soda)
Explanation: Note that not all estimated discounts yield a boost in sales.
End of explanation
import gurobipy as gu
model = gu.Model()
select_price = model.addVars(product_prices,vtype=gu.GRB.BINARY,name='X')
sales        = model.addVar(name='sales')
revenue      = model.addVar(name='revenue')
investment   = model.addVar(ub=max_investment, name='investment')
gusum = gu.quicksum
Explanation: Building a MIP model
$$
\begin{array}{l}
X_{b,p}\in{0,1}\quad\forall (b,p)\in Prod\
0 \leq sales\
0 \leq revenue\
0 \leq investment \leq max_{investment}
\end{array}
$$
End of explanation
model.addConstr(sales == select_price.prod(forecast_sales.Sales), name='sales')
model.addConstr(revenue == gusum(forecast_sales.Sales[b,p] * p * 
                                 select_price[b,p] for b,p in product_prices), 
                name='revenue')
model.addConstr(investment == 
                gusum(max(0,forecast_sales.Sales[b,p] - 
                            forecast_sales.Sales[b,normal_price[b]]) *
                        normal_price[b] * select_price[b,p] 
                        for b,p in product_prices),
                name='investment')
model.update()
Explanation: $$ 
sales = \sum\left(f_{b,p} X_{b,p}:{(b,p)\in Prod}\right)\
revenue = \sum\left(f_{b,p} p X_{b,p}:{(b,p)\in Prod} \right)\
investment = \sum\left(\left(f_{b,p} - f_{b,p_o}\right)+ p_o X{b,p}:{(b,p)\in Prod} \right)
$$
End of explanation
model.addConstrs((select_price.sum(b,'*') == 1 for b in soda), name='OnePrice')
model.addConstrs((gusum(select_price[b,p] for b,p in product_prices if 
                        soda_family[b] == f and p != normal_price[b] ) 
                  <= max_prom[f] for f in family),
                 name='MaxProm')
model.update()
Explanation: $$
\sum\left(X_{b,p}:{(b,p)\in Prod} \right)=1\quad\forall b\in Soda\
\sum\left(X_{b,p}:{(b,p)\in Prod, p\neq p_o, T(b)=t} \right) \leq max_t
$$
End of explanation
model.setObjective(sales, sense=gu.GRB.MAXIMIZE)
model.optimize()
model.status == gu.GRB.OPTIMAL
Explanation: Optimize and results
End of explanation
sales.X, revenue.X, investment.X
price_selections = {"Product":[], "Price":[], "Is Discount":[], "Family":[]}
for b, p in product_prices:
    if abs(select_price[b,p].X -1) < 0.0001: # i.e. almost one
        price_selections["Product"].append(b)
        price_selections["Price"].append(p)
        price_selections["Is Discount"].append(p < normal_price[b])
        price_selections["Family"].append(soda_family[b])
(DataFrame(price_selections).set_index("Product")
 [["Price", "Is Discount", "Family"]].sort_values("Family"))
Explanation: Only the paranoid survive
Carefully sanity check the solution.
End of explanation
simulated_KPI = {'Sales':[],'Revenue':[],'Investment':[]}
Z = select_price
num_infeas = 0
for i in range(100):
    np.random.seed(i)
    fitted = RandomForestRegressor(n_estimators=100,
                                   n_jobs=4).fit(y=df_hist["Sales"],
                                                 X=df_hist.drop("Sales", axis=1))
    forecast = df_superbowl_original[['Product', 'Cost Per Unit']].copy()
    forecast["Sales"] = fitted.predict(df_superbowl)
    forecast = forecast.set_index(['Product','Cost Per Unit'])
    sales, revenue, investment = 0, 0, 0
    for b,p in product_prices:
        sales   += forecast.Sales[b,p] * Z[b,p].X
        revenue += forecast.Sales[b,p] * p * Z[b,p].X
        investment += (max(0,forecast.Sales[b,p] - 
                             forecast.Sales[b,normal_price[b]]) * 
                       normal_price[b] * Z[b,p].X)
    if investment > max_investment:
        num_infeas += 1
    simulated_KPI['Sales'].append(sales)
    simulated_KPI['Revenue'].append(revenue)
    simulated_KPI['Investment'].append(investment)
data = {'Sales','Revenue','Investment'}
color=iter(cm.rainbow(np.linspace(0,1,3)))
for t in data:
    plt.figure(figsize=(7,4),dpi=300)
    plt.hist(simulated_KPI[t],50,normed=1,color=next(color), alpha=0.75)
    plt.ylabel('Probability')
    plt.xlabel(t)
    plt.grid(True)
    plt.show()
num_infeas
Explanation: Create a range of predictions to simulate the behavior of our solution under a range of conditions.
End of explanation
<END_TASK>