<jupyter_start><jupyter_text>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 ChallengeWe 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<jupyter_code>import pandas
df_hist = pandas.read_excel("soda_sales_historical_data.xlsx")
df_hist[:5]
df_hist.shape<jupyter_output><empty_output><jupyter_text>Convert categorical columns to numeric<jupyter_code>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]<jupyter_output><empty_output><jupyter_text>Picking the right predictor algorithm is of upmost importanceHence we examine our choices here in great deal.<jupyter_code>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 = []<jupyter_output><empty_output><jupyter_text>Resist the temptation to overfit!Instead, split the samples into train, test subsections.<jupyter_code>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)]<jupyter_output><empty_output><jupyter_text>By repeatedly splitting, training, and testing, you can create a realistic simulation of prediction accuracy.<jupyter_code>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"])<jupyter_output><empty_output><jupyter_text>Now lets visualize our results<jupyter_code>%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()<jupyter_output><empty_output><jupyter_text>More visualizations<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>Even more visualizations<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>Visualizations only take us so farLet's use `model_selection.cross_val_score` to automate the process of train/test split based assessment. (K-Fold Cross Validation)<jupyter_code>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")<jupyter_output><empty_output><jupyter_text>Now use the complete historical table to create a predictor object with the best algorithm<jupyter_code>fitted = (experiments["Objects"]
          [experiments["Algorithm"].index("Big Random Forest")]().
          fit(y=df_hist["Sales"], X=df_hist.drop("Sales", axis=1)))<jupyter_output><empty_output><jupyter_text>We will make predictions for the "Sales"-less table of current data<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Note that the current data table might have less categorical range than the historical data.<jupyter_code>for fld in set(df_hist.columns).difference(df_superbowl.columns, {"Sales"}):
    assert fld.startswith("dmy_")
    df_superbowl[fld] = 0<jupyter_output><empty_output><jupyter_text>**Take care!!** `sklearn` [has no concept of columns](https://github.com/scikit-learn/scikit-learn/issues/7242). We make sure that the `df_superbowl` columns are ordered consistently with the `df_hist` independent column sub-matrix.<jupyter_code>df_superbowl = df_superbowl[list(df_hist.drop("Sales", axis=1).columns)]
predicted = fitted.predict(df_superbowl)<jupyter_output><empty_output><jupyter_text>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<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Note that not all estimated discounts yield a boost in sales.<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>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}$$<jupyter_code>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<jupyter_output><empty_output><jupyter_text>$$ 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)$$<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>$$\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$$<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>Optimize and results<jupyter_code>model.setObjective(sales, sense=gu.GRB.MAXIMIZE)
model.optimize()
model.status == gu.GRB.OPTIMAL<jupyter_output><empty_output><jupyter_text>Only the paranoid surviveCarefully sanity check the solution.<jupyter_code>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"))<jupyter_output><empty_output><jupyter_text>Create a range of predictions to simulate the behavior of our solution under a range of conditions.<jupyter_code>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<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13905303313636025, 'fraction_numerical': 0.007853990575211309, 'mean_word_length': 4.188203337213815, 'pattern_counts': {'":': 6, '<': 89, '<?xml version=': 0, '>': 87, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '388329', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 4603, 'n_tokens_neox': 4382, 'n_words': 888, 'path': 'examples/expert_section/ml_soda_promotion/soda_promotion.ipynb', 'repo_name': 'adampkehoe/ticdat', 'size': '13369'}