Re: VEGAS algorithm for monte carlo integration in LuaJIT2

[Francesco Abbate <francesco.bbt@...>]

I'm sorry about the noise but Lesley have found an error in the Lua
benchmark. The seed used was different and the Lua version was doing
more iterations to terminate.

Of course the benchmark is more meaningful if exactly the same seed is used.

With the correction in the benchmark plus a rewrite with a more direct
coding style the LuaJIT2 results are much better and it is now
*faster* than C. Here some figures:

LuaJIT2 -jon (new benchmark): 0m2.991s
LuaJIT2 -joff(new benchmark): 2m31.857s
C (gcc -O2): 0m3.514s

So that LuaJIT2 is actually ~ 15% faster then equivalent optimized C
code. The gain due to the JIT compilation can be quantified up to a
factor of ~ 50x which is really amazing.

The patch, to be chained to the previous LuaJIT2 patch, is in
attachment. May be someone else can check to verify this figures on a
different system.

So, I'm sorry for the noise but I'm glad to see that, once again
LuaJIT2 is very competitive even on the most demanding numerical
algorithms.

I take the opportunity also to announce that now GSL Shell is able to
save the plot in SVG format. This result is already consolidated in
the gsl-shell-2 branch and I'm quite proud of the good quality of the
SVG output.

Another important contribution was made by Benjamin von Ardenne, he
kindly accepted to reimplement the module about special functions
using the FFI interface. His implementation does also extend the
previous implementation to include all the special functions provided
by the GSL library.

I've also implemented also other useful features for plots, the
support for user-defined axis labels (useful for bar plots or similar)
and the support for plot legends. These latter feature are not for the
moment in the main branch but they should be released soon to be
integrated with the 2.1 release.

Best regards,
Francesco

From a71bd172bb8594ed5262cd1c12cf706f73b95e91 Mon Sep 17 00:00:00 2001
From: Francesco Abbate <francesco.bbt@gmail.com>
Date: Sun, 22 Jan 2012 19:56:46 +0100
Subject: [PATCH] Cleanup in VEGAS benchmark files

---
 vegas-bench.c   |    2 +-
 vegas-bench.lua |   38 +++++++++++++++++---------------------
 2 files changed, 18 insertions(+), 22 deletions(-)

diff --git a/vegas-bench.c b/vegas-bench.c
index 3641533..cc45ac7 100644
--- a/vegas-bench.c
+++ b/vegas-bench.c
@@ -59,5 +59,5 @@ main (void)
     display_results ("vegas final", res, err, i);
 
     gsl_monte_vegas_free (s);
-    return 1;
+    return 0;
     }
diff --git a/vegas-bench.lua b/vegas-bench.lua
index 4e08e6b..bd8588a 100644
--- a/vegas-bench.lua
+++ b/vegas-bench.lua
@@ -1,30 +1,26 @@
-
-local rng = require 'rng'
-
 local monte_vegas = dofile('vegas.lua')
+local rng = require 'rng'
+local n=9
+local lo,hi = 0,2
+local exact = n*(n+1)/2 * (hi^3 - lo^3)/3 * (hi-lo)^(n-1)
+local a,b={},{}
+for i=1,n do
+  a[i],b[i]=lo,hi
+end
+local calls = 1e6*n
+local r = rng.new('taus2')
+r:set(30776)
+local function f(x)
+  return 1*x[1]^2+2*x[2]^2+3*x[3]^2+4*x[4]^2+5*x[5]^2+6*x[6]^2+7*x[7]^2+8*x[8]^2+9*x[9]^2
+end
 
-local function testdim(n)
-  local lo,hi = 0,2
-  local exact = n*(n+1)/2 * (hi^3 - lo^3)/3 * (hi-lo)^(n-1)
-  local t={}
-  local a,b={},{}
-  for i=1,n do
-    t[i]=string.format("%s*x[%s]^2",i,i)
-    a[i],b[i]=lo,hi
-  end
-  local s=table.concat(t,"+")
-  io.write("Integrating ",s,"\nExact integral = ",exact,"\n")
-  local calls = 1e6*n
-  local r = rng.new('taus2')
-  local result,sigma,runs,cont = monte_vegas(loadstring("return function(x) return "..s.." end")(),a,b,calls,r)
-  io.write( string.format([[
+local result,sigma,runs,cont = monte_vegas(f,a,b,calls,r)
+io.write( string.format([[
 ==================
 result = %.6f
 sigma  = %.6f
 exact  = %.6f
 error  = %.6f = %.2g sigma
-i      = %d
+i      = %d 
 ]] ,result,sigma,exact, result - exact,  math.abs(result - exact)/sigma, runs))
-end
 
-    testdim(9)
-- 
1.7.5.4