1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
open Gsl_math
let _ =
Gsl_error.init ()
let exact =
let e = Gsl_sf.gamma(0.25) ** 4. /. (4. *. pi ** 3.) in
Printf.printf "computing exact: %.9f\n" e ;
e
let g =
let a = 1. /. (pi *. pi *. pi) in
fun x ->
a /. (1. -. cos x.(0) *. cos x.(1) *. cos x.(2))
let display_results title
{ Gsl_fun.res=result; Gsl_fun.err=error } =
Printf.printf "%s ==================\n" title ;
Printf.printf "result = % .6f\n" result ;
Printf.printf "sigma = % .6f\n" error ;
Printf.printf "exact = % .6f\n" exact ;
Printf.printf "error = % .6f = %.1g sigma\n"
(result -. exact)
((abs_float (result -. exact)) /. error)
let compute rng =
let lo = [| 0.; 0.; 0.; |] in
let up = [| pi; pi; pi |] in
(* let gslfun = Gsl_monte.register g 3 in*)
let gslfun = g in
let calls = 500000 in
begin
(* let s = Gsl_monte.alloc_plain 3 in*)
(* let res = Gsl_monte.integrate_plain gslfun ~lo ~up calls rng s in*)
(* Gsl_monte.free_plain s ;*)
let res = Gsl_monte.integrate
Gsl_monte.PLAIN gslfun ~lo ~up calls rng in
display_results "PLAIN" res ;
print_newline ()
end ;
begin
(* let s = Gsl_monte.alloc_miser 3 in*)
(* let res = Gsl_monte.integrate_miser gslfun ~lo ~up calls rng s in*)
(* Gsl_monte.free_miser s ;*)
let res = Gsl_monte.integrate
Gsl_monte.MISER gslfun ~lo ~up calls rng in
display_results "MISER" res ;
print_newline ()
end ;
begin
let state = Gsl_monte.alloc_vegas 3 in
let params = Gsl_monte.get_vegas_params state in
let oc = open_out "truc" in
Gsl_monte.set_vegas_params state { params with
Gsl_monte.verbose = 0 ;
Gsl_monte.ostream = Some oc } ;
let res = Gsl_monte.integrate_vegas gslfun ~lo ~up 10000 rng state in
display_results "VEGAS warm-up" res ;
Printf.printf "converging...\n" ; flush stdout ;
let rec proc () =
let { Gsl_fun.res=result; Gsl_fun.err=err } as res =
Gsl_monte.integrate_vegas gslfun ~lo ~up (calls / 5) rng state in
let { Gsl_monte.chisq = chisq } = Gsl_monte.get_vegas_info state in
Printf.printf "result = % .6f sigma = % .6f chisq/dof = %.1f\n"
result err chisq ;
flush stdout ;
if (abs_float (chisq -. 1.)) > 0.5
then proc ()
else res
in
let res_final = proc () in
display_results "VEGAS final" res_final ;
Gsl_monte.free_vegas state ;
close_out oc
end
let _ =
Gsl_rng.env_setup ();
let rng = Gsl_rng.make (Gsl_rng.default ()) in
Printf.printf "using %s RNG\n" (Gsl_rng.name rng) ;
print_newline () ;
compute rng
|
