First commit

network/acceptance_rejection.py

@@ -0,0 +1,22 @@
+from network.network import (get_number_of_activated_antenna,
+                             is_valid_network)
+
+
+def _get_random_networks(box_size, random):
+    '''yields random networks, broken or not'''
+    number_of_square = box_size * box_size
+    while True:
+        yield random.randint(0, 2 ** number_of_square - 1)
+
+
+def get_random_networks(box_size, random):
+    for network in _get_random_networks(box_size, random):
+        if is_valid_network(network, box_size):
+            yield network
+
+
+def estimate_number_of_activated_antenna(n, box_size, random):
+    estimate = 0
+    for _, network in zip(range(n), get_random_networks(box_size, random)):
+        estimate += get_number_of_activated_antenna(network)
+    return estimate / n

network/brute_force.py

@@ -0,0 +1,27 @@
+from network.network import is_valid_network, get_number_of_activated_antenna
+
+def _get_all_possible_network(box_size):
+    '''returns all networks, broken or not'''
+    return range( 2**(box_size*box_size) )
+
+
+def get_all_possible_network(box_size):
+    if box_size < 1:
+        return
+    for network in _get_all_possible_network(box_size):
+        if is_valid_network(network, box_size):
+            yield network
+
+
+def _get_random_networks(n, box_size, random):
+    all_networks = list(get_all_possible_network(box_size))
+    n = min(n, len(all_networks))
+    return random.sample(all_networks, n), n
+
+
+def estimate_number_of_activated_antenna(n, box_size, random):
+    networks, n = _get_random_networks(n, box_size, random)
+    estimate = 0
+    for network in networks:
+        estimate += get_number_of_activated_antenna(network)
+    return estimate / n

network/gibbs.py

@@ -0,0 +1,36 @@
+from network.network import (get_number_of_activated_antenna,
+                             can_antenna_be_turned_on,
+                             turn_off_antenna,
+                             turn_on_antenna)
+
+def _flip_coin(random):
+    return random.random() < 0.5
+
+def create_new_network(network, box_size, random):
+    number_of_square = box_size * box_size
+    for antenna in range(number_of_square):
+        if can_antenna_be_turned_on(antenna, network, box_size) and _flip_coin(random):
+            network = turn_on_antenna(antenna, network)
+        else:
+            network = turn_off_antenna(antenna, network)
+    return network
+
+
+def get_pseudo_random_networks(box_size, random):
+    network = 0
+    while True:
+        network = create_new_network(network, box_size, random)
+        yield network
+
+def get_random_networks(box_size, random):
+    time_between_sample = 20
+    for i, network in enumerate(get_pseudo_random_networks(box_size, random), 1):
+        if i % time_between_sample == 0:
+            yield network
+
+
+def estimate_number_of_activated_antenna(n, box_size, random):
+    estimate = 0
+    for _, network in zip(range(n), get_random_networks(box_size, random)):
+        estimate += get_number_of_activated_antenna(network)
+    return estimate / n

network/importance_sampling.py

@@ -0,0 +1,92 @@
+from random import Random
+from network.network import create_network, get_number_of_activated_antenna
+
+
+class AntennaGenerator(object):
+
+    def __init__(self, box_size, random):
+        self.box_size = box_size
+        self._antennas = list(range(box_size * box_size))
+        self.random = random
+
+    def get_random_antenna(self):
+        idx, prob = self.get_random_index_and_its_probability()
+        if idx == -1:
+            self._antennas = []
+            return -1, prob
+        choosen_antenna = self._antennas[idx]
+        self._remove_antenna(idx)
+        return choosen_antenna, prob
+
+
+class SuccesiveAntennaGenerator(AntennaGenerator):
+    def _remove_antenna(self, idx):
+        ant = self._antennas[idx]
+        if ant == -1:
+            self._antennas = []
+            return ant
+        if (ant + 1) % self.box_size:
+            self._antennas = self._antennas[idx + 2:]
+        else:
+            self._antennas = self._antennas[idx + 1:]
+
+        if ant + self.box_size < self.box_size * self.box_size:
+            self._antennas.remove(ant + self.box_size)
+
+    def get_random_index_and_its_probability(self):
+        idx = self.random.randint(-1, len(self._antennas)-1)
+        prob = (len(self._antennas) + 1) ** -1
+        return idx, prob
+
+
+class BadAntennaGenerator(SuccesiveAntennaGenerator):
+    '''this antenna generator does not return 'real' probabilities
+    and thus cannot be used in importance_sampling'''
+    def _remove_antenna(self, idx):
+        ant = self._antennas[idx]
+        to_remove = {ant, ant + self.box_size, ant - self.box_size}
+        if (ant + 1) % self.box_size:
+            to_remove.add(ant + 1)
+
+        if ant % self.box_size:
+            to_remove.add(ant - 1)
+        self._antennas = [a for a in self._antennas if a not in to_remove]
+
+
+class BetterSuccesiveAntennaGenerator(SuccesiveAntennaGenerator):
+
+    def _get_random_index(self, max_number_inclusive):
+        if max_number_inclusive < 1:
+            return self.random.randint(-1, max_number_inclusive)
+        mean_ = max_number_inclusive/2
+        x = int(self.random.expovariate(mean_**-1))
+        if x > max_number_inclusive:
+            return -1
+        return x
+
+    def get_random_index_and_its_probability(self):
+        idx = self._get_random_index(len(self._antennas)-1)
+        prob = (len(self._antennas) + 1) ** -1
+        return idx, prob
+
+
+def get_random_network(box_size, random, Generator):
+    antenna_generator = Generator(box_size, random)
+    antennas = []
+    antenna, network_prob = antenna_generator.get_random_antenna()
+
+    while antenna != -1:
+        antennas.append(antenna)
+        antenna, p = antenna_generator.get_random_antenna()
+        network_prob *= p
+
+    return create_network(antennas), network_prob
+
+
+def estimate_number_of_activated_antenna(n, box_size, random, Generator):
+    ret, sum_of_probs = 0, 0
+    for i in range(n):
+        network, p = get_random_network(box_size, random, Generator)
+        ret += get_number_of_activated_antenna(network) / p
+        sum_of_probs += p**-1
+    return ret / sum_of_probs

network/network.py

@@ -0,0 +1,97 @@
+from collections import deque
+
+
+def _get_print_symbol(network, position):
+    black_square = '\u25A0'
+    white_square = '\u25A1'
+    activated_antennas = set(get_activated_antennas(network))
+    if position in activated_antennas:
+        return black_square
+    return white_square
+
+
+def pp_network(network, box_size):
+    ret = ""
+    for line in range(box_size):
+        for column in range(box_size):
+            symbol =  _get_print_symbol(network, line * box_size + column)
+            ret += symbol + '-'
+        ret = ret[:-1] + '\n'
+    return ret
+
+
+def get_activated_antennas(network):
+    i = 0
+    while network > 0:
+        if network % 2:
+            yield i
+        network = network >> 1
+        i += 1
+
+
+def create_network(activated_antennas):
+    return sum(2**activated_antenna for activated_antenna in activated_antennas)
+
+
+def _are_two_consicutive_horizontaly_antennas_activated(network, box_size):
+    predecessor = -2
+    for antenna in get_activated_antennas(network):
+        if antenna == predecessor + 1:
+            return True
+        if (antenna + 1) % box_size:
+            predecessor = antenna
+    return False
+
+
+def _are_two_consicutive_verticaly_antennas_activated(network, box_size):
+    invalid_entennas = set()
+    for antenna in get_activated_antennas(network):
+        if antenna in invalid_entennas:
+            return True
+        invalid_entennas.add(antenna + box_size)
+
+    return False
+
+
+def is_valid_network(network, box_size):
+    if _are_two_consicutive_horizontaly_antennas_activated(network, box_size):
+        return False
+
+    if _are_two_consicutive_verticaly_antennas_activated(network, box_size):
+        return False
+    return True
+
+
+def has_antenna(antenna, network):
+    return (network >> antenna) % 2
+
+
+def turn_on_antenna(antenna, network): # error with borders
+    if has_antenna(antenna, network):
+        return network
+    return network + 2**antenna
+
+
+def turn_off_antenna(antenna, network):
+    if has_antenna(antenna, network):
+        return network - 2**antenna
+    return network
+
+
+def can_antenna_be_turned_on(antenna, network, box_size):
+    for ant in get_activated_antennas(network):
+        if ant == antenna + 1 and ant % box_size:
+            return False
+        if ant == antenna - 1 and antenna % box_size:
+            return False
+        if ant == antenna + box_size:
+            return False
+        if ant == antenna - box_size:
+            return False
+        if ant > antenna + box_size:
+            return True
+    return True
+
+
+def get_number_of_activated_antenna(network):
+    return len(list(get_activated_antennas(network)))