Tcl Library Source Code

Check-in [8fdd9b5200]
Login
Bounty program for improvements to Tcl and certain Tcl packages.
Tcl 2019 Conference, Houston/TX, US, Nov 4-8
Send your abstracts to [email protected]
or submit via the online form by Sep 9.

Many hyperlinks are disabled.
Use anonymous login to enable hyperlinks.

Overview
Comment:Add a new package "quasirandom" for generating quasirandom numbers
Timelines: family | ancestors | descendants | both | trunk
Files: files | file ages | folders
SHA3-256: 8fdd9b52003df88c7ff463b1c951996e79c4648aee0ad7e4e358f4067cadb83a
User & Date: arjenmarkus 2019-04-23 19:46:34
Context
2019-04-26
12:43
coroutine properly quote coroutine name check-in: 3bea76f022 user: pooryorick tags: trunk
2019-04-23
19:46
Add a new package "quasirandom" for generating quasirandom numbers check-in: 8fdd9b5200 user: arjenmarkus tags: trunk
18:52
Solve a small problem with the math::stats proc (it did not correctly calculate the mean if the given numbers were all integers; now in the correct branch) check-in: 4c651a5ee7 user: arjenmarkus tags: trunk
Changes
Hide Diffs Side-by-Side Diffs Ignore Whitespace Patch

Changes to modules/math/ChangeLog.

            1  +2019-04-23  Arjen Markus <[email protected]>
            2  +	* quasirandom.tcl: New package - generate quasi-random numbers (for instance for estimating multidimensional integrals)
            3  +	* quasirandom.test: Tests for the new package
            4  +	* quasirandom.man: Documentation for the new package
            5  +	* pkgIndex.tcl: Add the new package
            6  +
     1      7   2019-04-18  Arjen Markus <[email protected]>
     2      8   	* misc.tcl: Add double() to calculation of mean and standard deviation in proc stats (ticket 0a030f850d4e3fc05da98aa954a6ec1b16e655d9)
     3      9   	* math.test: Correct the outcome of the test for stats (consequence of ticket 0a030f850d4e3fc05da98aa954a6ec1b16e655d9)
     4     10   
     5     11   2018-08-04  Arjen Markus <[email protected]>
     6     12   	* statistics.tcl: Source stat_wasserstein.tcl and stat_logit.tcl - for new commands
     7     13   	* statistics.test: Add corresponding tests

Changes to modules/math/pkgIndex.tcl.

    30     30   package ifneeded math::decimal           1.0.3 [list source [file join $dir decimal.tcl]]
    31     31   package ifneeded math::geometry          1.3.0 [list source [file join $dir geometry.tcl]]
    32     32   package ifneeded math::trig              1.0   [list source [file join $dir trig.tcl]]
    33     33   
    34     34   if {![package vsatisfies [package require Tcl] 8.6]} {return}
    35     35   package ifneeded math::exact             1.0.1 [list source [file join $dir exact.tcl]]
    36     36   package ifneeded math::PCA               1.0   [list source [file join $dir pca.tcl]]
           37  +package ifneeded math::quasirandom       1.0   [list source [file join $dir quasirandom.tcl]]

Added modules/math/quasirandom.man.

            1  +[vset VERSION 1]
            2  +[manpage_begin math::quasirandom n [vset VERSION]]
            3  +[keywords {quasi-random}]
            4  +[keywords mathematics]
            5  +[moddesc {Tcl Math Library}]
            6  +[titledesc {Quasi-random points for integration and Monte Carlo type methods}]
            7  +[category  Mathematics]
            8  +[require Tcl 8.6]
            9  +[require math::quasirandom [vset VERSION]]
           10  +[description]
           11  +[para]
           12  +
           13  +In many applications pseudo-random numbers and pseudo-random points in a (limited)
           14  +sample space play an important role. For instance in any type of Monte Carlo simulation.
           15  +Pseudo-random numbers, however, may be too random and as a consequence a large
           16  +number of data points is required to reduce the error or fluctuation in the results
           17  +to the desired value.
           18  +[para]
           19  +
           20  +Quasi-random numbers can be used as an alternative: instead of "completely" arbitrary
           21  +points, points are generated that are diverse enough to cover the entire sample space
           22  +in a more or less uniform way. As a consequence convergence to the limit can be
           23  +much faster, when such quasi-random numbers are well-chosen.
           24  +[para]
           25  +
           26  +The package defines a [term class] "qrpoint" that creates a command to generate
           27  +quasi-random points in 1, 2 or more dimensions. The command can either generate
           28  +separate points, so that they can be used in a user-defined algorithm or use these
           29  +points to calculate integrals of functions defined over 1, 2 or more dimensions.
           30  +It also holds several other common algorithms. (NOTE: these are not implemented yet)
           31  +[para]
           32  +One particular characteristic of the generators is that there are no tuning parameters
           33  +involved, which makes the use particularly simple.
           34  +
           35  +
           36  +[section "COMMANDS"]
           37  +A quasi-random point generator is created using the [term qrpoint] class:
           38  +
           39  +[list_begin definitions]
           40  +
           41  +[call [cmd "::math::quasirandom::qrpoint create"] [arg NAME] [arg DIM] [opt ARGS]]
           42  +This command takes the following arguments:
           43  +
           44  +[list_begin arguments]
           45  +[arg_def string NAME] The name of the command to be created (alternatively: the [term new] subcommand
           46  +will generate a unique name)
           47  +[arg_def integer/string DIM] The number of dimensions or one of: "circle", "disk", "sphere" or "ball"
           48  +[arg_def strings ARGS] Zero or more key-value pairs. The supported options are:
           49  +
           50  +[list_begin itemized]
           51  +[item] [term {-start index}]: The index for the next point to be generated (default: 1)
           52  +[item] [term {-evaluations number}]: The number of evaluations to be used by default (default: 100)
           53  +[list_end]
           54  +
           55  +[list_end]
           56  +
           57  +[list_end]
           58  +
           59  +The points that are returned lie in the hyperblock [lb]0,1[lb]^n (n the number of dimensions)
           60  +or on the unit circle, within the unit disk, on the unit sphere or within the unit ball.
           61  +[para]
           62  +
           63  +Each generator supports the following subcommands:
           64  +[list_begin definitions]
           65  +
           66  +[call [cmd "gen next"]]
           67  +Return the coordinates of the next quasi-random point
           68  +[nl]
           69  +
           70  +[call [cmd "gen set-start"] [arg index]]
           71  +Reset the index for the next quasi-random point. This is useful to control which list of points is returned.
           72  +Returns the new or the current value, if no value is given.
           73  +[nl]
           74  +
           75  +[call [cmd "gen set-evaluations"] [arg number]]
           76  +Reset the default number of evaluations in compound algorithms. Note that the actual number is the
           77  +smallest 4-fold larger or equal to the given number. (The 4-fold plays a role in the detailed integration
           78  +routine.)
           79  +[nl]
           80  +
           81  +[call [cmd "gen integral"] [arg func] [arg minmax] [arg args]]
           82  +Calculate the integral of the given function over the block (or the circle, sphere etc.)
           83  +
           84  +[list_begin arguments]
           85  +[arg_def string func] The name of the function to be integrated
           86  +
           87  +[arg_def list minmax] List of pairs of minimum and maximum coordinates. This can be used to
           88  +map the quasi-random coordinates to the desired hyper-block.
           89  +[nl]
           90  +If the space is a circle, disk etc. then this argument should be a single value, the radius.
           91  +The circle, disk, etc. is centred at the origin. If this is not what is required, then a coordinate
           92  +transformation should be made within the function.
           93  +
           94  +[arg_def strings args] Zero or more key-value pairs. The following options are supported:
           95  +[list_begin itemized]
           96  +[item] [term {-evaluations number}]: The number of evaluations to be used. If not specified use the
           97  +default of the generator object.
           98  +[list_end]
           99  +
          100  +[list_end]
          101  +
          102  +[list_end]
          103  +
          104  +[section TODO]
          105  +Implement other algorithms and variants
          106  +[para]
          107  +Implement more unit tests.
          108  +[para]
          109  +Comparison to pseudo-random numbers for integration.
          110  +
          111  +
          112  +[section References]
          113  +
          114  +Various algorithms exist for generating quasi-random numbers. The generators created in this package are based on:
          115  +[uri http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/]
          116  +
          117  +[manpage_end]

Added modules/math/quasirandom.tcl.

            1  +# quasirandom.tcl --
            2  +#     Generate quasi-random points in n dimensions and provide simple
            3  +#     methods to evaluate an integral
            4  +#
            5  +#     Note: provide a OO-style interface
            6  +#
            7  +#     TODO: integral-detailed, minimum, maximum
            8  +#
            9  +#     Based on the blog "The Unreasonable Effectiveness of Quasirandom Sequences" by Martin Roberts,
           10  +#     http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/
           11  +#
           12  +
           13  +package provide math::quasirandom 1.0
           14  +
           15  +namespace eval ::math::quasirandom {
           16  +
           17  +# qrpoints --
           18  +#     Create the class
           19  +#
           20  +::oo::class create qrpoints {
           21  +
           22  +    # constructor --
           23  +    #     Construct a new instance of the qrpoints class
           24  +    #
           25  +    # Arguments:
           26  +    #     dim             Number of dimensions, or one of: circle, disk, sphere, ball
           27  +    #     args            Zero or more key-value pairs:
           28  +    #                     -start       - start the generation with the given multiplier (integer)
           29  +    #                     -evaluations - default number of evaluations for the integration
           30  +    #                     (possibly others as well)
           31  +    #
           32  +    constructor {dimin args} {
           33  +        my variable dim
           34  +        my variable coord_factors
           35  +        my variable step
           36  +        my variable evaluations
           37  +        my variable use_radius
           38  +        my variable effective_dim
           39  +
           40  +        if { ( ![string is integer -strict $dimin] || $dimin <= 0 ) && $dimin ni {circle disk sphere ball} } {
           41  +            return -code error "The dimension argument should be a positive integer value or one of circle, disk, sphere or ball"
           42  +        }
           43  +
           44  +        set use_radius 1
           45  +        switch -- $dimin {
           46  +            "circle" {
           47  +                set dim 1
           48  +                set effective_dim 2
           49  +                ::oo::objdefine [self] {
           50  +                    forward next   my CircleNext
           51  +                    forward Volume my CircleVolume
           52  +                }
           53  +            }
           54  +            "disk" {
           55  +                set dim 2
           56  +                set effective_dim 2
           57  +                ::oo::objdefine [self] {
           58  +                    forward next   my DiskNext
           59  +                    forward Volume my DiskVolume
           60  +                }
           61  +            }
           62  +            "sphere" {
           63  +                set dim 2
           64  +                set effective_dim 3
           65  +                ::oo::objdefine [self] {
           66  +                    forward next   my SphereNext
           67  +                    forward Volume my SphereVolume
           68  +                }
           69  +            }
           70  +            "ball" {
           71  +                set dim 3
           72  +                set effective_dim 3
           73  +                ::oo::objdefine [self] {
           74  +                    forward next   my BallNext
           75  +                    forward Volume my BallVolume
           76  +                }
           77  +            }
           78  +            default {
           79  +                set dim $dimin
           80  +                set use_radius 0
           81  +                ::oo::objdefine [self] {
           82  +                    forward next   my PlainNext
           83  +                    forward Volume my PlainVolume
           84  +                }
           85  +            }
           86  +        }
           87  +
           88  +        set step        1
           89  +        set evaluations 100
           90  +
           91  +        set coord_factors [::math::quasirandom::CoordFactors $dim]
           92  +
           93  +        foreach {key value} $args {
           94  +            switch -- $key {
           95  +            "-start" {
           96  +
           97  +                 my set-step $value
           98  +            }
           99  +            "-evaluations" {
          100  +                 if { ![string is -strict integer $value] || $value <= 0 } {
          101  +                     return -code error "The value for the option $key should be a positive integer value"
          102  +                 }
          103  +
          104  +                 my set-evaluations $value
          105  +            }
          106  +            default {
          107  +                return -code error "Unknown option: $key -- value: $value"
          108  +            }
          109  +            }
          110  +        }
          111  +    }
          112  +
          113  +    # PlainNext --
          114  +    #     Generate the next point - for a hyperblock
          115  +    #
          116  +    method PlainNext {} {
          117  +        my variable step
          118  +        my variable coord_factors
          119  +
          120  +        set coords {}
          121  +        foreach f $coord_factors {
          122  +            lappend coords [expr {fmod( $f * $step, 1.0 )}]
          123  +        }
          124  +
          125  +        incr step
          126  +
          127  +        return $coords
          128  +    }
          129  +
          130  +    # PlainVolume --
          131  +    #     Calculate the volume of a hyperblock
          132  +    #
          133  +    # Arguments:
          134  +    #     minmax              List of minimum and maximum per dimension
          135  +    #
          136  +    # Returns:
          137  +    #     The volume
          138  +    #
          139  +    method PlainVolume {minmax} {
          140  +        set volume 1.0
          141  +        foreach range $minmax {
          142  +            lassign $range xmin xmax
          143  +            set volume [expr {$volume * ($xmax-$xmin)}]
          144  +        }
          145  +        return $volume
          146  +    }
          147  +
          148  +    # CircleNext --
          149  +    #     Generate the next point on a unit circle
          150  +    #
          151  +    method CircleNext {} {
          152  +
          153  +        set f      [lindex [my PlainNext] 0]
          154  +        set rad    [expr {2.0 * acos(-1.0) * $f}]
          155  +
          156  +        set coords [list [expr {cos($rad)}] [expr {sin($rad)}]]
          157  +
          158  +        return $coords
          159  +    }
          160  +
          161  +    # CircleVolume --
          162  +    #     Calculate the "volume" of the unit circle
          163  +    #
          164  +    # Arguments:
          165  +    #     radius        Radius of the circle
          166  +    #
          167  +    method CircleVolume {radius} {
          168  +         return [expr {$radius * 2.0*cos(-1.0)}]
          169  +    }
          170  +
          171  +    # DiskNext --
          172  +    #     Generate the next point on a unit disk
          173  +    #
          174  +    method DiskNext {} {
          175  +
          176  +        while {1} {
          177  +            set coords [my PlainNext]
          178  +
          179  +            lassign $coords x y
          180  +
          181  +            if { hypot($x-0.5,$y-0.5) <= 0.25 } {
          182  +                set coords [list [expr {2.0*$x-1.0}] [expr {2.0*$y-1.0}]]
          183  +                break
          184  +            }
          185  +        }
          186  +        return $coords
          187  +    }
          188  +
          189  +    # DiskVolume --
          190  +    #     Calculate the "volume" of the unit disk
          191  +    #
          192  +    # Arguments:
          193  +    #     radius        Radius of the disk
          194  +    #
          195  +    method DiskVolume {radius} {
          196  +         return [expr {$radius**2 * cos(-1.0)}]
          197  +    }
          198  +
          199  +    # BallNext --
          200  +    #     Generate the next point on a unit ball
          201  +    #
          202  +    method BallNext {} {
          203  +
          204  +        while {1} {
          205  +            set coords [my PlainNext]
          206  +
          207  +            lassign $coords x y z
          208  +
          209  +            set r [expr {($x-0.5)**2 + ($y-0.5)**2 + ($z-0.5)**2}]
          210  +            if { $r <= 0.25 } {
          211  +                set coords [list [expr {2.0*$x-1.0}] [expr {2.0*$y-1.0}] [expr {2.0*$z-1.0}]]
          212  +                break
          213  +            }
          214  +        }
          215  +
          216  +        return $coords
          217  +    }
          218  +
          219  +    # BallVolume --
          220  +    #     Calculate the volume of the unit ball
          221  +    #
          222  +    # Arguments:
          223  +    #     radius        Radius of the ball
          224  +    #
          225  +    method BallVolume {radius} {
          226  +         return [expr {4.0/3.0 * $radius**3 * cos(-1.0)}]
          227  +    }
          228  +
          229  +    # SphereNext --
          230  +    #     Generate the next point on a unit sphere
          231  +    #
          232  +    method SphereNext {} {
          233  +
          234  +        set coords [my PlainNext]
          235  +
          236  +        lassign $coords u v
          237  +
          238  +        set phi    [expr {2.0 * acos(-1.0) * $v}]
          239  +        set lambda [expr {acos(2.0 * $u - 1.0) + 0.5 * acos(-1.0)}]
          240  +
          241  +        set x      [expr {cos($lambda) * cos($phi)}]
          242  +        set y      [expr {cos($lambda) * sin($phi)}]
          243  +        set z      [expr {sin($lambda)}]
          244  +
          245  +        return [list $x $y $z]
          246  +    }
          247  +
          248  +    # SphereVolume --
          249  +    #     Calculate the "volume" of the unit sphere
          250  +    #
          251  +    # Arguments:
          252  +    #     radius        Radius of the sphere
          253  +    #
          254  +    method SphereVolume {radius} {
          255  +         return [expr {4.0 * $radius**2 * cos(-1.0)}]
          256  +    }
          257  +
          258  +    # set-step --
          259  +    #     Set the first step to be used
          260  +    #
          261  +    method set-step {{value ""}} {
          262  +        my variable step
          263  +
          264  +        if { $value eq "" } {
          265  +            return $step
          266  +        }
          267  +
          268  +        if { ![string is integer -strict $value] } {
          269  +            return -code error "The value for the option $key should be an integer value"
          270  +        }
          271  +
          272  +        set step [expr {int($value)}]
          273  +    }
          274  +
          275  +    # set-evaluations --
          276  +    #     Set the number of evaluations for integration
          277  +    #
          278  +    method set-evaluations {{value ""}} {
          279  +        my variable evaluations
          280  +
          281  +        if { $value eq "" } {
          282  +            return $evaluations
          283  +        }
          284  +
          285  +        if { ![string is integer -strict $value] || $value <= 0 } {
          286  +            return -code error "The value for the option $key should be a positive integer value"
          287  +        }
          288  +
          289  +        set evaluations [expr {4*int(($value+3)/4)}]  ;# Make sure it is a 4-fold
          290  +    }
          291  +
          292  +    # integral --
          293  +    #     Evaluate the integral of a function over a given (rectangular) domain
          294  +    #
          295  +    # Arguments:
          296  +    #     func              Function to be integrated
          297  +    #     minmax            List of minimum and maximum bounds for each coordinate
          298  +    #     args              Key-value pair: number of evaluations
          299  +    #
          300  +    # Returns:
          301  +    #     Estimate of the integral based on "evaluations" evaluations
          302  +    #     Note: no error estimate
          303  +    #
          304  +    method integral {func minmax args} {
          305  +        my variable dim
          306  +        my variable step
          307  +        my variable coord_factors
          308  +        my variable evaluations
          309  +        my variable use_radius
          310  +        my variable effective_dim
          311  +
          312  +        set evals $evaluations
          313  +
          314  +        foreach {key value} $args {
          315  +            switch -- $key {
          316  +            "-evaluations" {
          317  +                 if { ![string is integer -strict $value] || $value <= 0 } {
          318  +                     return -code error "The value for the option $key should be a positive integer value"
          319  +                 }
          320  +
          321  +                 set evals $value ;# Local only!
          322  +            }
          323  +            default {
          324  +                return -code error "Unknown option: $key -- value: $value"
          325  +            }
          326  +            }
          327  +        }
          328  +
          329  +        if { ! $use_radius } {
          330  +            if { [llength $minmax] != $dim } {
          331  +                return -code error "The number of ranges (minmax) should be equal to the dimension ($dim)"
          332  +            } else {
          333  +                set volume [my Volume $minmax]
          334  +            }
          335  +        } else {
          336  +            if { ! [string is double $minmax] } {
          337  +                return -code error "For a circle, disk, sphere or ball only the radius should be given"
          338  +            } else {
          339  +                set radius $minmax
          340  +                set minmax [lrepeat $effective_dim [list 0.0 $radius]]
          341  +                set volume [my Volume $radius]
          342  +            }
          343  +        }
          344  +
          345  +        set sum 0.0
          346  +
          347  +        for {set i 0} {$i < $evals} {incr i} {
          348  +            set coords {}
          349  +            foreach c [my next] range $minmax {
          350  +                lassign $range xmin xmax
          351  +                lappend coords [expr {$xmin + ($xmax-$xmin) * $c}]
          352  +            }
          353  +            set sum [expr {$sum + [$func $coords]}]
          354  +        }
          355  +
          356  +        return [expr {$sum * $volume / $evals}]
          357  +    }
          358  +
          359  +    # integral-detailed --
          360  +    #     Evaluate the integral of a function over a given (rectangular) domain
          361  +    #     and provide detailed information
          362  +    #
          363  +    # Arguments:
          364  +    #     func              Function to be integrated
          365  +    #     minmax            List of minimum and maximum bounds for each coordinate
          366  +    #     args              Key-value pair: number of evaluations
          367  +    #
          368  +    # Returns:
          369  +    #     Dictionary of:
          370  +    #     -estimate value     - estimate of the integral
          371  +    #     -evaluations number - total number of evaluations
          372  +    #     -error value        - estimate of the error
          373  +    #     -rawvalues list     - list of raw values obtained for the integral
          374  +    #
          375  +    method integral-detailed {func minmax args} {
          376  +        my variable evaluations
          377  +
          378  +        set evals $evaluations
          379  +
          380  +        foreach {key value} $args {
          381  +            switch -- $key {
          382  +            "-evaluations" {
          383  +                 if { ![string is integer -strict $value] || $value <= 0 } {
          384  +                     return -code error "The value for the option $key should be a positive integer value"
          385  +                 }
          386  +
          387  +                 set evals $value ;# Local only!
          388  +            }
          389  +            default {
          390  +                return -code error "Unknown option: $key -- value: $value"
          391  +            }
          392  +            }
          393  +        }
          394  +
          395  +        lappend args -evaluations [expr {($evals+3)/4}]
          396  +
          397  +        for {set i 0} {$i < 4} {incr i} {
          398  +            lappend rawvalues [my integral $func $minmax {*}$args]
          399  +        }
          400  +
          401  +        set sum   0.0
          402  +        set sqsum 0.0
          403  +
          404  +        foreach value $rawvalues {
          405  +            set sum   [expr {$sum + $value}]
          406  +            set sqsum [expr {$sqsum + $value**2}]
          407  +        }
          408  +
          409  +        set stdev [expr {sqrt(($sqsum - $sum**2/4.0)/3.0)}]
          410  +        set sum   [expr {$sum / 4.0}]
          411  +                                            # Standard error of mean
          412  +        return [dict create -estimate $sum -error [expr {$stdev/2.0}] -rawvalues $rawvalues -evaluations [expr {4*(($evals+3)/4)}]]
          413  +    }
          414  +
          415  +} ;# End of class
          416  +
          417  +} ;# End of namespace eval
          418  +
          419  +# CoordFactors --
          420  +#     Determine the factors for the coordinates
          421  +#
          422  +# Arguments:
          423  +#     dim         Number of dimensions
          424  +#
          425  +proc ::math::quasirandom::CoordFactors {dim} {
          426  +    set n [expr {$dim + 1}]
          427  +
          428  +    set f 1.0
          429  +    for {set i 0} {$i < 10} {incr i} {
          430  +        set f [expr {$f - ($f**$n-$f-1.0) / ($n*$f**($n-1)-1.0)}]
          431  +    }
          432  +
          433  +    set factors {}
          434  +    set af      1.0
          435  +
          436  +    for {set i 0} {$i < $dim} {incr i} {
          437  +        set af [expr {$af/$f}]
          438  +        lappend factors $af
          439  +    }
          440  +
          441  +    return $factors
          442  +}
          443  +
          444  +# End of code for package
          445  +
          446  +# --------------------------------------------
          447  +# test --
          448  +#
          449  +
          450  +if {0} {
          451  +
          452  +::math::quasirandom::qrpoints create square 2
          453  +
          454  +puts [square next]
          455  +puts [square next]
          456  +puts [square next]
          457  +
          458  +
          459  +proc f {coords} {
          460  +    lassign $coords x y
          461  +
          462  +    expr {$x**2+$y**2}
          463  +}
          464  +
          465  +proc g {coords} {
          466  +    lassign $coords x y
          467  +
          468  +    expr {(1.0-cos($x))**2 * (1.0-cos($y))**2}
          469  +}
          470  +
          471  +# Print four estimates - should not deviate too much from 10.0
          472  +puts [square integral f {{0 1} {0 3}}]
          473  +puts [square integral f {{0 1} {0 3}}]
          474  +puts [square integral f {{0 1} {0 3}}]
          475  +puts [square integral f {{0 1} {0 3}}]
          476  +
          477  +# Print a sequence of estimates - should converge to (3pi/2)**2
          478  +foreach n {20 40 100 300 1000} {
          479  +    square set-evaluations $n
          480  +
          481  +    puts "$n: [square integral g [list [list 0.0 [expr {acos(-1)}]] [list 0.0 [expr {acos(-1)}]]]]"
          482  +}
          483  +
          484  +
          485  +::math::quasirandom::qrpoints create block 3
          486  +puts [block next]
          487  +
          488  +puts "Circle ..."
          489  +::math::quasirandom::qrpoints create circle circle
          490  +puts [circle next]
          491  +puts [circle next]
          492  +puts [circle next]
          493  +
          494  +# Test values for CoordFactors
          495  +# dim = 1: 1.6180339887498948482045...
          496  +# dim = 2: 1.3247179572447460259609...
          497  +# dim = 3: 1.2207440846057594753616...
          498  +
          499  +set f [::math::quasirandom::CoordFactors 1]
          500  +puts 1.6180339887498948482045...
          501  +puts [expr {1.0/$f}]
          502  +
          503  +set f [lindex [::math::quasirandom::CoordFactors 2] 0]
          504  +puts 1.3247179572447460259609...
          505  +puts [expr {1.0/$f}]
          506  +
          507  +set f [lindex [::math::quasirandom::CoordFactors 3] 0]
          508  +puts 1.2207440846057594753616...
          509  +puts [expr {1.0/$f}]
          510  +}

Added modules/math/quasirandom.test.

            1  +# quasirandom.test --
            2  +#     Tests for the quasi-random numbers package
            3  +#
            4  +package require tcltest
            5  +namespace import ::tcltest::test
            6  +
            7  +source quasirandom.tcl
            8  +
            9  +#
           10  +# Functions for integration tests
           11  +#
           12  +proc const {coords} {
           13  +    return 1.0
           14  +}
           15  +
           16  +proc fx {coords} {
           17  +    set x [lindex $coords 0]
           18  +    return $x
           19  +}
           20  +
           21  +proc fy {coords} {
           22  +    set y [lindex $coords 1]
           23  +    return $y
           24  +}
           25  +
           26  +proc fz {coords} {
           27  +    set z [lindex $coords 2]
           28  +    return $z
           29  +}
           30  +
           31  +proc fxyz4 {coords} {
           32  +    lassign $coords x y z
           33  +    return [expr {($x*$y*$z)**4}]
           34  +}
           35  +
           36  +#
           37  +# Auxiliary proc
           38  +#
           39  +proc equalCoords {coords1 coords2} {
           40  +    set equal 1
           41  +    foreach c1 $coords1 c2 $coords2 {
           42  +        if { $c1 != $c2 } {
           43  +            set equal 0
           44  +            break
           45  +        }
           46  +    }
           47  +    return $equal
           48  +}
           49  +
           50  +#
           51  +# Create and register (in that order!) custom matching procedures
           52  +#
           53  +proc matchTolerant { expected actual } {
           54  +   set match 1
           55  +   foreach a $actual e $expected {
           56  +       if { $e != 0.0 } {
           57  +           if { abs($e-$a)>1.0e-7*abs($e) &&
           58  +                abs($e-$a)>1.0e-7*abs($a)     } {
           59  +               set match 0
           60  +               break
           61  +           }
           62  +       } else {
           63  +           if { abs($a) > 1.0e-7 } {
           64  +               set match 0
           65  +           }
           66  +       }
           67  +   }
           68  +   return $match
           69  +}
           70  +proc matchOnePercent { expected actual } {
           71  +   set match 1
           72  +   foreach a $actual e $expected {
           73  +       if { $e != 0.0 } {
           74  +           if { abs($e-$a)>1.0e-2*abs($e) &&
           75  +                abs($e-$a)>1.0e-2*abs($a)     } {
           76  +               set match 0
           77  +               break
           78  +           }
           79  +       } else {
           80  +           if { abs($a) > 1.0e-2 } {
           81  +               set match 0
           82  +           }
           83  +       }
           84  +   }
           85  +   return $match
           86  +}
           87  +
           88  +::tcltest::customMatch tolerant matchTolerant
           89  +::tcltest::customMatch error1percent matchOnePercent
           90  +::tcltest::customMatch equal equalCoords
           91  +
           92  +
           93  +#
           94  +# Testing CoordFactors: the basis of the algorithm
           95  +# Note: exact matching
           96  +#
           97  +test "Quasirandom-0.1" "Check basic factor for 1 dimension" -body {
           98  +    set f [::math::quasirandom::CoordFactors 1]
           99  +    return [expr {1.0/$f}]
          100  +} -result 1.618033988749895
          101  +
          102  +test "Quasirandom-0.2" "Check basic factor for 2 dimensions" -body {
          103  +    set f [lindex [::math::quasirandom::CoordFactors 2] 0]
          104  +    return [expr {1.0/$f}]
          105  +} -result 1.324717957244746
          106  +
          107  +test "Quasirandom-0.3" "Check basic factor for 3 dimensions" -body {
          108  +    set f [lindex [::math::quasirandom::CoordFactors 3] 0]
          109  +    return [expr {1.0/$f}]
          110  +} -result 1.2207440846057596
          111  +
          112  +test "Quasirandom-0.4" "Check number of factors for 10 dimensions" -body {
          113  +    return [llength [::math::quasirandom::CoordFactors 10]]
          114  +} -result 10
          115  +
          116  +#
          117  +# Basic interface to the qrpoints class
          118  +#
          119  +test "Quasirandom-1.0" "Simple QR generator for two dimensions" -match tolerant -body {
          120  +    ::math::quasirandom::qrpoints create simple 2
          121  +
          122  +    return [simple next]
          123  +} -result {0.7548776662466927 0.5698402909980532} -cleanup {simple destroy}
          124  +
          125  +test "Quasirandom-1.1" "Simple QR generator - negative dimension" -body {
          126  +    ::math::quasirandom::qrpoints create simple -1
          127  +} -returnCodes {error} -result {The dimension argument should be a positive integer value or one of circle, disk, sphere or ball}
          128  +
          129  +test "Quasirandom-1.2" "Simple QR generator - set start" -body {
          130  +    ::math::quasirandom::qrpoints create simple  2
          131  +    ::math::quasirandom::qrpoints create simple2 2 -start 2
          132  +
          133  +    simple next
          134  +    set coords  [simple next]
          135  +
          136  +    set coords2 [simple2 next]  ;# Should be equal to the second point for the [simple] generator
          137  +
          138  +    equalCoords $coords $coords2
          139  +} -result 1 -cleanup {simple destroy; simple2 destroy}
          140  +
          141  +#
          142  +# Test simple methods
          143  +#
          144  +test "Quasirandom-2.1" "set-step sets and returns the value" -match equal -body {
          145  +    ::math::quasirandom::qrpoints create simple 2
          146  +
          147  +    simple set-step 100
          148  +} -result 100 -cleanup {simple destroy}
          149  +
          150  +test "Quasirandom-2.2" "set-evaluations sets and returns the value" -match equal -body {
          151  +    ::math::quasirandom::qrpoints create simple 2
          152  +
          153  +    simple set-evaluations 100
          154  +} -result 100 -cleanup {simple destroy}
          155  +
          156  +test "Quasirandom-2.3" "set-step returns the value" -match equal -body {
          157  +    ::math::quasirandom::qrpoints create simple 2
          158  +
          159  +    simple set-step 100
          160  +    simple set-step
          161  +} -result 100 -cleanup {simple destroy}
          162  +
          163  +test "Quasirandom-2.4" "set-evaluations returns the value" -match equal -body {
          164  +    ::math::quasirandom::qrpoints create simple 2
          165  +
          166  +    simple set-evaluations 100
          167  +    simple set-evaluations
          168  +} -result 100 -cleanup {simple destroy}
          169  +
          170  +#
          171  +# Test of bounds on points
          172  +#
          173  +test "Quasirandom-3.1" "Points should fall within block" -body {
          174  +    ::math::quasirandom::qrpoints create simple 10
          175  +
          176  +    set correct_bound 1
          177  +
          178  +    for {set i 0} {$i < 100} {incr i} {
          179  +        set coords [simple next]
          180  +
          181  +        foreach c $coords {
          182  +            if { $c < 0.0 || $c > 1.0 } {
          183  +                set correct_bound 0
          184  +                break
          185  +            }
          186  +        }
          187  +    }
          188  +
          189  +    return $correct_bound
          190  +} -result 1 -cleanup {simple destroy}
          191  +
          192  +test "Quasirandom-3.2" "Points should fall on a circle" -match tolerant -body {
          193  +    ::math::quasirandom::qrpoints create simple circle
          194  +
          195  +    set correct_bound 1
          196  +    set radii {}
          197  +
          198  +    for {set i 0} {$i < 100} {incr i} {
          199  +        set coords [simple next]
          200  +
          201  +        lassign $coords x y
          202  +        lappend radii [expr {hypot($x,$y)}]
          203  +    }
          204  +
          205  +    return $radii
          206  +} -result [lrepeat 100 1.0] -cleanup {simple destroy}
          207  +
          208  +test "Quasirandom-3.3" "Points should fall within a disk" -match equal -body {
          209  +    ::math::quasirandom::qrpoints create simple disk
          210  +
          211  +    set correct_bounds {}
          212  +    for {set i 0} {$i < 100} {incr i} {
          213  +        set coords [simple next]
          214  +
          215  +        lassign $coords x y
          216  +        lappend correct_bounds [expr {hypot($x,$y) <= 1.0}]
          217  +    }
          218  +
          219  +    return $correct_bounds
          220  +} -result [lrepeat 100 1] -cleanup {simple destroy}
          221  +
          222  +test "Quasirandom-3.4" "Points should fall on a sphere" -match tolerant -body {
          223  +    ::math::quasirandom::qrpoints create simple sphere
          224  +
          225  +    set correct_bound 1
          226  +    set radii {}
          227  +
          228  +    for {set i 0} {$i < 100} {incr i} {
          229  +        set coords [simple next]
          230  +
          231  +        lassign $coords x y z
          232  +        lappend radii [expr {sqrt($x**2 + $y**2 + $z**2)}]
          233  +    }
          234  +
          235  +    return $radii
          236  +} -result [lrepeat 100 1.0] -cleanup {simple destroy}
          237  +
          238  +test "Quasirandom-3.5" "Points should fall within a ball" -match equal -body {
          239  +    ::math::quasirandom::qrpoints create simple ball
          240  +
          241  +    set correct_bounds {}
          242  +    for {set i 0} {$i < 100} {incr i} {
          243  +        set coords [simple next]
          244  +
          245  +        lassign $coords x y
          246  +        lappend correct_bounds [expr {sqrt($x**2 + $y**2 + $z**2) <= 1.0}]
          247  +    }
          248  +
          249  +    return $correct_bounds
          250  +} -result [lrepeat 100 1] -cleanup {simple destroy}
          251  +
          252  +
          253  +
          254  +
          255  +#
          256  +# Test of integral methods
          257  +#
          258  +# Integrating a constant function means the result is the volume
          259  +#
          260  +test "Quasirandom-4.1" "Integrate constant function - volume = 1" -match tolerant -body {
          261  +    ::math::quasirandom::qrpoints create simple 3
          262  +
          263  +    set result [simple integral const {{0.0 1.0} {0.0 1.0} {0.0 1.0}}]
          264  +
          265  +} -result 1.0 -cleanup {simple destroy}
          266  +
          267  +test "Quasirandom-4.2" "Integrate constant function - volume = 8" -match tolerant -body {
          268  +    ::math::quasirandom::qrpoints create simple 3
          269  +
          270  +    set result [simple integral const {{0.0 2.0} {0.0 2.0} {0.0 2.0}}]
          271  +
          272  +} -result 8.0 -cleanup {simple destroy}
          273  +
          274  +test "Quasirandom-4.3" "Integrate constant function - circle" -match tolerant -body {
          275  +    ::math::quasirandom::qrpoints create simple circle
          276  +
          277  +    set result [simple integral const 2.0]
          278  +
          279  +} -result [expr {2.0 * 2.0 * cos(-1.0)}] -cleanup {simple destroy}
          280  +
          281  +test "Quasirandom-4.3" "Integrate constant function - disk" -match tolerant -body {
          282  +    ::math::quasirandom::qrpoints create simple disk
          283  +
          284  +    set result [simple integral const 2.0]
          285  +
          286  +} -result [expr {2.0**2 * cos(-1.0)}] -cleanup {simple destroy}
          287  +
          288  +test "Quasirandom-4.4" "Integrate constant function - sphere" -match tolerant -body {
          289  +    ::math::quasirandom::qrpoints create simple sphere
          290  +
          291  +    set result [simple integral const 2.0]
          292  +
          293  +} -result [expr {4.0 * 2.0**2 * cos(-1.0)}] -cleanup {simple destroy}
          294  +
          295  +test "Quasirandom-4.5" "Integrate constant function - ball" -match tolerant -body {
          296  +    ::math::quasirandom::qrpoints create simple ball
          297  +
          298  +    set result [simple integral const 2.0]
          299  +
          300  +} -result [expr {4.0/3.0 * 2.0**3 * cos(-1.0)}] -cleanup {simple destroy}
          301  +
          302  +# We do not use too many evaluations ... error less than 1%
          303  +test "Quasirandom-4.6" "Integrate linear function (x, y, z)" -match error1percent -body {
          304  +    ::math::quasirandom::qrpoints create simple 3
          305  +
          306  +    set result [list [simple integral fx {{0.0 1.0} {0.0 1.0} {0.0 1.0}}] \
          307  +                     [simple integral fy {{0.0 1.0} {0.0 1.0} {0.0 1.0}}] \
          308  +                     [simple integral fz {{0.0 1.0} {0.0 1.0} {0.0 1.0}}] ]
          309  +
          310  +} -result {0.5 0.5 0.5} -cleanup {simple destroy}
          311  +
          312  +#
          313  +# The function varies "sharply", so we need more evaluations
          314  +#
          315  +test "Quasirandom-4.7" "Integrate (xyz)**4" -match error1percent -body {
          316  +    ::math::quasirandom::qrpoints create simple 3
          317  +
          318  +    # Exact answer is 1/125
          319  +    set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 1000]
          320  +
          321  +} -result 0.0080 -cleanup {simple destroy}
          322  +
          323  +
          324  +#
          325  +# Detailed integration: provides error estimates but also an indication that
          326  +# the values can differ quite a bit
          327  +#
          328  +test "Quasirandom-5.1" "Integrate constant function with details - volume = 1" -match tolerant -body {
          329  +    ::math::quasirandom::qrpoints create simple 3
          330  +
          331  +    set result [simple integral-detailed const {{0.0 1.0} {0.0 1.0} {0.0 1.0}}]
          332  +
          333  +    set rawvalues [dict get $result -rawvalues]
          334  +
          335  +} -result {1.0 1.0 1.0 1.0} -cleanup {simple destroy}
          336  +
          337  +
          338  +test "Quasirandom-5.2" "Integrate linear function with details - volume = 1" -match tolerant -body {
          339  +    ::math::quasirandom::qrpoints create simple 3
          340  +
          341  +    set result [simple integral-detailed fx {{0.0 1.0} {0.0 1.0} {0.0 1.0}}]
          342  +
          343  +    set rawvalues [dict get $result -rawvalues]
          344  +
          345  +} -result {0.48924267415013695 0.48855550905424594 0.5278683439583554 0.48718117886246404} -cleanup {simple destroy}
          346  +
          347  +
          348  +test "Quasirandom-5.3" "Integrate (xyz)**4 with details - volume = 1" -match tolerant -body {
          349  +    ::math::quasirandom::qrpoints create simple 3
          350  +
          351  +    set result [simple integral-detailed fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}}]
          352  +
          353  +    set rawvalues [dict get $result -rawvalues]
          354  +
          355  +} -result {0.0022115062627913935 0.009840104253511376 0.014937934937801888 0.007838969739655276} -cleanup {simple destroy}
          356  +
          357  +
          358  +# TODO:
          359  +# - func in different namespace
          360  +# - implement detailed integration and test the details
          361  +# - implement minimization
          362  +
          363  +#
          364  +# Hm, the less than 1% error in the above test is a coincidence. The error is more
          365  +# likely to be 10%.
          366  +#
          367  +if {0} {
          368  +::math::quasirandom::qrpoints create simple 3
          369  +# Exact answer is 1/125
          370  +set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 100]
          371  +puts "fxyz4: $result"
          372  +simple set-step 0
          373  +set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 1000]
          374  +puts "fxyz4: $result"
          375  +set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 1000]
          376  +puts "fxyz4: $result"
          377  +set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 1000]
          378  +puts "fxyz4: $result"
          379  +
          380  +package require math::statistics
          381  +set samples {}
          382  +for {set trial 0} {$trial < 10} {incr trial} {
          383  +    set sum 0.0
          384  +
          385  +    for {set p 0} {$p < 100} {incr p} {
          386  +        set x   [expr {rand()}]
          387  +        set y   [expr {rand()}]
          388  +        set z   [expr {rand()}]
          389  +        set sum [expr {$sum + [fxyz4 [list $x $y $z]]}]
          390  +    }
          391  +
          392  +    puts "Trial $trial: [expr {$sum/100.0}]"
          393  +
          394  +    lappend samples [expr {$sum/100.0}]
          395  +}
          396  +
          397  +puts "MonteCarlo (100):"
          398  +puts [::math::statistics::mean $samples]
          399  +puts [::math::statistics::stdev $samples]
          400  +
          401  +set samples {}
          402  +for {set trial 0} {$trial < 10} {incr trial} {
          403  +    set sum 0.0
          404  +
          405  +    for {set p 0} {$p < 1000} {incr p} {
          406  +        set x   [expr {rand()}]
          407  +        set y   [expr {rand()}]
          408  +        set z   [expr {rand()}]
          409  +        set sum [expr {$sum + [fxyz4 [list $x $y $z]]}]
          410  +    }
          411  +
          412  +    puts "Trial $trial: [expr {$sum/1000.0}]"
          413  +
          414  +    lappend samples [expr {$sum/1000.0}]
          415  +}
          416  +
          417  +puts "MonteCarlo (1000):"
          418  +puts [::math::statistics::mean $samples]
          419  +puts [::math::statistics::stdev $samples]
          420  +
          421  +set samples {}
          422  +for {set trial 0} {$trial < 10} {incr trial} {
          423  +    set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 100]
          424  +
          425  +    lappend samples $result
          426  +}
          427  +
          428  +puts "Quasi-random (100):"
          429  +puts [::math::statistics::mean $samples]
          430  +puts [::math::statistics::stdev $samples]
          431  +
          432  +set samples {}
          433  +for {set trial 0} {$trial < 10} {incr trial} {
          434  +    set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 1000]
          435  +
          436  +    lappend samples $result
          437  +}
          438  +
          439  +puts "Quasi-random (1000):"
          440  +puts [::math::statistics::mean $samples]
          441  +puts [::math::statistics::stdev $samples]
          442  +
          443  +
          444  +puts [simple integral-detailed fx {{0.0 1.0} {0.0 1.0} {0.0 1.0}}]
          445  +}