# let's start with a zero matrix
A = [
    [0,0,0],
    [0,0,0],
    [0,0,0]
]

# we will randomly pick the entries and then make normalization for each column
# it will be easier to iteratively construct the columns
# you may notice that each column is a probabilistic state
from random import randrange

normalization_factor = [0,0,0] # the normalization factor of each column may be different

for j in range(3): # each column is iteratively constructed
    normalization_factor[j] = 0
    while normalization_factor[j]==0: # the normalization factor cannot be zero
        for i in range(3):
            A[i][j] = randrange(101) # pick a random value between 0 and 100
            normalization_factor[j] += A[i][j]

# let's print matrix A before the normalization
# the entries are between 0 and 100
print("matrix A before normalization:")
for i in range(3):
    print(A[i])


# let's normalize each column
for j in range(3):
    for i in range(3):
        A[i][j] /= normalization_factor[j]  # shorter form of A[i][j] = A[i][j] / normalization_factor[j]
        
# let's print matrix A after the normalization
print() # print an empty line
print("matrix A after normalization:")
for i in range(3):
    print(A[i]) 

print()
print("the column summations must be 1")

sum = [0,0,0]
for j in range(3):
    for i in range(3):
        sum[j] += A[i][j]

print(sum) 

# Asja's probabilistic operator (coin flip protocol)
A = [
    [0.6,0.3],
    [0.4,0.7]
]

# one-step evolution of Asja's probabilistic operator on a given probabilistic state
def asja(prelist):
    newlist=[0,0]
    for i in range(2): # for each row
        for j in range(2): # for each column
            newlist[i] = newlist[i] + prelist[j] * A[i][j] # summation of pairwise multiplication 
    return newlist # return the new state
            
# initial state            
state = [1,0]  

# after one step
state = asja(state)
print("after one step, the state is",state)

# the new state is correct

# let's check one more step
state = asja(state)
print("after two steps, the state is",state)

# this is also correct
#
# then, let's evolve the system for more steps
for i in [3,6,9,12,24,48,96]:
    # start from the initial state
    state = [1,0]
    for t in range(i): # apply asja t times
        state = asja(state)
    # print the result
    print(state,"after",(t+1),"steps")

def evolve(Op,state):
    newstate=[]
    for i in range(len(Op)): # for each row
        # we calculate the corresponding entry of the new state
        newstate.append(0) # we set this value to 0 for the initialization
        for j in range(len(state)): # for each element in state
            newstate[i] = newstate[i] + Op[i][j] * state[j] # summation of pairwise multiplications
    return newstate # return the new probabilistic state

# test the function

# operator for the test
A = [
    [0.4,0.6,0],
    [0.2,0.1,0.7],
    [0.4,0.3,0.3]
]

# state for test
v = [0.1,0.3,0.6]

newstate = evolve(A,v)
print(newstate)

for step in [5,10,20,40]:
    new_state = [0.1,0.3,0.6] # initial state
    for i in range(step):
        new_state = evolve(A,new_state)
    print(new_state)

# change the initial state

for step in [5,10,20,40]:
    new_state = [1,0,0] # initial state
    for i in range(step):
        new_state = evolve(A,new_state)
    print(new_state)

# for random number generation
from random import randrange

# we will use evolve function
def evolve(Op,state):
    newstate=[]
    for i in range(len(Op)): # for each row
        newstate.append(0)
        for j in range(len(state)): # for each element in state
            newstate[i] = newstate[i] + Op[i][j] * state[j] # summation of pairwise multiplications
    return newstate # return the new probabilistic state

# the initial state
state = [0.5, 0, 0.5, 0]

# probabilistic operator for symbol a
A = [
    [0.5,  0, 0, 0],
    [0.25, 1, 0, 0],
    [0,    0, 1, 0],
    [0.25, 0, 0, 1]
]

# probabilistic operator for symbol b
B = [
    [1, 0, 0,    0],
    [0, 1, 0.25, 0],
    [0, 0, 0.5,  0],
    [0, 0, 0.25, 1]
]

#
# your solution is here
#

length = 40
total = 50
# total = 1000 # we will also test our code for 1000 strings

# we will check 5 cases
# let's use a list 
cases = [0,0,0,0,0]

for i in range(total): # total number of strings
    Na = 0
    Nb = 0
    string = ""
    state = [0.5, 0, 0.5, 0]
    for j in range(length): # generate random string
        if randrange(2) == 0: 
            Na = Na + 1 # new symbol is a
            string = string + "a"
            state = evolve(A,state) # update the probabilistic state by A
        else:
            Nb = Nb + 1 # new symbol is b
            string = string + "b"
            state = evolve(B,state) # update the probabilistic state by B
    # now we have the final state
    p0 = state[0] + state[1] # the probabilities of being in 00 and 01
    p1 = state[2] + state[3] # the probabilities of being in 10 and 11
    #print() # print an empty line
    print("(Na-Nb) is",Na-Nb,"and","(p0-p1) is",p0-p1)
    # let's check possible different cases
    
    # start with the case in which both are nonzero
    # then their multiplication is nonzero
    # let's check the sign of their multiplication
    if (Na-Nb) * (p0-p1) < 0: 
        print("they have opposite sign")
        cases[0] = cases[0] + 1
    elif (Na-Nb) * (p0-p1) > 0: 
        print("they have the same sign")
        cases[1] = cases[1] + 1
        
    #  one of them should be zero
    elif (Na-Nb) == 0:
        if (p0-p1) == 0: 
            print("both are zero")
            cases[2] = cases[2] + 1
        else: 
            print("(Na-Nb) is zero, but (p0-p1) is nonzero")
            cases[3] = cases[3] + 1
    elif (p0-p1) == 0: 
        print("(Na-Nb) is nonzero, while (p0-p1) is zero")
        cases[4] = cases[4] + 1
        
# check the case(s) that are observed and the case(s) that are not observed
print() # print an empty line
print(cases)        

