Algae light response - Calculate the level of oxygen at a give light level

algae_light_response.py

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)