Yes, this is crap. I just fixed that this morning and glad you sent out
this email so I don't have to explain it all. I've attached the amended
file method.c (which for the others belongs in
ccpnmr1.0/c/ccpnmr/analysis). You need to re-compile. I'm not sure what
planet I was on when those relevant lines got commented out.
In any case I would recommend against using this method (the "parabolic
fit") but would instead recommend one of the box sum methods (again, this
can be set in Crosspeaks -> Peak Find Parameters under "Other parameters",
the "volume method"). This time around I left the "parabolic fit" as the
default, but I'm going to switch that around next time.
Wayne
On Thu, 2 Jun 2005, Bruce D. Ray wrote:
> >Right, I think I know what is causing the problem. The new peaks are
> >coming out with infinite volume (you can check in the seleced peaks
> >dialog). And for some obscure reason which I can't fathom (especially
> >since I wrote that code) the point peak search algorithm is using the
> >volume to scale the search region. I'm going to change that but the real
> >point is that the volume should not be coming out infinite. And I don't
> >know yet why it is in 1.0.3 on the Mac (I cannot reproduce this on Linux).
> >It's quite possibly a stupid missing initialisation (or similar) in the C
> >world.
>
>
> I believe the volume determination is in method.c the diff from analysis1.02
> of which reads:
>
> 127a128
> > /*
> 128a130,131
> > */
> > float volume = 1, b;
> 172a176
> > /*
> 173a178
> > */
>
>
> Now, in method.c for analysis1.03, the relevant code sections read:
>
> 128 /*
> 129 float volume = 1, b, c;
> 130 */
> 131 float volume = 1, b;
>
>
> 176 /*
> 177 c = fit_center3(vm, y, vp, &b, TRUE);
> 178 */
> 179 /*
> 180 c = fit_center3(ym[i], y, yp[i], &b, TRUE);
> 181 */
> 182 /*
> 183 a = y / (float) exp((double) (- b * c * c));
> 184 volume *= a;
> 185 */
> 186 /* 13 Aug 2004: don't get what the hell this was doing
> 187 volume += b * c * c;
> 188 think it should be below instead
> 189 */
> 190 /* below should be fit from y = a exp -b(x-c)^2, with x = -1, 0, 1
> 191 do not know what proper x scale is so overall result is out */
> 192 /* integral if y is (a sqrt(pi/b)) but just ignore a and use y */
> 193 volume *= sqrt(PI / (double) ABS(b));
>
>
> This would be all well and good if I could just find someplace in the function
> fit_volume_gaussian3_method above this piece of code where b was set to some
> value. As it is, I'm not quite clear on what the value of the divisor, ABS(b)
> is. If this is working under Linux, then I must be missing something because
> the entire diff of method.c is so simple.
>
>
> Sincerely,
>
>
> --
> Bruce D. Ray, Ph.D.
> Associate Scientist, and Operations Director
> NMR Center
> IUPUI
> Physics Dept.
> 402 N. Blackford St.
> Indianapolis, IN 46202-3273
>
/*
======================COPYRIGHT/LICENSE START==========================
method.c: Part of the CcpNmr Analysis program
Copyright (C) 2004 Wayne Boucher and Tim Stevens (University of Cambridge)
=======================================================================
This file contains reserved and/or proprietary information
belonging to the author and/or organisation holding the copyright.
It may not be used, distributed, modified, transmitted, stored,
or in any way accessed, except by members or employees of the CCPN,
and by these people only until 31 December 2005 and in accordance with
the guidelines of the CCPN.
A copy of this license can be found in ../../../license/CCPN.license.
======================COPYRIGHT/LICENSE END============================
for further information, please contact :
- CCPN website (http://www.ccpn.ac.uk/)
- email: [log in to unmask]
- contact the authors: [log in to unmask], [log in to unmask]
=======================================================================
If you are using this software for academic purposes, we suggest
quoting the following references:
===========================REFERENCE START=============================
R. Fogh, J. Ionides, E. Ulrich, W. Boucher, W. Vranken, J.P. Linge, M.
Habeck, W. Rieping, T.N. Bhat, J. Westbrook, K. Henrick, G. Gilliland,
H. Berman, J. Thornton, M. Nilges, J. Markley and E. Laue (2002). The
CCPN project: An interim report on a data model for the NMR community
(Progress report). Nature Struct. Biol. 9, 416-418.
Wim F. Vranken, Wayne Boucher, Tim J. Stevens, Rasmus
H. Fogh, Anne Pajon, Miguel Llinas, Eldon L. Ulrich, John L. Markley, John
Ionides and Ernest D. Laue. The CCPN Data Model for NMR Spectroscopy:
Development of a Software Pipeline. Accepted by Proteins (2004).
===========================REFERENCE END===============================
*/
#include "method.h"
#define SMALL_VALUE (1.0e-4)
/* u, v, w are three y values in a row for points x-1, x and x+1 */
/* returns the offset of the center from x */
/* implicitly assumes that if v > 0 then v >= u,w and if v < 0 then v <= u,w */
/* b is used for volume calculation only */
static float fit_center3(float u, float v, float w, float *b, Bool use_log)
{
float c, d;
Bool is_positive;
if (v > 0)
is_positive = TRUE;
else
is_positive = FALSE;
if (use_log)
{
if (!is_positive)
{
u = -u;
v = -v;
w = -w;
}
if ((u > 0) && (w > 0))
{
u = (float) log((double) u);
v = (float) log((double) v);
w = (float) log((double) w);
}
/* otherwise just do non-log fit */
}
d = 0.5 * ABS(2*v - u - w);
if (is_positive)
*b = d;
else
*b = -d;
if (d > SMALL_VALUE)
{
c = 0.25 * ABS(w - u) / d;
if (is_positive)
{
if (w < u)
c = -c;
}
else
{
if (w > u)
c = -c;
}
c = MAX(c, -0.499);
c = MIN(c, 0.499);
}
else
{
c = 0;
}
return c;
}
#define FUDGE_MULT 0.01
float fit_volume_gaussian3_method(int ndim, float y, float *ym, float *yp)
{
int i;
/*
float volume, a, b, c;
*/
/*
float volume, b, c;
*/
float volume = 1, b, c;
float vm, v, vp;
/* change to log way of calculating volume to try and exclude infinities */
/* 13 Aug 2004: do not know what the heck the below was trying to do
so change back to not using logs, only problem here is overall scale */
for (i = 0; i < ndim; i++)
{
vm = ym[i];
v = y;
vp = yp[i];
/*
printf(" i = %d, ym = %3.2f, y = %3.2f, yp = %3.2f\n", i, vm, y, vp);
*/
if (v > 0)
{
if ((vm > v) || (vp > v))
{
y = MAX(vm, vp);
vm = vp = v;
}
else
{
vm = MAX(vm, FUDGE_MULT*y);
vp = MAX(vp, FUDGE_MULT*y);
}
}
else
{
if ((vm < v) || (vp < v))
{
y = MIN(vm, vp);
vm = vp = v;
}
else
{
vm = MIN(vm, FUDGE_MULT*y);
vp = MIN(vp, FUDGE_MULT*y);
}
}
/*
printf(" i = %d, vm = %3.2f, y = %3.2f, vp = %3.2f\n", i, vm, y, vp);
*/
c = fit_center3(vm, y, vp, &b, TRUE);
/*
c = fit_center3(ym[i], y, yp[i], &b, TRUE);
*/
/*
a = y / (float) exp((double) (- b * c * c));
volume *= a;
*/
/* 13 Aug 2004: don't get what the hell this was doing
volume += b * c * c;
think it should be below instead
*/
/* below should be fit from y = a exp -b(x-c)^2, with x = -1, 0, 1
do not know what proper x scale is so overall result is out */
/* integral if y is (a sqrt(pi/b)) but just ignore a and use y */
volume *= sqrt(PI / (double) ABS(b));
/*
printf(" i = %d, b = %3.2f, volume = %3.2f\n", i, b, volume);
*/
}
/*
volume = pow((double) volume, (double) 1.0/ndim);
*/
/* below was used before 13 Aug 2004
volume /= ndim;
volume += (float) log((double) ABS(y));
/ protect against exp overflow and underflow /
volume = (float) exp((double) MAX(-21.0, MIN(volume, 21.0)));
volume *= pow(sqrt(PI), (double) ndim);
*/
volume *= y;
/*
printf("at end y = %3.2f, volume = %3.2f\n", y, volume);
*/
return volume;
}
void fit_center_parabolic_method(int ndim, float *center, float y,
float *ym, float *yp)
{
int i;
float b;
for (i = 0; i < ndim; i++)
center[i] = fit_center3(ym[i], y, yp[i], &b, FALSE);
}
void fit_center_gaussian3_method(int ndim, float *center, float y,
float *ym, float *yp)
{
int i;
float b;
for (i = 0; i < ndim; i++)
center[i] = fit_center3(ym[i], y, yp[i], &b, TRUE);
}
|