Mating Sandbox.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#! /usr/bin/env python\n",
    "\n",
    "\"\"\"\n",
    "Simulation to analyze mating dynamics of cowbirds (or more generally, monogamy in the absence of parental care) \n",
    "written by Ammon Perkes (perkes.ammon@gmail.com) at University of Pennsylvania\n",
    "2016\n",
    "\"\"\" \n",
    "\n",
    "import sys, os\n",
    "import numpy as np\n",
    "from  matplotlib import pyplot as plt\n",
    "from matplotlib.widgets import Slider\n",
    "import plotly as py\n",
    "import networkx as nx               # Necessary for all the network analysis and plotting\n",
    "#from plotly.graph_objs import *    # For playing with nice network graphs, not necessary\n",
    "\n",
    "#External .py dependencies, contatining Strategies and Stats\n",
    "import SimStrategies,SimStats\n",
    "\n",
    "\n",
    "# Default parameters. These can be defined invdividually as well through the parameter class. \n",
    "\n",
    "N_MALES = 6\n",
    "N_FEMALES = 6\n",
    "RES_M = 1.0       # Resource limitation for a male\n",
    "RES_F = 1.0       # Resource limitation for a female\n",
    "RESS_M = [RES_M] * N_MALES     # Resource limitation vector for all males\n",
    "RESS_F = [RES_F] * N_FEMALES   # Resource limitation vector for all females\n",
    "Q_MALE = 1.0      # Default Quality for a male\n",
    "Q_FEMALE = 1.0    # Default Quality for a female\n",
    "Q_MALES = [Q_MALE] * N_MALES     # Quality vector for maless\n",
    "Q_FEMALES = [Q_FEMALE] * N_FEMALES # Quality vector for females \n",
    "TURNS = 1000       # Number of turns per trial\n",
    "TRIALS = 10        # Number of trials per simulation\n",
    "STRAT_M = 4        # Default strategy for males\n",
    "STRAT_F = 7        # Default strategy for females \n",
    "STRATS_M = [STRAT_M] * N_MALES       # Strategy vector for males\n",
    "STRATS_F = [STRAT_F] * N_FEMALES     # Strategy vector for females\n",
    "ALPHA = 1.0        # Alpha value scales the adjacency matrix \n",
    "ALPHAS = [ALPHA] * N_FEMALES         # Alpha vector for females\n",
    "BETA = 1.0         # scales the investment -> benefit function\n",
    "IOTA = 1.5         # scales the investment to adjacency function\n",
    "KAPPA = 0.0        # defines adjacency kost\n",
    "OMEGA = 0.5        # maximum cost of crowding\n",
    "PHI = 1.0          # Fleeing multiplier when determining adjacency\n",
    "RHO = 0.0          # base risk of individual association\n",
    "TAU = 2.0          # scales the rate at which crowding cost increases\n",
    "SHIFT = .1         # amount by which investment is shifted\n",
    "DECAY = .01         # decay constant \n",
    "BUMP = .01         # random encounter rate\n",
    "MIN_ADJ = .01      # Minimum adjacency\n",
    "MAX_ADJ = 1.0      # Maximum adjacency\n",
    "\n",
    "## Functions that define how matrices covert, these are two of the central mechanisms of the simulation\n",
    "def adjacency_to_reward(history, turn, params):\n",
    "    ## This provides males and females with a reward function on which to base decisions\n",
    "    # Noteably, this is different from the original implementation, in which the reward *is* the female response\n",
    "    # Because it needs to update immediately (and really everything should...) it needs the current turn also\n",
    "    \n",
    "    alpha = params.alpha # Adjacency benefit multiplier. \n",
    "    beta = params.beta # Investment benefit multiplier. \n",
    "    omega = params.omega ## This defines the maximum cost of crowding\n",
    "    tau = params.tau  ## This defines how quickly crowding cost decreases\n",
    "    kappa = params.kappa ## This defines base adjacency cost\n",
    "    rho = params.rho ## This defines risk from individual (also a cost)\n",
    "    n_birds = history.n_males + history.n_females\n",
    "    r_function = params.r_function\n",
    "    reward = np.zeros(np.shape(history.reward_matrix2[0]))\n",
    "\n",
    "    #current_adjacency = history.adjacency_matrix[history.current_turn]   # This way it gets the most recent adjacency\n",
    "    current_adjacency = turn.adjacency # This way it gets the most recent adjacency\n",
    "    current_investment = turn.invest2\n",
    "    base_adjacency_cost = kappa  #\n",
    "    individual_cost = rho # This is an invariant risk of associating with any individual ()\n",
    "    \n",
    "    # Presumably, costs/benefits should be different for males and females                  \n",
    "    if r_function == 0:\n",
    "        for i in range(n_birds):\n",
    "            for j in range(n_birds):\n",
    "                crowding_cost = omega / (tau ** (n_birds - current_adjacency[:,j].sum()) + .0001)\n",
    "                adjacency_cost = current_adjacency[i,j] * base_adjacency_cost + individual_cost\n",
    "                #adjacency_benefit = (current_adjacency[i,j] / (current_adjacency[:,j].sum() + .0001)) * alpha\n",
    "                adjacency_benefit = current_adjacency[i,j] ** alpha\n",
    "                investment_benefit = current_investment[j,i] ** beta\n",
    "                reward[i,j] = (adjacency_benefit * investment_benefit * history.quality_vector[j]\n",
    "                                   - crowding_cost - adjacency_cost)\n",
    "            reward[i,i] = 0\n",
    "    \n",
    "    #If desired, set reward for same-sex pairs to 0\n",
    "    #reward[0:n_males,0:n_males] = 0\n",
    "    #reward[n_males:,n_males:] = 0\n",
    "      \n",
    "    return reward\n",
    "\n",
    "# Translate investment into adjacency\n",
    "def investment_to_adjacency(history, params):\n",
    "    # Here, investment cuts a fraction of the distance, using new investment matrix\n",
    "    \n",
    "    iota = params.iota ## This is a parameter that determines how investment scales\n",
    "    phi = params.phi ## This is a parameter that determines how investment in fleeing scales (get it? phleeing)\n",
    "    shift_constant = params.shift_constant ## This scales down investment to make everything more gradual\n",
    "    decay_constant = params.decay_constant\n",
    "    min_adjacency = params.min_adjacency\n",
    "    max_adjacency = params.max_adjacency\n",
    "    previous_adjacency = history.adjacency_matrix[history.current_turn - 1]\n",
    "    a_function = params.a_function\n",
    "    ### The following line adds slippage, this is important to avoid crowding when no one is fleeing.\n",
    "    if a_function == 0:\n",
    "        encounter_rate = params.bump_constant\n",
    "        random_encounter = (1 - previous_adjacency) * encounter_rate\n",
    "        previous_adjacency = previous_adjacency - decay_constant * previous_adjacency \n",
    "\n",
    "        previous_investment = history.invest_matrix2[history.current_turn - 1]\n",
    "\n",
    "        shift_matrix = previous_investment ** iota\n",
    "        shift_matrix[previous_investment < 0] = - (phi + (1 - phi) * previous_investment[previous_investment < 0] ** iota) \n",
    "        current_adjacency = previous_adjacency + (previous_adjacency) * shift_matrix * shift_constant + random_encounter\n",
    "        current_adjacency[current_adjacency < min_adjacency] = min_adjacency\n",
    "        current_adjacency[current_adjacency > max_adjacency] = max_adjacency\n",
    "    return current_adjacency\n",
    "\n",
    "# Parameter object containing all necessary parameters. \n",
    "class Parameters(object):\n",
    "    def __init__(self, r_function = 0, a_function = 0, n_turns = TURNS, n_trials = TRIALS,\n",
    "                 strat_f = STRAT_F, strat_m = STRAT_M, strats_f = None, strats_m = None, \n",
    "                 res_m = RES_M, res_f = RES_F, ress_m = None, ress_f = None,\n",
    "                 n_males = N_MALES, n_females = N_FEMALES,\n",
    "                 q_male = Q_MALE, q_female = Q_FEMALE, q_males = None, q_females = None,\n",
    "                 alpha = ALPHA, alphas = None, beta = BETA, iota = IOTA, kappa = KAPPA, omega = OMEGA, phi = PHI, rho = RHO, tau = TAU,\n",
    "                 shift_constant = SHIFT, decay_constant = DECAY, bump_constant = BUMP, min_adjacency = MIN_ADJ, max_adjacency = MAX_ADJ):\n",
    "        self.r_function = r_function\n",
    "        self.a_function = a_function\n",
    "        self.alpha = alpha\n",
    "        self.aphas = alphas\n",
    "        self.beta = beta\n",
    "        self.iota = iota\n",
    "        self.kappa = kappa\n",
    "        self.omega = omega\n",
    "        self.phi = phi\n",
    "        self.rho = rho\n",
    "        self.tau = tau\n",
    "        self.shift_constant = shift_constant\n",
    "        self.decay_constant = decay_constant\n",
    "        self.bump_constant = bump_constant\n",
    "        self.min_adjacency = min_adjacency\n",
    "        self.max_adjacency = max_adjacency\n",
    "        self.strat_m = strat_m\n",
    "        self.strat_f = strat_f\n",
    "        self.strats_m = strats_m\n",
    "        self.strats_f = strats_f\n",
    "        self.res_m = res_m    # Resource limitation for a male\n",
    "        self.res_f = res_f    # Resource limitations for a female (NOTE: This might be unnecessary also....)\n",
    "        self.ress_m = ress_m\n",
    "        self.ress_f = ress_f\n",
    "        self.n_males = n_males    # Number of males per trial \n",
    "        self.n_females = n_females # Number of females per trial \n",
    "        self.q_male = q_male\n",
    "        self.q_female = q_female\n",
    "        self.q_males = q_males\n",
    "        self.q_females = q_females\n",
    "        self.n_turns = n_turns # Number of turns per trial\n",
    "        self.n_trials = n_trials  # Number of trials per simulation\n",
    "        \n",
    "        if alphas == None:\n",
    "            self.alphas = [self.alpha] * self.n_females\n",
    "        if strats_f == None:\n",
    "            self.strats_f = [self.strat_f] * self.n_females\n",
    "        if strats_m == None:\n",
    "            self.strats_m = [self.strat_m] * self.n_males\n",
    "        if ress_f == None:\n",
    "            self.ress_f = [self.res_f] * self.n_females\n",
    "        if ress_m == None:\n",
    "            self.ress_m = [self.res_m] * self.n_males\n",
    "        if self.q_females == None:\n",
    "            self.q_females = [self.q_female] * self.n_females\n",
    "        if self.q_males == None:\n",
    "            self.q_males = [self.q_male] * self.n_males\n",
    "\n",
    "# Object containing all the birds. Cute, eh? \n",
    "class Aviary(object):\n",
    "    def __init__(self, params = Parameters()):\n",
    "# Initialize some parameters:\n",
    "        self.params = params\n",
    "        self.n_males = params.n_males\n",
    "        self.n_females = params.n_females\n",
    "        self.strats_m = params.strats_m\n",
    "        self.strats_f = params.strats_f\n",
    "        self.ress_m = params.ress_m\n",
    "        self.ress_f = params.ress_f\n",
    "        self.q_males = params.q_males\n",
    "        self.q_females = params.q_females\n",
    "# Build the male and female lists in the aviary. \n",
    "        self.males = [Male_bird(num, params.strats_m[num], params.ress_m[num], params.q_males[num]) \n",
    "                      for num in range(params.n_males)]\n",
    "        self.females = [Female_bird(num, params.strats_f[num], params.ress_f[num], params.alphas[num], params.q_females[num]) \n",
    "                        for num in range(params.n_females)]\n",
    "\n",
    "        ## This is a big important function. It is how the aviary of birds responds every turn\n",
    "    def respond(self,history):\n",
    "        # Initialize Turn\n",
    "        turn = Turn(history.current_turn,history.n_males,history.n_females)\n",
    "            \n",
    "        # Save matrices to the turn\n",
    "        #turn.invest = self.mrespond(history)\n",
    "        ## mrespond & frespond also \n",
    "        #turn.reward = self.frespond(history)\n",
    "        turn.adjacency = self.update_adjacency2(history, turn)\n",
    "        turn.invest2 = self.update_invest2(history, turn)\n",
    "        turn.reward2 = self.update_reward2(history, turn)\n",
    "        return turn\n",
    "\n",
    "    def mrespond(self,history):\n",
    "        invest = np.zeros([self.n_males,self.n_females])\n",
    "        for m in range(self.n_males):\n",
    "            invest[m] = self.males[m].respond(history)[m]\n",
    "        return invest\n",
    "\n",
    "    def frespond(self,history):\n",
    "        reward = np.zeros([self.n_males,self.n_females])\n",
    "        for f in range(self.n_females):\n",
    "            reward[:,f] = self.females[f].respond(history)[:,f]   \n",
    "        return reward\n",
    "    \n",
    "    def update_adjacency(self,history, turn):\n",
    "        ## Here, investment adds to adjacency, with a hard cap at 1, and optional slippage\n",
    "        slippage = DECAY\n",
    "        min_adjacency = MIN_ADJ\n",
    "        previous_adjacency = history.adjacency_matrix[history.current_turn - 1]\n",
    "        previous_investment = history.invest_matrix[history.current_turn - 1]\n",
    "        investment_adj = get_adjacency(previous_investment)\n",
    "        current_adjacency = previous_adjacency + investment_adj\n",
    "        \n",
    "        current_adjacency = current_adjacency - slippage\n",
    "        #current_adjacency = current_adjacency * (investment_adj + 1) #update based on investment\n",
    "        \n",
    "        current_adjacency[current_adjacency <= min_adjacency] = min_adjacency #limit max distance\n",
    "        current_adjacency[current_adjacency > 1.0] = 1.0\n",
    "        return current_adjacency\n",
    "    \n",
    "    def update_adjacency1(self,history, turn):\n",
    "        # Here investment cuts a fraction of the distance using old investment matrix\n",
    "        # This means full inveestment will always get you to 1, but it also means you get diminishing returns\n",
    "        params = self.params\n",
    "        alpha = params.alpha ## This is how investment scales up or down\n",
    "        iota = params.iota ## This is a parameter that determines how investment scales\n",
    "        phi = params.phi ## This is a parameter that determines how ivnestment in fleeing scales (get it? phleeing)\n",
    "        previous_adjacency = history.adjacency_matrix[history.current_turn - 1]\n",
    "        \n",
    "        ### The following line adds slippage, this is important to avoid crowding when no one is fleeing.\n",
    "        previous_adjacency = previous_adjacency - previous_adjacency * .01 \n",
    "        previous_investment = history.invest_matrix[history.current_turn - 1]\n",
    "        investment_adj = get_adjacency(previous_investment)\n",
    "        shift_matrix = investment_adj ** alpha\n",
    "        shift_matrix[investment_adj < 0] = - (phi + (1 - phi) * investment_adj[investment_adj < 0] ** iota) \n",
    "        current_adjacency = previous_adjacency + (1 - previous_adjacency) * shift_matrix\n",
    "        return current_adjacency\n",
    "        \n",
    "    def update_adjacency2(self,history,turn):\n",
    "        \n",
    "        current_adjacency = investment_to_adjacency(history, self.params)\n",
    "        return current_adjacency\n",
    "    \n",
    "    def update_invest1(self,history,turn):\n",
    "        ## This provides both males and females an opportunity to return investment\n",
    "        previous_invest = history.invest_matrix2[turn.n - 1]\n",
    "        male_invest = np.zeros(np.shape(previous_invest))\n",
    "        female_invest = np.zeros(np.shape(previous_invest))\n",
    "        for m in range(history.n_males):\n",
    "            male_invest[m] = self.males[m].respond2(history)[m]\n",
    "        for f in range(history.n_females):\n",
    "            female_invest[f + history.n_males] = self.females[f].respond2(history)[f + history.n_males]\n",
    "        new_invest = male_invest + female_invest\n",
    "        return new_invest\n",
    "\n",
    "    def update_invest2(self,history,turn):\n",
    "        ## This provides both males and females an opportunity to return investment and competition\n",
    "        ## For it to work, I need to make sure everything works cleanly\n",
    "        previous_invest = history.invest_matrix2[turn.n - 1]\n",
    "        male_invest = np.zeros(np.shape(previous_invest))\n",
    "        female_invest = np.zeros(np.shape(previous_invest))\n",
    "        history.invest_matrix2[turn.n] = previous_invest \n",
    "        for m in range(history.n_males):\n",
    "            history.invest_matrix2[turn.n] = self.males[m].respond2(history)\n",
    "        for f in range(history.n_females):\n",
    "            history.invest_matrix2[turn.n] = self.females[f].respond2(history)\n",
    "        #new_invest = male_invest + female_invest\n",
    "        return history.invest_matrix2[turn.n]\n",
    "    \n",
    "    def update_reward2(self,history,turn):\n",
    "        new_reward2 = adjacency_to_reward(history,turn, self.params)\n",
    "        return new_reward2\n",
    "        \n",
    "        \n",
    "class Male_bird(object):\n",
    "    def __init__(self, num, strategy = 0, resources = 1, quality = 1.0): #Convention: Males have even strategies, females odd. \n",
    "        self.num = num\n",
    "        self.strategy = strategy\n",
    "        self.resources = resources\n",
    "        self.quality = quality\n",
    "\n",
    "##      Seed self investment, and normalize for resources\n",
    "#   Functions to adjust and get info. \n",
    "### NOTE: This is where males apply strategy, strategies are saved externally\n",
    "    def respond(self,history):\n",
    "        new_investment = SimStrategies.choose(self.strategy, self.resources, history, self.num)\n",
    "        return new_investment \n",
    "    \n",
    "    def respond2(self,history):\n",
    "        new_investment = SimStrategies.choose('M0', self.resources, history, self.num)\n",
    "        return new_investment\n",
    "\n",
    "#Class of female cowbirds: \n",
    "#Includes: resources, response matrix, reward matrix\n",
    "#Also includes choice function determining response\n",
    "class Female_bird(object):\n",
    "    def __init__(self, num, strategy = 1, resources = 1, alpha = ALPHA, quality = 1.0):\n",
    "        self.num = num\n",
    "        self.strategy = strategy\n",
    "        self.resources = resources\n",
    "        self.alpha = alpha\n",
    "        self.quality = quality\n",
    "##      Seed self  investment, NOTE that for females this is which males are investing in them, not vice versa\n",
    "#   Functions to adjust and get info. \n",
    "    def respond(self,history):\n",
    "### NOTE: This is where strategy is applied\n",
    "        resources = self.resources\n",
    "        new_response = SimStrategies.choose(self.strategy,self.resources, history, self.num, self.alpha)\n",
    "        return new_response \n",
    "    \n",
    "    def respond2(self,history):\n",
    "        new_investment = SimStrategies.choose('F0', self.resources, history, self.num)\n",
    "        return new_investment\n",
    "\n",
    "#Array type object containing cowbirds, allowing cleverness: \n",
    "#Keep history of investment, reward for both males and females\n",
    "\n",
    "#NOTE: I want to change this so that we have a full matrix of investment, where birds can invest in themselves, in same-sex birds, with negative investment, etc. This is a trivial change as it relates to history, but it will end up being quite a lot of changes down the line. \n",
    "class History(object):\n",
    "    def __init__(self, params = Parameters()):\n",
    "## Initialize the matrix for the whole sim (save a bunch of time)\n",
    "        n_males = params.n_males\n",
    "        n_females = params.n_females\n",
    "        n_turns = params.n_turns\n",
    "        self.invest_matrix = np.zeros([n_turns,n_males,n_females])\n",
    "        self.reward_matrix = np.zeros([n_turns,n_males,n_females])\n",
    "        self.adjacency_matrix = np.zeros([n_turns,n_males + n_females, n_males + n_females])\n",
    "        self.adjacency_matrix[0,:n_males,n_males:] = np.random.rand(n_males,n_females)\n",
    "        self.invest_matrix2 = np.zeros(np.shape(self.adjacency_matrix))\n",
    "        self.reward_matrix2 = np.zeros(np.shape(self.adjacency_matrix))       \n",
    "        self.n_turns = params.n_turns\n",
    "        self.n_males = params.n_males\n",
    "        self.n_females = params.n_females\n",
    "        self.current_turn = 0\n",
    "        self.quality_males = params.q_males\n",
    "        self.quality_females = params.q_females\n",
    "        self.quality_vector = np.append(params.q_males,params.q_females)\n",
    "        self.params = params\n",
    "\n",
    "    def record(self,turn): \n",
    "        #self.current_turn = turn.n\n",
    "        #self.invest_matrix[turn.n] = turn.invest\n",
    "        #self.reward_matrix[turn.n] = turn.reward\n",
    "        self.adjacency_matrix[turn.n] = turn.adjacency\n",
    "        self.invest_matrix2[turn.n] = turn.invest2\n",
    "        self.reward_matrix2[turn.n] = turn.reward2\n",
    "        self.advance()\n",
    "        \n",
    "    def initialize(self,initial_conditions = None): #set first investment conditions\n",
    "        self.invest_matrix[0] = np.random.random((self.n_males,self.n_females))\n",
    "        self.invest_matrix[0] = self.invest_matrix[0] / self.invest_matrix[0].sum(1) \n",
    "\n",
    "        self.reward_matrix[0] = np.random.random([self.n_males,self.n_females])\n",
    "        self.reward_matrix[0] = self.reward_matrix[0] / self.reward_matrix[0].sum(0)\n",
    "        \n",
    "        self.invest_matrix2[0] = np.random.random([self.n_males + self.n_females, self.n_males + self.n_females]) * (1 - np.identity(self.n_males + self.n_females))\n",
    "        self.invest_matrix2[0] = self.invest_matrix2[0] / self.invest_matrix2[0].sum(1)\n",
    "        self.adjacency_matrix[0] = np.random.random(np.shape(self.invest_matrix2[0])) * (1 - np.identity(self.n_males + self.n_females))\n",
    "        self.reward_matrix2[0] = self.update_reward_hist()\n",
    "        \n",
    "        if initial_conditions == None:\n",
    "            self.invest_matrix[0] = self.invest_matrix[0] / self.invest_matrix[0].sum(1).reshape(self.n_males,1)\n",
    "        else:\n",
    "            self.invest_matrix[0] = self.invest_matrix[0] * initial_conditions / self.invest_matrix[0].sum(1).reshape(self.n_males,1)\n",
    "\n",
    "        self.reward_matrix[-1] = self.reward_matrix[0] # This is for when the f_respond checks for previous\n",
    "        self.invest_matrix[-1] = self.invest_matrix[0]\n",
    "        self.invest_matrix2[-1] = self.invest_matrix2[0]\n",
    "        self.reward_matrix2[-1] = self.reward_matrix2[0]\n",
    "        self.adjacency_matrix[-1] = self.adjacency_matrix[0]\n",
    "        \n",
    "        \n",
    "    def advance(self):\n",
    "        self.current_turn = self.current_turn + 1\n",
    "    \n",
    "    def update_reward_hist(self):\n",
    "        ## This provides males and females with a reward function on which to base decisions\n",
    "        # Noteably, this is different from the original implementation, in which the reward *is* the female response\n",
    "        # Because it needs to update immediately (and really everything should...) it needs the current turn also\n",
    "        alpha = ALPHA # This should be a different parameter. Not everything is alpha\n",
    "        \n",
    "        new_reward2 = np.empty_like(self.reward_matrix2[0])\n",
    "        #current_adjacency = history.adjacency_matrix[history.current_turn]\n",
    "        current_adjacency = self.adjacency_matrix[self.current_turn]   # This way it gets the most recent adjacency\n",
    "        n_birds = self.n_females + self.n_males\n",
    "        omega = 0.5 ## This defines the maximum cost of crowding\n",
    "        tau = 2.0  ## This defines how quickly crowding changes\n",
    "        base_adjacency_cost = 0.0\n",
    "        individual_risk = 0.0\n",
    "        for i in range(n_birds):\n",
    "            for j in range(n_birds):\n",
    "                crowding_cost = omega / tau ** (1 - current_adjacency[:,j].sum())\n",
    "                adjacency_cost = current_adjacency[i,j] * base_adjacency_cost\n",
    "                new_reward2[i,j] = (current_adjacency[i,j] ** alpha * self.quality_vector[j] \n",
    "                                   - crowding_cost - base_adjacency_cost - individual_risk)\n",
    "            new_reward2[i,i] = 0\n",
    "        return new_reward2         \n",
    "\n",
    "#Object containing turn, I don't actually use this at all\n",
    "class Turn(object):\n",
    "    def __init__(self, n, n_males, n_females, previous_turn = None):\n",
    "## Initialize based on last turn if desired, otherwise start blank\n",
    "        self.n = n\n",
    "        if previous_turn != None:\n",
    "            self.invest = previous_turn.invest\n",
    "            self.reward = previous_turn.reward\n",
    "            self.adjacency = previous_turn.adjacency\n",
    "            self.invest2 = previous_turn.invest2\n",
    "            self.reward2 = previous_turn.reward2\n",
    "        else:\n",
    "            self.invest = np.zeros([n_males,n_females])\n",
    "            self.reward = np.zeros([n_males,n_females])\n",
    "            self.adjacency = np.zeros([n_males+n_females,n_females+n_males])\n",
    "            self.invest2 = np.zeros(np.shape(self.adjacency))\n",
    "            self.reward2 = np.zeros(np.shape(self.adjacency))\n",
    "    def change_invest(self,male,female,amount):\n",
    "        self.invest[male,female] += amount\n",
    "    def change_reward(self,male,female,amount):\n",
    "        self.reward[male,female] += amount\n",
    "    def set_invest(self,male,female,amount):\n",
    "        self.invest[male,female] = amount\n",
    "    def set_reward(self,male,female,amount):\n",
    "        self.reward[male,female] = amount\n",
    "\n",
    "\n",
    "#Functions determining success: \n",
    "#This is not technically important for the simulation, but it determines which strategy is best, which is important\n",
    "\n",
    "def female_success(params):\n",
    "\n",
    "    return success\n",
    "\n",
    "def male_success(params):\n",
    "\n",
    "    return success\n",
    "\n",
    "def run_trial(params = Parameters(), initial_conditions = None):\n",
    "    n_males = params.n_males\n",
    "    n_females = params.n_females\n",
    "    n_turns = params.n_turns\n",
    "## Initialize full record...\n",
    "## initialize history\n",
    "    history = History(params)\n",
    "    history.initialize(initial_conditions)\n",
    "# Build an aviary using the given parameters\n",
    "    aviary = Aviary(params)\n",
    "# Initialize initial investment (based on male resources, if specified)\n",
    "    history.initialize()\n",
    "# Get first female response:\n",
    "    history.reward_matrix[0] = aviary.frespond(history)\n",
    "    history.advance()\n",
    "# For every turn, calculate response and record it in history.\n",
    "    for t in range(n_turns-1):\n",
    "        turn = aviary.respond(history)\n",
    "        history.record(turn)\n",
    "    return history\n",
    "       \n",
    "def run_simulation(trials = TRIALS, n_turns = TURNS, \n",
    "                   n_males = N_MALES, n_females = N_FEMALES, \n",
    "                   strat_males = None, strat_females = None, \n",
    "                   ress_males = None, ress_females = None, alpha = None, alphas = ALPHAS, params = Parameters()):\n",
    "    if alpha == None:\n",
    "        pass\n",
    "    else:\n",
    "        alphas = [alpha] * n_females\n",
    "    record = [0 for tr in range(trials)]\n",
    "    for tr in range(trials):\n",
    "        history = run_trial(params)\n",
    "        record[tr] = SimStats.get_stats(history) \n",
    "# For tidiness, stats is saved in a seperate file\n",
    "    return record\n",
    "\n",
    "# Function to plot the response curve of a male or female bird (which is contained in their class)\n",
    "def plot_response(bird):\n",
    "    strategy = bird.strategy\n",
    "    function = SimStrategies.function(strategy)\n",
    "    plt.plot(function)\n",
    "    plt.show()\n",
    "#NOTE: This is probably wrong, go back and check it.\n",
    "\n",
    "#Function plotting history and outcome in interesting ways\n",
    "def show_his_stats(history):\n",
    "    stats = SimStats.get_stats(history)\n",
    "    for stat in stats:\n",
    "        print stat\n",
    "        #stat.print_stat()\n",
    "        #stat.plot_stat()\n",
    "        raw_input('press enter to continue...')\n",
    "\n",
    "#Function for plotting a full simulation\n",
    "def show_record_stats(record):\n",
    "    r_stats = rec_stats(record)\n",
    "    for stat in r_stats:\n",
    "        stat.print_stat()\n",
    "        stat.plot_stat()\n",
    "        raw_input('press enter to continue...')\n",
    "\n",
    "# Mini function to get birds\n",
    "def get_birds():\n",
    "    print \"Default males: \" + str(N_MALES)\n",
    "    print \"Default females \" + str(N_FEMALES)\n",
    "    choice = raw_input(\"Would you like to change these? [y/n]\")\n",
    "    if choice[0].lower() == 'y':\n",
    "        try:\n",
    "            n_males = raw_input('How many males? ')\n",
    "            n_females = raw_input('How many females? ')\n",
    "            n_males = int(n_males.strip())\n",
    "            n_females = int(n_males.strip())\n",
    "        except:\n",
    "            print 'bad input, try again:'\n",
    "            return get_birds()\n",
    "        return(n_females,n_females)\n",
    "    else:\n",
    "        return (N_MALES,N_FEMALES)\n",
    "\n",
    "## Mini function to get strategy. Strategy will always be an array of ints matching the number of birds.\n",
    "#  If a single number is given, it fills the list (2 -> [2,2,2,2,2]\n",
    "#  If it matches teh number of birds, it will just be one to one, otherwise, any leftovers will be filled by default.\n",
    "def get_strategy(n_males = N_MALES, n_females = N_FEMALES):\n",
    "    print \"Default strategy: 1\"\n",
    "    print \"Enter either a single strategy (2) or multiple strategies, seperated by commas (2,3,4,5)\"\n",
    "    print \"If there are more birds than strategies, the rest will be filled by default.\"\n",
    "    print \"If there are more strategies than birds, the last will be cut off.\"\n",
    "    print \"....\"\n",
    "    print \"N Males: \" + str(n_males)\n",
    "    print \"What strategies would you like to use?\"\n",
    "# get input: \n",
    "    m_strat = raw_input()\n",
    "# clean it up and put it in the right form.\n",
    "    m_strat = m_strat.replace(' ','') \n",
    "    m_strat = m_strat.split(',')\n",
    "    m_strat = map(int, m_strat)\n",
    "# If a single number was given, it fills the list\n",
    "    if len(m_strat) == 1:\n",
    "        m_strat = m_strat * n_males\n",
    "# If too few numbers were given, the rest are filled by default\n",
    "    elif len(m_strat) < n_males:\n",
    "        for i in range(len(m_strat),n_males):\n",
    "            m_strat[i] = 1\n",
    "# if the right number of numbers is given, well done\n",
    "    elif len(m_strat) == n_males:\n",
    "        pass\n",
    "# IF too many numbers are given, the last ones are chopped off.\n",
    "    elif len(m_strat) > n_males:\n",
    "        del m_strat[n_males:]\n",
    "    print \"....\"\n",
    "    print \"N Females: \" + str(n_females)\n",
    "    print \"What strategies would you like to use?\"\n",
    "# get input: \n",
    "    f_strat = raw_input()\n",
    "# clean it up and put it in the right form.\n",
    "    f_strat = f_strat.replace(' ','') \n",
    "    f_strat = f_strat.split(',')\n",
    "    f_strat = map(int, f_strat)\n",
    "    if len(f_strat) == 1:\n",
    "        f_strat = f_strat * n_females\n",
    "    elif len(f_strat) < n_females:\n",
    "        for i in range(len(f_strat),n_females):\n",
    "            f_strat[i] = 1\n",
    "    elif len(f_strat) == n_females:\n",
    "        pass\n",
    "    elif len(f_strat) > n_females:\n",
    "        del f_strat[n_females:]\n",
    "    return m_strat,f_strat\n",
    "\n",
    "def get_resources(n_males = N_MALES, n_females = N_FEMALES):\n",
    "    print \"This will work just like strategy.\"\n",
    "    print \"Default males resources: \" + str(RES_M)\n",
    "    print 'What resources would you like the males to have'\n",
    "    m_res = raw_input()\n",
    "    m_res = m_res.replace(' ','')\n",
    "    m_res = m_res.split(',')\n",
    "    m_res = map(int, m_res)\n",
    "    if len(m_res) == 1:\n",
    "        m_res = m_res * n_males\n",
    "    elif len(m_res) < n_males:\n",
    "        for i in range(len(m_res), n_males):\n",
    "            m_res[i] = 1\n",
    "    elif len(m_res) == n_males:\n",
    "        pass\n",
    "    elif len(m_res) > n_males:\n",
    "        del m_res[n_males:]\n",
    "    print \"....\"\n",
    "    print \"N Females: \" + str(n_females)\n",
    "    print 'What resources would you like the females to have?'\n",
    "# get input: \n",
    "    f_res = raw_input()\n",
    "# clean it up and put it in the right form.\n",
    "    f_res = f_res.replace(' ','') \n",
    "    f_res = f_res.split(',')\n",
    "    f_res = map(int, f_res)\n",
    "    if len(f_res) == 1:\n",
    "        f_res = f_res * n_females\n",
    "    elif len(f_res) < n_females:\n",
    "        for i in range(len(f_res),n_females):\n",
    "            f_res[i] = 1\n",
    "    elif len(f_res) == n_females:\n",
    "        pass\n",
    "    elif len(f_res) > n_females:\n",
    "        del f_res[n_females:]\n",
    "    return m_res,f_res\n",
    "\n",
    "# function to set up a run custon simulations\n",
    "def build_simulation():\n",
    "    print \"This will help you build a simulation. Enter all values as integers (i.e., 11)\"\n",
    "    n_turns = raw_input(\"How many turns per trial would you like? \")\n",
    "    n_turns = int(turns.strip())\n",
    "    trials = raw_input(\"How many trials in the simulation would you like? \")\n",
    "    trials = int(trials.strip())\n",
    "    n_males,n_females = get_birds()\n",
    "    stat_males,strat_females = get_strategy()\n",
    "    res_males,res_females = get_res()\n",
    "    run_simulation(trials, n_turns, n_males, n_females, strat_males, strat_females, res_males, res_females)\n",
    "\n",
    "### Plotting Function ###\n",
    "# Plot history (with a fancy slider)\n",
    "def plot_history(history):\n",
    "    # generate a five layer data \n",
    "    data = history.invest_matrix\n",
    "    # current layer index start with the first layer \n",
    "    idx = 0\n",
    "\n",
    "    # figure axis setup \n",
    "    fig, ax = plt.subplots()\n",
    "    fig.subplots_adjust(bottom=0.15)\n",
    "\n",
    "    # display initial image \n",
    "    im_h = ax.imshow(data[idx, :, :], vmin=0., vmax=1., cmap='hot', interpolation='nearest')\n",
    "\n",
    "    # setup a slider axis and the Slider\n",
    "    ax_depth = plt.axes([0.23, 0.02, 0.56, 0.04])\n",
    "    slider_depth = Slider(ax_depth, 'depth', 0, data.shape[0]-1, valinit=idx)\n",
    "\n",
    "    # update the figure with a change on the slider \n",
    "    def update_depth(val):\n",
    "        idx = int(round(slider_depth.val))\n",
    "        im_h.set_data(data[idx, :, :])\n",
    "\n",
    "    slider_depth.on_changed(update_depth)\n",
    "\n",
    "    plt.show()  \n",
    "    \n",
    "def plot_history2(history_matrix):\n",
    "    # generate a five layer data \n",
    "    data = history_matrix\n",
    "    # current layer index start with the first layer \n",
    "    idx = 0\n",
    "\n",
    "    # figure axis setup \n",
    "    fig, ax = plt.subplots()\n",
    "    fig.subplots_adjust(bottom=0.15)\n",
    "\n",
    "    # display initial image \n",
    "    im_h = ax.imshow(data[idx, :, :], vmin=-1.0, vmax=1., cmap='bwr', interpolation='nearest')\n",
    "\n",
    "    # setup a slider axis and the Slider\n",
    "    ax_depth = plt.axes([0.23, 0.02, 0.56, 0.04])\n",
    "    slider_depth = Slider(ax_depth, 'depth', 0, data.shape[0]-1, valinit=idx)\n",
    "\n",
    "    # update the figure with a change on the slider \n",
    "    def update_depth(val):\n",
    "        idx = int(round(slider_depth.val))\n",
    "        im_h.set_data(data[idx, :, :])\n",
    "\n",
    "    slider_depth.on_changed(update_depth)\n",
    "\n",
    "    plt.show()\n",
    "    \n",
    "### Several Functions for Network Plotting ###\n",
    "def plot_network_progression_inv(history):\n",
    "    turns = history.current_turn\n",
    "    seed_turn = 0\n",
    "    adj_rounded = get_adj_rounded(history.invest_matrix[seed_turn])\n",
    "\n",
    "    G = nx.from_numpy_matrix(adj_rounded)\n",
    "    G = nx.convert_node_labels_to_integers(G, first_label=0, ordering='default', label_attribute=None)\n",
    "\n",
    "    pos = nx.spring_layout(G)\n",
    "    #pos = nx.spectral_layout(G)\n",
    "    #pos = nx.circular_layout(G)\n",
    "    #pos = nx.shell_layout(G)\n",
    "    #pos = nx.kamada_kawai_layout(G)\n",
    "    \n",
    "    step = 30\n",
    "    wait_time = .4\n",
    "    for t in range(0,turns,step):\n",
    "        plt.cla()\n",
    "        adjacency = get_adjacency(history.invest_matrix[t]) * 4\n",
    "        adjacency.dtype = [('weight','float')]\n",
    "        G = nx.from_numpy_matrix(adjacency,parallel_edges=False)\n",
    "        pos = nx.spring_layout(G, pos=pos, fixed=None)\n",
    "        color_map = []\n",
    "        for node in G:\n",
    "            if int(node) < history.n_males:\n",
    "                color_map.append('blue')\n",
    "            else: color_map.append('red')    \n",
    "        edges = G.edges()\n",
    "        weights = [G[u][v]['weight'] for u,v in edges]\n",
    "        nx.draw(G, pos, node_color = color_map, with_labels = True, edges=edges, width=weights, node_size=450)\n",
    "        plt.pause(wait_time)\n",
    "        plt.draw()\n",
    "\n",
    "def plot_network_progression(history):\n",
    "    turns = history.current_turn\n",
    "    seed_turn = 0\n",
    "    adj_rounded = np.round(history.adjacency_matrix[seed_turn] + .4) * 4\n",
    "    G = nx.from_numpy_matrix(adj_rounded)\n",
    "    #G = nx.convert_node_labels_to_integers(G, first_label=0, ordering='default', label_attribute=None)\n",
    "\n",
    "    pos = nx.spring_layout(G)\n",
    "    #pos = nx.spectral_layout(G)\n",
    "    #pos = nx.circular_layout(G)\n",
    "    #pos = nx.shell_layout(G)\n",
    "    #pos = nx.kamada_kawai_layout(G)\n",
    "    \n",
    "    #turns = 160\n",
    "    step = turns / 20\n",
    "    wait_time = .5\n",
    "    for t in range(0,turns,step):\n",
    "        plt.cla()\n",
    "        adj_rounded = history.adjacency_matrix[t]\n",
    "        adj_rounded[adj_rounded < .1] = 0.0\n",
    "        adj_rounded.dtype = [('weight','float')]\n",
    "        G = nx.from_numpy_matrix(adj_rounded,parallel_edges=False)\n",
    "        pos = nx.spring_layout(G, pos=pos, fixed=None)\n",
    "        color_map = []\n",
    "        for node in G:\n",
    "            if int(node) < history.n_males:\n",
    "                color_map.append('blue')\n",
    "            else: color_map.append('red')    \n",
    "        edges = G.edges()\n",
    "        weights = [G[u][v]['weight'] for u,v in edges]\n",
    "        nx.draw(G, pos, node_color = color_map, with_labels = True, edges=edges, width=weights, node_size=450)\n",
    "        plt.pause(wait_time)\n",
    "        plt.draw()\n",
    "        \n",
    "        \n",
    "def get_adj_rounded(investment):\n",
    "    n_males,n_females = np.shape(investment)\n",
    "    size = n_males + n_females\n",
    "    adjacency = np.zeros([size,size])\n",
    "    adjacency[0:n_males,n_males:] = investment\n",
    "    adjacency[n_males:,:n_females] = np.transpose(investment)\n",
    "    adj_rounded = np.round(adjacency + .4)\n",
    "    return adj_rounded\n",
    "\n",
    "def get_adjacency(investment):\n",
    "    n_males,n_females = np.shape(investment)\n",
    "    size = n_males + n_females\n",
    "    adjacency = np.zeros([size,size])\n",
    "    adjacency[0:n_males,n_males:] = investment\n",
    "    adjacency[n_males:,:n_females] = np.transpose(investment)\n",
    "    return adjacency\n",
    "    \n",
    "def plot_network(investment):\n",
    "    plt.cla()\n",
    "    adjacency = get_adjacency(investment) * 4\n",
    "    adjacency.dtype = [('weight','float')]\n",
    "    G = nx.from_numpy_matrix(adjacency,parallel_edges=False)\n",
    "    # Set layount: \n",
    "    pos = nx.spring_layout(G)\n",
    "    color_map = []\n",
    "    n_males, n_females = np.shape(investment)\n",
    "    for node in G:\n",
    "        if int(node) < n_males:\n",
    "            color_map.append('blue')\n",
    "        else: color_map.append('red')    \n",
    "\n",
    "    edges = G.edges()\n",
    "    weights = [G[u][v]['weight'] for u,v in edges]\n",
    "    nx.draw(G, pos, node_color = color_map, with_labels = True, edges=edges, width=weights)\n",
    "\n",
    "    plt.show()\n",
    "        \n",
    "# List of Menu options\n",
    "def print_menu():\n",
    "    print \"[0] - Run Default Simulation\"    \n",
    "    print \"[1] - Run Default Single Trial\"\n",
    "    print \"[2] - Run Custom Simulation (set parameters)\"\n",
    "    print \"[9] - Quit\"\n",
    "\n",
    "#Menu allowing you to set the parameters (based on sys.argv)\n",
    "def menu():\n",
    "    choice = 0\n",
    "    while choice != 9:\n",
    "        print_menu()\n",
    "        choice = int(raw_input(\"Select option (enter a number):  \"))\n",
    "        if choice == 0:\n",
    "            run_simulation()\n",
    "        elif choice == 1:\n",
    "            run_trial()\n",
    "        elif choice == 2:\n",
    "            build_simulation()\n",
    "        elif choice == 9:\n",
    "            print \"How about a nice game of chess?\"\n",
    "    return choice\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    #menu()\n",
    "    \n",
    "    #history = run_trial()\n",
    "    #run_simulation(alphas = [1.5]*N_FEMALES)\n",
    "    #plot_history(history)\n",
    "    #plot_network_progression(history)\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No strategy: 7 found, quitting\n"
     ]
    },
    {
     "ename": "SystemExit",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[0;31mSystemExit\u001b[0m\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "'level' is an invalid keyword argument for this function",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/Users/Ammon/anaconda/lib/python2.7/site-packages/IPython/core/interactiveshell.pyc\u001b[0m in \u001b[0;36mrun_code\u001b[0;34m(self, code_obj, result)\u001b[0m\n\u001b[1;32m   2875\u001b[0m                 \u001b[0mresult\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0merror_in_exec\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2876\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshowtraceback\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexception_only\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2877\u001b[0;31m             \u001b[0mwarn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"To exit: use 'exit', 'quit', or Ctrl-D.\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlevel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2878\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcustom_exceptions\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2879\u001b[0m             \u001b[0metype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtb\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexc_info\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: 'level' is an invalid keyword argument for this function"
     ]
    }
   ],
   "source": [
    "# Fairly Default trial: \n",
    "params = Parameters(alpha = 1.5, bump_constant=0)\n",
    "history = run_trial(params)\n",
    "#plot_history(history)\n",
    "plot_history2(history.invest_matrix2)\n",
    "plot_history2(history.reward_matrix2)\n",
    "plot_history2(history.adjacency_matrix)\n",
    "plot_network_progression(history)\n",
    "\n",
    "#record = run_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Trial where investment is everthing: \n",
    "my_params = Parameters()\n",
    "my_params.alpha = 0\n",
    "my_params.beta = 1.5\n",
    "my_params.decay_constant = 0.001\n",
    "history = run_trial(params = my_params)\n",
    "plot_history2(history.invest_matrix2)\n",
    "plot_history2(history.reward_matrix2)\n",
    "plot_history2(history.adjacency_matrix)\n",
    "plot_network_progression(history)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Check effects of male variation\n",
    "\n",
    "plt.cla()\n",
    "span_list = []\n",
    "monos_means = []\n",
    "monos_stds = []\n",
    "\n",
    "\n",
    "for d in range(5):\n",
    "    n_males = N_MALES\n",
    "    span = d * .2\n",
    "    step = span / float(n_males)\n",
    "    res_max = 1.0 + span / 2\n",
    "    res_min = 1.0 - span / 2\n",
    "    ress_m = []\n",
    "    for m in range(n_males):\n",
    "        res_m = res_min + step * m\n",
    "        ress_m.append(res_m)\n",
    "    record = run_simulation(ress_males=ress_m, alpha=2.0, n_turns=1000, trials=10) \n",
    "    monos = []\n",
    "    for r in record:\n",
    "        monos.append(r['monogamy'])\n",
    "    span_list.append(span)\n",
    "    monos_means.append(np.mean(monos))\n",
    "    monos_stds.append(np.std(monos))\n",
    "\n",
    "plt.errorbar(span_list,monos_means, monos_stds)\n",
    "xmin, xmax = (0,2)\n",
    "ymin = min(monos_means) - 2 * max(monos_stds) - 1\n",
    "ymax = max(monos_means) + 2 * max(monos_stds) + 1\n",
    "\n",
    "plt.axis([xmin, xmax, ymin, ymax])\n",
    "\n",
    "plt.show()        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def crowding(invest):\n",
    "    fsums = invest.sum(0) ## Take the sum along the vertical axis\n",
    "    norm_fsums = fsums - 1   #Normalize to 1 (assumes equal resources for females...)\n",
    "    norm_sum = norm_fsums.sum()     #sum of norms\n",
    "    n_norm_sum = norm_sum / len(fsums)  # Normalize again for number of females\n",
    "    return n_norm_sum\n",
    "\n",
    "def get_waste(history):\n",
    "    waste = []\n",
    "    n_females = history.n_females\n",
    "    for r in history.reward_matrix:\n",
    "        lost_reward = n_females - sum(sum(r))\n",
    "        waste.append(lost_reward)\n",
    "    #plt.plot(waste)\n",
    "    #plt.show()\n",
    "    return waste\n",
    "\n",
    "def final_waste(history):\n",
    "    n_females = history.n_females\n",
    "    lost_reward = n_females - sum(sum(history.reward_matrix[-1]))\n",
    "    return lost_reward\n",
    "\n",
    "def get_all_degrees(history):\n",
    "    all_degrees = []\n",
    "    for i in history.invest_matrix():\n",
    "        adjacency = get_adjacency(history.invest_matrix[-1])\n",
    "        adjacency.dtype = [('weight','float')]\n",
    "        G = nx.from_numpy_matrix(adjacency,parallel_edges=False)\n",
    "        all_degrees.append(get_degrees(G))\n",
    "    return all_degrees\n",
    "\n",
    "def get_degrees_final(history):\n",
    "    adjacency = get_adjacency(history.invest_matrix[-1])\n",
    "    adjacency.dtype = [('weight','float')]\n",
    "    G = nx.from_numpy_matrix(adjacency,parallel_edges=False)\n",
    "    final_degrees = get_degrees(G)\n",
    "    return final_degrees\n",
    "    \n",
    "def get_deviation_final(history):\n",
    "    final_deviation = []\n",
    "    adjacency = get_adjacency(history.invest_matrix[-1])\n",
    "    adjacency.dtype = [('weight','float')]\n",
    "    G = nx.from_numpy_matrix(adjacency,parallel_edges=False)\n",
    "    return get_deviation(G)\n",
    "\n",
    "def get_degrees(G):\n",
    "    degrees = []\n",
    "    for n in G:\n",
    "        degrees.append(G.degree(n)\n",
    "    return degrees\n",
    "\n",
    "def get_deviation(G):\n",
    "    degrees = get_degrees(G)\n",
    "    deviation = 0\n",
    "    for d in degrees:\n",
    "        deviation += abs(d-1)\n",
    "    return deviation\n",
    "\n",
    "def get_deviations(history):\n",
    "    deviations = []\n",
    "    for i in history.invest_matrix():\n",
    "        adjacency = get_adjacency(i)\n",
    "        adjacency.dtype = [('weight','float')]\n",
    "        G = nx.from_numpy_matrix(adjacency,parallel_edges=False)\n",
    "        deviations.append(deviation(G))\n",
    "    return deviations\n",
    "\n",
    "def final_fornication(history):\n",
    "    final_fornication = 0\n",
    "    adjacency = get_adjacency(history.invest_matrix[-1])\n",
    "    adjacency.dtype = [('weight','float')]\n",
    "    G = nx.from_numpy_matrix(adjacency,parallel_edges=False)\n",
    "    final_fornication = fornication.append(deviation(G))\n",
    "    return final_fornication\n",
    "    \n",
    "def get_adjacency(investment):\n",
    "    n_males,n_females = np.shape(investment)\n",
    "    size = n_males + n_females\n",
    "    adjacency = np.zeros([size,size])\n",
    "    adjacency[0:n_males,n_males:] = investment\n",
    "    adjacency[n_males:,:n_females] = np.transpose(investment)\n",
    "    return adjacency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/site-packages/numpy/core/_methods.py:59: RuntimeWarning:\n",
      "\n",
      "Mean of empty slice.\n",
      "\n",
      "/usr/local/lib/python2.7/site-packages/numpy/core/_methods.py:82: RuntimeWarning:\n",
      "\n",
      "Degrees of freedom <= 0 for slice\n",
      "\n"
     ]
    }
   ],
   "source": [
    "## Check the effect of alpha\n",
    "\n",
    "plt.cla()\n",
    "alpha_list = []\n",
    "dev_means = []\n",
    "dev_stds = []\n",
    "mono = []\n",
    "\n",
    "for a in range(0,5):\n",
    "    al = 1.0 + a * .2\n",
    "    record = run_simulation(trials=15, n_turns=800, alpha=al)\n",
    "    devs = []\n",
    "    alpha_list.append(al)\n",
    "    for r in record:\n",
    "        mono.append(r['deviation'])\n",
    "\n",
    "    dev_means.append(np.mean(devs))\n",
    "    dev_stds.append(np.std(devs))\n",
    "\n",
    "plt.errorbar(alpha_list,dev_means, dev_stds)\n",
    "xmin, xmax = (0,2)\n",
    "ymin = min(dev_means) - 2 * max(dev_stds)\n",
    "ymax = max(dev_means) + 2 * max(dev_stds)\n",
    "\n",
    "plt.axis([xmin, xmax, ymin, ymax])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.cla()\n",
    "plt.errorbar(span_list,dev_means, dev_stds)\n",
    "xmin, xmax = (0,2)\n",
    "ymin = min(dev_means) - 2 * max(dev_stds) - 1\n",
    "ymax = max(dev_means) + 2 * max(dev_stds) + 1\n",
    "\n",
    "plt.axis([xmin, xmax, ymin, ymax])\n",
    "\n",
    "plt.show() \n",
    "\n",
    "history.adjacency_matrix"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}