#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Last update, 16/01/2018
Continuous Replicator Dynamic in Asexual model
@author: Shota Shibasaki
"""

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

#parameters
a1=1
b1=1.5
e1=0.0001
a2=0.5
b2=1.2
e2=0.0001
h=0.5#effect of forming smaller fruiting bodies
theta=0.5#prob of macrocyst formation
e0=0.000001

def Matrix(a2,b2,theta):
   """
   Defining patoff matrixes
   A1: payoff matrix for fruiting body formation
   A2: payoff matrix for macrocyst formation
   A: patoff matrix generated by the mixture of A1 and A2.
      see Eq [2] in the main text for more detail.
   Parameters
   a2: benefit of cmutual ooperation in macrocyst formation
   b2: benefit of exploitation for defectors in macrocyst formation
   theta: prob. of macrocyst fromation
      
   """
   A1=np.array([[a1,h*a1,0],\
                 [e0,e0,e0],\
                 [b1,h*e1,e1]])
    
   A2=np.array([[a2,0,a2],\
                 [b2,e2,b2],\
                 [a2,0,a2]])
   A=(1-theta)*A1+(theta)*A2
   return A

def func(x, t, a2,b2,theta):
    """
    Continuous replicator dynamics
    x: array of the frequency of the three strategy
    t: time
    a2, b2, theta: defined in function Matrix
    """
    #x is array and A is matrix
    A=Matrix(a2,b2,theta)
    phi=np.dot(x,A)
    phi=np.dot(phi,x.T)#mean fitness in the population
    #calculate replicator dynamics
    f0=np.dot(A[0,:],x.T)-phi
    f1=np.dot(A[1,:],x.T)-phi
    f2=np.dot(A[2,:],x.T)-phi     
    f=np.array([x[0]*f0,x[1]*f1,x[2]*f2])
    return(f)

def triplot(sol,col):
    """
    Mapping the phase space into the triangl plain
    """
    TP=np.array([1/2,np.sqrt(3)/2])#C
    LS=np.array([0,0])#D_M
    RS=np.array([1,0])#D_F
    orbX=sol[:,0]*TP[0]+sol[:,1]*LS[0]+sol[:,2]*RS[0]
    orbY=sol[:,0]*TP[1]+sol[:,1]*LS[1]+sol[:,2]*RS[1]
    #plot the dynamics
    plt.xlim(-0.2,1.2)
    plt.ylim(-0.1,1.2*TP[1])
    plt.plot([LS[0],TP[0]],[LS[1],TP[1]],linewidth=4,color="k")
    plt.plot([RS[0],TP[0]],[RS[1],TP[1]],linewidth=4,color="k")
    plt.plot([RS[0],LS[0]],[RS[1],LS[1]],linewidth=4,color="k")
    plt.text(-0.15,0,'$D^M$',fontsize=20,)
    plt.text(1.05,0,'$D^F$',fontsize=20,)
    plt.text(0.475,0.9,'$C$',fontsize=20,)
    plt.plot(orbX,orbY,linewidth=2,color=col)
    plt.tick_params(labelbottom='off',labelleft="off")

def main():
    """
    Simulate the continuous replicator dymnamics in asexual model
    from the three different initial conditions
    """
    T=1000# end the simulation at T
    init=np.array([10,1,1])#initial condition 1
    init=init/sum(init)#normalize
    t=np.linspace(0,T,100*T+1)#step size
    #using odeint function for integration
    sol=odeint(func, init, t, args=(a2,b2,theta))
    triplot(sol,"c")
    
    init=np.array([1,10,1])#initial condition 2
    init=init/sum(init)
    t=np.linspace(0,T,100*T+1)
    sol=odeint(func, init, t, args=(a2,b2,theta))
    triplot(sol,"m")
    
    init=np.array([1,1,12]) #initial condition 3
    init=init/sum(init)

    t=np.linspace(0,T,100*T+1)
    sol=odeint(func, init, t, args=(a2,b2,theta))
    triplot(sol,"y")
    
    """
    Check the inteirior equilibrium G
    """
    GX=sol[-1,0]*0.5+sol[-1,1]*0+sol[-1,2]*1
    GY=sol[-1,0]*np.sqrt(3)/2
    plt.plot(GX,GY,'k*',markersize=15)
    plt.text(GX+0.05,GY,'G',fontsize=20)
    print(GX,GY) # coordinate of interior equilibrium G
    plt.savefig("Cont_rep_asexual.pdf")
    plt.show()
    
    
if __name__ == '__main__':
    main()