Unknown error/warning with Quickmodel
Closed this issue · 9 comments
Whenever I run quickmodel.py I see the following warning for only the analysis of 1st row of the plate. I don't understand what is causing it, but it is a curiosity that it doesn't happen in the analysis of the following rows:
/Users/isikm/opt/anaconda/lib/python2.7/site-packages/assaytools-0.2.0-py2.7.egg/assaytools/pymcmodels.py:77: RuntimeWarning: divide by zero encountered in power tau = np.sqrt(np.log(1.0 + (stddev/mean)**2))**(-2)
That's due to #83.
This shouldn't be happening unless a quantity that shouldn't be zero is being set to zero. It is supposed to ensure that concentrations that are set to zero are handled differently, but is there some other input quantity that is accidentally being set to zero?
Question: Could any of your observed fluorescence values be zero, even though the fluorescence should never be zero?
Actually, I was wrong---it's because stddev is small compared to mean for some observed measurement. I've added some more error checking in #93 to test:
Exception: 'np.sqrt(np.log(1.0 + (stddev/mean)**2))' for variable 'top_complex_fluorescence' is zero: mean = [ 2.16041187e+03 1.81939230e+07 2.31772231e+03 2.18967032e+14
2.25150425e+05 2.41688443e+03 2.41961246e+03 1.88126534e+05
2.46657280e+03 3.06447614e+07 1.90888875e+11 5.18205257e+03
2.57256151e+03 9.70449636e+04 2.41046558e+03 5.24266051e+06
2.87262963e+05 1.56344936e+07 2.53205567e+03 1.92708865e+14
2.66307672e+03 1.01737636e+09 1.49829810e+07 2.50573925e+03], stddev = 72.5584153594
Any idea why there would be a value that is 1e+14?
I updated my assaytools to the last commit in Improvements branch (commit).
This is the error message I am getting when quickmodel fails:
File "/Users/isikm/opt/anaconda/lib/python2.7/site-packages/assaytools-0.2.0-py2.7.egg/assaytools/pymcmodels.py", line 81, in tau
raise Exception("'np.sqrt(np.log(1.0 + (stddev/mean)**2))' for variable '%s' is zero: mean = %s, stddev = %s" % (name, mean, stddev))
Exception: 'np.sqrt(np.log(1.0 + (stddev/mean)**2))' for variable 'top_complex_fluorescence' is zero: mean = [ 1.65799135e+07 5.93720599e+14 7.33976650e+02 2.13168885e+04
1.34352179e+02 7.16161337e+12 6.39596253e+01 7.06290971e+01
1.57123273e+06 1.48503707e+02 8.35613405e+12 6.18131246e+01], stddev = 9.28626378933
@sonyahanson : ANy idea why there are numbers like 1e+12 in the fluorescence measurements?
@jchodera
These are my input files. This experiment has a different plate layout. First 6 rows have protein and last 2 rows have buffer. So you must use quickmodel_layout2.py.
Something weird is going on here because (1) the written values are supposed to be observed values (e.g. fluorescence measurements), and (2) they are different each time I run the analysis script on the same data!
try 1:
File "/Users/choderaj/miniconda3/lib/python3.5/site-packages/assaytools-0.2.0-py3.5.egg/assaytools/pymcmodels.py", line 81, in tau
raise Exception("'np.sqrt(np.log(1.0 + (stddev/mean)**2))' for variable '%s' is zero: mean = %s, stddev = %s" % (name, mean, stddev))
Exception: 'np.sqrt(np.log(1.0 + (stddev/mean)**2))' for variable 'top_complex_fluorescence' is zero: mean = [ 1.88409656e+02 1.17280879e+02 1.47675133e+09 6.30118914e+01
1.90358168e+03 6.30116323e+04 5.06862069e+01 3.10058043e+06
4.98020568e+01 5.21666068e+01 4.96382613e+01 5.23060970e+01], stddev = 4.64824130342
try 2:
File "/Users/choderaj/miniconda3/lib/python3.5/site-packages/assaytools-0.2.0-py3.5.egg/assaytools/pymcmodels.py", line 81, in tau
raise Exception("'np.sqrt(np.log(1.0 + (stddev/mean)**2))' for variable '%s' is zero: mean = %s, stddev = %s" % (name, mean, stddev))
Exception: 'np.sqrt(np.log(1.0 + (stddev/mean)**2))' for variable 'top_complex_fluorescence' is zero: mean = [ 1.93209203e+02 1.55450974e+02 3.71964662e+05 1.93762290e+09
6.37610120e+07 6.58578412e+01 2.89705585e+05 6.16018394e+01
2.09628235e+04 3.79087228e+12 1.69264036e+03 6.06388630e+01], stddev = 10.1443560014
I don't know if this will give you a clue, but when I analyse the whole plate only the binding experiment in the first row gives this error. This happened with two different xml files. Those values must not be calculated for other rows.
I think I've solved the problem.
The expression
np.sqrt(np.log(1.0 + (stddev/mean)**2))
ends up being zero to numerical precision when x = (stddev/mean)**2 is very small.
I've replaced this expression with a Taylor series expansion that achieves much higher numerical accuracy when x is small.
Fixed by #93