%%javascript
IPython.OutputArea.auto_scroll_threshold = 9999;

%matplotlib inline

import re
import numpy as np

"""
Explanation: Let's consider some data for a fluorescence assay prepared with a Tecan HP D300 digital dispenser.
A summary of the experimental procedure:


Two protein assay plates were prepared with p38 and Src plus  Bosutinib, using rows A-E and columns 1-12, varying protein concentration from row to row (high to low where E is just Buffer)  and Bosutinib concentration from column to column (high to low where 12 has no Bosutinib). 
Bosutinib concentrations started at 20 uM and went down by ½ log as calcuated by the D300 dispenser. Src concentrations were ~1, 0.5, 0.25, and 0.125 uM.  p38 concentrations were 10, 5, 2.5, 1.25 uM.
End of explanation
"""


# Load data
filename = 'data/12022014-019-src-280scan-reboot-bottom.txt' # data filename
from assaytools import platereader
[SRC_280, SRC_280_x, SRC_280_x_num] = platereader.read_emission_spectra_text(filename)

## Plot the emission spectrum for 280 nm excitation

from assaytools import plots
plots.plot_scans(SRC_280, SRC_280_x, SRC_280_x_num)

"""
Explanation: Load fluorescence emission scans
End of explanation
"""


# Find desired emission wavelength to fit.
wavelength_to_fit = '480'
index = SRC_280_x_num.index(wavelength_to_fit)

# Specify row to process
row_index = 0 # row 0, 1, 2, or 3
Pstated_rows = np.array([1, 0.5, 0.25, 0.125]) * 1e-6 # protein concentrations, M

# Extract fluorescence intensities at desired wavelength.
bottom_complex_fluorescence = np.array([ float(SRC_280[i,index]) for i in range(row_index*12+1,row_index*12+13) ])
bottom_ligand_fluorescence = np.array([ float(SRC_280[i,index]) for i in range(49,61) ])

# Stated concentrations of protein and ligand.
Pstated = Pstated_rows[row_index] * np.ones([12],np.float64) # protein concentration, M
Lstated = 20.0e-6 / np.array([10**(float(i)/2.0) for i in range(12)])

# Assay configuration details
assay_volume = 75e-6 # well volume, L
well_area = 0.1586 # well area, cm^2 # half-area wells were used here
path_length = assay_volume / well_area

# Uncertainties in protein and ligand concentrations.
dPstated = 0.10 * Pstated # protein concentration uncertainty
dLstated = 0.08 * Lstated # ligand concentraiton uncertainty (due to gravimetric preparation and HP D300 dispensing)

"""
Explanation: Find desired emission wavelength to fit.
End of explanation
"""


# Create the pymc model
from assaytools import pymcmodels
pymc_model = pymcmodels.make_model(Pstated, dPstated, Lstated, dLstated,
               bottom_complex_fluorescence=bottom_complex_fluorescence,
               bottom_ligand_fluorescence=bottom_ligand_fluorescence,
               use_primary_inner_filter_correction=True, 
               assay_volume=assay_volume, well_area=well_area)

"""
Explanation: Create pymc model
End of explanation
"""


figure = plots.plot_measurements(Lstated, Pstated, pymc_model)

# Find the maximum a posteriori fit (will only be local optimum, and several cycles are needed for reasonable fit)
map = pymcmodels.map_fit(pymc_model)

figure = plots.plot_measurements(Lstated, Pstated, pymc_model, map=map)

# Sample the model posterior with MCMC
mcmc = pymcmodels.run_mcmc(pymc_model)

plots.plot_mcmc_results(Lstated, Pstated, path_length, mcmc)

figure = plots.plot_measurements(Lstated, Pstated, pymc_model, mcmc=mcmc)

# Fit maximum a posteriori estimate one more time
map = pymcmodels.map_fit(pymc_model)

pymcmodels.show_summary(pymc_model, map, mcmc)

"""
Explanation: Plot fluorescence measurements of complex and ligand as a function of ligand concentration
End of explanation
"""


