Celular automata tester.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "#%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "import matplotlib\n",
    "\n",
    "#modify some matplotlib parameters to manage the images for illustrator\n",
    "matplotlib.rcParams['pdf.fonttype'] = 42\n",
    "matplotlib.rcParams['ps.fonttype'] = 42"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a initial grid\n",
    "grid_size = 500    #evetually define as an input\n",
    "grid = np.zeros((grid_size,grid_size))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the initial pattern\n",
    "grid[int(grid_size/2), int(grid_size/2)]=2\n",
    "grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_grid(grid_array):\n",
    "    plt.figure()\n",
    "    plt.imshow(grid, cmap=plt.cm.gray)\n",
    "    plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_grid(grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define rules\n",
    "\n",
    "def rule1(grid):\n",
    "    # 2 = growing cell\n",
    "    # 1 = stationary cell\n",
    "    # 0 = empty space\n",
    "    \n",
    "    #this rule makes each cell divide only one time per step\n",
    "    new_grid = grid   #to separate the evaluation from the actualization\n",
    "    #g_ones = grid[grid == 1] \n",
    "    g_index =  np.nonzero(grid == 2) # growth index = where cell value == 2\n",
    "    #g_mask = np.array([[0,1,2],[3,8,4],[5,6,7]]) #position 8 is the center and not selectable\n",
    "    \n",
    "    #shuffle the elements in order to eliminate the bias of the computation order\n",
    "    \n",
    "    arr = np.arange(g_index[0].shape[0])\n",
    "    np.random.shuffle(arr)\n",
    "    \n",
    "    #for i in range(len(g_index[0])):  #go through every position with a \"growing cell\" ---> value = 2\n",
    "    for i in arr:  #go through every position with a \"growing cell\" ---> value = 2\n",
    "        \n",
    "        mask_pos = np.array([[-1,-1],[0,-1],[1,-1],[-1,0],[1,0],[-1,1],[0,1],[1,1]])\n",
    "        remove_list = []  #initialize a list to store \"neighbours spaces\"( = mask_pos) ocuped by cells\n",
    "        \n",
    "        # go through every position surrounding position --> eigh posibilities\n",
    "        \n",
    "        for j in range(len(mask_pos)):\n",
    "\n",
    "            m = g_index[0][i] + mask_pos[j][0]\n",
    "            n = g_index[1][i] + mask_pos[j][1]\n",
    "\n",
    "            if grid[m,n] !=0 :  #make a list with the positions which are not empty places\n",
    "                \n",
    "                remove_list.append(j)\n",
    "\n",
    "        new_mask_pos = np.delete(mask_pos, remove_list, 0)        \n",
    "\n",
    "        # to exit when there is not a surrounding empty position\n",
    "        l = len(new_mask_pos)\n",
    "        if l > 1:\n",
    "            r_pos = np.random.randint(l)   # a random number between [0,len[\n",
    "            new_pos = new_mask_pos[r_pos]\n",
    "\n",
    "            #m = g_ones_index[0][i] + mask_pos[new_pos][0]\n",
    "            #n = g_ones_index[1][i] + mask_pos[new_pos][1]\n",
    "            m = g_index[0][i] + new_pos[0]\n",
    "            n = g_index[1][i] + new_pos[1]\n",
    "            grid[m,n] = 2\n",
    "\n",
    "        elif l == 1:\n",
    "            \n",
    "            new_pos = new_mask_pos[0]\n",
    "\n",
    "            m = g_index[0][i] + new_pos[0]\n",
    "            n = g_index[1][i] + new_pos[1]\n",
    "            \n",
    "            grid[m,n] = 2\n",
    "\n",
    "        else: #when len(new_mask_pos) == 0\n",
    "            m = g_index[0][i]\n",
    "            n = g_index[1][i]\n",
    "            grid[m,n] = 1     # then, that position will not be evaluated again\n",
    "            \n",
    "    return(grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define rules\n",
    "\n",
    "def rule2(grid):\n",
    "    # 2 = growing cell\n",
    "    # 1 = stationary cell\n",
    "    # 0 = empty space\n",
    "    \n",
    "    #this rule makes each cell divide only one time per step\n",
    "    #g_ones = grid[grid == 1] \n",
    "    g_index =  np.nonzero(grid == 2) # growth index = where cell value == 2\n",
    "    #g_mask = np.array([[0,1,2],[3,8,4],[5,6,7]]) #position 8 is the center and not selectable\n",
    "    \n",
    "    #choose a random cell in the dividing state\n",
    "    cell_pos = int(np.random.rand(1)[0]*g_index[0].shape[0])\n",
    "    \n",
    "    mask_pos = np.array([[-1,-1],[0,-1],[1,-1],[-1,0],[1,0],[-1,1],[0,1],[1,1]])\n",
    "    remove_list = []  #initialize a list to store \"neighbours spaces\"( = mask_pos) ocuped by cells\n",
    "\n",
    "    # go through every position surrounding position --> eigh posibilities\n",
    "\n",
    "    for j in range(len(mask_pos)):\n",
    "\n",
    "        m = g_index[0][cell_pos] + mask_pos[j][0]\n",
    "        n = g_index[1][cell_pos] + mask_pos[j][1]\n",
    "\n",
    "        if grid[m,n] !=0 :  #make a list with the positions which are not empty places\n",
    "\n",
    "            remove_list.append(j)\n",
    "\n",
    "    new_mask_pos = np.delete(mask_pos, remove_list, 0)        \n",
    "\n",
    "    # to exit when there is not a surrounding empty position\n",
    "    l = len(new_mask_pos)\n",
    "    if l > 1:\n",
    "        r_pos = np.random.randint(l)   # a random number between [0,len[\n",
    "        new_pos = new_mask_pos[r_pos]\n",
    "\n",
    "        #m = g_ones_index[0][i] + mask_pos[new_pos][0]\n",
    "        #n = g_ones_index[1][i] + mask_pos[new_pos][1]\n",
    "        m = g_index[0][cell_pos] + new_pos[0]\n",
    "        n = g_index[1][cell_pos] + new_pos[1]\n",
    "        grid[m,n] = 2\n",
    "\n",
    "    elif l == 1:\n",
    "\n",
    "        new_pos = new_mask_pos[0]\n",
    "\n",
    "        m = g_index[0][cell_pos] + new_pos[0]\n",
    "        n = g_index[1][cell_pos] + new_pos[1]\n",
    "\n",
    "        grid[m,n] = 2\n",
    "\n",
    "    else: #when len(new_mask_pos) == 0\n",
    "        m = g_index[0][cell_pos]\n",
    "        n = g_index[1][cell_pos]\n",
    "        grid[m,n] = 1     # then, that position will not be evaluated again\n",
    "            \n",
    "    return(grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define rules\n",
    "\n",
    "def rule3(grid):\n",
    "    # 2 = growing cell\n",
    "    # 1 = stationary cell\n",
    "    # 0 = empty space\n",
    "    \n",
    "    #this rule makes each cell divide only one time per step\n",
    "    #g_ones = grid[grid == 1] \n",
    "    g_index =  np.nonzero(grid == 2) # growth index = where cell value == 2\n",
    "    #g_mask = np.array([[0,1,2],[3,8,4],[5,6,7]]) #position 8 is the center and not selectable\n",
    "    \n",
    "    #choose a random cell in the dividing state\n",
    "    index_pos = int(np.random.rand(1)[0]*g_index[0].shape[0])\n",
    "    m = g_index[0][index_pos]\n",
    "    n = g_index[1][index_pos]\n",
    "    #define the neighborhood\n",
    "    nb = grid[m-1:m+2,n-1:n+2]  #nb = neighborhood\n",
    "    \n",
    "    \n",
    "    #define the free spaces in nb\n",
    "    fs = np.where(nb == 0)\n",
    "    \n",
    "    if fs[0].shape[0] > 0: # go if there are a place in the neighbour\n",
    "\n",
    "        if len(fs[0]) == 1:\n",
    "            grid[m,n] = 1     # then, that position will not be evaluated again\n",
    "\n",
    "         #grown over an empty position\n",
    "        new_pos = int(np.random.rand(1)[0]*fs[0].shape[0]) #new pos in the neighbour matrix\n",
    "        m_new = m + fs[0][new_pos] - 1  #-1 to convert [0 1 2] to [-1 0 1]\n",
    "        n_new = n + fs[1][new_pos] - 1\n",
    "        grid[m_new, n_new] = 2\n",
    "    \n",
    "        #save the cell grid index positions\n",
    "        cell_index = [m,n]\n",
    "        ncell_index = [m_new, n_new]\n",
    "    else:\n",
    "        grid[m,n] = 1\n",
    "    \n",
    "    return(grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#to show rule 3\n",
    "time1 = time.clock()\n",
    "\n",
    "# create a initial grid\n",
    "grid_size = 150    #evetually define as an input\n",
    "grid = np.zeros((grid_size,grid_size))\n",
    "\n",
    "# define the initial pattern\n",
    "#grid[int(grid_size/2), int(grid_size/2)]=2\n",
    "grid = initial_pattern(grid, 0)\n",
    "\n",
    "#perform the loop\n",
    "steps = 10000\n",
    "#sleep_time = 0.01\n",
    "\n",
    "#to save the last figure\n",
    "filename = 'Segregation\\\\image_%03d.jpg'\n",
    "#filename = os.path.join(fpath, 'image_%05d.jpg')\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.show()\n",
    "show_grid(grid)  #show the initial grid\n",
    "plt.title('step 0')\n",
    "\n",
    "every = 100\n",
    "count = 0\n",
    "\n",
    "for i  in range(steps):\n",
    "    #time.sleep(sleep_time)    \n",
    "    grid = rule3(grid)\n",
    "    #plt.imshow(grid, cmap=plt.cm.gray)\n",
    "    #plt.title('step '+ str(i+1))\n",
    "    #fig.canvas.draw()\n",
    "    \n",
    "    if i%every == 0 or i == steps-1:\n",
    "        count += 1\n",
    "        plt.title('step '+ str(i+1))\n",
    "        plt.imshow(grid, cmap=plt.cm.gray)\n",
    "        #plt.savefig(str(filename) + \".pdf\", transparent=True)\n",
    "        plt.savefig(filename%(count))#, transparent=True)\n",
    "#plt.imshow(grid, cmap=plt.cm.gray)\n",
    "\n",
    "elapsed = time.clock() - time1\n",
    "\n",
    "print(elapsed)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def select_cell(grid):\n",
    "    # 2 = growing cell\n",
    "    # 1 = stationary cell\n",
    "    # 0 = empty space\n",
    "    \n",
    "    #this rule makes each cell divide only one time per step\n",
    "    g_index =  np.nonzero(grid == 2) # growth index = where cell value == 2\n",
    "    \n",
    "    #choose a random cell in the dividing state\n",
    "    index_pos = int(np.random.rand(1)[0]*g_index[0].shape[0])\n",
    "    m = g_index[0][index_pos]\n",
    "    n = g_index[1][index_pos]\n",
    "\n",
    "    \n",
    "    #save the cell grid index positions\n",
    "    cell_index = [m,n]\n",
    "        \n",
    "    return(cell_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_nbhd(grid, cell_index):\n",
    "    #chek free spaces in the neighbourhood\n",
    "    # fs: array\n",
    "    #    index of free spaces in the neighborhood\n",
    "    m = cell_index[0]\n",
    "    n = cell_index[1]\n",
    "    \n",
    "    #define the neighborhood\n",
    "    nb = grid[m-1:m+2,n-1:n+2]  #nb = neighborhood\n",
    "    \n",
    "    #define the free spaces in nb\n",
    "    fs = np.where(nb == 0)\n",
    "    \n",
    "    return(fs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nb_prob(grid, cell_index, prob_dist = 'contact_linear'):\n",
    "    # assign division probabilities based on empty space cell contacts\n",
    "    \n",
    "    # prob_dist: uniform - contact_linear - contact_exp\n",
    "    #    contact linear is the default or if another thing is written\n",
    "    \n",
    "    # fs: array\n",
    "    #    index of free spaces in the neighborhood\n",
    "    # return\n",
    "    # prob: list\n",
    "    #    list with the [0,1] probability partition limit of each free space\n",
    "    #    e.g. prob = [0.23, 0.81, 1] --> second cell has bigger probability\n",
    "    \n",
    "    m = cell_index[0]\n",
    "    n = cell_index[1]\n",
    "    \n",
    "    #define neighborhood\n",
    "    nb = grid[m-1:m+2,n-1:n+2]  #nb = neighborhood\n",
    "    \n",
    "    #define the free spaces in nb\n",
    "    fs = np.where(nb == 0)\n",
    "    #define cell spaces in bn\n",
    "    cs = np.where(nb != 0)\n",
    "    \n",
    "    fs_num = len(fs[0])\n",
    "    prob = np.zeros(fs_num)\n",
    "    contacts = np.zeros(fs_num)\n",
    "    \n",
    "    if prob_dist != 'uniform':\n",
    "        \n",
    "        # if prob_dist is something different from the options, contact_linear is the default\n",
    "        for i in range(fs_num):\n",
    "            mg = m + fs[0][i] - 1  #-1 to convert [0 1 2] to [-1 0 1]\n",
    "            ng = n + fs[1][i] - 1\n",
    "\n",
    "            i_nb = grid[mg-1:mg+2,ng-1:ng+2] # i position neighborhood\n",
    "\n",
    "            occup = np.where(i_nb != 0)\n",
    "            contacts[i] = len(occup[0])  #save the number of contacts of this position\n",
    "\n",
    "        if prob_dist == 'contact_exp':\n",
    "            contacts = np.exp(contacts)\n",
    "            \n",
    "    else:\n",
    "        contacts = np.ones(fs_num)  #assign uniform values\n",
    "        \n",
    "    total = sum(contacts)\n",
    "\n",
    "    prob[0] = (contacts[0]/total)\n",
    "\n",
    "    for i in range(1,fs_num):\n",
    "        prob[i] = prob[i-1]+contacts[i]/total\n",
    "        \n",
    "    return(prob)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cell_divition_uniform(grid, cell_index, fs):\n",
    "    # uniform neighborhood divition probability\n",
    "    #fs: free neighborhood spaces\n",
    "    \n",
    "    m = cell_index[0]\n",
    "    n = cell_index[1]\n",
    "    \n",
    "    if len(fs[0]) == 1:\n",
    "        grid[m,n] = 1     # then, that position will not divide again\n",
    "\n",
    "    #grown over an empty position\n",
    "    #new_pos = int(np.random.rand(1)[0]*fs[0].shape[0])\n",
    "    new_pos = int(np.random.rand(1)[0]*fs[0].shape[0]) #new pos in the neighbour matrix\n",
    "    m_new = m + fs[0][new_pos] - 1  #-1 to convert [0 1 2] to [-1 0 1]\n",
    "    n_new = n + fs[1][new_pos] - 1\n",
    "    grid[m_new, n_new] = 2     # crates the new cell\n",
    "        \n",
    "    ncell_index = [m_new ,n_new]\n",
    "    \n",
    "    return(grid, ncell_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cell_divition(grid, cell_index, fs, fs_proba):\n",
    "    #fs: free neighborhood spaces\n",
    "    #fs_proba: free spaces growth probabilities\n",
    "    \n",
    "    m = cell_index[0]\n",
    "    n = cell_index[1]\n",
    "    \n",
    "    if len(fs[0]) == 1:\n",
    "        grid[m,n] = 1     # then, that position will not divide again\n",
    "\n",
    "    #grown over an empty position\n",
    "\n",
    "    rand_val = np.random.rand(1)[0]\n",
    "    # find the first position which is bigger than rand_val\n",
    "    new_pos = np.where( (fs_proba > rand_val) == True )[0][0]  #new pos in the neighbour matrix\n",
    "   \n",
    "    m_new = m + fs[0][new_pos] - 1  #-1 to convert [0 1 2] to [-1 0 1]\n",
    "    n_new = n + fs[1][new_pos] - 1\n",
    "    grid[m_new, n_new] = 2     # crates the new cell\n",
    "        \n",
    "    ncell_index = [m_new ,n_new]\n",
    "    \n",
    "    return(grid, ncell_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "im_grid = np.zeros((100,100,4))\n",
    "im_grid[:,:,0] = np.ones((100,100))*0\n",
    "im_grid[:,:,1] = np.ones((100,100))*80\n",
    "im_grid[:,:,2] = np.ones((100,100))*0\n",
    "im_grid[:,:,3] = np.ones((100,100))*1\n",
    "plt.imshow(im_grid)\n",
    "plt.figure()\n",
    "im_grid[:,:,1] = np.ones((100,100))*1\n",
    "plt.imshow(im_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initial_plasmids(grid, pattern_num = 0, num_plas = 2, max_copy = 4):\n",
    "    # grid: initial grid\n",
    "    c_index =  np.nonzero(grid) # c_index, \n",
    "    #cell_number = c_index[0].shape[0]\n",
    "    \n",
    "    gs = grid.shape\n",
    "    pattern = np.zeros((gs[0],gs[1],max_copy)) #initialize the pattern array\n",
    "    \n",
    "    # add different patterns\n",
    "    if pattern_num == 0:  #random plasmid pattern\n",
    "        for i in range(c_index[0].shape[0]): #assign a random plasmid pattern to each cell position\n",
    "            pattern[c_index[0][i],c_index[1][i],:] = ((num_plas +1 )*np.random.rand(max_copy)).astype(int) \n",
    "                                                            #num_plas +1 to add \"no-plasmid\" state\n",
    "            \n",
    "    elif pattern_num == 1:\n",
    "        pattern = np.ones((grid.shape))\n",
    "        \n",
    "    return(pattern)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def role_divideFlag(plasmids):\n",
    "    #plasmids: cell plasmids vector\n",
    "    \n",
    "    max_plasmids = plasmids.shape[0]\n",
    "    num_plasmids = np.nonzero(plasmids)[0].shape[0]\n",
    "    divisor = max_plasmids*1.1  #arbitrary defined to make division(max_plasmids number) < 1\n",
    "       \n",
    "    # make a cuadratic function of probabilities\n",
    "    probability = (num_plasmids/divisor)**2 \n",
    "    #if a cell has no plasmids --> will not divide\n",
    "    \n",
    "    if np.random.rand(1) < probability:\n",
    "        return(1) # divide\n",
    "    else:\n",
    "        return(0) # not divide\n",
    "    \n",
    "    #Probability tables\n",
    "    #plasmid_nums = np.arange(max_plasmids +1)\n",
    "    #probability = (plasmid_nums/divisor)**2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_image(grid, plasgrid):\n",
    "    im_s = plasgrid.shape\n",
    "    aux_imR = np.zeros((im_s[0],im_s[1],im_s[2]))\n",
    "    aux_imG = np.zeros((im_s[0],im_s[1],im_s[2]))\n",
    "    \n",
    "    for i in range(im_s[2]):\n",
    "        aux_imR[:,:,i] = 1*(plasgrid[:,:,i]==1)\n",
    "        aux_imG[:,:,i] = 1*(plasgrid[:,:,i]==2) \n",
    "    \n",
    "    aux_imR = np.sum(aux_imR,axis=2)\n",
    "    aux_imG = np.sum(aux_imG,axis=2)\n",
    "    aux_transparency = 0.5*(grid[:,:]==1) + 1*(grid[:,:]==2)\n",
    "    \n",
    "    # create the image\n",
    "    im_grid = np.zeros((im_s[0],im_s[1],im_s[2]))\n",
    "    im_grid[:,:,0] = np.multiply(np.ones((im_s[0],im_s[1])),aux_imR)\n",
    "    im_grid[:,:,1] = np.multiply(np.ones((im_s[0],im_s[1])),aux_imG)\n",
    "    #im_grid[:,:,2] = np.ones((100,100))*250\n",
    "    im_grid[:,:,3] = np.multiply(np.ones((im_s[0],im_s[1])),aux_transparency)\n",
    "     # stationary cell -> transparency = 0.5)\n",
    "    \n",
    "    return(im_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_image2(grid, plasgrid):\n",
    "    im_s = plasgrid.shape\n",
    "    aux_imR = np.zeros((im_s[0],im_s[1],im_s[2]))\n",
    "    aux_imG = np.zeros((im_s[0],im_s[1],im_s[2]))\n",
    "    \n",
    "    for i in range(im_s[2]):\n",
    "        aux_imR[:,:,i] = 1*(plasgrid[:,:,i]==1)\n",
    "        aux_imG[:,:,i] = 1*(plasgrid[:,:,i]==2) \n",
    "    \n",
    "    aux_imR = np.multiply(1*(np.sum(aux_imR,axis=2)>0),50*(grid[:,:]==1)) + 1*(np.sum(aux_imR,axis=2)>0)\n",
    "    aux_imG = np.multiply(1*(np.sum(aux_imG,axis=2)>0),50*(grid[:,:]==1)) + 1*(np.sum(aux_imG,axis=2)>0)\n",
    "    \n",
    "    # create the image\n",
    "    im_grid = np.zeros((im_s[0],im_s[1],3))\n",
    "    im_grid[:,:,0] = np.multiply(np.ones((im_s[0],im_s[1])),aux_imR)\n",
    "    im_grid[:,:,1] = np.multiply(np.ones((im_s[0],im_s[1])),aux_imG)\n",
    "    #im_grid[:,:,2] = np.ones((100,100))*250\n",
    "    \n",
    "    return(im_grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plasmid_gProb(g_ratio=[1,1], p_types = [1,2]):\n",
    "    #define a growtn probability (=velocity) based on the plasmids\n",
    "    #g_ratio: ratio of growth rate between genotypes (i.e. plasmids)\n",
    "    #p_types: plasmids types or labels\n",
    "    \n",
    "    #built the probability class vector\n",
    "    cat_len = len(g_ratio)\n",
    "    probs = np.zeros(cat_len)\n",
    "    \n",
    "    denominator = sum(g_ratio)\n",
    "\n",
    "    probs[0] = g_ratio[0]/denominator\n",
    "\n",
    "    for i in range(1,cat_len):\n",
    "        probs[i] = probs[i-1]+g_ratio[i]/denominator\n",
    "    \n",
    "    return(probs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plasm_g_test(plasmids,probs):\n",
    "        #perform the probability test\n",
    "    rand_val = np.random.rand(1)\n",
    "    pos = np.where( (probs > rand_val) == True )[0][0]\n",
    "    \n",
    "    ptype = pos + 1\n",
    "    found = np.where(plasmids == ptype)[0]\n",
    "    growth = False\n",
    "    \n",
    "    if found.size>0:\n",
    "        growth = True\n",
    "    return(growth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cell_ratio(plasmgrid, ptype = [1,2]):\n",
    "    \n",
    "    c_num_plasm = np.sum(plasmgrid>0, axis=2) #number of plasmids in each grid\n",
    "\n",
    "    plasm_sum = np.sum(plasmgrid, axis = 2)\n",
    "    divition = np.divide(plasm_sum,c_num_plasm)\n",
    "    \n",
    "    #total = np.sum(np.isnan(divition) == False, axis = (0,1)) #it include cells with mix plasmids\n",
    "    found = np.zeros(len(ptype))\n",
    "    total = 0\n",
    "    for i in range(len(ptype)):\n",
    "        found[i] = len(np.where(divition == ptype[i])[0])\n",
    "        total += found[i]\n",
    "    \n",
    "    ratio = found[0]/total\n",
    "        \n",
    "    return(ratio)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "count=0\n",
    "plasmids = np.ones(4)*2\n",
    "for i in range(1000):\n",
    "    plas_probs = plasmid_gProb(g_ratio= [1,2])\n",
    "    ifG = plasm_g_test(plasmids, plas_probs)\n",
    "    if ifG == True:\n",
    "        count+=1\n",
    "print(count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#main\n",
    "sim_num = 1\n",
    "all_ratios = []\n",
    "for j in range(sim_num):\n",
    "    time1 = time.clock()\n",
    "    # create a initial empty grid\n",
    "    grid_size = 1000    #evetually define as an input\n",
    "    grid = np.zeros((grid_size,grid_size))\n",
    "\n",
    "    # define the initial grid and plasmid pattern\n",
    "    grid = initial_pattern(grid, 1)\n",
    "#    show_grid(grid)\n",
    "    plasm_grid = initial_plasmids(grid)\n",
    "\n",
    "    # Show the initial state\n",
    "#    im_grid = create_image(grid, plasm_grid)\n",
    "#    plt.imshow(im_grid)\n",
    "\n",
    "    #perform the loop\n",
    "    steps = 100000\n",
    "    #sleep_time = 0.01\n",
    "\n",
    "    #to save the last figure\n",
    "    #filename = 'null'\n",
    "#    filename = 'Segregation\\\\ratios\\\\image_%03d.jpg'\n",
    "    filename = 'Seg_1ratio.jpg'\n",
    "#    fig = plt.figure()\n",
    "#    plt.title('step 0')\n",
    "\n",
    "    #this two lines to save secuential images\n",
    "    every = 100   #every how many save images\n",
    "    count = 0\n",
    "\n",
    "    #define plasmid grwoth ratio\n",
    "    plas_probs = plasmid_gProb(g_ratio= [1,1])\n",
    "    # g_ratio = [2,1] --> plasmid 1 divide twice fast than plasmid 2\n",
    "    ratios = []  #to store the cell type ratios\n",
    "    for i  in range(steps):\n",
    "\n",
    "        #select a random growing cell\n",
    "        cell_pos = select_cell(grid)\n",
    "        free_nb = check_nbhd(grid, cell_pos)\n",
    "\n",
    "        if free_nb[0].shape[0] > 0: # go if there is a place in the neighborhood\n",
    "\n",
    "            plasmids = plasm_grid[cell_pos[0], cell_pos[1],:] #get its plasmids\n",
    "            c_growth = plasm_g_test(plasmids, plas_probs)\n",
    "\n",
    "            if c_growth == True: \n",
    "\n",
    "                # tal vez poner esto arriba y chequear que los plasmidos no estene en cero\n",
    "\n",
    "                #update its plasmids and cell state, n:new\n",
    "                n_plasmids, n_state = plasmid_update(plasmids, cell_pos)\n",
    "\n",
    "                plasm_grid[cell_pos[0], cell_pos[1],:] = n_plasmids\n",
    "                grid[cell_pos[0], cell_pos[1]] = n_state\n",
    "                 #state will not be evaluated before role_divide\n",
    "                 #role_divide function shouldn´t allow divition of that cell\n",
    "\n",
    "                divide_flag = role_divideFlag(n_plasmids)\n",
    "\n",
    "                #perform the divition if flag changed\n",
    "                if divide_flag != 0:\n",
    "                    #assign a cell to a new position\n",
    "                    free_proba = nb_prob(grid, cell_pos, prob_dist = 'contact_exp')\n",
    "                    grid, nCell_pos = cell_divition(grid, cell_pos, free_nb, free_proba) \n",
    "\n",
    "                    #split the mother plasmids\n",
    "                    m_plasmids, c_plasmids = divide_plasmids(n_plasmids)\n",
    "\n",
    "                    #assign mother and child plasmids\n",
    "                    plasm_grid[cell_pos[0], cell_pos[1],:] = m_plasmids\n",
    "                    plasm_grid[nCell_pos[0], nCell_pos[1],:] = c_plasmids\n",
    "        else:\n",
    "            grid[cell_pos[0],cell_pos[1]] = 1\n",
    "\n",
    "        #save cell type ratios\n",
    "        if i%every == 0:\n",
    "            ratios.append(cell_ratio(plasm_grid))\n",
    "        \n",
    "        #Plot the result\n",
    "        if i == steps-1:\n",
    "        #if i%every == 0 or i == steps-1:\n",
    "\n",
    "        #    count += 1\n",
    "\n",
    "            plt.title('step '+ str(i+1))\n",
    "            im_grid = create_image2(grid, plasm_grid)\n",
    "            plt.imshow(im_grid)\n",
    "        #    #fig.canvas.draw()\n",
    "        #    #plt.savefig(str(filename) + \".pdf\", transparent=True)\n",
    "            #plt.savefig(filename%(count), transparent=True)\n",
    "        #    plt.savefig(filename%(j), transparent=True)\n",
    "            plt.savefig(filename, transparent=True)\n",
    "            \n",
    "    all_ratios.append(np.asarray(ratios))\n",
    "    \n",
    "    elapsed = time.clock() - time1\n",
    "    print(elapsed)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_ratio = 0\n",
    "plt.figure()\n",
    "for i in range(len(all_ratios)):\n",
    "    plt.plot(all_ratios[i])\n",
    "    mean_ratio += all_ratios[i][-1]\n",
    "plt.show()\n",
    "mean_ratio= mean_ratio/len(all_ratios)\n",
    "print(mean_ratio)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio11=all_ratios\n",
    "ratio11_mean = mean_ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "for i in range(len(ratio11)):\n",
    "    plt.plot(ratio11[i])\n",
    "\n",
    "plt.title('growth ratio 1:1')\n",
    "plt.ylabel('cell type ratio')\n",
    "plt.xlabel('check point step number')\n",
    "plt.show()\n",
    "plt.savefig('ratio11', transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio32=all_ratios\n",
    "ratio32_mean = mean_ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio43=all_ratios\n",
    "ratio43_mean = mean_ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "for i in range(len(ratio43)):\n",
    "    plt.plot(ratio43[i])\n",
    "\n",
    "plt.title('growth ratio 4:3')\n",
    "plt.ylabel('cell type ratio')\n",
    "plt.xlabel('check point step number')\n",
    "plt.show()\n",
    "plt.savefig('ratio43', transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio109=all_ratios\n",
    "ratio109_mean = mean_ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_ratios = [ratio11_mean,ratio109_mean,ratio43_mean,ratio32_mean]\n",
    "expected = [1/2, 10/19, 4/7, 3/5]\n",
    "plt.figure()\n",
    "\n",
    "plt.plot([1, 2, 3, 4], mean_ratios, 'bo', label= 'observed ratio')\n",
    "plt.plot([1, 2, 3, 4], expected, 'ro', label = 'growth ratio')\n",
    "plt.ylabel('cell type ratio')\n",
    "plt.legend()\n",
    "plt.xticks([])\n",
    "#plt.xticks(['1:1','10:9','4:3','3:2'])\n",
    "plt.xticks([1,2,3,4], ['1:1','10:9','4:3','3:2'])\n",
    "plt.show()\n",
    "\n",
    "plt.savefig('ratios_obs_exp', transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_ratios = [ratio11_mean,ratio109_mean,ratio43_mean,ratio32_mean]\n",
    "expected = [1/2, 10/19, 4/7, 3/5]\n",
    "plt.figure()\n",
    "\n",
    "plt.plot(expected, mean_ratios, 'bo')\n",
    "plt.xlabel\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plasmid_update(plasmids, pos_index):\n",
    "    #plasmids: vector with plasmids. e.g [0,1,1,0,2]\n",
    "    state = 2 # cell state = growing state\n",
    "    \n",
    "    plasmids_pos = np.nonzero(plasmids)\n",
    "    empty_pos = np.where(plasmids == 0)\n",
    "    \n",
    "    num_plas = plasmids_pos[0].shape[0]\n",
    "    \n",
    "    if num_plas == 0:\n",
    "        #it means no plasmid in the cell\n",
    "        state = 1   #to not evaluate this cell in the loop again\n",
    "        \n",
    "    elif num_plas == plasmids.shape[0]:\n",
    "        #it means all plasmids positions are full\n",
    "        return(plasmids, state)\n",
    "    \n",
    "    else:\n",
    "        copied_pos = np.random.randint(num_plas)\n",
    "        plasmids[empty_pos[0][0]] = plasmids[plasmids_pos[0][copied_pos]]\n",
    "        #copy the plasmid in the first free space\n",
    "    \n",
    "    return(plasmids, state)   \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def divide_plasmids(plasmids):\n",
    "    #plasmids: cell plasmids\n",
    "    p_size = plasmids.size\n",
    "    \n",
    "    mother_p = np.zeros(p_size)\n",
    "    child_p = np.zeros(p_size)\n",
    "    \n",
    "    np.random.shuffle(plasmids) #shuffle the plasmids\n",
    "    \n",
    "    if (p_size & 1) == 1: #odd case\n",
    "        \n",
    "        #sum a random value to choose which cell keep more plasmids\n",
    "        rand_val = np.random.rand(1)\n",
    "        half_p = int(p_size/2 + rand_val) \n",
    "        \n",
    "    else:  #even case\n",
    "        half_p = int(p_size/2)\n",
    "    \n",
    "    mother_p[:half_p] = plasmids[:half_p]\n",
    "    child_p[half_p:]= plasmids[half_p:] \n",
    "    \n",
    "    return(mother_p, child_p)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initial_pattern(grid, pattern_num):\n",
    "    \n",
    "    pattern = {}  #initiate initial pattern dictionary\n",
    "    \n",
    "    # add different patterns\n",
    "    pattern[0] = np.array([[2]])\n",
    "    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]])\n",
    "    pattern[2] = np.ones((2,35))*2\n",
    "    \n",
    "    #make elements which are not in the border to be = 1\n",
    "    fixed_pat = pattern[pattern_num]\n",
    "    \n",
    "    \n",
    "    #put the pattern in the grid\n",
    "    gs = grid.shape\n",
    "    \n",
    "    m0 = int(gs[0]/2)\n",
    "    n0 = int(gs[1]/2)\n",
    "    \n",
    "    ps = fixed_pat.shape\n",
    "    mpm = int(ps[0]/2)\n",
    "    npm = int(ps[1]/2)\n",
    "            \n",
    "    for i in range(ps[0]):\n",
    "        for j in range(ps[1]):\n",
    "            m = m0 + (i - mpm)\n",
    "            n = n0 + (j - npm)\n",
    "            \n",
    "            grid[m,n] = fixed_pat[i,j]\n",
    "\n",
    "\n",
    "    return(grid)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#perform the loop\n",
    "steps = 50\n",
    "sleep_time = 0.1\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.show()\n",
    "show_grid(grid)  #show the initial grid\n",
    "plt.title('step 0')\n",
    "\n",
    "for i  in range(steps):\n",
    "    time.sleep(sleep_time)    \n",
    "    grid = rule1(grid)\n",
    "    plt.imshow(grid, cmap=plt.cm.gray)\n",
    "    plt.title('step '+ str(i+1))\n",
    "    fig.canvas.draw()\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#perform the loop\n",
    "steps = 10\n",
    "sleep_time = 0.5\n",
    "\n",
    "show_grid(grid)  #show the initial grid\n",
    "plt.figure(0)\n",
    "for i  in range(steps):\n",
    "    time.sleep(sleep_time)    \n",
    "    grid = rule1(grid)\n",
    "    show_grid(grid)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#to show\n",
    "time1 = time.clock()\n",
    "\n",
    "# create a initial grid\n",
    "grid_size = 150    #evetually define as an input\n",
    "grid = np.zeros((grid_size,grid_size))\n",
    "\n",
    "# define the initial pattern\n",
    "#grid[int(grid_size/2), int(grid_size/2)]=2\n",
    "grid = initial_pattern(grid, 2)\n",
    "\n",
    "#perform the loop\n",
    "steps = 10000\n",
    "#sleep_time = 0.01\n",
    "\n",
    "#to save the last figure\n",
    "filename = 'null'\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.show()\n",
    "show_grid(grid)  #show the initial grid\n",
    "plt.title('step 0')\n",
    "\n",
    "for i  in range(steps):\n",
    "    #time.sleep(sleep_time)    \n",
    "    #grid = rule1(grid)\n",
    "    grid = rule3(grid)\n",
    "    #plt.imshow(grid, cmap=plt.cm.gray)\n",
    "    plt.title('step '+ str(i+1))\n",
    "    #fig.canvas.draw()\n",
    "    \n",
    "    if i == steps-1:\n",
    "        plt.savefig(str(filename) + \".pdf\", transparent=True)\n",
    "plt.imshow(grid, cmap=plt.cm.gray)\n",
    "\n",
    "elapsed = time.clock() - time1\n",
    "\n",
    "print(elapsed)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#to show\n",
    "time1 = time.clock()\n",
    "# create a initial grid\n",
    "grid_size = 150    #evetually define as an input\n",
    "grid = np.zeros((grid_size,grid_size))\n",
    "\n",
    "# define the initial pattern\n",
    "grid = initial_pattern(grid, 0)\n",
    "\n",
    "#perform the loop\n",
    "steps = 15000\n",
    "#sleep_time = 0.01\n",
    "\n",
    "#to save the last figure\n",
    "filename = 'null'\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.title('step 0')\n",
    "\n",
    "for i  in range(steps):\n",
    "    #time.sleep(sleep_time)    \n",
    "    #grid = rule1(grid)\n",
    "    grid = rule2(grid)\n",
    "    #plt.imshow(grid, cmap=plt.cm.gray)\n",
    "    plt.title('step '+ str(i+1))\n",
    "    \n",
    "    \n",
    "    if i == steps-1:\n",
    "        \n",
    "        plt.imshow(grid, cmap=plt.cm.gray)\n",
    "        #fig.canvas.draw()\n",
    "        plt.savefig(str(filename) + \".pdf\", transparent=True)\n",
    "        \n",
    "\n",
    "elapsed = time.clock() - time1\n",
    "print(elapsed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask[1,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vals = np.array([[-1,-1],[0,-1],[1,-1],[-1,0],[1,0],[-1,1],[0,1],[1,1],[0,0]])\n",
    "vals[4][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos = np.nonzero(grid == 0)\n",
    "print(pos[0][10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(pos[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid[grid == 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make automata loop\n",
    "\n",
    "time_steps = 100\n",
    "\n",
    "for t in time_steps:\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.randint(8) # a random number between [0,8[ --> eigh posibilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid[0,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time1 = time.clock()\n",
    "arr = np.random.rand(5000,5000)\n",
    "elapsed = time.clock() - time1\n",
    "\n",
    "print(elapsed)  # it takes 0.21186613101770035 seg for me\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.amin(np.nonzero(grid == 2)[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#hacer clases\n",
    "# to not make:\n",
    "# cell = [px,py,vx,vy]\n",
    "# cell[0]+= cell[2]*deltaT\n",
    "\n",
    "\n",
    "# make classes is the same in general terms (in computation time, etc), but is a more recomended paradigma\n",
    "# because is more clear and organized.\n",
    "\n",
    "class cell:\n",
    "    def __init__(self,px,py): #primer parametro siempre es self -->que hace referencia al objeto cell\n",
    "        self.px = px\n",
    "        self.py = py\n",
    "        self.vx = 0\n",
    "        self.vy = 0\n",
    "        \n",
    "    # la gracia es que puedo definir funciones dentro de la clase    \n",
    "    def mover(self,vx,vy,t):  #siempre el primero es self!\n",
    "        self.px += vx*t\n",
    "        self.py += vy*t\n",
    "        \n",
    "celula = cell(-1,1)\n",
    "\n",
    "print(celula.px) #--> arroja -1\n",
    "ceula.px +=1\n",
    "print(celula.px) #--> arroja 0\n",
    "\n",
    "celula.mover(-1,1,1)\n",
    "\n",
    "lista = []\n",
    "for n in range(10):\n",
    "    lista.append(cell(n/10,n/10))\n",
    "    \n",
    "# or to make some computations faster when it is a big number of elements (>1000)\n",
    "\n",
    "class sim:\n",
    "    def __init__(self, ncells):\n",
    "        positions = np.array((ncells,2))\n",
    "\n",
    "#Also you can create an inheritance class\n",
    "class bacterium(cell):   #that means bacterium \"is a\" cell\n",
    "    def infect():   #then you can add new functions to the class bacterium\n",
    "        \n",
    "# todos los objetos en python pueden guardarse con el packete pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}