In [1]:

    

from __future__ import print_function
import sdm as sdmlib
from sdm import utils
import numpy as np
import matplotlib.mlab as mlab
from IPython.display import clear_output
import matplotlib.pyplot as plt


    

In [2]:

    

bits, radius = 1000, 451
#bits, radius = 256, 103

p = utils.calculate_probability(bits, radius)

sample = 1000000
scanner_type = sdmlib.SDM_SCANNER_OPENCL

address_space = sdmlib.AddressSpace.init_random(bits, sample)
address_space.opencl_init();


    

In [3]:

    

p


    

    
Out[3]:





0.0010718500489241473

In [4]:

    

counter = []


    

In [5]:

    

for i in range(100):
    if i%100 == 0:
        clear_output(wait=True)
        print(i)
    bs = sdmlib.Bitstring.init_random(bits)
    result = address_space.scan_opencl2(bs, radius)
    for idx in result:
        bs2 = address_space.get_bitstring(idx)
        counter.append(bs.distance_to(bs2))


    

    




0

In [6]:

    

from math import factorial
comb = lambda a, b: factorial(a) // factorial(b) // factorial(a-b)
acc = [0]
for i in range(radius+1):
    acc.append(acc[-1] + comb(bits, i))
acc = acc[1:]


    

In [7]:

    

x = []
y = []
for d in range(0, radius+1):
    x.append(d + 0.5)
    y.append(((1.0*comb(bits, d)) / acc[radius]))


    

In [8]:

    

plt.figure(figsize=(8, 6), dpi=100)
plt.hist(counter, bins=range(0, radius+2), density=True)
plt.plot(x, y, 'r', linewidth=2.0)
plt.xlim(400, radius+5)
plt.ylabel('Probability')
plt.xlabel('Distance from the center')
plt.show()


    

In [9]:

    

from math import log
I = [-log(p)/log(2) if p > 0 else 0 for p in y]


    

In [10]:

    

plt.figure(figsize=(8, 6), dpi=100)
plt.plot(x, I, 'r', linewidth=2.0)
plt.xlim(420, radius)
plt.ylim(0, 15)
plt.ylabel('Self information')
plt.xlabel('Distance to the center')
plt.grid()
plt.show()


    

In [11]:

    

plt.figure(figsize=(8, 6), dpi=100)
plt.plot(x, [int(p) for p in I], 'r', linewidth=2.0)
plt.xlim(420, radius)
plt.ylim(0, 15)
plt.ylabel('Self information')
plt.xlabel('Distance to the center')
plt.grid()
plt.show()


    

In [12]:

    

I[-1]


    

    
Out[12]:





2.3753165837275243

In [13]:

    

1.0*int(I[420]) / int(I[-1])


    

    
Out[13]:





6.5

In [14]:

    

repr([int(p) for p in I])


    

    
Out[14]:





'[990, 980, 971, 962, 954, 947, 939, 932, 925, 918, 912, 905, 899, 893, 887, 881, 875, 869, 863, 857, 852, 846, 841, 835, 830, 825, 819, 814, 809, 804, 799, 794, 789, 784, 779, 775, 770, 765, 761, 756, 751, 747, 742, 738, 733, 729, 725, 720, 716, 712, 707, 703, 699, 695, 691, 687, 682, 678, 674, 670, 666, 662, 659, 655, 651, 647, 643, 639, 636, 632, 628, 624, 621, 617, 613, 610, 606, 603, 599, 595, 592, 588, 585, 581, 578, 575, 571, 568, 564, 561, 558, 554, 551, 548, 544, 541, 538, 535, 532, 528, 525, 522, 519, 516, 513, 510, 507, 503, 500, 497, 494, 491, 488, 485, 482, 479, 477, 474, 471, 468, 465, 462, 459, 456, 454, 451, 448, 445, 442, 440, 437, 434, 431, 429, 426, 423, 421, 418, 415, 413, 410, 408, 405, 402, 400, 397, 395, 392, 390, 387, 385, 382, 380, 377, 375, 372, 370, 367, 365, 363, 360, 358, 355, 353, 351, 348, 346, 344, 341, 339, 337, 335, 332, 330, 328, 326, 323, 321, 319, 317, 314, 312, 310, 308, 306, 304, 302, 299, 297, 295, 293, 291, 289, 287, 285, 283, 281, 279, 277, 275, 273, 271, 269, 267, 265, 263, 261, 259, 257, 255, 253, 251, 249, 247, 246, 244, 242, 240, 238, 236, 235, 233, 231, 229, 227, 225, 224, 222, 220, 218, 217, 215, 213, 211, 210, 208, 206, 205, 203, 201, 200, 198, 196, 195, 193, 191, 190, 188, 187, 185, 183, 182, 180, 179, 177, 176, 174, 173, 171, 170, 168, 167, 165, 164, 162, 161, 159, 158, 156, 155, 153, 152, 150, 149, 148, 146, 145, 143, 142, 141, 139, 138, 137, 135, 134, 133, 131, 130, 129, 127, 126, 125, 124, 122, 121, 120, 118, 117, 116, 115, 114, 112, 111, 110, 109, 108, 106, 105, 104, 103, 102, 101, 99, 98, 97, 96, 95, 94, 93, 92, 90, 89, 88, 87, 86, 85, 84, 83, 82, 81, 80, 79, 78, 77, 76, 75, 74, 73, 72, 71, 70, 69, 68, 67, 66, 65, 64, 64, 63, 62, 61, 60, 59, 58, 57, 56, 56, 55, 54, 53, 52, 51, 51, 50, 49, 48, 47, 47, 46, 45, 44, 43, 43, 42, 41, 40, 40, 39, 38, 38, 37, 36, 35, 35, 34, 33, 33, 32, 31, 31, 30, 29, 29, 28, 28, 27, 26, 26, 25, 25, 24, 23, 23, 22, 22, 21, 21, 20, 19, 19, 18, 18, 17, 17, 16, 16, 15, 15, 14, 14, 13, 13, 13, 12, 12, 11, 11, 10, 10, 10, 9, 9, 8, 8, 8, 7, 7, 6, 6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 3, 2, 2, 2]'

In [ ]: