<jupyter_start><jupyter_text>Luria-Delbrück estimator corrected for death probability according to a continuous time model (Section 3.5 and section 3.5.3)This notebook is a Python implementation of the simulations described in section 3.5 and 3.5.3.It has been used _anytree_ library, from GitHub, to build the tree data structure that represent a generational tree.Repository: https://github.com/c0fec0de/anytreeDocumentation: https://anytree.readthedocs.io/<jupyter_code># import
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.stats import binom
from anytree import Node, RenderTree

# Check that used machine supports int64
# max int64 ~ 9*10^18
if(np.random.randint(0, high = 9223372036854775808)):
    print("Ok, this machine supports int64!")<jupyter_output>Ok, this machine supports int64!<jupyter_text>Implementation of the Luria-Delbrück stochastic experiment Import from *LD.py*.<jupyter_code>from LD import LD<jupyter_output><empty_output><jupyter_text>Number of cells with alive progeny varying extant cells' number and death probabilitySection 3.5, figure 10.<jupyter_code># !!! REALLY TIME CONSUMING if death probability goes to 0.5 !!!
# To reduce computational time (a lot) delete death_prob = 0.48 from death_probs

# Test for (extant, count_with_alive_progeny) == (extant, attempts)
# For different minimum sizes we plot 
# (extant, count_with_alive_progeny) varying death probability

# Note that extinct trees are discarded
print("Note that extinct trees are discarded!")

# Repeat simulations N times to get mean and devSt
N=1000

#variable parameters
death_probs = [0. , 0.1, 0.2, 0.3, 0.4, 0.48]
min_sizes = [32, 64, 128]

# results lists
# mean value
count_with_alive_progeny = []
# standard deviation
count_with_alive_progeny_std = []

# set seed to have reproducible results
np.random.seed(0)

# simulate
for min_size in min_sizes:
    print("Min size:", min_size)
    #results for this death probability
    this_min_size_results = []
    this_min_size_results_std = []
    
    for death_prob in death_probs:
        print("    death probability:", death_prob, end = " ")
        # sigle runs results
        count_alive = np.array([])
        
        for _ in range(N):
            if _%100 == 0:
                #show "progress bar"
                print(".", end=" ")
            if _ == N-1:
                print()
            tree = LD(bases = 1E+8, death_prob = death_prob, min_size = min_size)
            tree.count_with_alive_progeny()
            count_alive = np.append(count_alive, tree._count_with_alive_progeny)
            
        #save this death probability result
        this_min_size_results.append(np.mean(count_alive))
        this_min_size_results_std.append(np.std(count_alive))
        
        print("        Mean count_with_alive_progeny:",this_min_size_results[-1])
        print("        Std count_with_alive_progeny:",this_min_size_results_std[-1])
    
    count_with_alive_progeny.append(this_min_size_results)
    count_with_alive_progeny_std.append(this_min_size_results_std)
        
# plot results
fig = plt.figure(figsize=(10,7))
ax = fig.add_subplot(1, 1, 1)
#plt.title("Ruolo della ploidia\n"+str(N)+" simulazioni con soglia 1/32")
ax.set_xlabel('d/(b+d)', fontsize=20)
ax.set_ylabel('Cellule con progenie viva', fontsize=20)
ax.yaxis.offsetText.set_fontsize(20)

# set x y range
ax.set_xlim(-0.05, 0.55)
ax.set_ylim(0, 1300)

# set ticks
ax.tick_params(which='major', width=1.0, labelsize=20)
ax.tick_params(which='major', length=10, labelsize=20)

for i in range(len(min_sizes)):
    ax.errorbar(death_probs, count_with_alive_progeny[i], yerr=count_with_alive_progeny_std[i],
         marker ='o', capsize=6, capthick=3, lw=3,fmt=' ',markersize=10,
         label='min size '+str(min_sizes[i]))
    
ax.legend(fontsize=20, title='Legenda',title_fontsize=20)
plt.tight_layout()
ax.grid()

#plt.savefig("attempts_death_prob.pdf")
plt.show()

#print(count_with_alive_progeny)
#print(count_with_alive_progeny_std)<jupyter_output>Note that extinct trees are discarded!
Min size: 32
    death probability: 0.0 . . . . . . . . . . 
        Mean count_with_alive_progeny: 62.0
        Std count_with_alive_progeny: 0.0
    death probability: 0.1 . . . . . . . . . . 
        Mean count_with_alive_progeny: 93.84
        Std count_with_alive_progeny: 14.482278826206876
    death probability: 0.2 . . . . . . . . . . 
        Mean count_with_alive_progeny: 101.593
        Std count_with_alive_progeny: 14.970816644391848
    death probability: 0.3 . . . . . . . . . . 
        Mean count_with_alive_progeny: 115.889
        Std count_with_alive_progeny: 15.255644168634769
    death probability: 0.4 . . . . . . . . . . 
        Mean count_with_alive_progeny: 147.686
        Std count_with_alive_progeny: 21.115809338029173
    death probability: 0.48 . . . . . . . . . . 
        Mean count_with_alive_progeny: 230.377
        Std count_with_alive_progeny: 46.603271033265464
Min size: 64
    death probability: 0.0 . . . . . . . .[...]<jupyter_text>Test the analytic estimation of the number of cells with alive progeny given by the continuous time model Here trees are built given minimum final size, $N_\text{extant}$.Section 3.5.3, figure 11.<jupyter_code># Test for (extant, count_with_alive_progeny) = (extant, attempts)
# For different minimum sizes we plot 
# (mean generation number, count_with_alive_progeny) with fixed death probability

# Note that extinct trees are discarded
print("Note that extinct trees are discarded!")

#variable parameter
min_sizes = [32, 64, 128]

#fixed death probability
death_prob = 0.1
print("Fixed death probability:", death_prob)

# Repeat simulations N times to get mean and devSt
N=1000

# results lists
mean_count = []
count_std = []
mean_gen = []
gen_std = []

# set seed to have reproducible results
np.random.seed(0)
# simulate
for min_size in min_sizes:
       
    # sigle runs results
    count = np.array([])
    gen = np.array([])
        
    for _ in range(N):
        tree = LD(bases = 1E+8, death_prob = death_prob, min_size = min_size, accept_extinct=False)
        tree.count_with_alive_progeny()
        count = np.append(count, tree._count_with_alive_progeny)
        gen = np.append(gen, tree._num_layers)
    
    #save this min size result
    mean_count.append(np.mean(count))
    count_std.append(np.std(count))
    mean_gen.append(np.mean(gen))
    gen_std.append(np.std(gen))   
    
#analytic estimation
estimated = []
for t in mean_gen:
    attempts = 0.
    # mean number of cell with at least one alive daughter
    # from generation zero to the second-last one
    attempts = (1-math.pow(death_prob,2))*1. * (math.pow(2*(1-death_prob),t-1) -1) / (math.log(2*(1-death_prob)))
    # we add also the mean number of cells in the last layer 
    attempts += 1. * math.pow(2*(1-death_prob), t)
    estimated.append(attempts)
    

#plot results
fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('Numero di generazioni', fontsize=40)
ax.set_ylabel('Cellule con progenie viva', fontsize=40)
ax.yaxis.offsetText.set_fontsize(40)

### set x y range
#ax.set_xlim(5, 10)
#ax.set_ylim(50, 460)

## set ticks
ax.tick_params(which='major', width=2.0, labelsize=40)
ax.tick_params(which='major', length=20, labelsize=40)


for i in range(len(min_sizes)):
    ax.errorbar(mean_gen[i], mean_count[i], yerr=count_std[i], xerr=gen_std[i],
         marker ='o', capsize=12, capthick=6, lw=6,fmt=' ',markersize=20,
         label='min size '+str(min_sizes[i]))
    
ax.scatter(mean_gen, estimated, marker = 'x', s=800, label = 'Stima analitica', color = 'tab:red',zorder=10)

ax.legend(fontsize=40, title='Legenda',title_fontsize=40)
plt.tight_layout()
plt.grid()
plt.title("d/(b+d) = "+str(death_prob), fontsize=40)

#plt.savefig("attempts_size_integral_0.1.pdf")
plt.show()

print("Mean count:")
print(mean_count)
print("Std count:")
print(count_std)
print("Mean gen:")
print(mean_gen)
print("Std gen:")
print(gen_std)
print("Estimated:")
print(estimated)
# Test for (extant, count_with_alive_progeny) = (extant, attempts)
# For different minimum sizes we plot 
# (mean generation number, count_with_alive_progeny) with fixed death probability

# Note that extinct trees are discarded
print("Note that extinct trees are discarded!")

#variable parameter
min_sizes = [32, 64, 128]

#fixed death probability
death_prob = 0.2
print("Fixed death probability:", death_prob)

# Repeat simulations N times to get mean and devSt
N=1000

# results lists
mean_count = []
count_std = []
mean_gen = []
gen_std = []

# set seed to have reproducible results
np.random.seed(0)

# simulate
for min_size in min_sizes:
       
    # sigle runs results
    count = np.array([])
    gen = np.array([])
        
    for _ in range(N):
        tree = LD(bases = 1E+8, death_prob = death_prob, min_size = min_size, accept_extinct=False)
        tree.count_with_alive_progeny()
        count = np.append(count, tree._count_with_alive_progeny)
        gen = np.append(gen, tree._num_layers)
    
    #save this min size result
    mean_count.append(np.mean(count))
    count_std.append(np.std(count))
    mean_gen.append(np.mean(gen))
    gen_std.append(np.std(gen))   
    
#analytic estimation
estimated = []
for t in mean_gen:
    attempts = 0.
    # mean number of cell with at least one alive daughter
    # from generation zero to the second-last one
    attempts = (1-math.pow(death_prob,2))*1. * (math.pow(2*(1-death_prob),t-1) -1) / (math.log(2*(1-death_prob)))
    # we add also the mean number of cells in the last layer 
    attempts += 1. * math.pow(2*(1-death_prob), t)
    estimated.append(attempts)
    

#plot results
fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('Numero di generazioni', fontsize=40)
ax.set_ylabel('Cellule con progenie viva', fontsize=40)
ax.yaxis.offsetText.set_fontsize(40)

### set x y range
#ax.set_xlim(6, 12.5)
#ax.set_ylim(50, 500)

## set ticks
ax.tick_params(which='major', width=2.0, labelsize=40)
ax.tick_params(which='major', length=20, labelsize=40)


for i in range(len(min_sizes)):
    ax.errorbar(mean_gen[i], mean_count[i], yerr=count_std[i], xerr=gen_std[i],
         marker ='o', capsize=12, capthick=6, lw=6,fmt=' ',markersize=20,
         label='min size '+str(min_sizes[i]))
    
ax.scatter(mean_gen, estimated, marker = 'x', s=800, label = 'Stima analitica', color = 'tab:red',zorder=10)

ax.legend(fontsize=40, title='Legenda',title_fontsize=40)
plt.tight_layout()
plt.grid()
plt.title("d/(b+d) = "+str(death_prob), fontsize=40)

#plt.savefig("attempts_size_integral_0.2.pdf")
plt.show()

print("Mean count:")
print(mean_count)
print("Std count:")
print(count_std)
print("Mean gen:")
print(mean_gen)
print("Std gen:")
print(gen_std)
print("Estimated:")
print(estimated)
# Test for (extant, count_with_alive_progeny) = (extant, attempts)
# For different minimum sizes we plot 
# (mean generation number, count_with_alive_progeny) with fixed death probability

# Note that extinct trees are discarded
print("Note that extinct trees are discarded!")

#variable parameter
min_sizes = [32, 64, 128]

#fixed death probability
death_prob = 0.3
print("Fixed death probability:", death_prob)

# Repeat simulations N times to get mean and devSt
N=1000

# results lists
mean_count = []
count_std = []
mean_gen = []
gen_std = []

# set seed to have reproducible results
np.random.seed(0)

# simulate
for min_size in min_sizes:
       
    # sigle runs results
    count = np.array([])
    gen = np.array([])
        
    for _ in range(N):
        tree = LD(bases = 1E+8, death_prob = death_prob, min_size = min_size, accept_extinct=False)
        tree.count_with_alive_progeny()
        count = np.append(count, tree._count_with_alive_progeny)
        gen = np.append(gen, tree._num_layers)
    
    #save this min size result
    mean_count.append(np.mean(count))
    count_std.append(np.std(count))
    mean_gen.append(np.mean(gen))
    gen_std.append(np.std(gen))   
    
#analytic estimation
estimated = []
for t in mean_gen:
    attempts = 0.
    # mean number of cell with at least one alive daughter
    # from generation zero to the second-last one
    attempts = (1-math.pow(death_prob,2))*1. * (math.pow(2*(1-death_prob),t-1) -1) / (math.log(2*(1-death_prob)))
    # we add also the mean number of cells in the last layer 
    attempts += 1. * math.pow(2*(1-death_prob), t)
    estimated.append(attempts)
    

#plot results
fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('Numero di generazioni', fontsize=40)
ax.set_ylabel('Cellule con progenie viva', fontsize=40)
ax.yaxis.offsetText.set_fontsize(40)

### set x y range
#ax.set_xlim(8, 18)
#ax.set_ylim(80, 550)

## set ticks
ax.tick_params(which='major', width=2.0, labelsize=40)
ax.tick_params(which='major', length=20, labelsize=40)


for i in range(len(min_sizes)):
    ax.errorbar(mean_gen[i], mean_count[i], yerr=count_std[i], xerr=gen_std[i],
         marker ='o', capsize=12, capthick=6, lw=6,fmt=' ',markersize=20,
         label='min size '+str(min_sizes[i]))
    
ax.scatter(mean_gen, estimated, marker = 'x', s=800, label = 'Stima analitica', color = 'tab:red',zorder=10)

ax.legend(fontsize=40, title='Legenda',title_fontsize=40)
plt.tight_layout()
plt.grid()
plt.title("d/(b+d) = "+str(death_prob), fontsize=40)
#plt.savefig("attempts_size_integral_0.3.pdf")
plt.show()

print("Mean count:")
print(mean_count)
print("Std count:")
print(count_std)
print("Mean gen:")
print(mean_gen)
print("Std gen:")
print(gen_std)
print("Estimated:")
print(estimated)
# Test for (extant, count_with_alive_progeny) = (extant, attempts)
# For different minimum sizes we plot 
# (mean generation number, count_with_alive_progeny) with fixed death probability

# Note that extinct trees are discarded
print("Note that extinct trees are discarded!")

#variable parameter
min_sizes = [32, 64, 128]

#fixed death probability
death_prob = 0.4
print("Fixed death probability:", death_prob)

# Repeat simulations N times to get mean and devSt
N=1000

# results lists
mean_count = []
count_std = []
mean_gen = []
gen_std = []

# set seed to have reproducible results
np.random.seed(0)

# simulate
for min_size in min_sizes:
       
    # sigle runs results
    count = np.array([])
    gen = np.array([])
        
    for _ in range(N):
        tree = LD(bases = 1E+8, death_prob = death_prob, min_size = min_size, accept_extinct=False)
        tree.count_with_alive_progeny()
        count = np.append(count, tree._count_with_alive_progeny)
        gen = np.append(gen, tree._num_layers)
    
    #save this min size result
    mean_count.append(np.mean(count))
    count_std.append(np.std(count))
    mean_gen.append(np.mean(gen))
    gen_std.append(np.std(gen))   
    
#analytic estimation
estimated = []
for t in mean_gen:
    attempts = 0.
    # mean number of cell with at least one alive daughter
    # from generation zero to the second-last one
    attempts = (1-math.pow(death_prob,2))*1. * (math.pow(2*(1-death_prob),t-1) -1) / (math.log(2*(1-death_prob)))
    # we add also the mean number of cells in the last layer 
    attempts += 1. * math.pow(2*(1-death_prob), t)
    estimated.append(attempts)
    

#plot results
fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('Numero di generazioni', fontsize=40)
ax.set_ylabel('Cellule con progenie viva', fontsize=40)
ax.yaxis.offsetText.set_fontsize(40)

### set x y range
#ax.set_xlim(12.3, 33)
#ax.set_ylim(100, 700)

## set ticks
ax.tick_params(which='major', width=2.0, labelsize=40)
ax.tick_params(which='major', length=20, labelsize=40)


for i in range(len(min_sizes)):
    ax.errorbar(mean_gen[i], mean_count[i], yerr=count_std[i], xerr=gen_std[i],
         marker ='o', capsize=12, capthick=6, lw=6,fmt=' ',markersize=20,
         label='min size '+str(min_sizes[i]))
    
ax.scatter(mean_gen, estimated, marker = 'x', s=800, label = 'Stima analitica', color = 'tab:red',zorder=10)

ax.legend(fontsize=40, title='Legenda',title_fontsize=40)
plt.tight_layout()
plt.grid()
plt.title("d/(b+d) = "+str(death_prob), fontsize=40)
#plt.savefig("attempts_size_integral_0.4.pdf")
plt.show()

print("Mean count:")
print(mean_count)
print("Std count:")
print(count_std)
print("Mean gen:")
print(mean_gen)
print("Std gen:")
print(gen_std)
print("Estimated:")
print(estimated)<jupyter_output>Note that extinct trees are discarded!
Fixed death probability: 0.4<jupyter_text>Here trees are built until a fixed number of generations.Section 3.5.3, figure 12.<jupyter_code># Instead of using a minimum size to build the tree,
# a test with fixed numeber of generations is made here.

# fixed number of generations (time max)
max_gen = 10
# Number of simulations
N = 1000
# Death probability
death_prob = 0.1

print("Extinct trees are accepted!")
    
mean_count = np.zeros(max_gen)
std_count = np.zeros(max_gen)
estimated = np.zeros(max_gen)
    
for gen in range(max_gen) :
        
    this_gen_this_prob_count = []
        
    for _ in range (N):
        tree = LD(bases = 1E+8, death_prob = death_prob, gen=gen+1, accept_extinct=True)
        tree.count_with_alive_progeny()
        this_gen_this_prob_count.append(tree._count_with_alive_progeny)
            
    mean_count[gen] = np.mean(this_gen_this_prob_count)
    std_count[gen] = np.std(this_gen_this_prob_count)
        
    # Analytic estimation
    # In formulas t = gen, but here gen starts from zero
    t = gen + 1
    attempts = 0.
    # mean number of cell with at least one alive daughter
    # from generation zero to the second-last one
    attempts = (1-math.pow(death_prob,2))*1. * (math.pow(2*(1.-death_prob),t-1) -1) / (math.log(2*(1-death_prob)))
    # we add also the mean number of cells in the last layer 
    attempts += 1. * math.pow(2*(1-death_prob), t)
    estimated[gen] = attempts

               
#plot results
#plot results
fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('Numero di generazioni', fontsize=40)
ax.set_ylabel('Cellule con progenie viva', fontsize=40)
ax.yaxis.offsetText.set_fontsize(40)

### set x y range
#ax.set_xlim(0.5, 10.5)
#ax.set_ylim(-20, 1150)

## set ticksprobability
ax.tick_params(which='major', width=2.0, labelsize=40)
ax.tick_params(which='major', length=20, labelsize=40)

ax.errorbar(np.arange(1, max_gen+1), mean_count, yerr=std_count,
         marker ='o', c='tab:blue', capsize=12, capthick=6, lw=6,fmt=' ',markersize=20,
         label='Simulazione')
    
plt.scatter(np.arange(1, max_gen+1), estimated, marker = 'x', s=800, label = 'Stima analitica', color = 'tab:red',zorder=10)

ax.legend(fontsize=40, title='Legenda',title_fontsize=40)
plt.tight_layout()
plt.grid()
plt.title("d/(b+d) = "+str(death_prob), fontsize=40)


#plt.savefig("attempts_gen_integral_0.1.pdf")
plt.show()

print("Mean count:")
print(mean_count)
print("Std count:")
print(std_count)
print("Estimated:")
print(estimated)
    
# Instead of using a minimum size to build the tree,
# a test with fixed numeber of generations is made here.

# fixed number of generations (time max)
max_gen = 10
# Number of simulations
N = 1000
# Death probability
death_prob = 0.2

print("Extinct trees are accepted!")
    
mean_count = np.zeros(max_gen)
std_count = np.zeros(max_gen)
estimated = np.zeros(max_gen)
    
for gen in range(max_gen) :
        
    this_gen_this_prob_count = []
        
    for _ in range (N):
        tree = LD(bases = 1E+8, death_prob = death_prob, gen=gen+1, accept_extinct=True)
        tree.count_with_alive_progeny()
        this_gen_this_prob_count.append(tree._count_with_alive_progeny)
            
    mean_count[gen] = np.mean(this_gen_this_prob_count)
    std_count[gen] = np.std(this_gen_this_prob_count)
        
    # Analytic estimation
    # In formulas t = gen, but here gen starts from zero
    t = gen + 1
    attempts = 0.
    # mean number of cell with at least one alive daughter
    # from generation zero to the second-last one
    attempts = (1-math.pow(death_prob,2))*1. * (math.pow(2*(1.-death_prob),t-1) -1) / (math.log(2*(1-death_prob)))
    # we add also the mean number of cells in the last layer 
    attempts += 1. * math.pow(2*(1-death_prob), t)
    estimated[gen] = attempts

               
#plot results
fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('Numero di generazioni', fontsize=40)
ax.set_ylabel('Cellule con progenie viva', fontsize=40)
ax.yaxis.offsetText.set_fontsize(40)

### set x y range
#ax.set_xlim(0.5, 10.5)
#ax.set_ylim(-20, 500)

## set ticks
ax.tick_params(which='major', width=2.0, labelsize=40)
ax.tick_params(which='major', length=20, labelsize=40)

ax.errorbar(np.arange(1, max_gen+1), mean_count, yerr=std_count,
         marker ='o', c='tab:blue', capsize=12, capthick=6, lw=6,fmt=' ',markersize=20,
         label='Simulazione')
    
plt.scatter(np.arange(1, max_gen+1), estimated, marker = 'x', s=800, label = 'Stima analitica', color = 'tab:red',zorder=10)

ax.legend(fontsize=40, title='Legenda',title_fontsize=40)
plt.tight_layout()
plt.grid()
plt.title("d/(b+d) = "+str(death_prob), fontsize=40)
#plt.savefig("attempts_gen_integral_0.2.pdf")
plt.show()

print("Mean count:")
print(mean_count)
print("Std count:")
print(std_count)
print("Estimated:")
print(estimated)
    
# Instead of using a minimum size to build the tree,
# a test with fixed numeber of generations is made here.

# fixed number of generations (time max)
max_gen = 10
# Number of simulations
N = 1000
# Death probability
death_prob = 0.3

print("Extinct trees are accepted!")
    
mean_count = np.zeros(max_gen)
std_count = np.zeros(max_gen)
estimated = np.zeros(max_gen)
    
for gen in range(max_gen) :
        
    this_gen_this_prob_count = []
        
    for _ in range (N):
        tree = LD(bases = 1E+8, death_prob = death_prob, gen=gen+1, accept_extinct=True)
        tree.count_with_alive_progeny()
        this_gen_this_prob_count.append(tree._count_with_alive_progeny)
            
    mean_count[gen] = np.mean(this_gen_this_prob_count)
    std_count[gen] = np.std(this_gen_this_prob_count)
        
    # Analytic estimation
    # In formulas t = gen, but here gen starts from zero
    t = gen + 1
    attempts = 0.
    # mean number of cell with at least one alive daughter
    # from generation zero to the second-last one
    attempts = (1-math.pow(death_prob,2))*1. * (math.pow(2*(1.-death_prob),t-1) -1) / (math.log(2*(1-death_prob)))
    # we add also the mean number of cells in the last layer 
    attempts += 1. * math.pow(2*(1-death_prob), t)
    estimated[gen] = attempts

               
#plot results
fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('Numero di generazioni', fontsize=40)
ax.set_ylabel('Cellule con progenie viva', fontsize=40)
ax.yaxis.offsetText.set_fontsize(40)

### set x y range
#ax.set_xlim(0.5, 10.5)
#ax.set_ylim(-20, 170)

## set ticks
ax.tick_params(which='major', width=2.0, labelsize=40)
ax.tick_params(which='major', length=20, labelsize=40)

ax.errorbar(np.arange(1, max_gen+1), mean_count, yerr=std_count,
         marker ='o', c='tab:blue', capsize=12, capthick=6, lw=6,fmt=' ',markersize=20,
         label='Simulazione')
    
plt.scatter(np.arange(1, max_gen+1), estimated, marker = 'x', s=800, label = 'Stima analitica', color = 'tab:red',zorder=10)

ax.legend(fontsize=40, title='Legenda',title_fontsize=40)
plt.tight_layout()
plt.grid()
plt.title("d/(b+d) = "+str(death_prob), fontsize=40)

#plt.savefig("attempts_gen_integral_0.3.pdf")
plt.show()

print("Mean count:")
print(mean_count)
print("Std count:")
print(std_count)
print("Estimated:")
print(estimated)
    
# Instead of using a minimum size to build the tree,
# a test with fixed numeber of generations is made here.

# fixed number of generations (time max)
max_gen = 10
# Number of simulations
N = 1000
# Death probability
death_prob = 0.4

print("Extinct trees are accepted!")
    
mean_count = np.zeros(max_gen)
std_count = np.zeros(max_gen)
estimated = np.zeros(max_gen)
    
for gen in range(max_gen) :
        
    this_gen_this_prob_count = []
        
    for _ in range (N):
        tree = LD(bases = 1E+8, death_prob = death_prob, gen=gen+1, accept_extinct=True)
        tree.count_with_alive_progeny()
        this_gen_this_prob_count.append(tree._count_with_alive_progeny)
            
    mean_count[gen] = np.mean(this_gen_this_prob_count)
    std_count[gen] = np.std(this_gen_this_prob_count)
        
    # Analytic estimation
    # In formulas t = gen, but here gen starts from zero
    t = gen + 1
    attempts = 0.
    # mean number of cell with at least one alive daughter
    # from generation zero to the second-last one
    attempts = (1-math.pow(death_prob,2))*1. * (math.pow(2*(1.-death_prob),t-1) -1) / (math.log(2*(1-death_prob)))
    # we add also the mean number of cells in the last layer 
    attempts += 1. * math.pow(2*(1-death_prob), t)
    estimated[gen] = attempts

               
#plot results
fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('Numero di generazioni', fontsize=40)
ax.set_ylabel('Cellule con progenie viva', fontsize=40)
ax.yaxis.offsetText.set_fontsize(40)

### set x y range
#ax.set_xlim(0.5, 10.5)
#ax.set_ylim(-8, 55)

## set ticks
ax.tick_params(which='major', width=2.0, labelsize=40)
ax.tick_params(which='major', length=20, labelsize=40)

ax.errorbar(np.arange(1, max_gen+1), mean_count, yerr=std_count,
         marker ='o', c='tab:blue', capsize=12, capthick=6, lw=6,fmt=' ',markersize=20,
         label='Simulazione')
    
plt.scatter(np.arange(1, max_gen+1), estimated, marker = 'x', s=800, label = 'Stima analitica', color = 'tab:red',zorder=10)

ax.legend(fontsize=40, title='Legenda',title_fontsize=40)
plt.tight_layout()
plt.grid()
plt.title("d/(b+d) = "+str(death_prob), fontsize=40)

#plt.savefig("attempts_gen_integral_0.4.pdf")
plt.show()

print("Mean count:")
print(mean_count)
print("Std count:")
print(std_count)
print("Estimated:")
print(estimated)<jupyter_output>Extinct trees are accepted!<jupyter_text>Test LD estimator corrected by the continuos time modelSection 3.5.3, figure 13.<jupyter_code># Repeat simulations N times to get mean estimated mut rate and devSt

N = 10000 # number of simulation for each death probability
# N = 1000, time ~ 15 s
# N = 10^4, time ~ 150 s

# variable parameter
death_probabilities = np.linspace(0., 0.4, 9)

# results lists
# mean values
est_mu_no_dead = list([]) # using known number of cell with alive progeny
est_mu_no_dead_estimated = list([]) # using analytic estimate 
# standard deviations
est_mu_no_dead_std = list([])
est_mu_no_dead_estimated_std = list([])

# set seed to have reproducible results
np.random.seed(0)

# simulate
for death_prob in death_probabilities:
    # sigle runs results
    mu_no_dead = np.array([])
    mu_no_dead_estimated = np.array([])

    # perform N single runs
    for _ in np.arange(N):
        tree = LD(bases = 1E+8, death_prob = death_prob, min_size = 32)
        tree.test_LD_estimator_no_dead()
        tree.test_LD_estimator_no_dead_estimated()
        
        #save single run results
        mu_no_dead = np.append(mu_no_dead, tree._mu_est_no_dead)
        mu_no_dead_estimated = np.append(mu_no_dead_estimated, tree._mu_est_no_dead_estimated)
        
    # save results from N runs    
    est_mu_no_dead.append(np.mean(mu_no_dead))
    est_mu_no_dead_std.append(np.std(mu_no_dead))
    
    est_mu_no_dead_estimated.append(np.mean(mu_no_dead_estimated))
    est_mu_no_dead_estimated_std.append(np.std(mu_no_dead_estimated))    

    
#plot results
fig = plt.figure(figsize=(10,7))
ax = fig.add_subplot(1, 1, 1)
#plt.title('Stimatore corretto analiticamente\n'+str(N)+' simulazioni ogni punto')
ax.set_xlabel(r'$d/(b+d)$', fontsize=20)
ax.set_ylabel(r'$\mu (gen)^{-1}$', fontsize=20)
ax.yaxis.offsetText.set_fontsize(20)

ax.set_ylim(1.5*10**-9, 5.75*10**-9)

## set ticks
ax.tick_params(which='major', width=1.0, labelsize=20)
ax.tick_params(which='major', length=10, labelsize=20)

ax.errorbar(death_probabilities, est_mu_no_dead_estimated, yerr=est_mu_no_dead_estimated_std,
         marker ='o', c='tab:orange', capsize=6, capthick=3, lw=3,fmt=' ',markersize=10,
         label='Stimatore analitico')

plt.errorbar(death_probabilities, est_mu_no_dead, yerr=est_mu_no_dead_std,
         marker ='o', c='tab:blue', capsize=6, capthick=3, lw=3,fmt=' ',markersize=10,
         label='Stimatore corretto')



ax.axhline(y=tree._mu,  label='Reale', c='tab:gray',lw=3, ls=':')
ax.legend(fontsize=20, title='Legenda',title_fontsize=20)
plt.tight_layout()

#plt.savefig('dead_cells_estimated_correction_integral.pdf')
plt.show()     

# print results
print("death probabilities")
print(death_probabilities)
print("\nestimated mu no-dead")
print(est_mu_no_dead)
print("\nestimated mu no-dead std")
print(est_mu_no_dead_std)
print("\nestimated mu no-dead estimated")
print(est_mu_no_dead_estimated)
print("\nestimated mu no-dead estimated std")
print(est_mu_no_dead_estimated_std)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup