# quasirandom.tcl --
# Generate quasi-random points in n dimensions and provide simple
# methods to evaluate an integral
#
# Note: provide a OO-style interface
#
# TODO: integral-detailed, minimum, maximum
#
# Based on the blog "The Unreasonable Effectiveness of Quasirandom Sequences" by Martin Roberts,
# http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/
#
package require Tcl 8.5 9
package require TclOO
package provide math::quasirandom 1.1
namespace eval ::math::quasirandom {
# qrpoints --
# Create the class
#
::oo::class create qrpoints {
# constructor --
# Construct a new instance of the qrpoints class
#
# Arguments:
# dim Number of dimensions, or one of: circle, disk, sphere, ball
# args Zero or more key-value pairs:
# -start - start the generation with the given multiplier (integer)
# -evaluations - default number of evaluations for the integration
# (possibly others as well)
#
constructor {dimin args} {
my variable dim
my variable coord_factors
my variable step
my variable evaluations
my variable use_radius
my variable effective_dim
##nagelfar ignore
if { ( ![string is integer -strict $dimin] || $dimin <= 0 ) && $dimin ni {circle disk sphere ball} } {
return -code error "The dimension argument should be a positive integer value or one of circle, disk, sphere or ball"
}
set use_radius 1
switch -- $dimin {
"circle" {
set dim 1
set effective_dim 2
::oo::objdefine [self] {
forward next my CircleNext
forward Volume my CircleVolume
}
}
"disk" {
set dim 2
set effective_dim 2
::oo::objdefine [self] {
forward next my DiskNext
forward Volume my DiskVolume
}
}
"sphere" {
set dim 2
set effective_dim 3
::oo::objdefine [self] {
forward next my SphereNext
forward Volume my SphereVolume
}
}
"ball" {
set dim 3
set effective_dim 3
::oo::objdefine [self] {
forward next my BallNext
forward Volume my BallVolume
}
}
default {
set dim $dimin
set use_radius 0
::oo::objdefine [self] {
forward next my PlainNext
forward Volume my PlainVolume
}
}
}
set step 1
set evaluations 100
set coord_factors [::math::quasirandom::CoordFactors $dim]
foreach {key value} $args {
switch -- $key {
"-start" {
my set-step $value
}
"-evaluations" {
if { ![string is -strict integer $value] || $value <= 0 } {
return -code error "The value for the option $key should be a positive integer value"
}
my set-evaluations $value
}
default {
return -code error "Unknown option: $key -- value: $value"
}
}
}
}
# PlainNext --
# Generate the next point - for a hyperblock
#
method PlainNext {} {
my variable step
my variable coord_factors
set coords {}
foreach f $coord_factors {
lappend coords [expr {fmod( $f * $step, 1.0 )}]
}
incr step
return $coords
}
# PlainVolume --
# Calculate the volume of a hyperblock
#
# Arguments:
# minmax List of minimum and maximum per dimension
#
# Returns:
# The volume
#
method PlainVolume {minmax} {
set volume 1.0
foreach range $minmax {
lassign $range xmin xmax
set volume [expr {$volume * ($xmax-$xmin)}]
}
return $volume
}
# CircleNext --
# Generate the next point on a unit circle
#
method CircleNext {} {
set f [lindex [my PlainNext] 0]
set rad [expr {2.0 * acos(-1.0) * $f}]
set coords [list [expr {cos($rad)}] [expr {sin($rad)}]]
return $coords
}
# CircleVolume --
# Calculate the "volume" of the unit circle
#
# Arguments:
# radius Radius of the circle
#
method CircleVolume {radius} {
return [expr {$radius * 2.0*cos(-1.0)}]
}
# DiskNext --
# Generate the next point on a unit disk
#
method DiskNext {} {
while {1} {
set coords [my PlainNext]
lassign $coords x y
if { hypot($x-0.5,$y-0.5) <= 0.25 } {
set coords [list [expr {2.0*$x-1.0}] [expr {2.0*$y-1.0}]]
break
}
}
return $coords
}
# DiskVolume --
# Calculate the "volume" of the unit disk
#
# Arguments:
# radius Radius of the disk
#
method DiskVolume {radius} {
return [expr {$radius**2 * cos(-1.0)}]
}
# BallNext --
# Generate the next point on a unit ball
#
method BallNext {} {
while {1} {
set coords [my PlainNext]
lassign $coords x y z
set r [expr {($x-0.5)**2 + ($y-0.5)**2 + ($z-0.5)**2}]
if { $r <= 0.25 } {
set coords [list [expr {2.0*$x-1.0}] [expr {2.0*$y-1.0}] [expr {2.0*$z-1.0}]]
break
}
}
return $coords
}
# BallVolume --
# Calculate the volume of the unit ball
#
# Arguments:
# radius Radius of the ball
#
method BallVolume {radius} {
return [expr {4.0/3.0 * $radius**3 * cos(-1.0)}]
}
# SphereNext --
# Generate the next point on a unit sphere
#
method SphereNext {} {
set coords [my PlainNext]
lassign $coords u v
set phi [expr {2.0 * acos(-1.0) * $v}]
set lambda [expr {acos(2.0 * $u - 1.0) + 0.5 * acos(-1.0)}]
set x [expr {cos($lambda) * cos($phi)}]
set y [expr {cos($lambda) * sin($phi)}]
set z [expr {sin($lambda)}]
return [list $x $y $z]
}
# SphereVolume --
# Calculate the "volume" of the unit sphere
#
# Arguments:
# radius Radius of the sphere
#
method SphereVolume {radius} {
return [expr {4.0 * $radius**2 * cos(-1.0)}]
}
# set-step --
# Set the first step to be used
#
method set-step {{value ""}} {
my variable step
if { $value eq "" } {
return $step
}
##nagelfar ignore
if { ![string is integer -strict $value] } {
return -code error "The value for the option $key should be an integer value"
}
set step [expr {int($value)}]
}
# set-evaluations --
# Set the number of evaluations for integration
#
method set-evaluations {{value ""}} {
my variable evaluations
if { $value eq "" } {
return $evaluations
}
##nagelfar ignore
if { ![string is integer -strict $value] || $value <= 0 } {
return -code error "The value for the option $key should be a positive integer value"
}
set evaluations [expr {4*int(($value+3)/4)}] ;# Make sure it is a 4-fold
}
# integral --
# Evaluate the integral of a function over a given (rectangular) domain
#
# Arguments:
# func Function to be integrated
# minmax List of minimum and maximum bounds for each coordinate
# args Key-value pair: number of evaluations
#
# Returns:
# Estimate of the integral based on "evaluations" evaluations
# Note: no error estimate
#
method integral {func minmax args} {
my variable dim
my variable step
my variable coord_factors
my variable evaluations
my variable use_radius
my variable effective_dim
set evals $evaluations
set func [uplevel 1 [list namespace which -command $func]]
foreach {key value} $args {
switch -- $key {
"-evaluations" {
##nagelfar ignore
if { ![string is integer -strict $value] || $value <= 0 } {
return -code error "The value for the option $key should be a positive integer value"
}
set evals $value ;# Local only!
}
default {
return -code error "Unknown option: $key -- value: $value"
}
}
}
if { ! $use_radius } {
if { [llength $minmax] != $dim } {
return -code error "The number of ranges (minmax) should be equal to the dimension ($dim)"
} else {
set volume [my Volume $minmax]
}
} else {
if { ! [string is double $minmax] } {
return -code error "For a circle, disk, sphere or ball only the radius should be given"
} else {
set radius $minmax
set minmax [lrepeat $effective_dim [list 0.0 $radius]]
set volume [my Volume $radius]
}
}
set sum 0.0
for {set i 0} {$i < $evals} {incr i} {
set coords {}
foreach c [my next] range $minmax {
lassign $range xmin xmax
lappend coords [expr {$xmin + ($xmax-$xmin) * $c}]
}
set sum [expr {$sum + [$func $coords]}]
}
return [expr {$sum * $volume / $evals}]
}
# integral-detailed --
# Evaluate the integral of a function over a given (rectangular) domain
# and provide detailed information
#
# Arguments:
# func Function to be integrated
# minmax List of minimum and maximum bounds for each coordinate
# args Key-value pair: number of evaluations
#
# Returns:
# Dictionary of:
# -estimate value - estimate of the integral
# -evaluations number - total number of evaluations
# -error value - estimate of the error
# -rawvalues list - list of raw values obtained for the integral
#
method integral-detailed {func minmax args} {
my variable evaluations
set evals $evaluations
set func [uplevel 1 [list namespace which -command $func]]
foreach {key value} $args {
switch -- $key {
"-evaluations" {
##nagelfar ignore
if { ![string is integer -strict $value] || $value <= 0 } {
return -code error "The value for the option $key should be a positive integer value"
}
set evals $value ;# Local only!
}
default {
return -code error "Unknown option: $key -- value: $value"
}
}
}
lappend args -evaluations [expr {($evals+3)/4}]
for {set i 0} {$i < 4} {incr i} {
lappend rawvalues [my integral $func $minmax {*}$args]
}
set sum 0.0
set sqsum 0.0
foreach value $rawvalues {
set sum [expr {$sum + $value}]
set sqsum [expr {$sqsum + $value**2}]
}
set stdev [expr {sqrt(($sqsum - $sum**2/4.0)/3.0)}]
set sum [expr {$sum / 4.0}]
# Standard error of mean
return [dict create -estimate $sum -error [expr {$stdev/2.0}] -rawvalues $rawvalues -evaluations [expr {4*(($evals+3)/4)}]]
}
} ;# End of class
} ;# End of namespace eval
# CoordFactors --
# Determine the factors for the coordinates
#
# Arguments:
# dim Number of dimensions
#
proc ::math::quasirandom::CoordFactors {dim} {
set n [expr {$dim + 1}]
set f 1.0
for {set i 0} {$i < 10} {incr i} {
set f [expr {$f - ($f**$n-$f-1.0) / ($n*$f**($n-1)-1.0)}]
}
set factors {}
set af 1.0
for {set i 0} {$i < $dim} {incr i} {
set af [expr {$af/$f}]
lappend factors $af
}
return $factors
}
# End of code for package
# --------------------------------------------
# test --
#
if {0} {
::math::quasirandom::qrpoints create square 2
puts [square next]
puts [square next]
puts [square next]
proc f {coords} {
lassign $coords x y
expr {$x**2+$y**2}
}
proc g {coords} {
lassign $coords x y
expr {(1.0-cos($x))**2 * (1.0-cos($y))**2}
}
# Print four estimates - should not deviate too much from 10.0
puts [square integral f {{0 1} {0 3}}]
puts [square integral f {{0 1} {0 3}}]
puts [square integral f {{0 1} {0 3}}]
puts [square integral f {{0 1} {0 3}}]
# Print a sequence of estimates - should converge to (3pi/2)**2
foreach n {20 40 100 300 1000} {
square set-evaluations $n
puts "$n: [square integral g [list [list 0.0 [expr {acos(-1)}]] [list 0.0 [expr {acos(-1)}]]]]"
}
::math::quasirandom::qrpoints create block 3
puts [block next]
puts "Circle ..."
::math::quasirandom::qrpoints create circle circle
puts [circle next]
puts [circle next]
puts [circle next]
# Test values for CoordFactors
# dim = 1: 1.6180339887498948482045...
# dim = 2: 1.3247179572447460259609...
# dim = 3: 1.2207440846057594753616...
set f [::math::quasirandom::CoordFactors 1]
puts 1.6180339887498948482045...
puts [expr {1.0/$f}]
set f [lindex [::math::quasirandom::CoordFactors 2] 0]
puts 1.3247179572447460259609...
puts [expr {1.0/$f}]
set f [lindex [::math::quasirandom::CoordFactors 3] 0]
puts 1.2207440846057594753616...
puts [expr {1.0/$f}]
}
