diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..c1d6aa5 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +.venv +.vscode diff --git a/figure_one/ancestors_number.py b/figure_one/ancestors_number.py new file mode 100644 index 0000000..2feed8f --- /dev/null +++ b/figure_one/ancestors_number.py @@ -0,0 +1,130 @@ +import msprime +import concurrent +import pandas as pd +import matplotlib.pyplot as plt +import time +import numpy as np + +sample_size = 10000 +r = 1e-9 +Ne = 1e6 +L = 1e6 +filename = 'lineages.csv' + +models = {'Hudson':'Hudson', + 'SMCK(1)': msprime.SMCK(1), + } + +def csv(x): + return ",".join(map(str, x)) + "\n" + +def save(name): + plt.tight_layout() + plt.savefig(f"figures/{name}.pdf") + +def get_ancestors_after_n_generations(params): + model = params[0]; generations = params[1] + + model_class = models[model] + sim = msprime.ancestry._parse_simulate( + sample_size=sample_size, + recombination_rate=r, + Ne=Ne, + length=L, + model=model_class, + end_time=generations, + ) + sim.run() + ancestors = len(sim.ancestors) + finish_time = sim.time + return [Ne, L, r, sample_size, str(model), generations, finish_time, ancestors] + +def generate_data(gens_list=[1,10,100,1000,10000,100000,100000,1000000, None], reps=25): + tasks = [] + for gen in gens_list: + for model in models.keys(): + for rep in range(reps): + tasks.append((model, gen)) + + with open(filename, "w") as f: + f.write(csv(['N', 'L', 'r','num_samples', 'model', 'gens', 'finish_time', 'ancestors'])) + + + with concurrent.futures.ProcessPoolExecutor(max_workers=12) as executor: + # Submit all tasks and get futures + future_results = {executor.submit(get_ancestors_after_n_generations, params): params for params in tasks} + + # Process results as they complete + for future in concurrent.futures.as_completed(future_results): + try: + result = future.result() + f.write(csv(result)) + f.flush() + except Exception as exc: + params = future_results[future] + print(f"Task {params} generated an exception: {exc}") + +def plot(): + df = pd.read_csv(filename) + + models = sorted(df['model'].unique()) + + + + # Define a colormap for sample sizes + cmap = plt.get_cmap("viridis").resampled(len(models)) + + # Create the plot + plt.figure(figsize=(10, 6)) + + for i, model in enumerate(models): + subset = df[df['model'] == model] + subset = subset[subset['ancestors'] != 0] + grouped_data = subset.groupby('gens')['ancestors'].mean().reset_index() + grouped_data = grouped_data.sort_values('gens') + xs = list(grouped_data['gens']) + ys = list(grouped_data['ancestors']) + + # deal with simulations to the end of time + subset = df[df['gens'].isna()] + grouped_data = ( + subset[subset['model'] == model] + .groupby('model')[['ancestors', 'finish_time']] + .mean() + .reset_index() + ) + grouped_data = grouped_data.sort_values('finish_time') + + xs += list(grouped_data['finish_time']) + ys += list(grouped_data['ancestors']) + + + plt.plot(xs, ys, linestyle='-', color=cmap(i), + linewidth=2, marker='o', markersize=5, label=f'{model}', alpha=0.5) + + + plt.xscale('log') + #plt.yscale('log') + + # Add labels and title + plt.xlabel('generations') + plt.title(f'Number of Ancestors after n generations. Ne:{Ne}, L:{L}, r:{r}, num_samples:{sample_size}') + + # Add legend + plt.legend() + + # Add grid for better readability + plt.grid(True, which="both", ls="--", alpha=0.3) + + # Save the plot + plt.tight_layout() + save(filename.split('.')[0]) + plt.close() + #plt.show() + +if __name__ == "__main__": + gens_list= np.logspace(0, 4, num=8, dtype=int).tolist() + \ + np.logspace(4, 6, num=6, dtype=int).tolist() +\ + np.logspace(6, 7, num=5, dtype=int).tolist() + [None] + generate_data(gens_list=gens_list, reps=25) + plot() diff --git a/figure_one/archive/draw-figure.py b/figure_one/archive/draw-figure.py new file mode 100644 index 0000000..c19daf4 --- /dev/null +++ b/figure_one/archive/draw-figure.py @@ -0,0 +1,41 @@ +import msprime +import numpy as np +import matplotlib.pyplot as plt +import concurrent.futures + +# Parameters +sample_size = 10 +sequence_length = 1e6 # total length of genome +recombination_rate = 1e-2 +num_replicates = 100 + +def simulate_tmrcas(model): + ts = msprime.sim_ancestry( + samples=sample_size, + sequence_length=sequence_length, + recombination_rate=recombination_rate, + model=model, + ) + # Return the TMRCA for all trees in this replicate + return [tree.time(tree.root) for tree in ts.trees()] + +def parallel_tmrcas(model): + tmrcas = [] + + with concurrent.futures.ProcessPoolExecutor() as executor: + results = executor.map(simulate_tmrcas, [model] * num_replicates) + for tmrcas_per_replicate in results: + tmrcas.extend(tmrcas_per_replicate) + return tmrcas + +# Standard model +tmrcas_standard = parallel_tmrcas(msprime.StandardCoalescent()) # standard coalescent +print(f"Standard model finished") +# SMC' model +tmrcas_smck = parallel_tmrcas(msprime.SmcKApproxCoalescent()) + +# Box plot +plt.boxplot([tmrcas_standard, tmrcas_smck], labels=["Standard", "SMC'"], showfliers=False) +plt.ylabel("TMRCA") +plt.title("Distribution of Coalescence Times") +plt.show() \ No newline at end of file diff --git a/figure_one/archive/generate-data.py b/figure_one/archive/generate-data.py new file mode 100644 index 0000000..ad0eb2c --- /dev/null +++ b/figure_one/archive/generate-data.py @@ -0,0 +1,82 @@ +import msprime +import numpy as np + +# Parameters +sample_size = 10 +sequence_length = 1e5 # total length of genome +recombination_rate = 1e-3 +mutation_rate = 0 # not needed here +distance = 1000 # distance between sites (in bp) +num_replicates = 1000 # to average over simulations + +def average_coalescence_time_at_distance_d(ts, d): + times = [] + positions = np.arange(0, ts.sequence_length - d, d) + for pos in positions: + left_ts = ts.at(pos) + right_ts = ts.at(pos + d) + + for tree in left_ts.trees(): + if tree in right_ts: + print(f"Tree {tree} is in both left and right trees") + + # Ensure both positions are covered by trees + if left_tree is None or right_tree is None: + continue + + # Pick pairs of samples + for i in range(0, ts.num_samples, 2): + if i + 1 >= ts.num_samples: + break + n1, n2 = i, i + 1 + + # Get TMRCA at the two positions + t1 = left_tree.tmrca(n1, n2) + t2 = right_tree.tmrca(n1, n2) + + # Average TMRCA for the pair of positions + times.append((t1 + t2) / 2) + + return np.mean(times) if times else np.nan + +def _average_coalescence_time_at_distance_d(ts, d): + times = [] + positions = np.arange(0, ts.sequence_length - d, d) + for pos in positions: + left_tree = ts.at(pos) + right_tree = ts.at(pos + d) + + # Ensure both positions are covered by trees + if left_tree is None or right_tree is None: + continue + + # Pick pairs of samples + for i in range(0, ts.num_samples, 2): + if i + 1 >= ts.num_samples: + break + n1, n2 = i, i + 1 + + # Get TMRCA at the two positions + t1 = left_tree.tmrca(n1, n2) + t2 = right_tree.tmrca(n1, n2) + + # Average TMRCA for the pair of positions + times.append((t1 + t2) / 2) + + return np.mean(times) if times else np.nan + +# Run simulation +avg_tmrcas = [] +for _ in range(num_replicates): + ts = msprime.sim_ancestry( + samples=sample_size, + sequence_length=sequence_length, + recombination_rate=recombination_rate, + random_seed=None + ) + avg_tmrcas.append(average_coalescence_time_at_distance_d(ts, distance)) + +# Filter out None values +avg_tmrcas = [tmrca for tmrca in avg_tmrcas if tmrca is not None] +overall_average = np.nanmean(avg_tmrcas) +print(f"Average coalescence time for site pairs {distance} bp apart: {overall_average}") diff --git a/figure_one/speed.py b/figure_one/speed.py new file mode 100644 index 0000000..acfdf49 --- /dev/null +++ b/figure_one/speed.py @@ -0,0 +1,269 @@ +''' +blue two panels never used +''' + +import msprime +import concurrent +import pandas as pd +import matplotlib.pyplot as plt +import time +import numpy as np +import ast +import warnings +warnings.filterwarnings("ignore") + +max_workers=14 +filename = f'speed.csv' +replicates = 25 +Ne = 1e6 + +seq_len_dor = 25400000 +recombination_rate = 1.045e-8 #2.40463e-08 +sample_sizes = [2, 4, 10, 100, 1000] + +lengths = np.logspace(1, 7, num=7, dtype=int) +lengths = np.append(lengths, seq_len_dor) +shortened_lengths = lengths[lengths <= 1e6] +models = {'Hudson':'Hudson', + 'SMC(500k)': msprime.SMCK(500000), + 'SMC(1)': msprime.SMCK(1), + 'SMC(0)': msprime.SMCK(0) + } + +models_for_sample_size = {'Hudson':'Hudson', + 'SMC(1)': msprime.SMCK(1)} + +neL = Ne * lengths +shortened_neL = Ne * shortened_lengths +drosophila_neL = Ne * seq_len_dor +human_neL = 10**4 * 248956422 + +def csv(x): + return ",".join(map(str, x)) + "\n" + +def save(name): + plt.tight_layout() + #plt.savefig(f"figures/{name}.png") + plt.savefig(f"figures/{name}.pdf") + +def get_exc_time(params): + model = params[0]; length = params[1]; sample_size = params[2] + model_class = models[model] + start_time = time.time() + ts = msprime.sim_ancestry( + samples=sample_size, + ploidy=2, + sequence_length=length, + recombination_rate=recombination_rate, + population_size=Ne, + model=model_class + ) + ex_time = time.time() - start_time + + return [Ne, length, recombination_rate, sample_size, str(model), ex_time] + +def generate_data(): + + tasks = [] + + #loop for speed test per model and per length + for model in models: + for replicate in range(replicates): + if model not in ['SMC(1)', 'SMC(0)']: + allowed_L = shortened_lengths + else: + allowed_L = lengths + + for length in allowed_L: + + if model in models_for_sample_size: + for sample_size in sample_sizes: + tasks.append((model, length, sample_size)) + else: + tasks.append((model, length, 2)) + + with open(filename, "w") as f: + f.write(csv(['N', 'L', 'r','num_samples', 'model', 'ex_time'])) + + with concurrent.futures.ProcessPoolExecutor(max_workers=max_workers) as executor: + # Submit all tasks and get futures + future_results = {executor.submit(get_exc_time, params): params for params in tasks} + + # Process results as they complete + for future in concurrent.futures.as_completed(future_results): + try: + result = future.result() + f.write(csv(result)) + f.flush() + except Exception as exc: + params = future_results[future] + print(f"Task {params} generated an exception: {exc}") + +def format_time_label(time) -> str: + """Format time in seconds to a more readable string.""" + time = float(time) + if time < 60: + return f'{time:.1f}s' + elif time < 3600: + return f'{time/60:.1f}m' + elif time < 86400: + return f'{time/3600:.1f}h' + else: + return f'{time/86400:.1f}d' + +def text_on_plot_right(ax, y, text, color): + ax.text(1.02, y, text, + transform=ax.get_yaxis_transform(), + fontweight='bold', + va='center', fontsize=9, color=color) + +def text_on_plot_top(ax, x, text, color): + ax.text(x, 1.005, text, + transform=ax.get_xaxis_transform(), + #fontweight='bold', + ha='center', va='bottom', fontsize=9, color=color) + +def plot_speed(infile=filename): + + df = pd.read_csv(infile) + df = df[df['num_samples']==2] + + df_avg = df.groupby(['model', 'L']).mean().reset_index() + #shortened_neL = Ne * lengths[lengths <= 1e6] + hudson_times = df_avg[df_avg['model']=='Hudson']['ex_time'].to_numpy() + fit_times = np.polyfit(shortened_neL, hudson_times, 2) + fit_fn = np.poly1d(fit_times) + np.log10(fit_fn(neL[-1])) + fitted_line = fit_fn(neL[1:]) + + fig, ax = plt.subplots() + + for model in models_for_sample_size: + model_times = df_avg[df_avg['model']==model]['ex_time'].to_numpy() + xs = neL[:len(model_times)] + line = ax.plot(xs, model_times, marker='o', label=model) + + final_time = model_times[-1] + + if len(model_times) == len(lengths): + time_label = format_time_label(final_time) + text_on_plot_right(ax, final_time, time_label, line[0].get_color()) + + ax.plot(neL[1:], fitted_line, linestyle='--', color='gray', label='Quadratic fit (Hudson)') + final_fitted_time = fitted_line[-1] + fitted_time_label = format_time_label(final_fitted_time) + + text_on_plot_right(ax, final_fitted_time, fitted_time_label, 'gray') + + ax.axvline(x=drosophila_neL, color='green', linestyle=':', linewidth=3) + text_on_plot_top(ax, drosophila_neL, "Drosophila\n(chrom 3R)", "green") + ax.axvline(x=human_neL, color='purple', linestyle=':', linewidth=3) + text_on_plot_top(ax, human_neL, "Human\n(chrom 1)", "purple") + + ax.set_xscale('log') + ax.set_yscale('log') + ax.set_xlabel('Population-scaled Sequence Length (Ne * L)') + #ax.set_ylabel('Execution Time (seconds)') + #ax.set_title('SMC(k) vs Hudson Execution Time') + ax.set_title('Execution Time (seconds)', ha='right') + ax.legend() + #plt.grid(True, which="both", ls="--") + + plt.tight_layout() + plt.subplots_adjust(right=0.90) # Make room for the labels on the right + + save('speed') + plt.clf() + +def plot_speed_per_model(infile=filename): + + df = pd.read_csv(infile) + df = df[df['num_samples']==2] + + df_avg = df.groupby(['model', 'L']).mean().reset_index() + + hudson_times = df_avg[df_avg['model']=='Hudson']['ex_time'].to_numpy() + fit_times = np.polyfit(shortened_neL, hudson_times, 2) + fit_fn = np.poly1d(fit_times) + np.log10(fit_fn(neL[-1])) + fitted_line = fit_fn(neL[1:]) + + fig, ax = plt.subplots() + + for model in models: + model_times = df_avg[df_avg['model']==model]['ex_time'].to_numpy() + xs = neL[:len(model_times)] + line = ax.plot(xs, model_times, marker='o', label=model) + + final_time = model_times[-1] + if len(model_times) == len(lengths): + time_label = format_time_label(final_time) + text_on_plot_right(ax, final_time, time_label, line[0].get_color()) + + ax.plot(neL[1:], fitted_line, linestyle='--', color='gray', label='Quadratic fit (Hudson)') + final_fitted_time = fitted_line[-1] + fitted_time_label = format_time_label(final_fitted_time) + + text_on_plot_right(ax, final_fitted_time, fitted_time_label, 'gray') + + ax.axvline(x=drosophila_neL, color='green', linestyle=':', linewidth=3) + text_on_plot_top(ax, drosophila_neL, "Drosophila\n(chrom 3R)", "green") + ax.axvline(x=human_neL, color='purple', linestyle=':', linewidth=3) + text_on_plot_top(ax, human_neL, "Human\n(chrom 1)", "purple") + + ax.set_xscale('log') + ax.set_yscale('log') + ax.set_xlabel('Population-scaled Sequence Length (Ne * L)') + #ax.set_ylabel('Execution Time (seconds)') + #ax.set_title('SMC(k) vs Hudson Execution Time') + ax.set_title('Execution Time (seconds)', ha='right') + + ax.legend() + #plt.grid(True, which="both", ls="--") + + plt.tight_layout() + plt.subplots_adjust(right=0.90) # Make room for the labels on the right + + save('speed_per_model') + plt.clf() + +def plot_speed_per_sample_size(infile=filename): + df = pd.read_csv(infile) + df = df[df['model'].isin(models_for_sample_size.keys())] + + df_avg = df.groupby(['model', 'L', 'num_samples']).mean().reset_index() + + fig, ax = plt.subplots() + for model in models_for_sample_size: + if model not in ['SMC(1)']: continue + for sample_size in sample_sizes: + model_times = df_avg[(df_avg['model']==model) & (df_avg['num_samples']==sample_size)]['ex_time'].to_numpy() + xs = neL[:len(model_times)] + line = ax.plot(xs, model_times, marker='o', label=f"{model}, n={sample_size}") + + if len(model_times) == len(lengths): + final_time = model_times[-1] + time_label = format_time_label(final_time) + text_on_plot_right(ax, final_time, time_label, line[0].get_color()) + + ax.axvline(x=drosophila_neL, color='black', linestyle=':', linewidth=3) + text_on_plot_top(ax, drosophila_neL, "Drosophila\n(chrom 3R)", "black") + ax.axvline(x=human_neL, color='grey', linestyle=':', linewidth=3) + text_on_plot_top(ax, human_neL, "Human\n(chrom 1)", "grey") + + ax.set_xscale('log') + ax.set_yscale('log') + ax.set_xlabel('Population-scaled Sequence Length (Ne * L)') + #plt.ylabel('Execution Time (seconds)') + ax.set_title('Execution Time (seconds)', ha='right') + ax.legend() + #plt.grid(True, which="both", ls="--") + save('speed_per_sample_size') + plt.clf() + + +if __name__ == "__main__": + #generate_data() + plot_speed() + plot_speed_per_model() + plot_speed_per_sample_size() diff --git a/figure_one/tscompare_plot.py b/figure_one/tscompare_plot.py new file mode 100644 index 0000000..bce4699 --- /dev/null +++ b/figure_one/tscompare_plot.py @@ -0,0 +1,225 @@ +import msprime +import concurrent +import pandas as pd +import matplotlib.pyplot as plt +import seaborn as sns +import time +import numpy as np +import ast +import tscompare +import warnings +warnings.filterwarnings("ignore") + +max_workers=14 +filename = f'diff_tscompare.csv' +replicates = 1000 +Ne = 1e6 +sample_size = 2 +recombination_rate =1.045e-8 +seq_len = 1e5 +models = {'Hudson':'Hudson', + #'k=500k': msprime.SMCK(500000), + #'k=100k': msprime.SMCK(100000), + 'SMC(1)': msprime.SMCK(1), + 'SMC(0)': msprime.SMCK(0), + #'k=10': msprime.SMCK(10), + #'k=100': msprime.SMCK(100), + 'SMC(1kb)': msprime.SMCK(1000), + 'SMC(10kb)': msprime.SMCK(10000), + } +#models_ordered = ['Hudson', 'k=0', 'k=1','k=10', 'k=100','k=1k', 'k=10k', 'k=100k', 'k=500k'] +models_ordered = ['Hudson', 'SMC(0)', 'SMC(1)','SMC(1kb)', 'SMC(10kb)'] + + + +def csv(x): + return ",".join(map(str, x)) + "\n" + +def save(name): + plt.tight_layout() + #plt.savefig(f"figures/{name}.png") + plt.savefig(f"figures/{name}.pdf") + +def get_exc_time(params): + print(params) + model = params[0] + model_class = models[model] + start_time = time.time() + ts_hudson = msprime.sim_ancestry( + samples=sample_size, + ploidy=2, + sequence_length=seq_len, + recombination_rate=recombination_rate, + population_size=Ne, + coalescing_segments_only=False, + ) + ts_model = msprime.sim_ancestry( + samples=sample_size, + ploidy=2, + sequence_length=seq_len, + recombination_rate=recombination_rate, + population_size=Ne, + model=model_class, + coalescing_segments_only=False, + ) + node_times_matched, _span, best_id = tscompare.match_node_ages(ts_hudson, ts_model) + all_smc_node_times = np.array([n.time for n in ts_model.nodes()]) + node_times_smc = all_smc_node_times[best_id] + + mask = ~np.isnan(node_times_matched) & ~np.isnan(node_times_smc) + node_times_matched = node_times_matched[mask] + node_times_smc = node_times_smc[mask] + assert (node_times_matched == node_times_smc).all() + + node_times_hudson = np.array([n.time for n in ts_hudson.nodes()]) + node_times_hudson = node_times_hudson[mask] + assert len(node_times_hudson) == len(node_times_matched) + diff = np.log(1 + node_times_hudson[sample_size:]) - np.log(1+ node_times_matched[sample_size:]) + masked_span = _span[mask][sample_size:] + diff_by_span = (diff * masked_span)/seq_len + rmse = np.sqrt(np.mean(diff_by_span**2)) + + + x= np.sqrt( + np.sum(diff ** 2 * masked_span) + / np.sum(masked_span) + ) + + dis = tscompare.haplotype_arf(ts_hudson, ts_model) + return [Ne, seq_len, recombination_rate, sample_size, str(model), dis.arf, dis.tpr, dis.matched_span[0], dis.matched_span[1], dis.rmse, f'\"{diff.tolist()}\"', rmse] + +def generate_data(): + + tasks = [] + + #loop for speed test per model and per length + for model in models: + for replicate in range(replicates): + tasks.append((model,)) + + with open(filename, "w") as f: + f.write(csv(['N', 'L', 'r','num_samples', 'model', 'arf', 'tpr', 'matched_span', 'inverse_matched_span', 'rmse', 'time_diffs', 'calc_rmse'])) + + with concurrent.futures.ProcessPoolExecutor(max_workers=max_workers) as executor: + # Submit all tasks and get futures + future_results = {executor.submit(get_exc_time, params): params for params in tasks} + + # Process results as they complete + for future in concurrent.futures.as_completed(future_results): + try: + result = future.result() + f.write(csv(result)) + f.flush() + except Exception as exc: + params = future_results[future] + print(f"Task {params} generated an exception: {exc}") + raise exc + +def plot(infile=filename): + df = pd.read_csv(infile) + + df['model'] = pd.Categorical(df['model'], categories=models_ordered, ordered=True) + df['time_diffs'] = df['time_diffs'].apply(ast.literal_eval) + + #a box plot of arf per model + plt.figure(figsize=(8,6)) + ax = plt.subplot(1,1,1) + #df.boxplot(column='arf', by='model', ax=ax) + sns.violinplot(data=df, x='model', y='arf', ax=ax, alpha=0.95, palette='Set3') + mean_hudson = df.loc[df['model'] == 'Hudson', 'arf'].median() + ax.axhline(mean_hudson, color='red', linestyle='--') + plt.title('TSCompare Robinson-Foulds relative dissimilarity') + plt.suptitle('') + plt.ylabel('Average RF Distance (ARF)') + save('tscompare_accuracy_comparison') + + #a box plot of tpr per model + plt.figure(figsize=(8,6)) + ax = plt.subplot(1,1,1) + sns.violinplot(data=df, x='model', y='tpr', ax=ax, alpha=0.95, palette='Set3') + mean_hudson = df.loc[df['model'] == 'Hudson', 'tpr'].median() + ax.axhline(mean_hudson, color='red', linestyle='--') + plt.title('TSCompare true proportion represented Comparison') + plt.suptitle('') + plt.ylabel('true proportion represented (TPR)') + save('tscompare_tpr_comparison') + + #a box plot of rmse per model + plt.figure(figsize=(8,6)) + ax = plt.subplot(1,1,1) + sns.violinplot(data=df, x='model', y='rmse', ax=ax, alpha=0.95, palette='Set3') + mean_hudson = df.loc[df['model'] == 'Hudson', 'rmse'].median() + ax.axhline(mean_hudson, color='red', linestyle='--') + plt.title('TSCompare RMSE Comparison') + plt.suptitle('') + plt.ylabel('Root Mean Square Error (RMSE)') + save('tscompare_rmse_comparison') + + df['sum_matched_span'] = df['matched_span'] + df['inverse_matched_span'] + #a box plot of matched span length per model + plt.figure(figsize=(8,6)) + ax = plt.subplot(1,1,1) + sns.violinplot(data=df, x='model', y='sum_matched_span', ax=ax, alpha=0.95, palette='Set3') + mean_hudson = df.loc[df['model'] == 'Hudson', 'sum_matched_span'].median() + ax.axhline(mean_hudson, color='red', linestyle='--') + plt.title('TSCompare Matched Span Length Comparison') + plt.suptitle('') + plt.ylabel('Matched Span Length + inverse_match') + save('tscompare_sum_matched_span_length_comparison') + + #a box plot of inverse matched span length per model + plt.figure(figsize=(8,6)) + ax = plt.subplot(1,1,1) + + mean_hudson = df.loc[df['model'] == 'Hudson', 'inverse_matched_span'].median() + df_norm = df.copy() + df_norm['inverse_matched_span_norm'] = df_norm['inverse_matched_span'] / mean_hudson + df_norm = df_norm[df_norm['model'] != 'Hudson'] + df_norm['model'] = df_norm['model'].cat.remove_categories('Hudson') + + sns.violinplot(data=df_norm, x='model', y='inverse_matched_span_norm', ax=ax, alpha=0.95, palette='Set3') + ax.axhline(1, color='red', linestyle='--') + plt.title('Normalised similarity ($\\it{tscompare}$ matched span)', loc='left', fontsize=16) + plt.suptitle('') + plt.ylabel('') + plt.xlabel('') + ax.tick_params(labelsize=14) + save('tscompare_inverse_matched_span_length_comparison') + + plt.figure(figsize=(8,6)) + ax = plt.subplot(1,1,1) + mean_hudson = df.loc[df['model'] == 'Hudson', 'matched_span'].median() + ax.axhline(mean_hudson, color='red', linestyle='--') + plt.title('TSCompare Matched Span Length Comparison') + plt.suptitle('') + plt.ylabel('Matched Span Length') + save('tscompare_matched_span_length_comparison') + + grouped = df.groupby('model').agg({'time_diffs': lambda x: sum(x, [])}).reset_index() + plt.figure(figsize=(8,6)) + ax = plt.subplot(1,1,1) + data_to_plot = grouped['time_diffs'].tolist() + labels = grouped['model'].tolist() + sns.violinplot(data_to_plot, ax=ax, alpha=0.95, palette='Set3') + plt.xticks(ticks=range(len(labels)), labels=labels) + + plt.title('TSCompare time differences per node') + plt.suptitle('') + plt.ylabel('Time Differences') + save('tscompare_time_differences_comparison') + + #a box plot of rmse per model + plt.figure(figsize=(8,6)) + ax = plt.subplot(1,1,1) + sns.violinplot(data=df, x='model', y='calc_rmse', ax=ax, alpha=0.95, palette='Set3') + mean_hudson = df.loc[df['model'] == 'Hudson', 'calc_rmse'].median() + ax.axhline(mean_hudson, color='red', linestyle='--') + plt.title('TSCompare RMSE Comparison') + plt.suptitle('') + plt.ylabel('Root Mean Square Error (RMSE)') + save('tscompare_calc_rmse_comparison') + +if __name__ == "__main__": + + #generate_data() + plot() diff --git a/figure_one/unary_nodes_average_hulls.py b/figure_one/unary_nodes_average_hulls.py new file mode 100644 index 0000000..131d4c9 --- /dev/null +++ b/figure_one/unary_nodes_average_hulls.py @@ -0,0 +1,324 @@ +import msprime +import concurrent +import pandas as pd +import matplotlib.pyplot as plt +from matplotlib.patches import Patch +from matplotlib.legend_handler import HandlerTuple +import time +import numpy as np +import ast +import warnings +warnings.filterwarnings("ignore") + +sample_size = 2 +r = 1e-8 +Ne = 1e6 +L = 1e6 +max_workers=3 +filename = f'average_hulls_{sample_size}samples_segments.csv' +models = {'CwR':'Hudson', + 'SMC(500kb)': msprime.SMCK(500000), + 'SMC(1)': msprime.SMCK(1), + 'SMC(0)': msprime.SMCK(0) + } + +rs = [1e-11, 1e-10, 1e-9] + + +def csv(x): + return ",".join(map(str, x)) + "\n" + +def save(name): + plt.tight_layout() + #plt.savefig(f"figures/{name}.png") + plt.savefig(f"figures/{name}.pdf") + +def get_hulls_after_n_generations(params): + model = params[0]; _sample_size = params[1]; _r = params[2] + model_class = models[model] + + ts = msprime.sim_ancestry( + samples=_sample_size, + recombination_rate=_r, + population_size=Ne, + sequence_length=L, + model=model_class, + #additional_nodes=(msprime.NodeType.COMMON_ANCESTOR), + coalescing_segments_only=False, + stop_at_local_mrca=False + + ) + num_trees = ts.num_trees + l1 = np.zeros(ts.num_nodes+1) + l2 = np.zeros(ts.num_nodes+1) + is_root = np.zeros(ts.num_nodes, dtype=bool) + for tree in ts.trees(): + l1[(tree.num_children_array == 1)] += tree.span + l2[(tree.num_children_array == 2)] += tree.span + root = tree.root + is_root[root] = True + + assert np.all(ts.samples() == np.arange(ts.num_samples)) + + total_span = l1 + l2 + + start = np.array([-1] * ts.num_nodes) + not_started = np.ones(ts.num_nodes, dtype=bool) + not_started[:ts.sample_size] = False + + for tree in ts.trees(): + tree_nodes = np.zeros(ts.num_nodes, dtype=bool) + tree_nodes[tree.preorder()] = True + + tree_nodes_did_not_start = np.zeros(ts.num_nodes, dtype=bool) + tree_nodes_did_not_start[not_started & tree_nodes] = True + + start[tree_nodes_did_not_start] = tree.interval[0] + not_started[tree_nodes] = False + if not (not_started[ts.sample_size:].any()): + break + + end = np.array([-1] * ts.num_nodes) + not_ended = np.ones(ts.num_nodes, dtype=bool) + not_ended[:ts.sample_size] = False + + for tree in reversed(ts.trees()): + tree_nodes = np.zeros(ts.num_nodes, dtype=bool) + tree_nodes[tree.preorder()] = True + + tree_nodes_did_not_end = np.zeros(ts.num_nodes, dtype=bool) + tree_nodes_did_not_end[not_ended & tree_nodes] = True + + end[tree_nodes_did_not_end] = tree.interval[1] + not_ended[tree_nodes] = False + if not (not_ended[ts.sample_size:].any()): + break + + youngest_root = np.where(is_root)[0][0] + hulls = end - start + df = pd.DataFrame({'node': np.arange(ts.num_nodes), 'hull': hulls, 'l1':l1[:-1], 'l2':l2[:-1], 'total_span': total_span[:-1]}) + keep_rows = np.ones(ts.num_nodes, dtype=bool) + keep_rows[ts.samples()] = False + keep_rows[is_root] = False + #keep_rows[np.arange(youngest_root, ts.num_nodes)] = False + + + df_spans = pd.DataFrame(0,index=np.arange(ts.num_nodes)[keep_rows], columns=np.arange(ts.num_trees)) + for i, tree in enumerate(ts.trees()): + tree_nodes = np.zeros(ts.num_nodes, dtype=bool) + tree_nodes[tree.preorder()] = True + tree_nodes = tree_nodes[keep_rows] + df_spans.loc[tree_nodes, i] = tree.interval[1] - tree.interval[0] + + + all_segments = [] + no_of_segments = [] + for u in df_spans.index: + spans = np.array(df_spans.loc[u]) + segments = np.split(spans, np.where(spans == 0)[0]) + #drop segments that only has zeroes + segments = [seg for seg in segments if not np.all(seg == 0)] + #drop 0s from each segment + segments = [seg[seg != 0] for seg in segments] + all_segments.extend([sum(seg) for seg in segments]) + no_of_segments.append(len(segments)) + + avg_all_segments = np.mean(all_segments) + df = df[keep_rows] + avg_hulls = df['hull'].mean() + avg_l1 = df['l1'].mean() + avg_l2 = df['l2'].mean() + trapped = df['hull'] - (df['l1'] + df['l2']) + avg_trapped = trapped.mean() + avg_no_of_segments = np.mean(no_of_segments) + + oldest_node_time = ts.nodes()[-1].time + + return [Ne, L, _r, _sample_size, str(model), str(avg_hulls), + str(avg_l1), str(avg_l2), str(avg_trapped), str(avg_all_segments), + str(avg_no_of_segments), str(oldest_node_time), str(num_trees)] + +def generate_data(replicates=5): + + tasks = [] + + for _r in rs: + for model in models: + for replicate in range(replicates): + tasks.append((model, sample_size, _r)) + + timeout_seconds = 3 * 60 # 3 minutes in seconds + + with open(filename, "w") as f: + f.write(csv(['N', 'L', 'r','num_samples', 'model', 'avg_hulls', 'avg_l1', 'avg_l2', 'avg_trapped', 'avg_adj_hap_len', 'no_of_segments', 'oldest_node_time', 'num_trees'])) + + with concurrent.futures.ProcessPoolExecutor(max_workers=max_workers) as executor: + # Submit all tasks and get futures + future_results = {executor.submit(get_hulls_after_n_generations, params): params for params in tasks} + + # Process results as they complete + for future in concurrent.futures.as_completed(future_results): + try: + result = future.result(timeout=timeout_seconds) + f.write(csv(result)) + f.flush() + except concurrent.futures.TimeoutError: + params = future_results[future] + print(f"Task {params} exceeded the timeout of {timeout_seconds} seconds.") + except Exception as exc: + params = future_results[future] + print(f"Task {params} generated an exception: {exc}") + +def plot_single_parameter(param='num_trees', log_f=False, title=None): + # Read CSV + import ast + + df = pd.read_csv(filename) + rs = sorted(df['r'].unique()) + x = np.arange(len(rs)) + + ne = df['N'][0] + assert np.all(df['N'] == ne), "N values are not consistent in the data." + + #models = sorted(df['model'].unique()) + cmap = plt.get_cmap('tab10') + + fig, ax = plt.subplots(1, 1, figsize=(10, 10), sharex=True) + + width = 0.8 / len(models) # space bars for each model at the same r + + # Panel 1: avg_l1 + for i, model in enumerate(models): + df_model = df[df['model'] == model] + grouped_in = df_model.groupby('r').apply(lambda x: x[param].mean()).reset_index(name=f'mean_{param}') + + color = cmap(i) + + offsets = x + i * width - (width * len(models)) / 2 + ax.bar(offsets, grouped_in[f'mean_{param}'], width, alpha=0.8, color=color, label=model) + r_to_x = dict(zip(rs, x)) + + scatter_x = df_model['r'].map(r_to_x) + i * width - (width * len(models)) / 2 + + ax.scatter(scatter_x, + df_model[param], + color=color, + alpha=0.6, + s=20, + linewidth=0.3, + zorder=3) + if title is not None: + plt.title(title, loc='left', fontsize=20) + else: + plt.title(f'{param}, sample size {sample_size}, Ne {ne}, L {L}', loc='left', fontsize=20) + + human_rho = ((1e-8)*4*1e4) + plt.grid(True) + plt.legend(fontsize=16) + ax.set_xlabel(r'Normalised recombination rate ($\rho / \rho_{\mathrm{human}}$)', fontsize=16) + x_ticks = [f"{(i*4*ne)/human_rho:.2g}" for i in rs] + ax.set_xticks(x) + ax.set_xticklabels(x_ticks) + ax.tick_params(labelsize=16) + + if log_f: + plt.yscale('log') + + out_file_name = f'single_param{filename.split(".")[0]}_{param}_log' + else: + out_file_name = f'single_param{filename.split(".")[0]}_{param}' + + + plt.tight_layout() + save(out_file_name) + plt.close() + +def plot_stacked_bars(infile=filename): + df = pd.read_csv(infile) + + # Filter on sample size + df_samples = df[df['num_samples'] == sample_size].copy() + + + # Group and aggregate + grouped = df_samples.groupby(['r', 'model']).agg({ + 'avg_hulls': 'mean', + 'avg_l1': 'mean', + 'avg_l2': 'mean', + 'avg_trapped': 'mean', + 'no_of_segments': 'mean' + }).reset_index() + + '''grouped['trapped'] = grouped.apply( + lambda row: np.array(row['hulls']) - (np.array(row['l1']) + np.array(row['l2'])), + axis=1 + )''' + grouped['avg_hulls'] = grouped['avg_hulls'].apply(lambda x: (x) / L) + grouped['avg_l1'] = grouped['avg_l1'].apply(lambda x: (x) / L) + grouped['avg_l2'] = grouped['avg_l2'].apply(lambda x: (x) / L) + grouped['avg_trapped'] = grouped['avg_trapped'].apply(lambda x: (x) / L) + # Keep no_of_segments unnormalized for the label + segments_for_label = grouped['no_of_segments'].copy() + grouped['no_of_segments'] = grouped['no_of_segments'].apply(lambda x: (x) / L) + + #models = sorted(grouped['model'].unique()) + rs = sorted(grouped['r'].unique()) + x = np.arange(len(rs)) + ne = df['N'][0] + assert np.all(df['N'] == ne), "N values are not consistent in the data." + #x = [i* 4 * ne for i in x] + + width = 0.8 / len(models) # space bars for each model at the same r + fig, ax = plt.subplots(figsize=(12, 6)) + cmap = plt.get_cmap('tab10') + + legend_handles = [] + legend_labels = [] + for i, model in enumerate(models): + color = cmap(i) + model_data = grouped[grouped['model'] == model] + model_segments = segments_for_label[grouped['model'] == model] + offsets = x + i * width - (width * len(models)) / 2 + ax.bar(offsets, model_data['avg_l1'], width, alpha=0.6) + ax.bar(offsets, model_data['avg_l2'], width, bottom=model_data['avg_l1'], color=color, alpha=0.9) + ax.bar(offsets, model_data['avg_trapped'], width, + bottom=model_data['avg_l1'] + model_data['avg_l2'], + color='gray', alpha=0.5) + + unary_patch = Patch(facecolor=color, alpha=0.6) + binary_patch = Patch(facecolor=color, alpha=1.0) + legend_handles.append((binary_patch, unary_patch)) + legend_labels.append(f"{model} binary / unary") + + legend_handles.append(Patch(facecolor='grey', alpha=1.0)) + legend_labels.append("Trapped material") + ax.legend(legend_handles, legend_labels, + handler_map={tuple: HandlerTuple(ndivide=None)}, + fontsize=16) + + human_rho = ((1e-8)*4*1e4) + x_ticks = [f"{(i*4*ne)/human_rho:.2g}" for i in rs] + ax.set_xticks(x) + ax.set_xticklabels(x_ticks) + ax.set_xlabel(r'Normalised recombination rate ($\rho / \rho_{\mathrm{human}}$)', fontsize=18) + ax.tick_params(labelsize=16) + + #ax.set_yscale('log') + ax.set_title('Normalised length (per L)', loc='left', fontsize=18) + #ax.set_title(f'Stacked bar of l1, l2, and trapped material per r\nSample size {sample_size}, Ne {Ne}, L {L}') + #ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + #ax.legend(fontsize=13) + ax.grid(True) + + plt.tight_layout() + save(f'{infile.split(".")[0]}_stacked_bar') + plt.close() + +if __name__ == "__main__": + #generate_data(replicates=10) + plot_stacked_bars() + #for param in ['num_trees']: + #plot_single_parameter(param=param) + param = 'num_trees' + title = "Number of trees making the ARG" + plot_single_parameter(param=param, log_f=True, title=title) diff --git a/figure_three/demes.yaml b/figure_three/demes.yaml new file mode 100644 index 0000000..0269f7f --- /dev/null +++ b/figure_three/demes.yaml @@ -0,0 +1,24 @@ +description: + 3 population IM model with migration from Iphiclides feisthamelii into I. podalirius. +time_units: generations +demes: + - name: ancestral + epochs: + - {start_size: 1.15E6, end_time: 2.18E6} + - name: IF + ancestors: [ancestral] + epochs: + - start_size: 4.83E5 + - name: IP + ancestors: [ancestral] + epochs: + - start_size: 3.77E5 +migrations: + - source: IF + dest: IP + rate: 4.73E-08 + + - source: IP + dest: IF + rate: 1.553e-06 + start_time: 275 \ No newline at end of file diff --git a/figure_three/nb.ipynb b/figure_three/nb.ipynb new file mode 100644 index 0000000..724d02e --- /dev/null +++ b/figure_three/nb.ipynb @@ -0,0 +1,223 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 4, + "id": "7f3e3d35", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0sAAAHWCAYAAACmIxLCAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAU/NJREFUeJzt3Xl8VPW9//H3ObNlIQkECGtABISqbBW1ICooLS7VUltqXYr7bb1QF6y2qK11K9brRq/Wpa1Sb0u1Vmt/LlURxRW1gAtoK6IIVFmLELLOcs7vj3NmMpMcNCeZZCbJ6/l4nEcyk5mTb77ny+HzPt/vzBi2bdsCAAAAAGQwc90AAAAAAMhHhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8NCtw9KLL76oE044QQMHDpRhGHr00Ud978O2bd10003ab7/9FIlENGjQIF1//fXZbywAAACADhXMdQNyqaamRuPGjdPZZ5+tk046qVX7uPDCC/XMM8/opptu0pgxY7Rz507t3Lkzyy0FAAAA0NEM27btXDciHxiGob/+9a+aOXNm6r6GhgZdccUV+tOf/qRdu3bpwAMP1C9/+UtNnTpVkvTPf/5TY8eO1Zo1azRq1KjcNBwAAABAu+jWy/C+yNy5c7V8+XI98MADeueddzRr1iwdc8wx+uCDDyRJjz32mPbdd189/vjjGjZsmPbZZx+de+65zCwBAAAAXQBhaS82btyo++67Tw899JAOP/xwDR8+XD/60Y80ZcoU3XfffZKkjz76SBs2bNBDDz2k+++/X4sWLdLKlSv17W9/O8etBwAAANBW3fo1S59n9erVSiQS2m+//TLub2hoUO/evSVJlmWpoaFB999/f+pxv/vd73TQQQfp/fffZ2keAAAA0IkRlvaiurpagUBAK1euVCAQyPhZjx49JEkDBgxQMBjMCFRf+tKXJDkzU4QlAAAAoPMiLO3FhAkTlEgktG3bNh1++OGejznssMMUj8f14Ycfavjw4ZKktWvXSpKGDh3aYW0FAAAAkH3d+t3wqqurtW7dOklOOLrllls0bdo0lZeXa8iQITr99NP1yiuv6Oabb9aECRO0fft2LV26VGPHjtXxxx8vy7J08MEHq0ePHrrttttkWZbmzJmj0tJSPfPMMzn+6wAAAAC0RbcOS8uWLdO0adOa3X/GGWdo0aJFisViuu6663T//ffrk08+UZ8+ffSVr3xFV199tcaMGSNJ+vTTT/XDH/5QzzzzjIqLi3Xsscfq5ptvVnl5eUf/OQAAAACyqFuHJQAAAADYG946HAAAAAA8EJYAAAAAwEO3C0u2bauqqkqsPgQAAADwebrdW4fv2bNHZWVl2r17t0pLS3PdHKBDbauq1w//9KbGDCpTcaTb/fMHALSTmoa4iiNBnXboEFWUFuS6OUDWdLuZJaA727anQa+v36maaDzXTQEAdCE10bgWLv1A2/Y05LopQFYRlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADzkNCwtWLBABx98sEpKSlRRUaGZM2fq/fff/9znLFq0SIZhZGwFBQUd1GIAAAAA3UVOw9ILL7ygOXPm6LXXXtOSJUsUi8X0ta99TTU1NZ/7vNLSUm3evDm1bdiwoYNaDAAAAKC7CObylz/11FMZtxctWqSKigqtXLlSRxxxxF6fZxiG+vfv397NAwDkQGEooHGVZSoKB2UYuW5NbsUTttbvqNHGnbW5bgoAdEt59Zql3bt3S5LKy8s/93HV1dUaOnSoKisr9Y1vfEPvvvtuRzQPANDOAqahg4b2UnGEoCRJwYChkf16qH8py80BIBfyJixZlqWLLrpIhx12mA488MC9Pm7UqFG699579be//U1/+MMfZFmWJk+erH//+9+ej29oaFBVVVXGBgDIT8XhgCKhvPmvKW+UF4dz3QQA6Jby5n+kOXPmaM2aNXrggQc+93GTJk3S7NmzNX78eB155JF65JFH1LdvX919992ej1+wYIHKyspSW2VlZXs0HwCQBQGze08n/WP5yxpX2UtV7kqLpO7eLwCQK3kRlubOnavHH39czz//vAYPHuzruaFQSBMmTNC6des8fz5//nzt3r07tW3atCkbTQYAQJ9s2qhxlb30r3dX57opAIB2kNOwZNu25s6dq7/+9a967rnnNGzYMN/7SCQSWr16tQYMGOD580gkotLS0owNAICOFItGc90EAEAr5DQszZkzR3/4wx+0ePFilZSUaMuWLdqyZYvq6upSj5k9e7bmz5+fun3NNdfomWee0UcffaRVq1bp9NNP14YNG3Tuuefm4k8AALSzV55/VmecdIymHDBUR4zZV3PPPFmbPl4vqXFm59m/P6ZzvnOCDh05ULO+NkVvr3wjYx9v/uM1nTPr6zp05EBNOXAf/eC0b6lq1y5Jzmtmf3f7LTp28jgdMmKAZn1tipY88bfUc6t27dL8H56nqeNG6JARA3TC4Qfp0Qf/KEk6bvI4SdLJxxyhcZW9dM6sr0uSfnrxf+uic07Tb351k6Yf9CWdOPVgSdJjDz+gU46bpkmjK3XUl0fpJ3PP1X92bG/X/gMAtF5O3zr8zjvvlCRNnTo14/777rtPZ555piRp48aNMs3GTPfZZ5/pvPPO05YtW9SrVy8ddNBBevXVV7X//vt3VLMBAB2orq5W3ztvjvYbfYBqa2v065t+oYvPO11/fvql1GNuv/E6zbvyGg0ZNly333idfjL3XD320ioFg0H9693V+q9TZmrmd07TZVcvUCAQ1D+Wv6SElZAk/e72W/TEXx/Slb+4RUOHDdfK11/V5Rd+X73K+2jipMN0+03X68MP3tcd9z+knuW9tenjj1Rf71zU++NjS3XaCUfrnj89quH7jVYo1PhGDK+/8qKKS0p01+JHUvfFY3HNufRy7bPvSO38z3bddM0V+tm8/9Yd9z/UQb0JAPAjp2HJtu0vfMyyZcsybt9666269dZb26lFAIB8M/24EzNuX33z7Zo6boQ+XPsvFRX3kCTN/v5cHXH0DEnS+fN+opOOnqRNH3+kYSP206I7F2r/seN1xS9uTu1jxKgvSZKiDQ367e236p4//VXjDjpEkjR46D568x+v6S9/vE8TJx2mLZ/8W6MPGKsDxk2QJA2qHJLaT6/efSRJZb3K1aeiX0Y7C4uK9PMbf6VQuDFAffO7p6e+Hzx0H/346l/q1K8fpdqa6tTfAgDIHzkNSwAAfJEN6z/Ur2/6hVa/tVK7du6UZVmSpC2f/lv7jhwtSdpv9AGpx/etcD60fOeO7Ro2Yj+9/94affX4b3jue+PHH6m+rlbfP/WkjPtjsahGHzBWkvSd2Wfrkv86Q/9a87YmHTFN02Ycr/ETD/3Cdo8cvX9GUJKk9955S3feeoPWvrdGVbt3p/6WzZ/8W8P3G92S7gAAdCDCEgAgr11w1ikaMKhSV/1yofr26y/LsvSt6ZMVi8ZSjwmGQo1PcD/N1nJXL0QK9v6BrrW1NZKk2xc9qIr+mW8UFI44QWfKtK/q76+9o5efW6LlLz2v//ruTJ18xrm65KfXfm67CwuLmv2u80//liYfeZQW/Ooe9erdR5s/+bfOP/1bisVie9kLACCX8uKtwwEA8LLrs536+MMP9F8XXKJDpxypfUeOUtXuXb72MXL0AXr9lRc8fzZ85CiFIxFt/nSThgzbN2PrP7DxoyzKe/fRibNO0YJf3aNLf/4LPbz495Kcj6+QJCuR+MJ2fLzuA+36bKcunH+VvnzoZA0bsZ92/oc3dwCAfMbMEgAgb5WW9VTPXuX6y+Lfq09Ff23+9N9auOBqX/s4Z+7F+vZXD9P1l1+iWd87S6FQWP949SV99esz1au8t874r7m66eorZFu2Jhz8FVXvqdKbK15Xjx4lOnHWKbrjpl9o/zHjNXy/0YpGG/Ti0qc1bMR+kqTyPn1VUFCoV5Y9q34DBiociaiktMyzHf0HDVYoHNaf7rtHs04/W+vef0/3LLypzX0EAGg/zCwBAPKWaZr65R2/0z9Xv61vfXWybrr6cs274hpf+9hn3xG66w+PaO0/1+i0E6brezO/puef+bsCAed64ZxLr9B/XXipfnfHrZp51KE6/3vf1ktLn0m9kUMoFNavfnmNZn1tis7+9vEKmAH98o7fSZKCwaB+fM0N+ssfF2n6xC/ponNO22s7ynv30bU336Fnnvibvnn0V3Tvr2/TvCv9/S0AgI5l2C15S7oupKqqSmVlZdq9ezcfUItuZ80nu/X1/31ZpxxSqYqSvb+OA8iVXkUhfXlor1w3I+9sq2rQ6k9257oZwF5t21OvP72xSY//cIoOHOQ9uwp0RswsAQAAAIAHwhIAAAAAeCAsAQDyRkPcynUT8lJD/IvfbQ8AkH2EJQBA3qiNJrRpZ22um5FX6qIJbfgPfQIAucBbhwMA8srardXaXRdXj0hAhvsBs91VNG5pa1U9M24AkCOEJQBA3tlaVa+tuW4EAKDbYxkeAAAAAHggLAEAAACAB8ISAAAAAHggLAEAAACAB8ISAAAAAHggLAEAAACAB8ISAAAAAHggLAEAAACAB8ISAAAAAHggLAEAAACAB8ISAAAAAHggLAEAAACAB8ISAAAAAHggLAEAAACAB8ISAAAAAHggLAEAAACAB8ISAAAAAHggLAEAAACAB8ISAAAAAHggLAEAAACAB8ISAAAAAHggLAEAAACAB8ISAAAAAHggLAEAAACAB8ISAAAAAHgI5roBAAAA6Fz6lkS0T+9imYa05pPduW4O0G4ISwAAAGixyl6F+vbESv2nukGF4YAmj+ijW5e8n+tmAe2CsAQAAIAWG9mvRP/cXKVjF74kSXpt/tG6+KujZNmSYeS4cUCW8ZolAAAAtNjm3XUa1a8kdfuom5fp7hc+1JxpI3TAwLIctgzIPsISAAAAWuyfm/fINA1dNH2kCkMBzTposG58+n0ddO0S3ffKetm2rYRl5bqZQFawDC+X9myR/nK2NHCCFO6R69agGxjWENefRu7W2KIyFQd9/vOPVjvjdOJZUkn/9mkgAKDTiMYtlReHdfzYgXrj48/0z81Vuvqx91RWGNKMA/rrg93vaMmGJSoOFee6qcgDNbEaFYeKNWu/Wepb1DfXzWkxZpZyac8WacMrUkN1rluCbqI4EtSkfXurONKK6yQN1dILNzjjFgDQrVU3xPXrZR/qk111+s7dy9UQT+iqE/bX/54yQdP376fiSFB/fv/Pqo3V5rqpyBO1sVrd+fad2l63PddN8YWZJQAAAPjy+1c/bvz+7EN0+Ig+Ms3Gd3f4ePfHeuyjx/Ttkd/OQeuA7CEsAQAAwJezDxum48YM0LA+3kvsauPMKKFrYBkeAAAAfCkMB/YalICuhLAEAAAAAB4ISwAAAADggbAEAAAAAB4ISwAAAADggbAEAAAAAB4ISwAAAADggc9ZAoDWsCzJikmJmGQnJNuSbDvzMYYhyZAMUzIDzlfDlIyAe9vw3DVayLadfrfi7lf3OMh2j4WdeUyM5LEISWbQ3bhm6Mm2nf604o3j20ooo0+T49sMOv0acPuVce0t2Z/JcZq+SY39lnGOCKSdO/KkX+t3S+8+KpXvKw0+WAoV5LpFQLsiLAFoOyshNexxCqVggRTI41OLZUnxOineICWiUrzeLQJ9Mgyn6A6EnILGq5ix7cxC07aaFJ6tYBhSICIFwk4/h4qkYKR1++pItu32eYMUjzYGTSvetv0aZmOxnh5IZTQ/HpYl2XHJqnGL1njzgNuS3xcqdPo93CM/w5ZlSbEaKVojxeoai3E/ksHSCDTv36Tk+I7XO7+jtcfTDErhYqdfg4X52adJiZgUrZaitc5Y9jN+kmM148KJ2dinthtEbTegWom0r604hlLjOTkYzt65onqbtP4Fad0Sad1SacJpUp+Rbd8vkKfyuKIB0KlYbtFUu6NtQcAMOcWEGXRDSFrBmyzOkkVEIua/iDAMpyALRpxit6hPfoe7pizLKdISUSkRdwqXRNT/fgwjrY/dQjg5C9a0z2WnFW6tDRmGG/AiTuFmFnfOmQgr4QSQWK0z1n0Vy0aTGYPgXkK2lTYLEfP/78kwpFCxM76L+7rHNo8lYk6wq6+S4ltbF2DN9P5scvEiOQOZvFjR2jEsucGuh1RU7oSQfB+7iZhzXo5HpZrtzgULv1LnZTfo1e1y9hOrk2q2ObeBLqwTVQgA8pYZkIp7t30/yaVtqeU/aUuqUle6Q43FZiCU/4VgtpmmZBY6V+ElSa3sdyu5fK3pkqAmy9dM98p3MNJY4Hfn5WtmQIr0cDZV+HuuZTWZMYg372/DkIyg099msPHiQb4X5W0RCEmFPZ3Nr9RSzLRzRtNlsabp9Gn6TFl3GcMBd/a7LRNKtu1ckEkuISzsKfUaKvUcKu13jPM90IURlgDkD9OUzE6wpKwrME3JDOe6Fd2LaUpyX9uD7EifrRPjuV0Yhrt8zz03F5RK06/pXDPyQBt0g8sqAAAAyBqCEroRwhIAAAAAeCAsAQAAAIAHwhIAAAAAeCAsAQAAAIAHwhIAAAAAeCAsAQAAAICHnIalBQsW6OCDD1ZJSYkqKio0c+ZMvf/++1/4vIceekijR49WQUGBxowZoyeffLIDWgsAAACgO8lpWHrhhRc0Z84cvfbaa1qyZIlisZi+9rWvqaamZq/PefXVV3XKKafonHPO0ZtvvqmZM2dq5syZWrNmTQe2HAAAAEBXl9NPFXvqqacybi9atEgVFRVauXKljjjiCM/nLFy4UMccc4wuvfRSSdK1116rJUuW6Pbbb9ddd93V7m0GAAAA0D3k1WuWdu/eLUkqLy/f62OWL1+u6dOnZ9w3Y8YMLV++vF3bBgAAAKB7yenMUjrLsnTRRRfpsMMO04EHHrjXx23ZskX9+vXLuK9fv37asmWL5+MbGhrU0NCQul1VVZWdBgMAAADo0vJmZmnOnDlas2aNHnjggazud8GCBSorK0ttlZWVWd0/AAAAgK4pL8LS3Llz9fjjj+v555/X4MGDP/ex/fv319atWzPu27p1q/r37+/5+Pnz52v37t2pbdOmTVlrNwAAAICuK6fL8Gzb1g9/+EP99a9/1bJlyzRs2LAvfM6kSZO0dOlSXXTRRan7lixZokmTJnk+PhKJKBKJZKvJAOCwLMmKSYmYZCck25JsO/MxhiHJkAxTMgPOV8OUjIB728hJ07sM23b63Yq7X93jINs9FnbmMTGSxyIkmUF3y4trhvnHtp3+tOKN49tKKKNPk+PbDDr9GnD7lXHtLdmfyXGavkmN/ZZxjgiknTvypF+thPTpm1LpQGcDurichqU5c+Zo8eLF+tvf/qaSkpLU647KyspUWFgoSZo9e7YGDRqkBQsWSJIuvPBCHXnkkbr55pt1/PHH64EHHtCKFSt0zz335OzvALo9KyE17HEKpWCBFMibl0M2Z1lSvE6KN0iJqBSvd4tAnwzDKboDIaeg8SpmbDuz0LStJoVnKxiGFIhIgbDTz6EiKdgJLgjZttvnDVI82hg0rXjb9muYjcV6eiCV0fx4WJZkxyWrxi1a480Dbkt+X6jQ6fdwj/wMW5YlxWqkaI0Uq2ssxv1IBksj0Lx/k5LjO17v/I7WHk8zKIWLnX4NFuZnnyYlYlK0WorWOmPZz/hJjtWMCydmY5/abhC13YBqJdK+tuIYSo3n5GA4e+eKnR9J//itc/4bepg0crpU3Kft+wXyVE4rmjvvvFOSNHXq1Iz777vvPp155pmSpI0bN8pMO3FOnjxZixcv1pVXXqnLL79cI0eO1KOPPvq5bwoBoANYbtFUu6NtQcAMOcWEGXRDSFrBmyzOkkVEIua/iDAMpyALRpxit6hPfoe7pizLKdISUSkRl6q3Od/7ZRhpfewWwslZsKZ9LjutcGttyDDcgBdxCjezuHPORFgJJ4DEap2x7qtYNprMGAT3ErKttFmImP9/T4YhhYqd8V3c1z22eSwRc4JdfZUU39q6AGum92eTixfJGcjkxYrWjmHJDXY9pKJyJ4Tk+9hNxJzzcjwq1Wx3Llj4lTovu0Gvbpdzcaxhj9RQJZUNkoYdnvWmA/ki58vwvsiyZcua3Tdr1izNmjWrHVoEoFXMgFTcu+37SS5tSy3/SVtSlbrSHWosNgOh/C8Es800JbPQuQovSWplv1vJ5WtNlwQ1Wb5mule+g5HGAr87L18zA1Kkh7Opwt9zLavJjEG8eX8bhmQEnf42g40XD/K9KG+LQEgq7OlsfqWWYqadM5ouizVNp0/TZ8q6yxgOuLPfbZlQsm3ngkxyCWEgLEVKpJIB0qjjpEEHZa25QD7qRJdTAXR5pimZnWBJWVdgmpIZznUruhfTlOS+tgfZkT5bJ8ZzuzAMd/mee24eOE6aOt+ZtQwX5bRpQEcgLAEAAKDleg3NdQuADtMN5qABAAAAwD/CEgAAAAB4ICwBAAAAgAfCEgAAAAB4ICwBAAAAgAfCEgAAAAB4ICwBAAAAgAfCEgAAAAB4ICwBAAAAgAfCEgAAAAB4ICwBAAAAgIeg3yesX79eL730kjZs2KDa2lr17dtXEyZM0KRJk1RQUNAebQSA/GNZkhWTEjHJTki2Jdl25mMMQ5IhGaZkBpyvhikZAfe2kZOmdxm27fS7FXe/usdBtnss7MxjYiSPRUgyg+7GNUNPtu30pxVvHN9WQhl9mhzfZtDp14Dbr4xrb8n+TI7T9E1q7LeMc0Qg7dyRJ/1av1t691GpfF9p8MFSiNoPXVuLw9If//hHLVy4UCtWrFC/fv00cOBAFRYWaufOnfrwww9VUFCg0047TT/+8Y81dOjQ9mwzgHxjJaSGPU6hFCyQAr6vw3Qcy5LidVK8QUpEpXi9WwT6ZBhO0R0IOQWNVzFj25mFpm01KTxbwTCkQEQKhJ1+DhVJwUjr9tWRbNvt8wYpHm0Mmla8bfs1zMZiPT2Qymh+PCxLsuOSVeMWrfHmAbclvy9U6PR7uEd+hi3LkmI1UrRGitU1FuN+JIOlEWjev0nJ8R2vd35Ha4+nGZTCxU6/Bgvzs0+TEjEpWi1Fa52x7Gf8JMdqxoUTs7FPbTeI2m5AtRJpX1txDKXGc3IwnL1zRfU2af0L0rol0rql0oTTpD4j275fIE+1qKKZMGGCwuGwzjzzTD388MOqrKzM+HlDQ4OWL1+uBx54QBMnTtSvf/1rzZo1q10aDCBPWW7RVLujbUHADDnFhBl0Q0hawZsszpJFRCLmv4gwDKcgC0acYreoT36Hu6YsyynSElEpEXcKl0TU/34MI62P3UI4OQvWtM9lpxVurQ0ZhhvwIk7hZhZ3zpkIK+EEkFitM9Z9FctGkxmD4F5CtpU2CxHz/+/JMKRQsTO+i/u6xzaPJWJOsKuvkuJbWxdgzfT+bHLxIjkDmbxY0doxLLnBrodUVO6EkHwfu4mYc16OR6Wa7c4FC79S52U36NXtcvYTq5Nqtjm3gS6sRRXCDTfcoBkzZuz155FIRFOnTtXUqVN1/fXX6+OPP85W+wB0BmZAKu7d9v0kl7allv+kLalKXekONRabgVD+F4LZZpqSWehchZcktbLfreTytaZLgposXzPdK9/BSGOB352Xr5kBKdLD2VTh77mW1WTGIN68vw1DMoJOf5vBxosH+V6Ut0UgJBX2dDa/Uksx084ZTZfFmqbTp+kzZd1lDAfc2e+2TCjZtnNBJrmEsLCn1Guo1HOotN8xzvdAF9aisPR5Qamp3r17q3fvLBRNALof05TMTrCkrCswTckM57oV3YtpSnJf24PsSJ+tE+O5XRiGu3zPPTcXlErTr+lcM/JAG/i+rLJq1SqtXr06dftvf/ubZs6cqcsvv1zRaCuWggAAAKDzICihG/Edlr7//e9r7dq1kqSPPvpI3/3ud1VUVKSHHnpIl112WdYbCAAAAAC54DssrV27VuPHj5ckPfTQQzriiCO0ePFiLVq0SA8//HC22wcAAAAAOeE7LNm2Lcty3n3q2Wef1XHHHSdJqqys1I4dO7LbOgAAAADIEd9haeLEibruuuv0f//3f3rhhRd0/PHHS3I+rLZfv35ZbyAAAAAA5ILvsHTbbbdp1apVmjt3rq644gqNGDFCkvSXv/xFkydPznoDAQAAACAXfL+dydixYzPeDS/pf/7nfxQIdLPPOwEAAADQZbX6vR+j0ai2bduWev1S0pAhQ9rcKAAAAADINd9hae3atTrnnHP06quvZtxv27YMw1Aikcha4wAAAAAgV3yHpbPOOkvBYFCPP/64BgwYIMMw2qNdAAAAAJBTvsPSW2+9pZUrV2r06NHt0R4AAAAAyAu+3w1v//335/OUAAAAAHR5vsPSL3/5S1122WVatmyZ/vOf/6iqqipjAwAAAICuwPcyvOnTp0uSjj766Iz7eYMHAAAAAF2J77D0/PPPt0c7AAAAACCv+A5LRx55ZHu0AwAAAADySqs+lHbXrl363e9+p3/+85+SpAMOOEBnn322ysrKsto4AAAAAMgV32/wsGLFCg0fPly33nqrdu7cqZ07d+qWW27R8OHDtWrVqvZoIwAAAAB0ON8zSxdffLFOPPFE/eY3v1Ew6Dw9Ho/r3HPP1UUXXaQXX3wx640EAAAAgI7mOyytWLEiIyhJUjAY1GWXXaaJEydmtXEAAAAAkCu+l+GVlpZq48aNze7ftGmTSkpKstIoAAAAAMg132Hp5JNP1jnnnKMHH3xQmzZt0qZNm/TAAw/o3HPP1SmnnNIebQQAAACADud7Gd5NN90kwzA0e/ZsxeNxSVIoFNL555+vG264IesNBAAAAIBc8B2WwuGwFi5cqAULFujDDz+UJA0fPlxFRUVZbxwAAAAA5EqrPmdJkoqKijRmzJhstgUAAAAA8kaLwtJJJ52kRYsWqbS0VCeddNLnPvaRRx7JSsMAAAAAIJdaFJbKyspkGIYk593wkt8DAAAAQFfVorB03333pb5ftGhRe7UFAAAAAPKG77cOP+qoo7Rr165m91dVVemoo47KRpsAAAAAIOd8h6Vly5YpGo02u7++vl4vvfRSVhoFAAAAALnW4nfDe+edd1Lfv/fee9qyZUvqdiKR0FNPPaVBgwZlt3UAAAAAkCMtDkvjx4+XYRgyDMNzuV1hYaH+93//N6uNA4C8ZVmSFZMSMclOSLYl2XbmYwxDkiEZpmQGnK+GKRkB9zZvltMmtu30uxV3v7rHQbZ7LOzMY2Ikj0VIMoPu5nuBRfdg205/WvHG8W0llNGnyfFtBp1+Dbj9yrj2luzP5DhN36TGfss4RwTSzh150q8Ne6R/Pi712kca9GUpGMl1i4B21eKwtH79etm2rX333VdvvPGG+vbtm/pZOBxWRUWFAoFAuzQSQJ6zEs5/oGZQChZIgVZ/hFv7sywpXifFG6REVIrXu0WgT4bhFN2BkFPQeBUztp1ZaNpWk8KzFQxDCkSkQNjp51BR5yhWbNvt8wYpHm0Mmla8bfs1zMZiPT2Qymh+PCxLsuOSVeMWrfHmAbclvy9U6PR7uEd+hi3LkmI1UrRGitU1FuN+JIOlEWjev0nJ8R2vd35Ha4+nGZTCxU6/Bgvzs0+TEjEpWi1Fa52x7Gf8JMdqxoUTs7FPbTeI2m5AtRJpX1txDKXGc3IwnL1zxZ4t0rolTl/02U8af6rUe3jb9wvkqRZXNEOHDpUkWVYr/8EC6Nost2iq3dG2IGCGnGLCDLohJK3gTRZnySIiEfNfRBiGU5AFI06xW9Qnv8NdU5blFGmJqJSIS9XbnO/9Moy0PnYL4eQsWNM+l51WuLU2ZBhuwIs4hZtZ3DlnIqyEE0Bitc5Y91UsG01mDIJ7CdlW2ixEzP+/J8OQQsXO+C7u6x7bPJaIOcGuvkqKb21dgDXT+7PJxYvkDGTyYkVrx7DkBrseUlG5E0LyfewmYs55OR6VarY7Fyz8Sp2X3aBXt8vZT6xWqvpUqtlBWEKX1uoK4b333tPGjRubvdnDiSee2OZGAehkzIBU3Lvt+0kubUst/0lbUpW60h1qLDYDofwvBLPNNCWz0LkKL0lqZb9byeVrTZcENVm+ZrpXvoORxgK/Oy9fMwNSpIezqcLfcy2ryYxBvHl/G4ZkBJ3+NoONFw/yvShvi0BIKuzpbH6llmKmnTOaLos1TadP02fKussYDriz322ZULJt54JMcglhQalUNlgqq5RGHSeVD8tac4F85DssffTRR/rmN7+p1atXyzAM2e7JKPlBtYlEK68oA4BpSmYnWFLWFZimZIZz3YruxTQlua/tQXakz9aJ8dwuDMNdvueemwtKpa9e6y5B7sIhHnD5vqxy4YUXatiwYdq2bZuKior07rvv6sUXX9TEiRO1bNmydmgiAAAA8kYwTFBCt+F7Zmn58uV67rnn1KdPH5mmKdM0NWXKFC1YsEAXXHCB3nzzzfZoJwAAAAB0KN8zS4lEQiUlJZKkPn366NNPP5XkvAHE+++/n93WAQAAAECO+J5ZOvDAA/X2229r2LBhOvTQQ3XjjTcqHA7rnnvu0b777tsebQQAAACADuc7LF155ZWqqamRJF1zzTX6+te/rsMPP1y9e/fWgw8+mPUGAgAAAEAu+A5LM2bMSH0/YsQI/etf/9LOnTvVq1ev1DviAQAAAEBn5+s1S7FYTMFgUGvWrMm4v7y8nKAEAAAAoEvxFZZCoZCGDBnCZykBAAAA6PJ8vxveFVdcocsvv1w7d+5sj/YAAAAAQF7w/Zql22+/XevWrdPAgQM1dOhQFRcXZ/x81apVWWscAAAAAOSK77A0c+bMdmgGAAAAAOQX32Hpqquuao92AAAAAEBe8f2aJUnatWuXfvvb32r+/Pmp1y6tWrVKn3zyia/9vPjiizrhhBM0cOBAGYahRx999HMfv2zZMhmG0WzbsmVLa/4MAAAAANgr3zNL77zzjqZPn66ysjJ9/PHHOu+881ReXq5HHnlEGzdu1P3339/ifdXU1GjcuHE6++yzddJJJ7X4ee+//75KS0tTtysqKnz9DQAAAADwRXyHpXnz5unMM8/UjTfeqJKSktT9xx13nE499VRf+zr22GN17LHH+m2CKioq1LNnT9/PAwAAAICW8r0M7x//+Ie+//3vN7t/0KBBHbYcbvz48RowYIC++tWv6pVXXvncxzY0NKiqqipjAwAAAIAv4jssRSIRz8Cxdu1a9e3bNyuN2psBAwborrvu0sMPP6yHH35YlZWVmjp16ue+XfmCBQtUVlaW2iorK9u1jQAAAAC6Bt9h6cQTT9Q111yjWCwmSTIMQxs3btSPf/xjfetb38p6A9ONGjVK3//+93XQQQdp8uTJuvfeezV58mTdeuute33O/PnztXv37tS2adOmdm0jAAAAgK7Bd1i6+eabVV1drYqKCtXV1enII4/UiBEjVFJSouuvv7492vi5DjnkEK1bt26vP49EIiotLc3YAAAAAOCL+H6Dh7KyMi1ZskQvv/yy3nnnHVVXV+vLX/6ypk+f3h7t+0JvvfWWBgwYkJPfDQAAAKDr8h2WkqZMmaIpU6a06ZdXV1dnzAqtX79eb731lsrLyzVkyBDNnz9fn3zySertyG+77TYNGzZMBxxwgOrr6/Xb3/5Wzz33nJ555pk2tQMAAAAAmmpVWFq6dKmWLl2qbdu2ybKsjJ/de++9Ld7PihUrNG3atNTtefPmSZLOOOMMLVq0SJs3b9bGjRtTP49Go7rkkkv0ySefqKioSGPHjtWzzz6bsQ8AAAAAyAbfYenqq6/WNddco4kTJ2rAgAEyDKPVv3zq1KmybXuvP1+0aFHG7csuu0yXXXZZq38fAAAAALSU77B01113adGiRfre977XHu0BAAAAgLzg+93wotGoJk+e3B5tAQAAAIC84TssnXvuuVq8eHF7tAUAAAAA8obvZXj19fW655579Oyzz2rs2LEKhUIZP7/llluy1jgAAAAAyBXfYemdd97R+PHjJUlr1qzJ+Flb3uwBAAAAAPKJ77D0/PPPt0c7AAAAACCv+H7NUtK6dev09NNPq66uTpI+9y3AAQAAAKCz8R2W/vOf/+joo4/Wfvvtp+OOO06bN2+WJJ1zzjm65JJLst5AAAAAAMgF32Hp4osvVigU0saNG1VUVJS6/+STT9ZTTz2V1cYBAAAAQK74fs3SM888o6efflqDBw/OuH/kyJHasGFD1hoGAAAAALnke2appqYmY0YpaefOnYpEIllpFAAAAADkmu+wdPjhh+v+++9P3TYMQ5Zl6cYbb9S0adOy2jgAAAAAyBXfy/BuvPFGHX300VqxYoWi0aguu+wyvfvuu9q5c6deeeWV9mgjAAAAAHQ43zNLBx54oNauXaspU6boG9/4hmpqanTSSSfpzTff1PDhw9ujjQAAAADQ4XzPLElSWVmZrrjiimy3BQAAAADyhu+w9M4773jebxiGCgoKNGTIEN7oAQAAAECn5zssjR8/XoZhSJJs25ak1G1JCoVCOvnkk3X33XeroKAgS80EAAAAgI7l+zVLf/3rXzVy5Ejdc889evvtt/X222/rnnvu0ahRo7R48WL97ne/03PPPacrr7yyPdoLAAAAAB3C98zS9ddfr4ULF2rGjBmp+8aMGaPBgwfrpz/9qd544w0VFxfrkksu0U033ZTVxgIAAABAR/E9s7R69WoNHTq02f1Dhw7V6tWrJTlL9TZv3tz21gEAAABAjvgOS6NHj9YNN9ygaDSaui8Wi+mGG27Q6NGjJUmffPKJ+vXrl71WAgAAAEAH870M74477tCJJ56owYMHa+zYsZKc2aZEIqHHH39ckvTRRx/pv//7v7PbUgAAAADoQL7D0uTJk7V+/Xr98Y9/1Nq1ayVJs2bN0qmnnqqSkhJJ0ve+973sthIAAAAAOlirPpS2pKREP/jBD7LdFgAAAADIGy16zdJrr73W4h3W1tbq3XffbXWDAAAAACAftCgsfe9739OMGTP00EMPqaamxvMx7733ni6//HINHz5cK1euzGojAQAAAKCjtWgZ3nvvvac777xTV155pU499VTtt99+GjhwoAoKCvTZZ5/pX//6l6qrq/XNb35TzzzzjMaMGdPe7QYAAACAdtWisBQKhXTBBRfoggsu0IoVK/Tyyy9rw4YNqqur07hx43TxxRdr2rRpKi8vb+/2AgAAAECH8P0GDxMnTtTEiRPboy0AAAAAkDd8fygtAAAAAHQHhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPrQpLL7zwgk444QSNGDFCI0aM0IknnqiXXnop220DAAAAgJzxHZb+8Ic/aPr06SoqKkp99lJhYaGOPvpoLV68uD3aCAAAAAAdzvfnLF1//fW68cYbdfHFF6fuu+CCC3TLLbfo2muv1amnnprVBgIAAABALvieWfroo490wgknNLv/xBNP1Pr167PSKAAAAADINd9hqbKyUkuXLm12/7PPPqvKysqsNAoAAAAAcs33MrxLLrlEF1xwgd566y1NnjxZkvTKK69o0aJFWrhwYdYbCAAAAAC54DssnX/++erfv79uvvlm/fnPf5YkfelLX9KDDz6ob3zjG1lvIAAAAADkgu+wJEnf/OY39c1vfjPbbQEAAACAvMGH0gIAAACAB98zS7169ZJhGM3uNwxDBQUFGjFihM4880ydddZZWWkgAAAAAOSC77D0s5/9TNdff72OPfZYHXLIIZKkN954Q0899ZTmzJmj9evX6/zzz1c8Htd5552X9QYDAAAAQEfwHZZefvllXXfddfrBD36Qcf/dd9+tZ555Rg8//LDGjh2rX/3qV4QlAAAAAJ2W79csPf3005o+fXqz+48++mg9/fTTkqTjjjtOH330UdtbBwAAAAA54ntmqby8XI899pguvvjijPsfe+wxlZeXS5JqampUUlKSnRYCQD6yLMmKSYmYZCck25JsO/MxhiHJkAxTMgPOV8OUjIB7u/nrP+GDbTv9bsXdr+5xkO0eCzvzmBjJYxGSzKC78T5Hnmzb6U8r3ji+rYQy+jQ5vs2g068Bt18Z196S/Zkcp+mb1NhvGeeIQNq5I0/6tWGP9P7fpV77SAMnOMcd6MJ8h6Wf/vSnOv/88/X888+nXrP0j3/8Q08++aTuuusuSdKSJUt05JFHZrelAPKXlXD+AzWDUrBACrTqUwk6hmVJ8Top3iAlolK83i0CfTIMp+gOhJyCxquYse3MQtO2mhSerWAYUiAiBcJOP4eKpGCkdfvqSLbt9nmDFI82Bk0r3rb9GmZjsZ4eSGU0Px6WJdlxyapxi9Z484Dbkt8XKnT6PdwjP8OWZUmxGilaI8XqGotxP5LB0gg079+k5PiO1zu/o7XH0wxK4WKnX4OF+dmnSYmYFK2WorXOWPYzfpJjNePCidnYp7YbRG03oFqJtK+tOIZS4zk5GM7euWLPFmntU067+o6Wxn5H6j287fsF8pTviua8887T/vvvr9tvv12PPPKIJGnUqFF64YUXNHnyZEnSJZdckt1WAsh/lls01e5oWxAwQ04xYQbdEJJW8CaLs2QRkYj5LyIMwynIghGn2C3qk9/hrinLcoq0RFRKxKXqbc73fhlGWh+7hXByFqxpn8tOK9xaGzIMN+BFnMLNLO6cMxFWwgkgsVpnrPsqlo0mMwbBvYRsK20WIub/35NhSKFiZ3wX93WPbR5LxJxgV18lxbe2LsCa6f3Z5OJFcgYyebGitWNYcoNdD6mo3Akh+T52EzHnvByPSjXbnQsWfqXOy27Qq9vl7CdaK332sVS9nbCELq1VFcJhhx2mww47LNttAdBZmQGpuHfb95Nc2pZa/pO2pCp1pTvUWGwGQvlfCGabaUpmoXMVXpLUyn63ksvXmi4JarJ8zXSvfAcjjQV+d16+ZgakSA9nU4W/51pWkxmDePP+NgzJCDr9bQYbLx7ke1HeFoGQVNjT2fxKLcVMO2c0XRZrmk6fps+UdZcxHHBnv9syoWTbzgWZ5BLCSKlUOlAqq5RGHSv1HpG15gL5qE2XU+vr6xWNZl7RLC0tbVODAHRjpimZnWBJWVdgmpIZznUruhfTlGTyGo9sSp+tE+O5XRiGu3zPPTcXlEpfvda9gNKFQzzg8n1Zpba2VnPnzlVFRYWKi4vVq1evjA0AAABdWKgTLEEEssR3WLr00kv13HPP6c4771QkEtFvf/tbXX311Ro4cKDuv//+9mgjAAAAAHQ438vwHnvsMd1///2aOnWqzjrrLB1++OEaMWKEhg4dqj/+8Y867bTT2qOdAAAAANChfM8s7dy5U/vuu68k5/VJO3fulCRNmTJFL774YnZbBwAAAAA54jss7bvvvlq/fr0kafTo0frzn/8syZlx6tmzZ1YbBwAAAAC54jssnXXWWXr77bclST/5yU90xx13qKCgQBdffLEuvfTSrDcQAAAAAHLB92uWLr744tT306dP17/+9S+tXLlSI0aM0NixY7PaOAAAAADIlTZ/bP3QoUM1dOjQbLQFAAAAAPJGN/j4agAAAADwj7AEAAAAAB4ISwAAAADggbAEAAAAAB5aFZY+/PBDXXnllTrllFO0bds2SdLf//53vfvuu1ltHAAAAADkiu+w9MILL2jMmDF6/fXX9cgjj6i6ulqS9Pbbb+uqq67yta8XX3xRJ5xwggYOHCjDMPToo49+4XOWLVumL3/5y4pEIhoxYoQWLVrk908AAAAAgC/kOyz95Cc/0XXXXaclS5YoHA6n7j/qqKP02muv+dpXTU2Nxo0bpzvuuKNFj1+/fr2OP/54TZs2TW+99ZYuuuginXvuuXr66ad9/V4AAAAA+CK+P2dp9erVWrx4cbP7KyoqtGPHDl/7OvbYY3Xssce2+PF33XWXhg0bpptvvlmS9KUvfUkvv/yybr31Vs2YMcPX7wYAAACAz+N7Zqlnz57avHlzs/vffPNNDRo0KCuN2pvly5dr+vTpGffNmDFDy5cv3+tzGhoaVFVVlbEBAAAAwBfxHZa++93v6sc//rG2bNkiwzBkWZZeeeUV/ehHP9Ls2bPbo40pW7ZsUb9+/TLu69evn6qqqlRXV+f5nAULFqisrCy1VVZWtmsbAQAAAHQNvsPSL37xC40ePVqVlZWqrq7W/vvvryOOOEKTJ0/WlVde2R5tbJP58+dr9+7dqW3Tpk25bhIAAACATsD3a5bC4bB+85vf6Kc//anWrFmj6upqTZgwQSNHjmyP9mXo37+/tm7dmnHf1q1bVVpaqsLCQs/nRCIRRSKRdm8bAAAAgK7Fd1hKGjJkiIYMGZLNtnyhSZMm6cknn8y4b8mSJZo0aVKHtgMAAABA1+c7LNm2rb/85S96/vnntW3bNlmWlfHzRx55pMX7qq6u1rp161K3169fr7feekvl5eUaMmSI5s+fr08++UT333+/JOkHP/iBbr/9dl122WU6++yz9dxzz+nPf/6znnjiCb9/BgAAAAB8Lt9h6aKLLtLdd9+tadOmqV+/fjIMo9W/fMWKFZo2bVrq9rx58yRJZ5xxhhYtWqTNmzdr48aNqZ8PGzZMTzzxhC6++GItXLhQgwcP1m9/+1veNhwAAABA1vkOS//3f/+nRx55RMcdd1ybf/nUqVNl2/Zef75o0SLP57z55ptt/t0AAAAA8Hl8vxteWVmZ9t133/ZoCwAAAADkDd9h6ec//7muvvrqvX6uEQAAAAB0Bb6X4X3nO9/Rn/70J1VUVGifffZRKBTK+PmqVauy1jgAAAAAyBXfYemMM87QypUrdfrpp7f5DR4AAAAAIF/5DktPPPGEnn76aU2ZMqU92gMAAAAAecH3a5YqKytVWlraHm0BAAAAgLzhOyzdfPPNuuyyy/Txxx+3Q3MAAAAAID/4XoZ3+umnq7a2VsOHD1dRUVGzN3jYuXNn1hoHAAAAALniOyzddttt7dAMAAAAAMgvrXo3PAAAAADo6loUlqqqqlJv6lBVVfW5j+XNHwAAAAB0BS0KS7169dLmzZtVUVGhnj17en62km3bMgxDiUQi640EAAAAgI7WorD03HPPqby8XJL0/PPPt2uDAAAAACAftCgsHXnkkanvhw0bpsrKymazS7Zta9OmTdltHQAAAADkiO83eBg2bFhqSV66nTt3atiwYSzDA9A9WJZkxaRETLITkm1Jtp35GMOQZEiGKZkB56thSkbAvd18STN8sG2n3624+9U9DrLdY2FnHhMjeSxCkhl0N98fN9g92LbTn1a8cXxbCWX0aXJ8m0GnXwNuvzKuvSX7MzlO0zepsd8yzhGBtHNHnvRrIi59ukoqHSiVDc51a4B25zssJV+b1FR1dbUKCgqy0igAnYyVkBr2OIVSsEAK+D61dBzLkuJ1UrxBSkSleL1bBPpkGE7RHQg5BY1XMWPbmYWmbTUpPFvBMKRARAqEnX4OFUnBSOv21ZFs2+3zBikebQyaVrxt+zXMxmI9PZDKaH48LEuy45JV4xat8eYBtyW/L1To9Hu4R36GLcuSYjVStEaK1TUW434kg6URaN6/ScnxHa93fkdrj6cZlMLFTr8GC/OzT5MSMSlaLUVrnbHsZ/wkx2rGhROzsU9tN4jabkC1EmlfW3EMpcZzcjCcvXPFZ+ulFfc6+x4ySRr5NalH37bvF8hTLa5o5s2bJ0kyDEM//elPVVRUlPpZIpHQ66+/rvHjx2e9gQA6Ccstmmp3tC0ImCGnmDCDbghJK3iTxVmyiEjE/BcRhuEUZMGIU+wW9cnvcNeUZTlFWiLqXOGt3uZ875dhpPWxWwgnZ8Ga9rnstMKttSHDcANexCnczOLOORNhJZwAEqt1xrqvYtloMmMQ3EvIttJmIWL+/z0ZhhQqdsZ3cV/32OaxRMwJdvVVUnxr6wKsmd6fTS5eJGcgkxcrWjuGJTfY9ZCKyp0Qku9jNxFzzsvxqFSz3blg4VfqvOwGvbpdzsWx+irna88hhCV0aS2uEN58801JzszS6tWrFQ6HUz8Lh8MaN26cfvSjH2W/hQDynxmQinu3fT/JpW2p5T9pS6pSV7pDjcVmIJT/hWC2maZkFjpX4SVJrex3K7l8remSoCbL10z3yncw0ljgd+fla2ZAivRwNlV84cMzWFaTGYN48/42DMkIOv1tBhsvHuR7Ud4WgZBU2NPZ/EotxUw7ZzRdFmuaTp+mz5R1lzEccGe/2zKhZNvOBZnkEsJgRIqUSKWDpFHHSYMnZq25QD5qcVhKvgveWWedpYULF/J5SgCyzzQlsxMsKesKTFMyw1/8OGSPaUpyX9uD7EifrRPjuV0Yhrt8zz03DxgrTbtCKu6TdtEG6Lp8rz2577772qMdAAAA6Ax6Vua6BUCH6QZz0AAAAADgH2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAAADAA2EJAAAAADwQlgAAANBysbrGD/0Fujjfn7MEAJBkWZIVkxIxyU5IttW8eDAMSYZkmM6HZhqmu7kfomkYOWl6l2HbTr9bcferexxku8fCzjwmRvJYhCQz6G5cM/Rk205/WvHG8W0llNGnyfFtBp1+Dbj9yrj2luzP5DhN36TGfss4RwTSzh150q9Vm6Xlt0tlldKIo6XeI/KnbUA7ICwBaDsrITXscQqlYIEUyONTi2VJ8Top3iAlolK83i0CfTIMp+gOhJyCxquYse3MQtO2mhSerWAYUiAiBcJOP4eKpGCkdfvqSLbt9nmDFI82Bk0r3rb9GmZjsZ4eSGU0Px6WJdlxyapxi9a4/6vjhimFCp1+D/fIz7BlWVKsRorWuDMAlv99JIOlEWjev0nJ8R2vd35Ha4+nGZTCxU6/Bgvzs0+TEjEpWi1Fa52x7Gf8JMdqxoUTs7FPbTeI2m5AtRJpX1txDKXGc3IwnL1zRbTaCUw71kqfvikddKY0dFLb9wvkqTyuaAB0KpZbNNXuaFsQMENOMWEG3RCSVvAmi7NkEZGI+S8iDMMpyIIRp9gt6pPf4a4py3KKtERUSsSl6m3O934ZRlofu4VwchasaZ/LTivcWhsyDDfgRZzCzSzunDMRVsIJILFaZ6z7KpaNJjMGwb2EbCttFiLm/9+TYUihYmd8F/d1j20eS8ScYFdfJcW3ti7Amun92eTiRXIGMnmxorVjWHKDXQ+pqNwJIfk+dhMx57wcj0o1250LFn6lzstu0Kvb5ZxzzKDUax+ppF+2Ww3klU5UIQDIW2ZAKu7d9v0kl7allv+kLalKXekONRabgVD+F4LZZpqSWehchZcktbLfreTytaZLgposXzPdK9/BSGOB352Xr5kBKdLD2VTh77mW1WTGIN68vw1DMoJOf5vBxosH+V6Ut0UgJBX2dDa/Uksx084ZTZfFmqbTp+kzZd1lDAfc2e+2TCjZthOOkksISwdK+31NKhssDT2sc8xsA21AWAKQP0xTMvmPt0OYpmSGc92K7sU0Jbmv7UF2pM/WifHcLgzDDUTuubmgVPry7Jw2CehI3eCyCgAAAAD4R1gCAAAAAA+EJQAAAADwQFgCAAAAAA+EJQAAAADwQFgCAAAAAA+EJQAAAADwQFgCAAAAAA+EJQAAAADwQFgCAAAAAA+EJQAAAADwQFgCAAAAAA+EJQAAAADwQFgCAAAAAA+EJQAAAADwQFgCAAAAAA+EJQAAAADwQFgCAAAAAA+EJQAAAADwQFgCAAAAAA+EJQAAAADwQFgCAAAAAA+EJQAAALRctFay7Vy3AugQwVw3AAA6JcuSrJiUiEl2QrKt5sWDYUgyJMOUzIDz1TAlI+DeNnLS9C7Dtp1+t+LuV/c4yHaPhZ15TIzksQhJZtDduGboybad/rTijePbSiijT5Pj2ww6/Rpw+5Vx7S3Zn8lxmr5Jjf2WcY4IpJ078qRfqz6Vlt8ulVVKw4+W+ozMn7YB7YCwBKDtrITUsMcplIIFUiCPTy2WJcXrpHiDlIhK8Xq3CPTJMJyiOxByChqvYsa2MwtN22pSeLaCYUiBiBQIO/0cKpKCkdbtqyPZttvnDVI82hg0rXjb9muYjcV6eiCV0fx4WJZkxyWrxi1a4/6vjhumFCp0+j3cIz/DlmVJsRopWiPF6hqLcT+SwdIINO/fpOT4jtc7v6O1x9MMSuFip1+DhfnZp0mJmBStdmZWEg3+xk9yrGZcODEb+9R2g6jtBlQrkfa1FcdQajwnB8PZO1dEa6SqLdKOD6RP35IOOlMaOqnt+wXyVB5XNAA6Fcstmmp3tC0ImCGnmDCDbghJK3iTxVmyiEjE/BcRhuEUZMGIU+wW9cnvcNeUZTlFWiIqJeJS9Tbne78MI62P3UI4OQvWtM9lpxVurQ0ZhhvwIk7hZhZ3zpkIK+EEkFitM9Z9FctGkxmD4F5CtpU2CxHz/+/JMKRQsTO+i/u6xzaPJWJOAV5fJcW3ti7Amun92eTiRXIGMnmxorVjWHKDXQ+pqNwJIfk+dhMx57wcj0o1250LFn6lzstu0Kvb5ZxzzJBUPkwq6Z/1ZgP5pBNVCADylhmQinu3fT/JpW2p5T9pS6pSV7pDjcVmIJT/hWC2maZkFjpX4SVJrex3K7l8remSoCbL10z3yncw0ljgd+fla2ZAivRwNlX4e65lNZkxiDfvb8OQjKDT32aw8eJBvhflbREISYU9nc2v1FLMtHNG02Wxpun0afpMWXcZwwF39rstE0q27YSj5BLC0kHSqGOcr0Mnd46ZbaANCEsA8odpSib/8XYI05TMcK5b0b2YpiT3tT3IjvTZOjGe24VhuIHIPTcXlEoTTs9pk4CO1A0uqwAAAACAf4QlAAAAAPBAWAIAAAAAD4QlAAAAAPBAWAIAAAAAD3kRlu644w7ts88+Kigo0KGHHqo33nhjr49dtGiRDMPI2AoKCjqwtQAAAAC6g5yHpQcffFDz5s3TVVddpVWrVmncuHGaMWOGtm3bttfnlJaWavPmzaltw4YNHdhiAAAAAN1BzsPSLbfcovPOO09nnXWW9t9/f911110qKirSvffeu9fnGIah/v37p7Z+/fp1YIsBAAAAdAc5DUvRaFQrV67U9OnTU/eZpqnp06dr+fLle31edXW1hg4dqsrKSn3jG9/Qu+++2xHNBQAAANCN5DQs7dixQ4lEotnMUL9+/bRlyxbP54waNUr33nuv/va3v+kPf/iDLMvS5MmT9e9//9vz8Q0NDaqqqsrYAAAAAOCL5HwZnl+TJk3S7NmzNX78eB155JF65JFH1LdvX919992ej1+wYIHKyspSW2VlZQe3GAAAAEBnlNOw1KdPHwUCAW3dujXj/q1bt6p///4t2kcoFNKECRO0bt06z5/Pnz9fu3fvTm2bNm1qc7sBAAAAdH05DUvhcFgHHXSQli5dmrrPsiwtXbpUkyZNatE+EomEVq9erQEDBnj+PBKJqLS0NGMDAAAAgC8SzHUD5s2bpzPOOEMTJ07UIYccottuu001NTU666yzJEmzZ8/WoEGDtGDBAknSNddco6985SsaMWKEdu3apf/5n//Rhg0bdO655+byzwAAAADQxeQ8LJ188snavn27fvazn2nLli0aP368nnrqqdSbPmzcuFGm2TgB9tlnn+m8887Tli1b1KtXLx100EF69dVXtf/+++fqTwAAAADQBeU8LEnS3LlzNXfuXM+fLVu2LOP2rbfeqltvvbUDWgUAAACgO+t074YHAAAAAB2BsAQAAAAAHghLAAAAAOCBsAQAAAAAHghLAAAAAOCBsAQAAAAAHghLAAAAAOCBsAQAAAAAHghLAAAAAOCBsAQAAAAAHghLAAAAAOCBsAQAAAAAHghLAAAAAOCBsAQAAAAAHghLAAAAAOCBsAQAAAAAHghLAAAAAOAhmOsGAECnZFmSFZMSMclOSLYl2XbmYwxDkiEZpmQGnK+GKRkB97aRk6Z3Gbbt9LsVd7+6x0G2eyzszGNiJI9FSDKD7sY1Q0+27fSnFW8c31ZCGX2aHN9m0OnXgNuvjGtvyf5MjtP0TWrst4xzRCDt3JEn/dqwR/rXE1LPodKgg6RgONctAtoVYQlA21kJ5z9QMygFC6RAHp9aLEuK10nxBikRleL1bhHok2E4RXcg5BQ0XsWMbWcWmrbVpPBsBcOQAhEpEHb6OVQkBSOt21dHsm23zxukeLQxaFrxtu3XMBuL9fRAKqP58bAsyY5LVo1btMabB9yW/L5QodPv4R75GbYsS4rVSNEaKVbXWIz7kQyWRqB5/yYlx3e83vkdrT2eZlAKFzv9GizMzz5NSsSkaLUUrXXGsp/xkxyrGRdOzMY+td0garsB1UqkfW3FMZQaz8nBcPbOFXu2SB884/RFn1HS+FOk3sPbvl8gT+VxRQOgU7Hcoql2R9uCgBlyigkz6IaQtII3WZwli4hEzH8RYRhOQRaMOMVuUZ/8DndNWZZTpCWiUiIuVW9zvvfLMNL62C2Ek7NgTftcdlrh1tqQYbgBL+IUbmZx55yJsBJOAInVOmPdV7FsNJkxCO4lZFtpsxAx//+eDEMKFTvju7ive2zzWCLmBLv6Kim+tXUB1kzvzyYXL5IzkMmLFa0dw5Ib7HpIReVOCMn3sZuIOefleFSq2e5csPArdV52g17dLmc/sVqp6t9SzQ7CErq0TlQhAMhbZkAq7t32/SSXtqWW/6QtqUpd6Q41FpuBUP4XgtlmmpJZ6FyFlyS1st+t5PK1pkuCmixfM90r38FIY4HfnZevmQEp0sPZVOHvuZbVZMYg3ry/DUMygk5/m8HGiwf5XpS3RSAkFfZ0Nr9SSzHTzhlNl8WaptOn6TNl3WUMB9zZ77ZMKNm2c0EmuYSwoFQqGySVDZFGHSuV75u15gL5iLAEIH+YpmR2giVlXYFpSiavNehQpinJfW0PsiN9tk6M53ZhGO7yPffcXFAqffVaZ6a4K4d4wEVYAgAAQMt1htdJAlnSDeagAQAAAMA/whIAAAAAeCAsAQAAAIAHwhIAAAAAeCAsAQAAAIAHwhIAAAAAeCAsAQAAAIAHwhIAAAAAeCAsAQAAAIAHwhIAAAAAeAjmugEA0ClZlmTFpERMshOSbUm2nfkYw5BkSIYpmQHnq2FKRsC9beSk6V2GbTv9bsXdr+5xkO0eCzvzmBjJYxGSzKC7cc3Qk207/WnFG8e3lVBGnybHtxl0+jXg9ivj2luyP5PjNH2TGvst4xwRSDt35Em/RmukD5ZIPSul/uOkAKUkujZGOIC2sxJSwx6nUAoW5Pd/npYlxeukeIOUiErxercI9MkwnKI7EHIKGq9ixrYzC03balJ4toJhSIGIFAg7/RwqkoKR1u2rI9m22+cNUjzaGDSteNv2a5iNxXp6IJXR/HhYlmTHJavGLVrjzQNuS35fqNDp93CP/AxbliXFapyiNlbXWIz7kQyWRqB5/yYlx3e83vkdrT2eZlAKFzv9GizMzz5NSsSkaLUUrXXGsp/xkxyrGRdOzMY+td0garsB1UqkfW3FMZQaz8nBcPbOFVWfSv98zGlbxf7SmFlS+bC27xfIU3lc0QDoVCy3aKrd0bYgYIacYsIMuiEkreBNFmfJIiIR819EGIZTkAUjTrFb1Ce/w11TluUUaYmolIhL1duc7/0yjLQ+dgvh5CxY0z6XnVa4tTZkGG7AiziFm1ncOWcirIQTQGK1zlj3VSwbTWYMgnsJ2VbaLETM/78nw5BCxc74Lu7rHts8log5wa6+SopvbV2ANdP7s8nFi+QMZPJiRWvHsOQGux5SUbkTQvJ97CZiznk5HpVqtjsXLPxKnZfdoFe3y9lntEbasU6q2kxYQpfWiSoEAHnLDEjFvdu+n+TSttTyn7QlVakr3aHGYjMQyv9CMNtMUzILnavwkqRW9ruVXL7WdElQk+VrpnvlOxhpLPC78/I1MyBFejibKvw917KazBjEm/e3YUhG0OlvM9h48SDfi/K2CISkwp7O5ldqKWbaOaPpsljTdPo0faasu4zhgDv73ZYJJdt2LsgklxBGSqSS/lLJQGnUsVLFl7LWXCAfEZYA5A/TlMxOsKSsKzBNyQznuhXdi2lKcl/bg+xIn60T47ldGIa7fM89NxeUSl+9xlnW15VDPOAiLAEAAKDlwsW5bgHQYbrBHDQAAAAA+EdYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8JAXYemOO+7QPvvso4KCAh166KF64403PvfxDz30kEaPHq2CggKNGTNGTz75ZAe1FAAAAEB3kfOw9OCDD2revHm66qqrtGrVKo0bN04zZszQtm3bPB//6quv6pRTTtE555yjN998UzNnztTMmTO1Zs2aDm45AAAAgK4smOsG3HLLLTrvvPN01llnSZLuuusuPfHEE7r33nv1k5/8pNnjFy5cqGOOOUaXXnqpJOnaa6/VkiVLdPvtt+uuu+7q0LYDcCz/93JduOxC1SXqNLh4sL476rvqGemp/kX9NaZ8jIqKinLdxE5hR+0OfbjzQ+1s2KlttdsUt+KqidcoYSUkSQEzoOJgsQqDhepV0Ev9i/urd2Fv9Qn1oY/3ora2VptrN2v9nvXaVrtNu6O7FbNiGX0aMkMKm2GVRcpUUVShknAJ/doKtbW12hHboc17NqfGcG28dq/9XRgsVHGoWMWhYpWESjS8fLj6FPXJ8V/ROW3cvVEf7/5Yn9V/pl0NuzL6PdnnRcGi1HmZvgZaLqdhKRqNauXKlZo/f37qPtM0NX36dC1fvtzzOcuXL9e8efMy7psxY4YeffRRz8c3NDSooaEhdbuqqqrtDQeQcu2r1+rPH/w5dfvTmk/1wWcfqDZRq5pojaqj1UookfEcy7YkSaax98nt8ki5KksrNbB4oA4ecLD277N/+/wBWVJbW6tXt76qf2z9h9Z+tlZ18boWPS/ZF5IUMkLqEe6hwmChioJFCpkhFYQKFDScU3V9ol7/qf2PYlZMtfFaVceq1ZBoUMyK+WprQAH1LeqrUeWjdGTlkdqnaJ+8DgUbd2/UG5vf0Kqtq7ShakOz8fRFIoGISkIlKg2XqiBUoMJAoYJmUHEr7gRSu0b1sXrtie9RfaxeUSvqu18DCmhY2TCNrxivQwYcoiFlQ/z+mTlTW1urj2s/1obdG/Rp9afaWLVR2+u2a3fDbt99nSzKC4OF6hHqoXAwrJAZSo3h9P5O2AnVJmoVjUdVn6hv1ViWpLJwmQaVDNKg4kEaWjpU4/uN71RB4MMdH+r1ra9rw54N2l63XVurt2b0u2Vbn3uulJx+LwmVqCBUoKJAkYrDxQoaQQXNoBriDaqxa/TJnk8yzssxO/aF+/USCUTUr6ifeoR76MwDzuxUYx1oDcO2bTtXv/zTTz/VoEGD9Oqrr2rSpEmp+y+77DK98MILev3115s9JxwO6/e//71OOeWU1H2//vWvdfXVV2vr1q3NHv/zn/9cV199dbP7d+/erdLS0iz9JQAAAAC6mpy/Zqm9zZ8/X7t3705tu3bt0rZt21RSUpLrpgEAAADIYzldhtenTx8FAoFmM0Jbt25V//79PZ/Tv39/X4+PRCKKRCLZaTAAAACAbiOnM0vhcFgHHXSQli5dmrrPsiwtXbo0Y1leukmTJmU8XpKWLFmy18cDAAAAQGvk/N3w5s2bpzPOOEMTJ07UIYccottuu001NTWpd8ebPXu2Bg0apAULFkiSLrzwQh155JG6+eabdfzxx+uBBx7QihUrdM899+TyzwAAAADQxeQ8LJ188snavn27fvazn2nLli0aP368nnrqKfXr10+StHHjRplm4wTY5MmTtXjxYl155ZW6/PLLNXLkSD366KM68MADc/UnAAAAAOiCcvpueEA2nHnmmdq1a5ceffRRnXnmmfr973/f7DEffPCBRowYkYPWAQCA7mRvdUkoFNKQIUM0e/ZsXX755QoGcz5ngRbgKKHLOeaYY3Tfffdl3Ne3b98ctQYAAHRnybqkoaFBTz75pObMmaNQKJTxOaPIX13+rcPR/UQiEfXv3z9jCwQCuW4WAADohpJ1ydChQ3X++edr+vTp+n//7//lulloIcISAAAA0EEKCwsVjUZz3Qy0EGEJXc7jjz+uHj16pLZZs2blukkAAKCbs21bzz77rJ5++mkdddRRuW4OWojXLKHLmTZtmu68887U7eLi4hy2BgAAdGfJi7ixWEyWZenUU0/Vz3/+81w3Cy1EWEKXU1xczDvfAQCAvJC8iBsOhzVw4EDeBa+T4WgBAAAA7YSLuJ0br1kCAAAAAA+EJQAAAADwYNi2bee6EQAAAACQb5hZAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8EBYAgAAAAAPhCUAAAAA8PD/AZLQomUdsQpPAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import demes\n", + "import demesdraw\n", + "\n", + "graph = demes.load(\"demes.yaml\")\n", + "demesdraw.tubes(graph)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "896122bf", + "metadata": {}, + "outputs": [], + "source": [ + "import msprime\n", + "\n", + "graph = demes.load(\"demes.yaml\")\n", + "demography = msprime.Demography.from_demes(graph)\n", + "\n", + "ts1 = msprime.sim_ancestry(samples={\"IF\": 2}, demography=demography, random_seed=12)\n", + "ts1 = msprime.sim_mutations(ts, rate=1e-6, random_seed=12)\n", + "\n", + "ts1.dump(\"ts.trees\")\n", + "\n", + "ts2 = msprime.sim_ancestry(samples={\"IF\": 2}, demography=demography, random_seed=12)\n", + "ts2 = msprime.sim_mutations(ts, rate=1e-6, random_seed=12)\n", + "ts2.dump(\"test.trees\")\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "2e1c4862", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[32m2025-06-19 10:23:47.727\u001b[0m | \u001b[33m\u001b[1mWARNING \u001b[0m | \u001b[36mphlash.mcmc\u001b[0m:\u001b[36m_check_jax_gpu\u001b[0m:\u001b[36m23\u001b[0m - \u001b[33m\u001b[1mDetected that Jax is not running on GPU; you appear to have CPU-mode Jax installed. Performance may be improved by installing Jax-GPU instead. For installation instructions visit:\n", + "\n", + "\thttps://github.com/google/jax?tab=readme-ov-file#installation\n", + "\u001b[0m\n", + "\u001b[32m2025-06-19 10:23:47.731\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mphlash.mcmc\u001b[0m:\u001b[36mfit\u001b[0m:\u001b[36m88\u001b[0m - \u001b[1mLoading data\u001b[0m\n", + "\u001b[32m2025-06-19 10:23:47.733\u001b[0m | \u001b[33m\u001b[1mWARNING \u001b[0m | \u001b[36mphlash.data\u001b[0m:\u001b[36minit_mcmc_data\u001b[0m:\u001b[36m521\u001b[0m - \u001b[33m\u001b[1mThe chunk size is 0, which is less than 10 times the overlap (500).\u001b[0m\n", + "100%|██████████| 2.00/2.00 [00:01<00:00, 1.55bp/s]\n", + "\u001b[32m2025-06-19 10:23:49.309\u001b[0m | \u001b[34m\u001b[1mDEBUG \u001b[0m | \u001b[36mphlash.mcmc\u001b[0m:\u001b[36mfit\u001b[0m:\u001b[36m122\u001b[0m - \u001b[34m\u001b[1mMinibatch size: 1\u001b[0m\n", + "\u001b[32m2025-06-19 10:23:49.311\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mphlash.mcmc\u001b[0m:\u001b[36mfit\u001b[0m:\u001b[36m155\u001b[0m - \u001b[1mScaled mutation rate Θ=nan\u001b[0m\n", + "\u001b[32m2025-06-19 10:23:49.322\u001b[0m | \u001b[33m\u001b[1mWARNING \u001b[0m | \u001b[36mphlash.kernel\u001b[0m:\u001b[36mget_kernel\u001b[0m:\u001b[36m15\u001b[0m - \u001b[33m\u001b[1mError when loading GPU code, falling back on pure JAX implmentation. This will be **much slower**. Error was: No module named 'nvidia'\u001b[0m\n" + ] + }, + { + "data": { + "text/html": [ + " \n", + " " + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": {}, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Fitting model: 0%| | 0/1000 [00:04 \u001b[39m\u001b[32m8\u001b[39m results = \u001b[43mphlash\u001b[49m\u001b[43m.\u001b[49m\u001b[43mfit\u001b[49m\u001b[43m(\u001b[49m\n\u001b[32m 9\u001b[39m \u001b[43m \u001b[49m\u001b[43mdata\u001b[49m\u001b[43m=\u001b[49m\u001b[43m[\u001b[49m\u001b[43msample\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m 10\u001b[39m \u001b[43m \u001b[49m\u001b[43mmutation_rate\u001b[49m\u001b[43m=\u001b[49m\u001b[32;43m1e-6\u001b[39;49m\u001b[43m)\u001b[49m\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/Documents/SMCk/.venv/lib/python3.12/site-packages/phlash/mcmc.py:284\u001b[39m, in \u001b[36mfit\u001b[39m\u001b[34m(***failed resolving arguments***)\u001b[39m\n\u001b[32m 281\u001b[39m \u001b[38;5;28;01massert\u001b[39;00m jnp.isfinite(x).all()\n\u001b[32m 282\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m x\n\u001b[32m--> \u001b[39m\u001b[32m284\u001b[39m state = \u001b[43mjax\u001b[49m\u001b[43m.\u001b[49m\u001b[43mtree\u001b[49m\u001b[43m.\u001b[49m\u001b[43mmap\u001b[49m\u001b[43m(\u001b[49m\u001b[43mf\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mstate1\u001b[49m\u001b[43m)\u001b[49m\n\u001b[32m 285\u001b[39m _particles = state.particles\n\u001b[32m 286\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m test_data \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;129;01mand\u001b[39;00m i % \u001b[32m10\u001b[39m == \u001b[32m0\u001b[39m:\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/Documents/SMCk/.venv/lib/python3.12/site-packages/jax/_src/tree.py:155\u001b[39m, in \u001b[36mmap\u001b[39m\u001b[34m(f, tree, is_leaf, *rest)\u001b[39m\n\u001b[32m 115\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34mmap\u001b[39m(f: Callable[..., Any],\n\u001b[32m 116\u001b[39m tree: Any,\n\u001b[32m 117\u001b[39m *rest: Any,\n\u001b[32m 118\u001b[39m is_leaf: Callable[[Any], \u001b[38;5;28mbool\u001b[39m] | \u001b[38;5;28;01mNone\u001b[39;00m = \u001b[38;5;28;01mNone\u001b[39;00m) -> Any:\n\u001b[32m 119\u001b[39m \u001b[38;5;250m \u001b[39m\u001b[33;03m\"\"\"Maps a multi-input function over pytree args to produce a new pytree.\u001b[39;00m\n\u001b[32m 120\u001b[39m \n\u001b[32m 121\u001b[39m \u001b[33;03m Args:\u001b[39;00m\n\u001b[32m (...)\u001b[39m\u001b[32m 153\u001b[39m \u001b[33;03m - :func:`jax.tree.reduce`\u001b[39;00m\n\u001b[32m 154\u001b[39m \u001b[33;03m \"\"\"\u001b[39;00m\n\u001b[32m--> \u001b[39m\u001b[32m155\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mtree_util\u001b[49m\u001b[43m.\u001b[49m\u001b[43mtree_map\u001b[49m\u001b[43m(\u001b[49m\u001b[43mf\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtree\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m*\u001b[49m\u001b[43mrest\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mis_leaf\u001b[49m\u001b[43m=\u001b[49m\u001b[43mis_leaf\u001b[49m\u001b[43m)\u001b[49m\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/Documents/SMCk/.venv/lib/python3.12/site-packages/jax/_src/tree_util.py:361\u001b[39m, in \u001b[36mtree_map\u001b[39m\u001b[34m(f, tree, is_leaf, *rest)\u001b[39m\n\u001b[32m 359\u001b[39m leaves, treedef = tree_flatten(tree, is_leaf)\n\u001b[32m 360\u001b[39m all_leaves = [leaves] + [treedef.flatten_up_to(r) \u001b[38;5;28;01mfor\u001b[39;00m r \u001b[38;5;129;01min\u001b[39;00m rest]\n\u001b[32m--> \u001b[39m\u001b[32m361\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mtreedef\u001b[49m\u001b[43m.\u001b[49m\u001b[43munflatten\u001b[49m\u001b[43m(\u001b[49m\u001b[43mf\u001b[49m\u001b[43m(\u001b[49m\u001b[43m*\u001b[49m\u001b[43mxs\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mxs\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[38;5;28;43mzip\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m*\u001b[49m\u001b[43mall_leaves\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/Documents/SMCk/.venv/lib/python3.12/site-packages/jax/_src/tree_util.py:361\u001b[39m, in \u001b[36m\u001b[39m\u001b[34m(.0)\u001b[39m\n\u001b[32m 359\u001b[39m leaves, treedef = tree_flatten(tree, is_leaf)\n\u001b[32m 360\u001b[39m all_leaves = [leaves] + [treedef.flatten_up_to(r) \u001b[38;5;28;01mfor\u001b[39;00m r \u001b[38;5;129;01min\u001b[39;00m rest]\n\u001b[32m--> \u001b[39m\u001b[32m361\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m treedef.unflatten(\u001b[43mf\u001b[49m\u001b[43m(\u001b[49m\u001b[43m*\u001b[49m\u001b[43mxs\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mfor\u001b[39;00m xs \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mzip\u001b[39m(*all_leaves))\n", + "\u001b[36mFile \u001b[39m\u001b[32m~/Documents/SMCk/.venv/lib/python3.12/site-packages/phlash/mcmc.py:281\u001b[39m, in \u001b[36mfit..f\u001b[39m\u001b[34m(x)\u001b[39m\n\u001b[32m 280\u001b[39m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34mf\u001b[39m(x):\n\u001b[32m--> \u001b[39m\u001b[32m281\u001b[39m \u001b[38;5;28;01massert\u001b[39;00m jnp.isfinite(x).all()\n\u001b[32m 282\u001b[39m \u001b[38;5;28;01mreturn\u001b[39;00m x\n", + "\u001b[31mAssertionError\u001b[39m: " + ] + } + ], + "source": [ + "import phlash\n", + "\n", + "sample = phlash.contig(\"ts.trees\", samples=[(0, 1)])\n", + "test = phlash.contig(\"test.trees\", samples=[(0, 1)])\n", + "\n", + "#results = phlash.fit([sample], test_data=test)\n", + "\n", + "results = phlash.fit(\n", + " data=[sample],\n", + " mutation_rate=1e-6)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}