Calibration models

[bas-rustenburg]

Actually, do you mean the mechanical heat, or a different term that only depends on the fraction x? Trying to remember what we discussed last week.

I suppose the mechanical heat only depends on the volume of the injection... maybe I could model that as

$H_{mech} * v_{inj}/V_0?

[jchodera]

I think the mechanical heat should be just volume dependent: `$H{mech} * v{inj}$, without the $V_0$.

But you're right, the H_0 * x term may not have physical meaning---this form is essentially a Taylor expansion for Q(x), which may be all we need right now.

[bas-rustenburg]

Not sure if your implementation of the model above is correct. In the plot you give, the largest heat effect is for the first (tiny) injection. Perhaps you are computing the mole fraction x incorrectly? The mole fraction of titrant should be really tiny for the first injection.

I've played with my mole fraction function a bit, but I think it is correct. https://gist.github.com/anonymous/40a3e254938d074700df

[bas-rustenburg]

Adding a H_0 * x term does not improve the model much however.

[bas-rustenburg]

I've played with my mole fraction function a bit, but I think it is correct.

That said, does not mean the computed mole fraction is correct. Input might still be off. Checking the model.

[bas-rustenburg]

log scale is base 10

Initial mole fraction is definitely much smaller.

[jchodera]

Can you point to the code lines in question? The mole fraction plot looks right, but I don't understand why this leads to the first model heat q_n being the largest heat, despite it being the smallest change in mole fraction. That seems wrong.

[bas-rustenburg]

Heat is computed here, molefraction happens right above it. https://github.com/bas-rustenburg/bayesian-itc/blob/561b18db51ba4abc3dd0f76d2a16faef56f4620f/bitc/heats.py#L198-L202

I tried plugging in various different initial injection sizes into the model, plotting the expected heats (while fixing other parameters). I find the behavior quite striking...

See this ipython notebook where I play with the setup. https://gist.github.com/b0c93e29a63970193f1a

[bas-rustenburg]

Just noticed, havent pushed the model using dH * x yet, so that model just uses x * (1-x) * chi *.... let me push that.

[bas-rustenburg]

New code: https://github.com/choderalab/bayesian-itc/blob/2466cc442cbb4d1e12794bf091035986d8291560/bitc/heats.py#L251-L255

[bas-rustenburg]

I think I found a bug: https://github.com/choderalab/bayesian-itc/blob/2466cc442cbb4d1e12794bf091035986d8291560/bitc/heats.py#L255

- q_n[n-1]

This is wrong.. should have been

- (q_n[n-1]  q_n[n-2]  q_n[n-3] ....  q_n[0])

[bas-rustenburg]

Or something like that.. unless I'm mistaken.. since

is the cumulative heat...

[jchodera]

Maybe it's less confusing if you define the heat function Q_n[n] explicitly and compute

q_n[n] = Q_n[n] - Q_n[n-1] + H_mech

?

[bas-rustenburg]

Might be a good idea.

[bas-rustenburg]

This is already looking better, but now it just shifts the curve:

Implementation:

# Compute expected injection heats.
    # q_n_model[n] is the expected heat from injection n
    q_n = numpy.zeros([N])
    Q_n = 0.0
    # Instantaneous injection model (perfusion)
    # first injection
    # From units of cal/mole to ucal
    # page 149 of Atkins physical chemistry, 8th edition, eq. 5.30 H^E = n* chi *RT * x_a * x_b

    for n in range(N):
        # subsequent injections
        # From units of cal/mole to ucal
        Q_old = Q_n
        Q_n += 1.e9 * V0 * (DeltaH * ((X_frac[n])) + (buffer_concentration + Xn[n]) * chi * kt * (X_mfrac[n])) +  H_0 
        q_n[n] = Q_n - Q_old
    return q_n

Still has the largest heat effect though..

[bas-rustenburg]

Oh.. I'm probably treating H_mech wrong now.

[jchodera]

Is H_0 the same as H_mech? It seems to be really large. Your distribution says it is -4, but I'm not sure what units those are.

[jchodera]

Also, for this:

(buffer_concentration + Xn[n]) * chi * kt * (X_mfrac[n]))

Do you really mean to mix Xn[n] and X_mfrac[n]?

[bas-rustenburg]

Model shape looks a lot better now, although I feel that the heats plateau off too quickly still. May still be a bug somewhere, or maybe this is the best the model can do..

[pgrinaway]

great fool request

[jchodera]

Nicely done.

[bas-rustenburg]

Yep, still a bug somewhere, My molefractions look like: [ 3.78391307e-08 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00] .

Really happy to have pycharm debugger tools, would have taken a while to figure this out otherwise.

[bas-rustenburg]

The concentrations are already messed up.. weird..

Xn = [ 1.94288347e-06 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]

de:bug:ing

[bas-rustenburg]

Wrong indentation (without error) due to copy pasting between files.. sometimes I miss ruby...

[bas-rustenburg]

Edit: itchy clicking finger... here's the explanation.

Here is the data from the current model. Distributions look a bit weird. The H_mech now has units of ucal/liter so its excessively large. I think I will tweak the units a little bit.

# Compute expected injection heats.
    # q_n_model[n] is the expected heat from injection n
    q_n = numpy.zeros([N])
    Q_n = 0.0
    # Instantaneous injection model (perfusion)
    # first injection
    # From units of cal/mole to ucal
    # page 149 of Atkins physical chemistry, 8th edition, eq. 5.30 H^E = chi *nk_bT * x_a * x_b
    for n in range(N):
        # subsequent injections
        # From units of cal/mole to ucal
        Q_old = Q_n
        Q_n = 1.e9 * (DeltaH * (X_frac[n]) + chi * nkt * (X_mfrac[n]))
        q_n[n] = Q_n - Q_old + DeltaVn[n] * H_0
    return q_n

Model is not perfect yet.

[bas-rustenburg]

Hurray, more bugs!

my prior for H_mech is not defined correctly. Somehow a tuple ends up getting wrapped in a quantity..with nonsensical units too..

lower=(<Quantity(-3906835.54998, 'microcal / cal')>,)
upper=1552754.52936 microcal / liter
value=-3906835.54998 microcal / liter

[bas-rustenburg]

That bug is fixed. Seems that I introduced it when changing the units of H_mech.

Didn't change much.

[jchodera]

@bas-rustenburg : Can you take a look and see if this should be cleaned up and merged, or has this diverged so much from the current codebase that we should close the PR?

[bas-rustenburg]

From the looks of it, it has diverged too much. I also think in the current state has bugs and might be cleaner to revisit with a clean implementation if we end up needing this. I will leave the branch intact so we may salvage some stuff from it if we pick this up in the future, but the PR should be closed.