#
# your solution is here
#
# 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) 

# operator B
B = [
    [0.4,0.6,0],
    [0.2,0.1,0.7],
    [0.4,0.3,0.3]
]

# the current state
v = [0.1,0.3,0.6]

#
# your solution is here
#
# operator B
B = [
    [0.4,0.6,0],
    [0.2,0.1,0.7],
    [0.4,0.3,0.3]
]

# the current state
v = [0.1,0.3,0.6]

newstate = []

index = 0

for row in B:
    newstate.append(0)
    for i in range(len(row)):
        newstate[index] = newstate[index] + row[i] * v[i]
    index = index + 1
    
print(newstate)

# %%writefile linear_evolve.py

def linear_evolve(operator,state):
    newstate=[]
    for i in range(len(operator)): # 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] + operator[i][j] * state[j] # summation of pairwise multiplications
    return newstate # return the new probabilistic state

# test the function

# operator for the test
B = [
    [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 = linear_evolve(B,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 = linear_evolve(B,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 = linear_evolve(B,new_state)
    print(new_state)