From 1881c36e1f8a940dfd24e6a82f1279098c6c1025 Mon Sep 17 00:00:00 2001 From: Colin Campbell <44878672+ColinECampbell@users.noreply.github.com> Date: Sat, 21 Aug 2021 19:10:11 -0400 Subject: [PATCH] Add files via upload Initial commit --- analysis.py | 927 +++++++++++++++++++++++++++++++++++++++++ create_graphs.py | 35 ++ find_attractors.py | 115 ++++++ function_library.py | 985 ++++++++++++++++++++++++++++++++++++++++++++ invasions_random.py | 251 +++++++++++ invasions_whole.py | 234 +++++++++++ 6 files changed, 2547 insertions(+) create mode 100644 analysis.py create mode 100644 create_graphs.py create mode 100644 find_attractors.py create mode 100644 function_library.py create mode 100644 invasions_random.py create mode 100644 invasions_whole.py diff --git a/analysis.py b/analysis.py new file mode 100644 index 0000000..b5d7cea --- /dev/null +++ b/analysis.py @@ -0,0 +1,927 @@ +""" +Here we analyze the output of invasions_whole.py and invasions_random.py: +we have attractors that correspond to starting from a combination of states. + + +Master workflow: + 1. What are the species interaction graphs? (create_graphs.py) + 2. What are the stable communities in these graphs, if we segregate them + in to two "half communities" to enforce isolated communities that + nonetheless have inter-community interactions for when invasions occur? + (find_attractors.py) + 3. What happens if two stable communities are combined and allowed to + restabilize? (invasions_whole.py and invasions_random.py) + 4. What do the results from (3) look like? (THIS FILE) + +The file is structured so the pre-processing near the top (to about line 360) +is saved to file. The top can then be commented out and the final analysis +can be run to load those results, generate figures, run stats, etc. + +Some pre-publication figures are included for benefit of anyone who gets +this far and would like to start playing around with their own analysis. + +Python 3.x + +Colin Campbell +campbeco@mountunion.edu +""" + +import pickle +import glob +import numpy as np +import function_library as nv +import matplotlib.pyplot as plt +import time +from scipy.stats import spearmanr, ks_2samp + +TIME0 = time.time() + + +# Run analysis =============================================================== + +# # Load in graphs +# #Local +# f = open(r'graphlist.data','rb') + +# graphs = pickle.load(f) +# f.close() + +# # Stable community invasions -------------------- +# s_overlap = [] # compare source state (stable community + invaders) to stable state +# s_start_size = [] # how many species in the start state? +# s_end_size = [] # how many species in the end attractor (SS or LC superset)? +# s_delta_size = [] +# s_nestedness = [] +# s_connectance = [] +# s_invasive_species = [] +# s_invasive_species_end = [] +# s_nestedness_final = [] +# s_connectance_final = [] +# s_native_survive = [] +# s_dbl_positive = [] +# s_exception_count = 0 + + +# for fname in glob.glob(r'attractors_comb_whole_[0-9]*.data'): +# # Open the file for this batch of interaction networks +# f = open(fname,'rb') +# data = pickle.load(f) +# f.close() + +# PBS_ITERATOR = int(fname[fname.rfind('_')+1:fname.rfind('.')]) + +# TIME1 = time.time() +# print('Whole invasion, starting file {0} at {1:.2f} minutes elapsed.'.format(PBS_ITERATOR,(TIME1-TIME0)/60)) + +# for graph_nmbr in data.keys(): +# graph_number = graph_nmbr - 10*PBS_ITERATOR + +# # Load up this graph so we can look at species-species interactions below +# graph = graphs[graph_number] + +# for counter in data[graph_nmbr].keys(): +# # pull out parent states and final state +# start_AB = data[graph_nmbr][counter][0]['start_state'][0] +# start_BA = data[graph_nmbr][counter][1]['start_state'][0] + +# end_AB = data[graph_nmbr][counter][0]['stable_state'] +# end_BA = data[graph_nmbr][counter][1]['stable_state'] + +# #parents could be SS or (more likely) LC, so form superset +# if type(end_AB[0]) == list: +# end_AB = nv.superset(end_AB[0]) +# # end_AB = end_AB[0][0] # alternative: choose just 1 state to consider +# else: +# end_AB = end_AB[0] +# if type(end_BA[0]) == list: +# end_BA = nv.superset(end_BA[0]) +# # end_BA = end_BA[0][0] +# else: +# end_BA = end_BA[0] + +# # track fraction of edges that are double positive (mutually beneficial) +# edge_count_AB = 0 +# edge_count_BA = 0 +# dbl_positive_count_AB = 0 +# dbl_positive_count_BA = 0 +# for edge in graph.edges(): + +# # need to translate from e.g. 'po_1' to position in the state list +# # indices 0-49 are plants 0-49 and indices 50-199 are pollinators 0-149 + +# # first node +# if edge[0][:2] == 'pl': edge0_index = int(edge[0][3:]) +# else: edge0_index = int(edge[0][3:])+50 + +# # second node +# if edge[1][:2] == 'pl': edge1_index = int(edge[1][3:]) +# else: edge1_index = int(edge[1][3:])+50 + +# # we only care about edges between present species +# if start_AB[edge0_index] == '1' and start_AB[edge1_index] == '1': +# edge_count_AB += 1 +# # and track seperately the edges that are also double positive +# if graph[edge[0]][edge[1]]['data'] == 1.0 and graph[edge[1]][edge[0]]['data'] == 1.0: dbl_positive_count_AB += 1 + +# # repeat for BA invasion +# if start_BA[edge0_index] == '1' and start_BA[edge1_index] == '1': +# edge_count_BA += 1 +# if graph[edge[0]][edge[1]]['data'] == 1.0 and graph[edge[1]][edge[0]]['data'] == 1.0: dbl_positive_count_BA += 1 + +# # once we've gone through all the edges, record the fraction of double positive edges +# # Because start_state can be a single state from a LC, occasionally we pick up +# # a state with no edges whatsoever and the ratio here gives a divide-by-zero error. Record as 0 but +# # keep track of how many +# if edge_count_AB == 0: +# s_exception_count += 1 +# s_dbl_positive.append(0) +# else: +# s_dbl_positive.append(dbl_positive_count_AB/edge_count_AB) +# if edge_count_BA == 0: +# s_exception_count += 1 +# s_dbl_positive.append(0) +# else: s_dbl_positive.append(dbl_positive_count_BA/edge_count_BA) + + +# # track # of invader species +# # start +# A_invaders = start_AB[:25]+start_AB[50:125] +# B_invaders = start_BA[25:50]+start_BA[125:200] +# s_invasive_species.append(A_invaders.count('1')) +# s_invasive_species.append(B_invaders.count('1')) + +# # end +# A_invaders_end = end_AB[:25]+end_AB[50:125] +# B_invaders_end = end_BA[25:50]+end_BA[125:200] +# s_invasive_species_end.append(A_invaders_end.count('1')) +# s_invasive_species_end.append(B_invaders_end.count('1')) + +# # compare the starting state (source + invaders) to final stable state +# s_overlap.append(nv.overlap(end_AB,start_AB)) +# s_overlap.append(nv.overlap(end_BA,start_BA)) + +# s_start_size.append(data[graph_nmbr][counter][0]['start_state'][0].count('1')) +# s_start_size.append(data[graph_nmbr][counter][1]['start_state'][0].count('1')) + +# s_end_size.append(end_AB.count('1')) +# s_end_size.append(end_BA.count('1')) + +# s_delta_size.append(100*(s_end_size[-2] - start_AB.count('1'))/start_AB.count('1')) +# s_delta_size.append(100*(s_end_size[-1] - start_BA.count('1'))/start_BA.count('1')) + +# s_nestedness.append(data[graph_nmbr][counter][0]['SN'][0][0]) +# s_nestedness.append(data[graph_nmbr][counter][1]['SN'][0][0]) + +# s_connectance.append(data[graph_nmbr][counter][0]['C'][0][0]) +# s_connectance.append(data[graph_nmbr][counter][1]['C'][0][0]) + +# s_nestedness_final.append(data[graph_nmbr][counter][0]['SN'][0][-1]) +# s_nestedness_final.append(data[graph_nmbr][counter][1]['SN'][0][-1]) + +# s_connectance_final.append(data[graph_nmbr][counter][0]['C'][0][-1]) +# s_connectance_final.append(data[graph_nmbr][counter][1]['C'][0][-1]) + +# # track % of native species that make it to final state +# A_native = start_AB[25:50]+start_AB[125:200] +# A_end = end_AB[25:50]+end_AB[125:200] + +# B_native = start_BA[:25] + start_BA[:125] +# B_end = end_BA[:25]+end_BA[:125] + +# A_native_cnt = A_native.count('1') +# A_end_cnt = 0 +# for i in range(len(A_native)): +# if A_native[i] == '1' and A_end[i] == '1': A_end_cnt += 1 +# s_native_survive.append(A_end_cnt/A_native_cnt) + +# B_native_cnt = B_native.count('1') +# B_end_cnt = 0 +# for i in range(len(B_native)): +# if B_native[i] == '1' and B_end[i] == '1': B_end_cnt += 1 +# s_native_survive.append(B_end_cnt/B_native_cnt) + + +# print('s exception fraction = {0:.5g}'.format(s_exception_count/(s_exception_count + len(s_dbl_positive)))) + +# TIME1 = time.time() +# print('Halfway done in {0:.2f} minutes elapsed.'.format((TIME1-TIME0)/60)) + +# # Random community invasions -------------------- +# r_overlap = [] # compare source state (stable community + invaders) to stable state +# r_start_size = [] # how many species in the start state? +# r_end_size = [] # how many species in the end attractor (SS or LC superset)? +# r_delta_size = [] +# r_nestedness = [] +# r_connectance = [] +# r_invasive_species = [] +# r_invasive_species_end = [] +# r_nestedness_final = [] +# r_connectance_final = [] +# r_native_survive = [] +# r_dbl_positive = [] +# r_exception_count = 0 + +# for fname in glob.glob(r'attractors_comb_random_[0-9]*.data'): +# # Open the file for this batch of interaction networks +# f = open(fname,'rb') +# data = pickle.load(f) +# f.close() + +# PBS_ITERATOR = int(fname[fname.rfind('_')+1:fname.rfind('.')]) + +# TIME1 = time.time() +# print('Random invasion, starting file {0} at {1:.2f} minutes elapsed.'.format(PBS_ITERATOR,(TIME1-TIME0)/60)) + +# for graph_nmbr in data.keys(): +# graph_number = graph_nmbr - 10*PBS_ITERATOR + +# # Load up this graph so we can look at species-species interactions below +# graph = graphs[graph_number] + +# for counter in data[graph_nmbr].keys(): +# # pull out parent states and final state +# start_AB = data[graph_nmbr][counter][0]['start_state'][0] +# start_BA = data[graph_nmbr][counter][1]['start_state'][0] + +# end_AB = data[graph_nmbr][counter][0]['stable_state'] +# end_BA = data[graph_nmbr][counter][1]['stable_state'] + +# #parents could be SS or (more likely) LC, so form superset +# if type(end_AB[0]) == list: +# end_AB = nv.superset(end_AB[0]) +# # end_AB = end_AB[0][0] +# else: +# end_AB = end_AB[0] +# if type(end_BA[0]) == list: +# end_BA = nv.superset(end_BA[0]) +# # end_BA = end_BA[0][0] +# else: +# end_BA = end_BA[0] + +# # track fraction of edges that are double positive (mutually beneficial) +# edge_count_AB = 0 +# edge_count_BA = 0 +# dbl_positive_count_AB = 0 +# dbl_positive_count_BA = 0 +# for edge in graph.edges(): + +# # need to translate from e.g. 'po_1' to position in the state list +# # indices 0-49 are plants 0-49 and indices 50-199 are pollinators 0-149 + +# # first node +# if edge[0][:2] == 'pl': edge0_index = int(edge[0][3:]) +# else: edge0_index = int(edge[0][3:])+50 + +# # second node +# if edge[1][:2] == 'pl': edge1_index = int(edge[1][3:]) +# else: edge1_index = int(edge[1][3:])+50 + +# # we only care about edges between present species +# if start_AB[edge0_index] == '1' and start_AB[edge1_index] == '1': +# edge_count_AB += 1 +# # and track seperately the edges that are also double positive +# if graph[edge[0]][edge[1]]['data'] == 1.0 and graph[edge[1]][edge[0]]['data'] == 1.0: dbl_positive_count_AB += 1 + +# # repeat for BA invasion +# if start_BA[edge0_index] == '1' and start_BA[edge1_index] == '1': +# edge_count_BA += 1 +# if graph[edge[0]][edge[1]]['data'] == 1.0 and graph[edge[1]][edge[0]]['data'] == 1.0: dbl_positive_count_BA += 1 + +# # once we've gone through all the edges, record the fraction of double positive edges +# # Because start_state can be a single state from a LC, occasionally we pick up +# # a state with no edges whatsoever and the ratio here gives a divide-by-zero error. Just skip +# # and keep track of how many of these there are +# if edge_count_AB == 0: +# r_exception_count += 1 +# r_dbl_positive.append(0) +# else: +# r_dbl_positive.append(dbl_positive_count_AB/edge_count_AB) +# if edge_count_BA == 0: +# r_exception_count += 1 +# r_dbl_positive.append(0) +# else: r_dbl_positive.append(dbl_positive_count_BA/edge_count_BA) + +# # track # of invader species +# # start +# A_invaders = start_AB[:25]+start_AB[50:125] +# B_invaders = start_BA[25:50]+start_BA[125:200] +# r_invasive_species.append(A_invaders.count('1')) +# r_invasive_species.append(B_invaders.count('1')) + +# # end +# A_invaders_end = end_AB[:25]+end_AB[50:125] +# B_invaders_end = end_BA[25:50]+end_BA[125:200] +# r_invasive_species_end.append(A_invaders_end.count('1')) +# r_invasive_species_end.append(B_invaders_end.count('1')) + +# # compare the starting state (source + invaders) to final stable state +# r_overlap.append(nv.overlap(end_AB,start_AB)) +# r_overlap.append(nv.overlap(end_BA,start_BA)) + +# r_start_size.append(data[graph_nmbr][counter][0]['start_state'][0].count('1')) +# r_start_size.append(data[graph_nmbr][counter][1]['start_state'][0].count('1')) + +# r_end_size.append(end_AB.count('1')) +# r_end_size.append(end_BA.count('1')) + +# r_delta_size.append(100*(r_end_size[-2] - start_AB.count('1'))/start_AB.count('1')) +# r_delta_size.append(100*(r_end_size[-1] - start_BA.count('1'))/start_BA.count('1')) + +# r_nestedness.append(data[graph_nmbr][counter][0]['SN'][0][0]) +# r_nestedness.append(data[graph_nmbr][counter][1]['SN'][0][0]) + +# r_connectance.append(data[graph_nmbr][counter][0]['C'][0][0]) +# r_connectance.append(data[graph_nmbr][counter][1]['C'][0][0]) + +# r_nestedness_final.append(data[graph_nmbr][counter][0]['SN'][0][-1]) +# r_nestedness_final.append(data[graph_nmbr][counter][1]['SN'][0][-1]) + +# r_connectance_final.append(data[graph_nmbr][counter][0]['C'][0][-1]) +# r_connectance_final.append(data[graph_nmbr][counter][1]['C'][0][-1]) + +# # track % of native species that make it to final state +# A_native = start_AB[25:50]+start_AB[125:200] +# A_end = end_AB[25:50]+end_AB[125:200] + +# B_native = start_BA[:25] + start_BA[:125] +# B_end = end_BA[:25]+end_BA[:125] + +# A_native_cnt = A_native.count('1') +# A_end_cnt = 0 +# for i in range(len(A_native)): +# if A_native[i] == '1' and A_end[i] == '1': A_end_cnt += 1 +# r_native_survive.append(A_end_cnt/A_native_cnt) + +# B_native_cnt = B_native.count('1') +# B_end_cnt = 0 +# for i in range(len(B_native)): +# if B_native[i] == '1' and B_end[i] == '1': B_end_cnt += 1 +# r_native_survive.append(B_end_cnt/B_native_cnt) + + +# print('r exception fraction = {0:.5g}'.format(r_exception_count/(r_exception_count + len(r_dbl_positive)))) + +# TIME1 = time.time() +# print('Completed in {0:.2f} minutes elapsed.'.format((TIME1-TIME0)/60)) + +# End analysis =============================================================== + +# Dump or Load Above Analysis ================================================ +# Optionally dump this data to file (then load below to speed up figure tweaks) +# data = [r_overlap, r_start_size, r_end_size, r_delta_size, r_nestedness, r_connectance, r_invasive_species, r_nestedness_final, r_connectance_final, r_native_survive, r_dbl_positive, r_invasive_species_end, +# s_overlap, s_start_size, s_end_size, s_delta_size, s_nestedness, s_connectance, s_invasive_species, s_nestedness_final, s_connectance_final, s_native_survive, s_dbl_positive, s_invasive_species_end] +# f = open(r'attractors_analysis_data.data','wb') +# pickle.dump(data,f) +# f.close() + +#Load the above (using destination LC superset) +f = open(r'attractors_analysis_data.data','rb') +data = pickle.load(f) +f.close() + +r_overlap, r_start_size, r_end_size, r_delta_size, r_nestedness, r_connectance, r_invasive_species, r_nestedness_final, r_connectance_final, r_native_survive, r_dbl_positive, r_invasive_species_end,\ + s_overlap, s_start_size, s_end_size, s_delta_size, s_nestedness, s_connectance, s_invasive_species, s_nestedness_final, s_connectance_final, s_native_survive, s_dbl_positive, s_invasive_species_end = data + +# End Dump or Load ============================================================ + + +# Data Processing / Figure Generation ========================================= + +# # Jaccard Index --------------------------------- +plt.figure(figsize=(8,4)) +ax_J1 = plt.subplot(121) +h_J1 = ax_J1.hist2d(s_start_size,s_overlap,bins=40,cmap = 'Purples', cmin = 0, vmin=1, vmax = 6000) +ax_J1.set_xlabel('Amalgamated Community # Species') +ax_J1.set_ylabel('Amalgamated vs. Stable Community Jaccard Index') +ax_J1.set_title('Whole Community Invasion') +# plt.colorbar(h1[3],ax=ax1) + +ax_J2 = plt.subplot(122) +h_J2 = ax_J2.hist2d(r_start_size,r_overlap,bins=40,cmap = 'Purples', cmin = 0, vmin=1, vmax = 6000) +ax_J2.set_xlabel('Amalgamated Community # Species') +ax_J2.set_title('Random Species Invasion') + +cax_J = plt.axes([0.9,0.146,0.01,0.912-0.146]) +plt.colorbar(h_J2[3],cax=cax_J) + +plt.subplots_adjust(top=0.912, bottom=0.146, left=0.083, right=0.881, hspace=0.2, wspace=0.139) + +# Stats +print('Jaccard:') +print('\tStable Invaders average: {0:.3f}'.format(sum(s_overlap)/len(s_overlap))) +print('\tRandom Invaders average: {0:.3f}'.format(sum(r_overlap)/len(r_overlap))) + +# # # Size Change ----------------------------------- +plt.figure(figsize=(8,4)) +ax_S1 = plt.subplot(121) +h_S1 = ax_S1.hist2d(s_start_size,s_end_size,bins=40,cmap = 'summer', cmin=1) +ax_S1.set_aspect('equal') +ax_S1.plot([0,80],[0,80],'k-') +ax_S1.set_xlim(0,80) +ax_S1.set_ylim(0,80) + +ax_S1.set_xlabel('Amalgamated Community # Species') +ax_S1.set_ylabel('Stable Community # Species') + +ax_S1.hist2d(s_start_size,s_end_size,bins=40,range = ([0,80],[0,80]), cmap = 'Greens', cmin = 0, vmin=0,vmax=15000) +ax_S1.set_aspect('equal') +ax_S1.plot([0,80],[0,80],'k-') +ax_S1.set_xlim(0,80) +ax_S1.set_ylim(0,80) +ax_S1.set_title('Whole Community Invasion') + +ax_S2 = plt.subplot(122) +ax_S2.set_xlabel('Amalgamated Community # Species') +ax_S2.set_aspect('equal') +ax_S2.plot([0,80],[0,80],'k-') +ax_S2.set_xlim(0,80) +ax_S2.set_ylim(0,80) + +h_S2 = ax_S2.hist2d(r_start_size,r_end_size,bins=40, range = ([0,80],[0,80]), cmap = 'Greens', cmin = 0, vmin=0,vmax=15000) +ax_S2.set_aspect('equal') +ax_S2.plot([0,80],[0,80],'k-') +ax_S2.set_xlim(0,80) +ax_S2.set_ylim(0,80) +ax_S2.set_title('Random Species Invasion') + +cax_S = plt.axes([0.9,0.146,0.01,0.912-0.146]) +plt.colorbar(h_S2[3],cax=cax_S) + +plt.subplots_adjust(top=0.912, bottom=0.146, left=0.083, right=0.881, hspace=0.2, wspace=0.139) + +# Stats +s_size_compare = [] +r_size_compare = [] +for i in range(len(s_start_size)): + s_size_compare.append(100*(s_end_size[i]-s_start_size[i])/s_start_size[i]) +for i in range(len(r_start_size)): + r_size_compare.append(100*(r_end_size[i]-r_start_size[i])/r_start_size[i]) + +print('Size comparison:') +print('\tStable Invaders average % change: {0:.3f}'.format(sum(s_size_compare)/len(s_size_compare))) +print('\tRandom Invaders average % change: {0:.3f}'.format(sum(r_size_compare)/len(r_size_compare))) + +# # # Connectance and Nestedness vs size -------------------- +# plt.figure(figsize=(8,8)) +# ax_CN1 = plt.subplot(221) +# h_CN1 = ax_CN1.hist2d(s_nestedness,s_end_size,bins=40,cmap = 'summer', cmin=1, vmin = 1, vmax = 10000) +# ax_CN1.set_xlabel('Amalgamated Community Nestedness') +# ax_CN1.set_ylabel('Stable Community # Species') +# ax_CN1.set_title('Whole Community Invasion') +# ax_CN1.set_ylim(0,80) +# ax_CN1.set_xlim(0,7) + +# ax_CN2 = plt.subplot(222) +# h_CN2 = ax_CN2.hist2d(r_nestedness,r_end_size,bins=40,cmap = 'summer', cmin=1, vmin = 1, vmax = 10000 ) +# ax_CN2.set_xlabel('Amalgamated Community Nestedness') +# ax_CN2.set_title('Random Species Invasion') +# ax_CN2.set_ylim(0,80) +# ax_CN2.set_xlim(0,7) + +# cax_CN12 = plt.axes([0.9,0.4825+0.073,0.01,0.91*(0.956-0.073)/2]) +# plt.colorbar(h_CN2[3],cax=cax_CN12) + + +# ax_CN3 = plt.subplot(223) +# h_CN3 = ax_CN3.hist2d(s_connectance,s_end_size,bins=40,cmap = 'autumn', cmin=1, vmin = 1, vmax = 50000) +# ax_CN3.set_xlabel('Amalgamated Community Connectance') +# ax_CN3.set_ylabel('Stable Community # Species') +# ax_CN3.set_ylim(0,80) +# ax_CN3.set_xlim(0,1) + +# ax_CN4 = plt.subplot(224) +# h_CN4 = ax_CN4.hist2d(r_connectance,r_end_size,bins=40,cmap = 'autumn', cmin=1, vmin = 1, vmax = 50000) +# ax_CN4.set_xlabel('Amalgamated Community Connectance') +# ax_CN4.set_ylim(0,80) +# ax_CN3.set_xlim(0,1) + +# cax_CN34 = plt.axes([0.9,0.073,0.01,0.91*(0.956-0.073)/2]) +# plt.colorbar(h_CN4[3],cax=cax_CN34) + +# plt.subplots_adjust(top=0.956, bottom=0.073, left=0.103, right=0.887, hspace=0.198, wspace=0.24) + +# # Stats +CN_s1 = spearmanr(s_nestedness, s_end_size) +CN_s2 = spearmanr(r_nestedness, r_end_size) +CN_s3 = spearmanr(s_connectance, s_end_size) +CN_s4 = spearmanr(r_connectance, r_end_size) + +print('Nestedness vs. Size Spearman correlation, pval:') +print('\tWhole: {0:.3g} \t {1:.3g}'.format(CN_s1[0],CN_s1[1])) +print('\tRandom: {0:.3g} \t {1:.3g}'.format(CN_s2[0],CN_s2[1])) + +print('Connectance vs. Size Spearman correlation, pval:') +print('\tWhole: {0:.3g} \t {1:.3g}'.format(CN_s3[0],CN_s3[1])) +print('\tRandom: {0:.3g} \t {1:.3g}'.format(CN_s4[0],CN_s4[1])) + + +# # # # Connectance and Nestedness pre vs post ------ +# plt.figure(figsize=(8,8)) +# plt.subplots_adjust(top=0.956, bottom=0.073, left=0.083, right=0.887, hspace=0.198, wspace=0.178) +# ax_CNb1 = plt.subplot(221) +# h_CNb1 = ax_CNb1.hist2d(s_nestedness,s_nestedness_final,bins=40,cmap = 'summer', cmin=1, vmin = 1, vmax = 20000) +# ax_CNb1.plot([0,0],[7,7],'k-',lw=2) + +# ax_CNb1.set_xlabel('Amalgamated Community Nestedness') +# ax_CNb1.set_ylabel('Stable Community Nestedness') +# ax_CNb1.set_title('Whole Community Invasion') +# ax_CNb1.set_ylim(0,7) +# ax_CNb1.set_xlim(0,7) + +# ax_CNb2 = plt.subplot(222) +# h_CNb2 = ax_CNb2.hist2d(r_nestedness,r_nestedness_final,bins=40,cmap = 'summer', cmin=1, vmin = 1, vmax = 20000 ) +# ax_CNb2.plot([0,0],[7,7],'k-',lw=2) +# ax_CNb2.plot([0,0],[7,7],'k-',lw=2) +# ax_CNb2.set_xlabel('Amalgamated Community Nestedness') +# ax_CNb2.set_title('Random Species Invasion') +# ax_CNb2.set_ylim(0,7) +# ax_CNb2.set_xlim(0,7) + +# box = ax_CNb2.get_position() #box.bounds will give (x0, y0, width, height) +# cax_CNb2 = plt.axes([0.9,box.bounds[1],0.01,box.bounds[3]]) +# plt.colorbar(h_CNb2[3],cax=cax_CNb2) + +# ax_CNb3 = plt.subplot(223) +# ax_CNb3.plot([0,0],[1,1],'k-') +# h_CNb3 = ax_CNb3.hist2d(s_connectance,s_connectance_final,bins=40,cmap = 'autumn', cmin=1, vmin = 1, vmax = 90000) +# ax_CNb3.set_xlabel('Amalgamated Community Connectance') +# ax_CNb3.set_ylabel('Stable Community Connectance') +# # ax_CNb3.set_ylim(-100,1200) + +# ax_CNb4 = plt.subplot(224) +# h_CNb4 = ax_CNb4.hist2d(r_connectance,r_connectance_final,bins=40,cmap = 'autumn', cmin=1, vmin = 1, vmax = 90000) +# ax_CNb4.plot([0,0],[1,1],'k-') +# ax_CNb4.set_xlabel('Amalgamated Community Connectance') +# # ax_CNb4.set_ylim(-100,1200) + +# box = ax_CNb4.get_position() #box.bounds will give (x0, y0, width, height) +# cax_CNb4 = plt.axes([0.9,box.bounds[1],0.01,box.bounds[3]]) +# plt.colorbar(h_CNb4[3],cax=cax_CNb4) + +print('Nestedness pre vs. post:') +print('\tWhole\t pre: {0:.3g} \t post: {1:.3g}'.format(sum(s_nestedness)/len(s_nestedness),sum(s_nestedness_final)/len(s_nestedness_final))) +print('\2 sample KS test: {0:.3g}'.format(ks_2samp(s_nestedness,s_nestedness_final)[1])) +print('\tRandom\t pre: {0:.3g} \t post: {1:.3g}'.format(sum(r_nestedness)/len(r_nestedness),sum(r_nestedness_final)/len(r_nestedness_final))) +print('\2 sample KS test: {0:.3g}'.format(ks_2samp(r_nestedness,r_nestedness_final)[1])) + +print('Connectance pre vs. post:') +print('\tWhole\t pre: {0:.3g} \t post: {1:.3g}'.format(sum(s_connectance)/len(s_connectance),sum(s_connectance_final)/len(s_connectance_final))) +print('\2 sample KS test: {0:.3g}'.format(ks_2samp(s_connectance,s_connectance_final)[1])) +print('\tRandom\t pre: {0:.3g} \t post: {1:.3g}'.format(sum(r_connectance)/len(r_connectance),sum(r_connectance_final)/len(r_connectance_final))) +print('\2 sample KS test: {0:.3g}'.format(ks_2samp(r_connectance,r_connectance_final)[1])) + + +# # Size (detailed) ------------------------------- +# fig_SD = plt.figure(figsize=(8,9)) +# plt.subplots_adjust(top=0.977,bottom=0.065,left=0.1,right=0.89,hspace=0.276,wspace=0.22) + +# # top two panels for overall # species vs final % species +# ax_SD1 = plt.subplot(321) +# h_SD1 = ax_SD1.hist2d(s_start_size,s_delta_size,bins=40,cmap = 'summer', cmin=1, vmin = 1, vmax = 35000) +# ax_SD1.set_xlabel('Amalgamated State # Species') +# ax_SD1.set_xlim(0,60) +# ax_SD1.set_ylim(-100,1200) +# # ax_SD1.set_ylabel('Amalgamated to Child Change in # Species (%)') + +# ax_SD2 = plt.subplot(322) +# h_SD2 = ax_SD2.hist2d(r_start_size,r_delta_size,bins=40,cmap = 'summer', cmin=1, vmin = 1, vmax = 35000) +# ax_SD2.set_xlabel('Amalgamated State # Species') +# ax_SD2.set_xlim(0,60) +# ax_SD2.set_ylim(-100,1200) + +# box = ax_SD2.get_position() #box.bounds will give (x0, y0, width, height) +# cax_SD12 = plt.axes([0.9,box.bounds[1],0.01,box.bounds[3]]) +# plt.colorbar(h_SD2[3],cax=cax_SD12) + +# # middle two panels change horizontal axes to # invasive species +# ax_SD3 = plt.subplot(323) +# h_SD3 = ax_SD3.hist2d(s_invasive_species,s_delta_size,bins=40,cmap = 'autumn', cmin=1, vmin = 1, vmax = 35000) +# ax_SD3.set_xlabel('# Invasive Species') +# ax_SD3.set_xlim(0,60) +# ax_SD3.set_ylabel('Amalgamated to Stable Change in # Species (%)') +# ax_SD3.set_ylim(-100,1200) + +# ax_SD4 = plt.subplot(324) +# h_SD4 = ax_SD4.hist2d(r_invasive_species,r_delta_size,bins=40,cmap = 'autumn', cmin=1, vmin = 1, vmax = 35000) +# ax_SD4.set_xlabel('# Invasive Species') +# ax_SD4.set_xlim(0,60) +# ax_SD4.set_ylim(-100,1200) + +# box = ax_SD4.get_position() #box.bounds will give (x0, y0, width, height) +# cax_SD34 = plt.axes([0.9,box.bounds[1],0.01,box.bounds[3]]) +# plt.colorbar(h_SD4[3],cax=cax_SD34) + +# # bottom two panels change horizontal axes to # host species +s_xvals = [s_start_size[i] - s_invasive_species[i] for i in range(len(s_start_size))] +r_xvals = [r_start_size[i] - r_invasive_species[i] for i in range(len(r_start_size))] + +# ax_SD5 = plt.subplot(325) +# h_SD5 = ax_SD5.hist2d(s_xvals,s_delta_size,bins=40,cmap = 'winter', cmin=1, vmin = 1, vmax = 35000) +# ax_SD5.set_xlabel('# Invaded Species') +# ax_SD5.set_xlim(0,60) +# # ax_SD5.set_ylabel('Amalgamated to Child Change in # Species (%)') +# ax_SD5.set_ylim(-100,1200) + +# ax_SD6 = plt.subplot(326) +# h_SD6 = ax_SD6.hist2d(r_xvals,r_delta_size,bins=40,cmap = 'winter', cmin=1, vmin = 1, vmax = 35000) +# ax_SD6.set_xlabel('# Invaded Species') +# ax_SD6.set_xlim(0,60) +# ax_SD6.set_ylim(-100,1200) + +# box = ax_SD6.get_position() #box.bounds will give (x0, y0, width, height) +# cax_SD56 = plt.axes([0.9,box.bounds[1],0.01,box.bounds[3]]) +# plt.colorbar(h_SD6[3],cax=cax_SD56) + +# # Stats +SD_s1 = spearmanr(s_start_size, s_delta_size) +SD_s2 = spearmanr(r_start_size, r_delta_size) +SD_s3 = spearmanr(s_invasive_species, s_delta_size) +SD_s4 = spearmanr(r_invasive_species, r_delta_size) +SD_s5 = spearmanr(s_xvals, s_delta_size) +SD_s6 = spearmanr(r_xvals, r_delta_size) + +print('Size (detailed) Spearman correlation, pval:') +print('\tpanel 1: {0:.3g} \t {1:.3g}\tpanel 2: {2:.3g} \t {3:.3g}'.format(SD_s1[0],SD_s1[1],SD_s2[0],SD_s2[1])) +print('\tpanel 3: {0:.3g} \t {1:.3g}\tpanel 4: {2:.3g} \t {3:.3g}'.format(SD_s3[0],SD_s3[1],SD_s4[0],SD_s4[1])) +print('\tpanel 5: {0:.3g} \t {1:.3g}\tpanel 6: {2:.3g} \t {3:.3g}'.format(SD_s5[0],SD_s5[1],SD_s6[0],SD_s6[1])) + + +# # Size (detailed v2) -------------------------- +# # same vertical axis as earlier Size figure +# fig_SD = plt.figure(figsize=(8,9)) +# plt.subplots_adjust(top=0.977,bottom=0.065,left=0.1,right=0.89,hspace=0.276,wspace=0.22) + +# # top two panels for overall # species vs final % species +# ax_SD1 = plt.subplot(321) +# h_SD1 = ax_SD1.hist2d(s_start_size,s_end_size,bins=40,cmap = 'summer', cmin=1, vmin = 1, vmax = 15000) +# ax_SD1.set_xlabel('Amalgamated State # Species') +# ax_SD1.set_xlim(0,80) +# ax_SD1.set_ylim(0,80) +# ax_SD1.plot([0,80],[0,80],'k-') +# # ax_SD1.set_ylabel('Amalgamated to Child Change in # Species (%)') + +# ax_SD2 = plt.subplot(322) +# h_SD2 = ax_SD2.hist2d(r_start_size,r_end_size,bins=40,cmap = 'summer', cmin=1, vmin = 1, vmax = 15000) +# ax_SD2.set_xlabel('Amalgamated State # Species') +# ax_SD2.set_xlim(0,80) +# ax_SD2.set_ylim(0,80) +# ax_SD2.plot([0,80],[0,80],'k-') + +# box = ax_SD2.get_position() #box.bounds will give (x0, y0, width, height) +# cax_SD12 = plt.axes([0.9,box.bounds[1],0.01,box.bounds[3]]) +# plt.colorbar(h_SD2[3],cax=cax_SD12) + +# # middle two panels change horizontal axes to # invasive species +# ax_SD3 = plt.subplot(323) +# h_SD3 = ax_SD3.hist2d(s_invasive_species,s_end_size,bins=40,cmap = 'autumn', cmin=1, vmin = 1, vmax = 8000) +# ax_SD3.set_xlabel('# Invasive Species') +# ax_SD3.set_xlim(0,80) +# ax_SD3.set_ylabel('Stable Community # Species') +# ax_SD3.set_ylim(0,80) +# ax_SD3.plot([0,80],[0,80],'k-') + +# ax_SD4 = plt.subplot(324) +# h_SD4 = ax_SD4.hist2d(r_invasive_species,r_end_size,bins=40,cmap = 'autumn', cmin=1, vmin = 1, vmax = 8000) +# ax_SD4.set_xlabel('# Invasive Species') +# ax_SD4.set_xlim(0,80) +# ax_SD4.set_ylim(0,80) +# ax_SD4.plot([0,80],[0,80],'k-') + +# box = ax_SD4.get_position() #box.bounds will give (x0, y0, width, height) +# cax_SD34 = plt.axes([0.9,box.bounds[1],0.01,box.bounds[3]]) +# plt.colorbar(h_SD4[3],cax=cax_SD34) + +# # bottom two panels change horizontal axes to # host species +s_xvals = [s_start_size[i] - s_invasive_species[i] for i in range(len(s_start_size))] +r_xvals = [r_start_size[i] - r_invasive_species[i] for i in range(len(r_start_size))] + +# ax_SD5 = plt.subplot(325) +# h_SD5 = ax_SD5.hist2d(s_xvals,s_end_size,bins=40,cmap = 'winter', cmin=1, vmin = 1, vmax = 8000) +# ax_SD5.set_xlabel('# Invaded Species') +# ax_SD5.set_xlim(0,80) +# # ax_SD5.set_ylabel('Amalgamated to Child Change in # Species (%)') +# ax_SD5.set_ylim(0,80) +# ax_SD5.plot([0,80],[0,80],'k-') + +# ax_SD6 = plt.subplot(326) +# h_SD6 = ax_SD6.hist2d(r_xvals,r_end_size,bins=40,cmap = 'winter', cmin=1, vmin = 1, vmax = 8000) +# ax_SD6.set_xlabel('# Invaded Species') +# ax_SD6.set_xlim(0,80) +# ax_SD6.set_ylim(0,80) +# ax_SD6.plot([0,80],[0,80],'k-') + +# box = ax_SD6.get_position() #box.bounds will give (x0, y0, width, height) +# cax_SD56 = plt.axes([0.9,box.bounds[1],0.01,box.bounds[3]]) +# plt.colorbar(h_SD6[3],cax=cax_SD56) + +# # Stats +SD_s1 = spearmanr(s_start_size, s_end_size) +SD_s2 = spearmanr(r_start_size, r_end_size) +SD_s3 = spearmanr(s_invasive_species, s_end_size) +SD_s4 = spearmanr(r_invasive_species, r_end_size) +SD_s5 = spearmanr(s_xvals, s_end_size) +SD_s6 = spearmanr(r_xvals, r_end_size) + +print('Size (detailed v2) Spearman correlation, pval:') +print('\tpanel 1: {0:.3g} \t {1:.3g}\tpanel 2: {2:.3g} \t {3:.3g}'.format(SD_s1[0],SD_s1[1],SD_s2[0],SD_s2[1])) +print('\tpanel 3: {0:.3g} \t {1:.3g}\tpanel 4: {2:.3g} \t {3:.3g}'.format(SD_s3[0],SD_s3[1],SD_s4[0],SD_s4[1])) +print('\tpanel 5: {0:.3g} \t {1:.3g}\tpanel 6: {2:.3g} \t {3:.3g}'.format(SD_s5[0],SD_s5[1],SD_s6[0],SD_s6[1])) + + + +# % of species that survive from native community +# plt.figure() +# ax_ss1 = plt.subplot(121) +# ax_ss1.hist(s_native_survive) +# ax_ss1.set_xlabel('Whole') + +# ax_ss2 = plt.subplot(122) +# ax_ss2.hist(r_native_survive) +# ax_ss2.set_xlabel('Random') + +print('Native species survival:') +print('\t Whole Community: {0:.3g}'.format(sum(s_native_survive)/len(s_native_survive))) +print('\t Random Community: {0:.3g}'.format(sum(r_native_survive)/len(r_native_survive))) + +# # Fraction of double-positive edges vs. size of stable community +# fig_DPS = plt.figure() +# # plt.subplots_adjust(top=0.977,bottom=0.065,left=0.1,right=0.89,hspace=0.276,wspace=0.22) + + +# ax_DPS1 = plt.subplot(121) +# h_DPS1 = ax_DPS1.hist2d([100*x for x in s_dbl_positive],s_end_size,bins=40,cmap = 'summer', cmin=1, vmin = 1, vmax = 15000) + +# ax_DPS1.set_xlabel('Mutualistic Interactions (%)') +# ax_DPS1.set_ylabel('Stable Community # Species') +# ax_DPS1.set_title('Whole Community Invasion') +# # ax_DPS1.set_xlim(0,80) +# # ax_DPS1.set_ylim(0,80) + +# ax_DPS2 = plt.subplot(122) +# h_DPS2 = ax_DPS2.hist2d([100*x for x in r_dbl_positive],r_end_size,bins=40,cmap = 'summer', cmin=1, vmin = 1, vmax = 15000) + +# ax_DPS2.set_xlabel('Mutualistic Interactions (%)') +# ax_DPS2.set_title('Random Species Invasion') + +# # ax_DPS1.set_xlim(0,80) +# # ax_DPS1.set_ylim(0,80) + +DPS_s1 = spearmanr(s_dbl_positive, s_end_size) +DPS_s2 = spearmanr(r_dbl_positive, r_end_size) + +print('Mutualisms vs. end size; Spearman correlation, pval:') +print('\tWhole Community: {0:.3g}'.format(DPS_s1[0])) +print('\Random Species: {0:.3g}'.format(DPS_s2[0])) + +DPJ_s1 = spearmanr(s_dbl_positive, s_overlap) +DPJ_s2 = spearmanr(r_dbl_positive, r_overlap) + +print('Mutualisms vs. Jaccard; Spearman correlation, pval:') +print('\tWhole Community: {0:.3g}'.format(DPS_s1[0])) +print('\Random Species: {0:.3g}'.format(DPS_s2[0])) + +# Balance of native vs. invasive species pre vs post ------ +s_balance_pre = [] +s_balance_post = [] +s_exception_count = 0 +for i in range(len(s_invasive_species)): + native_start = s_start_size[i] - s_invasive_species[i] + inv_start = s_invasive_species[i] + + native_end = s_end_size[i] - s_invasive_species_end[i] + inv_end = s_invasive_species_end[i] + + + if native_start == 0 or inv_start == 0 or native_end == 0 or inv_end == 0: + s_exception_count += 1 + continue + + if native_start >= inv_start: + # just report the ratio as 1, 1.5, 2, or whatever, but shift "1" to value + # of "0" on the axis for ease of plotting + s_balance_pre.append(native_start/inv_start-1) + else: + # transform the ratio so we get equal spacing for "majority inv" and + # "majority native"; so native/invader = 0.5 is flipped and recorded + # as 2, then shifted to -1 so 0.5 is the same distance from 1:1 (plotted at 0) + # as a ratio of 2 + s_balance_pre.append(-inv_start/native_start + 1) + + # repeat for ending + if native_end >= inv_end: + s_balance_post.append(native_end/inv_end-1) + else: + s_balance_post.append(-inv_end/native_end + 1) + +print('Species balance, whole community exception count = {0} vs. {1} plotted.'\ + .format(s_exception_count,len(s_balance_pre) + len(s_balance_post))) + +# Repeat all of above for random invasion +r_balance_pre = [] +r_balance_post = [] +r_exception_count = 0 +for i in range(len(r_invasive_species)): + native_start = r_start_size[i] - r_invasive_species[i] + inv_start = r_invasive_species[i] + + native_end = r_end_size[i] - r_invasive_species_end[i] + inv_end = r_invasive_species_end[i] + + + if native_start == 0 or inv_start == 0 or native_end == 0 or inv_end == 0: + s_exception_count += 1 + continue + + if native_start >= inv_start: + # just report the ratio as 1, 1.5, 2, or whatever, centered so 1:1 is at 0 + r_balance_pre.append(native_start/inv_start - 1) + else: + # transform the ratio so we get equal spacing for "majority inv" and + # "majority native"; so native/invader = 0.5 is flipped and recorded + # as 2, then shifted to -1 so 0.5 is the same distance from 1 as a ratio + # of 2 + r_balance_pre.append(-inv_start/native_start + 1) + + # repeat for ending + if native_end >= inv_end: + r_balance_post.append(native_end/inv_end - 1) + else: + r_balance_post.append(-inv_end/native_end + 1) + +print('Species balance, random community exception count = {0} vs. {1} plotted.'\ + .format(r_exception_count,len(r_balance_pre) + len(r_balance_post))) + +# Now plot +plt.figure(figsize=(8,4)) +ax_bal1 = plt.subplot(121) +h_bal1 = ax_bal1.hist2d(s_balance_pre,s_balance_post,bins=40,cmap = 'Reds', range = [[-6,6],[-6,6]], vmin=1, vmax = 12000) + +ax_bal1.plot([-6,6],[-6,6], ls='-', c='k', lw=0.5) +ax_bal1.plot([-6,6],[0,0], ls='-', c='k', lw=0.5) +ax_bal1.plot([0,0],[-6,6], ls='-', c='k', lw=0.5) + +ax_bal1.set_xlabel('Amalgamated Community\nNative Species:Invasive Species') +ax_bal1.set_ylabel('Final Stable Community\nNative Species:Invasive Species') + +ax_bal1.set_xticks([-6,-4,-2,0,2,4,6]) +ax_bal1.set_xticklabels(['1:7','1:5','1:3','1:1','3:1','5:1','7:1'],fontsize='small') + +ax_bal1.set_yticks([-6,-4,-2,0,2,4,6]) +ax_bal1.set_yticklabels(['1:7','1:5','1:3','1:1','3:1','5:1','7:1'],fontsize='small') + +ax_bal1.set_title('Whole Community Invasion') +# plt.colorbar(h1[3],ax=ax1) + +ax_bal2 = plt.subplot(122) +h_bal2 = ax_bal2.hist2d(r_balance_pre,r_balance_post,bins=40,cmap = 'Reds', range = [[-6,6],[-6,30]], vmin=1, vmax = 12000) + +ax_bal2.plot([-6,30],[-6,30], ls='-', c='k', lw=0.5) +ax_bal2.plot([-6,6],[0,0], ls='-', c='k', lw=0.5) +ax_bal2.plot([0,0],[-6,30], ls='-', c='k', lw=0.5) + +ax_bal2.set_xlabel('Amalgamated Community\nNative Species:Invasive Species') +# ax_bal2.set_ylabel('Final Stable Community\nNative Species:Invasive Species') + +ax_bal2.set_xticks([-6,-4,-2,0,2,4,6]) +ax_bal2.set_xticklabels(['1:7','1:5','1:3','1:1','3:1','5:1','7:1'],fontsize='small') + +ax_bal2.set_yticks([-6,0,6,13,20,27]) +ax_bal2.set_yticklabels(['1:7','1:1','7:1','14:1','21:1','28:1'],fontsize='small') + +ax_bal2.set_title('Random Species Invasion') + +cax_bal = plt.axes([0.92,0.16,0.01,0.91-0.16]) +plt.colorbar(h_bal2[3],cax=cax_bal) + +plt.subplots_adjust(top=0.91, +bottom=0.16, +left=0.11, +right=0.9, +hspace=0.2, +wspace=0.25) + +# report correlation b/w dbl positive edges and # of invasive species +DP_s_inv = spearmanr(s_dbl_positive, s_invasive_species) +DP_r_inv = spearmanr(r_dbl_positive, r_invasive_species) + +print('Mutualisms vs. Nmbr Invasive Species Spearman correlation, pval:') +print('\tWhole: {0:.3g} \t {1:.3g}'.format(DP_s_inv[0],DP_s_inv[1])) +print('\tRandom: {0:.3g} \t {1:.3g}'.format(DP_r_inv[0],DP_r_inv[1])) + +# Repeat for # native species, but need to generate that list first +s_native = [] +r_native = [] +for i in range(len(r_invasive_species)): + s_native.append(s_start_size[i] - s_invasive_species[i]) + r_native.append(r_start_size[i] - r_invasive_species[i]) + +DP_s_nat = spearmanr(s_dbl_positive, s_native) +DP_r_nat = spearmanr(r_dbl_positive, r_native) + +print('Mutualisms vs. Nmbr Native Species Spearman correlation, pval:') +print('\tWhole: {0:.3g} \t {1:.3g}'.format(DP_s_nat[0],DP_s_nat[1])) +print('\tRandom: {0:.3g} \t {1:.3g}'.format(DP_r_nat[0],DP_r_nat[1])) + + +# Final draw command ---------------------------- +plt.show() diff --git a/create_graphs.py b/create_graphs.py new file mode 100644 index 0000000..1eaaa1b --- /dev/null +++ b/create_graphs.py @@ -0,0 +1,35 @@ +''' +Generates given # of species interaction graphs. + +Python 3.x + +Colin Campbell +campbeco@mountunion.edu +''' + +import function_library as nv +import numpy as np +import pickle +import time + +def graph_generator(x): + graph = nv.form_network(plants=50,pollinators=150) + if x%10==0: + T2 = time.time() + print('Generated graph {0} at {1:.2f} minutes elapsed.'.format(x,(T2-T1)/60.0)) + return graph + + +T1 = time.time() + +graphs=[graph_generator(i) for i in range(1000)] + +#Save graphs +f=open('graphlist.data','wb') +pickle.dump(graphs,f) +f.close() + +T2 = time.time() + +print('{0:.2f} minutes for algorithm to elapse.'.format(((T2-T1)/60))) + diff --git a/find_attractors.py b/find_attractors.py new file mode 100644 index 0000000..741aad6 --- /dev/null +++ b/find_attractors.py @@ -0,0 +1,115 @@ +''' +Here we consider our 1000 interaction graphs. For each, we: + 1. choose 100 random starting communities + 2. fix half of the plants and half of the pollinators to be OFF + 3. Allow each to dynamically evolve to a stable community, C + 4. return to steps 2-3 with the other half of the species fixed OFF + +For each iteration, we record: + 1. (list) spectral nestedness at each time step + 2. (list) # of present species at each time step + 3. (list) connectance at each time step + 4. (int) the time step where the graph initially stabilizes + + +Master workflow: + 1. What are the species interaction graphs? (create_graphs.py) + 2. What are the stable communities in these graphs, if we segregate them + in to two "half communities" to enforce isolated communities that + nonetheless have inter-community interactions for when invasions occur? + (THIS FILE) + +The file is structured to run either on a local machine or submitted to a +High-Performance Computing center. See comments. +(This research was performed via the Ohio Supercomputer Center, osc.edu) + +Python 3.x + +Colin Campbell +campbeco@mountunion.edu +''' + +import function_library as nv +import networkx as nx +import numpy as np +import pickle +import time +import os + +TIME1 = time.time() +#PBS_ITERATOR=int(os.environ["PBS_ARRAYID"]) #We use this variable to slice up a submission across an array on a HPC. Use: qsub -t 0-99 pbs_script.sh +PBS_ITERATOR=0 #Or dummy variable for local debugging + +# Load graphs from memory ------------------------------------------------------ +#Local +f = open(r'graphlist.data','rb') +# OSC +#f = open(r'/path/to/graphlist.data','rb') +graphs = pickle.load(f) +f.close() + +graphs=graphs[10*PBS_ITERATOR:(10*PBS_ITERATOR)+10] + +N = nx.number_of_nodes(graphs[0]) +nodes = nv.node_sorter(graphs[0].nodes()) +data = {} + +# Perform analysis ------------------------------------------------------------- +for graph_index in range(len(graphs)): + TIME2 = time.time() + print('Beginning graph {0} at {1:.2f} minutes elapsed.'.format(graph_index,(TIME2-TIME1)/60)) + data[10*PBS_ITERATOR+graph_index]={} + for start_index in range(100): + # if start_index % 10 == 0: + # TIME2 = time.time() + # print('\tBeginning iteration {0} at {1:.2f} minutes elapsed.'.format(start_index,(TIME2-TIME1)/60)) + #print('Start index = {0}'.format(start_index)) + #Initialize graph and find initial attractor + #inner dictionaries that will hold the data for this iteration + inner_data_A={'start_state':[],'stable_state':[],'sys_state':[],'SN':[],'C':[],'size':[],'benchmark':0,'adj':[]} + inner_data_B={'start_state':[],'stable_state':[],'sys_state':[],'SN':[],'C':[],'size':[],'benchmark':0,'adj':[]} + start_state = ''.join(['0' if np.random.randint(2) == 0 else '1' for x in range(N)]) #choose a random state for this graph + + # trial A with first half of the species forced OFF; trial B with second half + # (assume 50 plants, 150 pollinators, so 25 and 75 of each, respectively) + off_list_A = list(range(25))+list(range(50,125)) + off_list_B = list(range(25,50))+list(range(125,200)) + + # execute trial and store data for run A + G = graphs[graph_index].copy() + G,sys_state = nv.force_state(G,start_state) + stable_state,sys_state,SN_vals,C_vals,size_vals,adj1 = nv.check_sufficiency_mod(G,sys_state,posweight=4,ind=off_list_A,TRANS_TYPE=1, FORCE_EXTINCT=True,extra_data=True) #determine the resulting stable configuration + + inner_data_A['start_state'].append(start_state) + inner_data_A['stable_state'].append(stable_state) + inner_data_A['sys_state'].append(sys_state) + inner_data_A['SN'].append(SN_vals) + inner_data_A['C'].append(C_vals) + inner_data_A['size'].append(size_vals) + inner_data_A['benchmark']=len(size_vals) + + # execute trial and store data for run B + G = graphs[graph_index].copy() + G,sys_state = nv.force_state(G,start_state) + stable_state,sys_state,SN_vals,C_vals,size_vals,adj1 = nv.check_sufficiency_mod(G,sys_state,posweight=4,ind=off_list_B,TRANS_TYPE=1, FORCE_EXTINCT=True,extra_data=True) #determine the resulting stable configuration + + inner_data_B['start_state'].append(start_state) + inner_data_B['stable_state'].append(stable_state) + inner_data_B['sys_state'].append(sys_state) + inner_data_B['SN'].append(SN_vals) + inner_data_B['C'].append(C_vals) + inner_data_B['size'].append(size_vals) + inner_data_B['benchmark']=len(size_vals) + + #Finally, append this data to the master dictionary + data[10*PBS_ITERATOR+graph_index][start_index]=[inner_data_A.copy(),inner_data_B.copy()] + +# Write to file ---------------------------------------------------------------- +f = open(r'attractors_'+str(PBS_ITERATOR)+'.data','wb') +#f = open(r'/path/to/attractors_'+str(PBS_ITERATOR)+'.data','wb') +pickle.dump(data,f) +f.close() + +TIME2 = time.time() +print('Script completed in {0:.2f} minutes.'.format(((TIME2-TIME1)/60))) + diff --git a/function_library.py b/function_library.py new file mode 100644 index 0000000..2136f21 --- /dev/null +++ b/function_library.py @@ -0,0 +1,985 @@ +''' +This module contains the 'master functions' used to generate and analyze +networks in this project. + +Python 3.x + +Colin Campbell +campbeco@mountunion.edu +''' + +import numpy as np +import scipy.special as sp +import math + +import networkx as nx +from networkx.algorithms import bipartite as bp # needed for configuration_model + + +def form_network(plants=10,pollinators=10,\ +po_mean=3.5,po_stdev=3.272565,po_skewness=1.43,\ +pl_mean=3.5,pl_stdev=2.579218,pl_skewness=.73,\ +po_decay=1.273077,po_cutoff=.328301,\ +pl_decay=.876667,pl_cutoff=.141464): + ''' + OVERVIEW: + Forms a bipartite network of plants and pollinators. Plants and pollinators + are: + 1. Assigned characteristic proboscis/stem lengths from a skewed normal + distribution. + 2. Assigned edges (inter-type only) via an exponentially cut off power law + distribution. + The edges are then made bidirectional, with weights according to the general + schematic: + plants prefer prob<=stem + pollinators prefer prob>=stem + Edge weights are assigned values of -1, or 1. A plant->pollinator edge is + assigned -1 if probstem (good), or 1 if + prob=stem (preferred). We consider prob=stem if the % difference is <=10%. + The percent difference b/w x,y, is defined as 100*abs(x-y)/avg(x,y). + + + INPUT PARAMETERS: + * pollinators, plants = # of ea. type of species + + * mean, stdev, and skewness = mean, stdev, and skewness of distribution from + which the species' characteristic lengths will be drawn. alpha = 0 for + normal dist, - (+) for left (right) skew. + + The decay function used for assigning edges b/w plants and pollinators is: + p(x) = x^(-decay) * exp(-cutoff*x) + * decay = characteristic decay rate for power law + * cutoff = exponential cutoff term + + Default parameter values are drawn from: + * Stang et al, Annals of Botany 103: 1459-1469, 2009 (for length + distributions) + * Jordano et al, Ecology Letters 6: 69-81, 2003 (for decay function) + + PROCESS: + 1. Generate range of values from the specified distributions. + 2. Assign proboscis/stem lengths according to skewed normal distribution. + 3. Assign edges b/w pollinators and plants according to cut off exponential + distribution. + 4. Make edges bidirectional, and assign weights according to a 3-valued + fitting function. + ''' + + #1. Generate range of values from specified distributions + def cdf_skew(X,sgn): + ''' + Returns the cumulative density function of the skew normal distribution. + ''' + def dens(X,sgn): + ''' + Returns the probability density function of the skew normal distribution. + See http://en.wikipedia.org/wiki/Skew_normal_distribution for details. + sp.ndtr returns the area under the standard Gaussian probability density + function, integrated from minus infinity to (arg). + ''' + if sgn=='po': + #Compute scaling parameters from supplied information (see word doc for derivation) + po_skew = min([0.99,po_skewness]) # For lognormal the upper limit is 1; must assume sampling bias if value exceeds + arg1=2*po_skew/(4-math.pi) + delta=math.sqrt(math.pi/2.0)*arg1**(1/3)/((1+arg1**(2/3))**.5) + omega=math.sqrt((math.pi*po_stdev*po_stdev)/(math.pi-2*delta*delta)) + alpha=delta/(math.sqrt(1-delta*delta)) + epsilon=po_mean-omega*delta*math.sqrt(2/math.pi) + + arg2=(X-epsilon)/omega + Y = (2./po_stdev)*np.exp(-(arg2)**2/2)/np.sqrt(2*np.pi) + Y *= sp.ndtr(alpha*arg2) #The max value should be 1, but for some reason it is near 1.41424 + return Y #This is a normalization issue we account for below (I haven't tracked down the source) + elif sgn=='pl': + #Compute scaling parameters from supplied information (see word doc for derivation) + pl_skew = min([0.99,pl_skewness]) # For lognormal the upper limit is 1; must assume sampling bias if value exceeds + arg1=2*pl_skew/(4-math.pi) + delta=math.sqrt(math.pi/2.0)*arg1**(1/3)/((1+arg1**(2/3))**.5) + omega=math.sqrt((math.pi*pl_stdev*pl_stdev)/(math.pi-2*delta*delta)) + alpha=delta/(math.sqrt(1-delta*delta)) + epsilon=pl_mean-omega*delta*math.sqrt(2/math.pi) + + arg2=(X-epsilon)/omega + Y = (2./pl_stdev)*np.exp(-(arg2)**2/2)/np.sqrt(2*np.pi) + Y *= sp.ndtr(alpha*arg2) #The max value should be 1, but for some reason it is near 1.41424 + return Y #This is a normalization issue we account for below (I haven't tracked down the source) + + if sgn=='pl': + Y_pl = dens(X,'pl') + Y_pl = np.cumsum(Y_pl)*(X[1]-X[0]) + return Y_pl + elif sgn=='po': + Y_po = dens(X,'po') + Y_po = np.cumsum(Y_po)*(X[1]-X[0]) + return Y_po + + def cutoff_fn(sgn): + '''For power law/exponential cutoff: + We generate this with a standard exponential distribution, and force + the cutoff behavior by accepting each generated value with probability + p=(x/xmin)^(-cutoff). + Source for approach: + http://arxiv.org/PS_cache/arxiv/pdf/0706/0706.1062v2.pdf (pg.40) + + we assume xmin=1. (i.e. the lower cutoff for range where power law + applies is 1)''' + + if sgn=='pl': + decay,cutoff=pl_decay,pl_cutoff + elif sgn=='po': + decay,cutoff=po_decay,po_cutoff + + #Form exponential values + Y_cutoff=np.random.exponential(scale=1.0/cutoff,size=10000) + def cutoff_filter(x): + p=x**(-decay) + comparison=np.random.ranf() + return p>comparison + + #Remove some according to power law test + Y_cutoff=list(filter(cutoff_filter,Y_cutoff)) + + #Fill in empty values + temp=np.array([]) + while len(temp)<(10000-len(Y_cutoff)): + x=np.random.exponential(scale=1.0/cutoff) + p=x**(-decay) + comparison=np.random.ranf() + if p>comparison: + temp=np.append(temp,x) + + Y_cutoff=np.append(Y_cutoff,temp) + return Y_cutoff + + #Form x range for distributions + X_pl = np.arange(pl_mean-50, pl_mean+100, 0.1) #These ranges characterize an appropriate range for the default values + X_po = np.arange(po_mean-50, po_mean+100, 0.1) #If nondefault values are used, care should be taken that the leftmost... + #...value ~0 and the rightmost ~1 + #Form skew normal distribution + ''' + Procedure for generating values from a skew normal distribution: + 1. Form cdf of skew distribution w/ given params. + - Interpret values as CDF(x)=y; y=cdf value for input x + 2. Take a value, x1, from uniform distribution[0,1] + 3. Find largest value x2 s.t. CDF(x2)<=x1 + 4. x2 is the output from the end distribution. + ''' + Y_pl,Y_po = cdf_skew(X_pl,'pl'),cdf_skew(X_po,'po') + Y_pl_skew,Y_po_skew=[],[] + Y_pl_max=max(Y_pl) #The normalization factors alluded to in the cdf function defn + Y_po_max=max(Y_po) + for i in Y_pl: + x1=np.random.ranf()*(Y_pl_max) #If cdf gave a max value of 1 + x2=X_pl[Y_pl.searchsorted(x1)-1] # we'd just call np.random.ranf() here + if x2>0: + Y_pl_skew+=[x2] + for i in Y_po: + x1=np.random.ranf()*(Y_po_max) #If cdf gave a max value of 1 + x2=X_po[Y_po.searchsorted(x1)-1] #we'd just call np.random.ranf() here + if x2>0: + Y_po_skew+=[x2] + + #Form cutoff power law distribution + pl_cutoff=cutoff_fn('pl') + po_cutoff=cutoff_fn('po') + + #2. Assign proboscis/stem lengths + po_lengths=np.arange(0,pollinators) + pl_lengths=np.arange(0,plants) + + def assign_lengths_pl(x): + index=np.random.randint(0,len(Y_pl_skew)) + return Y_pl_skew[index] + + def assign_lengths_po(x): + index=np.random.randint(0,len(Y_po_skew)) + return Y_po_skew[index] + + po_lengths=list(map(assign_lengths_po,po_lengths)) + pl_lengths=list(map(assign_lengths_pl,pl_lengths)) + + #3. Assign edges b/w pollinators and plants according to cut off exponential distribution. + def obtain_sequence(n,cutoff): + #Obtains degree distributions from the cutoff powerlaw distribution + out=[] + for i in range(0,n): + loc=np.random.randint(0,np.size(cutoff)) + out+=[int(math.ceil(cutoff[loc]))] #ceiling function rounds up values to nearest int (avoids 0 assigned degree) + return out + + #Obtain degree sequences from exponential cut off degree distribution + aseq=obtain_sequence(pollinators,po_cutoff) + bseq=obtain_sequence(plants,pl_cutoff) + #Their sums must be equal in order to form graph + while sum(aseq)!=sum(bseq): + aseq=obtain_sequence(pollinators,po_cutoff) + bseq=obtain_sequence(plants,pl_cutoff) + #Form bipartite graph (networkx documentation claims using Graph() may lead to inexact degree sequences) + G=bp.configuration_model(aseq, bseq,create_using=nx.Graph()) + + #Relabel the nodes... + names={} + for i in range(0,pollinators): + names[i]='po_'+str(i) + for i in range(0,plants): + names[i+pollinators]='pl_'+str(i) + + G=nx.relabel_nodes(G,names) + G=G.to_directed() #now every edge a-b is duplicated & directed (a->b and b->a) + + #4. Assign edge weights according to fitting function + #First, assign nodes their type & lengths from po_lengths and pl_lengths + #(better storage of values now that the network is actually formed) + counter=0 + for i in po_lengths: + G.nodes['po_'+str(counter)]['type'] = 'pollinator' + G.nodes['po_'+str(counter)]['length'] = i + counter+=1 + counter=0 + for i in pl_lengths: + G.nodes['pl_'+str(counter)]['type'] = 'plant' + G.nodes['pl_'+str(counter)]['length'] = i + counter+=1 + + #Then assign edge weight according to the lengths + for i in G.edges(): + a,b=G.nodes[i[0]]['length'],G.nodes[i[1]]['length'] #prob/stem lengths + diff=100.*abs(a-b)/((a+b)/2.0) + #If they're w/in the specified range, both edges (a->b and b->a) are assigned 2 + if diff<=10.0: + G[i[0]][i[1]]['data']=1. + G[i[1]][i[0]]['data']=1. + elif alength2 + G[i[1]][i[0]]['data']=1. + G[i[0]][i[1]]['data']=-1. + + #Finally, output the graph + return G + +def force_state(Graph,state_string,mod=False): + ''' + This function forces the network into the state specified by a binary + string (e.g. '001110101'), where the first position applies to pl_0 and + the last to po_m (m=# of pollinators). + + It returns the modified graph with node property 'state' assigned + according to state_string. + + It is called internally in other functions (e.g. steady_states, + active_edges). + + the "mod" parameter was added to accomodate transient_analysis, when we + have removed nodes in a nonuniform manner. leaving mod=False is faster, + but mod=True allows the function to work appropriately in this case. + ''' + if mod==False: + #We iterate first over all plant nodes + sc=0 #state counter + nc=0 #node counter + while 1==1: + node='pl_'+str(nc) + if Graph.has_node(node): + Graph.nodes[node]['state']=state_string[sc] + nc+=1 + sc+=1 + else: + break + + #And then over pollinator nodes + nc=0 + while 1==1: + node='po_'+str(nc) + if Graph.has_node(node): + Graph.nodes[node]['state']=state_string[sc] + nc+=1 + sc+=1 + else: + break + + else: + def mapper(x): + try: + return int(x[-2:]) #for e.g. pl_15 + except ValueError: + return int(x[-1]) #for e.g. pl_5 + #Determine maximum label value + labels=Graph.nodes() + labels=list(map(mapper,labels)) + nc=max(labels)+1 + sc=0 + for i in range(nc): + node='pl_'+str(i) + if Graph.has_node(node): + Graph.nodes[node]['state']=state_string[sc] + sc+=1 + + #And then over pollinator nodes + for i in range(nc): + node='po_'+str(i) + if Graph.has_node(node): + Graph.nodes[node]['state']=state_string[sc] + sc+=1 + + #After we've iterated over all nodes, we also should have iterated + #over all entries in the state string + if len(state_string)!=sc: + print(sc, len(state_string),len(Graph.nodes())) + print(state_string) + print('error in assigning node states in function \'force_state\'!') + + #We also want a dictionary of the states + sys_state={} + for i in Graph.nodes(): + sys_state[i]={} + sys_state[i][0]=int(Graph.nodes[i]['state']) + + return Graph, sys_state + +def update_network(Graph,sys_state,threshold=1,posweight=1,mode='sync',negweight=-1,negspecial=-2,posspecial=4,negnodes=[],extra_data=False): + ''' + This function performs one timestep on the network. For an edge to be active + after the update, the total incoming edge weights from active nodes must be + >= threshold. "posweight" allows for the weight of positive edges to be + effectively increased, which allows for another variation on update rules + (instead of just varying the threshold). + + mode is 'sync' by default, which updates all the nodes. The only other mode + currently supported is 'omit1', which randomly selects 1 node to not update. + + negweight is the weight of a negative edges (recent addition). + + negnodes are nodes which have their weights modified to values indicated + by 'posspecial' and 'negspecial.' (negnodes is a deprecated name; it + originally only handled modifications to negative edges.) + + It returns the updated network and an updated dictionary of states. + + This function should not normally need to be called by the end user. It is + called internally in other functions, which repeatedly update the network + under some conditions. + + if extra_data==True, the fn performs analysis of spectral nestedness, + connectance, and size at every time step. + ''' + def edge_filter(x): + #used for filtering list of edges to give only edges ending at target + return x[1]==i + def state_filter(x): + #used for filtering list of edges to give only edges with active source node + return sys_state[x[0]][counter-1]==1 or sys_state[x[0]][counter-1]=='1' + def weight_map(x): + #maps list of edges into list of their weights (-1/1), modified to posweight, negweight, and negspecial + if x[0] in negnodes and x[2]['data']==-1: + return negspecial #negative edges originating from negnode list are changed from negweight to negspecial + if x[0] in negnodes and x[2]['data']==1: + return posspecial #positive edges originated from negnode list are changed from posweight to posspecial + if x[2]['data']==-1: + return negweight #negative edges are changed from -1 to negweight + return posweight #positives are increased from 1 to posweight + + counter=list(sys_state.values()) + counter=len(counter[0].values()) #determines which round of updating we're on + + edgelist=Graph.edges(data=True) #stores all edges + if mode=='omit1': + nodelist=Graph.nodes() #if we are to omit updating 1 node, make a list of the nodes + omit_node=np.random.randint(len(nodelist)) #choose 1 index to omit + for i in Graph.nodes(): + if mode=='omit1': + if nodelist.index(i)==omit_node: + sys_state[i][counter]=Graph.nodes[i]['state'] #don't update graph node dictionary, but make sure we update the state dictionary with the prev. state + continue #don't update this node if we're in omit1 mode + edges=list(filter(edge_filter,edgelist)) #edges is list of edges ending at i + edges=list(filter(state_filter,edges)) #edges is list of edges with active source (and ending at i) + edges=list(map(weight_map,edges)) #edges is list of weights of edges with active source nodes that end at i + weight=sum(edges) + if weight never). Setting equal to 'inf' makes the changes + permanent throughout the simulation. + + Note also that FORCE_EXTINCT can refer to forced introductions, as well. + + In this project we additionally output summary statistics of every network + configuration during the dynamics. + + ---- + This function takes an initialized graph and its dictionary of states and + updates it 'in stasis' (w/o externally introduced species) until it reaches + a steady state, limit cycle, or 100 iterations have occurred. + + threshold is what the total edge weight of incoming edges with active source + nodes must meet or exceed for a node to remain active after that update + round. + + posweight is the weight of a single positive edge. + + ind is an int that corresponds to the position in the state string that will + be forced to remain ON or OFF (based on TRANS_TYPE; 0 for 0->1 perturbations + and 1 for 1->0 perturbations.) ind is False if no such forcing is to occur. + + It returns which condition is met and the sys_state dictionary. + ''' + SN_out,C_out,size_out,adj_out=[],[],[],[] + if steplimit>2**20: + steplimit=2**20 + for i in range(steplimit): #Counting initialization, yields total of steplimit+1 states + #force a particular species to be ON or OFF prior to first update if we're forcing a particular type of perturbation + if i==0 and (ind!=False or type(ind)==int): #(0==False returns True, so account for that) + if TRANS_TYPE==1: #what species will be perturbed to (opposite of what they are perturbed from, i.e. TRANS_TYPE) + target=0 + else: + target=1 + nodes=Graph.nodes() + nodes=node_sorter(nodes) + end=sys_state[nodes[ind]].keys() + sys_state[nodes[ind]][end[-1]]=target + #perform update once we've made the perturbation + if extra_data: + Graph,sys_state,stats=update_network(Graph=Graph,sys_state=sys_state,threshold=threshold,posweight=posweight,mode=mode,negweight=negweight,negspecial=negspecial,posspecial=posspecial,negnodes=negnodes,extra_data=True) + SN_out+=[stats[0]] + C_out+=[stats[1]] + size_out+=[stats[2]] + adj_out+=[stats[3]] + else: + Graph,sys_state=update_network(Graph=Graph,sys_state=sys_state,threshold=threshold,posweight=posweight,mode=mode,negweight=negweight,negspecial=negspecial,posspecial=posspecial,negnodes=negnodes) + #force the specified node to remain ON/OFF if FORCE_EXTINCT==True + if (ind!=False or type(ind)==int) and FORCE_EXTINCT==True: + if TRANS_LIMIT=='inf' or i1 perturbations + and 1 for 1->0 perturbations.) ind is False if no such forcing is to occur. + + It returns which condition is met and the sys_state dictionary. + ''' + if TRANS_TYPE==1: + target=0 + else: + target=1 + SN_out,C_out,size_out,adj_out=[],[],[],[] + if steplimit>2**20: + steplimit=2**20 + for i in range(steplimit): #Counting initialization, yields total of steplimit+1 states + #force a particular species to be ON or OFF prior to first update if we're forcing a particular type of perturbation + if i==0 and (ind!=False or type(ind)==int): #(0==False returns True, so account for that) + if TRANS_TYPE==1: #what species will be perturbed to (opposite of what they are perturbed from, i.e. TRANS_TYPE) + target=0 + else: + target=1 + nodes=Graph.nodes() + nodes=node_sorter(nodes) + for j in ind: + end=sorted(sys_state[nodes[j]]) # (edited for Python 3) sorted list of dict keys: 0,1,2,... to current time step counter value + sys_state[nodes[j]][end[-1]]=target + #perform update + if extra_data: + Graph,sys_state,stats=update_network(Graph=Graph,sys_state=sys_state,threshold=threshold,posweight=posweight,mode=mode,negweight=negweight,negspecial=negspecial,posspecial=posspecial,negnodes=negnodes,extra_data=True) + SN_out+=[stats[0]] + C_out+=[stats[1]] + size_out+=[stats[2]] + adj_out+=[stats[3]] + else: + Graph,sys_state=update_network(Graph=Graph,sys_state=sys_state,threshold=threshold,posweight=posweight,mode=mode,negweight=negweight,negspecial=negspecial,posspecial=posspecial,negnodes=negnodes) + #force the specified node to remain ON/OFF if FORCE_EXTINCT==True + if ind!=False and FORCE_EXTINCT==True: + if TRANS_LIMIT=='inf' or i<=TRANS_LIMIT: #Here, only perform forcing if we're at a tolerable time step (always if TRANS_LIMIT = 0 , otherwise only for specified #) + for val in ind: + end=sorted(sys_state[nodes[val]]) + sys_state[nodes[val]][end[-1]]=target + + cond=check_steady(Graph=Graph,sys_state=sys_state,mode=mode) + if cond!='chaotic': #Stop evolving the system if we've reached a steady state or limit cycle + break + if extra_data: + return cond, sys_state, SN_out, C_out, size_out,adj_out + else: + return cond, sys_state + +def node_sorter(x): + ''' + Since our nodes are strings, just calling nodes.sort() will put e.g. 'pl_5' + after 'pl_10' - this function sorts them properly and returns them in a + list. + ''' + + def node_plucker(y): + #filters out nodes not of appropriate length + return len(y)==mn + + pl=[] + po=[] + for i in x: #split po and pl into two lists + if i.count('l')!=0: + pl+=[i] + else: + po+=[i] + + sorter=[] + while len(pl)>0: + mn=min(list(map(len,pl))) #minimum length (i.e. 1 digit, 2 digit,...) + small=list(filter(node_plucker,pl)) #small now has the nodes of the minimum length + pl=list(set(pl)-set(small)) #take all of 'small' out of 'pl' + small.sort() #order the nodes in small + sorter+=small #add these to sorter + + while len(po)>0: + mn=min(list(map(len,po))) #repeat process for po + small=list(filter(node_plucker,po)) #should really rewrite this in more condensed fashion + po=list(set(po)-set(small)) + small.sort() + sorter+=small + + return sorter + + + +def active_edges(Graph, states): + ''' + This function takes as input a graph and list of state strings (e.g. as in + the output of the function steady_states; in form ['00111','11000',...]), + and returns a list of the edges that are active (i.e. source and target + nodes are ON) and their edge weights. + ''' + + def state_filter(x): + #used for filtering list of edges to give only edges with active source and sink nodes + return Graph.nodes[x[0]]['state']=='1' and Graph.nodes[x[1]]['state']=='1' + + edgelist=Graph.edges(data=True) + for state in states: + print('Considering state vector:',state) + print('The active edges are:') + G,gbg=force_state(Graph,state) #here we don't care about the returned state dictionary + edges=list(filter(state_filter,edgelist)) #edges holds only edges with active source and target node + + for i in edges: + print(i[0],'->',i[1],'\t weight=',i[2]['data']) + +def perturb_state(Graph,states,steplimit=True,threshold=1, posweight=1, report=False,negweight=-1,negspecial=-2,negnodes=[]): + ''' + This function takes as input a graph and list of state strings (e.g. as in + the output of the function steady_states; in form ['00111','11000',...]), + and for each node in each state: + 1. Turns it from ON to OFF or vice versa + 2. Allows the dynamics to evolve until either... + - the state goes to a steady state, OR + - the state reaches a limit cycle, OR + - the state evolves the maximum number of times (i.e. goes to a complex + limit cycle). + A dictionary of the form + {state:[(perturbation1,result,stepcount),(perturbation2,result,stepcount)...],...} + is returned, where pertubations are each of the possible OFF-ON changes to + each string in 'states', and the results are the strings of the steady + states, or an int of the limit cycle length, or 'chaotic' otherwise. + stepcount is the number of steps it takes to reach the limit cycle or + steady state. + + NOTE: steplimit defaults to the maximum possible (2^N), but may be + overridden with any positive integer value (as a float or int). + ''' + def state_modifier(state,i): + #takes a state string and forces the value at index i to be 1 + out=state[0:i]+'1'+state[i+1:] + return out + + if report == True: + print('Executing perturb_state...') + out={} #keys will be entries in 'states'. Each will have 1 value; a list of (modified state,condition) tuples + N=len(states[0]) #len of state vector (i.e. # of nodes) + if steplimit==True: + steplimit=2**N #The maximum number of steps we can take on the transition graph for each considered steady state + count=0 #for printing status updates + for state in states: #iterate over all input states + if report==True: + print(100.0*count/len(states),'% complete') + out_list=[] #Will be stored as a value for this state in 'out' dictionary + for i in range(N): #iterate over all nodes in this state + if state[i]=='0': #If this node is OFF... + new_state=state[0:i]+'1'+state[i+1:] #Make a copy where it is ON + else: #Otherwise it is ON... + new_state=state[0:i]+'0'+state[i+1:] #So we make a copy where it is OFF + G, sys_state=force_state(Graph,new_state) #Initialize the graph to this modified state + #Allow the system to advance in time; cond is int of limit + #cycle length, 'chaotic', or the steady state string vector e.g. '01101..' + cond,sys_state=check_sufficiency(G,sys_state,int(steplimit),threshold=threshold, posweight=posweight,negweight=negweight,negspecial=negspecial,negnodes=negnodes) + state_count=len(sys_state['pl_0'].keys()) #The number of updates the network went through + if type(cond)==list: #If we enter a limit cycle we report the # of steps it takes to get TO the limit cycle + state_count-=(1+len(cond)) #(subtract additional 1 b/c #edges=#states-1) + elif type(cond)==str: #If we hit a steady state, we report the # of steps it takes to reach the SS + state_count-=2 #subtract 1 b/c #edges=#states-1; subtract additional 1 b/c the state dict has the steady state 2x + out_list+=[(new_state,cond,state_count)] #Store these results + out[state]=out_list #Store final set of perturbed states and the result of their evolution + count+=1 + if report==True: + print('Perturb_state finished.') + + if report==True: #If desired, we can output the results in an easy to read format + print('--------------') + for i in out.keys(): + print('Reporting perturbations on state %s'%i) + n=len(out[i]) + for j in range(n): + #Different formatting for chaotic, limit cycle, or steady state (if steplimit=True, chaotic will never be called) + if type(out[i][j][1])==list: + print('Perturbed state %s \t takes \t %d \t steps to go to a cycle of length: %s'%(out[i][j][0],out[i][j][2],out[i][j][1])) + elif out[i][j][1]=='chaotic': + print('Perturbed state %s \t takes \t %d \t steps and fails to reach a limit cycle or steady state.'%(steplimit,out[i][j][0])) + else: + print('Perturbed state %s \t takes \t %d \t steps to go to a steady state: %s'%(out[i][j][0],out[i][j][2],out[i][j][1])) + print('--------------') + + + return out + +def lc_comparison(LC): + ''' + This function compares all the states in a given limit cycle, and determines + how many of the nodes share the same state: x=(# of agreement/# of bits) + The average of x over all pairs of states in the limit cycle is calculated. + + LC is a LIST of state STRINGS in a + particular limit cycle. As such, 'x' is calculated for each entry in LC; + the average value is returned (so, "what is the average homogeneity of the + states in a limit cycle in this network?"). + ''' + + vals=[] + n=len(LC) + m=len(LC[0]) #The # of nodes in each state + + for i in range(n): + for j in range(i+1,n): #Compare every pair of states (order doesn't matter) + s1=LC[i] + s2=LC[j] + counter=0 + for k in range(m): #For every bit position... + if s1[k]==s2[k]: + counter+=1 #...determine if the bits from each state agree. + vals+=[1.0*counter/m] #Store the rating for each pair... + + out=np.mean(vals) + return out #...and return the average result. + + + +def lc_lc_comparison(lc1,lc2): + ''' + This function takes two limit cycles and determines an average similarity + between them. The process is to determine the average expression for each + species among a given limit cycle's states, and then make a comparison among + the LCs by abs(LC1_exp-LC2_exp). + + We normalize the results such that an end value of 0 corresponds to exactly + equal expression levels for all species, and 1 corresponds to exactly + opposite expression levels for all species. + ''' + + m=len(lc1[0]) #the number of species + n1,n2=len(lc1),len(lc2) + + avg1,avg2=[0]*m,[0]*m #holds the average expression (0 vs 1) among the LC states for each input LC + + #determine average expression levels in the LCs + for i in range(n1): + for j in range(i+1,n1): + s1,s2=lc1[i],lc1[j] + for k in range(m): + if s1[k]==s2[k]: + avg1[k]+=1.0/(math.factorial(n1)/(2*math.factorial(n1-2))) #each term is one of the possible state comparison sets (n1 choose 2) + for i in range(n2): + for j in range(i+1,n2): + s1,s2=lc2[i],lc2[j] + for k in range(m): + if s1[k]==s2[k]: + avg2[k]+=1.0/(math.factorial(n1)/(2*math.factorial(n1-2))) + + #now compare those averages + out=0 + for i in range(m): + out+=(1.0/m)*abs(avg1[i]-avg2[i]) + + return out + +def ss_lc_comparison(SS,LC): + ''' + This function takes a steady state and limit cycle and returns an average + similarity between them. An average expression level for the lc is + determined (as in lc_lc_comparision(), and compared to the ss expression + levels. + ''' + n=len(LC) + m=len(LC[0]) #The # of nodes in each state + avg=[0]*m + + for i in range(n): + for j in range(i+1,n): #Compare every pair of states (order doesn't matter) + s1=LC[i] + s2=LC[j] + for k in range(m): + if s1[k]==s2[k]: + avg[k]+=1.0/(math.factorial(n)/(2*math.factorial(n-2))) + + #compare this average LC expression level to the SS... + out=0 + for i in range(m): + out+=(1.0/m)*abs(avg[i]-int(SS[i])) + + return out + +def ss_ss_comparison(SS1,SS2): + ''' + This function takes two steady states and performs a bitwise comparison + of their expression values. The average difference is returned, so a value + of 0 corresponds to exact agreement (same SS) and a value of 1 correponds + to exact disagreement (e.g. 000 and 111). + ''' + n=len(SS1) + out=0 + for i in range(n): + out+=(1.0/n)*abs(int(SS1[i])-int(SS2[i])) + + return out + +def steady_states(Graph, threshold=1,report=False): + ''' + This function searches through all the states in a network, and determines + how many are steady states. + + It outputs as text the # of steady states, and if report=True, it prints + the states in binary according to the format: + pl_0,pl_1,...pl_n, po_0,...po_m + (where there are n, m plants, pollinators). + + NOTE: It generates an error for very large networks. This function is + not currently used in "official" analysis (i.e. in random_graph_analysis.py). + ''' + + def update(): + ''' + This function begins updating the states of the network. If any node's + new state is different from its last, the procedure stops and 'False' is + returned. + + If each node keeps its old state, 'True' is returned. + ''' + + def edge_filter(x): + #used for filtering list of edges to give only edges ending at target + return x[1]==i + def state_filter(x): + #used for filtering list of edges to give only edges with active source node + return Graph.nodes[x[0]]['state']=='1' + def weight_map(x): + #maps list of edges into list of their weights + return x[2]['data'] + + out=True #returned variable. Becomes False if a state is updated to a new value (ON->OFF or OFF->ON) + + edgelist=Graph.edges(data=True) #stores all edges + for i in Graph.nodes(): + edges=list(filter(edge_filter,edgelist)) #edges is list of edges ending at i + edges=list(filter(state_filter,edges)) #edges is list of edges with active source (and ending at i) + edges=list(map(weight_map,edges)) #edges is list of weights of edges with active source nodes that end at i + weight=sum(edges) + if weightb->c->a same as b->c->a->b + # do this by comparing to a separate list of LCs as sets + # but also store one copy of LC as usual list of strings for further analysis + x = set(stable_state_A[0]) + if x not in LC_set_list_A: + LC_set_list_A.append(x) + LC_list_A.append([stable_state_A[0],SN_A,C_A,size_A]) + elif type(stable_state_A[0])==str and stable_state_A[0] not in SS_list_A and stable_state_A[0].count('1') > 0: + # capture steady states here + # other possibility (outside of this if/elif block) is a repeat attractor + # or 'chaotic' if none is found -- very unlikely + # or 'all dead' -- possible but not of interest here + SS_list_A.append([stable_state_A[0],SN_A,C_A,size_A]) + + # Repeat For B... + # If this is an attractor we haven't seen yet, add it to the appropriate list + if type(stable_state_B[0]) == list: + # we have a LC; need to count a->b->c->a same as b->c->a->b + # do this by comparing to a separate list of LCs as sets + # but also store one copy of LC as usual list of strings for further analysis + x = set(stable_state_B[0]) + if x not in LC_set_list_B: + LC_set_list_B.append(x) + LC_list_B.append([stable_state_B[0],SN_B,C_B,size_B]) + elif type(stable_state_B[0])==str and stable_state_B[0] not in SS_list_B and stable_state_B[0].count('1') > 0: + # capture steady states here + # other possibility (outside of this if/elif block) is a repeat attractor + # or 'chaotic' if none is found -- very unlikely + # or 'all dead' -- possible but not of interest here + SS_list_B.append([stable_state_B[0],SN_B,C_B,size_B]) + + + # Now start the network at a state corresponding to the combination of + # two of the above attractors (A into B, then B into A) and allow + # dynamics to proceed; record final state as in find_attractors.py --- + + counter = 0 + counter_max = (len(SS_list_A) + len(LC_list_A))*(len(SS_list_B) + len(LC_list_B)) # nested for loop + + counter_p = 0 + for i,j in product(SS_list_A + LC_list_A,SS_list_B+LC_list_B): + # periodically report on progress + if 100*counter/counter_max > counter_p: + TIME1 = time.time() + print('{0:.2f}% through combination analysis (total = {1:.2g}) at {2:.2f} minutes elapsed.'.format(counter_p,counter_max,(TIME1-TIME0)/60.0)) + counter_p += 10 + + inner_data_AB = run_invasion(i[0],j[0],host='B') # A attractor i invades B network in state j + inner_data_BA = run_invasion(j[0],i[0],host='A') # B attractor j invades A network in state i + + counter += 1 + + #Finally, append this data to the master dictionary + out_data[10*PBS_ITERATOR+graph_index][counter]=[inner_data_AB, inner_data_BA,{'SN':[i[1],j[1]],'C':[i[2],j[2]],'size':[i[2],j[2]]}] + +# Write to file ---------------------------------------------------------------- +f = open(r'attractors_comb_random_'+str(PBS_ITERATOR)+'.data','wb') +#f = open(r'/to/path/attractors_comb_random_'+str(PBS_ITERATOR)+'.data','wb') +#pickle.dump(out_data,f) +#f.close() + +TIME1 = time.time() +print('Script completed in {0:.2f} minutes.'.format(((TIME1-TIME0)/60))) diff --git a/invasions_whole.py b/invasions_whole.py new file mode 100644 index 0000000..e1480ea --- /dev/null +++ b/invasions_whole.py @@ -0,0 +1,234 @@ +""" +Once we have determined the attractors of the "half interaction graphs", we +here start the system in a state that corresponds to an attractor from half A +invading an attractor from half B (and vice versa -- exhaustive combinations). + +We allow the system to re-stabilize and capture the resulting attractor along +with standard network properties (e.g. nestedness). + +Master workflow: + 1. What are the species interaction graphs? (create_graphs.py) + 2. What are the stable communities in these graphs, if we segregate them + in to two "half communities" to enforce isolated communities that + nonetheless have inter-community interactions for when invasions occur? + (find_attractors.py) + 3. What happens if two stable communities are combined and allowed to + restabilize? (THIS FILE and invasions_random.py) + +The file is structured to run either on a local machine or submitted to a +High-Performance Computing center. See comments. +(This research was performed via the Ohio Supercomputer Center, osc.edu) + +Python 3.x + +Colin Campbell +campbeco@mountunion.edu +""" + +import pickle +import glob +#import numpy as np +import function_library as nv +from itertools import product +import time +import random +import networkx as nx +import os +import sys + +# Function definitions -------------------------------------------------------- +def run_invasion(invader_state,host_state,host): + ''' + Considers two "half interaction graphs": an invading community + ('invader_state') from one half invades the other half, which is in state + 'host_state' to start. Parameter 'host' indicates which half (A=0,B=1) + is the host community. + + The system is started from invader_state + host_state and all species + in the host community + invading species are allowed to interact in + subsequent dynamics. The stable state is identified and key information + is stored and returned in a dict. + ''' + + #inner dictionary that will hold the data for this iteration + inner_data={'start_state':[],'stable_state':[],'sys_state':[], + 'SN':[],'C':[],'size':[],'benchmark':0,'adj':[],'i':[], + 'j':[],'start_i':[],'start_j':[]} + + + G = graphs[start].copy() + + #If either state is a LC, choose a constituent state at random + if type(invader_state) == str: start_inv = invader_state[:] + else: + invader_state_copy = invader_state.copy() + random.shuffle(invader_state_copy) + start_inv = invader_state_copy.pop() + + if type(host_state) == str: start_host = host_state[:] + else: + host_state_copy = host_state.copy() + random.shuffle(host_state_copy) + start_host = host_state_copy.pop() + + # now combine these starting states + start_state = ''.join(['0' if start_inv[x] == '0' and start_host[x] == '0' else '1' for x in range(N)]) + + # Whichever half is hosting, the species in the -other- half are by default + # excluded from participating in the dynamics... + if host == 'A': + off_list = list(range(25,50))+list(range(125,200)) + elif host == 'B': + off_list = list(range(25))+list(range(50,125)) + + #... with the exception of those species that are invading from the other + # half. + for position,val in enumerate(start_state): + if val == '1' and position in off_list: + off_list.remove(position) + + # execute trial and store data + G,sys_state = nv.force_state(G,start_state) + stable_state,sys_state,SN_vals,C_vals,size_vals,adj1 = nv.check_sufficiency_mod(G,sys_state,posweight=4,ind=off_list,TRANS_TYPE=1, FORCE_EXTINCT=True,extra_data=True) #determine the resulting stable configuration + + inner_data['start_state'].append(start_state) + inner_data['stable_state'].append(stable_state) + inner_data['sys_state'].append(sys_state) + inner_data['SN'].append(SN_vals) + inner_data['C'].append(C_vals) + inner_data['size'].append(size_vals) + inner_data['benchmark']=len(size_vals) + + return inner_data + + +# End function definitions ---------------------------------------------------- + + +TIME0 = time.time() + +# Which set of interaction graphs are we considering? ------------------------- +#PBS_ITERATOR=int(os.environ["PBS_ARRAYID"]) #We use this variable to slice up the job over 100 jobs on OSC. Use: qsub -t 0-99 pbs_script.sh +PBS_ITERATOR=10 #Or dummy variable for local debugging + +# Optionally do not overwrite existing output files +# if os.path.exists('/to/path/attractors_comb_whole_'+str(PBS_ITERATOR)+'.data'): +# sys.exit("Output file already exists.") + +# Load graphs from memory ----------------------------------------------------- +#Local +f = open(r'/graphlist.data','rb') +# OSC +#f = open(r'/to/path/graphlist.data','rb') +graphs = pickle.load(f) +f.close() + +graphs=graphs[10*PBS_ITERATOR:(10*PBS_ITERATOR)+10] +N = nx.number_of_nodes(graphs[0]) + +# Local vs OSC two lines below +for fname in glob.glob(r'attractors_'+str(PBS_ITERATOR)+'.data'): +# for fname in glob.glob(r'/to/path/attractors_'+str(PBS_ITERATOR)+'.data'): + # Open the file for this batch of 10 interaction networks ---------------- + f = open(fname,'rb') + data = pickle.load(f) + f.close() + + ''' + structure of 'data' is as follows: data[graph_index][start_index]=[inner_data for A, inner_data for B] + where: + graph_index is from 0 to 999 (which interaction graph?) + start_index is from 0 to 99 (which 'choose a starting state' trial) + and 'inner_data' is the inner dictionary for community half A and B; + they notably have keys 'stable state' and 'size'. + ''' + + out_data = {} + + + for start,graph_index in enumerate(data.keys()): + print('Starting graph...') + out_data[10*PBS_ITERATOR+graph_index]={} + + LC_list_A, LC_list_B = [], [] + LC_set_list_A, LC_set_list_B = [], [] + SS_list_A, SS_list_B = [], [] + + # First just identify the unique attractors for half A, B ------------- + for start_index in data[graph_index].keys(): + stable_state_A = data[graph_index][start_index][0]['stable_state'] + size_A = data[graph_index][start_index][0]['size'] + SN_A = data[graph_index][start_index][0]['SN'][0][-1] + C_A = data[graph_index][start_index][0]['C'][0][-1] + + stable_state_B = data[graph_index][start_index][1]['stable_state'] + size_B = data[graph_index][start_index][1]['size'] + SN_B = data[graph_index][start_index][1]['SN'][0][-1] + C_B = data[graph_index][start_index][1]['C'][0][-1] + + # For A... + # If this is an attractor we haven't seen yet, add it to the appropriate list + if type(stable_state_A[0]) == list: + # we have a LC; need to count a->b->c->a same as b->c->a->b + # do this by comparing to a separate list of LCs as sets + # but also store one copy of LC as usual list of strings for further analysis + x = set(stable_state_A[0]) + if x not in LC_set_list_A: + LC_set_list_A.append(x) + LC_list_A.append([stable_state_A[0],SN_A,C_A,size_A]) + elif type(stable_state_A[0])==str and stable_state_A[0] not in SS_list_A and stable_state_A[0].count('1') > 0: + # capture steady states here + # other possibility (outside of this if/elif block) is a repeat attractor + # or 'chaotic' if none is found -- very unlikely + # or 'all dead' -- possible but not of interest here + SS_list_A.append([stable_state_A[0],SN_A,C_A,size_A]) + + # Repeat For B... + # If this is an attractor we haven't seen yet, add it to the appropriate list + if type(stable_state_B[0]) == list: + # we have a LC; need to count a->b->c->a same as b->c->a->b + # do this by comparing to a separate list of LCs as sets + # but also store one copy of LC as usual list of strings for further analysis + x = set(stable_state_B[0]) + if x not in LC_set_list_B: + LC_set_list_B.append(x) + LC_list_B.append([stable_state_B[0],SN_B,C_B,size_B]) + elif type(stable_state_B[0])==str and stable_state_B[0] not in SS_list_B and stable_state_B[0].count('1') > 0: + # capture steady states here + # other possibility (outside of this if/elif block) is a repeat attractor + # or 'chaotic' if none is found -- very unlikely + # or 'all dead' -- possible but not of interest here + SS_list_B.append([stable_state_B[0],SN_B,C_B,size_B]) + + + # Now start the network at a state corresponding to the combination of + # two of the above attractors (A into B, then B into A) and allow + # dynamics to proceed; record final state as in find_attractors.py --- + + counter = 0 + counter_max = (len(SS_list_A) + len(LC_list_A))*(len(SS_list_B) + len(LC_list_B)) # nested for loop + + counter_p = 0 + for i,j in product(SS_list_A + LC_list_A,SS_list_B+LC_list_B): + # periodically report on progress + if 100*counter/counter_max > counter_p: + TIME1 = time.time() + print('{0:.2f}% through combination analysis (total = {1:.2g}) at {2:.2f} minutes elapsed.'.format(counter_p,counter_max,(TIME1-TIME0)/60.0)) + counter_p += 10 + + inner_data_AB = run_invasion(i[0],j[0],host='B') # A attractor i invades B network in state j + inner_data_BA = run_invasion(j[0],i[0],host='A') # B attractor j invades A network in state i + + counter += 1 + + #Finally, append this data to the master dictionary + out_data[10*PBS_ITERATOR+graph_index][counter]=[inner_data_AB, inner_data_BA,{'SN':[i[1],j[1]],'C':[i[2],j[2]],'size':[i[2],j[2]]}] + +# Write to file ---------------------------------------------------------------- +f = open(r'attractors_comb_whole_'+str(PBS_ITERATOR)+'.data','wb') +#f = open(r'/to/path/attractors_comb_whole_'+str(PBS_ITERATOR)+'.data','wb') +#pickle.dump(out_data,f) +#f.close() + +TIME1 = time.time() +print('Script completed in {0:.2f} minutes.'.format(((TIME1-TIME0)/60)))