<jupyter_start><jupyter_text>Design of Experiments in Python<NAME>    https://github.com/hpanderson/dexpy-pymntos    Agenda* What is Design of Experiments?* Design an experiment with `dexpy` to improve the office coffee* Measure the power of our experiment* Model the taste test results using `statsmodels`* Visualize the data with `seaborn` and `matplotlib` to find the best pot of coffee  What is Design of Experiments (DoE)?A systematic series of tests, in which purposeful changes are made to input factors, so that you may identify causes for significant changes in the output responses.NIST has a [nice primer on DOE](http://www.itl.nist.gov/div898/handbook/pri/section1/pri1.htm).<jupyter_code>dot = Digraph(comment='Design of Experiments')
dot.body.extend(['rankdir=LR', 'size="10,10"'])
dot.node_attr.update(shape='rectangle', style='filled', fontsize='20', fontname="helvetica")

dot.node('X', 'Controllable Factors', color='mediumseagreen', width='3')
dot.node('Z', 'Noise Factors', color='indianred2', width='3')
dot.node('P', 'Process', color='lightblue', height='1.25', width='3')
dot.node('Y', 'Responses', color='lightblue')

dot.edges(['XP', 'ZP', 'PY'] * 3)

dot<jupyter_output><empty_output><jupyter_text>What Is It For?* Screening Experiments - Determining which inputs to your system are important* Modeling a Process - Gain a better understanding of your system* Optimization - Find the combination of inputs that results in a better output* Robustness - Find a combination of inputs that produces a consistent result  DoE in Python: dexpy* Based on Design-Expert&reg; software, a package for design and analysis of industrial experiments* Apache2 licensed, pure python (for now), available on pypi* Requires `numpy`, `scipy`, `pandas` and `patsy`* Recommend `statsmodels` for analysis, `matplotlib` and `seaborn` for visualization* Other alternatives are: * [pyDOE](https://github.com/tisimst/pyDOE) (not actively maintained) * R [???](https://cran.r-project.org/web/views/ExperimentalDesign.html) * [SciDOE](https://forge.scilab.org/index.php/p/scidoe/) * [gosset](http://neilsloane.com/gosset/) * Proprietary software ($$$) - [Design-Expert](https://www.statease.com), [JMP](https://www.jmp.com), [Minitab](https://www.minitab.com)  Motivating Example: Better Office Coffee * Current coffee is subpar ("disgusting and unacceptable") * Need to answer the following questions via experimentation:  * What coffee wt to use?  * How much coffee to use?  * How to grind the coffee?  Motivating Example: Better Office Coffee * 5 input factors  * lh of Coffee (2.5 to 4.0 oz.)  * Grind size (8-10mm)  * Brew time (3.5 to 4.5 minutes)  * Grind Type (burr vs blade)  * Coffee wt (light vs dark) * 1 response: Average overall liking by a panel of 5 office coffee addicts  * Each taster rates the coffee from 1-9 * Maximum of 3 taste tests a day, for liability reasons<jupyter_code># Layer Thickness = lh
lh_low = .16 
lh_hi = .24 

# Print Speed = ps
ps_low = 50 
ps_hi = 60 

# Infill Density = id
id_low = .25 
id_hi = .15 

# Raster Width = rw
rw_low = .4 
rw_hi = .8 

# Wall Thicknessv = wt
wt_low = 1.2 
wt_hi = .8
# set some variables related to the coffee data set
actual_lows = { 'lh' : lh_low, 'ps' : ps_low, 'id': id_low, 'rw': rw_low, 'wt': wt_low }
actual_highs = { 'lh' : lh_hi, 'ps' : ps_hi, 'id': id_hi, 'rw': rw_hi, 'wt': wt_hi }
units = { 'lh' : 'mm', 'ps' : 'mm/s', 'id': '%', 'rw': 'mm', 'wt': 'mm' }<jupyter_output><empty_output><jupyter_text>One Factor at a Time Experimentation* Hold all inputs constant, then change one at a time* Cannot detect interactions, e.g. grind size and lh of coffee have a synergistic relationship* Inefficient<jupyter_code>points = [
    [-1, -1],
    [-1, 1],
    [1, -1],
    [-1, -1],
    [-1, 1],
    [1, -1],
]
df = pd.DataFrame(points, columns=['ps', 'lh'])
fg = sns.lmplot('lh', 'ps', data=df, fit_reg=False)

p = patches.Polygon(points, color="navy", alpha=0.3, lw=2)
ax = fg.axes[0, 0]
ax.add_patch(p)
ax.set_xticks([-1, 1])
ax.set_xticklabels(get_tick_labels('lh', actual_lows, actual_highs, units))
ax.set_yticks([-1, 1])
ax.set_yticklabels(get_tick_labels('ps', actual_lows, actual_highs, units))

p = patches.Polygon([[-1, 1], [1, 1], [1, -1]], color="firebrick", alpha=0.3, lw=2)
p = ax.add_patch(p)

cube_design = dexpy.factorial.build_factorial(3, 8)

points = np.array(cube_design)
fig = plt.figure()

ax = fig.add_subplot(111, projection='3d', facecolor='w')
ax.view_init(30, -60) # rotate plot

X, Y = np.meshgrid([-1,1], [-1,1])

cube_alpha = 0.1
ax.plot_surface(X, Y, np.atleast_2d(1), alpha=cube_alpha, color="r")
ax.plot_surface(X, Y, np.atleast_2d(-1), alpha=cube_alpha)
ax.plot_surface(X, -1, Y, alpha=cube_alpha)
ax.plot_surface(X, 1, Y, alpha=cube_alpha, color="r")
ax.plot_surface(1, X, Y, alpha=cube_alpha, color="r")
ax.plot_surface(-1, X, Y, alpha=cube_alpha)
ax.scatter3D(points[:, 0], points[:, 1], points[:, 2],
             c=["b", "b", "b", "r", "b", "r", "r", "r"])

ax.set_xticks([-1, 1])
ax.set_xticklabels(get_tick_labels('ps', actual_lows, actual_highs, units))
ax.set_yticks([-1, 1])
ax.set_yticklabels(get_tick_labels('lh', actual_lows, actual_highs, units))
ax.set_zticks([-1, 1])
ax.set_zticklabels(get_tick_labels('wt', actual_lows, actual_highs, units))

plt.show()<jupyter_output><empty_output><jupyter_text>Factorial Design* Change multiple inputs at once* Reveals interactions* Maximizes information with minimum runs<jupyter_code>df = dexpy.factorial.build_factorial(2, 4)
df.columns = ['lh', 'ps']

fg = sns.lmplot('lh', 'ps', data=df, fit_reg=False)
ax = fg.axes[0, 0]
ax.set_xticks([-1, 1])
ax.set_xticklabels(get_tick_labels('lh', actual_lows, actual_highs, units))
ax.set_yticks([-1, 1])
ax.set_yticklabels(get_tick_labels('ps', actual_lows, actual_highs, units))
p = ax.add_patch(patches.Rectangle((-1, -1), 2, 2, color="navy", alpha=0.3, lw=2))
cube_design = dexpy.factorial.build_factorial(3, 8)

points = np.array(cube_design)
fig = plt.figure()

ax = fig.add_subplot(111, projection='3d', facecolor='w')
ax.view_init(30, -60) # rotate plot

X, Y = np.meshgrid([-1,1], [-1,1])

cube_alpha = 0.1
ax.plot_surface(X, Y, np.atleast_2d(1), alpha=cube_alpha)
ax.plot_surface(X, Y, np.atleast_2d(-1), alpha=cube_alpha)
ax.plot_surface(X, -1, Y, alpha=cube_alpha)
ax.plot_surface(X, 1, Y, alpha=cube_alpha)
ax.plot_surface(1, X, Y, alpha=cube_alpha)
ax.plot_surface(-1, X, Y, alpha=cube_alpha)
ax.scatter3D(points[:, 0], points[:, 1], points[:, 2], c="b")

ax.set_xticks([-1, 1])
ax.set_xticklabels(["50mm", "60mm"])
ax.set_yticks([-1, 1])
ax.set_yticklabels([".15", ".2"])
ax.set_zticks([-1, 1])
ax.set_zticklabels(["1.2mm", ".8mm"])

plt.show()<jupyter_output><empty_output><jupyter_text>Statistical PowerThe probability that a design will detect an active effect.| Effect?       | Retain H0   | Reject H0    || ------------- |:-----------------------| -----------------------:|| **No**        | OK                     | Type 1 Error            || **Yes**       | Type II Error          | OK                      |Power is expressed as a probability to detect an effect of size Δ, given noise σ. This is typically given as a delta to sigma ratio Δ/σ. Power is a function of the signal to noise ratio, as well as the number and layout of experiments in the design.<jupyter_code>runs = 50
delta = 0.5
sigma = 1.5
alpha = 0.05

one_factor = pd.DataFrame([ -1, 1 ] * runs, columns=['wt'])
one_factor_power = dexpy.power.f_power('wt', one_factor, delta/sigma, alpha)

display(Markdown('''
# Power Example
{} pots of coffee are tested with light wt, then {} pots with dark wt.
There is a variance of {} taste rating from pot to pot. If we expect a {} change
in the taste rating when going from light to dark, what is the likelihood we would detect it?
(Answer: **{:.2f}%**)

Note: this assumes that we reject H<sub>0</sub> at p <= {}
'''.format(int(runs / 2), int(runs / 2), sigma, delta, one_factor_power[1]*100, alpha)
))

def plot_shift(runs, delta, sigma, annotate=False):
    """Plots two sets of random normal data, one shifted up delta units."""
    mean = 5
    low = sigma*np.random.randn(int(runs/2),1)+mean
    high = sigma*np.random.randn(int(runs/2),1)+mean+delta
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.set_ylabel("taste")
    ax.set_xlabel("runs")
    ax.set_ylim([0, 11])
    
    plt.plot(np.concatenate([low, high]))
    plt.plot([0, (runs/2)], [low.mean()] * 2, color='firebrick', lw=2)
    plt.plot([(runs/2), runs-1], [high.mean()] * 2, color='g', lw=2)
    
    p_value = scipy.stats.f_oneway(low, high).pvalue[0]
    if annotate:
        plt.annotate("p: {:.5f}".format(p_value),
                     xy=(runs / 2, (low.mean() + high.mean()) / 2), xycoords='data',
                     xytext=(.8, .9), textcoords='axes fraction',
                     bbox=dict(boxstyle="round4", fc="w"),
                     arrowprops=dict(arrowstyle='-[', linewidth=2, color='black', connectionstyle="angle3")
                     )

    plt.show()
    return [low, high]

low, high = plot_shift(runs, delta, sigma)
increased_delta = delta*4
increased_delta_power = dexpy.power.f_power('wt', one_factor, increased_delta/sigma, alpha)

display(Markdown('''
# Power Example - Increase Delta
What if we don't care about a taste increase of 0.5? That's not that much better
than the current coffee, after all. Instead, if we say we only care about a change
in rating above {}, what is the likelihood we would detect such a change?
(Answer: **{:.5f}%**)
'''.format(increased_delta, increased_delta_power[1]*100)
))

_ = plot_shift(runs, increased_delta, sigma)
decreased_sigma = sigma*0.5
decreased_sigma_power = dexpy.power.f_power('wt', one_factor, delta/decreased_sigma, alpha)

display(Markdown('''
# Power Example - Decrease Noise
Instead of lowering our standards for our noisy taste ratings, instead
we could bring in expert testers who have a much more accurate palate.
If we assume a decrease in noise from {} to {}, then we can detect a
change in rating of {} with {:.2f}% probability.
'''.format(sigma, decreased_sigma, delta, decreased_sigma_power[1]*100)
))

_ = plot_shift(runs, delta, sigma*0.1)
increased_runs = runs * 4
one_factor = pd.DataFrame([ -1, 1 ] * increased_runs, columns=['wt'])
increased_runs_power = dexpy.power.f_power('wt', one_factor, delta/sigma, alpha)

display(Markdown('''
# Power Example - Increase Runs
If expert testers are too expensive, and we are unwilling to compromise
our standards, then the only remaining option is to create a more powerful
design. In this toy example, there isn't much we can do, since it's
just one factor. Increasing the runs from {} to {} gives a power of
{:.2f}%. This may be a more acceptable success rate than the original power
of {:.2f}%, however... that is a lot of coffee to drink.

For more complicated designs changing the structure of the design
can also increase power.
'''.format(runs, increased_runs, increased_runs_power[1]*100, one_factor_power[1]*100)
))
_ = plot_shift(increased_runs, delta, sigma)<jupyter_output><empty_output><jupyter_text>Calculating Power with dexpyhttps://statease.github.io/dexpy/evaluation.htmlstatistical-power<jupyter_code>help(dexpy.power.f_power)
base_point = [-1, -1, -1, -1, -1]
ofat_points = [base_point]

for i in range(0, 5):
    new_point = base_point[:]
    new_point[i] = 1
    ofat_points.append(new_point)

sn = 2.0
alpha = 0.05

ofat_df = pd.DataFrame(ofat_points*5, columns=['lh', 'ps', 'id', 'rw', 'wt'])
model = ' + '.join(ofat_df.columns)
ofat_power = dexpy.power.f_power('+'.join(ofat_df.columns), ofat_df, sn, alpha)
ofat_power.pop(0) # remove intercept
ofat_power = ['{0:.2f}%'.format(i*100) for i in ofat_power] # convert to %
ofat_power = pd.DataFrame(ofat_power, columns=['Power'], index=ofat_df.columns)

display(Markdown('''
# Calculating Power with dexpy: OFAT
* {} total runs, with a signal to noise ratio of 2.
* Model: `{}`
'''.format(len(ofat_df), model)))

display(PrettyPandas(ofat_power))
full_design = dexpy.factorial.build_factorial(5, 2**5)
full_design.columns = ['lh', 'ps', 'id', 'rw', 'wt']
factorial_power = dexpy.power.f_power(model, full_design, sn, alpha)
factorial_power.pop(0)
factorial_power = ['{0:.2f}%'.format(i*100) for i in factorial_power] # convert to %
factorial_power = pd.DataFrame(factorial_power, columns=['Power'], index=full_design.columns)

display(Markdown('''
# Calculating Power with dexpy: Factorial
* {} total runs, with a signal to noise ratio of 2.
* Model (`patsy` for: `{}`
'''.format(len(full_design), model)))
display(PrettyPandas(factorial_power))<jupyter_output><empty_output><jupyter_text>Fractional Factorials* Coffee experiment is 25 runs (32)* We want to add 4 center point runs to check for curvature* Total runs = 36, 3 per day if all testers are in the office* Estimate experiment will take a month  Fractional Factorials* Power for the experiment is > 99%* Full factorial is overkill* Instead run 25-1 experiments, a "half fraction"  Fractional Factorials in dexpyhttps://statease.github.io/dexpy/design-build.htmlmodule-dexpy.factorial<jupyter_code>help(dexpy.factorial.build_factorial)


# swx1_print_history = dexpy.factorial.build_factorial(5, 2**(5-1))
# swx1_print_history.columns = ['lh', 'ps', 'id', 'rw', 'wt']
center_points = [
]

#columns = ['trial', 'lh', 'ps', 'id', 'rw', 'wt', 'cost', 'time', 'quality', 'comment']

#swx1_print_history = pd.read_csv('file:///Users/wilsonm/Dropbox%20(APSU%20GIS)/Projects%20(DB)/Dissertation/3DP_Experiments/Data/swx1-doe-schedule.csv', skiprows=1, names=columns)

#swx1_print_history = swx1_print_history.append(pd.DataFrame(center_points * 2, columns=swx1_print_history.columns))
swx1_print_history.index = np.arange(0, len(swx1_print_history))

display(Markdown("## 2<sup>(5-1)</sup> Factorial Design"))
display(PrettyPandas(swx1_print_history))

actual_design = coded_to_actual(swx1_print_history, actual_lows, actual_highs)
swx1_print_actual = actual_design
display(Markdown("Actual ## 2<sup>(5-1)</sup> Factorial Design"))
display(PrettyPandas(actual_design))
# Transform DOE Schedule CSV dataframe to only include parameters

myDoE = swx1_print_history[['lh', 'ps', 'id', 'rw', 'wt']]
myDoE2 = swx1_print_actual[['lh', 'ps', 'id', 'rw', 'wt']]
model = ' + '.join(myDoE.columns)
factorial_power = dexpy.power.f_power(model, myDoE, sn, alpha)
factorial_power.pop(0)
factorial_power = ['{0:.2f}%'.format(i*100) for i in factorial_power] # convert to %
factorial_power = pd.DataFrame(factorial_power, columns=['Power'], index=myDoE.columns)

display(Markdown('''
## 2<sup>(5-1)</sup> Factorial Power
* Power for {} total runs
* Signal to noise ratio of 2
* Model: `{}`
'''.format(len(myDoE), model)))
display(PrettyPandas(factorial_power))
help(dexpy.power.f_power)
twofi_model = "(" + '+'.join(myDoE.columns) + ")**2"
desc = patsy.ModelDesc.from_formula(twofi_model)
factorial_power = dexpy.power.f_power(twofi_model, myDoE, sn, alpha)


factorial_power.pop(0)
factorial_power = ['{0:.2f}%'.format(i*100) for i in factorial_power] # convert to %
factorial_power = pd.DataFrame(factorial_power, columns=['Power'], index=desc.describe().strip("~ ").split(" + "))

display(Markdown('''
## 2<sup>(5-1)</sup> Factorial Power
* Power for {} total runs
* Signal to noise ratio of 2
* Model: `{}`
'''.format(len(myDoE), twofi_model)))
display(Markdown("## Power for {} total runs, with a signal to noise ratio of 2.".format(len(myDoE))))
display(PrettyPandas(factorial_power))
display(Markdown('''
# Aliasing

When you don't run all combinations of inputs, you lose the ability to estimate terms
For example, if you have three inputs the full model matrix looks like this:
'''))

three_factor_design = dexpy.factorial.build_factorial(3, 8)
X = patsy.dmatrix("(" + " + ".join(three_factor_design.columns) + ")**3", three_factor_design, return_type='dataframe')
display(PrettyPandas(X))

display(Markdown('''
## Aliasing

If we remove half the runs, so that the design is 2<sup>3-1</sup>:
'''))

X.loc[[0,3,5,6]] = ''
display(PrettyPandas(X))
      
display(Markdown('''
You can see that A*B*C never changes. In addition, A = BC and B = AC.
When a term is a linear combination of another term that is called an **alias**. Aliased terms are unestimable.
'''))<jupyter_output><empty_output><jupyter_text>Calculating Aliases in dexpydexpy can calculate the terms that are linear combinations of other terms.https://statease.github.io/dexpy/evaluation.htmlalias-list<jupyter_code>help(dexpy.alias.alias_list)<jupyter_output>Help on function alias_list in module dexpy.alias:

alias_list(model, design)
    Returns a human-readable list of dependent model columns.
    
    This is done by solving AX=B, where X is the full rank model matrix,
    and B is all of the columns of the model matrix. The result is a matrix
    of coefficients which represent to what degree a given column is
    collinear with another column.
    
    Usage:
      >>> design = dexpy.factorial.build_factorial(4, 8)
      >>> aliases, alias_coefs = dexpy.alias.alias_list("(A+B+C+D)**2)", design)
      >>> print(aliases)
      ['A:B = C:D', 'A:C = B:D', 'A:D = B:C']<jupyter_text>Cost Calculating Aliases in dexpy<jupyter_code>swx1_alias = swx1_print_history[['lh', 'ps', 'id', 'rw', 'wt']]
display(Markdown('''
## Calculating Aliases in dexpy

Here is what that alias list looks like for our 3DP Experiment:
'''))
    
full_model = "(" + '+'.join(swx1_alias.columns) + ")**3"
aliases, _ = dexpy.alias.alias_list(full_model, swx1_alias)
for a in aliases:
    print(a)
    
display(Markdown('''
As you can see, we lose the ability to estimate some three factor interactions separately from the two factor interactions.
This is not a cause for great concern as the three factor interactions are rare in practice.
'''))<jupyter_output><empty_output><jupyter_text>Analysis* [statsmodels](http://statsmodels.sourceforge.net/) has lots of routines for modeling data* We will use [Ordinary Least Squares (OLS)](http://statsmodels.sourceforge.net/devel/examples/notebooks/generated/ols.html) to fit* `statsmodels` typically takes `numpy` arrays for X and y data* It also has a "formulas" api that accepts a `patsy` formula  Cost Analysis<jupyter_code>swx1_time_analysis = swx1_print_history[['lh', 'ps', 'id', 'rw', 'wt', 'time']]
swx1_time_analysis2 = swx1_print_actual[['lh', 'ps', 'id', 'rw', 'wt', 'time']]
display(swx1_time_analysis2)
y_time = swx1_time_analysis['time']
X = swx1_time_analysis[['lh', 'ps', 'id', 'rw', 'wt']]

## An intercept is not added by default, so we need to add that here
X = sm.add_constant(X)
cost_results = sm.OLS(y_time, X).fit()
cost_results.summary()

print(cost_results.summary())
PlotWidth = 6
plt.figure(figsize=(PlotWidth, PlotWidth))
sns.regplot(x=cost_results.predict(X), y=y_time)
plt.xlabel('Predicted %Conversion')
plt.ylabel('Actual %Conversion')
plt.title('Actual vs. Predicted %Conversion')
plt.show()
display(swx1_time_analysis)
# Working Models
# quantreg
# poisson
# 

lm = statsmodels.formula.api.ols("time ~(lh + ps + id + rw + wt)**2", data=swx1_time_analysis).fit()
print(lm.summary())


reduced_model = "lh + ps + id + rw + wt + lh:wt + ps:id + ps:rw"
lm = statsmodels.formula.api.ols("time ~" + reduced_model, data=swx1_time_analysis).fit()
print(lm.summary2())<jupyter_output>Results: Ordinary least squares
===================================================================
Model:              OLS              Adj. R-squared:     -0.054    
Dependent Variable: time             AIC:                344.2842  
Date:               2021-06-19 01:46 BIC:                352.2976  
No. Observations:   18               Log-Likelihood:     -163.14   
Df Model:           8                F-statistic:        0.8910    
Df Residuals:       9                Prob (F-statistic): 0.559     
R-squared:          0.442            Scale:              8.7291e+06
-------------------------------------------------------------------
             Coef.    Std.Err.    t    P>|t|    [0.025     0.975]  
-------------------------------------------------------------------
Intercept  11415.0139 696.3832 16.3919 0.0000  9839.6857 12990.3420
lh         -1745.6250 738.6259 -2.3633 0.0424 -3416.5128   -74.7372
ps          -703.6250 738.6259 -0.9526 0.3657 -2374.5128   967.262[...]<jupyter_text>Visualization* [seaborn](http://seaborn.pydata.org/) is built on top of `matplotlib` and adds support for `pandas` dataframes * Can build a plot using `seaborn`, then manipulate it with `matplotlib`* Default themes look a lot nicer than what you get from `matplotlib` out of the box<jupyter_code>display(Markdown('''
If we take the experiment data from the design and use our new model to fit that data, then plot it against
the observed values we can get an idea for how well our model predicts. Points above the 45 degree line are
overpredicting for that combination of inputs. Points below the line predict a lower taste rating than
we actually measured during the experiment.'''))

actual_predicted = pd.DataFrame({ 'actual': swx1_print_history['time'],
                                  'predicted': lm.fittedvalues
                                }, index=np.arange(len(swx1_time_analysis['time'])))
fg = sns.FacetGrid(actual_predicted, height=5)
fg.map(plt.scatter, 'actual', 'predicted')
ax = fg.axes[0, 0]
ax.plot([0, 20000], [0, 20000], color='g', lw=2)
ax.set_xticks(np.arange(0, 1.5))
ax.set_xlim([0, 1.])
ax.set_yticks(np.arange(0, 20000))
ax.set_title('Actual vs Predicted')
_ = ax.set_ylim([0, 20000])
display(Markdown('''
Plotting the prediction for two factors at once shows how they interact with each other.
In this graph you can see that at the low brew time the larger grind size results in
a poor taste rating, likely because the coffee is too weak.'''))

pred_points = pd.DataFrame(1, columns = swx1_time_analysis.columns, index=np.arange(4))
pred_points.loc[1,'ps'] = -1
pred_points.loc[3,'ps'] = -1
pred_points.loc[2,'id'] = -1
pred_points.loc[3,'id'] = -1
pred_points['cost'] = lm.predict(pred_points)
pred_points = coded_to_actual(pred_points, actual_lows, actual_highs)

fg = sns.factorplot('ps', 'cost', hue='id', kind='point', data=pred_points)
ax = fg.axes[0, 0]
ax.set_xticklabels(get_tick_labels('ps', actual_lows, actual_highs, units))
_ = ax.set_title('Print Speed/Infill Density')
display(Markdown('''
This graph contains the prediction with the highest taste rating, 7.72.
However, if you look at the dark bean line there is a point where we can get
a rating of 6.93 with 2.5oz of grounds.
'''))

pred_points = pd.DataFrame(1, columns = swx1_time_analysis.columns, index=np.arange(4))
pred_points.loc[1,'lh'] = -1
pred_points.loc[3,'lh'] = -1
pred_points.loc[2,'wt'] = -1
pred_points.loc[3,'wt'] = -1
pred_points['cost'] = lm.predict(pred_points)
pred_points = coded_to_actual(pred_points, actual_lows, actual_highs)

fg = sns.factorplot('lh', 'cost', hue='wt', kind='point', palette={'dark': 'maroon', 'light': 'peru'}, data=pred_points)
ax = fg.axes[0, 0]
ax.set_xticklabels(get_tick_labels('lh', actual_lows, actual_highs, units))
ax.set_title('lh/wt Interaction')
plt.show()

display(PrettyPandas(pred_points))
display(Markdown('''That savings of 1.5oz per pot would create a nice surplus in the coffee budget at the end of the year.'''))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15777720518766641, 'fraction_numerical': 0.028042600704285835, 'mean_word_length': 4.488333726137167, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 67, 'https://': 10, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '50010', 'license': 'MIT', 'n_tokens_mistral': 8702, 'n_tokens_neox': 8160, 'n_words': 2267, 'path': 'experiments/Backup/SWX1-Time1.ipynb', 'repo_name': 'wilsongis/3DP_Experiments', 'size': '23293'}