On Wed, 2006-08-09 at 15:40 +1000, Edward d'Auvergne wrote:
Here is my proposed solution to bug #6503
(https://gna.org/bugs/?func=detailitem&item_id=6503). I have used the
current 1.2 line as a reference and presented the diffs. Feel free to
propose better solutions.
The first thing is the catching of the NaN and Inf after the grid
search and minimisation. I propose to do the test within the
'self.minimise()' function of 'specific_fns/model_free.py' rather than
in the 'maths_fns/' or 'minimise/' code. Here is the diff:
===================================================================
--- specific_fns/model_free.py (revision 2527)
+++ specific_fns/model_free.py (working copy)
@@ -2374,6 +2374,14 @@
self.g_count = self.g_count + gc
self.h_count = self.h_count + hc
+ # Catch infinite chi-squared values.
+ if self.func == float('inf'):
+ raise RelaxInfError
+
+ # Catch chi-squared values of NaN.
+ if isnan(self.func):
+ raise RelaxNaNError
+
# Scaling.
if scaling:
self.param_vector =
matrixmultiply(self.scaling_matrix, self.param_vector)
===================================================================
As the function 'self.grid_search()' executes 'self.minimise()' this
single change will cover both the grid search and minimisation. The
function 'isnan()' currently non-existent but could be implemented as
Gary previously mentioned by checking the hex bit pattern of NaN.
If isnan() is to be implimented in this way, then a similar isinf() is
by far the best way to go. float('inf') is platform dependent, and
certainly will not work under Windows (I think I'm right in saying
windows has no valid string representation of inf or nan that would work
in that context).
This avoids fpconst.py and the switch from Numeric to Numpy. The
question is whether to use the RelaxErrors or set some warning? Would
this be alleviated by Gary's proposal of saving the program state just
prior to throwing the error? For example using the function
'self.save()' within 'generic_fns/state.py', printing the error
message, then quitting.
Saving state as part of the exception handling is a nice idea, but has
its limitations. Because the error could be thrown before, during or
after some change to the program state, the state which is saved will be
undefined and possibly inconsistent. As such a saved state may be of
limited use to most users. It would be a very valuable tool for
debugging, though.
The next change is the addition of new RelaxErrors:
===================================================================
--- errors.py (revision 2527)
+++ errors.py (working copy)
@@ -462,3 +462,17 @@
class RelaxInvalidColourError(BaseError):
def __init__(self, colour):
self.text = "The colour " + `colour` + " is invalid."
+
+
+ # Value errors.
+ ###############
+
+ # Infinity.
+ class RelaxInfError(BaseError):
+ def __init__(self, name):
+ self.text = "The invalid " + name + " floating point
value of infinity has occurred."
+
+ # NaN (Not a Number).
+ class RelaxNaNError(BaseError):
+ def __init__(self, name):
+ self.text = "The invalid " + name + " floating point
value of NaN (Not a Number) has occurred."
===================================================================
And as a preventative measure, the catching of zero length XH bond vectors:
===================================================================
--- generic_fns/pdb.py (revision 2527)
+++ generic_fns/pdb.py (working copy)
@@ -156,7 +156,11 @@
print "The PDB file " + `self.file_path` + "
cannot be found, no structures will be loaded."
return
+ # Test that the nuclei have been correctly set.
+ if self.heteronuc == self.proton:
+ raise RelaxError, "The proton and heteronucleus are set
to the same atom."
+
# Data creation.
################
@@ -256,10 +260,19 @@
# Get the heteronucleus position.
posX = pdb_res.atoms[self.heteronuc].position.array
- # Calculate the normalised vector.
+ # Calculate the XH bond vector.
vector = posH - posX
-
self.relax.data.res[self.run][j].xh_vect.append(vector /
sqrt(dot(vector, vector)))
+ # Normalisation factor.
+ norm_factor = sqrt(dot(vector, vector))
+
+ # Test for zero length.
+ if norm_factor == 0.0:
+ raise RelaxError, "The XH bond vector for
residue " + `self.relax.data.res[self.run][j].num` + " is of zero
length."
+
+ # Calculate the normalised vector.
+
self.relax.data.res[self.run][j].xh_vect.append(vector / norm_factor)
+
# Print out.
if self.print_flag:
if num_str > 1:
===================================================================
Finally, forcing termination of the backtracking line search (to break
the infinite loop):
===================================================================
--- minimise/line_search/backtrack.py (revision 2527)
+++ minimise/line_search/backtrack.py (working copy)
@@ -23,7 +23,7 @@
from Numeric import dot
-def backtrack(func, args, x, f, g, p, a_init=1.0, rho=0.5, c=1e-4):
+def backtrack(func, args, x, f, g, p, a_init=1.0, rho=0.5, c=1e-4,
max_iter=500):
"""Backtracking line search.
Procedure 3.1, page 41, from 'Numerical Optimization' by Jorge
Nocedal and Stephen J. Wright,
@@ -63,8 +63,9 @@
# Initialise values.
a = a_init
f_count = 0
+ i = 0
- while 1:
+ while i <= max_iter:
fk = apply(func, (x + a*p,)+args)
f_count = f_count + 1
@@ -73,3 +74,6 @@
return a, f_count
else:
a = rho*a
+
+ # Increment the counter.
+ i = i + 1
===================================================================
What do you think? I think I covered most of the issues. I've made a
branch of the 1.2 line located at
svn://svn.gna.org/svn/relax/branches/nan_catch_test which has all
these changes. Feel free to play around and modify this branch.
I'll have a bit more of a look, but it seems like a good fix.
Chris
Edward
_______________________________________________
relax (http://nmr-relax.com)
This is the relax-devel mailing list
relax-devel@xxxxxxx
To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-devel
Re: Proposed solution to bug #6503.