%pylab inline
%load_ext autoreload
%autoreload 2

import SimulationInstance
import SimulationMain
import SimulationInstanceVectorized

SternJudging = [[1, 0 ],
                [0, 1]]
SimpleStanding = [[1, 1],
                  [0, 1]]
Shunning = [[1, 0],
            [0, 1]]
ImageScore = [[1, 1],
              [0, 0]]

def SSP(population_size, socialnorm):
    runs = 2
    generations = 3*np.power(10, 5)
    mutation_rate = np.power(10*population_size, -1)

    execution_error = 0.08
    reputation_assignment_error = 0.01
    private_assessment_error = 0.01
    reputation_update_probability = 0.2
    randomseed = np.random.randint(999999)
    cost = 1
    benefit = 5
    coop_index = SimulationInstanceVectorized.run_instance(runs, generations, population_size,
        mutation_rate, execution_error, reputation_assignment_error, private_assessment_error,
        reputation_update_probability, randomseed, socialnorm,
        cost, benefit)
    return coop_index

"""
Explanation: Santos, Santos, Pacheco Results
The results achieved in the simulations by SSP (Santos Santos Pachecho) in their paper Social Norms of Cooperation in Small-Scale Societies can be seen in the following graph.
<img src="http://journals.plos.org/ploscompbiol/article/figure/image?size=large&id=info:doi/10.1371/journal.pcbi.1004709.g001"/>
One can see an almost constant cooperation index for population sizes greater than 50 for Stern Judging, Shunning, and Image Score. While Simple Standing increases steadily over the range 50-300.
Below we will be performing the same simulation and plotting the results in order to verify and compare the simulations results and those of the paper.
First however let us define the social norms we will use as well as a function to pass varying population sizes and social norms for retrieving the results.
End of explanation
"""


population_sizes = [5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120]

stern_judging_results = [0.508516473,0.775682093,0.826424982,0.760320876,0.697826156,0.643724905,0.561129066,0.594354744,0.652819772,0.598121451,0.527160869,0.531084403,0.608492096]

stern_judging_results = np.array(stern_judging_results)
# results_temp_x2
#x = stern_judging_results[:,0]
#y = stern_judging_results[:,1]
plt.plot(population_sizes, stern_judging_results, 'k', population_sizes, stern_judging_results, 'r^', linewidth=1.5)
plt.ylim((-0.01, 1.1))
plt.xlim((-0.01, 125))
plt.xticks([10*i for i in range(13)])
plt.yticks([0.1*i for i in range(11)])
plt.ylabel('Cooperation Index')
plt.xlabel('Population Size')
plt.grid(True)
plt.show()

"""
Explanation: Stern Judging
Now let us retrieve the results for populations in the ranges 10-150, with simulations at every multiple of 10, and then simulations at multiples of 25 from 150-300. 
Let's begin with the social norm of Stern Judging.
End of explanation
"""


coop_index = [0.508516473,0.775682093,0.826424982,0.760320876,0.697826156,0.643724905,0.561129066,0.594354744,0.652819772,0.598121451,0.527160869,0.531084403,0.608492096]
population_sizes_temp = [5,10,20,30,40,50,60,70,80,90,100,110,120]
zero_proportion = [0.472993542,0.135284757,0.039967674,0.127199045,0.183304844,0.212918715,0.317619392,0.277463455,0.217788767,0.277047708,0.345624097,0.345643889,0.257104097]
coop_index_wo_zeros = [0.964804944,0.896984755,0.860765941,0.869928644,0.853073938,0.814867559,0.814542487,0.786577321,0.770007229,0.814391772,0.795729616,0.797411836,0.811981669]
line_standard, = plt.plot(population_sizes_temp, coop_index, linewidth=1.5, label='Coop Index')
line_zero_prop, = plt.plot(population_sizes_temp, zero_proportion, label='Zero Proportion', marker='^')
line_wo_zeros, = plt.plot(population_sizes_temp, coop_index_wo_zeros, label='Coop Index w/o zeros', marker='o')
plt.ylim((-0.01, 1.05))
plt.xlim((-0.1, 125))
plt.legend(handles=[line_standard, line_wo_zeros, line_zero_prop], loc='upper left', bbox_to_anchor=(1,1))
plt.grid(b=True, which='both')
plt.rcParams["figure.figsize"][0] = 9
plt.rcParams["figure.figsize"][1] = 5
plt.xticks([10*i for i in range(13)])
plt.yticks([0.1*i for i in range(11)])
plt.ylabel('Cooperation Index/Proportion')
plt.xlabel('Population Size')
plt.show()

"""
Explanation: The results of the simulation as shown above do not seem to correlate with those of the paper by santos santos pacheco. However what seems to be the case is that the results of a number of runs seem to be either zero or a value close to 0.8. This means that some proportion of the runs seem to be fixated by defection, unable to recover to a stable cooperative population. While we cannot simply remove the zero values, if we take into account both the cooperation index average with zeros taken into account and the average without, we see the following.
We also look at the proportion of runs which have result in a cooperation index of zero against the total number of runs:
End of explanation
"""


coop_index = [0.55671953,0.717854758,0.764312844,0.748989805,0.707472638,0.834710665,0.838990698,0.767615598,0.708101936,0.709383435,0.671574256,0.730423794,0.606060629]
population_sizes_temp = [5,10,20,30,40,50,60,70,80,90,100,110,120]
zero_proportion = [0.286346651,0.172696458,0.145584712,0.158310288,0.186614808,0.074261587,0.076428718,0.129578862,0.171840962,0.171456635,0.185223429,0.147639872,0.235549487]
coop_index_wo_zeros = [0.780095131,0.86768886,0.89447719,0.889273361,0.869099708,0.900894275,0.907260109,0.878697814,0.849517584,0.854871055,0.820141081,0.852480787,0.781452704]
line_standard, = plt.plot(population_sizes_temp, coop_index, linewidth=1.5, label='Coop Index')
line_zero_prop, = plt.plot(population_sizes_temp, zero_proportion, label='Zero Proportion', marker='^')
line_wo_zeros, = plt.plot(population_sizes_temp, coop_index_wo_zeros, label='Coop Index w/o zeros', marker='o')
plt.ylim((-0.01, 1.05))
plt.xlim((-0.1, 125))
plt.legend(handles=[line_standard, line_wo_zeros, line_zero_prop], loc='upper left', bbox_to_anchor=(1,1))
plt.grid(b=True, which='both')
plt.rcParams["figure.figsize"][0] = 9
plt.rcParams["figure.figsize"][1] = 5
plt.xticks([10*i for i in range(13)])
plt.yticks([0.1*i for i in range(11)])
plt.ylabel('Cooperation Index/Proportion')
plt.xlabel('Population Size')
plt.show()

"""
Explanation: While initially one may come to tje conclusion that the zero proportion, which takes the form above of 1 minus the cooperation index, to be merely the defection index, however the way this is calculated is not by looking at the defection index, but moreso the number of runs which have a cooperation index of zero (calculated as less than 10^-5). This leads one to believe that the distribution of the cooperation index is zero inflated. WHere we see that for all cooperation indexes, it is either a value distributed around ~0.8 or it is equal to zero with certain probability.
Stern judging with high mutation and low execution error
The following data comes from a series of runs where the mutation rate is set to Z^-1 where Z is the population size, and where the execution error (epsilon) is 0.01 compared to the 0.08 value used in the previous simulations.
End of explanation
"""


coop_index = [0.270753407,0.33029135,0.365173906,0.255667253,0.215493488,0.15236222,0.150375082,0.123795391,0.145404647,0.098041625,0.079193873,0.064620633,0.09713768]
population_sizes_temp = [5,10,20,30,40,50,60,70,80,90,100,110,120]
zero_proportion = [0.68759684,0.601973681,0.549900174,0.677613316,0.717450052,0.779782604,0.783594149,0.813146215,0.788442622,0.842165365,0.860564097,0.87703309,0.839349549]
coop_index_wo_zeros = [0.866315918,0.827693345,0.806530326,0.736667922,0.66959725,0.50391085,0.386810394,0.337837888,0.419821748,0.354295807,0.292116211,0.234500168,0.241363372]
line_standard, = plt.plot(population_sizes_temp, coop_index, linewidth=1.5, label='Coop Index')
line_zero_prop, = plt.plot(population_sizes_temp, zero_proportion, label='Zero Proportion', marker='^')
line_wo_zeros, = plt.plot(population_sizes_temp, coop_index_wo_zeros, label='Coop Index w/o zeros', marker='o')
plt.ylim((-0.01, 1.05))
plt.xlim((-0.1, 125))
plt.legend(handles=[line_standard, line_wo_zeros, line_zero_prop], loc='upper left', bbox_to_anchor=(1,1))
plt.grid(b=True, which='both')
plt.rcParams["figure.figsize"][0] = 9
plt.rcParams["figure.figsize"][1] = 5
plt.xticks([10*i for i in range(13)])
plt.yticks([0.1*i for i in range(11)])
plt.ylabel('Cooperation Index/Proportion')
plt.xlabel('Population Size')
plt.show()

"""
Explanation: Simple Standing
End of explanation
"""


population_sizes = [5, 12, 25, 35, 50, 70]

# Simple Standing
coop_index_SS_5 = [0.30300328890107,0.3122170267297286, 0.3694477572139156, 0.32269225071203966,0.34660439483109196,0.35753302790010333,0.34472536935321985,0.33506662987640695,]
coop_index_SS_12 = [    0.389835217625291,    0.503741053593408,    0.414803410665924,    0.373144341121016,    0.404481667833337,    0.409381477004441,    0.411456629366117,    0.417190360092897]
coop_index_SS_25 = [    0.575895970725817,    0.593555191685895,    0.346345244241645,    0.586595733317530,    0.435658189104404,    0.510932066415060,    0.522842995714799,    0.464933493087011]
coop_index_SS_35 = [    0.6024419644486968,    0.5611443005344464,    0.6177017545265285,    0.5946591774980875,    0.5203268341659013,    0.5146393535395558,    0.5037492489579137,    0.5102253043114068]
coop_index_SS_50 = [    0.625971801377042,    0.728209162811714,    0.526364328756773,    0.663999642961640,    0.428279452701676,    0.557728084918833,    0.696717869325442,    0.541151875932873]
coop_index_SS_70 = [    0.6142549683828524]
coop_index_SS = [np.average(coop_index_SS_5), np.average(coop_index_SS_12), np.average(coop_index_SS_25), np.average(coop_index_SS_35), np.average(coop_index_SS_50), np.average(coop_index_SS_70)]

# Stern Judging
coop_index_SJ_5 = [0.6415401538524217,0.6446142578305846,0.6410012259304664,0.6680960299852963,0.6918668402475565,0.6640123852310564,0.6769155811688191,0.6492853677235841,]
coop_index_SJ_12 = [0.7954687062112575,0.8120924624823713,0.8133768743896582,0.8168074297908923,0.8086529789359302,0.8090324331448141,0.8007545571928044,0.8254361654528761,]
coop_index_SJ_25 = [0.828478139548665, 0.8283964836906803,0.8284225578613255,0.8284125859725406,0.8284379690295149,0.8284406979285108,0.8284578476216579,0.82846780617706,]
coop_index_SJ_35 = [0.8286760075679797,0.8286752539748635,0.8286820007874557,0.8286810253492246,0.8286692010423924,0.8286773564013956,0.8286702717190714,0.8286818004843707,]
coop_index_SJ_50 = [0.8289147343138308,0.8289050916855247,0.8289118700439632,0.8289093383679199,0.8289100281852896,0.8289160067259508,0.8289107137433842,0.8289103609642998,]
coop_index_SJ_70 = [0.8340885668850333,0.834185941885406, 0.8340846263499717,0.8341202176521099,0.8341425765400396,0.8341473538307335,0.8340992111946139]
coop_index_SJ = [np.average(coop_index_SJ_5), np.average(coop_index_SJ_12), np.average(coop_index_SJ_25), np.average(coop_index_SJ_35), np.average(coop_index_SJ_50), np.average(coop_index_SJ_70)]

# Zero Norm
coop_index_ZERO_5 = [0.019335418939590454,0.02019631862640381,0.02082008123397827,0.02134975790977478,0.022004514932632446,0.02099481225013733,0.020933866500854492,0.020770877599716187,]
coop_index_ZERO_12 = [0.010089797529085499,0.010102098243367979,0.00979763145213295,0.010382693418411838,0.009933781647750958,0.010198327380636389,0.010313705251377474,0.01053330804517383,]
coop_index_ZERO_25 = [0.008694934425683058,0.008677666954322052,0.00968594105904096, 0.008913428568253957,0.008520019703776681,0.008646592024414906,0.007169612480046222,0.008183021418597846,]
coop_index_ZERO_35 = [0.00755965195940635, 0.008507915108680297,0.00806827062322453, 0.007357088850895971,0.007357534808707909,0.006010654914746094,0.007110214164814367,0.005497417186201219,]
coop_index_ZERO_50 = [0.006506138356067986,0.004333089752184267,0.005001832784565869,0.005762029303812237,0.005873097678139549,0.005450347408694107,0.005294959610961522,0.005296161652450565,]
coop_index_ZERO_70 = [0.0009228635878207811,0.0006507234238636115,0.0006376077663370713,0.0004388977710386438,0.0007317914197910949,0.00239054263795588,0.0007066109140949628,0.0007344184431828324,]
coop_index_ZERO = [np.average(coop_index_ZERO_5), np.average(coop_index_ZERO_12), np.average(coop_index_ZERO_25), np.average(coop_index_ZERO_35), np.average(coop_index_ZERO_50), np.average(coop_index_ZERO_70)]

# Image Scoring
coop_index_IS_5 =[0.13157886266708374,0.12797850370407104,0.12999996542930603,0.1389734148979187, 0.1363530158996582, 0.12969133257865906,0.139630526304245, 0.1444832682609558, ]
coop_index_IS_12 =[0.07933607058626967,0.06608765153348305,0.07403314135216962,0.07385864546411901,0.07386893035303423,0.0728666632664483, 0.07001568498607608,0.07498191672829241,]
coop_index_IS_25 =[0.051464263200974816,0.05746115207834103,0.05665677206211875,0.055269822480824, 0.049308401882969935,0.0571610102259638, 0.060479542516824136,0.05478295697391472,]
coop_index_IS_35 =[0.04855178866497691,0.04499552001364177,0.04594197854549877,0.045567311247033344,0.04700417219988192,0.04882126547558642,0.04925220849268452,0.05192522821891026,]
coop_index_IS_50 =[0.049498805921503415,0.045692120526372317,0.04503519772295453,0.04016098504068533,0.03694476127689362,0.04080260482446231,0.04433218178585847,0.04906098281845488,]
coop_index_IS_70 =[0.0728822786061512,0.08704757792335796,0.06320217416341833,0.000601568897288247,]
coop_index_IS = [np.average(coop_index_IS_5), np.average(coop_index_IS_12), np.average(coop_index_IS_25), np.average(coop_index_IS_35), np.average(coop_index_IS_50), np.average(coop_index_IS_70)]


line_SS, = plt.plot(population_sizes, coop_index_SS, label='Simple Standing', marker='o', ms=8, linewidth=2)
line_SJ, = plt.plot(population_sizes, coop_index_SJ, label='Stern Judging', marker='o', ms=8, linewidth=2)
line_ZERO, = plt.plot(population_sizes, coop_index_ZERO, label='Uniform Zero', marker='o', ms=8, linewidth=2)
line_IS, = plt.plot(population_sizes, coop_index_IS, 'b--', label='Image Scoring', marker='s', ms=8, linewidth=2)
plt.ylim((-0.01, 1.05))
plt.xlim((-0.1, 125))
plt.legend(handles=[line_SS, line_SJ, line_ZERO, line_IS], loc='upper left', bbox_to_anchor=(1,1))
plt.grid(b=True, which='both')
plt.rcParams["figure.figsize"][0] = 9
plt.rcParams["figure.figsize"][1] = 5
plt.xticks([10*i for i in range(13)])
plt.yticks([0.1*i for i in range(11)])
plt.ylabel('Cooperation Index')
plt.xlabel('Population size')
plt.show()

"""
Explanation: Complete rerun of results with 1...Z loop in each generation
For a simulation with parameters:
R = 8, Z = 20, G = 3x10^3, 48 seconds
R = 8, Z = 20, G = 3x10^4, 432 seconds
R = 8, Z = 40, G = 3x10^3, 92 seconds
Estimated
Z = 20
Estimated R = 8, G = 3x10^5, ~4320 seconds (1.2 hours)
Estimated R = 32 (96//3), G = 3x10^5, ~5 hours
Z = 40
Estimated R = 8, G = 3x10^5, ~9200 seconds (2.5 hours)
Estimated R = 32 (96//3), G = 3x10^5, ~10 hours
Z = 80
Estimated R = 32 (96//3), G = 3x10^5, ~20 hours
Z = 120
Estimated R = 32 (96//3), G = 3x10^5, ~30 hours
Results
Social Norm Simple Standing
End of explanation
"""


