VAMP: properly handle IEEE exceptional values produced by the sampling function

[Juergen Reuter]

Using the input file (much the same as for #232)

process proc3 = g, g -> g, g
compile
sqrts = 7 TeV
integrate (proc3) { iterations = 10:10000 }

produces

| Loading process library 'processes'
| Process 'proc3': updating configuration
sqrts =    7000.0000000000000     
| Integrating process 'proc3'
| Generating phase space configuration ...
| ... found 5 phase space channels.
| Wrote phase-space configuration file 'proc3.phs'.
Warning: No cuts have been defined.
| Using partonic energy as event scale.
| iterations = 10:10000
| Creating grids
| 10000 calls, 5 channels, 2 dimensions, 20 bins, stratified = T
|=============================================================================|
| It      Calls  Integral[fb]  Error[fb]   Err[%]    Acc  Eff[%]   Chi2 N[It] |
|=============================================================================|
   1      10000  7.2964767E+11  4.13E+11   56.55   56.55*   0.07
   2      10000  7.3067373E+11  4.13E+11   56.47   56.47*   0.05
   3      10000  7.3072496E+11  4.13E+11   56.47   56.47*   0.05
   4      10000  7.3083000E+11  4.13E+11   56.46   56.46*   0.03
   5      10000  7.3083016E+11  4.13E+11   56.46   56.46*   0.05
   6      10000  7.3093205E+11  4.13E+11   56.45   56.45*   0.00
   7      10000  7.3093353E+11  4.13E+11   56.45   56.45*   0.01
   8      10000  7.3093361E+11  4.13E+11   56.45   56.45*   0.04
   9      10000  7.3093397E+11  4.13E+11   56.45   56.45*   0.01
  10      10000  7.3093398E+11  4.13E+11   56.45   56.45*   0.06
|=============================================================================|
  10     100000  7.3073737E+11  1.30E+11   17.86   56.47    0.06    0.00  10
|=============================================================================|
| Process 'proc3': time estimate for 10000 unweighted events = 0h 18m 12s
|-----------------------------------------------------------------------------|
| Contributing Helicity Combinations:     6 of    16
| Threshold: amp / avg >  0.22E-05 =  0.10E+11 * epsilon()

OK, there are no cuts, so it is clear that this academically partonic process cannot integrate very well, but at least this should be documented somehow.

[Juergen Reuter]

This should be rechecked with the NAG compiler. Basically the reason for it is that there are no more any magical cuts as in WHIZARD 1. But nevertheless it should not run into NaNs?. And actually here is an example with NaNs?:

| Loading process library 'processes'
| Process 'proc3': updating configuration
sqrts =    7000.0000000000000     
| Integrating process 'proc3'
| Generating phase space configuration ...
| ... found 5 phase space channels.
| Wrote phase-space configuration file 'proc3.phs'.
Warning: No cuts have been defined.
| Using partonic energy as event scale.
| iterations = 10:10000
| Creating grids
| 10000 calls, 5 channels, 2 dimensions, 20 bins, stratified = T
|=============================================================================|
| It      Calls  Integral[fb]  Error[fb]   Err[%]    Acc  Eff[%]   Chi2 N[It] |
|=============================================================================|
   1      10000  4.9068761E+11  1.55E+11   31.61   31.61*   0.10
   2      10000  4.9245669E+11  1.55E+11   31.50   31.50*   0.05
   3      10000  4.9538494E+11  1.55E+11   31.31   31.31*   0.04
   4      10000  4.9549916E+11  1.55E+11   31.30   31.30*   0.06
   5      10000  4.9549917E+11  1.55E+11   31.30   31.30*   0.05
   6      10000  4.9550026E+11  1.55E+11   31.30   31.30*   0.00
   7      10000  4.9550113E+11  1.55E+11   31.30   31.30*   0.00
   8      10000  4.9550119E+11  1.55E+11   31.30   31.30*   0.01
   9      10000  4.9550119E+11  1.55E+11   31.30   31.30*   0.06
  10      10000  4.9550119E+11  1.55E+11   31.30   31.30*    NaN
|=============================================================================|
  10     100000  4.9470337E+11  4.90E+10    9.91   31.35     NaN    0.00  10
|=============================================================================|
| Process 'proc3':
|    time estimate for 10000 unweighted events = 0h 14m 14s
|-----------------------------------------------------------------------------|

[Juergen Reuter]

The NaN comes from the vamp_grids%grids%mu data type. I will debug it once VAMP is there in its final form.

[Juergen Reuter]

If there is no comment until 1800 GST+0100 about this, this ticket will be moved to 2.0.1 and the trunk closed w.r.t. to this task.

[ohl]

Sorry, this will have to be fixed in milestone:2.0.1 (where #2 is due anyway).

[ohl]

Currently, I can NOT reproduce this with gfortran/64bit.

[Juergen Reuter]

Try around a little bit: it is easy to reproduce, e.g. set the iterations to 30:1000.

[ohl]

I see.

[Juergen Reuter]

TO, maybe if you have a look at the calculation of the efficiency, you could already have in mind the possibility of negative weights (i.e. indefinite integration functions).

[ohl]

The NaNs in question are not generated by VAMP, but are produced by the function being sampled. As mentioned above, this is caused by missing cuts and integrating a singular region.

There is no reason to catch these NaNs in the efficiency calculation.

At the moment they are reported as a negative weight (since "NaN > 0" is false per IEEE standard), but WHIZARD doesn't check the corresponding exception. Once the optional IEEE modules of Fortran 2003 are widely available, I can check portably and reliably for NaNs and !Infinity and handle both properly. I could test for "x != x" but that's a hack.

I leave this open for milestone:v2.1.

[Juergen Reuter]

Do we check for IEEE functions in the configure, actually? So we could screen over the status within the compiler world. As TO started to implement IEEE checks in VAMP, maybe we could boost this ticket a little bit.

[ohl]

wait for proper IEEE/2003 implementations ...

[Juergen Reuter]

-ffpe-trap=invalid,zero,overflow is the gfortran compiler to do this. This should be mentioned in the manual.

[Juergen Reuter]

It has been agreed upon to do mention this in the manual, but not to set it as a default. Will do that soon.

[Juergen Reuter]

Part of the manual in r5736. Closing.