← Back to branch summary

~ubuntu-branches/debian/sid/gsl/sid

« back to all changes in this revision

Viewing changes to monte/vegas.c

  • Committer: Bazaar Package Importer
  • Author(s): Dirk Eddelbuettel
  • Date: 2008-12-16 06:17:55 UTC
  • mfrom: (1.3.2 upstream) (3.1.15 jaunty)
  • Revision ID: james.westby@ubuntu.com-20081216061755-9la7p0qwrhopk8pl
Tags: 1.12+dfsg-1
* New upstream version released today

* doc/*: As before, removed the 'non-free' documentation to create a 
  source package that complies with Debian's interpretation of what is free. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
monte/vegas.c
71
71
#define NEW_COORD(s,i) ((s)->xin[(i)])
72
72
#define VALUE(s,i,j) ((s)->d[(i)*(s)->dim + (j)])
73
73
 
 
74
#define USE_ORIGINAL_CHISQ_FORMULA 0
 
75
 
74
76
/* predeclare functions */
75
77
 
76
78
typedef int coord;
149
151
      state->chi_sum = 0;
150
152
      state->it_num = 1;
151
153
      state->samples = 0;
 
154
      state->chisq = 0;
152
155
    }
153
156
 
154
157
  if (state->stage <= 2)
212
215
  cum_int = 0.0;
213
216
  cum_sig = 0.0;
214
217
 
215
 
  state->chisq = 0.0;
216
 
 
217
218
  for (it = 0; it < state->iterations; it++)
218
219
    {
219
220
      double intgrl = 0.0, intgrl_sq = 0.0;
220
 
      double sig = 0.0;
221
 
      double wgt;
 
221
      double tss = 0.0;
 
222
      double wgt, var, sig;
222
223
      size_t calls_per_box = state->calls_per_box;
223
224
      double jacbin = state->jac;
224
225
      double *x = state->x;
231
232
      
232
233
      do
233
234
        {
234
 
          double m = 0, q = 0;
 
235
          volatile double m = 0, q = 0;
235
236
          double f_sq_sum = 0.0;
236
237
 
237
238
          for (k = 0; k < calls_per_box; k++)
238
239
            {
239
 
              double fval, bin_vol;
 
240
              volatile double fval;
 
241
              double bin_vol;
240
242
 
241
243
              random_point (x, bin, &bin_vol, state->box, xl, xu, state, r);
242
244
 
243
245
              fval = jacbin * bin_vol * GSL_MONTE_FN_EVAL (f, x);
244
246
 
245
 
              /* recurrence for mean and variance */
 
247
              /* recurrence for mean and variance (sum of squares) */
246
248
 
247
249
              {
248
250
                double d = fval - m;
259
261
 
260
262
          intgrl += m * calls_per_box;
261
263
 
262
 
          f_sq_sum = q * calls_per_box ;
 
264
          f_sq_sum = q * calls_per_box;
263
265
 
264
 
          sig += f_sq_sum ;
 
266
          tss += f_sq_sum;
265
267
 
266
268
          if (state->mode == GSL_VEGAS_MODE_STRATIFIED)
267
269
            {
272
274
 
273
275
      /* Compute final results for this iteration   */
274
276
 
275
 
      sig = sig / (calls_per_box - 1.0)  ;
 
277
      var = tss / (calls_per_box - 1.0)  ;
276
278
 
277
 
      if (sig > 0) 
 
279
      if (var > 0) 
278
280
        {
279
 
          wgt = 1.0 / sig;
 
281
          wgt = 1.0 / var;
280
282
        }
281
283
      else if (state->sum_wgts > 0) 
282
284
        {
289
291
        
290
292
     intgrl_sq = intgrl * intgrl;
291
293
 
 
294
     sig = sqrt (var);
 
295
 
292
296
     state->result = intgrl;
293
 
     state->sigma  = sqrt(sig);
 
297
     state->sigma  = sig;
294
298
 
295
299
     if (wgt > 0.0)
296
300
       {
 
301
         double sum_wgts = state->sum_wgts;
 
302
         double wtd_int_sum = state->wtd_int_sum;
 
303
         double m = (sum_wgts > 0) ? (wtd_int_sum / sum_wgts) : 0;
 
304
         double q = intgrl - m;
 
305
 
297
306
         state->samples++ ;
298
307
         state->sum_wgts += wgt;
299
308
         state->wtd_int_sum += intgrl * wgt;
302
311
         cum_int = state->wtd_int_sum / state->sum_wgts;
303
312
         cum_sig = sqrt (1 / state->sum_wgts);
304
313
 
 
314
#if USE_ORIGINAL_CHISQ_FORMULA
 
315
/* This is the chisq formula from the original Lepage paper.  It
 
316
   computes the variance from <x^2> - <x>^2 and can suffer from
 
317
   catastrophic cancellations, e.g. returning negative chisq. */
305
318
         if (state->samples > 1)
306
319
           {
307
320
             state->chisq = (state->chi_sum - state->wtd_int_sum * cum_int) /
308
321
               (state->samples - 1.0);
309
322
           }
 
323
#else
 
324
/* The new formula below computes exactly the same quantity as above
 
325
   but using a stable recurrence */
 
326
         if (state->samples == 1) {
 
327
           state->chisq = 0;
 
328
         } else {
 
329
           state->chisq *= (state->samples - 2.0);
 
330
           state->chisq += (wgt / (1 + (wgt / sum_wgts))) * q * q;
 
331
           state->chisq /= (state->samples - 1.0);
 
332
         }
 
333
#endif
310
334
       }
311
335
     else
312
336
       {
318
342
      if (state->verbose >= 0)
319
343
        {
320
344
          print_res (state,
321
 
                     state->it_num, intgrl, sqrt (sig), cum_int, cum_sig,
 
345
                     state->it_num, intgrl, sig, cum_int, cum_sig,
322
346
                     state->chisq);
323
347
          if (it + 1 == state->iterations && state->verbose > 0)
324
348
            {