CALCULATIONS_FOR_SIR_THRESHOLDS_s2p5_GAMMA_2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pylab import *\n",
    "import random as rd\n",
    "import numpy as np\n",
    "import networkx as nx\n",
    "import csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def rectangular_geometric_lattice(F,C,r):\n",
    "    \n",
    "    a = ((C - 1.0) / (F - 1.0))**0.5\n",
    "    d = 1.0/(a * (F-1))\n",
    "\n",
    "    g = nx.Graph()\n",
    "    \n",
    "    pos_x = 0.0\n",
    "    pos_y = 0.0\n",
    "    \n",
    "    for i in range(0,int(F*C)):\n",
    "        \n",
    "        if i == 0:\n",
    "\n",
    "            g.add_node(i)\n",
    "            g.node[i]['x'] = pos_x\n",
    "            g.node[i]['y'] = pos_y\n",
    "            \n",
    "        else:\n",
    "            \n",
    "            g.add_node(i)\n",
    "            \n",
    "            if i % C == 0.0:\n",
    "                \n",
    "                pos_y += d\n",
    "                pos_x = 0.0\n",
    "                \n",
    "                g.node[i]['x'] = pos_x\n",
    "                g.node[i]['y'] = pos_y\n",
    "                \n",
    "            else:\n",
    "                \n",
    "                pos_x += d\n",
    "                \n",
    "                g.node[i]['x'] = pos_x\n",
    "                g.node[i]['y'] = pos_y\n",
    "\n",
    "\n",
    "    for i in g.nodes_iter():\n",
    "\n",
    "        for j in g.nodes_iter():\n",
    "\n",
    "            if i != j and g.has_edge(i,j)==False and\\\n",
    "            (g.node[i]['x'] - g.node[j]['x'])**2.0 + (g.node[i]['y']-g.node[j]['y'])**2.0 <= r**2.0:\n",
    "\n",
    "                g.add_edge(i,j)\n",
    "                \n",
    "    return g\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#ARRANGEMENTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def monocrop(F,C,fraction_infected_t0,inf_border):\n",
    "    \n",
    "    crop_inmune = []\n",
    "    crop_suscep = []\n",
    "    \n",
    "    for i in range(0,int(F*C)):\n",
    "            \n",
    "        crop_suscep.append(i)\n",
    "            \n",
    "    if inf_border == False:\n",
    "        \n",
    "        crop_infect = rd.sample(crop_suscep,int(fraction_infected_t0*len(crop_suscep)))\n",
    "        \n",
    "    else:\n",
    "        \n",
    "        borde_rec = []\n",
    "        \n",
    "        for i in range(0,int(C)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(int((F-1)*C),int(C*F)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(1,int(F-1)):\n",
    "            \n",
    "            borde_rec.append(int(C*i))\n",
    "            borde_rec.append(int(C*i)+int(C-1))\n",
    "            \n",
    "        borde_sus = [i for i in borde_rec if i not in crop_inmune]\n",
    "        crop_infect = rd.sample(borde_sus,int(1))\n",
    "    \n",
    "    return crop_inmune, crop_suscep, crop_infect\n",
    "\n",
    "def strips(F,C,fraction_infected_t0,inf_border): #Pensado para F=20, C= 50\n",
    "    \n",
    "    plantas = range(0,int(F*C))\n",
    "    \n",
    "    crop_inmune = []\n",
    "    \n",
    "    filas = [0,3,4,7,8,11,12,15,16,19]\n",
    "\n",
    "    for i in filas:\n",
    "            \n",
    "        for j in range(int(i*C),int(((i+1)*C))):\n",
    "\n",
    "            crop_inmune.append(j)\n",
    "        \n",
    "    crop_suscep = [x for x in plantas if x not in crop_inmune]\n",
    "    \n",
    "    if inf_border == False:\n",
    "        \n",
    "        crop_infect = rd.sample(crop_suscep,int(fraction_infected_t0*len(crop_suscep)))\n",
    "        \n",
    "    else:\n",
    "        \n",
    "        borde_rec = []\n",
    "        \n",
    "        for i in range(0,int(C)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(int((F-1)*C),int(C*F)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(1,int(F-1)):\n",
    "            \n",
    "            borde_rec.append(int(C*i))\n",
    "            borde_rec.append(int(C*i)+int(C-1))\n",
    "            \n",
    "        borde_sus = [i for i in borde_rec if i not in crop_inmune]\n",
    "        crop_infect = rd.sample(borde_sus,int(1))\n",
    "    \n",
    "    return crop_inmune, crop_suscep, crop_infect\n",
    "\n",
    "def random(F,C,fraction_infected_t0,inf_border):\n",
    "    \n",
    "    plantas = range(0,int(F*C))\n",
    "    \n",
    "    crop_inmune = rd.sample(plantas,int(0.5*F*C))\n",
    "    crop_suscep = [i for i in plantas if i not in crop_inmune]\n",
    "    \n",
    "    if inf_border == False:\n",
    "        \n",
    "        crop_infect = rd.sample(crop_suscep,int(fraction_infected_t0*len(crop_suscep)))\n",
    "        \n",
    "    else:\n",
    "        \n",
    "        borde_rec = []\n",
    "        \n",
    "        for i in range(0,int(C)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(int((F-1)*C),int(C*F)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(1,int(F-1)):\n",
    "            \n",
    "            borde_rec.append(int(C*i))\n",
    "            borde_rec.append(int(C*i)+int(C-1))\n",
    "            \n",
    "        borde_sus = [i for i in borde_rec if i not in crop_inmune]\n",
    "        crop_infect = rd.sample(borde_sus,int(1))\n",
    "    \n",
    "    return crop_inmune, crop_suscep, crop_infect\n",
    "\n",
    "def row(F,C,fraction_infected_t0,inf_border):\n",
    "    \n",
    "    a = ((C - 1.0) / (F - 1.0))**0.5\n",
    "    d = 1.0/(a * (F-1))\n",
    "    \n",
    "    crop_inmune = []\n",
    "    crop_suscep = []\n",
    "    \n",
    "    for i in range(0,int(F)):\n",
    "            \n",
    "        if 1.0*int(i) % 2 == 0:\n",
    "            \n",
    "            crop_suscep += range(i*C,(i+1)*int(C))\n",
    "            ##print i,1.0*int(C) % 2, 'sus',crop_suscep\n",
    "            \n",
    "        else:\n",
    "            \n",
    "            crop_inmune += range(i*C, (i+1)*int(C))\n",
    "            ##print i,1.0*int(C) % 2, 'inm',crop_inmune\n",
    "    \n",
    "    if inf_border == False:\n",
    "        \n",
    "        crop_infect = rd.sample(crop_suscep,int(fraction_infected_t0*len(crop_suscep)))\n",
    "        \n",
    "    else:\n",
    "        \n",
    "        borde_rec = []\n",
    "        \n",
    "        for i in range(0,int(C)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(int((F-1)*C),int(C*F)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(1,int(F-1)):\n",
    "            \n",
    "            borde_rec.append(int(C*i))\n",
    "            borde_rec.append(int(C*i)+int(C-1))\n",
    "            \n",
    "        borde_sus = [i for i in borde_rec if i not in crop_inmune]\n",
    "        crop_infect = rd.sample(borde_sus,int(1))\n",
    "    \n",
    "    return crop_inmune, crop_suscep, crop_infect\n",
    "\n",
    "def columna(F,C,fraction_infected_t0,inf_border):\n",
    "    \n",
    "   \n",
    "    crop_inmune = []\n",
    "    crop_suscep = []\n",
    "    \n",
    "    for i in range(0,int(C*F)):\n",
    "            \n",
    "        if 1.0*int(i) % 2 == 0:\n",
    "            \n",
    "            crop_suscep.append(i)\n",
    "            ##print i,1.0*int(C) % 2, 'sus',crop_suscep\n",
    "            \n",
    "        else:\n",
    "            \n",
    "            crop_inmune.append(i)\n",
    "            ##print i,1.0*int(C) % 2, 'inm',crop_inmune\n",
    "    \n",
    "    if inf_border == False:\n",
    "        \n",
    "        crop_infect = rd.sample(crop_suscep,int(fraction_infected_t0*len(crop_suscep)))\n",
    "        \n",
    "    else:\n",
    "        \n",
    "        borde_rec = []\n",
    "        \n",
    "        for i in range(0,int(C)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(int((F-1)*C),int(C*F)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(1,int(F-1)):\n",
    "            \n",
    "            borde_rec.append(int(C*i))\n",
    "            borde_rec.append(int(C*i)+int(C-1))\n",
    "            \n",
    "        borde_sus = [i for i in borde_rec if i not in crop_inmune]\n",
    "        crop_infect = rd.sample(borde_sus,int(1))\n",
    "    \n",
    "    return crop_inmune, crop_suscep, crop_infect\n",
    "\n",
    "def chessboard(F,C,fraction_infected_t0,inf_border):\n",
    "   \n",
    "    crop_inmune = []\n",
    "    crop_suscep = []\n",
    "    \n",
    "    for i in range(0,int(F)):\n",
    "            \n",
    "        for j in range(0,int(C)):\n",
    "            \n",
    "            if i % 2 == 0:\n",
    "                \n",
    "                if j % 2 == 0:\n",
    "            \n",
    "                    crop_suscep.append(i*C+j)\n",
    "                \n",
    "                else:\n",
    "                    \n",
    "                    crop_inmune.append(i*C+j)\n",
    "                \n",
    "            else:\n",
    "                \n",
    "                if j % 2 == 0:\n",
    "            \n",
    "                    crop_inmune.append(i*C+j)\n",
    "                \n",
    "                else:\n",
    "                    \n",
    "                    crop_suscep.append(i*C+j)\n",
    "    \n",
    "    if inf_border == False:\n",
    "        \n",
    "        crop_infect = rd.sample(crop_suscep,int(fraction_infected_t0*len(crop_suscep)))\n",
    "        \n",
    "    else:\n",
    "        \n",
    "        borde_rec = []\n",
    "        \n",
    "        for i in range(0,int(C)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(int((F-1)*C),int(C*F)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(1,int(F-1)):\n",
    "            \n",
    "            borde_rec.append(int(C*i))\n",
    "            borde_rec.append(int(C*i)+int(C-1))\n",
    "            \n",
    "        borde_sus = [i for i in borde_rec if i not in crop_inmune]\n",
    "        crop_infect = rd.sample(borde_sus,int(1))\n",
    "    \n",
    "    return crop_inmune, crop_suscep, crop_infect\n",
    "\n",
    "def patches(F,C,fraction_infected_t0,Fp,Cp,inf_border):\n",
    "   \n",
    "    crop_inmune = []\n",
    "    crop_suscep = []\n",
    "    \n",
    "    lista_Fp = [i for i in range(0,int(F)) if i % Fp == 0]\n",
    "    lista_Cp = [i for i in range(0,int(C)) if i % Cp == 0]\n",
    "    \n",
    "    lista_Fp_tipo = []\n",
    "    lista_Cp_tipo = []\n",
    "    \n",
    "    for i in range(0,len(lista_Fp)):\n",
    "        \n",
    "        if i%2==0:\n",
    "            \n",
    "            lista_Fp_tipo.append(1.0) #1 es susceptible\n",
    "            \n",
    "        else:\n",
    "            \n",
    "            lista_Fp_tipo.append(0.0) #0 es inmune\n",
    "    \n",
    "    for i in range(0,len(lista_Cp)):\n",
    "        \n",
    "        if i%2 == 0:\n",
    "            \n",
    "            lista_Cp_tipo.append(1.0) #1 es susceptible\n",
    "            \n",
    "        else:\n",
    "            \n",
    "            lista_Cp_tipo.append(0.0) #0 es inmune\n",
    "\n",
    "    for i in range(0,len(lista_Cp)):\n",
    "        \n",
    "        tipo_x = lista_Cp_tipo[i]\n",
    "        \n",
    "        for x in range(lista_Cp[i],lista_Cp[i]+Cp):\n",
    "            \n",
    "            for j in range(0,len(lista_Fp)):\n",
    "                \n",
    "                tipo_y = lista_Fp_tipo[j]\n",
    "                \n",
    "                for y in range(lista_Fp[j],lista_Fp[j]+Fp):\n",
    "                    \n",
    "                    if tipo_y==tipo_x:\n",
    "                        \n",
    "                        crop_suscep.append(y*C+x)\n",
    "                        ##print x,y,y*F+x,'sus'\n",
    "                        \n",
    "                    else:\n",
    "                        \n",
    "                        crop_inmune.append(y*C+x)\n",
    "                        ##print x,y,y*F+x,'inm'\n",
    "   \n",
    "    if inf_border == False:\n",
    "        \n",
    "        crop_infect = rd.sample(crop_suscep,int(fraction_infected_t0*len(crop_suscep)))\n",
    "        \n",
    "    else:\n",
    "        \n",
    "        borde_rec = []\n",
    "        \n",
    "        for i in range(0,int(C)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(int((F-1)*C),int(C*F)):\n",
    "            \n",
    "            borde_rec.append(i)\n",
    "        \n",
    "        for i in range(1,int(F-1)):\n",
    "            \n",
    "            borde_rec.append(int(C*i))\n",
    "            borde_rec.append(int(C*i)+int(C-1))\n",
    "            \n",
    "        borde_sus = [i for i in borde_rec if i not in crop_inmune]\n",
    "        crop_infect = rd.sample(borde_sus,int(1))\n",
    "    \n",
    "    return crop_inmune, crop_suscep, crop_infect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# SIR DYNAMIC WITHOUT LONG-RANGE INTERACTIONS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def q_i_SIN(nodo,lista_p_i):\n",
    "    \n",
    "    global g,beta_cons\n",
    "\n",
    "    aux = 1.0\n",
    "    \n",
    "    for i in g.neighbors_iter(nodo):\n",
    "        \n",
    "        aux *= (1.0 - beta_cons*lista_p_i[i])\n",
    "        \n",
    "    return 1.0 - aux\n",
    "\n",
    "def Markov_evolution_SIR_SIN_cuad(crop_inmune, crop_suscep, crop_infect):\n",
    "    \n",
    "    global g,beta_cons,mu_cons,tmax,lista_infectados_ini,\\\n",
    "    lista_recuperados_ini,lista_recuperados_final\n",
    "    \n",
    "    N0 = len(g.nodes())\n",
    "    \n",
    "    lista_s_i = [] #susceptibles\n",
    "    lista_p_i = [] #infectados\n",
    "    lista_rho_i = [] #recuperados\n",
    "    \n",
    "    for i in g.nodes_iter():\n",
    "\n",
    "        lista_s_i.append(1.0) #susceptibles\n",
    "        lista_p_i.append(0.0) #infectados\n",
    "        lista_rho_i.append(0.0) #recuperados\n",
    "        \n",
    "    N0 = len(g.nodes())\n",
    "    \n",
    "    cant_inmunes_t0 = len(crop_inmune)\n",
    "    lista_infectados_ini = crop_infect[:]\n",
    "    lista_recuperados_ini = crop_inmune[:]\n",
    "    \n",
    "    ##print 'inmunes:',1.0*len(crop_inmune)/N0,'infectados:',1.0*len(crop_infect)/len(crop_suscep)\n",
    "\n",
    "    for i in crop_inmune:\n",
    "        \n",
    "        lista_s_i[i] =  0.0\n",
    "        lista_p_i[i] = 0.0\n",
    "        lista_rho_i[i] = 1.0\n",
    "        \n",
    "    for i in crop_infect :\n",
    "        \n",
    "        lista_s_i[i] =  0.0\n",
    "        lista_p_i[i] = 1.0\n",
    "        lista_rho_i[i] = 0.0\n",
    "\n",
    "    \n",
    "    cantidad_inf.append(1.0*sum(lista_p_i)/(N0-cant_inmunes_t0))\n",
    "    cantidad_sus.append(1.0*sum(lista_s_i)/(N0-cant_inmunes_t0))\n",
    "    cantidad_rec.append(1.0*(sum(lista_rho_i)-cant_inmunes_t0)/(N0-cant_inmunes_t0))\n",
    "\n",
    "    \n",
    "    t = 0\n",
    "    terminar = False\n",
    "    \n",
    "    while t <= tmax and terminar==False:\n",
    "    \n",
    "        nuevo_s_i = []\n",
    "        nuevo_p_i = []\n",
    "        nuevo_rho_i = []\n",
    "\n",
    "        for a in g.nodes_iter():\n",
    "\n",
    "            \n",
    "            if lista_rho_i[a] == 1.0:\n",
    "                \n",
    "                nuevo_s_i.append(0.0)\n",
    "                nuevo_p_i.append(0.0)\n",
    "                nuevo_rho_i.append(1.0)\n",
    "                \n",
    "            else:  \n",
    "\n",
    "                q_i_nodo = q_i_SIN(a,lista_p_i)\n",
    "\n",
    "                aux_nuevo_p_i = lista_p_i[a]*(1.0 - mu_cons) + (1.0 - lista_p_i[a] - lista_rho_i[a])*q_i_nodo\n",
    "                aux_nuevo_rho_i = lista_rho_i[a] + mu_cons*lista_p_i[a]\n",
    "\n",
    "                nuevo_s_i.append(1.0 - aux_nuevo_p_i - aux_nuevo_rho_i)\n",
    "                nuevo_p_i.append(aux_nuevo_p_i)\n",
    "                nuevo_rho_i.append(aux_nuevo_rho_i)\n",
    "        \n",
    "        if max(nuevo_p_i) < 1e-6:\n",
    "\n",
    "            terminar = True\n",
    "\n",
    "        lista_s_i = nuevo_s_i[:]\n",
    "        lista_p_i = nuevo_p_i[:]\n",
    "        lista_rho_i = nuevo_rho_i[:]\n",
    "\n",
    "        \n",
    "        cantidad_inf.append(1.0*sum(lista_p_i)/(N0-cant_inmunes_t0))\n",
    "        cantidad_sus.append(1.0*sum(lista_s_i)/(N0-cant_inmunes_t0))\n",
    "        cantidad_rec.append(1.0*(sum(lista_rho_i)-cant_inmunes_t0)/(N0-cant_inmunes_t0))\n",
    "\n",
    "        t += 1\n",
    "        \n",
    "    \n",
    "    ##print 't_sim:',t\n",
    "    \n",
    "    lista_recuperados_final=[i for i in g.nodes_iter() if lista_rho_i[i] > lista_s_i[i]]\n",
    "    \n",
    "    return 1.0*(sum(lista_rho_i)-cant_inmunes_t0)/(N0-cant_inmunes_t0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def matriz_distancia():\n",
    "    \n",
    "    global g\n",
    "    \n",
    "    dic_shortest = nx.shortest_path(g)\n",
    "\n",
    "    M_nueva = []\n",
    "\n",
    "    for i in range(0,len(g.nodes())):\n",
    "\n",
    "        fila_aux = []\n",
    "\n",
    "        fila = dic_shortest.get(i)\n",
    "\n",
    "        for k in range(0,len(g.nodes())):\n",
    "\n",
    "            fila_aux.append(len(fila.get(k))-1.0)\n",
    "\n",
    "        M_nueva.append(fila_aux)\n",
    "    \n",
    "    return M_nueva\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def matriz_alargamiento():\n",
    "    \n",
    "    global g,crop_inmune\n",
    "    \n",
    "    dic_shortest = nx.shortest_path(g)\n",
    "\n",
    "    M_nueva = []\n",
    "\n",
    "    for i in range(0,len(g.nodes())):\n",
    "\n",
    "        fila_aux = []\n",
    "\n",
    "        fila = dic_shortest.get(i)\n",
    "\n",
    "        for k in range(0,len(g.nodes())):\n",
    "            \n",
    "            camino_ik = fila.get(k)\n",
    "            \n",
    "            aux = 0.0\n",
    "            \n",
    "            for p in camino_ik:\n",
    "                \n",
    "                if p in crop_inmune:\n",
    "                    \n",
    "                    aux += 1.0\n",
    "\n",
    "            fila_aux.append(aux)\n",
    "\n",
    "        M_nueva.append(fila_aux)\n",
    "    \n",
    "    return M_nueva"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# SIR DYNAMIC for a given value of gamma (i.e., the trap strength)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def q_i_MELLIN_SEXY3_mod(nodo,lista_p_i,s,gamma):\n",
    "    \n",
    "    global g,beta_cons,M_distancia,M_alargamiento,crop_suscep\n",
    "\n",
    "    aux = 1.0\n",
    "    \n",
    "    for i in crop_suscep:\n",
    "\n",
    "        if i != nodo:\n",
    "            \n",
    "            if s >= 100 and M_alargamiento[nodo][i] == 0.0:\n",
    "\n",
    "                aux *= (1.0 - 1.0*beta_cons*lista_p_i[i])\n",
    "                \n",
    "            elif s < 100 and M_alargamiento[nodo][i] == 0.0:\n",
    "\n",
    "                aux *= (1.0 - 1.0*beta_cons*lista_p_i[i]/((M_distancia[nodo][i])**s))\n",
    "                \n",
    "            else:\n",
    "                \n",
    "                aux *= (1.0 - 1.0*beta_cons*lista_p_i[i]/((M_distancia[nodo][i])**(gamma*s)))\n",
    " \n",
    "    return (1.0 - aux)\n",
    "    \n",
    "\n",
    "def Markov_evolution_SIR_MELLIN_SEXY_cuad3_time_mod(s,crop_inmune, crop_suscep, crop_infect):\n",
    "    \n",
    "    global g,beta_cons,mu_cons,tmax,lista_infectados_ini,\\\n",
    "    lista_recuperados_ini,lista_recuperados_final,M_distancia,M_alargamiento,rho_est,gamma\n",
    "    \n",
    "    N0 = len(g.nodes())\n",
    "    \n",
    "    lista_s_i = [] #susceptibles\n",
    "    lista_p_i = [] #infectados\n",
    "    lista_rho_i = [] #recuperados\n",
    "    \n",
    "    for i in g.nodes_iter():\n",
    "\n",
    "        lista_s_i.append(1.0) #susceptibles\n",
    "        lista_p_i.append(0.0) #infectados\n",
    "        lista_rho_i.append(0.0) #recuperados\n",
    "        \n",
    "    N0 = len(g.nodes())\n",
    "    \n",
    "    cant_inmunes_t0 = len(crop_inmune)\n",
    "    lista_infectados_ini = crop_infect[:]\n",
    "    lista_recuperados_ini = crop_inmune[:]\n",
    "    \n",
    "    ##print 'inmunes:',1.0*len(crop_inmune)/N0,'infectados:',1.0*len(crop_infect)/len(crop_suscep)\n",
    "\n",
    "\n",
    "    for i in crop_inmune:\n",
    "        \n",
    "        lista_s_i[i] =  0.0\n",
    "        lista_p_i[i] = 0.0\n",
    "        lista_rho_i[i] = 1.0\n",
    "        \n",
    "    for i in crop_infect :\n",
    "        \n",
    "        lista_s_i[i] =  0.0\n",
    "        lista_p_i[i] = 1.0\n",
    "        lista_rho_i[i] = 0.0\n",
    "\n",
    "    \n",
    "    cantidad_inf.append(1.0*sum(lista_p_i)/(N0-cant_inmunes_t0))\n",
    "    cantidad_sus.append(1.0*sum(lista_s_i)/(N0-cant_inmunes_t0))\n",
    "    cantidad_rec.append(1.0*(sum(lista_rho_i)-cant_inmunes_t0)/(N0-cant_inmunes_t0))\n",
    "    \n",
    "    \n",
    "    t = 0\n",
    "    terminar = False\n",
    "    \n",
    "    while t <= tmax and terminar==False:\n",
    "    \n",
    "        nuevo_s_i = []\n",
    "        nuevo_p_i = []\n",
    "        nuevo_rho_i = []\n",
    "\n",
    "        for a in g.nodes_iter():\n",
    "            \n",
    "            if g.node[a]['buscar'] == True: #lista_p_i[a] < lista_rho_i[a]:\n",
    "                \n",
    "                g.node[a]['time'] = t\n",
    "                \n",
    "                if 1.0 - lista_p_i[a] - lista_rho_i[a] < lista_rho_i[a]:\n",
    "                \n",
    "                    g.node[a]['buscar'] = False\n",
    "            \n",
    "            if s >= 100:\n",
    "            \n",
    "                q_i_nodo = q_i_SIN(a,lista_p_i)\n",
    "                \n",
    "            else:\n",
    "                \n",
    "                q_i_nodo = q_i_MELLIN_SEXY3_mod(a,lista_p_i,s,gamma) #q_i_MELLIN(a,lista_p_i,s)\n",
    "\n",
    "            aux_nuevo_p_i = lista_p_i[a]*(1.0 - mu_cons) + (1.0 - lista_p_i[a] - lista_rho_i[a])*q_i_nodo\n",
    "            aux_nuevo_rho_i = lista_rho_i[a] + mu_cons*lista_p_i[a]\n",
    "\n",
    "            nuevo_s_i.append(1.0 - aux_nuevo_p_i - aux_nuevo_rho_i)\n",
    "            nuevo_p_i.append(aux_nuevo_p_i)\n",
    "            nuevo_rho_i.append(aux_nuevo_rho_i)\n",
    "\n",
    "        \n",
    "        if max(nuevo_p_i) < 1e-6:\n",
    "\n",
    "            terminar = True\n",
    "\n",
    "        lista_s_i = nuevo_s_i[:]\n",
    "        lista_p_i = nuevo_p_i[:]\n",
    "        lista_rho_i = nuevo_rho_i[:]\n",
    "        \n",
    "        cantidad_inf.append(1.0*sum(lista_p_i)/(N0-cant_inmunes_t0))\n",
    "        cantidad_sus.append(1.0*sum(lista_s_i)/(N0-cant_inmunes_t0))\n",
    "        cantidad_rec.append(1.0*(sum(lista_rho_i)-cant_inmunes_t0)/(N0-cant_inmunes_t0))\n",
    "\n",
    "        t += 1\n",
    "        \n",
    "    \n",
    "    ##print 't_sim:',t\n",
    "    \n",
    "    return 1.0*(sum(lista_rho_i)-cant_inmunes_t0)/(N0-cant_inmunes_t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#DATA GENERATION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "F=20\n",
    "C=50\n",
    "\n",
    "Fp=5\n",
    "Cp=5\n",
    "\n",
    "a_2 = (C - 1.0) / (F - 1.0)\n",
    "d_2 = 1.0/(a_2 * (F-1)*(F-1))\n",
    "\n",
    "r=1.002*(d_2**0.5)\n",
    "\n",
    "g = rectangular_geometric_lattice(F,C,r)\n",
    "M_distancia = matriz_distancia()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "fraction_infected_t0 = 0.01\n",
    "mu_cons = 0.5\n",
    "tmax = 10000\n",
    "\n",
    "gamma = 2.0\n",
    "\n",
    "mell = 2.5\n",
    "\n",
    "lista_beta = np.logspace(log10(0.02), log10(1.0), num=50, endpoint=True, base=10.0)\n",
    "\n",
    "\n",
    "lista_resul_MELLIN_fila = []\n",
    "\n",
    "lista_resul_MELLIN_columna = []\n",
    "\n",
    "lista_resul_MELLIN_alternados = []\n",
    "\n",
    "lista_resul_MELLIN_aleatorio = []\n",
    "\n",
    "lista_resul_MELLIN_patches = []\n",
    "\n",
    "lista_resul_MELLIN_strips = []\n",
    "\n",
    "lista_resul_MELLIN_monocrop = []\n",
    "\n",
    "##############################################################################################\n",
    "#############################################################################################\n",
    "\n",
    "crop_inmune, crop_suscep, crop_infect = row(F,C,fraction_infected_t0,True)\n",
    "M_alargamiento = matriz_alargamiento()\n",
    "\n",
    "for beta_cons in lista_beta:\n",
    "\n",
    "    \n",
    "    cantidad_inf = []\n",
    "    cantidad_sus = []\n",
    "    cantidad_rec = []\n",
    "    lista_infectados_ini = []\n",
    "    lista_recuperados_ini = []\n",
    "    lista_recuperados_final = []\n",
    "    \n",
    "    crop_inmune, crop_suscep, crop_infect = row(F,C,fraction_infected_t0,True)\n",
    "    \n",
    "    for i in g.nodes_iter():\n",
    "        \n",
    "        g.node[i]['time'] = 0\n",
    "        g.node[i]['buscar'] = True\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    lista_resul_MELLIN_fila.append(Markov_evolution_SIR_MELLIN_SEXY_cuad3_time_mod(mell,crop_inmune, crop_suscep, crop_infect))\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    name =  'SIR_THRESHOLDS_ROW_MELLIN_SEXY_s_2p5_F20_C50_mu_0p5_phi_0p01_SEED_BORDE_beta_0p5_ATENUADO.txt'\n",
    "    res = [lista_resul_MELLIN_fila]\n",
    "    csvfile = name\n",
    "\n",
    "    with open(csvfile, \"w\") as output:\n",
    "        writer = csv.writer(output, lineterminator='\\n')\n",
    "        writer.writerows(res)\n",
    "    \n",
    "    \n",
    "##############################################################################################################\n",
    "##############################################################################################################\n",
    "\n",
    "crop_inmune, crop_suscep, crop_infect = columna(F,C,fraction_infected_t0,True)\n",
    "M_alargamiento = matriz_alargamiento()\n",
    "\n",
    "for beta_cons in lista_beta:\n",
    "\n",
    "    cantidad_inf = []\n",
    "    cantidad_sus = []\n",
    "    cantidad_rec = []\n",
    "    lista_infectados_ini = []\n",
    "    lista_recuperados_ini = []\n",
    "    lista_recuperados_final = []\n",
    "    \n",
    "    crop_inmune, crop_suscep, crop_infect = columna(F,C,fraction_infected_t0,True)\n",
    "    \n",
    "    for i in g.nodes_iter():\n",
    "        \n",
    "        g.node[i]['time'] = 0\n",
    "        g.node[i]['buscar'] = True\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    lista_resul_MELLIN_columna.append(Markov_evolution_SIR_MELLIN_SEXY_cuad3_time_mod(mell,crop_inmune, crop_suscep, crop_infect))\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    name =  'SIR_THRESHOLDS_COLUMNA_MELLIN_SEXY_s_2p5_F20_C50_mu_0p5_phi_0p01_SEED_BORDE_beta_0p5_ATENUADO.txt'\n",
    "    res = [lista_resul_MELLIN_columna]\n",
    "    csvfile = name\n",
    "\n",
    "    with open(csvfile, \"w\") as output:\n",
    "        writer = csv.writer(output, lineterminator='\\n')\n",
    "        writer.writerows(res)\n",
    "    \n",
    "    \n",
    "##############################################################################################################\n",
    "##############################################################################################################\n",
    "crop_inmune, crop_suscep, crop_infect = chessboard(F,C,fraction_infected_t0,True)\n",
    "M_alargamiento = matriz_alargamiento()\n",
    "\n",
    "for beta_cons in lista_beta:\n",
    "    \n",
    "    cantidad_inf = []\n",
    "    cantidad_sus = []\n",
    "    cantidad_rec = []\n",
    "    lista_infectados_ini = []\n",
    "    lista_recuperados_ini = []\n",
    "    lista_recuperados_final = []\n",
    "    \n",
    "    crop_inmune, crop_suscep, crop_infect = chessboard(F,C,fraction_infected_t0,True)\n",
    "    \n",
    "    for i in g.nodes_iter():\n",
    "        \n",
    "        g.node[i]['time'] = 0\n",
    "        g.node[i]['buscar'] = True\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    lista_resul_MELLIN_alternados.append(Markov_evolution_SIR_MELLIN_SEXY_cuad3_time_mod(mell,crop_inmune, crop_suscep, crop_infect))\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    name =  'SIR_THRESHOLDS_ALTERNADOS_MELLIN_SEXY_s_2p5_F20_C50_mu_0p5_phi_0p01_SEED_BORDE_beta_0p5_ATENUADO.txt'\n",
    "    res = [lista_resul_MELLIN_alternados]\n",
    "    csvfile = name\n",
    "\n",
    "    with open(csvfile, \"w\") as output:\n",
    "        writer = csv.writer(output, lineterminator='\\n')\n",
    "        writer.writerows(res)\n",
    "    \n",
    "##############################################################################################################\n",
    "##############################################################################################################\n",
    "crop_inmune, crop_suscep, crop_infect = patches(F,C,fraction_infected_t0,Fp,Cp,True)\n",
    "M_alargamiento = matriz_alargamiento()\n",
    "\n",
    "for beta_cons in lista_beta:\n",
    "\n",
    "    \n",
    "    cantidad_inf = []\n",
    "    cantidad_sus = []\n",
    "    cantidad_rec = []\n",
    "    lista_infectados_ini = []\n",
    "    lista_recuperados_ini = []\n",
    "    lista_recuperados_final = []\n",
    "    \n",
    "    crop_inmune, crop_suscep, crop_infect = patches(F,C,fraction_infected_t0,Fp,Cp,True)\n",
    "    \n",
    "    \n",
    "    for i in g.nodes_iter():\n",
    "        \n",
    "        g.node[i]['time'] = 0\n",
    "        g.node[i]['buscar'] = True\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    lista_resul_MELLIN_patches.append(Markov_evolution_SIR_MELLIN_SEXY_cuad3_time_mod(mell,crop_inmune, crop_suscep, crop_infect))\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    name =  'SIR_THRESHOLDS_PATCHES_MELLIN_SEXY_s_2p5_F20_C50_mu_0p5_phi_0p01_SEED_BORDE_beta_0p5_ATENUADO.txt'\n",
    "    res = [lista_resul_MELLIN_patches]\n",
    "    csvfile = name\n",
    "\n",
    "    with open(csvfile, \"w\") as output:\n",
    "        writer = csv.writer(output, lineterminator='\\n')\n",
    "        writer.writerows(res)\n",
    "    \n",
    "##############################################################################################################\n",
    "##############################################################################################################\n",
    "crop_inmune, crop_suscep, crop_infect = strips(F,C,fraction_infected_t0,True)\n",
    "M_alargamiento = matriz_alargamiento()\n",
    "\n",
    "for beta_cons in lista_beta:\n",
    "\n",
    "    \n",
    "    cantidad_inf = []\n",
    "    cantidad_sus = []\n",
    "    cantidad_rec = []\n",
    "    lista_infectados_ini = []\n",
    "    lista_recuperados_ini = []\n",
    "    lista_recuperados_final = []\n",
    "    \n",
    "    crop_inmune, crop_suscep, crop_infect = strips(F,C,fraction_infected_t0,True)\n",
    "    \n",
    "    for i in g.nodes_iter():\n",
    "        \n",
    "        g.node[i]['time'] = 0\n",
    "        g.node[i]['buscar'] = True\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    lista_resul_MELLIN_strips.append(Markov_evolution_SIR_MELLIN_SEXY_cuad3_time_mod(mell,crop_inmune, crop_suscep, crop_infect))\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    name =  'SIR_THRESHOLDS_STRIPS_MELLIN_SEXY_s_2p5_F20_C50_mu_0p5_phi_0p01_SEED_BORDE_beta_0p5_ATENUADO.txt'\n",
    "    res = [lista_resul_MELLIN_strips]\n",
    "    csvfile = name\n",
    "\n",
    "    with open(csvfile, \"w\") as output:\n",
    "        writer = csv.writer(output, lineterminator='\\n')\n",
    "        writer.writerows(res)\n",
    "    \n",
    "\n",
    "##############################################################################################################\n",
    "##############################################################################################################\n",
    "crop_inmune, crop_suscep, crop_infect = monocrop(F,C,fraction_infected_t0,True)\n",
    "M_alargamiento = matriz_alargamiento()\n",
    "\n",
    "for beta_cons in lista_beta:\n",
    "\n",
    "    \n",
    "    cantidad_inf = []\n",
    "    cantidad_sus = []\n",
    "    cantidad_rec = []\n",
    "    lista_infectados_ini = []\n",
    "    lista_recuperados_ini = []\n",
    "    lista_recuperados_final = []\n",
    "    \n",
    "    crop_inmune, crop_suscep, crop_infect = monocrop(F,C,fraction_infected_t0,True)\n",
    "    \n",
    "    for i in g.nodes_iter():\n",
    "        \n",
    "        g.node[i]['time'] = 0\n",
    "        g.node[i]['buscar'] = True\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    lista_resul_MELLIN_monocrop.append(Markov_evolution_SIR_MELLIN_SEXY_cuad3_time_mod(mell,crop_inmune, crop_suscep, crop_infect))\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    name =  'SIR_THRESHOLDS_MONOCROP_MELLIN_SEXY_s_2p5_F20_C50_mu_0p5_phi_0p01_SEED_BORDE_beta_0p5_ATENUADO.txt'\n",
    "    res = [lista_resul_MELLIN_monocrop]\n",
    "    csvfile = name\n",
    "\n",
    "    with open(csvfile, \"w\") as output:\n",
    "        writer = csv.writer(output, lineterminator='\\n')\n",
    "        writer.writerows(res)\n",
    "        \n",
    "##############################################################################################################\n",
    "##############################################################################################################\n",
    "\n",
    "for beta_cons in lista_beta:\n",
    "    \n",
    "    \n",
    "    cantidad_inf = []\n",
    "    cantidad_sus = []\n",
    "    cantidad_rec = []\n",
    "    lista_infectados_ini = []\n",
    "    lista_recuperados_ini = []\n",
    "    lista_recuperados_final = []\n",
    "    \n",
    "    crop_inmune, crop_suscep, crop_infect = random(F,C,fraction_infected_t0,True)\n",
    "    M_alargamiento = matriz_alargamiento()\n",
    "    \n",
    "    for i in g.nodes_iter():\n",
    "        \n",
    "        g.node[i]['time'] = 0\n",
    "        g.node[i]['buscar'] = True\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    lista_resul_MELLIN_aleatorio.append(Markov_evolution_SIR_MELLIN_SEXY_cuad3_time_mod(mell,crop_inmune, crop_suscep, crop_infect))\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    name =  'SIR_THRESHOLDS_ALEATORIO_MELLIN_SEXY_s_2p5_F20_C50_mu_0p5_phi_0p01_SEED_BORDE_beta_0p5_ATENUADO.txt'\n",
    "    res = [lista_resul_MELLIN_aleatorio]\n",
    "    csvfile = name\n",
    "\n",
    "    with open(csvfile, \"w\") as output:\n",
    "        writer = csv.writer(output, lineterminator='\\n')\n",
    "        writer.writerows(res) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}