Subject: Re: [nmr_sparky] modifying relax.py
From: Tom Goddard
Date: Feb 6, 2009

Previous: 556

Hi Asli,

The code assumes that the second fitting parameter (b if you have
a,b,c) is the exponent coefficient for purposes of computing the time
constant. That assumption is built into the calculate_decay_time()
routine. The output text us determined by the peak_line() routine. The
optimization uses a gradient descent method to find the local minimum of
the sum of squares from the fit curve to the data points. The initial
guess has to be close to the desired answer if there are multiple minima.

Tom

Asli Ertekin wrote:

Hi,
Thanks for the help. Now I can get nice fits. I changed the model to
a*exp(b*t)+c, and it understood this equation better, for some reason.
Now I get good fits, however this depends a lot on the initial
guesses, especially for the exponential factor. I could not implement
the command you wrote as t_values = [t for t,h in th_values]. I wrote
it as t_values = [t for t,h in th_pairs]. Some resides had good fit
some did not. I dont know how to correct for that. When I start with
initial guess 0.5 I get good fits as judged by my eyes.
Have you ever encountered such problems. The fact that the fits are so
much dependent on tiny changes in initial guesses makes me suspicious.
Another thing: Now the results table lists either 2 parameters in
T-decay column or it lists 0, even though there is a fit for the
corresponding residue. I could not see where the output is written, so
I cannot manipulate it. Do you have any recommendations for that??

Thanks a lot for your help.

Asli


On Wed, Feb 4, 2009 at 4:50 PM, Tom Goddard goddard@...
mailto:goddard@... wrote:

Hi Asli,

I have modified relax.py for other functions in the past and it
does work. One problem with your code is that the
parameter_count() routine should return 3, not 2 since you have 3
parameters a, b, c. The other thing I would be concerned with is
providing good initial values. Maybe max(h_values) and
min(h_values) would be better initial guesses for a,b. And for c
which you start at 1e-6 it might be better to use max(t_values)
with t_values = [t for t,h in th_values].

Tom


ertekinler wrote:


Hi,
I want to use relax.py (rh command) for T1 calculation.
However the
a*exp(b*t) model is not exactly right for T1, so I want to
modify the
python code for this. I changed where the model is defined and
where
(I believe) the initial guesses are defined as follows:

class a_exp_bx_function:

def parameter_count(self):
return 2
def value(self, params, x):
a, b ,c = params
# Multiply overflows in fit error calculation cause core dump on
DEC Alpha
cx = pyutil.bound_value(c*x, -50, 50)
e = math.exp(cx)
return a-(a-b)*e
def parameter_derivatives(self, params, x):
a, b, c = params
e = math.exp(c*x)
return (1-e,e,-a*x*e+b*x*e)

def fit_exponential(self, th_pairs):

fit_th_pairs = filter(lambda th: th[1] != None, th_pairs)
if len(fit_th_pairs) 2:
return None

exp_func = a_exp_bx_function()

h_values = map(lambda th: th[1], fit_th_pairs)
h_middle = .5 * (max(h_values) + min(h_values))
params = (h_middle, h_middle,0)
trials = self.error_estimate_trials

fit = curvefit.least_squares_fit(fit_th_pairs, exp_func,
params, 1e-6)
fit.estimate_parameter_error(fit.point_sd, trials)

self.calculate_decay_time(fit)

return fit

However, when I run this code in sparky I get a straight line (I
guess) at h_middle. I dont know anything about Python so I
can not
figure out where the problem is.
Is there anyone who had tried this before.
And maybe more important question is, can this code be modified to
give reasonable fit for y=a-(a-b)*exp(ct) kind of model, which
cannot
be linearized by taking simple log.

Thanks

Asli Ertekin