Skip to content

Instantly share code, notes, and snippets.

@ryanjdillon
Created March 13, 2017 14:51
Show Gist options
  • Save ryanjdillon/83f035853319c433c3710c9cf1d9705f to your computer and use it in GitHub Desktop.
Save ryanjdillon/83f035853319c433c3710c9cf1d9705f to your computer and use it in GitHub Desktop.
Algae light response - Calculate the level of oxygen at a give light level
import seaborn
def light_response(light, EC50, hillslop=2.5):
'''Typical response function for EC50 tests'''
import numpy
# http://www.sigmaplot.co.uk/splot/products/sigmaplot/productuses/prod-uses43.php
# Response to concentration
dose_min = light.min()
dose_max = light.max()
o2 = dose_min + ((dose_max - dose_min)/(1+(light/EC50)**-hillslope))
return o2
def nearest(a, value):
'''Find the element in array `a` nearest in value to `value`'''
a_near = min(a, key=lambda x: abs(x - value))
idx_near = numpy.where(a == a_near)[0][0]
return a_near, idx_near
def solve_for_light(light_fit, o2_fit, o2_val):
'''Get nearest light value from given o2 value'''
o2_ans, o2_idx = nearest(o2_fit, o2_val)
light_ans = light_fit[o2_idx]
print('Light at {} O2 level: {}'.format(o2_val, light_ans))
return light_ans, o2_ans
def plot_response(light_data, o2_data, light_fit, o2_fit, light_ans, o2_ans):
'''Plot the data, the fit, and the solution'''
import matplotlib.pyplot as plt
plt.plot(light_data, o2_data, label='data')
plt.plot(light_fit, o2_fit, label='fit')
plt.scatter(light_ans, o2_ans, label='{} o2'.format(o2_ans))
plt.legend()
plt.show()
return None
def quadratic_method(light_data, o2_data, o2_val, n_fit):
'''Use a quadratic fit of the data to estimate light at `o2_val`'''
# Make a 3rd order fit of the curve
fit_coeffs = numpy.polyfit(light_data, o2_data, deg=2)
# Make a range of light values to calculate an o2 value from the fit
light_fit = range(n_fit)
# Cacluate the o2 values from the fit
o2_fit = numpy.polyval(fit_coeffs, light_fit)
light_ans, o2_ans = solve_for_light(light_fit, o2_fit, o2_val)
plot_response(light_data, o2_data, light_fit, o2_fit, light_ans, o2_ans)
return None
def scipy_optimize_method(light_data, o2_data, o2_val, n_fit):
'''Use curf_fit() to fit the data to estimate light at `o2_val`.
Note
----
This assumes that you know the function whose parameters you want
to fit
'''
import scipy.optimize
popt, pcov = scipy.optimize.curve_fit(light_response, light_data, o2_data)
o2_fit = light_response(light_data, *popt)
# Make a range of light values to calculate an o2 value from the fit
light_fit = range(n_fit)
light_ans, o2_ans = solve_for_light(light_fit, o2_fit, o2_val)
plot_response(light_data, o2_data, light_fit, o2_fit, light_ans, o2_ans)
return None
if __name__ == '__main__':
import numpy
# Params for fake data
n = 1000
EC50 = 0.2*n
hillslope = 2.5
# Generate some fake data
light_data = numpy.arange(n)
o2 = light_response(light_data, EC50, hillslop=2.5)
o2_data = o2 + numpy.random.normal(size=len(light_data))*(0.1*light_data.max())
o2_data[o2_data < 0] = 0
# Value of O2 to solve for light with
o2_val = 400
# Number of points to generate from fits, higher gives better accuracy
n_fit = 1000
# Try both methods with plots
quadratic_method(light_data, o2_data, o2_val, n_fit)
scipy_optimize_method(light_data, o2_data, o2_val, n_fit)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment