We want to have, at the end of the day, a function that takes the rates in some fashion and returns the likelihood for the CH82 model. This is exactly what the code below does.
In [1]:
from dcpyps import samples, dataset, dcio, scalcslib as scl
from numpy import log, average, abs, array, any, all
mechanism = samples.CH82()
tres = 1e-4
tcrit = 4e-3
filename = "../data/CH82.scn"
ioffset, nint, calfac, header = dcio.scn_read_header(filename)
tint, iampl, iprops = dcio.scn_read_data(filename, ioffset, nint, calfac)
rec1 = dataset.SCRecord(filename, header, tint, iampl, iprops)
# Impose resolution, get open/shut times and bursts.
rec1.impose_resolution(tres)
rec1.get_open_shut_periods()
rec1.get_bursts(tcrit)
blength = rec1.get_burst_length_list()
openings = rec1.get_openings_burst_list()
# not sure impose_resolution works...
bursts = [u for u in rec1.bursts if all(abs(array(u)*1e-3) > tres)]
# Prepare parameter dict for simplex
opts = {}
opts['mec'] = mechanism
opts['conc'] = 50e-9
opts['tres'] = tres
opts['tcrit'] = tcrit
opts['isCHS'] = True
opts['data'] = bursts
# MAXIMUM LIKELIHOOD FIT.
def dcpyps(x):
start_lik, th = scl.HJClik(log(x), opts)
return start_lik
theta = mechanism.theta()
We want to create a function which does the exact same thing as the one above, but which uses internally the new dcprogs. In order to do that, we have dc-pyps translate its own input into a Q-matrix.
In [12]:
from numpy import sum
from dcprogs import read_idealized_bursts
from dcprogs.likelihood import QMatrix, Log10Likelihood, MissedEventsG
from dcprogs.likelihood.optimization import reduce_likelihood
name = 'CH82'
tau = 1e-4
tcrit = 4e-3
nopen = 2
bursts = read_idealized_bursts(name, tau=tau, tcrit=tcrit)
print str(bursts[2])
likelihood = Log10Likelihood(bursts, nopen, tau, tcrit)
def dcprogs(x):
mechanism.theta_unsqueeze(x)
mechanism.set_eff('c', opts['conc'])
return -likelihood(mechanism.Q) * log(10)
The code below compares the result from both codes for somewhat random inputs.
In [9]:
from numpy import random
x = theta + 10.0*random.uniform(low=-1, high=1, size=8)
print(repr(x))
print("DCPYPS {0}".format(dcpyps(x)))
print("DCPROGS {0}".format(dcprogs(x)))