Go deh!

Mainly Tech projects on Python and Electronic Design Automation.

Sunday, June 21, 2020

Successive Approximation, Hashes in high dimensions; Geohashing.


Successive Approximation

Given a named variable and the range in which its value must exist, generates successive binary digits which can be used to converge on the value successive binary digits approximating the value to ever greater precision.

Binary search for a value.

Starting with the original lower and upper range boundaries for the named variable, the value is compared to the range midpoint: if its lower than the midpoint then the range maximum is updated to become the midpoint as the value is in this half; otherwise the range minimum is updated with the midpoint to narrow therange containing the value by a half on each round. If the lower half of the range is used then generate a zero otherwise generate a one for this "round" in the binary search.

[1]

There is a TUTOR mode to turn on and before giving the encoder a whirl.

[2]
Successively encoding 0.6666666666666666 between [-2, 4]: is <= 1.0 is > -0.5 is > 0.25 is > 0.625 is <= 0.8125 is <= 0.71875 is <= 0.671875 is > 0.6484375 is > 0.66015625 Resulting bit stream: [0, 1, 1, 1, 0, 0, 0, 1, 1]. 
([0, 1, 1, 1, 0, 0, 0, 1, 1], 0.66015625, 0.671875)

You should be able to see how the encoder successively approaches the value.

Lets try some corner cases:

[3]
Successively encoding 0.5 between [0, 1]: is <= 0.5 is > 0.25 is > 0.375 is > 0.4375 is > 0.46875 is > 0.484375 is > 0.4921875 is > 0.49609375 Resulting bit stream: [0, 1, 1, 1, 1, 1, 1, 1]. 
([0, 1, 1, 1, 1, 1, 1, 1], 0.49609375, 0.5)
[4]
Successively encoding 0 between [0, 1]: is <= 0.5 is <= 0.25 is <= 0.125 is <= 0.0625 is <= 0.03125 is <= 0.015625 is <= 0.0078125 is <= 0.00390625 Resulting bit stream: [0, 0, 0, 0, 0, 0, 0, 0]. 
([0, 0, 0, 0, 0, 0, 0, 0], 0, 0.00390625)
[5]
Successively encoding 1 between [0, 1]: is > 0.5 is > 0.75 is > 0.875 is > 0.9375 is > 0.96875 is > 0.984375 is > 0.9921875 is > 0.99609375 Resulting bit stream: [1, 1, 1, 1, 1, 1, 1, 1]. 
([1, 1, 1, 1, 1, 1, 1, 1], 0.99609375, 1)

Bitstream decoding

With the encoders bitstream and theinitial range, the initial value can be recovered to ever greater accuracy as more bits are processed. The decoder takes the same decisions as the encoder when successively bisecting the initial range.

[6]

Lets encode a value then decode the whole of its bitstream

[7]
Successively encoding 0.6666666666666666 between [-2, 4]: is <= 1.0 is > -0.5 is > 0.25 is > 0.625 is <= 0.8125 is <= 0.71875 is <= 0.671875 is > 0.6484375 is > 0.66015625 is > 0.666015625 is <= 0.6689453125 is <= 0.66748046875 is <= 0.666748046875 is > 0.6663818359375 is > 0.66656494140625 Resulting bit stream: [0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1]. Successively decoding [0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1] between [-2, 4]: in range [-2, 1.0] in range [-0.5, 1.0] in range [0.25, 1.0] in range [0.625, 1.0] in range [0.625, 0.8125] in range [0.625, 0.71875] in range [0.625, 0.671875] in range [0.6484375, 0.671875] in range [0.66015625, 0.671875] in range [0.666015625, 0.671875] in range [0.666015625, 0.6689453125] in range [0.666015625, 0.66748046875] in range [0.666015625, 0.666748046875] in range [0.6663818359375, 0.666748046875] in range [0.66656494140625, 0.666748046875] Resulting value: 0.666656494140625 ±9.1552734375e-05.

Truncating binary digits

The full 15 bits of the encoding above give a result that is correct when rounded to 5 decimal digits. If a smaller number of bits (from the left), are decoded then the value should be reproduced to a lesser level of precision. Lets try:

[8]
# Decode the first 10 bits. Successively decoding [0, 1, 1, 1, 0, 0, 0, 1, 1, 1] between [-2, 4]: in range [-2, 1.0] in range [-0.5, 1.0] in range [0.25, 1.0] in range [0.625, 1.0] in range [0.625, 0.8125] in range [0.625, 0.71875] in range [0.625, 0.671875] in range [0.6484375, 0.671875] in range [0.66015625, 0.671875] in range [0.666015625, 0.671875] Resulting value: 0.6689453125 ±0.0029296875. # Decode only the first 5 bits. Successively decoding [0, 1, 1, 1, 0] between [-2, 4]: in range [-2, 1.0] in range [-0.5, 1.0] in range [0.25, 1.0] in range [0.625, 1.0] in range [0.625, 0.8125] Resulting value: 0.71875 ±0.09375.

Higher dimensions

What if we have more than one named variables values to encode in a similar fashion - preserving the ability to truncate on the right, (or generate more bits), to affect precision. Well lets give each variable, which we can think of as a "dimension" a name, a range, and all the dimensions a constant ordering.

If we assume the same precision for each dimension and interleave the bitstream for the dimensions in encoding, then the combined interleaved bitstream could be truncated as before ; de-interleaved, then individually decoded back to generate values for each dimension.

[9]
[10]
Encode then decode (-1.5, 2.7) -1.5001220703125 ±0.00012207 2.700042724609375 ±0.000152588

Again, we can truncate and take the leftmost 14 bits of the encoding to decode to less precision in both values

[11]
-1.53125 ±0.03125) 2.6953125 ±0.0390625)

A closer look at interleaving.

Create separate bitstreams for two dimensions and compare to that of the combined.

[12]
Dimensions: (Dim(name='s', rangemin=-4, rangemax=4), Dim(name='t', rangemin=-5, rangemax=5)) Values: (-1.5, 2.7)
[13]
Successively encoding -1.5 between [-4, 4]: is <= 0.0 is > -2.0 is <= -1.0 is <= -1.5 is > -1.75 is > -1.625 is > -1.5625 Resulting bit stream: [0, 1, 0, 0, 1, 1, 1].
[14]
Successively encoding 2.7 between [-5, 5]: is > 0.0 is > 2.5 is <= 3.75 is <= 3.125 is <= 2.8125 is > 2.65625 is <= 2.734375 Resulting bit stream: [1, 1, 0, 0, 0, 1, 0].
[15]
[16]
[0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0]
[17]
[0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0]
[18]
Successively encoding -1.5 between [-4, 4]: is <= 0.0 is > -2.0 is <= -1.0 is <= -1.5 is > -1.75 is > -1.625 is > -1.5625 Resulting bit stream: [0, 1, 0, 0, 1, 1, 1]. Successively encoding 2.7 between [-5, 5]: is > 0.0 is > 2.5 is <= 3.75 is <= 3.125 is <= 2.8125 is > 2.65625 is <= 2.734375 Resulting bit stream: [1, 1, 0, 0, 0, 1, 0]. OK

On precision.

The code works to the precision of the Python datatypes passed into it. The ranges so far have used Pythons 64 bit floating point datatype for range and value. If decimals or fractions are used then much more precision could be achieved with the same code due to duck typing.

Human readability

[19]
[0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0]

These strings of binary bits are hard for us to accurately remember and transcribe. Traditionally we have used octal, base-8; and hexadecimal, base 16 number bases to get shorter, more memorable, representations for binary numbers.

I am going to use abase 32 numbering scheme to encode the bits. Every five successive bits will become one of 32 different characters. The characters used are the digits zero to nine then 22 lowercase ASCII letters missing out certain letters, (a, i, l, and o), that could be mistaken for numbers.

Users will more naturally express precision in terms of the number of characters in the representation and so the underlying bit precision will be in multiples of 5-bits.

[20]
[21]
[22]
Dimensions: (Dim(name='s', rangemin=-4, rangemax=4), Dim(name='t', rangemin=-5, rangemax=5)) Values: (-1.5, 2.7)
[23]
Values encode to 'f2x'
[24]
Then decode to: -1.515625 ±0.015625 2.6953125 ±0.0390625

Geohashes

Geohashes are a convenient way of encoding geographical coordinates as a character string.

I have used the same textual mapping so that the following dimensional data will reproduce Geohash en/decoding.

[25]

Unfortunately Geohashes, whilst specifying coordinates in the order latitude then longitude; alternate bits from longitude and then latidude - i.e. in reverse, so adjustments are needed.

[26]
[27]
'gcpue5'
[28]
Geohash 'gcpue5' decodes to: latitude : 51.43524169921875 ±0.00274658 longitude : -0.2142333984375 ±0.00549316

External sites for testing Geohashes

You can enter the coordinates or the geohash above into this site and find that it points to the "Wimpledon All England Lawn Tennis and Croquet Club".

End Notes.

This blog post was written because the WIkipedia entry and a web search did nt turn up a good description of Geohashes for me to finish the Python Rosetta Code entry initially. After more reading I finally completed the Python entry.

I thought about generalising the method to more than two dimensions and also wanted to highlight the swithing of order in geohashing: lat/long normal coord order vs long/lat used when alternating the bits from the binary search.

Another generalisation might be to use n-ary search for a value producing a base-n "bit"; but I could not think of a need, and this post is already long.

END.

Not saved yet

python3 | launching

Posted by Paddy3118 at 4:50 pm 0 comments

Monday, April 06, 2020

Spin the table: Solution!

An answer to a puzzle set by Matt Parker on his YouTube channel:



The Puzzle:


  1. A Circular table with positions 1..7
  2. Delegates numbered 1..7 arranged arbitrarily around the table.
  3. Delegates start with only one match to a correct seating position.
  4. Is it ALWAYS possible to spin the table and get two or more matches?
Find arrangements where it is not possible to spin the table to get more than one match.

The Model.

Let's number the positions on the table with integers 1 .. 7; and also give the delegates numbers 1 .. 7. If the delegate number equals the table number then that constitutes a seating match.

Storing the table order as a fixed list of [1, 2, ..7], the all the rotations of the table  can be modelled by the additional six rotations of that initial list found by popping the last number from the last rotation and appending to the front of the list, i.e:
[[1, 2, .. 7],
 [7, 1, .. 6],
 [6, 7, 1, .. 5],
 ..
 [2, 3, .. 7, 1]]

This is stored in variable spun.

The permutations of delegate numbers 1 .. 7 are all the possible orderings of delegates.
Every "base" perm will have six other perms which are just other rotations of the "base" ordering . If only the permutations generated with an initial 1 are used, (for example), then we can assume that any of the other six rotations of the permutation will also be an answer.

The Code.



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
# -*- coding: utf-8 -*-
"""
Spin the table.

Created on Wed Mar 25 19:59:49 2020

@author: Paddy3118
"""

from itertools import permutations


table = list(range(1, 8))
spun = []                   # spinning table positions
for i in range(7):
    spun.append(table)
    table = table[-1:] + table[:-1]     # rotate once

def match_count(table, delegates):
    return sum(1 for t, d in zip(table, delegates)
               if t == d)

def max_matches(delegates):
    return max(match_count(table, delegates) for table in spun)

def spin_the_table():
    return [d for d in permutations(table)
            if d[0] == 1 and match_count(table, d) == 1 and max_matches(d) < 2]


if __name__ == '__main__':
    answer = spin_the_table()
    print(f"Found {len(answer)} initial arrangements of the delegates:")
    for a in answer:
        print(' ', str(a)[1:-1])
    print("The 6 rotations of each individual answer above are also solutions")

Answer:



Found 19 initial arrangements of the delegates:
  1, 3, 5, 7, 2, 4, 6
  1, 3, 6, 2, 7, 5, 4
  1, 3, 7, 6, 4, 2, 5
  1, 4, 2, 7, 6, 3, 5
  1, 4, 6, 3, 2, 7, 5
  1, 4, 7, 2, 6, 5, 3
  1, 4, 7, 3, 6, 2, 5
  1, 4, 7, 5, 3, 2, 6
  1, 5, 2, 6, 3, 7, 4
  1, 5, 4, 2, 7, 3, 6
  1, 5, 7, 3, 6, 4, 2
  1, 6, 2, 5, 7, 4, 3
  1, 6, 4, 2, 7, 5, 3
  1, 6, 4, 3, 7, 2, 5
  1, 6, 4, 7, 3, 5, 2
  1, 6, 5, 2, 4, 7, 3
  1, 7, 4, 6, 2, 5, 3
  1, 7, 5, 3, 6, 2, 4
  1, 7, 6, 5, 4, 3, 2
The 6 rotations of each individual answer above are also solutions


Matt Parker's Solution Video: (He mentions my solution at time 05:00).

END.

Posted by Paddy3118 at 5:35 pm 0 comments

Saturday, January 25, 2020

Sharing another way?

Tickled!

Sharing came up in something I was reading that got me revisiting the Wikipedia page on the Thue Morse Sequence. Tucked away to the right is a static image of ones and zeroes with the caption:
 "When counting in binary, the digit sum modulo 2 is the Thue-Morse sequence"
Now I had already blogged about the Thue-Morse sequence with respect to sharing fairly:

If you had two people, Asia and Bain chosing from different cash bounties then if Asia went first she'd choose the largest cash bounty, then Bain would chose the next larget, ... if they alternate like this then Asia, who went first, would always get more than Bain; furthermore, Asias total accumulated bounty would increase over Bains over time like this:





In my earlier blog entry: Sharing Fairly, I looked into the Thue-Morse sequence for who gets to pick next and showed that using that sequence, and itsextension for more than two people sharing, the person with the most accumulated wealth switches over time and does not diverge like in the simple taking of turns shown above.

I showed how using a Thue-Morse sequence for multiple people sharing gave a "fairer" distribution of wealth, for example  this graph of accumulated wealth for three people:


Old Code

My code to generate the Thue Morse sequence for sharing between multiple people was this function:


def thue_morse(persons='ABC', mx=20, tutor=True)

Its algorithm works, and was straight forward to describe as a natural extension of the two-person case, but could generate more terms than requested that are then trimmed.

Back to the tickle

 "When counting in binary, the digit sum modulo 2 is the Thue-Morse sequence"
Now that method might be coded to only give the number of terms required. I might code it as a generator.  But as it stands it applies only to two people (binary).

I thought that an equivalent statement for many people might be:
"When counting base b, the digit sum modulo b is the Thue-Morse sequence of fairer sharing between b people"
I set out to code it.

Code Journey

"When counting base b" ... "digit sum" ...

I am going to need to express an integer in different bases, and sum the digits.

Base changing

Normally when expressing an integer in a different base, b,  it is for printout and the routine returns characters: The columns in the string denote, by position, the multiples of successive powers ob b.
For example:
    "123" base 10 == 3 * 10**0 + 2 * 10**1 + 1 * 10**2

If the base goes beyond 10 such as in hexadecimal, then letters A... are used to represent numbers beyond 9 as single characters.

We don't want that!

We need to sum the digits so in the _basechange_int function a list is returned of integers to make for easier summation later.
Integer 123 base 10 would return the list  [1, 2, 3]

People representation

I am used to generating the Thue-Morse sequence as a string ofupper-case characters where each individual characters corresponds to one persons turn in sharing. as in my original, I will code my new function to have a persons parameter which is a string of different characters, each representing a person.
The number of persons becomes the base, b, of my altered  statement above.
The algorithm will generate an int that will then need to be represented as a corresponding character from the persons string.

Generate

I decided to code the algorithm as a generator of successive terms from the given persons. islice is my friend :-)


Thue-Morse sequence generator from digit counts




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
# -*- coding: utf-8 -*-
"""
Created on Sat Jan 18 08:34:29 2020

@author: Paddy3118
"""

#%%

from itertools import count, islice

def _basechange_int(num, b):
    """
    Return list of ints representing positive num in base b

    >>> b = 3
    >>> print(b, [_basechange_int(num, b) for num in range(11)])
    3 [[0], [1], [2], [1, 0], [1, 1], [1, 2], [2, 0], [2, 1], [2, 2], [1, 0, 0], [1, 0, 1]]
    >>>
    """
    if num == 0:
        return [0]
    result = []
    while num != 0:
        num, d = divmod(num, b)
        result.append(d)
    return result[::-1]

def _sumdigits_int(numlist, b):
    "Return the sum of numbers in numlist modulo b, in base b"
    return sum(numlist) % b

def fairshare(persons='ABC'):
    """
    Sequence of fairer sharing between persons
    We have:
      "When counting in binary, the digit sum modulo 2 is the Thue-Morse sequence"

    This computes:
        When counting base b, the digit sum modulo b is the Thue-Morse sequence
        of fairer sharing between b people
    """

    base = len(persons)     # number base
    num2person = dict(zip(range(base), persons))
    for i in count():
        yield num2person[_sumdigits_int(_basechange_int(i, base), base)]


if __name__ == '__main__':
    PERSONS = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
    for b in [10, 2, 3, 5, 7]:
        #print(b, [_basechange_int(num, b) for num in range(11)])
        print(b, ''.join(islice(fairshare(PERSONS[:b]), 44)))


Example runs:



In [24]: b = 3; print(b, ''.join(islice(fairshare(PERSONS[:b]), 44)))
3 ABCBCACABBCACABABCCABABCBCABCACABABCCABABCBC

In [25]: b = 2; print(b, ''.join(islice(fairshare(PERSONS[:b]), 44)))
2 ABBABAABBAABABBABAABABBAABBABAABBAABABBAABBA

In [26]:

Checks

I can't do the maths to formally proove it, but I have compared outputs between the old and the new and they generate the same sequences.

P.S:

I decided to start the "Fairshare between two and more" task on Rosetta Code based on this.

Posted by Paddy3118 at 11:44 am 1 comments

Sunday, November 10, 2019

Quasi-random Sobol sequence



In which I revisit low-discrepancy series, Hunt for an understanding of the Sobol sequence, then translate and optimise a C++ program into Python.

Re-visit VdC

I had left a mental note to get back to the Van der Corput sequence after creating the Rosetta Code task some years ago. Van der Corput is a type of low-discrepancy or quasi random sequence often used in Financial simulations; and with ever larger digital systems on automotive System on Chips, I keep a look-out for Verification techniques to add to the toolbox.

I decided to investigate the Sobol Sequence, understand the maths and write some Python!

Sobol Maths

I followed references; I googled; I read; and I just could not gain an understanding of what is happening, enough to write some code. What code I initially browsed gave me only snippets - Those snippets were bit twiddling so I coded those up for fun.

Fun bit twiddles

Highest 1 bit in binary number

def binhi1(b):
    "highest 1 bit of binary b. (== ceil(log2(b)))"
    pos = 0
    while b > 0:
        pos, b = pos + 1, b // 2
    return pos

Example 1
for n in range(18):
    print(f"{n:>5b} {binhi1(n)}")
    
    0 0
    1 1
   10 2
   11 2
  100 3
  101 3
  110 3
  111 3
 1000 4
 1001 4
 1010 4
 1011 4
 1100 4
 1101 4
 1110 4
 1111 4
10000 5
10001 5

Lowest 0 bit in binary number

def binlo0(b):
    "lowest 0 bit from right of binary b. Always >=1"
    pos, b2 = 1, b // 2
    while b2 * 2 != b:
        pos, b, b2 = pos + 1, b2, b2 // 2
    return pos

Example 2

for n in range(18):
    print(f"{n:<05b binlo0="" font="" n="">)
00000 1
00001 2
00010 1
00011 3
00100 1
00101 2
00110 1
00111 4
01000 1
01001 2
01010 1
01011 3
01100 1
01101 2
01110 1
01111 5
10000 1
10001 2


Code to translate

One of the Sobol references was this one by S. Joe and F. Y. Kuo. They had written academic papers on Sobol and had generated constants for the use of the sequence in up to 21201 dimensions. They also had a C++ program with BSD-style license and, crucially, sample outputs from their C++ code that could aid me in writing  Python conversion of their code.

I downloaded their C++ and their latest set of direction numbers and took a look. Eager to code something, I coded up a reader for their directions file:

Initial dirctions file reader


def direction_file_reader(fname="Sobol_new-joe-kuo-6.21201"):
    """\
    Reads direction file of format:
        d       s       a       m_i     
        2       1       0       1 
        3       2       1       1 3 
        4       3       1       1 3 1 
        5       3       2       1 1 1 
        6       4       1       1 1 3 3
        ...
    """
    with open(fname) as f:
        headers = f.readline().strip().split()
        assert headers == ['d', 's', 'a', 'm_i']
        directions = []
        for line in f:
            d, s, a, *m = [int(cell) for cell in line.strip().split()]
            directions.append((d, s, a, tuple([0] + m)))
        return directions

The direction file was three single columns then all the other columns made up the list m.

C++ code direct translation.

I decided to translate the C++ code near verbatim to Python, then run the example and ensure that the Python gave the same output as stated for the C++ program.

I cut-n-pasted sections of the C++ code as comments/unused triple-quoted strings into my Python then slowly translated from C++.
The C++ used several single-character, uppercase, variable names. I took them to be references to their papers or maybe earlier Fortran programs and either kept them, or made note of their re-naming.

I finally got my first working Python code which is below, but without the C++ code comments I followed:


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
def sobol_points(N, D):
    """From: C++ Source of Frances Y. Kuo "sobol.cc" """
    pow2_32 = 2**32

    L = binhi1(N)
    C = [binlo0(i) for i in range(N)]
    points= [[0 for j in range(D)] for i in range(N)]
    V = [1 << (32- i) for i in range(L+1)]

    X = [0]
    for i in range(1, N):
        X.append(X[i-1] ^ V[C[i-1]])
        points[i][0] = X[i]/pow2_32;
    del V, X
    
    directions = direction_file_reader()
    for j in range(1, D):
        d, s, a, m = directions[j-1]
        #print(d, s, a, m)
        
        if L <= s:
            V = [0] + [m[i] << (32 - i) for i in range(1, L+1)]
        else:
            V = [0] + [m[i] << (32 - i) for i in range(1, s+1)]
            for i in range(s+1, L+1):
                V.append(V[i-s] ^ (V[i-s] >> s))
                for k in range(1, s):
                    V[i] ^= (((a >> (s-1-k)) & 1) * V[i-k])
        X = [0]
        for i in range(1, N):
            X.append(X[i-1] ^ V[C[i-1]])
            points[i][j] = X[i]/pow2_32
    
    return points

Issues with the directly translated Python/C++:

  • The calculation ov V[] in lines 21 to 28 I thought should become a function.
  • All points are computed and returned at once. It would be better as a generator function.
  • The calculation of X[] in lines 29 to 31 seems like a very inefficient loop within a loop. Furthermore, if turned into a generator; Iwould need to unscramble the loops to try and recalculate X for each point without another loop.

Final Python generator function version


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
# -*- coding: utf-8 -*-
"""sobol_point_gen.py

 Usage:
    sobol_point_gen.py N D <dimensions_file>
    sobol_point_gen.py -h | --help

 Options:
     -h --help     Show this screen.     

 Arguments:
     N  Number of points to generate (<2**32)
     D  Number of dimensions for each point. (D <= max D in dimensions_file)
     dimensions_file    See: https://web.maths.unsw.edu.au/~fkuo/sobol/

 Created from an original c++ source by Frances Y. Kuo and Stephen Joe

Created on Fri Nov  1 09:13:12 2019

@author: Paddy3118
"""

from functools import lru_cache

#%%
def binhi1(b):
    "highest 1 bit of binary b. (== ceil(log2(b)))"
    pos = 0
    while b > 0:
        pos, b = pos + 1, b // 2
    return pos

def binlo0(b):
    "lowest 0 bit from right of binary b. Always >=1"
    pos, b2 = 1, b // 2
    while b2 * 2 != b:
        pos, b, b2 = pos + 1, b2, b2 // 2
    return pos

#%%
def direction_file_reader(d, fname="Sobol_new-joe-kuo-6.21201"):
    """\
    Reads d+header lines from direction file of format:
        d       s       a       m_i     
        2       1       0       1 
        3       2       1       1 3 
        4       3       1       1 3 1 
        5       3       2       1 1 1 
        6       4       1       1 1 3 3
        7       4       4       1 3 5 13 
        8       5       2       1 1 5 5 17 
        9       5       4       1 1 5 5 5 
        10      5       7       1 1 7 11 19 
        11      5       11      1 1 5 1 1 
        12      5       13      1 1 1 3 11 
        13      5       14      1 3 5 5 31 
        14      6       1       1 3 3 9 7 49 
        15      6       13      1 1 1 15 21 21 
        16      6       16      1 3 1 13 27 49 
        17      6       19      1 1 1 15 7 5 
        18      6       22      1 3 1 15 13 25 
        19      6       25      1 1 5 5 19 61 
        20      7       1       1 3 7 11 23 15 103 
        21      7       4       1 3 7 13 13 15 69 
        22      7       7       1 1 3 13 7 35 63 
        23      7       8       1 3 5 9 1 25 53 
        24      7       14      1 3 1 13 9 35 107 
        25      7       19      1 3 1 5 27 61 31 
        26      7       21      1 1 5 11 19 41 61 
        27      7       28      1 3 5 3 3 13 69 
        28      7       31      1 1 7 13 1 19 1 
        29      7       32      1 3 7 5 13 19 59 
        30      7       37      1 1 3 9 25 29 41 
        31      7       41      1 3 5 13 23 1 55 
        32      7       42      1 3 7 3 13 59 17 
        ...
    """
    with open(fname) as f:
        headers = f.readline().strip().split()
        assert headers == ['d', 's', 'a', 'm_i']
        directions = []
        for i, line in zip(range(d), f):
            d, s, a, *m = [int(cell) for cell in line.strip().split()]
            directions.append((d, s, a, tuple([0] + m)))
        return directions

#%%
@lru_cache(128)
def calc_dnum(s, a, m, mxbits):                      # V, direction numbers
    if mxbits <= s:
        dirnum = [0] + [m[i] << (32 - i) for i in range(1, mxbits+1)]
    else:
        dirnum = [0] + [m[i] << (32 - i) for i in range(1, s+1)]
        for i in range(s+1, mxbits+1):
            dirnum.append(dirnum[i-s] ^ (dirnum[i-s] >> s))
            for k in range(1, s):
                dirnum[i] ^= (((a >> (s-1-k)) & 1) * dirnum[i-k])
    return dirnum

def sobol_point_gen(N, D, fname="Sobol_new-joe-kuo-6.21201"):
    pow2_32 = 2**32
    mxbits = binhi1(N)                              # L, max bits needed
    dnum = [1 << (32- i) for i in range(mxbits+1)]  # V, direction numbers
    directions = direction_file_reader(D, fname)
    xd = [0 for j in range(D)]                      # X
    pt = [0 for j in range(D)]                      # POINTS
    yield pt
    for n in range(1, N):
        nfirst0 = binlo0(n-1)                       # C
        pt = [0 for j in range(D)]
        xd[0] = xd[0] ^ dnum[nfirst0]
        pt[0] = xd[0] / pow2_32
        for jj in range(1, D):
            _, s, a, m = directions[jj-1]
            dnum2 = calc_dnum(s, a, m, mxbits)      # V, direction numbers
            xd[jj] = xd[jj] ^ dnum2[nfirst0]
            pt[jj] = xd[jj]/pow2_32
        yield pt

#%%
if __name__ == '__main__':
    first_sobol_3D_points = [
         [0, 0, 0],
         [0.5, 0.5, 0.5],
         [0.75, 0.25, 0.25],
         [0.25, 0.75, 0.75],
         [0.375, 0.375, 0.625],
         [0.875, 0.875, 0.125],
         [0.625, 0.125, 0.875],
         [0.125, 0.625, 0.375],
         [0.1875, 0.3125, 0.9375],
         [0.6875, 0.8125, 0.4375]]
    generated = list(sobol_point_gen(10, 3))
    assert generated == first_sobol_3D_points


Function def direction_file_reader is changed to only read as much of thedirections file as is necessary.
Function def sobol_point_gen is now a generator function in which successive calls to next give successive points. The loops are swapped inside with an outer loop yielding a point every iteration, and an inner loop over dimensions > 1.
Function def calc_dnum is cached. Initial suspicions of a lot of recalculation proved true:

calc_dnum cache stats


In [1]: g = list(sobol_point_gen(10000, 3))
In [2]: calc_dnum.cache_info()
Out[2]: CacheInfo(hits=20012, misses=4, maxsize=128, currsize=4)



Posted by Paddy3118 at 11:45 am 1 comments
Subscribe to: Posts (Atom)

Followers

Subscribe Now: google

Add to Google Reader or Homepage

Go deh too!

whos.amung.us

Blog Archive