EcoEvo_agriculture_MULTISEQ.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Eco-evolutionary dynamics in agriculture: a study in crop rotations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#!/usr/bin/env python\n",
    "\n",
    "from __future__ import division\n",
    "from subprocess import call\n",
    "from scipy import *\n",
    "from scipy.linalg import eigh as largest_eigh\n",
    "from scipy.integrate import odeint\n",
    "import scipy.integrate as spi\n",
    "import numpy as np\n",
    "from itertools import combinations\n",
    "from itertools import repeat\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import io\n",
    "import csv\n",
    "import sys, os\n",
    "import pandas as pd\n",
    "from pylab import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## For crop rotation cash yield and soil quality\n",
    "\n",
    "q0 = 1 # init soil quality\n",
    "y0 = 0 # init cash yield\n",
    "K = 2 #carrying capacity\n",
    "beta_h1 = -1.5  #soil contribution h1\n",
    "beta_h1_r = -1.5 #soil contribution h1_r\n",
    "beta_h2 = 1 #soil contribution h2\n",
    "y_h1 = 1 #cash contribution h1\n",
    "y_h1_r = 0.9 #cash contribution h1_r\n",
    "y_h2 = 0 #cash contribution h2\n",
    "\n",
    "\n",
    "\n",
    "## For eco-evo dynamics\n",
    "epsilon = 0.5 # pathogen clearance\n",
    "k = 0.1 # cost of resistance\n",
    "mu = 0.1 # transition probability\n",
    "a = 0.5 # GfG proportion\n",
    "sigma = 0.04 # pathogen virulence of p1\n",
    "H_dens = 50 # host density\n",
    "dp = 0.5 # pathogen deathrate\n",
    "\n",
    "\n",
    "\n",
    "## Time (needed for simulation)\n",
    "t_step = 0.005\n",
    "t_last = int(1/t_step)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Infection environment (pathogen yes/no, time of infection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pathogen true? 1 = true, 0 = false1\n",
      "infection time?0\n"
     ]
    }
   ],
   "source": [
    "pathogen_true = 1 #int(raw_input(\"pathogen true? 1 = true, 0 = false\")) #presence or absence of pathogen\n",
    "infection_time = 0 #int(raw_input(\"infection time?\")) # time the pathogen infects\n",
    "\n",
    "df = pd.DataFrame() #data frame for output results\n",
    "df2 = pd.DataFrame() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sequence population (create all possible patterns of a certain length)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "l = 10  #length\n",
    "num_seq = 2**l   #population size\n",
    "\n",
    "cash = ['h1'] * l\n",
    "cover = ['h2'] * l\n",
    "\n",
    "list_seq = [] # population\n",
    "\n",
    "def replacements(cash, cover):\n",
    "    for count in range(len(cash)):\n",
    "        for positions in combinations(range(len(cash)), count):\n",
    "            yield [cover[i] if i in positions else name for i, name in enumerate(cash)]\n",
    "    yield cover\n",
    "\n",
    "for result in replacements(cash,cover):\n",
    "    list_seq.append(result)\n",
    "    \n",
    "    \n",
    "# t_step = 0.005\n",
    "# t_last = int(1/t_step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## define a binary sequence (1,0) for each crop type\n",
    "\n",
    "def get_h_seq(seq):\n",
    "    crop_types = ['h1','h2','h1_r']\n",
    "    seq_elements = seq\n",
    "    binaries = {k:[int(k in v) for v in seq_elements]for k in crop_types}\n",
    "    h1_seq,h2_seq,h1_r_seq = binaries['h1'],binaries['h2'],binaries['h1_r']\n",
    "    \n",
    "    return h1_seq,h2_seq,h1_r_seq\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transition and fitness matrices (section 2.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def transition_matrix():\n",
    "    d = np.array([\n",
    "    [1 - 2*mu ,mu,0,0,0,mu,0,0,0,0],\n",
    "    [mu,1 - 3*mu,mu,0,0,0,mu,0,0,0],\n",
    "    [0,mu,1 - 3*mu,mu,0,0,0,mu,0,0],\n",
    "    [0,0,mu,1 - 3*mu,mu,0,0,0,mu,0],\n",
    "    [0,0,0,mu,1 - 2*mu,0,0,0,0,mu],\n",
    "    [mu,0,0,0,0,1 - 2*mu,mu,0,0,0],\n",
    "    [0,mu,0,0,0,mu,1 - 3*mu,mu,0,0],\n",
    "    [0,0,mu,0,0,0,mu,1 - 3*mu,mu,0],\n",
    "    [0,0,0,mu,0,0,0,mu,1 - 3*mu,mu],\n",
    "    [0,0,0,0,mu,0,0,0,mu,1 - 2*mu]\n",
    "    ])\n",
    "    # d = np.identity(10) ### in case we don't want any transition between types\n",
    "    return(d)\n",
    "\n",
    "q_matrix = transition_matrix()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fitness_matrix():\n",
    "    w_matrix = np.array([[1,0,0],[1.25,0,0],[1.5,0,0],[1.75,0,0],[2,0,0],\n",
    "    [a*(1 - a*k),1 - a*k,0],[a*(1.25 - a*k),1.25 - a*k,0],[a*(1.5 - a*k),1.5 - a*k,0],[a*(1.75 - a*k),1.75 - a*k,0],\n",
    "    [a*(2 - a*k),2 - a*k,0]])\n",
    "\n",
    "    return(w_matrix)\n",
    "\n",
    "w_matrix = fitness_matrix()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eco-Evo dynamics (section 2.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def infection(h1_seq,h1_r_seq,h2_seq):\n",
    "    def HP_evol(t,init): ## equations for host-pathogen coevolution\n",
    "        h = np.array(init[0:3],dtype='float16')\n",
    "        p = np.array(init[3:13],dtype='float16')\n",
    "        hdot = np.zeros(3)\n",
    "        pdot = np.zeros(10)\n",
    "        \n",
    "        for h_type in range(0,3):    \n",
    "            hdot[h_type] = h[h_type] * np.sum((-sigma * np.multiply(w_matrix[:,h_type], p,dtype='float16')))\n",
    "\n",
    "        for p_type in range(0,10):\n",
    "            pdot[p_type] = np.sum((np.multiply(q_matrix[:,p_type], p,dtype='float16') * np.sum((sigma * np.multiply(w_matrix[p_type,:], h,dtype='float16'))))) - p[p_type]*dp\n",
    "        \n",
    "        ydot = np.append(hdot,pdot)\n",
    "        return(ydot)\n",
    "    \n",
    "    \n",
    "    t_i = 0 #init time (start of season)\n",
    "    t_f = 1 #final time (end of season)\n",
    "    general_ts = []\n",
    "    general_zs = []\n",
    "    z0 = np.zeros(13)\n",
    "    for i in range(0,l):\n",
    "\n",
    "      ## intiation of host density in each season\n",
    "    \n",
    "      h1 = H_dens * (h1_seq[i]) \n",
    "      h1_r = H_dens * (h1_r_seq[i])\n",
    "      h2 = H_dens * (h2_seq[i])\n",
    "    \n",
    "      \n",
    "      z0[0:3] = [h1,h1_r,h2]\n",
    "\n",
    "      if t_i == infection_time and pathogen_true == 1: ## introducing the pathogen at infection time with a density of 1\n",
    "        p1 = 1\n",
    "        z0[3] = p1\n",
    "\n",
    "      ## ODE solver\n",
    "    \n",
    "      t_end = 1\n",
    "      t_start = 0.\n",
    "      t_step = 0.005\n",
    "      t_interval = np.arange(t_start, t_end, t_step)\n",
    "\n",
    "      myode =  spi.ode(HP_evol)\n",
    "      myode.set_integrator('lsoda',nsteps = 500000000)\n",
    "      myode.set_initial_value(z0,t_start)\n",
    "\n",
    "      ts = [0]\n",
    "      zs = [z0]  \n",
    "\n",
    "\n",
    "      while myode.successful() and myode.t < t_end:\n",
    "         myode.integrate(myode.t + t_step)\n",
    "         ts.append(myode.t)\n",
    "         zs.append(myode.y)\n",
    "\n",
    "      z_final = zs[t_last]\n",
    "      #\n",
    "      h1_f = z_final[0]\n",
    "      h1_r_f = z_final[1]\n",
    "      h2_f = z_final[2]\n",
    "      \n",
    "      p = np.zeros(10)\n",
    "      for n in range(0,10):\n",
    "        p[n] = epsilon * z_final[(n + 3)] ## the pathogen for the next season depends on epsilon, or pathogen clearance\n",
    "      \n",
    "      z0[3:13] = p\n",
    "  \n",
    "      t_extra = np.full((t_last + 1,),i)\n",
    "      general_ts.append(ts + t_extra)\n",
    "      general_zs.append(zs)\n",
    "\n",
    "      t_i = t_i + 1  #discrete time increase\n",
    "      t_f = t_f + 1\n",
    "    return general_ts,general_zs\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Soil quality and cash yield in crop rotation patterns (section 2.1 and 2.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22.5105737509\n",
      "42.2727701109\n",
      "0.532507656628\n",
      "(10,)\n",
      "()\n",
      "22.5331140364\n",
      "42.1815807214\n",
      "0.534193210664\n",
      "(10,)\n",
      "()\n",
      "22.8807020778\n",
      "42.1815464342\n",
      "0.542433931706\n",
      "(10,)\n",
      "()\n",
      "22.9857914194\n",
      "42.1073700616\n",
      "0.545885230681\n",
      "(10,)\n",
      "()\n",
      "22.8774804559\n",
      "41.8274496819\n",
      "0.546948968438\n",
      "(10,)\n",
      "()\n",
      "22.4737396447\n",
      "41.1400026512\n",
      "0.546274628012\n",
      "(10,)\n",
      "()\n",
      "21.7256524989\n",
      "39.8335882959\n",
      "0.545410379239\n",
      "(10,)\n",
      "()\n",
      "20.3607835205\n",
      "37.3992633279\n",
      "0.544416699921\n",
      "(10,)\n",
      "()\n",
      "17.7986545195\n",
      "32.8550557236\n",
      "0.54173259267\n",
      "(10,)\n",
      "()\n",
      "13.2015867511\n",
      "24.7963062055\n",
      "0.532401344042\n",
      "(10,)\n",
      "()\n",
      "6.32861383172\n",
      "12.7959342247\n",
      "0.494580053366\n",
      "(10,)\n",
      "()\n",
      "24.2161368932\n",
      "42.0520852555\n",
      "0.575860548794\n",
      "(10,)\n",
      "()\n",
      "23.1553531031\n",
      "41.8644745565\n",
      "0.553102680696\n",
      "(10,)\n",
      "()\n",
      "22.8460470461\n",
      "41.697724583\n",
      "0.547896732366\n",
      "(10,)\n",
      "()\n",
      "22.6598867153\n",
      "41.2316765733\n",
      "0.549574710478\n",
      "(10,)\n",
      "()\n",
      "21.8997688639\n",
      "39.8607333341\n",
      "0.549407073882\n",
      "(10,)\n",
      "()\n",
      "20.3920702735\n",
      "37.2655992426\n",
      "0.547208972562\n",
      "(10,)\n",
      "()\n",
      "17.7702991512\n",
      "32.6786884674\n",
      "0.543788627532\n",
      "(10,)\n",
      "()\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-11-e2d6d684af2e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     58\u001b[0m     \u001b[0mh1_seq\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mh1_r_seq\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mh2_seq\u001b[0m\u001b[0;34m=\u001b[0m \u001b[0mget_h_seq\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     59\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 60\u001b[0;31m     \u001b[0mts\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mzs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minfection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mh1_seq\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mh1_r_seq\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mh2_seq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     61\u001b[0m     \u001b[0mdeltas\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgetdeltas\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseq\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mzs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     62\u001b[0m     \u001b[0minit_soil\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfull\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mq0\u001b[0m\u001b[0;34m)\u001b[0m    \u001b[0;31m#init soil is 1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-8-e1825c7ed73d>\u001b[0m in \u001b[0;36minfection\u001b[0;34m(h1_seq, h1_r_seq, h2_seq)\u001b[0m\n\u001b[1;32m     52\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     53\u001b[0m       \u001b[0;32mwhile\u001b[0m \u001b[0mmyode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msuccessful\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mmyode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mt_end\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 54\u001b[0;31m          \u001b[0mmyode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mintegrate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmyode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mt_step\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     55\u001b[0m          \u001b[0mts\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmyode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     56\u001b[0m          \u001b[0mzs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmyode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/bargues/anaconda2/lib/python2.7/site-packages/scipy/integrate/_ode.pyc\u001b[0m in \u001b[0;36mintegrate\u001b[0;34m(self, t, step, relax)\u001b[0m\n\u001b[1;32m    406\u001b[0m             self._y, self.t = mth(self.f, self.jac or (lambda: None),\n\u001b[1;32m    407\u001b[0m                                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_y\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 408\u001b[0;31m                                 self.f_params, self.jac_params)\n\u001b[0m\u001b[1;32m    409\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mSystemError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    410\u001b[0m             \u001b[0;31m# f2py issue with tuple returns, see ticket 1187.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/bargues/anaconda2/lib/python2.7/site-packages/scipy/integrate/_ode.pyc\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, f, jac, y0, t0, t1, f_params, jac_params)\u001b[0m\n\u001b[1;32m   1207\u001b[0m         \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcall_args\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1208\u001b[0m                \u001b[0;34m[\u001b[0m\u001b[0mjac\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcall_args\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mf_params\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mjac_params\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1209\u001b[0;31m         \u001b[0my1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mistate\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrunner\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1210\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mistate\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1211\u001b[0m             warnings.warn('lsoda: ' +\n",
      "\u001b[0;32m<ipython-input-8-e1825c7ed73d>\u001b[0m in \u001b[0;36mHP_evol\u001b[0;34m(t, init)\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mp_type\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m             \u001b[0mpdot\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mp_type\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mq_matrix\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mp_type\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'float16'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msigma\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mw_matrix\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mp_type\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mh\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'float16'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mp\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mp_type\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mdp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     13\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m         \u001b[0mydot\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhdot\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpdot\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/bargues/anaconda2/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc\u001b[0m in \u001b[0;36msum\u001b[0;34m(a, axis, dtype, out, keepdims)\u001b[0m\n\u001b[1;32m   1722\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1723\u001b[0m         return _methods._sum(a, axis=axis, dtype=dtype,\n\u001b[0;32m-> 1724\u001b[0;31m                             out=out, keepdims=keepdims)\n\u001b[0m\u001b[1;32m   1725\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1726\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mproduct\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/bargues/anaconda2/lib/python2.7/site-packages/numpy/core/_methods.pyc\u001b[0m in \u001b[0;36m_sum\u001b[0;34m(a, axis, dtype, out, keepdims)\u001b[0m\n\u001b[1;32m     30\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     31\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_sum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 32\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mumr_sum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     33\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     34\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_prod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "def getdeltas(seq,zs):  ## section 2.3, delta as the ratio of host at the end of season and initial host density\n",
    "    zs = np.array(zs)\n",
    "    zz = zs[:,t_last]\n",
    "    zz_all = []\n",
    "    for j in range(0,l):\n",
    "        zz_1 = zz[j]\n",
    "        zz_2 = zz_1[0:3]\n",
    "        zz_all.append(zz_2)\n",
    "    zz_all = np.array(zz_all)\n",
    "    x = np.sum(zz_all,axis = 1)    \n",
    "    deltas = x/H_dens \n",
    "    return deltas\n",
    "\n",
    "\n",
    "\n",
    "def fitness_function(init_cash,init_soil,seq,deltas): ## fitness of sequences\n",
    "# \n",
    "    #print(\"deltas\",deltas)\n",
    "\n",
    "    soil_vector = []\n",
    "    cash_vector = []\n",
    "    n_h2 = 0 #n of consecutive cover crops\n",
    "    n_h1_hr = 0 #n of consecutive cash crops\n",
    "    \n",
    "    soil_seq = np.zeros(l + 1) #vector that will have soil changes\n",
    "    soil_seq[0] = init_soil\n",
    "\n",
    "    cash_seq = np.zeros(l + 1) #vector that will have cash changes\n",
    "    cash_seq[0] = init_cash\n",
    "\n",
    "    for t in range(0,l):\n",
    "        d = deltas[t]\n",
    "        \n",
    "        if seq[t] == 'h2': #soil quality and cash yield for cv crops\n",
    "            n_h2 = n_h2 + 1\n",
    "            n_h1_hr = 0\n",
    "            soil_seq[t + 1] = max(0.01, min(1.99,((K * soil_seq[t] * math.exp(beta_h2)) / (K + soil_seq[t] * (math.exp(beta_h2) - 1)))))\n",
    "            cash_seq[t + 1] = cash_seq[t] + y_h2 * d * soil_seq[t]\n",
    "        elif seq[t]=='h1': \n",
    "            n_h1_hr = n_h1_hr + 1\n",
    "            n_h2 = 0\n",
    "            soil_seq[t + 1] = max(0.01, min(1.99,((K * soil_seq[t] * math.exp(beta_h1)) / (K + soil_seq[t] * (math.exp(beta_h1) - 1)))))\n",
    "            cash_seq[t + 1] = cash_seq[t] + y_h1 * d * soil_seq[t]\n",
    "        elif seq[t]=='hr': \n",
    "            n_h1_hr = n_h1_hr + 1\n",
    "            n_h2 = 0\n",
    "            soil_seq[t + 1] = max(0.01, min(1.99,((K * soil_seq[t] * math.exp(beta_hr)) / (K + soil_seq[t] * (math.exp(beta_hr) - 1)))))\n",
    "            cash_seq[t + 1] = cash_seq[t] + y_hr * d * soil_seq[t]\n",
    "\n",
    "    fitness = cash_seq[l]\n",
    "    return cash_seq,soil_seq,fitness\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "for j in range(0,num_seq):\n",
    "    seq = list_seq[j]\n",
    "    h1_seq,h1_r_seq,h2_seq= get_h_seq(seq)\n",
    "    \n",
    "    ts,zs = infection(h1_seq,h1_r_seq,h2_seq)\n",
    "    deltas = getdeltas(seq,zs)\n",
    "    init_soil = np.full(1,q0)    #init soil is 1\n",
    "    init_cash = np.full(1,y0)   #init cash is 0\n",
    "    cash_seq,soil_seq,fitness = fitness_function(init_cash,init_soil,seq,deltas)\n",
    "    z = np.array(zs).T\n",
    "\n",
    "    h1,h1_r,h2,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10 = z\n",
    "    \n",
    "    \n",
    "    p_sum1 = p1 + p2 + p3 + p4 + p5\n",
    "    final_p5 = p5[-1,-1] \n",
    "    final_pathogen = p_sum1[-1,-1]    \n",
    "    freq_p5 = final_p5 / final_pathogen  \n",
    "    \n",
    "    \n",
    "    soil_seq = soil_seq[1:11]\n",
    "    cash_seq = cash_seq[1:11]\n",
    "\n",
    "    \n",
    "    df['h1_seq'] = np.append(np.nan,h1_seq)\n",
    "    df['h1_r_seq'] = np.append(np.nan,h1_r_seq)\n",
    "    df['h2_seq'] = np.append(np.nan, h2_seq)\n",
    "    df['deltas'] = np.append(np.nan,deltas)\n",
    "    df['soil'] = np.append(np.nan,soil_seq)\n",
    "    df['cash'] = np.append(np.nan,cash_seq)\n",
    "    df['freq p5'] = np.append(h1_seq,freq_p5)\n",
    "    \n",
    "    df2['seq and fitness'] = np.append(h1_seq,cash_seq[9])\n",
    "    \n",
    "    df.to_csv('out{0}.csv'.format(j+1))\n",
    "    df2.to_csv('h1_seq{0}.csv'.format(j+1), index = False)\n"
   ]
  },
  {
   "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}