[0.9.8] Regression with mpfit test about gdl HOT 19 CLOSED

Interesting. On the same code IDL8 gives [-1.19891 , 47.1908 , 1.69077] !

from gdl.

And what CHI-SQUARE?

from gdl.

Ah I see, mpfit is from Markwardt.
Versions may differ. Strange.

IDL> test_mpfit
% Compiled module: TEST_MPFIT.
% Compiled module: MPFIT.
Iter 1 CHI-SQUARE = 2102302.5 DOF = 98
P(0) = 10.0000
P(1) = 45.0000
P(2) = 10.0000
Iter 2 CHI-SQUARE = 10848.010 DOF = 98
P(0) = -0.276637
P(1) = 44.1631
P(2) = 10.9445
Iter 3 CHI-SQUARE = 10380.288 DOF = 98
P(0) = -0.177669
P(1) = 81.9140
P(2) = -25.1137
(...)
Iter 15 CHI-SQUARE = 9812.0020 DOF = 98
P(0) = -1.19995
P(1) = 47.1943
P(2) = 1.69246
Iter 16 CHI-SQUARE = 9812.0000 DOF = 98
P(0) = -1.19891
P(1) = 47.1908
P(2) = 1.69077
Iter 16 CHI-SQUARE = 9812.0000 DOF = 98
P(0) = -1.19891
P(1) = 47.1908
P(2) = 1.69077
-1.19891 47.1908 1.69077
28.2524 31.9683 -4.26773
% Variable is undefined: .
% Execution halted at: $MAIN$
% Program caused arithmetic error: Floating underflow

from gdl.

% Variable is undefined: .

Is ".5" not a valid floating point number in IDL (but in GDL)?

from gdl.

No, not at all. There is something fussy here. The test_mpfit.pro in the GDL testsuite gives no errors.

from gdl.

Is there any progress with this? I am worrying to ignore this, since this seems to be a real regression.

from gdl.

As test_mpfit was not passing with idl6 and idl8 I suspected a trick. Besides, as many of the instrument packages that work OK with GDL do use mpfit, such a problem would have been reported almost immediately.

The point is, as "y" was written, it was no more a noisy gaussian. And the noise figure (err) was not representative of the variance of y. And was negative at some points. Note that the difference between 0.9.7 and 0.9.8 is that randomn() has changed and is now the same as in IDL.
Well, I believe that the following procedure, for which the test is OK both in IDL8+ and GDL, can be used for debian tests:

FUNCTION fitfunc, p, x=x, y=y, err=err
  gauss = p[0]+p[1]*exp(-(x-p[2])^2/(2*p[3]^2))
  return, (y - gauss)/err
END
PRO test_mpfit
  x = findgen(101)
  err = replicate(3.0,101) ; error estimate
  y = 25.*exp(-(x-32.)^2/(2*5.^2))+12.0+err*randomn(33,101)
  start_params = [12.,10.,45.,10.]
  functargs = {x:x, y:y, err:err}
  params = mpfit("fitfunc",start_params,functargs=functargs)
  ref = [12.032, 24.0875, 32.2197, 5.096]
  print, params, ref
  w=where(abs(params-ref) gt 0.5, count)
  if (count gt 0)  then exit, status=1
  exit, status=0
END

from gdl.

I discovered this post today, and will give a look ASAP

from gdl.

I am totally lost here ! the code in Debian (with same name that code in GDL testsuite) changed several times and at the same time the RANDOM code changed in GDL :(

GDL in the range June 2017 up to Sep 2018 was OK with current version mpfit 1.85+2017.01.03-2

now, new GDL RANDOM with fast DSFMT don't give same random values than IDL, then is not normal the test in Debian fails. GDL in the range Dec 2017 up to Sep 2018 was OK

from gdl.

warning : IDL 8.2.2 starts to have the Mersenne

and to have the Debian test OK, you should do :

export GDL_NO_DSFMT=1
then start current GDL version :)

from gdl.

Hi,
I see the mpfit test as too naive. The small number of noisy samples does not insure the test success.
It fails statistically 2 times over 10 with IDL6 and IDL8. This is quite evident as there are 100 points with variance 9, and we want for fit 4 parameters, that is similar of taking the rms of 25 values of variance 9 for any point and hoping to find 0.0 \pm 0.5 everytime. Of course GDL, with or without dSFMT, exhibits the same behaviour.
I will provide a more statistically robust test, but my question is, what do the Debian guys want to measure with this test? If it is a TAT test to check the results of MPFIT, then it is not OK also.

from gdl.

I would just like to see that the package works somehow. The test was provided from someone (sorry, can't trace it back in the moment); I agree that it is not very stable.
If you have something that works better -- just with a small constant array --, I would happily include it.

from gdl.

here is the new test:

FUNCTION fitfunc, p, x=x, y=y, err=err
  gauss = p[0]+p[1]*exp(-(x-p[2])^2/(2*p[3]^2))
  return, (y - gauss)/err
END
; this test checks that, statistically, mpfit fits a gaussian in a
; noisy distribution. The gaussian individual parameters MUST be
; correctly sampled by the distribution (e.g., the width must not bee
; too small wrt the number of points etc. In order not to depend on a
; particular noisy sample, we bootstrap a bit with 10 trials. Also, we
; do it on 2 sample lengths, short and long. Short sample is too few
; to insure a good convergence of mpfit for all random number
; realizations, so a 3-over-10 trial erros is tolerated. Long sample
; should be OK in approximately all cases. Failure of this test is
; still marginally possible, in this case: retry.
PRO test_mpfit
  np=[101,10001]
  ecount=fltarr(2)
  ; test
  for i=0,1 do begin
     npoints=np[i]
     x = findgen(npoints)
; true gaussian parameters
     gpar = [12.,32., npoints/2., npoints/16.]
; the function  
     y0 = gpar[1]*exp(-(x-gpar[2])^2/(2*gpar[3]^2))+gpar[0]
; the noisy function  
     for j=0,9 do begin
        err = 3*randomn(seed,npoints) ; random error
        y = y0+err
        ;plot,x,y,psym=3
        functargs = {x:x, y:y, err:err}
; blind start parameters
        est_base=mean(y);
        est_height=max(y)-est_base
        w=where(y-est_base gt est_height/3, count)
        if (count le 0) then message,"error estimating initial parameters"
        est_center=mean(x[w])
        dummy=moment(x[w],mdev=est_width) 
        start_params = [est_base,est_height,est_center, est_width]
;        print,start_params
        params = mpfit("fitfunc",start_params,functargs=functargs,/quiet)
        y = params[1]*exp(-(x-params[2])^2/(2*params[3]^2))+params[0]
        ;oplot,x,y
        print, abs(params-gpar)
        w=where(abs(params-gpar) gt 0.5, num)
;        if (num gt 0) then stop ; if you want to check
     end
     ecount[i]+=num
  end
  if ecount[0] lt 3 then ecount[0]=0 ; accept 3 error on 10 due to 3*randomn stats.
  w=where(ecount gt 0, errcount) 
  if (errcount gt 0) then message,/info,"test failed! bad!"
  if (errcount gt 0) then exit, status=1
  exit, status=0
END

from gdl.

Hi @olebole , the test above is probably not sufficiently naive anymore.
Was the Debian test intended to test solely MPFIT? in an anti-regression test apporach?
In this case the random values err can be provided as an array in the original test (unreadable but efficient).
Then the only variations on the result will be the solely dependent on mpfit (or, the arithmetics in GDL).
Need that?

from gdl.

Yes, something like that. I just want to see that mpfit can be called somehow from gdl and does not completely fail. That can be really naive.

from gdl.

The following fixes the "random" part of the test. IDL8 and GDL give the same result. Any change will be due to a change in MPFIT.

FUNCTION fitfunc, p, x=x, y=y, err=err
  gauss = p[0]+p[1]*exp(-(x-p[2])^2/(2*p[3]^2))
  return, (y - gauss)/err
END
PRO test_mpfit
  x = findgen(101)
  err = [2.45583,-3.25605, 5.87131,-1.73504, 2.27937,-0.48254,-1.51753,-0.93256, 1.12209, 1.65635, 1.28370,-3.57939, 4.03734,-0.17061, 1.32193,-0.30438,-0.21467,-2.46095,-2.34313,-0.88627, 1.51226,-2.06973,-3.97503,-0.07323,-5.06270,-3.67407, 3.76492,-0.55431, 2.76809, 2.23146,-2.30901,-1.99149,-4.16457, 2.09245,-0.71644,-3.27360, 0.84205,-4.13405, 5.02517, 3.07811,-2.83230, 3.39967, 0.24383, 0.85002, 6.02679, 0.94512, 0.88598,-6.16427,-1.20659,-0.38692,-0.12944,-2.47544,-2.79558, 0.39396, 1.71503,-0.87366,-3.12669,-1.27370, 3.74370,-2.54965,-1.25783, 0.62415,-0.27057,-0.97599,-3.04811,-0.21875,-2.03241, 1.32219,-0.41634, 4.63014, 1.16504,-0.75943, 3.12185,-2.25148, 3.99848, 1.31917, 1.54279, 1.44990, 1.52935, 0.00041,-2.96998,-0.87376, 4.16183, 1.99984, 1.48010,-3.42435, 1.37966,-0.21252,-0.74505,-1.38282,-0.92772, 5.10892, 2.45020,-0.64636, 2.03258,-0.77305, 0.05296,-7.53868, 5.39158,-0.78637,-1.50639] ; on shot error estimate
  y = 25.*exp(-(x-32.)^2/(2*5.^2))+12.0+err
  start_params = [12.,10.,45.,10.]
  functargs = {x:x, y:y, err:err}
  params = mpfit("fitfunc",start_params,functargs=functargs,/quiet)
  ref = [12.0004, 24.4652, 32.1227, 5.07643]
  print, params, ref, abs(ref-params)
  w=where(abs(params-ref) gt 0.0001, count)
  if (count gt 0)  then exit, status=1
  exit, status=0
END

from gdl.

With the GDL version in unstable, I get the error:

% Compiled module: TEST_MPFIT.
% Compiled module: MPFIT.
      12.0004      24.4658      32.1226      5.07634
      12.0004      24.4652      32.1227      5.07643
  5.72205e-06  0.000638962  9.91821e-05  8.67844e-05

Shall I just adjust the reference? Or increase the tolerance?

from gdl.

One could increase the tolerance of course. Previoulsy it was 0.5 😄
Let's put it a 10-2. Incidentally, your GDL values are the same as mine for an unoptimized version of GDL.
With the -DCMAKE_CXX_FLAGS_RELEASE='-O3 ' the results are idl-compatible.
Does this mean that the Debian version is not optimized?

from gdl.

OK, I increased the tolerance and uploaded the new version. Thanks!
GDL is compiled with -O2 which is the Debian default setting. See logfile. We could change that, but I don't see a big win here.

from gdl.