<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Load Data
Step1:  Graphs of population distribution and trajectories for base parameter set
Step2:  Given a delta_f of .03 and a K of 1 million, we expect it would take about 450 generations for a new mutation increasing fitness by one step to fix after it invades the population
Step3:  Varying the probability of beneficial mutations
Step4:  Varying the probability of accurate mutations
Step5:  Varying $\Delta f$
Step6:  Varying K
The population size changed due to a bug in all the simulations where N = 10^10. Because of this, it is necessary to cut off the last set of simulations for K=10^10
<ASSISTANT_TASK:>
Python Code:
fbvary_results = []
for file in glob.glob('runs/evolved_mu_f_b_vary?replicate?datetime.datetime(2019, 5, *).hdf5'):
    try:
        fbvary_results.append(popev.PopulationReader(file))
    except OSError:
        pass
favary_results = []
for file in glob.glob('runs/evolved_mu_f_a_vary?replicate?datetime.datetime(2019, 5, *).hdf5'):
    try:
        favary_results.append(popev.PopulationReader(file))
    except OSError:
        pass
delta_f_vary_results = []
for file in glob.glob('runs/evolved_mu_delta_f_vary?replicate?datetime.datetime(2019, 5, *).hdf5'):
    try:
        delta_f_vary_results.append(popev.PopulationReader(file))
    except OSError:
        pass
    
K_vary_results = []
for file in glob.glob('runs/evolve_mu_K_vary?replicate?datetime.datetime(2019, 6, *).hdf5'):
    try:
        K_vary_results.append(popev.PopulationReader(file))
    except OSError:
        pass
len(fbvary_results[0])
cd Thesis_Data_and_Figures/
Explanation: Load Data
End of explanation
basepop = fbvary_results[50]
print('The parameters of the base population are delta_f: {},'
      ' M: {}, P_b: {}, P_a: {}, P_mu: {}, K: {}'.format(basepop.delta_fitness, basepop.mu_multiple,
                                                        basepop.fraction_beneficial, basepop.fraction_accurate,
                                                        basepop.fraction_mu2mu, basepop.pop_cap))
Explanation: Graphs of population distribution and trajectories for base parameter set
End of explanation
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')
logNbarplot(ax, basepop(15000))
plt.savefig('traveling_wave_full_distribution.pdf')
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
logNvsfplot(ax, basepop(15000))
plt.savefig('traveling_wave_f_distribution.pdf')
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
logNvsmuplot(ax, basepop(15000))
plt.savefig('traveling_wave_mu_distribution.pdf')
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
plt.plot(basepop.mode_fitness[:10000], marker='');
ax.set_xlabel('generation', fontsize=32);
ax.set_ylabel('mode of population fitness, $f_{mode}$', fontsize=28);
plt.savefig('traveling_wave_f_evolution.pdf')
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
plt.semilogy(basepop.mode_mutation_rate[:], marker='');
ax.minorticks_off()
ax.get_yaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
ax.set_yticks(np.unique(basepop.mode_mutation_rate[:]))
ax.set_xlabel('generation', fontsize=32);
ax.set_ylabel('mode of population\nmutation rate, $\mu_{mode}$', fontsize=28);
plt.savefig('traveling_wave_mu_evolution.pdf')
Explanation: Given a delta_f of .03 and a K of 1 million, we expect it would take about 450 generations for a new mutation increasing fitness by one step to fix after it invades the population
End of explanation
fbs = [pop.fraction_beneficial for pop in fbvary_results]
dfdts = [mean_df_dt(pop) for pop in fbvary_results]
fbs, mean_dfdts, std_dfdts = mean_and_std_by(fbs, dfdts)
fig = plt.figure(figsize=(10,10))
plt.errorbar(fbs, mean_dfdts, 2*std_dfdts, linestyle='');
plt.xscale('log')
plt.yscale('log')
plt.xlabel('probability of a beneficial mutation $P_{b}$', fontsize=28);
plt.ylabel(r'rate of fitness increase $\frac{df_{mean}}{dt}$', fontsize=28);
plt.savefig('df_dt_varying_P_b.pdf')
fbs = [pop.fraction_beneficial for pop in fbvary_results]
mu_avs = [np.mean(pop.mean_mutation_rate[:]) for pop in fbvary_results]
fbs, mean_mu_avs, std_mu_avs = mean_and_std_by(fbs, mu_avs)
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
ax.errorbar(fbs, mean_mu_avs, 2*std_mu_avs, linestyle='')
ax.set_xscale('log')
ax.set_yscale('log')
ax.minorticks_off()
ax.get_yaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
ax.set_yticks([.00125, .0025, .005, .01, .02, .04])
ax.set_xlabel('probability of a beneficial mutation $P_{b}$', fontsize=28);
ax.set_ylabel('mean mutation rate, $\mu$,\n(over population and time)', fontsize=28);
plt.savefig('mean_mu_varying_P_b.pdf')
Explanation: Varying the probability of beneficial mutations
End of explanation
fas = [pop.fraction_accurate for pop in favary_results]
dfdts = [mean_df_dt(pop) for pop in favary_results]
fas, mean_dfdts, std_dfdts = mean_and_std_by(fas, dfdts)
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
ax.errorbar(fas, mean_dfdts, 2*std_dfdts, linestyle='')
ax.set_xscale('log')
ax.set_xlabel('probability of an accuracy increasing mutation $P_{a}$', fontsize=28);
ax.set_ylabel(r'rate of fitness increase $\frac{df_{mean}}{dt}$', fontsize=28);
plt.savefig('df_dt_varying_P_a.pdf')
fas = [pop.fraction_accurate for pop in favary_results]
mu_avs = [np.mean(pop.mean_mutation_rate[:]) for pop in favary_results]
fas, mean_mu_avs, std_mu_avs = mean_and_std_by(fas, mu_avs)
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
ax.errorbar(fas, mean_mu_avs, 2*std_mu_avs, linestyle='')
ax.set_xscale('log')
ax.set_xlabel('probability of an accuracy increasing mutation $P_{a}$', fontsize=28);
ax.set_ylabel('mean mutation rate, $\mu$,\n(over population and time)', fontsize=28);
plt.savefig('mean_mu_varying_P_a.pdf')
Explanation: Varying the probability of accurate mutations
End of explanation
delta_fs = [pop.delta_fitness for pop in delta_f_vary_results]
dfdts = [mean_df_dt(pop) for pop in delta_f_vary_results]
delta_fs, mean_dfdts, std_dfdts = mean_and_std_by(delta_fs, dfdts)
fig = plt.figure(figsize=(10,10))
plt.errorbar(delta_fs, mean_dfdts, 2*std_dfdts, linestyle='', label='simulations')
def parabola(x, a):
    return a*x**2
popt, pcov = spopt.curve_fit(parabola,
                             np.array([pop.delta_fitness for pop in delta_f_vary_results]),
                             np.array([mean_df_dt(pop) for pop in delta_f_vary_results]))
effs = np.linspace(.005,.05)
plt.plot(effs, parabola(effs,*popt), 'r', marker='', label='parabolic fit')
plt.xlabel('size of changes in fitness, $\Delta f$', fontsize=28);
plt.ylabel(r'rate of fitness increase $\frac{df_{mean}}{dt}$', fontsize=28);
plt.legend()
plt.savefig('df_dt_varying_delta_f.pdf')
popt
delta_fs = [pop.delta_fitness for pop in delta_f_vary_results]
mu_avs = [np.mean(pop.mean_mutation_rate[:]) for pop in delta_f_vary_results]
delta_fs, mean_mu_avs, std_mu_avs = mean_and_std_by(delta_fs, mu_avs)
fig = plt.figure(figsize=(10,10))
plt.errorbar(delta_fs, mean_mu_avs, 2*std_mu_avs, linestyle='')
plt.xlabel('size of changes in fitness, $\Delta f$', fontsize=28);
plt.ylabel('mean mutation rate, $\mu$,\n(over population and time)', fontsize=28);
plt.savefig('mean_mu_varying_delta_f.pdf')
Explanation: Varying $\Delta f$
End of explanation
Ks = [pop.pop_cap for pop in K_vary_results]
dfdts = [mean_df_dt(pop) for pop in K_vary_results]
Ks, mean_dfdts, std_dfdts = mean_and_std_by(Ks, dfdts)
fig = plt.figure(figsize=(10,10))
plt.errorbar(Ks[:-1], mean_dfdts[:-1], 2*std_dfdts[:-1], linestyle='')
plt.xscale('log')
plt.xlabel('population size, $K$', fontsize=28);
plt.ylabel(r'rate of fitness increase $\frac{df_{mean}}{dt}$', fontsize=28);
plt.savefig('df_dt_varying_K.pdf')
Ks = [pop.pop_cap for pop in K_vary_results]
mu_avs = [np.mean(pop.mean_mutation_rate[:]) for pop in K_vary_results]
Ks, mean_mu_avs, std_mu_avs = mean_and_std_by(Ks, mu_avs)
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
ax.errorbar(Ks[:-1], mean_mu_avs[:-1], 2*std_mu_avs[:-1], linestyle='')
ax.set_xscale('log')
ax.set_xlabel('population size, $K$', fontsize=28);
ax.set_ylabel('mean mutation rate, $\mu$,\n(over population and time)', fontsize=28);
plt.savefig('mean_mu_varying_K.pdf')
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
ax.set_ylim(-.1,.06)
ax.set_yticks(np.linspace(-.1,.06,9))
ax.set_xlim(.01/4,.04)
ax.set_xscale('log')
ax.minorticks_off()
ax.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
ax.set_xticks([.01/4, .01/2, .01, .02, .04, .08])
ax.grid(True)
ax.text(.01, 0, '823', ha='center', va='center', fontsize=32, color='red')
ax.text(.01, -.02, '110', ha='center', va='center', fontsize=32, color='red')
ax.text(.02, 0, '58', ha='center', va='center', fontsize=32, color='red')
ax.text(.02, -.02, '23', ha='center', va='center', fontsize=32, color='red')
ax.set_xlabel('mutation rate $\mu$', fontsize=32)
ax.set_ylabel('fitness $f$', fontsize=32)
plt.savefig('example_fitness_and_mu_layout.pdf')
Explanation: Varying K
The population size changed due to a bug in all the simulations where N = 10^10. Because of this, it is necessary to cut off the last set of simulations for K=10^10
End of explanation
<END_TASK>