Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Autoprotocol improvements #88

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

[WIP] Autoprotocol improvements #88

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0c70384
Backport improvements from pymcmodels to assaytools version. Lognorma…
jchodera Apr 21, 2017
ebc2a3b
Updates to autoprotocol analysis
jchodera Apr 22, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
85 changes: 56 additions & 29 deletions AssayTools/analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,30 +44,53 @@
# SUBROUTINES
#=============================================================================================

def LogNormalWrapper(name, mean, stddev):
def LogNormalWrapper(name, mean, stddev, observed=False, value=None, log_trace=True, trace=True):
"""
Create a PyMC Normal stochastic, automatically converting parameters to be appropriate for a `LogNormal` distribution.
Note that the resulting distribution is Normal, not LogNormal.
This is appropriate if we want to describe the log of a volume or a fluorescence reading.

# TODO: Automatically generate corresponding deterministic
# TODO: Squeeze out mean=0 entries from stochastic to allow for zero concentrations

Parameters
----------
mean : float
Mean of exp(X), where X is lognormal variate.
stddev : float
Standard deviation of exp(X), where X is lognormal variate.
observed : bool, optional, False
If True, the values are observed
value : float, optional, default=None
The observed values; observed must be True if provided.

Returns
-------
stochastic : pymc.Normal
Normal stochastic for random variable X
Stochastic version of log quantity
deterministic
Deterministic version of real quantity

"""
# Compute parameters of lognormal distribution
mu = np.log(mean**2 / np.sqrt(stddev**2 + mean**2))
tau = np.sqrt(np.log(1.0 + (stddev/mean)**2))**(-2)
stochastic = pymc.Normal(name, mu=mu, tau=tau)
return stochastic
@pymc.deterministic
def mu(mean=mean, stddev=stddev):
return np.log(mean**2 / np.sqrt(stddev**2 + mean**2))
@pymc.deterministic
def tau(mean=mean, stddev=stddev):
return np.sqrt(np.log(1.0 + (stddev/mean)**2))**(-2)

# Create stochastic
logname = 'log ' + name
stochastic = pymc.Normal(logname, mu=mu, tau=tau, trace=log_trace, observed=observed, value=value)

# Create deterministic
if observed:
deterministic = value
else:
deterministic = pymc.Lambda(name, lambda log_quantity=stochastic: np.exp(log_quantity), trace=trace)

return stochastic, deterministic

def wellname(well):
"""
Expand Down Expand Up @@ -156,8 +179,9 @@ def _create_solutions_model(self):
for solution in self.solutions.values():
if solution.species is None:
continue # skip buffers or pure solvents
name = 'log concentration of %s' % solution.shortname
self.model[name] = LogNormalWrapper(name, mean=solution.concentration.to_base_units().m, stddev=solution.uncertainty.to_base_units().m)
name = 'concentration of %s' % solution.shortname
self.model['log ' + name], self.model[name] = LogNormalWrapper(name, mean=solution.concentration.to_base_units().m, stddev=solution.uncertainty.to_base_units().m)
self.parameter_names['solution concentrations'].append('log ' + name)
self.parameter_names['solution concentrations'].append(name)

def _identify_ligand_names(self):
Expand Down Expand Up @@ -192,26 +216,19 @@ def _create_dispensing_model(self):
log_volumes = list() # log volumes are in Liters
for component in well.properties['contents']:
name = 'volume of %s dispensed into well %s' % (component, wellname(well))
logname = 'log %s' % name
logname = 'log ' + name
(volume, error) = well.properties['contents'][component]
log_volume_dispensed = LogNormalWrapper(logname, mean=volume.to_base_units().m, stddev=error.to_base_units().m)
self.model[logname] = log_volume_dispensed
#self.parameter_names['dispensed volumes'].append(logname)
log_volumes.append(log_volume_dispensed)
# Store real (non-log) value
self.model[name] = pymc.Lambda(name, lambda log_quantity=self.model[logname]: np.exp(log_quantity) / volume_unit.to_base_units().m, trace=True)
self.model[logname], self.model[name] = LogNormalWrapper(name, mean=volume.to_base_units().m, stddev=error.to_base_units().m)
log_volumes.append(self.model[logname])
self.parameter_names['dispensed volumes'].append(logname)
self.parameter_names['dispensed volumes'].append(name)

# Total volume in well
name = 'volume of well %s' % wellname(well)
logname = 'log %s' % name
@pymc.deterministic(name=logname)
def log_total_volume(log_volumes=log_volumes):
return logsumexp(log_volumes)
self.model[logname] = log_total_volume
#self.parameter_names['well volumes'].append(logname)
# Store real (non-log) value
self.model[name] = pymc.Lambda(name, lambda log_quantity=self.model[logname]: np.exp(log_quantity) / volume_unit.to_base_units().m, trace=True)
logname = 'log ' + name
self.model[logname] = pymc.Lambda(logname, lambda log_volumes=log_volumes: logsumexp(log_volumes))
self.model[name] = pymc.Lambda(name, lambda log_total_volume=self.model[logname]: np.exp(log_total_volume), trace=True)
self.parameter_names['well volumes'].append(logname)
self.parameter_names['well volumes'].append(name)

# Total concentration of all species involving each component in well
Expand All @@ -231,7 +248,7 @@ def log_total_volume(log_volumes=log_volumes):
def log_total_concentration(log_solution_concentration=log_solution_concentration, log_solution_volume=log_solution_volume, log_total_volume=log_total_volume):
return log_solution_concentration + log_solution_volume - log_total_volume
self.model[logname] = log_total_concentration
#self.parameter_names['well concentrations'].append(logname)
self.parameter_names['well concentrations'].append(logname)
# Store real (non-log) value
self.model[name] = pymc.Lambda(name, lambda log_quantity=self.model[logname]: np.exp(log_quantity) / concentration_unit.to_base_units().m, trace=True)
self.parameter_names['well concentrations'].append(name)
Expand Down Expand Up @@ -638,6 +655,10 @@ def _create_fluorescence_model(self):
self.model[logname] = pymc.Uniform(logname, lower=MIN_LOG_FLUORESCENCE_UNCERTAINTY, upper=MAX_LOG_FLUORESCENCE_UNCERTAINTY, value=0)
#self.parameter_names['fluorescence'].append(name)

name = 'top fluorescence uncertainty'
self.model[name] = pymc.Lambda(name, lambda log_sigma=self.model[logname] : np.exp(log_sigma))
#self.parameter_names['fluorescence'].append(name)

name = 'top fluorescence precision'
self.model[name] = pymc.Lambda(name, lambda log_sigma=self.model[logname] : np.exp(-2*log_sigma))
self.parameter_names['fluorescence'].append(name)
Expand All @@ -647,6 +668,10 @@ def _create_fluorescence_model(self):
self.model[name] = pymc.Uniform(name, lower=MIN_LOG_FLUORESCENCE_INTENSITY, upper=MAX_LOG_FLUORESCENCE_INTENSITY, value=0)
self.parameter_names['fluorescence'].append(name)

name = 'bottom fluorescence uncertainty'
self.model[name] = pymc.Lambda(name, lambda log_sigma=self.model[logname] : np.exp(log_sigma))
#self.parameter_names['fluorescence'].append(name)

logname = 'bottom fluorescence log uncertainty'
self.model[logname] = pymc.Uniform(logname, lower=MIN_LOG_FLUORESCENCE_UNCERTAINTY, upper=MAX_LOG_FLUORESCENCE_UNCERTAINTY, value=0)
#self.parameter_names['fluorescence'].append(name)
Expand Down Expand Up @@ -728,7 +753,9 @@ def fluorescence_model(log_well_volume=log_well_volume,

measured_fluorescence = measurements['fluorescence'][(excitation_wavelength, emission_wavelength, geometry)]
name = 'measured %s fluorescence of well %s at excitation wavelength %s and emission wavelength %s' % (geometry, wellname(well), excitation_wavelength, emission_wavelength)
self.model[name] = pymc.Normal(name, mu=fluorescence_model, tau=self.model['%s fluorescence precision' % geometry], observed=True, value=measured_fluorescence)
logname = 'log ' + name
self.model[logname], self.model[name] = LogNormalWrapper(name, mean=fluorescence_model, stddev=self.model['%s fluorescence uncertainty' % geometry], observed=True, value=measured_fluorescence)
self.parameter_names['fluorescence'].append(logname)
self.parameter_names['fluorescence'].append(name)

def map_fit(self):
Expand Down Expand Up @@ -781,15 +808,15 @@ def run_mcmc(self, dbfilename='output'):
# TODO: Allow
mcmc = pymc.MCMC(self.model, db='sqlite', name='output', verbose=True)

nthin = 10
nburn = nthin*100
niter = nthin*100
nthin = 20
nburn = nthin*0
niter = nthin*10000

# Specify initial parameter standard deviations to apply to specific classes of parameters
keywords = {
'concentration' : 0.1,
'affinity' : 0.1,
'volume' : 0.01,
'volume' : 0.1,
}

print('Assigning initial guesses for Metropolis step method proposal standard deviations:')
Expand Down