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 tclconference@googlegroups.com
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 Unified Diffs Ignore Whitespace Patch

Changes to modules/math/ChangeLog.







1
2
3
4
5
6
7






2019-04-18  Arjen Markus <arjenmarkus@users.sourceforge.net>
	* misc.tcl: Add double() to calculation of mean and standard deviation in proc stats (ticket 0a030f850d4e3fc05da98aa954a6ec1b16e655d9)
	* math.test: Correct the outcome of the test for stats (consequence of ticket 0a030f850d4e3fc05da98aa954a6ec1b16e655d9)

2018-08-04  Arjen Markus <arjenmarkus@users.sourceforge.net>
	* statistics.tcl: Source stat_wasserstein.tcl and stat_logit.tcl - for new commands
	* statistics.test: Add corresponding tests
>
>
>
>
>
>







1
2
3
4
5
6
7
8
9
10
11
12
13
2019-04-23  Arjen Markus <arjenmarkus@users.sourceforge.net>
	* quasirandom.tcl: New package - generate quasi-random numbers (for instance for estimating multidimensional integrals)
	* quasirandom.test: Tests for the new package
	* quasirandom.man: Documentation for the new package
	* pkgIndex.tcl: Add the new package

2019-04-18  Arjen Markus <arjenmarkus@users.sourceforge.net>
	* misc.tcl: Add double() to calculation of mean and standard deviation in proc stats (ticket 0a030f850d4e3fc05da98aa954a6ec1b16e655d9)
	* math.test: Correct the outcome of the test for stats (consequence of ticket 0a030f850d4e3fc05da98aa954a6ec1b16e655d9)

2018-08-04  Arjen Markus <arjenmarkus@users.sourceforge.net>
	* statistics.tcl: Source stat_wasserstein.tcl and stat_logit.tcl - for new commands
	* statistics.test: Add corresponding tests

Changes to modules/math/pkgIndex.tcl.

30
31
32
33
34
35
36

package ifneeded math::decimal           1.0.3 [list source [file join $dir decimal.tcl]]
package ifneeded math::geometry          1.3.0 [list source [file join $dir geometry.tcl]]
package ifneeded math::trig              1.0   [list source [file join $dir trig.tcl]]

if {![package vsatisfies [package require Tcl] 8.6]} {return}
package ifneeded math::exact             1.0.1 [list source [file join $dir exact.tcl]]
package ifneeded math::PCA               1.0   [list source [file join $dir pca.tcl]]








>
30
31
32
33
34
35
36
37
package ifneeded math::decimal           1.0.3 [list source [file join $dir decimal.tcl]]
package ifneeded math::geometry          1.3.0 [list source [file join $dir geometry.tcl]]
package ifneeded math::trig              1.0   [list source [file join $dir trig.tcl]]

if {![package vsatisfies [package require Tcl] 8.6]} {return}
package ifneeded math::exact             1.0.1 [list source [file join $dir exact.tcl]]
package ifneeded math::PCA               1.0   [list source [file join $dir pca.tcl]]
package ifneeded math::quasirandom       1.0   [list source [file join $dir quasirandom.tcl]]

Added modules/math/quasirandom.man.











































































































































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
[vset VERSION 1]
[manpage_begin math::quasirandom n [vset VERSION]]
[keywords {quasi-random}]
[keywords mathematics]
[moddesc {Tcl Math Library}]
[titledesc {Quasi-random points for integration and Monte Carlo type methods}]
[category  Mathematics]
[require Tcl 8.6]
[require math::quasirandom [vset VERSION]]
[description]
[para]

In many applications pseudo-random numbers and pseudo-random points in a (limited)
sample space play an important role. For instance in any type of Monte Carlo simulation.
Pseudo-random numbers, however, may be too random and as a consequence a large
number of data points is required to reduce the error or fluctuation in the results
to the desired value.
[para]

Quasi-random numbers can be used as an alternative: instead of "completely" arbitrary
points, points are generated that are diverse enough to cover the entire sample space
in a more or less uniform way. As a consequence convergence to the limit can be
much faster, when such quasi-random numbers are well-chosen.
[para]

The package defines a [term class] "qrpoint" that creates a command to generate
quasi-random points in 1, 2 or more dimensions. The command can either generate
separate points, so that they can be used in a user-defined algorithm or use these
points to calculate integrals of functions defined over 1, 2 or more dimensions.
It also holds several other common algorithms. (NOTE: these are not implemented yet)
[para]
One particular characteristic of the generators is that there are no tuning parameters
involved, which makes the use particularly simple.


[section "COMMANDS"]
A quasi-random point generator is created using the [term qrpoint] class:

[list_begin definitions]

[call [cmd "::math::quasirandom::qrpoint create"] [arg NAME] [arg DIM] [opt ARGS]]
This command takes the following arguments:

[list_begin arguments]
[arg_def string NAME] The name of the command to be created (alternatively: the [term new] subcommand
will generate a unique name)
[arg_def integer/string DIM] The number of dimensions or one of: "circle", "disk", "sphere" or "ball"
[arg_def strings ARGS] Zero or more key-value pairs. The supported options are:

[list_begin itemized]
[item] [term {-start index}]: The index for the next point to be generated (default: 1)
[item] [term {-evaluations number}]: The number of evaluations to be used by default (default: 100)
[list_end]

[list_end]

[list_end]

The points that are returned lie in the hyperblock [lb]0,1[lb]^n (n the number of dimensions)
or on the unit circle, within the unit disk, on the unit sphere or within the unit ball.
[para]

Each generator supports the following subcommands:
[list_begin definitions]

[call [cmd "gen next"]]
Return the coordinates of the next quasi-random point
[nl]

[call [cmd "gen set-start"] [arg index]]
Reset the index for the next quasi-random point. This is useful to control which list of points is returned.
Returns the new or the current value, if no value is given.
[nl]

[call [cmd "gen set-evaluations"] [arg number]]
Reset the default number of evaluations in compound algorithms. Note that the actual number is the
smallest 4-fold larger or equal to the given number. (The 4-fold plays a role in the detailed integration
routine.)
[nl]

[call [cmd "gen integral"] [arg func] [arg minmax] [arg args]]
Calculate the integral of the given function over the block (or the circle, sphere etc.)

[list_begin arguments]
[arg_def string func] The name of the function to be integrated

[arg_def list minmax] List of pairs of minimum and maximum coordinates. This can be used to
map the quasi-random coordinates to the desired hyper-block.
[nl]
If the space is a circle, disk etc. then this argument should be a single value, the radius.
The circle, disk, etc. is centred at the origin. If this is not what is required, then a coordinate
transformation should be made within the function.

[arg_def strings args] Zero or more key-value pairs. The following options are supported:
[list_begin itemized]
[item] [term {-evaluations number}]: The number of evaluations to be used. If not specified use the
default of the generator object.
[list_end]

[list_end]

[list_end]

[section TODO]
Implement other algorithms and variants
[para]
Implement more unit tests.
[para]
Comparison to pseudo-random numbers for integration.


[section References]

Various algorithms exist for generating quasi-random numbers. The generators created in this package are based on:
[uri http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/]

[manpage_end]

Added modules/math/quasirandom.tcl.





























































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
# 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 provide math::quasirandom 1.0

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

        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
        }

        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
        }

        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

        foreach {key value} $args {
            switch -- $key {
            "-evaluations" {
                 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

        foreach {key value} $args {
            switch -- $key {
            "-evaluations" {
                 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}]
}

Added modules/math/quasirandom.test.



























































































































































































































































































































































































































































































































































































































































































































































































































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
# quasirandom.test --
#     Tests for the quasi-random numbers package
#
package require tcltest
namespace import ::tcltest::test

source quasirandom.tcl

#
# Functions for integration tests
#
proc const {coords} {
    return 1.0
}

proc fx {coords} {
    set x [lindex $coords 0]
    return $x
}

proc fy {coords} {
    set y [lindex $coords 1]
    return $y
}

proc fz {coords} {
    set z [lindex $coords 2]
    return $z
}

proc fxyz4 {coords} {
    lassign $coords x y z
    return [expr {($x*$y*$z)**4}]
}

#
# Auxiliary proc
#
proc equalCoords {coords1 coords2} {
    set equal 1
    foreach c1 $coords1 c2 $coords2 {
        if { $c1 != $c2 } {
            set equal 0
            break
        }
    }
    return $equal
}

#
# Create and register (in that order!) custom matching procedures
#
proc matchTolerant { expected actual } {
   set match 1
   foreach a $actual e $expected {
       if { $e != 0.0 } {
           if { abs($e-$a)>1.0e-7*abs($e) &&
                abs($e-$a)>1.0e-7*abs($a)     } {
               set match 0
               break
           }
       } else {
           if { abs($a) > 1.0e-7 } {
               set match 0
           }
       }
   }
   return $match
}
proc matchOnePercent { expected actual } {
   set match 1
   foreach a $actual e $expected {
       if { $e != 0.0 } {
           if { abs($e-$a)>1.0e-2*abs($e) &&
                abs($e-$a)>1.0e-2*abs($a)     } {
               set match 0
               break
           }
       } else {
           if { abs($a) > 1.0e-2 } {
               set match 0
           }
       }
   }
   return $match
}

::tcltest::customMatch tolerant matchTolerant
::tcltest::customMatch error1percent matchOnePercent
::tcltest::customMatch equal equalCoords


#
# Testing CoordFactors: the basis of the algorithm
# Note: exact matching
#
test "Quasirandom-0.1" "Check basic factor for 1 dimension" -body {
    set f [::math::quasirandom::CoordFactors 1]
    return [expr {1.0/$f}]
} -result 1.618033988749895

test "Quasirandom-0.2" "Check basic factor for 2 dimensions" -body {
    set f [lindex [::math::quasirandom::CoordFactors 2] 0]
    return [expr {1.0/$f}]
} -result 1.324717957244746

test "Quasirandom-0.3" "Check basic factor for 3 dimensions" -body {
    set f [lindex [::math::quasirandom::CoordFactors 3] 0]
    return [expr {1.0/$f}]
} -result 1.2207440846057596

test "Quasirandom-0.4" "Check number of factors for 10 dimensions" -body {
    return [llength [::math::quasirandom::CoordFactors 10]]
} -result 10

#
# Basic interface to the qrpoints class
#
test "Quasirandom-1.0" "Simple QR generator for two dimensions" -match tolerant -body {
    ::math::quasirandom::qrpoints create simple 2

    return [simple next]
} -result {0.7548776662466927 0.5698402909980532} -cleanup {simple destroy}

test "Quasirandom-1.1" "Simple QR generator - negative dimension" -body {
    ::math::quasirandom::qrpoints create simple -1
} -returnCodes {error} -result {The dimension argument should be a positive integer value or one of circle, disk, sphere or ball}

test "Quasirandom-1.2" "Simple QR generator - set start" -body {
    ::math::quasirandom::qrpoints create simple  2
    ::math::quasirandom::qrpoints create simple2 2 -start 2

    simple next
    set coords  [simple next]

    set coords2 [simple2 next]  ;# Should be equal to the second point for the [simple] generator

    equalCoords $coords $coords2
} -result 1 -cleanup {simple destroy; simple2 destroy}

#
# Test simple methods
#
test "Quasirandom-2.1" "set-step sets and returns the value" -match equal -body {
    ::math::quasirandom::qrpoints create simple 2

    simple set-step 100
} -result 100 -cleanup {simple destroy}

test "Quasirandom-2.2" "set-evaluations sets and returns the value" -match equal -body {
    ::math::quasirandom::qrpoints create simple 2

    simple set-evaluations 100
} -result 100 -cleanup {simple destroy}

test "Quasirandom-2.3" "set-step returns the value" -match equal -body {
    ::math::quasirandom::qrpoints create simple 2

    simple set-step 100
    simple set-step
} -result 100 -cleanup {simple destroy}

test "Quasirandom-2.4" "set-evaluations returns the value" -match equal -body {
    ::math::quasirandom::qrpoints create simple 2

    simple set-evaluations 100
    simple set-evaluations
} -result 100 -cleanup {simple destroy}

#
# Test of bounds on points
#
test "Quasirandom-3.1" "Points should fall within block" -body {
    ::math::quasirandom::qrpoints create simple 10

    set correct_bound 1

    for {set i 0} {$i < 100} {incr i} {
        set coords [simple next]

        foreach c $coords {
            if { $c < 0.0 || $c > 1.0 } {
                set correct_bound 0
                break
            }
        }
    }

    return $correct_bound
} -result 1 -cleanup {simple destroy}

test "Quasirandom-3.2" "Points should fall on a circle" -match tolerant -body {
    ::math::quasirandom::qrpoints create simple circle

    set correct_bound 1
    set radii {}

    for {set i 0} {$i < 100} {incr i} {
        set coords [simple next]

        lassign $coords x y
        lappend radii [expr {hypot($x,$y)}]
    }

    return $radii
} -result [lrepeat 100 1.0] -cleanup {simple destroy}

test "Quasirandom-3.3" "Points should fall within a disk" -match equal -body {
    ::math::quasirandom::qrpoints create simple disk

    set correct_bounds {}
    for {set i 0} {$i < 100} {incr i} {
        set coords [simple next]

        lassign $coords x y
        lappend correct_bounds [expr {hypot($x,$y) <= 1.0}]
    }

    return $correct_bounds
} -result [lrepeat 100 1] -cleanup {simple destroy}

test "Quasirandom-3.4" "Points should fall on a sphere" -match tolerant -body {
    ::math::quasirandom::qrpoints create simple sphere

    set correct_bound 1
    set radii {}

    for {set i 0} {$i < 100} {incr i} {
        set coords [simple next]

        lassign $coords x y z
        lappend radii [expr {sqrt($x**2 + $y**2 + $z**2)}]
    }

    return $radii
} -result [lrepeat 100 1.0] -cleanup {simple destroy}

test "Quasirandom-3.5" "Points should fall within a ball" -match equal -body {
    ::math::quasirandom::qrpoints create simple ball

    set correct_bounds {}
    for {set i 0} {$i < 100} {incr i} {
        set coords [simple next]

        lassign $coords x y
        lappend correct_bounds [expr {sqrt($x**2 + $y**2 + $z**2) <= 1.0}]
    }

    return $correct_bounds
} -result [lrepeat 100 1] -cleanup {simple destroy}




#
# Test of integral methods
#
# Integrating a constant function means the result is the volume
#
test "Quasirandom-4.1" "Integrate constant function - volume = 1" -match tolerant -body {
    ::math::quasirandom::qrpoints create simple 3

    set result [simple integral const {{0.0 1.0} {0.0 1.0} {0.0 1.0}}]

} -result 1.0 -cleanup {simple destroy}

test "Quasirandom-4.2" "Integrate constant function - volume = 8" -match tolerant -body {
    ::math::quasirandom::qrpoints create simple 3

    set result [simple integral const {{0.0 2.0} {0.0 2.0} {0.0 2.0}}]

} -result 8.0 -cleanup {simple destroy}

test "Quasirandom-4.3" "Integrate constant function - circle" -match tolerant -body {
    ::math::quasirandom::qrpoints create simple circle

    set result [simple integral const 2.0]

} -result [expr {2.0 * 2.0 * cos(-1.0)}] -cleanup {simple destroy}

test "Quasirandom-4.3" "Integrate constant function - disk" -match tolerant -body {
    ::math::quasirandom::qrpoints create simple disk

    set result [simple integral const 2.0]

} -result [expr {2.0**2 * cos(-1.0)}] -cleanup {simple destroy}

test "Quasirandom-4.4" "Integrate constant function - sphere" -match tolerant -body {
    ::math::quasirandom::qrpoints create simple sphere

    set result [simple integral const 2.0]

} -result [expr {4.0 * 2.0**2 * cos(-1.0)}] -cleanup {simple destroy}

test "Quasirandom-4.5" "Integrate constant function - ball" -match tolerant -body {
    ::math::quasirandom::qrpoints create simple ball

    set result [simple integral const 2.0]

} -result [expr {4.0/3.0 * 2.0**3 * cos(-1.0)}] -cleanup {simple destroy}

# We do not use too many evaluations ... error less than 1%
test "Quasirandom-4.6" "Integrate linear function (x, y, z)" -match error1percent -body {
    ::math::quasirandom::qrpoints create simple 3

    set result [list [simple integral fx {{0.0 1.0} {0.0 1.0} {0.0 1.0}}] \
                     [simple integral fy {{0.0 1.0} {0.0 1.0} {0.0 1.0}}] \
                     [simple integral fz {{0.0 1.0} {0.0 1.0} {0.0 1.0}}] ]

} -result {0.5 0.5 0.5} -cleanup {simple destroy}

#
# The function varies "sharply", so we need more evaluations
#
test "Quasirandom-4.7" "Integrate (xyz)**4" -match error1percent -body {
    ::math::quasirandom::qrpoints create simple 3

    # Exact answer is 1/125
    set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 1000]

} -result 0.0080 -cleanup {simple destroy}


#
# Detailed integration: provides error estimates but also an indication that
# the values can differ quite a bit
#
test "Quasirandom-5.1" "Integrate constant function with details - volume = 1" -match tolerant -body {
    ::math::quasirandom::qrpoints create simple 3

    set result [simple integral-detailed const {{0.0 1.0} {0.0 1.0} {0.0 1.0}}]

    set rawvalues [dict get $result -rawvalues]

} -result {1.0 1.0 1.0 1.0} -cleanup {simple destroy}


test "Quasirandom-5.2" "Integrate linear function with details - volume = 1" -match tolerant -body {
    ::math::quasirandom::qrpoints create simple 3

    set result [simple integral-detailed fx {{0.0 1.0} {0.0 1.0} {0.0 1.0}}]

    set rawvalues [dict get $result -rawvalues]

} -result {0.48924267415013695 0.48855550905424594 0.5278683439583554 0.48718117886246404} -cleanup {simple destroy}


test "Quasirandom-5.3" "Integrate (xyz)**4 with details - volume = 1" -match tolerant -body {
    ::math::quasirandom::qrpoints create simple 3

    set result [simple integral-detailed fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}}]

    set rawvalues [dict get $result -rawvalues]

} -result {0.0022115062627913935 0.009840104253511376 0.014937934937801888 0.007838969739655276} -cleanup {simple destroy}


# TODO:
# - func in different namespace
# - implement detailed integration and test the details
# - implement minimization

#
# Hm, the less than 1% error in the above test is a coincidence. The error is more
# likely to be 10%.
#
if {0} {
::math::quasirandom::qrpoints create simple 3
# Exact answer is 1/125
set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 100]
puts "fxyz4: $result"
simple set-step 0
set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 1000]
puts "fxyz4: $result"
set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 1000]
puts "fxyz4: $result"
set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 1000]
puts "fxyz4: $result"

package require math::statistics
set samples {}
for {set trial 0} {$trial < 10} {incr trial} {
    set sum 0.0

    for {set p 0} {$p < 100} {incr p} {
        set x   [expr {rand()}]
        set y   [expr {rand()}]
        set z   [expr {rand()}]
        set sum [expr {$sum + [fxyz4 [list $x $y $z]]}]
    }

    puts "Trial $trial: [expr {$sum/100.0}]"

    lappend samples [expr {$sum/100.0}]
}

puts "MonteCarlo (100):"
puts [::math::statistics::mean $samples]
puts [::math::statistics::stdev $samples]

set samples {}
for {set trial 0} {$trial < 10} {incr trial} {
    set sum 0.0

    for {set p 0} {$p < 1000} {incr p} {
        set x   [expr {rand()}]
        set y   [expr {rand()}]
        set z   [expr {rand()}]
        set sum [expr {$sum + [fxyz4 [list $x $y $z]]}]
    }

    puts "Trial $trial: [expr {$sum/1000.0}]"

    lappend samples [expr {$sum/1000.0}]
}

puts "MonteCarlo (1000):"
puts [::math::statistics::mean $samples]
puts [::math::statistics::stdev $samples]

set samples {}
for {set trial 0} {$trial < 10} {incr trial} {
    set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 100]

    lappend samples $result
}

puts "Quasi-random (100):"
puts [::math::statistics::mean $samples]
puts [::math::statistics::stdev $samples]

set samples {}
for {set trial 0} {$trial < 10} {incr trial} {
    set result [simple integral fxyz4 {{0.0 1.0} {0.0 1.0} {0.0 1.0}} -evaluations 1000]

    lappend samples $result
}

puts "Quasi-random (1000):"
puts [::math::statistics::mean $samples]
puts [::math::statistics::stdev $samples]


puts [simple integral-detailed fx {{0.0 1.0} {0.0 1.0} {0.0 1.0}}]
}