︠c54e2089-deaf-40d2-94a1-1577acd999a9︠
# See http://www.packtpub.com/article/what-can-you-do-with-sagemath

# First we want to correlate the concentration of bacteria in the liquid to the optical density
# so we load the data (this could come from a file).
import numpy
CFU = numpy.array([2.60E+08, 3.14E+08, 3.70E+08, 4.62E+08,
    8.56E+08, 1.39E+09, 1.84E+09])
optical_density = numpy.array([0.083, 0.125, 0.213, 0.234,
    0.604, 1.092, 1.141])
OD_vs_CFU = zip(optical_density, CFU)
︠5c8c2b0f-31f4-4047-9828-d7b5b0dd5aa0︠
# Plot the data
data_plot = scatter_plot(OD_vs_CFU, markersize=20,
     facecolor='red', axes_labels=['OD at 600nm', 'CFU/ml'])
show(data_plot)
︠732dcf4e-b27a-4a3a-ba0a-9163272e8283︠
# Looks linear, lets try and fit a straight line

# This is the function we will try and fit to the data
var('OD, slope, intercept')
def standard_curve(OD, slope, intercept):
   return OD * slope + intercept

# Now we call SAGE to do the fit
std_params = find_fit(OD_vs_CFU, standard_curve,
    parameters=[slope, intercept],
    variables=[OD], initial_guess=[1e9, 3e8],
    solution_dict = True)
# And print the retuls
for param, value in std_params.iteritems():
     print(str(param) + ' = %e' % value)
︠79f10f9b-a0a9-46f2-8e60-dc7432f501df︠
# We can plot the result
fit_plot = plot(standard_curve(OD, std_params[slope],
     std_params[intercept]), (OD, 0, 1.2))
# And show it together with the data
show(data_plot + fit_plot)
︠682f80b6-b483-4369-b6f7-bc1bb52da81d︠
# We will use the intercept = 1.237283e+08 slope = 1.324714e+09 as the coorelation between number
# of bacteria and the optical density.  So now it is time to check out the real data, which is measurements
# of optical density over time
sample_times = numpy.array([0, 20, 40, 60, 80, 100, 120,
    140, 160, 180, 200, 220, 240, 280, 360, 380, 400, 420,
    440, 460, 500, 520, 540, 560, 580, 600, 620, 640, 660,
    680, 700, 720, 760, 1240, 1440, 1460, 1500, 1560])

OD_readings = numpy.array([0.083, 0.087, 0.116, 0.119, 0.122,
    0.123, 0.125, 0.131, 0.138, 0.142, 0.158, 0.177, 0.213,
    0.234, 0.424, 0.604, 0.674, 0.726, 0.758, 0.828, 0.919,
    0.996, 1.024, 1.066, 1.092, 1.107, 1.113, 1.116, 1.12,
    1.129, 1.132, 1.135, 1.141, 1.109, 1.004, 0.984, 0.972, 0.952])

# Now convert OD_readings to number of bacteria using the slope and intercept computed above
concentrations = standard_curve(OD_readings, std_params[slope],
     std_params[intercept])

# The final data is a combination of sample times and concentrations
exp_data = zip(sample_times, concentrations)
︠6c4c12aa-2e31-40b5-832a-2f9f10bdd006︠
# We can plot the data
data_plot = scatter_plot(exp_data, markersize=20, facecolor='red',
     axes_labels=['time (sec)', 'CFU/ml'])
show(data_plot)
︠f54ad300-ed04-4b73-bd9d-26ae2fc23fc4︠
# Now we need to guess a type of function to attempt a fit.  Here is where SAGE doesn't help
# us.  This example seems to be fit by http://en.wikipedia.org/wiki/Gompertz_function but this
# is where you would try several different fit functions.

# Here is the function we guess with the parameters we want SAGE to find
var('t, max_rate, lag_time, y_max, y0')
def gompertz(t, max_rate, lag_time, y_max, y0):
    """Define a growth model based upon the Gompertz growth curve"""
    return y0 + (y_max - y0) * numpy.exp(-numpy.exp(1.0 +
    max_rate * numpy.exp(1) * (lag_time - t) / (y_max - y0)))
# Estimated parameter values for initial guess
max_rate_est = (1.4e9 - 5e8)/200.0
lag_time_est = 100
y_max_est = 1.7e9
y0_est = 2e8
︠53603e3d-fd7f-43c7-b2e0-4d8d754e5f67︠
# Now do the fit
gompertz_params = find_fit(exp_data, gompertz,
    parameters=[max_rate, lag_time, y_max, y0],
    variables=[t],
    initial_guess=[max_rate_est, lag_time_est, y_max_est, y0_est],
    solution_dict = True)

for param,value in gompertz_params.iteritems():
    print(str(param) + ' = %e' % value)
︠e6f9fa2a-7276-45e0-a565-735ff00d79b0︠
# And the plot
gompertz_model_plot = plot(gompertz(t, gompertz_params[max_rate],
   gompertz_params[lag_time], gompertz_params[y_max],
   gompertz_params[y0]), (t, 0, sample_times.max()))

show(gompertz_model_plot + data_plot)
︠eb00c7bb-d909-4340-8292-e2fa36add1cb︠