HW1 for EE5393: Circuits, Computation and Biology

cell_tower-save.pickle

cnumpy.core.multiarray
_reconstruct
p0
(cnumpy
ndarray
p1
(I0
tp2
S'b'
p3
tp4
Rp5
(I1
(I21
I4
tp6
cnumpy
dtype
p7
(S'f8'
p8
I0
I1
tp9
Rp10
(I3
S'<'
p11
NNNI-1
I-1
I0
tp12
bI00
S'\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x004@\x10-\xd6\x83-\xfek?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x00>@k\x18\x0f\x0e\\\xcc\xa2?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x80A@\x13\xa6^\x843\x02\xb8?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x80B@\x88CcjU\xb0\xc0?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x00D@f\xe9\xbc{S\x9a\xd0?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x00E@\xf8<ic\xe0\xae\xd4?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x80F@\xb5\xcf\x82\x82#b\xdc?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x00H@\xf2\xea\x1c\xa0\x0e\xb6\xe0?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x00I@\xdac\xdag\xb7\xea\xe4?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x00J@\x0e\x14\xc7\xb0\xfc\x94\xe2?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x80K@\x95\x1f.\x8bZ2\xe8?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x80L@\xe6\x96\x92^`\x94\xea?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x00N@\xfb\xa6\x06B=*\xeb?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x00O@\xa3$P2\xee\xab\xec?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00@P@\n\xf2\xfe\x1a\xec\x8f\xee?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x80Q@\xb5\xce1\xb2\xde\xb8\xef?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x00T@H{\xc3\x99\x04#\xef?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x80V@\xa6\xdd\xf9;\x85\xfd\xef?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x00Y@K\x13\xd3\r\xe2\xff\xef?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x00i@8\xb8\xbe\xb5\x9b\xf7\xef?\x00\x00\x00\x00\x00\x00Y@\x00\x00\x00\x00\x00\x00y@\x00\x00\x00\x00\x00\x00y@\xf5\xde\xccI\x0b\xf0\xef?'
p13
tp14
b.

hw1.py

from __future__ import division
from pylab import *
from sympy import init_printing, init_session
from sympy import var, simplify, Poly, pprint
from scipy.special import comb
from sympy import lambdify # eval our function numerically
import pickle


"""
README:
    This code is not well commented. This code is for development only; not
    nearly as well formatted or as well presented as the notebook. Additionally,
    it doesn't give the exact same results; the notebook results are better
    thought out etc. Go with that, but if you want to run it...
"""

__author__ = "Scott Sievert [email protected]"
__class__ = "EE5393: Circuits, Computation and Biology"
__prof__ = "Marc Riedel"

def prob_1():
    def bern(n, i, t): return comb(n, i) * t**i * (1-t)**(n-i)
    def calc_B(n, debug=False):
        """ calculates the Berstein basis for polynomials of degree n """
        B = zeros((n+1,n+1))
        for i in arange(n+1):
            B_in = comb(n, i) * t**i * (1-t)**(n-i)
            if debug: pprint(B_in)
            B_in = B_in.expand()
            a = Poly(B_in, t)
            B[i, :] = a.all_coeffs()
        return B

    def convert_to_berstein(x, t):
        """ represents a polynomial in the berstein basis. assumes x(t)=sum_i c_i t^i
        return the polynomials of the function.

        [1 2 3] corresponds to t^2 + 2*t^1 + 3*t^0

        (index i corresponds to B_{i+1, n})
        """
        for m in arange(0, 10):
            # define our n
            xz = concatenate((zeros(m), x))
            n = xz.shape[0] - 1

            B = calc_B(n)
            b = inv(B).dot(xz)
            if (b > -1).all(): break
        g = 0
        for i in arange(n+1):
            g += b[i] * bern(n, i, t)
        # returns the symbolic expression and the coeffs in front of each term
        return g, b

    t = var('t')
    # (1 - 1*t + 2*t^2 - 2*t^3 + 3*t^4 - 3*t^5) / 11
    c = array([1, -1, 2, -2, 3, -3]) * 1 / 11

    g, b = convert_to_berstein(c[::-1], t)

    # checking they are equal
    f = 0
    for i in arange(c.shape[0]): f += c[i] * pow(t, i) * 1.0

    g = g.expand()
    print f==g # within machine epsilon (ie, 0.999999*t**2)
    pprint(f)
    pprint(g)

    def eval_poly(t):
        """ Evaluate the polynomial using the circuit model shown on paper (mux with
        sum block as adder)
        """
        N = 1e4
        X = zeros((5, N))
        for i in arange(5):
            X[i, :] = around(rand(N) - 1/2 + t)
        # which input line do we want to select?
        select = sum(X, axis=0)
        imp = 0
        for i in select:
            imp += b[i]
        imp /= N
        return imp

    x = []
    t_values = linspace(0, 1)
    for t in t_values:
        x += [eval_poly(t)]

    t = var('t')
    true_f = lambdify(t, f)(linspace(0, 1))

    figure()
    plot(t_values, true_f, label='Ground truth')
    plot(t_values, x, 'o--', label='Simulation')
    legend(loc='best')
    show()

    figure()
    plot(t_values, abs(true_f - x), label='Error')
    legend(loc='best')
    title('Error, $|x - \\widehat{x}|$')
    show()

def prob_2():
    """
    My final result:
        x | y | z | P_e_c | P_e_s   
        0 | 0 | 0 | 0.179 | 0.252
        1 | 0 | 0 | 0.180 | 0.747
        0 | 1 | 0 | 0.180 | 0.771
        1 | 1 | 0 | 0.820 | 0.165
        0 | 0 | 1 | 0.819 | 0.748
        1 | 0 | 1 | 0.180 | 0.748
        0 | 1 | 1 | 0.180 | 0.229
        1 | 1 | 1 | 0.179 | 0.166
    
    It seems that a formula for this is

    (!x means `not x`. This operation is not x!)
    """
    def print_in_binary(a):
        print np.array(map(bin, a.flatten())).reshape(a.shape)
    def circuit_output(x, y, z, failure_prob=0.0):
        # & and, ^ xor, | or
        g1 = corrupt(x & y, failure_prob);
        g2 = corrupt(x ^ y, failure_prob)
        g3 = corrupt(g2 & z, failure_prob)
        g4 = corrupt(g1 | g3, failure_prob)
        g5 = corrupt(z ^ g2, failure_prob)
        c, s = g4, g5
        return c, s
    def corrupt(x, p):
        if type(x)==ndarray:
            for i in arange(x.shape[0]):
                if rand() < p: x[i] = 1 - x[i]
            return x
        else:
            if rand() < p: return 1-x
            else: return x

    N = 8
    n = arange(N, dtype=int)
    x = (1 & n) // 1 # 0 1 0 1 0 1 0 1
    y = (2 & n) // 2 # 0 0 1 1 0 0 1 1
    z = (4 & n) // 4 # 0 0 0 0 1 1 1 1

    M = 1e5
    epsilon = 0.05
    print "Given the inputs x, y, z and the error rate (%0.2f) what's the probability of getting the incorrect result for c or s? (probability of error in s/c or P_e_s/P_e_c" % (epsilon)
    print "x | y | z | P_e_c | P_e_s   "
    for i in arange(N):
        x_m = ones(M, dtype=int) * x[i] # x_m for x_multiple
        y_m = ones(M, dtype=int) * y[i]
        z_m = ones(M, dtype=int) * y[i]

        c_correct, s_correct = circuit_output(x[i], y[i], z[i])
        c_corrupt, s_corrupt = circuit_output(x_m, y_m, z_m, failure_prob=epsilon)

        c_wrong = argwhere(c_correct != c_corrupt).shape[0]
        s_wrong = argwhere(s_correct != s_corrupt).shape[0]

        P_e_s = c_wrong / M # Prob_error_s
        P_e_c = s_wrong / M # Prob_error_c
        print "%d | %d | %d | %0.3f | %0.3f" % (x[i], y[i], z[i], 1*P_e_c, 1*P_e_s)
def prob_3():
    n = arange(8, dtype=int)
    a = (n & 1) // 1
    b = (n & 2) // 2
    c = (n & 4) // 4
    d = (n & 8) // 8
    e = (n & 16) // 16

    n = lambda(x): 1-x
    f = (a&b&c) ^ (n(a)&n(b)&n(c)) ^ (b&n(c)) ^ (a&n(b))
    print "a ", a
    print "b ", b
    print "c ", c
    print "f ", f

def prob_3():
    def xor_bar(k, dtype=int):
        if k==1: return array([[-1]], dtype=dtype)
        else:
            lhs = xor(k-1)
            rhs = xor_bar(k-1)
            bottom = hstack((lhs, rhs))
            top = k * ones_like(bottom[0])
            top[top.shape[0]/2:] *= -1
            return vstack((top, bottom))
    def xor(k, dtype=dtype):
        if k==1: return array([[1]], dtype=dtype)
        else:
            lhs = xor_bar(k-1)
            rhs = xor(k-1)
            bottom = hstack((lhs, rhs))
            top = k * ones_like(bottom[0])
            top[top.shape[0]/2:] *= -1
            return vstack((top, bottom))
    def part_d():
        print "xor_3:\n", xor(3)
        print "xor_5:\n", xor(5)
    def part_e():
        """
        Because the block diagram for the AND function is relatively simple and
        we have 
            ab (+) bc (+) cd ...
        we can call x=ab, y=bc so we can use part (d). This function is just a
        fancy way of 

        1 corresponds to x_1 and -1 corresponds to \\bar{x}_1

        We expand our array to hold the values for the block diagram (a, b). The
        block for ab is 2x1, NOT(ab) is 1x2. We can implement this with a 2x2
        block and have the row/cols (accordingly) be the same.
        """
        x = xor(4, dtype=object)
        N = x.shape[0] # number of 2x2 blocks
        x = repeat(x, 2, axis=0)
        x = repeat(x, 2, axis=1)

        # what pairs do we have included? This represents x_1*x_2
        pairs = array([[1, 2], [2, 3], [3, 4], [5, 1]])
        for i in arange(x.shape[0]/2):
            for j in arange(x.shape[1]/2):
                block = x[2*i:2*(i+1), 2*j:2*(j+1)]
                idx = int(abs(block[0, 0])) - 1
                if sum(block) < 0:
                    block[:, 0] = -pairs[idx, 0]
                    block[:, 1] = -pairs[idx, 1]
                elif sum(block) > 0:
                    block[0, :] = pairs[idx, 0]
                    block[1, :] = pairs[idx, 1]
                else: assert False, "Something went terribly wrong"
                x[2*i:2*(i+1), 2*j:2*(j+1)] = block
        print x
    part_d()
    part_e()


def part_4a():
    def calculate_connectedness(N, M, P, SEED=42):
        x = arange(N*M*10) % N

        # to keep the same neighbors
        seed(42)
        np.random.shuffle(x)

        seed(int(SEED))
        for i in arange(x.shape[0]):
            if rand() > P: x[i] = -1
        x = reshape(x[:N*M], (N, M))

        # x is the array of neighbors
        # now see if we're connected from nodes 0 to N-1
        # the neighbors are at x[0]
        # bounce around randomly and see if we reach N-1
        for m in arange(M):
            next_neigh = m
            connected = False
            for its in arange(4*N):
                next_neigh = x[next_neigh, int(rand() * M)]
                if next_neigh == N-1: 
                    connected = True
                    break
        return connected
    def calc_prob(P, N_PROB=5e2):
        c = 0
        if P==0: return 0
        for i in arange(N_PROB):
            if calculate_connectedness(N, M, P, SEED=i): c += 1
        c /= N_PROB
        return c

    N = 100 # nodes
    M = 10 # neighbors
    P = 0.7 # prob connected to neighbor

    MODEL = 'load'
    if MODEL=='train':
        N_PROB = 25
        probs = linspace(0, 1, num=N_PROB)
        p2 = zeros_like(probs)
        for i, prop in enumerate(probs):
            p2[i] = calc_prob(prop)
            print prop, p2[i]
        pickle.dump(probs, open('prop.pickle', 'wb'))
        pickle.dump(p2, open('p2.pickle', 'wb'))
    elif MODEL=='load':
        probs = pickle.load(open('prop.pickle', 'rb'))
        p2 =   pickle.load(open('p2.pickle', 'rb'))

        p = linspace(0, 1)
        def f(p): return 1 / (1 + exp(-7*(p-0.55)))

        figure()
        plot(probs, p2)
        plot(p, f(p))
        show()


def part_4b():
    def neighoring_pixel_in_patch(im, x, y, value=2/3):
        in_patch = False
        idx_values = [[x,y-1], [x,y+1],       \
                      [x-1, y], [x+1,y],      \
                      [x+1, y+1], [x+1, y-1], \
                      [x-1, y+1], [x-1, y-1]  \
                      #[x-2, y+2], [x+2, y], [y, x+1],
                      #[x-2, y], [x-1, y-1], [x-2, y+2]\
                      ]
        for idx in idx_values:
            if im[idx[0], idx[1]] == value:
                in_patch = True
                break
        return in_patch
    def make_cell_map(N, R, Q, CONNECTED=2/3, PWR=1/3, N_R=1024):
        cell_map = zeros((N_R+1, N_R+1))

        # making the circle
        r = arange(-N_R/2, N_R/2)
        r_x, r_y = meshgrid(r, r)
        i = argwhere(r_x**2 + r_y**2 < R**2)
        cell_map[i[:,0], i[:,1]] = 1

        xc, yc = N_R/2, N_R/2
        # putting N random users down
        shuffle(i)
        users = i[0:N, :]
        cell_map[users[:, 0], users[:, 1]] = PWR
        for u in users:
            u -= N_R/2 # center the input
            i = argwhere((r_x-u[0])**2 + (r_y-u[1])**2 < Q**2)
            cell_map[i[:, 0], i[:, 1]] = PWR

        cell_map[yc, xc] = CONNECTED
        yc, xc = yc, xc
        rect = zeros_like(cell_map)
        rs = arange(0, N_R/2)
        rs = repeat(rs, 3)
        for r in rs:
            if False and r%100==0: print r
            rect = zeros_like(cell_map)
            rect[yc-r-1:yc+r+1,   xc-r] = 1
            rect[yc-r-1:yc+r+1,   xc+r] = 1
            rect[yc-r, xc-r-1:xc+r+1] = 1
            rect[yc+r, xc-r-1:xc+r+1] = 1
            i = argwhere(rect==1)
            for y, x in i:
                if cell_map[y, x]==PWR and neighoring_pixel_in_patch(cell_map, y, x):
                    cell_map[y, x] = CONNECTED
        connected = argwhere(cell_map == CONNECTED).shape[0]
        pwr = argwhere(cell_map == PWR).shape[0]
        p_connected = connected / (pwr + connected)
        return cell_map, p_connected
    def test_cell_map(N, R, Q, iterations=10):
        Ps = []
        result = []
        for i in arange(iterations):
            cell_map, p_connected = make_cell_map(N, R, Q, N_R=N_R)
            result += [[N, R, Q, p_connected]]
        result = asarray(result)
        return mean(result[:, 3])

    N, R, Q  = 100, 400, 50
    N_R = 1024
    MODEL = 'load'
    if MODEL=='train':
        params = array([
            [100, 400, 20, 0],
            [100, 400, 30, 0],
            [100, 400, 35, 0],
            [100, 400, 37, 0],
            [100, 400, 40, 0],
            [100, 400, 42, 0],
            [100, 400, 45, 0],
            [100, 400, 48, 0],
            [100, 400, 50, 0],
            [100, 400, 52, 0],
            [100, 400, 55, 0],
            [100, 400, 57, 0],
            [100, 400, 60, 0],
            [100, 400, 62, 0],
            [100, 400, 65, 0],
            [100, 400, 70, 0],
            [100, 400, 80, 0],
            [100, 400, 90, 0],
            [100, 400, 100, 0],
            [100, 400, 200, 0],
            [100, 400, 400, 0],
            ], dtype=float)
        for i in arange(params.shape[0]):
            N, R, Q, p = params[i]
            print "##### ", N, R, Q
            p = test_cell_map(N, R, Q, iterations=40)
            params[i, -1] = p
            pickle.dump(params, open("cell_tower.pickle", "wb"))

        pickle.dump(params, open("cell_tower.pickle", "wb"))
    elif MODEL=='test':
        cell_map, p = make_cell_map(N, R, Q, N_R=N_R)
        figure()
        imshow(cell_map)
        savefig('cell_map.png', dpi=300)
        show()
    elif MODEL=='load':
        results = pickle.load(open("cell_tower-save.pickle", "rb"))

        r = linspace(0, 1)
        p_r = 1 / (1 + exp(-60*(r-0.12)))

        figure()
        title('Probabilty of connectedness vs radius size, N=%d' % N)
        plot(results[:, 2] / R, results[:, -1], 'o-')
        plot(r, p_r, 'r', linewidth=2)
        xlabel('Radius (as a fraction of cell tower radius)')
        ylabel('Probability of connectedness')
        ylim([0, 1.1])
        show()
if __name__ == "__main__":
    prob_1()
    prob_2()
    prob_3()
    part_4a()
    part_4b()

p2.pickle

cnumpy.core.multiarray
_reconstruct
p0
(cnumpy
ndarray
p1
(I0
tp2
S'b'
p3
tp4
Rp5
(I1
(I25
tp6
cnumpy
dtype
p7
(S'f8'
p8
I0
I1
tp9
Rp10
(I3
S'<'
p11
NNNI-1
I-1
I0
tp12
bI00
S'\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\xfc\xa9\xf1\xd2Mb\x80?\xfc\xa9\xf1\xd2Mb\x90?\x9c\xc4 \xb0rh\xa1?\x19\x04V\x0e-\xb2\xad?\xaa\xf1\xd2Mb\x10\xb8?L7\x89A`\xe5\xc0?#\xdb\xf9~j\xbc\xc4?\x19\x04V\x0e-\xb2\xcd?\xe3\xa5\x9b\xc4 \xb0\xd2?b\x10X9\xb4\xc8\xd6?\x8d\x97n\x12\x83\xc0\xda?\xac\x1cZd;\xdf\xdf?H\xe1z\x14\xaeG\xe1?%\x06\x81\x95C\x8b\xe4?\x19\x04V\x0e-\xb2\xe5?V\x0e-\xb2\x9d\xef\xe7?\xf8S\xe3\xa5\x9b\xc4\xe8?B`\xe5\xd0"\xdb\xe9?\xdb\xf9~j\xbct\xeb?\xcd\xcc\xcc\xcc\xcc\xcc\xec?\x1f\x85\xebQ\xb8\x1e\xed?j\xbct\x93\x18\x04\xee?\x17\xd9\xce\xf7S\xe3\xed?'
p13
tp14
b.

prop.pickle

cnumpy.core.multiarray
_reconstruct
p0
(cnumpy
ndarray
p1
(I0
tp2
S'b'
p3
tp4
Rp5
(I1
(I25
tp6
cnumpy
dtype
p7
(S'f8'
p8
I0
I1
tp9
Rp10
(I3
S'<'
p11
NNNI-1
I-1
I0
tp12
bI00
S'\x00\x00\x00\x00\x00\x00\x00\x00UUUUUU\xa5?UUUUUU\xb5?\x00\x00\x00\x00\x00\x00\xc0?UUUUUU\xc5?\xaa\xaa\xaa\xaa\xaa\xaa\xca?\x00\x00\x00\x00\x00\x00\xd0?\xaa\xaa\xaa\xaa\xaa\xaa\xd2?UUUUUU\xd5?\x00\x00\x00\x00\x00\x00\xd8?\xaa\xaa\xaa\xaa\xaa\xaa\xda?UUUUUU\xdd?\x00\x00\x00\x00\x00\x00\xe0?UUUUUU\xe1?\xaa\xaa\xaa\xaa\xaa\xaa\xe2?\x00\x00\x00\x00\x00\x00\xe4?UUUUUU\xe5?\xaa\xaa\xaa\xaa\xaa\xaa\xe6?\x00\x00\x00\x00\x00\x00\xe8?UUUUUU\xe9?\xaa\xaa\xaa\xaa\xaa\xaa\xea?\x00\x00\x00\x00\x00\x00\xec?UUUUUU\xed?\xaa\xaa\xaa\xaa\xaa\xaa\xee?\x00\x00\x00\x00\x00\x00\xf0?'
p13
tp14
b.