Re: r25483 - /trunk/test_suite/system_tests/relax_disp.py -- September 01, 2014 - 11:24

Hi,

I just ran the following script:

-----
from math import exp
from numpy import array

times = array([ 0.7,  1. ,  0.8,  0.4,  0.9])
I = array([ 476.76174875,  372.43328777,  454.20339981,  656.87936253,
 419.16726341])
errors = array([  9.48032653,  11.34093541,   9.35149017,
10.84867928,  12.17590736])

R = - 500.
I0 = 1000.
params = [R, I0]

I_bc = []
chi2 = 0.0
for i in range(len(I)):
    I_bc.append(I0*exp(-R*times[i]))
    chi2 += (I[i] - I_bc[i])**2 / errors[i]**2
print(chi2)
-----

And I see:

-----
aaa.py:16: RuntimeWarning: overflow encountered in double_scalars
  chi2 += (I[i] - I_bc[i])**2 / errors[i]**2
inf
-----

However optimisation will never reach this R = -500 value.  When R =
-5, the chi2 value is 270307925.086.  When -0.1, chi2 = 16974.5262241.
Therefore any optimisation algorithm will never let R go up to such a
high negative value to start with.  This is not an issue that needs to
be worried about.

Regards,

Edward



On 31 August 2014 14:33,  <tlinnet@xxxxxxxxxxxxx> wrote:
> Author: tlinnet
> Date: Sun Aug 31 14:33:33 2014
> New Revision: 25483
>
> URL: http://svn.gna.org/viewcvs/relax?rev=25483&view=rev
> Log:
> Inserted systemtest Relax_disp.test_finite_value, to illustrate the return 
> of inf from C-code exponential, when R is negative.
>
> This can be an issue, if minfx takes a wrong step when no constraints are 
> implemented.
>
> bug #22552(https://gna.org/bugs/index.php?22552): Chi2 value returned from 
> C-code Curve-fitting return 0.0 for wrong parameters -> Expected influence 
> on Monte-Carlo sim.
>
> Modified:
>     trunk/test_suite/system_tests/relax_disp.py
>
> Modified: trunk/test_suite/system_tests/relax_disp.py
> URL: 
> http://svn.gna.org/viewcvs/relax/trunk/test_suite/system_tests/relax_disp.py?rev=25483&r1=25482&r2=25483&view=diff
> ==============================================================================
> --- trunk/test_suite/system_tests/relax_disp.py (original)
> +++ trunk/test_suite/system_tests/relax_disp.py Sun Aug 31 14:33:33 2014
> @@ -44,6 +44,15 @@
>  from status import Status; status = Status()
>  from test_suite.system_tests.base_classes import SystemTestCase
>
> +# C modules.
> +if dep_check.C_module_exp_fn:
> +    from specific_analyses.relax_fit.optimisation import func_wrapper, 
> dfunc_wrapper, d2func_wrapper
> +    from target_functions.relax_fit import jacobian, jacobian_chi2, setup
> +    # Call the python wrapper function to help with list to numpy array 
> conversion.
> +    func = func_wrapper
> +    dfunc = dfunc_wrapper
> +    d2func = d2func_wrapper
> +
>
>  class Relax_disp(SystemTestCase):
>      """Class for testing various aspects specific to relaxation dispersion 
> curve-fitting."""
> @@ -65,7 +74,8 @@
>                  
> "test_bug_21344_sparse_time_spinlock_acquired_r1rho_fail_relax_disp",
>                  "test_estimate_r2eff_err",
>                  "test_estimate_r2eff_err_auto",
> -                "test_estimate_r2eff_err_methods"
> +                "test_estimate_r2eff_err_methods",
> +                "test_finite_value",
>                  "test_exp_fit",
>                  "test_m61_exp_data_to_m61",
>                  "test_r1rho_kjaergaard_auto",
> @@ -3288,6 +3298,26 @@
>          self.assert_('test' not in cdp.clustering)
>          self.assertEqual(cdp.clustering['free spins'], [':2@N'])
>          self.assertEqual(cdp.clustering['cluster'], [':1@N', ':3@N'])
> +
> +
> +    def test_finite_value(self):
> +        """Test return from C code, when parameters are wrong.  This can 
> happen, if minfx takes a wrong step."""
> +
> +        times = array([ 0.7,  1. ,  0.8,  0.4,  0.9])
> +        I = array([ 476.76174875,  372.43328777,  454.20339981,  
> 656.87936253,  419.16726341])
> +        errors = array([  9.48032653,  11.34093541,   9.35149017,  
> 10.84867928,  12.17590736])
> +
> +        scaling_list = [1.0, 1.0]
> +        setup(num_params=2, num_times=len(times), values=I, sd=errors, 
> relax_times=times, scaling_matrix=scaling_list)
> +
> +        R = - 500.
> +        I0 = 1000.
> +        params = [R, I0]
> +
> +        chi2 = func(params)
> +
> +        print("The chi2 value returned from C-code for R=%3.2f and 
> I0=%3.2f, then chi2=%3.2f"%(R, I0, chi2))
> +        self.assertNotEqual(chi2, inf)
>
>
>      def test_hansen_catia_input(self):
>
>
> _______________________________________________
> relax (http://www.nmr-relax.com)
>
> This is the relax-commits mailing list
> relax-commits@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-commits