# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

import numpy as np
import matplotlib.pyplot as plt
import time
# #%matplotlib inline

# +
# %matplotlib notebook

import matplotlib

#modify some matplotlib parameters to manage the images for illustrator
matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.rcParams['ps.fonttype'] = 42
# -

# create a initial grid
grid_size = 500    #evetually define as an input
grid = np.zeros((grid_size,grid_size))

# define the initial pattern
grid[int(grid_size/2), int(grid_size/2)]=2
grid


def show_grid(grid_array):
    plt.figure()
    plt.imshow(grid, cmap=plt.cm.gray)
    plt.show()
    


show_grid(grid)


# +
# Define rules

def rule1(grid):
    # 2 = growing cell
    # 1 = stationary cell
    # 0 = empty space
    
    #this rule makes each cell divide only one time per step
    new_grid = grid   #to separate the evaluation from the actualization
    #g_ones = grid[grid == 1] 
    g_index =  np.nonzero(grid == 2) # growth index = where cell value == 2
    #g_mask = np.array([[0,1,2],[3,8,4],[5,6,7]]) #position 8 is the center and not selectable
    
    #shuffle the elements in order to eliminate the bias of the computation order
    
    arr = np.arange(g_index[0].shape[0])
    np.random.shuffle(arr)
    
    #for i in range(len(g_index[0])):  #go through every position with a "growing cell" ---> value = 2
    for i in arr:  #go through every position with a "growing cell" ---> value = 2
        
        mask_pos = np.array([[-1,-1],[0,-1],[1,-1],[-1,0],[1,0],[-1,1],[0,1],[1,1]])
        remove_list = []  #initialize a list to store "neighbours spaces"( = mask_pos) ocuped by cells
        
        # go through every position surrounding position --> eigh posibilities
        
        for j in range(len(mask_pos)):

            m = g_index[0][i] + mask_pos[j][0]
            n = g_index[1][i] + mask_pos[j][1]

            if grid[m,n] !=0 :  #make a list with the positions which are not empty places
                
                remove_list.append(j)

        new_mask_pos = np.delete(mask_pos, remove_list, 0)        

        # to exit when there is not a surrounding empty position
        l = len(new_mask_pos)
        if l > 1:
            r_pos = np.random.randint(l)   # a random number between [0,len[
            new_pos = new_mask_pos[r_pos]

            #m = g_ones_index[0][i] + mask_pos[new_pos][0]
            #n = g_ones_index[1][i] + mask_pos[new_pos][1]
            m = g_index[0][i] + new_pos[0]
            n = g_index[1][i] + new_pos[1]
            grid[m,n] = 2

        elif l == 1:
            
            new_pos = new_mask_pos[0]

            m = g_index[0][i] + new_pos[0]
            n = g_index[1][i] + new_pos[1]
            
            grid[m,n] = 2

        else: #when len(new_mask_pos) == 0
            m = g_index[0][i]
            n = g_index[1][i]
            grid[m,n] = 1     # then, that position will not be evaluated again
            
    return(grid)


# +
# Define rules

def rule2(grid):
    # 2 = growing cell
    # 1 = stationary cell
    # 0 = empty space
    
    #this rule makes each cell divide only one time per step
    #g_ones = grid[grid == 1] 
    g_index =  np.nonzero(grid == 2) # growth index = where cell value == 2
    #g_mask = np.array([[0,1,2],[3,8,4],[5,6,7]]) #position 8 is the center and not selectable
    
    #choose a random cell in the dividing state
    cell_pos = int(np.random.rand(1)[0]*g_index[0].shape[0])
    
    mask_pos = np.array([[-1,-1],[0,-1],[1,-1],[-1,0],[1,0],[-1,1],[0,1],[1,1]])
    remove_list = []  #initialize a list to store "neighbours spaces"( = mask_pos) ocuped by cells

    # go through every position surrounding position --> eigh posibilities

    for j in range(len(mask_pos)):

        m = g_index[0][cell_pos] + mask_pos[j][0]
        n = g_index[1][cell_pos] + mask_pos[j][1]

        if grid[m,n] !=0 :  #make a list with the positions which are not empty places

            remove_list.append(j)

    new_mask_pos = np.delete(mask_pos, remove_list, 0)        

    # to exit when there is not a surrounding empty position
    l = len(new_mask_pos)
    if l > 1:
        r_pos = np.random.randint(l)   # a random number between [0,len[
        new_pos = new_mask_pos[r_pos]

        #m = g_ones_index[0][i] + mask_pos[new_pos][0]
        #n = g_ones_index[1][i] + mask_pos[new_pos][1]
        m = g_index[0][cell_pos] + new_pos[0]
        n = g_index[1][cell_pos] + new_pos[1]
        grid[m,n] = 2

    elif l == 1:

        new_pos = new_mask_pos[0]

        m = g_index[0][cell_pos] + new_pos[0]
        n = g_index[1][cell_pos] + new_pos[1]

        grid[m,n] = 2

    else: #when len(new_mask_pos) == 0
        m = g_index[0][cell_pos]
        n = g_index[1][cell_pos]
        grid[m,n] = 1     # then, that position will not be evaluated again
            
    return(grid)


# +
# Define rules

def rule3(grid):
    # 2 = growing cell
    # 1 = stationary cell
    # 0 = empty space
    
    #this rule makes each cell divide only one time per step
    #g_ones = grid[grid == 1] 
    g_index =  np.nonzero(grid == 2) # growth index = where cell value == 2
    #g_mask = np.array([[0,1,2],[3,8,4],[5,6,7]]) #position 8 is the center and not selectable
    
    #choose a random cell in the dividing state
    index_pos = int(np.random.rand(1)[0]*g_index[0].shape[0])
    m = g_index[0][index_pos]
    n = g_index[1][index_pos]
    #define the neighborhood
    nb = grid[m-1:m+2,n-1:n+2]  #nb = neighborhood
    
    
    #define the free spaces in nb
    fs = np.where(nb == 0)
    
    if fs[0].shape[0] > 0: # go if there are a place in the neighbour

        if len(fs[0]) == 1:
            grid[m,n] = 1     # then, that position will not be evaluated again

         #grown over an empty position
        new_pos = int(np.random.rand(1)[0]*fs[0].shape[0]) #new pos in the neighbour matrix
        m_new = m + fs[0][new_pos] - 1  #-1 to convert [0 1 2] to [-1 0 1]
        n_new = n + fs[1][new_pos] - 1
        grid[m_new, n_new] = 2
    
        #save the cell grid index positions
        cell_index = [m,n]
        ncell_index = [m_new, n_new]
    else:
        grid[m,n] = 1
    
    return(grid)

# +
#to show rule 3
time1 = time.clock()

# create a initial grid
grid_size = 150    #evetually define as an input
grid = np.zeros((grid_size,grid_size))

# define the initial pattern
#grid[int(grid_size/2), int(grid_size/2)]=2
grid = initial_pattern(grid, 0)

#perform the loop
steps = 10000
#sleep_time = 0.01

#to save the last figure
filename = 'Segregation\\image_%03d.jpg'
#filename = os.path.join(fpath, 'image_%05d.jpg')

fig = plt.figure()
fig.show()
show_grid(grid)  #show the initial grid
plt.title('step 0')

every = 100
count = 0

for i  in range(steps):
    #time.sleep(sleep_time)    
    grid = rule3(grid)
    #plt.imshow(grid, cmap=plt.cm.gray)
    #plt.title('step '+ str(i+1))
    #fig.canvas.draw()
    
    if i%every == 0 or i == steps-1:
        count += 1
        plt.title('step '+ str(i+1))
        plt.imshow(grid, cmap=plt.cm.gray)
        #plt.savefig(str(filename) + ".pdf", transparent=True)
        plt.savefig(filename%(count))#, transparent=True)
#plt.imshow(grid, cmap=plt.cm.gray)

elapsed = time.clock() - time1

print(elapsed)  


# -

def select_cell(grid):
    # 2 = growing cell
    # 1 = stationary cell
    # 0 = empty space
    
    #this rule makes each cell divide only one time per step
    g_index =  np.nonzero(grid == 2) # growth index = where cell value == 2
    
    #choose a random cell in the dividing state
    index_pos = int(np.random.rand(1)[0]*g_index[0].shape[0])
    m = g_index[0][index_pos]
    n = g_index[1][index_pos]

    
    #save the cell grid index positions
    cell_index = [m,n]
        
    return(cell_index)


def check_nbhd(grid, cell_index):
    #chek free spaces in the neighbourhood
    # fs: array
    #    index of free spaces in the neighborhood
    m = cell_index[0]
    n = cell_index[1]
    
    #define the neighborhood
    nb = grid[m-1:m+2,n-1:n+2]  #nb = neighborhood
    
    #define the free spaces in nb
    fs = np.where(nb == 0)
    
    return(fs)


def nb_prob(grid, cell_index, prob_dist = 'contact_linear'):
    # assign division probabilities based on empty space cell contacts
    
    # prob_dist: uniform - contact_linear - contact_exp
    #    contact linear is the default or if another thing is written
    
    # fs: array
    #    index of free spaces in the neighborhood
    # return
    # prob: list
    #    list with the [0,1] probability partition limit of each free space
    #    e.g. prob = [0.23, 0.81, 1] --> second cell has bigger probability
    
    m = cell_index[0]
    n = cell_index[1]
    
    #define neighborhood
    nb = grid[m-1:m+2,n-1:n+2]  #nb = neighborhood
    
    #define the free spaces in nb
    fs = np.where(nb == 0)
    #define cell spaces in bn
    cs = np.where(nb != 0)
    
    fs_num = len(fs[0])
    prob = np.zeros(fs_num)
    contacts = np.zeros(fs_num)
    
    if prob_dist != 'uniform':
        
        # if prob_dist is something different from the options, contact_linear is the default
        for i in range(fs_num):
            mg = m + fs[0][i] - 1  #-1 to convert [0 1 2] to [-1 0 1]
            ng = n + fs[1][i] - 1

            i_nb = grid[mg-1:mg+2,ng-1:ng+2] # i position neighborhood

            occup = np.where(i_nb != 0)
            contacts[i] = len(occup[0])  #save the number of contacts of this position

        if prob_dist == 'contact_exp':
            contacts = np.exp(contacts)
            
    else:
        contacts = np.ones(fs_num)  #assign uniform values
        
    total = sum(contacts)

    prob[0] = (contacts[0]/total)

    for i in range(1,fs_num):
        prob[i] = prob[i-1]+contacts[i]/total
        
    return(prob)



def cell_divition_uniform(grid, cell_index, fs):
    # uniform neighborhood divition probability
    #fs: free neighborhood spaces
    
    m = cell_index[0]
    n = cell_index[1]
    
    if len(fs[0]) == 1:
        grid[m,n] = 1     # then, that position will not divide again

    #grown over an empty position
    #new_pos = int(np.random.rand(1)[0]*fs[0].shape[0])
    new_pos = int(np.random.rand(1)[0]*fs[0].shape[0]) #new pos in the neighbour matrix
    m_new = m + fs[0][new_pos] - 1  #-1 to convert [0 1 2] to [-1 0 1]
    n_new = n + fs[1][new_pos] - 1
    grid[m_new, n_new] = 2     # crates the new cell
        
    ncell_index = [m_new ,n_new]
    
    return(grid, ncell_index)


def cell_divition(grid, cell_index, fs, fs_proba):
    #fs: free neighborhood spaces
    #fs_proba: free spaces growth probabilities
    
    m = cell_index[0]
    n = cell_index[1]
    
    if len(fs[0]) == 1:
        grid[m,n] = 1     # then, that position will not divide again

    #grown over an empty position

    rand_val = np.random.rand(1)[0]
    # find the first position which is bigger than rand_val
    new_pos = np.where( (fs_proba > rand_val) == True )[0][0]  #new pos in the neighbour matrix
   
    m_new = m + fs[0][new_pos] - 1  #-1 to convert [0 1 2] to [-1 0 1]
    n_new = n + fs[1][new_pos] - 1
    grid[m_new, n_new] = 2     # crates the new cell
        
    ncell_index = [m_new ,n_new]
    
    return(grid, ncell_index)


plt.figure()
im_grid = np.zeros((100,100,4))
im_grid[:,:,0] = np.ones((100,100))*0
im_grid[:,:,1] = np.ones((100,100))*80
im_grid[:,:,2] = np.ones((100,100))*0
im_grid[:,:,3] = np.ones((100,100))*1
plt.imshow(im_grid)
plt.figure()
im_grid[:,:,1] = np.ones((100,100))*1
plt.imshow(im_grid)


def initial_plasmids(grid, pattern_num = 0, num_plas = 2, max_copy = 4):
    # grid: initial grid
    c_index =  np.nonzero(grid) # c_index, 
    #cell_number = c_index[0].shape[0]
    
    gs = grid.shape
    pattern = np.zeros((gs[0],gs[1],max_copy)) #initialize the pattern array
    
    # add different patterns
    if pattern_num == 0:  #random plasmid pattern
        for i in range(c_index[0].shape[0]): #assign a random plasmid pattern to each cell position
            pattern[c_index[0][i],c_index[1][i],:] = ((num_plas +1 )*np.random.rand(max_copy)).astype(int) 
                                                            #num_plas +1 to add "no-plasmid" state
            
    elif pattern_num == 1:
        pattern = np.ones((grid.shape))
        
    return(pattern)
    


def role_divideFlag(plasmids):
    #plasmids: cell plasmids vector
    
    max_plasmids = plasmids.shape[0]
    num_plasmids = np.nonzero(plasmids)[0].shape[0]
    divisor = max_plasmids*1.1  #arbitrary defined to make division(max_plasmids number) < 1
       
    # make a cuadratic function of probabilities
    probability = (num_plasmids/divisor)**2 
    #if a cell has no plasmids --> will not divide
    
    if np.random.rand(1) < probability:
        return(1) # divide
    else:
        return(0) # not divide
    
    #Probability tables
    #plasmid_nums = np.arange(max_plasmids +1)
    #probability = (plasmid_nums/divisor)**2 


def create_image(grid, plasgrid):
    im_s = plasgrid.shape
    aux_imR = np.zeros((im_s[0],im_s[1],im_s[2]))
    aux_imG = np.zeros((im_s[0],im_s[1],im_s[2]))
    
    for i in range(im_s[2]):
        aux_imR[:,:,i] = 1*(plasgrid[:,:,i]==1)
        aux_imG[:,:,i] = 1*(plasgrid[:,:,i]==2) 
    
    aux_imR = np.sum(aux_imR,axis=2)
    aux_imG = np.sum(aux_imG,axis=2)
    aux_transparency = 0.5*(grid[:,:]==1) + 1*(grid[:,:]==2)
    
    # create the image
    im_grid = np.zeros((im_s[0],im_s[1],im_s[2]))
    im_grid[:,:,0] = np.multiply(np.ones((im_s[0],im_s[1])),aux_imR)
    im_grid[:,:,1] = np.multiply(np.ones((im_s[0],im_s[1])),aux_imG)
    #im_grid[:,:,2] = np.ones((100,100))*250
    im_grid[:,:,3] = np.multiply(np.ones((im_s[0],im_s[1])),aux_transparency)
     # stationary cell -> transparency = 0.5)
    
    return(im_grid)


def create_image2(grid, plasgrid):
    im_s = plasgrid.shape
    aux_imR = np.zeros((im_s[0],im_s[1],im_s[2]))
    aux_imG = np.zeros((im_s[0],im_s[1],im_s[2]))
    
    for i in range(im_s[2]):
        aux_imR[:,:,i] = 1*(plasgrid[:,:,i]==1)
        aux_imG[:,:,i] = 1*(plasgrid[:,:,i]==2) 
    
    aux_imR = np.multiply(1*(np.sum(aux_imR,axis=2)>0),50*(grid[:,:]==1)) + 1*(np.sum(aux_imR,axis=2)>0)
    aux_imG = np.multiply(1*(np.sum(aux_imG,axis=2)>0),50*(grid[:,:]==1)) + 1*(np.sum(aux_imG,axis=2)>0)
    
    # create the image
    im_grid = np.zeros((im_s[0],im_s[1],3))
    im_grid[:,:,0] = np.multiply(np.ones((im_s[0],im_s[1])),aux_imR)
    im_grid[:,:,1] = np.multiply(np.ones((im_s[0],im_s[1])),aux_imG)
    #im_grid[:,:,2] = np.ones((100,100))*250
    
    return(im_grid)


def plasmid_gProb(g_ratio=[1,1], p_types = [1,2]):
    #define a growtn probability (=velocity) based on the plasmids
    #g_ratio: ratio of growth rate between genotypes (i.e. plasmids)
    #p_types: plasmids types or labels
    
    #built the probability class vector
    cat_len = len(g_ratio)
    probs = np.zeros(cat_len)
    
    denominator = sum(g_ratio)

    probs[0] = g_ratio[0]/denominator

    for i in range(1,cat_len):
        probs[i] = probs[i-1]+g_ratio[i]/denominator
    
    return(probs)


def plasm_g_test(plasmids,probs):
        #perform the probability test
    rand_val = np.random.rand(1)
    pos = np.where( (probs > rand_val) == True )[0][0]
    
    ptype = pos + 1
    found = np.where(plasmids == ptype)[0]
    growth = False
    
    if found.size>0:
        growth = True
    return(growth)


def cell_ratio(plasmgrid, ptype = [1,2]):
    
    c_num_plasm = np.sum(plasmgrid>0, axis=2) #number of plasmids in each grid

    plasm_sum = np.sum(plasmgrid, axis = 2)
    divition = np.divide(plasm_sum,c_num_plasm)
    
    #total = np.sum(np.isnan(divition) == False, axis = (0,1)) #it include cells with mix plasmids
    found = np.zeros(len(ptype))
    total = 0
    for i in range(len(ptype)):
        found[i] = len(np.where(divition == ptype[i])[0])
        total += found[i]
    
    ratio = found[0]/total
        
    return(ratio)


count=0
plasmids = np.ones(4)*2
for i in range(1000):
    plas_probs = plasmid_gProb(g_ratio= [1,2])
    ifG = plasm_g_test(plasmids, plas_probs)
    if ifG == True:
        count+=1
print(count)

#main
sim_num = 1
all_ratios = []
for j in range(sim_num):
    time1 = time.clock()
    # create a initial empty grid
    grid_size = 1000    #evetually define as an input
    grid = np.zeros((grid_size,grid_size))

    # define the initial grid and plasmid pattern
    grid = initial_pattern(grid, 1)
#    show_grid(grid)
    plasm_grid = initial_plasmids(grid)

    # Show the initial state
#    im_grid = create_image(grid, plasm_grid)
#    plt.imshow(im_grid)

    #perform the loop
    steps = 100000
    #sleep_time = 0.01

    #to save the last figure
    #filename = 'null'
#    filename = 'Segregation\\ratios\\image_%03d.jpg'
    filename = 'Seg_1ratio.jpg'
#    fig = plt.figure()
#    plt.title('step 0')

    #this two lines to save secuential images
    every = 100   #every how many save images
    count = 0

    #define plasmid grwoth ratio
    plas_probs = plasmid_gProb(g_ratio= [1,1])
    # g_ratio = [2,1] --> plasmid 1 divide twice fast than plasmid 2
    ratios = []  #to store the cell type ratios
    for i  in range(steps):

        #select a random growing cell
        cell_pos = select_cell(grid)
        free_nb = check_nbhd(grid, cell_pos)

        if free_nb[0].shape[0] > 0: # go if there is a place in the neighborhood

            plasmids = plasm_grid[cell_pos[0], cell_pos[1],:] #get its plasmids
            c_growth = plasm_g_test(plasmids, plas_probs)

            if c_growth == True: 

                # tal vez poner esto arriba y chequear que los plasmidos no estene en cero

                #update its plasmids and cell state, n:new
                n_plasmids, n_state = plasmid_update(plasmids, cell_pos)

                plasm_grid[cell_pos[0], cell_pos[1],:] = n_plasmids
                grid[cell_pos[0], cell_pos[1]] = n_state
                 #state will not be evaluated before role_divide
                 #role_divide function shouldn´t allow divition of that cell

                divide_flag = role_divideFlag(n_plasmids)

                #perform the divition if flag changed
                if divide_flag != 0:
                    #assign a cell to a new position
                    free_proba = nb_prob(grid, cell_pos, prob_dist = 'contact_exp')
                    grid, nCell_pos = cell_divition(grid, cell_pos, free_nb, free_proba) 

                    #split the mother plasmids
                    m_plasmids, c_plasmids = divide_plasmids(n_plasmids)

                    #assign mother and child plasmids
                    plasm_grid[cell_pos[0], cell_pos[1],:] = m_plasmids
                    plasm_grid[nCell_pos[0], nCell_pos[1],:] = c_plasmids
        else:
            grid[cell_pos[0],cell_pos[1]] = 1

        #save cell type ratios
        if i%every == 0:
            ratios.append(cell_ratio(plasm_grid))
        
        #Plot the result
        if i == steps-1:
        #if i%every == 0 or i == steps-1:

        #    count += 1

            plt.title('step '+ str(i+1))
            im_grid = create_image2(grid, plasm_grid)
            plt.imshow(im_grid)
        #    #fig.canvas.draw()
        #    #plt.savefig(str(filename) + ".pdf", transparent=True)
            #plt.savefig(filename%(count), transparent=True)
        #    plt.savefig(filename%(j), transparent=True)
            plt.savefig(filename, transparent=True)
            
    all_ratios.append(np.asarray(ratios))
    
    elapsed = time.clock() - time1
    print(elapsed)


mean_ratio = 0
plt.figure()
for i in range(len(all_ratios)):
    plt.plot(all_ratios[i])
    mean_ratio += all_ratios[i][-1]
plt.show()
mean_ratio= mean_ratio/len(all_ratios)
print(mean_ratio)

ratio11=all_ratios
ratio11_mean = mean_ratio

# +
plt.figure()
for i in range(len(ratio11)):
    plt.plot(ratio11[i])

plt.title('growth ratio 1:1')
plt.ylabel('cell type ratio')
plt.xlabel('check point step number')
plt.show()
plt.savefig('ratio11', transparent=True)
# -

ratio32=all_ratios
ratio32_mean = mean_ratio

ratio43=all_ratios
ratio43_mean = mean_ratio

# +
plt.figure()
for i in range(len(ratio43)):
    plt.plot(ratio43[i])

plt.title('growth ratio 4:3')
plt.ylabel('cell type ratio')
plt.xlabel('check point step number')
plt.show()
plt.savefig('ratio43', transparent=True)
# -

ratio109=all_ratios
ratio109_mean = mean_ratio

# +
mean_ratios = [ratio11_mean,ratio109_mean,ratio43_mean,ratio32_mean]
expected = [1/2, 10/19, 4/7, 3/5]
plt.figure()

plt.plot([1, 2, 3, 4], mean_ratios, 'bo', label= 'observed ratio')
plt.plot([1, 2, 3, 4], expected, 'ro', label = 'growth ratio')
plt.ylabel('cell type ratio')
plt.legend()
plt.xticks([])
#plt.xticks(['1:1','10:9','4:3','3:2'])
plt.xticks([1,2,3,4], ['1:1','10:9','4:3','3:2'])
plt.show()

plt.savefig('ratios_obs_exp', transparent=True)

# +
mean_ratios = [ratio11_mean,ratio109_mean,ratio43_mean,ratio32_mean]
expected = [1/2, 10/19, 4/7, 3/5]
plt.figure()

plt.plot(expected, mean_ratios, 'bo')
plt.xlabel
plt.show()


# -

def plasmid_update(plasmids, pos_index):
    #plasmids: vector with plasmids. e.g [0,1,1,0,2]
    state = 2 # cell state = growing state
    
    plasmids_pos = np.nonzero(plasmids)
    empty_pos = np.where(plasmids == 0)
    
    num_plas = plasmids_pos[0].shape[0]
    
    if num_plas == 0:
        #it means no plasmid in the cell
        state = 1   #to not evaluate this cell in the loop again
        
    elif num_plas == plasmids.shape[0]:
        #it means all plasmids positions are full
        return(plasmids, state)
    
    else:
        copied_pos = np.random.randint(num_plas)
        plasmids[empty_pos[0][0]] = plasmids[plasmids_pos[0][copied_pos]]
        #copy the plasmid in the first free space
    
    return(plasmids, state)   
    


def divide_plasmids(plasmids):
    #plasmids: cell plasmids
    p_size = plasmids.size
    
    mother_p = np.zeros(p_size)
    child_p = np.zeros(p_size)
    
    np.random.shuffle(plasmids) #shuffle the plasmids
    
    if (p_size & 1) == 1: #odd case
        
        #sum a random value to choose which cell keep more plasmids
        rand_val = np.random.rand(1)
        half_p = int(p_size/2 + rand_val) 
        
    else:  #even case
        half_p = int(p_size/2)
    
    mother_p[:half_p] = plasmids[:half_p]
    child_p[half_p:]= plasmids[half_p:] 
    
    return(mother_p, child_p)
    


def initial_pattern(grid, pattern_num):
    
    pattern = {}  #initiate initial pattern dictionary
    
    # add different patterns
    pattern[0] = np.array([[2]])
    pattern[1] = np.array([[0, 0, 2, 0, 0],[0,2,2,2,0],[2,2,1,2,2],[0,2,2,2,0],[0,0,2,0,0]])
    pattern[2] = np.ones((2,35))*2
    
    #make elements which are not in the border to be = 1
    fixed_pat = pattern[pattern_num]
    
    
    #put the pattern in the grid
    gs = grid.shape
    
    m0 = int(gs[0]/2)
    n0 = int(gs[1]/2)
    
    ps = fixed_pat.shape
    mpm = int(ps[0]/2)
    npm = int(ps[1]/2)
            
    for i in range(ps[0]):
        for j in range(ps[1]):
            m = m0 + (i - mpm)
            n = n0 + (j - npm)
            
            grid[m,n] = fixed_pat[i,j]


    return(grid)
    

# +
#perform the loop
steps = 50
sleep_time = 0.1

fig = plt.figure()
fig.show()
show_grid(grid)  #show the initial grid
plt.title('step 0')

for i  in range(steps):
    time.sleep(sleep_time)    
    grid = rule1(grid)
    plt.imshow(grid, cmap=plt.cm.gray)
    plt.title('step '+ str(i+1))
    fig.canvas.draw()
  
# -

plt.figure()
plt.imshow(grid)

# +
#perform the loop
steps = 10
sleep_time = 0.5

show_grid(grid)  #show the initial grid
plt.figure(0)
for i  in range(steps):
    time.sleep(sleep_time)    
    grid = rule1(grid)
    show_grid(grid)
    

# +
#to show
time1 = time.clock()

# create a initial grid
grid_size = 150    #evetually define as an input
grid = np.zeros((grid_size,grid_size))

# define the initial pattern
#grid[int(grid_size/2), int(grid_size/2)]=2
grid = initial_pattern(grid, 2)

#perform the loop
steps = 10000
#sleep_time = 0.01

#to save the last figure
filename = 'null'

fig = plt.figure()
fig.show()
show_grid(grid)  #show the initial grid
plt.title('step 0')

for i  in range(steps):
    #time.sleep(sleep_time)    
    #grid = rule1(grid)
    grid = rule3(grid)
    #plt.imshow(grid, cmap=plt.cm.gray)
    plt.title('step '+ str(i+1))
    #fig.canvas.draw()
    
    if i == steps-1:
        plt.savefig(str(filename) + ".pdf", transparent=True)
plt.imshow(grid, cmap=plt.cm.gray)

elapsed = time.clock() - time1

print(elapsed)  

# +
#to show
time1 = time.clock()
# create a initial grid
grid_size = 150    #evetually define as an input
grid = np.zeros((grid_size,grid_size))

# define the initial pattern
grid = initial_pattern(grid, 0)

#perform the loop
steps = 15000
#sleep_time = 0.01

#to save the last figure
filename = 'null'

fig = plt.figure()
plt.title('step 0')

for i  in range(steps):
    #time.sleep(sleep_time)    
    #grid = rule1(grid)
    grid = rule2(grid)
    #plt.imshow(grid, cmap=plt.cm.gray)
    plt.title('step '+ str(i+1))
    
    
    if i == steps-1:
        
        plt.imshow(grid, cmap=plt.cm.gray)
        #fig.canvas.draw()
        plt.savefig(str(filename) + ".pdf", transparent=True)
        

elapsed = time.clock() - time1
print(elapsed)
# -

mask[1,2]

vals = np.array([[-1,-1],[0,-1],[1,-1],[-1,0],[1,0],[-1,1],[0,1],[1,1],[0,0]])
vals[4][1]

pos = np.nonzero(grid == 0)
print(pos[0][10])

len(pos[1])

grid[grid == 1]

# +
#make automata loop

time_steps = 100

for t in time_steps:
    
# -

np.random.randint(8) # a random number between [0,8[ --> eigh posibilities

grid[0,0]

# +
time1 = time.clock()
arr = np.random.rand(5000,5000)
elapsed = time.clock() - time1

print(elapsed)  # it takes 0.21186613101770035 seg for me



# -

np.amin(np.nonzero(grid == 2)[0])

# +
#hacer clases
# to not make:
# cell = [px,py,vx,vy]
# cell[0]+= cell[2]*deltaT


# make classes is the same in general terms (in computation time, etc), but is a more recomended paradigma
# because is more clear and organized.

class cell:
    def __init__(self,px,py): #primer parametro siempre es self -->que hace referencia al objeto cell
        self.px = px
        self.py = py
        self.vx = 0
        self.vy = 0
        
    # la gracia es que puedo definir funciones dentro de la clase    
    def mover(self,vx,vy,t):  #siempre el primero es self!
        self.px += vx*t
        self.py += vy*t
        
celula = cell(-1,1)

print(celula.px) #--> arroja -1
ceula.px +=1
print(celula.px) #--> arroja 0

celula.mover(-1,1,1)

lista = []
for n in range(10):
    lista.append(cell(n/10,n/10))
    
# or to make some computations faster when it is a big number of elements (>1000)

class sim:
    def __init__(self, ncells):
        positions = np.array((ncells,2))

#Also you can create an inheritance class
class bacterium(cell):   #that means bacterium "is a" cell
    def infect():   #then you can add new functions to the class bacterium
        
# todos los objetos en python pueden guardarse con el packete pickle
# -



StarcoderdataJupyterScriptsDedupFiltered